# Fourier analysis perspective for sufficient dimension reduction problem

A theory of sufficient dimension reduction (SDR) is developed from an optimizational perspective. In our formulation of the problem, instead of dealing with raw data, we assume that our ground truth includes a mapping f: R^n→ R^m and a probability distribution function p over R^n, both given analytically. We formulate SDR as a problem of finding a function g: R^k→ R^m and a matrix P∈ R^k× n such that E_ x∼ p( x)| f( x) - g(P x)|^2 is minimal. It turns out that the latter problem allows a reformulation in the dual space, i.e. instead of searching for g(P x) we suggest searching for its Fourier transform. First, we characterize all tempered distributions that can serve as the Fourier transform of such functions. The reformulation in the dual space can be interpreted as a problem of finding a k-dimensional linear subspace S and a tempered distribution t supported in S such that t is "close" in a certain sense to the Fourier transform of f. Instead of optimizing over generalized functions with a k-dimensional support, we suggest minimizing over ordinary functions but with an additional term R that penalizes a strong distortion of the support from any k-dimensional linear subspace. For a specific case of R, we develop an algorithm that can be formulated for functions given in the initial form as well as for their Fourier transforms. Eventually, we report results of numerical experiments with a discretized version of the latter algorithm.

## Authors

• 13 publications
03/12/2019

### Dimension reduction as an optimization problem over a set of generalized functions

Classical dimension reduction problem can be loosely formulated as a pro...
02/24/2022

### Numerical reconstruction from the Fourier transform on the ball using prolate spheroidal wave functions

We implement numerically formulas of [Isaev, Novikov, arXiv:2107.07882] ...
07/17/2019

### Interesting Open Problem Related to Complexity of Computing the Fourier Transform and Group Theory

The Fourier Transform is one of the most important linear transformation...
08/05/2020

### A fast algorithm for the electromagnetic scattering from a large rectangular cavity in three dimensions

The paper is concerned with the three-dimensional electromagnetic scatte...
01/22/2021

### An in-place truncated Fourier transform

We show that simple modifications to van der Hoeven's forward and invers...
03/30/2018

### An efficient high dimensional quantum Schur transform

The Schur transform is a unitary operator that block diagonalizes the ac...
01/22/2020

### Substitution of subspace collections with nonorthogonal subspaces to accelerate Fast Fourier Transform methods applied to conducting composites

We show the power of the algebra of subspace collections developed in Ch...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The dimensionality reduction

is an important problem in data science that has many facets and non-equivalent formulations coming from different contexts, either purely mathematical or appearing in applications. The classical one was first formulated in the work of R. Fisher

[3] and currently known as

. Subsequently, the idea of principal components was applied to more general frameworks, giving birth to new branches of statistics/machine learning such as the manifold learning (e.g. the nonlinear dimensionality reduction) and the sufficient dimension reduction. In the manifold learning formulation (which is the direct generalization of the classical) we are usually given a finite number of points in

(sampled according to some unknown distribution) and our goal is to find a “low-dimensional” geometric structure that approximates “the support” of the distribution and satisfies some additional properties such as smoothness, low complexity etc.

Unlike the latter formulations, in the sufficient dimension reduction (sometimes called the supervised dimension reduction), we are given a finite number of pairs

, also generated according to some unknown joint distribution

, and our goal is to find vectors (where ) such that symbolically:

 y⊥⊥x|wT1x,⋯,wTkx

The latter means that an output is conditionally independent of , given . Or, that conditional distribution is the same as .

Of course, the latter formulation can hardly be solved if we do not make any assumptions on the joint distribution, or more specifically on the conditional distribution . A standard assumption is the following semi-parametric discriminative model:

 y=g(wT1x,⋯,wTkx)+\boldmathε (1)

where is a Gaussian noise with and . The function is an unknown smooth function. Then, the function is called the regression function.

There are 3 major methods to estimate parameters of model

1: (1) sliced inverse regression [6],[2]; (2) methods based on an analysis of gradient and Hessian of the regression function [7], [13], [8]

; (3) methods based on combining local classifiers

[4], [11].

Probably, the closest to ours is the second approach. Let us briefly outline its idea for . According to that approach we first recover the regression function and estimate the distribution from our data

