JP6466863B2 - 最適化装置、最適化方法、およびプログラム - Google Patents

Description

本発明は、統計的技術に関し、特に、出力変数を元に、当該出力変数とある程度の相関を持つように入力変数を最適化する技術に関する。

目的音を強調する技術として、ウィナーフィルタに基づく音強調を説明する。時間周波数領域での観測音は以下のように近似できる。
Ｘ_ω，τ＝Ｓ_ω，τ＋Ｎ_ω，τ （１）
ここでω＝｛１，２，・・・，Ω｝とτ＝｛１，２，・・・，Ｆ｝は時間と周波数のインデックス、Ｓ_ω，τは目的音、Ｎ_ω，τは雑音である。ここで音源からマイクロホンまでの伝達特性は問題の簡単のために無視した。さらに、目的音と雑音は無相関であると仮定し、目的音のパワースペクトル密度（ＰＳＤ）をφ_Ｓ，ω,τ＝｜Ｓ_ω,τ｜^２、雑音のＰＳＤをφ_Ｎ，ω,τ＝｜Ｎ_ω,τ｜^２としたとき、目的音を抽出するウィナーフィルタは以下のように近似できる。

ここでξ_ω,τ＝φ_Ｓ，ω,τ／φ_Ｎ，ω,τは事前ＳＮＲを表す。入力されたＸ_ω，τにウィナーフィルタを乗ずることで、目的音Ｙ_ω，τが抽出される（ウィナーフィルタリング）。
Ｙ_ω，τ＝Ｇ_ω，τＸ_ω，τ （３）
式（２）（３）より、雑音下で目的音だけをクリアに収音するためには、目的音と雑音のＰＳＤ φ_Ｓ，ω,τ, φ_Ｎ，ω,τか、事前ＳＮＲξ_ω,τを正確に推定すればよいことが分かる。

従来の雑音下で目的音を強調する技術では、混合ガウスモデル（ＧＭＭ）（例えば、非特許文献１等参照）やディープニューラルネットワーク（ＤＮＮ）（例えば、非特許文献２等参照）などの音響特徴量を用いた音強調技術が代表的である。音源のモデル化に基づくウィナーフィルタ設計は、観測信号からの音響特徴量の抽出と、音響特徴量を事前に学習した統計モデルを用いて事前ＳＮＲ等にマッピングする２つの処理から成る。これらの方法の性能を高めるには、入力された音響特徴量と事前ＳＮＲが強い（非線形な）相関を持つことが必要である。音響特徴量と事前ＳＮＲがいかなる相関も持たない場合、柔軟で洗練されたマッピング法を用いても、収音性能は向上しない。すなわち、音源のモデル化に基づくウィナーフィルタ設計を達成するためには、事前ＳＮＲを正確に推定できる、有効な音響特徴量を選択する必要がある。

Ｄ次元の音響特徴量をｆ_τ＝（ｆ_１,τ，・・・，ｆ_Ｄ,τ）^Ｔ、推定したい事前ＳＮＲをξ_τ、（・）の転置を（・）^Ｔと記述する。ただしξ_τは全周波数ビンの事前ＳＮＲを並べた物でもよいし、フィルタバンクごとの事前ＳＮＲを並べた物でもよいし、ある一つの周波数ビンやフィルタバンクの事前ＳＮＲでもよい。

音響特徴量の選択法として、特徴選択という枠組みがある。これは、大量の音響特徴量の候補の中から目的音の強調に有効な特徴量だけを取り出すものである。ここではＱ（＞Ｄ）次元の音響特徴量の候補ｇ_τの中から、目的音の強調に有効なＤ個の音響特徴量ｆ_τだけを用いて事前ＳＮＲ ξ_τを推定する。特徴選択の手続きは、選択行列Ａ：Ｒ^Ｑ→Ｒ^Ｄ，Ｑ＞Ｄを用いて以下のように表現できる。
ｆ_τ＝Ａｇ_τ （４）
ただし、ｇ_τは音響特徴量の候補を要素とするＱ次元のベクトルであり、ｆ_τは目的音の強調に有効なＤ個の音響特徴量を要素とするＱ次元のベクトルである。選択行列Ａの各行は、１つの要素だけが正の値を持ち、それ以外の要素の値は０となる。つまり音響特徴量の選択問題は、選択行列Ａの最適化問題である。

従来の選択行列Ａの最適化手法の一つに、音響特徴量ｆ_τと事前ＳＮＲ ξ_τとの相互情報量を最大化するように選択行列Ａを最適化するものがある（例えば、非特許文献３等参照）。しかし、この手法で相互情報量を計算するためには、同時分布ｐ（ξ_τ,Ａｇ_τ）や周辺分布ｐ（ξ_τ），ｐ（Ａｇ_τ）が既知である必要がある。多くの場合、これらの分布は未知であり、何らかの形で推定したり近似したりしなくてはならない。非特許文献３では同時分布ｐ（ξ_τ,Ａｇ_τ）をＧＭＭで近似表現し、選択行列Ａと同時分布ｐ（ξ_τ），ｐ（Ａｇ_τ）を一般化ＥＭアルゴリズムで同時最適化しているが、同時分布を十分に近似できず音質が劣化する。

他の従来手法として、再生核ヒルベルト空間上での相互共分散作用素を評価することで相互情報量を計算し、選択行列Ａを最適化する“カーネル次元圧縮”が提案されている（例えば、非特許文献４等参照）。

M. Fujimoto, et al., "Frame-wise model re-estimation method based on Gaussian pruning with weight normalization for noise robust voice activity detection,"Speech communication, vol. 54, pp.229-244, 2012. A. Narayanan, et al., "Ideal ratio mask estimation using deep neural networks for robust speech recognition," in Proc. ICASSP, 2013. 小泉悠馬ほか,"競技音を抽出するための特徴選択と音源強調の統合的アプローチの検討," 音講論（秋）, 2015. K. Fukumizu, et al., "Dimension Reduction for Supervised Learning with Reproducing Kernel Hilbert Space," Journal of Machine Learning Research, vol.5, pp.73-99, 2004.

非特許文献４の手法では、選択行列Ａの最適化が組み合わせ最適化になり、選択行列Ａの各行の設定に、組み合わせ最適化ないしランダムサーチを用いなくてはならない。つまり、全ての組み合わせで相互共分散作用素を評価し、その中で相互共分散作用素が最大となる音響特徴量の組み合わせを選択しなければならない。そのため、音響特徴量の候補の次元が大きくなるにつれ評価が困難になる。

このような問題は、事前ＳＮＲを元に、当該事前ＳＮＲとある程度の相関を持つように音響特徴量を最適化する場合に限られたものではない。何らかの出力変数（出力情報）を元に、当該出力変数（出力情報）とある程度の相関を持つように入力変数を最適化する場合に共通する問題である。

本発明の課題は、出力変数を元に、当該出力変数とある程度の相関を持つように入力変数を最適化する際の演算量を削減することである。

