Speed up estimate_s_rss

[aksarkar]

Pre-compute some quantities outside optim to speed up the diagnostic.

Profiling code

devtools::load_all()
data(N3finemapping)
ss <- compute_suff_stat(N3finemapping$X, N3finemapping$Y[,1])
Rprof()
s <- estimate_s_rss(ss$Xty / ss$n, ss$XtX / ss$n, ss$n)
Rprof(NULL)
print(summaryRprof()$by.total)

On revision f63f9cb

                  total.time total.pct self.time self.pct
"estimate_s_rss"        2.84    100.00      0.00     0.00
"optim"                 1.94     68.31      0.00     0.00
"f"                     1.92     67.61      0.00     0.00
"fn"                    1.92     67.61      0.00     0.00
"optimize"              1.92     67.61      0.00     0.00
".External2"            1.90     66.90      0.00     0.00
"<Anonymous>"           1.90     66.90      0.00     0.00
"t.default"             1.10     38.73      1.10    38.73
"standardGeneric"       1.10     38.73      0.00     0.00
"t"                     1.10     38.73      0.00     0.00
"eigen"                 0.64     22.54      0.60    21.13
"%*%"                   0.56     19.72      0.56    19.72
"*"                     0.26      9.15      0.26     9.15
"/"                     0.26      9.15      0.26     9.15
"is.finite"             0.02      0.70      0.02     0.70
"lazyLoadDBfetch"       0.02      0.70      0.02     0.70
"structure"             0.02      0.70      0.02     0.70

New code

                 total.time total.pct self.time self.pct
"estimate_s_rss"       1.00       100      0.00        0
"eigen"                0.70        70      0.64       64
"/"                    0.28        28      0.28       28
"structure"            0.04         4      0.04        4
"is.finite"            0.02         2      0.02        2
"sum"                  0.02         2      0.02        2
".External2"           0.02         2      0.00        0
"<Anonymous>"          0.02         2      0.00        0
"f"                    0.02         2      0.00        0
"fn"                   0.02         2      0.00        0
"optim"                0.02         2      0.00        0
"optimize"             0.02         2      0.00        0

[gaow]

Thanks @aksarkar it makes complete sense. Do all unit test pass?

[aksarkar]

Yes, all tests passed on my machine.

> devtools::test()

ℹ Testing susieR
✔ | F W  S  OK | Context
✔ |         18 | test_compute_tf.R
✔ |          2 | test_compute_colsds.R
✔ |          2 | test_Eloglik_rss.R
✔ |          2 | test_Eloglik.R
✔ |          2 | test_est_V.R
✔ |          2 | test_get_ER2_rss.R
✔ |          2 | test_get_ER2.R
✔ |          4 | test_get_pip.R
✔ |          8 | test_get_samples.R [3.7s]
✔ |          2 | test_init.R
✔ |          6 | test_intercept_standardize.R [4.7s]
✔ |         34 | test_null_weight.R [4.5s]
✔ |         45 | test_prior_weights.R [5.8s]
✔ |          2 | test_SER_posterior_e_loglik_rss.R
✔ |          2 | test_SER_posterior_e_loglik.R
✔ |          6 | test_sparse_multiplication.R
✔ |          4 | test_sparse_set_X_attributes.R
✔ |         30 | test_susie_beta_se.R [2.1s]
✔ |          2 | test_susie_get_cs.R
✔ |          2 | test_get_objective.R
✔ |          2 | test_get_objective_rss.R
✔ |      2   0 | test_susie_rss.R
✔ |        151 | test_susie_trendfilter.R [1.8s]
✔ |         30 | test_susie_XtX_Xty.R
⠏ |          0 | test_update_each_effect.R
══ Results ═══════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════
Duration: 25.9 s

── Skipped tests (2) ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
• susie_rss_lambda is not currently being used (2): test_susie_rss.R:5:3, test_susie_rss.R:21:3

[ FAIL 0 | WARN 0 | SKIP 2 | PASS 360 ]

[pcarbo]

Looks good, @aksarkar. I can confirm that the tests pass, and the result was the same in the few tests that I ran.

@gaow Shall I go ahead and merge this PR?

[gaow]

@aksarkar @pcarbo that's great! Okay I merged it. Thanks a lot @aksarkar !