Algorithms and Hardness for Subspace Approximation

Introduction

Large data sets that arise in data mining, machine learning, statistics and computational geometry problems are naturally modeled as sets of points in a high-dimensional Euclidean space. Even though these points live in a high-dimensional space, in practice they are observed to have low intrinsic dimension and it is an algorithmic challenge to capture their underlying low-dimensional structure. The subspace approximation problem described below generalizes several problems formulated in this context.

We describe below the special cases of the subspace approximation problem which have been studied previously and the known results about them.

For small $k$: Bădoiu, Har-Peled and Indyk gave a $(1+\varepsilon)$-approximation algorithm running in polynomial time for the minimum enclosing cylinder problem (equivalent to Subspace(1,$\infty$)), which was further extended by Har-Peled and Varadarajan to Subspace($k$,$\infty$) for constant $k$.

For $p\neq 2$, we do not know any suitable generalization of SVD, and therefore, have no exact characterization of the optimal subspace. The approximation techniques used so far to overcome this are: (i) coresets and sampling-based techniques: which give nearly optimal approximations but only for small or constant $k$ and $p$. (ii) convex relaxations and rounding: which give somewhat sub-optimal approximations mostly for large values of $k$; the only exception is the result of Varadarajan, Venkatesh, Ye and Zhang which works for any $k$ (but only for $p=\infty$).

In this paper, we study the problem Subspace($k$,$p$) for $p<\infty$, about which little is known in general. One motivation for doing so is that often the case $p<\infty$ gives significantly better approximation guarantees and requires somewhat different techniques to analyze than $p=\infty$. This is evident in the work for subspace approximation for small $k$ ( and for $p<\infty$ versus and for $p=\infty$) and in the work on regression ( and versus the $p=\infty$ case which is solvable by fixed dimensional linear programming). Also, in the study of hardness of approximation, the case $p=\infty$ can often be reduced to a discrete problem; while the case $p<\infty$ is inherently of a more continuous nature, and requires somewhat different techniques.

We also investigate the hardness of approximation for Subspace($k$,$p$). We give a reduction from the Unique Label Cover problem of Khot to the problem of approximating Subspace($n-1$,$p$) within a factor $\gamma_{p}$ (which can trivially be extended to a reduction to Subspace($k$,$p$) for $k=n^{\Omega(1)}$). The reduction is related to the ones used for similar geometric problems in , and . However, an interesting difference here in comparison to usual reductions is that we use a different (real-valued) encoding of the assignment to Unique Label Cover (in terms of the Fourier coefficients of the long-code instead of the truth table) which is more natural in our context. This may also be useful for other problems of a continuous nature.

Very recently, our techniques were also extended by Guruswami et al. to give a reduction from the Label Cover problem (without assuming the uniqueness property) to approximating Subspace($n-1$,$p$) within a factor of $\gamma_{p}$. This proves an unconditional NP-hardness for the latter problem.

Other related problems

In this special case, using Grothendieck’s inequality and a technique by Alon and Naor , one can get $O(1)$-approximation. Moreover, in this case, the above problem is also equivalent to finding diameters of convex bodies given by $\left\lVert Ax\right\rVert_{p}\leq 1$ and computing $p\mapsto 2$ norm of the matrix $A^{-1}$.

Preliminaries and Notation

2 Bernoulli and Gaussian Random Variables

A Bernoulli random variable is a discrete random variable taking values in $\{-1,1\}$ with probability $\nicefrac{{1}}{{2}}$ each. A standard normal random variables (or $1$-dimensional Gaussian) is a continuous random variable with probability density function $1/\sqrt{2\pi}\cdot\exp(-x^{2}/2)$. We use $\gamma_{p}$ to denote the $p^{th}$ moment of $N(0,1)$,

We shall require both upper and lower bounds on moments of a sum of Bernoulli random variables by the moment of an appropriate Gaussian. The following upper bound is one direction of the Khintchine inequality (see ) well-known in functional analysis.

The following version of the reverse direction, when all $c_{i}$’s are much smaller than $\left\lVert\bf c\right\rVert$, can be derived using the Berry-Esseen Theorem (as in ). A proof of the statement below appears in (as Lemma 2.5).

Technical Overview

In this section we describe our results and give a general outline of the sections that follow.

