Julia ordschur() エラーとトラブルシューティング【解決策まとめ】

実シューア分解とは

まず、実シューア分解について簡単に説明します。正方行列 A の実シューア分解は、次のような形になります。

A=QTQT

ここで、Q は直交行列（つまり、QTQ=QQT=I、I は単位行列）であり、T は準三角行列（upper quasi-triangular matrix）です。準三角行列とは、対角成分が 1times1 のブロック（実数の固有値に対応）または 2times2 のブロック（複素共役な固有値のペアに対応）である上三角ブロック行列のことです。

LinearAlgebra.ordschur() 関数の役割

LinearAlgebra.ordschur() 関数は、標準の LinearAlgebra.schur() 関数によって得られる実シューア分解に加えて、固有値（または対応するシューアブロック）を指定した順序で並べ替える機能を提供します。これにより、特定の固有値に対応する不変部分空間（invariant subspace）を抽出したり、固有値の性質に基づいて行列の構造を分析したりすることが可能になります。

関数の使い方

LinearAlgebra.ordschur() 関数の基本的な使い方は以下の通りです。

using LinearAlgebra

A = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0]

# 固有値の順序を指定しない場合（これは LinearAlgebra.schur() と同様の結果になります）
F = schur(A)
println("シューア分解 T:\n", F.T)
println("直交行列 Q:\n", F.Q)

# 固有値の順序を指定する場合
# 順序付け関数（predicate function）を定義します。
# 例えば、実部が正の固有値を先頭に並べる場合：
positive_real_part(w) = real(w) > 0

# ordschur! 関数は行列を上書きするので、コピーを作成することが推奨されます。
Acopy = copy(A)
F_ordered = ordschur!(Acopy, select=positive_real_part)
println("\n順序付けられたシューア分解 T (実部が正の固有値が先頭):\n", F_ordered.T)
println("順序付けられた直交行列 Q:\n", F_ordered.Q)

# 固有値の絶対値の大きさで順序付ける例：
by_magnitude(w) = abs(w)

Acopy2 = copy(A)
F_ordered_magnitude = ordschur!(Acopy2, select=by_magnitude)
println("\n順序付けられたシューア分解 T (絶対値が大きい固有値が先頭):\n", F_ordered_magnitude.T)
println("順序付けられた直交行列 Q:\n", F_ordered_magnitude.Q)

引数と戻り値

• 戻り値
  

  • Schur 型のオブジェクト。このオブジェクトは以下のフィールドを持ちます。
    • T: 順序付けられた準三角シューア行列。
    • Q: 順序付けに対応する直交行列。
    • values: シューア行列の対角成分（1times1 ブロックの実数値）または 2times2 ブロックから計算される固有値。順序付けられています。
• • A: 分解したい正方行列。
  • select: 省略可能な引数で、固有値（またはシューアブロックに対応する代表値）に基づいて順序を決定するための関数（predicate function）です。この関数は、複素数または 2times2 の実数行列を受け取り、true または false を返します。true を返す固有値（またはブロック）が、シューア形式の先頭に移動します。select! という破壊的なバージョンも存在し、これは select 関数が状態を持つ場合に利用できます。

重要な点

• ordschur() の非破壊的なバージョンは、Julia の比較的新しいバージョンで導入された可能性があります。古いバージョンでは ordschur! のみが提供されている場合があります。
• select 関数は、固有値そのものだけでなく、2times2 の実数ブロックの代表値（例えば、ブロックの最初の固有値）に対しても適用されます。
• ordschur! は、入力行列 A を直接変更する破壊的な関数であることに注意してください。元の行列を保持したい場合は、事前に copy(A) などでコピーを作成してください。

---

引数の型の不一致

• トラブルシューティング
  
  • 入力行列が正方行列であることを確認してください (size(A, 1) == size(A, 2))。
  • 行列の要素の型を確認してください (eltype(A))。実シューア分解は実数または複素数の要素を持つ行列に対して定義されます。通常は Float64 が用いられます。必要に応じて convert(Matrix{Float64}, A) などで型を変換してください。
• 原因
  ordschur! (または ordschur) に渡された引数の型が期待される型と異なる場合に発生します。例えば、正方行列でない行列を渡したり、要素の型が Float64 などの浮動小数点数でない行列を渡したりするとエラーになります。
• エラー
  MethodError: no method matching ordschur!(...)

select 関数の定義ミス

• トラブルシューティング
  
  • select 関数が、LinearAlgebra.eigvals(A) の要素と同じ型の引数を受け取るように定義されているか確認してください。実数の固有値の場合は Real 型、複素数の固有値の場合は Complex 型です。2×2 ブロックの場合は、そのブロック自体が渡されます。
  • select 関数が必ず true または false の Bool 型の値を返すように定義されているか確認してください。
  • 関数のロジックが、意図した順序付けの条件を満たしているか再確認してください。テスト用の固有値やブロックで select 関数を手動で評価してみるのも有効です。
• 原因
  select 引数に渡す関数が、期待される引数（固有値または 2×2 の実数ブロック）を受け取らない、または Bool 型の値を返さない場合に発生します。
• エラー
  MethodError (具体的なエラーメッセージは select 関数の定義による)

ordschur! による元の行列の変更

• トラブルシューティング
  
  • 元の行列を保持したい場合は、ordschur! を呼び出す前に copy(A) で行列のコピーを作成し、そのコピーに対して ordschur! を適用してください。
  • 非破壊的な ordschur 関数が利用できる Julia のバージョンを使用している場合は、そちらを使用することを検討してください。
• 問題
  ordschur! は破壊的な関数であり、入力として与えた行列 A の内容が分解の結果で上書きされます。元の行列を保持しておきたい場合に、意図せず内容が変更されてしまうことがあります。

順序付けが期待通りにならない

• トラブルシューティング
  
  • select 関数の条件が、意図した固有値（またはブロック）を正しく識別しているか再度確認してください。
  • 複数の固有値が select 関数の条件を満たす場合、それらの相対的な順序は保証されないことがあります。より厳密な順序付けが必要な場合は、より複雑な select 関数を定義するか、分解後に結果をさらに処理することを検討してください。
  • 複素共役な固有値は、2×2 の実数ブロックとして扱われます。select 関数がこれらのブロックの代表値（通常はブロックの左上の要素の固有値）に基づいて評価されることを理解しておいてください。
• 問題
  select 関数を定義しても、固有値の順序が期待通りにならない場合があります。

数値的な安定性の問題

• トラブルシューティング
  
  • 入力行列の条件数 (cond(A)) を評価し、非常に大きい場合は注意が必要です。
  • 問題によっては、特異値分解 (SVD) など、より安定した他の分解方法を検討する必要があるかもしれません。
• 問題
  極端に条件数の悪い行列に対して ordschur! を適用すると、数値的な不安定性により精度が低下したり、エラーが発生したりする可能性があります。

Julia のバージョンによる違い

• トラブルシューティング
  
  • 使用している Julia のバージョンを確認し (versioninfo())、公式ドキュメントでそのバージョンにおける LinearAlgebra.ordschur の仕様を確認してください。
• 問題
  Julia のバージョンによっては、ordschur 関数の挙動や利用可能なオプションが異なる場合があります。

• ドキュメントを参照する
  Julia の公式ドキュメントは、関数の詳細な仕様や使用例を提供しています。
• 中間結果の確認
  分解の結果 (F.T, F.Q, F.values) を出力して、順序付けがどのように行われているかを確認します。
• 簡単な例で試す
  問題のある行列の代わりに、既知の固有値を持つ小さな行列で ordschur! を試してみることで、select 関数のロジックなどを検証できます。
• エラーメッセージをよく読む
  Julia のエラーメッセージは、問題の原因の手がかりとなる重要な情報を含んでいます。

