用maxima學數值分析-特徵值和特徵向量

高中曾學過用特徵方程式來求特徵值和特徵向量，但在實務上反而用power法、逆power法或QR法。
我會採用"線代啟示錄"網站的內容，並搭配maxima程式來介紹各種求特徵值和特徵向量的方法。
https://ccjou.wordpress.com

估計特徵值範圍的Gershgorin圓
https://ccjou.wordpress.com/2010 ... 的-gershgorin-圓/

A=−4−11032115，圓心(−40)，半徑0+1=1圓心(30)，半徑−1+1=2圓心(50)，半徑1+2=3
矩陣A和轉置矩陣AT有相同的特徵值
AT=−401−131125，圓心(−40)，半徑−1+1=2圓心(30)，半徑0+2=2圓心(50)，半徑1+1=2
矩陣A的特徵值為-4.144、2.390、5.754，都在Gershgorin圓內。


要先載入implicit_plot才能使用implicit_plot指令
(%i1)　load(implicit_plot);
(%o1)　C:\maxima-5.38.0\share\maxima\5.38.0_dirty\share\contrib\implicit_plot.lisp

矩陣A
(%i2)
A:matrix([-4,0,1],
              [-1,3,1],
              [1,2,5]);
(%o2)　−4−11032115

GershgorinCircle副程式
(%i3)
GershgorinCircle(A):=block
([circles:[],centers:"",xrange:[],yrange:[]],
for i:1 thru length(A) do
   (r:sum(if j#i then abs(A[i,j]) else 0,j,1,length(A[1])),
    circles:append(circles,[(x-A[i,i])^2+y^2=r^2]),
    centers:concat(centers,"set label ",i," at ",A[i,i],",0 point;"),
    xrange:append(xrange,[A[i,i]-r,A[i,i]+r]),
    yrange:append(yrange,[r,-r]),
    print("圓心(",A[i,i],"0),半徑",r)
  ),
  implicit_plot(circles,
                     [x,lmin(xrange)-1,lmax(xrange)+1],
                     [y,lmin(yrange)-1,lmax(yrange)+1],
                     same_xy,
                     [gnuplot_preamble,centers])
)$

計算矩陣A的Gershgorin Circle
(%i4)　GershgorinCircle(A);
圓心(−40),半徑1
圓心(30),半徑2
圓心(50),半徑3
(%o4)　done


計算矩陣A^t的Gershgorin Circle
(%i5)　GershgorinCircle(transpose(A));
圓心(−40),半徑2
圓心(30),半徑2
圓心(50),半徑2
(%o5)　done


[ 本帖最後由 bugmens 於 2016-7-28 09:37 AM 編輯 ]

Power迭代法
https://ccjou.wordpress.com/2010/11/02/power-遞迴法/

GooglePageRank就是用Power法來計算各網頁的PageRank值
https://math.pro/db/viewthread.php?tid=2299&page=1#pid13825

矩陣A=−4−11032115絕對值的特徵值由大到小排列57542−4144423902
－－－－－－－－－－－－－－－－－－－－－－－－－－－
虛擬碼
1.給定一初始向量u0
2.對於k=012，直至(k)收斂止，計算
3.uk+1=AukAuk
4.(k+1)=uTk+1Auk+1
－－－－－－－－－－－－－－－－－－－－－－－－－－－－
計算A=−4−11032115最大的絕對值特徵值和特徵向量
[計算過程]
設定初始向量u0=111

第1輪：
Au0=−4−11032115111=−338
特徵向量u1=Au0Au0=−338(−3)2+32+82=−033130331308835
特徵值1=uT1Au1=−033130331308835−4−11032115−033130331308835=41951

第2輪：
Au_1=\left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right] \left[ \matrix{-0.3313 \cr 0.3313 \cr 0.8835} \right]=\left[ \matrix{2.2087 \cr 2.2087 \cr 4.7488} \right]
特徵向量\displaystyle u_2=\frac{Au_1}{\Vert\;Au_1\Vert\;}=\frac{\left[ \matrix{2.2087 \cr 2.2087 \cr 4.7488}\right]}{\sqrt{2.2087^2+2.2087^2+4.7488^2}}=\left[ \matrix{0.3886 \cr 0.3886 \cr 0.8355} \right]
特徵值\lambda_2=u_2^T A u_2=\left[\matrix{0.3886&0.3886&0.8355}\right] \left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right] \left[\matrix{0.3886 \cr 0.3886 \cr 0.8355}\right]=4.8112

第3輪：
Au_2=\left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right] \left[ \matrix{0.3886 \cr 0.3886 \cr 0.8355} \right]=\left[ \matrix{-0.7189 \cr 1.6127 \cr 5.3433} \right]
特徵向量 \displaystyle u_3=\frac{Au_2}{\Vert\;Au_2\Vert\;}=\frac{\left[ \matrix{-0.7189 \cr 1.6127 \cr 5.3433}\right]}{\sqrt{(-0.7189)^2+1.6127^2+5.3433^2}}=\left[ \matrix{-0.1278 \cr 0.2866 \cr 0.9495} \right]
特徵值\lambda_3=u_3^T A u_3=\left[\matrix{-0.1278&0.2866&0.9495}\right] \left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right] \left[\matrix{-0.1278 \cr 0.2866 \cr 0.9495}\right]=5.2992

反覆計算直到最大絕對值特徵值\lambda=5.7542收斂為止


設定為數值運算
(%i1)　numer:true;
(%o1)　true

Power迭代法副程式
(%i2)
Power(A,epsilon):=block
([lambda1:0,lambda2:0,
   u:genmatrix(lambda([i,j],1),length(A),1)
],
do
  (norm:sqrt((A.u).(A.u)),
   u: (A.u)/norm,
   lambda1:transpose(u).A.u,
    print("特徵值",lambda1,",特徵向量",u),
    if abs(lambda1-lambda2)<epsilon then
      (return([lambda1,u])
      ),
    lambda2:lambda1
   )
)$

矩陣A
(%i3)
A:matrix([-4,0,1],
              [-1,3,1],
              [1,2,5]);
(%o3)　 \left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5}\right]

計算矩陣A的最大絕對值特徵值和特徵向量,誤差10^-5
(%i4)　Power(A,10^-5);
特徵值4.195121951219512,特徵向量 \left[ \matrix{-0.3312945782245396 \cr 0.3312945782245396 \cr 0.8834522085987723} \right]
特徵值4.811249528123821,特徵向量 \left[ \matrix{0.3885876702894004 \cr 0.3885876702894004 \cr 0.835463491122211} \right]
特徵值5.299190539203416,特徵向量 \left[ \matrix{-0.127751129016673 \cr 0.2865768569833476 \cr 0.9495016345833805} \right]
特徵值5.481375464225344,特徵向量 \left[ \matrix{0.2548162255461269 \cr 0.3379477601214591 \cr 0.9060132463862292} \right]
特徵值5.630841932875022,特徵向量 \left[ \matrix{-0.01983355591413924 \cr 0.2915954709929835 \cr 0.9563360870301721} \right]
特徵值5.675246645900037,特徵向量 \left[ \matrix{0.1801013767454748 \cr 0.3218782372414542 \cr 0.9294933536527943} \right]
特徵值5.723067597718835,特徵向量 \left[ \matrix{0.03644148428537539 \cr 0.2989084211151836 \cr 0.9535857454940834} \right]
特徵值5.730714847076354,特徵向量 \left[ \matrix{0.1403550033047277 \cr 0.3151515480071531 \cr 0.9386053349710008} \right]
特徵值5.747206056797989,特徵向量 \left[ \matrix{0.06562511526259228 \cr 0.3033809467502935 \cr 0.9506068300804844} \right]
特徵值5.7468154332089,特徵向量 \left[ \matrix{0.1195466125233259 \cr 0.3118719274730653 \cr 0.9425733437183723} \right]
特徵值5.753083430150924,特徵向量 \left[ \matrix{0.08074385491236928 \cr 0.3057786038135667 \cr 0.9486726913660559} \right]
特徵值5.751683756373387,特徵向量 \left[ \matrix{0.1087140210130338 \cr 0.310187221869381 \cr 0.9444390658079183} \right]
特徵值5.754312277364633,特徵向量 \left[ \matrix{0.08857834565554373 \cr 0.3070250780781634 \cr 0.9475703024641626} \right]
特徵值5.75326102960935,特徵向量 \left[ \matrix{0.1030864250836731 \cr 0.3093107842420867 \cr 0.9453570900538133} \right]
特徵值5.754455685847601,特徵向量 \left[ \matrix{0.09263978964425018 \cr 0.3076692887809196 \cr 0.9469727969248706} \right]
特徵值5.753822420970558,特徵向量 \left[ \matrix{0.1001652563476691 \cr 0.308854253924219 \cr 0.9458202637148985} \right]
特徵值5.754396491271145,特徵向量 \left[ \matrix{0.09474585840858986 \cr 0.3080022900102694 \cr 0.9466561211246931} \right]
特徵值5.754044259194526,特徵向量 \left[ \matrix{0.098649468141091 \cr 0.308616716056323 \cr 0.9460570833760988} \right]
特徵值5.754329734340657,特徵向量 \left[ \matrix{0.09583814464185311 \cr 0.3081746223192885 \cr 0.9464900697788494} \right]
特徵值5.754140620372628,特徵向量 \left[ \matrix{0.09786305388692572 \cr 0.3084932665854396 \cr 0.9461790143283508} \right]
特徵值5.754285386809937,特徵向量 \left[ \matrix{0.09640469653336128 \cr 0.3082638919733257 \cr 0.946403459097529} \right]
特徵值5.754185562769438,特徵向量 \left[ \matrix{0.09745508102955758 \cr 0.3084291606593783 \cr 0.9462420198006817} \right]
特徵值5.754259764059531,特徵向量 \left[ \matrix{0.09669857231455173 \cr 0.3083101629526002 \cr 0.9463584043756721} \right]
特徵值5.75420751750356,特徵向量 \left[ \matrix{0.09724344241992158 \cr 0.3083958871361926 \cr 0.9462746375676587} \right]
特徵值5.754245767596055,特徵向量 \left[ \matrix{0.09685101227284877 \cr 0.3083341551256965 \cr 0.9463349989325384} \right]
特徵值5.754218540433258,特徵向量 \left[ \matrix{0.09713365556656234 \cr 0.308378621525002 \cr 0.9462915400354248} \right]
特徵值5.754238317320978,特徵向量 \left[ \matrix{0.09693008733225696 \cr 0.308346597897462 \cr 0.9463228485748517} \right]
特徵值5.754224159959818,特徵向量 \left[ \matrix{0.0970767045344243 \cr 0.3083696637309211 \cr 0.9463003032479778} \right]
特徵值5.754234401503306,特徵向量 \left[ \matrix{0.09697110615386674 \cr 0.3083530516383455 \cr 0.946316543296489} \right]
特徵值5.754227048509759,特徵向量 \left[ \matrix{0.0970471617982491 \cr 0.3083650165891291 \cr 0.9463048477794515} \right]
(%o4)　[\; 5.754227048509759, \left[ \matrix{0.0970471617982491 \cr 0.3083650165891291 \cr 0.9463048477794515} \right] ]\;

[ 本帖最後由 bugmens 於 2016-8-4 07:07 AM 編輯 ]

Inverse Power迭代法