A convex relaxation of Subspace($k$,$p$) is then obtained by optimizing over arbitrary positive semidefinite matrices $X$ and replacing the requirement that the matrix have rank $n-k$ by a condition on the trace of $X$ (see Figure 1). This is similar to the relaxations used in . The problem then reduces to giving a “rounding algorithm” which reduces the rank of the matrix $X$ (which might be as large as $n$) to $n-k$, and achieves a good approximation of the objective value of the convex program.

In keeping with the intuition that the singular vectors of $ZZ^{T}$ span the orthogonal complement of $V$, our algorithm looks at the singular vectors of the matrix $X$ obtained by solving the convex relaxation. It then divides the singular vectors into $n-k$ “bins”, and constructs one vector for each bin by taking a random linear combination of vectors within each bin.

Our algorithm described in Section 4 achieves an approximation ratio of $\gamma_{q}\cdot(2-\frac{1}{n-k})^{\nicefrac{{1}}{{2}}}$ for Subspace($k$,$p$), where $q=2\cdot\lceil p/2\rceil$. (See Theorem 4.4.)

We remark that the problem of obtaining low-rank solutions to a semidefinite program was also considered by , and was addressed by simply taking random (chosen according to a Gaussian) linear combinations of the singular vectors of the relevant matrix. However, in their case, they were’ only interested in satisfying the constraints, with an error depending inversely on the rank parameter. In our case, we require a rank $n-k$ positive semidefinite matrix, all of whose eigenvalues are exactly 1. Since the only constraint enforcing this is a constraint on the trace of the matrix, even a small multiplicative error in satisfying the constraint can make some singular values quite small. To resolve this, we proceed by dividing the singular vectors in various bins and take Bernoulli linear combinations, do directly generate the orthogonal singular vectors.

2 A gap instance

3 Unique-Games hardness

In Section 6, we describe a reduction from Unique Label Cover to the problem of approximating Subspace($n-1$,$p$) within a factor better than $\gamma_{p}$ (for a constant $p\geq 1$). By a trivial reduction from Subspace($n-1$,$p$) to Subspace($k$,$p$) for any $k=n^{\Omega(1)}$, this gives the hardness of approximating Subspace($k$,$p$) better than $\gamma_{p}$, assuming the Unique Games Conjecture.

To understand the intuition for the reduction, let us consider the simpler problem of testing whether a given function $f:\{-1,1\}^{R}\to\{-1,1\}$ is a “dictator” i.e. $f(x_{1},\ldots,x_{R})=x_{i}$ for some $i\in[R]$, which is a useful primitive in such reductions. The problem is to design an instance $\cal I$ of Subspace($n-1$,$p$) and interpret the description of $f$ as a solution to $\cal I$. The required property is that if $f$ is a dictator then the corresponding subspace fits the points in $\cal I$ with small error. On the other hand, if $f$ is “far from being a dictator”, the error is required to be larger by a factor of $\gamma_{p}$.

Approximation Algorithm via Convex Programming

To relax the minimization problem for Subspace($k$,$p$) to a convex problem, we rewrite the distances $\left\lVert a_{i}^{T}Z\right\rVert$ in the objective as $(a_{i}^{T}ZZ^{T}a_{i})^{\nicefrac{{1}}{{2}}}$. Noting that $ZZ^{T}$ is a positive semidefinite matrix of rank $n-k$, we get the following natural relaxation similar to the one used in .

Note that this relaxation removes the constraint on the rank and relaxes the constraint on the length of the individual vectors $Z^{(j)}$ to the trace of entire matrix $X$. Also, the objective function is written as $\left(\sum_{i}\left|a_{i}^{T}Xa_{i}\right|^{\nicefrac{{p}}{{2}}}\right)^{\nicefrac{{1}}{{p}}}$ which is not convex. However, for solving the convex program, we can work with $\sum_{i}\left|a_{i}^{T}Xa_{i}\right|^{\nicefrac{{p}}{{2}}}$, which is convex for $p\geq 2$.

In Figure 2, we give a “rounding algorithm” for the relaxation. Note that the problem here is not really to round the solution to an integer solution as with most convex relaxations, but instead to reduce the rank of the solution to the program, while obtaining a good approximation of the objective.

