Reducing the dimensionality of data using tempered distributions

TL;DR

This paper introduces a novel approach to unsupervised dimension reduction by reformulating it with tempered distributions, establishing a connection with sufficient dimension reduction, and proposing algorithms for optimization and practical implementation.

Contribution

It presents a new formulation of dimension reduction using tempered distributions, linking it to SDR, and develops algorithms including an efficient approximate method for practical use.

Findings

Algorithms successfully reduce data dimensionality in synthetic and real datasets.

The approach effectively approximates empirical distributions with low-dimensional support.

The method demonstrates comparable or improved performance over existing techniques.

Abstract

We reformulate unsupervised dimension reduction problem (UDR) in the language of tempered distributions, i.e. as a problem of approximating an empirical probability density function by another tempered distribution, supported in a $k$-dimensional subspace. We show that this task is connected with another classical problem of data science -- the sufficient dimension reduction problem (SDR). In fact, an algorithm for the first problem induces an algorithm for the second and vice versa. In order to reduce an optimization problem over distributions to an optimization problem over ordinary functions we introduce a nonnegative penalty function that ``forces'' the support of the model distribution to be $k$-dimensional. Then we present an algorithm for the minimization of the penalized objective, based on the infinite-dimensional low-rank optimization, which we call the alternating scheme. Also, we design an efficient approximate algorithm for a special case of the problem, where the distance between the empirical distribution and the model distribution is measured by Maximum Mean Discrepancy defined by a Mercer kernel of a certain type. We test our methods on four examples (three UDR and one SDR) using synthetic data and standard datasets.

Tables

Table 1. Table 1: The cross-validated accuracies/R 2 of KNN on 2 or 3-dimensional input representations.

	PCA		SIR		KDR		AS 1
Dimension $k$	2	3	2	3	2	3	2	3
Heart (acc)	79.80	79.46	82.49	81.82	86.33	88.77	81.48	83.50
Breast (acc)	93.46	93.65	97.30	96.73	93.13	95.95	97.88	97.69
Ionosphere (acc)	80.29	86.57	89.14	89.43	83.43	86.29	88.29	90.57
Diabetes (R2)	25.34	28.72	43.47	43.61	41.82	44.30	43.07	44.48
Boston (R2)	56.42	67.12	76.03	74.29	77.88	79.97	73.21	77.88
Wine (R2)	93.91	94.12	98.68	99.24	98.30	96.02	97.10	96.93

Table 2. Table 2: Original images (the first row) and their projections to 1-dimensional subspaces computed by PCA (the second row) and computed by Gauss (the third row).

Table 3. Table 3: Original image, projected image and the difference between them (Gauss).

Table 4. Table 4: Original image, its projection, and their grayscale difference (Gauss).

Table 5. Table 5: Computed backgrounds are almost identical, though noise is more often in PCA’s output.

	Input images
	A background computed by PCA
	A background computed by Gauss

Table 6. Table 6: Measures of consistency with the ground truth background image for the SBMnet dataset.

Method	AGE	pEPs	pCEPS	MSSSIM	PSNR	CQM
MSCL (SOTA) [46]	5.9547	0.0524	0.0171	0.9410	30.8952	31.7049
BRTF [47]	9.5385	0.1140	0.0876	0.9621	28.4655	29.3246
GoDec [48]	11.5934	0.1584	0.0974	0.8854	24.9954	25.9955
PCA	9.3774	0.0904	0.0522	0.9027	26.1549	27.5052
Kurtosis	9.3509	0.0936	0.0544	0.9032	26.1475	27.468
Gauss ($a=0.2$)	9.251	0.09	0.0521	0.9025	26.1649	27.464
Gauss ($a=1.0$)	8.8679	0.0876	0.05	0.9049	26.7391	28.0609
Gauss ($a=5.0$)	8.85	0.0876	0.0497	0.9045	26.7254	28.0586
Gauss ($a=25.0$)	8.8781	0.0886	0.0509	0.9065	27.0038	28.3913
Laplace ($a=0.2$)	9.0745	0.089	0.0511	0.9032	26.3269	27.619
Laplace ($a=1.0$)	8.9428	0.088	0.0505	0.904	26.5121	27.819
Laplace ($a=5.0$)	8.8424	0.0873	0.0498	0.906	26.8228	28.1728
Poisson ($a=0.04$)	9.2481	0.0905	0.0523	0.9021	26.173	27.4645
Poisson ($a=1.0$)	9.2483	0.0906	0.0523	0.9022	26.173	27.4644
Poisson ($a=25.0$)	9.2481	0.0906	0.0523	0.9021	26.173	27.4646

Equations

$$p_{\mathrm{emp}}(\mathbf{x})=\frac{1}{N}\sum _{i=1}^N{\delta }^n(x-x_i)$$

$$I\left(q\right)=D\left(p_{\mathrm{emp}},q\right)\to \underset{q\in {\mathcal{G}}_k}{\min }$$

$$R(q)=\sum _{i=k+1}^n{\lambda }_i(M_q)$$

$$I\left(q\right)+\lambda \Vert M_q{\Vert }_{n-k}\to \underset{q}{\min }$$

$$\begin{array}{r}{q}_{l+1}=\arg \underset{q}{\min }I\left(q\right)+\lambda \Vert S_q-S_l{\Vert }_{\ast }^2,\\ S_{l+1}=\arg \underset{S:\mathrm{rank}S=k}{\min }\Vert S_{q_{l+1}}-S{\Vert }_{\ast }^2.\end{array}$$

$$⟨T_f\otimes {\delta }^{n-k},\varphi ⟩={\int }_{{\mathbb{R}}^k}f({\mathbf{x}}_{1:k})\varphi ({\mathbf{x}}_{1:k},0_{n-k})d{\mathbf{x}}_{1:k}$$

$${\mathcal{G}}_k^{\prime }=\{(f\otimes {\delta }^{n-k}{)}_U|f\in {\mathcal{S}}^{\prime }({\mathbb{R}}^k),U\in \mathcal{O}(n)\},$$

$${\mathcal{G}}_k=\left\{(T_f\otimes {\delta }^{n-k}{)}_U|f\in \mathcal{S}({\mathbb{R}}^k),U\in \mathcal{O}(n)\right\},$$

$$\begin{array}{r}{\mathcal{F}}_k=\{T_r\mid r(\mathbf{x})=f(U\mathbf{x}),f\in \mathcal{S}({\mathbb{R}}^k),U\in {\mathbb{R}}^{k\times n},\mathrm{rank}(U)=k\}\end{array}$$

$$\begin{array}{r}⟨\mathcal{F}[g],\varphi ⟩=⟨g,\mathcal{F}[\varphi ]⟩=⟨T_f\otimes {\delta }^{n-k},{\int }_{{\mathbb{R}}^n}\varphi (\mathbf{y})e^{-\mathrm{i}{\mathbf{x}}^T\mathbf{y}}d\mathbf{y}⟩=\\ ⟨T_f,{\int }_{{\mathbb{R}}^n}\varphi (\mathbf{y})e^{-\mathrm{i}{\mathbf{x}}_{1:k}^T{\mathbf{y}}_{1:k}}d\mathbf{y}⟩={\int }_{{\mathbb{R}}^{n+k}}f({\mathbf{x}}_{1:k})\varphi (\mathbf{y})e^{-\mathrm{i}{\mathbf{x}}_{1:k}^T{\mathbf{y}}_{1:k}}d\mathbf{y}d{\mathbf{x}}_{1:k}=\\ {\int }_{{\mathbb{R}}^n}\hat{f}({\mathbf{y}}_{1:k})\varphi (\mathbf{y})d\mathbf{y}=⟨T_r,\varphi ⟩.\end{array}$$

$$\begin{array}{r}⟨\mathcal{F}[g_U],\varphi (\mathbf{x})⟩=⟨g_U,\mathcal{F}[\varphi ](\mathbf{x})⟩=⟨g,\mathcal{F}[\varphi ](U^T\mathbf{x})⟩=\\ ⟨g,\mathcal{F}[\varphi (U^T\mathbf{x})]⟩=⟨\mathcal{F}[g],\varphi (U^T\mathbf{x})⟩=⟨(\mathcal{F}[g]{)}_U,\varphi (\mathbf{x})⟩.\end{array}$$

$$\begin{array}{r}\left(\mathcal{F}[g_U]\right)=(T_r{)}_U=T_{r^{\prime }}\end{array}$$

$$\mathcal{Q}=\{g(U_k\mathbf{x})|g\in \mathcal{S}({\mathbb{R}}^k),U_k\in {\mathbb{R}}^{k\times n},U_k{U}_k^T=I_k\}$$

$$\begin{array}{r}\dim {\mathrm{span}}_{\mathbb{R}}\{x_{1}T,x_{2}T,\cdots ,x_{n}T\}\le k.\end{array}$$

$$\begin{array}{r}{\mathcal{P}}_k=\{\mu \in \mathcal{P}|\exists {\mathbf{v}}_1,\cdots ,{\mathbf{v}}_k\in {\mathbb{R}}^n,\forall A\in \mathcal{B}({\mathbb{R}}^n):\mu (A)=\mu (A\cap \mathrm{span}({\mathbf{v}}_1,\cdots ,{\mathbf{v}}_k))\}\end{array}$$

$$\mu \in \mathcal{P}\stackrel{\varphi }{\to }{\mathbb{E}}_{\mathbf{y}\sim \mu }K(\mathbf{x},\mathbf{y})=\int K(\mathbf{x},\mathbf{y})d\mu (\mathbf{y}).$$

$$d_{\text{textsc}MMD}(\mu ,\nu )=\Vert \varphi (\mu )-\varphi (\nu ){\Vert }_{{\mathcal{H}}_K}.$$

$$\underset{\nu \in {\mathcal{P}}_k}{\min }{d}_{\text{textsc}MMD}({\mu }_{\mathrm{data}},\nu )=\underset{\nu \in {\mathcal{P}}_k}{\min }\Vert \varphi ({\mu }_{\mathrm{data}})-\varphi (\nu ){\Vert }_{{\mathcal{H}}_K}$$

$$d_{\text{textsc}MMD}(\mu ,\nu )=\Vert \psi (\mu )-\psi (\nu ){\Vert }_{L_2({\mathbb{R}}^n)},$$

$$d_{\text{textsc}HM}(\mu ,\nu {)}^2=\sum _{s=1}^4\frac{{\lambda }_s}{n^s}\sum _{1\le i_1,\cdots ,i_s\le n}(m_{i_1\cdots i_s}-n_{i_1\cdots i_s}{)}^2$$

$$\underset{\nu \in {\mathcal{P}}_k}{\min }{d}_{\text{textsc}HM}({\mu }_{\mathrm{data}},\nu )$$

$$W_p(\mu ,\nu )=(\underset{\pi \in \Pi (\mu ,\nu )}{\inf }\int \Vert \mathbf{x}-\mathbf{y}{\Vert }^{p}d\pi {)}^{1/p}$$

$$\underset{\nu \in {\mathcal{P}}_k}{\min }{W}_p({\mu }_{\mathrm{data}},\nu )$$

$$\underset{f\in {\mathcal{F}}_k}{\inf }J(f)$$

$$\begin{array}{r}J(f)=\frac{1}{N}\sum _{i=1}^N{\mathbb{E}}_{\mathbf{ϵ}\sim N(0,{\upsilon }^2{I}_n)}|y_i-f({\mathbf{x}}_i+\mathbf{ϵ}){|}^2\end{array}$$

$$\begin{array}{r}J(f)=\frac{1}{N}\sum _{i=1}^N{\mathbb{E}}_{\mathbf{ϵ}\sim N(0,{\upsilon }^2{I}_n)}H\left(y_i,\frac{e^{f({\mathbf{x}}_i+\mathbf{ϵ})}}{1+e^{f({\mathbf{x}}_i+\mathbf{ϵ})}}\right)\end{array}$$

$$⟨f|M|g⟩={\iint }_{{\mathbb{R}}^n\times {\mathbb{R}}^n}f(\mathbf{x}{)}^{\ast }M(\mathbf{x},\mathbf{y})g(\mathbf{y})d\mathbf{x}d\mathbf{y}.$$

$$⟨f|M|g⟩\le \underset{\mathbf{x},\mathbf{y}}{\sup }|M(\mathbf{x},\mathbf{y})|\cdot \Vert f{\Vert }_{L_1}\Vert g{\Vert }_{L_1}<\infty .$$

$${\left[\begin{array}{c}⟨x_{i}f|M|x_{j}f⟩\end{array}\right]}_{1\le i,j\le n}$$

$$R(f)=\Vert M_f{\Vert }_{n-k}.$$

Full text

In order to reduce an optimization problem over distributions to an optimization problem over ordinary functions we introduce a nonnegative penalty function that “forces” the support of the model distribution to be $k$-dimensional. Then we present an algorithm for the minimization of the penalized objective, based on the infinite-dimensional low-rank optimization, which we call the alternating scheme. Also, we design an efficient approximate algorithm for a special case of the problem, where the distance between the empirical distribution and the model distribution is measured by Maximum Mean Discrepancy defined by a Mercer kernel of a certain type. We test our methods on four examples (three UDR and one SDR) using synthetic data and standard datasets.

Keywords: linear dimensionality reduction, sufficient dimension reduction, alternating scheme, tempered distribution.

1 Introduction

Linear dimension reduction (LDR) is a family of problems in data science that includes principal component analysis, factor analysis, linear multidimensional scaling, Fisher’s linear discriminant analysis, canonical correlations analysis, sufficient dimension reduction (SDR), maximum autocorrelation factors, slow feature analysis and more. In unsupervised dimension reduction (UDR) we are given a finite number of points in ${\mathbb{R}}^{n}$ (sampled according to some unknown distribution) and the goal is to find a “low-dimensional” manifold (e.g. an affine or a linear subspace) that approximates “the support” of the distribution. UDR, historically, was approached by linear methods and, therefore, has developed into a set of standard tools in data science. Though non-linear dimensionality reduction (a.k.a. the manifold learning) techniques gained a wide popularity in modern research, the potential of linear methods is far from being exhausted. For high-dimensional datasets, due to the phenomenon of concentration of measure [1], LDR often can give us an interpretable and low-dimensional representation of data. The linearity of a projection operator is a restrictive property that allows avoiding the overfitting in the dimension reduction (which is a key problem for the manifold learning techniques).

The LDR study field currently achieved a saturation level at which unifying frameworks for the problem become of special interest. One of such frameworks, that covers many cases of LDR, frames LDR as the optimization task over matrix manifolds such as the Stiefel manifold [2]. Elements of the Stiefel manifold $V_{k}({\mathbb{R}}^{n})$ are orthogonal $k$-frames $O\in{\mathbb{R}}^{n\times k}$ whose column space is the $k$-dimensional space onto which a dataset is projected. Different loss functions on $V_{k}({\mathbb{R}}^{n})$ define different versions of LDR. Table 1 of [2] lists fourteen common LDR techniques (such as principal component analysis, multi-dimensional scaling, linear discriminant analysis etc), nine of which are formulated over Stiefel manifolds. Such a general treatment allows to approach all LDR problems by a single algorithm, i.e. by an adaptation of the gradient descent method to Stiefel manifolds [3]. This adaptation consists of a series of projected gradient steps where a common gradient descent is followed by a projection onto a Stiefel manifold, which is equivalent to the computation of a singular value decomposition of a current point. Note that the Stiefel manifold is a non-convex set and even minimizing a convex function on such a manifold is an NP-hard task, in general. Although, in applications, the projected gradient descent method demonstrated a relatively fast convergence to good quality solutions.

The paper’s main contribution is a development of a novel view of LDR. First, an argument over which we search in an optimization task is not a $k$-frame, but a tempered distribution (which is a generalization of a probabilistic distribution) that is concentrated on a $k$-dimensional linear subspace of ${\mathbb{R}}^{n}$. Thus, an argument has a more complex structure, it includes not just a $k$-dimensional subspace, but also a distribution on that subspace. The justification of our optimization framework uses the theory of generalized functions, or tempered distributions [4, 5]. An important generalized function that cannot be represented as an ordinary function is the Dirac delta function, denoted $\delta$, and $\delta^{n}$ denotes its $n$-dimensional version.

This more general formulation allows us to analyze new types of objectives for LDR. In Section 4 we list four examples of such objectives that, to our knowledge, have not been considered in the LDR field so far. A notable specifics of such objectives is that, even for a fixed $k$-dimensional subspace $\mathcal{L}$, finding an optimal distribution supported in $\mathcal{L}$ is a non-trivial optimization task. In other words, our problems can not be simply reduced to the previous formalisms based on the Stiefel manifold, or the Grassmannian [6, 7].

Let us briefly describe an optimization problem that motivates the new formalism. Any dataset $\{{\mathbf{x}}_{i}\}_{i=1}^{N}\subseteq{\mathbb{R}}^{n}$ naturally corresponds to the distribution

[TABLE]

which, with some abuse of terminology, can be called the empirical probability density function. Based on that, UDR can be understood as a task whose goal is to approximate $p_{{\rm emp}}({\mathbf{x}})$ by $q({\mathbf{x}})$, where $q({\mathbf{x}})$ is a distribution whose density is supported in a $k$-dimensional linear subspace $\mathcal{L}\subseteq{\mathbb{R}}^{n}$. Note that a function whose density is supported in some low-dimensional subset of ${\mathbb{R}}^{n}$ is not an ordinary function. An exact definition of a set of such distributions, denoted by $\mathcal{G}_{k}$, is given in Section 3. To formulate an optimization task we additionally need a loss $D(p_{{\rm emp}},q)$ that measures the distance between the ground truth $p_{{\rm emp}}$ and a distribution $q$, that we search for. Thus, in our approach, the UDR problem is defined as

[TABLE]

under the condition that $q({\mathbf{x}})$ has a $k$-dimensional support. In most of our statements we do not consider any specific loss functions, though in our basic examples we deal with the Maximum Mean Discrepancy distance or the Wasserstein distance.

The UDR and SDR. Within our formalism the sufficient dimension reduction problem is tightly connected with the UDR problem. In the SDR, given supervised data, the goal is to find the so called effective subspace, defined by its orthogonal basis (or, a $k$-frame) $\{{\mathbf{w}}_{1},\cdots,{\mathbf{w}}_{k}\}\subseteq{\mathbb{R}}^{n}$ , such that the regression function can be searched in the form $g({\mathbf{w}}^{T}_{1}{\mathbf{x}},\cdots,{\mathbf{w}}^{T}_{k}{\mathbf{x}})$. In literature, these functions are known under different names, e.g. functions with low effective dimensionality [8], functions with active subspaces [9] and multi-ridge functions [10, 11]. In [12] it was shown that a method originally developed for the SDR can be turned into a UDR method, i.e. applied to unsupervised data, by simply setting an output to be equal to an input. In such methods for the SDR problem as the Sliced Inverse Regression [13], the Principal Hessian Direction [14], the Sliced Average Variance Estimation [15], an effective subspace is recovered from the Singular Value Decomposition applied to a certain matrix that is constructed from a training set in a straightforward way. Other methods, such as the Principal Fitted Components [16], the Likelihood Acquired Direction [17], the Kernel Dimensionality Reduction [18], are based on analytic expressions measuring the affinity of a $k$-dimensional subspace to the effective subspace. The second type of methods reduce the SDR problem to an optimization problem over the Stiefel manifold, or the Grassmanian. For other methods we refer to a tutorial on SDR methods [19]. Again, an important aspect of all these methods is that, given a fixed effective subspace, the regression function that predicts an output variable has a relatively straightforward structure and is not optimized by any additional supervised learning procedure. The key novelty that our framework brings to the SDR is that we suggest to search for an effective subspace and a regression function in a joint manner.

The key observation of our analysis, stated in Theorem 2, is that a class of functions of the form $g({\mathbf{w}}^{T}_{1}{\mathbf{x}},\cdots,{\mathbf{w}}^{T}_{k}{\mathbf{x}})$ can be characterized as functions whose Fourier transform is supported in the corresponding effective subspace. In other words, functions with an effective dimensionality $k$ are dual to $\mathcal{G}_{k}$ under the Fourier transform. Three examples of UDR problems that we give in Section 4 are cast as (2), whereas in the fourth example we formulate SDR as an optimization task with the search space dual to that of UDR (to distinguish our formulation from a general SDR problem we call it an SDR with optimized regression function). Thus, all four examples can be studied within our optimization framework.

Besides the problem setup we also suggest a general algorithm that tackles it. The basic idea of that algorithm, which we call the alternating scheme, instead of optimizing over $\mathcal{G}_{k}$, to optimize over ordinary functions with a penalty added to an objective that forces the ordinary function’s support to be low-dimensional.

The penalty based reformulation. The starting point of our approach is to reduce the task (2) to the minimization of $I\left(q\right)+\lambda R(q)$ over ordinary functions $q$. We define the penalty function $R(q)$ in such a way that forcing $R(q)$ to be small is equivalent to forcing “the support” of $q$ to be $k$-dimensional. Our definition of $R$ is based on using a positive definite kernel $M:{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}\to{\mathbb{C}}$.

First we note that $M$ defines a billinear form on pairs of (possibly, generalized) functions by $\langle f|M|g\rangle=\int_{{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}}f({\mathbf{x}})^{\ast}M({\mathbf{x}},{\mathbf{y}})g({\mathbf{y}})d{\mathbf{x}}d{\mathbf{y}}$. On a properly defined space of (generalized) functions, the billinear form $\langle\cdot|M|\cdot\rangle$ is the hermitian inner product, using which one can define distances and other geometrical notions on that space. Note that if $f$ and $g$ are probability density functions and $M$ is real-valued, the corresponding distance function, i.e. ${\rm dist}_{M}(f,g)=(\langle f|M|f\rangle-\langle g|M|g\rangle-2\langle f|M|g\rangle)^{1/2}$ coincides with the maximum mean discrepancy metric [20]. We define $R(q)$ as

[TABLE]

where $M_{q}={\rm Re}\begin{bmatrix}\langle x_{i}q({\mathbf{x}})|M|x_{j}q({\mathbf{x}})\rangle\end{bmatrix}_{i,j=\overline{1,n}}$ and $\lambda_{1}(M_{q})\geq\lambda_{2}(M_{q})\geq\cdots$ are ordered eigenvalues of the matrix $M_{q}$. The sum of all but first $k$ eigenvalues of a positive semidefinite matrix $A$ is a well-known penalty function, denoted by $\|A\|_{n-k}$ and called a Ky Fan $n-k$-antinorm. Applications of the Ky Fan $n-k$-antinorm to low-rank optimization problems can be found in [21, 22, 23, 24] and its properties are studied in [25].

