For a computational representation theory program, I have to calculate the trace of $A^+B$ for a fixed $A$ and many different $B$, where $A^+$ is the pseudoinverse of $A$ ($A,B$ have full column rank and the same image). Currently, I compute the LU decomposition of $A$ and use this to solve $Ax_i=b_i$ for each column $b_i$ of $B$, and use $\mathrm{tr}(A^+B)=\sum_i x_{ii}$.

A factor of 2 or so speedup could be achieved by halting backsolve once $x_{ii}$ has been computed rather than computing all of $x_i$. However, this is not enough to make a difference in my code and would require reimplementing UMFPACK methods, which would probably give me slower code in the end because my implementation wouldn't be as good.

Is there a (hopefully asymptotically) faster approach? For example, is it possible to compute $x_{ii}$ asymptotically faster than $x_i$?

In case it matters, all matrices and methods in question are sparse.

  • $\begingroup$ Are $A,B$ square matrices? If $A$ has full column rank and is square, then $A^+=A^{-1}$. $\endgroup$ Dec 2 '13 at 3:11
  • $\begingroup$ @WolfgangBangerth No, they are not. $\endgroup$ Dec 2 '13 at 3:18
  • $\begingroup$ Does $A$ have more rows or more columns? Is it of full row or column rank? You've stated that the matrices are sparse- are they small enough that it would be reasonably possible to store (a presumably dense) $A^{+}$ if it could reasonably be computed? How many $B$ matrices are there compared with the dimension of $A$? $\endgroup$ Dec 2 '13 at 5:39
  • $\begingroup$ @BrianBorchers $A$ and $B$ have the same dimensions, more rows than columns and full column rank. I could probably feasibly store $A^+$, it wouldn't be much larger than $5,000\times10,000$. The number of $B$ matrices varies, but is very large: it is roughly $$\sum_{k=1}^n p(k)$$ where $p(k)$ is the number of partitions of $k$, for values of $n$ up to $15$ or so, which gives $683$. The dimension of $A$ is generally in the thousands, with something like twice as many rows as columns. $\endgroup$ Dec 2 '13 at 5:46

Since $A$ has full column rank,


You could handle this computation by computing the Cholesky factorization of $(A^{T}A)=R^{T}R$, then compute $(A^{T}A)^{-1}=R^{-1}R^{-T}$.