. The former can be done by solving the supervised learning problem using any suitable model, e.g. by neural networks, and the latter is typically done by assuming that

and estimating the parameters

of the multivariate normal distribution. At the second stage we no longer need our data and treat

as the ground truth. Since, for recovered it is natural to expect that , then a natural way to reconstruct vectors is to set them equal to first principal components of the matrix , where is a Hessian matrix of at point .

In our paper we also assume that is an already given ground truth, though unlike the previous approach, we formulate the main problem optimizationally, i.e. our goal is to find

 Ex∼p(x)∣∣f(x)−g(wT1x,⋯,wTkx)∣∣2→ming,w1,⋯,wk (2)

It is easy to see that the latter corresponds to the maximum likelihood approach to estimating of the parameters . Since is an infinite-dimensional object, we analyse it by the tools of functional analysis, specifically using a theory of tempered distributions. The key observation of our analysis, stated in theorem 3.3 of section 3, is that a class of functions of the form can be characterized as those functions whose Fourier transform is supported in a -dimensional linear subspace. Instead of optimizing over generalized functions with a -dimensional support, we suggest minimizing over ordinary functions given in a generic form but with an additional constraint. In order to force their support to be -dimensional, in section 4 we introduce a class of penalty functions such that large values of indicate a strong distortion of the support from any -dimensional linear subspace. For a specific case of , in section 5 we develop an algorithm for our problem that can be formulated for functions given in the frequency coordinate form as well as in the initial coordinate form. The last section is dedicated to experiments on synthetic data.

## 2 Preliminaries

Throughout the paper we will use common terminology and notations from functional analysis. The Schwartz space of functions, denoted , is a space of infinitely differentiable functions such that , and equipped with a standard topology, which is complete and metrizable. A cartesian power is a set of vector-valued functions, i.e. if and only if .

By the tempered distribution we understand an element from the dual space, . The Fourier and inverse Fourier transforms are first defined as operators by:

 F[f](\boldmathξ)=1√2πn∫Rnf(x)e−i\boldmathξTxdx,f∈S(Rn)
 F−1[f](x)=1√2πn∫Rnf(\boldmathξ)ei\boldmathξTxd% \boldmathξ,f∈S(Rn)

and then extended to continuous bijective linear operators by the rule: . The Fourier transform can be applied component-wise to objects from the cartesian power which we will also call the tempered distributions.

If a function is such that for any then it induces a tuple , where .

For a measure , by we denote the Hilbert space of functions from to , square-integrable w.r.t , with the inner product: . The induced norm is then . A space (i.e. when is Lebesgue measure) can be embedded into , i.e. , where corresponds to a tempered distribution . Therefore, Fourier transform can be defined on and we will use the fact that is a unitary operator.

For the convolution is defined as . For , the convolution is defined as a tempered distribution such that:

 ψ∗T[ϕ]=T[~ψ∗ϕ]∀ϕ∈S(Rn)

where and the multiplication is defined by:

 (ψT)[ϕ]=T[ψϕ]

Both operations can be extended to the case when by applying them to every component of .

A set of infinitely differentiable functions with a compact support in is denoted as . The Sobolev -norm on for is defined as . The Sobolev space is a the completion of w.r.t. the norm .

For a matrix the Frobenius norm is .

## 3 Problem formulation

Let

be a probability density function such that

. The probability density function defines the Hilbert space , i.e. where . We are also given a real-valued function from which can be given in an arbitrary form, keeping in mind the case of

defined by a feed-forward neural network. Our goal is to approximate

in the following form (for fixed in advance):

 f(x)≈g(wT1x,⋯,wTkx)

where is an arbitrary function from and .

###### Theorem 3.1

For , we have and .

###### Proof (Proof of theorem 3.1)

It is enough to prove the theorem for . W.l.o.g. we can assume that are linearly independent. If they are linearly dependent and, e.g. , then we define . It is easy to see that and we reduced to the case of theorem for .

If and

is an invertible matrix, then

. Indeed, if we denote and , then:

 xα11⋯xαnn∂β1x1⋯∂βnxns(Ax)=(bT1y)α1⋯(bTny)αn⋅(aT1∂y)β1⋯(aTn∂y)βns(y)

and after opening all the brackets we will obtain a finite sum of expressions of the kind that is bounded. In fact, we proved that Schwartz class is invariant under invertible linear change of variables.