原本Power法求矩陣的最大絕對值特徵值，若將原來的矩陣改成反矩陣A^{-1}，可以求該矩陣的最小絕對值特徵值
矩陣A=\left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right]絕對值的特徵值由大到小排列|\;5.7542|\;>|\;-4.1444|\;>|\;2.3902|\;
－－－－－－－－－－－－－－－－－－－－－－－－－－－
虛擬碼
1.給定一初始向量u_0
2.對於k=0,1,2,\ldots，直至\lambda(k)收斂止，計算
3.\displaystyle u_{k+1}=\frac{A^{-1}u_k}{\Vert\; A^{-1}u_k \Vert\;}
4.\displaystyle \lambda(k+1)=u_{k+1}^T A u_{k+1}
－－－－－－－－－－－－－－－－－－－－－－－－－－－
計算A=\left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right]最小的絕對值特徵值和特徵向量
[計算過程]
設定初始向量 u_0=\left[ \matrix{1 \cr 1 \cr 1} \right]
反方陣 A^{-1}=\left[ \matrix{\displaystyle -\frac{13}{57}&-\frac{2}{57}&\frac{1}{19}\cr -\frac{2}{19}&\frac{7}{19}&-\frac{1}{19}\cr \frac{5}{57}&-\frac{8}{57}&\frac{4}{19}} \right]=\left[ \matrix{-0.2281&-0.0351&0.0526 \cr -0.1053&0.3684&-0.0526 \cr 0.0877&-0.1404&0.2105} \right]

第1輪：
A^{-1}u_0=\left[ \matrix{-0.2281&-0.0351&0.0526 \cr -0.1053&0.3684&-0.0526 \cr 0.0877&-0.1404&0.2105} \right] \left[ \matrix{1 \cr 1 \cr 1}\right]=\left[ \matrix{-0.2105 \cr 0.2105 \cr 0.1579} \right]
特徵向量 \displaystyle u_1=\frac{A^{-1}u_0}{\Vert\;A^{-1}u_0\Vert\;}=\frac{\left[ \matrix{-0.2105 \cr 0.2105 \cr 0.1579} \right]}{\sqrt{(-0.2105)^2+0.2105^2+0.1579^2}}=\left[ \matrix{-0.6247 \cr 0.6247 \cr 0.4685} \right]
特徵值 \displaystyle \lambda_1=u_1^T A u_1=[\; \matrix{-0.6247 & 0.6247 & 0.4685} ]\; \left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right] \left[ \matrix{-0.6247 \cr 0.6247 \cr 0.4685} \right]=1.3902

第2輪：
A^{-1}u_1=\left[ \matrix{-0.2281&-0.0351&0.0526 \cr -0.1053&0.3684&-0.0526 \cr 0.0877&-0.1404&0.2105} \right] \left[ \matrix{-0.6247 \cr 0.6247 \cr 0.4685} \right]=\left[ \matrix{0.1452 \cr 0.2712 \cr -0.0438} \right]
特徵向量\displaystyle u_2=\frac{A^{-1}u_1}{\Vert\; A^{-1}u_1 \Vert\;}=\frac{\left[ \matrix{0.1451 \cr 0.2712 \cr -0.0438} \right]}{\sqrt{0.1452^2+0.2712^2+(-0.0438)^2}}=\left[ \matrix{0.4673 \cr 0.8728 \cr -0.1411} \right]
特徵值 \displaystyle \lambda_2=u_2^T A u_2=[\; \matrix{0.4673 & 0.8728 & -0.1411} ]\; \left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right] \left[ \matrix{0.4673 \cr 0.8728 \cr -0.1411} \right]=0.6025

第3輪：
A^{-1}u_2=\left[ \matrix{-0.2281&-0.0351&0.0526 \cr -0.1053&0.3684&-0.0526 \cr 0.0877&-0.1404&0.2105} \right]\left[ \matrix{0.4673 \cr 0.8728 \cr -0.1411} \right]=\left[ \matrix{-0.1446 \cr 0.2798 \cr -0.1112} \right]
特徵向量\displaystyle u_3=\frac{A^{-1}u_2}{\Vert\; A^{-1}u_2 \Vert\;}=\frac{\left[ \matrix{-0.1446 \cr 0.2798 \cr -0.1112} \right]}{\sqrt{(-0.1446)^2+0.2798^2+(-0.1112)^2}}=\left[ \matrix{-0.4330 \cr 0.8377 \cr -0.3329} \right]
特徵值 \displaystyle \lambda_3=u_3^T A u_3=\left[ \matrix{-0.4330 & 0.8377 & -0.3329} \right] \left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right] \left[ \matrix{-0.4330 \cr 0.8377 \cr -0.3329} \right]=1.7238

反覆計算直到最小絕對值特徵值\lambda=2.3902收斂為止


設定為數值運算
(%i1)　numer:true;
(%o1)　true

InversePower迭代法副程式
(%i2)
InversePower(A,epsilon):=block
([lambda1:0,lambda2:0,
   u:genmatrix(lambda([i,j],1),length(A),1),
   invA:invert(A)
],
do
  (norm:sqrt((invA.u).(invA.u)),
   u: (invA.u)/norm,
   lambda1:transpose(u).A.u,
    print("特徵值",lambda1,",特徵向量",u),
    if abs(lambda1-lambda2)<epsilon then
      (return([lambda1,u])
      ),
    lambda2:lambda1
   )
)$

矩陣A
(%i3)
A:matrix([-4,0,1],
              [-1,3,1],
              [1,2,5]);
(%o3)　 \left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right]

計算矩陣A的最小絕對值特徵值和特徵向量,誤差10^-5
(%i4)　InversePower(A,10^-5);
特徵值1.390243902439024,特徵向量 \left[ \matrix{-0.6246950475544242 \cr 0.6246950475544242 \cr 0.4685212856658182} \right]
特徵值0.6025182651950873,特徵向量 \left[ \matrix{0.4672551492229502 \cr 0.8727973542089068 \cr -0.1410581582559849} \right]
特徵值1.723810898162325,特徵向量 \left[ \matrix{-0.4329600692855059 \cr 0.8376735458154866 \cr -0.3329393473970894} \right]
特徵值1.898404099984411,特徵向量 \left[ \matrix{0.1183560105491442 \cr 0.8488240781640499 \cr -0.5152567700630842} \right]
特徵值2.353621222544481,特徵向量 \left[ \matrix{-0.2088230074124524 \cr 0.8148882500740622 \cr -0.5406940830695776} \right]
特徵值2.290289620832044,特徵向量 \left[ \matrix{-0.02198003227500158 \cr 0.8178754811107383 \cr -0.5749752825809717} \right]
特徵值2.404381834227433,特徵向量 \left[ \matrix{-0.1304823266239602 \cr 0.8076149441908717 \cr -0.5750934396760006} \right]
特徵值2.366004441947428,特徵向量 \left[ \matrix{-0.06838457207541901 \cr 0.8096852683579864 \cr -0.5828664654156365} \right]
特徵值2.398226522212817,特徵向量 \left[ \matrix{-0.1043686777721781 \cr 0.8067637967691548 \cr -0.5815834895547717} \right]
特徵值2.383319198106618,特徵向量 \left[ \matrix{-0.0837054173929339 \cr 0.8077016380546078 \cr -0.5836192825661075} \right]
特徵值2.393280270604373,特徵向量 \left[ \matrix{-0.09565793262222155 \cr 0.8068381181023898 \cr -0.5829766814405501} \right]
特徵值2.388045766621258,特徵向量 \left[ \matrix{-0.08878099624356632 \cr 0.8071987795210824 \cr -0.5835649621470388} \right]
特徵值2.391257973753078,特徵向量 \left[ \matrix{-0.09275383808600854 \cr 0.8069324494868957 \cr -0.5833153070899075} \right]
特徵值2.389479539742537,特徵向量 \left[ \matrix{-0.09046555403219637 \cr 0.8070614399676388 \cr -0.5834962002027171} \right]
特徵值2.390533916277717,特徵向量 \left[ \matrix{-0.09178648747560744 \cr 0.8069766901476567 \cr -0.5834071153793234} \right]
特徵值2.389937075673304,特徵向量 \left[ \matrix{-0.09102519898238785 \cr 0.8070212218484542 \cr -0.5834647895429893} \right]
特徵值2.390285726718177,特徵向量 \left[ \matrix{-0.09146446943282166 \cr 0.8069937196208362 \cr -0.5834341328066945} \right]
特徵值2.390086420379639,特徵向量 \left[ \matrix{-0.0912112245795963 \cr 0.8070088200522667 \cr -0.5834528917303767} \right]
特徵值2.390202073842328,特徵向量 \left[ \matrix{-0.09135731484464309 \cr 0.8069997936976503 \cr -0.5834425198734888} \right]
特徵值2.390135660125875,特徵向量 \left[ \matrix{-0.09127307748649945 \cr 0.8070048669335477 \cr -0.5834486867512093} \right]
特徵值2.390174078419733,特徵向量 \left[ \matrix{-0.0913216657128404 \cr 0.8070018858881337 \cr -0.5834452069769944} \right]
特徵值2.390151969029395,特徵向量 \left[ \matrix{-0.09129364662248557 \cr 0.8070035821341262 \cr -0.5834472456949792} \right]
特徵值2.390164739323511,特徵向量 \left[ \matrix{-0.09130980700744958 \cr 0.8070025943037633 \cr -0.5834460831398715} \right]
特徵值2.390157382341876,特徵向量 \left[ \matrix{-0.09130048744691514 \cr 0.8070031600260253 \cr -0.5834467590962906} \right]
(%o4)　[2.390157382341876,\left[ \matrix{-0.09130048744691514 \cr 0.8070031600260253 \cr -0.5834467590962906} \right] ]

[ 本帖最後由 bugmens 於 2016-8-5 05:26 PM 編輯 ]

Shift Inverse Power迭代法

要求其他的特徵值，將原來的矩陣改成 (A-\alpha I)^{-1} ，可以求離\alpha最接近的特徵值。

－－－－－－－－－－－－－－－－－－－－－－－－－－－
虛擬碼
1.給定一初始向量u_0，偏移值\alpha
2.對於k=0,1,2,\ldots，直至\lambda(k)收斂止，計算
3.\displaystyle u_{k+1}=\frac{(A-\alpha I)^{-1}u_k}{\Vert\; (A-\alpha I)^{-1}u_k \Vert\;}
4.\displaystyle \lambda(k+1)=u_{k+1}^T A u_{k+1}
－－－－－－－－－－－－－－－－－－－－－－－－－－－
給定\alpha值，計算A=\left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right]離\alpha最接近的特徵值和特徵向量
[計算過程]
設定初始向量 u_0=\left[ \matrix{1 \cr 1 \cr 1} \right]
設定偏移值\alpha=-4(由Gershgorin Circle定理可知有個特徵值在圓心(-4,0)半徑1的圓內)
反方陣(A-\alpha I)^{-1}=\left[ \matrix{0&0&1 \cr -1&7&1 \cr 1&2&9} \right]^{-1}=\left[ \matrix{\displaystyle -\frac{61}{9}&-\frac{2}{9}&\frac{7}{9} \cr -\frac{10}{9}&\frac{1}{9}&\frac{1}{9} \cr 1&0&0} \right]=\left[ \matrix{-6.7778&-0.2222&0.7778 \cr -1.1111&0.1111&0.1111 \cr 1&0&0} \right]

