非線形逆問題における局所解の探索手法

土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
増本清（島根大）
非線形逆問題における局所解の探索手法
１非線形逆問題の定式化
ここでは，微分方程式などで記述される様々な現象に対して，より一般的な逆問題を考
える．物性値などのパラメータ k ( x, y, z, t ) と状態量 h( x, y, z, t ) とを含む，微分方程式
R(h, k , x, y, z, t )  0 と付随する条件（初期・境界条件）を時空間で離散化することにより，
   
(1)
R( h , k )  0
なる連立方程式が得られる．
ただし，
 
 
 

R  ( R1 (h , k ), R2 (h , k ),  , RI (h , k ))T

h  (h1 , h2 ,  , hI )T ：離散化した各節点における状態量（I は時空間の総節点数）

k  (k1 , k2 ,  , k M )T ：システムを構成する物性値などのパラメータ（逆問題においては未知
パラメータになる．M は未知パラメータの次元）
とする．また，上付き添え字 T は行列の転置を表す．


k が与えられたとき，(1)式を満たす h を求める問題が順問題に相当する．順問題が解け



るとき，それを解いて得られた h は k の関数とみなせるので，このときの h を，
  
h  H (k )
(2)
と書くことにする．
 
 
 
一方，(1)式を満たす h ，k の内，目的関数 e(h , k ) を最小とする h ，k を求める問題（以下，
問題 A とよぶ）が逆問題に相当する．
 
   
e(h , k )  min . （制約条件： R(h , k )  0 ）
［問題 A］
これは(1)式を制約条件とする非線形最適化問題とみなすことができる．最小化すべき目的
 

関数 e(h , k ) も k のみの関数として，
 

 
e(h , k )  e( H (k ), k )  E (k )
(3)
と書ける．このとき，逆問題は，

E (k )  min .
［問題 B］


となる k を求める制約条件無しの最小化問題（問題 B）に帰着できる．ただし， k 自体に不
等式制約条件などがある場合は，制約条件付きの非線形最小化問題となるが，ここでは簡

単のため制約条件無しで説明する．問題 B が非線形問題のとき，与えられた k に対して，
土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
増本清（島根大）

目的関数 E (k ) の勾配


 T
  E (k ) E (k )

E
(
k
)

Ek (k )  
,
,  ,


k

k

k
1
2
M 

(4)
が求められれば，後述する勾配法（最急降下法，準ニュートン法など）を用いることによ
り問題 B の解を得ることができる．そのとき，解は，
 
Ek (k )  0
(5)
を満たす．問題 A と問題 B は実質的には同じである．

 
また，目的関数が計算値 H i (k ) （ H (k ) の i 番目の要素）と観測値 hobs,i の重み付き残差二
乗和で表される場合，特に非線形最小二乗問題（問題 C）と呼ばれ，この形に特化した非線
形最小二乗法（Gauss-Newton 法，Marquardt 法など）を利用できる．

I


E (k )   wi H i (k )  hobs,i
i 1
  min .
2
［問題 C］
ここで，wi は節点 i が観測点のときの重みを表し，
hobs,i は i が観測点のときの観測値を表す．
ただし，i が観測点でないときは wi=0 とし，hobs,i は不要とする．
多くの逆問題は非線形最小二乗問題として定式化されるので非線形最小二乗法を適用で
きる．しかし，大規模問題で adjoint 法（後述）を適
用する場合など，より一般的な非線形最適化問題に
開始
対する非線形最適化手法を用いた方が効率の良い場

k 0 の入力
合がある．
以下，これらの問題を数値的に解く手法について
述べる．その後，特に大規模問題に対応するための
j 0
効率的解法として，adjoint 法による勾配計算法とそ

k i の計算
の具体形，準ニュートン法におけるメモリ節約法に
ついて説明する．



k j 1  k j  k j
j  j 1
2 局所最適化の方法
収束？
ｙ
ここでは非線形問題を線形近似の繰り返しにより

