統計数理（2022）
第 70 巻 第 2 号 195–208©2022 統計数理研究所

特集「データ同化の方法」［研究ノート］

共役ベクトルと BFGS 公式を用いた
解析誤差共分散行列の導出
丹羽 洋介1,2 ・藤井 陽介2,3（受付 2022 年 1 月 4 日；改訂 6 月 30 日；採択 6 月 30 日）

要旨

データ同化や逆解析の問題において，4 次元変分法は有効な手法であるが，最適解の誤差が自動的に得られないという欠点がある．本稿は Niwa and Fujii（2020）で示された，4 次元変分法を用いた場合の解析誤差共分散行列の推定手法について解説を行った．本手法では，最適解の探索アルゴリズムとして用いられる BFGS 公式を用いた準ニュートン法を採用しているが，厳密な直線探索やアンサンブル法，直交化を新たに導入することで，BFGS 公式の計算に必要なベクトルの共役性を保ちつつ数を増やして解析誤差共分散行列の推定精度を向上させている．
本稿では，BFGS 公式において前処理を施した場合に観測と同じ数の反復回数（または BFGS
公式で必要なベクトルのペア数）で解析的に解析誤差共分散行列が求められることを解説するのアルゴリズムを詳述する．さらに，本手法を CO2 の逆解析とともに，Niwa and Fujii（2020）問題に適用した際の結果についても紹介する．
キーワード：解析誤差共分散行列，データ同化，逆解析，4 次元変分法，BFGS 公式，
準ニュートン法．
1.

はじめに

4 次元変分法（アジョイント法）は問題規模の大きいデータ同化・逆解析の分野において有用な手法であり，広く用いられている手法である（例えば Sasaki, 1969）．しかしながら，アンサ
ンブル・カルマンフィルターなどの手法とは異なり，最適解の誤差（解析誤差）
が自動的には得られないという欠点がある．本論文では 4 次元変分法を用いた場合に，どのように解析誤差を得ることができるか，一つの手法を紹介する．
観測データを y o とし，あるパラメータ x からこの観測データに対応する値を計算するモデ
ル演算を M (x) と定義する．ここで，x, y o はそれぞれ，解析対象とする時空間内のあらゆる時間・場所の制御変数，観測データを含むものとし，要素数は n, m とする．また，モデル演
算子 M (.) は時間発展演算子と観測演算子が含まれる．4 次元変分法では，x の初期推定値 xpriが与えられたとき，条件付確率密度関数 p(x|xpri , y o ) が最大となる解 xpos を求める（maximum
a posteriori estimate: MAP 推定）．

1

国立環境研究所 地球システム領域：〒 305–8506 茨城県つくば市小野川 16–2

2

気象庁気象研究所 全球大気海洋研究部：〒 305–0052 茨城県つくば市長峰 1–1

3

統計数理研究所：〒190–8562 東京都立川市緑町 10–3

196

(1.1)

統計数理

第 70 巻 第 2 号 2022

xpos = arg maxn p(x|xpri , y o )x∈R

ここで，ベイズの定理を用いると，(1.2)

p(x|xpri , y o ) =p(xpri , y o |x)p(x)p(xpri , y o )であり，さらに，パラメータ x が与えられたときに xpri と y o が独立に生起するとすると，(1.3)

p(x|xpri , y o ) =p(xpri |x)p(y o |x)p(x)p(xpri , y o )が得られる．式（1.3）の分母は x に依存しないので式（1.1）は(1.4)

xpos = arg maxn [p(xpri |x)p(y o |x)p(x)]x∈Rと変形することができる．p(x) は初期推定値も観測もないときの x の確率密度関数であり，一般に得ることは容易ではないため，ここでは一様分布と仮定すると，(1.5)

xpos = arg maxn [p(xpri |x)p(y o |x)]x∈Rとなり，MAP 推定は最尤推定（maximum likelihood estimate）に帰着する（露木・川畑, 2008）．
4 次元変分法では式（1.5）
をもとにして定義した評価関数の勾配をアジョイントモデルとよばれるモデルを用いて計算し，反復しながらその勾配がゼロとなる最適解を探索する．ここでは簡単のため，p(xpri |x) と p(y o |x) がガウス分布であると仮定すると，評価関数は一般的に対数で
定義され，
1
1
(1.6)
J(x) = (x − xpri )T B −1 (x − xpri ) + (M (x) − y o )T R−1 (M (x) − y o )22と表すことができる．4 次元変分法では，この評価関数 J を最小とする解を求める．なお，Bおよび R はそれぞれ初期推定誤差とモデル - 観測間の差の共分散行列を表している．さらに，
M (.) が十分に線形であるとみなせる場合には，式（1.6）は接線形モデルを表す行列 M を用いて以下のように書き直すことができる．
1
1
(1.7)
J(δx) = δxT B −1 δx + (M δx − do )T R−1 (M δx − do )2
2
なお，δx := x − xpri であり do := y o − M (xpri ) である．式（1.7）の評価関数が最小となる解は(1.8)

δxpos = arg minn J(δx) = (B −1 + M T R−1 M )−1 M T R−1 doδx∈Rであり，このとき，評価関数の勾配 ∇δx J は(1.9)∇δx J = B −1 δx + M T R−1 (M δx − do ) = 0を満たす．問題の規模が大きい場合，式（1.8）を解析的に解くことは不可能なため，実際の 4 次元変分法の計算では M T で表されるアジョイントモデルを用いて評価関数の勾配を計算し，最適解の探索を行う．この時，最適解の誤差共分散行列（解析誤差共分散行列）
Pa は必要なく，同時に求められることはないが，実際には Pa は評価関数のヘッセ行列−1+ M T R−1 M ）の逆行列と等しく，（A = ∇δx ∇T
δx J(δx) = B(1.10)

Pa = A−1 = (B −1 + M T R−1 M )−1と表すことができる．問題規模が大きい場合には式（1.8）と同様に式（1.10）も解析的に解くこと

共役ベクトルと BFGS 公式を用いた解析誤差共分散行列の導出

197

は不可能なため，Pa を推定する手法が必要となってくる．なお，式（1.7）（1.10）では，M (.) で–はなく線形のモデル M が使われていることに注意が必要である．モデルが線形であるために，評価関数（式（1.7））
が 2 次関数であり，極小値が一つしか存在せず，ヘッセ行列が定数行列となっている．したがって，解析誤差共分散行列はモデルが線形である場合にしか定義することができない．本稿では，この解析誤差共分散行列を推定する手法を提案しているが，その適用範囲は，モデルが線形であり，かつ，p(xpri |x) と p(y|x) がガウス分布と仮定した線形問題に限っていることに注意されたい．
これまで 4 次元変分法をフレームワークとした Pa の推定には様々な手法が提案されており，
法，モンテ・カルロ法，ランチョス法，Broyden-Fletcher-Goldfarbランダム化（randomization）
法などがある（Rabier and Courtier, 1992; Fisher and Courtier, 1995; Chevallier
Shanno（BFGS）
et al., 2007）．他にも second-order adjoint 法を用いる手法も存在するが（Ito et al., 2016），制
御変数の最適化計算では用いない second-order adjoint モデルを新たに開発する必要があり，モデルが大規模で複雑である場合には開発コストがかかる．ランダム化法やモンテ・カルロ法はその確率論的な性質からサンプル数に精度が大きく依存し，特に共分散（Pa の非対角要素）の推定精度に難があることが指摘されている（Fisher and Courtier, 1995）．一方で，ランチョス
法および BFGS 法は，それぞれ，4 次元変分法の最適解探索のアルゴリズムとして用いられる
共役勾配法や準ニュートン法と密接に関係，またはアルゴリズムの一部として用いられているため，最適解探索のアルゴリズムを応用することで Pa を推定することが可能であるが，その．
推定精度を上げるためには，さらなる追加の手法が必要となってくる（Bousserez et al., 2015）
そこで筆者らは，制御変数の最適化に用いる既存の BFGS 法を使用しつつ，実装が簡易なベクトルの共役化やアンサンブルを新たに導入することにより，高精度に Pa を導出する手法を．本論文では，その手法について解説を行う．なお，これから開発した（Niwa and Fujii, 2020）の記述においては，簡単のために δx の δ は省略するとともに，慣例的な表記に従って，ベク
トル d は制御変数の探索方向，ベクトル y は評価関数の勾配のインクリメントとして表すこととする．
2.

BFGS 公式

4 次元変分法で用いられる最適解探索のアルゴリズムの一つとして準ニュートン法がある．
この準ニュートン法は，ある k 回目の反復計算における制御変数の更新式(2.1)

xk = xk−1 + αk dk−1

で必要な探索ベクトル dk−1 を求める際に，ヘッセ行列の近似逆行列を用いる手法である（なお，αk はその探索幅で J(xk−1 + αk dk−1 ) が最小となるように決められる）．具体的には，探索ベクトルはヘッセ行列の近似逆行列 Hk と評価関数の勾配ベクトル gk を使って，(2.2)

dk = −Hk gkとして決定される．この時に必要な Hk の求め方には様々な手法が存在するが，その一つに
BFGS 公式があり，以下のような式で書くことができる．
(2.3)

T
Hk = (I − ρk pk ykT )Hk−1 (I − ρk yk pTk ) + ρk pk pk

なお，ここで，pk := xk − xk−1 , yk := gk − gk−1 , ρk := 1/(ykT pk ) であり，I は単位行列を表す．
また，最初の反復で必要となる H0 は任意の正定値対称行列を用いることができる．実際の（2.3）の計算では，y および p のベクトルのみを使った再帰的な計算（two-loop recursion（2.2）と統計数理

198

第 70 巻 第 2 号 2022

scheme）
を用いることで，Hk をメモリに保持することなく dk を計算することができる（例え．
ば Nocedal and Wright, 2006）
ここで注目すべき点は，Hk がヘッセ行列の近似逆行列であるということから，Pa の近似行列と見なすこともできるということである．したがって，式（2.2）で gk を (1, 0, . . . , 0)T ,(0, 1, . . . , 0)T , . . . , (0, 0, . . . , 1)T といったベクトルに置き換えてやれば，Pa のそれぞれの列を推定することが可能となる．しかし，その推定精度は Hk の計算に必要なベクトル p, y や初期推では，アンサンブル手法を使ってベクトル (p, y)
定行列 H0 に依存する．Niwa and Fujii（2020）のペア数を増やしつつ，(p, y) が本来満たすべき共役性を達成することで Pa の推定精度が大幅に向上することを実証した．以下の章ではその手法について解説していく．
3.

共役性であるとき，共役性が満たされているとは，以下問題が線形（評価関数が下に凸な 2 次関数）の式が成り立っていることを意味する．
(3.1)

pT
i Apj = 0,

for all i = j

または，線形の場合には Apj = yj であることから
(3.2)

pT
i yj = 0,

for all i = j

BFGS 公式を用いた反復計算において，この共役性が満たされていれば，制御変数の数 (= n)の回数の反復で解が収束して解析解と等しくなり，また，そこでの Hk もヘッセ行列の逆行列と等しくなる（Hn = A−1 ）ことが知られている（Nocedal and Wright, 2006）．逆にいえば，共役性が満たされていなければ，n 以上の回数が解の収束に必要となる．
実際には前処理を行うことで，この収束速度は速めることができる．最適化計算において前の手法がある．
処理は一般的に使われるものであるが，最もよく使われるものに Lorenc（1988）の前処理では，制御変数がこの Lorenc（1988）(3.3)

x̃ = B −1/2 xと変換される．この時，変換されたヘッセ行列は(3.4)

Ã = I + B 1/2 M T R−1 M B 1/2 = I + XX T
X := B 1/2 M T R−1/2と書くことができる．通常のデータ同化または逆解析の問題では，n > m であるため，このヘッセ行列は 1 以外の固有値はたかだか m 個しか存在せず，収束性を速めることができる．な
おこの前処理をした後の BFGS 公式において，そのヘッセ近似逆行列の初期値 H̃0 を最もシン（H̃0 = I ），前処理前の BFGS 公式において H0 = B
プルな正定値対称行列である I とした場合とした場合と同等となる．したがって，通常の BFGS 公式において H0 = B とすれば，陰に式（3.3）の前処理をしたこととなり，収束を速めることができる．このような前処理をした BFGS
公式においては，観測数 m の数だけの共役性を満たす (y, p) のベクトルのペアさえあれば，それらによって得られるヘッセ近似逆行列 Hm が解析的に得られる行列，すなわち Pa と等しくなる．
このことは以下のように証明することができる．まず，共役性が満たされていれば，変換
後も
(3.5)

T
p̃T
i ỹj = p̃i Ãp̃j = 0,

for all i = j

共役ベクトルと BFGS 公式を用いた解析誤差共分散行列の導出

199であるため，
T(3.6) H̃k ỹj = [(I − ρ̃k p̃k ỹkT )H̃k−1 (I − ρ̃k ỹk p̃Tk ) + ρ̃k p̃k p̃k ]ỹj

= (I − ρ˜k p̃k ỹkT )H̃k−1 ỹjTT= (I − ρ˜k p̃k ỹkT )(I − ρ̃k−1 p̃k−1 ỹk−1) · · · (I − ρ̃j+1 p̃j+1 ỹj+1)H̃j ỹj ,for j ≤ kが導かれる．ここで，ỹk := B 1/2 yk , p̃k := B −1/2 pk , ρ̃k := 1/(ỹkT p̃k ) である．また，(3.7)

TH̃j ỹj = [(I − ρ̃j p̃j ỹjT )H̃j−1 (I − ρ̃j ỹj p̃Tj ) + ρ̃j p̃j p̃j ]ỹj = p̃j ,for j ≤ k

により，式（3.6）
はH̃k ỹj = p̃j ,(3.8)

for j ≤ kとなり，シンプルな関係性が得られる（Dennis, Jr. and Moré, 1977; Fisher and Courtier, 1995）．
さらに X の固有ベクトル {u1 , u2 , . . . , um } のいずれとも直交する任意のベクトルを ū とす
ると，
(3.9)