Thus, we reduce the task (2) to

[TABLE]

over ordinary functions. An analysis that we make in Subsection 5.3 of Section 5 (based on theory of tempered distributions) shows that if the kernel $M$ is chosen from a class of so called proper kernels and the solution of (2) satisfies certain regularity conditions, the solution of the task (4) for $\lambda\to+\infty$ will approach the solution of (2).

The alternating scheme. The task (4) can be understood as an infinite dimensional low-rank optimization task in which the penalty term forces the matrix $M_{q}$ to be of rank $k$. In Section 6 we prove that $M_{q}=S_{q}S_{q}^{\dagger}$ where $S_{q}$ is a linear operator between a suitable space $\mathcal{H}$ and ${\mathbb{R}}^{n}$ that itself depends on $q$ linearly, and this automatically gives us that $R(q)=\min_{S}\|S_{q}-S\|^{2}_{\ast}$ where the minimum is taken over all operators between $\mathcal{H}$ and ${\mathbb{R}}^{n}$ of rank $k$ and $\|\cdot\|_{\ast}$ is a suitable norm on the space of bounded linear operators from $\mathcal{H}$ to ${\mathbb{R}}^{n}$.

Then, a natural idea to solve the task (4) is to present it as the joint minimum $\min\limits_{q}\min\limits_{S:{\rm rank\,}S=k}I\left(q\right)+\lambda\|S_{q}-S\|^{2}_{\ast}$ and to minimize the objective over $q$ and over $S$ of rank $k$ in an alternating fashion, i.e.

[TABLE]

This algorithm, called the alternating scheme, is suitable for a practical implementation due to the fact that the second step of it, i.e. the optimization over $S$ of rank $k$, is solvable analytically. In fact, $S_{l+1}$ is the Singular Value Decomposition of $S_{q_{l+1}}$ truncated at $k$-th term. In E we give an algorithm whose every step is equivalent to a corresponding step of the alternating algorithm, but it operates on Fourier transforms of functions rather than on functions of initial coordinates. Numerical specifications of the alternating scheme for different special cases of UDR/SDR problems are given in G,  H,  I and J. In Section 8 we describe results of our experiments with the alternating scheme that we conducted for various synthetic and practical datasets. As a result we conclude that the alternating scheme is a practical algorithm that can be applied to datasets of moderate size. For the SDR tasks its performance is comparable with classical algorithms.

An approximate algorithm for a special case. For a special case of the task (2), where the distance function is the Maximum Mean Discrepancy with the kernel of the form $K({\mathbf{x}},{\mathbf{y}})=({\mathbf{x}}\cdot{\mathbf{y}})H({\mathbf{x}},{\mathbf{y}})$ and $H$ is itself a Mercer kernel, we develop an approximate algorithm that can be applied to large datasets. In Section 7 we demonstrate that a solution with provable approximation ratios is given by the following simple procedure: given a dataset $\{{\mathbf{x}}_{i}\}_{i=1}^{N}$ we build a data matrix $X=[{\mathbf{x}}_{1},\cdots,{\mathbf{x}}_{N}]$, a Gram matrix $G=[H({\mathbf{x}}_{i},{\mathbf{x}}_{j})]$, and output first $k$ principal components of the matrix $XGX^{T}$. This algorithm is tested on Yale B dataset for the shadow/black removal and SBMnet datasets for the background modeling. In both applications, our approximate algorithm showed a performance comparable to the performance of other low-rank approximation algorithms.

The structure of the paper is as follows. In Section 2 we give some notations and define standard notions from functional analysis that we use throughout the paper. In Section 3 we formally define the search space in Problem 2, denoted $\mathcal{G}_{k}$, and an image of $\mathcal{G}_{k}$ under the Fourier transform, denoted $\mathcal{F}_{k}$. In Section 4 we formulate some UDR/SDR problems as optimization tasks over $\mathcal{G}_{k}$/$\mathcal{F}_{k}$. Instead of searching directly in a set of generalized functions, $\mathcal{G}_{k}$, in Section 5 we describe how we substitute an ordinary function for a distribution in the optimization task at the expence of adding a new penalty term to its objective, $\lambda R(f)$. Using a kernel $M({\mathbf{x}},{\mathbf{y}})$, Theorem 4 characterizes generalized $g\in\mathcal{G}_{k}$ as such $g$ for which the matrix of properly defined integrals $M_{g}={\rm Re}\begin{bmatrix}\iint_{{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}}x_{i}y_{j}g({\mathbf{x}})^{\ast}M({\mathbf{x}},{\mathbf{y}})g({\mathbf{y}})d{\mathbf{x}}d{\mathbf{y}}\end{bmatrix}_{i,j=\overline{1,n}}$ is of rank $k$. In Section 6 we suggest a method for solving $\min_{\phi}I(\phi)+\lambda R(\phi)$ which we call the alternating scheme. In Section 7 we describe a simple approximate algorithm for the task (2) in a special subcase of the Maximum Mean Discrepancy distance and prove some theoretical guarantees on the approximation ratio of this algorithm. Section 8 is dedicated to experiments with the alternating scheme on synthetic and real world data and with the approximate algorithm on the shadow/black removal and the background modeling applications. Proofs of all theorems and lemmas are given after their formulations, or can be found in the appendix.

1.1 Related work

As was already mentioned, another unifying framework for LDR tasks is suggested by [2] in which the basic search space is the Stiefel manifold $V_{k}({\mathbb{R}}^{n})$. The main advantage of the Stiefel manifold over ${\mathcal{G}}_{k}$ is that its elements are finite-dimensional. Because a distribution from ${\mathcal{G}}_{k}$ is an infinite-dimensional object, an optimization over ${\mathcal{G}}_{k}$ requires additional constructions to turn it into a finite-dimensional task. Both an optimization over ${\mathcal{G}}_{k}$ and over $V_{k}({\mathbb{R}}^{n})$ is typically hard: for a final point, at best one can guarantee that it is a local extremum. Promising aspects of ${\mathcal{G}}_{k}$ are: a) ${\mathcal{G}}_{k}$ allows to formulate a new class of objectives naturally on it, b) local extrema on ${\mathcal{G}}_{k}$ substantially differ from local extrema on $V_{k}({\mathbb{R}}^{n})$, because a local search over ${\mathcal{G}}_{k}$ uses more degrees of freedom.

There is plenty of literature on the SDR problem some of which was already mentioned. In [26] the Fourier transform was applied for estimating the effective subspace in SDR, implicitly using an analog of Theorem 2. The closest to ours is a recent approach of [27], where an effective subspace was computed in a two step process. First, given supervised data, a regression function was trained in the form of a neural network (with a general architecture), then the obtained regression function was approximated by another neural network with a bottleneck architecture (by which a low effective dimensionality is guaranteed by construction). Like in this approach, we train a regression function as a neural network, though we search for it and an effective subspace jointly. In our approach, it is a regularization term $R(f)$ that forces the neural network to have a low effective dimensionality.

Using Ky Fan $k$-antinorm as a regularizer for the matrix completion problem has been suggested by [21] and further developed in [22, 23, 24]. Unlike this chain of works, we formulate an infinite-dimensional task and our regularizer $R(f)=\|M_{f}\|_{n-k}$ is a sum of smallest $n-k$ squared singular values of the infinite-dimensional operator $S_{f}$ where $S_{f}$ depends on $f$ linearly and $M_{f}=S_{f}S_{f}^{\dagger}$. Thus, our algorithms are substantially different from algorithms designed within the latter approach. The idea of alternating two basic stages, convex optimization and SVD, is ubiquitous in low-rank optimization, see e.g. [28, 29].

2 Preliminaries and notations

Throughout this paper we use standard terminology and notation from functional analysis. For details one can address the textbook on the theory of distributions [30]. The Schwartz space, denoted by $\mathcal{S}({\mathbb{R}}^{n})$, is a space of infinitely differentiable functions $f:{\mathbb{R}}^{n}\rightarrow{\mathbb{C}}$ such that $\forall\alpha,\beta\in{\mathbb{N}}^{n},\sup_{{\mathbf{x}}\in\mathbb{R}^{n}}$ $|{\mathbf{x}}^{\alpha}D^{\beta}f({\mathbf{x}})|<\infty$, and equipped with standard topology. Its dual space is denoted by $\mathcal{S^{\prime}}({\mathbb{R}}^{n})$ and is equipped with weak topology. For a tempered distribution $T\in\mathcal{S^{\prime}}({\mathbb{R}}^{n})$ and $\phi\in\mathcal{S}({\mathbb{R}}^{n})$, $\langle T,\phi\rangle$ denotes $T(\phi)$. Thus, for a sequence $\{f_{s}\}\subseteq\mathcal{S^{\prime}}({\mathbb{R}}^{n})$ and $f\in\mathcal{S^{\prime}}({\mathbb{R}}^{n})$, $\lim_{s\rightarrow\infty}f_{s}=f$ (or, $f_{s}\rightarrow^{\ast}f$) means that $\lim_{s\rightarrow\infty}\langle f_{s},\phi\rangle=\langle f,\phi\rangle$ for any $\phi\in\mathcal{S}({\mathbb{R}}^{n})$. For a sequence $\{f_{s}\}^{\infty}_{s=1}\subseteq{\mathcal{S}}^{\prime}({\mathbb{R}}^{n})$, $\mathop{\rm Lim}\limits_{s\rightarrow\infty}f_{s}$ denotes a set of points $f\in{\mathcal{S}}^{\prime}({\mathbb{R}}^{n})$, such that there exists a growing sequence $\{s_{i}\}\subseteq{\mathbb{N}}$ and $\lim_{i\rightarrow\infty}f_{s_{i}}=f$. The Fourier and inverse Fourier transforms are denoted by $\mathcal{F},\mathcal{F}^{-1}:\mathcal{S^{\prime}}({\mathbb{R}}^{n})\rightarrow\mathcal{S^{\prime}}({\mathbb{R}}^{n})$. For brevity, we denote $\mathcal{F}[f]$ by $\hat{f}$. If all required conditions are satisfied, an integrable $f:{\mathbb{R}}^{n}\rightarrow{\mathbb{C}}$ (or, a Borel measure $\mu$ on ${\mathbb{R}}^{n}$) is used as the tempered distribution $T_{f}$ (or, $T_{\mu}$) where $\langle T_{f},\phi\rangle=\int_{{\mathbb{R}}^{n}}f({\mathbf{x}})\phi({\mathbf{x}})d{\mathbf{x}}$ (or, $\langle T_{\mu},\phi\rangle=\int_{{\mathbb{R}}^{n}}\phi({\mathbf{x}})d\mu$). For $\Omega\subseteq\mathcal{S}({\mathbb{R}}^{n})$, $\overline{\Omega}$ denotes the sequential closure of $\Omega$ with respect to standard topology of $\mathcal{S}({\mathbb{R}}^{n})$. For $\Omega\subseteq\mathcal{S^{\prime}}({\mathbb{R}}^{n})$, $\overline{\Omega}^{\ast}$ denotes the sequential closure of $\Omega$ with respect to weak topology of $\mathcal{S^{\prime}}({\mathbb{R}}^{n})$. For $\psi\in\mathcal{S}({\mathbb{R}}^{n}),T\in\mathcal{S^{\prime}}({\mathbb{R}}^{n})$, the convolution is defined as a tempered distribution $\psi\ast T$ such that $\langle\psi\ast T,\phi\rangle=\langle T,\tilde{\psi}\ast\phi\rangle$ where $\tilde{\psi}({\mathbf{x}})=\psi(-{\mathbf{x}})$. If $T\in\mathcal{S^{\prime}}({\mathbb{R}}^{n})$ and a function $\psi$ is such that $\psi({\mathbf{x}})\phi({\mathbf{x}})\in\mathcal{S}({\mathbb{R}}^{n})$ whenever $\phi({\mathbf{x}})\in\mathcal{S}({\mathbb{R}}^{n})$, then the multiplication $\psi T$ is defined by $\langle\psi T,\phi\rangle=\langle T,\psi\phi\rangle$. Given a measure $\mu$, by $L_{2,\mu}({\mathbb{R}}^{n})$ we denote the complex $L_{2}$-space with the inner product $\langle u,v\rangle_{L_{2,\mu}}=\int u({\mathbf{x}})^{\ast}v({\mathbf{x}})d\mu$. The induced norm is then $\|{\mathbf{u}}\|_{L_{2,\mu}}=\sqrt{\langle{\mathbf{u}},{\mathbf{u}}\rangle_{L_{2,\mu}}}$. If $d\mu=p({\mathbf{x}})d{\mathbf{x}}$, then $L_{2,\mu}$ is denoted by $L_{2,p}$. A set of infinitely differentiable functions in ${\mathbb{R}}^{n}$ is denoted by $C^{\infty}({\mathbb{R}}^{n})$. A set of infinitely differentiable functions with compact support in ${\mathbb{R}}^{n}$ is denoted by $C_{c}^{\infty}({\mathbb{R}}^{n})$. If $T$ is a topological space, then a subset $S\subseteq T$ is said to be dense in $T$ if the sequential closure of $S$ is equal to $T$. For a square matrix $A$, ${\rm Tr}(A)$ denotes its trace and for an arbitrary matrix, $\|A\|_{F}\overset{def}{=}\sqrt{{\rm Tr}(A^{T}A)}$. The identity matrix of size $n$ is denoted by $I_{n}$. The notation $f\propto g$ means $f=cg$ where $c$ is some universal constant.

3 Basic function classes

To formalize distributions supported in a $k$-dimensional subspace, we need a number of standard definitions [31]. For $\phi_{1}\in{\mathcal{S}}({\mathbb{R}}^{k})$ and $\phi_{2}\in{\mathcal{S}}({\mathbb{R}}^{n-k})$, their tensor product is the function $\phi_{1}\otimes\phi_{2}\in{\mathcal{S}}({\mathbb{R}}^{n})$ such that $(\phi_{1}\otimes\phi_{2})({\mathbf{x}},{\mathbf{y}})=\phi_{1}({\mathbf{x}})\phi_{2}({\mathbf{y}})$. The span of $\{\phi_{1}\otimes\phi_{2}|\phi_{1}\in{\mathcal{S}}({\mathbb{R}}^{k}),\phi_{2}\in{\mathcal{S}}({\mathbb{R}}^{n-k})\}$, denoted by ${\mathcal{S}}({\mathbb{R}}^{k})\otimes{\mathcal{S}}({\mathbb{R}}^{n-k})$, is called the tensor product of ${\mathcal{S}}({\mathbb{R}}^{k})$ and ${\mathcal{S}}({\mathbb{R}}^{n-k})$. For $g_{1}\in{\mathcal{S}}^{\prime}({\mathbb{R}}^{k})$ and $g_{2}\in{\mathcal{S}}^{\prime}({\mathbb{R}}^{n-k})$, their tensor product is defined by the following rule: $\langle g_{1}\otimes g_{2},\phi_{1}\otimes\phi_{2}\rangle=\langle g_{1},\phi_{1}\rangle\langle g_{2},\phi_{2}\rangle$ for any $\phi_{1}\in{\mathcal{S}}({\mathbb{R}}^{k}),\phi_{2}\in{\mathcal{S}}({\mathbb{R}}^{n-k})$. Since $\overline{{\mathcal{S}}({\mathbb{R}}^{k})\otimes{\mathcal{S}}({\mathbb{R}}^{n-k})}={\mathcal{S}}({\mathbb{R}}^{n})$, there is only one distribution $g_{1}\otimes g_{2}\in{\mathcal{S}}^{\prime}({\mathbb{R}}^{n})$ that satisfies the identity.

An example of a generalized function, whose density is concentrated in a $k$-dimensional subspace, is any distribution that can be represented as g\otimes\delta^{n-k}\overset{def}{=}g\otimes\underbrace{\delta\otimes\cdots\otimes\delta}_{\text{n-k times}} where $g\in{\mathcal{S}}^{\prime}({\mathbb{R}}^{k})$. If $g=T_{f}$, where $f:{\mathbb{R}}^{k}\rightarrow{\mathbb{R}}$ is an ordinary function, then $g\otimes\delta^{n-k}$ can be understood as a generalized function whose density is concentrated in a subspace $\{{\mathbf{x}}\in{\mathbb{R}}^{n}|x_{i}=0,i>k\}$ and equals $f({\mathbf{x}}_{1:k})$. It can be shown that the distribution acts on $\phi\in{\mathcal{S}}({\mathbb{R}}^{n})$ in the following way:

[TABLE]

Now to generalize the latter definition to any $k$-dimensional subspace we have to introduce a change of variables in tempered distributions.

Let $g\in{\mathcal{S}}^{\prime}({\mathbb{R}}^{n})$ and $U\in{\mathbb{R}}^{n\times n}$ be an orthogonal matrix, i.e. $U^{T}U=I_{n}$. Then, $g_{U}\in{\mathcal{S}}^{\prime}({\mathbb{R}}^{n})$ is defined by the rule: $\langle g_{U},\phi\rangle=\langle g,\psi\rangle$ where $\psi({\mathbf{x}})=\phi(U^{T}{\mathbf{x}})$. If $g=T_{f}$, the latter definition gives $g_{U}=T_{f^{\prime}}$ where $f^{\prime}({\mathbf{x}})=f(U{\mathbf{x}})$. Now, we define classes of tempered distributions:

[TABLE]

[TABLE]

and

[TABLE]

where ${\mathcal{O}}(n)=\{U\in{\mathbb{R}}^{n\times n}\mid U^{T}U=I_{n}\}$. The first two classes are related as:

Theorem 1**.**

$\mathcal{G}^{\prime}_{k}=\overline{\mathcal{G}_{k}}^{\ast}$.

The last two classes are isomorphic under the Fourier transform.

Theorem 2**.**

$\mathcal{F}[\mathcal{G}_{k}]=\mathcal{F}_{k}$* and $\mathcal{F}^{-1}[\mathcal{F}_{k}]=\mathcal{G}_{k}$.*

Proof.

Let us prove first that if $g=T_{f}\otimes\delta^{n-k}$, then $\mathcal{F}[g]=T_{r}$, where $r({\mathbf{x}})=\hat{f}({\mathbf{x}}_{1:k}),{\mathbf{x}}\in{\mathbb{R}}^{n}$. For that we have to prove that $\langle\mathcal{F}[g],\phi\rangle=\langle T_{r},\phi\rangle$ for any $\phi\in{\mathcal{S}}({\mathbb{R}}^{n})$. Indeed,

[TABLE]

Let us calculate the image of $\mathcal{G}_{k}$ under the Fourier transform. It is easy to see that for any $g\in{\mathcal{S}}^{\prime}({\mathbb{R}}^{n}),\phi\in{\mathcal{S}}({\mathbb{R}}^{n})$ and orthogonal $U\in{\mathbb{R}}^{n\times n}$ we have:

[TABLE]

Therefore, $\mathcal{F}[g_{U}]=(\mathcal{F}[g])_{U}$. Thus, if $g=T_{f}\otimes\delta^{n-k}$, then

[TABLE]

where $r^{\prime}({\mathbf{x}})=r(U{\mathbf{x}})=\hat{f}(U_{k}{\mathbf{x}})$ and $U_{k}\in{\mathbb{R}}^{k\times n}$ is a matrix consisting of first $k$ rows of $U$. Thus, $T_{r^{\prime}}\in\mathcal{F}_{k}$.

Let us show that by varying $f\in{\mathcal{S}}({\mathbb{R}}^{k})$ and $U$ in the expression $\hat{f}(U_{k}{\mathbf{x}})$ we can obtain any function from $\mathcal{F}_{k}$. For this, it is enough to show that $\mathcal{F}_{k}$ is equivalent to the following set of functions:

[TABLE]

The fact $\mathcal{Q}\subseteq\mathcal{F}_{k}$ is obvious. Let us now prove that $\mathcal{Q}\supseteq\{g(P{\mathbf{x}})|g\in{\mathcal{S}}({\mathbb{R}}^{k}),P\in{\mathbb{R}}^{k\times n},\operatorname{rank}P=k\}=\mathcal{F}_{k}$. Indeed, if $f({\mathbf{x}})=g(P{\mathbf{x}})$, then $f({\mathbf{x}})=g^{\prime}(U_{k}{\mathbf{x}})$ where $U_{k}=(PP^{T})^{-1/2}P$ and $g^{\prime}({\mathbf{y}})=g((PP^{T})^{1/2}{\mathbf{y}})$. By construction, $U_{k}U_{k}^{T}=I_{k}$ and $g^{\prime}\in{\mathcal{S}}({\mathbb{R}}^{k})$. Thus, $\mathcal{Q}=\mathcal{F}_{k}$.

Therefore, $\mathcal{F}[\mathcal{G}_{k}]=\mathcal{F}_{k}$, and from the bijectivity of the Fourier transform we obtain $\mathcal{F}^{-1}[\mathcal{F}_{k}]=\mathcal{G}_{k}$. ∎

For any collection $f_{1},\cdots,f_{l}\in{\mathcal{S}}^{\prime}({\mathbb{R}}^{n})$, ${\rm span}_{{\mathbb{R}}}\{f_{i}\}^{l}_{1}$ denotes $\{\sum_{i=1}^{l}\lambda_{i}f_{i}|\lambda_{i}\in{\mathbb{R}}\}\subseteq{\mathcal{S}}^{\prime}({\mathbb{R}}^{n})$, which is a linear space over ${\mathbb{R}}$. The set $\mathcal{G}^{\prime}_{k}$ has the following simple characterization:

Theorem 3**.**

For any $T\in{\mathcal{S}}^{\prime}({\mathbb{R}}^{n})$, $T\in\mathcal{G}^{\prime}_{k}$ if and only if

[TABLE]

Informally, the theorem holds because any linear dependency $\alpha_{1}x_{1}T+\cdots+\alpha_{n}x_{n}T=0$ over ${\mathbb{R}}$ implies that if $\alpha_{1}x_{1}+\cdots+\alpha_{n}x_{n}\neq 0$, then $T=0$. This is equivalent to a statement that the support of $T$ is concentrated on a subspace $\alpha_{1}x_{1}+\cdots+\alpha_{n}x_{n}=0$. If $\dim{\rm span}_{{\mathbb{R}}}\{x_{1}T,x_{2}T,\cdots,x_{n}T\}\leq k$, then one can find $n-k$ such dependencies, which means that the support of $T$ is $k$-dimensional.