本発明では、インデックスｔでの出力変数ξ_ｔおよびＱ個の入力変数の候補を要素とするベクトルｇ_ｔに対し、出力変数ξ_ｔとベクトルＡｇ_ｔとの相関の高さを表すコスト関数の値が大きくなるように、Ａ^ＴＡの対角成分に対応するベクトルａを更新し、更新されたベクトルａから選択行列Ａの要素を得て出力する。ただし、Ｑ＞Ｄ≧１であり、（・）^Ｔが（・）の転置であり、Ａがベクトルｇ_ｔのＤ個の要素に応じたＤ個の要素からなるベクトルＡｇ_ｔを得るためのＤ行Ｑ列の選択行列であり、ξ_ｔはベクトルｇ_ｔの少なくとも一部の要素と相関を持つ。

以上により、Ａを直接最適化するよりも、出力変数を元に、当該出力変数にある程度の相関を持つように入力変数を最適化する際の演算量を削減できる。

図１は実施形態の最適化装置の全体構成を例示したブロック図である。 図２は実施形態の更新処理部の構成を例示したブロック図である。 図３は実施形態の更新処理部の処理を説明するためのフロー図である。 図４Ａは実施形態の更新量計算部の構成を例示したブロック図である。図４Ｂは実施形態の更新量計算の処理を説明するためのフロー図である。

以下、本発明の実施形態を説明する。
［理論］
まず数学的な理論を説明し、その後で図面を用いて本発明の実施形態を説明する。
ここでは時間周波数領域での観測音の事前ＳＮＲ ξ_ｔ（出力変数）を元に、事前ＳＮＲ ξ_ｔとある程度の相関を持つように音響特徴量（入力変数）を最適化する選択行列Ａを得る場合を説明する。本形態の特徴点は以下の通りである。
（１）選択行列とガウシアンカーネルの特性を利用することにより、組み合わせ最適化を非線形最適化に置き換えた点。
（２）最適化に「確率的最急降下法」を導入し、全学習データを適切なサイズのミニバッチごとに分割してグラム行列の逆行列計算を近似することで、高速に最適化が可能になった点。

まずガウシアンカーネルは以下のように定義される．

ただし、ｋ_ｓ（ξ_τ，ξ_τ’）は各時間インデックスτ，τ’での事前ＳＮＲ ξ_τ，ξ_τ’に対応するガウシアンカーネルを表し、ｋ_ｇ（Ａｇ_τ，Ａｇ_τ’）は各時間インデックスτ，τ’でのＡｇ_τ，Ａｇ_τ’に対応するガウシアンカーネルを表す。ｇ_τ＝（ｇ_１,τ，・・・，ｇ_Ｑ,τ）^Ｔは時間インデックスτでのＱ個の音響特徴量の候補ｇ_ｑ,τ（ただし、ｑ＝１，・・・，Ｑ）を要素としたＱ次元ベクトルであり、ｇ_τ’＝（ｇ_１,τ’，・・・，ｇ_Ｑ,τ’）^Ｔは時間インデックスτ’でのＱ個の音響特徴量の候補ｇ_ｑ,τ’を要素としたＱ次元ベクトルである。ξ_τはベクトルｇ_τの少なくとも一部の要素と相関を持ち、ξ_τ’はベクトルｇ_τ’の少なくとも一部の要素と相関を持つ。ＡはＤ行Ｑ列の選択行列である。選択行列Ａの各行は、１つの要素だけが正の値を持ち、それ以外の要素の値は０となる。ｆ_τ＝Ａｇ_τによって、Ｄ個の音響特徴量ｇ_ｄ,τに対応するＤ個の音響特徴量ｆ_ｄ,τを要素としたＤ次元ベクトルｆ_τ＝（ｆ_１,τ，・・・，ｆ_Ｄ,τ）^Ｔが得られる。また、ｆ_τ’＝Ａｇ_τ’の演算によって、Ｄ個の音響特徴量ｇ_ｄ,τ’に対応するＤ個の音響特徴量ｆ_ｄ,τ’を要素としたＤ次元ベクトルｆ_τ’＝（ｆ_１,τ，・・・，ｆ_Ｄ,τ’）^Ｔが得られる。ＱおよびＤはＱ＞Ｄ≧１を満たす整数であり、例えばＤ≧２である。ｅｘｐ（・）は（・）の指数関数を表し、（・）^Ｔは（・）の転置を表す。

式（５）（６）を用いて計算されるグラム行列は以下となる。

ただし、これらのグラム行列は時間区間［１，・・・，Ｆ］（所定集合）内の各時間インデックスτ＝１，・・・，Ｆ，τ’＝１，・・・，Ｆでのガウシアンカーネルｋ_ｓ（ξ_τ，ξ_τ’）およびｋ_ｇ（Ａｇ_τ，Ａｇ_τ’）に対応するものである。Ｆは１以上の整数であり、例えばＦ≧２である。

カーネルを用いて計算される相互共分散作用素Σ_ｓｓ｜ｇは、グラム行列を用いて以下のように計算できる。
Σ_ｓｓ｜ｇ＝Σ_ｇｇ−Σ_ｓｇΣ_ｇｇ ^−１Σ_ｇｓ （９）
ただし、
Σ_ｓｓ＝Ｋ_ｓＫ_ｓ （１０）
Σ_ｓｇ＝Ｋ_ｓＫ_ｇ （１１）
Σ_ｇｓ＝Ｋ_ｇＫ_ｓ （１２）
Σ_ｇｇ＝Ｋ_ｇＫ_ｇ （１３）
である。Ｋ_ｓおよびＫ_ｇは以下のように計算できる中心化グラム行列である。
Ｋ_ｓ＝ＰＧ_ｓＰ （１４）
Ｋ_ｇ＝ＰＧ_ｇＰ （１５）
ただし、

であり、１_Ｆ＝（１，・・・，１）^Ｔ∈Ｒ^Ｆ（Ｆ次元のベクトル）であり、Ｉ_ＦはＦ×Ｆの単位行列である。

二次モーメントまでで分布系が特定できるということは、再生核ヒルベルト空間で各要素とその条件付分布がガウス分布で表現できることに等しい。ゆえに、ガウス分布のエントロピーの性質から、相互共分散作用素Σ_ｓｓ｜ｇの大きさ（例えば、行列式や負のトレース）を最大化することで、相互情報量を最大化できる。

ところで行列Ａが選択行列の場合、Ａ^ＴＡ∈Ｒ^Ｑ×Ｑは、Ａで選択される音響特徴量の重みに対応する対角成分のみに正の値を持つ特殊な対角行列となる。すると音響特徴量に対応するガウシアンカーネルｋ_ｇ（Ａｇ_τ，Ａｇ_τ’）は以下のように変形できる。

ここでａ_ｑ ^２は行列Ａ^ＴＡのｑ対角要素番目の対角要素である。するとｋ_ｇ（Ａｇ_τ，Ａｇ_τ’）はベクトルａ＝√ｄｉａｇ［Ａ^ＴＡ］に関して微分可能になるため、行列Ａの最適化を行列Ａ^ＴＡの最適化に置き換えれば、非線形最適化問題として解ける。ただし、ｄｉａｇ［Ａ^ＴＡ］は行列Ａ^ＴＡの対角成分を要素とするベクトルを表し、√ｄｉａｇ［Ａ^ＴＡ］は行列Ａ^ＴＡの対角成分の平方根を要素とするベクトルを表す。つまり、相互共分散作用素Σ_ｓｓ｜ｇの大きさの最大化を、選択行列Ａに対してではなくベクトルａに対して行うことで選択行列Ａの最適化が容易になる。相互共分散作用素Σ_ｓｓ｜ｇの最大化に有効な音響特徴量に対応するベクトルａの要素はその絶対値が大きくなり、不要な音響特徴量に対応するベクトルａの要素は０に縮退していく。以降では、相互共分散作用素Σ_ｓｓ｜ｇの大きさをベクトルａに対して最大化する。

相互共分散作用素Σ_ｓｓ｜ｇの大きさは、Σ_ｓｓ｜ｇの行列式や負のトレースで求められるが、ここでは負のトレースをコスト関数（出力変数ξ_ｔとベクトルＡｇ_ｔとの相関の高さを表すコスト関数）として用いた計算方法を説明する。また計算量削減のため、負のトレースを以下のように近似計算する。

ただし、Ｔｒ（・）は（・）のトレースを表す。式（１８）を最大化するための更新式を導出する。式（１８）の最大化は勾配法で行う。勾配法には何を用いてもよいが、更新の収束を速めるために、以下ではAdaDeltaによる実装を説明する。AdaDeltaによるａの更新式は以下となる。

ｓ←γｓ＋（１−γ）ν^２ （２１）
ａ←ａ＋ν （２２）
ただし、式（１９）〜（２２）の更新式におけるベクトルの累乗や除算などの演算は、各要素ごとに行われる。すなわち、式（１９）〜（２２）を要素ごとに書くと以下のようになる。

ｓ_ｑ←γｓ_ｑ＋（１−γ）ν_ｑ ^２ （２５）
ａ_ｑ←ａ_ｑ＋ν_ｑ （２６）
なお、γは０以上１未満の定数であり、εは整数の定数である。「α_１←α_２」はα_２の結果をα_１とする（α_２を新たなα_１とする）ことを意味する。

勾配ベクトル▽ａは以下のように計算できる。

ただしＫ_{τ，τ’,ｑ}は、スペースの関係上、ｋ_ｇ（Ａｇ_τ，Ａｇ_τ’）を単にｋ_ｇ（τ，τ’）と表記し、以下のように表される。

ガウシアンカーネルｋ_ｇ（Ａｇ_τ，Ａｇ_τ’）の偏微分は以下のようになる。

また学習データに対応する時間インデックスの総数（例えば、総フレーム数）がＨの時、Ｆ＝Ｈとした式（２８）中の（Ｋ_ｇ＋εＩ_Ｆ）はＲ^Ｈ×Ｈの対称行列となる。これが学習データの増加により逆行列の計算が困難になる部分である。そこで本形態では、本特許では「確率的最急降下法」にならい、全学習データを適切なサイズのミニバッチにランダム分割し、式（２７）の評価を段階的に行うことで、この問題を回避する。

また、式（２２）の更新では、ベクトルの要素が完全に０になることは稀である。この場合には更新前後での要素の変動が大きくなり、更新が不安定となる場合もある。そこで更新の安定性のために、式（２２）による更新毎に以下のソフトスレッショルディングを行ってもよい。

ただし、βは正則化パラメータ（正値）である。これは、コスト関数（目的関数）にＬ_１正規化項を付与して最適化するのに等しい。

以上のように、相互共分散作用素Σ_ｓｓ｜ｇの大きさを選択行列Ａについて最大化する問題を、行列Ａ^ＴＡの対角成分に対応するベクトルａについて最大化する問題に置き換えることで、選択行列Ａの最適化が容易になる。また、最適化に「確率的最急降下法」を導入し、全学習データを適切なサイズのミニバッチごとに分割してミニバッチごとに逆行列を計算する（グラム行列の逆行列計算を近似することに相当する）ことで演算量を削減し、最適化を高速化できる。なお、ミニバッチサイズは、大きい方が精度はよいものの計算コストが大きくなることから、実装する装置のメモリの大きさや計算能力等を元に、事前に求めておく。

［実施形態］
次に、図面を用いて本形態を詳細に説明する。
＜構成＞
図１に例示するように、本形態の最適化装置１は、記憶部１０１，１０２，１０７，１０９，１１０、周波数領域変換部１０３，１０４、重畳部１０５、事前ＳＮＲ計算部１０８、更新処理部１２０、および出力部１３０を有する。図２に例示するように、更新処理部１２０は、正規化部１２１、初期化部１２２、分割部１２３、更新部１２４、収束判定部１２５、および生成部１２６を有する。更新部１２４は、行列生成部１２４１、更新量計算部１２４２、ベクトル更新部１２４３、およびミニバッチ判定部１２４４を有する。図４Ａに例示するように、更新量計算部１２４２は、更新部１２４２ａ〜１２４２ｃを有する。最適化装置１は、例えば、ＣＰＵ（central processing unit）等のプロセッサ（ハードウェア・プロセッサ）およびＲＡＭ（random-access memory）・ＲＯＭ（read-only memory）等のメモリ等を備える汎用または専用のコンピュータが所定のプログラムを実行することで構成される装置である。このコンピュータは１個のプロセッサやメモリを備えていてもよいし、複数個のプロセッサやメモリを備えていてもよい。このプログラムはコンピュータにインストールされてもよいし、予めＲＯＭ等に記録されていてもよい。また、ＣＰＵのようにプログラムが読み込まれることで機能構成を実現する電子回路（circuitry）ではなく、プログラムを用いることなく処理機能を実現する電子回路を用いて一部またはすべての処理部が構成されてもよい。また、１個の装置を構成する電子回路が複数のＣＰＵを含んでいてもよい。

＜処理＞
次に、本形態の処理を説明する。
≪学習データ≫
目的音の学習データｓ_ｍと雑音の学習データｎ_ｍの時間波形を用意する。ただし、ｍ＝１，・・・，Ｍであり、Ｍは正整数である。ここでサンプリングレートや量子化ビット数は任意であるが、たとえばサンプリングレートを４８ｋＨｚ，量子化ビット数を１６ｂｉｔなどに設定できる。目的音の学習データｓ_ｍは記憶部１０１に格納され、雑音の学習データｎ_ｍは記憶部１０２に格納される（図１）。

≪周波数領域への変換≫
周波数領域変換部１０３，１０４が、ぞれぞれ、記憶部１０１，１０２から読み込んだ目的音と雑音の学習データｓ_ｍ，ｎ_ｍを短時間フーリエ変換（ＳＴＦＴ）などを用いて周波数領域に変換し、目的音の周波数領域信号Ｓ_ω，ｔおよび雑音の周波数領域信号Ｎ_ω，ｔを得て出力する。例えば、フーリエ変換長は１０２４点（サンプリング周波数４８ｋＨｚで約２２ｍｓ），シフト長は５１２点（サンプリング周波数４８ｋＨｚで約１１ｍｓ）などに設定できる。なお、ω＝｛１，２，・・・，Ω｝とｔ＝｛１，２，・・・，Ｆ｝は時間と周波数のインデックスである。ΩおよびＦはそれぞれ正の整数である。

≪重畳≫
重畳部１０５はＳ_ω，ｔおよびＮ_ω，ｔを入力とし、観測信号を模擬的に設計するために、Ｓ_ω，ｔとＮ_ω，ｔを重畳し、時間周波数領域での観測音Ｘ_ω，ｔ＝Ｓ_ω，ｔ＋Ｎ_ω，ｔを得て出力する。

≪音響特徴量候補の抽出≫
音響特徴量候補抽出部１０６は、観測音Ｘ_ω，ｔを入力とし、観測音Ｘ_ω，ｔから時間インデックスｔごとにＱ個の音響特徴量（入力変数）の候補ｇ_ｑ,ｔ（ただし、ｑ＝１，・・・，Ｑ、Ｑ≧２）を抽出し、それらを要素とするＱ次元のベクトルｇ_ｔ＝（ｇ_１,ｔ，・・・，ｇ_Ｑ,ｔ）^Ｔを出力する。候補として用いる音響特徴量は任意であるが、例えば４８次元のメル周波数ケプストラム係数（ＭＦＣＣ）ならびにその一階差分と二階差分、および、４８次元のメルフィルタバンク出力（ＭＦＢＯ）ならびにその一階差分と二階差分などを用いることができる。また、学習データの観測に用いたマイクロホンの個数が複数である場合、ビームフォーミングを行って、方向別にＭＦＣＣやＭＦＢＯを求めることもできる。その他にも、スペクトルフラックスやスペクトルセントロイドなど、Ｑ＝５１２程度の様々な音響特徴量を用いることができる。Ｑ次元のベクトルｇ_ｔは記憶部１０７に格納される。なお、ベクトルｇ_ｔが上記選択行列更新アルゴリズムの入力変数に相当する。

≪事前ＳＮＲの計算≫
事前ＳＮＲ計算部１０８は、Ｓ_ω，ｔおよびＮ_ω，ｔを入力とし、これらから事前ＳＮＲ ξ_ｔ（出力変数）を計算して出力する。例えば、事前ＳＮＲ計算部１０８は、φ_Ｓ，ω,ｔ＝｜Ｓ_ω,ｔ｜^２、φ_Ｎ，ω,ｔ＝｜Ｎ_ω,ｔ｜^２とし、各周波数インデックスωに対応する事前ＳＮＲ ξ_ω,ｔ＝φ_Ｓ，ω,ｔ／φ_Ｎ，ω,ｔからなる列（ξ_１,ｔ，・・・，ξ_Ω,ｔ）を事前ＳＮＲ ξ_ｔとしてもよいし、フィルタバンクごとの事前ＳＮＲを並べたものを事前ＳＮＲ ξ_ｔとしてもよいし、ある一つの周波数インデックスωやフィルタバンクの事前ＳＮＲを事前ＳＮＲ ξ_ｔとしてもよい。ξ_ｔ＝（ξ_１,ｔ，・・・，ξ_Ω,ｔ）の場合、フーリエ変換長が大きいと事前SNRの次元Ωも大きくなるため、演算結果をメルフィルタバンクで圧縮してもよい。メルフィルタバンクの個数はたとえば３２程度に設定できる。事前ＳＮＲ ξ_ｔは記憶部１０９に格納される。なお、事前ＳＮＲ ξ_ｔが上記選択行列更新アルゴリズムの出力変数に相当する。

≪パラメータ≫
以下の定数のパラメータが設定され、記憶部１１０に格納される。
カーネルパラメータ：式（５）のカーネルパラメータσ^２はチューニングして決定すべきであるが、例えば２．０×１０^-２程度に設定できる。
勾配法パラメータ：式（１９）〜（２２）の勾配法パラメータγ，εは例えば、γ＝０．９，ε＝１０^-５に設定できる。
ミニバッチサイズ：ミニバッチサイズＢは学習データの総フレーム数Ｈに応じて変更すべきだが例えばＢ＝２０４８に設定できる。
総フレーム数Ｈ：総フレーム数Ｈは任意であるが、本形態ではＨ＞Ｂである。

≪更新処理≫
更新処理部１２０は、ベクトルｇ_ｔ、事前ＳＮＲ ξ_ｔ、およびパラメータσ^２，γ，ε，Ｂ，Ｈを入力とし、ξ_ｔとＡｇ_ｔとの相関の高さを表すコスト関数の値（関数値、スコア）が大きくなるように、Ａ^ＴＡの対角成分に対応するベクトルａを更新し、更新されたベクトルａから選択行列Ａの要素を得る。

すなわち、本形態の更新処理部１２０は、ξ_ｔを入力として「所定集合」に属する各インデックスτ，τ’でのガウシアンカーネルｋ_ｓ（ξ_τ，ξ_τ’）に対応する中心化グラム行列Ｋ_ｓを得、ｇ_ｔを入力として選択行列Ａを変数としたインデックスτ，τ’でのガウシアンカーネルｋ_ｇ（Ａｇ_τ，Ａｇ_τ’）に対応する中心化グラム行列Ｋ_ｇを得る。さらに、更新処理部１２０は、Σ_ｓｓ＝Ｋ_ｓＫ_ｓ，Σ_ｓｇ＝Ｋ_ｓＫ_ｇ，Σ_ｇｓ＝Ｋ_ｇＫ_ｓ，Σ_ｇｇ＝Ｋ_ｇＫ_ｇとした相互共分散作用素Σ_ｓｓ｜ｇ＝Σ_ｇｇ−Σ_ｓｇΣ_ｇｇ ^−１Σ_ｇｓの大きさが大きくなるように、Ａ^ＴＡの対角成分に対応するベクトルａを更新する。さらに更新処理部１２０は、更新されたベクトルａから選択行列Ａの要素を得て出力する。

特に本形態では、更新処理部１２０は、複数のミニバッチ（部分集合）のそれぞれを「所定集合」として中心化グラム行列Ｋ_ｓおよび中心化グラム行列Ｋ_ｇを得、複数のミニバッチのそれぞれでベクトルａを更新し、複数のミニバッチで更新されたベクトルａから選択行列Ａの要素を得て出力する。これらの処理の詳細は後述する。

