Source code for GPy.util.choleskies

# Copyright James Hensman and Max Zwiessele 2014, 2015
# Licensed under the BSD 3-clause license (see LICENSE.txt)

import numpy as np
from . import linalg
from .config import config

    from . import choleskies_cython
    use_choleskies_cython = config.getboolean('cython', 'working')
except ImportError:
    print('warning in choleskies: failed to import cython module: falling back to numpy')
    use_choleskies_cython = False

[docs]def safe_root(N): i = np.sqrt(N) j = int(i) if i != j: raise ValueError("N is not square!") return j
def _flat_to_triang_pure(flat_mat): N, D = flat_mat.shape M = (-1 + safe_root(8*N+1))//2 ret = np.zeros((D, M, M)) for d in range(D): count = 0 for m in range(M): for mm in range(m+1): ret[d,m, mm] = flat_mat[count, d]; count = count+1 return ret def _flat_to_triang_cython(flat_mat): N, D = flat_mat.shape M = (-1 + safe_root(8*N+1))//2 return choleskies_cython.flat_to_triang(flat_mat, M) def _triang_to_flat_pure(L): D, _, M = L.shape N = M*(M+1)//2 flat = np.empty((N, D)) for d in range(D): count = 0; for m in range(M): for mm in range(m+1): flat[count,d] = L[d, m, mm] count = count +1 return flat def _triang_to_flat_cython(L): return choleskies_cython.triang_to_flat(L) def _backprop_gradient_pure(dL, L): """ Given the derivative of an objective fn with respect to the cholesky L, compute the derivate with respect to the original matrix K, defined as K = LL^T where L was obtained by Cholesky decomposition """ dL_dK = np.tril(dL).copy() N = L.shape[0] for k in range(N - 1, -1, -1): for j in range(k + 1, N): for i in range(j, N): dL_dK[i, k] -= dL_dK[i, j] * L[j, k] dL_dK[j, k] -= dL_dK[i, j] * L[i, k] for j in range(k + 1, N): dL_dK[j, k] /= L[k, k] dL_dK[k, k] -= L[j, k] * dL_dK[j, k] dL_dK[k, k] /= (2 * L[k, k]) return dL_dK
[docs]def triang_to_cov(L): return np.dstack([[:,:,i], L[:,:,i].T) for i in range(L.shape[-1])])
[docs]def multiple_dpotri(Ls): return np.array([linalg.dpotri(np.asfortranarray(Ls[i]), lower=1)[0] for i in range(Ls.shape[0])])
[docs]def indexes_to_fix_for_low_rank(rank, size): """ Work out which indexes of the flatteneed array should be fixed if we want the cholesky to represent a low rank matrix """ #first we'll work out what to keep, and the do the set difference. #here are the indexes of the first column, which are the triangular numbers n = np.arange(size) triangulars = (n**2 + n) / 2 keep = [] for i in range(rank): keep.append(triangulars[i:] + i) #add the diagonal keep.append(triangulars[1:]-1) keep.append((size**2 + size)/2 -1)# the very last element keep = np.hstack(keep) return np.setdiff1d(np.arange((size**2+size)/2), keep)
if use_choleskies_cython: triang_to_flat = _triang_to_flat_cython flat_to_triang = _flat_to_triang_cython backprop_gradient = choleskies_cython.backprop_gradient_par_c else: backprop_gradient = _backprop_gradient_pure triang_to_flat = _triang_to_flat_pure flat_to_triang = _flat_to_triang_pure