If you can afford to store $A^{+}$ (which is probably dense), then once you have $A^{+}$, computing the trace of $A^{+}B$ takes $O(mn)$ time for each $B$ matrix (and a lot less if $B$ is sparse- you haven't said whether the $B$ matrices are sparse.)

If you can't afford to store $A^{+}$ but can store $(A^{T}A)^{-1}$ then it may be helpful to precompute the inverse of $A^{T}A$.

  • $\begingroup$ Thanks for the answer. As I said, all matrices in question are sparse; this includes $B$. I in fact suspect that $A^+$ would be sparse, since the density of $A$ is $\approx .001$. $\endgroup$ Dec 2 '13 at 6:08
  • 1
    $\begingroup$ Wouldn't you want to solve $R^{T}X = A^{T}$, then $RA^{+} = X$ instead of doing explicit inversion and multiplication, to avoid ill-conditioning and loss of accuracy? Furthermore, sparse QR decomposition of $A$ (say, using SuiteSparseQR, by the same author of UMFPACK) will produce a sparse $R$, and it should be possible to use sparse LU decompositions to solve these linear systems. Such an approach should save memory compared to explicit inversion. $\endgroup$ Dec 2 '13 at 9:42
  • $\begingroup$ The Sparse QR methods require you to have the right hand sides for all of the least squares problems in place when you compute the "Q-less" QR factorization. It seems unlikely that all of the B matrices could be kept in storage. $\endgroup$ Dec 2 '13 at 15:06
  • $\begingroup$ Why would you need $B$ to do a QR decomposition of $A$? I am not talking about solving a linear system with QR directly; $R$ just happens to be the coefficient matrix of the linear systems I mentioned. $\endgroup$ Dec 2 '13 at 18:11
  • $\begingroup$ Sparse QR factorization methods avoid forming the (big and dense) Q matrix by keeping track of the orthogonal transformations that produce the QR factorization. Storing this information would require as much storage as $Q$, but if we apply the same transformations to the right hand side of the least squares problem as $R$ is computed, then we can avoid this storage cost. You can do this for multiple right hand sides, but only if you have all of the right hand sides in hand as you perform these transformations. $\endgroup$ Dec 2 '13 at 21:01

I agree with BrianBorchers' general approach. To summarize my discussion with him in the comment thread below his answer, I suggest some modifications to his approach:

To recap the gist of Brian's answer, the quickest way to calculate the trace of $A^{+}B$ for many values of $B$ is to precompute $A^{+}$, and then perform only the matrix-vector products needed for the diagonal elements of $A^{+}B$, and sum these diagonal elements. In the limit of large numbers of $B$, the expense of those matrix-vector products will dominate the cost of forming $A^{+}$. $A$ is also sparse, so you could exploit sparsity for speed, too. Memory becomes a concern when $A$ is large. You mentioned that $A$ is no larger than $5000 \times 10000$, so with double precision arithmetic, it could require up to 400 MB to store $A$. That's a lot, but storing it as a dense matrix is possible; also storing each $B$ as a dense matrix would not be possible.

Brian suggests calculating $A^{+}$ as $(A^{T}A)^{-1}A = A^{+}$, since $A$ has full column rank. We can rearrange this problem into $(A^{T}A)A^{+} = A$. This immediately suggests three approaches, all based on the normal equations of linear regression (i.e., least squares):

1) Use sparse Cholesky factorization of $A^{T}A$ to find a sparse matrix $R$ such that $R^{T}R = A^{T}A$; this step can be carried out using a call to CHOLMOD. Call $RA^{+} = X$. Solve the equation $R^{T}X = A$ for $X$ using a sparse LU factorization such as in UMFPACK; $R^{T}$ is lower triangular, so this solve should be quick, and only use symbolic factorization and some pivoting. Then solve $RA^{+} = X$ for $A^{+}$, again using a sparse LU factorization. $R$ is upper triangular, and it might be possible to reuse some of the information from the previous LU factorization in this solve.

2) Use a sparse QR factorization to find $R$ -- it's the same matrix as the one above. SuiteSparseQR is one implementation of sparse QR by Tim Davis, author of UMFPACK and CHOLMOD. This part saves you a matrix multiply (which doesn't really matter in the grand scheme of things), and you can either use the two sparse LU factorization approach above, or treat it as a least squares problem directly.

3) Use an SVD factorization of $A$, and the formula in this Wikipedia article.

Of these three approaches, the sparse QR decomposition and SVD approach are what I used in my thesis. In your case, I wouldn't necessarily recommend the SVD approach because I'm not aware of a good recent implementation of sparse SVD; if one is available, it's worth trying. Cholesky factorization is supposed to do a better job of exploiting sparsity than QR, and is probably fastest. QR is generally used in practice. SVD has the best numerical stability, followed by Cholesky and QR.

I assume in Brian's answer that when he suggests calculating $(A^{T}A)^{-1} = (R^{T}R)^{-1}$, he probably really means, "solve the system $(R^{T}R)A^{+} = A$". Solving a linear system by calculating the inverse explicitly and then multiplying is conventionally considered bad practice. Evidence suggests it's really not that bad for accuracy, but the inverse of a sparse matrix is frequently dense, so inverting the matrix is likely to require more storage and more operations than a sparse factorization. Furthermore, $R$ is better conditioned than $R^{T}R$, so recasting the linear system $(R^{T}R)A^{+} = A$ as two separate linear systems is at least as accurate as solving the single linear system; solving it as two separate systems is the conventional advice given in least squares problems.

Brian also rightly pointed out that there are approaches to calculate $A^{+}B$ directly using QR factorizations, and that these approaches would be unfavorable because they would require more computation and probably more memory.


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.