Writing Real Symmetric Matrices as Linear Combination of Rank One Symmetric Terms $uu^T$

Question

Given a real symmetric matrix $M$, ostensibly of "low rank", efficiently find an expression $M = \sum \alpha_i u_i u_i^T$ using the number of terms rank($M$).

A 2011 StackOverflow Question Dense Cholesky Update in Python asked about doing low rank updates to Cholesky decompositions, but was answered largely with ways to do rank 1 updates in Python.

Naturally a rank one symmetric term would have the form $\alpha u u^T$, and if we have a rank two symmetric update of the form $uv^T + vu^T$, then we can rewrite it:

$$ uv^T + vu^T = \frac{1}{2} ((u+v)(u+v)^T - (u-v)(u-v)^T) $$

but this does not seem to lead to an efficient general method for similarly expanding an arbitrary symmetric real matrix into as many such terms as the rank requires.

Existence of such an expresssion is guaranteed by eigenbasis expansion:

$$ M = \sum \alpha_i u_i u_i^T $$

where the $u_i$ are an orthonormal basis of eigenvectors with corresponding eigenvalues $\alpha_i$ (and the limitation to as many terms as the rank requires is met by ignoring zero eigenvalues). But a full eigenbasis expansion would call for a disproportionate amount of work relative to the Cholesky updates.

We could easily expand a rank $k$ matrix as a sum of $k$ not necessarily symmetric rank one terms and use the symmetry of the sum to rewrite it as $2k$ symmetric rank one terms:

$$ \sum u_i v_i^T = \frac{1}{2} \sum (u_i v_i^T + v_i u_i^T) = \frac{1}{4} \sum ((u_i+v_i)(u_i+v_i)^T - (u_i-v_i)(u_i-v_i)^T) $$

But this uses twice the number of rank one updates as ought to be necessary.

Answer 1

Performing $k$ steps of a pivoted (rank-revealing) LDL^T factorization should give you such an approximation. This is more obvious when you recognize that an $LDL^T$ factorization is in the form of your symmetric expansion, where each column of $L$ is one of your basis vectors, say $u_i$.

Pivoted $LDL^T$ factorizations yield:

$$ P A P^T = L D L^T, $$

and so

$$A = (P^T L) D (P^T L)^T $$

provides a symmetric decomposition. You can stop the factorization process once the remaining portion of the matrix is sufficiently small in your favorite norm. The total amount of work for a rank-$k$ approximation should be $O(k n^2)$.

See Gu and Miranian's paper Strong rank-revealing Cholesky factorization for a discussion of this issue in the SPD regime.

Answer 2

Even for an indefinite real symmetric matrix, an exact symmetric expansion can be found in $O(kn^2)$ operations, where $n\times n$ matrix $M$ has rank $k$.

Proof is by induction and constructive. Because our steps reduce $n$ either by 1 or 2, base cases $n=0$ (trivial sum) and $n=1$ (scalar) are adduced by inspection.

Case 0: The simplest case is where a row and column "pair" (same index) are all zero. After any needed row and column swaps:

$$ M = \begin{bmatrix} 0 & 0 \\ 0 & C \end{bmatrix} $$

It is evident that applying the induction hypothesis to $(n-1)\times (n-1)$ symmetric $C$, modified only by tacking a zero entry to the front of each vector used, produces the required expansion using precisely rank($M$) = rank($C$) terms.

Case 1: Otherwise, if some diagonal entry of $M$ is nonzero, then a Cholesky-like step as in Jack Poulson's Answer reduces the matrix size by 1. After distributing out said nonzero diagonal entry, which may be positioned first by suitable row/column exchanges, we have:

$$ M = \begin{bmatrix} 1 & u^T \\ u & B \end{bmatrix} = \begin{bmatrix} 1 \\ u \end{bmatrix} \begin{bmatrix} 1 & u^T \end{bmatrix} + \begin{bmatrix} 0 & 0 \\ 0 & C \end{bmatrix} $$

where $C = B - uu^T$. After applying Case 0 for the last term (which tacks zero on as leading entries), the nonzero leading entry in the first term makes it linearly independent from the terms in the symmetric expansion of the latter. Thus rank($M$) is one more than rank($C$), and the combined expansion has the right number of terms.

Case 2: Finally suppose all rows and columns contain a nonzero entry but have zero on the diagonal. By distributing out the scalar and swapping rows and columns suitably, we can arrange:

$$ M = \begin{bmatrix} 0 & u^T & 1 \\ u & B & v \\ 1 & v^T & 0 \end{bmatrix} = \begin{bmatrix} 0 \\ u \\ 1 \end{bmatrix} \begin{bmatrix} 1 & v^T & 0 \end{bmatrix} + \begin{bmatrix} 1 \\ v \\ 0 \end{bmatrix} \begin{bmatrix} 0 & u^T & 1 \end{bmatrix} + \begin{bmatrix} 0 & 0 & 0 \\ 0 & C & 0 \\ 0 & 0 & 0 \end{bmatrix} $$

where $C = B - (uv^T + vu^T)$. Since as noted in the Question:

$$ uv^T + vu^T = \frac{1}{2} ((u+v)(u+v)^T - (u-v)(u-v)^T) $$

applying the induction hypothesis to $C$ and tacking on first and last zero coordinates to each vector in its symmetric expansion gives terms that may be summed with the two above. Since $uv^T$ has first row zero and last row not, and conversely for $vu^T$, these two terms are independent in combination with the terms from $C$ (which have both first and last rows zero). Thus rank($M$) is two more than rank($C$), and again the combined expansion has the right number of terms.

Work in each reduction step is $O(n^2)$, so provided we halt (in exact arithmetic) when the reduced block $C$ is recognized to be zero, $O(kn^2)$ is the total complexity.