WOLFRAM言語チュートリアル

硬さ検出

概要

多くの微分方程式には,刻み幅を制限し,それにより陽的な解法の効率をも制限するある種の硬さというものが存在する.

この問題を回避するために,多数の陰的メソッドが長年に渡り開発されてきた.

同じ刻み幅の場合,陰的メソッドは内的線形代数に関連するオーバーヘッドのために,陽的メソッドよりも効率がはるかに劣ることがある.

このコストは,ある特定の領域において陰的メソッドはかなり大きい刻み幅を取ることができるということにより埋め合せられる.

ランタイム時に硬さを自動的に検出し,必要に応じて適切なメソッドに切り替えるユーザフレンドリーなコードを提供するため,数々の試みが行われてきた.

ここでは,コードに硬さ検出を自動的に装備するために提案されてきた多数のストラテジーを紹介する.

硬さ検出をどのようにNDSolveに実装するかを記述するために,行列の優勢固有値の推定問題に特に注目してみる.

ストラテジーの効率は数値例題で例示する.

初期化

定義済みの例題とユーティリティ関数を含むパッケージをロードする.
In[1]:=
Click for copyable input

はじめに

初期値問題の数値解を考える:

硬さとは問題,解法,初期条件,局所的誤差許容度の組合せである.

硬さはメソッドが取ることのできる刻み幅を制約するため,陽的解法の効率を制約する.

硬さは線の方法による偏微分方程式の数値解法だけでなく,多数の実践的な系で発生する.

例題

ファン・デル・ポル(van der Pol)振動子は,非線形減衰を受けた非保存系の振動子であり,常微分方程式の硬い系の例である:

このとき ;

初期条件

を考慮し,区間 上で解く.

メソッドではデフォルトで1ペアの外挿メソッドが使われる.

  • 陽的修正中点法(グラッグ法),二重調和数列2, 4, 6,
  • 線形陰的オイラー法,部分調和数列2, 3, 4,

以下のようにしてパッケージから問題をロードする.
In[4]:=
Click for copyable input
硬くないメソッドを使って系を数値的に解く.
硬さが発生したら切り替わるメソッドを使って系を解く.
In[6]:=
Click for copyable input
次は2つの解の成分のプロットである.尖った頂点(青)は大きさほぼ450くらいまで広がっているので切り取られている.
In[7]:=
Click for copyable input
Out[7]=

硬さは急速な過渡解に従う領域で起ることが多い.

以下は時間に対する刻み幅をプロットしたものである.
In[8]:=
Click for copyable input
Out[8]=

問題は,解が急速に変化するとき,局所的確度が優勢な問題であるため硬いソルバを使う点はほとんどないということである.

効率を考えると,局所的確度(安定性ではない)が重要である領域を,メソッドが自動的に検出できたら便利であろう.

線形安定

線形安定理論は,ダールキスト(Dahlquist)のスカラー線形テスト方程式の研究から派生する:

これは初期問題(1)を研究するための簡約化されたモデルである.

安定とは下の

(ここで であり は有理安定関数)を得るために(2)に適用されたメソッドを分析することである.

絶対安定の境界は,領域を考慮することで得られる:

前進オイラー法

前進オイラー法:

を(3)に適用すると,以下が求まる:

陰影の付いた領域は不安定を表し,ここで TemplateBox[{{R, (, z, )}}, Abs]>1である.

In[9]:=
Click for copyable input
Out[9]=

線形安定境界(LSB)は負の実軸との交点と取られることが多い.

前進オイラー法ではである

固有値 では,線形安定の必要条件は,刻み幅が非常にゆるい制約の を満足する必要があるということを意味する.

しかし固有値 のとき,線形安定の必要条件は,刻み幅が非常に厳しい制約である を満足する必要があるということを意味する.

例題

以下の例題は,硬い系を解くのに陽的ルンゲ・クッタ法を使ったときの刻み幅シーケンスに対する硬さの効果を示している.

