1-10-12
時系列事前分布モデルとスペクトル基底の同時適応を用いた
バイノーラル音源分離の実験的評価 ∗◎室田勇騎 (奈良先端大)，北村 大地 (総研大)，小山 翔一，猿渡 洋 (東京大学)，
中村 哲 (奈良先端大)1 はじめに
近年，信号処理の分野において音源分離技術が盛んに研究されている．この技術は，3D オーディオシステ
ムの実現や自動採譜への応用など，様々な可能性を秘めている．本稿では特にバイノーラル形式の音源分離について述べる．バイノーラル信号における音源分離
では，音の臨場感を保つために，分離音の定位や残響などは保持されていることが望ましい．このような分離を行うためには，頭部伝達関数 (head related transferfunction: HRTF) などの両耳情報を利用することが考えられる．しかし各ユーザの HRTF 情報は基本的に未知であるため，本 HRTF の利用は実用上において現実的ではない．そこで著者らは，事前分布パラメータ推定を用いたゲイン共通型一般化平均二乗誤差最小化短時間振幅スペクトル推定法 (minimum mean-square error
shot-time amplitude estimator: MMSE-STSA) によるバイノーラル音楽信号分離手法を提案している [1]. この手法は，教師あり非負値行列因子分解 (supervisednonnegative matrix factorization: SNMF)[2] による動的な妨害音推定，高次統計量を用いた事前分布パラメータ推定手法，そして両耳で用いるゲイン関数を共通化するゲイン共通型一般化 MMSE-STSA 推定器の 3 処理から構成されている．この手法では，両耳観測データから分離対象音の時系列情報である事前分布モデルを左右個別に推定し，分離に利用する．これにより，HRTF 等の両耳情報を用いることなく，精度のよい分離を行うことができる．
一方で前述の手法では，時系列方向の適応に比べ，周波数方向の適応は十分に行われていなかった．本
手法では，まず分離対象音のスペクトル基底情報をSNMF によって事前に学習し，それを用いて観測音から妨害音成分推定を行う．しかし，事前学習によって得られた教師成分と実際に観測される目的音信号との差異により，推定精度が劣化する可能性がある．
この問題を解決するために，観測音信号に合わせて教師基底を変形し分離を行うことが考えられる．このような手法は基底変形型 SNMF[3] として提案されているが，教師基底の変形と信号の分離を同時に行うため最適化が難しく分離精度が良くないという問題点があった．そこで本研究では，この教師基底を中間的な音源分離音 (MMSE-STSA 推定器出力) に適応させることにより，時系列事前分布モデルとスペクトル基底を反復的に同時適応する手法を提案する．その後，評価実験を行うことで提案手法の有用性を確かめる．

