# Column reordering of Cholesky using Givens rotations

**URL:** <https://discourse.julialang.org/t/column-reordering-of-cholesky-using-givens-rotations/323>\
**Category:** Numerics\
**Tags:** linearalgebra\
**Created:** [November 15, 2016, 3:42pm UTC](https://discourse.julialang.org/t/column-reordering-of-cholesky-using-givens-rotations/323 "2016-11-15T15:42:34Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![dmbates](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dmbates/32/44_2.png) [@dmbates](https://discourse.julialang.org/u/dmbates)\
**Post date:** [November 15, 2016, 3:42pm UTC](https://discourse.julialang.org/t/column-reordering-of-cholesky-using-givens-rotations/323/1 "2016-11-15T15:42:34Z")

</div>

`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](https://en.wikipedia.org/wiki/Variance_inflation_factor) for linear models. I wanted to check before writing such a function if similar capabilities have already been written.

---

<div class="post-metadata">

**Author:** ![andreasnoack](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasnoack/32/27_2.png) [@andreasnoack](https://discourse.julialang.org/u/andreasnoack)\
**Post date:** [December 12, 2016, 7:26pm UTC](https://discourse.julialang.org/t/column-reordering-of-cholesky-using-givens-rotations/323/2 "2016-12-12T19:26:53Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dmbates](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dmbates/32/44_2.png) [@dmbates](https://discourse.julialang.org/u/dmbates)\
**Post date:** [December 12, 2016, 8:24pm UTC](https://discourse.julialang.org/t/column-reordering-of-cholesky-using-givens-rotations/323/3 "2016-12-12T20:24:16Z")

</div>

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)

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

---

<div class="post-metadata">

**Author:** ![andreasnoack](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasnoack/32/27_2.png) [@andreasnoack](https://discourse.julialang.org/u/andreasnoack)\
**Post date:** [December 13, 2016, 10:00pm UTC](https://discourse.julialang.org/t/column-reordering-of-cholesky-using-givens-rotations/323/4 "2016-12-13T22:00:10Z")

</div>

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

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

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