import numpy as np
import matplotlib.pyplot as plt
import scipy
from scipy import linalg
It turns out that a proper permutation in rows (or columns) is sufficient for LU factorization. LU factorization with partial pivoting (LUP) refers often to LU factorization with row permutations only:
$PA=LU$
where
It turns out that all square matrices can be factorized in this form Pavel Okunev, Charles R. Johnson 1997
and the factorization is numerically stable in practice.
This makes LUP decomposition a useful technique in practice.
Look at discussion and python 2 code
A = np.array([1,4,-1,2]).reshape(2,2)
scipy.linalg.lu(A)
(array([[1., 0.],
[0., 1.]]),
array([[1., 0.],
[1., 1.]]),
array([[1., 1.],
[0., 1.]]))
B = [-1,3,0,1,0,-4,2,0,5]
B = np.array(B).reshape(3,3)
B
array([[-1, 3, 0],
[ 1, 0, -4],
[ 2, 0, 5]])
P,L,U = scipy.linalg.lu(B)
P @ L @ U
array([[-1., 3., 0.],
[ 1., 0., -4.],
[ 2., 0., 5.]])
P
array([[0., 1., 0.],
[0., 0., 1.],
[1., 0., 0.]])
np.linalg.inv(P) @ B
array([[ 2., 0., 5.],
[-1., 3., 0.],
[ 1., 0., -4.]])
np.linalg.solve(B,[1,2,3])
array([ 1.69230769, 0.8974359 , -0.07692308])
np.linalg.inv(P) @ [1,2,3]
array([3., 1., 2.])
Given a system of linear equations in matrix form
$ A\mathbf {x} =\mathbf {b} $ we want to solve the equation for x, given A and b.
Suppose we have already obtained the $LUP$ decomposition of $A$ such that
$PA=LU$ , so $LU\mathbf {x} =P\mathbf {b}$.
In this case the solution is done in two logical steps:
In both cases we are dealing with triangular matrices (L and U), which can be solved directly by forward and backward substitution without using the Gaussian elimination process (however we do need this process or equivalent to compute the LU decomposition itself).
a = 1
A = np.array([a,1,1,2]).reshape(2,2)
import scipy
scipy.linalg.lu(A)
(array([[1., 0.],
[0., 1.]]),
array([[1., 0.],
[1., 1.]]),
array([[1., 1.],
[0., 1.]]))
a = 0
A = np.array([a,1,1,2]).reshape(2,2)
P, L , U = scipy.linalg.lu(A)
P
array([[0., 1.],
[1., 0.]])
P @ U, U
(array([[0., 1.],
[1., 2.]]),
array([[1., 2.],
[0., 1.]]))
a = .5
A = np.array([a,1,1,2]).reshape(2,2)
P, L , U = scipy.linalg.lu(A)
L, U
(array([[1. , 0. ],
[0.5, 1. ]]),
array([[1., 2.],
[0., 0.]]))
In linear algebra, the Cholesky decomposition is a decomposition of a matrix $\mathbf {A}$ which is
into the product of a lower triangular matrix and its conjugate transpose ie
$\mathbf {A} =\mathbf {LL} ^{*}$
The recipe for the coefficients of $\mathbf {L}$ is :
A = np.array([1,2,1,2,13, -1,1,-1, 3]).reshape(3,3)
A
array([[ 1, 2, 1],
[ 2, 13, -1],
[ 1, -1, 3]])
import scipy.linalg
scipy.linalg.cholesky(A).T
array([[ 1., 0., 0.],
[ 2., 3., 0.],
[ 1., -1., 1.]])
np.array(cholesky(A))
array([[ 1., 0., 0.],
[ 2., 3., 0.],
[ 1., -1., 1.]])
def cholesky(A):
'''Performs a Cholesky decomposition
A, which must be a symmetric and positive definite matrix.
returns L = lower variant triangular matrix
such that A = L L^*'''
n = len(A)
# Initialise L as the zero matrix
L = np.zeros((n,n))
for i in range(n):
# under the diagonal
for k in range(i):
# LaTeX: l_{ik} = \frac{1}{l_{kk}} \left( a_{ik} - \sum^{k-1}_{j=1} l_{ij} l_{kj} \right)
L[i,k] = ( (A[i,k] - L[i,:] @ L[k,:]) / L[k,k] )
#on the diagonal
# LaTeX: l_{kk} = \sqrt{ a_{kk} - \sum^{k-1}_{j=1} l^2_{kj}}
L[i,i] = sqrt(A[i,i] - L[i,:] @ L[i,:])
return L
L = cholesky(A)
L, L @ L.T
(array([[ 1., 0., 0.],
[ 2., 3., 0.],
[ 1., -1., 1.]]),
array([[ 1., 2., 1.],
[ 2., 13., -1.],
[ 1., -1., 3.]]))
There are a several different algorithms for calculating the matrices $Q$ and $R$.
We will outline the method of Householder Reflections, which is known to be more numerically stable than the alternative Gramm-Schmidt method.
For more details see QR Decomposition using Householder Reflections.
A Householder Reflection is a linear transformation that enables a vector to be reflected through a plane or hyperplane. Essentially, we use this method because we want to create an upper triangular matrix, $R$. The householder reflection is able to carry out this vector reflection such that all but one of the coordinates disappears. The matrix $Q$ will be built up as a sequence of matrix multiplications that eliminate each coordinate in turn, up to the rank of the matrix $A$.
First step create the vector $\mathbb{x}$, which is the $k$-th column of the matrix $A$, for step $k$.
We define $\alpha = -sgn(\mathbb{x}_k)(||\mathbb{x}||)$. The norm $||\cdot||$ used here is the Euclidean norm.
Given the first column vector of the identity matrix, $I$ of equal size to $A$, $\mathbb{e}_1 = (1,0,...,0)^T$, we create the vector $\mathbb{u}$:
\begin{eqnarray*} \mathbb{u} = \mathbb{x} + \alpha \mathbb{e}_1 \end{eqnarray*}Once we have the vector $\mathbb{u}$, we need to convert it to a unit vector, which we denote as $\mathbb{v}$:
\begin{eqnarray*} \mathbb{v} = \mathbb{u}/||\mathbb{u}|| \end{eqnarray*}Second step form the matrix of the Householder transformation $Q$ associated to $\mathbb{v}$ :
\begin{eqnarray*} Q = I - 2 \mathbb{v} \mathbb{v}^T \end{eqnarray*}Third step $Q$ is now an $m\times m$ Householder matrix, with $Q\mathbb{x} = \left( \alpha, 0, ..., 0\right)^T$. We will use $Q$ to transform $A$ to upper triangular form, giving us the matrix $R$.
Write $Q_k$ for $Q$ at the $k$th tep and, since $k=0$ in this first step, we have $Q_0$ as our first Householder matrix. Muliplying by $A$ gives us:
\begin{eqnarray*} Q_0A = \begin{bmatrix} \alpha_1&\star&\dots&\star\\ 0 & & & \\ \vdots & & A' & \\ 0 & & & \end{bmatrix} \end{eqnarray*}The whole process is recursive and we repeat the 3 steps above for the minor matrix $A'$, which will give a second Householder matrix $Q'_1$. Now we have to "pad out" this minor matrix with elements from the identity matrix such that we can consistently multiply the Householder matrices together. Hence, we define $Q_k$ as the block matrix:
\begin{eqnarray*} Q_k = \begin{pmatrix} I_{k-1} & 0\\ 0 & Q_k'\end{pmatrix} \end{eqnarray*}Once we have carried out $t$ iterations of this process we have $R$ as an upper triangular matrix:
\begin{eqnarray*} R = Q_t ... Q_2 Q_1 A \end{eqnarray*}$Q$ is then fully defined as the multiplication of the transposes of each $Q_k$:
\begin{eqnarray*} Q = Q^T_1 Q^T_2 ... Q^T_t \end{eqnarray*}This gives $A=QR$, the QR Decomposition of $A$.
def QR_householder(A):
'''Performs a Householder Reflections based QR becomposition of the
matrix A an np.array
Returns
- Q, an orthogonal matrix
- R upper triangular matrix
such that A = QR.
'''
n = A.shape[0]
# Set R equal to A, and Q to the identity matrix of the same size
R = A
Q = np.identity(n)
# The Householder procedure
for k in range(n-1): # reduce the index by 1 to skip the 1x1 matrix
# get the vectors x, e and the scalar alpha
x = R[k:,k]
e_0 = np.identity(n-k)[0]
alpha = -np.sign(x[0]) * np.linalg.norm(x)
u = x + alpha*e_0
v = u/np.linalg.norm(u)
# matrix of the reflection x -> x - 2<v,x>v
Householder_matrix = np.identity(n-k) - 2*np.array([ v[i]*v for i in range(n-k)])
# pad out the matrix to match A.shape()
Q_k = np.identity(n)
Q_k[k:,k:] = Householder_matrix
Q = Q_k @ Q
R = Q_k @ R
# Q is the inverse of the product of the Q_k
# we need to take the transpose/inverse now
return Q.T, R
A = np.array([[12, -51, 4],
[6, 167, -68],
[-4, 24, -41]])
Q, R = QR_householder(A)
Q, R
(array([[ 0.85714286, 0.39428571, -0.33142857],
[ 0.42857143, -0.90285714, 0.03428571],
[-0.28571429, -0.17142857, -0.94285714]]),
array([[ 1.40000000e+01, 2.10000000e+01, -1.40000000e+01],
[-5.50670620e-16, -1.75000000e+02, 7.00000000e+01],
[ 3.01980663e-16, -3.63797881e-14, 3.50000000e+01]]))
def QR(A):
'''Performs a Householder Reflections based QR becomposition of the
matrix A an np.array
Returns
- Q, an orthogonal matrix
- R upper triangular matrix
such that A = QR.
'''
n = A.shape[0]
# base case 1x1 matrix do nothing
if n == 1 : return [1], A
R = A.copy()
# get the vectors x, e and the scalar alpha
x = R[:,0]
e_0 = np.identity(n)[0]
alpha = -np.sign(x[0]) * np.linalg.norm(x)
u = x + alpha*e_0
v = u/np.linalg.norm(u)
# matrix of the reflection x -> x - 2<v,x>v
Q = np.identity(n) - 2*np.array([ v[i]*v for i in range(n)])
R = Q @ R
# do the recursion
Q1, R1 = QR(R[1:,1:])
# copy the results into Q, R
Q[1:,1:] = Q[1:,1:] @ Q1
R[1:,1:] = R1
return Q, R
A = np.array([[12, -51, 4],
[6, 167, -68],
[-4, 24, -41]])
QR(A)
(array([[ 0.85714286, 0.39428571, -0.33142857],
[ 0.42857143, -0.90285714, 0.03428571],
[-0.28571429, -0.17142857, -0.94285714]]),
array([[ 1.40000000e+01, 2.10000000e+01, -1.40000000e+01],
[-5.50670620e-16, -1.75000000e+02, 7.00000000e+01],
[ 3.01980663e-16, -3.63797881e-14, 3.50000000e+01]]))
A = np.identity(3)
A[2:,2:].shape
(1, 1)
import numpy as np
from tabulate import tabulate
# A is a square random matrix of size n
n = 5
A = np.random.rand(n, n)
print("A=")
print(tabulate(A))
def eigen_qr_simple(A, iterations=500000):
Ak = np.copy(A)
n = A.shape[0]
QQ = np.eye(n)
for k in range(iterations):
Q, R = np.linalg.qr(Ak)
Ak = R @ Q
QQ = QQ @ Q
# we "peek" into the structure of matrix A from time to time
# to see how it looks
if k%100000 == 0:
print("A",k,"=")
print(tabulate(Ak))
print("\n")
return Ak, QQ
# We call the function
eigen_qr_simple(A)
# We compare our results with the official numpy algorithm
print(np.linalg.eigvals(A))
A= --------- -------- --------- --------- --------- 0.201161 0.867271 0.812846 0.478103 0.756584 0.821249 0.752742 0.272163 0.784753 0.780868 0.975179 0.562201 0.272817 0.931913 0.663289 0.0999714 0.389918 0.0610292 0.0613589 0.217314 0.39985 0.159648 0.318877 0.549889 0.0273645 --------- -------- --------- --------- --------- A 0 = ---------- ----------- ---------- ---------- ---------- 1.68787 -1.48075 0.156866 -0.900978 -0.726793 -1.16421 -0.0989761 -0.0168406 0.437746 0.0419396 0.308498 0.00277441 -0.147396 0.194095 -0.296607 0.0446465 0.0444515 0.120628 -0.0920303 -0.0155305 0.0386219 0.0246236 0.0823022 -0.0842476 -0.0340297 ---------- ----------- ---------- ---------- ---------- A 100000 = ------- --------- --------- ------------- ---------- 2.37559 -0.191106 -0.907508 -0.857174 0.273345 0 -0.798147 -0.107658 -0.245126 -0.069745 0 0 -0.212483 0.160252 0.246612 0 0 0 -0.104604 -0.28233 0 0 0 4.94066e-324 0.0550875 ------- --------- --------- ------------- ---------- A 200000 = ------- --------- --------- ------------- ---------- 2.37559 -0.191106 -0.907508 -0.857174 0.273345 0 -0.798147 -0.107658 -0.245126 -0.069745 0 0 -0.212483 0.160252 0.246612 0 0 0 -0.104604 -0.28233 0 0 0 4.94066e-324 0.0550875 ------- --------- --------- ------------- ---------- A 300000 = ------- --------- --------- ------------- ---------- 2.37559 -0.191106 -0.907508 -0.857174 0.273345 0 -0.798147 -0.107658 -0.245126 -0.069745 0 0 -0.212483 0.160252 0.246612 0 0 0 -0.104604 -0.28233 0 0 0 4.94066e-324 0.0550875 ------- --------- --------- ------------- ---------- A 400000 = ------- --------- --------- ------------- ---------- 2.37559 -0.191106 -0.907508 -0.857174 0.273345 0 -0.798147 -0.107658 -0.245126 -0.069745 0 0 -0.212483 0.160252 0.246612 0 0 0 -0.104604 -0.28233 0 0 0 4.94066e-324 0.0550875 ------- --------- --------- ------------- ---------- [ 2.3755882 -0.79814693 0.0550875 -0.10460353 -0.2124826 ]
def pivot_matrix(M):
"""Returns the pivoting matrix for M, used in Doolittle's method."""
m = len(M)
# Create an identity matrix, with floating point values
id_mat = [[float(i == j) for i in range(m)] for j in range(m)]
# Rearrange the identity matrix such that the largest element of
# each column of M is placed on the diagonal of of M
for j in range(m):
row = max(range(j, m), key=lambda i: abs(M[i][j]))
if j != row:
# Swap the rows
id_mat[j], id_mat[row] = id_mat[row], id_mat[j]
return id_mat
M = np.array([0,1,1,1]).reshape(2,2)
pivot_matrix(M)
[[0.0, 1.0], [1.0, 0.0]]
def pivot_m(M):
"""Returns the pivoting matrix for M,
used in Doolittle's method."""
m = len(M)
P = np.identity(m, dtype=np.int32)
# Rearrange P such that the largest element of
# each column of M is placed on the diagonal of of M
for j in range(m):
row = max(range(j, m), key = lambda i: abs(M[i,j]))
if j != row:
# Swap the rows
P[[j,row]] = P[[row,j]]
return P
pivot_m(M)
array([[0, 1],
[1, 0]], dtype=int32)