JP5323860B2 - 混合された信号を複数の成分信号に分離する方法 - Google Patents

Description

本発明は、信号処理の一般的な技術分野に関するものであり、より具体的には、ブラインドソース分離を用いてソース信号の混合体からソース信号を分離する一般的な技術分野に関するものである。

本発明は、数多くの分野に応用できる。特に、脳波図、無線通信の分野における信号処理、さらには質量スペクトル処理にこれを応用することができる。

いずれも所与の瞬間にセンサーネットワークによって記録された複数の混合体から統計的に独立した複数のソースを分離すること（空間分散）からなるソース分離の問題は、信号処理において比較的新しい問題である。

観察された物理的現象が、複数の独立した基本的現象の積み重ねであり、この積み重ねが瞬間線形結合（または混合体）によって表現されることは、頻繁にあることである。

この現象の観察に由来する測定データは、このとき、各々の独立した基本的現象からの寄与の瞬間線形混合体である。

この場合、複数のセンサーを設置することによって、これらの寄与の複数の異なる瞬間線形混合体を測定することが可能である。例えば、脳波検査においては、患者の頭皮上にさまざまな電極が設置される。これらの電極によって、患者の脳活動を測定することができる。

各電極は、脳のさまざまなソースが発出する電気信号の混合体を測定し、異なる電極が異なる混合体を測定する。ソースの混合体は、ソースがその中で異なる割合でまたは異なる重みをもって出現する場合に、異なるものであると考えられている。
この場合、次のものについてのその他の情報、すなわち、
−ソースについて、および、
−これらのソース由来の信号がどのように混合されているか
は全く知られていない。

このとき、第一段階において混合体の重みを識別するために、そして第二段階において測定値内で各ソースの寄与を分離するために、ソース分離が使用される。こうして、混合体に妨害されることなく異なる寄与の個別分析が可能となる。

現在、ソース分離方法には二つのタイプが存在する。

第一のタイプの方法は、予備ステップ、いわゆる「白色化（ｗｈｉｔｅｎｉｎｇ）」ステップを使用し、これにより、混合体の重みを識別するために後続する計算ステップを容易にすることができる。しかしながら、この白色化ステップは、各ソースの寄与の分離時点で補正不能なエラーを誘発する。したがって、このタイプの方法はあまり正確ではない。

第二のタイプの方法は、白色化ステップを一切必要としない。しかしながら、この第二のタイプの方法は、効率および計算コストに関する問題を誘発する。

本発明の一つの目的は、混合された信号を複数の成分信号に分離するための方法を提案することにあり、これにより前述のタイプの方法の欠点を少なくとも一つ克服することができる。

この目的で、少なくとも二つのセンサーにより記録された混合された信号（ｘ_１（ｔ），．．．，ｘ_Ｎ（ｔ））を複数の成分信号（ｓ_１（ｔ），．．．，ｓ_Ｍ（ｔ））へと分離するための方法が提供されるものであり、該方法においては、混合された信号（ｘ_１（ｔ），．．．，ｘ_ｎ（ｔ））は独立したソースに由来するＭ個の成分信号（ｓ_１（ｔ），．．．，ｓ_Ｍ（ｔ））のＮ個の線形結合に対応しており、ソースに由来する成分信号（ｓ_１（ｔ）、．．．、ｓ_Ｍ（ｔ））のベクトルを混合行列（Ａ）に乗じた積として各混合された信号を記すことが可能であり、該方法はさらに結合対角化ステップを含み、このステップでは一つまたは複数のギブンス回転行列（Ｇｉｖｅｎｓ ｒｏｔａｔｉｏｎ ｍａｔｒｉｃｅｓ）（Ｇ（θ））、および、一つまたは複数の「ハイパボリック（双曲線）回転」行列（ｈｙｐｅｒｂｏｌｉｃ ｒｏｔａｔｉｏｎ ｍａｔｒｉｃｅｓ）（Ｈ（φ））に対応する複数の回転行列の積から混合行列Ａが推定される。

本発明の範囲内で「ハイパボリック回転行列」というのは、ハイパボリック関数（例えばｃｏｓｈ、ｓｉｎｈ）である項を含み、かつ１に等しい行列式を有する回転行列を意味する。

かくして、本発明によると、高速で信頼性の高い単一の最適化ステップにより混合行列を識別できる。単一のステップしか存在しないという事実によって、後述の文献［１］、［１ｂ］、［２］の場合などの２ステップの方法（白色化と、それに続く結合対角化）に比べてさらに正確な結果が得られる。実際、第一のステップの不可避なエラーを第二のステップにより補正することはできない。さらに、提案される方法では、特にロバストな、そして計算能力の観点からみてさほどコストが高くない代数アルゴリズムが使用される（これは、２×２サイズの「三角法」および「ハイパボリック」回転行列の積である）。

白色化ステップを備える方法は、基本直交行列（ｅｌｅｍｅｎｔａｒｙ ｏｒｔｈｏｇｏｎａｌ ｍａｔｒｉｃｅｓ）の積により直交混合行列（ｏｒｔｈｏｇｏｎａｌ ｍｉｘｔｕｒｅ ｍａｔｒｉｘ）を確立する（ギブンス回転）。

本発明に関して言うと、これは１に等しい行列式を有する行列を乗算することにより、１に等しい行列式を有する混合行列を構築することを可能にする。直交行列群（正の行列式＋１を有する）は、１に等しい行列式を有する行列群の中に含み入れられる。したがって、混合行列を計算するステップは、混合体が当初受けていた制約を制限することのできる本発明の方法によって容易になる。最初の「白色化」ステップはもはや不可欠ではなく、それでも「代数」法のロバスト性および効率は保たれる。

結果として、本発明の方法は、白色化ステップを含む方法で得られた混合行列よりも混合体の物理的現実に近い混合行列の計算を可能にする。

本発明に係る方法の好ましいものの非限定的な態様は、以下の通りである。

−ギブンス回転行列が、θの三角関数に対応する要素を含み、ハイパボリック回転行列がφのハイパボリック関数に対応する要素を含む；

−対角化ステップには、右にはギブンス回転行列とハイパボリック回転行列の積、そして左にはこの積の転置行列により推定行列族（Ｒ：_ｋ）を乗算するためのステップ、そして、この乗算の結果としてのｎ次の行列の非対角項を最小化するためにθ_ｎおよびφ_ｎの値を計算するためのステップとが含まれている；
なお、本明細書においては、明細書の記載として使用できる文字の制約を考慮し、

上記〔数１〕を、テキスト上、「（Ｒ：_ｋ）」と表記する。

−結果として得られたｎ次の行列に対して乗算ステップが反復される；

−ｃｏｓθおよびｃｏｓｈφがそれぞれ１に向かうギブンス角度θおよび「ハイパボリック角度」φになるまで乗算ステップが反復される；

−推定行列Ｒ：_ｋは、（推定エラーの範囲内で）混合行列と対角行列と転置された混合行列との積に等しく（Ｒ_ｋ＝ＡＤ_ｋＡ^Ｔ）、これらは例えば、キュムラント（ｃｕｍｕｌａｎｔｓ）行列、または異なる周派数帯域内の共分散行列、または異なる時間推移を有する相互相関行列、または異なる時間的間隔にわたって推定された共分散行列である；

−ギブンス回転行列は、対角線上のｃｏｓθ要素、対角線より下のｓｉｎθ要素および対角線より上の−ｓｉｎθ要素を含み、「ハイパボリック回転」行列は対角線上のｃｏｓｈφ要素、対角線より下のｓｉｎｈφおよび対角線より上のｓｉｎｈφ要素を含む；

−ギブンス回転行列は、

タイプのものであり、「ハイパボリック回転」行列は、

タイプのものである。

本発明は同様に、少なくとも二つのセンサーによって記録された混合された信号を複数の成分信号へと分離するための装置に関するものであり、各混合された信号は独立したソースに由来する成分信号の瞬間線形結合に対応し、ソースに由来する成分信号のベクトルを混合行列に乗じた積として各混合された信号を記すことができるものであり、該装置は、さらに結合対角化用手段を備え、該手段では一つまたは複数のギブンス回転行列および一つまたは複数のハイパボリック回転行列に対応する複数の回転行列の積から混合行列が推定されることを特徴とする。

本発明は同様に、コンピュータで使用することのできる媒体に記録されたプログラムコード命令を含むコンピュータプログラム製品にも関するものであり、該プログラム製品は上述の方法を応用するための命令を含むことを特徴とする。

本発明のその他の特徴、目的および利点は以下の記述からさらに明らかとなるが、この記述は、純粋に例証であって限定的なものではなく、添付図面を参照に読まれるべきものである。

本発明にしたがった方法の実施態様を示す。 混合された信号を分離するための三つの異なる方法（Ｊ−Ｄｉ、ＦＦＤｉａｇおよびＬＵＪ１Ｄ）について、サイズＮ＝５０のＫ＝５０の行列のセットを用いた、５０回の独立した試行での、多くのスイープ（ｓｗｅｅｐｓ）に応じて、対角性に向かう収束を示す図である。 混合された信号を分離するための三つの異なる方法（Ｊ−Ｄｉ、ＦＦＤｉａｇおよびＬＵＪ１Ｄ）について、サイズＮ＝５０のＫ＝５０の行列のセットを用いた５０回の独立した試行での、多くのスイープに応じて、分離基準（ｓｅｐａｒａｔｉｏｎ ｃｒｉｔｅｒｉｏｎ）の収束を示す図である。 Ｋ＝Ｎ（Ｎ＋１）／２＝２１０ 四次キュムラントスライスＮ＝２０およびσ＝０．１である場合の、５０回の独立した試行での、多くのスイープに応じて、分離基準の収束を示す図である。 本発明に係る混合された信号を分離するためのシステムの実施態様を表示したものである。

当業者にとって公知の以下の代数的性質が、以降で使用されている：
−二つの行列ＡおよびＢの積の転置：
（ＡＢ）^Ｔ＝Ｂ^ＴＡ^Ｔ
−二つの行列ＡおよびＢの積のインバース：
（ＡＢ）^−１＝Ｂ^−１Ａ^−１
−ギブンス回転：
Ｇ（θ）^−１＝Ｇ（θ）^Ｔ＝Ｇ（−θ）
−ハイパボリック回転
Ｈ（φ）^−１＝Ｈ（−φ）
−ハイパボリック回転の対称性：
Ｈ（φ）＝Ｈ（φ）^Ｔ

本発明に係る方法の一つの実施態様についてここで図１を参考にしながらより詳細に記述する。

混合体を独立した成分に分解することは、混合体の中で相互に隠れている現象の各々の解釈、分析または処理を容易にすることから有利である。

これらの現象は、例えば以下のものであってよい：
−脳波信号（ＥＥＧ）に寄与するさまざまな（脳および筋肉）電気現象、
−狭周波数帯域のさまざまな無線送信機、
−溶液中に存在する、異なる蛍光色素からの発光（ｌｉｇｈｔ ｅｍｉｓｓｉｏｎ）、
−同じ分子の異なる混合体の質量スペクトル。

異なる寄与を分離することにより、このとき次のことが可能である。
−筋肉寄生（ｍｕｓｃｕｌａｒ ｐａｒａｓｉｔｅ）を免れたＥＥＧ信号を研究すること（明滅（ｂｌｉｎｋ）、ＥＣＧなど）；
−各高周波ソースからの信号を別々に復調すること；
−各蛍光色素の濃度を（例えばサイトメトリによって）推定すること、または、
−質量分析法を通して研究された媒質内に存在する、異なる分子を識別すること。

ソースの信号は、時間（ＥＥＧ、高周波数）、波長（蛍光性）、質量（質量スペクトル）などであってもよい可変数量に左右される。

以下では、ソース由来の信号が時間に左右されるという点を考慮に入れながら、本発明に係る方法を提示するが、これによってこの方法の全体的内容が限定されることはない。

ソースの分離の問題は、以下の要領で数学的にモデル化される。
Ｍ個のソースの各々は、他のソースが発出する信号とは独立した形で信号ｓ_ｍ（ｔ）を発出する。
これらのＭ個の寄与の（所与の瞬間ｔにおける）混合体は、Ｎ個の全く異なる混合体を提供するＮ個のセンサーのネットワーク（ここでＮ＞ＭまたはＮ＝Ｍ）により観察される。
したがって、ｘ（ｔ）と記される瞬間ｔにおいて観察されたデータは、多次元である。

所与の瞬間におけるこれらのＮ回の観察のベクトルは、行列Ａ（Ｎ行、Ｍ列）つまりいわゆる「混合行列」とソース信号のベクトルｓ（ｔ）との積として記してよい。すなわち

先の等式を再公式化することで測定値ｘ（ｔ）に対する異なるソースの寄与ｓ_ｍ（ｔ）を示してもよい。

本発明に係る方法によると、混合行列Ａに対するいかなる先験的情報もなくこの混合行列を推定することができ（この行列は単に非特異行列であることだけが仮定されている）、このときソース信号が推定される。

ソース信号は独立していると仮定され、その非ガウス性（ｎｏｎ−Ｇａｕｓｓｉａｎ ｃｈａｒａｃｔｅｒ）、そのスペクトル差、またはその非定常性などのその統計的特性が利用される。

第一のステップ１００においては、Ｋ個の推定された行列Ｒ：_ｋの族が、センサーネットワークにより記録されたデータから（１とＫの間に含まれるｋについて）計算される。
これらの推定された行列Ｒ：_ｋは、以下の構造を有する行列Ｒ_ｋと（推定エラー内まで）等しい：

なお式中、Ｄ_ｋは対角行列であり、Ａ^Ｔは転置された混合行列Ａである。当然のことながら、当業者であれば、行列Ｄ_ｋがブロック状に対角であってよいということを理解すると思われる。

等式の第３項は、ｒａｎｋが１のＭ個の行列の加重和（ａ_ｍａ^Ｔ _ｍ）としての、Ｒ_ｋ個の行列を示す第２項の再公式化にすぎない。項Ｄ_ｋ（ｍ，ｍ）ａ_ｍａ^Ｔ _ｍは、行列Ｒ_ｋに対するソースｍの寄与である。

これらの状況に応じて、これらの推定行列Ｒ：_ｋは、３次または４次のキュムラントの行列（ＪＡＤＥ方法［１］、［１ｂ］）、異なる周波数帯域内の共分散行列［３］、異なる時間推移の相互相関行列（ＳＯＢＩ方法［２］）または異なる時間的間隔にわたり推定された共分散行列［３］であってよい。当然のことながら、これらの推定行列は、当業者にとって公知のその他のタイプのものであってよい。

第二のステップ２００においては、推定行列の結合対角化が、以下でより詳細に記述されているように実施される。

本発明に係る方法は、Ａ行列をＲ：_１，．．．，Ｒ：_ｋ，．．．Ｒ：_Ｋ行列から計算することを可能にする。

行列Ｂは、すべての行列ＢＲ_ｋ：Ｂ^Ｔが可能なかぎり「対角」となるように求められる：
ＢＲ：_ｋＢ^Ｔ＝Ｄ_ｋ
Ｂは、以下の形の基本Ｎ×Ｎ行列の積として構築される：

ここで、ｉ＜ｊは、１とＮの間の任意の二つのインデックスであり、Ｉ_Ｐは、ｐ×ｐサイズの恒等行列を表わし、下記〔数８〕は、下記〔数１１〕の等式を満たすような、下記〔数９〕タイプのギブンス行列（なお式中θは、ギブンス回転の角度である）と〔数１０〕タイプのハイパボリック回転行列（なお式中φはハイパボリック回転角度である）との積である。

したがって、先行技術の対角化方法との差異は、行列Ｂ（ひいては混合行列Ａ）が、回転行列と一つまたは複数のギブンス回転行列および一つまたは複数のハイパボリック回転行列との積から推定されるという点にある。

「ハイパボリック回転行列」とは、φのハイパボリック関数に対応する要素を含む行列のことである。

より詳細には、推定行列Ｒ：_ｋの族には以下のものが乗算される：
−各推定行列Ｒ：_ｋの右側に、ギブンス回転行列とハイパボリック回転行列との積。当業者にとっては、該方法の等価のインプリメンテーション（ｅｑｕｉｖａｌｅｎｔ ｉｍｐｌｅｍｅｎｔａｔｉｏｎ）が回転の反転（ｉｎｖｅｒｓｉｏｎ）からもたらされるであろうことは明白である（最初にギブンス回転で次にハイパボリック回転ではなく、ハイパボリック回転で次にギブンス回転）。
−そして、各推定行列Ｒ：_ｋの左側に、この積の転置行列：すなわち

となる。ここで

である。

この積の結果としてもたらされる行列、下記〔数１４〕の非対角項を最小化する角度θ_０およびφ_０が次に求められる。

この演算は、ギブンスおよびハイパボリック回転行列の要素を、結果として得られる行列の以下の要素へと移動させることによって、結果として得られた行列について反復される：

ここで

であって、インデックスｉおよびｊの全ての対について同様である。

これらの演算は、ｃｏｓθおよびｃｏｓｈφがそれぞれに１に向かうような、全てのギブンス回転角度θおよび全てのハイパボリック回転角度φとなるまで、インデックスｉおよびｊの全ての対について複数回反復される。

当然のことながら、当業者であれば、（該方法の大域収束性を意味する）他の停止基準が可能であるということを知っている。例えば、対角要素の全体的増加または非対角要素の全体的減少を監視して、反復停止を決定してもよい。

行列Ｍ（ｉ，ｊ，θ，φ）^Ｔ×Ｒ：_ｋ×Ｍ（ｉ，ｊ，θ，φ）の非対角項を全体的に最小化する角度θおよびφは、以下の要領で得ることができる：

−行列Ｃを以下のように構築する。

−Ｃ^Ｔ×Ｃおよび対角行列（−１，１，１）の一般固有分解（ｇｅｎｅｒａｌｉｚｅｄ ｅｉｇｅｎｄｅｃｏｍｐｏｓｉｔｉｏｎ）を計算する。固有ベクトルはｚ＞０および−ｘ^２＋ｙ^２＋ｚ^２＝１となるような正の最小の一般固有値のν＝（ｘ，ｙ．ｚ）^Ｔとして記される。

を計算する。

を計算する。

問題は以下の等式を解くことのみにあることから、その他の等価の式も存在することを当業者であれば認識する：

上述の式は、該方法の収束の終りでｘが０に向かう場合に優れた数値的安定性を有するという利点をもつ。

コンピュータ計算の有限精度に適応させるためには、その他のマイナーな修正が有用であるかもしれない。例えば、正の最小の固有値は、理論的には実際上ごくわずかに負であるかもしれない。それでも、それを保持することが必要である。これらの点は、当業者にとって明白である。

行列Ｃからの角度θおよびφの計算のより完全な説明を以下で記述する。

論理的進行は以下の通りである。

ｋ＝１，．．．，Ｋについて、結果として得られる行列の積（Ｒ’_ｋと記す）はＲ’_ｋ＝Ｍ（ｉ，ｊ，θ，φ）^Ｔ×Ｒ：_ｋ×Ｍ（ｉ，ｊ，θ，φ）と計算される：

θのコサインとθのサインの積ならびにφのハイパボリックコサイン（ｈｙｐｅｒｂｏｌｉｃ ｃｏｓｉｎｅｓ）とφのハイパボリックサイン（ｈｙｐｅｒｂｏｌｉｃ ｓｉｎｅｓ）の積が出現し、倍角（２θ、２φ）のコサインおよびサインへと変換され、これらは項のグルーピングの時点で、以下の等式を提供する：

ここで、

このとき、結果として得られる行列の積Ｍ（ｉ，ｊ，θ，φ）^Ｔ×Ｒ：_ｋ×Ｍ（ｉ，ｊ，θ，φ）の非対角要素（ｉ，ｊ）の二乗の和は二次形式ν（θ，φ）^ＴＣ^ＴＣν（θ，φ）に等しい。

ベクトルν（θ，φ）がν（θ，φ）^ＴＪν（θ，φ）＝１（なおＪ＝ｄｉａｇ（−１，１，１））により定義される双曲面（ｈｙｐｅｒｂｏｌｏｉｄ）を記述する、という点に留意されたい。

制約ν（θ，φ）^ＴＪν（θ，φ〉＝１の下で二次形式ν（θ，φ）^ＴＣ^ＴＣν（θ，φ）を最大化するベクトルν（θ，φ）が求められる。このベクトルが必然的に行列対（Ｃ^ＴＣ，Ｊ）の三つの一般固有ベクトルのうちの一つであるということは周知である。

この一般固有分解の研究は、これが（Ｃ^ＴＣ，Ｊ）の正の最小の一般固有値に対応する一般固有値ν＝（ｘｙｚ）^Ｔであることを示している：

最終的に、積Ｍ（ｉ，ｊ，θ，φ）^Ｔ×Ｒ：_ｋ×Ｍ（ｉ，ｊ，θ，φ）を計算するために、θおよびφの最適値は不要である。その（従来のおよびハイパボリックの）コサインおよびサインの値のみが必要とされ、これは例えば以下の式から得られる：

したがって該方法は、１とＮの間のｉ＜ｊについて、収束まで数回、
Ｍ（ｉ，ｊ，θ，φ）^Ｔ×Ｒ_ｋ：×Ｍ（ｉ，ｊ，θ，φ）、
の積の反復を含む。

混合行列Ａはこのとき

のインバース（ｉｎｖｅｒｓｅ）に等しい。混合行列Ａがひとたび決定されると、ソース信号を分離することが可能になる。

Ｍａｔｌａｂ言語での例示的なインプランテーション（ｉｍｐｌａｎｔａｔｉｏｎ）を以下に示す：
ｆｕｎｃｔｉｏｎ［Ａ，ｎｕｍｂｅｒ＿ｏｆ＿ｌｏｏｐｓ，ｍａｔｒｉｃｅｓ］＝Ｊ＿Ｄｉ［ｍａｔｒｉｃｅｓ，ｔｈｒｅｓｈｏｌｄ］