≪出力≫
更新処理で得られた選択行列Ａが出力される。任意の装置は、記憶部１０７に格納されたベクトルｇ_ｔと選択行列Ａを用いてｆ_ｔ＝Ａｇ_ｔを計算することで、事前ＳＮＲの推定に有効な音響特徴量を得ることができる。

＜更新処理の詳細＞
図２から図４を用い、更新処理部１２０が行う更新処理の詳細を説明する。
《入力変数の正規化》
まず、正規化部１２１がｇ_ｔ，ξ_ｔ，Ｈを入力とし、以下のようにｇ_ｔ，ξ_ｔを正規化する。

ただし、式（３６）から式（３８）は、ξ_ｔを複数の周波数ビンやフィルタバンクの事前ＳＮＲを並べたものの場合は、各要素について実行する。式（３５）のように更新されたｇ_ｑ，ｔからなる新たなｇ_ｔ＝（ｇ_１，ｔ，・・・，ｇ_Ｑ，ｔ）および新たなξ_ｔは分割部１２３に送られる（ステップＳ１２１）。

《選択行列とAdaDelta更新係数の初期化》
初期化部１２２は、Ｑ次元のベクトルａ＝（ａ_１，…，ａ_Ｑ），ｒ＝（ｒ_１，…，ｒ_Ｑ），ｓ＝（ｓ_１，…，ｓ_Ｑ）を初期化する。初期値は任意だが、例えばａ＝σ^２１_Ｑ，ｒ＝１_Ｑ，ｓ＝０×１_Ｑなどに初期化できる。ただし、１_Ｑ＝（１，・・・，１）^Ｔ∈Ｒ^Ｑである。初期化されたベクトルａ，ｒ，ｓは分割部１２３に送られる（ステップＳ１２２）。

《入出力変数のミニバッチ分割》
分割部１２３は、正規化部１２１から出力された入力変数である（ｇ_１，・・・，ｇ_Ｈ）および出力変数である（ξ_１，・・・，ξ_Ｈ）をランダムにＢ個ずつのサブセットに分割する。各サブセットをミニバッチと呼ぶ。１つのミニバッチは連続する複数の時間インデックスに対応していてもよいし、隣り合わない時間インデックスに対応していてもよい。すなわち、処理対象の時間インデックスの区間［１，・・・，Ｈ］（処理対象区間）が複数の部分集合であるミニバッチに区分されればよい。例えば分割部１２３は、（ｇ_１，・・・，ｇ_Ｈ）をミニバッチ（ｇ_１，・・・，ｇ_Ｂ），（ｇ_Ｂ＋１，・・・，ｇ_２Ｂ），・・・，（ｇ_{Ｈ−Ｂ＋１}，・・・，ｇ_Ｈ）に分割し、（ξ_１，・・・，ξ_Ｈ）をミニバッチ（ξ_１，・・・，ξ_Ｂ），（ξ_Ｂ＋１，・・・，ξ_２Ｂ），・・・，（ξ_{Ｈ−Ｂ＋１}，・・・，ξ_Ｈ）に分割する（ステップＳ１２３）。

《中心化グラム行列の計算》
行列生成部１２４１は、（ｇ_１，・・・，ｇ_Ｈ）の１つのミニバッチ、（ξ_１，・・・，ξ_Ｈ）の１つのミニバッチ、およびσ^２が入力される。入力される（ｇ_１，・・・，ｇ_Ｈ）のミニバッチと（ξ_１，・・・，ξ_Ｈ）のミニバッチとは同じ時間インデックスに対応する。説明の便宜上、ステップＳ１２４１〜Ｓ１２４４では、処理対象として入力された、（ｇ_１，・・・，ｇ_Ｈ）のミニバッチを（ｇ_１，・・・，ｇ_Ｆ）と読み替え、（ξ_１，・・・，ξ_Ｈ）のミニバッチを（ξ_１，・・・，ξ_Ｆ）と読み替えて説明する。

行列生成部１２４１は、ミニバッチ（ξ_１，・・・，ξ_Ｆ）および（ｇ_１，・・・，ｇ_Ｆ）について、前述した式（１４）および（１５）に従って、選択行列Ａを変数とした中心化グラム行列Ｋ_ｓおよびＫ_ｇを得て出力する。すなわち、行列生成部１２４１は、ξ_ｔを入力としてミニバッチ（所定集合）の各インデックスτ＝１，・・・，Ｆ，τ’＝１，・・・，Ｆでのガウシアンカーネルｋ_ｓ（ξ_τ，ξ_τ’）に対応する中心化グラム行列Ｋ_ｓを得、ベクトルｇ_ｔを入力として選択行列Ａを変数としたインデックスτ＝１，・・・，Ｆ，τ’＝１，・・・，Ｆでのガウシアンカーネルｋ_ｇ（Ａｇ_τ，Ａｇ_τ’）に対応する中心化グラム行列Ｋ_ｇを得て出力する（ステップＳ１２４２）。

《更新量の計算》
更新量計算部１２４２は、γ，ε，Ｋ_ｓ，Ｋ_ｇを入力とし、前述の式（１９）〜（２１）（すなわち、式（２３）〜（２５））を計算してベクトルνおよびベクトルｓを更新する（ステップＳ１２４２）。ベクトルνはベクトルａの更新量を表す。すなわち、更新量計算部１２４２の更新部１２４２ａが式（１９）（すなわち、式（２３））に従ってベクトルｒを更新し（ステップＳ１２４２ａ）、更新部１２４２ｂが式（２０）（すなわち、式（２４））に従ってベクトルνを更新し（ステップＳ１２４２ｂ）、更新部１２４２ｃが式（２１）（すなわち、式（２５））に従ってベクトルｓを更新する（ステップＳ１２４２ｃ）。更新されたｒ，ｓは図示していないメモリに格納される。

《ベクトルａの更新》
ベクトル更新部１２４３は、ステップＳ１２４２ｂで更新されたベクトルνを用い、式（２２）（すなわち、式（２６））に従ってベクトルａを更新する。あるいは、更新が不安定となる場合、式（２２）に代えて以下のように慣性項付きの更新を行ってもよい。
ａ←ａ＋｛ζν＋（１−ζ）ν’｝ （３９）
ただし、ζは０＜ζ＜１の慣性項であり、例えばζ＝０．１に設定できる。またν’は、ステップＳ１２４２の更新前（１回の更新前）のνである。また、ａはＡ^ＴＡの対角項の各値の平方根を要素としたベクトルであり、その要素は０以上であるべきである。しかし、上記の方法で更新されたベクトルνの値によっては更新されたベクトルａが負の要素を持つ場合もある。そのため、ベクトルａ＝（ａ_１，・・・，ａ_Ｑ）の更新のたびに、その各要素ａ_ｑ（ただし、ｑ＝１，・・・，Ｑ）に対して以下の処理を行ってもよい。
ａ_ｑ＝ｍａｘ（ａ_ｑ，０） （４０）
ただし、ｍａｘ（ａ_ｑ，０）は、ａ_ｑと０とのうち大きい方の値を意味する。ただし、ａ_ｑ＝０の場合には、ｍａｘ（ａ_ｑ，０）＝０とする。あるいは、更新の安定のために前述の式（３２）によるソフトスレッショルディングを行ってもよい。すなわち、ベクトル更新部１２４３は、ベクトルａと、ベクトルνに応じたＱ次元ベクトルとを加算したＱ次元ベクトルに応じた新たなベクトルａを得ればよい。更新されたａは図示していないメモリに格納される（ステップＳ１２４３）。