Let ${\mathcal{B}}({\mathbb{R}}^{n})$ denote the Borel sigma-algebra on ${\mathbb{R}}^{n}$ and $\mathcal{P}$ denote a set of all Borel probability measures on ${\mathbb{R}}^{n}$. Let us now define

[TABLE]

i.e. $\mathcal{P}_{k}$ is a set of probability measures with all probability concentrated in some subspace ${\rm span}({\mathbf{v}}_{1},\cdots,{\mathbf{v}}_{k})$ whose dimension is not greater than $k$. It is easy to see that $T_{\mu}\in\mathcal{G}^{\prime}_{k}$ for any $\mu\in\mathcal{P}_{k}$.

4 Examples of LDR formulations

Maximum mean discrepancy PCA (MMD-PCA) Let $K:{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}\to{\mathbb{R}}$ be a continuous Mercer kernel, and $\mathcal{H}_{K}$ be a reproducing kernel Hilbert space (RKHS) defined by $K$. The kernel $K({\mathbf{x}},{\mathbf{y}})$ defines the so-called kernel embedding of probability measures $\phi$ [32]:

[TABLE]

The Maximum Mean Discrepancy (MMD) distance [20] is defined as the distance induced by metrics on $\mathcal{H}_{K}$, i.e. for two probability measures $\mu,\nu\in\mathcal{P}$,

[TABLE]

Let ${\mathbf{x}}_{1},\cdots,{\mathbf{x}}_{N}\in{\mathbb{R}}^{n}$ be the dataset of points. This dataset defines the empirical probabilistic measure $\mu_{\rm{data}}$ that corresponds to the tempered distribution $T_{\mu_{\rm{data}}}=\frac{1}{N}\sum_{i=1}^{N}\delta^{n}({\mathbf{x}}-{\mathbf{x}}_{i})$. We shall study a method concurrent to PCA that is based on solving the following problem:

[TABLE]

i.e. we shall attempt to approximate the empirical probabilistic measure $\mu_{\rm{data}}$ with another probabilistic measure $\nu$ which is supported in some $k$-dimensional subspace of ${\mathbb{R}}^{n}$. To our knowledge, the task (17) has not been yet considered in the research field of LDR.

Example 1** (Gaussian MMD-PCA).**

Let $k({\mathbf{x}})=G^{n}_{h}({\mathbf{x}})$ where $G^{n}_{h}({\mathbf{x}})=\frac{e^{-\frac{\|{\mathbf{x}}\|^{2}}{2h^{2}}}}{(2\pi h^{2})^{n/2}}$ is the radial Gaussian kernel on ${\mathbb{R}}^{n}$ and $K({\mathbf{x}},{\mathbf{y}})=(k\ast k)({\mathbf{x}}-{\mathbf{y}})=G^{n}_{2h}({\mathbf{x}})$. For such a kernel, we have

[TABLE]

where $\psi(\mu)=\int k({\mathbf{x}}-{\mathbf{y}})d\mu({\mathbf{y}})$ is just a smoothing of the distribution $\mu$ via the Weierstrass trasform.

In this example, as $h\to+0$, the optimal measure $\nu^{\ast}=\arg\min_{\nu\in\mathcal{P}_{k}}\|\psi(\mu_{\rm{data}})-\psi(\nu)\|_{L_{2}({\mathbb{R}}^{n})}$ is supported in a $k$-dimensional subspace that contains the largest possible number of points from $\{{\mathbf{x}}_{1},\cdots,{\mathbf{x}}_{N}\}$. Khachiyan demonstrated [33] that the following problem is NP-hard: given $\{{\mathbf{x}}_{1},\cdots,{\mathbf{x}}_{N}\}\subseteq{\mathbb{R}}^{n}$, find an $n-1$-dimensional subspace of ${\mathbb{R}}^{n}$ that contains at least $(1-\varepsilon)(1-\frac{1}{n})N$ points from the dataset. This indicates that in the regime $h\to+0$, the task (17) is NP-hard. In other words, it is unlikely that the task admits an efficient algorithm, in general. In G we describe an algorithm for the Gaussian MMD-PCA.

The higher moments PCA (HM-MMD-PCA) Another natural approach to measuring the similarity of two distributions is based on the difference between moments:

[TABLE]

where $m_{i_{1}\cdots i_{s}}={\mathbb{E}}_{{\mathbf{X}}\sim\mu}\left[{\mathbf{X}}[{i_{1}}]\cdots{\mathbf{X}}[{i_{s}}]\right]$ and $n_{i_{1}\cdots i_{s}}={\mathbb{E}}_{{\mathbf{X}}\sim\nu}\left[{\mathbf{X}}[{i_{1}}]\cdots{\mathbf{X}}[{i_{s}}]\right]$ are corresponding moments. The positive parameters $\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$, $\lambda_{4}$ are chosen to fix the relative importance of the mean, the co-variance, the co-skewness and the co-kurtosis.

Thus, we will be interested in the following optimization task (analogous to 17):

[TABLE]

If we set $\lambda_{2}=1$ and $\lambda_{1}=\lambda_{3}=\lambda_{4}=0$, then the solution of the task (20) coinsides with the solution of the classical PCA. Let us briefly demonstrate that. Let $X=[{\mathbf{x}}_{1},\cdots,{\mathbf{x}}_{N}]$ be the data matrix, $Y=[{\mathbf{y}}_{1},\cdots,{\mathbf{y}}_{N}]$ be the SVD of $X$ truncated at $k$-th term, $\sigma_{i}(X)$ be an $i$th singular value of $X$. By $\mu_{\rm{pca}}$ we denote a probabilistic measure concentrated in points $\{{\mathbf{y}}_{i}\}_{i=1}^{N}$. In that case we have $d_{\textsc{HM}}(\mu_{\rm{data}},\mu_{\rm{pca}})^{2}=\|\frac{1}{N}XX^{T}-\frac{1}{N}YY^{T}\|_{F}^{2}=\frac{1}{N^{2}}\sum_{i=k+1}^{\min\{N,n\}}\sigma_{i}^{4}(X)$. But for any $\nu\in\mathcal{P}_{k}$ the covariance matrix ${\rm cov}(\nu)=[{\mathbb{E}}_{{\mathbf{x}}\sim\nu}x_{i}x_{j}]_{i,j\in[n]}$ is of rank $k$. Therefore, by Eckart-Young-Mirsky’s theorem, we have $d_{\textsc{HM}}(\mu_{\rm{data}},\nu)^{2}=\|\frac{1}{N}XX^{T}-{\rm cov}(\nu)\|^{2}_{F}\geq\frac{1}{N^{2}}\sum_{i=k+1}^{\min\{N,n\}}\sigma_{i}^{4}(X)$. Thus, the minimum of $d_{\textsc{HM}}(\mu_{\rm{data}},\nu)$ is attained at $\nu=\mu_{\rm{pca}}$.

Thus, the task (20) can be considered as a direct generalization of PCA that takes into account higher moments. Note that the distance based on higher moments is a special case of maximum mean discrepance metric, where $K({\mathbf{x}},{\mathbf{y}})=\sum_{s=1}^{4}\frac{\lambda_{s}}{n^{s}}({\mathbf{x}}\cdot{\mathbf{y}})^{s}$. That is why we denote the task as HM-MMD-PCA. In Section 7 we prove that there is an efficient 2-approximating algorithm for the HM-MMD-PCA. In H we additionally describe another algorithm for the HM-MMD-PCA based on a generic alternating scheme.

Wasserstein distance PCA (WD-PCA) Another significant distance between probability measures with the origins in the transport theory is the Wasserstein distance (see [34]).

Let $({\mathbb{R}}^{n},\|\cdot\|)$ be a Banach space and $p\geq 1$. Between any two Borel probability measures $\mu,\nu$ on ${\mathbb{R}}^{n}$ with $\int\|{\mathbf{x}}\|^{p}d\mu<\infty$ and $\int\|{\mathbf{x}}\|^{p}d\nu<\infty$ the $p$th Wasserstein distance is:

[TABLE]

where $\Pi(\mu,\nu)$ is a set of all couplings of $\mu$ and $\nu$. The Wasserstein distance defines another version of LDR problem:

[TABLE]

In the B one can find proofs that in the case of $l_{1}$ norm $\|{\mathbf{x}}\|=\sum_{i}|x_{i}|$ and $p=1$, the task (22) corresponds to the well-studied robust PCA problem [35]. If, instead of the $l_{1}$-norm, we use the $l_{2}$-norm and set $p=1$, this leads to another well-studied task, which is known as the outlier pursuit problem [36, 37]. In the case of the $l_{2}$-norm and a general $p\geq 1$ we obtain the $l_{p}$ subspace approximation problem [38, 39]. Note that, except for the $l_{2}$ subspace approximation problem, all these problems are NP-hard. In I we describe an algorithm for the WD-PCA in the case of $l_{2}$-norm and $p=1$.

Sufficient dimension reduction with optimized regression function (SDR-ORF). Given a labeled dataset $\{({\mathbf{x}}_{i},y_{i})\}_{i=1}^{N}$ where ${\mathbf{x}}_{i}\in{\mathbb{R}}^{n},y_{i}\in{\mathcal{C}}$ (${\mathcal{C}}$ is a finite set of classes for a classification, or ${\mathbb{R}}$ for a regression problem), the sufficient dimension reduction problem can be informally described as a problem of finding vectors ${\mathbf{w}}_{1},\cdots,{\mathbf{w}}_{k}\in{\mathbb{R}}^{n}$ such that conditional distributions satisfy $p(y|{\mathbf{w}}^{T}_{1}{\mathbf{x}},\cdots,{\mathbf{w}}^{T}_{k}{\mathbf{x}})\approx p(y|{\mathbf{x}})$ (possibly, under some additional assumptions on the form of $p(y|{\mathbf{x}})$).

We formulate the SDR-ORF problem as an optimization task:

[TABLE]

The object $f:{\mathbb{R}}^{n}\rightarrow{\mathbb{R}}$ is a smooth real-valued function. We assume that $f$ is a candidate for the regression function and $J(f)$ is a cost function that values how strongly $f$ fits in this role. In practice for the regression case and for the binary classification case with 0-1 outputs we use the following cost functions correspondingly:

[TABLE]

and

[TABLE]

where $H(y,p)=-y\log p-(1-y)\log(1-p)$ and $\upsilon>0$ is a parameter.

By requiring $f\in\mathcal{F}_{k}$, we assume that the regression function $f$ satisfies (for $k$ fixed in advance): $f({\mathbf{x}})=g({\mathbf{w}}^{T}_{1}{\mathbf{x}},\cdots,{\mathbf{w}}^{T}_{k}{\mathbf{x}})$, where ${\mathbf{w}}_{1},\cdots,{\mathbf{w}}_{k}\in{\mathbb{R}}^{n}$. Thus, given an input ${\mathbf{x}}$, an output of $f$ depends on the projection of ${\mathbf{x}}$ onto ${\rm span}({\mathbf{w}}_{1},\cdots,{\mathbf{w}}_{k})$. The set ${\rm span}({\mathbf{w}}_{1},\cdots,{\mathbf{w}}_{k})$ is called the effective subspace. In J we describe an algorithm for the SDR-ORF problem.

5 Reduction of the optimization problem to ordinary functions

The central problem that our paper addresses is the optimization of an objective function over $\mathcal{G}^{\prime}_{k}$? In this section we suggest an approach based on penalty functions and kernels.

5.1 The definition of the penalty function

In this subsection we introduce a penalty function $R(f)$. Let $M:{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}\to{\mathbb{C}}$111Throughout the paper the kernel that induces the MMD distance is denoted by $K$ and the kernel that is used to define a penalty is denoted by $M$. be some bounded function such that $[M({\mathbf{z}}_{i},{\mathbf{z}}_{j})]_{i,j\in[x]}$ is a positive semidefinite matrix for any $\{{\mathbf{z}}_{i}\}_{i\in[x]}\subseteq{\mathbb{R}}^{n},x\in{\mathbb{N}}$. For $f,g:{\mathbb{R}}^{n}\rightarrow{\mathbb{C}}$ let us denote

[TABLE]

For $f,g\in L_{1}({\mathbb{R}}^{n})$,

[TABLE]

For general $f,g\in{\mathcal{S}}^{\prime}({\mathbb{R}}^{n})$ the expression $\langle f|M|g\rangle$ is defined if there are $f_{\epsilon},g_{\epsilon}\in L_{1}({\mathbb{R}}^{n})$ such that $T_{f_{\epsilon}}=f\ast G^{n}_{\epsilon}$, $T_{g_{\epsilon}}=g\ast G^{n}_{\epsilon}$ and $\lim_{\epsilon\rightarrow 0}\langle f_{\epsilon}|M|g_{\epsilon}\rangle=A<\infty$. Then, $\langle f|M|g\rangle\overset{def}{=}A$. For example, for continuous $M$ we have $\langle\delta^{n}|M|\delta^{n}\rangle=M(0,0)$.

One can build a Gram matrix from the collection of functions $\{x_{i}f\}_{i=1}^{n}$, $\begin{bmatrix}\langle x_{i}f|M|x_{j}f\rangle\end{bmatrix}_{1\leq i,j\leq n}$. Let us denote a real part of the Gram matrix

[TABLE]

by $M_{f}$.

Theorem 3 concludes, from $f\in\mathcal{G}_{k}$, that $\dim{\rm span}_{\mathbb{R}}\{x_{1}f,x_{2}f,\cdots,x_{n}f\}\leq k$.

Theorem 4**.**

Let $M({\mathbf{x}},{\mathbf{y}})$ be a bounded Lipschitz function. If $f=(T_{g}\otimes\delta^{n-k})_{U}\in\mathcal{G}^{\prime}_{k}$ is such that $\{x_{i}g\}_{i=1}^{k}\subseteq L_{1}({\mathbb{R}}^{k})$, then $\langle x_{i}f|M|x_{j}f\rangle$ is defined and $\operatorname{rank}M_{f}\leq k$.

Definition 1**.**

Let $A\in{\mathbb{R}}^{n\times n}$ be a positive semidefinite matrix with eigenvalues $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{n}$ (with counting multiplicities). Then, the Ky Fan $k$–anti-norm of $A$ is $\|A\|_{k}=\sum_{i=1}^{k}\lambda_{n+1-k}$.

Let

[TABLE]

By construction, by penalizing the value of $R(f)$, we enforce $M_{f}$ to be close to some matrix of rank $k$. Equivalently, we enforce a real part of the Gram matrix of $\{x_{i}f\}_{i\in[n]}$ to be of close to a rank $k$ matrix. By Theorem 3, the condition $\dim{\rm span}_{\mathbb{R}}\left\{x_{1}f,x_{2}f,\cdots,x_{n}f\right\}\leq k$ implies $f\in\mathcal{G}^{\prime}_{k}$, therefore, we enforce $f$ to be close to some function from $\mathcal{G}^{\prime}_{k}$. In the next section we will justify the latter informal logic by reducing the optimization over $\mathcal{G}^{\prime}_{k}$ to the optimization over ordinary functions with the penalty function $R(f)$.

5.2 Proper kernels

For a function $M({\mathbf{x}},{\mathbf{y}}):{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}\rightarrow{\mathbb{C}}$, let us denote by ${\rm O}_{M}$ a linear operator between ${\rm Dom}({\rm O}_{M})$ and $L_{2}({\mathbb{R}}^{n})$ given by ${\rm O}_{M}[f]=\int_{{\mathbb{R}}^{n}}M({\mathbf{x}},{\mathbf{y}})f({\mathbf{y}})d{\mathbf{y}}$ where ${\rm Dom}({\rm O}_{M})=\{f\in L_{2}({\mathbb{R}}^{n})\mid{\rm O}_{M}[f]\in L_{2}({\mathbb{R}}^{n})\}$. For any operator $O$ between spaces $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$, we denote its range by ${\rm Range\,}[O]=\{O(x)|x\in\mathcal{H}_{1}\}$.

Definition 2**.**

The function $M({\mathbf{x}},{\mathbf{y}}):{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}\rightarrow{\mathbb{C}}$ is called the proper kernel if and only if

${\rm O}_{M}:L_{2}({\mathbb{R}}^{n})\to L_{2}({\mathbb{R}}^{n})$* is a properly defined, strictly positive and self-adjoint operator,* 2. 2.

$\max_{{\mathbf{x}},{\mathbf{y}}}|M({\mathbf{x}},{\mathbf{y}})|<\infty$, 3. 3.

$\overline{{\rm Range\,}[{\rm O}_{M}]\cap{\mathcal{S}}({\mathbb{R}}^{n})}={\mathcal{S}}({\mathbb{R}}^{n})$.

Note that the latter definition implies that $M({\mathbf{y}},{\mathbf{x}})=M({\mathbf{x}},{\mathbf{y}})^{\ast}$ (modulo some null set) and $\langle f,{\rm O}_{M}[f]\rangle_{L_{2}({\mathbb{R}}^{n})}>0,\forall f\in L_{2}({\mathbb{R}}^{n}),f\neq{\mathbf{0}}$.

Example 2**.**

The Gaussian kernel is of special interest in applications: $M({\mathbf{x}},{\mathbf{y}})=G^{n}_{\sigma}({\mathbf{x}}-{\mathbf{y}}).$

It is captured by the following lemma:

Lemma 1**.**

If $\zeta,\hat{\zeta}\in C({\mathbb{R}}^{n})$ are bounded, $\forall{\mathbf{x}}\,\,\hat{\zeta}({\mathbf{x}})>0$, then $M({\mathbf{x}},{\mathbf{y}})=\zeta({\mathbf{x}}-{\mathbf{y}})$ is a proper kernel.

Proof.

Verification of the first three conditions is easy, so we only check the fourth condition. Let us denote linear operators $C_{\zeta}[f]=\zeta\ast f$ and $O_{g}[f]({\mathbf{x}})=g({\mathbf{x}})f({\mathbf{x}})$. Then we have $\mathcal{F}[C_{\zeta}[L_{2}({\mathbb{R}}^{n})]]=O_{\hat{\zeta}}[L_{2}({\mathbb{R}}^{n})]\supseteq C_{c}^{\infty}({\mathbb{R}}^{n})$. Therefore, ${\rm Range\,}[{\rm O}_{M}]=C_{\zeta}[L_{2}({\mathbb{R}}^{n})]\supseteq\mathcal{F}^{-1}[C_{c}^{\infty}({\mathbb{R}}^{n})]$. Since $C_{c}^{\infty}({\mathbb{R}}^{n})$ is dense in ${\mathcal{S}}({\mathbb{R}}^{n})$, then $\mathcal{F}^{-1}[C_{c}^{\infty}({\mathbb{R}}^{n})]$ also has this property. Thus, $\overline{{\rm Range\,}[{\rm O}_{M}]\cap{\mathcal{S}}({\mathbb{R}}^{n})}={\mathcal{S}}({\mathbb{R}}^{n})$. ∎

Besides the Gaussian kernel the lemma also captures a case of the Laplace kernel $\zeta({\mathbf{x}})=e^{-|{\mathbf{x}}|}$. It is well-known that the Fourier tranform of the Laplace kernel is the Poisson kernel: $\hat{\zeta}({\mathbf{x}})=\frac{c_{n}}{(1+|{\mathbf{x}}|^{2})^{\frac{n+1}{2}}}$ (which is also proper).

For $I:\mathcal{G}^{\prime}_{k}\cup{\mathcal{S}}({\mathbb{R}}^{n})\rightarrow{\mathbb{R}}^{+}$, it is natural to reduce the optimization task over tempered distributions

[TABLE]

to an optimization task over ordinary functions with a penalty term $R$,

[TABLE]

where we assume that the set of functions $\mathfrak{F}$ is rich enough to approximate weakly solutions of (29), i.e. $\overline{\mathfrak{F}}^{\ast}\supset\mathcal{G}^{\prime}_{k}$. Since we cannot guarantee that the minimum in (30) is attainable, we substitute it by infimum. For this reduction to be effective it is desirable to have the following property: if a sequence $\{f_{n}\}\subset\mathfrak{F}$ is such that $I(f_{n})+\lambda_{n}R(f_{n})-\inf\limits_{f\in\mathfrak{F}}\left(I(f)+\lambda_{n}R(f)\right)\rightarrow+0$ for $\lambda_{n}\mathop{\rightarrow}\limits^{n\rightarrow\infty}+\infty$ (i.e. $\{f_{n}\}$ solves (30) for arbitrarily large values of the regularization parameter), then there exists a growing subsequence $\{n_{k}\}$ such that $T_{f_{n_{k}}}\rightarrow^{\ast}T$ (weakly) where $T$ is a solution of (29).

We make a thorough theoretical analysis of the case $\mathfrak{F}={\mathcal{S}}({\mathbb{R}}^{n})$. If to formulate in a simplified way, for the last property to hold, the sequence ${\rm Tr}(M_{f_{n}})$ should be bounded. Details on the conditions under which this reduction holds can be found in the following subsection.

5.3 Regular solutions and reduction theorems for $\mathfrak{F}={\mathcal{S}}({\mathbb{R}}^{n})$

For a sequence $\{f_{s}\}^{\infty}_{s=1}\subseteq{\mathcal{S}}^{\prime}({\mathbb{R}}^{n})$, $\mathop{\rm Lim}\limits_{s\rightarrow\infty}f_{s}$ denotes a set of points $f\in{\mathcal{S}}^{\prime}({\mathbb{R}}^{n})$, such that there exists a growing sequence $\{s_{i}\}\subseteq{\mathbb{N}}$ and $\lim_{i\rightarrow\infty}f_{s_{i}}=f$.

For $I:\mathcal{G}^{\prime}_{k}\cup{\mathcal{S}}({\mathbb{R}}^{n})\rightarrow{\mathbb{R}}^{+}$, it is natural to reduce the optimization task (29) to an optimization task over ordinary functions with a penalty term (30). To have an equivalence between (29) and (30) we need to assume that $I$’s behaviour when approaching $f\in\mathcal{G}^{\prime}_{k}$ from a set ${\mathcal{S}}({\mathbb{R}}^{n})$ is continuous, i.e. for any sequence $\{f_{i}\}\subseteq{\mathcal{S}}({\mathbb{R}}^{n})$ such that $T_{f_{i}}\rightarrow^{\ast}f\in\mathcal{G}^{\prime}_{k}$, we have $\lim_{i\rightarrow\infty}I(T_{f_{i}})=I(f)$.

