[WIP]Replace Matrix with NMatrix for high performance

[lokeshh]

No description provided.

[lokeshh]

I just replaced all instances of Matrix with NMatrix in logistic regression and instead of getting a speed up, the benchmark suggests there's 130% increase in time when using NMatrix compared to Matrix.

[agisga]

Not sure if this check is necessary. Computing the determinant is very expensive. I think we can continue without it, then methods like invert or solve will complain later on anyway, if h is (nearly) singular.

[lokeshh]

@agarie I tried all the suggestions, there's some improvement but not that much. Here's the benchmark that I used and the result I got.

[agisga]

Interesting... I'll take a look sometime this week, and let you know if I can figure something out...

[agisga]

@lokeshh I just had another look at your benchmark code, and I noticed that you generate data from a perfectly linear equation. When you generate the data according to a statistical linear model, there is always a random noise term added to the linear equation. That is, add another vector of (small) random numbers to line 10 (preferably sampled from the Normal distribution).

I think the problem why NMatrix versions is slower can lie therein in your example. The data is just too perfectly linear (because there isn't any noise), and therefore the algorithm does not have to work hard at all, which means that most time is spend on internal data preparation rather than number crunching.
Once you add a vector of noise term, you may also experiment with the variance of the noise, which you can control by scaling the noise with a scalar.

Anyway, would be great if you can try it out. It's just one little change in line 10 of your code.

[lokeshh]

@agarie I changed the last line to

df[:y] = df.a * 3 + df.b * 5 + df.c * -2 + df.d * 2 + Statsample::Shorthand.rnorm(n)

and Result was almost the same. For Matrix it took 1.2 seconds while for NMatrix(lapack) it took 4 seconds.

[agisga]

@lokeshh Thanks for trying it out! I guess the problem lies somewhere else.
I will try to find some time and write a GLM solver with nmatrix-lapacke from scratch to understand the problem better.

[agisga]

Okay. So, I have written a simple gradient descent algorithm from scratch to fit GLMs with NMatrix.
Here is the (quick and dirty) code: https://github.com/agisga/gradient_descent/blob/master/examples/logistic_regression.rb
The benchmarks don't look good at all 😕 :

       user     system      total        real
statsample-glm  0.090000   0.000000   0.090000 (  0.097168)
GD without backtracking  6.400000   0.010000   6.410000 (  6.397935)
GD with backtracking 18.040000   0.000000  18.040000 ( 18.048868)

which I think is mostly due to the gradient descent method, which requires many more iterations than the algorithm in statsample-glm to converge:

------------------------
(1) Gradient descent without backtracking:
Step size: 0.0036953813884534462
Number of iterations: 258
Estimated optimum:
3.242609900192333
5.192787416360552
-1.953170820366363
1.9075632959144226
------------------------
(2) Gradient descent with backtracking:
Number of iterations: 201
Estimated optimum:
3.2075997913676497
5.136035342275983
-1.9315688832370708
1.8864649165825929
------------------------
(3) Statsample-glm algorithm:
Number of iterations: 7
Estimated optimum:
3.2707387315649643
5.238379138664602
-1.9705222484841016
1.9245196488076983

NB: @lokeshh I was wrong with my comment above about adding random noise to y. Since the response vector of a logistic regression model should be binary, the proper way to generate y is actually

p = df.a * 3 + df.b * 5 + df.c * -2 + df.d * 2
p.map! { |i| 1/(1+Math::exp(-i)) }
y = p.map do |prob|
  unif_draw = rand
  unif_draw < prob ? 1 : 0
end
df[:y] = Daru::Vector.new(y)