This notebook presents an experimental implementation of the LU decomposition. It is not meant for use in the wild. Its purpose is to provide a tool to actually see LU decomposition as it arises during Gaussian elimination. It uses lists of lists, rather than NumPy arrays, so that you can experiment with objects of various precision - like Decimal.
lu_step¶lu_step is the main function allowing you to see the action. It accepts four mandatory inputs:
LLPartial: The multipliers below the diagonal of ones in current L matrixUpartial: The current U matrix obtained by Gaussian eliminationrow_col: The current step or position along the diagonal that we're woring onpivot: The pivot position you want to use.In addition, there are two inputs with default values:
show: Set to True to print the current resultswidth: A parameter to tweak the pretty printing.Here's the function:
import pprint as pp
from copy import deepcopy
def lu_step(Lpartial, Upartial, row_col, pivot, show = True, width = 50):
L = deepcopy(Lpartial) # Need to deepcopy these to avoid the contents
U = deepcopy(Upartial) # of the input from being overwritten.
dim = len(U)
for i in range(dim):
L[i][i] = 0
multipliers = []
row,col = row_col,row_col
perm = list(range(dim))
temp = perm[pivot]
perm[pivot] = perm[row]
perm[row] = temp
U = [U[i] for i in perm]
L = [L[i] for i in perm]
pivot_term = U[row][col]
for i in range(row+1,dim):
multiplier = -U[i][col]/pivot_term
multipliers.append(multiplier)
for j in range(dim):
U[i][j] = U[i][j]+multiplier*U[row][j]
for i in range(row+1,dim):
L[i][col] = -multipliers[i-(row+1)]
for i in range(dim):
L[i][i] = 1
if show == True:
pp.pprint('-------')
pp.pprint('Step ' + str(row_col))
pp.pprint('-------')
pp.pprint('L = ')
pp.pprint(L, width = width)
pp.pprint('-------')
pp.pprint('U = ')
pp.pprint(U, width = width)
return L,U
And here's how to use it. Note that we are using NumPy simply for the convenience in some spots, like the setup of an identity matrix. It is not used in the computations.
import numpy as np
A = [
[10,-7,0],
[-3,2,6],
[5,-1,5]
]
Linit = np.identity((len(A))).tolist()
L,U = lu_step(Linit, A, 0, 0, width = 40)
Note that the resulting L matrix has the additive inverses of the multipliers used in Gaussian elimination below the ones on the diagonal and the resulting U matrix contains the result of the Gaussian elimination. Next, we can use the L and U computed to this point to take the next step. We also need to indicate that we've moved from row_col zero to row_col one and we'll need to indicate where we want the pivot to be. The pivot should be at least as large as row_col but not larger than the dimension of the matrix. Looking at the current state of U, the best choice of a pivot should be 2 - right? Let's see what happens if we use naive Gaussian elimination and use 1 as the pivot.
L,U = lu_step(L, U, 1, 1, width = 40)
np.array(L).dot(np.array(U))
It appears to have worked!
Here's an example that (I think) shows what goes wrong when the wrong pivot is chosen.
from decimal import Decimal, getcontext
context = getcontext()
context.prec = 5
n = 8
A = [[10**(-n), 1], [1,1]]
dim = len(A)
for i in range(dim):
for j in range(dim):
A[i][j] = Decimal(A[i][j])
zero = [[Decimal(0) for i in range(dim)] for j in range(dim)]
L,U = lu_step([[1,0],[0,1]],A,0,0)
L = np.array(L)
U = np.array(U)
L.dot(U)
The function lu_decomp allows you to visualize the whole decompostion, rather than just a single step. Of course, the choice of the pivot should be automated.
def choose_pivot(M,k):
row,col = k,k
dim = len(M)
indexed_elements = [(abs(M[i][col]),i) for i in range(row,dim)]
pivot = max(indexed_elements)[1]
return pivot
def lu_decomp(M, show = True, width = 40):
U = deepcopy(M)
dim = len(U)
U = [[float(U[i][j]) for j in range(dim)] for i in range(dim)]
L = np.identity(dim).tolist()
perm = list(range(dim))
for k in range(dim):
pivot = choose_pivot(U,k)
temp = perm[pivot]
perm[pivot] = perm[k]
perm[k] = temp
L,U = lu_step(L, U, k, pivot, show=show, width=width)
if show == True:
pp.pprint('perm = ' + str(perm))
return L,U,np.identity(dim)[perm].transpose().tolist()
from numpy.random import seed, randint
seed(1)
A = randint(-8,8,(4,4))
A
L,U,P = lu_decomp(A, width = 60)
L = np.array(L)
U = np.array(U)
P = np.array(P)
P.dot(L.dot(U))