 AMath 584, Autumn Quarter, 2011 Applied Linear Algebra and Introductory Numerical Analysis

Python¶

Demo 1¶

The Python version looks a bit more complicated because older versions of NumPy do not provide commands rank or null. These are computed here by using the SVD.

Note that in NumPy the rank of an array generally refers to the number of indices the array takes (e.g. 1 for a vector, 2 for a matrix), not to the matrix rank. More recent versions have a function numpy.linalg.matrix_rank that computes the matrix rank (using the SVD).

[demo1.py]

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 from numpy import * from numpy.linalg import * caseno = 1 if caseno == 1: A = matrix('[1,2,1; 2,4,1; 3,6,1; 4,8,1]') elif caseno == 2: A = matrix('[1+1j, 1; 2+2j, 2; 3+3j, 3]') print "\nA = " print A print "\nThe adjoint of A is:" print A.H print "\nThe matrix (adjoint of A) * A is:" print A.H * A U,s,Vh = svd(A) # singular value decomposition tolerance = 1.e-12 # set smaller singular values to 0 rankA = sum(s > tolerance) # number of nonzero singular values print "\nThe rank of A is %g" % rankA # Select columns of Vh corresponding to non-zero singular values: Vh_zero_s = compress(s < tolerance, Vh, axis=0) null_space = Vh_zero_s.H print "\nBasis for the null space of A is:" print null_space x = null_space print "\nx = " print x b = A*x print "\nb = " print b print "\n2-norm of A is %g" % norm(A) print "1-norm of A is %g" % norm(A,1) print "max-norm of A is %g" % norm(A,inf)

The SVD in Python¶

Note that Python differs from Matlab in the following ways:

• The s returned is a vector of singular values rather than a matrix.
• The matrix Vh returned is the conjugate transpose of V, whereas Matlab returns V.

Recommender systems¶

A more complicated example than shown in Matlab.

[recommend2.py]

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 """ To use, start Ipython with the -pylab flag: \$ ipython -pylab and then In[1}: run recommend2.py """ from pylab import * from matplotlib import cm # colormaps A = matrix("""[-2 0 -1 1 0 2 0; 0 -1 -2 0 1 0 2; -2 1 -2 2 0 1 2; -2 0 1 1 2 1 0; 2 2 0 1 -2 0 -1; 0 1 2 2 -1 -1 0; 1 1 2 0 0 0 -2; 0 0 0 0 2 0 0]""") (U, s, Vh) = svd(A) figure(1) clf() pcolor(flipud(array(A)), cmap=cm.RdYlBu, edgecolors='k') xlabel('Movies') ylabel('People') title('Original data') colorbar() A1 = U[:,0] * s * Vh[0,:]; figure(2) clf() pcolor(flipud(array(A1)), cmap=cm.RdYlBu, edgecolors='k') xlabel('Movies') ylabel('People') title('Rank 1 approximation') colorbar()

Rayleigh Quotient Iteration¶

[rqiter.py]

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 """ Sample code for Rayleigh quotient iteration applied to shift matrix. """ from pylab import * A = array([[0,1,0,0],[0,0,1,0],[0,0,0,1],[1,0,0,0]]) print "\nA = " print A lam,X = eig(A) print "\nEigenvalues: ",lam v = randn(4) print "\nStarting guess v: ",v print " " mu = 1. + 1j print "k = %s, mu = %21.16f + %21.16fj" % (0,mu.real,mu.imag) for k in range(1,10): B = A - mu*eye(4) try: w = solve(B,v) except: print "Matrix is singular! Converged solution" break v = w / norm(w,2) mu = dot(conj(v.T), dot(A,v)) print "k = %s, mu = %21.16f + %21.16fj" % (k,mu.real,mu.imag)

Gives output like:

In : run rqiter.py

A =
[[0 1 0 0]
[0 0 1 0]
[0 0 0 1]
[1 0 0 0]]

Eigenvalues:  [ -1.00000000e+00+0.j  -5.55111512e-17+1.j
-5.55111512e-17-1.j
1.00000000e+00+0.j]

Starting guess v:  [ 1.41419139  0.85394083  0.04821873 -0.04396223]

k = 0, mu =    1.0000000000000000 +    1.0000000000000000j
k = 1, mu =    0.6006347079781388 +    0.2528188603994022j
k = 2, mu =    0.8765977589170861 +    0.0999099526208390j
k = 3, mu =    0.9978231303333175 +    0.0018522523017547j
k = 4, mu =    0.9999999910960606 +    0.0000000076154739j
k = 5, mu =    1.0000000000000000 +   -0.0000000000000000j
k = 6, mu =    1.0000000000000000 +    0.0000000000000000j
Matrix is singular!   Converged solution