(q(t),p(t)) (q(t),p(t)) q = E[q], p = E[p], v 1 = E[(q q) 2 ], v 2 = E[(q q)(p p)], v 3 = E[(p p) 2 ]. ( t ). v 2 (t) q(t) p(t) E ξ(t), q(t), p(t) ()

解説：特集

カルマンフィルタを中心とした状態推定の理論から応用まで

連続時間カルマンフィルタと量子状態推定

山 本 直 樹

＊

＊慶應義塾大学理工学部物理情報工学科 神奈川県横浜市港北区日吉 3–14–1 ＊Department of Applied Physics and Physico-Informatics, Keio University,

3-14-1 Hiyoshi, Kohoku, Yokohama, Kanagawa, Japan ＊E-mail: yamamoto@appi.keio.ac.jp

キーワード：連続時間カルマンフィルタ (continuous-time Kalman filter), 量 子系 (quantum systems). JL 0009/17/5609–06622017 SICEC

1

.

はじめに

本稿の第一の趣旨は，カルマンフィルタの連続時間版 の紹介である．実は，フィルタの実装はデジタルで行う 事がほとんどで，かつ，離散時間のフィルタ理論はシン プルな数学展開が可能なので，実務的・教育的側面から 見ると離散系に関するガイダンスの方が重宝されそうで ある．しかし，物理系や，金融などの高頻度な不規則現 象が起こる系については，連続時間でモデリングするこ とが多い．この場合，離散近似せず，連続時間のままフィ ルタを作る方が自然である．種々のテクニカルな問題

(

サ ンプリング時刻の恣意性，離散近似系の発散など

)

が生じ ないからである．何より，微積分の確率版と言える確率 解析が使えることの恩恵は大きい．とは言っても，連続 時間のフィルタ理論は，数学的に厳密にやろうとすると 難しい．本稿では，そのような数学的難所を見つつかわ しつつ，考え方を伝えたい．第二の趣旨は，量子論でカ ルマンフィルタリングをすることができる，という事実 の紹介である．カルマンフィルタは工学では不動の地位 を築いているが，実は，物理とくに量子論においてはま だ新参者である．しかしすでにその有用性は大きなイン パクトをもって認められており，しかも，その概念が量 子論の思想にまで影響を与えつつある．カルマンフィル タの重要性は，工学の枠を超えてますます高まっている のである．そういう話を，ざっくばらんに書いてみたい．

2

.

簡単な例

—

調和振動子の動的推定

2.1

モデリングの仕方 まず，百聞は一見に如かずということで，具体例から 出発しよう．基本中の基本素子，調和振動子である．た とえば，図

1 (a)

に示すような，バネにつながっていて

1

次元的に振動する質点をイメージするとよい．質点の重 さを

M ,

周波数を

ω

とする．これがさらに大きさ

γ > 0

図1 (a)ふつうの調和振動子 (b)量子調和振動子 の抵抗力とランダムな力

Ξ(t)

を受けるとき，質点の位 置

Q(t)

はつぎの運動方程式に従う：

M

d

2

_Q(t)

dt

2

+ Mγ

dQ(t)

dt

+ Mω

2

_{Q(t) = Ξ(t)}

ここで運動量

P = M dQ/dt

を導入し，さらにメートル の次元をもつ量

a > 0

によって

Q = aq, P = aM ωp

と 無次元化した量

(q, p)

を考えると，これは

d

dt



q

p



=



0

ω

−ω −γ

 

q

p



+



0

1



ξ ( 1 )

なる

1

階の微分方程式に従う．

ω

と

γ

は，頑張れば正確 な数値を実験的に特定できる．つまり，既知の量として よい．ノイズ

ξ(t) = Ξ(t)/aM ω

のモデル化は難しい． よくやるのは，スペクトルが全周波数で一様である白色 ノイズを仮定する（注1），というものである．これは数学 的には

E[ξ(t)] = 0, E[ξ(t)ξ(s)] = ¯nδ(t − s)

( 2 )

を満たす確率変数である．

ξ(t)

はガウス分布に従ってお り，

E

はこの分布に関する平均を表わす

(Expectation).

確率論が必要になるのか，と嫌気がさすかもしれないが 深く考えなくても良い．要は，大きさ

¯n

で，振動子よ りはるかに速く変化するノイズが加わる，ということで ある

(

たとえば，質点が

1

往復する間に

1000

回質点を つっつく

).

一応追記しておくと，

δ(t − s)

は

t = s

のと きゼロで

t = s

のとき無限大をとるデルタ関数であり，

(2)

式は「

ξ(t)

は任意の異なる

2

時点で独立である」こ とを意味する．ついでに言うと，デルタ関数は直観的に は便利だが数学的取り扱いには注意が必要で，そこで確 率微分方程式に基づく精緻な理論が必要云々となるのだ が，これが，連続時間系のフィルタ理論が離散系のそれ と比べて難しくなってしまう原因である． さて，無事に

(?)

振動子のモデリングができた．ここ で，この振動子が複雑な機械系

(

自動車のエンジン等

)

の 一部をなしていて，直接モニターできないとする．では 問題．変数

(q(t), p(t))

はどんな値をとっているか？しか （注1）_{もちろん，現実のノイズはある周波数帯だけで効いてくる「有色ノ} イズ」である．この帯域を無限に延ばすのである．なお，白色ノイ ズをローパスフィルタに入れて有色ノイズをモデル化してもよい．

しいま

(q(t), p(t))

も確率変数なので，これらは「この くらいの値をとるだろう」としか答えようがない．そこ で，上で確率については深く考えなくて良いと言ったば かりで恐縮だが，

(q(t), p(t))

の平均と分散：

¯q = E[q], ¯p = E[p], ¯v

1

= E[(q − ¯q)

2

],

¯v

2

= E[(q − ¯q)(p − ¯p)], ¯v

3

= E[(p − ¯p)

2

].

を考える

(

時刻

t

は省略した

). ¯v

2

(t)

は，

q(t)

と

p(t)

の 相関を表わす共分散である．また，ここでの

E

は，

ξ(t),

q(t), p(t)

が従う確率分布

(

つまり同時分布

)

についての 平均である．これら統計量の動きがわかれば，変数のお およその振る舞いは把握できる．実際，

(2)

式から

d¯

q

dt

= ω¯p,

d¯

p

dt

= −ω¯q− γ ¯p,

d¯

v1

dt

= 2ω¯v

2,

d¯

v2

dt

= ω(¯v

3

− ¯v1

) − γ¯v

2,

d¯

v3

dt

= −2ω¯v

2

− 2γ¯v3

+ ¯n

( 3 )

が出てくる．導出法については後述する．これは連立の

1

階線形微分方程式なので解けて，そして時間を大きく 取ると，変数はつぎの値に収束する：

¯q(∞) = ¯p(∞) = 0, ¯v

1

(∞) = ¯v

3

(∞) =

_2γ

¯n

,

¯v

2

(∞) = 0.

つまり，

(q(t), p(t))

は原点まわりを大きさ

¯n/2γ

でふ らついている（注2）．これより，抵抗力

γ

が小さくなる と，エネルギーが散逸しにくくなるので質点のふらつき が大きくなることがわかる．とくに

γ → 0

とすると，

¯v

1

(∞), ¯v

3

(∞) → ∞

となる．パチンコ玉のような小抵 抗の物体をランダムにつっつきまわすと，その物体がど こにあるのかさっぱりわからなくなる，ということであ る

(

バネ付き質点の場合，実際にはそうはならないが

).

2.2

フィルタリングの方法 さて，この機械系を精密に制御したいとする．どうす ればよいか？ひとつの方法として，金銭的・技術的ハー ドルは上がるが，質点の付近に

CCD

カメラ等を取り付 けて，質点の位置をモニターし続ければよい

(

図

1 (a)

参 照

)

．そうすれば，質点の運動量についての情報も得られ るかもしれないし，ひいては，その情報を使って質点を 自由自在にフィードバック制御できるかもしれない（注3）_． そのような，「質点の位置がモニターできている」という 状況をつぎのようにモデリングする．すなわち， （注2）_{¯v1+ ¯v3 は質点のエネルギーを表わす．なのでこの式は，系の定常} 状態におけるエネルギーは，系を揺らす力の大きさ ¯n と系からエネ ルギーを散逸させる力の大きさ γ の釣り合いで決まる，という「揺 動散逸定理」の一形態を表わす． （注3）_{質点を原点付近にとどまらせたいなら，つぎのようなフィードバッ} ク制御を行えばよい：質点の運動量が正の値をとるなら負方向に力 を加え，質点の運動量が負の値をとるなら正方向に力を加える．

