Lab 4B: Numerical Linear Algebra
Lab 4B: Numerical Linear Algebra
Due: Sunday, March 26 th at 11:59PM
In this lab you will first explore the numerical stability of several variations of LU factorization.
1. LU factorization with no pivoting will decompose a matrix into a lower triangular matrix
L, and upper triangular matrix U, such that
A = LU .
Pivot selection: The pivot at the ith stage will be A i,i
Stability and efficiency: This method is so unstable that it can only be used in special
circumstances.
In Matlab: Matlab always uses partial pivoting. The command [L, U] = lu(A) might seem
to ignore partial pivoting, but that is wrong. You will have to write your own routine for
LU factoring. The program lutx.m from Moler’s book is available on OWL. Moler’s program
uses partial pivoting. On line 16 of lutx , the pivot is selected. Read Matlab help about the
find command, and replace max with an appropriate command.
2. LU factorization with partial pivoting will also compute a decomposition of a matrix A
into a lower triangular matrix L and upper triangular matrix U, except this time rows can
be interchanged. The row interchanges are recorded in a permutation matrix P. This gives a
factorization
PA = LU .
Pivot selection: The pivot at the ith stage will be the largest element in the ith column on
or below the ith row.
Stability and efficiency: This is an efficient method that is usually stable but in rare
instances can be unstable.
In Matlab: [L, U, P] = lu(A)
3. LU factorization with complete pivoting allows both row and columns to be interchanged.
Two permutation matrices record these interchanges. P records the row interchanges and Q
records column interchanges. The factorization looks like
PAQ = LU .
Pivot selection: The largest (in magnitude) element of corner submatrix (row and column
≥ i) is selected to be the pivot.
Lab 4B: Numerical Linear Algebra 2
Stability and efficiency: This method is stable but more expensive than the other methods.
Finding the largest entry in a matrix can be a slow, memory intensive operation, especially
for large matrices.
In Matlab: Matlab only offers complete pivoting for matrices declared as sparse, which are
not considered in this course. If allowed, the command is [L, U, P, Q] = lucp(Asparse) .
For this lab, see OWL for the lucp.m file to obtain a complete-pivoting factoring.
Problem 1: LU factorization
In this problem you will compare the backward error of each of the three methods. To do this you
will generate a large random sample of matrices and compute the backward error of the result from
each method.
The random matrices you use should have the entries uniformly sampled from [−1,1]. The rand
function in Matlab will be helpful for this. Since the differences in the stability of each of the
methods is most noticeable when the entries vary over many orders of magnitude, you should
replace a random entry in each random matrix with 10 10 . The randi function will be helpful here.
The backward error from each of the 4 methods is:
• No pivoting: ||LU − A|| 2
• Partial pivoting: ||LU − PA|| 2
• Complete pivoting: ||LU − PAQ|| 2
For 200 n × n matrices, compute the average backward error for each of the methods. You will
need to do this 20 times for n ranging from 1 to 20. Plot the resulting averages. You should have a
plot where the x axis is n (ranging from 1 to 20) and the y axis is the backward error. Don’t forget
axis labels, a title and a legend. Comment on this plot in your report. Is it what you expected?
Problem 2: SVD and condition numbers
In this problem, you will check the operation of the svd command, and explore the condition of
linear solving, i.e., solving Ax = b.
Eigenvectors and values and the SVD
We shall use plots to explore eigenvectors geometrically.
• Create a random 2 × 2 matrix, A.
• Use the eig command to calculate the eigenvectors and eigenvalues. If you get complex values,
try another matrix.
• Plot the two eigenvectors.
• On the same axes, plot a unit circle.
Lab 4B: Numerical Linear Algebra 3
• Transform the circle and the eigenvectors using the matrix A.
• Plot the transformed circle and eigenvectors.
• Comment on what you observe.
• Some textbooks like to draw the original figure and the transformed one on the same axes,
so combine both figures into one.
• Hint: You may appreciate your figures more if you use the axis equal command.
Now explore the SVD.
• Use the same matrix as before. Now use the Matlab svd command to obtain SVD factors of
A.
• Verify the properties of the SVD, including multiplying the factors to get A, the orthonormal
properties of some matrices, the relation between the singular values and the eigenvalues of
another matrix.
• Draw the equivalent diagrams for the SVD that you did for the eigenvectors. Thus draw a
circle and the singular vectors V , then the transformed circle and the singular vectors U.
• Hint: You may like to use the singular values to adjust the lengths of the singular vectors.
Finally, we construct a problem for which the condition number is realized.
• Create a random matrix A that is 20 × 20, at least.
• Create a random vector b that is 20 × 1.
• Obtain the SVD of A.
• Set b equal to the first column of u from the matlab SVD.
• Set δb equal to 10 −3 of the last column of u.
• Solve Ax = b and Aδx = δb, and calculate the relative errors.
• Use the Matlab cond command to obtain the condition number.
• Comment on the result, including an explanation of how this choice of b and δb leads to the
results obtained.
Problem 3: Condition number for eigenvalues
Here we explore the theory on eigenvalue condition number.
We shall use two test matrices:
A =
5 9 9
−72 −79 −72
90 90 80
B =
2 −1 4
−1 1 −3
4 −3 5
.
We first check how Matlab functions implement the textbook definitions.
Lab 4B: Numerical Linear Algebra 4
• Read Moler’s discussion of eigen-condition number, and the Matlab help on condeig .
• For the test matrices, use the eig command to obtain the eigensystems. Hint: eig can obtain
both left and right eigenvectors.
• Use condeig to obtain condition numbers for eigenvalues and compare them against condition
numbers calculated from the definitions in §10.6 in Moler’s book.
Now perturb the test matrices:
A 1 =
5 9.01 9
−72 −79 −72
90 90 80
A 2 =
5 9 9
−72 −79.01 −72
90 90 80
A 3 =
5 9 9
−72 −79 −72
90 90.01 80
and similarly for B. Compute the new eigenvalues in each case and see how accurate the analysis
in Moler is.
Submission Requirements
All files should be submitted on OWL on the Lab 4B submission. Only .m , .mlx , and .pdf files
will be accepted for grading. Submit all files needed to reproduce the results you found in the lab.
Your report should provide a discussion of how you solved the problems and your interpretation
of the results. Your discussion should be able to convince the marker that you have a thorough
understanding of the problems and results.
If you submit your report as a Matlab Live Script you will receive 2 bonus marks.
Due: Sunday, March 26 th at 11:59PM
In this lab you will first explore the numerical stability of several variations of LU factorization.
1. LU factorization with no pivoting will decompose a matrix into a lower triangular matrix
L, and upper triangular matrix U, such that
A = LU .
Pivot selection: The pivot at the ith stage will be A i,i
Stability and efficiency: This method is so unstable that it can only be used in special
circumstances.
In Matlab: Matlab always uses partial pivoting. The command [L, U] = lu(A) might seem
to ignore partial pivoting, but that is wrong. You will have to write your own routine for
LU factoring. The program lutx.m from Moler’s book is available on OWL. Moler’s program
uses partial pivoting. On line 16 of lutx , the pivot is selected. Read Matlab help about the
find command, and replace max with an appropriate command.
2. LU factorization with partial pivoting will also compute a decomposition of a matrix A
into a lower triangular matrix L and upper triangular matrix U, except this time rows can
be interchanged. The row interchanges are recorded in a permutation matrix P. This gives a
factorization
PA = LU .
Pivot selection: The pivot at the ith stage will be the largest element in the ith column on
or below the ith row.
Stability and efficiency: This is an efficient method that is usually stable but in rare
instances can be unstable.
In Matlab: [L, U, P] = lu(A)
3. LU factorization with complete pivoting allows both row and columns to be interchanged.
Two permutation matrices record these interchanges. P records the row interchanges and Q
records column interchanges. The factorization looks like
PAQ = LU .
Pivot selection: The largest (in magnitude) element of corner submatrix (row and column
≥ i) is selected to be the pivot.
Lab 4B: Numerical Linear Algebra 2
Stability and efficiency: This method is stable but more expensive than the other methods.
Finding the largest entry in a matrix can be a slow, memory intensive operation, especially
for large matrices.
In Matlab: Matlab only offers complete pivoting for matrices declared as sparse, which are
not considered in this course. If allowed, the command is [L, U, P, Q] = lucp(Asparse) .
For this lab, see OWL for the lucp.m file to obtain a complete-pivoting factoring.
Problem 1: LU factorization
In this problem you will compare the backward error of each of the three methods. To do this you
will generate a large random sample of matrices and compute the backward error of the result from
each method.
The random matrices you use should have the entries uniformly sampled from [−1,1]. The rand
function in Matlab will be helpful for this. Since the differences in the stability of each of the
methods is most noticeable when the entries vary over many orders of magnitude, you should
replace a random entry in each random matrix with 10 10 . The randi function will be helpful here.
The backward error from each of the 4 methods is:
• No pivoting: ||LU − A|| 2
• Partial pivoting: ||LU − PA|| 2
• Complete pivoting: ||LU − PAQ|| 2
For 200 n × n matrices, compute the average backward error for each of the methods. You will
need to do this 20 times for n ranging from 1 to 20. Plot the resulting averages. You should have a
plot where the x axis is n (ranging from 1 to 20) and the y axis is the backward error. Don’t forget
axis labels, a title and a legend. Comment on this plot in your report. Is it what you expected?
Problem 2: SVD and condition numbers
In this problem, you will check the operation of the svd command, and explore the condition of
linear solving, i.e., solving Ax = b.
Eigenvectors and values and the SVD
We shall use plots to explore eigenvectors geometrically.
• Create a random 2 × 2 matrix, A.
• Use the eig command to calculate the eigenvectors and eigenvalues. If you get complex values,
try another matrix.
• Plot the two eigenvectors.
• On the same axes, plot a unit circle.
Lab 4B: Numerical Linear Algebra 3
• Transform the circle and the eigenvectors using the matrix A.
• Plot the transformed circle and eigenvectors.
• Comment on what you observe.
• Some textbooks like to draw the original figure and the transformed one on the same axes,
so combine both figures into one.
• Hint: You may appreciate your figures more if you use the axis equal command.
Now explore the SVD.
• Use the same matrix as before. Now use the Matlab svd command to obtain SVD factors of
A.
• Verify the properties of the SVD, including multiplying the factors to get A, the orthonormal
properties of some matrices, the relation between the singular values and the eigenvalues of
another matrix.
• Draw the equivalent diagrams for the SVD that you did for the eigenvectors. Thus draw a
circle and the singular vectors V , then the transformed circle and the singular vectors U.
• Hint: You may like to use the singular values to adjust the lengths of the singular vectors.
Finally, we construct a problem for which the condition number is realized.
• Create a random matrix A that is 20 × 20, at least.
• Create a random vector b that is 20 × 1.
• Obtain the SVD of A.
• Set b equal to the first column of u from the matlab SVD.
• Set δb equal to 10 −3 of the last column of u.
• Solve Ax = b and Aδx = δb, and calculate the relative errors.
• Use the Matlab cond command to obtain the condition number.
• Comment on the result, including an explanation of how this choice of b and δb leads to the
results obtained.
Problem 3: Condition number for eigenvalues
Here we explore the theory on eigenvalue condition number.
We shall use two test matrices:
A =
5 9 9
−72 −79 −72
90 90 80
B =
2 −1 4
−1 1 −3
4 −3 5
.
We first check how Matlab functions implement the textbook definitions.
Lab 4B: Numerical Linear Algebra 4
• Read Moler’s discussion of eigen-condition number, and the Matlab help on condeig .
• For the test matrices, use the eig command to obtain the eigensystems. Hint: eig can obtain
both left and right eigenvectors.
• Use condeig to obtain condition numbers for eigenvalues and compare them against condition
numbers calculated from the definitions in §10.6 in Moler’s book.
Now perturb the test matrices:
A 1 =
5 9.01 9
−72 −79 −72
90 90 80
A 2 =
5 9 9
−72 −79.01 −72
90 90 80
A 3 =
5 9 9
−72 −79 −72
90 90.01 80
and similarly for B. Compute the new eigenvalues in each case and see how accurate the analysis
in Moler is.
Submission Requirements
All files should be submitted on OWL on the Lab 4B submission. Only .m , .mlx , and .pdf files
will be accepted for grading. Submit all files needed to reproduce the results you found in the lab.
Your report should provide a discussion of how you solved the problems and your interpretation
of the results. Your discussion should be able to convince the marker that you have a thorough
understanding of the problems and results.
If you submit your report as a Matlab Live Script you will receive 2 bonus marks.
2018-03-17
In this lab you will first explore the numerical stability of several variations of LU factorization.