本帖最后由 chlh_jd 于 2013-1-4 01:12 编辑
求非奇异实数矩阵A特征值与特征向量常用的求解方法有: 幂法,反幂法,QR迭代法, Givens变换收缩法,Jocobi法,同步求解法。
幂法用来求第一个特征值和特征向量,后面将给出第一次计算范例(详细步骤可在百度文库找到)。幂法收敛速度通常比较慢,如果提高计算速度就要损失精度,当前后两次迭代得到的特征值较为接近时,收敛很慢,因此还衍生出了加入瑞利(Rayleigh)商加速法、原位平移法等。幂法+收缩法就可以求全部特征值和特征向量,主要步骤为:
1.对矩阵A第一次求出特征值λ1和第1个特征向量x1;
2. 然后用0代替λ1所在行,将矩阵降阶,再次使用幂法,就可以得到第2个特征值λ2与特征向量x2 ;
3.依次类推,最终得到所有的特征值λ。
4. 对并列矩阵[ (λiI-A)^t | I ] 采用高斯消元法转为上三角矩阵,左阵0向量或相邻向量线性相关的向量对应的右阵向量即为特征向量,依次代入所有特征值求出所有特征向量。
- ;;;---------------------------------------幂法-----------------------------------;;;
- (defun PowerEigen (mat eps / v0 f0 vk f1)
- ;; 求矩阵第一特征值与特征向量求解
- ;; Get First EigenValue and it's EigenVector
- ;; mat 为方阵 ; eps 为精度、该精度不宜太小否收敛很慢, 宜为1e-6~1e-8
- ;; by GSLS(SS) 2012-12-24
- ;; (setq a '((2. 1. 0.) (1. 3. 1.) (0. 1. 2.))) -->(4. (0.5 1.0 0.5)) 即4为其中一个特征值,该特征值对应的特征向量为(0.5 1.0 0.5)
- ;; (setq a '((7. 3. -2.) (3. 4. 1.) (-2. -1. 3.))) --> (9.27492 (1.0 0.493399 -0.39736))
- ;; (setq a '((4. 1. 0.) (1. 0. -1.) (1. 1. -4.))) --> (4.20303 (1.0 0.20303 0.146657))
- ;; (setq a '((1. 1. 0.5) (1. 1. 0.25) (0.5 0.25 2.))) --> (2.53652587 (0.74822116 0.64966115 1.0))
- ;; (PowerEigen a 1e-9)
- (setq v0 (mapcar (function (lambda (x) 1.0)) (car mat))
- f0 1.0
- vk (mxv mat v0)
- f1 (apply (function max) (mapcar (function abs) vk))
- vk (pt* vk (1/ f1)))
- (while
- (not
- (and
- (<= (abs (- f1 f0)) eps)
- (<=
- ((lambda (x)
- (sqrt (apply (function +) (mapcar (function *) x x))))
- (pt- vk v0))
- eps)))
- (setq v0 vk
- f0 f1
- vk (mxv mat v0)
- f1 (apply (function max) (mapcar (function abs) vk))
- vk (pt* vk (1/ f1))))
- (cons f1 vk))
- ;;
- (defun Power->EV (mat eps / m v v1 eig res)
- ;; Get all EigenValues
- ;; by GSLS(SS) 2013-1-3
- ;; (PowerEigenValues a)
- (while (cadr mat)
- (setq eig (powereigen mat eps)
- res (cons eig res))
- (if res
- (setq v (car mat)
- v1 (cdar res)
- v1 (pt* v1 (1/ (car v1)))
- m (vtxv v1 v)
- mat ([-] mat m)
- mat (mapcar (function cdr) (cdr mat))
- ))
- )
- (reverse (mapcar (function car) (cons (car mat) res)))
- )
反幂法用来求最小特征值和特征向量,其迭代过程与幂法雷同。
QR迭代法是20世纪60年代产生的一种求中小规模的矩阵的特征值与特征向量的方法,先求得准确的矩阵模,将该模与叠代中间矩阵的主迹(脚标相同的元素积)做比较,如果在容差范围内,则停止迭代,否则进行下一步迭代。QR法的主要步骤:
R——矩阵A的上三角矩阵,A=QR
将A做Householder(镜像变换,由著名的数值分析专家Householder在1958年为讨论矩阵特征值问题而提出来的,在Grandyang的百度空间有详细的公式和算法说明http://hi.baidu.com/grandyang/item/e024400bf3c879ce915718e2)变换,得到Q1和R1,
新矩阵A1=R1Q1,如果A1的主迹等于|A|,停止迭代,否则继续
-----(镜像变换)---->Q2、R2
A2=R2Q2,如果A2的主迹等于|A|,停止迭代,否则继续
-----(镜像变换)---->Q3、R3
A3=R3Q3 ....
...
Ak=RkQk
...
An=RnQn
QR法经过N次迭代总能得到收敛,收敛速度依赖于相邻特征值的比值,为了加速收敛可以引入原点平移法,此时迭代公式改为Ak-sk*I=RkQk+sk*I , sk为位移量。求得的所有特征值后,采用幂法的第4步求特征向量。
QR法可以把求解精度设得很高,如1e-10,所需增加的计算量并不是太多(一个5x5方阵,计算精度由1e-8提高到1e-10其迭代次数也仅是由24次提高到28次),当然由于VLISP自身的计算精度并不是很高(如乘*、除/、幂expt等函数),精度也不适合设得太高,建议采用1e-10。值得一提的是QR法同时还可以用来求超秩矩阵(超定线性方程组的系数)的特征值。
- ;;;----------------------------------Householder变换------------------------------;;;
- (defun householder (mat / imat i m n d a b u hk q)
- ;; Use Householder transformation to Get QR matrix for QR-decomposition
- ;; (householder mat)
- ;; (setq mat '((2. 3.) (1. 2.) (5. 7.)))
- ;; (setq mat '((1. 2. -3.) (1. 3. 10.) (1. -1. 6.)))
- (setq i 0
- n (length (car mat))
- m (length mat)
- imat ([i] m)
- q imat)
- (while (< i n)
- (setq d (mapcar (function (lambda (x) (nth i x))) mat))
- (repeat i
- (setq d (cdr d)))
- (setq a (* (sign (car d))
- (sqrt (apply (function +)
- (mapcar (function (lambda (x)
- (* x x)))
- d))))
- b (* a (+ a (car d)))
- d (cons (+ (car d) a)
- (cdr d)))
- (repeat i
- (setq d (cons 0. d)))
- (setq u (vtxv d d))
- (if (not(zerop b))
- (setq u ([k*] (1/ b) u)))
- (setq hk ([-] imat u)
- q ([*] q hk)
- mat ([*] hk mat)
- i (1+ i)))
- (list q mat))
- ;;;------------------------------------QR分解迭代---------------------------------;;;
- ;;
- (defun QR->EV (mat eps / d i b t1)
- ;; QR-decomposition
- ;; by GSLS(SS) 2013-1-2
- (setq t1 (getvar "millisecs")
- d (gc:determ mat)
- i 0);_Set the maximum number of iterations to prevent stack overflow
- (while (and (not (equal d ([tr] mat) eps)) (< i 1000))
- (setq b (householder mat)
- mat ([*] (cadr b) (car b))
- i (1+ i)))
- (princ (strcat "\nQR-decomposition of iterations is : "
- (rtos i 2)
- " times ."))
- (princ (strcat "\nTotel time of QR-decomposition : "
- (rtos (- (getvar "millisecs") t1) 2 0)
- " m.s."))
- mat)
Givens变换(这是20世纪90年代的新算法,又称平面旋转变换) 是一种重要的正交交换, 它具有计算量小, 数值稳定的优点;欲将一个向量许多分量化为0采用Householder变换;欲将一个向量某个分量化为0采用Givens变换。Givens变换向量解和迭代法求解特征值和特征向量,主要步骤是:
1.先将矩阵A做右旋转相似变换,依次将从左到右的各列下三角元素变换为0,在这个过程中,下三角元素会逐步趋于0,因此需要多次迭代,将矩阵转换为上三角矩阵S;则各个对角线上的元素即为所有特征值。
2.之后可以用幂法的第4步来求解特征向量。
- ;;--------------------------------Givens变换法-----------------------------------;;;
- (defun Givens (mat Q / i j aij aii ajj aji ang a b m n0 n r Q1)
- ;;(Givens a)
- ;; Givens 变换求上三角阵
- (setq i (length mat)
- m0 i
- n0 (length (car mat))
- j -1)
- (while (< j (1- n0))
- (setq i m0
- j (1+ j))
- (while (> (setq i (1- i)) j)
- (setq aij (nth j (nth i mat)))
- (if (not (equal aij 0.0 1e-14))
- (progn
- (setq aii (nth i (nth i mat))
- ajj (nth j (nth j mat))
- aji (nth i (nth j mat)))
- (if (= aii ajj)
- (if (> aij 0)
- (setq ang _pi4)
- (if (< aij 0)
- (setq ang (- _pi4))))
- (setq ang (/ (atan (/ aij (- ajj aii) 0.5)) 2.))
- )
- (setq ang (max (min ang _pi4) (- _pi4)))
- (setq m m0
- Q1 nil)
- (while (<= 0 (setq m (1- m)))
- (setq n n0
- r nil)
- (while (<= 0 (setq n (1- n)))
- (setq r (cons (cond ((and (= m i) (= n j))
- (- (sin ang)))
- ((and (= m j) (= n i))
- (sin ang))
- ((and (= m i) (= n i))
- (cos ang))
- ((and (= m j) (= n j))
- (cos ang))
- ((= m n)
- 1.)
- (0.0))
- r)))
- (setq Q1 (cons r Q1))
- )
- (setq Q ([*] Q1 Q))
- (setq mat ([*] ([*] Q1 mat) ([trp] Q1)))
- )
- )
- )
- )
- (list mat Q))
- ;;-------------------
- (defun Givens->EV (mat eps / d i b t1 q )
- ;; Givens method main function
- ;; by GSLS(SS) 2013-1-4
- ;; (setq a '((4.0 2.0 6.0 8.0 2.0) (2.0 10.0 9.0 13.0 7.0) (6.0 9.0 14.0 20.0 12.0) (8.0 13.0 20.0 54.0 35.0) (2.0 7.0 12.0 35.0 43.0)))
- ;; (Givens->EV a 1e-12) --> (0.0660769 4.18334 18.2577 94.9829 7.50996)
- (setq t1 (getvar "millisecs")
- d (gc:determ mat)
- q ([i] (length (car mat)))
- i 0 ;_Set the maximum number of iterations to prevent stack overflow
- )
- (while (and (not (equal d ([tr] mat) eps)) (< i 1000))
- (setq b (Givens mat Q)
- mat (car b)
- q (cadr b))
- (setq i (1+ i)))
- (print mat)
- (princ (strcat "\nJacobi iterations is : "
- (rtos i 2)
- " times ."))
- (princ (strcat "\nTotel time of Jacobi solving : "
- (rtos (- (getvar "millisecs") t1) 2 0)
- " m.s."))
- (setq i -1)
- (mapcar (function (lambda (r1 )
- (nth (setq i (1+ i)) r1)))
- mat)
- )
仅求全部特征值而言,同精度下比较,Givens变换的迭代次数仅为QR分解迭代的1/5左右,总计算量约为QR法的一半,效率高了很多。 当然,Givens变换和Jacobi方法结合起来,就可以同时求出特征值和特征向量;但是这里编程起来比较麻烦,需要逐次迭代将从大到小得非对角线上元素采用Givens变换消0;并同时记录变换矩阵的乘积Q=Q*Qi,上三角矩阵S主对角线上的特征值对应的变换矩阵积Q的各列为特征向量;这种方法也可以算是同步求解法之一,据我所知,该法是目前应用最为广泛的方法。
同步求解法首先得列并阵,行式[ (λiI-A)^t | I ]和列式[ λiI-A / I ]都可以;对于行式采用初等行变换将左阵转换为上三角矩阵,主对角线上的元素即为特征值,特征所对应右阵的行即为特征向量;对于列式采用初等列变换,将上阵转换为下三角方阵,主对角线上的元素即为特征值,其对应下阵的各列即为特征向量。该法最为精准,但是对于大型矩阵的求解,要书写计算式就麻烦了。同步求解法,还有其他的方法来实现,如Jacobi叠代法、行列互换法。
下面是特征值与特征向量求解的主要函数,请下载附件文件做测试。
非常感谢Robert和HighFlyingBird斑竹的热心帮助,同时也致谢ElpanovEvgeniy的高斯消元函数、Gile的矩阵高斯求模函数。
希望大家多多指导
- (setq a '((2. 1. 0.) (1. 3. 1.) (0. 1. 2.))) ([E] a)--> ((4.0 1.0 2.0 1.0) (2.0 -1.0 0.0 1.0) (1.0 1.0 -1.0 1.0)) (setq a '((7. 3. -2.) (3. 4. 1.) (-2. -1. 3.)))([E] a)--> ((9.27492 -2.51661 -1.24169 1.0) (3.0 -1.0 2.0 1.0) (1.72508 2.51661 -3.75831 1.0)) (setq a '((4. 1. 0.) (1. 0. -1.) (1. 1. -4.)))([E] a)--> ((4.20303 6.81864 1.38439 1.0) (-3.7601 -0.0354877 0.275388 1.0) (-0.442931 -1.03315 4.59022 1.0)) (setq a '((4.0 2.0 6.0 8.0 2.0) (2.0 10.0 9.0 13.0 7.0) (6.0 9.0 14.0 20.0 12.0) (8.0 13.0 20.0 54.0 35.0) (2.0 7.0 12.0 35.0 43.0)))([E] a)--> ((94.9829 0.171168 0.328923 0.502966 1.23721 1.0) (18.2577 -0.391138 -0.491509 -0.584425 -0.385896 1.0) (7.50996 -0.156679 1.91471 1.45052 -1.88531 1.0)(4.18334 4.56872 -4.5105 3.56226 -1.68936 1.0) (0.0660769 6.41687 3.11375 -4.83309 -0.559052 1.0))
|