diff --git a/QRFactorization/QRFactors.m b/QRFactorization/QRFactors.m new file mode 100644 index 0000000..a1a95c7 --- /dev/null +++ b/QRFactorization/QRFactors.m @@ -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 diff --git a/QRFactorization/QRFactors.py b/QRFactorization/QRFactors.py new file mode 100644 index 0000000..f14986e --- /dev/null +++ b/QRFactorization/QRFactors.py @@ -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 diff --git a/QRFactorization/TestQR.m b/QRFactorization/TestQR.m new file mode 100644 index 0000000..e76c3bf --- /dev/null +++ b/QRFactorization/TestQR.m @@ -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 diff --git a/QRFactorization/TestQR.py b/QRFactorization/TestQR.py new file mode 100755 index 0000000..a7b7786 --- /dev/null +++ b/QRFactorization/TestQR.py @@ -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