Rで線形代数:ジョルダン分解

行列の各要素

• 行列 J (ジョルダン標準形):
  • J は、対角線上にジョルダンブロックと呼ばれる正方行列が並んだブロック対角行列です。
  • ジョルダンブロックとは、以下のような形の行列です。
    • 対角成分には、ある一つの固有値 λ が並びます。
    • そのすぐ一つ上の対角線（超対角）には 1 が並びます。
    • それ以外の要素はすべて 0 です。

  ジョルダンブロックの例:

  • 1×1 ブロック: [λ] (固有値そのもの)
  • 2×2 ブロック: [[λ, 1], [0, λ]]
  • 3×3 ブロック: [[λ, 1, 0], [0, λ, 1], [0, 0, λ]]

  ジョルダン行列 J の例:

J = [ [λ₁,  1,  0 ],  [ 0,  0,  0 ],  [ 0,  0,  0 ] ]
    [ [ 0, λ₁,  1 ],  [ 0,  0,  0 ],  [ 0,  0,  0 ] ]
    [ [ 0,  0,  λ₁],  [ 0,  0,  0 ],  [ 0,  0,  0 ] ]
    [ [ 0,  0,  0 ],  [λ₂,  1     ],  [ 0,  0,  0 ] ]
    [ [ 0,  0,  0 ],  [ 0,  λ₂    ],  [ 0,  0,  0 ] ]
    [ [ 0,  0,  0 ],  [ 0,  0,  0 ],  [λ₃         ] ]

この例では、固有値 λ₁ に対応する 3×3 のジョルダンブロック、λ₂ に対応する 2×2 のブロック、λ₃ に対応する 1×1 のブロックが並んでいます。

なお、もし行列が対角化可能であれば、すべてのジョルダンブロックが 1×1 となり、J は対角行列 D と一致します。

• 行列 P (一般化固有ベクトル行列):
  • P の列は、通常の固有ベクトルだけでなく、一般化固有ベクトルを含みます。
  • 一般化固有ベクトルとは、(A - λI)v = 0 を満たす通常の固有ベクトル v を拡張したもので、(A - λI)ᵏ v = 0 となるようなベクトルです。これらは「鎖（チェーン）」を形成し、その鎖の長さがジョルダンブロックのサイズに対応します。

1. テスト用の非対角化可能行列を定義

A <- matrix(c(4, 1, -1, 2), nrow = 2, byrow = TRUE)

# A - λI = [[4, 1], [-1, 2]] - λ * [[1, 0], [0, 1]]
#        = [[4, 1], [-1, 2]] - [[λ, 0], [0, λ]]
#        = [[4 - λ,   1   ],
#           [ -1  , 2 - λ]]

# det(A - λI) = (4 - λ) * (2 - λ) - (1) * (-1)
#             = (8 - 4λ - 2λ + λ²) + 1
#             = λ² - 6λ + 8 + 1
#             = λ² - 6λ + 9
#             = (λ - 3)(λ - 3) = 0
lambda <- 3

cat("--- ジョルダン分解の計算プロセス実演 ---\n")
cat("元の行列 A:\n")
print(A)
cat("重複する固有値 λ:", lambda, "\n\n")

--- ジョルダン分解の計算プロセス実演 ---
元の行列 A:
     [,1] [,2]
[1,]    4    1
[2,]   -1    2
重複する固有値 λ: 3 

2. 固有ベクトルを求める

I <- diag(1, 2)
A_minus_lambda_I <- A - lambda * I
cat("A - λI:\n")
print(A_minus_lambda_I)

A - λI:
     [,1] [,2]
[1,]    1    1
[2,]   -1   -1

方程式 (A - λI)v₁ = 0 を満たす固有ベクトル v₁ を探します((A - λI)の零空間の基底として)。

v₁ を [x, y] とすると、方程式は以下のようになります。

[ 1  1 ] [x] = [0]
[-1 -1 ] [y]   [0]

これを連立方程式として書き下すと、

1.   1*x + 1*y = 0 →  x + y = 0
2.  -1*x - 1*y = 0 → -x - y = 0

となります。

この2つの式は、式2に-1を掛けると式1になるため、本質的に全く同じです。

したがって、解くべき方程式は、実質的に x + y = 0 のみとなります。

この式は、y = -x と書き換えることができますので、つまり「固有ベクトル v₁ = [x, y] の2番目の要素 y は、1番目の要素 x の符号を反転させたものである」となります。

この条件を満たす 方向が同じ（または真逆）のベクトルは

• x = 1 なら v₁ = [1, -1]
• x = -1 なら v₁ = [-1, 1]
• x = 10 なら v₁ = [10, -10]
• x = 0.707 なら v₁ = [0.707, -0.707]

と、無数に存在します が、ここでは以降の計算を単純にするために、

v₁ = c(-1, 1)

とします。

もちろん v₁ = c(1, -1) や v₁ = c(2, -2) 等を選んで計算を続けても、最終的には正しいジョルダン分解にたどり着きます（ただし、途中の v₂ や P の数値は変わります）。

v1 <- matrix(c(-1, 1), ncol = 1)
cat("固有ベクトル v1 (手動設定):\n")
print(v1)

固有ベクトル v1 (手動設定):
     [,1]
[1,]   -1
[2,]    1

3. 一般化固有ベクトルを求める

(A - λI)v2 = v1 を満たす v2 を探します。

[ 1  1 ] [x] = [-1]
[-1 -1 ] [y]   [ 1]

連立方程式の形に書き直しますと

1. ( 1 * x) +  (1 * y)
2. (-1 * x) + (-1 * y)

この連立方程式の結果が、右辺のベクトル [-1, 1] と等しくならなければなりません。

1.  1x + 1y = -1 →  x + y = -1
2. -1x - 1y = 1  → -x - y = 1

式2の両辺に -1 を掛けますと 式1と全く同じになりますので、制約条件は x + y = -1 のみです。

よって、連立方程式は解が無数に存在する劣決定系であるため、solve() や qr.solve() ではエラーを返します。

そこで、ここでは連立方程式満たす解の一つを選択し、x = 0 とすると y = -1 となるため、v2 = [0, -1] を採用します。

v2 <- matrix(c(0, -1), ncol = 1)
cat("一般化固有ベクトル v2 (手動選択):\n")
print(v2)

# 検算: (A - λI)v2 が v1 と等しいか確認
cat("\n検算 (A - λI)v2:\n")
print(A_minus_lambda_I %*% v2)
cat("v1と等しいか:", all.equal(A_minus_lambda_I %*% v2, v1), "\n")

一般化固有ベクトル v2 (手動選択):
     [,1]
[1,]    0
[2,]   -1

検算 (A - λI)v2:
     [,1]
[1,]   -1
[2,]    1
v1と等しいか: TRUE 

4. 行列 P と J を組み立てる

# Pは固有ベクトルと一般化固有ベクトルを列に持つ
P <- cbind(v1, v2)
J <- matrix(c(lambda, 1, 0, lambda), nrow = 2, byrow = TRUE)

cat("--- 分解結果 ---\n")
cat("遷移行列 P = [v1, v2]:\n")
print(P)
cat("\nジョルダン標準形 J:\n")
print(J)

--- 分解結果 ---
遷移行列 P = [v1, v2]:
     [,1] [,2]
[1,]   -1    0
[2,]    1   -1

ジョルダン標準形 J:
     [,1] [,2]
[1,]    3    1
[2,]    0    3

5. 最終検算: A = P J P⁻¹ が成り立つか確認

# Pはv1とv2が線形独立であるため正則であり、solve()が使える
P_inv <- solve(P)
A_reconstructed <- P %*% J %*% P_inv

cat("--- 最終検算 ---\n")
cat("P %*% J %*% P⁻¹ の結果:\n")
print(A_reconstructed)
cat("\n復元された行列は元の行列と等しいか:", all.equal(A, A_reconstructed), "\n")

--- 最終検算 ---
P %*% J %*% P⁻¹ の結果:
     [,1] [,2]
[1,]    4    1
[2,]   -1    2

復元された行列は元の行列と等しいか: TRUE 

以上です。