以下の系は化学反応をモデル化する.
In[10]:=
Click for copyable input
この系は,組込みの硬さ検出を無効にすることで解かれる.
In[11]:=
Click for copyable input
安定領域に達すると,刻み幅シーケンスは振動し始める.
In[12]:=
Click for copyable input
Out[12]=
  • 多数のステップの拒否は,パフォーマンスにマイナスのインパクトを与える.
    • 多数のステップが取られると,計算された解の確度に逆の影響を及ぼす.
    組込みの検出機能は,硬さが発生したときの位置を見付けるのに効果的である.

後退オイラー法

後退オイラー法

を(4)に適用すると

が得られる.このメソッドは左半分の平面には無条件で安定である.

In[14]:=
Click for copyable input
Out[14]=

つまり,安定を保つために刻み幅には制約が必要ないということである.

欠点は,方程式の陰的な系をそれぞれの積分ステップで解かなければならないということである.

タイプの不感性

タイプ不感性ソルバは各ステップで効率よく硬さを認識し反応するため,問題のタイプ(変化)に敏感でない.

このクラスで最もよく構築されたソルバの一つがLSODA [H83],[P83]である.

CVODE等,後の世代のLSODAはもはや固さ検出デバイスを装備しない.それは,後で説明するが,LSODAは優勢固有値を推測するためにノルム限界を使うが,その限界が非常に不正確であり得るためである.

A(α)安定のBDFメソッドの低次数は,高確度の系あるいは優勢固有値に大きい虚部のある系をとくのにLSODAおよびCVODEがあまり適していないことを意味する.線形の陰的スキームの外挿に基づくメソッド等,別のメソッドではそのようなことは問題にならない.

硬さ検出については1980年代,1990年代にスタンドアロンのFORTRANコードを使って多数の研究がされた.

それ以降新しい線形代数技術や効率的なソフトウェアが利用できるようになり,それがWolfram言語ですぐにアクセスできるようになっている.

硬さは一時的な現象であり得るため,硬さのないことを検出することも同様に重要である[S77],[B90].

"StiffnessTest"メソッドオプション

硬さのないソルバから硬い祖父場に切り替えるために使用できるアプローチがいくつかある.

直接の推定

硬さを検出する便利な方法に,問題のヤコビアン の優勢固有値を直接推定するというものがある([S77],[P83],[S83],[S84a],[S84c],[R87],[HW96]を参照).

そのような推定は数値積分の二次的結果として利用できることが多いので,比較的安価である.

がヤコビアンの優勢固有値に対応する固有ベクトルの近似を表し,が十分小さい場合,平均値定理により,主固有値のよい近似は以下のようになる:

リチャードソン(Richardson)の補外法は,ある種の陽的ルンゲ・クッタ(RungeKutta)法で行うように,この公式の量を算出する一連の微調整を提供する.

コストはせいぜい2回の関数評価ですむが,少なくともそのうちの1回は数値積分の二次的結果として利用できるので比較的安価である.

が線形安定領域を表すとする線形安定領域と負の実軸との交点である.

は硬さを検出するためにメソッドの線形安定領域と比較することのできる推定値を与える:

ここで は安全因子である.

解説

メソッドにはオプションがあり,これは与えられた問題に,指定されたAccuracyGoalおよびPrecisionGoalで適用されたメソッドが硬いかどうかを見分けるために使われる.

メソッドオプションの自体は(5)の弱形式を実装する多数のオプションを受け入れる.ここでテストは決まった回数だけ失敗することが許される.

その理由は,問題の中にはある領域だけがほどほどに硬いだけということがあり,明示的な積分法はまだ有効であるためである.

"NonstiffTest"メソッドオプション

メソッドにはオプションがある.これは硬いメソッドから硬くないメソッドへと切り替えるために使われる.

次の設定がオプション に許される.

  • NoneまたはFalse(テストをしない)

硬くないソルバへの切換え

硬いメソッドと独立のアプローチが使われる.

ヤコビアンJ(あるいは近似)があるとき,以下のうちの一つを計算する:

ノルム限界:

スペクトル半径: rho(J)=max TemplateBox[{{lambda, _, i}}, Abs]

優勢固有値

線形代数の手法の多くは,1つの問題を高確度までとくことに焦点を当てている.

硬さ検出では,1桁か2桁までの解を持つ問題の連続が妥当である.

