Approximation of Generalized Ridge Functions in High Dimensions

TL;DR

This paper introduces algorithms for approximating generalized ridge functions, especially linear-sleeve functions, in high-dimensional spaces, with proven error bounds and numerical comparisons.

Contribution

The paper presents novel algorithms for approximating linear-sleeve functions and extends these methods to general sleeve functions, including error analysis and practical implementation.

Findings

Algorithms effectively approximate linear-sleeve functions.

Error bounds are established for the proposed algorithms.

Numerical experiments demonstrate the algorithms' performance.

Abstract

This paper studies the approximation of generalized ridge functions, namely of functions which are constant along some submanifolds of $\mathbb{R}^N$. We introduce the notion of linear-sleeve functions, whose function values only depend on the distance to some unknown linear subspace $L$. We propose two effective algorithms to approximate linear-sleeve functions $f(x)=g(\text{dist}(x,L)^2)$, when both the linear subspace $L\subset \mathbb{R}^N$ and the function $g\in C^s[0,1]$ are unknown. We will prove error bounds for both algorithms and provide an extensive numerical comparison of both. We further propose an approach of how to apply these algorithms to capture general sleeve functions, which are constant along some lower dimensional submanifolds.

Equations

$${\mathbb{R}}^N\supseteq \Omega \ni x\mapsto f(x)=g(Ax),$$

$$f(x)=g(⟨x,a⟩).$$

$$f(x)=g(\mathrm{dist}(x,M{)}^2),$$

$$f(x)=g(\mathrm{dist}(x,L{)}^2).$$

$$\Vert L-\tilde{L}{\Vert }_{\text{HS}}\le C(1+K{)}^{N-d}\sqrt{N-d}{h}^s,$$

$$\Vert L-\tilde{L}{\Vert }_{\text{HS}}\le \tilde{C}\sqrt{N-d}{M}^{-1},$$

$$\Vert g{\Vert }_{C^s}:=\Vert g{\Vert }_{C^s[0,1]}:=|g^{(k)}{|}_{\text{Lip}(s-k)}+\sum _{j=0}^k\Vert g^{(j)}{\Vert }_{C[0,1]},$$

$$|g{|}_{\text{Lip}(\beta )}:=\underset{x\ne y}{\sup }\frac{\left|g(x)-g(y)\right|}{{\left|x-y\right|}^{\beta }}.$$

$$\#\{i:x_i\ge \epsilon \}\le M{\epsilon }^{-1/p},\epsilon >0.$$

$$\Vert \mathrm{sign}(\lambda )\frac{\tilde{x}}{\Vert \tilde{x}\Vert }-x\Vert \le \frac{2\Vert \tilde{x}-\lambda x\Vert }{\Vert \tilde{x}\Vert }.$$

$$\Vert P{\Vert }_{\text{HS}}=\sqrt{\sum _{i=1}^N\Vert P{e}_i{\Vert }_2^2}.$$

$$\mathrm{dist}(x,P):=\underset{y\in P}{\min }\Vert x-y{\Vert }_2.$$

$$\Vert g-Q_{h}g{\Vert }_{C[0,1]}\le C|g{|}_{C^s[0,1]}{h}^s,$$

$$f(x)=g(\mathrm{dist}(x,L{)}^2)=g(\Vert P_{P}x{\Vert }_2^2),$$

$$\left\{x\in {\mathbb{R}}^N:\mathrm{dist}(x,L{)}^2=\mathrm{dist}(x_0,L{)}^2\right\}$$

$$\nabla f(x_N)=2g^{\prime }(\Vert P{x}_N{\Vert }_2^2)P{x}_N,$$

$$h^{N-1}(\hat{x}):=f(V\left[\begin{array}{c}\hat{x}\\ 0\end{array}\right])=f^{N-1}(V\left[\begin{array}{c}\hat{x}\\ 0\end{array}\right]),$$

$$\left[\begin{array}{c}\nabla h^{N-1}(\hat{x})\\ 0\end{array}\right]=V_{N-1}^T\nabla f(x),$$

$$\nabla f^{N-1}(x_{N-1})=V{V}_{N-1}^T\nabla f(x_{N-1})=2g^{\prime }(\Vert P{x}_{N-1}{\Vert }_2^2)V{V}_{N-1}^{T}Px.$$

$$u_{N-1}=\nabla f^{N-1}(x_{N-1})/\Vert \nabla f^{N-1}(x_{N-1}){\Vert }_2$$

$${\nabla }_{h}f(x)={\left[\frac{f(x+h{e}_i)-f(x)}{h}\right]}_{i=1}^N.$$

$${\nabla }_h{\tilde{f}}^i({\tilde{x}}_i)={\left[\frac{{\tilde{f}}^i({\tilde{x}}_i+h{e}_k)-{\tilde{f}}^i({\tilde{x}}_i)}{h}\right]}_{k=1}^N.$$

$$\Vert L-\tilde{L}{\Vert }_{\text{HS}}\le C\hat{C}(1+\hat{C}K{)}^{N-d}\sqrt{N-d}{h}^{s/2(s-1{)}^i},$$

$$K=4|g^{\prime }{|}_{C^s}+2c_3+\max \{\Vert \nabla f(x_i)\Vert :i=N,\dots ,d+1\}$$

$$\hat{C}=2/\min \{\Vert \nabla f^{N-i}(x_i){\Vert }_2:i=N,\dots ,d+1\}$$

$$C^2=h^{2-s}{c}_3+4|g^{\prime }{|}_C^s{h}^{s-1}(2+h{)}^{2(s-1)}+c_3|g^{\prime }{|}_{C^s}(2+h{)}^{s-1},$$

$$\Vert L-\tilde{L}{\Vert }_{\text{HS}}\lesssim (1+\tilde{C}{)}^{N-d}\sqrt{N-d}h,$$

$$f(x_t)=g(\Vert P{x}_t{\Vert }_2^2)=g(\frac{t^2}{\Vert Px{\Vert }_2^2}\Vert Px{\Vert }_2^2)=g(t^2).$$

$$\Vert {\tilde{u}}_N-u_N{\Vert }_2\le \frac{2C\sqrt{N-d}{h}^{s/2}}{\nabla f({\tilde{x}}_N)}=:S_0,$$

$${\nabla }_{h}f(x{)}_i$$

$$=g^{\prime }({\xi }_{i,h})\frac{\Vert P(x+h{e}_i){\Vert }_2^2-\Vert Px{\Vert }_2^2}{h}$$

$$=g^{\prime }({\xi }_{i,h})\frac{{\sum }_{j=d+1}^N⟨x+h{e}_i,u_j{⟩}^2-⟨x,u_j{⟩}^2}{h}$$

Full text

Key words. Ridge Functions, Function Approximation, Big Data, High Dimensions, Active Variables, Active Subspaces, Optimization over Grassmannian Manifolds

1 Introduction

Nowadays we are living in a world where the acquisition, analysis and storage of big data play a major role. Digital communication, medical imaging, seismology and cosmology are only a few examples, which show the necessity to handle massive data sets. Usually data is modeled as functions $f:X\rightarrow Y$, where $X$ can be $\mathbb{R}^{N}$ or a general curved surface. In particular the approximation of such functions from point queries, when $N$ is very large, is an important field. Such problems arise, for example, in learning theory [24], in modeling physical and biological systems [17], as well as neural networks [3] and in parametric and stochastic PDEs [5].

Because of the so-called curse of dimensionality, a notion introduced in 1961 by Richard Bellman [1], the handling of functions in many variables is an ambitious task. Namely, functions on $\mathbb{R}^{N}$ with smoothness of order $s$ can in general be recovered with an accuracy of at most $n^{-s/N}$, applying $n$-dimensional spaces of linear or nonlinear approximation. Thus, the learning of functions depending on a large number of variables is particularly difficult even with smoothness assumptions on $f$ [8, 9, 20]. Certainly, we need to impose additional structure on $f$ to achieve efficient learning [23, 22, 4, 14, 16].