Let us introduce the notion of a regular solution both for (29) and (30). Let

[TABLE]

Definition 3**.**

Any $f\in{\rm Arg}\min\limits_{f\in\mathcal{G}^{\prime}_{k}}I(f)\bigcap\mathcal{B}_{k}$ is called a regular solution of (29).

In other words, $\mathcal{B}_{k}$ formalizes a set of distributions from $\mathcal{G}^{\prime}_{k}$, that can be approached through sequences $\{f_{i}\}\subseteq\mathcal{G}_{k}$, for which $\operatorname{Tr}(M_{f_{i}})$ does not blow up. Obviously, $\mathcal{G}_{k}\subseteq\mathcal{B}_{k}\subseteq\mathcal{G}^{\prime}_{k}$. In applications, regular solutions include all ${\rm Arg}\min\limits_{f\in\mathcal{G}^{\prime}_{k}}I(f)$ if we choose the kernel $M$ correctly. This regularity is important for a reduction to the penalty form (30), because when approaching a non-regular solution we are unable to guarantee a bounded behaviour of $M_{f}$ (and of $R(f)$).

Definition 4**.**

A sequence $\{f_{i}\}_{1}^{\infty}\subseteq{\mathcal{S}}({\mathbb{R}}^{n})$ is said to solve (30) if

[TABLE]

where $\epsilon_{i}\rightarrow+0$ and $\lambda_{i}\rightarrow+\infty,i\rightarrow+\infty$. If, additionally, $\operatorname{Tr}(M_{f_{i}})$ is bounded, then $\{f_{i}\}_{1}^{\infty}$ is said to solve (30) regularly.

Let us define

[TABLE]

Theorem 5**.**

If $M$ is a proper kernel, then ${\rm rsol\,}(I(f),R(f))\subseteq{\rm Arg}\min\limits_{f\in\mathcal{G}^{\prime}_{k}}I(f)$.

Theorem 6**.**

If $M$ is a proper kernel and ${\rm rsol\,}(I(f),R(f))\neq\emptyset$, then

[TABLE]

Theorem 7** (Reduction theorem).**

If $M$ is a proper kernel, ${\rm Arg}\min\limits_{f\in\mathcal{G}^{\prime}_{k}}I(f)\subseteq\mathcal{B}_{k}$ and ${\rm rsol\,}(I(f),R(f))\neq\emptyset$, then ${\rm rsol\,}(I(f),R(f))={\rm Arg}\min\limits_{f\in\mathcal{G}^{\prime}_{k}}I(f)$.

Suppose that we now solve a sequence of problems (30) and find $\{f_{s}\}_{1}^{\infty}$. According to Theorems 5 and 6, the following are potential scenarios:

(1) $\operatorname{Tr}(M_{f_{s}})$ blows up and the convergence is not guaranteed. This situation can be avoided by controlling $\operatorname{Tr}(M_{f})$ in an optimization process. In practice, when $f$ has a parameterized form, this can be done by bounding parameters.

If $\operatorname{Tr}(M_{f_{s}})$ does not blow up, we still have two subcases:

(2.1) $\mathop{\rm Lim}\limits_{s\rightarrow\infty}T_{f_{s}}\neq\emptyset$. This implies a positive outcome to approach (30) to the optimization problem, Problem (29).

(2.2) $\mathop{\rm Lim}\limits_{s\rightarrow\infty}T_{f_{s}}=\emptyset$. This exotic situation can happen only if a sequence $T_{f_{s}}$ leaves any sequentially compact subset of ${\mathcal{S}}^{\prime}({\mathbb{R}}^{n})$. Bounding parameters also tackles this case.

Let us now concentrate on the task (30) and describe the alternating scheme for its solution.

6 The alternating scheme

We will concentrate on problem (30). It is known [25] that the Ky Fan anti-norm is a concave function, i.e. $R(\phi)=\|M_{\phi}\|_{n-k}$ depends on $M_{\phi}$ in a concave way. It can be shown that the dependence of $R(\phi)$ on $\phi$ is both non-convex and non-concave, i.e. we deal with a non-convex optimization task.

Let $\mathcal{B}(H_{1},H_{2})$ denote a set of bounded linear operators between Hilbert spaces $H_{1}$ and $H_{2}$. For $O\in\mathcal{B}(H_{1},H_{2})$ the rank of $O$ is defined as $\dim{\mathcal{R}}(O)$. Let $L^{r}_{2}({\mathbb{R}}^{n})$ be the Hilbert space (over ${\mathbb{R}}$) of real-valued functions from $L_{2}({\mathbb{R}}^{n})$ (i.e. the real-valued $L_{2}$-space) and $L^{\ast}_{2}({\mathbb{R}}^{n})=L^{r}_{2}({\mathbb{R}}^{n})\times L^{r}_{2}({\mathbb{R}}^{n})$. The space $L^{\ast}_{2}({\mathbb{R}}^{n})$ is equivalent to $L_{2}({\mathbb{R}}^{n})$ treated as a linear space over ${\mathbb{R}}$. Below we do not distinguish $[\phi_{1},\phi_{2}]\in L^{\ast}_{2}({\mathbb{R}}^{n})$ and $\phi_{1}+{\rm i}\phi_{2}\in L_{2}({\mathbb{R}}^{n})$. It is easy to see that any $O\in\mathcal{B}(L^{\ast}_{2}({\mathbb{R}}^{n}),{\mathbb{R}}^{n})$ can be given by formula:

[TABLE]

i.e. $O\in\mathcal{B}(L^{\ast}_{2}({\mathbb{R}}^{n}),{\mathbb{R}}^{n})$ can be identified with a vector of functions $O=\begin{bmatrix}O_{i}\end{bmatrix}_{i=\overline{1,n}},O_{i}\in L_{2}({\mathbb{R}}^{n})$ and the Hilbert–Schmidt norm on $\mathcal{B}(L^{\ast}_{2}({\mathbb{R}}^{n}),{\mathbb{R}}^{n})$ (i.e. $\sqrt{{\rm Tr\,}O^{\dagger}O}$) is

[TABLE]

Recall that for a Mercer kernel $M$, ${\rm O}_{M}[\phi]({\mathbf{x}})=\int_{{\mathbb{R}}^{n}}M({\mathbf{x}},{\mathbf{y}})\phi({\mathbf{y}})d{\mathbf{y}}$ is a positive operator whose domain is ${\rm Dom}({\rm O}_{M})=\{f\in L_{2}({\mathbb{R}}^{n})\mid{\rm O}_{M}[f]\in L_{2}({\mathbb{R}}^{n})\}$ and range is a subset of $L_{2}({\mathbb{R}}^{n})$. If we assume that ${\rm Dom}({\rm O}_{M})$ is dense in $L_{2}({\mathbb{R}}^{n})$, then its adjoint ${\rm O}^{\dagger}_{M}$ and the square root $\sqrt{{\rm O}_{M}}:{\rm Dom}({\rm O}_{M})\to L_{2}({\mathbb{R}}^{n})$ can be properly defined [40]. Thus, ${\rm O}_{M}$ is self-adjoint. For any complex-valued function $f$ such that ${\rm Tr\,}M_{f}<\infty$ let us introduce a linear operator $S_{f}:L^{\ast}_{2}({\mathbb{R}}^{n})\rightarrow{\mathbb{R}}^{n}$ by the following rule:

[TABLE]

i.e. $(S_{f})_{i}=\sqrt{{\rm O}_{M}}[x_{i}f({\mathbf{x}})],i=\overline{1,n}$. In the latter definition the expression $\langle\sqrt{{\rm O}_{M}}[x_{i}f({\mathbf{x}})],\phi\rangle_{L_{2}({\mathbb{R}}^{n})}$ is finite due to $\langle\sqrt{{\rm O}_{M}}[x_{i}f({\mathbf{x}})],\sqrt{{\rm O}_{M}}[x_{i}f({\mathbf{x}})]\rangle_{L_{2}({\mathbb{R}}^{n})}=(M_{f})_{ii}<\infty$ and the Cauchy-Schwarz inequality.

Theorem 8**.**

Let $M$ be a Mercer kernel such that ${\rm Dom}({\rm O}_{M})$ is dense in $L_{2}({\mathbb{R}}^{n})$ and $\,\operatorname{Tr}M_{f}<\infty$. Then, $S_{f}\in\mathcal{B}(L^{\ast}_{2}({\mathbb{R}}^{n}),{\mathbb{R}}^{n})$ and $S_{f}S^{\dagger}_{f}=M_{f}$. Moreover,

[TABLE]

and the minimum is attained at $S=P_{f}S_{f}$ where $P_{f}=\sum_{i=1}^{k}{\mathbf{u}}_{i}{\mathbf{u}}_{i}^{\dagger}$ and $\{{\mathbf{u}}_{i}\}^{k}_{1}$ are unit eigenvectors of $M_{f}$ corresponding to the $k$ largest eigenvalues (counting multiplicities).

Proof.

The boundedness of $S_{f}$ follows from the Cauchy-Schwarz inequality:

[TABLE]

and therefore:

[TABLE]

Thus, we have checked that $S_{f}$ is bounded.

By definition, $S^{\dagger}_{f}:{\mathbb{R}}^{n}\rightarrow L^{r}_{2}({\mathbb{R}}^{n})\times L^{r}_{2}({\mathbb{R}}^{n})$ and $\langle{\mathbf{u}},S_{f}[\phi_{1},\phi_{2}]\rangle=\langle S^{\dagger}_{f}[{\mathbf{u}}],[\phi_{1},\phi_{2}]\rangle$, ${\mathbf{u}}\in{\mathbb{R}}^{n}$, $[\phi_{1},\phi_{2}]\in L^{r}_{2}({\mathbb{R}}^{n})\times L^{r}_{2}({\mathbb{R}}^{n})$. Let us denote $f_{1}={\rm Re\,}f,f_{2}={\rm Im\,}f$. It is easy to see that the following operator satisfies the latter identity:

[TABLE]

Since the adjoint is unique, then $S^{\dagger}_{f}=O$. Let us calculate $S_{f}S^{\dagger}_{f}$:

[TABLE]

Thus, $S_{f}S^{\dagger}_{f}=M_{f}$. Since $\operatorname{Tr}S_{f}S^{\dagger}_{f}<\infty$ and $\|S^{\dagger}_{f}[{\mathbf{u}}]\|^{2}\leq\langle{\mathbf{u}},M_{f}{\mathbf{u}}\rangle$, we obtain $S^{\dagger}_{f}$ is a bounded operator.

Let ${\mathbf{u}}_{1},\cdots{\mathbf{u}}_{n}$ be orthonormal eigenvectors of $\mathcal{M}_{f}=S_{f}S^{\dagger}_{f}$ and $\lambda_{1}\geq\cdots\geq\lambda_{n^{\prime}}>0$ be corresponding nonzero eigenvalues. For $\sigma_{i}=\sqrt{\lambda_{i}}$ let us define ${\mathbf{v}}_{i}=\frac{S_{f}^{\dagger}[{\mathbf{u}}_{i}]}{\sigma_{i}}$. Vector ${\mathbf{v}}_{i}$ corresponds to a pair of functions

[TABLE]

It is easy to see that ${\mathbf{v}}_{1},\cdots{\mathbf{v}}_{n^{\prime}}$ is an orthonormal basis in $\operatorname{Im}S_{f}^{\dagger}$, and $S_{f}^{\dagger}$ can be expanded in the following way:

[TABLE]

and therefore, SVD for $S_{f}$ is

[TABLE]

By the Eckart-Young-Mirsky theorem (see Theorem 4.4.7 from [41]), an optimal $S$ in $\min\limits_{S\in\mathcal{B}(L^{\ast}_{2}({\mathbb{R}}^{n}),{\mathbb{R}}^{n}),\operatorname{rank}S\leq k}\|S_{f}-S\|_{\ast}^{2}$ is defined by a truncation of SVD for $S_{f}$ at $k$th term, i.e.

[TABLE]

where $P_{f}=\sum_{i=1}^{k}{\mathbf{u}}_{i}{\mathbf{u}}_{i}^{\dagger}$ is a projection operator to first $k$ principal components of $\mathcal{M}_{f}$. Moreover, $\|S_{f}-P_{f}S_{f}\|^{2}=\sum_{i=k+1}^{n^{\prime}}\sigma^{2}_{i}=\|M_{f}\|_{n-k}=R(f)$. ∎

Given the new representation $R(f)=\min\limits_{S\in\mathcal{B}(L^{\ast}_{2}({\mathbb{R}}^{n}),{\mathbb{R}}^{n}),\operatorname{rank}S\leq k}\|S_{f}-S\|_{\ast}^{2}$ we have

[TABLE]

Thus, it is natural to view the Task (30) as a minimization of $I(\phi)+\lambda\|S_{\phi}-S\|_{\ast}^{2}$ over two objects: $f\in\mathfrak{F}$ and $S\in\mathcal{B}(L^{\ast}_{2}({\mathbb{R}}^{n}),{\mathbb{R}}^{n}):\operatorname{rank}S\leq k$. The simplest approach to minimize a function over two arguments is to optimize alternatingly, i.e. first over $f$, and then over $S:\operatorname{rank}S\leq k$, and so on. Theorem 8 gives that the minimization over $S$ is equivalent to the truncation of ${\rm SVD}(S_{f})$ at the $k$-th term. This idea, that we dub the alternating scheme (AS), is described in Algorithm 1.

The alternating algorithm 1 allows for a reformulation in the dual space. By this we mean that in Algorithm 1 we substitute $\widehat{\phi}_{t}$ for the original $\phi_{t}$. If the primal Algorithm 1 deals with operators $S_{\phi},S_{\phi_{t-1}}$, the dual version deals with vectors of functions $\sqrt{{\widehat{G}_{\sigma}}}\frac{\partial{\widehat{\phi}}}{\partial{\mathbf{x}}},\sqrt{{\widehat{G}_{\sigma}}}\frac{\partial{\widehat{\phi}}_{t-1}}{\partial{\mathbf{x}}}$. Details of the dual algorithm can be found in E.

The objective $I(f)+\lambda R(f)$ can have many local minima due to the effect of the penalty term $R(f)$. Therefore, the Alternating Scheme 1 is strongly dependant on the initialization step. One of such initialization procedures for the task (17) is described in the next section.

7 An approximate algorithm for the MMD-PCA

Let us analyze the task (17) in the case where $K({\mathbf{x}},{\mathbf{y}})=({\mathbf{x}}\cdot{\mathbf{y}})H({\mathbf{x}},{\mathbf{y}})$ and $H$ is a Mercer kernel (by construction, $K$ is also a Mercer kernel). In this section we demonstrate that, given a distribution $f({\mathbf{x}})=\frac{1}{N}\sum_{i=1}^{N}\delta^{n}({\mathbf{x}}-{\mathbf{x}}_{i})$, a good guess for a $k$-dimensional space in which an optimal solution is supported is a span of the first $k$ principal components of $H_{f}$ (see the Algorithm 2).

A specifics of this type of kernels is that the MMD distance (induced by $K$) till an optimal solution of the MMD-PCA task is bounded below by the Ky Fan $k$-antinorm of $H_{f}$, as shown in the following theorem.

Theorem 9**.**

Let $K({\mathbf{x}},{\mathbf{y}})=({\mathbf{x}}\cdot{\mathbf{y}})H({\mathbf{x}},{\mathbf{y}})$ where $H:{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}\to{\mathbb{R}}$ is a Mercer kernel, ${\rm Dom}({\rm O}_{K})$ is dense in $L_{2}({\mathbb{R}}^{n})$ and $f$ is such that ${\rm Tr\,}H_{f}<\infty$. Then,

[TABLE]

where $\|g\|_{K}^{2}=\langle g|K|g\rangle$, $\lambda_{1}\geq\cdots\geq\lambda_{n}$ are eigenvalues (counting multiplicities) of $H_{f}$.

Sketch..

Let us apply Theorem 8 to the kernel $H$ and the function $f$. Recall that an element $O\in\mathcal{B}(L^{\ast}_{2}({\mathbb{R}}^{n}),{\mathbb{R}}^{n})$ can be identified with a vector-function $[O_{i}]_{i=1}^{n},O_{i}\in L_{2}({\mathbb{R}}^{n})$ where $O[\phi]_{i}={\rm Re\,}\langle O_{i},\phi\rangle_{L_{2}({\mathbb{R}}^{n})}$. Since $S_{f}$ corresponds to $[\sqrt{{\rm O}_{H}}[x_{i}f({\mathbf{x}})]]_{i=1}^{n}$, the representation of Theorem 8 gives us

[TABLE]

The restriction $\operatorname{rank}S\leq k$ is equivalent to ${\rm dim\,}\mathcal{L}_{S}\leq k$ where $\mathcal{L}_{S}=\{\sum_{i=1}^{n}\xi_{i}S_{i}({\mathbf{x}})\mid\xi_{i}\in{\mathbb{R}}\}$ and $S$ corresponds to $[S_{i}]_{i=1}^{n}$.

Let $f^{\prime}\in\mathcal{G}_{k}$. By Theorem 3 we have ${\rm dim\,}{\rm span}_{\mathbb{R}}(\{x_{1}f^{\prime},\cdots,x_{n}f^{\prime}\})\leq k$. By Theorem 4, $\langle x_{i}f^{\prime}|H|x_{i}f^{\prime}\rangle$ is finite, therefore $\sqrt{{\rm O}_{H}}[x_{i}f^{\prime}({\mathbf{x}})]$ can be properly defined and is in $L_{2}({\mathbb{R}}^{n})$. Therefore,

[TABLE]

For any $f^{\prime}\in\mathcal{G}_{k}$ one can set $S_{i}=\sqrt{{\rm O}_{H}}[x_{i}f^{\prime}({\mathbf{x}})]$ and search over all possible $f^{\prime}\in\mathcal{G}_{k}$ in the minimization operator. Thus,

[TABLE]

∎

Let $\mu_{\rm data}$ be a uniform distribution over $\{{\mathbf{x}}_{i}\}_{i=1}^{N}$ and $d_{\rm MMD}$ be the MMD distance induced by $K({\mathbf{x}},{\mathbf{y}})=({\mathbf{x}}\cdot{\mathbf{y}})H({\mathbf{x}},{\mathbf{y}})$. The last theorem can be applied to a smoothed empirical distribution $f_{\varepsilon}({\mathbf{x}})=\frac{1}{N}\sum_{i=1}^{N}G^{n}_{\varepsilon}({\mathbf{x}}-{\mathbf{x}}_{i})$ and then, we can send $\varepsilon\to 0$. All the more, the inequality will be satisfied if we search over $\mu\in\mathcal{P}_{k}$ due to $T_{\mu}\in\mathcal{G}_{k}$. Thus,

[TABLE]

where $\lambda_{1}\geq\cdots\geq\lambda_{n}$ are eigenvalues (counting multiplicities) of $H_{f}$, $f({\mathbf{x}})=\frac{1}{N}\sum_{i=1}^{N}\delta^{n}({\mathbf{x}}-{\mathbf{x}}_{i})$. Thus, $(\sum_{i=k+1}^{n}\lambda_{i})^{1/2}$ is a lower bound of the solution of (17).

For such an important practical case as the HM-MMD-PCA, a multiple of the square root of the Ky Fan $k$-antinorm of $H_{f}$ is also an upper bound.

Theorem 10**.**

Let $H({\mathbf{x}},{\mathbf{y}})=P({\mathbf{x}}\cdot{\mathbf{y}})$ where $P(x)=c_{0}+c_{1}x+\cdots+c_{l-1}x^{l-1}$, $c_{i}\geq 0,i\in[l-1]$, $f({\mathbf{x}})=\frac{1}{N}\sum_{i=1}^{N}\delta^{n}({\mathbf{x}}-{\mathbf{x}}_{i})=T_{\mu_{\rm data}}$ and $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{n}$ are eigenvalues of $H_{f}$. Then,

[TABLE]

The following corollary is straightforward from the last theorem.

Corollary 1**.**

A $2$-approximating solution of the task (20) can be efficiently found by the Algorithm 2.

For the case when $H$ is the Gaussian kernel the situation is slightly trickier.

Theorem 11**.**

Let $H({\mathbf{x}},{\mathbf{y}})=e^{-\frac{\sigma^{2}\|{\mathbf{x}}-{\mathbf{y}}\|^{2}}{2}}$, $m_{i}=\int_{{\mathbb{R}}^{n}}\|{\mathbf{x}}\|^{i}|f({\mathbf{x}})|d{\mathbf{x}}<\infty,i\in[4]$ and $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{n}$ are eigenvalues of $H_{f}$. Then,

[TABLE]

where $M={\mathcal{O}}(m_{2}+\sqrt{2}\sigma\sqrt{m_{4}m_{2}}+\sqrt{2}\sigma m_{3})$.

An analogous theorem can be proved for $H({\mathbf{x}},{\mathbf{y}})=(1+\sigma^{2}\|{\mathbf{x}}-{\mathbf{y}}\|^{2})^{-\frac{n+1}{2}}$, i.e. the Poisson kernel.

8 Experiments

The alternating scheme 1 is a general optimization method that needs to be specified for every optimization task. We designed numerical specifications of the alternating scheme 1 for all 4 optimization tasks: (17), (20), (22) and (23) and made experiments with all of them. Details of the algorithms, i.e. numerical methods to minimize over $\phi$ and calculate $M_{\phi_{t}}$, can be found in Appendix. Note that for WD-PCA (22) we exploit the alternating scheme in the initial form (i.e. 1), and for MMD-PCA (17), HM-MMD-PCA (20) and SDR-ORF (23) we use the dual version of the scheme.

Behaviour of the Gaussian MMD-PCA for small $h$. We studied the difference in the behavior of PCA and a solution of (17), for the distance function induced by the kernel $K({\mathbf{x}},{\mathbf{y}})=\frac{1}{\sqrt{(8\pi h^{2})^{n}}}e^{-\frac{\|{\mathbf{x}}\|^{2}}{8h^{2}}}$, obtained by the alternating scheme 1 (AS for MMD-PCA), for the case when $h$ is small compared to the standard deviation of features. Experiments show that they are sharply different when data points are sampled along a low-dimensional manifold ${\mathfrak{M}}$, which is bent globally, goes through the origin $O$ and has a large curvature at $O$ (see Fig. 1(a)). Since generated points do not lie on an affine subspace, the global nature of PCA makes it hard to interprete principal directions.

