Dimension reduction as an optimization problem over a set of generalized functions
Authors: Rustem Takhanov
Venue:  (2021)
Keywords: Unsupervised dimension reduction
sufficient dimension reduction
alternating scheme
Fourier transform
maximum mean discrepancy
Wasserstein distance
positive definite functions
Bochner’s theorem

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 $p_{\rm{emp}}({\mathbf x})$ by another tempered distribution $q({\mathbf x})$ whose support is in a $k$-dimensional subspace. Thus, our problem is reduced to the minimization of the distance between $q$ and $p_{\rm{emp}}$, $D(q, p_{\rm{emp}})$, over a pertinent set of generalized functions.

This infinite-dimensional formulation allows to establish a connection with another classical problem of data science --- the sufficient dimension reduction problem (SDR). Thus, 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 $R(f)$ that ``forces'' the support of $f$ to be $k$-dimensional.
Then we present an algorithm for minimization of $I(f)+\lambda R(f)$, based on the idea of two-step iterative computation, briefly described as a) an adaptation to real data and to fake data sampled around a $k$-dimensional subspace found at a previous iteration, b) calculation of a new $k$-dimensional subspace. We demonstrate the method on 4 examples (3 UDR and 1 SDR) using synthetic data and standard datasets.

Under review as a conference paper at ICLR 2021
DIMENSION REDUCTION AS AN OPTIMIZATION PROBLEM OVER A SET OF GENERALIZED FUNCTIONS
Anonymous authors Paper under double-blind review
ABSTRACT
We reformulate unsupervised dimension reduction problem (UDR) in the lan- guage of tempered distributions, i.e. as a problem of approximating an empirical probability density function Pemp (x) by another tempered distribution g(x) whose support in k-dimensional subspace. Thus, our problem is reduced to the min- imization of the distance between q and Pemp- D(g, Pemp). over pertinent set of generalized functions.
This infiiitediimensional formulation allows to establish connection with an- other classical problem of data science the sufficient dimension reduction prob- lem (SDR). Thus, algorithm for the first problem induces algorithm for the second and vice versa. In order to reduce an optimization problem over distribu- tions an optimization problem over ordinary functions ative penalty function R(ff) that "forces" the support of f to be k-dimensional. Then we present an algorithm for minimization of I [()) XR(S), based on the idea of two-step iterative computation, briefly described as a) an adaptation to real data and fake data sampled around a k-dimensional subspace found at a pre- vious iteration, b) calculation of a new k-dimensional subspace. We demonstrate the method on 4 examples (3 UDR and SDR) using synthetic data and standard datasets.
introduce nonneg-
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 dimensionality reduction (SDR), maximum au- tocorrelation factors, slow feature analysis and more. In unsupervised dimension reduction (UDR) we are given a finite number of points in I²" (sampled according to some unknown distribution) and the goal is to find "low-dimensional"" affine (or linear) subspace that approximates "the sup- port" the distribution. The study field currently achieved a saturation level at which unifying frameworks for the problem become of special interest Cunningham & Chahramani (2015). An approach that we present in that paper is based on the theory of generalized functions, or tempered distributions (Soboleff] (1936); Schwwarzz (1949). An important generalized function that cannot be represented as an ordinary function is the Dirac delta function, denoted 8, and d" denotes its n- dimensional version.
Any dataset {x:] S R2 naturally corresponds to the distribution Pemp (x) 1 Σi 6x x;) 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 Pemp (x) by gxx), where (xx) is distribution whose density supported in d--dimensional affine subspace A C R... Note that function whose density is supported in some cow-dimensional subset of I2 is not an ordinary function. Exact definitions of such distributions can be found in Section 3 To formulate an optimization task we additionally need loss D ([pemp: g) that measures the distance hetween the ground truth Pemp and a distribution q, that we search for. Thus, in our approach, the UDR problem defined as:
N
I (g) D (pemp. q) min 9
(1)
under the condition that q(x) has k-dimensional support.

---

Under review as a conference paper at ICLR 2021
The SDR problem is tightly connected with the UDR problem. In SDR, given supervised data, the goal is to find the so called effective subspace, defined by its basis vectors (W1. w...... w S R", such that the regression function can be searched in the form g( g(wixx... w x). In Wang et al. (211)) it was shown how a method originally developed for SDR can be turned into an UDR method, i.e. applied unsupervised data, by simply setting
output to be equal to an input. The key observation of our analysis, stated in Theorem 2 is that a class of functions of the form g(w] x..,, w/'x) can be characterized as functions whose Fourier transform is supported in the corresponding effective subspace. In Section 4 we give 3 examples of UDR problems that we cast as [and in the fourth example we formulate SDR as an optimization task with the search space dual to that of UDR. Thus, all 4 examples can be studied within the same optimization framework.
The structure of the paper as follows: in Section B we formally define the search space in Problem [ denoted GR. and an image of 9R under the Fourier transform, denoted Fg. In- stead of searching directly in set of generalized functions, Gu. in Section [ we describe how we substitute ordinary function for distribution in the optimization task at the expence of adding a new penalty term to its objective, AR((). Using a gaussian kernel M (x, y). Theorem characterizes generalized g E Gj such g for which the matrix of properly defined integrals Mg- = Re [J] "xy,g(x) (x.ygg(5 (erllyydxddy..... I.==.m is of rank k. We define R(f) as a squared
Frobenius distance from MM to the closest matrix of rank k. In Section 6 we suggest a method for solving ming (()) + (R()) which we call the alternating scheme. Section is dedicated to experiments with the alternating scheme on synthetic data and standard datasets.
2 PRELIMINARIES AND NOTATIONS
Throughout this paper we use standard terminology and notation from functional analysis. For exact definitions one can address the textbook on the theory of distributions Friedlander & Joshi (1998). The Schwartz space of functions and its dual space are denoted by S(R") and (RR) correspondingly. For a tempered distribution T' E 8' ('(R") and Φ E ((R"), (T, ø) denotes 7 (Ø)) The Fourier and inverse Fourier transforms are denoted by F S' S'(R")
- (R""). For brevity,
denote Fif] by f. If all required conditions are satisfied, an integrable f : R" - C (or, Borel measure μ on IR") is used as the tempered distribution T (or, Tw) where (Ty = SR- f(x)x(x)dx
(or, (Ty. 0) = Sgn ((x)dz). For C (R"),, Ω2 denotes the sequential closure of 52 with respect weak topology of S' ('")"). By La (R)) we denote the L2-space with the inner product: (),e),,
[(x)"" (xx)xxx For & E S(RR), ) + E S (I²",, their convolution and multiplication are denoted by & & and øú correspondingly. For 91 E (IR and 92 E (R)"")),, 91 @ 92 E ('()") denotes their tensor product. For square matrix A, Tr(A) denotes its trace and for arbitrary matrix,
||A||E def ATr(ATA). Identity matrix of size is denoted by In.
3 BASIC FUNCTION CLASSES
An example of a generalized function, whose density is concentrated in a k-dimensional subspace, is any distribution that can be represented as X - def
X & where g E (.*
n k times
If TT, where ʃ R + R is an ordinary function, then g @ 8"- can be understood generalized function whose density is concentrated in a subspace (x E RR = 0, > k} and equals fx It can be shown that the distribution acts à E S(R") in the following way:
0) = S/e fxx1 ((x1==)(((x1+ fn-x)dx1..
(T
Now to generalize the latter definition to any k-dimensional subspace we have to introduce a change of variables in tempered distributions.
Let E (RR)) and U E Ra** be an orthogonal matrix, i.e. UTU = I,. Then, gu E S(R) is (g, t) where (x() = o(*xx), x). If = Ty, the latter definition gives gu = Ty. where f'(x)= = (Ux). Now, we define classes of tempered distributions:
defined by the rule: (gu,()
9 = (( s 8t Dols E S' S'RRR),U Raxn UT =,,}
2

---

Under review as a conference paper at ICLR 2021
9x = ( S R DUUn, =
(3)
F& = [0x) ra)())))
Rex" rank(U))
(4)
The first two classes related as: Theorem 1. 9% Si
The last two classes are isomorphic under the Fourier transform. Theorem 2. F = FR and F [F& = Gz.
For any collection fi.... which is a linear space over IR. The set Gj has the following simple characterization:
8' (IR"), spang [f:]l denotes - AcHiAe E R] S S'RR"),
Theorem 3. For any E ('((")) T Gj if and only if dim spang { xTT,
In T S
Informally, the theorem holds because any linear dependency aj ... an 0 over I implies that a æ + 2n, * 0, then T 0. This equivalent to a statement that the support Tf is concentrated on a subspace 01 +.. t-. +,, = O. If dim spang [x1 222.,... ,.,,,,,,] S k, then one can find n l such dependencies, which means that the support of is k-dimensional.
and P denote set of all Borel probability mea- Let ((R") denote the Borel aamna-l R" sures on R". Let us now define
Ps {u PlEv1, v& E R., AA E B(R"): (aA) = ( span(vi, v&))) (5)
Le. Ps is set of probability measures with all probability concentrated in some subspace span(vi,, "** V&.) whose dimension is not greater than k. It is easy to see that Tç E 9j for any J Pk.
EXAMPLES OF LDR FORMULATIONS
UDR: Maximum mean discrepancy (MMD) Let k(x) 4(2x
247
be the radial gaussian
kernel on R". The kernel k(x) defines the so-called kernel embedding of probability measures Muandet al. 2017)
/ 1(xx yyddu(y)
p 0(M) k = Ey-nul(x- y)
The Maximum Mean Discrepancy (MMD) distance Gretton al. (2012) is defined as the distance induced by metrics on L2 (R"),, i.e. for two probability measures n1, v EVE P:
dMMD (p, v) [1Φ(1) +(v) ((2))))
Let XI, XN E IR" the dataset of points. This dataset defines the empirical probabilistic measure pdata that corresponds to the tempered distribution Tp = + > 8" - x,). We
shall study method concurrent PCA that is based on solving the following problem:
I(v) dMMDD (Paata v) ll@(paata) - <(v)llea(R>) min NEPx
(6)
ie. we shall attempt to approximate the empirical probabilistie measure jlata with another proba- bilistic measure which is supported in some k-dimensional subspace of IR".
UDR: Distance based on higher moments (HM) It is well-known that maximum mean discrep- ancy measures the similarity between characteristic functions of two probability distributions in the O (H) neighboourhoo the origin. Another approach to measure the similarity of two distributions based on the difference between moments:
dHM(U, v/²
Σ
(mij Tii...., )2
ISh.""
where mi.-. Ex~ X Xiiel and ni,--i, Ex [(ii X[ix]] are corresponding moments. The positive parameters dj. A2. Ag. Xa are chosen to fix the relative importance of the mean, the co-variance, the co-skewness and the co-kurtosis.
3

---

Under review as a conference paper at ICLR 2021
Thus, will be interested in the following optimization task (analogous 6o
duM((aatta v) min
MEPx
UDR: Wasserstein distance (WD) Another important distance between probability measures that has the origins in the transport theory is the Wasserstein distance Villani (2008).
Let (R", I|- II) be a Banach space. Between any two Borel probability measures μu, on R" with |||] do and [|x|||e < the Wasserstein distance is:
/ ||xdT
w(,,v)=
inf xEll(e
y/|dT
where [[(,, v) is set of couplings y and v. The Wasserstein distance defines another version of L.RR problem:
W (Hudata: v) + min DEP
(8)
In the appendix B one can find proofs that in the case of L1 norm |
Σ; lz the task
corresponds to the well-studied robust PCA problem Candès et al. (2011). If, instead of the L1- norm, we the L2 L2-norm, this leads to another well-studied task, which is known the outlier pursuit problem Xu al. (2010).
Sufficient dimension reduction (SDR) Given labeled dataset {(x,, 1 where x; E R" :/ (C a finite set of classes for a classification,, or R for a regression problem), the suffi- cient dimension reduction problem
be informally described a problem of finding vectors W1,' Wa E IR" such that p(y) [ x,, ,"," wfx) Ra p(y(x) (possibly, under some additional as- sumptions on the form of p(y x)).
We formulate the SDR problem as an optimization task:
inf J(f)
(9)
SEF
The object f R -> R is a smooth real-valued function. We assume that f is a candidate for the regression function and (J)) is cost function that values how strongly 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:
J(f) = B f(x; c
N
mm
N
es(x)te)
J(f)
Σ
E E
yi, ex(tte))
N
where H(,.)) y log (1 y) log(1 p) and v > 0 is parameter.
By requiring Fk. we assume that the regression function I satisfies (for & fixed in advance):
f(x)
where Wi,'.., W. E R". Thus, given an input x, an output of f depends on the projection of onto span(w (W.............................. w&). The set span W1, ... Wi) is called the effective subspace. Note that the way we defined the SDR's objective J(f) for the regression and the classification cases is not unique. There definitions that has the same form [ but deal with the conditional distribution p(y(x) x) as an argument, instead of the regression function.
5 REDUCTION OF THE OPTIMIZATION PROBLEM TO ORDINARY FUNCTIONS
The central problem that our paper addresses is how to minimize an objective function over 9j (or ?Pt)? In this section describe approach based on penalty functions and kernels.
Let us assume for simplicity that M is the gaussian kernel, ie. M(x = Gx - y) where
2x]2
(x)
20 Besides the gaussian kernel our theory also captures many other kernels,
V
4

---

Under review as a conference paper at ICLR 2021
1 x
including cases of the Abel Kernel
o" and the Fourier tranform of the Abel Kernel, the
Poisson kernel: Go x0²+ "
For f.g S(R") let us denote:
(f) def [| ((x)))((x)))))) S max (x. y)ll IIlle, < oo
x.y
For general f.g E (1²)) the expression (f is defined if 3ffe 9. E S(R") such that Ty, = f * G., Tg. = g G and (f(l)) 5 A. Then, (f| def M)),0).
A. For example, (6" (M|8))) =
Theorem 3 concludes, from f E Gs. that dim spang {21 f.i.ffffff............... ,2,f....,.ff I k. Using the kernel M. one can build the Gram matrix from the collection of distributions, [fzf fllssse For ff any f E S' let us denote a real part of the Gram matrix (Re (I.f( xy f filsssn by My (if is defined).
Theorem 4. f E Gk. then (faf\M[z,f) is defined and rank M, S k.
Definition 1. Let R be positive semidefinite matrix with eigenvalues 2 2 A2 2 .... dn (with counting multiplicities). Then, the Ky Fan k-anti-norn of A is ||A|| ls = Ei= Ant1-k.
Let R()) ||M Theorem tells us that that for f E GR. R(S)
For ordinary f. the Eckart- Eckart-Young-Miin Young-Mirsky theorem gives us R(f) minn rank ASA |I All Thus,
by penalizing the value of R(f), we enforce M, to be close some matrix of rank For I 9 USR()) + IRt, it is natural to reduce the optimization task over tempered distributions
I(f) + min SECL
(10)
to an optimization task over functions with a penalty term R:
+ M| I() XR() + inf fES(B")
(11)
Details on the conditions, under which this reduction holds, can be found in the appendix D Let us now concentrate on the task [1 and describe the alternating scheme for its solution.
THE ALTERNATING SCHEME
We will concentrate problem 11 It is known Hiai (013) that the Ky Fan anti-norm is concave function, i.c. R(ő) = |M depends on Mr in a concave way. It can be shown that the dependence of R(()) on á is both non-convex and non-concave, i.c. deal with non-convex optimization task.
The kernel M(x, y) R" x R" O(MIffI = Jgm (x..)))x() its range as R(O0 = (0(z){ E H₂}. Let B(H, H2) denote a set of bounded linear oper- ators between Hilbert spaces H1 and H2. For O
induces a linear operator from L2(R) L2(R)): For any operator O between spaces Hi and H2. we denote
H(H,,H2) the rank of is defined as dim R().. Let LJ(R") be the Hilbert space (over IR) of real-valued functions from Le (KR" and L2 (R.) = LF( L (R"). The space (R") is equivalent (R(") treated linear space over R. Below we do not distinguish (01 Ø2) (FR)) and Φ1 + idz E L2RR).. It is easy to see
that any O E B(L (R²), R²) can be given by formula:
[[o], Re (O,, Φ) LR(R),, O₂ (2 " I,
Le. O E B(L;(R"), R") can be identified with a vector of functions O = [01l O2 E L2(R))
and the Hilbert-Schmidt norm on B( [.3 (R"), R²) (i.e. vTrOtO) is:
||0||. =
110:112 (R))
(12)
5

---

Under review as a conference paper at ICLR 2021
Recall that for the kernel M, (MM) is positive and self-adjoint. Since ((M) is also bounded, then the square root vO(M) can be correctly defined Rudin (1991). For any complex-valuet function let us introduce a linear operator S L3 + K" by the following rule:
Sy[ø]; = Re (x; v VO(M) (2()) i.e. (Sy): (x
Theorem 5. If Tr My < ca, then Sy E (L) R2) and Sy = My. Moreover
min R() BEB(L1(R"))EE Jaank ssk
||| - -
and the minimum attained at S S where P==E: u.u! and {u,} ! are unit eigenvectors of M, corresponding to the k largest eigenvalues (counting multiplicities).
Given the new representation R() (EB(T.((")))) R). rank SSE
min
||S} - -11 it is natural to view the
Task 11 as a minimization of T(Ø) A1|Sa 511 over two objects: & and S E (L) B rank
minimize a function over two arguments is optimize alternatingly, i.e. first over o, and then over S rank S < ki and so on. Theorem 5 gives that the minimization over S is equivalent to the truncation of SVD(S,) at the k-th term. This idea, that we dub the alternating scheme, is described in Algorithm
The simplest approach
Algorithm The alternating scheme for I1
Po + 0, S. { 0 for t = 1, .............................. T do
Φ: - arg min (0) + 211s, - Pt-1SA- (minimizing over 6)
Calculate Mo. and find [vi]] s.t. Ms. V, Av, 2 ... An
P Ei=1 v.v? (Truncated SVD(S,,) is P.S..)
Output: Vi,-. Vk
The alternating scheme | allows for reformulation in the dual space. By this we mean that in Scheme I
we substitute Φ for the original Φt. If the primal Scheme
I deals with operators Sg, So.-.>
the dual version deals with vectors of functions V G. can be found in the appendix F
0e Ox
VGo
08,
Details of the dual algorithm
Ux
EXPERIMENTS
The alternating scheme is general optimization method which needs to be specified for every optimization task. We designed numerical specifications of the alternating scheme [ for all 4 opti- mization tasks: 7 g and Jnd and made experiments with all of them. Details of the algorithms, i.e. numerical methods to minimize over & and calculate Ma., can be found in the appendix (G H |
and D). Note that for Wasserstein distance minimization we exploit the alternating scheme the initial form (i.c. []. and for MMD 6 7 and SDR 9 we use the dual version of the scheme.
Behaviour of MMD for small h. We studied the difference in the behaviour of PCA and solution of 6 obtained by the alternating scheme [ (MMD), 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 low-dimensional manifold Mt, which is bent globally, goes through the origin O and has a large curvature O. Because PCA global method and points do not lie on an affine subspace, interpreting principal directions is not straightforward.
and generate points in the following way: points x1, X2, Xx [-10, 10 are sampled uniformly, after calculation of y/ = f(x;) we add some noise: 2; = [x,, y/ + C, G; ~ Nf(0, 0.01I,).. Both PCA and MMD are applied to the dataset (first 3 pictures
We select a smooth function f 1B-11 + R, such that f(0)
Figure Tap. As we see, MMD, unlike PCA, tries to catch ideal alignments of points rather that searching for a global alignment of points (which can be non- existent). This property of MMD makes it promising tool for the calculation of the tangent space

---

Under review as a conference paper at ICLR 2021
(a) Visualization of outputs of the PCA and MMD methods. MMD (green line) tends select subcollection of points that sharply aligns along the main direction, whereas the first principal component (red line) could be a result of averaging over different directions in the data.
(b) Left plot: || P 0.05, 20.0, case 0.05, À 20.0, case 0.05, À
100.0, case I, 0.05, À = 100.0, case II. 8 0. 100.0, case I, 8 0.1, A 100.0, case II.
Right plot: P - P|[r as a function of In az MMD, HM, WD.
to data manifold at given point. Fourth picture shows that when we have equally important directions data such that the first principal direction of PCA is between them (red line), and we set k = I, then MMD (green line) always chooses one of those directions.
Experiments with outlier detection (MMD, HM, Wasserstein distance). Following the experi- ment setup of Xu et al. (2010), we choose parameters N = = 400, 5 = 0.05(0.1), = 10 and generate random matrices A E (1 B E Ba*k whose entries are iid /O, 1). Then, to the columns of the matrix BAT B=xN((16) (whose rank is < k) are concatenated with the columns of the matrix E RaxN6; X = concat(BA" C) E RaxN The entries C are either üid as (0, (case I) or N6 copies of the same vector whose entries are iid as NN0 (case II). Let X = [x,, xw]. i.c. columns of X are the data points. Thus, NV - 8) columns of BAT lie in k-dimensional subspace of R and N6 columns of C are outliers, and solutions of tasks or for this dataset are expected to be supported in a column space of BAT.
After every iteration (step of the alternating scheme calculate the Frobenius distance between the projection operator P. of 1 and the projection operator P to the column space of BAT i.e. ||, PP.-P[[r. For the task 8 the dependence of P| P|]r for different values of parameters 8 and À is shown in Figure iterations are enough approach the optimal subspace.
lb For tasks 6 7 the behaviour of the alternating scheme is similar,
Besides the speed of convergence we were also interested in how ||P* Pl where P lim; Pt is the final projection operator (e.g. P2o in practice), depends on the parameter of the kernel M = G.. It is natural to expect the quality of the solution p to degrade as + too (this corresponds to M(x.y) + 0), and, less trivially, as a + (this corresponds to M(x.y) + 80 (x-x))
Experiments with the sufficient dimension reduction. We made experiments on the standard datasets, Heart Breast Cancer lonosphere, Diabetes, Boston house prices and Wine quality First we applied Sliced Inverse Regression algorithm (SIR) [ (191) 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- dimee-iimesiioal representations of input points. In the last step we applied 10 nearest neighbours algorithm (KNN) to predict outputs (based on reduced inputs) 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 Fukumizu al. (2004) and the alternating scheme T adapted for SDR.
7

---

Under review as a conference paper at ICLR 2021
We experimented with the dual version of algorithm setting (after the data was standardized) the kernel's 's parameter o = 0. and À = 10.0. Details of its numerical implementation can be found in the appendixpi In the table one can see the obtained test set accuracy the classification tasks and R² on the regression tasks. As we see from the table [ after reducing the dimension of an input to k 2, we still able obtain good accuracy of prediction test set.
Method
PCA
STR
KDR
AS
Dataset
Dimension k
2
3
2
3
2
3
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
lonosphere (acc)
80.29
86.57
89.14
89.43
83.43
86.29
88.29
90.57
Diabetes (R²)
25.34
28.72
43.47
43.61
41,82
44.30
43.07
44.48
Boston (R²)
56.42
67.12
76.03
74.29
77.88
79.97
73.21
77.88
Wine (R3)
93.91
94.12
98.68
99.24
98.30
96.02
97.10
96.93
Table I: The cross-validated accuracies/R² KNN on or 3-dimensional input representations.
The code is available on github to facilitate the reproducibility of our results.
RELATED WORK
We present optimization framework in which the search space is g%-- or Py- Another unifying framework for LDR tasks is suggested by Cunningham & Ghahramani (2015) in which the basic search space is the Stiefel manifold S(n, k). The main disadvantage of using Jr instead of the Stiefel manifold, is that its infinite number of dimensions requires special procedure to turn an op- timization into a mmite-dimensional task. Both an optimization over G/ and over S(n, k) is typically hard: for final point, at best one can guarantee that it is a local extremum. Promising aspects of 9% are: a) 9% allows formulate naturally a new class of objectives on it, b) local extrema on 9 substantially differ from local extrema on S(n, k). because local search over 9% uses more degrees of freedom.
Using Ky-Fan k-antinorm a regularizer for the matrix completion problem has been suggested by Hu et al. (2013) and further developed in Oh et al. (016)) Liu et al. (201)); Hong et al. (2016). Unlike this chain of works, formulate an
task. Also, regularizer R(f)
||M is a sum of smallest n k squared singular values of the operator S where Sy depends on
linearly. The idea of alternating two basic stages, the convex optimization and SVD, is ubiquitous in low-rank optimization, see e.g. Mazumder et al. (2010); Hastie et al. (2015).
Zhu & Zeng (2006) applied the Fourier transform for estimating the effective subspace in SDR, implicitely using an analog of Theorem
CONCLUSIONS
We develope 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, a, and the penalty parameter X. We believe that the MMD/HM/WD methods for UDR could be used as alternative to PCA in study fields in which data demonstrate "heavy-tailed" and "non-gaussian" behaviour, such as financial applications. Also, our formulation of SDR is free from any assumptions on the distribution of input-output pairs, which makes it an alternative to other methods the efficient subspace estimation. More detailed report on these topies is a subject of future research.
Since the role of the parameter a is similar to that of the bandwidth in the kernel density estimation, we Silverman's 's rule of thumb to set N-1/6644)
8

---

Under review as a conference paper at ICLR 2021
REFERENCES
Martin Arjovsky, Soumith Chintala, and Léon Bottou. Wasserstein generative adversarial networks. In Doina Precup and Yee Whye Tch (eds.), Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pp. 214-223, International Convention Centre, Sydney, Australia, 06-11 Aug 2017. PMLR. URL http://
proceedings. mlm
html
A. R. Barron. Universal approximation bounds for superpositions of sigmoidal function. IEEE Transactions on Information Theory, 39(3):930-94 May 1993. ISSN 0018-9448 doi: 10.1 109/ 18.256500.
S. Bochner. Vorlesungen über Fouriersche Integrale: von S. Bochner. Mathematik und ihre Anwen- dungen in Monographien und Lehrbüchern. Akad. Verl.-Ges., 1932.
Emmanuel J. Candès, Xiaodong Li, Yi Ma, and John Wright. Robust principal component analysis?
J. ACM, 58(3), June 2011. ISSN 0004-5411. doi: 10. 00.11445/9773992119003953 URL https:
//doi. 11w...111199999 htt.:///1709999999909999 1970395
John P. Cunningham and Zoubin Ghahramani. Linear dimensionality reduction: Survey, insights, and generalizations. Journal of Machine Learning Research, 16(89):2859-2900, 2015, URL
http://jmlr.ox/ www.ooogggppp......o //v//vvv5. http:///////uunngnaaaaaaaaaaa
html
F.G. Friedlander and M.S. Joshi. Introduction to the Theory of Distributions. Cambridge University Press, 1998. ISBN 9780521649711.
Kenji Fukumizu, Francis R. Bach, and Michael Jordan. Dimensionality reduction for supervised learning with reproducing kernel hilbert spaces. Mach. Learn. Res., 5:73-99, December 2004. ISSN 1532-4435.
A. Gretton, K. Borgwardt, M. Rasch, B. Schölkopf, and A. Smola. A kernel two-sample test. Journal of Machine Learning Research, 13:723-773, March 2012.
Ishaan Gulrajani, Faruk Ahmed, Martin Arjovsky, Vincent Dumoulin, and Aaron Courville. Im- proved training of wasserstein gans. In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (eds.), Advances in Neural Information Processing Systems 30, pp. 5767-5777. Curran Associates, Inc., 2017. URL http: /ttp:///ppppr..... nips cc/paper/ 7159- 1mp rwverddddddooooooooooooooo
of www.erstei......... gans.pa........
Trevor Hastie, Rahul Mazumder, Jason D. Lee, and Reza Zadeh. Matrix completion and low-rank svd via fast alternating least squares. Journal of Machine Learning Research, 16(104):3377-3402 2015. http: http:///ww......./////////////aa
Fumio Hiai. Concavity of certain matrix trace and norm functions. Linear Algebra and its
Applications. 439(5):1568 1589 2013. ISSN 0024-3795. doi: http:://doi.org//..000///
laa 22200013..000 URL http:/ www. www.::seeeecceicc....c...cm com/ http:///eetttccccc/ccc///////rre cctmmm/ccceecceaatile/ppp..
S0024379513002882
Bin Hong, Long Wei, Yao Hu, Deng Cai, and Xiaofei He. Online robust principal component analysis via truncated nuclear norm regularization. Neurocomputing, 175:216 222, 2016. ISSN 0925-2312. doi: https////laiiir/////oo.o////ooop
URL http: //tp:w//www
sciencedirect. com/ wwtp:///ww...neee......////////a
T. Hsing and R. Eubank. Theoretical Foundations of Functional Data Analysis, with an Intro- duction to Linear Operators. Wiley Series in Probability and Statistics. Wiley, 2015. ISBN
9780470016916.
Y. Hu, D. Zhang. J. Ye, X. Li, and X. He. Fast and accurate matrix completion via truncated nuclear norm regularization. IEEE Transactions on Pattern Analysis and Machine Intelligence, 35(9);
2117-2130, 2013.
Ker- Ker-Chau Li. Sliced inverse regression for dimension reduction. Journal of the American Statisti- cal Association, 86(414):316-32 1991. ISSN 01621459. URL http:// / / www. . jstor. org/ sttable/229055

---

Under review as a conference paper at ICLR 2021
Q. Liu, Z. Lai, Z. Zhou, F. Kuang, and Z. Jin. A truncated nuclear norm regularization method based on weighted residual error for matrix completion. IEEE Transactions on Image Processing, 25 (1):316-330, 2016.
Rahul Mazumder, Trevor Hastie, and Robert Tibshirani. Spectral regularization algorithms for learn- ing large incomplete matrices. Journal of Machine Learning Research, 11080))2287-2322, 2010.
URL http:h/tmmmmi. http://mmmmmrr.ro/ppaeeeoee wwtppppppp.pppoppopotooppppw /vlas/ss/lii/// htal///mmmdderrreeeeeeee html
Krikamol Muandet, Kenji Fukumizu, Bharath Sriperumbudur, and Bernhard Schölkopf. Kernel mean embedding of distributions: A review and beyond. Foundations and Trends@ in Machine
Learning, 10(1-2):1-141, 2017. ISSN 1935-8237.
doi: 10. 101151122000000600 URL http:
77dx. www.dooi..... org 0. 1561 22200000600
T. Oh, Y. Tai, J. Bazin, H. Kim, and I. S. Kweon. Partial sum minimization of singular values in robust pea: Algorithm and applications. IEEE Transactions on Pattern Analysis and Machine Intelligence. 38(4):744-758, 2016.
Henning Petzka, Asja Fischer, and Denis Lukovnikov. On the regularization of wasserstein GANs. In International Conference on Learning Representations, 2018. URL https:// openreview wwwpeeeiwwwwwweeeeeeee www.eet///oooe fereet/foormmeee Ld-BLLYEMC
Walter Rudin. Functional analysis. International Series in Pure and Applied Mathematics. McGraw- Hill Inc., New York, second edition, 1991. ISBN 0-07-054236-8.
Laurent Schwartz. Théorie des distributions et transformation de fourier. Analyse Harmonique, Colloques Internationaux du C.N.R.S., 15:1-8, 1949.
S. Soboleff. Méthode nouvelle à resoudre le problème de Cauchy pour les équations linéaires hy- perboliques normales. Rec. Math. Moscou, n. Ser,, 1:39-71, 1936.
C. Villani. Optimal Transport: Old and New. Grundlehren der mathematischen Wissenschaften. Springer Berlin Heidelberg, 2008. ISBN 9783540710509, URL htttht:s///oook google.
kz/books 2 5tkC
Meihong Wang, Fei Sha, and Michael Jordan. Unsupervised kernel dimension re- duction. In J. D. Lafferty, C. K. I. Williams, J. Shawe-Taylor, R. S. Zemel, and A. Culotta (eds.), Advances in Neural Information Processing Systems 23, pp. 2379- 2387. Curran Associates, Inc., 2010. URL http::/papeer: nips. cips.cc/paper/ wtrp-uunsseriiseddeeeeeee -iw.ennl-ddimssiieeeee reduction.
Xiang Wei, Zixia Liu, Liqiang Wang, and Boqing Gong. Improving the improved training of wasserstein GANs. In International Conference on Learning Representations, 2018. URL
https://openreviie net meet/FoouumlddSS35SCC
Huan Xu, Constantine Caramanis, and Sujay Sanghavi. Robust pea via outlier pursuit. In J. D. Laf- ferty, C. K. I. Williams, J. Shawe-Taylor, R. S. Zemel, and A. Culotta (eds.), Advances Neural Information Processing Systems 23, pp. 2504. Curran Associates, Inc., 2010. URL http: 7/ //papers.... nips wwtppcpppppeeeeteeccceccccc..e http:///oo...oooooooooooaaaaaaa- -
pdf
Yu Zhu and Peng Zeng. Fourier methods for estimating the central subspace and the central mean subspace in regression. Journal of the American Statistical Association, 163844) 2006. doi: 10. I
A PROOFS FOR SECTION 3
A. PROOF OF THBOREM 1 GIVEN FOR COMPLETENESS
Proof. The inclusion 0% S Th follows from well well- fact that S(R*) is dense in S' ('(RR)). Le. for any E ('RI"') one can always find a sequence {ft} C S(RR such that Ty, - f. Therefore, for any (f 8 50= E G there is sequence ((T,, @ kn A) S Go such that (T,, @ 8" k) 4 (ʃ ).. Thus, 3% C Gz.
10

---

Under review as a conference paper at ICLR 2021
Since 9 S 9% prove 9 9% it is enough to show that g is sequentially closed. We need a simple fact from a theory of distributions.
Lemma 1. IfT, - T and o; + 6. then (Ti, di) + (T, Φ).
Proof of Lemma. Schwartz space S(R") is Fréchet space, therefore the Banach-Steinhaus theo- rem applies S' ). Since T; - T, have sup, [T;: d) theorem, applied to set [Ti] we obtain for any >> 0, there is neighbour-
for any & E S(R"). From the
Banach-Steinhhe
hood U of S(R") such that [(T,, <)] whenever Φ U. Thus, KT,, 0z Φ) < for large
enough i. From that we conclude that (T), 0i) + o).
For any T ('(R") and & ES(R"=)), let us define T E SR( as (TV o) = (T
Suppose that {fil ( { are such that (fr @ 8"~*)u; +' f. We need to prove that I E Gz- Since set of orthogonal matrices is compact, then one can always find subsequence (Un,) } such that Un, - U. Since (fn, @ 8m- Ju,, - and Φ(Un,x) + @(Ux) (for any fixed
S(R")), using lemma | we obtain:
(Fn, @ 8" ②) = {(Sn, @ 20 8mt Un, o(Un,))) x)) - (I,((Ux)) = (Foz,((x))
Thus, we have:
Sn, 2 5m~ -* fUT
From the last we see that Sn, -" f0r where & is such that (()) = 1. Therefore, Sor für 608
and f = (frr @ 8" 8"-k)o E Ik.
A.2 PROOF OF THEBREM 2
Proof. Let us prove first that if g = -Ty 2 5n then
Flg]=
where r(x) f(x11)), E R". For that we have to prove that (F(g). = (T Φ) for any Φ ((R).. Indeed,
(F(g).0) = Flo) (T, 80 øly)e oly)e-x"xds xdy)))
(TB, / dly)e M ix JRIK [(xxx) i dddxxxx dydx1=
/ f(yn))olydd = (Tr,0)
JRn
Let us calculate the image of Is under the Fourier transform. easy to see that for any g S' S'(RR), S(R") and orthogonal U E Raxx we have:
(Flou). d(x)) (ou,
(g. F(o((())
(g,F( x))) (F(g), pTUT x)) ((Flg)u,, o(x))
Therefore, Flgol = (Flg)v. Thus, if g Ty 8 8- then
(Flgul) (Tr) = T
where rx) = r(Ux) = f(Uxx) and Up kx is a matrix consisting of first rows of U. Thus, T. E Fk.
Let us show that by varying f S(R*) and U in the expression f(Uxx) we can obtain any function from Ft. For this it enough to show that F& is equivalent to the following set of functions:
º= = (o(U,x)le E S(RR), ), U, E R _U = [s)
The fact Q Fg is obvious. Let now prove that Q 2 (s(Px)lg E S(R*), P Rkxn rank P = k) Fr. Indeed, if fxx) o( Px), then f(x) g (Uxx) where Uk
(ppT)-1/22 and
d'(y)= = ( By construction, U.U. In and g' E S(RR). Thus, Q = Fr.
Therefore, Flgg] = Fr, and from the bijectivity of the Fourier transform we obtain JF F Gk.
11

---

Under review as a conference paper at ICLR 2021
A.3 PROOF OF THEOREM 3
Proof of Theorem (+)) Let us prove that from T = ( 2 )u, E S'(R)),U" = In
follows that dim spang [I1T, z2I.... .IIT. S k.
It is casy to see that lf @ 8-=k]= = i > k. IU [u, unl", then for > k we have 0 (z. f 8 "I)u u x(f @ 8n~*)o = uf XT.
Thus, we have n k: orthogonal vectors, uk+1,-- Un. such that [a] ... Jn Ju; 0. Using standard linear algebra we obtain there are at most k distributions 2,, T, ... 3ig T. S that form a basis of spank {2,T])
To prove the second part of theorem need the following lemma. Lemma 2. If' E SRR))) is such that y:T 0 for any i k, then T E lk.
Proof of lemma. Recall from functional analysis, for f S'(R"), the tempered distribution 0/ is Dé
dz. [J, dz, Once the Fourier transform is applied, our lemma's dual version is equivalent to the following formulation: if öf = 0, > k, then f E FR Let us prove it in this formulation.
aL defined by the condition (āx, Φ) -
Sw:
A of infinitely differentiable functions with a compact support is denoted by C (RR). Suppose S(R") and p E C (R) chosen in such way that J p(yx)dy: 1, suppy p C [A, B]. Let us define:
(1x-i18))d S p(u)dd: I
r(x)
((x-in v:)dy:
It is easy to see that for any a E Nm- over æ, is present):
EN 1M-1, , N we have (at least one derivative
08.148T7
88,8 (0,8(((x) p(xx)
1 (xx 1.k-,,,))duu y1)dul
x
x
Ox"
.xx
1+80
0x~,0t ,xx""
0 6(x)
29' p(x:)
0
J px Mi dy:
0x 3 0x,
dx
θx=;
The terms x" x 00.00 6(x) and n 0² p(x.) vx_,dn
are bounded by the definition of S(R"), C. (R). The
oz."
2 (1xx)))11) dy: is a consequence of the inequality (which holds because boundedness of - xa Ox"
S(R")): xx" soxx 19t1) I
1t47 1
dx".
Analogously (not single derivative over z; is present):
0
I 00 vs) dy: - ag' 1 p(yi)dy; L xa
2
100(x Ui)
naa
xa,
dyfi
Ox
dx
dx
22 plyi)dyi)
00 00 days n2f P plyi)dy:
10
0 (xx dys
ad (1 1
(1-
xa;
0x
Ox
8 p(x i. ox_,,
The second term is 0 when 2; S A. It also bounded when v; A because xx"
C² and:
(1+
0 (xx(x.))))) dy:
S | f
O'
x²
(1 + H 1 dy
Ox"
The latter is bounded because lim [2/]'" s (1
dy, = 0.
12

---

Under review as a conference paper at ICLR 2021
The first term is when 2; 2 B and it bounded for 2; < B:
Fe
2 (x
C²
dyj
S |0il"
dy
0x_11
+ y )2))a'+
The latter is also bounded, since lim 7-20
|2/| S (1+1,) ) dy: = 0.
Thus, xa r(x) is bounded and S(R"). Therefore
85 Om.
0 implies:
Ox
U, Or: - 0 [| = s(p(x;)
(xx-i, yi ddy
Since this sequence of arguments works for any i > k. we can apply them sequentially to initial Φ S(R") w.r.t. 1k+1,,..., Vn and obtain for any Pk+ Pr. E (R) such that S Pi Pi(yiddy:
f(ø] flps+1 ( (I() [ Rm
o(x1 Xd
Moreover, since C (R) is dense in S(R), assume that Pk+1, 1.. Pr. S(R). For the inverse Fourier transform 7 = [-] the latter condition becomes equivalent to:
(T. (T, ( ... p'a (In ((t)))(X::k)) Ox+1s)))
for any Phrth E S(R) such that (0) = Leet us define pf (x:) = e-
It is easy to check
that g & 8a-k where g E (R (9. vb) (T, |8 mmdl (x1x)))) for E S(R*). Le.
-
T E 9 and lemma is proved.
Proof of Theorem 3 (=). If dim spang (x. x2T, T. .I,,,,,,, ,InT} k. then
dim(v RR" [2) In T' 0} 2
Thus, there exist least n k orthonormal vectors VA+1, Vn. such that [ryT. ... 2,, 0. Therefore, [x T, ... T = (v: x)'x""==00
Let us complete Vk+1.-.. Vn to form an orthonormal basis of R": Vi Vn. Let define matrix V [V1. vn.. It is easy to see that:
(vv]x)))
(v.Vx)TV x.IV
Since for i > &o we have (vF x)T 0, then x,TV 0. Using lemma obtain Tv G Therefore, (Tv)vr = T E Gt- Theorem proved.
B STRUCTURE OF WD
Recall that (R", ||- []) is a Banach space. Now, let us consider optimization problem: for given X TRx. solve
||X Lll
min rank(L)gk
(13)
where ||
is extended RaxN by !||s),... sn]ll
Σx 1ss.ll
The following simple theorem shows that the two tasks are connected, so that the solution of one directly leads to the solution of another.
Theorem 6. Given data points (xI,".., let X = (x,,"" [x]] RaXN Then,
1
min W (Mdatas v) N
min
X--YII
PEP
YEREXN (a)
Moreover, minver, W (Pdatas v) is attained on v", where v²t a uniform distribution over {y:) and "y,.... .yn] arg miny eRxx "enkkN..nn(""" (onk((Y)sk X|X VII.
13

---

Under review as a conference paper at ICLR 2021
Proof. Let prove first that inf uEP, W (Pdatas p1) S x - Y+1 where
Y [y,, yN E arg min VERaxN yeRoxx, (rank(Y)
||- - YII
Let ı be uniform distribution over y:)) and μ be uniform distribution over {ys Since π E Muaatas j²),, we obtain W (Pdata: v=)) S N > ||x y:ll = & Y The support of j²t is k-dimensional, because vank(Y*) S k. Thus, have , Pk and infueP, W (Pdata: p.) I W (Hdatas:P p*) S & - YI Now, if we prove the inverse inequal- ity, i.e. infNEPe W(Paata (4) 2 MIX 1l, this will imply that inf, W (Maata: P)) + -| and therefore, infuEP, W W(Hdata ()) W ((ddatas ((") This will in the end give us u E aarg inginf W (Hdaaa: 11).
Let {u be such that 1l Pp and W (Pdatas (4) inf W (Hddtta: H1) O. Let L, denote k-dimensional support of #4 and Pt is a projection operator onto Lt-
Let j4 be a uniform distribution over (Pix1,".. Pxxn) i.e. easy see that W (nt (dasa) I W (Pt: Pdata). because je and fr share the same k-dimensional support Lr but the "transportation of a mass" concentrated in point x; of the empirical distribution Hemp can be most optimally done by just moving we have pEP, W (Haatas (4) I W (Pdata: (7) I W (ddata: Mc), and therefore, W(Hdata: m4) inlREp, W (Haataa p)) - O.
E A]. is
Pix; (i.e. to the closest point on Lj). Thus,
Since set of projection operators is compact, one can always extract a subsequence {Pr, 1.00 such that Pt. - P. is easy to see that Jr, -> μ (i.e. , - 0) where / is uniform distribution over {Px1,.... PxN]. For that distribution have W W(ddtta: p²*)
lim W (Udatas (i,)) unf W W(ddttas ut). Thus, the infinum is attained on n²"
is easy to see that W (Adata: p**) = W("", Mdata) = & |X - PX||. Since rank (PX) I & we obtain W (Adatas 1=*))= 2 N milly san (aanl ||X - Y|. This completes the proof.
Note that the case of L1 norm |x||
= >, |as| in the task 13 corresponds to the well-studied robust PCA problem [andès et al. (001T) If, instead of the -nonom, we use the Ly-norm, this leads to another task:
|| - L1,,2 +
min
(14)
rank(L)gk
where ||5||1.2 Es VE Bt. which known as the outlier pursuit problem Xu al. (2010)
C PROPER KERNELS AND PROOF OF THEOREM 4
C.I PROPER KERNELS
In the main part of the paper we assume M to be a gaussian kernel, though the theory can be applied to more general case of the called proper kernels.
Recall that for any operator O between spaces H. and H2 we denote its range as R[O] = (0())) e H.J. For Q2 C S((R",, Q denotes the sequential closure of Ω with respect to natural topology of S(R"). A set of continuous functions in 1n is denoted by C(R*). A set of infinitely differentiable functions with compact support in I" is denoted as C (RR)
Definition 2. The function M Mx, y) R x R. -> C is called the proper kernel if and only if
O(M)[[]= Ig.. (x. fy) is a linear operator from L2(")) to L2 (RR),
= ((ly,))MM(x,,
(f:((M))1)))2(77) > 0, Vf E L2(R)), f 00
maxxy (M(x)) <0.
R[O(M)) ((R))) S(R)).
14

---

Under review as a conference paper at ICLR 2021
The gaussian kernel (x,,) = (" y). which is of special interest from application-oriented oriented perspective, is captured by the following lemma:
Lemma 3. If S. C(R²) are bounded, Vx {(x) > 0, then Mxx = Ç(x y) is a proper kernel.
Proof. Verification of the first four conditions is easy, we only check the fifth condition. Let us denote linear operators Cx[f] = {+ ff and O,[f[(x) = fxx) Then we have FIC₂ (2 = 0g[42( 2 Ce (I Therefore, R[O(M)] = Cc(L2( 2 F "CC ]]. Since dense in (IR"), then F [C2 (R²") also has this property. L.e. R[O(M) (RR)) S(RR")
Besides the gaussian kernel the lemma also captures a case of the Abel Kernel ((x) = e 1x is well-known that the Fourier tranform of the Abel Kernel is the Poisson kernel: {(x) (1+1xP)) = (which is also proper).
C.2 PROOF OF THEOREM 4
We will prove a more general statement:
Theorem 7. Let yx. be a proper kernel and, additionally, a Lipschitr,function If E Gv, then (m f) is defined and rank M k.
Proof. Let us first show that (f]M)g) is defined for all f.g E Gk. Note that for any f = (Ta &
8n-k)u Ju E 96 we have
Tj. (Ta 2 8" 8"= G. ((Ta G.)) & +t
Let us denote a, * G and b. = b G.: It is easy to see that
fe = (a. (x1:k (x (x)+1m)) E S(R")
From a well-known property of the Weierstrass transform we have:
||Sellr, ||ac|lr, || Malle
From this we obtain for any f = (Ta @ )u, = (Th @ 8"~*)v E Gsi
[[f[|M max [M (x, y)l IIfelle max M M y) Melle,11611.. <0 x.y
x.y
Thus, L.elM/9s is defined and:
(fll /
(x(x1x)) (X)+)))))()))))) V yb)) (y
R XR"
/ Me F1:: (V1nk ddx
I2 xIl
where
Me Mx Vi.k)
G (X4+m)))) V. (y
XR
Let UR, VA E Rxn be matrices that comprise the first & rows of U, V correspondingly and k zero rows below. Also, let L denote Lipschitz constant for M such that [M y) M(x', x') S L(Xx x + y y')). For the function M. (X y V. we have:
MM V
M( x, V y
G - (Xe+1m)) (M <M xX -M0M-MMUXXX. Vi2 v)) -h-* (Y Im )dx Ion dy k torl
2
15

---

Under review as a conference paper at ICLR 2021
S Ll
C. (x (I(U (1))x/ 1(V yl) G. (y=="(yxtmmm
R2n-2k
L
Gm- (xx++mm) (xx+1mml ( 1m
22n-26
2T.
Xl+1im C ( (XI-+1mmx) 2Te
Rn-k
Rn-a
Thus, there exists bounded (x = M(U x, VfT) such that:
Me (x VI.k 3 M (x Yiik) in Loo (I22k)
Further we assume that e >>>> 0 is small enough, so that M, (x.x, SC= CC-2mam 2mxx IM Now we have:
/ 22 xB
M. Yra)) - a" /
R2 XBI
/ Yix (d)
[ R XB=
a" (M(NA
( S
R$XR X
C'laelle, - bll2, + 1116112,110 - aLLL, + (( IIMe - Mille
is well-known (e.g. see Theorem 2.25) that |la. alle,. ||be bll2, 0, l|aelle, S ||all,, and ||M Milex + O. Thus, lim -o (felM (g)) exists and (f] is defined.
Let now prove that rank My S k. Since E Gu. then f (T, 8m-*)er where is an orthogonal matrix and U [W1, wn). It is easy to see that:
(e.f)Ma,,f) ((=:forr x, U. (=)
(w."x s0 8" 80~ X,'"""UU"" U²
('y)|w/] y)[ w 60 8m
Let us now denote = [ui,... u.] E kx submatrix of U in which only first k rows of U are present. Then, the latter integral equal to:
M(VT X V VT (d uf Bu,
where
[ (< MM M(V X1:k, VJ (1ysk)) B
the Gram matrix of the collection (x.g(x1:))1 C S(R*).
Obviously, rank M, rank [Re ui Bu, Buylissn rank V (Re B)V S rank V = k.
D GENERAL THEORY OF THE REDUCTION FOR SECTION 5
For sequence ( C S'(R"), Lim Ss denotes set of points E S'(RR",, such that there
exists growing sequence (st) C N and tim Ss, f-
16

---

Under review as a conference paper at ICLR 2021
D.1 REGULAR SOLUTIONS AND REDUCTION THEOREMS
For : 9% S + Rt, it is natural to reduce the optimization task 10 to an optimization task over ordinary functions with penalty term [T] To have equivalence between T0 and [[ we need to assume that I's behaviour when approaching f E 9k from a set ((R)) is continuous, L.e. for any sequence {ft] S S(R") such that Ty. + f E Gj+ we have lim ( I(f).
[
Let us introduce the notion of a regular solution both for 10 and 11
Let
B₁ U if GR Tr(M)) CC C20
Definition 3. Any f Arg min I(f) is called a regular solution of 10 fec
of distributions from Gh that can be approached through quences (fi S Gs. for which Tr(M,, ) does not blow up. Obviously, Gs C B, S 9-- In ap- plications, regular solutions include all Arg min T(J) if we choose the kernel M correctly. This FEG
In other words, By formalizes
regularity is important for reduction to the penalty form 11 because when approaching non- regular solution we are unable to guarantee bounded behaviour of M (and R(f)).
Definition 4. sequence {fi]i C S(R") is said to solve 17 if:
I(fi) + X,R(f)) - inf I(f) + X,(()) +
(15)
fES(
where G, - +0 and >; +00, - toc. If. additionally, {)) bounded, then {f:} is said
to solve 77 regularly.
Let us define
rsol (I(),
U
Lim Ty:
(AH r. solves (11)
Theorem 8. M is proper kernel, then rsol (I(f), R(f)) S Arg min I(f).
SE9:
Theorem 9. If M is a proper kernel and rsol ([ff),((f)) * b. then Arg min I(f) B. S rsol ((ff),(f)))
fE%
Theorem 10 (Reduction theorem). If M is a proper kernel, Arg min I()) rsol (((f),(()) # 9, then rsol (I(f), R(f)) Arg min I(f). fE9
C
B₁ and
feo
Suppose that we now solve sequence of problems and 9 the following are potential scenarios:
11 and find {f According to Theorems
(1) Tr(M,,) blows up and the convergence is not guaranteed. This situation can be avoided by controlling Tr(M,) in an optimization process. practice when f has a parameterized form, this can be done by bounding parameters.
Tr(M,). does not blow up, we still have two subcases:
(2.1) Lim Ty. # 0. This implies positive outcome to approach 11 to the optimization problem, Problem 10
(2.2) Lim T,, 0. This exotic situation can happen only if sequence Ty, leaves any sequentially compact subset of S' (RR). Bounding parameters also tackles this case.
D.2 PROOFS THEOREM AND 9
For any f = (Ti @ 8"
Ju E 9b and a > 0, let us define fo as:
Tj. = (Ti 3 6t 8--h)o = G = (Ti, 8 T -
fa (1x
17

---

Under review as a conference paper at ICLR 2021
1 1 G
We have Ts, -+" (Ti by as
Lemma For any E Gb. to | 2 for any (i.j) & [L..,,]'. and
(aPae M[ oc, for any (i,5) {1, (1,... k)²
Proof. W.l.o.g. we can assume that f = Ti @ 8m-k e S(R*). If k,j S k we have:
(I.folMazsfa) (2T02)n //
Isttin!
Tillee
2n
(x(x1x) Mx,
20
La (
k+m
22T02"" 20op""=- -
Lo (x11)) P(x)dx 20
where P(x) = Sun V2ng=
Uy (x,
Lo (yI:)ddy. Using the Hólder inequality we
obtain:
Mktin 2x
RiE
La (x1x (xxx))llaa
/2TT22
1
24C
20
(R) I L Plle_
2ma²
Since |M y)l y for some 7, we have:
P(()) sill
2TTO2 n Uje
20
1
Vt' 1im 20
711 V2ma2" n-k
|1,, (R- M (ya))le,, (Vlusla(y
Thus,
18 1 Ree
Ma S
20
V2T02"
Using [oll -11 of [| (1ll (R+) and proceed:
+0 0. l
(B*))
|l (I(A))
+0 0, we see the boundedness
s C'|l
- I&e
2u
(R=,(R
v V
18 20
It is easy to see that V2na2" n- Tipe (rafalM IM (Mlll (ryJo) Similarly, The entries of the main & x & minor Mo are bounded, because:
(BI
O as
0, therefore
+ O.
-ifii k. prove that (rafo MRrsSo)
Tr Ms (2m02)n-K J/e
Mar1 20
Iy 1n
x-
lo (x1x))
20
l (y
(2 /1e xR
18 1+ vy Iin
(IxIa [x
20a
I (x1x)) (Y S
(iisk - Viskla (V + 1ll (n k)o²
18zlo 11 7116,112, (n k)o²
Again, using ||4, [|,,l,,,(R))) - |1l1,, t o. 1155bolle, - " 0, we obtain the boundedness of RHS.
18

---

Under review as a conference paper at ICLR 2021
Corollary 1. For any Gg. limo-o R(fo) 0.
Proof. W.I.o.g. we can assume that f = Ti 7
"o E ((R)). By lemma, all entries of Mr except
those the main k minor approach 0 as - 0. This means that
lim Q(f)== = a-+0
where Q(fo) - - (rifo (iffo.) Let Vi. Vn be unit eigenvectors of Mr correspond- ing to the eigenvalues A 2 ... 2 An. p= Ei= e then
R(fo)= = Σ A=
Xips S - Q(fa)) Σ
min AE P. 1=1
k=++1
1=1
Since R(fo)= S Q(fo), we obtain li R(fo) = -.
D.2.1 PROOF THEOREM 8
Proof. Suppose that sequence (fi) S S((R)) regularly solves (7) and T Lim fi. W.l.o.g.
we can assume that Tç +* T and Tr(Ms;) is bounded and I(f:) + X,R(s, S 1-100
inf I()+ +
SES(Ra)
X.R(f) Edi 6, 0. Below we use continuity of I and corollary !
inf I(f) + +R(() & inf inf I(fo) A,R(fo) ( S fES(Rn)
inf lim I(f=) + A,f(f) inf I(f) fEin
fES a20
feEx a-+0
from which we conclude that R( inf I(f) and, therefore, R(fi) '230 0.
SEGx
For each let us define P1 as the projection operator to a subspace spanned by first principal com- ponents of the matrix /M i.c.
P
where vi...., va. are orthonormal eigenvectors that correspond to k largest eigenvalues of v MM, From the Eekart Yoautt-Youn---rssse theorem we see that R(fi) IIVM, P Since set of all projection operators {P E Ruxx p2 P, PT = P} is a compact subset of R" one always find projection operator P P= E v.v." and growing subsequence (4,) such that I|Pr Plle + 0 as + oc. Thus, for the subsequence {fi,} } we have:
IV M PV V Mr 11. 11 Mr. Pt, M₁, Pt, AL PP PM. Il
IV Mn. PI. Mr. lE + +II., Pilelly Mr.
(R(f.)))))))))))o) + |P,, PlIE Tr(M) IFF 1I
PIlF
and using the boundedness of TT(M,,) we obtain I PyMM 1 +
0, let complete Vi,..., Va an orthonormal basis Vi, Vn and make the change of variables y/s = v x. Let us denote V = [Vi,..., v.. and let V = [Wi.... w..]. Then, after that change of variables any function fxx corresponds to f'(y) f(Vy) and the kernel M corresponds to (y. = M(Vy, Vy'). If we apply that change variables the integral expression of Me f), will obtain:
Since
(I = (wTyl' [W] = wT Kuv f' f') Wj
Re (a.f) = w [ee (y M' f'yl W;
Le. M, VM,.V. Mr. - M,V. Note that P VITV where Iht is diagonal matrix whose main & x 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 yt /M,, V =
V /VT M V, we obtain:
Mr, - PPMM, I= / V MA,V Vll==
V
19

---

Under review as a conference paper at ICLR 2021
IVVTMMV V V VVV V Mn Vllf = MYH - 1h M lle
Thus, the property I - PVM,, F + 0 implies that:
Re (wi (T,MM' 1,, - 0, IF IF k
Moreover, for i = j we have Re (w. |M 1y, fi,)) (v. MM,,,, It is easy to see that after the change of variables we still have fi, + Tv. Since fi E S(R"), we have yrfi, E S(R") and, therefore, y:fi, E L2(I)). Let us treat now M' as an operator O(M') (2( (Rx",, (MM) = Je- M'(x, y) Let us take any function E L2(R²) such that v = (0) E ((R"). Since O(M") a strictly positive self-adjoint operator, by the Cauchy-
Schwarz inequality, we obtain:
O(M' & V (u.JI, (o. ((M))
Therefore, for any vs E RIO(M') n S(R") and i > k we have lim, (V.fi, 1) = lim (fi, m.b) = O. Since S. +' Tv we obtain (Tv. vep) = (w:Tv, vb) = for any & E ((((R)))S(R")) But the denseness of "R(O(M"))SS(EE in S(R") implies that T. = O.
Using lemma 2 and (Tv)rr T we obtain T E G%- Thus, proved that T₁ T E 9%- Since I(f+) & I(fi) + R(f.)
& inf fEot
I(f) + and I is continuous, we finally get that ((T) -
inf i..,, E Arg min T(ʃ). fEe
fest
D.2.2 PROOF OF THEOREM 9
Proof. Suppose f E Arg min I(f)B i.e. f E Be and I(f*) = min I(f). Since f E B %59
then there exists sequence {st C Ge such that Te +" f and Tr Mdt oc.
Let us define s E ((RR)) as
T Tg. G
Since lim R(s))) = 0 (lemma 4h there exists 0; > 0, such that R(s)) < whenever < S or. Also, by definition Tr M, l Tr Mag. Therefore, there exists of 0, such that Tr Mg; Tr M.+ +1 whenever 0 a o
we set of min o }}. then a sequence (s } S S(R)) satisfies:
lim R(5;;) ) 1-000
Tr Ma 00
and (using lemma]
Toi; - l
Due to the continuity of I we have
lim ]
1-000
Now set fi= s'a A; (R())
and obtain the needed sequence:
lim T(A:) = lim T(J+) + X,R(J)) = T("*), lim \; = +00 ++000
{0000
+4300
where Tr Mr is bounded. remains to check that our sequence regularly solves (7), i.e.
lim inf +(f)+A,
T(f*) (this will imply lim +(, inf T(f)+
SES(RM)
SES(R")
X,R(f) 0). The inequality in one direction is obvious.
inf T()) + XRR()) & inf inf T(fa) + A, & inf lim T(fa) + Xe = fE9. 020
FES(2")
SE9x 0-+
20

---

Under review as a conference paper at ICLR 2021
inf I(f) I(f*)
Let us prove the inverse inequality.
Since rsol (I(f), R(f)) # 0, there exists {fit C S(R") such that:
((.)) + R.((()) - inf 1()) ( ++++)((f) Es t00 SES(R*)
lim d t00, lim E = 0, Or <0 1++00
and + F. From theorem 5 we obtain a E Arg min I(f). SELL
One can always find subset {da.) S {Ao} such that da, Ais dd, oo and obtain:
inf SES(R2)
T(ʃ) + A() 2 inf SES(R()
T(ʃ) + Xa R()) z
I(fa.) Na, R(Fa.) cd. 2 T(fa,) cdi
Therefore,
lim inf 2-000 fES(R")
T()) + AR() 2 lim I(fa) Ed, =I(a)= inf I(f)= ((}*)
10000
fEG,
This proves that {fr} regularly solves (7) and lim f: = i... f E rsol R(f))
E PROOF OF THEOREM 5
Again we will prove a more general statement.
Theorem 11. If M a proper and a real-valued kernel, ((M) is bounded and Tr My < cc. then Sy E B(L)(R"), and S Mf. Moreover
min R(f) SEB(LG rank S5k
||s - S112
and the minimum is attained at S = S where Pj -i uau, and {u;} are unit eigenvectors
of Mj corresponding to the k largest eigenvalues (counting multiplicities).
Proof. The boundedness of Sy follows from the Cauchy- Baanhy-Sccwarziine inequality
[5)()) |R (v (OM v; f) (VO(M) () Vo(M))
and therefore:
| < Tr ||S Σ
Lee. we have checked that Sy is bounded.
By definition, S n + L5(R") x L)(R") and (u. Syl01, Φ2)) (s)(u), @2)). u
R.a [11, Φ2) I5(R") [T()).. Let us denote fi Ref.fz Im easy to see that the following operator satisfies the latter identity:
Ouu (x)( u), (x(M)] ul]
Since the adjoint is unique, then 57 O. Let us calculate SyS}:
VO(M) VO(M) S4 (aaf1e))(y00000110000)))11111 ((1fa(x)) (ONM)]/ONN)))) (e2.f ottp:///.w...aaaa/./.aa.........
(x. (
21

---

Under review as a conference paper at ICLR 2021
Σj- (x
[Re (a. Myu
E (nn.f (xM) u]}
Thus, SyS} My. Since Tr S,S} 00 and || (u, My we obtain sj is bounded operator.
Let U1,--0, Un be orthonormal eigenvectors of M, = SyS} and A 2 ... 2 An' > 0 be corre-
Vector v; corresponds to S/[u sponding nonzero eigenvalues. For o; V let us define oi pair of functions:
Vi (VO(M) u.), (x) u,]] E L5(R") x (İ(R))
easy see that Vi, Var is orthonormal basis in Im St. and f can be expanded in the following way:
n
S) =
v.v
i=
and therefore, SVD for Sy is:
v S =
i=1
By the Eckart- Young Mousky theorem (see Theorem 4.4.7 from & Eubank (2015)), an optimal || S11: is defined by truncation SVD for Sy at kth term, i.e.:
in
min
SEB(L)((")."").ann S5k
S Σ u.v P,S
(16)
where P, -E u is projection operator to first & principal components of My. Moreover, ||S P/ Ei=k+1 of ||M R(f).
F THE ALTERNATING SCHEME IN THE DUAL SPACE FOR M (xx = C(X - y)
When M(x, y) {(x - y), the alternating scheme allows for reformulation in the dual space. By this that in Scheme we substitute Oe for the original Φv. If the primal Scheme 1 deals with operators Sz, S the dual version deals with vectors of functions / (àà substitution is based on the following simple fact:
- The
Dx ax
Theorem 12. If (x, y) <(x - v), E C(R²) and Vx ((x) O, then there exist constants Cı
and such that ||S Pt- 112 = cill 11 P
dd.-1
||2 (R.)) and (aaf)
8x
ai ai 02. (a,, dz, - (R²)
Proof. Let f _2a C be such that ||==f(le2(R)) < 00.
((M)[ C* + (0( ~ Eه + F {VoLM))w)} ~ VcB +
Sy[p], = Re (r,f, Vo(M) ~ Re (F {oif) F (VO()))1))
Ri of
( Rc (i 05 o0 ())
dr,
22

---

Under review as a conference paper at ICLR 2021
Since Sy(b)= = Re ((S)),,)) ~ Re ((S))s, 18), we obtain
(8)) KV of
(17)
Or.
where x: is constant.
Let us now introduce a vector of functions Vy = [(S)),, .. (5))ml E L2(R"). Using 17 we
obtain (Sy), = =k/ā af 0ar
and therefore:
V, - vi=mVi
of
0x
Thus, the expression ||S PI-1Só,- 11 in the alternating scheme can be rewritten
/
VIVe - Po-1Von [izee ~ 1
- PL-INV
1
0x
8x
-01101111 ll2 (122)
0
-P--1-
0x
0x
The matrix My can also be calculated from f using the following identity:
aî
_ 0s
ae as
(at. = (a:f. (xzf))
Ox: Ozz" Ox::" dx,"" L (R*)
Let us introduce a function i such that I(f) = I(f). Then, we see that all steps in Scheme can be performed with Ot rather than with Φs, using the algorithm 2
Informally, the dual algorithm works as follows: at each iteration t we compute function år adapt- ing it to data (the term f())) and adapting its gradient field to the rank reduced gradient field of the previous 6.-1. For a sufficiently large T, it will converge and Φr = Ôr Then, the second aár tx
term the last step will be approximately equal to a|| Dér: ix
llz I12 (R=) enforcing Pr-1
BdT Pr-1 8x
BØr for random x |K|] σx
?_ Thus, gradients
aer lie in a k-dimensional subspace 8x
col Pr-1. This last property is a characteristic property of functions from Fx.
Algorithm 2 The alternating scheme in the dual space.
Po 0, Φ0 0
for t = 1, 1....................... T do
d. arg min ī(0) + X|l 31 P--
0d,
ll2 (R
ix
Calculate Mr Re < aà. Dē:
azx omz 11 ( (3n)
Find [vr]T s.t. M,,. A. A1 2 2 An
P Ei V.V.
Output: VI,".,
G A NUMERICAL ALTERNATING SCHEME FOR MMD
G.1 STRUCTURE OF F(Px]
From theorems and 2 C In fact, famous theorem of Bochner 1932 gives that the Fourier transform of any positive finite Borel measure a continuous positive definite function.
That is, if f F(P), then for any distinct ",,"", TRth the matrix [f(y; yy)l,,=IIm is
23

---

Under review as a conference paper at ICLR 2021
positive semidefinite. Since (k(R)) 1, we additionally have f(0) = 1. Let PDF denote the set of all continuous positive definite functions on IR" and
M& = {f E PDF V& E RR 9 R + s.t. f(x) - glvl x..,,,,.,, v.....vv xxx 1) (18)
Thus, the following characterization of F[Px] becomes evident. Theorem 13. F(PA] = Ma.
G.2 THE DUAL FORM OF MMD
Let us define another gaussian kernel 7(x) e-- 1/ Flk]. Let Patsa (x) denote the charac-
teristic function of the random vector Xdata ~ Pdaata By definition, Pdata(x) Ele
N E ei Thus, Pdata [Mduta and Mdata x F(paata].
Using the isometry property of the Fourier transform for L2 (RR) and the convolution theorem, we that:
dMDD dMMD(H, ||k μ h (ll x (F]ll())) - F(x(()))
Thus, from Theorem 13 we obtain that the task66la is equivalent to:
Pdata - 2 21l22,2(80) + min geMMz
(19)
G.3 ALGORITIMM FOR MMD
Let II I& - +00} and M2 Fk > {1, +00} be simple penalty functions:
(Ø) i, Φ Px and II, (Ø) oc, otherwise
Mz if & M₁ and M₃ (Ø) DO, otherwise
Then, the task is equivalent to:
T(Ø) BhADD (Jdata: (T) inf
9E9i
From the result of the previous section we see that if (Φ)) = I(δ), then:
[(क) Pdata
011 2, M
Thus, the Algorithm 3 is adaptation of Algorithm to MMD.
Algorithm 3 The alternating scheme in the dual space for MMD
Po 0, 40 < 0 for 1,... do
gt -- arg min SE. 7(x)) (x) - g(x) + A AZz ((x)] 0q ge-Ms
Pa-
Sqc-1
112dx
8x
Fx
Oyc 8qc Calculate Mr (az dzy L2 (RR)
Find {v:}T s.t. M.V; \vv; X An Pt - Σi= V.VV" Output: C = span(Vi.."" VR)
If the function Pdata is real-valued, then only real-valued functions can appear the Algorithm 3 This assumption can be satisfied by adding reflections of initial points to the dataset (after it was centered).
At step 1, we search over given in the following parameterized form:
(x)
Σ accos(wf (wfo x)
(20)
24

---

Under review as a conference paper at ICLR 2021
where a; 0 and Ei= a2 1. In our implementation, set [a=]j= softmax( (uul)ii.. and u, 's are unconstrained. The number neurons in a single layer neural network with cosine is hyperparameter. Let us denote parameters {w; km by 8. It is easy
activation function, to see the function Qo is positive definite. Moreover, using Theorem from [Barron (1993). it can be shown that set of all such functions, i.e. the convex hull of {cos(w" x) RR", J. is dense in set of real-valued functions from My. Though this parameterization is quite natural, finding architectures with more expressive power in space of real-valued positive definite functions is an open problem.
Now, to minimize
(0) = / (Pxtta)d dxx + A / {(x))
Ogo 0x
0gu.-1 dx
P:-1
1|xx
with stochastic gradient descent methods (in our case, the Adam optimizer) we need to have an unbiased estimator of
040 Voll
0g0,- (z')l
V (0) OC Egnn? Volpdata (z))) z
0x P
(z')
0x
where ~ f denotes that the random vector is sampled according to the probability density (x) function Jam [xx)xxx Thus, natural estimator of the gradient is:
(g)().)
Ogo, (E.))
V Pdata - q0(z))
V Vell
Po-1 δx
8x
1=1
where {2.)/=1 72 and
ym ( ,dd Ç.
The last important issue with the practical numerical algorithm is the calculation of M, at step 2. It easy to that:
dqt M₁ = E x- 0x Ox
(x) da
(x)²
In practice we sample XI. Xi ? and estimate M, as follows:
1 Og:.
Oq.
(xi)
(x
Ox
Ox
The details of the numerical algorithm are given below ?
In all our experiments with MMD we set
Algorithm The numerical algorithm for MMD. Hyperparameters: ,. h, m, a, 32, nn.
Po + 0, 0, + 0 for 1 = 1,....777 do while 0 has not converged do
Sample {2 mm ind 72
Sample {E, 1m Rid
Σn Ipiata (z:) - qo (2:))=++ Em
@qa(E.)) dx
P.-
age, (E.))
Ux
0
Adam VoL, O, a, B1, (2)
0,
Sample (xilk=1 widd ?
Calculate Mr = 1 _!
Ago, (x.)) 0gn, (x.))
8x
8x
Find [v s.t. MAV, \, \1 2 2 \n
P Ei V.V!
Output: " Vk
25

---

Under review as a conference paper at ICLR 2021
H A NUMERICAL ALTERNATING SCHEME FOR HM
H. THE DUAL FORM OF HM
Due to well-known relationship between moments of the probability measure f and its character-
istic function p. i.e. imig
0'p(0) the task 7 is equivalent to: dæ., -u.,
8 Paata (0)
07q(0) 2 min dai Ox. gEMs
Σ
(21)
Jri dz,,
I5h, Sn
Note that the maximum mean discrepancy distance and the distance based on higher moments are substantially different. Indeed, even we set / large value (which makes & 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 21 it clear that dμM (Hdata/ 2) is degenerate case of weighted Sobolev norm between characteristic functions of Mdata and v.
H.2 ALGORIHTMS FOR HM
Analogously to the case of MMD we see that the task 7 is equivalent to: (H Φ) (lx(d)
[()) =
+ inf 6ESx
and
A
Or Pdata (0) 2e 00)(0)) ar. ... -O,,
[(0)
Σ
Σ
(x)
n2
Dz, ... ..,r,,
Sn
Thus, the Algorithm 5
is an adaptation of Algorithm to HM.
Algorithm 5 The alternating scheme in the dual space for HM
Po 0, 40 <- 0
for ,, do
qt - arg min E1 gEMi
de Eish Sn De. Ue,
Patta
0'q(0) D2,, "0
2 + A fg {(x) 80
mml 12
P. 0q 1
11
Buc Soc, dz, L,, (Bx)
2 Caleulate Mj
3 Find [v.}1 s.t. M.V. = X d 2 ... An
4 P. Ei= V.V? Output: C = span(Vv..".. Vk)
Again, as in numerical algorithm for MMD, at step 1, we search over given in the form 20 The objective of step can be represented as:
8 (paata dz,, Ōx;,
qo (0) 12 XE 0x
oqo (z') P.-
Oqe__
(0)) = Σ A.E (1)))))))))))
(z') 11)
0x
where (1,, n) is the discrete uniform distribution over {1.... ... n.. descent methods we need to have an unbiased estimator of (( which equal to:
). To apply the stochastic gradient
8 (Paata - (0) 2 + ĀE
Oqo
oqa.-1
A. Eg ((
(z') P
(z²)113
Or,, .. Oxy
Ox
0x
Thus, a natural estimator of the gradient is:
mi
8 (Paata - (0)
ma
Oqo(E:))
Oqe. (E.))
+
P.-
Ir 1f 1. ... 05
Ox
8x
where {a[s, é, .1. lowing Algorithm
I. rid u(1, n) and {m iid Č. Overall, we obtain the fol-
26

---

Under review as a conference paper at ICLR 2021
Algorithm
The numerical
algorithm
for
HM.
Hyperparameters:
X, ( 1,4: m1, m2, a, 81, 82, nn.
Po 0,00 4- 0 for t = 1, 1,, , do
while / has not converged do
Sample {als, i, ==
Rid U(1, n)
I.
Sample {5:) 12
.d
L
Σ. mi
i=1
8 (Pdata (0)
(E.) /x
ma
P 0gu (E,))
ix
0 Adam VoL,0,a, 82)
0g
Sample txil/=1 rid ?
Oga, (x.)) 0ge, (xi)) 0x
Calculate Mr
8x
Find {v:]T s.t. M.V; A.vi, A1 An
P Ei=1 V.V.
Output: V1
A NUMERICAL ALTERNATING SCHEME FOR WD
By Theorem 6. the task 13 equívalent min W(u, pdata), to the following task:
[ + inf
dE
where I(T,) W(,, (data) if μ Ps and [()) oc, if otherwise. The alternating scheme 1 is designed to solve the penalty form of the problem, i.e.
(()) + AR(o) + min
$ES(2*)
which is equivalent
( (data) + AR(
min (E5)
where (R S S(R") is set of Schwartz functions that can serve as pdf: o(x) > O. Jan o(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
1.1 HOw TO MINIMIZE OVER Φ?
In the case WD, the minimization step of the alternating scheme makes the following:
Φ - arg min (6, Mdata) + Al|S - Pt-1S (E5)
(22)
where Sy = /O(M)
For numerical implementation of that step we need choose some family of functions that is dense in (PP) (or, rich enough to approach the solution n²)). Following the tradition of GAN research let us assume that the family is given in the following form
H { is pdf of random vector (z( p(z), E Q)
(23)
where {gol0 E θ} a parameterized family of smooth functions (usually, neural network) and Ballowing/Arrovsky
p(z) is some fixed distribution (usually, the gaussian distribution). Following Arjovsky et al.) (201)). we make the assumption] [] numerical algorithm we need an access a procedure that samples according to dh ()), not the function itself.
Assumption 1. llger (z') go (z)ll < L(e, z ( 0|| + ||z - z||) where
E Ez~p(1) I(0, x) +o0
II C S(R") not satisfied, then choose H. (No G. O) for very small
27

---

Under review as a conference paper at ICLR 2021
Thus, instead of solving 22 we solve:
Φr arg min W(6, (ddaa) + A11s Pt- 11 PEH
taking into account that di-1 E H.
The Kantorovich-Rubinster duality theorem gives us that:
W (0o, (ddaa)
max Exc-pass [f(x) R ( (()o())
fill
which turns 22 into the following minimax task:
arg min max Ex Ex--pan Alls, Pt-1So._1l? (24) PE fillll
In practice, we choose a family of functions L = {fw e E W) and internal maximization is made over w W with an additional penalty term that penalizes violation of the Lipschitz condition: Vx ||ʃx|| < 1.
family of minimax algorithms for the minimization of W ( (omm)))mm)) was developed in series of papers Arjovsky et al. (2017): Gulrajani et al.] (2017); Wei et al. (2018). The standard minimax scheme that gained popularity GAN literature iterates two steps: a) Niter times make a gradient ascent over E W, b) make a gradient descent over 0. The task 24 can be viewed as Wasser- stein GAN with an additional regularization term (T where T(0) 1|S
P 11 To
adapt these algorithms to the minimization of our function, we only need to have an unbiased esti- mator of the gradient ar 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 fw)) can be set in standard fashion
choose Petzka al. (2018)'s version, in which the term max (0. || f (Ex + {))(())) 1]2 enforces Lipschitz condition (see step (*) of the Algo- rithm
Algorithm 7 Numerical algorithm for WD. We use M(x, = l1x-y112 and default values of 10, A 100, neritic 5, 40, 10000n, 0.00001, 0.5, 82 0.9
Po + 0, 0 10
for t , do
Minimax realization of min W (coo Memp) + (T (*):
0
while 0 has not converged do for 1, 1.... Deritic do Discriminator updates w
Sample {2, 1T {2/) p(z)
for (z,)) + E., 2(0.2.2))
L レ 0 - Adam VoL, O, a. 32,
is defined in equation 25
0,
Realization of step (**);
Sample {2,)/en [1 p(z)
M,
Eis 90. (z,)go, (z MM M (z,) go. (z,)) 12 s.t. M.,; A, d 2 ...
Find {v:}]
Am
Pt 4 E/ v.v.
Output: VI Vk
1.2 HOw TO ESTIMATE 30
8T AND M ?
Another important aspect of the numerical algorithm is the complexity of estimating the matrix Moo, at step (**). The following theorem shows that we only need to sample p a sufficient number of times to estimate dø
8T and M
Theorem 14. the pdf of the random vector (z), p(z), then
28

---

Under review as a conference paper at ICLR 2021
0T
0E(0, 2, z')
Exxx'p
00
aB
Mè. Bu,x' ~p (z)go (z')² MM (g0(z))gozz)) go
where
[(0, x,,') = (ge go (z))) (z(())
(25)
((go(z) P1-190. (g) (z'))
and RHS is well-defined.
1,2.1 DEFINITION OF H
Specifically, for robust PCA/outlier pursuit applications, we define (x(x) as a probability density function of the random vector a + b, where a, b are independent and a is the i-th column of matrix 01 E RaxN (where ~ (1(1, N) is sampled uniformly from {1, (1,..., N.)) b g0, (c), с (O,,, In) and g0, : :R" + RR" is a neural network with weights 02. Thus, 0 (0,, O2),
can be checked that H. defined in this way, satisfies the Assumption [ We specifically introduce the random vector a here because, according to Theorem 6 the ultimate solution of the problem corresponds to 01 Y and b 0. This guarantees that the solution is approachable from set H.
1.3 PROOF OF THEOREM
We need to following lemma
Lemma 5. ||S PS = Exy~o(x x) + Exy-p(x Py)M(x,y) 2xx (x,
Proof of lemma.
||S 100--///|11110000000021-------00
- - (PX)))))dxx)(((sss(()))
Exy M) + Exyrv (x) - 2F 2Ix (x (()))))x)))))))
Proof of theorem [4
Using lemma we have:
T(0) = Ex 0o (x y) + Ex,y (x P.-
2 2xx (x (x,
Ex (g0 gn (z) ) (2) go (z())))
F Bz,zz'~pp 90. (z) Pt 190
(z), (2'))
2E 2E, x. up (0u
(z((z)) (z) P. 199.., (z) )) (2), 90. (z'))
The second term does depend on 0. Therefore,
0T a 00
E
00
where
E(0,2, (go(z)) go((z))) (go go(())) 2((g(z)) P2-190,-1 (z')) (60)))))))))) (2))
29

---

Under review as a conference paper at ICLR 2021
If Zzp (0( well-defined (the proof of sufficiency of that condition similar the proof
of Theorem 3 from Arjovsky et al. (2017)), then, using Leibniz integral rule, we obtain:
0= 2, xz)
E xx E(0, = Exxx' mp
00
00
The fact that
Mè. (px (go go()))
obvious from the definition Møy Ex xy" y).
NUMERICAL ALTERNATING SCHEME FOR SDR
For binary classification case, given labeled dataset ((x; H x, E R" y: EC, C {0.1) we formulate the sufficient dimension reduction problem as the minimization task:
J( = (2 == Me,
+ min fEF
where L(c,y) = - log(v) (1 - log(1 (og(1-y) y)
We apply the alternating scheme in the dual space (Algorithm 2 to this task. We set M(x <(x - y), where & is strictly positive probability density function. A numerical version of the scheme is given below (Algorithm 8)
At every iteration = 1, of the Algorithm 2 we solve the task (in our case i = J):
Oá
d
a arg min ((0) + XlI
Pa-
0x (12
In numerical version of the algorithm we assume that given as neural network fo. i.. our task becomes:
05e
(E) - P.-1 0f0.-1 0f0--1 (2)112
0g - arg min (Se) + λE
E
0x
0x
The gradient of the function (0) = ((fo + XI
ale (5)
P.-
al6,11 8x
(E)112 equals:
g~:
ex
@+( 00 Paaa))
θ [c, fo + c)) + ÅE g-
05o
Ofe... (E)
(E) P--1 Ox
0e
ag
Ôx
That is why VoL (given Adam optimizer in the gradient descent loop) in the Algorithm is an Thus, the "while loop" we find optimal Φz = for-
unbiased estimator of 04(8)
According to Algorithm that
the next goal to estimate M, = Ro
De: ad, (ar; Oxy 1.2,
easy to
aô. (x)
ao. (x)" E Ex 8x
Ofe.
Of..
M, = E x 8x
(x)
(x)
Ox
8x
From the last we see that the matrix Mj can be estimated by sampling X ~ a sufficient number of times (the parameter I in our algorithm). All the rest is identical to Algorithm 2
The regression version of the algorithm can be obtained by setting L(c, c') = (c c²)2. Implemen- tations for different databases be found at github
30

---

Under review as a conference paper at ICLR 2021
Algorithm The numerical alternating scheme for SDR. We use v = 1.0, {(x) Gö (x) and
default values of = 10, m = 50, m' = 100, 1 30000, 0.0001, 31 = =0.5. 82 0.9
Po 0, 0o 0
for f = ................................ , do
while 9 has not converged do
Sample [(z c))] Paata
Sample [= v=,)
Sample {c,
1m ;= ~
Σ. Mc,,folz + Ei))
En'
0f.(E.))
P-1
0fa, (E,))
2
m
8x
8x
0
Adam VoL, 0, a, B1, (2)
02
Sample (xille1 ~ ?
Calculate M, = + 12t Find {v:]T s.t. M,V, X, d 2 ... ...2 An
0fe (x.1) 0fa, (x.1) 8x
8x
P E V.V!
Output: VI," Vk
31


Reviews:
[{"review_id": "boK0P0c-ADa", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper considers the unsupervised dimension reduction problem. That is, given a set of points in R^n, find a low dimensional affine subspace that approximate the support of the distribution that generated the points. More specifically, the paper considers the empirical probability density function p_emp of a dataset which is the average of \\delta^n(x-x_i), where x_i's are points of the dataset and \\delta^n is the n-dimensional version of the Dirac function. Then the goal is to find a distribution q such that its density is supported in a k-dimensional affine space and it minimized a certain loss D(p_emp,q), where D is a measure of distance between two distributions.\nThe paper then presents 4 examples of problems that can be formulated in this framework: 1) maximum mean discrepancy; 2) distance based on the higher moments; 3) Wasserstein distance; and 4) sufficient dimension reduction.\nFinally, the paper proposes an alternating optimization scheme to solve this optimization problem and presents experiments that compare the accuracy of the proposed method with other dimensionality reduction methods like PCA.\n\nI think the paper is not well-motivated, and it is not clear what are the novelties of the paper. Please explicitly state what are the contributions of the paper. The experiments are also very inconclusive. Table 1 reports the accuracy of KNN on 2 and 3 dimensions. First, I think it is better to report the reconstruction error of PCA and other methods instead of this. Moreover, it is better to test the projection on a bit higher dimensions as well. For example what happens for k=10?", "extra_scores": {}, "extra_text": {}}, {"review_id": "OgIiL2D70bU", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper considers the unsupervised dimension reduction problem, in which one is given a finite number of points in R^n drawn from some distribution and wants to find a low-dimensional distribution that approximates the underlying unknown distribution. The paper proposes to solve a minimization problem min_f { I(f) + lambda*R(f) }, where I(f) is the distance between f and the empirical distribution of the data and R(f) is a penalty term that tries to force f to be low-dimensional. When the penalty is defined via a kernel function, R(f) admits a highly tractable form and, for a given f, can be solved by an SVD. This leads to an alternating scheme to minimize I(f) + lambda*R(f), which is the main contribution of the paper. The algorithm is iterative. In each iteration, it finds the minimizer f with respect to the distance to the empirical distribution and the penalty in terms of a k-dimensional subspace found in the previous iteration, then finds a new k-dimensional space by minimizing the penalty. The paper then conducts experiments by applying this general theory to specific distance functions and kernels.\n\nStrengths:\n- Theory stated in general normed spaces, clean and neat\n- Algorithm is conceptually clear and simple\n- Experiments seem to confirm that the proposed alternating scheme is a serious competitor with the existing Sliced Inverse Regression and Kernel Dimensionality Reduction algorithms.\n\nWeaknesses:\n- It is not clear how one should choose a finite-dimensional domain of phi (which is defined in an infinite-dimensional space). Although the authors have specified some domains for specific problems in the Appendix, in general, it is not clear how to choose the domain phi and how to solve the optimization problem in Line 3 of Algorithm 1. It is also not clear when to use the primal and when to use the dual form of Algorithm 1. Some discussions are expected to be included in the main body.\n- No convergence analysis of Algorithm 1?\n- It would be better to compare the running time with the existing algorithms, too.\n\nMinor points:\n- page 2, end of Section 2, “Identity matrix” -> “The identity matrix”\n- page 5, two lines above Theorem 4: “a real part” -> “the real part”\n- page 5 and 6, “task 11”, “problem 11”, etc -> “11” should be “(11)” so that readers know it refers to Equation (11).\n- page 6, the paragraph below Algorithm 1, “Scheme 1” -> “Algorithm 1”?", "extra_scores": {}, "extra_text": {}}, {"review_id": "5jfOUrYDwEr", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Summary: In the dimension reduction problem, we are given a set of high-dimensional points and would like to embed them in a lower dimensional space so as to preserve relevant properties. This paper proposes a certain optimization framework for dimension reduction, and suggests to solve it by alternating optimization. The problem is formulated in terms of a general distance measure between distributions, and the alternating optimization scheme requires instantiating a certain optimization step for the specific choice of distance measure. The paper instantiates it for four such choices. The method is implemented and some experiments are presented.\n\nComments:\nMuch of the paper consists of highly technical mathematical derivations, and unfortunately I do not have the background to assess this content. As a matter of presentation, it seems to me an unfortunate choice to not even define the basic notions on which the paper is based (Schwartz space, tempered distributions, generalized functions etc); instead, the paper just refers to textbooks. I do not believe the current manuscript is widely accessible to ICLR audience, as it seems to require a rather specialized background in functional analysis. An introductory section or appendix defining the basic terms and facts is standard in such cases, and would go a long way in making the paper accessible and self-contained.\n\nThe experiments do not show a clear advantage of the proposed methods; in fact, when compared with two other baselines, each of the three methods gets the best result on exactly one third of the experiments in Table 1, and the differences between them are generally small and doubtedly significant. Is there a further message in these results that point to some advantages of the proposed method?\n\nThe real datasets used are mentioned by name and with a link to their UCI source. This seems insufficient; please include a written out reference or link for them, as hidden hyperlinks are not visible and are of course lost when the paper is printed. It would also help to specify the parameters of the datasets (number of points, ambient dimension etc), which are needed in order to put your empirical results in context.\n\nI could not find a discussion of the running time / scalability of the proposed method compared to PCA and the other baselines. It would help if the authors could comment on that (both asymptotically and in practice) as it is relevant to assessing their claim that their method is suitable for use instead of PCA. \n\nAs a final remark on presentation, please use the \\citep command where appropriate. Currently most of the references are unbracketed and make the text difficult to read (e,g., \"the problem becomes of special interest Cunningham & Ghahramani (2015)\", and many other such instances).\n\nConclusion: The mathematical content of the manuscript is largely opaque to me. The experimental results are not outstanding, though perhaps the new approach presented here has some conceptual merit that could justify acceptance (this is currently difficult for me to judge). Given the diverse spectrum of ICLR target audiences, I would advise revising the presentation to be friendlier to the general ML community.\n\nPost-rebuttal update: I thank the authors for their response. The points raised above were largely acknowledged and the authors chose to not revise the manuscript, so my assessment remains the same.\n", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: Overall, there were significant concerns about the motivation and experiments in this paper, and these were thought not to merit acceptance on their own. Because of this, the reviewers started discussing the theory to see if that would justify acceptance. The reviewers were not able to find a clear advantage over existing approaches, nor sufficient motivation; also the presentation was found to be largely inaccessible. In the rebuttal there was a brief mentioning of background and possible implications, but they were hard to assess and the paper itself did not have such context nor was updated to have such context. For a future version, one recommendation could be to focus significantly more on context, motivation, and improvements over prior work. Also, making the paper more self-contained could help. 
Decision: Reject