单步QR方法计算上Hessenberg矩阵特征值的Python实现
taoyuliang
|
1#
taoyuliang 发表于 2008-11-15 00:33
单步QR方法计算上Hessenberg矩阵特征值的Python实现刚开通博客，把今晚写的一个Python程序传了上来。 这个类实现了用单步QR分解来计算上Hessenberg矩阵的全部特征值。QR方法是一种计算一般矩阵（中小型矩阵）全部特征值问最有效的方法之一。在很多数值分析的书本上都有该算法的说明，所以在这里不慢慢解释。 其中，在import中代入的几个模块Matrix，Househoulder和Hessenberg，均是用Python实现的几个类，在后续几篇博客中，我将慢慢介绍这些类。[color="#804040"] [color="#804040"] 1 [color="#0000ff"]#!/usr/bin/python [color="#804040"] 2 [color="#804040"] 3 [color="#0000ff"]################################################################ [color="#804040"] 4 [color="#0000ff"]# Filename : eigenvalueQR.py [color="#804040"] 5 [color="#0000ff"]# Compute the eigenvalue of a matrix with QR method [color="#804040"] 6 [color="#0000ff"]################################################################ [color="#804040"] 7 [color="#804040"] 8 [color="#0000ff"]# import [color="#804040"] 9 [color="#a020f0"]from copy [color="#a020f0"]import deepcopy [color="#804040"] 10 [color="#a020f0"]from matrix [color="#a020f0"]import Matrix [color="#804040"] 11 [color="#a020f0"]from householder [color="#a020f0"]import Householder [color="#804040"] 12 [color="#a020f0"]from hessenberg [color="#a020f0"]import Hessenberg [color="#804040"] 13 [color="#804040"] 14 [color="#0000ff"]################################################################ [color="#804040"] 15 [color="#0000ff"]#============= CLASS EIGENVALUEQR DEFINITION BEGIN ============# [color="#804040"] 16 [color="#0000ff"]################################################################ [color="#804040"] 17 [color="#804040"]class [color="#008080"]EigenvalueQR : [color="#804040"] 18 """[color="#ff00ff"]Compute the eigenvalue of a matrix with QR decomposition""" [color="#804040"] 19 [color="#804040"] 20 [color="#804040"]def [color="#008080"]__init__( self, matrix, eps = None ) : [color="#804040"] 21 """[color="#ff00ff"]Constructor of class""" [color="#804040"] 22 [color="#804040"] 23 [color="#0000ff"]# matrix is must be a squre matrix [color="#804040"] 24 [color="#804040"]if [color="#804040"]not matrix.isSquareMatrix() : [color="#804040"] 25 [color="#804040"]raise Exception, "[color="#ff00ff"]Square matrix is required" [color="#804040"] 26 [color="#804040"]return [color="#804040"] 27 [color="#804040"] 28 self.__n = matrix.getdimRow() [color="#804040"] 29 [color="#804040"] 30 [color="#0000ff"]# define the precision [color="#804040"] 31 [color="#804040"]if eps : [color="#804040"] 32 self.__eps = eps [color="#804040"] 33 [color="#804040"]else : [color="#804040"] 34 self.__eps = 0.0001 [color="#804040"] 35 [color="#804040"] 36 [color="#0000ff"]# Judge if matrix is a upper hessenberg matrix [color="#804040"] 37 [color="#0000ff"]# If not, implement the hessenberg transfromation [color="#804040"] 38 flag = self._isHessenbergMatrix( matrix ) [color="#804040"] 39 [color="#0000ff"]# is not a irreducible hessenberg matrix [color="#804040"] 40 [color="#804040"]if flag == 1 : [color="#804040"] 41 [color="#804040"]raise Exception, "[color="#ff00ff"]Must be irreducible Hessenberg Matrix" [color="#804040"] 42 [color="#804040"]return [color="#804040"] 43 [color="#0000ff"]# is not a upper hessenberg matrix [color="#804040"] 44 [color="#804040"]elif flag == 0 : [color="#804040"] 45 hb = Hessenberg( matrix ) [color="#804040"] 46 self.__H = hb.getH() [color="#804040"] 47 [color="#0000ff"]# is a irreducible upper hessenberg matrix [color="#804040"] 48 [color="#804040"]else : [color="#804040"] 49 self.__H = deepcopy( matrix ) [color="#804040"] 50 [color="#804040"] 51 [color="#0000ff"]# now the H is a irreducible upper hessenberg matrix [color="#804040"] 52 [color="#0000ff"]# we can implement QR decomposition [color="#804040"] 53 [color="#0000ff"]# control the condition for termination [color="#804040"] 54 iterator = self.__n - 1 [color="#804040"] 55 subH = self.__H [color="#804040"] 56 self.__eigenvalues = [] [color="#804040"] 57 [color="#804040"] 58 [color="#804040"]while iterator : [color="#804040"] 59 [color="#0000ff"]# one step QR mwthod [color="#804040"] 60 eigenv, subH = self._onestepQR( subH ) [color="#804040"] 61 self.__eigenvalues.append( eigenv ) [color="#804040"] 62 iterator -= 1 [color="#804040"] 63 [color="#804040"] 64 [color="#0000ff"]# get the last eigenvalue [color="#804040"] 65 self.__eigenvalues.append( subH[0][0] ) [color="#804040"] 66 [color="#804040"] 67 [color="#804040"] 68 [color="#0000ff"]#-----------------------------------------------------# [color="#804040"] 69 [color="#0000ff"]# Define utility methods of class [color="#804040"] 70 [color="#0000ff"]#-----------------------------------------------------# [color="#804040"] 71 [color="#804040"]def [color="#008080"]_isHessenbergMatrix( self, A ) : [color="#804040"] 72 """[color="#ff00ff"]Judge if A is a hessenberg matrix""" [color="#804040"] 73 [color="#804040"] 74 dim = A.getdimRow() [color="#804040"] 75 [color="#804040"] 76 [color="#804040"]for i [color="#804040"]in range( 2, dim ) : [color="#804040"] 77 [color="#804040"]for j [color="#804040"]in range( i-1 ) : [color="#804040"] 78 [color="#0000ff"]# must be hessenberg matrix [color="#804040"] 79 [color="#804040"]if A[j] != 0 : [color="#804040"] 80 [color="#804040"]return 0 [color="#804040"] 81 [color="#804040"] 82 [color="#804040"]for i [color="#804040"]in range( 1, dim ) : [color="#804040"] 83 [color="#0000ff"]# must be irreducible hessenberg matrix [color="#804040"] 84 [color="#804040"]if A[i-1] == 0 : [color="#804040"] 85 [color="#804040"]return 1 [color="#804040"] 86 [color="#804040"] 87 [color="#804040"]return 2 [color="#804040"] 88 [color="#804040"] 89 [color="#804040"]def [color="#008080"]_onestepQR( self, matrix ) : [color="#804040"] 90 """[color="#ff00ff"]Compute the upper hessenberg matrix's eigenvalue [color="#804040"] 91 [color="#ff00ff"] with onw step QR method""" [color="#804040"] 92 [color="#804040"] 93 [color="#0000ff"]#------------------------------------------------- [color="#804040"] 94 [color="#0000ff"]# Argorithm : [color="#804040"] 95 [color="#0000ff"]#------------------------------------------------- [color="#804040"] 96 [color="#0000ff"]# H[k] - s[k] * I = Q[k] * R[k] [color="#804040"] 97 [color="#0000ff"]# H[k+1] = R[k] * Q[k] + s[k] * I [color="#804040"] 98 [color="#0000ff"]#------------------------------------------------- [color="#804040"] 99 H = deepcopy(matrix) [color="#804040"]100 [color="#804040"]101 dim = H.getdimRow() [color="#804040"]102 n = dim - 1 [color="#804040"]103 [color="#804040"]104 [color="#804040"]while 1 : [color="#804040"]105 [color="#0000ff"]# construct H[k] - s[k] * I [color="#804040"]106 [color="#0000ff"]# get s[k] [color="#804040"]107 sk = H[n][n] [color="#804040"]108 offset = self._genDiagMatrix( dim, sk ) [color="#804040"]109 H -= offset [color="#804040"]110 [color="#804040"]111 [color="#0000ff"]# implement QR decomposition for H [color="#804040"]112 hh = Householder( H ) [color="#804040"]113 Q = hh.getQ() [color="#804040"]114 R = hh.getR() [color="#804040"]115 [color="#804040"]116 [color="#0000ff"]# update H : H[k+1] = R * Q + s[k] * I [color="#804040"]117 H = R * Q + offset [color="#804040"]118 [color="#804040"]119 [color="#0000ff"]# [color="#804040"]120 [color="#804040"]if H[n][n-1] |