Linear Systems and LU Decomposition
Chapter 2, Linear Systems,
Mainly LU Decomposition
Linear Systems
In matrix notation:
A x = b
See, e.g., “Matrix Analysis and Applied Linear Algebra”, Carl D. Meyer,
SIAM, 2000. Or “Matrix Computations”, Golub and van Loan.
Existence and Solutions
•
M
=
N,
for nonsingular square matrix, a
unique solution exists
•
M
<
N
, more variables than equations –
infinitely many solutions, need to find the
null space of
A
,
A
x
= 0 (e.g. SVD)
•
M
>
N
, one can ask for a least-square
solution, e.g., from the normal equation
(
A
T
A
)
x
= (
A
T
b
)
Some Concepts in Linear Algebra
•
Vector space, linear independence, and
dimension
•
Null space of a matrix
A
: the set of all
x
such that
Ax
= 0
•
Range of
A
: the set of all
Ax
•
Rank of
A
: max number of linearly
independent rows or columns
•
Rank of
A
+ Null space Dim of
A
= number
of columns of
A
Computation Tasks and Pitfalls
•
Solve
A x
=
b
•
Find
A
-1
•
Compute det(
A
)
•
Round off error makes the system
singular, or numerical instability makes the
answer wrong.
Gauss Elimination
•
Basic facts about linear equations
–
Interchanging any two rows of
A
and
b
does
not change the solution
x
–
Replace a row by a linear combination of itself
and any other row does not change
x
–
Interchange column permutes the solution
•
Example of Gauss elimination
•
Pivoting
•
LU equivalent to Gauss elimination
LU Decomposition
L . U
=
A
Thus
A x
= (
LU
)
x
=
L
(
U
x) =
b
Let
y
=
U x
, then solve
y
in
L y
=
b
by
forward substitution
Solve
x
in
U x
=
y
by backward substitution
Crout’s Algorithm
•
Set
for all
i
•
For each
j
= 1, 2, 3, …,
N
(a) for
i
= 1, 2, …,
j
(b) for
i
=
j
+
1
,
j
+2, …,
N
a
ij
for A
for
L
for
U
A = LU
L y = b
U x = y
Ax = LUx = b
Order of Update in Crout’s
Algorithm
LU Decomposition (without
pivoting)
Determining the entries
in an order scanning the
column from top to
bottom and left to right.
L
and
U
are stored in
the same place as the
original
A
.
Pivoting
•
Pivoting is essential for stability
•
Interchange rows of
A
to get largest
β
ii
,
this results
LU
=
P A
.
•
Implicit pivoting (when comparing for the
biggest elements in a column, use the
normalized one so that the largest
coefficient in an equation is 1)
ludcmp(…)
import math
def ludcmp(a, n, indx, d):
d[0] = 1.0
vv=indx.copy()
for i in range(0,n):
big = 0.0
for j in range(0,n):
temp = math.fabs(a[i][j])
if (temp > big):
big = temp
vv[i]=1.0/big
Indentations are
important in Python
Run Crout’s algorithm
for j in range(0,n):
big = 0.0
for i in range(0,n):
if(i<j):
l = i
else:
l = j
sum = a[i][j]
for k in range(0,l):
sum -= a[i][k]*a[k][j]
a[i][j] = sum
if(i>=j):
dum = vv[i]*math.fabs(sum)
if (dum >= big):
big = dum
imax = i
if (j != imax):
dum = a[imax]
a[imax] = a[j]
a[j] = dum
d[0] = - d[0]
vv[imax] = vv[j]
indx[j] = imax
dum = 1.0/a[j][j]
for i in range(j+1,n):
a[i][j] *= dum
Need very carefully
pay attention to
indentation!
Computational Complexity
•
How many basic steps does the LU
decomposition program take?
•
Big O( … ) notation & asymptotic
performance
•
O(
N
3
)
Compute
A
-1
•
Let
B
=
A
-1
then
AB
=
I
(
I
is identity matrix)
or
A
[
b
1
,
b
2
, …,
b
N
] = [e
1
, e
2
, …,e
N
]
or
A
b
j
= e
j
for
j
= 1, 2, …,
N
where
b
j
is the
j
-th column of
B
.
•
I.e., to compute A
-1
, we solve a linear
system
N
times, each with a unit vector e
j
.
Compute det(
A
)
•
Definition of determinant
•
Properties of determinant (anti-symmetry,
Laplace expansion, etc)
•
Since det(
LU
) = det(
L
) det(
U
), thus
det(
A
) = det(
U
) =
Use numpy and scipy
•
numpy: define
N
-dimensional arrays of
various types
•
scipy: efficient numerical routines such as
integration, interpolation, linear algebra,
and statistics covered in this course. scipy
need numpy in most cases.
Solve A
x
=
b
in scipy
import numpy as np
from scipy import linalg
# A is 2x2 matrix
A = np.array([[1.03,2.0],[3.5,4.1]])
# b is the right-side vector
b = np.array([[5.0],[6.1]])
# x is solution
x = np.linalg.solve(A,b)
Use LAPACK
•
Lapack is a free, high quality linear
algebra solver package (downloadable at
www.netlib.org/lapack/
). Much more
sophisticated than NR routines.
•
Written in Fortran 90
•
Calling inside Python made possible
through scipy
What Lapack can do?
•
Solution of linear systems,
Ax
=
b
•
Least-square problem, min ||
Ax
-
b
||
2
•
Singular value decomposition,
A
=
UsV
T
•
Eigenvalue problems,
Ax
=
x
An Example for Lapack Fortran
Routine
SUBROUTINE DSYEVR(JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL, M,
W, Z, LDZ, ISUPPZ, WORK, LWORK, IWORK, LIWORK, INFO)
DSYEVR computes selected eigenvalues and, optionally, eigenvectors of a real symmetric
matrix A. Eigenvalues and eigenvectors can be selected by specifying either a range
of values or a range of indices for the desired eigenvalues.
JOBZ (input) CHARACTER*1 = 'N': Compute eigenvalues only; = 'V': Compute
eigenvalues and eigenvectors.
RANGE (input) CHARACTER*1 = 'A': all eigenvalues will be found. = 'V': all eigenvalues
in the half-open interval (VL,VU] will be found. = 'I': the IL-th through IU-th eigenvalues
will be found. For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and DSTEIN are
called
UPLO (input) CHARACTER*1 = 'U': Upper triangle of A is stored; = 'L': Lower triangle of A
is stored.
N (input) INTEGER The order of the matrix A. N >= 0.
A (input/output) DOUBLE PRECISION array, dimension (LDA, N) On entry, the symmetric
matrix A. If UPLO = 'U', the leading N-by-N upper triangular part of A contains the
upper triangular part of the matrix A. If UPLO = 'L', the leading N-by-N lower triangular
part of A contains the lower triangular part of the matrix A. On exit, the lower triangle (if
UPLO='L') or the upper triangle (if UPLO='U') of A, including the diagonal, is destroyed.
skip
Other matrix decompositions
•
Cholesky:
A
=
LL
T
,
A
is real symmetric,
positive definite.
•
QR factorization:
A=QR, Q
is orthogonal
Q
T
Q=I, R
is upper triangular.
•
Schur:
U
†
AU=
upper triangular with
eigenvalues
λ
on the diagonals.
•
Singular value decomposition:
A=U
Σ
V
†
,
U
and
V
are unitary
U
†
U=I, V
†
V=I
,
Σ
is
nonnegative diagonal.
Reading, Reference
•
Read NR, Chap. 2
•
M. T. Heath, “Scientific Computing: an
introductory survey”
•
For a more thorough treatment on
numerical linear algebra computation
problems, see G H Golub & C F van Loan,
“Matrix Computations”, as well as Horn &
Johnson, “Matrix Analysis.”
P
r
o
b
l
e
m
s
f
o
r
w
e
e
k
3
/
c
h
a
p
t
e
r
2
(
d
u
e
T
h
u
1
2
S
e
p
2
0
2
4
)
1
.
C
o
n
s
i
d
e
r
t
h
e
f
o
l
l
o
w
i
n
g
l
i
n
e
a
r
e
q
u
a
t
i
o
n
(a) Solve the system (by hand) with LU decomposition, following Crout’s
algorithm (slide 8, NR p.44) exactly, without pivoting.
(b) Find the inverse of the 3
3 matrix, using LU decomposition.
(c) Find the determinant of the 3
3 matrix, using LU decomposition.
2
.
C
a
n
w
e
d
o
L
U
d
e
c
o
m
p
o
s
i
t
i
o
n
f
o
r
a
l
l
s
q
u
a
r
e
,
r
e
a
l
m
a
t
r
i
c
e
s
?
W
h
a
t
i
s
t
h
e
c
o
n
d
i
t
i
o
n
f
o
r
t
h
e
e
x
i
s
t
e
n
c
e
o
f
a
n
L
U
d
e
c
o
m
p
o
s
i
t
i
o
n
?
T
h
e
a
n
s
w
e
r
s
d
e
p
e
n
d
i
f
p
i
v
o
t
i
n
g
i
s
a
p
p
l
i
e
d
o
r
n
o
t
.
Chapter number refers to “Numerical Recipes” chapters.
Explore the fundamental concepts of linear algebra, including matrix notation, existence of solutions, vector spaces, computation tasks, and LU decomposition techniques. Learn about Gauss elimination, Crout's algorithm, and how to solve linear systems efficiently using LU decomposition.
Download Presentation
Please find below an Image/Link to download the presentation.
The content on the website is provided AS IS for your information and personal use only. It may not be sold, licensed, or shared on other websites without obtaining consent from the author. Download presentation by click this link. If you encounter any issues during the download, it is possible that the publisher has removed the file from their server.
E N D
Presentation Transcript
Chapter 2, Linear Systems, Mainly LU Decomposition
Linear Systems In matrix notation: A x = b See, e.g., Matrix Analysis and Applied Linear Algebra , Carl D. Meyer, SIAM, 2000. Or Matrix Computations , Golub and van Loan.
Existence and Solutions M = N, for nonsingular square matrix, a unique solution exists M < N, more variables than equations infinitely many solutions, need to find the null space of A, Ax = 0 (e.g. SVD) M > N, one can ask for a least-square solution, e.g., from the normal equation (ATA) x = (AT b)
Some Concepts in Linear Algebra Vector space, linear independence, and dimension Null space of a matrix A: the set of all x such that Ax = 0 Range of A: the set of all Ax Rank of A: max number of linearly independent rows or columns Rank of A + Null space Dim of A = number of columns of A
Computation Tasks and Pitfalls Solve A x = b Find A-1 Compute det(A) Cramer s rule ??=det(??) det(?) Round off error makes the system singular, or numerical instability makes the answer wrong.
Gauss Elimination Basic facts about linear equations Interchanging any two rows of A and b does not change the solution x Replace a row by a linear combination of itself and any other row does not change x Interchange column permutes the solution Example of Gauss elimination Pivoting LU equivalent to Gauss elimination
LU Decomposition 1 = ii L . U = A Thus A x = (LU) x = L (U x) = b Let y = U x, then solve y in L y = b by forward substitution Solve x in U x = y by backward substitution
Crouts Algorithm = 1 Set For each j= 1, 2, 3, , N (a) for i= 1, 2, , j a ij ij = for all i aij for A ii for L for U 1 i A = LU ik kj L y = b U x = y Ax = LUx = b = 1 k (b) for i = j +1, j +2, , N 1 jj 1 j k = a ij ij ik kj = 1
Order of Update in Crouts Algorithm ???= 1
LU Decomposition (without pivoting) 2 1 3 0 5 1 0 1 0 0 1 11 0 0 12 13 = = = 2 1 2 A LU 21 22 0 23 1 31 32 33 Determining the entries in an order scanning the column from top to bottom and left to right. L and U are stored in the same place as the original A. 1 0 1 1/8 1 0 0 2 0 0 1 4 0 5 4 = 1 1/ 2 4
Pivoting Pivoting is essential for stability Interchange rows of A to get largest ii, this results LU = P A. Implicit pivoting (when comparing for the biggest elements in a column, use the normalized one so that the largest coefficient in an equation is 1)
ludcmp() import math def ludcmp(a, n, indx, d): d[0] = 1.0 vv=indx.copy() for i in range(0,n): big = 0.0 for j in range(0,n): temp = math.fabs(a[i][j]) if (temp > big): big = temp vv[i]=1.0/big Indentations are important in Python
Run Crouts algorithm for j in range(0,n): big = 0.0 for i in range(0,n): if(i<j): l = i else: l = j sum = a[i][j] for k in range(0,l): sum -= a[i][k]*a[k][j] a[i][j] = sum
if(i>=j): dum = vv[i]*math.fabs(sum) if (dum >= big): big = dum imax = i if (j != imax): dum = a[imax] a[imax] = a[j] a[j] = dum d[0] = - d[0] vv[imax] = vv[j] indx[j] = imax dum = 1.0/a[j][j] for i in range(j+1,n): a[i][j] *= dum Need very carefully pay attention to indentation!
Computational Complexity How many basic steps does the LU decomposition program take? Big O( ) notation & asymptotic performance O(N3)
Compute A-1 Let B = A-1 then AB = I (I is identity matrix) or A [b1, b2, ,bN] = [e1, e2, ,eN] or Abj = ej for j= 1, 2, , N where bj is the j-th column of B. I.e., to compute A-1, we solve a linear system N times, each with a unit vector ej.
Compute det(A) Definition of determinant = ( 1) P det( ) A a a a 1, 2, , 1 i i n i 2 n P Properties of determinant (anti-symmetry, Laplace expansion, etc) Since det(LU) = det(L) det(U), thus det(A) = det(U) = 1 j N = jj
Use numpy and scipy numpy: define N-dimensional arrays of various types scipy: efficient numerical routines such as integration, interpolation, linear algebra, and statistics covered in this course. scipy need numpy in most cases.
Solve Ax =b in scipy import numpy as np from scipy import linalg # A is 2x2 matrix A = np.array([[1.03,2.0],[3.5,4.1]]) # b is the right-side vector b = np.array([[5.0],[6.1]]) # x is solution x = np.linalg.solve(A,b)
Use LAPACK Lapack is a free, high quality linear algebra solver package (downloadable at www.netlib.org/lapack/). Much more sophisticated than NR routines. Written in Fortran 90 Calling inside Python made possible through scipy
What Lapack can do? Solution of linear systems, Ax = b Least-square problem, min ||Ax-b||2 Singular value decomposition, A = UsVT Eigenvalue problems, Ax = x
An Example for Lapack Fortran Routine SUBROUTINE DSYEVR(JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, ISUPPZ, WORK, LWORK, IWORK, LIWORK, INFO) DSYEVR computes selected eigenvalues and, optionally, eigenvectors of a real symmetric matrix A. Eigenvalues and eigenvectors can be selected by specifying either a range of values or a range of indices for the desired eigenvalues. JOBZ (input) CHARACTER*1 = 'N': Compute eigenvalues only; = 'V': Compute eigenvalues and eigenvectors. RANGE (input) CHARACTER*1 = 'A': all eigenvalues will be found. = 'V': all eigenvalues in the half-open interval (VL,VU] will be found. = 'I': the IL-th through IU-th eigenvalues will be found. For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and DSTEIN are called UPLO (input) CHARACTER*1 = 'U': Upper triangle of A is stored; = 'L': Lower triangle of A is stored. N (input) INTEGER The order of the matrix A. N >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA, N) On entry, the symmetric matrix A. If UPLO = 'U', the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A. If UPLO = 'L', the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A. On exit, the lower triangle (if UPLO='L') or the upper triangle (if UPLO='U') of A, including the diagonal, is destroyed. skip
Other matrix decompositions Cholesky: A = LLT, A is real symmetric, positive definite. QR factorization: A=QR, Q is orthogonal QTQ=I, R is upper triangular. Schur: U AU= upper triangular with eigenvalues on the diagonals. Singular value decomposition: A=U V , U and V are unitary U U=I, V V=I , is nonnegative diagonal.
Reading, Reference Read NR, Chap. 2 M. T. Heath, Scientific Computing: an introductory survey For a more thorough treatment on numerical linear algebra computation problems, see G H Golub & C F van Loan, Matrix Computations , as well as Horn & Johnson, Matrix Analysis.
Problems for week 3/chapter 2 (due Thu 12 Sep 2024) 1. Consider the following linear equation 1 2 5 0 1 2 x x x 1 = 1 2 3 6 0 1 2 3 (a) Solve the system (by hand) with LU decomposition, following Crout s algorithm (slide 8, NR p.44) exactly, without pivoting. (b) Find the inverse of the 3 3 matrix, using LU decomposition. (c) Find the determinant of the 3 3 matrix, using LU decomposition. 2. Can we do LU decomposition for all square, real matrices? What is the condition for the existence of an LU decomposition? The answers depend if pivoting is applied or not.