k の出力
数値的に解く方法を示す．局所最小化アルゴリズム
の基本的な流れは Fig.1 に示すとおりである．すなわ


ち，未知ベクトル k の初期推定ベクトル k 0 を与えて，

N
終了

何らかの方法で k j に対する修正ベクトル k j を求め，




k j 1  k j  k j により k j を更新する手続きを収束する

まで繰り返し（j は反復回数），最終の k を解とする
Fig.1
局所最適化アルゴリズム
の基本的な流れ図
土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
増本清（島根大）

ものである．計算効率等は，修正量 k i の求め方の違いにより異なる．本節では以下，修正

ベクトル k j の主要な求め方として，問題 B に対しては(1)最急降下法(2)Newton 法(5)準ニュ
ートン法，問題 C に対しては(3)Gauss-Newton 法(4)Marquardt 法を示してから，(6)その他の
手法に関する概要と(7)制約条件付き問題に対する対応について簡単なコメントを述べる．
（１）最急降下法

E の勾配 E ( Ek ) を用いて k を求める手法を勾配法という．最も簡単な勾配法が最急降
下法であり，


k  E (k )
(6)

により修正ベクトル k を計算する．修正ベクトルの方向ベクトルを  E として算出してか
ら，この方向で  ( 0) を求める． を求める作業は，半直線上の極小点を探索する（直線探
索という）もので，他の最適化法でも使用される．直線探索法には様々な方法（黄金分割
法，フィボナッチ法，ニュートン法など）があるが詳細は文献（例えば嘉納（1987）1））を
参照されたい．最急降下法は単純だが非線形性の強い問題などでは収束性（速度と安定性）
が悪いため，多くの場合，次に示す Newton 法を基礎とした手法が用いられる．
（２）ニュ－トン（Newton）法
Newton 法は多くの局所最小問題の数値的解法の基礎となる．元来，Newton 法は非線形方
  
 

程式 f ( x )  0 を，一階微分値 f (x ) （  f x ，ヤコビ行列という）を用いて線形近似の繰り返
しにより解く方法である．局所最適化問題（問題 B,C）では E  0 （問題 B,C の解の必要
条件）を Newton 法で解くので，E の 2 階微分値  2 E （  Ekk ，ヘッセ行列という，通常は
対称行列）を用いて，


k  2 E 1  E (k )
(7)

とする．初期推定ベクトル k 0 が解に近いとき収束性が良い．しかし，  2 E （またはその逆

行列）の評価が必要であること， k 0 が解から遠いとき収束性が悪い場合があることが弱点
となる．
問題 C に Newton 法を適用すると，勾配ベクトル E の m 番目の要素 Em  Ek m  は，
Em  2 wi ( H i  hobs,i )
i
H i
km
(8)
となり，ヘッセ行列  2 E の m 行 n 列の要素  2 Enm  Ek m k n  は，
 H
H i
2 Hi
H i  hobs,i 
 2 Emn  2  i wi

kn km kn
i  k m

となる．
(9)
土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
増本清（島根大）
Newton 法の弱点を補強するために様々な方法が開発されており，いくつかの主要な方法

を以下に示す．なお，Newton 法は勾配 E のみを使って k を求めることができない点で通
常の勾配法とは異なる．
（３）Gauss-Newton 法
問題 C の最小二乗問題に対して，Newton 法における(9)式の第 2 項を省略して，m 行 n 列
の要素が
 H
H i 
Gmn  2  i wi

kn 
i  k m
(10)
となる行列 G を  2 E の代わりに使用する方法を Gauss-Newton 法という．(9)，(10)式からわ

かるように，線形に近く  2 H / kmkn が無視できる場合や， k j が解に近く H i (k j )  hobs,i が小さ

