In [181]:

    

import numpy as np

N = 1000
J = 8
alpha = 1.0+np.random.rand(J//2) + 1.0j * np.random.rand(J//2)
alpha = np.append(alpha, np.conj(alpha))
beta = 1.0+np.random.rand(J//2) + 1.0j * np.random.rand(J//2)
beta = np.append(beta, np.conj(beta))
t = np.sort(np.random.uniform(-100.0, 100.0, N))
y0 = np.sin(t)
u = alpha*np.exp(-beta * t[:, None])
v = np.exp(beta * t[:, None])
diag = 0.1 + np.sum(alpha).real + np.zeros(N)

K = np.sum(alpha*np.exp(-beta*np.abs(t[:, None] - t[None, :])[:, :, None]), axis=-1)
K[np.diag_indices_from(K)] = diag

K0 = np.tril(np.dot(u, v.T), -1) + np.triu(np.dot(v, u.T), 1)
K0[np.diag_indices_from(K0)] = diag
print("Semiseparable error: {0}".format(np.max(np.abs(K - K0))))

# Cholesky method
dt = np.diff(t)
phi = np.exp(-beta * dt[:, None])
D = np.empty(N, dtype=alpha.dtype)
X = np.empty((N, J), dtype=alpha.dtype)

# Explicit first step
# D[0] = np.sqrt(diag[0])
D[0] = diag[0]
X[0] = 1.0 / D[0]
S = X[0][:, None] * X[0][None, :] * D[0]

# Then the rest
for n in range(1, N):
    St = phi[n-1][:, None] * phi[n-1][None, :] * S
#     D[n] = np.sqrt(diag[n] - np.sum(alpha[None, :] * alpha[:, None] * St))
    D[n] = diag[n] - np.sum(alpha[None, :] * alpha[:, None] * St)
    X[n] = (1.0 - np.sum(alpha[None, :] * St, axis=1)) / D[n]
    S = St + X[n][:, None] * X[n][None, :] * D[n]

# Check factorization
L = np.tril(np.dot(u, (v*X).T), -1)
L[np.diag_indices_from(L)] = 1.0
# L[np.diag_indices_from(L)] = D

print("Cholesky error: {0}".format(np.max(np.abs(np.dot(L, np.dot(np.diag(D), L.T)) - K))))
print(np.sum(np.log(D)).real - np.linalg.slogdet(K)[1])


    

    




Semiseparable error: 1.3500315867439023e-13
Cholesky error: 2.2026890678642093e-13
-1.36424205266e-12

In [182]:

    

D_orig = np.array(D)


    

In [183]:

    

y = np.array(y0)
z = np.empty(N, dtype=alpha.dtype)
z[0] = y[0] / D[0]
f = 0.0
for n in range(1, N):
    f = phi[n-1] * (f + alpha * X[n-1] * z[n-1] * D[n-1]) 
    z[n] = (y[n] - np.sum(f)) / D[n]
print("Forward sub error: {0}".format(np.max(np.abs(z - np.linalg.solve(L, y) / D))))

y = np.array(z)  # / D
z = np.empty(N, dtype=alpha.dtype)
z[-1] = y[-1]  # / D[-1]
f = 0.0
for n in range(N-2, -1, -1):
    f = phi[n] * (f + alpha * z[n+1]) 
    z[n] = (y[n] - np.sum(f * X[n]))  # / D[n]
print("Backward sub error: {0}".format(np.max(np.abs(z - np.linalg.solve(L.T, y)))))

print("Full solve error: {0}".format(np.max(np.abs(np.linalg.solve(K, y0) - z))))


    

    




Forward sub error: 2.1902621458465837e-14
Backward sub error: 4.05925762865106e-16
Full solve error: 2.1331253615756953e-15

In [ ]:

    




    

In [208]:

    

a = 2*alpha[:J//2].real
b = 2*alpha[:J//2].imag
c = beta[:J//2].real
d = beta[:J//2].imag

tt = t[:, None]
u1 = a*np.exp(-c*tt)*np.cos(d*tt) + b*np.exp(-c*tt)*np.sin(d*tt)
u2 = a*np.exp(-c*tt)*np.sin(d*tt) - b*np.exp(-c*tt)*np.cos(d*tt)
v1 = np.exp(c*tt)*np.cos(d*tt)
v2 = np.exp(c*tt)*np.sin(d*tt)

ut1 = a*np.cos(d*tt) + b*np.sin(d*tt)
ut2 = a*np.sin(d*tt) - b*np.cos(d*tt)
vt1 = np.cos(d*tt)
vt2 = np.sin(d*tt)
dt = np.diff(t)[:, None]
phi = np.exp(-c * dt)

K2 = np.tril(np.dot(u1, v1.T), -1) + np.triu(np.dot(v1, u1.T), 1)
K2 += np.tril(np.dot(u2, v2.T), -1) + np.triu(np.dot(v2, u2.T), 1)
K2[np.diag_indices_from(K2)] = diag
print("Semiseparable error: {0}".format(np.max(np.abs(K.real - K2))))

D = np.empty(N)
X1 = np.empty((N, J//2))
X2 = np.empty((N, J//2))

D[0] = diag[0]
X1[0] = vt1[0] / D[0]
X2[0] = vt2[0] / D[0]
S11 = X1[0, None, :] * X1[0, :, None] * D[0]
S12 = X1[0, :, None] * X2[0, None, :] * D[0]
S22 = X2[0, None, :] * X2[0, :, None] * D[0]

for n in range(1, N):
    S11 *= phi[n-1, :, None] * phi[n-1, None, :]
    S12 *= phi[n-1, :, None] * phi[n-1, None, :]
    S22 *= phi[n-1, :, None] * phi[n-1, None, :]
    
    D[n] = diag[n]
    D[n] -= np.sum(S11 * ut1[n, None, :] * ut1[n, :, None])
    D[n] -= np.sum(S22 * ut2[n, None, :] * ut2[n, :, None])
    D[n] -= 2*np.sum(S12 * ut1[n, :, None] * ut2[n, None, :])
    
    X1[n] = vt1[n]
    X1[n] -= np.sum(S11 * ut1[n, None, :], axis=1)
    X1[n] -= np.sum(S12 * ut2[n, None, :], axis=1)
    X1[n] /= D[n]
    
    X2[n] = vt2[n]
    X2[n] -= np.sum(S22 * ut2[n, None, :], axis=1)
    X2[n] -= np.sum(S12 * ut1[n, :, None], axis=0)
    X2[n] /= D[n]

    S11 += X1[n, None, :] * X1[n, :, None] * D[n]
    S12 += X1[n, :, None] * X2[n, None, :] * D[n]
    S22 += X2[n, None, :] * X2[n, :, None] * D[n]
    
L = np.tril(np.dot(u1, (X1 * np.exp(c*tt)).T), -1) + np.tril(np.dot(u2, (X2 * np.exp(c*tt)).T), -1)
L[np.diag_indices_from(L)] = 1.0

print("Cholesky error: {0}".format(np.max(np.abs(np.dot(L, np.dot(np.diag(D), L.T)) - K2))))
print("Log determinant error: {0}".format(np.sum(np.log(D.astype(np.complex))).real - np.linalg.slogdet(K2)[1]))


    

    




Semiseparable error: 1.3677947663381929e-13
Cholesky error: 2.2382096176443156e-13
Log determinant error: -1.8189894035458565e-12

In [211]:

    

y = np.array(y0)
z = np.empty(N)
z[0] = y[0]
f1 = 0.0
f2 = 0.0
for n in range(1, N):
    f1 = phi[n-1] * (f1 + X1[n-1] * z[n-1])
    f2 = phi[n-1] * (f2 + X2[n-1] * z[n-1])
    z[n] = (y[n] - np.dot(ut1[n], f1) - np.dot(ut2[n], f2))
print("Forward sub error: {0}".format(np.max(np.abs(z - np.linalg.solve(L, y)))))

y = np.array(z) / D
z = np.empty(N)
z[-1] = y[-1]
f1 = 0.0
f2 = 0.0
for n in range(N-2, -1, -1):
    f1 = phi[n] * (f1 + ut1[n+1] * z[n+1])
    f2 = phi[n] * (f2 + ut2[n+1] * z[n+1])
    z[n] = (y[n] - np.dot(X1[n], f1) - np.dot(X2[n], f2))
print("Backward sub error: {0}".format(np.max(np.abs(z - np.linalg.solve(L.T, y)))))

print("Full solve error: {0}".format(np.max(np.abs(np.linalg.solve(K, y0) - z))))


    

    




Forward sub error: 9.825473767932635e-15
Backward sub error: 4.49293380277993e-16
Full solve error: 2.0982564644111967e-15

In [212]:

    

z = np.array(y0)
y = np.empty(N)
y[-1] = z[-1]
f1 = 0.0
f2 = 0.0
for n in range(N-2, -1, -1):
    f1 = phi[n] * (f1 + ut1[n+1] * z[n+1])
    f2 = phi[n] * (f2 + ut2[n+1] * z[n+1])
    y[n] = (z[n] + np.dot(X1[n], f1) + np.dot(X2[n], f2))
print("Backward dot error: {0}".format(np.max(np.abs(y - np.dot(L.T, z)))))

z = np.array(y) * D
y = np.empty(N)
y[0] = z[0]
f1 = 0.0
f2 = 0.0
for n in range(1, N):
    f1 = phi[n-1] * (f1 + X1[n-1] * z[n-1])
    f2 = phi[n-1] * (f2 + X2[n-1] * z[n-1])
    y[n] = (z[n] + np.dot(ut1[n], f1) + np.dot(ut2[n], f2))
print("Forward dot error: {0}".format(np.max(np.abs(y - np.dot(L, z)))))

print("Full dot error: {0}".format(np.max(np.abs(np.dot(K, y0) - y))))


    

    




Backward dot error: 3.4638958368304884e-14
Forward dot error: 5.826450433232822e-13
Full dot error: 9.947598300641403e-14

In [230]:

    

z = np.array(y0)
y = np.empty(N)
y[0] = diag[0]*z[0]
f1 = 0.0
f2 = 0.0
for n in range(1, N):
    f1 = phi[n-1] * (f1 + vt1[n-1] * z[n-1])
    f2 = phi[n-1] * (f2 + vt2[n-1] * z[n-1])
    y[n] = (diag[n] * z[n] + np.dot(ut1[n], f1) + np.dot(ut2[n], f2))

mat = np.tril(np.dot(u1, v1.T), -1) + np.tril(np.dot(u2, v2.T), -1)
mat[np.diag_indices_from(mat)] = diag
print("Forward dot error: {0}".format(np.max(np.abs(np.dot(mat, y0) - y))))

f1 = 0.0
f2 = 0.0
for n in range(N-2, -1, -1):
    f1 = phi[n] * (f1 + ut1[n+1] * z[n+1])
    f2 = phi[n] * (f2 + ut2[n+1] * z[n+1])
    y[n] += np.dot(vt1[n], f1) + np.dot(vt2[n], f2)
    
print("Full dot error: {0}".format(np.max(np.abs(np.dot(K, y0) - y))))


    

    




Forward dot error: 3.268496584496461e-13
Full dot error: 8.526512829121202e-14

In [231]:

    

z = np.array(y0)
y = np.empty(N)
y[-1] = diag[-1]*z[-1]
f1 = 0.0
f2 = 0.0
for n in range(N-2, -1, -1):
    f1 = phi[n] * (f1 + ut1[n+1] * z[n+1])
    f2 = phi[n] * (f2 + ut2[n+1] * z[n+1])
    y[n] = diag[n] * z[n] + np.dot(vt1[n], f1) + np.dot(vt2[n], f2)
    
mat = np.triu(np.dot(v1, u1.T), 1) + np.triu(np.dot(v2, u2.T), 1)
mat[np.diag_indices_from(mat)] = diag
print("Forward dot error: {0}".format(np.max(np.abs(np.dot(mat, y0) - y))))

f1 = 0.0
f2 = 0.0
for n in range(1, N):
    f1 = phi[n-1] * (f1 + vt1[n-1] * z[n-1])
    f2 = phi[n-1] * (f2 + vt2[n-1] * z[n-1])
    y[n] += np.dot(ut1[n], f1) + np.dot(ut2[n], f2)

print("Full dot error: {0}".format(np.max(np.abs(np.dot(K, y0) - y))))


    

    




Forward dot error: 4.405364961712621e-13
Full dot error: 8.526512829121202e-14

In [ ]: