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

            # 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
    return CC

but the simple one-liner seems to be faster.