23 123
發新話題
打印

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

推到噗浪
推到臉書
3-2.減少QR法運算-對稱矩陣的Wilkinson Shift
若矩陣\(A\)是對稱矩陣,則計算過程可以再簡化。
矩陣\(A\)的右下角\(2\times2\)子矩陣\( \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}&b_{n-1}\cr b_{n-1}&a_n} \right] \)代替
\( \left[ \matrix{a_{n-1}-x & b_{n-1}\cr b_{n-1}&a_n-x} \right]=0 \),特徵方程式\( x^2-(a_{n-1}+a_n)x+(a_{n-1}a_n-b_{n-1}^2)=0 \)
特徵值\( \displaystyle \lambda_{1,2}=\frac{a_{n-1}+a_n}{2}\pm \sqrt{\left( \frac{a_{n-1}-a_{n}}{2}\right)^2+b_{n-1}^2} \)
令\( \displaystyle \delta=\frac{a_{n-1}-a_{n}}{2} \)
\( \displaystyle \lambda_{1,2}=a_n+\delta\pm \sqrt{\delta+b_{n-1}^2} \)
\( \lambda_{1,2}=a_n+x_{1,2} \),其中\( x_{1,2} \)是方程式\( x^2-2\delta x-b_{n-1}^2=0 \)的根
設\( x_1 \)是方程式比較大的根,\( \left| x_1 \right|=sign(\delta)(\left| \delta \right|+\sqrt{\delta^2+b_{n-1}^2})=\delta+sign(\delta)\sqrt{\delta^2+b_{n-1}^2} \)
利用根與係數的關係\( x_1 \cdot x_2=-b_{n-1}^2 \),求比較小的根\( \displaystyle x_2=\frac{-b_{n-1}^2}{\delta+sign(\delta)\sqrt{\delta^2+b_{n-1}^2}} \)
位移值\( \displaystyle \alpha=a_n-\frac{b_{n-1}^2}{\delta+sign(\delta)\sqrt{\delta^2+b_{n-1}^2}} \)


當\( a_{n-1}=a_n \),\( \delta=0 \)時
\( \lambda_{1,2}=a_n \pm \left| b_{n-1} \right| \)
位移值\( \alpha=a_n-\left| b_{n-1} \right| \)

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


以簡單的整數為例子
矩陣\(A\)的右下角\(2\times2\)子矩陣\( \left[ \matrix{a_{n-1}&b_{n-1}\cr b_{n-1}&a_n} \right]=\left[ \matrix{1&2 \cr 2&3} \right] \)
\( \left[ \matrix{1-x&2 \cr 2 &3-x} \right]=0 \),特徵方程式\( x^2-4x-1=0 \)
特徵值\( \lambda_{1,2}=2 \pm \sqrt{5} \)
\( 2+\sqrt{5} \)比較接近\( A[n,n]=3 \)
位移值\( \alpha=2+\sqrt{5} \)


代公式
\( \displaystyle \delta=\frac{a_{n-1}-a_{n}}{2}=\frac{1-3}{2}=-1 \)
\( sign(\delta)=-1 \)
位移值\( \displaystyle \alpha=a_n-\frac{b_{n-1}^2}{\delta+sign(\delta)\sqrt{\delta^2+b_{n-1}^2}}=3-\frac{2^2}{-1+(-1)\sqrt{(-1)^2+2^2}}=2+\sqrt{5} \)


\(a_{n-1}=a_{n}\)的例子
矩陣\(A\)的右下角\( 2 \times 2 \)子矩陣\( \left[ \matrix{a_{n-1} & b_{n-1}\cr b_{n-1}&a_n} \right]=\left[ \matrix{1&2 \cr 2&1} \right] \)
\( \left[ \matrix{1-x& 2 \cr 2 &1-x} \right]=0 \),特徵方程式\(x^2-2x-3=0\)
特徵值\(\lambda_{1,2}=-1,3\)
\(-1,3\)和\(A[n,n]=1\)距離相同
位移值\(\alpha=-1\)或\(\alpha=3\)


代公式
\( \displaystyle \delta=\frac{a_{n-1}-a_{n}}{2}=0 \)
\( \lambda_{1,2}=a_n \pm \left| b_{n-1} \right|=-1,3 \)
位移值\( \alpha=a_n-\left| b_{n-1} \right|=-1 \)


位移值\(\alpha\)的maxima程式
delta: (A[n-1,n-1]-A[n,n])/2,
if delta=0 then
  (alpha:A[n,n]-abs(A[n-1,n])
  )
else
  (alpha:float(A[n,n]-A[n-1,n]^2/(delta+signum(delta)*sqrt(delta^2+A[n-1,n]^2)))
  ),



要先載入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

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

對稱矩陣的WilkinsonShift
(%i3)
SymmetricWilkinsonShift(A):=block
([subA,delta,alpha],   
subA:genmatrix(lambda([i,j],A[n-2+i,n-2+j]),2,2),
print("A的右下角2*2矩陣",subA),
delta: (A[n-1,n-1]-A[n,n])/2,
if delta=0 then
  (alpha:A[n,n]-abs(A[n-1,n])
  )
else
  (alpha:float(A[n,n]-A[n-1,n]^2/(delta+signum(delta)*sqrt(delta^2+A[n-1,n]^2)))
  ),
print("計算A-",alpha,"*I的QR分解=QR"),
[Q,R]:dgeqrf(A-alpha*I),
print("計算A=RQ+",alpha,"*I=",A:R.Q+alpha*I)
)$


具有WilkinsonShift的QR法副程式
(%i4)
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,"輪"),
   A:SymmetricWilkinsonShift(A),
   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("-----------")
  )
)$


對稱矩陣
(%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) QRwithWilkinsonShift(A,10^-5);
第1輪
\(A\)的右下角2*2矩陣\( \left[ \matrix{
374/225 & 901/450 \cr
901/450 & 2099/900} \right] \)
計算\( A-4.02728*I \)的QR分解\( =QR \)
計算\( A=RQ+4.02728*I=\left[ \matrix{
1.62809 & 2.80636 & 0.881862 & -0.487506 & 0.72242 & 0.0113621 \cr
2.80636 & 3.91942 & -0.110835 & 0.32049 & -0.475787 & 0.00735244 \cr
0.881862 & -0.110835 & 5.07773 & 0.0843398 & -0.127243 & 0.0369233 \cr
-0.487506 & 0.32049 & 0.0843398 & 2.89772 & 0.153266 & -0.0269407 \cr
0.72242 & -0.475787 & -0.127243 & 0.153266 & 0.467693 & -0.0602476 \cr
0.0113621 & 0.00735244 & 0.0369233 & -0.0269407 & -0.0602476 & 3.99934} \right] \)

特徵值\( [1.62809,3.91942,5.07773,2.89772,0.467693,3.99934] \)
-----------
第2輪
\(A\)的右下角2*2矩陣\( \left[ \matrix{
0.467693 & -0.0602476 \cr
-0.0602476 & 3.99934} \right] \)
計算\( A-4.00037*I \)的QR分解\( =QR \)
計算\( A=RQ+4.00037*I=\left[ \matrix{
-0.229037 & -1.82178 & -0.262006 & -0.138638 & 0.613611 & -1.08914*10^{-6} \cr
-1.82178 & 5.40498 & 0.165062 & -0.0407398 & 0.180314 & 2.27869*10^{-6} \cr
-0.262006 & 0.165062 & 5.05008 & -0.00508227 & 0.022494 & 1.13371*10^{-5} \cr
-0.138638 & -0.0407398 & -0.00508227 & 2.99446 & 0.0245179 & 9.86431*10^{-6} \cr
0.613611 & 0.180314 & 0.022494 & 0.0245179 & 0.769524 & 6.97715*10^{-6} \cr
-1.08914*10^{-6} & 2.27869*10^{-6} & 1.13371*10^{-5} & 9.86431*10^{-6} & 6.97715*10^{-6} & 4.0} \right] \)

特徵值\( [-0.229037,5.40498,5.05008,2.99446,0.769524,4.0] \)
-----------
第3輪
\(A\)的右下角2*2矩陣\( \left[ \matrix{
0.769524 & 6.97715*10^{-6} \cr
6.97715*10^{-6} & 4.0} \right] \)
計算\( A-4.0*I \)的QR分解\( =QR \)
計算\( A=RQ+4.0*I=\left[ \matrix{
-0.818759 & 0.830858 & 0.0586549 & -0.0309183 & 0.41086 & -2.45934*10^{-16} \cr
0.830858 & 5.88358 & 0.115146 & 0.00377153 & -0.0501182 & 5.74485*10^{-16} \cr
0.0586549 & 0.115146 & 5.01463 & 2.32382*10^{-4} & -0.00308803 & 2.84037*10^{-16} \cr
-0.0309183 & 0.00377153 & 2.32382*10^{-4} & 2.99975 & 0.00328852 & -3.02836*10^{-16} \cr
0.41086 & -0.0501182 & -0.00308803 & 0.00328852 & 0.910799 & -2.75267*10^{-16} \cr
5.58974*10^{-18} & 2.99876*10^{-17} & 2.54508*10^{-16} & -2.35496*10^{-16} & -5.41112*10^{-17} & 4.0} \right] \)

