Skip to content

Commit

Permalink
Added demo codes for QR factorization
Browse files Browse the repository at this point in the history
  • Loading branch information
@drreynolds
drreynolds committed Jan 26, 2021
1 parent 6d426ba commit bbf1e06
Show file tree
Hide file tree
Showing 4 changed files with 316 additions and 0 deletions.
66 changes: 66 additions & 0 deletions QRFactorization/QRFactors.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,66 @@
function [Q,R] = QRFactors(A)
% usage: [Q,R] = QRFactors(A)
%
% Function to compute the QR factorization of a (possibly rank-deficient)
% 'thin' matrix A (m x n, with m >=n) using Householder reflection matrices.
%
% Input: A - thin matrix
% Outputs: Q - orthogonal matrix
% R - essentially upper triangular matrix, i.e. R = [ Rhat ]
% [ 0 ]
% with Rhat an (n x n) upper-triangular matrix
%
% Daniel R. Reynolds
% SMU Mathematics
% Math 4315

% get dimensions of A
[m,n] = size(A);

% initialize results
Q = eye(m);
R = A;

% determine elimination extent
if (m==n)
kend = n-1;
else
kend = n;
end

% iterate over columns
for k=1:kend

% extract subvector from diagonal down and compute norm
z = R(k:m,k);
v = -z;
v(1) = -sign(z(1))*norm(z) - z(1);
vnorm = norm(v);

% if subvector has norm zero, continue to next column
if (vnorm < eps)
continue
end

% compute u = v/||v||;
% the Householder matrix is then Qk = I-2*u*u'
u = v/vnorm;

% update rows k through m of R
for j=1:n
utR = 2*u'*R(k:m,j);
R(k:m,j) = R(k:m,j) - u*utR;
end

% update rows k through m of Q
for j=1:m
utQ = 2*u'*Q(k:m,j);
Q(k:m,j) = Q(k:m,j) - u*utQ;
end

end

% transpose Q before return
Q = Q';

% end function
69 changes: 69 additions & 0 deletions QRFactorization/QRFactors.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,69 @@
# QRFactors.py
#
# Daniel R. Reynolds
# SMU Mathematics
# Math 4315

# imports
import numpy

def QRFactors(A):
"""
usage: Q, R = QRFactors(A)
Function to compute the QR factorization of a (possibly rank-deficient)
'thin' matrix A (m x n, with m >=n) using Householder reflection matrices.
Input: A - thin matrix
Outputs: Q - orthogonal matrix
R - "upper triangular" matrix, i.e. R = [ Rhat ]
[ 0 ]
with Rhat an (n x n) upper-triangular matrix
"""

# get dimensions of A
m, n = numpy.shape(A)

# initialize results
Q = numpy.identity(m)
R = A.copy()

# determine elimination extent
if (m==n):
kend = n-1
else:
kend = n

# iterate over columns
for k in range(kend):

# extract subvector from diagonal down and compute norm
z = R[k:m,k]
v = -z;
v[0] = -numpy.sign(z[0])*numpy.linalg.norm(z) - z[0];
vnorm = numpy.linalg.norm(v)

# if subvector has norm zero, continue to next column
if (vnorm < numpy.finfo(float).eps):
continue

# compute u = u = v/||v||;
# the Householder matrix is then Qk = I-2*u*u'
u = v/vnorm

# update rows k through m of R
for j in range(n):
utR = 2 * numpy.transpose(u) @ R[k:m, j]
R[k:m, j] -= u*utR

# update rows k through m of Q
for j in range(m):
utQ = 2 * numpy.transpose(u) @ Q[k:m, j]
Q[k:m, j] -= u*utQ

# transpose Q before return
Q = numpy.transpose(Q)

return [Q, R]

# end function
91 changes: 91 additions & 0 deletions QRFactorization/TestQR.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,91 @@
% Script to test QRFactors on a variety of matrices.
%
% Daniel R. Reynolds
% SMU Mathematics
% Math 4315
clear

% set matrix sizes for tests
nvals = [50, 100, 200, 400];

% full-rank square matrix tests
for n = nvals

fprintf('Testing with full-rank square matrix of dimension %i\n',n);

% create the matrix
I = eye(n);
A = rand(n,n) + I;

% call QRFactors
[Q,R] = QRFactors(A);