いとき，ニュートン法の良い近似となる．Gauss-Newton 法では  2 E は不要だが， H k （ヤコ
ビ行列）を評価する必要がある．
なお，後に示す adjoint 法を適用してヤコビ行列を求める方法もあるが，計算効率はあま
り上がらない（観測点数によっては効率的な場合もある）ので，大規模問題に対しては，
Gauss-Newton 法よりも，準ニュートン法と adjoint 法の組み合わせの方が効率的である．
（４）Marquardt 法

初期推定ベクトル ko が解から遠く，非線形性が強いとき，Gauss-Newton 法では良い収束
性を示さないことがある．これを改善するために，(10)式の G の代わりに
G  G  I
(11)
を用いる方法を Marquardt 法という．ただし，ここでは I は単位行列とし， ≧0 である．
=0 のときが Gauss-Newton 法であり，解から遠いときは  を大きくすることにより最急降下
法に近い方法をとる．これにより大域的収束性（初期推定値が解から遠くても収束する性
質）を得ることができ，収束性が改善する．これは Gauss-Newton 法と最急降下法を組み合
わせた方法とみなすことができる．解に近づくにつれて  を小さくしていくことにより，効
率のよい収束性を獲得している．さらに，  の決め方を工夫することにより効率を上げる
方法（修正 Marquardt 法）がある（詳細は中川・小柳(1982）2)を参照のこと）．
（５）準ニュートン法(quasi-Newton method)
Newton 法におけるヘッセ行列を評価せずに，大域的収束性を保障する方法が準ニュート

ン法である．準ニュートン法は，(7)式の  2 E (k j ) の代わりに，適当な M×M 正定値対称行列
土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
増本清（島根大）
Fj（あるいはその逆行列 Bj）を用いる．Fj は，更新公式により更新されるヘッセ行列の近似
行列である．公式導出過程については文献（例えば，今野・山下（1978）2），八巻ほか(1991)3)）
を参照されたい．準ニュートン法の特徴として以下の点が挙げられる．


①  2 E (k j ) の近似行列 Fj を目的関数の一階微分値 E (k j ) で構成する．
② 大域的収束性を得るために Fj の正定値性が保証されている．
③ ヘッセ行列の近似をするために，セカント条件（準ニュートン条件ともいう）




E (k j 1 )  E (k j )  B j (k j 1  k j )
(12)
を満たすように公式が作られている．
なお，準ニュートン法では，近似ヘッセ行列を用いるため，収束性を保証するために，
修正ベクトル方向の直線探索を行うのが一般的である．さらに，大規模問題に対応するた
めの方法（直線探索を最小限に抑える SSVM 法，メモリ節約アルゴリズム）については後
述する．
（６）その他の手法について
上記の手法の他にも数多くの解法があるが，後学のために主要なものを挙げておく．
・Powell のハイブリッド法
問題 C の解法である．関数の微分値（ヤコビ行列）の計算を行わずに済ませる Powell の
最小二乗法を基礎として，様々な手法を組み合わせている意味でこの名がついている．
Gauss-Newton 法と最急降下法の折衷となっている点で Marquardt 法に似ているが，反復とと
もに計算される関数値を用いてヤコビ行列を逐次修正していく．修正ベクトルの大きさな
どに応じた場合分けを必要とする（信頼領域法）など，独特のパラメータがあるが，比較
的収束性が良いとされている（詳細は中川・小柳(1982）2)を参照のこと)．
・共役勾配法
準ニュートン法と並んで，勾配法の一つの有力な解法である．修正ベクトルの方向ベク
トルとして過去の修正ベクトルと共役な方向を選ぶことにより，2 次関数においては有限回
の反復で収束することを保証する方法である．前処理付きなど様々な発展手法があり，線
型方程式の解法としても使用される．
・信頼領域法
直線探索の代わりに，反復毎の修正ベクトルの大きさを信頼領域の考え方で制限するこ
とにより大域的収束性を得る方法である．
・シンプレックス法と内点法
これらは非線形最適化の手法ではないが，線形計画問題の解法として有名なので挙げて
おく．線形計画法の解法であるシンプレックス法の計算量を大幅に削減したアルゴリズム
として内点法が生まれた．内点法は非線形最適化の考え方が使われているが，線形問題を
解くための方法である．
土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
増本清（島根大）
（７）制約条件付き問題への対応について
なお，ここで示した方法により，大抵の無制約条件非線形問題に対応できると思われる
が，さらに発展的な方法については文献を参照されたい．また，制約条件付き問題への対
応については，
・ペナルティ法
・拡張ラグランジェ法
・逐次二次計画法
・変数変換による方法
などがある．実務においては制約条件が必要になることも多いと思われる．最も簡単な制
約条件に対する方法として，ペナルティ法と変数変換により無制約条件化する方法につい
て補足しておこう．

