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)
.