Lecture 3
Numerical methods
System of linear equations
Lecture 3
OUTLINE
1. System of linear equations (linear system)
2. Gaussian elimination
3. Pivoting
4. LU decomposition
5. Cholesky decomposition
6. Effect of round-off errors
7. Conditionality
8. System of homogeneous linear equations
System of linear equations (linear system)
Basic definitions
x1a11 + + xn a1n = b1
x1a 21 + + xn a 2 n = b2
(1)
x1a m1 + + xn a mn = bm
We can assign two matrixes to the system (1):
- coefficient matrix - augmented matrix
æ a11 a12 a1n ö÷ æ a11 a12 a1n b1 ö÷
çç ÷÷ çç ÷÷
çç a 21 a 22 a 2n ÷÷ çç a 21 a 22 a 2n b2 ÷÷
A = çç ÷÷ A ¢ = çç ÷÷
çç ÷÷ çç ÷÷
çç ÷÷ çç ÷÷
èa m1 a m2 a mn ø èa m1 a m2 a mn bm ø
System of linear equations (linear system)
Denote
æ x1 ö æ b1 ö
çç ÷÷ çç ÷÷
çç x2 ÷÷ çç b2 ÷÷
x = çç ÷÷÷ b = çç ÷÷÷
çç ÷÷ çç ÷÷
çç ÷÷ çç ÷÷
è xn ø èbm ø
Then the linear system (1) reads
æ a11 a1n ÷öçæ x1 ÷ö çæ b1 ÷ö
ççç ÷ç ÷÷ ç ÷÷
÷
çç a 21 a 2n ÷÷çç x2 ÷÷ çç b2 ÷÷
÷çç ÷÷ = çç ÷÷
ççç ÷
÷÷çç ÷÷ ç ÷÷
çèça m1 ÷÷çç ÷÷ ççç ÷÷
a mn øè xn ø èbm ø
Ax = b
System of linear equations (linear system)
0
0
If b= then the system (1) is called homogeneous
0
0
0
if b then it is called nonhomogeneous
0
System of linear equations (linear system)
Definition We say that the n-tuple (r1, …, rn) is solution of system (1),
if Ar = b, where
r1
r= r2 .
r
n
Note
System (1) is consistent,
if it has at least one solution;
otherwise it is inconsistent.
System of linear equations (linear system)
Direct methods
attempt to solve the problem
by a finite sequence of operations.
In the absence of rounding errors,
direct methods would deliver an exact solution.
Iterative methods
delivers only approximate solution.
The number of steps depends
on required precision.
Example: Find the solution of linear system in real numbers.
x-2y+4z+ t =-6
2x+3y- z+2t =13 | (2)-2*(1)
2x+5y+ z+ t = 8 | (3)-2*(1)
3x+ y+3z+ t = 1 | (4)-3*(1)
x-2y+4z+ t =-6
7y-9z =25
9y-7z- t =20 | (3)-9/7*(2)
7y-9z-2t =19 | (4)-7/7*(2)
x-2y+ 4z+ t = -6
7y -9z = 25
32/7z- t =-85/7
-2t = -6
x = 1, y = 1, z = - 2, t = 3
Elementary row operations
Definition: Elementary row operations (ERO)
on a matrix are:
– row switching ↔
– row multiplication →
– row addition →
Row equivalence of matrices
Definition: We say that matrix A of type m x n is
row equivalent with the matrix B of type m x n
if it is possible to transform A into B
by a sequence of elementary row operations.
Upper triangular matrix
Definition: We say that squared matrix U is
upper triangular (or right triangular)
if all the entries below the main diagonal are zero.
Row equivalence of matrices
Theorem: Each matrix is row equivalent
with some upper triangular matrix.
Rank of a matrix
Definition: The rank of matrix A is
the number of nonzero rows
of triangular matrix equivalent
to the matrix A,
we will denote it h(A).
or
The rank of A is the maximal number
of linearly independent rows of A.
Consequence: Row equivalent matrices have the same rank.
Example: Find the solution of linear system in real numbers.
x-2y+4z+ t =-6 1 2 4 1 6
2x+3y- z+2t =13 /(2)-2*(1) 2 3 1 2 13
2x+5y+ z+ t = 8 /(3)-2*(1) 2 5 1 1 8
3x+ y+3z+ t = 1 /(4)-3*(1) 3 1 3 1 1
æ1 - 2 4 1 -6 ö÷
x-2y+4z+ t = -6 çç ÷÷
7y-9z = 25 çç0 7 -9 0 25÷÷
çç ÷÷
9y-7z- t = 20 /(3)-9/7*(2) çç0 9 -7 -1 20÷÷
çç ÷÷
7y-9z-2t = 19 /(4)-7/7*(2) è 0 7 -9 - 2 19ø
æ1 - 2 4 1 -6 ö
x-2y+ 4z+ t = -6 çç ÷÷
7y- 9z = 25 ççç0 7 -9 0 25 ÷÷÷
32/7z- t =-85/7 çç0 0 32 / 7 -1 -85/7÷÷÷
çç ÷÷
-2t = -6 çè0 0 0 -2 -6 ÷ø
x = 1, y = 1, z = - 2, t = 3
Gaussian elimination
We showed a principle of
Gaussian elimination (GE):
Using the finite number of elementary row operations
on the coefficient matrix we get
the upper triangular matrix
and
using back substitution we calculate the vector of unknowns.
Gaussian elimination
We showed a principle of
Gaussian elimination (GE):
Using the finite number of elementary row operations
on the coefficient matrix we get
the upper triangular matrix
and
using back substitution we calculate the vector of unknowns.
Forward elimination:
Let
A(0) = A, b(0) = b
elements of matrix A(0) will be aij(0) = aij
elements of vector b(0) will be bi(0) = bi
Gaussian elimination
algorithm of GE
Gaussian elimination
algorithm of GE
multiplicator assures
that
element (i,k) of matrix A(k)
will be zero
Gaussian elimination
algorithm of GE
leading coefficient
(pivot)
has to be non-zero
Gaussian elimination
Example:
Solve the linear system
on the hypothetical computer,
which works in decimal system
with five digits mantissa.
The exact solution is
Gaussian elimination
Gaussian elimination
Gaussian elimination – partial pivoting
Partial pivoting
is the modification of GE which assures
that the absolute values of multiplicators
are less than or equal to 1.
The algorithm selects the entry
with largest absolute value
from the column of the matrix
that is currently being considered
as the pivot element.
Gaussian elimination – partial pivoting
Gaussian elimination – complete (maximal) pivoting
Pivoting
Partial or complete pivoting?
usually partial pivoting is enough
Partial pivoting is generally sufficient to adequately reduce round-off error.
Complete pivoting is usually not necessary
to ensure numerical stability and,
due to the additional cost
of searching for the maximal element,
the improvement in numerical stability
is typically outweighed
by its reduced efficiency
for all but the smallest matrices.
Gaussian elimination
Computational efficiency
2 3
The number of arithmetic operations in forward elimination: » n
3
The number of arithmetic operations in back substitution : » n2
LU decomposition (factorization)
x-2y+4z+ t =-6 1 2 4 1 6
2x+3y- z+2t =13 /(2)-2*(1) 2 3 1 2 13
2x+5y+ z+ t = 8 /(3)-2*(1) 2 5 1 1 8
3x+ y+3z+ t = 1 /(4)-3*(1) 3 1 3 1 1
æ1 - 2 4 1 -6 ö÷
x-2y+4z+ t = -6 çç ÷÷
7y-9z = 25 çç0 7 -9 0 25÷÷
çç ÷÷
9y-7z- t = 20 /(3)-9/7*(2) çç0 9 -7 -1 20÷÷
çç ÷÷
7y-9z-2t = 19 /(4)-7/7*(2) è 0 7 -9 - 2 19ø
æ1 - 2 4 1 -6 ö
x-2y+ 4z+ t = -6 çç ÷÷
7y- 9z = 25 ççç0 7 -9 0 25 ÷÷÷
32/7z- t =-85/7 çç0 0 32 / 7 -1 -85/7÷÷÷
çç ÷÷
-2t = -6 çè0 0 0 -2 -6 ÷ø
x = 1, y = 1, z = - 2, t = 3
LU decomposition (factorization)
x-2y+4z+ t =-6 1 2 4 1 6
2x+3y- z+2t =13 /(2)-2*(1) 2 3 1 2 13
2x+5y+ z+ t = 8 /(3)-2*(1) 2 5 1 1 8
3x+ y+3z+ t = 1 /(4)-3*(1) 3 1 3 1 1
æ 1 -2 4 1 -6 ö÷
x-2y+4z+ t = -6 çç ÷÷
7y-9z = 25 çç2 7 -9 0 25÷÷
çç ÷÷
9y-7z- t = 20 /(3)-9/7*(2) çç2 9 -7 -1 20÷÷
çç ÷÷
7y-9z-2t = 19 /(4)-7/7*(2) è 3 7 -9 - 2 19ø
æ 1 -2 4 1 -6 ö
x-2y+ 4z+ t = -6 çç ÷÷
7y- 9z = 25 çç2 7 -9 0 25 ÷÷÷
çç ÷÷
32/7z- t =-85/7 çç2 0 32 / 7 -1 -85/7÷÷
-2t = -6 çç ÷÷
è3 0 0 -2 -6 ø
x = 1, y = 1, z = - 2, t = 3
LU decomposition (factorization)
x-2y+4z+ t =-6 1 2 4 1 6
2x+3y- z+2t =13 /(2)-2*(1) 2 3 1 2 13
2x+5y+ z+ t = 8 /(3)-2*(1) 2 5 1 1 8
3x+ y+3z+ t = 1 /(4)-3*(1) 3 1 3 1 1
æ 1 -2 4 1 -6 ö÷
x-2y+4z+ t = -6 çç ÷÷
7y-9z = 25 çç2 7 -9 0 25÷÷
çç ÷÷
9y-7z- t = 20 /(3)-9/7*(2) çç2 9 -7 -1 20÷÷
çç ÷÷
7y-9z-2t = 19 /(4)-7/7*(2) è 3 7 -9 - 2 19ø
æ 1 -2 4 1 -6 ö
x-2y+ 4z+ t = -6 çç ÷÷
7y- 9z = 25 çç2 7 -9 0 25 ÷÷÷
çç ÷÷
32/7z- t =-85/7 çç2 9 / 7 32 / 7 -1 -85/7÷÷
-2t = -6 çç ÷÷
è3 1 0 -2 -6 ø
x = 1, y = 1, z = - 2, t = 3
LU decomposition (factorization)
x-2y+4z+ t =-6 1 2 4 1 6
2x+3y- z+2t =13 /(2)-2*(1) 2 3 1 2 13
2x+5y+ z+ t = 8 /(3)-2*(1) 2 5 1 1 8
3x+ y+3z+ t = 1 /(4)-3*(1) 3 1 3 1 1
æ 1 -2 4 1 -6 ö÷
x-2y+4z+ t = -6 çç ÷÷
7y-9z = 25 çç2 7 -9 0 25÷÷
çç ÷÷
9y-7z- t = 20 /(3)-9/7*(2) çç2 9 -7 -1 20÷÷
çç ÷÷
7y-9z-2t = 19 /(4)-7/7*(2) è 3 7 -9 - 2 19ø
æ 1 -2 4 1 -6 ö
x-2y+ 4z+ t = -6 çç ÷÷
7y- 9z = 25 çç2 7 -9 0 25 ÷÷÷
çç ÷÷
32/7z- t =-85/7 çç2 9 / 7 32 / 7 -1 -85/7÷÷
-2t = -6 çç ÷÷
è3 1 0 -2 -6 ø
x = 1, y = 1, z = - 2, t = 3
LU decomposition (factorization)
Multiplicators mij from the forward elimination are placed
into lower triangular matrix L)
LU decomposition (factorization)
We obtained LU decomposition of matrix A
(so called Doolittle version – ones on the diagonal of matrix L)
because
A = L⋅ U
æ1 -2 4 1ö æ1 0 0 0ö æ1 -2 4 1 ö
çç ÷÷ ç ÷÷ ç ÷÷
ç
çç2 3 -1 2÷÷ çç2 1 0 0÷÷ çç0 7 ç ÷÷
÷÷ = ç ÷÷⋅ ç - 9 0 ÷÷
çç
çç2 5 1 1÷÷÷ çç2 9 / 7 1 0÷÷÷ çç0 0 32 / 7 -1÷÷÷
çç ÷÷ ççç ÷÷ ççç ÷÷
3 1 3 1ø è3 1 0 1ø è0 0
è 0 -2
ø
A L U
Doolittle algorithm
makes the LU decomposition
using the half number
of arithmetic operations
1 3
» n
3
LU decomposition (factorization)
LU decomposition of A
is useful for solving of sequences of tasks
if the new right-hand-side bi could be estimated only after
we solve the previous linear systems
for k < i.
To solve the system means to solve
LU decomposition with partial pivoting
LU decomposition with partial pivoting
is standard routine of each library of codes.
Input parameter is matrix A.
Output are three matrices:
L, U and permutation matrix P
where
LU = PA.
Permutation matrix comes form identity matrix
after row exchange.
Solution of Ax=b is given by solution of
(show this)
LU decomposition with partial pivoting
Instead of matrix P we work with
the vector of row permutations p.
At the beginning we put
and we just exchange the elements of vector p.
Gaussian elimination
Definition: We say that square matrix A
is diagonally dominant if
Gaussian elimination
Definition: Symmetric matrix A is positive-definite if
for any non-zero vector x it holds
xT ⋅ A ⋅ x > 0
If A is regular then
AT ⋅ A
is positive-definite.
Gaussian elimination
Sylvester’s criterion: Symmetric matrix A
is positive-definite, if and only if the leading
principal minors are positive, i.e. if
Gaussian elimination
Gaussian elimination
is numerically stable for
diagonally dominant
or
positive-definite matrices.
Cholesky decomposition
Cholesky decomposition theorem: There is only one way, how the
decompose the positive-definite matrix A into
A = L ⋅ LT
where L is the lower triangular matrix.
Cholesky decomposition
Cholesky decomposition theorem: There is only one way, how the
decompose the positive-definite matrix A into
A = L ⋅ LT
where L is the lower triangular matrix.
Cholesky algorithm:
1 3
The number of arithmetic operations: » n
3
The effect of round-off errors
Example:
On the hypothetical computer working decimal system
with three digits mantissa solve the following system
The effect of round-off errors
GE with partial pivoting
assures the small residua
However, the small residuum does not give
an estimation of the error of solution
Conditionality
In order to evaluate the conditionality of problem
„to find solution x of system Ax = b“
we need to asses
the effect of change of A a b
to the solution x.
Conditionality
Norm of the vector
class of vector norms lp,
The most often we use p=1, p=2, or p
Manhattan Euclid Chebyshev
Conditionality
Properties
1. av = | a | v , (absolute homogeneity or absolute scalability)
2. u+v ≤ u v (triangle inequality or subadditivity)
3. If v = 0 then v is the zero vector (separates points)
Conditionality
Condition number of a matrix
Let
where maximum and minimum is considered for all nonzero x.
Ratio M/m is called condition number of matrix A:
Conditionality
Lets consider linear system
Ax = b
and other system, which could be obtained by changing the r.h.s.:
, or
we can consider as an error of b
and
corresponding error of solution x.
Conditionality
Because then
from the definition of M and m
it follows that
so for
where is a residuum.
Condition number of matrix acts as
an amplifier of relative error !
It is possible to show, that also changes in matrix A
leads to the similar effects on results.
Conditionality
Norm of matrix
the number M is also known as the norm of matrix
and also it holds
Equivalently we could write the condition number of matrix as
Conditionality
“Entrywise” norm of matrix
1
æ n n pöp
÷÷
A = vec (A) p = çççåå aij ÷÷
p
çè j=1 j=1 ÷ø
n n
A F
= åå ij
a 2
Frobenius norm (consistent with l2)
i=1 j =1
A ¥
= max aij max norm
i, j
Example: Find the solution of linear system in real numbers.
x+2y+ z- t = 1
2x+3y- z+2t = 3
4x+7y+ z = 5
5x+7y-4z+7t =10
Example: Find the solution of linear system in real numbers.
1 2 1 1 1
x+2y+ z- t = 1
2x+3y- z+2t = 3 2 3 1 2 3
4x+7y+ z = 5 4 7 1 0 5
5x+7y-4z+7t =10 5 7 -4 7 10
1 2 1 1 1
0 1 3 4 1
0 1 3 4 1
0 3 9 12 5
x+2y+ z- t = 1 1 2 1 1 1
-y-3z+4t = 1 0 1 3 4 1
0z+0t = 2 0 0 0 0 2
0t = 0 0 0 0 0 0
system has no solution
Another example: Find the solution of linear system in rational numbers.
x-2y+4z+ t =-6
2x+3y- z+2t =13
3x+ y+3z+3t = 7
x+5y-5z+ t =19
Another example: Find the solution of linear system in rational numbers.
1 2 4 1 -6
x-2y+4z+ t =-6
2x+3y- z+2t =13 2 3 1 2 13
3x+ y+3z+3t = 7 3 1 3 3 7
x+5y-5z+ t =19 1 5 -5 1 19
1 2 4 1 -6
0 7 9 0 25
0 7 9 0 25
0 7 9 0 25
x-2y+4z+ t =-6 1 2 4 1 - 6
7y-9z+0t =25 0 7 9 0 25
0z+0t = 0 0 0 0 0 0
0t = 0 0 0 0 0 0
Another example: Find the solution of linear system in rational numbers.
1 2 4 1 - 6
x-2y+4z+ t =-6
7y-9z+0t =25 0 7 9 0 25
0z+0t = 0 0 0 0 0 0
0t = 0 0 0 0 0 0
System has infinite number of solutions
t = p, z = q, p,q Q
25 9 8 10
y q x qp
7 7 7 7
8 10 25 9
S q p, q, q, p , p, q Q
7 7 7 7
2 34
For example: if q = 1, p = 0, we get particular solution , ,1,0
7 7
System of linear equations
Frobenius theorem: System of nonhomogeneous equations
has solution
if and only if,
the rank of coefficient matrix
is equal to
the rank of augmented matrix.
Consequence 1: If h(A) = h(A´) = n (n is the number of unknowns),
then the system has just one solution.
Consequence 2: If h(A) = h(A´) < n (n is the number of unknowns),
then the system has infinite number of solutions
and
n-h unknowns could be freely chosen.
Consequence 3: If h(A) h(A´) then the system has no solution.
System of homogeneous linear equations
System of homogeneous linear equations
x1a11 x n a1n 0
x1a 21 x n a 2n 0
x1a m1 x n a mn 0
The system has always solution because h(A)=h(A´)
Theorem: System of homogeneous equations has
nontrivial solution
if and only if h(A) n.
(trivial solution is (0,0,...,0) )
Example: Find the solution of linear system in rational numbers.
x-2y+4z+ t =0 1 2 4 1 0
2x+3y- z+2t =0 2 3 1 2 0
3x+ y+3z+3t =0 3 1 3 3 0
x+5y-5z+ t =0 1 5 -5 1 0
1 2 4 1 0
0 7 9 0 0
0 7 9 0 0
0 7 9 0 0
x-2y+4z+ t =0 1 2 4 1 0
7y-9z+0t =0 0 7 9 0 0
0z+0t =0 0 0 0 0 0
0 0 0 0 0
0t =0
Example: Find the solution of linear system in rational numbers.
x-2y+4z+ t =0
7y-9z+0t =0
0z+0t =0
0t =0
t = p, z = q, p,q Q
9
y q
7
10
x qp
7
10 9
S q p, q, q, p , p, q Q
7 7
Summary of notions
• homogeneous and nonhomogeneous system of linear
equations
• coefficient matrix and augmented matrix
• elementary row operations
• row equivalent matrices
• rank of a matrix
• upper and lower triangular matrix
• Gaussian elimination (pivoting)
• LU decomposition (Doolittle method)
• Cholesky decomposition
• Frobenius theorem
• solution of homogeneous linear systems
