Choosing subset of vectors to approximate a subspace

Question

Suppose I have a high-dimensional vector space $X$, a subspace $V \subset X$, and a collection of $n$ vectors $\{x_i\}_{i=1}^n \subset X$.

My question is: How can I choose a small collection $k < n$ of the vectors $x_i$ so that the span of this smaller collection "well-approximates" the subspace $V$?

The notion of "well-approximation" is intentionally left vague since, although it's intuitive that some subspaces approximate each other better than others, it's not clear to me the best way to introduce definitions that make this precise.

For concreteness, in my scenario the sizes of the various objects are of the following orders $dim(X)\approx 10000$, $dim(V)\approx 20$, $n\approx 5000$, and $k$ can be varied but has a target of $k \approx 100$.

It seems like this should be well studied, but I'm having trouble finding the right terms to search for. In particular, the subject of "subspace approximation" appears to deal with the opposite problem of choosing a subspace to approximate vectors, and the topic of "basis selection" appear to be interested with choosing linear combinations of basis vectors that make certain things sparse - both very different problems from this (as far as I can tell).

Edit: some clarifications based on discussion below

• The dimension of the space $X$ is larger than the number of candidate basis vectors $x_i$, and the subspace $V$ does not necessairily lie in the span of the $x_i$'s.
• As an illustrative example of where it might be useful to consider more basis vectors than the dimension of the space being approximated, consider the following situation: $X=\mathbb{R}^4$, $V=span((1,0,0,0))$, $x_1=(1,1,\epsilon,0)$, $x_2=(1,-1,\epsilon,0)$, $x_3=(0,0,0,1)$. It would be useful to choose 2 vectors $x_1$ and $x_2$, even though the space to be approximated, $V$, has dimension 1.
• Or in 3D, consider the situation in the following picture. You can approximate the space 1D $V$ perfectly with 3 vectors $x_1,x_2,x_3$, very well with 2 vectors $x_1,x_2$, and poorly with only one vector.
• One possible measure of how well a candidate space $\tilde V$ approximates the target space $V$ would be the expected value of the size of the projection of a random unit vector in $V$ onto $\tilde V$. Ie, for a uniformly distributed random unit vector $v \in V$, maximize $\mathbb{E}||\Pi_{\tilde V} v||$. If the approximation is exact this will be 1, otherwise it will be less than 1. Other definitions of "well approximation" may be better, this is just the first thing I thought of.

Answer 1

Form the matrix $X:=[x_1 ... x_N]$ and the matrix $V$ whose columns define the subspace of interest. Thus I use $X$ and $V$ to denote matrices rather than their column spaces.

An appropriate goal seems to be the solution of the problem $XU \approx V$ with $U$ being of reasonable size and having at most $k$ nonzero rows. Without the restriction on the number of rows, the problem is both under- and overdetermined, and can be represented as the least squares problem
$~~~~~~~~~~$ (1) $~~~~$ minimize $\|XU-V\|_F^2$ in the Frobenius norm, subject to a suitable regularization.
In the present case, the restriction on the number of rows serves to regularize the problem. This problem can be solved as follows.

Compute an orthogonal factorization $V=:Q_0\pmatrix{R_0\\0}$ with $Q_0$ a product of reflections, and conformally partition $Y:=Q_0^TX=:\pmatrix{Y_1\\Y_2}$. Then the normal equations reduce to solving
$~~~~~~~~~~$ (2) $~~~~$ $Y U = \pmatrix{R_0\\0}~~~~~~~~$.
The permutation $P_1$ of an orthogonal factorization $Y_1=Q_1R_1P_1$ with column pivoting defines a set of $dim~ V$ indices of columns of $X$ needed so that the projected vectors form a good basis of $V$.

Permute $Y$ accordingly and denote the result again by $Y$. Conformally partition $Y=\pmatrix{Y_{11}& Y_{12}\\Y_{21}& Y_{22}}$. By construction, $Y_{11}$ is nonsingular (and usually well-conditioned).

