|
|
23#
大 中
小 發表於 2017-11-4 09:54 只看該作者
3-4.減少QR法運算-Hessenberg矩陣的Implicit Double Shift
前面文章提到 https://math.pro/db/viewthread.php?tid=2561&page=2#pid17755
取矩陣 A 的右下角 2 \times 2 子矩陣 \left[ \matrix{A[n-1,n-1]& A[n-1,n]\cr A[n,n-1]&A[n,n]} \right] ，
計算特徵方程式 \left| \matrix{A[n-1,n-1]-x& A[n-1,n] \cr A[n,n-1]&A[n,n]-x} \right|=0 ，
x^2-(A[n-1,n-1]+A[n,n])x+(A[n-1,n-1]A[n,n]-A[n-1,n]A[n,n-1])=0
判別式 D=(A[n-1,n-1]+A[n,n])^2-4(A[n-1,n-1]A[n,n]-A[n-1,n]A[n,n-1])
D=(A[n-1,n-1]-A[n,n])^2+4A[n-1,n]A[n,n-1]
D\ge 0 特徵值為實數， D<0 特徵值為複數
(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 矩陣
}
－－－－－－－－－－－－－－－－－－－－－－－－－－－
參考資料 https://books.google.com.tw/book ... e&q&f=false
要先載入diag才能使用diag指令
(%i1) load(diag);
(%o1) C:\maxima-5.40.0\share\maxima\5.40.0\share\contrib\diag.mac
轉換成Householder矩陣副程式
(%i2)
Householder(X):=block
([n,e1,normx,v,normv,P_hat],
n:length(X),
if X=zeromatrix(n,1) then
(print("X為0向量,H=",ident(n)),
return(ident(n))
),
if X[1,1]>=0 then sign:+1 else sign:-1,
print("X向量=",X,"第1元=",X[1,1]," , sign(X_1)=",sign,",長度∥X∥=sqrt(",X,".",X,")=",normx:ratsimp(sqrt(X.X))),
print("第1元是1,其餘為0的向量e1=",e1:genmatrix(lambda([i,j],if i=1 then 1 else 0),n,1)),
print("X+sign(X_1)*∥X∥*e1=",X,"+",sign,"*",normx,"*",e1,"=",v:X+sign*normx*e1),
print("長度∥X+sign(X_1)*∥X∥*e1∥=sqrt(",v,".",v,")=",normv:ratsimp(sqrt(v.v))),
print("v=",v,"/",normv,"=",v:ratsimp(v/normv)),
print("P_hat=I-2*v.v^T=",ident(n),"-2*",v,".",transpose(v),"=",P_hat:ratsimp(ident(n)-2*v.transpose(v))),
return(P_hat)
)$
針對Hessenberg矩陣Implicit Double Shift副程式
(%i3)
HessenbergImplicitDoubleShift(H):=block
([n,t,d,x,y,z,k,X,P_hat_k],
n:length(H),
print("t=H[",n-1,",",n-1,"]+H[",n,",",n,"]=",H[n-1,n-1],"+",H[n,n],"=",t:H[n-1,n-1]+H[n,n]),
print("d=H[",n-1,",",n-1,"]*H[",n,",",n,"]-H[",n,",",n-1,"]*H[",n-1,",",n,"]=",
H[n-1,n-1],"*",H[n,n],"-",H[n,n-1],"*",H[n-1,n],"=",d:H[n-1,n-1]*H[n,n]-H[n,n-1]*H[n-1,n]),
print("x=H[1,1]^2-t*H[1,1]+d+H[1,2]*H[2,1]=("
,H[1,1],")^2-",t,"*",H[1,1],"+",d,"+",H[1,2],"*",H[2,1],"=",x:H[1,1]^2-t*H[1,1]+d+H[1,2]*H[2,1]),
print("y=H[2,1]*(H[1,1]+H[2,2]-t)=",H[2,1],"*(",H[1,1],"+",H[2,2],"-",t,")=",y:H[2,1]*(H[1,1]+H[2,2]-t)),
print("z=H[2,1]*H[3,2]=",H[2,1],"*",H[3,2],"=",z:H[2,1]*H[3,2]),
print("取向量",matrix(["x"],["y"],["z"]),"=",matrix([x],[y],[z])),
print("------------------------"),
for k:0 thru n-3 do
(print("k=",k),
print("取向量",X:matrix([x],[y],[z]),"找一個3*3Householder矩陣P_hat_",k,",使得",P_hat_k:Householder(X),X,"=",P_hat_k.X),
if k=0 then
(P_k:diag([P_hat_k,ident(n-k-3)]))
else if k=n-3 then
(P_k:diag([ident(n-3),P_hat_k]))
else
(P_k:diag([ident(k),P_hat_k,ident(n-k-3)])),
print("設P_",k,"=",P_k,",計算H_",k,"=P_",k,"^T.H_",k-1,".P_",k,"=",H:ratsimp(transpose(P_k).H.P_k)),
x:H[k+2,k+1],
y:H[k+3,k+1],
if k<n-3 then
(z:H[k+4,k+1],
print("取向量",matrix([concat("H",k,"[",k+2,",",k+1,"]")],[concat("H",k,"[",k+3,",",k+1,"]")],[concat("H",k,"[",k+4,",",k+1,"]")]),"=",matrix([x],[y],[z]))
)
else
(print("取向量",matrix([concat("H",k,"[",k+2,",",k+1,"]")],[concat("H",k,"[",k+3,",",k+1,"]")]),"=",matrix([x],[y]))
),
print("------------------------")
),
k:n-2,
print("取向量",X:matrix([x],[y]),"找一個2*2Householder矩陣P_hat_",k,",使得",P_hat_k:Householder(X),X,"=",P_hat_k.X),
print("設P_",k,"=",P_k:diag([ident(n-2),P_hat_k]),",計算H",k,"=P_",k,"^T.H_",k-1,".P_",k,"=",H:ratsimp(transpose(P_k).H.P_k))
)$
Hessenberg矩陣，各元數字經過設計， \left[ \matrix{x \cr y \cr z} \right] 向量長度均為完全平方數，處理後的 Hessenberg 矩陣各元最多是分數，不會出現根號數字
(%i4)
H:matrix([-1,-1,2/5,-1/4,11/8,31/80],
[1,-2,2,11/5,-11/5,-21/5],
[0,-2,1,2,2,-10],
[0,0,3,-3,0,10],
[0,0,0,3,-2,-2],
[0,0,0,0,1,-2]);
(%o4) \left[ \matrix{\displaystyle
-1& -1 & \frac{2}{5} & -\frac{1}{4} & \frac{11}{8} & \frac{31}{80} \cr
1 & -2 & 2 & \frac{11}{5} & -\frac{11}{5} & -\frac{21}{5} \cr
0& -2 & 1 & 2 & 2 & -10 \cr
0& 0 & 3 & -3 & 0 & 10 \cr
0& 0 & 0 & 3 & -2 & -2 \cr
0& 0 & 0 & 0 & 1 & -2} \right]
執行Hessenberg矩陣的Implicit Double Shift
(%i5) HessenbergImplicitDoubleShift(H);
t=H[5,5]+H[6,6]=-2+-2=-4
d=H[5,5]*H[6,6]-H[6,5]*H[5,6]=-2*-2-1*-2=6
x=H[1,1]^2-t*H[1,1]+d+H[1,2]*H[2,1]=(-1)^2--4*-1+6+-1*1=2
y=H[2,1]*(H[1,1]+H[2,2]-t)=1*(-1+-2--4 )=1
z=H[2,1]*H[3,2]=1*-2=-2
取向量 \left[ \matrix{x \cr y \cr z} \right]=\left[ \matrix{2 \cr 1 \cr -2} \right]
------------------------
k=0
X 向量 =\left[ \matrix{2 \cr 1 \cr -2} \right] 第1元 =2, sign(X_1)=1,長度 \Vert\; X \Vert\;=sqrt(\left[ \matrix{2 \cr 1 \cr -2} \right] \cdot \left[ \matrix{2 \cr 1 \cr -2} \right])=3
第1元是1,其餘為0的向量 e1=\left[ \matrix{1 \cr 0 \cr 0} \right]
X+sign(X_1)*\Vert\; X \Vert\;*e1=\left[ \matrix{2 \cr 1 \cr -2} \right]+1*3*\left[ \matrix{1 \cr 0 \cr 0} \right]=\left[ \matrix{5 \cr 1 \cr -2} \right]
長度 \Vert\; X+sign(X_1)*\Vert\; X \Vert\;*e1 \Vert\;=sqrt(\left[ \matrix{5 \cr 1 \cr -2} \right] \cdot \left[ \matrix{5 \cr 1 \cr -2} \right])=\sqrt{30}
v=\left[ \matrix{5 \cr 1 \cr -2} \right]/ \sqrt{30}=\left[ \matrix{\displaystyle \frac{5}{\sqrt{30}}\cr \frac{1}{\sqrt{30}}\cr -\frac{2}{\sqrt{30}}} \right]
\hat{P}=I-2*v.v^T=\left[ \matrix{1&0&0 \cr 0&1&0 \cr 0&0&1} \right]-2*\left[ \matrix{\displaystyle \frac{5}{\sqrt{30}} \cr \frac{1}{\sqrt{30}} \cr -\frac{2}{\sqrt{30}}} \right] \cdot \left[ \matrix{\displaystyle \frac{5}{\sqrt{30}}& \frac{1}{\sqrt{30}}& -\frac{2}{\sqrt{30}}} \right]=\left[ \matrix{\displaystyle -\frac{2}{3}&-\frac{1}{3}&\frac{2}{3} \cr
-\frac{1}{3}&\frac{14}{15}&\frac{2}{15} \cr
\frac{2}{3}&\frac{2}{15}&\frac{11}{15}} \right]
取向量 \left[ \matrix{2 \cr 1 \cr -2} \right] 找一個3*3 Householder 矩陣 \hat{P_0} ,使得 \left[ \matrix{\displaystyle -\frac{2}{3}&-\frac{1}{3}&\frac{2}{3} \cr
-\frac{1}{3}&\frac{14}{15}&\frac{2}{15} \cr
\frac{2}{3}&\frac{2}{15}&\frac{11}{15}} \right] \left[ \matrix{2 \cr 1 \cr -2} \right]=\left[ \matrix{-3 \cr 0 \cr 0} \right]
設 P_0=\left[ \matrix{\displaystyle -\frac{2}{3}&-\frac{1}{3}&\frac{2}{3}&0&0&0 \cr
-\frac{1}{3}&\frac{14}{15}&\frac{2}{15}&0&0&0 \cr
\frac{2}{3}&\frac{2}{15}&\frac{11}{15}&0&0&0 \cr
0&0&0&1&0&0 \cr
0&0&0&0&1&0 \cr
0&0&0&0&0&1 } \right] ,計算 H_0=P_0^T.H_{-1}.P_0=\left[ \matrix{\displaystyle
-\frac{2}{5}&-\frac{11}{75}&\frac{2}{75}&\frac{23}{30}&\frac{23}{20}&-\frac{221}{40} \cr
1&-\frac{139}{75}&\frac{148}{75}&\frac{721}{300}&-\frac{449}{200}&-\frac{2153}{400} \cr
2&-\frac{142}{75}&\frac{19}{75}&\frac{239}{150}&\frac{209}{100}&-\frac{1527}{200} \cr
2&\frac{2}{5}&\frac{11}{5}&-3&0&10 \cr
0&0&0&3&-2&-2 \cr
0&0&0&0&1&-2 \cr} \right]
取向量 \left[ \matrix{H0[2,1]\cr H0[3,1]\cr H0[4,1]} \right]=\left[ \matrix{1 \cr 2 \cr 2} \right]
------------------------
k=1
X 向量 =\left[ \matrix{1 \cr 2 \cr 2} \right] 第1元 =1, sign(X_1)=1,長度 \Vert\; X \Vert\;=sqrt(\left[ \matrix{1 \cr 2 \cr 2} \right] \cdot \left[ \matrix{1 \cr 2 \cr 2} \right])=3
第1元是1,其餘為0的向量 e1=\left[ \matrix{1 \cr 0 \cr 0} \right]
X+sign(X_1)*\Vert\; X \Vert\;*e1=\left[ \matrix{1 \cr 2 \cr 2} \right]+1*3*\left[ \matrix{1 \cr 0 \cr 0} \right]=\left[ \matrix{4 \cr 2 \cr 2} \right]
長度 \Vert\; X+sign(X_1)*\Vert\; X \Vert\;*e1 \Vert\;=sqrt(\left[ \matrix{4 \cr 2 \cr 2} \right] \cdot \left[ \matrix{4 \cr 2 \cr 2} \right])=2\sqrt{6}
v=\left[ \matrix{4 \cr 2 \cr 2} \right]/ 2\sqrt{6}=\left[ \matrix{\displaystyle \frac{2}{\sqrt{6}}\cr \frac{1}{\sqrt{6}}\cr \frac{1}{\sqrt{6}}} \right]
\hat{P}=I-2*v.v^T=\left[ \matrix{1&0&0 \cr 0&1&0 \cr 0&0&1} \right]-2*\left[ \matrix{\displaystyle \frac{2}{\sqrt{6}} \cr \frac{1}{\sqrt{6}} \cr \frac{1}{\sqrt{6}}} \right] \cdot \left[ \matrix{\displaystyle \frac{2}{\sqrt{6}}& \frac{1}{\sqrt{6}}& \frac{1}{\sqrt{6}}} \right]=\left[ \matrix{\displaystyle -\frac{1}{3}&-\frac{2}{3}& -\frac{2}{3} \cr
-\frac{2}{3}&\frac{2}{3}&-\frac{1}{3} \cr
-\frac{2}{3}&-\frac{1}{3}&\frac{2}{3}} \right]
取向量 \left[ \matrix{1 \cr 2 \cr 2} \right] 找一個3*3 Householder 矩陣 \hat{P_1} ,使得 \left[ \matrix{\displaystyle -\frac{1}{3}&-\frac{2}{3}&-\frac{2}{3} \cr
-\frac{2}{3}&\frac{2}{3}&-\frac{1}{3} \cr
-\frac{2}{3}&-\frac{1}{3}&\frac{2}{3}} \right] \left[ \matrix{1 \cr 2 \cr 2} \right]=\left[ \matrix{-3 \cr 0 \cr 0} \right]
設 P_1=\left[ \matrix{\displaystyle 1&0&0&0&0&0 \cr
0&-\frac{1}{3}&-\frac{2}{3}&-\frac{2}{3}&0&0 \cr
0&-\frac{2}{3}&\frac{2}{3}&-\frac{1}{3}&0&0 \cr
0&-\frac{2}{3}&-\frac{1}{3}&\frac{2}{3}&0&0 \cr
0&0&0&0&1&0 \cr
0&0&0&0&0&1 } \right] ,計算 H_1=P_1^T.H_0.P_1=\left[ \matrix{\displaystyle
-\frac{2}{5}&-\frac{12}{25}&-\frac{7}{50}&\frac{3}{5}&\frac{23}{20}&-\frac{221}{40} \cr
-3&\frac{9}{10}&-\frac{53}{20}&-\frac{11}{50}&-\frac{129}{200}&\frac{87}{400} \cr
0&1&-\frac{13}{10}&\frac{26}{25}&\frac{289}{100}&-\frac{967}{200} \cr
0&2&0&-\frac{21}{5}&\frac{4}{5}&\frac{64}{5} \cr
0&-2&-1&2&-2&-2 \cr
0&0&0&0&1&-2 \cr} \right]
取向量 \left[ \matrix{H1[3,2]\cr H1[4,2]\cr H1[5,2]} \right]=\left[ \matrix{1 \cr 2 \cr -2} \right]
------------------------
k=2
X 向量 =\left[ \matrix{1 \cr 2 \cr -2} \right] 第1元 =1, sign(X_1)=1,長度 \Vert\; X \Vert\;=sqrt(\left[ \matrix{1 \cr 2 \cr -2} \right] \cdot \left[ \matrix{1 \cr 2 \cr -2} \right])=3
第1元是1,其餘為0的向量 e1=\left[ \matrix{1 \cr 0 \cr 0} \right]
X+sign(X_1)*\Vert\; X \Vert\;*e1=\left[ \matrix{1 \cr 2 \cr -2} \right]+1*3*\left[ \matrix{1 \cr 0 \cr 0} \right]=\left[ \matrix{4 \cr 2 \cr -2} \right]
長度 \Vert\; X+sign(X_1)*\Vert\; X \Vert\;*e1 \Vert\;=sqrt(\left[ \matrix{4 \cr 2 \cr -2} \right] \cdot \left[ \matrix{4 \cr 2 \cr -2} \right])=2\sqrt{6}
v=\left[ \matrix{4 \cr 2 \cr -2} \right]/ 2\sqrt{6}=\left[ \matrix{\displaystyle \frac{2}{\sqrt{6}}\cr \frac{1}{\sqrt{6}}\cr -\frac{1}{\sqrt{6}}} \right]
\hat{P}=I-2*v.v^T=\left[ \matrix{1&0&0 \cr 0&1&0 \cr 0&0&1} \right]-2*\left[ \matrix{\displaystyle \frac{2}{\sqrt{6}} \cr \frac{1}{\sqrt{6}} \cr -\frac{1}{\sqrt{6}}} \right] \cdot \left[ \matrix{\displaystyle \frac{2}{\sqrt{6}}& \frac{1}{\sqrt{6}}& -\frac{1}{\sqrt{6}}} \right]=\left[ \matrix{\displaystyle -\frac{1}{3}&-\frac{2}{3}& \frac{2}{3} \cr
-\frac{2}{3}&\frac{2}{3}&\frac{1}{3} \cr
\frac{2}{3}&\frac{1}{3}&\frac{2}{3}} \right]
取向量 \left[ \matrix{1 \cr 2 \cr -2} \right] 找一個3*3 Householder 矩陣 \hat{P_2} ,使得 \left[ \matrix{\displaystyle -\frac{1}{3}&-\frac{2}{3}&\frac{2}{3} \cr
-\frac{2}{3}&\frac{2}{3}&\frac{1}{3} \cr
\frac{2}{3}&\frac{1}{3}&\frac{2}{3}} \right] \left[ \matrix{1 \cr 2 \cr -2} \right]=\left[ \matrix{-3 \cr 0 \cr 0} \right]
設 P_2=\left[ \matrix{\displaystyle 1&0&0&0&0&0 \cr
0&1&0&0&0&0 \cr
0&0&-\frac{1}{3}&-\frac{2}{3}&\frac{2}{3}&0 \cr
0&0&-\frac{2}{3}&\frac{2}{3}&\frac{1}{3}&0 \cr
0&0&\frac{2}{3}&\frac{1}{3}&\frac{2}{3}&0 \cr
0&0&0&0&0&1 } \right] ,計算 H_2=P_2^T.H_1.P_2=\left[ \matrix{\displaystyle
-\frac{2}{5}&-\frac{12}{25}&\frac{31}{75}&\frac{263}{300}&\frac{131}{150}&-\frac{221}{40} \cr
-3&\frac{9}{10}&\frac{3}{5}&\frac{281}{200}&-\frac{227}{100}&\frac{87}{400} \cr
0&-3&-\frac{13}{3}&\frac{521}{300}&-\frac{39}{50}&-\frac{1651}{200} \cr
0&0&\frac{1}{3}&-\frac{439}{150}&-\frac{49}{25}&\frac{1109}{100} \cr
0&0&\frac{2}{3}&\frac{259}{150}&-\frac{6}{25}&-\frac{29}{100} \cr
0&0&\frac{2}{3}&\frac{1}{3}&\frac{2}{3}&-2 \cr} \right]
取向量 \left[ \matrix{H2[4,3]\cr H2[5,3]\cr H2[6,3]} \right]=\left[ \matrix{\displaystyle \frac{1}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right]
------------------------
k=3
X 向量 =\left[ \matrix{\displaystyle \frac{1}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right] 第1元 =\displaystyle \frac{1}{3}, sign(X_1)=1,長度 \Vert\; X \Vert\;=sqrt(\left[ \matrix{\displaystyle \frac{1}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right] \cdot \left[ \matrix{\displaystyle \frac{1}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right])=1
第1元是1,其餘為0的向量 e1=\left[ \matrix{1 \cr 0 \cr 0} \right]
X+sign(X_1)*\Vert\; X \Vert\;*e1=\left[ \matrix{\displaystyle \frac{1}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right]+1*1*\left[ \matrix{1 \cr 0 \cr 0} \right]=\left[ \matrix{\displaystyle \frac{4}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right]
長度 \Vert\; X+sign(X_1)*\Vert\; X \Vert\;*e1 \Vert\;=sqrt(\left[ \matrix{\displaystyle \frac{4}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right] \cdot \left[ \matrix{\displaystyle \frac{4}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right])=\frac{2^{3/2}}{\sqrt{3}}
v=\left[ \matrix{\displaystyle \frac{4}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right]/ \frac{2^{3/2}}{\sqrt{3}}=\left[ \matrix{\displaystyle \frac{\sqrt{2}}{\sqrt{3}}\cr \frac{1}{\sqrt{2}\sqrt{3}}\cr \frac{1}{\sqrt{2}\sqrt{3}}} \right]
\hat{P}=I-2*v.v^T=\left[ \matrix{1&0&0 \cr 0&1&0 \cr 0&0&1} \right]-2*\left[ \matrix{\displaystyle \frac{\sqrt{2}}{\sqrt{3}} \cr \frac{1}{\sqrt{2}\sqrt{3}} \cr \frac{1}{\sqrt{2}\sqrt{3}}} \right] \cdot \left[ \matrix{\displaystyle \frac{\sqrt{2}}{\sqrt{3}}& \frac{1}{\sqrt{2}\sqrt{3}}& \frac{1}{\sqrt{2}\sqrt{3}}} \right]=\left[ \matrix{\displaystyle -\frac{1}{3}&-\frac{2}{3}& -\frac{2}{3} \cr
-\frac{2}{3}&\frac{2}{3}&-\frac{1}{3} \cr
-\frac{2}{3}&-\frac{1}{3}&\frac{2}{3}} \right]
取向量 \left[ \matrix{\displaystyle \frac{1}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right] 找一個3*3 Householder 矩陣 \hat{P_3} ,使得 \left[ \matrix{\displaystyle -\frac{1}{3}&-\frac{2}{3}&-\frac{2}{3} \cr
-\frac{2}{3}&\frac{2}{3}&-\frac{1}{3} \cr
-\frac{2}{3}&-\frac{1}{3}&\frac{2}{3}} \right] \left[ \matrix{\displaystyle \frac{1}{3} \cr \frac{2}{3} \cr \frac{2}{3}} \right]=\left[ \matrix{-1 \cr 0 \cr 0} \right]
設 P_3=\left[ \matrix{\displaystyle 1&0&0&0&0&0 \cr
0&1&0&0&0&0 \cr
0&0&1&0&0&0 \cr
0&0&0&-\frac{1}{3}&-\frac{2}{3}&-\frac{2}{3} \cr
0&0&0&-\frac{2}{3}&\frac{2}{3}&-\frac{1}{3} \cr
0&0&0&-\frac{2}{3}&-\frac{1}{3}&\frac{2}{3} } \right] ,計算 H_3=P_3^T.H_2.P_3=\left[ \matrix{\displaystyle
-\frac{2}{5}&-\frac{12}{25}&\frac{31}{75}&\frac{632}{225}&\frac{3311}{1800}&-\frac{4103}{900} \cr
-3&\frac{9}{10}&\frac{3}{5}&\frac{9}{10}&-\frac{1009}{400}&-\frac{7}{200} \cr
0&-3&-\frac{13}{3}&\frac{49}{9}&\frac{1933}{1800}&-\frac{5761}{900} \cr
0&0&-1&\frac{4}{3}&\frac{1111}{900}&-\frac{587}{450} \cr
0&0&0&3&\frac{209}{225}&-\frac{1556}{225} \cr
0&0&0&4&\frac{2729}{900}&-\frac{3343}{450} \cr} \right]
取向量 \left[ \matrix{H3[5,4]\cr H3[6,4]} \right]=\left[ \matrix{\displaystyle 3 \cr 4} \right]
------------------------
X 向量 =\left[ \matrix{3 \cr 4} \right] 第1元 =3, sign(X_1)=1,長度 \Vert\; X \Vert\;=sqrt(\left[ \matrix{3 \cr 4} \right] \cdot \left[ \matrix{3 \cr 4} \right])=5
第1元是1,其餘為0的向量 e1=\left[ \matrix{1 \cr 0} \right]
X+sign(X_1)*\Vert\; X \Vert\;*e1=\left[ \matrix{3 \cr 4} \right]+1*5*\left[ \matrix{1 \cr 0} \right]=\left[ \matrix{8 \cr 4 } \right]
長度 \Vert\; X+sign(X_1)*\Vert\; X \Vert\;*e1 \Vert\;=sqrt(\left[ \matrix{8 \cr 4} \right] \cdot \left[ \matrix{8 \cr 4} \right])=4\sqrt{5}
v=\left[ \matrix{8 \cr 4} \right]/ 4\sqrt{5}=\left[ \matrix{\displaystyle \frac{2}{\sqrt{5}}\cr \frac{1}{\sqrt{5}}} \right]
\hat{P}=I-2*v.v^T=\left[ \matrix{1&0 \cr 0&1} \right]-2*\left[ \matrix{\displaystyle \frac{2}{\sqrt{5}} \cr \frac{1}{\sqrt{5}}} \right] \cdot \left[ \matrix{\displaystyle \frac{2}{\sqrt{5}}& \frac{1}{\sqrt{5}}} \right]=\left[ \matrix{\displaystyle -\frac{3}{5}&-\frac{4}{5} \cr
-\frac{4}{5}&\frac{3}{5}} \right]
取向量 \left[ \matrix{3 \cr 4} \right] 找一個2*2 Householder 矩陣 \hat{P_4} ,使得 \left[ \matrix{\displaystyle -\frac{3}{5}&-\frac{4}{5} \cr
-\frac{4}{5}&\frac{3}{5}} \right] \left[ \matrix{3 \cr 4} \right]=\left[ \matrix{-5 \cr 0 \cr} \right]
設 P_4=\left[ \matrix{\displaystyle 1&0&0&0&0&0 \cr
0&1&0&0&0&0 \cr
0&0&1&0&0&0 \cr
0&0&0&1&0&0 \cr
0&0&0&0&-\frac{3}{5}&-\frac{4}{5} \cr
0&0&0&0&-\frac{4}{5}&\frac{3}{5} } \right] ,計算 H_4=P_4^T.H_3.P_4=\left[ \matrix{\displaystyle
-\frac{2}{5}&-\frac{12}{25}&\frac{31}{75}&\frac{632}{225}&\frac{22891}{9000}&-\frac{18931}{4500} \cr
-3&\frac{9}{10}&\frac{3}{5}&\frac{9}{10}&\frac{3083}{2000}&\frac{1997}{1000} \cr
0&-3&-\frac{13}{3}&\frac{49}{9}&\frac{40289}{9000}&-\frac{21149}{4500} \cr
0&0&-1&\frac{4}{3}&\frac{1363}{4500}&-\frac{3983}{2250} \cr
0&0&0&-5&-\frac{35348}{5625}&\frac{47486}{5625} \cr
0&0&0&0&-\frac{33881}{22500}&-\frac{2429}{11250} \cr} \right]
(%o5) \left[ \matrix{\displaystyle
-\frac{2}{5}&-\frac{12}{25}&\frac{31}{75}&\frac{632}{225}&\frac{22891}{9000}&-\frac{18931}{4500} \cr
-3&\frac{9}{10}&\frac{3}{5}&\frac{9}{10}&\frac{3083}{2000}&\frac{1997}{1000} \cr
0&-3&-\frac{13}{3}&\frac{49}{9}&\frac{40289}{9000}&-\frac{21149}{4500} \cr
0&0&-1&\frac{4}{3}&\frac{1363}{4500}&-\frac{3983}{2250} \cr
0&0&0&-5&-\frac{35348}{5625}&\frac{47486}{5625} \cr
0&0&0&0&-\frac{33881}{22500}&-\frac{2429}{11250} \cr} \right]
|
