An Infinite Dimensional Analysis of Kernel Principal Components

TL;DR

This paper develops a new kernel PCA framework for nonlinear data dimension reduction, extending classical PCA with manifold and feature space transforms, and provides theoretical insights into optimal kernel choices.

Contribution

Introduces a novel kernel PCA method with theoretical analysis for optimal Gaussian kernel selection and extends probabilistic Karhunen-Loève transforms to nonlinear settings.

Findings

Proves new theorems for data-dimension reduction.

Identifies conditions for optimal Gaussian kernel choice.

Enhances digital image representation and compression.

Abstract

We study non-linear data-dimension reduction. We are motivated by the classical linear framework of Principal Component Analysis. In nonlinear case, we introduce instead a new kernel-Principal Component Analysis, manifold and feature space transforms. Our results extend earlier work for probabilistic Karhunen-Lo\`eve transforms on compression of wavelet images. Our object is algorithms for optimization, selection of efficient bases, or components, which serve to minimize entropy and error; and hence to improve digital representation of images, and hence of optimal storage, and transmission. We prove several new theorems for data-dimension reduction. Moreover, with the use of frames in Hilbert space, and a new Hilbert-Schmidt analysis, we identify when a choice of Gaussian kernel is optimal.

Equations

$$X_A=\sum _{n=0}^{\infty }\left({\int }_A{f}_{n}d\mu \right)Z_n,A\in {\mathcal{B}}_{fin}$$

$$\sum _{i,j=1}^N\overline{c_i}{c}_{j}K\left(x_i,x_j\right)\ge 0$$

$$K\left(x,y\right)={⟨\Phi \left(x\right),\Phi \left(y\right)⟩}_{\mathcal{H}\left(K\right)}.$$

$$f\left(x\right)={⟨K_x,f⟩}_{\mathcal{H}\left(K\right)},$$

$$\text{span}\left\{K_x:=K\left(\cdot ,x\right)\right\}$$

$${⟨\sum c_i{K}_{x_i},\sum d_j{K}_{x_j}⟩}_{\mathcal{H}\left(K\right)}:=\sum \overline{c_i}{d}_{j}K\left(x_i,x_j\right).$$

$${\Vert X\Vert }_{HS}^2=\sum {\left|x_{ij}\right|}^2=tr\left(X^{\ast }X\right),$$

$$⟨X,Y⟩=tr\left(X^{\ast }Y\right).$$

$$HS\left(\mathcal{H}\right)\ni \left|u⟩⟨v\right|⟶u\otimes \overline{v}\in \mathcal{H}\otimes \overline{\mathcal{H}},$$

$${\Vert \left|u⟩⟨v\right|\Vert }_{HS}^2={\Vert u\Vert }_{\mathcal{H}}^2{\Vert v\Vert }_{\mathcal{H}}^2,$$

$$\mathcal{H}\left(K\right)\ni f\stackrel{T}{\to }{\left(f\left(x\right)\right)}_{x\in X}\in L^2\left(\mu \right)$$

$$Q_{\psi ,\beta }\left(f\right)={\Vert \psi -Tf\Vert }_{L^2\left(\mu \right)}^2+\beta {\Vert f\Vert }_{\mathcal{H}\left(K\right)}^2$$

$$L^2\left(\mu \right)\stackrel{T^{\ast }}{\to }\mathcal{H}\left(K\right)$$

$${⟨F,T^{\ast }\psi ⟩}_{\mathcal{H}\left(K\right)}={⟨Tf,\psi ⟩}_{L^2\left(\mu \right)}={\int }_X\overline{f\left(s\right)}\psi \left(s\right)d\mu \left(s\right).$$

$$\underset{f\in \mathcal{H}\left(K\right)}{\inf }{Q}_{\psi ,\beta }\left(f\right)$$

$$F={\left(\beta I+T^{\ast }T\right)}^{-1}{T}^{\ast }\psi$$

$$\frac{d}{d\varepsilon}\big{|}_{\varepsilon=0}Q_{\psi,\beta}\left(f_{\varepsilon}\right)=2\Re\left\langle h,\left(\beta I+T^{*}T\right)F-T^{*}\psi\right\rangle_{\mathcal{H}\left(K\right)}=0$$

$$Q\left(f\right)=\sum _{i=1}^m\underset{\text{least square}}{{\left|f\left(t_i\right)-y_i\right|}^2}+\beta \underset{\text{penalty form}}{{\Vert f\Vert }_{\mathcal{H}\left(K\right)}^2},f\in \mathcal{H}\left(K\right).$$

$$\begin{array}{rl}T:\mathcal{H}\left(K\right)& ⟶{\mathbb{R}}^m\simeq l_m^2,\text{mbox}and\\ T^{\ast }:l_m^2& ⟶\mathcal{H}\left(K\right)\end{array}$$

$$Tf$$

$$T^{\ast }y$$

$${⟨T^{\ast }y,f⟩}_{\mathcal{H}\left(K\right)}={⟨y,Tf⟩}_{l_m^2},\forall f\in \mathcal{H}\left(K\right),\forall y\in l_m^2;$$

$$\mu =\sum _{i=1}^m{\delta }_{t_i}={\delta }_{\left\{t_i\right\}}$$

$$T^{\ast }Tf$$

$$T{T}^{\ast }\vec{y}$$

$$K_m={\left(K\left(t_i,t_j\right)\right)}_{i,j=1}^m=\left(\begin{array}{cccc}K\left(t_1,t_1\right)& \cdots & \cdots & K\left(t_1,t_m\right)\\ K\left(t_2,t_1\right)& \cdots & \cdots & K\left(t_2,t_m\right)\\ ⋮& & & \\ K\left(t_m,t_1\right)& & & K\left(t_m,t_m\right)\end{array}\right).$$

$$Q_{\beta ,\left\{t_i\right\},\left\{y_i\right\}}\left(f\right)=\sum _{i=1}^m{\left|y_i-f\left(t_i\right)\right|}^2+\beta {\Vert f\Vert }_{\mathcal{H}\left(K\right)}^2,f\in \mathcal{H}\left(K\right).$$

$$F\left(\cdot\right)=\sum_{i=1}^{m}\left(K_{m}+\beta I_{m}\right)_{i}^{-1}K\left(\cdot,t_{i}\right)\>\mbox{on $X$;}$$

$$\beta F+T^{\ast }TF=T^{\ast }y,$$

$$\beta F\left(\cdot \right)=\sum _{i=1}^m\left(y_i-F\left(t_i\right)\right)K\left(\cdot ,t_i\right)$$

$$\beta \vec{F}=K_m\left(\vec{y}-\vec{F}\right)$$

Full text

1. 1 Introduction

2. 1.1 Reproducing kernels and operator theory

3. 2 Kernel PCA

4. 2.1 Application to Optimization

5. 2.2 The Case of Gaussian Fields

6. 2.3 Optimization and Frames

7. 2.4 The Dual Problem

8. 2.5 Feature Selection

9. 2.6 Dynamic PCA

1. 2.1 Spectral clustering and linear decision boundary via Gaussian kernels and KPCA.
2. 2.2 Non-linear data, and detection of rotation angles via KPCA.

1. Introduction

Recently a number of new features of principal component analysis (PCA) have lead to exciting and new improved dimension reduction (DR). See e.g., [BN03, GGB18, JHZW19, GJ06, Bis06, Bis13, ZBB04, AFS18, VD16]. In general DR refers to the process of reducing the number of random variables under consideration in such areas as machine learning, statistics, and information theory. Within machine learning, it involves both the steps of feature selection and feature extraction. In the present paper, we shall consider linear as well as non-linear data models. The linear case arises naturally in principal component analysis (PCA). See [Son08, JS07]. Here one starts with a linear mapping of the given data into a suitable lower-dimensional space. However, this must be done in such a way that the variance of the data in the low-dimensional representation is maximized. As for the variance, we study both covariance, and the correlation matrix for the underlying data. The eigenvectors that correspond to the largest eigenvalues (the principal components) are then used in a construction of a large fraction of the initial variance, i.e., that which corresponds to the original data. The first few eigenvectors are typically interpreted in terms of the large-scale physical behavior of a particular system; and will retain the most important variance features.

In nonlinear settings, principal component analysis can still be adapted, but now by means of suitable kernel tricks. In some applications, instead of starting from a fixed kernel, the optimization will instead try to learn, or adapt, for example with the use of semidefinite programs [WGY10, BM03]. The most prominent such a technique is known as maximum variance unfolding (MVU) [ACP13, SBXD06].

The aim of this paper is to set the stage for further developments of the linking of the mathematics of kernel-PCA (KPCA) with real-life non-linear data-sets, coming from a host of new developments in the current applied literature. A key focus for our work is a study of two complementary directions: on one side, the mathematics of kernel theory, and on the other, empirical (here, numerical experiments starting with databases from real-world situations). For the purpose, we develop (Section 2 ) particular kernel tools with view to KPCA via statistical simulation, as well as corresponding analyses of kernel based choices of feature spaces and feature maps.

In recent years, the subject of kernel-principal component analysis (KPCA), and its applications, has been extensively studied, and made progress in diverse directions. There is a variety of themes covered in the literature dealing with applications of kernel-PCA tools is vast, each of some relevance to our present theme. A short user guide: the paper [GGB18] deals with sensitivity issues, [JHZW19] direct integration algorithms, [GJ06] image processing, [Bis06, Bis13, AFS18] pattern recognition and machine learning, [ZBB04] statistical properties, and finally [VD16] on selection of efficient parameterizations. In addition to earlier papers by the co-authors [Son08, JS07], we also include here a partial list of other relevant and current citations; see e.g., the papers mentioned above, [LHN18, Raj18, WGLP19, RS99, GGB18, DWGC18, THH19, LLY*+*19, MMP19, CLS*+*19, GK19].

1.1. Reproducing kernels and operator theory

In our present paper we stress, and motivate, the most general framework for choices of positive definite kernels, their corresponding RKHSs, and families of compatible measure spaces. These considerations encompass more general contexts for Principal Component Analysis, going beyond classical cases of kernel-based PCA in the literature. In particular, this very general and non-linear framework for PCA data-dimension reduction goes beyond existing, more classical, considerations admitting easy realizations in $\mathbb{R}^{n}$, classical kernels, and classical Fourier tools. Indeed, there are existing parallel considerations, in the literature involving only kernel analyses for such special cases.

However, existing non-linear realizations of big and un-structured data motivate choices from a much wider family of RKHSs, and measure spaces. This further entails infinite-dimensional stochastic analysis frameworks which do not admit realizations in $\mathbb{R}^{n}$, nor classical PDE/Fourier techniques. Examples of uses, in the literature, of general RKHSs in stochastic approaches include online learning with Markov sampling, data mining (e.g., in finance, intelligence, hidden data, tele-communication, energy, unstructured information integrated with traditional structured data), as well as a diverse host of different investigations, see e.g., [ACP13, AHR*+*20, AJS14, Bis06, BM03, BTA04, CS02, CWG19, JT16a, JT18a, JT18b, Mok07, SZ09a, SZ09b, ZBB04, JT17, Agg18, Hog20]. Moreover “big data” further dictate diverse, large, infinite-dimensional, and extremely non-linear data-sets which do not lend themselves to feature analysis with choices of anyone of the simple kernels with easy realizations in $\mathbb{R}^{n}$. Examples of this includes many public available datasets on Global Health Observatory, U.S. Census Bureau, NIH Data Sharing Repositories, Kaggle, etc. Such cases of course dictate new and diverse choices from a much wider family of kernels and RKHSs. Encoding of meaningful programs for corresponding non-linear data-dimension reduction algorithms, while feasible, is a gigantic task, going far beyond the scope of our present paper. It is the subject of future projects involving the present authors, as well as others, see e.g., the papers cited above.

In addition, we call attention to the 2018 paper [YTLM*+*18] by Yamada et al. It is motivated by existing and new biological data. For the purpose, the authors introduce a novel Hilbert-Schmidt independence criterion. They apply it to specific high dimensional, and large-scale datasets, as arising in ultra-high dimensional biological data. The authors obtain classification into phenotypes, module expression in human prostate cancer patients, and detection of enzymes in protein structures. The 2019 paper [CLW*+*19] by Chen et al deals with a new class of kernel based nonlinear clustering algorithms. They study locality structures in machine learning. Also motivated by specific large-scale datasets is the 2020 paper [UTK20]. This paper introduces new hyper-parameter initialization tools for kernel-based regression analysis. Yet different but related kernel based technology, and specific data sets, are further covered in the papers [FZ05], [WSS04], [MMR*+*01], and [DSR*+*21].

Motivated by related applications, we focus on non-traditional positive definite kernels and the associated RKHSs, in particular, when the family of feature spaces may be chosen in the form $L^{2}\left(\mu\right)$. This raises the question of which measures $\mu$ are right for a particular kernel $K$ and its RKHS.

In the use of probabilistic tools in big data, such as Monte Carlo simulation (see e.g., (1.1) below), a main issue is reduction of the required sample size; and still maintain (i) acceptable small noise in output, (ii) small error term estimates, as well as (iii) reduction of computation time. These issues are especially relevant for modern image processing. As a consequence, more traditional kernel tools must therefore be adapted. Both our present paper, and an earlier one [JT22], serve to do this by shifting the focus in the design of optimization algorithms to that of making choices of measurable partitions, as opposed to working directly with points.

For instance, when $K=K^{\left(\mu\right)}$ is defined on a $\sigma$-finite measure space $\left(V,\mathscr{B},\mu\right)$, where $K^{\left(\mu\right)}\left(A,B\right)=\mu\left(A\cap B\right)$, for all $A,B$ in $\mathscr{B}_{fin}=\left\{A\in\mathscr{B}:\mu\left(A\right)<\infty\right\}$, the associated RKHS $\mathscr{H}(K^{\left(\mu\right)})$ is shown to be a Hilbert space of signed finite measures $F$ on $\left(V,\mathscr{B}\right)$, such that $F\ll\mu$ and $\left\|F\right\|_{\mathscr{H}(K^{\left(\mu\right)})}=\left\|dF/d\mu\right\|_{L^{2}\left(\mu\right)}$; see [JT22, Thms 3.2, 3.4]. Moreover, in this measurable setting, the classical Karhunen-Loève decomposition takes the following form instead:

[TABLE]

where $X_{A}$ is a centered Gaussian process indexed by $\mathscr{B}_{fin}$, $\left\{f_{n}\right\}_{n\in\mathbb{N}_{0}}$ is an orthonormal basis (ONB) in $L^{2}\left(\mu\right)$, and $\left\{Z_{n}\right\}_{n\in\mathbb{N}_{0}}$ is an i.i.d. $N\left(0,1\right)$ system. These new kernels encompass realistic cases when random measurements (random variables) might not be performed at points, but instead at sample sets that are selected from suitable choices of sigma-algebras. For more details on the more general context of Parseval frames in the measure category, see 2.2 below, and [JT22, Cor 3.13 and sect. 5].

Main results. The presentation of our main results is organized as follows: 2 covers multiple interrelated tools, playing a key role in our main results. The first of these tools is our use of kernels (also called reproducing kernels) and their associated Hilbert spaces, called reproducing kernel Hilbert spaces (RKHSs.) The use of these RKHSs yields powerful tools, and they allow explicit realizations of a variety of embeddings of non-linear structures into linear ones (Hilbert spaces). This is useful in turn for our solution to optimization questions, see 2.14 and 2.16. Our results are tested on datasets, see Figures 2.1 and 2.2.

Organization. The goal of this paper is to extend our previous results on Karhunen-Loève transform to a nonlinear setting by means of kernel-principal component analysis (KPCA). This paper is organized as follows:

In 2, we show our main results using KPCA on nonlinear data. See 2.9. The latter is for rank-1 projections, but is then extended to PCA selection, and algorithms, in subsequent results. Our focus is the non-linear case. Indeed, the focus in 2 is nonlinear data dimension reduction, and the corresponding kernels, the core of our paper. In particular, our 2.11 deals with kernel PCA for nonlinear data dimension reduction (see Examples 2.17 and 2.18). In the remaining part of 2, we turn to the case of Gaussian kernels.

KPCA have found wide applications. The focus of our current paper is to present a framework of operator in suitable Hilbert space, and an associated spectral theory. Even though the applications we include here are presented in finite dimensional setting, most of our results extend to infinite dimensional spaces as well. Nonetheless, for use in recursion schemes, the finite-dimensional case is most relevant.

Our main results deal with algorithms for optimization in maximal variance, and dimension reduction-problems, from PCA. The points where our results go beyond the earlier literature includes the following list of four closely interrelated items: (i) Our design and use of kernel tricks (specific reproducing kernel Hilbert spaces (RKHSs) integrated into new designs of PCA-tools for image-analysis, thus serving as practical tools in dimension-reduction algorithms; (ii) Our combination of RKHS-tools with spectral theory and Hilbert space frame-estimates; — this serves in turn to make more precise both specific numerical PCA-algorithms, and their error estimates; (iii) Our use of RKHSs in creating explicit and practical embeddings of non-linear data into suitable linear spaces (feature spaces); — applications to optimization (2.1); and (iv) A new dynamical PCA-analysis (2.6)

Our results in 2 include 2.9, 2.11, and 2.16; each of which are formulated and proved in a general setting of kernel analysis; hence in a non-linear framework of feature selection.

Our main purpose is the identification and proof of rigorous results from infinite-dimensional analysis which are needed in the solution of optimization problems in the general framework of Principal Component Analysis with an emphasis on the role of PCA-algorithms for selection of features in infinite dimensions. Main points in the paper include our infinite-dimensional optimization results, Theorems 2.5, 2.6, 2.11, Remark 2.15, and Corollaries 2.16, and 2.20.

2. Kernel PCA

Comparison of KPCA to standard PCA (brief sketch.) The standard PCA always finds linear principal components. It serves to represent a given large data set into a suitable choice of principal components for lower dimension. But PCA will fail to find good representative directions when applied to most non-linear data of interest. For non-linear and more tricky data sets, one turns instead to Kernel PCA (KPCA). It has the following advantages and features (not a complete list): (i) KPCA will perform PCA but in a new space. (ii) It uses diverse kernel tricks in order to still find principal components, but they will be in a different space (typically in a higher dimensional space.) (iii) KPCA will use Hilbert space geometry in order to still find new directions based on families of kernel-matrices. Via geometric algorithms, KPCA will then extract corresponding eigenvalues (corresponding to a suitable number of observations.) However, (iv) the computational complexity for KPCA, dictated by extraction of principal components, will take more time, as compared to standard PCA. The list of papers on KPCA and applications to diverse sets of non-linear data is long; some recent ones are [HKF*+*20, ZSS19, Law05, AH11, HLMS04, Nai17, FCB].

In this section, we make precise, and show that the choice of Gaussian kernel is “optimal” for KPCA in capturing maximal variance (i.e., maximal variability of data). There is an extensive literature on PCA and KPCA. Our approach and emphasis are different. For instance, the paper [RT01] deals with certain Least Squares Regression problems with the use of a particular (small) class of Reproducing tools involving Kernel Hilbert Space. While there are connections to our present themes, our present results are much more general, offering in particular specific spectral theoretic, and stochastic tools, penalty terms, especially relevant for digital images, and for dynamical PCA. The focus of [SSM98] is the design of certain Nonlinear Component Analysis as a particular class of Kernel Eigenvalue Problems. By comparison, our present results and applications have a different focus, use different tools both as regards to infinite-dimensional analysis, optimization, and Gaussian fields; as we outline below. This paper [KW71] deals with a particular class of splines (numerical analysis). This in turn entails a class of kernels. But the focus of [KW71] is quite different from ours. The paper [WMd97] deals with matrix algorithms close in spirit to our present paper, but our focus and applications are global in nature, and are very different as we outline below. The book [CST00] deals with kernel methods and some of their applications to support vector machines. This is an exciting application of feature selections, and learning algorithms, but the aim and the tools developed in this book are quite different from our present focus.

This allows us to obtain extensions of some PCA-results from [JS07]. We shall focus here on their use in kernel PCA (KPCA).

PCA is used in data dimension reduction on linear case. However, this cannot be done on nonlinear case and thus kernel principal component analysis (KPCA) is used for nonlinear dimension reduction. See, e.g., [SZ09a, SZ09b, SZ07, SY06, PS03, CS02].

Standard PCA is effective at identifying linear subspaces carrying the greatest variance in a data set. However, this method is not able to detect nonlinear submanifolds. A popular technique to tackle the latter case is kernel PCA. It first maps data into a higher dimensional space $\mathcal{H}\left(K\right)$, and performs PCA there. Here $\mathcal{H}\left(K\right)$ is the reproducing kernel Hilbert space (RKHS) associated with a given positive definite kernel $K$. (Recall that an RKHS is a Hilbert space $\mathcal{H}$ of functions defined on a set $X\left(\neq\emptyset\right)$, such that the evaluation functionals $\mathcal{H}\ni f\rightarrow f\left(x\right)\in\mathbb{C}$ are bounded for all $x\in X$. See also 2.2 below.) The mapping in this context presumably sends a nonlinear submanifold in the input space to a linear subspace in $\mathcal{H}\left(K\right)$. For example, in classification problems, a kernel is usually chosen so that the mapped data can be separated by a linear decision boundary in $\mathcal{H}\left(K\right)$ (see 2.1).

Remark 2.1*.*

It would be intriguing to compare Smale’s Dimension Reduction algorithm from [SZ09a] with ours. The two approaches are along a different lines of development.

The approach in Belkin’s paper [BN03] is popular in current Machine Learning research. Both our results and those of Belkin et al aim for dimension reduction algorithms. Other methods exist, which constitute variants of KPCA, but with different choices of kernels, and with the Laplacian eigenmap (LE) as one of them. (See also [CWG19, VVQCR*+*19, TF19, SGS*+*19].) For recent developments on graph Laplacians, and Perron-Frobenius eigenfunctions as principal components, we refer to e.g., [BJ02].

As applications of the theory presented in the last section, 2, we include a discussion of new non-linear, and real-life data/simulated examples. Moreover, inside the sections, we offer explanations for our use of RKHS-theory, as applied in new KPCA based Gaussian process simulations. In particular, we outline how they serve to produce optimal choices for detection of maximal variance.

The optimization issues in the present section refer to details in subsections 2.3 and 2.4; see especially 2.9, 2.10 and 2.11, all adapted to the PCA issues at hand. Our notation for solutions to optimization problems is “argmax” (standard terminology in the optimization literature.) The role of kernels enters via associated choices of* feature maps* $\Phi$, see (2.45), (2.46), (2.47), (2.48), (2.49), and subsection 2.5. Furthermore, our choices of frames, global frame-operators, and their adjoints (2.14), then allow a* duality approach* which in turn helps us gain insight into algorithms for the optimizers.

By a feature we mean an individual measurable property or characteristic. One chooses informative, discriminating and independent features to be used in algorithmic constructions in machine learning, in pattern recognition, and regression. Vector spaces associated with particular features are called the feature spaces. In KPCA, the choices of feature spaces come about from specified kernels (assumed positive definite, p.d.), see 2.2 below. If $K$ is a fixed p.d. kernel defined on general non-linear sets of data, a natural choice of feature space is then in the form of a Hilbert space is called the reproducing kernel Hilbert space (RKHS), written $\mathscr{H}\left(K\right)$, see 2.3 below.

Definition 2.2**.**

Let $S$ be a set. A positive definite (p.d.) kernel on $S$ is a function $K:S\times S\rightarrow\mathbb{C}$, such that

[TABLE]

for all $\left\{x_{i}\right\}_{i=1}^{N}\subset S$, $\left\{c_{i}\right\}_{i=1}^{N}\subset\mathbb{C}$, and $N\in\mathbb{N}$.

Given a p.d. kernel as in (2.1), there exists a reproducing kernel Hilbert space (RKHS) $\mathcal{H}\left(K\right)$ and a mapping $\Phi:S\rightarrow\mathcal{H}\left(K\right)$ such that

[TABLE]

The function $\Phi$ in (2.2) is called a feature map for the problem.

Moreover, the following reproducing property holds:

[TABLE]

for all $f\in\mathcal{H}\left(K\right)$, and $x\in S$.

While the theory of positive definite kernels and their RKHSs have served for a long time as powerful tools in diverse areas of pure and applied mathematics, the particular applications of them we need here are of a more recent vintage; for background literature, see e.g., [AJ15, BTA04, JT15, JT18a, JT18b]. The theory of frames (non-orthogonal expansions) also plays an essential role in our present approach to PCA, and the list of relevant papers for background includes [HKLW07, HWW05, CCEL15, CCK13].

Feature maps. In the subsequent discussion, reference to kernel will always mean positive definite kernel, see (2.1). For each choice of kernel, defined on data sets (typically non-linear configurations), there will then be associated feature maps. Generally, a feature map $\Phi$ is simply a function which maps sets of data configurations (non-linear) into some choice of feature space, a Hilbert space, say $\mathscr{L}$. Typically the inner product in $\mathscr{L}$ will serve to model correlation, or regression numbers, or other measurements corresponding to other features of interest. Given a kernel $K$, a natural choice of feature space $\mathscr{L}$ is the RKHS defined directly from $K$, see (2.2), but there are other possibilities. In a different choice for $\Phi$, the range of $\Phi$, and the inner product on the right-hand side in (2.2) might be different. The main logic in the use of kernels in machine learning is that it yields representations of learning algorithm for data sets that will then become more amenable to regression analysis, detection, learning, and to classification. Analysis is aided with the use of powerful tools from the theory of linear operators in Hilbert space; see e.g., [HKLW07, JT18b, SZ07]. We stress that the phrase feature map is broad, and that a wide variety of functions $\Phi$ may serve as feature maps. But the main use of them relates to suitable choices of kernels. Support Vector Machines (and other kernel-based methods) make use both implicit and explicit feature maps. This leads to remapping of data and will allow non-linearly separable data sets to get linear representations. The feature space leads to effective separation of data via suitable choices of hyperplanes in higher dimension. But reaching these dimensions might be computationally expensive because bad choices of feature mappings might require many computations.

Remark 2.3*.*

$\mathcal{H}\left(K\right)$ may be chosen as the Hilbert completion of

[TABLE]

with respect to the $\mathcal{H}\left(K\right)$-inner product

[TABLE]

Initially the LHS in formula (2.5) only refers to finite linear combinations. Hence, the vector space (2.4) becomes a pre-Hilbert space. (By pre-Hilbert space we mean an inner product space that is not complete.) The RKHS $\mathcal{H}(K)$ itself then results from the standard Hilbert completion. It is this Hilbert space we will use in our subsequent study of optimization problems, and in our KPCA-dimension reduction. Sections 2.1–2.3 deal with separate issues of kernel-optimization. Before turning to these, however, we will first introduce a setting of Hilbert-Schmidt operators. This will play a crucial role in the formulation of our main result, 2.11 in 2.3.

Recall that a data set $\left\{x_{j}\right\}_{j=1}^{n}$, $x_{j}\in\mathbb{C}^{m}$, may be viewed as an $m\times n$ matrix $X$, where $x_{j}$ is the $j^{th}$ column vector. Here, $m$ is the number of features, and $n$ the number of sample points. Set

[TABLE]

where $\left\|\cdot\right\|_{HS}$ in (2.6) denotes the Hilbert-Schmidt norm.

Remark 2.4*.*

Let $\mathcal{H}$ be a Hilbert space, and let $HS\left(\mathcal{H}\right)$ be the Hilbert-Schmidt operators $\mathcal{H}\xrightarrow{\;X\;}\mathcal{H}$ with inner product

[TABLE]

Then the two Hilbert spaces $HS\left(\mathcal{H}\right)$, and $\mathcal{H}\otimes\overline{\mathcal{H}}$ (tensor-product), are naturally isometrically isomorphic via

[TABLE]

Indeed,

[TABLE]

and the assertion follows from isometric extension of (2.8).

2.1. Application to Optimization

One of the more recent applications of kernels and the associated reproducing kernel Hilbert spaces (RKHS) is to optimization, also called kernel-optimization. See [YLTL18, LLL11]. In the context of machine learning, it refers to training-data and feature spaces. In the context of numerical analysis, a popular version of the method is used to produce splines from sample points; and to create best spline-fits. In statistics, there are analogous optimization problems going by the names “least-square fitting,” and “maximum-likelihood” estimation. In the latter instance, the object to be determined is a suitable probability distribution which makes “most likely” the occurrence of some data which arises from experiments, or from testing.

What these methods have in common is a minimization (or a max problem) involving a “quadratic” expression $Q$ with two terms. The first in $Q$ measures a suitable $L^{2}\left(\mu\right)$-square applied to a difference of a measurement and a “best fit.” The latter will then to be chosen from a number of suitable reproducing kernel Hilbert spaces (RKHS). The choice of kernel and RKHS will serve to select desirable features. So we will minimize a quantity $Q$ which is the sum of two terms as follows: (i) a $L^{2}$-square applied to a difference, and (ii) a penalty term which is a RKHS norm-squared. (See eq. (2.10).) In the application to determination of splines, the penalty term may be a suitable Sobolev normed-square; i.e., $L^{2}$ norm-squared applied to a chosen number of derivatives. Hence non-differentiable choices will be “penalized.”

In all of the cases, discussed above, there will be a good choice of (i) and (ii), and we show that there is then an explicit formula for the optimal solution; see eq (2.14) in 2.5 below.

Let $X$ be a set, and let $K:X\times X\longrightarrow\mathbb{C}$ be a positive definite (p.d.) kernel. Let $\mathcal{H}\left(K\right)$ be the corresponding reproducing kernel Hilbert space (RKHS). Let $\mathscr{B}$ be a sigma-algebra of subsets of $X$, and let $\mu$ be a positive measure on the corresponding measure space $\left(X,\mathscr{B}\right)$. We assume that $\mu$ is sigma-finite. We shall further assume that the associated operator $T$ given by

[TABLE]

is densely defined and closable.

Fix $\beta>0$, and $\psi\in L^{2}\left(\mu\right)$, and set

[TABLE]

defined for $f\in\mathcal{H}\left(K\right)$ , or in the dense subspace $dom\left(T\right)$ where $T$ is the operator in (2.9). Let

[TABLE]

be the corresponding adjoint operator, i.e.,

[TABLE]

Our present RKHS/$L^{2}\left(\mu\right)$ framework is close to that of [SS16]. But we have included the results we need in our present framework, and for our purpose.

Theorem 2.5**.**

Let $K$, $\mu$, $\psi$, $\beta$ be as specified above; then the optimization problem

[TABLE]

has a unique solution $F$ in $\mathcal{H}\left(K\right)$, it is

[TABLE]

where the operator $T$ and $T^{*}$ are as specified in (2.9)-(2.12).

Proof.

We fix $F$, and assign $f_{\varepsilon}:=F+\varepsilon h$ where $h$ varies in the dense domain $dom\left(T\right)$ from (2.9). For the derivative \frac{d}{d\varepsilon}\big{|}_{\varepsilon=0} we then have:

[TABLE]

for all $h$ in a dense subspace in $\mathcal{H}\left(K\right)$. The desired conclusion follows. Note that convexity of the function $Q_{\psi,\beta}$ makes these conditions sufficient as well. ∎

As already mentioned, the general conclusion in the present 2.5 connects to both many uses of kernel techniques in classical analysis, as well as to new and exciting applications. In the second group, we mention the following approach to Neural Networks [YO21], leading in turn to encoding of deep regularization into training of inner layers which make up specific Neural Network (NN) constructions, with the use of kernel flows. For details, we refer to [YO21] and the cited literature for details. But here we only stress that, as a step in particular NN-designs, one makes use of an recursive iteration of substitution into a fixed kernel (the kernel resulting from the substitution iterations is called a warped kernel); and then, at each step, there will be an RKHS-optimization of a term that arises as a special case of (2.13), see [YO21, eq (1.2)].

Least-square Optimization

To help readers put our present results on feature space and RKHS penalty-term into context, we now specialize the optimization formula from 2.5 to the problem of minimize a “quadratic” quantity $Q$. It is still the sum of two individual terms: (i) a $L^{2}$-square applied to a difference, and (ii) a penalty term which is the RKHS norm-squared. But the least-square term in (i) will simply be a sum of a finite number of squares of differences; hence “least-squares.” As an application, we then get an easy formula (2.6) for the optimal solution.

Let $K$ be a positive definite kernel on $X\times X$ where $X$ is an arbitrary set, and let $\mathcal{H}\left(K\right)$ be the corresponding reproducing kernel Hilbert space (RKHS). Let $m\in\mathbb{N}$, and consider sample points:

$\left\{t_{j}\right\}_{j=1}^{m}$ as a finite subset in $X$, and

$\left\{y_{i}\right\}_{i=1}^{m}$ as a finite subset in $\mathbb{R}$, or equivalently, a point in $\mathbb{R}^{m}$.

Fix $\beta>0$, and consider $Q=Q_{\left(\beta,t,y\right)}$, defined by

[TABLE]

We introduce the associated dual pair of operators as follows:

[TABLE]

where

[TABLE]

for all $\vec{y}=\left(y_{i}\right)\in\mathbb{R}^{m}$.

Note that the duality then takes the following form:

[TABLE]

consistent with (2.12).

Applying 2.5 to the counting measure

[TABLE]

for the set of sample points $\left\{t_{i}\right\}_{i=1}^{m}$, we get the two formulas:

[TABLE]

where $K_{m}$ denotes the $m\times m$ matrix

[TABLE]

Theorem 2.6**.**

Let $K$, $X$, $\left\{t_{i}\right\}_{i=1}^{m}$, and $\left\{y_{i}\right\}_{i=1}^{m}$ be as above, and let $K_{m}$ be the induced sample matrix (2.22). Fix $\beta>0$; consider the optimization problem with

[TABLE]

Then the unique solution to (2.23) is given by

[TABLE]

i.e., $F=\arg\min Q$ on $\mathcal{H}\left(K\right)$.

Proof.

From 2.5, we get that the unique solution $F\in\mathcal{H}\left(K\right)$ is given by:

[TABLE]

and by (2.20)-(2.21), we further get

[TABLE]

where the dot $\cdot\phantom{}$ refers to a free variable in $X$. An evaluation of (2.25) on the sample points yields:

[TABLE]

where $\vec{F}:=\left(F\left(t_{i}\right)\right)_{i=1}^{m}$, and $\vec{y}=\left(y_{i}\right)_{i=1}^{m}$. Hence

[TABLE]

Now substitute (2.27) into (2.26), and the desired conclusion in the theorem follows. We used the matrix identity

[TABLE]

∎

2.2. The Case of Gaussian Fields

For a number of applications, it will be convenient to consider general stochastic processes $\left(X_{s}\right)$ indexed by $s\in S$, where $S$ is merely a set; so not a priori equipped with any additional structure. Consideration of stochastic processes will always assume some fixed probability space, $\left(\Omega,\mathscr{A},\mathbb{P}\right)$ where $\Omega$ is a set of sample points; $\mathscr{A}$ is a $\sigma$-algebra of events, fixed at the outset; and $\mathbb{P}$ is a probability measure defined on $\mathscr{A}$. A given process $\left(X_{s}\right)_{s\in S}$ is then said to be Gaussian and centered iff (Def.) for all choice of finite subsets of $S$ ($s_{1},s_{2},\cdots,s_{N}$), then the system of random variables $\left\{X_{s_{i}}\right\}_{i=1}^{N}$ is jointly Gaussian, i.e., the joint distribution of $\left\{X_{s_{i}}\right\}_{i=1}^{N}$ on $\mathbb{R}^{N}$ is the Gaussian $g_{N}\left(x_{1},\cdots,x_{N}\right)$ which has mean zero, and covariance matrix

[TABLE]

so for $x=\left(x_{1},\cdots,x_{N}\right)\in\mathbb{R}^{N}$,

[TABLE]

If $A_{N}\subset\mathbb{R}^{N}$ is a Borel set, then

[TABLE]

holds. Note we consider the joint distributions for all finite subsets of $S$.

Let $S$ be a set, and let $\left\{X_{s}\right\}_{s\in S}$ be a Gaussian process with $\mathbb{E}\left(X_{s}\right)=0$, $\forall s\in S$; and with

[TABLE]

as its covariance kernel. Finally, let $\mathcal{H}\left(K\right)$ be the corresponding RKHS.

Then the following general results hold (see e.g., [AJL11, AJ12, AJS14, AJ15, JT15, JT16a, JT16b, AJL17, JT18a, JT18b]):

1. (i)

Every positive definite kernel $S\times S\xrightarrow{\;K\;}\mathbb{C}$ arises as in (2.31) from some Gaussian process $\left\{X_{s}\right\}_{s\in S}$. 2. (ii)

Assume $\mathcal{H}\left(K\right)$ is separable; then we have a representation $\left\{f_{n}\right\}_{n\in\mathbb{N}}$ for a system of functions $f_{n}:S\rightarrow\mathbb{C}$, $n\in\mathbb{N}$,

[TABLE]

absolutely convergent on $S\times S$. 3. (iii)

A system $\left\{f_{n}\right\}_{n\in\mathbb{N}}$ satisfies (ii) if and only if it forms a Parseval frame in $\mathcal{H}\left(K\right)$. 4. (iv)

Given (2.32), then, for every sequence of independent identically distributed (i.i.d.) Gaussian system \left\{Z_{n}\right\}_{n\in\mathbb{N}},$$Z_{n}\sim N\left(0,1\right), i.e., each $Z_{n}$ is a standard Gaussian random variable, $\mathbb{E}\left(Z_{n}\right)=0$, $\mathbb{E}\left(Z_{n}Z_{m}\right)=\delta_{n,m}$; the representation

[TABLE]

is valid in $L^{2}$ of the underlying probability space, and $\mathbb{E}\left(X_{s}X_{t}\right)=K\left(s,t\right)$ on $S\times S$.

Note 2.7*.*

When a fixed Gaussian process $(X_{s})$ is given, then the associated decomposition (2.33) is called a Karhunen-Loève (KL) transform for $X_{s}$. The conclusion from (i)–(iv) above is that there is a direct connection between the two KL transforms the relatively better known KL-transforms for positive definite kernels ([JS07] Theorem 4.15).

The object in principal component analysis (PCA) is to find optimal representations; and to select from them the “leading terms”, the principal components.

Remark 2.8*.*

The present general kernel framework (RKHSs and Gaussian processes) encompasses the special case we outlined in [JS07] Example 3.1. By way of comparison, note that the particular positive definite kernel in the latter example is only a special case of the present ones, see (2.32) and (2.33). These types of kernels are often referred to as the case of Mercer kernels; see also [SY06, SZ07, SZ09a, SZ09b]. A Mercer kernel is continuous, and it defines a trace class operator, as illustrated in the example. This latter feature in turn leads to a well defined “top part of the spectrum.” And this then allows us to select the principal components; i.e., the maximally correlated variables. We shall show, in 2.11 below, that there is an alternative approach to principal components which applies to the general class of positive definite kernels, and so goes far beyond the case of Mercer kernels.

2.3. Optimization and Frames

Fix a p.d. kernel $K$ on $S:=\mathbb{C}^{m}$, i.e., a functional $K:\mathbb{C}^{m}\times\mathbb{C}^{m}\rightarrow\mathbb{C}$ satisfying (2.1); and let $\mathcal{H}\left(K\right)$ be the associated RKHS. In PCA, one solves the quadratic optimization problem:

[TABLE]

where $Q$ is a selfadjoint projection.

Kernel PCA, by contrast, solves a similar problem in $\mathcal{H}\left(K\right)$:

[TABLE]

where

[TABLE]

It is understood that $\left\|\cdot\right\|_{HS}$ as in (2.35) refers to the Hilbert-Schmidt class in $B\left(\mathcal{H}\left(K\right)\right)$.

A Finite Frame

Let $X$, $K$, and $\Phi$ be as above. Then $\left(\Phi\left(x_{j}\right)\right)_{j=1}^{n}$ is a finite frame whose span $\mathcal{H}_{\Phi}$ is a closed subspace in $\mathcal{H}\left(K\right)$.

Set $L:\mathcal{H}\left(K\right)\rightarrow\mathbb{C}^{n}$ by

[TABLE]

where $\delta_{j}$ is the standard ONB in $\mathbb{C}^{m}$. The adjoint $L^{*}:\mathbb{C}^{n}\rightarrow\mathcal{H}\left(K\right)$ is given by

[TABLE]

It follows that

[TABLE]

Lemma 2.9**.**

Let $L^{*}L$ be as in (2.39), and let

[TABLE]

be the corresponding spectral representation, with $\lambda_{1}^{2}\geq\lambda_{2}^{2}\geq\cdots$. Then

[TABLE]

Equivalently, the best rank-1 approximation to $G$ is

[TABLE]

Remark 2.10*.*

Note that the conclusion of the lemma yields a solution the optimization problem we introduced above. Indeed, in the statement of the lemma (see (2.40)) we use the standard notation argmax for the data which realizes a particular optimization. In the present case, we are maximizing a certain quadratic expression over the unit-ball in the Hilbert-Schmidt operators. The Hilbert-Schmidt norm is designated with the subscript HS. Part of the conclusion of the lemma asserts that the maximum, as specified in (2.40), is attained for a definite rank-one operator.

Proof of 2.9.

Note that

[TABLE]

Let $w$ be a unit vector in $\mathcal{H}\left(K\right)$, and set $Q_{1}=\left|w\left\rangle\right\langle w\right|$; then

[TABLE]

Since $\sum_{j}\left|\left\langle w,\phi_{j}\right\rangle\right|^{2}=\left\|w\right\|^{2}=1$, the r.h.s. of (2.41) is a convex combination of $\lambda_{j}^{2}$’s; therefore,

[TABLE]

and equality holds if and only if $\left|\left\langle w,\phi_{1}\right\rangle\right|=1$, and $\left\langle w,\phi_{j}\right\rangle=0$, for $j>1$.

Note that, given (2.41), the conclusion of the lemma follows directly from comparison of positive series. ∎

2.9 can be applied inductively which yields the best rank-1 approximation at each iteration. In fact, the result holds more generally; see 2.11 below.

Our present RKHS/$L^{2}(\mu)$ framework is close to that of [SS16]. But we have included the results we need in our present framework, and for our purpose.

Theorem 2.11**.**

(a) Let $T:\mathcal{H}\rightarrow\mathcal{H}$ be a Hilbert Schmidt operator, and let

[TABLE]

with $\lambda_{1}^{2}\geq\lambda_{2}^{2}\geq\cdots$, and $\lambda_{j}\rightarrow 0$ as $j\rightarrow\infty$. Then

[TABLE]

Note that $\left\|Q_{n}T\right\|_{HS}^{2}=tr\left(Q_{n}G\right)$, where $G:=TT^{*}$. See [JS07]. (b) Recall the assumption that $T$ is Hilbert Schmidt is equivalent to $T^{*}T$ trace class. We then apply the Spectral Theorem to $T^{*}T$ yielding an orthonormal basis of eigenvectors for $T^{*}T$, and a convergent representation for the infinite sum in (2.42).

Proof.

Let $Q_{n}=\sum_{j=1}^{n}P_{j}$, where each $P_{j}$ is rank-1, and $P_{i}P_{j}=0$, for $i\neq j$. Then,

[TABLE]

Let $\mathcal{L}$ be the corresponding Lagrangian, i.e.,

[TABLE]

It follows that

[TABLE]

Hence $P_{k}$ is a spectral projection of $TT^{*}$. The conclusion of the theorem follows from this.

Note that one may use the usual variational directional derivative argument applied to a sesquilinear form $\left(x,x\right)\mapsto\left\langle x,Tx\right\rangle_{\mathscr{H}}$ in any Hilbert space $\mathscr{H}$, and obtain

[TABLE]

In particular, if $T$ is self-adjoint, then

[TABLE]

Especially, the Hilbert space in (2.43) is the space of Hilbert-Schmidt operators. (This argument is a key step in the proof of the spectral theorem, i.e., taking the variational directional derivative of a bounded bilinear form in Hilbert space.)

The precise meaning of (and the justification for) the generalized gradient assertion (2.43) is included in 2.12 below. See also the following citations of sources for the underlying operator theory. ∎

The variational arguments in the proof of 2.11, used in the present setting of Hilbert space, and the variety of projections, can be justified with standard tools from operator theory, including the spectral theorem; see e.g., [JT17, ch.2], [JT21, ch.1], [Con90, ch.2] and also [AK06, BJ02, HKLW07, SZ09a]. We have included more details below:

Lemma 2.12**.**

Suppose $T$ is a Hilbert-Schmidt (HS) operator in a separable Hilbert space $\mathscr{H}$. For the directional derivative $\nabla_{T}:=\frac{\partial}{\partial T}$, we have

[TABLE]

(Here, $\left|x\left\rangle\right\langle x\right|$ is the rank-1 operator $\mathscr{H}\ni u\mapsto\left\langle x,u\right\rangle x$ using Dirac’s notation.)

Proof.

Let $Q=\left|x\left\rangle\right\langle x\right|$. Note that $f_{Q}\left(T\right):=\left\langle x,Tx\right\rangle=\text{trace}\left(TQ\right)$. In fact, if $\left\{e_{n}\right\}_{1}^{\infty}$ is an orthonormal basis in $\mathscr{H}$, then

[TABLE]

But

[TABLE]

where $\left\langle\cdot,\cdot\right\rangle_{HS}$ denotes the Hilbert-Schmidt inner product.

Now we use a general fact about inner products:

[TABLE]

where $\mathscr{H}S=$ the Hilbert-Schmidt operators, viewed as a Hilbert space, relative to its trace-inner product (2.44). That is, the Schwarz inequality is applied to the $\left\langle\cdot,\cdot\right\rangle_{HS}$-inner product. One obtains maximum when $T\sim Q$. In particular, if $Q=\left|x\left\rangle\right\langle x\right|$, then $Q^{2}=\left\|x\right\|^{2}Q$ and $\text{trace}\left(Q^{2}\right)=\left\|x\right\|_{\mathscr{H}}^{4}$.

Therefore, it follows from $f_{Q}\left(T\right)=\left\langle Q,T\right\rangle_{HS}$ that

[TABLE]

which is identified with $Q=\left|x\left\rangle\right\langle x\right|$, since

[TABLE]

∎

Remark 2.13*.*

Our use of frame analysis serves as the technical tools for “error estimates.” More specifically, frames are designed such that, for a given problem involving finite-dimensional PCA computations, the frame estimates are designed to allow estimation of the corresponding “error terms.” (See [JS07].)

For our present considerations of optimization questions, it may be of interest to consider a comparison between the following two settings, one general, and the other special: On the one hand, there is a variety of (i) general calculus of variation issues in infinite-dimensional contexts; and on the other, (ii) particular optimization questions in the restricted framework of specific choices of pairs of Hilbert spaces (still infinite-dimensional). We note that the context for (i) is wider than that of (ii). On the other hand (as outlined below), for the specialized framework of (ii) considered here, the use of standard Hilbert space geometry, and a little operator theory, offer much simplification.

Detail: Our present focus is (ii). By contrast, in consideration of (i) one must address separately such technical issues as (a) functional-derivatives, and (b) existence of extremal (optimal) solutions. By contrast, the smaller class of optimization questions from (ii) lend themselves directly to computation in Hilbert space. Here we note that application of Hilbert space geometry will then simplify the two issues (a) and (b) in (i). The point we make for our present optimization question (i.e., a specific choice of a Hilbert space context (ii)), is that a “natural” choice of a pair of Hilbert spaces will be a combination of the two: an appropriate RKHS, and a compatible $L^{2}\left(\mu\right)$-Hilbert space. This then allows for relatively simple solutions.

2.4. The Dual Problem

Fix a data set $X=\left(x_{j}\right)_{j=1}^{n}$, $x_{j}\in\mathbb{C}^{m}$. Let $\Phi:X\rightarrow\mathcal{H}\left(K\right)$ be the feature map in (2.36), i.e.,

[TABLE]

Let $L$, $L^{*}$ be the analysis and synthesis operators from (2.37)-(2.38), and

[TABLE]

be the frame operator in (2.39). In particular,

[TABLE]

In view of [JS07] and 2.11, the KL basis for $L^{*}L$ contains the principal directions carrying the greatest variance in $\Phi\left(X\right)$. In applications, it is more convenient to first find the KL basis of $LL^{*}$ instead, where

[TABLE]

as an $n\times n$ matrix in $\mathbb{C}^{n}$; see (2.2). (By general theory, if $A:\mathcal{H}\rightarrow\mathcal{H}$ is a linear operator in a Hilbert space with dense domain, then $\sigma\left(A^{*}A\right)\backslash\left\{0\right\}=\sigma\left(AA^{*}\right)\backslash\left\{0\right\}$.)

Theorem 2.14**.**

Set $A:\mathbb{C}^{n}\rightarrow\mathcal{H}\left(K\right)$ by

[TABLE]

and extend linearly, where $\left(\delta_{j}\right)_{j=1}^{n}$ denotes the standard basis in $\mathbb{C}^{n}$. Then the adjoint operator $A^{*}:\mathcal{H}\left(K\right)\rightarrow\mathbb{C}^{n}$ is

[TABLE]

That is, $A=L^{*}$ and $A^{*}=L$.

Proof.

Let $v\in\mathbb{C}^{n}$, and $h\in\mathcal{H}$, then

[TABLE]

and the assertions follows. ∎

Hence $LL^{*}:\mathbb{C}^{n}\rightarrow\mathbb{C}^{n}$ is the Gramian matrix in $\mathbb{C}^{n}$ given by

[TABLE]

see (2.47).

By the singular value decomposition, $L^{*}=WDU^{*}$, so that

[TABLE]

where $D=\text{diag}\left(\lambda_{j}\right)$ consists of the non-negative eigenvalues of $\sqrt{LL^{*}}$. Therefore,

[TABLE]

Note that $W=\left(w_{j}\right)$ is the KL basis that diagonalizes $L^{*}L$ as in (2.46), i.e.,

[TABLE]

It also follows from (2.51)–(2.52), that

[TABLE]

Remark 2.15*.*

In the above discussion, $\Phi\left(X\right)$ may be centered by removing its mean. Specifically, let

[TABLE]

be the projection onto $\text{span}\left\{\mathbbm{1}\right\}^{\perp}$, where $\mathbbm{1}$ denotes the constant vector $\begin{bmatrix}1&\cdots&1\end{bmatrix}$. Then

[TABLE]

and so

[TABLE]

The effect of $J$ in (2.56) is to exclude the eigenspace of the Gramian $\left(K\left(x_{i},x_{j}\right)\right)_{i,j=1}^{n}$ spanned by the constant eigenvector.

In what follows, we shall always assume $\Phi\left(X\right)$ is centered as in (2.55)-(2.56).

2.5. Feature Selection

Feature selection, also called variable selection, or attribute selection, is a procedure for automatic selection of those attributes in data sets which are most relevant to particular predictive modeling problems. Which features should one use in designs of predictive models? This is a difficult question that requires detailed knowledge of the problem at hand. The aim is algorithmic designs which automatically select those features from prescribed data, which are most useful, or most relevant, for the particular problem. The process is called feature selection. A central premise of feature selection is that the input data will contain features that are either redundant or irrelevant, and can therefore be removed. The use of sample correlations in the process is based in turn on the following principle: A particular relevant feature might be redundant, in the presence of some other relevant feature, with which it is strongly correlated.

Our present purpose is not a systematic treatment of feature selection, but merely to identify how our present tools suggest recursive algorithms in the general area. With this in mind we now consider the following setup:

Let $x\in\mathbb{C}^{m}$ be a test example. The image $\Phi\left(x\right)$ under the feature map can be projected onto the principal directions in $\mathcal{H}\left(K\right)$, via

[TABLE]

The mapping $x\mapsto WW^{*}\Phi\left(x\right)$ is in general nonlinear. See Examples 2.17 and 2.18, and Figures 2.1 and 2.2 below. We now turn to our main result which is the following corollary.

Corollary 2.16**.**

Let $L^{*}=WDU^{*}$ be as above, assuming $D$ is full rank. For all $x\in\mathbb{C}^{m}$, the coefficients of the projection $WW^{*}\Phi\left(x\right)$ are

[TABLE]

Proof.

By (2.54), $W^{*}=D^{-1}U^{*}L$, so that

[TABLE]

∎

Example 2.17** (Spectral Clustering, see 2.1).**

This is included as an instructive example. The data set $X$ has two classes: Class “0” consists of $867$ points uniformly distributed in $\{x\in\mathbb{R}^{3}\mathrel{;}0.6<\left\|x\right\|^{2}<1\}$; class “1” consists of 126 points in the open ball $\{x\in\mathbb{R}^{3}\mathrel{;}\left\|x\right\|^{2}<0.2\}$. Hence $X$ has dimension $3\times 993$, and each column $x_{j}$ of $X$ corresponds to a sample point.

We choose the Gaussian kernel $K\left(x,y\right)=e^{-\left\|x-y\right\|^{2}/\sigma}$ with $\sigma=0.05$, so that $x_{j}$ is embedded into the associated RKHS by

[TABLE]

By projecting $\Phi\left(x_{j}\right)$ onto the first two principal components, each sample point $x_{j}\in\mathbb{R}^{3}$ has a 2D representation via the mapping

[TABLE]

See the general formula in 2.16; note that $n=2$ in the present example.

As shown in 2.1b, the two clusters are linearly separable in $\mathcal{H}\left(K\right)$.

Example 2.18** (Dimension Reduction, see 2.2).**

Let $X$ be a collection of 100 grayscale $256\times 256$ images of an ellipse, rotated successively by $\pi/100$. 2.2a shows 6 sample images corresponding to different rotation angles. The images are unrolled as column vectors, thus $X$ has dimension $65536\times 100$.

This data set may be viewed as 1D submanifold embedded in $\mathbb{R}^{65536}$, i.e., it has only one degree of freedom, the rotation angle. For dimension reduction, KPCA will ideally extract this information, and each image is then represented by a single projection coefficient. We choose the Gaussian kernel with $\sigma=300$. In 2.2b, there are 4 subplots consisting of the projections onto the first, second, third, and fourth principal directions. The rotation angle is encoded in e.g. PC 1. As a consequence, the dimension of the data set is reduced from $65536\times 100$ to $1\times 100$.

In particular, the projection coefficients onto the $k^{th}$ principal direction (see 2.2b, PC1 – PC4) are proportional to the $k^{th}$ eigenvector of the centered Gramian $G:=JK\left(x_{i},x_{j}\right)J$; see (2.56) and 2.16.

2.6. Dynamic PCA

In some applications, the tool of Dynamic PCA (DPCA), see e.g., [HAI*+*07] serves as a useful tool for high-dimensional and time-dependent data. The idea is that in favorable cases, with the use of DPCA, one can arrange that the input matrices can be augmented by addition of time-lagged values of the variables under consideration. In summary, the method is based on the behavior of the eigenvalues of the lagged autocorrelation, and partial autocorrelation, matrices.

Our present 2.11 states that when a system of PCA eigenvalues $\left\{\lambda_{i}\right\}$ is given, then we have a algorithmic solution to the corresponding optimization question (see (2.34)-(2.35)). We now turn to a formula for generating PCA features as a limit of a certain iteration of operators. The family of operators $T$ discussed below is a generalization of the operators from 2.3 above.

When principal component analysis (PCA) is used for problems involving stochastic (statistical) processes, for example in monitoring applications, it often relies on an implicit, but unrealistic, assumption that data involved are time independent. This however is evidently unrealistic, as work with industrial data shows: In a host of diverse applications, one is faced with serial correlations. Hence, the literature over recent years (see e.g., [AHR*+*20, KOL20, Sha20, GS20, Sha19]) has witnessed a host of new approaches to PCA, going by what is now called dynamic PCA (DPCA) methods. DPCA-models involve a variety of stochastic and dynamical systems-tools, each one serving as a remedy for plain-vanilla PCA applied to data as they appear in: (i) in extremely high-dimensions, (ii) in time-dependent data-sets, and (iii) in digital image processing. A main feature of DPCA is time-feedback; — hence in DPCA, the input-matrix is dynamically augmented, thus taking the form of a transfer operator;— in the simplest models, it operates by dynamic additions of time-lagged values of the variables.

Proposition 2.19**.**

Suppose $T:\mathcal{H}\rightarrow\mathcal{H}$ is compact, and $T^{*}T$ has simple spectrum, i.e.,

[TABLE]

with $\lambda_{1}^{2}>\lambda_{2}^{2}>\cdots>0$; $\lambda_{j}^{2}\rightarrow 0$ as $j\rightarrow\infty$. Setting $S=\sqrt{T^{*}T}$, then,

[TABLE]

Moreover, given $\lambda_{1}$, then

[TABLE]

where convergence in (2.59) is w.r.t. the norm topology of $B\left(\mathcal{H}\right)$.

Proof.

Note that

[TABLE]

Now, given $\lambda_{1}$, we have

[TABLE]

∎

Corollary 2.20**.**

The system of PCA eigenvalues $\left\{\lambda_{j}\right\}$ can be obtained inductively as follows: Let $S=\sqrt{T^{*}T}$ as before, and set $Q_{k}=\sum_{j=1}^{k}\left|u_{j}\left\rangle\right\langle u_{j}\right|$. Then

[TABLE]

and

[TABLE]

Proof.

One checks that

[TABLE]

and so the assertion follows from 2.19. ∎

Example 2.21**.**

Consider the covariance function of standard Brownian motion $B_{t}$, $t\in[0,\infty)$, i.e., a Gaussian process $\left\{B_{t}\right\}$ with mean zero and covariance function

[TABLE]

Let $F_{N}=\left\{x_{1},x_{2},\ldots,x_{N}\right\}$ be a finite subsets of $V$, such that

[TABLE]

and let

[TABLE]

Lemma 2.22**.**

Let $K_{N}$ be as in (2.61). Then

1. (i)

The determinant of $K_{N}$ is given by

[TABLE] 2. (ii)

$K_{N}$* assumes the LU decomposition*

[TABLE]

where

[TABLE]

Proof.

For details, see e.g., [JT15]. ∎

Example 2.23**.**

Fix $N$, then the top eigenvalue of $K:=K_{N}$ can be extracted by the method from 2.19. For instance, if $N=3$ and $F_{3}=\left\{1,2,3\right\}$, we have

[TABLE]

Let $e_{1}=\left[1,0,0\right]$, then

[TABLE]

Standard numerical algorithm returns

[TABLE]

Acknowledgement*.*

The co-authors thank the following colleagues for helpful and enlightening discussions: Professors Daniel Alpay, Sergii Bezuglyi, Ilwoo Cho, Wayne Polyzou, David Stewart, Eric S. Weber, and members in the Math Physics seminar at The University of Iowa. The present work was started with discussions at the NSF-CBMS Conference, Harmonic Analysis: Smooth and Non-Smooth, held at Iowa State University, June 4-8, 2018. Jorgensen was the main speaker. We are grateful to the organizers, especially to Prof Eric Weber, and to the NSF for financial support. Sooran Kang was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (#NRF-2017R1D1A1B03034697 and #NRF-2020R1F1A1A01076072).