Ex.19 の続き
Ex.19 をめげずに続けてみましょう。
次は $(2,3)$-成分をターゲットにして、
再び $(1,2)$-成分、$(1,3)$-成分、$\cdots$ と
３周回してみると、
めでたく対角化に成功します：
\begin{align}
A &=\mat{rrr}{ 1.0000000000 & 2.0000000000 & 3.0000000000 \\
2.0000000000 & 5.0000000000 & 4.0000000000 \\
3.0000000000 & 4.0000000000 & 7.0000000000 \\
}
\end{align}
------------------------------ round 1 ------------------------------
\begin{align}
R_1 &= R(1, 2; -0.3926990817)\\
A_1 &= {}^tR_1\, A\, R_1\\
&=\mat{rrr}{ 0.1715728753 & -0.0000000000 & 1.2409048681 \\
-0.0000000000 & 5.8284271247 & 4.8435684271 \\
1.2409048681 & 4.8435684271 & 7.0000000000 \\
}
\end{align}
\begin{align}
R_2 &= R(1, 3; -0.1743043552)\\
A_2 &= {}^tR_2\, A_1\, R_2 \\
&=\mat{rrr}{ -0.0469396931 & -0.8399865337 & 0.0000000000 \\
-0.8399865337 & 5.8284271247 & 4.7701758596 \\
0.0000000000 & 4.7701758596 & 7.2185125684 \\
}
\end{align}
\begin{align}
R_3 &= R(2, 3; -0.7130543098)\\
A_3 &= {}^tR_3\, A_2\, R_3 \\
&=\mat{rrr}{ -0.0469396931 & -0.6353384637 & -0.5494746703 \\
-0.6353384637 & 1.7029240147 & 0.0000000000 \\
-0.5494746703 & -0.0000000000 & 11.3440156784 \\
}
\end{align}
----------------------------- round 2 ------------------------------
\begin{align}
R_4 &= R(1, 2; 0.3140332911)\\
A_4 &= {}^tR_4\, A_3\, R_4 \\
&=\mat{rrr}{ -0.2532851914 & 0.0000000000 & -0.5226028516 \\
0.0000000000 & 1.9092695130 & 0.1697311779 \\
-0.5226028516 & 0.1697311779 & 11.3440156784 \\
}
\end{align}
\begin{align}
R_5 &= R(1, 3; 0.0449410398)\\
A_5 &= {}^tR_5\, A_4\, R_5 \\
&=\mat{rrr}{ -0.2767873315 & 0.0076253282 & 0.0000000000 \\
0.0076253282 & 1.9092695130 & 0.1695598040 \\
0.0000000000 & 0.1695598040 & 11.3675178185 \\
}
\end{align}
\begin{align}
R_6 &= R(2, 3; -0.0179195128)\\
A_6 &= {}^tR_6\, A_5\, R_6 \\
&=\mat{rrr}{ -0.2767873315 & 0.0076241040 & 0.0001366349 \\
0.0076241040 & 1.9062307587 & -0.0000000000 \\
0.0001366349 & -0.0000000000 & 11.3705565729 \\
}
\end{align}
------------------------------ round 3 ------------------------------
\begin{align}
R_7 &= R(1, 2; -0.0034924035)\\
A_7 &= {}^tR_7\, A_6\, R_7 \\
&=\mat{rrr}{ -0.2768139581 & 0.0000000000 & 0.0001366340 \\
0.0000000000 & 1.9062573852 & 0.0000004772 \\
0.0001366340 & 0.0000004772 & 11.3705565729 \\
}
\end{align}
\begin{align}
R_8 &= R(1, 3; -0.0000117309)\\
A_8 &= {}^tR_8\, A_7\, R_8 \\
&=\mat{rrr}{ -0.2768139597 & -0.0000000000 & -0.0000000000 \\
-0.0000000000 & 1.9062573852 & 0.0000004772 \\
0.0000000000 & 0.0000004772 & 11.3705565745 \\
}
\end{align}
\begin{align}
\require{color}
R_9 &= R(2, 3; -0.0000000504)\\
A_9 &= {}^tR_9\, A_8\, R_9 \\
&=\mat{rrr}{ -0.2768139597 & \textcolor{red}{-0.0000000000} & \textcolor{red}{-0.0000000000} \\
\textcolor{red}{-0.0000000000} & 1.9062573852 & \textcolor{red}{0.0000000000} \\
\textcolor{red}{0.0000000000} & \textcolor{red}{-0.0000000000} & 11.3705565745 \\
}
\end{align}
すなわち
\begin{align}
P &= R_1\,R_2\,R_3\,R_4\,R_5\,R_6\,R_7\,R_8\,R_9 \\
&=\mat{rrr}{ 0.9385567220 & -0.1080624304 & 0.3277709425 \\
-0.1070043104 & 0.8118025918 & 0.5740441007 \\
-0.3281179013 & -0.5738458531 & 0.7503596336 \\
}\\
\end{align}
を用いて
\begin{align}
P^{-1}\,A\,P
&={}^tP\,A\,P \\
&=\mat{rrr}{ -0.2768139597 & -0.0000000000 & 0.0000000000 \\
-0.0000000000 & 1.9062573852 & -0.0000000000 \\
0.0000000000 & -0.0000000000 & 11.3705565745 \\
}
\end{align}
と対角化できました。
シリアル・ヤコビ法
上述の方法はシリアル・ヤコビ法 ( Serial Jacobi Method ) と呼ばれます。
Alg.20 ( シリアル・ヤコビ法 )
入力： $n$ 次実対称行列 $A$
出力： $A$ の固有値・固有ベクトル $(\lambda_i, \vvv_i)$ ( $i =1,2,\cdots,n$ )
$P=E$
while (継続条件)
for $p = 1$ to $n$
for $q = p+1$ to $n$
$\theta=\left\{
\begin{array}{lll}
\dps{\frac{1}{2}\arctan\left(\frac{2a_{pq}}{a_{pp}-a_{qq}}\right)} & \mbox{ if} & a_{pp} \neq a_{qq} \\
\dps{\frac{\pi}{4}} & \mbox{ if} & a_{pp} = a_{qq}
\end{array}
\right.$
$R=R(p,q;\theta)$
(新 $A)={}^tR\,A\,R$
(新 $P)=P\,R$
$\lambda_i = a_{ii}$, $\vvv_i=(P$ の第 $i$ 列ベクトル) ( $i =1,2,\cdots,n$ ) を出力
継続条件は、非対角要素の絶対値の最大値が許容誤差より大きいとか、
反復回数で決めます。
古典的ヤコビ法
シリアル・ヤコビ法のように順番にターゲットを決めるのではなく、
絶対値最大の非対角要素をターゲットにする方法を古典的ヤコビ法と呼びます。
Alg.21 ( 古典的ヤコビ法 )
入力： $n$ 次実対称行列 $A$
出力： $A$ の固有値・固有ベクトル $(\lambda_i, \vvv_i)$ ( $i =1,2,\cdots,n$ )
$P=E$
while (継続条件)
絶対値最大の非対角要素を $a_{pq}$ ( $p \lt q$ ) とする
$\theta=\left\{
\begin{array}{lll}
\dps{\frac{1}{2}\arctan\left(\frac{2a_{pq}}{a_{pp}-a_{qq}}\right)} & \mbox{ if} & a_{pp} \neq a_{qq} \\
\dps{\frac{\pi}{4}} & \mbox{ if} & a_{pp} = a_{qq}
\end{array}
\right.$
$R=R(p,q;\theta)$
(新 $A)={}^tR\,A\,R$
(新 $P)=P\,R$
$\lambda_i = a_{ii}$, $\vvv_i=(P$ の第 $i$ 列ベクトル) ( $i =1,2,\cdots,n$ ) を出力
この方法はシリアル・ヤコビ法より効率よく非対角成分を小さくできますが、
行列サイズ $n$ が大きくなると絶対値最大の非対角要素を選択するのに時間が掛かるのが難点です。
閾値(しきいち)ヤコビ法
シリアル・ヤコビ法と古典的ヤコビ法のいいとこ取りで、
設定した閾値より大きい絶対値を持つ非対角要素だけをターゲットにする方法を閾値ヤコビ法 ( Threshold Jacobi Method ) と呼びます。
非対角要素の選択に掛かる時間を節約し、かつ、効率よく変形することができます。
Alg.22 ( 閾値ヤコビ法 )
入力： $n$ 次実対称行列 $A$
出力： $A$ の固有値・固有ベクトル $(\lambda_i, \vvv_i)$ ( $i =1,2,\cdots,n$ )
$P=E$
閾値 $\varepsilon$ を設定
while (継続条件)
for $p = 1$ to $n$
for $q = p+1$ to $n$
if $|\,a_{pq}\,| \gt \varepsilon$
$\theta=\left\{
\begin{array}{lll}
\dps{\frac{1}{2}\arctan\left(\frac{2a_{pq}}{a_{pp}-a_{qq}}\right)} & \mbox{ if} & a_{pp} \neq a_{qq} \\
\dps{\frac{\pi}{4}} & \mbox{ if} & a_{pp} = a_{qq}
\end{array}
\right.$
$R=R(p,q;\theta)$
(新 $A)={}^tR\,A\,R$
(新 $P)=P\,R$
$\varepsilon = \varepsilon / 10$
$\lambda_i = a_{ii}$, $\vvv_i=(P$ の第 $i$ 列ベクトル) ( $i =1,2,\cdots,n$ ) を出力
塩田は研究では閾値ヤコビ法を使っています。