y(t) =

√

λq(t) + ζ(t)

( 4 )

なる時系列データ

y(t)

がリアルタイムで入手できてい る，とする．

ζ(t)

はカメラの精度を反映する測定ノイズ で，

ξ(t)

とは独立で，大きさ

1

の白色ノイズとする

(

ゆ えに

E[ζ(t)] = 0, E[ζ(t)ζ(s)] = δ(t − s)

を満たす

). λ

は信号

q(t)

とノイズ

ζ(t)

の比を表わすパラメータ

(

信 号対ノイズ比：

S/N

比

)

で，これが大きいほど位置が上 手くモニターできている． さて，問題を記述しよう．いま，振動子の情報を含んで いるデータ

y(t)

がリアルタイムで得られている．やりた いことは，これを用いて質点の変数

(q(t), p(t))

を「推定 する」ことである．具体的には，

y(t)

を用いて，全時刻

t

で，

(q(t), p(t))

に「近い」確率変数

(ˆq(t), ˆp(t))

をつく る．どうやってそのような変数をつくればよいか？すぐ 思いつく案として，データ

y(t)

を

S/N

比で割った量は ほぼ

q(t)

だろうということで

ˆq(t) = y(t)/

√

λ

とする． では

p(t)

ˆ

をどうするかと言うと，これはモデルを通して 動的に推定する．つまり，

(1)

式から

(q(t), p(t))

の関係 はわかっているのだから，思い切ってノイズ

ξ(t)

を無視 して，推定値

(ˆq(t), ˆp(t))

を

dˆ

p(t)/dt = −ω ˆ

q(t) − γ ˆ

p(t)

で関係付ける．これは，

dˆ

p(t)

dt

= −γ ˆp(t) −

ω

√

λ

y(t)

( 5 )

と表わせる．

y(t)

は時々刻々手に入るのだから，以上の アルゴリズムを用いて

(ˆq(t), ˆp(t))

をリアルタイムで更 新していけばよい．このやり方は実はさほど良い方法で はないのだが，しかしともあれ「モデルベースの，動的 な推定器をつくる」という方針が得られた．こう言うと 簡単に聞こえるが，この方針こそが，カルマンが発見し た推定理論におけるパラダイムシフトだったのである． ではどうやって動的推定器をつくれば良いのか．望ま しい推定値は，変数

(q(t), p(t))

に最も近いもの，つまり 各時刻でつぎの推定誤差を最小とするものが良いだろう：

v1

(t) = E[(q(t) − ˆq(t))

2

],

v3

(t) = E[(p(t) − ˆp(t))

2

].

ここでの

E

は，確率変数

ξ, ζ, q, p

が従う同時確率 分布についての平均を表わす．問題は，この最適推定値

(ˆq(t), ˆp(t))

を時々刻々更新する式

((5)

式のような，

y(t)

で駆動されるダイナミクス

)

を求めることである．導出 は相当に非自明ではあるが，答えはつぎである：

dˆ

q

dt

= −λv

1

ˆq+ ωˆp + λv

1y,

dˆ

p

dt

= −(ω + λv

2

)ˆq − γ ˆp + λv

2y,

dv1

dt

= 2ωv

2

− λv

2 1

,

dv2

dt

= ω(v

3

− v1

) − γv

2

− λv1

v2,

dv3

dt

= −2ωv

2

− 2γv3

+ ¯n − λv

2 2

( 6 )

これがカルマンフィルタである．ノイジーなデータ

y

か ら，有用な情報である

(ˆq, ˆp)

を「濾過」する働きをもつ． なお，

v2

= E[(q − ˆq)(p − ˆp)]

は

(ˆq, ˆp)

間の相関を表 わす．

(6)

式のうち，推定誤差

(v

1, v2, v3

)

に関する式はリ カッチ方程式と呼ばれており，非線形方程式ではあるが 多くの性質が判明している．まず重要なことに，リカッチ 方程式はデータ

y(t)

に依存しない．データと関係なく， これくらいの精度で推定できますよ，と理論予測できる 訳である．また，いま推定誤差が収束する条件ははっき りしており，これは

γ = 0

というケースも含んでいる． このとき前述のとおり，振動子はノイズでつっつかれて どこか遠くに吹っ飛んでしまうこともあるのだが，カル マンフィルタはそれをキチンと追いかけてくれるのであ る！なお

(v

1, v3

)

の式を見ると，これは平均ダイナミクス

(3)

の

(¯v

1, ¯

v3

)

に関する式の右辺から

λv

2 1

, λv

22 という非 負量を引いたものであることがわかる．つまり

(v

1, v3

)

の増加率は

(¯v

1, ¯

v3

)

のそれと比べて小さいので，結果，

v1

≤ ¯v1

, v

3

≤ ¯v3

となる．これは，データ

y(t)

を用い ている分，振動子に関する不確かさが減少したことを意 味する．なお，この不確かさの減少度合いを表わすパラ メータ

(S/N

比

)

を

λ = 0

とする，つまりデータが全く 得られないとすると，カルマンフィルタ

(6)

は平均ダイ ナミクス

(3)

と一致する．これは自然である．

2.3

数値例 カルマンフィルタがどれほど上手くシステムの変数を 追跡できるのか，数値例で確認しておこう（注4）．パラメー タは

ω = 1, λ = 1, γ = 0.01, ¯

n = 1

とする．図

2

は，振 動子の変数

(q, p),

平均値

(¯q, ¯p),

推定値

(ˆq, ˆp)

の時間変 化を示したものである．初期値はいずれの場合も

(2, 0)

で，また初期分散値は

(¯v

1, ¯

v2, ¯

v3

), (v

1, v2, v3

)

のどちら も

(0.1, 0, 0.1)

とする．

(¯q, ¯p)

と

(ˆq, ˆp)

のまわりのシェー ドは，誤差の標準偏差

(√¯v

1,

√

¯v

3

)

と

(√v

1

,

√

v3

)

を表 わし，つまりたとえば

¯q

を太さ

2√¯v

1 の線としてプロッ トしている．図からわかるとおり，推定値

(ˆq, ˆp)

は真値

(q, p)

を良く追跡できている．一方，平均値

(¯q, ¯p)

は真 値の動きの傾向は捉えているものの，追跡とは言い難い． また，

(ˆq, ˆp)

まわりのシェードの方が，

(¯q, ¯p)

のそれより 小さいことがわかる．前述したとおり，データを用いる事 で対象系に関する不確かさが減るのである．とくに真値

(q, p)

はほぼ

(ˆq, ˆp)

まわりのシェードに入っており，つ まり理論予測の誤差範囲内に真値がある，ということも 確認できる．なお，データ

y(t)

は

q(t)

を陽に含んでい （注4）_{Matlabで実装した．数値計算の仕方については文献}1)_{を参照され} たい． 図2 振動子の変数(q, p), 平均値(¯q, ¯p), 推定値(ˆq, ˆp)の 時間変化 q(t) は確率揺らぎを直接的には受けていないので，その動きは滑らかに見える． るが，

p(t)

を陽には含んでいない．しかしそれでも

p(t)

が推定できているのは，ひとえにカルマンフィルタがモ デルベースの動的推定機構であるからである．すなわち，

p(t)

は，ダイナミクスのモデルを通して推定しているの である． 最 後 に ，重 要 な 点 を 注 意 し た い ．第 一 に ，平 均 値

(¯q(t), ¯p(t))

はリアルタイムで得られる量ではなく，実 験をたくさん行って，その後で

(

ヒストグラムを作る等 して

)

算出するオフラインの予測値である．当然，これを 使ってフィードバック制御などはできない．他方，推定 値

(ˆq(t), ˆp(t))

は，データ

y(t)

を用いてリアルタイムで 算出するオンラインの予測値である．ゆえに「いま，こ の瞬間に振動子はどこにいて，どのくらいの速さで動い ているか」について検討がつけられるので，したがって フィードバック制御ができるのである．

3

.

一般論らしきこと

連続時間カルマンフィルタの導出は，数学的に真面目 にやろうとすると難しい．そこで，ここでは直観的な方法 でそれを示そうと思う．連続時間系における数学的難所 を知ることにも意義があるだろう．難所の解決には確率 微分方程式が必要である．たとえば1), 2)_{を参照されたい．} つぎのシステムを考えよう：

dx(t)

dt

= Ax(t) + Bu(t) + w(t),

y(t) = Cx(t) + v(t)

( 7 )

x, y, w, v

は確率変数を成分とするベクトルで，

A, B, C

は行列である．また，前章の例では導入していなかった が，