数値的離散化

で,部分区間のいくつかにある行列のシーケンス を考えてみる:

連続行列のスペクトルはステップからステップでゆっくりと変化することが多い.

目的は下のような行列の連続 の優勢固有値(の限界)を推定する方法を探すことである:

  • 硬いソルバの各ステップで,線形代数で実行される操作よりも低いコスト.
  • ソルバのステップごとの特性を考慮する.

NormBound

優勢固有値の限界を得るための簡単で効率的な方法は,ヤコビアンTemplateBox[{J, p}, Norm2] のノルム(ここで通常 あるいは )を使うことである.

この方法には,硬いソルバで実施される作業よりも小さい複雑性 がある.

これはLSODAで使われるアプローチである.

  • 密な行列のノルム限界は多く見積もりすぎ,境界は次元が増えるにつれ悪くなる.
    • ノルム限界は非常に大きい次元の疎行列あるいは帯行列では狭いものになる.

    オプションの設定 TemplateBox[{J, 1}, Norm2]TemplateBox[{J, infty}, Norm2]を計算し,その2つのうちの小さい方を返す.

例題

以下のヤコビ行列は,硬いソルバを使ったファン・デル・ポルの数値解で現れる.
In[15]:=
Click for copyable input
ノルムに基づく限界は,スペクトル半径を大きさの次数より多く見積もる.
In[16]:=
Click for copyable input
Out[16]=

固有値の直接計算

小さい問題(例えば )では,優勢固有値を直接計算した方が効率がよいことがある.

  • エルミート行列はLAPACK関数xgeevを使う.
    • 一般行列はLAPACK関数xsyevrを使う.

    オプションの設定は,Eigenvaluesと同じLAPACKルーチンを使って の優勢固有値を計算する.

    大きい問題では,固有値の直接計算のコストは である.これは硬いソルバでの線形代数のコストと比較するとひどく高いものである.

    この目的のために,多数の反復スキームが実装されている.これらは小さい部分空間での優勢固有空間を近似し,小さい問題では密な固有値メソッドを使って効率的に動作する.

ベキ乗法

Shampineはヤコビアンの優勢固有値を推定するためにベキ乗法を使うことを提唱した[S91].

ベキ乗法はあまり高く評価されている方法ではないが,Google社のページランキングに使用されているため,非常に関心が高まっている.

ベキ乗法は次の場合に使うことができる.

  • n 個の線形独立固有ベクトルを持つ(対角化可能)
  • 固有ベクトルがとして大きさで順序付けることができる
  •  lambda_1A の優勢固有値である.

解説

出発ベクトルを与えると は以下を計算する.

優勢固有値の近似を計算するために,レイリー商が使われる:

実際には,近似固有ベクトルは各ステップでスケールされる:

プロパティ

ベキ乗法は以下の割合で線形に収束する:

これは遅い.

特に,このメソッドは優勢固有値の複素共役ペアを持つ行列に適用されると収束しない.

一般化

ベキ乗法は,同一モジュール固有値の問題を克服するために適用することができる(例:NAPACK).

しかし,個の修正は一般にクラスタ化された固有値の遅い収束率という問題は回避しない.

ベキ乗法を一般化するための主なアプローチが2つある:

  • 小さい次元から中程度の次元では部分空間反復
    • 大きい時限ではアーノルディ反復

    これらのメソッドはかなり異なる動作をするが,共通で最適化できる多数の中心的コンポーネントがある.

    部分空間およびクリロフ反復は O(n2 m) 回の操作を必要とする.

    これらは 行列を 行列(ここで )に投影する.

    小さい行列は,優勢固有空間を表し,近似は密な固有値ルーチンを使用する.

部分空間反復

部分空間(同時)反復は,各ステップで m 個のベクトルに作用することによりベキ乗法の考えを一般化する.

ベクトル の正規直交集合から始める.ここで通常 である.

行列 との積から:

すべてのベクトルが の同じ優勢固有ベクトル の倍数へと収束するのを妨げるため,正規直交化される:

正規直交化のステップは,行列の積と比較すると高価である.