1.1 Ridge Functions

One popular approach to break the curse of dimensionality is to consider ridge functions of the form

[TABLE]

where $A\in\mathbb{R}^{m\times N}$, with $m$ considerably smaller than $N$, is usually called ridge matrix and $g\in C^{s}(\mathbb{R}^{m})$, $1\!\leq\!s\leq\!2$, is called the ridge profile. The requirement for the function to have at least one derivative is essential. In fact, it was shown in [18] that ridge functions need to have a first derivative uniformly bounded away from zero in the origin in order to reduce the complexity of the approximation task.

For particular choices of $A$ different approaches to successfully learn ridge functions have been investigated. For example, if $A$ is of the form $A^{T}=[e_{i_{1}},\dots,e_{i_{m}}]$, for $e_{i_{k}}\in\mathbb{R}^{N}$ being the canonical unit vectors and $i_{k}\in\{1,\dots,N\}$, $f$ can be rewritten as a function which depends only on a few variables, i.e., $f(x_{1},\dots,x_{N})=g(x_{i_{1}},\dots,x_{i_{m}})$. An approach to recover the active variables and approximating the ridge profile $g$ has been given in [10]. It was shown that by adaptive sampling we can obtain similar estimates as if the active coordinates $i_{1},\dots,i_{m}$ are known to us.

Another special case of (1) is to assume that $m=1$ and that the matrix $A$ is, therefore, a vector, usually called ridge vector and denoted by $A^{T}=:a$. In this case, $f$ is of the form

[TABLE]

The recovery of such ridge functions from point queries was first considered by Cohen, Daubechies, DeVore, Kerkyacharian, and Picard in [4] for ridge functions with a positive ridge vector. It was shown that the accuracy of their method is close to the approximation rate of one-dimensional continuous functions.

However, the algorithm from [4] does not apply to arbitrary ridge vectors. In [23, 14, 16] new algorithms were introduced to waive the assumption of a positive ridge vector. The main idea of the algorithm in [16] is to approximate the gradient of $f$ by divided differences, exploiting the fact that the gradient of $f$ is some scalar multiple of the ridge vector. The accuracy of the approximation of the gradient is determined by the choice of the step size in the computation of the divided differences, whereas the number of sampling points is fixed.

The approach by Fornasier, Schnass and Vybiral [14] is rather based on compressed sensing and applies to (1) very generally. Thus, not the gradient but the directional derivatives of $f$ were approximated at a certain number of random points in random directions. However, especially for the methods in [14, 4], the authors need a restrictive assumption to use compressed sensing techniques. That is, the ridge vector $a$ can be well-approximated by a sparse subset of its coefficients. In [23] this assumption could be removed by leveraging the Dantzig selector [2] to recover an approximation of $a$.

However, the structure assumption on $f$ to be a ridge function can be very restrictive. If we, for example, consider a sensor network, where we have a certain number of sensors, say $N$, which measure the moisture, temperature and pressure to forecast forest fire, the aim is to compute the risk of fire by a function $f:\mathbb{R}^{3N}\rightarrow\mathbb{R}$ depending on the measurements of the sensors. It is then very unlikely that the combination of measurements which yield the same risk of fire lie on a $3N-1$-dimensional hyperplane, since also parameters like topography and vegetation influence the prediction. Much more likely is the assumption that these combinations lie on a lower dimensional manifold.

1.2 Sleeve Functions

To allow for the recovery of more general functions, which are constant along some lower-dimensional submanifolds, we will introduce the notion of sleeve functions and as a special case of linear-sleeve functions. Within this paper we will then investigate and analyze algorithms to capture linear-sleeve functions and we will propose a technique to apply these methods to general sleeve functions.

Definition 1.1

Let $g\in C^{s}[-r,r]$, $r\in\mathbb{R}_{+}$, $M\subset\mathbb{R}^{N}$ a $d$-dimensional smooth submanifold of $\mathbb{R}^{N}$, and $\text{tub}_{r}(M):=\{x\in\mathbb{R}^{N}:\operatorname{dist}(x,M)<r\}$, then we call $f:\text{tub}_{r}(M)\rightarrow\mathbb{R}$ a sleeve function if we can rewrite $f$ in terms of $g$ by

[TABLE]

for $x\in\text{tub}_{r}(M)$. In the case where $M$ is a linear subspace, we call $f$ a linear-sleeve function and denote $L:=M$, to emphasize the special case, i.e., we write:

[TABLE]

The need to restrict $g$ to a bounded domain is twofold; on the one hand, if we wish to recover $g$ from a finite number of sampling points, we do need this restriction and on the other hand, it is useful for the approximation task to have a unique mapping $x\mapsto x_{0}$ with $\operatorname{dist}(x,M)=\operatorname{dist}(x,x_{0})$. Thus, in the case of $M$ being a linear subspace, $r$ can be chosen arbitrarily, where in the general case $r$ is chosen to be the radius of a non-self-intersecting tube around M. For an illustration of linear-sleeve functions we refer to Figure 2.

Note that the notion of sleeve functions is indeed a generalization of ridge function, thus, if $L$ is an $N-1$-dimensional subspace, we can rewrite $\operatorname{dist}(x,L)^{2}=\langle x,a\rangle^{2}$, where $a$ is the normal vector of $L$. Also note that this formulation is very different from the one introduced in [14]. Indeed, if $f$ is of the form $f=g(A\cdot)$, the level sets are linear subspaces, whereas this is not true for linear-sleeve functions (cf. Figure 2).

Furthermore, observe that even by separating the approximation task in approximating $g$ and $M$, we cannot simply use manifold learning algorithms to approximate $M$, since manifold learning algorithms (cf. e.g. [25, 7]) usually assume that we can sample from the manifold. However, we need to reconstruct the level sets (or at least one, namely $M$) without knowing in advance to which level set the sampling points belong; actually it is very likely that all sampling points belong to different level sets.

1.3 Our Contribution

Our work studies the approximation of linear-sleeve function of the form $f(x)=g(\operatorname{dist}(x,L))$ for $x\in\text{tub}_{1}(L)$, where $g\in C^{s}[0,1]$ and $L\subset R^{N}$ is a $d$-dimensional subspace of $\mathbb{R}^{N}$. We will provide and analyze two different algorithms to capture linear sleeve functions from point queries. Our main contributions can be summarized as follows.

• •

Adaptive Algorithm. The first algorithm, to which we refer to as ATPE, is based on the fact that the gradient of $f$ in some $x\in\mathbb{R}^{N}$ is perpendicular to the level set of $f$ in $x$. We will show that the restriction of $f$ to the plane, which is perpendicular to the gradient and which is the tangent plane of the corresponding level set, is again a linear-sleeve function. We will then argue that applying the same fact iteratively to the restrictions of $f$, the tangent plane computed in the $N-d$-th step gives a reconstruction of $L$. In ATPE will then substitute the gradient by divided differences, because we cannot compute the gradient by point queries of $f$.

• •

Optimization Algorithm The second algorithm, to which we refer to as OGM, is based on a minimization over the Grassmannian manifold. Namely, it will define an objective function, whose minimizer is $L$. However, we will see that we cannot define this objective functions using only point samples of $f$. In OGM we therefore approximate this objective function by an objective function whose minimizer $\tilde{L}$ will be proven to be close to $L$.

• •

Error Bounds. Those two algorithms are of a rather different nature. Whereas the approximation success of the first algorithm depends only on the error of the gradient approximation by divided differences, the success of the second algorithms depends on the error of the approximation of the objective function. The first main theorem states that the error of the approximation of the $d$-dimensional subspace $L$ using ATPE can be bounded by

[TABLE]