例えば，未知ベクトル k  (k1 , k2 ,  , k M )T に対して，
km  0
(13)
という制約がある場合，次のように目的関数にペナルティ項を加える方法がある．
E'  E   
m
1
k m2
(14)
ここで，   0 である．また，例えば次の変数変換
X m  log km
(15)
により，未知パラメータを X に変換すれば，無制約条件の最適化問題に変換できる．その
際，変換した未知パラメータに関する微分値は合成関数の微分で求めればよい．たとえば
(15)式の変数変換の場合は，
E
E dk m
E



 km
X m k m dX m k m
(16)
となる．

3 adjoint 法による目的関数 E (k ) の勾配の計算法
本節では，目的関数の勾配 E k を効率的に算出する方法を示す．この方法は adjoint method
（随伴法）あるいは adjoint state method （随伴関数法）として知られる方法であり，変分法
により説明することができる．しかしここでは，実用性を考慮して，離散化した問題（問
題 A）に限定することにより，adjoint 法の代わりにラグランジュ乗数法の考え方を用いて
説明することにする．より実際的な複雑で大規模な問題では順問題として離散化した問題
を扱うことが大半であることを考慮すると，変分法で得られた式を離散化して解くよりも，
はじめから支配方程式を離散化した問題 A に限定してから導いた式を利用する方が数値解
析上の精度が高いと考えられるので実用的である．なお，以下の説明を理解するためには，
土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
増本清（島根大）
偏微分法やベクトル・行列の掛け算といった程度の数学的基礎知識があれば十分であろう
と思われる．


状態量 h と同じ I 次元（時空間節点総数次元）のベクトル  を導入してラグランジュ関数
Lを
  
     
L(h , k ,  )  e(h , k )    R(h , k )
(17)

とおく．導入したベクトル   (1 , 1 ,  , I ) は，いわゆる随伴ベクトルあるいは乗数ベクト
ルなどと呼ばれるものである（横ベクトルであることに注意すること）．なお，(17)式の右
  

辺第 2 項は  と R(h , k ) の内積を表している．
ここで，
  
Lh  eh    Rh  0
(18)

（ここでは Rh は i 行 j 列の成分が
Ri
の I×I 行列とする．ヤコビ行列の転置に相当．）
h j
を要請することにより，
  

Ek (k )  Lk (h , k ,  )
(19)

が成り立つことが示される（証明は後述）．これを利用すると，与えられた k に対する勾

配 Ek (k ) は，次の手順により効率的に算出することができる．なお，(18)式は I 個の方程式
から成ることは明らかであろう．


（手順 1）与えられた k を用いて(1)式を解いて， h を求める．





（手順 2） k と h を用いて，(18)式を満たす  を求める．

（手順 3） k と h と  を用いて(19)式より E k を求める．
なお，(1)式を解く（手順 1）が順解析である．(18)式は随伴方程式と呼ばれ，これを解く
（手順 2）は(1)式と同等以下の計算量で解ける上に，(1)式を解く際に用いた手続き（サブ
ルーチン）と共通する部分が多い．したがって，(1)式を解くプログラムができていれば，
(18)式を解くためのプログラム作成はそれほど難しくない．また，（手順 3）の計算量は（手
順 2,3）に比べると無視できるほど小さく，プログラム作成も容易である．






つぎに， k が与えられたとき， h と k が (1)式を満たし，かつ， h と k と  が(18)式を満た
しているならば，(19)式が成立することを示す．
（仮定）(1)式，(18)式，（結論）(19)式
（証明）


h と k が (1)式を常に満たすとき，
土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
増本清（島根大）
 
 
 
dR  Rh  dh  Rk  dk  0

が成り立つ（ R の各成分について全微分が 0 ということ．）ので，
 

Rh  H k  Rk  0
(20)
が成り立つ（右辺は全成分 0 の I×M 行列）．したがって，

Ek  ek  eh  H k

（E の定義（(3)式）より明らか）

  
 ek     Rh  H k


 
 ek     Rh  H k
 
（∵(18)式より eh    Rh ）

 
 ek    Rk



（∵(20)式より  Rh  H k  Rk ）
 Lk
（L の定義（(17)式）より明らか）
（証明終わり）


 
いわゆるラグランジュ乗数法では，目的関数 e(h , k ) を最小化するような h と k を求めるた



めに， h ， k ，  に関する連立方程式
   
L (h , k ,  )  0
(21)
   
Lh (h , k ,  )  0
(22)
   
Lk (h , k ,  )  0
(23)

を解く．ここでは，(1)式が(21)式に，(18)式が(22)式に相当する．勾配 E k が 0 になること（す
なわち(5)式）が(23)式に相当する．
以上をまとめると，本節で示した勾配 Ek の計算法を準ニュートン法などの勾配法と組み


合わせて，勾配 E k が 0 となるような k を求めれば逆問題が解けたことになる．このことは，
ラグランジュ乗数法における連立方程式(21)，(22)，(23)を解いたのと同じ結果になる．ま
た，勾配計算に要する計算量は，未知パラメータ数 M に関わりなく順問題を解く計算量の
2 倍程度以下となる．
土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
4 adjoint 法の適用例 ―
増本清（島根大）
具体的な式展開の例
ここでは例として，次の放物型微分方程式
S
h   h 
 K   Q
t x  x 
(24)
が支配する現象について具体的な式展開を示そう．問題 A の目的関数 e および(1), (18),(19)
式の具体的な例を示せば十分であろう．
(24)式を時空間について差分法（陰解法）で離散化すると，例えば次のような差分方程式
 S K i 1 2 K i 1 2  n S i n1 K i 1 2 n K i 1 2 n
  hi   hi 
Rin   i 

 hi 1 
 hi 1  Qin  0
2
2 
2
2
t
x 
x 
 t x  x  
(25)
が得られる．ただし，時空間のステップ幅 Δt，Δx は，簡単のため時間ステップ n および空
間格子点番号 i によらず一定とした．また，初期条件 hi0 は既知とし，境界条件は閉境界，す
なわち
K1 2  K I 1 2  0
(26)
の場合について考えることにする．これにより(25)式は，
i  1,2,  , I
（I：空間格子点数）
n  1,2,  , N
（N：時間ステップ数）
で成立することになる．なお，(25)式は初期・境界条件を内包していることに注意されたい．
この(25)式が(1)式に相当する．なお，前節における I は，本節の I×N に相当することにも
注意されたい．
n
次に目的関数は，対象となる 1 次元対象領域内の観測点における観測値 hobs
と(25)式を解
,i
いて得られる計算値 hin の重み付き残差 2 乗和

n
E   win hin  hobs
,i

2
(27)
i ,n
とする．ただし，時空間ステップ（i,n）が観測点に当たらない（ hin が観測点に対応してい
ない）ときは win  0 とする．また，未知パラメータとして，各隣接格子点間の K（拡散係数
や伝導率に相当する物性値）からなるベクトル

K  K 3 2 , K 5 2 ,  , K i1 2 ,  , K I 1 2 
(28)

を例として考える．逆問題は(27)式の目的関数 E を最小にする(28)式のベクトル K を求める
土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
増本清（島根大）
非線形最小問題となり，これが問題 B の目的関数に相当する．
ラグランジュ関数 L を
L  E   in  Rin
(29)
i ,n
とおくことにより，随伴方程式
b
L
E
b Ra

  a  n  0
hin hin a ,b
hi
(30)
の具体形は，
K i 1 2 n
 Si K i 1 2 K i 1 2  n Si n1 K i 1 2 n
n
 
  i   i 

 i 1 
 i 1  2win hin  hobs
,i  0
2
2 
2
2
t
x 
x 
 t x  x  


(31)
となる．また，
iN 1  0
（ i  1,2,  , I ）
(32)
とした．これにより，(26)式の境界条件も併せて考慮すると，(31)式は i=1，,i=I，n=N にお
いても成立することになる．この(31)式が(18)式の具体形の例となる．
これより，(31)式は終末条件(32)式から，逆時間方向に順解析（（25）式）と同じ係数行
列（非対称の場合は転置行列）の方程式を時間ステップ毎に解くことにより， in
（ i  1,2,  , I ， n  1,2,  , N ）を求められることがわかる．
次に勾配計算式の具体形は，
Ran
L
E

  na 
K i 1 2 K i 1 2 a ,n
K i 1 2

Rin
Rin1 
 0    in 
 in1 

K i 1 2
K i 1 2 
n 
 n  in hin1  hin 

   i 1

x
x 
n 
(33)
となる．これが(19)式の勾配計算式に相当する．
より複雑なモデル，多様な未知量，多様な目的関数形などについても同様の手順で式展
開を行うことにより，問題に応じた adjoint 法による具体的な勾配計算式が得られる．
土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
増本清（島根大）
5. 大規模無制約問題に対する準ニュートン法
非線形最小問題として定式化された逆問題（問題 B）を，目的関数の勾配（(4)式）を繰
り返し計算して解く方法を勾配法という．大規模問題についても adjoint 法により計算効率
のよい勾配計算が可能になることはすでに示した．ここでは，得られた勾配を用いて，未
知パラメータ数の多い大規模問題を解く方法として，勾配法の一種である準ニュートン法
の内 SSVM 法(self scaling variable metric method)5)について，計算速度と計算機記憶容量の両
面で効率のよい方法を示す 6)7)．
（１）SSVM 法
準ニュートン法では，修正ベクトルの方向ベクトル（  B j E ）の大きさを調節するため
に順解析計算の繰り返しを要する直線探索（1 変数関数の最小点の探索）が必要だが，一部
の更新公式では目的関数値の減少（Ej+1<Ej）を保証するだけの粗い直線探索でも十分収束す
ることがわかっている．さらに，スケーリングファクターωj をヘッセ逆行列 Bj に反復毎に
乗じることにより直線探索をほとんど行わない次の SSVM 法公式族が一層計算効率の良い
方法である．
 
 T

B j y j y Tj B j
 T    T  s j s j

B j 1   j B j   T    j y j B j y j v j v j   T  
 s y
y j Bj y j

 j j


(34)
ただし，


wj
sj

v j  T   T 
y j wj s j y j
(35)

 

s Tj y j
s Tj E (k j )
 
 j  1   j   T    j
y j wj
E (k j )T w j
(36)
 j ,  j  0,1
(37)

y j  E j 1  E j (38)

 
s j  k j 1  k j
(39)