u

は制御入力である．フィードバック制御をしたけ れば，

u

を

y

の関数として設計する．

w, v

は前章同様，

平均ゼロの白色ノイズのベクトルであり，つぎを満たす：

E[w(t)w(s)



_{] = Nδ(t − s),}

E[v(t)v(s)



_{] = Mδ(t − s),}

E[w(t)v(s)



] = Lδ(t − s)

( 8 )

N, M, L

は適当なサイズのある行列である（注5）．

L = 0

のとき，

w(t), v(t)

に相関があることになる．前章の調 和振動子系が，この一般系に含まれている事を確認され たい． さて，目標はデータ

y(t)

を用いて

x(t)

に最も近い推 定値ベクトル

ˆx(t)

を獲得することである．前章で明らか になったように，これを，モデルベースで，動的に実行す る．とくに，

ˆx(t)

の更新則としてつぎを採用する（注6）：

dˆ

x(t)

dt

= Aˆx(t) + Bu(t) + K(t)[y(t) − C ˆx(t)]

( 9 )

u(t)

は観測者が設計する量なので，

(7)

，

(9)

式で共通に 取れる事に注意しよう．

(9)

式に現れる

K(t)

を上手く 選びたい．そこで，誤差ベクトル

e(t) = x(t) − ˆ

x(t)

を 考え，これの「大きさ」を小さくすることを考える．具 体的には，共分散行列

V (t) = E[e(t)e(t)



]

を小さくす る．「行列を小さくする」ことの意味が曖昧だが，気にせ ず先に進もう．まず，誤差ベクトルは

de(t)

dt

= (A − K(t)C)e(t) + w − K(t)v

に従う（注7）．これに基づいて

V (t)

の変化則も求める． そのために，

e(t)

のテーラー展開：

e(t + Δt) ≈ e(t) +

de(t)

dt

Δt

を思い出そう．

Δt

は時間の微小量である．これより，

d[e(t)e(t)



]

dt

≈

e(t + Δt)e(t + Δt)



− e(t)e(t)



Δt

= de

dt

e



_{+ e}



de

dt

_

+ de

dt



_de

dt

_

Δt

と計算できる．

Δt → 0

とやれば，普通の積の微分則に 見える．ただ，いま

de/dt

に

(

デルタ関数で特徴付けら （注5）_{具体的には，N, M は正定対称行列である．} （注6）_{dˆx/dt = F ˆx + Bu + Ky から出発しよう．これに「不偏性」の条件} E[x(t)] = E[ˆx(t)] ∀t を課すと，この式が出てくる． （注7）_{推定対象が (7) 式で w = 0, v = 0 とした確定系の場合，e → 0 と} なれば良いので，基本的には K は A − KC が安定であれば何で もよい．より発展的には，A − KC が望ましい極をもつように K を選ぶ．((A, C) 可観測であれば任意の極配置が可能である．)この ような推定器をオブザーバという．つまり，オブザーバは統計的仮 定を置かず極配置により K を決める．一方，カルマンフィルタは ノイズに統計的仮定を置き，推定誤差を最小にするという規範で K を決める． れる

)

瞬間パルス的な入力

w, v

が含まれているので，こ こは注意が必要である．雑で恐縮だが，

(8)

式を

E[w(t)w(t)



] = N/Δt, E[v(t)v(t)



] = M/Δt,

E[w(t)v(t)



] = L/Δt

と読み替える．

Δt → 0

とすると，たしかに右辺はすべ て無限大になる．これにより，

V (t)

の時間変化がつぎの ように計算できる：

dV (t)

dt

= E



_d[e(t)e(t)



_]

dt



= E



(A − KC)ee



+ (w − Kv)e



+ ee



_{(A − KC)}



_{+ e(w − Kv)}



+ (w − Kv)(w − Kv)



Δt



= (A − KC)V + V (A − KC)



+ N − KL



_{− LK}



_{+ KMK}

 普通の積の微分則を使うと，最後の行が出てこないので ある．この辺りが連続時間の確率解析ならではの結果で， 同時に，難しいところである．さて，上式はさらに

dV

dt

= AV + V A



_{+ N}

− (V C



_{+ L)M}

−1

_{(V C}



_{+ L)}



+ [K − (V C



+ L)M

−1

]M[K

− (V C



+ L)M

−1

]

 と変形できる．決定変数

