linear-solvers

線形ソルバー

目標

数値シミュレーションで生じる線形系に対して、ソルバーを選択し、条件の悪さを評価し、収束を診断するための汎用的なワークフローを提供します。

要件

• Python 3.10+
• NumPy、SciPy（行列演算用）
• 依存関係については個々のスクリプトを参照してください

収集する入力

入力	説明	例
行列サイズ	系の次元	n = 1000000
疎性	非ゼロ要素の割合	0.01%
対称性	A = Aᵀ ですか？	はい
正定値性	A は正定値ですか？	はい（SPD）
条件	推定条件数	10⁶

選択の指針

ソルバー選択フローチャート

Is matrix dense and small enough to factor in memory (dense float64
storage n²·8 bytes < ~2 GB, i.e. n ≲ 16000)?
├── YES → Use direct solver (Cholesky/LDLᵀ/LU by symmetry)
└── NO → Is matrix symmetric?
    ├── YES → Is it positive definite?
    │   ├── YES → Use CG with AMG/IC preconditioner
    │   └── NO → Use MINRES
    └── NO → Is it nearly symmetric?
        ├── YES → Use BiCGSTAB
        └── NO → Use GMRES with ILU/AMG

クイックリファレンス

行列の種類	ソルバー	前処理
SPD、疎行列	CG	AMG、IC
対称不定値	MINRES	SPD 前処理（SSOR、対称ブロック対角、または SPD 部分に対する AMG）
非対称	GMRES、BiCGSTAB	ILU、AMG
密行列	LU、Cholesky	なし
鞍点問題	Schur 補行列、Uzawa	ブロック前処理

スクリプトの出力（JSON フィールド）

スクリプト	主な出力
scripts/solver_selector.py	recommended、alternatives、notes
scripts/convergence_diagnostics.py	rate、asymptotic_rate、stagnation、recommended_action
scripts/sparsity_stats.py	nnz、density、bandwidth、symmetry
scripts/preconditioner_advisor.py	suggested、notes
scripts/scaling_equilibration.py	row_scale、col_scale、notes
scripts/residual_norms.py	residual_norms、relative_norms、converged

ワークフロー

1. 行列の特性を把握する - 対称性、正定値性、疎性
2. 疎性を分析する - scripts/sparsity_stats.py を実行します
3. ソルバーを選択する - scripts/solver_selector.py を実行します
4. 前処理を選択する - scripts/preconditioner_advisor.py を実行します
5. スケーリングを適用する - 条件が悪い場合は scripts/scaling_equilibration.py を使用します
6. 収束を監視する - scripts/convergence_diagnostics.py を使用します
7. 問題を診断する - scripts/residual_norms.py で残差の履歴を確認します

対話形式のワークフロー例

ユーザー：GMRES ソルバーが 50 回の反復後に停滞しています。残差は 1e-3 まで減少しますが、その後は改善しません。

エージェントのワークフロー：

1. 収束を診断します：

   python3 scripts/convergence_diagnostics.py --residuals 1,0.1,0.01,0.005,0.003,0.002,0.002,0.002 --json
   

2. 前処理に関する推奨事項を確認します：

   python3 scripts/preconditioner_advisor.py --matrix-type nonsymmetric --sparse --ill-conditioned --json
   

3. 推奨事項：リスタートパラメーターを増やす、高い k を指定した ILU(k) を試す、または AMG に切り替えます。

解く前のチェックリスト

• 行列の対称性／正定値性を確認する
• サイズと疎性に基づいて直接法と反復法のどちらを使うか決める
• 物理スケールに応じて残差の許容誤差を設定する
• 行列構造に適した前処理を選択する
• 必要に応じてスケーリング／平衡化を適用する
• 収束を追跡し、停滞が起きた場合は調整する

CLI の例

# Analyze sparsity pattern
python3 scripts/sparsity_stats.py --matrix A.npy --json

# Select solver for SPD sparse system
python3 scripts/solver_selector.py --symmetric --positive-definite --sparse --size 1000000 --json