Write conformally $U = \pmatrix{U_1\\U_2}$, solve the first block of the normal equations (2) for $U_1$, and substitute it into the second block of equations. This gives a smaller least squares problem of the form (1). Thus we can proceed recursively until the desired number of vectors has been determined.

Alternatively, one can solve the optimization problem
$~~~~~~~~~~$ (3) $~~~~$ minimize $\|XU-V\|_F^2+\lambda\sum_i\|U_{i:}\|_2$.
with a suitable regularization parameter $\lambda>0$. The lack of the square in the penalty terms tends to force unneeded rows of $U$ to zero. One starts with a fairly large value of $\lambda$ (forcing few nonzero rows) and uses the result as a starting point for a problem with smaller $\lambda$ (if $U$ was too sparse) or with larger $\lambda$ (if $U$ was not sparse enough).

Problem (3) is nonsmooth but convex, and there are a number of recent efficient (''first order convex'') methods for handling these at the scale given in the original problem.

Edit: The modern era of efficient convex first order methods started with Nesterov's paper ''Smooth minimization of nonsmooth functions'' http://www.springerlink.com/content/c7wheayx86u1vfwx/ . Some of the later work done in this area is dedicated to group regularization (as problem (3) is called). See http://scholar.google.at/scholar?q=group&btnG=&hl=en&as_sdt=2005&sciodt=0%2C5&cites=14895856039931556754&scipsc=1

Answer 2

I would generate vectors $y_i = \Pi_V x_i, i=1\ldots n$ where $\Pi_V$ is the projector into the subspace $V$. Not all of these $y_i$ will be mutually linearly independent, but you can of course apply the Gram-Schmidt orthogonalization procedure to the $y_i$ which will produce a set of mutually orthogonal vectors $z_i,i=1\ldots n'$ where $n'<\textrm{min}\{n,\textrm{dim}(V)\}$. These vectors $z_i$ are then the best basis for $V$ you can find from the subspace spanned by the $x_i$.

Alternatively, if you don't want to apply the Gram-Schmidt procedure, you can consider pairs of vectors $x_i,x_j$ and eliminate one of them if the angle $\theta_{ij}$ they form in $V$, with $\cos \theta_{ij} = \frac{(\Pi_V x_i) \cdot (\Pi_V x_j)}{\|\Pi_V x_i\| \|\Pi_V x_j\|}$, is too small. You also need to eliminate vectors that are mostly orthogonal to $V$, i.e. for which $\|\Pi_V x_i\|$ is below a certain threshold.

Answer 3

Update

The following greedy algorithm should do what you want reasonably fast. $\DeclareMathOperator{\argmax}{arg\,max} $

Assume the original vectors are normalized. Compute $$ y_1 = \argmax_{x_i} \lVert \Pi_V x_i \rVert . $$ Now, given $Y_j = [y_1|\dotsb| y_j]$, let $$y_{j+1} = \argmax_{x_i \notin Y_j} \lVert (1 - \Pi_{Y_j}) \Pi_V x_i \rVert . $$ At each step, this algorithm selects the $y_j$ that "fills" the most of the remaining part of $V$.

Note that even though $\hat X = [x_1| x_2| \dotsb | x_n]$ may be full rank (implying that all $x_i$ are linearly independent, $\Pi_V \hat X$ has rank at most $\dim(V)$. This means that if the above algorithm selects $\dim(V)$ vectors, then $V$ has been spanned. Given that you don't care about basis conditioning, there is no use in considering more than $\dim(V)$ vectors.

Original

How is the subspace $V$ defined? In terms of data consisting of vectors that live in $V$? You state that $\dim V \approx 20$, but you are willing to use $k \approx 100$ vectors to "approximately span" it? Are you really insisting on using a subset of the vectors $x_i$ or would linear combinations do?

Have you looked at principle component analysis? Computationally, this will involve computing a partial SVD of the data matrix (the vectors mentioned above). For your problem size, the SVD can be computed directly, but if the size grows much, you can approximate it using iterative methods (e.g. in SLEPc) or the randomized methods discussed in this question. (PCA is not quite what you asked for, but there was enough ambiguity in the question that it may be what you want.)