K(t)

は最後の項だけに含まれ ている．一方で，

M

が正定

(

ゆえに

M

−1 も正定

)

であ ることから，この最後の項は非負である．なので，一般 の

K(t)

についてはこの項は

dV (t)/dt

が増加すること に寄与してしまうが，唯一，

K(t) = (V (t)C



+ L)M

−1 のときに限って，この最後の項がゼロになる．つまり，こ れが

V (t)

の増加を最も抑える，すなわち

V (t)

を最も小 さくする

K(t)

の選び方である．このようにして，

ˆx(t)

の更新則と誤差共分散行列の変化則が求まった：

dˆ

x

dt

= Aˆx + Bu + (V C



_{+ L)M}

−1

_{(y − C ˆx)}

dV

dt

= AV + V A



_{+ N}

− (V C



_{+ L)M}

−1

_{(V C}



_{+ L)}

 これがカルマンフィルタである．行列を前章のものに合 わせると

(6)

式が出てくる．また，この式から，情報が 全く得られないときの平均ダイナミクスの式も得られる． すなわち，

C = 0, L = 0

とすれば

d¯

x

dt

= A¯x + Bu,

d ¯

V

dt

= A ¯

V + ¯

V A



_{+ N}

(10)

となる．

(3)

式はここから出てくる．この場合，

u(t)

は

y(t)

の関数ではなく，事前設計する制御入力

(

オープン ループ制御

)

であることに注意されたい．

4

.

量子カルマンフィルタ

4.1

雑談

30

∼

40

年くらい前まで，スタンダードな量子力学は 「見ない」理論だった．先に，調和振動子が複雑な機械系 に組み込まれていて直接モニターできないという状況を 考えた．量子ではほとんどこういったことになっていて， 典型的には，原子集団や超伝導回路が冷却器などの中に 閉じ込められていて，差し当たり「見ない」のである．こ のとき，前述のとおり，平均過程を考えるしかない．ふ つうの調和振動子が相手の場合，

(3)

式を調べる．量子 でも調和振動子はきわめて重要なシステムで，実際，

(3)

式とほぼ同じ方程式が解析の対象となる．これはエーレ ンフェスト方程式と呼ばれる．また，平均過程を考えて いるので，当然，背後に確率分布が存在する．いわゆる シュレーディンガー方程式とは，この確率分布

(

らしき もの

)

の時間発展則のことである．さて，何かしらの実 験結果を得るためには，やはり「見る」必要がある．そ こで，ある時刻

(

とくに定常になっているような時刻

)

で 系を測定する．たとえば冷却器を開けて，漏れ出る発光 を調べる．この測定結果は，たとえば原子がたくさんあ るなどの理由で，ほぼ理論で予測できる平均値と一致す る

(

大数の法則！

)

．要するに，従来は，量子力学は平均 過程の力学だった．「測定」は統計量を算出するための道 具という以上の意味をもたず，またもたせようとした研 究は「変わり種研究」という不幸な扱われ方をされる傾 向もあった． しかし，量子情報科学が状況を一変させた．量子系に 対して巧妙な測定を行うと，古典ではあり得ない現象が 起こることが判明したのである．その極端なものが量子 計算である．この事情で，測定を積極的に行うあるいは デザインする研究が，突如，変わり種研究から王道研究 に格上げされた．そうなると，「時間連続で測定し続ける」 ための方法論も欲しくなる．そして実際，そのような理 論が急速に整備され，連続測定された量子系の時間発展 方程式がある直観的な方法で導かれた．この時間発展則 は，平均過程のそれとは全く異なる．つまり，エーレン フェスト方程式ないしシュレーディンガー方程式ではな い．この式は，確率的に変動する，リアルタイムで入手で きる時系列データを入力する方程式になっている．導出 が直観的ゆえその物理的意味付けが不明であったが，と もあれ量子軌道方程式，あるいは確率シュレーディンガー 方程式と名前が付けられた．しかしあるとき，物理学者 たちは，これがカルマンフィルタと同じ形式になってい ることに気付いた．フィルタとは何だ？カルマンとは何 者だ？また，過去の「変わり種研究」の中に，量子連続 測定とフィルタの関係を調べたものがあった．ベラフキ ンの研究である．彼がやったことは何なのか？ ここで，話が筆者の個人的な回想に変わる．ベラフキ ンの理論は重要そうではあるが，数学的に難解すぎ，誰 もフォローできない代物であった．