％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％
％
％ Ｊ＿Ｄｉ．ｍ ｗｒｉｔｔｅｎ ｔｈｅ １９ｔｈ ｏｆ ｄｅｃｅｍｂｅｒ ２００７ ｂｙ Ａｎｔｏｉｎｅ Ｓｏｕｌｏｕｍｉａｃ．
％
％ Ｔｈｉｓ ｓｏｆｔｗａｒｅ ｃｏｍｐｕｔｅｓ ｔｈｅ Ｊｏｉｎｔ Ｄｉａｇｏｎａｌｉｓａｔｉｏｎ ｏｆ ａ ｓｅｔ ｏｆ ｍａｔｒｉｃｅｓ．
％
％ ｉｎｐｕｔ：

％ ｍａｔｒｉｃｅｓ： ａｒｒａｙ ｏｆ ｍａｔｒｉｃｅｓ ｔｏ ｂｅ ｊｏｉｎｔｌｙ ｄｉａｇｏｎａｌｉｚｅｄ
％ ｔｈｒｅｓｈｏｌｄ： ｔｈｒｅｓｈｏｌｄ ｆｏｒ ｓｔｏｐ ｃｒｉｔｅｒｉｏｎ ｏｎ ｍａｘｉｍａｌ ｓｉｎｅｓ
％
％ ｏｕｔｐｕｔｓ：
％ Ａ： ｅｓｔｉｍａｔｅｄ ｍｉｘｔｕｒｅ ｍａｔｒｉｘ
％ ｎｕｍｂｅｒ＿ｏｆ＿ｌｏｏｐｓ：ｎｕｍｂｅｒ ｏｆ ｌｏｏｐｓ ｎｅｃｅｓｓａｒｙ ｆｏｒ ｃｏｎｖｅｒｇｅｎｃｅ
％ ｍａｔｒｉｃｅｓ：ａｒｒａｙ ｏｆ ｍａｔｒｉｃｅｓ ａｆｔｅｒ ｊｏｉｎｔ ｄｉａｇｏｎａｌｉｚａｔｉｏｎ
％
％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％
％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％

％ ｉｎｉｔｉａｌｉｚａｔｉｏｎ
ｍａｔ ｓｉｚｅ ＝ ｓｉｚｅ（ｍａｔｒｉｃｅｓ，１）；
ｎｂ ｍａｔ ＝ ｓｉｚｅ（ｍａｔｒｉｃｅｓ，３）；
Ａ ＝ ｅｙｅ（ｍａｔ＿ｓｉｚｅ）；
ｇｏ＿ｏｎ ＝ １；
ｎｕｍｂｅｒ＿ｏｆ＿ｌｏｏｐｓ ＝ −１；

％ ｊｏｉｎｔ ｄｙａｄｉｚａｔｉｏｎ ｌｏｏｐｓ
ｗｈｉｌｅ ｇｏ＿ｏｎ

ｆｏｒ ｎ＝ ２：ｍａｔ＿ｓｉｚｅ
ｆｏｒ ｍ＝ １：ｎ−１

％ ｉｎｉｔｉａｌｉｚａｔｉｏｎ
ｓ ｍａｘ ＝ ０．０；
ｓｈ ｍａｘ ＝ ０．０；

％ Ｇｉｖｅｎｓ ａｎｄ ｈｙｐｅｒｂｏｌｉｃ ａｎｇｌｅｓ ｃａｌｃｕｌａｔｉｏｎ
Ｃ＝［（ｓｑｕｅｅｚｅ（ｍａｔｒｉｃｅｓ（ｍ，ｍ，：））＋ｓｑｕｅｅｚｅ（ｍａｔｒｉｃｅｓ（ｎ，ｎ，：）））／２．．．
（ｓｑｕｅｅｚｅ（ｍａｔｒｉｃｅｓ（ｍ，ｍ，：））−ｓｑｕｅｅｚｅ（ｍａｔｒｉｃｅｓ（ｎ，ｎ，：）））／２．．．
ｓｑｕｅｅｚｅ（ｍａｔｒｉｃｅｓ（ｍ，ｎ，：））］；
［Ｖ，Ｄ］＝ｅｉｇ（Ｃ’^＊Ｃ，ｄｉａｇ（［−１ １ １］））；
％ ｎｏｒｍａｌｉｓａｔｉｏｎ（−１ １ １ ） ｏｆ ｇｅｎｅｒａｌｉｚｅｄ ｅｉｇｅｎｖｅｃｔｏｒｓ
Ｖ＝Ｖ^＊ｄｉａｇ（１．／ｓｑｒｔ（ａｂｓ（ｄｉａｇ（Ｖ’^＊ｄｉａｇ（［−１ １１］）^＊Ｖ））））；
％ ｃｏｍｐｕｔａｔｉｏｎ ｏｆ ｔｈｅ ｍｅｄｉｕｍ ｇｅｎｅｒａｌｉｚｅｄ ｅｉｇｅｎｖａｌｕｅ
［ｓｏｒｔＤ，ｉｎｄｅｘ］＝ ｓｏｒｔ（ｄｉａｇ（Ｄ））；
ｉｎｄ＿ｍｉｎ ＝ ｉｎｄｅｘ（２）；
ｖ＝Ｖ（：，ｉｎｄ＿ｍｉｎ）^＊ｓｉｇｎ（Ｖ（３，ｉｎｄ＿ｍｉｎ））；

％ ｃｏｍｐｕｔａｔｉｏｎ ｏｆ ｔｒｉｇｏｎｏｍｅｔｒｉｃ ａｎｄ ｈｙｐｅｒｂｏｌｉｃ ｃｏｓｉｎｅｓ ａｎｄ ｓｉｎｅｓ
ｃｈ ＝ ｓｑｒｔ（（１＋ｓｑｒｔ（１＋ｖ（１）^∧２））／２）；
ｓｈ ＝ ｖ（１）／（２^＊ｃｈ）；
ｃ ＝ ｓｑｒｔ（（１＋ｖ（３）／ｓｑｒｔ（１＋ｖ（１）^∧２））／２）；
ｓ ＝ −ｖ（２）／（２^＊ｃ^＊ｓｑｒｔ（１＋ｖ（１）^∧２））；
％ ｍａｘｉｍａｌ ｓｉｎｅ ａｎｄ ｓｉｎｈ ｃｏｍｐｕｔａｔｉｏｎ ｆｏｒ ｓｔｏｐ ｔｅｓｔ
ｓ ｍａｘ ＝ ｍａｘ（ｓ＿ｍａｘ，ａｂｓ（ｓ））；
ｓｈ＿ｍａｘ ＝ ｍａｘ（ｓｈ＿ｍａｘ，ａｂｓ（ｓｈ））；

％ ｐｒｏｄｕｃｔ ｏｆ Ｇｉｖｅｎｓ ａｎｄ ｈｙｐｅｒｂｏｌｉｃ ｒｏｔａｔｉｏｎｓ
ｒｍ１１ ＝ ｃ^＊ｃｈ − ｓ^＊ｓｈ；
ｒｍ１２ ＝ ｃ^＊ｓｈ − ｓ^＊ｃｈ；
ｒｍ２１ ＝ ｃ^＊ｓｈ ＋ ｓ^＊ｃｈ；
ｒｍ２２ ＝ ｃ^＊ｃｈ ＋ ｓ^＊ｓｈ；

％ ｍａｔｒｉｃｅｓ Ｇｉｖｅｎｓ ｕｐｄａｔｅ
ｈ＿ｓｌｉｃｅ１ ＝ ｓｑｕｅｅｚｅ（ｍａｔｒｉｃｅｓ（ｍ，：，：））；
ｈ＿ｓｌｉｃｅ２ ＝ ｓｑｕｅｅｚｅ（ｍａｔｒｉｃｅｓ（ｎ，：，：））；
ｍａｔｒｉｃｅｓ（ｍ，：，：）＝ ｒｍ１１^＊ｈ＿ｓｌｉｃｅ１ ＋ ｒｍ２１^＊ｈ＿ｓｌｉｃｅ２；
ｍａｔｒｉｃｅｓ（ｎ，：，：）＝ ｒｍ１２^＊ｈ＿ｓｌｉｃｅ１ ＋ ｒｍ２２^＊ｈ＿ｓｌｉｃｅ２；
ｖ＿ｓｌｉｃｅ１ ＝ ｓｑｕｅｅｚｅ（ｍａｔｒｉｃｅｓ（：，ｍ，：））；
ｖ＿ｓｌｉｃｅ２ ＝ ｓｑｕｅｅｚｅ（ｍａｔｒｉｃｅｓ（：，ｎ，：））；
ｍａｔｒｉｃｅｓ（：，ｍ，：）＝ ｒｍ１１^＊ｖ＿ｓｌｉｃｅ１ ＋ ｒｍ２１^＊ｖ＿ｓｌｉｃｅ２；
ｍａｔｒｉｃｅｓ（：，ｎ，：）＝ ｒｍ１２^＊ｖ＿ｓｌｉｃｅ１ ＋ ｒｍ２２^＊ｖ＿ｓｌｉｃｅ２；

％ ｄｙａｄｓ Ｇｉｖｅｎｓ ｕｐｄａｔｅ
ｃｏｌ１ ＝ Ａ（：，ｍ）；
ｃｏｌ２ ＝ Ａ（：，ｎ）；
Ａ（：，ｍ）＝ ｒｍ２２^＊ｃｏｌ１ − ｒｍ１２^＊ｃｏｌ２；
Ａ（：，ｎ）＝ −ｒｍ２１^＊ｃｏｌ１ ＋ ｒｍ１１^＊ｃｏｌ２；

ｅｎｄ
ｅｎｄ
ｎｕｍｂｅｒ＿ｏｆ＿ｌｏｏｐｓ＝ｎｕｍｂｅｒ＿ｏｆ＿ｌｏｏｐｓ＋１；

％ ｒｅｎｏｒｍａｌｉｚａｔｉｏｎ ｏｆ ｔｈｅ Ａ ｃｏｌｕｍｎｓ（ｅｑｕａｌ ｎｏｒｍ）
ｃｏｌ＿ｎｏｒｍ ＝ ｓｑｒｔ（ｓｕｍ（Ａ．^∧２，１））；
ｄ ＝（ｐｒｏｄ（ｃｏｌ＿ｎｏｒｍ）．^∧（１／ｍａｔ＿ｓｉｚｅ））．／ｃｏｌ＿ｎｏｒｍ；
Ａ ＝ Ａ．^＊ｒｅｐｍａｔ（ｄ，ｍａｔ＿ｓｉｚｅ，１）；
ｍａｔｒｉｃｅｓ ＝ ｍａｔｒｉｃｅｓ．^＊ｒｅｐｍａｔ（（１．／ｄ）’^＊（１．／ｄ），［１ １ ｎｂ＿ｍａｔ］）；

％ ｓｔｏｐ ｔｅｓｔ
ｉｆ ｓ＿ｍａｘ＜ｔｈｒｅｓｈｏｌｄ ＆ ｓｈ＿ｍａｘ＜ｔｈｒｅｓｈｏｌｄ
ｇｏ＿ｏｎ＝０；
ｄｉｓｐ（［’ｓｔｏｐ ｗｉｔｈ ｍａｘ ｓｉｎｕｓ ａｎｄ ｓｉｎｈ＝’ ｎｕｍ２ｓｔｒ（［ｓ＿ｍａｘ ｓｈ＿ｍａｘ］）］）
ｅｎｄ
ｉｆ ｎｕｍｂｅｒ＿ｏｆ＿ｌｏｏｐｓ＞１００
ｇｏ＿ｏｎ ＝ ０；
ｗａｒｎｉｎｇ（’ｍａｘｉｍｕｍ ｎｕｍｂｅｒ ｏｆ ｌｏｏｐｓ ｈａｓ ｂｅｅｎ ｒｅａｃｈｅｄ’）
ｅｎｄ
ｅｎｄ