レイリー・リッツ(Rayleigh-Ritz)の投影法

入力:行列 A とベクトルの正規直交集合 V

  • レイリー商 を計算する
    • シューア分解 を計算する

    行列 S は小さい次元 を持つ.

    シューア分解は S in R^(m x m) のとき,擬似上三角行列 を使って実数計算で求めることができる.

収束

部分空間(同時)反復は,各ステップで 個のベクトルに作用することによりベキ乗法の考えを一般化する.

SRRITは割合

で線形に収束する.特に,優勢固有値の割合は

となる.したがって,優勢固有値にのみ関心があるとしても例えば 以上を取るとよいことがある.

誤差制御

優勢固有値の連続した近似に対する相対的な誤差検定は

である.これは収束が遅いときは満足されないので十分ではない.

あるいはならば, 番目の列は一意に決められない.

SRRITで使われる残差検定は

であり,ここで

番目の列, 番目の列である.

これは同一モジュラの固有値にも使えるので便利である.

次数付きのシューア分解を使用するため,上三角行列 の最初の列の位置が検定される.

実装

部分空間反復の実装にはいくつかある.

  • チェビシェフの加速を伴う部分空間反復 [S84b], [DS93]
    • シューア,レイリー・リッツ反復 ([BS97] と [SLEPc05])

    での使用のための実装は以下に基づく.

    • シューア,レイリー・リッツ反復 [BS97]

    「SRRITの魅力的な機能はモノトニックな一貫性を表示することである.つまり収束の誤差許容範囲は計算された残差の大きさのように減少する」[LS96].

    SRRITは,最大のモジュロの固有値が左上の項目に現れる,次数付きのシューア分解を利用する.

    再正規直交化の修正グラム・シュミット(Gram-Schmidt)は を形成するために使われる.これはハウスホルダー(householder)変換より速い.

    次の積分ステップ における反復を開始するために,積分時間 における近似優勢部分空間 が使われる:

KrylovIteration

が,指定された部分空間 の直行基底をなすような 行列 があるとする:

レイリー・リッツ法は, の計算と関連する固有問題 を解くことからなる.

もとの問題 の近似固有ペアはそれぞれリッツ値およびリッツベクトルと呼ばれる および を満足する.

この方法は部分空間 の不変量部分空間を近似するときに最も効果がある.

この方法は, が行列 および与えられた初期ベクトル に関連付けられたクリロフ部分空間に等しい場合に効果的である:

解説

アーノルディ法はクリロフ部分空間の直交基底を計算し, のときの投影された 行列 を生成するクリロフベースの投影アルゴリズムである.

入力:行列 ,ステップ数 ,ノルム1の初期ベクトル

出力:であり

アーノルディ法の場合, は簡約されないヘッセンベルグ形式(追加の非零部分対角要素を含む上三角)を持つ.

正規直交化は通常グラム・シュミット法で実行される.

このアルゴリズムで計算された数量は以下を満足する:

残差 は不変量部分空間への近さを示し,関連付けられたノルム は計算されたリッツペアの正確さを示す:

再出発

リッツペアは,初期ベクトル が所望の固有値の方向に富んでいれば迅速に収束する.

そうでない場合,作業およびメモリが過剰になるのを避けるために再出発が必要になる.

再出発には,特に以下のようないくつかのストラテジーがある:

  • 明示的再出発新しい出発ベクトルはリッツベクトルの部分集合の線形結合である.
    • 暗示的再出発新しい出発ベクトルはアーノルディ法と暗示的にシフトされたQR法を組み合せたもので形成される.

    明示的再出発は比較的実装が簡単であるが,暗示的再出発は大きい問題に必要な固有情報を保存しているため,より効率的である.しかしながら,暗示的再出発は数値的に安定した方法では実装が難しい.

    実装がずっと簡単であるが暗示的再出発と同じくらい効果的である別の方法は,クリロフ・シューア法である[S01].

実装

多数のソフトウェア実装が利用できる:

    の実装はクリロフ・シューア反復に基づいている.

自動ストラテジー

