cholesky.py

"""
Functions for updating Cholesky decompositions and speeding up least-square projections
"""

import numpy as np
from scipy.linalg import cho_solve, solve_triangular


def cho_rk1_update(L, x, inplace_L=True, inplace_x=False):
    """ Cholesky rank one update (in-place)
    """
    if not inplace_L:
        L = np.copy(L)
    if not inplace_x:
        x = np.copy(x)
    n = x.size
    for k in range(n):
        r = np.sqrt(L[k, k]**2 + x[k]**2)
        c = r / L[k, k]
        s = x[k] / L[k, k]
        L[k, k] = r
        if k < n - 1:
            L[k+1:n, k] = (L[k+1:n, k] + s * x[k+1:n]) / c
            x[k+1:n] = c * x[k+1:n] - s * L[k+1:n, k]
    if not inplace_L:
        return L


def cho_rk1_downdate(L, x, inplace_L=True, inplace_x=False):
    """ Cholesky rank one downdate (in-place)
    """
    if not inplace_L:
        L = np.copy(L)
    if not inplace_x:
        x = np.copy(x)
    n = x.size
    for k in range(n):
        r = np.sqrt(L[k, k]**2 - x[k]**2)
        c = r / L[k, k]
        s = x[k] / L[k, k]
        L[k, k] = r
        if k < n - 1:
            L[k+1:n, k] = (L[k+1:n, k] - s * x[k+1:n]) / c
            x[k+1:n] = c * x[k+1:n] - s * L[k+1:n, k]
    if not inplace_L:
        return L


def chol_1d_update_old(L, v):
    """
    Cholesky 1-dimension update

    Update a Cholesky decomposition of the form $A = L \times L^H$ by
    computing the Cholesky decomposition of $(A, v[:-1]; v^H)$, i.e., by adding
    v as a new row and column to A
    """
    n = L.shape[0]
    assert L.shape[1] == n
    assert v.shape[0] == n + 1
    assert v.ndim == 1
    if L is None:
        return np.array([[np.sqrt(np.real(v[0]))]])
    L_new = np.zeros((n+1, n+1))
    L_new[:n, :n] = L
    for j in range(n):
        L_new[-1, j] = (v[j] - np.vdot(L_new[j, :j], L_new[-1, :j])) / L_new[j, j]
    L_new[-1, -1] = np.sqrt(v[-1] - np.sum(np.abs(L_new[-1, :-1])**2))
    return L_new


def chol_1d_update(L, v):
    """
    Cholesky 1-dimension update

    Update a Cholesky decomposition of the form $A = L \times L^H$ by
    computing the Cholesky decomposition of $(A, v[:-1]; v^H)$, i.e., by adding
    v as a new row and column to A
    """
    n = L.shape[0]
    assert L.shape[1] == n
    assert v.shape[0] == n + 1
    assert v.ndim == 1
    if L is None:
        return np.array([[np.sqrt(np.real(v[0]))]])
    L_new = np.zeros((n+1, n+1))
    L_new[:n, :n] = L
    if n != 0:
        L_new[-1, :-1] = solve_triangular(L, v[:-1], lower=True, check_finite=False)
    L_new[-1, -1] = np.sqrt(v[-1] - np.sum(np.abs(L_new[-1, :-1])**2))
    return L_new


def chol_1d_downdate(L, i):
    """
    Cholesky 1-dimension downdate

    Update a Cholesky decomposition of the form $A = L \times L^H$ by
    computing the Cholesky decomposition of the matrix obtained from $A$ by
    removing row $i$ and column $i$.
    """
    L_new = np.delete(np.delete(L, i, axis=0), i, axis=1)
    cho_rk1_update(L_new[i:, i:], L[i+1:, i], inplace_L=True, inplace_x=False)
    return L_new


def chol_1d_downdate_inplace_version(L, atom_indice):
    """
    Cholesky 1-dimension downdate
    """
    # https://en.wikipedia.org/w/index.php?title=Cholesky_decomposition#Adding_and_removing_rows_and_columns
    n = L.shape[0]
    assert L.shape[1] == n

    L_new = np.delete(L, atom_indice, axis=0)
    x = L_new[:, atom_indice]
    L_new = np.delete(L_new, atom_indice, axis=1)
    cho_rk1_update_inplace_version(L_new, x, atom_indice)
    return L_new


def cho_rk1_update_inplace_version(L, x, atom_indice):
    """
    Cholesky rank-1 update
    """
    # https://en.wikipedia.org/w/index.php?title=Cholesky_decomposition#Rank-one_update
    n = x.size
    for k in range(atom_indice, n):
        r = np.sqrt(L[k, k]**2 + x[k]**2)
        c = r / L[k, k]
        s = x[k] / L[k, k]
        L[k, k] = r
        if k < n - 1:
            L[k+1:n, k] = (L[k+1:n, k] + s * x[k+1:n]) / c
            x[k+1:n] = c * x[k+1:n] - s * L[k+1:n, k]