We shall show the algorithm outputs a matrix $Z$ of rank $n-k$ which achieves an approximation ratio of $\gamma_{p}\cdot\sqrt{2-\nicefrac{{1}}{{n-k}}}$ in expectation, for even integers $p\geq 2$. An approximation guarantee for other values of $p$ can be obtained via Jensen’s inequality. We state the dependence on $n-k$ precisely as we shall be interested in the case $n-k=1$. For notational convenience, we shall use $\alpha_{n,k}$ to denote the quantity $\sqrt{2-\nicefrac{{1}}{{n-k}}}$ in the rest of this section.

It is clear that the columns of the matrix $Z$ given by the algorithm form an orthonormal set since they are all in the span of distinct eigenvectors of $X$, and are normalized to have length 1. However, this assumes that the lengths of the vectors $y_{j}$ are nonzero. Since a vector $y_{j}$ is a weighted sum of orthogonal vectors, $\left\lVert y_{j}\right\rVert^{2}=\sum_{t\in S_{j}}\lambda_{t}$. The following claim gives a lower bound on this quantity which is also useful in bounding the approximation ratio.

Let $S_{1},\ldots,S_{n-k}$ be the partition constructed by the algorithm in step 2. Then

Proof: Let ${j_{0}}\stackrel{{\scriptstyle\textup{def}}}{{=}}\operatorname{argmin}_{j}\left\{\sum_{t\in S_{j}}\lambda_{t}\right\}$ and let $s^{*}\stackrel{{\scriptstyle\textup{def}}}{{=}}\sum_{t\in S_{j_{0}}}\lambda_{t}$. Let $Q\stackrel{{\scriptstyle\textup{def}}}{{=}}\{j_{0}\}\cup\{j~{}|~{}|S_{j}|>1\}$. Note that the algorithm ensures that $|S_{j}|>0$ for all $j$ but in $T$ we discard the singleton sets. We will show that $s^{*}\geq 1/\left(2-\nicefrac{{1}}{{|Q|}}\right)$, which will prove the claim since $|Q|\leq n-k$.

We argue that for each $j\in Q$, $j\neq j_{0}$, $\sum_{t\in S_{j}}\lambda_{t}~{}\leq~{}2s^{*}$. To see this, let $t_{j}$ be the maximal index in $S_{j}$. At step $t=t_{j}$, $t_{j}$ was added to set $S_{j}$ and not to the set $S_{j_{0}}$. Hence,

Also, there exists at least one $t_{0}\in S_{j_{0}}$ such that $t_{0}<t_{j}$. This is because $S_{j}$ was non-empty at step $t_{j}$ (otherwise it would be a singleton). But then $\lambda_{t_{j}}\leq\lambda_{t_{0}}\leq s^{*}$ and, hence, $\sum_{t\in S_{j}}\lambda_{t}\leq 2s^{*}$.

Finally, we note that for each $j\notin Q$, $S_{j}$ contains exactly one element $t$, the eigenvalue $\lambda_{t}$ corresponding to which is at most 1. Thus,

The following lemma proves the required approximation guarantee for the expected $p^{th}$ moment of the distance a single point $a_{i}$ from the orthogonal complement of the column span of $Z$.

Let $X$ be the solution of the convex relaxation and let $Z$ be the matrix returned by the algorithm. Also, let $p$ be even. Then, for each $i\in[m]$

Proof: We can expand $\left\lVert a_{i}^{T}Z\right\rVert$, using $W_{j}$ to denote $\left\langle a_{i},Z^{(j)}\right\rangle$, as

Note that the $W_{j}$-s are independent random variables since each $W_{j}$ only depends on $b_{t}$ such that $t\in S_{j}$, and the sets are disjoint. Using the multinomial expansion and the fact that $p$ is even, the above can be written as

The following claim then finishes the proof.

For each $j$, let $D_{j}$ denote $\sum_{t\in S_{j}}\lambda_{t}\left\langle a_{i},x_{t}\right\rangle^{2}$ and let $\Lambda_{j}$ denote $\sum_{t\in S_{j}}\lambda_{t}$. Using the above claim we get that

An approximation guarantee for other values of $p$ can be obtained via a standard application of Jensen’s Inequality. We state the dependence on $n-k$ precisely as we shall be interested in the case $n-k=1$ in the later sections. Notice that the approximation factor is $\gamma_{q}$, where $q=2\cdot\lceil p/2\rceil$, in the case $n-k=1$, and thus matches the integrality gap and unique-games hardness that appear in the later sections.