We select a smooth function $f:{\mathbb{R}}^{n-1}\rightarrow{\mathbb{R}}$, such that $f({\mathbf{0}})=0$ and generate points in the following way: points ${\mathbf{x}}_{1},{\mathbf{x}}_{2},\cdots,{\mathbf{x}}_{N}\sim[-10,10]^{n-1}$ are sampled uniformly, after calculation of $y_{i}=f({\mathbf{x}}_{i})$ we add some noise: ${\mathbf{z}}_{i}=({\mathbf{x}}_{i},y_{i})+\boldsymbol{\epsilon}_{i},\boldsymbol{\epsilon}_{i}\sim{\mathcal{N}}(0,0.01I_{n})$. Both PCA and MMD-PCA are applied to the dataset (first 3 pictures on Figure 1(a)). As we see, MMD-PCA, unlike PCA, tries to catch ideal alignments of points rather that searching for a global alignment of points (which is non-existent). This property of MMD-PCA makes it a promising tool for a calculation of the tangent space to a data manifold at a given point. Fourth picture shows that when we have 2 equally important directions in data such that the first principal direction of PCA is between them (red line), and we set $k=1$, then MMD-PCA (green line) always chooses one of those directions. These experimental results are aligned with the theoretical observation given in Example 1, in which we show that the Gaussian MMD-PCA task for $h\to 0+$ is equivalent to finding a $k$-dimensional subspace that contains as many points of a dataset as possible. Thus, the Gaussian MMD-PCA can be considered as a method that can be potentially used to tackle the latter NP-hard problem. Some informal discussion of this problem can be found in [42].

Experiments with outlier detection (MMD-PCA, HM-MMD-PCA, WD-PCA). Following the experiment setup of [37], we choose parameters $N=n=400,\updelta=0.05$ (or $0.1$), $k=10$ and generate random matrices $A\in{\mathbb{R}}^{N(1-\updelta)\times k},B\in{\mathbb{R}}^{n\times k}$ whose entries are iid as ${\mathcal{N}}(0,1)$. Then, columns of the matrix $BA^{T}\in{\mathbb{R}}^{n\times N(1-\updelta)}$ (whose rank is $\leq k$) are concatenated with columns of the matrix $C\in{\mathbb{R}}^{n\times N\updelta}$: $X={\rm concat}(BA^{T},C)\in{\mathbb{R}}^{n\times N}$. The entries in $C$ are either iid as ${\mathcal{N}}(0,1)$ (case I) or $N\updelta$ copies of the same vector whose entries are iid as ${\mathcal{N}}(0,1)$ (case II). Let $X=[{\mathbf{x}}_{1},\cdots,{\mathbf{x}}_{N}]$, i.e. columns of $X$ are the data points. Thus, $N(1-\updelta)$ columns of $BA^{T}$ lie in a $k$-dimensional subspace of ${\mathbb{R}}^{n}$ and $N\updelta$ columns of $C$ are outliers, and solutions of tasks (17), (20) or (22) for this dataset are expected to be supported in a column space of $BA^{T}$.

After every iteration (step $t$ of the alternating scheme 1) we calculate the Frobenius distance between the projection operator $P_{t}$ of 1 and the projection operator $P$ to the column space of $BA^{T}$, i.e. $\|P_{t}-P\|_{F}$. For the task (22), the dependence of $\|P_{t}-P\|_{F}$ on $t$ for different values of parameters $\updelta$ and $\lambda$ is shown in Figure 1(b). For tasks (17), (20) the behaviour of the alternating scheme is similar, 7 iterations are enough to approach the optimal subspace.

One of main goals of this experimental setup was to study how the kernel $M$, that defines the regularizer $R(f)$ by equation(28), affects the quality of a solution. Besides the speed of convergence we were interested in how $\|P^{\ast}-P\|_{F}$, where $P^{\ast}=\lim_{t\rightarrow\infty}P_{t}$ is the final projection operator (e.g. $P_{20}$ in practice), depends on the parameter ${\sigma}$ of the kernel $M=G^{n}_{\sigma}$ (bandwidth). It is natural to expect the quality of the solution $P^{\ast}$ to degrade as ${\sigma}\rightarrow+\infty$ (this corresponds to $M({\mathbf{x}},{\mathbf{y}})\rightarrow 0$), and, less trivially, as ${\sigma}\rightarrow 0$ (this corresponds to $M({\mathbf{x}},{\mathbf{y}})\rightarrow\delta^{n}({\mathbf{x}}-{\mathbf{y}})$). As the right plot on Figure 1(b) shows, for the HM-MMD-PCA, the solution $P^{\ast}$ is close to the correct $P$ if the bandwidth $\sigma$ is in interval $[e^{-2},e^{3}]$ and it degrades beyond that interval. For the Gaussian MMD-PCA the degrading occurs beyond $[e^{1.3},e^{3}]$. For the WD-PCA the interval for $\sigma$ is sligtly narrower than $[e^{1.3},e^{3}]$. Our numerical specification of the alternating scheme for WD-PCA involves training regularized Generative Adversarial Network (see for details I) and are based on numerically unstable algorithms for the Wasserstein distance minimization. Finding numerical specifications for WD-PCA with a more stable behavior is a future work.

Experiments with SDR-ORF. We made experiments on the standard datasets, Heart, Breast Cancer, Ionosphere, Diabetes, Boston House Prices and Wine Quality. First, we applied the Sliced Inverse Regression algorithm (SIR) [13] to the training set and calculated the effective subspace for $k=2,3$. All points were projected onto that space and we obtained two- or three-dimensional representations of input points. In the last step we applied the ten nearest neighbors algorithm (KNN) to predict outputs (based on reduced inputs) on the test set (for the regression case, the 10-KNN regression was used). The same scheme was repeated with PCA, Kernel Dimensionality Reduction (KDR) algorithm [18] and the alternating scheme 1 (AS) adapted for the SDR-ORF.

We experimented with the dual version of algorithm 1, setting (after the data was standardized) the kernel’s parameter ${\sigma}=0.8$222Since the role of the parameter ${\sigma}$ is similar to that of the bandwidth in the kernel density estimation, we use Silverman’s rule of thumb to set ${\sigma}=N^{-1/(n+4)}$. and $\lambda=10.0$. Details of its numerical implementation can be found in J. In the table 1 one can see the obtained test set accuracy on the classification tasks and R2 on the regression tasks. As we see from the table 1, after reducing the dimension of an input to $k=2,3$, we are still able to obtain good accuracy of prediction on a test set and the AS for the SDR-ORF is competitive in comparison with other methods. Note that all listed datasets are of moderate size and our Python scripts managed to compute an effective subspace in 3-5 minutes on a PC with GTX Titan X (Pascal), Intel Core i7-7700K (4.20 GHz), 64 GB RAM.

Experiments with the shadow/black removal. We made experiments with Yale B dataset [43], which is a popular benchmark for testing robust versions of PCA. That dataset contains images of 28 human subjects under 9 poses and 64 illumination conditions. Test images used in the experiments are cropped and re-sized to 168x192 images, making the dimensionality of every image 32256. Thus, each human subject corresponds to a collection of 32256-dimensional vectors that lie on some low-dimensional subspace $\mathcal{L}$ of ${\mathbb{R}}^{32256}$. We search for this subspace, assuming that its dimension is either 1 or 5, using PCA and the Algorithm 2 with the kernel $K({\mathbf{x}},{\mathbf{y}})=({\mathbf{x}}\cdot{\mathbf{y}})e^{-\frac{\|{\mathbf{x}}-{\mathbf{y}}\|^{2}}{n}}$ (which we simply call Gauss). Our experiments showed that behaviour of PCA and Gauss are quite similar if the dimension of $\mathcal{L}$ is 5, though Gauss removes more shadows and preserves more details of an original image if the dimension of $\mathcal{L}$ is 1 (see Figures 2 and 3). A processing of each human subject by Gauss takes seconds on Google Colab.

Experiments with the background modeling. For these experiments we use the dataset for testing background estimation algorithms SBMnet [44]. The dataset contains frames of short videos and the frame of a background for each video (so called the ground truth). Spatial resolutions of the videos vary from 240x240 to 800x600. Thus, a collection of frames of every video is a set of high-dimensional vectors (with a dimension up to 480000) that, again, lie on a low dimensional subspace $\mathcal{L}$. We assume that the dimension of $\mathcal{L}$ is 5. We recover $\mathcal{L}$ using PCA and the Algorithm 2 for kernels $K({\mathbf{x}},{\mathbf{y}})=\sum_{i=1}^{4}(\frac{{\mathbf{x}}\cdot{\mathbf{y}}}{n})^{i}$, $K({\mathbf{x}},{\mathbf{y}})=({\mathbf{x}}\cdot{\mathbf{y}})e^{-\frac{a\|{\mathbf{x}}-{\mathbf{y}}\|^{2}}{n}}$, $K({\mathbf{x}},{\mathbf{y}})=({\mathbf{x}}\cdot{\mathbf{y}})e^{-\frac{a\|{\mathbf{x}}-{\mathbf{y}}\|}{\sqrt{n}}}$ and $K({\mathbf{x}},{\mathbf{y}})=({\mathbf{x}}\cdot{\mathbf{y}})(1+\frac{a\|{\mathbf{x}}-{\mathbf{y}}\|^{2}}{n})^{-\frac{n+1}{2}}$ (which we simply call Kurtosis, Gauss, Laplace and Poisson respectively). Recall that, according to corollary 1, this algorithm is 2-approximating for Kurtosis. Subsequently, we compute the median of the vectors, projected onto $\mathcal{L}$, and define the latter to be the recovered background image (see Figure 4). Measures of consistency with the ground truth backgrounds are then calculated using Python scripts downloaded from [45]. Six measures are used: AGE (average of the gray-level absolute difference between a ground truth image and a computed background image), pEPs (percentage of pixels in a computed background whose value differs from the value of the corresponding pixel in a ground truth by more than a threshold), pCEPS (percentage of pixels whose 4-connected neighbors are also error pixels), MSSSIM (estimate of the perceived visual distortion), PSNR (Peak-Signal-to-Noise-Ratio, or $10\log_{10}(\frac{(L-1)^{2}}{MSE})$ where $L$ is the maximum number of grey levels and MSE is the mean squared error between a ground truth and a computed background images), CQM (Color image Quality Measure). Codes that compute listed metrics can be found in [45]. As shown on Table 6, experiments again demonstrated very similar behavior of PCA, Kurtosis, Gauss, Laplace and Poisson with very close accuracies of the background reconstruction. Best measures of consistency with the ground truth images were achieved for Gauss ($a=5.0,25.0$) and Laplace ($a=5.0$). For a comparison with other methods, we also give accuracies of other methods based on low-rank approximation and an accuracy of a state-of-the-art method that was specifically tailored for that task [46]. On figure 5 one can see that background images computed by PCA and Gauss are almost identical, though Gauss is less likely than PCA to add local artefacts, such as blurs, noise etc.

The processing of the whole SBMnet dataset using PCA/Kurtosis/Gauss/ Laplace takes approximately the same time — 25 minutes on a cluster equipped with Intel Xeon Platinum 8168 Processors (33M Cache, 2.70 GHz) and 1TB RAM. The code is available on github to facilitate the reproducibility of our results.

Scalability of algorithms. A major practical limitation of the alternating scheme 1 comes from the fact that it involves an optimization over a set of functions $\mathfrak{F}$, which in applications is either a feedforward neural network (as in our specifications of the AS for SDR-ORF, MMD-PCA, HM-MMD-PCA) or a generative neural network (WD-PCA). A speed of optimization is also strongly dependant on the objective’s landscape. Thus, for large scale datasets, with a dimension of vectors $\gg 10^{3}$, and a sophisticated structure of a regression function (SDR-ORF) or a data distribution (MMD-PCA, HM-MMD-PCA, WD-PCA), the alternating scheme is substantially slower in comparison with other popular methods (such as PCA for the UDR, or SIR/KDR for the SDR).

In the special case of MMD-PCA (that includes HM-MMD-PCA), the approximate algorithm 2 can be used as a surrogate of the alternating scheme. It requires the same time as PCA and can be applied to datasets with a dimension of vectors $\sim 10^{6}-10^{7}$. Also, the Algorithm 2 can be used for an initialization of the alternating scheme.

9 Conclusions

We develop a new optimization framework for LDR problems. The alternating scheme for the optimization task demonstrates both the computational efficiency and the applicability to real-world data. The algorithm performs quite stably when we vary most of the hyperparameters, though it crucially depends on two parameters, the bandwidth of the “smoothing” kernel $M$, $\sigma$, and the penalty parameter $\lambda$. We believe that the MMD-PCA/HM-MMD-PCA/WD-PCA methods for UDR could be used as an alternative to PCA in study fields in which data demonstrate “heavy-tailed” and “non-Gaussian” behavior, such as financial applications or computer vision. Also, our formulation of SDR-ORF is free from any assumptions on the distribution of input-output pairs, which makes it an alternative to other methods of efficient subspace estimation. A more detailed report on these topics is a subject of future research.

Appendix A Proofs for section 3

A.1 Proof of Theorem 1: given for completeness

Proof.

The inclusion $\mathcal{G}^{\prime}_{k}\subseteq\overline{\mathcal{G}_{k}}^{\ast}$ follows from a well-known fact that ${\mathcal{S}}({\mathbb{R}}^{k})$ is dense in ${\mathcal{S}}^{\prime}({\mathbb{R}}^{k})$. I.e. for any $f\in{\mathcal{S}}^{\prime}({\mathbb{R}}^{k})$ one can always find a sequence $\{f_{i}\}\subseteq{\mathcal{S}}({\mathbb{R}}^{k})$ such that $T_{f_{i}}\rightarrow^{\ast}f$. Therefore, for any $(f\otimes\delta^{n-k})_{U}\in\mathcal{G}^{\prime}_{k}$ there is a sequence $\{(T_{f_{i}}\otimes\delta^{n-k})_{U}\}\subseteq\mathcal{G}_{k}$ such that $(T_{f_{i}}\otimes\delta^{n-k})_{U}\rightarrow^{\ast}(f\otimes\delta^{n-k})_{U}$. Thus, $\mathcal{G}^{\prime}_{k}\subseteq\overline{\mathcal{G}_{k}}^{\ast}$.

Since $\mathcal{G}_{k}\subseteq\mathcal{G}^{\prime}_{k}$, to prove $\mathcal{G}^{\prime}_{k}=\overline{\mathcal{G}_{k}}^{\ast}$ it is enough to show that $\mathcal{G}^{\prime}_{k}$ is sequentially closed.

We need a simple fact from a theory of distributions.

Lemma 2**.**

If $T_{i}\rightarrow^{\ast}T$ and $\phi_{i}\rightarrow\phi$, then $\langle T_{i},\phi_{i}\rangle\rightarrow\langle T,\phi\rangle$.

Proof of Lemma.

Schwartz space ${\mathcal{S}}({\mathbb{R}}^{n})$ is a Fréchet space, therefore the Banach-Steinhaus theorem applies to ${\mathcal{S}}^{\prime}({\mathbb{R}}^{n})$. Since $T_{i}\rightarrow^{\ast}T$, we have $\sup_{i}|\langle T_{i},\phi\rangle|<\infty$ for any $\phi\in{\mathcal{S}}({\mathbb{R}}^{n})$. From the Banach-Steinhaus theorem, applied to a set $\{T_{i}\}_{1}^{\infty}$, we obtain for any $\epsilon>0$, there is a neighbourhood $U$ of ${\mathbf{0}}\in{\mathcal{S}}({\mathbb{R}}^{n})$ such that $|\langle T_{i},\phi\rangle|<\epsilon$ whenever $\phi\in U$. Thus, $|\langle T_{i},\phi_{i}-\phi\rangle|<\epsilon$ for a large enough $i$. From that we conclude that $\langle T_{i},\phi_{i}\rangle\rightarrow\langle T,\phi\rangle$. ∎

For any $T\in{\mathcal{S}}^{\prime}({\mathbb{R}}^{n})$ and $\psi\in{\mathcal{S}}({\mathbb{R}}^{n-k})$, let us define $T^{\psi}\in{\mathcal{S}}^{\prime}({\mathbb{R}}^{k})$ as $\langle T^{\psi},\phi\rangle=\langle T,\phi\otimes\psi\rangle$.

Suppose that $\{f_{i}\}_{1}^{\infty}\subseteq{\mathcal{S}}^{\prime}({\mathbb{R}}^{k})$, $\{U_{i}\}_{1}^{\infty}$ are such that $(f_{i}\otimes\delta^{n-k})_{U_{i}}\rightarrow^{\ast}f$. We need to prove that $f\in\mathcal{G}^{\prime}_{k}$. Since a set of orthogonal matrices is compact, then one can always find a subsequence $\{U_{n_{i}}\}$ such that $U_{n_{i}}\rightarrow U$. Since $(f_{n_{i}}\otimes\delta^{n-k})_{U_{n_{i}}}\rightarrow^{\ast}f$ and $\phi(U_{n_{i}}{\mathbf{x}})\rightarrow\phi(U{\mathbf{x}})$ (for any fixed $\phi\in{\mathcal{S}}({\mathbb{R}}^{n})$), using lemma 2 we obtain:

[TABLE]

Thus, we have $f_{n_{i}}\otimes\delta^{n-k}\rightarrow^{\ast}f_{U^{T}}$. From the last we see that $f_{n_{i}}\rightarrow^{\ast}f^{\psi}_{U^{T}}$ where $\psi$ is such that $\psi({\mathbf{0}})=1$. Therefore, $f_{U^{T}}=f^{\psi}_{U^{T}}\otimes\delta^{n-k}$ and $f=(f^{\psi}_{U^{T}}\otimes\delta^{n-k})_{U}\in\mathcal{G}^{\prime}_{k}$. ∎

A.2 Proof of Theorem 3

Proof of Theorem 3 ($\Rightarrow$).

Let us prove that from $T=(f\otimes\delta^{n-k})_{U}$, $f\in{\mathcal{S}}^{\prime}({\mathbb{R}}^{k})$, $U^{T}U=I_{n}$ it follows that $\dim{\rm span}_{\mathbb{R}}\{x_{1}T,x_{2}T,\cdots,x_{n}T\}\leq k$.

It is easy to see that $x_{i}[f\otimes\delta^{n-k}]=0$ if $i>k$. If $U=[{\mathbf{u}}_{1},\cdots,{\mathbf{u}}_{n}]^{T}$, then for $i>k$ we have $0=(x_{i}[f\otimes\delta^{n-k}])_{U}={\mathbf{u}}^{T}_{i}{\mathbf{x}}(f\otimes\delta^{n-k})_{U}={\mathbf{u}}^{T}_{i}{\mathbf{x}}T$.

Thus, we have $n-k$ orthogonal vectors, ${\mathbf{u}}_{k+1},\cdots,{\mathbf{u}}_{n}$, such that

[TABLE]

Using standard linear algebra we obtain there are at most $k^{\prime}$ distributions $x_{i_{1}}T,\cdots,x_{i_{k^{\prime}}}T,k^{\prime}\leq k$ that form a basis of ${\rm span}_{{\mathbb{R}}}\{x_{i}T\}^{n}_{1}$. ∎

To prove the second part of theorem we need the following lemma.

Lemma 3**.**

If $T\in{\mathcal{S}}^{\prime}({\mathbb{R}}^{n})$ is such that $y_{i}T=0$ for any $i>k$, then $T\in\mathcal{G}^{\prime}_{k}$.

Proof of lemma.

Recall from functional analysis, for $f\in{\mathcal{S}}^{\prime}({\mathbb{R}}^{n})$, the tempered distribution $\frac{\partial f}{\partial x_{i}}$ is defined by the condition $\langle\frac{\partial f}{\partial x_{i}},\phi\rangle=-\langle f,\frac{\partial\phi}{\partial x_{i}}\rangle$. Once the Fourier transform is applied, our lemma’s dual version is equivalent to the following formulation: if $\frac{\partial f}{\partial x_{i}}=0,i>k$, then $f\in\overline{\mathcal{F}_{k}}^{\ast}$. Let us prove it in this formulation.

Recall that a set of infinitely differentiable functions with a compact support is denoted by $C^{\infty}_{c}({\mathbb{R}})$. Suppose $\phi\in{\mathcal{S}}({\mathbb{R}}^{n})$ and $p\in C^{\infty}_{c}({\mathbb{R}})$ are chosen in such a way that $\int^{\infty}_{-\infty}p(y_{i})dy_{i}=1$, $\textbf{supp\,}p\subseteq[A,B]$. Let us define:

[TABLE]

It is easy to see that for any $\alpha\in{\mathbb{N}}^{n-1},\alpha^{\prime}\in{\mathbb{N}},\beta\in{\mathbb{N}}^{n-1},\beta^{\prime}\in{\mathbb{N}}$ we have (at least one derivative over $x_{i}$ is present):

[TABLE]

The terms ${\mathbf{x}}^{\alpha}_{-i}x^{\alpha^{\prime}}_{i}\frac{\partial^{\beta,\beta^{\prime}}\phi({\mathbf{x}})}{\partial{\mathbf{x}}^{\beta}_{-i}\partial x^{\beta^{\prime}}_{i}}$ and $x^{\alpha^{\prime}}_{i}\frac{\partial^{\beta^{\prime}}p(x_{i})}{\partial x^{\beta^{\prime}}_{i}}$ are bounded by the definition of ${\mathcal{S}}({\mathbb{R}}^{n}),C^{\infty}_{c}({\mathbb{R}})$. The boundedness of $\int^{\infty}_{-\infty}{\mathbf{x}}^{\alpha}_{-i}\frac{\partial^{\beta}\phi({\mathbf{x}}_{-i},y_{i})}{\partial{\mathbf{x}}^{\beta}_{-i}}dy_{i}$ is a consequence of the inequality $|{\mathbf{x}}^{\alpha}_{-i}\frac{\partial^{\beta}\phi({\mathbf{x}}_{-i},y_{i})}{\partial{\mathbf{x}}^{\beta}_{-i}}|\leq\frac{C}{1+y_{i}^{2}}$ (which holds because $\phi\in{\mathcal{S}}({\mathbb{R}}^{n})$).

Analogously (when not a single derivative over $x_{i}$ is present):

[TABLE]

