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 )
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 )
import flame
import laff as laff
def QR_Gram_Schmidt_unb(A, Q, R):
AL, AR = flame.part_1x2(A, \
0, 'LEFT')
QL, QR = flame.part_1x2(Q, \
0, 'LEFT')
RTL, RTR, \
RBL, RBR = flame.part_2x2(R, \
0, 0, 'TL')
while AL.shape[1] < A.shape[1]:
A0, a1, A2 = flame.repart_1x2_to_1x3(AL, AR, \
1, 'RIGHT')
Q0, q1, Q2 = flame.repart_1x2_to_1x3(QL, QR, \
1, 'RIGHT')
R00, r01, R02, \
r10t, rho11, r12t, \
R20, r21, R22 = flame.repart_2x2_to_3x3(RTL, RTR, \
RBL, RBR, \
1, 1, 'BR')
#------------------------------------------------------------#
laff.gemv( 'Transpose', 1.0, Q0, a1, 0.0, r01 )
laff.copy( a1, q1 )
laff.gemv( 'No transpose', -1.0, Q0, r01, 1.0, q1 )
rho11[:,:] = laff.norm2( q1 )
laff.invscal( rho11, q1 )
#------------------------------------------------------------#
AL, AR = flame.cont_with_1x3_to_1x2(A0, a1, A2, \
'LEFT')
QL, QR = flame.cont_with_1x3_to_1x2(Q0, q1, Q2, \
'LEFT')
RTL, RTR, \
RBL, RBR = flame.cont_with_3x3_to_2x2(R00, r01, R02, \
r10t, rho11, r12t, \
R20, r21, R22, \
'TL')
flame.merge_1x2(QL, QR, Q)
flame.merge_2x2(RTL, RTR, \
RBL, RBR, R)
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 )
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 )
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 )
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 )
This should equal, approximately, the zero matrix.
Now, let's check if the columns of Q are mutually orthogonal:
print( np.transpose( Q ) * Q )
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)