I have an algorithm that uses only the first k rows and k columns of a Cholesky factor L (i.e., the kth leading principal submatrix L = L[1:k, 1:k]). Currently I just extract it after the computation of the whole Cholesky factor, but is there any way to have cholesky from LinearAlgebra compute only these elements given some positive definite matrix A?
Moreover, is there a way within Julia to efficiently compute the dot product of the kth row of the Cholesky factor with another vector of length k, or at least compute just the kth row efficiently? Or is the optimal way just taking this principal submatrix above and extracting L[k, 1:k]?
Yup, just compute the Cholesky factor of the submatrix. For example, here is the 5x5 leading principal submatrix of the Cholesky factor of the 100x100 SPD matrix A = rand(100,100); A = A'A, computed the slow way and the fast way:
julia> cholesky(A).L[1:5,1:5] # slow way: form whole Cholesky factor, then take submatrix
5×5 Matrix{Float64}:
5.81234 0.0 0.0 0.0 0.0
4.50388 3.68507 0.0 0.0 0.0
4.1697 1.13275 3.27478 0.0 0.0
3.77458 1.3239 1.19011 3.14419 0.0
4.65957 1.30715 1.30752 0.63263 2.95488
julia> cholesky(@view A[1:5,1:5]).L # fast way: form Cholesky factor from submatrix of A
5×5 LowerTriangular{Float64, Matrix{Float64}}:
5.81234 ⋅ ⋅ ⋅ ⋅
4.50388 3.68507 ⋅ ⋅ ⋅
4.1697 1.13275 3.27478 ⋅ ⋅
3.77458 1.3239 1.19011 3.14419 ⋅
4.65957 1.30715 1.30752 0.63263 2.95488
This is a consequence of how Cholesky factorization is defined: a leading principal submatrix of the Cholesky factor of A is always equivalent to a Cholesky factor of the corresponding leading principal submatrix of A.
From above, just do @views dot(cholesky(A[1:k,1:k]).L[k,1:k], vec).