I have a somewhat simple linear problem.

I have data $D$ (a vector with a few million elements), the parameter vector $X$ (a couple of thousands elements) and the design matrix $A$ which is extremely sparse.

To find $X$ I am minimizing the $L_2$ norm $||D-A X||_2$ using and iterative sparse chi-square solver (LSMR). The problem however is that the $A$ operator is essentially a convolution operator, therefore the deconvolution is not well determined and I am getting various ringing artefacts caused by the deconvolution (the covariance matrix of the solution is essentially degenerate/not positive-definite).

Because of this I actually do not care about $X$ but a reduced parameter space $Y=BX$ where $Y$ is a smaller dimensionality vector than $X$ ($B$ is some matrix). I also know that $Y$ should be much better constrained by the data. And I am trying to see if there is a way of solving for $Y$ without first solving for $X$. ($B$ is also sparse and I do know it). I can still solve first for X then transform to Y, but because the linear system is poorly conditioned, the number of iterations needed to get X right is too large.


I don't know any way around computing $X$, but I can give you a few suggestions for computing $X$: sparse QR and sparse SVD.

If A were dense and a reasonable size, the QR factorization would be the standard approach. The analytic solution to your problem is $X = (A^T A)^{-1} A^T D$. If you factor $A = Q R$, then the solution becomes:

\begin{array}{} X & = & (A^T A)^{-1} A^T D \\ & = & ((Q R)^T (Q R))^{-1} (QR)^T D \\ & = & (R^T Q^T Q R)^{-1} R^T Q^T D \\ & = & (R^T R)^{-1} R^T Q^T D \\ & = & R^{-1} R^{-T} R^T Q^T D \\ & = & R^{-1} Q^T D \end{array}

The problem with using the QR factorization on a large, sparse system is fill-in. If $A$ is $m \times n$ then $Q$ is $m \times n$ and dense. In your case, $m \approx 10^6$, so forming $Q$ is out of the question.

This is where Givens rotations come in to play. I won't go in to tremendous detail about Givens rotations, but here's the general idea. You basically transform $A$ in to $R$ by applying a series of Givens rotations. Each Givens rotation will zero out exactly one element of $A$. You only need to zero out the elements below the diagonal. Unfortunately, in your case, almost all of the elements are below the diagonal. So if your matrix has $nnz(A)$ nonzeros, you'll need close to $nnz(A)$ Givens rotations compute R. On the plus side, Givens rotations are very simple and easy to apply. Still, this may be intractable for your problem; I can't say.

But here's where the magic comes in to play: You NEVER need to explicitly form $Q$. Instead, you can represent $Q$ as the product of all the Givens rotations: $Q = G_1 G_2 ... G_k$. And thus $Q^T D = G_k^T G_{k-1}^T ... G_1^T D$. To be sure, this approach requires no small amount of bookkeeping. And it might be computationally intractable. But it is a sure-fire solution to your problem.

The other approach is a sparse SVD. Instead of factoring $A = QR$, we now factor $A = U \Sigma V^T$. Then we have:

\begin{array}{} X & = & (A^T A)^{-1} A^T D \\ & = & ((U \Sigma V^T)^T (U \Sigma V^T))^{-1} (U \Sigma V^T)^T D \\ & = & (V \Sigma U^T U \Sigma V^T)^{-1} V \Sigma U^T D \\ & = & (V \Sigma^2 V^T)^{-1} V \Sigma U^T D \\ & = & V \Sigma^{-2} V^T V \Sigma U^T D \\ & = & V^T \Sigma^{-1} U^T D \end{array}

Unfortunately, we have the same problem as before. $U$ is $m \times n$ and dense, so we can't even store it, let alone compute it! But there's still something we can do. Suppose that, instead of computing the full SVD, we compute only the $r$ largest singular values and corresponding singular vectors (there are pretty good algorithms for doing this). The result is the optimal rank $r$ approximation $A \approx U_r \Sigma_r V_r^T$ Now we can approximate the $X$ using the reduced SVD.

$$X \approx V_r^T \Sigma_r^{-1} U_r^T D$$

Of course, this is not an exact answer. However, you already said that $A$ is not full-rank. So even if you can only compute a few singular values of $A$, that might be enough to get a good approximation. This is the case when $A$ is dominated by a handful of very large singular values. So the approximation might turn out to be quite good. Might be worth a try.

  • $\begingroup$ Golub & Van Loan call the "The Sparse QR Challenge" and it is discussed in pp. 606-608 $\endgroup$ – Reinstate Monica Aug 29 at 21:24
  • $\begingroup$ Oh yeah, I completely forgot to mention column reordering. You can also reorder the columns of A before factorizing it to help prevent fill-in. This is more of an art than a science. There are several algorithms to estimate the optimal reordering and they're all pretty much graph partitioning algorithms: nested dissection, approximate minimum degree, Dulmage-Mendelsohn and Cuthill-McKee are some of the better-known algorithms. $\endgroup$ – Bill Woessner Aug 30 at 3:29
  • $\begingroup$ On suggestion from @Sycorax I checked the new edition of Golub & Van Loan. They mention another approach called the semi-normal equations. This allows you to do away with the $Q$ matrix completely at the cost of an extra triangular solve. Your solution becomes $X = R^{-1} R^{-T} A^T D$. This is not as numerically stable as the traditional QR approach, but you can use an iterative refinement to improve the precision of the solution. $\endgroup$ – Bill Woessner Aug 30 at 15:53
  • $\begingroup$ Thank you for the suggestions. I didn't have time over last few days. But I will seen if I can implement some of them. $\endgroup$ – sega_sai Sep 3 at 20:08

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.