UW AMath

AMath 584, Autumn Quarter, 2011

Applied Linear Algebra and

Introductory Numerical Analysis

Table Of Contents

Python

Demo 1

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).

[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[0] * 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 [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