Problem 8.2: Write a MATLAB function [Q,R] mgs(A) (see next lecture) that computes a reduced QR factorization A=QR of an m x n matrix A with m >= n using modifed Gram-Schmidt orthogonalization. The output variables are a martix Q (within complex domain m x n) with orthonormal columns and a triangular matrix R (within complex domain n x n).

Solution:

Matlab Code - QR by Modified Gram-Schmidt process algorithm

function [Q,R] = MGS(A)
% Usage: [Q,R] = MGS(A)
%
% This routine numerically approximates a thin QR decomposition, using 
% BLAS - 2 (Basic Linear Algebra Subroutine - 2), Modified Gran-Schmidt. 
% The routine also calculates the factorization and orthogonality errors
% accrued by the decomposition algrotihm.
%
% Inputs:  A      Randomly generated a m-by-n sized Matrix 
%                 
% Outputs: Q      The orthogonal matrix: Q
%          R      The upper triangular matrix: R
%

[~,n] = size(A);
Q = A;
R = zeros(n,n);

for j = 1:n
    R(j,j) = norm(Q(:,j));
    if (R(j,j) == 0)
        error('linearly dependent columns');
    end
    Q(:,j) = Q(:,j)/R(j,j);
    R(j,j+1:n) = Q(:,j)'*Q(:,j+1:n);
    Q(:,j+1:n) = Q(:,j+1:n) - Q(:,j)*R(j,j+1:n);
end
% Matlab Code #1
n=10; m=15;
A=randn(n,m);
[Q,R] = MGS(A);

Result: In this case we implemented a BLAS-2 (Basic Linear Algebra Subroutines-2) Modified Gram-Schmidt routine to approximate a thin QR decomposition. We constructed the algorithm around the previous implementation, and tested the decomposition routine on a pseudo-random normally distributed number sampled matrix of size 10 x 15, whom was generated by MATLAB’s built-in function randn(). The result of the simulation was a bit unexpected, as the decomposition error for this small matrix ranged from 10^(-14) to ϵ (Machine epsilon), an exceptional result. The results of this experiment should have been expected, as the Modified Gram-Schmidt process constructed in MGS substantially contributed to the routine’s conditioning and stability; resulting in significant less error (noise accumulation during computation).

Problem 10.2:

  1. Write a MATLAB function [W,R] house(A) that computes an implicit representation of a full QR factorization A=QR of an m x n matrix A with m >= n using Householder reflections. The output variables are a lower-triangular matrix W (within complex domain m x n) whose columns are the vectors V(k) defining the successive Householder reflections, and a triangular matrix R (within complex domain n x n).

Solution:

**Result:** 
    **(A)** In this experiment we construct a routine to approximate a full QR decomposition A=QR of a matrix A of size m x n (m≥n) using the Householder reflection method. But instead of constructing the corresponding unitary and orthogonal matrix Q; we instead constructed a matrix W to capture all the Householder transformations that were performed upon matrix A. The modification needed for the routine to capture the household reflections was quite simple, as we only needed to create an additional matrix to capture the householder transformation for each iteration of the for loop. At the end of this capturing process the W matrix that was constructed was a lower triangular matrix of size 10 x 10. Note, we attained this result through a simulation upon another matrix A which was generated through randn() function.

Matlab Code - QR by Householder reflection

function [W,R] = house(A)
% Usage: [W,R] = house(A)
%
% This routine numerically approximates a full QR Decomposition of A by
% Householder reflection.
%
% Inputs:  A      Randomly generated a m-by-n sized Matrix
%                 
% Outputs: Q      The orthogonal matrix: Q
%          R      The upper triangular matrix: R
%

[m,n] = size(A);
W=zeros(m,n);