《ミニバッチ判定》
ミニバッチ判定部１２４４は、ステップＳ１２３で分割されたすべてのミニバッチについてステップＳ１２４１〜Ｓ１２４３の処理を終了したかを判断する（ステップＳ１２４４）。ステップＳ１２４１〜Ｓ１２４３の処理を実行していないミニバッチが存在する場合にはステップＳ１２４１に戻り、これらの処理が実行されてないミニバッチについて処理を行う。一方、すべてのミニバッチについてステップＳ１２４１〜Ｓ１２４３の処理をしていた場合、ステップＳ１２５の収束判定に進む。

《収束判定》
収束判定部２１５は、終了条件を満たしたか否かを判定する（ステップＳ１２５）。例えば、収束判定部２１５は、ステップＳ１２３，Ｓ１２４（Ｓ１２４１〜Ｓ１２４４）の処理の繰り返し回数が一定数以上となった場合に終了条件を満たしたと判断してもよいし、全ミニバッチに対する更新処理（ステップＳ１２４）前後でのベクトルａの変化量のノルムが一定値以下となった場合に終了条件を満たしたと判断してもよい。終了条件を満たしていないと判断した場合にはステップＳ１２３に戻って処理をやり直す。終了条件を満たしたと判断した場合には、ステップＳ１２６の選択行列の生成に進む。

《選択行列の生成》
生成部１２６は、ステップＳ１２４で更新されたベクトルａを入力とし、このベクトルａから選択行列Ａの要素を得て出力する。生成部１２６は、例えば、ベクトルａ＝（ａ_１，・・・，ａ_Ｑ）の要素ａ_ｑのうち閾値α_ｔｈ以上の要素数をＤ’とし、閾値α_ｔｈ以上の要素の次元をｑ_ｄ（ただし、ｄ＝１，・・・，Ｄ’）とし、選択行列Ａのｄ行ｑ_ｄ列の要素Ａ（ｄ，ｑ_ｄ）を正値（例えば１）とし、他の要素を０とする。なお、閾値α_ｔｈは正値であり、Ｄ’が所定値以上（例えば、１以上）となるように設定されてもよいし、予め定められていてもよい。具体的には以下の手順によって選択行列Ａを生成できる。
１．生成部１２６は、Ａを要素がすべて０のＤ’×Ｑ行列として初期化する。
２．生成部１２６は、ｄ＝１，・・・，Ｄ’に対し，以下の処理を実行する。
３．生成部１２６は、ａのｄ番目の閾値α_ｔｈ以上の要素の次元をｑ_ｄとして、ｄ行ｑ_ｄ列の要素をＡ（ｄ，ｑ_ｄ）＝１とする。

有効な要素数Ｄを固定したい場合、生成部１２６は、例えば、ベクトルａ＝（ａ_１，・・・，ａ_Ｑ）の要素ａ_ｑのうち大きいものからＤ個を選び、それらの次元をｑ_ｄ（ただし、ｄ＝１，・・・，Ｄ）とし、選択行列Ａのｄ行ｑ_ｄ列の要素Ａ（ｄ，ｑ_ｄ）を正値（例えば１）とし、他の要素を０としてもよい。具体的には以下の手順によって選択行列Ａを生成できる。
１．生成部１２６は、Ａを要素がすべて０のＤ×Ｑ行列として初期化する。
２．生成部１２６は、ｄ＝１からＤに対し、以下の処理を実行する。
３．生成部１２６は、ａの要素を降順に並び替え、ａのｄ番目の要素に対応する次元をｑ_ｄとして、ｄ行ｑ_ｄ列の要素をＡ（ｄ，ｑ_ｄ）＝１とする。

＜本形態の特徴＞
以上のように、本形態では、相互共分散作用素の大きさ（コスト関数の大きさ）を最大化するＡを探索する問題を、相互共分散作用素の大きさを最大化する行列Ａ^ＴＡの対角成分に対応するベクトルａを探索する問題に置き換えた。これによって偏微分が可能となり、非線形最適化問題として解を求めることができる。その結果、従来よりも演算量を削減することができる。

すなわち、相互共分散作用素の大きさを最大化するＡを探索する問題、すなわち、選択行列Ａの各行でどの音響特徴量をアクティブにするかの探索は、組み合わせ最適化ないしランダムサーチを用いなくてはならなかった。つまり、全てのパターンで相互共分散作用素の大きさを評価し、その中で相互共分散作用素が最大となる音響特徴量の組み合わせを選択しなければならなかった。そのため、音響特徴量の候補の次元が大きくなるにつれ、評価が困難になった。一般に音響特徴量は、複数の変数の組み合わせで音の性質を説明する。例えば音声認識で広く用いられるＭＦＣＣは、スペクトル包絡を２４個程度の変数で説明する。また、音量の計算で用いられるフィルタバンク分析は、全帯域の周波数を３２程度のフィルタで分析する。突発性を調べるためには、その一階差分（Δ特徴量）や二階差分（ΔΔ特徴量）を計算する。このように、音響特徴量の候補は数多くあり、その次元数は非常に大きい。この大量の音響特徴量の候補から最適な組み合わせを選ぶとき、たとえばＱ＝５１２次元の候補から，Ｄ＝４８次元の音響特徴量を選択する場合、その組み合わせ数は５１２Ｃ４８≒９．２×１０^６７という膨大な数になり、その探索は事実上困難である。これに対し、本形態の手法では非線形最適化問題として解を求めることができるため、探索に必要な演算量を大幅に削減できる。

さらに、本形態ではミニバッチごとに段階的にベクトルａを更新するため、式（２８）中の（Ｋ_ｇ＋εＩ_Ｆ）の行列サイズを小さくでき、その逆行列（Ｋ_ｇ＋εＩ_Ｆ）^−１の演算量を小さくできる。

すなわち、一般的なカーネル法を用いて相互情報量を計算するとき、計算途中でグラム行列呼ばれる行列の逆行列の評価が必要となる。学習データの総フレーム数がθの時、グラム行列はＲ_θ×θの対称行列となる。音データを、例えば分析窓３２ｍｓ，シフト幅１６ｍｓで分析するとき、フレーム数は１秒間で６０フレームとなる。音データを用いた機械学習の学習データ量は数十時間となることが多く、総フレーム数θは簡単に１万を越えてしまう．その逆行列の計算は事実上困難であり、音データの学習にカーネル法を用いることは困難であった。これに対し、本形態ではミニバッチごとに段階的にベクトルａを更新するため、総フレーム数が増加した場合でも演算量の増加を大幅に抑制できる。

