Support faster Cholesky Decomposition of Toeplitz matrices

[mike-lawrence]

I recently discovered that cholesky decomposition of toeplitz matrices (useful for GPs evaluated on a grid) can be done with O(n^2) rather than the standard cholesky O(n^3). I found an implementation here, where an academic has posted C, C++, Fortran, MATLAB & Python code (also mirrored in separate repositories here). I took a look at the python and it wasn't particularly optimized (used loops instead of vectorized alternatives), so in case it's useful here's the vectorized version I came up with:

import numpy as np

def cholesky_toeplitz(r):
    """
    Cholesky decomposition of a Toeplitz matrix.
    :param r: the first row of the Toeplitz matrix
    :return: the lower triangular matrix L
    """
    n = len(r)
    l = np.zeros([n, n])
    l[:, 0] = r
    g = np.zeros([2, n])
    g[0, 1:] = r[0:(n-1)]
    g[1, 1:] = r[1:]
    for i in range ( 1, n ):
        rho = - g[1,i] / g[0,i]
        gam = np.sqrt ( ( 1.0 - rho ) * ( 1.0 + rho ) )
        alf = g[0, i:n]
        bet = g[1, i:n]
        temp_g_0 = (alf + rho * bet) / gam
        temp_g_1 = (rho * alf + bet) / gam
        g[0, i:n] = temp_g_0
        g[1, i:n] = temp_g_1
        l[i:n, i] = g[0, i:n]
        g[0, i+1:n] = g[0, i:n-1]
        g[0,i] = 0.0
    return(l)

And in some light testing, it does seem to outperform scipy.linalg.cholesky() when n is large (& n.b. scipy uses LAPACK, so my comparison was native python vs compiled fortran).

[mike-lawrence]

Oh, and Aki just made me aware that the case of decomposing the correlation matrix (one can add the amplitude parameter later by multiplication), it's a symmetric circulant matrix that can be decomposed using FFT methods in O(n*log(n)).

[spinkney]

Here's it written in Stan

      matrix cholesky_toeplitz(vector r) {
        int n = num_elements(r);
        matrix[n, n] L;
        int nm1 = n - 1;
        array[2] row_vector[nm1] g = {r[1:nm1]', r[2:n]'};
        
        L[:, 1] = r; 
        L[1, 2:] = rep_row_vector(0, nm1);
        
        int j = 0;
            
        for (i in 2:n) {
            j += 1;
            real rho = -g[2, j] / g[1, j];
            real gam = sqrt((1 - rho) * (1 + rho));
            row_vector[n - i + 1] alf = g[1, j:nm1];
            row_vector[n - i + 1] bet = g[2, j:nm1];
            g[1, j:nm1] = (alf + rho * bet) / gam;
            g[2, j:nm1] = (rho * alf + bet) / gam;

            L[i:n, i] = to_vector(g[1, j:nm1]);
            
            if (i < n) {
                g[1, i:nm1] = g[1, j:nm1 - 1];
                L[i, i + 1:n] = L[1, 2:n - i + 1];
            }
            g[1, j] = 0;
        }
        return L;
    }