wj  Bj y j
(40)
土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
増本清（島根大）
以上からわかるように，準ニュートン法の反復公式をこのまま使用する場合は，対称行
列であるヘッセ行列の逆行列の近似行列 Bj（以下，ヘッセ逆行列と略す）を反復計算毎に
記憶しておく必要があることがわかる．このため，本手法のコード化に際して，膨大な未
知パラメータ数に対応するために，次節で述べる記憶容量節約アルゴリズムを適用するこ
とができる．
（２）記憶容量節約法による SSVM 公式の変形 6)
SSVM 法におけるヘッセ逆行列の計算機記憶容量はほぼ M2 に比例するため，未知パラメ
ータ数の多い大規模問題ではこれを節約する必要が生じる．ここでは，SSVM 公式でこの記
憶容量を節約する方法を示す 6)．これにより，計算速度を損なわずに，この部分の記憶領域
をアルゴリズムの修正により M（未知パラメータ数）に比例する容量にすることができる．
SSVM 公式(34)は次のように変形できる．
B j 1   j B j  1j s j sTj  2j v j vTj  3j w j wTj
(41)
ただし，
1
s yj
(42)
2j   j j yTj w j
(43)
1j 
3j 
T
j
j
(44)
y tj w j
計算開始時に与える正値対称行列 B0 として単位行列 I を与えるものとして，(41)式の漸
化式を解くと，ヘッセ逆行列 Bj は(45)式のように表すことができる．
j 1
j 1
j 1
r 0
r 0
r 0
B j  aI   br sr srT   cr vr vrT   d r wr wrT
(45)
ただし，
j 1
a   l
(46)
l 0
j 1
br  1r  l
(47)
l  r 1
j 1
cr  2r  l
(48)
l  r 1
j 1
d r  3r  l
l  r 1
(49)
土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
増本清（島根大）
(45)式の右辺の形からわかるように，Bj はベクトルのダイアデック（縦ベクトル×横ベクト
ル）と単位行列の線形和で構成されている．これにより，未知パラメータの修正ベクトル

の方向ベクトル d j は，

d j   B j E j
j 1


 j 1 
 j 1

 
 aE j   br srT E j sr   cr vrT E j vr   d r wrT E j wr 
r 0
r 0
r 0








(50)
として計算できることがわかる．この(50)式による計算では，ベクトルの内積計算を先に行
うために未知パラメータ次元の行列を記憶する必要がない．
さらに，
(41)式の右辺の Bj に(45)
式を代入することにより，
j 1
 

 j 1 
 j 1

 
w j  ay j   br srT y j sr   cr vrT y j vr   d r wrT y j wr 
r 0
r 0
r 0








(51)
とすれば全アルゴリズムの中から，Bjとダイアデックを消去することができる．これにより
対称行列Bjないしダイアデックからなる行列の記憶領域を必要としないアルゴリズムを導
  
くことができた．ただし，反復ステップごとに算出される3つのベクトル s , v , w （3M個の数
値)をすべて記憶することが必要になる．また，ここで示したようにSSVM法においては反
復毎に変化するスケーリングファクターωjの存在があるために，a,br,cr,drの計算を反復毎に
行うことが必要になるが，この分の計算量は小さくできる6)．
参考文献
1)
嘉納秀明（1987）：システムの最適理論と最適化，コロナ社
2)
中川徹，小柳義夫（1982）：最小二乗法による実験データ解析－プログラム SALS，UP 応用数学選
書７，東京大学出版会．
3)
今野
浩，山下
浩（1978）：非線形計画法，日科技連．
4)
八巻直一，宮田雅智，本郷茂，高橋悟，矢部博，内田智史（1991）:パソコン FORTRAN 版非線形最適
化プログラミング，ASNOP 研究会編，日刊工業新聞社．
5)
Oren, S.S. and D.G. Luenberger (1974):Self-scaling variable metric (SSVM) algorithm, Part1: Criteria and
sufficient conditions for scaling and class of algorithms, Management Science, 20, 845-862.
6)
増本
清（2004）：地下水理逆解析における計算機記憶容量節約アルゴリズム，応用力学論文集，７，
191-199．
7)
西垣
誠ほか（2009）：地下水のトレーサー試験，技報堂出版．

Download Report