---

例1: 基本的な実シューア分解と固有値の表示

この例では、まずランダムな実数行列を作成し、その実シューア分解を計算します。そして、得られたシューア行列の対角成分（または 2times2 ブロックから計算される固有値）を表示します。

using LinearAlgebra
using Random

# 4x4 のランダムな実数行列を作成
rng = MersenneTwister(1234) # 乱数生成器を固定して結果を再現可能にする
A = randn(rng, 4, 4)

# 実シューア分解を計算
F = schur(A)

println("元の行列 A:\n", A)
println("\nシューア行列 T:\n", F.T)
println("\n直交行列 Q:\n", F.Q)
println("\n固有値:\n", F.values)

このコードでは、schur() 関数を使って基本的な実シューア分解を行っています。ordschur() との違いは、この関数では固有値の順序付けは行われない点です。F.values には、シューア行列の対角ブロックから計算された固有値が格納されています。実数の固有値はそのまま、複素共役な固有値のペアは 2times2 のブロックに対応します。

例2: select 関数を用いた固有値の実部による順序付け

この例では、ordschur! 関数と select 関数を使って、固有値の実部が正であるものをシューア行列の先頭に移動させます。

using LinearAlgebra
using Random

rng = MersenneTwister(5678)
A = randn(rng, 4, 4)

# 固有値の実部が正であるかを判定する関数
is_positive_real_part(w) = real(w) > 0

# 行列のコピーを作成（ordschur! は破壊的な関数）
Acopy = copy(A)

# 順序付けられた実シューア分解を計算
F_ordered = ordschur!(Acopy, select=is_positive_real_part)

println("元の行列 A:\n", A)
println("\n順序付けられたシューア行列 T (実部が正の固有値が先頭):\n", F_ordered.T)
println("\n順序付けられた直交行列 Q:\n", F_ordered.Q)
println("\n順序付けられた固有値:\n", F_ordered.values)

ここでは、is_positive_real_part(w) という関数を select 引数に渡しています。この関数は、与えられた固有値 w の実部が正であれば true を、そうでなければ false を返します。ordschur! は、この条件を満たす固有値に対応するシューアブロックを先頭に移動するように分解を行います。

例3: select 関数を用いた固有値の絶対値による順序付け

この例では、固有値の絶対値が大きい順にシューア行列を並べ替えます。

using LinearAlgebra
using Random

rng = MersenneTwister(9012)
A = randn(rng, 4, 4)

# 固有値の絶対値が大きいかを判定する関数（ここでは常に true を返すようにして、
# 比較関数を select に渡す方法を示します）
function by_magnitude(w)
    return abs(w)
end

# 比較関数を select に直接渡すことはできません。
# select 関数は Bool 型を返す必要があります。
# したがって、少し工夫が必要です。

# eigenvalues 関数で固有値を事前に計算し、その順序に基づいて select 関数を作成します。
eigenvalues = eigvals(A)
sorted_indices = sortperm(eigenvalues, by=abs, rev=true) # 絶対値の大きい順にソートしたインデックス

function select_by_magnitude(w)
    # 与えられた固有値 w が、ソートされた固有値の先頭のグループに
    # 対応するかどうかを判定するロジックをここに記述する必要があります。
    # これは少し複雑になるため、より直接的な方法として、
    # ordschur! のドキュメントで推奨されている方法（predicate function）
    # に戻ります。

    # より簡単な例として、絶対値が 1 より大きい固有値を先頭にする場合：
    return abs(w) > 1.0
end

Acopy2 = copy(A)
F_ordered_magnitude = ordschur!(Acopy2, select=select_by_magnitude)

println("元の行列 A:\n", A)
println("\n順序付けられたシューア行列 T (絶対値が 1 より大きい固有値が先頭):\n", F_ordered_magnitude.T)
println("\n順序付けられた直交行列 Q:\n", F_ordered_magnitude.Q)
println("\n順序付けられた固有値:\n", F_ordered_magnitude.values)

