diff options
author | David Luevano <55825613+luevano@users.noreply.github.com> | 2019-12-08 21:54:19 -0700 |
---|---|---|
committer | David Luevano <55825613+luevano@users.noreply.github.com> | 2019-12-08 21:54:19 -0700 |
commit | 7c8b391dddd65ff58417df80dfc79cadddf4f4e6 (patch) | |
tree | 0bad908f5b40848addc3307de4307b1dae4c6e7e /cholesky_solve.py |
First commit
Diffstat (limited to 'cholesky_solve.py')
-rw-r--r-- | cholesky_solve.py | 42 |
1 files changed, 42 insertions, 0 deletions
diff --git a/cholesky_solve.py b/cholesky_solve.py new file mode 100644 index 000000000..b2addc01c --- /dev/null +++ b/cholesky_solve.py @@ -0,0 +1,42 @@ +import numpy as np +from numpy.linalg import cholesky + + +def cholesky_solve(K, y): + """ + Applies Cholesky decomposition to obtain the 'alpha coeficients'. + K: kernel. + y: known parameters. + """ + # The initial mathematical problem is to solve Ka=y. + + # First, add a small lambda value. + K[np.diag_indices_from(K)] += 1e-8 + + # Get the Cholesky decomposition of the kernel. + L = cholesky(K) + size = len(L) + + # Solve Lx=y for x. + x = np.zeros(size) + x[0] = y[0] / L[0, 0] + for i in range(1, size): + temp_sum = 0.0 + for j in range(i): + temp_sum += L[i, j] * x[j] + x[i] = (y[i] - temp_sum) / L[i, i] + + # Now, solve LTa=x for a. + L2 = L.T + a = np.zeros(size) + a_ms = size - 1 + a[a_ms] = x[a_ms] / L2[a_ms, a_ms] + # Because of the form of L2 (upper triangular matriz), an inversion of + # range() needs to be done. + for i in range(0, a_ms)[::-1]: + temp_sum = 0.0 + for j in range(i, size)[::-1]: + temp_sum += L2[i, j] * a[j] + a[i] = (x[i] - temp_sum) / L2[i, i] + + return a |