Let $X$ be the solution of the convex relaxation and let $Z$ be the matrix returned by the algorithm. Let $p\geq 1$ and let $q=2\cdot\lceil\nicefrac{{p}}{{2}}\rceil$ be the smallest even integer such that $q\geq p$. Then,

Proof: (Proof of Theorem 4.4) By the concavity of the function $f(u)=u^{\nicefrac{{1}}{{p}}}$ and Jensen’s Inequality we have that

and by linearity it suffices to consider a single term of the summation. Another application of Jensen’s (using $p\leq q$) and Lemma 4.2 give that

which completes the proof of the theorem.

Our results are stated in terms of the expected approximation ratio achieved by the algorithm. However, one can get arbitrarily close to this ratio with high probability, simply by considering few independent runs of the algorithm and picking the best solution. In particular, one can achieve an approximation guarantee $(1+\varepsilon)\cdot\gamma_{q}\cdot\sqrt{2-(\nicefrac{{1}}{{n-k}})}$ with probability $1-p_{e}$, by using $O(\nicefrac{{1}}{{\varepsilon}}\cdot\log(\nicefrac{{1}}{{p_{e}}}))$ runs.

A Gap Instance for the Convex Relaxation

Here we describe an instance of Subspace($n-1$,$p$) such that the value of any valid solution (which is of rank $1$) is at least $\gamma_{p}$ times the value of the convex relaxation. Note that approximation ratio of the algorithm for the case $n-k=1$ (and even $p$) is exactly $\gamma_{p}$ and hence this shows that our analysis is optimal for this case.

Proof: We first consider the value of the LHS. By the rotational invariance of the Gaussian measure, the value is equal for all $z$ and we can restrict ourselves to $z=e_{1}$.

In comparison, the optimum of the convex relaxation can be upper bounded by using the matrix $X=\nicefrac{{1}}{{n}}\cdot I$.

2 Discretizing the gap example

A discrete analog of the above, i.e., picking sufficiently many samples from the same distribution, gives us our final integrality gap (or “rank gap”) example.

The theorem can be proved by using the continuous gap instance, and concentration bounds for the samples $a_{1},\ldots,a_{m}$. We defer a full proof to the appendix.

Unique-Games Hardness

We shall show a reduction to subspace approximation problem from the Unique Label Cover problem defined below.

An instance of Unique Label Cover with alphabet size $R$ is specified as a bipartite graph ${{\mathcal{U}}=(V,W,E)}$ with a set of permutations $\{\pi_{vw}:[R]\to[R]\}_{(v,w)\in E}$. A labeling ${\cal L}:V\cup W\to[R]$ is said to satisfy an edge $(v,w)$ if ${\cal L}(w)=\pi_{vw}({\cal L}(v))$. We denote by $\mathsf{val}({\cal U})$ the maximum fraction of edges satisfied by any labeling ${\cal L}$.

The Unique Games Conjecture proposed by Khot in conjectures the hardness of distinguishing between the cases when the optimum to the above problem is very close to 1 and when it is very close to 0. This conjecture is an important complexity assumption as several approximation problems have been shown to be at least as hard as deciding if a given instance ${\cal U}$ of Unique Label Cover problem has $\mathsf{val}({\cal U})>1-\varepsilon$ or $\mathsf{val}({\cal U})<\delta$ for appropriate positive constants $\varepsilon$ and $\delta$.

Given any constants $\varepsilon,\delta>0$, there is an integer $R$ such that it is NP-hard to decide if for given an instance ${{\mathcal{U}}=(V,W,E)}$ of Unique Label Cover with alphabet size $R$, $\mathsf{val}({\cal U})\geq 1-\varepsilon$ or $\mathsf{val}({\cal U})\leq\delta$.

2 Reduction from Unique Label Cover

Norms for functions are defined as usual (over the uniform probability measure). Note that $\left\lVert f_{b}\right\rVert_{2}^{2}=\left\lVert{\bf b}\right\rVert_{2}^{2}$. When the exponent in the norm is unspecified, $\left\lVert\cdot\right\rVert$ denotes $\left\lVert\cdot\right\rVert_{2}$.

