讲解MATLAB程序设计、辅导data编程、MATLAB语言编程调试 辅导留学生 Statistics统计、回归、迭代|辅导留学生 Statistics统计、回
- 首页 >> Java编程 MATLAB PROJECT 4
Please read the Instructions located on the Assignments page prior
to working on the Project.
BEGIN with creating a Live Script Project4.
Note: All exercises in this project have to be completed in the Live Script using the Live
Editor. Please refer to the MATLAB video that explains how to use the Live Script:
https://www.mathworks.com/videos/using-the-live-editor-117940.html?s_tid=srchtitle
The final script has to be generated by exporting the Live Script to PDF.
Each exercise has to begin with the line
Exercise#
You should also mark down the parts such as (a), (b), (c), and etc. This makes grading easier.
Important: we use the default format short for the numbers in all exercises unless it is
specified otherwise. We do not employ format rat since it may cause problems with
running the codes and displaying matrices in the Live Script. If format long had been used
in an exercise, please make sure to return to the default format in the next exercise.
Part I. Evaluating Eigenvalues
Exercise 1 (4 points): Difficulty: Easy
In this exercise, you will calculate the eigenvalues of a matrix by using various methods and
functions – one of the functions, quer, was created earlier in Exercise 4 of Project 3, and the
other functions will be MATLAB built-in functions.
NOTE: The function quer from Project 3 can be used for approximating eigenvalues by
consecutive iterations, but not for all matrices the process converges. For computing
eigenvalues in MATLAB, especially complex eigenvalues, it is advisable to use a MATLAB
built-in function eig().
**Create a function that begins with
function [q,e,c,S]=eigenval(A)
format
The input is a square matrix A. The outputs will be specified below when you proceed with
the programming.
**First Method: your function will calculate the eigenvalues by the QR factorization. Include
in your code for eigenval those lines of the function quer that deal with the case when m=n
and make the following adjustments:
(1) Increase the accuracy from p=7 to p=12 and output and display the number of iterations k
needed to archive this accuracy:
fprintf('the number of iterations is %i\n',k)
(2) Output the vector of the “sorted” diagonal entries of the matrix B, which is the last
iteration in the sequence, and assign it to the vector q. Also, convert to a zero the entries of the
vector q that are in the range of 10^(-7) from a zero. The outputs for this part should be:
2
q=sort(diag(B));
disp('The eigenvalues of A from QR factorization are:')
q=closetozeroroundoff(q,7)
Note: Function sort, applied to a real vector, re-arranges its entries in the ascending order.
**Second Method: you will calculate the eigenvalues of A using a built-in MATLAB function
eig().The outputs for this part will be:
e=eig(A);
e=sort(e,'ascend','ComparisonMethod','real');
disp('The eigenvalues of A from a MATLAB function are:')
e=closetozeroroundoff(e,7)
**Third Method: there is a built-in MATLAB function poly(A) that outputs the vector of the
coefficients of the characteristic polynomial of A written in the descending order according to
the degree. Output the vector
S=poly(A);
and, then, convert to zero the entries of S that are in the range of 10^(-7) from 0 – assign this
vector to S again (do not display S).
Next, output the characteristic polynomial of A in a symbolic form as indicated below:
disp('the MATLAB characteristic polynomial of A is:')
Q=poly2sym(S)
Continue working with the vector S: use the MATLAB function roots(S) to output the
vector of zeros of the polynomial whose coefficients form the vector S. “Sort” this vectors in
the ascending order with respect to the real part and assign it to c. Since c is the vector of
zeros of the characteristic polynomial, it is also the vector of the eigenvalues of A. Your
outputs for this part will be as below:
c=roots(S);
c=sort(c,'ascend','ComparisonMethod','real');
disp('The eigenvalues from the MATLAB characteristic polynomial are:')
c=closetozeroroundoff(c,7)
**Then, we will need to see how close the vectors q, e, and c are to each other.
Output and display the three 2-norms:
eq=norm(e-q)
cq=norm(c-q)
ec=norm(e-c)
**Next, we will use the vector of the eigenvalues e to build a polynomial whose zeros are the
entries of e, that is, to construct the characteristic polynomial of A. We will run a MATLAB
function poly(e) that returns the vector of the coefficients of the polynomial given the vector
e of its zeros. Your outputs here will be the characteristic polynomial R of the matrix A,
written in a symbolic form, and a message (see below):
disp('the constructed characteristic polynomial of A is:')
R=poly2sym(poly(e))
**Finally, you will chick if the symbolic polynomials Q and R are the same: run a logical
comparison statement and, if Q and R match, code a message:
disp('Yes! Q and R are the same!')
If they don’t, a message could be something like:
disp('What?! Q and R do not match?')
3
**Print the functions closetozeroroundoff, jord, and eigenval in your Live Script.
Note: function jord was created in Exercise 1 of Project 1.
**Run the function eigenval(A); as indicated below:
%(a)
A=[3 3;0 3]
[q,e,c,S]=eigenval(A);
%(b)
A=jord(5,4)
[q,e,c,S]=eigenval(A);
%(c)
A=ones(5);
[q,e,c,S]=eigenval(A);
%(d)
A=tril(magic(4))
[q,e,c,S]=eigenval(A);
%(e)
A=triu(magic(4),-1)
[q,e,c,S]=eigenval(A);
%(f)
A=[4 3 0 0; 7 10 3 0;0 3 1 5;0 0 6 7]
[q,e,c,S]=eigenval(A);
% Analyze the outputs eq, cq, ec after running the function eigenval on the matrices (a)-
(f) and write a comment whether the QR factorization gives a good approximation of the
eigenvalues of these matrices.
BONUS! (2 points)
Insert Section Break and run the part below as a separate Section.
**Input the matrix and run the function:
A=[4 0 0 0; 1 3 0 0; 0 -1 3 0; 0 0 5 4]
[q,e,c,S]=eigenval(A);
(It may take some time for the outputs to be displayed – please be patient!)
The last output message should indicate that the polynomials Q and R are not equal. You will
need to find out how close to each other are S and poly(e), which are the vectors of
coefficients of the polynomials Q and R, respectively. Proceed as follows:
** Initialize
p=1;
and set up a “while” loop in your Live Script using the variable p and the function
closetozeroroundoff(S-poly(e),p)
to output the largest value of p (p is a positive integer) for which the function
closetozeroroundoff outputs the zero vector.
Code a message that will display the number of the matching decimals in the vectors of the
coefficients of polynomials Q and R.
4
Part II. Eigenvectors & Diagonalization
Exercise 2 (5 points) Difficulty: Hard
In this exercise, you will work with the eigenvalues of an n n × matrix A. First, you will find
all eigenvalues of A. Then, you will consider the distinct eigenvalues and find orthonormal
bases for the corresponding eigenspaces and the dimensions of the eigenspaces. Next, you will
check if A is diagonalizable by applying the general theory. If a matrix A is diagonalizable,
the output has to contain an invertible matrix P and a diagonal matrix D, such that, 1 A PDP− = ,
or, equivalently, AP PD = and P is invertible. For diagonalizable matrices, you will run a
built-in MATLAB function, which also performs diagonalization, and you will compare its
outputs with the outputs P and D of your function.
**Create a function in MATLAB that begins with
function [L,P,D]=eigen(A)
format
[~,n]=size(A);
P=[];
D=[];
Your function [L,P,D]=eigen(A) will have a set of commands that generates the following
outputs for an n n × input matrix A.
Part 1. Output the vector L of all eigenvalues of A.
To do that, you will go through several steps. Please suppress all intermediate outputs and
display only the final output vector L.
**Output a row vector L of all eigenvalues, where each eigenvalue repeats as many times as
its multiplicity. A basic MATLAB command for this part
L=eig(A);
returns a column vector of all eigenvalues of A, and we use the function
L=transpose(L);
to get a row vector.
The input matrices have real eigenvalues; however, it is possible that the MATLAB command
eig() outputs them as complex numbers with negligibly small imaginary parts – to eliminate
this discrepancy, we will need to run the function
L=real(L);
on the vector L.
Then, we will “sort” the entries of L using the MATLAB command
L=sort(L);
To output vector L correctly, you will go through two more steps:
(1) In your code, you will need to ensure that the multiple eigenvalues show as the same
number. Use the function closetozeroroundoff with p = 7 to compare the eigenvalues and
assign the same value to the ones that are within 10^(-7) tolerance of each other. To perform
this task, you can go through the sequence of the sorted eigenvalues and, if there is a set of
eigenvalues that are within 10^(-7) from each other, you will assign the first eigenvalue that
comes across to all others in that set.
Important Note: here we work with the eigenvalues that are stored in MATLAB – we do not
round off any eigenvalues on this step.
5
(2) A singular matrix A has a zero eigenvalue; however, when running your code, a zero
eigenvalue may not show as a zero due to round-off errors in MATLAB.
If a matrix A is not invertible, you will need to run the function closetozeroroundoff with
p = 7 on the vector of eigenvalues L to ensure that the zero eigenvalues show as a zero, and
assign the output to L again. To check whether a matrix A is invertible or not, please use a
MATLAB command rank().
After performing all of the tasks outlined above, you will display your final output L with a
message:
fprintf('all eigenvalues of A are\n')
display(L)
Please make sure that L is a row vector of the “sorted” real eigenvalues of A where all
multiple eigenvalues are equal to each other and the zero eigenvalues (for singular matrices)
show as a zero.
Part 2. Construct an orthonormal basis W for each eigenspace.
**Output and display with a message a row vector M of all distinct eigenvalues (no repetitions
are allowed). The MATLAB command unique() can be used here.
For each distinct eigenvalue M(i), your function has to do the following:
**Find the multiplicity, m(i), of the eigenvalue M(i). Output it with the message:
fprintf('eigenvalue %d has multiplicity %i\n',M(i),m(i))
**Output an orthonormal basis W for the eigenspace for an eigenvalue M(i). Display it with
the message:
fprintf('a basis for the eigenspace for lambda=%d is:\n',M(i))
display(W)
Note: to output an orthonormal basis for the null space of a matrix, use a MATLAB command
null().
**Calculate the dimension d(i) of W. Output it with the message:
fprintf('dimension of eigenspace for lambda = %d is %i\n',M(i),d(i))
Part 3. Construct a Diagonalization when possible.
For help with this part, please review the material of Lecture 24.
**First, determine if A is diagonalizable: for each eigenvalue M(i), you will compare the
multiplicity, m(i), with the dimension of the corresponding eigenspace, d(i).
**If A is not diagonalizable, output a corresponding message and terminate the program – the
empty outputs for P and D will stay.
**If A is diagonalizable, output a corresponding message, and your function will continue
with constructing a diagonalization.
Output and display a matrix P constructed by combining the bases W for the eigenspaces.
Output and display the diagonal matrix D with the corresponding eigenvalues on its main
diagonal.
Verify that your diagonalization is correct, that is, both conditions hold: AP=PD and P is
invertible. To verify the first condition, you will need to use the function
closetozeroroundoff with p = 7. To verify the second condition, please use the command
rank().
6
If both conditions hold, output a message: 'Great! I got a diagonalization!'
Otherwise, output a message: 'Oops! I got a bug in my code!' and terminate the
program.
Part 4. Compare your outputs with the outputs of a built-in MATLAB function.
**If the diagonalization is confirmed, your code will continue with the following task.
There is a MATLAB built-in function that runs as
[U,V]=eig(A)
which, for a diagonalizable matrix A, generates a diagonal matrix V, with the eigenvalues of
A on its main diagonal, and an n n × invertible matrix U of the corresponding eigenvectors,
such that AU=UV. Place this function in your code and display the outputs U and V.
**Compare the output matrix U with the matrix P: write a set of commands that would check
on the following case:
P and U have either the same columns or the columns match up to a scalar (-1), where
the order of the columns does not matter.
If it is the case, output a message 'Columns of P and U are the same or match up to
scalar (-1)'.
If it is not the case (which is possible), the output message could be 'There is no specific
match among the columns of P and U'
Note: You will need to use the function closetozeroroundoff with p = 7 when comparing
the columns of the matrices P and U.
**Verify that the matrices D and V have the same diagonal elements (the order does not count
either). If it is the case, output a message 'The diagonal elements of D and V match'. If
it is not the case, output something like: 'That cannot be true!'
Hint: To perform this task, you can “sort” the diagonal elements of V and compare them with
the vector L using the function closetozeroroundoff with p = 7.
This is the end of the function eigen.
BONUS! (1 point)
Theory: A matrix A is said to be orthogonally diagonalizable if there exists an orthogonal
matrix P ( 1 T P P − = ) and a diagonal matrix D, such that, 1 A PDP− = , or equivalently,
T A PDP = .
Theorem: An n n × matrix A is orthogonally diagonalizable if and only if A is symmetric.
(For more information please read Section 7.1 of the textbook.)
Add a set of commands to your function eigen for a diagonalizable matrix A that, first, will
check if A is symmetric.
**If A is not symmetric, output a message:
disp('diagonalizable matrix A is not orthogonally diagonalizable')
**If A is symmetric, output a message
disp('matrix A is symmetric')
and verify the orthogonal diagonalization: check if the matrix P is orthogonal (you will need
to use here the function closetozeroroundoff with p = 7).
If the orthogonal diagonalization is confirmed, output a message:
7
disp('the orthogonal diagonalization is confirmed')
otherwise, output something like:
disp('Wow! A symmetric matrix is not orthogonally diagonalizable?!')
This is the end of your function eigen with the BONUS part if included.
**Print the functions eigen, closetozeroroundoff, jord in the Live Script.
Note: the function jord was created in Exercise 1 of Project 1.
**Run the function [L,P,D]=eigen(A); as indicated below:
%(a)
A=[3 3; 0 3]
[L,P,D]=eigen(A);
%(b)
A=[2 4 3;-4 -6 -3;3 3 1]
[L,P,D]=eigen(A);
%(c)
A=[4 0 1 0; 0 4 0 1; 1 0 4 0; 0 1 0 4]
[L,P,D]=eigen(A);
%(d)
A=jord(5,4);
[L,P,D]=eigen(A);
%(e)
A=ones(4)
[L,P,D]=eigen(A);
%(f)
A=[4 1 3 1;1 4 1 3;3 1 4 1;1 3 1 4]
[L,P,D]=eigen(A);
%(g)
A=[3 1 1;1 3 1;1 1 3]
[L,P,D]=eigen(A);
%(h)
A=magic(4)
[L,P,D]=eigen(A);
%(k)
A=magic(5)
[L,P,D]=eigen(A);
%(l)
A=hilb(5)
[L,P,D]=eigen(A);
Part III. Orthogonal Projections & Least-Squares Solutions
Exercise 3 (4 points) Difficulty: Moderate
In this exercise, you will create a function proj(A,b)which will work with the projection of a
vector b onto the Column Space of an m n × matrix A.
For a help with this exercise, please refer to Lectures 28-31.
8
Theory: When the columns of A are linearly independent, a vector p is the orthogonal
projection of a vector b∉ Col A onto the Col A if and only if p x = Aˆ , where xˆ is the unique
least-squares solution of Ax b = , or equivalently, the unique solution of the normal equations
T T AA A x b = .
Also, if p is the orthogonal projection of a vector b onto Col A, then there exists a unique
vector z, such that, z bp = − and z is orthogonal to Col A.
**Your program should allow a possibility that the columns of A are not linearly independent.
In order for the algorithm to work, we will use the function shrink() and generate a new
matrix, which will be assigned to A, whose columns form a basis for Col A. The function
shrink() was created in Exercise 2 of Project 3 and you should have it in your Current
Folder in MATLAB. If not, please see the code below:
function B=shrink(A)
[~,pivot]=rref(A);
B=A(:,pivot);
end
**Create your function which begins with the lines:
function [p,z]=proj(A,b)
format
A=shrink(A);
m=size(A,1);
p=[];
z=[];
The inputs are an m n × matrix A and a column vector b. Your outputs p and z will be either
empty (when the sizes of A and b do not match) or the projection of b onto the Col A and the
component of b orthogonal to the Col A, respectively.
**First, your function has to check whether the input vector b has exactly m entries, where m
is the number of rows in A. If it doesn’t, the program terminates with a message 'No
solution: sizes of A and b disagree', and the empty outputs for p and z will stay.
If b has exactly m entries, we proceed to the next step:
Determine whether it is the case that b∈ Col A – use here a MATLAB function rank().
**If b∈ Col A, output a message 'b is in Col A' and make assignments to the vectors p
and z based on the general theory (no computations are needed). Display the vectors p and z.
After that, the program terminates.
**If b∉ Col A, your function will continue with the following tasks:
Determine whether it is the case that b is orthogonal to Col A. Due to the round-off errors in
MATLAB computations, you will need to use the function closetozeroroundoff with p = 7
to check if this condition holds. If b is orthogonal to Col A, output a message 'b is
orthogonal to Col A' and make assignments to the vectors p and z based on the general
theory (no computations are needed either). Display the vectors p and z. After that, the
program terminates.
If b is not orthogonal to Col A, your function will continue:
9
**Find the least-squares solution x as the solution of the normal equations by using the
inverse matrix or the backslash operator, \ . Output it with a message:
fprintf('the least squares solution of inconsistent system Ax=b is')
display(x)
**Then, output a vector (do not display it)
x1=A\b;
and verify by using the function closetozeroroundoff with p=12 that the vectors x and x1
match within the given precision. If your code confirms it, display a message:
disp('A\b returns the least-squares solution of inconsistent system Ax=b')
**Next, output the projection p by using the least-squares solution x. Display p with the
message:
disp('the projection of b onto Col A is')
display(p)
**Then, output vector z, which is the component of b orthogonal to the Col A. Verify in your
code whether the vector z is, indeed, orthogonal to the Col A: use the function
closetozeroroundoff with p=7 here. If your code confirms it, display z with a message:
disp('the component of b orthogonal to Col A is')
diplay(z)
Otherwise, an output message could be:
disp('Oops! Is there a bug in my code?')
and the program terminates.
**Finally, use the vector z to compute the distance, d, from b to Col A. Output d with the
message:
fprintf('the distance from b to Col A is %i',d)
Hint: use here a MATLAB built-in function norm().
This is the end of your function proj.
**Print the functions closetozeroroundoff, shrink, proj in your Live Script.
**Run the function [p,z]=proj(A,b); as indicated below:
%(a)
A=magic(4), b=(1:4)'
[p,z]=proj(A,b);
%(b)
A=magic(6), b=A(:,6)
[p,z]=proj(A,b);
%(c)
A=magic(6), b=(1:5)'
[p,z]=proj(A,b);
%(d)
A=magic(5), b = rand(5,1)
[p,z]=proj(A,b);
%(e)
A=ones(4); A(:)=1:16, b=[1;0;1;0]
[p,z]=proj(A,b);
%(f)
B=ones(4); B(:)=1:16; A=null(B,'r'), b=ones(4,1)
[p,z]=proj(A,b);
10
% Analyze the outputs in part (d) and write a comment that would explain a reason why a
random vector b belongs to the Col A.
% For part (f), based on the input matrix A and the outputs, state to which vector space the
vector b has to belong.
**You will need to run some command in your Live Script to support your responses. Supply
them with the corresponding output messages and/or comments.
Exercise 4 (4 points) Difficulty: Moderate
In this exercise, we will work with the matrix equation Ax b = , where the columns of A are
linearly independent. We find unique “exact” solution for a consistent system (b∈Col A) and
unique least-squares solution for an inconsistent system (b∉Col A).
Please read the section Theory carefully. You may also find it helpful to review the material
of Lectures 28-31 prior to working on this exercise.
Theory: When solving in MATLAB a general linear system Ax b = with a matrix A whose
columns are linearly independent, we can always find the solution of a consistent system or
the least-squares solution of an inconsistent system by using the MATLAB backslash
operator, \ .
When matrix A has orthonormal columns, we can find the “exact” solution of a consistent
system using the Decomposition Theorem (Lecture 28) or the least squares-solution of an
inconsistent system using the Statement in Lecture 31 – both, the Theorem and the Statement,
are particular cases of the more general Orthogonal Decomposition Theorem (Lecture 29) and
both employ the same formulas for computing the entries of the solutions.
Recall: An m n × matrix A has an orthonormal set of columns if and only if A’*A=eye(n).
If a matrix A does not have orthonormal columns, we can, first, generate an orthonormal basis
U for Col A, then, find the projection b1 of the vector b onto the Col A, and, finally, calculate
the “exact” or the least-squares solution of a system Ax=b as the solution of a consistent
system Ax=b1. This method may be effective for calculating the least-squares solution, and,
for a purpose of an exercise, it can be also extended to a calculation of the “exact” solution,
which would unify an approach to solving a linear system in general.
Note: when b is in the Col A, that is, the system Ax=b is consistent, the projection of b onto
the Col A is the vector b itself, thus, the system Ax=b1 becomes Ax=b.
Recall: The projection b1 of b onto the Col A can be computed by the formula b1=U*U’*b,
where U is a matrix whose columns form an orthonormal basis for Col A.
Moreover, if x1 is the solution or the least-squares solution of the system Ax=b, that is,
A*x1=b1, we can compute the least-squares error, n1, of approximation of the vector b by
the elements of the Col A, where n1=norm(b-A*x1), or equivalently, n1=norm(b-b1).
Geometrically, n1 is the distance from b to the Col A (see Exercise 3 of this Project).
**Create a function in MATLAB that begins with
function [x1,x2]=solveall(A,b)
format compact
format long
[m,n]=size(A);
11
The inputs are an m n × matrix A, whose columns are linearly independent, and a vector m b∈ . The two outputs x1 and x2 are either two solutions of a consistent system or two
least-squares solutions of an inconsistent system calculated in two different ways.
**Continue your function with a conditional statement that has to determine if b∈Col A or
not, that is, whether the system is consistent or not – use the MATLAB function rank() in
this part. If the system is consistent, output a message:
disp('the system is consistent - find the "exact" solution')
otherwise, output a message:
disp('the system is inconsistent - find the least-squares solution')
**Next, calculate the solution or the least-squares solution of the system Ax=b, vector x1, by
using the backslash operator, \ . Output and display it with a message as indicated below:
x1=A\b;
disp('the solution calculated by the backslash operator is')
display(x1)
Notice that the type of the solution has been stated in the previous message.
**Then, determine whether the matrix A has orthonormal columns – you will need to use here
the function closetozeroroundoff() with p = 7.
If matrix A has orthonormal columns, output a message:
disp('A has orthonormal columns')
and find the solution (or the least-squares solution) x2 using the Orthogonal Decomposition
Theorem (Lecture 29), which is more general form of the Decomposition Theorem (Lecture
28) and the Statement in Lecture 31 (see Theory above).
Hint: You can either use a “for loop” to calculate the entries of the vector x2, one by one, or
you can employ the transpose of A to output x2 at once by a MATLAB operation (second way
is preferable – remember, we work with a matrix whose columns form an orthonormal set).
Display the solution x2 with a message:
disp('solution calculated by the Orthogonal Decomposition Theorem is')
display(x2)
**If matrix A does not have orthonormal columns, output a message:
disp('A does not have orthonormal columns')
and find the solution (or the least-squares solution) of Ax b = by, first, calculating the
projection of the vector b onto the Col A. Proceed as follows:
**Run the MATLAB function orth() to output and display a matrix U whose columns form
an orthonormal basis for the Col A supplied with a message:
disp('an orthonormal basis for Col A is')
U=orth(A)
**Next, using the created orthonormal basis U for Col A, calculate the projection b1 of the
vector b onto the Col A (see Theory above). Display vector b1 with a message:
disp('the projection of b onto Col A is')
display(b1)
**Then, output the solution (or the least-squares solution) x2 as the “exact” solution of the
equation Ax = b1 (use the backslash operator). Display x2 with the message:
disp('the solution calculated using the projection b1 is')
display(x2)
12
**Verify that the solution x2 is calculated correctly: use a conditional statement and the
function closetozeroroundoff()with p = 12 that you will run on the vector x1-x2. If your
code confirms that the solutions calculated by two different methods match, code a message:
disp('solutions x1 and x2 are sufficiently close to each other')
Otherwise, your message could be
disp('Check the code!')
and the program terminates.
**If the code passes the check point above, your function will continue with a few more tasks:
Use the solution x1 to output the least-squares error n1 of approximation of b by the
elements of the Col A (see Theory above). Display n1 with a message:
disp('least-squares error of approximation of b by elements of Col A is')
display(n1)
**Next, input a vector
x=rand(n,1);
and compute the error of approximation, n2, of the vector b by a vector A*x of the Col A (for
a random vector x). Output and display it as below:
n2=norm(b-A*x);
disp('an error of approximation of b by A*x of Col A for a random x is')
display(n2)
This is the end of your function solveall.
**Print the functions closetozeroroundoff, shrink, solveall in your Live Script.
Note: the function shrink was created in Project 2. We use it here to input some matrices. If
you do not have it in your Current Folder, please see the code in Exercise3 of this Project.
**Run the function [x1,x2]=solveall(A,b); as indicated below:
%(a)
A=magic(4); b=A(:,4), A=orth(A)
[x1,x2]=solveall(A,b);
%(b)
A=magic(5), b=rand(5,1)
[x1,x2]=solveall(A,b);
%(c)
A=magic(4); A=shrink(A), b=ones(4,1)
[x1,x2]=solveall(A,b);
%(d)
A=magic(4); A=shrink(A), b=rand(4,1)
[x1,x2]=solveall(A,b);
%(e)
A=magic(6); A=orth(A), b=rand(6,1)
[x1,x2]=solveall(A,b);
% Analyze the outputs n1 and n2 for both consistent and inconsistent systems. Based on the
geometrical meaning of the least-squares error n1 explain the difference between the forms of
the outputs n1 and n2 for consistent systems vs. inconsistent systems.
% Compare n1 and n2 for each inconsistent system and write a comment whether the leastsquares
solution, x1 (or x2) may, indeed, minimize the distance between the vector b and the
vectors Ax of the Col A (which has to be true according to the general theory).
13
Part IV. Orthogonalization by Gram-Schmidt with Normalization
According to the theory, we can create an orthonormal basis for the column space of a matrix
by using the Gram-Schmidt process. However, a direct computer realization of the GramSchmidt
process involves round-off errors that build up when the vectors are calculated, one
by one. In this exercise, we will be programming a Modified Gram-Schmidt process, which
involves normalization, and we may need to run the algorithm several times until the required
accuracy is archived.
Exercise 5 (5 points) Difficulty: Hard
**Create a function that begins with the following lines:
function U=orthonorm(A)
format compact
format long
The input is an m x n matrix A. The final output U will be a matrix whose columns form an
orthonormal basis for Col A satisfying the required accuracy.
We begin with generating a basis X for the Col A, which will be represented by the matrix X.
In general, the columns of A span Col A, but they may not be linearly independent, thus, it is
possible that not all columns of A will be in a basis X for Col A.
**If the set of columns of A is linearly independent, we assign to X matrix A itself and
display a basis X with a message:
disp('a basis X for Col A is formed by all columns of A')
X=A
**Otherwise, we “shrink” the set of columns of A into a basis for Col A and assign it to X. In
this case, we display both a basis X and the indexes of the pivot columns of A as below:
disp('a basis X for Col A is formed by the columns of A with indexes')
X=shrink_1(A)
where the function shrink_1 is a modification of the function shrink constructed in the way
that, along with a basis for the Col A, it also displays the indexes of the pivot columns of A
(see Exercise 2 of Project 3). You will need to create and save the function shrink_1 in your
Current Folder in MATLAB.
**Next, assign (do not display):
U=X;
n=size(U,2);
and continue with a check if the columns of U (or X) already form an orthonormal basis for
Col A with the accuracy p=15, that is, if the function
closetozeroroundoff(U'*U-eye(n),15)
returns the zero matrix.
If it is the case, display U with a message as below:
disp('a required orthonormal basis U for Col A is X')
display(U)
and your code terminates here.
If the columns of U do not form the required orthonormal basis for Col A, we proceed with
orthogonalization by Gram-Schmidt along with normalization. It is possible that running this
process just one time will not create an orthonormal basis U with the required accuracy
14
(p=15). In this case, we will need to re-run the process on the created matrix U and repeat it as
many times as needed until the required accuracy is archived. We proceed as follows.
**Initialize
count=1;
and set up an exterior “while” loop using the variable count that will keep track on how many
times we need to run the process of orthonormalization.
**Within the exterior loop, set up a “for” loop and use the formulas for Gram-Schmidt
orthogonalization (Lecture 30) to recalculate the columns of the matrix U, one by one, making
the following adjustments: after completing a consecutive column, you will need to normalize
it (use a MATLAB function norm()).
Hint: here we need to employ MATLAB operations such as matrix-vector multiplic
Please read the Instructions located on the Assignments page prior
to working on the Project.
BEGIN with creating a Live Script Project4.
Note: All exercises in this project have to be completed in the Live Script using the Live
Editor. Please refer to the MATLAB video that explains how to use the Live Script:
https://www.mathworks.com/videos/using-the-live-editor-117940.html?s_tid=srchtitle
The final script has to be generated by exporting the Live Script to PDF.
Each exercise has to begin with the line
Exercise#
You should also mark down the parts such as (a), (b), (c), and etc. This makes grading easier.
Important: we use the default format short for the numbers in all exercises unless it is
specified otherwise. We do not employ format rat since it may cause problems with
running the codes and displaying matrices in the Live Script. If format long had been used
in an exercise, please make sure to return to the default format in the next exercise.
Part I. Evaluating Eigenvalues
Exercise 1 (4 points): Difficulty: Easy
In this exercise, you will calculate the eigenvalues of a matrix by using various methods and
functions – one of the functions, quer, was created earlier in Exercise 4 of Project 3, and the
other functions will be MATLAB built-in functions.
NOTE: The function quer from Project 3 can be used for approximating eigenvalues by
consecutive iterations, but not for all matrices the process converges. For computing
eigenvalues in MATLAB, especially complex eigenvalues, it is advisable to use a MATLAB
built-in function eig().
**Create a function that begins with
function [q,e,c,S]=eigenval(A)
format
The input is a square matrix A. The outputs will be specified below when you proceed with
the programming.
**First Method: your function will calculate the eigenvalues by the QR factorization. Include
in your code for eigenval those lines of the function quer that deal with the case when m=n
and make the following adjustments:
(1) Increase the accuracy from p=7 to p=12 and output and display the number of iterations k
needed to archive this accuracy:
fprintf('the number of iterations is %i\n',k)
(2) Output the vector of the “sorted” diagonal entries of the matrix B, which is the last
iteration in the sequence, and assign it to the vector q. Also, convert to a zero the entries of the
vector q that are in the range of 10^(-7) from a zero. The outputs for this part should be:
2
q=sort(diag(B));
disp('The eigenvalues of A from QR factorization are:')
q=closetozeroroundoff(q,7)
Note: Function sort, applied to a real vector, re-arranges its entries in the ascending order.
**Second Method: you will calculate the eigenvalues of A using a built-in MATLAB function
eig().The outputs for this part will be:
e=eig(A);
e=sort(e,'ascend','ComparisonMethod','real');
disp('The eigenvalues of A from a MATLAB function are:')
e=closetozeroroundoff(e,7)
**Third Method: there is a built-in MATLAB function poly(A) that outputs the vector of the
coefficients of the characteristic polynomial of A written in the descending order according to
the degree. Output the vector
S=poly(A);
and, then, convert to zero the entries of S that are in the range of 10^(-7) from 0 – assign this
vector to S again (do not display S).
Next, output the characteristic polynomial of A in a symbolic form as indicated below:
disp('the MATLAB characteristic polynomial of A is:')
Q=poly2sym(S)
Continue working with the vector S: use the MATLAB function roots(S) to output the
vector of zeros of the polynomial whose coefficients form the vector S. “Sort” this vectors in
the ascending order with respect to the real part and assign it to c. Since c is the vector of
zeros of the characteristic polynomial, it is also the vector of the eigenvalues of A. Your
outputs for this part will be as below:
c=roots(S);
c=sort(c,'ascend','ComparisonMethod','real');
disp('The eigenvalues from the MATLAB characteristic polynomial are:')
c=closetozeroroundoff(c,7)
**Then, we will need to see how close the vectors q, e, and c are to each other.
Output and display the three 2-norms:
eq=norm(e-q)
cq=norm(c-q)
ec=norm(e-c)
**Next, we will use the vector of the eigenvalues e to build a polynomial whose zeros are the
entries of e, that is, to construct the characteristic polynomial of A. We will run a MATLAB
function poly(e) that returns the vector of the coefficients of the polynomial given the vector
e of its zeros. Your outputs here will be the characteristic polynomial R of the matrix A,
written in a symbolic form, and a message (see below):
disp('the constructed characteristic polynomial of A is:')
R=poly2sym(poly(e))
**Finally, you will chick if the symbolic polynomials Q and R are the same: run a logical
comparison statement and, if Q and R match, code a message:
disp('Yes! Q and R are the same!')
If they don’t, a message could be something like:
disp('What?! Q and R do not match?')
3
**Print the functions closetozeroroundoff, jord, and eigenval in your Live Script.
Note: function jord was created in Exercise 1 of Project 1.
**Run the function eigenval(A); as indicated below:
%(a)
A=[3 3;0 3]
[q,e,c,S]=eigenval(A);
%(b)
A=jord(5,4)
[q,e,c,S]=eigenval(A);
%(c)
A=ones(5);
[q,e,c,S]=eigenval(A);
%(d)
A=tril(magic(4))
[q,e,c,S]=eigenval(A);
%(e)
A=triu(magic(4),-1)
[q,e,c,S]=eigenval(A);
%(f)
A=[4 3 0 0; 7 10 3 0;0 3 1 5;0 0 6 7]
[q,e,c,S]=eigenval(A);
% Analyze the outputs eq, cq, ec after running the function eigenval on the matrices (a)-
(f) and write a comment whether the QR factorization gives a good approximation of the
eigenvalues of these matrices.
BONUS! (2 points)
Insert Section Break and run the part below as a separate Section.
**Input the matrix and run the function:
A=[4 0 0 0; 1 3 0 0; 0 -1 3 0; 0 0 5 4]
[q,e,c,S]=eigenval(A);
(It may take some time for the outputs to be displayed – please be patient!)
The last output message should indicate that the polynomials Q and R are not equal. You will
need to find out how close to each other are S and poly(e), which are the vectors of
coefficients of the polynomials Q and R, respectively. Proceed as follows:
** Initialize
p=1;
and set up a “while” loop in your Live Script using the variable p and the function
closetozeroroundoff(S-poly(e),p)
to output the largest value of p (p is a positive integer) for which the function
closetozeroroundoff outputs the zero vector.
Code a message that will display the number of the matching decimals in the vectors of the
coefficients of polynomials Q and R.
4
Part II. Eigenvectors & Diagonalization
Exercise 2 (5 points) Difficulty: Hard
In this exercise, you will work with the eigenvalues of an n n × matrix A. First, you will find
all eigenvalues of A. Then, you will consider the distinct eigenvalues and find orthonormal
bases for the corresponding eigenspaces and the dimensions of the eigenspaces. Next, you will
check if A is diagonalizable by applying the general theory. If a matrix A is diagonalizable,
the output has to contain an invertible matrix P and a diagonal matrix D, such that, 1 A PDP− = ,
or, equivalently, AP PD = and P is invertible. For diagonalizable matrices, you will run a
built-in MATLAB function, which also performs diagonalization, and you will compare its
outputs with the outputs P and D of your function.
**Create a function in MATLAB that begins with
function [L,P,D]=eigen(A)
format
[~,n]=size(A);
P=[];
D=[];
Your function [L,P,D]=eigen(A) will have a set of commands that generates the following
outputs for an n n × input matrix A.
Part 1. Output the vector L of all eigenvalues of A.
To do that, you will go through several steps. Please suppress all intermediate outputs and
display only the final output vector L.
**Output a row vector L of all eigenvalues, where each eigenvalue repeats as many times as
its multiplicity. A basic MATLAB command for this part
L=eig(A);
returns a column vector of all eigenvalues of A, and we use the function
L=transpose(L);
to get a row vector.
The input matrices have real eigenvalues; however, it is possible that the MATLAB command
eig() outputs them as complex numbers with negligibly small imaginary parts – to eliminate
this discrepancy, we will need to run the function
L=real(L);
on the vector L.
Then, we will “sort” the entries of L using the MATLAB command
L=sort(L);
To output vector L correctly, you will go through two more steps:
(1) In your code, you will need to ensure that the multiple eigenvalues show as the same
number. Use the function closetozeroroundoff with p = 7 to compare the eigenvalues and
assign the same value to the ones that are within 10^(-7) tolerance of each other. To perform
this task, you can go through the sequence of the sorted eigenvalues and, if there is a set of
eigenvalues that are within 10^(-7) from each other, you will assign the first eigenvalue that
comes across to all others in that set.
Important Note: here we work with the eigenvalues that are stored in MATLAB – we do not
round off any eigenvalues on this step.
5
(2) A singular matrix A has a zero eigenvalue; however, when running your code, a zero
eigenvalue may not show as a zero due to round-off errors in MATLAB.
If a matrix A is not invertible, you will need to run the function closetozeroroundoff with
p = 7 on the vector of eigenvalues L to ensure that the zero eigenvalues show as a zero, and
assign the output to L again. To check whether a matrix A is invertible or not, please use a
MATLAB command rank().
After performing all of the tasks outlined above, you will display your final output L with a
message:
fprintf('all eigenvalues of A are\n')
display(L)
Please make sure that L is a row vector of the “sorted” real eigenvalues of A where all
multiple eigenvalues are equal to each other and the zero eigenvalues (for singular matrices)
show as a zero.
Part 2. Construct an orthonormal basis W for each eigenspace.
**Output and display with a message a row vector M of all distinct eigenvalues (no repetitions
are allowed). The MATLAB command unique() can be used here.
For each distinct eigenvalue M(i), your function has to do the following:
**Find the multiplicity, m(i), of the eigenvalue M(i). Output it with the message:
fprintf('eigenvalue %d has multiplicity %i\n',M(i),m(i))
**Output an orthonormal basis W for the eigenspace for an eigenvalue M(i). Display it with
the message:
fprintf('a basis for the eigenspace for lambda=%d is:\n',M(i))
display(W)
Note: to output an orthonormal basis for the null space of a matrix, use a MATLAB command
null().
**Calculate the dimension d(i) of W. Output it with the message:
fprintf('dimension of eigenspace for lambda = %d is %i\n',M(i),d(i))
Part 3. Construct a Diagonalization when possible.
For help with this part, please review the material of Lecture 24.
**First, determine if A is diagonalizable: for each eigenvalue M(i), you will compare the
multiplicity, m(i), with the dimension of the corresponding eigenspace, d(i).
**If A is not diagonalizable, output a corresponding message and terminate the program – the
empty outputs for P and D will stay.
**If A is diagonalizable, output a corresponding message, and your function will continue
with constructing a diagonalization.
Output and display a matrix P constructed by combining the bases W for the eigenspaces.
Output and display the diagonal matrix D with the corresponding eigenvalues on its main
diagonal.
Verify that your diagonalization is correct, that is, both conditions hold: AP=PD and P is
invertible. To verify the first condition, you will need to use the function
closetozeroroundoff with p = 7. To verify the second condition, please use the command
rank().
6
If both conditions hold, output a message: 'Great! I got a diagonalization!'
Otherwise, output a message: 'Oops! I got a bug in my code!' and terminate the
program.
Part 4. Compare your outputs with the outputs of a built-in MATLAB function.
**If the diagonalization is confirmed, your code will continue with the following task.
There is a MATLAB built-in function that runs as
[U,V]=eig(A)
which, for a diagonalizable matrix A, generates a diagonal matrix V, with the eigenvalues of
A on its main diagonal, and an n n × invertible matrix U of the corresponding eigenvectors,
such that AU=UV. Place this function in your code and display the outputs U and V.
**Compare the output matrix U with the matrix P: write a set of commands that would check
on the following case:
P and U have either the same columns or the columns match up to a scalar (-1), where
the order of the columns does not matter.
If it is the case, output a message 'Columns of P and U are the same or match up to
scalar (-1)'.
If it is not the case (which is possible), the output message could be 'There is no specific
match among the columns of P and U'
Note: You will need to use the function closetozeroroundoff with p = 7 when comparing
the columns of the matrices P and U.
**Verify that the matrices D and V have the same diagonal elements (the order does not count
either). If it is the case, output a message 'The diagonal elements of D and V match'. If
it is not the case, output something like: 'That cannot be true!'
Hint: To perform this task, you can “sort” the diagonal elements of V and compare them with
the vector L using the function closetozeroroundoff with p = 7.
This is the end of the function eigen.
BONUS! (1 point)
Theory: A matrix A is said to be orthogonally diagonalizable if there exists an orthogonal
matrix P ( 1 T P P − = ) and a diagonal matrix D, such that, 1 A PDP− = , or equivalently,
T A PDP = .
Theorem: An n n × matrix A is orthogonally diagonalizable if and only if A is symmetric.
(For more information please read Section 7.1 of the textbook.)
Add a set of commands to your function eigen for a diagonalizable matrix A that, first, will
check if A is symmetric.
**If A is not symmetric, output a message:
disp('diagonalizable matrix A is not orthogonally diagonalizable')
**If A is symmetric, output a message
disp('matrix A is symmetric')
and verify the orthogonal diagonalization: check if the matrix P is orthogonal (you will need
to use here the function closetozeroroundoff with p = 7).
If the orthogonal diagonalization is confirmed, output a message:
7
disp('the orthogonal diagonalization is confirmed')
otherwise, output something like:
disp('Wow! A symmetric matrix is not orthogonally diagonalizable?!')
This is the end of your function eigen with the BONUS part if included.
**Print the functions eigen, closetozeroroundoff, jord in the Live Script.
Note: the function jord was created in Exercise 1 of Project 1.
**Run the function [L,P,D]=eigen(A); as indicated below:
%(a)
A=[3 3; 0 3]
[L,P,D]=eigen(A);
%(b)
A=[2 4 3;-4 -6 -3;3 3 1]
[L,P,D]=eigen(A);
%(c)
A=[4 0 1 0; 0 4 0 1; 1 0 4 0; 0 1 0 4]
[L,P,D]=eigen(A);
%(d)
A=jord(5,4);
[L,P,D]=eigen(A);
%(e)
A=ones(4)
[L,P,D]=eigen(A);
%(f)
A=[4 1 3 1;1 4 1 3;3 1 4 1;1 3 1 4]
[L,P,D]=eigen(A);
%(g)
A=[3 1 1;1 3 1;1 1 3]
[L,P,D]=eigen(A);
%(h)
A=magic(4)
[L,P,D]=eigen(A);
%(k)
A=magic(5)
[L,P,D]=eigen(A);
%(l)
A=hilb(5)
[L,P,D]=eigen(A);
Part III. Orthogonal Projections & Least-Squares Solutions
Exercise 3 (4 points) Difficulty: Moderate
In this exercise, you will create a function proj(A,b)which will work with the projection of a
vector b onto the Column Space of an m n × matrix A.
For a help with this exercise, please refer to Lectures 28-31.
8
Theory: When the columns of A are linearly independent, a vector p is the orthogonal
projection of a vector b∉ Col A onto the Col A if and only if p x = Aˆ , where xˆ is the unique
least-squares solution of Ax b = , or equivalently, the unique solution of the normal equations
T T AA A x b = .
Also, if p is the orthogonal projection of a vector b onto Col A, then there exists a unique
vector z, such that, z bp = − and z is orthogonal to Col A.
**Your program should allow a possibility that the columns of A are not linearly independent.
In order for the algorithm to work, we will use the function shrink() and generate a new
matrix, which will be assigned to A, whose columns form a basis for Col A. The function
shrink() was created in Exercise 2 of Project 3 and you should have it in your Current
Folder in MATLAB. If not, please see the code below:
function B=shrink(A)
[~,pivot]=rref(A);
B=A(:,pivot);
end
**Create your function which begins with the lines:
function [p,z]=proj(A,b)
format
A=shrink(A);
m=size(A,1);
p=[];
z=[];
The inputs are an m n × matrix A and a column vector b. Your outputs p and z will be either
empty (when the sizes of A and b do not match) or the projection of b onto the Col A and the
component of b orthogonal to the Col A, respectively.
**First, your function has to check whether the input vector b has exactly m entries, where m
is the number of rows in A. If it doesn’t, the program terminates with a message 'No
solution: sizes of A and b disagree', and the empty outputs for p and z will stay.
If b has exactly m entries, we proceed to the next step:
Determine whether it is the case that b∈ Col A – use here a MATLAB function rank().
**If b∈ Col A, output a message 'b is in Col A' and make assignments to the vectors p
and z based on the general theory (no computations are needed). Display the vectors p and z.
After that, the program terminates.
**If b∉ Col A, your function will continue with the following tasks:
Determine whether it is the case that b is orthogonal to Col A. Due to the round-off errors in
MATLAB computations, you will need to use the function closetozeroroundoff with p = 7
to check if this condition holds. If b is orthogonal to Col A, output a message 'b is
orthogonal to Col A' and make assignments to the vectors p and z based on the general
theory (no computations are needed either). Display the vectors p and z. After that, the
program terminates.
If b is not orthogonal to Col A, your function will continue:
9
**Find the least-squares solution x as the solution of the normal equations by using the
inverse matrix or the backslash operator, \ . Output it with a message:
fprintf('the least squares solution of inconsistent system Ax=b is')
display(x)
**Then, output a vector (do not display it)
x1=A\b;
and verify by using the function closetozeroroundoff with p=12 that the vectors x and x1
match within the given precision. If your code confirms it, display a message:
disp('A\b returns the least-squares solution of inconsistent system Ax=b')
**Next, output the projection p by using the least-squares solution x. Display p with the
message:
disp('the projection of b onto Col A is')
display(p)
**Then, output vector z, which is the component of b orthogonal to the Col A. Verify in your
code whether the vector z is, indeed, orthogonal to the Col A: use the function
closetozeroroundoff with p=7 here. If your code confirms it, display z with a message:
disp('the component of b orthogonal to Col A is')
diplay(z)
Otherwise, an output message could be:
disp('Oops! Is there a bug in my code?')
and the program terminates.
**Finally, use the vector z to compute the distance, d, from b to Col A. Output d with the
message:
fprintf('the distance from b to Col A is %i',d)
Hint: use here a MATLAB built-in function norm().
This is the end of your function proj.
**Print the functions closetozeroroundoff, shrink, proj in your Live Script.
**Run the function [p,z]=proj(A,b); as indicated below:
%(a)
A=magic(4), b=(1:4)'
[p,z]=proj(A,b);
%(b)
A=magic(6), b=A(:,6)
[p,z]=proj(A,b);
%(c)
A=magic(6), b=(1:5)'
[p,z]=proj(A,b);
%(d)
A=magic(5), b = rand(5,1)
[p,z]=proj(A,b);
%(e)
A=ones(4); A(:)=1:16, b=[1;0;1;0]
[p,z]=proj(A,b);
%(f)
B=ones(4); B(:)=1:16; A=null(B,'r'), b=ones(4,1)
[p,z]=proj(A,b);
10
% Analyze the outputs in part (d) and write a comment that would explain a reason why a
random vector b belongs to the Col A.
% For part (f), based on the input matrix A and the outputs, state to which vector space the
vector b has to belong.
**You will need to run some command in your Live Script to support your responses. Supply
them with the corresponding output messages and/or comments.
Exercise 4 (4 points) Difficulty: Moderate
In this exercise, we will work with the matrix equation Ax b = , where the columns of A are
linearly independent. We find unique “exact” solution for a consistent system (b∈Col A) and
unique least-squares solution for an inconsistent system (b∉Col A).
Please read the section Theory carefully. You may also find it helpful to review the material
of Lectures 28-31 prior to working on this exercise.
Theory: When solving in MATLAB a general linear system Ax b = with a matrix A whose
columns are linearly independent, we can always find the solution of a consistent system or
the least-squares solution of an inconsistent system by using the MATLAB backslash
operator, \ .
When matrix A has orthonormal columns, we can find the “exact” solution of a consistent
system using the Decomposition Theorem (Lecture 28) or the least squares-solution of an
inconsistent system using the Statement in Lecture 31 – both, the Theorem and the Statement,
are particular cases of the more general Orthogonal Decomposition Theorem (Lecture 29) and
both employ the same formulas for computing the entries of the solutions.
Recall: An m n × matrix A has an orthonormal set of columns if and only if A’*A=eye(n).
If a matrix A does not have orthonormal columns, we can, first, generate an orthonormal basis
U for Col A, then, find the projection b1 of the vector b onto the Col A, and, finally, calculate
the “exact” or the least-squares solution of a system Ax=b as the solution of a consistent
system Ax=b1. This method may be effective for calculating the least-squares solution, and,
for a purpose of an exercise, it can be also extended to a calculation of the “exact” solution,
which would unify an approach to solving a linear system in general.
Note: when b is in the Col A, that is, the system Ax=b is consistent, the projection of b onto
the Col A is the vector b itself, thus, the system Ax=b1 becomes Ax=b.
Recall: The projection b1 of b onto the Col A can be computed by the formula b1=U*U’*b,
where U is a matrix whose columns form an orthonormal basis for Col A.
Moreover, if x1 is the solution or the least-squares solution of the system Ax=b, that is,
A*x1=b1, we can compute the least-squares error, n1, of approximation of the vector b by
the elements of the Col A, where n1=norm(b-A*x1), or equivalently, n1=norm(b-b1).
Geometrically, n1 is the distance from b to the Col A (see Exercise 3 of this Project).
**Create a function in MATLAB that begins with
function [x1,x2]=solveall(A,b)
format compact
format long
[m,n]=size(A);
11
The inputs are an m n × matrix A, whose columns are linearly independent, and a vector m b∈ . The two outputs x1 and x2 are either two solutions of a consistent system or two
least-squares solutions of an inconsistent system calculated in two different ways.
**Continue your function with a conditional statement that has to determine if b∈Col A or
not, that is, whether the system is consistent or not – use the MATLAB function rank() in
this part. If the system is consistent, output a message:
disp('the system is consistent - find the "exact" solution')
otherwise, output a message:
disp('the system is inconsistent - find the least-squares solution')
**Next, calculate the solution or the least-squares solution of the system Ax=b, vector x1, by
using the backslash operator, \ . Output and display it with a message as indicated below:
x1=A\b;
disp('the solution calculated by the backslash operator is')
display(x1)
Notice that the type of the solution has been stated in the previous message.
**Then, determine whether the matrix A has orthonormal columns – you will need to use here
the function closetozeroroundoff() with p = 7.
If matrix A has orthonormal columns, output a message:
disp('A has orthonormal columns')
and find the solution (or the least-squares solution) x2 using the Orthogonal Decomposition
Theorem (Lecture 29), which is more general form of the Decomposition Theorem (Lecture
28) and the Statement in Lecture 31 (see Theory above).
Hint: You can either use a “for loop” to calculate the entries of the vector x2, one by one, or
you can employ the transpose of A to output x2 at once by a MATLAB operation (second way
is preferable – remember, we work with a matrix whose columns form an orthonormal set).
Display the solution x2 with a message:
disp('solution calculated by the Orthogonal Decomposition Theorem is')
display(x2)
**If matrix A does not have orthonormal columns, output a message:
disp('A does not have orthonormal columns')
and find the solution (or the least-squares solution) of Ax b = by, first, calculating the
projection of the vector b onto the Col A. Proceed as follows:
**Run the MATLAB function orth() to output and display a matrix U whose columns form
an orthonormal basis for the Col A supplied with a message:
disp('an orthonormal basis for Col A is')
U=orth(A)
**Next, using the created orthonormal basis U for Col A, calculate the projection b1 of the
vector b onto the Col A (see Theory above). Display vector b1 with a message:
disp('the projection of b onto Col A is')
display(b1)
**Then, output the solution (or the least-squares solution) x2 as the “exact” solution of the
equation Ax = b1 (use the backslash operator). Display x2 with the message:
disp('the solution calculated using the projection b1 is')
display(x2)
12
**Verify that the solution x2 is calculated correctly: use a conditional statement and the
function closetozeroroundoff()with p = 12 that you will run on the vector x1-x2. If your
code confirms that the solutions calculated by two different methods match, code a message:
disp('solutions x1 and x2 are sufficiently close to each other')
Otherwise, your message could be
disp('Check the code!')
and the program terminates.
**If the code passes the check point above, your function will continue with a few more tasks:
Use the solution x1 to output the least-squares error n1 of approximation of b by the
elements of the Col A (see Theory above). Display n1 with a message:
disp('least-squares error of approximation of b by elements of Col A is')
display(n1)
**Next, input a vector
x=rand(n,1);
and compute the error of approximation, n2, of the vector b by a vector A*x of the Col A (for
a random vector x). Output and display it as below:
n2=norm(b-A*x);
disp('an error of approximation of b by A*x of Col A for a random x is')
display(n2)
This is the end of your function solveall.
**Print the functions closetozeroroundoff, shrink, solveall in your Live Script.
Note: the function shrink was created in Project 2. We use it here to input some matrices. If
you do not have it in your Current Folder, please see the code in Exercise3 of this Project.
**Run the function [x1,x2]=solveall(A,b); as indicated below:
%(a)
A=magic(4); b=A(:,4), A=orth(A)
[x1,x2]=solveall(A,b);
%(b)
A=magic(5), b=rand(5,1)
[x1,x2]=solveall(A,b);
%(c)
A=magic(4); A=shrink(A), b=ones(4,1)
[x1,x2]=solveall(A,b);
%(d)
A=magic(4); A=shrink(A), b=rand(4,1)
[x1,x2]=solveall(A,b);
%(e)
A=magic(6); A=orth(A), b=rand(6,1)
[x1,x2]=solveall(A,b);
% Analyze the outputs n1 and n2 for both consistent and inconsistent systems. Based on the
geometrical meaning of the least-squares error n1 explain the difference between the forms of
the outputs n1 and n2 for consistent systems vs. inconsistent systems.
% Compare n1 and n2 for each inconsistent system and write a comment whether the leastsquares
solution, x1 (or x2) may, indeed, minimize the distance between the vector b and the
vectors Ax of the Col A (which has to be true according to the general theory).
13
Part IV. Orthogonalization by Gram-Schmidt with Normalization
According to the theory, we can create an orthonormal basis for the column space of a matrix
by using the Gram-Schmidt process. However, a direct computer realization of the GramSchmidt
process involves round-off errors that build up when the vectors are calculated, one
by one. In this exercise, we will be programming a Modified Gram-Schmidt process, which
involves normalization, and we may need to run the algorithm several times until the required
accuracy is archived.
Exercise 5 (5 points) Difficulty: Hard
**Create a function that begins with the following lines:
function U=orthonorm(A)
format compact
format long
The input is an m x n matrix A. The final output U will be a matrix whose columns form an
orthonormal basis for Col A satisfying the required accuracy.
We begin with generating a basis X for the Col A, which will be represented by the matrix X.
In general, the columns of A span Col A, but they may not be linearly independent, thus, it is
possible that not all columns of A will be in a basis X for Col A.
**If the set of columns of A is linearly independent, we assign to X matrix A itself and
display a basis X with a message:
disp('a basis X for Col A is formed by all columns of A')
X=A
**Otherwise, we “shrink” the set of columns of A into a basis for Col A and assign it to X. In
this case, we display both a basis X and the indexes of the pivot columns of A as below:
disp('a basis X for Col A is formed by the columns of A with indexes')
X=shrink_1(A)
where the function shrink_1 is a modification of the function shrink constructed in the way
that, along with a basis for the Col A, it also displays the indexes of the pivot columns of A
(see Exercise 2 of Project 3). You will need to create and save the function shrink_1 in your
Current Folder in MATLAB.
**Next, assign (do not display):
U=X;
n=size(U,2);
and continue with a check if the columns of U (or X) already form an orthonormal basis for
Col A with the accuracy p=15, that is, if the function
closetozeroroundoff(U'*U-eye(n),15)
returns the zero matrix.
If it is the case, display U with a message as below:
disp('a required orthonormal basis U for Col A is X')
display(U)
and your code terminates here.
If the columns of U do not form the required orthonormal basis for Col A, we proceed with
orthogonalization by Gram-Schmidt along with normalization. It is possible that running this
process just one time will not create an orthonormal basis U with the required accuracy
14
(p=15). In this case, we will need to re-run the process on the created matrix U and repeat it as
many times as needed until the required accuracy is archived. We proceed as follows.
**Initialize
count=1;
and set up an exterior “while” loop using the variable count that will keep track on how many
times we need to run the process of orthonormalization.
**Within the exterior loop, set up a “for” loop and use the formulas for Gram-Schmidt
orthogonalization (Lecture 30) to recalculate the columns of the matrix U, one by one, making
the following adjustments: after completing a consecutive column, you will need to normalize
it (use a MATLAB function norm()).
Hint: here we need to employ MATLAB operations such as matrix-vector multiplic