# Column reordering of Cholesky using Givens rotations

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

I donâ€™t know why I havenâ€™t seen this question earlier. I think that `LinAlg.lowrankupdate!` and `LinAlg.lowrankdowndate!` might do this.

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.

1 Like

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
``````

but the simple one-liner seems to be faster.