18 12
發新話題
打印

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

推到噗浪
推到臉書
1.求QR分解的方法--Gram-Schmidt

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



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


本文使用的範例
(%i2)
A:matrix([-4,0,1],
              [-1,3,1],
              [1,2,5]);

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

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

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

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

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

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

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

[ 本帖最後由 bugmens 於 2016-10-29 03:37 PM 編輯 ]

TOP

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

TOP

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

TOP

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

TOP

假設矩陣\(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 編輯 ]

TOP

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] \)

TOP

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

TOP

拿之前的例子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 編輯 ]

TOP

 18 12
發新話題