where $\tilde{L}$ is the approximation of $L$, $h$ can be chosen arbitrarily small but fixed, $C,K$ are some positive constants and the number of function evaluations is given by $(N+1)(N-d)$. For the second main theorem we prove that using OGM the approximation error is given by

[TABLE]

where $M$ is the number of function samples and $\tilde{C}$ a constant only depending polynomial on the dimension of the space. Note that OGM, differently to ATPE, yields a reconstruction error which decreases with the number of sampling points and is not constrained by a fixed number of sampling points and is, therefore, advantageous. However, the first algorithm is more promising to apply also to the manifold case.

• •

Impact on the Approximation of General Sleeve Functions. The next step would be to find algorithms to recover general sleeve functions of the type (3). Due to the fact, that we would need to optimize over all possible $d$-dimensional submanifolds, to approximate general sleeve functions in a similar way as proposed by OGM, we anticipate that an adaptation of ATPE is more promising.

We believe that one can also use gradient approximations to capture general sleeve functions of the type (3). Roughly said, we propose to use the gradients to compute samples from the manifold. More precisely, knowing the gradient of $f$ at some point $x$, again would enable us to approximate the sleeve profile $g$ and, under additional assumption, we could use the direction given by the gradient and the value of $f$ in $x$ to translate $x$ to the manifold. A careful estimation of the distribution of the translated sample points should then enable us to apply manifold learning algorithms (e.g [25]) to estimate the manifold $M$.

The paper is organized as follows: After introducing some preliminaries, we will present and analyze ATPE in Section 3. In Section 4 we will introduce and analyze OGM. The consideration will be completed in Section 5 by some promising numerical results.

2 Preliminaries

To put our results in a precise setting, we introduce the class $\mathcal{LR}(s)$ of all linear-sleeve functions $f(x):=g(\operatorname{dist}(x,L)^{2})$, where $g\in C^{s}[0,1]$ and $L\subset\mathbb{R}^{N}$ a subspace. We use the following norm, subsequently referred to as Hölder norm, on $C^{s}$. For $k<s\leq k+1$, with $k\in\mathbb{N}$, we define

[TABLE]

where $g^{(j)}$ denotes the $j$-th derivative of $g$, and, for $0<\beta\leq 1$, we set

[TABLE]

Note that we call $g$ Lipschitz continuous if $|g|_{\text{Lip}(1)}$ is bounded. We then call $|g|_{\text{Lip}(1)}$ Lipschitz constant or Lipschitz norm of $g$. If we want to highlight the dimension of the vector space, we sometimes write $\|\cdot\|_{\ell_{p}^{N}}:=\|\cdot\|_{p}$ for the $\ell_{p}$ norm of a vector, for $p=1,2$. The weak $\ell_{p}$ norm of a vector $x\in\mathbb{R}^{N}$ is the smallest constant $M$, such that

[TABLE]

We further recall the following useful property of any norm on $\mathbb{R}^{N}$.

Lemma 2.1** **([16])

Let $\|\cdot\|$ be any norm on $\mathbb{R}^{N}$ and $x\in\mathbb{R}^{N}$ with $\|x\|=1$, $\tilde{x}\in\mathbb{R}^{N}\setminus\{0\}$ and $\lambda\in\mathbb{R}$. Then

[TABLE]

We will denote the $i$-th canonical unit vector, with a one in the $i$-th coordinate and zero elsewhere, by $e_{i}$. The Grassmannian manifold of all $d$-dimensional subspaces of $\mathbb{R}^{N}$ is denoted by $G(d,N)$ and for an orthogonal projection $P:\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$, the operator norm is given by the Hilbert-Schmidt norm

[TABLE]

The orthogonal complement of a subspace $P$ is denoted by $P^{\perp}$ and the distance of a vector $x\in\mathbb{R}^{N}$ to a subspace, respective subset, $P\subset\mathbb{R}^{N}$ is defined by

[TABLE]

In the sequel, for two quantities $A,B\in\mathbb{R}$, which may depend on several parameters, we shall write $A\lesssim B$, if there exists a constant $C>0$ such that $A\leq CB$, uniformly in the parameters. If the converse inequality holds true, we write $A\gtrsim B$ and if both inequalities hold, we shall write $A\asymp B$.

Finally, we want to recall some approximation properties of functions in $C[0,1]$. For two given integers $S>1$ and $M\geq 2$, we consider the space $\mathcal{S}_{h,S}$, $h:=1/M$, of piecewise polynomials of degree $S-1$ with equally spaced knots at the points $ih$, $i=1,\dots,M-1$, and having continuous derivatives of order $S-2$. It is well-known (cf. e.g. [8]) that there is a class of linear operators $Q_{h}$ which maps $C[0,1]$ into $\mathcal{S}_{h,S}$. These operators are usually called quasi-interpolants. For a function $g\in C[0,1]$, the application of a quasi-interpolant only depends on the values of $g$ at the points $ih$, $i=0,\dots,M$. Furthermore, we can choose the operator $Q_{h}$ to fulfill the following property: For all $g\in C^{s}([0,1])$

[TABLE]

with $C$, a constant depending only on $S$ [8].

3 An Adaptive Algorithm Estimating Tangent Planes

An obvious approach to approximate sleeve functions of the form (3) is to apply the methods to recover classical ridge functions of the form (2), i.e., to approximate a linear sleeve function by a classical ridge function. Of course, this method can only provide good approximation results if the sleeve function is close to a classical ridge function in a certain sense, cf. [15]. However, it would be more convenient to approximate a function of the form (3) by an estimator of the same form.

We first observe that we can rewrite a linear-sleeve function as

[TABLE]

where $P_{P}$ is the orthogonal projection to the $(N-d)$-dimensional subspace $P$ orthogonal to $L$. For simplicity we will denote the orthogonal projection of a vector $x\in\mathbb{R}^{N}$ to a subspace $H=\operatorname{\operatorname{span}}\{u_{1},\dots,u_{d}\}\in G(d,N)$ by $Hx:=P_{H}(x)$. This notation relates to the matrix representation of an orthogonal projection given by $H=\sum_{i=1}^{d}u_{i}u_{i}^{T}$.

The algorithm, we will introduce in this section, will, similarly as in [4, 16], exploit the fact that we can estimate the tangent plane in some $x_{0}\in\mathbb{R}^{N}$ of the $(N-1)$-dimensional submanifold

[TABLE]

as the unique hyperplane which is perpendicular to the gradient of $f$ in $x_{0}$. We will then show that the function $f$ restricted to this tangent plane is again of the form (4). Of course, we cannot compute the gradient by sampling the function; in the subsequent proposed algorithm we, therefore, approximate the gradient by divided differences.

3.1 The Algorithm

As mentioned before the idea of the first algorithm is to use the fact that the gradient of the linear-sleeve function $f$ in some $x_{0}$ is perpendicular to the level sets and that the restriction of $f$ to the corresponding tangent plane is again a linear sleeve function (cf. Figure 3). We will further show that applying this fact iteratively to restrictions of $f$ finds after $N-d$ steps the wanted subspace $L$. To prove this statement we introduce an adaptive algorithm to exactly recover the subspace $L$ by computing gradients of $f$ (see Algorithm 1) and show in Theorem 3.1 that this algorithm can recover the subspace $L$ exactly. Thus, our first aim is to show that the system, which is formed by the gradients of the restrictions of $f$, forms a basis for $P=L^{\perp}$. Observe that an essential step in this algorithm is to compute gradients of $f$ and is therefore not useful to approximate $f$ from point samples. It only serves as an auxiliary tool to introduce the first main Algorithm 2.

