# 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:** 1
### Author:  [@Johncowk](https://discourse.julialang.org/u/Johncowk)
#### Post date: [November 18, 2019, 5:21pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/1 "2019-11-18T17:21:15Z")
Hello,
I need an efficient way to performs linear regression because I have to fit several segments in a big for loops which takes time to achieve. In the past, their was a `linreg` function which has been depreciated. Currently I am using `polyfit` to do a polynomial fit of order 1, for each segment.
Is there a more efficient way to do this ?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [November 18, 2019, 5:31pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/2 "2019-11-18T17:31:32Z")
The `\` operation is for linear regression. Just build the matrix.
---
### Author:  [@Johncowk](https://discourse.julialang.org/u/Johncowk)
#### Post date: [November 18, 2019, 5:33pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/3 "2019-11-18T17:33:23Z")
Oh I see. I’ll try it thanks !
---
### Author:  [@Johncowk](https://discourse.julialang.org/u/Johncowk)
#### Post date: [November 18, 2019, 8:02pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/4 "2019-11-18T20:02:44Z")
Ok so i didn’t understand this as good as i thought i did. Considering that I have two 1D arrays and I want to find the best linear curve to find one element of y given the corresponding element of x, how do I do it with `\` ?
---
### Author:  [@fipelle](https://discourse.julialang.org/u/fipelle)
#### Post date: [November 18, 2019, 8:14pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/5 "2019-11-18T20:14:18Z")
You can estimate the coefficients via OLS: `β̂ = (x'*x)\x'*y`.
The estimated regression line is `ŷ=x*β̂`.
---
### Author:  [@MatFi](https://discourse.julialang.org/u/MatFi)
#### Post date: [November 18, 2019, 8:32pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/6 "2019-11-18T20:32:53Z")
OK say you have the datapoints (x1, y1) (x2, y2) and (x3, y3).
in this case you could simply describe it via `M*v=d` with :
```julia
M =[1 x1
1 x2
1 x3]
v= [a
b]
d=[y1
y2
y3]
```
`v` contains then the coefficients of your regression so that `y = a+b*x`.
@fipelle’s answer is a shortcut for that process. `v` is then obtained with `v=d\M`.
---
### Author:  [@Johncowk](https://discourse.julialang.org/u/Johncowk)
#### Post date: [November 18, 2019, 8:48pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/8 "2019-11-18T20:48:39Z")
Super clear thanks.
To come back to the original issue, it looks like doing the linear regression this way is slower than with `polyfit` ?
with polyfit :
```julia
0.058074 seconds (48.54 k allocations: 2.364 MiB)
```
with matrix inversion :
```julia
0.061431 seconds (65.20 k allocations: 3.137 MiB)
```
---
### Author:  [@MatFi](https://discourse.julialang.org/u/MatFi)
#### Post date: [November 18, 2019, 9:13pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/9 "2019-11-18T21:13:29Z")
Are you benchmarking in global scope? I wonder because of the allocs. Could type instability play a role here?
---
### Author:  [@Johncowk](https://discourse.julialang.org/u/Johncowk)
#### Post date: [November 18, 2019, 9:16pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/10 "2019-11-18T21:16:06Z")
What do you mean by global scope ?
---
### Author:  [@MatFi](https://discourse.julialang.org/u/MatFi)
#### Post date: [November 18, 2019, 9:25pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/11 "2019-11-18T21:25:41Z")
I think [Performance tips](https://docs.julialang.org/en/v1/manual/performance-tips/index.html#Measure-performance-with-%5B@time%5D(@ref)-and-pay-attention-to-memory-allocation-1) gives the answer to that quastion.
---
### Author:  [@Paul\_Soderlind](https://discourse.julialang.org/u/Paul_Soderlind)
#### Post date: [November 18, 2019, 9:45pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/12 "2019-11-18T21:45:06Z")
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.
---
### Author:  [@mkborregaard](https://discourse.julialang.org/u/mkborregaard)
#### Post date: [November 18, 2019, 10:09pm UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/13 "2019-11-18T22:09:36Z")
I find it rather dissatisfying that `linreg` was replaced, with the message of “just do `(x'*x)\x'*y`”, like “it’s only 5 more characters and much more explicit”. Right but a small percentage of people might not have the formula for linear regression memorized at all times.
---
### Author:  [@Johncowk](https://discourse.julialang.org/u/Johncowk)
#### Post date: [November 19, 2019, 12:22am UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/14 "2019-11-19T00:22:01Z")
To answer your question : no I was not running my benchmark on global scope 😕
---
### Author:  [@oxinabox](https://discourse.julialang.org/u/oxinabox)
#### Post date: [November 19, 2019, 12:41am UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/15 "2019-11-19T00:41:04Z")
> [@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.
Something about making it positive semidefiniate so that `\` will use cholesky ?
(Don’t quote me on that)
---
### Author:  [@fipelle](https://discourse.julialang.org/u/fipelle)
#### Post date: [November 19, 2019, 12:46am UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/16 "2019-11-19T00:46:26Z")
> [@mkborregaard](#):
>
> Right but a small percentage of people might not have the formula for linear regression memorized at all times
This is perfectly reasonable. In my previous post I was trying to expand on the use of `\`. If you do not use regressions or modelling very often it is more convenient to use something like `linreg`.
However, I think that a technical audience would find somewhat harder to memorise the name of a function for simple regressions - plus with small changes in the formula you can also do ridge 🙂
---
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [November 19, 2019, 3:13am UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/17 "2019-11-19T03:13:14Z")
I totally agree, except that the function should be called `linear_regression` or similar, rather than `linreg`, and should probably be in a package like `StatsBase.jl`.
---
### Author:  [@Nosferican](https://discourse.julialang.org/u/Nosferican)
#### Post date: [November 20, 2019, 6:44am UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/18 "2019-11-20T06:44:21Z")
Using the normal equation will likely use Cholesky decomposition while the `\` operator uses QR decomposition. Cholesky is more efficient which is why it may be faster. QR decomposition is more stable and precise. There are a few variants of Cholesky decomposition implemented in Julia (`LL'`, `LDL'`, Bunch-Kaufman, and their pivoted/sparse variants). Using `factorize` on a dense symmetric/hermitian matrix will use the appropriate one.
```julia
using LinearAlgebra
F = Hermitian(A'A) |> factorize
x = F \ (A'b) # Solves a linear equation system Ax = b
```
Most Cholesky implementations are O(n^3) while most implementations of QR decomposition (e.g., Householder reflections) are O(mn^2) which means if the linear system is overly specified (e.g., many more linearly rows than columns in A) then Cholesky will be faster.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [November 20, 2019, 7:23am UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/19 "2019-11-20T07:23:03Z")
> [@mkborregaard](#):
>
> with the message of "just do `(x'*x)\x'*y` ", like “it’s only 5 more characters and much more explicit”.
Can you please link the discussion with that message? I am rather surprised that this was suggested, one should never use the textbook OLS formula for actual calculation (except for some special cases, eg 1 covariate); it is very easy to run into ill-conditioning in practice.
The polyalgorithm for `\` picks a very common and reasonable solution, based on the QR decomposition. For “small” problems, it is not uncommon to use SVD, for large datasets an iterative algorithm is often used.
Since the original problem of this discussion seems to be 1D, just using the specific formula without forming the matrices may be fastest:
> **[Ordinary least squares | Simple linear regression model](https://en.wikipedia.org/wiki/Ordinary_least_squares#Simple_linear_regression_model)**
>
> If the data matrix X contains only two variables, a constant and a scalar regressor xi, then this is called the "simple regression model". This case is often considered in the beginner statistics classes, as it provides much simpler formulas even suitable for manual calculation. The parameters are commonly denoted as (α, β):
---
### Author:  [@mkborregaard](https://discourse.julialang.org/u/mkborregaard)
#### Post date: [November 20, 2019, 9:44am UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/20 "2019-11-20T09:44:12Z")
I was paraphrasing, not arguing that someone exactly said that. But that is the gist I get - every time someone asks “why was linreg removed” the answer is “just use `\`”, and then the asker goes “huh” because of course that ignores the intercept.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [November 20, 2019, 9:58am UTC](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/21 "2019-11-20T09:58:09Z")
I don’t think I have seen people asking for `linreg` for more than a year now, but I have a very different impression: when someone asks about linreg or linear regression in general, generally they are just directed to the appropriate package, eg:
> [@0.7: how to call linreg?](https://discourse.julialang.org/t/0-7-how-to-call-linreg/13105):
>
> I did add StatsBase as instructed, but still fail: $ julia \_ \_ \_ \_(\_)\_ | A fresh approach to technical computing (\_) | (\_) (\_) | Documentation: https://docs.julialang.org \_ \_ \_| |\_ \_\_ \_ | Type "?" for help, "]?" for Pkg help. | | | | | | |/ \_` | | | | |\_| | | | (\_| | | Version 0.7.0 (2018-08-08 06:46 UTC) \_/ |\_\_'\_|\_|\_|\_\_'\_| | Official http://julialang.org/ release |\_\_/ | x86\_64-pc-linux-gnu julia\> x=[1:50;] 50-eleme…
> [@Simple linear regresion in Julia? (and a comment on the stability of user interfaces)](https://discourse.julialang.org/t/simple-linear-regresion-in-julia-and-a-comment-on-the-stability-of-user-interfaces/20867):
>
> I have two arrays of the same dimension, and I want to do linear regression (least squares) between them. After searching in google, It seems that Julia used to have a very useful simple function linreg, for doing that But it seems that it has been removed in 1.1 (or moved to some package?) How can I replace it? Please help me! Just one comment: a thing that I really don’t like about Julia is that everything got deprecated every time. This makes impossible for new newcomers to pick a exa…
There are now excellent packages for various implementations of linear regression and more general models. They should be the first choice, there is little need to implement OLS in Julia for most users now.
[Next page](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232.md?page=2)