Board logo

標題: 用maxima學數值分析-特徵值和特徵向量 [打印本頁]

作者: bugmens    時間: 2016-7-27 14:27     標題: 用maxima學數值分析-特徵值和特徵向量

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

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

\(A=\left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right]\),\( \matrix{圓心(-4,0),半徑|\;0|\;+|\;1|\;=1 \cr 圓心(3,0),半徑|\;-1|\;+|\;1|\;=2 \cr 圓心(5,0),半徑|\;1|\;+|\;2|\;=3} \)
矩陣\(A\)和轉置矩陣\(A^T\)有相同的特徵值
\(A^T=\left[ \matrix{-4&-1&1 \cr 0&3&2 \cr 1&1&5} \right]\),\( \matrix{圓心(-4,0),半徑|\;-1|\;+|\;1|\;=2 \cr 圓心(3,0),半徑|\;0|\;+|\;2|\;=2 \cr 圓心(5,0),半徑|\;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) \(\left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right]\)

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);
圓心\((-4,0)\),半徑1
圓心\((3,0)\),半徑2
圓心\((5,0)\),半徑3
(%o4) done


計算矩陣A^t的Gershgorin Circle
(%i5) GershgorinCircle(transpose(A));
圓心\((-4,0)\),半徑2
圓心\((3,0)\),半徑2
圓心\((5,0)\),半徑2
(%o5) done


[ 本帖最後由 bugmens 於 2016-7-28 09:37 AM 編輯 ]
作者: bugmens    時間: 2016-7-31 17:39

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

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

矩陣\(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{Au_k}{\Vert\; Au_k \Vert\;}\)
4.\(\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]\)

第1輪:
\(Au_0=\left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5}\right] \left[ \matrix{1 \cr 1 \cr 1} \right]=\left[ \matrix{-3 \cr 3 \cr 8} \right] \)
特徵向量\(\displaystyle u_1=\frac{Au_0}{\Vert\;Au_0\Vert\;}=\frac{\left[ \matrix{-3\cr 3 \cr 8}\right]}{\sqrt{(-3)^2+3^2+8^2}}=\left[ \matrix{-0.3313 \cr 0.3313 \cr 0.8835} \right] \)
特徵值\( \lambda_1=u_1^T A u_1=\left[ \matrix{-0.3313 & 0.3313 & 0.8835} \right] \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]=4.1951 \)

第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 編輯 ]
作者: bugmens    時間: 2016-8-4 07:08

\(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 編輯 ]
作者: bugmens    時間: 2016-8-6 08:46

\(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 編輯 ]
作者: bugmens    時間: 2016-8-11 17:34

\(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次

作者: bugmens    時間: 2016-9-15 10:54

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 編輯 ]
作者: bugmens    時間: 2016-9-17 11:23

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 編輯 ]
作者: bugmens    時間: 2016-9-25 17:09

但不是每個矩陣都可以用\(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 編輯 ]
作者: bugmens    時間: 2016-9-27 09:24

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] \)
作者: bugmens    時間: 2016-10-1 00:44

1.求QR分解的方法--Gram-Schmidt

詳細方法請見線代啟示錄的"Gram-Schmidt 正交化與 QR 分解"
https://ccjou.wordpress.com/2010 ... %E5%88%86%E8%A7%A3/



利用Gram-Schmidt正交化的QR分解
(%i1)
QRbyGramSchmidt(A):=block
([columns:length(A[1]),rank:rank(A),x,y,normy,q,Q,R],
x:create_list(col(A,i),i,1,columns),
y:makelist(0,columns),
normy:makelist(0,columns),
print("1.正交化"),
for i:1 thru rank do
   (print("y",i,"=x",i,"=",y[ i ]:x[ i ]),
    for j:1 thru i-1 do
      (print("y",i,"=y",i,"-(y",j,"^T.x",i,")/(|y",j,"|^2)y",j),
       print("=",y[ i ],"-(",transpose(y[j]),".",col(A,i),")/(",y[j],".",y[j],")",y[j]),
       numerator:transpose(y[j]).col(A,i),
       denominator:y[j].y[j],
       print("=",y[ i ],"-",numerator,"/",denominator,y[j]),
       print("=",y[ i ]:y[ i ]-numerator/denominator*y[j])
      ),
    normy[ i ]:sqrt(y[ i ].y[ i ]),
    print("|y",i,"|=sqrt(",y[ i ],".",y[ i ],")=",normy[ i ]),
    print("------------------------")
   ),
   print("2.單位化得正交矩陣Q"),
   q:makelist(0,rank),
   for i:1 thru rank do
     (q[ i ]:ratsimp(y[ i ]/normy[ i ]),
      print("q",i,"=y",i,"/|y",i,"|=",y[ i ],"/",normy[ i ],"=",q[ i ])
     ),
    print("Q=",Q:transpose(lreduce(append,map(transpose,q)))),
    print("------------------------"),
    print("3.上三角矩陣R"),
    R:zeromatrix(rank,columns),
    for i:1 thru rank do
      (R[i,i]:normy[ i ],
       print("R[",i,",",i,"]=|y",i,"|=",R[i,i]),
       for j:i+1 thru columns do
         (R[i,j]:ratsimp(transpose(q[ i ]).x[j]),
          print("R[",i,",",j,"]=q",i,"^Tx",j=transpose(q[ i ]),".",x[j],"=",R[i,j])
         )
       ),
    print("R=",R),
    return([Q,R])
)$


本文使用的範例
(%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] \)

執行QR分解
(%i3) [Q,R]: QRbyGramSchmidt(A);
1.正交化
\(y1=x1=\left[ \matrix{-4 \cr -1 \cr 1} \right] \)
\( |\; y1 |\;=sqrt(\left[ \matrix{-4 \cr -1 \cr 1} \right]\cdot \left[ \matrix{-4 \cr -1 \cr 1} \right])=3 \sqrt{2} \)
------------------------
\( y2=x2=\left[ \matrix{0 \cr 3 \cr 2} \right] \)
\( y2=y2-(y1^T \cdot x2)/(|\; y1 |\;^2)y1 \)
\( =\left[ \matrix{0 \cr 3 \cr 2} \right]-(\left[ \matrix{-4&-1&1} \right]\cdot \left[ \matrix{0 \cr 3 \cr 2} \right])/(\left[ \matrix{-4 \cr -1 \cr 1} \right]\cdot \left[ \matrix{-4 \cr -1 \cr 1} \right])\left[ \matrix{-4 \cr -1 \cr 1} \right] \)
\( =\left[ \matrix{0 \cr 3 \cr 2}\right]- -1/18 \left[ \matrix{-4 \cr -1 \cr 1} \right] \)
\( =\left[ \matrix{\displaystyle -\frac{2}{9}\cr \frac{53}{18}\cr \frac{37}{18}} \right] \)
\( \displaystyle |\; y2 |\;=sqrt(\left[ \matrix{\displaystyle -\frac{2}{9}\cr \frac{53}{18}\cr \frac{37}{18}} \right]\cdot \left[ \matrix{\displaystyle -\frac{2}{9}\cr \frac{53}{18}\cr \frac{37}{18}} \right])=\frac{\sqrt{233}}{3 \sqrt{2}} \)
------------------------
\( y3=x3=\left[ \matrix{1 \cr 1 \cr 5} \right] \)
\( y3=y3-(y1^T \cdot x3)/(|\; y1 |\;^2)y1 \)
\( =\left[ \matrix{1 \cr 1 \cr 5} \right]-(\left[ \matrix{-4&-1&1} \right]\cdot \left[ \matrix{1 \cr 1 \cr 5} \right])/(\left[ \matrix{-4 \cr -1 \cr 1} \right]\cdot \left[ \matrix{-4 \cr -1 \cr 1} \right])\left[ \matrix{-4 \cr -1 \cr 1} \right] \)
\( =\left[ \matrix{1 \cr 1 \cr 5} \right]-0/18\left[ \matrix{-4 \cr -1 \cr 1} \right] \)
\( =\left[ \matrix{1 \cr 1 \cr 5} \right] \)
\( y3=y3-(y2^T \cdot x3)/(|\; y2 |\;^2)y2 \)
\( =\left[ \matrix{1 \cr 1 \cr 5} \right]-(\left[ \matrix{\displaystyle -\frac{2}{9}&\frac{53}{18}&\frac{37}{18}} \right]\cdot \left[ \matrix{1 \cr 1 \cr 5} \right])/(\left[ \matrix{\displaystyle -\frac{2}{9}\cr \frac{53}{18} \cr \frac{37}{18}} \right]\cdot \left[ \matrix{\displaystyle -\frac{2}{9}\cr \frac{53}{18} \cr \frac{37}{18}} \right])\left[ \matrix{\displaystyle -\frac{2}{9}\cr \frac{53}{18} \cr \frac{37}{18}} \right] \)
\( \displaystyle =\left[ \matrix{1 \cr 1 \cr 5} \right]-13/\frac{233}{18}\left[ \matrix{\displaystyle -\frac{2}{9}\cr \frac{53}{18} \cr \frac{37}{18}} \right] \)
\( =\left[ \matrix{\displaystyle \frac{285}{233}\cr -\frac{456}{233} \cr \frac{684}{233}} \right] \)
\( \displaystyle |\; y3 |\;=sqrt(\left[ \matrix{\displaystyle \frac{285}{233}\cr -\frac{456}{233} \cr \frac{684}{233}} \right]\cdot \left[ \matrix{\displaystyle \frac{285}{233}\cr -\frac{456}{233} \cr \frac{684}{233}} \right])=\frac{57}{\sqrt{233}} \)
------------------------
2.單位化得正交矩陣\(Q\)
\( q1=y1/|\; y1 |\;=\left[ \matrix{-4 \cr -1 \cr 1} \right]/3\sqrt{2}=\left[ \matrix{\displaystyle -\frac{2^{3/2}}{3}\cr -\frac{1}{3 \sqrt{2}}\cr \frac{1}{3\sqrt{2}}} \right] \)
\( \displaystyle q2=y2/|\; y2 |\;=\left[ \matrix{\displaystyle -\frac{2}{9}\cr \frac{53}{18} \cr \frac{37}{18}} \right] / \frac{\sqrt{233}}{3\sqrt{2}}=\left[ \matrix{\displaystyle -\frac{2^{3/2}}{3 \sqrt{233}} \cr \frac{53}{3 \sqrt{2} \sqrt{233}} \cr \frac{37}{3 \sqrt{2} \sqrt{233}}} \right] \)
\( \displaystyle q3=y3/|\; y3 |\;=\left[ \matrix{\displaystyle \frac{285}{233}\cr -\frac{456}{233} \cr \frac{684}{233}}\right]/ \frac{57}{\sqrt{233}}=\left[ \matrix{\displaystyle \frac{5}{\sqrt{233}}\cr -\frac{8}{\sqrt{233}}\cr \frac{12}{\sqrt{233}}} \right] \)
\( Q=\left[ \matrix{\displaystyle -\frac{2^{3/2}}{3}&-\frac{2^{3/2}}{3\sqrt{233}}&\frac{5}{\sqrt{233}}\cr -\frac{1}{3\sqrt{2}}&\frac{53}{3\sqrt{2}\sqrt{233}}&-\frac{8}{\sqrt{233}}\cr \frac{1}{3\sqrt{2}}&\frac{37}{3 \sqrt{2}\sqrt{233}}&\frac{12}{\sqrt{233}}} \right] \)
------------------------
3.上三角矩陣\(R\)
\( R\left[1,1 \right]=|\; y1 |\;=3 \sqrt{2} \)
\( \displaystyle R\left[1,2 \right]=q1^Tx2=\left[ \matrix{\displaystyle -\frac{2^{3/2}}{3}&-\frac{1}{3\sqrt{2}}&\frac{1}{3\sqrt{2}}} \right]\cdot \left[ \matrix{0 \cr 3 \cr 2} \right]=-\frac{1}{3\sqrt{2}} \)
\( R\left[1,3 \right]=q1^Tx3=\left[ \matrix{\displaystyle -\frac{2^{3/2}}{3}&-\frac{1}{3\sqrt{2}}&\frac{1}{3\sqrt{2}}} \right]\cdot \left[ \matrix{1 \cr 1 \cr 5} \right]=0\)
\( \displaystyle R\left[2,2 \right]=|\; y2 |\;=\frac{\sqrt{233}}{3\sqrt{2}} \)
\( \displaystyle R \left[ 2,3 \right]=q2^Tx3=\left[ \matrix{\displaystyle -\frac{2^{3/2}}{3 \sqrt{233}}&\frac{53}{3\sqrt{2}\sqrt{233}}&\frac{37}{3\sqrt{2}\sqrt{233}}} \right]\cdot \left[ \matrix{1 \cr 1 \cr 5} \right]=\frac{39\sqrt{2}}{\sqrt{233}} \)
\( \displaystyle R \left[ 3,3 \right]=|\; y3 |\;=\frac{57}{\sqrt{233}} \)
\( R=\left[ \matrix{\displaystyle 3\sqrt{2}&-\frac{1}{3\sqrt{2}}&0\cr 0&\frac{\sqrt{233}}{3\sqrt{2}}&\frac{39\sqrt{2}}{\sqrt{233}}\cr 0&0&\frac{57}{\sqrt{233}}} \right] \)
(%o3) \( [\; \left[ \matrix{\displaystyle -\frac{2^{3/2}}{3}&-\frac{2^{3/2}}{3\sqrt{233}}&\frac{5}{\sqrt{233}}\cr -\frac{1}{3\sqrt{2}}&\frac{53}{3\sqrt{2}\sqrt{233}}&-\frac{8}{\sqrt{233}}\cr \frac{1}{3\sqrt{2}}&\frac{37}{3 \sqrt{2}\sqrt{233}}&\frac{12}{\sqrt{233}}} \right],\left[ \matrix{\displaystyle 3\sqrt{2}&-\frac{1}{3\sqrt{2}}&0\cr 0&\frac{\sqrt{233}}{3\sqrt{2}}&\frac{39\sqrt{2}}{\sqrt{233}}\cr 0&0&\frac{57}{\sqrt{233}}} \right] ]\; \)

Q矩陣和R矩陣相乘得到矩陣A
(%i4) Q.R;
(%o4) \( \left[ \matrix{-4&0&1\cr -1&3&1 \cr 1&2&5} \right] \)

線代啟示錄文章中的範例
(%i5)
A:matrix([1,2,-1],
              [1,-1,2],
              [-1,1,1],
              [1,-1,2]);

(%o5) \( \left[ \matrix{1&2&-1\cr 1&-1&2\cr -1&1&1\cr 1&-1&2} \right] \)

執行QR分解
(%i6) [Q,R]: QRbyGramSchmidt(A);
1.正交化
\( y1=x1=\left[ \matrix{1\cr 1 \cr -1 \cr 1} \right] \)
\( |\; y1 |\;=sqrt(\left[ \matrix{1\cr 1\cr -1\cr 1} \right]\cdot \left[ \matrix{1\cr 1\cr -1\cr 1} \right])=2 \)
------------------------
\( y2=x2=\left[ \matrix{2\cr -1\cr 1\cr -1} \right] \)
\( y2=y2-(y1^T \cdot x2)/(|\; y1 |\;^2)y1 \)
\( =\left[ \matrix{2\cr -1\cr 1\cr -1} \right]-(\left[ \matrix{1&1&-1&1} \right]\cdot \left[ \matrix{2\cr -1\cr 1\cr -1} \right])/(\left[ \matrix{1\cr 1\cr -1\cr 1} \right]\cdot \left[ \matrix{1\cr 1\cr -1\cr 1} \right])\left[ \matrix{1\cr 1\cr -1\cr 1} \right] \)
\( =\left[ \matrix{2\cr -1\cr 1\cr -1} \right]- -1/4 \left[ \matrix{1 \cr 1\cr -1\cr 1} \right] \)
\( =\left[ \matrix{\displaystyle \frac{9}{4}\cr -\frac{3}{4}\cr \frac{3}{4}\cr -\frac{3}{4}} \right] \)
\( \displaystyle |\; y2 |\;=sqrt(\left[ \matrix{\displaystyle \frac{9}{4}\cr -\frac{3}{4}\cr \frac{3}{4}\cr -\frac{3}{4}} \right] \cdot \left[ \matrix{\displaystyle \frac{9}{4}\cr -\frac{3}{4}\cr \frac{3}{4}\cr -\frac{3}{4}} \right])=\frac{3^{3/2}}{2} \)
------------------------
\( y3=x3=\left[ \matrix{-1 \cr 2 \cr 1 \cr 2} \right] \)
\( y3=y3-(y1^T \cdot x3)/(|\; y1 |\;^2)y1 \)
\( =\left[ \matrix{-1 \cr 2 \cr 1 \cr 2} \right]-(\left[ \matrix{1&1&-1&1} \right]\cdot \left[ \matrix{-1 \cr 2 \cr 1 \cr 2} \right])/(\left[ \matrix{1 \cr 1 \cr -1 \cr 1} \right] \cdot \left[ \matrix{1 \cr 1 \cr -1 \cr 1} \right])\left[ \matrix{1 \cr 1 \cr -1 \cr 1} \right] \)
\( =\left[ \matrix{-1\cr 2 \cr 1 \cr 2} \right]-2/4 \left[ \matrix{1 \cr 1 \cr -1 \cr 1} \right] \)
\( =\left[ \matrix{\displaystyle -\frac{3}{2}\cr \frac{3}{2}\cr \frac{3}{2}\cr \frac{3}{2}} \right] \)
\( y3=y3-(y2^T \cdot x3)/(|\; y2 |\;^2)y2 \)
\( =\left[ \matrix{\displaystyle -\frac{3}{2}\cr \frac{3}{2}\cr \frac{3}{2}\cr \frac{3}{2}} \right]-(\left[ \matrix{\frac{9}{4}&-\frac{3}{4}&\frac{3}{4}&-\frac{3}{4}} \right]\cdot \left[ \matrix{-1\cr 2\cr 1\cr 2} \right])/(\left[ \matrix{\displaystyle \frac{9}{4}\cr -\frac{3}{4}\cr \frac{3}{4}\cr -\frac{3}{4}} \right] \cdot \left[ \matrix{\displaystyle \frac{9}{4}\cr -\frac{3}{4}\cr \frac{3}{4}\cr -\frac{3}{4}} \right])\left[ \matrix{\displaystyle \frac{9}{4}\cr -\frac{3}{4}\cr \frac{3}{4}\cr -\frac{3}{4}} \right] \)
\( \displaystyle =\left[ \matrix{\displaystyle -\frac{3}{2}\cr \frac{3}{2}\cr \frac{3}{2}\cr \frac{3}{2}} \right]--\frac{9}{2}/\frac{27}{4}\left[ \matrix{\displaystyle \frac{9}{4}\cr -\frac{3}{4}\cr \frac{3}{4}\cr -\frac{3}{4}} \right] \)
\( =\left[ \matrix{0\cr 1\cr 2\cr 1} \right] \)
\( |\; y3 |\;=sqrt(\left[ \matrix{0\cr 1\cr 2\cr 1} \right] \cdot \left[ \matrix{0\cr 1\cr 2\cr 1} \right])=\sqrt{6} \)
------------------------
2.單位化得正交矩陣\(Q\)
\( q1=y1/ |\; y1 |\;=\left[ \matrix{1 \cr 1 \cr -1 \cr 1} \right]/2=\left[ \matrix{\displaystyle \frac{1}{2}\cr \frac{1}{2}\cr -\frac{1}{2}\cr\frac{1}{2}} \right] \)
\( \displaystyle q2=y2/ |\; y2 |\;=\left[ \matrix{\displaystyle \frac{9}{4}\cr -\frac{3}{4}\cr\frac{3}{4}\cr -\frac{3}{4}} \right]/ \frac{3^{3/2}}{2}=\left[ \matrix{\displaystyle \frac{\sqrt{3}}{2}\cr -\frac{1}{2\sqrt{3}}\cr\frac{1}{2\sqrt{3}}\cr -\frac{1}{2\sqrt{3}}} \right] \)
\( q3=y3/ |\; y3 |\;=\left[ \matrix{0\cr 1\cr 2\cr 1} \right]/ \sqrt{6}=\left[ \matrix{\displaystyle 0\cr\frac{1}{\sqrt{6}}\cr\frac{2}{\sqrt{6}}\cr\frac{1}{\sqrt{6}}} \right] \)
\( Q=\left[ \matrix{\displaystyle \frac{1}{2}&\frac{\sqrt{3}}{2}&0\cr \frac{1}{2}&-\frac{1}{2\sqrt{3}}&\frac{1}{\sqrt{6}}\cr -\frac{1}{2}&\frac{1}{2\sqrt{3}}&\frac{2}{\sqrt{6}}\cr \frac{1}{2}&-\frac{1}{2\sqrt{3}}&\frac{1}{\sqrt{6}}} \right] \)
------------------------
3.上三角矩陣\(R\)
\( R[1,1]=|\; y1 |\;=2 \)
\( \displaystyle R[1,2]=q1^Tx2=\left[ \matrix{\displaystyle \frac{1}{2}&\frac{1}{2}&-\frac{1}{2}&\frac{1}{2}} \right]\cdot \left[ \matrix{2\cr -1\cr 1\cr -1} \right]=-\frac{1}{2} \)
\( R[1,3]=q1^Tx3=\left[ \matrix{\displaystyle \frac{1}{2}&\frac{1}{2}&-\frac{1}{2}&\frac{1}{2}} \right] \cdot \left[ \matrix{-1\cr 2\cr 1\cr 2} \right]=1 \)
\( \displaystyle R[2,2]=|\; y2 |\;=\frac{3^{3/2}}{2} \)
\( R[2,3]=q2^Tx3=\left[ \matrix{\displaystyle \frac{\sqrt{3}}{2}& -\frac{1}{2 \sqrt{3}}& \frac{1}{2 \sqrt{3}}& -\frac{1}{2\sqrt{3}}} \right] \cdot \left[ \matrix{-1\cr 2 \cr 1 \cr 2} \right]=-\sqrt{3} \)
\( R[3,3]=|\; y3 |\;=\sqrt{6} \)
\( R=\left[ \matrix{\displaystyle 2&-\frac{1}{2}&1\cr 0&\frac{3^{3/2}}{2}&-\sqrt{3} \cr 0&0&\sqrt{6}} \right] \)
(%o6) \( [\left[ \matrix{\displaystyle \frac{1}{2}&\frac{\sqrt{3}}{2}&0\cr \frac{1}{2}&-\frac{1}{2\sqrt{3}}&\frac{1}{\sqrt{6}}\cr -\frac{1}{2}&\frac{1}{2\sqrt{3}}&\frac{2}{\sqrt{6}}\cr \frac{1}{2}&-\frac{1}{2\sqrt{3}}&\frac{1}{\sqrt{6}}} \right],\left[ \matrix{\displaystyle 2&-\frac{1}{2}&1\cr 0&\frac{3^{3/2}}{2}&-\sqrt{3} \cr 0&0&\sqrt{6}} \right]] \)

Q矩陣和R矩陣相乘得到矩陣A
(%i7) Q.R;
(%o7) \( \left[ \matrix{1&2&-1\cr 1&-1&2\cr -1&1&1\cr 1&-1&2} \right] \)

[ 本帖最後由 bugmens 於 2016-10-29 03:37 PM 編輯 ]
作者: bugmens    時間: 2016-10-8 23:11

2.求QR分解的方法--Householder

詳細方法請見線代啟示錄的"Householder 變換於 QR 分解的應用"
https://ccjou.wordpress.com/2011 ... qr-分解的應用/



要先載入diag才能使用diag指令
(%i1) load(diag);
(%o1) C:\maxima-5.38.1\share\maxima\5.38.1_5_gdf93b7b_dirty\share\contrib\diag.mac

執行1/(√3-√2)類型的分母有理化
(%i2) algebraic : true;
(%o2) true

利用Householder的QR分解
(%i3)
QRbyHouseholder(A):=block
([columns:length(A[1]),rows:length(A),rank:rank(A),x,e1,v,Hhat,Ahat,I,H,Q,R,list],
print("1.初始值"),
Q:ident(rows),
R:A,
Ahat:A,
print("Q=",Q,",R=",R,",A^=",Ahat),
for i:1 thru rank do
  (print("2.計算H矩陣"),
   x:col(Ahat,1),
   e1:zeromatrix(rows+1-i,1),
   e1[1,1]:1,
   delta:ratsimp(sqrt(x.x)),
   v:x-delta*e1,
   print("A矩陣第1行x=",x,",第1元為1其餘為0的行向量e1=",e1,",x的長度δ=",delta),
   if sqrt(v.v)=0 then
     (print("x-δe1=",x,"-",delta,"*",e1,"=",v,"為0向量,結束")
     )
   else/*若x-δe1不為0向量*/
     (print("v=(x-δe1)/|x-δe1|=",x,"-",delta,"*",e1,"/sqrt(",v,".",v,")=",v,"/",sqrt(v.v),"=",v:ratsimp(v/sqrt(v.v))),
      I:ident(length(Ahat)),
      Hhat:ratsimp(I-2*v.transpose(v)),
      print("H^=I-2*v.v^T=",I,"-2*",v,".",transpose(v),"=",Hhat),
      if i=1 then/*第一次的H矩陣對角線不用補1,所以H=H^*/
        (H:Hhat,
         print("H矩陣就是H^=",H)
        )
      else/*第2次之後H矩陣對角線要補上1,讓H矩陣大小和A矩陣相同*/
        (H:diag([ident(i-1),Hhat]),
         print("H^對角線補上1形成H矩陣=",H)
        ),   
      print("3,計算Q矩陣和R矩陣"),
      print("Q=Q.H=",Q,".",H,"=",Q:ratsimp(Q.H)),
      print("R=H.R=",H,".",R,"=",R:ratsimp(H.R)),
      if i<rank then   
        (print("取R矩陣右下角得到下一個A^=",Ahat:ratsimp(submatrix(1,Hhat.Ahat,1)))),
      print("------------------------")
     )
  ),
  if rows#rank then/*刪除R矩陣0列和Q矩陣對應的行*/
    (list:makelist(i,i,rank(A)+1,rows),
     print("4-1.刪除Q矩陣的第",list,"行"),
     print(Q,"=>",Q:apply(submatrix,cons(Q,list))),
     print("4-2,刪除R矩陣的第",list,"列"),
     print(R,"=>",R:apply(submatrix,reverse(cons(R,list))))
    ),
  for i:1 thru rank do/*讓R矩陣對角線元素為正*/
    (if R[i,i]<0 then
       (print("R[",i,",",i,"]=",R[i,i],"<0"),
        print("R矩陣第",i,"列",R[ i ],",Q矩陣第",i,"行",col(Q,i),"乘上-1倍"),
        R[ i ]:-R[ i ],
        Q:columnop(Q,i,i,2)
       )
    ),
  return([Q,R])
)$


線代啟示錄文章中的範例
(%i4)
A:matrix([0,-15,14],
              [4,32,2],
              [3,-1,4]);

(%o4) \( \left[ \matrix{0&-15&14\cr 4&32&2\cr 3&-1&4} \right] \)

執行QR分解               
(%i5) [Q,R]: QRbyHouseholder(A);
1.初始值
\( Q=\left[ \matrix{1&0&0\cr 0&1&0 \cr 0&0&1} \right],R=\left[ \matrix{0&-15&14\cr 4&32&2 \cr 3&-1&4} \right],\hat{A}=\left[ \matrix{0&-15&14\cr 4&32&2 \cr 3&-1&4} \right] \)
2.計算\(H\)矩陣
\(A\)矩陣第1行\(x=\left[ \matrix{0\cr 4\cr 3}\right] \),第1元為1其餘為0的行向量\( e1=\left[ \matrix{1 \cr 0\cr 0} \right],x \)的長度\( \delta=5 \)
\( v=(x-\delta e1)/|\; x-\delta e1 |\;=\left[ \matrix{\displaystyle 0\cr 4\cr 3} \right]-5 *\left[ \matrix{1\cr 0\cr 0} \right]/sqrt(\left[ \matrix{-5\cr 4\cr 3} \right] \cdot \left[ \matrix{-5\cr 4\cr 3} \right])=\left[ \matrix{-5\cr 4\cr 3} \right]/5\sqrt{2}=\left[ \matrix{\displaystyle -\frac{1}{\sqrt{2}}\cr \frac{2^{3/2}}{5}\cr \frac{3}{5\sqrt{2}}} \right] \)
\( \hat{H}=I-2*v \cdot v^T=\left[ \matrix{1&0&0\cr 0&1&0 \cr 0&0&1} \right]-2*\left[ \matrix{\displaystyle -\frac{1}{\sqrt{2}}\cr \frac{2^{3/2}}{5}\cr \frac{3}{5\sqrt{2}}}\right] \cdot \left[ \matrix{\displaystyle -\frac{1}{\sqrt{2}}&\frac{2^{3/2}}{5}&\frac{3}{5\sqrt{2}}} \right]=\left[ \matrix{\displaystyle 0&\frac{4}{5}&\frac{3}{5}\cr \frac{4}{5}&\frac{9}{25}&-\frac{12}{25} \cr \frac{3}{5}&-\frac{12}{25}&\frac{16}{25}} \right] \)
\(H\)矩陣就是\( \hat{H}=\left[\matrix{\displaystyle 0&\frac{4}{5}&\frac{3}{5}\cr \frac{4}{5}&\frac{9}{25}&-\frac{12}{25} \cr \frac{3}{5}&-\frac{12}{25}&\frac{16}{25}} \right] \)
3.計算\(Q\)矩陣和\(R\)矩陣
\( Q=Q \cdot H=\left[ \matrix{1&0&0\cr 0&1&0 \cr 0&0&1} \right]\cdot \left[ \matrix{\displaystyle 0&\frac{4}{5}&\frac{3}{5}\cr \frac{4}{5}&\frac{9}{25}&-\frac{12}{25} \cr \frac{3}{5}&-\frac{12}{25}&\frac{16}{25}} \right]=\left[ \matrix{\displaystyle 0&\frac{4}{5}&\frac{3}{5}\cr \frac{4}{5}&\frac{9}{25}&-\frac{12}{25} \cr \frac{3}{5}&-\frac{12}{25}&\frac{16}{25}} \right] \)
\( R=H \cdot R=\left[ \matrix{\displaystyle 0&\frac{4}{5}&\frac{3}{5}\cr \frac{4}{5}&\frac{9}{25}&-\frac{12}{25} \cr \frac{3}{5}&-\frac{12}{25}&\frac{16}{25}} \right] \cdot \left[ \matrix{0&-15&14\cr 4&32&2 \cr 3&-1&4}\right]=\left[ \matrix{5&25&4\cr 0&0&10 \cr 0&-25&10} \right] \)
取\(R\)矩陣右下角得到下一個\( \hat{A}=\left[ \matrix{0&10\cr -25&10\cr } \right]\)
------------------------
2.計算\(H\)矩陣
\( A矩陣第1行x=\left[ \matrix{0 \cr -25} \right] \),第1元為1其餘為0的行向量\( e1=\left[ \matrix{1 \cr 0} \right],x \)的長度\( \delta=25 \)
\( v=(x-\delta e1)/|\; x-\delta e1 |\;=\left[ \matrix{0 \cr -25} \right]-25*\left[ \matrix{1 \cr 0} \right]/sqrt(\left[ \matrix{-25 \cr -25} \right] \cdot \left[ \matrix{-25 \cr -25} \right])=\left[ \matrix{-25 \cr -25} \right]/25 \sqrt{2}=\left[ \matrix{\displaystyle -\frac{1}{\sqrt{2}} \cr -\frac{1}{\sqrt{2}}} \right] \)
\( \hat{H}=I-2*v \cdot v^T=\left[ \matrix{1&0\cr 0&1} \right]-2*\left[ \matrix{\displaystyle -\frac{1}{\sqrt{2}}\cr -\frac{1}{\sqrt{2}}} \right] \cdot \left[ \matrix{\displaystyle -\frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}}} \right]=\left[ \matrix{0&-1 \cr -1&0} \right] \)
\( \hat{H} \)對角線補上1形成\(H\)矩陣\( =\left[ \matrix{1&0&0\cr 0&0&-1 \cr 0&-1&0} \right] \)
3.計算\(Q\)矩陣和\(R\)矩陣
\( Q=Q \cdot H=\left[ \matrix{\displaystyle 0&\frac{4}{5}&\frac{3}{5}\cr \frac{4}{5}&\frac{9}{25}&-\frac{12}{25} \cr \frac{3}{5}&-\frac{12}{25}&\frac{16}{25}} \right]\cdot \left[ \matrix{1&0&0\cr 0&0&-1 \cr 0&-1&0} \right]=\left[ \matrix{\displaystyle 0&-\frac{3}{5}&-\frac{4}{5}\cr \frac{4}{5}&\frac{12}{25}&-\frac{9}{25} \cr \frac{3}{5}&-\frac{16}{25}&\frac{12}{25}} \right] \)
\( R=H \cdot R=\left[ \matrix{1&0&0\cr 0&0&-1 \cr 0&-1&0} \right]\cdot \left[ \matrix{5&25&4\cr 0&0&10 \cr 0&-25&10} \right]=\left[ \matrix{5&25&4\cr 0&25&-10 \cr 0&0&-10} \right] \)
取\(R\)矩陣右下角得到下一個\( \hat{A}=[-10] \)
------------------------
2.計算\(H\)矩陣
\(A\)矩陣第1行\(x=\left[ -10 \right]\),第1元為1其餘為0的行向量\(e1=\left[ 1 \right]\),\(x\)的長度\( \delta=10 \)
\( v=(x-\delta e1)/ |\; x-\delta e1 |\;=[-10]-10*[1]/sqrt([-20]\cdot [-20])=[-20]/20=[-1] \)
\( \hat{H}=I-2*v \cdot v^T=[1]-2*[-1]\cdot [-1]=[-1] \)
\( \hat{H} \)對角線補上1形成\(H\)矩陣\(=\left[ \matrix{1&0&0 \cr 0&1&0 \cr 0&0&-1} \right]\)
3.計算\(Q\)矩陣和\(R\)矩陣
\( Q=Q \cdot H=\left[ \matrix{\displaystyle 0&-\frac{3}{5}&-\frac{4}{5} \cr \frac{4}{5}&\frac{12}{25}&-\frac{9}{25} \cr\frac{3}{5}&-\frac{16}{25}&\frac{12}{25}} \right] \cdot \left[ \matrix{1&0&0\cr 0&1&0\cr 0&0&-1} \right]=\left[ \matrix{\displaystyle 0&-\frac{3}{5}&\frac{4}{5} \cr \frac{4}{5}&\frac{12}{25}&\frac{9}{25} \cr\frac{3}{5}&-\frac{16}{25}&-\frac{12}{25}} \right] \)
\( R=H \cdot R=\left[ \matrix{1&0&0\cr 0&1&0\cr 0&0&-1} \right] \cdot \left[ \matrix{5&25&4\cr 0&25&-10\cr 0&0&-10} \right]=\left[ \matrix{5&25&4\cr 0&25&-10\cr 0&0&10} \right] \)
(%o5) \( [\; \left[ \matrix{\displaystyle 0&-\frac{3}{5}&\frac{4}{5}\cr \frac{4}{5}&\frac{12}{25}&\frac{9}{25} \cr \frac{3}{5}&-\frac{16}{25}&-\frac{12}{25}} \right],\left[ \matrix{5&25&4\cr 0&25&-10 \cr 0&0&10} \right] ]\; \)

Q矩陣和R矩陣相乘得到矩陣A
(%i6) ratsimp(Q.R);
(%o6) \( \left[ \matrix{0&-15&14\cr 4&32&2 \cr 3&-1&4} \right] \)

和前一篇用GramSchmidt方法比較
(%i7)
A:matrix([1,2,-1],
              [1,-1,2],
              [-1,1,1],
              [1,-1,2]);

(%o7) \( \left[ \matrix{1&2&-1\cr 1&-1&2 \cr -1&1&1 \cr 1&-1&2} \right] \)

執行QR分解
(%i8) [Q,R]: QRbyHouseholder(A);
………
計算過程太繁雜不列出來,請自行執行程式。
………

(%o8) \( [\; \left[ \matrix{\displaystyle \frac{1}{2}&\frac{\sqrt{3}}{2}&0 \cr
\frac{1}{2}&-\frac{1}{2\sqrt{3}}&\frac{1}{\sqrt{6}}\cr
-\frac{1}{2}&\frac{1}{2\sqrt{3}}&\frac{\sqrt{6}}{3} \cr
\frac{1}{2}&-\frac{1}{2\sqrt{3}}&\frac{1}{\sqrt{6}} } \right],
\left[ \matrix{\displaystyle 2&-\frac{1}{2}&1\cr 0&\frac{3^{3/2}}{2}&-\sqrt{3}\cr 0&0&\sqrt{6}} \right] ]\; \)

Q矩陣和R矩陣相乘得到矩陣A
(%i9) ratsimp(Q.R);
(%o9) \( \left[ \matrix{1&2&-1\cr 1&-1&2\cr -1&1&1\cr 1&-1&2} \right] \)



註:若沒加algebraic : true,執行結果為
\( [ \left[ \matrix{\displaystyle
\frac{1}{2}&-\frac{3^{3/2}-3}{2\sqrt{3}-6}&0 \cr
\frac{1}{2}&\frac{\sqrt{3}-1}{2 \sqrt{3}-6}&-\frac{(313496\sqrt{3}-542991)\sqrt{6}-767444 \sqrt{3}+1329252}{(767444\sqrt{3}-1329252)\sqrt{6}-626992 \cdot 3^{3/2}+3257946} \cr
-\frac{1}{2}&-\frac{\sqrt{3}-1}{2\sqrt{3}-6}&-\frac{(22508\sqrt{3}-38985)\sqrt{6}-55100\sqrt{3}+95436}{(27550\sqrt{3}-47718)\sqrt{6}-22508 \cdot 3^{3/2}+116955} \cr
\frac{1}{2}&\frac{\sqrt{3}-1}{2\sqrt{3}-6}&-\frac{(22508\sqrt{3}-38985)\sqrt{6}-55100\sqrt{3}+95436}{(55100\sqrt{3}-95436)\sqrt{6}-45016 \cdot 3^{3/2}+233910}} \right], \)
\( \left[ \matrix{\displaystyle 2&-\frac{1}{2}&1 \cr 0&-\frac{3^{5/2}-9}{2\sqrt{3}-6}&\frac{3^{3/2}-3}{\sqrt{3}-3} \cr 0&0&-\frac{(22508 \cdot 3^{3/2}-116955)\sqrt{6}-55100 \cdot 3^{3/2}+286308}{(27550\sqrt{3}-47718)\sqrt{6}-22508 \cdot 3^{3/2}+116955}} \right] \)

[ 本帖最後由 bugmens 於 2016-11-12 08:43 AM 編輯 ]
作者: bugmens    時間: 2016-10-10 18:38

3.求QR分解的方法--GivensRotation

詳細方法請見線代啟示錄的"Givens 旋轉於 QR 分解的應用"
https://ccjou.wordpress.com/2010 ... qr-分解的應用/



利用Givens旋轉的QR分解
(%i1)
QRbyGivensRotation(A):=block
([columns:length(A[1]),rows:length(A),rank:rank(A),n,r,c,s,Qji,list,Q,R],
n:max(columns,rows),
print("1.初始值"),
Q:ident(rows),
R:copymatrix(A),
print("Q=",Q,",R=",R),
for i:1 thru min(columns,rows) do
    (for j:i+1 thru rows do
       (if R[j,i]#0 then/*若R[j,i]已經是0就不處理*/
          (print("2.以R[",i,",",i,"]=",R[i,i],"消掉R[",j,",",i,"]=",R[j,i],"為0的旋轉矩陣Q",j,i),
           print("r=sqrt(R[",i,",",i,"]^2+R[",j,",",i,"]^2)=",r:sqrt(R[i,i]^2+R[j,i]^2)),
           c:R[i,i]/r,
           s:-R[j,i]/r,
           print("cosφ=R[",i,",",i,"]/r=",R[i,i],"/",r,"=",c),
           print("sinφ=-R[",j,",",i,"]/r=-",R[j,i],"/",r,"=",s),
           Qji:ident(rows),
           Qji[i,i]:c,    Qji[i,j]:-s,
           Qji[j,i]:s,    Qji[j,j]:c,
           print("Q",j,i,"[",i,",",i,"]=cosφ=",c,",    Q",j,i,"[",i,",",j,"]=-sinφ=",-s),
           print("Q",j,i,"[",j,",",i,"]=sinφ=",s,",     Q",j,i,"[",j,",",j,"]=cosφ=",c),
           print("旋轉矩陣Q",j,i,"=",Qji),
           print("3.更新Q矩陣和R矩陣"),
           print("Q=Q.Q",j,i,"^T=",Q,".",transpose(Qji),"=",Q:ratsimp(Q.transpose(Qji))),
           print("R=Q",j,i,".R=",Qji,".",R,"=",R:ratsimp(Qji.R)),
           print("----------------------")
          )
       )
    ),
if rows#rank then
  (list:makelist(i,i,rank(A)+1,rows),
   print("4-1.刪除Q矩陣的第",list,"行"),
   print(Q,"=>",Q:apply(submatrix,cons(Q,list))),
   print("4-2.刪除R矩陣的第",list,"列"),
   print(R,"=>",R:apply(submatrix,reverse(cons(R,list))))
  ),
for i:1 thru rank do
  (if R[i,i]<0 then
    (print("R[",i,",",i,"]=",R[i,i],"<0"),
     print("R矩陣第",i,"列",R[ i ],",Q矩陣第",i,"行",col(Q,i),"乘上-1倍"),
     R[ i ]:-R[ i ],
     Q:columnop(Q,i,i,2)
    )
  ),
return(ratsimp([Q,R]))
)$


線代啟示錄文章中的範例
(%i2)
A:matrix([0,-15,14],
              [4,32,2],
              [3,-1,4]);

(%o2) \( \left[ \matrix{0&-15&14\cr 4&32&2\cr 3&-1&4} \right] \)

執行QR分解
(%i3) [Q,R]: QRbyGivensRotation(A);
1.初始值
\( Q=\left[ \matrix{1&0&0\cr 0&1&0\cr 0&0&1} \right],R=\left[ \matrix{0&-15&14\cr 4&32&2\cr 3&-1&4} \right] \)
2.以\(R[1,1]=0\)消掉\(R[2,1]=4\)為0的旋轉矩陣\(Q21\)
\(r=sqrt(R[1,1]^2+R[2,1]^2)=4\)
\(cosφ=R[1,1]/r=0/4=0\)
\(sinφ=-R[2,1]/r=-4/4=-1\)
\(Q21[1,1]=cosφ=0\),    \(Q21[1,2]=-sinφ=1\)
\(Q21[2,1]=sinφ=-1\),   \(Q21[2,2]=cosφ=0\)
旋轉矩陣\(Q21=\left[ \matrix{0&1&0\cr -1&0&0\cr 0&0&1} \right]\)
3.更新\(Q\)矩陣和\(R\)矩陣
\( Q=Q \cdot Q21^T=\left[ \matrix{1&0&0\cr 0&1&0\cr 0&0&1} \right] \cdot \left[ \matrix{0&-1&0\cr 1&0&0\cr 0&0&1} \right]=\left[ \matrix{0&-1&0\cr 1&0&0\cr 0&0&1} \right] \)
\( R=Q21 \cdot R=\left[ \matrix{0&1&0\cr -1&0&0\cr 0&0&1} \right] \cdot \left[ \matrix{0&-15&14\cr 4&32&2\cr 3&-1&4} \right]=\left[ \matrix{4&32&2\cr 0&15&-14\cr 3&-1&4} \right] \)
------------------------
2.以\(R[1,1]=4\)消掉\(R[3,1]=3\)為0的旋轉矩陣\(Q31\)
\(r=sqrt(R[1,1]^2+R[3,1]^2)=5\)
\( \displaystyle cosφ=R[1,1]/r=4/5=\frac{4}{5} \)
\( \displaystyle sinφ=-R[3,1]/r=-3/5=-\frac{3}{5} \)
\( \displaystyle Q31[1,1]=cosφ=\frac{4}{5}\),    \(\displaystyle Q31[1,3]=-sinφ=\frac{3}{5}\)
\( \displaystyle Q31[3,1]=sinφ=-\frac{3}{5}\),     \(\displaystyle Q31[3,3]=cosφ=\frac{4}{5} \)
旋轉矩陣\(Q31=\left[\matrix{\displaystyle \frac{4}{5}&0&\frac{3}{5}\cr 0&1&0\cr -\frac{3}{5}&0&\frac{4}{5}} \right]\)
3.更新\(Q\)矩陣和\(R\)矩陣
\( Q=Q \cdot Q31^T=\left[ \matrix{0&-1&0\cr 1&0&0\cr 0&0&1} \right] \cdot \left[ \matrix{\displaystyle \frac{4}{5}&0&-\frac{3}{5}\cr 0&1&0\cr \frac{3}{5}&0&\frac{4}{5}} \right]=\left[ \matrix{\displaystyle 0&-1&0\cr \frac{4}{5}&0&-\frac{3}{5}\cr \frac{3}{5}&0&\frac{4}{5}} \right] \)
\(R=Q31 \cdot R=\left[ \matrix{\displaystyle \frac{4}{5}&0&\frac{3}{5}\cr 0&1&0\cr -\frac{3}{5}&0&\frac{4}{5}} \right] \cdot \left[ \matrix{4&32&2\cr 0&15&-14\cr 3&-1&4} \right]=\left[ \matrix{5&25&4\cr 0&15&-14\cr 0&-20&2} \right] \)
------------------------
2.以\(R[2,2]=15\)消掉\(R[3,2]=-20\)為0的旋轉矩陣\(Q32\)
\(r=sqrt(R[2,2]^2+R[3,2]^2)=25\)
\( \displaystyle cosφ=R[2,2]/r=15/25=\frac{3}{5} \)
\( \displaystyle sinφ=-R[3,2]/r=--20/25=\frac{4}{5} \)
\( \displaystyle Q32[2,2]=cosφ=\frac{3}{5}\),    \( \displaystyle Q32[2,3]=-sinφ=-\frac{4}{5} \)
\( \displaystyle Q32[3,2]=sinφ=\frac{4}{5}\),    \( \displaystyle Q32[3,3]=cosφ=\frac{3}{5} \)
旋轉矩陣\( Q32=\left[ \matrix{\displaystyle 1&0&0\cr 0&\frac{3}{5}&-\frac{4}{5}\cr 0&\frac{4}{5}&\frac{3}{5}} \right] \)
3.更新\(Q\)矩陣和\(R\)矩陣
\(Q=Q \cdot Q32^T=\left[ \matrix{\displaystyle 0&-1&0\cr \frac{4}{5}&0&-\frac{3}{5}\cr \frac{3}{5}&0&\frac{4}{5}} \right]\cdot \left[ \matrix{\displaystyle 1&0&0\cr 0&\frac{3}{5}&\frac{4}{5}\cr 0&-\frac{4}{5}&\frac{3}{5}} \right]=\left[ \matrix{\displaystyle 0&-\frac{3}{5}&-\frac{4}{5}\cr \frac{4}{5}&\frac{12}{25}&-\frac{9}{25}\cr \frac{3}{5}&-\frac{16}{25}&\frac{12}{25}} \right] \)
\( R=Q32 \cdot R=\left[ \matrix{\displaystyle 1&0&0\cr 0&\frac{3}{5}&-\frac{4}{5}\cr 0&\frac{4}{5}&\frac{3}{5}} \right] \cdot \left[ \matrix{5&25&4\cr 0&15&-14\cr 0&-20&2} \right]=\left[ \matrix{5&25&4\cr 0&25&-10\cr 0&0&-10} \right] \)
------------------------
\( R[3,3]=-10<0 \)
\(R\)矩陣第3列\( [\matrix{0,0,-10}] \),\(Q\)矩陣第3行\( \left[ \matrix{\displaystyle -\frac{4}{5}\cr -\frac{9}{25}\cr \frac{12}{25}} \right] \)乘上\(-1\)倍
0 errors, 0 warnings
(%o3) \( [\; \left[ \matrix{\displaystyle 0&-\frac{3}{5}&\frac{4}{5}\cr \frac{4}{5}&\frac{12}{25}&\frac{9}{25}\cr \frac{3}{5}&-\frac{16}{25}&-\frac{12}{25}} \right],\left[ \matrix{5&25&4\cr 0&25&-10\cr 0&0&10} \right] ]\; \)

Q矩陣和R矩陣相乘得到矩陣A
(%i4) Q.R;
(%o4) \( \left[ \matrix{0&-15&14\cr 4&32&2\cr 3&-1&4} \right] \)

和前一篇用GramSchmidt方法比較
(%i5)
A:matrix([1,2,-1],
              [1,-1,2],
              [-1,1,1],
              [1,-1,2]);

(%o5) \( \left[ \matrix{1&2&-1\cr 1&-1&2\cr -1&1&1 \cr 1&-1&2} \right] \)

執行QR分解
(%i6) [Q,R]: QRbyGivensRotation(A);
1.初始值
\( Q=\left[ \matrix{1&0&0&0\cr 0&1&0&0\cr 0&0&1&0\cr 0&0&0&1} \right],R=\left[ \matrix{1&2&-1\cr 1&-1&2\cr -1&1&1 \cr 1&-1&2} \right] \)
2.以\(R[1,1]=1\)消掉\(R[2,1]=1\)為0的旋轉矩陣\(Q21\)
\( r=sqrt(R[1,1]^2+R[2,1]^2)=\sqrt{2} \)
\( \displaystyle \displaystyle cosφ=R[1,1]/r=1/ \sqrt{2}=\frac{1}{\sqrt{2}} \)
\( \displaystyle sinφ=-R[2,1]/r=-1/ \sqrt{2}=-\frac{1}{\sqrt{2}} \)
\( \displaystyle Q21[1,1]=cosφ=\frac{1}{\sqrt{2}}\),    \( \displaystyle Q21[1,2]=-sinφ=\frac{1}{\sqrt{2}} \)
\( \displaystyle Q21[2,1]=sinφ=-\frac{1}{\sqrt{2}}\),    \( \displaystyle Q21[2,2]=cosφ=\frac{1}{\sqrt{2}} \)
旋轉矩陣\( Q21=\left[ \matrix{\displaystyle \frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}&0&0\cr -\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}&0&0\cr 0&0&1&0 \cr 0&0&0&1} \right] \)
3.更新\(Q\)矩陣和\(R\)矩陣
\( Q=Q \cdot Q21^T=\left[ \matrix{1&0&0&0\cr 0&1&0&0\cr 0&0&1&0 \cr 0&0&0&1} \right] \cdot \left[ \matrix{\displaystyle \frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{2}}&0&0\cr \frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}&0&0\cr 0&0&1&0 \cr 0&0&0&1} \right]=\left[ \matrix{\displaystyle \frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{2}}&0&0\cr \frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}&0&0\cr 0&0&1&0 \cr 0&0&0&1} \right] \)
\( R=Q21 \cdot R=\left[ \matrix{\displaystyle \frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}&0&0\cr -\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}&0&0\cr 0&0&1&0 \cr 0&0&0&1} \right] \cdot \left[ \matrix{1&2&-1\cr 1&-1&2\cr -1&1&1 \cr 1&-1&2} \right]=\left[ \matrix{\displaystyle \sqrt{2}&\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}\cr 0&-\frac{3}{\sqrt{2}}&\frac{3}{\sqrt{2}}\cr -1&1&1 \cr 1&-1&2} \right] \)
------------------------
2.以\(R[1,1]=\sqrt{2}\)消掉\(R[3,1]=-1\)為0的旋轉矩陣\(Q31\)
\( r=sqrt(R[1,1]^2+R[3,1]^2)=\sqrt{3} \)
\( \displaystyle cosφ=R[1,1]/r=\sqrt{2}/\sqrt{3}=\frac{\sqrt{2}}{\sqrt{3}} \)
\( \displaystyle sinφ=-R[3,1]/r=--1/ \sqrt{3}=\frac{1}{\sqrt{3}} \)
\( \displaystyle Q31[1,1]=cosφ=\frac{\sqrt{2}}{\sqrt{3}}\),    \( \displaystyle Q31[1,3]=-sinφ=-\frac{1}{\sqrt{3}} \)
\( \displaystyle Q31[3,1]=sinφ=\frac{1}{\sqrt{3}}\),     \( \displaystyle Q31[3,3]=cosφ=\frac{\sqrt{2}}{\sqrt{3}} \)
旋轉矩陣\(Q31=\left[ \matrix{\displaystyle \frac{\sqrt{2}}{\sqrt{3}}&0&-\frac{1}{\sqrt{3}}&0\cr 0&1&0&0\cr \frac{1}{\sqrt{3}}&0&\frac{\sqrt{2}}{\sqrt{3}}&0\cr 0&0&0&1} \right] \)
3.更新\(Q\)矩陣和\(R\)矩陣
\( Q=Q.Q31^T=\left[ \matrix{\displaystyle \frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{2}}&0&0\cr \frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}&0&0\cr 0&0&1&0\cr 0&0&0&1} \right] \cdot \left[ \matrix{\displaystyle \frac{\sqrt{2}}{\sqrt{3}}&0&\frac{1}{\sqrt{3}}&0\cr 0&1&0&0 \cr -\frac{1}{\sqrt{3}}&0&\frac{\sqrt{2}}{\sqrt{3}}&0 \cr 0&0&0&1} \right]=\left[ \matrix{\displaystyle \frac{1}{\sqrt{3}}&-\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}\sqrt{3}}&0\cr \frac{1}{\sqrt{3}}&\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}\sqrt{3}}&0\cr -\frac{1}{\sqrt{3}}&0&\frac{\sqrt{2}}{\sqrt{3}}&0\cr 0&0&0&1} \right]\)
\( R=Q31 \cdot R=\left[ \matrix{\displaystyle \frac{\sqrt{2}}{\sqrt{3}}&0&-\frac{1}{\sqrt{3}}&0\cr 0&1&0&0\cr \frac{1}{\sqrt{3}}&0&\frac{\sqrt{2}}{\sqrt{3}}&0\cr 0&0&0&1} \right] \cdot \left[ \matrix{\displaystyle \sqrt{2}&\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}\cr 0&-\frac{3}{\sqrt{2}}&\frac{3}{\sqrt{2}}\cr -1&1&1\cr 1&-1&2} \right]=\left[ \matrix{\displaystyle \sqrt{3}&0&0\cr 0&-\frac{3}{\sqrt{2}}&\frac{3}{\sqrt{2}}\cr 0&\frac{\sqrt{3}}{\sqrt{2}}&\frac{\sqrt{3}}{\sqrt{2}}\cr 1&-1&2} \right] \)
------------------------
2.以\(R[1,1]=\sqrt{3}\)消掉\(R[4,1]=1\)為0的旋轉矩陣\(Q41\)
\( r=sqrt(R[1,1]^2+R[4,1]^2)=2 \)
\( \displaystyle cosφ=R[1,1]/r=\sqrt{3}/2=\frac{\sqrt{3}}{2} \)
\( \displaystyle sinφ=-R[4,1]/r=-1/2=-\frac{1}{2} \)
\( \displaystyle Q41[1,1]=cosφ=\frac{\sqrt{3}}{2}\),    \( \displaystyle Q41[1,4]=-sinφ=\frac{1}{2} \)
\( \displaystyle Q41[4,1]=sinφ=-\frac{1}{2}\),     \( \displaystyle Q41[4,4]=cosφ=\frac{\sqrt{3}}{2} \)
旋轉矩陣\(Q41=\left[ \matrix{\displaystyle \frac{\sqrt{3}}{2}&0&0&\frac{1}{2}\cr 0&1&0&0\cr 0&0&1&0\cr -\frac{1}{2}&0&0&\frac{\sqrt{3}}{2}} \right]\)
3.更新\(Q\)矩陣和\(R\)矩陣
\( Q=Q \cdot Q41^T=\left[ \matrix{\displaystyle \frac{1}{\sqrt{3}}&-\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}\sqrt{3}}&0\cr \frac{1}{\sqrt{3}}&\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}\sqrt{3}}&0\cr -\frac{1}{\sqrt{3}}&0&\frac{\sqrt{2}}{\sqrt{3}}&0\cr 0&0&0&1} \right] \cdot \left[ \matrix{\displaystyle \frac{\sqrt{3}}{2}&0&0&-\frac{1}{2}\cr 0&1&0&0\cr 0&0&1&0\cr \frac{1}{2}&0&0&\frac{\sqrt{3}}{2}} \right]=\left[ \matrix{\displaystyle \frac{1}{2}&-\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}\sqrt{3}}&-\frac{1}{2\sqrt{3}}\cr \frac{1}{2}&\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}\sqrt{3}}&-\frac{1}{2\sqrt{3}}\cr -\frac{1}{2}&0&\frac{\sqrt{2}}{\sqrt{3}}&\frac{1}{2\sqrt{3}}\cr \frac{1}{2}&0&0&\frac{\sqrt{3}}{2}} \right]\)
\( R=Q41 \cdot R=\left[ \matrix{\displaystyle \frac{\sqrt{3}}{2}&0&0&\frac{1}{2}\cr 0&1&0&0\cr 0&0&1&0\cr -\frac{1}{2}&0&0&\frac{\sqrt{3}}{2}} \right] \cdot \left[ \matrix{\displaystyle \sqrt{3}&0&0\cr 0&-\frac{3}{\sqrt{2}}&\frac{3}{\sqrt{2}}\cr 0&\frac{\sqrt{3}}{\sqrt{2}}&\frac{\sqrt{3}}{\sqrt{2}}\cr 1&-1&2} \right]=\left[ \matrix{\displaystyle 2&-\frac{1}{2}&1\cr 0&-\frac{3}{\sqrt{2}}&\frac{3}{\sqrt{2}}\cr 0&\frac{\sqrt{3}}{\sqrt{2}}&\frac{\sqrt{3}}{\sqrt{2}}\cr 0&-\frac{\sqrt{3}}{2}&\sqrt{3}} \right] \)
------------------------
2.以\( \displaystyle R[2,2]=-\frac{3}{\sqrt{2}}\)消掉\(\displaystyle R[3,2]=\frac{\sqrt{3}}{\sqrt{2}} \)為0的旋轉矩陣\(Q32\)
\( r=sqrt(R[2,2]^2+R[3,2]^2)=\sqrt{6} \)
\( \displaystyle cosφ=R[2,2]/r=-\frac{3}{\sqrt{2}}/ \sqrt{6}=-\frac{3}{\sqrt{2}\sqrt{6}} \)
\( \displaystyle sinφ=-R[3,2]/r=-\frac{\sqrt{3}}{\sqrt{2}}/ \sqrt{6}=-\frac{\sqrt{3}}{\sqrt{2}\sqrt{6}} \)
\( \displaystyle Q32[2,2]=cosφ=-\frac{3}{\sqrt{2}\sqrt{6}}\),    \(\displaystyle Q32[2,3]=-sinφ=\frac{\sqrt{3}}{\sqrt{2}\sqrt{6}} \)
\( \displaystyle Q32[3,2]=sinφ=-\frac{\sqrt{3}}{\sqrt{2}\sqrt{6}}\),     \( \displaystyle Q32[3,3]=cosφ=-\frac{3}{\sqrt{2}\sqrt{6}} \)
旋轉矩陣\( Q32=\left[ \matrix{\displaystyle 1&0&0&0\cr 0&-\frac{3}{\sqrt{2}\sqrt{6}}&\frac{\sqrt{3}}{\sqrt{2}\sqrt{6}}&0\cr 0&-\frac{\sqrt{3}}{\sqrt{2}\sqrt{6}}&-\frac{3}{\sqrt{2}\sqrt{6}}&0\cr 0&0&0&1} \right] \)
3.更新\(Q\)矩陣和\(R\)矩陣
\( Q=Q \cdot Q32^T=\left[ \matrix{\displaystyle \frac{1}{2}&-\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}\sqrt{3}}&-\frac{1}{2\sqrt{3}}\cr \frac{1}{2}&\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}\sqrt{3}}&-\frac{1}{2\sqrt{3}}\cr -\frac{1}{2}&0&\frac{\sqrt{2}}{\sqrt{3}}&\frac{1}{2\sqrt{3}}\cr \frac{1}{2}&0&0&\frac{\sqrt{3}}{2}} \right] \cdot \left[ \matrix{\displaystyle 1&0&0&0\cr 0&-\frac{3}{\sqrt{2}\sqrt{6}}&-\frac{\sqrt{3}}{\sqrt{2}\sqrt{6}}&0 \cr 0&\frac{\sqrt{3}}{\sqrt{2}\sqrt{6}}&-\frac{3}{\sqrt{2}\sqrt{6}}&0\cr 0&0&0&1} \right]=\left[ \matrix{\displaystyle \frac{1}{2}&\frac{2}{\sqrt{6}}&0&-\frac{1}{2\sqrt{3}}\cr \frac{1}{2}&-\frac{1}{\sqrt{6}}&-\frac{\sqrt{3}}{\sqrt{6}}&-\frac{1}{2\sqrt{3}}\cr -\frac{1}{2}&\frac{1}{\sqrt{6}}&-\frac{\sqrt{3}}{\sqrt{6}}&\frac{1}{2\sqrt{3}}\cr \frac{1}{2}&0&0&\frac{\sqrt{3}}{2}} \right] \)
\( R=Q32 \cdot R=\left[ \matrix{\displaystyle 1&0&0&0\cr 0&-\frac{3}{\sqrt{2}\sqrt{6}}&\frac{\sqrt{3}}{\sqrt{2}\sqrt{6}}&0\cr 0&-\frac{\sqrt{3}}{\sqrt{2}\sqrt{6}}&-\frac{3}{\sqrt{2}\sqrt{6}}&0\cr 0&0&0&1} \right] \cdot \left[ \matrix{\displaystyle 2&-\frac{1}{2}&1\cr 0&-\frac{3}{\sqrt{2}}&\frac{3}{\sqrt{2}}\cr 0&\frac{\sqrt{3}}{\sqrt{2}}&\frac{\sqrt{3}}{\sqrt{2}}\cr 0&-\frac{\sqrt{3}}{2}&\sqrt{3}} \right]=\left[ \matrix{\displaystyle 2&-\frac{1}{2}&1\cr 0&\sqrt{6}&-\frac{3}{\sqrt{6}}\cr 0&0&-\frac{3^{3/2}}{\sqrt{6}}\cr 0&-\frac{\sqrt{3}}{2}&\sqrt{3}} \right] \)
------------------------
2.以\(R[2,2]=\sqrt{6}\)消掉\( \displaystyle R[4,2]=-\frac{\sqrt{3}}{{2}}\)為0的旋轉矩陣\(Q42\)
\( \displaystyle r=sqrt(R[2,2]^2+R[4,2]^2)=\frac{3^{3/2}}{2} \)
\( \displaystyle cosφ=R[2,2]/r=\sqrt{6}/\frac{3^{3/2}}{2}=\frac{2 \sqrt{6}}{3^{3/2}} \)
\( \displaystyle sinφ=-R[4,2]/r=--\frac{\sqrt{3}}{2}/ \frac{3^{3/2}}{2}=\frac{1}{3} \)
\( \displaystyle Q42[2,2]=cosφ=\frac{2 \sqrt{6}}{3^{3/2}}\),    \(\displaystyle Q42[2,4]=-sinφ=-\frac{1}{3} \)
\( \displaystyle Q42[4,2]=sinφ=\frac{1}{3}\),     \( \displaystyle Q42[4,4]=cosφ=\frac{2 \sqrt{6}}{3^{3/2}}\)
旋轉矩陣\( Q42=\left[ \matrix{\displaystyle 1&0&0&0\cr 0&\frac{2\sqrt{6}}{3^{3/2}}&0&-\frac{1}{3}\cr 0&0&1&0\cr 0&\frac{1}{3}&0&\frac{2\sqrt{6}}{3^{3/2}}} \right] \)
3.更新\(Q\)矩陣和\(R\)矩陣
\( Q=Q \cdot Q42^T=\left[ \matrix{\displaystyle \frac{1}{2}&\frac{2}{\sqrt{6}}&0&-\frac{1}{2\sqrt{3}}\cr \frac{1}{2}&-\frac{1}{\sqrt{6}}&-\frac{\sqrt{3}}{\sqrt{6}}&-\frac{1}{2\sqrt{3}}\cr -\frac{1}{2}&\frac{1}{\sqrt{6}}&-\frac{\sqrt{3}}{\sqrt{6}}&\frac{1}{2\sqrt{3}}\cr \frac{1}{2}&0&0&\frac{\sqrt{3}}{2}} \right] \cdot \left[ \matrix{\displaystyle 1&0&0&0\cr 0&\frac{2\sqrt{6}}{3^{3/2}}&0&\frac{1}{3}\cr 0&0&1&0\cr 0&-\frac{1}{3}&0&\frac{2 \sqrt{6}}{3^{3/2}}} \right]=\left[ \matrix{\displaystyle \frac{1}{2}&\frac{\sqrt{3}}{2}&0&0\cr \frac{1}{2}&-\frac{1}{2\sqrt{3}}&-\frac{\sqrt{3}}{\sqrt{6}}&-\frac{1}{\sqrt{6}}\cr -\frac{1}{2}&\frac{1}{2\sqrt{3}}&-\frac{\sqrt{3}}{\sqrt{6}}&\frac{1}{\sqrt{6}}\cr \frac{1}{2}&-\frac{1}{2\sqrt{3}}&0&\frac{\sqrt{6}}{3}} \right] \)
\( R=Q42 \cdot R=\left[ \matrix{\displaystyle 1&0&0&0\cr 0&\frac{2\sqrt{6}}{3^{3/2}}&0&-\frac{1}{3}\cr 0&0&1&0 \cr 0&\frac{1}{3}&0&\frac{2\sqrt{6}}{3^{3/2}}} \right]\cdot \left[ \matrix{\displaystyle 2&-\frac{1}{2}&1\cr 0&\sqrt{6}&-\frac{3}{\sqrt{6}}\cr 0&0&-\frac{3^{3/2}}{\sqrt{6}} \cr 0&-\frac{\sqrt{3}}{2}&\sqrt{3}} \right]=\left[ \matrix{\displaystyle 2&-\frac{1}{2}&1\cr 0&\frac{3^{3/2}}{2}&-\sqrt{3} \cr 0&0&-\frac{3^{3/2}}{\sqrt{6}} \cr 0&0&\frac{3}{\sqrt{6}}} \right] \)
------------------------
2.以\(\displaystyle R[3,3]=-\frac{3^{3/2}}{\sqrt{6}}\)消掉\(\displaystyle R[4,3]=\frac{3}{\sqrt{6}}\)為0的旋轉矩陣\(Q43\)
\(r=sqrt(R[3,3]^2+R[4,3]^2)=\sqrt{6} \)
\( \displaystyle cosφ=R[3,3]/r=-\frac{3^{3/2}}{\sqrt{6}}/ \sqrt{6}=-\frac{3^{3/2}}{6} \)
\( \displaystyle sinφ=-R[4,3]/r=-\frac{3}{\sqrt{6}}/ \sqrt{6}=-\frac{1}{2} \)
\( \displaystyle Q43[3,3]=cosφ=-\frac{3^{3/2}}{6}\),    \( \displaystyle Q43[3,4]=-sinφ=\frac{1}{2}\)
\( \displaystyle Q43[4,3]=sinφ=-\frac{1}{2}\),     \( \displaystyle Q43[4,4]=cosφ=-\frac{3^{3/2}}{6}\)
旋轉矩陣\(Q43=\left[ \matrix{\displaystyle 1&0&0&0\cr 0&1&0&0\cr 0&0&-\frac{3^{3/2}}{6}&\frac{1}{2} \cr 0&0&-\frac{1}{2}&-\frac{3^{3/2}}{6}} \right] \)
3.更新\(Q\)矩陣和\(R\)矩陣
\( Q=Q \cdot Q43^T=\left[ \matrix{\displaystyle \frac{1}{2}&\frac{\sqrt{3}}{2}&0&0\cr \frac{1}{2}&-\frac{1}{2\sqrt{3}}&-\frac{\sqrt{3}}{\sqrt{6}}&-\frac{1}{\sqrt{6}}\cr -\frac{1}{2}&\frac{1}{2\sqrt{3}}&-\frac{\sqrt{3}}{\sqrt{6}}&\frac{1}{\sqrt{6}}\cr \frac{1}{2}&-\frac{1}{2\sqrt{3}}&0&\frac{\sqrt{6}}{3}} \right] \cdot \left[ \matrix{\displaystyle 1&0&0&0\cr 0&1&0&0\cr 0&0&-\frac{3^{3/2}}{6}&-\frac{1}{2}\cr 0&0&\frac{1}{2}&-\frac{3^{3/2}}{6}} \right]=\left[ \matrix{\displaystyle \frac{1}{2}&\frac{\sqrt{3}}{2}&0&0\cr \frac{1}{2}&-\frac{1}{2\sqrt{3}}&\frac{1}{\sqrt{6}}&\frac{\sqrt{6}}{2\sqrt{3}}\cr -\frac{1}{2}&\frac{1}{2\sqrt{3}}&\frac{2}{\sqrt{6}}&0\cr \frac{1}{2}&-\frac{1}{2\sqrt{3}}&\frac{1}{\sqrt{6}}&-\frac{\sqrt{3}}{\sqrt{6}}} \right] \)
\( R=Q43\cdot R=\left[ \matrix{\displaystyle 1&0&0&0\cr 0&1&0&0\cr 0&0&-\frac{3^{3/2}}{6}&\frac{1}{2} \cr 0&0&-\frac{1}{2}&-\frac{3^{3/2}}{6}} \right] \cdot \left[ \matrix{\displaystyle 2&-\frac{1}{2}&1\cr 0&\frac{3^{3/2}}{2}&-\sqrt{3} \cr 0&0&-\frac{3^{3/2}}{\sqrt{6}} \cr 0&0&\frac{3}{\sqrt{6}}} \right]=\left[ \matrix{\displaystyle 2&-\frac{1}{2}&1 \cr 0&\frac{3^{3/2}}{2}&-\sqrt{3} \cr 0&0&\sqrt{6} \cr 0&0&0} \right] \)
------------------------
4-1.刪除\(Q\)矩陣的第[4]行
\( \left[ \matrix{\displaystyle
\frac{1}{2}&\frac{\sqrt{3}}{2}&0&0\cr
\frac{1}{2}&-\frac{1}{2 \sqrt{3}}&\frac{1}{\sqrt{6}}&\frac{\sqrt{6}}{2 \sqrt{3}}\cr
-\frac{1}{2}&\frac{1}{2\sqrt{3}}&\frac{2}{\sqrt{6}}&0\cr
\frac{1}{2}&-\frac{1}{2\sqrt{3}}&\frac{1}{\sqrt{6}}&-\frac{\sqrt{3}}{\sqrt{6}}} \right] \)=>\( \left[ \matrix{\displaystyle
\frac{1}{2}&\frac{\sqrt{3}}{2}& 0\cr
\frac{1}{2}&-\frac{1}{2\sqrt{3}}&\frac{1}{\sqrt{6}}\cr
-\frac{1}{2}&\frac{1}{2\sqrt{3}}&\frac{2}{\sqrt{6}}\cr
\frac{1}{2}&-\frac{1}{2\sqrt{3}}&\frac{1}{\sqrt{6}} } \right] \)
4-2.刪除\(R\)矩陣的第[4]列
\( \left[ \matrix{\displaystyle 2&-\frac{1}{2}&1\cr 0&\frac{3^{3/2}}{2}&-\sqrt{3}\cr 0&0&\sqrt{6} \cr 0&0&0} \right] \)=>\( \left[ \matrix{\displaystyle 2&-\frac{1}{2}&1\cr 0&\frac{3^{3/2}}{2}&-\sqrt{3}\cr 0&0&\sqrt{6}} \right] \)
(%o6) \( [\; \left[ \matrix{\displaystyle \frac{1}{2}&\frac{\sqrt{3}}{2}&0\cr \frac{1}{2}&-\frac{1}{2\sqrt{3}}&\frac{1}{\sqrt{6}}\cr -\frac{1}{2}&\frac{1}{2\sqrt{3}}&\frac{2}{\sqrt{6}}\cr \frac{1}{2}&-\frac{1}{2\sqrt{3}}&\frac{1}{\sqrt{6}}} \right],\left[ \matrix{\displaystyle 2&-\frac{1}{2}&1 \cr 0&\frac{3^{3/2}}{2}&-\sqrt{3} \cr 0&0&\sqrt{6}} \right] ]\; \)

Q矩陣和R矩陣相乘得到矩陣A
(%i7) Q.R;
(%o7) \( \left[ \matrix{1&2&-1\cr 1&-1&2\cr -1&1&1\cr 1&-1&2}\right] \)

[ 本帖最後由 bugmens 於 2017-3-5 12:03 編輯 ]
作者: bugmens    時間: 2016-11-12 09:03

https://ccjou.wordpress.com/2011 ... 用/#comment-6583
如果限定\(R\)的主對角元為正值,則非奇異方陣與full col rank矩陣有唯一的\(QR\)分解,不過後者的形式為
\( A=QR=[ \matrix{Q_1&Q_2} ] \left[ \matrix{R_1 \cr 0} \right]=Q_1 R_1 \),\(Q_1\)和\(R_1\)是唯一的,但\(Q_2\)則否。


\( \left[ \matrix{0&-15&14\cr 4&32&2\cr 3&-1&4} \right]=\left[ \matrix{\displaystyle 0&-\frac{3}{5}&\frac{4}{5}\cr \frac{4}{5}&\frac{12}{25}&\frac{9}{25} \cr \frac{3}{5}&-\frac{16}{25}&-\frac{12}{25}} \right]\left[ \matrix{5&25&4\cr 0&25&-10 \cr 0&0&10} \right] \)
測試Maple和Mathematica的執行結果\(R\)矩陣對角線元素為正值。

\( \left[ \matrix{0&-15&14\cr 4&32&2\cr 3&-1&4} \right]=\left[ \matrix{\displaystyle 0&-\frac{3}{5}&-\frac{4}{5}\cr \frac{4}{5}&\frac{12}{25}&-\frac{9}{25} \cr \frac{3}{5}&-\frac{16}{25}&\frac{12}{25}} \right]\left[ \matrix{5&25&4\cr 0&25&-10 \cr 0&0&-10} \right] \)
測試MATLAB和maxima的dgeqrf指令執行結果\(R\)矩陣右下角元素為\(-10\)。

----------------------------------

\( \left[ \matrix{1&2&-1\cr 1&-1&2 \cr -1&1&1 \cr 1&-1&2} \right]=\left[ \matrix{\displaystyle \frac{1}{2}&\frac{\sqrt{3}}{2}&0 \cr \frac{1}{2}&-\frac{1}{2\sqrt{3}}&\frac{1}{\sqrt{6}}\cr -\frac{1}{2}&\frac{1}{2\sqrt{3}}&\frac{\sqrt{6}}{3} \cr \frac{1}{2}&-\frac{1}{2\sqrt{3}}&\frac{1}{\sqrt{6}} } \right] \left[ \matrix{\displaystyle 2&-\frac{1}{2}&1\cr 0&\frac{3^{3/2}}{2}&-\sqrt{3}\cr 0&0&\sqrt{6}} \right] \)
測試Maple和Mathematica會將\(R\)矩陣的0列和\(Q\)矩陣對應的列刪除。

\( \left[ \matrix{1&2&-1\cr 1&-1&2 \cr -1&1&1 \cr 1&-1&2} \right]= \left[ \matrix{\displaystyle \frac{1}{2}&\frac{\sqrt{3}}{2}&0&0\cr \frac{1}{2}&-\frac{1}{2 \sqrt{3}}&\frac{1}{\sqrt{6}}&\frac{\sqrt{6}}{2 \sqrt{3}}\cr -\frac{1}{2}&\frac{1}{2\sqrt{3}}&\frac{2}{\sqrt{6}}&0\cr \frac{1}{2}&-\frac{1}{2\sqrt{3}}&\frac{1}{\sqrt{6}}&-\frac{\sqrt{3}}{\sqrt{6}}} \right] \left[ \matrix{\displaystyle 2&-\frac{1}{2}&1\cr 0&\frac{3^{3/2}}{2}&-\sqrt{3}\cr 0&0&\sqrt{6} \cr 0&0&0} \right] \)
測試Mablab和maxima的dgeqrf指令沒有刪除\(R\)矩陣的0列和\(Q\)矩陣對應的列。

所以我修改前三篇的maxima程式碼,以符合Maple和Mathematica的執行結果。



Maple執行結果


Mathematica執行結果


MATLAB執行結果


LAPACK的dgeqrf執行結果
作者: bugmens    時間: 2017-1-31 09:25

假設矩陣\(A\)的特徵值\(\lambda_1,\lambda_2,\ldots,\lambda_n\)滿足\( |\; \lambda_1 |\;>|\; \lambda_2 |\;>\ldots>|\; \lambda_n |\;>0 \)。QR法的收斂速度在於兩個相鄰特徵值\( |\; \lambda_i |\; >|\; \lambda_{i+1} |\; \)的商\( \displaystyle \Bigg\vert\; \frac{\lambda_{i+1}}{\lambda_{i}} \Bigg\vert\; \)。

設矩陣\(A\)有6個特徵值\( \lambda_1=6,\lambda_2=5,\lambda_3=4,\lambda_4=3,\lambda_5=-1.01,\lambda_6=1 \),其中兩個特徵值\( |\; \lambda_5 |\;=|\;-1.01|\; ,|\; \lambda_6|\; =|\;1|\;\)很接近,QR法會執行很多次。
以下介紹如何給定特徵值來建構出對應的矩陣\(A\)
(1)相伴矩陣
最小多項式\( p(t)=t^n+a_{n-1}t^{n-1}+\ldots+a_1 t+a_0 \)
對應的相伴矩陣為\( C=\left[ \matrix{0&0&\ldots&0&-a_0 \cr 1&0&\ldots&0&-a_1 \cr \vdots&\vdots&\ddots&\vdots&\vdots \cr 0&\ldots&1&0&-a_{n-2} \cr0&0&\ldots&1&-a_{n-1}} \right] \)

詳細內容請參閱線代啟示錄的"多項式的相伴矩陣"
https://ccjou.wordpress.com/2011/01/17/多項式的相伴矩陣/


指定特徵值
\( \lambda_1=6,\lambda_2=5,\lambda_3=4,\lambda_4=3,\lambda_5=-1.01,\lambda_6=1 \)
特徵方程式
\( p(x)=(x-6)(x-5)(x-4)(x-3)(x+1.01)(x-1) \)
  \( =x^6-17.99x^5+117.81x^4-322.63x^3+236.39x^2+349.02x-363.6 \)
相伴矩陣
\( C=\left[ \matrix{0&0&0&0&0&363.6 \cr 1&0&0&0&0&-349.02 \cr 0&1&0&0&0&-236.39 \cr 0&0&1&0&0&322.63 \cr 0&0&0&1&0&-117.81 \cr 0&0&0&0&1&17.99} \right] \)


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

不顯示小數轉換成分數的過程
例如rat: replaced -4.25 by -17/4 = -4.25

(%i2) ratprint:false;
(%o2) false

給定特徵值產生相伴矩陣副程式
(%i3)
CompanionMatrix(eig):=block
([n:length(eig),px],
print("特徵值λ=",eig),
print("特徵方程式p(x)=",px:product((x-eig[ i ]),i,1,n)),
print("多項式展開p(x)=",px:expand(px)),
h[i,j]:=if j=n then (-ratcoef(px,x,i-1))
           else if i-j=1 then (1)
           else (0),
print("相伴矩陣C=",genmatrix(h,n,n))
)$


特徵值6,5,4,3,-1.01,1
(%i4) CompanionMatrix([6,5,4,3,-1.01,1]);
特徵值\( \lambda=\left[ 6,5,4,3,-1.01,1 \right] \)
特徵方程式\( p(x)=(x-6)(x-5)(x-4)(x-3)(x-1)(x+1.01) \)
多項式展開\( p(x)=x^6-17.99x^5+117.81x^4-322.63x^3+236.39x^2+349.02x-363.6 \)
相伴矩陣\( C=\left[ \matrix{0&0&0&0&0&363.6 \cr 1&0&0&0&0&-349.02 \cr 0&1&0&0&0&-236.39 \cr 0&0&1&0&0&322.63 \cr 0&0&0&1&0&-117.81 \cr 0&0&0&0&1&17.99} \right] \)
(%o4) \( \left[ \matrix{0&0&0&0&0&363.6 \cr 1&0&0&0&0&-349.02 \cr 0&1&0&0&0&-236.39 \cr 0&0&1&0&0&322.63 \cr 0&0&0&1&0&-117.81 \cr 0&0&0&0&1&17.99} \right] \)



(2)實對稱矩陣可正交對角化
令\(A\)為一個\(n \times n\)階實對稱矩陣且\( \lambda_1,\ldots,\lambda_n \)為\(A\)的特徵值。所謂正交對角化是指存在一個實正交矩陣(orthogonal matrix)\(Q\),\(Q^T=Q^{-1}\),使得\( Q^T AQ=\Lambda=diag(\lambda_1,\ldots,\lambda_n) \),其中\(Q\)的行向量(column vector)是\(A\)的特徵向量。
以上內容節錄線代啟示錄的"實對稱矩陣可正交對角化的證明"
https://ccjou.wordpress.com/2011 ... 對角化的證明/

利用\(\displaystyle Q=I-2 \frac{vv^T}{v^Tv}\)可以構造出正交矩陣\(Q\)
https://zh.wikipedia.org/wiki/ ... C.E5.8F.98.E6.8D.A2


給定特徵值產生對稱矩陣副程式
(%i1)
OrthDiag(eig):=block
([n:length(eig),Q,D],
print("向量v=",v:genmatrix(lambda([i,j],1),n,1)),
print("正交矩陣Q=I-2(vv^T)/(v^Tv)"),
print("Q=",ident(n),"-2",v,transpose(v),"/",transpose(v),v),
print("Q=",ident(n),"-2",v.transpose(v),"/",transpose(v).v),
print("Q=",Q:ident(n)-2*v.transpose(v)/transpose(v).v),
print("對角矩陣D=",D:genmatrix(lambda([i,j],if i=j then eig[ i ] else 0),n,n)),
print("對稱矩陣A=Q^TDQ"),
print("A=",transpose(Q),D,Q),
print("A=",transpose(Q).D.Q)
)$


特徵值6,5,4,3,-1.01,1
(%i2) OrthDiag([6,5,4,3,-101/100,1]);
向量\( v=\left[ \matrix{1 \cr 1 \cr 1 \cr 1 \cr 1 \cr 1} \right] \)
正交矩陣\(Q=I-2(vv^T)/(v^Tv)\)
\( Q=\left[ \matrix{1&0&0&0&0&0 \cr 0&1&0&0&0&0 \cr 0&0&1&0&0&0 \cr 0&0&0&1&0&0 \cr 0&0&0&0&1&0 \cr 0&0&0&0&0&1} \right]-2 \left[ \matrix{1 \cr 1 \cr 1 \cr 1 \cr 1 \cr 1} \right] \left[ \matrix{1&1&1&1&1&1} \right] / \left[ \matrix{1&1&1&1&1&1} \right] \left[ \matrix{1 \cr 1 \cr 1 \cr 1 \cr 1 \cr 1} \right] \)
\( Q=\left[ \matrix{1&0&0&0&0&0 \cr 0&1&0&0&0&0 \cr 0&0&1&0&0&0 \cr 0&0&0&1&0&0 \cr 0&0&0&0&1&0 \cr 0&0&0&0&0&1} \right]-2 \left[ \matrix{1&1&1&1&1&1 \cr 1&1&1&1&1&1 \cr 1&1&1&1&1&1 \cr 1&1&1&1&1&1 \cr 1&1&1&1&1&1 \cr 1&1&1&1&1&1} \right] /6 \)
\( Q=\left[ \matrix{\displaystyle \frac{2}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3} \cr -\frac{1}{3}&\frac{2}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3} \cr -\frac{1}{3}&-\frac{1}{3}&\frac{2}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3} \cr -\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&\frac{2}{3}&-\frac{1}{3}&-\frac{1}{3} \cr -\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&\frac{2}{3}&-\frac{1}{3} \cr -\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&\frac{2}{3}} \right] \)
對角矩陣\( \displaystyle D=\left[ \matrix{\displaystyle 6&0&0&0&0&0\cr 0&5&0&0&0&0\cr 0&0&4&0&0&0\cr 0&0&0&3&0&0\cr 0&0&0&0&-\frac{101}{100}&0\cr 0&0&0&0&0&1} \right] \)
對稱矩陣\( A=Q^TDQ \)
\( A=\left[ \matrix{\displaystyle \frac{2}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3} \cr -\frac{1}{3}&\frac{2}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3} \cr -\frac{1}{3}&-\frac{1}{3}&\frac{2}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3} \cr -\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&\frac{2}{3}&-\frac{1}{3}&-\frac{1}{3} \cr -\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&\frac{2}{3}&-\frac{1}{3} \cr -\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&\frac{2}{3}} \right]
\left[ \matrix{\displaystyle 6&0&0&0&0&0\cr 0&5&0&0&0&0\cr 0&0&4&0&0&0\cr 0&0&0&3&0&0\cr 0&0&0&0&-\frac{101}{100}&0\cr 0&0&0&0&0&1} \right]
\left[ \matrix{\displaystyle \frac{2}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3} \cr -\frac{1}{3}&\frac{2}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3} \cr -\frac{1}{3}&-\frac{1}{3}&\frac{2}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3} \cr -\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&\frac{2}{3}&-\frac{1}{3}&-\frac{1}{3} \cr -\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&\frac{2}{3}&-\frac{1}{3} \cr -\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&-\frac{1}{3}&\frac{2}{3}} \right] \)
\( A=\left[ \matrix{ \displaystyle
\frac{3599}{900}&-\frac{1501}{900}&-\frac{1201}{900}&-\frac{901}{900}&\frac{151}{450}&-\frac{301}{900}\cr
-\frac{1501}{900}&\frac{3299}{900}&-\frac{901}{900}&-\frac{601}{900}&\frac{301}{450}&-\frac{1}{900}\cr
-\frac{1201}{900}&-\frac{901}{900}&\frac{2999}{900}&-\frac{301}{900}&\frac{451}{450}&\frac{299}{900}\cr
-\frac{901}{900}&-\frac{601}{900}&-\frac{301}{900}&\frac{2699}{900}&\frac{601}{450}&\frac{599}{900}\cr
\frac{151}{450}&\frac{301}{450}&\frac{451}{450}&\frac{601}{450}&\frac{374}{450}&\frac{901}{450}\cr
-\frac{301}{900}&-\frac{1}{900}&\frac{299}{900}&\frac{599}{900}&\frac{901}{450}&\frac{2099}{900}} \right] \)

(%o2) \( \left[ \matrix{ \displaystyle
\frac{3599}{900}&-\frac{1501}{900}&-\frac{1201}{900}&-\frac{901}{900}&\frac{151}{450}&-\frac{301}{900}\cr
-\frac{1501}{900}&\frac{3299}{900}&-\frac{901}{900}&-\frac{601}{900}&\frac{301}{450}&-\frac{1}{900}\cr
-\frac{1201}{900}&-\frac{901}{900}&\frac{2999}{900}&-\frac{301}{900}&\frac{451}{450}&\frac{299}{900}\cr
-\frac{901}{900}&-\frac{601}{900}&-\frac{301}{900}&\frac{2699}{900}&\frac{601}{450}&\frac{599}{900}\cr
\frac{151}{450}&\frac{301}{450}&\frac{451}{450}&\frac{601}{450}&\frac{374}{450}&\frac{901}{450}\cr
-\frac{301}{900}&-\frac{1}{900}&\frac{299}{900}&\frac{599}{900}&\frac{901}{450}&\frac{2099}{900}} \right] \)



以QR(A,10^-5);測試執行次數
http://math.pro/db/viewthread.php?tid=2561&page=1#pid16149

(1)相伴矩陣
A:matrix([0,0,0,0,0,363.6],
         [1,0,0,0,0,-349.02],
         [0,1,0,0,0,-236.39],
         [0,0,1,0,0,322.63],
         [0,0,0,1,0,-117.81],
         [0,0,0,0,1,17.99]);
(2)實對稱矩陣
A:matrix([3599/900,-1501/900,-1201/900,-901/900,151/450,-301/900],
[-1501/900,3299/900,-901/900,-601/900,301/450,-1/900],
[-1201/900,-901/900,2999/900,-301/900,451/450,299/900],
[-901/900,-601/900,-301/900,2699/900,601/450,599/900],
[151/450,301/450,451/450,601/450,374/225,901/450],
[-301/900,-1/900,299/900,599/900,901/450,2099/900]);

原始的QR法

571次

348次



[ 本帖最後由 bugmens 於 2017-1-31 09:37 編輯 ]
作者: bugmens    時間: 2017-3-3 23:44

1.減少QR法運算-化簡成Hessenberg矩陣

在執行QR法之前,會將矩陣化簡成\( A=\left[ \matrix{*&*&*&*&* \cr *&*&*&*&* \cr 0&*&*&*&* \cr 0&0&*&*&* \cr 0&0&0&*&*} \right] \)讓左下角元為零,減少QR法運算量。

若對稱矩陣化簡後變成三對角矩陣\( A=\left[ \matrix{*&*&0&0&0 \cr *&*&*&0&0 \cr 0&*&*&*&0 \cr 0&0&*&*&* \cr 0&0&0&*&*} \right] \),減少QR法運算次數。

詳細方法請見線代啟示錄
特殊矩陣 (19):Hessenberg 矩陣
https://ccjou.wordpress.com/2013 ... hessenberg-矩陣/

Hessenberg化簡還需要Householder變換,詳細方法請見線代啟示錄
特殊矩陣 (4):Householder 矩陣
https://ccjou.wordpress.com/2009 ... householder-矩陣/




虛擬碼
Householder(x)
{
\( e1 \)為第1元為1,其餘為0的向量
定義\( sign(x_1)=+1 \)若\( x_1 \ge 0 \)
  \( sign(x_1)=-1 \)若\( x_1 <0 \)
\( \displaystyle v=\frac{x+sign(x_1)\Vert\; x \Vert\; e1}{\Vert\; x+sign(x_1)\Vert\; x \Vert\; e1 \Vert\;} \),其中\( \Vert\; * \Vert\; \)為向量長度
Householder矩陣\(H=I-2vv^T \)
\( return(H) \)
}

Hessenberg(A)
{
\(n\)為矩陣\(A\)行的個數
for i=1 ~ (n-2)
 {擷取\( A \)矩陣第\(i\)列末\(n-i\)元形成\(x\)向量
  將\(x\)向量代入\(Householder(x)\)副程式得到Householder矩陣\(H\)
  構造\( P=\left[ \matrix{I&0 \cr 0&H} \right] \),其中\(I\)為單位矩陣
  計算\(A=PAP\),得到新的矩陣A
 }
若\(A\)有整列的0向量,調整位置放在最後一列
\( return(A) \)
}



要先載入diag才能使用diag指令
(%i1) load(diag);
(%o1) C:\maxima-5.39.0\share\maxima\5.39.0_2_g5a49f11_dirty\share\contrib\diag.mac

轉換成Householder矩陣副程式
(%i2)
Householder(x):=block
([n,e1,normx,v,normv],
n:length(x),
if x=zeromatrix(n,1) then
  (print("x為0向量,H=",ident(n)),
   return(ident(n))
  ),
if x[1,1]>=0 then sign:+1 else sign:-1,/*maxima內建指令signum(0)=0會導致錯誤*/
print("x向量第1元=",x[1,1]," , sign(x_1)=",sign,",長度∥x∥=sqrt(",x,".",x,")=",normx:ratsimp(sqrt(x.x))),
print("第1元是1,其餘為0的向量e1=",e1:genmatrix(lambda([i,j],if i=1 then 1 else 0),n,1)),
print("x+sign(x_1)*∥x∥*e1=",x,"+",sign,"*",normx,"*",e1,"=",v:x+sign*normx*e1),
print("長度∥x+sign(x_1)*∥x∥*e1∥=sqrt(",v,".",v,")=",normv:ratsimp(sqrt(v.v))),
print("v=",v,"/",normv,"=",v:ratsimp(v/normv)),
print("H=I-2*v.v^T=",ident(n),"-2*",v,".",transpose(v),"=",H:ratsimp(ident(n)-2*v.transpose(v))),
return(H)
)$


轉換成Hessenberg矩陣副程式
(%i3)
Hessenberg(A):=block
([n,H,P],
n:length(A),
for i:1 thru n-2 do
  (print("擷取A=",A,"第",i,"行末",n-i,"元,x=",x:genmatrix(lambda([ii,jj],A[i+ii,i]),n-i,1)),
    print("Householder矩陣,H=",H:Householder(x)),
    print("P=",matrix([I,0],[0,"H"]),"=",P:diag([ident(i),H])),
    print("A=PAP=",P,".",A,".",P,"=",A:ratsimp(P.A.P)),
    print("------------------------")
   ),
zerorow:zeromatrix(1,n),
for i:1 thru n-1 do
  (if matrix(A[ i ])=zerorow then
     (print(A,"第",i,"列為0向量,放在最後一行",A:copymatrix(append(submatrix(i,A),zerorow)))
     )
  ),
return(A)
)$


非對稱矩陣A
(%i4)
A:matrix([1,3,2,1],
         [-2,-2,-2,1],
         [2,-2,2,-2],
         [-1,0,2,0]);


(%o4) \( \left[ \matrix{1&3&2&1 \cr -2&-2&-2&1 \cr 2&-2&2&-2 \cr-1&0&2&0} \right] \)

轉換成Hessenberg矩陣
(%5) Hessenberg(A);
擷取\( A=\left[ \matrix{1&3&2&1\cr -2&-2&-2&1\cr 2&-2&2&-2\cr -1&0&2&0} \right] \)第1行末3元,\( x=\left[ \matrix{-2 \cr 2 \cr -1} \right] \)
\(x\)向量第1元\(=-2\),\(sign(x_1)=-1\),長度\( \Vert\; x \Vert\;=sqrt( \left[ \matrix{-2 \cr 2 \cr -1} \right] \cdot \left[ \matrix{-2 \cr 2 \cr -1} \right] )=3 \)
第1元是1,其餘為0的向量\( e1=\left[ \matrix{1\cr 0 \cr 0} \right] \)
\( x+sign(x_1)*\Vert\; x \Vert\;*e1=\left[ \matrix{-2 \cr 2 \cr -1} \right]+-1*3*\left[ \matrix{1\cr 0 \cr 0} \right]=\left[ \matrix{-5 \cr 2 \cr -1} \right] \)
長度\( \Vert\; x+sign(x_1)*\Vert\; x \Vert\;*e1 \Vert\;=sqrt(\left[ \matrix{-5 \cr 2 \cr -1} \right] \cdot \left[ \matrix{-5 \cr 2 \cr -1} \right])=\sqrt{30} \)
\( v=\left[ \matrix{-5 \cr 2 \cr -1} \right]/ \sqrt{30}=\left[ \matrix{\displaystyle -\frac{5}{\sqrt{30}} \cr \frac{2}{\sqrt{30}} \cr -\frac{1}{\sqrt{30}}} \right] \)
\( H=I-2*v \cdot v^T=\left[ \matrix{1&0&0\cr 0&1&0\cr 0&0&1} \right]-2 *\left[ \matrix{\displaystyle -\frac{5}{\sqrt{30}} \cr \frac{2}{\sqrt{30}} \cr -\frac{1}{\sqrt{30}}} \right] \cdot \left[ \matrix{\displaystyle -\frac{5}{\sqrt{30}} & \frac{2}{\sqrt{30}} & -\frac{1}{\sqrt{30}}} \right]=\left[ \matrix{\displaystyle -\frac{2}{3}&\frac{2}{3}&-\frac{1}{3} \cr \frac{2}{3}&\frac{11}{15}&\frac{2}{15} \cr -\frac{1}{3}&\frac{2}{15}&\frac{14}{15}} \right] \)
Householder矩陣,\( H=\left[ \matrix{\displaystyle -\frac{2}{3}&\frac{2}{3}&-\frac{1}{3} \cr \frac{2}{3}&\frac{11}{15}&\frac{2}{15} \cr -\frac{1}{3}&\frac{2}{15}&\frac{14}{15}} \right] \)
\( P=\left[ \matrix{I&0 \cr 0&H} \right]=\left[ \matrix{\displaystyle 1&0&0&0 \cr 0&-\frac{2}{3}&\frac{2}{3}&-\frac{1}{3} \cr 0&\frac{2}{3}&\frac{11}{15}&\frac{2}{15} \cr 0&-\frac{1}{3}&\frac{2}{15}&\frac{14}{15}} \right] \)
\( A=PAP=\left[ \matrix{\displaystyle 1&0&0&0 \cr 0&-\frac{2}{3}&\frac{2}{3}&-\frac{1}{3} \cr 0&\frac{2}{3}&\frac{11}{15}&\frac{2}{15} \cr 0&-\frac{1}{3}&\frac{2}{15}&\frac{14}{15}} \right] \cdot \left[ \matrix{\displaystyle 1&3&2&1\cr -2&-2&-2&1\cr 2&-2&2&-2\cr -1&0&2&0} \right] \cdot \left[ \matrix{\displaystyle 1&0&0&0 \cr 0&-\frac{2}{3}&\frac{2}{3}&-\frac{1}{3} \cr 0&\frac{2}{3}&\frac{11}{15}&\frac{2}{15} \cr 0&-\frac{1}{3}&\frac{2}{15}&\frac{14}{15}} \right]=\left[ \matrix{\displaystyle 1&-1&\frac{18}{5}&\frac{1}{5}\cr 3&2&\frac{6}{5}&-\frac{8}{5}\cr 0&\frac{12}{5}&-\frac{42}{25}&\frac{6}{25}\cr 0&\frac{9}{5}&\frac{56}{25}&-\frac{8}{25}} \right] \)
---------------
擷取\( A=\left[ \matrix{\displaystyle 1&-1&\frac{18}{5}&\frac{1}{5}\cr 3&2&\frac{6}{5}&-\frac{8}{5}\cr 0&\frac{12}{5}&-\frac{42}{25}&\frac{6}{25}\cr 0&\frac{9}{5}&\frac{56}{25}&-\frac{8}{25}} \right] \)第2行末2元,\( x=\left[ \matrix{\displaystyle \frac{12}{5} \cr \frac{9}{5}} \right] \)
\(x\)向量第1元\(\displaystyle =\frac{12}{5}\),\( sign(x_1)=1 \),長度\( \Vert\; x \Vert\;=sqrt(\left[ \matrix{\displaystyle \frac{12}{5} \cr \frac{9}{5}} \right] \cdot \left[ \matrix{\displaystyle \frac{12}{5} \cr \frac{9}{5}} \right])=3 \)
第1元是1,其餘為0的向量\(e1=\left[ \matrix{1 \cr 0} \right] \)
\( x+sign(x_1)*\Vert\; x \Vert\;*e1=\left[ \matrix{\displaystyle \frac{12}{5} \cr \frac{9}{5}} \right]+1*3*\left[ \matrix{1 \cr 0} \right]=\left[ \matrix{\displaystyle \frac{27}{5} \cr \frac{9}{5}} \right] \)
長度\( \displaystyle \Vert\; x+sign(x_1)*\Vert\; x \Vert\;*e1 \Vert\;=sqrt(\left[ \matrix{\displaystyle \frac{27}{5} \cr \frac{9}{5}} \right] \cdot \left[ \matrix{\displaystyle \frac{27}{5} \cr \frac{9}{5}} \right])=\frac{9\sqrt{2}}{\sqrt{5}} \)
\( \displaystyle v=\left[ \matrix{\displaystyle \frac{27}{5}\cr \frac{9}{5}} \right]/ \frac{9\sqrt{2}}{\sqrt{5}}=\left[ \matrix{\displaystyle \frac{3}{\sqrt{2}\sqrt{5}} \cr \frac{1}{\sqrt{2}\sqrt{5}}} \right] \)
\( H=I-2*v \cdot v^T=\left[ \matrix{1&0 \cr 0&1} \right]-2*\left[ \matrix{\displaystyle \frac{3}{\sqrt{2}\sqrt{5}}\cr \frac{1}{\sqrt{2}\sqrt{5}}} \right] \cdot \left[ \matrix{\displaystyle \frac{3}{\sqrt{2}\sqrt{5}}& \frac{1}{\sqrt{2}\sqrt{5}}} \right]=\left[ \matrix{\displaystyle -\frac{4}{5}&-\frac{3}{5}\cr -\frac{3}{5}&\frac{4}{5}} \right] \)
Householder矩陣,\( H=\left[ \matrix{\displaystyle -\frac{4}{5}&-\frac{3}{5}\cr -\frac{3}{5}&\frac{4}{5}} \right] \)
\( P=\left[ \matrix{I&0 \cr 0&H} \right]=\left[ \matrix{\displaystyle 1&0&0&0 \cr 0&1&0&0 \cr 0&0&-\frac{4}{5}&-\frac{3}{5}\cr 0&0&-\frac{3}{5}&\frac{4}{5}} \right] \)
\( A=PAP=\left[ \matrix{\displaystyle 1&0&0&0 \cr 0&1&0&0 \cr 0&0&-\frac{4}{5}&-\frac{3}{5}\cr 0&0&-\frac{3}{5}&\frac{4}{5}} \right] \cdot \left[ \matrix{\displaystyle 1&-1&\frac{18}{5}&\frac{1}{5}\cr 3&2&\frac{6}{5}&-\frac{8}{5}\cr 0&\frac{12}{5}&-\frac{42}{25}&\frac{6}{25}\cr 0&\frac{9}{5}&\frac{56}{25}&-\frac{8}{25}} \right] \cdot \left[ \matrix{\displaystyle 1&0&0&0 \cr 0&1&0&0 \cr 0&0&-\frac{4}{5}&-\frac{3}{5}\cr 0&0&-\frac{3}{5}&\frac{4}{5}} \right]=\left[ \matrix{\displaystyle 1&-1&-3&-2\cr 3&2&0&-2\cr 0&-3&0&0\cr 0&0&-2&-2} \right] \)
---------------
(%o5) \( \left[ \matrix{1&-1&-3&-2\cr 3&2&0&-2\cr0&-3&0&0\cr 0&0&-2&-2} \right] \)

矩陣A為線代啟示錄"特殊矩陣 (4):Householder 矩陣"文章中的範例
(%i6)
A:matrix([4,1,-2,2],
              [1,2,0,1],
              [-2,0,3,-2],
              [2,1,-2,-1]);

(%o6) \( \left[ \matrix{4&1&-2&2 \cr 1&2&0&1 \cr -2&0&3&-2 \cr 2&1&-2&-1} \right] \)

轉換成Hessenberg矩陣
(%i7) Hessenberg(A);
擷取\( A=\left[ \matrix{4&1&-2&2\cr 1&2&0&1\cr -2&0&3&-2\cr 2&1&-2&-1} \right] \)第1行末3元,\( x=\left[ \matrix{1 \cr -2 \cr 2} \right] \)
\(x\)向量第1元\(=1\),\(sign(x_1)=1\),長度\( \Vert\; x \Vert\;=sqrt(\left[ \matrix{1 \cr -2 \cr 2} \right].\left[ \matrix{1 \cr -2 \cr 2} \right])=3 \)
第1元是1,其餘為0的向量\(e1=\left[ \matrix{1 \cr 0 \cr 0} \right] \)
\( x+sign(x_1)* \Vert\; x \Vert\;*e1=\left[ \matrix{1 \cr -2 \cr 2} \right]+1*3*\left[ \matrix{1 \cr 0 \cr 0} \right]=\left[ \matrix{4 \cr -2 \cr 2} \right] \)
長度\( \Vert\; x+sign(x_1)*\Vert\; x \Vert\; *e1 \Vert\;=sqrt(\left[ \matrix{4 \cr -2 \cr 2} \right].\left[ \matrix{4 \cr -2 \cr 2} \right])=2 \sqrt{6} \)
\( v=\left[ \matrix{4 \cr -2 \cr 2} \right]/2\sqrt{6}=\left[ \matrix{\displaystyle \frac{2}{\sqrt{6}} \cr -\frac{1}{\sqrt{6}} \cr \frac{1}{\sqrt{6}}} \right] \)
\( H=I-2*v.v^T=\left[ \matrix{1&0&0\cr 0&1&0\cr 0&0&1} \right]-2*\left[ \matrix{\displaystyle \frac{2}{\sqrt{6}} \cr -\frac{1}{\sqrt{6}} \cr \frac{1}{\sqrt{6}}} \right] \cdot \left[ \matrix{\displaystyle \frac{2}{\sqrt{6}}&-\frac{1}{\sqrt{6}}&\frac{1}{\sqrt{6}}} \right]=
\left[ \matrix{\displaystyle
-\frac{1}{3}&\frac{2}{3}&-\frac{2}{3} \cr
\frac{2}{3}&\frac{2}{3}&\frac{1}{3} \cr
-\frac{2}{3}&\frac{1}{3}&\frac{2}{3}} \right] \)
Householder矩陣,\( H=\left[ \matrix{\displaystyle
-\frac{1}{3}&\frac{2}{3}&-\frac{2}{3} \cr
\frac{2}{3}&\frac{2}{3}&\frac{1}{3} \cr
-\frac{2}{3}&\frac{1}{3}&\frac{2}{3}} \right] \)
\( P=\left[ \matrix{I&0 \cr 0&H} \right]=
\left[ \matrix{\displaystyle
1&0&0&0 \cr
0&-\frac{1}{3}&\frac{2}{3}&-\frac{2}{3} \cr
0&\frac{2}{3}&\frac{2}{3}&\frac{1}{3} \cr
0&-\frac{2}{3}&\frac{1}{3}&\frac{2}{3}} \right] \)
\( A=PAP=\left[ \matrix{\displaystyle
1&0&0&0 \cr
0&-\frac{1}{3}&\frac{2}{3}&-\frac{2}{3} \cr
0&\frac{2}{3}&\frac{2}{3}&\frac{1}{3} \cr
0&-\frac{2}{3}&\frac{1}{3}&\frac{2}{3}} \right] \cdot
\left[ \matrix{4&1&-2&2\cr 1&2&0&1\cr -2&0&3&-2\cr 2&1&-2&-1} \right] \cdot
\left[ \matrix{\displaystyle
1&0&0&0 \cr
0&-\frac{1}{3}&\frac{2}{3}&-\frac{2}{3} \cr
0&\frac{2}{3}&\frac{2}{3}&\frac{1}{3} \cr
0&-\frac{2}{3}&\frac{1}{3}&\frac{2}{3}} \right]=
\left[ \matrix{\displaystyle
4&-3&0&0 \cr
-3&\frac{10}{3}&1&\frac{4}{3} \cr
0&1&\frac{5}{3}&-\frac{4}{3} \cr
0&\frac{4}{3}&-\frac{4}{3}&-1} \right] \)
---------------
擷取\( A=\left[ \matrix{\displaystyle
4&-3&0&0 \cr
-3&\frac{10}{3}&1&\frac{4}{3} \cr
0&1&\frac{5}{3}&-\frac{4}{3} \cr
0&\frac{4}{3}&-\frac{4}{3}&-1} \right] \)第2行末2元,\( x=\left[ \matrix{\displaystyle 1 \cr \frac{4}{3}} \right] \)
\(x\)向量第1元=1,\(sign(x_1)=1\),長度\( \Vert\; x \Vert\;=sqrt(\left[ \matrix{\displaystyle 1 \cr \frac{4}{3}} \right] \cdot \left[ \matrix{\displaystyle 1 \cr \frac{4}{3}} \right])=\displaystyle \frac{5}{3} \)
第1元是1,其餘為0的向量\(e1=\left[ \matrix{1 \cr 0} \right] \)
\( \displaystyle x+sign(x_1)*\Vert\; x \Vert\;*e1=\left[ \matrix{\displaystyle 1 \cr \frac{4}{3}} \right]+1*\frac{5}{3}*\left[ \matrix{1 \cr 0} \right]=\left[ \matrix{\displaystyle \frac{8}{3} \cr \frac{4}{3}} \right] \)
長度\( \Vert\; x+sign(x_1)*\Vert\; x \Vert\;*e1 \Vert\;=sqrt(\left[ \matrix{\displaystyle \frac{8}{3} \cr \frac{4}{3}} \right] \cdot \left[ \matrix{\displaystyle \frac{8}{3} \cr \frac{4}{3}} \right])=\displaystyle \frac{4 \sqrt{5}}{3} \)
\( \displaystyle v=\left[ \matrix{\displaystyle \frac{8}{3} \cr \frac{4}{3}} \right] / \frac{4 \sqrt{5}}{3}=\left[ \matrix{\displaystyle \frac{2}{\sqrt{5}} \cr \frac{1}{\sqrt{5}}} \right] \)
\( H=I-2*v \cdot v^T=\left[ \matrix{1&0 \cr 0&1} \right]-2* \left[ \matrix{\displaystyle \frac{2}{\sqrt{5}} \cr \frac{1}{\sqrt{5}}} \right] \cdot \left[ \matrix{\displaystyle \frac{2}{\sqrt{5}} & \frac{1}{\sqrt{5}}} \right]=\left[ \matrix{\displaystyle -\frac{3}{5}&-\frac{4}{5}\cr -\frac{4}{5}&\frac{3}{5}} \right] \)
Householder矩陣,\( H=\left[ \matrix{\displaystyle -\frac{3}{5}&-\frac{4}{5}\cr-\frac{4}{5}&\frac{3}{5}} \right] \)
\( P=\left[ \matrix{I&0 \cr 0&H} \right]=
\left[ \matrix{\displaystyle 1&0&0&0\cr 0&1&0&0 \cr 0&0&-\frac{3}{5}&-\frac{4}{5}\cr 0&0&-\frac{4}{5}&\frac{3}{5}} \right] \)
\( A=PAP=\left[ \matrix{\displaystyle 1&0&0&0\cr 0&1&0&0 \cr 0&0&-\frac{3}{5}&-\frac{4}{5}\cr 0&0&-\frac{4}{5}&\frac{3}{5}} \right] \cdot
\left[ \matrix{\displaystyle 4&-3&0&0\cr -3&\frac{10}{3}&1&\frac{4}{3} \cr 0&1&\frac{5}{3}&-\frac{4}{3}\cr 0&\frac{4}{3}&-\frac{4}{3}&-1} \right]
\left[ \matrix{\displaystyle 1&0&0&0\cr 0&1&0&0 \cr 0&0&-\frac{3}{5}&-\frac{4}{5}\cr 0&0&-\frac{4}{5}&\frac{3}{5}} \right]=
\left[ \matrix{\displaystyle 4&-3&0&0\cr -3&\frac{10}{3}&-\frac{5}{3}&0 \cr 0&-\frac{5}{3}&-\frac{33}{25}&\frac{68}{75}\cr 0&0&\frac{68}{75}&\frac{149}{75}} \right] \)
---------------
(%o7) \( \left[ \matrix{\displaystyle 4&-3&0&0\cr -3&\frac{10}{3}&-\frac{5}{3}&0 \cr 0&-\frac{5}{3}&-\frac{33}{25}&\frac{68}{75}\cr 0&0&\frac{68}{75}&\frac{149}{75}} \right] \)
作者: bugmens    時間: 2017-3-17 21:56

Hessenberg矩陣已經具有許多零元,我們只要消去主對角下標元即可獲得QR分解的上三角矩陣\(R\)。
我將前一篇文章的Hessenberg矩陣的QR分解過程列出來,讓網友可以更加了解線代啟示錄文章內容。

QRbyGivensRotation(A)副程式取自
http://math.pro/db/viewthread.php?tid=2561&page=2#pid16180

線代啟示錄"特殊矩陣 (19):Hessenberg 矩陣"的"QR分解"
https://ccjou.wordpress.com/2013 ... hessenberg-矩陣/



利用Givens旋轉的QR分解
(%i1)
QRbyGivensRotation(A):=block
([columns:length(A[1]),rows:length(A),rank:rank(A),n,r,c,s,Qji,list,Q,R],
n:max(columns,rows),
print("1.初始值"),
Q:ident(rows),
R:copymatrix(A),
print("Q=",Q,",R=",R),
for i:1 thru min(columns,rows) do
    (for j:i+1 thru rows do
       (if R[j,i]#0 then/*若R[j,i]已經是0就不處理*/
          (print("2.以R[",i,",",i,"]=",R[i,i],"消掉R[",j,",",i,"]=",R[j,i],"為0的旋轉矩陣Q",j,i),
           print("r=sqrt(R[",i,",",i,"]^2+R[",j,",",i,"]^2)=",r:sqrt(R[i,i]^2+R[j,i]^2)),
           c:R[i,i]/r,
           s:-R[j,i]/r,
           print("cosφ=R[",i,",",i,"]/r=",R[i,i],"/",r,"=",c),
           print("sinφ=-R[",j,",",i,"]/r=-",R[j,i],"/",r,"=",s),
           Qji:ident(rows),
           Qji[i,i]:c,    Qji[i,j]:-s,
           Qji[j,i]:s,    Qji[j,j]:c,
           print("Q",j,i,"[",i,",",i,"]=cosφ=",c,",    Q",j,i,"[",i,",",j,"]=-sinφ=",-s),
           print("Q",j,i,"[",j,",",i,"]=sinφ=",s,",     Q",j,i,"[",j,",",j,"]=cosφ=",c),
           print("旋轉矩陣Q",j,i,"=",Qji),
           print("3.更新Q矩陣和R矩陣"),
           print("Q=Q.Q",j,i,"^T=",Q,".",transpose(Qji),"=",Q:ratsimp(Q.transpose(Qji))),
           print("R=Q",j,i,".R=",Qji,".",R,"=",R:ratsimp(Qji.R)),
           print("----------------------")
          )
       )
    ),
if rows#rank then
  (list:makelist(i,i,rank(A)+1,rows),
   print("4-1.刪除Q矩陣的第",list,"行"),
   print(Q,"=>",Q:apply(submatrix,cons(Q,list))),
   print("4-2.刪除R矩陣的第",list,"列"),
   print(R,"=>",R:apply(submatrix,reverse(cons(R,list))))
  ),
for i:1 thru rank do
  (if R[i,i]<0 then
    (print("R[",i,",",i,"]=",R[i,i],"<0"),
     print("R矩陣第",i,"列",R[ i ],",Q矩陣第",i,"行",col(Q,i),"乘上-1倍"),
     R[ i ]:-R[ i ],
     Q:columnop(Q,i,i,2)
    )
  ),
return(ratsimp([Q,R]))
)$


前一篇文章的Hessenberg矩陣
(%i2)
A:matrix([1,-1,-3,-2],
              [3,2,0,-2],
              [0,-3,0,0],
              [0,0,-2,-2]);

(%o2) \( \left[ \matrix{1&-1&-3&-2\cr 3&2&0&-2\cr 0&-3&0&0\cr 0&0&-2&-2} \right] \)

(%i3) [Q,R]: QRbyGivensRotation(A);
1.初始值
\( Q=\left[ \matrix{1&0&0&0\cr 0&1&0&0\cr 0&0&1&0\cr 0&0&0&1} \right],R=\left[ \matrix{1&-1&-3&-2\cr 3&2&0&-2\cr 0&-3&0&0\cr 0&0&-2&-2} \right] \)
2.以\( R[1,1]=1 \)消掉\( R[2,1]=3 \)為0的旋轉矩陣\(Q21\)
\( r=sqrt(R[1,1]^2+R[2,1]^2)=\sqrt{10} \)
\( \displaystyle cos \phi=R[1,1]/r=1/ \sqrt{10}=\frac{1}{\sqrt{10}} \)
\( \displaystyle sin \phi=-R[2,1]/r=-3/ \sqrt{10}=-\frac{3}{\sqrt{10}} \)
\( \displaystyle Q21[1,1]=cos \phi=\frac{1}{\sqrt{10}},Q21[1,2]=-sin \phi=\frac{3}{\sqrt{10}} \)
\( \displaystyle Q21[2,1]=sin \phi=-\frac{3}{\sqrt{10}},Q21[2,2]=cos \phi=\frac{1}{\sqrt{10}} \)
旋轉矩陣\( Q21=\left[ \matrix{\displaystyle \frac{1}{\sqrt{10}}&\frac{3}{\sqrt{10}}&0&0\cr -\frac{3}{\sqrt{10}}&\frac{1}{\sqrt{10}}&0&0\cr 0&0&1&0\cr 0&0&0&1} \right] \)
3.更新\(Q\)矩陣和\(R\)矩陣
\( Q=Q \cdot Q21^T=\left[ \matrix{\displaystyle 1&0&0&0\cr 0&1&0&0\cr 0&0&1&0\cr 0&0&0&1} \right] \cdot \left[ \matrix{\displaystyle \frac{1}{\sqrt{10}}&-\frac{3}{\sqrt{10}}&0&0\cr \frac{3}{\sqrt{10}}&\frac{1}{\sqrt{10}}&0&0\cr 0&0&1&0\cr 0&0&0&1} \right]=\left[ \matrix{\displaystyle \frac{1}{\sqrt{10}}&-\frac{3}{\sqrt{10}}&0&0\cr \frac{3}{\sqrt{10}}&\frac{1}{\sqrt{10}}&0&0\cr 0&0&1&0\cr 0&0&0&1} \right]\)
\( R=Q21 \cdot R=\left[ \matrix{\displaystyle \frac{1}{\sqrt{10}}&\frac{3}{\sqrt{10}}&0&0\cr -\frac{3}{\sqrt{10}}&\frac{1}{\sqrt{10}}&0&0\cr 0&0&1&0\cr 0&0&0&1} \right] \cdot \left[ \matrix{1&-1&-3&-2\cr 3&2&0&-2\cr 0&-3&0&0\cr 0&0&-2&-2} \right]=\left[ \matrix{\displaystyle \sqrt{10}&\frac{5}{\sqrt{10}}&-\frac{3}{\sqrt{10}}&-\frac{8}{\sqrt{10}}\cr 0&\frac{5}{\sqrt{10}}&\frac{9}{\sqrt{10}}&\frac{4}{\sqrt{10}}\cr 0&-3&0&0\cr 0&0&-2&-2} \right] \)
----------------------
2.以\( \displaystyle R[2,2]=\frac{5}{\sqrt{10}}\)消掉\( R[3,2]=-3 \)為0的旋轉矩陣\(Q32\)
\( \displaystyle r=sqrt(R[2,2]^2+R[3,2]^2)=\frac{\sqrt{23}}{\sqrt{2}}\)
\( \displaystyle cos \phi=R[2,2]/r=\frac{5}{\sqrt{10}}/ \frac{\sqrt{23}}{\sqrt{2}}=\frac{5\sqrt{2}}{\sqrt{10}\sqrt{23}} \)
\( \displaystyle sin \phi=-R[3,2]/r=--3/ \frac{\sqrt{23}}{\sqrt{2}}=\frac{3\sqrt{2}}{\sqrt{23}} \)
\( \displaystyle Q32[2,2]=cos \phi=\frac{5\sqrt{2}}{\sqrt{10}\sqrt{23}},Q32[2,3]=-sin \phi=-\frac{3\sqrt{2}}{\sqrt{23}} \)
\( \displaystyle Q32[3,2]=sin \phi=\frac{3\sqrt{2}}{\sqrt{23}},Q32[3,3]=cos \phi=\frac{5\sqrt{2}}{\sqrt{10}\sqrt{23}} \)
旋轉矩陣\(Q32=\left[ \matrix{\displaystyle 1&0&0&0\cr 0&\frac{5\sqrt{2}}{\sqrt{10}\sqrt{23}}&-\frac{3\sqrt{2}}{\sqrt{23}}&0\cr 0&\frac{3\sqrt{2}}{\sqrt{23}}&\frac{5\sqrt{2}}{\sqrt{10}\sqrt{23}}&0\cr 0&0&0&1} \right] \)
更新\(Q\)矩陣和\(R\)矩陣
\( Q=Q \cdot Q32^T=\left[ \matrix{\displaystyle \frac{1}{\sqrt{10}}&-\frac{3}{\sqrt{10}}&0&0\cr \frac{3}{\sqrt{10}}&\frac{1}{\sqrt{10}}&0&0\cr 0&0&1&0\cr 0&0&0&1} \right] \cdot \left[ \matrix{\displaystyle 1&0&0&0\cr 0&\frac{5\sqrt{2}}{\sqrt{10}\sqrt{23}}&\frac{3\sqrt{2}}{\sqrt{23}}&0\cr 0&-\frac{3\sqrt{2}}{\sqrt{23}}&\frac{5\sqrt{2}}{\sqrt{10}\sqrt{23}}&0\cr 0&0&0&1} \right]=\left[ \matrix{\displaystyle \frac{1}{\sqrt{10}}&-\frac{3}{\sqrt{2}\sqrt{23}}&-\frac{9\sqrt{2}}{\sqrt{10}\sqrt{23}}&0\cr \frac{3}{\sqrt{10}}&\frac{1}{\sqrt{2}\sqrt{23}}&\frac{3\sqrt{2}}{\sqrt{10}\sqrt{23}}&0\cr 0&-\frac{3\sqrt{2}}{\sqrt{23}}&\frac{5\sqrt{2}}{\sqrt{10}\sqrt{23}}&0\cr 0&0&0&1} \right] \)
\( R=Q32 \cdot R=\left[ \matrix{\displaystyle 1&0&0&0\cr 0&\frac{5\sqrt{2}}{\sqrt{10}\sqrt{23}}&-\frac{3\sqrt{2}}{\sqrt{23}}&0\cr 0&\frac{3\sqrt{2}}{\sqrt{23}}&\frac{5\sqrt{2}}{\sqrt{10}\sqrt{23}}&0\cr 0&0&0&1} \right] \cdot \left[ \matrix{\displaystyle \sqrt{10}&\frac{5}{\sqrt{10}}&-\frac{3}{\sqrt{10}}&-\frac{8}{\sqrt{10}}\cr 0&\frac{5}{\sqrt{10}}&\frac{9}{\sqrt{10}}&\frac{4}{\sqrt{10}}\cr 0&-3&0&0\cr 0&0&-2&-2} \right]=\left[ \matrix{\displaystyle \sqrt{10}&\frac{5}{\sqrt{10}}&-\frac{3}{\sqrt{10}}&-\frac{8}{\sqrt{10}}\cr 0&\frac{\sqrt{23}}{\sqrt{2}}&\frac{9}{\sqrt{2}\sqrt{23}}&\frac{2^{3/2}}{\sqrt{23}}\cr 0&0&\frac{27\sqrt{2}}{\sqrt{10}\sqrt{23}}&\frac{3 \cdot 2^{5/2}}{\sqrt{10}\sqrt{23}}\cr 0&0&-2&-2} \right] \)
----------------------
2.以\( \displaystyle R[3,3]=\frac{27 \sqrt{2}}{\sqrt{10}\sqrt{23}} \)消掉\(R[4,3]=-2\)為0的旋轉矩陣\(Q43\)
\( \displaystyle r=sqrt(R[3,3]^2+R[4,3]^2)=\frac{\sqrt{1189}}{\sqrt{115}} \)
\( \displaystyle cos \phi=R[3,3]/r=\frac{27\sqrt{2}}{\sqrt{10}\sqrt{23}}/ \frac{\sqrt{1189}}{\sqrt{115}}=\frac{27 \sqrt{2}\sqrt{115}}{\sqrt{10}\sqrt{23}\sqrt{1189}} \)
\( \displaystyle sin \phi=-R[4,3]/r=--2/ \frac{\sqrt{1189}}{\sqrt{115}}=\frac{2 \sqrt{115}}{\sqrt{1189}} \)
\( \displaystyle Q43[3,3]=cos \phi=\frac{27\sqrt{2}\sqrt{115}}{\sqrt{10}\sqrt{23}\sqrt{1189}},Q43[3,4]=-sin \phi=-\frac{2\sqrt{115}}{\sqrt{1189}} \)
\( \displaystyle Q43[4,3]=sin \phi=\frac{2\sqrt{115}}{\sqrt{1189}},Q43[4,4]=cos \phi=\frac{27\sqrt{2}\sqrt{115}}{\sqrt{10}\sqrt{23}\sqrt{1189}} \)
旋轉矩陣\( Q43=\left[ \matrix{\displaystyle 1&0&0&0\cr 0&1&0&0\cr 0&0&\frac{27\sqrt{2}\sqrt{115}}{\sqrt{10}\sqrt{23}\sqrt{1189}}&-\frac{2\sqrt{115}}{\sqrt{1189}}\cr 0&0&\frac{2\sqrt{115}}{\sqrt{1189}}&\frac{27\sqrt{2}\sqrt{115}}{\sqrt{10}\sqrt{23}\sqrt{1189}}} \right] \)
3.更新\(Q\)矩陣和\(R\)矩陣
\( Q=Q.Q43^T=\left[ \matrix{\displaystyle \frac{1}{\sqrt{10}}&-\frac{3}{\sqrt{2}\sqrt{23}}&-\frac{9\sqrt{2}}{\sqrt{10}\sqrt{23}}&0\cr \frac{3}{\sqrt{10}}&\frac{1}{\sqrt{2}\sqrt{23}}&\frac{3\sqrt{2}}{\sqrt{10}\sqrt{23}}&0\cr 0&-\frac{3\sqrt{2}}{\sqrt{23}}&\frac{5\sqrt{2}}{\sqrt{10}\sqrt{23}}&0\cr 0&0&0&1} \right] \cdot \left[ \matrix{\displaystyle 1&0&0&0\cr 0&1&0&0\cr 0&0&\frac{27\sqrt{2}\sqrt{115}}{\sqrt{10}\sqrt{23}\sqrt{1189}}&\frac{2\sqrt{115}}{\sqrt{1189}}\cr 0&0&-\frac{2\sqrt{115}}{\sqrt{1189}}&\frac{27\sqrt{2}\sqrt{115}}{\sqrt{10}\sqrt{23}\sqrt{1189}}} \right]=\left[ \matrix{\displaystyle \frac{1}{\sqrt{10}}&-\frac{3}{\sqrt{2}\sqrt{23}}&-\frac{243}{\sqrt{115}\sqrt{1189}}&-\frac{9 \cdot 2^{3/2}\sqrt{115}}{\sqrt{10}\sqrt{23}\sqrt{1189}}\cr \frac{3}{\sqrt{10}}&\frac{1}{\sqrt{2}\sqrt{23}}&\frac{81}{\sqrt{115}\sqrt{1189}}&\frac{3 \cdot 2^{3/2}\sqrt{115}}{\sqrt{10}\sqrt{23}\sqrt{1189}}\cr 0&-\frac{3\sqrt{2}}{\sqrt{23}}&\frac{27\sqrt{115}}{23\sqrt{1189}}&\frac{5 \cdot 2^{3/2}\sqrt{115}}{\sqrt{10}\sqrt{23}\sqrt{1189}}\cr 0&0&-\frac{2\sqrt{115}}{\sqrt{1189}}&\frac{27\sqrt{2}\sqrt{115}}{\sqrt{10}\sqrt{23}\sqrt{1189}}} \right] \)
\( R=Q43 \cdot R=\left[ \matrix{\displaystyle 1&0&0&0\cr 0&1&0&0\cr 0&0&\frac{27\sqrt{2}\sqrt{115}}{\sqrt{10}\sqrt{23}\sqrt{1189}}&-\frac{2\sqrt{115}}{\sqrt{1189}}\cr 0&0&\frac{2\sqrt{115}}{\sqrt{1189}}&\frac{27\sqrt{2}\sqrt{115}}{\sqrt{10}\sqrt{23}\sqrt{1189}}} \right] \cdot \left[ \matrix{\displaystyle \sqrt{10}&\frac{5}{\sqrt{10}}&-\frac{3}{\sqrt{10}}&-\frac{8}{\sqrt{10}}\cr 0&\frac{\sqrt{23}}{\sqrt{2}}&\frac{9}{\sqrt{2}\sqrt{23}}&\frac{2^{3/2}}{\sqrt{23}}\cr 0&0&\frac{27\sqrt{2}}{\sqrt{10}\sqrt{23}}&\frac{3 \cdot 2^{5/2}}{\sqrt{10}\sqrt{23}}\cr 0&0&-2&-2} \right]=\left[ \matrix{\displaystyle \sqrt{10}&\frac{5}{\sqrt{10}}&-\frac{3}{\sqrt{10}}&-\frac{8}{\sqrt{10}}\cr 0&\frac{\sqrt{23}}{\sqrt{2}}&\frac{9}{\sqrt{2}\sqrt{23}}&\frac{2^{3/2}}{\sqrt{23}}\cr 0&0&\frac{\sqrt{1189}}{\sqrt{115}}&\frac{784}{\sqrt{115}\sqrt{1189}}\cr 0&0&0&-\frac{3\sqrt{2}\sqrt{10}\sqrt{115}}{\sqrt{23}\sqrt{1189}}} \right] \)
----------------------
\( \displaystyle R[4,4]=-\frac{3\sqrt{2}\sqrt{10}\sqrt{115}}{\sqrt{23}\sqrt{1189}}<0 \)
\( R \)矩陣第4列\( \displaystyle [0,0,0,-\frac{3\sqrt{2}\sqrt{10}\sqrt{115}}{\sqrt{23}\sqrt{1189}}],Q \)矩陣第4行\( \left[ \matrix{\displaystyle -\frac{9 \cdot 2^{3/2}\sqrt{115}}{\sqrt{10}\sqrt{23}\sqrt{1189}} \cr\frac{3 \cdot 2^{3/2}\sqrt{115}}{\sqrt{10}\sqrt{23}\sqrt{1189}} \cr\frac{5 \cdot 2^{3/2}\sqrt{115}}{\sqrt{10}\sqrt{23}\sqrt{1189}} \cr\frac{27\sqrt{2}\sqrt{115}}{\sqrt{10}\sqrt{23}\sqrt{1189}}}\right] \)乘上\(-1\)倍
0 errors, 0 warnings
(%o3) \( [\left[ \matrix{\displaystyle \frac{1}{\sqrt{10}}&-\frac{3}{\sqrt{2}\sqrt{23}}&-\frac{243}{\sqrt{115}\sqrt{1189}}&\frac{9\sqrt{2}\sqrt{10}\sqrt{115}}{5\sqrt{23}\sqrt{1189}}\cr \frac{3}{\sqrt{10}}&\frac{1}{\sqrt{2}\sqrt{23}}&\frac{81}{\sqrt{115}\sqrt{1189}}&-\frac{3\sqrt{2}\sqrt{10}\sqrt{115}}{5\sqrt{23}\sqrt{1189}}\cr 0&-\frac{3\sqrt{2}}{\sqrt{23}}&\frac{27\sqrt{115}}{23\sqrt{1189}}&-\frac{\sqrt{2}\sqrt{10}\sqrt{115}}{\sqrt{23}\sqrt{1189}}\cr 0&0&-\frac{2\sqrt{115}}{\sqrt{1189}}&-\frac{27\sqrt{10}\sqrt{115}}{5\sqrt{2}\sqrt{23}\sqrt{1189}}} \right] ,\left[ \matrix{\displaystyle \sqrt{10}&\frac{5}{\sqrt{10}}&-\frac{3}{\sqrt{10}}&-\frac{8}{\sqrt{10}}\cr 0&\frac{\sqrt{23}}{\sqrt{2}}&\frac{9}{\sqrt{2}\sqrt{23}}&\frac{2^{3/2}}{\sqrt{23}}\cr 0&0&\frac{\sqrt{1189}}{\sqrt{115}}&\frac{784}{\sqrt{115}\sqrt{1189}}\cr 0&0&0&\frac{3\sqrt{2}\sqrt{10}\sqrt{115}}{\sqrt{23}\sqrt{1189}}} \right]] \)

\(R.Q\)相乘後仍是Hessenberg矩陣
(%i4) ratsimp(R.Q);
(%o4) \( \left[ \matrix{\displaystyle
\frac{5}{2}&-\frac{7\sqrt{10}}{5 \cdot 2^{3/2}\sqrt{23}}&-\frac{59\sqrt{10}}{\sqrt{115}\sqrt{1189}}&\frac{99 \cdot 2^{3/2}\sqrt{115}}{5\sqrt{23}\sqrt{1189}}\cr
\frac{3\sqrt{23}}{\sqrt{2}\sqrt{10}}&-\frac{31}{46}&\frac{1079\sqrt{2}}{\sqrt{23}\sqrt{115}\sqrt{1189}}&-\frac{168\sqrt{10}}{\sqrt{115}\sqrt{1189}}\cr
0&-\frac{3\sqrt{2}\sqrt{1189}}{\sqrt{23}\sqrt{115}}&-\frac{3961}{27347}&-\frac{16529\sqrt{2}\sqrt{10}}{5945\sqrt{23}}\cr
0&0&-\frac{15 \cdot 2^{3/2}\sqrt{10}\sqrt{23}}{1189}&-\frac{810}{1189}} \right] \)

前一篇文章的三對角矩陣
(%i5)
A:matrix([4,-3,0,0],
              [-3,10/3,-5/3,0],
              [0,-5/3,-33/25,68/75],
              [0,0,68/75,149/75]);

(%o5) \( \left[ \matrix{\displaystyle 4&-3&0&0\cr -3&\frac{10}{3}&-\frac{5}{3}&0\cr 0&-\frac{5}{3}&-\frac{33}{25}&\frac{68}{75}\cr 0&0&\frac{68}{75}&\frac{149}{75}} \right] \)

(%i6) [Q,R]: QRbyGivensRotation(A);
1.初始值
\( Q=\left[ \matrix{1&0&0&0\cr 0&1&0&0\cr 0&0&1&0\cr 0&0&0&1} \right],R=\left[ \matrix{\displaystyle 4&-3&0&0\cr -3&\frac{10}{3}&-\frac{5}{3}&0\cr 0&-\frac{5}{3}&-\frac{33}{25}&\frac{68}{75}\cr 0&0&\frac{68}{75}&\frac{149}{75}} \right] \)
2.以\(R[1,1]=4\)消掉\(R[2,1]=-3\)為0的旋轉矩陣\(Q21\)
\( r=sqrt(R[1,1]^2+R[2,1]^2)=5 \)
\( \displaystyle cos \phi=R[1,1]/r=4/5=\frac{4}{5} \)
\( \displaystyle sin \phi=-R[2,1]/r=--3/5=\frac{3}{5} \)
\( \displaystyle Q21[1,1]=cos \phi=\frac{4}{5},Q21[1,2]=-sin \phi=-\frac{3}{5} \)
\( \displaystyle Q21[2,1]=sin \phi=\frac{3}{5},Q21[2,2]=cos \phi=\frac{4}{5} \)
旋轉矩陣\(Q21=\left[ \matrix{\displaystyle \frac{4}{5}&-\frac{3}{5}&0&0\cr \frac{3}{5}&\frac{4}{5}&0&0\cr 0&0&1&0\cr 0&0&0&1} \right] \)
3.更新\(Q\)矩陣和\(R\)矩陣
\( Q=Q \cdot Q21^T=\left[ \matrix{\displaystyle 1&0&0&0\cr 0&1&0&0\cr 0&0&1&0\cr 0&0&0&1} \right] \cdot \left[ \matrix{\displaystyle \frac{4}{5}&\frac{3}{5}&0&0\cr -\frac{3}{5}&\frac{4}{5}&0&0\cr 0&0&1&0\cr 0&0&0&1} \right]=\left[ \matrix{\displaystyle \frac{4}{5}&\frac{3}{5}&0&0\cr -\frac{3}{5}&\frac{4}{5}&0&0\cr 0&0&1&0\cr 0&0&0&1} \right] \)
\( R=Q21 \cdot R=\left[ \matrix{\displaystyle \frac{4}{5}&-\frac{3}{5}&0&0\cr \frac{3}{5}&\frac{4}{5}&0&0\cr 0&0&1&0\cr 0&0&0&1} \right] \cdot \left[ \matrix{\displaystyle 4&-3&0&0\cr -3&\frac{10}{3}&-\frac{5}{3}&0\cr 0&-\frac{5}{3}&-\frac{33}{25}&\frac{68}{75}\cr 0&0&\frac{68}{75}&\frac{149}{75}} \right]=\left[ \matrix{\displaystyle 5&-\frac{22}{5}&1&0\cr 0&\frac{13}{15}&-\frac{4}{3}&0\cr 0&-\frac{5}{3}&-\frac{33}{25}&\frac{68}{75}\cr 0&0&\frac{68}{75}&\frac{149}{75}} \right] \)
----------------------
2.以\( \displaystyle R[2,2]=\frac{13}{15} \)消掉\( \displaystyle R[3,2]=-\frac{5}{3} \)為0的旋轉矩陣\( Q32 \)
\( \displaystyle r=sqrt(R[2,2]^2+R[3,2]^2)=\frac{\sqrt{794}}{15} \)
\( \displaystyle cos \phi=R[2,2]/r=\frac{13}{15}/ \frac{\sqrt{794}}{15}=\frac{13}{\sqrt{794}} \)
\( \displaystyle sin \phi=-R[3,2]/r=--\frac{5}{3}/ \frac{\sqrt{794}}{15}=\frac{25}{\sqrt{794}} \)
\( \displaystyle Q32[2,2]=cos \phi=\frac{13}{\sqrt{794}},Q32[2,3]=-sin \phi=-\frac{25}{\sqrt{794}} \)
\( \displaystyle Q32[3,2]=sin \phi=\frac{25}{\sqrt{794}},Q32[3,3]=cos \phi=\frac{13}{\sqrt{794}} \)
旋轉矩陣\( Q32=\left[ \matrix{\displaystyle 1&0&0&0\cr 0&\frac{13}{\sqrt{794}}&-\frac{25}{\sqrt{794}}&0\cr 0&\frac{25}{\sqrt{794}}&\frac{13}{\sqrt{794}}&0\cr 0&0&0&1} \right] \)
3.更新\(Q\)矩陣和\(R\)矩陣
\( Q=Q \cdot Q32^T=\left[ \matrix{\displaystyle \frac{4}{5}&\frac{3}{5}&0&0\cr -\frac{3}{5}&\frac{4}{5}&0&0\cr 0&0&1&0\cr 0&0&0&1} \right] \cdot \left[ \matrix{\displaystyle 1&0&0&0\cr 0&\frac{13}{\sqrt{794}}&\frac{25}{\sqrt{794}}&0\cr 0&-\frac{25}{\sqrt{794}}&\frac{13}{\sqrt{794}}&0\cr 0&0&0&1} \right]=\left[ \matrix{\displaystyle \frac{4}{5}&\frac{39}{5\sqrt{794}}&\frac{15}{\sqrt{794}}&0\cr -\frac{3}{5}&\frac{52}{5\sqrt{794}}&\frac{20}{\sqrt{794}}&0\cr 0&-\frac{25}{\sqrt{794}}&\frac{13}{\sqrt{794}}&0\cr 0&0&0&1} \right] \)
\( R=Q32 \cdot R=\left[ \matrix{\displaystyle 1&0&0&0\cr 0&\frac{13}{\sqrt{794}}&-\frac{25}{\sqrt{794}}&0\cr 0&\frac{25}{\sqrt{794}}&\frac{13}{\sqrt{794}}&0\cr 0&0&0&1} \right] \cdot \left[ \matrix{\displaystyle 5&-\frac{22}{5}&1&0\cr 0&\frac{13}{15}&-\frac{4}{3}&0\cr 0&-\frac{5}{3}&-\frac{33}{25}&\frac{68}{75}\cr 0&0&\frac{68}{75}&\frac{149}{75}} \right]=\left[ \matrix{\displaystyle 5&-\frac{22}{5}&1&0\cr 0&\frac{\sqrt{794}}{15}&\frac{47}{3\sqrt{794}}&-\frac{68}{3\sqrt{794}}\cr 0&0&-\frac{3787}{75\sqrt{794}}&\frac{884}{75\sqrt{794}}\cr 0&0&\frac{68}{75}&\frac{149}{75}} \right] \)
----------------------
2.以\( \displaystyle R[3,3]=-\frac{3787}{75\sqrt{794}} \)消掉\( \displaystyle R[4,3]=\frac{68}{75} \)為0的旋轉矩陣\(Q43\)
\( \displaystyle r=sqrt(R[3,3]^2+R[4,3]^2)=\frac{\sqrt{80057}}{5\sqrt{794}} \)
\( \displaystyle cos \phi=R[3,3]/r=-\frac{3787}{75\sqrt{794}}/ \frac{\sqrt{80057}}{5\sqrt{794}}=-\frac{3787}{15\sqrt{80057}} \)
\( \displaystyle sin \phi=-R[4,3]/r=-\frac{68}{75}/ \frac{\sqrt{80057}}{5\sqrt{794}}=-\frac{68\sqrt{794}}{15\sqrt{80057}} \)
\( \displaystyle Q43[3,3]=cos \phi=-\frac{3787}{15\sqrt{80057}},Q43[3,4]=-sin \phi=\frac{68\sqrt{794}}{15\sqrt{80057}} \)
\( \displaystyle Q43[4,3]=sin \phi=-\frac{68\sqrt{794}}{15\sqrt{80057}},Q43[4,4]=cos \phi=-\frac{3787}{15\sqrt{80057}} \)
旋轉矩陣\(Q43=\left[ \matrix{\displaystyle 1&0&0&0\cr 0&1&0&0\cr 0&0&-\frac{3787}{15\sqrt{80057}}&\frac{68\sqrt{794}}{15\sqrt{80057}}\cr 0&0&-\frac{68\sqrt{794}}{15\sqrt{80057}}&-\frac{3787}{15\sqrt{80057}}} \right] \)
3.更新\(Q\)矩陣和\(R\)矩陣
\( Q=Q \cdot Q43^T=\left[ \matrix{\displaystyle \frac{4}{5}&\frac{39}{5\sqrt{794}}&\frac{15}{\sqrt{794}}&0\cr -\frac{3}{5}&\frac{52}{5\sqrt{794}}&\frac{20}{\sqrt{794}}&0\cr 0&-\frac{25}{\sqrt{794}}&\frac{13}{\sqrt{794}}&0\cr 0&0&0&1} \right] \cdot \left[ \matrix{\displaystyle 1&0&0&0\cr 0&1&0&0\cr 0&0&-\frac{3787}{15\sqrt{80057}}&-\frac{68\sqrt{794}}{15\sqrt{80057}}\cr 0&0&\frac{68\sqrt{794}}{15\sqrt{80057}}&-\frac{3787}{15\sqrt{80057}}} \right]=\left[ \matrix{\displaystyle \frac{4}{5}&\frac{39}{5\sqrt{794}}&-\frac{3787}{\sqrt{794}\sqrt{80057}}&-\frac{68}{\sqrt{80057}}\cr -\frac{3}{5}&\frac{52}{5\sqrt{794}}&-\frac{15148}{3\sqrt{794}\sqrt{80057}}&-\frac{272}{3\sqrt{80057}}\cr 0&-\frac{25}{\sqrt{794}}&-\frac{49231}{15\sqrt{794}\sqrt{80057}}&-\frac{884}{15\sqrt{80057}}\cr 0&0&\frac{68\sqrt{794}}{15\sqrt{80057}}&-\frac{3787}{15\sqrt{80057}}} \right] \)
\( R=Q43 \cdot R=\left[ \matrix{\displaystyle 1&0&0&0\cr 0&1&0&0\cr 0&0&-\frac{3787}{15\sqrt{80057}}&\frac{68\sqrt{794}}{15\sqrt{80057}}\cr 0&0&-\frac{68\sqrt{794}}{15\sqrt{80057}}&-\frac{3787}{15\sqrt{80057}}} \right] \cdot \left[ \matrix{\displaystyle 5&-\frac{22}{5}&1&0\cr 0&\frac{\sqrt{794}}{15}&\frac{47}{3\sqrt{794}}&-\frac{68}{3\sqrt{794}}\cr 0&0&-\frac{3787}{75\sqrt{794}}&\frac{884}{75\sqrt{794}}\cr 0&0&\frac{68}{75}&\frac{149}{75}} \right]=\left[ \matrix{\displaystyle 5&-\frac{22}{5}&1&0\cr 0&\frac{\sqrt{794}}{15}&\frac{47}{3\sqrt{794}}&-\frac{68}{3\sqrt{794}}\cr 0&0&\frac{\sqrt{80057}}{5\sqrt{794}}&\frac{20876}{5\sqrt{794}\sqrt{80057}}\cr 0&0&0&-\frac{555}{\sqrt{80057}}} \right] \)
----------------------
\( \displaystyle R[4,4]=-\frac{555}{\sqrt{80057}}<0 \)
\( R \)矩陣第4列\( \displaystyle [0,0,0,-\frac{555}{\sqrt{80057}}],Q \)矩陣第4行\( \left[ \matrix{\displaystyle -\frac{68}{\sqrt{80057}}\cr -\frac{272}{3\sqrt{80057}}\cr -\frac{884}{15\sqrt{80057}}\cr -\frac{3787}{15\sqrt{80057}}} \right] \)乘上\(-1\)倍
(%o6) \( [\left[ \matrix{\displaystyle \frac{4}{5}&\frac{39}{5\sqrt{794}}&-\frac{3787}{\sqrt{794}\sqrt{80057}}&\frac{68}{\sqrt{80057}}\cr -\frac{3}{5}&\frac{52}{5\sqrt{794}}&-\frac{15148}{3\sqrt{794}\sqrt{80057}}&\frac{272}{3\sqrt{80057}}\cr 0&-\frac{25}{\sqrt{794}}&-\frac{49231}{15\sqrt{794}\sqrt{80057}}&\frac{884}{15\sqrt{80057}}\cr 0&0&\frac{68\sqrt{794}}{15\sqrt{80057}}&\frac{3787}{15\sqrt{80057}}} \right],\left[ \matrix{\displaystyle 5&-\frac{22}{5}&1&0\cr 0&\frac{\sqrt{794}}{15}&\frac{47}{3\sqrt{794}}&-\frac{68}{3\sqrt{794}}\cr 0&0&\frac{\sqrt{80057}}{5\sqrt{794}}&\frac{20876}{5\sqrt{794}\sqrt{80057}}\cr 0&0&0&\frac{555}{\sqrt{80057}}} \right]] \)

\(R.Q\)相乘後仍是三對角矩陣
(%i7) ratsimp(R.Q);
(%o7) \( \left[ \matrix{\displaystyle
\frac{166}{25}&-\frac{\sqrt{794}}{25}&0&0\cr
-\frac{\sqrt{794}}{25}&\frac{3971}{19850}&-\frac{5\sqrt{80057}}{794}&0\cr
0&-\frac{5\sqrt{80057}}{794}&-\frac{37521989}{63565258}&\frac{2516\sqrt{794}}{80057}\cr
0&0&\frac{2516\sqrt{794}}{80057}&\frac{140119}{80057}} \right] \)

[ 本帖最後由 bugmens 於 2017-3-21 19:58 編輯 ]




歡迎光臨 Math Pro 數學補給站 (https://math.pro/db/) 論壇程式使用 Discuz! 6.1.0