As mentioned before this algorithm iteratively computes restrictions of $f$, such that after $N-d$ steps the restriction of $f$ will be exactly defined on $L$. Assume we have computed the tangent plane $T^{i}$ and the restriction $f^{i}$ of $f$ the subspace $T^{i}\subset\mathbb{R}^{N}$, for $i=N,\dots,d+1$. The algorithm then chooses uniformly at random a point $x_{i}\in T^{i}$ and computes the gradient of $f^{i}$ in $x_{i}$ (1.). It then normalizes this gradient (2.) and determines the subspace $T^{i-1}$ which is orthogonal to this gradient in $T^{i}$ (3.). Finally we restrict $f^{i}$ to $T^{i-1}$ and repeat this procedure until $i=d+1$. The following theorem states that ATPC indeed succeeds to recover $L$.

Theorem 3.1

Let $f(x)=g(\operatorname{dist}(x,L)^{2})=g(\|Px\|_{2}^{2})$, $g\in C^{s}[0,1]$, $1<s\leq 2$, for $x\in\text{tub}_{1}(L)$ and the linear subspace $P^{\perp}=L\subset\mathbb{R}^{N}$. Compute $T^{d}$ as proposed in Algorithm 1. Then $L$ coincides with $T^{d}$ almost surely (e.g. due to the Lebesque measure on $\mathbb{R}^{N}$).

Proof.

We write $L$ as $L=\operatorname{\operatorname{span}}\left\{u_{1},\dots,u_{d}\right\}$ where $\left\{u_{1},\dots,u_{d},u_{d+1},\dots,u_{N}\right\}$ is an orthonormal basis of $\mathbb{R}^{N}$, and let $V=\left[u_{1}\dots u_{d}\,u_{d+1}\dots u_{N}\right]$ be the corresponding matrix. We begin by computing the gradient of $f$ and obtain

[TABLE]

which is obviously perpendicular to $L$ if $Px_{N}\neq 0$, which is almost surely true.

Due to the orthogonality of $\nabla{f}(x_{N})$ to $L$, we can assume that $u_{N}\!=\!\nabla{f}(x_{N})/\|\nabla{f}(x_{N})\|_{2}$ and then by definition $T^{N-1}=(\operatorname{\operatorname{span}}\{u_{N}\})^{\perp}=\operatorname{\operatorname{span}}\{u_{1},\dots,u_{N-1}\}$. We now define $f^{N-1}$ to be the restriction of $f$ to $T^{N-1}$ and $h^{N-1}:\mathbb{R}^{N-1}\rightarrow\mathbb{R}$ by

[TABLE]

for $\hat{x}\in\mathbb{R}^{N-1}$. Then $\nabla{h^{N-1}}^{T}(\hat{x})=\nabla{f}(x)^{T}\hat{V}_{N-1}$, where $\hat{V}_{i}:=\left[u_{1}\dots u_{i}\right]$, $i=1,\dots,N$, and $x=:V(\hat{x},0)^{T}\in T^{N-1}$. Thus, the gradient of $h^{N-1}$, considered as a vector in $\mathbb{R}^{N}$ is given by

[TABLE]

where $V_{i}:=\left[u_{1}\dots u_{i}0\dots 0\right]$. We conclude that for some $x_{N-1}\in T^{N-1}$ chosen uniformly at random, we have

[TABLE]

A straightforward computation shows that $VV_{N-1}^{T}Px$ is the projection of $Px$ to $T^{N-1}$. Thus, it is obvious that for some $x_{N-1}\in T^{N-1}$, chosen as in Algorithm 1, $\nabla{f}^{N-1}(x_{N-1})$ is perpendicular to $L$, i.e., if $Px_{N-1}\neq 0$, which is almost surely the case.

Further, $\nabla{f}^{N-1}(x_{N-1})$ is also orthogonal to $u_{N}$, since it lies in $T^{N-1}$. Therefore, we set

[TABLE]

for some $x_{N-1}\in T^{N-1}$ and $T^{N-2}:=\operatorname{\operatorname{span}}\left\{u_{1},\dots,u_{N-2}\right\}$. Note that again $Px_{N-1}\neq 0$ holds almost surely. We repeat this procedure until we get a basis $\left\{u_{d+1},\dots,u_{N}\right\}$ of $P$, which yields the desired space $L=P^{\perp}$. ∎

The previous theorem shows that, if we could compute the gradient in $(N-d)$ points, we would be able to recover the space $L$, respectively its orthogonal complement $P$, exactly. However, ATPC is not based on sampling the function $f$, since gradients cannot be computed exactly using point queries. Thus, we can only approximate the gradients by computing the divided differences

[TABLE]

We adapt the first step in ATPC by substituting the computation of the gradients by computing divided difference and propose Algorithm 2, to which we will refer as ATPE, for the approximation task.

The described procedure, of course, cannot find the correct plane $P$. However, it is able to compute a good approximation of $P$, where the approximation error depends on the choice of $h$. For reasons of clarity, the proof of the next theorem is moved to the next subsection.

Theorem 3.2

Let $f$ be a linear-sleeve function of the form (4), i.e., $f\in\mathcal{LR}(s)$ for some $s\in(1,2]$. Assume that the derivative of $g\in C^{s}([0,1])$ is bounded by some positive constants $c_{2},c_{3}$. By sampling the function $f$ at $(N-d)(N+1)$ appropriate points, ATPE constructs almost surely an approximation of $L$ by a subspace $\tilde{L}\subset\mathbb{R}^{N}$, such that the error is bounded by

[TABLE]

for some arbitrarily small $h>0$, where

[TABLE]

which are constants depending only on the Hölder norm and bounds of $g^{\prime}$.

In particular if $g^{\prime}$ is Lipschitz continuous, i.e. $s=2$, it holds

[TABLE]

with $\tilde{C}=\hat{C}K$.

It then only remains to recover the ridge profile $g$. The estimation of $g$ is rather straightforward. As the gradient gives the direction in which $f$ changes, f becomes a one-dimensional function in the direction of the gradient. Hence, we can estimate $g$ with well-known numerical methods. Indeed, we have already seen that the gradient of $f$ in some point $x$ is given by $\nabla f(x)=g^{\prime}(\|Px\|_{2}^{2})Px$, i.e., the normalized direction is $a:=Px/\|Px\|_{2}$. Setting $x_{t}:=ta$ yields

[TABLE]

However, similar to the algorithm in [16], this algorithm uses a fixed number of samples and the estimation cannot be improved by taking more samples. We, therefore, also aim for an algorithm which yields a reconstruction whose error decreases with the number of sampling points (cf. Section 4). Also note that due to the adaptive character of ATPE the reconstruction error increases for smaller values of the Lipschitz continuity $s$.

To complete this subsection, we want to remark that we can perform a similar, but slightly worse, error analysis for the case that $f$ is of the form $f(x)=g(\|Px\|_{2})$, whereas before we considered sleeve functions of the form $f(x)=g(\|Px\|^{2}_{2})$. See Remark 1 in the next subsection for a short explanation.

3.2 Proof of Theorem 3.2

As mentioned above, the idea is to approximate the gradients of $f=f^{N}$ and $f^{i}$ for $i=d+1,\dots,N-1$. Since we need $N+1$ samples for each gradient approximation, we need $(N-d)(N+1)$ samples altogether. We already know from Theorem 3.1 that the subspace $P$ can be written in terms of the gradients of the restrictions of $f$ and $f$ itself. Hence, we assume $P=L^{\perp}=\operatorname{\operatorname{span}}\left\{u_{d+1}\dots,u_{N}\right\}$, where the $u_{i}$’s are given as stated in Theorem 3.1. We split the proof by establishing several lemmata.

Lemma 3.3

Under the assumptions of Theorem 3.2 and with the choice of the $\tilde{u}_{i}$’s, $i=d+1,\dots,N$, as proposed in Algorithm 2, we have

[TABLE]

where $\tilde{x}_{N}\in\mathbb{R}^{N}$ is chosen uniformly at random $u_{N}=\frac{\nabla f(\tilde{x}_{N})}{\|\nabla f(\tilde{x}_{N})\|_{2}}$.