計算の複雑性は以下のような近似値を有する：
スイープ数×２ＫＮ^３
なお式中、
−スイープ数は、インデックス対（ｉ，ｊ）全体について反復が行われる回数であり、
−ＫはＲ：_ｋ行列の数であり、
−ＮはＲ：_ｋ行列のサイズである。

本発明および既存のアルゴリズム例えばＦＦＤｉａｇ［４］、ＤＮＪＤ［６］、ＬＵＪ１Ｄ［７］などは、２ＫＮ^３に等しくスイープすることにより等価の計算コストを有する。本発明の利点は、スイープ数にある。より具体的には、当該方法によると、先行技術の公知の方法と比べてスイープ数を削減してよい。

一例として、本発明に係る方法（以下「Ｊ−Ｄｉ」と呼ぶ）とＦＦＤｉａｇ（Ｚｉｅｈｅ［４］）およびＬＵＪ１Ｄ（Ａｆｓａｒｉ［７］）方法の比較は以下の通りに提示される。すなわちこの比較は、先行技術の方法が、本発明の方法のものの少なくとも二倍のスイープ数を有するという事実を示す。

その上、先行技術の一部の方法は、特定の行列セットに特化されている（Ｐｈａｍ［３］：正値定符号行列、Ｖｏｌｌｇｒａｆ［５］：少なくとも一つの正値定符号行列を含むセット）。本発明に係る方法はあらゆるタイプの行列セットの処理を可能にする。

図５は本発明に係る混合された信号を分離するためのシステムの実施態様を示す。このシステムは、例えばプロセッサーのような処理手段１０（一つまたは複数のギブンス回転行列および一つまたは複数のハイパボリック回転行列に対応する回転行列の積から混合行列Ａを推定する結合対角化用）、例えば表示スクリーンのような表示手段２０、そして、例えばキーボードやマウスのような入力手段３０を備える。システムは、処理すべき混合された信号（ｘ_１（ｔ），．．．，ｘ_Ｎ（ｔ））を受信するため通信手段に接続されている。

読者は、本明細書の教示から実質的に逸脱することなく上述に記載の本発明に対し数多くの修正を加えることができると認識するものである。

例えば、ギブンス回転（Ｇ（θ））とハイパボリック回転（Ｈ（φ））行列の積は、ハイパボリック回転行列にギブンス回転行列を乗じたもの、または、ギブンス回転行列にハイパボリック回転行列を乗じたもののいずれかに対応することができる。

結果として、このタイプのすべての修正は、添付の請求の範囲に包含されるものと意図される。

「ＡＮＮＥＸ」
添付の図面により、上述の方法をよりよく理解することができる。先に用いた一部の表記法は、以下では異なる可能性がある。いずれにせよ、表記法の変更は読者の理解を容易にする目的で特定されるものである。

Ｉ．要約
独立成分分析およびブラインドソース分離（ｂｌｉｎｄ ｓｏｕｒｃｅ ｓｅｐａｒａｔｉｏｎ）のアプリケーション用として、行列セットの非直交結合対角化を計算するための新しいアルゴリズムが提案されている。

このアルゴリズムは、直交結合対角化用の、固有行列の同時近似対角化（Ｊｏｉｎｔ Ａｐｐｒｏｘｉｍａｔｅ ｄｉａｇｏｎａｌｉｚａｔｉｏｎ ｏｆ ｅｉｇｅｎｍａｔｒｉｃｅｓ）（ＪＡＤＥ）方法において最初に提案されたヤコビのようなアルゴリズムの拡張である。改善は、主として、直交混合行列の代りに等ノルム（ｅｑｕａｌ ｎｏｒｍ）の列および行列式１の混合行列を計算することからなる。この標的行列は、ギブンス回転のみならずハイパボリック回転および対角行列の連続的乗算によって反復的に構築される。合成データについて評価されたアルゴリズム性能は、収束速度および複雑性の視点から既存の方法と比較して有利である。

ＩＩ．序
結合対角化アルゴリズムは、ブラインドソース分離（ＢＳＳ）および独立成分分析（ＩＣＡ）のアプリケーションにおいて使用される。

これらのアプリケーションでは、瞬時に混合された独立した信号（ソース）を分離すること、すなわち、混合体のみを観察する行列Ａによる一時的コンボリューション（ｔｅｍｐｏｒａｌ ｃｏｎｖｏｌｕｔｉｏｎ）無しで、かつ、ソースの統計的独立性以外の先行情報無しで、独立した信号（ソース）を分離することが目標とされる。

［１］、［１ｂ］、［２］、［３］、［４］、［５］、［１１］などの数多くのＢＳＳ方法は、Ｎ×ＮサイズのＫ個の対称行列Ｒ_１，．．．，Ｒ_ｋ，．．．，Ｒ_ＫのセットＲの結合対角化ステップを含む。手元のアプリケーションおよびソース信号の対応する予想特性によると、考慮対象の行列Ｒ_ｋは、例えば異なる時間窓で、三次または四次（［１］，［１ｂ］）以上のキュムラントスライスに基づいて推定される共分散行列、または時間推移を伴う相互相関行列［２］などである。これらの行列はすべて、共通の構造を共有する：

なお式中、Ａは混合行列であり、Ｄ_ｋは対角行列である。これらの条件下で、

［数２９］により定義される分離行列は、各Ｒ_ｋを対角化する特性、すなわち１≦ｋ≦ＫについてＢＲ_ｋＢ^Ｔ＝Ｄ_ｋにより特徴づけられる。弱い仮説の下で、置換および不定スケーリング（ｓｃａｌｉｎｇ ｉｎｄｅｔｅｒｍｉｎａｔｉｏｎ）までは、Ｒ_ｋ行列の結合対角化手順により分離行列また同等に混合行列を推定するのにこの特性で充分である。実際には、Ｒ_ｋ行列は、一部の推定エラーによって破損され、Ａ（またはＢ）はセットＲの近似結合対角化を介して効率良く推定される。

その上、例えばＪＡＤＥ［１］、［１ｂ］またはＳＯＢＩ［２］のような最初の結合対角化技術には、Ｒ内の一つの行列の、一般的には（ノイズ無しの）混合されたソースの共分散行列の、正定性が必要とされる。これに関連して、このとき、混合されたソースを予め白色化することが可能であり、こうして直交（または複雑なケースではユニタリ）混合行列を計算する上での結合対角化の問題が削減される。残念ながら、この白色化ステップは、混合されたソースの共分散行列内の推定エラーに起因して決して正確ではない。この予備的エラーは、後続する直交結合対角化（ＯＪＤ）により補正できない。結果として、数多くの製作者が、予備白色化を回避するための非直交結合対角化（ＮＯＪＤ）方法を開発した。これらの最適化方法は、正値定符号行列のセット［３］、またはＱＤｉａｇのような一つの正値定符号行列を含むセット［１２］、［５］または最終的に任意の行列セット［１１］、［４］、［７］に対応している。

ＪＡＤＦなどのヤコビ技術の優れた収束特性を享受する代数アプローチに基づいたＢＳＳおよびＩＣＡアプリケーション用の任意の行列セットのための新たな非直交結合対角化アルゴリズムが提案されている。この結果を得るために、混合行列Ａは同じく基本行列の乗算によって構築される。しかしながら、Ａは直交でないことから、下位群として直交行列群を含む新しい乗法群、およびギブンス回転を一般化する基本行列の新しいセットを定義することが必要である。これは、ギブンス回転、ハイパボリック回転および対角行列の連続的乗算により等ノルムの列を有し１に等しい行列式の混合行列を反復的に構築することによって達成できる。

以下の第三節では、混合行列の列ノルムおよび行列式の制約およびそれがＮＯＪＤ問題の不確定にもたらす帰結について吟味する。

第四節では、行列式１の行列をいかにしてギブンス回転、ハイパボリック回転および対角行列の積として生成することができるかを示す。

Ｊｏｉｎｔ−Ｄｉａｇｏｎａｌｉｚａｔｉｏｎ（結合対角化）を略してＪ−Ｄｉと呼ばれる新しい非直交近似結合対角化アルゴリズムは第五節で記述されており、その性能は、第六節内で合成データについて評価される。特に、例えば２００×２００サイズの２００個の行列などの大きい行列の大きいセットについてのシミュレーションが提供されている。これは、電極数が１２８個または２５６個に達する場合の脳波図などの利用分野にとって有利である。

最後に、［６］においてＷａｎｇ、ＬｉｕおよびＺｈａｎｇによって提案された最近のアルゴリズムとの短い複雑性比較が付加されている。

ＩＩＩ．１に等しい行列式および等ノルム列の行列
ＩＩＩ．Ａ．特性

モデル（１）内の任意の正則正方混合行列Ａが１に等しい行列式および等ノルムの列を有するものと仮定される。実際のところ、（１）の中に正則対角行列Ｄ_Ａを挿入することは常に可能である。

行列Ｄ_Ａ ^−１Ｄ_ｋＤ_Ａ ^−１は明らかに対角であり、Ｄ_Ａは、ＡＤ_Ａが所要の行列式および列ノルム特性を確認するような形で調整可能である。

を考慮するだけで充分である。

ここで、ａ_ｎはＡのｎ番目の列であり、｜｜ａ_ｎ｜｜はそのユークリッド・ノルムである。混合行列に対するこれらの制約は、数値的精度にとって有用な一種の正規化を提供する。特に、ＡＤ_Ａの条件数はＡの条件数よりも小さいことが多いが、これは通則ではない。

ＩＩＩ．Ｂ．混合行列の一意性条件および不確定
［９］、［１０］に記述されているように、Ｄ_ｋの対角線（ＤはＫ×Ｎであり、その（ｋ，ｎ）エントリはＤ_ｋのｎ番目の対角エントリである）の垂直連結（ｖｅｒｔｉｃａｌ ｃｏｎｃａｔｅｎａｔｉｏｎ）の結果もたらされる行列Ｄ内に共線列（ｃｏｌｉｎｅａｒ ｃｏｌｕｍｎｓ）が全く存在しない場合、そしてその場合にのみＮＯＪＤ問題に対する一意解が存在する。Ρ_ｎｍが、Ｄのｎ番目およびｍ番目の列の角度のコサインの絶対値であり、

である場合には、一意性条件を単に

として再公式化できる。一意性のモジュラスと呼ばれるこのパラメータρは、ρが１に近いほどＮＯＪＤはより難しいということを明白な規則として、ＮＯＪＤ問題の難度を示すものでもある。

ＮＯＪＤ問題の「一意解」が、混合行列の列の古典的置換およびスケール不確定までのみ「一意的である」ということが理解される。（本発明の）単位行列式および等しい列ノルム特性は、これらの不確定性をわずかに変える。実際のところ、置換不確定は存在し続けるがスケール不確定が符号不確定（ｓｉｇｎ ｉｎｄｅｔｅｒｍｉｎａｔｉｏｎ）まで削減されることは容易に示すことができる。より厳密に言うと、偶数の混合行列の列に−１を乗じることができる。

