# GLM is slow on large datasets. Using OnlineStats for regressions? MixedModels?

**URL:** <https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695>\
**Category:** Performance\
**Tags:** glm\
**Created:** [November 19, 2018, 12:49am UTC](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695 "2018-11-19T00:49:40Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Juan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juan/32/7657_2.png) [@Juan](https://discourse.julialang.org/u/Juan)\
**Post date:** [November 19, 2018, 12:49am UTC](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695/1 "2018-11-19T00:49:40Z")

</div>

Hello.

I’m planning to move from R to Julia and doing some tests about how to properly deal with large datasets and do simple tasks like regressions or survival analysis.

I’ve done a benchmark with R (microbenchmark) for the regressions.

> N ← 3000  
> x1 ← rep(1:N, N)  
> x2 ← rep(1:N, each = N)  
> x3 ← sqrt(rep(1:N^2))  
> x1x2 ← x1_x2  
> gg ← rep(1:5, each=N^2/5)  
> y ← 1-2_x1+3_x2+0.5_x1x2+rnorm(N^2)+x3\*rnorm(N^2)  
> dat ← data.frame(y,x1,x2,x1x2,x3, gg)  
> dat2 ← cbind(1,x1,x2,x1x2,x3,gg)

> lm(y ~ x1 + x2 + x1:x2 + x3, data = dat)  
> coef(.lm.fit(dat2, y))  
> lmfit(dat2, y)$be

Results:

> lm 3.72s  
> .lm.fit 1.66s  
> Rfast lmfit 0.67s

Now I’ve tried to reproduce something similar on Julia with GLM and FixedEffectModels. This is my first day and maybe I need to improve many things.

> N=3000  
> x1 = repeat(1:N, outer=N)  
> x2 = repeat(1:N, inner=N)  
> x3 = sqrt.(repeat(1:N^2))  
> x1x2 = x1 .\* x2  
> gg = repeat(1:5,inner=trunc(Int,N^2/5))  
> y = 1 .- 2_x1 + 3_x2 + 0.5\*x1x2 + rand(N^2) + x3.\*rand(N^2)  
> data = DataFrame(x1=x1, x2=x2, x3=x3, x1x2=x1x2, y=y,gg=gg)  
> categorical!(data, :gg)

These are the results after several runs: (@time)

> fit(LinearModel, @formula(y ~ x1+x2+x3+x1x2), data)  
> reg(data, @model(y ~ x1 + x2 + x3 + x1x2 ))  
> GLM.jl 2.4s  
> FixedEffectModels.jl 2.1s

They are both slightly faster than lm() but much slower than Rfast.  
I’m using 4 threads on Julia and I have left R by default, I guess it uses 4 threads as well for these packages.  
How can I make Julia regressions faster?  
I will be using them with more complex models.

I’m not confident using the package FixedEffectModels because it has an option to specify “fixed effects” but I don’t understand the difference between that effects and the any other regular variable you use, both continuous or categorical. It’s not supposed to have random effects, then I don’t understand the distinction.

I’ve also tried to fit a random effects model with MixedModels.jl

> fit(LinearMixedModel, @formula(y ~x1+x2+x3+x1x2 + (1|gg)), data)

It just takes 2.3s, almost the same than the GLM model. Maybe because it’s very simple or because the package is very optimized.  
lme4 can’t do it with this large dataset. With smaller ones lme4 it takes 6 times more time, and mgcv bam lies in between.

If I wanted to work with larger datasets or more complex models none of these solutions work. R has a Revoscaler package but it doesn’t do regressions with random effects.  
Julia also complains about the memory and MixedEffects says “ERROR: LinearAlgebra.RankDeficientException(1)”.

I think with Julia I have to try with OnlineStats.jl.  
How can I run the former regressions (simple and random effects) using OnlineStats or JuliaDB or other solution?  
For example for N=30000

---

<div class="post-metadata">

**Author:** ![Yifan\_Liu](https://avatars.discourse-cdn.com/v4/letter/y/4da419/32.png) [@Yifan\_Liu](https://discourse.julialang.org/u/Yifan_Liu)\
**Post date:** [November 19, 2018, 1:44am UTC](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695/2 "2018-11-19T01:44:32Z")

</div>

Have you tried XGBoost?

---

<div class="post-metadata">

**Author:** ![Juan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juan/32/7657_2.png) [@Juan](https://discourse.julialang.org/u/Juan)\
**Post date:** [November 19, 2018, 1:54am UTC](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695/3 "2018-11-19T01:54:32Z")

</div>

I didn’t know it, I’m going to try.  
Do you know how to easily translate my example to XGBoost?

I’m also trying Alistair.jl but it doesn’t compile with Julia 1.0.

I’ve also seen LsqFit.jl but I can’t figure out how to use it with my data.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [November 19, 2018, 6:06am UTC](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695/4 "2018-11-19T06:06:31Z")

</div>

Do XGBoost and GLM have anything in common except for them doing regression? GLM fits linear models whereas XGBoost fits boosted regression trees?

---

<div class="post-metadata">

**Author:** ![matthieu](https://avatars.discourse-cdn.com/v4/letter/m/da6949/32.png) [@matthieu](https://discourse.julialang.org/u/matthieu)\
**Post date:** [November 19, 2018, 2:23pm UTC](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695/5 "2018-11-19T14:23:40Z")

</div>

The Rfast package is faster because it accepts a matrix instead of a dataframe, and does not compute any standard error. Just use `\` in Julia to do this:

```julia
N=3000
x1 = repeat(1:N, outer=N)
x2 = repeat(1:N, inner=N)
x3 = sqrt.(repeat(1:N^2))
x1x2 = x1 .* x2
gg = repeat(1:5,inner=trunc(Int,N^2/5))
y = 1 .- 2x1 + 3x2 + 0.5*x1x2 + rand(N^2) + x3.*rand(N^2)
dat2 = hcat(ones(length(x1)), x1, x2, x1x2, x3, gg)
@time dat2 \ y
#> 0.773727 seconds (84 allocations: 480.724 MiB, 2.50% gc time)
```

---

<div class="post-metadata">

**Author:** ![Juan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juan/32/7657_2.png) [@Juan](https://discourse.julialang.org/u/Juan)\
**Post date:** [November 19, 2018, 4:49pm UTC](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695/6 "2018-11-19T16:49:46Z")

</div>

Nice.  
I guess you mean:

> hcat(ones(length(x1)), x1, x2, x1x2, x3)

without the gg, because I didn’t include it on the simple models without random effects.  
On my computer the speed of your **dat2 \ y** now is 1.4s, faster than GLM but still much slower than Rfast. And I expected much more speed from Julia.  
The good news is it also uses much less memory (1/4).  
And to be fair if we wanted to use this process many times on a loop we should also compare the time needed to create the hcat(ones(…)) matrix with the time to create the dataframe.

How can I further accelerate the speed of your division method?  
Maybe using a different kind of structure or container for the data?  
Using a different algorithm or using special macros?

Maybe the speed on R has also to do with it using MKL.

---

<div class="post-metadata">

**Author:** ![Juan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juan/32/7657_2.png) [@Juan](https://discourse.julialang.org/u/Juan)\
**Post date:** [November 19, 2018, 5:35pm UTC](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695/7 "2018-11-19T17:35:33Z")

</div>

I’ve found a much faster solution for the regression with fixed effects,

> using MultivariateStats  
> llsq(dat2, y; bias=false)

It’s just 0.21 seconds, though it doesn’t give you any information about variances, p-values…

The documentation it’s not clear but I’ve found it accepts the same matrix structure, the first column are ones, the other are the input variables.

I hope Alistairs.jl gets updated and we can try it on Julia 1.0

---

<div class="post-metadata">

**Author:** ![matthieu](https://avatars.discourse-cdn.com/v4/letter/m/da6949/32.png) [@matthieu](https://discourse.julialang.org/u/matthieu)\
**Post date:** [November 19, 2018, 6:37pm UTC](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695/8 "2018-11-19T18:37:22Z")

</div>

All these benchmarks are about matrix division / multiplication etc. Julia is not better than matlab or python for this, since all these languages call BLAS.

---

<div class="post-metadata">

**Author:** ![matthieu](https://avatars.discourse-cdn.com/v4/letter/m/da6949/32.png) [@matthieu](https://discourse.julialang.org/u/matthieu)\
**Post date:** [November 19, 2018, 6:46pm UTC](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695/9 "2018-11-19T18:46:52Z")

</div>

Btw, if you look at the source code of `llsq`, in the background, `llsq` basically does

```julia
cholesky!(Symmetric(dat2' * dat2)) \ (dat2' *y)

```

that is, it solves the least square problem using a cholesky decomposition of the matrix dat2’dat2.  
When you do

```julia
dat2 \ y

```

it solves the least squares problem using a QR decomposition of the matrix dat2. The QR method is usually slower but more precise (especially when the system is highly collinear). That is why it is the default — in most cases though, the cholesky method works well.

---

<div class="post-metadata">

**Author:** ![Juan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juan/32/7657_2.png) [@Juan](https://discourse.julialang.org/u/Juan)\
**Post date:** [November 24, 2018, 1:59am UTC](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695/10 "2018-11-24T01:59:50Z")

</div>

I’ve seen the code from Alistair.jl, it doesn’t work on Julia but it’s supposed to be fast.

> <https://github.com/giob1994/Alistair.jl/blob/master/src/linregress.jl>

and they are using:  
`(X' * X) \ (X' * Y)`

What’s the advantage of this method instead of X \ Y?

---

<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 24, 2018, 6:55am UTC](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695/11 "2018-11-24T06:55:11Z")

</div>

> [@Juan](#):
>
> it doesn’t work on Julia but it’s supposed to be fast

Writing code that does not work but is supposed to be fast must be a niche industry.

> [@Juan](#):
>
> What’s the advantage of this method instead of X \ Y?

None. But as a compensation, it comes with plenty of disadvantages (numerical stability).

---

<div class="post-metadata">

**Author:** ![Juan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juan/32/7657_2.png) [@Juan](https://discourse.julialang.org/u/Juan)\
**Post date:** [November 25, 2018, 3:01am UTC](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695/12 "2018-11-25T03:01:19Z")

</div>

I’ve being trying and it seems that

`(dat2' * dat2) \ (dat2' * y)`  
gives the same coefficients than  
`lm(@formula(y ~ x1+x2+x3+x1x2), data)`

Then the method should be right or at least GLM uses the same method.

On the other hand  
`dat2 \ y`  
gives coefficients slightly different and it’s much slower.

(dat2’ \* dat2) \ (dat2’ \* y) is faster because the products produce small matrices and the inverse is easier to calculate than calculating the inverse of a large matrix.

Then what should I trust?

---

<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:** [November 25, 2018, 3:57am UTC](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695/13 "2018-11-25T03:57:45Z")

</div>

`(X'X)\(X'y)` is the canonical way that [OLS](https://en.wikipedia.org/wiki/Ordinary_least_squares) is estimated.

But what are you trying to do? In my experience, linear regression without standard errors is not very useful.

---

<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 25, 2018, 8:39am UTC](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695/14 "2018-11-25T08:39:23Z")

</div>

> [@Juan](#):
>
> gives the same coefficients than

Mathemathically, of course it does. To see the error, you would need an ill-conditioned matrix.

> [@Juan](#):
>
> Then what should I trust?

Any reasonable text on numerical linear algebra?

---

<div class="post-metadata">

**Author:** ![andreasnoack](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasnoack/32/27_2.png) [@andreasnoack](https://discourse.julialang.org/u/andreasnoack)\
**Post date:** [November 25, 2018, 2:04pm UTC](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695/15 "2018-11-25T14:04:23Z")

</div>

You can pass the model matrix directly to `fit` instead of using a `Formula` and `DataFrame`. That will save you the costly computation of the model matrix. I.e.

```julia
julia> @time ft = fit(LinearModel, @formula(y ~ x1+x2+x3+x1x2), data);
  1.254742 seconds (540 allocations: 1.745 GiB, 37.68% gc time)

julia> @time fit(LinearModel, ft.mm.m, data.y);
  0.551531 seconds (26 allocations: 823.976 MiB, 41.03% gc time)

```

We are already using the Cholesky.

---

<div class="post-metadata">

**Author:** ![Juan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juan/32/7657_2.png) [@Juan](https://discourse.julialang.org/u/Juan)\
**Post date:** [November 25, 2018, 5:01pm UTC](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695/16 "2018-11-25T17:01:46Z")

</div>

How do you create that “ft.mm.m” without first running the first line?  
I mean, if I want to use your second syntax but also need the first one I’ll end spending more time. If you already have the first why do you need the second?

Do you mean GLM.jl always uses Cholesky or the latter syntax with the model matrix uses Cholesky?  
@matthieu mentioned QR is more stable and precise.

---

<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:** [November 25, 2018, 7:31pm UTC](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695/17 "2018-11-25T19:31:02Z")

</div>

`Matrix(df)` will work, or you can dig into the unexported functions in `StatsModels` that are used to construct a model matrix.

---

<div class="post-metadata">

**Author:** ![Juan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juan/32/7657_2.png) [@Juan](https://discourse.julialang.org/u/Juan)\
**Post date:** [November 25, 2018, 8:37pm UTC](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695/18 "2018-11-25T20:37:23Z")

</div>

Does it need to have the same structure than the dataframe or the same than dat2 = hcat(ones(length(x1)), x1, x2, x1x2, x3) ?

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [November 25, 2018, 9:39pm UTC](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695/19 "2018-11-25T21:39:56Z")

</div>

I’m not a numeric analyst, but won’t the condition number of `(X' * X) \ (X' * Y)` be the square of that of `X \ Y`? Theoretically, they should give the same result, but for ill-conditioned problems, the first formulation is horrible in comparison? With the original problem being to find the LS solution of `X*b = Y`, a standard procedure would be to do the QR factorization of X. I assume that `X \ Y` is doing that.

If you want to test this for yourself: splitting up the QR factorization, `X = Q*R = [Q1 Q2]*[R1;R2]` with `Q1` having `r = rank(X)` columns and using `Q'*Q = I`, leads to two sets of equations, [ahem… fixing this: `R1*X` and `R2*X` should of course be `R1*b` and `R2*b`…] `R1*b = Q1' * Y` and `R2*b = Q2' * Y`. The first is trivial to solve as `R1` is upper triangular. For the second, in principle `R2 = 0`, so one has `0 = Q2'*Y`: if the model `X*b` is a perfect description of the data `Y`, then `Y` should lie in the nullpace of `Q2'`, hence `Q2'*Y = 0`. With an imperfect model (e.g., noisy data), `Q2'*Y` contains information about the squared error; this equation is ditched.

Of course, in doing the QR factorization, one doesn’t need `Q` and `R` – it suffices to keep `Q1` and `R1` – which normally are much smaller than `Q` and `R`.

For many problems, `X` is of full rank, and then the algorithm is straightforward. If `X` has some rank loss, it is necessary to find the rank. This is sometimes done using costly `SVD` computations. There is a “rank revealing” QR factorization which orders the main diagonal elements of `R`, but I don’t think that algorithms is available in Julia (???). There is a MATLAB version of it, but apparently the MATLAB implementation is copyrighted.

---

<div class="post-metadata">

**Author:** ![Nosferican](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nosferican/32/9275_2.png) [@Nosferican](https://discourse.julialang.org/u/Nosferican)\
**Post date:** [November 26, 2018, 12:20am UTC](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695/20 "2018-11-26T00:20:22Z")

</div>

```julia
using LinearAlgebra
X = hcat(ones(100), rand(100, 2))
X = hcat(X, 1.5 * X[:,2])
F = qr(X, Val(true))
good = abs.(diag(F.R)) .> √eps()
@views if ~all(good)
    X = X[:,good[invperm(F.jpvt)]]
    F = qr(X, Val(true))
end

```

You can use `F.jpvt` to get the full rank matrix.

[Next page](https://discourse.julialang.org/t/glm-is-slow-on-large-datasets-using-onlinestats-for-regressions-mixedmodels/17695.md?page=2)