Math Pro 數學補給站 » 數學軟體 » 用maxima學數值分析-特徵值和特徵向量

bugmens 發表於 2017-7-19 22:31

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

參考資料
[url]https://books.google.com.tw/books?id=QX7HBAAAQBAJ&lpg=PA447&ots=PR8DfwRHZk&dq=%22an%20even%20better%20approximation%20the%20eigenvalue%22&hl=zh-TW&pg=PA447#v=onepage&q&f=false[/url]


以簡單的整數為例子
矩陣\(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)))
  ),



[color=green]要先載入lapack才能使用dgeqrf指令[/color]
[color=red](%i1)[/color]　[color=blue]load(lapack);[/color]
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
[color=red](%o1)[/color]　[i]C:\maxima-5.40.0\share\maxima\5.40.0\share\lapack\lapack.mac[/i]

[color=green]設定小數點底下第6位四捨五入[/color]
[color=red](%i2)[/color]　[color=blue]fpprintprec:6;[/color]
[color=red](%o2)[/color]　6

[color=green]對稱矩陣的WilkinsonShift[/color]
[color=red](%i3)[/color]
[color=blue]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)
)$[/color]

[color=green]具有WilkinsonShift的QR法副程式[/color]
[color=red](%i4)[/color]
[color=blue]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("－－－－－－－－－－－")
  )
)$[/color]

[color=green]對稱矩陣[/color]
[color=red](%i5)[/color]
[color=blue]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]);[/color]
[color=red](%o5)[/color]　\( \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] \)

[color=green]計算矩陣A全部的特徵值,誤差10^-5[/color]
[color=red](%i6)[/color]　[color=blue]QRwithWilkinsonShift(A,10^-5);[/color]
第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
[color=red](%o6)[/color]　\( [-1.01,6.0,5.0,3.0,0.999997,4.0] \)

[color=green]三對角矩陣[/color]
[color=red](%i7)[/color]
[color=blue]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]);[/color]
[color=red](%o7)[/color]　\( \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] \)

[color=green]計算矩陣A全部的特徵值,誤差10^-5[/color]
[color=red](%i8)[/color]　[color=blue]QRwithWilkinsonShift(A,10^-5);[/color]
第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
[color=red](%o8)[/color]　\( [-1.01,1.00001,5.99999,3.32957,4.67043,4.0] \)

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

bugmens 發表於 2017-10-22 22:03

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\)舉例說明
將對稱矩陣先轉換成三對角矩陣
[url]https://math.pro/db/viewthread.php?tid=2561&page=2#pid16538[/url]
\( 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)\)，處理後仍是三對角矩陣
}
－－－－－－－－－－－－－－－－－－－－－－－－－－－
參考資料[url]https://books.google.com.tw/books?id=X5YfsuCWpxMC&lpg=PA462&dq=%22WE%20ARE%20THUS%20IN%20A%20POSITION%20TO%20APPLY%20THE%20IMPLICIT%22&hl=zh-TW&pg=PA463#v=onepage&q&f=false[/url]

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

別的書採用另一種旋轉矩陣
若旋轉矩陣採用\( 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 [url]https://books.google.com.tw/books?id=X5YfsuCWpxMC&lpg=PA462&dq=%22WE%20ARE%20THUS%20IN%20A%20POSITION%20TO%20APPLY%20THE%20IMPLICIT%22&hl=zh-TW&pg=PA463#v=onepage&q&f=false[/url]
Handbook of Linear Algebra [url]https://books.google.com.tw/books?id=kATOBQAAQBAJ&lpg=SA42-PA1&ots=GAktZeCEnB&dq=Tridiagonal%20Implicit%20Shift&hl=zh-TW&pg=SA42-PA10#v=onepage&q&f=false[/url]




[color=green]要先載入diag才能使用diag指令[/color]
[color=red](%i1)[/color]　[color=blue]load(diag);[/color]
[color=red](%o1)[/color]　[i]C:\maxima-5.41.0\share\maxima\5.41.0\share\contrib\diag.mac[/i]

[color=green]針對三對角矩陣Implicit Shift副程式[/color]
[color=red](%i2)[/color]
[color=blue]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)
)$[/color]

[color=green]三對角矩陣[/color]
[color=red](%i3)[/color]
[color=blue]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]);[/color]
[color=red](%o3)[/color]　\(\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] \)

[color=green]執行三對角矩陣的Implicit Shift[/color]
[color=red](%i4)[/color]　[color=blue]TridiagonalImplicitShift(T);[/color]
\( \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] \)
[color=red](%o4)[/color]　\( \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] \)

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

bugmens 發表於 2017-11-4 09:54

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

取矩陣\( A \)的右下角\( 2 \times 2 \)子矩陣\( \left[ \matrix{A[n-1,n-1]& A[n-1,n]\cr A[n,n-1]&A[n,n]} \right] \)，
計算特徵方程式\( \left| \matrix{A[n-1,n-1]-x& A[n-1,n] \cr A[n,n-1]&A[n,n]-x} \right|=0 \)，

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

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

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

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

(1)若特徵值\( \sigma_1,\sigma_2 \)為實數，直接計算\( (A-\sigma_1I)(A- \sigma_2I) \)

(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 \)矩陣
}
－－－－－－－－－－－－－－－－－－－－－－－－－－－
參考資料[url]https://books.google.com.tw/books?id=1V9PbyYGZIIC&lpg=PA327&dq=%22implicit%20DOUBLE%22%20SHIFT%20QR%20%20householder&hl=zh-TW&pg=PA330#v=onepage&q&f=false[/url]



[color=green]要先載入diag才能使用diag指令[/color]
[color=red](%i1)[/color]　[color=blue]load(diag);[/color]
[color=red](%o1)[/color]　[i]C:\maxima-5.40.0\share\maxima\5.40.0\share\contrib\diag.mac[/i]

[color=green]轉換成Householder矩陣副程式[/color]
[color=red](%i2)[/color]
[color=blue]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,
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)
)$[/color]

[color=green]針對Hessenberg矩陣Implicit Double Shift副程式[/color]
[color=red](%i3)[/color]
[color=blue]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))
)$[/color]

[color=green]Hessenberg矩陣，各元數字經過設計，\( \left[ \matrix{x \cr y \cr z} \right] \)向量長度均為完全平方數，處理後的\( Hessenberg \)矩陣各元最多是分數，不會出現根號數字[/color]
[color=red](%i4)[/color]
[color=blue]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]);[/color]
[color=red](%o4)[/color]　\( \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] \)

