import numpy as np
import numpy.linalg as la
n = 3
np.random.seed(15)
A = np.round(5*np.random.randn(n, n))
A
Initialize L
and U
with zeros:
L = np.zeros((n,n))
U = np.zeros((n,n))
Set U
to be the first row of A
:
U[0,:] = A[0,:]
U
Compute the first column of L
:
L[:,0] = A[:,0]/U[0,0]
L
Compare what we have to A
:
print(A)
print(L@U)
Perform the Schur complement update and store the result in A1
:
A1 = A - L @ U
A1
Take the second row of U
to be the second row of A1
:
U[1,1:] = A1[1,1:]
U
We can now compute the next column of L
:
L[1:,1] = A1[1:,1]/U[1,1]
L
And finally, compute the bottom right elements of L
and U
U[2,2] = A1[2,2] - L[2,1]*U[1,2]
L[2,2] = 1.0
print(L)
print(U)
print(A)
print(L@U)
Implement the general LU factorization algorithm
n = 4
A = np.random.random((n,n))
L = np.zeros((n,n))
U = np.zeros((n,n))
M = A.copy()
for i in range(n):
U[i,i:] = M[i,i:]
L[i:,i] = M[i:,i]/U[i,i]
M[i+1:,i+1:] -= np.outer(L[i+1:,i:i+1],U[i:i+1,i+1:])
print(L)
print(U)
print(A-L@U)
When a divisor (U_ii) becomes zero, while nonzeros remain in the trailing matrix, the LU factorization does not exist. Further, if U_ii is small, the magnitude of elements in L and U can blow up, increasing round-off error. Partial pivoting circumvents both problems.
n = 4
A = np.random.random((n,n))
L = np.zeros((n,n))
U = np.zeros((n,n))
P = np.eye(n,n)
M = A.copy()
def perm(i,j):
P = np.eye(n,n)
P[i,i] = 0
P[j,j] = 0
P[i,j] = 1
P[j,i] = 1
return P
for i in range(n):
imax = np.argmax(np.abs(M[i:,i]))
Pi = perm(i,i+imax)
P = Pi @ P
L = Pi @ L
M = Pi @ M
U[i,i:] = M[i,i:]
L[i:,i] = M[i:,i]/U[i,i]
M[i+1:,i+1:] -= np.outer(L[i+1:,i:i+1],U[i:i+1,i+1:])
Check residual error for the permuted A, i.e., PA-LU. What property do off-diagonal entries in L satisfy?
print(P@A - L @ U)
L
P