From 7c8b391dddd65ff58417df80dfc79cadddf4f4e6 Mon Sep 17 00:00:00 2001
From: David Luevano <55825613+luevano@users.noreply.github.com>
Date: Sun, 8 Dec 2019 21:54:19 -0700
Subject: First commit

---
 cholesky_solve.py | 42 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 42 insertions(+)
 create mode 100644 cholesky_solve.py

(limited to 'cholesky_solve.py')

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
-- 
cgit v1.2.3-70-g09d2