The second term is 0 when $x_{i}\leq A$. It is also bounded when $x_{i}>A$ because $|{\mathbf{x}}^{\alpha}_{-i}\frac{\partial^{\beta}\phi({\mathbf{x}}_{-i},y_{i})}{\partial{\mathbf{x}}^{\beta}_{-i}}|\leq\frac{C^{\prime}}{(1+y^{2}_{i})^{\alpha^{\prime}+1}}$. Therefore,

[TABLE]

The latter is bounded because $\lim_{x_{i}\rightarrow+\infty}|x_{i}|^{\alpha^{\prime}}\int^{\infty}_{x_{i}}\frac{C^{\prime}}{(1+y^{2}_{i})^{\alpha^{\prime}+1}}dy_{i}=0$.

The first term is 0 when $x_{i}\geq B$ and for $x_{i}<B$ it satisfies:

[TABLE]

The latter is also bounded, since $\lim_{x_{i}\rightarrow-\infty}|x_{i}|^{\alpha^{\prime}}\int^{x_{i}}_{-\infty}\frac{C^{\prime}}{(1+y^{2}_{i})^{\alpha^{\prime}+1}}dy_{i}=0$.

Thus, ${\mathbf{x}}^{\alpha}\frac{\partial^{\beta}r({\mathbf{x}})}{\partial{\mathbf{x}}^{\beta}}$ is bounded and $r\in{\mathcal{S}}({\mathbb{R}}^{n})$. Therefore, $\frac{\partial f}{\partial x_{i}}=0$ implies:

[TABLE]

Since this sequence of arguments works for any $i>k$, we can apply them sequentially to initial $\phi\in{\mathcal{S}}({\mathbb{R}}^{n})$ w.r.t. $x_{k+1},...,x_{n}$. Thus, for any $p_{k+1},...,p_{n}\in C_{c}({\mathbb{R}})$ such that $\int^{\infty}_{-\infty}p_{i}(y_{i})dy_{i}=1$ we obtain:

[TABLE]

Moreover, since $C^{\infty}_{c}({\mathbb{R}})$ is dense in ${\mathcal{S}}({\mathbb{R}})$, we can assume that $p_{k+1},...,p_{n}\in{\mathcal{S}}({\mathbb{R}})$. For the inverse Fourier transform $T=\mathcal{F}^{-1}[f]$ the latter condition becomes equivalent to:

[TABLE]

for any $p^{\prime}_{k+1},...,p^{\prime}_{n}\in{\mathcal{S}}({\mathbb{R}})$ such that $p^{\prime}_{i}(0)=1$. Let us define $p^{\prime}_{i}(x_{i})=e^{-x^{2}_{i}}$. It is easy to check that $T=g\otimes\delta^{n-k}$ where $g\in{\mathcal{S}}^{\prime}({\mathbb{R}}^{k}),\langle g,\psi\rangle=\langle T,e^{-|{\mathbf{x}}_{k+1:n}|^{2}}\psi({\mathbf{x}}_{1:k})\rangle$ for $\psi\in{\mathcal{S}}({\mathbb{R}}^{k})$. Thus, $T\in\mathcal{G}^{\prime}_{k}$ and lemma is proved. ∎

Proof of Theorem 3 ($\Leftarrow$).

If $\dim{\rm span}_{\mathbb{R}}\{x_{1}T,x_{2}T,\cdots,x_{n}T\}\leq k$, then

[TABLE]

Thus, there exist at least $n-k$ orthonormal vectors ${\mathbf{v}}_{k+1},\cdots,{\mathbf{v}}_{n}$, such that $[x_{1}T,\cdots,x_{n}T]{\mathbf{v}}_{i}=0$. Therefore, $[x_{1}T,\cdots,x_{n}T]{\mathbf{v}}_{i}=({\mathbf{v}}^{T}_{i}{\mathbf{x}})T=0$.

Let us complete ${\mathbf{v}}_{k+1},\cdots,{\mathbf{v}}_{n}$ to form an orthonormal basis of ${\mathbb{R}}^{n}$: ${\mathbf{v}}_{1},\cdots,{\mathbf{v}}_{n}$. Let us define a matrix $V=[{\mathbf{v}}_{1},\cdots,{\mathbf{v}}_{n}]$. It is easy to see that:

[TABLE]

Since for $i>k$ we have $({\mathbf{v}}^{T}_{i}{\mathbf{x}})T=0$, then $x_{i}T_{V}=0$. Using lemma 3 we obtain $T_{V}\in\mathcal{G}^{\prime}_{k}$. Therefore, $(T_{V})_{V^{T}}=T\in\mathcal{G}^{\prime}_{k}$. Theorem proved. ∎

Appendix B Structure of WD-PCA

Recall that $({\mathbb{R}}^{n},\|\cdot\|)$ is a Banach space and $p\geq 1$. Now, let us consider an optimization problem: for a given $X\in{\mathbb{R}}^{n\times N}$ solve

[TABLE]

where $\|\cdot\|_{p}$ is a norm on ${\mathbb{R}}^{n\times N}$ defined by $\|[{\mathbf{s}}_{1},\cdots,{\mathbf{s}}_{N}]\|\overset{def}{=}(\sum_{i=1}^{N}\|{\mathbf{s}}_{i}\|^{p})^{1/p}$.

The following simple theorem shows that the two tasks are connected so that the solution of one directly leads to another’s solution.

Theorem 12**.**

Given data points $\{{\mathbf{x}}_{1},\cdots,{\mathbf{x}}_{N}\}$, let $X=[{\mathbf{x}}_{1},\cdots,{\mathbf{x}}_{N}]\in{\mathbb{R}}^{n\times N}$. Then,

[TABLE]

Moreover, $\min_{\nu\in{\mathcal{P}}_{k}}W_{p}(\mu_{\rm{data}},\nu)$ is attained on $\nu^{\ast}$, where $\nu^{\ast}$ is a uniform distribution over $\{{\mathbf{y}}_{i}\}_{i=1}^{N}$ and $[{\mathbf{y}}_{1},\cdots,{\mathbf{y}}_{N}]\in\arg\min_{Y\in{\mathbb{R}}^{n\times N},{\rm rank}(Y)\leq k}{\|X-Y\|_{p}}$.

Proof.

Let us prove first that $\inf_{\mu\in{\mathcal{P}}_{k}}W_{p}(\mu_{\rm{data}},\mu)\leq\frac{1}{N}\|X-Y^{\ast}\|_{p}$ where

[TABLE]

Let $\pi$ be a uniform distribution over $\{({\mathbf{x}}_{i},{\mathbf{y}}_{i})\}_{i=1}^{N}$ and $\mu^{\ast}$ be a uniform distribution over $\{{\mathbf{y}}_{i}\}_{i=1}^{N}$. Since $\pi\in\Pi(\mu_{{\rm data}},\mu^{\ast})$, we obtain $W_{p}(\mu_{{\rm data}},\mu^{\ast})\leq(\frac{1}{N}\sum_{i=1}^{N}\|{\mathbf{x}}_{i}-{\mathbf{y}}_{i}\|^{p})^{1/p}=\frac{1}{N^{p}}\|X-Y^{\ast}\|_{p}$. The support of $\mu^{\ast}$ is $k$-dimensional, because ${\rm rank}(Y^{\ast})\leq k$. Thus, we have $\mu^{\ast}\in{\mathcal{P}}_{k}$ and $\inf_{\mu\in{\mathcal{P}}_{k}}W_{p}(\mu_{\rm{data}},\mu)\leq W_{p}(\mu_{{\rm data}},\mu^{\ast})\leq\frac{1}{N^{p}}\|X-Y^{\ast}\|_{p}$. Now, if we prove the inverse inequality, i.e. $\inf_{\mu\in{\mathcal{P}}_{k}}W_{p}(\mu_{\rm{data}},\mu)\geq\frac{1}{N^{p}}\|X-Y^{\ast}\|_{p}$, this will imply that $\inf_{\mu\in{\mathcal{P}}_{k}}W_{p}(\mu_{\rm{data}},\mu)=\frac{1}{N^{p}}\|X-Y^{\ast}\|_{p}$ and therefore, $\inf_{\mu\in{\mathcal{P}}_{k}}W_{p}(\mu_{\rm{data}},\mu)=W_{p}(\mu_{{\rm data}},\mu^{\ast})$. This will in the end give us $\mu^{\ast}\in\arg\inf_{\mu\in{\mathcal{P}}_{k}}W_{p}(\mu_{\rm{data}},\mu)$.

Let $\{\mu_{t}\}_{1}^{\infty}$ be such that $\mu_{t}\in{\mathcal{P}}_{k}$ and $W_{p}(\mu_{\rm{data}},\mu_{t})-\inf_{\mu\in{\mathcal{P}}_{k}}W_{p}(\mu_{\rm{data}},\mu)\rightarrow 0$. Let $L_{t}$ denote a $k$-dimensional support of $\mu_{t}$ and $P_{t}$ is a projection operator onto $L_{t}$.

Let $\mu^{\ast}_{t}$ be a uniform distribution over $\{P_{t}{\mathbf{x}}_{1},\cdots,P_{t}{\mathbf{x}}_{N}\}$, i.e. $\mu^{\ast}_{t}(A)=\frac{1}{N}\sum_{i=1}^{N}[P_{t}{\mathbf{x}}_{i}\in A]$. It is easy to see that $W_{p}(\mu^{\ast}_{t},\mu_{{\rm data}})\leq W_{p}(\mu_{t},\mu_{{\rm data}})$, because $\mu^{\ast}_{t}$ and $\mu_{t}$ share the same $k$-dimensional support $L_{t}$, but the “transportation of a mass” concentrated in point ${\mathbf{x}}_{i}$ of the empirical distribution $\mu_{{\rm emp}}$ can be most optimally done by just moving it to $P_{t}{\mathbf{x}}_{i}$ (i.e. to the closest point on $L_{t}$). Thus, we have $\inf_{\mu\in{\mathcal{P}}_{k}}W_{p}(\mu_{\rm{data}},\mu)\leq W_{p}(\mu_{\rm{data}},\mu^{\ast}_{t})\leq W_{p}(\mu_{\rm{data}},\mu_{t})$, and therefore, $W_{p}(\mu_{\rm{data}},\mu^{\ast}_{t})-\inf_{\mu\in{\mathcal{P}}_{k}}W_{p}(\mu_{\rm{data}},\mu)\rightarrow 0$.

Since a set of projection operators is compact, one can always extract a subsequence $\{P_{t_{s}}\}_{s=1}^{\infty}$, such that $P_{t_{s}}\rightarrow P$. It is easy to see that $\mu^{\ast}_{t_{s}}\rightarrow\mu^{\ast\ast}$ (i.e. $W_{p}(\mu^{\ast}_{t_{s}},\mu^{\ast\ast})\rightarrow 0$) where $\mu^{\ast\ast}$ is a uniform distribution over $\{P{\mathbf{x}}_{1},\cdots,P{\mathbf{x}}_{N}\}$. For that distribution we have

[TABLE]

Thus, the infinum is attained on $\mu^{\ast\ast}$.

It is easy to see that $W_{p}(\mu_{{\rm data}},{\mu^{\ast\ast}})=\frac{1}{N^{p}}\|X-PX\|_{p}$. Since ${\rm rank}(PX)\leq k$ we obtain $W_{p}(\mu_{\rm{data}},\mu^{\ast\ast})\geq\frac{1}{N^{p}}\min_{Y\in{\mathbb{R}}^{n\times N},{\rm rank}(Y)\leq k}{\|X-Y\|_{p}}=\|X-Y^{\ast}\|_{p}$. This completes the proof. ∎

Note that in the case of $l_{1}$ norm and $p=1$, i.e. $\|{\mathbf{x}}\|=\sum_{i}|x_{i}|$, the task 66 corresponds to the well-studied robust PCA problem [35]. If, instead of the $l_{1}$-norm, we use the $l_{2}$-norm and $p\geq 1$, this leads to another task:

[TABLE]

where $\|[{\mathbf{s}}_{1},\cdots,{\mathbf{s}}_{N}]\|_{p,2}=(\sum_{i=1}^{N}\|{\mathbf{s}}_{i}\|_{2}^{p})^{1/p}$. This task has many applications in mathematics and is known as the subspace approximation problem [39] .

Appendix C Proper kernels and proof of Theorem 4

C.1 Proof of Theorem 4

Let us first show that $\langle f|M|g\rangle$ is defined for any $f=(T_{a}\otimes\delta^{n-k})_{U}\in\mathcal{G}_{k}$ and $g=(T_{b}\otimes\delta^{n-k})_{V}\in\mathcal{G}_{k}$ where $a,b\in L_{1}({\mathbb{R}}^{k})$. We have

[TABLE]

Let us denote $a_{\epsilon}=a\ast G_{\epsilon}^{k}$ and $b_{\epsilon}=b\ast G_{\epsilon}^{k}$. It is easy to see that

[TABLE]

From a well-known property of the Weierstrass transform we have

[TABLE]

From this we obtain that

[TABLE]

Thus, $\langle f_{\epsilon}|M|g_{\epsilon}\rangle$ is properly defined and

[TABLE]

where

[TABLE]

Let $U_{k},V_{k}\in{\mathbb{R}}^{n\times n}$ be matrices that comprise the first $k$ rows of $U,V$ correspondingly and $n-k$ zero rows below. Also, let $L$ denote the Lipschitz constant for $M$ such that $|M({\mathbf{x}},{\mathbf{y}})-M({\mathbf{x}}^{\prime},{\mathbf{y}}^{\prime})|\leq L(|{\mathbf{x}}-{\mathbf{x}}^{\prime}|+|{\mathbf{y}}-{\mathbf{y}}^{\prime}|)$. For the function $M_{\epsilon}({\mathbf{x}}_{1:k},{\mathbf{y}}_{1:k})$ we have:

[TABLE]

Thus, there exists bounded $\tilde{M}({\mathbf{x}}_{1:k},{\mathbf{y}}_{1:k})=M(U_{k}^{T}{\mathbf{x}},V_{k}^{T}{\mathbf{y}})$ such that

[TABLE]

Further we assume that $\epsilon>0$ is small enough, so that $M_{\epsilon}({\mathbf{x}}_{1:k},{\mathbf{y}}_{1:k})\leq C=2\max|M({\mathbf{x}},{\mathbf{y}})|$. Now we have:

[TABLE]

It is well-known (e.g. see Theorem 2.25 from [49]) that $\|a_{\epsilon}-a\|_{L_{p}}$, $\|b_{\epsilon}-b\|_{L_{p}}\rightarrow 0$, $\|a_{\epsilon}\|_{L_{1}}\leq\|a\|_{L_{1}}$ and $\|M_{\epsilon}-\tilde{M}\|_{L_{\infty}}\rightarrow 0$. Thus, $\lim_{\epsilon\rightarrow 0}\langle f_{\epsilon}|M|g_{\epsilon}\rangle$ exists and $\langle f|M|g\rangle$ is defined.

Let us now prove that $\operatorname{rank}M_{f}\leq k$. The function $f\in\mathcal{G}^{\prime}_{k}$ is such that $f=(T_{g}\otimes\delta^{n-k})_{U}$ where $\{x_{i}g\}_{i=1}^{k}\subseteq L_{1}({\mathbb{R}}^{k})$ and $U=\begin{bmatrix}{\mathbf{w}}_{1},\cdots,{\mathbf{w}}_{n}\end{bmatrix}$ is an orthogonal matrix. By construction,

[TABLE]

Let us now denote $V=\begin{bmatrix}{\mathbf{u}}_{1},\cdots,{\mathbf{u}}_{n}\end{bmatrix}\in{\mathbb{R}}^{k\times n}$ a submatrix of $U$ in which only first $k$ rows of $U$ are present. Then, the latter integral is equal to

[TABLE]

where

[TABLE]

is the Gram matrix of the collection $\{x_{i}g({\mathbf{x}}_{1:k})\}_{i=1}^{k}\subseteq L_{1}({\mathbb{R}}^{k})$.

Obviously, $\operatorname{rank}M_{f}=\operatorname{rank}\begin{bmatrix}{\rm Re\,}{\mathbf{u}}^{T}_{i}B{\mathbf{u}}_{j}\end{bmatrix}_{1\leq i,j\leq n}=\operatorname{rank}V^{T}({\rm Re\,}B)V\leq\operatorname{rank}V=k$.

Appendix D Proofs of Theorem 5 and 6

For any $f=(T_{l}\otimes\delta^{n-k})_{U}\in\mathcal{G}_{k}$ and $\sigma>0$, let us define $f_{\sigma}$ as

[TABLE]

We have $T_{f_{\sigma}}\rightarrow^{\ast}(T_{l}\otimes\delta^{n-k})_{U}$ as $\sigma\rightarrow+0$.

Lemma 4**.**

For any $f\in\mathcal{G}_{k}$, $\lim_{\sigma\rightarrow+0}\langle x_{i}f_{\sigma}|M|x_{j}f_{\sigma}\rangle=0$, for any $(i,j)\notin\{1,...,k\}^{2}$, and $\sup_{\sigma\in[0,1]}\langle x_{i}f_{\sigma}|M|x_{j}f_{\sigma}\rangle<\infty$, for any $(i,j)\in\{1,...,k\}^{2}$.

Proof.

W.l.o.g. we can assume that $f=T_{l}\otimes\delta^{n-k},l\in{\mathcal{S}}({\mathbb{R}}^{k})$. If $i>k,j\leq k$ we have

[TABLE]

where $P({\mathbf{x}})=\int_{{\mathbb{R}}^{n}}\frac{1}{\sqrt{2\pi\sigma^{2}}^{n-k}}y_{j}M({\mathbf{x}},{\mathbf{y}})e^{-\frac{|{\mathbf{y}}_{k+1:n}|^{2}}{2\sigma^{2}}}l_{\sigma}({\mathbf{y}}_{1:k})d{\mathbf{y}}$. Using the Hölder inequality we obtain

[TABLE]

Since $|M({\mathbf{x}},{\mathbf{y}})|\leq\gamma$ for some $\gamma$, we have

[TABLE]

Thus,

[TABLE]

Using $\|l_{\sigma}\|_{L_{1}({\mathbb{R}}^{k})}-\|l\|_{L_{1}({\mathbb{R}}^{k})}\mathop{\rightarrow}\limits^{\sigma\rightarrow+0}0$, $\|y_{j}l_{\sigma}\|_{L_{1}({\mathbb{R}}^{k})}-\|y_{j}l\|_{L_{1}({\mathbb{R}}^{k})}\mathop{\rightarrow}\limits^{\sigma\rightarrow+0}0$, we see the boundedness of $\|l_{\sigma}\|_{L_{1}({\mathbb{R}}^{k})}\gamma\|y_{j}l_{\sigma}\|_{L_{1}({\mathbb{R}}^{k})}$ and proceed

[TABLE]

It is easy to see that $\|\frac{1}{\sqrt{2\pi\sigma^{2}}^{n-k}}x_{i}e^{-\frac{|{\mathbf{x}}_{k+1:n}|^{2}}{2\sigma^{2}}}\|_{L_{1}({\mathbb{R}}^{n-k})}\rightarrow 0$ as $\sigma\rightarrow 0$, therefore $\langle x_{i}f_{\sigma}|M|x_{j}f_{\sigma}\rangle\rightarrow 0$.

Similarly, we can prove that $\langle x_{i}f_{\sigma}|M|x_{j}f_{\sigma}\rangle\rightarrow 0$ if $i,j>k$.

The entries of the main $k\times k$ minor $[\langle x_{i}f_{\sigma}|M|x_{j}f_{\sigma}\rangle]_{1\leq i,j\leq k}$ are bounded, due to

[TABLE]

Again, using $\|l_{\sigma}\|_{L_{1}({\mathbb{R}}^{k})}-\|l\|_{L_{1}({\mathbb{R}}^{k})}\mathop{\rightarrow}\limits^{\sigma\rightarrow+0}0$, $\|y_{j}l_{\sigma}\|_{L_{1}({\mathbb{R}}^{k})}-\|y_{j}l\|_{L_{1}({\mathbb{R}}^{k})}\mathop{\rightarrow}\limits^{\sigma\rightarrow+0}0$, we obtain the boundedness of RHS. ∎

Corollary 2**.**

For any $f\in\mathcal{G}_{k}$, $\lim_{\sigma\rightarrow 0}R(f_{\sigma})=0$.

Proof.

W.l.o.g. we can assume that $f=T_{l}\otimes\delta^{n-k},l\in{\mathcal{S}}({\mathbb{R}}^{k})$. By lemma, all entries of $M_{f_{\sigma}}$ except those of the main $k\times k$ minor approach [math] as $\sigma\rightarrow 0$. This means that $\lim_{\sigma\rightarrow+0}Q(f_{\sigma})=0$, where $Q(f_{\sigma})=\sum_{i=k+1}^{n}\langle x_{i}f_{\sigma}|M|x_{i}f_{\sigma}\rangle$. Let ${\mathbf{v}}_{1},\cdots,{\mathbf{v}}_{n}$ be unit eigenvectors of $M_{f_{\sigma}}$ corresponding to the eigenvalues $\lambda_{1}\geq\cdots\geq\lambda_{n}$, $P=\sum_{i=k+1}^{n}{\mathbf{e}}_{i}{\mathbf{e}}^{T}_{i}$, then

[TABLE]

Since $R(f_{\sigma})\leq Q(f_{\sigma})$, we obtain $\lim_{\sigma\rightarrow 0}R(f_{\sigma})=0$. ∎

D.0.1 Proof of Theorem 5

Proof.

Suppose that a sequence $\{f_{i}\}^{\infty}_{s=1}\subseteq{\mathcal{S}}({\mathbb{R}}^{n})$ regularly solves (7) and $T\in\mathop{\rm Lim}\limits_{i\rightarrow\infty}f_{i}$. W.l.o.g. we can assume that $T_{f_{i}}\rightarrow^{\ast}T$ and $\operatorname{Tr}(M_{f_{i}})$ is bounded and $I(f_{i})+\lambda_{i}R(f_{i})\leq\inf\limits_{f\in{\mathcal{S}}({\mathbb{R}}^{n})}I(f)+\lambda_{i}R(f)+\epsilon_{i},\epsilon_{i}\rightarrow 0$. Below we use continuity of $I$ and corollary 2:

[TABLE]

from which we conclude that $\lambda_{i}R(f_{i})\leq\inf\limits_{f\in\mathcal{G}_{k}}I(f)+\epsilon_{i}$ and, therefore, $R(f_{i})\mathop{\rightarrow}\limits^{i\rightarrow\infty}0$.

