JuliaのSchur分解をマスター：固有値問題からデータ解析まで活用
Julia言語のLinearAlgebra.Schurは、行列のシューア分解 (Schur decomposition) を計算するための機能を提供します。線形代数における重要な行列分解の一つで、特に非対称行列の固有値解析において非常に有用です。
シューア分解とは？
任意の正方行列 A は、以下の形に分解できます。
A=QTQ∗
ここで、
- T は上三角行列です。
- Q はユニタリ行列です。つまり、Q∗Q=I (I は単位行列) を満たします。実数の場合は、ユニタリ行列は直交行列 (QTQ=I) になります。
この分解をシューア分解と呼びます。
JuliaのLinearAlgebraモジュールに含まれるschur()関数は、このシューア分解を実行します。
使い方
基本的な使い方は以下の通りです。
using LinearAlgebra
A = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0] # 分解したい行列
S = schur(A)
schur(A)の結果として返されるSは、LinearAlgebra.Schur型のオブジェクトです。このオブジェクトは、以下のプロパティを持っています。
S.values: 行列 A の固有値のリスト (これは上三角行列 T の対角成分に対応します)S.Z: シューア分解によって得られたユニタリ行列 Q (Juliaでは慣習的にZと表記されることがあります)S.T: シューア分解によって得られた上三角行列 T
したがって、A≈S.Z∗S.T∗S.Z′ (ここで'は共役転置を表す) の関係が成り立ちます。
シューア分解のメリット
- 固有値の計算: 上三角行列 T の対角成分が元の行列 A の固有値になります。これは、固有値を数値的に安定して求める上で非常に重要な性質です。
- 安定性: 一般的に、非対称行列の固有値や固有ベクトルを直接計算するよりも、シューア分解を介して行う方が数値的に安定しています。特に、行列が対角化可能でない場合でも、シューア分解は常に存在します。
- 一般化固有値問題:
schur(A, B)のように2つの行列を与えると、一般化シューア分解 (generalized Schur decomposition) を計算することもできます。これは、Ax=λBx のような一般化固有値問題を解く際に利用されます。
実装の詳細
JuliaのLinearAlgebra.schurは、通常、高性能な数値線形代数ライブラリであるLAPACK (Linear Algebra PACKage) のルーチンを内部で利用しています。これにより、効率的かつ高精度な計算が実現されています。
例
using LinearAlgebra
# 複素数の行列の場合
A_complex = [1.0+2im 3.0-1im; 0.5+0.5im -2.0-3im]
S_complex = schur(A_complex)
println("上三角行列 T_complex:")
display(S_complex.T)
println("\nユニタリ行列 Z_complex:")
display(S_complex.Z)
println("\n固有値 values_complex:")
display(S_complex.values)
# 検証: Z * T * Z' が元の行列 A に戻ることを確認
println("\nZ * T * Z' の計算結果:")
display(S_complex.Z * S_complex.T * S_complex.Z')
println("\n元の行列 A_complex:")
display(A_complex)
入力行列の要件
schur関数は正方行列に対してのみ定義されます。非正方行列を与えるとエラーになります。
よくあるエラー
using LinearAlgebra
A = rand(3, 4) # 3x4 の非正方行列
schur(A)
エラーメッセージの例
DimensionMismatch("matrix must be square")
トラブルシューティング
入力行列が正方行列であることを確認してください。
データ型 (精度) の選択
schur関数は、入力行列の要素の型に基づいて計算の精度を決定します。
ComplexF64(倍精度複素数) は、複素行列の場合に使用します。Float32(単精度浮動小数点数) は、メモリ使用量を抑えたい場合や、精度がそれほど重要でない場合に利用できますが、精度問題が発生しやすくなります。Float64(倍精度浮動小数点数) が最も一般的で推奨されます。
よくある問題
Float32を使用した場合に、期待通りの精度が得られない、あるいは微小な数値誤差によって結果が不安定になることがあります。
トラブルシューティング
可能な限り、Float64またはComplexF64を使用することをお勧めします。特に、大規模な行列や、数値的に悪条件な行列を扱う場合は重要です。
# 精度が問題になる可能性のある例 (Float32)
A_float32 = Float32.(rand(3, 3))
S_float32 = schur(A_float32)
# 一般的には Float64 を使う
A_float64 = rand(3, 3)
S_float64 = schur(A_float64)
数値的な安定性 (悪条件な行列)
シューア分解は数値的に安定なアルゴリズムですが、元の行列が非常に「悪条件 (ill-conditioned)」である場合、結果の精度が低下したり、予期せぬ挙動を示すことがあります。悪条件とは、行列の微小な摂動が結果に大きな影響を与えることを意味します。
よくある問題
- 複数の固有値が非常に近い、または同じである行列: これらの固有値に対応する不変部分空間の計算が困難になることがあります。
- ほぼ特異な行列: 行列式がゼロに近い行列は、数値的に不安定になる傾向があります。
- 非常に大きな値と小さな値が混在する行列: スケーリングが適切でないと、小さい値が浮動小数点誤差で失われやすくなります。
トラブルシューティング
- 高精度演算の使用:
BigFloatを使用することで、より高い精度で計算できますが、計算コストが大幅に増加します。 - Perturbation (摂動): ごくわずかなランダムな摂動を加えることで、特定の特異ケースを回避できる場合があります。ただし、これは問題の性質を歪める可能性もあるため、慎重に行う必要があります。
- 条件数の確認:
cond(A)を使って行列の条件数を確認します。条件数が非常に大きい（例: 1015以上）場合、結果の信頼性は低くなる可能性があります。 - スケーリング: 必要に応じて、行列を事前にスケーリングして、要素の値の範囲を均一にすることを検討してください。
メモリ不足エラー (OutOfMemoryError)
大規模な行列のシューア分解は、大量のメモリを消費します。
よくあるエラー
OutOfMemoryError()
トラブルシューティング
- 計算環境の変更: より多くのメモリを持つマシンや、クラウドコンピューティングなどの環境を利用することを検討します。
- システムメモリの増強: 物理的なメモリが不足している場合は、増設を検討します。
- メモリ使用量の削減:
- 可能であれば、より小さいデータ型（
Float32）を使用する（ただし精度トレードオフに注意）。 - 行列が疎行列の場合、シューア分解は通常、密行列に変換されるため、疎行列の特性を活かせる別のアルゴリズム（例: 繰り返し法）を検討する必要があるかもしれません。
- 可能であれば、より小さいデータ型（
- 行列のサイズを確認: 処理している行列が大きすぎないか確認してください。
LAPACKからのエラー
LinearAlgebra.Schurは内部でBLAS/LAPACKライブラリを使用しているため、それらのライブラリからエラーが返されることがあります。エラーメッセージはしばしばLAPACKのルーチン名（例: _GEES、_GGEVなど）を含んでいます。
よくあるエラー
ArgumentError: invalid argument #X to LAPACK call
トラブルシューティング
これは、入力行列がLAPACKルーチンが期待する形式を満たしていない場合に発生することがあります。例えば、NaN（Not a Number）やInf（Infinity）などの非数値が含まれている場合です。
- NaN/Infの発生源を特定: これらの非数値がどこで生成されたかを追跡し、その計算プロセスを修正する必要があります。
- 入力行列の検査:
any(isnan, A)やany(isinf, A)を使って、行列に非数値が含まれていないか確認してください。
schur関数は、入力行列が実数であるか複素数であるかによって、異なる分解を行います。実数行列に対しては、通常、実シューア分解 (Real Schur decomposition) が計算され、対角ブロックに 2×2 のブロック（複素共役な固有値に対応）が含まれることがあります。複素行列に対しては、必ず上三角行列 T が得られます。
よくある混乱
実数行列を与えたのに、T が完全に上三角にならない（2×2 ブロックが現れる）ことに戸惑うことがあります。
トラブルシューティング
これは期待される挙動です。もし常に完全に上三角のTが必要な場合は、事前にComplexF64.(A)として行列を複素数型に変換してからschurを実行してください。
# 実数行列のシューア分解 (2x2ブロックが現れる可能性あり)
A_real = rand(4, 4)
S_real = schur(A_real)
display(S_real.T)
# 複素数行列として分解 (常に上三角)
A_complex = ComplexF64.(A_real)
S_complex = schur(A_complex)
display(S_complex.T)
例 1: 基本的なシューア分解
最も基本的な使い方です。実数行列のシューア分解を行います。
using LinearAlgebra
using Printf # 出力フォーマットのために使用
println("--- 例 1: 基本的なシューア分解 ---")
# 分解する実数行列を定義
A = [1.0 2.0 3.0;
4.0 5.0 6.0;
7.0 8.0 9.0]
println("元の行列 A:")
display(A)
# シューア分解を実行
S = schur(A)
println("\nシューア分解の結果:")
println("型: ", typeof(S)) # LinearAlgebra.Schur{Float64, Matrix{Float64}} のような型
# 各成分にアクセス
T = S.T # 上三角行列 (または準上三角行列)
Z = S.Z # ユニタリ行列 (または直交行列)
eigenvalues = S.values # 固有値
println("\n上三角行列 T (または準上三角行列):")
display(T)
println("\nユニタリ行列 Z (または直交行列):")
display(Z)
println("\n固有値 (S.values):")
display(eigenvalues)
# 検証: Z * T * Z' が A に戻ることを確認
# Z' は共役転置 (実数の場合は単なる転置)
A_reconstructed = Z * T * Z'
println("\nZ * T * Z' (再構築された行列 A):")
display(A_reconstructed)
# 誤差の確認
error_matrix = A - A_reconstructed
println("\n元の行列 A と再構築された行列の差 (誤差):")
display(error_matrix)
@printf "最大絶対誤差: %.2e\n" maximum(abs.(error_matrix))
解説
S.valuesは、行列 A の固有値のリストです。これらはTの対角要素、または 2×2 ブロックから計算される複素数となります。S.Zは直交行列であり、ZTZ=I を満たします。S.Tは、実数行列の場合、厳密な意味での上三角行列ではなく、対角ブロックに 2×2 のブロック（複素共役な固有値に対応）を持つ準上三角行列となることがあります。schur(A)はLinearAlgebra.Schur型のオブジェクトを返します。
例 2: 複素数行列のシューア分解
複素数行列を扱う場合、T は常に厳密な上三角行列になります。
using LinearAlgebra
using Printf
println("\n--- 例 2: 複素数行列のシューア分解 ---")
# 複素数行列を定義
B = [1.0+2.0im 3.0-1.0im;
0.5+0.5im -2.0-3.0im]
println("元の行列 B:")
display(B)
# シューア分解を実行
S_complex = schur(B)
T_complex = S_complex.T
Z_complex = S_complex.Z
eigenvalues_complex = S_complex.values
println("\n上三角行列 T_complex:")
display(T_complex)
println("\nユニタリ行列 Z_complex:")
display(Z_complex)
println("\n固有値 (S_complex.values):")
display(eigenvalues_complex)
# 検証: Z * T * Z' が B に戻ることを確認 (Z' は共役転置)
B_reconstructed = Z_complex * T_complex * Z_complex'
println("\nZ_complex * T_complex * Z_complex' (再構築された行列 B):")
display(B_reconstructed)
error_matrix_complex = B - B_reconstructed
println("\n元の行列 B と再構築された行列の差 (誤差):")
display(error_matrix_complex)
@printf "最大絶対誤差: %.2e\n" maximum(abs.(error_matrix_complex))
解説
S.Zはユニタリ行列であり、Z∗Z=I を満たします (Z∗ は共役転置)。- 複素数行列の場合、
schurは常にTを厳密な上三角行列として返します。
例 3: 特定の対角ブロックをソートする (schur with select argument)
シューア分解では、T の対角要素（固有値）の順序は保証されません。特定の固有値を特定の順序に並べ替えたい場合は、select 引数を使用できます。これは、固有値を並べ替えるための論理関数 (述語) を受け取ります。
using LinearAlgebra
using Printf
println("\n--- 例 3: 特定の対角ブロックをソートする ---")
# 例として、異なる固有値を持つ行列
C = [1.0 1.0; 0.0 2.0] # 固有値は 1.0 と 2.0
println("元の行列 C:")
display(C)
# デフォルトのシューア分解 (順序は不定)
S_default = schur(C)
println("\nデフォルトの T (S_default.T):")
display(S_default.T) # 固有値の順序はランダムになりうる
# 固有値を昇順にソートする例
# 固有値 `λ` を引数にとり、並べ替えたいかどうかを `true` または `false` で返す関数を定義
# ここでは、絶対値が小さい固有値を優先するようにソート
select_func_abs_asc = λ -> abs(λ) < 1.5 # 1.5 より小さい絶対値を持つ固有値を優先
S_sorted = schur(C, select_func_abs_asc)
println("\n絶対値が小さい固有値を優先してソートされた T (S_sorted.T):")
display(S_sorted.T)
println("ソートされた固有値 (S_sorted.values):")
display(S_sorted.values)
# 検証
C_reconstructed = S_sorted.Z * S_sorted.T * S_sorted.Z'
error_matrix_sorted = C - C_reconstructed
@printf "\n最大絶対誤差 (ソート後): %.2e\n" maximum(abs.(error_matrix_sorted))
# 固有値を降順にソートする例 (大きい固有値を優先)
select_func_desc = λ -> real(λ) > 1.5 # 実部が 1.5 より大きい固有値を優先
S_desc = schur(C, select_func_desc)
println("\n実部が大きい固有値を優先してソートされた T (S_desc.T):")
display(S_desc.T)
println("ソートされた固有値 (S_desc.values):")
display(S_desc.values)
解説
- この機能は、特定の固有値に対応する不変部分空間（invariant subspace）を抽出したい場合に非常に便利です。
schur(A, select_func)は、select_funcがtrueを返す固有値を、Tの左上に集めるように分解します。
2つの行列 A と B に対する一般化シューア分解は、Ax=λBx の形式の一般化固有値問題を解くために使用されます。分解は A=QSZ∗ と B=QTZ∗ の形になります。ここで、Q と Z はユニタリ行列、S と T は上三角行列です。
using LinearAlgebra
using Printf
println("\n--- 例 4: 一般化シューア分解 ---")
# 2つの行列 A と B を定義
A_gen = [1.0 2.0; 3.0 4.0]
B_gen = [5.0 6.0; 7.0 8.0]
println("元の行列 A_gen:")
display(A_gen)
println("\n元の行列 B_gen:")
display(B_gen)
# 一般化シューア分解を実行
GS = schur(A_gen, B_gen)
# 各成分にアクセス
S_gen = GS.S # 上三角行列 S
T_gen = GS.T # 上三角行列 T
Q_gen = GS.Q # ユニタリ行列 Q
Z_gen = GS.Z # ユニタリ行列 Z
generalized_eigenvalues = GS.alpha ./ GS.beta # 一般化固有値 (alpha / beta)
println("\n上三角行列 S_gen:")
display(S_gen)
println("\n上三角行列 T_gen:")
display(T_gen)
println("\nユニタリ行列 Q_gen:")
display(Q_gen)
println("\nユニタリ行列 Z_gen:")
display(Z_gen)
println("\n一般化固有値 (alpha ./ beta):")
display(generalized_eigenvalues)
# 検証: Q * S * Z' が A に戻ることを確認
A_reconstructed_gen = Q_gen * S_gen * Z_gen'
println("\nQ_gen * S_gen * Z_gen' (再構築された行列 A_gen):")
display(A_reconstructed_gen)
# 検証: Q * T * Z' が B に戻ることを確認
B_reconstructed_gen = Q_gen * T_gen * Z_gen'
println("\nQ_gen * T_gen * Z_gen' (再構築された行列 B_gen):")
display(B_reconstructed_gen)
error_A_gen = A_gen - A_reconstructed_gen
error_B_gen = B_gen - B_reconstructed_gen
@printf "\nA_gen の最大絶対誤差: %.2e\n" maximum(abs.(error_A_gen))
@printf "B_gen の最大絶対誤差: %.2e\n" maximum(abs.(error_B_gen))
GS.alphaとGS.betaは、一般化固有値 λ=α/β を定義するスカラーのベクトルです。betaがゼロに近い場合、対応する固有値は無限大に近づきます。GS.S,GS.T,GS.Q,GS.Zはそれぞれ、分解された行列の成分です。schur(A, B)はLinearAlgebra.GeneralizedSchur型のオブジェクトを返します。
固有値分解 (eigen / eigvals / eigvecs)
シューア分解と最も密接に関連しているのは、固有値分解 (Eigendecomposition) です。これは、行列 A を A=VΛV−1 の形に分解します。ここで、V は固有ベクトルを列にもつ行列、Λ は固有値を対角成分にもつ対角行列です。
LinearAlgebra.eigen 関数
- 関連関数
eigvals(A): 固有値のみを計算します。eigvecs(A): 固有ベクトルのみを計算します。
- 使用例
using LinearAlgebra A = [1.0 2.0; 3.0 4.0] E = eigen(A) println("--- 固有値分解 (eigen) ---") println("固有値 (E.values):") display(E.values) println("\n固有ベクトル行列 (E.vectors):") display(E.vectors) # 検証: A * E.vectors ≈ E.vectors * Diagonal(E.values) println("\nA * E.vectors:") display(A * E.vectors) println("\nE.vectors * Diagonal(E.values):") display(E.vectors * Diagonal(E.values)) - シューア分解との違い
- シューア分解は、行列が対角化可能でない場合でも常に存在し、結果として上三角行列（または準上三角行列）T を返します。
- 固有値分解は、行列が対角化可能な場合にのみ存在し、V が逆行列をもつことを保証します。
- 数値的安定性の観点から、
schurは非対称行列の固有値を安定して計算するためにしばしば推奨されます。特に、固有空間が複雑な場合（例：重複固有値がある場合など）。
- 目的
行列のすべての固有値と固有ベクトルを計算します。
特異値分解 (svd)
特異値分解 (Singular Value Decomposition, SVD) は、任意の m×n 行列 A を A=UΣV∗ の形に分解します。ここで、U と V はユニタリ行列、Σ は特異値を対角成分にもつ対角行列です。
- 使用例
using LinearAlgebra A = [1.0 2.0 3.0; 4.0 5.0 6.0] # 非正方行列 F = svd(A) println("\n--- 特異値分解 (svd) ---") println("ユニタリ行列 U (F.U):") display(F.U) println("\n特異値 (F.S):") display(F.S) println("\nユニタリ行列 V (F.V):") display(F.V) # 検証: U * Diagonal(S) * V' が A に戻ることを確認 A_reconstructed_svd = F.U * Diagonal(F.S) * F.V' println("\nU * Diagonal(S) * V' (再構築された行列 A):") display(A_reconstructed_svd) - シューア分解との違い
- SVDは正方行列に限定されず、任意の行列に適用できます。
- SVDの特異値は常に非負の実数です。固有値は複素数になることがあります。
- シューア分解は固有値と固有ベクトル（不変部分空間）に関連していますが、SVDは行列の「作用」を理解するのに役立ちます。
- 目的
行列のランク、条件数、近似などの情報を提供し、データ解析（主成分分析など）や連立一次方程式の解法に広く用いられます。
ヘッセンベルグ分解 (hessenberg)
ヘッセンベルグ分解は、正方行列 A を A=QHQ∗ の形に分解します。ここで、Q はユニタリ行列、H はヘッセンベルグ行列（準上三角行列に似ていますが、主対角線の下の最初の対角線以外はゼロ）です。
- 使用例
using LinearAlgebra A = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0] H = hessenberg(A) println("\n--- ヘッセンベルグ分解 (hessenberg) ---") println("ヘッセンベルグ行列 H (H.H):") display(H.H) println("\nユニタリ行列 Q (H.Q):") display(H.Q) # 検証: Q * H * Q' が A に戻ることを確認 A_reconstructed_hessenberg = H.Q * H.H * H.Q' println("\nQ * H * Q' (再構築された行列 A):") display(A_reconstructed_hessenberg) - シューア分解との違い
- ヘッセンベルグ分解はシューア分解よりも計算コストが低く、多くの固有値アルゴリズムの最初のステップとして利用されます。
- シューア分解は最終的な上三角形式を提供しますが、ヘッセンベルグ分解は中間形式です。
- 目的
シューア分解や固有値分解の予備ステップとして使用されることがよくあります。ヘッセンベルグ行列は疎な構造を持つため、固有値計算が効率的になります。
繰り返し法 (KrylovKit.jlなど)
LinearAlgebra.Schurやeigenは、行列のすべての固有値/固有ベクトル（またはシューア分解）を計算する「直接法」に分類されます。しかし、非常に大規模な疎行列の場合、これらの直接法はメモリや計算時間の点で非現実的になることがあります。このような場合、繰り返し法 (Iterative methods) が代替として用いられます。
-
他のパッケージ
Arpack.jl: 疎行列の固有値問題を解くための別の人気のあるパッケージです。
-
KrylovKit.jlのschursolve関数
KrylovKit.jlには、部分シューア分解を計算するschursolve関数もあります。これは、特定の数の極値固有値に対応するシューア分解を計算したい場合に便利です。using KrylovKit using LinearAlgebra A = rand(100, 100) # 大規模な行列を想定 # A = sprand(1000, 1000, 0.01) # 疎行列の例 println("\n--- KrylovKit.jl の schursolve (部分シューア分解) ---") # 最大絶対値を持つ3つの固有値に対応する部分シューア分解を計算 vals, vecs, info = schursolve(A, 3, :LM) # :LM は Largest Magnitude を意味する println("選択された固有値 (vals):") display(vals) println("\nシュールベクトル (vecs の最初の3つの列):") # vecs はシュールベクトル (Q) の列をベクトルとして持つリスト # 通常、この vecs は Q の最初のいくつかの列に対応し、完全な Q ではない display(hcat(vecs...)) # ベクトルのリストを結合して行列にする println("\n収束情報 (info):") display(info) -
シューア分解との違い
- 通常、
schurは密行列に対して最適化されていますが、繰り返し法は疎行列に優れています。 - 繰り返し法は近似解を返しますが、直接法は（浮動小数点精度内で）正確な解を返します。
- 通常、
-
Juliaのパッケージ
KrylovKit.jlは、Krylov部分空間法（アーノルディ法、ランチョス法など）を実装しており、疎行列や行列の明示的な形式が利用できない（行列-ベクトル積のみが利用可能）場合に特に有用です。 -
目的
行列のごく一部の固有値（例：最大または最小の固有値、特定の値に近い固有値）のみが必要な場合に効率的です。
- 固有値問題の前処理として、より効率的な計算が必要な場合
LinearAlgebra.hessenberg(ヘッセンベルグ分解) を使用します。 - 行列の構造解析、ランク、近似などが目的の場合
LinearAlgebra.svd(特異値分解) を使用します。 - 行列が大規模で疎行列であり、一部の固有値/固有ベクトルのみが必要な場合
KrylovKit.jlやArpack.jlなどの繰り返し法を使用します。 - すべての固有値と固有ベクトルが必要な場合
LinearAlgebra.eigenを使用します。 - 完全なシューア分解が必要な場合
LinearAlgebra.schurを使用します。