Updating Cholesky Factorization submatrix

Specific Domains Numerics
,
1

I have a positive semidefinite matrix A \in \mathbb{R}^{n \times n} and its Cholesky factorization is A = L L^T. I would like to update the Cholesky factorization by removing the last row and column of L.

What is the correct way to update the Cholesky factorization object returned by cholesky() in this way? I’d like to pass the updated Cholesky factorization object to linear system solvers later in my program. I imagine this is easy to do since we are just throwing data away; but I can’t figure out how to access the Cholesky factorization object in the right way.

Here’s an example of what I’d want.

using LinearAlgebra

# create a positive definite matrix
A = [2 0 1; 0 3 0; 1 0 2] 

# obtain cholesky factorization
C = cholesky(A)
L = C.L 

# I want to update C to have this lower triangular matrix
L_sub = L[1:end-1,1:end-1]

# solve this linear system, but with updated Cholesky
x = (L_sub' * L_sub) \ ones(2) # answer = (1/2, 1/3)
1 Like
2

This is pretty hacky, but I don’t see a better way. Using your C above:

C2 = Cholesky(C.factors[1:end-1, 1:end-1], :U, 0)

This uses the direct Cholesky constructor rather than cholesky, which actually performs the factorization. It needs the factor, whether it is upper (default) or lower triangular, and the return code where 0 is success. I found the constructors here.

1 Like
3

Thanks for your response, @jkbest2! I’ll give this a try in my code. The only thing I’m worried about is the overhead I’d incur by writing C.factors again in memory.

I’ll keep this question open in case anyone has other suggestions. I’ll also take a look through the constructors link you sent.

4

This is getting pretty far outside my experience, but maybe you could avoid allocations by using a view into the original C.factors? It might cause problems down the line because the view wouldn’t be strided so couldn’t be passed to eg BLAS. Not sure about that though.

5

Have you looked at
https://github.com/JuliaLang/julia/issues/2929
and
https://github.com/mpf/QRupdate.jl
?

2 Likes
6

I had known about lowrankupate et al., but missed that they could be used to remove rows/columns. Thanks!

7

I have implemented this functionality in my QP solver QPDAS.jl, see https://github.com/mfalt/QPDAS.jl/blob/master/src/choleskySpecial.jl

It actually sets any row and column to identity, and you are then able to reverse the update.

For your example:

import QPDAS, LinearAlgebra
using LinearAlgebra

A = [2 0 1; 0 3 0; 1 0 2] 
C = cholesky(A)
C2 = QPDAS.CholeskySpecial(C)
# Remove last row/column (actually sets it to identity)
QPDAS.deleterowcol!(C2, 3)
(C2\[ones(2),randn()])[1:2] # answer [1/2,1/3]

Edit: I am not sure how it will handle the semi-definite case, I have only used it for positive-definite factorizations.

1 Like
8

It is a late response :). For a case of positive-definite symmetric matrix you may use this approach:
Algorithms for updating Cholesky factorization
BR,
Igor