The short version is: I have a sparse matrix A like
julia> m1.R[1,2]
6040x3706 sparse matrix with 1000209 Float64 entries:
[1 , 1] = 0.0874371
[6 , 1] = 0.0765784
[8 , 1] = 0.0558497
[9 , 1] = 0.0635299
[10 , 1] = 0.0333554
[18 , 1] = 0.0381539
[19 , 1] = 0.041645
[21 , 1] = 0.126607
[23 , 1] = 0.0382153
[26 , 1] = 0.0333964
[28 , 1] = 0.0632488
[34 , 1] = 0.0515869
[36 , 1] = 0.035613
[38 , 1] = 0.0652987
[44 , 1] = 0.0476856
[45 , 1] = 0.0386538
[48 , 1] = 0.0273836
[49 , 1] = 0.0629714
[51 , 1] = 0.0989613
[56 , 1] = 0.0786409
[60 , 1] = 0.0770788
[65 , 1] = 0.0596687
[68 , 1] = 0.0760876
[73 , 1] = 0.041645
[75 , 1] = 0.0499972
⋮
[4751, 3706] = 0.113348
[4790, 3706] = 0.0343391
[4802, 3706] = 0.0264789
[4816, 3706] = 0.0346104
[4823, 3706] = 0.0308447
[4831, 3706] = 0.0431582
[4834, 3706] = 0.039785
[4858, 3706] = 0.0616373
[4939, 3706] = 0.0501356
[5049, 3706] = 0.0478057
[5074, 3706] = 0.0243015
[5087, 3706] = 0.0387177
[5100, 3706] = 0.0229637
[5205, 3706] = 0.0353664
[5304, 3706] = 0.104785
[5333, 3706] = 0.0231531
[5359, 3706] = 0.0296257
[5405, 3706] = 0.0608765
[5475, 3706] = 0.0343839
[5602, 3706] = 0.0742148
[5682, 3706] = 0.0271598
[5812, 3706] = 0.0224472
[5831, 3706] = 0.0192226
[5837, 3706] = 0.0289602
[5927, 3706] = 0.045151
[5998, 3706] = 0.0525206
and a dense matrix B, which in this case starts off as diagonal of size
3706, but stored as a full dense matrix. I want to downdate B as
B -= A'A
For each iteration in an optimization of a criterion I modify the values of
the non-zeros in A, go back to the original value of B and do the downdate
again.
At present I am forming the product A'A with a general sparse matrix
product at each iteration. Thus each iteration bears the cost of
allocating and freeing storage and the cost of determining the location of
the nonzeros in B. For the particular model fit
julia> @time fit!(lmm(Rating ~ 1 + (1|Subject) + (1|Item), df))
49.765195 seconds (36.12 M allocations: 6.503 GB, 2.80% gc time)
Linear mixed model fit by maximum likelihood
logLik: -1331986.005811, deviance: 2663972.011622, AIC: 2663980.011622,
BIC: 2664027.274500
Variance components:
Variance Std.Dev.
Subject 0.12985209 0.36034996
Item 0.36879693 0.60728653
Residual 0.81390867 0.90216887
Number of obs: 1000209; levels of grouping factors: 6040, 3706
Fixed-effects parameters:
Estimate Std.Error z value
(Intercept) 3.33902 0.0114624 291.302
where Subject, Item and Rating are from the GroupLens 1M data
set http://files.grouplens.org/datasets/movielens/ml-1m.zip and lmm is in
the MixedModels package, about 2/3 of the time is spent in this downdate
operation.
I have tried a couple of different approaches, such as using the dot
product of sparse vectors as in the v0.5.0-dev file
base/sparse/sparsevector.jl and it takes much longer than does forming A'A,
with its associated overhead. That is good news, in a sense, because it
means the sparse matrix product and sparse matrix transpose are performing
well. But I still have the nagging suspicion that I should be able to do
better because I know the product will be dense and I know I am working
with A'A. Tim Davis has a function to evaluate AA' in CHOLMOD but I don't
want to read his GPL or LGPL code to see how it is done.