Proof.

First, let us estimate the error between $\nabla f(x)$ and $\nabla_{h}f(x)$. We can compute

[TABLE]

for some $\xi_{i,h}$ between $\|Px\|_{2}^{2}$ and $\|P(x+he_{i})\|_{2}^{2}$, where $u_{ji}$ denotes the $i$-th entry of the vector $u_{j}$. We then estimate

[TABLE]

where we used the fact that $u$ is a unit vector and that, therefore, all its entries (in absolute value) and the entries of its projection are smaller than or equal to one. Thus, the error which we obtain by approximating the gradient can be estimated as

[TABLE]

To estimate those terms we take the following inequality into account:

[TABLE]

where we used in the second as well as in the last step that $\{u_{j}\}_{j=1}^{d}$ forms an orthonormal system, which spans $H$, so that $\sum_{j=1}^{d}u_{ji}^{2}\leq 1$ for each $i=1,\dots,N$ and $\sum_{i=1}^{N}u_{ji}^{2}=1$ for each $j=1,\dots,d$. Now the desired estimates follow immediately:

[TABLE]

Thus, we can find a constant $C>0$, independent of the dimensions $d$ and $N$, such that

[TABLE]

For the exact choice of $C$ we refer to Theorem 3.2. Hence, applying Lemma 2.1, with $\lambda=1/\|\nabla_{h}f(x)\|_{2}$, proves the claim. ∎

Next, we use the approximation of the gradient to approximate the tangent plane $T^{N-1}$ at $x$ with $\tilde{T}^{N-1}=\operatorname{\operatorname{span}}\{\nabla_{h}f(x)^{\perp}\}$. The approximation error is then, of course, given by (21). Further, we let $f^{N-1}$ and $\tilde{f}^{N-1}$ be the restriction of $f$ to $T^{N-1}$ and $\tilde{T}^{N-1}$, respectively.

Again we want to compute the column vectors $u_{i}$ of $V$, $i=N,\dots,d+1$, step by step as the normalized gradients of $f$, $f^{i}$. But instead of computing the gradient of $f^{j}$ we can only approximate it through an approximation of the gradient of $\tilde{f}^{j}$. Thus, we iteratively set the columns $\tilde{u}_{i}$, $i=N,\dots,d+1$ of $\tilde{V}$ as the normalized approximated gradients of $\tilde{f}_{i}$. The error of the approximation in each step can then be estimated by means of the following lemmata, in particular, by means of Lemma 3.4 for the first step and Lemma 3.7 for the $i$-th step.

Before stating these lemmata, let us recall the definition of the matrices $V,\tilde{V},V_{i}$ and $\tilde{V}_{i}$, $i=N-1,\dots,d+1$. Let

[TABLE]

for $i=N,\dots,d+1$, where $x_{i}=VV_{i}^{T}\tilde{x}_{i}$ and $\tilde{x}_{i}$ as well as $f^{i}$ and $\tilde{f}^{i}$ are chosen as proposed in the Algorithms ATPC and ATPE. Then $\{u_{d+1},\dots,u_{N}\}$ as well as $\{\tilde{u}_{d+1},\dots,\tilde{u}_{N}\}$ form orthonormal systems and $u_{i},\tilde{u}_{i}$, $i=1,\dots,d$, are chosen that the whole systems $\{u_{1},\dots,u_{N}\}$ and $\{\tilde{u}_{1},\dots,\tilde{u}_{N}\}$ form an orthonormal basis for $\mathbb{R}^{N}$. We can now define

[TABLE]

Lemma 3.4

With the same assumptions and choices as in Theorem 3.2 and Lemma 3.3, we have

[TABLE]

We first have to prove the following lemma:

Lemma 3.5

With the same assumptions and choices as in Theorem 3.2 and Lemma 3.3, let $x:=VV_{N-1}^{T}\tilde{x}$, for some $\tilde{x}\in\tilde{T}^{N-1}$. We then have

[TABLE]

Proof.

We write $\tilde{x}=\sum_{i=1}^{N-1}\tilde{x}_{i}\tilde{u}_{i}$ with $\tilde{x}_{i}\in\mathbb{R}$ and $|\tilde{x}_{i}|\leq 1$. Then we compute

[TABLE]

In the step before the last we used that $\{\tilde{u}_{1},\dots,\tilde{u}_{N}\}$ is an orthonormal basis (according to Step 2. and 3. in the ATPE Algorithm) and in the last step additionally that $1-\langle\tilde{u}_{N},u_{N}\rangle^{2}=(1+\langle\tilde{u}_{N},u_{N}\rangle)(1-\langle\tilde{u}_{N},u_{N}\rangle)\leq 2(1-\langle\tilde{u}_{N},u_{N}\rangle)$. By observing that

[TABLE]

we deduce the claim. ∎

We are now able to prove the error of the first step of our algorithm ATPE, i.e., Lemma 3.4.

Proof of Lemma 3.4.

For simplicity we will write $\tilde{x}$ and $x$ instead of $\tilde{x}_{N-1}$ and $x_{N-1}$. Note that it is straightforward to show (compare to Equation (16)) that $\nabla f^{i}(x)=VV_{i}^{T}\nabla f(x)$ and $\nabla_{h}f^{i}(x)=\tilde{V}\tilde{V}_{i}\nabla_{h}f$. We, hence, obtain that

[TABLE]

where we applied Lemma 2.1 with $\lambda=\|\nabla_{h}\tilde{f}^{N-1}(\tilde{x})\|_{2}$ in the second step. Now, we can estimate

[TABLE]

as well as

[TABLE]

Finally, we estimate $T_{2}$ by

[TABLE]

which proves the lemma. ∎

One can now prove similar estimates as in Lemma 3.4 for $\|u_{i}-\tilde{u}_{i}\|_{2}$, $i=d+1,\dots,N-2$. However, we first need to prove a more general version of Lemma 3.5.

Lemma 3.6

With the same assumptions and choices as in Theorem 3.2 and Lemma 3.3, we have for $\tilde{x}_{i}\in\tilde{T}^{i}$ and $x=VV^{T}_{i}\tilde{x}_{i}$, $i=N-2,\dots,d+1$, that

[TABLE]

This inequality in turn yields the desired generalization of Lemma 3.4:

Lemma 3.7

With the same assumption and choices as in Theorem 3.2 and Lemma 3.3, for $i=0,\dots,d$, we have

[TABLE]

for all $i=0,\dots,N-d-1$.

Proof.

Using the same methods as in the proof of Lemma 3.4 and using Lemma 3.6, one can show

[TABLE]

where the $S_{j}$, $j=1,\dots,d$ are recursively defined. The claim then follows by induction. ∎

Putting the conclusions of the previous lemmata together and observing that

[TABLE]

finishes the proof of Theorem 3.2. ∎

As mentioned in the end of the last subsection, we can perform a similar, but slightly worse, error analysis for the case that $f$ is of the form $f(x)=g(\operatorname{dist}(x,L))=g(\|Px\|_{2})$. Indeed, we can estimate the approximation error of the gradient in the following way:

Remark 1

To approximate a function of the form $f(x)=g(\operatorname{dist}(x,L))=g(\|Px\|_{2})$, for $x\in\mathbb{R}^{N}$, $g\in C^{1}([0,1])$ and $g^{\prime}$ Lipschitz continuous, we can utilize the following observation to obtain a worse approximation result than for linear sleeve functions of the form (4). Namely, we can rewrite the $i$-th entry of the divided difference of $\|Px\|_{2}$ as:

[TABLE]

We then obtain the following estimate:

[TABLE]

And, hence, we have:

[TABLE]

If $\|Px\|_{2}$ is small, this upper bound can become large. Fortunately, in the case of linear sleeve functions of the form (4), we are not constrained by this term (cf. Equation (22)).

