MixedModels.jl icon indicating copy to clipboard operation
MixedModels.jl copied to clipboard

Extend rankUpdate! method to allow for upper tri

Open dmbates opened this issue 3 years ago • 2 comments

  • Although this commit allows for rankUpdate! of the upper triangle it turns out that this is much slower than updating the lower triangle.
julia> using BenchmarkTools, LinearAlgebra, MixedModels
[ Info: Precompiling MixedModels [ff71e718-51f3-5ec2-a782-8ffcbfa3c316]

julia> const contrasts = Dict{Symbol,Any}(:subj => Grouping(), :s => Grouping(), :d => Grouping(), :dept => Grouping(), :service => EffectsCoding());

julia> const form = @formula(y ~ 1 + service + (1|s) + (1|d) + (1|dept));

julia> insteval = MixedModels.dataset(:insteval);

julia> iem01 = fit(MixedModel, form, insteval; contrasts)
Minimizing 115 	 Time: 0:00:01 (16.15 ms/it)
Linear mixed model fit by maximum likelihood
 y ~ 1 + service + (1 | s) + (1 | d) + (1 | dept)
    logLik     -2 logLik       AIC         AICc          BIC     
 -118860.8844  237721.7688  237733.7688  237733.7699  237788.9926

Variance components:
            Column    Variance  Std.Dev. 
s        (Intercept)  0.1059724 0.3255340
d        (Intercept)  0.2652055 0.5149811
dept     (Intercept)  0.0061676 0.0785341
Residual              1.3864886 1.1774925
 Number of obs: 73421; levels of grouping factors: 2972, 1128, 14

  Fixed-effects parameters:
─────────────────────────────────────────────────────
                  Coef.  Std. Error       z  Pr(>|z|)
─────────────────────────────────────────────────────
(Intercept)   3.23629    0.0281513   114.96    <1e-99
service: Y   -0.0462943  0.00669159   -6.92    <1e-11
─────────────────────────────────────────────────────

julia> @time fit(MixedModel, form, insteval; contrasts);
Minimizing 115 	 Time: 0:00:00 ( 8.68 ms/it)
  1.082987 seconds (310.67 k allocations: 52.726 MiB, 2.81% gc time)

julia> const C = zeros(1128, 1128);

julia> const A = iem01.L[2];

julia> @benchmark MixedModels.rankUpdate!(CC, AA, 1.0, 1.0) setup=(CC=Symmetric(fill!(C, 0.0), :L); AA=A)
BenchmarkTools.Trial: 2676 samples with 1 evaluation.
 Range (min … max):  1.583 ms …  2.572 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.621 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.629 ms ± 42.664 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

        ▃▃▄▄▇▇█▆▄▄▂                                           
  ▁▁▂▃▄████████████▇▇▅▅▅▅▄▄▃▃▃▃▃▂▂▂▂▂▂▂▂▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
  1.58 ms        Histogram: frequency by time        1.75 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark MixedModels.rankUpdate!(CC, AA, 1.0, 1.0) setup=(CC=Symmetric(fill!(C, 0.0), :U); AA=sparse(transpose(A)))
BenchmarkTools.Trial: 26 samples with 1 evaluation.
 Range (min … max):  193.801 ms … 194.178 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     194.037 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   194.027 ms ± 100.684 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▁   ▁▁                   ▁ ███ ▁  ▁▁    ▁▁█ ▁ ▁  ▁█    ▁  █ ▁  
  █▁▁▁██▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁███▁█▁▁██▁▁▁▁███▁█▁█▁▁██▁▁▁▁█▁▁█▁█ ▁
  194 ms           Histogram: frequency by time          194 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.
  • This particular modification is not helpful but it is probably worthwhile using this as a starting point for me to review and modernize much of the linear algebra code in the package.
  • In particular, I would plan to use Base.require_one_base_indexing and axes liberally.

dmbates avatar Aug 15 '22 17:08 dmbates

Codecov Report

Base: 94.17% // Head: 94.24% // Increases project coverage by +0.06% :tada:

Coverage data is based on head (cf45767) compared to base (04c2724). Patch coverage: 100.00% of modified lines in pull request are covered.

:exclamation: Current head cf45767 differs from pull request most recent head 8358103. Consider uploading reports for the commit 8358103 to get more accurate results

Additional details and impacted files
@@            Coverage Diff             @@
##             main     #636      +/-   ##
==========================================
+ Coverage   94.17%   94.24%   +0.06%     
==========================================
  Files          29       29              
  Lines        2766     2797      +31     
==========================================
+ Hits         2605     2636      +31     
  Misses        161      161              
Impacted Files Coverage Δ
src/linalg/cholUnblocked.jl 100.00% <100.00%> (ø)
src/linalg/logdet.jl 100.00% <100.00%> (ø)
src/linalg/rankUpdate.jl 98.19% <100.00%> (+0.69%) :arrow_up:

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

:umbrella: View full report at Codecov.
:loudspeaker: Do you have feedback about the report comment? Let us know in this issue.

codecov[bot] avatar Aug 15 '22 17:08 codecov[bot]

xref: #579

palday avatar Aug 16 '22 18:08 palday

@palday Should this PR be merged with #637? If so, could you do rebase one or the other - I still bear the scars of trying to rebase the julia repository when I was just learning git.

dmbates avatar Aug 27 '22 14:08 dmbates

#637 pulled in

palday avatar Aug 29 '22 02:08 palday