# How to approx. decompose a sym. p.d. matrix M into X'X?

M: pxp symmetric p.d. matrix with unit diagonals

n: number much smaller than p

Want a nonrandom nxp matrix X such that X'X is close to M element-wise. If n gets larger, hopefully difference gets smaller.

I have a no so good method. Get the eigendecomposition M = D'VD, and rearrange the diagonals of V in decreasing order and so the rows of D. Take only the first n rows of D as X. This method works only if the eigenvalues decrease very fast.

Any input is appreciated. Thanks a lot!

-

Unfortunately, you can't do any better than your "no so good method." It's a standard result that this is the best (as measured by the Frobenius norm) rank $n$ approximation to $M$.

-
Thanks. Frobenius norm may not be the good measure for me, since I care more about the diagonals. –  Peter May 4 '11 at 4:04
Then you'll need to make your question more precise by deciding exactly what it is you want to optimize... –  Brian Borchers May 4 '11 at 4:24
Incidentally, the same result holds for the Euclidean induced norm and for any orthogonal-invariant norm (magic words: Eckart-Young-Mirsky theorem). –  Federico Poloni May 4 '11 at 6:53
Thanks. I do not know what norm is appropriate, but at least some random approximations seem "better" than this one. For example, generate a random nxp matrix from $N(0, 1/\sqrt{n})$ and multiply by $L$ on the right, where $L'L=M$ . The element-wise convergence is ensured by the law of large numbers. But for the eigen-decomposition method, if the eigenvalues are all around 1, $X'X$ are too far from $M$. –  Peter May 4 '11 at 7:47
What do you use as a measure of "distance" between $M$ and $X′X$? If you take the Euclidean or Frobenius norm, then there must be something wrong in your computations since the Eckart-Young-Mirsky theorem says that the method you suggested in the question gives the best bound achievable. If you use another distance, let us know and maybe we can provide a better solution. –  Federico Poloni May 4 '11 at 11:32

Not a new solution but an elaboration from a computational viewpoint.

From a practical point of view one can solve your problem as follows:

Firstly, we do not care about the diagonal values being $1$. Secondly, we enlarge the problem to real symmetric matrices with no strictly negative eigenvalues.

Such a matrix $A$ of size $p\times p$ is then given by $A=\sum_i \lambda_i v_i^tv_i$ where $v_1,\dots,v_p$ is an orthogonal basis of eigenvectors with eigenvalues $\lambda_1\geq \lambda_2\geq\dots$ of $A$.

Suppose there is a fast way for computing an eigenvector $v_1$ (of norm $1$) associated to the largest eigenvalue $\lambda_1$ of $A$. (One can for instance consider the projective limit of $A^k v$ for $v\in\mathbb R^p$ a generic vector.)

The largest eigenvalue $\lambda_1$ is then given by $$\lambda_1=\sum_{i,j}A_{i,j}{P(1)}_{i,j}=\sum_{i,j}A_{i,j}{(v_1)}_i{(v_1)}_j$$ where $P(1)=v_1^t v_1$ is the orthogonal rank one projector associated to $v_1$.

Replacing $A$ by $A-\lambda_1P_1$ and iterating one gets $P_2,P_3,\dots$ and $\lambda_2,\lambda_3,\dots$.

Your solution is then given by $\sum_{i=1}^p\lambda_iP(i)=\sum_{i=1}^p\lambda_iv_i^tv_i$.

-
Sorry for the teX problem which I was unable to solve. –  Roland Bacher May 4 '11 at 11:15
This idea is widely studied in numerical linear algebra; unfortunately, often it does not give computationally sound results after the first few eigenvectors since orthogonality w.r.t. the first $v_i$'s is a property that is easily lost numerically. Better solutions are the orthogonal iteration and in general Krylov subspace methods, which try to compute all the eigenvalues at the same time rather than one after the other. –  Federico Poloni May 4 '11 at 11:31