4 An Optimization Algorithm on the Grassmannian Manifold

We will now reformulate the given approximation problem as an optimization problem over the Grassmannian manifold $G(d,N)$. This reformulation allows us to develop an algorithm which yields a reconstruction whose error decreases with the number of sampling points. Remember that the previous algorithm needed a fixed number of sampling points and the error has decreased with the step size $h$ in the computation of the divided differences. We again use the following notation for $f$

[TABLE]

where the operator $P$ denotes the orthogonal projection $P_{P}$ to the subspace $P\subset\mathbb{R}^{N}$ orthogonal to $L$. In the sequel of this section, we will, for brevity, assume that $\dim P=d$ where we assumed $\dim P=N-d=N-\dim L$ in the last section.

4.1 The Algorithm

Let us assume without loss of generality that $g$ is not the constant function. Otherwise we do not need to find the subspace, since every subspace can be used to represent $g$ as stated. We first define for each $H\in G(d,N)$ a function $f_{H}$ as linear-sleeve function, namely by

[TABLE]

The main idea of the algorithm then uses the fact that $f_{P}=f$ and $f_{H}\neq f$ for $H\neq P$ and, thus, that $P$ is the unique minimizer of

[TABLE]

Unfortunately, we cannot express this objective function in terms of sampling the input function $f$. On the one hand, we, therefore, need to replace the integral by a finite sum and on the other hand, the definition of $f_{H}$ is not clear, since we do not know $g$ in advance.

However, note that we can easily recover $g$ by sampling $f$ in some random direction $\hat{\theta}$. Indeed, it holds for $t\in\mathbb{R}$ that $f(t\hat{\theta})=g(t\|P\hat{\theta}\|_{2}^{2})$ and, since $\hat{\theta}$ is almost surely not contained in the orthogonal complement of $P$, $g$ is up to the constant $\|P\hat{\theta}\|_{2}^{2}$ uniquely determined by $f(\cdot\hat{\theta})$. Hence, if we knew $\|P\hat{\theta}\|_{2}^{2}$ approximately, sampling $f$ at $ih\hat{\theta}$, $i=1,\dots,M$, where $h\in(0,1)$ is the step size, gave an approximation to $g$. Namely, with $g_{\hat{\theta}}=g(\cdot\|P\hat{\theta}\|_{2}^{2})$, we let $\hat{g}_{\hat{\theta}}^{M}$ be the approximation of $g_{\hat{\theta}}$ from the sampling points $(ih,f(ih\hat{\theta}))$. We can, then, set $\hat{f}_{H}^{M}(x)=\hat{g}^{M}_{\theta}({\hat{p}_{\hat{\theta}}}^{-1}{\|Hx\|_{2}^{2}})$, where $\hat{p}_{\hat{\theta}}$ is an approximation of $\|P\hat{\theta}\|_{2}^{2}$.

One possibility to choose $\hat{\theta}$ such that we know $\|P\hat{\theta}\|_{2}^{2}$ approximately, is to choose $\hat{\theta}$ as the approximation of the normalized gradient of $f$ in some random direction $\eta$. Indeed, the normalized gradient of $f$ is given by

[TABLE]

Therefore, we have almost surely (e.g., with respect to the Lebesque measure) $\|P\theta\|_{2}=1$. Thus, we choose

[TABLE]

and let $\hat{g}_{\hat{\theta}}^{M}=Q_{h}g_{\hat{\theta}}$ be the approximation of $g_{\hat{\theta}}$, with $Q_{h}$ as introduced in the preliminaries (see Equation (13)).

The only remaining task is now to substitute the integral by a finite sum. Hence, we aim to define the objective function

[TABLE]

for some $x_{1},\dots,x_{n}\in\mathbb{R}^{N}$ such that $P$ is the unique minimizer of

[TABLE]

to ensure that the minimizer of $\hat{F}^{M}$ is a good approximation of $P$. Certainly, $P$ can only be the unique minimizer of $F$ if it uniquely minimizes the function

[TABLE]

It is, therefore, necessary to find $x_{1},\dots,x_{n}$ such that the map

[TABLE]

is injective. This problem is known as projection retrieval [12] and discussed in the next subsection.

The proposed procedure is summarized in Algorithm 3, to which we refer to as OGM.

We will now be able to prove the following main result in Subsection 4.3.

Theorem 4.1

Let $f$ be a linear-sleeve function of the form (4), i.e., $f\in\mathcal{LR}(s)$, $s\in(1,2]$, and $f=g(\operatorname{dist}(x,L)^{2})=g(\|Px\|_{2}^{2})$ for some $d$-dimensional subspace $P\subset\mathbb{R}^{N}$. Suppose that the derivative of $g$ is bounded from below by some positive constant, and let $\hat{P}:=\operatorname*{argmin}_{H\in{G(N-d,N)}}\hat{\mathcal{F}}^{M}(H)$. Then, we have

[TABLE]

almost surely (e.g., with respect to Haar measure), with a constant depending only polynomial on the dimension of the space. In particular, if $f\in\mathcal{LR}(2)$, then

[TABLE]

Note that the statement holds indeed only almost surely, since we have to ensure that $P\eta\neq 0$.

The first step is to find measurements $\{x_{i}\}_{i=1}^{n}$ which ensure that the map defined in (57) is injective. Then we can turn to the error analysis and the proof of Theorem 4.1.

4.2 Projection Retrieval

To find sampling points which ensure that the objective function has a unique minimizer, we consider the necessary case where $g$ is the identity and $\mathcal{F}$ is, therefore, given by $\mathcal{F}(H)=\sum_{i=1}^{n}\left(\|Px_{i}\|_{2}-\|Hx_{i}\|_{2}\right)^{2}$. Thus, $P$ is the unique minimizer of $\mathcal{F}(H)$ if and only if the sampling points $x_{i}$, $i=1,\dots,n$, determine $P$ uniquely, i.e., if the map

[TABLE]

is injective. We can show the following theorem.

Theorem 4.2

For every $P\in G(d,N)$ the quantities

[TABLE]

uniquely determine $P$.

Proof.

Let $P,H\subset\mathbb{R}^{N}$ be two subspaces and let $\{u_{1},\dots,u_{d}\}$ be an orthonormal basis for $P$ and $\{v_{1},\dots,v_{d}\}$ an orthonormal basis for $H$. Further, suppose that $\|Px_{i,k}\|_{2}^{2}=\|Hx_{i,k}\|_{2}^{2}$ for $i=1,\dots,N$, $k=i,\dots,N$. For $i=k$, we obtain

[TABLE]

This shows that the entries of the diagonal of the projection matrices corresponding to $P$ and $H$ coincide. For the case $i\neq k$, we can compute

[TABLE]

Therefore, using the knowledge from the case where $i$ equals $k$, this equation gives

[TABLE]

Thus, due to the symmetry of a real projection matrix, both projection matrices coincide, since the left-hand side equals the $(i,k)$-th entry of the projection matrix corresponding to $H$ and the right-hand side to $P$, respectively. Therefore, we can conclude that $H=P$. ∎

We see that we need $(N^{2}+N)/2$ sampling points to ensure injectivity, if we choose them as suggested by the last theorem. However, we believe a smaller number of sample points should be sufficient. Indeed, we can ensure that a fewer number of sampling points are sufficient to ensure almost surely injectivity and, therefore, almost surely a unique minimizer. For this, we adapt Theorem 4 in [12] to the real case and deduce that to ensure almost injectivity, we can require only $(d+1)(N+d/2)$ points.

Theorem 4.3** **([12])

Draw a random subspace $P$ uniformly with respect to the Haar measure from the Grassmannian manifold $G(d,N)$. Then the quantities

[TABLE]

for $i\in\{1,\dots,N\}$ and $k\in\{i+1,\dots,d\}$, uniquely determine $P$ a with probability of $1$.

