社区 › 几何算法 › 实对称矩阵的特征值与特征向量求解

chlh_jd 发表于 2013-1-4 00:56:09

实对称矩阵的特征值与特征向量求解

本帖最后由 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 / mvv1 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 ( 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 ( (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)
i0);_Set the maximum number of iterations to prevent stack overflow
(while (and (not (equal d ( 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) ( 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( (length (car mat)))
i0 ;_Set the maximum number of iterations to prevent stack overflow
)
(while (and (not (equal d ( 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叠代法、行列互换法。

下面是特征值与特征向量求解的主要函数，请下载附件文件做测试。
;;;-------------------------------------------------------------------------------;;;
;;;                         Get EigenValues and EigenVectors                      ;;;
;;;-------------------------------------------------------------------------------;;;
(defun (mat / em i ai aii el vl res)
;; get Eigenvalues and Eigenvectors
;;
;; First Use QR-decomposition to get all EigenValues ,
;; Then Use Gaussian-elimination to get all EigenVectors .
;;
;; mat ---- a nonsingular matrix
;;
;; returns a list of all Eigenvalues and Eigenvectors , like '((Value1 Vector1) (Value2 Vector2) ... (Valuen Vectorn))
;;
;; e.g.
;; (setq a '((2. 1. 0.) (1. 3. 1.) (0. 1. 2.))) ;_( 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.))) ;_( 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.))) ;_( 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))
;;******if use Power method , this one will too slow .******
;;
;; (setq a '((1. 1. 0.5) (1. 1. 0.25) (0.5 0.25 2.))) ;_( a)--> ((2.53653 0.748221 0.649661 1.0) (1.48012 -0.63687 -0.805775 1.0) (-0.0166473 -7.69468 7.32278 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)))
;;          ;_( 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))
;;
;; (setq a '((2.0 3.0 8.0)(1.0 0.0 5.0) (-1.0 5.0 0.0))) ;_( a)--> ((5.0 4.58333 1.91667 1.0)(-4.19258 -0.807418 -1.0 1.0)(1.19258 -6.19258 -1.0 1.0))
;;******if use Power method , this one will too slow .******
;;
;; (setq a '((-3. 1. -1.) (-7. 5. -1.) (-6. 6. -2.))) ;_( a)--> ((4.0 0.0 1.0 1.0) (-2.0 -1.0 -1.0 0.0) (-2.0 -1.0 -1.0 0.0))
;;******if use Power method , this one will too slow .******
;;      
;;      get:   (setq vl ( a) v (nth 1 vl))
;;               
;;      checkEg: (equal (gc:determ ([-] ([λ] (length (car a)) (car v)) a)) 0.0 1e-8)
;;               
;;      checkEv(equal (mxv a (cdr v)) (mxv ([λ] (length a) (car v)) (cdr v)) 1e-8)
;;
;; by GSLS(SS) 2013-1-3
;;
;; Any suggestions are welcome , and please post them to chlh_jd@126.com or 56454300@qq.com
;;
;; Thanks Robert for method suggestions .
;; Thanks HighFlyingBird's help .
;; Thanks ElpanovEvgeniy's Gaussian-elimination function , it's really a pretty codes .;;----------- QR method------------
(setq Em (QR->EV mat 1e-12));_1e--10 The tolerance between the matrix modules with the product of it's Trace ,
                              ;_You can set a larger value like 1e-8 to reduce computations .
    (setq i 0)
    (while (and (setq ai (nth i Em)) (setq aii (nth i ai)))
    (setq el (cons aii el)
   i (1+ i)))                        
;;-----------------------------------

;;----------- Power method ----------
;_(setq el (Power->EV mat 1e-9))
;;-----------------------------------

(foreach ael;_Get all EigenVectors
    (setq vl (EigenValue->EigenVectors a mat))
    (if vl
      (foreach cvl
(setq res (cons (cons a c) res)))
      (setq res (cons (cons a nil) res))))

;;----------- QR method--------------
(mapcar (function(lambda (x)
       (mapcar (function (lambda (y)
      (if (equal y 0.0 1e-9)
      0.0
      y))) x)))
      res)
;;-----------------------------------;;----------- Power method ----------
;|(reverse (mapcar (function(lambda (x)
       (mapcar (function (lambda (y)
      (if (equal y 0.0 1e-9)
      0.0
      y))) x)))
      res))|;
;;-----------------------------------

)
;;
(defun EigenValue->EigenVectors(λ mat / n imat ires)
;_(setq mat a λ -2. ires nil)
;_(if (equal λ 4. 1e-4)(princ))
;_(princ (rtos λ 2 8))
(setq n    (length mat)
imat ( n)
mat( ([-] ([λ] n λ) mat))
mat(mapcar (function (lambda (x y)
   (append x y)))
       mat
       imat);_Creat matrix [(λI-A)^t | I]
mat(eea:gauss mat);_Gaussian elimination
imat nil)
(setq
    mat (mapcar (function (lambda (x / y)
       (setq x (reverse x))
       (repeat n
         (setq y (cons (car x) y)
      x (cdr x)))
       (setq imat (cons (reverse y) imat))
       x))
mat)
    )
(setq mat (reverse mat))
;; Check ZeroVector of Left Matrix .
(while (<= (apply (function +)
      (mapcar (function *) (car mat) (car mat)))
      5e-7)
    (setq mat(cdr mat)
   ires (cons (car imat) ires)
   imat (cdr imat)))
;; determine Linealy Independent
(while (and (cadr mat) (not (vslip (car mat) (cadr mat))))
    (setq
      ires (cons
      (vrefzero (car imat) (cadr imat) (car mat) (cadr mat))
      ires)
      imat (cdr imat)
      mat(cdr mat)))
;; Result deal ...
(mapcar
    (function reverse)
    (mapcar (function (lambda (b)
   (if (equal b 0.0 5e-7)
   0.0
   b)))
   ires))
)
;;;-------------------------------------高斯消元法--------------------------------;;;
;;;
(defun eea:gauss(lst)
;; Gauss eliminationby
;; by ElpanovEvgeniy
;; http://www.theswamp.org/index.php?action=post;quote=163485;topic=13505.13505.msg163485#msg163485
;; Implementation Gaussian elimination
;; For text:
;;1x+2y+3z=2
;;10x+1y+8z=17
;;7z+2y=5
;; (eea:gauss '((1.0 2.0 3.0) (10.0 1.0 8.0) (0.0 2.0 7.0)))
;; =>
;; ((1.0 2.0 3.0 2.0) (0.0 -19.0 -22.0 -3.0) (0.0 0.0 4.68421 4.68421))
;;(gauss lst)
(if (car lst)
    (if ;|(zerop (caar lst))|;(equal (caar lst) 0.0 1e-8) ;_here changed for
      (if (vl-every ;|(function zerop)|;
   (function (lambda (x) (equal x 0.0 1e-8))) ;_here changed for
   (mapcar (function car) lst))
(if (cdr lst)
   (cons
   (car lst)
   (mapcar
       (function (lambda (x) (cons 0. x)))
       (eea:gauss (mapcar (function cdr) (cdr lst)))
       ) ;_mapcar
   ) ;_cons
   lst
   ) ;_if
(eea:gauss
   (cons
   (mapcar
       (function +)
       (car lst)
       (car (vl-remove-if
       (function (lambda (x) (zerop (car x))))
       (cdr lst)))
       ) ;_mapcar
   (cdr lst)
   ) ;_cons
   ) ;_gauss
) ;_if
      (cons
(car lst)
(mapcar
   (function (lambda (x) (cons 0. x)))
   (eea:gauss
   (mapcar
       (function
(lambda (x / i)
    (if (and (equal (car x) 0.0 1e-6)
      (equal (caar lst) 0.0 1e-6)) ;_here add '0 / 0 = 1.'
      (setq i 1.)
      (setq i (/ (car x) (caar lst))))
    ;|(setq i (/ (car x) (caar lst)))|;
    (mapcar
      (function -)
      (cdr x)
      (mapcar (function (lambda (x1) (* x1 i)))
       (cdar lst))
      );_mapcar
    );_lambda
);_function
       (cdr lst)
       );_mapcar
   );_test
   );_mapcar
);_cons
      );_if
    );_if
);_defun
非常感谢Robert和HighFlyingBird斑竹的热心帮助，同时也致谢ElpanovEvgeniy的高斯消元函数、Gile的矩阵高斯求模函数。
希望大家多多指导
(setq a '((2. 1. 0.) (1. 3. 1.) (0. 1. 2.))) ( 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.)))( 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.)))( 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)))( 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))

qjchen 发表于 2022-6-10 17:29:35

landsat99 发表于 2022-6-8 18:46
eigen类的lsp实现，难得见到，赞。数值计算在commonlisp或scheme这类lisp中有什么优势特点吗？不了解lisp， ...
很久之前写过一些 矩阵求逆什么的 求方程解的，估计改改也能算特征值和特征向量
http://www.theswamp.org/index.php?topic=13505.msg163149#msg163149



但是可能是水平问题，就写得很辛苦，而autolisp的效率也比较低。所以个人觉得假如计算简单的有限元，几百根杆件以内的，是没有什么问题，非线性的可能就比较麻烦了。可能用objectarx会好些，3d3s就是这么做的。

不知道纯lisp会否好些，作为一个和fortran一样古老的语言活到现在，肯定有它厉害的地方。

帖子中的Evgeniy就是一个大高手，好像用common Lisp写过一个系统。

landsat99 发表于 2022-6-13 12:00:55

qjchen 发表于 2022-6-10 17:29
很久之前写过一些 矩阵求逆什么的 求方程解的，估计改改也能算特征值和特征向量
http://www.theswamp.or ...
资深大神！！ 深耕lsp十几年，大赞一个

这个矩阵运算的UI方式 太独特新颖了，创意很受启发。。。相当精彩，果断收藏！
lsp果然神一般的存在于CAD

landsat99 发表于 2022-6-8 18:46:03

eigen类的lsp实现，难得见到，赞。数值计算在commonlisp或scheme这类lisp中有什么优势特点吗？不了解lisp，，请大神科普下

chlh_jd 发表于 2013-1-4 19:56:54

本帖最后由 chlh_jd 于 2013-1-5 16:23 编辑

修改了下Jacobi函数，可以同时得到矩阵的特征值与特征向量，效率比QR法高了非常多；
使用的是依次对下三角元素消0叠代，每次迭代都必须进行下三角检查0检查，如果与0不在误差范围内则必须进行矩阵计算；
一般迭代收敛到1e-10需要n次迭代，n为方阵行数；每次迭代的最大矩阵乘法运算为1.5(n-1)^2次，但总的矩阵乘法次数有望在
(n-1)^2--->n=∞时为3(n-1)^2次以内；算法的复杂程度低于QR法（后续）；且可以改进Givens一次变换的元素数量，设置一个阀值，仅对于元素大于阀值的所有下三角元素进行矩阵运算，依次减小阀值，直到阀值足够小。

(defun Givens(mat Q / i j aij aii ajj aji ang a b m n0 n r Q1)
;;(Givens a)
;; Givens 变换求上三角阵
;; by GSLS(SS) 2013-1-4
(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 mm0
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) ( Q1)))
   )
)
      )
    )
(list mat Q))
;;-------------------------------------------------------------------------------;;;
(defun Jacobi(mat eps / d i b t1 q)
;; Jacobi method
;; by GSLS(SS) 2013-1-4
;;
;; (setq a '((2. 1. 0.) (1. 3. 1.) (0. 1. 2.)))--> ((1.0 0.57735 -0.57735 0.57735)(4.0 0.408248 0.816497 0.408248)(2.0 -0.707107 -2.59319e-010 0.707107))
;;
;; (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))) -->
;;      -->((0.0660769 0.738289 0.358251 -0.556068 -0.0643214 0.115054)(4.18334 -0.60114 0.59348 -0.468712 0.222282 -0.131577)(18.2577 0.284885 0.35799 0.425666 0.281067 -0.72835)
;;         (94.9829 0.100149 0.192449 0.29428 0.723881 0.585089) (7.50996 -0.0487036 0.595184 0.450894 -0.586047 0.310849))
;;
;; (Jacobi a 1e-12)
(setq t1 (getvar "millisecs")
d(gc:determ mat)
q( (length (car mat)))
i0 ;_Set the maximum number of iterations to prevent stack overflow
)
(while (and (not (equal d ( mat) eps)) (< i 1000))
    (setq b   (Givens mat Q)
   mat (car b)
   q   (cadr b))
    (setq i (1+ i)))
(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 (a b)
      (cons a b)))
   (mapcar (function (lambda (r1)
         (nth (setq i (1+ i)) r1)))
    mat)
   q)
)

chlh_jd 发表于 2013-1-5 16:22:06

本帖最后由 chlh_jd 于 2013-1-5 16:49 编辑

改进QR法函数代码，可以同时求出特征值与特征向量；与Jacobi法相比，要达到准确的特征向量解，需要把求模精度提高1倍，否则求到的特征向量会反号，尽管如此得到的特征向量的精度任然小于Jacobi法，这是因为Jacobi法采用的Givens平面旋转变换比Householder正交变换的数值计算来得稳定。Householder每运算一遍要作3n次矩阵乘法，即QR法每迭代一次需要做3n次矩阵乘法，对于N维方阵，QR法要迭代5~6N次才能得到稳定向量解，则QR法需要共需要15~18n^2次矩阵乘法；与Jacobi法相比，
15~18n^2/ =10~12n倍，如果Jacobi法临近下限3(n-2)^2次矩阵迭代，则比值变为5~6n^2。可见QR法的效率明显低于Jacobi法。
当然QR法也可以仅把方阵变成海森伯格(Hessenberg)上三角矩阵，这样来减少计算量：Ak=QkRk,Ak+1=RkQk。
;;;----------------------------------Householder变换------------------------------;;;
(defun householder(q mat / i m n d a b u hk)
;; 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))
)
(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 ( (1/ b) u)))
    (setq hk([-] imat u)
   q   ([*] q hk)
   mat ([*] ([*] hk mat) ( hk))
   i   (1+ i)))
(list q mat))
;;;------------------------------------QR分解迭代---------------------------------;;;
;;
(defun QR->EV(mat eps / imat q d i b t1)
;; QR-decomposition
;; by GSLS(SS) 2013-1-2
;; (setq a '((2. 1. 0.) (1. 3. 1.) (0. 1. 2.)))
;; (QR->EV a 5e-11)
(setq t1   (getvar "millisecs")
d    (gc:determ mat)
imat ( (length (car mat)))
q    imat
i    0) ;_Set the maximum number of iterations to prevent stack overflow
(while (and (not (equal d ( mat) eps)) (< i 1000))
    (setq b   (householder q mat)
   mat (cadr b)
   q   (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."))
(setq i -1)
(mapcar (function (lambda (a b)
      (cons a b)))
   (mapcar (function (lambda (r1)
         (nth (setq i (1+ i)) r1)))
    mat)
   ( q))
)

dwg001 发表于 2013-1-6 19:25:39

矩阵方面很好的资料，学习了。

自贡黄明儒 发表于 2013-1-8 08:01:58

东西很好，就是太深了。有没有简单实用的东西呢？

chlh_jd 发表于 2013-1-9 01:35:17

矩阵特征值与特征向量的应用在现代科学里面已经非常广泛了：图像问题，趋势问题，矩阵计算机...
朱小平老师有一本书叫《实对称矩阵的拟特征值理论与应用》（内容摘要）
这里有介绍2个应用例子：天气Markov链的稳定状态、人口流动问题；
这里介绍了3个应用例子：动物繁殖的规律问题、商品的市场占有率问题、常染色体遗传问题

但是在ACAD上具体应用，我还没发现；我最近一直在想，二维边界映射问题是否可以通过特征向量矩阵来求解，如果可以那么任意边界1范围内的图形就可以变换到任意边界2内了。

mj0000 发表于 2013-5-22 12:31:34

太深了，完全看不懂

pengfei2010 发表于 2013-6-22 09:32:37

矩阵很强大，，，

monkey_w 发表于 2018-3-23 09:27:06

向老师学习，感谢分享

zhenz02 发表于 2019-7-26 13:37:56

感谢楼主分享，收藏了

查看完整版本: 实对称矩阵的特征值与特征向量求解