# Choosing subset of vectors to approximate a subspace

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.

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

• Is this the recent first order convex paper you're referring to: Nesterov, Y. E., Smooth minimization of nonsmooth functions ? It seems to solve these sort of problems well and is highly cited. – Nick Alger Sep 11 '12 at 19:39
• This is just the tip of an iceberg. The Nesterov paper does not directly address this problem (but a more general class), but some of the later work done in this area is dedicated to group regularization (as the above problem is called). See scholar.google.at/… – Arnold Neumaier Sep 12 '12 at 9:50
• Do you think these Nesterov methods would be better than Bregman optimization? Eg, caam.rice.edu/~optimization/L1/bregman – Nick Alger Sep 12 '12 at 18:43
• Note that the regularization term must be a sum of 2-norms to have the desired group regularization effect! The papers on the Rice site you mention don't address the group regularization problem. – Arnold Neumaier Sep 13 '12 at 9:08
• Thank you for the help. Group regularization looks like what I want, but it will take me a while to read and understand these papers before I can say for sure. – Nick Alger Sep 14 '12 at 6:16

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.

• The first idea is out because I need a strict subset of the original vectors and not linear combinations. The second idea seems promising, but I forsee a problem - if the vectors have a very large component perpendicular to $V$, then they might be unsuitable even if they satisfy the conditions youve listed. Perhaps normalizing first would fix this problem though.. – Nick Alger Sep 8 '12 at 21:08
• Yes, throw out the ones that are almost perpendicular first (these are the ones for which $\|\Pi_V x_i\| \ll \|x_i\|$) and only compute the angles for the remaining ones. – Wolfgang Bangerth Sep 9 '12 at 2:38
• Thinking about this more generally, the goal is basically to balance the goals of maximizing the projections of the vectors onto the space, and minimizing the projected angles between the vectors. – Nick Alger Sep 9 '12 at 3:05

### 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.)

• It's very important that it's a strict subset of the vectors and not linear combinations. However if there were linear combinations of only $k=100$ vectors that well-approximated the space, then choosing that subset of the vectors would, of course, suffice. – Nick Alger Sep 8 '12 at 21:04
• Okay, we know that only 20 vectors are needed to span the space. Presumably your collection of $n$ vectors is sufficient to span the space. Now do you want (a) the best-conditioned collection of $k=20$ vectors that spans $V$, (b) $k=20$ vectors minimizing some norm of their component orthogonal to $V$, or (c) something else? – Jed Brown Sep 8 '12 at 21:12
• The original vectors probably don't span the space exactly, and conditioning is not important. The goal is something along the lines of the following: let $\tilde V$ be the approximation of $V$, and let $x$ be a uniformly chosen random unit vector on $V$. We want to maximize $\mathbb{E} ||\Pi_{\tilde V} x||$, the expected value of the projection of $x$ onto $\tilde V$. There may also be other better characterizations of "goodness of approximation", but that's what I've thought of so far. – Nick Alger Sep 8 '12 at 21:23
• I updated my answer. Please put the information from your comments into your question. If you are having trouble coming up with a useful measure of quality, you should explain what the problem background and what you are going to use the selected vectors for. Note that it is very rare for basis conditioning to be unimportant for approximation problems. – Jed Brown Sep 8 '12 at 22:00
• The algorithm in the updated post will probably do well, but some of the post misinterprets the problem. $X \neq span(x1,x2,...,x_n)$, and the small space $V$ does not lie in $span(x_1,...,x_n)$. As an illustrative example, consider $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)$. Here it is fruitful to choose two vectors $x_1$ and $x_2$ to approximate a 1-dimensional space $V$. – Nick Alger Sep 8 '12 at 23:16