Thus, if we complete with to form a basis in , and make the change of variables , then from we obtain a function . It remains to prove that this function is also in .

For any the expression will be a sum if terms each of them being bounded.

Eventually, we note that and therefore .

If we choose the squared error as the loss function, then we come to the following optimizational problem:

 (3)

The problem is non-convex and the minimum is taken over infinite-dimensional object. Let us reveal the structure of the objective:

 ||f(x)−g(wT1x,⋯,wTkx)||Lm2,p=||√p(x)f(x)−√p(x)g(wT1x,⋯,wTkx)||Lm2

We can apply Fourier transform to our functions, taking into account that Fourier transform is unitary on .

 ||√p(x)f(x)−√p(x)g(⋯)||Lm2=||F[√p(x)f(x)]−F[√p(x)g(⋯)]||Lm2

Let us denote . The following statement is an application of the convolution theorem to our case:

###### Theorem 3.2

If and , then and

 Tk=1√2πnγ∗F[Tl]
###### Proof (Proof of theorem 3.2)

W.l.o.g. we again assume that . For we have:

 ||l||L∞(Rn)≤||g||L∞(Rn)<∞

I.e. . Unfortunately, is not a rapidly decreasing function, because , in general, defines a nonempty affine subspace and ’s value on the whole subspace will be constant . Therefore, the Fourier transform of is not necessarily an ordinary function.

Since ,

 Tl[ϕ]=∫Rnl(x)ϕ(x)dx,ϕ∈S(Rn)

is a continuous operator (i.e. a tempered distribution), therefore is also a tempered distribution.

By definition . Let us prove that

 Tk=1√2πnγ∗F[Tl]

Since , there exists a sequence of functions , such that

 Tϕn→Tl,n→∞\textscor∀ϕ∈S(Rn),∫Rnϕn(x)ϕ(x)dx→∫Rnl(x)ϕ(x)dx

The latter follows from the well-known fact that is dense in .

It is easy to see that

 √pTϕn→√pTl,n→∞\textscor∀ψ∈S(Rn),∫Rn√p(x)ϕn(x)ψ(x)dx→∫Rn√p(x)l(x)ψ(x)dx

because we can set in the former expression.

The convolution theorem states that for any 2 functions we have:

 F[uv]=1√2πnF[u]∗F[v],F[uTv]=1√2πnF[u]∗F[Tv]

Therefore:

 F[√pTϕn]=1√2πnγ∗F[Tϕn]

Since is a continuous operator, then and in . In order to obtain the needed result it remains to show that the convolution operator , is also continuous.

By definition where . I.e. we have to show that if

 Ti→T\textscor∀ϕ∈S(Rn),Ti[ϕ]→T[ϕ]

then

 γ∗Ti→γ∗T\textscor∀ψ∈S(Rn),Ti[~γ∗ψ]→T[~γ∗ψ]

The latter is obvious if we can set in the former expression. Thus, theorem proved.

The basic phenomenon behind our approach to optimization of (3) is the following statement:

###### Theorem 3.3

A function can be represented as if and only if there is an orthonormal basis such that:

 F[Tl]=r(aT1x,⋯,aTk′x)n∏i=k′+1δ(aTix),r∈Sm(Rk) (4)

where – Dirac’s delta function. Moreover, .

###### Proof (Sketch of the proof of theorem 3.3)

W.l.o.g. we can assume that and are linearly independent. A rigorous proof of the theorem would require a carefull checking of certain integral identitites. Instead we will present a sketch of the proof at the level of strictness common to theoretical physics papers.

() We also can assume that are orthonormal. Indeed, after every redefinition of given by the rule we get the same function if we simultaneously transform to . By making such redefinitions, we can always orthogonolize by Gramm-Schmidt process with a subsequent scaling of ’s arguments.

Let us complete with to form an orthonormal basis in and set:

Then in the Fourier transform formula we will make the change of variables , , :

 F[l](\boldmathξ)=1√2πn∫Rng(wT1x,⋯,wTkx)e−i\boldmathξTxdx=1√2πn∫Rng(y1)e−i\boldmathξTQ[y1y2]dy1dy2==1√2πn∫Rng(y1)e−i(QT1\boldmathξ)Ty1−i(QT2% \boldmathξ)Ty2dy1dy2=1√2πn∫Rkg(y1)e−i(QT1\boldmathξ)Ty1dy1⋅⋅∫Rn−ke−i(QT2\boldmathξ% )Ty2dy2=√2πn−kF[g](QT1\boldmathξ)δn−k(QT2\boldmathξ)

where . Here we used that . Thus, we obtain the needed representation.

() Suppose that:

 F[l]=r(aT1x,⋯,aTk′x)n∏i=k′+1δ(aTix)

Using inverse Fourier transform we get:

 l(\boldmathξ)=F−1[F[l]](\boldmathξ)=1√2πn∫Rnr(aT1x,⋯,aTk′x)n∏i=k′+1δ(aTix)eixT\boldmathξdx

After the change of variables , where

 O=[a1,⋯,an]

we get:

 l(\boldmathξ)=1√2πn∫Rnr(y1:k′)n∏i=k′+1δ(yi)ei∑ni=1yiaTi\boldmathξdy1:n=1√2πn∫Rnr(y1:k′)ei∑k′i=1yiaTi\boldmathξdy1:k′=1√2πn−k~g(aT1% \boldmathξ,⋯,aTk\boldmathξ)

where .

Substantively, the theorem claims that if the function’s value depends only on the projection of an argument on , then frequencies from the spectrum of such function are all in .

###### Definition 1

A set of tempered distributions of the form (4) is denoted as and called a set of functions with -dimensional support.

Thus, our problem becomes equivalent to:

 ||f′−k||Lm2→minTk=γ∗g,g∈Gk

For simplicity of our notation, let us use and interchangeably (from the context it is always clear what we mean). Thus, our problem is:

 ||f′−γ∗g||Lm2→ming∈Gk (5)

Note that if we would restrict to be any ordinary function, the latter problem is known in the theory of inverse problems. E.g., in a case when , a problem of finding such that is known as the deconvolution of gaussian kernel, and has many applications in mathematical physics [9], [10], [12]. But with our type of restriction, besides that we cannot guarantee that the minimum is attainable on a function from , the set itself does not suit as a good optimization space as it lacks obvious metrics, completeness properties etc.

Instead of minimization over tempered distributions we will relax the property that the support of the function is strictly -dimensional, reducing the problem to optimization over ordinary functions:

 ||f′−γ∗g||Lm2→ming:R(g)≤ϵ

where is a penalty term that penalizes if “the dimensionality of its support is greater than ”. In the next section we describe one natural approach to construct such a penalty term .

## 4 Penalty function

Let be a continuous function such that and . Let us consider a set of functions:

 LI={g:Rn→Cm|∫RnI(g(x))dx<∞}

We believe that practically the most interesting case is . Since , we will correspond to the finite measure function (induced by the density ):

 μg(A)=∫AI(g(x))dx,A⊆Rn,μg(Rn)<∞

on the -algebra of Lebesgue measurable sets. Any finite measure induces the probability measure via the normalization: . We will call a finite measure on a -dimensional measure if there is a -dimensional linear subspace such that .

In the previous section we proved that our problem (3) can be reduced to optimization task (5) over functions with -dimensional support. As we have already pointed out, (as well as ) lacks standard metrics on it, so we need to devise a certain way to measure a distance from an ordinary function to a set . If is an ordinary function, then its support cannot be strictly -dimensional. It is natural to define a distance till as , for a proper distance function on measures. It turns out that -dimensional measures can be characterized in a very simple way:

###### Theorem 4.1

Let be a finite measure on such that . The measure is -dimensional if and only if

 rank(M)≤k

where .

###### Proof

Let be i.i.d. random vectors sampled according to and

. A natural estimator for the matrix of second moments

is:

 1NN∑i=1xixTi=1NXTX

where .

This estimator is consistent, i.e.:

 limN→∞P[||1NXTX−M||F>ϵ]=0,

If we denote , then the latter can be shown after analysis of: . Indeed,

are i.i.d. random variables with finite second moment

. Therefore, by weak law of large numbers:

I.e and therefore:

 limN→∞P[||1NXTX−M||F>ϵ]=0.

() Now suppose that

. I.e. we can find orthonormal vectors

such that . Since , then