［変形例等］
なお、本発明は上述の実施形態に限定されるものではない。本発明は、例えば、音響信号強調フィルタ（出力変数）を元に、複数の音響特徴量の候補（入力変数）の中から音響信号強調フィルタとの相関が小さな音響特徴量を除外するための選択行列を求めるために利用できる。しかしながら、本発明の用途はこれに限定されない。すなわち、本発明は、出力変数を元に、当該出力変数とある程度の相関を持つように入力変数を最適化する用途、例えば、相互情報量最大化に基づき、ある入力変数からそれに対応する出力変数を推定する際に、出力変数を推定するために有効な情報を残すように、入力変数のサブセットを選択する用途であれば、どのような用途にも適用できる。そのため、この使用用途は音源強調や事前ＳＮＲの推定に限らない。つまり入力変数は音響特徴量（音響信号の特徴量）に限らず、画像やセンサなどのなんらかのデータから特徴抽出した結果（画像信号の特徴量やセンサ信号の特徴量）であってもよいし、音響信号、画像信号、センサ信号、位置座標などの生データであってもよい。すなわち、入力変数が、音響信号、画像信号、センサ信号、位置座標、その他の時系列データ（例えば、生データ）や、それらの特徴量を含んでもよい。同様に出力変数も事前ＳＮＲに限らず、画像の属するクラスを表す変数でもよいし、元信号の振幅スペクトル｜Ｓ_ω，τ｜であってもよい。すなわち、出力変数が音響信号、画像信号、またはセンサ信号、位置座標、その他の時系列データの何れかに対応する情報を含んでもよい。その他、入力変数が時系列データ以外のデータまたはその特徴量を含んでもよく、出力変数が時系列データ以外のデータに対応する情報を含んでもよい。すなわち、入力変数および出力変数が時系列の情報でなくてもよく、本発明は時系列信号以外にも適用可能である。

例えば、以下のようにして、画像信号に対する判別を行うことができる。この場合の最適化装置は、記憶部１０１，１０２、周波数領域変換部１０３，１０４、重畳部１０５、事前ＳＮＲ計算部１０８、音響特徴量候補抽出部１０６に代えて、学習データである画像信号を格納する記憶部、および、画像特徴量の候補を抽出する画像特徴量候補抽出部を含む。学習データである画像信号は、判別対象の画像に対応するものと、それ以外の画像に対応するものを含む時系列信号である。時間インデックスｔでの画像信号が判別対象の画像である場合、その画像信号に出力変数ξ_ｔ＝１が対応付けられている。一方、時間インデックスｔでの画像信号が判別対象ではない画像である場合、その画像信号に出力変数ξ_ｔ＝０が対応付けられている。画像特徴量候補抽出部は、入力された学習データである画像信号から、時間インデックスｔごとにＱ個の画像特徴量の候補を抽出し、それらを要素とするＱ次元のベクトルをｇ_ｔとして記憶部１０７に格納する。画像特徴量の候補は、例えば、コーナー検出やＳＨＩＦＴ（Scale-Invariant Feature Transform）特徴などを用いて、入力された学習データに含まれる画像信号から抽出した特徴量である。また、各出力変数ξ_ｔは記憶部１０９に格納される。それ以外は上述の実施形態で説明した通りである。

また上記の実施形態では、コスト関数として相互共分散作用素の負のトレースを用いたが、コスト関数として相互共分散作用素の行列式を用いてもよい。その他、入力変数の候補と出力変数との相関の強さを示すコスト関数であればどのようなものを用いてもよい。また、上記の実施形態では、相互共分散作用素の大きさを最大化するベクトルａを選択した。しかしながら、相互共分散作用素の大きさを大きくするようにベクトルａを更新して得られる値であれば、最終的に選択されるベクトルａが相互共分散作用素の大きさを最大化しなくてもよい。その他のコスト関数を用いる場合も同様である。

上記の実施形態では、ａ＝√ｄｉａｇ［Ａ^ＴＡ］としたが、Ａ^ＴＡの対角成分に対応するベクトルであればどのようなものをａとしてもよい。例えば、Ａ^ＴＡの対角成分そのものをベクトルａの要素にしてもよいし、Ａ^ＴＡの対角成分やその平方根の関数値をベクトルａの要素にしてもよい。

また、学習データに対応する時間インデックスの総数Ｈが小さい場合には、ミニバッチごとではなく、すべての時間インデックスｔ＝１，・・・，Ｈについてまとめて更新処理を行ってもよい。この場合には、分割部１２３、ミニバッチ判定部１２４４、ステップＳ１２３、Ｓ１２４４は不要であり、ステップＳ１２２の後、Ｆ＝Ｈとして、Ｓ１２４１〜Ｓ１２４３の処理を実行し、その後ステップＳ１２５の処理を行えばよい。

また実施形態ではインデックスｔがフレーム番号などの時間インデックスである場合を例示した。しかしながら、インデックスｔがフレーム番号以外の時間インデックスであってもよいし、周波数ビンに対応する周波数インデックスであってもよいし、その他のインデックスであってもよい。

上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。

上述の構成をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。このプログラムをコンピュータで実行することにより、上記処理機能がコンピュータ上で実現される。この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体の例は、非一時的な（non-transitory）記録媒体である。このような記録媒体の例は、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等である。

このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。

このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。処理の実行時、このコンピュータは、自己の記憶装置に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。

上記実施形態では、コンピュータ上で所定のプログラムを実行させて本装置の処理機能が実現されたが、これらの処理機能の少なくとも一部がハードウェアで実現されてもよい。

１ 最適化装置
１２０ 更新処理部

Claims (6)