2004

年に学位をとり ポスドクとして渡米したばかりの私も，この理論の上辺 だけを使って研究をしていた．本質的理解を伴わない研 究を行っているとフラストレーションがたまる．しかし いつの世も勇者はいるもので，彼の名はボウテンという． 彼はその才能をもってベラフキン論文の解読を行い，そ の本質を見抜いた．ベラフキンがやったことは，まさに カルマンフィルタそのものだったのである．つまり，量 子軌道方程式とか確率シュレーディンガー方程式と呼ば れていた方程式は，量子力学的変数を

(

最小

2

乗平均の 意味で

)

最もよく近似する確率変数の時間発展則にほか ならなかったのである．こう書くと簡単そうに聞こえる が，量子力学は変数が交換しない云々という古典とは劇 的に違う体系で記述されるので，量子のフィルタリング は全くもって非自明である．この非自明な理論について， ボウテンは皆が理解できるような懇切丁寧な論文5)_を書 いた．私は幸運だった．ボウテンは私と同年齢で，同じ 大学・同じ研究室に同じ年にポスドクとして渡米し，そ して彼は私に量子フィルターのレクチャーをマンツーマ ンで行ってくれたのである．すべてが氷解した．フラス トレーションの解消どころではない．恥ずかしながら当 時私は通常の

(

量子ではない

)

フィルタ理論も良くわかっ ていなかったが，それの本質も同時に理解できた．そし てカルマンフィルタがいかに強力で革新的な方法論なの かを知った．今に至るまでの私の研究の多くは，この米 国での経験に基づいている．大げさに言うと，カルマン フィルタとその量子版が，私の人生を決めたのである．

4.2

量子力学的調和振動子の動的推定 調子に乗って書いているうちに紙数が尽きかけてきた． 看板に偽り有りと言われないように，最後に，簡単に量 子カルマンフィルタの事例を紹介したい．一般論に興味 がある方は3), 4)_{を参照されたい．} 対象は図

1 (b)

で示す調和振動子を極端に小さくした ものである．質点はナノグラムスケールで，

1

次元方向 の移動距離はフェムトメートルのオーダである．ここま で小さいと，

CCD

カメラ等の分解能では質点の位置を モニターすることはできない．そこで，振動方向にレー ザを照射して，跳ね返ってきた光を測定する．

(

図で示す ように，しばしば，部分透過鏡を置いて光を往復・増幅 して感度を上げることをする．

)

光の変化は時間連続で検 出する事ができるので，これで，質点の動きを時間連続

測定するわけである．このシステムのダイナミクスは次 式で与えられる：

d

dt



q

p



=



0

ω

−ω −γ

 

q

p



+



0

√

λζ1

+

√

2γξ



y =

√

λq + ζ2

これは古典の調和振動子のダイナミクスと良く似ている． どこが量子なのか？実は，変数が

q(t)p(t) − p(t)q(t) = i

という関係式をつねに満たすのである．

(q, p) = (3, −1)

などとするとこの関係式は成立せず，なのでこれは行列 のような「交換しない」量である．ほかにも相違点があ る．ダイナミクスのノイズ項に，

√

λζ1

という，測定の

S/N

比を係数とするノイズ項が加わっている．これは，

S/N

比を上げるとダイナミクスに加わるノイズも大きく なることを意味し，この意味で反作用ノイズと呼ばれる． 一方で

ζ2

はショットノイズと呼ばれ，これら

(ζ

1, ζ2

)

は， 上述した，測定用の入射レーザ光の揺らぎを表わす．重要 なことに，これらは同時に小さくする事はできない．不 確定性関係

E(ζ

₁2

)E(ζ

₂2

) ≥ ¯h

2

/4

を満たすのである！

(E

は量子的な意味での平均である．以下同様．

)

また，一般 には

L = 0

である．最後に，

ξ

は質点を構成する分子 の熱的揺らぎを表わし，

(2)

式と同じ統計的性質をもつ．

¯n

は熱揺らぎの大きさを表わす．また，散逸の大きさ

γ

が，この熱ノイズに起因することが見て取れる． さて，

y(t)

は跳ね返ってきた光を時間連続的に測定し て得られる時系列データである．これは，式からわかる とおり質点の

“

位置

” q(t)

に関する情報を含んでいる．な ので，