特徵值\( [-0.818759,5.88358,5.01463,2.99975,0.910799,4.0] \)
-----------

省略第4輪到第10輪執行過程

-----------
第11輪
\(A\)的右下角2*2矩陣\( \left[ \matrix{
0.999929 & -2.81267*10^{-16} \cr
0.0 & 4.0} \right] \)
計算\( A-4.0*I \)的QR分解\( =QR \)
計算\( A=RQ+4.0*I=\left[ \matrix{
-1.00997 & 5.64188*10^{-4} & 1.55186*10^{-7} & -8.21073*10^{-8} & 0.00715985 & -2.63617*10^{-16} \cr
5.64188*10^{-4} & 6.0 & 4.8144*10^{-4} & 6.60776*10^{-12} & -5.7625*10^{-7} & 5.09567*10^{-16} \cr
1.55186*10^{-7} & 4.8144*10^{-4} & 5.0 & 1.41015*10^{-15} & -1.09379*10^{-10} & -3.51744*10^{-17} \cr
-8.21073*10^{-8} & 6.6083*10^{-12} & 1.23425*10^{-15} & 3.0 & 1.33342*10^{-10} & -6.4829*10^{-17} \cr
0.00715985 & -5.7625*10^{-7} & -1.09378*10^{-10} & 1.33342*10^{-10} & 0.999974 & -2.81897*10^{-16} \cr
0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 4.0} \right] \)

特徵值\( [-1.00997,6.0,5.0,3.0,0.999974,4.0] \)
-----------
第12輪
\(A\)的右下角2*2矩陣\( \left[ \matrix{
0.999974 & -2.81897*10^{-16} \cr
0.0 & 4.0} \right] \)
計算\( A-4.0*I \)的QR分解\( =QR \)
計算\( A=RQ+4.0*I=\left[ \matrix{
-1.00999 & -2.25226*10^{-4} & -3.09753*10^{-8} & -1.63887*10^{-8} & -0.00428737 & 2.63271*10^{-16} \cr
-2.25226*10^{-4} & 6.0 & 2.4072*10^{-4} & -5.27096*10^{-13} & -1.3775*10^{-7} & -5.09529*10^{-16} \cr
-3.09753*10^{-8} & 2.4072*10^{-4} & 5.0 & 4.69149*10^{-16} & 7.14311*10^{-11} & 3.5297*10^{-17} \cr
-1.63887*10^{-8} & -5.26559*10^{-13} & 2.93119*10^{-16} & 3.0 & 2.22622*10^{-11} & 6.4829*10^{-17} \cr
-0.00428737 & -1.3775*10^{-7} & 7.14309*10^{-11} & 2.22623*10^{-11} & 0.999991 & -2.82273*10^{-16} \cr
0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 4.0} \right] \)

特徵值\( [-1.00999,6.0,5.0,3.0,0.999991,4.0] \)
-----------
第13輪
\(A\)的右下角2*2矩陣\( \left[ \matrix{
0.999991 & -2.82273*10^{-16} \cr
0.0 & 4.0} \right] \)
計算\( A-4.0*I \)的QR分解\( =QR \)
計算\( A=RQ+4.0*I=\left[ \matrix{
-1.01 & 8.99109*10^{-5} & 6.18255*10^{-9} & -3.27108*10^{-9} & 0.00256729 & -2.63006*10^{-16} \cr
8.99109*10^{-5} & 6.0 & 1.2036*10^{-4} & 4.14181*10^{-14} & -3.29283*10^{-8} & 5.09536*10^{-16} \cr
6.18255*10^{-9} & 1.2036*10^{-4} & 5.0 & -1.77861*10^{-16} & 2.6205*10^{-10} & -3.53584*10^{-17} \cr
-3.27108*10^{-9} & 4.19552*10^{-14} & -3.53956*10^{-16} & 3.0 & -1.1726*10^{-10} & -6.4829*10^{-17} \cr
0.00256729 & -3.29283*10^{-8} & 2.6205*10^{-10} & -1.1726*10^{-10} & 0.999997 & -2.82498*10^{-16} \cr
0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 4.0} \right] \)

特徵值\( [-1.01,6.0,5.0,3.0,0.999997,4.0] \)
執行次數13
(%o6) \( [-1.01,6.0,5.0,3.0,0.999997,4.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) QRwithWilkinsonShift(A,10^-5);
第1輪
\(A\)的右下角2*2矩陣\( \left[ \matrix{
3.60161 & -0.72545\cr
-0.72545 & 3.9793} \right] \)

計算\( A-4.54008*I\)的QR分解\(=QR\)
計算\(A=RQ+4.54008*I=\left[ \matrix{
0.765738 & -2.52898 & -2.22045*10^{-16} & 1.66533*10^{-16} & 1.11022*10^{-16} & -2.22045*10^{-16}\cr
-2.52898 & 3.9155 & 1.8145 & 1.11022*10^{-16} & 5.55112*10^{-17} & -1.11022*10^{-16}\cr
0.0 & 1.8145 & 1.78123 & -1.07006 & 1.66533*10^{-16} & -4.996*10^{-16}\cr
0.0 & 0.0 & -1.07006 & 2.8407 & -0.672845 & -2.22045*10^{-16}\cr
0.0 & 0.0 & 0.0 & -0.672845 & 4.37745 & -0.527042\cr
0.0 & 0.0 & 0.0 & 0.0 & -0.527042 & 4.30939} \right] \)

特徵值\( [0.765738,3.9155,1.78123,2.8407,4.37745,4.30939] \)
-----------
第2輪
\(A\)的右下角2*2矩陣\( \left[ \matrix{
4.37745 & -0.527042\cr
-0.527042 & 4.30939} \right] \)
計算\( A-3.81528*I\)的QR分解\(=QR\)
計算\(A=RQ+3.81528*I=\left[ \matrix{
-0.436081 & 1.58352 & -1.41226*10^{-16} & 2.11044*10^{-16} & -1.5642*10^{-18} & 2.62182*10^{-16}\cr
1.58352 & 4.72567 & -1.92785 & 1.32422*10^{-16} & 1.76772*10^{-16} & 2.95571*10^{-16}\cr
0.0 & -1.92785 & 1.74151 & -0.364505 & 2.62527*10^{-16} & 4.26581*10^{-16}\cr
0.0 & 0.0 & -0.364505 & 3.28272 & 0.738764 & 1.11022*10^{-16}\cr
0.0 & 0.0 & 0.0 & 0.738764 & 4.6657 & 0.123502\cr
0.0 & 0.0 & 0.0 & 0.0 & 0.123502 & 4.01048} \right] \)

特徵值\( [-0.436081,4.72567,1.74151,3.28272,4.6657,4.01048] \)
-----------
第3輪
\(A\)的右下角2*2矩陣\( \left[ \matrix{
4.6657 & 0.123502\cr
0.123502 & 4.01048} \right] \)
計算\( A-3.98797*I\)的QR分解\(=QR\)
計算\(A=RQ+3.98797*I=\left[ \matrix{
-0.854732 & -0.770321 & -4.23756*10^{-16} & 1.93093*10^{-16} & 1.10126*10^{-16} & -1.34706*10^{-16}\cr
-0.770321 & 4.36891 & 2.32945 & 1.44982*10^{-17} & -2.19437*10^{-16} & 1.36724*10^{-16}\cr
0.0 & 2.32945 & 2.47901 & -0.131906 & 4.97574*10^{-16} & -4.90434*10^{-16}\cr
0.0 & 0.0 & -0.131906 & 3.3263 & -0.743098 & -1.66533*10^{-16}\cr
0.0 & 0.0 & 0.0 & -0.743098 & 4.67051 & -0.00148842\cr
0.0 & 0.0 & 0.0 & 0.0 & -0.00148842 & 4.0} \right] \)

特徵值\( [-0.854732,4.36891,2.47901,3.3263,4.67051,4.0] \)
-----------

省略第4輪到第18輪執行過程

-----------
第19輪
\(A\)的右下角2*2矩陣\( \left[ \matrix{
4.67044 & -1.3167*10^{-16}\cr
0.0 & 4.0} \right] \)
計算\( A-4.0*I\)的QR分解\(=QR\)
計算\(A=RQ+4.0*I=\left[ \matrix{
-1.01 & 1.22187*10^{-4} & 6.31425*10^{-16} & 1.99522*10^{-16} & -1.43708*10^{-17} & -3.9486*10^{-17}\cr
1.22187*10^{-4} & 1.00003 & 0.0114425 & 1.0112*10^{-16} & -4.9102*10^{-16} & -5.11954*10^{-16}\cr
0.0 & 0.0114425 & 5.99997 & 1.10762*10^{-6} & 1.03762*10^{-16} & 1.39068*10^{-16}\cr
0.0 & 0.0 & 1.10762*10^{-6} & 3.32956 & 0.741963 & -1.46328*10^{-16}\cr
0.0 & 0.0 & 0.0 & 0.741963 & 4.67043 & -3.44542*10^{-17}\cr
0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 4.0} \right] \)

特徵值\( [-1.01,1.00003,5.99997,3.32956,4.67043,4.0] \)
-----------
第20輪
\(A\)的右下角2*2矩陣\( \left[ \matrix{
4.67043 & -3.44542*10^{-17}\cr
0.0 & 4.0} \right] \)
計算\( A-4.0*I\)的QR分解\(=QR\)
計算\(A=RQ+4.0*I=\left[ \matrix{
-1.01 & 7.31656*10^{-5} & 6.29722*10^{-16} & 1.44419*10^{-16} & -1.3841*10^{-16} & 3.94735*10^{-17}\cr
7.31656*10^{-5} & 1.00001 & -0.00762834 & 4.32693*10^{-16} & 2.54117*10^{-16} & 5.12481*10^{-16}\cr
0.0 & -0.00762834 & 5.99999 & -5.53806*10^{-7} & 1.35945*10^{-17} & -1.37114*10^{-16}\cr
0.0 & 0.0 & -5.53806*10^{-7} & 3.32957 & 0.741964 & 7.25401*10^{-17}\cr
0.0 & 0.0 & 0.0 & 0.741964 & 4.67043 & -1.3167*10^{-16}\cr
0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 4.0} \right] \)

特徵值\( [-1.01,1.00001,5.99999,3.32957,4.67043,4.0] \)
-----------
第21輪
\(A\)的右下角2*2矩陣\( \left[ \matrix{
4.67043 & -1.3167*10^{-16}\cr
0.0 & 4.0} \right] \)
計算\( A-4.0*I\)的QR分解\(=QR\)
計算\(A=RQ+4.0*I=\left[ \matrix{
-1.01 & 4.38117*10^{-5} & 6.30883*10^{-16} & 1.99518*10^{-16} & -1.43522*10^{-17} & -3.9466*10^{-17}\cr
4.38117*10^{-5} & 1.00001 & 0.00508558 & 1.01268*10^{-16} & -4.91152*10^{-16} & -5.12132*10^{-16}\cr
0.0 & 0.00508558 & 5.99999 & 2.76901*10^{-7} & 1.03138*10^{-16} & 1.38417*10^{-16}\cr
0.0 & 0.0 & 2.76901*10^{-7} & 3.32957 & 0.741966 & -1.46329*10^{-16}\cr
0.0 & 0.0 & 0.0 & 0.741966 & 4.67043 & -3.44539*10^{-17}\cr
0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 4.0} \right] \)
特徵值\( [-1.01,1.00001,5.99999,3.32957,4.67043,4.0] \)
執行次數21
(%o8) \( [-1.01,1.00001,5.99999,3.32957,4.67043,4.0] \)

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

TOP

3-3.減少QR法運算-三對角矩陣的Implicit Shift
前一篇文章所述對稱矩陣偏移值的計算可以簡化成\( \displaystyle \alpha=a_n-\frac{b_{n-1}^2}{\delta+sign(\delta)\sqrt{\delta^2+b_{n-1}^2}} \)
再計算\(A-\alpha I\)的\(QR\)分解
計算\(RQ+\alpha I\)得到下一輪的\(A\)矩陣
這是\(Explicit\) \(Single-Shift\) \(QR\)疊代方法

這篇文章介紹\( Implicit \) \(Shift\)疊代方法
以\(n=6\)舉例說明
將對稱矩陣先轉換成三對角矩陣
https://math.pro/db/viewthread.php?tid=2561&page=2#pid16538
\( A=\left[ \matrix{\times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times} \right] \Rightarrow T=\left[ \matrix{
\times & \times & 0 & 0 & 0 & 0 \cr
\times & \times & \times & 0 & 0 & 0 \cr
0 & \times & \times & \times & 0 & 0 \cr
0 & 0 & \times & \times & \times & 0 \cr
0 & 0 & 0 & \times & \times & \times \cr
0 & 0 & 0 & 0 & \times & \times} \right] \)
找一個\(2 \times 2\)旋轉矩陣\(G\),使得\( \left[ \matrix{cos\theta & -sin\theta \cr sin\theta&cos\theta} \right] \left[ \matrix{T[1,1]-\alpha \cr T[2,1]} \right]=\left[ \matrix{* \cr 0} \right] \)
設\(G_1=diag(\left[ \matrix{cos\theta& -sin\theta \cr sin\theta & cos\theta} \right],I_4)\),計算\( T_1=G_1 T G_1^T=\left[ \matrix{
\times & \times & + & 0 & 0 & 0 \cr
\times & \times & \times & 0 & 0 & 0 \cr
+ & \times & \times & \times & 0 & 0 \cr
0 & 0 & \times & \times & \times & 0 \cr
0 & 0 & 0 & \times & \times & \times \cr
0 & 0 & 0 & 0 & \times & \times} \right] \) 其中\( \matrix{&+ \cr +&} \)形成凸起(bulge)
接下來將凸起往右下角推
取\( \left[ \matrix{T_1[2,1]\cr T_1[3,1]} \right]=\left[ \matrix{\times \cr +} \right] \)找一個\(2 \times 2\)旋轉矩陣,使得\( \left[ \matrix{cos \theta & -sin \theta \cr sin \theta & cos \theta} \right] \left[ \matrix{\times \cr +} \right]=\left[ \matrix{* \cr 0} \right] \)
設\(G_2=diag(I_1,\left[ \matrix{cos\theta& -sin\theta \cr sin\theta & cos\theta} \right],I_3)\),計算\( T_2=G_2 T_1 G_2^T=\left[ \matrix{
\times & \times & 0 & 0 & 0 & 0 \cr
\times & \times & \times & + & 0 & 0 \cr
0 & \times & \times & \times & 0 & 0 \cr
0 & + & \times & \times & \times & 0 \cr
0 & 0 & 0 & \times & \times & \times \cr
0 & 0 & 0 & 0 & \times & \times} \right] \) 將凸起\( \matrix{&+ \cr +&} \)往右下角推
取\( \left[ \matrix{T_2[3,2]\cr T_2[4,2]} \right]=\left[ \matrix{\times \cr +} \right] \)找一個\(2 \times 2\)旋轉矩陣,使得\( \left[ \matrix{cos \theta & -sin \theta \cr sin \theta & cos \theta} \right] \left[ \matrix{\times \cr +} \right]=\left[ \matrix{* \cr 0} \right] \)
設\(G_3=diag(I_2,\left[ \matrix{cos\theta& -sin\theta \cr sin\theta & cos\theta} \right],I_2)\),計算\( T_3=G_3 T_2 G_3^T=\left[ \matrix{
\times & \times & 0 & 0 & 0 & 0 \cr
\times & \times & \times & 0 & 0 & 0 \cr
0 & \times & \times & \times & + & 0 \cr
0 & 0 & \times & \times & \times & 0 \cr
0 & 0 & + & \times & \times & \times \cr
0 & 0 & 0 & 0 & \times & \times} \right] \) 將凸起\( \matrix{&+ \cr +&} \)往右下角推
取\( \left[ \matrix{T_3[4,3]\cr T_3[5,3]} \right]=\left[ \matrix{\times \cr +} \right] \)找一個\(2 \times 2\)旋轉矩陣,使得\( \left[ \matrix{cos \theta & -sin \theta \cr sin \theta & cos \theta} \right] \left[ \matrix{\times \cr +} \right]=\left[ \matrix{* \cr 0} \right] \)
設\(G_4=diag(I_3,\left[ \matrix{cos\theta& -sin\theta \cr sin\theta & cos\theta} \right],I_1)\),計算\( T_4=G_4 T_3 G_4^T=\left[ \matrix{
\times & \times & 0 & 0 & 0 & 0 \cr
\times & \times & \times & 0 & 0 & 0 \cr
0 & \times & \times & \times & 0 & 0 \cr
0 & 0 & \times & \times & \times & + \cr
0 & 0 & 0 & \times & \times & \times \cr
0 & 0 & 0 & + & \times & \times} \right] \) 將凸起\( \matrix{&+ \cr +&} \)往右下角推
取\( \left[ \matrix{T_4[5,4]\cr T_4[6,4]} \right]=\left[ \matrix{\times \cr +} \right] \)找一個\(2 \times 2\)旋轉矩陣,使得\( \left[ \matrix{cos \theta & -sin \theta \cr sin \theta & cos \theta} \right] \left[ \matrix{\times \cr +} \right]=\left[ \matrix{* \cr 0} \right] \)
設\(G_5=diag(I_4,\left[ \matrix{cos\theta& -sin\theta \cr sin\theta & cos\theta} \right])\),計算\( T_5=G_5 T_4 G_5^T=\left[ \matrix{
\times & \times & 0 & 0 & 0 & 0 \cr
\times & \times & \times & 0 & 0 & 0 \cr
0 & \times & \times & \times & 0 & 0 \cr
0 & 0 & \times & \times & \times & 0 \cr
0 & 0 & 0 & \times & \times & \times \cr
0 & 0 & 0 & 0 & \times & \times} \right] \) 仍是三對角矩陣


---------------------------
虛擬碼
TridiagonalImplicitShift(T)
{
\( \displaystyle \delta=\frac{T[n-1,n-1]-T[n,n]}{2} \)
偏移值\( \displaystyle \alpha=T[n,n]-\frac{T[n,n-1]^2}{\delta+sign(\delta)\sqrt{\delta^2+T[n,n-1]^2}} \)
\( x=T[1,1]-\alpha \)
\( y=T[2,1] \)
for k=1 ~ (n-1)
  {\( \displaystyle cos \theta=\frac{x}{\sqrt{x^2+y^2}} \)
   \( \displaystyle sin \theta=-\frac{y}{\sqrt{x^2+y^2}} \)
   旋轉矩陣\( G=\left[ \matrix{cos \theta & -sin \theta \cr sin \theta & cos \theta} \right] \)
   \( G_k=diag(I_{k-1},G,I_{n-k-1}) \)
   計算\( T=G_k \cdot T \cdot G_k^T \)
  }
\(return(T)\),處理後仍是三對角矩陣
}
---------------------------
參考資料https://books.google.com.tw/book ... e&q&f=false

註:本文章採用第一種旋轉矩陣
若旋轉矩陣採用\( G=\left[ \matrix{cos \theta & -sin \theta \cr sin \theta & cos \theta} \right] \),則下一輪三對角矩陣計算方式為\(T_k=G_k.T_{k-1}.G_k^T\)
WIKI https://en.wikipedia.org/wiki/Givens_rotation
線代啟示錄 https://ccjou.wordpress.com/2010 ... qr-分解的應用/

別的書採用另一種旋轉矩陣
若旋轉矩陣採用\( G=\left[ \matrix{cos \theta & sin \theta \cr -sin \theta & cos \theta} \right] \),則下一輪三對角矩陣計算方式為\(T_k=G_k^T.T_{k-1}.G_k\)
Matrix Computations https://books.google.com.tw/book ... e&q&f=false
Handbook of Linear Algebra https://books.google.com.tw/book ... e&q&f=false




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

針對三對角矩陣Implicit Shift副程式
(%i2)
TridiagonalImplicitShift(T):=block
([n:length(T),delta,alpha,x,y,r,cos,sin],
print("δ=(T[",n-1,",",n-1,"]-T[",n,",",n,"])/2=(",T[n-1,n-1],"-",T[n,n],")/2=",delta: (T[n-1,n-1]-T[n,n])/2),
if delta=0 then
  (print("偏移值α=T[",n,",",n,"]-|T[",n-1,",",n,"]|=",T[n,n],"-|",T[n-1,n],"|=",alpha:T[n,n]-abs(T[n-1,n])))
else
  (print("偏移值α=T[",n,",",n,"]-T[",n,",",n-1,"]^2/(δ+sign(δ)*sqrt(δ^2+T[",n,",",n-1,"]^2))"),
   print("偏移值α=",T[n,n],"-(",T[n,n-1],")^2/(",delta,"+",signum(delta),"*sqrt((",delta,")^2+(",T[n,n-1],")^2))=",
   alpha:T[n,n]-T[n-1,n]^2/(delta+signum(delta)*sqrt(delta^2+T[n-1,n]^2)))
  ),
print("取",matrix(["T[1,1]-α"],["T[2,1]"]),"=",matrix(["x"],["y"]),"=",matrix([x:T[1,1]-alpha],[y:T[2,1]])),
for k:1 thru n-1 do
   (print("------------------------"),
    print("k=",k),
    print("r=sqrt(x^2+y^2)=",r:sqrt(x^2+y^2)),
    print("cosθ=x/r=",cos:x/r," , sinθ=-y/r=",sin:-y/r),
    G:ident(n),
    G[k,k]:cos,     G[k,k+1]:-sin,  
    G[k+1,k]:sin, G[k+1,k+1]:cos,
    print("取",matrix([x],[y]),"找一個2*2旋轉矩陣,使得",matrix([cos,-sin],[sin,cos]),matrix([x],[y]),"=",matrix([cos,-sin],[sin,cos]).matrix([x],[y])),
    print("設G_",k,"=",G,",計算T_",k,"=G_",k,".T_",k-1,".G_",k,"^T=",T:ratsimp(G.T.transpose(G))),
    if k<n-1 then
      (x:T[k+1,k],
       y:T[k+2,k],
       print("取",matrix([concat("T",k,"[",k+1,",",k,"]")],[concat("T",k,"[",k+2,",",k,"]")]),"=",matrix(["x"],["y"]),"=",matrix([x],[y]))
      )
   ),
return(T)
)$


三對角矩陣
(%i3)
T:matrix([1,4,0,0,0,0],
              [4,25/12,1,0,0,0],
              [0,1,-31/20,1,0,0],
              [0,0,1,1/20,1,0],
              [0,0,0,1,2,2],
              [0,0,0,0,2,-1]);

(%o3) \(\left[ \matrix{\displaystyle 1&4&0&0&0&0 \cr 4&\frac{25}{12}&1&0&0&0 \cr 0&1&-\frac{31}{20}&1&0&0 \cr 0&0&1&\frac{1}{20}&1&0 \cr 0&0&0&1&2&2 \cr 0&0&0&0&2&-1} \right] \)

執行三對角矩陣的Implicit Shift
(%i4) TridiagonalImplicitShift(T);
\( \displaystyle \delta=(T[5,5]-T[6,6])/2=(2--1)/2=\frac{3}{2} \)
偏移值\( \alpha=T[6,6]-T[6,5]^2/(\delta+sign(\delta)*sqrt(\delta^2+T[6,5]^2)) \)
偏移值\( \displaystyle \alpha=-1-(2)^2/(\frac{3}{2}+1*sqrt((\frac{3}{2})^2+(2)^2))=-2 \)
取\( \left[ \matrix{T[1,1]-\alpha \cr T[2,1]} \right]=\left[ \matrix{x \cr y} \right]=\left[ \matrix{3 \cr 4} \right] \)
------------------------
\( k=1 \)
\( r=sqrt(x^2+y^2)=5 \)
\( \displaystyle cos \theta=x/r=\frac{3}{5},sin \theta=-y/r=-\frac{4}{5} \)
取\( \left[ \matrix{3 \cr 4} \right] \)找一個2*2旋轉矩陣,使得\( \left[ \matrix{\displaystyle \frac{3}{5}&\frac{4}{5}\cr -\frac{4}{5}&\frac{3}{5}} \right] \left[ \matrix{3 \cr 4} \right]=\left[ \matrix{5 \cr 0} \right] \)
設\( G_1=\left[ \matrix{\displaystyle \frac{3}{5}&\frac{4}{5}&0&0&0&0 \cr -\frac{4}{5}&\frac{3}{5}&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] \),計算\( T_1=G_1 \cdot T_0 \cdot G_1^T=\left[ \matrix{\displaystyle
\frac{83}{15}&-\frac{3}{5}&\frac{4}{5}&0&0&0 \cr
-\frac{3}{5}&-\frac{49}{20}&\frac{3}{5}&0&0&0 \cr
\frac{4}{5}&\frac{3}{5}&-\frac{31}{20}&1&0&0 \cr
0&0&1&\frac{1}{20}&1&0 \cr
0&0&0&1&2&2 \cr
0&0&0&0&2&-1} \right] \)
取\( \left[ \matrix{T1[2,1] \cr T1[3,1]} \right]=\left[ \matrix{x \cr y} \right]=\left[ \matrix{\displaystyle -\frac{3}{5} \cr \frac{4}{5}} \right] \)
------------------------
\( k=2 \)
\( r=sqrt(x^2+y^2)=1 \)
\( \displaystyle cos \theta=x/r=-\frac{3}{5},sin \theta=-y/r=-\frac{4}{5} \)
取\( \left[ \matrix{\displaystyle -\frac{3}{5}\cr \frac{4}{5}} \right] \)找一個2*2旋轉矩陣,使得\( \left[ \matrix{\displaystyle -\frac{3}{5}&\frac{4}{5}\cr -\frac{4}{5}&-\frac{3}{5}} \right] \left[ \matrix{\displaystyle -\frac{3}{5}\cr \frac{4}{5}} \right]=\left[ \matrix{1 \cr 0}\right] \)
設\( G_2=\left[ \matrix{\displaystyle 1&0&0&0&0&0 \cr 0&-\frac{3}{5}&\frac{4}{5}&0&0&0 \cr 0&-\frac{4}{5}&-\frac{3}{5}&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] \),計算\( T_2=G_2 \cdot T_1 \cdot G_2^T=\left[ \matrix{\displaystyle
\frac{83}{15}&1&0&0&0&0 \cr
1&-\frac{49}{20}&-\frac{3}{5}&\frac{4}{5}&0&0 \cr
0&-\frac{3}{5}&-\frac{31}{20}&-\frac{3}{5}&0&0 \cr
0&\frac{4}{5}&-\frac{3}{5}&\frac{1}{20}&1&0 \cr
0&0&0&1&2&2 \cr
0&0&0&0&2&-1} \right] \)
取\( \left[ \matrix{T2[3,2] \cr T2[4,2]} \right]=\left[ \matrix{x \cr y} \right]=\left[ \matrix{\displaystyle -\frac{3}{5} \cr \frac
{4}{5}} \right] \)
------------------------
\( k=3 \)
\( r=sqrt(x^2+y^2)=1 \)
\( \displaystyle cos \theta=x/r=-\frac{3}{5},sin \theta=-y/r=-\frac{4}{5} \)
取\( \left[ \matrix{\displaystyle -\frac{3}{5} \cr \frac{4}{5}} \right] \)找一個2*2旋轉矩陣,使得\( \left[ \matrix{\displaystyle -\frac{3}{5}&\frac{4}{5}\cr -\frac{4}{5}&-\frac{3}{5}} \right] \left[ \matrix{\displaystyle -\frac{3}{5} \cr \frac{4}{5}} \right]=\left[ \matrix{1 \cr 0} \right] \)
設\( G_3=\left[ \matrix{\displaystyle 1&0&0&0&0&0 \cr 0&1&0&0&0&0 \cr 0&0&-\frac{3}{5}&\frac{4}{5}&0&0 \cr 0&0&-\frac{4}{5}&-\frac{3}{5}&0&0 \cr 0&0&0&0&1&0 \cr 0&0&0&0&0&1} \right] \),計算\( T_3=G_3 \cdot T_2 \cdot G_3^T=\left[ \matrix{\displaystyle
\frac{83}{15}&1&0&0&0&0 \cr
1&-\frac{49}{20}&1&0&0&0 \cr
0&1&\frac{1}{20}&-\frac{3}{5}&\frac{4}{5}&0 \cr
0&0&-\frac{3}{5}&-\frac{31}{20}&-\frac{3}{5}&0 \cr
0&0&\frac{4}{5}&-\frac{3}{5}&2&2 \cr
0&0&0&0&2&-1} \right] \)
取\( \left[ \matrix{T3[4,3] \cr T3[5,3]} \right]=\left[ \matrix{x \cr y} \right]=\left[ \matrix{\displaystyle -\frac{3}{5} \cr \frac{4}{5}} \right] \)
------------------------
\( k=4 \)
\( r=sqrt(x^2+y^2)=1 \)
\( \displaystyle cos \theta=x/r=-\frac{3}{5},sin \theta=-y/r=-\frac{4}{5} \)
取\( \left[ \matrix{\displaystyle -\frac{3}{5}\cr \frac{4}{5}} \right] \)找一個2*2旋轉矩陣,使得\( \left[ \matrix{\displaystyle -\frac{3}{5}&\frac{4}{5}\cr -\frac{4}{5}&-\frac{3}{5}} \right] \left[ \matrix{\displaystyle -\frac{3}{5} \cr \frac{4}{5}} \right]=\left[ \matrix{1 \cr 0} \right] \)
設\( G_4=\left[ \matrix{\displaystyle 1&0&0&0&0&0 \cr 0&1&0&0&0&0 \cr 0&0&1&0&0&0 \cr 0&0&0&-\frac{3}{5}&\frac{4}{5}&0 \cr 0&0&0&-\frac{4}{5}&-\frac{3}{5}&0 \cr 0&0&0&0&0&1} \right] \),計算\( T_4=G_4 \cdot T_3 \cdot G_4^T=\left[ \matrix{\displaystyle
\frac{83}{15}&1&0&0&0&0 \cr
1&-\frac{49}{20}&1&0&0&0 \cr
0&1&\frac{1}{20}&1&0&0 \cr
0&0&1&\frac{649}{500}&-\frac{192}{125}&\frac{8}{5} \cr
0&0&0&-\frac{192}{125}&-\frac{106}{125}&-\frac{6}{5} \cr
0&0&0&\frac{8}{5}&-\frac{6}{5}&-1} \right] \)
取\( \left[ \matrix{T4[5,4] \cr T4[6,4]} \right]=\left[ \matrix{x \cr y} \right]=\left[ \matrix{\displaystyle -\frac{192}{125} \cr \frac{8}{5}} \right] \)
------------------------
\( k=5 \)
\( \displaystyle r=sqrt(x^2+y^2)=\frac{8\sqrt{1201}}{125} \)
\( \displaystyle cos \theta=x/r=-\frac{24}{\sqrt{1201}},sin \theta=-y/r=-\frac{25}{\sqrt{1201}} \)
取\( \left[ \matrix{\displaystyle -\frac{192}{125} \cr \frac{8}{5}} \right] \)找一個2*2旋轉矩陣,使得\(\left[ \matrix{\displaystyle -\frac{24}{\sqrt{1201}}&\frac{25}{\sqrt{1201}}\cr -\frac{25}{\sqrt{1201}}&-\frac{24}{\sqrt{1201}}} \right]
\left[ \matrix{\displaystyle -\frac{192}{125} \cr \frac{8}{5}} \right]=\left[ \matrix{\displaystyle \frac{8\sqrt{1201}}{125} \cr 0} \right] \)
設\( G_5=\left[ \matrix{\displaystyle 1&0&0&0&0&0 \cr 0&1&0&0&0&0 \cr 0&0&1&0&0&0 \cr 0&0&0&0&-\frac{24}{\sqrt{1201}}&\frac{25}{\sqrt{1201}} \cr 0&0&0&0&-\frac{25}{\sqrt{1201}}&-\frac{24}{\sqrt{1201}}} \right] \),計算\( T_5=G_5 \cdot T_4 \cdot G_5^T=\left[ \matrix{\displaystyle
\frac{83}{15}&1&0&0&0&0 \cr
1&-\frac{49}{20}&1&0&0&0 \cr
0&1&\frac{1}{20}&1&0&0 \cr
0&0&1&\frac{649}{500}&\frac{8\sqrt{1201}}{125}&0 \cr
0&0&0&\frac{8\sqrt{1201}}{125}&\frac{40819}{150125}&\frac{150}{1201} \cr
0&0&0&0&\frac{150}{1201}&-\frac{2546}{1201}} \right] \)
(%o4) \( \left[ \matrix{\displaystyle
\frac{83}{15}&1&0&0&0&0 \cr
1&-\frac{49}{20}&1&0&0&0 \cr
0&1&\frac{1}{20}&1&0&0 \cr
0&0&1&\frac{649}{500}&\frac{8\sqrt{1201}}{125}&0 \cr
0&0&0&\frac{8\sqrt{1201}}{125}&\frac{40819}{150125}&\frac{150}{1201} \cr
0&0&0&0&\frac{150}{1201}&-\frac{2546}{1201}} \right] \)

[ 本帖最後由 bugmens 於 2017-11-4 09:55 編輯 ]

TOP

3-4.減少QR法運算-Hessenberg矩陣的Implicit Double Shift
前面文章提到https://math.pro/db/viewthread.php?tid=2561&page=2#pid17755

取矩陣\( 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 \)特徵值為複數

(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\)矩陣

但\(A^2\)矩陣相乘是\( O(n^3) \)的運算,可以先將\(A\)矩陣化簡成\(Hessenberg\)矩陣再進行\( Implicit\) \(Double\) \(Shift \),只需要\(O(n^2)\)的運算。
以\(n=6\)舉例說明
將\(A\)矩陣轉換成\(Hessenberg\)矩陣
\( A=\left[ \matrix{\times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times} \right] \Rightarrow H=\left[ \matrix{\times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times \cr
0 & \times & \times & \times & \times & \times \cr
0 & 0 & \times & \times & \times & \times \cr
0 & 0 & 0 & \times & \times & \times \cr
0 & 0 & 0 & 0 & \times & \times} \right] \)
計算\( H^2-2Re(\sigma)H+\left| \sigma \right|^2 I \)的第一列\( n_1=\left[ \matrix{x \cr y \cr z \cr 0 \cr 0 \cr 0} \right] \),前三個元素非零,其餘元素為0
取向量\( \left[ \matrix{x \cr y \cr z} \right] \)找一個\( 3\times 3 Householder \)矩陣\( \hat{P_0} \),使得\( \hat{P_0}\left[ \matrix{x \cr y \cr z} \right]=\left[ \matrix{* \cr 0 \cr 0} \right] \)
設\( P_0=diag(\hat{P_0},I_3) \),計算\( H_0=P_0^T H P_0 =\left[ \matrix{\times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times \cr
+ & \times & \times & \times & \times & \times \cr
+ & + & \times & \times & \times & \times \cr
0 & 0 & 0 & \times & \times & \times \cr
0 & 0 & 0 & 0 & \times & \times} \right] \) 其中\( \matrix{+& \cr +& +} \)形成凸起(bulge)

接下來將凸起往右下角推
取向量\( \left[ \matrix{H_0[2,1] \cr H_0[3,1] \cr H_0[4,1]} \right]=\left[ \matrix{\times \cr + \cr +} \right] \)找一個\( 3\times 3 Householder \)矩陣\( \hat{P_1} \),使得\( \hat{P_1}\left[ \matrix{\times \cr + \cr +} \right]=\left[ \matrix{* \cr 0 \cr 0} \right] \)
設\( P_1=diag(I_1,\hat{P_1},I_2) \),計算\( H_1=P_1^T H_0 P_1 =\left[ \matrix{\times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times \cr
0 & \times & \times & \times & \times & \times \cr
0 & + & \times & \times & \times & \times \cr
0 & + & + & \times & \times & \times \cr
0 & 0 & 0 & 0 & \times & \times} \right] \) 將凸起\( \matrix{+& \cr +& +} \)往右下角推

取向量\( \left[ \matrix{H_1[3,2] \cr H_1[4,2] \cr H_1[5,2]} \right]=\left[ \matrix{\times \cr + \cr +} \right] \)找一個\( 3\times 3 Householder \)矩陣\( \hat{P_2} \),使得\( \hat{P_2}\left[ \matrix{\times \cr + \cr +} \right]=\left[ \matrix{* \cr 0 \cr 0} \right] \)
設\( P_2=diag(I_2,\hat{P_2},I_1) \),計算\( H_2=P_2^T H_1 P_2 =\left[ \matrix{\times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times \cr
0 & \times & \times & \times & \times & \times \cr
0 & 0 & \times & \times & \times & \times \cr
0 & 0 & + & \times & \times & \times \cr
0 & 0 & + & + & \times & \times} \right] \) 將凸起\( \matrix{+& \cr +& +} \)往右下角推

取向量\( \left[ \matrix{H_2[4,3] \cr H_2[5,3] \cr H_2[6,3]} \right]=\left[ \matrix{\times \cr + \cr +} \right] \)找一個\( 3\times 3 Householder \)矩陣\( \hat{P_3} \),使得\( \hat{P_3}\left[ \matrix{\times \cr + \cr +} \right]=\left[ \matrix{* \cr 0 \cr 0} \right] \)
設\( P_3=diag(I_3,\hat{P_3}) \),計算\( H_3=P_3^T H_2 P_3 =\left[ \matrix{\times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times \cr
0 & \times & \times & \times & \times & \times \cr
0 & 0 & \times & \times & \times & \times \cr
0 & 0 & 0 & \times & \times & \times \cr
0 & 0 & 0 & + & \times & \times} \right] \) 將凸起+往右下角推

取向量\( \left[ \matrix{H_3[5,4] \cr H_3[6,4]} \right]=\left[ \matrix{\times \cr +} \right] \)找一個\( 2\times 2 Householder \)矩陣\( \hat{P_4} \),使得\( \hat{P_4}\left[ \matrix{\times \cr +} \right]=\left[ \matrix{* \cr 0} \right] \)
設\( P_4=diag(I_4,\hat{P_4}) \),計算\( H_4=P_4^T H_3 P_4=\left[ \matrix{\times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times \cr
0 & \times & \times & \times & \times & \times \cr
0 & 0 & \times & \times & \times & \times \cr
0 & 0 & 0 & \times & \times & \times \cr
0 & 0 & 0 & 0 & \times & \times} \right] \) 仍是\(Hessenberg\)矩陣


---------------------------
虛擬碼
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) \)
}

HessenbergImplicitDoubleShift(H)
{
\( t=H[n-1,n-1]+H[n,n] \)
\( d=H[n-1,n-1]H[n,n]-H[n,n-1]H[n-1,n] \)
\( \cases{x=H[1,1]^2-tH[1,1]+d+H[1,2]H[2,1] \cr y=H[2,1](H[1,1]+H[2,2]-t) \cr z=H[2,1]H[3,2]} \),計算\( H^2-tH+dI \)的第一列的前三元\( \left[ \matrix{x \cr y \cr z} \right] \)
for k=1 ~ (n-3)
 {\( \hat{P_k}=Householder(\left[ \matrix{x \cr y \cr z} \right]) \),找一個\( 3 \times 3 Householder \)矩陣\( \hat{P_k} \),使得\( \hat{P_k}\left[ \matrix{x \cr y \cr z} \right]=\left[ \matrix{* \cr 0 \cr 0} \right] \),*為非零數字
  \( P_k=diag(I_k,\hat{P_k},I_{n-k-3}) \)
  \( H=P_k^T H P_k \)
  \( x=H[k+2,k+1] \)
  \( y=H[k+3,k+1] \)
  if \(k<n-3\) then
  \( \{\; z=H[k+4,k+1] \}\; \)
 }
\( \hat{P}_{n-2}=Householder(\left[ \matrix{x \cr y} \right]) \),找一個\( 2 \times 2 Householder \)矩陣\( \hat{P}_{n-2} \),使得\( \hat{P}_{n-2}\left[ \matrix{x \cr y} \right]=\left[ \matrix{* \cr 0} \right] \),*為非零數字
\( P_{n-2}=diag(I_{n-2},\hat{P}_{n-2}) \)
\( H=P_{n-2}^T H P_{n-2} \)
\( return(H) \),處理後仍是\( Hessenberg \)矩陣
}
---------------------------
參考資料https://books.google.com.tw/book ... e&q&f=false



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

轉換成Householder矩陣副程式
(%i2)
Householder(X):=block
([n,e1,normx,v,normv,P_hat],
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向量=",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("P_hat=I-2*v.v^T=",ident(n),"-2*",v,".",transpose(v),"=",P_hat:ratsimp(ident(n)-2*v.transpose(v))),
return(P_hat)
)$


針對Hessenberg矩陣Implicit Double Shift副程式
(%i3)
HessenbergImplicitDoubleShift(H):=block
([n,t,d,x,y,z,k,X,P_hat_k],
n:length(H),
print("t=H[",n-1,",",n-1,"]+H[",n,",",n,"]=",H[n-1,n-1],"+",H[n,n],"=",t:H[n-1,n-1]+H[n,n]),
print("d=H[",n-1,",",n-1,"]*H[",n,",",n,"]-H[",n,",",n-1,"]*H[",n-1,",",n,"]=",
         H[n-1,n-1],"*",H[n,n],"-",H[n,n-1],"*",H[n-1,n],"=",d:H[n-1,n-1]*H[n,n]-H[n,n-1]*H[n-1,n]),
print("x=H[1,1]^2-t*H[1,1]+d+H[1,2]*H[2,1]=("
            ,H[1,1],")^2-",t,"*",H[1,1],"+",d,"+",H[1,2],"*",H[2,1],"=",x:H[1,1]^2-t*H[1,1]+d+H[1,2]*H[2,1]),
print("y=H[2,1]*(H[1,1]+H[2,2]-t)=",H[2,1],"*(",H[1,1],"+",H[2,2],"-",t,")=",y:H[2,1]*(H[1,1]+H[2,2]-t)),
print("z=H[2,1]*H[3,2]=",H[2,1],"*",H[3,2],"=",z:H[2,1]*H[3,2]),
print("取向量",matrix(["x"],["y"],["z"]),"=",matrix([x],[y],[z])),
print("------------------------"),
for k:0 thru n-3 do
  (print("k=",k),
   print("取向量",X:matrix([x],[y],[z]),"找一個3*3Householder矩陣P_hat_",k,",使得",P_hat_k:Householder(X),X,"=",P_hat_k.X),
   if k=0 then
      (P_k:diag([P_hat_k,ident(n-k-3)]))
    else if k=n-3 then
      (P_k:diag([ident(n-3),P_hat_k]))
    else
      (P_k:diag([ident(k),P_hat_k,ident(n-k-3)])),
   print("設P_",k,"=",P_k,",計算H_",k,"=P_",k,"^T.H_",k-1,".P_",k,"=",H:ratsimp(transpose(P_k).H.P_k)),
   x:H[k+2,k+1],
   y:H[k+3,k+1],
   if k<n-3 then
     (z:H[k+4,k+1],
      print("取向量",matrix([concat("H",k,"[",k+2,",",k+1,"]")],[concat("H",k,"[",k+3,",",k+1,"]")],[concat("H",k,"[",k+4,",",k+1,"]")]),"=",matrix([x],[y],[z]))
     )
   else
     (print("取向量",matrix([concat("H",k,"[",k+2,",",k+1,"]")],[concat("H",k,"[",k+3,",",k+1,"]")]),"=",matrix([x],[y]))
     ),
   print("------------------------")
  ),
  k:n-2,
  print("取向量",X:matrix([x],[y]),"找一個2*2Householder矩陣P_hat_",k,",使得",P_hat_k:Householder(X),X,"=",P_hat_k.X),
  print("設P_",k,"=",P_k:diag([ident(n-2),P_hat_k]),",計算H",k,"=P_",k,"^T.H_",k-1,".P_",k,"=",H:ratsimp(transpose(P_k).H.P_k))
)$