ＩＶ．特別な線形群全体にわたる基本行列の新たな部類
この節の最終目標は、反復した乗算によりセットＡ（Ｒ，Ｎ）および群ＳＬ（Ｒ，Ｎ）を生成する基本行列セットを提案することにある。

ＩＶ．Ａ．Ａ（Ｒ，２）内の一意的分解
Ａ（Ｒ，２）（ｄｅｔ（Ａ）＝１および｜｜ａ_１｜｜＝｜｜ａ_２｜｜）内の行列Ａおよびその極分解

を考慮し、式中Ｇは直交であり、Ｈは対称的正である。ＡおよびＨの行列式は両方共正であることから、Ｇの行列式はプラス１に等しく、Ｇは必然的に以下の〔数３７〕でギブンス回転Ｇ（θ）である。

ここで、同じく行列式がプラス１に等しいＨが実際にはハイパボリック回転であることを示す。

Ｈのエントリをｅ、ｆおよびｇと表示する。

Ｈ＝Ｇ（−θ）Ａは、Ａと同様、等ノルムの列を有する。したがって、ｅ^２＋ｆ^２＝ｆ^２＋ｇ^２であり、Ｈは構成上正であり、したがってｅ＞０、ｇ＞０そしてｅ＝ｇである。Ｈの行列式は１に等しく、このときｅ^２−ｆ^２＝１であり、ｅ＝ｃｏｓｈ（φ）およびｆ＝ｓｉｎｈ（φ）であるような一意的角度φが存在する。以下の〔数３９〕で最終的にＨ＝Ｈ（φ）である。

結論として、Ａ（Ｒ，２）の任意の行列をギブンスおよびハイパボリック回転の一意的積に分解することができる。

ＩＶ．Ｂ． ＳＬ（Ｒ，２）中の複数解
ＳＬ（Ｒ，２）の任意の行列Ｆは、Ａ（Ｒ、２）内の行列Ａに２×２直交行列Ｄ（λ）を乗じた積へと、一意的に分解することができる。ここで、

である。ｆ_１、ｆ_２はＦの列であり、｜｜ｆ_１｜｜、｜｜ｆ_２｜｜はそれらのノルムである。上の節で示した通り、ギブンスおよびハイパボリック回転の積へと行列Ａを分解することができる。

最終的に、二つの一意的角度、すなわち−пと＋пの間のθおよびＲ内のφそして、一意の正の実数λが存在し、次のような関係となる：

さらなる分解を

の特異値分解（ｓｉｎｇｕｌａｒ ｖａｌｕｅ ｄｅｃｏｍｐｏｓｉｔｉｏｎすなわちＳＶＤ）から導出することができる。

必然的に、ｄｅｔ（Ｓ）＝１およびｄｅｔ（Ｕ）＝ｄｅｔ（Ｖ）＝±１である。ただし、ｄｅｔ（Ｕ）＝ｄｅｔ（Ｖ）＝＋１となるように、ＵとＶの第一列（または行）の符号を交換することがつねに可能であるという点に留意されたい。このことは、ＵおよびＶの両方をギブンス回転すなわちＵ＝Ｇ（θ_１）およびＶ＝Ｇ（θ_２）として選択できることを証明している。明らかに、Ｓ＝Ｄ（ｓ_１）であり、ここでｓ_１＞０は第一のＳ特異値である。以下の式が最終的に得られる。

さらに、任意のハイパボリック回転Ｈ（φ）の固有値分解（ＥＶＤ）から、以下が得られる：

そして結果として

これらの分解は、ＳＬ（Ｒ，２）の任意の行列を、ギブンス回転と対角行列、またはギブンス回転とハイパボリック回転の乗算によって生成できるということを証明している。さらに、対角行列を使用して、ＳＬ（Ｒ，２）中の任意の行列の列を正規化し、かくしてＡ（Ｒ，２）の行列を得ることができる。これらの結果はＮ次元まで拡張できる。

ＩＶ．Ｃ． ＳＬ（Ｒ，Ｎ）への拡張
以下の表記法が最初に導入される。

任意のｉ＜ｊ∈｛１，．．．，Ｎ｝Ｇ（θ，ｉ，ｊ）は古典的ギブンス回転

を表わし、Ｈ（φ，ｉ，ｊ）は、対応するハイパボリック回転

であり、Ｄ（λ，ｉ，ｊ）は対角行列

を表わす。

Ｎ＝２のケースときわめて類似して、ＳＬ（Ｒ，Ｎ）中の、行列Ｆの以下の〔数４９〕のＳＶＤが使用される。

直交行列ＵおよびＶの行列式は、ＵおよびＶの第一列の符号を交換することによりプラス１に等しく選択できる。行列式１（ＳＯ（Ｒ，Ｎ）中）の全ての直交行列をギブンス回転の積の中で分解できるということは周知である。構成上、Ｓは対角正定値であり、その行列式はプラス１に等しい（特異値ｓ_１，ｓ_２，．．．，ｓ_Ｎの積は１に等しい）。同様に以下を示すこともできる。

さらに、二次元の場合のように、Ｄ（λ，ｉ，ｊ）行列をギブンス回転およびハイパボリック回転の積に変換することができる。最終的に、同じ結論がＳＬ（Ｒ，２）内と同様ＳＬ（Ｒ，Ｎ）内にもあてはまる。すなわち、ＳＬ（Ｒ，Ｎ）内の全ての行列を、ギブンス回転と対角行列またはギブンス回転とハイパボリック回転を乗算することによって生成することができる。対角行列を用いてＳＬ（Ｒ，Ｎ）内の任意の行列の列を正規化し、Ａ（Ｒ，Ｎ）内の行列を得ることができる。

次節では、上述の三つの基本行列のパラメータをＮＯＪＤコンテクスト内で容易に最適化できるということを見ていく。

Ｖ． Ｊ−ＤＩアルゴリズム
本節では、Ｒ_１，…Ｒ_ｋ，…Ｒ_Ｋの正確な行列の推定に対応するセットを結合対角化することにより、基本行列Ｇ（θ，ｉ，ｊ）、Ｈ（φ，ｉ，ｊ）およびＤ（λ，ｉ，ｊ）の積として、１に等しい行列式および等ノルムの列の未知の行列Ａをどのように構築するかについて記述されている。

より厳密には、行列Ｇ（θ，ｉ，ｊ）およびＨ（φ，ｉ，ｊ）は、非対角エントリ（ｏｆｆ ｄｉａｇｏｎａｌ ｅｎｔｒｉｅｓ）を最小化するために使用され、一方でＤ（λ，ｉ，ｊ）行列は列を正規化するために使用される。

Ｖ．Ａ． （ｉ，ｊ）非対角エントリを最小化するためのパラメータθおよびφの最適化
ｋ＝１，…，Ｋについて、以下の要領でＲ：_ｋを更新することによって得られた行列をＲ’_ｋと表わす。

Ｒ’_ｋ＝Ｈ^Ｔ（φ，ｉ，ｊ）×Ｇ^Ｔ（θ，ｉ，ｊ）×Ｒ：_ｋ×Ｇ（θ，ｉ，ｊ）×Ｈ（φ，ｉ，ｊ） （１７）

Ｋ個の行列、Ｒ’_ｋのＲ’_ｋ［ｉ，ｊ］で表わされる（ｉ，ｊ）−エントリの二乗の和を最小化することが提案される：すなわち

この数量を最小化することは、例えばＪＡＤＥにおけるようにグローバルコスト関数（ｇｌｏｂａｌ ｃｏｓｔ ｆｕｎｃｔｉｏｎ）または非対角エントリの全てを最小化することとは等価ではない。それでも、以下の性能評価は、この反復的ローカル最小化が実際に行列セットＲを結合対角化することを示すものである。一定の操作の後、次のことがわかる：

なお式中

Ｒ’_ｋ中のＧおよびＨの係数の二乗が倍角（三角法角度およびハイパボリック角度）へと変換されるという点に留意されたい。Ｒ’_ｋ［ｉ，ｊ］の二乗和は、ν（θ，φ）^ＴＣ^ＴＣν（θ，φ）に等しい。したがって（１８）における最小化は次のようになる：

なお式中、ベクトルν（θ，φ）は、ユークリッド・ノルムではなくハイパボリックノルムで正規化される。実際のところ、−ｖ_１（θ，φ）^２＋ｖ_２（θ，φ）^２＋ｖ_３（θ，φ）^２＝１、すなわち、Ｊ＝ｄｉａｇ（−１，１，１）としてν（θ，φ）^ＴＪν（θ，φ）＝１となる。

ν（θ，φ）^ＴＪν（θ，φ）＝１で定義されるハイパボリック全体が−п／２≦θ≦п／２そしてＲ内のφという状態で、ν（θ，φ）により網羅されていることを証明することができる。喚言すると、二次制約（ｑｕａｄｒａｔｉｃ ｃｏｎｓｔｒａｉｎｔ）ｖ^ＴＪｖ＝１の下で二次数量（ｑｕａｄｒａｔｉｃ ｑｕａｎｔｉｔｙ）ｖ^ＴＣ^ＴＣｖの最小化が達成される。解ｖは（Ｃ^ＴＣ，Ｊ）の一般固有ベクトルであることがわかっている。

すなわちｖは、下記〔数５５〕となるようなものである。

この一般固有分解について以下で研究する。
その列が一般固有ベクトルである３×３行列をＶと表示し、（Ｃ^ＴＣ，Ｊ）の一般固有値の対角行列をΔと表わすと、Ｃ^ＴＣＶ＝ＪＶΔであり、したがってＶ^ＴＣ^ＴＣＶ＝Ｖ^ＴＪＶΔである。Ｖ^ＴＪＶは対角であり、Ｊのような一つの負の対角エントリと二つの正の対角エントリ（シルベスタの慣性の法則（Ｓｙｌｖｅｓｔｅｒ ｌａｗ ｏｆ ｉｎｅｒｔｉａ）を有することがわかっている。Ｃ^ＴＣは正であることから、ΔはＪと同じ慣性（ｉｎｅｒｔｉａ）を有し、次のものが存在する：
・負の一般固有値と負のＪ−ｎｏｒｍ ｖ^ＴＪｖ＝−１を有する一つの一般固有ベクトル、および
・正の一般固有値と正のＪ−ｎｏｒｍ ｖ^ＴＪｖ＝＋１を有する二つの一般固有ベクトル。

我々が求めている一般固有ベクトルは、ｖ^ＴＪｖが＋１に等しく、ｖ^ＴＣ^ＴＣｖが最小であるようなものであり、このことはすなわちそれが（Ｃ^ＴＣ，Ｊ）の最小の正一般固有値に対応することを意味している。

要約すると、最適なギブンスおよびハイパボリック角θおよびφは、以下の等式により定義される：

なお式中、ｖ＝（ｘ，ｙ，ｚ）^Ｔは、最小の正固有値の（Ｃ^ＴＣ，Ｊ）の一般固有ベクトルである。

固有ベクトルがその符号次第で決定されることから、ｚ＞０すなわちｃｏｓ２θが正となるように、ｖを選択することが常に可能である。その上、θをθ＋пで置換した場合でも、式（１７）は不変である。