For each $l$, let us define $P_{l}$ as the projection operator to a subspace spanned by first principal components of the matrix $\sqrt{M_{f_{l}}}$, i.e.

[TABLE]

where ${\mathbf{v}}^{l}_{1},...,{\mathbf{v}}^{l}_{k}$ are orthonormal eigenvectors that correspond to $k$ largest eigenvalues of $\sqrt{M_{f_{l}}}$. From the Eckart-Young-Mirsky theorem we see that $R(f_{l})=\|\sqrt{M_{f_{l}}}-P_{l}\sqrt{M_{f_{l}}}\|^{2}_{F}$. Since a set of all projection operators $\{P\in{\mathbb{R}}^{n\times n}|P^{2}=P,P^{T}=P\}$ is a compact subset of ${\mathbb{R}}^{n^{2}}$, one can always find a projection operator $P=\sum_{i=1}^{k}{\mathbf{v}}_{i}{\mathbf{v}}^{T}_{i}$ and a growing subsequence $\{l_{s}\}$ such that $\|P_{l_{s}}-P\|_{F}\rightarrow 0$ as $s\rightarrow\infty$. Thus, for the subsequence $\{f_{l_{s}}\}$ we have

[TABLE]

and using the boundedness of $\operatorname{Tr}(M_{f_{s}})$ we obtain $\|\sqrt{M_{f_{l_{s}}}}-P\sqrt{M_{f_{l_{s}}}}\|_{F}\rightarrow 0$.

Since $\|\sqrt{M_{f_{l_{s}}}}-P\sqrt{M_{f_{l_{s}}}}\|_{F}\rightarrow 0$, let us complete ${\mathbf{v}}_{1},...,{\mathbf{v}}_{k}$ to an orthonormal basis ${\mathbf{v}}_{1},...,{\mathbf{v}}_{n}$ and make the change of variables $y_{i}={\mathbf{v}}^{T}_{i}{\mathbf{x}}$. Let us denote $V=\begin{bmatrix}{\mathbf{v}}_{1},...,{\mathbf{v}}_{n}\end{bmatrix}$ and let $V^{T}=\begin{bmatrix}{\mathbf{w}}_{1},...,{\mathbf{w}}_{n}\end{bmatrix}$. Then, after that change of variables any function $f({\mathbf{x}})$ corresponds to $f^{\prime}({\mathbf{y}})=f(V{\mathbf{y}})$ and the kernel $M$ corresponds to $M^{\prime}({\mathbf{y}},{\mathbf{y}}^{\prime})=M(V{\mathbf{y}},V{\mathbf{y}}^{\prime})$. After we apply that change of variables in the integral expression of $\langle x_{i}f|M|x_{j}f\rangle$, we obtain

[TABLE]

I.e. $M_{f}=VM^{\prime}_{f^{\prime}}V^{T}$, or $M^{\prime}_{f^{\prime}}=V^{T}M_{f}V$. Note that $P=VI_{n}^{k}V^{T}$ where $I^{k}_{n}$ is a diagonal matrix whose main $k\times k$ minor is the identity matrix, and all other entries are zeros. Using the fact that the Frobenius norm of orthogonally similar matrices are equal and the identity $V^{T}\sqrt{M_{f_{l_{s}}}}V=\sqrt{V^{T}M_{f_{l_{s}}}V}$, we obtain

[TABLE]

Thus, the property $\|\sqrt{M_{f_{l_{s}}}}-P\sqrt{M_{f_{l_{s}}}}\|_{F}\rightarrow 0$ implies

[TABLE]

Moreover, for $i=j$ we have ${\rm Re\,}\langle y_{i}f^{\prime}_{l_{s}}|M^{\prime}|y_{j}f^{\prime}_{l_{s}}\rangle=\langle y_{i}f^{\prime}_{l_{s}}|M^{\prime}|y_{j}f^{\prime}_{l_{s}}\rangle$. It is easy to see that after the change of variables we still have $f^{\prime}_{l_{s}}\rightarrow^{\ast}T_{V}$. Since $f^{\prime}_{l_{s}}\in{\mathcal{S}}({\mathbb{R}}^{n})$, we have $y_{i}f^{\prime}_{l_{s}}\in{\mathcal{S}}({\mathbb{R}}^{n})$ and, therefore, $y_{i}f^{\prime}_{l_{s}}\in L_{2}({\mathbb{R}}^{n})$. Let us treat now $M^{\prime}$ as an operator ${\rm O}_{M^{\prime}}:L_{2}({\mathbb{R}}^{n})\rightarrow L_{2}({\mathbb{R}}^{n}),{\rm O}_{M^{\prime}}[f]({\mathbf{x}})=\int_{{\mathbb{R}}^{n}}M^{\prime}({\mathbf{x}},{\mathbf{y}})f({\mathbf{y}})d{\mathbf{y}}$. Let us take any function $\phi\in L_{2}({\mathbb{R}}^{n})$ such that $\psi={\rm O}_{M^{\prime}}[\phi]\in{\mathcal{S}}({\mathbb{R}}^{n})$. Since ${\rm O}_{M^{\prime}}$ is a strictly positive self-adjoint operator, by the Cauchy-Schwarz inequality, we obtain

[TABLE]

Therefore, for any $\psi\in{\rm Range\,}[{\rm O}_{M^{\prime}}]\cap{\mathcal{S}}({\mathbb{R}}^{n})$ and $i>k$ we have $\lim_{s\rightarrow\infty}\langle y_{i}f^{\prime}_{l_{s}},\psi\rangle=\lim_{s\rightarrow\infty}\langle f^{\prime}_{l_{s}},y_{i}\psi\rangle=0$. Since $f^{\prime}_{l_{s}}\rightarrow^{\ast}T_{V}$ we obtain $\langle T_{V},y_{i}\psi\rangle=\langle y_{i}T_{V},\psi\rangle=0$ for any $\psi\in{\rm Range\,}[{\rm O}_{M^{\prime}}]\cap{\mathcal{S}}({\mathbb{R}}^{n})$. But the denseness of ${\rm Range\,}[{\rm O}_{M^{\prime}}]\cap{\mathcal{S}}({\mathbb{R}}^{n})$ in ${\mathcal{S}}({\mathbb{R}}^{n})$ implies that $y_{i}T_{V}=0$.

Using lemma 3 and $(T_{V})_{V^{T}}=T$ we obtain $T\in\mathcal{G}^{\prime}_{k}$. Thus, we proved that $T_{f_{i}}\rightarrow T\in\mathcal{G}^{\prime}_{k}$.

Since $I(f_{i})\leq I(f_{i})+\lambda_{i}R(f_{i})\leq\inf\limits_{f\in\mathcal{G}^{\prime}_{k}}I(f)+\epsilon_{i}$ and $I$ is continuous, we finally get that $I(T)\leq\inf\limits_{f\in\mathcal{G}^{\prime}_{k}}I(f)$, i.e. $T\in{\rm Arg}\min\limits_{f\in\mathcal{G}^{\prime}_{k}}I(f)$. ∎

D.0.2 Proof of Theorem 6

Proof.

Suppose $f^{\ast}\in{\rm Arg}\min\limits_{f\in\mathcal{G}^{\prime}_{k}}I(f)\bigcap\mathcal{B}_{k}$, i.e. $f^{\ast}\in\mathcal{B}_{k}$ and $I(f^{\ast})=\min\limits_{f\in\mathcal{G}^{\prime}_{k}}I(f)$. Since $f^{\ast}\in\mathcal{B}_{k}$, then there exists a sequence $\{s^{i}\}\subseteq\mathcal{G}_{k}$ such that $T_{s^{i}}\rightarrow^{\ast}f^{\ast}$ and $\sup\limits_{i}{\rm Tr\,}M_{s^{i}}<\infty$.

Let us define $s^{i}_{\sigma}\in{\mathcal{S}}({\mathbb{R}}^{n})$ as $T_{s^{i}_{\sigma}}=T_{s^{i}}\ast G_{\sigma}^{n}$. Since $\lim_{\sigma\rightarrow 0}R(s^{i}_{\sigma})=0$ (lemma 4), there exists $\sigma_{i}>0$, such that $R(s^{i}_{\sigma})<\frac{1}{i}$ whenever $0<\sigma\leq\sigma_{i}$. Also, by definition ${\rm Tr\,}M_{s^{i}}=\lim_{\sigma\rightarrow 0}{\rm Tr\,}M_{s^{i}_{\sigma}}$. Therefore, there exists $\sigma^{\prime}_{i}>0$, such that ${\rm Tr\,}M_{s^{i}_{\sigma}}<{\rm Tr\,}M_{s^{i}}+1$ whenever $0<\sigma\leq\sigma^{\prime}_{i}$.

If we set $\sigma_{i}^{\ast}=\min\{\sigma_{i},\sigma^{\prime}_{i},\frac{1}{i}\}$, then a sequence $\{s^{i}_{\sigma^{\ast}_{i}}\}\subseteq{\mathcal{S}}({\mathbb{R}}^{n})$ satisfies

[TABLE]

and (using lemma 2) $T_{s^{i}_{\sigma^{\ast}_{i}}}\rightarrow^{\ast}f^{\ast}$.

Due to the continuity of $I$ we have $\lim_{i\rightarrow\infty}I(s^{i}_{\sigma^{\ast}_{i}})=I(f^{\ast})$. Now we set $f_{i}=s^{i}_{\sigma^{\ast}_{i}}$, $\lambda_{i}=\frac{1}{\sqrt{R(f_{i})}}$ and we obtain the needed sequence:

[TABLE]

where ${\rm Tr\,}M_{f_{i}}$ is bounded. It remains to check that our sequence regularly solves (7), i.e. $\lim_{i\rightarrow\infty}\inf\limits_{f\in{\mathcal{S}}({\mathbb{R}}^{n})}I(f)+\lambda_{i}R(f)=I(f^{\ast})$ (this will imply $\lim_{i\rightarrow\infty}I(f_{i})+\lambda_{i}R(f_{i})-\inf\limits_{f\in{\mathcal{S}}({\mathbb{R}}^{n})}I(f)+\lambda_{i}R(f)=0$). The inequality in one direction is obvious,

[TABLE]

Let us prove the inverse inequality.

Since ${\rm rsol\,}(I(f),R(f))\neq\emptyset$, there exists $\{\tilde{f}_{i}\}\subseteq{\mathcal{S}}({\mathbb{R}}^{n})$ such that

[TABLE]

and $a=\lim_{i\rightarrow+\infty}T_{\tilde{f}_{i}}$. From theorem 5 we obtain $a\in{\rm Arg}\min\limits_{f\in\mathcal{G}^{\prime}_{k}}I(f)$.

One can always find a subset $\{\tilde{\lambda}_{d_{i}}\}\subseteq\{\tilde{\lambda}_{i}\}$ such that $\tilde{\lambda}_{d_{i}}<\lambda_{i}$, $\tilde{\lambda}_{d_{i}}\rightarrow\infty$ and obtain

[TABLE]

Therefore,

[TABLE]

This proves that $\{{f_{i}}\}$ regularly solves (7) and $\lim_{i\rightarrow\infty}f_{i}=f^{\ast}$ i.e. $f^{\ast}\in{\rm rsol\,}(I(f),R(f))$. ∎

Appendix E The alternating scheme in the dual space for $M({\mathbf{x}},{\mathbf{y}})=\zeta({\mathbf{x}}-{\mathbf{y}})$

When $M({\mathbf{x}},{\mathbf{y}})=\zeta({\mathbf{x}}-{\mathbf{y}})$, the alternating scheme 1 allows for a reformulation in the dual space. By this we mean that in Scheme 1 we substitute $\hat{\phi}_{t}$ for the original $\phi_{t}$. If the primal Scheme 1 deals with operators $S_{\phi},S_{\phi_{t-1}}$, the dual version deals with vectors of functions $\sqrt{{\hat{\zeta}}}\frac{\partial{\hat{\phi}}}{\partial{\mathbf{x}}},\sqrt{{\hat{\zeta}}}\frac{\partial{\hat{\phi}}_{t-1}}{\partial{\mathbf{x}}}$. The substitution is based on the following simple fact:

Theorem 13**.**

If $M({\mathbf{x}},{\mathbf{y}})=\zeta({\mathbf{x}}-{\mathbf{y}}),\zeta,\hat{\zeta}\in C({\mathbb{R}}^{n})$ and $\forall{\mathbf{x}}\,\,\hat{\zeta}({\mathbf{x}})>0$, then there exist constants $c_{1}$ and $c_{2}$ such that $\|S_{\phi}-P_{t-1}S_{\phi_{t-1}}\|_{\ast}^{2}=c_{1}\|\,\,\|\frac{\partial{\hat{\phi}}}{\partial{\mathbf{x}}}-P_{t-1}\frac{\partial{\hat{\phi}}_{t-1}}{\partial{\mathbf{x}}}\|_{2}\,\,\|^{2}_{L_{2,\hat{\zeta}}({\mathbb{R}}^{n})}$ and $\langle x_{i}f|M|x_{j}f\rangle=c_{2}\langle\frac{\partial{\hat{f}}}{\partial x_{i}},\frac{\partial{\hat{f}}}{\partial x_{j}}\rangle_{L_{2,\hat{\zeta}}({\mathbb{R}}^{n})}$

Proof.

Let $f:{\mathbb{R}}^{n}\rightarrow{\mathbb{C}}$ be such that $\|x_{i}f\|_{L_{2}({\mathbb{R}}^{n})}<\infty$.

[TABLE]

Since $S_{f}[\psi]_{i}={\rm Re\,}\langle(S_{f})_{i},\psi\rangle\propto{\rm Re\,}\langle\widehat{(S_{f})_{i}},\hat{\psi}\rangle$, we obtain

[TABLE]

where $\kappa$ is a constant.

Let us now introduce a vector of functions $V_{f}=\begin{bmatrix}(S_{f})_{1},\cdots,(S_{f})_{n}\end{bmatrix}^{T}\in L^{n}_{2}({\mathbb{R}}^{n})$. Using 105 we obtain $\widehat{(S_{f})_{i}}=\kappa\sqrt{\hat{\zeta}}\frac{\partial{\hat{f}}}{\partial x_{i}}$, and therefore ${\widehat{V}_{f}}=\kappa\sqrt{\hat{\zeta}}\frac{\partial{\hat{f}}}{\partial{\mathbf{x}}}$. Thus, the expression $\|S_{\phi}-P_{t-1}S_{\phi_{t-1}}\|_{\ast}^{2}$ in the alternating scheme can be rewritten as

[TABLE]

The matrix $M_{f}$ can also be calculated from $\hat{f}$ using the following identity:

[TABLE]

∎

Let us introduce a function $\tilde{I}$ such that $\tilde{I}(f)=I(\hat{f})$. Then, we see that all steps in Scheme 1 can be performed with ${\hat{\phi}_{t}}$ rather than with $\phi_{t}$, using the algorithm 3.

Informally, the dual algorithm works as follows: at each iteration $t$ we compute a function ${\hat{\phi}}_{t}$ adapting it to data (the term $\tilde{I}({\hat{\phi}})$) and adapting its gradient field to the rank reduced gradient field of the previous ${\hat{\phi}}_{t-1}$. For a sufficiently large $T$, it will converge and ${\hat{\phi}}_{T}\approx{\hat{\phi}}_{T-1}$. Then, the second term in the last step will be approximately equal to $\lambda\|\,\,\|\frac{\partial{\hat{\phi}}_{T}}{\partial{\mathbf{x}}}-P_{T-1}\frac{\partial{\hat{\phi}}_{T}}{\partial{\mathbf{x}}}\|_{2}\,\,\|^{2}_{L_{2,\hat{\zeta}}({\mathbb{R}}^{n})}$, enforcing $\frac{\partial{\hat{\phi}}_{T}}{\partial{\mathbf{x}}}\approx P_{T-1}\frac{\partial{\hat{\phi}}_{T}}{\partial{\mathbf{x}}}$ for random ${\mathbf{x}}\sim\frac{\hat{\zeta}}{\|\hat{\zeta}\|_{L_{1}}}$. Thus, gradients $\frac{\partial{\hat{\phi}}_{T}}{\partial{\mathbf{x}}}$ lie in a $k$-dimensional subspace ${\rm col\,}P_{T-1}$. This last property is a characteristic property of functions from ${\mathcal{F}}_{k}$.

Absolutely analogously to the Algorithm 3, one can construct a dual algorithm that deals with inverse Fourier transforms of functions, i.e. with $\mathcal{F}^{-1}[\phi]$, $\mathcal{F}^{-1}[\phi_{t}]$, $M_{t}=\begin{bmatrix}{\rm Re\,}\langle\frac{\partial{\mathcal{F}^{-1}[\phi_{t}]}}{\partial x_{i}},\frac{\partial{\mathcal{F}^{-1}[\phi_{t}]}}{\partial x_{j}}\rangle_{L_{2,\mathcal{F}^{-1}[\zeta]}({\mathbb{R}}^{n})}\end{bmatrix}$ etc. This version of the dual alternating scheme will be used for designing numerical algorithms for the Gaussian MMD-PCA and HM-MMD-PCA.

Appendix F Proofs for Section 7

Proof of Theorem 10..

Let $X=[{\mathbf{x}}_{1},\cdots,{\mathbf{x}}_{N}]$. Note that $H_{f}=XSX^{T}$ where $S=[P({\mathbf{x}}_{i}\cdot{\mathbf{x}}_{j})]_{i,j\in[N]}$. For any $U\in\mathcal{O}(n)$ and we have

[TABLE]

Let $U=[{\mathbf{u}}_{1},\cdots,{\mathbf{u}}_{n}]$ where $\{{\mathbf{u}}_{i}\}_{i=1}^{n}$ are eigenvectors such that $H_{f}{\mathbf{u}}_{i}=\lambda_{i}{\mathbf{u}}_{i}$. Then, the rotated distribution $f_{U}$ is such that $H_{f_{U}}$ is diagonal. Note that $\inf_{\nu\in\mathcal{P}_{k}}\|f_{U}-T_{\nu}\|_{K}=\inf_{\nu\in\mathcal{P}_{k}}\|f-T_{\nu}\|_{K}$.

Therefore, w.l.o.g. we can assume that principal components of $H_{f}$ are ${\mathbf{e}}_{1},\cdots,{\mathbf{e}}_{n}$ and $H_{f}{\mathbf{e}}_{i}=\lambda_{i}{\mathbf{e}}_{i},i\in[n]$, where $\{{\mathbf{e}}_{i}\}_{i=1}^{n}$ is a canonical basis in ${\mathbb{R}}^{n}$ and $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{n}$ are eigenvalues of $H_{f}$. From the latter we conclude that $\langle x_{i}f|H|x_{j}f\rangle=\lambda_{i}\delta_{ij}$. Using $P({\mathbf{x}}\cdot{\mathbf{y}})=\sum_{j=0}^{l-1}c_{j}\sum_{\alpha\in({\mathbb{N}}\cup\{0\})^{n}:|\alpha|=j}\frac{j!}{\alpha!}{\mathbf{x}}^{\alpha}{\mathbf{y}}^{\alpha}$, we obtain

[TABLE]

For an input distribution $f({\mathbf{x}})=\frac{1}{N}\sum_{i=1}^{N}\delta^{n}({\mathbf{x}}-{\mathbf{x}}_{i})$, let us denote $f^{\prime}({\mathbf{x}})=\frac{1}{N}\sum_{i=1}^{N}\delta^{k}({\mathbf{x}}_{1:k}-({\mathbf{x}}_{i})_{1:k})\otimes\delta^{n-k}({\mathbf{x}}_{k+1:n})$, where ${\mathbf{x}}=[{\mathbf{x}}_{1:k},{\mathbf{x}}_{k+1:n}]$ and ${\mathbf{z}}_{1:k}\in{\mathbb{R}}^{k}$ equals the first $k$ components of ${\mathbf{z}}\in{\mathbb{R}}^{n}$. By construction,

[TABLE]

for $i_{1},\cdots,i_{s}\in[k]$ and

[TABLE]

whenever $i_{j}\in[n]\setminus[k]$ for at least one $j\in[s]$. Therefore,

[TABLE]

The second sum $F_{2}$ equals $\sum_{i=k+1}^{n}\lambda_{i}$. Let us compare the first sum, $F_{1}$, with the second, $F_{2}$. $F_{1}$ is a sum of positive factors of $({\mathbb{E}}_{{\mathbf{x}}\sim f}{\mathbf{x}}^{\alpha})^{2}$ where $\alpha_{1:k}\neq 0,\alpha_{k+1:n}\neq 0$. Let $e_{i}\in({\mathbb{N}}\cup\{0\})^{n}$ be an $i$th canonical unit vector and $\alpha_{i}$ denote an $i$th component of $\alpha$. The coefficient in front of $({\mathbb{E}}_{{\mathbf{x}}\sim f}{\mathbf{x}}^{\alpha})^{2}$ in $F_{1}$ equals $A_{\alpha}=c_{|\alpha|-1}(|\alpha|-1)!\sum_{i\in[k]:\alpha_{i}>0}\frac{1}{(\alpha-e_{i})!}=\frac{c_{|\alpha|-1}(|\alpha|-1)!}{\alpha!}|\alpha_{1:k}|$ and the coefficient in front of $({\mathbb{E}}_{{\mathbf{x}}\sim f}{\mathbf{x}}^{\alpha})^{2}$ in $F_{2}$ equals $B_{\alpha}=c_{|\alpha|-1}(|\alpha|-1)!\sum_{i\in[n]\setminus[k]:\alpha_{i}>0}\frac{1}{(\alpha-e_{i})!}=\frac{c_{|\alpha|-1}(|\alpha|-1)!}{\alpha!}|\alpha_{k+1:n}|$. Since $|\alpha_{k+1:n}|\geq 1$ and $|\alpha_{1:k}|\leq l-1$, we conclude that $A_{\alpha}\leq(l-1)B_{\alpha}$. Therefore, $F_{1}\leq(l-1)F_{2}$.

Thus, overall we have

[TABLE]

From $\|T_{\mu}-T_{\nu}\|^{2}_{K}=d_{\rm MMD}(\mu,\nu)^{2}$ the statement of theorem directly follows. ∎

In our proof of Theorem 11. we will need the following classical theorem.

Theorem 14** (The Gaussian Poincaré inequality).**

Let $g:{\mathbb{R}}^{n}\to{\mathbb{R}}$ be a smooth function, then

[TABLE]

Proof of Theorem 11..

