# 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:** 20\
**Page:** 2

<div class="post-metadata">

**Author:** ![Paul\_Soderlind](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paul_soderlind/32/1753_2.png) [@Paul\_Soderlind](https://discourse.julialang.org/u/Paul_Soderlind)\
**Post date:** [November 20, 2019, 12:35pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/22 "2019-11-20T12:35:59Z")

</div>

more details on this are found in, for instance, Golub&van Loan (3rd ed). I like the very last sentence of section 5.3: “At the very minimum, this discussion should convince you how difficult it can be to choose the ‘right’ algorithm!”

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [November 21, 2019, 3:27pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/23 "2019-11-21T15:27:12Z")

</div>

Could you point out which those packages are, please?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [November 21, 2019, 3:37pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/24 "2019-11-21T15:37:14Z")

</div>

I have only used

> **[GitHub - JuliaStats/GLM.jl: Generalized linear models in Julia](https://github.com/JuliaStats/GLM.jl)**
>
> Generalized linear models in Julia. Contribute to JuliaStats/GLM.jl development by creating an account on GitHub.

and

> **[GitHub - joshday/OnlineStats.jl: ⚡ Single-pass algorithms for statistics](https://github.com/joshday/OnlineStats.jl)**
>
> ⚡ Single-pass algorithms for statistics. Contribute to joshday/OnlineStats.jl development by creating an account on GitHub.

in practice.

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [November 21, 2019, 3:47pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/25 "2019-11-21T15:47:48Z")

</div>

I use polynomials.jl and lsqfit.jl. Not caring about speed I’ve never benchmarked.

Example:

```julia
using Polynomials
p = polyfit(x,y,1) # fit first order polynomial
coeffs(p)

```

---

<div class="post-metadata">

**Author:** ![dmbates](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dmbates/32/44_2.png) [@dmbates](https://discourse.julialang.org/u/dmbates)\
**Post date:** [November 21, 2019, 4:24pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/26 "2019-11-21T16:24:08Z")

</div>

I describe a number of different ways of performing a simple linear regression in

[https://github.com/dmbates/CopenhagenEcon/blob/master/jmd/03-LinearAlgebra.jmd](https://github.com/dmbates/CopenhagenEcon/blob/master/jmd/03-LinearAlgebra.jmd)

As shown in the `README.md` file for the repository an effective way of running the .jmd file is to use `Weave.convert_doc` to create a Jupyter notebook and run that. (By the way, this notebook also shows using the R `ggplot2` graphics package from Julia through `RCall`.).

It is difficult to benchmark a simple calculation like this but I suspect that the method of augmenting the model matrix with the response vector and using a Cholesky factorization will be competitive. It has the advantage of also providing the sum of squared residuals without needing to evaluate the residuals.

---

<div class="post-metadata">

**Author:** ![markushhh](https://avatars.discourse-cdn.com/v4/letter/m/f1d935/32.png) [@markushhh](https://discourse.julialang.org/u/markushhh)\
**Post date:** [November 25, 2019, 8:50am UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/27 "2019-11-25T08:50:11Z")

</div>

what about?

```julia
using Random

N = 10000
x = rand(N)
X = [ones(N) x]
y = 10 .+ x .* 0.3

function linreg1(y, X)
    β_hat = (X' * X) \ X' * y
    return(β_hat)
end

function linreg2(y, X)
    β_hat = X \ y
    return(β_hat)
end

using GLM
GLM.fit(LinearModel, X, y, true)

using DataFrames, GLM
data = DataFrame(X = x, Y = y)
lm(@formula(Y ~ X), data)

using BenchmarkTools
@benchmark linreg1(y, X)
@benchmark linreg2(y, X)
@benchmark GLM.fit(LinearModel, X, y, true)
@benchmark lm(@formula(Y ~ X), data)

```

and the results are

lingerg1:

```julia
BenchmarkTools.Trial: 
  memory estimate: 156.78 KiB
  allocs estimate: 8
  --------------
  minimum time: 213.900 μs (0.00% GC)
  median time: 289.351 μs (0.00% GC)
  mean time: 552.780 μs (4.66% GC)
  maximum time: 76.192 ms (0.00% GC)
  --------------
  samples: 8836
  evals/sample: 1

```

linreg2:

```julia
BenchmarkTools.Trial: 
  memory estimate: 303.69 KiB
  allocs estimate: 45
  --------------
  minimum time: 118.300 μs (0.00% GC)
  median time: 170.300 μs (0.00% GC)
  mean time: 344.598 μs (17.20% GC)
  maximum time: 57.911 ms (0.00% GC)
  --------------
  samples: 10000
  evals/sample: 1

```

fit():

```julia
BenchmarkTools.Trial: 
  memory estimate: 470.45 KiB
  allocs estimate: 26
  --------------
  minimum time: 133.700 μs (0.00% GC)
  median time: 189.299 μs (0.00% GC)
  mean time: 428.768 μs (17.15% GC)
  maximum time: 45.171 ms (0.00% GC)
  --------------
  samples: 10000
  evals/sample: 1

```

lm():

```julia
BenchmarkTools.Trial: 
  memory estimate: 1.08 MiB
  allocs estimate: 186
  --------------
  minimum time: 546.600 μs (0.00% GC)
  median time: 839.500 μs (0.00% GC)
  mean time: 1.917 ms (10.94% GC)
  maximum time: 127.190 ms (0.00% GC)
  --------------
  samples: 2588
  evals/sample: 1

```

---

<div class="post-metadata">

**Author:** ![dmbates](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dmbates/32/44_2.png) [@dmbates](https://discourse.julialang.org/u/dmbates)\
**Post date:** [November 25, 2019, 7:34pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/28 "2019-11-25T19:34:47Z")

</div>

First you should decide what you want the function to be. Here you are assuming that the model matrix `X` has been externally created. I believe the original suggestion was to have two vectors, `x` and `y`, as the arguments. In that case it is best to go old school and create `X'X` and `X'y` without ever creating X. I would say that this knowledge was from back in the old days except that, astonishingly, you could probably get these formulas from a recent introductory stats text.

```julia
julia> linreg1(x::AbstractVector{T}, y::AbstractVector{T}) where {T<:AbstractFloat} = [ones(length(x)) x]\y
linreg1 (generic function with 1 method)

julia> linreg1(x, y)
2-element Array{Float64,1}:
 0.4986028246024461
 0.29993048137652306

julia> @benchmark linreg1($x, $y)
BenchmarkTools.Trial: 
  memory estimate: 538.31 KiB
  allocs estimate: 52
  --------------
  minimum time: 135.517 μs (0.00% GC)
  median time: 140.943 μs (0.00% GC)
  mean time: 150.165 μs (3.99% GC)
  maximum time: 832.710 μs (78.93% GC)
  --------------
  samples: 10000
  evals/sample: 1

julia> function linreg2(x::AbstractVector{T}, y::AbstractVector{T}) where {T<:AbstractFloat}
           X = [ones(length(x)) x]
           ldiv!(cholesky!(Symmetric(X'X, :U)), X'y)
       end
linreg2 (generic function with 1 method)

julia> linreg2(x, y)
2-element Array{Float64,1}:
 0.4986028246024445
 0.29993048137652634

julia> @benchmark linreg2($x, $y)
BenchmarkTools.Trial: 
  memory estimate: 234.94 KiB
  allocs estimate: 13
  --------------
  minimum time: 67.063 μs (0.00% GC)
  median time: 70.251 μs (0.00% GC)
  mean time: 74.033 μs (3.30% GC)
  maximum time: 652.841 μs (86.41% GC)
  --------------
  samples: 10000
  evals/sample: 1

julia> function linreg3(x::AbstractVector{T}, y::AbstractVector{T}) where {T<:AbstractFloat}
           (N = length(x)) == length(y) || throw(DimensionMismatch())
           ldiv!(cholesky!(Symmetric([T(N) sum(x); zero(T) sum(abs2, x)], :U)), [sum(y), dot(x, y)])
       end
linreg3 (generic function with 1 method)

julia> linreg3(x, y)
2-element Array{Float64,1}:
 0.49860282460244537
 0.299930481376524

julia> @benchmark linreg3($x, $y)
BenchmarkTools.Trial: 
  memory estimate: 368 bytes
  allocs estimate: 9
  --------------
  minimum time: 10.041 μs (0.00% GC)
  median time: 10.260 μs (0.00% GC)
  mean time: 10.473 μs (0.00% GC)
  maximum time: 105.643 μs (0.00% GC)
  --------------
  samples: 10000
  evals/sample: 1

```

If you look closely the big differences are in the amount of storage that is allocated.

---

<div class="post-metadata">

**Author:** ![mkborregaard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkborregaard/32/556_2.png) [@mkborregaard](https://discourse.julialang.org/u/mkborregaard)\
**Post date:** [November 25, 2019, 7:56pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/29 "2019-11-25T19:56:31Z")

</div>

Sounds like a pretty strong case for having an efficient `linreg` function in the Statistics stdlib.

---

<div class="post-metadata">

**Author:** ![markushhh](https://avatars.discourse-cdn.com/v4/letter/m/f1d935/32.png) [@markushhh](https://discourse.julialang.org/u/markushhh)\
**Post date:** [November 25, 2019, 8:53pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/30 "2019-11-25T20:53:51Z")

</div>

Wouldn’t you always have to create X externally in a for loop?

---

<div class="post-metadata">

**Author:** ![dmbates](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dmbates/32/44_2.png) [@dmbates](https://discourse.julialang.org/u/dmbates)\
**Post date:** [November 25, 2019, 9:21pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/31 "2019-11-25T21:21:54Z")

</div>

I’m not sure I understand the question. If you have the vector `x` and you want `X'X` and `X'y` you don’t need to allocate and populate the N by 2 matrix `X`. That’s what the section in the code with

```julia
Symmetric([T(N) sum(x); zero(T) sum(abs2, x)], :U)

```

and

```julia
[sum(y); dot(x, y)]

```

is doing.

Once upon a time (I took my first statistics course in 1968) we needed to learn this kind of short cut because we were doing the calculations by hand. But these days you would need to be doing a very large number of such calculations before saving a few hundred microseconds would be of importance.

And, from the viewpoint of a statistician, there is a lot of information in addition to the coefficient estimates that you would want.

---

<div class="post-metadata">

**Author:** ![markushhh](https://avatars.discourse-cdn.com/v4/letter/m/f1d935/32.png) [@markushhh](https://discourse.julialang.org/u/markushhh)\
**Post date:** [November 25, 2019, 10:07pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/32 "2019-11-25T22:07:11Z")

</div>

I just learned so much from just those two lines!  
Thank you! I didn’t even know there was such a thing…

---

<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.

---

<div class="post-metadata">

**Author:** ![brett\_knoss](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/brett_knoss/32/13050_2.png) [@brett\_knoss](https://discourse.julialang.org/u/brett_knoss)\
**Post date:** [March 17, 2020, 4:51am UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/34 "2020-03-17T04:51:57Z")

</div>

does x’ mean the derivative of x? and is x independent, while y is dependent?

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [March 17, 2020, 5:38am UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/35 "2020-03-17T05:38:37Z")

</div>

x’ means transpose of x in this context.

If we were talking about odes it’d probably mean x’(t)

---

<div class="post-metadata">

**Author:** ![brett\_knoss](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/brett_knoss/32/13050_2.png) [@brett\_knoss](https://discourse.julialang.org/u/brett_knoss)\
**Post date:** [March 17, 2020, 6:11am UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/36 "2020-03-17T06:11:12Z")

</div>

that makes sense.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [August 12, 2020, 12:35am UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/37 "2020-08-12T00:35:00Z")

</div>

Sorry to bring this discussion back. But I think that many times the user is just searching for a function which given two vectors x and y, returns the coefficients of the linear fit and the Pearson correlation coefficient, and perhaps a third vector with the fitted Y values in the same x positions, such that plotting a linear regression over some data is direct. Something like:

```julia
x = [...]
y = [...]
L = linreg(x,y)
scatter(x,y)
plot!(x,L.Y,label="R^2=$(L.R2)")

```

Actually I think one such alternative would be nice to exist in the Plots package perhaps.

I mean, that is what I was searching for now, and it was easier to just plot what I wanted somewhere else at the end. A suite of the most simple regressions (linear, square, log, perhaps) would be nice accompanying Plots for a fast visual inspection of data.

Ps. I know this is three-lines of code anyway. But I assure you that most people in the world doing linear regressions have no idea of how they are done.

Ps2. Noticed now that this discussion is the “Performance” section, so that was not the case here. But this is one of the threads one reaches searching for a simple way to plot a linear regression.

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [August 12, 2020, 2:19am UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/38 "2020-08-12T02:19:39Z")

</div>

You can do exactly that already. I will agree that it isn’t as “in your face” or easy to find as it could be though.

```julia
using GLM
x = [...]
y = [...]
reg = lm(x,y)
scatter(x,y)
plot!(x,predict(reg),label="R^2=$(r2(reg))")

```

also `Plots.jl` or `StatsPlots.jl` (I can’t remember which one) provides a direct option in the `scatter` command to put the regression line on the plot for you.

---

<div class="post-metadata">

**Author:** ![BeastyBlacksmith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/beastyblacksmith/32/4741_2.png) [@BeastyBlacksmith](https://discourse.julialang.org/u/BeastyBlacksmith)\
**Post date:** [August 12, 2020, 3:11pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/39 "2020-08-12T15:11:18Z")

</div>

The keyword is `smooth = true`. However, it doesn’t provide any analytics beyond the curve itself.

---

<div class="post-metadata">

**Author:** ![mkborregaard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkborregaard/32/556_2.png) [@mkborregaard](https://discourse.julialang.org/u/mkborregaard)\
**Post date:** [August 12, 2020, 6:17pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/40 "2020-08-12T18:17:54Z")

</div>

see also [https://github.com/JuliaPlots/StatsPlots.jl/pull/293](https://github.com/JuliaPlots/StatsPlots.jl/pull/293)

---

<div class="post-metadata">

**Author:** ![piever](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/piever/32/1815_2.png) [@piever](https://discourse.julialang.org/u/piever)\
**Post date:** [August 12, 2020, 7:41pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/41 "2020-08-12T19:41:06Z")

</div>

This (linear regression with confidence interval) is also implemented in AlgebraOfGraphics with the [linear](http://juliaplots.org/AlgebraOfGraphics.jl/stable/analyses/#AlgebraOfGraphics.linear) analysis.

[Previous page](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232.md?page=1)

[Next page](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232.md?page=3)