したがってθは、ｃｏｓθ＞０となるように−п／２と−п／２の間で選択できる。
こうして、ｃｏｓθとｃｏｓ２θの両方が正であると仮定することができ、これは、−п／４≦θ≦п／４を仮定することと等価である。

角度θおよびφの明示的表現が不要であるという点に留意されたい。更新（１７）をインプリメンテーションするにはｃｏｓｈφ、ｓｉｎφ、ｃｏｓθおよびｓｉｎθの計算のみが必要である。単純な三角法により（２２）の以下の解が直ちに得られる：

この段階で、ギブンス回転Ｇ（θ，ｉ，ｊ）およびハイパボリック回転Ｈ（φ，ｉ，ｊ）を用いて所与のインデックスの対１≦ｉ＜ｊ≦ＮについてＲ’_ｋ［ｉ，ｊ］を全体的に最小化するための方法が知られる。これらの最適な角度θおよびφがひとたび計算されたならば、セットＲの全ての行列は（１７）により更新され、これは、Ｒ_ｋの行または列ｉおよびｊのわずか二つの線形結合を計算することによって達成することができる。

Ｎ×Ｎ恒等行列で初期化された混合行列は、ＡＲ_ｋＡ^Ｔ＝Ａ’Ｒ’_ｋＡ’^Ｔを保つために以下の式によって更新される。

更新の後、Ａの列のノルムはもはや等しくなく、行列Ｄ（ｉ，ｊ，λ）で再正規化され得る。次節では、この正規化を各ステップで行なう必要のないことが説明されている。Ｒ_ｋおよびＡの更新は、ヤコビアルゴリズム［８］の場合と同様に、１≦ｉ＜ｊ≦Ｎで各々のインデックスの対（ｉ，ｊ）について反復される。Ｎ（Ｎ−１）／２つのインデックスの対にわたる完全な更新は、スイープと呼ばれる。古典的なヤコビアルゴリズムと類似して、完全な収束に至るまで複数のスイープが一般に必要である。

Ｖ．Ｂ． インプリメンテーションの詳細
その他の回転シーケンスが考えられる。ギブンス回転は直交であるため、これらの回転はハイパボリック回転に比べ発生させる丸めノイズが少ない。したがって、例えば非直交ハイパボリック回転を使用する前に、純粋なギブンス回転更新の複数のスイープ（すなわちＪＡＤＥスイープ）により結合対角化を開始することができると思われる。この事前準備は、たとえ悪コンディションの混合行列についてさえ（最高１０００００の条件数）、我々の数値的実験において不要であった。数値的精度の喪失は恐らく、大きい角度φでのハイパボリック回転に結びつけられる：これらの事象も同様に監視可能であると思われる。

ＪＡＤＥまたはＳＯＢＩ予備白色化は不要であるが、それでもＲが正値定符号行列を含む場合、数値精度を改善するために使用することができる。ＢＳＳアプリケーションにおいては、ソースの数がセンサーの数よりも少ないという予備的知識が知られている場合、わずかな主成分の部分空間へのデータの削減を明らかに計算することができる。

式（２３）、（２４）、（２５）、（２６）は、（２２）の考えられる唯一の反転公式（ｉｎｖｅｒｓｉｏｎ ｆｏｒｍｕｌａｅ）ではない。これらは、ｖ１およびｖ２が０まで収束する時のその優れた数値的精度のために選択された。より良い選択肢が存在する可能性もある。収束に近い場合、行列ＣおよびＣ^ＴＣは特異点に近い。このとき、（Ｃ^ＴＣ，Ｊ）の一般固有分解は数値的に困難であり得、時として標的一般固有ベクトルはもはや、より小さい正の一般固有値に対応せず、非常に小さい負の固有値に対応する。より大きな正の固有値が誤って選択されるのを回避するため、より小さい正の固有値の代りにメジアン一般固有値に付随する固有ベクトルを選択することが好ましい。

以下のシミュレーションでインプリメンテーションが行われる。固有値符号が正しくない場合（＋、＋、−）に更新を放棄することまたは標的固有ベクトルのみを計算するための具体的手順を開発することといったその他の解決法を調査することができると思われる。

シミュレーションにおいては、Ａの列の正規化が各ステップにおいて必要でないが、一回のスイープにつき一回のみ実施してよいということが観察されている。それは、以下の式で、右側では対角行列Ｄ_ＡをＡに乗算し、両側からＤ^−１ _ＡをＲ_ｋに乗算する（式（２）および（３）と同様）ことによって簡単に達成される。

なお式中、ａ_ｎはＡのｎ番目の列であり、｜｜ａ_ｎ｜｜はそのユークリッド・ノルムである。

各スイープにおけるこの正規化（２ＫＮ^２）のコストは、Ｒ_ｋ（２ＫＮ^３）の更新（１７）に比較して無視できるものである。シミュレーションを考慮すると、本発明はなおも、列のノルムのいかなる正規化もなく収束する。

非対角項のノルムを監視することおよびその最小値での安定化を検出することといった複数の停止基準を想像することができる。例えば当該ケースにおいては、スイープの全ての三角法およびハイパボリックサインが非常に小さい場合、反復を停止することが選択される（Ｍａｔｌａｂの二倍精度で正確に同時対数化することが可能な行列および計算のためには、１０^−２０〜１０^−２５の間の閾値が適切であると思われる：この値は、正確に結合対角化することが可能でないセットについては１０^−１０〜１０^−２０の間で増大されなければならない）。

Ｖ．Ｃ． Ｊ−Ｄｉアルゴリズムのまとめ
次節の数値的シミュレーションはＪ−Ｄｉの以下のバージョンをインプリメンテーションする。
要件：Ｋ個の対称Ｎ×Ｎ行列Ｒ１，…ＲＫおよび閾値Τ

該アルゴリズムの計算コストは、主としてインデックスｉおよびｊの各対についてＲ_ｋ行列を更新するコスト、すなわち（１７）のコストである。ＧおよびＨの疎性とＲ_ｋの対称性を用いることで、この更新のコストを４ＫＮフロップ（ｆｌｏｐｓ）に削減することができる。Ｎ（Ｎ−１）／２個のインデックスの対の全体にわたる一回のスイープのためのコストは、このときおおよそ２ＫＮ^３である。収束に至るまでに行なうべきスイープの数は、次節で記述する通り、手元のＮＯＪＤ問題のサイズと難度により左右される。

Ｖ．Ｄ． Ａｆｓａｒｉアルゴリズムとの簡潔な比較
Ｊ−Ｄｉと［７］で提案されている四つのアルゴリズムとの間にあるいくつかの類似性は、強調するに値する。
−両方の方法が共に、平面変換（ｐｌａｎａｒ ｔｒａｎｓｆｏｒｍａｔｉｏｎ）の積として、対角化手段を構築する、すなわち、ＱＲＪ１ＤおよびＱＲＪ２Ｄにおけるギブンス回転および三角行列、またはＬＵＪ１ＤおよびＬＵＪ２Ｄ方法における三角行列のみ、そしてＪ−Ｄｉの場合のギブンス回転およびハイパボリック回転。全てのケースにおいて、これらの変換の行列式は１であり、それらの積は特殊線形群にまたがる。
−ＬＵＪ１Ｄ、ＬＵＪ２Ｄ、ＱＲＪ１ＤおよびＱＲＪ２Ｄは、グローバルコスト関数を最小化し、Ｊ−Ｄｉは、各々の現行のインデックスの対（ｉ；ｊ）のローカル費用関数を最小化する。
−各々のインデックスの対の二つの角度はＪ−Ｄｉにより同時に最適化され、一方で、二つの最適化はＬＵＪ１Ｄ、ＬＵＪ２Ｄ、ＱＲＪ１ＤおよびＱＲＪ２Ｄにおいて独立したものである。
−（必然的に）限定されたシミュレーションを考慮すると、ＬＵＪ１Ｄ、ＬＵＪ２Ｄ、ＱＲＪ１ＤおよびＱＲＪ２Ｄは、一次的に収束し、Ｊ−Ｄｉは二次的に収束する。
−スイープ一回あたりの計算コストは（二つの（ｉ；ｊ）−更新をグループにした場合）同一である。

ＬＵＪ１Ｄは、特定のシミュレーションされたコンテクスト内で比較的速いと思われることから、次節において比較のために選択される。

ＶＩ． 数値的シミュレーションについての性能評価
ＦＦＤｉａｇ［４］およびＬＵＪ１Ｄ［７］という二つの終りのＮＯＪＤアルゴリズムがＲの任意の行列の正値符号性（ｐｏｓｉｔｉｖｅ ｄｅｆｉｎｉｔｅｎｅｓｓ）を必要としないことから、これら二つのアルゴリズムと本発明に係る方法（すなわちＪ−Ｄｉ）を比較する。

ＦＦＤｉａｇアルゴリズムの１の反復（以下でスイープと呼ぶ）は、一つのＪ−ＤｉまたはＬＵＪ１Ｄスイープと比べてほぼ同じ複雑性を有する（Ｒを更新するために２ＫＮ^３）。したがって、唯一比較されるのは、収束に必要なスイープ数だけである。

それでも、一つのスイープ内で演算は異なる形で構造化されている。すなわち、Ｊ−ＤｉおよびＬＵＪ１ＤはＲの非常に疎な更新を必要とし、一方でＦＦＤｉａｇは、完全な高密度の更新のみを必要とする。

この差がおそらく、Ｍａｔｌａｂのようなインタプリタ型言語にインプリメンテーションされた場合のＪ−ＤｉまたはＬＵＪ１Ｄスイープに比べ、ＦＦＤｉａｇのスイープをわずかに速いものにしている。ただし、例えば、Ｃ言語によるインプリメンテーションは類似の実行時間を示すはずである。

二つの異なる性能の基準が定義され、最初に正確に対角化可能な行列のセットについて測定され、第二にほぼ対角化可能な行列のセットについて測定される。シミュレーションは、Ｍａｔｌａｂの二倍精度で実行される。

ＶＩ．Ａ． 非直交結合対角化についての性能の基準
以下では結合対角化の質を測定する二つの方法が用いられる。最初の基準は、Ｒの非対角エントリの二乗の和として定義される非対角ノルムとＲの対角エントリの二乗の和として定義される対角ノルムの比率を計算することから成る。この非標準正規化は、対角ノルムが最初のスイープの間一定でないことから有用である。

第二の基準は、混合行列の推定の質と、対応する分離性能を査定するのに必要である。実際のところ、下記の〔数６５〕の推定された混合行列のインバースは、Ａによって混合されたソースを分離するのに（付加ノイズが無い場合に）最適な空間フィルタである。なお、本明細書においては、明細書の記載として使用できる文字の制約を考慮し、

上記〔数６５〕を、テキスト上、「Ａ：」と表記する。

積Ａ：^−１Ａが対角行列Ｄの置換行列Π倍に近いほど、分離は優れたものとなる。したがって、

と、置換行列と対角行列の最も近い積ΠＤとの間の距離を計算することによって、分離の質を測定するのが自然である。

ＢＳＳにおいては、たとえそれがΠＤまでの完全な距離でなくても、