[color=green]執行\(Hessenberg\)矩陣的\(Implicit\) \(Double\) \(Shift\)[/color]
[color=red](%i5)[/color]　[color=blue]HessenbergImplicitDoubleShift(H);[/color]
\( 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]\)
[color=red](%o5)[/color]　\( \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 發表於 2018-11-11 12:02

3-5.減少QR法運算-Hessenberg矩陣的Implicit Multishift
前一篇\( Double Shift \)文章取矩陣\(H\)的右下角\(2 \times 2\)子矩陣得到2個特徵值\( \mu_1,\mu_2 \)，以\((H-\mu_1I)(H-\mu_2I)\)來計算\(QR\)分解。
而\(Multishift\)則是取矩陣\(H\)的右下角\(k \times k\)子矩陣得到\(k\)個特徵值\(\mu_1,\mu_2,\ldots,\mu_{k}\)，以\((H-\mu_1I)(H-\mu_2I)...(H-\mu_{k}I)\)來計算\(QR\)分解。
\(k\)個特徵值中若\(\mu\)為實數則直接計算\((H-\mu I)\)，若\(\mu,\overline{\mu}\)互為共軛複數，則計算\((H-\mu I)(H-\overline{\mu}I)=H^2-2Re(\mu)H+|\; \mu |\;^2I\)以避免複數運算。

以\(n=9\)舉例說明
將\(A\)矩陣轉換成\(Hessenberg\)矩陣
\( A=\left[ \matrix{\times & \times & \times & \times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times & \times & \times & \times} \right] \Rightarrow H=\left[ \matrix{\times & \times & \times & \times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times & \times & \times & \times \cr
0 & \times & \times & \times & \times & \times & \times & \times & \times \cr
0 & 0 & \times & \times & \times & \times & \times & \times & \times \cr
0 & 0 & 0 & \times & \times & \times & \times & \times & \times \cr
0 & 0 & 0 & 0 & \times & \times & \times & \times & \times \cr
0 & 0 & 0 & 0 & 0 & \times & \times & \times & \times \cr
0 & 0 & 0 & 0 & 0 & 0 & \times & \times & \times \cr
0 & 0 & 0 & 0 & 0 & 0 & 0 & \times & \times} \right] \)

取\(H\)右下角\(4 \times 4\)子矩陣\(T=\left[ \matrix{\times & \times & \times & \times \cr \times & \times & \times & \times \cr 0 & \times & \times & \times \cr 0 & 0 & \times & \times} \right]\)為例
計算\(det(T-\mu I)=0\)，得到4個特徵值\(\mu=\mu_1,\mu_2,\mu_3,\mu_4\)
計算\((H-\mu_1I)(H-\mu_2I)(H-\mu_3I)(H-\mu_4I)e1=\left[ \matrix{\times \cr \times \cr \times \cr \times \cr \times \cr 0 \cr 0 \cr 0 \cr 0} \right]\)
其中\(e1\)為第一元為1，其餘為0的向量，換句話說就是取\( (H-\mu_1I)(H-\mu_2I)(H-\mu_3I)(H-\mu_4I) \)相乘後的第1行

取前5個非0元形成向量\( \left[ \matrix{\times \cr \times \cr \times \cr \times \cr \times} \right] \)找一個\( 5\times 5 Householder \)矩陣\( \hat{P_0} \)，使得\( \hat{P_0}\left[ \matrix{\times \cr \times \cr \times \cr \times \cr \times} \right]=\left[ \matrix{* \cr 0 \cr 0 \cr 0 \cr 0} \right] \)
設\( P_0=diag(\hat{P_0},I_4) \)，計算\( H_0=P_0^T H P_0 =\left[\matrix{
\times & \times & \times & \times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times & \times & \times & \times \cr
+ & \times & \times & \times & \times & \times & \times & \times & \times \cr
+ & + & \times & \times & \times & \times & \times & \times & \times \cr
+ & + & + & \times & \times & \times & \times & \times & \times \cr
+ & + & + & + & \times & \times & \times & \times & \times \cr
0 & 0 & 0 & 0 & 0 & \times & \times & \times & \times \cr
0 & 0 & 0 & 0 & 0 & 0 & \times & \times & \times \cr
0 & 0 & 0 & 0 & 0 & 0 & 0 & \times & \times} \right] \)　其中\( \matrix{+& & & \cr + & + & & \cr + &+&+&\cr + & + & +&+} \)形成凸起(bulge)

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

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

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

取向量\( \left[ \matrix{H_3[5,4] \cr H_3[6,4] \cr H_3[7,4] \cr H_3[8,4] \cr H_2[9,4]} \right]=\left[ \matrix{\times \cr + \cr + \cr + \cr +} \right] \)找一個\( 5\times 5 Householder \)矩陣\( \hat{P_4} \)，使得\( \hat{P_4}\left[ \matrix{\times \cr + \cr + \cr + \cr +} \right]=\left[ \matrix{* \cr 0 \cr 0 \cr 0 \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 & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times & \times & \times & \times \cr
0 & \times & \times & \times & \times & \times & \times & \times & \times \cr
0 & 0 & \times & \times & \times & \times & \times & \times & \times \cr
0 & 0 & 0 & \times & \times & \times & \times & \times & \times \cr
0 & 0 & 0 & 0 & \times & \times & \times & \times & \times \cr
0 & 0 & 0 & 0 & + & \times & \times & \times & \times \cr
0 & 0 & 0 & 0 & + & + & \times & \times & \times \cr
0 & 0 & 0 & 0 & + & + & + & \times & \times} \right] \)　將凸起\( \matrix{+& & \cr +&+& \cr +&+&+} \)往右下角推

取向量\( \left[ \matrix{H_4[6,5] \cr H_4[7,5] \cr H_4[8,5] \cr H_4[9,5]} \right]=\left[ \matrix{\times \cr + \cr + \cr +} \right] \)找一個\( 4\times 4 Householder \)矩陣\( \hat{P_5} \)，使得\( \hat{P_5}\left[ \matrix{\times \cr + \cr + \cr +} \right]=\left[ \matrix{* \cr 0 \cr 0 \cr 0} \right] \)
設\( P_5=diag(I_5,\hat{P_5}) \)，計算\( H_5=P_5^T H_4 P_5 =\left[\matrix{
\times & \times & \times & \times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times & \times & \times & \times \cr
0 & \times & \times & \times & \times & \times & \times & \times & \times \cr
0 & 0 & \times & \times & \times & \times & \times & \times & \times \cr
0 & 0 & 0 & \times & \times & \times & \times & \times & \times \cr
0 & 0 & 0 & 0 & \times & \times & \times & \times & \times \cr
0 & 0 & 0 & 0 & 0 & \times & \times & \times & \times \cr
0 & 0 & 0 & 0 & 0 & + & \times & \times & \times \cr
0 & 0 & 0 & 0 & 0 & + & + & \times & \times} \right] \)　將凸起\( \matrix{+& \cr +&+} \)往右下角推

取向量\( \left[ \matrix{H_5[7,6] \cr H_5[8,6] \cr H_5[9,6]} \right]=\left[ \matrix{\times \cr + \cr +} \right] \)找一個\( 3\times 3 Householder \)矩陣\( \hat{P_6} \)，使得\( \hat{P_6}\left[ \matrix{\times \cr + \cr +} \right]=\left[ \matrix{* \cr 0 \cr 0} \right] \)
設\( P_6=diag(I_6,\hat{P_6}) \)，計算\( H_6=P_6^T H_5 P_6 =\left[\matrix{
\times & \times & \times & \times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times & \times & \times & \times \cr
0 & \times & \times & \times & \times & \times & \times & \times & \times \cr
0 & 0 & \times & \times & \times & \times & \times & \times & \times \cr
0 & 0 & 0 & \times & \times & \times & \times & \times & \times \cr
0 & 0 & 0 & 0 & \times & \times & \times & \times & \times \cr
0 & 0 & 0 & 0 & 0 & \times & \times & \times & \times \cr
0 & 0 & 0 & 0 & 0 & 0 & \times & \times & \times \cr
0 & 0 & 0 & 0 & 0 & 0 & + & \times & \times} \right] \)　將凸起\(+\)往右下角推

取向量\( \left[ \matrix{H_6[8,7] \cr H_6[9,7]} \right]=\left[ \matrix{\times \cr +} \right] \)找一個\( 2\times 2 Householder \)矩陣\( \hat{P_7} \)，使得\( \hat{P_7}\left[ \matrix{\times \cr +} \right]=\left[ \matrix{* \cr 0} \right] \)
設\( P_7=diag(I_7,\hat{P_7}) \)，計算\( H_7=P_7^T H_6 P_7 =\left[\matrix{
\times & \times & \times & \times & \times & \times & \times & \times & \times \cr
\times & \times & \times & \times & \times & \times & \times & \times & \times \cr
0 & \times & \times & \times & \times & \times & \times & \times & \times \cr
0 & 0 & \times & \times & \times & \times & \times & \times & \times \cr
0 & 0 & 0 & \times & \times & \times & \times & \times & \times \cr
0 & 0 & 0 & 0 & \times & \times & \times & \times & \times \cr
0 & 0 & 0 & 0 & 0 & \times & \times & \times & \times \cr
0 & 0 & 0 & 0 & 0 & 0 & \times & \times & \times \cr
0 & 0 & 0 & 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) \)
}

theFirstColumnofΠk(H,Bulge)
{
\(H\)為\(n \times n\)的\(Hessenberg\)矩陣，取右下角\(Bulge \times Bulge\)子矩陣\(T\)
計算子矩陣\(T\)的特徵值\(\lambda=\mu_1,\mu_2,\ldots,\mu_{Bulge}\)
若特徵值\(\mu\)為實數則計算\((H-\mu I)\)
若特徵值\(\mu,\overline{\mu}\)為共軛複數則計算\((H-\mu I)(H-\overline{\mu}I)=H^2-2Re(\mu)H+|\;\mu|\;^2I\)
\(Column1=(H-\mu_1I)(H-\mu_2I)\ldots(H-\mu_{Bulge}I)e_1\)，其中\(e_1=\left[\matrix{1 \cr 0 \cr \vdots \cr 0} \right]\)
換句話說取\( \displaystyle \prod_{i=1}^{Bulge}(H-\mu_i I)\)的第一行向量
取\(Column1\)前\(Bulge+1\)個非零0元得到向量\(X\)
}

HessenbergImplicitMultishift(H,Bulge)
{
計算\( \displaystyle \prod_{i=1}^{Bulge}(H-\mu_i I)\)的第一行向量
theFirstColumnofΠk(H,Bulge)
for k=0 ~ (n-2)
  {\(\hat{P_k}=Householder(X)\)，找一個\(Householder\)矩陣\(\hat{P_k}\)，使得\(\hat{P_k}X=\left[ \matrix{* \cr 0 \cr \vdots \cr 0} \right]\)，*為非零數字
   \(P_k=diag(I_k,\hat{P_k},I_{n-Bulge-1-k})\)
   \(H=P_k^THP_k\)
   if \(k=n-2\) then，已將凸起完全消除
     {\(return(H)\)，\(H\)仍是\(Hessenberg\)矩陣
     }
   若\(k=0 \sim (n-Bulge-2)\)
     {\(Xrows=Bulge+1\)}
   若\(k=(n-Bulge-1) \sim(n-3)\)
     {\(Xrows=n-k-1\)}
   取\(X=\left[ \matrix{H[k+1+1,k+1] \cr H[k+1+2,k+1] \cr \vdots \cr H[k+1+Xrows,k+1]} \right]\)做為下一次Householder(X)使用
  }
}
－－－－－－－－－－－－－－－－－－－－－－－－－－－
參考資料
Z. Bai and J.W. Demmel(1989). "On a Block Implementation of Hessenberg Multishift QR Iteration," Int. J. High Speed Comput. 1,97-112.
[url]http://www.netlib.org/lapack/lawnspdf/lawn08.pdf[/url]


[color=green]要先載入diag才能使用diag指令[/color]
[color=red](%i1)[/color]　[color=blue]load(diag);[/color]
[color=red](%o1)[/color]　[i]C:\maxima-5.41.0\share\maxima\5.41.0\share\contrib\diag.mac[/i]

[color=green]轉換成Householder矩陣副程式[/color]
[color=red](%i2)[/color]
[color=blue]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,
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)
)$[/color]

[color=green]計算(H-u1I)(H-u2I)...(H-ukI)第一行[/color]
[color=red](%i3)[/color]
[color=blue]theFirstColumnofΠk(H,Bulge):=block
([n:length(H),T,equation,column1,mu,eigenvalue,realpart,cabs,HH,X],
print("取H矩陣右下角",Bulge,"*",Bulge,"子矩陣T=",
       T:genmatrix(lambda([i,j],H[n-Bulge+i,n-Bulge+j]),Bulge,Bulge)),
print("計算T的特徵值"),
print("det(T-μ*I)=0,",equation:T-%mu*ident(Bulge),"=0"),
print("展開行列式",equation:factor(determinant(equation)),"=0"),
print("特徵值μ=",eigenvalue:map(rhs,solve(equation,%mu))),
print("第1元為1,其餘為0的向量e1=column1=",column1:genmatrix(lambda([i,j],if i=1 then 1 else 0),n,1)),
for i:1 thru Bulge do
    (mu:eigenvalue[ i ],
     if numberp(mu)=true then
        (print("當特徵值為實數μ=",mu),
         print("計算column1=(H-",mu,"I).column1=",column1: (H-mu*ident(n)).column1)
        )
     else
        (realpart:realpart(mu),
         cabs:cabs(mu)^2,
         print("當特徵值為共軛複數μ=",eigenvalue[ i ],",",eigenvalue[i+1],",Re(μ)=",realpart,",|μ|^2=",cabs),
         print("column1=(H^2-2*(",realpart,")H+",cabs,"I).column1"),
         HH:H.H-2*realpart*H+cabs*ident(n),
         print("column1=",HH,".",column1,"=",column1:HH.column1),
         i:i+1
        )
    ),
print("column1=(H-μ1I)(H-μ2I)...(H-μ",Bulge,"I)e1=",column1),
print("擷取第1行前",Bulge+1,"元數字X=",X:genmatrix(lambda([i,j],column1[i,1]),Bulge+1,1)),
return(X)
)$[/color]

