import numpy # # QR decomposition via classical Gram-Schmidt process # def gramschmidt_classical( A ): n = len( A ) Q = numpy.zeros( (n, n) ) R = numpy.zeros( (n, n) ) for i in range( 0, n ): ai = A[:,i] qi = A[:,i] for j in range(0,i): qj = Q[:,j] R[j,i] = ai.dot( qj ) #take dot product with original column. qi = qi - R[j,i] * qj R[i,i] = numpy.sqrt( qi.dot( qi ) ) Q[:,i] = qi / R[i,i] return (Q,R) # # QR decomposition via modified Gram-Schmidt process # def gramschmidt_modified(A): n = len( A ) Q = numpy.zeros( (n, n) ) R = numpy.zeros( (n, n) ) for i in range( 0, n ): ai = A[:,i] qi = ai for j in range(0,i): qj = Q[:,j] R[j,i] = qi.dot( qj ) #take dot product with current partially orthogonalized column. qi = qi - R[j,i] * qj R[i,i] = numpy.sqrt( qi.dot( qi ) ) Q[:,i] = qi / R[i,i] return (Q,R) # # Given a matrix A and a QR decomposition algorithm 'method' # compute the decomposition A = Q R # and return the Frobenius norms of # A - Q R # I - Q^t Q # def test_decomposition( A, method ): N = len( A ) (Q,R) = method( A ) A_diff = A - numpy.matmul( Q, R ) Q_diff = numpy.identity( N ) - numpy.matmul( numpy.transpose( Q ), Q ) A_error = numpy.linalg.norm( A_diff, ord='fro' ) Q_error = numpy.linalg.norm( Q_diff, ord='fro' ) return (A_error,Q_error) # # Test the QR decomposition via classical and modified Gram-Schmidt process # with random (Gaussian) matrices. # Output the errors ratios and their averages. # def run_tests( N, T ): # Arrays for the error terms of the QR decomposition methods ... A_errors_classic = numpy.zeros( (T) ) Q_errors_classic = numpy.zeros( (T) ) A_errors_modified = numpy.zeros( (T) ) Q_errors_modified = numpy.zeros( (T) ) # ... and their ratios. A_errors_ratios = numpy.zeros( (T) ) Q_errors_ratios = numpy.zeros( (T) ) for t in range(0,T): # Draw a random matrix. A = numpy.random.normal( size = (N,N) ) # Uncomment the code line below to exacerbate # the conditioning of the matrix A. # #A = numpy.matmul(A,A) # Compute the error terms of the two decomposition methods. (A_c,Q_c) = test_decomposition( A, gramschmidt_classical ) (A_m,Q_m) = test_decomposition( A, gramschmidt_modified ) # Save the error terms and print them A_errors_classic[t] = A_c Q_errors_classic[t] = Q_c A_errors_modified[t] = A_m Q_errors_modified[t] = Q_m A_errors_ratios[t] = A_c / A_m Q_errors_ratios[t] = Q_c / Q_m print( A_c/A_m, Q_c/Q_m ) # The for-loop having finished, we print the averages of the error ratios. print numpy.mean( A_errors_ratios ) print numpy.mean( Q_errors_ratios ) N = 300 T = 100 run_tests(N,T)