Note that the change in the second index set in comparison to [12] is due to taking the symmetry of a real projection matrix into account. And the change in the number of necessary measurements is due to some small typo in [12], since for this proposed procedure we need to compute not only the first $d$ columns of the projection matrix, but also all its diagonal entries. However, also in [12] it is proven that the first $d$ columns of the projection matrix $P$ determine the corresponding subspace almost surely uniquely. Thus, it would be desirable to find points which allow us to directly determine the entries of the first $d$ columns, without computing all diagonal entries. This would deduce the number of necessary measurements to $Nd$. In the case $d=1$, the following result by Fickus, Mixon, Nelson, Wang [13], tells us that $N+1$ measurements are sufficient to ensure almost injectivity, which is almost the conjectured number $Nd$.

Theorem 4.4** **([13])

Consider $\Phi=\left\{\phi_{i}\right\}_{i=1}^{n}\subset\mathbb{R}^{N}$ and the intensity measurement mapping $A:\mathbb{R}^{N}/{\pm 1}\rightarrow\mathbb{R}^{n}$ defined by $(A(x))(i):=|\langle x,\phi_{i}\rangle|^{2}$ . Suppose each $\phi_{i}$ is nonzero. Then $A$ is almost surely injective if and only if $\Phi$ spans $\mathbb{R}^{N}$ and $\text{rank }\Phi_{S}+\text{rank }\Phi_{S^{C}}>N$ for each nonempty proper subset $S\subset\{1,\dots,N\}$.

This shows that $A$ cannot be almost injective if $n<N+1$. Moreover, for the case $n=N+1$, it is almost injective if and only if $\Phi$ is full spark, which means that every size-$N$ collection of vectors of $\Phi$ is linearly independent. We remark that this result does not stand in contradiction to the above mentioned conjecture that $Nd$ measurements are sufficient. Indeed, in our case we want to recover the subspace and in this sense we can interpret the condition $1=\|v\|_{2}=\langle v,v\rangle$, for some vector $v\in\mathbb{R}^{N}$ of a orthonormal basis of this vector, as an additional measurement. Of course, it is not generally true that every size-$N$ subcollection of $\left\{e_{1},\dots,e_{N},v\right\}$ forms a spanning set, because $v$ could have zero entries. However, if we consider $v$ as embedded in $\mathbb{R}^{m}$, where $\{i_{1},\dots,i_{m}\}$ are the indices of the nonzero-entries of $v$, then every size-$m$ subcollection of $\left\{e_{i_{1}},\dots,e_{i_{m}},v\right\}$ forms a spanning set for $\mathbb{R}^{m}$. Thus by Theorem 12 in [13], the entries $i_{1},\dots,i_{m}$ of $v$ are uniquely determined by these sampling points with a probability of $1$. That the other entries are equal to zero is already determined by the other sampling points.

As indicated by this discussion, we suspect that $Nd$ sampling points are sufficient to ensure almost surely injectivity of $A$. And, indeed, we are able to prove the following theorem, which states that even a fewer number of sampling points are sufficient, although, we do not directly determine the first $d$ columns of the projection matrix.

Theorem 4.5

Draw a random subspace $P$ uniformly with respect to the Haar measure from the Grassmannian manifold $G(d,N)$. Then the quantities

[TABLE]

for a randomly chosen vector $x\in\mathbb{R}^{N}$, $i\in\{1,\dots,N-1\}$, $j\in\{1,\dots,d\}$ and $k\in\{j+1,\dots,N\}$, uniquely determine $P$ with a probability of $1$.

Proof.

We start with the same argumentation as in [12], which tells us that the first $d$ columns of the projection matrix $P$ are linearly independent for almost every $P$, that the diagonal entries of the projection matrix $P$ are given by

[TABLE]

and that the other entries can be computed as

[TABLE]

We further note that we only need to observe $N-1$ diagonal entries to determine all $N$ diagonal entries of $P$. Indeed, if $u_{1},\dots,u_{d}$ is any orthonormal system which spans $P$, it holds that

[TABLE]

and, thus, $P_{NN}=d-\sum_{i=1}^{N-1}P_{ii}$. Hence, by observing $\|Pe_{i}\|_{2}^{2}$, $i=1,\dots,N-1$, as well as $\|P(e_{j}+e_{k})\|_{2}^{2}$, $j=1,\dots,d$ and $k=j+1,\dots,N$, we can recover the first $d-1$ columns of the projection matrix $P$.

We now claim that there exist only finitely many projection matrices with the same first $d-1$ columns. This would show that there are only finitely many subspaces $H$ which yield the same measurements as $P$ for the stated collection of quantities and that we, therefore, can almost surely uniquely recover $P$ by an additional random measurement.

Thus, suppose we already know the first $d-1$ columns of the projection matrix. These are clearly linearly independent in the same event, where the first $d$ columns of $P$ are linearly independent. We can now use the fact that for a projection matrix $P$, it has to hold that $PP=P$ and, hence, that each column of $P$ lies in the span of the corresponding subspace. Applying Gram-Schmidt orthonormalization to the first $d-1$ columns, which are linearly independent, gives an orthonormal system of $d-1$ vectors, which we denote by $u_{1},\dots,u_{d-1}$. Thus, there is only one unknown basis vector, denoted by $u_{d}$, for the subspace $P$ left. However, the measurements we have already taken determine the entries of this vector uniquely in absolute value. Indeed, we have

[TABLE]

which is equivalent to

[TABLE]

Note that we already know the right-hand side of the last equation. This shows that there are indeed only finitely many subspaces which produce the same measurements as $P$. Hence, taking some random measurement in addition, yields the desired almost injectivity.∎

The following corollary shows that if the dimension of the subspace is $d>N/2$, we can choose the same measurements as if the dimension would be $N-d$. So we can further deduce the number of measurements. For example, if the dimension $d=N-1$, we can apply Theorem 4.4 and find that $N$ measurements are sufficient to ensure almost injectivity of $A$.

Corollary 4.6

With the same choice of measurements as in Theorem 4.5 we can uniquely determine a randomly drawn subspace $P\in G(N-d,N)$ with a probability of $1$, i.e., the measurements

[TABLE]

for some random vector $x\in\mathbb{R}^{N}$, $i\in\{1,\dots,N-1\}$, $j\in\{1,\dots,d\}$ and $k\in\{j+1,\dots,N\}$, uniquely determine $P$ with probability $1$.

Proof.

For every $y\in\mathbb{R}^{N}$ it holds $\|Py\|_{2}^{2}=\|Hy\|_{2}^{2}$ if and only if $\|P^{\perp}y\|_{2}^{2}=\|y\|_{2}^{2}-\|Py\|_{2}^{2}=\|y\|_{2}^{2}-\|Hy\|_{2}^{2}=\|H^{\perp}y\|_{2}^{2}$, and therefore, we can apply the results of the above theorem. ∎

4.3 Proof of Theorem 4.1

In the last subsection we have shown that we need less than $Nd$ sampling points to ensure that the objective function defined in (57), where $g$ is assumed to be the identity, has almost surely a unique minimizer. However, for ease of computation we will use the sampling points proposed in Theorem 4.2. We first show that the measurements given in Theorem 4.2 also ensure a unique minimizer of $F_{H}$. For this purpose we introduce a bijective mapping

[TABLE]

and set $n=N(N+1)/2$ as well as

[TABLE]

Lemma 4.7

Suppose that $g$ fulfills the assumption of Theorem 4.1. Then $P$ is the unique minimizer of

[TABLE]

where $x_{i}$, $i=1,\dots,n$, are defined as above.

Proof.

