RENSx アルゴリズム解説
RENSx algorithm: a versatile source localization algorithm forvarious types of biomagnetic imaging
株式会社

シグナルアナリシス

関原謙介

はじめに

1

RENSx アルゴリズムは脳磁界計測を始めとして，脊髄磁場計測，末梢神経磁場計測あるいは心磁界計測など，生体磁場計測のいろいろな分野において有用な信号源推定アルゴリズムである．RENSx
アルゴリズムはスパースベイズラーニング（sparse Bayes learning, SBL）をベースにしたアルゴリズ
ムであるが，生体磁気計測データからの信号源再構成問題に対応するため，SBL を土台としていくつかの工夫・改良が施されている．一般的な SBL アルゴリズムについては，ベイズ信号処理 [1] の第 5 章
やシグナルアナリシス社技術資料 [2] に詳しい解説があるので，このノートは，このアルゴリズムの特に信号源再構成問題に特化した部分を解説する．このノートにおける説明は，ほぼ Electromagnetic

Brain Imaging (EBI)[3] の第 4 章の内容に沿ったものである．

2

種々の定義とデータモデル

2.1

データベクトルとソースベクトル

m 番目のセンサーの時刻 t における出力を ym (t) で表し，センサーアレイ全体の出力を列ベクトル
y(t) = y1 (t)y2 (t)
..
.








(1)

yM (t)で表す．ここで，M はセンサーアレイに含まれるセンサー数である．3 次元空間における位置をベ
クトル r = (x, y, z) で表す．位置 r に計測の対象となっている信号を発している信号源が存在すると仮定する．信号源は時間的に変動する 3 次元ベクトルであり，位置 r における時刻 t の信号源（ソース）ベクトルを s(r) と表し，


sx (r, t)s(r, t) =  sy (r, t) sz (r, t)とする．ここで sx (r, t), sy (r, t), sz (r, t) はソースベクトルの x, y, z 成分である．

1

(2)

2.2

センサーリードフィールド（センサー感度ベクトル）

次にセンサーリードフィールドを定義する．位置 r に単位強度の信号源が存在すると仮定し，この
x信号源が，x 方向を向いている場合のこの単位強度の信号源による m 番目のセンサーの出力を lm(r)yとする．同様に，この信号源が y 方向を向いている場合の出力を lm(r), z 方向を向いている場合のz
出力を lm(r) とする．すると，3 次元の行ベクトル lm (r),xyzlm (r) = [lm(r), lm(r), lm(r)](3)

が m 番目のセンサーの位置 r における感度を表す．この lm (r) を用いて，M × 3 の行列を
l1 (r)
..

L(r) = .
lM (r)
(4)と定義すれば，この L(r) がセンサーアレイ全体の位置 r における感度を表す．この行列 L(r) はリー
ドフィールド行列と呼ばれる．このリードフィールド行列を用いて，観測ベクトル y と空間に分布した信号源ベクトル s(r) の関係はZy=

L(r)s(r) dr(5)Ωと表現できる． ここで，右辺の積分は信号源が存在する可能性のある 3 次元の領域 Ω について行う．
この Ω を信号源空間 (source space) と呼ぶ．

2.3

離散ソースモデル

式 (5) に示される積分方程式を，信号源空間をボクセル分割により離散化し，線型方程式に変換する．各ボクセルの座標を r 1 , r 2 , . . . , r N とする．ここでボクセル上に存在する離散的な信号源で信号源分布が近似できるとして，信号源分布をデルタ関数を用いてs(r, t) =NXs(r j , t)δ(r − r j )j=1と表せば，これを式 (5) に代入して

y(t) =N
X

L(r j )s(r j , t)

(6)

j=1

を得る．ここで以下に定義される M × 3N 行列

F = [L(r 1 ), L(r 2 ), . . . , L(r N )]2(7)

を導入する．この F はボクセルリードフィールド行列と呼ばれる．さらに，3N × 1 の列ベクトルx(t) をs(r 1 , t)
 s(r 2 , t) 
x(t) = 
..



.

(8)

s(r N , t)と定義する．この x(t) はすべてのボクセル位置の信号源ベクトルを含む列ベクトルであり，式 (2) で定義されたソースベクトルと区別するためボクセルソースベクトルあるいはボクセルベクトルと呼ぶ．すると，式 (6) はs(r 1 , t)
 s(r 2 , t) 
 = F x(t)y(t) = [L(r 1 ), L(r 2 ), . . . , L(r N )] 
..



.

(9)

s(r N , t)と表すことができる．さらに，観測結果には，信号成分に加法的にノイズ ε が重畳すると仮定すれ
ば，観測データのモデルはy(t) = F x(t) + ε(10)となる．ノイズベクトル ε は M 次元の列ベクトルで ε = [ε1 , ε2 , . . . , εM ]T であり，j 番目の要素 εj
は j 番目の観測データ yj (t) に重畳するノイズを表す．式 (10) が観測データを表す線形離散モデルである．
データが時刻 t1 , t2 , . . . , tK で計測される時系列データであるとする．表記法を簡単にするため
時刻 tk で計測されたデータベクトルを y k （y(tk ) = y k ）で表し，時刻 tk のソースベクトルを xk（x(tk ) = xk ）と表す．すると，式 (10) は

y k = F xk + ε(11)と表される．時系列データの全体 y 1 , y 2 , . . . , y K を y で表す．
（つまり，y と書いてあれば y 1 , y 2 , . . . , y Kを意味する．
）同様に，x1 , x2 , . . . , xK を x と表す．後の議論のため行列 F の j 番目の列ベクトルを

f j で表す．すなわち，f j を用いるとF = [f 1 , f 2 , . . . , f 3N ]
である．

3

3

スパースベイズ学習

3.1

確率モデル

式 (11) で表されたのモデルのもとで，観測データ y から未知な信号源分布 x を推定する問題を考える．ノイズの確率分布として平均ゼロ，共分散行列 Σ ε の多次元正規分布p(ε) = N (ε|0, Σ ε )(12)を仮定する．式 (11) のデータモデルと式 (12) のノイズの確率分布からデータ y k の条件付き確率分
布 p(y k |xk ) がp(y k |xk ) = N (y k |F xk , Σ ε )