Given an instance ${{\mathcal{U}}=(V,W,E)}$ of Unique Label Cover we output the following instance of subspace approximation, for a suitable constant $B$ to be determined later:

The following claim shows that the optimum of the subspace approximation problem is low when the Unique Label Cover is instance is highly satisfiable.

If $\mathsf{val}({\cal U})\geq 1-\varepsilon$, then $(\mathsf{opt})^{p}\leq 1+\varepsilon\cdot B\cdot 2^{p}.$

Note that $\left\lVert f_{\pi_{w_{2}v}(b_{w_{2}})}-f_{\pi_{w_{1}v}(b_{w_{1}})}\right\rVert_{p}^{p}$ equals $2^{p-1}$ if $\pi_{w_{1}v}({\cal L}(w_{1}))\neq\pi_{w_{2}v}({\cal L}(w_{2}))$ and 0 otherwise. Hence,

Combining the two bounds above gives $(\mathsf{opt})^{p}\leq 1+\varepsilon\cdot B\cdot 2^{p}$.

Soundness

For the soundness, we need to prove that if $\mathsf{val}({\cal U})\leq\delta$, then $\mathsf{opt}\geq\gamma_{p}^{p}\cdot(1-\nu)$ where $\nu$ is a small constant depending on $\varepsilon$ and $\delta$. We first make some simple observations about the optimal solution.

For any optimal solution $\{{\bf b}_{w}\}_{w\in W}$ to the above instance of Subspace($n-1$,$p$), it must be true that

We show that if $\mathsf{val}({\cal U})\leq\delta$, then in fact the first term itself is approximately $\gamma_{p}^{p}$. As is standard in Unique Games based reductions, the proof proceeds by arguing separately about the “high-influence” and “low-influence” cases. However, since the inputs for our problem are not in the form of a long-code but the vectors ${\bf b}$, we will use $\max_{i\in R}\{\left|b_{i}\right|/\left\lVert{\bf b}\right\rVert\}$ as a substitute for influence of the $i^{th}$ variable on the function $f_{b}$.

For the vertices $v\in V$ where the functions $f_{b_{v}}$ have no influential coordinates, the Central Limit Theorem shows that $\left\lVert f_{b_{v}}\right\rVert_{p}$ is very close to $\gamma_{p}$. We then show that the contribution of the remaining vertices to the objective function is small.

Below, we define $S_{1}$ to be the set of vertices corresponding to low influence functions and divide the remaining vertices into three cases which we shall analyze separately. The parameters $\tau,\beta\in(0,1/2)$ will be chosen later.

Since $f_{b_{v}}(x_{1},\ldots,x_{R})=b_{v,1}\cdot x_{1}+\cdots+b_{v,R}\cdot x_{R}$ is a linear function of Bernoulli variables, Claim 2.3 gives that

as $\left\lVert{\bf b}_{w}\right\rVert=\left\lVert\pi_{wv}({\bf b}_{w})\right\rVert$ and that ${\bf b}_{v}$ is the mean of $\pi_{wv}({\bf b}_{w})$. Now, since $\left\lVert{\bf b}\right\rVert=\left\lVert f_{b}\right\rVert$, we get that

Proof: Consider a vertex $v\in S_{3}$. Since we know that $v\notin S_{2}$, we get that

Using this we can again say that $\left\lVert{\bf b}_{v}-\pi_{wv}({\bf b}_{w})\right\rVert$ must be large on average and, hence, derive a bound on the measure of $S_{3}$.

As in the previous claim, we use this to conclude that

Proof: Since $v\notin S_{1}\cup S_{2}\cup S_{3}$, we know that

Construct a labeling for ${\cal U}$ by assigning to each $v\in V$, the special label $i$ as above, and to each $w\in W$, a random label $j$ satisfying $\left|b_{w,j}\right|\geq\nicefrac{{\tau}}{{4}}\cdot\left\lVert{\bf b}_{w}\right\rVert$. For, $w\in W$ when no such $j$ exists or for $v\notin S_{4}$, we fix a label arbitrarily.

Note that there can be at most $\nicefrac{{16}}{{\tau^{2}}}$ choices of $j$ satisfying $\left|b_{w,j}\right|\geq\nicefrac{{\tau}}{{4}}\cdot\left\lVert{\bf b}_{w}\right\rVert$. By the condition on $i$, we know that, in expectation, the labeling satisfies $\nicefrac{{1}}{{4}}\cdot\nicefrac{{\tau^{2}}}{{16}}$ fraction of the edges incident on a $v\in S_{4}$. Since the fraction of edges satisfied overall is at most $\delta$, we get that