設定では,下のようなメソッドの融合が使われる.

  • では,固有値直接計算が使われる.のいずれが使われるかは,どちらのメソッドがアクティブかに依存する.
  • では,デフォルトの基底サイズが の部分空間反復が使われる.この方法が成功すると,次の積分ステップでメソッドを開始するのに結果となった規定が使われる.
  • 部分空間反復が最大反復回数後に収束に失敗すると,デフォルト基底サイズ のクリロフ方を開始するために,優勢ベクトルが使われる.後続する積分ステップでは,前のステップからの結果のベクトルで開始したクリロフ法が使われる.
  • クリロフ反復が最大反復回数の後に収束に失敗すると,現在のステップにノルム限界が使われる.次の積分ステップでは引き続きクリロフ反復を使おうとする.
  • これらは高価ではないので,部分空間あるいはクリロフ反復が使われ,そのうちの絶対値の小さいほうが使われると,ノルム限界は常に計算される.

ステップの却下

評価時間のキャッシュは,優勢固有値推定が却下されたステップに対して再計算されないことを確実にする.

硬さ検出は以下の理由で,却下されたステップに対しても実行される.

  • ステップの却下は,安定領域付近で作業しているときに硬くないソルバに対してよく起る
  • ステップの却下は,早い過渡値を解決するときに硬いソルバに対してよく起る

反復メソッドオプション

の反復メソッドには,変更することのできるオプションがある:

In[20]:=
Click for copyable input
Out[20]=
In[21]:=
Click for copyable input
Out[21]=

デフォルトの誤差許容範囲はただ1つの正しい数値を目標としているが,特に連続するステップで数回の反復が成功した後,実質的にもっと多くの正確な値を得ることがある.

反復の回数を制限するデフォルトの値は以下のようなものである:

  • 部分空間反復
    • クリロフ反復

    これらの値が大きすぎると,収束の失敗は非常にコストの高いものとなる.

    難しい問題では,ステップ全体で収束作業を分けた方がよい.このメソッドは前のステップからの基底ベクトルを効率よく調整するので,後続のステップで収束が成功する可能性もある.

レイテンシと切換え

メソッドが継続的に硬いメソッドと硬くないメソッドとを切り換えようとするサイクルを避けるために,ある形でレイテンシを含ませることが重要である.

およびのオプションはこの目的で使用される.

デフォルト設定では,切換えが非常に反応が早い.これは1ステップの積分メソッドに適している.

  • は硬くないメソッドのステップの最後に実行される.オプションのどちらかの値に達したら,ステップの却下が起り,ステップは硬いメソッドで再計算される.
  • は先手を取るものである.これは前のステップのヤコビアン行列を使って,硬いソルバでステップを取る前に実行される.

例題

ヴァンデルポール

例題の系を選ぶ.
In[17]:=
Click for copyable input

StiffnessTest

この系は,与えれたメソッドとのデフォルトオプション設定での積分が成功する.
In[18]:=
Click for copyable input
Out[18]=
積分が長くなると放棄され,硬さ検出条件が満足されないとメッセージが出力される.
In[19]:=
Click for copyable input
Out[19]=
単位安全因子を使い,硬さの失敗が1回だけ許されることを指定すると,厳密な検定が効果的に与えられる.この指定ではネストされたメソッドオプションのシンタックスが使われている.
In[20]:=
Click for copyable input
Out[20]=

NonstiffTest

小さい系の場合は,直接の固有値計算が使われる.

この例題は,全体的な硬さ切換えフレームワークが期待通りに動いているというより検定の例である.

硬いメソッドと硬くないメソッドとの切換えおよび取られるステップを監視する関数を設定する.硬いソルバと硬くないソルバのデータは,"Sow"に対して異なるタグを使うことで別々のリストに入れられる.
In[21]:=
Click for copyable input
系を解き,メソッド切換えのデータを収集する.
In[23]:=
Click for copyable input
明示的ソルバ(青)と暗示的ソルバ(赤)を使って,取られる刻み幅をプロットする.
In[25]:=
Click for copyable input
Out[25]=
取られた硬くないステップおよび硬いステップの数(却下されたステップを含む)を計算する.
In[26]:=
Click for copyable input
Out[26]=

