# How can we solve the normal equations with limited memory?

I was asked this open ended question in an interview once:

How would you find a solution to the normal equations with limited memory?

Unlike Solving sparse least squares system with limited memory, the matrix $$A$$ is not necessarily sparse. It could be dense.

The question was asked in a machine learning / stats context, so I gave the answer that instead of solving it exactly (e.g., with QR), you can solve it using stochastic or batch gradient descent and obtain an approximate solution.

In a linear algebra and scientific computing context, is there any way to solve a dense least squares system with limited memory?

• What do you mean by 'limited memory' here exactly? If $A$ is $m\times n$, do $mn$ real numbers fit in RAM? Do $\min(m,n)^2$? Do $n$? Feb 28 at 18:45
• @FedericoPoloni That was sort of open ended during the interview. But limited definitely meant that we could not store all of $A$ in RAM. But I think it would be reasonable to assume that a single column or row could fit in RAM. Feb 28 at 19:12
• There are so-called limited memory Krylov subspace methods which can be used as direct solvers if you iterate them $n$ iterations where $n$ is the number of equations. Of course, I neglected the round-off issues due to computer arithmetic, so you would stop earlier than that. Also they might be looking for an answer like L-BFGS but I don't know how you would apply that to a general normal equation. Feb 28 at 21:34
• There are situations in which $A$ is sparse but $A^{T}A$ is dense. In those situations, iterative methods for least squares such as LSQR can be used to solve the least squares problem without requiring additional storage. Feb 28 at 23:48
• @user5965026: It’s a good question too. Pretty common when things get serious not to have enough memory to solve the system of interest. Or if there is, there eventually won’t be. 🙂 Mar 1 at 0:03

That might also have been a trick question. Let's say you want to solve the normal equations for $$Ax=b$$, i.e., $$(A^T A) x = A^T b$$. Let's assume for a moment that the questioner meant that $$A$$ is actually already stored in memory, so we know that that much memory is already available. Let's also assume that $$A$$ is tall and narrow (more specifically, has fewer columns than rows) -- because otherwise solving the normal equations.

Then, if the matrix is dense, then so is $$B=A^TA$$ and furthermore $$B$$ is smaller than $$A$$ and so at the cost of at most a factor of two in memory, $$B$$ can be stored if $$A$$ can be stored.

If $$A$$ is sparse, then it could be that $$B$$ is substantially more expensive to store, though that depends on the sparsity pattern of $$A$$. If that were the case, then we can solve the normal equations with Conjugate Gradients (CG) because $$B$$ is a symmetric and (hopefully) positive definite matrix. In CG, all you need is the ability to multiply a vector by the matrix $$B$$, which can be implemented with minimal additional memory by writing $$Bv=(A^TA)v=A^T(Av)$$, i.e., you just have to do two subsequent matrix-vector products -- it is never necessary to form the product matrix $$B$$.

• Ahh good catch. I never thought about that $B$ would necessarily take less memory to store since it'd be a smaller matrix than $A$ (assuming $A$ isn't sparse). Couldn't you also use CG if $A$ isn't sparse? Not sure if you were implying that CG only applies if $A$ is sparse or not. Mar 1 at 0:42
• Yes, of course you can use CG for any kind of $A$. But if $A$ is not sparse, it's probably cheaper to either form $B$ right away, or at least do something like a $QR$ decomposition of $A$ for use in solving the normal equations. Mar 1 at 3:49
• Ah I see. When you say "cheaper" in this context, do you mean computationally cheaper in terms of time, but not in terms of RAM usage? Mar 1 at 20:54
• Maybe both. These days, memory and CPU time are highly correlated. You can do so many more computations whereas memory transfer is the bottleneck. (You could do ~25 FP operations in the time necessary to load the data for just one of these operations.) So if $A$ is not sparse but large, it's expensive to do matrix-vector products because you have to get so much memory onto the chip. You might be better off not doing that many times over, but just forming $B$ once (now smaller -- might fit into the cache) or using an efficient QR decomposition that works on matrix blocks. Mar 2 at 2:28

So called “matrix free” methods, relying primarily on the ability to perform multiplication of the matrix by a vector, lend themselves nicely to iterative techniques such as GMRES. The matrix itself might be on disk, but portions retrieved selectively to compute the necessary matrix-vector products.

• But in this situation, the entire matrix still has to be stored somewhere right? Feb 28 at 17:14
• Not really. Either in memory (most convenient), or stored on disk (a bit less convenient), or able to have portions, for example blocks, generated on demand, or having some way of computing $Ax$. As an example you can set up a finite differences solution for the steady state heat distribution problem and solve it (say Gauss-Seidel) without ever explicitly creating a matrix. Not the same as the least squares problem you have, but for illustration. Feb 28 at 17:32
• I see. Why is this called a "matrix free" method if it still involves the matrix? Feb 28 at 17:38
• Probably worthwhile having a look at en.wikipedia.org/wiki/… Feb 28 at 17:40
• If $A^T A$ is positive definite, it is very simple to implement a least squares conjugate gradient solver. It will require two matrix vector products per iteration, but no more memory than is required to store $A$ Feb 28 at 23:26

Among approximate techniques:

• Gradient descent should do a decent job, given the limitations.
• Randomized SVD is an effective technique. It is quick to implement, and there are ready-to-use error bounds (see e.g. https://doi.org/10.1137/090771806 for a review of that area).

Among techniques that produce an exact solution:

• There is research on out-of-core QR and SVD computation; essentially the main idea for QR is computing the QR of thinner 'panels' of $$A$$ beforehand and then merging them (but then of course the details get very technical).
• If your problem is well-conditioned, and if $$\min(m,n)^2$$ fits in memory, you can just use the normal equations solution $$x=(A^TA)^{-1}(A^Tb)$$, which you can compute with a single pass on your out-of-core data. This is by far the simplest solution in terms of implementation (followed by randomized SVD), but its stability is only conditional and problem-dependent.
• What is the significance of the square on the $\min(m,n)$? And also, for least squares isn't the smaller dimension always the number of columns of $A$, otherwise $A^TA$ would necessarily be singular since $A$ would not be full rank if the number of columns is greater than the number of rows. Mar 1 at 0:39
• @user5965026 Yes, in this formulation $A$ is a tall thin matrix and $\min(m,n)=n$. However, people in machine learning often work with the transposed problem, so I preferred to use $\min(m,n)$ above in a comment for clarity, and then I copied it here below as well without thinking that in this post I am also using that explicit formula for $x$ that works only in one case. Anyhow, that $\min(m,n)^2$ is the size of $A^TA$; if that does not fit in memory it is trickier to assemble it on disk. Mar 1 at 7:11