# Accelerating linear methods

**URL:** <https://discourse.julialang.org/t/accelerating-linear-methods/99685>\
**Category:** Machine Learning\
**Tags:** question\
**Created:** [May 31, 2023, 7:14pm UTC](https://discourse.julialang.org/t/accelerating-linear-methods/99685 "2023-05-31T19:14:06Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![Leticia-maria](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leticia-maria/32/30981_2.png) [@Leticia-maria](https://discourse.julialang.org/u/Leticia-maria)\
**Post date:** [May 31, 2023, 7:14pm UTC](https://discourse.julialang.org/t/accelerating-linear-methods/99685/1 "2023-05-31T19:14:06Z")

</div>

Dear Julia friends, I am benchmarking Julia for linear regression models and trying to optimize that for high-dimensional spaces. Any suggestions on the code below to reduce allocations and speed up the code? Thank you so much

```julia
using Random
using GLM
using DataFrames
using BenchmarkTools

function benchmark(dimensions::Int)
    Random.seed!(1234)
    X = randn(100000, dimensions)
    β = randn(dimensions)
    y = X * β + 0.1 * randn(100000)

    # Avoiding column name generation
    colnames = [Symbol("x$i") for i in 1:dimensions]
    push!(colnames, :y)
    data = DataFrame(hcat(X, y), colnames) # Creating DataFrame in one step to avoid extra allocations

    formula = @inbounds Term(:y) ~ sum(Term(Symbol("x$i")) for i in 1:dimensions)

    GLM.lm(formula, data)
end

for dims in [10, 100, 1000, 10000]
    println("Warming up for dimension $dims")
    benchmark(dims) # Warm-up run

    println("Benchmarking for dimension $dims")
    @btime benchmark($dims) # Using BenchmarkTools.@btime for more accurate timing and allocation measurement
end

```

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [May 31, 2023, 7:51pm UTC](https://discourse.julialang.org/t/accelerating-linear-methods/99685/2 "2023-05-31T19:51:42Z")

</div>

I just profiled the code with `@profview benchmark(1000)` and on my computer, this call spent roughly

- 50% of the time in `lm`, which is split pretty evenly between `ModelMatrix`, `ModelFrame` and `fit`
- 30% of the time in `DataFrame` (constructing `data`)
- 20% of the time in `randn` (generating `X`)

One would assume the `fit` method does the actual work, everything else is just moving stuff around in memory. What is the context of this function? Will you call it several times? Can you reuse `X` or `data`?

---

<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:** [May 31, 2023, 7:55pm UTC](https://discourse.julialang.org/t/accelerating-linear-methods/99685/3 "2023-05-31T19:55:27Z")

</div>

Have you considered just skipping the dataframe and GLM stuff and just doing `fit = X \ y` directly? That’s all it’s doing in the end.

(Reducing it to `X \ y` lets you further optimize things by customizing the method of the linear solve, doing things in-place, etcetera. It may also trim off a bunch of overhead.)

---

<div class="post-metadata">

**Author:** ![jar1](https://avatars.discourse-cdn.com/v4/letter/j/c0e974/32.png) [@jar1](https://discourse.julialang.org/u/jar1)\
**Post date:** [May 31, 2023, 7:58pm UTC](https://discourse.julialang.org/t/accelerating-linear-methods/99685/4 "2023-05-31T19:58:50Z")

</div>

`lm` also calculates other values like t-statistics etc. It would be cool to have more control over the solver though.

---

<div class="post-metadata">

**Author:** ![pdeffebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pdeffebach/32/10320_2.png) [@pdeffebach](https://discourse.julialang.org/u/pdeffebach)\
**Post date:** [May 31, 2023, 7:59pm UTC](https://discourse.julialang.org/t/accelerating-linear-methods/99685/5 "2023-05-31T19:59:17Z")

</div>

You might want to try [FixedEffectModels.jl](https://github.com/FixedEffects/FixedEffectModels.jl), which also does linear regression but may be faster.

Also, your benchmark includes creating the `DataFrame`. You should split the data creation and estimation process into separate pieces to best analyze the performance of the estimation.

---

<div class="post-metadata">

**Author:** ![Leticia-maria](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leticia-maria/32/30981_2.png) [@Leticia-maria](https://discourse.julialang.org/u/Leticia-maria)\
**Post date:** [May 31, 2023, 11:48pm UTC](https://discourse.julialang.org/t/accelerating-linear-methods/99685/6 "2023-05-31T23:48:28Z")

</div>

Yes, I am actually calling it several times (like thousands), so any small improvement is critical

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [June 1, 2023, 12:10am UTC](https://discourse.julialang.org/t/accelerating-linear-methods/99685/7 "2023-06-01T00:10:35Z")

</div>

`fit` or `X \ y` complexity also scale differently from other operations (cubic in `dimension`, IIRC), and therefore testing `benchmark(10_000)` should spend most time there and optimizing the rest will be less crucial.  
Having said this, for smaller problems it is always nice to reduce overhead.  
One would want to minimize copying of the data (as just generating it takes 20% with dimension 1\_000). Creating the DataFrame with `view`s to columns might help, by adding `copycols=false` like so:

```julia
data = DataFrame(hcat(X, y), colnames; copycols=false)

```

and maybe to avoid `hcat` copying, using some in-place versions:

```julia
function benchmark(dimensions::Int)
    Random.seed!(1234)
    D = Matrix{Float64}(undef, 100_000, dimensions+1)
    randn!(@view D[:,1:end-1])
    β = randn(dimensions)
    @views D[:,end] .= D[:,1:end-1] * β + 0.1 * randn(100000)

    # Avoiding column name generation
    colnames = [Symbol("x$i") for i in 1:dimensions]
    push!(colnames, :y)
    data = DataFrame(D, colnames; copycols=false) # Creating DataFrame in one step to avoid extra allocations
...

```

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 1, 2023, 5:03am UTC](https://discourse.julialang.org/t/accelerating-linear-methods/99685/8 "2023-06-01T05:03:57Z")

</div>

Do you always have the same X? or the same y? Do they at least keep the same dimensions?

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [June 1, 2023, 5:38am UTC](https://discourse.julialang.org/t/accelerating-linear-methods/99685/9 "2023-06-01T05:38:11Z")

</div>

This discussion was had before, much of it is likely unchanged:

> [@Efficient way of doing linear regression](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232):
>
> 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 ?

I agree with Peter, if you are fitting models with high dimensional fixed effects look at FixedEffectModels which is designed for this and also has GPU support.

---

<div class="post-metadata">

**Author:** ![Ajay\_Shah](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ajay_shah/32/10520_2.png) [@Ajay\_Shah](https://discourse.julialang.org/u/Ajay_Shah)\
**Post date:** [June 2, 2023, 6:50am UTC](https://discourse.julialang.org/t/accelerating-linear-methods/99685/10 "2023-06-02T06:50:28Z")

</div>

The R world has a nice separation between a small lm.fit() which is fast and the full lm() which does more things. We should think similarly.

Mousum Dutta has been building a QR-decomposition based engine for GLM.jl which is more precise. Here also, we see the tradeoffs between a fast and precise method. It’s good to have both in the implementation and in the function design.

---

<div class="post-metadata">

**Author:** ![Leticia-maria](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leticia-maria/32/30981_2.png) [@Leticia-maria](https://discourse.julialang.org/u/Leticia-maria)\
**Post date:** [June 3, 2023, 11:15pm UTC](https://discourse.julialang.org/t/accelerating-linear-methods/99685/11 "2023-06-03T23:15:38Z")

</div>

They have the same dimension: I have posted a more detailed code here: [Error when processing data from file for linear regression](https://discourse.julialang.org/t/error-when-processing-data-from-file-for-linear-regression/99824)

---

<div class="post-metadata">

**Author:** ![Leticia-maria](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leticia-maria/32/30981_2.png) [@Leticia-maria](https://discourse.julialang.org/u/Leticia-maria)\
**Post date:** [June 3, 2023, 11:15pm UTC](https://discourse.julialang.org/t/accelerating-linear-methods/99685/12 "2023-06-03T23:15:57Z")

</div>

That would be great, we would have more control