この例では、固有値の絶対値に基づいて順序付けを行おうとしましたが、select 関数は Bool 型を返す必要があるため、直接的な比較関数を渡すことはできません。代わりに、絶対値が特定の閾値（ここでは 1.0）より大きい固有値に対応するブロックを先頭に移動する select_by_magnitude 関数を定義しています。

より複雑な順序付け（例えば、絶対値の大きい順に厳密に並べる）を行う場合は、ordschur! を直接使用するのではなく、まず eigvals() で固有値を計算し、それらをソートした順序に基づいて、適切な変換行列を構築し、元のシューア分解に適用するなどの追加のステップが必要になる場合があります。しかし、ordschur! の select 引数は、特定の性質を持つ固有値をグループ化するのに非常に便利です。

例4: 2times2 ブロックに対する select 関数の使用

実シューア分解では、複素共役な固有値のペアは 2times2 の実数ブロックとして現れます。select 関数は、これらの 2times2 ブロックに対しても適用されます。select 関数に渡されるのは、この 2times2 の行列です。

using LinearAlgebra

# 複素共役な固有値を持つ可能性のある 2x2 行列
A = [0.0 1.0; -2.0 0.0]

function select_2x2_block(block)
    # ここでは単純に、左上の要素が正である $2 \times 2$ ブロックを先頭に移動する例
    return block[1, 1] > 0
end

Acopy3 = copy(A)
F_ordered_block = ordschur!(Acopy3, select=select_2x2_block)

println("元の行列 A:\n", A)
println("\n順序付けられたシューア行列 T:\n", F_ordered_block.T)
println("\n直交行列 Q:\n", F_ordered_block.Q)
println("\n固有値:\n", F_ordered_block.values)

---

LinearAlgebra.schur() と手動での並べ替え

ordschur() の主な機能は、実シューア分解の計算と同時に、指定した条件に基づいて固有値（またはシューアブロック）を並べ替えることです。したがって、代替案として、まず LinearAlgebra.schur() で標準の実シューア分解を計算し、その後、得られたシューア行列 T と直交行列 Q を手動で並べ替えるという方法が考えられます。

このアプローチの利点は、並べ替えのロジックを完全に自分で制御できることです。例えば、複数の条件を組み合わせて複雑な順序付けを行いたい場合や、並べ替えの過程で追加の計算を行いたい場合に有効です。

using LinearAlgebra

A = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0]
F = schur(A)
T = F.T
Q = F.Q
eigenvalues = eigvals(A) # 必要に応じて固有値を計算

# ここで、T と Q を eigenvalues の順序に基づいて手動で並べ替えるロジックを実装します。
# これは、T の対角ブロックと Q の対応する列を入れ替える操作になります。

# (具体的な並べ替えのコードは、目的とする順序付けの条件によって大きく異なります。)

# 例: eigenvalues の実部が大きい順に並べ替える（概念的なコード）
sorted_indices = sortperm(real.(eigenvalues), rev=true)
T_ordered = zeros(size(T))
Q_ordered = zeros(size(Q))

current_row = 1
for index in sorted_indices
    # eigenvalues[index] に対応する T のブロックと Q の列を取り出し、
    # T_ordered と Q_ordered の適切な位置に配置する
    # (1x1 ブロックの場合と 2x2 ブロックの場合を考慮する必要があります)
    # ... (複雑な処理) ...
end

# 結果の T_ordered と Q_ordered を使用する
# println("手動で並べ替えた T:\n", T_ordered)
# println("手動で並べ替えた Q:\n", Q_ordered)

この方法は柔軟性が高い反面、実装が複雑になる可能性があります。特に、2times2 のブロックを扱う場合や、直交性の維持に注意が必要です。

固有値分解 (eigen() または eigvals()) に基づく方法

