Column reordering of Cholesky using Givens rotations

Specific Domains Numerics
1

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.

2

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

3

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
4

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.