[color=green]針對Hessenberg矩陣Implicit Multishift副程式[/color]
[color=red](%i4)[/color]
[color=blue]HessenbergImplicitMultishift(H,X):=block
([n:length(H),Bulge:length(X)-1,P_hat_k,P_k,Xrows:length(X)],
for k:0 thru n-2 do
  (print("k=",k),
   print("取向量",X,"找一個Householder矩陣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-Bulge-1)]))
   else if k>n-Bulge-2 then
     (P_k:diag([ident(k),P_hat_k]))
   else
     (P_k:diag([ident(k),P_hat_k,ident(n-k-Bulge-1)])),
   print("設P_",k,"=",P_k),
   print("計算H_",k,"=P_",k,"^T.H_",k-1,".P_",k,"=",H:ratsimp(transpose(P_k).H.P_k)),
   if k=n-2 then
     (numer:false,
      return(H)
     ),
   if k+Bulge>n-2 then
     (Xrows:n-k-1),
   print("取向量",genmatrix(lambda([i,j],concat("H",k,"[",i+k+1,",",k+1,"]")),Xrows,1),"=",
                   X:genmatrix(lambda([i,j],H[i+k+1,k+1]),Xrows,1)),
   print("------------------------"),
   if numer=false and ratnump(sqrt(X.X))=false then
     (print("向量X長度=",sqrt(X.X),"不是完全平方數,改用浮點數計算"),
      numer:true,
      ratprint:false
     )
  )
)$[/color]

[color=green]Hessenberg矩陣[/color]
[color=red](%i5)[/color]
[color=blue]H:matrix([-2,1,-1, 1,2,-2/3,-3,0,29/6],
             [-1,-1,1,2,-2,1,2/3,-5/2,-23/6],
             [0,-1,1,2,-1,1,6,-8/3,-65/6],
             [0,0,-1,-2,0,3,-1,7/2,10/3],
             [0,0,0,-1,-1,2,-1,-1,-4],
             [0,0,0,0,-2,1,2,2,-1],
             [0,0,0,0,0,-1,-2,0,1],
             [0,0,0,0,0,0,1,-1,-2],
             [0,0,0,0,0,0,0,1,-1] );[/color]
[color=red](%o5)[/color]　\( \left[ \matrix{\displaystyle -2&1&-1&1&2&-\frac{2}{3}&-3&0&\frac{29}{6}\cr
-1&-1&1&2&-2&1&\frac{2}{3}&-\frac{5}{2}&-\frac{23}{6}\cr
0&-1&1&2&-1&1&6&-\frac{8}{3}&-\frac{65}{6}\cr
0&0&-1&-2&0&3&-1&\frac{7}{2}&\frac{10}{3}\cr
0&0&0&-1&-1&2&-1&-1&-4\cr
0&0&0&0&-2&1&2&2&-1\cr
0&0&0&0&0&-1&-2&0&1\cr
0&0&0&0&0&0&1&-1&-2\cr
0&0&0&0&0&0&0&1&-1} \right] \)

[color=green]取H右下角4*4矩陣計算特徵值，計算(H-u1I)(H-u2I)...(H-u4I)第一行[/color]
[color=red](%i6)[/color]　[color=blue]X:theFirstColumnofΠk(H,4);[/color]
取\(H\)矩陣右下角\(4*4\)子矩陣\(T=\left[ \matrix{1&2&2&-1 \cr-1&-2&0&1 \cr0&1&-1&-2 \cr0&0&1&-1}\right] \)
計算\(T\)的特徵值
\(det(T-\mu*I)=0,\left[ \matrix{1-\mu&2&2&-1 \cr-1&-\mu-2&0&1 \cr0&1&-\mu-1&-2 \cr0&0&1&-\mu-1}\right] \)=0
展開行列式\((\mu^2+\mu+1)(\mu^2+2\mu+2)=0\)
特徵值\( \displaystyle \mu=[-\%i-1,\%i-1,-\frac{\sqrt{3}\%i+1}{2},\frac{\sqrt{3}\%i-1}{2} \)
第1元為1,其餘為0的向量\(e1=column1=\left[ \matrix{1 \cr0 \cr0 \cr0 \cr0 \cr0 \cr0 \cr0 \cr0}\right] \)
當特徵值為共軛複數\(\mu=-\%i-1,\%i-1,Re(\mu)=-1,|\; \mu |\;^2=2 \)
\(column1=(H^2-2*(-1)H+2I).column1\)
\(column1=\left[ \matrix{ \displaystyle 1&0&-1&-4&-\frac{5}{3}&\frac{28}{3}&-\frac{11}{3}&\frac{31}{6}&-\frac{29}{6} \cr1&-1&1&1&-5&5&\frac{47}{6}&\frac{9}{2}&\frac{11}{3} \cr1&-2&2&1&-2&1&\frac{11}{3}&-\frac{11}{3}&\frac{19}{6} \cr0&1&-1&0&-5&3&\frac{11}{2}&\frac{17}{2}&-\frac{7}{2} \cr0&0&1&1&-3&2&5&-\frac{7}{2}&-\frac{13}{3} \cr0&0&0&2&-4&-1&6&5&4 \cr0&0&0&0&2&-1&0&-1&0 \cr0&0&0&0&0&-1&-1&-1&1 \cr0&0&0&0&0&0&1&0&-1}\right].
\left[ \matrix{1 \cr0 \cr0 \cr0 \cr0 \cr0 \cr0 \cr0 \cr0}\right]=\left[ \matrix{ \displaystyle 1 \cr1 \cr1 \cr0 \cr0 \cr0 \cr0 \cr0 \cr0}\right] \)
當特徵值為共軛複數\( \displaystyle \mu=-\frac{\sqrt{3}\%i+1}{2},\frac{\sqrt{3}\%i-1}{2},Re(\mu)=-\frac{1}{2},|\; \mu |\;^2=1 \)
\(column1=(\displaystyle H^2-2*(-\frac{1}{2})H+1I).column1\)
\(column1=\left[ \matrix{ \displaystyle 2&-1&0&-5&-\frac{11}{3}&10&-\frac{2}{3}&\frac{31}{6}&-\frac{29}{3} \cr2&-1&0&-1&-3&4&\frac{43}{6}&7&\frac{15}{2} \cr1&-1&0&-1&-1&0&-\frac{7}{3}&-1&14 \cr0&1&0&1&-5&0&\frac{13}{2}&5&-\frac{41}{6} \cr0&0&1&2&-3&0&6&-\frac{5}{2}&-\frac{1}{3} \cr0&0&0&2&-2&-3&4&3&5 \cr0&0&0&0&2&0&1&-1&-1 \cr0&0&0&0&0&-1&-2&-1&3 \cr0&0&0&0&0&0&1&-1&-1}\right].
\left[ \matrix{ \displaystyle 1 \cr1 \cr1 \cr0 \cr0 \cr0 \cr0 \cr0 \cr0}\right]=
\left[ \matrix{ \displaystyle 1 \cr1 \cr0 \cr1 \cr1 \cr0 \cr0 \cr0 \cr0}\right] \)
\(column1=(H-\mu 1I)(H-\mu 2I)...(H-\mu 4I)e1=\left[ \matrix{ \displaystyle 1 \cr1 \cr0 \cr1 \cr1 \cr0 \cr0 \cr0 \cr0}\right] \)
擷取第1行前5元數字\(X=\left[ \matrix{ \displaystyle 1 \cr1 \cr0 \cr1 \cr1}\right]\)
[color=red](%o6)[/color]　\(\left[ \matrix{1 \cr1 \cr0 \cr1 \cr1}\right] \)

[color=red](%i7)[/color]　[color=blue]HessenbergImplicitMultishift(H,X);[/color]
\(k=0\)
\(X\)向量\(=\left[ \matrix{\displaystyle 1\cr 1\cr 0\cr 1\cr 1}\right]\)第1元\(=1,sign(X_1)=1\),長度\( \Vert\;X \Vert\;=sqrt(\left[ \matrix{\displaystyle 1\cr 1\cr 0\cr 1\cr 1}\right].\left[ \matrix{\displaystyle 1\cr 1\cr 0\cr 1\cr 1}\right] )=2\)
第1元是1,其餘為0的向量\( e1=\left[ \matrix{\displaystyle 1\cr 0\cr 0\cr 0\cr 0}\right] \)
\(X+sign(X_1)*\Vert\; X \Vert\;*e1=\left[ \matrix{\displaystyle 1\cr 1\cr 0\cr 1\cr 1}\right]+1*2*\left[ \matrix{\displaystyle 1\cr 0\cr 0\cr 0\cr 0}\right]=\left[ \matrix{\displaystyle 3\cr 1\cr 0\cr 1\cr 1}\right] \)
長度\(\Vert\; X+sign(X_1)*\Vert\; X \Vert\;*e1 \Vert\;=sqrt(\left[ \matrix{\displaystyle 3\cr 1\cr 0\cr 1\cr 1}\right].\left[ \matrix{\displaystyle 3\cr 1\cr 0\cr 1\cr 1}\right] )=2\sqrt{3} \)
\(v=\left[ \matrix{\displaystyle 3\cr 1\cr 0\cr 1\cr 1}\right]/2\sqrt{3}=\left[\matrix{\displaystyle \frac{\sqrt{3}}{2}\cr \frac{1}{2\sqrt{3}}\cr 0\cr \frac{1}{2\sqrt{3}}\cr \frac{1}{2\sqrt{3}}}\right] \)
\( \hat{P}=I-2*v.v^T=\left[ \matrix{\displaystyle 1&0&0&0&0\cr 0&1&0&0&0\cr 0&0&1&0&0\cr 0&0&0&1&0\cr 0&0&0&0&1}\right]-2*\left[ \matrix{\displaystyle \frac{\sqrt{3}}{2}\cr \frac{1}{2\sqrt{3}}\cr 0\cr \frac{1}{2\sqrt{3}}\cr \frac{1}{2\sqrt{3}}}\right].\left[ \matrix{\displaystyle \frac{\sqrt{3}}{2}&\frac{1}{2\sqrt{3}}&0&\frac{1}{2\sqrt{3}}&\frac{1}{2\sqrt{3}}}\right]=\left[ \matrix{\displaystyle -\frac{1}{2}&-\frac{1}{2}&0&-\frac{1}{2}&-\frac{1}{2}\cr -\frac{1}{2}&\frac{5}{6}&0&-\frac{1}{6}&-\frac{1}{6}\cr 0&0&1&0&0\cr -\frac{1}{2}&-\frac{1}{6}&0&\frac{5}{6}&-\frac{1}{6}\cr -\frac{1}{2}&-\frac{1}{6}&0&-\frac{1}{6}&\frac{5}{6}}\right] \)
取向量\( \left[ \matrix{\displaystyle 1\cr 1\cr 0\cr 1\cr 1}\right] \)找一個\(Householder\)矩陣\( \hat{P}_0\),使得\(\left[ \matrix{\displaystyle -\frac{1}{2}&-\frac{1}{2}&0&-\frac{1}{2}&-\frac{1}{2}\cr -\frac{1}{2}&\frac{5}{6}&0&-\frac{1}{6}&-\frac{1}{6}\cr 0&0&1&0&0\cr -\frac{1}{2}&-\frac{1}{6}&0&\frac{5}{6}&-\frac{1}{6}\cr -\frac{1}{2}&-\frac{1}{6}&0&-\frac{1}{6}&\frac{5}{6}}\right] \left[ \matrix{1\cr 1\cr 0\cr 1\cr 1}\right]=\left[ \matrix{\displaystyle -2\cr 0\cr 0\cr 0\cr 0}\right] \)
設\(P_0=\left[ \matrix{\displaystyle -\frac{1}{2}&-\frac{1}{2}&0&-\frac{1}{2}&-\frac{1}{2}&0&0&0&0\cr -\frac{1}{2}&\frac{5}{6}&0&-\frac{1}{6}&-\frac{1}{6}&0&0&0&0\cr 0&0&1&0&0&0&0&0&0\cr -\frac{1}{2}&-\frac{1}{6}&0&\frac{5}{6}&-\frac{1}{6}&0&0&0&0\cr -\frac{1}{2}&-\frac{1}{6}&0&-\frac{1}{6}&\frac{5}{6}&0&0&0&0\cr 0&0&0&0&0&1&0&0&0\cr 0&0&0&0&0&0&1&0&0\cr 0&0&0&0&0&0&0&1&0\cr 0&0&0&0&0&0&0&0&1}\right] \)
計算\(H_0=P_0^T.H_{-1}.P_0=\left[ \matrix{\displaystyle -1&-\frac{5}{6}&\frac{1}{2}&-\frac{5}{6}&-\frac{1}{3}&-\frac{8}{3}&\frac{13}{6}&0&-\frac{1}{6}\cr
1&-\frac{19}{18}&\frac{3}{2}&\frac{35}{18}&-\frac{20}{9}&\frac{1}{3}&\frac{43}{18}&-\frac{5}{2}&-\frac{11}{2}\cr
0&-1&1&2&-1&1&6&-\frac{8}{3}&-\frac{65}{6}\cr
1&-\frac{7}{18}&-\frac{1}{2}&-\frac{43}{18}&-\frac{5}{9}&\frac{7}{3}&\frac{13}{18}&\frac{7}{2}&\frac{5}{3}\cr
1&-\frac{7}{18}&\frac{1}{2}&-\frac{25}{18}&-\frac{14}{9}&\frac{4}{3}&\frac{13}{18}&-1&-\frac{17}{3}\cr 1&\frac{1}{3}&0&\frac{1}{3}&-\frac{5}{3}&1&2&2&-1\cr 0&0&0&0&0&-1&-2&0&1\cr 0&0&0&0&0&0&1&-1&-2\cr 0&0&0&0&0&0&0&1&-1}\right] \)
取向量\( \left[ \matrix{\displaystyle H0[2,1]\cr H0[3,1]\cr H0[4,1]\cr H0[5,1]\cr H0[6,1]}\right]=\left[ \matrix{\displaystyle 1\cr 0\cr 1\cr 1\cr 1}\right] \)
------------------------
\(k=1\)
\(X\)向量\(=\left[ \matrix{\displaystyle 1\cr 0\cr 1\cr 1\cr 1}\right] \)第1元\(=1,sign(X_1)=1\),長度\( \Vert\;X \Vert\;=sqrt(\left[ \matrix{\displaystyle 1\cr 0\cr 1\cr 1\cr 1}\right].\left[ \matrix{\displaystyle 1\cr 0\cr 1\cr 1\cr 1}\right] )=2\)
第1元是1,其餘為0的向量\( e1=\left[ \matrix{\displaystyle 1\cr 0\cr 0\cr 0\cr 0}\right] \)
\(X+sign(X_1)*\Vert\; X \Vert\;*e1=\left[ \matrix{\displaystyle 1\cr 0\cr 1\cr 1\cr 1}\right]+1*2*\left[ \matrix{\displaystyle 1\cr 0\cr 0\cr 0\cr 0}\right]=\left[ \matrix{\displaystyle 3\cr 0\cr 1\cr 1\cr 1}\right] \)
長度\(\Vert\; X+sign(X_1)*\Vert\; X \Vert\;*e1 \Vert\;=sqrt(\left[ \matrix{\displaystyle 3\cr 0\cr 1\cr 1\cr 1}\right].\left[ \matrix{\displaystyle 3\cr 0\cr 1\cr 1\cr 1}\right] )=2\sqrt{3}\)
\(v=\left[ \matrix{\displaystyle 3\cr 0\cr 1\cr 1\cr 1}\right]/2\sqrt{3}=\left[ \matrix{\displaystyle \frac{\sqrt{3}}{2}\cr 0\cr \frac{1}{2\sqrt{3}}\cr \frac{1}{2\sqrt{3}}\cr \frac{1}{2\sqrt{3}}}\right] \)
\(\hat{P}=I-2*v.v^T=\left[ \matrix{\displaystyle 1&0&0&0&0\cr 0&1&0&0&0\cr 0&0&1&0&0\cr 0&0&0&1&0\cr 0&0&0&0&1}\right]-2*\left[ \matrix{\displaystyle \frac{\sqrt{3}}{2}\cr 0\cr \frac{1}{2\sqrt{3}}\cr \frac{1}{2\sqrt{3}}\cr \frac{1}{2\sqrt{3}}}\right].\left[ \matrix{\displaystyle \frac{\sqrt{3}}{2}&0&\frac{1}{2\sqrt{3}}&\frac{1}{2\sqrt{3}}&\frac{1}{2\sqrt{3}}}\right]=\left[\matrix{\displaystyle-\frac{1}{2}&0&-\frac{1}{2}&-\frac{1}{2}&-\frac{1}{2}\cr 0&1&0&0&0\cr -\frac{1}{2}&0&\frac{5}{6}&-\frac{1}{6}&-\frac{1}{6}\cr -\frac{1}{2}&0&-\frac{1}{6}&\frac{5}{6}&-\frac{1}{6}\cr -\frac{1}{2}&0&-\frac{1}{6}&-\frac{1}{6}&\frac{5}{6}}\right] \)
取向量\( \left[ \matrix{\displaystyle 1\cr 0\cr 1\cr 1\cr 1}\right] \)找一個\(Householder\)矩陣\(\hat{P}_1\),使得\(\left[ \matrix{\displaystyle -\frac{1}{2}&0&-\frac{1}{2}&-\frac{1}{2}&-\frac{1}{2}\cr 0&1&0&0&0\cr -\frac{1}{2}&0&\frac{5}{6}&-\frac{1}{6}&-\frac{1}{6}\cr -\frac{1}{2}&0&-\frac{1}{6}&\frac{5}{6}&-\frac{1}{6}\cr -\frac{1}{2}&0&-\frac{1}{6}&-\frac{1}{6}&\frac{5}{6}}\right] \left[ \matrix{1\cr 0\cr 1\cr 1\cr 1}\right]=\left[ \matrix{\displaystyle -2\cr 0\cr 0\cr 0\cr 0}\right] \)
設\(P_1=\left[ \matrix{\displaystyle 1&0&0&0&0&0&0&0&0\cr 0&-\frac{1}{2}&0&-\frac{1}{2}&-\frac{1}{2}&-\frac{1}{2}&0&0&0\cr 0&0&1&0&0&0&0&0&0\cr 0&-\frac{1}{2}&0&\frac{5}{6}&-\frac{1}{6}&-\frac{1}{6}&0&0&0\cr 0&-\frac{1}{2}&0&-\frac{1}{6}&\frac{5}{6}&-\frac{1}{6}&0&0&0\cr 0&-\frac{1}{2}&0&-\frac{1}{6}&-\frac{1}{6}&\frac{5}{6}&0&0&0\cr 0&0&0&0&0&0&1&0&0\cr 0&0&0&0&0&0&0&1&0\cr 0&0&0&0&0&0&0&0&1}\right] \)
計算\(H_1=P_1^T.H_0.P_1=\left[ \matrix{\displaystyle -1&\frac{7}{3}&\frac{1}{2}&\frac{2}{9}&\frac{13}{18}&-\frac{29}{18}&\frac{13}{6}&0&-\frac{1}{6}\cr -2&-1&-\frac{3}{4}&\frac{1}{6}&\frac{29}{12}&-\frac{37}{12}&-\frac{35}{12}&-1&\frac{21}{4}\cr 0&-\frac{1}{2}&1&\frac{13}{6}&-\frac{5}{6}&\frac{7}{6}&6&-\frac{8}{3}&-\frac{65}{6}\cr 0&0&-\frac{5}{4}&-\frac{463}{162}&\frac{361}{324}&\frac{427}{324}&-\frac{113}{108}&4&\frac{21}{4}\cr 0&\frac{1}{2}&-\frac{1}{4}&-\frac{137}{81}&\frac{91}{324}&\frac{157}{324}&-\frac{113}{108}&-\frac{1}{2}&-\frac{25}{12}\cr 0&-\frac{1}{2}&-\frac{3}{4}&-\frac{44}{81}&-\frac{131}{324}&-\frac{137}{324}&\frac{25}{108}&\frac{5}{2}&\frac{31}{12}\cr 0&\frac{1}{2}&0&\frac{1}{6}&\frac{1}{6}&-\frac{5}{6}&-2&0&1\cr 0&0&0&0&0&0&1&-1&-2\cr 0&0&0&0&0&0&0&1&-1}\right] \)
取向量\( \left[ \matrix{\displaystyle H1[3,2]\cr H1[4,2]\cr H1[5,2]\cr H1[6,2]\cr H1[7,2]}\right]=\left[ \matrix{\displaystyle -\frac{1}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}}\right] \)
------------------------
\(k=2\)
\(X\)向量\(=\left[ \matrix{\displaystyle -\frac{1}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}}\right] \)第1元\(\displaystyle =-\frac{1}{2},sign(X_1)=-1\),長度\( \Vert\;X \Vert\;=sqrt(\left[ \matrix{\displaystyle -\frac{1}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}}\right].\left[ \matrix{\displaystyle -\frac{1}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}}\right] )=1\)
第1元是1,其餘為0的向量\( e1=\left[ \matrix{\displaystyle 1\cr 0\cr 0\cr 0\cr 0}\right] \)
\(X+sign(X_1)*\Vert\; X \Vert\;*e1=\left[ \matrix{\displaystyle -\frac{1}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}}\right]+-1*1*\left[ \matrix{\displaystyle 1\cr 0\cr 0\cr 0\cr 0}\right]=\left[ \matrix{\displaystyle -\frac{3}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}}\right] \)
長度\(\Vert\; X+sign(X_1)*\Vert\; X \Vert\;*e1 \Vert\;=sqrt(\left[ \matrix{\displaystyle -\frac{3}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}}\right].\left[ \matrix{\displaystyle -\frac{3}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}}\right] )=\sqrt{3}\)
\(v=\left[ \matrix{\displaystyle -\frac{3}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}}\right]/\sqrt{3}=\left[ \matrix{\displaystyle -\frac{\sqrt{3}}{2}\cr 0\cr \frac{1}{2\sqrt{3}}\cr -\frac{1}{2\sqrt{3}}\cr \frac{1}{2\sqrt{3}}}\right] \)
\(\hat{P}=I-2*v.v^T=\left[ \matrix{\displaystyle 1&0&0&0&0\cr 0&1&0&0&0\cr 0&0&1&0&0\cr 0&0&0&1&0\cr 0&0&0&0&1}\right]-2*\left[ \matrix{\displaystyle -\frac{\sqrt{3}}{2}\cr 0\cr \frac{1}{2\sqrt{3}}\cr -\frac{1}{2\sqrt{3}}\cr \frac{1}{2\sqrt{3}}}\right].\left[ \matrix{\displaystyle -\frac{\sqrt{3}}{2}&0&\frac{1}{2\sqrt{3}}&-\frac{1}{2\sqrt{3}}&\frac{1}{2\sqrt{3}}}\right]=\left[\matrix{\displaystyle -\frac{1}{2}&0&\frac{1}{2}&-\frac{1}{2}&\frac{1}{2}\cr 0&1&0&0&0\cr \frac{1}{2}&0&\frac{5}{6}&\frac{1}{6}&-\frac{1}{6}\cr -\frac{1}{2}&0&\frac{1}{6}&\frac{5}{6}&\frac{1}{6}\cr \frac{1}{2}&0&-\frac{1}{6}&\frac{1}{6}&\frac{5}{6}}\right] \)
取向量\( \left[ \matrix{\displaystyle -\frac{1}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}}\right]\)找一個\(Householder\)矩陣\( \hat{P}_2\),使得\(\left[ \matrix{\displaystyle -\frac{1}{2}&0&\frac{1}{2}&-\frac{1}{2}&\frac{1}{2}\cr 0&1&0&0&0\cr \frac{1}{2}&0&\frac{5}{6}&\frac{1}{6}&-\frac{1}{6}\cr -\frac{1}{2}&0&\frac{1}{6}&\frac{5}{6}&\frac{1}{6}\cr \frac{1}{2}&0&-\frac{1}{6}&\frac{1}{6}&\frac{5}{6}}\right] \left[ \matrix{\displaystyle -\frac{1}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}}\right]=\left[ \matrix{\displaystyle 1\cr 0\cr 0\cr 0\cr 0}\right] \)
設\(P_2=\left[ \matrix{\displaystyle 1&0&0&0&0&0&0&0&0\cr 0&1&0&0&0&0&0&0&0\cr 0&0&-\frac{1}{2}&0&\frac{1}{2}&-\frac{1}{2}&\frac{1}{2}&0&0\cr 0&0&0&1&0&0&0&0&0\cr 0&0&\frac{1}{2}&0&\frac{5}{6}&\frac{1}{6}&-\frac{1}{6}&0&0\cr 0&0&-\frac{1}{2}&0&\frac{1}{6}&\frac{5}{6}&\frac{1}{6}&0&0\cr 0&0&\frac{1}{2}&0&-\frac{1}{6}&\frac{1}{6}&\frac{5}{6}&0&0\cr 0&0&0&0&0&0&0&1&0\cr 0&0&0&0&0&0&0&0&1}\right] \)
計算\(H_2=P_2^T.H_1.P_2=\left[ \matrix{\displaystyle -1&\frac{7}{3}&2&\frac{2}{9}&\frac{2}{9}&-\frac{10}{9}&\frac{5}{3}&0&-\frac{1}{6}\cr -2&-1&\frac{5}{3}&\frac{1}{6}&\frac{29}{18}&-\frac{41}{18}&-\frac{67}{18}&-1&\frac{21}{4}\cr 0&1&-\frac{3}{2}&-\frac{85}{54}&\frac{34}{27}&-\frac{26}{27}&-\frac{38}{9}&-\frac{1}{6}&\frac{43}{12}\cr
0&0&0&-\frac{463}{162}&\frac{113}{162}&\frac{281}{162}&-\frac{79}{54}&4&\frac{21}{4}\cr 0&0&\frac{1}{2}&-\frac{4}{9}&-\frac{7}{18}&\frac{7}{6}&\frac{43}{18}&-\frac{4}{3}&-\frac{62}{9}\cr 0&0&-\frac{1}{2}&-\frac{145}{81}&-\frac{11}{162}&-\frac{125}{162}&-\frac{191}{54}&\frac{10}{3}&\frac{133}{18}\cr 0&0&\frac{1}{2}&\frac{229}{162}&-\frac{34}{81}&-\frac{19}{81}&\frac{41}{27}&-\frac{5}{6}&-\frac{137}{36}\cr 0&0&\frac{1}{2}&0&-\frac{1}{6}&\frac{1}{6}&\frac{5}{6}&-1&-2\cr 0&0&0&0&0&0&0&1&-1}\right] \)
取向量\( \left[ \matrix{\displaystyle H2[4,3]\cr H2[5,3]\cr H2[6,3]\cr H2[7,3]\cr H2[8,3]}\right]=\left[ \matrix{\displaystyle 0\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}\cr \frac{1}{2}}\right] \)
------------------------
\(k=3\)
\(X\)向量\(=\left[ \matrix{\displaystyle 0\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}\cr \frac{1}{2}}\right]\)第1元\(=0,sign(X_1)=1\),長度\( \Vert\;X \Vert\;=sqrt(\left[ \matrix{\displaystyle 0\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}\cr \frac{1}{2}}\right].\left[ \matrix{\displaystyle 0\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}\cr \frac{1}{2}}\right] )=1 \)
第1元是1,其餘為0的向量\( e1=\left[ \matrix{\displaystyle 1\cr 0\cr 0\cr 0\cr 0}\right] \)
\(X+sign(X_1)*\Vert\; X \Vert\;*e1=\left[ \matrix{\displaystyle 0\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}\cr \frac{1}{2}}\right]+1*1*\left[ \matrix{\displaystyle 1\cr 0\cr 0\cr 0\cr 0}\right]=\left[ \matrix{\displaystyle 1\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}\cr \frac{1}{2}}\right] \)
長度\(\Vert\; X+sign(X_1)*\Vert\; X \Vert\;*e1 \Vert\;=sqrt(\left[ \matrix{\displaystyle 1\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}\cr \frac{1}{2}}\right].\left[ \matrix{\displaystyle 1\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}\cr \frac{1}{2}}\right] )=\sqrt{2} \)
\(v=\left[ \matrix{\displaystyle 1\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}\cr \frac{1}{2}}\right]/\sqrt{2}=\left[ \matrix{\displaystyle \frac{1}{\sqrt{2}}\cr \frac{1}{2^{3/2}}\cr -\frac{1}{2^{3/2}}\cr \frac{1}{2^{3/2}}\cr \frac{1}{2^{3/2}}}\right] \)
\(\hat{P}=I-2*v.v^T=\left[ \matrix{\displaystyle 1&0&0&0&0\cr 0&1&0&0&0\cr 0&0&1&0&0\cr 0&0&0&1&0\cr 0&0&0&0&1}\right]-2*\left[ \matrix{\displaystyle \frac{1}{\sqrt{2}}\cr \frac{1}{2^{3/2}}\cr -\frac{1}{2^{3/2}}\cr \frac{1}{2^{3/2}}\cr \frac{1}{2^{3/2}}}\right].\left[ \matrix{\displaystyle \frac{1}{\sqrt{2}}&\frac{1}{2^{3/2}}&-\frac{1}{2^{3/2}}&\frac{1}{2^{3/2}}&\frac{1}{2^{3/2}}}\right]=\left[\matrix{\displaystyle 0&-\frac{1}{2}&\frac{1}{2}&-\frac{1}{2}&-\frac{1}{2}\cr -\frac{1}{2}&\frac{3}{4}&\frac{1}{4}&-\frac{1}{4}&-\frac{1}{4}\cr \frac{1}{2}&\frac{1}{4}&\frac{3}{4}&\frac{1}{4}&\frac{1}{4}\cr -\frac{1}{2}&-\frac{1}{4}&\frac{1}{4}&\frac{3}{4}&-\frac{1}{4}\cr -\frac{1}{2}&-\frac{1}{4}&\frac{1}{4}&-\frac{1}{4}&\frac{3}{4}}\right] \)
取向量\( \left[ \matrix{\displaystyle 0\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}\cr \frac{1}{2}}\right]\)找一個\(Householder\)矩陣\( \hat{P}_3\),使得\(\left[ \matrix{\displaystyle 0&-\frac{1}{2}&\frac{1}{2}&-\frac{1}{2}&-\frac{1}{2}\cr -\frac{1}{2}&\frac{3}{4}&\frac{1}{4}&-\frac{1}{4}&-\frac{1}{4}\cr \frac{1}{2}&\frac{1}{4}&\frac{3}{4}&\frac{1}{4}&\frac{1}{4}\cr -\frac{1}{2}&-\frac{1}{4}&\frac{1}{4}&\frac{3}{4}&-\frac{1}{4}\cr -\frac{1}{2}&-\frac{1}{4}&\frac{1}{4}&-\frac{1}{4}&\frac{3}{4}}\right] \left[ \matrix{\displaystyle 0\cr \frac{1}{2}\cr -\frac{1}{2}\cr \frac{1}{2}\cr \frac{1}{2}}\right]=\left[ \matrix{\displaystyle -1\cr 0\cr 0\cr 0\cr 0}\right] \)
設\(P_3=\left[ \matrix{\displaystyle 1&0&0&0&0&0&0&0&0\cr 0&1&0&0&0&0&0&0&0\cr 0&0&1&0&0&0&0&0&0\cr 0&0&0&0&-\frac{1}{2}&\frac{1}{2}&-\frac{1}{2}&-\frac{1}{2}&0\cr 0&0&0&-\frac{1}{2}&\frac{3}{4}&\frac{1}{4}&-\frac{1}{4}&-\frac{1}{4}&0\cr 0&0&0&\frac{1}{2}&\frac{1}{4}&\frac{3}{4}&\frac{1}{4}&\frac{1}{4}&0\cr 0&0&0&-\frac{1}{2}&-\frac{1}{4}&\frac{1}{4}&\frac{3}{4}&-\frac{1}{4}&0\cr 0&0&0&-\frac{1}{2}&-\frac{1}{4}&\frac{1}{4}&-\frac{1}{4}&\frac{3}{4}&0\cr 0&0&0&0&0&0&0&0&1}\right] \)
計算\(H_3=P_3^T.H_2.P_3=\left[ \matrix{\displaystyle -1&\frac{7}{3}&2&-\frac{3}{2}&-\frac{23}{36}&-\frac{1}{4}&\frac{29}{36}&-\frac{31}{36}&-\frac{1}{6}\cr -2&-1&\frac{5}{3}&\frac{5}{12}&\frac{125}{72}&-\frac{173}{72}&-\frac{259}{72}&-\frac{7}{8}&\frac{21}{4}\cr 0&1&-\frac{3}{2}&\frac{13}{12}&\frac{559}{216}&-\frac{55}{24}&-\frac{625}{216}&\frac{251}{216}&\frac{43}{12}\cr 0&0&-1&-\frac{1}{4}&\frac{55}{54}&-\frac{3}{2}&-\frac{193}{54}&\frac{103}{27}&\frac{241}{24}\cr0&0&0&\frac{1}{2}&-\frac{529}{1296}&-\frac{13}{48}&\frac{1495}{1296}&-\frac{2081}{1296}&-\frac{647}{144}\cr0&0&0&-\frac{1}{2}&\frac{1385}{1296}&-\frac{65}{144}&-\frac{1535}{1296}&\frac{6121}{1296}&\frac{719}{144}\cr 0&0&0&0&-\frac{233}{144}&-\frac{71}{144}&-\frac{43}{48}&-\frac{329}{144}&-\frac{203}{144}\cr 0&0&0&\frac{1}{2}&-\frac{529}{1296}&-\frac{151}{144}&-\frac{809}{1296}&-\frac{1937}{1296}&\frac{19}{48}\cr 0&0&0&-\frac{1}{2}&-\frac{1}{4}&\frac{1}{4}&-\frac{1}{4}&\frac{3}{4}&-1}\right] \)
取向量\( \left[ \matrix{\displaystyle H3[5,4]\cr H3[6,4]\cr H3[7,4]\cr H3[8,4]\cr H3[9,4]}\right]=\left[ \matrix{\displaystyle \frac{1}{2}\cr -\frac{1}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}}\right] \)
------------------------
\(k=4\)
\(X\)向量\(=\left[ \matrix{\displaystyle \frac{1}{2}\cr -\frac{1}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}}\right]\)第1元\(\displaystyle=\frac{1}{2},sign(X_1)=1\),長度\( \Vert\;X \Vert\;=sqrt(\left[ \matrix{\displaystyle \frac{1}{2}\cr -\frac{1}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}}\right].\left[ \matrix{\displaystyle \frac{1}{2}\cr -\frac{1}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}}\right] )=1\)
第1元是1,其餘為0的向量\( e1=\left[ \matrix{\displaystyle 1\cr 0\cr 0\cr 0\cr 0}\right] \)
\(X+sign(X_1)*\Vert\; X \Vert\;*e1=\left[ \matrix{\displaystyle \frac{1}{2}\cr -\frac{1}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}}\right]+1*1*\left[ \matrix{\displaystyle 1\cr 0\cr 0\cr 0\cr 0}\right]=\left[ \matrix{\displaystyle \frac{3}{2}\cr -\frac{1}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}}\right] \)
長度\(\Vert\; X+sign(X_1)*\Vert\; X \Vert\;*e1 \Vert\;=sqrt(\left[ \matrix{\displaystyle \frac{3}{2}\cr -\frac{1}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}}\right].\left[ \matrix{\displaystyle \frac{3}{2}\cr -\frac{1}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}}\right] )=\sqrt{3}\)
\(v=\left[ \matrix{\displaystyle \frac{3}{2}\cr -\frac{1}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}}\right]/\sqrt{3}=\left[ \matrix{\displaystyle \frac{\sqrt{3}}{2}\cr -\frac{1}{2\sqrt{3}}\cr 0\cr \frac{1}{2\sqrt{3}}\cr -\frac{1}{2\sqrt{3}}}\right] \)
\(\hat{P}=I-2*v.v^T=\left[ \matrix{\displaystyle 1&0&0&0&0\cr 0&1&0&0&0\cr 0&0&1&0&0\cr 0&0&0&1&0\cr 0&0&0&0&1}\right]-2*\left[ \matrix{\displaystyle \frac{\sqrt{3}}{2}\cr -\frac{1}{2\sqrt{3}}\cr 0\cr \frac{1}{2\sqrt{3}}\cr -\frac{1}{2\sqrt{3}}}\right].\left[ \matrix{\displaystyle \frac{\sqrt{3}}{2}&-\frac{1}{2\sqrt{3}}&0&\frac{1}{2\sqrt{3}}&-\frac{1}{2\sqrt{3}}}\right]=\left[ \matrix{\displaystyle -\frac{1}{2}&\frac{1}{2}&0&-\frac{1}{2}&\frac{1}{2}\cr \frac{1}{2}&\frac{5}{6}&0&\frac{1}{6}&-\frac{1}{6}\cr 0&0&1&0&0\cr -\frac{1}{2}&\frac{1}{6}&0&\frac{5}{6}&\frac{1}{6}\cr \frac{1}{2}&-\frac{1}{6}&0&\frac{1}{6}&\frac{5}{6}}\right] \)
取向量\( \left[ \matrix{\displaystyle \frac{1}{2}\cr -\frac{1}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}}\right] \)找一個\(Householder\)矩陣\( \hat{P}_4\),使得\(\left[ \matrix{\displaystyle -\frac{1}{2}&\frac{1}{2}&0&-\frac{1}{2}&\frac{1}{2}\cr \frac{1}{2}&\frac{5}{6}&0&\frac{1}{6}&-\frac{1}{6}\cr 0&0&1&0&0\cr -\frac{1}{2}&\frac{1}{6}&0&\frac{5}{6}&\frac{1}{6}\cr \frac{1}{2}&-\frac{1}{6}&0&\frac{1}{6}&\frac{5}{6}}\right] \left[ \matrix{\displaystyle \frac{1}{2}\cr -\frac{1}{2}\cr 0\cr \frac{1}{2}\cr -\frac{1}{2}}\right]=\left[ \matrix{\displaystyle -1\cr 0\cr 0\cr 0\cr 0}\right] \)
設\(P_4=\left[ \matrix{\displaystyle 1&0&0&0&0&0&0&0&0\cr 0&1&0&0&0&0&0&0&0\cr 0&0&1&0&0&0&0&0&0\cr 0&0&0&1&0&0&0&0&0\cr 0&0&0&0&-\frac{1}{2}&\frac{1}{2}&0&-\frac{1}{2}&\frac{1}{2}\cr 0&0&0&0&\frac{1}{2}&\frac{5}{6}&0&\frac{1}{6}&-\frac{1}{6}\cr 0&0&0&0&0&0&1&0&0\cr 0&0&0&0&-\frac{1}{2}&\frac{1}{6}&0&\frac{5}{6}&\frac{1}{6}\cr 0&0&0&0&\frac{1}{2}&-\frac{1}{6}&0&\frac{1}{6}&\frac{5}{6}}\right] \)
計算\(H_4=P_4^T.H_3.P_4=\left[ \matrix{\displaystyle -1&\frac{7}{3}&2&-\frac{3}{2}&\frac{13}{24}&-\frac{139}{216}&\frac{29}{36}&-\frac{101}{216}&-\frac{121}{216}\cr-2&-1&\frac{5}{3}&\frac{5}{12}&\frac{143}{144}&-\frac{931}{432}&-\frac{259}{72}&-\frac{485}{432}&\frac{2375}{432}\cr 0&1&-\frac{3}{2}&\frac{13}{12}&-\frac{59}{48}&-\frac{1321}{1296}&-\frac{625}{216}&-\frac{143}{1296}&\frac{6293}{1296}\cr 0&0&-1&-\frac{1}{4}&\frac{89}{48}&-\frac{2305}{1296}&-\frac{193}{54}&\frac{5305}{1296}&\frac{12653}{1296}\cr 0&0&0&-1&-\frac{1}{4}&\frac{3557}{3888}&-\frac{2545}{2592}&\frac{15283}{3888}&\frac{17111}{3888}\cr 0&0&0&0&-1&-\frac{679}{3888}&-\frac{1225}{2592}&\frac{8575}{3888}&\frac{10499}{3888}\cr 0&0&0&0&1&-\frac{295}{216}&-\frac{43}{48}&-\frac{305}{216}&-\frac{493}{216}\cr 0&0&0&0&1&-\frac{4297}{3888}&-\frac{3463}{2592}&\frac{3121}{3888}&\frac{11309}{3888}\cr 0&0&0&0&-1&\frac{113}{1296}&\frac{133}{288}&-\frac{1721}{1296}&-\frac{4837}{1296}}\right] \)
取向量\( \left[ \matrix{\displaystyle H4[6,5]\cr H4[7,5]\cr H4[8,5]\cr H4[9,5]}\right]=\left[ \matrix{\displaystyle -1\cr 1\cr 1\cr -1}\right] \)
------------------------
\(k=5\)
\(X\)向量\(=\left[ \matrix{-1\cr 1\cr 1\cr -1}\right] \)第1元\(=-1,sign(X_1)=-1\),長度\( \Vert\;X \Vert\;=sqrt(\left[ \matrix{-1\cr 1\cr 1\cr -1}\right].\left[ \matrix{-1\cr 1\cr 1\cr -1}\right] )=2\)
第1元是1,其餘為0的向量\( e1=\left[ \matrix{1\cr 0\cr 0\cr 0}\right] \)
\(X+sign(X_1)*\Vert\; X \Vert\;*e1=\left[ \matrix{-1\cr 1\cr 1\cr -1}\right]+-1*2*\left[ \matrix{\displaystyle 1\cr 0\cr 0\cr 0}\right]=\left[ \matrix{-3\cr 1\cr 1\cr -1}\right] \)
長度\(\Vert\; X+sign(X_1)*\Vert\; X \Vert\;*e1 \Vert\;=sqrt(\left[ \matrix{-3\cr 1\cr 1\cr -1}\right].\left[ \matrix{-3\cr 1\cr 1\cr -1}\right] )=2\sqrt{3}\)
\(v=\left[ \matrix{-3\cr 1\cr 1\cr -1}\right]/2\sqrt{3}=\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] \)
\(\hat{P}=I-2*v.v^T=\left[ \matrix{\displaystyle 1&0&0&0\cr 0&1&0&0\cr 0&0&1&0\cr 0&0&0&1}\right]-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].\left[ \matrix{\displaystyle -\frac{\sqrt{3}}{2}&\frac{1}{2\sqrt{3}}&\frac{1}{2\sqrt{3}}&-\frac{1}{2\sqrt{3}}}\right]=\left[ \matrix{\displaystyle -\frac{1}{2}&\frac{1}{2}&\frac{1}{2}&-\frac{1}{2}\cr \frac{1}{2}&\frac{5}{6}&-\frac{1}{6}&\frac{1}{6}\cr \frac{1}{2}&-\frac{1}{6}&\frac{5}{6}&\frac{1}{6}\cr -\frac{1}{2}&\frac{1}{6}&\frac{1}{6}&\frac{5}{6}}\right] \)
取向量\( \left[ \matrix{-1\cr 1\cr 1\cr -1}\right] \)找一個\(Householder\)矩陣\( \hat{P}_5\),使得\(\left[ \matrix{\displaystyle -\frac{1}{2}&\frac{1}{2}&\frac{1}{2}&-\frac{1}{2}\cr \frac{1}{2}&\frac{5}{6}&-\frac{1}{6}&\frac{1}{6}\cr \frac{1}{2}&-\frac{1}{6}&\frac{5}{6}&\frac{1}{6}\cr -\frac{1}{2}&\frac{1}{6}&\frac{1}{6}&\frac{5}{6}}\right] \left[ \matrix{-1\cr 1\cr 1\cr -1}\right]=\left[ \matrix{\displaystyle 2\cr 0\cr 0\cr 0}\right] \)
設\(P_5=\left[ \matrix{\displaystyle 1&0&0&0&0&0&0&0&0\cr 0&1&0&0&0&0&0&0&0\cr 0&0&1&0&0&0&0&0&0\cr 0&0&0&1&0&0&0&0&0\cr 0&0&0&0&1&0&0&0&0\cr 0&0&0&0&0&-\frac{1}{2}&\frac{1}{2}&\frac{1}{2}&-\frac{1}{2}\cr 0&0&0&0&0&\frac{1}{2}&\frac{5}{6}&-\frac{1}{6}&\frac{1}{6}\cr 0&0&0&0&0&\frac{1}{2}&-\frac{1}{6}&\frac{5}{6}&\frac{1}{6}\cr 0&0&0&0&0&-\frac{1}{2}&\frac{1}{6}&\frac{1}{6}&\frac{5}{6}}\right] \)
計算\(H_5=P_5^T.H_4.P_5=\left[ \matrix{\displaystyle -1&\frac{7}{3}&2&-\frac{3}{2}&\frac{13}{24}&\frac{37}{48}&\frac{433}{1296}&-\frac{1217}{1296}&-\frac{115}{1296}\cr -2&-1&\frac{5}{3}&\frac{5}{12}&\frac{143}{144}&-\frac{129}{32}&-\frac{7703}{2592}&-\frac{1289}{2592}&\frac{12629}{2592}\cr 0&1&-\frac{3}{2}&\frac{13}{12}&-\frac{59}{48}&-\frac{985}{288}&-\frac{16277}{7776}&\frac{5365}{7776}&\frac{31535}{7776}\cr 0&0&-1&-\frac{1}{4}&\frac{89}{48}&-\frac{1075}{288}&-\frac{22727}{7776}&\frac{36895}{7776}&\frac{70853}{7776}\cr 0&0&0&-1&-\frac{1}{4}&-\frac{2045}{1728}&-\frac{13177}{46656}&\frac{216029}{46656}&\frac{172699}{46656}\cr 0&0&0&0&2&-\frac{859}{1152}&-\frac{39167}{31104}&-\frac{27749}{31104}&\frac{30413}{31104}\cr 0&0&0&0&0&\frac{8155}{10368}&-\frac{360449}{279936}&-\frac{289307}{279936}&-\frac{295117}{279936}\cr 0&0&0&0&0&-\frac{10907}{10368}&-\frac{287807}{279936}&\frac{526555}{279936}&\frac{962189}{279936}\cr 0&0&0&0&0&\frac{13199}{10368}&-\frac{72541}{279936}&-\frac{789199}{279936}&-\frac{1077113}{279936}}\right] \)
取向量\( \left[ \matrix{\displaystyle H5[7,6]\cr H5[8,6]\cr H5[9,6]}\right]=\left[ \matrix{\displaystyle \frac{8155}{10368}\cr -\frac{10907}{10368}\cr \frac{13199}{10368}}\right] \)
------------------------
向量\(X\)長度\(\displaystyle=\frac{5\sqrt{1598579}}{3456}\)不是完全平方數,改用浮點數計算
\(k=6\)
\(X\)向量\(=\left[ \matrix{\displaystyle \frac{8155}{10368}\cr -\frac{10907}{10368}\cr \frac{13199}{10368}}\right]\)第1元\(=0.7865547839506173,sign(X_1)=1\),長度\( \Vert\;X \Vert\;=sqrt(\left[ \matrix{\displaystyle \frac{8155}{10368}\cr -\frac{10907}{10368}\cr \frac{13199}{10368}}\right].\left[ \matrix{\displaystyle \frac{8155}{10368}\cr -\frac{10907}{10368}\cr \frac{13199}{10368}}\right] )=1.829208969513448\)
第1元是1,其餘為0的向量\( e1=\left[ \matrix{1\cr 0\cr 0}\right] \)
\(X+sign(X_1)*\Vert\; X \Vert\;*e1=\left[ \matrix{\displaystyle \frac{8155}{10368}\cr -\frac{10907}{10368}\cr \frac{13199}{10368}}\right]+1*1.829208969513448*\left[ \matrix{\displaystyle 1\cr 0\cr 0}\right]=\left[ \matrix{\displaystyle 2.615763753464065\cr -1.05198688271605\cr 1.273051697530864}\right] \)
長度\(\Vert\; X+sign(X_1)*\Vert\; X \Vert\;*e1 \Vert\;=sqrt(\left[ \matrix{2.615763753464065\cr -1.05198688271605\cr 1.273051697530864}\right].\left[ \matrix{2.615763753464065\cr -1.05198688271605\cr 1.273051697530864}\right] )=3.093470064495413 \)
\(v=\left[ \matrix{2.615763753464065\cr -1.05198688271605\cr 1.273051697530864}\right]/3.093470064495413=\left[ \matrix{0.8455759063214779\cr -0.3400669347959709\cr 0.4115286946339066}\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{0.8455759063214779\cr -0.3400669347959709\cr 0.4115286946339066}\right].\left[ \matrix{ 0.8455759063214779&-0.3400669347959709&0.4115286946339066}\right]=\)
\(\left[ \matrix{ -0.4299972267027781&0.5751048132001405&-0.6959574978847209\cr 0.5751048132001405&0.7687089597169466&0.2798946035294795\cr -0.6959574978847209&0.2798946035294795&0.661288266985825}\right] \)
取向量\( \left[ \matrix{\displaystyle \frac{8155}{10368}\cr -\frac{10907}{10368}\cr \frac{13199}{10368}}\right] \)找一個\(Householder\)矩陣\( \hat{P}_6\),使得\(\left[ \matrix{ -0.4299972267027781&0.5751048132001405&-0.6959574978847209\cr 0.5751048132001405&0.7687089597169466&0.2798946035294795\cr -0.6959574978847209&0.2798946035294795&0.661288266985825}\right] \left[ \matrix{\displaystyle \frac{8155}{10368}\cr -\frac{10907}{10368}\cr \frac{13199}{10368}}\right]=\left[ \matrix{-1.829208969513457\cr 7.771561172376096*10^{-16}\cr -2.886579864025407*10^{-15}}\right] \)
設\(P_6=\left[ \matrix{1&0&0&0&0&0&0&0&0\cr 0&1&0&0&0&0&0&0&0\cr 0&0&1&0&0&0&0&0&0\cr 0&0&0&1&0&0&0&0&0\cr 0&0&0&0&1&0&0&0&0\cr 0&0&0&0&0&1&0&0&0\cr 0&0&0&0&0&0&-0.4299972267027781&0.5751048132001405&-0.6959574978847209\cr 0&0&0&0&0&0&0.5751048132001405&0.7687089597169466&0.2798946035294795\cr 0&0&0&0&0&0&-0.6959574978847209&0.2798946035294795&0.661288266985825}\right] \)
計算\(H_6=P_6^T.H_5.P_6=\left[ \matrix{ -1&2.333333333333333&2&-1.5&0.5416666666666666&0.7708333333333334&-0.6219569788349768&-0.5545418975815991&-0.5540350924250236\cr -2&-1&1.666666666666667&0.4166666666666667&0.9930555555555556&-4.03125&-2.39903113723365&-0.7276656161967713&5.151074840285687\cr 0&1&-1.5&1.083333333333333&-1.229166666666667&-3.420138888888889&-1.525529515680916&0.4616279381076218&4.331720709928399\cr 0&0&-1&-0.25&1.854166666666667&-3.732638888888889&-2.355920465963594&4.516780134134694&9.387602241820074\cr 0&0&0&-1&-0.25&-1.183449074074074&0.2082096026209324&4.4329301669687&3.940325053062458\cr 0&0&0&0&2&-0.7456597222222222&-0.6521038280089702&-1.136304030499402&1.273261703644264\cr 0&0&0&0&0&-1.829208969513453&-1.610172351137944&4.19317878751438&4.254779928150215\cr 0&0&0&0&0&7.7715611723761*10^{-16}&0.0579988841413669&-0.605664190520036&1.767500658619938\cr 0&0&0&0&0&-2.886579864025409*10^{-15}&0.04700389329356565&-0.4617357301648267&-1.038503736119814}\right] \)
取向量\( \left[ \matrix{H6[8,7]\cr H6[9,7]}\right]=\left[ \matrix{\displaystyle 0.0579988841413669\cr 0.04700389329356565}\right] \)
------------------------
\(k=7\)
\(X\)向量\(=\left[ \matrix{\displaystyle 0.0579988841413669\cr 0.04700389329356565}\right]\)第1元\(=0.0579988841413669,sign(X_1)=1\),長度\( \Vert\;X \Vert\;=sqrt(\left[ \matrix{0.0579988841413669\cr 0.04700389329356565}\right].\left[ \matrix{0.0579988841413669\cr 0.04700389329356565}\right] )=0.07465411272258612 \)
第1元是1,其餘為0的向量\( e1=\left[ \matrix{1\cr 0}\right] \)
\(X+sign(X_1)*\Vert\; X \Vert\;*e1=\left[ \matrix{0.0579988841413669\cr 0.04700389329356565}\right]+1*0.07465411272258612*\left[ \matrix{1\cr 0}\right]=\left[ \matrix{0.132652996863953\cr 0.04700389329356565}\right] \)
長度\(\Vert\; X+sign(X_1)*\Vert\; X \Vert\;*e1 \Vert\;=sqrt(\left[ \matrix{0.132652996863953\cr 0.04700389329356565}\right].\left[ \matrix{0.132652996863953\cr 0.04700389329356565}\right] )=0.1407344434093547 \)
\(v=\left[ \matrix{\displaystyle 0.132652996863953\cr 0.04700389329356565}\right]/0.1407344434093547=\left[ \matrix{0.9425766262357311\cr 0.3339899754095393}\right] \)
\(\hat{P}=I-2*v.v^T=\left[ \matrix{\displaystyle 1&0\cr 0&1}\right]-2*\left[ \matrix{0.9425766262357311\cr 0.3339899754095393}\right].\left[ \matrix{0.9425766262357311&0.3339899754095393}\right]=\left[\matrix{-0.7769013926518663&-0.6296222884361566\cr -0.6296222884361566&0.7769013926518692}\right] \)
取向量\( \left[ \matrix{0.0579988841413669\cr 0.04700389329356565}\right]\)找一個\(Householder\)矩陣\( \hat{P}_7\),使得\(\left[ \matrix{-0.7769013926518663&-0.6296222884361566\cr -0.6296222884361566&0.7769013926518692}\right] \left[ \matrix{0.0579988841413669\cr 0.04700389329356565}\right]=\left[ \matrix{-0.0746541127225859\cr 6.245004513516506*10^{-17}}\right] \)
設\(P_7=\left[ \matrix{1&0&0&0&0&0&0&0&0\cr 0&1&0&0&0&0&0&0&0\cr 0&0&1&0&0&0&0&0&0\cr 0&0&0&1&0&0&0&0&0\cr 0&0&0&0&1&0&0&0&0\cr 0&0&0&0&0&1&0&0&0\cr 0&0&0&0&0&0&1&0&0\cr 0&0&0&0&0&0&0&-0.7769013926518663&-0.6296222884361566\cr 0&0&0&0&0&0&0&-0.6296222884361566&0.7769013926518692}\right] \)
計算\(H_7=P_7^T.H_6.P_7=\left[ \matrix{-1&2.333333333333333&2&-1.5&0.5416666666666666&0.7708333333333334&-0.6219569788349768&0.7796572152815339&-0.08127869629395262\cr -2&-1&1.666666666666667&0.4166666666666667&0.9930555555555556&-4.03125&-2.39903113723365&-2.677907098238435&4.460031707558076\cr 0&1&-1.5&1.083333333333333&-1.229166666666667&-3.420138888888889&-1.525529515680916&-3.085987294254227&3.074668613324933\cr 0&0&-1&-0.25&1.854166666666667&-3.732638888888889&-2.355920465963594&-9.419736382934685&4.449375810914967\cr 0&0&0&-1&-0.25&-1.183449074074074&0.2082096026209324&-5.924866097337959&0.2701723850207669\cr 0&0&0&0&2&-0.7456597222222222&-0.6521038280089702&0.08112223614429317&1.704641134813782\cr 0&0&0&0&0&-1.829208969513453&-1.610172351137944&-5.936590714812356&0.6654256275903683\cr 0&0&0&0&0&1.213681349951592*10^{-15}&-0.07465411272258583&-0.1385315695498772&-1.038138425003079\cr 0&0&0&0&0&-2.731902729369484*10^{-15}&6.245004513516506*10^{-17}&1.191097963781674&-1.505636357089966}\right] \)
[color=red](%o7)[/color]　\( \left[\matrix{-1&2.333333333333333&2&-1.5&0.5416666666666666&0.7708333333333334&-0.6219569788349768&0.7796572152815339&-0.08127869629395262\cr -2&-1&1.666666666666667&0.4166666666666667&0.9930555555555556&-4.03125&-2.39903113723365&-2.677907098238435&4.460031707558076\cr 0&1&-1.5&1.083333333333333&-1.229166666666667&-3.420138888888889&-1.525529515680916&-3.085987294254227&3.074668613324933\cr 0&0&-1&-0.25&1.854166666666667&-3.732638888888889&-2.355920465963594&-9.419736382934685&4.449375810914967\cr 0&0&0&-1&-0.25&-1.183449074074074&0.2082096026209324&-5.924866097337959&0.2701723850207669\cr 0&0&0&0&2&-0.7456597222222222&-0.6521038280089702&0.08112223614429317&1.704641134813782\cr 0&0&0&0&0&-1.829208969513453&-1.610172351137944&-5.936590714812356&0.6654256275903683\cr 0&0&0&0&0&1.213681349951592*10^{-15}&-0.07465411272258583&-0.1385315695498772&-1.038138425003079\cr 0&0&0&0&0&-2.731902729369484*10^{-15}&6.245004513516506*10^{-17}&1.191097963781674&-1.505636357089966}\right]\)

bugmens 發表於 2019-8-3 04:45

上一篇文章找\(H\)右下角\(4 \times 4\)子矩陣\(T=\left[ \matrix{\times & \times & \times & \times \cr \times & \times & \times & \times \cr 0 & \times & \times & \times \cr 0 & 0 & \times & \times} \right]\)
計算\(det(T-\mu I)=0\)，得到4個特徵值\(\mu=\mu_1,\mu_2,\mu_3,\mu_4\)
計算\((H-\mu_1 I)(H-\mu_2 I)(H-\mu_3 I)(H-\mu_4 I)e1=\left[\matrix{\times \cr \times \cr \times \cr \times \cr \times \cr 0 \cr 0 \cr 0 \cr 0} \right]\)
其中\(e1\)為第一元為1，其餘為0的向量，換句話說就是取\((H−\mu_1I)(H−\mu_2I)(H−\mu_3I)(H−\mu_4I)\)相乘後的第1行

但實際上不需要計算子矩陣\(T\)的特徵值，可利用遞迴式求得偏移向量。[table=98%][tr][td]論文遞迴式[/td][td]符號解釋[/td][/tr]
[tr][td]\(w_1=e_1\)，\(z_1=h_1-t_{11}e_1\)，
\(w_i=t_{ii-1}^{-1}z_{i-1}\)
\(\displaystyle z_i=Hw_i-\sum_{j=1}^i t_{ji}w_j\)　　\(i=2,\ldots,m.\)[/td][td]\(e_1\)為第一元為1其餘為0的向量
\(h_1\)為\(H\)矩陣的第一行
\(t_{ii-1}\)為子矩陣\(T\)第\(i\)列,第\(i-1\)行數字
\(m\)為子矩陣\(T\)的大小[/td][/tr][/table]最後得到的\(z_m\)就是\((H-\mu_1 I)(H-\mu_2 I)(H-\mu_3 I)(H-\mu_4 I)\)相乘後第一行，論文中稱為\(shift\) \(vector\)。

