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.

  • 1
    $\begingroup$ What is wrong with the eigenvector expansion? You know that if the matrix is low-rank, then most of the eigenvalues are zero and all you have to find are the few nonzero eigenvalues and corresponding eigenvectors! $\endgroup$ Jun 24, 2013 at 22:06
  • $\begingroup$ @WolfgangBangerth: The eigenvector expansion is mathematically perfect, but I fear the cost of computation destroys any efficiency one might obtain from using Cholesky updates (over refactoring). The goal is a cost comparable to the updates, i.e. $O(kn^2)$, where $M$ is $n \times n$ and rank($M$) = $k$. $\endgroup$
    – hardmath
    Jun 25, 2013 at 1:06

2 Answers 2


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.

  • $\begingroup$ Thanks for the reference. Doing Cholesky on the incremental matrix agrees with my intuition, but I wasn't sure if an "exact" expansion is feasible. $\endgroup$
    – hardmath
    Jun 25, 2013 at 1:10
  • 1
    $\begingroup$ $L$ is unit-diagonal, so its columns must be linearly independent...My comment on approximations was referring to the more general, and significantly harder, case where the matrix is only approximately low-rank. $\endgroup$ Jun 29, 2013 at 0:30

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.