In [1]:

    

%matplotlib inline
%config InlineBackend.figure_format = "retina"

from __future__ import print_function

from matplotlib import rcParams
rcParams["savefig.dpi"] = 100
rcParams["figure.dpi"] = 100
rcParams["font.size"] = 20


    

In [2]:

    

import numpy as np
import matplotlib.pyplot as plt


    

In [3]:

    

N = 101
x = np.sort(np.random.uniform(0, 10, N))
tau = x[:, None] - x[None, :]
K = np.exp(-0.5 * tau**2) * np.cos(2*np.pi*np.abs(tau))
K[np.diag_indices_from(K)] += 0.1
# K = K[:N//2, N//2:]


    

In [41]:

    

def low_rank_approx(K,
                    start_row=0, n_rows=None,
                    start_col=0, n_cols=None,
                    tol=1.234e-5, max_rank=None,
                    random=None):
    if random is None:
        random = np.random.RandomState()
    
    n_rows = K.shape[0] - start_row if n_rows is None else n_rows
    n_cols = K.shape[1] - start_col if n_cols is None else n_cols
    if max_rank is None:
        max_rank = min(n_rows, n_cols)
    U = np.empty((n_rows, max_rank))
    V = np.empty((n_cols, max_rank))

    rank = 0
    index = np.arange(n_rows)
    norm = 0.0
    finished = False
    rows = slice(start_row, start_row+n_rows)
    cols = slice(start_col, start_col+n_cols)
    
    while True:
        while True:
            if not len(index):
                finished = True
                break
            k = random.randint(len(index))
            i = index[k]
            index[k] = index[-1]
            index = index[:-1]
            
            row = (K[start_row+i, cols] -
                   np.dot(U[i, :rank], V[:, :rank].T))
            j = np.argmax(np.abs(row))
            if np.abs(row[j]) > 1e-14:
                break
                        
        if finished:
            if not rank:
                print("failed")
                if n_cols <= n_rows:
                    return K[rows, cols], np.eye(n_cols)
                return np.eye(n_rows), K[rows, cols].T
            break
        
        row /= row[j]
        V[:, rank] = row
        
        U[:, rank] = col = (K[rows, start_col+j] -
                            np.dot(U[:, :rank], V[j, :rank]))
        
        rank += 1
        if rank >= max_rank:
            break
    
        rowcol_norm = np.dot(row, row) * np.dot(col, col)
        if rowcol_norm < tol**2 * norm:
            break
        
        norm += rowcol_norm
        s = np.dot(U[:, :rank-1].T, col)
        norm += 2 * np.sum(np.abs(s))
        s = np.dot(V[:, :rank-1].T, row)
        norm += 2 * np.sum(np.abs(s))
    
    return U[:, :rank], V[:, :rank]

