Overwrite Cholesky factorization

General Usage
1

I am trying to figure out the best way to overwrite/update a Cholesky factorization. In short, I have a symmetric matrix K which is periodically updated with new data. Every time I update K I want to recompute its Cholesky factorization (call it C). So what I want is a method like:

cholesky!(C::Cholesky, K::Matrix)

But as far as I can tell this doesn’t exist? I do see there is a cholesky!(F::Factor, A::Matrix), but I’m not sure if or how to use that.

I think the snippet below works? But I’m not sure and it feels a bit dangerous… Is there a better alternative?

K = # recompute symmetric matrix K
copyto!(C.U.data, K)
cholesky!(C.U.data)

The function below gives a slightly more expanded sketch of what I’m using this for, in case that’s helpful…

"""
    update_kernel!(K, C, X)

Overwrite `K` with kernel matrix computed between
columns of `X`. Also compute the Cholesky factorization
of `K` and put the result into `C`.
"""
function update_kernel!(K::Matrix, C::Cholesky, X::Matrix)
    
    # Compute kernel matrix, overwrite K.
    d, n = size(X)
    for (i, j) in IterTools.product(1:n, 1:n)
        K[i, j] = exp(-sum((X[:, i] .- X[:, j]).^2))  # RBF kernel, for example...
    end

    # Compute Cholesky factorization of K, overwrite C.
    cholesky!(C, K) # <--- how to do this?
end
2

AFAIK generally, there is no algorithm that is better than simply refactoring, other than rank-1 updates.

3

Sorry, I just meant to ask how to best avoid allocating new memory while doing this. I guess the answer might be that the memory allocation is a small cost compared with computing the factorization itself…

4

I think it would be

cholesky!(Symmetric(copyto!(C.factors, K), :U))

The matrix that actually stores the data for the Cholesky object is the factors field, either in the upper triangle or lower triangle according to the uplo field. If you copy the contents of K to C.factors and then force it to be considered as Symmetric stored in the upper triangle (or lower, if that works better for you) then the cholesky! method will decompose it in place.

2 Likes