from __future__ import division, absolute_import, print_function
import numpy as np
def form_a(n, m):
a = np.zeros((n, n))
for i in range(1, n+1):
for j in range(1, n+1):
if i == j:
a[i-1][j-1] = n + m**2 + j/m + i/n
else:
a[i-1][j-1] = (i + j)/(m + n)
return a
def form_b(n):
b = np.zeros(n)
for i in range(1,n+1):
b[i-1] = 200 + 50*i
return b
def p_print(x):
return str(x.tolist()).replace(']','').replace('[', '').replace(',','')
def cholesky(A):
n = len(A)
L = np.zeros((n, n))
for i in range(n):
for j in range(i + 1):
t_sum = 0
for k in range(j):
t_sum += L[i][k] * L[j][k]
if i == j:
L[i][j] = np.sqrt(A[i][i] - t_sum)
else:
L[i][j] = 1.0 / L[j][j] * (A[i][j] - t_sum)
return L
def solve_ls(L, b):
n = len(L)
y = np.matrix([0.0 for i in range(n)]).T
for i in range(n):
y[i, 0] = (b[i] - L[i, :i] * y[:i, 0]) / L[i, i]
L = L.T
x = np.matrix([0.0 for i in range(n)]).T
for i in range(n - 1, -1, -1):
x[i, 0] = (y[i] - L[i, i:] * x[i:, 0]) / L[i, i]
return p_print(x)
if __name__ == "__main__":
n = 20
m = 8
print(solve_ls(cholesky(form_a(n, m)), form_b(n)))
