Solving sparse least squares system with limited memory

This was a question on a past final that we can't figure out. Take the least squares system

$$\min_x ||Ax-b||_2\, ,$$

where $A\in\mathbb{R}^{mxn}$, $m<n$, and A is full rank. A has $\mathcal{O}(n)$ non-zero entries. Assume you only have $\mathcal{O}(n)$ computer memory.

1. Suggest an algorithm for approximating the solution to the least squares problem
2. What is the expected number of iterations and overall computational cost?
3. What is the expected error in the solution x? (Provide an appropriate definition of error for this problem)

We've only learned a little bit about iterative methods, just like Rayleigh Quotient, Inverse Iteration, Power Iteration, and Conjugate Gradient. None of these seem to help. Any ideas would be much appreciated.

Edit: I'm thinking now, what if we take a random mm columns so that a square matrix can be formed, and then use CG on that problem? And then do this with various random columns to find average solutions?

• hint: the $A$ matrix is sparse (it has ${\cal O} (n) \ll n \times m$ nonzero entries, i.e., a constant number of nonzeros per column). you'll want to exploit that sparsity in your iterative least squares algorithm. – GoHokies Dec 17 '17 at 15:48
• Yes, that's why I'm reading up on algorithms suited to sparse matrices. Currently I'm looking at CG since we discussed it in class, but that's only for symmetric pos-def. matrices. I'm thinking now about taking $A^TA$ and using that (solving the normal equations I guess, $A^TAx = A^Tb$) – wobertson Dec 17 '17 at 15:54
• yes, but keep in mind that you don't have enough memory to store all elements of $A^T A$. what you will need with CG is a matrix-vector product routine that implements $(A^T A) v$ for a given vector $v$ and requires only ${\cal O}(n)$ memory space. can you see how such a routine could be implemented? – GoHokies Dec 17 '17 at 15:59
• Actually $A^TA$ is non-singular because $m<n$, but $AA^T$ could work, and we could first store the elements of $x = A^Ty$ and use these for $Ax=b$, is this what you're talking about? – wobertson Dec 17 '17 at 16:51
• @GoHokies -- I don't think that $A^T A$ has $O(n^2)$ entries. It should still be a sparse matrix. – Wolfgang Bangerth Dec 17 '17 at 21:15

Your call for $O(N)$ immediately triggers the multigrid methods in my mind. Multigrid techniques give algorithms that solve sparse linear systems $Ax = b$ of $N$ unknowns with $O(N)$ computational and storage cost for large classes of problems. Thus it is also very suitable for your problem.

The steps of the general method follow:

• Compute residuals
• Restrict the residual to coarse grid
• Solve the coarse system
• Prolongate the error to a finer grid
• Update the approximation

All of these schemes are quite well studied and many approaches exist. Multigrids are also an active research area. Check this out for an intro.

• This is an exercise -- it is very likely the instructor had something else in mind. – Federico Poloni Jan 3 '18 at 21:48