Juliaのtgsen!()：パフォーマンス最適化とLAPACK直接利用のメリット
2025-05-27
一般化された固有値問題は、次の形式で表されます。 Ax=λBx
ここで、A と B は行列、x は一般化された固有ベクトル、λ は一般化された固有値です。
tgsen!()は、LAPACKのD/Z TGSEN (またはそのシングル精度版) に対応しており、以下の処理を実行します:
-
一般化されたシュール分解の計算: 入力された2つの行列 A と B に対して、一般化されたシュール分解 A=QSZH および B=QTZH を計算します。ここで、Q と Z はユニタリー行列（または直交行列）、S と T は上三角行列です。S と T の対角要素から一般化された固有値を得ることができます。
-
（オプション）シュール形式の並べ替え: 一般化された固有値を、ユーザーが指定した条件に基づいて並べ替えることができます。例えば、特定の一般化された固有値のペアをシュール行列の左上隅に移動させることで、部分空間の抽出や、特定のモードの解析に利用できます。
-
（オプション）一般化された固有ベクトルの計算: 必要であれば、一般化された固有ベクトルを行列 Q と Z の列として計算します。
-
（オプション）条件数の計算: 一般化された固有値や部分空間の条件数を計算することができます。これは、計算結果の精度や安定性を評価するのに役立ちます。
関数の名前にある!は、Juliaの慣習として、入力引数をインプレース（in-place）で変更する関数であることを示しています。つまり、tgsen!()は、入力された行列 A や B の内容を計算結果で上書きする可能性があります。
なぜLAPACKを直接呼び出すのか？
JuliaのLinearAlgebraモジュールは、多くの線形代数関数を提供していますが、それらの多くは内部的にLAPACKやBLASといった最適化された低レベルライブラリを呼び出しています。LinearAlgebra.LAPACK.tgsen!()のようにLAPACKの名前空間を直接使用する場合、通常は以下のような目的があります。
- 特定のLAPACK関数の機能を利用:
LinearAlgebraモジュールが高レベルなインターフェースとして提供していない、特定のLAPACK関数が持つユニークな機能を利用したい場合。 - パフォーマンスの最適化: 非常にパフォーマンスが重要な場面で、Juliaの高レベルなラッパーを介さずに、直接LAPACKの関数呼び出しを行うことで、わずかなオーバーヘッドを削減したい場合。
- より細かい制御: LAPACK関数が提供する詳細なオプション（例えば、どの情報を計算するか、どのように結果を並べ替えるかなど）にアクセスしたい場合。
tgsen!()はLAPACKのD/Z TGSENルーチンに対応しており、エラーが発生した場合は主に LAPACKException がスローされます。この例外には、LAPACKルーチンから返された INFO コードが含まれており、これがエラーの原因を特定する鍵となります。
LAPACKException(i): 負の INFO コード (i<0)
- トラブルシューティング:
- 引数の確認:
tgsen!()のドキュメントを丁寧に読み、すべての引数の型、次元、許容される値が正しいかを確認してください。特に、行列のサイズや文字引数の指定が正しいか注意深く見ます。 NaNやInfのチェック: 入力行列AとBにNaNやInfが含まれていないか確認します。any(isnan, A)やany(isinf, B)などでチェックできます。- ワークスペースの不足:
tgsen!()は内部的にワークスペースを必要としますが、Juliaのラッパーが通常は自動的に適切なサイズのワークスペースを割り当てます。しかし、非常に特殊なケースや、手動でワークスペースを渡す場合は、そのサイズが不足していないか確認します。
- 引数の確認:
- 考えられる原因:
- 不正な入力引数:
tgsen!()に渡された行列の次元が不適切、オプションの文字引数（例:jobvl,jobvrなど）が不正な文字である、ワークスペース配列のサイズが小さすぎる、といったケースです。 - Fortranのインデックスとのずれ: Juliaのインデックス（1から始まる）とFortranのインデックスが一致しない場合（ただし、
tgsen!はJuliaのラッパーなので、通常は開発者がこのずれを吸収しています）。 - 未定義の要素: 入力行列に
NaN(Not a Number) やInf(Infinity) が含まれている場合。
- 不正な入力引数:
- エラーの意味:
INFO = -iの場合、ルーチンの$i$番目の引数が不正な値を持っていたことを示します。LAPACKルーチンは、Fortranの呼び出し規約に基づいて引数を数えるため、1から数えられます。
- トラブルシューティング:
- 入力行列の性質を調べる:
- 行列の条件数 (
cond(A),cond(B)) を調べて、非常に大きい場合は数値的に不安定な問題を示唆します。 AとBが特異であるか (det(A),det(B)がゼロに近いか) を確認します。AとBの一般化された固有値 (eigvals(A, B)) を計算し、その分布や重複度を確認します。
- 行列の条件数 (
- スケール調整: 行列の要素が極端に大きい/小さい場合は、適切なスケール調整（例: 行列を定数で割るなど）を試みてください。
- 異なるアルゴリズムの検討:
tgsen!()が失敗する場合、より堅牢な（ただし、場合によっては遅い）別のアルゴリズムや、高レベルな関数（例:schur(A, B)）で対応できるか検討します。Juliaのschur(A, B)は、内部でより安定したLAPACKルーチンを呼び出している可能性があります。 select引数の調整: 固有値の並べ替えに問題がある場合は、select引数の指定を見直すか、より緩やかな条件を設定できないか検討します。- エラーメッセージの確認: スタックトレースと
LAPACKException(i)のiの値を確認し、LAPACKのドキュメントでそのINFOコードが何を意味するかを調べます（LAPACKのドキュメントはオンラインで公開されています）。
- 入力行列の性質を調べる:
- 考えられる原因:
- 特異な、あるいは数値的に不安定な行列: 入力行列 A や B が特異、または特異に近い（条件数が非常に大きい）場合、あるいは固有値が非常に近い値を持つ場合、数値的な問題によりアルゴリズムが収束しないことがあります。
- 極端な値の範囲: 行列の要素が非常に小さい値や非常に大きい値を持つ場合、浮動小数点数の精度限界に達し、計算が不安定になることがあります。
- 特定の固有値クラスターの並べ替えの問題:
select引数を使用して固有値を並べ替える際に、数値的に区別が難しい固有値クラスターを指定した場合、並べ替えがうまくいかないことがあります。 - バグ: 非常に稀ですが、LAPACKライブラリ自体のバグである可能性もゼロではありません（ただし、これは非常に安定したライブラリです）。
- エラーの意味:
INFO = iの場合、アルゴリズムが収束しなかったか、計算上の問題が発生したことを示します。具体的な意味はLAPACKルーチンによって異なりますが、tgsenの場合、一般化された固有値の計算や並べ替えの際に収束しなかったことを意味することが多いです。
- 特定のマシンでのみ発生する問題: まれに、特定のハードウェアやOS環境でのみ発生する問題があります。これは、数値演算の細かな差異や、BLAS/LAPACKの実装の違いに起因する可能性があります。
- トラブルシューティング: 別のマシンで同じコードを実行してみて、問題が再現するか確認します。Juliaのバージョンやパッケージの依存関係を統一して試すことも重要です。
- BLAS/LAPACKライブラリの設定問題: JuliaはデフォルトでOpenBLASを使用しますが、MKLなどの他のBLAS/LAPACKライブラリに切り替えることができます。ライブラリのバージョンや設定によっては、互換性の問題やパフォーマンスの問題が発生することがあります。
- トラブルシューティング: 使用しているBLAS/LAPACKライブラリが正しく設定されているか確認します。
versioninfo()でBLASとLAPACKのバージョンを確認できます。MKL.jlなどを使用している場合は、その設定も確認してください。Juliaのアップデートによって、BLAS/LAPACKライブラリのリンク方法が変更されることもあります。
- トラブルシューティング: 使用しているBLAS/LAPACKライブラリが正しく設定されているか確認します。
- メモリ不足: 非常に大きな行列を扱う場合、メモリ不足が発生する可能性があります。
- トラブルシューティング: より多くのメモリを確保するか、行列のサイズを小さくできないか検討します。大規模な問題では、疎行列を利用したソルバー（
Arpack.jlなど）を検討することも有効です。
- トラブルシューティング: より多くのメモリを確保するか、行列のサイズを小さくできないか検討します。大規模な問題では、疎行列を利用したソルバー（
トラブルシューティングの一般的なヒント:
- LAPACKドキュメントの参照: 最も確実なのは、エラーメッセージの
INFOコードに対応するLAPACKルーチン（D/Z TGSEN）の公式ドキュメントを参照することです。これにより、エラーの具体的な意味と、それに対するLAPACK開発者の推奨する対処法がわかります。 - 型に注意:
tgsen!()のような低レベルな関数は、入力データの型（Float64、ComplexF64など）に非常に敏感です。意図しない型が渡されていないか確認します。 - Julia Discourseなどのコミュニティの活用: Juliaの公式フォーラム（Julia Discourse）では、多くのユーザーが同様の問題に遭遇しており、解決策やアドバイスが得られる可能性があります。エラーメッセージやスタックトレースを正確に共有することで、より的確な助言が得られます。
- 最小限の再現可能な例 (MWE) の作成: エラーが発生するコードの中で、問題の箇所だけを切り出して最小限のコードにすることで、原因の特定がしやすくなります。
tgsen!()の主な役割は、以下の一般化された固有値問題を解く際に使用される、行列ペア (A,B) の一般化されたシュール分解 (Generalized Schur Decomposition) を計算し、さらに操作することです。
AX=lambdaBX
ここで、A と B は正方行列、X は一般化された固有ベクトル、$ \lambda $ は一般化された固有値です。
以下に、LinearAlgebra.LAPACK.tgsen!()の使用例をいくつか示します。
例1: 基本的な一般化されたシュール分解の計算
この例では、2つの行列 A と B を定義し、tgsen!()を使って一般化されたシュール分解を計算します。
using LinearAlgebra
# 行列 A と B を定義
A = [1.0 2.0; 3.0 4.0]
B = [5.0 6.0; 7.0 8.0]
# tgsen! は引数をインプレースで変更するため、コピーを作成
A_copy = copy(A)
B_copy = copy(B)
# tgsen! の呼び出し
# この例では、jobq='N', jobz='N' で Q と Z は計算しない
# select は全て false (並べ替えなし)
# m, pl, pr, dif は出力引数で、初期化が必要
# work, lwork, iwork, liwork はワークスペース
# INFO は LAPACK の戻り値
m = Ref{BlasInt}(0)
pl = Ref{Float64}(0.0)
pr = Ref{Float64}(0.0)
dif = Vector{Float64}(undef, 2) # dif(1) = difu, dif(2) = difl
# LAPACK のドキュメントに基づいて、適切なワークスペースサイズを決定する必要がある
# 通常、lwork = -1 を指定して最初の呼び出しで最適な lwork を取得し、
# その後、そのサイズでワークスペースを確保して再度呼び出す
# 簡単のためにここでは大きめのサイズで初期化
lwork = BlasInt(max(1, 4 * size(A, 1) + 16)) # 推奨サイズ
work = Vector{Float64}(undef, lwork)
liwork = BlasInt(max(1, size(A, 1) + 6)) # 推奨サイズ
iwork = Vector{BlasInt}(undef, liwork)
info = Ref{BlasInt}(0)
# 一般化されたシュール分解を計算
# 返り値: INFO
try
LinearAlgebra.LAPACK.tgsen!(
'N', # jobq: Q を計算するか ('N': 計算しない)
'N', # jobz: Z を計算するか ('N': 計算しない)
fill(false, size(A, 1)), # select: 固有値を並べ替えるための論理配列
size(A, 1), # n: 行列の次数
A_copy, # A: 入力行列 (シュール形式 S に上書きされる)
size(A, 1), # lda: A の先行次元
B_copy, # B: 入力行列 (シュール形式 T に上書きされる)
size(B, 1), # ldb: B の先行次元
zeros(ComplexF64, size(A,1)), # alpha: 一般化された固有値の分子部 (複素数)
zeros(Float664, size(A,1)), # beta: 一般化された固有値の分母部 (実数)
zeros(Float64, size(A, 1), size(A, 1)), # Q: ユニタリー変換行列 Q (jobq='V'の場合に更新)
size(A, 1), # ldq: Q の先行次元
zeros(Float64, size(A, 1), size(A, 1)), # Z: ユニタリー変換行列 Z (jobz='V'の場合に更新)
size(A, 1), # ldz: Z の先行次元
m, # m: 選択された固有値の数
pl, # pl: 選択された部分空間の射影のノルム
pr, # pr: 選択された部分空間の射影のノルム
dif, # dif: 選択された部分空間と残りの部分空間との間の距離 (条件数)
work, lwork, # work, lwork: ワークスペース
iwork, liwork, # iwork, liwork: 整数型ワークスペース
info # info: LAPACK の戻り値
)
if info[] != 0
@warn "tgsen! returned info = $(info[]). This might indicate a numerical issue or an invalid argument."
end
println("Original A:\n", A)
println("Original B:\n", B)
println("\nGeneralized Schur Form S (from A_copy):\n", A_copy)
println("Generalized Schur Form T (from B_copy):\n", B_copy)
println("Selected eigenvalues count (m): ", m[])
println("Condition number (dif): ", dif)
# 一般化された固有値の計算
alpha = zeros(ComplexF64, size(A,1))
beta = zeros(Float64, size(A,1))
# tgsen!の引数リストを再確認して、alphaとbetaが返されるように修正
# LAPACKのCTGSEN/DTGSENの仕様では、alpha/betaは呼び出し時に引数として渡され、
# その配列に結果が書き込まれる
# 上記のtgsen!呼び出しはalpha, betaの部分がゼロ配列になっていたので修正
# JuliaのLinearAlgebra.LAPACK.tgsen!の実際の引数リストを確認
# Juliaのドキュメントに直接の記載がない場合、LAPACKのFortranルーチンのドキュメントを参照
# 通常、JuliaのLAPACKインターフェースは、これらの引数を直接提供します。
# ここでは仮に引数に追加して記述しますが、実際のJuliaのソースコードで引数リストを確認することが重要です。
# (Re-run tgsen! with alpha/beta as output arguments if not already done)
# Since alpha and beta were passed as zeros(ComplexF64, N) and zeros(Float64, N)
# they should be updated by tgsen! if the routine itself computes them.
# The example above is conceptually correct in showing where they would be passed.
# To demonstrate, let's manually compute them from S and T:
println("\nGeneralized Eigenvalues (λ = S[i,i] / T[i,i]):")
for i in 1:size(A, 1)
# S[i,i] と T[i,i] は A_copy[i,i] と B_copy[i,i]
# 複素数対応のため ComplexF64 にキャスト
eigenvalue = ComplexF64(A_copy[i,i]) / B_copy[i,i]
println("λ_$i = ", eigenvalue)
end
catch e
if isa(e, LAPACKException)
println("LAPACKException occurred: INFO = $(e.info)")
else
rethrow(e)
end
end
解説:
infoは、LAPACKルーチンからの戻り値で、0であれば成功、負の値は引数のエラー、正の値は収束しないなどの問題を示します。work,lwork,iwork,liworkは、計算に必要なワークスペースです。LAPACK関数の多くは、lwork = -1を指定して一度呼び出すことで、最適なワークスペースサイズをwork[1]に返します。その後、そのサイズでワークスペースを確保して再度呼び出すのが一般的な使用法です。上記の例では簡単のため、推奨される大きさに初期化しています。pl,pr,difは、選択された部分空間の条件数に関連する情報です。dif[1](difu) は一般化された固有ベクトル空間の条件数、dif[2](difl) は一般化された固有空間の条件数を示します。mは、select配列によって選択された固有値の数です。alpha,betaは、計算された一般化された固有値の分子と分母を格納する配列です。一般化された固有値は lambda=alpha/beta となります。lda,ldb,ldq,ldzは、それぞれの行列の先行次元（leading dimension）です。Juliaでは通常、行列の行数と同じです。nは行列の次元です。jobqとjobzは、ユニタリー変換行列QとZを計算するかどうかを制御します。'N'は計算しない、'V'は計算するという意味です。tgsen!()は、入力行列AとBをインプレースで一般化されたシュール形式SとTに変換します。そのため、元の行列を保持したい場合はcopy()で複製してから渡します。
この例では、jobq='V'とjobz='V'を指定して変換行列QとZも計算し、さらにselect関数を使って特定の固有値を並べ替えます。
using LinearAlgebra
# 行列 A と B を定義
A = [1.0 2.0; 3.0 4.0]
B = [5.0 6.0; 7.0 8.0]
# tgsen! はインプレースで変更するためコピーを作成
A_s = copy(A)
B_s = copy(B)
Q_mat = Matrix{Float64}(I, size(A, 1), size(A, 1)) # 単位行列で初期化
Z_mat = Matrix{Float64}(I, size(A, 1), size(A, 1)) # 単位行列で初期化
N = size(A, 1)
# 一般化された固有値を並べ替えるための選択関数
# 例えば、実部の絶対値が0.5より大きい固有値を選択
# (alpha_r, alpha_i, beta) の形式で受け取る
select_func = (αr, αi, β) -> (abs(αr / β) > 0.5)
# select 配列の初期化（この配列が select_func によって埋められる）
# LAPACK tgsen の select 引数は、論理配列（Bool[]）として渡され、
# 内部的にこの配列に基づいて並べ替えが実行される。
# select_func は Julia の tgsen の高レベルインターフェースで使うものであり、
# LAPACK.tgsen! を直接使う場合は、Fortranのコールバック関数や
# 内部で select_func に基づいて予め select 論理配列を生成しておく必要がある。
# ここでは例として select_arr を手動で設定します。
# 例えば、一つ目の固有値を選択する場合:
select_arr = fill(false, N)
# 真の固有値を知っている場合に、その固有値に対応するインデックスをtrueにする
# tgsen! を実行する前に、A と B の一般化された固有値を計算して、
# どの固有値を "選択" するかを決定するのが現実的
# ここでは、例として最初の固有値を選択すると仮定
select_arr[1] = true
m = Ref{BlasInt}(0)
pl = Ref{Float64}(0.0)
pr = Ref{Float64}(0.0)
dif = Vector{Float64}(undef, 2)
lwork = BlasInt(max(1, 4 * N + 16))
work = Vector{Float64}(undef, lwork)
liwork = BlasInt(max(1, N + 6))
iwork = Vector{BlasInt}(undef, liwork)
info = Ref{BlasInt}(0)
alpha_out = zeros(ComplexF64, N) # 出力用 alpha
beta_out = zeros(Float64, N) # 出力用 beta
try
LinearAlgebra.LAPACK.tgsen!(
'V', # jobq: Q を計算する
'V', # jobz: Z を計算する
select_arr, # select: 固有値を並べ替えるための論理配列
N, # n: 行列の次数
A_s, # A: 入力行列 (S に上書き)
N, # lda
B_s, # B: 入力行列 (T に上書き)
N, # ldb
alpha_out, # alpha: 一般化された固有値の分子部
beta_out, # beta: 一般化された固有値の分母部
Q_mat, # Q: ユニタリー変換行列 Q
N, # ldq
Z_mat, # Z: ユニタリー変換行列 Z
N, # ldz
m, # m: 選択された固有値の数
pl, # pl
pr, # pr
dif, # dif
work, lwork,
iwork, liwork,
info
)
if info[] != 0
@warn "tgsen! returned info = $(info[]). This might indicate a numerical issue or an invalid argument."
end
println("\n--- With Q, Z, and selection ---")
println("Original A:\n", A)
println("Original B:\n", B)
println("\nGeneralized Schur Form S (from A_s):\n", A_s)
println("Generalized Schur Form T (from B_s):\n", B_s)
println("\nQ (left transformation matrix):\n", Q_mat)
println("Z (right transformation matrix):\n", Z_mat)
println("Selected eigenvalues count (m): ", m[])
println("Condition number (dif): ", dif)
# 選択された固有値
println("\nCalculated Generalized Eigenvalues (alpha/beta):")
for i in 1:N
eigenvalue = alpha_out[i] / beta_out[i]
println("λ_$i = ", eigenvalue)
end
# 検証: A = Q * S * Z' および B = Q * T * Z'
println("\nVerification (Q * S * Z'):")
println("A_reconstructed = \n", Q_mat * A_s * Z_mat')
println("B_reconstructed = \n", Q_mat * B_s * Z_mat')
println("\nOriginal A compared to reconstructed A (should be close):")
println(A ≈ Q_mat * A_s * Z_mat') # 約等しいか確認
println("\nOriginal B compared to reconstructed B (should be close):")
println(B ≈ Q_mat * B_s * Z_mat') # 約等しいか確認
catch e
if isa(e, LAPACKException)
println("LAPACKException occurred: INFO = $(e.info)")
else
rethrow(e)
end
end
重要な注意点:
- LAPACKException: LAPACKルーチンは、エラーが発生すると
infoを0以外の値に設定します。Juliaでは、これがLAPACKExceptionとしてスローされます。トラブルシューティングのために、この例外を捕捉し、e.infoの値を確認することが重要です。 select引数: LAPACKのtgsenのselect引数は、並べ替えたい固有値のブロックに対応するブール配列です。JuliaのLinearAlgebra.LAPACK.tgsen!を直接呼び出す場合、このselect配列を事前に適切に設定する必要があります。これは、通常、ユーザーが事前に一般化された固有値を計算し、どの固有値をグループ化したいかを決定してから、それに応じてselect配列を構築することを意味します。- 複素数行列: 上記の例は実数行列 (
Float64) 用ですが、複素数行列 (ComplexF64) を扱う場合は、対応するLAPACKの複素数ルーチン（Z TGSEN）が呼び出されます。引数の型は自動的に推論されますが、alphaのような出力配列はComplexF64で初期化する必要があります。 - ワークスペースの管理: 上記の例では、
lworkとliworkを固定の大きさに初期化していますが、これは常に最適な方法ではありません。理想的には、最初にlwork = -1(またはliwork = -1) でtgsen!()を呼び出し、work[1](またはiwork[1]) に返される最適なワークスペースサイズを取得し、そのサイズで配列を割り当ててから再度tgsen!()を呼び出すべきです。これにより、メモリの効率的な使用とパフォーマンスの最適化が図られます。
以下に、tgsen!() の代替となる主なプログラミング手法と、それぞれの利用シナリオを説明します。
LinearAlgebra.schur() および LinearAlgebra.schurfact() (推奨)
LinearAlgebra モジュールには、一般化されたシュール分解を計算するための高レベルな関数が用意されています。
-
schur!(A, B): これはschur(A, B)のインプレース版です。入力行列AとBを直接シュール形式に変換し、メモリ割り当てを減らすことができます。using LinearAlgebra A = [1.0 2.0; 3.0 4.0] B = [5.0 6.0; 7.0 8.0] # インプレース版を呼び出すため、AとBが変更される F = schur!(A, B) println("\n--- schur!() (in-place) ---") println("Schur form S (F.S, now original A):\n", A) # AはF.Sと同じ内容になる println("Schur form T (F.T, now original B):\n", B) # BはF.Tと同じ内容になる println("Left transformation matrix Q (F.Q):\n", F.Q) println("Right transformation matrix Z (F.Z):\n", F.Z)メリット: メモリ効率が高い（特に大規模行列の場合）。 注意: 元の行列
AとBは破壊的に変更されます。
LinearAlgebra.eigvals() (一般化された固有値のみが必要な場合)
もし一般化されたシュール分解自体は必要なく、単に一般化された固有値だけが必要な場合は、eigvals(A, B) が最も簡潔な方法です。
using LinearAlgebra
A = [1.0 2.0; 3.0 4.0]
B = [5.0 6.0; 7.0 8.0]
# 一般化された固有値を直接計算
gen_eigenvalues = eigvals(A, B)
println("\nGeneralized Eigenvalues using eigvals(A, B):\n", gen_eigenvalues)
メリット:
- 効率的: 必要最小限の計算を実行します。
- シンプル: 最もシンプルなインターフェースで、必要な情報（固有値）のみを返します。
LinearAlgebra.eigen() (一般化された固有値と固有ベクトルが必要な場合)
一般化された固有値とそれに対応する一般化された固有ベクトルの両方が必要な場合は、eigen(A, B) を使用します。これも Eigen オブジェクトを返します。
using LinearAlgebra
A = [1.0 2.0; 3.0 4.0]
B = [5.0 6.0; 7.0 8.0]
# 一般化された固有値と固有ベクトルを計算
E = eigen(A, B)
println("\n--- Generalized Eigen-decomposition using eigen(A, B) ---")
println("Eigenvalues (E.values):\n", E.values)
println("Eigenvectors (E.vectors):\n", E.vectors) # 行列の各列が固有ベクトル
# 検証: A * E.vectors ≈ B * E.vectors * Diagonal(E.values)
println("\nVerification (A * X ≈ B * X * Λ):")
println("Left side (A * E.vectors):\n", A * E.vectors)
println("Right side (B * E.vectors * Diagonal(E.values)):\n", B * E.vectors * Diagonal(E.values))
println("Are they approximately equal?: ", (A * E.vectors) ≈ (B * E.vectors * Diagonal(E.values)))
メリット:
- 使いやすい:
tgsen!()のような低レベルな複雑さを避けることができます。
上記の関数は密行列（Dense Matrix）に対して非常に効率的ですが、非常に大規模な疎行列（Sparse Matrix）の場合には、より特化したパッケージを使用することを検討します。
-
Arpack.jl: 疎行列の一部の（最大値、最小値など、端に位置する）固有値と固有ベクトルを計算するためのパッケージです。Krylov法の実装が含まれています。# using Pkg; Pkg.add("Arpack") using SparseArrays, LinearAlgebra, Arpack n = 1000 A_sparse = sprand(n, n, 0.01) # 疎行列A B_sparse = I + sprand(n, n, 0.01) # 疎行列B (対角要素を持つ) # 最大マグニチュードの5つの一般化された固有値と固有ベクトルを計算 # 実際には、LinearMaps.jl などと組み合わせて線形作用素として渡すことが多い # eigs は標準固有値問題向けなので、一般化には少し工夫が必要 # 例: Ax = λBx => B⁻¹Ax = λx # ただし、B⁻¹ を明示的に計算するのは非推奨なので、線形作用素で表現 # 厳密には、Arpack.eigs は直接 Ax = λBx 形式をサポートしていない場合があるため、 # JacobiDavidson.jl のような専門のパッケージがより適しています。 # しかし、簡単な例として、B が逆行列を持つ場合に標準問題に変換する例を示します。 # 通常は A*x = λ*B*x の形式で、AとBを直接渡すAPIを探すか、 # 内部的にシフト＆反転変換を適用して標準問題に変換する形で行われます。 # 例えば、Bが正定値の場合、Cholesky分解 B = L*L' を用いて # (L⁻¹ * A * L⁻') * y = λ * y, x = L⁻' * y # という標準固有値問題に変換することが可能です。 # ここでは、簡略化のため、BがinvertibleでAとBが密行列に変換されると仮定した場合の # 概念的な例にとどめます。 # eigs(A, B) のような直接のAPIは、通常の `eigvals` や `eigen` と異なり、 # `Arpack.eigs` では直接提供されていないことが多いです。 # 通常は、A と B から線形作用素を構成して渡します。 # 例: `eigs(A, B, nev=k, which=:LM)` は LinearAlgebra.eigen の疎行列版とは異なる。 # Sparse-specific libraries often implement their own generalized eigenvalue solvers.
ほとんどの場合、LinearAlgebra.LAPACK.tgsen!() を直接使用するよりも、LinearAlgebra.schur(A, B) または LinearAlgebra.eigen(A, B) を使用することを強く推奨します。これらはより安全で使いやすく、かつパフォーマンスも最適化されています。