Let $\nu$ denote $10\tau\cdot\left(\log(\nicefrac{{1}}{{\tau}})\right)^{\nicefrac{{p}}{{2}}}$. Using these estimates, we can now prove the soundness of the reduction.

If $\mathsf{val}({\cal U})<\delta$, then for the reduction with parameters $B,\tau$ and $\beta=\tau^{2}$

We lower bound ${\mathds{1}}_{\{S_{1}\}}(v)$ by $1-{\mathds{1}}_{\{S_{2}\}}(v)-{\mathds{1}}_{\{S_{3}\}}(v)-{\mathds{1}}_{\{S_{4}\}}(v)$. Claims 6.5 and 6.6 and give bounds on the first two terms (with $\beta=\tau^{2}$).

We bound the third term using Claim 6.7 and Hölder’s inequality

where the last bound used that since $\mathsf{opt}\leq\gamma_{p}$ (see Claim 6.4), we must have

Combining the bounds for the above three terms proves the lemma.

For a small constant $\eta$ such that $\eta(\log(\nicefrac{{1}}{{\eta}}))^{\nicefrac{{p}}{{2}}}<\nicefrac{{2^{-\nicefrac{{p}}{{2}}}}}{{50}}$, choosing parameters as

in Lemma 6.8 would imply that $\mathsf{opt}\leq 1+\eta$ in the completeness case and $\mathsf{opt}\geq\gamma_{p}\cdot(1-\eta)$ in the soundness case. This gives the following theorem.

For any $p\geq 2$ and sufficiently small constant $\eta$, there exist constants $\varepsilon,\delta>0$ and a reduction from Unique Label Cover to Subspace($n-1$,$p$) such that if $\mathsf{val}({\cal U})$ is the fraction of edges satisfiable in the given instance of Unique Label Cover and $\mathsf{opt}$ is the optimum of the instance of Subspace($n-1$,$p$), then

Acknowledgments

We thank Kasturi Varadarajan for initiating the work on this problem by suggesting that we generalize the algorithm of and for generous help with an early draft of this paper on which he was offered a co-authorship (but later opted out). MT would also like to thank David Steurer for helpful discussions. We thank the anonymous reviewers of this manuscript for their suggestions and references, and for pointing out an error in the previous proof of Lemma 6.8.

References

Appendix A Proof of Theorem 5.2

as long as we choose $m$ large enough so that

Hence, choosing $m>\nicefrac{{5}}{{\varepsilon^{2}}}$, we have

On the other hand to analyze the value of the corresponding convex relaxation, we use

Choosing $m>\nicefrac{{9}}{{\eta^{2}}}$, we get

Therefore, the convex relaxation satisfies

Appendix B NP-hardness of Subspace Approximation

In this section, we show unconditionally that the problem Subspace($n-1$,$p$) is NP-hard, for $p>2$, using a reduction from the Min-Uncut problem on graphs. Such a result was also obtained independently by Gibson and Xiao (personal communication).

Min-Uncut problem: Given a graph $G=(V,E)$, find a bipartition of its vertices $V=S\cup T$ that minimizes the number of edges with both endpoints on the same side of the bipartition.

where $N$ is an integer polynomially large in $n$ and $m$ which will be chosen later.

No case: Otherwise, for any bipartition the Min-Uncut has at least $t+1$ edges, i.e., for any $x\in\{-1,1\}^{n}$ we have $\sum_{ij\in E}(x_{i}+x_{j})^{p}\geq(t+1)2^{p}$. Now divide the sphere of radius $\sqrt{n}$ into two parts as follows:

where $\varepsilon<\nicefrac{{1}}{{p}}\cdot(m+1)$. For any $y\in T$,

Case $1$: $\left|y_{i}\right|=1+\varepsilon_{i}\geq 1+\varepsilon$ for some $i$. Then,

Case $2$: $\left|y_{i}\right|=1-\varepsilon_{i}\leq 1-\varepsilon$ for some $i$. Then,

Therefore, using the same analysis as in the previous case, we get

Using the above property of $y\in T$, we get