Linpack, the predecessor of LAPACK, provided update/downdate functions for Cholesky decompositions (dchud, dchdd) as well as a function, dchex, to restore a Cholesky factor to triangular form after a left or right circular shift of columns.

I plan on using a dchex-like scheme based on the givens function to evaluate variance inflation factors for linear models. I wanted to check before writing such a function if similar capabilities have already been written.

Your timing is off by one day if you wanted the name of the movie. (It says â€ś27 DAYS LATERâ€ť in the browser version.)

I actually donâ€™t want update or downdate, I want the cyclic permutation of the columns. I know how to do it, I just havenâ€™t written it yet. Some people I was working with used a hideously inefficient function to calculate the VIFs so, of course, I had to spend a while thinking about how to do it efficiently.

If you extract a k by k R from the factorization, either Cholesky or QR, then the VIF for the kâ€™th column is (assuming that the first column is for the (Intercept) term - otherwise VIF doesnâ€™t make a whole lot of sense)

sumabs2(view(R, 2:k, k)) / abs2(R[k,k])

(The proof is left as an exercise for the reader.)

To get the other VIFs just rearrange the columns and restore R to upper triangular via Givens rotations. It works best if you use a cyclic shift of a trailing subset of the columns.

Okay, I get it now. It is actually pretty neat that you can do that. However, Iâ€™m wondering how much of a speedup it is relative to something like

X = [ones(n) randn(n,k)]
F = cholfact(X'X)
diag(inv(F)) .* sum(abs2, UpperTriangular(F.factors), 1)[:]

which will probably traverse memory in more friendly patterns. To test it, I wrote this implementation which can only compute a single rotation given a Cholesky

function _circshiftchol!(C::Cholesky)
n = size(C, 1)
CC = C.factors
if C.uplo == 'U'
# Loop over rows of CC
@inbounds for i in 1:n - 1
# Set active column
k = i + 1
# Compute givens rotation
G, r = givens(CC[i, k], CC[i + 1, k], 1, 2)
# Set value in active column
CC[i, k] = r
CC[i + 1, k] = 0 # not strictly required
# Apply rotation to columns to the right
@simd for j in i + 2:n
cij = CC[i , j]
ci1j = CC[i + 1, j]
CC[i, j] = G.c*cij + G.s*ci1j
CC[i + 1, j] = -G.s*cij + G.c*ci1j
end
# Apply rotation to left column (where lower element is zero)
cij = CC[i, 1]
CC[i , 1] = G.c*cij
CC[i + 1, 1] = -G.s*cij
end
end
return CC
end