import numpy as np
from scipy.linalg import lu

---

Gauss, LU, Cholesky

http://homepages.warwick.ac.uk/~ecsgaj/matrixAlgSlidesC.pdf

https://www.quantstart.com/articles/Cholesky-Decomposition-in-Python-and-NumPy/

---

invert UT matrix

• I rewrote this using recursion below
• it's 3 or 4 lines max lines if you do it properly

https://stackoverflow.com/questions/420612/is-there-around-a-straightforward-way-to-invert-a-triangular-upper-or-lower-ma

n = size(A,1);
B = zeros(n);
for i=1:n
    B(i,i) = 1/A(i,i);
    for j=1:i-1
        s = 0;
        for k=j:i-1
            s = s + A(i,k)*B(k,j);
        end
        B(i,j) = -s*B(i,i);
    end
end

---

Code for a Gaussian elimination

Solves a linear system

$A \vec{x} = \vec{y}$

def gaussElim(A,y):
    a = A.copy()
    b = y.copy()
    n = len(b)
    # Elimination phase
    for k in range(0,n-1):
        for i in range(k+1,n):
            if a[i,k] != 0.0:
                #if not null define λ
                lam = a [i,k]/a[k,k]
                #we calculate the new row of the matrix
                a[i,k+1:] = a[i,k+1:] - lam*a[k,k+1:]
                #we update vector b
                b[i] = b[i] - lam*b[k]
                
    # backward substitution
    for k in range(n-1,-1,-1):
        b[k] = (b[k] - a[k,k+1:] @ b[k+1:])/a[k,k]
    
    return a,b

#initial coefficients
a=np.array([[1.0,1.0,1.0],[1.0,-1.0,-1.0],[1.0,-2.0,3.0]])
b=np.array([1.0,1.0,-5.0])

U,x = gaussElim(a,b)
#print A transformed for check
print(a)
print("x =\n",x)
#det = np.prod(np.diagonal(a)) #determinant
#print("\ndet =",det)
#check result and numerical precision
print("\nCheck result: [A]{x} - b =\n",a @ x - b)

[[ 1.  1.  1.]
 [ 1. -1. -1.]
 [ 1. -2.  3.]]
x =
 [ 1.   1.2 -1.2]

Check result: [A]{x} - b =
 [2.22044605e-16 0.00000000e+00 0.00000000e+00]

---

Basic code for LU

def LU(X):
    def inv(U):
        if U.shape[0] == 1: return U
        X = U.copy()
        # upper block is calculated recursively
        X[:-1,:-1] = inv(U[:-1,:-1])
        X[-1,:-1] = -U[-1,:-1] @ X[:-1,:-1]
        return X
    
    #gaussian pivot
    n = len(X)
    U = np.copy(X)
    L = np.identity(n)
    for k in range(0,n-1):
        for i in range(k+1,n):
            if U[i,k] == 0.0:
                raise ValueError
            #should fix this to permute
            #if not null define λ
            lam = U [i,k]/U[k,k]
            #we calculate the new row of the matrix
            U[i] -= lam*U[k]
            #we are calculating the inverse of L
            L[i] -= lam*L[k]
    
    return inv(L), U

a = np.random.randint(1, high = 5, size=16).reshape(4,4).astype(float)

L, U = LU(a)
L @ U - a

array([[0., 0., 0.],
       [0., 0., 0.],
       [0., 0., 0.]])

---

Exo 1

Numerical stability

https://dspace.mit.edu/bitstream/handle/1721.1/75282/18-335j-fall-2006/contents/lecture-notes/lec12.pdf

def M(t):
    return np.array([t,1,1,1]).reshape(2,2)

B = np.array([1,2])

M(1) 

array([[1, 1],
       [1, 1]])

t = .1**20
A = M(t)
A

array([[1.e-20, 1.e+00],
       [1.e+00, 1.e+00]])

A[1] = A[1] - A[0]/A[0,0]
A

array([[ 1.e-20,  1.e+00],
       [ 0.e+00, -1.e+20]])

M(2) @ np.linalg.inv(M(2)) @ B

array([[ 1.e-16,  1.e+00],
       [ 0.e+00, -1.e+16]])

P,L,U = lu(A)

L@U

array([[ 1.e-20,  1.e+00],
       [ 0.e+00, -1.e+20]])

L, U

(array([[1., 0.],
        [0., 1.]]),
 array([[ 1.e-20,  1.e+00],
        [ 0.e+00, -1.e+20]]))

---

Exo 2

A = np.array([int(x) for x in '-2 3 0 1 0 -4 2 0 5'.split()]).reshape(3,3)

P,L,U = lu(A)

L,U

(array([[ 1. ,  0. ,  0. ],
        [-1. ,  1. ,  0. ],
        [-0.5,  0.5,  1. ]]),
 array([[-2. ,  3. ,  0. ],
        [ 0. ,  3. ,  5. ],
        [ 0. ,  0. , -6.5]]))