(13)

と求まる．
ベイズ学習ではこれから推定する未知量 x に対する事前分布を導入する．スパースベイズ学習では事前分布として，要素 xj (t) に対して分散が νj で，要素間で独立な正規分布：

p(xj (t)) = N (xj (t)|0, νj )(14)を仮定する．ベクトル確率変数 xk （xk ∈ R3N ×1 ）に対する確率分布はp(xk ) = N (xk |0, Υ )(15)となる．この事前分布の共分散行列 Υ （Υ ∈ R3N ×3N ）はベクトル ν （ν ∈ R3N ×1 ）を ν =[ν1 , ν2 , . . . , ν3N ]T と定義して，
Υ = diag(ν)(16)と表される対角行列となる．つまり Υ はこの ν の要素を対角成分として持つ対角行列である．

3.2

事後確率分布の導出

ボクセルソースベクトル xk のもっとも確からしい推定解は，事後分布 p(xk |y k ) を最大とする xk
として求めることができる（MAP（maximum a posteriori）推定解）．事前分布（式 (15)）と尤度（式 (13)）が正規分布の場合，事後分布 p(xk |y k ) も正規分布で与えられる．したがって，事後分布を平均 x̄k ，精度行列（共分散行列の逆行列）Γ の正規分布

p(xk |y k ) = N (xk |x̄k , Γ −1 )(17)

Γ = Υ −1 + F T Σ −1
ε F

(18)

とすれば，精度行列と平均は

x̄k = Γ

−1

F

4

T

Σ −1
ε yk

(19)

として求まる1 ．したがって，事後分布 p(xk |y k ) を最大とする解，すなわち，未知量 xk の MAP 推
定解は，正規分布では期待値で確率最大となるのであるから，式 (19) で与えられる x̄k である．
式 (19) を等価な別の表現に変換する．逆行列公式（ベイズ信号処理 (A.51)）を用いて，
TT−1x̄k = (Υ −1 + F T Σ ε−1 F )−1 F T Σ −1
y k = Υ F T Σ −1
ε y k = Υ F (F Υ F + Σ ε )
y yk

(20)

が成り立つ．ここで，

Σy = Σε + F Υ F T(21)

であり，この Σ y はモデルデータ共分散行列と呼ぶ．さらに，後の議論のために式 (20) を重み行列を用いて

x̄k = Ωy k

(22)と書く．つまり，事後平均 x̄k は観測データ y k の線形推定解として表され，x̄k を推定する重み Ω はΩ = Υ F T Σ −1y(23)で与えられる．
以上の議論から明らかであるが，未知量 xk のベイズ推定解（MAP 推定解）x̄k を計算するにはハイパーパラメータ Υ （すなわちその対角成分である ν ）とノイズ共分散行列 Σ ε が必要である．

RENSx アルゴリズムではこれらを観測データ y から推定する．以下にまず，ハイパーパラメータ Υの推定ついて説明する．

3.3

ハイパーパラメータ ν の学習

3.3.1

周辺尤度の最大化による ν の学習

A. 周辺尤度の導出
b は時系列データの周辺尤度
ハイパーパラメータ Υ ，すなわちその対角成分 ν の最適推定解 µlog p(y|ν) = log p(y 1 , . . . , y K |ν) =KXlog p(y k |ν)

k=1

を最大とする ν として求めることができる．ここで，log p(y k |ν) は k 番目のタイムサンプルにおける観測データ y k に対する周辺尤度である．この周辺尤度 p(y k |ν) を求めるにはいくつかの方法がある．最も直截的な方法は

p(y k |ν) =ZZp(y k , xk |ν)dxk =p(y k |xk )p(xk |ν)dxk(24)から求めるものである．式 (24) において p(y k |xk ) は式 (13) で与えられ，p(xk |ν) はで与えられるのでこれらを代入し，右辺の積分を計算すれば p(y k |ν) を求めることができる．この正攻法による周辺
1 ベイズ信号処理の 3.4 節を参照のこと．

5

尤度の計算は「ベイズ信号処理」5.3 節に詳しく説明してある．ただし，この積分の実行にはかなりな計算を必要とする．この積分を計算する方法よりも若干計算の楽な別の方法について「信号処理のための線形代数入門」[4] の 8.2.3 項に説明がある．ここでは，周辺尤度 p(y k |ν) のさらに簡単な計算
法を示す．
ここでは，

¤
1
1£
T −1
log p(y k |ν) = − log |Σ y | −(y k − F x̄k )T Σ −1x̄ε (y k − F x̄k ) + x̄ Υ
2
2

(25)

と求まることを証明する．式 (12) に示すように

ε ∼ N (ε|0, Σ ε )を仮定している，また，式 (11) と (15) が成り立つので，

F xk ∼ N (F xk |0, F Υ F T )(26)が成り立つ．確率変数 y k は独立な正規分布に従う 2 つの確率変数 F xk と ε の和で表される．独立な正規分布の和はやはり正規分布であり，

y k = F xk + ε ∼ N (y k |0, F Υ F T + Σ ε ) = N (y k |0, Σ y )である2 ．したがって，

(27)

1
1
log p(y k ) = − log |Σ y | − y Tk Σ −1
y yk
2
2

(28)

−1
T −1
T
x̄k
y Tk Σ −1
y y k = (y k − F x̄k ) Σ ε (y k − F x̄k ) + x̄k Υ(29)を得る．さらに，の関係3 を用いれば結局，時系列データに対する周辺尤度として
K
X

¤
K
1 X£
T −1
log |Σ y | −(y k − F x̄k )T Σ −1x̄k
ε (y k − F x̄k ) + x̄k Υ
2
2
k=1
k=1
(30)
を得る．ここで，Σ y は式 (21) に示すモデルデータ共分散行列である．
log p(y|ν) =K

log p(y k |ν) = −B. ハイパーパラメータ更新式の導出
b は式 (30) においてハイパーパラメータ ν の最適な推定解 ν∂
log p(y|ν) = 0∂ν
とおいて求まる ν である．この微分計算は「ベイズ信号処理」5.4 節で詳しく述べられていおり，事後共分散行列 Σ を Σ = Γ −1 と定義して，Σ の {j, j} 要素を Σj,j と書けば，

