Thought 主成分分析

目次

FrontPage
・Windows11
・Chromebook
・Random
・機械学習
・さくらVPS
・Fedora13
・SuSe10
・Docker
・Ansible
・Java
・Scala
・Python
・Ruby
・Lisp
・Computer
・GIS
・HTML
・Culture
・Link

訪問者

total: 7317
today: 4
yesterday: 3
now: 13

更新

最新の10件

人気の10件

MenuBar

主成分分析は何をするものか? †

• こんな入力データがあったとする

  \(X_{1} = (x_{11}, x_{21}, ... ,x_{I1}) \\ X_{2} = (x_{12}, x_{22}, ... ,x_{I2}) \\ \vdots \\ X_{N} = (x_{1N}, x_{2N}, ... ,x_{IN}) \)

• 主成分分析は、X の特徴をよく表すような Z (=主成分) を作ること

  \(z_{i} = (w_{1} \times x_{1i}) + (w_{2} \times x_{2i}) + ... + (w_{I} \times x_{Ii}) \)

• 具体的には、z の分散が最も大きくなるように (w₁, w₂, ... , w_I) を決めるのが主成分分析

考え方 †

• z の分散を計算する

  \(\sigma _{z} ^{2} = \frac{1}{N} \sum_{i=1}^{N} (z_{i} - \overline{z}) ^{2} \\ = \frac{1}{N} \sum_{i=1}^{N} ( w_{1} (x_{1i} - \overline{x_{1}}) + w_{2} (x_{2i} - \overline{x_{2}}) + ... + w_{I} (x_{Ii} - \overline{x_{I}}) ) ^{2} \\ = w_{1}^{2} \sigma _{1}^{2} + 2w_{1}w_{2} \sigma _{12} + 2w_{1}w_{3} \sigma _{13} + 2w_{1}w_{4} \sigma _{14} + \dots + 2w_{1}w_{I} \sigma _{1I} \\ + w_{2}^{2} \sigma _{2}^{2} + 2w_{2}w_{3} \sigma _{23} + 2w_{2}w_{4} \sigma _{24} + \dots + 2w_{2}w_{I} \sigma _{2I} \\ + w_{3}^{2} \sigma _{3}^{2} + 2w_{3}w_{4} \sigma _{34} + \dots + 2w_{3}w_{I} \sigma _{3I} \\ \vdots \\ + w_{I}^{2} \sigma _{I}^{2} \)

  (σ_ij は x_i と x_j の共分散)
• σ_z² を最大化したい ⇒ ラグランジュの未定乗数法 を使う

  \(maximize : f(w_1,w_2,...,w_I) = \sigma_z^2\\ subject~to : g(w_1,w_2,...,w_I) = w_1^2 + w_2^2 + ... + w_I^2 - 1 = 0\\ とすると\\ F(w_1,w_2,...,w_I,\lambda) = f() - \lambda g() \\ について \\ \\ \frac{\sigma F}{\sigma w_1} = 2 \sigma_{1}^{2} w_{1} + 2 \sigma_{12} w_{2} + 2 \sigma_{13} w_{3} + ... + 2 \sigma_{1I} w_{I} - \lambda 2 w_{1} = 0\\ \frac{\sigma F}{\sigma w_2} = 2 \sigma_{12} w_{1} + 2 \sigma_{2}^{2} w_{2} + 2 \sigma_{23} w_{3} + ... + 2 \sigma_{2I} w_{I} - \lambda 2 w_{2} = 0\\ \vdots \\ \frac{\sigma F}{\sigma w_I} = 2 \sigma_{1I} w_{1} + 2 \sigma_{2I} w_{2} + 2 \sigma_{3I} w_{3} + ... + 2 \sigma_{I}^2 w_{I} - \lambda 2 w_{I} = 0\\ を満たす (w_{1}, w_{2}, ... , w_{I}) が f() を最大化する \)

  整理すると

  \(\begin{bmatrix} \sigma_1^2 & \sigma_{12} & \dots & \sigma_{1I} \\ \sigma_{12} & \sigma_2^2 & \dots & \sigma_{2I} \\ \vdots & \vdots & \ddots & \vdots \\ \sigma_{1I} & \sigma_{2I} & \dots & \sigma_I^2 \\ \end{bmatrix} \begin{bmatrix} w_{1} \\ w_{2} \\ \vdots \\ w_{I} \end{bmatrix} = \lambda \begin{bmatrix} w_{1} \\ w_{2} \\ \vdots \\ w_{I} \end{bmatrix} \\ を満たす (w_{1}, w_{2}, ... , w_{I}, \lambda) を見つければいい \)

• これは、固有値問題
  • 解は I 個。固有値問題の解法の一つである QR 法 を使うと

    \(X = \begin{bmatrix} \overrightarrow{W_{1}} & \overrightarrow{W_{2}} & \dots & \overrightarrow{W_{I}} \end{bmatrix} \\ (\overrightarrow{W_{i}} = (w_{i1}, w_{i2}, \dots , w_{iI})^T) \\ \\ \Lambda = \begin{bmatrix} \lambda_{1} & \lambda_{2} & \dots & \lambda_{I} \end{bmatrix} \)

    形式で全ての解が求まる
  • \(\overrightarrow{W_{i}}\) 同士は一次独立
  • λ = σ_z² である
    

    \(\sigma _{z} ^{2} = (w_{1}^{2} \sigma _{1}^{2} + w_{1}w_{2} \sigma _{12} + \dots + w_{1}w_{I} \sigma _{1I}) + (w_{1}w_{2} \sigma _{12} + w_{2}^{2} \sigma _{1}^{2} + \dots + w_{2}w_{I} \sigma _{2I}) + \dots + (w_{1}w_{I} \sigma _{1I} + w_{2}w_{I} \sigma _{2I} + \dots + w_{I}^{2} \sigma _{I}^{2}) \\ = w_{1} (w_{1} \sigma _{1}^{2} + w_{2} \sigma _{12} + \dots + w_{I} \sigma _{1I}) + w_{2} (w_{1} \sigma _{12} + w_{2} \sigma _{1}^{2} + \dots + w_{I} \sigma _{2I}) + \dots + w_{I} (w_{1} \sigma _{1I} + w_{2} \sigma _{2I} + \dots + w_{I} \sigma _{I}^{2}) \\ = w_{1} (\lambda w_{1}) + w_{2} (\lambda w_{2}) + \dots + w_{I} (\lambda w_{I}) \\ = \lambda (w_{1}^2 + w_{2}^2 + \dots + w_{I}^2) \\ 制約条件より g() = w_{1}^2 + w_{2}^2 + \dots + w_{I}^2 - 1 = 0 \\ \therefore \sigma _{z} ^{2} = \lambda \)

  • i 番目の解の寄与率 (contribution)

    \(c_{i} = \frac{\lambda_{i}}{\lambda_{1} + \lambda_{2} + \dots + \lambda_{I}} \)

    i 番目の解で元データの分散の内どれくらいを実現できるかを示す指標
  • 累積寄与率 (cumulative contribution)

    \(s_{i} = \sum_{j=1}^i c_{j} \)

    解を寄与率順にソートしたとき 1 〜 i 番目までの解で元データの分散の内どれくらいを実現できるかを示す指標
    上位 20 % の解で、元の分散の 99.9 % を表すことができるなら、分析すべきデータ項目数を元データの項目数の 20 %にすることができる。

ラグランジュの未定乗数法 †

• 最大値が知りたい関数 f(x,y) があるとき
• x,y の適当な制約条件 g(x,y) = 0 を設定すると
• F(x,y,λ) = f(x,y) - λg(x,y) について

  \(\frac{\sigma F}{\sigma x} = \frac{\sigma F}{\sigma y} = \frac{\sigma F}{\sigma \lambda} = 0 \)

  となるような (x,y,λ) が 制約条件 g(x,y) = 0 のもとで f(x,y) を最大化する
• これは何をやってるの?
  F を x、y それぞれで偏微分して 0 になるということは、

  \(\frac{\sigma f}{\sigma x} = \lambda\frac{\sigma g}{\sigma x} \\ \frac{\sigma f}{\sigma y} = \lambda\frac{\sigma g}{\sigma y} \\ \)

  つまり、f の勾配と g の勾配の方向が同じ地点を探している。
  ⇒ f の等高線と g の接点を探している
  
• パラメータ数が増えても考え方は同じ
• 普通、制約条件には、各パラメータの自乗和 = 1 が使われる
  • ２次元の場合は円になる
  • 制約関数の勾配方向が進行方向でないのは変曲点だけ、円、球、... なら表面が全部変曲点

固有値問題 †

一般的な固有値問題 †

• 固有値問題
  正方行列 A に対して

  \(A \overrightarrow{x} = \lambda \overrightarrow{x} \\ \begin{bmatrix} a_{11} & a_{12} & \dots & a_{1N} \\ a_{21} & a_{22} & \dots & a_{2N} \\ \vdots & \vdots & \ddots & \vdots \\ a_{N1} & a_{N2} & \dots & a_{NN} \\ \end{bmatrix} \begin{bmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{N} \end{bmatrix} = \lambda \begin{bmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{N} \end{bmatrix} \\ \)

  となるような、固有値 λ、固有ベクトル \(\overrightarrow{x}\) を見つけること。
  
   
  幾何学的には、写像 A を行っても方向が変わらないベクトル \(\overrightarrow{x}\) を見つけること。
• 固有方程式

  \(A \overrightarrow{x} = \lambda \overrightarrow{x} \\ (\lambda I - A) \overrightarrow{x} = 0 \\ もしも (\lambda I - A) に逆行列が存在すると \\ (\lambda I - A)^{-1} (\lambda I - A) \overrightarrow{x} = (\lambda I - A)^{-1} 0 \\ \overrightarrow{x} = 0 となって固有ベクトルではなくなる。これはおかしい。 \\ \therefore (\lambda I - A) に逆行列は存在しない \Rightarrow (\lambda I - A) は正則でない \\ \Rightarrow det(\lambda I - A) = 0 (固有方程式) \)

  • 固有方程式はλのN次方程式になるので、これを解けば N 個の λ を求めることが出来て、λ_i に対する \(\overrightarrow{x}_i\) を求めればよい
  • でも、行列式 (\(det(\lambda I - A)\) (別の書き方では \(|\lambda I - A|\))) を展開できるのは、せいぜい 3 次まで。それ以上の高次になると、固有値・固有ベクトルの計算を現実的な時間でするのには、もうチョット工夫が必要。
• 固有値の性質
  • 転置行列の固有値は同じ (固有ベクトルは変わる)

    \(A x = \lambda x \\ A^{T} x' = \lambda x' \)

  • A^k の固有値は λ^k、固有ベクトルは x

    \(A x = \lambda x \\ A^{k} x = \lambda^{k} x \)

  • A^-1 の固有値は 1/λ、固有ベクトルは x

    \(A x = \lambda x \\ A^{-1} x = \frac{1}{\lambda} x \)

• 相似変換
  • A の相似変換

    \(X^{-1}AX = A' (Xは正則行列) \)

  • 相似変換をしても固有値λは変わらない

    \(eigenvalue(A) = det(\lambda I - A) = det(X^{-1})det(\lambda I - A)det(X) = det(X^{-1}(\lambda I - A)X) = det(\lambda I - X^{-1}AX) = det(\lambda I - A') = eigenvalue(A') \)

  • 固有ベクトル \(\overrightarrow{x}\) は \(X^{-1} \overrightarrow{x}\) になる。
    A の固有ベクトル \(\overrightarrow{x}\)、A'の固有ベクトル \(\overrightarrow{x'}\) について

    \(A' \overrightarrow{x'} = \lambda \overrightarrow{x'} \\ X^{-1}AX \overrightarrow{x'} = \lambda \overrightarrow{x'} \\ XX^{-1}AX \overrightarrow{x'} = \lambda X \overrightarrow{x'} \\ A (X \overrightarrow{x'}) = \lambda (X \overrightarrow{x'}) \\ \therefore \overrightarrow{x} = X \overrightarrow{x'} \\ \overrightarrow{x'} = X^{-1} \overrightarrow{x} \)

• 行列の対角化
  • 正方行列 A の固有値 λ_i、固有ベクトル \(\overrightarrow{x_{i}}\) について

    \(P = [\overrightarrow{x_{1}}~\overrightarrow{x_{2}} \dots \overrightarrow{x_{n}}] \\ \Lambda = \begin{bmatrix} \lambda_1 & 0 & \dots & 0 \\ 0 & \lambda_2 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \lambda_n \end{bmatrix} \)

    とすると

    \(AP = [A\overrightarrow{x_{1}}~A\overrightarrow{x_{2}} \dots A\overrightarrow{x_{n}}] \\ = [\lambda_{1}\overrightarrow{x_{1}}~\lambda_{2}\overrightarrow{x_{2}} \dots \lambda_{n}\overrightarrow{x_{n}}] \\ = [\overrightarrow{x_{1}}~\overrightarrow{x_{2}} \dots \overrightarrow{x_{n}}] \begin{bmatrix} \lambda_1 & 0 & \dots & 0 \\ 0 & \lambda_2 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \lambda_n \end{bmatrix} \\ = P\Lambda \\ \therefore P^{-1}AP = \Lambda \)

  • Aの対角化 : \(P^{-1}AP = \Lambda\)
  • Aの標準系 : \(\Lambda\)
  • 対角化は、変換結果が標準型になる特殊な相似変換
    ⇒ 対角化された行列 Λ の固有値は自明 (対角成分) && 対角化された行列 Λ の固有値は元の行列 A の固有値と同じ
    ⇒ 固有値問題を解くことは、対角化することと同じ
    ⇒ 固有値問題を解くアルゴリズムは、何らかの相似変換を繰り返して対角化を行う方向性で考えられている。「何らかの相似変換」が腕の見せ所

    \(X_{1}^{-1}AX_{1} \\ X_{2}^{-1}(X_{1}^{-1}AX_{1})X_{2} \\ \vdots \\ X_{\infty}^{-1}( \dots (X_{2}^{-1}(X_{1}^{-1}AX_{1})X_{2}) \dots )X_{\infty} \approx P^{-1}AP = \Lambda \)

対称行列の固有値問題 †

• 対称行列

  \(a_{ij} = a_{ji}\\ A^T = A \)

• この性質により一般的な固有値問題よりも簡単に計算できる
  ⇒ 主成分分析は相関行列の固有値問題なので、対称行列の固有値問題を考えればいい (σ₁₂ = σ₂₁)
• A が対称行列のとき、固有ベクトルはそれぞれ直交する(一次独立になる)

  \(固有値 \lambda_{i}, \lambda_{j}\\ 固有ベクトル x_{i}, x_{j} \\ について \\ \lambda_{j} (x_{i} \bullet x_{j}) = x_{i}^T (\lambda_{j} x_{j}) = x_{i}^T (A x_{j}) = x_{i}^T A^T x_{j} = (A x_{i})^T x_{j} = \lambda_{i} (x_{i} \bullet x_{j}) \\ \therefore \lambda_{i} (x_{i} \bullet x_{j}) - \lambda_{j} (x_{i} \bullet x_{j}) = (\lambda_{i} - \lambda_{j})(x_{i} \bullet x_{j}) = 0\\ ここで、\lambda_{i} \neq \lambda_{j} だから x_{i} \bullet x_{j} = 0 \)

• A が対称行列のとき、正規直交行列で対角化できる
  正方行列 A の固有値 λ_i、固有ベクトル \(\overrightarrow{x_{i}}\) について

  \(U = [\frac{\overrightarrow{x_{1}}}{|\overrightarrow{x_{1}}|}~\frac{\overrightarrow{x_{2}}}{|\overrightarrow{x_{2}}|} \dots \frac{\overrightarrow{x_{n}}}{|\overrightarrow{x_{n}}|}] \\ \Lambda = \begin{bmatrix} \lambda_1 & 0 & \dots & 0 \\ 0 & \lambda_1 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \lambda_n \end{bmatrix} \)

  とすると、さっき計算したように固有ベクトル同士は直行するから、U は正規直交行列 (構成するベクトル同士が直交して、各々のノルムは1)

  \(AU = [A \frac{\overrightarrow{x_{1}}}{|\overrightarrow{x_{1}}|}~A \frac{\overrightarrow{x_{2}}}{|\overrightarrow{x_{2}}|} \dots A \frac{\overrightarrow{x_{n}}}{|\overrightarrow{x_{n}}|}] \\ = [\lambda_{1} \frac{\overrightarrow{x_{1}}}{|\overrightarrow{x_{1}}|}~\lambda_{2} \frac{\overrightarrow{x_{2}}}{|\overrightarrow{x_{2}}|} \dots \lambda_{n} \frac{\overrightarrow{x_{n}}}{|\overrightarrow{x_{n}}|}] \\ = [\frac{\overrightarrow{x_{1}}}{|\overrightarrow{x_{1}}|}~\frac{\overrightarrow{x_{2}}}{|\overrightarrow{x_{2}}|} \dots \frac{\overrightarrow{x_{n}}}{|\overrightarrow{x_{n}}|}] \begin{bmatrix} \lambda_1 & 0 & \dots & 0 \\ 0 & \lambda_1 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \lambda_n \end{bmatrix} \\ = U\Lambda \\ \therefore U^{-1}AU = \Lambda \)

  正規直交行列は、構成するベクトル同士が直交して、各々のベクトルのノルムが１なので \(U^{T}U=I\) つまり \(U^{-1}=U^T\)

  \(\therefore U^{-1}AU = U^{T}AU = \Lambda \)

固有値から固有ベクトルを求める †

• 相似変換で、行列を簡略化して固有値を求めることができることは分かる
• で、どのようにして固有ベクトルを求めるか?
• 正直にやると

  \(A \overrightarrow{x} = \lambda_{k} \overrightarrow{x} \)

  のN次不定方程式を k=1,2,...,N まで N個解けばいい ⇒ (°Д°)ﾒﾝﾄﾞｸｾ(°Д°)ﾒﾝﾄﾞｸｾ(°Д°)ﾒﾝﾄﾞｸｾ
• 固有方程式を使う

  \((A-\lambda_{1}I)(A-\lambda_{2}I)\dots(A-\lambda_{N}I) = 0 \\ (A-\lambda_{1}I)M_1 = (A-\lambda_{1}I) \begin{bmatrix} \overrightarrow{v_1} & \overrightarrow{v_2} & \dots & \overrightarrow{v_N} \end{bmatrix} = 0 \\ \therefore (A-\lambda_{1}I) \overrightarrow{v_1} = 0 \\ A \overrightarrow{v_1} = \lambda_{1}\overrightarrow{v_1} \)

  ⇒ 固有値 λ₁ に対応する固有ベクトル(の一つ)は、v₁ になる
• つまり固有値 λ_k に対応する固有ベクトルは

  \(\prod_{i=1}^{N (i \neq k)} (A-\lambda_{i}I) \)

  の任意の列ベクトル
• N次元空間での固有ベクトルの方向だけわかればいいので、一つわかればいい
  (実数解があるなら M の列ベクトルの方向は全部同じはず)

QR法 †

• QR分解

  \(A = QR \)

  A を Q(直交行列) と R(上三角行列) の積の形として表すこと
   
  • 幾何学的には、線形変換を回転(Rotate)、拡大・縮小(Extension/Reduction)、歪み(Skew) に分解している

    \(\begin{bmatrix}a & b \\c & d \end{bmatrix} \begin{bmatrix}x \\y \end{bmatrix} = \begin{bmatrix}cos\theta & -sin\theta \\ sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix}p & 0 \\0 & q \end{bmatrix} \begin{bmatrix}1 & r/p \\0 & 1 \end{bmatrix} \begin{bmatrix}x \\y \end{bmatrix} \)

    
  • QR分解のやり方
    • グラム・シュミットの直交化
    • ハウスホールダー変換
• QR法

  \(A_{1} = Q_{1}R_{1} \\ R_{1}Q_{1} = A_{2} \\ \vdots \\ A_{k} = Q_{k}R_{k} \\ R_{k}Q_{k} = A_{k+1} \\ \vdots \\ A_{\infty} = \begin{bmatrix} a_{11} & a_{12} & a_{13} & \dots & a_{1N} \\ 0 & a_{22} & a_{23} & \dots & a_{2N} \\ 0 & 0 & a_{33} & \dots & a_{3N} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & a_{NN} \\ \end{bmatrix} \)

  A_∞ は、上三角行列になり、その対角成分が固有値になる
• 一回分の変換 \(A_{k} \Rightarrow A_{k+1}\) は相似変換 (固有値 λ は変わらない)

  \(A_{k+1} = R_{k}Q_{k} \\ Q_{k}A_{k+1} = Q_{k}R_{k}Q_{k} \\ Q_{k}A_{k+1} = A_{k}Q_{k} \\ A_{k+1} = Q_{k}^{-1}A_{k}Q_{k} \\ \)

  ⇒ QR法は、相似変換を繰り返して (固有値λが同じ) 上三角行列を作っている
   
• なぜ、上三角行列の対角成分が固有値になるのか?
  固有値 λ は、固有方程式

  \(det(\lambda I - A_{\infty}) = 0 \)

  の解。

  \(det(\lambda I - A_{\infty}) \\ = det( \begin{bmatrix} \lambda & 0 & 0 & \dots & 0 \\ 0 & \lambda & 0 & \dots & 0 \\ 0 & 0 & \lambda & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & \lambda \\ \end{bmatrix} - \begin{bmatrix} a_{11} & a_{12} & a_{13} & \dots & a_{1N} \\ 0 & a_{22} & a_{23} & \dots & a_{2N} \\ 0 & 0 & a_{33} & \dots & a_{3N} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & a_{NN} \\ \end{bmatrix} ) \\ = det( \begin{bmatrix} \lambda - a_{11} & -a_{12} & -a_{13} & \dots & -a_{1N} \\ 0 & \lambda - a_{22} & -a_{23} & \dots & -a_{2N} \\ 0 & 0 & \lambda - a_{33} & \dots & -a_{3N} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & \lambda - a_{NN} \\ \end{bmatrix} ) \\ = (\lambda - a_{11})(\lambda - a_{22})(\lambda - a_{33}) \dots (\lambda - a_{NN}) = 0 \\ ※ 三角行列の行列式の値は、対角成分の積 \\ \therefore \lambda = a_{11}, a_{22}, a_{33} \dots a_{NN} \)

  ⇒ とにかく、行列 A を三角行列に相似変換できれば固有値λが計算できる
   
• なぜ、QR法は上三角行列に収束するのか?
  → QR法が上三角行列に収束する証明

グラム・シュミットの直交化 †

• QR分解 の一手法
• Gram-Schmidt orthonormalization
• QR分解

  \(A = QR \)

  A を Q(直交行列) と R(上三角行列) の積の形として表すこと
   
• グラム・シュミットの直交化の考え方

  \(A = QR \)

  をベクトルの集まり A を互いに直交するベクトルの集まり Q に変換する。そして、その変換行列が R と考える。
  つまり、

  \(A = \begin{bmatrix} \overrightarrow{a_{1}} & \overrightarrow{a_{2}} & \dots & \overrightarrow{a_{N}} \end{bmatrix} \\ Q = \begin{bmatrix} \overrightarrow{q_{1}} & \overrightarrow{q_{2}} & \dots & \overrightarrow{q_{N}} \end{bmatrix} \\ (ただし、\overrightarrow{q_{1}},\overrightarrow{q_{2}},\dots,\overrightarrow{q_{N}} は互いに直交) \\ R は、A = QR とするための適当な変換行列 \)

  いま、

  \(\overrightarrow{q_{1}},\dots,\overrightarrow{q_{i}} \)

  まで、互いに直行するベクトルができているとすると、\(\overrightarrow{q_{i+1}}\) は、

  \(\overrightarrow{q_{1}},\dots,\overrightarrow{q_{i}} \)

  からなる超平面に \(\overrightarrow{a_{i+1}}\) から垂線を下ろしたものになる。
  このようにして、i=1 から N まで順に、直行するベクトル qi を増やしていく
  • i=1
    

    \(\overrightarrow{t} + \overrightarrow{q_{2}} = \overrightarrow{a_{2}}\\ \therefore \overrightarrow{q_{2}} = \overrightarrow{a_{2}} - \overrightarrow{t} \)

    ここで、

    \(\overrightarrow{t} = \frac{|\overrightarrow{a_{2}}|cos\theta} {|\overrightarrow{q_{1}}|} \overrightarrow{q_{1}} = \frac{|\overrightarrow{q_{1}}||\overrightarrow{a_{2}}|cos\theta} {|\overrightarrow{q_{1}}||\overrightarrow{q_{1}}|cos0} \overrightarrow{q_{1}} = \frac{\overrightarrow{q_{1}} \bullet \overrightarrow{a_{2}}} {\overrightarrow{q_{1}} \bullet \overrightarrow{q_{1}}} \overrightarrow{q_{1}} \\ \)

    \(\therefore \overrightarrow{q_{2}} = \overrightarrow{a_{2}} - \frac {\overrightarrow{q_{1}} \bullet \overrightarrow{a_{2}}} {\overrightarrow{q_{1}} \bullet \overrightarrow{q_{1}}} \overrightarrow{q_{1}} \)

  • i=2
    
    i=1の時と同様に考えると

    \(\therefore \overrightarrow{q_{3}} = \overrightarrow{a_{3}} - \frac {\overrightarrow{q_{1}} \bullet \overrightarrow{a_{3}}} {\overrightarrow{q_{1}} \bullet \overrightarrow{q_{1}}} \overrightarrow{q_{1}} - \frac {\overrightarrow{q_{2}} \bullet \overrightarrow{a_{3}}} {\overrightarrow{q_{2}} \bullet \overrightarrow{q_{2}}} \overrightarrow{q_{2}} \)

  • i=k

    \(\therefore \overrightarrow{q_{k+1}} = \overrightarrow{a_{k+1}} - \frac {\overrightarrow{q_{1}} \bullet \overrightarrow{a_{k+1}}} {\overrightarrow{q_{1}} \bullet \overrightarrow{q_{1}}} \overrightarrow{q_{1}} - \frac {\overrightarrow{q_{2}} \bullet \overrightarrow{a_{k+1}}} {\overrightarrow{q_{2}} \bullet \overrightarrow{q_{2}}} \overrightarrow{q_{2}} - \dots - \frac {\overrightarrow{q_{k}} \bullet \overrightarrow{a_{k+1}}} {\overrightarrow{q_{k}} \bullet \overrightarrow{q_{k}}} \overrightarrow{q_{k}} \)

• 行列形式に書き換えると

  \(\begin{bmatrix} \overrightarrow{a_{1}} & \overrightarrow{a_{2}} & \overrightarrow{a_{3}} & \dots & \overrightarrow{a_{N}} \end{bmatrix} = \begin{bmatrix} \overrightarrow{q_{1}} & \overrightarrow{q_{2}} & \overrightarrow{q_{3}} & \dots & \overrightarrow{q_{N}} \end{bmatrix} \begin{bmatrix} 1 & \frac {\overrightarrow{q_{1}} \bullet \overrightarrow{a_{2}}} {\overrightarrow{q_{1}} \bullet \overrightarrow{q_{1}}} & \frac {\overrightarrow{q_{1}} \bullet \overrightarrow{a_{3}}} {\overrightarrow{q_{1}} \bullet \overrightarrow{q_{1}}} & \dots & \frac {\overrightarrow{q_{1}} \bullet \overrightarrow{a_{N}}} {\overrightarrow{q_{1}} \bullet \overrightarrow{q_{1}}} \\ 0 & 1 & \frac {\overrightarrow{q_{2}} \bullet \overrightarrow{a_{3}}} {\overrightarrow{q_{2}} \bullet \overrightarrow{q_{2}}} & \dots & \frac {\overrightarrow{q_{2}} \bullet \overrightarrow{a_{N}}} {\overrightarrow{q_{2}} \bullet \overrightarrow{q_{2}}} \\ 0 & 0 & 1 & \dots & \frac {\overrightarrow{q_{3}} \bullet \overrightarrow{a_{N}}} {\overrightarrow{q_{3}} \bullet \overrightarrow{q_{3}}} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots \ & 1 \end{bmatrix} \)

  A=QR (Q は直交行列、R は上三角行列)
   
• QR分解の方法として理解しやすいけど、逐次処理なので並列化できないし、計算誤差もたまりやすい → 計算機を使う場合には次のハウスホールダー変換を使ったほうがいい

ハウスホールダー変換 †

• Householder transformation
• ハウスホールダー変換の考え方
  変換元の行列 A をベクトルの集まりと考える

  \(A = \begin{bmatrix} \overrightarrow{a_{1}} & \overrightarrow{a_{2}} & \dots & \overrightarrow{a_{N}} \end{bmatrix} \)

  これらのベクトルを原点を含むある超平面で反転させて、\(\overrightarrow{a_{i}}\) の 第 j 要素を 0 にする変換を考える
  ⇒ そういう変換は、各ベクトルのノルムや相対的な関係は変わらないはずなので、固有値λが変わらない相似変換になるはず
  ⇒ 0 を原点を含む超平面で反転させても 0 なので、いったん 0 になった要素は 0 のまま。
  そこで、このような変換を次々使えば上三角行列を作ることができるはず
  ⇒ QR分解に使えんじゃね?
   
  いま a を a' に変更して、第2要素 (y要素) を 0 にすることを考える
  

  \(\overrightarrow{v} = \overrightarrow{a} - \overrightarrow{a'} \\ \overrightarrow{u} = \frac{\overrightarrow{v}}{|\overrightarrow{v}|} \)

  とすると、

  \(\overrightarrow{a'} = \overrightarrow{a} - \overrightarrow{v} = \overrightarrow{a} - 2 (|\overrightarrow{a}| cos \theta)\overrightarrow{u} = \overrightarrow{a} - 2 (\overrightarrow{a} \bullet \overrightarrow{u})\overrightarrow{u} \\ = \overrightarrow{a} - 2 \begin{bmatrix} u_{1}^2 & u_{1}u_{2} & \dots & u_{1}u_{N} \\ u_{1}u_{2} & u_{2}^2 & \dots & u_{2}u_{N} \\ \vdots & \vdots & \ddots & \vdots \\ u_{1}u_{N} & u_{2}u_{N} & \dots & u_{N}^2 \end{bmatrix} \overrightarrow{a} \\ = \overrightarrow{a} - 2 uu^T \overrightarrow{a} = (E - 2 uu^T) \overrightarrow{a} = H(u) \overrightarrow{a} \)

• ハウスホールダー変換
  • 適当なベクトル v について

    \(\overrightarrow{u} = \frac{\overrightarrow{v}}{|\overrightarrow{v}|} \)

    としたとき

    \(A' = H(u)A \)

    をハウスホルダー変換と呼ぶ。変換行列

    \(H(u) = \begin{bmatrix} 1-2u_{1}^2 & -2u_{1}u_{2} & \dots & -2u_{1}u_{N} \\ -2u_{1}u_{2} & 1-2u_{2}^2 & \dots & -2u_{2}u_{N} \\ \vdots & \vdots & \ddots & \vdots \\ -2u_{1}u_{N} & -2u_{2}u_{N} & \dots & 1-2u_{N}^2 \end{bmatrix}\\ \)

    はハウスホルダー行列。
     
  • ハウスホルダー行列は、正規直交行列
    それぞれのベクトルは直行する

    \(h_{i} \bullet h_{j} = - 2 u_{i} u_{j} (1 - 2 u_{i}^2) - 2 u_{j} u_{i} (1 - 2 u_{j}^2) + 4 u_{i} u_{j} \sum_{k=1}^{N(k \neq i,j)} u_k^2 \\ = 4 u_{i} u_{j} \sum_{k=1}^{N} u_k^2 - 4 u_{i} u_{j} \\ = 4 u_{i} u_{j} | \overrightarrow{u}| - 4 u_{i} u_{j} \\ = 4 u_{i} u_{j} - 4 u_{i} u_{j} = 0 \)

    ベクトルのノルムは、どれも１

    \(h_{i} \bullet h_{i} = (1 - 2 u_{i}^2)^2 + 4 u_{j}^2 \sum_{k=1}^{N(k \neq i)} u_k^2 \\ = 4 u_{i}^2 \sum_{k=1}^{N} u_k^2 - 4 u_{i}^2 + 1 \\ = 4 u_{i}^2 | \overrightarrow{u}| - 4 u_{i}^2 + 1 \\ = 4 u_{i}^2 - 4 u_{i}^2 + 1 = 1 \)

ハウスホールダー変換によるQR法 †

• Step 1

  \(A = \begin{bmatrix} a_{11} & a_{12} & a_{13} & \dots & a_{1N} \\ a_{21} & a_{22} & a_{23} & \dots & a_{2N} \\ a_{31} & a_{32} & a_{33} & \dots & a_{3N} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ a_{N1} & a_{N2} & a_{N3} & \dots & a_{NN} \\ \end{bmatrix} \)

  の１列目の２行目以降を 0 にする。
  ⇔

  \(\overrightarrow{a_1} = [a_{11},a_{21},a_{31},\dots,a_{N1}]^T \)

  をノルムが同じで、1行目以外が 0 の

  \(\overrightarrow{a'_1} = [\sqrt{a_{11}^2 + a_{21}^2 + a_{31}^2 + \dots + a_{N1}^2},0,0,\dots,0]^T \)

  にする。つまり、

  \(\overrightarrow{v_1} = \overrightarrow{a_1} - \overrightarrow{a'_1}\\ \overrightarrow{u_1} = \frac{\overrightarrow{v_1}}{|\overrightarrow{v_1}|} \)

  この u₁ を使って

  \(H(u_1) A = A' = \begin{bmatrix} a'_{11} & a'_{12} & a'_{13} & \dots & a'_{1N} \\ 0 & a'_{22} & a'_{23} & \dots & a'_{2N} \\ 0 & a'_{32} & a'_{33} & \dots & a'_{3N} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & a'_{N2} & a'_{N3} & \dots & a'_{NN} \\ \end{bmatrix} \)

• Step 2

  \(A' = \begin{bmatrix} a'_{11} & a'_{12} & a'_{13} & \dots & a'_{1N} \\ 0 & a'_{22} & a'_{23} & \dots & a'_{2N} \\ 0 & a'_{32} & a'_{33} & \dots & a'_{3N} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & a'_{N2} & a'_{N3} & \dots & a'_{NN} \\ \end{bmatrix} \)

  の２列目の３行目以降を 0 にする。
  ⇔

  \(\overrightarrow{a'_2} = [a'_{12},a'_{22},a'_{32},\dots,a'_{N2}]^T \)

  をノルムが同じで、1行目は変えずに、3行目以降が 0 の

  \(\overrightarrow{a''_2} = [a'_{12},\sqrt{{a'}_{22}^2 + {a'}_{32}^2 + \dots + {a'}_{N2}^2},0,\dots,0]^T \)

  にする。つまり、

  \(\overrightarrow{v_2} = \overrightarrow{a'_2} - \overrightarrow{a''_2}\\ \overrightarrow{u_2} = \frac{\overrightarrow{v_2}}{|\overrightarrow{v_2}|} \)

  この u₂ を使って

  \(H(u_2)H(u_1) A = H(u_2)A' = \begin{bmatrix} a''_{11} & a''_{12} & a''_{13} & \dots & a''_{1N} \\ 0 & a''_{22} & a''_{23} & \dots & a''_{2N} \\ 0 & 0 & a''_{33} & \dots & a''_{3N} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & a''_{N3} & \dots & a''_{NN} \\ \end{bmatrix} \)

  【補足】

  \(\overrightarrow{a'_2} \Rightarrow \overrightarrow{a''_2} \)

  の変換は、第一要素を変えないので、

  \(\overrightarrow{u_2} = [0,p,q,\dots,z]^T \)

  となるはず。すると

  \(H(u_2)A' = \begin{bmatrix} 1 & 0 & 0 & \dots & 0 \\ 0 & h_{22} & h_{23} & \dots & h_{2N} \\ 0 & h_{32} & h_{33} & \dots & h_{3N} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & h_{N2} & h_{N3} & \dots & h_{NN} \\ \end{bmatrix} \begin{bmatrix} a'_{11} & a'_{12} & a'_{13} & \dots & a'_{1N} \\ 0 & a'_{22} & a'_{23} & \dots & a'_{2N} \\ 0 & a'_{32} & a'_{33} & \dots & a'_{3N} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & a'_{N2} & a'_{N3} & \dots & a'_{NN} \\ \end{bmatrix} \)

  となって、Step 1 で整理した A' の 1 列目を変えない
• ...
• Step N
  最終的に

  \(H(u_N) \dots H(u_2)H(u_1) A = \begin{bmatrix} r_{11} & r_{12} & r_{13} & \dots & r_{1N} \\ 0 & r_{22} & r_{23} & \dots & r_{2N} \\ 0 & 0 & r_{33} & \dots & r_{3N} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & r_{NN} \\ \end{bmatrix} \\ \)

  と、上三角行列に変換できる。\(Q=H(u_N) \dots H(u_2)H(u_1)\) 最後にできる上三角行列を R とすると

  \(QA=R \\ A=Q^{-1}R \)

  Q は、直交行列の積なので直交行列だから \(Q^{-1} = Q^{T} (\because Q^{T}Q=I)\)

  \(A=Q^{T}R \)

  と表せる。これは QT 分解 //

ハウスホールダー変換によるヘッセンベルグ行列への変換 †

• 目的
  • QT分解を上三角行列に収束するまで何回もやるのは大変なので、ハウスホールダー変換でヘッセンベルグ行列まで相似変換してやる。
  • ヘッセンベルグ行列は、上三角行列の対角要素のひとつ下まで要素が入っている行列

    \(\begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} & a_{15} \\ a_{21} & a_{22} & a_{23} & a_{24} & a_{25} \\ 0 & a_{32} & a_{33} & a_{34} & a_{35} \\ 0 & 0 & a_{43} & a_{44} & a_{45} \\ 0 & 0 & 0 & a_{54} & a_{55} \end{bmatrix} \\ \)

  • 元が対称行列なら、三重対角行列になる

    \(\begin{bmatrix} a_{11} & a_{12} & 0 & 0 & 0 \\ a_{21} & a_{22} & a_{23} & 0 & 0 \\ 0 & a_{32} & a_{33} & a_{34} & 0 \\ 0 & 0 & a_{43} & a_{44} & a_{45} \\ 0 & 0 & 0 & a_{54} & a_{55} \end{bmatrix} \\ \)

  • ヘッセンベルグ行列または三重対角行列まで相似変換してからQR法を用いる
• 相似変換

  \(H(u)^{-1}AH(u) = A' \)

  は相似変換。H(u) は、正規直交行列だから

  \(H(u)=H(u)^T=H(u)^{-1} \)

  したがって

  \(H(u)^{-1}AH(u) = H(u)AH(u) = A' \)

• 対称行列の相似変換
  
   
  A が対称行列の時

  \(H(u)AH(u) = A'H(u) = A'' \)

  を計算する。

  \(\begin{bmatrix} u_{11} & u_{12} & u_{13} \\ u_{12} & u_{22} & u_{23} \\ u_{13} & u_{23} & u_{33} \\ \end{bmatrix} \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{12} & a_{22} & a_{23} \\ a_{13} & a_{23} & a_{33} \\ \end{bmatrix} = \begin{bmatrix} a'_{11} & a'_{12} & a'_{13} \\ a'_{21} & a'_{22} & a'_{23} \\ a'_{31} & a'_{32} & a'_{33} \end{bmatrix} \\ a'_{rc} = \sum_{i=1}^3 u_{ri}a_{ic} \\ \begin{bmatrix} a'_{11} & a'_{12} & a'_{13} \\ a'_{21} & a'_{22} & a'_{23} \\ a'_{31} & a'_{32} & a'_{33} \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & u_{13} \\ u_{21} & u_{22} & u_{23} \\ u_{31} & u_{32} & u_{33} \\ \end{bmatrix} = \begin{bmatrix} a''_{11} & a''_{12} & a''_{13} \\ a''_{21} & a''_{22} & a''_{23} \\ a''_{31} & a''_{32} & a''_{33} \end{bmatrix} \\ a''_{rc} = \sum_{j=1}^3 a'_{rj}u_{jc} = \sum_{j=1}^3 (\sum_{i=1}^3 u_{jc}u_{ri}a_{ij}) \\ a''_{cr} = \sum_{j=1}^3 a'_{cj}u_{jr} = \sum_{j=1}^3 (\sum_{i=1}^3 u_{jr}u_{ci}a_{ij}) \\ \)

  ここで、ハウスホルダー行列は対称行列なので、u_ij=u_ji、
  変換対象 A が対称行列の時の話なので a_ij=a_ji

  \(a''_{rc} = \sum_{j=1}^3 (\sum_{i=1}^3 u_{jc}u_{ri}a_{ij}) \\ a''_{cr} = \sum_{j=1}^3 (\sum_{i=1}^3 u_{ic}u_{rj}a_{ij}) \\ \therefore a''_{rc} = a''_{cr} \)

  つまり、対称行列をハウスホルダー行列 H(u) で相似変換すると、対称行列になる。
  
   
  ⇒ 行列 A を相似変換して左下が 0 なヘッセンベルグ行列にしようとした時、元の行列 A が対称行列なら自然に三重対角行列になってしまう

• Step 1
  A の 1 列目の 3 行目以降を 0 にするような相似変換を考える。

  \(H(u_1) = \left[ \begin{array}{c:c} 1 & \overrightarrow{0}^t \\ \hdashline % \overrightarrow{0} & Q \end{array} \right] \)

  とすると

  \(H(u_1) A H(u_1)= \left[ \begin{array}{c:c} 1 & \overrightarrow{0}^t \\ \hdashline % \overrightarrow{0} & Q \end{array} \right] \left[ \begin{array}{c:c} a_{11} & \overrightarrow{h}^t \\ \hdashline % \overrightarrow{v} & C \end{array} \right] \left[ \begin{array}{c:c} 1 & \overrightarrow{0}^t \\ \hdashline % \overrightarrow{0} & Q \end{array} \right] = \left[ \begin{array}{c:c} a_{11} & \overrightarrow{h}^t \\ \hdashline % Q\overrightarrow{v} & QC \end{array} \right] \left[ \begin{array}{c:c} 1 & \overrightarrow{0}^t \\ \hdashline % \overrightarrow{0} & Q \end{array} \right] = \left[ \begin{array}{c:c} a_{11} & \overrightarrow{h}^tQ \\ \hdashline % Q\overrightarrow{v} & QCQ \end{array} \right] \)

  つまり Qv の2行目以降が 0 になるような Q を見つけばよい
  ⇔
  

  \(\overrightarrow{a_1} = [a_{11},a_{21},a_{31},\dots,a_{N1}]^T \)

  をノルムが同じで、1行目は変えずに、3行目以降が 0 の

  \(\overrightarrow{a'_1} = [a_{11},\sqrt{a_{21}^2 + a_{31}^2 + \dots + a_{N1}^2},0,\dots,0]^T \)

  にするハウスホルダー変換。
  この変換で

  \(A_{2} = H(u_1) A_{1} H(u_1) = \left[ \begin{array}{cc:c} a_{11} & {a'}_{12} & \dots \\ {a'}_{21} & {a'}_{22} & \dots \\ \hdashline % \overrightarrow{0} & \overrightarrow{v_2} & C_{2} \end{array} \right] = \left[ \begin{array}{c:c} {\alpha}_{2} & H_{2} \\ \hdashline % V_{2} & C_{2} \end{array} \right] \)

  となる
• Step 2
  A の 2 列目の 3 行目以降を 0 にするような相似変換を考える。

  \(H(u_1) = \left[ \begin{array}{c:c} I & 0 \\ \hdashline % 0 & Q_{2} \end{array} \right] \)

  とすると

  \(H(u_2) A_2 H(u_2) = \left[ \begin{array}{c:c} I & 0 \\ \hdashline % 0 & Q_{2} \end{array} \right] \left[ \begin{array}{c:c} {\alpha}_{2} & H_{2} \\ \hdashline % V_{2} & C_{2} \end{array} \right] \left[ \begin{array}{c:c} I & 0 \\ \hdashline % 0 & Q_{2} \end{array} \right] = \left[ \begin{array}{c:c} {\alpha}_{2} & H_{2}Q_{2} \\ \hdashline % Q_{2}V_{2} & Q_{2}C_{2}Q_{2} \end{array} \right] \\ \\ Q_{2}V_{2} = Q_{2} \begin{bmatrix} \overrightarrow{0} & \overrightarrow{v_2} \end{bmatrix} = \begin{bmatrix} \overrightarrow{0} & Q_{2}\overrightarrow{v_2} \end{bmatrix} \)

  つまり Q₂V₂ の2行目以降が 0 になるような Q₂ を見つけばよい
  ⇔
  

  \(\overrightarrow{a''_2} = [a''_{12},a''_{22},a''_{32},\dots,a''_{N2}]^T \)

  をノルムが同じで、1行目、2行目は変えずに、4行目以降が 0 の

  \(\overrightarrow{a'''_2} = [a''_{12},a''_{22},\sqrt{a_{32}^2 + \dots + a_{N2}^2},0,\dots,0]^T \)

  にするハウスホルダー変換。
  この変換で

  \(A_{3} = H(u_2) A_{2} H(u_2) = \left[ \begin{array}{ccc:c} a_{11} & {a'}_{12} & {a''}_{13} & \dots \\ {a'}_{21} & {a'}_{22} & {a''}_{23} & \dots \\ 0 & {a''}_{32} & {a''}_{33} & \dots \\ \hdashline % \overrightarrow{0} & \overrightarrow{0} & \overrightarrow{v_3} & C_3 \end{array} \right] \)

  となる
• ・・・
• Step N-2
  最終的に

  \(H(u_{N-2}) A_{N-2} H(u_{N-2}) = \begin{bmatrix} h_{11} & h_{12} & h_{13} & \dots & h_{1~N-2} & h_{1~N-1} & h_{1N} \\ h_{21} & h_{22} & h_{23} & \dots & h_{2~N-2} & h_{2~N-1}& h_{2N} \\ 0 & h_{23} & h_{33} & \dots & h_{3~N-2} & h_{3~N-1}& h_{3N} \\ 0 & 0 & h_{43} & \dots & h_{4~N-2} & h_{4~N-1}& h_{4N} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & h_{N-3~N-2} & h_{N-3~N-1} & h_{N-3~N} \\ 0 & 0 & 0 & \dots & h_{N-2~N-2} & h_{N-2~N-1} & h_{N-2~N} \\ 0 & 0 & 0 & \dots & h_{N-1~N-2} & h_{N-1~N-1} & h_{N-1~N} \\ 0 & 0 & 0 & \dots & 0 & h_{N~N-1} & h_{N~N} \\ \end{bmatrix} \\ \)

  とヘッセンベルグ行列になる。

QR法が上三角行列に収束する証明 †

• Step 1
  • QR 法

    \(A_k = Q_k R_k \\ A_{k+1} = R_k Q_k \)

  • 記号の定義

    \(P_k = Q_1 Q_2 \dots Q_{k-1} Q_{k} \\ U_k = R_{k} R_{k-1} \dots R_2 R_1 \)

    ここで、Q は正則行列、R は上三角行列
• Step 2

  \(A_k = Q_k R_k = R_{k-1} Q_{k-1} \\ = (Q_{k-1}^{-1} Q_{k-1}) R_{k-1} Q_{k-1} \\ = Q_{k-1}^{-1} (Q_{k-1} R_{k-1}) Q_{k-1} \\ = Q_{k-1}^{-1} A_{k-1} Q_{k-1} \\ = (Q_1 Q_2 \dots Q_{k-1})^{-1} A (Q_1 Q_2 \dots Q_{k-1}) \\ = P_{k-1}^{-1} A P_{k-1} \\ \\ \therefore P_{k-1} A_k = A P_{k-1} \)

• Step3

  \(P_k U_k \\ = Q_1 Q_2 \dots Q_{k-1} (Q_{k} R_{k}) R_{k-1} \dots R_2 R_1 \\ = Q_1 Q_2 \dots Q_{k-1} A_k R_{k-1} \dots R_2 R_1 \\ = P_{k-1} A_k U_{k-1} \\ \\ \therefore P_k U_k = P_{k-1} A_k U_{k-1} \)

• Step4
  Step2,3より

  \(P_k U_k = A P_{k-1} U_{k-1} \\ P_2 U_2 = A P_1 U_1 = A Q_1 R_1 = A^2 \\ P_3 U_3 = A P_2 U_2 = A^3 \\ \vdots \\ P_k U_k = A^k \\ \\ \therefore A^k = Q_1 Q_2 \dots Q_{k-1} Q_{k} R_{k} R_{k-1} \dots R_2 R_1 \)

• Step 5
  話は変わって、A を対角化する X と対角化された行列 Λ を考える

  \(X^{-1}AX = \Lambda = \begin{bmatrix} \lambda_1 & 0 & \dots & 0 & 0 \\ 0 & \lambda_2 & \dots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & \lambda_{n-1} & 0 \\ 0 & 0 & \dots & 0 & \lambda_{n} \end{bmatrix} \\ ( \lambda_1 > \lambda_2 > \dots > \lambda_{n-1} > \lambda_{n} ) \)

  ここで、

  \(X^{-1}AX = \Lambda \\ X(X^{-1}AX)X^{-1} = X \Lambda X^{-1} \\ \therefore A = X \Lambda X^{-1} \)

  したがって

  \(A^k \\ = (X \Lambda X^{-1})(X \Lambda X^{-1}) \dots (X \Lambda X^{-1})(X \Lambda X^{-1}) = X \Lambda (X^{-1}X) \Lambda (X^{-1} \dots X) \Lambda (X^{-1}X) \Lambda X^{-1} = X \Lambda^k X^{-1} \)

  ここで、X の QR 分解と X^-1 の LU 分解を

  \(X = Q_X R_X \\ X^{-1} = L_X U_X \)

  とすると

  \(A^k \\ = X \Lambda^k X^{-1} \\ = Q_X R_X \Lambda^k L_X U_X \\ = Q_X R_X \Lambda^k L_X ((\Lambda^{-1})^k \Lambda^k) U_X \\ = Q_X R_X (\Lambda^k L_X (\Lambda^{-1})^k) (\Lambda^k U_X) \)

• Step 6
  Step 5 で求めた式の右辺の一部を計算する

  \(\Lambda^k L_X (\Lambda^{-1})^k \\ = \begin{bmatrix} \lambda_1^k & 0 & \dots & 0 & 0 \\ 0 & \lambda_2^k & \dots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & \lambda_{n-1}^k & 0 \\ 0 & 0 & \dots & 0 & \lambda_{n}^k \end{bmatrix} \begin{bmatrix} 1 & 0 & \dots & 0 & 0 \\ l_{ij} & 1 & \dots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ l_{ij} & l_{ij} & \dots & 1 & 0 \\ l_{ij} & l_{ij} & \dots & l_{ij} & 1 \end{bmatrix} \begin{bmatrix} \frac{1}{\lambda_1^k} & 0 & \dots & 0 & 0 \\ 0 & \frac{1}{\lambda_2^k} & \dots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & \frac{1}{\lambda_{n-1}^k} & 0 \\ 0 & 0 & \dots & 0 & \frac{1}{\lambda_n^k} \end{bmatrix} = \begin{bmatrix} 1 & 0 & \dots & 0 & 0 \\ l_{ij}(\frac{\lambda_i}{\lambda_j})^k & 1 & \dots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ l_{ij}(\frac{\lambda_i}{\lambda_j})^k & l_{ij}(\frac{\lambda_i}{\lambda_j})^k & \dots & 1 & 0 \\ l_{ij}(\frac{\lambda_i}{\lambda_j})^k & l_{ij}(\frac{\lambda_i}{\lambda_j})^k & \dots &l_{ij}(\frac{\lambda_i}{\lambda_j})^k & 1 \end{bmatrix} \)

  ここで k→∞ のとき λ₁ > λ₂ > … > λ_n-1 > λ_n だから、上式の左下 (行 i > 列 j) にでてくる

  \(\lim_{k \rightarrow \infty} (\frac{\lambda_i}{\lambda_j})^k = 0 \)

  従って k→∞ のとき

  \(\Lambda^k L_X (\Lambda^{-1})^k = I \)

• Step 7
  Step 6 より Step 5 の式は k→∞ のとき

  \(A^k = Q_X R_X (\Lambda^k L_X (\Lambda^{-1})^k) (\Lambda^k U_X) = Q_X R_X I \Lambda^k U_X = Q_X R_X \Lambda^k U_X \)

  ここで、\(R_X \Lambda^k U_X\) は上三角行列になり、QR分解は一意なので、Step4 で求めた A^k の QR 分解と比較すると

  \(P_k = Q_1 Q_2 \dots Q_{k-1} Q_{k} = Q_X \\ U_k = R_{k} R_{k-1} \dots R_2 R_1 = R_X \Lambda^k U_X \)

• Step 8
  ここで Step2 に立ち返って

  \(A_k = P_{k-1}^{-1} A P_{k-1} \\ = Q_X^{-1} A Q_X \\ = Q_X^{-1} X X^{-1} A X X^{-1} Q_X \\ = Q_X^{-1} X \Lambda X^{-1} Q_X \\ = Q_X^{-1} Q_X R_X \Lambda (Q_X R_X)^{-1} Q_X \\ = R_X \Lambda R_X^{-1} Q_X^{-1} Q_X \\ = R_X \Lambda R_X^{-1} \)

  ここで、R_X は上三角行列で、上三角行列の逆行列も上三角行列だから R_X^-1 も上三角行列なので、
  

  \(R_X R_X^{-1} = I \\ \begin{bmatrix} r_{11} & r_{12} & r_{13} & \dots \\ 0 & r_{22} & r_{23} & \dots \\ 0 & 0 & r_{33} & \dots \\ \vdots & \vdots & \vdots & \vdots \end{bmatrix} \begin{bmatrix} {r'}_{11} & {r'}_{12} & {r'}_{13} & \dots \\ 0 & {r'}_{22} & {r'}_{23} & \dots \\ 0 & 0 & {r'}_{33} & \dots \\ \vdots & \vdots & \vdots & \vdots \end{bmatrix} = \begin{bmatrix} r_{11}{r'}_{11} & {r''}_{12} & {r''}_{13} & \dots \\ 0 & r_{22}{r'}_{22} & {r''}_{23} & \dots \\ 0 & 0 & r_{33}{r'}_{33} & \dots \\ \vdots & \vdots & \vdots & \vdots \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & \dots \\ 0 & 1 & 0 & \dots \\ 0 & 0 & 1 & \dots \\ \vdots & \vdots & \vdots & \vdots \end{bmatrix} \\ \therefore r_{ii}{r'}_{ii} = 1 \)

  よって、

  \(A_k = R_X \Lambda R_X^{-1} = \\ \begin{bmatrix} r_{11} & r_{12} & r_{13} & \dots \\ 0 & r_{22} & r_{23} & \dots \\ 0 & 0 & r_{33} & \dots \\ \vdots & \vdots & \vdots & \vdots \end{bmatrix} \begin{bmatrix} \lambda_1 & 0 & 0 & \dots \\ 0 & \lambda_2 & 0 & \dots \\ 0 & 0 & \lambda_3 & \dots \\ \vdots & \vdots & \vdots & \vdots \end{bmatrix} \begin{bmatrix} {r'}_{11} & {r'}_{12} & {r'}_{13} & \dots \\ 0 & {r'}_{22} & {r'}_{23} & \dots \\ 0 & 0 & {r'}_{33} & \dots \\ \vdots & \vdots & \vdots & \vdots \end{bmatrix} \\ = \begin{bmatrix} \lambda_1 r_{11} & \lambda_2 r_{12} & \lambda_3 r_{13} & \dots \\ 0 & \lambda_2 r_{22} & \lambda_3 r_{23} & \dots \\ 0 & 0 & \lambda_3 r_{33} & \dots \\ \vdots & \vdots & \vdots & \vdots \end{bmatrix} \begin{bmatrix} {r'}_{11} & {r'}_{12} & {r'}_{13} & \dots \\ 0 & {r'}_{22} & {r'}_{23} & \dots \\ 0 & 0 & {r'}_{33} & \dots \\ \vdots & \vdots & \vdots & \vdots \end{bmatrix} \\ = \begin{bmatrix} \lambda_1 r_{11} {r'}_{11} & {r'''}_{12} & {r'''}_{13} & \dots \\ 0 & \lambda_2 r_{22} {r'}_{12} & {r'''}_{23} & \dots \\ 0 & 0 & \lambda_3 r_{33} {r'}_{33} & \dots \\ \vdots & \vdots & \vdots & \vdots \end{bmatrix} = \begin{bmatrix} \lambda_1 & {r'''}_{12} & {r'''}_{13} & \dots \\ 0 & \lambda_2 & {r'''}_{23} & \dots \\ 0 & 0 & \lambda_3 & \dots \\ \vdots & \vdots & \vdots & \vdots \end{bmatrix} \)

  つまり、A_k は k→∞ のとき上三角行列になる。また、対角成分には固有値が大きい順に左上から左下に向かって並ぶ。

Culture