參考資料
A. Dubrulle, G. Golub, A multishift QR iteration without computation of the shifts, Numer. Algorithms 7 (1994) 173-181.
[url]https://www.hpl.hp.com/techreports/93/HPL-93-54.pdf[/url]



[color=green]使用遞迴式得到(H-u1I)(H-u2I)...(H-ukI)第一行[/color]
[color=red](%i1)[/color]
[color=blue]theFirstColumnofΠk(H,Bulge):=block
([n:length(H),T,e1,W,Z,X],
powerdisp : true,/*改變變數顯示順序*/
print("取H矩陣右下角",Bulge,"*",Bulge,"子矩陣T=",
       T:genmatrix(lambda([i,j],H[n-Bulge+i,n-Bulge+j]),Bulge,Bulge)),
e1:genmatrix(lambda([i,j],if i=1 then 1 else 0),n,1),
W:create_list(0,i,1,Bulge),
Z:create_list(0,i,1,Bulge),
print("第1元為1,其餘為0的向量",w[1],"=",e[1],"=",W[1]:e1),
print(z[1],"=",h[1],"-",t[1,1],e[1],"=",col(H,1),"-",T[1,1],e1,"=",Z[1]:col(H,1)-T[1,1]*e1),
for i:2 thru Bulge do
    (print(w[ i ],"=1/",t[i,i-1],z[i-1],"=1/",T[i,i-1],Z[i-1],"=",W[ i ]:1/T[i,i-1]*Z[i-1]),
     print(z[ i ],"=H",w[ i ],"-(",apply("+",create_list(t[j,i]*w[j],j,1,i)),")"),
     printList:[z[ i ],"=",H.W[ i ],"-("],/*先將要印的東西放在printList*/
     for j:1 thru i do
        (if T[j,i]>=0 then
            (printList:append(printList,["+",T[j,i],W[j]]))/*T[j,i]是正的,加上+號*/
         else
            (printList:append(printList,[T[j,i],W[j]]))/*T[j,i]是負的,原本就有-號,不需要再加*/
        ),
     apply(print,append(printList,[")=",Z[ i ]:H.W[ i ]-sum(T[j,i]*W[j],j,1,i)]))/*再用apply(print,)將全部內容印在同一行*/
    ),
print("擷取",z[Bulge],"前",Bulge+1,"元數字X=",X:genmatrix(lambda([i,j],Z[Bulge][i,j]),Bulge+1,1)),
return(X)
)$[/color]

