单步QR方法计算上Hessenberg矩阵特征值的Python实现
时间：2008-11-15 来源：taoyuliang
刚开通博客，把今晚写的一个Python程序传了上来。
这个类实现了用单步QR分解来计算上Hessenberg矩阵的全部特征值。QR方法是一种计算一般矩阵（中小型矩阵）全部特征值问最有效的方法之一。在很多数值分析的书本上都有该算法的说明，所以在这里不慢慢解释。
其中，在import中代入的几个模块Matrix，Househoulder和Hessenberg，均是用Python实现的几个类，在后续几篇博客中，我将慢慢介绍这些类。
1 #!/usr/bin/python
2
3 ################################################################
4 # Filename : eigenvalueQR.py
5 # Compute the eigenvalue of a matrix with QR method
6 ################################################################
7
8 # import
9 from copy import deepcopy
10 from matrix import Matrix
11 from householder import Householder
12 from hessenberg import Hessenberg
13
14 ################################################################
15 #============= CLASS EIGENVALUEQR DEFINITION BEGIN ============#
16 ################################################################
17 class EigenvalueQR :
18 """Compute the eigenvalue of a matrix with QR decomposition"""
19
20 def __init__( self, matrix, eps = None ) :
21 """Constructor of class"""
22
23 # matrix is must be a squre matrix
24 if not matrix.isSquareMatrix() :
25 raise Exception, "Square matrix is required"
26 return
27
28 self.__n = matrix.getdimRow()
29
30 # define the precision
31 if eps :
32 self.__eps = eps
33 else :
34 self.__eps = 0.0001
35
36 # Judge if matrix is a upper hessenberg matrix
37 # If not, implement the hessenberg transfromation
38 flag = self._isHessenbergMatrix( matrix )
39 # is not a irreducible hessenberg matrix
40 if flag == 1 :
41 raise Exception, "Must be irreducible Hessenberg Matrix"
42 return
43 # is not a upper hessenberg matrix
44 elif flag == 0 :
45 hb = Hessenberg( matrix )
46 self.__H = hb.getH()
47 # is a irreducible upper hessenberg matrix
48 else :
49 self.__H = deepcopy( matrix )
50
51 # now the H is a irreducible upper hessenberg matrix
52 # we can implement QR decomposition
53 # control the condition for termination
54 iterator = self.__n - 1
55 subH = self.__H
56 self.__eigenvalues = []
57
58 while iterator :
59 # one step QR mwthod
60 eigenv, subH = self._onestepQR( subH )
61 self.__eigenvalues.append( eigenv )
62 iterator -= 1
63
64 # get the last eigenvalue
65 self.__eigenvalues.append( subH[0][0] )
66
67
68 #-----------------------------------------------------#
69 # Define utility methods of class
70 #-----------------------------------------------------#
71 def _isHessenbergMatrix( self, A ) :
72 """Judge if A is a hessenberg matrix"""
73
74 dim = A.getdimRow()
75
76 for i in range( 2, dim ) :
77 for j in range( i-1 ) :
78 # must be hessenberg matrix
79 if A[i][j] != 0 :
80 return 0
81
82 for i in range( 1, dim ) :
83 # must be irreducible hessenberg matrix
84 if A[i][i-1] == 0 :
85 return 1
86
87 return 2
88
89 def _onestepQR( self, matrix ) :
90 """Compute the upper hessenberg matrix's eigenvalue
91 with onw step QR method"""
92
93 #-------------------------------------------------
94 # Argorithm :
95 #-------------------------------------------------
96 # H[k] - s[k] * I = Q[k] * R[k]
97 # H[k+1] = R[k] * Q[k] + s[k] * I
98 #-------------------------------------------------
99 H = deepcopy(matrix)
100
101 dim = H.getdimRow()
102 n = dim - 1
103
104 while 1 :
105 # construct H[k] - s[k] * I
106 # get s[k]
107 sk = H[n][n]
108 offset = self._genDiagMatrix( dim, sk )
109 H -= offset
110
111 # implement QR decomposition for H
112 hh = Householder( H )
113 Q = hh.getQ()
114 R = hh.getR()
115
116 # update H : H[k+1] = R * Q + s[k] * I
117 H = R * Q + offset
118
119 #
120 if H[n][n-1] < self.__eps :
121 print "H:\n", H
122 break
123 else :
124 print "H:\n", H
125
126 return ( H[n][n], H.getSubMatrix(1,n,1,n) )
127
128
129 def _genDiagMatrix( self, dim, k = None ) :
130 """Generate a diag matrix. If k is None, then
131 generate a unit matrix"""
132
133 I = Matrix( dim, dim )
134 cons = 1.0
135 if k :
136 cons = k
137
138 for i in range( dim ) :
139 I[i][i] = cons
140
141 return I
142
143
144 #-----------------------------------------------------#
145 # Define user interface methods of class
146 #-----------------------------------------------------#
147 def getEigenvalues( self ) :
148 """Return the list of all eigenvalues"""
149
150 return self.__eigenvalues
151
152 def getDimension( self ) :
153 """Return the dimension of matrix"""
154
155 return self.__n
156
157 def getEpsilon( self ) :
158 """Return the required precision epsilon"""
159
160 return self.__eps
161
162
163 ################################################################
164 #============== CLASS EIGENVALUEQR DEFINITION END =============#
165 ################################################################
166 def main() :
167 """Test"""
168
169 data = [ [ 2, 1, 0],
170 [ 1, 3, 1],
171 [ 0, 1, 4] ]
172
173 H = Matrix( 3, 3, data )
174 print "Original matrix: "
175 print H
176 print
177
178 qr = EigenvalueQR( H )
179 values = qr.getEigenvalues()
180 print "Eigenvalues : "
181 for i in range( len(values) ) :
182 print "%8f" % values[i]
183
184 print
185
186 ################################################################
187 if __name__ == '__main__' :
188 main()
189
这个类实现了用单步QR分解来计算上Hessenberg矩阵的全部特征值。QR方法是一种计算一般矩阵（中小型矩阵）全部特征值问最有效的方法之一。在很多数值分析的书本上都有该算法的说明，所以在这里不慢慢解释。
其中，在import中代入的几个模块Matrix，Househoulder和Hessenberg，均是用Python实现的几个类，在后续几篇博客中，我将慢慢介绍这些类。
1 #!/usr/bin/python
2
3 ################################################################
4 # Filename : eigenvalueQR.py
5 # Compute the eigenvalue of a matrix with QR method
6 ################################################################
7
8 # import
9 from copy import deepcopy
10 from matrix import Matrix
11 from householder import Householder
12 from hessenberg import Hessenberg
13
14 ################################################################
15 #============= CLASS EIGENVALUEQR DEFINITION BEGIN ============#
16 ################################################################
17 class EigenvalueQR :
18 """Compute the eigenvalue of a matrix with QR decomposition"""
19
20 def __init__( self, matrix, eps = None ) :
21 """Constructor of class"""
22
23 # matrix is must be a squre matrix
24 if not matrix.isSquareMatrix() :
25 raise Exception, "Square matrix is required"
26 return
27
28 self.__n = matrix.getdimRow()
29
30 # define the precision
31 if eps :
32 self.__eps = eps
33 else :
34 self.__eps = 0.0001
35
36 # Judge if matrix is a upper hessenberg matrix
37 # If not, implement the hessenberg transfromation
38 flag = self._isHessenbergMatrix( matrix )
39 # is not a irreducible hessenberg matrix
40 if flag == 1 :
41 raise Exception, "Must be irreducible Hessenberg Matrix"
42 return
43 # is not a upper hessenberg matrix
44 elif flag == 0 :
45 hb = Hessenberg( matrix )
46 self.__H = hb.getH()
47 # is a irreducible upper hessenberg matrix
48 else :
49 self.__H = deepcopy( matrix )
50
51 # now the H is a irreducible upper hessenberg matrix
52 # we can implement QR decomposition
53 # control the condition for termination
54 iterator = self.__n - 1
55 subH = self.__H
56 self.__eigenvalues = []
57
58 while iterator :
59 # one step QR mwthod
60 eigenv, subH = self._onestepQR( subH )
61 self.__eigenvalues.append( eigenv )
62 iterator -= 1
63
64 # get the last eigenvalue
65 self.__eigenvalues.append( subH[0][0] )
66
67
68 #-----------------------------------------------------#
69 # Define utility methods of class
70 #-----------------------------------------------------#
71 def _isHessenbergMatrix( self, A ) :
72 """Judge if A is a hessenberg matrix"""
73
74 dim = A.getdimRow()
75
76 for i in range( 2, dim ) :
77 for j in range( i-1 ) :
78 # must be hessenberg matrix
79 if A[i][j] != 0 :
80 return 0
81
82 for i in range( 1, dim ) :
83 # must be irreducible hessenberg matrix
84 if A[i][i-1] == 0 :
85 return 1
86
87 return 2
88
89 def _onestepQR( self, matrix ) :
90 """Compute the upper hessenberg matrix's eigenvalue
91 with onw step QR method"""
92
93 #-------------------------------------------------
94 # Argorithm :
95 #-------------------------------------------------
96 # H[k] - s[k] * I = Q[k] * R[k]
97 # H[k+1] = R[k] * Q[k] + s[k] * I
98 #-------------------------------------------------
99 H = deepcopy(matrix)
100
101 dim = H.getdimRow()
102 n = dim - 1
103
104 while 1 :
105 # construct H[k] - s[k] * I
106 # get s[k]
107 sk = H[n][n]
108 offset = self._genDiagMatrix( dim, sk )
109 H -= offset
110
111 # implement QR decomposition for H
112 hh = Householder( H )
113 Q = hh.getQ()
114 R = hh.getR()
115
116 # update H : H[k+1] = R * Q + s[k] * I
117 H = R * Q + offset
118
119 #
120 if H[n][n-1] < self.__eps :
121 print "H:\n", H
122 break
123 else :
124 print "H:\n", H
125
126 return ( H[n][n], H.getSubMatrix(1,n,1,n) )
127
128
129 def _genDiagMatrix( self, dim, k = None ) :
130 """Generate a diag matrix. If k is None, then
131 generate a unit matrix"""
132
133 I = Matrix( dim, dim )
134 cons = 1.0
135 if k :
136 cons = k
137
138 for i in range( dim ) :
139 I[i][i] = cons
140
141 return I
142
143
144 #-----------------------------------------------------#
145 # Define user interface methods of class
146 #-----------------------------------------------------#
147 def getEigenvalues( self ) :
148 """Return the list of all eigenvalues"""
149
150 return self.__eigenvalues
151
152 def getDimension( self ) :
153 """Return the dimension of matrix"""
154
155 return self.__n
156
157 def getEpsilon( self ) :
158 """Return the required precision epsilon"""
159
160 return self.__eps
161
162
163 ################################################################
164 #============== CLASS EIGENVALUEQR DEFINITION END =============#
165 ################################################################
166 def main() :
167 """Test"""
168
169 data = [ [ 2, 1, 0],
170 [ 1, 3, 1],
171 [ 0, 1, 4] ]
172
173 H = Matrix( 3, 3, data )
174 print "Original matrix: "
175 print H
176 print
177
178 qr = EigenvalueQR( H )
179 values = qr.getEigenvalues()
180 print "Eigenvalues : "
181 for i in range( len(values) ) :
182 print "%8f" % values[i]
183
184 print
185
186 ################################################################
187 if __name__ == '__main__' :
188 main()
189
相关阅读 更多 +
