(57)【要約】  本発明は、独立成分の瞬間線形混合体で構成された信号セットを分解する方法に関するものである。該方法は、一つまたは複数の三角法回転行列と一つまたは複数のハイパボリック回転行列との積を用いることによって、混合行列Ａを推定することからなるステップ（２００）を備えている。【選択図】図１
2011509037.tif 000002
【特許請求の範囲】
【請求項１】
少なくとも二つのセンサーにより記録された混合された信号（ｘ１（ｔ），．．．，ｘＮ（ｔ））を複数の成分信号（ｓ１（ｔ），．．．，ｓＭ（ｔ））へと分離するための方法であり、各混合された信号Ｘｎ（ｔ）は独立したソースに由来する成分信号（ｓ１（ｔ），．．．，ｓＭ（ｔ））の線形結合に対応しており、ソースに由来する成分信号（ｓ１（ｔ），．．．，ｓＭ（ｔ））のベクトルを混合行列（Ａ）に乗じた積として各混合された信号を記すことが可能であり、
該方法は、さらに結合対角化ステップ（３００）を含み、このステップでは一つまたは複数のギブンス回転行列（Ｇ（θ））および一つまたは複数のハイパボリック回転行列（Ｈ（φ））に対応する複数の回転行列の積から混合行列（Ａ）が推定されることを特徴とする方法。
【請求項２】
ギブンス回転行列が、θの三角関数に対応する要素を含み、ハイパボリック回転行列がφのハイパボリック関数に対応する要素を含むことを特徴とする、請求項１に記載の方法。
【請求項３】
対角化ステップには、一方にはギブンス回転行列とハイパボリック回転行列の積、そして他方にはこの積の転置行列による推定された行列の族
【数１】
2011509037.tif 000073
の乗算と、この乗算の結果としてのｎ次の行列の非対角項を最小化しうるθｎおよびφｎ値の推定
とが含まれることを特徴とする、請求項２に記載の方法。
【請求項４】
結果として得られたｎ次の行列に対して乗算ステップが反復されることを特徴とする、請求項３に記載の方法。
【請求項５】
ｃｏｓθｎおよびｃｏｓｈφｎがそれぞれ１に向かうか、対応するサインがゼロに向かうギブンス角度θｎおよびハイパボリック角度φｎになるまで乗算ステップが反復されることを特徴とする、請求項４に記載の方法。
【請求項６】
推定された行列が、混合行列と任意の対角行列による転置された混合行列との積に等しく、これらが例えば、キュムラント行列、または異なる周派数帯域内の共分散行列、または異なる時間推移を有する相互相関行列、または異なる時間的間隔にわたって推定された共分散行列であることを特徴とする、請求項１〜５のいずれか一つに記載の方法。
【請求項７】
ギブンス回転行列が、対角線上のｃｏｓθ項、対角線より下のｓｉｎθ項および対角線より上の−ｓｉｎθ項を含み、ハイパボリック回転行列が対角線上のｃｏｓｈφ要素、対角線より下のｓｉｎｈφおよび対角線より上のｓｉｎｈφ要素を含むことを特徴とする、請求項１〜６のいずれか一つに記載の方法。
【請求項８】
ギブンス回転行列が、
【数２】
2011509037.tif 000074
タイプのものであり、ハイパボリック回転行列が、
【数３】
2011509037.tif 000075
タイプのものであることを特徴とする、請求項１〜７のいずれか一つに記載の方法。
【請求項９】
少なくとも二つのセンサーによって記録された混合された信号（ｘ１（ｔ），．．．，ｘＮ（ｔ））を複数の成分信号（ｓ１（ｔ），．．．，ｓＭ（ｔ））へと分離するための装置であり、各混合された信号Ｘｎ（ｔ）は独立したソースに由来する成分信号（ｓ１（ｔ），．．．，ｓＭ（ｔ））の線形結合に対応し、これらのソースに由来する成分信号のベクトルを混合行列（Ａ）に乗じた積として各混合された信号を記すことができるものであり、
該装置は、さらに結合対角化用手段を備え、該手段では一つまたは複数のギブンス回転行列および一つまたは複数のハイパボリック回転行列に対応する複数の回転行列の積から混合行列Ａが推定されることを特徴とする、装置。
【請求項１０】
コンピュータ内で使用することのできる媒体に記録されたプログラムコード命令を含むコンピュータプログラム製品であり、請求項１〜８のいずれか一つに記載の方法を適用するための命令を含むことを特徴とする、コンピュータプログラム製品。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、信号処理の一般的な技術分野に関するものであり、より具体的には、ブラインドソース分離を用いてソース信号の混合体からソース信号を分離する一般的な技術分野に関するものである。
【０００２】
本発明は、数多くの分野に応用できる。特に、脳波図、無線通信の分野における信号処理、さらには質量スペクトル処理にこれを応用することができる。
【背景技術】
【０００３】
いずれも所与の瞬間にセンサーネットワークによって記録された複数の混合体から統計的に独立した複数のソースを分離すること（空間分散）からなるソース分離の問題は、信号処理において比較的新しい問題である。
【０００４】
観察された物理的現象が、複数の独立した基本的現象の積み重ねであり、この積み重ねが瞬間線形結合（または混合体）によって表現されることは、頻繁にあることである。
【０００５】
この現象の観察に由来する測定データは、このとき、各々の独立した基本的現象からの寄与の瞬間線形混合体である。
【０００６】
この場合、複数のセンサーを設置することによって、これらの寄与の複数の異なる瞬間線形混合体を測定することが可能である。例えば、脳波検査においては、患者の頭皮上にさまざまな電極が設置される。これらの電極によって、患者の脳活動を測定することができる。
【０００７】
各電極は、脳のさまざまなソースが発出する電気信号の混合体を測定し、異なる電極が異なる混合体を測定する。ソースの混合体は、ソースがその中で異なる割合でまたは異なる重みをもって出現する場合に、異なるものであると考えられている。
この場合、次のものについてのその他の情報、すなわち、
−ソースについて、および、
−これらのソース由来の信号がどのように混合されているか
は全く知られていない。
【０００８】
このとき、第一段階において混合体の重みを識別するために、そして第二段階において測定値内で各ソースの寄与を分離するために、ソース分離が使用される。こうして、混合体に妨害されることなく異なる寄与の個別分析が可能となる。
【０００９】
現在、ソース分離方法には二つのタイプが存在する。
【００１０】
第一のタイプの方法は、予備ステップ、いわゆる「白色化（ｗｈｉｔｅｎｉｎｇ）」ステップを使用し、これにより、混合体の重みを識別するために後続する計算ステップを容易にすることができる。しかしながら、この白色化ステップは、各ソースの寄与の分離時点で補正不能なエラーを誘発する。したがって、このタイプの方法はあまり正確ではない。
【００１１】
第二のタイプの方法は、白色化ステップを一切必要としない。しかしながら、この第二のタイプの方法は、効率および計算コストに関する問題を誘発する。
【発明の概要】
【発明が解決しようとする課題】
【００１２】
本発明の一つの目的は、混合された信号を複数の成分信号に分離するための方法を提案することにあり、これにより前述のタイプの方法の欠点を少なくとも一つ克服することができる。
【課題を解決するための手段】
【００１３】
この目的で、少なくとも二つのセンサーにより記録された混合された信号（ｘ１（ｔ），．．．，ｘＮ（ｔ））を複数の成分信号（ｓ１（ｔ），．．．，ｓＭ（ｔ））へと分離するための方法が提供されるものであり、該方法においては、混合された信号（ｘ１（ｔ），．．．，ｘｎ（ｔ））は独立したソースに由来するＭ個の成分信号（ｓ１（ｔ），．．．，ｓＭ（ｔ））のＮ個の線形結合に対応しており、ソースに由来する成分信号（ｓ１（ｔ）、．．．、ｓＭ（ｔ））のベクトルを混合行列（Ａ）に乗じた積として各混合された信号を記すことが可能であり、該方法はさらに結合対角化ステップを含み、このステップでは一つまたは複数のギブンス回転行列（Ｇｉｖｅｎｓ  ｒｏｔａｔｉｏｎ  ｍａｔｒｉｃｅｓ）（Ｇ（θ））、および、一つまたは複数の「ハイパボリック（双曲線）回転」行列（ｈｙｐｅｒｂｏｌｉｃ  ｒｏｔａｔｉｏｎ  ｍａｔｒｉｃｅｓ）（Ｈ（φ））に対応する複数の回転行列の積から混合行列Ａが推定される。
【００１４】
本発明の範囲内で「ハイパボリック回転行列」というのは、ハイパボリック関数（例えばｃｏｓｈ、ｓｉｎｈ）である項を含み、かつ１に等しい行列式を有する回転行列を意味する。
【００１５】
かくして、本発明によると、高速で信頼性の高い単一の最適化ステップにより混合行列を識別できる。単一のステップしか存在しないという事実によって、後述の文献［１］、［１ｂ］、［２］の場合などの２ステップの方法（白色化と、それに続く結合対角化）に比べてさらに正確な結果が得られる。実際、第一のステップの不可避なエラーを第二のステップにより補正することはできない。さらに、提案される方法では、特にロバストな、そして計算能力の観点からみてさほどコストが高くない代数アルゴリズムが使用される（これは、２×２サイズの「三角法」および「ハイパボリック」回転行列の積である）。
【００１６】
白色化ステップを備える方法は、基本直交行列（ｅｌｅｍｅｎｔａｒｙ  ｏｒｔｈｏｇｏｎａｌ  ｍａｔｒｉｃｅｓ）の積により直交混合行列（ｏｒｔｈｏｇｏｎａｌ  ｍｉｘｔｕｒｅ  ｍａｔｒｉｘ）を確立する（ギブンス回転）。
【００１７】
本発明に関して言うと、これは１に等しい行列式を有する行列を乗算することにより、１に等しい行列式を有する混合行列を構築することを可能にする。直交行列群（正の行列式＋１を有する）は、１に等しい行列式を有する行列群の中に含み入れられる。したがって、混合行列を計算するステップは、混合体が当初受けていた制約を制限することのできる本発明の方法によって容易になる。最初の「白色化」ステップはもはや不可欠ではなく、それでも「代数」法のロバスト性および効率は保たれる。
【００１８】
結果として、本発明の方法は、白色化ステップを含む方法で得られた混合行列よりも混合体の物理的現実に近い混合行列の計算を可能にする。
【００１９】
本発明に係る方法の好ましいものの非限定的な態様は、以下の通りである。
【００２０】
−ギブンス回転行列が、θの三角関数に対応する要素を含み、ハイパボリック回転行列がφのハイパボリック関数に対応する要素を含む；
【００２１】
−対角化ステップには、右にはギブンス回転行列とハイパボリック回転行列の積、そして左にはこの積の転置行列により推定行列族（Ｒ：ｋ）を乗算するためのステップ、そして、この乗算の結果としてのｎ次の行列の非対角項を最小化するためにθｎおよびφｎの値を計算するためのステップとが含まれている；
なお、本明細書においては、明細書の記載として使用できる文字の制約を考慮し、
【００２２】
【数１】
2011509037.tif 000003
【００２３】
上記〔数１〕を、テキスト上、「（Ｒ：ｋ）」と表記する。
【００２４】
−結果として得られたｎ次の行列に対して乗算ステップが反復される；
【００２５】
−ｃｏｓθおよびｃｏｓｈφがそれぞれ１に向かうギブンス角度θおよび「ハイパボリック角度」φになるまで乗算ステップが反復される；
【００２６】
−推定行列Ｒ：ｋは、（推定エラーの範囲内で）混合行列と対角行列と転置された混合行列との積に等しく（Ｒｋ＝ＡＤｋＡＴ）、これらは例えば、キュムラント（ｃｕｍｕｌａｎｔｓ）行列、または異なる周派数帯域内の共分散行列、または異なる時間推移を有する相互相関行列、または異なる時間的間隔にわたって推定された共分散行列である；
【００２７】
−ギブンス回転行列は、対角線上のｃｏｓθ要素、対角線より下のｓｉｎθ要素および対角線より上の−ｓｉｎθ要素を含み、「ハイパボリック回転」行列は対角線上のｃｏｓｈφ要素、対角線より下のｓｉｎｈφおよび対角線より上のｓｉｎｈφ要素を含む；
【００２８】
−ギブンス回転行列は、
【００２９】
【数２】
2011509037.tif 000004
【００３０】
タイプのものであり、「ハイパボリック回転」行列は、
【００３１】
【数３】
2011509037.tif 000005
【００３２】
タイプのものである。
【００３３】
本発明は同様に、少なくとも二つのセンサーによって記録された混合された信号を複数の成分信号へと分離するための装置に関するものであり、各混合された信号は独立したソースに由来する成分信号の瞬間線形結合に対応し、ソースに由来する成分信号のベクトルを混合行列に乗じた積として各混合された信号を記すことができるものであり、該装置は、さらに結合対角化用手段を備え、該手段では一つまたは複数のギブンス回転行列および一つまたは複数のハイパボリック回転行列に対応する複数の回転行列の積から混合行列が推定されることを特徴とする。
【００３４】
本発明は同様に、コンピュータで使用することのできる媒体に記録されたプログラムコード命令を含むコンピュータプログラム製品にも関するものであり、該プログラム製品は上述の方法を応用するための命令を含むことを特徴とする。
【００３５】
本発明のその他の特徴、目的および利点は以下の記述からさらに明らかとなるが、この記述は、純粋に例証であって限定的なものではなく、添付図面を参照に読まれるべきものである。
【図面の簡単な説明】
【００３６】
【図１】本発明にしたがった方法の実施態様を示す。
【図２】混合された信号を分離するための三つの異なる方法（Ｊ−Ｄｉ、ＦＦＤｉａｇおよびＬＵＪ１Ｄ）について、サイズＮ＝５０のＫ＝５０の行列のセットを用いた、５０回の独立した試行での、多くのスイープ（ｓｗｅｅｐｓ）に応じて、対角性に向かう収束を示す図である。
【図３】混合された信号を分離するための三つの異なる方法（Ｊ−Ｄｉ、ＦＦＤｉａｇおよびＬＵＪ１Ｄ）について、サイズＮ＝５０のＫ＝５０の行列のセットを用いた５０回の独立した試行での、多くのスイープに応じて、分離基準（ｓｅｐａｒａｔｉｏｎ  ｃｒｉｔｅｒｉｏｎ）の収束を示す図である。
【図４】Ｋ＝Ｎ（Ｎ＋１）／２＝２１０  四次キュムラントスライスＮ＝２０およびσ＝０．１である場合の、５０回の独立した試行での、多くのスイープに応じて、分離基準の収束を示す図である。
【図５】本発明に係る混合された信号を分離するためのシステムの実施態様を表示したものである。
【発明を実施するための形態】
【００３７】
当業者にとって公知の以下の代数的性質が、以降で使用されている：
−二つの行列ＡおよびＢの積の転置：
（ＡＢ）Ｔ＝ＢＴＡＴ
−二つの行列ＡおよびＢの積のインバース：
（ＡＢ）−１＝Ｂ−１Ａ−１
−ギブンス回転：
Ｇ（θ）−１＝Ｇ（θ）Ｔ＝Ｇ（−θ）
−ハイパボリック回転
Ｈ（φ）−１＝Ｈ（−φ）
−ハイパボリック回転の対称性：
Ｈ（φ）＝Ｈ（φ）Ｔ
【００３８】
本発明に係る方法の一つの実施態様についてここで図１を参考にしながらより詳細に記述する。
【００３９】
混合体を独立した成分に分解することは、混合体の中で相互に隠れている現象の各々の解釈、分析または処理を容易にすることから有利である。
【００４０】
これらの現象は、例えば以下のものであってよい：
−脳波信号（ＥＥＧ）に寄与するさまざまな（脳および筋肉）電気現象、
−狭周波数帯域のさまざまな無線送信機、
−溶液中に存在する、異なる蛍光色素からの発光（ｌｉｇｈｔ  ｅｍｉｓｓｉｏｎ）、
−同じ分子の異なる混合体の質量スペクトル。
【００４１】
異なる寄与を分離することにより、このとき次のことが可能である。
−筋肉寄生（ｍｕｓｃｕｌａｒ  ｐａｒａｓｉｔｅ）を免れたＥＥＧ信号を研究すること（明滅（ｂｌｉｎｋ）、ＥＣＧなど）；
−各高周波ソースからの信号を別々に復調すること；
−各蛍光色素の濃度を（例えばサイトメトリによって）推定すること、または、
−質量分析法を通して研究された媒質内に存在する、異なる分子を識別すること。
【００４２】
ソースの信号は、時間（ＥＥＧ、高周波数）、波長（蛍光性）、質量（質量スペクトル）などであってもよい可変数量に左右される。
【００４３】
以下では、ソース由来の信号が時間に左右されるという点を考慮に入れながら、本発明に係る方法を提示するが、これによってこの方法の全体的内容が限定されることはない。
【００４４】
ソースの分離の問題は、以下の要領で数学的にモデル化される。
Ｍ個のソースの各々は、他のソースが発出する信号とは独立した形で信号ｓｍ（ｔ）を発出する。
これらのＭ個の寄与の（所与の瞬間ｔにおける）混合体は、Ｎ個の全く異なる混合体を提供するＮ個のセンサーのネットワーク（ここでＮ＞ＭまたはＮ＝Ｍ）により観察される。
したがって、ｘ（ｔ）と記される瞬間ｔにおいて観察されたデータは、多次元である。
【００４５】
所与の瞬間におけるこれらのＮ回の観察のベクトルは、行列Ａ（Ｎ行、Ｍ列）つまりいわゆる「混合行列」とソース信号のベクトルｓ（ｔ）との積として記してよい。すなわち
【００４６】
【数４】
2011509037.tif 000006
【００４７】
先の等式を再公式化することで測定値ｘ（ｔ）に対する異なるソースの寄与ｓｍ（ｔ）を示してもよい。
【００４８】
【数５】
2011509037.tif 000007
【００４９】
本発明に係る方法によると、混合行列Ａに対するいかなる先験的情報もなくこの混合行列を推定することができ（この行列は単に非特異行列であることだけが仮定されている）、このときソース信号が推定される。
【００５０】
ソース信号は独立していると仮定され、その非ガウス性（ｎｏｎ−Ｇａｕｓｓｉａｎ  ｃｈａｒａｃｔｅｒ）、そのスペクトル差、またはその非定常性などのその統計的特性が利用される。
【００５１】
第一のステップ１００においては、Ｋ個の推定された行列Ｒ：ｋの族が、センサーネットワークにより記録されたデータから（１とＫの間に含まれるｋについて）計算される。
これらの推定された行列Ｒ：ｋは、以下の構造を有する行列Ｒｋと（推定エラー内まで）等しい：
【００５２】
【数６】
2011509037.tif 000008
【００５３】
なお式中、Ｄｋは対角行列であり、ＡＴは転置された混合行列Ａである。当然のことながら、当業者であれば、行列Ｄｋがブロック状に対角であってよいということを理解すると思われる。
【００５４】
等式の第３項は、ｒａｎｋが１のＭ個の行列の加重和（ａｍａＴｍ）としての、Ｒｋ個の行列を示す第２項の再公式化にすぎない。項Ｄｋ（ｍ，ｍ）ａｍａＴｍは、行列Ｒｋに対するソースｍの寄与である。
【００５５】
これらの状況に応じて、これらの推定行列Ｒ：ｋは、３次または４次のキュムラントの行列（ＪＡＤＥ方法［１］、［１ｂ］）、異なる周波数帯域内の共分散行列［３］、異なる時間推移の相互相関行列（ＳＯＢＩ方法［２］）または異なる時間的間隔にわたり推定された共分散行列［３］であってよい。当然のことながら、これらの推定行列は、当業者にとって公知のその他のタイプのものであってよい。
【００５６】
第二のステップ２００においては、推定行列の結合対角化が、以下でより詳細に記述されているように実施される。
【００５７】
本発明に係る方法は、Ａ行列をＲ：１，．．．，Ｒ：ｋ，．．．Ｒ：Ｋ行列から計算することを可能にする。
【００５８】
行列Ｂは、すべての行列ＢＲｋ：ＢＴが可能なかぎり「対角」となるように求められる：
ＢＲ：ｋＢＴ＝Ｄｋ
Ｂは、以下の形の基本Ｎ×Ｎ行列の積として構築される：
【００５９】
【数７】
2011509037.tif 000009
【００６０】
ここで、ｉ＜ｊは、１とＮの間の任意の二つのインデックスであり、ＩＰは、ｐ×ｐサイズの恒等行列を表わし、下記〔数８〕は、下記〔数１１〕の等式を満たすような、下記〔数９〕タイプのギブンス行列（なお式中θは、ギブンス回転の角度である）と〔数１０〕タイプのハイパボリック回転行列（なお式中φはハイパボリック回転角度である）との積である。
【００６１】
【数８】
2011509037.tif 000010
【００６２】
【数９】
2011509037.tif 000011
【００６３】
【数１０】
2011509037.tif 000012
【００６４】
【数１１】
2011509037.tif 000013
【００６５】
したがって、先行技術の対角化方法との差異は、行列Ｂ（ひいては混合行列Ａ）が、回転行列と一つまたは複数のギブンス回転行列および一つまたは複数のハイパボリック回転行列との積から推定されるという点にある。
【００６６】
「ハイパボリック回転行列」とは、φのハイパボリック関数に対応する要素を含む行列のことである。
【００６７】
より詳細には、推定行列Ｒ：ｋの族には以下のものが乗算される：
−各推定行列Ｒ：ｋの右側に、ギブンス回転行列とハイパボリック回転行列との積。当業者にとっては、該方法の等価のインプリメンテーション（ｅｑｕｉｖａｌｅｎｔ  ｉｍｐｌｅｍｅｎｔａｔｉｏｎ）が回転の反転（ｉｎｖｅｒｓｉｏｎ）からもたらされるであろうことは明白である（最初にギブンス回転で次にハイパボリック回転ではなく、ハイパボリック回転で次にギブンス回転）。
−そして、各推定行列Ｒ：ｋの左側に、この積の転置行列：すなわち
【００６８】
【数１２】
2011509037.tif 000014
【００６９】
となる。ここで
【００７０】
【数１３】
2011509037.tif 000015
【００７１】
である。
【００７２】
この積の結果としてもたらされる行列、下記〔数１４〕の非対角項を最小化する角度θ０およびφ０が次に求められる。
【００７３】
【数１４】
2011509037.tif 000016
【００７４】
この演算は、ギブンスおよびハイパボリック回転行列の要素を、結果として得られる行列の以下の要素へと移動させることによって、結果として得られた行列について反復される：
【００７５】
【数１５】
2011509037.tif 000017
【００７６】
ここで
【００７７】
【数１６】
2011509037.tif 000018
【００７８】
であって、インデックスｉおよびｊの全ての対について同様である。
【００７９】
【数１７】
2011509037.tif 000019
【００８０】
これらの演算は、ｃｏｓθおよびｃｏｓｈφがそれぞれに１に向かうような、全てのギブンス回転角度θおよび全てのハイパボリック回転角度φとなるまで、インデックスｉおよびｊの全ての対について複数回反復される。
【００８１】
当然のことながら、当業者であれば、（該方法の大域収束性を意味する）他の停止基準が可能であるということを知っている。例えば、対角要素の全体的増加または非対角要素の全体的減少を監視して、反復停止を決定してもよい。
【００８２】
行列Ｍ（ｉ，ｊ，θ，φ）Ｔ×Ｒ：ｋ×Ｍ（ｉ，ｊ，θ，φ）の非対角項を全体的に最小化する角度θおよびφは、以下の要領で得ることができる：
【００８３】
−行列Ｃを以下のように構築する。
【００８４】
【数１８】
2011509037.tif 000020
【００８５】
−ＣＴ×Ｃおよび対角行列（−１，１，１）の一般固有分解（ｇｅｎｅｒａｌｉｚｅｄ  ｅｉｇｅｎｄｅｃｏｍｐｏｓｉｔｉｏｎ）を計算する。固有ベクトルはｚ＞０および−ｘ２＋ｙ２＋ｚ２＝１となるような正の最小の一般固有値のν＝（ｘ，ｙ．ｚ）Ｔとして記される。
【００８６】
【数１９】
2011509037.tif 000021
【００８７】
を計算する。
【００８８】
【数２０】
2011509037.tif 000022
【００８９】
を計算する。
【００９０】
問題は以下の等式を解くことのみにあることから、その他の等価の式も存在することを当業者であれば認識する：
【００９１】
【数２１】
2011509037.tif 000023
【００９２】
上述の式は、該方法の収束の終りでｘが０に向かう場合に優れた数値的安定性を有するという利点をもつ。
【００９３】
コンピュータ計算の有限精度に適応させるためには、その他のマイナーな修正が有用であるかもしれない。例えば、正の最小の固有値は、理論的には実際上ごくわずかに負であるかもしれない。それでも、それを保持することが必要である。これらの点は、当業者にとって明白である。
【００９４】
行列Ｃからの角度θおよびφの計算のより完全な説明を以下で記述する。
【００９５】
論理的進行は以下の通りである。
【００９６】
ｋ＝１，．．．，Ｋについて、結果として得られる行列の積（Ｒ’ｋと記す）はＲ’ｋ＝Ｍ（ｉ，ｊ，θ，φ）Ｔ×Ｒ：ｋ×Ｍ（ｉ，ｊ，θ，φ）と計算される：
【００９７】
【数２２】
2011509037.tif 000024
【００９８】
θのコサインとθのサインの積ならびにφのハイパボリックコサイン（ｈｙｐｅｒｂｏｌｉｃ  ｃｏｓｉｎｅｓ）とφのハイパボリックサイン（ｈｙｐｅｒｂｏｌｉｃ  ｓｉｎｅｓ）の積が出現し、倍角（２θ、２φ）のコサインおよびサインへと変換され、これらは項のグルーピングの時点で、以下の等式を提供する：
【００９９】
【数２３】
2011509037.tif 000025
【０１００】
ここで、
【０１０１】
【数２４】
2011509037.tif 000026
【０１０２】
このとき、結果として得られる行列の積Ｍ（ｉ，ｊ，θ，φ）Ｔ×Ｒ：ｋ×Ｍ（ｉ，ｊ，θ，φ）の非対角要素（ｉ，ｊ）の二乗の和は二次形式ν（θ，φ）ＴＣＴＣν（θ，φ）に等しい。
【０１０３】
ベクトルν（θ，φ）がν（θ，φ）ＴＪν（θ，φ）＝１（なおＪ＝ｄｉａｇ（−１，１，１））により定義される双曲面（ｈｙｐｅｒｂｏｌｏｉｄ）を記述する、という点に留意されたい。
【０１０４】
制約ν（θ，φ）ＴＪν（θ，φ〉＝１の下で二次形式ν（θ，φ）ＴＣＴＣν（θ，φ）を最大化するベクトルν（θ，φ）が求められる。このベクトルが必然的に行列対（ＣＴＣ，Ｊ）の三つの一般固有ベクトルのうちの一つであるということは周知である。
【０１０５】
この一般固有分解の研究は、これが（ＣＴＣ，Ｊ）の正の最小の一般固有値に対応する一般固有値ν＝（ｘｙｚ）Ｔであることを示している：
【０１０６】
【数２５】
2011509037.tif 000027
【０１０７】
最終的に、積Ｍ（ｉ，ｊ，θ，φ）Ｔ×Ｒ：ｋ×Ｍ（ｉ，ｊ，θ，φ）を計算するために、θおよびφの最適値は不要である。その（従来のおよびハイパボリックの）コサインおよびサインの値のみが必要とされ、これは例えば以下の式から得られる：
【０１０８】
【数２６】
2011509037.tif 000028
【０１０９】
したがって該方法は、１とＮの間のｉ＜ｊについて、収束まで数回、
Ｍ（ｉ，ｊ，θ，φ）Ｔ×Ｒｋ：×Ｍ（ｉ，ｊ，θ，φ）、
の積の反復を含む。
【０１１０】
混合行列Ａはこのとき
【０１１１】
【数２７】
2011509037.tif 000029
【０１１２】
のインバース（ｉｎｖｅｒｓｅ）に等しい。混合行列Ａがひとたび決定されると、ソース信号を分離することが可能になる。
【０１１３】
Ｍａｔｌａｂ言語での例示的なインプランテーション（ｉｍｐｌａｎｔａｔｉｏｎ）を以下に示す：
ｆｕｎｃｔｉｏｎ［Ａ，ｎｕｍｂｅｒ＿ｏｆ＿ｌｏｏｐｓ，ｍａｔｒｉｃｅｓ］＝Ｊ＿Ｄｉ［ｍａｔｒｉｃｅｓ，ｔｈｒｅｓｈｏｌｄ］
％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％
％
％  Ｊ＿Ｄｉ．ｍ  ｗｒｉｔｔｅｎ  ｔｈｅ  １９ｔｈ  ｏｆ  ｄｅｃｅｍｂｅｒ  ２００７  ｂｙ  Ａｎｔｏｉｎｅ  Ｓｏｕｌｏｕｍｉａｃ．
％
％  Ｔｈｉｓ  ｓｏｆｔｗａｒｅ  ｃｏｍｐｕｔｅｓ  ｔｈｅ  Ｊｏｉｎｔ  Ｄｉａｇｏｎａｌｉｓａｔｉｏｎ  ｏｆ  ａ  ｓｅｔ  ｏｆ  ｍａｔｒｉｃｅｓ．
％
％  ｉｎｐｕｔ：
％  ｍａｔｒｉｃｅｓ：  ａｒｒａｙ  ｏｆ  ｍａｔｒｉｃｅｓ  ｔｏ  ｂｅ  ｊｏｉｎｔｌｙ  ｄｉａｇｏｎａｌｉｚｅｄ
％  ｔｈｒｅｓｈｏｌｄ：  ｔｈｒｅｓｈｏｌｄ  ｆｏｒ  ｓｔｏｐ  ｃｒｉｔｅｒｉｏｎ  ｏｎ  ｍａｘｉｍａｌ  ｓｉｎｅｓ
％
％  ｏｕｔｐｕｔｓ：
％  Ａ：  ｅｓｔｉｍａｔｅｄ  ｍｉｘｔｕｒｅ  ｍａｔｒｉｘ
％  ｎｕｍｂｅｒ＿ｏｆ＿ｌｏｏｐｓ：ｎｕｍｂｅｒ  ｏｆ  ｌｏｏｐｓ  ｎｅｃｅｓｓａｒｙ  ｆｏｒ  ｃｏｎｖｅｒｇｅｎｃｅ
％  ｍａｔｒｉｃｅｓ：ａｒｒａｙ  ｏｆ  ｍａｔｒｉｃｅｓ  ａｆｔｅｒ  ｊｏｉｎｔ  ｄｉａｇｏｎａｌｉｚａｔｉｏｎ
％
％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％
％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％％
％  ｉｎｉｔｉａｌｉｚａｔｉｏｎ
ｍａｔ ｓｉｚｅ  ＝  ｓｉｚｅ（ｍａｔｒｉｃｅｓ，１）；
ｎｂ ｍａｔ  ＝  ｓｉｚｅ（ｍａｔｒｉｃｅｓ，３）；
Ａ  ＝  ｅｙｅ（ｍａｔ＿ｓｉｚｅ）；
ｇｏ＿ｏｎ  ＝  １；
ｎｕｍｂｅｒ＿ｏｆ＿ｌｏｏｐｓ  ＝  −１；
％  ｊｏｉｎｔ  ｄｙａｄｉｚａｔｉｏｎ  ｌｏｏｐｓ
ｗｈｉｌｅ  ｇｏ＿ｏｎ
ｆｏｒ  ｎ＝  ２：ｍａｔ＿ｓｉｚｅ
ｆｏｒ  ｍ＝  １：ｎ−１
％  ｉｎｉｔｉａｌｉｚａｔｉｏｎ
ｓ ｍａｘ  ＝  ０．０；
ｓｈ ｍａｘ  ＝  ０．０；
％  Ｇｉｖｅｎｓ  ａｎｄ  ｈｙｐｅｒｂｏｌｉｃ  ａｎｇｌｅｓ  ｃａｌｃｕｌａｔｉｏｎ
Ｃ＝［（ｓｑｕｅｅｚｅ（ｍａｔｒｉｃｅｓ（ｍ，ｍ，：））＋ｓｑｕｅｅｚｅ（ｍａｔｒｉｃｅｓ（ｎ，ｎ，：）））／２．．．
（ｓｑｕｅｅｚｅ（ｍａｔｒｉｃｅｓ（ｍ，ｍ，：））−ｓｑｕｅｅｚｅ（ｍａｔｒｉｃｅｓ（ｎ，ｎ，：）））／２．．．
ｓｑｕｅｅｚｅ（ｍａｔｒｉｃｅｓ（ｍ，ｎ，：））］；
［Ｖ，Ｄ］＝ｅｉｇ（Ｃ’＊Ｃ，ｄｉａｇ（［−１  １  １］））；
％  ｎｏｒｍａｌｉｓａｔｉｏｎ（−１  １  １  ）  ｏｆ  ｇｅｎｅｒａｌｉｚｅｄ  ｅｉｇｅｎｖｅｃｔｏｒｓ
Ｖ＝Ｖ＊ｄｉａｇ（１．／ｓｑｒｔ（ａｂｓ（ｄｉａｇ（Ｖ’＊ｄｉａｇ（［−１  １１］）＊Ｖ））））；
％  ｃｏｍｐｕｔａｔｉｏｎ  ｏｆ  ｔｈｅ  ｍｅｄｉｕｍ  ｇｅｎｅｒａｌｉｚｅｄ  ｅｉｇｅｎｖａｌｕｅ
［ｓｏｒｔＤ，ｉｎｄｅｘ］＝  ｓｏｒｔ（ｄｉａｇ（Ｄ））；
ｉｎｄ＿ｍｉｎ  ＝  ｉｎｄｅｘ（２）；
ｖ＝Ｖ（：，ｉｎｄ＿ｍｉｎ）＊ｓｉｇｎ（Ｖ（３，ｉｎｄ＿ｍｉｎ））；
％  ｃｏｍｐｕｔａｔｉｏｎ  ｏｆ  ｔｒｉｇｏｎｏｍｅｔｒｉｃ  ａｎｄ  ｈｙｐｅｒｂｏｌｉｃ  ｃｏｓｉｎｅｓ  ａｎｄ  ｓｉｎｅｓ
ｃｈ  ＝  ｓｑｒｔ（（１＋ｓｑｒｔ（１＋ｖ（１）∧２））／２）；
ｓｈ  ＝  ｖ（１）／（２＊ｃｈ）；
ｃ  ＝  ｓｑｒｔ（（１＋ｖ（３）／ｓｑｒｔ（１＋ｖ（１）∧２））／２）；
ｓ  ＝  −ｖ（２）／（２＊ｃ＊ｓｑｒｔ（１＋ｖ（１）∧２））；
％  ｍａｘｉｍａｌ  ｓｉｎｅ  ａｎｄ  ｓｉｎｈ  ｃｏｍｐｕｔａｔｉｏｎ  ｆｏｒ  ｓｔｏｐ  ｔｅｓｔ
ｓ ｍａｘ  ＝  ｍａｘ（ｓ＿ｍａｘ，ａｂｓ（ｓ））；
ｓｈ＿ｍａｘ  ＝  ｍａｘ（ｓｈ＿ｍａｘ，ａｂｓ（ｓｈ））；
％  ｐｒｏｄｕｃｔ  ｏｆ  Ｇｉｖｅｎｓ  ａｎｄ  ｈｙｐｅｒｂｏｌｉｃ  ｒｏｔａｔｉｏｎｓ
ｒｍ１１  ＝  ｃ＊ｃｈ  −  ｓ＊ｓｈ；
ｒｍ１２  ＝  ｃ＊ｓｈ  −  ｓ＊ｃｈ；
ｒｍ２１  ＝  ｃ＊ｓｈ  ＋  ｓ＊ｃｈ；
ｒｍ２２  ＝  ｃ＊ｃｈ  ＋  ｓ＊ｓｈ；
％  ｍａｔｒｉｃｅｓ  Ｇｉｖｅｎｓ  ｕｐｄａｔｅ
ｈ＿ｓｌｉｃｅ１  ＝  ｓｑｕｅｅｚｅ（ｍａｔｒｉｃｅｓ（ｍ，：，：））；
ｈ＿ｓｌｉｃｅ２  ＝  ｓｑｕｅｅｚｅ（ｍａｔｒｉｃｅｓ（ｎ，：，：））；
ｍａｔｒｉｃｅｓ（ｍ，：，：）＝  ｒｍ１１＊ｈ＿ｓｌｉｃｅ１  ＋  ｒｍ２１＊ｈ＿ｓｌｉｃｅ２；
ｍａｔｒｉｃｅｓ（ｎ，：，：）＝  ｒｍ１２＊ｈ＿ｓｌｉｃｅ１  ＋  ｒｍ２２＊ｈ＿ｓｌｉｃｅ２；
ｖ＿ｓｌｉｃｅ１  ＝  ｓｑｕｅｅｚｅ（ｍａｔｒｉｃｅｓ（：，ｍ，：））；
ｖ＿ｓｌｉｃｅ２  ＝  ｓｑｕｅｅｚｅ（ｍａｔｒｉｃｅｓ（：，ｎ，：））；
ｍａｔｒｉｃｅｓ（：，ｍ，：）＝  ｒｍ１１＊ｖ＿ｓｌｉｃｅ１  ＋  ｒｍ２１＊ｖ＿ｓｌｉｃｅ２；
ｍａｔｒｉｃｅｓ（：，ｎ，：）＝  ｒｍ１２＊ｖ＿ｓｌｉｃｅ１  ＋  ｒｍ２２＊ｖ＿ｓｌｉｃｅ２；
％  ｄｙａｄｓ  Ｇｉｖｅｎｓ  ｕｐｄａｔｅ
ｃｏｌ１  ＝  Ａ（：，ｍ）；
ｃｏｌ２  ＝  Ａ（：，ｎ）；
Ａ（：，ｍ）＝  ｒｍ２２＊ｃｏｌ１  −  ｒｍ１２＊ｃｏｌ２；
Ａ（：，ｎ）＝  −ｒｍ２１＊ｃｏｌ１  ＋  ｒｍ１１＊ｃｏｌ２；
ｅｎｄ
ｅｎｄ
ｎｕｍｂｅｒ＿ｏｆ＿ｌｏｏｐｓ＝ｎｕｍｂｅｒ＿ｏｆ＿ｌｏｏｐｓ＋１；
％  ｒｅｎｏｒｍａｌｉｚａｔｉｏｎ  ｏｆ  ｔｈｅ  Ａ  ｃｏｌｕｍｎｓ（ｅｑｕａｌ  ｎｏｒｍ）
ｃｏｌ＿ｎｏｒｍ  ＝  ｓｑｒｔ（ｓｕｍ（Ａ．∧２，１））；
ｄ  ＝（ｐｒｏｄ（ｃｏｌ＿ｎｏｒｍ）．∧（１／ｍａｔ＿ｓｉｚｅ））．／ｃｏｌ＿ｎｏｒｍ；
Ａ  ＝  Ａ．＊ｒｅｐｍａｔ（ｄ，ｍａｔ＿ｓｉｚｅ，１）；
ｍａｔｒｉｃｅｓ  ＝  ｍａｔｒｉｃｅｓ．＊ｒｅｐｍａｔ（（１．／ｄ）’＊（１．／ｄ），［１  １  ｎｂ＿ｍａｔ］）；
％  ｓｔｏｐ  ｔｅｓｔ
ｉｆ  ｓ＿ｍａｘ＜ｔｈｒｅｓｈｏｌｄ  ＆  ｓｈ＿ｍａｘ＜ｔｈｒｅｓｈｏｌｄ
ｇｏ＿ｏｎ＝０；
ｄｉｓｐ（［’ｓｔｏｐ  ｗｉｔｈ  ｍａｘ  ｓｉｎｕｓ  ａｎｄ  ｓｉｎｈ＝’  ｎｕｍ２ｓｔｒ（［ｓ＿ｍａｘ  ｓｈ＿ｍａｘ］）］）
ｅｎｄ
ｉｆ  ｎｕｍｂｅｒ＿ｏｆ＿ｌｏｏｐｓ＞１００
ｇｏ＿ｏｎ  ＝  ０；
ｗａｒｎｉｎｇ（’ｍａｘｉｍｕｍ  ｎｕｍｂｅｒ  ｏｆ  ｌｏｏｐｓ  ｈａｓ  ｂｅｅｎ  ｒｅａｃｈｅｄ’）
ｅｎｄ
ｅｎｄ
【０１１４】
計算の複雑性は以下のような近似値を有する：
スイープ数×２ＫＮ３
なお式中、
−スイープ数は、インデックス対（ｉ，ｊ）全体について反復が行われる回数であり、
−ＫはＲ：ｋ行列の数であり、
−ＮはＲ：ｋ行列のサイズである。
【０１１５】
本発明および既存のアルゴリズム例えばＦＦＤｉａｇ［４］、ＤＮＪＤ［６］、ＬＵＪ１Ｄ［７］などは、２ＫＮ３に等しくスイープすることにより等価の計算コストを有する。本発明の利点は、スイープ数にある。より具体的には、当該方法によると、先行技術の公知の方法と比べてスイープ数を削減してよい。
【０１１６】
一例として、本発明に係る方法（以下「Ｊ−Ｄｉ」と呼ぶ）とＦＦＤｉａｇ（Ｚｉｅｈｅ［４］）およびＬＵＪ１Ｄ（Ａｆｓａｒｉ［７］）方法の比較は以下の通りに提示される。すなわちこの比較は、先行技術の方法が、本発明の方法のものの少なくとも二倍のスイープ数を有するという事実を示す。
【０１１７】
その上、先行技術の一部の方法は、特定の行列セットに特化されている（Ｐｈａｍ［３］：正値定符号行列、Ｖｏｌｌｇｒａｆ［５］：少なくとも一つの正値定符号行列を含むセット）。本発明に係る方法はあらゆるタイプの行列セットの処理を可能にする。
【０１１８】
図５は本発明に係る混合された信号を分離するためのシステムの実施態様を示す。このシステムは、例えばプロセッサーのような処理手段１０（一つまたは複数のギブンス回転行列および一つまたは複数のハイパボリック回転行列に対応する回転行列の積から混合行列Ａを推定する結合対角化用）、例えば表示スクリーンのような表示手段２０、そして、例えばキーボードやマウスのような入力手段３０を備える。システムは、処理すべき混合された信号（ｘ１（ｔ），．．．，ｘＮ（ｔ））を受信するため通信手段に接続されている。
【０１１９】
読者は、本明細書の教示から実質的に逸脱することなく上述に記載の本発明に対し数多くの修正を加えることができると認識するものである。
【０１２０】
例えば、ギブンス回転（Ｇ（θ））とハイパボリック回転（Ｈ（φ））行列の積は、ハイパボリック回転行列にギブンス回転行列を乗じたもの、または、ギブンス回転行列にハイパボリック回転行列を乗じたもののいずれかに対応することができる。
【０１２１】
結果として、このタイプのすべての修正は、添付の請求の範囲に包含されるものと意図される。
【０１２２】
「ＡＮＮＥＸ」
添付の図面により、上述の方法をよりよく理解することができる。先に用いた一部の表記法は、以下では異なる可能性がある。いずれにせよ、表記法の変更は読者の理解を容易にする目的で特定されるものである。
【０１２３】
Ｉ．要約
独立成分分析およびブラインドソース分離（ｂｌｉｎｄ  ｓｏｕｒｃｅ  ｓｅｐａｒａｔｉｏｎ）のアプリケーション用として、行列セットの非直交結合対角化を計算するための新しいアルゴリズムが提案されている。
【０１２４】
このアルゴリズムは、直交結合対角化用の、固有行列の同時近似対角化（Ｊｏｉｎｔ  Ａｐｐｒｏｘｉｍａｔｅ  ｄｉａｇｏｎａｌｉｚａｔｉｏｎ  ｏｆ  ｅｉｇｅｎｍａｔｒｉｃｅｓ）（ＪＡＤＥ）方法において最初に提案されたヤコビのようなアルゴリズムの拡張である。改善は、主として、直交混合行列の代りに等ノルム（ｅｑｕａｌ  ｎｏｒｍ）の列および行列式１の混合行列を計算することからなる。この標的行列は、ギブンス回転のみならずハイパボリック回転および対角行列の連続的乗算によって反復的に構築される。合成データについて評価されたアルゴリズム性能は、収束速度および複雑性の視点から既存の方法と比較して有利である。
【０１２５】
ＩＩ．序
結合対角化アルゴリズムは、ブラインドソース分離（ＢＳＳ）および独立成分分析（ＩＣＡ）のアプリケーションにおいて使用される。
【０１２６】
これらのアプリケーションでは、瞬時に混合された独立した信号（ソース）を分離すること、すなわち、混合体のみを観察する行列Ａによる一時的コンボリューション（ｔｅｍｐｏｒａｌ  ｃｏｎｖｏｌｕｔｉｏｎ）無しで、かつ、ソースの統計的独立性以外の先行情報無しで、独立した信号（ソース）を分離することが目標とされる。
【０１２７】
［１］、［１ｂ］、［２］、［３］、［４］、［５］、［１１］などの数多くのＢＳＳ方法は、Ｎ×ＮサイズのＫ個の対称行列Ｒ１，．．．，Ｒｋ，．．．，ＲＫのセットＲの結合対角化ステップを含む。手元のアプリケーションおよびソース信号の対応する予想特性によると、考慮対象の行列Ｒｋは、例えば異なる時間窓で、三次または四次（［１］，［１ｂ］）以上のキュムラントスライスに基づいて推定される共分散行列、または時間推移を伴う相互相関行列［２］などである。これらの行列はすべて、共通の構造を共有する：
【０１２８】
【数２８】
2011509037.tif 000030
【０１２９】
なお式中、Ａは混合行列であり、Ｄｋは対角行列である。これらの条件下で、
【０１３０】
【数２９】
2011509037.tif 000031
【０１３１】
［数２９］により定義される分離行列は、各Ｒｋを対角化する特性、すなわち１≦ｋ≦ＫについてＢＲｋＢＴ＝Ｄｋにより特徴づけられる。弱い仮説の下で、置換および不定スケーリング（ｓｃａｌｉｎｇ  ｉｎｄｅｔｅｒｍｉｎａｔｉｏｎ）までは、Ｒｋ行列の結合対角化手順により分離行列また同等に混合行列を推定するのにこの特性で充分である。実際には、Ｒｋ行列は、一部の推定エラーによって破損され、Ａ（またはＢ）はセットＲの近似結合対角化を介して効率良く推定される。
【０１３２】
その上、例えばＪＡＤＥ［１］、［１ｂ］またはＳＯＢＩ［２］のような最初の結合対角化技術には、Ｒ内の一つの行列の、一般的には（ノイズ無しの）混合されたソースの共分散行列の、正定性が必要とされる。これに関連して、このとき、混合されたソースを予め白色化することが可能であり、こうして直交（または複雑なケースではユニタリ）混合行列を計算する上での結合対角化の問題が削減される。残念ながら、この白色化ステップは、混合されたソースの共分散行列内の推定エラーに起因して決して正確ではない。この予備的エラーは、後続する直交結合対角化（ＯＪＤ）により補正できない。結果として、数多くの製作者が、予備白色化を回避するための非直交結合対角化（ＮＯＪＤ）方法を開発した。これらの最適化方法は、正値定符号行列のセット［３］、またはＱＤｉａｇのような一つの正値定符号行列を含むセット［１２］、［５］または最終的に任意の行列セット［１１］、［４］、［７］に対応している。
【０１３３】
ＪＡＤＦなどのヤコビ技術の優れた収束特性を享受する代数アプローチに基づいたＢＳＳおよびＩＣＡアプリケーション用の任意の行列セットのための新たな非直交結合対角化アルゴリズムが提案されている。この結果を得るために、混合行列Ａは同じく基本行列の乗算によって構築される。しかしながら、Ａは直交でないことから、下位群として直交行列群を含む新しい乗法群、およびギブンス回転を一般化する基本行列の新しいセットを定義することが必要である。これは、ギブンス回転、ハイパボリック回転および対角行列の連続的乗算により等ノルムの列を有し１に等しい行列式の混合行列を反復的に構築することによって達成できる。
【０１３４】
以下の第三節では、混合行列の列ノルムおよび行列式の制約およびそれがＮＯＪＤ問題の不確定にもたらす帰結について吟味する。
【０１３５】
第四節では、行列式１の行列をいかにしてギブンス回転、ハイパボリック回転および対角行列の積として生成することができるかを示す。
【０１３６】
Ｊｏｉｎｔ−Ｄｉａｇｏｎａｌｉｚａｔｉｏｎ（結合対角化）を略してＪ−Ｄｉと呼ばれる新しい非直交近似結合対角化アルゴリズムは第五節で記述されており、その性能は、第六節内で合成データについて評価される。特に、例えば２００×２００サイズの２００個の行列などの大きい行列の大きいセットについてのシミュレーションが提供されている。これは、電極数が１２８個または２５６個に達する場合の脳波図などの利用分野にとって有利である。
【０１３７】
最後に、［６］においてＷａｎｇ、ＬｉｕおよびＺｈａｎｇによって提案された最近のアルゴリズムとの短い複雑性比較が付加されている。
【０１３８】
ＩＩＩ．１に等しい行列式および等ノルム列の行列
ＩＩＩ．Ａ．特性
【０１３９】
モデル（１）内の任意の正則正方混合行列Ａが１に等しい行列式および等ノルムの列を有するものと仮定される。実際のところ、（１）の中に正則対角行列ＤＡを挿入することは常に可能である。
【０１４０】
【数３０】
2011509037.tif 000032
【０１４１】
【数３１】
2011509037.tif 000033
【０１４２】
行列ＤＡ−１ＤｋＤＡ−１は明らかに対角であり、ＤＡは、ＡＤＡが所要の行列式および列ノルム特性を確認するような形で調整可能である。
【０１４３】
【数３２】
2011509037.tif 000034
【０１４４】
【数３３】
2011509037.tif 000035
【０１４５】
を考慮するだけで充分である。
【０１４６】
ここで、ａｎはＡのｎ番目の列であり、｜｜ａｎ｜｜はそのユークリッド・ノルムである。混合行列に対するこれらの制約は、数値的精度にとって有用な一種の正規化を提供する。特に、ＡＤＡの条件数はＡの条件数よりも小さいことが多いが、これは通則ではない。
【０１４７】
ＩＩＩ．Ｂ．混合行列の一意性条件および不確定
［９］、［１０］に記述されているように、Ｄｋの対角線（ＤはＫ×Ｎであり、その（ｋ，ｎ）エントリはＤｋのｎ番目の対角エントリである）の垂直連結（ｖｅｒｔｉｃａｌ  ｃｏｎｃａｔｅｎａｔｉｏｎ）の結果もたらされる行列Ｄ内に共線列（ｃｏｌｉｎｅａｒ  ｃｏｌｕｍｎｓ）が全く存在しない場合、そしてその場合にのみＮＯＪＤ問題に対する一意解が存在する。Ρｎｍが、Ｄのｎ番目およびｍ番目の列の角度のコサインの絶対値であり、
【０１４８】
【数３４】
2011509037.tif 000036
【０１４９】
である場合には、一意性条件を単に
【０１５０】
【数３５】
2011509037.tif 000037
【０１５１】
として再公式化できる。一意性のモジュラスと呼ばれるこのパラメータρは、ρが１に近いほどＮＯＪＤはより難しいということを明白な規則として、ＮＯＪＤ問題の難度を示すものでもある。
【０１５２】
ＮＯＪＤ問題の「一意解」が、混合行列の列の古典的置換およびスケール不確定までのみ「一意的である」ということが理解される。（本発明の）単位行列式および等しい列ノルム特性は、これらの不確定性をわずかに変える。実際のところ、置換不確定は存在し続けるがスケール不確定が符号不確定（ｓｉｇｎ  ｉｎｄｅｔｅｒｍｉｎａｔｉｏｎ）まで削減されることは容易に示すことができる。より厳密に言うと、偶数の混合行列の列に−１を乗じることができる。
【０１５３】
ＩＶ．特別な線形群全体にわたる基本行列の新たな部類
この節の最終目標は、反復した乗算によりセットＡ（Ｒ，Ｎ）および群ＳＬ（Ｒ，Ｎ）を生成する基本行列セットを提案することにある。
【０１５４】
ＩＶ．Ａ．Ａ（Ｒ，２）内の一意的分解
Ａ（Ｒ，２）（ｄｅｔ（Ａ）＝１および｜｜ａ１｜｜＝｜｜ａ２｜｜）内の行列Ａおよびその極分解
【０１５５】
【数３６】
2011509037.tif 000038
【０１５６】
を考慮し、式中Ｇは直交であり、Ｈは対称的正である。ＡおよびＨの行列式は両方共正であることから、Ｇの行列式はプラス１に等しく、Ｇは必然的に以下の〔数３７〕でギブンス回転Ｇ（θ）である。     
【０１５７】
【数３７】
2011509037.tif 000039
【０１５８】
ここで、同じく行列式がプラス１に等しいＨが実際にはハイパボリック回転であることを示す。
【０１５９】
Ｈのエントリをｅ、ｆおよびｇと表示する。
【０１６０】
【数３８】
2011509037.tif 000040
【０１６１】
Ｈ＝Ｇ（−θ）Ａは、Ａと同様、等ノルムの列を有する。したがって、ｅ２＋ｆ２＝ｆ２＋ｇ２であり、Ｈは構成上正であり、したがってｅ＞０、ｇ＞０そしてｅ＝ｇである。Ｈの行列式は１に等しく、このときｅ２−ｆ２＝１であり、ｅ＝ｃｏｓｈ（φ）およびｆ＝ｓｉｎｈ（φ）であるような一意的角度φが存在する。以下の〔数３９〕で最終的にＨ＝Ｈ（φ）である。
【０１６２】
【数３９】
2011509037.tif 000041
【０１６３】
結論として、Ａ（Ｒ，２）の任意の行列をギブンスおよびハイパボリック回転の一意的積に分解することができる。
【０１６４】
ＩＶ．Ｂ．  ＳＬ（Ｒ，２）中の複数解
ＳＬ（Ｒ，２）の任意の行列Ｆは、Ａ（Ｒ、２）内の行列Ａに２×２直交行列Ｄ（λ）を乗じた積へと、一意的に分解することができる。ここで、
【０１６５】
【数４０】
2011509037.tif 000042
【０１６６】
である。ｆ１、ｆ２はＦの列であり、｜｜ｆ１｜｜、｜｜ｆ２｜｜はそれらのノルムである。上の節で示した通り、ギブンスおよびハイパボリック回転の積へと行列Ａを分解することができる。
【０１６７】
最終的に、二つの一意的角度、すなわち−пと＋пの間のθおよびＲ内のφそして、一意の正の実数λが存在し、次のような関係となる：
【０１６８】
【数４１】
2011509037.tif 000043
【０１６９】
さらなる分解を
【０１７０】
【数４２】
2011509037.tif 000044
【０１７１】
の特異値分解（ｓｉｎｇｕｌａｒ  ｖａｌｕｅ  ｄｅｃｏｍｐｏｓｉｔｉｏｎすなわちＳＶＤ）から導出することができる。
【０１７２】
必然的に、ｄｅｔ（Ｓ）＝１およびｄｅｔ（Ｕ）＝ｄｅｔ（Ｖ）＝±１である。ただし、ｄｅｔ（Ｕ）＝ｄｅｔ（Ｖ）＝＋１となるように、ＵとＶの第一列（または行）の符号を交換することがつねに可能であるという点に留意されたい。このことは、ＵおよびＶの両方をギブンス回転すなわちＵ＝Ｇ（θ１）およびＶ＝Ｇ（θ２）として選択できることを証明している。明らかに、Ｓ＝Ｄ（ｓ１）であり、ここでｓ１＞０は第一のＳ特異値である。以下の式が最終的に得られる。
【０１７３】
【数４３】
2011509037.tif 000045
【０１７４】
さらに、任意のハイパボリック回転Ｈ（φ）の固有値分解（ＥＶＤ）から、以下が得られる：
【０１７５】
【数４４】
2011509037.tif 000046
【０１７６】
そして結果として
【０１７７】
【数４５】
2011509037.tif 000047
【０１７８】
これらの分解は、ＳＬ（Ｒ，２）の任意の行列を、ギブンス回転と対角行列、またはギブンス回転とハイパボリック回転の乗算によって生成できるということを証明している。さらに、対角行列を使用して、ＳＬ（Ｒ，２）中の任意の行列の列を正規化し、かくしてＡ（Ｒ，２）の行列を得ることができる。これらの結果はＮ次元まで拡張できる。
【０１７９】
ＩＶ．Ｃ．  ＳＬ（Ｒ，Ｎ）への拡張
以下の表記法が最初に導入される。
【０１８０】
任意のｉ＜ｊ∈｛１，．．．，Ｎ｝Ｇ（θ，ｉ，ｊ）は古典的ギブンス回転
【０１８１】
【数４６】
2011509037.tif 000048
【０１８２】
を表わし、Ｈ（φ，ｉ，ｊ）は、対応するハイパボリック回転
【０１８３】
【数４７】
2011509037.tif 000049
【０１８４】
であり、Ｄ（λ，ｉ，ｊ）は対角行列
【０１８５】
【数４８】
2011509037.tif 000050
を表わす。
【０１８６】
Ｎ＝２のケースときわめて類似して、ＳＬ（Ｒ，Ｎ）中の、行列Ｆの以下の〔数４９〕のＳＶＤが使用される。
【０１８７】
【数４９】
2011509037.tif 000051
【０１８８】
直交行列ＵおよびＶの行列式は、ＵおよびＶの第一列の符号を交換することによりプラス１に等しく選択できる。行列式１（ＳＯ（Ｒ，Ｎ）中）の全ての直交行列をギブンス回転の積の中で分解できるということは周知である。構成上、Ｓは対角正定値であり、その行列式はプラス１に等しい（特異値ｓ１，ｓ２，．．．，ｓＮの積は１に等しい）。同様に以下を示すこともできる。
【０１８９】
【数５０】
2011509037.tif 000052
【０１９０】
さらに、二次元の場合のように、Ｄ（λ，ｉ，ｊ）行列をギブンス回転およびハイパボリック回転の積に変換することができる。最終的に、同じ結論がＳＬ（Ｒ，２）内と同様ＳＬ（Ｒ，Ｎ）内にもあてはまる。すなわち、ＳＬ（Ｒ，Ｎ）内の全ての行列を、ギブンス回転と対角行列またはギブンス回転とハイパボリック回転を乗算することによって生成することができる。対角行列を用いてＳＬ（Ｒ，Ｎ）内の任意の行列の列を正規化し、Ａ（Ｒ，Ｎ）内の行列を得ることができる。
【０１９１】
次節では、上述の三つの基本行列のパラメータをＮＯＪＤコンテクスト内で容易に最適化できるということを見ていく。
【０１９２】
Ｖ．  Ｊ−ＤＩアルゴリズム
本節では、Ｒ１，…Ｒｋ，…ＲＫの正確な行列の推定に対応するセットを結合対角化することにより、基本行列Ｇ（θ，ｉ，ｊ）、Ｈ（φ，ｉ，ｊ）およびＤ（λ，ｉ，ｊ）の積として、１に等しい行列式および等ノルムの列の未知の行列Ａをどのように構築するかについて記述されている。
【０１９３】
より厳密には、行列Ｇ（θ，ｉ，ｊ）およびＨ（φ，ｉ，ｊ）は、非対角エントリ（ｏｆｆ  ｄｉａｇｏｎａｌ  ｅｎｔｒｉｅｓ）を最小化するために使用され、一方でＤ（λ，ｉ，ｊ）行列は列を正規化するために使用される。
【０１９４】
Ｖ．Ａ．  （ｉ，ｊ）非対角エントリを最小化するためのパラメータθおよびφの最適化
ｋ＝１，…，Ｋについて、以下の要領でＲ：ｋを更新することによって得られた行列をＲ’ｋと表わす。
【０１９５】
Ｒ’ｋ＝ＨＴ（φ，ｉ，ｊ）×ＧＴ（θ，ｉ，ｊ）×Ｒ：ｋ×Ｇ（θ，ｉ，ｊ）×Ｈ（φ，ｉ，ｊ）  （１７）
【０１９６】
Ｋ個の行列、Ｒ’ｋのＲ’ｋ［ｉ，ｊ］で表わされる（ｉ，ｊ）−エントリの二乗の和を最小化することが提案される：すなわち
【０１９７】
【数５１】
2011509037.tif 000053
【０１９８】
この数量を最小化することは、例えばＪＡＤＥにおけるようにグローバルコスト関数（ｇｌｏｂａｌ  ｃｏｓｔ  ｆｕｎｃｔｉｏｎ）または非対角エントリの全てを最小化することとは等価ではない。それでも、以下の性能評価は、この反復的ローカル最小化が実際に行列セットＲを結合対角化することを示すものである。一定の操作の後、次のことがわかる：
【０１９９】
【数５２】
2011509037.tif 000054
【０２００】
なお式中
【０２０１】
【数５３】
2011509037.tif 000055
【０２０２】
Ｒ’ｋ中のＧおよびＨの係数の二乗が倍角（三角法角度およびハイパボリック角度）へと変換されるという点に留意されたい。Ｒ’ｋ［ｉ，ｊ］の二乗和は、ν（θ，φ）ＴＣＴＣν（θ，φ）に等しい。したがって（１８）における最小化は次のようになる：
【０２０３】
【数５４】
2011509037.tif 000056
【０２０４】
なお式中、ベクトルν（θ，φ）は、ユークリッド・ノルムではなくハイパボリックノルムで正規化される。実際のところ、−ｖ１（θ，φ）２＋ｖ２（θ，φ）２＋ｖ３（θ，φ）２＝１、すなわち、Ｊ＝ｄｉａｇ（−１，１，１）としてν（θ，φ）ＴＪν（θ，φ）＝１となる。
【０２０５】
ν（θ，φ）ＴＪν（θ，φ）＝１で定義されるハイパボリック全体が−п／２≦θ≦п／２そしてＲ内のφという状態で、ν（θ，φ）により網羅されていることを証明することができる。喚言すると、二次制約（ｑｕａｄｒａｔｉｃ  ｃｏｎｓｔｒａｉｎｔ）ｖＴＪｖ＝１の下で二次数量（ｑｕａｄｒａｔｉｃ  ｑｕａｎｔｉｔｙ）ｖＴＣＴＣｖの最小化が達成される。解ｖは（ＣＴＣ，Ｊ）の一般固有ベクトルであることがわかっている。
【０２０６】
すなわちｖは、下記〔数５５〕となるようなものである。
【０２０７】
【数５５】
2011509037.tif 000057
【０２０８】
この一般固有分解について以下で研究する。
その列が一般固有ベクトルである３×３行列をＶと表示し、（ＣＴＣ，Ｊ）の一般固有値の対角行列をΔと表わすと、ＣＴＣＶ＝ＪＶΔであり、したがってＶＴＣＴＣＶ＝ＶＴＪＶΔである。ＶＴＪＶは対角であり、Ｊのような一つの負の対角エントリと二つの正の対角エントリ（シルベスタの慣性の法則（Ｓｙｌｖｅｓｔｅｒ  ｌａｗ  ｏｆ  ｉｎｅｒｔｉａ）を有することがわかっている。ＣＴＣは正であることから、ΔはＪと同じ慣性（ｉｎｅｒｔｉａ）を有し、次のものが存在する：
・負の一般固有値と負のＪ−ｎｏｒｍ  ｖＴＪｖ＝−１を有する一つの一般固有ベクトル、および
・正の一般固有値と正のＪ−ｎｏｒｍ  ｖＴＪｖ＝＋１を有する二つの一般固有ベクトル。
【０２０９】
我々が求めている一般固有ベクトルは、ｖＴＪｖが＋１に等しく、ｖＴＣＴＣｖが最小であるようなものであり、このことはすなわちそれが（ＣＴＣ，Ｊ）の最小の正一般固有値に対応することを意味している。
【０２１０】
要約すると、最適なギブンスおよびハイパボリック角θおよびφは、以下の等式により定義される：
【０２１１】
【数５６】
2011509037.tif 000058
【０２１２】
なお式中、ｖ＝（ｘ，ｙ，ｚ）Ｔは、最小の正固有値の（ＣＴＣ，Ｊ）の一般固有ベクトルである。
【０２１３】
固有ベクトルがその符号次第で決定されることから、ｚ＞０すなわちｃｏｓ２θが正となるように、ｖを選択することが常に可能である。その上、θをθ＋пで置換した場合でも、式（１７）は不変である。
【０２１４】
したがってθは、ｃｏｓθ＞０となるように−п／２と−п／２の間で選択できる。
こうして、ｃｏｓθとｃｏｓ２θの両方が正であると仮定することができ、これは、−п／４≦θ≦п／４を仮定することと等価である。
【０２１５】
角度θおよびφの明示的表現が不要であるという点に留意されたい。更新（１７）をインプリメンテーションするにはｃｏｓｈφ、ｓｉｎφ、ｃｏｓθおよびｓｉｎθの計算のみが必要である。単純な三角法により（２２）の以下の解が直ちに得られる：
【０２１６】
【数５７】
2011509037.tif 000059
【０２１７】
【数５８】
2011509037.tif 000060
【０２１８】
【数５９】
2011509037.tif 000061
【０２１９】
【数６０】
2011509037.tif 000062
【０２２０】
この段階で、ギブンス回転Ｇ（θ，ｉ，ｊ）およびハイパボリック回転Ｈ（φ，ｉ，ｊ）を用いて所与のインデックスの対１≦ｉ＜ｊ≦ＮについてＲ’ｋ［ｉ，ｊ］を全体的に最小化するための方法が知られる。これらの最適な角度θおよびφがひとたび計算されたならば、セットＲの全ての行列は（１７）により更新され、これは、Ｒｋの行または列ｉおよびｊのわずか二つの線形結合を計算することによって達成することができる。
【０２２１】
Ｎ×Ｎ恒等行列で初期化された混合行列は、ＡＲｋＡＴ＝Ａ’Ｒ’ｋＡ’Ｔを保つために以下の式によって更新される。
【０２２２】
【数６１】
2011509037.tif 000063
【０２２３】
【数６２】
2011509037.tif 000064
【０２２４】
更新の後、Ａの列のノルムはもはや等しくなく、行列Ｄ（ｉ，ｊ，λ）で再正規化され得る。次節では、この正規化を各ステップで行なう必要のないことが説明されている。ＲｋおよびＡの更新は、ヤコビアルゴリズム［８］の場合と同様に、１≦ｉ＜ｊ≦Ｎで各々のインデックスの対（ｉ，ｊ）について反復される。Ｎ（Ｎ−１）／２つのインデックスの対にわたる完全な更新は、スイープと呼ばれる。古典的なヤコビアルゴリズムと類似して、完全な収束に至るまで複数のスイープが一般に必要である。
【０２２５】
Ｖ．Ｂ．  インプリメンテーションの詳細
その他の回転シーケンスが考えられる。ギブンス回転は直交であるため、これらの回転はハイパボリック回転に比べ発生させる丸めノイズが少ない。したがって、例えば非直交ハイパボリック回転を使用する前に、純粋なギブンス回転更新の複数のスイープ（すなわちＪＡＤＥスイープ）により結合対角化を開始することができると思われる。この事前準備は、たとえ悪コンディションの混合行列についてさえ（最高１０００００の条件数）、我々の数値的実験において不要であった。数値的精度の喪失は恐らく、大きい角度φでのハイパボリック回転に結びつけられる：これらの事象も同様に監視可能であると思われる。
【０２２６】
ＪＡＤＥまたはＳＯＢＩ予備白色化は不要であるが、それでもＲが正値定符号行列を含む場合、数値精度を改善するために使用することができる。ＢＳＳアプリケーションにおいては、ソースの数がセンサーの数よりも少ないという予備的知識が知られている場合、わずかな主成分の部分空間へのデータの削減を明らかに計算することができる。
【０２２７】
式（２３）、（２４）、（２５）、（２６）は、（２２）の考えられる唯一の反転公式（ｉｎｖｅｒｓｉｏｎ  ｆｏｒｍｕｌａｅ）ではない。これらは、ｖ１およびｖ２が０まで収束する時のその優れた数値的精度のために選択された。より良い選択肢が存在する可能性もある。収束に近い場合、行列ＣおよびＣＴＣは特異点に近い。このとき、（ＣＴＣ，Ｊ）の一般固有分解は数値的に困難であり得、時として標的一般固有ベクトルはもはや、より小さい正の一般固有値に対応せず、非常に小さい負の固有値に対応する。より大きな正の固有値が誤って選択されるのを回避するため、より小さい正の固有値の代りにメジアン一般固有値に付随する固有ベクトルを選択することが好ましい。
【０２２８】
以下のシミュレーションでインプリメンテーションが行われる。固有値符号が正しくない場合（＋、＋、−）に更新を放棄することまたは標的固有ベクトルのみを計算するための具体的手順を開発することといったその他の解決法を調査することができると思われる。
【０２２９】
シミュレーションにおいては、Ａの列の正規化が各ステップにおいて必要でないが、一回のスイープにつき一回のみ実施してよいということが観察されている。それは、以下の式で、右側では対角行列ＤＡをＡに乗算し、両側からＤ−１ＡをＲｋに乗算する（式（２）および（３）と同様）ことによって簡単に達成される。
【０２３０】
【数６３】
2011509037.tif 000065
【０２３１】
なお式中、ａｎはＡのｎ番目の列であり、｜｜ａｎ｜｜はそのユークリッド・ノルムである。
【０２３２】
各スイープにおけるこの正規化（２ＫＮ２）のコストは、Ｒｋ（２ＫＮ３）の更新（１７）に比較して無視できるものである。シミュレーションを考慮すると、本発明はなおも、列のノルムのいかなる正規化もなく収束する。
【０２３３】
非対角項のノルムを監視することおよびその最小値での安定化を検出することといった複数の停止基準を想像することができる。例えば当該ケースにおいては、スイープの全ての三角法およびハイパボリックサインが非常に小さい場合、反復を停止することが選択される（Ｍａｔｌａｂの二倍精度で正確に同時対数化することが可能な行列および計算のためには、１０−２０〜１０−２５の間の閾値が適切であると思われる：この値は、正確に結合対角化することが可能でないセットについては１０−１０〜１０−２０の間で増大されなければならない）。
【０２３４】
Ｖ．Ｃ．  Ｊ−Ｄｉアルゴリズムのまとめ
次節の数値的シミュレーションはＪ−Ｄｉの以下のバージョンをインプリメンテーションする。
要件：Ｋ個の対称Ｎ×Ｎ行列Ｒ１，…ＲＫおよび閾値Τ
【０２３５】
【数６４】
2011509037.tif 000066
【０２３６】
該アルゴリズムの計算コストは、主としてインデックスｉおよびｊの各対についてＲｋ行列を更新するコスト、すなわち（１７）のコストである。ＧおよびＨの疎性とＲｋの対称性を用いることで、この更新のコストを４ＫＮフロップ（ｆｌｏｐｓ）に削減することができる。Ｎ（Ｎ−１）／２個のインデックスの対の全体にわたる一回のスイープのためのコストは、このときおおよそ２ＫＮ３である。収束に至るまでに行なうべきスイープの数は、次節で記述する通り、手元のＮＯＪＤ問題のサイズと難度により左右される。
【０２３７】
Ｖ．Ｄ．  Ａｆｓａｒｉアルゴリズムとの簡潔な比較
Ｊ−Ｄｉと［７］で提案されている四つのアルゴリズムとの間にあるいくつかの類似性は、強調するに値する。
−両方の方法が共に、平面変換（ｐｌａｎａｒ  ｔｒａｎｓｆｏｒｍａｔｉｏｎ）の積として、対角化手段を構築する、すなわち、ＱＲＪ１ＤおよびＱＲＪ２Ｄにおけるギブンス回転および三角行列、またはＬＵＪ１ＤおよびＬＵＪ２Ｄ方法における三角行列のみ、そしてＪ−Ｄｉの場合のギブンス回転およびハイパボリック回転。全てのケースにおいて、これらの変換の行列式は１であり、それらの積は特殊線形群にまたがる。
−ＬＵＪ１Ｄ、ＬＵＪ２Ｄ、ＱＲＪ１ＤおよびＱＲＪ２Ｄは、グローバルコスト関数を最小化し、Ｊ−Ｄｉは、各々の現行のインデックスの対（ｉ；ｊ）のローカル費用関数を最小化する。
−各々のインデックスの対の二つの角度はＪ−Ｄｉにより同時に最適化され、一方で、二つの最適化はＬＵＪ１Ｄ、ＬＵＪ２Ｄ、ＱＲＪ１ＤおよびＱＲＪ２Ｄにおいて独立したものである。
−（必然的に）限定されたシミュレーションを考慮すると、ＬＵＪ１Ｄ、ＬＵＪ２Ｄ、ＱＲＪ１ＤおよびＱＲＪ２Ｄは、一次的に収束し、Ｊ−Ｄｉは二次的に収束する。
−スイープ一回あたりの計算コストは（二つの（ｉ；ｊ）−更新をグループにした場合）同一である。
ＬＵＪ１Ｄは、特定のシミュレーションされたコンテクスト内で比較的速いと思われることから、次節において比較のために選択される。
【０２３８】
ＶＩ．  数値的シミュレーションについての性能評価
ＦＦＤｉａｇ［４］およびＬＵＪ１Ｄ［７］という二つの終りのＮＯＪＤアルゴリズムがＲの任意の行列の正値符号性（ｐｏｓｉｔｉｖｅ  ｄｅｆｉｎｉｔｅｎｅｓｓ）を必要としないことから、これら二つのアルゴリズムと本発明に係る方法（すなわちＪ−Ｄｉ）を比較する。
【０２３９】
ＦＦＤｉａｇアルゴリズムの１の反復（以下でスイープと呼ぶ）は、一つのＪ−ＤｉまたはＬＵＪ１Ｄスイープと比べてほぼ同じ複雑性を有する（Ｒを更新するために２ＫＮ３）。したがって、唯一比較されるのは、収束に必要なスイープ数だけである。
【０２４０】
それでも、一つのスイープ内で演算は異なる形で構造化されている。すなわち、Ｊ−ＤｉおよびＬＵＪ１ＤはＲの非常に疎な更新を必要とし、一方でＦＦＤｉａｇは、完全な高密度の更新のみを必要とする。
【０２４１】
この差がおそらく、Ｍａｔｌａｂのようなインタプリタ型言語にインプリメンテーションされた場合のＪ−ＤｉまたはＬＵＪ１Ｄスイープに比べ、ＦＦＤｉａｇのスイープをわずかに速いものにしている。ただし、例えば、Ｃ言語によるインプリメンテーションは類似の実行時間を示すはずである。
【０２４２】
二つの異なる性能の基準が定義され、最初に正確に対角化可能な行列のセットについて測定され、第二にほぼ対角化可能な行列のセットについて測定される。シミュレーションは、Ｍａｔｌａｂの二倍精度で実行される。
【０２４３】
ＶＩ．Ａ．  非直交結合対角化についての性能の基準
以下では結合対角化の質を測定する二つの方法が用いられる。最初の基準は、Ｒの非対角エントリの二乗の和として定義される非対角ノルムとＲの対角エントリの二乗の和として定義される対角ノルムの比率を計算することから成る。この非標準正規化は、対角ノルムが最初のスイープの間一定でないことから有用である。
【０２４４】
第二の基準は、混合行列の推定の質と、対応する分離性能を査定するのに必要である。実際のところ、下記の〔数６５〕の推定された混合行列のインバースは、Ａによって混合されたソースを分離するのに（付加ノイズが無い場合に）最適な空間フィルタである。なお、本明細書においては、明細書の記載として使用できる文字の制約を考慮し、
【０２４５】
【数６５】
2011509037.tif 000067
【０２４６】
上記〔数６５〕を、テキスト上、「Ａ：」と表記する。
【０２４７】
積Ａ：−１Ａが対角行列Ｄの置換行列Π倍に近いほど、分離は優れたものとなる。したがって、
【０２４８】
【数６６】
2011509037.tif 000068
【０２４９】
と、置換行列と対角行列の最も近い積ΠＤとの間の距離を計算することによって、分離の質を測定するのが自然である。
【０２５０】
ＢＳＳにおいては、たとえそれがΠＤまでの完全な距離でなくても、
【０２５１】
【数６７】
2011509037.tif 000069
【０２５２】
によって定義される性能Ｉ（Ｐ）のインデックスを使用することが標準的である。正規化項１／２Ｎ（Ｎ−１）は、Ｐが対角エントリが１に等しく、非対角エントリがαに等しい行列である場合に、０≦α＜１においてＩ（Ｐ）がαに等しくなるような形で選択された。したがって、Ｉ（Ｐ）は、分離されたソースにおける相対的な干渉振幅（ａｍｐｌｉｔｕｄｅ  ｏｆ  ｉｎｔｅｒｆｅｒｅｎｃｅｓ）の桁（ｏｒｄｅｒ  ｏｆ  ｍａｇｎｉｔｕｄｅ）を提供する（クロストークインデックス（ｃｒｏｓｓ−ｔａｌｋ  ｉｎｄｅｘ））。
【０２５３】
ＶＩ．Ｂ．  正確に対角化可能な行列セット
混合行列エントリは独立して正規分布しており、その行列式は１に設定され、その列は、等式（３）を用いて同じノルムに正規化される。これらのランダム混合行列は一般的に、以下の図面で示されている通り、１０４または１０５未満の条件数を有し、それより多いことはまれである。Ｊ−Ｄｉアルゴリズムのより特定的な妥当性確認は、より悪条件の混合行列が関与するアプリケーションのために行なわれるべきである。Ｄｋの行列の対角エントリは、同じく、標準的な正規分布から選択される。このようなＤｋ生成方法が、ＫおよびＮのみに左右される一意性のモジュラス（ｍｏｄｕｌｕｓ  ｏｆ  ｕｎｉｑｕｅｎｅｓｓ）の分布を決定するという点に留意されたい。ＲｋはＲｋ＝ＡＤｋＡＴによって計算される。
【０２５４】
１）二次収束：
Ｊ−Ｄｉ、ＦＦＤｉａｇおよびＬＵＪ１Ｄアルゴリズムの対角性比は、収束に至るまで各スイープ後に計算され、図２には、Ｎ＝Ｋ＝５０および５０の独立した試行についてプロットされている。この場合、一意性のモジュラスは小さく、ほぼ常に０．４〜０．６の間にある。
【０２５５】
分離基準は、同じ条件下で図３にプロットされている。
【０２５６】
これらの結果は、ＦＦＤｉａｇおよびＪ−Ｄｉアルゴリズムが（少なくとも）二次収束を有すると思われるが、Ｊ−Ｄｉはより少ないスイープしか必要としないことを示している。最終的な一次収束にも関わらず、ＬＵＪ１Ｄは中間の性能を示す。
【０２５７】
２）Ｊ−Ｄｉ収束速度に対する行列およびセットサイズの影響
異なる混合行列および異なるセットサイズ（Ｎ＝Ｋ＝５、１０、２０、５０、１００および２００）について、Ｊ−ＤＩの収束に達するのに必要なスイープ数が計数され、１００回の独立した試行についてのその分布が、下表に示されている：
【０２５８】
【表１】
2011509037.tif 000070
【０２５９】
Ｊ−Ｄｉスイープ数は行列サイズと共に増大するものの、大きい行列（Ｎ＝２００）の大きなセット（Ｋ＝２００）においても、ほぼ１０未満にとどまることがわかる。
【０２６０】
ＶＩ．Ｃ．  正確に対角化できない行列セット：ブラインドソース分離に対する応用
【０２６１】
Ｊ−Ｄｉ、ＦＦＤｉａｇおよびＬＵＪ１Ｄ非直交結合対角化アルゴリズムの性能を、正確に対角化できない行列セットについて比較する。
【０２６２】
（共分散行列のように）正でもなく負でもない行列セットについて本発明に係る方法（Ｊ−Ｄｉ）を評価するため、四次キュミュラントスライスセットが選択される。
【０２６３】
古典的ＢＢＳモデルすなわちｘ＝Ａｓ＋ｎを考慮する。なお式中、ｘ、ｓおよびｎはそれぞれに観察された混合体のＮ次元ベクトル、ソースおよび付加的ガウスノイズ（ａｄｄｉｔｉｖｅ  Ｇａｕｓｓｉａｎ  ｎｏｉｓｅ）を表わし、ＡはＮ×Ｎの混合行列である。考慮対象の行列セットＲは、
【０２６４】
【数６８】
2011509037.tif 000071
【０２６５】
によって定義されるサイズＮ×ＮのＫ＝Ｎ（Ｎ＋１）／２個の行列Ｒｍｎを含んでいる。ここで、全てのｍおよびｎが１≦ｍ≦ｎ≦Ｎであり、式中ｘｎはｘのｎ番目のエントリを表わす。
【０２６６】
Ｎ＝２０ソースであり、したがってＫ＝Ｎ（Ｎ＋１）／２＝２１０の２０×２０サイズのキュムラントスライス（行列）であると考えられる。ソース信号のサンプルは独立しており、−√３と√３の間で一様分布して非ガウス、ゼロ平均（ｚｅｒｏ−ｍｅａｎ）かつユニティバリアンスな（ｕｎｉｔｙ−ｖａｒｉａｎｃｅ）信号を提供する。ゼロ平均ガウス付加ノイズ（ｚｅｒｏ−ｍｅａｎ  Ｇａｕｓｓｉａｎ  ｎｏｉｓｅ）σの標準偏差は、まずゼロに設定され、第二にσ＝０．１に設定される。キュムラントは、１０００００個の独立したサンプルについて推定される。これらのケースにおいて、行列数Ｋは充分大きいため、一意性のモジュラス（ｍｏｄｕｌｕｓ  ｏｆ  ｕｎｉｑｕｅｎｅｓｓ）は一般的に０．３〜０．６の間という小さいものとなる。このパラメータがシミュレーションの性能に影響を及ぼしそうにない理由がこれによって説明されており、したがってこれは提示されない。
【０２６７】
収束速度
対角性比および分離基準の収束は、図４にノイズを伴うスイープ数の関数としてプロットされている。
【０２６８】
全てのケースにおいて、本発明に係る方法（Ｊ−Ｄｉ）は、［７］の場合と同様に挙動すると思われるＦＦＤｉａｇおよびＬＵＪ１Ｄよりも速く収束する。Ｊ−Ｄｉは、ＢＳＳアプリケーションにとって重要である最終的な分離基準の観点からみて、ＦＦＤｉａｇおよびＬＵＪ１Ｄよりもわずかに性能が優れている。
【０２６９】
これらのシミュレーションは同様に、ＦＦＤｉａｇは二次でありＬＵＪ１Ｄは一次である異なる収束速度が、行列セットが正確に対角化可能でない場合の類似の性能と適合可能であることをも示している。
【０２７０】
ＶＩＩ．  ＷＡＮＧ．ＬＩＵおよびＺＨＡＮＧのアルゴリズムとの比較
Ｗａｎｇ、ＬｉｕおよびＺｈａｎｇ［６］（以下ＤＮＪＤ）のアルゴリズムとの比較について以下で論述する。提案されているアルゴリズムＤＮＪＤは、
【０２７１】
【数６９】
2011509037.tif 000072
【０２７２】
によって定義される行列Ｊ（θ，φ，ｉ，ｊ）の連続的乗算によって混合行列を構築する。
【０２７３】
この定義は明らかに、ギブンス回転およびＪＡＤＥのヤコビのようなアルゴリズムを一般化する試みである。ＤＮＪＤアルゴリズムが本発明と非常に異なるものであることは明白である。例えば、ＤＮＪＤは、本発明に係る方法において用いられる行列式制約もハイパボリック三角法も使用しない。
【０２７４】
ＤＮＪＤにおける主たる困難性は、等式（１６）の反転にある。Ｒｋ行列が正確に対角化できない場合に、非対角エントリの実際の最小化を保証しない近似式を使用することが必要となり、正確な解は提供されない。したがって、（各インデックスの対について）各反復毎に非対角エントリの進化を計算し、非対角エントリが増大した場合更新を放棄することが必要である。残念なことに、この確認の計算コスト４ＫＮは、アルゴリズム自体のコストと同等である（第ＩＩＩ節の終りの５３０３頁を参照のこと）。したがって、各スイープのコストは二倍となる。さらに、二つ一組の角度（θ，φ）が放棄される毎に、数多くの計算（８ＫＮフロップ）が費用関数の改善無く行なわれる。この無駄が、Ｊ−Ｄｉに比べてＤＮＪＤが用いるスイープ数の多さを説明し得ると思われる（［６］中の５３０３項、図１）（表１）。
【０２７５】
要するに、Ｊ−Ｄｉはスイープ毎の複雑性の観点から見てＤＮＪＤよりも二倍以上優れた性能を示す。
【０２７６】
ＶＩＩＩ．  結論
非直交結合対角化の新しい方法（Ｊ−Ｄｉ）が提示された。一つの実施態様において、この方法は、以下のステップを含む：
−特殊な直交群の中で構築する代わりに、列が等ノルムであり、また、特殊線形群に属する混合行列を構築するステップ；
−ギブンス回転だけでなく、ギブンス回転、ハイパボリック回転および対角行列を乗算するステップ。
【０２７７】
この新しい方法は、数多くの変形形態すなわち基本行列の乗算の次数、停止テスト、予備白色化の有無などを許容する。そのロバスト性および削減された複雑性により、標準コンピュータ上で最高Ｎ＝５００およびＫ＝６００までの大規模な問題を扱うことが可能になる。
【符号の説明】
【０２７８】
１０  処理手段
２０  表示手段
３０  入力手段
【先行技術文献】
【非特許文献】
【０２７９】
【非特許文献１】［１］Ｃｏｍｐａｒａｉｓｏｎ  ｄｅ  ｍｅｔｈｏｄｅｓ  ｄｅ  ｓｅｐａｒａｔｉｏｎｓ  ｄｅ  ｓｏｕｒｃｅｓ，Ｓｏｕｌｏｕｍｉａｃ  Ａ．，Ｃａｒｄｏｓｏ  Ｊ．−Ｆ．，Ｐｒｏｃｅｅｄｉｎｇｓ  ｏｆ  ＧＲＥＴＳＩ  １９９１，Ｊｕａｎ  ｌｅｓ  Ｐｉｎｓ，Ｓｅｐｔｅｍｂｒｅ  １９９１．
【非特許文献２】［１ｂ］Ｂｌｉｎｄ  ｂｅａｍｆｏｒｍｉｎｇ  ｆｏｒ  ｎｏｎ  Ｇａｕｓｓｉａｎ  ｓｉｇｎａｌｓ，Ｊｅａｎ−Ｆｒａｎｃｏｉｓ  Ｃａｒｄｏｓｏ  ａｎｄ  Ａｎｔｏｉｎｅ  Ｓｏｕｌｏｕｍｉａｃ，ｉｎ  ＩＥＥ  Ｐｒｏｃｅｅｄｉｎｇｓ−Ｆ，１４０（６）：３６２−３７０，Ｄｅｃｅｍｂｅｒ  １９９３．
【非特許文献３】［２］Ａ  ｂｌｉｎｄ  ｓｏｕｒｃｅ  ｓｅｐａｒａｔｉｏｎ  ｔｅｃｈｎｉｑｕｅ  ｕｓｉｎｇ  ｓｅｃｏｎｄ−ｏｒｄｅｒ  ｓｔａｔｉｓｔｉｃｓ，Ｂｅｌｏｕｃｈｒａｎｉ  Ａ．，ＡｂｅｄＭｅｒａｉｍ  Ｋ．，Ｃａｒｄｏｓｏ  Ｊ．−Ｆ．，Ｍｏｕｌｉｎｅｓ  Ｅ．，ＩＥＥＥ  ＴＲＡＮＳＡＣＴＩＯＮＳ  ＯＮ  ＳＩＧＮＡＬ  ＰＲＯＣＥＳＳＩＮＧ  ４５（２）：４３４−４４４，ＦＥＢ  １９９７．
【非特許文献４】［３］Ｊｏｉｎｔ  Ａｐｐｒｏｘｉｍａｔｅ  Ｄｉａｇｏｎａｌｉｚａｔｉｏｎ  ｏｆ  Ｐｏｓｉｔｉｖｅ  Ｄｅｆｉｎｉｔｅ  Ｈｅｒｍｉｔｉａｎ  Ｍａｔｒｉｃｅｓ，Ｐｈａｍ  Ｄ．Ｔ．，ＳＩＡＭ  Ｊ．ｏｎ  Ｍａｔｒｉｘ  Ａｎａｌ．ａｎｄ  Ａｐｐｌ．，２２（４）：１１３６−１１５２，２００１．
【非特許文献５】［４］Ａ  Ｆａｓｔ  Ａｌｇｏｒｉｔｈｍ  ｆｏｒ  Ｊｏｉｎｔ  Ｄｉａｇｏｎａｌｉｚａｔｉｏｎ  ｗｉｔｈ  Ｎｏｎ−ｏｒｔｈｏｇｏｎａｌ  Ｔｒａｎｓｆｏｒｍａｔｉｏｎｓ  ａｎｄ  ｉｔｓ  Ａｐｐｌｉｃａｔｉｏｎ  ｔｏ  Ｂｌｉｎｄ  Ｓｏｕｒｃｅ  Ｓｅｐａｒａｔｉｏｎ，Ｚｉｅｈｅ  Ａ．，Ｌａｓｋｏｖ  Ｐ．，Ｎｏｌｔｅ  Ｇ．，Ｍｕｌｌｅｒ  Ｋ．−Ｒ．，Ｊｏｕｒｎａｌ  ｏｆ  Ｍａｃｈｉｎｅ  Ｌｅａｒｎｉｎｇ  Ｒｅｓｅａｒｃｈ  ５（２００４）７７７−８００．
【非特許文献６】［５］Ｑｕａｄｒａｔｉｃ  Ｏｐｔｉｍｉｚａｔｉｏｎ  ｆｏｒ  Ｓｉｍｕｌｔａｎｅｏｕｓ  Ｍａｔｒｉｘ  Ｄｉａｇｏｎａｌｉｚａｔｉｏｎ，Ｖｏｌｌｇｒａｆ  Ｒ．，Ｏｂｅｒｍａｙｅｒ  Ｋ．，ＩＥＥＥ  Ｔｒａｎｓ，ｏｎ  Ｓｉｇｎａｌ  Ｐｒｏｃｅｓｓｉｎｇ，ｖｏｌ．５４，ｎｏ．９，ｐｐ．３２７０−３２７８，Ｓｅｐｔｅｍｂｅｒ  ２００６．
【非特許文献７】［６］Ｎｏｎｏｒｔｈｏｇｏｎａｌ  Ｊｏｉｎｔ  Ｄｉａｇｏｎａｌｉｚａｔｉｏｎ  Ａｌｇｏｒｉｔｈｍ  Ｂａｓｅｄ  ｏｎ  Ｔｒｉｇｏｎｏｍｅｔｒｉｃ  Ｐａｒａｍｅｔｅｒｉｚａｔｉｏｎ，Ｗａｎｇ  Ｆ．，Ｌｉｕ  Ｚ．，Ｚｈａｎｇ  Ｊ．，ＩＥＥＥ  Ｔｒａｎ．ｏｎ  Ｓｉｇｎａｌ  Ｐｒｏｃｅｓｓｉｎｇ，ｖｏｌ．５５，ｎｏ．１１，ｐｐ．５２９９−５３０８，Ｎｏｖｅｍｂｅｒ  ２００７．
【非特許文献８】［７］Ｓｉｍｐｌｅ  ＬＵ  ａｎｄ  ＱＲ  Ｂａｓｅｄ  Ｎｏｎ−Ｏｒｔｈｏｇｏｎａｌ  Ｍａｔｒｉｘ  Ｊｏｉｎｔ  Ｄｉａｇｏｎａｌｉｚａｔｉｏｎ，Ｂ．ＡＦＳＡＲＩ，Ｉｎｄｅｐｅｎｄｅｎｔ  Ｃｏｍｐｏｎｅｎｔ  Ａｎａｌｙｓｉｓ  ａｎｄ  Ｂｌｉｎｄ  Ｓｏｕｒｃｅ  Ｓｅｐａｒａｔｉｏｎ，Ｌｅｃｔｕｒｅ  Ｎｏｔｅｓ  ｉｎ  Ｃｏｍｐｕｔｅｒ  Ｓｃｉｅｎｃｅ，ｖｏｌ．３８８９，ｐｐ．１−７，２００６．
【非特許文献９】［８］Ｍａｔｒｉｘ  Ｃｏｍｐｕｔａｔｉｏｎｓ，Ｇ．Ｈ．ＧＯＬＵＢ，Ｃ．Ｆ．ＶＡＮ  ＬＯＡＮ，Ｂａｌｔｉｍｏｒｅ，Ｊｏｈｎｓ  Ｈｏｐｋｉｎｓ  Ｕｎｉｖｅｒｓｉｔｙ  Ｐｒｅｓｓ，３ｒｄ  ｅｄｉｔｉｏｎ，１９９６．
【非特許文献１０】［９］Ｗｈａｔ  Ｃａｎ  Ｍａｋｅ  Ｊｏｉｎｔ  Ｄｉａｇｏｎａｌｉｚａｔｉｏｎ  Ｄｉｆｆｉｃｕｌｔ？，Ｂ．ＡＦＳＡＲＩ，ＩＣＡＳＳＰ  ２００７，Ｖｏｌ．３，ｐｐ．ＩＩＩ−１３７７−１３８０，１５−２０  Ａｐｒｉｌ  ２００７．
【非特許文献１１】［１０］Ｓｅｎｓｉｔｉｖｉｔｙ  Ａｎａｌｙｓｉｓ  ｆｏｒ  ｔｈｅ  ｐｒｏｂｌｅｍ  ｏｆ  Ｎｏｎ−Ｏｒｔｈｏｇｏｎａｌ  Ｍａｔｒｉｘ  Ｊｏｉｎｔ  Ｄｉａｇｏｎａｌｉｚａｔｉｏｎ，Ｂ．ＡＦＳＡＲＩ，Ｓｕｂｍｉｔｔｅｄ  ｔｏ  ＳＩＡＭ  Ｊ．ｏｎ  Ｍａｔｒｉｘ  Ｃｏｍｐｕｔａｔｉｏｎｓ．
【非特許文献１２】［１１］Ｎｏｎ−Ｏｒｔｈｏｇｏｎａｌ  Ｊｏｉｎｔ  Ｄｉａｇｏｎａｌｉｚａｔｉｏｎ  ｉｎ  ｔｈｅ  Ｌｅａｓｔ−Ｓｑｕａｒｅｓ  Ｓｅｎｓｅ  ｗｉｔｈ  Ａｐｐｌｉｃａｔｉｏｎ  ｉｎ  Ｂｌｉｎｄ  Ｓｏｕｒｃｅ  Ｓｅｐａｒａｔｉｏｎ，Ａ．ＹＥＲＥＤＯＲ，ＩＥＥＥ  Ｔｒａｎｓ．ｏｎ  Ｓｉｇｎａｌ  Ｐｒｏｃｅｓｓｉｎｇ，ｖｏｌ．５０，ｎｏ．７，ｐｐ．１５４５−１５５３，Ｊｕｌｙ  ２００２．
【非特許文献１３】［１２］Ｓｕｒ  ｌａ  ｄｉａｇｏｎａｌｉｓａｔｉｏｎ  ｃｏｎｊｏｉｎｔｅ  ａｐｐｒｏｃｈｅｅ  ｐａｒ  ｕｎ  ｃｒｉｔｅｒｅ  ｄｅｓ  ｍｏｉｎｄｒｅｓ  ｃａｒｒｅｓ，Ｓ．ＤＥＧＥＲＩＮＥ，ｉｎ  Ｐｒｏｃｅｅｄｉｎｇｓ  ＸＶＩＩＩ  Ｃｏｌｌｏｑｕｅ  Ｇｒｅｔｓｉ  ２００１，ｎ２，ｐｐ．３１１−３１４，Ｔｏｕｌｏｕｓｅ，Ｆｒａｎｃｅ，Ｓｅｐｔ．２００１．
【図１】
2011509037.tif 000076
【図２】
2011509037.tif 000077
【図３】
2011509037.tif 000078
【図４】
2011509037.tif 000079
【図５】
2011509037.tif 000080
【国際調査報告】
2011509037.tif 000081
2011509037.tif 000082
2011509037.tif 000083
2011509037.tif 000084
2011509037.tif 000085
2011509037.tif 000086