% output results
fprintf(' ||I-Q^TQ|| = %g\n', norm(I-Q'*Q,2));
fprintf(' ||I-QQ^T|| = %g\n', norm(I-Q*Q',2));
fprintf(' ||A-QR|| = %g\n', norm(A-Q*R,2));
fprintf(' ||tril(R,-1)|| = %g\n', norm(tril(R,-1),2));

end

% full-rank thin matrix tests
for n = nvals

fprintf('Testing with full-rank rectangular matrix of dimension %ix%i\n',2*n,n);

% create the matrix
I = eye(2*n);
A = rand(2*n,n) + I(:,1:n);

% call QRFactors
[Q,R] = QRFactors(A);

% output results
fprintf(' ||I-Q^TQ|| = %g\n', norm(I-Q'*Q,2));
fprintf(' ||I-QQ^T|| = %g\n', norm(I-Q*Q',2));
fprintf(' ||A-QR|| = %g\n', norm(A-Q*R,2));
fprintf(' ||tril(R,-1)|| = %g\n', norm(tril(R,-1),2));

end

% rank-deficient square matrix tests
for n = nvals

fprintf('Testing with rank-deficient square matrix of dimension %i\n',n);

% create the matrix
I = eye(n);
A = rand(n,n) + I;
A(:,3) = 2*A(:,2);

% call QRFactors
[Q,R] = QRFactors(A);

% output results
fprintf(' ||I-Q^TQ|| = %g\n', norm(I-Q'*Q,2));
fprintf(' ||I-QQ^T|| = %g\n', norm(I-Q*Q',2));
fprintf(' ||A-QR|| = %g\n', norm(A-Q*R,2));
fprintf(' ||tril(R,-1)|| = %g\n', norm(tril(R,-1),2));

end

% rank-deficient thin matrix tests
for n = nvals

fprintf('Testing with rank-deficient rectangular matrix of dimension %ix%i\n',2*n,n);

% create the matrix
I = eye(2*n);
A = rand(2*n,n) + I(:,1:n);
A(:,3) = 2*A(:,2);

% call QRFactors
[Q,R] = QRFactors(A);

% output results
fprintf(' ||I-Q^TQ|| = %g\n', norm(I-Q'*Q,2));
fprintf(' ||I-QQ^T|| = %g\n', norm(I-Q*Q',2));
fprintf(' ||A-QR|| = %g\n', norm(A-Q*R,2));
fprintf(' ||tril(R,-1)|| = %g\n', norm(tril(R,-1),2));

end
90 changes: 90 additions & 0 deletions QRFactorization/TestQR.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,90 @@
#!/usr/bin/env python3
#
# Script to test QRFactors on a variety of matrices.
#
# Daniel R. Reynolds
# SMU Mathematics
# Math 4315

# imports
import numpy
from QRFactors import QRFactors

# set matrix sizes for tests
nvals = [50, 100, 200, 400]

# full-rank square matrix tests
for n in nvals:

print("Testing with full-rank square matrix of dimension ", n)

# create the matrix
I = numpy.eye(n)
A = numpy.random.rand(n,n) + I

# call QRFactors
Q, R = QRFactors(A)

# output results
print(" ||I-Q^TQ|| = ", numpy.linalg.norm(I-numpy.transpose(Q)@Q,2))
print(" ||I-QQ^T|| = ", numpy.linalg.norm(I-Q@numpy.transpose(Q),2))
print(" ||A-QR|| = ", numpy.linalg.norm(A-Q@R,2))
print(" ||tril(R,-1)|| = ", numpy.linalg.norm(numpy.tril(R,-1),2))

# full-rank rectangular matrix tests
for n in nvals:

print("Testing with full-rank rectangular matrix of dimension ", 2*n, "x", n)

# create the matrix
I = numpy.eye(2*n)
A = numpy.random.rand(2*n,n) + I[:,:n]

# call QRFactors
Q, R = QRFactors(A)

# output results
print(" ||I-Q^TQ|| = ", numpy.linalg.norm(I-numpy.transpose(Q)@Q,2))
print(" ||I-QQ^T|| = ", numpy.linalg.norm(I-Q@numpy.transpose(Q),2))
print(" ||A-QR|| = ", numpy.linalg.norm(A-Q@R,2))
print(" ||tril(R,-1)|| = ", numpy.linalg.norm(numpy.tril(R,-1),2))

# rank-deficient square matrix tests
for n in nvals:

print("Testing with rank-deficient square matrix of dimension ", n)

# create the matrix
I = numpy.eye(n)
A = numpy.random.rand(n,n) + I
A[:,2] = 2*A[:,1]

# call QRFactors
Q, R = QRFactors(A)

# output results
print(" ||I-Q^TQ|| = ", numpy.linalg.norm(I-numpy.transpose(Q)@Q,2))
print(" ||I-QQ^T|| = ", numpy.linalg.norm(I-Q@numpy.transpose(Q),2))
print(" ||A-QR|| = ", numpy.linalg.norm(A-Q@R,2))
print(" ||tril(R,-1)|| = ", numpy.linalg.norm(numpy.tril(R,-1),2))

# rank-deficient rectangular matrix tests
for n in nvals:

print("Testing with rank-deficient rectangular matrix of dimension ", 2*n, "x", n)

# create the matrix
I = numpy.eye(2*n)
A = numpy.random.rand(2*n,n) + I[:,:n]

# call QRFactors
Q, R = QRFactors(A)

# output results
print(" ||I-Q^TQ|| = ", numpy.linalg.norm(I-numpy.transpose(Q)@Q,2))
print(" ||I-QQ^T|| = ", numpy.linalg.norm(I-Q@numpy.transpose(Q),2))
print(" ||A-QR|| = ", numpy.linalg.norm(A-Q@R,2))
print(" ||tril(R,-1)|| = ", numpy.linalg.norm(numpy.tril(R,-1),2))


# end of script

0 comments on commit bbf1e06

Please sign in to comment.