y(t)

を使って

(q, p)

を時間連続的に推定すること ができる．これが，量子カルマンフィルタである．上述し たように

(q, p)

はいま「行列のようなもの」なので，こ れを推定することの意味が曖昧であるが，ともあれ古典 と同じく，

E[(q − ˆq)

2

]

と

E[(p− ˆp)

2

]

を最小にする

(ˆq, ˆp)

を求めることができる．なお，この

(ˆq, ˆp)

は普通のスカ ラー値確率変数である．これらが従う式，すなわち推定 値の更新則は

(6)

式と似たものになる．文献6)_{は，この量} 子カルマンフィルタの実験的実現とその性能について報 告したものである．この実験では

M = 1.76 × 10

−10

kg

という極小の質点を用いて，

ω = 8.03 MHz

の調和振動 を実現している．とてつもなく速い振動である．ほかの パラメータは

γ = 1.66 kHz, λ = 467 Hz, ¯

n = 2 × 10

5 と同定されている

(

文献における弱結合領域の場合

).

こ の設定で，まず，測定値を使わずに算出した平均値

(¯q, ¯p)

の不確かさの幅

(

おおよそ標準偏差

)

が，いずれも約

800

であることが示されている．そして，量子カルマンフィ （注8）_{無次元化する前の位置 Q = aq に換算すると，これは約 43.7 フェ} ムトメートルの推定誤差に相当する．なお，Q と P = aMωp は QP − P Q = i¯h を満たし，これより換算定数が a =¯h/M ω = 2.73 × 10−16mと決まる．¯h = 1.055 × 10−34はプランク定数． ルタが算出する推定値

(¯q, ¯p)

が，およそ

150

の不確かさ で変化していることが示されている（注8）_{．すなわち，カ} ルマンフィルタによって，不確かさが約

1/5

に減少する． なお，この

150

という大きさは熱ノイズの大きさ

¯n

に 起因するもので，これがなければ

1

のオーダまで下がる． つまり，実験ではこの零点振動と呼ばれる量の約

150

倍 の不確かさで推定ができているわけである． 以上のように，量子系の連続モニタリングとカルマン フィルタリングが，いまや実験的に可能である．こうな ると，量子フィードバック制御もできるようになる．と くに推定値をフィードバックして振動子を原点付近にと どまらせようとする方法が「フィードバック冷却」であ り，最近，この方法の実験的実装も報告されている

(

た とえば7)

)

_{．このような状況は，量子力学が「見ない理論」} であった時代では夢物語である．なお，その夢物語では 「デーモン」なるものが登場する．ミクロな粒子ひとつひ とつを「見て」，その情報を使って速い粒子と遅い粒子 を選別する．遅い粒子だけからなる集団は，もとの集団 より「冷えている」．なので上述のフィードバック冷却 は，デーモンの実現と見ることもできる．夢の実現の鍵 はカルマンフィルターであったわけである． （2017 年 6 月 1 日受付） 参 考 文 献 1） 大住 晃：確率システム入門，朝倉書店 (2002) 2） 山本直樹：確率微分方程式の基礎，計測と制御，50-11, 937/943 (2011) 3） 山本直樹：量子フィードバック制御の数理，数理科学，585, 21/27 (2012)

4） H.I. Nurdin and N. Yamamoto: Linear Dynamical Quantum Systems: Analysis, Synthesis, and Control, Springer (2017) 5） L. Bouten, R. van Handel, and M.R. James: A discrete

invita-tion to quantum filtering and feedback control,SIAM Review,

51, 239/316 (2009)

6） W. Wieczorek, et al.: Optimal state estimation for cavity op-tomechanical systems,Phys. Rev. Lett., 114, 223601 (2015) 7） D.J. Wilson, et al.: Measurement-based control of a

mechan-ical oscillator at its thermal decoherence rate, Nature, 524, 325 (2015) ［著 者 紹 介］ やま 山 もと本 なお直 樹 君き (正会員) 1999年東京大学工学部計数工学科卒業，2004 年 同大学大学院情報理工学系研究科博士課程修了．カ リフォルニア工科大学研究員，オーストラリア国立 大学研究員を経て 2008 年，慶應義塾大学理工学部 物理情報工学科専任講師．2011 年，同大学准教授． 量子制御理論，量子情報理論，量子光学の研究に従 事．博士 (情報理工学)．