Least Squares: Numerically, is solving normal equations okay for nice matrices?

I have to solve a least squares problem: $$x=\arg \min\|Ax-b\|$$ where $A$ is a $m\times n$ matrix, $m>n$, and $b\in\mathbb{R}^m$.

I always thought that doing this via QR factorization is better than solving the normal equations $$A^*Ax=A^*b$$ directly.

However, I am now in the situation where the columns of $A$ are close to orthogonal and my code is significantly faster (factor 10) when using normal equation. Just to be sure I am doing nothing wrong, is this a reasonable result?

($A$ has size $20000\times 2000$. Computing the Gramian matrix takes longer than solving the normal equation)

• Solving normal equations should take $\sim mn^2+\frac13 n^3$ flops (forming $A^*A$, then solving it with Cholesky factorization), and doing QR should take $2mn^2$ flops (for example: seas.ucla.edu/~vandenbe/103/lectures/qr.pdf), which is why a factor of 10 speedup is not expected. – Kirill Oct 25 '16 at 19:07
• I include the time to form $A^*A$ into the cost of "the normal equations". I am confused: the columns are close to orthogonal to each other, thus $A^*A$ has a very good condition. The condition numbers are at most 5. I am using python linalg.solve for the normal equations and scipy.lstsq for the "direct" appraoch (as far as I know this uses SVD by standard, but I also tried all the other LINPACK options that scipy offers) – Bananach Oct 25 '16 at 19:14
• Sorry, I misread that as "columns are close" (to equal), my mistake. – Kirill Oct 25 '16 at 19:16
• Another thought: is the matrix multiplication in $A^*A$ multithreaded? – Kirill Oct 26 '16 at 2:02
• I don't know. I am using Numpy's .dot() – Bananach Oct 26 '16 at 14:03

For solving the least squares problem in general as principal methods there are (matrix $A$ with full rank):

• solve the system of normal equations $A^{T}Ax = A^{T}b$
• use QR factorization
• use SVD decomposition

Generally speaking the QR factorization is a method with a good balancing between accuracy and computational cost.

Using the normal equations is possible, but it is often avoided due to condition number. The matrix of the system is $A^{T}A$, and it has got $$K_{A^{T}A} = (K_{A})^{2}$$ where $K$ is the condition number. So if your matrix has got a good $K$, the normal equations are OK.

Note that the matrix $A^{T}A$ is symmetric and definite positive (the eigenvalues are $\sigma^{2}_{i} > 0$ with $\sigma_{i}$ singular value of $A$) so to solve directly the system you can use the Cholesky factorization.

In a comment you state that the condition number of the system of the normal equentions is as small as $\kappa=5$. For such a well condtioned problem it might be a good idea to try the conjugate gradient iterative method. This method should converge fast, according to the well known conjugate gradient error estimate $$||x_*-x_m||_A \le 2 \left[ \frac{\sqrt \kappa -1}{\sqrt \kappa +1} \right]^m ||x_*-x_0||_A \, ,$$ where $A$ is not your $A$, but the matrix $A^*A$ from the normal equations (the one with condition number 5). Note that $\left[ \frac{\sqrt 5 -1}{\sqrt 5 +1} \right]^{40}=$1.9e-17 .

An advantage of the conjugate gradient method is, that there is no need to compute $A^* A$ explicitly, only matrix vector products of the form $q_i=A^* Ap_i$ are needed. These are computed as $q_i=A^* (Ap_i)$ .

See also Saad's book: Iterative methods for sparse linear systems, chapter 6 and 8.

• That's a great observation. However, CG ends up not being faster for me, I guess because pythons sparse module is not as optimized as the usual solver – Bananach Oct 26 '16 at 16:09
• You can even use the LSQR algorithm, which is mathematically equivalent to CG on the normal equations, but numerically more stable – Richard Zhang Nov 18 '16 at 9:34
• Yes, but with cond(A*A)=5, I wouldn't expect much difference between LSQR and CG. – wim Nov 18 '16 at 10:15

The other answers already give good advice, so I want to offer a guess as to the cause of the unexpected speedup.

Solving normal equations should take $\sim mn^2+13n^3$ flops (forming $A^*A$, then solving it with Cholesky factorization), and doing QR should take $2mn^2$ flops (for example: https://seas.ucla.edu/~vandenbe/103/lectures/qr.pdf), which makes us expect that normal equations would be 2 times faster.

Since you say that you use numpy.dot, and numpy.dot is automatically parallelized (https://scipy.github.io/old-wiki/pages/ParallelProgramming), it seems possible that if QR is not parallelized, you would see a speedup of $2\times \text{number-of-threads}$, so a speedup of 10 might then be reasonable. But this multithreading speedup is usually not mentioned because the issue with squaring condition number is more important.

To check, you'd need to find out what BLAS/LAPACK library your numpy is using (numpy.show_config), and, to make sure, set the appropriate environment variable (like OMP_NUM_THREADS) to 1 to check if it's really multithreading that causes this. There are already a number of questions on StackOverflow about how to do this.