1. Ｑ＞Ｄ≧１であり、Ｆが正整数であり、所定区間［１，…，Ｆ］に属するインデックスτ，τ’，ｔがτ＝１，・・・，Ｆ，τ’＝１，・・・，Ｆ，ｔ＝１，・・・，Ｆであり、（・）^Ｔが（・）の転置であり、ｇ_ｔがインデックスｔでのＱ個の入力変数の候補を要素とするベクトルであり、Ａが前記ベクトルｇ_ｔのＤ個の要素に応じたＤ個の要素からなるベクトルＡｇ_ｔを得るための選択行列であり、ξ_ｔが前記ベクトルｇ_ｔの少なくとも一部の要素と相関を持つ出力変数であり、
   前記出力変数ξ_ｔを入力として前記所定区間［１，…，Ｆ］に属する各インデックスτ＝１，・・・，Ｆ，τ’＝１，・・・，Ｆでのガウシアンカーネルｋ_ｓ（ξ_τ，ξ_τ’）に対応する中心化グラム行列Ｋ_ｓを得、前記ベクトルｇ_ｔを入力として前記選択行列Ａを変数とした前記インデックスτ＝１，・・・，Ｆ，τ’＝１，・・・，Ｆでのガウシアンカーネルｋ_ｇ（Ａｇ_τ，Ａｇ_τ’）に対応する中心化グラム行列Ｋ_ｇを得る行列生成部と、
   Σ_ｓｓ＝Ｋ_ｓＫ_ｓ，Σ_ｓｇ＝Ｋ_ｓＫ_ｇ，Σ_ｇｓ＝Ｋ_ｇＫ_ｓ，Σ_ｇｇ＝Ｋ_ｇＫ_ｇとした相互共分散作用素Σ_ｓｓ｜ｇ＝Σ_ｇｇ−Σ_ｓｇΣ_ｇｇ ^−１Σ_ｇｓの大きさが大きくなるように、Ａ^ＴＡの対角成分に対応するベクトルａを更新するベクトル更新部と、
   更新された前記ベクトルａから前記選択行列Ａの要素を得て出力する生成部と、
   を有し、
   

   

   であり、Ｉ _Ｆ がＦ×Ｆの単位行列であり、１ _Ｆ がＦ次元のベクトル（１，…，１） ^Ｔ であり、Ｐ＝Ｉ _Ｆ −（１／Ｆ）１ _Ｆ １ _Ｆ ^Ｔ であり、Ｋ _ｓ ＝ＰＧ _ｓ Ｐであり、Ｋ _ｇ ＝ＰＧ _ｇ Ｐである、最適化装置。
2. 請求項１の最適化装置であって、
   ＨがＨ＞Ｆを満たす正整数であり、前記出力変数ξ _ｔ の集合（ξ _１ ，・・・，ξ _Ｈ ）および前記ベクトルｇ _ｔ の集合（ｇ _１ ，・・・，ｇ _Ｈ ）が学習データであり、
   前記行列生成部は、前記学習データを構成する部分集合（ξ _１ ，・・・，ξ _Ｆ ）および（ｇ _１ ，・・・，ｇ _Ｆ ）ごとに前記中心化グラム行列Ｋ_ｓおよび前記中心化グラム行列Ｋ_ｇを得、
   前記ベクトル更新部は、前記部分集合（ξ _１ ，・・・，ξ _Ｆ ）および（ｇ _１ ，・・・，ｇ _Ｆ ）ごとに前記ベクトルａを更新し、
   前記生成部は、前記部分集合（ξ _１ ，・・・，ξ _Ｆ ）および（ｇ _１ ，・・・，ｇ _Ｆ ）ごとに更新された前記ベクトルａから前記選択行列Ａの要素を得て出力する最適化装置。
3. 請求項１または２の最適化装置であって、
   γが０以上１未満の定数であり、εが定数であり、Ｊ＝−Ｔｒ｛Ｋ_ｓ（Ｋ_ｇ＋εＩ_Ｆ）^−１｝であり、Ｔｒ（・）が（・）のトレースであり、
   Ｑ個の要素からなるベクトルｒ＝（ｒ_１，…，ｒ_Ｑ），ｓ＝（ｓ_１，…，ｓ_Ｑ）および前記ベクトルａ＝（ａ_１，…，ａ_Ｑ）を初期化する初期化部と、
   

   

   をｑ番目（ただし、ｑ＝１，…，Ｑ）の要素とするＱ次元のベクトルを新たな前記ベクトルｒとする第１更新部と、
   

   

   をｑ番目の要素とするＱ次元のベクトルを新たなベクトルνとする第２ベクトル更新部と、
   γｓ_ｑ＋（１−γ）ν_ｑ ^２をｑ番目の要素とするＱ次元のベクトルを新たな前記ベクトルｓとする第３更新部と、をさらに有し、
   前記ベクトル更新部は、前記ベクトルａと、前記ベクトルνに応じたＱ次元ベクトルと、を加算したＱ次元ベクトルに応じた新たな前記ベクトルａを得る第４ベクトル更新部と、
   を含む最適化装置。
4. 請求項１から３の何れかの最適化装置であって、
   前記入力変数が、音響信号、画像信号、センサ信号、音響信号の特徴量、画像信号の特徴量、またはセンサ信号の特徴量を含み、
   前記出力変数が、音響信号、画像信号、またはセンサ信号の何れかに対応する情報を含む、最適化装置。
5. Ｑ＞Ｄ≧１であり、Ｆが正整数であり、所定区間［１，…，Ｆ］に属するインデックスτ，τ’，ｔがτ＝１，・・・，Ｆ，τ’＝１，・・・，Ｆ，ｔ＝１，・・・，Ｆであり、（・）^Ｔが（・）の転置であり、ｇ_ｔがインデックスｔでのＱ個の入力変数の候補を要素とするベクトルであり、Ａが前記ベクトルｇ_ｔのＤ個の要素に応じたＤ個の要素からなるベクトルＡｇ_ｔを得るための選択行列であり、ξ_ｔが前記ベクトルｇ_ｔの少なくとも一部の要素と相関を持つ出力変数であり、
   行列生成部が、前記出力変数ξ_ｔを入力として前記所定区間［１，…，Ｆ］に属する各インデックスτ＝１，・・・，Ｆ，τ’＝１，・・・，Ｆでのガウシアンカーネルｋ_ｓ（ξ_τ，ξ_τ’）に対応する中心化グラム行列Ｋ_ｓを得、前記ベクトルｇ_ｔを入力として前記選択行列Ａを変数とした前記インデックスτ＝１，・・・，Ｆ，τ’＝１，・・・，Ｆでのガウシアンカーネルｋ_ｇ（Ａｇ_τ，Ａｇ_τ’）に対応する中心化グラム行列Ｋ_ｇを得る行列生成ステップと、
   ベクトル更新部が、Σ_ｓｓ＝Ｋ_ｓＫ_ｓ，Σ_ｓｇ＝Ｋ_ｓＫ_ｇ，Σ_ｇｓ＝Ｋ_ｇＫ_ｓ，Σ_ｇｇ＝Ｋ_ｇＫ_ｇとした相互共分散作用素Σ_ｓｓ｜ｇ＝Σ_ｇｇ−Σ_ｓｇΣ_ｇｇ ^−１Σ_ｇｓの大きさが大きくなるように、Ａ^ＴＡの対角成分に対応するベクトルａを更新するベクトル更新ステップと、
   生成部が、更新された前記ベクトルａから前記選択行列Ａの要素を得て出力する生成ステップと、
   を有し、
   

   

   であり、Ｉ _Ｆ がＦ×Ｆの単位行列であり、１ _Ｆ がＦ次元のベクトル（１，…，１） ^Ｔ であり、Ｐ＝Ｉ _Ｆ −（１／Ｆ）１ _Ｆ １ _Ｆ ^Ｔ であり、Ｋ _ｓ ＝ＰＧ _ｓ Ｐであり、Ｋ _ｇ ＝ＰＧ _ｇ Ｐである、最適化方法。
6. 請求項１から４の何れかの最適化装置としてコンピュータを機能させるためのプログラム。

Images