I would like to give some closure to the topic.
I rewrote the code to use regular matrices instead of sparse and I could leverage mul!(), ldiv!(), and @view very efficiently. I got tremendous speedup and much less allocations.
Here is the function itself with 4 implementations, where the first two use sparse matrices and the last two are without, only measured with @time
0.721347 seconds (944 allocations: 266.088 MiB, 41.04% gc time)
0.727664 seconds (914 allocations: 265.672 MiB, 40.43% gc time)
0.158879 seconds (744 allocations: 33.713 MiB, 1.27% gc time)
0.170166 seconds (722 allocations: 33.677 MiB, 2.59% gc time)
And here is the function as part of a larger model:
16.876194 seconds (391.23 k allocations: 8.699 GiB, 23.92% gc time, 1.13% compilation time)
14.811057 seconds (263.82 k allocations: 8.685 GiB, 16.73% gc time, 0.33% compilation time)
6.620201 seconds (260.26 k allocations: 4.154 GiB, 16.99% gc time, 0.61% compilation time)
6.404180 seconds (260.03 k allocations: 4.154 GiB, 12.55% gc time, 0.63% compilation time)
At least for my application, where I am doing a lot of basic linear algebra operations on the same matrices, even with very high sparsity (context dependent), sparse matrices appear not to be the optimal choice. The calculations that I am doing, for context of readers of this involve calculating
\mathbf{K}_\beta=\mathbf{V}_{\mathrm{Minn}}^{-1}+\mathbf{X}^{\prime}\left(\mathbf{I}_T \otimes \widehat{\Sigma}^{-1}\right) \mathbf{X}, and
\mathbf{C}_{\mathbf{K}_\beta}^{\prime} \backslash\left(\mathbf{C}_{\mathbf{K}_\beta} \backslash\left(\mathbf{V}_{\operatorname{Minn}}^{-1} \beta_{\operatorname{Minn}}+\mathbf{X}^{\prime}\left(\mathbf{I}_T \otimes \widehat{\Sigma}^{-1}\right) \mathbf{y}\right)\right),
where \mathbf{C}_{\mathbf{K}_\beta} is the lower Cholesky factor and \mathbf{C}_{\mathbf{K}_\beta}' is the upper Cholesky factor and the matrices are of the size, for example, (3555, 1140) for \mathbf{X} and \mathbf{K}_\beta is 1140 \times 1140. For some applications, as shown above, the density is around 5%, while for others some matrices are not sparse at all.