P@L@U

array([[-2.,  3.,  0.],
       [ 1.,  0., -4.],
       [ 2.,  0.,  5.]])

---

Determinants

the determinant of L is 1 so

det A = det P det U

np.linalg.det(P), np.linalg.det(U), np.linalg.det(A)

(-1.0, 38.99999999999999, -38.99999999999999)

some row operations

A[1] + .5*A[0], A[2] + A[0]

(array([ 0. ,  1.5, -4. ]), array([0, 3, 5]))

array([[-2.,  3.,  0.],
       [ 1.,  0., -4.],
       [ 2.,  0.,  5.]])

---

Sort the rows

...so that the ones with first coeff close to 0 are at the end

B = list(A)
B.sort(key=lambda x: -abs(x[0]))
B

[array([-2,  3,  0]), array([2, 0, 5]), array([ 1,  0, -4])]

---

P is trivial now

lu(B)

(array([[1., 0., 0.],
        [0., 1., 0.],
        [0., 0., 1.]]),
 array([[ 1. ,  0. ,  0. ],
        [-1. ,  1. ,  0. ],
        [-0.5,  0.5,  1. ]]),
 array([[-2. ,  3. ,  0. ],
        [ 0. ,  3. ,  5. ],
        [ 0. ,  0. , -6.5]]))

---

Solution to the system

C = np.zeros((3,4))
C[:,:-1] = B
C[:,-1] = [1,2,3]
C

array([[-2.,  3.,  0.,  1.],
       [ 2.,  0.,  5.,  2.],
       [ 1.,  0., -4.,  3.]])

C = list(C)
C.sort(key=lambda x: -abs(x[0]))
C = np.array(C)
C

array([[-2.,  3.,  0.,  1.],
       [ 2.,  0.,  5.,  2.],
       [ 1.,  0., -4.,  3.]])

L, U = LU(C[:,:-1])

L@U

array([[-2.,  3.,  0.],
       [ 2.,  0.,  5.],
       [ 1.,  0., -4.]])

def Linv_U(X):
    #gaussian pivot
    n = len(X)
    U = np.copy(X)
    L = np.identity(n)
    for k in range(0,n-1):
        for i in range(k+1,n):
            if U[i,k] == 0.0:
                raise ValueError
            #should fix this to permute
            #if not null define λ
            lam = U [i,k]/U[k,k]
            #we calculate the new row of the matrix
            U[i] -= lam*U[k]
            #we are calculating the inverse of L
            L[i] -= lam*L[k]
    
    return L, U

---

Hilbert matrix

N= 5
A = np.zeros((N,N))

for i in range(0,N):
    for j in range(0,N):
        A[i,j] = 1/(i+j+1)

A

array([[1.        , 0.5       , 0.33333333, 0.25      , 0.2       ],
       [0.5       , 0.33333333, 0.25      , 0.2       , 0.16666667],
       [0.33333333, 0.25      , 0.2       , 0.16666667, 0.14285714],
       [0.25      , 0.2       , 0.16666667, 0.14285714, 0.125     ],
       [0.2       , 0.16666667, 0.14285714, 0.125     , 0.11111111]])

---

X = np.random.randint(8,size= 9).reshape(3,3)
X = np.triu(X, k=0)
for i in range(X.shape[0]):
    X[i,i]= 1
    
X

array([[1, 2, 5],
       [0, 1, 4],
       [0, 0, 1]])

---

invert an upper triangular unipotent matrix by recursion

def ran_uni(N):
    X = np.random.randint(8,size= N*N).reshape(N,N)
    X = np.triu(X, k=0)
    for i in range(N):
        X[i,i]= 1
    return X


def inv_uni(U): 
    if U.shape[0] == 1: return U
    #should check its unimodular
    X = U.copy()
    # lower block is calculated recursively
    X[1:,1:] = inv_uni(U[1:,1:])
    # first line is a matrix multiplication
    X[0,1:] = -U[0,1:] @ X[1:,1:]
    
    return X

V = ran_uni(5)

V @ inv_uni(V)

array([[1, 0, 0, 0, 0],
       [0, 1, 0, 0, 0],
       [0, 0, 1, 0, 0],
       [0, 0, 0, 1, 0],
       [0, 0, 0, 0, 1]])

inv_uni(V)

array([[  1,  -7,  19, -98, -30],
       [  0,   1,  -3,  15,   5],
       [  0,   0,   1,  -5,  -2],
       [  0,   0,   0,   1,  -1],
       [  0,   0,   0,   0,   1]])

V

array([[1, 6, 4, 7, 5],
       [0, 1, 1, 7, 7],
       [0, 0, 1, 0, 1],
       [0, 0, 0, 1, 2],
       [0, 0, 0, 0, 1]])