# Efficient way of doing linear regression

**URL:** https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232
**Category:** Performance
**Tags:** regression
**Created:** [November 18, 2019, 5:21pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232 "2019-11-18T17:21:15Z")
**Posts on this page:** 1
**Showing post:** 33

<div class="post-metadata">

### Author: ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [January 14, 2020, 8:17pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/33 "2020-01-14T20:17:05Z")

</div>

> [@Paul\_Soderlind](#):
>
> I have recently noticed that `b=(x'x)\(x'y)` is often twice as fast as `b = x\y` , for instance, when `x is 1000x10` and `y is 1000x1` . So far, I have no good explanation for it.

Sorry to come late to this thread, but I’m surprised that no one seems to have given the correct explanation here for the performance differences between various methods to solve \min\_x \Vert Ax-b\Vert\_2.

- `x̂ = (A'A)\(A'b)` is the “normal equations” approach. (You could also use `cholesky(Hermitian(A'A))\(A'b)` to explicitly take advantage of the fact that `A'A` is SPD). This (and variations thereof) is the **fastest** approach, but it is also the **least accurate**. In particular, it squares the condition number of `A`, so it _doubles_ the number of digits you lose to roundoff errors and similar. You should only use this method if you know that you have a well-conditioned `A`.

- `x̂ = A \ b` is equivalent to `qr(A, Val{true}()) \ b` — it uses a pivoted QR factorization. This is slower than the normal-equations approach, but it is **much more accurate** for badly conditioned matrices because it doesn’t square the condition number.

- `x̂ = pinv(A) * b` uses the SVD of `A` to apply the pseudo-inverse. This is the **slowest** method, but it gives you the **most control over accuracy** for ill-conditioned matrices because you can specify a tolerance in `pinv` to regularize the problem by dropping noisy singular values (though a very similar result can also be accomplished by [Tikhonov regularization](https://en.wikipedia.org/wiki/Ridge_regression)). _Update_: using the pseudo-inverse by an explicit call to `pinv` may not be backwards-stable [(Liu & Barnett, 2016, and references therein)](https://doi.org/10.1016/j.jcp.2016.08.011) — it is better to use `x̂ = svd(A) \ b`, which computes \hat{x} = \hat{V} (\hat{\Sigma}^+ (\hat{U}^\* b)), though you [currently need](https://github.com/JuliaLang/LinearAlgebra.jl/pull/1387) to do the latter manually if you want to pass a tolerance.

QR is the default choice used for `A \ b` (in Julia and many other systems) because it is reliable without being too slow.

---

_[View the full topic](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232)._