# Get preconditioner recommendation
python3 scripts/preconditioner_advisor.py --matrix-type spd --sparse --json

# Diagnose convergence from residual history
python3 scripts/convergence_diagnostics.py --residuals 1,0.2,0.05,0.01 --json

# Apply scaling
python3 scripts/scaling_equilibration.py --matrix A.npy --symmetric --json

# Compute residual norms
python3 scripts/residual_norms.py --residual 1,0.1,0.01 --rhs 1,0,0 --json

エラー処理

エラー	原因	解決方法
Matrix file not found	無効なパス	ファイルが存在することを確認します
Matrix must be square	正方行列ではない入力	行列の次元を確認します
Residuals must be positive	無効な残差データ	入力形式を確認します

解釈の指針

収束率

convergence_diagnostics.py は 2 つの収束率を報告します。rate は全履歴にわたる各反復の残差比の平均で、asymptotic_rate は直近の短い区間における平均です。stagnation フラグは asymptotic_rate（> 0.95）に基づきます。停滞は履歴の末尾に関する性質であり、初期の急速な減少によって末尾の停滞が見えなくなることがあるためです。以下の状態を判断する際は、asymptotic_rate を参照してください。

漸近収束率	意味	対応
< 0.1	優秀	現在の設定が最適です
0.1 - 0.5	良好	ほとんどの問題で許容範囲です
0.5 - 0.95	低速	より良い前処理を検討してください
> 0.95	停滞	ソルバーまたは前処理を変更してください

停滞の診断

パターン	考えられる原因	対処法
残差が横ばい	前処理が不十分	前処理を改善する
振動する	ほぼ特異、または不定	行列を確認し、別のソルバーを試す
減衰が非常に遅い	悪条件	スケーリングを適用し、AMG を使用する

検証チェックリスト

記録された値を使って以下の項目をすべて確認するまでは、求解結果を信頼しないでください。「問題なさそう」という判断では不十分です。

• convergence_diagnostics.py の asymptotic_rate を記録し、停滞しきい値の 0.95 未満であること（理想的には < 0.5）を確認しました。履歴全体の rate が低いだけでは、末尾の横ばいは否定できません。
• residual_norms.py の相対残差を、絶対ノルムだけでなく、物理スケールに合わせた --rel-tol（デフォルトは 1e-6）と照合しました。スケーリングされていない RHS では、過小な rhs が収束を装えないように --require-both を指定しました。
• solver_selector.py の recommended が、sparsity_stats.py で記録された実際の行列特性（symmetry、既知であれば正定値性）と一致することを確認しました。たとえば、CG は対称かつ正定値の場合に限り、対称不定値には MINRES、非対称には GMRES/BiCGSTAB を使用します。
• 悪条件と判定された系については、scaling_equilibration.py を実行し、row_scale_max/row_scale_min と col_scale_max/col_scale_min を記録しました。対称行列では対称性を保つために --symmetric（D A D）を指定し、非対称の両側スケーリングでは row_scale を適用してから col_scale を適用しました。
• scaling_equilibration.py の sparsity_stats.py の notes/zero_rows/zero_cols を確認しました。ゼロ行またはゼロ列がある場合、その系は構造的に特異であり、スケール 1 へのフォールバックでは解決しません。
• preconditioner_advisor.py の前処理が、選択した Krylov 法に適用可能であることを確認しました。特に、MINRES の前処理は SPD でなければならず、不定値の不完全 LDLᵀ は無効です。

よくある落とし穴と誤った理由付け

