Extend rankUpdate! method to allow for upper tri
- 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_indexingandaxesliberally.
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.
xref: #579
@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.
#637 pulled in