Suppose $H\in G(d,N)$ is another minimizer. Then for all $i\in\left\{1,\dots,n\right\}$, we conclude $g(\|Px_{i}\|^{2}_{2})=g(\|Hx_{i}\|^{2}_{2})$. But since $g$ is injective, this implies $\|Px_{i}\|_{2}=\|Hx_{i}\|_{2}$ for all $i\in\left\{1,\dots,n\right\}$. Thus by the statements proved in subsection 4.2 we conclude that $P=H$. ∎

Lemma 4.8

Under the same assumption as in Theorem 4.1, we have

[TABLE]

where $\hat{\mathcal{F}}^{M}$ was introduced in Equation (55).

Proof.

We start by estimating $\hat{\mathcal{F}}^{M}(H)$ as

[TABLE]

where we used the triangle inequality for $\|\cdot\|_{\ell_{2}^{n}}$ in the second step. We can apply a similar argument to $\mathcal{F}(H)$ to derive $\mathcal{F}(H)\leq\hat{\mathcal{F}}^{M}(H)+\sqrt{\sum_{i=1}^{n}(f_{H}(x_{i})-\hat{f}^{M}_{H}(x_{i}))^{2}}$. This in turn yields the inequality

[TABLE]

We now split this inequality by

[TABLE]

For $T_{1}$, we estimate

[TABLE]

with a constant depending on $\|g^{\prime}\|_{\infty}$ and $\|\nabla f(\eta)\|_{2}$. Note that we used the estimate from Lemma 3.3.

Using Property (13), we can bound $T_{2}$ by

[TABLE]

where $C$ is a constant depending only on the degree of the interpolating polynomials. This proves the lemma. ∎

Theorem 4.9** **(Convergence)

Under the assumptions of Theorem 4.1, suppose that $\hat{P}$ is a minimizer of $\hat{\mathcal{F}}^{M}$. Then

[TABLE]

with a constant depending on the Hölder norm and bounds of $g^{\prime}$.

Proof.

Let $H_{0}$ be a minimizer of $\hat{\mathcal{F}}^{M}$. First note that

[TABLE]

because $\hat{\mathcal{F}}^{M}(H_{0})\leq\hat{F}^{M}(P)=\hat{\mathcal{F}}^{M}(P)-\mathcal{F}(P)\leq C\sqrt{d}h^{s/2}$, where we used the fact that $H$ is a minimizer in the first step, that $\mathcal{F}(P)=0$ in the second step and the statement of the last lemma in the third step. The stated bound then follows from $\mathcal{F}(H_{0})\leq\left|\mathcal{F}(H_{0})-\hat{\mathcal{F}}^{M}(H_{0})\right|+\left|\hat{\mathcal{F}}^{M}(H_{0})\right|\leq 2C\sqrt{d}h^{s/2}$. This yields

[TABLE]

Now define the matrix $\tilde{P}$ by

[TABLE]

and $\tilde{H}$ analogously by

[TABLE]

We further denote the matrix which only contains the diagonal entries of a matrix $P$, and is zero elsewhere, by $P_{D}$ and the matrix which is zero along the diagonal and coincides at all off-diagonal entries with $P$ by $P_{OD}$. We then have

[TABLE]

Note that the $\{ij\}$-th entry $P_{ij}$ of the projection matrix $P$ to the belonging subspace is given by

[TABLE]

Thus, defining

[TABLE]

gives

[TABLE]

Applying (85), we can, therefore, deduce

[TABLE]

which proves the claim. ∎

5 Numerical Results

In this section we investigate the performance of the approximation schemes presented in the last sections. Our algorithms separate the approximation task in approximating the one-dimensional function $g$ and the subspace $P$ independently. Consequently, the quality of the uniform approximation of $f$ by $\hat{f}$ is then bounded by the corresponding error between $g$ and $\hat{g}$ and the error between $P$ and $\hat{P}$. In what follows, we will only discuss the approximation error between $P$ and $\hat{P}$, because the approximation error between $g$ and $\hat{g}$ is well known, cf. [8] and Section 2.

We consider two different functions, one which fulfills all assumptions of Theorem 4.1, namely $g=\tanh$, and one which does not fulfill the assumption of a positive derivative, namely $g=\sin(5\cdot)$, which is not monotone on its domain. We further consider for the dimension of the subspace $d=1$ and $8$. For the dimension of the ambient space we choose $N=10$ and $50$. For each combination of $N$, $d$ and $g$ we ran $100$ experiments, where we drew a subspace $P\in G(d,N)$ uniformly at random. Note that the analysis of the algorithm makes heavy use of the monotonicity of $g$. However, the numerics show that we have comparably good results for the non-monotonic function $\sin(5\cdot)$.

5.1 Numerical Results for ATPE

The implementation of ATPE (Algorithm 2) is straightforward. However, to draw a random vector from some tangent plane, we used a method of the Matlab toolbox Manopt [19] to draw a random vector of $\mathbb{R}^{N}$ and projected it to the tangent plane. The results can be seen in Figure 4. They show as promised that the error of the approximation becomes arbitrarily small for all considered choices of $d,N$ and $g$, if we choose $h$ in computing the divided differences certainly small.

5.2 Numerical Results for OGM (Algorithm 3)

To solve the optimization problem (59), we leverage the freely available Matlab toolbox Manopt [19]. In particular, we imposed the manifold constraint by choosing the built-in grassmannfactory and we selected the built-in steepestdescent solver, as steepest descent is a well-known method to solve optimization problems. This solver requires both a cost function and the Euclidean gradient of the cost function as inputs:

[TABLE]

where $\text{{\boldmath H}}_{k}=\left[H_{1}x_{k}x_{k}\dots H_{N}x_{k}x_{k}\right]$ and $\texttt{interp1}(x,v,xq,\text{'spline'})$ returns interpolated values of a one-dimensional function at specific query points $xq$ using spline interpolation. The vector $x$ contains the sample points, and $v$ contains the corresponding function values. Note, that the Euclidean gradient ignores the manifold constraints.

Further observe, that even so OGM is shown to succeed to find an objective function whose minimizer is a suitable approximation of the wanted subspace, it is not obvious that the optimization algorithm can succeed to find this minimizer. In order to hope for this, we need to input a default subspace to the optimization algorithm which is not to far from the wanted one. In the following we choose those default subspaces uniformly at random in two different neighborhoods of the wanted subspace.

Figure 5 shows the results for both functions and the different choices of $N$ and $d$, if the default value for the optimization is a randomly chosen subspace in a distance of at most $\sqrt{2d(1-\cos(\pi/3))}$ to the unique minimum $P$ of the objective function (55), i.e., if the default value is an rotation of $P$ with an angle of at most $\pi/3$. The error is given in a logarithmic scale and the lines correspond to the different choices of $g$, $N$ and $d$. We see that we can recover all randomly drawn subspaces successfully, whenever the dimension is $d=1$ or whenever $g=\tanh$. This fits to our analysis, where the theorems hold true for injective functions and indeed $\sin(5\cdot)$ is not at all injective.

Figure 6 shows that for the case that the dimension is $d=8$ and and that we have $g=\sin(5\cdot)$, we can recover $95\%$ of randomly drawn subspaces if we ensure that the default value for the optimization is at most in a distance of $\sqrt{2d(1-\cos(\pi/4))}$ to $P$. Thus, we see that for a more carefully chosen default value, all subspaces can be recovered with a reasonable small error. Note that the case $d=1$ corresponds to the case of a usual ridge function, because if $\dim P=1$ we measure the distance to the $N-1$-dimensional subspace $P^{\perp}$.

Furthermore, as we have seen in Subsection 4.2 we can use fewer measurements to ensure almost injectivity. We then also have to adapt the convegrence analysis of $\hat{F}_{M}$.

Acknowledgements.

The author acknowledges support by the DFG Grant 1446/18 and the Berlin Mathematical School. In particular the author acknowledges Ingrid Daubechies, Gitta Kutyniok and Mauro Maggioni for helpful discussions.