Reducing the dimensionality of data using tempered distributions

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.


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 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 (R n ) are orthogonal k-frames O ∈ R n×k whose column space is the k-dimensional space onto which a dataset is projected. Different loss functions on V k (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 arXiv:1903.05083v2 [math.ST] 7 Nov 2022 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 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 δ, and δ 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 L, finding an optimal distribution supported in 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 {x i } N i=1 ⊆ R n naturally corresponds to the distribution 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 emp (x) by q(x), where q(x) is a distribution whose density is supported in a k-dimensional linear subspace L ⊆ R n . Note that a function whose density is supported in some low-dimensional subset of R n is not an ordinary function. An exact definition of a set of such distributions, denoted by G k , is given in Section 3. To formulate an optimization task we additionally need a loss D(p emp , q) that measures the distance between the ground truth p emp and a distribution q, that we search for. Thus, in our approach, the UDR problem is defined as I (q) = D (p emp , q) → min q∈G k (2) under the condition that q(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) {w 1 , · · · , w k } ⊆ R n , such that the regression function can be searched in the form g(w T 1 x, · · · , w T k 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(w T 1 x, · · · , w T k 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 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 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 (q) + λ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 : R n × R n → C.
First we note that M defines a billinear form on pairs of (possibly, generalized) functions by f |M |g = R n ×R n f (x) * M (x, y)g(y)dxdy. On a properly defined space of (generalized) functions, the billinear form ·|M |· 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. dist M (f, g) = ( f |M |f − g|M |g − 2 f |M |g ) 1/2 coincides with the maximum mean discrepancy metric [20]. We define R(q) as where M q = Re [ x i q(x)|M |x j q(x) ] i,j=1,n and λ 1 (M q ) ≥ λ 2 (M q ) ≥ · · · 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 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 λ → +∞ 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 where S q is a linear operator between a suitable space H and R n that itself depends on q linearly, and this automatically gives us that R(q) = min S S q − S 2 * where the minimum is taken over all operators between H and R n of rank k and · * is a suitable norm on the space of bounded linear operators from H to R n .
Then, a natural idea to solve the task (4) is to present it as the joint minimum min q min S:rank S=k I (q) + λ S q − S 2 * and to minimize the objective over q and over S of rank k in an alternating fashion, i.e. q l+1 = arg min 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(x, y) = (x · y)H(x, 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 {x i } N i=1 we build a data matrix X = [x 1 , · · · , x N ], a Gram matrix G = [H(x i , 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 G k , and an image of G k under the Fourier transform, denoted F k . In Section 4 we formulate some UDR/SDR problems as optimization tasks over G k /F k . Instead of searching directly in a set of generalized functions, 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, λR(f ). Using a kernel M (x, y), Theorem 4 characterizes generalized g ∈ G k as such g for which the matrix of properly defined integrals M g = Re R n ×R n x i y j g(x) * M (x, y)g(y)dxdy i,j=1,n is of rank k. In Section 6 we suggest a method for solving min φ I(φ) + λR(φ) 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.

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 (R n ). The main advantage of the Stiefel manifold over G k is that its elements are finite-dimensional. Because a distribution from G k is an infinite-dimensional object, an optimization over G k requires additional constructions to turn it into a finite-dimensional task. Both an optimization over G k and over V k (R n ) is typically hard: for a final point, at best one can guarantee that it is a local extremum. Promising aspects of G k are: a) G k allows to formulate a new class of objectives naturally on it, b) local extrema on G k substantially differ from local extrema on V k (R n ), because a local search over 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 . 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].

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 S(R n ), is a space of infinitely differentiable functions f : R n → C such that ∀α, β ∈ N n , sup x∈R n |x α D β f (x)| < ∞, and equipped with standard topology. Its dual space is denoted by S (R n ) and is equipped with weak topology. For a tempered distribution T ∈ S (R n ) and φ ∈ S(R n ), T, φ denotes T (φ). Thus, for a sequence {f s } ⊆ S (R n ) and f ∈ S (R n ), For Ω ⊆ S(R n ), Ω denotes the sequential closure of Ω with respect to standard topology of S(R n ). For Ω ⊆ S (R n ), Ω * denotes the sequential closure of Ω with respect to weak topology of S (R n ). For ψ ∈ S(R n ), T ∈ S (R n ), the convolution is defined as a tempered distribution ψ * T such that ψ * T, φ = T,ψ * φ whereψ(x) = ψ(−x). If T ∈ S (R n ) and a function ψ is such that ψ(x)φ(x) ∈ S(R n ) whenever φ(x) ∈ S(R n ), then the multiplication ψT is defined by ψT, φ = T, ψφ . Given a measure µ, by L 2,µ (R n ) we denote the complex L 2 -space with the inner product u, v L2,µ = u(x) * v(x)dµ. The induced norm is then u L2,µ = u, u L2,µ . If dµ = p(x)dx, then L 2,µ is denoted by L 2,p . A set of infinitely differentiable functions in R n is denoted by C ∞ (R n ). A set of infinitely differentiable functions with compact support in R n is denoted by C ∞ c (R n ). If T is a topological space, then a subset S ⊆ T is said to be dense in T if the sequential closure of S is equal to T . For a square matrix A, Tr(A) denotes its trace and for an arbitrary matrix, A F def = Tr(A T A). The identity matrix of size n is denoted by I n . The notation f ∝ g means f = cg where c is some universal constant.

Basic function classes
To formalize distributions supported in a k-dimensional subspace, we need a number of standard definitions [31]. For φ 1 ∈ S(R k ) and φ 2 ∈ S(R n−k ), their tensor product is the function , is called the tensor product of S(R k ) and S(R n−k ). For g 1 ∈ S (R k ) and g 2 ∈ S (R n−k ), their tensor product is defined by the following rule: there is only one distribution g 1 ⊗ g 2 ∈ S (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 an ordinary function, then g ⊗ δ n−k can be understood as a generalized function whose density is concentrated in a subspace {x ∈ R n |x i = 0, i > k} and equals f (x 1:k ). It can be shown that the distribution acts on φ ∈ S(R n ) in the following way: Now to generalize the latter definition to any k-dimensional subspace we have to introduce a change of variables in tempered distributions.
Let g ∈ S (R n ) and U ∈ R n×n be an orthogonal matrix, i.e. U T U = I n . Then, g U ∈ S (R n ) is defined by the rule: . Now, we define classes of tempered distributions: and where O(n) = {U ∈ R n×n | U T U = I n }. The first two classes are related as: The last two classes are isomorphic under the Fourier transform.
Proof. Let us prove first that if g = T f ⊗ δ n−k , then F[g] = T r , where r(x) =f (x 1:k ), x ∈ R n . For that we have to prove that F[g], φ = T r , φ for any φ ∈ S(R n ). Indeed, Let us calculate the image of G k under the Fourier transform. It is easy to see that for any g ∈ S (R n ), φ ∈ S(R n ) and orthogonal U ∈ R n×n we have: where r (x) = r(U x) =f (U k x) and U k ∈ R k×n is a matrix consisting of first k rows of U . Thus, T r ∈ F k .
Let us show that by varying f ∈ S(R k ) and U in the expressionf (U k x) we can obtain any function from F k . For this, it is enough to show that F k is equivalent to the following set of functions: where U k = (P P T ) −1/2 P and g (y) = g((P P T ) 1/2 y). By construction, U k U T k = I k and g ∈ S(R k ). Thus, Q = F k . Therefore, F[G k ] = F k , and from the bijectivity of the Fourier transform we obtain For any collection which is a linear space over R. The set G k has the following simple characterization: Informally, the theorem holds because any linear dependency α 1 x 1 T + · · · + α n x n T = 0 over R implies that if α 1 x 1 + · · · + α n x n = 0, then T = 0. This is equivalent to a statement that the support of T is concentrated on a subspace α 1 x 1 +· · ·+α n x n = 0. If dim span R {x 1 T, x 2 T, · · · , x n T } ≤ k, then one can find n−k such dependencies, which means that the support of T is k-dimensional.
Let B(R n ) denote the Borel sigma-algebra on R n and P denote a set of all Borel probability measures on R n . Let us now define i.e. P k is a set of probability measures with all probability concentrated in some subspace span(v 1 , · · · , v k ) whose dimension is not greater than k. It is easy to see that T µ ∈ G k for any µ ∈ P k .

Examples of LDR formulations
Maximum mean discrepancy PCA (MMD-PCA) Let K : R n × R n → R be a continuous Mercer kernel, and H K be a reproducing kernel Hilbert space (RKHS) defined by K. The kernel K(x, y) defines the so-called kernel embedding of probability measures φ [32]: The Maximum Mean Discrepancy (MMD) distance [20] is defined as the distance induced by metrics on H K , i.e. for two probability measures µ, ν ∈ P, Let x 1 , · · · , x N ∈ R n be the dataset of points. This dataset defines the empirical probabilistic measure µ data that corresponds to the tempered distribution We shall study a method concurrent to PCA that is based on solving the following problem: i.e. we shall attempt to approximate the empirical probabilistic measure µ data with another probabilistic measure ν which is supported in some k-dimensional subspace of R n . To our knowledge, the task (17) has not been yet considered in the research field of LDR.
2h 2 (2πh 2 ) n/2 is the radial Gaussian kernel on R n and K(x, y) = (k * k)(x − y) = G n 2h (x). For such a kernel, we have where ψ(µ) = k(x − y)dµ(y) is just a smoothing of the distribution µ via the Weierstrass trasform.
In this example, as h → +0, the optimal measure ν * = arg min ν∈P k ψ(µ data ) − ψ(ν) L2(R n ) is supported in a k-dimensional subspace that contains the largest possible number of points from {x 1 , · · · , x N }. Khachiyan demonstrated [33] that the following problem is NP-hard: given {x 1 , · · · , x N } ⊆ R n , find an n − 1-dimensional subspace of R n that contains at least (1 − ε)(1 − 1 n )N points from the dataset. This indicates that in the regime h → +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: where are corresponding moments. The positive parameters λ 1 , λ 2 , λ 3 , λ 4 are chosen to fix the relative importance of the mean, the co-variance, the coskewness and the co-kurtosis.
Thus, we will be interested in the following optimization task (analogous to 17): If we set λ 2 = 1 and λ 1 = λ 3 = λ 4 = 0, then the solution of the task (20) coinsides with the solution of the classical PCA. Let us briefly demonstrate that. Let X = [x 1 , · · · , x N ] be the data matrix, Y = [y 1 , · · · , y N ] be the SVD of X truncated at k-th term, σ i (X) be an ith singular value of X. By µ pca we denote a probabilistic measure concentrated in points is of rank k. Therefore, by Eckart-Young- 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(x, y) = 4 s=1 λs n s (x · 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 (R n , · ) be a Banach space and p ≥ 1. Between any two Borel probability measures µ, ν on R n with x p dµ < ∞ and x p dν < ∞ the pth Wasserstein distance is: where Π(µ, ν) is a set of all couplings of µ and ν. The Wasserstein distance defines another version of LDR problem: In the B one can find proofs that in the case of l 1 norm x = 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 ≥ 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
where x i ∈ R n , y i ∈ C (C is a finite set of classes for a classification, or R for a regression problem), the sufficient dimension reduction problem can be informally described as a problem of finding vectors w 1 , · · · , w k ∈ R n such that conditional distributions satisfy p(y|w T 1 x, · · · , w T k x) ≈ p(y|x) (possibly, under some additional assumptions on the form of p(y|x)).
We formulate the SDR-ORF problem as an optimization task: The object f : R n → 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: and By requiring f ∈ F k , we assume that the regression function f satisfies (for k fixed in advance): Thus, given an input x, an output of f depends on the projection of x onto span(w 1 , · · · , w k ). The set span(w 1 , · · · , w k ) is called the effective subspace. In J we describe an algorithm for the SDR-ORF problem.

Reduction of the optimization problem to ordinary functions
The central problem that our paper addresses is the optimization of an objective function over G k ? In this section we suggest an approach based on penalty functions and kernels.

The definition of the penalty function
In this subsection we introduce a penalty function R(f ). Let M : R n × R n → C 1 be some bounded function such that

One can build a Gram matrix from the collection of functions
Definition 1. Let A ∈ R n×n be a positive semidefinite matrix with eigenvalues λ 1 ≥ λ 2 ≥ · · · ≥ λ n (with counting multiplicities). Then, the Ky Fan Let 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∈[n] to be of close to a rank k matrix. By Theorem 3, the condition dim span R {x 1 f, x 2 f, · · · , x n f } ≤ k implies f ∈ G k , therefore, we enforce f to be close to some function from G k .
In the next section we will justify the latter informal logic by reducing the optimization over G k to the optimization over ordinary functions with the penalty function R(f ).

Proper kernels
For any operator O between spaces H 1 and H 2 , we denote its range by is a properly defined, strictly positive and self-adjoint operator, Note that the latter definition implies that M (y, It is captured by the following lemma: Proof. Verification of the first three conditions is easy, so we only check the fourth condition. Let us denote linear Besides the Gaussian kernel the lemma also captures a case of the Laplace kernel ζ(x) = e −|x| . It is well-known that the Fourier tranform of the Laplace kernel is the Poisson kernel: (which is also proper).
For I : G k ∪ S(R n ) → R + , it is natural to reduce the optimization task over tempered distributions to an optimization task over ordinary functions with a penalty term R, where we assume that the set of functions F is rich enough to approximate weakly solutions of (29), i.e. F * ⊃ G 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 } ⊂ F is such that (30) for arbitrarily large values of the regularization parameter), then there exists a growing subsequence {n k } such that T fn k → * T (weakly) where T is a solution of (29).
We make a thorough theoretical analysis of the case F = S(R n ). If to formulate in a simplified way, for the last property to hold, the sequence Tr(M fn ) should be bounded. Details on the conditions under which this reduction holds can be found in the following subsection.

Regular solutions and reduction theorems for
, Lim s→∞ f s denotes a set of points f ∈ S (R n ), such that there exists a growing sequence {s i } ⊆ N and lim i→∞ f si = f .
For I : G k ∪ S(R n ) → 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 ∈ G k from a set S(R n ) is continuous, i.e. for any sequence Let us introduce the notion of a regular solution both for (29) and (30). Let Definition 3. Any f ∈ Arg min In other words, B k formalizes a set of distributions from G k , that can be approached through sequences  (30), because when approaching a non-regular solution we are unable to guarantee a bounded behaviour of M f (and of R(f )). Arg min Theorem 7 (Reduction theorem). If M is a proper kernel, Arg min Suppose that we now solve a sequence of problems (30) and find {f s } ∞ 1 . According to Theorems 5 and 6, the following are potential scenarios: (1) Tr(M fs ) blows up and the convergence is not guaranteed. This situation can be avoided by controlling Tr(M f ) in an optimization process. In practice, when f has a parameterized form, this can be done by bounding parameters.
If Tr(M fs ) does not blow up, we still have two subcases: This implies a positive outcome to approach (30) to the optimization problem, Problem (29).
(2.2) Lim s→∞ T fs = ∅. This exotic situation can happen only if a sequence T fs leaves any sequentially compact subset of S (R n ). Bounding parameters also tackles this case.
Let us now concentrate on the task (30) and describe the alternating scheme for its solution.

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(φ) = M φ n−k depends on M φ in a concave way. It can be shown that the dependence of R(φ) on φ is both non-convex and non-concave, i.e. we deal with a non-convex optimization task.
Let B(H 1 , H 2 ) denote a set of bounded linear operators between Hilbert spaces H 1 and It is easy to see that any O ∈ B(L * 2 (R n ), R n ) can be given by formula: Recall that for a Mercer kernel can be properly defined [40]. Thus, O M is self-adjoint. For any complex-valued function f such that Tr M f < ∞ let us introduce a linear operator S f : L * 2 (R n ) → R n by the following rule: and the minimum is attained at 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: and therefore: Thus, we have checked that S f is bounded.
It is easy to see that the following operator satisfies the latter identity: Since the adjoint is unique, Thus, Let u 1 , · · · u n be orthonormal eigenvectors of M f = S f S † f and λ 1 ≥ · · · ≥ λ n > 0 be corresponding nonzero eigenvalues. For σi . Vector v i corresponds to a pair of functions It is easy to see that v 1 , · · · v n is an orthonormal basis in Im S † f , and S † f can be expanded in the following way: and therefore, SVD for S f is By the Eckart-Young-Mirsky theorem (see Theorem 4.4.7 from [41]), an optimal S in min where Given the new representation R(f ) = min Thus, it is natural to view the Task (30) as a minimization of I(φ) + λ S φ − S 2 * over two objects: f ∈ F and S ∈ B(L * 2 (R n ), R n ) : rank S ≤ k. The simplest approach to minimize a function over two arguments is to optimize alternatingly, i.e. first over f , and then over S : rank S ≤ k, and so on. Theorem 8 gives that the minimization over S is equivalent to the truncation of SVD(S f ) at the k-th term. This idea, that we dub the alternating scheme (AS), is described in Algorithm 1.

Algorithm 1
The alternating scheme (AS) for (30) The alternating algorithm 1 allows for a reformulation in the dual space. By this we mean that in Algorithm 1 we substitute φ t for the original φ t . If the primal Algorithm 1 deals with operators S φ , S φt−1 , the dual version deals with vectors of functions G σ ∂ φ ∂x , G σ ∂ φt−1 ∂x . Details of the dual algorithm can be found in E.
The objective I(f ) + λ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(x, y) = (x · y)H(x, 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 (x) = 1 N N i=1 δ n (x − 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).
Algorithm 2 An approximate algorithm for (17) where K(x, y) = (x · y)H(x, y) 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.

Sketch. Let us apply Theorem 8 to the kernel H and the function f . Recall that an element
, the representation of Theorem 8 gives us The restriction rank S ≤ k is equivalent to can be properly defined and is in L 2 (R n ). Therefore, For any f ∈ G k one can set and search over all possible f ∈ G k in the minimization operator. Thus, Let µ data be a uniform distribution over {x i } N i=1 and d MMD be the MMD distance induced by K(x, y) = (x · y)H(x, y). The last theorem can be applied to a smoothed empirical distribution f ε (x) = 1 N N i=1 G n ε (x − x i ) and then, we can send ε → 0. All the more, the inequality will be satisfied if we search over µ ∈ P k due to T µ ∈ G k . Thus, where λ 1 ≥ · · · ≥ λ n are eigenvalues (counting multiplicities) of Thus, ( n i=k+1 λ 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.
The following corollary is straightforward from the last theorem. For the case when H is the Gaussian kernel the situation is slightly trickier.

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)  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(x, y) = 1 √ (8πh 2 ) n e − 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 M, which is bent globally, goes through the origin O and has a large curvature at O (see Fig. 1a). 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 : R n−1 → R, such that f (0) = 0 and generate points in the following way: points x 1 , x 2 , · · · , x N ∼ [−10, 10] n−1 are sampled uniformly, after calculation of y i = f (x i ) we add some noise: z i = (x i , y i )+ i , i ∼ N (0, 0.01I n ). Both PCA and MMD-PCA are applied to the dataset (first 3 pictures on Figure 1a). 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 → 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, δ = 0.05 (or 0.1), k = 10 and generate random matrices A ∈ R N (1−δ)×k , B ∈ R n×k whose entries are iid as N (0, 1). Then, columns of the matrix BA T ∈ R n×N (1−δ) (whose rank is ≤ k) are concatenated with columns of the matrix C ∈ R n×N δ : X = concat(BA T , C) ∈ R n×N . The entries in C are either iid as N (0, 1) (case I) or N δ copies of the same vector whose entries are iid as N (0, 1) (case II). Let X = [x 1 , · · · , x N ], i.e. columns of X are the data points. Thus, N (1 − δ) columns of BA T lie in a k-dimensional subspace of R n and N δ 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 δ and λ is shown in Figure 1b. 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 * −P F , where P * = lim t→∞ P t is the final projection operator (e.g. P 20 in practice), depends on the parameter σ of the kernel M = G n σ (bandwidth). It is natural to expect the quality of the solution P * to degrade as σ → +∞ (this corresponds to M (x, y) → 0), and, less trivially, as σ → 0 (this corresponds to M (x, y) → δ n (x − y)). As the right plot on Figure 1b shows, for the HM-MMD-PCA, the solution P * is close to the correct P if the bandwidth σ 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 σ 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 σ = 0.8 2 and λ = 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 R 2 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 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 32256dimensional vectors that lie on some low-dimensional subspace L of 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(x, y) = (x · y)e − x−y 2 n (which we simply call Gauss). Our experiments showed that behaviour of PCA and Gauss are quite similar if the dimension of L is 5, though Gauss removes more shadows and preserves more details of an original image if the dimension of 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 L. We assume that the dimension of L is 5. We recover L using PCA and the Algorithm 2 for kernels K(x, y) = (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 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 ( (L−1) 2 M SE ) 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 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 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).  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 ∼ 10 6 − 10 7 . Also, the Algorithm 2 can be used for an initialization of the alternating scheme.

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 , σ, and the penalty parameter λ. 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.

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

A Proofs for section 3
A.1 Proof of Theorem 1: given for completeness Proof. The inclusion G k ⊆ G k * follows from a well-known fact that S(R k ) is dense in S (R k ). I.e. for any f ∈ S (R k ) one can always find a sequence {f i } ⊆ S(R k ) such that T fi → * f . Therefore, for any Since G k ⊆ G k , to prove G k = G k * it is enough to show that G k is sequentially closed.
We need a simple fact from a theory of distributions.
Proof of Lemma. Schwartz space S(R n ) is a Fréchet space, therefore the Banach-Steinhaus theorem applies to For any T ∈ S (R n ) and ψ ∈ S(R n−k ), let us define We need to prove that f ∈ G k . Since a set of orthogonal matrices is compact, then one can always find a subsequence {U ni } such that U ni → U . Since (f ni ⊗ δ n−k ) Un i → * f and φ(U ni x) → φ(U x) (for any fixed φ ∈ S(R n )), using lemma 2 we obtain: Thus, we have f ni ⊗ δ n−k → * f U T . From the last we see that f ni → * f ψ U T where ψ is such that ψ(0) = 1. Therefore,

A.2 Proof of Theorem 3
Proof of Theorem 3 (⇒). Let us prove that from T = (f ⊗ δ n−k ) U , f ∈ S (R k ), U T U = I n it follows that dim span R {x 1 T, x 2 T, · · · , x n T } ≤ k.
It is easy to see that Thus, we have n − k orthogonal vectors, u k+1 , · · · , u n , such that Using standard linear algebra we obtain there are at most k distributions x i1 T, · · · , x i k T, k ≤ k that form a basis of To prove the second part of theorem we need the following lemma. Lemma 3. If T ∈ S (R n ) is such that y i T = 0 for any i > k, then T ∈ G k .
Proof of lemma. Recall from functional analysis, for f ∈ S (R n ), the tempered distribution ∂f ∂xi is defined by the condition ∂f ∂xi , φ = − f, ∂φ ∂xi . Once the Fourier transform is applied, our lemma's dual version is equivalent to the following formulation: if ∂f ∂xi = 0, i > k, then f ∈ F k * . Let us prove it in this formulation.
Recall that a set of infinitely differentiable functions with a compact support is denoted by C ∞ c (R). Suppose φ ∈ S(R n ) and p ∈ C ∞ c (R) are chosen in such a way that It is easy to see that for any α ∈ N n−1 , α ∈ N, β ∈ N n−1 , β ∈ N we have (at least one derivative over x i is present): The terms Analogously (when not a single derivative over x i is present): The second term is 0 when The latter is bounded because The first term is 0 when x i ≥ B and for x i < B it satisfies: The latter is also bounded, since ∂x β is bounded and r ∈ S(R n ). Therefore, ∂f ∂xi = 0 implies: Since this sequence of arguments works for any i > k, we can apply them sequentially to initial φ ∈ S(R n ) w.r.t.
x k+1 , ..., x n . Thus, for any p k+1 , ..., p n ∈ C c (R) such that ∞ −∞ p i (y i )dy i = 1 we obtain: Moreover, since C ∞ c (R) is dense in S(R), we can assume that p k+1 , ..., p n ∈ S(R). For the inverse Fourier transform T = F −1 [f ] the latter condition becomes equivalent to: for any p k+1 , ..., p n ∈ S(R) such that p i (0) = 1. Let us define p i (x i ) = e −x 2 i . It is easy to check that T = g ⊗ δ n−k where g ∈ S (R k ), g, ψ = T, e −|x k+1:n | 2 ψ(x 1:k ) for ψ ∈ S(R k ). Thus, T ∈ G k and lemma is proved.
Thus, there exist at least n − k orthonormal vectors v k+1 , · · · , v n , such that [x 1 T, · · · , x n T ]v i = 0. Therefore, [x 1 T, · · · , x n T ]v i = (v T i x)T = 0. Let us complete v k+1 , · · · , v n to form an orthonormal basis of R n : v 1 , · · · , v n . Let us define a matrix V = [v 1 , · · · , v n ]. It is easy to see that: ( Since for i > k we have (v T i x)T = 0, then x i T V = 0. Using lemma 3 we obtain T V ∈ G k . Therefore, (T V ) V T = T ∈ G k . Theorem proved.

B Structure of WD-PCA
Recall that (R n , · ) is a Banach space and p ≥ 1. Now, let us consider an optimization problem: for a given 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 Proof. Let us prove first that inf µ∈P k W p (µ data , µ) Let π be a uniform distribution over {(x i , y i )} N i=1 and µ * be a uniform distribution over The support of µ * is kdimensional, because rank(Y * ) ≤ k. Thus, we have µ * ∈ P k and inf µ∈P k W p (µ data , µ) ≤ W p (µ data , µ * ) ≤ 1 N p X − Y * p . Now, if we prove the inverse inequality, i.e. inf µ∈P k W p (µ data , µ) ≥ 1 N p X − Y * p , this will imply that inf µ∈P k W p (µ data , µ) = 1 N p X − Y * p and therefore, inf µ∈P k W p (µ data , µ) = W p (µ data , µ * ). This will in the end give us µ * ∈ arg inf µ∈P k W p (µ data , µ).
It is easy to see that W p (µ * t , µ data ) ≤ W p (µ t , µ data ), because µ * t and µ t share the same k-dimensional support L t , but the "transportation of a mass" concentrated in point x i of the empirical distribution µ emp can be most optimally done by just moving it to P t x i (i.e. to the closest point on L t ). Thus, we have inf µ∈P k W p (µ data , µ) ≤ W p (µ data , µ * t ) ≤ W p (µ data , µ t ), and therefore, W p (µ data , µ * t ) − inf µ∈P k W p (µ data , µ) → 0. Since a set of projection operators is compact, one can always extract a subsequence {P ts } ∞ s=1 , such that P ts → P . It is easy to see that µ * ts → µ * * (i.e. W p (µ * ts , µ * * ) → 0) where µ * * is a uniform distribution over {P x 1 , · · · , P x N }. For that distribution we have Thus, the infinum is attained on µ * * .
It is easy to see that W p (µ data , µ * * ) = 1 N p X − P X p . Since rank(P X) ≤ k we obtain W p (µ data , µ * * ) ≥ This completes the proof.
A PREPRINT -NOVEMBER 8,2022 Note that in the case of l 1 norm and p = 1, i.e. x = 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 ≥ 1, this leads to another task: where [s 1 , · · · , s N ] p,2 = ( N i=1 s i p 2 ) 1/p . This task has many applications in mathematics and is known as the subspace approximation problem [39] . Let us first show that f |M |g is defined for any f = (T a ⊗ δ n−k ) U ∈ G k and g = ( We have Let us denote a = a * G k and b = b * G k . It is easy to see that From a well-known property of the Weierstrass transform we have From this we obtain that where M (x 1:k , y 1:k ) = R n−k ×R n−k G n−k (x k+1:n )M (U T x, V T y)G n−k (y k+1:n )dx k+1:n dy k+1:n .
Thus, there exists boundedM (x 1:k , y 1: Further we assume that > 0 is small enough, so that M (x 1:k , y 1:k ) ≤ C = 2 max |M (x, y)|. Now we have: It is well-known (e.g. see Theorem 2.25 from [49] Thus, lim →0 f |M |g exists and f |M |g is defined.
Let us now prove that rank is an orthogonal matrix. By construction, Let us now denote V = [u 1 , · · · , u n ] ∈ R k×n a submatrix of U in which only first k rows of U are present. Then, the latter integral is equal to R k ×R k u T i x 1:k y T 1:k u j g(x 1:k ) * M (V T x 1:k , V T y 1:k )g(y 1:k )dx 1:k dy 1: is the Gram matrix of the collection {x i g(
Proof. W.l.o.g. we can assume that f = T l ⊗ δ n−k , l ∈ S(R k ). By lemma, all entries of M fσ except those of the main k × k minor approach 0 as σ → 0. This means that lim σ→+0 Q( Let v 1 , · · · , v n be unit eigenvectors of M fσ corresponding to the eigenvalues λ 1 ≥ · · · ≥ λ n , P = n i=k+1 e i e T i , then Since R(f σ ) ≤ Q(f σ ), we obtain lim σ→0 R(f σ ) = 0.
D.0.1 Proof of Theorem 5 Proof. Suppose that a sequence {f i } ∞ s=1 ⊆ S(R n ) regularly solves (7) and T ∈ Lim i→∞ f i . W.l.o.g. we can assume that T fi → * T and Tr(M fi ) is bounded and Below we use continuity of I and corollary 2: from which we conclude that For each l, let us define P l as the projection operator to a subspace spanned by first principal components of the matrix M f l , i.e.
Let us treat now M as an operator O M : Since O M is a strictly positive self-adjoint operator, by the Cauchy-Schwarz inequality, we obtain Therefore, for any ψ ∈ Range [O M ] ∩ S(R n ) and i > k we have lim s→∞ y i f ls , ψ = lim s→∞ f ls , y i ψ = 0.
Using lemma 3 and (T V ) V T = T we obtain T ∈ G k . Thus, we proved that T fi → T ∈ G k . Since T ∈ Arg min

D.0.2 Proof of Theorem 6
Proof. Suppose f * ∈ Arg min and (using lemma 2) Due to the continuity of I we have lim i→∞ I(s i and we obtain the needed sequence: lim where Tr M fi is bounded. It remains to check that our sequence regularly solves (7), i.e. lim i→∞ inf f ∈S(R n ) Let us prove the inverse inequality.
Since rsol (I(f ), R(f )) = ∅, there exists {f i } ⊆ S(R n ) such that and a = lim i→+∞ Tf i . From theorem 5 we obtain a ∈ Arg min f ∈G k I(f ).
E The alternating scheme in the dual space for M (x, y) = ζ(x − y) When M (x, y) = ζ(x − y), the alternating scheme 1 allows for a reformulation in the dual space. By this we mean that in Scheme 1 we substituteφ t for the original φ t . If the primal Scheme 1 deals with operators S φ , S φt−1 , the dual version deals with vectors of functions ζ ∂φ ∂x , ζ ∂φt−1 ∂x . The substitution is based on the following simple fact: Theorem 13. If M (x, y) = ζ(x − y), ζ,ζ ∈ C(R n ) and ∀xζ(x) > 0, then there exist constants c 1 and c 2 such that where κ is a constant.

Let us now introduce a vector of functions
Using 105 we obtain (S f ) i = κ ζ ∂f ∂xi , and therefore V f = κ ζ ∂f ∂x . Thus, the expression S φ − P t−1 S φt−1 2 * in the alternating scheme can be rewritten as The matrix M f can also be calculated fromf using the following identity: Let us introduce a functionĨ such thatĨ(f ) = I(f ). Then, we see that all steps in Scheme 1 can be performed witĥ φ t rather than with φ t , using the algorithm 3.
Informally, the dual algorithm works as follows: at each iteration t we compute a functionφ t adapting it to data (the termĨ(φ)) and adapting its gradient field to the rank reduced gradient field of the previousφ t−1 . For a sufficiently large T , it will converge andφ T ≈φ T −1 . Then, the second term in the last step will be approximately equal to Thus, gradients ∂φ T ∂x lie in a k-dimensional subspace col P T −1 . This last property is a characteristic property of functions from 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 etc. This version of the dual alternating scheme will be used for designing numerical algorithms for the Gaussian MMD-PCA and HM-MMD-PCA.

F Proofs for Section 7
Proof of Theorem 10. Let X = [x 1 , · · · , x N ]. Note that H f = XSX T where S = [P (x i · x j )] i,j∈ [N ] . For any U ∈ O(n) and we have Algorithm 3 The alternating scheme in the dual space.
Therefore, w.l.o.g. we can assume that principal components of H f are e 1 , · · · , e n and H f e i = λ i e i , i ∈ [n], where {e i } n i=1 is a canonical basis in R n and λ 1 ≥ λ 2 ≥ · · · ≥ λ n are eigenvalues of H f . From the latter we conclude that For :n ] and z 1:k ∈ R k equals the first k components of z ∈ R n . By construction, The second sum F 2 equals n i=k+1 λ i . Let us compare the first sum, F 1 , with the second, F 2 . F 1 is a sum of positive factors of (E x∼f x α ) 2 where α 1:k = 0, α k+1:n = 0. Let e i ∈ (N ∪ {0}) n be an ith canonical unit vector and α i denote an ith component of α.
Thus, overall we have From T µ − T ν 2 K = d MMD (µ, ν) 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 : R n → R be a smooth function, then Proof of Theorem 11. Let us assume w.l.o.g. that principal components of H f are e 1 , · · · , e n and H f e i = λ i e i , i ∈ [n], where {e i } n i=1 is a canonical basis in R n and λ 1 ≥ λ 2 ≥ · · · ≥ λ n are eigenvalues of H f . From the latter we conclude that f (x) = E y 1:n−k ∼N (0,σ 2 I n−k )f (x 1:k , y 1:n−k ) . From the isometry of the Fourier transform, we have The latter expression decomposes into two terms. The first is The second is we finally obtain . By construction, f ∈ G k and it can be approached by elements of G k w.r.t. norm · K . The statement of theorem directly follows from this observation. From theorems 1 and 2, F −1 [P k ] ⊆ F k * . 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 ∈ F −1 [P], then for any distinct y 1 , · · · , y s ∈ R n the matrix [f (y i − y j )] i,j=1,n is positive semidefinite. Since µ(R n ) = 1, we additionally have f (0) = 1. Let PDF denote the set of all continuous positive definite functions on R n and Thus, the following characterization of F −1 [P k ] becomes evident.

G.2 The dual form of the Gaussian MMD-PCA
Recall that k(x) = G n h (x). Let us define another Gaussian kernel γ(x) = e − h 2 |x| 2 Let p data (x) denote the characteristic function of the random vector X data ∼ µ data . By definition, Using the isometry property of the inverse Fourier transform for L 2 (R n ) and the convolution theorem, we see that Thus, from Theorem 15 we obtain that the task 17 is equivalent to

G.3 Algorithms for the Gaussian MMD-PCA
Let Π k : G k → {1, +∞} and M k : F k → {1, +∞} be simple penalty functions: if φ ∈ M k and M k (φ) = ∞, otherwise Then, the task 17 is equivalent to From the result of the previous section we see that if I(φ) =Ĩ(φ), theñ Thus, the Algorithm 4 is an adaptation of Algorithm 3 to MMD-PCA.

Algorithm 4
The alternating scheme in the dual space for the Gaussian MMD-PCA If the function p 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: where α i > 0 and ,nn ) and u i 's are unconstrained. The number of neurons in a single layer neural network with a cosine activation function, nn, is a hyperparameter. Let us denote parameters {ω i , u i } nn i=1 by θ. It is easy to see the function q θ 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(ω T x)|ω ∈ R n }, is dense in a set of real-valued functions from 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 with stochastic gradient descent methods (in our case, the Adam optimizer) we need to have an unbiased estimator of where z ∼ f denotes that the random vector z is sampled according to the probability density function R n f (x)dx . Thus, a natural estimator of the gradient is where The last important issue with the practical numerical algorithm is the calculation of M t at step 2. By construction, In practice we sample χ 1 , · · · , χ l ∼ζ and estimate M t as The details of the numerical algorithm 5 are given below. In all our experiments with MMD-PCA we setζ = γ 2 .
which is equivalent to W (φ, µ data ) + λR(φ) → min where S p (R n ) ⊆ S(R n ) is a set of Schwartz functions that can serve as pdf: φ(x) ≥ 0, R n φ(x)dx = 1. A numerical version of the alternating scheme requires additional specifications on a) how to minimize over φ at step 1, and b) how to estimate M φt .

I.1 How to minimize over φ?
In the case of WD-PCA, the minimization step of the alternating scheme makes the following: For a numerical implementation of that step we need to choose some family of functions that is dense in S p (R n ) (or, rich enough to approach the solution µ * ). Following the tradition of GAN research let us assume that the family is given in the following form 3 : where {g θ |θ ∈ Θ} is a parameterized family of smooth functions (usually, a neural network) and p(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 φ θ (x), not the function itself.
Thus, instead of solving 144 we solve: taking into account that φ t−1 ∈ H.
The Kantorovich-Rubinstein duality theorem gives us that: which turns 144 into the following minimax task: In practice, we choose a family of functions L = {f w |w ∈ W} and internal maximization is made over w ∈ W with an additional penalty term that penalizes a violation of the Lipschitz condition: ∀x : f x ≤ 1.
A family of minimax algorithms for the minimization of W (φ θ , µ 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 iter times make a gradient ascent over w ∈ W, b) make a gradient descent over θ. The task 149 can be viewed as a Wasserstein GAN with an additional regularization term λT (θ) where T (θ) = S φ θ − P t−1 S φ θ t−1 2 . To adapt these algorithms to the minimization of our function, we only need to have an unbiased estimator of the gradient ∂T ∂θ . 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, ∂fw ∂x (ξx + (1 − ξ)g θ (z)) − 1} 2 enforces Lipschitz condition (see step (*) of the Algorithm 8).
Proof of lemma.