[color=green]Hessenberg矩陣[/color]
[color=red](%i2)[/color]
[color=blue]H:matrix([-2,1,-1, 1,2,-2/3,-3,0,29/6],
             [-1,-1,1,2,-2,1,2/3,-5/2,-23/6],
             [0,-1,1,2,-1,1,6,-8/3,-65/6],
             [0,0,-1,-2,0,3,-1,7/2,10/3],
             [0,0,0,-1,-1,2,-1,-1,-4],
             [0,0,0,0,-2,1,2,2,-1],
             [0,0,0,0,0,-1,-2,0,1],
             [0,0,0,0,0,0,1,-1,-2],
             [0,0,0,0,0,0,0,1,-1] );[/color]
[color=red](H)[/color]　\( \left[ \matrix{\displaystyle -2&1&-1&1&2&-\frac{2}{3}&-3&0&\frac{29}{6}\cr
-1&-1&1&2&-2&1&\frac{2}{3}&-\frac{5}{2}&-\frac{23}{6}\cr
0&-1&1&2&-1&1&6&-\frac{8}{3}&-\frac{65}{6}\cr
0&0&-1&-2&0&3&-1&\frac{7}{2}&\frac{10}{3}\cr
0&0&0&-1&-1&2&-1&-1&-4\cr
0&0&0&0&-2&1&2&2&-1\cr
0&0&0&0&0&-1&-2&0&1\cr
0&0&0&0&0&0&1&-1&-2\cr
0&0&0&0&0&0&0&1&-1} \right] \)