Hessenberg矩陣,各元數字經過設計,\( \left[ \matrix{x \cr y \cr z} \right] \)向量長度均為完全平方數,處理後的\( Hessenberg \)矩陣各元最多是分數,不會出現根號數字
(%i4)
H:matrix([-1,-1,2/5,-1/4,11/8,31/80],
             [1,-2,2,11/5,-11/5,-21/5],
             [0,-2,1,2,2,-10],
             [0,0,3,-3,0,10],
             [0,0,0,3,-2,-2],
             [0,0,0,0,1,-2]);

(%o4) \( \left[ \matrix{\displaystyle
-1& -1 & \frac{2}{5} & -\frac{1}{4} & \frac{11}{8} & \frac{31}{80} \cr
1 & -2 & 2 & \frac{11}{5} & -\frac{11}{5} & -\frac{21}{5} \cr
0& -2 & 1 & 2 & 2 & -10 \cr
0& 0 & 3 & -3 & 0 & 10 \cr
0& 0 & 0 & 3 & -2 & -2 \cr
0& 0 & 0 & 0 & 1 & -2} \right] \)

執行\(Hessenberg\)矩陣的\(Implicit\) \(Double\) \(Shift\)
(%i5) HessenbergImplicitDoubleShift(H);
\( t=H[5,5]+H[6,6]=-2+-2=-4 \)
\( d=H[5,5]*H[6,6]-H[6,5]*H[5,6]=-2*-2-1*-2=6 \)
\( x=H[1,1]^2-t*H[1,1]+d+H[1,2]*H[2,1]=(-1)^2--4*-1+6+-1*1=2 \)
\( y=H[2,1]*(H[1,1]+H[2,2]-t)=1*(-1+-2--4 )=1 \)
\( z=H[2,1]*H[3,2]=1*-2=-2 \)
取向量\( \left[ \matrix{x \cr y \cr z} \right]=\left[ \matrix{2 \cr 1 \cr -2} \right] \)
------------------------
\( k=0 \)
\( X \)向量\( =\left[ \matrix{2 \cr 1 \cr -2} \right] \)第1元\(=2\),\(sign(X_1)=1\),長度\( \Vert\; X \Vert\;=sqrt(\left[ \matrix{2 \cr 1 \cr -2} \right] \cdot \left[ \matrix{2 \cr 1 \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{2 \cr 1 \cr -2} \right]+1*3*\left[ \matrix{1 \cr 0 \cr 0} \right]=\left[ \matrix{5 \cr 1 \cr -2} \right] \)
長度\( \Vert\; X+sign(X_1)*\Vert\; X \Vert\;*e1 \Vert\;=sqrt(\left[ \matrix{5 \cr 1 \cr -2} \right] \cdot \left[ \matrix{5 \cr 1 \cr -2} \right])=\sqrt{30} \)
\( v=\left[ \matrix{5 \cr 1 \cr -2} \right]/ \sqrt{30}=\left[ \matrix{\displaystyle \frac{5}{\sqrt{30}}\cr \frac{1}{\sqrt{30}}\cr -\frac{2}{\sqrt{30}}} \right] \)
\( \hat{P}=I-2*v.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{1}{\sqrt{30}} \cr -\frac{2}{\sqrt{30}}} \right] \cdot \left[ \matrix{\displaystyle \frac{5}{\sqrt{30}}& \frac{1}{\sqrt{30}}& -\frac{2}{\sqrt{30}}} \right]=\left[ \matrix{\displaystyle -\frac{2}{3}&-\frac{1}{3}&\frac{2}{3} \cr
-\frac{1}{3}&\frac{14}{15}&\frac{2}{15} \cr
\frac{2}{3}&\frac{2}{15}&\frac{11}{15}} \right]\)
取向量\( \left[ \matrix{2 \cr 1 \cr -2} \right] \)找一個3*3\( Householder \)矩陣\( \hat{P_0} \),使得\( \left[ \matrix{\displaystyle -\frac{2}{3}&-\frac{1}{3}&\frac{2}{3} \cr
-\frac{1}{3}&\frac{14}{15}&\frac{2}{15} \cr
\frac{2}{3}&\frac{2}{15}&\frac{11}{15}} \right] \left[ \matrix{2 \cr 1 \cr -2} \right]=\left[ \matrix{-3 \cr 0 \cr 0} \right] \)
設\( P_0=\left[ \matrix{\displaystyle -\frac{2}{3}&-\frac{1}{3}&\frac{2}{3}&0&0&0 \cr
-\frac{1}{3}&\frac{14}{15}&\frac{2}{15}&0&0&0 \cr
\frac{2}{3}&\frac{2}{15}&\frac{11}{15}&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] \),計算\(H_0=P_0^T.H_{-1}.P_0=\left[ \matrix{\displaystyle
-\frac{2}{5}&-\frac{11}{75}&\frac{2}{75}&\frac{23}{30}&\frac{23}{20}&-\frac{221}{40} \cr
1&-\frac{139}{75}&\frac{148}{75}&\frac{721}{300}&-\frac{449}{200}&-\frac{2153}{400} \cr
2&-\frac{142}{75}&\frac{19}{75}&\frac{239}{150}&\frac{209}{100}&-\frac{1527}{200} \cr
2&\frac{2}{5}&\frac{11}{5}&-3&0&10 \cr
0&0&0&3&-2&-2 \cr
0&0&0&0&1&-2 \cr} \right]\)
取向量\( \left[ \matrix{H0[2,1]\cr H0[3,1]\cr H0[4,1]} \right]=\left[ \matrix{1 \cr 2 \cr 2} \right] \)
------------------------
\(k=1\)
\( X \)向量\( =\left[ \matrix{1 \cr 2 \cr 2} \right] \)第1元\(=1\),\(sign(X_1)=1\),長度\( \Vert\; X \Vert\;=sqrt(\left[ \matrix{1 \cr 2 \cr 2} \right] \cdot \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] \cdot \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] \)
\( \hat{P}=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]\)
取向量\( \left[ \matrix{1 \cr 2 \cr 2} \right] \)找一個3*3\( Householder \)矩陣\( \hat{P_1} \),使得\( \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] \left[ \matrix{1 \cr 2 \cr 2} \right]=\left[ \matrix{-3 \cr 0 \cr 0} \right] \)
設\( P_1=\left[ \matrix{\displaystyle 1&0&0&0&0&0 \cr
0&-\frac{1}{3}&-\frac{2}{3}&-\frac{2}{3}&0&0 \cr
0&-\frac{2}{3}&\frac{2}{3}&-\frac{1}{3}&0&0 \cr
0&-\frac{2}{3}&-\frac{1}{3}&\frac{2}{3}&0&0 \cr
0&0&0&0&1&0 \cr
0&0&0&0&0&1 } \right] \),計算\(H_1=P_1^T.H_0.P_1=\left[ \matrix{\displaystyle
-\frac{2}{5}&-\frac{12}{25}&-\frac{7}{50}&\frac{3}{5}&\frac{23}{20}&-\frac{221}{40} \cr
-3&\frac{9}{10}&-\frac{53}{20}&-\frac{11}{50}&-\frac{129}{200}&\frac{87}{400} \cr
0&1&-\frac{13}{10}&\frac{26}{25}&\frac{289}{100}&-\frac{967}{200} \cr
0&2&0&-\frac{21}{5}&\frac{4}{5}&\frac{64}{5} \cr
0&-2&-1&2&-2&-2 \cr
0&0&0&0&1&-2 \cr} \right]\)
取向量\( \left[ \matrix{H1[3,2]\cr H1[4,2]\cr H1[5,2]} \right]=\left[ \matrix{1 \cr 2 \cr -2} \right] \)
------------------------
\(k=2\)
\( X \)向量\( =\left[ \matrix{1 \cr 2 \cr -2} \right] \)第1元\(=1\),\(sign(X_1)=1\),長度\( \Vert\; X \Vert\;=sqrt(\left[ \matrix{1 \cr 2 \cr -2} \right] \cdot \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] \cdot \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] \)
\( \hat{P}=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]\)
取向量\( \left[ \matrix{1 \cr 2 \cr -2} \right] \)找一個3*3\( Householder \)矩陣\( \hat{P_2} \),使得\( \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] \left[ \matrix{1 \cr 2 \cr -2} \right]=\left[ \matrix{-3 \cr 0 \cr 0} \right] \)
設\( P_2=\left[ \matrix{\displaystyle 1&0&0&0&0&0 \cr
0&1&0&0&0&0 \cr
0&0&-\frac{1}{3}&-\frac{2}{3}&\frac{2}{3}&0 \cr
0&0&-\frac{2}{3}&\frac{2}{3}&\frac{1}{3}&0 \cr
0&0&\frac{2}{3}&\frac{1}{3}&\frac{2}{3}&0 \cr
0&0&0&0&0&1 } \right] \),計算\(H_2=P_2^T.H_1.P_2=\left[ \matrix{\displaystyle
-\frac{2}{5}&-\frac{12}{25}&\frac{31}{75}&\frac{263}{300}&\frac{131}{150}&-\frac{221}{40} \cr
-3&\frac{9}{10}&\frac{3}{5}&\frac{281}{200}&-\frac{227}{100}&\frac{87}{400} \cr
0&-3&-\frac{13}{3}&\frac{521}{300}&-\frac{39}{50}&-\frac{1651}{200} \cr
0&0&\frac{1}{3}&-\frac{439}{150}&-\frac{49}{25}&\frac{1109}{100} \cr
0&0&\frac{2}{3}&\frac{259}{150}&-\frac{6}{25}&-\frac{29}{100} \cr
0&0&\frac{2}{3}&\frac{1}{3}&\frac{2}{3}&-2 \cr} \right]\)
取向量\( \left[ \matrix{H2[4,3]\cr H2[5,3]\cr H2[6,3]} \right]=\left[ \matrix{\displaystyle \frac{1}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right] \)
------------------------
\(k=3\)
\( X \)向量\( =\left[ \matrix{\displaystyle \frac{1}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right] \)第1元\(=\displaystyle \frac{1}{3}\),\(sign(X_1)=1\),長度\( \Vert\; X \Vert\;=sqrt(\left[ \matrix{\displaystyle \frac{1}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right] \cdot \left[ \matrix{\displaystyle \frac{1}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right])=1 \)
第1元是1,其餘為0的向量\( e1=\left[ \matrix{1 \cr 0 \cr 0} \right] \)
\( X+sign(X_1)*\Vert\; X \Vert\;*e1=\left[ \matrix{\displaystyle \frac{1}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right]+1*1*\left[ \matrix{1 \cr 0 \cr 0} \right]=\left[ \matrix{\displaystyle \frac{4}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right] \)
長度\( \Vert\; X+sign(X_1)*\Vert\; X \Vert\;*e1 \Vert\;=sqrt(\left[ \matrix{\displaystyle \frac{4}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right] \cdot \left[ \matrix{\displaystyle \frac{4}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right])=\frac{2^{3/2}}{\sqrt{3}} \)
\( v=\left[ \matrix{\displaystyle \frac{4}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right]/ \frac{2^{3/2}}{\sqrt{3}}=\left[ \matrix{\displaystyle \frac{\sqrt{2}}{\sqrt{3}}\cr \frac{1}{\sqrt{2}\sqrt{3}}\cr \frac{1}{\sqrt{2}\sqrt{3}}} \right] \)
\( \hat{P}=I-2*v.v^T=\left[ \matrix{1&0&0 \cr 0&1&0 \cr 0&0&1} \right]-2*\left[ \matrix{\displaystyle \frac{\sqrt{2}}{\sqrt{3}} \cr \frac{1}{\sqrt{2}\sqrt{3}} \cr \frac{1}{\sqrt{2}\sqrt{3}}} \right] \cdot \left[ \matrix{\displaystyle \frac{\sqrt{2}}{\sqrt{3}}& \frac{1}{\sqrt{2}\sqrt{3}}& \frac{1}{\sqrt{2}\sqrt{3}}} \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]\)
取向量\( \left[ \matrix{\displaystyle \frac{1}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right] \)找一個3*3\( Householder \)矩陣\( \hat{P_3} \),使得\( \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] \left[ \matrix{\displaystyle \frac{1}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right]=\left[ \matrix{-1 \cr 0 \cr 0} \right] \)
設\( P_3=\left[ \matrix{\displaystyle 1&0&0&0&0&0 \cr
0&1&0&0&0&0 \cr
0&0&1&0&0&0 \cr
0&0&0&-\frac{1}{3}&-\frac{2}{3}&-\frac{2}{3} \cr
0&0&0&-\frac{2}{3}&\frac{2}{3}&-\frac{1}{3} \cr
0&0&0&-\frac{2}{3}&-\frac{1}{3}&\frac{2}{3} } \right] \),計算\(H_3=P_3^T.H_2.P_3=\left[ \matrix{\displaystyle
-\frac{2}{5}&-\frac{12}{25}&\frac{31}{75}&\frac{632}{225}&\frac{3311}{1800}&-\frac{4103}{900} \cr
-3&\frac{9}{10}&\frac{3}{5}&\frac{9}{10}&-\frac{1009}{400}&-\frac{7}{200} \cr
0&-3&-\frac{13}{3}&\frac{49}{9}&\frac{1933}{1800}&-\frac{5761}{900} \cr
0&0&-1&\frac{4}{3}&\frac{1111}{900}&-\frac{587}{450} \cr
0&0&0&3&\frac{209}{225}&-\frac{1556}{225} \cr
0&0&0&4&\frac{2729}{900}&-\frac{3343}{450} \cr} \right]\)
取向量\( \left[ \matrix{H3[5,4]\cr H3[6,4]} \right]=\left[ \matrix{\displaystyle 3 \cr 4} \right] \)
------------------------
\( X \)向量\( =\left[ \matrix{3 \cr 4} \right] \)第1元\(=3\),\(sign(X_1)=1\),長度\( \Vert\; X \Vert\;=sqrt(\left[ \matrix{3 \cr 4} \right] \cdot \left[ \matrix{3 \cr 4} \right])=5 \)
第1元是1,其餘為0的向量\( e1=\left[ \matrix{1 \cr 0} \right] \)
\( X+sign(X_1)*\Vert\; X \Vert\;*e1=\left[ \matrix{3 \cr 4} \right]+1*5*\left[ \matrix{1 \cr 0} \right]=\left[ \matrix{8 \cr 4 } \right] \)
長度\( \Vert\; X+sign(X_1)*\Vert\; X \Vert\;*e1 \Vert\;=sqrt(\left[ \matrix{8 \cr 4} \right] \cdot \left[ \matrix{8 \cr 4} \right])=4\sqrt{5} \)
\( v=\left[ \matrix{8 \cr 4} \right]/ 4\sqrt{5}=\left[ \matrix{\displaystyle \frac{2}{\sqrt{5}}\cr \frac{1}{\sqrt{5}}} \right] \)
\( \hat{P}=I-2*v.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]\)
取向量\( \left[ \matrix{3 \cr 4} \right] \)找一個2*2\( Householder \)矩陣\( \hat{P_4} \),使得\( \left[ \matrix{\displaystyle -\frac{3}{5}&-\frac{4}{5} \cr
-\frac{4}{5}&\frac{3}{5}} \right] \left[ \matrix{3 \cr 4} \right]=\left[ \matrix{-5 \cr 0 \cr} \right] \)
設\( P_4=\left[ \matrix{\displaystyle 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&-\frac{3}{5}&-\frac{4}{5} \cr
0&0&0&0&-\frac{4}{5}&\frac{3}{5} } \right] \),計算\(H_4=P_4^T.H_3.P_4=\left[ \matrix{\displaystyle
-\frac{2}{5}&-\frac{12}{25}&\frac{31}{75}&\frac{632}{225}&\frac{22891}{9000}&-\frac{18931}{4500} \cr
-3&\frac{9}{10}&\frac{3}{5}&\frac{9}{10}&\frac{3083}{2000}&\frac{1997}{1000} \cr
0&-3&-\frac{13}{3}&\frac{49}{9}&\frac{40289}{9000}&-\frac{21149}{4500} \cr
0&0&-1&\frac{4}{3}&\frac{1363}{4500}&-\frac{3983}{2250} \cr
0&0&0&-5&-\frac{35348}{5625}&\frac{47486}{5625} \cr
0&0&0&0&-\frac{33881}{22500}&-\frac{2429}{11250} \cr} \right]\)
(%o5) \( \left[ \matrix{\displaystyle
-\frac{2}{5}&-\frac{12}{25}&\frac{31}{75}&\frac{632}{225}&\frac{22891}{9000}&-\frac{18931}{4500} \cr
-3&\frac{9}{10}&\frac{3}{5}&\frac{9}{10}&\frac{3083}{2000}&\frac{1997}{1000} \cr
0&-3&-\frac{13}{3}&\frac{49}{9}&\frac{40289}{9000}&-\frac{21149}{4500} \cr
0&0&-1&\frac{4}{3}&\frac{1363}{4500}&-\frac{3983}{2250} \cr
0&0&0&-5&-\frac{35348}{5625}&\frac{47486}{5625} \cr
0&0&0&0&-\frac{33881}{22500}&-\frac{2429}{11250} \cr} \right]\)

[ 本帖最後由 bugmens 於 2017-11-4 09:55 編輯 ]

TOP

 23 123
發新話題