Let us assume w.l.o.g. that principal components of $H_{f}$ are ${\mathbf{e}}_{1},\cdots,{\mathbf{e}}_{n}$ and $H_{f}{\mathbf{e}}_{i}=\lambda_{i}{\mathbf{e}}_{i},i\in[n]$, where $\{{\mathbf{e}}_{i}\}_{i=1}^{n}$ is a canonical basis in ${\mathbb{R}}^{n}$ and $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{n}$ are eigenvalues of $H_{f}$. From the latter we conclude that $\langle x_{i}f|H|x_{j}f\rangle=\lambda_{i}\delta_{ij}$.

Note that $H({\mathbf{x}},{\mathbf{y}})=\mathcal{F}^{-1}[G_{\sigma}]({\mathbf{x}}-{\mathbf{y}})$. Let $\hat{f}=\mathcal{F}[f]$,

[TABLE]

and $f^{\prime}=\mathcal{F}^{-1}[\hat{f^{\prime}}]$ where ${\mathbf{x}}=[{\mathbf{x}}_{1:k},{\mathbf{x}}_{k+1:n}]$. From the isometry of the Fourier transform, we have

[TABLE]

The latter expression decomposes into two terms. The first is

[TABLE]

The second is

[TABLE]

The latter can be bounded using the Gaussian Poincaré inequality by

[TABLE]

After changing an order of summations one can bound the internal sum using integration by parts, i.e.

[TABLE]

Then, using the Cauchy–Schwarz inequality we bound the latter by

[TABLE]

The first term equals $(\frac{1}{C_{n}}\lambda_{j})^{1/2}$ and the second term, after making the inverse Fourier transform, is bounded by

[TABLE]

where $T_{i}=\mathcal{F}^{-1}[\frac{x_{i}^{2}G_{\sigma}}{\sigma^{2}}]=e^{-\sigma^{2}\|{\mathbf{x}}\|^{2}/2}(1-\sigma^{2}x_{i}^{2})$. Since $T_{i}({\mathbf{x}}-{\mathbf{y}})=H({\mathbf{x}},{\mathbf{y}})(1-\sigma^{2}x_{i}^{2}-\sigma^{2}y_{i}^{2}+2\sigma^{2}x_{i}y_{i})$, we have $\langle x_{i}x_{j}f|T_{i}({\mathbf{x}}-{\mathbf{y}})|y_{i}y_{j}f\rangle=\langle x_{i}x_{j}f|H|x_{i}x_{j}f\rangle-2\sigma^{2}\langle x^{3}_{i}x_{j}f|H|x_{i}x_{j}f\rangle+2\sigma^{2}\langle x^{2}_{i}x_{j}f|H|x^{2}_{i}x_{j}f\rangle$. After noting that

[TABLE]

we finally obtain

[TABLE]

where $M=\mathcal{O}(\int_{{\mathbb{R}}^{n}}\|{\mathbf{x}}\|^{2}|f({\mathbf{x}})|d{\mathbf{x}}+\sqrt{2}\sigma\sqrt{\int_{{\mathbb{R}}^{n}}\|{\mathbf{x}}\|^{4}|f({\mathbf{x}})|d{\mathbf{x}}\int_{{\mathbb{R}}^{n}}\|{\mathbf{x}}\|^{2}|f({\mathbf{x}})|d{\mathbf{x}}}+\sqrt{2}\sigma\int_{{\mathbb{R}}^{n}}\|{\mathbf{x}}\|^{3}|f({\mathbf{x}})|d{\mathbf{x}})$.

By construction, $f^{\prime}\in\mathcal{G}_{k}^{\prime}$ and it can be approached by elements of $\mathcal{G}_{k}$ w.r.t. norm $\|\cdot\|_{K}$. The statement of theorem directly follows from this observation. ∎

Appendix G A numerical alternating scheme for the Gaussian MMD-PCA

G.1 Structure of $\mathcal{F}^{-1}[\mathcal{P}_{k}]$

From theorems 1 and 2, $\mathcal{F}^{-1}[\mathcal{P}_{k}]\subseteq\overline{\mathcal{F}_{k}}^{\ast}$. In fact, Bochner’s theorem [50] gives us that the inverse Fourier transform of any positive finite Borel measure is a continuous positive definite function. That is, if $f\in\mathcal{F}^{-1}[\mathcal{P}]$, then for any distinct ${\mathbf{y}}_{1},\cdots,{\mathbf{y}}_{s}\in{\mathbb{R}}^{n}$ the matrix $[f({\mathbf{y}}_{i}-{\mathbf{y}}_{j})]_{i,j=\overline{1,n}}$ is positive semidefinite. Since $\mu({\mathbb{R}}^{n})=1$, we additionally have $f({\mathbf{0}})=1$. Let ${\rm PDF}$ denote the set of all continuous positive definite functions on ${\mathbb{R}}^{n}$ and

[TABLE]

Thus, the following characterization of $\mathcal{F}^{-1}[\mathcal{P}_{k}]$ becomes evident.

Theorem 15**.**

$\mathcal{F}^{-1}[\mathcal{P}_{k}]=\mathcal{M}_{k}$.

G.2 The dual form of the Gaussian MMD-PCA

Recall that $k({\mathbf{x}})=G_{h}^{n}({\mathbf{x}})$. Let us define another Gaussian kernel $\gamma({\mathbf{x}})=e^{-\frac{h^{2}|{\mathbf{x}}|^{2}}{2}}=\mathcal{F}^{-1}[k]$. Let $p_{\rm{data}}({\mathbf{x}})$ denote the characteristic function of the random vector ${\mathbf{X}}_{\rm{data}}\sim\mu_{\rm{data}}$. By definition, $p_{\rm{data}}({\mathbf{x}})={\mathbb{E}}[e^{{\rm i}{\mathbf{X}}_{\rm{data}}^{T}{\mathbf{x}}}]=\frac{1}{N}\sum_{i=1}^{N}e^{{\rm i}{\mathbf{x}}^{T}_{i}{\mathbf{x}}}$. Thus, $p_{\rm{data}}=\mathcal{F}^{-1}[\mu_{\rm{data}}]$ and $\mu_{\rm{data}}=\mathcal{F}[p_{\rm{data}}]$.

Using the isometry property of the inverse Fourier transform for $L_{2}({\mathbb{R}}^{n})$ and the convolution theorem, we see that

[TABLE]

Thus, from Theorem 15 we obtain that the task 17 is equivalent to

[TABLE]

where $L_{2,\gamma^{2}}({\mathbb{R}}^{n})\overset{def}{=}L_{2,\nu}({\mathbb{R}}^{n})$ with $d\nu=\gamma^{2}d{\mathbf{x}}$.

G.3 Algorithms for the Gaussian MMD-PCA

Let $\Pi_{k}:{\mathcal{G}_{k}}\rightarrow\{1,+\infty\}$ and ${\text{M}}_{k}:{\mathcal{F}_{k}}\rightarrow\{1,+\infty\}$ be simple penalty functions:

[TABLE]

[TABLE]

Then, the task 17 is equivalent to

[TABLE]

From the result of the previous section we see that if $I(\phi)=\tilde{I}({\hat{\phi}})$, then

[TABLE]

Thus, the Algorithm 4 is an adaptation of Algorithm 3 to MMD-PCA.

If the function $p_{\rm{data}}$ is real-valued, then only real-valued functions can appear in the Algorithm 4. This assumption can be satisfied by adding reflections of initial points to the dataset (after it was centered).

At step 1, we search over $q$ given in the following parameterized form:

[TABLE]

where $\alpha_{i}>0$ and $\sum_{i=1}^{{\rm nn}}\alpha_{i}=1$. In our implementation, we set $[\alpha_{i}]_{i=\overline{1,{\rm nn}}}={\rm softmax}([u_{i}]_{i=\overline{1,{\rm nn}}})$ and $u_{i}$’s are unconstrained. The number of neurons in a single layer neural network with a cosine activation function, ${\rm nn}$, is a hyperparameter. Let us denote parameters $\{\omega_{i},u_{i}\}_{i=1}^{\rm nn}$ by $\theta$. It is easy to see the function $q_{\theta}$ is positive definite. Moreover, using Theorem 2 from [51], it can be shown that a set of all such functions, i.e. the convex hull of $\{cos(\omega^{T}{\mathbf{x}})|\omega\in{\mathbb{R}}^{n}\}$, is dense in a set of real-valued functions from $\mathcal{M}_{k}$. Though this parameterization is quite natural, finding architectures with more expressive power in a space of real-valued positive definite functions is an open problem.

Now, to minimize

[TABLE]

with stochastic gradient descent methods (in our case, the Adam optimizer) we need to have an unbiased estimator of

[TABLE]

where ${\mathbf{z}}\sim f$ denotes that the random vector ${\mathbf{z}}$ is sampled according to the probability density function $\frac{f({\mathbf{x}})}{\int_{{\mathbb{R}}^{n}}f({\mathbf{x}})d{\mathbf{x}}}$. Thus, a natural estimator of the gradient is

[TABLE]

where $\{{\mathbf{z}}_{i}\}_{i=1}^{m}\sim^{iid}\gamma^{2}$ and \{\text{\boldmath\xi}_{i}\}_{i=1}^{m}\sim^{iid}{\hat{\zeta}}.

The last important issue with the practical numerical algorithm is the calculation of $M_{t}$ at step 2. By construction,

[TABLE]

In practice we sample ${\boldsymbol{\chi}}_{1},\cdots,{\boldsymbol{\chi}}_{l}\sim{\hat{\zeta}}$ and estimate $M_{t}$ as

[TABLE]

The details of the numerical algorithm 5 are given below. In all our experiments with MMD-PCA we set ${\hat{\zeta}}=\gamma^{2}$.

Appendix H A numerical alternating scheme for HM-MMD-PCA

H.1 The dual form of HM-MMD-PCA

Due to a well-known relationship between moments of the probability measure $\mu$ and its characteristic function $p$, i.e. ${\rm i}^{s}m_{i_{1}\cdots i_{s}}=\frac{\partial^{s}p({\mathbf{0}})}{\partial x_{i_{1}}\cdots\partial x_{i_{s}}}$, the task (20) is equivalent to

[TABLE]

Note that the maximum mean discrepancy distance for the Gaussian kernel and the distance based on higher moments are substantially different. Indeed, even if we set $h$ as a large value (which makes $\frac{1}{h}\approx 0$), the MMD distance, unlike the HM distance, neglects higher order derivatives of the characteristic functions in the neigbourhood of the origin. Moreover, from the dual form (135) it is clear that $d_{\rm HM}(\mu_{\rm data},\nu)$ is a degenerate case of a weighted Sobolev norm between characteristic functions of $\mu_{\rm data}$ and $\nu$.

H.2 Algorithms for HM-MMD-PCA

Analogously to the case of MMD-PCA we see that the task (20) is equivalent to:

[TABLE]

and

[TABLE]

Thus, the Algorithm 6 is an adaptation of Algorithm 3 to HM-MMD-PCA.

Again, as in a numerical algorithm for MMD-PCA, at step 1, we search over $q$ given in the form (129). The objective of step 1 can be represented as

[TABLE]

where ${\mathcal{U}}(1,n)$ is the discrete uniform distribution over $\{1,\cdots,n\}$. To apply the stochastic gradient descent methods we need to have an unbiased estimator of $\nabla_{\theta}\Phi(\theta)$ which is equal to

[TABLE]

Thus, a natural estimator of the gradient is:

[TABLE]

where $\{a[s,i,j]\}_{s=\overline{1,4},i=\overline{1,m_{1}},j=\overline{1,s}}\sim^{iid}{\mathcal{U}}(1,n)$ and \{\text{\boldmath\xi}_{i}\}_{i=1}^{m_{2}}\sim^{iid}{\hat{\zeta}}. Overall, we obtain the following Algorithm 7.

Appendix I A numerical alternating scheme for WD-PCA

Let us consider the case $p=1$ and denote $W(\mu,\nu)=W_{1}(\mu,\nu)$. By Theorem 12, the task 66 is equivalent to $\min_{\mu\in{\mathcal{P}}_{k}}W(\mu,\mu_{\rm{data}})$, or to the following task:

[TABLE]

where $I(T_{\mu})=W(\mu,\mu_{\rm{data}})$ if $\mu\in{\mathcal{P}_{k}}$ and $I(\phi)=\infty$, if otherwise. The alternating scheme 1 is designed to solve the penalty form of the problem, i.e.

[TABLE]

which is equivalent to

[TABLE]

where ${\mathcal{S}}_{p}({\mathbb{R}}^{n})\subseteq{\mathcal{S}}({\mathbb{R}}^{n})$ is a set of Schwartz functions that can serve as pdf: $\phi({\mathbf{x}})\geq 0$, $\int_{{\mathbb{R}}^{n}}\phi({\mathbf{x}})d{\mathbf{x}}=1$. A numerical version of the alternating scheme requires additional specifications on a) how to minimize over $\phi$ at step 1, and b) how to estimate $M_{\phi_{t}}$.

I.1 How to minimize over $\phi$?

In the case of WD-PCA, the minimization step of the alternating scheme makes the following:

[TABLE]

where $S_{f}=\sqrt{{\rm O}_{M}}[{\mathbf{x}}f({\mathbf{x}})]$.

For a numerical implementation of that step we need to choose some family of functions that is dense in ${\mathcal{S}}_{p}({\mathbb{R}}^{n})$ (or, rich enough to approach the solution $\mu^{\ast}$). Following the tradition of GAN research let us assume that the family is given in the following form333If $\mathcal{H}\subseteq{\mathcal{S}}({\mathbb{R}}^{n})$ is not satisfied, then we can choose $\mathcal{H}_{\epsilon}=\left\{\phi_{\theta}\ast G_{\epsilon}^{n}|\theta\in\Theta\right\}$ for a very small $\epsilon$.:

[TABLE]

where $\{g_{\theta}|\theta\in\Theta\}$ is a parameterized family of smooth functions (usually, a neural network) and $p({\mathbf{z}})$ is some fixed distribution (usually, the Gaussian distribution). Following [52], we make the assumption 1. In a numerical algorithm we need an access to a procedure that samples according to $\phi_{\theta}({\mathbf{x}})$, not the function itself.

Assumption 1**.**

$\|g_{\theta^{\prime}}({\mathbf{z}}^{\prime})-g_{\theta}({\mathbf{z}})\|\leq L(\theta,{\mathbf{z}})(\|\theta^{\prime}-\theta\|+\|{\mathbf{z}}^{\prime}-{\mathbf{z}}\|)$* where*

[TABLE]

Thus, instead of solving 144 we solve:

[TABLE]

taking into account that $\phi_{t-1}\in{\mathcal{H}}$.

The Kantorovich-Rubinstein duality theorem gives us that:

[TABLE]

which turns 144 into the following minimax task:

[TABLE]

In practice, we choose a family of functions $\mathcal{L}=\{f_{w}|w\in\mathcal{W}\}$ and internal maximization is made over $w\in\mathcal{W}$ with an additional penalty term that penalizes a violation of the Lipschitz condition: $\forall{\mathbf{x}}:\|f_{{\mathbf{x}}}\|\leq 1$.

A family of minimax algorithms for the minimization of $W(\phi_{\theta},\mu_{{\rm emp}})$ was developed in a series of papers [52, 53, 54]. The standard minimax scheme that gained popularity in GAN literature iterates two steps: a) $n_{{\rm iter}}$ times make a gradient ascent over $w\in\mathcal{W}$, b) make a gradient descent over $\theta$. The task 149 can be viewed as a Wasserstein GAN with an additional regularization term $\lambda T(\theta)$ where $T(\theta)=\|S_{\phi_{\theta}}-P_{t-1}S_{\phi_{\theta_{t-1}}}\|^{2}$. To adapt these algorithms to the minimization of our function, we only need to have an unbiased estimator of the gradient $\frac{\partial T}{\partial\theta}$. This estimator is needed for the generator to make its gradient descent step. The discriminator’s part of the algorithm (in which we maximize over Lipschitz functions $f_{w}$) can be set in a standard fashion — we choose [55]’s version, in which the term $\max\{0,\|\frac{\partial f_{w}}{\partial{\mathbf{x}}}(\xi{\mathbf{x}}+(1-\xi)g_{\theta}({\mathbf{z}}))\|-1\}^{2}$ enforces Lipschitz condition (see step (*) of the Algorithm 8).

I.2 How to estimate $\frac{\partial T}{\partial\theta}$ and $M_{\phi_{\theta_{t}}}$?

Another important aspect of the numerical algorithm is the complexity of estimating the matrix $M_{\phi_{\theta_{t}}}$ at step (**). The following theorem shows that we only need to sample ${\mathbf{z}}\sim p$ a sufficient number of times to estimate $\frac{\partial T}{\partial\theta}$ and $M_{\phi_{\theta_{t}}}$.

Theorem 16**.**

If $\phi_{\theta}$ is pdf of the random vector $g_{\theta}({\mathbf{z}})$, ${\mathbf{z}}\sim p({\mathbf{z}})$, then

[TABLE]

where

[TABLE]

and RHS is well-defined.

To prove the theorem we need the following lemma first.

Lemma 5**.**

$\|S_{\phi}-PS_{\psi}\|^{2}={\mathbb{E}}_{{\mathbf{x}},{\mathbf{y}}\sim\phi}({\mathbf{x}}\cdot{\mathbf{y}})M({\mathbf{x}},{\mathbf{y}})+{\mathbb{E}}_{{\mathbf{x}},{\mathbf{y}}\sim\psi}({\mathbf{x}}\cdot P{\mathbf{y}})M({\mathbf{x}},{\mathbf{y}})-2{\mathbb{E}}_{{\mathbf{x}}\sim\phi,{\mathbf{y}}\sim\psi}({\mathbf{x}}\cdot P{\mathbf{y}})M({\mathbf{x}},{\mathbf{y}})$.

Proof of lemma.

[TABLE]

∎

Proof of Theorem 16.

Using lemma 5 we have

[TABLE]

The second term does not depend on $\theta$. Therefore,

[TABLE]

where

[TABLE]

If ${\mathbb{E}}_{{\mathbf{z}},{\mathbf{z}}^{\prime}\sim p}\frac{\partial\Xi(\theta,{\mathbf{z}},{\mathbf{z}}^{\prime})}{\partial\theta}$ is well-defined (the proof of sufficiency of that condition is similar to the proof of Theorem 3 from [52]), then, using Leibniz integral rule, we obtain

[TABLE]

The fact that

[TABLE]

is obvious from the definition $M_{\phi_{\theta}}={\mathbb{E}}_{{\mathbf{x}},{\mathbf{y}}\sim\phi_{\theta}}{\mathbf{x}}{\mathbf{y}}^{T}M({\mathbf{x}},{\mathbf{y}})$. ∎

I.2.1 Definition of ${\mathcal{H}}$

Specifically, for robust PCA/outlier pursuit applications, we define $\phi_{\theta}({\mathbf{x}})$ as a probability density function of the random vector ${\mathbf{a}}+{\mathbf{b}}$, where ${\mathbf{a}}$, ${\mathbf{b}}$ are independent and ${\mathbf{a}}$ is the $i$-th column of matrix $\theta_{1}\in{\mathbb{R}}^{n\times N}$ (where $i\sim{\mathcal{U}}(1,N)$ is sampled uniformly from $\{1,\cdots,N\}$), ${\mathbf{b}}=g_{\theta_{2}}({\mathbf{c}})$, ${\mathbf{c}}\sim{\mathcal{N}}({\mathbf{0}},I_{n})$ and $g_{\theta_{2}}:{\mathbb{R}}^{n}\rightarrow{\mathbb{R}}^{n}$ is a neural network with weights $\theta_{2}$. Thus, $\theta=(\theta_{1},\theta_{2})$. It can be checked that ${\mathcal{H}}$, defined in this way, satisfies the Assumption 1. We specifically introduce the random vector ${\mathbf{a}}$ here because, according to Theorem 12, the ultimate solution of the problem corresponds to $\theta_{1}=Y$ and ${\mathbf{b}}={\mathbf{0}}$. This guarantees that the solution is approachable from set ${\mathcal{H}}$.

Appendix J A numerical alternating scheme for SDR-ORF

For a binary classification case, given a labeled dataset $\{({\mathbf{x}}_{i},y_{i})\}_{i=1}^{N}$, ${\mathbf{x}}_{i}\in{\mathbb{R}}^{n},y_{i}\in{\mathcal{C}}$, ${\mathcal{C}}=\{0,1\}$ we formulate the sufficient dimension reduction problem as the minimization task

[TABLE]

where $L(c,y)=-c\log(y)-(1-c)\log(1-y)$.

We apply the alternating scheme in the dual space (Algorithm 3) to this task. We set $M({\mathbf{x}},{\mathbf{y}})=\zeta({\mathbf{x}}-{\mathbf{y}})$, where $\hat{\zeta}$ is a strictly positive probability density function. A numerical version of the scheme is given below (Algorithm 9).

At every iteration $t=1,\cdots,T$ of the Algorithm 3 we solve the task (in our case $\tilde{I}=J$)

[TABLE]

In a numerical version of the algorithm we assume that ${\hat{\phi}}$ is given as a neural network $f_{\theta}$, i.e. our task becomes

[TABLE]

The gradient of the function \Phi(\theta)=J(f_{\theta})+\tilde{\lambda}{\mathbb{E}}_{\text{\boldmath\xi}\sim\hat{\zeta}}\|\frac{\partial f_{\theta}}{\partial{\mathbf{x}}}(\text{\boldmath\xi})-P_{t-1}\frac{\partial f_{\theta_{t-1}}}{\partial{\mathbf{x}}}(\text{\boldmath\xi})\|^{2} equals

[TABLE]

That is why $\nabla_{\theta}L$ (given to Adam optimizer in the gradient descent loop) in the Algorithm 9 is an unbiased estimator of $\frac{\partial\Phi(\theta)}{\partial\theta}$. Thus, in the “while loop” we find optimal ${\hat{\phi}}_{t}=f_{\theta_{t}}$.

According to Algorithm 3, the next goal is to estimate

[TABLE]

It is easy to see that

[TABLE]

From the last we see that the matrix $M_{t}$ can be estimated by sampling \text{\boldmath\chi}\sim\hat{\zeta} a sufficient number of times (the parameter $l$ in our algorithm). All the rest is identical to Algorithm 3.

The regression version of the algorithm can be obtained by setting $L(c,c^{\prime})=(c-c^{\prime})^{2}$. Implementations for different databases can be found at github.