第1輪：
(A-\alpha I)^{-1}u_0=\left[ \matrix{-6.7778&-0.2222&0.7778 \cr -1.1111&0.1111&0.1111 \cr 1&0&0} \right] \left[ \matrix{1 \cr 1 \cr 1} \right]=\left[ \matrix{-6.2222 \cr -0.8889 \cr 1} \right]
特徵向量\displaystyle u_1=\frac{(A-\alpha I)^{-1}u_0}{\Vert\; (A-\alpha I)^{-1}u_0 \Vert\;}=\frac{\left[ \matrix{-6.2222 \cr -0.8889 \cr 1} \right]}{\sqrt{(-6.2222)^2+(-0.8889)^2+1^2}}=\left[ \matrix{-0.9777 \cr -0.1397 \cr 0.1571} \right]
特徵值 \lambda_1=u_1^T A u_1=[\; \matrix{-0.9777 & -0.1397 & 0.1571} ]\; \left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right] \left[ \matrix{-0.9777 \cr -0.1397 \cr 0.1571} \right]=-4.1509

第2輪：
(A-\alpha I)^{-1}u_1=\left[ \matrix{-6.7778&-0.2222&0.7778 \cr -1.1111&0.1111&0.1111 \cr 1&0&0} \right] \left[ \matrix{-0.9777 \cr -0.1397 \cr 0.1571} \right]=\left[ \matrix{6.7796 \cr 1.0882 \cr -0.9777} \right]
特徵向量\displaystyle u_2=\frac{(A-\alpha I)^{-1}u_1}{\Vert\; (A-\alpha I)^{-1}u_1 \Vert\;}=\frac{\left[ \matrix{6.7796 \cr 1.0882 \cr -0.9777} \right]}{\sqrt{6.7796^2+1.0882^2+(-0.9777)^2}}=\left[ \matrix{0.9775 \cr 0.1569 \cr -0.1410} \right]
特徵值\lambda_2=u_2^T A u_2=\left[ \matrix{0.9775 & 0.1569 & -0.1410} \right] \left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right] \left[ \matrix{0.9775 \cr 0.1569 \cr -0.1410} \right]=-4.1441

第3輪：
(A-\alpha I)^{-1}u_2=\left[ \matrix{-6.7778&-0.2222&0.7778 \cr -1.1111&0.1111&0.1111 \cr 1&0&0} \right] \left[ \matrix{0.9775 \cr 0.1569 \cr -0.1410} \right]=\left[ \matrix{-6.7698 \cr -1.0843 \cr 0.9775} \right]
特徵向量\displaystyle u_3=\frac{\left[ \matrix{-6.7698 \cr -1.0843 \cr 0.9775} \right]}{\sqrt{(-6.7698)^2+(-1.0843)^2+0.9775^2}}=\left[ \matrix{-0.9775 \cr -0.1566 \cr 0.1411} \right]
特徵值\lambda_3=u_3^T A u_3=\left[ \matrix{-0.9775 & -0.1566 & 0.1411} \right] \left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right] \left[ \matrix{-0.9775 \cr -0.1566 \cr 0.1411} \right]=-4.1444

反覆計算直到特徵值\lambda=-4.1444收斂為止


設定為數值運算
(%i1)　numer:true;
(%o1)　true

ShiftInversePower迭代法副程式
(%i2)
ShiftInversePower(A,alpha,epsilon):=block
([lambda1:0,lambda2:0,
  u:genmatrix(lambda([i,j],1),length(A),1),
  ShiftInvA:invert(A-alpha*ident(length(A)))
],
do
  (norm:sqrt((ShiftInvA.u).(ShiftInvA.u)),
   u: (ShiftInvA.u)/norm,
   lambda1:transpose(u).A.u,
   print("特徵值",lambda1,",特徵向量",u),
   if abs(lambda1-lambda2)<epsilon then
     (return([lambda1,u])
     ),
     lambda2:lambda1
   )
)$

矩陣A
(%i3)
A:matrix([-4,0,1],
              [-1,3,1],
              [1,2,5]);
(%o3)　 \left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5}\right]

A[1,1]=-4,計算矩陣A離-4最近的特徵值,誤差10^-5
(%i4)　ShiftInversePower(A,A[1,1],10^-5);
特徵值-4.150868637610485,特徵向量 \left[ \matrix{-0.9776533929054586 \cr -0.1396647704150655 \cr 0.1571228667169487} \right]
特徵值-4.14414348452357,特徵向量 \left[ \matrix{0.9775024211666079 \cr 0.1569038221100049 \cr -0.1409617225373306} \right]
特徵值-4.144395503309593,特徵向量 \left[ \matrix{-0.9775285462559674 \cr -0.1565741860704146 \cr 0.1411469642290686} \right]
特徵值-4.144390102997756,特徵向量 \left[ \matrix{0.9775277105485285 \cr 0.1565807707392749 \cr -0.1411454474804167} \right]
(%o4)　 [\;-4.144390102997756,\left[ \matrix{0.9775277105485285 \cr 0.1565807707392749 \cr -0.1411454474804167} \right] ]\;

A[2,2]=3,計算矩陣A離3最近的特徵值,誤差10^-5
(%i5)　ShiftInversePower(A,A[2,2],10^-5);
特徵值3.4,特徵向量\left[ \matrix{0.0 \cr -0.447213595499958 \cr 0.8944271909999159} \right]
特徵值2.289558665231432,特徵向量\left[ \matrix{-0.06561787149247868 \cr 0.885841265148462 \cr -0.4593251004473507} \right]
特徵值2.426983132689005,特徵向量\left[ \matrix{0.0923650023965832 \cr -0.788287520453598 \cr 0.6083350157843926} \right]
特徵值2.382512315819198,特徵向量\left[ \matrix{-0.090698027506321 \cr 0.8111159042179215 \cr -0.5778103994661282} \right]
特徵值2.391869258878447,特徵向量\left[ \matrix{0.09140680193936326 \cr -0.8060912929597601 \cr 0.5846893397127062} \right]
特徵值2.389781526971133,特徵向量\left[ \matrix{-0.09127848203849248 \cr 0.8072048325603888 \cr -0.5831711558435508} \right]
特徵值2.390243815947408,特徵向量\left[ \matrix{0.09130930087416761 \cr -0.8069582576595488 \cr 0.5835074823589975} \right]
特徵值2.390141525165331,特徵向量\left[ \matrix{-0.09130268093133849 \cr 0.8070128481625035 \cr -0.583433015311436} \right]
特徵值2.390164183499367,特徵向量\left[ \matrix{0.09130416404587217 \cr -0.8070007600281967 \cr 0.5834495033349479} \right]
特徵值2.390159167449219,特徵向量\left[ \matrix{-0.09130383713123971 \cr 0.807003436501031 \cr -0.5834458525010168} \right]
(%o5)　 [\;2.390159167449219,\left[ \matrix{-0.09130383713123971 \cr 0.807003436501031 \cr -0.5834458525010168} \right] ]\;

A[3,3]=5,計算矩陣A離5最近的特徵值,誤差10^-5
(%i6)　ShiftInversePower(A,A[3,3],10^-5);
特徵值5.607038123167155,特徵向量 \left[ \matrix{0.05415303610738823 \cr 0.2166121444295529 \cr 0.9747546499329882} \right]
特徵值5.780069595037571,特徵向量 \left[ \matrix{0.09942088604805867 \cr 0.3337701174470541 \cr 0.9373969255959819} \right]
特徵值5.745800818656059,特徵向量 \left[ \matrix{0.09715673137405657 \cr 0.3011858672595754 \cr 0.9486029954157886} \right]
特徵值5.756608639911317,特徵向量 \left[ \matrix{0.09690361111870185 \cr 0.3104257059630139 \cr 0.9456455843652629} \right]
特徵值5.753536916223192,特徵向量 \left[ \matrix{0.09705322894473802 \cr 0.3077642101489163 \cr 0.9464997948773225} \right]
特徵值5.754430086933437,特徵向量 \left[ \matrix{0.09700387900903709 \cr 0.3085321381841197 \cr 0.9462548109070491} \right]
特徵值5.754172318128385,特徵向量 \left[ \matrix{0.09701866284758788 \cr 0.3083102704326569 \cr 0.946325607919921} \right]
特徵值5.754246849656881,特徵向量 \left[ \matrix{0.09701434887646511 \cr 0.3083743806347206 \cr 0.9463051608652613} \right]
特徵值5.754225312899597,特徵向量 \left[ \matrix{0.09701559914502492 \cr 0.3083558535212965 \cr 0.9463110699562194} \right]
特徵值5.754231537172271,特徵向量 \left[ \matrix{0.0970152375416282 \cr 0.3083612077029521 \cr 0.9463093623486551} \right]
(%o6)　 [\;5.754231537172271, \left[ \matrix{0.0970152375416282 \cr 0.3083612077029521 \cr 0.9463093623486551} \right] ]\;

[ 本帖最後由 bugmens 於 2016-8-11 04:22 PM 編輯 ]

ShiftInversePower迭代法+Rayleigh quotient

https://en.wikipedia.org/wiki/Rayleigh_quotient_iteration
前一篇(A-\alpha I)^{-1}的\alpha值是固定的，但Rayleigh quotient改成 \displaystyle \alpha=\frac{u^T A u}{u^T u}來加快收斂速度。


將特徵值問題的Au=\lambda u，(\lambda為特徵值，u為特徵向量)
代入 \displaystyle \alpha=\frac{u^T A u}{u^T u}=\frac{\lambda u^T u}{u^T u}=\lambda
也就是將(A-\alpha I)^{-1}替換成(A-\lambda I)^{-1}就可以加快收斂速度
－－－－－－－－－－－－－－－－－－－－－－－－－－－
虛擬碼
1.給定一初始向量u_0，初始特徵值\lambda_0
2.對於k=0,1,2,\ldots，直至\lambda(k)收斂止，計算
3.\displaystyle u_{k+1}=\frac{(A-\lambda_k I)^{-1}u_k}{\Vert\; (A-\lambda_k I)^{-1}u_k \Vert\;}
4.\lambda(k+1)=u_{k+1}^T A u_{k+1}
－－－－－－－－－－－－－－－－－－－－－－－－－－－
給定\lambda_0值，計算A=\left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right]離\lambda_0最接近的特徵值和特徵向量
[計算過程]
設定初始向量u_0=\left[ \matrix{1 \cr 1 \cr 1} \right]
設定初始特徵值\lambda_0=−4(由GershgorinCircle定理可知有個特徵值在圓心(−4,0)半徑1的圓內)

第1輪：
A-\lambda_0 I=\left[ \matrix{0&0&1 \cr -1&7&1 \cr 1&2&9} \right]
(A-\lambda_0 I)^{-1}u_0=\left[ \matrix{-6.7778&-0.2222&0.7778 \cr -1.1111&0.1111&0.1111 \cr 1&0&0} \right] \left[ \matrix{1 \cr 1 \cr 1} \right] = \left[ \matrix{-6.2222 \cr -0.8889 \cr 1} \right]
特徵向量 \displaystyle u_1=\frac{(A-\lambda_0 I)^{-1}u_0}{\Vert\; (A-\lambda_0 I)^{-1}u_0 \Vert\;}=\frac{\left[ \matrix{-6.2222 \cr -0.8889 \cr 1} \right]}{\sqrt{(-6.2222)^2+(-0.8889)^2+1^2}}=\left[ \matrix{-0.9777 \cr -0.1397 \cr 0.1571} \right]
特徵值 \lambda_1=u_1^T A u_1=\left[ \matrix{-0.9777&-0.1397&0.1571} \right] \left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right] \left[ \matrix{-0.9777 \cr -0.1397 \cr 0.1571} \right]=-4.1507

第2輪：
A-\lambda_1 I=\left[ \matrix{0.1509&0&1 \cr -1&7.1509&1 \cr 1&2&9.1509} \right]
(A-\lambda_1 I)^{-1}u_1=\left[ \matrix{151.1354&4.7649&-17.0367 \cr 24.1841&0.9067&-2.7419 \cr -21.8016&-0.7189&2.5703} \right] \left[ \matrix{-0.9777 \cr -0.1397 \cr 0.1571} \right]=\left[ \matrix{-151.1004 \cr -24.2011 \cr 21.8187} \right]
特徵向量 \displaystyle u_2=\frac{(A-\lambda_1 I)^{-1}u_1}{\Vert\; (A-\lambda_1 I)^{-1}u_1 \Vert\;}=\frac{\left[ \matrix{-151.1004 \cr -24.2011 \cr 21.8187} \right]}{\sqrt{(-151.1004)^2+(-24.2011)^2+21.8187^2}}=\left[ \matrix{-0.9775 \cr -0.1566 \cr 0.1412} \right]
特徵值 \lambda_2=u_1^T A u_1=\left[ \matrix{-0.9775 & -0.1566 & 0.1412} \right] \left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right] \left[ \matrix{-0.9775 \cr -0.1566 \cr 0.1412} \right]=-4.1444

第3輪：
A-\lambda_2 I=\left[ \matrix{0.1444&0&1 \cr -1&7.1444&1 \cr 1&2&9.1444} \right]
(A-\lambda_2 I)^{-1}u_2=\left[ \matrix{5.7856 \cdot 10^{10} & 1.8271 \cdot 10^9 & -6.5267 \cdot 10^9 \cr 9.2674 \cdot 10^9 & 2.9296 \cdot 10^8 & -1.0455 \cdot 10^9 \cr -8.3538 \cdot 10^9 & -2.6381 \cdot 10^8 & 9.7239 \cdot 10^8} \right] \left[ \matrix{-0.9776 \cr -0.1566 \cr 0.1411} \right]=\left[ \matrix{-5.7762 \cdot 10^{10} \cr -9.2525 \cdot 10^9 \cr 8.3404 \cdot 10^9}\right]
特徵向量 \displaystyle u_3=\frac{(A-\lambda_2 I)^{-1}u_2}{\Vert\; (A-\lambda_2 I)^{-1}u_2 \Vert\;}=\frac{\left[ \matrix{-5.7763 \cdot 10^{10} \cr -9.2525 \cdot 10^9 \cr 8.3404 \cdot 10^9} \right]}{\sqrt{(-5.7763 \cdot 10^{10})^2+(-9.2525 \cdot 10^9)^2+(8.3404 \cdot 10^9)^2}}=\left[ \matrix{-0.9775 \cr -0.1566 \cr 0.1411} \right]
特徵值 \lambda_3=u_2^T A u_2=\left[ \matrix{-0.9775 & -0.1566 & 0.1411} \right] \left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right] \left[ \matrix{-0.9775 \cr -0.1566 \cr 0.1411} \right]=-4.1444

反覆計算直到特徵值\lambda=-4.1444收斂為止


設定為數值運算
(%i1)　numer:true;
(%o1)　true

RayleighQuotientPower迭代法副程式
(%i2)
RayleighQuotientPower(A,lambda1,epsilon):=block
([lambda2:0,
  u:genmatrix(lambda([i,j],1),length(A),1)
],
do
  (y:invert(A-lambda1*ident(length(A))).u,
   norm:sqrt(y.y),
   u:y/norm,
   lambda1:transpose(u).A.u,
   print("特徵值",lambda1,",特徵向量",u),
   if abs(lambda1-lambda2)<epsilon then
     (return([lambda1,u])
     ),
   lambda2:lambda1
  )
)$

矩陣A
(%i3)
A:matrix([-4,0,1],
              [-1,3,1],
              [1,2,5]);
(%o3)　 \left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&2}\right]

A[1,1]=-4,計算矩陣A離-4最近的特徵值,誤差10^-5
(%i4)　RayleighQuotientPower(A,A[1,1],10^-5);
特徵值-4.150868637610485,特徵向量 \left[ \matrix{-0.9776533929054586 \cr -0.1396647704150655 \cr 0.1571228667169487} \right]
特徵值-4.144400989358965,特徵向量 \left[ \matrix{-0.9775288279606554 \cr -0.1565664321230263 \cr 0.141153614328975} \right]
特徵值-4.144390218552043,特徵向量 \left[ \matrix{-0.9775277325400289 \cr -0.1565806375486704 \cr 0.1411454429303385} \right]
特徵值-4.144390218535118,特徵向量 \left[ \matrix{-0.9775277325374503 \cr -0.1565806375693802 \cr 0.1411454429252224} \right]
(%o4)　 [\;-4.144390218535118,\left[ \matrix{-0.9775277325374503 \cr -0.1565806375693802 \cr 0.1411454429252224} \right] ]\;

A[2,2]=3,計算矩陣A離3最近的特徵值,誤差10^-5
(%i5)　RayleighQuotientPower(A,A[2,2],10^-5);
特徵值3.4,特徵向量 \left[ \matrix{0.0 \cr -0.447213595499958 \cr 0.8944271909999159} \right]
特徵值2.34480387797816,特徵向量 \left[ \matrix{-0.04418147155945099 \cr 0.9440107756536034 \cr -0.3269428895399374} \right]
特徵值2.386261415958715,特徵向量 \left[ \matrix{-0.09118695120399674 \cr 0.8091554911793659 \cr -0.5804759521499571} \right]
特徵值2.390155311149527,特徵向量 \left[ \matrix{-0.09130350144505062 \cr 0.8070054621306113 \cr -0.5834431032373528} \right]
特徵值2.390160076833697,特徵向量 \left[ \matrix{-0.09130389650666068 \cr 0.8070029513063456 \cr -0.5834465143143361} \right]
(%o5)　 [\; 2.390160076833697,\left[ \matrix{-0.09130389650666068 \cr 0.8070029513063456 \cr -0.5834465143143361} \right] ]\;

A[3,3]=5,計算矩陣A離5最近的特徵值,誤差10^-5
(%i6)　RayleighQuotientPower(A,A[3,3],10^-5);
特徵值5.607038123167155,特徵向量 \left[ \matrix{0.05415303610738823 \cr 0.2166121444295529 \cr 0.9747546499329882} \right]
特徵值5.75863661416414,特徵向量 \left[ \matrix{0.09749930430177767 \cr 0.3123865901727191 \cr 0.9449383598631875} \right]
特徵值5.754236049858843,特徵向量 \left[ \matrix{-0.09701530629188788 \cr -0.308365169743175 \cr -0.9463080642340285} \right]
特徵值5.754230141705235,特徵向量 \left[ \matrix{0.09701531863069972 \cr 0.3083600072972683 \cr 0.9463097451947816} \right]
(%o6)　 [\; 5.754230141705235,\left[ \matrix{0.09701531863069972 \cr 0.3083600072972683 \cr 0.9463097451947816} \right] ]\;



比較兩種方法所需的計算次數，發現採用Rayleigh quotient計算次數少於Shift Invers Power。

	\lambda=-4.1444	\lambda=2.3902	\lambda=5.7542
Shift Inverse Power 迭代法	4次	10次	10次
Rayleigh quotient 迭代法	4次	5次	4次

power法一次只能求一個特徵值，Rutishauser在1958年提出LR演算法，可以求出全部的特徵值。
LR演算法是基於LU分解，也就是將A分解成L和U兩個矩陣，其中L為下三角矩陣，U為上三角矩陣，再將U乘上L形成下一個矩陣A，如此反覆計算最後的矩陣A的對角線就是全部的特徵值。
關於LU分解請見https://ccjou.wordpress.com/2010/09/01/lu-分解/
－－－－－－－－－－－－－－－－－－－－－－－－－－－
虛擬碼
1.對於k=0,1,2,\ldots，直到 \lambda(1),\lambda(2),\ldots,\lambda(n) 收斂為止，計算
2.A_k=L_k U_k，計算A_k的LU分解。
3.A_{k+1}=U_k L_k，將L_k和U_k對調後相乘得到下一個矩陣A_{k+1}
4.特徵值 \lambda(1),\lambda(2),\ldots,\lambda(n) 在矩陣A_{k+1}的對角線上。
－－－－－－－－－－－－－－－－－－－－－－－－－－－
計算A=\left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right] 全部的特徵值。
[計算過程]
第1輪：
計算矩陣A的LU分解
\left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right]=\left[ \matrix{1&0&0 \cr 0.25&1&0 \cr -0.25&0.6667&1} \right] \left[ \matrix{-4&0&1 \cr 0&3.0&0.75 \cr 0&0&4.75} \right]
將L和U對調後相乘得到下一個矩陣A
\left[ \matrix{-4&0&1 \cr 0&3.0&0.75 \cr 0&0&4.75} \right] \left[ \matrix{1&0&0 \cr 0.25&1&0 \cr -0.25&0.6667&1} \right]=\left[ \matrix{-4.25&0.6667&1 \cr 0.5625&3.5&0.75 \cr -1.1875&3.1667&4.75} \right]
特徵值 \left[-4.25,3.5,4.75 \right]

第2輪：
計算矩陣A的LU分解
\left[ \matrix{-4.25&0.6667&1 \cr 0.5625&3.5&0.75 \cr -1.1875&3.1667&4.75} \right]=\left[ \matrix{1&0&0 \cr -0.1324&1&0 \cr 0.2794&0.8306&1} \right] \left[ \matrix{-4.25&0.6667&1 \cr 0&3.5882&0.8824 \cr 0&0&3.7377} \right]
將L和U對調後相乘得到下一個矩陣A
\left[ \matrix{-4.25&0.6667&1 \cr 0&3.5882&0.8824 \cr 0&0&3.7377} \right] \left[ \matrix{1&0&0 \cr -0.1324&1&0 \cr 0.2794&0.8306&1} \right]=\left[ \matrix{-4.0588&1.4973&1.0 \cr -0.2284&4.3211&0.8824 \cr 1.0444&3.1045&3.7377} \right]
特徵值 \left[ -4.0588,4.3211,3.7377 \right]

第3輪：
計算矩陣A的LU分解
\left[ \matrix{-4.0588&1.4973&1.0 \cr -0.2284&4.3211&0.8824 \cr 1.0444&3.1045&3.7377} \right]=\left[ \matrix{1&0&0 \cr 0.0563&1&0 \cr -0.2573&0.8237&1} \right] \left[ \matrix{-4.0588&1.4973&1.0 \cr 0&4.2369&0.8261 \cr 0&0&3.3146} \right]
將L和U對調後相乘得到下一個矩陣A
\left[ \matrix{-4.0588&1.4973&1.0 \cr 0&4.2369&0.8261 \cr 0&0&3.3146} \right] \left[ \matrix{1&0&0 \cr 0.0563&1&0 \cr -0.2573&0.8237&1} \right]=\left[ \matrix{-4.2319&2.3209&1.0 \cr 0.0253&4.9173&0.8261 \cr -0.8529&2.7301&3.3145} \right]
特徵值 \left[ -4.2319 , 4.9173 , 3.3146 \right]

反覆計算直到全部特徵值 \left[ 5.7542,-4.1444,2.3902 \right]


設定為數值運算
(%i1)　numer:true;
(%o1)　true

不顯示小數轉換成分數的過程
例如rat: replaced -4.25 by -17/4 = -4.25
(%i2)　ratprint:false;
(%o2)　false

LR迭代法副程式
/*print(...);*/將兩邊的/*和*/刪除後可以顯示計算過程
但顯示太多計算過程maxima程式會當掉,只好註解起來不顯示
(%i3)
LR(A,epsilon):=block
([count:0,P,L,U,n:length(A),eig1,eig2],
eig2:makelist(0,n),
do
   (count:count+1,
    /*print("第",count,"輪"),*/
    [P,L,U]:get_lu_factors(lu_factor(A)),
    /*print("計算矩陣A的LU分解",A,"="),*/
    /*print(L,"*",U),*/
    A:U.L,
    /*print("將L和U對調後相乘得到下一個矩陣A"),*/
    /*print(U,"*",L,"=",A),*/
    eig1:create_list(A[i,i],i,1,n),
    print("特徵值",eig1),
    if lmax(abs(eig1-eig2))<epsilon then
      (return(eig1)),
    eig2:eig1
  )
)$

矩陣A
(%i4)
A:matrix([-4,0,1],
              [-1,3,1],
              [1,2,5]);
(%o4)　 \left[ \matrix{-4&0&1\cr -1&3&1\cr 1&2&5} \right]

計算矩陣A全部的特徵值,誤差10^-5
(%i5)　LR(A,10^-5);
0 errors, 0 warnings
特徵值 \left[-4.25,3.5,4.75\right]
特徵值 \left[-4.058823529411764,4.321118611378978,3.737704918032787\right]
特徵值 \left[-4.231884057971015,4.91729873861925,3.314585319351763\right]
特徵值 \left[-4.044520547945205,5.313248850521489,2.731271697423716\right]
特徵值 \left[-4.27519051651143,5.669125227250055,2.606065289261376\right]
特徵值 \left[-3.97464844523668,5.528967112469076,2.445681332767604\right]
特徵值 \left[-4.386336456049431,5.951388036215438,2.434948419833993\right]
特徵值 \left[-3.82905992615734,5.432559904220157,2.396500021937184\right]
特徵值 \left[-4.616971218176369,6.21703651117396,2.39993470700241\right]
特徵值 \left[-3.556070568397674,5.165962621376724,2.39010794702095\right]
特徵值 \left[-5.09592634683869,6.703414031205933,2.392512315632756\right]
特徵值 \left[-3.07014305149453,4.680374938114763,2.389768113379766\right]
特徵值 \left[-6.15763906410022,7.766852181380823,2.390786882719396\right]
特徵值 \left[-2.263100148442174,3.873135110765942,2.389965037676232\right]
特徵值 \left[-8.927750729051786,10.53740902683053,2.390341702221251\right]
特徵值 \left[-1.061374420236201,2.671290372000583,2.390084048235619\right]
特徵值 \left[-20.85888707850496,22.46867122537871,2.390215853126253\right]
特徵值 \left[0.4665444106715815,1.143322725027585,2.390132864300833\right]
特徵值 \left[52.72558253627871,-51.11576036528454,2.390177829005824\right]
特徵值 \left[2.06213873325895,-0.4522894269100236,2.390150693651076\right]
特徵值 \left[13.17442195234688,-11.56458779576555,2.390165843418669\right]
特徵值 \left[3.419996923785756,-1.810153822408377,2.390156898622622\right]
特徵值 \left[8.58287969592947,-6.97304166698185,2.39016197105238\right]
特徵值 \left[4.388368293639632,-2.77852730349426,2.39015900985463\right]
特徵值 \left[7.044154778522382,-5.434315481286824,2.390160702764444\right]
特徵值 \left[4.995309976952969,-3.385469697199792,2.390159720246825\right]
特徵值 \left[6.383873007186085,-4.774033291505344,2.390160284319262\right]
特徵值 \left[5.34546797814534,-3.735627936085535,2.390159957940198\right]
特徵值 \left[6.071147185395128,-4.46130733112236,2.390160145727236\right]
特徵值 \left[5.537890814471758,-3.928050851714158,2.390160037242404\right]
特徵值 \left[5.916131959742552,-4.306292059474105,2.390160099731557\right]
特徵值 \left[5.640814092845454,-4.030974156506412,2.390160063660963\right]
特徵值 \left[5.8375586037984,-4.227718688248872,2.390160084450476\right]
特徵值 \left[5.695070878935839,-4.08523095139107,2.390160072455235\right]
特徵值 \left[5.797281242393553,-4.187441321764426,2.390160079370877\right]
特徵值 \left[5.723453549499145,-4.113613624880672,2.390160075381531\right]
特徵值 \left[5.776515673768994,-4.166675751450871,2.390160077681882\right]
特徵值 \left[5.738241274640624,-4.128401350995675,2.390160076355056\right]
特徵值 \left[5.765777948950188,-4.15593802607038,2.390160077120198\right]
特徵值 \left[5.745929679925635,-4.136089756604526,2.390160076678896\right]
特徵值 \left[5.760217049729156,-4.150377126662543,2.390160076933392\right]
特徵值 \left[5.749922650695337,-4.140082727481946,2.390160076786614\right]
特徵值 \left[5.757334865818812,-4.14749494269007,2.390160076871262\right]
特徵值 \left[5.751995220747789,-4.142155297570226,2.390160076822442\right]
特徵值 \left[5.7558404324125,-4.146000509263092,2.390160076850597\right]
特徵值 \left[5.753070680591324,-4.143230757425678,2.39016007683436\right]
特徵值 \left[5.755065392906378,-4.145225469750097,2.390160076843724\right]
特徵值 \left[5.753628652968832,-4.14378872980715,2.390160076838324\right]
特徵值 \left[5.754663399407986,-4.144823476249418,2.390160076841438\right]
特徵值 \left[5.753918118401036,-4.144078195240673,2.390160076839642\right]
特徵值 \left[5.754454883559434,-4.144614960400106,2.390160076840678\right]
特徵值 \left[5.75406828136153,-4.144228358201605,2.390160076840081\right]
特徵值 \left[5.754346722233892,-4.144506799074311,2.390160076840425\right]
特徵值 \left[5.754146178167928,-4.144306255008148,2.390160076840227\right]
特徵值 \left[5.754290615935278,-4.144450692775613,2.390160076840341\right]
特徵值 \left[5.754186586569208,-4.144346663409477,2.390160076840275\right]
特徵值 \left[5.754261511799679,-4.144421588639986,2.390160076840314\right]
特徵值 \left[5.754207548012657,-4.144367624852942,2.390160076840291\right]
特徵值 \left[5.754246414489054,-4.144406491329351,2.390160076840304\right]
特徵值 \left[5.75421842151387,-4.144378498354161,2.390160076840297\right]
特徵值 \left[5.754238582980808,-4.144398659821102,2.390160076840301\right]
特徵值 \left[5.75422406200547,-4.144384138845762,2.390160076840299\right]
特徵值 \left[5.75423452049635,-4.144394597336643,2.3901600768403\right]
特徵值 \left[5.75422698793699,-4.144387064777282,2.390160076840299\right]
(%o5)　 \left[5.75422698793699,-4.144387064777282,2.390160076840299\right]


參考資料
https://books.google.com.tw/book ... e&q&f=false

[ 本帖最後由 bugmens 於 2016-9-16 10:41 PM 編輯 ]

https://ccjou.wordpress.com/2010/09/01/lu-分解/
這篇利用maxima來計算LU分解的計算過程，可以觀察到其實是高斯消去法，L矩陣紀錄消去法化簡A矩陣的過程，而U矩陣則儲存化簡的結果。

這個程式無法處理對角線元素為0的情況，詳細內容請見"PA=LU分解"。
https://ccjou.wordpress.com/2012/04/13/palu-分解/



LU分解副程式
(%i1)
LU(A):=block
([n:length(A),U,L],
print("初始值L=",L:ident(n),",U=",U:A),
print("------------------"),
for i:1 thru n do
   (for j:i+1 thru n do
      (L[j,i]:U[j,i]/U[i,i],
       print("設L[",j,",",i,"]=U[",j,",",i,"]/U[",i,",",i,"]=",L[j,i],",L=",L),
       print("原本U=",U,",第",j,"列-",L[j,i],"*第",i,"列,得到U=",U:rowop(U,j,i,L[j,i])),
       print("------------------")
      )
   ),
return([L,U])
)$

本文使用的範例
(%i2)
A:matrix([-4,0,1],
              [-1,3,1],
              [1,2,5]);
(%o2)　 \left[ \matrix{-4&0&1\cr -1&3&1\cr 1&2&5} \right]

執行LU分解
(%i3)　[L,U]: LU(A);
初始值L=\left[ \matrix{1&0&0\cr 0&1&0\cr 0&0&1} \right],U=\left[ \matrix{-4&0&1\cr -1&3&1\cr 1&2&5} \right]
------------------
設 \displaystyle L [\; 2,1 ]\;=U[\; 2,1 ]\; / U[\; 1,1 ]\;=\frac{1}{4} ,L=\left[ \matrix{\displaystyle 1&0&0\cr \frac{1}{4}&1&0\cr 0&0&1} \right]
0 errors, 0 warnings
原本U=\left[ \matrix{-4&0&1\cr -1&3&1\cr 1&2&5} \right],第2列 \displaystyle -\frac{1}{4}* 第1列,得到U=\left[ \matrix{\displaystyle -4&0&1\cr 0&3&\frac{3}{4} \cr 1&2&5} \right]
------------------
設 \displaystyle L [\; 3,1 ]\;=U[\; 3,1 ]\; / U[\; 1,1 ]\;=-\frac{1}{4} ,L=\left[ \matrix{\displaystyle 1&0&0\cr \frac{1}{4}&1&0\cr -\frac{1}{4}&0&1} \right]
原本U=\left[ \matrix{\displaystyle -4&0&1\cr 0&3&\frac{3}{4}\cr 1&2&5} \right],第3列 \displaystyle --\frac{1}{4}* 第1列,得到U=\left[ \matrix{\displaystyle -4&0&1\cr 0&3&\frac{3}{4}\cr 0&2&\frac{21}{4}} \right]
------------------
設 \displaystyle L [\; 3,2 ]\;=U[\; 3,2 ]\; / U[\; 2,2 ]\;=\frac{2}{3} ,L=\left[ \matrix{\displaystyle 1&0&0 \cr \frac{1}{4}&1&0\cr -\frac{1}{4}&\frac{2}{3}&1} \right]
原本U=\left[ \matrix{\displaystyle -4&0&1\cr 0&3&\frac{3}{4}\cr 0&2&\frac{21}{4}} \right],第3列 \displaystyle -\frac{2}{3}* 第2列,得到U=\left[ \matrix{\displaystyle -4&0&1\cr 0&3&\frac{3}{4}\cr 0&0&\frac{19}{4}} \right]
------------------
(%o3)　 [\; \left[ \matrix{\displaystyle 1&0&0 \cr \frac{1}{4}&1&0\cr -\frac{1}{4}&\frac{2}{3}&1} \right],\left[ \matrix{\displaystyle -4&0&1 \cr 0&3&\frac{3}{4} \cr 0&0&\frac{19}{4}} \right] ]\;

利用maxima內建指令得到相同的答案
(%i4)　[P,L,U]:get_lu_factors(lu_factor(A));
(%o4)　 [\; \left[ \matrix{1&0&0\cr 0&1&0\cr 0&0&1} \right],\left[ \matrix{\displaystyle 1&0&0\cr \frac{1}{4}&1&0\cr -\frac{1}{4}&\frac{2}{3}&1} \right],\left[ \matrix{\displaystyle -4&0&1\cr 0&3&\frac{3}{4}\cr 0&0&\frac{19}{4}} \right] ]\;

