Unfortunately, this creates a dense matrix. (F.Q
is not stored explicitly, but as a composition of sparse or low-rank operations so that you can compute F.Q
times a vector efficiently. Trying to pull out a subset of the columns of Q
discards this structure.)
Instead, I would compute (F.Q' * b[F.prow])[1:rnk]
, which should be equivalent and efficient without constructing a dense matrix.
Note that this is also not a sparse matrix in general, since inv(R11)
is in general not sparse. In general you want to compute the action of matrices like this on vectors rather than computing the matrix itself. But it may be okay for your purposes since maybe you are interested only in the case where the rank is very low, so a dense R11
matrix inverse is acceptable?
I’ve just skimmed your code looking for common sparse-matrix gotchas; I haven’t really looked at the structure of your algorithm. Note that there are some papers on minimum-norm solutions with sparse matrices that may be helpful; for example, Solution of Sparse Underdetermined Systems of Linear Equations (1984)
See also Finding the least norm solution of least squares problem with a sparse matrix for iterative methods.