[color=green]取H右下角4*4子矩陣，使用遞迴式得到(H-u1I)(H-u2I)...(H-u4I)第一行[/color]
[color=red](%i3)[/color]　[color=blue]X:theFirstColumnofΠk(H,4);[/color]
取\(H\)矩陣右下角\(4*4\)子矩陣\(T=\left[\matrix{1&2&2&-1\cr -1&-2&0&1\cr 0&1&-1&-2\cr 0&0&1&-1}\right]\)
第1元為1,其餘為0的向量\(w_1=e_1=\left[\matrix{1\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0}\right]\)
\(z_1=h_1-t_{1,1}e_1=\left[\matrix{-2\cr -1\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0}\right]-1 \left[\matrix{1\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0}\right]=\left[\matrix{-3\cr -1\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0}\right]\)
\(w_2=1/t_{2,1}z_1=1/-1\left[\matrix{-3\cr -1\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0}\right]=\left[\matrix{3\cr 1\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0}\right]\)
\(z_2=Hw_2-(w_1t_{1,2}+w_2t_{2,2})\)
\(z_2=\left[\matrix{-5\cr -4\cr -1\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0}\right]-(+2\left[\matrix{1\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0}\right]-2\left[\matrix{3\cr 1\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0}\right])=\left[\matrix{-1\cr -2\cr -1\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0}\right]\)
\(w_3=1/t_{3,2}z_2=1/1 \left[\matrix{-1\cr -2\cr -1\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0}\right]=\left[\matrix{-1\cr -2\cr -1\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0}\right]\)
\(z_3=Hw_3-(w_1 t_{1,3}+w_2 t_{2,3}+w_3 t_{3,3})\)
\(z_3=\left[\matrix{1\cr 2\cr 1\cr 1\cr 0\cr 0\cr 0\cr 0\cr 0}\right]-(+2\left[\matrix{1\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0}\right]+0\left[\matrix{3\cr 1\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0}\right]-1\left[\matrix{-1\cr -2\cr -1\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0}\right])=\left[\matrix{-2\cr 0\cr 0\cr 1\cr 0\cr 0\cr 0\cr 0\cr 0}\right]\)
\(w_4=1/t_{4,3}z_3=1/1 \left[\matrix{-2\cr 0\cr 0\cr 1\cr 0\cr 0\cr 0\cr 0\cr 0}\right]=\left[\matrix{-2\cr 0\cr 0\cr 1\cr 0\cr 0\cr 0\cr 0\cr 0}\right]\)
\(z_4=Hw_4-(w_1 t_{1,4}+w_2 t_{2,4}+w_3 t_{3,4}+w_4 t_{4,4})\)
\(z_4=\left[\matrix{5\cr 4\cr 2\cr -2\cr -1\cr 0\cr 0\cr 0\cr 0}\right]-(-1\left[\matrix{1\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0}\right]+1\left[\matrix{3\cr 1\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0}\right]-2\left[\matrix{-1\cr -2\cr -1\cr 0\cr 0\cr 0\cr 0\cr 0\cr 0}\right]-1\left[\matrix{-2\cr 0\cr 0\cr 1\cr 0\cr 0\cr 0\cr 0\cr 0}\right])=\left[\matrix{-1\cr -1\cr 0\cr -1\cr -1\cr 0\cr 0\cr 0\cr 0}\right]\)
擷取\(z_4\)前5元數字\(X=\left[\matrix{-1\cr -1\cr 0\cr -1\cr -1}\right]\)
[color=red](X)[/color]　\(\left[\matrix{-1\cr -1\cr 0\cr -1\cr -1}\right]\)

查看完整版本: 用maxima學數值分析-特徵值和特徵向量