Sparse linear poorly constrained least-squares problem

Question

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.

Answer 1

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.