CUSP

神経インパルスメカニズムに対するcusp catastrophyモデル[Z72]:

ヴァンデルポールと組み合せるとCUSP系になる[HW96]:

ここで

であり,および である.

線の方法を使った拡散項の離散化が,次元の常微分方程式系を得るために使われる.

ヴァンデルポール系とは異なり,問題の大きさのため,固有値推定には反復メソッドが使われる.

刻み幅と次数選択

解く問題を選ぶ.
In[27]:=
Click for copyable input
使用されるメソッドのタイプと刻み幅を監視する関数を設定する.さらにメソッドの次数がツールチップとして含まれる.
In[28]:=
Click for copyable input
積分が進むにつれて,メソッドの次数のデータを収集する.
In[30]:=
Click for copyable input
明示的ソルバ(青)と暗示的ソルバ(赤)を使って,取られる刻み幅をプロットする.各ステップでのメソッドの次数は,ツールチップで示される.
In[32]:=
Click for copyable input
Out[32]=
取られる硬くないステップと硬いステップの合計数(却下されたステップを含む)を計算する.
In[39]:=
Click for copyable input
Out[39]=

ヤコビアンの例題

最初の2,3のヤコビアン行列を集める関数を定義する.
In[33]:=
Click for copyable input
In[35]:=
Click for copyable input

硬いメソッドへの切換えは0.00113425付近で起きており,硬くないことの最初の検定は次のステップ で起きている.

ヤコビアン のグラフィックスによる表示である.
In[37]:=
Click for copyable input
Out[37]=
計算する関数を定義し, , , の最初のいくつかの固有値とノルム限界を計算する関数を定義する.
In[38]:=
Click for copyable input
In[39]:=
Click for copyable input
Out[39]=

この例題では,ノルム限界は非常に急である.

KortewegdeVries

KortewegdeVriesの偏微分方程式は浅い水面の波動の数学的モデルである:

次の境界条件

を考え,区間 上で解く.

線の方法を使った離散化は192の常微分方程式系を形成するために使われる.

刻み幅

解く問題を選ぶ.
In[40]:=
Click for copyable input
LSODAで使用されるBackward Differentiation Formulaメソッドでは,この問題を解く際に困難が生じる.
In[41]:=
Click for copyable input
Out[41]=
このプロットで刻み幅が急激に減少しているのが分かる.
In[42]:=
Click for copyable input
Out[42]=
これとは対照的に,はほとんど積分ステップを必要としない線形後退オイラー法を使うよう,直ちに切り替える.
In[43]:=
Click for copyable input
Out[43]=
In[44]:=
Click for copyable input
Out[44]=

積分のはじめに硬いソルバが選ばれたら,外挿法は決して硬くないソルバに切り換えない.

従って,これは硬くないことの検出で最悪の場合の例だといえる.

それでも部分空間反復を使用するコストは積分全体の時間のほんの数パーセントに過ぎない.

硬くないメソッドへの切換えを無効にして時間を計算する.
In[45]:=
Click for copyable input
Out[45]=

ヤコビアンの例

前に定義した監視関数を使って最初のいくつかのヤコビアン行列のデータを集める.
In[46]:=
Click for copyable input
初期ヤコビアン のグラフによる表示である.
In[48]:=
Click for copyable input
Out[48]=
最初のいくつかの , , の固有値およびノルム境界を計算し表示する.
In[49]:=
Click for copyable input
Out[49]=

ノルム境界はわずかに過大評価であるが,より重要なのは,実部および虚部の相対的な大きさにてついて言及していないという点である.

オプションの要約

StiffnessTest

オプション名
デフォルト値
"MaxRepetitions"{3,5}硬さ検定(6)の失敗が許される連続および合計回数を指定する
"SafetyFactor"硬さ検定(7)の右辺と左辺に使う安全係数を指定する

メソッドオプションのオプション

NonstiffTest

オプション名
デフォルト値
"MaxRepetitions"{2,}硬さ検定(8)の失敗が許される連続および合計回数を指定する
"SafetyFactor"硬さ検定(9)の右辺と左辺に使う安全係数を指定する

メソッドオプションのオプション