安易な近道	誤りである理由／対処法
「平均 rate が低いから収束しました。」	rate は履歴全体の平均で、序盤の急速な低下に大きく左右されます。停滞は末尾の特性です。asymptotic_rate を確認し、0.95 未満であることを確かめてください。
「絶対残差が小さいから完了です。」	RHS が大きい場合やスケーリングされていない場合、絶対ノルムが小さくても意味がないことがあります。relative_norms / relative_value を、物理スケールに合わせた --rel-tol と照合してください。
「対称だから CG を使えばよい。」	CG には対称性と正定値性の両方が必要です。対称不定値行列には、SPD 前処理を用いた MINRES が必要です。CG を使うと破綻または停滞する可能性があります。選択前に正定値性を確認してください。
「大規模な系だから、直接分解すればよい。」	solver_selector.py は、密な float64 ストレージ（n²·8 bytes < ~2 GB、n ≈ 16384）の場合に密行列直接ソルバーを選択します。それを超えると密な Cholesky/LU は実行不可能なため、反復法を選択する必要があります。
「スケーリングとは、各行をその最大値で割るだけです。」	片側の行スケーリングでは均等化されません。非対称行列では、行スケーリング後の行列から col_scale を導出し、両方を適用してください。対称行列では対称 D A D スケーリングを使わないと、対称性が失われます。
「GMRES が停滞するので、反復回数を増やせばよい。」	末尾が横ばいなら、問題は反復回数ではなく前処理または再始動長です。前処理を強化する（ILU のフィルレベルを上げる／AMG を使用する）、再始動パラメーターを増やす、または別の手法に切り替えてください。

セキュリティ

入力の検証

• すべての数値入力（残差、許容誤差、行列要素）が有限数として検証されます
• カンマ区切りの残差／ベクトル入力は 100,000 件を上限とします
• solver_selector.py の --size パラメーターの上限は 100 億です
• --matrix-type は固定の許可リスト（spd、symmetric-indefinite、nonsymmetric）に照らして検証されます
• ブール値フラグ（--symmetric、--positive-definite、--sparse、--ill-conditioned）は、型安全な argparse フラグです

ファイルアクセス

• sparsity_stats.py と scaling_equilibration.py は、--matrix で指定された単一の行列ファイル（.npy 形式）を読み込みます
• 作成された .npy ファイルを介した任意コードの実行を防ぐため、np.load() は allow_pickle=False を指定して呼び出されます
• 行列ファイルは、解析の前にサイズが 500 MB を超えている場合、読み込みを拒否されます
• 行列の各次元に上限（100,000）を設け、メモリ枯渇を防ぎます
• その他のすべてのスクリプトは外部ファイルを読み込まず、入力は CLI 引数で指定します

ツールの制限

• Read：スクリプトのソース、参照資料、行列ファイルの確認に使用します
• Write：解析結果やソルバーの推奨事項の保存に使用します。書き込み先はユーザーの作業ディレクトリ内に限定されます
• Grep/Glob：関連ファイルの特定や参照箇所の検索に使用します
• このスキルの allowed-tools では Bash を除外しており、信頼できない行列ファイルや数値入力の処理中にエージェントが任意のコマンドを実行するのを防ぎます

安全対策

• eval()、exec()、動的コード生成は使用しません
• すべてのサブプロセス呼び出しで明示的な引数リストを使用し、shell=True は使用しません
• ツールの利用範囲を制限し（Bash なし）、エージェントが読み取り／書き込み操作のみを行えるようにしています
• JSON 出力モードでは、シェルで解釈される内容を含まない、構造化された解析可能な結果が生成されます

制限事項

• 大規模な密行列：直接法ソルバーではメモリ不足が発生する場合があります
• 非常に不定な行列：標準的な前処理では失敗する場合があります
• 鞍点問題：専用のブロック前処理が必要です

参考資料

• references/solver_decision_tree.md - 選択ロジック
• references/preconditioner_catalog.md - 前処理の選択肢
• references/convergence_patterns.md - 失敗の診断
• references/scaling_guidelines.md - 平衡化の指針

バージョン履歴

• v1.2.0（2026-06-23）：漸近的な停滞検出、密行列の実行可能性に基づくソルバー選択、鞍点問題／小規模密行列への直接法ソルバーの振り分け、両側の平衡化スケーリング、CG の反復回数上限表、ドキュメントと評価の整合性を修正
• v1.1.0（2024-12-24）：ドキュメント、選択ガイダンス、例を強化
• v1.0.0：6 つのソルバー解析スクリプトを含む初回リリース