2 関連技術
2.1 事前分布パラメータ推定を用いたゲイン共通型
一般化 MMSE-STSA 推定器によるバイノーラ
ル音楽信号分離
本章では著者らが提案した，事前分布パラメータ推定を用いたゲイン共通型一般化 MMSE-STSA 推定法によるバイノーラル音楽音源分離手法について述べる．Figure 1 に従来手法におけるバイノーラル信号モデルを示す．この手法では，左右の耳に到来する目的音信号振幅スペクトルの確率密度関数 (probability∗Fig. 1 Conventional and proposed binaural models.
density function: p.d.f.）の差を両耳情報として利用し分離を行うことにより，品質の良い分離音を得ることができる．具体的には，伝達関数が畳み込まれた目的音信号をカイ分布を用いた統計的な信号モデルで表すことで，従来の決定論的な HRTF の推定問題を統計パラメータの推定問題に変換する．そして観測デー
タから左右チャネル毎にこの統計パラメータを推定し分離に利用する．この手法では，SNMF による動的な妨害音推定，高次統計量を用いた閉形式の目的音事前分布パラメータの推定手法，及び音像定位の改善を目的とした両耳共通スペクトルゲインを用いている．次節では一般化 MMSE-STSA 推定器，SNMFによる妨害音推定について説明する （両耳共通スペクトルゲインの導出，目的音事前分布パラメータ推定法については文献 [1] を参照されたい）．
2.2 信号の混合モデル
本研究では，左耳と右耳に一個ずつの 2 マイクロホンで観測されたバイノーラル信号を考える．このとき時間周波数領域における観測音信号 x( f, τ) =[xL ( f, τ), xR ( f, τ)]T は，目的音信号 s( f, τ) と伝達関数 h( f ) = [hL ( f ), hR ( f )]T 及び妨害音信号 n( f, τ) =[nL ( f, τ), nR ( f, τ)]T を用いて以下で表されるとする．
x( f, τ) = h( f )s( f, τ) + n( f, τ)(1)
ここで， f は周波数ビン，τ は時間フレームを表し，下付き添え字 ∗ (∗ = {L, R}) はそれぞれ左耳と右耳での信号を表すとする．

2.3 一般化 MMSE-STSA 推定器 [4]一般化 MMSE-STSA 推定器は，目的音を表す事前分布の基で，真の目的音振幅スペクトルとその推定値との平均二乗誤差を最小化するスペクトルゲインを求める．このスペクトルゲインを観測信号スペクトルに乗算することで，推定目的音が得られ，音源分離が達成される．一般化 MMSE-STSA 推定器における推定目的音を s̃∗ ( f, τ) とすると，これは以下の式で与えられる．
s̃∗ ( f, τ) = G∗ ( f, τ)x∗ ( f, τ)(2)G∗ ( f, τ)√()1/βν∗ ( f, τ) Γ(ρ∗ +0.5) Φ(0.5−ρ∗ , 1, −ν∗ ( f, τ))··=γ∗ ( f, τ)Γ(ρ∗ )
Φ(1−ρ∗ , 1, −ν∗ ( f, τ))(3)“Experimental evaluation of simultaneous adaptation method of time-series a priori statistical model and spectral basis
for binaural source separation,” by Yuki Murota (NAIST), Daichi Kitamura (SOKENDAI), Shoichi Koyama, Hiroshi
Saruwatari (The University of Tokyo), and Satoshi Nakamura (NAIST).

日本音響学会講演論文集- 549 -

2015年3月

ここで Γ(·) はガンマ関数， Φ(a, b; k) = F 1 (a, b; k) は合流型超幾何関数，β は振幅圧縮パラメータを表す．また，ν∗ ( f, τ) は以下の式で表される()−1ν∗ ( f, τ) = γ̃∗ ( f, τ)ξ̃∗ ( f, τ) 1 + ξ̃∗ ( f, τ)(4)ξ̃∗ ( f, τ)，γ̃∗ ( f, τ) はそれぞれ事前 SNR，事後 SNR であり，以下の式で表される．
ξ̃∗ ( f, τ) = αγ̃∗ ( f, τ − 1)G2 ( f, τ) + (1 − α)max[γ∗ ( f, τ) − 1, 0](5)γ̃∗ ( f, τ) =|x∗ ( f, τ)|2 /Pñ∗ ( f )(6)

ここで，Pñ∗ ( f ) は推定妨害音のパワースペクトル，αは忘却係数を表す．
一般化 MMSE-STSA 推定器における目的音振幅ス
ペクトルの事前分布は，次のカイ分布で表される．

p(x) = 2ϕρ Γ(ρ)−1 x2ρ−1 exp(−ϕx2 )ρ
ϕ=E{|x|2 }

(7)
(8)

ここで p(x) は信号の振幅スペクトルの p.d.f. であり，
ρ は形状母数である．カイ分布において，ρ = 1 であるとき，これはレイリー分布と一致し，時間領域では信号が複素ガウス分布に従うと仮定している．また 1より小さい ρ では，優ガウス性の分布に従うと仮定している．
一般化 MMSE-STSA 推定器では， γ˜∗ ( f, τ) の計算に妨害音のパワースペクトルが必要となるが，妨害音が非定常である場合，動的な推定が必要となる．また，目的音の事前分布を最もよく表すための最適な形状母数についても推定する必要がある．

2.4 SNMF による妨害音推定
SNMF は，事前学習によって得られた教師基底F( f, k) を用いて観測信号のスペクトログラムを近似分解する手法である．SNMF による観測スペクトログラムの分解は以下の式で表される.
∑A∗ ( f, τ) = |x∗ ( f, τ)| ≈ k F∗ ( f, k)V∗ (k, τ)∑+ n H∗ ( f, n)U∗ (n, τ) (9)
ここで，F∗ ( f, k) は事前学習された目的音信号のスペクトルパターンを含む基底行列の要素値である．V∗ ( f, k)
は F∗ ( f, k) に対応するアクティベーション行列の要素値であり，各スペクトルパターンの時間強度変化を示す．H∗ ( f, n) 及び U∗ (n, t) はそれぞれ，目的音信号以外の成分を表すための基底行列及びアクティベーション
行列の要素値であり，F∗ ( f, k)，V∗ (k, t)，H∗ ( f, n)，及び U∗ (n, t) はいずれも非負の実数値となる．また, k はF∗ ( f, k) の基底のインデックスを表し，n は H∗ ( f, n) の基底のインデックスを表す．ここで，事前学習で得られる教師基底 F∗ ( f, k) は，目的音信号のみが含まれる教師信号を用いて構成される．この分解では F∗ ( f, k)を固定した状態で V∗ (k, τ)∑，H∗ ( f, n)，及び U∗ (n, τ) を求めるため，理想的には k F∗ ( f, k)V∑∗ (k, τ) は観測音信号に含まれる目的音信号の成分を， n H∗ ( f, n)U∗ (n, τ)は目的音信号以外の成分をそれぞれ表している．
∑こ れ よ∑り 得 ら れ たn H∗ ( f, n)U ∗ (n, τ) (又 は ，A∗ ( f, τ)− k F∗ ( f, k)V∗ (k, τ)) を推定妨害音の振幅スペクトルとして一般化 MMSE-STSA 推定器に用いる．
この手法は教師基底を用いるので，音楽信号のように非定常で動的な妨害音信号の推定も可能であり，∑( n H∗ ( f, n)U∗ (n, τ))2 は式 (6) の Pñ∗ ( f ) のよい推定値であると考えられる．

日本音響学会講演論文集

Fig. 2Block diagram of proposed method.

3 提案手法3.1 提案手法の概要
バイノーラル信号では，頭部回折などの影響によりスペクトル上での変形が生じ，左右チャネルでの目的音成分のスペクトル基底は事前学習されたスペクトル基底と異なる．よって，SNMF による妨害音信号の推定をそのまま用いると，事前学習によって得られた教師成分と実際の両耳の目的音成分との差異により，推定精度が劣化する可能性がある．これを解決するためには，例えば観測音信号に合わせて教師基底を変形し分離を行う基底変形型 SNMF を導入することが考えられる．
しかし，従来の基底変形型 SNMF では変形成分が各教師基底の各グリッドに対して求まるため，畳み込みで表されるような時不変の変形に対して過度な変形を許容してしまう．また教師基底の変形と信号分離を同時に行うため，最適化が困難である．そこで
本研究では，SNMF において一般化 MMSE-STSA 推定器の出力値を用いて教師基底の変形を行う手法を提案する．そして，この基底変形手法を従来手法である事前分布パラメータ推定を備えたゲイン共通型
一般化 MMSE-STSA 推定法と組み合わることを考える．Figure 2 に提案手法の概要を示す．
本手法では，観測音信号に対する周波数方向への適応を SNMF における教師基底の変形により行い，時間方向への適応を目的音信号の事前分布パラメータ推定により行う．具体的にはまず，従来手法により推定目的音信号を得る．次に，その出力値を用いて SNMFに用いる教師基底の変形を行う．そして，変形した教師基底を用いて事前分布パラメータ推定を備えたゲイン共通型一般化 MMSE-STSA 推定器を行う．最後にこれらの処理を任意の反復数 itnum だけ行い，その後得られる推定目的音を最終的な出力値とする．
3.2 一般化 MMSE-STSA 推定器の出力を用いた教師基底の適応
本節では，一般化 MMSE-STSA 推定器の出力を用いた教師基底の変形手法について述べる．提案手法
では式 (9) と同じく，目的音信号のみが含まれる教師信号を用いて教師基底 F を学習し，その後，前の反復により得られた推定目的音信号を用いて教師基底の変形成分を学習する．教師基底の変形成分の学習は次式のように行われる．
∑( m)
m m
|I∗m (ω, τ) ◦ s̃m∗ (ω, τ)| ≈ k I∗ (ω, τ) ◦ a∗ F ∗ (ω, k)V∗ (k, τ)(10)
ここで s̃m∗ ( f, τ) は前の反復において得られる一般化
MMSE-STSA 推定器出力の要素値を表す．m, ω はそれぞれ要素値が属するサブバンドのインデックス，サブバンド内における周波数インデックスを表す．
I∗ (ω, τ) は 1 と 0 の値のみを持つバイナリマスクであり，前の反復において式 (2) より計算されるゲイン関数 Gm∗ (ω, τ) が閾値 θ ( 0 ≤ θ ≤ 1 ) を超えているグリッドは 1 を，それ以外は 0 をとるように定める．ま- 550 -2015年3月

た am∗ は教師基底の変形を行うためのスペクトル重みを表し，◦ はアダマール積を表す．この学習では，サブバンド毎にスペクトル重み am∗ を学習し，教師基底に乗ずることで観測音信号中の目的音信号に適応させる．それにより頭部回折や部屋の伝達特性を表すことを期待している．また，推定目的音信号をそのまま用いると，残留している妨害音信号などの影響によりスペクトル重みの学習が適切に実行されない場合がある．そこでバイナリマスク I∗ (ω, τ) を用いることでスペクトル重みの学習に用いるデータの選抜を行っている．これにより，主に目的音信号のみが存在するグリッドを用いて学習を行うことでスペクトル重みの推定精度を向上させることを期待している．

3.3 目的関数と更新式
本節では，スペクトル重みの学習に用いる NMF の目的関数及び更新式の導出過程を示す．本研究では，一般化 KL ダイバージェンスを NMF の目的関数として用いる．一般化 KL ダイバージェンスは次式で定義される．
DKL (y∥x) = y(log y − log x) + x − y(11)これを用いると，スペクトル重みの学習における NMFの目的関数 J は以下の様に表される．
∑ ∑(J = m ω,τ I∗m (ω, τ)DKL | s̃m(ω, τ)| ∥∑ m m)∗aF(ω,k)V(k,τ)∗k ∗ ∗]∑ ∑ m [ mm= m ω,τ iω,τ,∗ −sω,τ,∗ log Rmω,τ,∗ + Rω,τ,∗ + C k(12)∑ m mRm=afv(13)k ∗ ω,k,∗ k,τ,∗ω,τ,∗()mmm
Ck = im(14)
ω,τ,∗ sω,τ,∗ log sω,τ,∗ − sω,τ,∗mmm
ここで，imω,τ,∗ , sω,τ,∗ , fω,k,∗ , vk,τ,∗ はそれぞれ I∗ (ω, τ),mm| s̃∗ (ω, τ)|, F∗ (ω, k), V∗ (k, τ) に対応する要素値である．
式 (12) を最小化する際に更新するものは am∗ 及び vk,τ,∗である．しかしこれらを解析的に導出すことは困難なため，式 (12) の上限を与える関数を用いた補助関数法より更新式を導出する．式 (12) の対数項に関して補助変数 βmk,ω,τ,∗ ≥ 0 と Jensen の不等式を用いるとJ の上限関数 J + は以下の様に与えられる．ただし∑mk βk,ω,τ,∗ = 1 とする．

J ≤ J+[()∑ ∑m ∑mm m= m ω,τ imω,τ,∗ −sω,τ,∗ k βk,ω,τ,∗ log a∗ fω,k,∗ vk,τ,∗]∑m′+ k am(15)∗ fω,k,∗ vk,τ,∗ + C k()′mmm
Ck = Ck + iω,τ,∗ sω,τ,∗ log βk,ω,τ,∗(16)
ここで，式 (15) の等号成立条件は以下で与えられる．
mam∗ fω,k,∗ vk,τ,∗∑βm=k,ω,τ,∗m m
k a∗ fω,k,∗ vk,τ,∗(17)

J + を各変数で偏微分し 0 とおいた式に，補助関数の等号成立条件を代入することで，各変数の乗法更新式を+m導出できる．まず am∗ の更新式を導出する．∂J /∂a∗ =0 より以下が成立する．
∑∑
 m ∑ βmk,ω,τ,∗m
 = 0 (18)
fv
−s+imω,τ,∗  ω,τ,∗ω,k,∗ k,τ,∗ maω,τk∗k

式 (18) に補助変数を代入し整理すると，am∗ の更新式が導出できる．

日本音響学会講演論文集∑mam∗ ← a∗∑ m∑ m mmm−1
ω,τ iω,τ,∗ sω,τ,∗ k ( fω,k,∗ vk,τ,∗ )( k a∗ fω,k,∗ vk,τ,∗ )
∑

ω,τ

(∑m
m
k iω,τ,∗ fω,k,∗ vk,τ,∗)(19)
同様にして vk,τ,∗ の更新式は以下で与えられる
∑
−1mm
m mm∑
mm,ω iω,τ,∗ sω,τ,∗ (a∗ fω,k,∗ )(a∗k fω,k,∗ vk,τ,∗ )()vk,τ,∗ ← vk,τ,∗∑mm∑
mm,τ iω,τ,∗ a∗k fω,k,∗(20)
本手法では，スペクトル重み am∗ の学習は他のサブバンドに依存しない．また，この時に更新されるアク
ティベーション vk,τ,∗ は学習後に破棄する．このようにして学習されたスペクトル重みを教師基底に乗ずる事により，教師基底の変形を行う．その後，変形された教師基底を用いて SNMF を行うことにより妨害音推定を行う．

4 評価実験
4.1 実験条件
提案手法の有用性を確認するために，評価実験を行った．まず，実験用データとして MIDI 音で作成した Oboe，Clarinet，Piano，Cello の 4 種類を用意し，Oboe を抽出対象音として選んだ ( それぞれの楽譜は文献 [5] を参照されたい )．観測音信号に関しては Oboe以外の 3 種類の音源データの中から 1 つを妨害音として選び, 目的音信号，妨害音信号共に −90◦ から 90◦の間で 15◦ の間隔で等パワーにて配置した．各音源の作成にあたっては，ダミーヘッドにより収録された
HRTF[6] を使用し，該当方位の HRTF を畳み込むこと
で両耳信号を作成した．また教師音信号として MIDIシーケンサで Oboe の人工信号を作成し，観測音信号とは別の正面方位の HRTF を畳み込んだものを使用する．この教師音信号は，観測音信号中の各楽器音の音域を含む範囲で，2 オクターブを半音階ずつ上昇する24 音で構成されている．また教師音信号及び，目的音信号はいずれもサンプリング周波数 44.1 kHz であり，FFT 点数 8192, 窓長 4096 点の矩形窓を 512 点でシフトさせてスペクトログラムを作成した．また事前学習における学習基底数は 100, 信号分離における妨害音信号の基底数は 50 とした．忘却係数 α については0.97，振幅圧縮パラメータ β については 1.0 とした．
また，客観評価値として BSS EVAL TOOLBOX[7] による source to distortion ratio (SDR) を用いる． SDRは目的音と妨害音の分離度合いと，一連の信号処理によって生じた目的音の線形ひずみ，非線形ひずみを考慮した推定目的音の品質を表すものであり，値が大きくなるほど精度が良いことを示す．
4.2 実験結果
評価実験では，提案手法と従来手法との比較実験を行う．比較対照法として，反復を行うことなく分離と教師基底の適応を同時に行い妨害音信号を推定する
SNMF と，他のゲイン共通型ポストフィルタとを組み合わせた手法を想定する．従来手法に用いるポストフィルタは，ゲイン共通型 WF (Equi-gain WF) [8], ゲ
イン共通型 MMSE-STSA 推定器 (Equi-gain MMSESTSA) [9], 最小ゲインベースのゲイン共通型 MMSESTSA 推定器 (Gain-min MMSE-STSA) [10], そして，2 章で述べたゲイン共通型一般化 MMSE-STSA 推定器 (Equi-gain Generalized MMSE-STSA) [1] を用いる．これらの従来手法に用いる SNMF については，付録を参照されたい．また，提案手法における反復数は4，閾値 θ は 0.8 とし，全てのサブバンドにおいてス
ペクトル重みの初期値を 1 とした．

- 551 -2015年3月

Equi-gain WF

Equi-gain MMSE-STSA
Equi-gain Generalized
MMSE-STSA

Gain-min MMSE-STSA
Proposed method[6] S. Shimada, K. Sugiyama, H. Hokari “HeadRelated Transfer Function,” Cyber Publishing
Center Press, Cyber Creative Institute, Japan, 2014
(in Japanese).

11SDR [dB]

10[7] E. Vincent, R. Gribonval, C. Fevotte, “Performance measurement in blind audio source separation,” IEEE Trans. Audio, Speech and Language
Processing, vol.14, no.4, pp.1462–1469, 2006.

98765-90-75-60-45-30-150153045607590Angle [deg]Fig. 3

Average SDRs for each method and each direction.

各手法における音源方位別の評価値を Fig. 3 に示す．本図より，ほぼすべての方位で提案手法が従来手法より良い結果を示していることが分かる．この結果
より，教師基底の変形と妨害音の分離を個別に行い，そして学習データの中間的な分離音を用いて教師基底の変形成分を学習することの有効性が確認できた．

5 まとめ
本研究では，SNMF に教師基底の変形成分を反復的に学習する手法を取り入れることにより，従来手法の事前分布パラメータ推定を用いたゲイン共通型一般化 MMSE-STSA 推定法の改良を行った．従来手法の問題点であるスペクトル基底の周波数方向への適応精度が提案手法により向上し，より精度の高い分離ができることを示した．

[8] K. Reindl, Y. Zheng, W. Kellermann, “Speech enhancement for binaural hearing aids based on blindsource separation,” Proc. 2010 International Symposium on Communications, Control and Signal
Processing (ISCCSP2010), 2010.
[9] H. Saruwatari, M. Go, R. Okamoto, K. Shikano,
H. Hosoi, “Binaural hearing aid using soundlocalization-preserved MMSE STSA estimatorwith ICA-based noise estimation,” Proc. 2010 International Workshop on Acoustic Echo and NoiseControl (IWAENC2010), 2010.
[10] A. H. Kamkar-Parsi, M. Bouchard, “Improved
noise power spectrum density estimation for binaural hearing aids operating in a diﬀuse noise field environment, IEEE Trans. Audio, Speech and Lang.
Process., vol.17, no.4, pp.521–533, 2009.

A 付録：基底変形型 SNMF

従来手法では，教師基底の変形と分離を同時に行う基底変形型 SNMF を用いて妨害音推定を行う．目
的関数 J は以下の様に表される．
∑ ∑
∑ m m
(J = m ω,τ DKL Amk a∗ F ∗ (ω, k)V∗ (k, τ)∗ (ω, τ) ∥References∑ ∑ (∑ m))
m+H∗ (ω, n)U∗ (n, τ) + µ m k,n ω F∗ (ω, k)H∗m (ω, n) 2
[1] 室田勇騎，北村大地，小山翔一，猿渡洋，中村哲，(21)“チャネル別事前分布推定と両耳共通スペクトル
ゲインを用いた定位保持型バイノーラル音源分
第二項は教師基底 F∗m (ω, k) とその他の基底 H∗m (ω, n)
離,” 第 29 回信号処理シンポジウム, pp. 486-491,が互いに無相関となる制約を与える直交化罰則項で2014.
あり，µ はその重みパラメータである．各変数の更新[2] D. Kitamura, H. Saruwatari, K. Yagi, K. Shikano,Y. Takahashi, K. Kondo, “Music signal separation
based on supervised nonnegative matrix factorization with orthogonality and maximum-divergencepenalties,” IEICE Trans. Fundamentals of Electronics, Communications and Computer Sciences,
vol.E97-A, no.5, pp.1113–1118, 2014.
[3] D. Kitamura,
H. Saruwatari,
K. Shikano,Y. Takahashi, K. Kondo, ”Music signal separation
by supervised nonnegative matrix factorization with basis deformation,” Proc. DSP2013,
T3P(C)-1, 2013.
[4] C. Breithaupt, M. Krawczyk, R. Martin, “Parameterized MMSE spectral magnitude estimation for
the enhancement of noisy speech,” Proc. IEEE International Conference on Acoustics, Speech andSignal Processing, pp.4037–4040, 2008.
[5] Y. Murota, D. Kitamura, S. Nakai, H. Saruwatari,
S. Nakamura,
K. Shikano,Y. Takahashi,
K. Kondo, “Music signal separation based onbayesian spectral amplitude estimator withautomatic target prior adaptation,” Proc. IEEE
International Conference on Acoustics, Speech andSignal Processing (ICASSP2014), pp.7490–7494,
2014.

日本音響学会講演論文集

式は以下の式で与えられる．
∑∑ m
Am Q mfv
mm ω,τ ω,τ,∗ ω,τ,∗ k ω,k,∗ k,τ,∗)(22)a∗ ← a∗∑ (∑ mω,τ
k fω,k,∗ vk,τ,∗∑m
m
m mm,ω Aω,τ,∗ Qω,τ,∗ (a∗ fω,k,∗ )()vk,τ,∗ ← vk,τ,∗(23)∑m∑mm,τ a∗k fω,k,∗∑mmm
ω,τ Aω,τ,∗ Qω,τ,∗ un,τ,∗mm∑ m ∑ mhω,n,∗ ← hω,n,∗ ∑ mm
τ un,τ,∗ + µ k fω,k,∗ ω′ fω′ ,k,∗ hω′ ,n,∗
(24)
∑mm
m
ω,τ Aω,τ,∗ Qω,τ,∗ hω,n,∗m
∑
um
(25)
n,τ,∗ ← un,τ,∗mm,τ hω,n,∗)−1( ∑∑ mm
m
m
(26)
Qm
ω,τ,∗ = a∗ k fω,k,∗ vk,τ,∗ + n hω,n,∗ un,τ,∗mmm
こ こ で ，Amω,τ,∗ , fω,k,∗ , vk,τ,∗ , hω,n,∗ , un,τ,∗ は そ れ ぞ れmmmA∗ (ω, τ), F∗ (ω, k), V∗ (k, τ), H∗ (ω, n), U∗ (n, τ) に対応する要素値である．

- 552 -

2015年3月