線代啟示錄文章中的範例
(%i5)
A:matrix([3,-1,2],
              [6,-1,5],
              [-9,7,3]);
(%o5)　 \left[ \matrix{3&-1&2\cr 6&-1&5\cr -9&7&3} \right]

(%i6)　[L,U]: LU(A);
初始值L=\left[ \matrix{1&0&0\cr 0&1&0\cr 0&0&1} \right],U=\left[ \matrix{3&-1&2\cr 6&-1&5\cr -9&7&3} \right]
------------------
設 \displaystyle L [\; 2,1 ]\;=U[\; 2,1 ]\; / U[\; 1,1 ]\;=2 ,L=\left[ \matrix{1&0&0\cr 2&1&0\cr 0&0&1} \right]
原本U=\left[ \matrix{3&-1&2\cr 6&-1&5\cr -9&7&3} \right],第2列-2*第1列,得到U=\left[ \matrix{3&-1&2\cr 0&1&1\cr -9&7&3} \right]
------------------
設 \displaystyle L [\; 3,1 ]\;=U[\; 3,1 ]\; / U[\; 1,1 ]\;=-3 ,L=\left[ \matrix{1&0&0\cr 2&1&0\cr -3&0&1} \right]
原本U=\left[ \matrix{3&-1&2\cr 0&1&1\cr -9&7&3} \right],第3列--3*第1列,得到U=\left[ \matrix{3&-1&2\cr 0&1&1\cr 0&4&9} \right]
------------------
設 \displaystyle L [\; 3,2 ]\;=U[\; 3,2 ]\; / U[\; 2,2 ]\;=4 ,L=\left[ \matrix{1&0&0\cr 2&1&0\cr -3&4&1} \right]
原本U=\left[ \matrix{3&-1&2\cr 0&1&1\cr 0&4&9} \right],第3列-4*第2列,得到U=\left[ \matrix{3&-1&2\cr 0&1&1\cr 0&0&5} \right]
------------------
(%o6)　 [\; \left[ \matrix{1&0&0\cr 2&1&0\cr -3&4&1} \right],\left[ \matrix{3&-1&2\cr 0&1&1\cr 0&0&5} \right] ]\;

利用maxima內建指令得到相同的答案
(%i7)　[P,L,U]:get_lu_factors(lu_factor(A));
(%o7)　 [\;\left[ \matrix{1&0&0\cr 0&1&0\cr 0&0&1} \right],\left[ \matrix{1&0&0\cr 2&1&0\cr -3&4&1} \right],\left[ \matrix{3&-1&2\cr 0&1&1\cr 0&0&5} \right] ]\;

[ 本帖最後由 bugmens 於 2016-9-17 11:24 AM 編輯 ]

但不是每個矩陣都可以用LR法求得全部的特徵值，William Edward Stephenson提供一個特徵值無法收斂的例子。
http://krex.k-state.edu/dspace/b ... e=1&isAllowed=y

在第21頁底下，該矩陣為 A_1=\left[ \matrix{1&-1&1\cr 4&6&-1 \cr 4&4&1} \right] ，特徵值 \matrix{\lambda_1=5 \cr \lambda_2=2 \cr \lambda_3=1}
其中X=\left[ \matrix{0&-1&-1\cr 1&1&1 \cr 1&0&1} \right] ， Y=X^{-1}=\left[ \matrix{1&1&0 \cr 0&1&-1 \cr -1&-1&1} \right] ，D=\left[ \matrix{5&0&0 \cr 0&2&0 \cr 0&0&1} \right] ，符合 YAX=D
因為矩陣X的leading principal minor為0，導致原本的LR法無法收斂。

原本的LR法
A_2=\left[ \matrix{1&-0.2&1 \cr 20&6&-5 \cr 4&0.8&1} \right] ， A_3=\left[ \matrix{1&-0.04&1 \cr 100&6&-25 \cr 4&0.16&1} \right] ， A_4=\left[ \matrix{1&-0.008&1\cr 500&6&-125 \cr 4&0.032&1} \right]
發現對角線都是 \left[ 1,6,1 \right] ，無法收斂到正確的特徵值 \left[ 5,2,1 \right] ，以下是利用maxima呈現 A_2,A_3,A_4 的計算過程


(%i1)
A:matrix([1,-1,1],
              [4,6,-1],
              [4,4,1]);
(%o1)　 \left[ \matrix{1&-1&1\cr 4&6&-1 \cr 4&4&1}\right]

(%i2)
for i:2 thru 4 do
  ([P,L,U]:get_lu_factors(lu_factor(A)),
    print("計算矩陣A的LU分解"),
    print(A,"=",L,"*",U),
    A:U.L,
    print("將L和U對調後相乘得到下一個矩陣A"),
    print(U,"*",L,"=",A),
    eig:create_list(A[i,i],i,1,length(A)),
    print("特徵值",eig),
    print("------------------")
  );
0 errors, 0 warnings
計算矩陣A的LU分解
\left[ \matrix{1&-1&1\cr 4&6&-1 \cr 4&4&1} \right] =\left[ \matrix{\displaystyle 1&0&0\cr 4&1&0 \cr 4&\frac{4}{5}&1} \right]*\left[ \matrix{1&-1&1\cr 0&10&-5\cr 0&0&1}\right]
將L和U對調後相乘得到下一個矩陣A
\left[ \matrix{1&-1&1\cr 0&10&-5\cr 0&0&1} \right] *\left[ \matrix{\displaystyle 1&0&0\cr 4&1&0 \cr 4&\frac{4}{5}&1}\right]=\left[ \matrix{\displaystyle 1&-\frac{1}{5}&1\cr 20&6&-5 \cr 4&\frac{4}{5}&1}\right]
特徵值 \left[ 1,6,1 \right]
------------------
計算矩陣A的LU分解
\left[ \matrix{\displaystyle 1&-\frac{1}{5}&1\cr 20&6&-5\cr 4&\frac{4}{5}&1}\right] =\left[\matrix{\displaystyle 1&0&0\cr 20&1&0 \cr 4&\frac{4}{25}&1}\right]*\left[ \matrix{\displaystyle 1&-\frac{1}{5}&1\cr 0&10&-25 \cr 0&0&1} \right]
將L和U對調後相乘得到下一個矩陣A
\left[ \matrix{\displaystyle 1&-\frac{1}{5}&1\cr 0&10&-25 \cr 0&0&1}\right]*\left[ \matrix{\displaystyle 1&0&0\cr 20&1&0 \cr 4&\frac{4}{25}&1}\right]=\left[ \matrix{\displaystyle 1&-\frac{1}{25}&1\cr 100&6&-25 \cr 4&\frac{4}{25}&1}\right]
特徵值 \left[ 1,6,1 \right]
------------------
計算矩陣A的LU分解
\left[ \matrix{\displaystyle 1&-\frac{1}{25}&1 \cr 100&6&-25 \cr 4&\frac{4}{25}&1}\right]= \left[ \matrix{\displaystyle 1&0&0\cr 100&1&0 \cr 4&\frac{4}{125}&1}\right] * \left[ \matrix{\displaystyle 1&-\frac{1}{25}&1\cr 0&10&-125 \cr 0&0&1}\right]  
將L和U對調後相乘得到下一個矩陣A
\left[ \matrix{\displaystyle 1&-\frac{1}{25}&1\cr 0&10&-125 \cr 0&0&1}\right]*\left[ \matrix{\displaystyle 1&0&0\cr 100&1&0 \cr 4&\frac{4}{125}&1}\right]=\left[ \matrix{\displaystyle 1&-\frac{1}{125}&1\cr 500&6&-125 \cr 4&\frac{4}{125}&1}\right]
特徵值 \left[ 1,6,1 \right]
------------------
(%o2)　done



改進的方法稱為LR with interchanges，也就是將其中兩列交換後再利用LR法計算
\left[ \matrix{1&-1&1 \cr 4&6&-1 \cr 4&4&1}\right] 第1列和第2列交換 \left[ \matrix{4&6&-1\cr 1&-1&1 \cr 4&4&1}\right] 第2列 \displaystyle -\frac{1}{4} 第1列 \left[ \matrix{4&6&-1 \cr 0&-2.5&1.25 \cr 4&4&1}\right]
\left[ \matrix{4&6&-1 \cr 0&-2.5&1.25 \cr 4&4&1}\right] 第3列 - 第1列 \left[ \matrix{4&6&-1 \cr 0&-2.5&1.25 \cr 0&-2&2}\right] 第3列-0.8第2列 \left[ \matrix{4&6&-1\cr 0&-2.5&1.25 \cr 0&0&1}\right]=R ，就是LU分解的U矩陣

\left[ \matrix{4&6&-1\cr 0&-2.5&1.25 \cr 0&0&1}\right] 第1行和第2行交換 \left[ \matrix{6&4&-1\cr -2.5&0&1.25 \cr 0&0&1}\right] 第1行 \displaystyle +\frac{1}{4} 第2行 \left[ \matrix{7&4&-1\cr -2.5&0&1.25 \cr 0&0&1}\right]
\left[ \matrix{7&4&-1\cr -2.5&0&1.25 \cr 0&0&1}\right] 第1行+第3行 \left[ \matrix{6&4&-1\cr -1.25&0&1.25 \cr 1&0&1}\right] 第2行+0.8第3行 \left[ \matrix{6&3.2&-1\cr -1.25&1&1.25 \cr 1&0.8&1}\right]=A_2
第一部份先將第1列和第2列交換後針對列做高斯消去法，將運算過程記錄下來於第二部分對行做相對應運算得到A_2。
接下來就不用再列交換和行交換，就用原本的LR法計算出
A_3=\left[ \matrix{5.167&3.040&-1.00\cr -0.174&1.833&1.042\cr 0.167&0.160&1.000} \right] ，A_4=\left[ \matrix{5.032&3.008&-1.00\cr -0.03&1.968&1.008\cr 0.032&0.032&1.000} \right]，全部計算過程利用maxima呈現。



設定為數值運算
(%i1)　numer:true;
(%o1)　true

不顯示小數轉換成分數的過程
例如rat: replaced -4.25 by -17/4 = -4.25
(%i2)　ratprint:false;
(%o2)　false

僅顯示小數5位
(%i3)　fpprintprec:5;
(%o3)　5

列交換的LR迭代法副程式
(%i4)
LRwithInterchanges(A,row1,row2,epsilon):=block
([count:0,P,L,U,n:length(A),eig1,eig2],
eig2:makelist(0,n),
do
   (count:count+1,
    print("第",count,"輪"),
    if count=1 then
      (print(A,"第",row1,"列和第",row2,"列交換",A:rowswap(A,row1,row2))),
    print("計算矩陣A的LU分解"),
    [P,L,U]:get_lu_factors(lu_factor(A)),
    print(A,"=",L,"*",U),
    if count=1 then
       (print("L=",L,"第",row1,"列和第",row2,"列交換",L:rowswap(L,row1,row2))),
     print("將L和U對調後相乘得到下一個矩陣A"),
     A:U.L,
     print(U,"*",L,"=",A),
     eig1:create_list(A[i,i],i,1,n),
     print("特徵值",eig1),
     print("------------------"),
     if lmax(abs(eig1-eig2))<epsilon then
        (return(eig1)),
     eig2:eig1
    )
)$

原本LR法特徵值無法收斂的矩陣A
(%i5)
A:matrix([1,-1,1],
              [4,6,-1],
              [4,4,1]);
(%o5)　 \left[ \matrix{1&-1&1 \cr 4&6&-1 \cr 4&4&1} \right]

指定第1列和第2列交換,計算矩陣A全部的特徵值,誤差10^-5
(%i6)　LRwithInterchanges(A,1,2,10^-5);
第1輪
0 errors, 0 warnings
\left[ \matrix{1&-1&1 \cr 4&6&-1 \cr 4&4&1}\right] 第1列和第2列交換 \left[ \matrix{4&6&-1\cr 1&-1&1 \cr 4&4&1}\right]
計算矩陣A的LU分解
\left[ \matrix{4&6&-1\cr 1&-1&1 \cr 4&4&1}\right]=\left[ \matrix{1&0&0\cr 0.25&1&0 \cr 1.0&0.8&1}\right] * \left[ \matrix{4&6&-1\cr 0&-2.5&1.25 \cr 0&0&1.0}\right]
L=\left[ \matrix{1&0&0\cr 0.25&1&0 \cr 1.0&0.8&1}\right] 第1列和第2列交換 \left[ \matrix{0.25&1&0\cr 1&0&0 \cr 1.0&0.8&1}\right]
將L和U對調後相乘得到下一個矩陣A
\left[ \matrix{4&6&-1\cr 0&-2.5&1.25 \cr 0&0&1.0}\right] * \left[ \matrix{0.25&1&0\cr 1&0&0 \cr 1.0&0.8&1}\right]=\left[ \matrix{6.0&3.2&-1\cr -1.25&1.0&1.25 \cr 1.0&0.8&1.0}\right]
特徵值 \left[6.0,1.0,1.0 \right]
------------------
第2輪
計算矩陣A的LU分解
\left[ \matrix{6.0&3.2&-1 \cr -1.25&1.0&1.25 \cr 1.0&0.8&1.0}\right]=\left[ \matrix{1&0&0\cr -0.20833&1&0 \cr 0.16667&0.16&1}\right] * \left[ \matrix{6.0&3.2&-1\cr 0&1.6667&1.0417 \cr 0&0&1.0}\right]
將L和U對調後相乘得到下一個矩陣A
\left[ \matrix{6.0&3.2&-1\cr 0&1.6667&1.0417 \cr 0&0&1.0}\right]*\left[ \matrix{1&0&0\cr -0.20833&1&0 \cr 0.16667&0.16&1}\right]=\left[ \matrix{5.1667&3.04&-1.0\cr -0.17361&1.8333&1.0417 \cr 0.16667&0.16&1.0}\right]
特徵值 \left[ 5.1667,1.8333,1.0 \right]
------------------
第3輪
計算矩陣A的LU分解
\left[ \matrix{5.1667&3.04&-1.0\cr -0.17361&1.8333&1.0417 \cr 0.16667&0.16&1.0}\right]=\left[ \matrix{1&0&0\cr -0.033602&1&0 \cr 0.032258&0.032&1}\right]*\left[ \matrix{5.1667&3.04&-1.0\cr 0&1.9355&1.0081 \cr 0&0&1.0}\right]
將L和U對調後相乘得到下一個矩陣A
\left[ \matrix{5.1667&3.04&-1.0\cr 0&1.9355&1.0081 \cr 0&0&1.0}\right] *\left[ \matrix{1&0&0\cr -0.033602&1&0 \cr 0.032258&0.032&1}\right]=\left[ \matrix{5.0323&3.008&-1.0\cr -0.032518&1.9677&1.0081 \cr 0.032258&0.032&1.0}\right]
特徵值 \left[5.0323,1.9677,1.0 \right]
------------------
第4輪
計算矩陣A的LU分解
\left[ \matrix{5.0323&3.008&-1.0\cr -0.032518&1.9677&1.0081 \cr 0.032258&0.032&1.0}\right]=\left[ \matrix{1&0&0\cr -0.006462&1&0 \cr 0.0064103&0.0064&1}\right] * \left[ \matrix{5.0323&3.008&-1.0 \cr 0&1.9872&1.0016 \cr 0&0&1.0}\right]
將L和U對調後相乘得到下一個矩陣A
\left[ \matrix{5.0323&3.008&-1.0\cr 0&1.9872&1.0016 \cr 0&0&1.0}\right]*\left[ \matrix{1&0&0\cr -0.006462&1&0 \cr 0.0064103&0.0064&1}\right]=\left[ \matrix{5.0064&3.0016&-1.0\cr -0.0064205&1.9936&1.0016 \cr 0.0064103&0.0064&1.0}\right]
特徵值 \left[5.0064,1.9936,1.0 \right]
------------------
第5輪
計算矩陣A的LU分解
\left[ \matrix{5.0064&3.0016&-1.0\cr -0.0064205&1.9936&1.0016 \cr 0.0064103&0.0064&1.0}\right]=\left[ \matrix{1&0&0\cr -0.0012825&1&0 \cr 0.0012804&0.00128&1}\right]*\left[ \matrix{5.0064&3.0016&-1.0\cr 0&1.9974&1.0003 \cr 0&0&1.0}\right]
將L和U對調後相乘得到下一個矩陣A
\left[ \matrix{5.0064&3.0016&-1.0\cr 0&1.9974&1.0003 \cr 0&0&1.0}\right]*\left[ \matrix{1&0&0\cr -0.0012825&1&0 \cr 0.0012804&0.00128&1}\right]=\left[ \matrix{5.0013&3.0003&-1.0\cr -0.0012808&1.9987&1.0003 \cr 0.0012804&0.00128&1.0}\right]
特徵值 \left[5.0013,1.9987,1.0 \right]
------------------
第6輪
計算矩陣A的LU分解
\left[ \matrix{5.0013&3.0003&-1.0\cr -0.0012808&1.9987&1.0003 \cr 0.0012804&0.00128&1.0}\right]=\left[ \matrix{1&0&0\cr -2.561 \cdot 10^{-4}&1&0 \cr 2.5602 \cdot 10^{-4}&2.56 \cdot 10^{-4}&1}\right]*\left[ \matrix{5.0013&3.0003&-1.0 \cr 0&1.9995&1.0001 \cr 0&0&1.0}\right]
將L和U對調後相乘得到下一個矩陣A
\left[ \matrix{5.0013&3.0003&-1.0\cr 0&1.9995&1.0001 \cr 0&0&1.0}\right]*\left[ \matrix{1&0&0\cr -2.561 \cdot 10^{-4}&1&0 \cr 2.5602 \cdot 10^{-4}&2.56 \cdot 10^{-4}&1}\right]=\left[ \matrix{5.0003&3.0001&-1.0 \cr -2.5603 \cdot 10^{-4}&1.9997&1.0001 \cr 2.5602 \cdot 10^{-4}&2.56 \cdot 10^{-4}&1.0}\right]
特徵值 \left[5.0003,1.9997,1.0 \right]
------------------
第7輪
計算矩陣A的LU分解
\left[ \matrix{5.0003&3.0001&-1.0\cr -2.5603 \cdot 10^{-4}&1.9997&1.0001 \cr 2.5602 \cdot 10^{-4}&2.56 \cdot 10^{-4}&1.0}\right]=\left[ \matrix{1&0&0\cr -5.1204 \cdot 10^{-5}&1&0 \cr 5.1201 \cdot 10^{-5}&5.12 \cdot 10^{-5}&1}\right] * \left[ \matrix{5.0003&3.0001&-1.0\cr 0&1.9999&1.0 \cr 0&0&1.0}\right]
將L和U對調後相乘得到下一個矩陣A
\left[ \matrix{5.0003&3.0001&-1.0\cr 0&1.9999&1.0 \cr 0&0&1.0}\right]*\left[ \matrix{1&0&0\cr -5.1204 \cdot 10^{-5}&1&0 \cr 5.1201 \cdot 10^{-5}&5.12 \cdot 10^{-5}&1}\right]=\left[ \matrix{5.0001&3.0&-1.0\cr -5.1201 \cdot 10^{-5}&1.9999&1.0 \cr 5.1201 \cdot 10^{-5}&5.12 \cdot 10^{-5}&1.0}\right]
特徵值 \left[5.0001,1.9999,1.0\right]
------------------
第8輪
計算矩陣A的LU分解
\left[ \matrix{5.0001&3.0&-1.0\cr -5.1201 \cdot 10^{-5}&1.9999&1.0 \cr 5.1201 \cdot 10^{-5}&5.12 \cdot 10^{-5}&1.0}\right]=\left[ \matrix{1&0&0\cr -1.024 \cdot 10^{-5}&1&0 \cr 1.024 \cdot 10^{-5}&1.024 \cdot 10^{-5}&1}\right]*\left[ \matrix{5.0001&3.0&-1.0\cr 0&2.0&1.0 \cr 0&0&1.0}\right]
將L和U對調後相乘得到下一個矩陣A
\left[ \matrix{5.0001&3.0&-1.0\cr 0&2.0&1.0 \cr 0&0&1.0}\right]*\left[ \matrix{1&0&0\cr -1.024 \cdot 10^{-5}&1&0 \cr 1.024 \cdot 10^{-5}&1.024 \cdot 10^{-5}&1}\right]=\left[ \matrix{5.0&3.0&-1.0\cr -1.024 \cdot 10^{-5}&2.0&1.0 \cr 1.024 \cdot 10^{-5}&1.024 \cdot 10^{-5}&1.0}\right]
特徵值 \left[5.0,2.0,1.0 \right]
------------------
第9輪
計算矩陣A的LU分解
\left[ \matrix{5.0&3.0&-1.0\cr -1.024 \cdot 10^{-5}&2.0&1.0 \cr 1.024 \cdot 10^{-5}&1.024 \cdot 10^{-5}&1.0}\right]=\left[ \matrix{1&0&0 \cr -2.048 \cdot 10^{-6}&1&0 \cr 2.048 \cdot 10^{-6}&2.048 \cdot 10^{-6}&1}\right]*\left[ \matrix{5.0&3.0&-1.0 \cr 0&2.0&1.0 \cr 0&0&1.0}\right]
將L和U對調後相乘得到下一個矩陣A
\left[ \matrix{5.0&3.0&-1.0\cr 0&2.0&1.0 \cr 0&0&1.0}\right]*\left[ \matrix{1&0&0\cr -2.048 \cdot 10^{-6}&1&0 \cr 2.048 \cdot 10^{-6}&2.048 \cdot 10^{-6}&1}\right]=\left[ \matrix{5.0&3.0&-1.0\cr -2.048 \cdot 10^{-6}&2.0&1.0 \cr 2.048 \cdot 10^{-6}&2.048 \cdot 10^{-6}&1.0}\right]
特徵值 \left[5.0,2.0,1.0 \right]
------------------
(%o6)　 \left[5.0,2.0,1.0 \right]

[ 本帖最後由 bugmens 於 2016-9-26 05:04 PM 編輯 ]

1961年，英國計算機科學家弗朗西斯 (John G.F. Francis) 和俄國數學家卡布諾夫斯凱雅 (Vera Kublanovskaya) 獨立提出了 QR演算法。經過五十年的演化，QR演算法逐漸成為當今特徵值和特徵向量數值計算的骨幹，並被視為20世紀最重要的十個演算法之一。另一方面，由於 LR演算法的穩定性不如QR演算法，今天已鮮少被提及。

節錄自線代啟示錄
https://ccjou.wordpress.com/2013/12/10/qr-演算法-上/

弗朗西斯1961年的論文
http://comjnl.oxfordjournals.org/content/4/3/265.full.pdf+html
－－－－－－－－－－－－－－－－－－－－－－－－－－－
虛擬碼
1.對於k=0,1,2,\ldots，直到 \lambda(1),\lambda(2),\ldots,\lambda(n) 收斂為止，計算
2.A_k=Q_k R_k，計算A_k的QR分解。
3.A_{k+1}=R_k Q_k，將Q_k和R_k對調後相乘得到下一個矩陣A_{k+1}
4.特徵值 \lambda(1),\lambda(2),\ldots,\lambda(n) 在矩陣A_{k+1}的對角線上。
－－－－－－－－－－－－－－－－－－－－－－－－－－－
計算 A=\left[ \matrix{-4&0&1\cr -1&3&1\cr 1&2&5} \right] 全部的特徵值。
[計算過程]
第1輪：
計算矩陣A的QR分解
\left[ \matrix{-4&0&1\cr -1&3&1\cr 1&2&5} \right]=\left[ \matrix{-0.9428&0.0618&0.3276\cr -0.2357&-0.8184&-0.5241\cr 0.2357&-0.5713&0.7861} \right]*\left[ \matrix{4.2426&-0.2357&2.2204 \cdot 10^{-16}\cr 0.0&-3.5978&-3.6133\cr 0.0&0.0&3.7342} \right]
將Q和R對調後相乘得到下一個矩陣A
\left[ \matrix{4.2426&-0.2357&2.2204 \cdot 10^{-16}\cr 0.0&-3.5978&-3.6133\cr 0.0&0.0&3.7342} \right]*\left[ \matrix{- 0.9428&0.0618&0.3276\cr -0.2357&-0.8184&-0.5241\cr 0.2357&-0.5713&0.7861} \right]=\left[ \matrix{-3.9444&0.4549&1.5133\cr - 0.0036&5.0088&-0.9549\cr 0.8802&-2.1335&2.9356} \right]
特徵值 \left[ -3.9444,5.0088,2.9356 \right]

第2輪：
計算矩陣A的QR分解
\left[ \matrix{-3.9444&0.4549&1.5133\cr -0.0036&5.0088&-0.9549\cr 0.8802&-2.1335&2.9356} \right]=\left[ \matrix{- 0.9760&0.0810&0.2022\cr -9.0056 \cdot 10^{-4}&-0.9298&0.3682\cr 0.2178&0.3592&0.9075} \right]*\left[ \matrix{4.0415&-0.9132&- 0.8367\cr 0.0&-5.3864&2.0648\cr 0.0&0.0&2.6184} \right]
將Q和R對調後相乘得到下一個矩陣A
\left[ \matrix{4.0415&-0.9132&-0.8367\cr 0.0&-5.3864&2.0648\cr 0.0&0.0&2.6184} \right]* \left[ \matrix{- 0.9760&0.0810&0.2022\cr -9.0056 \cdot 10^{-4}&-0.9298&0.3682\cr 0.2178&0.3592&0.9075} \right]= \left[ \matrix{-4.1259&0.8759&- 0.2785\cr 0.4545&5.7496&-0.1093\cr 0.5702&0.9404&2.3762} \right]
特徵值 \left[ -4.1259,5.7496,2.3762 \right]

第3輪：
計算矩陣A的QR分解
\left[ \matrix{-4.1259&0.8759&-0.2785\cr 0.4545&5.7496&-0.1093\cr 0.5702&0.9404&2.3762} \right]= \left[ \matrix{-0.9847&- 0.1302&0.1155\cr 0.1085&-0.9781&-0.1775\cr 0.1361&-0.1622&0.9773} \right]* \left[ \matrix{4.1898&-0.1107&0.5858\cr 0.0&- 5.8904&-0.2423\cr 0.0&0.0&2.3096} \right]
將Q和R對調後相乘得到下一個矩陣A
\left[ \matrix{4.1898&-0.1107&0.5858\cr 0.0&-5.8904&-0.2423\cr 0.0&0.0&2.3096} \right]*\left[ \matrix{-0.9847&- 0.1302&0.1155\cr 0.1085&-0.9781&-0.1775\cr 0.1361&-0.1622&0.9773} \right]=\left[ \matrix{-4.0581&-0.5321&1.0763\cr - 0.6720&5.8009&0.8085\cr 0.3143&-0.3746&2.2572} \right]
特徵值 \left[ -4.0581,5.8009,2.2572 \right]

反覆計算直到全部特徵值 \left[ 5.7542,-4.1444,2.3902 \right]


－－－－－－－－－－－－－－－－－－－－－－－－－－－－
這次程式碼需要安裝https://sourceforge.net/projects ... ows/5.38.1-Windows/
的maxima-clisp-sbcl-5.38.1.exe，若安裝maxima-sbcl-5.38.1-win64.exe在執行load(lapack);時出現以下的錯誤訊息。
While evaluating the form starting at line 13, column 0
  of #P"C:/Program Files (x86)/Maxima-sbcl-5.38.1/share/maxima/5.38.1/share/lapack/load-lapack.lisp":;
                                                                                                     ; compilation unit aborted
                                                                                                     ;   caught 1 fatal ERROR condition
;
; compilation unit aborted
;   caught 1 fatal ERROR condition
;
; compilation unit aborted
;   caught 1 fatal ERROR condition
loadfile: failed to load C:/Program Files (x86)/Maxima-sbcl-5.38.1/share/maxima/5.38.1/share/lapack/load-lapack.lisp
-- an error. To debug this try: debugmode(true);
－－－－－－－－－－－－－－－－－－－－－－－－－－－－


要先載入lapack才能使用dgeqrf指令
(%i1)　load(lapack);
;;  Loading file C:\Users\user\maxima\binary\5_38_1_5_gdf93b7b_dirty\clisp\2_49__2010_07_07___built_on_toshiba__192_168_43_3__\share\lapack\lapack-package.fas ...
;;  Loaded file C:\Users\user\maxima\binary\5_38_1_5_gdf93b7b_dirty\clisp\2_49__2010_07_07___built_on_toshiba__192_168_43_3__\share\lapack\lapack-package.fas
...有幾十行類似的訊息...
;;  Loading file C:\Users\user\maxima\binary\5_38_1_5_gdf93b7b_dirty\clisp\2_49__2010_07_07___built_on_toshiba__192_168_43_3__\share\lapack\dgemm.fas ...
;;  Loaded file C:\Users\user\maxima\binary\5_38_1_5_gdf93b7b_dirty\clisp\2_49__2010_07_07___built_on_toshiba__192_168_43_3__\share\lapack\dgemm.fas
0 errors, 0 warnings
(%o1)　C:\maxima-5.38.1\share\maxima\5.38.1_5_gdf93b7b_dirty\share\lapack\lapack.mac

QR分解副程式
(%i2)
QR(A,epsilon):=block
([count:0,Q,R,n:length(A),eig1,eig2],
eig2:makelist(0,n),
do
   (count:count+1,
    /*print("第",count,"輪"),*/
    [Q,R]:dgeqrf(A),
    /*print("計算矩陣A的QR分解",A,"="),*/
    /*print(Q,"*",R),*/
    A:R.Q,
    /*print("將Q和R對調後相乘得到下一個矩陣A"),*/
    /*print(R,"*",Q,"=",A),*/
    eig1:create_list(A[i,i],i,1,n),
    print("特徵值",eig1),
    if lmax(abs(eig1-eig2))<epsilon then
      (return(eig1)),
    eig2:eig1
  )
)$

矩陣A
(%i3)
A:matrix([-4,0,1],
              [-1,3,1],
              [1,2,5]);
(%o3)　 \left[ \matrix{-4&0&1\cr -1&3&1 \cr 1&2&5} \right]

計算矩陣A全部的特徵值,誤差10^-5
(%i4)　QR(A,10^-5);
特徵值 \left[-3.944444444444443,5.008822126847878,2.935622317596567\right]
特徵值 \left[-4.125850340136053,5.749602762836074,2.37624757729998\right]
特徵值 \left[-4.05812826971517,5.800914022263142,2.257214247452029\right]
特徵值 \left[-4.062424573627523,5.699337059058267,2.363087514569257\right]
特徵值 \left[-3.9888058478382,5.63243637676476,2.356369471073442\right]
特徵值 \left[-3.846431951098917,5.458191384059295,2.388240567039621\right]
特徵值 \left[-3.626635311452812,5.244350090320102,2.382285221132707\right]
特徵值 \left[-3.147542963553108,4.756572866807507,2.390970096745598\right]
特徵值 \left[-2.447555281262403,4.059406511538823,2.388148769723577\right]
特徵值 \left[-1.273077705628198,2.882396372616165,2.390681333012031\right]
特徵值 \left[0.1887543282710361,1.421650419107125,2.389595252621835\right]
特徵值 \left[1.792441558902001,-0.1828180203586464,2.39037646145664\right]
特徵值 \left[3.236224515952411,-1.626214507426892,2.389989991474478\right]
特徵值 \left[4.228710740781744,-2.618950141448976,2.390239400667228\right]
特徵值 \left[4.9356478091601,-3.325754357659397,2.390106548499294\right]
特徵值 \left[5.2817166444885,-3.671904365492856,2.390187721004352\right]
特徵值 \left[5.52698054546436,-3.917123340678473,2.390142795214107\right]
特徵值 \left[5.617321691451712,-4.007491179509199,2.39016948805748\right]
特徵值 \left[5.695768989139763,-4.085923407784471,2.3901544186447\right]
特徵值 \left[5.71433037843271,-4.104493622778735,2.390163244346017\right]
特徵值 \left[5.740366480751086,-4.130524691091749,2.390158210340657\right]
特徵值 \left[5.742069167509259,-4.13223030430456,2.390161136795295\right]
特徵值 \left[5.75149968365114,-4.141659142342601,2.390159458691456\right]
特徵值 \left[5.750232183663825,-4.140392614159591,2.39016043049576\right]
特徵值 \left[5.754015952686897,-4.144175824385144,2.390159871698241\right]
特徵值 \left[5.752779110658064,-4.142939305317683,2.390160194659612\right]
特徵值 \left[5.754442586633282,-4.144602595320878,2.39016000868759\right]
特徵值 \left[5.753645152532918,-4.143805268593751,2.390160116060828\right]
特徵值 \left[5.754427402570066,-4.144587456755889,2.390160054185819\right]
特徵值 \left[5.753971827867405,-4.144131917758443,2.390160089891033\right]
特徵值 \left[5.754355892472823,-4.14451596178042,2.390160069307593\right]
特徵值 \left[5.754108293244391,-4.14426837442645,2.390160081182056\right]
特徵值 \left[5.754301676088646,-4.144461750423915,2.390160074335267\right]
特徵值 \left[5.754170203599655,-4.144330281884226,2.390160078284571\right]
特徵值 \left[5.754268945170628,-4.144429021177805,2.390160076007176\right]
特徵值 \left[5.754199929325,-4.144360006645707,2.390160077320705\right]
特徵值 \left[5.754250726832598,-4.144410803395806,2.390160076563208\right]
特徵值 \left[5.754214706157051,-4.144374783157142,2.390160077000091\right]
特徵值 \left[5.754240942756682,-4.14440101950482,2.390160076748139\right]
特徵值 \left[5.754222198416984,-4.144382275310429,2.390160076893448\right]
特徵值 \left[5.754235777643633,-4.144395854453276,2.390160076809646\right]
特徵值 \left[5.754226038367267,-4.144386115225238,2.390160076857978\right]
(%o4)　 \left[5.754226038367267,-4.144386115225238,2.390160076857978 \right]