によって定義される性能Ｉ（Ｐ）のインデックスを使用することが標準的である。正規化項１／２Ｎ（Ｎ−１）は、Ｐが対角エントリが１に等しく、非対角エントリがαに等しい行列である場合に、０≦α＜１においてＩ（Ｐ）がαに等しくなるような形で選択された。したがって、Ｉ（Ｐ）は、分離されたソースにおける相対的な干渉振幅（ａｍｐｌｉｔｕｄｅ ｏｆ ｉｎｔｅｒｆｅｒｅｎｃｅｓ）の桁（ｏｒｄｅｒ ｏｆ ｍａｇｎｉｔｕｄｅ）を提供する（クロストークインデックス（ｃｒｏｓｓ−ｔａｌｋ ｉｎｄｅｘ））。

ＶＩ．Ｂ． 正確に対角化可能な行列セット
混合行列エントリは独立して正規分布しており、その行列式は１に設定され、その列は、等式（３）を用いて同じノルムに正規化される。これらのランダム混合行列は一般的に、以下の図面で示されている通り、１０^４または１０^５未満の条件数を有し、それより多いことはまれである。Ｊ−Ｄｉアルゴリズムのより特定的な妥当性確認は、より悪条件の混合行列が関与するアプリケーションのために行なわれるべきである。Ｄ_ｋの行列の対角エントリは、同じく、標準的な正規分布から選択される。このようなＤ_ｋ生成方法が、ＫおよびＮのみに左右される一意性のモジュラス（ｍｏｄｕｌｕｓ ｏｆ ｕｎｉｑｕｅｎｅｓｓ）の分布を決定するという点に留意されたい。Ｒ_ｋはＲ_ｋ＝ＡＤ_ｋＡ^Ｔによって計算される。

１）二次収束：
Ｊ−Ｄｉ、ＦＦＤｉａｇおよびＬＵＪ１Ｄアルゴリズムの対角性比は、収束に至るまで各スイープ後に計算され、図２には、Ｎ＝Ｋ＝５０および５０の独立した試行についてプロットされている。この場合、一意性のモジュラスは小さく、ほぼ常に０．４〜０．６の間にある。

分離基準は、同じ条件下で図３にプロットされている。

これらの結果は、ＦＦＤｉａｇおよびＪ−Ｄｉアルゴリズムが（少なくとも）二次収束を有すると思われるが、Ｊ−Ｄｉはより少ないスイープしか必要としないことを示している。最終的な一次収束にも関わらず、ＬＵＪ１Ｄは中間の性能を示す。

２）Ｊ−Ｄｉ収束速度に対する行列およびセットサイズの影響
異なる混合行列および異なるセットサイズ（Ｎ＝Ｋ＝５、１０、２０、５０、１００および２００）について、Ｊ−ＤＩの収束に達するのに必要なスイープ数が計数され、１００回の独立した試行についてのその分布が、下表に示されている：

Ｊ−Ｄｉスイープ数は行列サイズと共に増大するものの、大きい行列（Ｎ＝２００）の大きなセット（Ｋ＝２００）においても、ほぼ１０未満にとどまることがわかる。

ＶＩ．Ｃ． 正確に対角化できない行列セット：ブラインドソース分離に対する応用

Ｊ−Ｄｉ、ＦＦＤｉａｇおよびＬＵＪ１Ｄ非直交結合対角化アルゴリズムの性能を、正確に対角化できない行列セットについて比較する。

（共分散行列のように）正でもなく負でもない行列セットについて本発明に係る方法（Ｊ−Ｄｉ）を評価するため、四次キュミュラントスライスセットが選択される。

古典的ＢＢＳモデルすなわちｘ＝Ａｓ＋ｎを考慮する。なお式中、ｘ、ｓおよびｎはそれぞれに観察された混合体のＮ次元ベクトル、ソースおよび付加的ガウスノイズ（ａｄｄｉｔｉｖｅ Ｇａｕｓｓｉａｎ ｎｏｉｓｅ）を表わし、ＡはＮ×Ｎの混合行列である。考慮対象の行列セットＲは、

によって定義されるサイズＮ×ＮのＫ＝Ｎ（Ｎ＋１）／２個の行列Ｒ_ｍｎを含んでいる。ここで、全てのｍおよびｎが１≦ｍ≦ｎ≦Ｎであり、式中ｘ_ｎはｘのｎ番目のエントリを表わす。

Ｎ＝２０ソースであり、したがってＫ＝Ｎ（Ｎ＋１）／２＝２１０の２０×２０サイズのキュムラントスライス（行列）であると考えられる。ソース信号のサンプルは独立しており、−√３と√３の間で一様分布して非ガウス、ゼロ平均（ｚｅｒｏ−ｍｅａｎ）かつユニティバリアンスな（ｕｎｉｔｙ−ｖａｒｉａｎｃｅ）信号を提供する。ゼロ平均ガウス付加ノイズ（ｚｅｒｏ−ｍｅａｎ Ｇａｕｓｓｉａｎ ｎｏｉｓｅ）σの標準偏差は、まずゼロに設定され、第二にσ＝０．１に設定される。キュムラントは、１０００００個の独立したサンプルについて推定される。これらのケースにおいて、行列数Ｋは充分大きいため、一意性のモジュラス（ｍｏｄｕｌｕｓ ｏｆ ｕｎｉｑｕｅｎｅｓｓ）は一般的に０．３〜０．６の間という小さいものとなる。このパラメータがシミュレーションの性能に影響を及ぼしそうにない理由がこれによって説明されており、したがってこれは提示されない。

収束速度
対角性比および分離基準の収束は、図４にノイズを伴うスイープ数の関数としてプロットされている。

全てのケースにおいて、本発明に係る方法（Ｊ−Ｄｉ）は、［７］の場合と同様に挙動すると思われるＦＦＤｉａｇおよびＬＵＪ１Ｄよりも速く収束する。Ｊ−Ｄｉは、ＢＳＳアプリケーションにとって重要である最終的な分離基準の観点からみて、ＦＦＤｉａｇおよびＬＵＪ１Ｄよりもわずかに性能が優れている。

これらのシミュレーションは同様に、ＦＦＤｉａｇは二次でありＬＵＪ１Ｄは一次である異なる収束速度が、行列セットが正確に対角化可能でない場合の類似の性能と適合可能であることをも示している。

ＶＩＩ． ＷＡＮＧ．ＬＩＵおよびＺＨＡＮＧのアルゴリズムとの比較
Ｗａｎｇ、ＬｉｕおよびＺｈａｎｇ［６］（以下ＤＮＪＤ）のアルゴリズムとの比較について以下で論述する。提案されているアルゴリズムＤＮＪＤは、

によって定義される行列Ｊ（θ，φ，ｉ，ｊ）の連続的乗算によって混合行列を構築する。

この定義は明らかに、ギブンス回転およびＪＡＤＥのヤコビのようなアルゴリズムを一般化する試みである。ＤＮＪＤアルゴリズムが本発明と非常に異なるものであることは明白である。例えば、ＤＮＪＤは、本発明に係る方法において用いられる行列式制約もハイパボリック三角法も使用しない。

ＤＮＪＤにおける主たる困難性は、等式（１６）の反転にある。Ｒ_ｋ行列が正確に対角化できない場合に、非対角エントリの実際の最小化を保証しない近似式を使用することが必要となり、正確な解は提供されない。したがって、（各インデックスの対について）各反復毎に非対角エントリの進化を計算し、非対角エントリが増大した場合更新を放棄することが必要である。残念なことに、この確認の計算コスト４ＫＮは、アルゴリズム自体のコストと同等である（第ＩＩＩ節の終りの５３０３頁を参照のこと）。したがって、各スイープのコストは二倍となる。さらに、二つ一組の角度（θ，φ）が放棄される毎に、数多くの計算（８ＫＮフロップ）が費用関数の改善無く行なわれる。この無駄が、Ｊ−Ｄｉに比べてＤＮＪＤが用いるスイープ数の多さを説明し得ると思われる（［６］中の５３０３項、図１）（表１）。

要するに、Ｊ−Ｄｉはスイープ毎の複雑性の観点から見てＤＮＪＤよりも二倍以上優れた性能を示す。

ＶＩＩＩ． 結論
非直交結合対角化の新しい方法（Ｊ−Ｄｉ）が提示された。一つの実施態様において、この方法は、以下のステップを含む：
−特殊な直交群の中で構築する代わりに、列が等ノルムであり、また、特殊線形群に属する混合行列を構築するステップ；
−ギブンス回転だけでなく、ギブンス回転、ハイパボリック回転および対角行列を乗算するステップ。

この新しい方法は、数多くの変形形態すなわち基本行列の乗算の次数、停止テスト、予備白色化の有無などを許容する。そのロバスト性および削減された複雑性により、標準コンピュータ上で最高Ｎ＝５００およびＫ＝６００までの大規模な問題を扱うことが可能になる。

１０ 処理手段
２０ 表示手段
３０ 入力手段

