|
AMath 584, Autumn Quarter, 2011 Applied Linear Algebra and Introductory Numerical Analysis |
|
|
AMath 584, Autumn Quarter, 2011 Applied Linear Algebra and Introductory Numerical Analysis |
See also the Matlab version of Demo 1 for comparison.
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).
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)
|
Note that Python differs from Matlab in the following ways:
A more complicated example than shown in Matlab.
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[0] * Vh[0,:];
figure(2)
clf()
pcolor(flipud(array(A1)), cmap=cm.RdYlBu, edgecolors='k')
xlabel('Movies')
ylabel('People')
title('Rank 1 approximation')
colorbar()
|
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 [151]: 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