21 123
發新話題
打印

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

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


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

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


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


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


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


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


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



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

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

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


具有WilkinsonShift的QR法副程式
(%i4)
QRwithWilkinsonShift(A,epsilon):=block
([count:0,Q,R,I,n:length(A),subA,eig1,eig2,eigenvalue],
eig2:makelist(0,n),
I:ident(n),
do
  (count:count+1,
   print("第",count,"輪"),
   A:SymmetricWilkinsonShift(A),
   eig1:create_list(A[i,i],i,1,n),
   print("特徵值",eig1),
   if lmax(abs(eig1-eig2))<epsilon then
     (print("執行次數",count),
     return(eig1)
     ),
   eig2:eig1,
   print("-----------")
  )
)$


對稱矩陣
(%i5)
A:matrix([3599/900,-1501/900,-1201/900,-901/900,151/450,-301/900],
[-1501/900,3299/900,-901/900,-601/900,301/450,-1/900],
[-1201/900,-901/900,2999/900,-301/900,451/450,299/900],
[-901/900,-601/900,-301/900,2699/900,601/450,599/900],
[151/450,301/450,451/450,601/450,374/225,901/450],
[-301/900,-1/900,299/900,599/900,901/450,2099/900]);

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

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

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

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

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

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

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

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

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

特徵值\( [-1.01,6.0,5.0,3.0,0.999997,4.0] \)
執行次數13
(%o6) \( [-1.01,6.0,5.0,3.0,0.999997,4.0] \)

三對角矩陣
(%i7)
A:matrix([3.99889,2.40601,0,0,0,0],
   [2.40601,1.68454,-1.93868,0,0,0],
   [0,-1.93868,2.37077,-1.67177,0,0],
   [0,0,-1.67177,2.35489,1.31944,0],
   [0,0,0,1.31944,3.60161,-0.72545],
   [0,0,0,0,-0.72545,3.97930]);

(%o7) \( \left[ \matrix{3.99889&2.40601&0&0&0&0 \cr
2.40601&1.68454&-1.93868&0&0&0 \cr
0&-1.93868&2.37077&-1.67177&0&0 \cr
0&0&-1.67177&2.35489&1.31944&0 \cr
0&0&0&1.31944&3.60161&-0.725455 \cr
0&0&0&0&-0.725455&        3.9793} \right] \)

計算矩陣A全部的特徵值,誤差10^-5
(%i8) QRwithWilkinsonShift(A,10^-5);
第1輪
\(A\)的右下角2*2矩陣\( \left[ \matrix{
3.60161 & -0.72545\cr
-0.72545 & 3.9793} \right] \)

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

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

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

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

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

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

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

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

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

TOP

 21 123
發新話題