X∂log p(y|ν) = −KΣj,j + Kνj −x̄j (tk )
∂νj
K

k=1

2 「統計的信号処理」の 2.2 節参照．
3 証明は「信号処理のための線形代数入門」の問題 8.1 の解答を参照のこと

6

(31)

となる．ここで，x̄j (tk ) はベクトル xk の j 番目の要素である．したがって，上式の右辺をゼロと置くことにより νj の推定解として，

1 X
x̄j (tk )KKνbj = Σj,j +

(32)

k=1

を得る．
（式 (32) は EM アルゴリズムを用いて導出した νj の推定解と同じものとなっている．詳し
くは「ベイズ信号処理」の第 5 章参照．
）式 (32) は推定が劣決定の場合（特にその度合いが激しい場
合，すなわち M << N の場合）に収束が遅いことが知られている．次節では収束の早い ν の学習について述べる．

3.3.2

凸関数の性質を用いた ν の学習

A. 代替コスト関数の導出
式 (32) よりも収束の早い更新式を凸関数の性質を用いて導出する．
（「ベイズ信号処理」5.5 節にここで用いる手法についての解説がある．
）コスト関数を式 (30) の周辺尤度を用いて

¤
1 X£
T −1
x̄k (33)(y k − F x̄k )T Σ −1
ε (y k − F x̄k ) + x̄k ΥKK

F(ν) = −(2/K) log p(y|ν) = log |Σ y | +k=1と定義する．この F(ν) を最小とする ν が ν のベストな推定解である．この最小化に際して，log |Σ y |は変数 ν1 , ν2 , . . . , ν3N に関して上に凸な関数であるので，任意の ν に対して，z0 をある定数としてz T ν − z0 ≥ log |Σ y |(34)

