5
$\begingroup$

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)

Improve this question
$\endgroup$
5
  • 2
    $\begingroup$ 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. $\endgroup$
    – Kirill
    Oct 25, 2016 at 19:07
  • $\begingroup$ 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) $\endgroup$
    – Bananach
    Oct 25, 2016 at 19:14
  • $\begingroup$ Sorry, I misread that as "columns are close" (to equal), my mistake. $\endgroup$
    – Kirill
    Oct 25, 2016 at 19:16
  • $\begingroup$ Another thought: is the matrix multiplication in $A^*A$ multithreaded? $\endgroup$
    – Kirill
    Oct 26, 2016 at 2:02
  • $\begingroup$ I don't know. I am using Numpy's .dot() $\endgroup$
    – Bananach
    Oct 26, 2016 at 14:03

3 Answers 3

Reset to default
7
$\begingroup$

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.

Improve this answer
$\endgroup$
2
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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 $\endgroup$
    – Bananach
    Oct 26, 2016 at 16:09
  • $\begingroup$ You can even use the LSQR algorithm, which is mathematically equivalent to CG on the normal equations, but numerically more stable $\endgroup$
    – Richard Zhang
    Nov 18, 2016 at 9:34
  • $\begingroup$ Yes, but with cond(A*A)=5, I wouldn't expect much difference between LSQR and CG. $\endgroup$
    – wim
    Nov 18, 2016 at 10:15
1
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.