10. SVDアシスト ......................................................................................................................199
10.1 基本概念....................................................................................................................199
10.1.1 スーパーパラメータ...............................................................................................199
10.1.2 実装概要.................................................................................................................200
10.1.3 スーパーパラメータの数は?............................................................................202
10.1.4 利点と欠点.............................................................................................................202
10.2 SVDAPREP ..............................................................................................................203
10.2.1 PEST入力データセットの準備........................................................................203
10.2.2 SVDAPREPの実行........................................................................................204
10.2.3 SVDAPREPの機能........................................................................................206
10.2.4 新しいモデルバッチファイル...........................................................................207
10.3 PARCALCおよびPICALCユーティリティ.....................................................................208
10.3.1 PARCALC ........................................................................................................208
10.3.2 PICALC.............................................................................................................208
10.4 スーパーパラメータによるPEST制御ファイル..................................................................209 10.4.1 全般................................................................................................................209
10.4.2 スーパーパラメータ定義に関する重要な注意点................................................209
10.4.3 「SVDアシスト」セクション.............................................................................210
10.5 SVDアシストによる逆解析............................................................................................217
10.5.1 PESTの処理.......................................................................................................217
10.5.2 最適パラメータ推定..........................................................................................218
10.5.3 並列SVDアシスト逆解析..................................................................................218
10.6 最終的な留意事項.....................................................................................................219
10.6.1 JCO2JCOおよびその他のユーティリティプログラム............................................219
10.6.2 効率性................................................................................................................219
10.6.3 ヌル空間モンテカルロ....................................................................................220
- SVDアシスト
10.1 概念
10.1.1 スーパーパラメータ
PESTにおける「特異値分解アシスト(SVD-Assist)」を活用することで、高次元パラメータの逆解析において実行時間の大幅削減が可能となり、処理困難な問題も数値的に扱いやすくなります。
この手法では、特定の「スーパーパラメータ」と呼ばれる変数群を定義し、これらを直接推定します。
これらのスーパーパラメータは、全モデルパラメータを含む重み付きヤコビ行列の特異値分解に基づいて定義されます。
数学的には、スーパーパラメータとは、元のパラメータ空間を補正解空間の直交軸に射影したスカラー値です。
これらのスカラー値を推定し、それぞれに対応する軸方向の単位ベクトルに掛け合わせることで逆問題の解を得ます。その後、これらのベクトルを足し合わせて推定されたパラメータ集合を作成します。
理想的には、観測データとの良好な適合を実現するために必要なスーパーパラメータ数は、補正解空間の次元数に等しくなります。これは、モデルで実際に使用されるパラメータ数よりも小さい場合があります。しかし、これは補正データセットの情報量に基づいて最適に定義されます。したがって、これらはその情報を最も効率的に表現するものです。さらに、必要な最小限の数です。詳細についてはDoherty(2015)の6.2.6節を参照してください。
スーパーパラメータが定義されれば、PESTは各イテレーションでスーパーパラメータ数だけのモデル実行を行うだけで済みます。なぜなら、PESTはこれらのスーパーパラメータに対する有限差分導関数を計算するためです。
このため、SVD-assistにより高度なパラメータ化された逆解析の計算負荷が大幅に軽減できます。一方、バックグラウンドではPESTがスーパーパラメータの値をローカルパラメータに変換します。
スーパーパラメータについてユーザーが知る必要があるのは、それらがPESTの労働時間削減のための手段であり、モデル実行の効率性を高めるために使用されることだけです。実際の値は重要ではありません。
SVDアシストを逆問題解決手段として使用しても、Tikhonov正則化の使用は除外されません。実際、前章やDoherty(2015)で述べた理由により、問題が不安定な場合は常にTikhonov正則化の使用が推奨されます。これにより、専門知識に従った解が得られ、誤差分散の最小化が可能になります。そして、この解は後処理不確実性解析の最良の基盤となります。幸いにも、SVD-assistを使用しても、正則化関係をスーパーパラメータで表現する必要はありません。代わりに、元のパラメータで表現できます。PESTが自動的に変換を行います。(ただし、Tikhonov正則化を使わずにSVD-assistを適用することも可能です。これは必須ではありません。)
10.1.2 実装の概要
次の説明では、PESTによる正則化逆解析を使用していることを前提とします。したがって、基本パラメータPESTコントロールファイルのpesmodeは「estimate」または「regularization」に設定されている必要があります。ただし、SVD-assistはPEST予測解析やPareto操作の解決手段としても使用できますが、これはまれです。したがって、pesmodeは「estimate」「regularization」「prediction」「pareto」のいずれかに設定できます。
逆解析はSVDアシストなしでも可能です。
SVD-assistはあくまで解決戦略です。
使用するかどうかは数値的な利便性の問題です。
したがって、前章で提案された数値安定な正則化逆解析の実装方法に従うべきです。
• 可能な限り多くのパラメータを逆解析に含めることで、シミュレーションシステムの非一様性を反映させます。
これにより、逆解析の柔軟性が保たれ、非一様性が補正データの情報内容に従って出現します。また、補正後の不確実性解析がパラメータの非一意性が予測の不確実性に与える影響を評価できるようになります。
• 逆解析プロセスにTikhonov正則化を追加することで、すべてのパラメータにバックアップ位置を提供し、データが保証されている場合にそれらから逸脱する方法を実現します。
目標測定目的関数を適切な値に設定します。
これは初期の低速PEST実行、または過剰適合を抑制するためのレベルで行うことができます。
• PESTコントロールファイルに特異値分解またはLSQRセクションを含めます。これにより、逆解析プロセスの数値的安定性が保証されます。
• rlamda1を10、RLAMFACを-3に設定します。マーカートのλの高速移動は数値不安定性への最終防衛線であり、モデル線形性への適応も可能です。
• モデルの数値性能に懸念がある場合は、第4章で述べられているスプリットスロープ分析やモデル実行の許容機構などの他の支援ツールも考慮すべきです。
• 異なる観測グループ間の重み相関が正しく設定されていることを確認します。各グループが初期測定目的関数に同じ貢献をするようにします。(本手順書の第2部で説明されているPWTADJ1ユーティリティが役立ちます。)PEST入力データセットが準備できたら、PESTコントロールファイルの「control data」セクションでNOPTMAX変数を
-2に設定します。その後PESTを実行します。すると、PESTはローカルモデルパラメータに基づいてヤコビ行列を計算し、それをファイルに保存します。ファイル名はcase.jcoで、caseはPESTコントロールファイルのベース名です。このファイルの保存後、PESTは即座に実行を停止します。
次に、使用するスーパーパラメータ数を決定します(後述)。その後、SVDAPREPユーティリティを実行して、SVDアシスト逆解析のための新しいPESTコントロールファイルと補助ファイルを構築します。便宜上、このファイルをcase-svda.pstと呼びます。このファイルの詳細は必要ありません。
しかし、興味のある読者にはこのファイルの仕様が後で示されます。
その後、以下のコマンドで新しいPESTコントロールファイルをベースにしてPESTを実行します。
pest case-svda
任意のPESTバージョンが使用でき、並列PESTやBEOPESTも可能です。
SVDアシストの最初のイテレーションは非常に高速に行われます。なぜなら、PESTはcase.jcoに保存された基本パラメータ導関数からスーパーパラメータ導関数を計算するからです。
しかし、以降のイテレーションでは、前述の標準有限差分法によってスーパーパラメータ導関数を計算します。
SVDアシストによるモデル実行時、PESTは通常のPEST実行と同じ情報を画面に出力します。
特に、測定と正則化の目的関数、および各観測グループの目的関数への寄与が記録され、逆解析の進行状況を監視できます。
実行ログファイルにも同じ情報が記録されます。
同時に、case-svda.reiという名前の一時残差ファイルなども記録され、特定の残差の減少状況を監視できます。
SVDアシスト逆解析の場合、第5節で説明されている再開スイッチを使用して、通常通りPESTを停止・再開できます。
逆解析のどの段階でも、最適なパラメータ値はcase-svda.parファイルに保存されます。
ただし、これはスーパーパラメータの最適値を記録しているため、ユーザーにとってあまり有用ではありません。スーパーパラメータは「par1」「par2」などの名前で記録されています。
基本パラメータの最適値はcase.bpaファイルに保存されます(「Bpa」は「best parameter values」の略)。
このファイル名は基本PESTコントロールファイルのベース名ではなく、スーパーパラメータPESTコントロールファイルのベース名です。
SVDアシスト逆解析では、通常のPEST逆解析と同じ終了基準が使用されます。
しかし、既に十分な作業が完了していると感じられた場合は、PESTが「終了」を宣言するのを待つ必要はありません。
その場合は、通常通りPESTを手動で停止できます。
PESTが自ら終了したとしても、SVDアシスト逆解析では最適パラメータを使った最終モデル実行は行われません。
なぜなら、逆解析の進行に伴ってスーパーパラメータと基本パラメータの関係が変化する可能性があるためです。
そのため、初期イテレーションでの最適パラメータをそのまま使用すると、現在のスーパーパラメータ値から基本パラメータ値への変換が失われる可能性があります。
幸いにも、SVDアシスト逆解析から得られた最適パラメータを用いて手動でモデル実行することは簡単です。PARREPユーティリティ(本手順書第2部参照)を使用し、以下のコマンドで新しい基本PESTコントロールファイル(ここではcase-soln.pst)を作成します。その初期値は最適パラメータです。
parrep case.bpa case.pst case-soln.pst
その後、case-soln.pstファイルでNOPTMAXを0に設定します。PESTを実行します。これにより、最適パラメータを使ってモデルを一度だけ実行し、実行を停止します。あるいは、SUBREG1ユーティリティを使用してcase-solnから正則化を削除し、新しい雅可ビ行列を計算するためにNOPTMAXを-2に設定することもできます。PWTADJ2ユーティリティを使用して、新しいPESTコントロールファイルの重みを調整し、モデルと測定の不適合のランダム特性を反映させることも可能です。その後、本手順書第2部で説明されているGENLINPREDなどのユーティリティを使用して、補正後の線形パラメータと予測の不確実性解析を実行できます。
10.1.3 スーパーパラメータの数は?
使用するスーパーパラメータの数は?可能な限り多くの計算リソースを使用してください。
本手順書第2部で述べられているSUPCALCユーティリティは、特定の逆問題の解空間の最適次元を計算しようとします。これは基本パラメータで計算されたヤコビ行列が必要です。
しかし、SVD-assist PESTコントロールファイルに含まれるスーパーパラメータ数をsupcalcで算出された解空間次元に限定する必要はありません。
以下で述べるように、スーパーパラメータの定義は線形仮定に基づいており、多くのモデルはこの仮定を破る可能性があります。解空間とヌル空間の次元と構成はパラメータ値の変化に応じて変化します。推定されるスーパーパラメータが多いほど、この問題は小さくなります。最後に、使用するスーパーパラメータ数の唯一の考慮点は、利用可能な計算リソースです。もし十分な計算リソースがある場合、SVD-assistを数値手段として使用する必要はありません。
定義上、スーパーパラメータ数が解空間の次元を超える場合、スーパーパラメータ推定問題は不適定になります。通常これは問題ではありません。Tikhonov正則化が使用されている場合、逆問題は最小誤差分散解に向かいます。さらに、基本パラメータPESTコントロールファイルでSVDまたはLSQRが指定されている場合、SVDAPREPユーティリティ(以下参照)は、スーパーパラメータPESTコントロールファイルの解決メカニズムとしてそれらを指定します。したがって、スーパーパラメータ数に関係なく、SVDアシスト逆解析プロセスの数値的安定性は保証されます。
場合によっては、可用な計算リソースがSUPCALCが推奨する最小スーパーパラメータ数を満たすことができないこともあります。その場合は、可能な限り多くのスーパーパラメータを使用します。Doherty(2015)によれば、スーパーパラメータ定義に含まれるパラメータ簡略化は最適であり、これは最小限の「パラメータコンテナ」で補正データセットの情報を最大限に受け入れ可能にするからです。最終的に、環境モデリングと環境モデル補正は妥協の問題です。スーパーパラメータは非常に良い妥協点を持っています。
10.1.4 利点と欠点
SVDアシストの強みは、モデル実行時間の大幅な削減です。逆解析中の特徴パラメータ数に対して、補正解空間の次元が小さいほど、その節約効果は大きくなります。
SVDアシストの大きな欠点は、スーパーパラメータ定義が線形仮定に基づいていること、つまりモデルがその仮定に従わない可能性があることです。前述のように、解空間に必要な以上のスーパーパラメータを使用することで、非線形性に対する防御が可能です。
別の防御戦略として、SVDアシスト逆解析プロセスの複数イテレーション後にスーパーパラメータ定義プロセスを繰り返す方法があります。例えば、PESTが第三次イテレーションで測定目的関数の減少にほとんど進展がない場合を考えます。(まず、スーパーパラメータの移動に相対的または係数制限がかけられていないか、基本パラメータの多くが境界に達していないかを確認してください。)その後、PESTの実行を停止し、PARREPユーティリティを使用して、停止したSVDアシスト逆解析から得られた最適パラメータを持つ新しい基本パラメータPESTコントロールファイルを作成できます。新しいファイルでNOPTMAXを-2に設定してPESTを実行し、新しい基本パラメータヤコビ行列を計算します。その後、SVDAPREPを使用して新しいスーパーパラメータPESTコントロールファイルを構築します。このファイルを使用してPESTを実行し、SVDアシスト逆解析プロセスを再開します。
モデルの非線形性の影響を緩和するためのさらなる対策も取れます。通常通り、基本PESTコントロールファイルの初期パラメータ値は専門家の知識に基づく優先値であるべきです。そして、Tikhonov正則化の逆問題の解は、補正データセットに適合しつつ、優先パラメータセットから最も遠くないパラメータセットを求めます。
場合によっては、特に地下水や水庫モデリングでは、PESTコントロールファイル内のすべてのパラメータの優先値を決定することが難しいことがあります。しかし、SVDアシスト逆解析を良いスタートにする方法もあります。
例えば、マルチレイヤーモデルのパラメータを推定しており、各レイヤーに多くのパラメータ(数百個)が割り当てられているとします。最初のパラメータ推定実行では、各レイヤーのすべてのパラメータを単一のパラメータに束縛して、レイヤー全体を表すようにします。これにより、逆問題の次元がモデル内のレイヤー数に減少します。各レイヤーの初期値が同じであれば、この戦略はレイヤー属性の一様性を仮定しています。このようにして定義された少数のレイヤー固有パラメータが推定されると、パラメータを相互に分離し、Tikhonov正則化を導入し(レイヤー一様性パラメータ制約をカプセル化)、NOPTMAXを-2に設定してPESTを実行してすべてのモデルパラメータのヤコビ行列を計算できます。その後、上記のようにSVDアシスト逆解析を実行できます。
SVDアシスト逆解析を実行する際、基本パラメータ境界の強制は妥協が必要な領域です。PESTがスーパーパラメータを推定する際には、基本パラメータPESTコントロールファイルで設定された境界が適用されます。スーパーパラメータの更新ベクトルが基本パラメータの境界に達する場合、PESTは更新ベクトルを短くしてそのような事態を回避します。(これは測定目的関数のイテレーション内での改善を制限する可能性があります。)その後、PESTは基本パラメータを逆解析から削除し、境界に永続的に固定します。これには二つの欠点があります。まず、基本パラメータの数値が非常に高い場合、スーパーパラメータの再定義は数値的に計算コストが高く(つまり遅い)なります。第二の欠点は、境界付きパラメータがさらに調整される機会を失うことです。しかし、スーパーパラメータとその封入する基本パラメータを開放状態に保つには、この戦略以外に選択肢はありません。
したがって、可能であれば、基本パラメータPESTコントロールファイルの境界は広めに設定すべきです。
10.2 SVDAPREP
10.2.1 PEST入力データセットの準備
前述のように、SVDアシスト逆解析の最初のステップは、通常の逆解析と同じく、PEST入力データセットの構築です。これは、コントロールファイルを含むものであり、必要に応じてTikhonov正則化を含むかもしれません。必要な場合は、固定パラメータや束縛パラメータを含めることも可能です。単一のパラメータは対数変換または変換なしでも構いません。準備ができたら、PESTCHEKユーティリティを使用してPEST入力データセットを確認します。このPESTコントロールファイルは「基本パラメータPESTコントロールファイル」と呼ばれます。
その後、NOPTMAXを-2に設定してPESTを実行し、ヤコビ行列を計算します。
ヤコビ行列が計算された後、SVDAPREPユーティリティを実行して、スーパーパラメータに基づいてPESTデータセットを構築します。(もしSVD-assistを逆問題解決手段として使用しないことに変わりがあった場合、ヤコビ行列の計算に費やされたモデル実行時間は無駄ではありません。PESTコントロールファイルのNOPTMAXをより大きな数値に設定し、「/i」スイッチでPESTを起動します。前述のように、PESTは最初のイテレーションに使用するヤコビ行列の名前を尋ねます。)基本パラメータPESTコントロールファイルは他のPESTコントロールファイルと同様に使用できますが、SVDAPREPは特定の規則に従う必要があります。(そうでない場合は、実行を停止しエラーメッセージを表示します。)モデル実行コマンドは「.bat」で終わる必要があります。Windowsでは、これはバッチファイルを意味します。Unixではこの命名規則はスクリプトファイルを意味しませんが、SVDAPREPはそれに従う必要があります。
実際には、PEST実行の「モデル」はバッチファイルまたはスクリプトファイルである必要があります。なぜなら、SVDAPREPはそのファイルを複製し、その上部にいくつかの行を追加する必要があるからです。これらの行はPARCALCおよびPICALCユーティリティを実行するコマンドであり、スーパーパラメータを基本パラメータに変換し、基本パラメータ間の関係を示す事前情報方程式を評価します。これらのユーティリティはPATH環境変数で指定されたディレクトリ(つまりフォルダ)に存在する必要があり、OSがそれらを実行する際にどこにあるかを認識できるようにする必要があります。あるいは、PARCALCおよびPICALC実行ファイルをPEST作業ディレクトリにコピーする(Parallel PESTまたはBEOPESTを使用している場合は各エージェントの作業ディレクトリにコピーする)必要があります。
10.2.2 SVDAPREPの実行
SVDAPREPを画面プロンプトで実行します。多くのプロンプトには、デフォルト値を受け入れるキー入力で応答できます。
実行開始時に、SVDAPREPは基本パラメータPESTコントロールファイル名を尋ねます。
Enter name of existing PEST control file:
(「.pst」拡張子を省略した場合、SVDAPREPは自動的に追加します。)その後、SVDAPREPはこのファイルを読み込みます。また、case.jcoファイル(caseは基本パラメータPESTコントロールファイルのベース名)を検索します。これは上記の方法でPESTを実行して計算された基本パラメータヤコビ行列を含みます。次の質問は:
Use pre-defined super-parameter file? [y/n] ( if "no"):
この質問に「y」と答えた場合、SVDAPREPはスーパーパラメータファイル名を尋ねます。新しいPESTコントロールファイルにSVDA_EXTSUPER制御変数を1に設定し、この新しいPESTコントロールファイルにスーパーパラメータファイル名を記録します。このセクションの後半でスーパーパラメータPESTコントロールファイルの他の特徴について説明します。
外部でスーパーパラメータを定義しない場合(通常はその通り)、SVDAPREPの次のプロンプトは:
For computation of super parameters:
SVD on Q^(1/2)X - enter 1
if SVD on XtQX - enter 2
if LSQR without orthogonalisation - enter 3
if LSQR with orthogonalisation - enter 4
これらのプロンプトでは、Xは基本パラメータに関連するヤコビ行列(通常はJと表記)であり、正則化は除外されています。Qは観測重み行列です。LSQRはQ^(1/2)Jに適用されます。
オプション1が最も良いですが、基本パラメータが2500を超える場合はオプション3または4を選択します。パラメータ数が増加すると特異値分解は非常に遅くなるためです。注意:SVDAPREP自体はスーパーパラメータを計算しません。代わりに、PESTがどのように計算するかをスーパーパラメータPESTコントロールファイルで適切な制御変数を設定することで伝えます(後述)。PESTはSVDアシスト逆解析の開始時およびその後、基本パラメータが境界に達した場合にスーパーパラメータを計算します。
次にSVDAPREPは:
Enter number of super parameters to estimate:
適切な数値を入力します。前述の議論を参照してください。
SVDAPREPは次に:
Enter name of new super pest control file:
新しいPESTコントロールファイル名を入力します。(「.pst」拡張子を省略した場合、SVDAPREPが自動的に追加します。)この新しいPESTコントロールファイルはSVDアシストパラメータ推定に使用されます。次に、SVDAPREPは:
Enter offset for super parameters ( if 10):
SVDアシスト逆解析中、スーパーパラメータの初期値はゼロ(基本パラメータの摂動がゼロを示す)です。しかし、ゼロ値のパラメータはPESTに問題を引き起こす可能性があり、特にパラメータ変更制限を強制する場合です。したがって、これらのパラメータにオフセットを与えることで、値をゼロから離すのが望ましいです。ほとんどの場合、SVDアシストパラメータ推定の初期値は10が適切です。SVDAPREPがこれを受諾するようにします。
SVDAPREPはスーパーパラメータを「相対制限」に指定します。本手順書第4章で述べたように、この相対制限の実際の値はPEST変数RELPARMAXとして提供される必要があります。この変数は通常よりも大幅に低く設定する必要がありますが、あまり低すぎるとパラメータ更新が小さくなり使い物にならなくなります。ほとんどの場合、0.1が適切ですが、PESTの収束が遅い場合はこれを上げ、パラメータの振動や境界に早く達する場合は下げてください。SVDAPREPは、これから書き込まれる制御ファイルに記録するために、ユーザーにRELPARMAXの適切な値を要求します。