With this notebook, you will implement the QR factorization of a matrix $ A $ with linearly independent columns, producing matrices $ Q $ and $ R $ such that $ A = Q R $. The algorithm is equivalent to performing Gram-Schmidt orthogonalization on the vectors that are the columns of $ A $.
Be sure to make a copy!!!!
import numpy as np
import laff
import flame
A = np.matrix( ' 1., -1., 2;\
2., 1., -3;\
-1., 3., 2;\
0., -2., -1' )
Q = np.matrix( np.zeros( (4,3) ) )
R = np.matrix( np.zeros( (3,3) ) )
print( 'A = ' )
print( A )
A = [[ 1. -1. 2.] [ 2. 1. -3.] [-1. 3. 2.] [ 0. -2. -1.]]
Here is the algorithm:
Important: if you make a mistake, rerun ALL cells above the cell in which you were working, and then the one where you are working.
Create the routine
QR_Gram_Schmidt_unb
with the Spark webpage for the algorithm
Hints:
q1
. This will mean first copying a1
to q1
.
laff.copy ( x, y )
laff.gemv ( trans, alpha, A, x, beta, y )
laff.norm2 ( x )
rho11[:,:] = laff.norm2
.
laff.invscal ( alpha, x )
# Insert code here
Important: if you make a mistake, rerun ALL cells above the cell in which you were working, and then the one where you are working.
QR_Gram_Schmidt_unb( A, Q, R )
print( 'A = ' )
print( A )
print( 'Q = ' )
print( Q )
print( 'R = ' )
print( R )
print( 'Q * R - A:' )
print( Q * R - A )
A = [[ 1. -1. 2.] [ 2. 1. -3.] [-1. 3. 2.] [ 0. -2. -1.]] Q = [[ 0.40824829 -0.17609018 0.8820199 ] [ 0.81649658 0.44022545 -0.32318287] [-0.40824829 0.70436073 0.23565417] [ 0. -0.52827054 -0.24912013]] R = [[ 2.44948974 -0.81649658 -2.44948974] [ 0. 3.7859389 0.26413527] [ 0. 0. 3.45401687]] Q * R - A: [[ 0.00000000e+00 0.00000000e+00 -2.22044605e-16] [ 0.00000000e+00 0.00000000e+00 0.00000000e+00] [ 0.00000000e+00 0.00000000e+00 0.00000000e+00] [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]]
The result should be a 4 x 3 (approximately) zero matrix.
Now, let's check if the columns of $ Q $ are mutually orthogonal:
print( np.transpose( Q ) * Q )
[[ 1.00000000e+00 1.11022302e-16 1.66533454e-16] [ 1.11022302e-16 1.00000000e+00 2.77555756e-17] [ 1.52655666e-16 5.55111512e-17 1.00000000e+00]]
The above should approximately be a 3 x 3 identity matrix.
If you read the enrichment on the Gram-Schmidt method, you will find the following interesting...
import numpy as np
import laff
import flame
epsilon = 1.0e-7
A = np.matrix( ' 1., 1., 1.;\
0, 0., 0.;\
0., 0, 0.;\
0., 0., 0.' )
A[ 1,0 ] = epsilon
A[ 2,1 ] = epsilon
A[ 3,2 ] = epsilon
# This creates the matrix
# A = np.matrix( ' 1., 1., 1.;\
# epsilon, 0., 0.;\
# 0., epsilon, 0.;\
# 0., 0., epsilon' )
Q = np.matrix( np.zeros( (4,3) ) )
R = np.matrix( np.zeros( (3,3) ) )
print( 'A = ' )
print( A )
A = [[ 1.00000000e+00 1.00000000e+00 1.00000000e+00] [ 1.00000000e-07 0.00000000e+00 0.00000000e+00] [ 0.00000000e+00 1.00000000e-07 0.00000000e+00] [ 0.00000000e+00 0.00000000e+00 1.00000000e-07]]
QR_Gram_Schmidt_unb( A, Q, R )
print( 'A = ' )
print( A )
print( 'Q = ' )
print( Q )
print( 'R = ' )
print( R )
print( 'Q * R - A:' )
print( Q * R - A )
A = [[ 1.00000000e+00 1.00000000e+00 1.00000000e+00] [ 1.00000000e-07 0.00000000e+00 0.00000000e+00] [ 0.00000000e+00 1.00000000e-07 0.00000000e+00] [ 0.00000000e+00 0.00000000e+00 1.00000000e-07]] Q = [[ 1.00000000e+00 6.90840682e-08 4.07886015e-08] [ 1.00000000e-07 -7.07106781e-01 -4.17602698e-01] [ 0.00000000e+00 7.07106781e-01 -3.98821881e-01] [ 0.00000000e+00 0.00000000e+00 8.16424579e-01]] R = [[ 1.00000000e+00 1.00000000e+00 1.00000000e+00] [ 0.00000000e+00 1.41421356e-07 6.90840682e-08] [ 0.00000000e+00 0.00000000e+00 1.22485288e-07]] Q * R - A: [[ 0. 0. 0.] [ 0. 0. 0.] [ 0. 0. 0.] [ 0. 0. 0.]]
This should equal, approximately, the zero matrix.
Now, let's check if the columns of Q are mutually orthogonal:
print( np.transpose( Q ) * Q )
[[ 1.00000000e+00 -1.62660994e-09 -9.71668373e-10] [ -1.62660994e-09 1.00000000e+00 1.32800433e-02] [ -9.71668373e-10 1.32800433e-02 1.00000000e+00]]
You will notice that the resulting 3 x 3 matrix, which should (approximately) equal the identity matrix, has off-diagonal entries that do not equal zero at all.
One of them even equals 0.5 (when I tried)