for k = 1:n
    v=A(k:m,k);
    v(1)=v(1)+mysign(v(1))*norm(v);
    v=v/norm(v);
    A(k:m,k:n)=A(k:m,k:n)-2*v*(v'*A(k:m,k:n));
    W(k:m,k)=v;
end
R=triu(A);
end
% Matlab Code #2
n=10;
A=randn(n,n);
[W,R] = house(A)

  1. Write a MATLAB function Q=formQ(W) that takes the matrix W produced by house.m as input and generates a corresponding m x m orthogonal matrix Q.

    result (B) In the next part of this experiment we construct the orthogonal matrix Q from the householder transformation matrix W. Normally this line of executable code is implemented as the last line of iteration in a Householder processed QR decomposition, but due to the interest of the experimentation aims; it was implemented here as such. Nevertheless, the implementation of code is the same in operation, but indices and pattern of operation is quite different, as one would observe from the code below. The construction of the orthogonal matrix occurs from the row operations, instead of the normal column operations that is normally implemented due to the previous aforementioned implementation style discarded for the purposes of this experimentation. Nonetheless the new code construction previous computationally well-condition and stable and resulted in an orthogonal matrix with significantly negligible error, ranged from 10^(-15) to ϵ machine epsilon. All together this implementation of Householder reflection resulted with similar result to traditional academic and industry implementations of QR decomposition algorithms by householder reflection process.

Matlab Code - Used Householder transformation matrix to construct the orthogonal matrix Q

function [Q] = formQ(W)
% Usage: [Q] = formQ(W)
%
% This routine that uses the captured Householder transformations from
% matrix W to construct the orthogonal matrix Q.
%
% Inputs:  W      Lower-triangular matrix of Householder reflections
%                 
% Outputs: Q      The orthogonal matrix: Q
%

[m,n]=size(W);
Q = eye(m);

for j = n:-1:1
    Q(j:m,:)=Q(j:m,:)-2*W(j:m,j)*(W(j:m,j)'*Q(j:m,:));
end

end
% Matlab Code #3
Q = formQ(W)
[n,m]=size(A);
I = eye(n,m);
ErrorFact(1,1) = norm(A - (Q*R))
ErrorO(1,1) = norm(I - (Q'*Q))

Problem 10.3: Compute three reduced QR factorization of Z in MATLAB: by the Gram-Schmidt rountine mgs.m of Exercise 8.2, by Householder routines house.m and formQ.m of Excercise 10.2, and by MATLAB’s built-in command [Q,R]=qr(Z,0). Compare these three and comment on any differences you see.

Solution:

**Result:** Result: First difference observed among three is that MGS and MATLAB’s Economical QR decomposition algorithms both approximate a thin QR decomposition, while the house.m along accompanying formQ.m approximate a full QR decomposition. And as a result, both of the former algorithms form matrices of equal size (with reduce sized R matrices), while the later algorithm forms matrices equal to the matrix that it factorized. And a further implication is that QR by house.m must create two additional orthogonal vectors for orthogonal matrix Q, so that it an matrix operation between the Q and R can occur to attain the original matrix A. The last difference amongst the three QR decomposition algorithms is the sign of the elements if each of their vectors, and I speculate that it has to do with the type of process – such as GS, MGS, Householder, or Givens – that was used to construction that QR approximation.  
% Matlab Code #4
Z=[1 2 3;4 5 6;7 8 7;4 2 3;4 2 2];
[Q,R] = MGS(Z)
[W,R] = house(Z)
Q = formQ(W)
[Q,R]=qr(Z,0)
## Warning in matrix(unlist(Q3), nrow = 5, ncol = 3, byrow = FALSE): data length
## [25] is not a sub-multiple or multiple of the number of columns [3]

Problem 11.3: Take m=50, n=12. Using MATLAB’s linspace, define t to be the m-vector corresponding to linearly spaced grid points from 0 to 1.Using MATLAB’s vander.m and fliplr, define A to be the m x n matrix associated with least sqaures fitting on this grid by a polynomial of degree n-1. Take b to be the function cos(4t) evaluated on the grid. Now, calculate and print (to sixteen-digit precision) the least sqaures coefficient vector x by six methods:

Solution:

Note: I matrix A and vector b were too large for plot.matrix package. Thus I could not display either matrix without R screaming warnings, and preventing the kniting of the website. It is unfortunate, but something completely outside of my control.

**Result:** I noticed that Normal Equation had the greatest error from the six methods, and its coefficients were one magnitude off from the other methods. Noticed that mgs.m attained the correct magnitude, but had rounding error which separate it coefficients from the other methods; this was noticed when comparing mgs coefficients to the other methods coeifficients. Thus it was observed house.m, MATLAB's qr.m, MATLAB's \, MATLAB's SVD performed the similarly amongst each other, with slight differences with very minor rounding error in each of these method's coefficients.
The significant error by the Normal equation in solving this least square problem must be due to the inherit ill-condition-ness, and furthermore, instability that this algorithm possesses due to it pattern of operations. It is particularly not wise to perform Matrix-vector operations, as they are quite expensive, and also prone to rounding error. On the other hand Modified-Gram-Schmidt least square solution's error is due to the instability of minusing small numbers, and the disadvantage of forcing orthogonality amongst a set of vectors.
% Matlab Code #5
m=50; n=12;
t=linspace(0,1,m);

A=fliplr(vander(t));
A=A(:,1:n);
b=cos(4*t)';
## Warning in matrix(unlist(Error), nrow = 50, ncol = 1, byrow = FALSE): data
## length [6] is not a sub-multiple or multiple of the number of rows [50]
  1. Formation and solution of the normal equations, using MATLAB’S
% Matlab Code #6
[R,~] = Cholesky_Factorization(A'*A);
x1=R\(R'\(A'*b))
Error1 = norm((A*x1) - b)

  1. QR factorization computed by mgs.m (modified Gram-Schmid Exercise 8.2)
% Matlab Code #7
[Q,R] = MGS(A);
x2=R\(Q'*b)
Error2 = norm((A*x2) - b)

  1. QR factorization computed by house.m (Householder triangularization, Exercise 10.2)
% Matlab Code #8
[W1,R] = house(A);
Q = formQ(W1);
x3=R\(Q'*b)
Error3 = norm((A*x3) - b)

  1. QR factorization computed by MATLAB’s qr.m (also Householder triangularization)
% Matlab Code #9
[Q,R]=qr(A,0);
x4=R\(Q'*b)
Error4 = norm((A*x4) - b)

  1. x=A
% Matlab Code #10
x5=A\b
Error5 = norm((A*x5) - b)

  1. SVD, MATLAB’s svd.m
% Matlab Code #11
[U,S,V]=svd(A);
x6=V*(S\(U'*b))
Error6 = norm((A*x6) - b)

  1. The calculations above will produce six lists of twelve coefficients. In each list, shade with red the digits that appear to be wrong (affected by rounding error). Comment on what differences you observe. Do the normal equations exhibit instability? You do not have to explain your explain your obervations.
‘Screenshot of matrix X for the six x solution for the six methods to solve this Least Square’

‘Screenshot of matrix X for the six x solution for the six methods to solve this Least Square’

**Result:** I noticed that Normal Equation had the greatest error from the six methods, and its coefficients were one magnitude off from the other methods. Noticed that mgs.m attained the correct magnitude, but had rounding error which separate it coefficients from the other methods; this was noticed when comparing mgs coefficients to the other methods coeifficients. Thus it was observed house.m, MATLAB's qr.m, MATLAB's \, MATLAB's SVD performed the similarly amongst each other, with slight differences with very minor rounding error in each of these method's coefficients.
‘Kromia’

‘Kromia’

Below is the error for each of the six method used to solve this Least Sqaure problem. ‘Error for the Six Least Sqaure Solutions’

Extra Experiment: The extra experiment is to show the overal Factorization and Orthogonality error that each of the % accrued by all three QR decomposition algorithms tested here.

Result:

n=10:10:1000;
Factorization_Error=zeros(length(n),3); Orthogonality_Error=zeros(length(n),3);
FactError=zeros(length(n),1); OrthoError=zeros(length(n),1);

for i=1:length(n)
    A=randn(n(i),n(i));
    I = eye(n(i),n(i));
    [Q,R] = MGS(A);
    FactError(i)=norm(A - (Q*R));
    OrthoError(i)=norm(I - (Q'*Q));
end
Factorization_Error(:,1)=FactError(:);
Orthogonality_Error(:,1)=OrthoError(:);

for i=1:length(n)
    A=randn(n(i),n(i));
    I = eye(n(i),n(i));
    [W,R] = house(A);
    Q = formQ(W);
    FactError(i)=norm(A - (Q*R));
    OrthoError(i)=norm(I - (Q'*Q));
end
Factorization_Error(:,2)=FactError(:);
Orthogonality_Error(:,2)=OrthoError(:);

for i=1:length(n)
    A=randn(n(i),n(i));
    I = eye(n(i),n(i));
    [Q,R]=qr(A,0);
    FactError(i)=norm(A - (Q*R));
    OrthoError(i)=norm(I - (Q'*Q));
end
Factorization_Error(:,3)=FactError(:);
Orthogonality_Error(:,3)=OrthoError(:);

figure (1)
hold on
plot(n,Factorization_Error(:,1))
plot(n,Factorization_Error(:,2))
plot(n,Factorization_Error(:,3))
title('Factorization Error for QR Factorization Algorithms')
xlabel('size of Matrix A')
ylabel('Factorization Error')
legend({'MGS','House', 'qr'})
hold off


figure (2)
hold on
plot(n,Orthogonality_Error(:,1))
plot(n,Orthogonality_Error(:,2))
plot(n,Orthogonality_Error(:,3))
title('Orthogonality Error for QR factorization Algorithms')
xlabel('size of Matrix A')
ylabel('Orthogonality Error')
legend({'MGS','House', 'qr'})
hold off
‘Factorization Error’

‘Factorization Error’

‘Orthogonality Error’

‘Orthogonality Error’

Including Kromia (UwU)

It is Kromia. Because I wanted to include a cute picture at the end. ‘Kromia’