In [39]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In our last Jupyter Notebook we learned how to solve 1D equations. Now, we'll go to the next level and will learn how to solve not just one equation, but a system of linear equations. This is a set of $n$ equations involving $n$ variables wherein all the equations must be satisfied at the same time. You probably know how to solve small 2D systems with methods such as substitution and reduction, but in practical real-life situations it's very likely that you'll find problems of bigger dimensions. As usual, we'll present some useful methods for solving systems of linear equations below.
Firstly, we will study direct methods. They compute the analytic solution of the system (from here comes the name direct) limited only by the loss of numerical precision, because of the arithmetic operations performed by the computer. Their counterpart is the iterative methods, which calculate an approximate solution that evolves iteratively converging to the real solution.
Given the matrix $A \in \mathbb{R}^{n \times n}$ square and non singular, the main goal of this method involves finding a decomposition like $A = L U$ where $L,U \in \mathbb{R}^{n \times n}$ are lower and upper triangular matrices respectively.
The algorithm to perform this decomposition is basically a modified version of Gaussian Elimination. It basically iterates through the first $n-1$ columns, making $0$ all the entries below the main diagonal. This is accomplished by performing row operations.
In [40]:
def lu_decomp(A, show=False):
N,_ = A.shape
U = np.copy(A)
L = np.identity(N)
if show:
print('Initial matrices')
print('L = '); print(np.array_str(L, precision=2, suppress_small=True))
print('U = '); print(np.array_str(U, precision=2, suppress_small=True))
print('----------------------------------------')
#iterating through columns
for j in range(N-1):
#iterating through rows
for i in range(j+1,N):
L[i,j] = U[i,j]/U[j,j]
U[i] -= L[i,j]*U[j]
if show:
print('L = '); print(np.array_str(L, precision=2, suppress_small=True))
print('U = '); print(np.array_str(U, precision=2, suppress_small=True))
print('----------------------------------------')
return L,U
Once the decomposition is done, solving a linear system like $A x = b$ is straightforward:
$$A x = b \rightarrow L U x = b \ \ \text{ if we set } \ \ U x = c \rightarrow L c = b \ \ \text{ (solve for c) } \ \rightarrow U x = c$$and as you might know, solving lower and upper triangular systems can be easily performed by back-substitution and forward-subsitution respectively.
In [41]:
"""
Solves a linear system A x = b, where A is a
triangular (upper or lower) matrix
"""
def solve_triangular(A, b, upper=True):
n = b.shape[0]
x = np.zeros_like(b)
if upper==True:
#perform back-substitution
x[-1] = (1./A[-1,-1]) * b[-1]
for i in range(n-2, -1, -1):
x[i] = (1./A[i,i]) * (b[i] - np.sum(A[i,i+1:] * x[i+1:]))
else:
#perform forward-substitution
x[0] = (1./A[0,0]) * b[0]
for i in range(1,n):
x[i] = (1./A[i,i]) * (b[i] - np.sum(A[i,:i] * x[:i]))
return x
def solve_lu(A, b, show=False):
L,U = lu_decomp(A, show)
# L.c = b with c = U.x
c = solve_triangular(L, b, upper=False)
x = solve_triangular(U, c)
return x
Let's now try our implementations. We begin by creating a random 100$\times$100 linear system:
In [42]:
A = np.random.random((3,3))
b = np.ones(3)
and then we compute the solution with our LU solver, and aditionally with the NumPy solver which computes the solution using LAPACK routines.
In [43]:
lu_sol = solve_lu(A,b, show=True)
np_sol = np.linalg.solve(A,b)
in order to compare these huge vectors, we use the Euclidean metric as follows:
In [44]:
np.linalg.norm(lu_sol - np_sol)
Out[44]:
which is a very good result!
This method has two important facts to be noted:
L[i,j] = U[i,j]/U[j,j]. Here we divide an entry below the main diagonal by the pivot value. What happens if the pivot equals 0? How can we prevent that? Answer: PALU.As you might've noted previously, LU has a problem when a pivot has the value of $0$. To handle this problem, we add row permutations to the original LU algorithm. The procedure is as follows:
To keep track of all the permutations performed, we use the permutation matrix $P$. It's inicially an identity matrix which permutes its rows in the same way the algorithm does on the resulting matrix.
In [45]:
#permutation between rows i and j on matrix A
def row_perm(A, i, j):
tmp = np.copy(A[i])
A[i] = A[j]
A[j] = tmp
def palu_decomp(A, show=False):
N,_ = A.shape
P = np.identity(N)
L = np.zeros((N,N))
U = np.copy(A)
if show:
print('Initial matrices')
print('P = '); print(np.array_str(P, precision=2, suppress_small=True))
print('L = '); print(np.array_str(L, precision=2, suppress_small=True))
print('U = '); print(np.array_str(U, precision=2, suppress_small=True))
print('----------------------------------------')
#iterating through columns
for j in range(N-1):
#determine the new pivot
p_index = np.argmax(np.abs(U[j:,j]))
if p_index != 0:
row_perm(P, j, j+p_index)
row_perm(U, j, j+p_index)
row_perm(L, j, j+p_index)
if show:
print('A permutation has been made')
print('P = '); print(np.array_str(P, precision=2, suppress_small=True))
print('L = '); print(np.array_str(L, precision=2, suppress_small=True))
print('U = '); print(np.array_str(U, precision=2, suppress_small=True))
print('----------------------------------------')
#iterating through rows
for i in range(j+1,N):
L[i,j] = U[i,j]/U[j,j]
U[i] -= L[i,j]*U[j]
if show:
print('P = '); print(np.array_str(P, precision=2, suppress_small=True))
print('L = '); print(np.array_str(L, precision=2, suppress_small=True))
print('U = '); print(np.array_str(U, precision=2, suppress_small=True))
print('----------------------------------------')
np.fill_diagonal(L,1)
return P,L,U
The procedure to solve the system $Ax=b$ remains almost the same. We have to add the efect of the permutation matrix $P$:
$$A x = b \rightarrow P A x = P b \rightarrow L U x = b' \ \ \text{ if we set } \ \ U x = c \rightarrow L c = b' \ \ \text{ (solve for c) } \ \rightarrow U x = c$$
In [46]:
def solve_palu(A, b, show=False):
P,L,U = palu_decomp(A, show)
#A.x = b -> P.A.x = P.b = b'
b = np.dot(P,b)
# L.c = b' with c = U.x
c = solve_triangular(L, b, upper=False)
x = solve_triangular(U, c)
return x
Let's test this new method against the LU and NumPy solvers
In [47]:
palu_sol = solve_palu(A, b, show=True)
In [48]:
np.linalg.norm(palu_sol - lu_sol)
Out[48]:
In [49]:
np.linalg.norm(palu_sol - np_sol)
Out[49]:
Here are some questions about PALU:
This is another direct method only applicable to symmetric positive-definite matrices. In order to try this algorithm we have to create this kind of matrices. The next function generates random symmetric positive-definite matrices.
In [50]:
"""
Randomly generates an nxn symmetric positive-
definite matrix A.
"""
def generate_spd_matrix(n, flag=True):
if flag:
A = np.random.random((n,n))
# Constructing symmetry
A += A.T
# A = np.dot(A.T,A) # Another way
#symmetric+diagonally dominant -> symmetric positive-definite
deltas = 0.1*np.random.random(n)
row_sum = A.sum(axis=1)-np.diag(A)
np.fill_diagonal(A, row_sum+deltas)
else:
B = np.random.random((n,n))
# A way to make sure the quadratic form is greater or equal to zero:
# this means x^T*B^T\B*x >= ||B*x||, but if B is singular, it could be zero.
A = np.dot(B.T,B)
# To avoid a being singular, we just add a positive diagonal matrix
A = A + np.eye(n)
return A
Given a symmetric positive-definite matrix $A \in \mathbb{R}^{n \times n}$, the Cholesky decomposition is of the form $A =R^T R$, with $R$ being an upper triangular matrix. This method takes advantage of the properties of symmetric matrices, reaching approximately twice the efficiency of LU.
In [51]:
def cholesky_decomp(A, show=False):
N,_ = A.shape
A = np.copy(A)
R = np.zeros((N,N))
if show:
print('Initial matrix')
print('A = '); print(np.array_str(A, precision=2, suppress_small=True))
print('R = '); print(np.array_str(R, precision=2, suppress_small=True))
print('----------------------------------------')
for i in range(N):
R[i,i] = np.sqrt(A[i,i])
u = (1./R[i,i])*A[i,i+1:]
R[i,i+1:] = u
A[i+1:,i+1:] -= np.outer(u,u)
if show:
print('A = '); print(np.array_str(A, precision=2, suppress_small=True))
print('R = '); print(np.array_str(R, precision=2, suppress_small=True))
print('----------------------------------------')
return R
The solve stage remains the same as LU:
In [52]:
def solve_cholesky(A, b, show=False):
R = cholesky_decomp(A, show)
#R^T.R.x = b -> R^T.c = b with R.x = c
c = solve_triangular(R.T, b, upper=False)
x = solve_triangular(R, c)
return x
Now we test our implementation, comparing time execution with LU and PALU on two different linear systems
In [53]:
A = generate_spd_matrix(3)
b = np.ones(3)
In [54]:
b=np.array([4,2,0])
solve_cholesky(A, b, show=True)
Out[54]:
In [55]:
A = generate_spd_matrix(100)
b = np.ones(100)
In [56]:
%timeit solve_cholesky(A, b)
%timeit solve_lu(A, b)
%timeit solve_palu(A, b)
In [57]:
A = generate_spd_matrix(1000)
b = np.ones(1000)
In [58]:
%timeit solve_cholesky(A, b)
%timeit solve_lu(A, b)
%timeit solve_palu(A, b)
In [59]:
"""
Randomly generates an nxn strictly diagonally
dominant matrix A.
"""
def generate_dd_matrix(n):
A = np.random.random((n,n))
deltas = 0.1*np.random.random(n)
row_sum = A.sum(axis=1)-np.diag(A)
np.fill_diagonal(A, row_sum+deltas)
return A
"""
Computes relative error between each row on
X matrix and y vector.
"""
def error(X, y):
D = X-y
err = np.linalg.norm(D, axis=1, ord=np.inf)
return err
As before we will create a linear system $A x = b$, with $A$ as a diagonally dominant matrix, which is a sufficient condition for the methods we will study in this section converge
In [60]:
A = np.array([[3, -1, 0, 0, 0, 0.5],[-1, 3, -1, 0, 0.5, 0],[0, -1, 3, -1, 0, 0],[0, 0, -1, 3, -1, 0],
[0, 0.5, 0, -1, 3, -1],[0.5, 0, 0, 0, -1, 3]])
b = np.array([2.5, 1.5, 1., 1., 1.5, 2.5])
print ('A='); print (A)
print ('b='); print (b)
and find the solution $x$ through np.linalg.solve to use it as the reference solution-
In [61]:
np_sol = np.linalg.solve(A,b)
In [62]:
"""
Iterative methods implementations returns an array X
with the the solutions at each iteration
"""
def jacobi(A, b, n_iter=50):
n = A.shape[0]
#array with solutions
X = np.empty((n_iter, n))
#initial guess
X[0] = np.zeros(n)
#submatrices
D = np.diag(A)
Dinv = D**-1
R = A - np.diag(D) # R = (L+U)
for i in range(1, n_iter):
# X[i] = Dinv*(b - np.dot(R, X[i-1]))
# v1.12
ri = b - np.dot(A, X[i-1])
X[i] = X[i-1]+Dinv*ri
return X
def jacobi_M(A):
L = np.tril(A,-1)
U = np.triu(A,1)
D = np.diag(np.diag(A))
M = -np.dot(np.linalg.inv(D),L+U)
return M
# \mathbf{x}_{n+1}=M\,\mathbf{x}_{n}+\mathbf{c}
Now let's resolve the same linear system with Jacobi method!
In [63]:
jac_sol = jacobi(A,b, n_iter=50)
In [64]:
jac_err = error(jac_sol, np_sol)
it = np.linspace(1, 50, 50)
In [65]:
plt.figure(figsize=(12,6))
plt.semilogy(it, jac_err, marker='o', linestyle='--', color='b')
plt.grid(True)
plt.xlabel('Iterations')
plt.ylabel('Error')
plt.title('Infinity norm error for Jacobi method')
plt.show()
In [66]:
Mj = jacobi_M(A)
print(np.linalg.norm(Mj))
In [67]:
np.linalg.eigvals(Mj)
Out[67]:
In [68]:
np.abs(np.linalg.eigvals(Mj))
Out[68]:
In [69]:
np.max(np.abs(np.linalg.eigvals(Mj)))
Out[69]:
In [70]:
def gauss_seidel(A, b, n_iter=50):
n = A.shape[0]
#array with solutions
X = np.empty((n_iter, n))
#initial guess
X[0] = np.zeros(n)
#submatrices
R = np.tril(A) #R=(L+D)
U = A-R
for i in range(1, n_iter):
#X[i] = solve_triangular(R, b-np.dot(U, X[i-1]), upper=False)
# v1.11
X[i] = X[i-1]+solve_triangular(R, b-np.dot(A, X[i-1]), upper=False)
return X
def gauss_seidel_M(A):
L = np.tril(A,-1)
U = np.triu(A,1)
D = np.diag(np.diag(A))
M = -np.dot(np.linalg.inv(L+D),U)
return M
Now let's resolve the same linear system with Gauss-Seidel method!
In [71]:
gauss_sol = gauss_seidel(A,b)
In [72]:
gauss_err = error(gauss_sol, np_sol)
In [73]:
plt.figure(figsize=(12,6))
plt.semilogy(it, gauss_err, marker='o', linestyle='--', color='r')
plt.grid(True)
plt.xlabel('Iterations')
plt.ylabel('Error')
plt.title('Infinity norm error for Gauss method')
plt.show()
Here are some questions about Gauss-Seidel:
solve_triangular instead of np.linalg.solve or something similar?
In [74]:
Mgs = gauss_seidel_M(A)
print(np.linalg.norm(Mgs))
In [75]:
np.max(np.abs(np.linalg.eigvals(Mgs)))
Out[75]:
In [76]:
def sor(A, b, w=1.05, n_iter=50):
n = A.shape[0]
#array with solutions
X = np.empty((n_iter, n))
#initial guess
X[0] = np.zeros(n)
#submatrices
R = np.tril(A) #R=(L+D)
U = A-R
# v1.11
L = np.tril(A,-1)
D = np.diag(np.diag(A))
M = L+D/w
for i in range(1, n_iter):
#X_i = solve_triangular(R, b-np.dot(U, X[i-1]), upper=False)
#X[i] = w*X_i + (1-w)*X[i-1]
# v1.11
X[i] = X[i-1]+solve_triangular(M, b-np.dot(A, X[i-1]), upper=False)
return X
def sor_M(A,w=1.05):
L = np.tril(A,-1)
U = np.triu(A,1)
D = np.diag(np.diag(A))
M = np.dot(np.linalg.inv(w*L + D),((1-w)*D -w*U))
return M
Now let's resolve the same linear system with Jacobi method!
In [77]:
sor_sol = sor(A, b, w=1.15)
In [78]:
sor_err = error(sor_sol, np_sol)
In [79]:
plt.figure(figsize=(12,6))
plt.semilogy(it, sor_err, marker='o', linestyle='--', color='g')
plt.grid(True)
plt.xlabel('Iterations')
plt.ylabel('Error')
plt.title('Infinity norm error for SOR method')
plt.show()
In [80]:
Msor = sor_M(A)
print(np.linalg.norm(Msor))
In [81]:
np.max(np.abs(np.linalg.eigvals(Msor)))
Out[81]:
How can we choose a good value of $\omega$? Well there are some methods you could search, but for now we will try a naive way, i.e, computing the solution for a range $\omega \in [1,1.3]$ as follows:
In [82]:
n = 30 #width of subdivisions
sor_solutions = list()
for w in np.linspace(1., 1.3, n):
sor_solutions.append(sor(A, b, w, n_iter=5)[-1])
np.asarray(sor_solutions)
#now compute error solutions with each w
sor_errors = error(sor_solutions, np_sol)
w = np.linspace(1., 1.3, n)
as you can see, we compute the SOR solution with 5 iterations for each $\omega$ on the given range.
In [83]:
plt.figure(figsize=(12,6))
plt.semilogy(w, sor_errors, marker='o', linestyle='--', color='g')
plt.grid(True)
plt.xlabel('w')
plt.ylabel('Error')
plt.title('Infinity norm error after 5 steps of SOR as a function of w')
plt.show()
Here are some questions about SOR:
Let's see convergence plots all together
In [84]:
plt.figure(figsize=(12,6))
plt.semilogy(it, jac_err, marker='o', linestyle='--', color='b', label='Jacobi')
plt.semilogy(it, gauss_err, marker='o', linestyle='--', color='r', label='Gauss-Seidel')
plt.semilogy(it, sor_err, marker='o', linestyle='--', color='g', label='SOR')
plt.grid(True)
plt.xlabel('Iterations')
plt.ylabel('Error')
plt.title('Infinity norm error for all methods')
plt.legend(loc=0)
plt.show()
Now that you know how to solve systems of linear equations problem with these methods, let's try to answer a few questions!
$a)$ Find the values of $\alpha$ that make possible to do a LU descomposition of the following matrix:
$$ \begin{bmatrix} \alpha & 2 \\[0.3em] 1 & \alpha \end{bmatrix} $$$b)$- Let $A$ be the following matrix:
$$ A = \begin{bmatrix} 2 & 4 & 2 \\[0.3em] -1 & 1 & 2 \\[0.3em] -1 & -3 & -1 \end{bmatrix} $$Find the PALU descomposition of the matrix $A$.
Solve the system of equations $Ax = \[1 , \frac{1}{2}, \frac{1}{3}\]^T$.
$c)$ Considering this matrix:
$$ \begin{bmatrix} 1 & 1 & 0 \\[0.3em] 1 & 5 & 2 \\[0.3em] 0 & 2 & 3 \end{bmatrix} $$Find the LU descomposition.
Find the Cholesky descomposition.
Compare the efficiency of both methods.
$d)$ Use Jacobi, Gauss Seidel, and SOR to solve the following system of equations (number of iterations = 2):
$$2x + y = 3$$$$x + 2y = 2$$$e)$ Explain the pros and cons of using iterative methods instead of the direct ones.
In [85]:
from scipy.linalg import hilbert
In [86]:
N=20
errors=np.zeros(N+1)
kappas=np.zeros(N+1)
my_range=np.arange(5,N+1)
for n in my_range:
A=hilbert(n)
x_exact=np.ones(n)
b=np.dot(A,x_exact)
x=np.linalg.solve(A,b)
errors[n]=np.linalg.norm(x-x_exact,ord=np.inf)
kappas[n]=np.linalg.cond(A,2)
In [87]:
plt.figure(figsize=(12,6))
plt.semilogy(my_range, 10.**(-16+np.log10(kappas[my_range])), marker='o', linestyle='--', color='b',label='Estimated forward error')
plt.semilogy(my_range, errors[my_range], marker='o', linestyle='--', color='r',label='Forward error')
plt.grid(True)
plt.xlabel('n')
#plt.ylabel('$\kappa$ and errors')
plt.title('')
plt.legend(loc=0)
plt.show()
(1/Kappa(A))||b-A x_a||/||b|| <= ||x-x_a||/||x||<= Kappa(A)||b-A x_a||/||b||
In [88]:
n=200
A=hilbert(n)
x_exact=np.ones(n)
b=np.dot(A,x_exact)
x=np.linalg.solve(A,b)
kappa=np.linalg.cond(A,2)
print(np.log10(kappa))
In [89]:
np.linalg.norm(b-np.dot(A,x))
Out[89]:
In [90]:
kappa
Out[90]:
In [91]:
np.linalg.norm(x-x_exact)
Out[91]:
In [92]:
x_exact[0]
Out[92]:
In [93]:
x
Out[93]:
In [94]:
plt.plot(x,'.')
plt.show()
In [95]:
lu_decomp(np.array([[1e-20, 1],[1,2]]), show=True)
Out[95]:
In [ ]: