Source Code for Module SloppyCell.Subspaces

 1  import scipy 
 2   
 3 -def subspace_angles(A, B): 
 4      """ 
 5      Return the principle angles and vectors between two subspaces.  
 6   
 7      See Bjorck and Golub, "Numerical Methods for Computing Between Linear  
 8             Subspaces" _Mathematics_of_Computation_, 27(123), 1973, pp. 579-594 
 9      Or, for a more understandable exposition, Golub and Van Loan,  
10             _Matrix_Computations_ (3rd ed.) pp. 603-604 
11      """ 
12      A, B = scipy.asarray(A), scipy.asarray(B) 
13      if A.shape[0] != B.shape[0]: 
14          raise ValueError, 'Input subspaces must live in the same dimensional '\ 
15                  'space.' 
16   
17      # Get orthogonal bases for our two subspaces. 
18      QA, QB = scipy.linalg.orth(A), scipy.linalg.orth(B) 
19   
20      M = scipy.matrixmultiply(scipy.transpose(scipy.conjugate(QA)), QB) 
21      Y, C, Zh = scipy.linalg.svd(M) 
22   
23      U = scipy.matrixmultiply(QA, Y) 
24      V = scipy.matrixmultiply(QB, scipy.transpose(scipy.conjugate(Zh))) 
25      return scipy.arccos(C), U, V 
26   
27 -def test(): 
28      # This is the example from Bjorck and Golub (1973) 
29      m, p = 26, 13 
30   
31      A = scipy.zeros((m, p), scipy.Float) 
32      for col in range(p): 
33          A[col*(m/p):(col+1)*(m/p), col:col+1] = scipy.ones((m/p, 1)) 
34      A *= 1/scipy.sqrt(m/p) 
35   
36      B = scipy.zeros((m, p), scipy.Complex) 
37      for col in range(0, p): 
38          B[:, col:col+1] = scipy.transpose(scipy.mat((-1 + 2j*scipy.arange(m)/(m+1))**col)) 
39   
40      theta, U, V = subspace_angles(A, B) 
41      assert scipy.round(scipy.cos(theta), 4) == \ 
42              scipy.array([ 1.    ,  0.9982,  0.9981,  0.9903,  0.9899,  0.9765, 
43                           0.9628,  0.9415, 0.9176,  0.8701,  0.7637,  0.0608,   
44                           0.0156])  
45      assert scipy.round(scipy.sin(theta), 4) == \ 
46              scipy.array([ 0.    ,  0.0594,  0.0609,  0.1388,  0.1418,  0.2157, 
47                           0.2701,  0.337 , 0.3975,  0.4928,  0.6456,  0.9982,   
48                           0.9999]) 
49      print 'Comparison with Bjorck and Golub (1973) example passed.' 
50   
51      # Example from p 605 of Matrix Computations 
52      A = [[1,2],[3,4],[5,6]] 
53      B = [[1,5],[3,7],[5,-1]] 
54      theta, U, V = subspace_angles(A, B) 
55      assert scipy.round(scipy.cos(scipy.real(theta)), 3) ==\ 
56              scipy.array([1.000, 0.856]) 
57       
58 -def plot_test(): 
59      A = scipy.rand(3, 2) - 0.5 
60      B = scipy.rand(3, 1) - 0.5 
61      theta, U, V = subspace_angles(A, B) 
62       
63      import os, tempfile 
64      file_names = [] 
65      for vects in [A, B, U, V]: 
66          fd, name = tempfile.mkstemp() 
67          file_names.append(name) 
68          f = os.fdopen(fd, 'w') 
69          lines = [] 
70          for column in scipy.transpose(vects): 
71              column = scipy.asarray(column) 
72              normalized = column/scipy.sqrt(scipy.dot(column, column)) 
73              lines.append(' '.join([str(0) for elem in column])) 
74              lines.append(' '.join([str(elem) for elem in normalized])) 
75              lines.append('') 
76              lines.append('') 
77          f.write(os.linesep.join(lines)) 
78          f.close() 
79       
80      gnuplot_command = 'splot "%s" lw 10 t "A", "%s" lw 10 t "B", '\ 
81              '"%s" lw 8 t "U", "%s" lw 4 t "V"' % tuple(file_names) 
82   
83      print theta * 180/scipy.pi 
84      print gnuplot_command 
85