U, V = low_rank_approx(K, N//2, N-N//2, 0, N//2)


    

In [5]:

    

U.shape, V.shape, K[N//2:, :N//2].shape


    

    
Out[5]:





((51, 17), (50, 17), (51, 50))

In [6]:

    

plt.imshow(K[N//2:, :N//2], interpolation="nearest")


    

    
Out[6]:





<matplotlib.image.AxesImage at 0x10f1abe80>






    

In [7]:

    

plt.imshow(np.dot(U, V.T), interpolation="nearest")


    

    
Out[7]:





<matplotlib.image.AxesImage at 0x113a81da0>






    

In [8]:

    

plt.imshow((K[N//2:, :N//2] - np.dot(U, V.T)), interpolation="nearest")
plt.colorbar()


    

    
Out[8]:





<matplotlib.colorbar.Colorbar at 0x1133ed208>






    

In [ ]:

    




    

In [433]:

    

from scipy.linalg import lu_factor, lu_solve

class Node(object):
    
    def __init__(self,
                 K,               # The matrix
                 start=0,         # The starting index of this node
                 size=None,       # The width of the node
                 tol=1.234e-5,    # Tolerance for approximation
                 min_size=100,    # Minimum leaf width
                 left=None,       # Flag indicating direction of branch
                 parent=None,     # The parent node
                 random=None      # A random number generator
                ):
        self.left = left
        self.parent = parent
        size = K.shape[0] if size is None else size
        self.start = start
        self.size = size
        if size // 2 >= min_size:
            self.U = [None, None]
            self.V = [None, None]
            self.U[1], self.V[0] = low_rank_approx(
                K,
                start_row=start+size//2,
                n_rows=size-size//2,
                start_col=start,
                n_cols=size//2,
                tol=tol, random=random)
            
            # Make copies for applying the inverses of children
            self.U[0] = np.array(self.V[0])
            self.V[1] = np.array(self.U[1])
            self.rank = self.U[0].shape[1]

            self.is_leaf = False
            self.children = (
                Node(K, start, size//2,
                     tol, min_size, random=random,
                     parent=self, left=True),
                Node(K, start+size//2, size-size//2,
                     tol, min_size, random=random,
                     parent=self, left=False),
            )
            
        else:
            self.is_leaf = True
            self.children = tuple()
            self.S = np.array(K[start:start+size, start:start+size])
            
    def compute(self):
        for child in self.children:
            child.compute()
        self.compute_S()
        node = self.parent
        start = self.start
        ind = 0 if self.left else 1
        while node is not None:
            node.U[ind] = self._apply_inverse(node.U[ind], start)
            start = node.start
            ind = 0 if node.left else 1
            node = node.parent
    
    def compute_S(self):
        if not self.is_leaf:
            rank = self.rank
            self.S = np.eye(2*rank)
            self.S[:rank, rank:] = np.dot(self.V[1].T, self.U[1])
            self.S[rank:, :rank] = np.dot(self.V[0].T, self.U[0])

        self.S = lu_factor(self.S, overwrite_a=True)
    
    def __repr__(self):
        return ("Node({0.start}, {0.size})").format(self)
    
    def _apply_inverse(self, x, start):
        start = self.start - start
        if self.is_leaf:
            x[start:start+self.size] = lu_solve(self.S, x[start:start+self.size], overwrite_b=True)
            return x
    
        shape = list(x.shape)
        shape[0] = 2*self.rank
        temp = np.empty(shape)
        s1 = self.children[0].size
        s2 = self.children[1].size

        temp[:self.rank] = np.dot(self.V[1].T, x[start+s1:start+s1+s2])
        temp[self.rank:] = np.dot(self.V[0].T, x[start:start+s1])
        temp = lu_solve(self.S, temp, overwrite_b=True)

        x[start:start+s1] -= np.dot(self.U[0], temp[:self.rank])
        x[start+s1:start+s1+s2] -= np.dot(self.U[1], temp[self.rank:])
        return x
    
    def apply_inverse(self, x):
        for child in self.children:
            x = child.apply_inverse(x)
        return self._apply_inverse(x, 0)
        
    def log_determinant(self):
        logdet = np.sum(np.log(np.abs(np.diag(self.S[0]))))
        for child in self.children:
            logdet += child.log_determinant()
        return logdet


    

In [447]:

    

random = np.random.RandomState(42)
root = Node(K, random=random, min_size=5)
root.compute()


    

In [453]:

    

root.children[0].children


    

    
Out[453]:





(Node(0, 25), Node(25, 25))

In [449]:

    

root.log_determinant(), np.linalg.slogdet(K)[1]


    

    
Out[449]:





(-165.30233035504668, -165.3023383604226)

In [450]:

    

x0 = np.random.randn(N, 5)


    

In [451]:

    

x = np.array(x0)
root.apply_inverse(x) - np.linalg.solve(K, x0)


    

    
Out[451]:





array([[ -1.78608359e-06,  -3.27255662e-06,   1.24005590e-05,
          2.20071992e-07,   6.64210257e-06],
       [ -1.57363594e-05,  -1.10600048e-05,  -2.58190304e-05,
          6.53466572e-06,  -2.42662167e-05],
       [ -1.39367082e-05,  -9.12023448e-06,  -4.12748314e-05,
          7.02454800e-06,  -3.18533311e-05],
       [ -9.86236281e-06,  -8.27448694e-06,  -2.69307123e-05,
          1.57924312e-06,  -1.92856083e-05],
       [  8.29791446e-06,   7.64493485e-07,   1.91371415e-05,
         -6.67502317e-06,   1.92562417e-05],
       [  1.36779326e-05,   4.58143211e-06,   2.88335966e-05,
         -5.32764574e-06,   2.67338361e-05],
       [  1.88581488e-05,   8.93897989e-06,   3.63039512e-05,
         -1.81965036e-06,   3.18880901e-05],
       [  1.06177061e-05,   4.05441309e-06,   2.59217741e-05,
         -1.25362217e-06,   2.19066180e-05],
       [ -9.72044926e-07,  -5.66840206e-06,   1.55273507e-05,
         -9.25543623e-06,   1.39564488e-05],
       [ -2.85433315e-06,  -7.29192445e-06,   1.39748428e-05,
         -1.07169650e-05,   1.28264812e-05],
       [ -6.61823747e-05,  -3.56298763e-05,  -1.45437374e-04,
          3.88789098e-05,  -1.30603061e-04],
       [ -1.76862577e-05,  -1.77002433e-05,  -2.16902712e-05,
          9.62538696e-06,  -2.16973043e-05],
       [  5.31599273e-05,   1.12535059e-05,   1.33110939e-04,
         -3.42606065e-05,   1.24331478e-04],
       [  8.65405495e-06,  -4.40207004e-06,   5.50621445e-06,
         -1.75541670e-05,   2.20050719e-05],
       [ -1.82257859e-05,  -1.97280797e-05,  -6.71747646e-05,
          1.40313886e-06,  -4.28589429e-05],
       [ -6.62104088e-05,  -2.95303178e-05,  -1.64937883e-04,
          3.78109252e-05,  -1.48556670e-04],
       [  2.54727837e-05,   2.23348502e-05,   7.20642975e-05,
         -1.96925260e-05,   6.43394606e-05],
       [  1.36366557e-05,   2.01320536e-05,   1.01915436e-05,
         -1.45982504e-05,   1.36777044e-05],
       [  1.70312627e-06,   1.41639268e-05,  -2.05148784e-05,
         -7.83824687e-06,  -1.46965584e-05],
       [  6.42742369e-06,   2.67686444e-05,  -6.96026353e-05,
         -9.38852365e-06,  -4.54593212e-05],
       [  5.51659170e-05,   6.40070708e-05,   6.76089531e-05,
         -2.99442221e-05,   8.35161472e-05],
       [ -5.02696666e-06,   3.46051439e-05,  -1.64046715e-05,
          1.10570166e-05,  -3.53594067e-05],
       [ -2.55200188e-05,   2.21437225e-05,  -2.26098293e-05,
          3.07688507e-05,  -7.01772074e-05],
       [ -9.14315157e-06,   2.84004699e-05,   2.71773981e-05,
          3.04769330e-05,  -2.38987283e-05],
       [ -2.31100617e-05,   2.11874575e-05,   4.60345918e-05,
          1.83734536e-05,  -3.27420969e-06],
       [ -1.94589218e-05,   7.19682213e-06,   4.93749478e-06,
         -3.87761039e-06,  -2.74213811e-06],
       [  1.50131886e-04,   9.59031241e-05,  -2.49435006e-05,
         -7.17224487e-05,   6.46174211e-05],
       [ -1.07543359e-05,   7.40731696e-06,  -2.65109879e-06,
         -1.04284279e-05,  -1.08025565e-05],
       [ -8.73083381e-06,   3.52590971e-05,  -1.11144697e-05,
          2.53181812e-05,  -1.57719920e-05],
       [ -1.38292100e-05,   9.75772476e-05,   1.27009551e-04,
          1.28392538e-04,  -6.38523367e-05],
       [  1.59179781e-05,   2.11815385e-05,   9.76146252e-06,
         -3.75694207e-06,   1.24249012e-05],
       [  2.76683340e-05,   3.83298825e-05,  -1.64771406e-05,
         -1.42121399e-05,   2.02345333e-05],
       [  2.17721661e-05,   2.53362178e-05,   2.24771696e-05,
          3.15204374e-05,  -4.31167524e-06],
       [  1.84445550e-05,   1.67082278e-05,   2.16713412e-05,
          3.74645531e-05,  -4.01890214e-06],
       [  1.51599836e-05,  -6.97755537e-06,   2.61940017e-05,
          3.65928788e-05,  -6.45016107e-06],
       [ -7.08851224e-06,  -3.23937307e-05,  -1.06570388e-05,
         -7.48985421e-06,   1.20213235e-05],
       [  9.48258970e-05,   3.14713108e-05,  -7.61827868e-06,
          1.34439560e-05,  -5.71996390e-06],
       [  4.49629415e-06,   1.35997252e-05,  -5.21436807e-06,
         -6.31914696e-06,   1.18954795e-05],
       [  8.98326555e-06,   2.00782119e-05,  -1.25105802e-05,
          2.52122335e-06,   1.29232011e-05],
       [ -1.60796502e-05,  -8.74914104e-05,   7.55429970e-05,
          2.62127451e-06,  -7.56508666e-06],
       [  7.28938007e-06,   1.70631804e-05,  -9.31463750e-06,
         -7.79747262e-06,   1.04360726e-06],
       [  5.98075498e-06,   1.63850005e-05,  -7.76458248e-06,
         -1.03012829e-05,  -1.63220207e-06],
       [ -6.45130704e-06,   1.53760844e-05,   6.13171652e-06,
         -2.73315631e-05,  -1.91862293e-05],
       [  8.90097262e-05,   3.22852247e-05,   2.16427339e-05,
         -7.42081055e-05,   1.89242709e-05],
       [ -5.49743748e-05,   1.91472777e-05,   4.36413985e-05,
         -1.95606664e-05,  -3.82520547e-05],
       [  6.34311392e-06,   2.03723855e-05,   2.51161030e-05,
         -3.70880477e-05,  -1.17570017e-06],
       [ -1.64850573e-06,  -1.01129005e-05,   4.20448146e-06,
          8.90292440e-06,   3.85006043e-06],
       [  2.09892494e-05,  -3.15421325e-05,  -2.01609370e-05,
          3.31973522e-05,   1.95547303e-05],
       [  1.75713424e-05,  -2.93786827e-05,  -1.45702082e-05,
          3.03264188e-05,   2.07088247e-05],
       [  1.75482364e-05,  -2.93091088e-05,  -1.45419167e-05,
          3.02860953e-05,   2.06876046e-05],
       [  5.54104513e-05,   1.22017439e-05,  -8.99237692e-05,
          2.08009920e-05,   2.96419326e-05],
       [ -2.03388230e-04,   2.51412699e-04,   2.95273187e-04,
         -3.37230761e-04,  -1.99465858e-04],
       [ -9.78384958e-05,   1.08876847e-04,   1.41750250e-04,
         -1.46684523e-04,  -8.64335359e-05],
       [ -5.07618299e-05,   4.93581522e-05,   9.02395214e-05,
         -6.66891549e-05,  -1.05452707e-05],
       [  1.08042712e-04,  -1.27700154e-04,  -1.51405557e-04,
          1.65060468e-04,   8.34532913e-05],
       [  1.25007883e-04,  -1.43494916e-04,  -1.75572513e-04,
          1.86606478e-04,   9.37276674e-05],
       [ -3.31708495e-05,   7.46994082e-05,   4.05813337e-05,
         -7.63351690e-05,  -4.60186440e-05],
       [ -2.34735729e-04,   3.14860760e-04,   3.09794446e-04,
         -3.80153503e-04,  -1.98263039e-04],
       [ -1.67950503e-04,   2.28087587e-04,   2.14837629e-04,
         -2.71539251e-04,  -1.42309114e-04],
       [ -3.08671773e-06,   1.33794655e-05,  -2.11643356e-05,
         -4.46743433e-07,  -1.54500012e-06],
       [  2.09324017e-04,  -2.63583132e-04,  -3.27519188e-04,
          3.47566470e-04,   1.75428694e-04],
       [ -9.35090261e-07,   8.97328067e-05,   1.40052634e-05,
         -5.87027665e-05,  -1.49317102e-05],
       [ -1.84808924e-04,   1.99303347e-04,   1.69634969e-04,
         -2.56421216e-04,  -1.38184916e-04],
       [  4.93008028e-05,  -5.34377080e-05,  -1.11416177e-04,
          9.20152158e-05,   5.19251030e-05],
       [  5.21444061e-05,  -5.59692198e-05,  -1.02217163e-04,
          9.49134104e-05,   5.72720147e-05],
       [  2.06599311e-05,  -1.64717860e-05,  -4.49413695e-05,
          4.11113380e-05,   3.38612007e-05],
       [  2.01859294e-05,  -1.59116188e-05,  -4.39236608e-05,
          4.03586385e-05,   3.33197524e-05],
       [ -3.10960305e-05,   4.78992900e-05,   4.17659157e-05,
         -4.64587721e-05,  -7.05067829e-06],
       [ -1.28905514e-04,   1.69258270e-04,   2.00946154e-04,
         -2.12224975e-04,  -8.55166344e-05],
       [  4.23107312e-05,  -4.93761947e-05,   9.50243411e-05,
          6.91899687e-06,   2.97785768e-05],
       [  6.25641119e-05,  -7.90069367e-05,  -5.18838464e-07,
          7.80633199e-05,   5.85461063e-05],
       [ -7.89157820e-05,   9.54700932e-05,   1.76127121e-04,
         -1.40683719e-04,  -4.39631482e-05],
       [ -8.60234030e-05,   1.04250779e-04,   1.82535791e-04,
         -1.50376676e-04,  -4.82833227e-05],
       [  3.80597157e-05,  -4.81287953e-05,  -3.31687143e-05,
          7.35242639e-05,   6.14977146e-05],
       [ -9.64290900e-06,   1.97468112e-05,  -1.24066700e-04,
          5.86610902e-05,   7.37526098e-05],
       [  6.82637258e-05,  -2.66734535e-04,   8.24504222e-04,
         -9.32397798e-05,  -4.16023076e-04],
       [ -1.98615452e-05,   4.50173627e-05,  -1.49952278e-04,
          3.54752789e-05,   9.00448319e-05],
       [ -2.05672862e-05,   4.71943303e-05,  -1.57321602e-04,
          3.64050745e-05,   9.32012557e-05],
       [ -2.40199890e-05,   5.06512949e-05,  -1.76253531e-04,
          3.42797891e-05,   6.92968191e-05],
       [ -1.01105034e-05,   4.29175714e-06,  -2.40543422e-05,
          8.36908092e-06,  -2.52524872e-05],
       [  6.50260701e-05,  -2.13003381e-04,   7.42763513e-04,
         -1.24643071e-04,  -4.34944089e-04],
       [ -1.05620531e-06,  -1.67990165e-04,   4.18540274e-04,
          4.24179655e-05,  -3.66362035e-04],
       [ -3.07929857e-05,  -1.13324029e-04,   1.37308890e-04,
          1.19747529e-04,  -2.36026512e-04],
       [ -1.09787285e-05,  -9.14346457e-05,   7.17830633e-05,
          9.32660170e-05,  -1.81826580e-04],
       [ -1.03918627e-05,  -9.10145364e-05,   7.10258352e-05,
          9.22835261e-05,  -1.81032454e-04],
       [ -9.44037166e-06,  -9.03357669e-05,   6.98288671e-05,
          9.06801629e-05,  -1.79761964e-04],
       [  6.62441600e-05,  -2.98410041e-05,   2.62677501e-05,
         -6.63063653e-05,  -6.45218758e-05],
       [  1.00079830e-04,   2.93470040e-06,   6.97409657e-06,
         -1.44652543e-04,  -6.38491608e-06],
       [  1.39420654e-04,   4.45515119e-05,  -1.61929430e-05,
         -2.39711306e-04,   6.59623085e-05],
       [  1.79174125e-04,   9.42206518e-05,  -3.94491781e-05,
         -3.27073497e-04,   1.33618747e-04],
       [ -7.29691269e-04,  -6.57612054e-04,   3.76105226e-04,
          1.65799326e-03,  -1.22186976e-03],
       [  1.62853238e-04,   6.92221578e-05,  -3.50613052e-05,
         -3.11612004e-04,   1.25673732e-04],
       [  1.54269136e-04,   6.21199619e-05,  -3.39973072e-05,
         -2.90667378e-04,   1.12698707e-04],
       [  3.37537183e-05,  -3.16396975e-08,  -1.07422020e-05,
         -3.58974709e-05,  -6.67997992e-06],
       [ -3.03166688e-06,  -1.04276258e-05,   1.95635158e-06,
          2.39532606e-05,  -2.06142330e-05],
       [ -1.55972051e-04,  -1.43227454e-05,   7.75965882e-05,
          2.09490102e-04,   7.37616529e-06],
       [ -1.89937725e-04,  -6.82732312e-06,   9.84002728e-05,
          2.39355177e-04,   3.00840290e-05],
       [ -5.87283668e-04,   3.79256906e-04,   4.08499230e-04,
          2.20389771e-04,   8.89752014e-04],
       [ -1.69133332e-04,   1.21702111e-04,   1.05144297e-04,
          5.27545965e-05,   2.86830695e-04],
       [  4.27646813e-05,  -1.33334602e-05,  -4.56818135e-05,
         -2.15556567e-05,  -3.48483795e-05],
       [ -6.05265747e-04,   4.40601979e-04,   4.45386983e-04,
          1.60131729e-04,   1.01303033e-03]])

In [452]:

    

np.linalg.solve(K, x0)


    

    
Out[452]:





array([[ -1.49098277e+01,   7.65754720e+00,   1.25192136e+01,
          7.64545796e-01,   1.22362483e+01],
       [  1.73264167e+01,  -3.12039006e+00,  -2.00608811e+01,
         -6.00275288e+00,  -1.16559814e+01],
       [ -9.80813225e+00,  -1.20992008e+01,  -1.82160047e+00,
         -1.67143469e+00,   2.21550800e+01],
       [  1.77431158e+00,  -4.62247065e+00,   3.13637007e+00,
         -2.01924526e+00,  -8.23517286e+00],
       [  1.60298306e+01,   5.00124204e+00,   1.33704046e+01,
          1.10980256e+01,  -4.20627993e+00],
       [ -3.82900075e+00,   2.77115547e+00,   1.02196022e-02,
          1.06011391e+01,  -7.80628323e+00],
       [ -6.23062510e+00,  -1.63690047e+01,   2.74999510e+00,
          2.26600579e+00,   3.05889135e+00],
       [  6.92416062e+00,   1.58740771e+01,   2.84550091e+00,
          4.29720667e+00,   8.72652502e-02],
       [ -1.52389091e+01,   2.77249883e+00,  -2.56473498e+01,
         -8.15763323e+00,  -4.04124341e+00],
       [ -5.24129621e-01,   1.02581431e+01,   1.82878609e+00,
         -1.92704057e+01,   8.56080507e+00],
       [ -2.34763513e+00,  -1.47444588e+01,  -7.90714275e-01,
          7.83657381e-01,  -2.47728709e+00],
       [  1.39787408e+01,  -1.45460525e+01,   3.14555692e+00,
         -5.15007995e+00,   1.02987765e+01],
       [  1.87753794e+00,   5.54964911e+00,  -6.36022691e+00,
          1.18214137e+01,  -1.94674276e+00],
       [ -1.10433043e+01,   1.23462561e+00,  -4.41871982e+00,
         -4.00436770e+00,   1.06594029e+01],
       [ -6.37539306e+00,  -1.05005284e+01,  -3.41953118e+00,
         -1.23456880e+01,  -1.09683583e+01],
       [ -7.81773803e-01,  -5.90100519e+00,  -4.71539855e+00,
         -7.02314911e+00,   5.63508421e+00],
       [ -6.79607312e-01,  -2.02090415e+01,   3.54594827e+00,
         -1.26749122e+01,   4.20546525e+00],
       [ -7.33785953e+00,  -2.24347762e+00,   7.22213011e+00,
         -1.36601889e+00,   5.44595322e+00],
       [ -1.27688295e-01,  -5.85055517e+00,  -1.20199482e+01,
         -1.42414709e+01,   5.38882654e+00],
       [ -6.41296378e+00,  -1.15527750e+01,   7.35893493e+00,
         -8.26866722e-01,  -1.11153224e+01],
       [  8.43294772e-02,  -6.52328789e+00,   2.72351228e+01,
         -1.25402722e+01,  -3.49074478e+00],
       [ -1.62107945e+01,   5.85193761e+00,  -1.08024034e+01,
         -7.08098499e+00,   1.94506807e+01],
       [  1.64462324e+01,   5.25568829e+00,  -5.52927342e-01,
         -1.01673842e+01,   9.35176680e-01],
       [  2.64206277e+00,  -6.27259657e+00,   2.78331519e+00,
          8.29663417e+00,   9.00437287e+00],
       [  2.41194209e+00,  -2.31084922e+00,   2.78305709e+01,
         -6.19119038e+00,  -6.19817582e+00],
       [ -2.02349427e+00,   1.21632121e+01,  -1.55861630e+01,
          5.63228213e+00,   4.01016640e+00],
       [ -9.19380818e+00,  -1.62609568e+00,   3.35429085e-01,
         -1.12057495e+00,   2.39947635e+00],
       [  9.86799135e+00,  -6.64823754e+00,   2.42711560e+00,
         -6.10223206e+00,  -1.03232084e+01],
       [ -5.59479014e+00,   4.17072541e+00,   2.49701382e+01,
         -2.65799997e+00,   2.15504112e+00],
       [  1.00070214e+01,   2.68496620e+00,   1.66819728e+01,
         -1.96970116e+00,   1.95265402e+01],
       [ -5.75278407e-01,   9.65217507e+00,   1.39434355e+01,
         -1.01018322e+01,  -6.75705881e+00],
       [  5.86729955e+00,  -3.68862172e+00,   3.58468569e+00,
          9.56335409e+00,  -1.92818230e-01],
       [  6.88615954e-01,   9.90901964e+00,   4.31872786e+00,
         -1.23691538e+01,   3.72646785e+00],
       [ -1.13824990e+01,  -1.82421637e+01,   1.13372309e+01,
          1.11414191e+01,   1.10534571e+00],
       [  5.83014869e+00,  -4.88154401e-01,   1.74443303e+00,
         -5.00230968e+00,   1.02521818e+00],
       [  8.16001813e-01,  -1.26748269e+01,  -8.74758260e+00,
         -4.89085319e+00,   4.24975117e+00],
       [ -1.52401980e+01,  -8.32245263e+00,   2.75314904e+00,
          5.66799539e+00,  -6.23914768e+00],
       [  1.57285715e+00,  -6.15208347e+00,   8.96371219e+00,
         -9.93070788e+00,  -2.51386413e+00],
       [ -1.12588551e+01,   3.05341497e+00,  -1.48790606e+00,
          1.21141715e+01,   2.93738325e+00],
       [  1.60487188e+01,  -7.36902783e-01,  -5.66112449e+00,
         -9.10539610e-01,   8.17319135e+00],
       [  1.10123034e+01,  -5.45365450e+00,  -2.68704905e+00,
         -2.04960031e+00,  -1.02781128e+01],
       [ -1.19185050e+01,  -4.16321198e+00,   1.23445691e+01,
         -9.79701006e+00,   1.06787407e+01],
       [ -6.15841506e+00,  -1.00948779e+01,  -1.41633780e+01,
         -3.08161463e+00,  -9.82100437e-01],
       [ -2.58698215e+00,  -9.26186892e+00,   7.65809273e+00,
         -6.58916487e-01,  -1.35369830e+00],
       [  1.93881035e+00,  -5.85636278e-01,  -5.13177036e+00,
          2.30765235e+00,   3.30716265e+00],
       [  4.41553359e+00,  -4.50596258e+00,   2.80411431e+00,
          9.21826480e+00,  -5.68673545e+00],
       [ -2.56707451e+00,  -1.54466669e+01,   6.73945153e+00,
          2.63677988e+00,  -1.34780489e+00],
       [  8.99219168e+00,   1.16499608e+00,   8.10078684e+00,
         -9.34467283e+00,  -1.00075264e+01],
       [ -1.36289128e+00,   5.16413279e-01,  -1.03051487e+01,
          1.94498571e+01,  -2.77195306e+00],
       [  6.07100098e-01,   9.98501090e+00,  -8.73855633e+00,
         -4.02153422e+00,   8.84130707e+00],
       [ -5.34778956e+00,  -1.16592536e+01,   1.32908523e+00,
          3.44225170e+00,   1.92505555e+00],
       [  8.91259556e+00,   5.87543941e+00,   3.23442444e+00,
         -6.97635796e+00,  -2.76831332e+00],
       [  1.12379377e+01,   5.34640348e+00,   1.10720193e+01,
         -7.60222261e+00,   2.25934718e+00],
       [ -9.16142130e+00,  -1.38379906e+01,  -5.59282698e+00,
          2.24186541e+00,  -5.67751558e+00],
       [  3.32427512e+00,   4.20948467e+00,   2.92125264e+00,
         -1.47730067e+01,   5.25845851e+00],
       [ -6.83168489e+00,  -1.19650521e+01,   5.09111461e+00,
          1.79079784e+00,   2.60482420e+00],
       [  1.29910585e+01,   2.31148024e+00,  -4.82821899e+00,
          8.58116287e+00,  -5.36556963e+00],
       [ -1.27553634e+01,  -4.31694453e-01,   1.39881139e+00,
         -3.10516277e+00,  -1.01280969e+00],
       [  4.00448177e+00,  -1.52235836e+00,  -1.36390950e+01,
          8.24433941e-01,   1.42050618e+01],
       [  1.23611946e+01,   1.27009483e+01,  -1.54438607e-01,
         -9.51647495e+00,  -4.20286783e+00],
       [  7.15691531e+00,  -1.58071490e+01,   1.63507066e+01,
         -1.15302204e+01,   3.50593912e+00],
       [ -6.01855187e-02,   8.88819834e-01,  -4.04210490e+00,
         -1.46017192e+00,  -8.74334962e+00],
       [ -4.57077468e+00,   2.55310431e+00,   2.10984225e+00,
          5.46887832e+00,   1.15456112e+01],
       [  2.07341542e+00,   4.91963429e+00,   4.47713232e+00,
         -3.95084074e+00,  -1.05025062e+01],
       [  1.10955423e+01,   1.48696074e+01,  -4.60078257e+00,
         -3.24573237e+00,   1.17281324e+01],
       [ -7.30601297e+00,  -3.54564839e+00,  -5.81524796e+00,
          1.87972371e+00,   3.90595408e+00],
       [ -5.40424643e+00,  -5.49232508e+00,   1.28925570e+01,
         -7.48625545e+00,  -1.04102047e+01],
       [  1.39307393e+01,  -8.93744891e+00,  -1.17993733e+01,
         -1.54024348e+01,   3.36904264e+00],
       [ -1.15353283e+01,  -8.77044599e+00,  -1.46596001e+00,
          5.79206888e+00,   2.58081642e+00],
       [ -1.01897426e+00,  -1.33368664e+00,   1.89607650e+01,
          2.57746353e+00,  -1.54530599e+01],
       [ -1.22723753e+01,  -4.42298615e+00,   5.37786924e-01,
          7.93535768e-02,   8.26226256e-01],
       [  1.59469951e+01,  -3.93242420e+00,  -1.39354198e+01,
          6.24776612e+00,   1.19068040e+01],
       [ -1.24434590e+01,  -5.02778630e+00,  -5.23006820e+00,
         -1.81618472e+01,  -2.70843594e+00],
       [  6.23093349e+00,   1.03245016e+01,   3.83071433e+00,
         -4.08503747e+00,  -2.04476813e+01],
       [ -4.01540889e+00,  -1.70275217e+01,   1.67146705e+01,
          2.13117187e+01,  -1.63675634e+01],
       [  5.37087287e+00,  -1.14593133e+01,   6.90395464e+00,
          2.47341128e+00,   1.48986560e+01],
       [  6.11349813e+00,   7.12791740e+00,  -1.15005446e+01,
         -5.26841942e+00,  -4.41651515e-01],
       [  2.19485128e+00,   1.55550940e+01,   4.79028444e+00,
         -6.96100896e+00,   1.37100962e-01],
       [  2.38100821e-01,  -3.27140798e+00,   1.36533494e+01,
          5.92832141e+00,  -1.28838791e+00],
       [ -9.52045869e+00,  -1.04537740e+01,   7.04792204e+00,
         -2.05483210e+01,  -7.99875402e+00],
       [  3.06907633e+00,   3.90039579e+00,  -2.61286191e+01,
          1.38711600e+01,  -7.92809311e+00],
       [  4.27621429e+00,  -1.89343075e-02,   1.24146725e+01,
          4.80227789e+00,   3.51102918e+00],
       [  3.32695379e+00,  -1.61614912e+01,  -7.82264754e+00,
         -6.02616743e+00,  -6.36179361e+00],
       [  7.54545941e+00,  -1.08163634e+00,   1.64921194e+01,
          5.50746561e+00,   2.86673317e+00],
       [ -8.37499158e+00,   6.32300130e+00,   4.41429763e+00,
          1.21406625e+01,  -6.49115215e+00],
       [ -4.20885079e+00,   1.06487280e+01,  -9.16959575e+00,
         -1.34141879e+01,   2.87429432e+00],
       [  4.55322575e+00,  -5.42462185e+00,  -2.19389539e+01,
          5.53908390e-01,   1.09635314e+01],
       [ -7.20194792e+00,   7.37902892e+00,   5.66289037e+00,
          6.67981348e+00,  -6.07223740e+00],
       [ -4.87387486e-01,  -2.65615272e+00,  -2.09024551e+00,
          6.06901767e+00,  -1.32300301e+01],
       [  5.91711846e+00,  -7.65686734e-01,   1.20483017e+01,
         -8.33972875e+00,   1.82128435e+01],
       [ -9.94251062e+00,   6.82033052e+00,   7.41381711e+00,
          3.10017238e+00,   1.57997180e+01],
       [  6.25290905e+00,  -5.32363597e+00,  -1.55661208e+01,
          3.59354721e+00,  -9.10383684e-01],
       [ -4.47855347e+00,   1.47653984e+00,   1.39390671e+00,
          8.69877750e+00,  -1.30640960e+01],
       [ -8.27514907e+00,   2.84356424e+00,   7.64382670e+00,
         -7.27789607e+00,  -2.44438367e+00],
       [  4.20649706e-01,  -8.42124920e+00,  -1.49917818e+00,
          3.55170747e+00,   5.73842154e+00],
       [ -1.46509500e+01,   9.76900608e+00,  -2.42189843e+00,
         -3.96563919e-01,  -3.13424921e+00],
       [  8.04883303e+00,  -2.86628028e+00,   2.33993264e+00,
          1.89224236e+01,   8.66747889e+00],
       [ -1.46970381e+01,  -9.84860627e+00,   1.01740098e+01,
          1.72379030e+01,  -2.83726184e+00],
       [ -2.97748127e+00,   9.89478961e-01,  -1.77364564e+01,
         -4.07554252e+00,  -8.89218404e+00],
       [  8.32317168e+00,   2.94005077e+00,   1.19028909e+01,
          5.71255728e+00,   4.46370568e+00],
       [ -3.73537229e+00,  -3.22672069e+00,   1.47281284e+00,
          9.21443241e+00,  -6.89751749e+00]])

In [291]:

    

root.V[0][0]


    

    
Out[291]:





array([ -3.37478176e-14,   1.94577600e-11,  -7.24225978e-09,
        -3.00535730e-10,   2.35734213e-07,   7.01111040e-07,
         1.26845259e-06,  -2.38237772e-06,   1.65330121e-05,
         8.30605143e-06,   1.24076079e-05,  -2.13668341e-04,
         2.07632691e-05,   2.57517049e-04,   2.94201801e-03,
         3.63567685e-03,   4.54317683e-04])

In [90]:

    

root.compute_K()


    

In [91]:

    

np.linalg.slogdet(K)


    

    
Out[91]:





(1.0, -165.3023383604226)

In [92]:

    




    

In [94]:

    




    

In [15]:

    

root.children[0].children[0].children[0].children[0].log_determinant()


    

    
Out[15]:





-6.8880219216076322

In [331]:

    

np.sum(np.log(np.abs(np.diag(lu_factor(K)[0]))))


    

    
Out[331]:





-163.67731577167658

In [332]:

    

np.linalg.slogdet(K)


    

    
Out[332]:





(1.0, -163.67731577167658)

In [ ]: