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)　−4−11032115

執行QR分解
(%i3)　[Q,R]: QRbyGramSchmidt(A);
1.正交化
y1=x1=−4−11
y1=sqrt(−4−11−4−11)=32
------------------------
y2=x2=032
y2=y2−(y1Tx2)(y12)y1
=032−(−4−11032)(−4−11−4−11)−4−11
=032−−118−4−11
=−9218531837
y2=sqrt(−9218531837−9218531837)=32233
------------------------
y3=x3=115
y3=y3−(y1Tx3)(y12)y1
=115−(−4−11115)(−4−11−4−11)−4−11
=115−018−4−11
=115
y3=y3−(y2Tx3)(y22)y2
=115−(−9218531837115)(−9218531837−9218531837)−9218531837
=115−1318233−9218531837
=233285−233456233684
y3=sqrt(233285−233456233684233285−233456233684)=57233
------------------------
2.單位化得正交矩陣Q
q1=y1y1=−4−1132=−3232−132132
q2=y2y2=−921853183732233=−232323353322333732233
q3=y3y3=233285−23345623368457233=5233−823312233
Q=−3232−132132−2323233533223337322335233−823312233
------------------------
3.上三角矩陣R
R11=y1=32 
R12=q1Tx2=−3232−132132032=−132 
R13=q1Tx3=−3232−132132115=0 
R22=y2=32233 
R23=q2Tx3=−232323353322333732233115=233392 
R33=y3=57233 
R=3200−132322330023339257233
(%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 編輯 ]

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 編輯 ]

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 編輯 ]

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執行結果

假設矩陣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法	1142次	348次

[ 本帖最後由 bugmens 於 2017-5-28 11:57 編輯 ]

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]

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 編輯 ]

拿之前的例子https://math.pro/db/viewthread.php?tid=2561&page=2#pid16487
相伴矩陣 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] 也是Hessenberg矩陣，QR法執行次數和原來一樣。

實對稱矩陣 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] 化簡成三對角矩陣，QR法執行次數反而增加到444次。


Householder,Hessenberg副程式取自https://math.pro/db/viewthread.php?tid=2561&page=2#pid16538
只需要執行結果，不需要計算過程，所以去掉print指令。

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

轉換成Householder矩陣副程式(去掉print指令)
(%i2)
Householder(x):=block
([n,e1,normx,v,normv],
n:length(x),
if x=zeromatrix(n,1) then
  (return(ident(n))
  ),
if x[1,1]>=0 then sign:+1 else sign:-1,
normx:ratsimp(sqrt(x.x)),
e1:genmatrix(lambda([i,j],if i=1 then 1 else 0),n,1),
v:x+sign*normx*e1,
normv:ratsimp(sqrt(v.v)),
v:ratsimp(v/normv),
H:ratsimp(ident(n)-2*v.transpose(v)),
return(H)
)$

轉換成Hessenberg矩陣副程式(去掉print指令)
(%i3)
Hessenberg(A):=block
([n,H,P],
n:length(A),
for i:1 thru n-2 do
  (x:genmatrix(lambda([ii,jj],A[i+ii,i]),n-i,1),
   H:Householder(x),
   P:diag([ident(i),H]),
   A:ratsimp(P.A.P)
  ),
zerorow:zeromatrix(1,n),
for i:1 thru n-1 do
  (if matrix(A[ i ])=zerorow then
     (A:copymatrix(append(submatrix(i,A),zerorow))
     )
  ),
return(A)
)$

實對稱矩陣
(%i4)
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]);
(%o4)　 \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]

設定為數值運算
(%i5)　numer:true;
(%o5)　true

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

設定小數點底下第6位四捨五入
(%i7)　fpprintprec:6;
(%o7)　7

轉換成Hessenberg矩陣
(%i8)　Hessenberg(A);
(%o8)　 \left[ \matrix{3.99889&2.40601&1.09577\cdot 10^{-15}&1.52435\cdot 10^{-16}&        3.36937\cdot 10^{-15}&-1.96632\cdot 10^{-15} \cr 2.40601&1.68454&-1.93868&-1.46452\cdot 10^{-15}&4.19666\cdot 10^{-15}&-9.19618\cdot 10^{-15} \cr 1.09577\cdot 10^{-15}&-1.93868&2.37077&-1.67177&2.83246\cdot 10^{-16}&5.38814\cdot 10^{-16} \cr 1.52435\cdot 10^{-16}&-1.64384\cdot 10^{-15}&-1.67177&2.35489&1.31944&-6.66134\cdot 10^{-16} \cr 3.36937\cdot 10^{-15}&3.6443\cdot 10^{-15}&-8.64785\cdot 10^{-16}&1.31944&3.60161&-0.725455 \cr -1.96632\cdot 10^{-15}&-8.90842\cdot 10^{-15}&2.04396\cdot 10^{-16}&        -6.66134\cdot 10^{-16}&-0.725455&        3.9793} \right]


接近0的數字當成0，三對角矩陣為 \left[ \matrix{3.99889&2.40601&0&0&0&0 \cr 2.40601&1.68454&-1.93868&0&0&0 \cr 0&-1.93868&2.37077&-1.67177&0&0 \cr 0&0&-1.67177&2.35489&1.31944&0 \cr 0&0&0&1.31944&3.60161&-0.725455 \cr 0&0&0&0&-0.725455&        3.9793} \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法	1142次	348次
	(3)Hessenberg矩陣 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]);	(4)三對角矩陣 A:matrix([3.99889,2.40601,0,0,0,0],    [2.40601,1.68454,-1.93868,0,0,0],    [0,-1.93868,2.37077,-1.67177,0,0],    [0,0,-1.67177,2.35489,1.31944,0],    [0,0,0,1.31944,3.60161,-0.72545],    [0,0,0,0,-0.72545,3.97930]);
原始的QR法	1142次	444次

雖然化簡成Hessenberg矩陣可以減少QR分解的運算量，但QR法執行次數卻增加了。

[ 本帖最後由 bugmens 於 2017-5-28 14:46 編輯 ]

2.減少QR法運算-Rayleigh Quotient Shift

之前Power迭代法有介紹利用 A-\alpha I 來加快收斂速度。
https://math.pro/db/viewthread.php?tid=2561&page=1#pid16043
在QR迭代法則取 \alpha=A[n,n] 也就是矩陣的右下角元當作位移值。
－－－－－－－－－－－－－－－－－－－－－－－－－－－
虛擬碼
n為A矩陣的列數或行數
1.對 k=0,1,2,\ldots ，直到 \lambda(k) 收斂為止，計算
2. A_k-A[n,n]\cdot I=Q_kR_k ，計算 A_k-A[n,n]\cdot I 的QR分解
3. A_{k+1}=R_k Q_k+A[n,n]\cdot I ，將 Q_k 和R_k對調後相乘再加 A[n,n]\cdot I 得到下一個矩陣A_{k+1}
4.特徵值 \lambda(1),\lambda(2),\ldots,\lambda(n) 在矩陣 A_{k+1} 的對角線上。
－－－－－－－－－－－－－－－－－－－－－－－－－－－
計算 A=\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] 全部的特徵值。
[計算過程]
第一輪：
計算 A-A[n,n]\cdot I的QR分解
\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&\bbox[border:1px solid blue]{17.99}} \right]-17.99I=\left[ \matrix{-17.99&0&0&0&0&363.6\cr 1&-17.99&0&0&0&-349.02\cr 0&1&-17.99&0&0&-236.39\cr 0&0&1&-17.99&0&322.63\cr 0&0&0&1&-17.99&-117.81\cr 0&0&0&0&1&0} \right]=QR
將Q和R對調後相乘再加 A[n,n]\cdot I 得到下一個矩陣A
RQ+17.99 I= \left[ \matrix{-0.0554152&-0.00307557&-1.70959\cdot 10^{-4}&-9.50302\cdot 10^{-6}&-21.2568&-381.819\cr 0.998463&-1.70696\cdot 10^{-4}&-9.48832\cdot 10^{-6}&-5.27422\cdot 10^{-7}&17.4905&314.168\cr 0&0.999995&-5.2742\cdot 10^{-7}&-2.93174\cdot 10^{-8}&15.1092&271.394\cr 0&0&1.0&-1.62964\cdot 10^{-9}&-17.458&-313.584\cr 0&0&0&1.0&5.6338&101.195\cr 0&0&0&0&-0.310553&12.4118} \right]
特徵值 \left[ -0.0554152,-1.70696\cdot 10^{-4},-5.2742\cdot 10^{-7},-1.62964\cdot 10^{-9},5.6338,12.4118 \right]

第二輪：
計算A-A[n,n]\cdot I的QR分解
\left[ \matrix{-0.0554152&-0.00307557&-1.70959\cdot 10^{-4}&-9.50302\cdot 10^{-6}&-21.2568&-381.819\cr 0.998463&-1.70696\cdot 10^{-4}&-9.48832\cdot 10^{-6}&-5.27422\cdot 10^{-7}&17.4905&314.168\cr 0&0.999995&-5.2742\cdot 10^{-7}&-2.93174\cdot 10^{-8}&15.1092&271.394\cr 0&0&1&-1.62964\cdot 10^{-9}&-17.458&-313.584\cr 0&0&0&1&5.6338&101.195\cr 0&0&0&0&-0.310553&\bbox[border:1px solid blue]{12.4118}} \right]- 12.4118 \cdot I
= \left[ \matrix{-12.4672&-0.00307557&-1.70959\cdot 10^{-4}&-9.50302\cdot 10^{-6}&-21.2568&-381.819\cr 0.998463&-12.412&-9.48832\cdot 10^{-6}&-5.27422\cdot 10^{-7}&17.4905&314.168\cr 0&0.999995&-12.4118&-2.93174\cdot 10^{-8}&15.1092&271.394\cr 0&0&1&-12.4118&-17.458&-313.584\cr 0&0&0&1&-6.77799&101.195\cr 0&0&0&0&-0.310553&0} \right]=QR
將Q和R對調後相乘再加 A[n,n]\cdot I 得到下一個矩陣A
RQ+ 12.4118 I= \left[ \matrix{-0.134273&-0.0137569&-0.00127716&1.81955&-38.1301&404.512\cr 0.990944&-0.00186405&-1.73055\cdot 10^{-4}&-1.16579&24.4282&-259.152\cr 0&0.999904&-1.7894\cdot 10^{-5}&-1.42457&29.851&-316.681\cr 0&0&0.999999&1.32836&-27.8351&295.295\cr 0&0&0&0.649217&7.37399&-77.7816\cr 0&0&0&0&0.115242&9.4238} \right]
特徵值
\left[ -0.134273,-0.00186405,-1.7894\cdot 10^{-5},1.32836,7.37399,9.4238 \right]

第三輪：
計算A-A[n,n]\cdot I的QR分解
\left[ \matrix{-0.134273&-0.0137569&-0.00127716&1.81955&-38.1301&404.512\cr 0.990944&-0.00186405&-1.73055\cdot 10^{-4}&-1.16579&24.4282&-259.152\cr 0&0.999904&-1.7894\cdot 10^{-5}&-1.42457&29.851&-316.681\cr 0&0&0.999999&1.32836&-27.8351&295.295\cr 0&0&0&0.649217&7.37399&-77.7816\cr 0&0&0&0&0.115242&\bbox[border:1px solid blue]{9.4238}} \right]- 9.4238 \cdot I
= \left[ \matrix{-9.55808&-0.0137569&-0.00127716&1.81955&-38.1301&404.512\cr 0.990944&-9.42567&-1.73055\cdot 10^{-4}&-1.16579&24.4282&-259.152\cr 0&0.999904&-9.42382&-1.42457&29.851&-316.681\cr 0&0&0.999999&-8.09544&-27.8351&295.295\cr 0&0&0&0.649217&-2.04981&-77.7816\cr 0&0&0&0&0.115242&0.0} \right]=QR
將Q和R對調後相乘再加 A[n,n]\cdot I 得到下一個矩陣A
RQ+ 9.4238 I= \left[ \matrix{-0.233098&-0.0373848&-0.210022&5.10063&-52.6375&-427.73\cr 0.972453&-0.0089612&0.0852744&-2.15513&22.236&180.689\cr 0&0.999261&0.175089&-4.36603&45.0504&366.077\cr 0&0&0.874009&3.15576&-32.2784&-262.283\cr 0&0&0&0.312162&7.14176&57.2873\cr 0&0&0&0&-0.0484512&7.75945} \right]
特徵值 \left[ -0.233098,-0.0089612,0.175089,3.15576,7.14176,7.75945 \right]

反覆計算直到全部特徵值 \left[ -1.01,1,3,4,5,6 \right]



要先載入lapack才能使用dgeqrf指令
(%i1)   load(lapack);
0 errors, 0 warnings
0 errors, 0 warnings
0 errors, 0 warnings
0 errors, 0 warnings
0 errors, 0 warnings
0 errors, 0 warnings
0 errors, 0 warnings
0 errors, 0 warnings
0 errors, 0 warnings
(%o1)　C:\maxima-5.39.0\share\maxima\5.39.0_2_g5a49f11_dirty\share\lapack\lapack.mac

有Rayleigh位移的QR分解副程式
(%i2)
QRwithRayleighShift(A,epsilon):=block
([count:0,Q,R,I,n:length(A),eig1,eig2],
eig2:makelist(0,n),
I:ident(n),
do
  (count:count+1,
   [Q,R]:dgeqrf(A-A[n,n]*I),
   A:R.Q+A[n,n]*I,
   eig1:create_list(A[i,i],i,1,n),
   print("特徵值",eig1),
   if lmax(abs(eig1-eig2))<epsilon then
     (print("執行次數",count),
      return(eig1)
     ),
   eig2:eig1
   )
)$

Hessenberg矩陣
(%i3)
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]);
(%o3)　 \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]

計算矩陣A全部的特徵值,誤差10^-5
(%i4)　QRwithRayleighShift(A,10^-5);
特徵值 [-0.0554152121071958,-1.706957479434834 \cdot 10^{-4},-5.274195551407956 \cdot 10^{-7},-1.629643975320505 \cdot 10^{-9},5.633796457744232,12.41178997916011]
特徵值 [-0.1342728973892626,-0.001864054657636061,-1.789401310148264 \cdot 10^{-5},1.328361385039951,7.373990728556812,9.423802732463246]
特徵值 [-0.2330984496644373,-0.00896120187027627,0.1750887851007601,3.155759970239055,7.141761495529245,7.75944940066566]
特徵值 [-0.3419389276255256,-0.03277311167899644,0.7656578204246607,4.323032959391513,6.463362526448887,6.812658733039465]
特徵值 [-0.4516838784880166,-0.04859673999350367,1.632578188284038,4.666773060400056,5.895419068145634,6.295510301651798]
特徵值 [-0.5633422880002747,0.02944794977870835,2.355825335361045,4.60082766258795,5.50449211080626,6.062749229466317]
特徵值 [-0.6788276506780742,0.2239951458769376,2.749431534145878,4.436268986900575,5.255300885689792,6.003831098064898]
特徵值 [-0.7863651279227009,0.4527888848794159,2.906886108392509,4.294225869684157,5.12244896958134,6.000015295385284]
特徵值 [-0.8711542188771304,0.6444491049212555,2.960203806951765,4.196380448349329,5.060120858415357,6.00000000023943]
特徵值 [-0.9289699805732585,0.7799453397954821,2.977380864737738,4.131653967421146,5.02998980861915,5.99999999999975]
特徵值 [-0.9645416477980442,0.8674499800082502,2.98345617185446,4.088593345440747,5.015042150494844,5.99999999999975]
特徵值 [-0.9850455588070135,0.9212586824227689,2.986528792091727,4.059703590286082,5.007554494006693,5.99999999999975]
特徵值 [-0.996420748038723,0.9535030515319471,2.988896306535816,4.040228509856105,5.00379288011511,5.99999999999975]
特徵值 [-1.002605207530562,0.9725812027521856,2.991042001366524,4.027079206073314,5.001902797338794,5.99999999999975]
特徵值 [-1.005940275258888,0.9838100028962309,2.992972662236778,4.01820380067515,5.000953809450982,5.99999999999975]
特徵值 [-1.007739426659477,0.9904122174984673,2.994628810107148,4.012220616788556,5.000477782265562,5.99999999999975]
特徵值 [-1.008716781517061,0.994299741198426,2.995984283944792,4.008193554048924,5.000239202325174,5.99999999999975]
特徵值 [-1.009254404197364,0.9965954520976856,2.997051813148713,4.005487429353162,5.000119709598057,5.99999999999975]
特徵值 [-1.009555347544265,0.9979563073390505,2.997867521069139,4.003671627009723,5.000059892126607,5.99999999999975]
特徵值 [-1.009727505038192,0.9987665417094682,2.998476221504767,4.002454783021602,5.00002995880261,5.99999999999975]
特徵值 [-1.00982848582706,0.999251263181379,2.998922039936204,4.001640198987505,5.000014983722228,5.99999999999975]
特徵值 [-1.009889339001543,0.9995427357842255,2.999243733525915,4.001095376371471,5.000007493320185,5.99999999999975]
特徵值 [-1.009927029226854,0.9997189546363359,2.999473089254277,4.000731238185238,5.000003747151255,5.99999999999975]
特徵值 [-1.009950993097018,0.9998261013041958,2.999635017979005,4.000488000073412,5.00000187374066,5.99999999999975]
特徵值 [-1.009966595539297,0.9998916408332965,2.99974842403746,4.000325593743018,5.000000936925779,5.99999999999975]
特徵值 [-1.009976964039406,0.9999319830257809,2.9998273170086,4.000217195523692,5.000000468481591,5.99999999999975]
特徵值 [-1.009983971927121,0.9999569799022288,2.999881892978214,4.000144864799742,5.000000234247194,5.99999999999975]
特徵值 [-1.009988772855428,0.9999725762149927,2.999919468523729,4.000096610991085,5.00000011712588,5.99999999999975]
特徵值 [-1.009992096532979,0.9999823779538515,2.999945235207018,4.000064424808519,5.000000058563849,5.99999999999975]
特徵值 [-1.009994415923166,0.9999885845273688,2.999962843388455,4.000042958725227,5.000000029282374,5.99999999999975]
特徵值 [-1.009996044151435,0.9999925452102545,2.999974840673362,4.000028643626595,5.000000014641483,5.99999999999975]
特徵值 [-1.009997192214406,0.999995092806369,2.999982994075778,4.000019098011532,5.000000007320986,5.99999999999975]
執行次數 32
(%o4)　 [-1.009997192214406,0.999995092806369,2.999982994075778,4.000019098011532,5.000000007320986,5.99999999999975]

實對稱矩陣
(%i5)
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]);
(%o5)　 \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]

計算矩陣A全部的特徵值,誤差10^-5
(%i6)　QRwithRayleighShift(A,10^-5);
特徵值 [4.687463048832852,2.303233931699447,2.895757871614659,4.180419375119289,1.698867537928739,2.224258234805014]
特徵值 [5.0472697019595,1.567035460671847,3.220795257274317,4.165306321709514,1.524463222953819,2.465130035431002]
特徵值 [5.093280274517875,0.9153126236507791,3.918620820796786,4.066872770913875,1.093783248203376,2.902130261917306]
特徵值 [4.670843575256507,0.6452057194127008,4.657952826984372,4.017888431300075,0.9983637095017406,2.9997457375446]
特徵值 [3.93434857995697,1.130324432377812,4.921565821472631,4.004239302946366,0.9995218632488947,2.99999999999732]
特徵值 [3.001259108881703,2.002038670362651,4.985791453569853,4.001029727196395,0.9998810399893925,3.0]
特徵值 [1.989400649137261,3.001349949453158,4.999025079966701,4.000253914617052,0.9999704068258226,3.0]
特徵值 [1.057592981093276,3.931439082990065,5.000912219500666,4.000063077937745,0.999992638478242,3.0]
特徵值 [0.3196378154597554,4.669606730711599,5.000741560680952,4.000015724365658,0.9999981687820307,3.0]
特徵值 [-0.1980838105378795,5.187667036440498,5.000413303579961,4.000003926041754,0.9999995444756622,3.0]
特徵值 [-0.5311792385543539,5.520973849385642,5.00020452153543,4.000000980947106,0.9999998866861723,3.0]
特徵值 [-0.7337023068465403,5.723606010076015,5.000096079783825,4.000000245174043,0.9999999718126538,3.0]
特徵值 [-0.8526274607253512,5.842583415671053,5.0000439907795,4.000000061286531,0.9999999929882648,3.0]
特徵值 [-0.9210400359561284,5.911020150414187,5.000019871965285,4.000000015320857,0.999999998255797,3.0]
特徵值 [-0.9599296300952491,5.94992071499559,5.000008911703409,4.000000003830127,0.99999999956612,3.0]
特徵值 [-0.9818873306755895,5.971883349239707,5.000003980586289,4.000000000957522,0.999999999892069,3.0]
特徵值 [-0.9942375695225221,5.984235795227988,5.000001774082003,4.00000000023938,0.9999999999731504,3.0]
特徵值 [-1.001169057827592,5.991168268066513,5.000000789707915,4.000000000059845,0.99999999999332,3.0]
特徵值 [-1.005054596592398,5.995054245292558,5.000000351286541,4.00000000001496,0.9999999999983373,3.0]
特徵值 [-1.007231207609367,5.997231051402752,5.000000156203289,4.00000000000374,0.9999999999995857,3.0]
特徵值 [-1.008450042880667,5.998449973437261,5.000000069442574,4.000000000000934,0.9999999999998965,3.0]
特徵值 [-1.009132407625798,5.999132376757524,5.000000030868065,4.000000000000233,0.9999999999999738,3.0]
特徵值 [-1.009514383811725,5.999514370091365,5.00000001372031,4.000000000000057,0.9999999999999929,3.0]
特徵值 [-1.009728193294757,5.999728187196538,5.000000006098205,4.000000000000013,0.9999999999999978,3.0]
特徵值 [-1.009847867730331,5.999847865019944,5.000000002710385,4.000000000000003,0.9999999999999991,3.0]
特徵值 [-1.009914851057516,5.999914849852884,5.000000001204633,4.0,0.9999999999999996,3.0]
特徵值 [-1.009952342051302,5.999952341515906,5.000000000535396,3.999999999999999,0.9999999999999996,3.0]
特徵值 [-1.00997332586102,5.999973325623067,5.000000000237954,3.999999999999999,0.9999999999999996,3.0]
特徵值 [-1.009985070512432,5.999985070406677,5.000000000105756,3.999999999999999,0.9999999999999996,3.0]
特徵值 [-1.009991643987654,5.999991643940653,5.000000000047002,3.999999999999999,0.9999999999999996,3.0]
執行次數 30
(%o6)　 [-1.009991643987654,5.999991643940653,5.000000000047002,3.999999999999999,0.9999999999999996,3.0]

三對角矩陣
(%i7)
A:matrix([3.99889,2.40601,0,0,0,0],
   [2.40601,1.68454,-1.93868,0,0,0],
   [0,-1.93868,2.37077,-1.67177,0,0],
   [0,0,-1.67177,2.35489,1.31944,0],
   [0,0,0,1.31944,3.60161,-0.72545],
   [0,0,0,0,-0.72545,3.97930]);
(%o7)　 \left[ \matrix{3.99889&2.40601&0&0&0&0 \cr 2.40601&1.68454&-1.93868&0&0&0 \cr 0&-1.93868&2.37077&-1.67177&0&0 \cr 0&0&-1.67177&2.35489&1.31944&0 \cr 0&0&0&1.31944&3.60161&-0.725455 \cr 0&0&0&0&-0.725455&        3.9793} \right]

計算矩陣A全部的特徵值,誤差10^-5
(%i8)　QRwithRayleighShift(A,10^-5);
特徵值 [1.72387082004926,3.324540397328139,1.33606542554601,3.693330312406069,3.912210035607812,3.999983009062711]
特徵值 [-0.1957077674919923,4.850881309601665,1.378752309771135,4.000089214192345,3.955987137170426,3.99999779675642]
特徵值 [-0.8127503789372117,4.932944318289477,1.873316361922162,4.037006805484861,3.959485096484287,3.999997796756424]
特徵值 [-0.9592847557581403,4.283345647667558,2.665987856408768,4.040426288835855,3.959527166089535,3.999997796756424]
特徵值 [-0.9952758240675927,3.322858800886055,3.662340297890198,4.040634960426104,3.959443968108812,3.999997796756424]
特徵值 [-1.005287026228012,2.395061618114314,4.600198787928333,4.040618030497149,3.959410792931794,3.999997796756424]
特徵值 [-1.008400216221442,1.734394949564231,5.264000184602603,4.040608164880656,3.959399120417531,3.999997796756424]
特徵值 [-1.009439482743673,1.355559926441152,5.643880858248862,4.040607079913188,3.959393821384052,3.999997796756424]
特徵值 [-1.009799761617996,1.164579046933425,5.835223723597451,4.040608967826792,3.959390226503906,3.999997796756424]
特徵值 [-1.00992696494193,1.07452788620796,5.925402526712892,4.04061168529806,3.959387069966598,3.999997796756424]
特徵值 [-1.00997225893254,1.033407889174767,5.966567928871993,4.040614619747304,3.959384024382054,3.999997796756424]
特徵值 [-1.009988449151831,1.014907085280123,5.985084950921726,4.040617609505592,3.959381006687969,3.999997796756424]
特徵值 [-1.009994246304271,1.006638754385619,5.993359085969478,4.040620613208874,3.95937799598388,3.999997796756424]
特徵值 [-1.009996323664509,1.002954616134844,5.99704530333256,4.04062362041098,3.959374987029705,3.999997796756424]
特徵值 [-1.009997068325025,1.001315290801098,5.998685373765045,4.040626628488697,3.959371978513766,3.999997796756424]
特徵值 [-1.009997335300511,1.000586289686123,5.999414641965082,4.040629636784925,3.95936897010796,3.999997796756424]
特徵值 [-1.009997431023016,1.000262196765743,5.999738830635367,4.040632645135246,3.959365961730242,3.999997796756424]
特徵值 [-1.009997465344805,1.000118133210399,5.999882928519348,4.04063565349854,3.959362953360099,3.999997796756424]
特徵值 [-1.009997477651221,1.000054099107653,5.999946974930223,4.040638661864524,3.959359944992403,3.999997796756424]
特徵值 [-1.009997482063836,1.000025637814941,5.999975440635978,4.040641670230628,3.959356936625871,3.999997796756424]
特徵值 [-1.009997483646035,1.000012987820502,5.999988092212724,4.04064467859621,3.959353928260182,3.999997796756424]
特徵值 [-1.009997484213353,1.000007365422972,5.9999937151776,4.04064768696111,3.959350919895256,3.999997796756424]
執行次數 22
(%o8)　 [-1.009997484213353,1.000007365422972,5.9999937151776,4.04064768696111,3.959350919895256,3.999997796756424]

	(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法	1142次	348次
Rayleigh Quotient Shift	32次	30次
	(3)Hessenberg矩陣 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]);	(4)三對角矩陣 A:matrix([3.99889,2.40601,0,0,0,0],    [2.40601,1.68454,-1.93868,0,0,0],    [0,-1.93868,2.37077,-1.67177,0,0],    [0,0,-1.67177,2.35489,1.31944,0],    [0,0,0,1.31944,3.60161,-0.72545],    [0,0,0,0,-0.72545,3.97930]);
原始的QR法	1142次	444次
Rayleigh Quotient Shift	32次	22次

看得出來Rayleigh Quotient Shift大幅減少QR迭代次數。

[ 本帖最後由 bugmens 於 2017-7-13 19:20 編輯 ]

3-1.減少QR法運算-Wilkinson Shift
取矩陣 A 的右下角 2 \times 2 子矩陣 \left[ \matrix{A[n-1,n-1]& A[n-1,n]\cr A[n,n-1]&A[n,n]} \right] ，
計算特徵方程式 \left| \matrix{A[n-1,n-1]-x& A[n-1,n] \cr A[n,n-1]&A[n,n]-x} \right|=0 ，

x^2-(A[n-1,n-1]+A[n,n])x+(A[n-1,n-1]A[n,n]-A[n-1,n]A[n,n-1])=0

判別式 D=(A[n-1,n-1]+A[n,n])^2-4(A[n-1,n-1]A[n,n]-A[n-1,n]A[n,n-1])

D=(A[n-1,n-1]-A[n,n])^2+4A[n-1,n]A[n,n-1]

D\ge 0 特徵值為實數， D<0 特徵值為複數

(1)若特徵値 \sigma_1,\sigma_2 為實數，看哪個特徵値比較接近 A[n,n] 當作位移値 \alpha 。
if \left| \sigma_1-A[n,n] \right|<\left| \sigma_2-A[n,n] \right| then
　( \alpha=\sigma_1 )
else
　( \alpha=\sigma_2 )
計算 A-\alpha I 的QR分解
計算 RQ+\alpha I 得到A矩陣

(2)若特徵值 \sigma_1,\sigma_2 為複數，使用Implicit Double Shift避免複數運算
但 \sigma_1,\sigma_2 為共軛複數，改以 \sigma,\overline{\sigma} 表示
個別計算 A-\sigma I 和 A-\overline{\sigma}I 的QR分解都會遇到複數運算，故將兩次Shift運算結合一起變成實數矩陣再做QR分解
計算 (A-\sigma I)(A-\overline{\sigma}I)=A^2-(\sigma+\overline{\sigma})A+\sigma \overline{\sigma}I=A^2-2Re(\sigma)A+\left| \sigma \right|^2 I 為實數矩陣的QR分解
計算Q^T \cdot A \cdot Q得到A矩陣

參考資料：Implicit Double Shift QR Algorithm
https://books.google.com.tw/book ... ge&q&f=true



要先載入lapack才能使用dgeqrf指令
(%i1)　load(lapack);
0 errors, 0 warnings
0 errors, 0 warnings
0 errors, 0 warnings
0 errors, 0 warnings
0 errors, 0 warnings
0 errors, 0 warnings
0 errors, 0 warnings
0 errors, 0 warnings
0 errors, 0 warnings
(%o1)　C:\maxima-5.40.0\share\maxima\5.40.0\share\lapack\lapack.mac

Hessenberg矩陣
(%i2)
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]);
(%o2)　 \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]

設定小數點底下第6位四捨五入
(%i3)　fpprintprec:6;
(%o3)　6

右下角2*2矩陣特徵值為實數的Wilkinson Shift
(%i4)
WilkinsonShift(A,eig2by2):=block
(print("特徵值為實數",eig2by2),
if abs(eig2by2[1]-A[n,n])<abs(eig2by2[2]-A[n,n]) then
   (print("特徵值",eig2by2[1],"比較接近A[n,n]=",A[n,n]),
    alpha:eig2by2[1]
   )
else
   (print("特徵值",eig2by2[2],"比較靠近A[n,n]=",A[n,n]),
    alpha:eig2by2[2]
   ),
print("計算A-",alpha,"*I的QR分解=QR"),
[Q,R]:dgeqrf(A-alpha*I),
print("計算A=RQ+",alpha,"*I=",A:R.Q+alpha*I)
)$

右下角2*2矩陣特徵值為複數的Implicit Double Shift
(%i5)
ImplicitDoubleShift(A,eig2by2):=block
(print("特徵值為複數δ=",eig2by2),
print("計算A^2-2Re(δ)A+|δ|^2*I=A^2-2*",realpart(eig2by2[1]),"A+",cabs(eig2by2[1])^2,"I的QR分解"=QR),
[Q,R]:dgeqrf(A.A-2*realpart(eig2by2[1])*A+cabs(eig2by2[2])^2*I),
print("計算A=Q^T.A.Q=",A:transpose(Q).A.Q)
)$

具有WilkinsonShift的QR法副程式
(%i6)
QRwithWilkinsonShift(A,epsilon):=block
([count:0,Q,R,I,n:length(A),subA,eig1,eig2,eigenvalue],
eig2:makelist(0,n),
I:ident(n),
do
  (count:count+1,
   print("第",count,"輪"),
   subA:genmatrix(lambda([i,j],A[n-2+i,n-2+j]),2,2),
   print("A的右下角2*2矩陣",subA),
   print(subA-x*ident(2),"=0,特徵方程式",expand(determinant(subA-x*ident(2))),"=0"),
   eig2by2:dgeev(subA),/*特徵值*/
   if (subA[1,1]-subA[2,2])^2+4*subA[1,2]*subA[2,1]>=0 then/*特徵方程式判別式≧0*/
      (A:WilkinsonShift(A,eig2by2[1])/*特徵值為實數*/
      )
    else
      (A:ImplicitDoubleShift(A,eig2by2[1])/*特徵值為複數*/
      ),
    eig1:create_list(A[i,i],i,1,n),
    print("特徵值",eig1),
    if lmax(abs(eig1-eig2))<epsilon then
      (print("執行次數",count),
      return(eig1)
      ),
    eig2:eig1,
    print("－－－－－－－－－－－")
  )
)$

計算矩陣A全部的特徵值,誤差10^-5
(%i7)　QRwithWilkinsonShift(A,10^-5);
第1輪
A 的右下角2*2矩陣 \left[ \matrix{0 & -117.81 \cr 1 & 17.99} \right]
\left[ \matrix{-x & -117.81 \cr 1 & 17.99-x } \right] =0 ,特徵方程式x^2-17.99*x+117.81=0
特徵值為複數 \sigma=[6.07453*\%i+8.995,8.995-6.07453*\%i ]
計算 A^2-2Re(\sigma)A+|\sigma|^2*I=A^2-2*8.995A+117.81I 的QR分解 =QR
計算 A=Q^T.A.Q=\left[ \matrix{ -0.15048 & -0.0145152 & -9.53934*10^{-4} & 3.48119 & -21.2753 & -409.554 \cr 0.988613 & -0.00220941 & -1.45201*10^{-4} & -2.10275 & 12.8509 & 247.382 \cr -2.88398*10^{-17} & 0.999892 & -1.41689*10^{-5} & -2.78352 & 17.0114 & 327.473 \cr -1.21431*10^{-17} & -2.94903*10^{-17} & 1.0 & 2.48535 & -15.1892 & -292.395 \cr -1.0842*10^{-18} & -1.73472*10^{-18} & 2.52077*10^{-17} & 0.381684 & 3.84825 & 72.1532 \cr -2.96462*10^{-20} & -1.6263*10^{-19} & 0.0 & -2.1684*10^{-19} & -0.382882 & 11.8091 } \right]
特徵值 [-0.15048,-0.00220941,-1.41689*10^{-5},2.48535,3.84825,11.8091]
－－－－－－－－－－－
第2輪
A 的右下角2*2矩陣 \left[ \matrix{3.84825 & 72.1532 \cr -0.382882 & 11.8091 } \right]
\left[ \matrix{3.84825-x & 72.1532 \cr -0.382882 & 11.8091-x } \right] =0 ,特徵方程式x^2-15.6573*x+73.0706=0
特徵值為複數 \sigma=[3.43255*\%i+7.82867,7.82867-3.43255*\%i ]
計算 A^2-2Re(\sigma)A+|\sigma|^2*I=A^2-2*7.82867A+73.0706I 的QR分解 =QR
計算 A=Q^T.A.Q=\left[ \matrix{ -0.340459 & -0.0187684 & -0.846134 & 7.56341 & -44.0982 & 445.365 \cr 0.94026 & -0.034987 & 0.14736 & -1.37467 & 8.01098 & -80.9063 \cr -1.77764*10^{-17} & 0.94948 & 0.787709 & -6.90939 & 40.2584 & -406.585 \cr -1.04869*10^{-17} & -2.56953*10^{-17} & 0.571107 & 3.93811 & -22.4175 & 226.378 \cr -2.61695*10^{-19} & -3.25585*10^{-21} & -1.39862*10^{-17} & 0.155215 & 4.98351 & -48.8466 \cr -4.16098*10^{-20} & 1.22237*10^{-19} & 0.0 & 1.56125*10^{-17} & 0.136995 & 8.65612} \right]
特徵值 [-0.340459,-0.034987,0.787709,3.93811,4.98351,8.65612 ]
－－－－－－－－－－－
第3輪
A 的右下角2*2矩陣 \left[ \matrix{4.98351 & -48.8466 \cr0.136995 & 8.65612 } \right]
\left[ \matrix{4.98351-x & -48.8466 \cr 0.136995 & 8.65612-x } \right] =0 ,特徵方程式x^2-13.6396*x+49.8296=0
特徵值為複數 \sigma=[1.82201*\%i+6.81982,6.81982-1.82201*\%i ]
計算 A^2-2Re(\sigma)A+|\sigma|^2*I=A^2-2*6.81982A+49.8296I 的QR分解 =QR
計算 A=Q^T.A.Q=\left[ \matrix{ -0.550722 & 0.160314 & -2.0241 & 10.9002 & -59.6766 & -443.499 \cr 0.843302 & 0.0181879 & -0.597236 & 3.09386 & -16.9466 & -125.941 \cr -3.643*10^{-17} & 0.660459 & 2.00489 & -10.3831 & 56.7744 & 421.933 \cr -6.91901*10^{-18} & 7.53929*10^{-17} & 0.225077 & 4.30867 & -22.6477 & -168.297 \cr 3.45526*10^{-19} & 1.51358*10^{-18} & -7.81981*10^{-18} & 0.0606534 & 5.15209 & 36.991 \cr -1.26193*10^{-20} & 4.42551*10^{-20} & 1.0842*10^{-19} & 4.33681*10^{-18} & -0.0442203 & 7.05688 } \right]
特徵值 [-0.550722,0.0181879,2.00489,4.30867,5.15209,7.05688]
－－－－－－－－－－－
第4輪
A 的右下角2*2矩陣 \left[ \matrix{5.15209 & 36.991 \cr -0.0442203 & 7.05688 } \right]
\left[ \matrix{5.15209-x & 36.991 \cr -0.0442203 & 7.05688-x } \right] =0 ,特徵方程式x^2-12.209*x+37.9934=0
特徵值為複數 \sigma=[0.853639*\%i+6.10448,6.10448-0.853639*\%i]
計算 A^2-2Re(\sigma)A+|\sigma|^2*I=A^2-2*6.10448A+37.9934I 的QR分解 =QR
計算 A=Q^T.A.Q=\left[ \matrix{ -0.769656 & 0.352202 & -2.62459 & 11.8663 & -61.7088 & 395.945 \cr 0.659697 & 0.377261 & -2.01053 & 8.94582 & -46.53 & 298.552 \cr -5.61889*10^{-17} & 0.295707 & 2.71313 & -11.7142 & 60.8213 & -390.256 \cr -3.03024*10^{-18} & 1.56599*10^{-17} & 0.0767723 & 4.28338 & -21.1565 & 135.744 \cr -3.18495*10^{-20} & -1.6239*10^{-19} & -1.61032*10^{-18} & 0.0181863 & 5.08056 & -31.3718 \cr 4.51695*10^{-22} & -4.18436*10^{-21} & -3.38813*10^{-20} & 4.33681*10^{-19} & 0.011151 & 6.30532 } \right]
特徵值 [-0.769656,0.377261,2.71313,4.28338,5.08056,6.30532]
－－－－－－－－－－－
第5輪
A 的右下角2*2矩陣 \left[ \matrix{5.08056 & -31.3718 \cr0.011151 & 6.30532 } \right]
\left[ \matrix{5.08056-x & -31.3718 \cr 0.011151 & 6.30532-x } \right] =0 ,特徵方程式x^2-11.3859*x+32.3844=0
特徵值為實數 [5.53426,5.85162]
特徵值5.85162比較靠近 A[n,n]=6.30532
計算 A-5.85162*I 的QR分解 =QR
計算 A=RQ+5.85162*I=\left[ \matrix{ -0.85821 & 0.360367 & -2.65291 & 11.4251 & -60.8763 & -363.718 \cr 0.537404 & 0.579626 & -2.62892 & 11.1957 & -59.6618 & -356.461 \cr -3.22653*10^{-17} & 0.177737 & 2.86104 & -11.7597 & 62.5505 & 373.727 \cr -9.36941*10^{-19} & 4.63981*10^{-18} & 0.0431351 & 4.21815 & -21.2607 & -127.025 \cr -2.58762*10^{-21} & -5.26611*10^{-20} & -7.55784*10^{-19} & 0.00961362 & 5.23107 & 30.1214 \cr 7.24317*10^{-24} & -8.11741*10^{-23} & -1.063*10^{-21} & 3.12971*10^{-20} & 0.00122958 & 5.95832 } \right]
特徵值 [-0.85821,0.579626,2.86104,4.21815,5.23107,5.95832]
－－－－－－－－－－－
第6輪
A 的右下角2*2矩陣 \left[ \matrix{5.23107 & 30.1214 \cr 0.00122958 & 5.95832 } \right]
\left[ \matrix{5.23107-x & 30.1214 \cr 0.00122958 & 5.95832-x } \right] =0 ,特徵方程式x^2-11.1894*x+31.1313=0
特徵值為實數 [5.18328,6.00611]
特徵值6.00611比較靠近 A[n,n]=5.95832
計算 A-6.00611*I 的QR分解 =QR
計算 A=RQ+6.00611*I=\left[ \matrix{ -0.919309 & 0.32341 & -2.5702 & 10.7607 & -55.5116 & 334.865 \cr 0.420282 & 0.730591 & -3.06377 & 12.7235 & -65.6431 & 395.981 \cr -1.66882*10^{-17} & 0.106845 & 2.92217 & -11.6734 & 60.1084 & -362.601 \cr -2.71431*10^{-19} & 1.42958*10^{-18} & 0.0258404 & 4.16855 & -20.2876 & 122.381 \cr -3.60428*10^{-22} & -8.17873*10^{-21} & -1.87791*10^{-19} & 0.00434508 & 5.08827 & -29.5139 \cr -6.71049*10^{-27} & 9.51864*10^{-26} & 2.0198*10^{-24} & -1.14183*10^{-22} & -8.94925*10^{-6} & 5.99973 } \right]
特徵值 [-0.919309,0.730591,2.92217,4.16855,5.08827,5.99973 ]
－－－－－－－－－－－

省略第7輪到第27輪執行過程

－－－－－－－－－－－
第28輪
A 的右下角2*2矩陣 \left[ \matrix{5.0 & 28.9828 \cr 0.0 & 6.0 } \right]
\left[ \matrix{5.0-x & 28.9828 \cr 0.0 & 6.0-x} \right] =0 ,特徵方程式x^2-11.0*x+30.0=0
特徵值為實數 [5.0,6.0]
特徵值6.0比較靠近 A[n,n]=6.0
計算 A-6.0*I 的QR分解 =QR
計算 A=RQ+6.0*I=\left[ \matrix{ -1.01 & 0.0214797 & -1.96762 & 7.69405 & 38.9703 & 235.303 \cr 2.87689*10^{-4} & 0.999994 & -4.0136 & 15.6682 & 79.3635 & 479.196 \cr -1.78149*10^{-25} & 1.43282*10^{-6} & 2.99997 & -11.2221 & -56.7202 & -342.484 \cr -2.56298*10^{-31} & 1.61709*10^{-27} & 2.58843*10^{-6} & 4.00003 & 18.9798 & 114.601 \cr 8.39634*10^{-41} & 2.80728*10^{-36} & 4.40677*10^{-30} & -1.10018*10^{-9} & 5.0 & 28.9828 \cr 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 6.0} \right]
特徵值 [-1.01,0.999994,2.99997,4.00003,5.0,6.0]
－－－－－－－－－－－
第29輪
A的右下角2*2矩陣 \left[ \matrix{5.0 & 28.9828 \cr 0.0 & 6.0} \right]
\left[ \matrix{5.0-x & 28.9828 \cr 0.0 & 6.0-x} \right] =0 ,特徵方程式x^2-11.0*x+30.0=0
特徵值為實數 [5.0,6.0]
特徵值6.0比較靠近 A[n,n]=6.0
計算 A-6.0*I 的QR分解 =QR
計算 A=RQ+6.0*I=\left[ \matrix{-1.01 & 0.0213978 & -1.96746 & 7.6934 & -38.967 & -235.283 \cr 2.052*10^{-4} & 0.999996 & -4.01369 & 15.6686 & -79.3651 & -479.206 \cr -7.62419*10^{-26} & 8.59697*10^{-7} & 2.99998 & -11.2221 & 56.7202 & 342.484 \cr -7.31228*10^{-32} & 6.46829*10^{-28} & 1.72559*10^{-6} & 4.00002 & -18.9798 & -114.601 \cr -1.19777*10^{-41} & -5.61456*10^{-37} & -1.46891*10^{-30} & 5.50096*10^{-10} & 5.0 & 28.9828 \cr 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 6.0} \right]
特徵值 [-1.01,0.999996,2.99998,4.00002,5.0,6.0 ]
執行次數29
(%o7)　 [-1.01,0.999996,2.99998,4.00002,5.0,6.0]

[ 本帖最後由 bugmens 於 2017-7-19 07:15 編輯 ]