行列が対角化可能である場合（すべての固有ベクトルが線形独立である場合）、固有値分解 A=PDP−1 を利用することも考えられます。ここで、D は固有値を対角成分に持つ対角行列、P は対応する固有ベクトルを列とする行列です。

固有値分解を計算した後、固有値（D の対角成分）を任意の順序に並べ替え、それに対応するように固有ベクトル（P の列）も並べ替えることができます。ただし、一般の非対称行列の場合、P は直交行列とは限らないため、実シューア分解とは異なる分解になります。

using LinearAlgebra

A = [1.0 2.0; 3.0 1.0]
eigen_decomposition = eigen(A)
eigenvalues = eigen_decomposition.values
eigenvectors = eigen_decomposition.vectors

# 固有値を絶対値の大きい順に並べ替える
sorted_indices = sortperm(abs.(eigenvalues), rev=true)
eigenvalues_sorted = eigenvalues[sorted_indices]
eigenvectors_sorted = eigenvectors[:, sorted_indices]

# 並べ替えた固有値と固有ベクトルを使用して、必要に応じて行列を再構成できます。
# A_reconstructed = eigenvectors_sorted * Diagonal(eigenvalues_sorted) * inv(eigenvectors_sorted)

println("並べ替えた固有値:\n", eigenvalues_sorted)
println("並べ替えた固有ベクトル:\n", eigenvectors_sorted)

この方法は、対角化可能な行列に対しては比較的簡単ですが、実シューア分解とは異なる性質を持つため、目的によっては適切でない場合があります。また、非対称行列の場合、固有ベクトルは一般に直交しません。

外部ライブラリの利用

特定の高度な順序付けアルゴリズムや、数値的な安定性を重視する場合には、他の専門的なライブラリの利用を検討することもできます。Julia のエコシステムには、線形代数に関連する様々なパッケージが存在し、より特化した機能を提供している可能性があります。ただし、LinearAlgebra パッケージは Julia の標準ライブラリであり、多くの基本的な線形代数演算を効率的に実行できるため、通常はまず標準ライブラリの機能を検討することが推奨されます。

反復法による特定の固有空間の計算

特定の性質を持つ固有値に対応する不変部分空間（固有空間または一般固有空間）のみに関心がある場合は、反復法（例えば、べき乗法や逆べき乗法、部分空間反復法など）を適用して、直接それらの部分空間を計算する方法もあります。これらの方法は、行列全体を分解するよりも計算コストが低い場合があります。

using LinearAlgebra

A = [1.0 2.0; 3.0 1.0]

# 絶対値が最大の固有値と対応する固有ベクトルをべき乗法で求める（簡単な例）
function power_iteration(A, q0, num_iterations=100, tol=1e-6)
    q = q0 / norm(q0)
    λ_prev = Inf
    for i = 1:num_iterations
        v = A * q
        λ = dot(q, v)
        q = v / norm(v)
        if abs(λ - λ_prev) < tol
            return λ, q
        end
        λ_prev = λ
    end
    return λ, q
end

# 初期ベクトルをランダムに設定
q0 = randn(size(A, 1))
λ, v = power_iteration(A, q0)

println("絶対値が最大の固有値:", λ)
println("対応する固有ベクトル:", v)

この方法は、特定の固有値や固有空間に焦点を当てる場合に有効ですが、行列全体の構造を理解するには不向きです。

LinearAlgebra.ordschur() は、実シューア分解と固有値の順序付けを同時に行う便利な関数ですが、代替的な方法として以下のようなものが考えられます。

• 反復法
  特定の固有値や固有空間のみに関心がある場合に有効です。
• 外部ライブラリの利用
  特定の高度な機能が必要な場合に検討します。
• 固有値分解 (eigen())
  対角化可能な行列に対して有効ですが、実シューア分解とは異なります。
• LinearAlgebra.schur() と手動での並べ替え
  より柔軟な順序付けが可能ですが、実装が複雑になる場合があります。