を満たすベクトル z = [z1 , . . . , z3N ]T が必ず存在するという事実を用いる．すなわち，代替コスト関数をe z) = z T ν − z0 + 1
F(ν,
K

K
X
£
¤
T −1
(y k − F x̄k )T Σ −1x̄k
ε (y k − F x̄k ) + x̄k Υ

(35)

k=1

とすれば，

e z) ≥ F(ν)
F(ν,

(36)

e z) は真のコスト関数 F(ν) の上限を形成する．すなわち，代が成り立つので，代替コスト関数 F(ν,ee z) を最小化する ν は真の替コスト関数 F(ν, z) を補助変数 z と共に変数 ν で最小化すれば，F(ν,e z) を最小とする
b は代替コスト関数 F(ν,コスト関数をも最小化する．したがって，ν の最適推定値 νν として求めることができる．
B. ハイパーパラメータ ν の推定式の導出

7

b は，
実際に ν
"
#
K
1 X T −1Teb = argmin F(ν, z) = argmin z ν +νx̄k Υ x̄k
K
ν
ν
k=1

= argminνNX"
zj νj +

1KPK

2
k=1 x̄j (tk )

j=1

νj

#
(37)

として求める．再右辺を νj に関し微分してゼロと置けば，

s
νbj =

1
K

PK

2
k=1 x̄j (tk )

zj

(38)

を得る．この式に示す νbj が νj の推定解である．しかし上式には，補助変数 zj が入っている．補助
変数 zj の推定値は（「ベイズ信号処理」5.5.2 項を参照のこと．
）

zbj =として求まる．ここで，∂log |Σ y |
∂µj

·
¸
∂
−1 ∂
log |Σ y | = tr Σ y
Σy
∂µj
∂µj

(39)

(40)

であり，式 (21) から

∂
∂
Σy =
F Υ F T = f j f Tj∂µj∂µjであるので，hiTT−1
zbj = tr Σ −1
y f j f j = f j Σy f jを得る．式 (38) および (41) が ν の更新式である．式 (38) に (41) を代入してしてv
u 1 PK2u
k=1 x̄j (tk )νbj = t K T −1
f j Σy f j

(41)

(42)

として ν の更新式を表すこともできる．

3.4

ハイパーパラメータ・タイイング

前節で述べたアルゴリズムには問題がある．ボクセルソースベクトル x の j 番目の要素が j = 3×n+1と表せれば，これは n 番目のボクセルの x 成分の値であり，この場合，j = 3×n+2 および j = 3×n+3
ば，n 番目のボクセルの y および z 成分の値である．このようにソースベクトルのインデックス j は異なる 2 つのインデックス (ボクセルのインデックス n と同じボクセルにおけるソースベクトルの 3
つの成分を表すインデックス) の和で表されるが，これまで述べてきたスパースベイズアルゴリズム
はこの違いを考慮せず同等に扱うため，信号源の成分間でスパース制約が働いてしまい，信号源の向8きの推定に（その結果として信号源強度の推定にも）誤差を生じてしまう．したがって，このアルゴリズムによる再構成結果は誤差を含む．この実例は付録 B（Biomag2014 の発表ポスター）を参照されたい．
この問題に対する解決策は同じボクセル由来のハイパーパラメータ ν3n+1 ，ν3n+2 ，ν3n+3 を結び付
けて同じ 1 つの値を持つようにすることである．すなわち，それぞれの更新値 νb3n+1 ，νb3n+2 ，νb3n+3を用いて新しい更新値 µb3n+m （m = 1, 2, 3）を

v
u 3
u1 X
µ
b3n+m = tνb2
3 i=1 3n+i(43)とすることにより，n 番目のボクセルのハイパーパラメータの値を３つの成分に関して同じ値にするのである．これを Hyperparameter Tying4 と呼ぶ．
ハイパーパラメータを結合する，すなわち同じ値にしてしまうことによりパラメータ間のスパー
ス制約をコントロールできる．すなわち，タイイングしたパラメータの間ではスパースな制約は消失するが，タイイングしていないパラメータの間ではスパースな制約が保たれる．ここで議論している式 (43) を用いた場合には，各ボクセルのソースベクトルの成分間ではスパース制約は消失する
が，ボクセル間ではスパース制約はなお有効である5 ．Hyperparameter Tying がスパース制約に与える影響についてのコスト関数を使った議論は付録 B あるいは Electromagnetic Brain Imaging の4.9 節に説明がある．また，補助変数 zbj についても同様に
1X
zb3n+i3 i=1
3new
zb3n+m=(m = 1, 2, 3)(44)とすることによりタイイングを行う．
このハイパーパラメータ・タイイングの有効性に関するコスト関数を用いた議論と実験結果を発表
したポスター（Biomag2014,Halifax, Canada）を付録 B に示す．ポスター６ページ目に，2 次元データの場合におけるハイパーパラメータ・タイイングを行った場合と行わない場合のコスト関数の導出を示す．また，ポスター 7 ページ目にコスト関数のプロットを示す．このプロットから，ハイパーパ
ラメータ・タイイングを行なわない場合にはコスト関数が L0 -あるいは L1 -ノルムに類似のスパースな解を生ずるタイプの形状をしているが，ハイパーパラメータ・タイイングによってこの形状が L2 ノルムのコスト関数に類似したスパースな解を生じないタイプの形状に変化することが示される．

8 ページから 9 ページにかけて，ハイパーパラメータ・タイイングの有効性を示すコンピュータシ
ミュレーションとその結果を示す．9 ページの結果は，(a) 事前分布の共分散行列を対角行列とする
方法，(b) 事前分布の共分散行列を対角行列として，ハイパーパラメータ・タイイングを行う方法，(c) 事前分布の共分散行列を対角行列とする方法（A 節に述べる方法）を比較し，(a) では再構成結
果に大きな誤差を含むが，(b) であれば誤差の無い再構成を得ることができる．また，(c) でも誤差の無い再構成が可能であるが，計算時間が 12 倍（方法 (b) と比べて）かかることが述べられている．
4 ハイパーパラメータの結合とでも訳すべきか?
5 ちなみにベイズミニマムノルム法はこの極端な場合である．この方法では全ボクセルのハイパーパラメータを結合し，1つの値で表していると考えることができる．したがって，ベイズミニマムノルム法ではスパースな解は生じない．なお，ベイズミニマムノルム法については Electromagnetic Brain Imaging の 2.10 節，あるいは技術ノート（参考文献 [5]）を参照の
こと．

9

3.5

ハイパーパラメータのスムージング

SBL アルゴリズムは理論的にはシングルタイムサンプルのデータに対しても有効なアルゴリズムである．しかし，実際にはある程度の数のタイムサンプルを必要とする．この，ある程度の数のタイ
ムサンプルの必要性は式 (38) において明らかである．式 (38) において，ハイパーパラメータ νj の更新式は j 番目のボクセルにおける再構成パワーの測定タイムタイムサンプルでの平均値

1 X 2
x̄j (tk )
K
K

(45)

k=1

を用いて計算される．個々のタイムサンプルにおける値 x̄2j (tk ) は観測データに混入するノイズの影
響を受けるので，ある瞬時値ではなく，式 (45) に示すようにタイムタイムサンプルでの平均値を用
いて νj の更新を行うことにより，観測データに混入するノイズに対しより安定な更新を行うことができる．
脳磁界計測以外の生体磁気では，少ないタイムサンプルでの信号源再構成を行わなければならな
い場合も多い．例えば，脊髄や神経磁場の計測では信号源が高速で移動するため，通常はシングルタ
イムサンプルで信号源再構成を行う．ハイパーパラメータのスムージングは，このような少数タイム
サンプルを用いた再構成において少しでも安定した再構成を行うことを目的とする．

b
ハイパーパラメータスムージングは，式 (43) を用いて計算したハイパーパラメータの更新値 µ（= [bµ1 , . . . , µb3N ）を，あるスムージング重み関数 Θ を用いて，

b S = Θµ
b
µ

(46)

S

b を最終的な更新値として用いる．ハイパーパラメー
としてスムージングを行うことにより求めた µ
タスムージングを行うことにより，少数タイムサンプルの場合においてもノイズの影響を含む更新値
を周囲のボクセルにおける重み付き平均を計算することにより，更新値に含まれるノイズの影響を低減し，アルゴリズムの安定性の向上が期待できる．

4

ノイズの学習
ノイズ共分散行列をハイパーパラメータとみなせば，周辺尤度関数の最大化から，ノイズ共分散行列も推定することができる．ノイズ共分散行列の推定をノイズ学習（noise learning）と呼ぶ．ノイズ共分散行列の最適解は，周辺尤度を最大とするノイズ共分散行列として，すなわち∂∂ X
log p(y|ν, Σ ε ) =log p(y k |ν, Σ ε ) = 0
∂Σ ε
∂Σ ε
K

(47)

k=1

を満たす Σ ε として求めることができる．ただし，この計算のためにはコスト関数を直接最大化するよりも EM アルゴリズムを用いて平均データ尤度を最大化する方が計算が楽である．

10

ノイズ共分散行列 Σ ε ではなく精度行列 Λ（= Σ −1ε ）を用いて平均データ尤度を表すと，

Θ(Λ) =KX

E [log p(y k |xk , Λ)] +

k=1

K
X

E [log p(xk |Υ )]

(48)

k=1

ここで，E [·] は事後分布の期待値を表す．Λ の更新値は∂
Θ(Λ) = 0∂Λを満たす Λ である．式 (48) において右辺第 2 項は Λ を含まずこの計算には関係しない．したがって，

·
¸
K
X
∂
∂
Θ(Λ) =
E
log p(y k |xk , Λ)
∂Λ
∂Λ

(49)

k=1

であり，

1
1
log |Λ| − (y k − F xk )T Λ(y k − F xk )22である．式 (50) では計算に関係のない定数は書くのを省略した．ここで，であるので，

log p(y k |xk , Λ) =

(50)

∂
1
1
log p(y k |xk , Λ) = Λ−1 − (y k − F xk )(y k − F xk )T
∂Λ
2
2

(51)

K
¤
∂
K
1X £
E (y k − F xk )(y k − F xk )TΘ(Λ) = Λ−1 −
∂Λ
2
2

(52)

k=1

を得る．さらに，£¤E (y k − F xk )(y k − F xk )T = (y k − F x̄k )(y k − F x̄k )T + F Γ −1 F T(53)であるので，結局，Λ すなわち Σ ε の更新式として

bε = Λ
b −1 = 1
ΣKK
X(y k − F x̄k )(y k − F x̄k )T + F Γ −1 F T

(54)

k=1

を得る．
式 (54) をそのまま使うのではなく，計算しやすいように変形する．式 (54) の右辺第 1 項は式 (22)
を用いれば

1 X
1 X
(y k − F x̄k )(y k − F x̄k )T =(y k − F Ωy k )(y k − F Ωy k )T = ΠRyy Π T
K
K
K

K

k=1

k=1

(55)

である．ここでΠ = I − F Ω = I − F Υ F T Σ −1
y
Ryy =

1
K

K
X

y k y Tk

k=1

11

(56)
(57)

である．
式 (54) の右辺第 2 項については，まず，Γ −1（式 (18)）を以下のように変形する．Γ −1 は，逆行
列公式より，³´−1
Γ −1 = Υ −1 + F T Σ −1
= Υ − Υ F T (Λ + F Υ F T )−1 F Υ = Υ − Υ F T Σ −1
ε F
y FΥ

(58)

を得る．したがって，式 (54) の右辺第 2 項は³´
TTF Γ −1 F T = F Υ − Υ F T Σ −1
= F Υ F T − F Υ F T Σ −1
y FΥ Fy FΥF³´= I − F Υ F T Σ −1
F Υ F T = (I − F Ω) F Υ F T = ΠF Υ F T
y(59)となるので，結局，Σ ε の更新式 (54) はb ε = ΠRyy Π T + ΠF Υ F T
Σ(60)と表される．通常は ε を対角ノイズ（diagonal noise）と仮定しているので，Σ ε を対角行列として，
更新式は

h
i
b ε = diag ΠRyy Π T + ΠF Υ F T
Σ

(61)となる6 ．

5

SBL ビームフォーミング
SBL ビームフォーミングは SBL アルゴリズムから出力されたモデルデータ共分散行列 Σ y を用いてビームフォーマーを構成する信号源再構成手法である．従来のビームフォーマーは，サンプルデー
タ共分散行列 R を用いて構成した重み行列：hi−1
W (r) = R−1 L(r) LT (r)R−1 L(r )を用いて，位置 r におけるソースベクトルを
hi−1bs(r, t) = W T (r)y(t) = LT (r)R−1 L(r )LT (r)R−1 y(t)

(62)

(63)

として推定する．SBL ビームフォーミングは，サンプル共分散行列 R の代わりに，SBL アルゴリズ
ムの出力として求まるモデルデータ共分散行列 Σ y を用いて重み行列

i−1
h
T
−1
W (r) = Σ −1y L(r) L (r)Σ y L(r )(64)

を求め，位置 r におけるソースベクトルを

i−1hb
LT (r)Σ −1s(r, t) = W T (r)y(t) = LT (r)Σ −1
L(r)
y y(t)y(65)

6 ここで，行列を A として B = diag(A) とは，「A の対角成分を対角成分に持つ対角行列を B とする」の意味である．

12

ν

σ

xk

νj
νj

zj

σ

Σy

xk

xk

図 1: RENSx アルゴリズムの流れ．ベクトル σ はノイズ共分散行列 Σ ε の対角成分を要素に持つ列ベクトルであり，Σ ε = diag(σ) の関係がある．

と推定する．

SBL ビームフォーマーは次の 2 つの点で，従来のビームフォーマーよりも優れている．サンプル
共分散行列 R を精度良く求めるためには多数のタイムサンプル（すなわち大きな K の値）を必要とする．そのため，従来のビームフォーマーは多数のタイムサンプルを含むデータ時間窓を必要とす
る．これがためにビームフォーマーは多数のタイムサンプルを含むデータ時間窓を確保できない計測
データ，例えば誘発脳磁界データなどにはあまり用いられることはなかった．一方で，SBL アルゴリ
ズムによるモデルデータ共分散行列の推定ではこのような制約は存在せず，少数のタイムサンプルからもモデルデータ共分散行列の推定が可能である．
さらに，ビームフォーマーには複数の信号源活動に時間的な相関があると信号源推定結果に大きな誤差を生じることが知られている．これは，データサンプル共分散行列が信号源相関の影響を大きく受け，その信号およびノイズ部分空間が縮退してしまうためである．一方，SBL アルゴリズムにより求まるモデルデータ共分散行列は信号源活動の相関には影響を受けなため，信号源活動に大きな相関があっても信号源活動の情報を保ったモデルデータ共分散行列の推定が可能である．

13

付録 C に SBL ビームフォーマーに関する実験結果を発表したポスター（Biomag2016,Soeul,Korea）を示す．このポスターは特に SBL ビームフォーマーの multi-resolution capability に関するコンピュー
タシミュレーションの報告であり，low-resolution ボクセル（ボクセル間隔：15mm）でモデルデー
タ共分散行列を求め，high-resolution ボクセル（ボクセル間隔：2mm）でビームフォーマ再構成を行うことで，low-resolution ボクセルでの再構成よりはるかに高分解能の再構成結果を得ることができることを報告している．

6

RENSx アルゴリズム−まとめ
RENSx アルゴリズムを図 1 にまとめた．RENSx アルゴリズムはまず SBL アルゴリズムを実行する．この部分では，まずハイパーパラメータ ν と σ に初期値を与える．ここで，σ はノイズ共分散
行列 Σ ε の対角成分を要素に持つ列ベクトルであり，ノイズ共分散行列は Σ ε = diag(σ) として求める．次に，ν と Σ ε を用いて事後平均 x̄k （k = 1, . . . , K ）を求める．さらに，この x̄k を用いて，ν
と Σ ε を更新する．以上の手順をある条件が満たされるまで行う．SBL アルゴリズムの出力は事後
平均 x̄k （k = 1, . . . , K ）とモデルデータ共分散行列 Σ y （式 (21)）である．RENSx アルゴリズムにおいでは，SBL アルゴリズムの出力であるモデルデータ共分散行列 Σ y を用いて SBL ビームフォー
ミングを行う．アルゴリズムの最終的な出力は SBL から出力される SBL ボクセルタイムコース x̄k（k = 1, . . . , K ）と，ビームフォーマーの出力 x̄k （k = 1, . . . , K ）である．

7

参考文献
以下に SBL アルゴリズムの説明がある．

1. Electromagnetic Brain Imaging, Springer 2015.
2. ベイズ信号処理，共立出版 2015.
3. 信号処理のための線形代数入門，2019.
4. シグナルアナリシス社技術資料「スパースベイズ推定」
http://www.signalanalysis.jp/pdfs/SBL_note.pdf
5. シグナルアナリシス社技術資料「ベイズミニマムノルム推定」
http://www.signalanalysis.jp/pdfs/BMN_note.pdf

14

A

Appendix: マトリックス SBL アルゴリズム

A.1

確率モデルと学習測本ノートの 3 節で述べた SBL アルゴリズムは各ボクセルの共分散行列をスカラー倍の単位行列に仮定したおり，スカラータイプの SBL アルゴリズムと呼ばれる．これに対し，ボクセル共分散行列に特別な仮定を置かない SBL アルゴリズムを構築することも可能である．このアルゴリズムは（いわば）マトリックスタイプの SBL アルゴリズムであり，matrix champagne と呼ばれる．
このアルゴリズムでは各ボクセルにおける j 番目のソース s(r j , t) に対して事前分布を，p(sj (t)) = N (sj (t)|0, Υ j )(66)とする．簡単のため s(r j , t) を sj (t) と書いた．ここで，共分散行列 Υ j は 3 × 3 の行列である．各
ボクセルにおけるソースベクトル sj (t) は（固定の）任意の向きを取りえるので対角行列と仮定はできない．したがって，ボクセルベクトル xk に対する事前分布はp(xk |Υ ) =NY

N (sj (tk )|0, Υ j ) = N (xk |0, Υ )

(67)

j=1

で表される．ボクセルベクトルの事前共分散行列 Υ は 3N × 3N のブロック対角行列で



Υ1


 0
Υ =
 ..
 .
0

0

···

Υ2

···
..
.

0
0

···

0




0 
.. 
.
. 
ΥN

(68)

と表される．
この事前分布と式 (12) のノイズの仮定のもとで，事後分布の平均はやはりx̄k = Υ F T Σ −1
y yk

(69)

と求まる．個々のボクセルにおけるソースベクトルの事後平均はs̄j (tk ) = Υ j LTj Σ −1
y yk

(70)

である．ここで，モデルデータ共分散行列 Σ y は

Σy = Σε +NX
j=1である．ただし，Lj = L(r j ) と書いた．

15

Lj Υ j LTj(71)

式 (33) に対応したコスト関数はF(Υ ) = −(2/K) log p(y|Υ )



K
N
X
1 X

s̄Tj (tk )Υ −1= log |Σ y | +(y k − F x̄k )T Σ −1
ε (y k − F x̄k ) +j s̄j (tk )
K
j=1

(72)

k=1

である．ここで，式 (34) に対応して，
N
X³´
tr Z Tj Υ j − Z0 ≥ log |Σ y |

(73)

j=1

を満たす 3 × 3 の行列 Z j （j = 1, . . . , N ）が必ず存在する．したがって，代替コスト関数を以下のように定義すれば7



K
N
X
X
1
e , Z) =(y k − F x̄k )T Σ −1F(Υs̄Tj (tk )Υ −1
ε (y k − F x̄k ) +j s̄j (tk )
K
j=1
k=1

+

N
X

´
³
tr Z Tj Υ j − Z0

(74)

j=1

このコスト関数は必ず

e , Z) ≥ F(Υ )F(Υe , Z) を Υ と Z で最小化すれば，その Υ は必ず F(Υ ) も最小化する．
であり，F(Υハイパーパラメータ Z j と vU pj の更新式は以下のように導出される．
（これらの導出は Electro-

magnetic Brain Imaging の 4.8 節にある．
）ハイパーパラメータ Z j の更新式は∂log |Σ y |
∂Υ j

(75)

b j = LT Σ −1 Lj
Z
j
y

(76)

bj =
Z
から求まり，右辺の計算を行うとである．また，Υ j の更新式は

"

³
´
X
T
b j = argmin 1
s̄
(t
)
+
tr
Z
Υ
s̄Tj (tk )Υ −1
Υ
j k
j
j
j
K
Υj
K

#
(77)

k=1

から求まる．右辺の最小値を与える Z j は

"
"
#
#1/2
K
1 X
−1/2
1/2
1/2
−1/2
T
b
Υ j = ZjZjs̄j (tk )s̄j (tk ) Z j
Zj
K
k=1

7Z

1 , . . . , Z N をまとめて Z す．

16

(78)

と求まる．
ここで述べた matrix champagne はボクセル共分散構造に特別な仮定を置かない，言ってみれば理論的には正統な SBL アルゴリズムであるが，第 3.3.2 節で述べたスカラータイプの SBL アルゴリ
ズムに対して，明らかにアルゴリズムが複雑（スカラー計算が行列計算になる等）になり，その分計算時間を必要とする8 ．その一方で，（理論的な正当性はともかく）実用的にはスカラータイプの SBL
アルゴリズムに対して特段の利点は見当たらない．

A.2

Hyperparameter Tying のマトリックス SBL からの導出

マトリックスタイプの SBL はボクセルの共分散構造に何の仮定も置かないで導かれるアルゴリズムである．それでは，ボクセルの共分散行列をスカラーと仮定し，第 3.4 節で導入した Hyperparameter

Tying と等価な式が更新式として導かれることを示す．
まず，ボクセルベクトルの事前共分散行列 Υ は 3N × 3N のブロック対角行列で




Υ =



µ1 I 3

0

···

0
..
.
0

µ2 I 3

···
..
.

0
0

···

0




0 
.
..


.
µN I 3(79)と仮定する．ここで，I 3 は 3 × 3 の単位行列を表す．各ボクセルのソースベクトル sj (t) の事後平均
は，したがって，s̄j (tk ) = µj LTj Σ −1
y yk

(80)

となる．モデルデータ共分散行列は，この場合，

Σy = Σε +NXµj Lj LTj

(81)

j=1

である．
したがって，コスト関数はF(µ1 , . . . , µN ) == log |Σ y | +

K
1 X

K


(y k − F x̄k )T Σ −1
ε (y k − F x̄k ) +

N
X


T

µ−1
j s̄j (tk )s̄j (tk )

(82)

j=1

k=1

となる．式 (73) に Υ j = µj I 3 と Z j = zj I 3 を代入すれば，NX

zj µj − Z0 = z T µ − Z0 ≥ log |Σ y |j=1
8 以前の検討では計算時間はスカラーアルゴリズムに比べ 10 倍程度必要であった．

17

(83)

の関係を満たす．ここで，z = [z1 , . . . , zN ] および µ = [µ1 , . . . , µN ] とした．したがって，代替コス
ト関数を



K
N
X
X
1
−1
T
e
(y k − F x̄k )T Σ ε (y k − F x̄k ) +F(µ,z) = z T µ − Z0 +µ−1j s̄j (tk )s̄j (tk )
K
j=1

(84)

k=1

とすれば，やはり，eF(µ,z) ≥ F(ν)(85)

eが成り立つので，F(µ,z) を最小とする µ = [µ1 , . . . , µN ] は真のコスト関数 F(µ1 , . . . , µN ) を最小とする．

e
実際，F(µ,z) を最小とする µ は



K X
N
X
1
T
e

b = argmin F(µ,µz) = argmin z T µ +µ−1
j s̄j (tk )s̄j (tk )
K
µ
µ
j=1
k=1

として求める．

∂ e
1 X Ts̄j (tk )s̄j (tk ) = 0
F(µ, z) = zj − µ−2
j
∂µj
K
K

k=1

より，s̄j (tk ) = [s̄xj (tk ), s̄yj (tk ), s̄zj (tk )] とすれば，結局，

s
µ
bj =

1
K

s

PK

T
k=1 s̄j (tk )s̄j (tk )

zj

=

1
K

¤
PK £ xyz
222
k=1 s̄j (tk ) + s̄j (tk ) + s̄j (tk )
zj

(86)

e
を得る．また，F(µ,z) を最小とする z は式 (39) から求まる．式 (81) から∂
Σ y = Lj LTj∂µjであるので，式 (40) を用いて，
i
i
h
h
y T
x T
z T
T
T
−1
−1 x
−1 y
−1 z
L
Σ
zbj = tr Σ −1
L
L
=
tr
L
j = (lj ) Σ y lj + (lj ) Σ y lj + (lj ) Σ y lj
j j
y
j
y

(87)

£
¤
を得る．ここで，Lj = lxj , lyj , lzj とした．式 (86) に (87) を代入して，結局vuuνb = tj¤
PK £ xy1z222
k=1 s̄j (tk ) + s̄j (tk ) + s̄j (tk )
K
y T
x
z T
−1 y
−1 z
(lxj )T Σ −1
y lj + (lj ) Σ y lj + (lj ) Σ y lj

(88)と求まる．式 (88) がハイパーパラメータ νj の更新式である．
本節で述べた式 (86) は式 (43) と，式 (88) は式 (44) と同じものであり，結局，3.4 節で導入した

Hyperparameter tying を，Υ を式 (79) と仮定することにより，matrix タイプの SBL から導いたこ
とになる．

18

付録B

Hyperparameter tyingに関する発表
Biomag2014におけるポスター

Hyper‐parameter tying: a novel method ofcontrolling the sparsity in sparse Bayes source
imaging (Champagne) algorithm
Sekihara K. 1

Nagarajan S. S.2

1 Department of Systems Design and Engineering, Tokyo

Metropolitan University, Hino, Tokyo 191‐0065, Japan.
2Biomagnetic Imaging Laboratory, Department of Radiology,

University of California, San Francisco, 513 Parnassus Avenue,
S362, San Francisco, CA 94143, USASummary: Part I
The sparse Bayes source imaging (Champagne) algorithm has
proven to be effective for neuromagnetic imaging [1]. The method
produces sparse solutions because a sparsity constraint is
integrated in its cost function. However, there are some occasions
where the sparsity constraint rather leads to wrong solutions. One
such occasion is the estimation of voxel source orientations.
In Champagne algorithm, when the voxel source orientation is
known, a diagonal covariance matrix can be used for the Gaussianprior distribution, and this leads to a simple and efficient algorithm.
When the voxel source orientation is unknown, a block‐diagonalcovariance matrix must be used, which prevents from the sparsity
constraint being imposed on each component of a source
vector. However, the use of such a non‐diagonal covariance matrix
leads to an algorithm that is significantly complex andcomputationally expensive.

Summary: Part II
In this paper, we propose a simple method to overcome thisdifficulty, called the hyper‐parameter tying. The proposed
method can use a diagonal covariance matrix even when thesource orientation is unknown. However, in each step of updating
iteration, the source variance parameters of vector components
are tied together to correctly estimate the source orientation.
Since a diagonal covariance matrix is used, the algorithm keeps
the simplicity and computational efficiency. We discuss why thehyper‐parameter tying can control the sparsity by presenting thecost function analysis of the Champagne algorithm[2].
References
Reference[1] Wipf DP, Owen JP, Attias HT, Sekihara K, Nagarajan SS.,
Neuroimage. 2010; 49(1) :641‐55.
[2] Wipf DP, Nagarajan SS, Advances in Neural Information
Processing Systems, 2008

Sparse Bayes source imaging (Champagne) algorithm
Data model:y  Fx  εProbability model: p( y | x)  N ( y | Fx,  2 I)p( x)  N ( x | 0, ¡ )Covariance matrix of the prior distribution¡  diag[ , , , ]xˆ  arg min F ( x) : F ( x) 

1



2

1

2

y - Hx L2 regularized minimum‐norm method¡  diag[ 1, 2 , , N ] diag(n )n  [ 1, 2 , , N ]Tnˆ  arg min F (n ) : F (n )  log | Σ |  yT ΣyΣ   2 I  F¡ FT
xˆ  ¡ F T Σˆ 1 y
Champagne algorithmx2

When voxel orientation is unknown:

F  [ , lx (rj ), ly (rj ), lz (rj ) , ]
lead field for the j th voxel¡  diag[ 1, 2 , , N ]Leading to the shrinkage over the source
vector components.

¡  diag[¡1, ¡ 2 , , ¡ N ]¡ j : 3  3 submatrix

Leading to a computer‐intensive, slow
algorithm, called the matrix champagne.

Hyperparameter tying:  diag[ ,  3 j 1 , 3 j 2 , 3 j  3 , ]

hyperparameters for the j th voxel

Set the same update value such as ˆj  (1 / 3)m3 1 3 j mto  3 j 1, 3 j 2 , 3 j  3

Cost function for sparse Bayes source imaging (Champagne algorithm)xˆ  arg min F ( x) : F ( x) 

1



2

2

y - Hx   ( x)

( x)  min( xT ¡ 1 x  log | S |)ν

If we assume that the columns of H are orthogonal,Nlog | S |  log( 2   j )j 1

2‐dimensional case2

x 2j

(x 1, x 2 )  min  (  log( 2   j ))
1, 2 j 1 
j




x j2 | x j |2 | xj |222
  log   
x j  4  j 1 22
x j 2  4 2  | x j |



2

(x 1, x 2 )  min(x 12  x 22

No hyperparametertying

 2 log( 2   ))Hyperparametera  a 2  8a 2  
4a222  : (a  x 1  x 2 ) tying 2 log   
2
2


4



 a  8a  a

Plots of cost functions

x2
x1

x1

x1

x1

x2

(a) Champagne cost function withtying hyperparameters.
(b) Champagne cost function withno hyperparameter tying.
(c) Cost function for L2‐norm
regularization.
(d) Cost function for L1‐norm
regularization.

(x 1, x 2 )  x 12  x 22 : L2-norm regularization
(x 1, x 2 ) | x 1 |  | x 2 |: L1-norm regularization

Computer Simulation
Source time courses

Signal magnetic field• 275 CTF sensor array• signal to sensor noise ratio:16

Results of Computer Simulation
Block‐diagonal covariance used

Diagonal covariance used

Diagonal covariance with hyperparameter tying

Champagne with a block‐diagonal covarianceneeds twelve times longer computation time,compared to the champagne algorithm with adiagonal covariance matrix.

付録C

SBL beamformerに関する発表
Biomag2016におけるポスター

Multi‐resolution champagne beamformer
Sekihara K. 1,2

Nagarajan S. S.3

1 Tokyo Medical and Dental University, Tokyo, Japan
2 Signal Analysis Inc., Tokyo, Japan
3 University of California, San Francisco, USASummary
The sparse Bayesian source imaging (Champagne) algorithm has
proven to be effective for neuromagnetic imaging [1,2]. It is robust to
various types of noise and interference, and it can reconstructmultiple sources, even when they are highly correlated [2]. However,
one serious problem of the algorithm is that it is computationallyexpensive because of its iterative nature.
In this paper, we propose a novel algorithm, called multi‐resolution
champagne beamformer. In the proposed algorithm, the first step
estimates the model data covariance by using the sparse Bayesian
scheme, and the second step implements the adaptive beamformer
source imaging using the data covariance estimated in the first step[3].
The key point of the algorithm is to exploit our empirical findings thatthe quality of the estimated data covariance does not very muchdepend on the voxel size in the first step, while the voxel size primarily
determines the quality of the source reconstruction in the second
step.

Therefore, we can significantly shorten the computational time byusing a low‐resolution voxel grid in the first step but still obtain thehigh‐quality source images by using a high‐resolution grid in thesecond step.
Computer simulation verifies the performance of the algorithm. Unlike
the conventional adaptive beamformer, which uses a sample
covariance matrix, the proposed champagne beamformer can beapplied to a case where correlated sources exist and a case where only
a short time window with a small number of time sample is available.
The algorithm may be useful in non‐brain electrophysiological imaging
such as cardiac imaging in which problems of correlated sources
and/or a short time window are common.
References[1] Wipf DP et al. , Neuroimage. 2010; 49(1): 641‐55.
[2] Sekihara S and Nagarajan SS, Electromagnetic brain imaging: aBayesian perspective, Springer, 2015.
[3] Wipf DP and Nagarajan SS, Proceedings of the 24th international
conference on Machine learning (pp. 1023‐1030). ACM.

Sparse Bayes source imaging (Champagne) algorithmy  Fx  εData model:

Probability model: p( y | x)  N ( y | Fx,  2 I)p( x)  N ( x | 0, ¡ )The hyper parameter ¡ is obtained through maximizing
the marginal likelihood

log p ( y | ¡ )  

1
T
1
|
S
|
S
y
y  ,yy
2where the model data covariance is defined as

Sy   2 I  F ¡ F T
This maximization can be implemented in several ways, and the one
we have developed is described in [1,2].

Champagne beamformer
The champagne beamformer uses the model data covariance derived fromthe champagne algorithm[3]:Sy1l(r )w(r )  Tl (r )Sy1l(r )The conventional adaptive beamformer uses the sample data covariance R
derived from a data time window containing many time points[4]:R 1l(r )w(r )  Tl (r )R 1l(r )[4] Sekihara K and Nagarajan SS, Adaptive spatial filters for electromagnetic brain imaging.
Springer, 2008.

Data generation for computer simulation
Computer simulation geometry

Time courses assigned to thethree sources

Simulated sensor recordings

275 (CTF) whole head sensor array
was assumed for data generation.

The model data covariance was estimated with voxel grid of 1.5cm.
The beamformer scan was performed with voxel grid of 0.2cm.

Results of Computer Simulation with no array mismatch
Champagne beamformer images with 2mm voxels

Champagne images with 1.5cm voxels

SNR:8

SNR:4

coordinate of three sources: [0 ‐1 7], [0. 1 9], [0. 1 5] (cm)Results of Computer Simulation with array mismatch
Champagne beamformer images with 2mm voxels

Champagne images with 1.5cm voxels

SNR:8

SNR:4

coordinate of three sources: [0.5 -1.75 9.75], [0.5 2

9.75], [0.5 1.25 6] (cm)Reconstruction of perfectly correlated sources
Champagne beamformer images

Conventional beamformer images

no source
correlation

nearly perfectsource correlation

The first and the second (the upper two) sources are nearly perfectly correlated
in the bottom case.

Comparison of computational timechampagne

champagne beamforming

Condition # of iteration: 400# of time points: 1200
# of sensor channels: 275
reconstruction volume: 9cmX9cmX9 cm

Intel® Core i7 3.4GHz
32GB memory
64 bit windows 7
matlab 2014b