T
H̃k ū = [(I − ρ̃k p̃k ỹkT )H̃k−1 (I − ρ̃k ỹk p̃Tk ) + ρ̃k p̃k p̃k ]ū= (I − ρ̃k p̃k ỹkT )H̃k−1 ūT) · · · (I − ρ̃1 p̃1 ỹ1T )H̃0 ū
= (I − ρ̃k p̃k ỹkT )(I − ρ̃k−1 p̃k−1 ỹk−1= ūという関係が得られる．なお，ここで(3.10)ỹk ∈ span{u1 , u2 , . . . , um }, p̃k ∈ span{u1 , u2 , . . . , um },for any k ≤ m

を満たしており，(3.11)

span{u1 , u2 , . . . , um } = span{p̃1 , p̃2 , . . . , p̃m }であることから（これらの詳しい証明は Niwa and Fujii, 2020 を参照のこと），1 ≤ k ≤ m につTいて p̃Tk ū = 0 かつ ỹk ū = 0 であること，さらに H̃0 = I であることを用いた．
ここで，制御変数の空間の任意のベクトル b は式（3.11）により(3.12)b=m

ai p̃i + ū = P̃ a + ūi=1P̃ := (p̃1 , p̃2 , . . . , p̃m ), a := (a1 , a2 , . . . , am )Tと表すことができ，その b は以下の式を満たす．
(3.13)

H̃m Ãb = H̃m (Ỹ a + ū) = P̃ a + ū = b

ここで Ỹ := (ỹ1 , ỹ2 , . . . , ỹm ) であり，式
（3.4）
，
（3.8）
，
（3.9）
，また，ỹ = Ãp̃, X T ū = 0 であるはことを用いた．さらに式（3.13）(3.14)(H̃m Ã − I)b = 0と書き換えることができ，この式が任意のベクトル b に対して成り立つことから，H̃m = Ã−1であり，前処理前の行列においても(3.15)

Hm = A−1 = Pa

200

統計数理

第 70 巻 第 2 号 2022が成り立つ．以上のことから，共役性が保たれていれば，m 個の (y, p) のペアによって Pa が解析的に導出できることが示されたが，m 個だけあれば良いとはいえ，観測データの数は通常，可能な反復回数よりも大きい．したがって，(y, p) のベクトルのペアを如何に増やすか，また，それらの (y, p) のベクトルについて，どうやって共役性を保たせることができるか，これらの問題に対応した手法を次節で紹介する．
4.

共役性のある (y, p) ベクトルの導出

上記では解析誤差共分散行列 Pa の推定のために BFGS 公式を用いることを提案しているが，そこで必要な (y, p) の作成には必ずしも BFGS 公式を用いた準ニュートン法を用いる必然
では BFGS 公式を用いた準ニュートン法をベースと性はない．しかし，Niwa and Fujii（2020）
したアルゴリズム POpULar（Preconditioned Optimizing Utility for Large-dimensional analyse）（Fujii and Kamachi, 2003; Fujii, 2005）を用いて (y, p) のペアを作成しており，ここで，新たに厳密な直線探索（exact line search）
を導入して (y, p) の共役性が保たれるようにしている．さらに，(y, p) のペア数を増やすためにアンサンブルを導入し，アンサンブルメンバー間の (y, p)についても共役性が保たれるように直交化を行っている．
4.1 厳密な直線探索
問題が線形であり，評価関数が下に凸な二次関数である場合，式（2.1）において J(xk−1 +αk dk−1 ) が最小となる αk は解析的に
αk = −(4.1)

dT
k−1 gk−1dT
k−1 Adk−1と表すことができ，この αk を用いた直線探索を厳密な直線探索（exact line search）と呼ぶ．問
題が線形な場合，この厳密な直線探索が BFGS 公式を用いた準ニュートン法で行われれば，得られる (y, p) ベクトルが共役性を持つことが知られている（Dennis, Jr. and Moré, 1977）．しは A が基本的に未知であり，また，既知であったとしても，その行列のかしながら，式（4.1）
から容易に計算することはできない．しかし，新たな変数
サイズの大きさ（n × n ）(4.2)

fk = ek−1 + M T R−1 M dk−1を導入することにより，式
（4.1）
は
(4.3)

αk = −dT
k−1 gk−1dT
k−1 fkとなり，ベクトルのみで計算することが可能となる．なお，ここで，e := B −1 d である．式（4.2）式に M と M T があるため，f の計算には接線形モデルとアジョイントモデルによる計算が必要となっているが，一方で評価関数の勾配 g は(4.4)

gk = gk−1 + αk fkという式で更新することができるため，通常の最適化計算において必要であった g の導出における M および M T の計算は必要ではなくなり，全体的な計算量は厳密な直線探索を導入しても変わらないようになっている．
で提案されている．ただし，
実はこの f を用いた厳密な直線探索は Derber and Rosati（1989）
その最適化手法は準ニュートン法ではなく共役勾配法が用いられている．この共役勾配法での前処理を処しているものの，漸化式を用いることによって，は，先ほど述べた Lorenc（1988）

共役ベクトルと BFGS 公式を用いた解析誤差共分散行列の導出

201

前処理に必要な B 1/2 や B −1 の計算は必要ないというのがもう一つの特徴である．POpULar
は，非線形の場合にも対応できるように，この Derber and Rosati（1989）の手法を BFGS 公式
を用いた準ニュートン法に拡張したもので，直線探索も厳密なものではないものに置き換えられている．しかし，Niwa and Fujii（2020）はこの POpULar をベースとしつつも，直線探索を（Linear
上記の厳密なものに戻して共役性のある (y, p) ベクトルのペアを作成するようにした
POpULar）．そのアルゴリズムは以下のようになる．
Algorithm 1 Linear POpULar
x0 ← 0;g0 ← −M T R−1 d0 ;d0 ← −Bg0 ;e0 ← −g0 ;k ← 1;
repeatfk = ek−1 + M T R−1 M dk−1 ;xk ← xk−1 + αk dk−1 ,ここで，αk は厳密な直線探索（式 4.3）で得られる;
gk = gk−1 + αk fk ;dk ← −Hk gk ,ek ← B −1 dk ,ここで，two-loop recursion scheme を用いた BFGS 公式の計算により dk のみならず ek の計算も行う．
その際，B −1 の計算を行わずとも ek を計算することができる;k ← k + 1;
until convergence

なお，BFGS 公式において，最初のヘッセ行列の逆行列の推定（H0 ）に B を用いることでLorenc（1988）の前処理を導入しているが，上記のアルゴリズムで見て分かる通り，B 1/2 やB −1 の計算は必要ない．
問題が線形で厳密な直線探索が使われている場合，共役勾配法と BFGS 公式を用いた準ニュートン法は同一であることが知られているが（例えば Nazareth, 1979），共役勾配法は丸め
誤差の影響を受けやすいため（Golub and van Loan, 1996），安定して共役性の保たれた (y, p)のペアを得るためには，準ニュートン法を用いた方が適している．また，ここで注意しておきたいのが，オリジナルの POpULar の BFGS 公式で使われており，また，一般的にもよく使わでは用いていないという点である．メモリ節約とれているメモリ節約は Niwa and Fujii（2020）
は式（2.3）の H の更新式において，全ての (y, p) ベクトルを用いるのではなく，直近のあらか
じめ定めた数のペアのみを用いることであるが，ここでの目的は，共役性の保たれた (y, p) ベ
クトルのペアをできるだけ多く作るということであるため，わざわざベクトルの数を減らすメ
モリ節約は行わない．なお，近年の計算機性能の向上により，全ての (y, p) のデータを保持することは容易であり（必ずしもメモリに載せる必要はなく，ディスクに格納できれば良い），全ての (y, p) のデータを保持することに大きな障害は生じない．
最終的な反復回数を K ，モデル計算にかかるコストを C とすると，Algorithm 1 全体の計算量は O(nK 2 ) + O(KC) と表すことができる．第一項は反復回数の二乗に比例しているが，大気・海洋などの問題でモデル計算のコストが圧倒的に大きい場合には，第二項が卓越し，反復はた
回数に比例して計算量が増加する．また，理論的に必要とされる (y, p) の数（= 反復回数）
かだか m であるが，n  m であったとしても，相当の反復回数が必要となる場合がある．特

202

統計数理

第 70 巻 第 2 号 2022にモデル計算が大規模な場合には，アルゴリズムが逐次的であるために，実行可能な反復回数がかかりすぎる）．そこで次節では，実行可能な反に限界が生じる（計算時間（wall-clock time）
復回数で Algorithm 1 の計算を打ち切る一方，同様の計算を複数，同時に行うことで，効率よく (y, p) の数を増やすことを考える．
4.2 アンサンブルと直交化
BFGS 公式の式（2.3）
を見て分かる通り，H は (y, p) の順序には関係ない．したがって，
Ap = y である (y, p) のベクトルの数を何らかの方法で増やせばよく，そのため，Niwa andでは，準ニュートン法による最適化計算を同時に複数行うアンサンブル手法を用い
Fujii（2020）
た．先行研究である Bousserez et al.（2015）でも同様に，BFGS 公式とアンサンブルを組み合わせた手法を提案しているが，Niwa and Fujii（2020）では，さらに，得られた (y, p) ベクトルに直交化を施すプロセスを入れている．これは，それぞれのアンサンブルメンバーは独立に最適化計算を行うことから，メンバーをまたいでは (y, p) の共役性が保たれないためで，直交化は，これらのベクトルを変換して共役性を満たすようにする．この直交化は以下のようにして行う．
まず得られた全ての (y, p) ベクトルをそれぞれ並べて以下のような行列を作る．
(4.5)

P = (p1 , p2 , p3 , . . . , pN )(4.6)

Y = (y1 , y2 , y3 , . . . , yN )ここで N は全てのアンサンブルメンバーから得られた (y, p) ベクトルの数である．次に，Pと Y をかけて行列 Q を作成する．
(4.7)

Q = P TY

pT y = pT Ap であり A がそもそも対称行列であることから，Q も対称行列であり，次のように固有値分解することができる．
(4.8)

Q = U LU T

ここで U は固有ベクトルからなる直交行列，L は固有値を対角成分にもつ対角行列である．なお，N ≤ m であり，Q は N × N の次元である．m が十分に小さい場合には（実際に大気・海洋，式（4.8）の固有値分解は容易などのデータ同化・逆問題などでは，n  m であることが多い）に実行可能である．この固有値分解において，理論的には Q は正定値行列であるために，得られる固有値も全て正となるべきだが，実際の計算では計算誤差により非常に小さい，または負の固有値が計算されてしまう．そこで Niwa and Fujii（2020）では，計算を安定化させるため，最大の固有値と比して 1.0 × 10−6 より小さい固有値を除外している．その除外された後の固有ベクトル，固有値からなる行列をそれぞれ，U  , L とし，それらを使って P と Y を以下のように変換する．
(4.9)

P  = (p1 , p2 , p3 , . . . , pN  ) = P U  L−1/2(4.10) −1/2
Y  = (y1 , y2 , y3 , . . . , yN) = Y U L

ここで，N  は除外後の有効な固有値の数であり，U  は N × N  , L は N  × N  の次元となっている．なお，変換後でも(4.11)P T Y  = I

共役ベクトルと BFGS 公式を用いた解析誤差共分散行列の導出

203が成り立ち，共役性(4.12)

TpTi yj = pi Apj = 0,

for all i = jが満たされていることに注目されたい．この変換後の (y  , p ) を BFGS 公式の計算に用いるが，単一の反復計算で得られる数よりもはるかに多くの数で，かつ，共役性が保たれているため，Pa を高精度に推定することが可能となる．
なお，アンサンブルで得られた複数の解析値のばらつきを解析誤差とするのではなく，単純に (y, p) のベクトルのペア数を増やせれば良いので，アンサンブルメンバーの作成はシンプルで表せられるように，Pa は誤差共分散行列 B と R，またモデルの計算 Mにできる．式（1.10）にしかよらないため，これらの行列さえ共通であれば，特段，摂動の与え方に制約はない．そでは，単純にランダムに生成した do (= y o − M (xpri )) を入力とし
こで，Niwa and Fujii（2020）て用いている．
Niwa and Fujii（2020）の手法では，4.1 節の厳密な直線探索を使った準ニュートン法の計算をアンサンブルで行い，そこで得られた (y, p) ベクトルを直交化して BFGS 公式を計算している．最終的に直交化するため，理論的には各アンサンブルメンバーの計算で厳密な直線探索を用いる必然性はないが（Ap = y さえ満たされていれば良い），あらかじめ，厳密な直線探索
で各メンバー内の (y, p) の共役性が保たれていた方が安定的に Pa を推定できることを実際の実験で経験している．また，直線探索を厳密にすることで発生する追加の計算コストはないた
め，安全のためにも直線探索を厳密に行っておく方が良い．なお，以上のアルゴリズムでは全て問題が線形であることを前提としているので，モデル計算も線形であることが必須条件である．線形の現象を扱っているモデルであっても，精度向上のために，しばしば非線形の処理がなされる場合があるが，そういった非線形の部分はあらかじめ外しておく必要がある．
5.

CO2 逆解析問題への適用

上記で紹介したアルゴリズムの有用性を確認するため，Niwa and Fujii（2020）では，大気二濃度の観測データから地表面フラックスを推定するという CO2 逆解析問題に
酸化炭素（CO2 ）
おいて，上記アルゴリズムを使って解析誤差共分散行列 Pa の導出を試みた．用いた数値モデを参照されたい．
ルや実験設定の詳細は Niwa and Fujii（2020）
図 1 は BFGS 公式で推定した Pa が反復回数を重ねるごとに解析解に近づいていく様子を示（a），（b）
では一つの BFGS 公式による計算結果を示し，図 1（c），（d）は 50 メ
している．図 1
ンバーのアンサンブルを使った場合を示す．また，ここで解析解との差をノルムと Ueno andで定義された divergence の 2 通りで示すことで，メトリックによる違いも確
Nakamura（2016）
認している．図 1 でわかることは，反復を重ねるとともにほぼ単調に推定誤差が減少すること，また，アンサンブルを用いることで，誤差の減少速度を大幅に短縮していることがわかる（上段と下段で横軸のスケールが異なることに注意）．このことから，BFGS 公式で用いる (y, p) の数を増やすことで，Pa の推定精度を高められていることがわかる．さらに，ここでは，共役性のある (y, p) を用いた場合に加え（実線），共役性のない場合（点線）も示している．共役性のない場合では，通常，よく用いられる厳密ではない直線探索を使っており，また，直交化も行（以後同じ．詳なっていないが，精度を高めるために BFGS 公式の計算を 10 回繰り返している．この両者の違いから，共役性を保つことで推定精度
細は Niwa and Fujii, 2020 を参照のこと）
が高められていることがわかる．
図 2 は Pa の非対角成分から計算した，東アジアのある場所のフラックスに対する誤差相関の分布を示す．ここでは，50 メンバーを用いて 15 回の反復を行った結果について，共役性が

204

統計数理

第 70 巻 第 2 号 2022

図 1．CO2 逆解析問題において，BFGS 公式で導出した解析誤差共分散行列 Pa とその解析解，また，Ueno and Nakamura（2016）
で定義された divergenceとの差のノルム（a，c）（b，d）の反復回数との関係．上段は一つの BFGS 公式の計算によるもの，下段は 50
メンバーを用いた場合を示す．また，実線は (y, p) ベクトルに共役性がある場合，点線はない場合を示す．Niwa and Fujii（2020）の Fig. 6 および Fig. 9 より改変．

図 2．Pa から求めた東アジアのある場所のフラックスに対する誤差相関の分布．共役性のあ．Niwa andる (y, p)（a）とない (y, p)（b）から推定したものと解析的に求めたもの（c）
Fujii（2020）の Fig. 8 より改変．

共役ベクトルと BFGS 公式を用いた解析誤差共分散行列の導出

205

保たれた場合と保たれていない場合の 2 通りを示している．一般的には対角成分よりも非対角成分の方が導出が困難であるが，共役性が保たれている場合には，解析解の分布をより正確に再現していることが見て取れ，本手法の精度の高さがうかがえる．なお，この時，直交化を経て実際の BFGS 公式に用いた有効な (y, p) ベクトルのペア数は 125 個であり，一方で，用いた観測の数は 260 個である．3 節で，前処理をした BFGS 公式では，共役性のある (y, p) ベクト
ルを観測の数だけ用意すれば，解析的に得たものと同等の Pa が得られることを示したが，この実験においては，およそ半分の数でほぼ十分な精度で推定できていることがわかる．このことは，多くの観測が互いに相関していることに起因していると考えられる．
なお，実際の解析では真の解析誤差共分散行列を知り得ないが，図 1 で示すように，2 つのメトリックは単調に減少しているため，その減少幅が十分に小さくなったところでの推定値が解析解とほぼ同等のものであると期待できる．なお，図 2 のケースでは，解析解と同等の分布を得るには 50 アンサンブルのケースで 15 回の反復回数が必要であったことから，より収束が遅く見えるノルムで判断した方が良いことがわかる．
6.

まとめと今後の課題

本稿では，Niwa and Fujii（2020）による解析誤差共分散行列 Pa の新たな推定手法について解説を行った．この手法は，データ同化または逆解析において 4 次元変分法を用いた際に，その最適解の誤差情報を正確に得ることを目的としている．Niwa and Fujii（2020）の手法では，まず，収束速度を速めるために，ヘッセ行列の近似逆行列の初期値に B を用いる．さらに，
BFGS 公式の計算で必要な (y, p) ベクトルの共役性を保つために厳密な直線探査を用いる．また，計算時間を短縮するためにアンサンブルで (y, p) の数を増やす場合には，得られた全ての(y, p) が共役性を保つように直交化を行う．Niwa and Fujii（2020）
では，以上の手法を CO2 逆解析の問題に適用して，実際に高い精度で Pa を推定できることを実証した．
また，Niwa and Fujii（2020）では，前処理を施した場合には，観測と同じ数の共役性のある(y, p) ベクトルを用意すれば解析値と同等の Pa が得られることを示したが，実際の実験ではその観測数よりもはるかに少ない数の (y, p) ベクトルで，解析的に得られるものとほぼ同等のPa が得られていた．しかし，それでも必要とされる (y, p) ベクトルの数は大きく，問題サイズが大きくなればなるほど，その数も大きくなるため，アルゴリズムのさらなる効率化が望まれる（本稿で示している実験では，比較のための解析解が容易に得られるように，実際の問題よりもサイズを落としている）．特に，直交化後の有効な (y, p) ベクトルの数 N  の N に対する割合は，N の数が大きければ大きいほど小さくなってしまうことが，一つの問題としてある．
図 3 は，反復回数に対する N  /N の推移を示しているが，回数を重ねるごとに有効なベクトルの取得率が低下していることがわかる．また，アンサンブル数が大きい方が，取得率が低いこ
とも図 3 で示されており，この取得率の面にだけ注目すれば，単一の反復計算によるものが効率が良い．したがって，本手法を実際の問題に適用する際には，利用可能なトータルの計算資
源と計算時間との兼ね合いでアンサンブル数と個々の反復数を決める必要がある．この N  /Nの減少について，一つの原因として考えられるのは，y や p の有効桁が反復計算を重ねるにつれて減少することである（アンサンブル数を増やすと N  /N がより減少するのは，直交化で yと p をかけている（式（4.7））
ためかもしれない）が，現在のところ確実な原因の解明には至っていない．今後は，その原因を突き止め，N  /N を大きくする手法を開発することが，さらなる効率化の 1 つのアプローチであると考えられる．

統計数理

206

第 70 巻 第 2 号 2022

図 3．反復回数を増やしていった場合の N  /N の推移．N はアンサンブルで得られた (y, p)の合計の数．N  は直交化後の有効な (y, p) の数．アンサンブルメンバー数が 10（灰（黒）の場合を示す．Niwa and Fujii（2020）の Fig. 10 より改変．
色）と 50謝辞査読者から有益なコメントを頂きました．感謝申し上げます．本研究は統計数理研究所共同
研究プログラム（28-共研-2004, 29-共研-2007, 30-共研-2007, 2019-ISMCRP-2030, 2020-ISMCRPの助成を受けたものです．
2049, 2021-ISMCRP-2023）

参

考

文

献

Bousserez, N., Henze, D. K., Perkins, A., Bowman, K. W., Lee, M., Liu, J., Deng, F. and Jones, D. B. A.
(2015). Improved analysis-error covariance matrix for high-dimensional variational inversions: Application to source estimation using a 3D atmospheric transport model, Quarterly Journal of theRoyal Meteorological Society, 141(690), 1906–1921, http://dx.doi.org/10.1002/qj.2495.
Chevallier, F., Bréon, F.-M. and Rayner, P. J. (2007). Contribution of the Orbiting Carbon Observatory to the estimation of CO2 sources and sinks: Theoretical study in avariational data assimilation framework, Journal of Geophysical Research, 112(D09307),
http://dx.doi.org/10.1029/2006JD007375.
Dennis, Jr., J. E. and Moré, J. J. (1977). Quasi-Newton methods, motivation and theory, SIAM Review,
19(1), 46–89, http://dx.doi.org/10.1137/1019005.
Derber, J. and Rosati, A. (1989). A global oceanic data assimilation system, Journal of Physical
Oceanography, 19, 1333–1347.
Fisher, M. and Courtier, P. (1995). Estimating the covariance matrices of analysis and forecast error
in variational data assimilation, Technical Memorandum, No.220, ECMWF, Reading, UK.
Fujii, Y. (2005). Preconditioned Optimizing Utility for Large-dimensional analyses (POpULar), Journalof Oceanography, 61(1), 167–181, http://dx.doi.org/10.1007/s10872-005-0029-z.
Fujii, Y. and Kamachi, M. (2003). A nonlinear preconditioned quasi-Newton method without inversion of a first-guess covariance matrix in variational analyses, Tellus A, 55(5), 450–454,
http://dx.doi.org/10.1034/j.1600-0870.2003.00030.x.
Golub, G. H. and van Loan, C. F. (1996). Matrix Computations, 3rd ed., Johns Hopkins University
Press, Baltimore, Maryland.
Ito, S., Nagao, H., Yamanaka, A., Tsukada, Y., Koyama, T., Kano, M. and Inoue, J. (2016). Dataassimilation for massive autonomous systems based on a second-order adjoint method, PhysicalReview E, 94, p.043307, http://dx.doi.org/10.1103/PhysRevE.94.043307.

共役ベクトルと BFGS 公式を用いた解析誤差共分散行列の導出

207

Lorenc, A. C. (1988). Optimal nonlinear objective analysis, Quarterly Journal of the Royal Meteorological Society, 114(479), 205–240, http://dx.doi.org/10.1002/qj.49711447911.
Nazareth, L. (1979). A Relationship between the BFGS and Conjugate Gradient Algorithms and
Its Implications for New Algorithms, SIAM Journal on Numerical Analysis, 16(5), 794–800,
http://dx.doi.org/10.1137/0716059.
Niwa, Y. and Fujii, Y. (2020). A conjugate BFGS method for accurate estimation of a posterior error covariance matrix in a linear inverse problem, Quarterly Journal of the Royal MeteorologicalSociety, 146(732), 3118–3143, http://dx.doi.org/10.1002/qj.3838.
Nocedal, J. and Wright, S. J. (2006). Numerical Optimization, 2nd ed, Springer, New York.
Rabier, F. and Courtier, P. (1992). Four-dimensional assimilation in the presence of baroclinic instability, Quarterly Journal of the Royal Meteorological Society, 118(506), 649–672,
http://dx.doi.org/10.1002/qj.49711850604.
Sasaki, Y. (1969). Proposed inclusion of time variation terms, observational and theoretical, in numerical variational objective analysis, Journal of the Meteorological Society of Japan, 47(2), 115–124,
http://dx.doi.org/10.2151/jmsj1965.47.2_115.
露木 義, 川畑拓矢 編 (2008).『気象学におけるデータ同化』, 気象研究ノート 217 号, 日本気象学会,
東京.
Ueno, G. and Nakamura, N. (2016). Bayesian estimation of the observation-error covariance matrix inensemble-based filters, Quarterly Journal of the Royal Meteorological Society, 142(698), 2055–
2080, http://dx.doi.org/10.1002/qj.2803.

208

Proceedings of the Institute of Statistical Mathematics Vol. 70, No. 2, 195–208 (2022)

Estimation of a Posterior Error Covariance Matrix Using Conjugate
Vectors and the BFGS Formula
Yosuke Niwa1,2 and Yosuke Fujii2,3
1 Earth System Division, National Institute for Environmental Studies

2 Department of Atmosphere, Ocean, and Earth System Modeling Research,

Meteorological Research Institute, Japan Meteorological Agency
3 The Institute of Statistical Mathematics

A four-dimensional variational method is commonly used for data assimilation/inverse
problems. However, it cannot automatically provide a posterior error. Niwa and Fujii(2020) developed a technique to estimate a posterior error covariance matrix within theframework of the four-dimensional variational method. Their technique adopts a quasiNewton method with the BFGS formula, which is a conventional optimizing method. To
enhance the estimation accuracy of a posterior error covariance matrix, their techniquealso employs an exact line search, an ensemble method, and orthogonalization to increase
the number of conjugate vectors used in the BFGS formula.
This report explains the fact that with preconditioning in the BFGS formula, ananalytical posterior error covariance matrix can be obtained from the same number of iterations (or vector pairs used in the BFGS formula) as observations followed by a detailed
description of Niwa and Fujii’s technique (2020). Finally, the results obtained by applying
this technique to an inverse problem of atmospheric CO2 are demonstrated.

Key words: Posterior error covariance matrix, data assimilation, inverse analysis, four-dimensional variational
method, BFGS formula, quasi-Newton method.