［１］Ｃｏｍｐａｒａｉｓｏｎ ｄｅ ｍｅｔｈｏｄｅｓ ｄｅ ｓｅｐａｒａｔｉｏｎｓ ｄｅ ｓｏｕｒｃｅｓ，Ｓｏｕｌｏｕｍｉａｃ Ａ．，Ｃａｒｄｏｓｏ Ｊ．−Ｆ．，Ｐｒｏｃｅｅｄｉｎｇｓ ｏｆ ＧＲＥＴＳＩ １９９１，Ｊｕａｎ ｌｅｓ Ｐｉｎｓ，Ｓｅｐｔｅｍｂｒｅ １９９１． ［１ｂ］Ｂｌｉｎｄ ｂｅａｍｆｏｒｍｉｎｇ ｆｏｒ ｎｏｎ Ｇａｕｓｓｉａｎ ｓｉｇｎａｌｓ，Ｊｅａｎ−Ｆｒａｎｃｏｉｓ Ｃａｒｄｏｓｏ ａｎｄ Ａｎｔｏｉｎｅ Ｓｏｕｌｏｕｍｉａｃ，ｉｎ ＩＥＥ Ｐｒｏｃｅｅｄｉｎｇｓ−Ｆ，１４０（６）：３６２−３７０，Ｄｅｃｅｍｂｅｒ １９９３． ［２］Ａ ｂｌｉｎｄ ｓｏｕｒｃｅ ｓｅｐａｒａｔｉｏｎ ｔｅｃｈｎｉｑｕｅ ｕｓｉｎｇ ｓｅｃｏｎｄ−ｏｒｄｅｒ ｓｔａｔｉｓｔｉｃｓ，Ｂｅｌｏｕｃｈｒａｎｉ Ａ．，ＡｂｅｄＭｅｒａｉｍ Ｋ．，Ｃａｒｄｏｓｏ Ｊ．−Ｆ．，Ｍｏｕｌｉｎｅｓ Ｅ．，ＩＥＥＥ ＴＲＡＮＳＡＣＴＩＯＮＳ ＯＮ ＳＩＧＮＡＬ ＰＲＯＣＥＳＳＩＮＧ ４５（２）：４３４−４４４，ＦＥＢ １９９７． ［３］Ｊｏｉｎｔ Ａｐｐｒｏｘｉｍａｔｅ Ｄｉａｇｏｎａｌｉｚａｔｉｏｎ ｏｆ Ｐｏｓｉｔｉｖｅ Ｄｅｆｉｎｉｔｅ Ｈｅｒｍｉｔｉａｎ Ｍａｔｒｉｃｅｓ，Ｐｈａｍ Ｄ．Ｔ．，ＳＩＡＭ Ｊ．ｏｎ Ｍａｔｒｉｘ Ａｎａｌ．ａｎｄ Ａｐｐｌ．，２２（４）：１１３６−１１５２，２００１． ［４］Ａ Ｆａｓｔ Ａｌｇｏｒｉｔｈｍ ｆｏｒ Ｊｏｉｎｔ Ｄｉａｇｏｎａｌｉｚａｔｉｏｎ ｗｉｔｈ Ｎｏｎ−ｏｒｔｈｏｇｏｎａｌ Ｔｒａｎｓｆｏｒｍａｔｉｏｎｓ ａｎｄ ｉｔｓ Ａｐｐｌｉｃａｔｉｏｎ ｔｏ Ｂｌｉｎｄ Ｓｏｕｒｃｅ Ｓｅｐａｒａｔｉｏｎ，Ｚｉｅｈｅ Ａ．，Ｌａｓｋｏｖ Ｐ．，Ｎｏｌｔｅ Ｇ．，Ｍｕｌｌｅｒ Ｋ．−Ｒ．，Ｊｏｕｒｎａｌ ｏｆ Ｍａｃｈｉｎｅ Ｌｅａｒｎｉｎｇ Ｒｅｓｅａｒｃｈ ５（２００４）７７７−８００． ［５］Ｑｕａｄｒａｔｉｃ Ｏｐｔｉｍｉｚａｔｉｏｎ ｆｏｒ Ｓｉｍｕｌｔａｎｅｏｕｓ Ｍａｔｒｉｘ Ｄｉａｇｏｎａｌｉｚａｔｉｏｎ，Ｖｏｌｌｇｒａｆ Ｒ．，Ｏｂｅｒｍａｙｅｒ Ｋ．，ＩＥＥＥ Ｔｒａｎｓ，ｏｎ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｉｎｇ，ｖｏｌ．５４，ｎｏ．９，ｐｐ．３２７０−３２７８，Ｓｅｐｔｅｍｂｅｒ ２００６． ［６］Ｎｏｎｏｒｔｈｏｇｏｎａｌ Ｊｏｉｎｔ Ｄｉａｇｏｎａｌｉｚａｔｉｏｎ Ａｌｇｏｒｉｔｈｍ Ｂａｓｅｄ ｏｎ Ｔｒｉｇｏｎｏｍｅｔｒｉｃ Ｐａｒａｍｅｔｅｒｉｚａｔｉｏｎ，Ｗａｎｇ Ｆ．，Ｌｉｕ Ｚ．，Ｚｈａｎｇ Ｊ．，ＩＥＥＥ Ｔｒａｎ．ｏｎ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｉｎｇ，ｖｏｌ．５５，ｎｏ．１１，ｐｐ．５２９９−５３０８，Ｎｏｖｅｍｂｅｒ ２００７． ［７］Ｓｉｍｐｌｅ ＬＵ ａｎｄ ＱＲ Ｂａｓｅｄ Ｎｏｎ−Ｏｒｔｈｏｇｏｎａｌ Ｍａｔｒｉｘ Ｊｏｉｎｔ Ｄｉａｇｏｎａｌｉｚａｔｉｏｎ，Ｂ．ＡＦＳＡＲＩ，Ｉｎｄｅｐｅｎｄｅｎｔ Ｃｏｍｐｏｎｅｎｔ Ａｎａｌｙｓｉｓ ａｎｄ Ｂｌｉｎｄ Ｓｏｕｒｃｅ Ｓｅｐａｒａｔｉｏｎ，Ｌｅｃｔｕｒｅ Ｎｏｔｅｓ ｉｎ Ｃｏｍｐｕｔｅｒ Ｓｃｉｅｎｃｅ，ｖｏｌ．３８８９，ｐｐ．１−７，２００６． ［８］Ｍａｔｒｉｘ Ｃｏｍｐｕｔａｔｉｏｎｓ，Ｇ．Ｈ．ＧＯＬＵＢ，Ｃ．Ｆ．ＶＡＮ ＬＯＡＮ，Ｂａｌｔｉｍｏｒｅ，Ｊｏｈｎｓ Ｈｏｐｋｉｎｓ Ｕｎｉｖｅｒｓｉｔｙ Ｐｒｅｓｓ，３ｒｄ ｅｄｉｔｉｏｎ，１９９６． ［９］Ｗｈａｔ Ｃａｎ Ｍａｋｅ Ｊｏｉｎｔ Ｄｉａｇｏｎａｌｉｚａｔｉｏｎ Ｄｉｆｆｉｃｕｌｔ？，Ｂ．ＡＦＳＡＲＩ，ＩＣＡＳＳＰ ２００７，Ｖｏｌ．３，ｐｐ．ＩＩＩ−１３７７−１３８０，１５−２０ Ａｐｒｉｌ ２００７． ［１０］Ｓｅｎｓｉｔｉｖｉｔｙ Ａｎａｌｙｓｉｓ ｆｏｒ ｔｈｅ ｐｒｏｂｌｅｍ ｏｆ Ｎｏｎ−Ｏｒｔｈｏｇｏｎａｌ Ｍａｔｒｉｘ Ｊｏｉｎｔ Ｄｉａｇｏｎａｌｉｚａｔｉｏｎ，Ｂ．ＡＦＳＡＲＩ，Ｓｕｂｍｉｔｔｅｄ ｔｏ ＳＩＡＭ Ｊ．ｏｎ Ｍａｔｒｉｘ Ｃｏｍｐｕｔａｔｉｏｎｓ． ［１１］Ｎｏｎ−Ｏｒｔｈｏｇｏｎａｌ Ｊｏｉｎｔ Ｄｉａｇｏｎａｌｉｚａｔｉｏｎ ｉｎ ｔｈｅ Ｌｅａｓｔ−Ｓｑｕａｒｅｓ Ｓｅｎｓｅ ｗｉｔｈ Ａｐｐｌｉｃａｔｉｏｎ ｉｎ Ｂｌｉｎｄ Ｓｏｕｒｃｅ Ｓｅｐａｒａｔｉｏｎ，Ａ．ＹＥＲＥＤＯＲ，ＩＥＥＥ Ｔｒａｎｓ．ｏｎ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｉｎｇ，ｖｏｌ．５０，ｎｏ．７，ｐｐ．１５４５−１５５３，Ｊｕｌｙ ２００２． ［１２］Ｓｕｒ ｌａ ｄｉａｇｏｎａｌｉｓａｔｉｏｎ ｃｏｎｊｏｉｎｔｅ ａｐｐｒｏｃｈｅｅ ｐａｒ ｕｎ ｃｒｉｔｅｒｅ ｄｅｓ ｍｏｉｎｄｒｅｓ ｃａｒｒｅｓ，Ｓ．ＤＥＧＥＲＩＮＥ，ｉｎ Ｐｒｏｃｅｅｄｉｎｇｓ ＸＶＩＩＩ Ｃｏｌｌｏｑｕｅ Ｇｒｅｔｓｉ ２００１，ｎ２，ｐｐ．３１１−３１４，Ｔｏｕｌｏｕｓｅ，Ｆｒａｎｃｅ，Ｓｅｐｔ．２００１．

Claims (10)

1. 少なくとも二つのセンサーにより記録された混合された信号（ｘ_１（ｔ），．．．，ｘ_Ｎ（ｔ））を複数の成分信号（ｓ_１（ｔ），．．．，ｓ_Ｍ（ｔ））へと分離するための方法であり、各混合された信号Ｘ_ｎ（ｔ）は独立したソースに由来する成分信号（ｓ_１（ｔ），．．．，ｓ_Ｍ（ｔ））の線形結合に対応しており、ソースに由来する成分信号（ｓ_１（ｔ），．．．，ｓ_Ｍ（ｔ））のベクトルを混合行列（Ａ）に乗じた積として各混合された信号を記すことが可能であり、
   該方法は、さらに結合対角化ステップ（３００）を含み、このステップでは一つまたは複数のギブンス回転行列（Ｇ（θ））および一つまたは複数のハイパボリック回転行列（Ｈ（φ））に対応する複数の回転行列の積から混合行列（Ａ）が推定されることを特徴とする方法。
2. ギブンス回転行列が、θの三角関数に対応する要素を含み、ハイパボリック回転行列がφのハイパボリック関数に対応する要素を含むことを特徴とする、請求項１に記載の方法。
3. 対角化ステップには、一方にはギブンス回転行列とハイパボリック回転行列の積、そして他方にはこの積の転置行列による推定された行列の族
   

   

   の乗算と、この乗算の結果としてのｎ次の行列の非対角項を最小化しうるθ_ｎおよびφ_ｎ値の推定
   とが含まれることを特徴とする、請求項２に記載の方法。
4. 結果として得られたｎ次の行列に対して乗算ステップが反復されることを特徴とする、請求項３に記載の方法。
5. ｃｏｓθ_ｎおよびｃｏｓｈφ_ｎがそれぞれ１に向かうか、対応するサインがゼロに向かうギブンス角度θ_ｎおよびハイパボリック角度φ_ｎになるまで乗算ステップが反復されることを特徴とする、請求項４に記載の方法。
6. 推定された行列が、混合行列と任意の対角行列による転置された混合行列との積に等しく、これらが例えば、キュムラント行列、または異なる周派数帯域内の共分散行列、または異なる時間推移を有する相互相関行列、または異なる時間的間隔にわたって推定された共分散行列であることを特徴とする、請求項１〜５のいずれか一つに記載の方法。
7. ギブンス回転行列が、対角線上のｃｏｓθ項、対角線より下のｓｉｎθ項および対角線より上の−ｓｉｎθ項を含み、ハイパボリック回転行列が対角線上のｃｏｓｈφ要素、対角線より下のｓｉｎｈφおよび対角線より上のｓｉｎｈφ要素を含むことを特徴とする、請求項１〜６のいずれか一つに記載の方法。
8. ギブンス回転行列が、
   

   

   タイプのものであり、ハイパボリック回転行列が、
   

   

   タイプのものであることを特徴とする、請求項１〜７のいずれか一つに記載の方法。
9. 少なくとも二つのセンサーによって記録された混合された信号（ｘ_１（ｔ），．．．，ｘ_Ｎ（ｔ））を複数の成分信号（ｓ_１（ｔ），．．．，ｓ_Ｍ（ｔ））へと分離するための装置であり、各混合された信号Ｘ_ｎ（ｔ）は独立したソースに由来する成分信号（ｓ_１（ｔ），．．．，ｓ_Ｍ（ｔ））の線形結合に対応し、これらのソースに由来する成分信号のベクトルを混合行列（Ａ）に乗じた積として各混合された信号を記すことができるものであり、
   該装置は、さらに結合対角化用手段を備え、該手段では一つまたは複数のギブンス回転行列および一つまたは複数のハイパボリック回転行列に対応する複数の回転行列の積から混合行列Ａが推定されることを特徴とする、装置。
10. コンピュータに請求項１〜８のいずれか一つに記載の方法を実行させるためのプログラムを記録したコンピュータ読み取り可能な記録媒体。
    

Images