カーネル主成分分析方法、カーネル主成分分析装置、カーネル主成分分析プログラム
【課題】計算量を削減することができ、かつ、計算量の削減に伴う計算精度の低下が従来よりも小さなカーネル主成分分析方法を実現する。
【解決手段】本発明の主成分分析方法は、ｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝からｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝を選択する選択ステップと、一般化固有値問題ＫxyTＫxyｚ＝λＫyzを解くことによって、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出する設計ステップと、評価対象ベクトルｘの特徴ベクトルｙ＝［ｚ1，ｚ2，…，ｚr］TｘΦを算出する評価ステップとを含んでいる。
【解決手段】本発明の主成分分析方法は、ｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝からｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝を選択する選択ステップと、一般化固有値問題ＫxyTＫxyｚ＝λＫyzを解くことによって、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出する設計ステップと、評価対象ベクトルｘの特徴ベクトルｙ＝［ｚ1，ｚ2，…，ｚr］TｘΦを算出する評価ステップとを含んでいる。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、主成分分析に関するものであり、特に、カーネル主成分分析に関するものである。
【背景技術】
【０００２】
多変量解析のひとつとして、主成分分析（ＰＣＡ：principal component analysis）がある。主成分分析は、ｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝（ｄ次元ベクトル）に基づき、ｄ次元の評価対象ベクトルｘの特徴を表すｒ次元の特徴ベクトルｙを算出する方法であり、（１）ｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝から定義される分散共分散行列の固有値問題を解き、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出する設計ステップと、（２）ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を用いて評価対象ベクトルｘを特徴ベクトルｙに変換する評価ステップとにより構成することができる。得られる特徴ベクトルｙの各成分は、評価対象ベクトルｘの各主成分への射影を表し、「主成分得点」と称される。
【０００３】
この主成分分析の拡張として、カーネル主成分分析（ＫＰＣＡ：Kernel principal component analysis）がある（非特許文献１参照）。カーネル主成分分析は、特徴空間と呼ばれる再生核ヒルベルト空間Ｆ上での主成分分析であり、例えば（１）ｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝から定義されるカーネルグラム行列Ｋｘの固有値値問題を解き、ｒ個の固有ベクトル固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出する設計ステップと、（２）ｒ個の固有ベクトル固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を用いて評価対象ベクトルｘを特徴ベクトルｙに変換する評価ステップとにより構成することができる。ここで、カーネルグラム行列Ｋｘとは、ｋ（ｘ、ｙ）＝＜Φ（ｘ）｜Φ（ｙ）＞をＭｅｒｃｅｒのカーネル関数として、カーネル関数値ｋ（ｘi，ｘj）を第ｉｊ成分とするｎ×ｎ行列のことである。Φ（ｘ）は入力空間Ｒdから特徴空間Ｆへの写像を表し、＜・｜・＞は内積を表す。
【０００４】
カーネル主成分分析を用いることにより、非線形相関を考慮した多変量解析を行うことができる。しかしながら、ｎ×ｎ行列の固有値問題を解く必要があるため、計算量（計算の複雑さ）はｎ3に比例し、必要なメモリサイズはｎ2に比例する。したがって、現在のコンピュータの性能では、標本数ｎが数万から数十万を超えると、カーネル主成分分析を行うことが困難になる。このため、与えられたｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝の一部のみを用いたカーネル主成分分析しか行えないことがしばしばであった。
【０００５】
このような問題を考慮したカーネル主成分分析の改良としては、改良カーネル主成分分析（ＩＫＰＣＡ：Improved Kernel principal component analysis）（非特許文献２参照）や、スパースカーネル主成分分析（非特許文献３、４）などが知られている。
【先行技術文献】
【非特許文献】
【０００６】
【非特許文献１】B. Scholkopf and A. Smola and K.-R. Muller，“Nonlinear Component Analysis as a Kernel Eigenvalue Problem”，Neural computation，10，5，pp. 1299-1319，1998
【非特許文献２】Y. Xu, D. Zhang, F. Song, J. Yang, Z. Jing and M. Li，“Amethod for speeding up feature extraction based on KPCA”，Neurocomputing， pp. 1056-1061，2007
【非特許文献３】M. E. Tipping，“Sparse kernel principal component analysis”，Advances in Neural Information Processing Systems(NIPS)，13，pp. 633-639，2001
【非特許文献４】A. J. Smola, O. L. Mgngasarian and B. Scholkopf，“Sparse kernel feature analysis”，Technical report 99-04，Universityof Wisconsin，1999
【発明の概要】
【発明が解決しようとする課題】
【０００７】
しかしながら、上記従来の技術には、計算量を削減すると計算精度が大幅に低下するという問題があった。
【０００８】
すなわち、与えられたｎ個の標本ベクトルの一部のみを用いてカーネル主成分分析を行えば、計算量も標本数の減少に応じて当然に減少する。しかし、用いる標本ベクトルの数を減らすと、得られる特徴ベクトルの精度低下を招来する。例えば、１０００個の標本ベクトルが与えられたときに、５０個の標本ベクトルを用いたカーネル主成分分析を行うと、得られる特徴ベクトルの精度が大幅に低下する。
【０００９】
また、改良カーネル主成分分析は、適切な問題設定ではないため、近似が正確でないという問題を有している。また、スパースカーネル主成分分析は、評価ステップにおける計算量を削減することを目的とするものであり、設計ステップにおける計算量を削減することができないという問題を有している。
【００１０】
本発明は、上記の問題に鑑みてなされたものであり、その目的は、計算量を削減することができ、かつ、計算量の削減に伴う計算精度の低下が従来よりも小さなカーネル主成分分析方法を実現することにある。
【課題を解決するための手段】
【００１１】
本発明に係るカーネル主成分分析方法は、カーネル主成分分析装置を用いて、標本として与えられたｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝に基づき、評価対象として与えられたｄ次元ベクトルｘの特徴を表現するｒ次元ベクトルｙを算出するカーネル主成分分析方法である。上記課題を解決するために、本発明に係るカーネル主成分分析方法は、上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝からｍ個（ｍはｎ未満の自然数）のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を選択する選択ステップ、又は、上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝とは別にｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を取得する取得ステップと、カーネル関数値ｋ（ｙi，ｙj）を第ｉｊ成分とするｍ×ｍ行列をＫyとし、カーネル関数値ｋ（ｘi，ｙj）を第ｉｊ成分とするｎ×ｍ行列をＫxyとして、固有値λと固有ベクトルｚを求める一般化固有値問題ＫxyTＫxyｚ＝λＫyz、又は、正則化一般化固有値問題ＫxyTＫxy（ＫxyTＫxy＋μＫy）ＫxyTＫxyｚ＝λＫyzを解くことによって、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出する設計ステップと、上記ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を用いて、上記ｄ次元ベクトルｘを上記ｒ次元ベクトルｙに変換する評価ステップと、を含んでいる。ここで、ｋ（ｘ、ｙ）＝＜Φ（ｘ）｜Φ（ｙ）＞は、従来のカーネル主成分分析に用いられているものと同じMercerのカーネル関数であり、μ＞０は正則化パラメータである。
【００１２】
上記の構成によれば、解くべき固有値問題はｍ×ｍ行列（ＫxyTＫxy及びＫy）に関する一般化固有値問題になる（計算量ｍ3）。このため、ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝を用いた従来のカーネル主成分分析方法（計算量ｎ3）と比べて計算量を減少させることができる。しかも、本発明において解くべき一般化固有値問題には、（選択された）ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝に関する情報のみならず、（選択されたなかった）ｎ−ｍ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝＼｛ｙ1，ｙ2，…，ｙm｝に関する情報がｎ×ｍ行列Ｋxyを介して取り込まれている。したがって、（選択された）ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝のみを用いて従来のカーネル主成分分析を行う場合に比べて、より高い精度で特徴ベクトルｙを算出することができる。したがって、計算量を削減することができ、かつ、計算量の削減に伴う計算精度の低下を従来よりも小さくすることができるという効果を奏する。
【００１３】
なお、ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を用いて従来のカーネル主成分分析を行う場合に比べて、より高い精度の特徴ベクトルｙを算出することができることの、数学的な裏付けについては発明の詳細な説明の〔部分カーネル主成分分析の原理〕を、実験的な検証については発明の詳細な説明の〔数値実験による検証〕をそれぞれ参照されたい。
【００１４】
なお、上記設計ステップは、例えば、（１）上記ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝から、上記ｍ×ｍ行列Ｋyを算出するステップと、（２）上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝、及び、上記ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝から、上記ｎ×ｍ行列Ｋxyを算出するステップと、（３）上記一般化固有値問題ＫxyTＫxyｚ＝λＫyz、又は、ＫxyTＫxy（ＫxyTＫxy＋μＫy）ＫxyTＫxyｚ＝λＫyzを解くことによって、上記ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出するステップとを含んで構成することができる。
【００１５】
ただし、上記設計ステップは、カーネル関数値ｋ（ｙi，ｙj）を第ｉｊ成分とするｍ×ｍ行列をＫyとし、カーネル関数値ｋ（ｘi，ｙj）を第ｉｊ成分とするｎ×ｍ行列をＫxyとしたときに、ＫxyTＫxyｚ＝λＫyzにより表現される一般化固有値問題、又は、ＫxyTＫxy（ＫxyTＫxy＋μＫy）ＫxyTＫxyｚ＝λＫyzにより表現される正則化一般化固有値問題を解くことによって、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出するものであれば何でもよく、必ずしも、上記ｍ×ｍ行列Ｋy自体、及び、上記ｎ×ｍ行列Ｋxy自体の算出を要するものではない。
【００１６】
また、上記評価ステップは、例えば、（１）上記ｄ次元ベクトルｘ、及び、上記ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝から、カーネル関数値ｋ（ｘ，ｙj）を第ｊ成分とするｍ次元ベクトルｘΦを算出するステップと、（２）上記ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝、及び、上記ｍ次元ベクトルｘΦから、上記ｒ次元ベクトルｙ＝［ｚ1，ｚ2，…，ｚr］TｘΦを算出するステップと、を含んで構成することができる。
【００１７】
ただし、上記評価ステップは、カーネル関数値ｋ（ｘ，ｙj）を第ｊ成分とするｍ次元ベクトルをｘΦとし、上記ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を並べたｍ×ｒ行列［ｚ1，ｚ2，…，ｚr］をＺとしたときに、ＺTｘΦにより表現されるｒ次元ベクトルを算出するものであれば何でもよく、必ずしも、上記ｍ次元ベクトルをｘΦ自体、及び、上記ｍ×ｒ行列Ｚ自体の算出を要するものではない。
【００１８】
また、本発明に係るカーネル主成分分析方法は、上記評価ステップにて算出された上記ｒ次元ベクトルｙを他の装置に出力する出力ステップを含んでいてもよい。出力ステップにおける出力先となる装置は特に限定されるものではないが、例えば、上記ｒ次元ベクトルｙを表示装置（ディスプレイ等）、印刷装置（プリンタ等）、記録装置（光ディスクドライブ等）、送信装置（モデム等）を出力先とすることができる。表示装置／印刷装置を出力先とする場合には、上記ｒ次元ベクトルｙを数値列又はグラフとして表示又は印刷することができるし、記録装置を出力先とする場合には、上記ｒ次元ベクトルｙを記録媒体に記録することができるし、送信装置を出力先とする場合には、送信装置を介して上記ｒ次元ベクトルｙを他のコンピュータに送信することができる。これにより、ユーザが上記ｒ次元ベクトルｙを視覚的に認識したり、上記カーネル主成分分析装置または上記記録媒体から情報を読み出すことができる他装置が上記記録媒体から上記ｒ次元ベクトルｙを読み出して利用したり、他のコンピュータが上記送信装置から上記ｒ次元ベクトルｙを受信して利用したりすることができる。
【００１９】
本発明に係るカーネル主成分分析方法において、上記選択ステップは、上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝をｍ個のクラスタ｛Ｃ1，Ｃ2，…，Ｃm｝にクラスタリングするステップと、上記ｍ個のクラスタ｛Ｃ1，Ｃ2，…，Ｃm｝の各々から上記ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝の各々を選択するステップと、を含んでいる、ことが好ましい。
【００２０】
上記構成によれば、ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を選択するために要する計算量を徒に増加させることなく、特徴ベクトルｙを精度良く計算することができるよう、ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を選択することができる、という更なる効果を奏する。
【００２１】
本発明に係るカーネル主成分分析方法は、上記選択ステップにおいて、上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝をｍ個のクラスタ｛Ｃ1，Ｃ2，…，Ｃm｝にクラスタリングする際に、特徴空間上のノルムを使用する、ことが好ましい。
【００２２】
上記の構成によれば、特徴ベクトルｙを更に精度良く計算することができるよう、ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を選択することができる、という更なる効果を奏す。
【００２３】
本発明に係るカーネル主成分分析方法において、上記ベクトルｙの次元ｒは、上記ベクトルｘの次元ｄよりも小さい、ことが好ましい。
【００２４】
上記の構成によれば、伝送効率に優れたデータサイズの小さい特徴ベクトルｙを得ることができる。例えば、本発明を情報源符号化に適用する場合、送信側（装置）と受信側（装置）とが上記ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を予め共有することによって、送受信する情報量を削減することができる。
【００２５】
本発明に係るカーネル主成分分析方法において、上記ベクトルｙの次元ｒは、上記ベクトルｘの次元ｄよりも大きい、ことが好ましい。
【００２６】
上記の構成によれば、ノイズに強い冗長な特徴ベクトルｙを得ることができる。例えば、本発明を通信路符号化に適用する場合、送信側（装置）と受信側（装置）とが上記ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を予め共有することによって、雑音に対して頑健な情報を送受信することができる。
【００２７】
本発明に係るカーネル主成分分析装置は、標本として与えられたｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝に基づき、評価対象として与えられたｄ次元ベクトルｘの特徴を表現するｒ次元ベクトルｙを算出するカーネル主成分分析装置である。そして、上記課題を解決するために、上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝からｍ個（ｍはｎ未満の自然数）のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を選択する選択手段、又は、上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝とは別にｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を取得する取得手段と、カーネル関数値ｋ（ｙi，ｙj）を第ｉｊ成分とするｍ×ｍ行列をＫyとし、カーネル関数値ｋ（ｘi，ｙj）を第ｉｊ成分とするｎ×ｍ行列をＫxyとして、一般化固有値問題ＫxyTＫxyｚ＝λＫyz、又は、正則化一般化固有値問題ＫxyTＫxy（ＫxyTＫxy＋μＫy）ＫxyTＫxyｚ＝λＫyzを解くことによって、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出する設計手段と、上記ｄ次元ベクトルｘを上記ｒ次元ベクトルｙに変換する評価手段と、とを備えている。
【００２８】
上記の構成によれば、上記カーネル主成分分析方法と同様、計算量を削減することができ、かつ、計算量の削減に伴う計算精度の低下を従来よりも小さくすることができるという効果を奏する。
【００２９】
なお、本発明に係るカーネル主成分分析装置は、コンピュータによって実現してもよい。この場合、コンピュータを上記各手段として動作させることにより、上記カーネル主成分分析装置をコンピュータにおいて実現するプログラムも本発明の範疇に入る。
【００３０】
また、本発明に係るカーネル主成分分析方法は、雑音除去方法、パターン分類方法、パターン検出方法、パターン分類方法、及び、欠損値補間方法に利用することができ、これらの方法の一部として利用されるカーネル主成分分析方法についても本発明の範疇に入る。
【発明の効果】
【００３１】
本発明によれば、計算量を削減することができ、かつ、計算量の削減に伴う計算精度の低下を従来よりも小さくすることができるという効果を奏する。
【図面の簡単な説明】
【００３２】
【図１】本発明の実施形態を示すものであり、部分カーネル主成分分析装置のブロック図である。
【図２】本発明の実施形態を示すものであり、部分カーネル主成分分析装置として機能するコンピュータのブロック図である。
【図３】本発明の実施形態を示すものであり、部分カーネル主成分分析方法のフローチャートである。
【図４】本発明の実施形態を示すものであり、部分カーネル主成分分析方法に含まれる選択ステップの構成例を示すフローチャートである。
【図５】本発明の実施形態を示すものであり、部分カーネル主成分分析方法に含まれる設計ステップの構成例を示すフローチャートである。
【図６】本発明の実施形態を示すものであり、部分カーネル主成分分析方法に含まれる評価ステップの構成例を示すフローチャートである。
【図７】人工データを用いた、本発明のカーネル主成分分析方法（部分カーネル主成分分析法）の計算精度と、従来のカーネル主成分分析方法の計算精度とを比較したグラフである。
【図８】実データを用いた、本発明のカーネル主成分分析方法（部分カーネル主成分分析法）の計算精度と、従来のカーネル主成分分析方法の計算精度とを比較したグラフである。
【図９】実データを用いた、本発明のカーネル主成分分析方法（部分カーネル主成分分析法）の計算精度と、従来のカーネル主成分分析方法の計算精度とを比較したグラフである。
【図１０】実データを用いた、本発明のカーネル主成分分析方法（部分カーネル主成分分析法）の計算精度と、従来のカーネル主成分分析方法の計算精度とを比較したグラフである。
【図１１】Ｋ−ｍｅａｎｓと無作為抽出との比較を示したグラフである。
【発明を実施するための形態】
【００３３】
本発明に係るカーネル主成分分析方法は、発明者が見出した新規なカーネル主成分分析であり、部分カーネル主成分分析（ＳＫＰＣＡ：Subset Kernel Principal Component Analysis）と呼ぶことにしたものである。そこで、本発明に係るカーネル主成分分析方法を、以下では、部分カーネル主成分分析方法と呼称する。また、本発明に係るカーネル主成分分析装置、すなわち、部分カーネル主成分分析を行う装置を、以下では、部分カーネル主成分分析装置と呼称する。
【００３４】
〔実施形態〕
本発明の一実施形態について、図面に基づいて説明すれば以下のとおりである。
【００３５】
（部分カーネル主成分分析装置）
まず、部分カーネル主成分分析装置（以下、「ＳＫＰＣＡ装置」と略記）について、図１を参照して説明する。図１は、本実施形態に係るＳＫＰＣＡ装置１の構成を示したブロック図である。ＳＫＰＣＡ装置１は、図１に示したように、選択部１０と、設計部２０と、評価部３０とを備えている。
【００３６】
ＳＫＰＣＡ装置１に入力される入力データは、標本として与えられたｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝と、評価対象として与えられたｄ次元ベクトルｘとである。以下では、標本として与えられたｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝の各々を「標本ベクトル」と呼称し、評価対象として与えられたｄ次元ベクトルｘを「評価対象ベクトル」と呼称する。一方、ＳＫＰＣＡ装置１から出力される出力データは、評価対象ベクトルｘの特徴を表現するｒ次元ベクトルｙである。以下では、このｒ次元ベクトルｙを「特徴ベクトル」と呼称する。
【００３７】
選択部１０は、入力されたｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝からｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝を選択するための手段である（ｍはｎ未満の自然数）。選択されたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝は、後述する自己カーネルグラム行列計算部２１、相互カーネルグラム行列計算部２２、及び、経験カーネル写像ベクトル計算部３１がこれを参照することができるよう、不図示の記憶装置に格納される。選択部１０により実行される選択処理の詳細については、参照する図面を代えて後述する。
【００３８】
設計部２０は、一般化固有値問題ＫxyTＫxyｚ＝λＫyzを解くことによって、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出するための手段であり、例えば、図１に示したように、自己カーネルグラム行列計算部２１と、相互カーネルグラム行列計算部２２と、一般化固有値問題計算部２３とにより構成することができる。ここで、Ｋyは、カーネル関数値ｋ（ｙi，ｙj）を第ｉｊ成分とするｍ×ｍ行列であり、「自己カーネルグラム行列」と称される。また、Ｋxyは、カーネル関数値ｋ（ｘi，ｙj）を第ｉｊ成分とするｎ×ｍ行列であり、「相互カーネルグラム行列」と称される。
【００３９】
設計部２０を構成する各部が担う機能は、以下のとおりである。
【００４０】
自己カーネルグラム行列計算部２１は、選択部１０により選択されたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝から、自己カーネルグラム行列Ｋyを計算するための手段である。より具体的に言うと、自己カーネルグラム行列計算部２１は、選択された標本ベクトルｙi及びｙjを記憶装置から読み出し、これらを予め定められたカーネル関数ｋ（ｘ、ｙ）に代入することによって、ｍ×ｍ行列である自己カーネルグラム行列Ｋyの第ｉｊ成分（Ｋy）ij＝ｋ（ｙi，ｙj）を算出する（ｉ＝１，２，…，ｍ／ｊ＝１，２，…，ｍ）。算出された自己カーネルグラム行列Ｋyは、後述する一般化固有値問題計算部２３がこれを参照することができるよう、不図示の記憶装置に格納される。
【００４１】
相互カーネルグラム行列計算部２２は、入力されたｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝と、選択部１０により選択されたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝とから、相互カーネルグラム行列Ｋxyを算出するための手段である。より具体的に言うと、相互カーネルグラム行列計算部２２は、入力された標本ベクトルｘiと選択された標本ベクトルｙjとを記憶装置から読み出し、これらをカーネル関数ｋ（ｘ，ｙ）に代入することによって、ｎ×ｍ行列である相互カーネルグラム行列Ｋxyの第ｉｊ成分（Ｋxy）ij＝ｋ（ｘi，ｙj）を算出する（ｉ＝１，２，…，ｎ／ｊ＝１，２，…，ｍ）。算出された相互カーネルグラム行列Ｋxyは、後述する一般化固有値問題計算部２３がこれを参照することができるよう、不図示の記憶装置に格納される。
【００４２】
一般化固有値問題計算部２３は、自己カーネルグラム行列計算部２１により算出された自己カーネルグラム行列Ｋyと、相互カーネルグラム行列計算部２２により算出された相互カーネルグラム行列Ｋxyとにより定義される一般化固有値問題ＫxyTＫxyｚ＝λＫyzを解くことによって、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出するための手段である。
【００４３】
より具体的に言うと、一般化固有値問題計算部２３は、（１）相互カーネルグラム行列Ｋxyの転置行列ＫxyTを算出し、（２）相互カーネルグラム行列Ｋxyとその転置行列ＫxyTとの行列積ＫxyTＫxyを算出し、（３）一般化固有値問題ＫxyTＫxyｚ＝λＫyzを解き、（４）得られた固有ベクトル（ｍ次元の列ベクトル）のうち固有値の大きいものから順にｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を選択し、（５）選択したｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝の各々を、ｚiTＫyｚiのノルムが１になるように規格化する。規格化されたｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝は、後述する特徴ベクトル計算部３２がこれを参照することができるよう、不図示の記憶装置に格納される。
【００４４】
なお、一般化固有値問題には公知の数値解法が幾つも存在しており、コンピュータ上で一般化固有値問題を解くためサブルーチンも広く利用されている。例えば、ＭＡＴＬＡＢ（登録商標）に含まれているeigsやIntel Math Kernel Library（登録商標）に含まれているdsygvxなどは、その一例である。これらのサブルーチンを利用すれば、コンピュータを一般化固有値問題計算部２３として機能させることができる。
【００４５】
評価部３０は、評価対象ベクトルｘを特徴ベクトルｙ＝ＺTｘΦに変換するための手段であり、例えば、図１に示したように、経験カーネル写像ベクトル計算部３１と、特徴ベクトル計算部３２とにより構成することができる。ここで、ｘΦは、カーネル関数値ｋ（ｘ，ｙj）を第ｊ成分とするｍ次元ベクトルであり、「経験カーネル写像ベクトル」と称される。また、Ｚは、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を並べたｍ×ｒ行列［ｚ1，ｚ2，…，ｚr］である。
【００４６】
評価部３０を構成する各部が担う機能は、以下のとおりである。
【００４７】
経験カーネル写像ベクトル計算部３１は、入力された評価対象ベクトルｘと、選択部１０により選択されたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝とから、経験カーネル写像ベクトルｘΦを算出するための手段である。より具体的に言うと、経験カーネル写像ベクトル計算部３１は、入力された評価対象ベクトルｘと選択された標本ベクトルｙjとを記憶装置から読み出し、これらをカーネル関数ｋ（ｘ、ｙ）に代入することによって、ｍ次元の列ベクトルである経験カーネル写像ベクトルｘΦの第ｊ成分（ｘΦ）j＝ｋ（ｘ，ｙj）を算出する（ｊ＝１，２，…，ｍ）。算出された経験カーネル写像ベクトルｘΦは、後述する特徴ベクトル計算部３２がこれを参照することができるよう、不図示の記憶装置に格納される。
【００４８】
特徴ベクトル計算部３２は、一般化固有値問題計算部２３により算出されたｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝と、経験カーネル写像計算部２１により算出された経験カーネル写像ベクトルｘΦとから、特徴ベクトルｙ＝ＺTｘΦを算出するための手段である。より具体的に言うと、特徴ベクトル計算部３２は、規格化された固有ベクトルｚiと経験カーネル写像ベクトルｘΦとを記憶装置から読み出し、これらの内積を算出することによって、ｒ次元ベクトルである特徴ベクトルｙの第ｉ成分を算出する（ｉ＝１，２，…，ｒ）。あるいは、行列ＺTに経験カーネル写像ベクトルｘΦを右から乗ずることによって、特徴ベクトルｙを算出する。
【００４９】
なお、算出する特徴ベクトルｙの次元ｒは、入力される評価対象ベクトルｘの次元ｄより高く設定されていてもよいし、低く設定されていてもよい。前者の場合、ノイズに強い冗長な特徴ベクトルｙを得ることができ（冗長情報付加）、後者の場合、伝送効率に優れたデータサイズの小さい特徴ベクトルｙを得ることができる（情報圧縮）。
【００５０】
また、本実施形態においては、与えられたｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝と、ｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝から選択されたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝とから算出された自己カーネルグラム行列Ｋy、相互カーネルグラム行列Ｋxy、及び、経験カーネル写像ベクトルｘΦを用いて部分カーネル主成分分析を行う構成を示したが、本発明はこれに限定されるものではない。すなわち、後述する〔部分カーネル主成分分析の原理〕からも明らかなように、ｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝と、これとは別に与えられたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝とにより定義される自己カーネルグラム行列Ｋy、相互カーネルグラム行列Ｋxy、及び、経験カーネル写像ベクトルｘΦを用いて部分カーネル主成分分析を行う構成を採用してもよい。この場合、自己カーネルグラム行列計算部２１、相互カーネルグラム行列計算部２２、及び、経験カーネル写像ベクトル計算部３１が、ｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝と同様に、ｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝を直接外部から取得するようにすればよい。
【００５１】
また、本実施形態においては、一般化固有値問題ＫxyTＫxyｚ＝λＫyzを解くことによって、部分カーネル主成分分析を行う構成を示したが、本発明はこれに限定されるものではない。すなわち、一般化固有値問題ＫxyTＫxyｚ＝λＫyzに代えて、正則化一般化固有値問題ＫxyTＫxy（ＫxyTＫxy＋μＫy）ＫxyTＫxyｚ＝λＫyzを解くことによって、部分カーネル主成分分析を行う構成を採用してもよい。なお、正則化一般化固有値問題については、後述する〔部分カーネル主成分分析の原理〕を参照されたい。
【００５２】
（カーネル主成分分析装置の構成例）
ＳＫＰＣＡ装置１は、コンピュータ（電子計算機）を用いて構成することができる。図２は、ＳＫＰＣＡ装置１として利用可能なコンピュータ１００の構成を例示したブロック図である。
【００５３】
コンピュータ１００は、図２に示したように、バス１１０を介して互いに接続された演算装置１２０と、主記憶装置１３０と、補助記憶装置１４０と、入出力インタフェース１５０とを備えている。演算装置１２０として利用可能なデバイスとしては、ＣＰＵ（Central Processing Unit）を挙げることができる。また、主記憶装置１３０として利用可能なデバイスとしては、例えば、半導体ＲＡＭ（random access memory）を挙げることができる。また、補助記憶装置１４０として利用可能なデバイスとしては、例えば、ハードディスクドライブを挙げることができる。
【００５４】
入出力インタフェース１５０には、図２に示したように、入力装置２００及び出力装置３００が接続される。ユーザにより指定された数値を入力データ（標本ベクトル及び評価対象ベクトル）とする場合には、例えば、キーボードを入力装置２００として利用することができる。また、測定値を入力データとする場合には、例えば、その測定値に応じた測定器（センサ）を入力装置２００として利用することもできる。また、出力データ（特徴ベクトル）を数値列又はグラフとしてユーザに提示する場合、例えば、モニタやプリンタなどを出力装置３００として利用することができる。
【００５５】
なお、出力装置３００を介して出力データを出力する代わりに、記録媒体に出力データを記録するようにしてもよい。また、出力データを出力装置３００に出力する代わりに、他の装置に出力データを送信するようにしてもよい。
【００５６】
補助記憶装置１４０には、コンピュータ１００をＳＫＰＣＡ装置１として動作させるためのカーネル主成分分析プログラム（以下「ＫＰＣＡプログラム」と略記する）が格納されている。ＫＰＣＡプログラムは、選択プログラムと、自己カーネルグラム行列計算プログラムと、相互カーネルグラム行列計算プログラムと、一般化固有値問題計算プログラムと、経験カーネル写像ベクトル計算プログラムと、特徴ベクトル計算プログラムとを含んでいる。
【００５７】
演算装置１２０は、主記憶装置１３０上に展開された標本プログラムに含まれる命令を実行することによって、コンピュータ１００を選択部１０として機能させる。同様に、主記憶装置１３０上に展開された自己カーネルグラム行列計算プログラム、相互カーネルグラム行列計算プログラム、一般化固有値問題計算プログラム、経験カーネル写像ベクトル計算プログラム、及び、特徴ベクトル計算プログラムに含まれる命令を実行することによって、コンピュータ１００を自己カーネルグラム行列計算部２１、相互カーネルグラム行列計算部２２、一般化固有値問題計算部２３、経験カーネル写像ベクトル計算部３１、及び、特徴ベクトル計算部３２として機能させる。
【００５８】
ＫＰＣＡプログラムを実行する過程で生成される中間データ、すなわち、標本ベクトル｛ｙ1，ｙ2，…，ｙm｝、相互カーネルグラム行列Ｋxy、自己カーネルグラム行列Ｋy、相互カーネルグラム行列の転置行列ＫxyT、相互カーネルグラム行列Ｋxyとその転置行列ＫxyTとの行列積ＫxyTＫxy、固有ベクトル｛ｚ1，ｚ2，…，ｚr｝、及び、経験カーネル写像ベクトルｘΦは、演算装置１２０がこれを参照することができるよう主記憶装置１３０に格納される。ただし、主記憶装置１３０の記憶容量に制約がある場合には、これらの中間データを補助記憶装置１４０に格納するようにしてもよい。他のデータと比べてサイズの大きい相互カーネルグラム行列Ｋxyだけを補助記憶装置１４０に格納するような構成も有効である。
【００５９】
（カーネル主成分分析方法）
次に、部分カーネル主成分分析方法（以下、「ＳＫＰＣＡ方法」と略記する）について、図３〜図６を参照して説明する。
【００６０】
図３は、本実施形態に係るＳＫＰＣＡ方法Ｓ１の流れを示したフローチャートである。ＳＫＰＣＡ方法Ｓ１は、ＳＫＰＣＡ装置１を用いたＳＫＰＣＡ方法であり、図３に示したように、選択ステップＳ１０と、設計ステップＳ２０と、評価ステップＳ３０とを含んでいる。また、図３に示したように、標本ベクトル入力ステップＳin１と、評価対象ベクトル入力ステップＳin２と、特徴ベクトル出力ステップＳoutとを更に含んでいる。
【００６１】
標本ベクトル入力ステップＳin１においては、ＳＫＰＣＡ装置１が、ｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝の入力を受け付ける。換言すれば、ｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝を入力装置等の外部装置から取得する。
【００６２】
選択ステップＳ１０においては、ＳＫＰＣＡ装置１の選択部１０が、標本ベクトル入力ステップＳin１にて入力されたｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝からｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝を選択する。
【００６３】
設計ステップＳ２０においては、ＳＫＰＣＡ装置１の設計部２０が、一般化固有値問題ＫxyTＫxyｚ＝λＫyzを解くことによって、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出する。自己カーネルグラム行列Ｋy、及び、相互カーネルグラム行列Ｋxyは、上述したとおり、標本ベクトル入力ステップＳin１にて入力されたｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝、及び、選択ステップＳ１０にて選択されたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝に基づいて定義することができる行列である。
【００６４】
評価対象ベクトル入力ステップＳin２においては、ＳＫＰＣＡ装置１が、評価対象ベクトルｘの入力を受け付ける。換言すれば、評価対象ベクトルｘを入力装置等の外部装置から取得する。なお、評価対象ベクトル入力ステップＳin２は、評価ステップＳ３０より前に実行されればよく、設計ステップＳ２０より後に実行されることを要さない。すなわち、設計ステップＳ２０より前に、あるいは、選択ステップＳ１０より前に実行されてもよい。
【００６５】
評価ステップＳ３０においては、ＳＫＰＣＡ装置１の評価部３０が、設計ステップＳ２０にて算出されたｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝とを用いて、評価対象ベクトル入力ステップＳin２にて入力された評価対象ベクトルｘを特徴ベクトルｙ＝ＺTｘΦに変換する。
【００６６】
特徴ベクトル出力ステップＳoutにおいては、ＳＫＰＣＡ装置１が、評価ステップＳ３０にて算出された特徴ベクトルｙを出力する。例えば、モニタを用いて特徴ベクトルｙを数値列又はグラフとして表示したり、プリンタを用いて特徴ベクトルｙを数値又はグラフとして印刷したりする。なお、特徴ベクトルｙを表示したり印刷したりする代わりに、特徴ベクトルｙを記録媒体に記録したり他の装置に送信したりしてもよい。
【００６７】
なお、特徴ベクトルｙをグラフとして表示または印刷する場合、特徴ベクトルｙの次元が１であれば直線上の点としてこれをプロットすることができるし、特徴ベクトルｙの次元が２であれば平面上の点としてこれをプロットすることができるし、特徴ベクトルｙの次元が３であれば、（仮想）３次元空間上の点としてこれをプロットすることができる。この場合、評価対象ベクトルｘと合わせて表示すれば、評価対象ベクトルｘと特徴ベクトルｙとの関係も一目瞭然である。
【００６８】
図４は、選択ステップＳ１０の構成例を示したフローチャートである。選択ステップＳ１０は、例えば図４に示したように、クラスタリングステップＳ１１と、標本ベクトル選択ステップＳ１２とにより構成することができる。
【００６９】
クラスタリングステップＳ１１においては、選択部１０が、標本ベクトル入力ステップＳin１にて入力されたｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝をｍ個のクラスタ｛Ｃ1，Ｃ2，…，Ｃm｝にクラスタリングする。
【００７０】
標本ベクトル選択ステップＳ１２においては、選択部１０が、クラスタリングステップＳ１１にてクラスタリングされたｍ個のクラスタ｛Ｃ1，Ｃ2，…，Ｃm｝の各々からｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝の各々を選択する。
【００７１】
なお、クラスタリングステップＳ１１においては、Φ（ｘi）が属する特徴空間Ｆにおける距離を用いたクラスタリングを行うことができる。これにより、標本ベクトルｘiが属する空間Ｒdにおける距離を用いたクラスタリングを行う場合と比べて、より精度の高い特徴ベクトルを得ることができる。特徴空間Ｆにおける距離としては、例えば、‖Φ（ｘi）−Φ（ｘj）‖2＝ｋ（ｘi，ｘi）＋ｋ（ｘj，ｘj）−２ｋ（ｘi，ｘj）が挙げられる。
【００７２】
クラスタリングステップＳ１１において利用可能な公知のクラスタリング手法としては、階層化クラスタリングやK-meansなどを挙げることができる。例えば、K-meansを利用した場合、（１）初期クラスタ重心｛＜Ｃ1＞，＜Ｃ2＞，…，＜Ｃm＞｝をｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝から適当に選び出し、（２）各標本ベクトルｘｊを、特徴空間Ｆにおける距離‖Φ（ｘｊ）−＜Ｃｉ＞‖2が最小となるクラスタＣｉに割り振り、（３）各クラスタ重心＜Ｃｉ＞を、＜Ｃi＞＝（１／｜Ｃi｜）Σｘ∈ＣiΦ（ｘ）により更新する処理を、クラスタ重心｛＜Ｃ1＞，＜Ｃ2＞，…，＜Ｃm＞｝が更新されなくなるまで繰り返すことによって、ｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝をｍ個のクラスタ｛Ｃ1，Ｃ2，…，Ｃm｝にクラスタリングすることができる。
【００７３】
また、標本ベクトル選択ステップＳ１２においては、各クラスタＣiに属する標本ベクトルのうち、特徴空間ＦにおいてクラスタＣiの重心に最も近い標本ベクトル、すなわち、argminｘ∈Ｃi‖Φ（ｘ）−＜Ｃi＞‖2を取り出し、これを標本ベクトルｙiとすることが好ましい。あるいは、クラスタＣiに属する標本ベクトルのうち、入力空間Ｒdにおいてクラスタ重心Ｃiに最も近い標本ベクトル（いわゆるPre-image）を取り出し、これを標本ベクトルｙiとすることが好ましい。これにより、各クラスタから他の標本ベクトルを取り出す場合と比べて、より精度の高い特徴ベクトルを得ることができる。なお、Pre-imageについては、以下の文献を参照されたい。
【００７４】
S. Mika, B. Scholkopf and A. Smola，“Kernel PCA and denoising in feature space”，Advances in neural information processing systems (NIPS)，11，pp，536-542 (1999)。
【００７５】
B. Scholkopf，S. Mika，C. Burges，P. Knirsch，K.-R. Muller，G. Ratsch and A. Smola，“Input space vs. feature space in kernel-based methods”，IEEE Transactions on Neural Networks，10，5，pp. 1000-1017 (1999)。
【００７６】
M. Girolami，“Mercer kernel-based clustering in feature space”，IEEE trans. on neural networks，13，3，pp. 780-784(2002)。
【００７７】
なお、選択ステップＳ１０の構成は、図４に示したものに限らない。特に、ＳＫＰＣＡ方法Ｓ１では近似誤差あるいは残差を容易に算出することができるので、前向き選択や後向き選択などの標本選択の手法を容易に導入することができる。更に、無作為抽出されたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝を用いたＳＫＰＣＡ方法Ｓ１を繰り返し、所望の精度を有する特徴ベクトルｙを算出するような構成を採ることもできる。
【００７８】
図５は、設計ステップＳ２０の構成例を示したフローチャートである。設計ステップ２０は、例えば図５に示したように、自己カーネルグラム行列計算ステップＳ２１と、相互カーネルグラム行列計算ステップＳ２２と、一般化固有値問題計算ステップＳ２３とにより構成することができる。
【００７９】
自己カーネルグラム行列計算ステップＳ２１においては、評価部２０の自己カーネルグラム行列計算部２１が、選択ステップ１０にて選択されたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝から、自己カーネルグラム行列Ｋyを算出する。
【００８０】
相互カーネルグラム行列計算ステップＳ２２においては、評価部２０の相互カーネルグラム行列計算部２２が、標本ベクトル入力ステップＳin１にて入力されたｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝と、選択ステップ１０にて選択されたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝とから、相互カーネルグラム行列Ｋxyを算出する。
【００８１】
一般化固有値問題計算ステップＳ２３においては、評価部２０の一般化固有値問題計算部２３が、自己カーネルグラム行列計算ステップＳ２１にて算出された自己カーネルグラム行列Ｋyと、相互カーネルグラム行列計算ステップＳ２２にて算出された相互カーネルグラム行列Ｋxyとにより定義される一般化固有値問題ＫxyTＫxyｚ＝λＫyzを解くことによって、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出する。
【００８２】
なお、自己カーネルグラム行列計算ステップＳ２１及び相互カーネルグラム行列計算ステップＳ２２は、一般化固有値問題計算ステップＳ２３より以前に実行されればよく、これらの実行順序は図５に示したものに限定されない。すなわち、図５においては自己カーネルグラム行列計算ステップＳ２１を実行した後で相互カーネルグラム行列計算ステップＳ２２を実行するようにしているが、相互カーネルグラム行列計算ステップＳ２２を実行した後で自己カーネルグラム行列計算ステップＳ２１を実行するようにしてもよい。
【００８３】
図６は、評価ステップＳ３０の構成例を示したフローチャートである。評価ステップＳ３０は、例えば図６に示したように、経験カーネル写像ベクトル計算ステップＳ３１と、特徴ベクトル計算ステップＳ３２と、により構成することができる。
【００８４】
経験カーネル写像ベクトル計算ステップＳ３１においては、評価部３０の経験カーネル写像ベクトル計算部３１が、評価対象ベクトル取得ステップＳin２にて入力された評価対象ベクトルｘと、選択ステップ１０にて選択されたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝とから、経験カーネル写像ベクトルｘΦを算出する。
【００８５】
特徴ベクトル計算ステップＳ３２においては、評価部３０の特徴ベクトル計算部３２が、一般化固有値問題計算ステップＳ２３にて算出されたｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝と、経験カーネル写像計算ステップＳ２１にて算出された経験カーネル写像ベクトルｘΦとから、特徴ベクトルｙを算出する。
【００８６】
〔部分カーネル主成分分析の応用〕
部分カーネル主成分分析は、カーネル主成分分析と同様、雑音除去、パターン分類、パターン検出、欠損値補完、統計的推定、統計的予測、脳コンピュータインタフェース（ＢＣＩ）などに応用することができる。これらの問題に部分カーネル主成分分析を適用する仕方は、これらの問題に従来のカーネル主成分分析を適用する仕方と同様である。
【００８７】
例えば、部分カーネル主成分分析を用いた雑音除去方法は、従来、カーネル主成分分析により得られた特徴ベクトルを用いて行われていた逆変換写像の計算を、部分カーネル主成分分析により得られた特徴ベクトルを用いて行うようにすることにより実現することができる。逆変換写像の計算に部分カーネル主成分分析を用いた場合、従来のカーネル主成分分析を用いた場合と比べて、少ない計算量で効率の良い雑音除去を行うことができる。
【００８８】
パターン分類、パターン検出、及び、欠損値補間においては、例えば、以下のようにして部分主成分分析方法を用いることができる。
【００８９】
パターン分類は、ｐ個のクラス｛Ａ1，Ａ2，…，Ａp｝に分類されているｑ個の標本ベクトル｛ｓ1，ｓ2，…，ｓq｝基づき、対象ベクトルｔを分類すべきクラスを決定する方法である。（１）各クラスＡiに分類されているｎi個の標本ベクトル｛ｓi1，ｓi2，…，ｓini｝に基づく部分カーネル主成分分析を実行することによって、対象ベクトルｔの特徴ベクトルｕiを各j=1,…,pについて算出し、（２）特徴ベクトルｕjのノルムが最大となるクラスＡjに対象ベクトルｔを分類すれば、少ない計算量で精度良くパターン分類を行うことができる。
【００９０】
パターン検出は、或るパターンを有するｑ個の標本ベクトル｛ｓ1，ｓ2，…，ｓq｝に基づき、対象ベクトルｔがそのパターンを有するか否かを判定する方法である。（１）
ｑ個の標本ベクトル｛ｓ1，ｓ2，…，ｓq｝に基づく部分カーネル主成分分析を実行することによって、対象ベクトルｔの特徴ベクトルｕを算出し、（２）特徴ベクトルｕのノルムを予め定められた閾値Ｔｈと比較することによって、対象ベクトルｔが特定のパターンを有するか否かを判定すれば、少ない計算量で精度良くパターン検出を行うことができる。
【００９１】
欠損値補完は、ｑ個の標本ベクトル｛ｓ1，ｓ2，…，ｓq｝に基づき、第ｋ成分に欠損のある対象ベクトルｔに対する欠損値補完を行う方法である。（１）第ｋ成分を異なる値に設定したｐ個の対象ベクトル｛ｔ1，ｔ2，…，ｔp｝の各々について、ｑ個の標本ベクトル｛ｘ1，ｘ2，…，ｘq｝に基づく部分カーネル主成分分析を実行することによって、対象ベクトルｔiの特徴ベクトルｕiを算出し、（２）特徴ベクトルｕjのノルムが最大となる対象ベクトルｔjの第ｋ成分の値を欠損のある対象ベクトルｔの第ｋ成分の値とすれば、少ない計算量で精度良く欠損値補間を行うことができる。
【００９２】
なお、与えられたベクトルデータを対象ベクトルｔと看做し、与えられたベクトルデータにおける未知の値を対象ベクトルｔにおける欠損成分と看做せば、対象ベクトルｔの欠損値補完を行うことにより、ベクトルデータにおける未知の値を統計的に推定することができる（統計的推定）。また、時系列データを対象ベクトルｔと看做し、時系列データにおける未来の値を対象ベクトルｔにおける欠損成分と看做せば、対象ベクトルｔの欠損値補完を行うことにより、時系列データにおける未来の値を統計的に予測することができる（統計的予測）。
【００９３】
更に、部分カーネル主成分分析は、脳コンピュータインタフェース（ＢＣＩ）に応用することができる。脳信号の測定には、脳電図（ＥＥＧ）、脳磁図（ＭＥＧ）、近赤外光トポグラフィー（ＮＩＲＳ）、ｆＭＲＩなどが用いられるが、ここでは、ＥＥＧを用いたＢＣＩについて説明する。
【００９４】
ＥＥＧによって得られた脳信号は、通常、多チャンネルであり、かつ、雑音を多く含む。そこで、ＢＣＩの前処理として、アーチファクト除去、バンドパスフィルタ（ＢＰＦ）、ＣＳＰ（Common Spatial Patterns)などの処理が施される。その後、入力信号は適当な時間幅をもつ窓で区切られ、データは数値列で表される。この数値列を評価対象ベクトルｘとして部分カーネル主成分分析を利用する。この際、標本ベクトル｛ｘ1，ｘ2，…，ｘn｝は予め与えておく。そして、部分カーネル主成分分析により得られた特徴ベクトルｙを、識別器で判別することによって、ＢＣＩの結果を得る。
【００９５】
〔部分カーネル主成分分析の原理〕
次に、本発明の基礎となる部分カーネル主成分分析の原理について説明する。
【００９６】
１．カーネル主成分分析
本節では、簡略にＫＰＣＡについて示した後、Mercerカーネル（１）を用いたＫＰＣＡの特徴付けを行う。
【００９７】
【数１】
【００９８】
ｘ1，ｘ2，…，ｘnをｄ次元の標本とする。特徴空間Ｆへ変換された標本Φ（ｘ1），…，Φ（ｘn）の相関作用素ＲFは、（２）で与えられる。
【００９９】
【数２】
【０１００】
ここで，ケット-ブラ表記｜・＞＜・｜は｜ａ＞＜ｂ｜ｃ＝＜ｃ，ｂ＞ａを満たす作用素である。有限次元では転置・Tを用いた表記ａｂTと同値であるが、Ｇａｕｓｓｉａｎカーネルなどを用いた場合は、Ｆが無限次元となるためケット-ブラ表記を用いる。ＲFの固有値分解を（３）とおく。
【０１０１】
【数３】
【０１０２】
ここで固有値λiは降順に並んでいるとする。ｒ次元の固有空間への変換ＵKと射影作用素ＰKは、（４）（５）で与えられる。
【０１０３】
【数４】
【０１０４】
【数５】
【０１０５】
ここでｅriは、Ｒrのｉ番目の標準基底（ｉ番目の要素が１でそれ以外の要素が０のベクトル）を表す。
【０１０６】
通常，ＲFは非常に大きいか、無限次元から無限次元への作用素であるため固有値分解が難しい。そこで、以下の作用素Ｓを考える。
【０１０７】
【数６】
【０１０８】
Ｓの随伴作用素をＳ*と書くと、ＲF＝（１／ｎ）ＳＳ*の固有値、固有ベクトルは、Ｋx＝Ｓ*Ｓ＝ｎΣni=1λiｖiｖiT、ｕi＝（１／√λi）Ｓｖiの関係を満たす。これより、（７）（８）が得られる。
【０１０９】
【数７】
【０１１０】
【数８】
【０１１１】
ここで、（Ｋx）ij＝ｋ（ｘi，ｘj）はカーネルグラム行列と呼ばれる。入力ｘに対して、以下の出力を得る。
【０１１２】
【数９】
【０１１３】
【数１０】
【０１１４】
ここで、‖・‖はＥｕｃｌｉｄノルムまたはｌ2ノルムを表し、Ｓ*Φ（ｘ）＝［ｋ（ｘ，ｘ1），…，ｋ（ｘ，ｘｎ）］Tは経験カーネル写像と呼ばれるｎ次元ベクトルである。
【０１１５】
以上の議論より、ＵKあるいはＰKを求めるためには、ｎ×ｎの大きさの対称行列の固有値分解が必要であり，経験カーネル写像Ｓ*Φ（ｘ）を求めるためにはｎ回のカーネル関数の計算が必要となる。さらに入力ｘの評価のために、すべての標本を保存しておく必要がある。このため、ｎが大きくなると計算量、メモリ量が問題となることがある。特に標本数が数万から数十万を超えるような大規模な問題ではＫＰＣＡを適用することができない。
【０１１６】
ＫＰＣＡを拡張するために，ＫＰＣＡの特徴付けを行う。主成分分析は、ランク制約の下で近似誤差を最小とする行列と定義できる。
【０１１７】
【数１１】
【０１１８】
ここで、Ｅxは平均を表す。ＫＰＣＡについても同様の特徴付けができる。
【０１１９】
【数１２】
【０１２０】
ここで、Ｎ（Ａ）はＡの核空間、Ｒ（Ａ）はＡの値域あるいは像を表す。制約Ｎ（Ｘ）⊃Ｒ（Ｓ）⊥は特徴空間Ｆ中のＳが張る空間以外についても考慮する必要があるため必要となる。入力次元空間の問題では、標本数が充分にあれば全空間を張ることが多いため、核空間や値域の議論を省略することが多い。しかしながら、高次元の特徴空間では、標本が全空間を張らないことが多いため、これらを正確に扱うことが重要である。本明細書においてもこれらを正確に扱う。
【０１２１】
２. 部分カーネル主成分分析
前節で示した通り、ＫＰＣＡは特徴空間Ｆで、標本Φ（ｘ1），…，Φ（ｘn）が張るＲnと同型の空間上で問題が議論される。すなわち。標本Φ（ｘ1），…，Φ（ｘn）は基底を成す。ＫＰＣＡはすべての標本を基底に用いるため、固有値分解の大きさやカーネル関数の評価回数が大きい。そこで、基底を少数の標本に置き換えることにより、この問題を解決する。
【０１２２】
ｙ1，…，ｙm（ｍ＜ｎ）を基底のための標本とし、作用素（１１）を考える。
【０１２３】
【数１３】
【０１２４】
この基底の下で全標本を最良に低階数近似する作用素を部分カーネル主成分分析（ＳＫＰＣＡ：ＳｕｂｓｅｔＫＰＣＡ）と定義する。
【０１２５】
［定義１］（ＳｕｂｓｅｔＫＰＣＡ）以下の最適化問題の解を部分ＫＰＣＡと呼ぶ。
【０１２６】
【数１４】
【０１２７】
Ｋy＝Ｔ*Ｔ、Ｋxy＝Ｓ*Ｔとおき、Ａ+で作用素または行列のＭｏｏｒｅ−Ｐｅｎｒｏｓｅ一般逆、Ａ1/2で自己共役作用素あるいは対称行列の平方根Ａ1/2ＡT/2＝Ａを表す。Ａ1/2は一般に対称ではなく、一意ではない。行列平方根やＣｈｏｌｅｓｋｙ分解などを用いることができる。次にＳＫＰＣＡの解を示す。
【０１２８】
［定理１］（ＳＫＰＣＡの解）（Ｋy1/2）+ＫxyTＫxy（ＫyT/2）+の固有値分解を（１３）とおく。
【０１２９】
【数１５】
【０１３０】
ここで、ξi (ｉ＝１，…，ｍ)は降順に並んでいるものとする。Ｗ＝［ｗ1，…，ｗr］とおくと、最適化問題（１２）の解の１つは、以下で与えられる。
【０１３１】
【数１６】
【０１３２】
特にＲ（Ｔ）⊂Ｒ（Ｓ）であるような場合、すなわち｛ｙi｝mi=1が｛ｘi｝ni=1の部分集合で与えられるような場合、（１４）で与えられる。
【０１３３】
【数１７】
【０１３４】
証明は後述する。
【０１３５】
［系１］Ｒ（Ｔ）⊂Ｒ（Ｓ）ならばＰSは正射影である。
【０１３６】
（証明）明らかにＰS*＝ＰSである。（Ｋy1/2）+Ｋy（ＫyT/2）+は、Ｒ（Ｋy）への正射影であり、Ｒ（Ｋy）⊃Ｒ（Ｗ）であることに注意すると、ＰSＰS＝Ｔ（ＫyT/2）+ＷＷT（Ｋy1/2）+Ｋy（ＫyT/2）+ＷＷT（Ｋy1/2）+Ｔ*＝ＰSであり、ＰSは正射影である(証明終)。
【０１３７】
［系２］基底にｘ1，…，ｘnを用いるとＳＫＰＣＡはＫＰＣＡに一致する。
【０１３８】
明らかであるため，証明は省略する。
【０１３９】
次にＳＫＰＣＡに別の表現を与える。
【０１４０】
［補題１］Ａ，ＢをＲ（Ａ）⊂ Ｒ（Ｂ）を満たす非負定値の行列とする。このとき、（Ｂ1/2）+Ａ（ＢT/2）+の固有値λ、と固有ベクトルｖは、任意のスカラーαに関して、Ａｕ＝λＢｕ、ｕ＝α（BT/2）+ｖを満たす。さらに、（λ，u）がＡｕ＝λＢｕを満たすとき、任意のスカラーβに対して、（Ｂ1/2）+Ａ（ＢT/2）+ｖ＝λｖ、ｖ＝βＢT/2ｕ を満たす。
【０１４１】
［命題１］補題より、Ｒ（Ｔ）⊂ Ｒ（Ｓ）ならば、一般固有値問題ＫxyTＫxyｚ＝λＫyｚの固有値の上位ｒ個に対応する固有ベクトル｛ｚi｝ri=1を用いて、（１５）と表される。
【０１４２】
【数１８】
【０１４３】
ここで、Ｚ＝［ｚ1…ｚr］であり、ｚ1，…，ｚrは，＜ｚi｜Ｋyｚi＞＝１，（ｉ＝１，…，ｒ）となるように正規化されているとする。
【０１４４】
Ｒ（Ｔ）⊂Ｒ（Ｓ）ならば、ＳＫＰＣＡによるｒ次元への変換は、ＵS＝ＺTＴ*で与えられ、入力をｘとすると（１６）である。
【０１４５】
【数１９】
【０１４６】
以上の議論より、ＳＫＰＣＡは大きさｍの（一般）固有値問題に帰着する。評価時についてもｍ回のカーネル関数の計算で評価でき，ｍ個の標本を保存しておくだけでよい。これ以降、Ｒ（Ｔ）⊂Ｒ（Ｓ）である、すなわちｙ1，…，ｙmはｘ1，…，ｘnの部分集合で与えられているとする。
【０１４７】
ＩＫＰＣＡでは、変換作用素は、異なる考え方から導き出され、（１７）で与えられている。
【０１４８】
【数２０】
【０１４９】
非特許文献２ではＩＫＰＣＡは識別のための特徴抽出に（１８）の形で用いられている。
【０１５０】
【数２１】
【０１５１】
すなわち、ＳＫＰＣＡと比較すると各次元の要素が（１／√ξi）倍されている。この影響は後段の識別性能にはあまり影響しないと考えられるが、機械学習のための特徴抽出以外の応用においては大きな誤差が生じる。
【０１５２】
３．正則化
標本数が少ないとき、統計的推定量はしばしば標本に過適合（overfit）することがある。特に高次元中で問題を解くカーネルトリックは自由度が高いため過適合しやすい。線形逆問題や回帰問題（リッジ回帰）などにおいては、過適合を防ぐために（Tikhonov）正則化がよく用いられる。そこで、本手法にも過適合を防ぐために正則化を導入する。
【０１５３】
正則化は、パラメータベクトルあるいは行列のノルムを最小化することにより実現される。ＳＫＰＣＡの正則化を正則化ＳＫＰＣＡ（ＲＳＫＰＣＡ：ＲｅｇｕｌａｒｉｚｅｄＳＫＰＣＡ）と呼び、以下に定義を示す。
【０１５４】
［定義２］（正則化ＳＫＰＣＡ）
【０１５５】
【数２２】
【０１５６】
ここで、μ＞０は正則化パラメータ、‖・‖FはＦｒｏｂｅｎｉｕｓノルムを表す。目的関数の第２項は適合の度合いを表し、正則化パラメータμが適合度合いと近似誤差のトレードオフの調整を行う。
【０１５７】
正則化ＳＫＰＣＡの解を示す。
【０１５８】
［定理２］（正則化ＳＫＰＣＡの解）ＫxyTＫxy＋μＫyが正則であるとき、正則化ＳＫＰＣＡの解の１つは、（２０）で与えられる。
【０１５９】
【数２３】
【０１６０】
ここで、ｑiは固有値分解（２１）の固有ベクトルとして与えられる。
【０１６１】
【数２４】
【０１６２】
ここで、固有値σiは降順に並んでいるものとする。ＫxyTＫxy＋μＫyが特異の場合には、逆行列は一般逆行列に置き換わる。
【０１６３】
ＳＫＰＣＡと同様に、（ＫyT/2）+ｑiは、一般固有値問題ＫxyTＫxy（ＫxyTＫxy＋μＫy）-1ＫxyTＫxyｐ＝λＫyｐで計算することができる。
【０１６４】
［命題２］Ｒ（Ｔ）⊂Ｒ（Ｓ）であるとき、ＫxyTＫxy＋μＫyが正則となる必要十分条件は、Ｎ（Ｔ）＝｛０｝である。
【０１６５】
（証明）（ｉ）必要条件を示す。ＫxyTＫxyとＫyは非負定値であるため、どちらかが正定値であれば，ＫxyTＫxy＋μＫyは正則である。Ｎ（Ｔ）＝｛０｝であれば明らかにＫyは正則である。（ｉｉ）十分条件を示す。正則であれば明らかに正定値であり、任意のｘ∈Ｒm，ｘ≠０に対して＜ｘ｜（ＫxyTＫxy＋μＫy）ｘ＞＝‖Ｓ*Ｔｘ‖2＋μ‖Ｔｘ‖2が正であるため、Ｎ（Ｔ）＝｛０｝である（証明終）。
【０１６６】
［命題３］ＧａｕｓｓｉａｎカーネルＫ（ｘ，ｙ）＝ｅｘｐ（−ｃ‖ｘ−ｙ‖2）を用いた場合、ｙ1，…，ｙmが互いに異なる場合、Ｋyは正定値である。
【０１６７】
以上の命題より、ＫxyTＫxy＋μＫyはGaussianカーネルを用いた場合は、標本が重複しなければ正則であり、また、多項式カーネルｋ（ｘ，ｙ）＝＜ｘ，ｙ＋ｃ＞dの場合においても、ｄが充分に大きければＫyは正則となる。
【０１６８】
正則化パラメータμは交差検定（cross validation）などにより実験的に決定出来る。教師なし学習として用いる場合は、学習用標本とは別に検定用標本ｚ1，…，ｚlを用意し、その検定用標本の近似誤差Σli=1‖Φ（ｚi）−ＰRΦ（ｚi）‖2を最小化するようなμを探す。検定用標本と学習用標本を入れ替え、交差検定を行うことで最適な正則化パラメータが推定できる。
【０１６９】
〔数値実験による検証〕
次に、本発明の効果を検証するための数値実験について説明する。
【０１７０】
１．人工データによる実験
本実験においては、図７の（ａ）に示す２次元の１０００個のデータを利用した。図７の（ｂ）は、Gaussianカーネルｋ（ｘ，ｙ）＝ｅｘｐ（−０．１‖ｘ−ｙ‖2）を用いたＫＰＣＡで５次元の部分空間への射影ノルムの等高線を表す。図７の（ｃ）に、標本を無作為に５０個選び出し、ＫＰＣＡを行った結果を表す。簡便のため、これ以降、（ｂ）のように標本すべてを用いたＫＰＣＡを全体ＫＰＣＡと呼び、（ｃ）のように標本の一部を使って求めたＫＰＣＡを縮減ＫＰＣＡと呼ぶ。図７の（ｄ）は縮減ＫＰＣＡと同じ標本を基底に用いたＳＫＰＣＡ、図７の（ｅ）は縮減ＫＰＣＡと同じ標本を基底に用いたＲＳＫＰＣＡ（μ＝０．１）、図７の（ｆ）はＩＫＰＣＡを示す。
【０１７１】
定量的評価には２つの手法を用いた。１つ目の評価基準は、全標本に対する近似誤差（残差）である。
【０１７２】
【数２５】
【０１７３】
Ｘは、ＰKやＰSで置き換えられる。もう１つの評価基準は全体ＫＰＣＡからの距離である。ｍ次元中のｒ次元空間部分空間への正射影行列の全体の集合はGrassmann多様体（Grassmannian）と呼ばれ、その空間における距離の１つは射影行列の差のノルムで与えられる。ＫＰＣＡとＳＫＰＣＡも特徴空間Ｆ中のｒ次元部分空間への正射影作用素であるため、差のノルムで距離を定義する。ノルムにはFrobeniusノルムを用いた。全体ＫＰＣＡの正射影作用素をＰFKとおき、距離の２乗を全体ＫＰＣＡのノルムの２乗‖ＰFK‖F2で割ったものを正規化距離Ｄとして定義する。
【０１７４】
【数２６】
【０１７５】
図より、提案手法は全体ＫＰＣＡとほとんど変わらない等高線を描くことが分かる。一方、縮減ＫＰＣＡとＩＫＰＣＡの等高線はＫＰＣＡから大きくかけ離れたものとなっている。近似誤差Ｅと正規化距離Ｄにおいても、提案手法は全体ＫＰＣＡとほぼ変わらないのに対して、縮減ＫＰＣＡ、ＩＫＰＣＡはどちらも大きな値となっている。
【０１７６】
Ｉｎｔｅｌ Ｃｏｒｅ２ Ｄｕｏ（登録商標） ２．１６ＧＨｚを搭載したＭａｃＢｏｏｋ（登録商標）でＭＡＴＬＡＢ（登録商標）のＩＲＡＭを利用した固有値分解の関数ｅｉｇｓでは、１０００×１０００の対称行列の固有値分解で０．４秒、５０×５０の対称行列の固有値分解で０．０２秒であった。この実験の場合では、提案法と縮減ＫＰＣＡは、全体ＫＰＣＡの場合より２０倍高速に固有値分解を行うことができる。さらにカーネル関数の評価回数も２０ 分の１の回数で計算ができる。
【０１７７】
２．実データを用いた実験
UCI Machine learning repositoryにあるconcreteとhousingのデータセットを用いて実験を行った。concreteは９次元、１０３０個の標本があり、housingは１４次元、５０６個の標本がある。カーネル関数にはGaussianカーネルを用い、パラメータｃには、全要素の分散σ から、ｃ＝１／２σ2を求めて用いた。階数は入力次元と同じ次元で実験を行った。
【０１７８】
図８の（ａ）及び（ｂ）は全体ＫＰＣＡからの正規化距離を示す。ＩＫＰＣＡは、比較できない程度の大きな距離を示したため省略した。標本の部分集合は無作為抽出で行い、１００回の無作為抽出の平均、標準偏差、最小値を示す。図より、標本数が全体の２０％から１０％程度であっても提案手法は全体ＫＰＣＡから１％程度の誤差しか生じないことが確認できる。図９の（ａ）及び（ｂ）は式(２２)の誤差を全体ＫＰＣＡの誤差で割った正規化誤差を示す。標本数が２０％程度に減少した場合でも、正規化誤差の最小値は、全体ＫＰＣＡの誤差の１．０１倍程度であることが分かる。図１０の（ａ）及び（ｂ）は設計時の計算時間を示す。数値計算はIntel Math Kernel Library（登録商標）を用い、全体ＫＰＣＡ、縮減ＫＰＣＡの固有値分解にはdsyevx、ＳＫＰＣＡの一般固有値分解には、dsygvxを用いた。図より、基底に用いる標本数が全標本数の７０％以下の場合、提案手法の方が高速に設計を行うことができることがわかる。
【０１７９】
図１１はhousingにおける無作為抽出と特徴空間中Ｋ‐ｍｅａｎｓによる標本選択の比較である。Ｋ‐ｍｅａｎｓの方が平均的によい近似性能を示し、標準偏差も小さいことが確認できる。
【０１８０】
最後に、正則化の効果を示すためにhousingの部分集合として１００個（ｎ＝１００）、基底にｍ個の標本を用いて、この部分集合に対する近似を最小にするＳＫＰＣＡ及び、正則化ＳＫＰＣＡを求め、それぞれの全体ＫＰＣＡからの正規化距離を比較した。１００回の試行における正規化距離の平均値を表１に示す。実験結果より、正則化を用いた方がよい近似を与えることが確認できた。
【０１８１】
【表１】
【０１８２】
〔定理の証明〕
［補題２］制約条件Ｎ（Ｘ）⊃Ｒ（Ｔ）⊥は、ある作用素Ａ：Ｒm→Ｆを用いて、Ｘ＝ＡＴ*と表されることと同値である。
【０１８３】
【数２７】
【０１８４】
（証明）（ｉ）<＝はＲ（Ｔ）⊥＝Ｎ（Ｔ）⊂Ｎ（Ｘ）より明らかである。（ｉｉ）⇒の方向、すなわち、「Ｎ（Ｘ）⊃Ｒ（Ｔ）⊥のときＸ＝ＡＴ*となるＡが存在する。」を示す。Ａ＝Ｘ（Ｔ*）+とおけばＡＴ*＝Ｘ（Ｔ*）+Ｔ*＝ＸＰR(T)＝Ｘ（Ｉ−ＰR(T)⊥）である。ここで、Ｐχはχへの正射影を表す。Ｎ（Ｘ）⊃Ｒ（Ｔ）⊥より、ＸＰR(T)⊥＝０であり、ＡＴ*＝Ｘとなる。（ｉ），（ｉｉ）より補題を得る（証明終）。
【０１８５】
［補題３］制約条件Ｒ（Ｘ）⊂Ｒ（Ｔ）⊥は、ある作用素Ａ：Ｆ→Ｒmを用いて、Ｘ＝ＴＡと表されることと同値である。
【０１８６】
証明は、補題２と同様なので省略する。
【０１８７】
定理２でμ＝０とおけば定理１となるため、μ≧０として定理２のみを証明する。補題２，３よりＸ＝ＴＡＴ*とおける。作用素Ａに対するトレースを任意の正規直交基底ｕi，ｉ＝１，…を用いてＴｒａｃｅ（Ａ）＝Σi＜ｕi｜Ａｕi＞と定義する。最適化問題（１９）の目的関数をＪとおくと以下のとおりである。
【０１８８】
【数２８】
【０１８９】
Ｒ（ＫxyT）＝Ｒ（Ｔ*Ｓ）⊂Ｒ（Ｔ*）＝Ｒ（Ｋy）より、ＫxyT＝ＰR(Ky)ＫxyT＝ＫyＫy+ＫxyTであり、Ｋxy＝Ｋxy（ＫyT/2）+ＫyT/2である。また、Ｒ（ＫxyTＫxy＋μＫy）⊃Ｒ（ＫxyT）より、ＰR(KxyTKxy+μKy)ＫxyT＝（ＫxyTＫxy＋μＫy）1/2（（ＫxyTＫxy＋μＫy）1/2）+ＫxyT＝ＫxyTである。これより、以下のようになる。
【０１９０】
【数２９】
【０１９１】
ここでＣはＡに対する定数項を表す。式（２１） の固有値分解より、特異値分解を以下のように置く。
【０１９２】
【数３０】
【０１９３】
Schmidtの近似定理（Eckart−Youngの定理）より、以下を満たすＡが存在するときは、上式を満たすＡがＪを最小化する。
【０１９４】
【数３１】
【０１９５】
【数３２】
【０１９６】
以上より、式（Ａ．１）を満たすＡが存在し、その解の一つは、以下で与えられる。
【０１９７】
【数３３】
【０１９８】
証明終
【産業上の利用可能性】
【０１９９】
本発明は、カーネル主成分分析を用いたデータ解析一般に広く適用することができる。特徴抽出、雑音除去、パターン分類、パターン検出、欠損値補完、統計的推定、統計的予測などに好適に利用することができる。
【符号の説明】
【０２００】
１ ＳＫＰＣＡ装置（カーネル主成分分析装置）
１０ 選択部（選択手段）
２０ 設計部（設計手段）
２１ 自己カーネルグラム行列計算部
２２ 相互カーネルグラム行列計算部
２３ 一般化固有値問題計算部
３０ 評価部（評価手段）
３１ 経験カーネル写像ベクトル計算部
３２ 特徴ベクトル計算部
【技術分野】
【０００１】
本発明は、主成分分析に関するものであり、特に、カーネル主成分分析に関するものである。
【背景技術】
【０００２】
多変量解析のひとつとして、主成分分析（ＰＣＡ：principal component analysis）がある。主成分分析は、ｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝（ｄ次元ベクトル）に基づき、ｄ次元の評価対象ベクトルｘの特徴を表すｒ次元の特徴ベクトルｙを算出する方法であり、（１）ｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝から定義される分散共分散行列の固有値問題を解き、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出する設計ステップと、（２）ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を用いて評価対象ベクトルｘを特徴ベクトルｙに変換する評価ステップとにより構成することができる。得られる特徴ベクトルｙの各成分は、評価対象ベクトルｘの各主成分への射影を表し、「主成分得点」と称される。
【０００３】
この主成分分析の拡張として、カーネル主成分分析（ＫＰＣＡ：Kernel principal component analysis）がある（非特許文献１参照）。カーネル主成分分析は、特徴空間と呼ばれる再生核ヒルベルト空間Ｆ上での主成分分析であり、例えば（１）ｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝から定義されるカーネルグラム行列Ｋｘの固有値値問題を解き、ｒ個の固有ベクトル固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出する設計ステップと、（２）ｒ個の固有ベクトル固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を用いて評価対象ベクトルｘを特徴ベクトルｙに変換する評価ステップとにより構成することができる。ここで、カーネルグラム行列Ｋｘとは、ｋ（ｘ、ｙ）＝＜Φ（ｘ）｜Φ（ｙ）＞をＭｅｒｃｅｒのカーネル関数として、カーネル関数値ｋ（ｘi，ｘj）を第ｉｊ成分とするｎ×ｎ行列のことである。Φ（ｘ）は入力空間Ｒdから特徴空間Ｆへの写像を表し、＜・｜・＞は内積を表す。
【０００４】
カーネル主成分分析を用いることにより、非線形相関を考慮した多変量解析を行うことができる。しかしながら、ｎ×ｎ行列の固有値問題を解く必要があるため、計算量（計算の複雑さ）はｎ3に比例し、必要なメモリサイズはｎ2に比例する。したがって、現在のコンピュータの性能では、標本数ｎが数万から数十万を超えると、カーネル主成分分析を行うことが困難になる。このため、与えられたｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝の一部のみを用いたカーネル主成分分析しか行えないことがしばしばであった。
【０００５】
このような問題を考慮したカーネル主成分分析の改良としては、改良カーネル主成分分析（ＩＫＰＣＡ：Improved Kernel principal component analysis）（非特許文献２参照）や、スパースカーネル主成分分析（非特許文献３、４）などが知られている。
【先行技術文献】
【非特許文献】
【０００６】
【非特許文献１】B. Scholkopf and A. Smola and K.-R. Muller，“Nonlinear Component Analysis as a Kernel Eigenvalue Problem”，Neural computation，10，5，pp. 1299-1319，1998
【非特許文献２】Y. Xu, D. Zhang, F. Song, J. Yang, Z. Jing and M. Li，“Amethod for speeding up feature extraction based on KPCA”，Neurocomputing， pp. 1056-1061，2007
【非特許文献３】M. E. Tipping，“Sparse kernel principal component analysis”，Advances in Neural Information Processing Systems(NIPS)，13，pp. 633-639，2001
【非特許文献４】A. J. Smola, O. L. Mgngasarian and B. Scholkopf，“Sparse kernel feature analysis”，Technical report 99-04，Universityof Wisconsin，1999
【発明の概要】
【発明が解決しようとする課題】
【０００７】
しかしながら、上記従来の技術には、計算量を削減すると計算精度が大幅に低下するという問題があった。
【０００８】
すなわち、与えられたｎ個の標本ベクトルの一部のみを用いてカーネル主成分分析を行えば、計算量も標本数の減少に応じて当然に減少する。しかし、用いる標本ベクトルの数を減らすと、得られる特徴ベクトルの精度低下を招来する。例えば、１０００個の標本ベクトルが与えられたときに、５０個の標本ベクトルを用いたカーネル主成分分析を行うと、得られる特徴ベクトルの精度が大幅に低下する。
【０００９】
また、改良カーネル主成分分析は、適切な問題設定ではないため、近似が正確でないという問題を有している。また、スパースカーネル主成分分析は、評価ステップにおける計算量を削減することを目的とするものであり、設計ステップにおける計算量を削減することができないという問題を有している。
【００１０】
本発明は、上記の問題に鑑みてなされたものであり、その目的は、計算量を削減することができ、かつ、計算量の削減に伴う計算精度の低下が従来よりも小さなカーネル主成分分析方法を実現することにある。
【課題を解決するための手段】
【００１１】
本発明に係るカーネル主成分分析方法は、カーネル主成分分析装置を用いて、標本として与えられたｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝に基づき、評価対象として与えられたｄ次元ベクトルｘの特徴を表現するｒ次元ベクトルｙを算出するカーネル主成分分析方法である。上記課題を解決するために、本発明に係るカーネル主成分分析方法は、上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝からｍ個（ｍはｎ未満の自然数）のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を選択する選択ステップ、又は、上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝とは別にｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を取得する取得ステップと、カーネル関数値ｋ（ｙi，ｙj）を第ｉｊ成分とするｍ×ｍ行列をＫyとし、カーネル関数値ｋ（ｘi，ｙj）を第ｉｊ成分とするｎ×ｍ行列をＫxyとして、固有値λと固有ベクトルｚを求める一般化固有値問題ＫxyTＫxyｚ＝λＫyz、又は、正則化一般化固有値問題ＫxyTＫxy（ＫxyTＫxy＋μＫy）ＫxyTＫxyｚ＝λＫyzを解くことによって、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出する設計ステップと、上記ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を用いて、上記ｄ次元ベクトルｘを上記ｒ次元ベクトルｙに変換する評価ステップと、を含んでいる。ここで、ｋ（ｘ、ｙ）＝＜Φ（ｘ）｜Φ（ｙ）＞は、従来のカーネル主成分分析に用いられているものと同じMercerのカーネル関数であり、μ＞０は正則化パラメータである。
【００１２】
上記の構成によれば、解くべき固有値問題はｍ×ｍ行列（ＫxyTＫxy及びＫy）に関する一般化固有値問題になる（計算量ｍ3）。このため、ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝を用いた従来のカーネル主成分分析方法（計算量ｎ3）と比べて計算量を減少させることができる。しかも、本発明において解くべき一般化固有値問題には、（選択された）ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝に関する情報のみならず、（選択されたなかった）ｎ−ｍ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝＼｛ｙ1，ｙ2，…，ｙm｝に関する情報がｎ×ｍ行列Ｋxyを介して取り込まれている。したがって、（選択された）ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝のみを用いて従来のカーネル主成分分析を行う場合に比べて、より高い精度で特徴ベクトルｙを算出することができる。したがって、計算量を削減することができ、かつ、計算量の削減に伴う計算精度の低下を従来よりも小さくすることができるという効果を奏する。
【００１３】
なお、ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を用いて従来のカーネル主成分分析を行う場合に比べて、より高い精度の特徴ベクトルｙを算出することができることの、数学的な裏付けについては発明の詳細な説明の〔部分カーネル主成分分析の原理〕を、実験的な検証については発明の詳細な説明の〔数値実験による検証〕をそれぞれ参照されたい。
【００１４】
なお、上記設計ステップは、例えば、（１）上記ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝から、上記ｍ×ｍ行列Ｋyを算出するステップと、（２）上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝、及び、上記ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝から、上記ｎ×ｍ行列Ｋxyを算出するステップと、（３）上記一般化固有値問題ＫxyTＫxyｚ＝λＫyz、又は、ＫxyTＫxy（ＫxyTＫxy＋μＫy）ＫxyTＫxyｚ＝λＫyzを解くことによって、上記ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出するステップとを含んで構成することができる。
【００１５】
ただし、上記設計ステップは、カーネル関数値ｋ（ｙi，ｙj）を第ｉｊ成分とするｍ×ｍ行列をＫyとし、カーネル関数値ｋ（ｘi，ｙj）を第ｉｊ成分とするｎ×ｍ行列をＫxyとしたときに、ＫxyTＫxyｚ＝λＫyzにより表現される一般化固有値問題、又は、ＫxyTＫxy（ＫxyTＫxy＋μＫy）ＫxyTＫxyｚ＝λＫyzにより表現される正則化一般化固有値問題を解くことによって、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出するものであれば何でもよく、必ずしも、上記ｍ×ｍ行列Ｋy自体、及び、上記ｎ×ｍ行列Ｋxy自体の算出を要するものではない。
【００１６】
また、上記評価ステップは、例えば、（１）上記ｄ次元ベクトルｘ、及び、上記ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝から、カーネル関数値ｋ（ｘ，ｙj）を第ｊ成分とするｍ次元ベクトルｘΦを算出するステップと、（２）上記ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝、及び、上記ｍ次元ベクトルｘΦから、上記ｒ次元ベクトルｙ＝［ｚ1，ｚ2，…，ｚr］TｘΦを算出するステップと、を含んで構成することができる。
【００１７】
ただし、上記評価ステップは、カーネル関数値ｋ（ｘ，ｙj）を第ｊ成分とするｍ次元ベクトルをｘΦとし、上記ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を並べたｍ×ｒ行列［ｚ1，ｚ2，…，ｚr］をＺとしたときに、ＺTｘΦにより表現されるｒ次元ベクトルを算出するものであれば何でもよく、必ずしも、上記ｍ次元ベクトルをｘΦ自体、及び、上記ｍ×ｒ行列Ｚ自体の算出を要するものではない。
【００１８】
また、本発明に係るカーネル主成分分析方法は、上記評価ステップにて算出された上記ｒ次元ベクトルｙを他の装置に出力する出力ステップを含んでいてもよい。出力ステップにおける出力先となる装置は特に限定されるものではないが、例えば、上記ｒ次元ベクトルｙを表示装置（ディスプレイ等）、印刷装置（プリンタ等）、記録装置（光ディスクドライブ等）、送信装置（モデム等）を出力先とすることができる。表示装置／印刷装置を出力先とする場合には、上記ｒ次元ベクトルｙを数値列又はグラフとして表示又は印刷することができるし、記録装置を出力先とする場合には、上記ｒ次元ベクトルｙを記録媒体に記録することができるし、送信装置を出力先とする場合には、送信装置を介して上記ｒ次元ベクトルｙを他のコンピュータに送信することができる。これにより、ユーザが上記ｒ次元ベクトルｙを視覚的に認識したり、上記カーネル主成分分析装置または上記記録媒体から情報を読み出すことができる他装置が上記記録媒体から上記ｒ次元ベクトルｙを読み出して利用したり、他のコンピュータが上記送信装置から上記ｒ次元ベクトルｙを受信して利用したりすることができる。
【００１９】
本発明に係るカーネル主成分分析方法において、上記選択ステップは、上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝をｍ個のクラスタ｛Ｃ1，Ｃ2，…，Ｃm｝にクラスタリングするステップと、上記ｍ個のクラスタ｛Ｃ1，Ｃ2，…，Ｃm｝の各々から上記ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝の各々を選択するステップと、を含んでいる、ことが好ましい。
【００２０】
上記構成によれば、ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を選択するために要する計算量を徒に増加させることなく、特徴ベクトルｙを精度良く計算することができるよう、ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を選択することができる、という更なる効果を奏する。
【００２１】
本発明に係るカーネル主成分分析方法は、上記選択ステップにおいて、上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝をｍ個のクラスタ｛Ｃ1，Ｃ2，…，Ｃm｝にクラスタリングする際に、特徴空間上のノルムを使用する、ことが好ましい。
【００２２】
上記の構成によれば、特徴ベクトルｙを更に精度良く計算することができるよう、ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を選択することができる、という更なる効果を奏す。
【００２３】
本発明に係るカーネル主成分分析方法において、上記ベクトルｙの次元ｒは、上記ベクトルｘの次元ｄよりも小さい、ことが好ましい。
【００２４】
上記の構成によれば、伝送効率に優れたデータサイズの小さい特徴ベクトルｙを得ることができる。例えば、本発明を情報源符号化に適用する場合、送信側（装置）と受信側（装置）とが上記ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を予め共有することによって、送受信する情報量を削減することができる。
【００２５】
本発明に係るカーネル主成分分析方法において、上記ベクトルｙの次元ｒは、上記ベクトルｘの次元ｄよりも大きい、ことが好ましい。
【００２６】
上記の構成によれば、ノイズに強い冗長な特徴ベクトルｙを得ることができる。例えば、本発明を通信路符号化に適用する場合、送信側（装置）と受信側（装置）とが上記ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を予め共有することによって、雑音に対して頑健な情報を送受信することができる。
【００２７】
本発明に係るカーネル主成分分析装置は、標本として与えられたｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝に基づき、評価対象として与えられたｄ次元ベクトルｘの特徴を表現するｒ次元ベクトルｙを算出するカーネル主成分分析装置である。そして、上記課題を解決するために、上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝からｍ個（ｍはｎ未満の自然数）のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を選択する選択手段、又は、上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝とは別にｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を取得する取得手段と、カーネル関数値ｋ（ｙi，ｙj）を第ｉｊ成分とするｍ×ｍ行列をＫyとし、カーネル関数値ｋ（ｘi，ｙj）を第ｉｊ成分とするｎ×ｍ行列をＫxyとして、一般化固有値問題ＫxyTＫxyｚ＝λＫyz、又は、正則化一般化固有値問題ＫxyTＫxy（ＫxyTＫxy＋μＫy）ＫxyTＫxyｚ＝λＫyzを解くことによって、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出する設計手段と、上記ｄ次元ベクトルｘを上記ｒ次元ベクトルｙに変換する評価手段と、とを備えている。
【００２８】
上記の構成によれば、上記カーネル主成分分析方法と同様、計算量を削減することができ、かつ、計算量の削減に伴う計算精度の低下を従来よりも小さくすることができるという効果を奏する。
【００２９】
なお、本発明に係るカーネル主成分分析装置は、コンピュータによって実現してもよい。この場合、コンピュータを上記各手段として動作させることにより、上記カーネル主成分分析装置をコンピュータにおいて実現するプログラムも本発明の範疇に入る。
【００３０】
また、本発明に係るカーネル主成分分析方法は、雑音除去方法、パターン分類方法、パターン検出方法、パターン分類方法、及び、欠損値補間方法に利用することができ、これらの方法の一部として利用されるカーネル主成分分析方法についても本発明の範疇に入る。
【発明の効果】
【００３１】
本発明によれば、計算量を削減することができ、かつ、計算量の削減に伴う計算精度の低下を従来よりも小さくすることができるという効果を奏する。
【図面の簡単な説明】
【００３２】
【図１】本発明の実施形態を示すものであり、部分カーネル主成分分析装置のブロック図である。
【図２】本発明の実施形態を示すものであり、部分カーネル主成分分析装置として機能するコンピュータのブロック図である。
【図３】本発明の実施形態を示すものであり、部分カーネル主成分分析方法のフローチャートである。
【図４】本発明の実施形態を示すものであり、部分カーネル主成分分析方法に含まれる選択ステップの構成例を示すフローチャートである。
【図５】本発明の実施形態を示すものであり、部分カーネル主成分分析方法に含まれる設計ステップの構成例を示すフローチャートである。
【図６】本発明の実施形態を示すものであり、部分カーネル主成分分析方法に含まれる評価ステップの構成例を示すフローチャートである。
【図７】人工データを用いた、本発明のカーネル主成分分析方法（部分カーネル主成分分析法）の計算精度と、従来のカーネル主成分分析方法の計算精度とを比較したグラフである。
【図８】実データを用いた、本発明のカーネル主成分分析方法（部分カーネル主成分分析法）の計算精度と、従来のカーネル主成分分析方法の計算精度とを比較したグラフである。
【図９】実データを用いた、本発明のカーネル主成分分析方法（部分カーネル主成分分析法）の計算精度と、従来のカーネル主成分分析方法の計算精度とを比較したグラフである。
【図１０】実データを用いた、本発明のカーネル主成分分析方法（部分カーネル主成分分析法）の計算精度と、従来のカーネル主成分分析方法の計算精度とを比較したグラフである。
【図１１】Ｋ−ｍｅａｎｓと無作為抽出との比較を示したグラフである。
【発明を実施するための形態】
【００３３】
本発明に係るカーネル主成分分析方法は、発明者が見出した新規なカーネル主成分分析であり、部分カーネル主成分分析（ＳＫＰＣＡ：Subset Kernel Principal Component Analysis）と呼ぶことにしたものである。そこで、本発明に係るカーネル主成分分析方法を、以下では、部分カーネル主成分分析方法と呼称する。また、本発明に係るカーネル主成分分析装置、すなわち、部分カーネル主成分分析を行う装置を、以下では、部分カーネル主成分分析装置と呼称する。
【００３４】
〔実施形態〕
本発明の一実施形態について、図面に基づいて説明すれば以下のとおりである。
【００３５】
（部分カーネル主成分分析装置）
まず、部分カーネル主成分分析装置（以下、「ＳＫＰＣＡ装置」と略記）について、図１を参照して説明する。図１は、本実施形態に係るＳＫＰＣＡ装置１の構成を示したブロック図である。ＳＫＰＣＡ装置１は、図１に示したように、選択部１０と、設計部２０と、評価部３０とを備えている。
【００３６】
ＳＫＰＣＡ装置１に入力される入力データは、標本として与えられたｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝と、評価対象として与えられたｄ次元ベクトルｘとである。以下では、標本として与えられたｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝の各々を「標本ベクトル」と呼称し、評価対象として与えられたｄ次元ベクトルｘを「評価対象ベクトル」と呼称する。一方、ＳＫＰＣＡ装置１から出力される出力データは、評価対象ベクトルｘの特徴を表現するｒ次元ベクトルｙである。以下では、このｒ次元ベクトルｙを「特徴ベクトル」と呼称する。
【００３７】
選択部１０は、入力されたｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝からｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝を選択するための手段である（ｍはｎ未満の自然数）。選択されたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝は、後述する自己カーネルグラム行列計算部２１、相互カーネルグラム行列計算部２２、及び、経験カーネル写像ベクトル計算部３１がこれを参照することができるよう、不図示の記憶装置に格納される。選択部１０により実行される選択処理の詳細については、参照する図面を代えて後述する。
【００３８】
設計部２０は、一般化固有値問題ＫxyTＫxyｚ＝λＫyzを解くことによって、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出するための手段であり、例えば、図１に示したように、自己カーネルグラム行列計算部２１と、相互カーネルグラム行列計算部２２と、一般化固有値問題計算部２３とにより構成することができる。ここで、Ｋyは、カーネル関数値ｋ（ｙi，ｙj）を第ｉｊ成分とするｍ×ｍ行列であり、「自己カーネルグラム行列」と称される。また、Ｋxyは、カーネル関数値ｋ（ｘi，ｙj）を第ｉｊ成分とするｎ×ｍ行列であり、「相互カーネルグラム行列」と称される。
【００３９】
設計部２０を構成する各部が担う機能は、以下のとおりである。
【００４０】
自己カーネルグラム行列計算部２１は、選択部１０により選択されたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝から、自己カーネルグラム行列Ｋyを計算するための手段である。より具体的に言うと、自己カーネルグラム行列計算部２１は、選択された標本ベクトルｙi及びｙjを記憶装置から読み出し、これらを予め定められたカーネル関数ｋ（ｘ、ｙ）に代入することによって、ｍ×ｍ行列である自己カーネルグラム行列Ｋyの第ｉｊ成分（Ｋy）ij＝ｋ（ｙi，ｙj）を算出する（ｉ＝１，２，…，ｍ／ｊ＝１，２，…，ｍ）。算出された自己カーネルグラム行列Ｋyは、後述する一般化固有値問題計算部２３がこれを参照することができるよう、不図示の記憶装置に格納される。
【００４１】
相互カーネルグラム行列計算部２２は、入力されたｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝と、選択部１０により選択されたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝とから、相互カーネルグラム行列Ｋxyを算出するための手段である。より具体的に言うと、相互カーネルグラム行列計算部２２は、入力された標本ベクトルｘiと選択された標本ベクトルｙjとを記憶装置から読み出し、これらをカーネル関数ｋ（ｘ，ｙ）に代入することによって、ｎ×ｍ行列である相互カーネルグラム行列Ｋxyの第ｉｊ成分（Ｋxy）ij＝ｋ（ｘi，ｙj）を算出する（ｉ＝１，２，…，ｎ／ｊ＝１，２，…，ｍ）。算出された相互カーネルグラム行列Ｋxyは、後述する一般化固有値問題計算部２３がこれを参照することができるよう、不図示の記憶装置に格納される。
【００４２】
一般化固有値問題計算部２３は、自己カーネルグラム行列計算部２１により算出された自己カーネルグラム行列Ｋyと、相互カーネルグラム行列計算部２２により算出された相互カーネルグラム行列Ｋxyとにより定義される一般化固有値問題ＫxyTＫxyｚ＝λＫyzを解くことによって、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出するための手段である。
【００４３】
より具体的に言うと、一般化固有値問題計算部２３は、（１）相互カーネルグラム行列Ｋxyの転置行列ＫxyTを算出し、（２）相互カーネルグラム行列Ｋxyとその転置行列ＫxyTとの行列積ＫxyTＫxyを算出し、（３）一般化固有値問題ＫxyTＫxyｚ＝λＫyzを解き、（４）得られた固有ベクトル（ｍ次元の列ベクトル）のうち固有値の大きいものから順にｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を選択し、（５）選択したｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝の各々を、ｚiTＫyｚiのノルムが１になるように規格化する。規格化されたｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝は、後述する特徴ベクトル計算部３２がこれを参照することができるよう、不図示の記憶装置に格納される。
【００４４】
なお、一般化固有値問題には公知の数値解法が幾つも存在しており、コンピュータ上で一般化固有値問題を解くためサブルーチンも広く利用されている。例えば、ＭＡＴＬＡＢ（登録商標）に含まれているeigsやIntel Math Kernel Library（登録商標）に含まれているdsygvxなどは、その一例である。これらのサブルーチンを利用すれば、コンピュータを一般化固有値問題計算部２３として機能させることができる。
【００４５】
評価部３０は、評価対象ベクトルｘを特徴ベクトルｙ＝ＺTｘΦに変換するための手段であり、例えば、図１に示したように、経験カーネル写像ベクトル計算部３１と、特徴ベクトル計算部３２とにより構成することができる。ここで、ｘΦは、カーネル関数値ｋ（ｘ，ｙj）を第ｊ成分とするｍ次元ベクトルであり、「経験カーネル写像ベクトル」と称される。また、Ｚは、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を並べたｍ×ｒ行列［ｚ1，ｚ2，…，ｚr］である。
【００４６】
評価部３０を構成する各部が担う機能は、以下のとおりである。
【００４７】
経験カーネル写像ベクトル計算部３１は、入力された評価対象ベクトルｘと、選択部１０により選択されたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝とから、経験カーネル写像ベクトルｘΦを算出するための手段である。より具体的に言うと、経験カーネル写像ベクトル計算部３１は、入力された評価対象ベクトルｘと選択された標本ベクトルｙjとを記憶装置から読み出し、これらをカーネル関数ｋ（ｘ、ｙ）に代入することによって、ｍ次元の列ベクトルである経験カーネル写像ベクトルｘΦの第ｊ成分（ｘΦ）j＝ｋ（ｘ，ｙj）を算出する（ｊ＝１，２，…，ｍ）。算出された経験カーネル写像ベクトルｘΦは、後述する特徴ベクトル計算部３２がこれを参照することができるよう、不図示の記憶装置に格納される。
【００４８】
特徴ベクトル計算部３２は、一般化固有値問題計算部２３により算出されたｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝と、経験カーネル写像計算部２１により算出された経験カーネル写像ベクトルｘΦとから、特徴ベクトルｙ＝ＺTｘΦを算出するための手段である。より具体的に言うと、特徴ベクトル計算部３２は、規格化された固有ベクトルｚiと経験カーネル写像ベクトルｘΦとを記憶装置から読み出し、これらの内積を算出することによって、ｒ次元ベクトルである特徴ベクトルｙの第ｉ成分を算出する（ｉ＝１，２，…，ｒ）。あるいは、行列ＺTに経験カーネル写像ベクトルｘΦを右から乗ずることによって、特徴ベクトルｙを算出する。
【００４９】
なお、算出する特徴ベクトルｙの次元ｒは、入力される評価対象ベクトルｘの次元ｄより高く設定されていてもよいし、低く設定されていてもよい。前者の場合、ノイズに強い冗長な特徴ベクトルｙを得ることができ（冗長情報付加）、後者の場合、伝送効率に優れたデータサイズの小さい特徴ベクトルｙを得ることができる（情報圧縮）。
【００５０】
また、本実施形態においては、与えられたｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝と、ｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝から選択されたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝とから算出された自己カーネルグラム行列Ｋy、相互カーネルグラム行列Ｋxy、及び、経験カーネル写像ベクトルｘΦを用いて部分カーネル主成分分析を行う構成を示したが、本発明はこれに限定されるものではない。すなわち、後述する〔部分カーネル主成分分析の原理〕からも明らかなように、ｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝と、これとは別に与えられたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝とにより定義される自己カーネルグラム行列Ｋy、相互カーネルグラム行列Ｋxy、及び、経験カーネル写像ベクトルｘΦを用いて部分カーネル主成分分析を行う構成を採用してもよい。この場合、自己カーネルグラム行列計算部２１、相互カーネルグラム行列計算部２２、及び、経験カーネル写像ベクトル計算部３１が、ｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝と同様に、ｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝を直接外部から取得するようにすればよい。
【００５１】
また、本実施形態においては、一般化固有値問題ＫxyTＫxyｚ＝λＫyzを解くことによって、部分カーネル主成分分析を行う構成を示したが、本発明はこれに限定されるものではない。すなわち、一般化固有値問題ＫxyTＫxyｚ＝λＫyzに代えて、正則化一般化固有値問題ＫxyTＫxy（ＫxyTＫxy＋μＫy）ＫxyTＫxyｚ＝λＫyzを解くことによって、部分カーネル主成分分析を行う構成を採用してもよい。なお、正則化一般化固有値問題については、後述する〔部分カーネル主成分分析の原理〕を参照されたい。
【００５２】
（カーネル主成分分析装置の構成例）
ＳＫＰＣＡ装置１は、コンピュータ（電子計算機）を用いて構成することができる。図２は、ＳＫＰＣＡ装置１として利用可能なコンピュータ１００の構成を例示したブロック図である。
【００５３】
コンピュータ１００は、図２に示したように、バス１１０を介して互いに接続された演算装置１２０と、主記憶装置１３０と、補助記憶装置１４０と、入出力インタフェース１５０とを備えている。演算装置１２０として利用可能なデバイスとしては、ＣＰＵ（Central Processing Unit）を挙げることができる。また、主記憶装置１３０として利用可能なデバイスとしては、例えば、半導体ＲＡＭ（random access memory）を挙げることができる。また、補助記憶装置１４０として利用可能なデバイスとしては、例えば、ハードディスクドライブを挙げることができる。
【００５４】
入出力インタフェース１５０には、図２に示したように、入力装置２００及び出力装置３００が接続される。ユーザにより指定された数値を入力データ（標本ベクトル及び評価対象ベクトル）とする場合には、例えば、キーボードを入力装置２００として利用することができる。また、測定値を入力データとする場合には、例えば、その測定値に応じた測定器（センサ）を入力装置２００として利用することもできる。また、出力データ（特徴ベクトル）を数値列又はグラフとしてユーザに提示する場合、例えば、モニタやプリンタなどを出力装置３００として利用することができる。
【００５５】
なお、出力装置３００を介して出力データを出力する代わりに、記録媒体に出力データを記録するようにしてもよい。また、出力データを出力装置３００に出力する代わりに、他の装置に出力データを送信するようにしてもよい。
【００５６】
補助記憶装置１４０には、コンピュータ１００をＳＫＰＣＡ装置１として動作させるためのカーネル主成分分析プログラム（以下「ＫＰＣＡプログラム」と略記する）が格納されている。ＫＰＣＡプログラムは、選択プログラムと、自己カーネルグラム行列計算プログラムと、相互カーネルグラム行列計算プログラムと、一般化固有値問題計算プログラムと、経験カーネル写像ベクトル計算プログラムと、特徴ベクトル計算プログラムとを含んでいる。
【００５７】
演算装置１２０は、主記憶装置１３０上に展開された標本プログラムに含まれる命令を実行することによって、コンピュータ１００を選択部１０として機能させる。同様に、主記憶装置１３０上に展開された自己カーネルグラム行列計算プログラム、相互カーネルグラム行列計算プログラム、一般化固有値問題計算プログラム、経験カーネル写像ベクトル計算プログラム、及び、特徴ベクトル計算プログラムに含まれる命令を実行することによって、コンピュータ１００を自己カーネルグラム行列計算部２１、相互カーネルグラム行列計算部２２、一般化固有値問題計算部２３、経験カーネル写像ベクトル計算部３１、及び、特徴ベクトル計算部３２として機能させる。
【００５８】
ＫＰＣＡプログラムを実行する過程で生成される中間データ、すなわち、標本ベクトル｛ｙ1，ｙ2，…，ｙm｝、相互カーネルグラム行列Ｋxy、自己カーネルグラム行列Ｋy、相互カーネルグラム行列の転置行列ＫxyT、相互カーネルグラム行列Ｋxyとその転置行列ＫxyTとの行列積ＫxyTＫxy、固有ベクトル｛ｚ1，ｚ2，…，ｚr｝、及び、経験カーネル写像ベクトルｘΦは、演算装置１２０がこれを参照することができるよう主記憶装置１３０に格納される。ただし、主記憶装置１３０の記憶容量に制約がある場合には、これらの中間データを補助記憶装置１４０に格納するようにしてもよい。他のデータと比べてサイズの大きい相互カーネルグラム行列Ｋxyだけを補助記憶装置１４０に格納するような構成も有効である。
【００５９】
（カーネル主成分分析方法）
次に、部分カーネル主成分分析方法（以下、「ＳＫＰＣＡ方法」と略記する）について、図３〜図６を参照して説明する。
【００６０】
図３は、本実施形態に係るＳＫＰＣＡ方法Ｓ１の流れを示したフローチャートである。ＳＫＰＣＡ方法Ｓ１は、ＳＫＰＣＡ装置１を用いたＳＫＰＣＡ方法であり、図３に示したように、選択ステップＳ１０と、設計ステップＳ２０と、評価ステップＳ３０とを含んでいる。また、図３に示したように、標本ベクトル入力ステップＳin１と、評価対象ベクトル入力ステップＳin２と、特徴ベクトル出力ステップＳoutとを更に含んでいる。
【００６１】
標本ベクトル入力ステップＳin１においては、ＳＫＰＣＡ装置１が、ｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝の入力を受け付ける。換言すれば、ｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝を入力装置等の外部装置から取得する。
【００６２】
選択ステップＳ１０においては、ＳＫＰＣＡ装置１の選択部１０が、標本ベクトル入力ステップＳin１にて入力されたｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝からｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝を選択する。
【００６３】
設計ステップＳ２０においては、ＳＫＰＣＡ装置１の設計部２０が、一般化固有値問題ＫxyTＫxyｚ＝λＫyzを解くことによって、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出する。自己カーネルグラム行列Ｋy、及び、相互カーネルグラム行列Ｋxyは、上述したとおり、標本ベクトル入力ステップＳin１にて入力されたｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝、及び、選択ステップＳ１０にて選択されたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝に基づいて定義することができる行列である。
【００６４】
評価対象ベクトル入力ステップＳin２においては、ＳＫＰＣＡ装置１が、評価対象ベクトルｘの入力を受け付ける。換言すれば、評価対象ベクトルｘを入力装置等の外部装置から取得する。なお、評価対象ベクトル入力ステップＳin２は、評価ステップＳ３０より前に実行されればよく、設計ステップＳ２０より後に実行されることを要さない。すなわち、設計ステップＳ２０より前に、あるいは、選択ステップＳ１０より前に実行されてもよい。
【００６５】
評価ステップＳ３０においては、ＳＫＰＣＡ装置１の評価部３０が、設計ステップＳ２０にて算出されたｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝とを用いて、評価対象ベクトル入力ステップＳin２にて入力された評価対象ベクトルｘを特徴ベクトルｙ＝ＺTｘΦに変換する。
【００６６】
特徴ベクトル出力ステップＳoutにおいては、ＳＫＰＣＡ装置１が、評価ステップＳ３０にて算出された特徴ベクトルｙを出力する。例えば、モニタを用いて特徴ベクトルｙを数値列又はグラフとして表示したり、プリンタを用いて特徴ベクトルｙを数値又はグラフとして印刷したりする。なお、特徴ベクトルｙを表示したり印刷したりする代わりに、特徴ベクトルｙを記録媒体に記録したり他の装置に送信したりしてもよい。
【００６７】
なお、特徴ベクトルｙをグラフとして表示または印刷する場合、特徴ベクトルｙの次元が１であれば直線上の点としてこれをプロットすることができるし、特徴ベクトルｙの次元が２であれば平面上の点としてこれをプロットすることができるし、特徴ベクトルｙの次元が３であれば、（仮想）３次元空間上の点としてこれをプロットすることができる。この場合、評価対象ベクトルｘと合わせて表示すれば、評価対象ベクトルｘと特徴ベクトルｙとの関係も一目瞭然である。
【００６８】
図４は、選択ステップＳ１０の構成例を示したフローチャートである。選択ステップＳ１０は、例えば図４に示したように、クラスタリングステップＳ１１と、標本ベクトル選択ステップＳ１２とにより構成することができる。
【００６９】
クラスタリングステップＳ１１においては、選択部１０が、標本ベクトル入力ステップＳin１にて入力されたｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝をｍ個のクラスタ｛Ｃ1，Ｃ2，…，Ｃm｝にクラスタリングする。
【００７０】
標本ベクトル選択ステップＳ１２においては、選択部１０が、クラスタリングステップＳ１１にてクラスタリングされたｍ個のクラスタ｛Ｃ1，Ｃ2，…，Ｃm｝の各々からｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝の各々を選択する。
【００７１】
なお、クラスタリングステップＳ１１においては、Φ（ｘi）が属する特徴空間Ｆにおける距離を用いたクラスタリングを行うことができる。これにより、標本ベクトルｘiが属する空間Ｒdにおける距離を用いたクラスタリングを行う場合と比べて、より精度の高い特徴ベクトルを得ることができる。特徴空間Ｆにおける距離としては、例えば、‖Φ（ｘi）−Φ（ｘj）‖2＝ｋ（ｘi，ｘi）＋ｋ（ｘj，ｘj）−２ｋ（ｘi，ｘj）が挙げられる。
【００７２】
クラスタリングステップＳ１１において利用可能な公知のクラスタリング手法としては、階層化クラスタリングやK-meansなどを挙げることができる。例えば、K-meansを利用した場合、（１）初期クラスタ重心｛＜Ｃ1＞，＜Ｃ2＞，…，＜Ｃm＞｝をｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝から適当に選び出し、（２）各標本ベクトルｘｊを、特徴空間Ｆにおける距離‖Φ（ｘｊ）−＜Ｃｉ＞‖2が最小となるクラスタＣｉに割り振り、（３）各クラスタ重心＜Ｃｉ＞を、＜Ｃi＞＝（１／｜Ｃi｜）Σｘ∈ＣiΦ（ｘ）により更新する処理を、クラスタ重心｛＜Ｃ1＞，＜Ｃ2＞，…，＜Ｃm＞｝が更新されなくなるまで繰り返すことによって、ｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝をｍ個のクラスタ｛Ｃ1，Ｃ2，…，Ｃm｝にクラスタリングすることができる。
【００７３】
また、標本ベクトル選択ステップＳ１２においては、各クラスタＣiに属する標本ベクトルのうち、特徴空間ＦにおいてクラスタＣiの重心に最も近い標本ベクトル、すなわち、argminｘ∈Ｃi‖Φ（ｘ）−＜Ｃi＞‖2を取り出し、これを標本ベクトルｙiとすることが好ましい。あるいは、クラスタＣiに属する標本ベクトルのうち、入力空間Ｒdにおいてクラスタ重心Ｃiに最も近い標本ベクトル（いわゆるPre-image）を取り出し、これを標本ベクトルｙiとすることが好ましい。これにより、各クラスタから他の標本ベクトルを取り出す場合と比べて、より精度の高い特徴ベクトルを得ることができる。なお、Pre-imageについては、以下の文献を参照されたい。
【００７４】
S. Mika, B. Scholkopf and A. Smola，“Kernel PCA and denoising in feature space”，Advances in neural information processing systems (NIPS)，11，pp，536-542 (1999)。
【００７５】
B. Scholkopf，S. Mika，C. Burges，P. Knirsch，K.-R. Muller，G. Ratsch and A. Smola，“Input space vs. feature space in kernel-based methods”，IEEE Transactions on Neural Networks，10，5，pp. 1000-1017 (1999)。
【００７６】
M. Girolami，“Mercer kernel-based clustering in feature space”，IEEE trans. on neural networks，13，3，pp. 780-784(2002)。
【００７７】
なお、選択ステップＳ１０の構成は、図４に示したものに限らない。特に、ＳＫＰＣＡ方法Ｓ１では近似誤差あるいは残差を容易に算出することができるので、前向き選択や後向き選択などの標本選択の手法を容易に導入することができる。更に、無作為抽出されたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝を用いたＳＫＰＣＡ方法Ｓ１を繰り返し、所望の精度を有する特徴ベクトルｙを算出するような構成を採ることもできる。
【００７８】
図５は、設計ステップＳ２０の構成例を示したフローチャートである。設計ステップ２０は、例えば図５に示したように、自己カーネルグラム行列計算ステップＳ２１と、相互カーネルグラム行列計算ステップＳ２２と、一般化固有値問題計算ステップＳ２３とにより構成することができる。
【００７９】
自己カーネルグラム行列計算ステップＳ２１においては、評価部２０の自己カーネルグラム行列計算部２１が、選択ステップ１０にて選択されたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝から、自己カーネルグラム行列Ｋyを算出する。
【００８０】
相互カーネルグラム行列計算ステップＳ２２においては、評価部２０の相互カーネルグラム行列計算部２２が、標本ベクトル入力ステップＳin１にて入力されたｎ個の標本ベクトル｛ｘ1，ｘ2，…，ｘn｝と、選択ステップ１０にて選択されたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝とから、相互カーネルグラム行列Ｋxyを算出する。
【００８１】
一般化固有値問題計算ステップＳ２３においては、評価部２０の一般化固有値問題計算部２３が、自己カーネルグラム行列計算ステップＳ２１にて算出された自己カーネルグラム行列Ｋyと、相互カーネルグラム行列計算ステップＳ２２にて算出された相互カーネルグラム行列Ｋxyとにより定義される一般化固有値問題ＫxyTＫxyｚ＝λＫyzを解くことによって、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出する。
【００８２】
なお、自己カーネルグラム行列計算ステップＳ２１及び相互カーネルグラム行列計算ステップＳ２２は、一般化固有値問題計算ステップＳ２３より以前に実行されればよく、これらの実行順序は図５に示したものに限定されない。すなわち、図５においては自己カーネルグラム行列計算ステップＳ２１を実行した後で相互カーネルグラム行列計算ステップＳ２２を実行するようにしているが、相互カーネルグラム行列計算ステップＳ２２を実行した後で自己カーネルグラム行列計算ステップＳ２１を実行するようにしてもよい。
【００８３】
図６は、評価ステップＳ３０の構成例を示したフローチャートである。評価ステップＳ３０は、例えば図６に示したように、経験カーネル写像ベクトル計算ステップＳ３１と、特徴ベクトル計算ステップＳ３２と、により構成することができる。
【００８４】
経験カーネル写像ベクトル計算ステップＳ３１においては、評価部３０の経験カーネル写像ベクトル計算部３１が、評価対象ベクトル取得ステップＳin２にて入力された評価対象ベクトルｘと、選択ステップ１０にて選択されたｍ個の標本ベクトル｛ｙ1，ｙ2，…，ｙm｝とから、経験カーネル写像ベクトルｘΦを算出する。
【００８５】
特徴ベクトル計算ステップＳ３２においては、評価部３０の特徴ベクトル計算部３２が、一般化固有値問題計算ステップＳ２３にて算出されたｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝と、経験カーネル写像計算ステップＳ２１にて算出された経験カーネル写像ベクトルｘΦとから、特徴ベクトルｙを算出する。
【００８６】
〔部分カーネル主成分分析の応用〕
部分カーネル主成分分析は、カーネル主成分分析と同様、雑音除去、パターン分類、パターン検出、欠損値補完、統計的推定、統計的予測、脳コンピュータインタフェース（ＢＣＩ）などに応用することができる。これらの問題に部分カーネル主成分分析を適用する仕方は、これらの問題に従来のカーネル主成分分析を適用する仕方と同様である。
【００８７】
例えば、部分カーネル主成分分析を用いた雑音除去方法は、従来、カーネル主成分分析により得られた特徴ベクトルを用いて行われていた逆変換写像の計算を、部分カーネル主成分分析により得られた特徴ベクトルを用いて行うようにすることにより実現することができる。逆変換写像の計算に部分カーネル主成分分析を用いた場合、従来のカーネル主成分分析を用いた場合と比べて、少ない計算量で効率の良い雑音除去を行うことができる。
【００８８】
パターン分類、パターン検出、及び、欠損値補間においては、例えば、以下のようにして部分主成分分析方法を用いることができる。
【００８９】
パターン分類は、ｐ個のクラス｛Ａ1，Ａ2，…，Ａp｝に分類されているｑ個の標本ベクトル｛ｓ1，ｓ2，…，ｓq｝基づき、対象ベクトルｔを分類すべきクラスを決定する方法である。（１）各クラスＡiに分類されているｎi個の標本ベクトル｛ｓi1，ｓi2，…，ｓini｝に基づく部分カーネル主成分分析を実行することによって、対象ベクトルｔの特徴ベクトルｕiを各j=1,…,pについて算出し、（２）特徴ベクトルｕjのノルムが最大となるクラスＡjに対象ベクトルｔを分類すれば、少ない計算量で精度良くパターン分類を行うことができる。
【００９０】
パターン検出は、或るパターンを有するｑ個の標本ベクトル｛ｓ1，ｓ2，…，ｓq｝に基づき、対象ベクトルｔがそのパターンを有するか否かを判定する方法である。（１）
ｑ個の標本ベクトル｛ｓ1，ｓ2，…，ｓq｝に基づく部分カーネル主成分分析を実行することによって、対象ベクトルｔの特徴ベクトルｕを算出し、（２）特徴ベクトルｕのノルムを予め定められた閾値Ｔｈと比較することによって、対象ベクトルｔが特定のパターンを有するか否かを判定すれば、少ない計算量で精度良くパターン検出を行うことができる。
【００９１】
欠損値補完は、ｑ個の標本ベクトル｛ｓ1，ｓ2，…，ｓq｝に基づき、第ｋ成分に欠損のある対象ベクトルｔに対する欠損値補完を行う方法である。（１）第ｋ成分を異なる値に設定したｐ個の対象ベクトル｛ｔ1，ｔ2，…，ｔp｝の各々について、ｑ個の標本ベクトル｛ｘ1，ｘ2，…，ｘq｝に基づく部分カーネル主成分分析を実行することによって、対象ベクトルｔiの特徴ベクトルｕiを算出し、（２）特徴ベクトルｕjのノルムが最大となる対象ベクトルｔjの第ｋ成分の値を欠損のある対象ベクトルｔの第ｋ成分の値とすれば、少ない計算量で精度良く欠損値補間を行うことができる。
【００９２】
なお、与えられたベクトルデータを対象ベクトルｔと看做し、与えられたベクトルデータにおける未知の値を対象ベクトルｔにおける欠損成分と看做せば、対象ベクトルｔの欠損値補完を行うことにより、ベクトルデータにおける未知の値を統計的に推定することができる（統計的推定）。また、時系列データを対象ベクトルｔと看做し、時系列データにおける未来の値を対象ベクトルｔにおける欠損成分と看做せば、対象ベクトルｔの欠損値補完を行うことにより、時系列データにおける未来の値を統計的に予測することができる（統計的予測）。
【００９３】
更に、部分カーネル主成分分析は、脳コンピュータインタフェース（ＢＣＩ）に応用することができる。脳信号の測定には、脳電図（ＥＥＧ）、脳磁図（ＭＥＧ）、近赤外光トポグラフィー（ＮＩＲＳ）、ｆＭＲＩなどが用いられるが、ここでは、ＥＥＧを用いたＢＣＩについて説明する。
【００９４】
ＥＥＧによって得られた脳信号は、通常、多チャンネルであり、かつ、雑音を多く含む。そこで、ＢＣＩの前処理として、アーチファクト除去、バンドパスフィルタ（ＢＰＦ）、ＣＳＰ（Common Spatial Patterns)などの処理が施される。その後、入力信号は適当な時間幅をもつ窓で区切られ、データは数値列で表される。この数値列を評価対象ベクトルｘとして部分カーネル主成分分析を利用する。この際、標本ベクトル｛ｘ1，ｘ2，…，ｘn｝は予め与えておく。そして、部分カーネル主成分分析により得られた特徴ベクトルｙを、識別器で判別することによって、ＢＣＩの結果を得る。
【００９５】
〔部分カーネル主成分分析の原理〕
次に、本発明の基礎となる部分カーネル主成分分析の原理について説明する。
【００９６】
１．カーネル主成分分析
本節では、簡略にＫＰＣＡについて示した後、Mercerカーネル（１）を用いたＫＰＣＡの特徴付けを行う。
【００９７】
【数１】
【００９８】
ｘ1，ｘ2，…，ｘnをｄ次元の標本とする。特徴空間Ｆへ変換された標本Φ（ｘ1），…，Φ（ｘn）の相関作用素ＲFは、（２）で与えられる。
【００９９】
【数２】
【０１００】
ここで，ケット-ブラ表記｜・＞＜・｜は｜ａ＞＜ｂ｜ｃ＝＜ｃ，ｂ＞ａを満たす作用素である。有限次元では転置・Tを用いた表記ａｂTと同値であるが、Ｇａｕｓｓｉａｎカーネルなどを用いた場合は、Ｆが無限次元となるためケット-ブラ表記を用いる。ＲFの固有値分解を（３）とおく。
【０１０１】
【数３】
【０１０２】
ここで固有値λiは降順に並んでいるとする。ｒ次元の固有空間への変換ＵKと射影作用素ＰKは、（４）（５）で与えられる。
【０１０３】
【数４】
【０１０４】
【数５】
【０１０５】
ここでｅriは、Ｒrのｉ番目の標準基底（ｉ番目の要素が１でそれ以外の要素が０のベクトル）を表す。
【０１０６】
通常，ＲFは非常に大きいか、無限次元から無限次元への作用素であるため固有値分解が難しい。そこで、以下の作用素Ｓを考える。
【０１０７】
【数６】
【０１０８】
Ｓの随伴作用素をＳ*と書くと、ＲF＝（１／ｎ）ＳＳ*の固有値、固有ベクトルは、Ｋx＝Ｓ*Ｓ＝ｎΣni=1λiｖiｖiT、ｕi＝（１／√λi）Ｓｖiの関係を満たす。これより、（７）（８）が得られる。
【０１０９】
【数７】
【０１１０】
【数８】
【０１１１】
ここで、（Ｋx）ij＝ｋ（ｘi，ｘj）はカーネルグラム行列と呼ばれる。入力ｘに対して、以下の出力を得る。
【０１１２】
【数９】
【０１１３】
【数１０】
【０１１４】
ここで、‖・‖はＥｕｃｌｉｄノルムまたはｌ2ノルムを表し、Ｓ*Φ（ｘ）＝［ｋ（ｘ，ｘ1），…，ｋ（ｘ，ｘｎ）］Tは経験カーネル写像と呼ばれるｎ次元ベクトルである。
【０１１５】
以上の議論より、ＵKあるいはＰKを求めるためには、ｎ×ｎの大きさの対称行列の固有値分解が必要であり，経験カーネル写像Ｓ*Φ（ｘ）を求めるためにはｎ回のカーネル関数の計算が必要となる。さらに入力ｘの評価のために、すべての標本を保存しておく必要がある。このため、ｎが大きくなると計算量、メモリ量が問題となることがある。特に標本数が数万から数十万を超えるような大規模な問題ではＫＰＣＡを適用することができない。
【０１１６】
ＫＰＣＡを拡張するために，ＫＰＣＡの特徴付けを行う。主成分分析は、ランク制約の下で近似誤差を最小とする行列と定義できる。
【０１１７】
【数１１】
【０１１８】
ここで、Ｅxは平均を表す。ＫＰＣＡについても同様の特徴付けができる。
【０１１９】
【数１２】
【０１２０】
ここで、Ｎ（Ａ）はＡの核空間、Ｒ（Ａ）はＡの値域あるいは像を表す。制約Ｎ（Ｘ）⊃Ｒ（Ｓ）⊥は特徴空間Ｆ中のＳが張る空間以外についても考慮する必要があるため必要となる。入力次元空間の問題では、標本数が充分にあれば全空間を張ることが多いため、核空間や値域の議論を省略することが多い。しかしながら、高次元の特徴空間では、標本が全空間を張らないことが多いため、これらを正確に扱うことが重要である。本明細書においてもこれらを正確に扱う。
【０１２１】
２. 部分カーネル主成分分析
前節で示した通り、ＫＰＣＡは特徴空間Ｆで、標本Φ（ｘ1），…，Φ（ｘn）が張るＲnと同型の空間上で問題が議論される。すなわち。標本Φ（ｘ1），…，Φ（ｘn）は基底を成す。ＫＰＣＡはすべての標本を基底に用いるため、固有値分解の大きさやカーネル関数の評価回数が大きい。そこで、基底を少数の標本に置き換えることにより、この問題を解決する。
【０１２２】
ｙ1，…，ｙm（ｍ＜ｎ）を基底のための標本とし、作用素（１１）を考える。
【０１２３】
【数１３】
【０１２４】
この基底の下で全標本を最良に低階数近似する作用素を部分カーネル主成分分析（ＳＫＰＣＡ：ＳｕｂｓｅｔＫＰＣＡ）と定義する。
【０１２５】
［定義１］（ＳｕｂｓｅｔＫＰＣＡ）以下の最適化問題の解を部分ＫＰＣＡと呼ぶ。
【０１２６】
【数１４】
【０１２７】
Ｋy＝Ｔ*Ｔ、Ｋxy＝Ｓ*Ｔとおき、Ａ+で作用素または行列のＭｏｏｒｅ−Ｐｅｎｒｏｓｅ一般逆、Ａ1/2で自己共役作用素あるいは対称行列の平方根Ａ1/2ＡT/2＝Ａを表す。Ａ1/2は一般に対称ではなく、一意ではない。行列平方根やＣｈｏｌｅｓｋｙ分解などを用いることができる。次にＳＫＰＣＡの解を示す。
【０１２８】
［定理１］（ＳＫＰＣＡの解）（Ｋy1/2）+ＫxyTＫxy（ＫyT/2）+の固有値分解を（１３）とおく。
【０１２９】
【数１５】
【０１３０】
ここで、ξi (ｉ＝１，…，ｍ)は降順に並んでいるものとする。Ｗ＝［ｗ1，…，ｗr］とおくと、最適化問題（１２）の解の１つは、以下で与えられる。
【０１３１】
【数１６】
【０１３２】
特にＲ（Ｔ）⊂Ｒ（Ｓ）であるような場合、すなわち｛ｙi｝mi=1が｛ｘi｝ni=1の部分集合で与えられるような場合、（１４）で与えられる。
【０１３３】
【数１７】
【０１３４】
証明は後述する。
【０１３５】
［系１］Ｒ（Ｔ）⊂Ｒ（Ｓ）ならばＰSは正射影である。
【０１３６】
（証明）明らかにＰS*＝ＰSである。（Ｋy1/2）+Ｋy（ＫyT/2）+は、Ｒ（Ｋy）への正射影であり、Ｒ（Ｋy）⊃Ｒ（Ｗ）であることに注意すると、ＰSＰS＝Ｔ（ＫyT/2）+ＷＷT（Ｋy1/2）+Ｋy（ＫyT/2）+ＷＷT（Ｋy1/2）+Ｔ*＝ＰSであり、ＰSは正射影である(証明終)。
【０１３７】
［系２］基底にｘ1，…，ｘnを用いるとＳＫＰＣＡはＫＰＣＡに一致する。
【０１３８】
明らかであるため，証明は省略する。
【０１３９】
次にＳＫＰＣＡに別の表現を与える。
【０１４０】
［補題１］Ａ，ＢをＲ（Ａ）⊂ Ｒ（Ｂ）を満たす非負定値の行列とする。このとき、（Ｂ1/2）+Ａ（ＢT/2）+の固有値λ、と固有ベクトルｖは、任意のスカラーαに関して、Ａｕ＝λＢｕ、ｕ＝α（BT/2）+ｖを満たす。さらに、（λ，u）がＡｕ＝λＢｕを満たすとき、任意のスカラーβに対して、（Ｂ1/2）+Ａ（ＢT/2）+ｖ＝λｖ、ｖ＝βＢT/2ｕ を満たす。
【０１４１】
［命題１］補題より、Ｒ（Ｔ）⊂ Ｒ（Ｓ）ならば、一般固有値問題ＫxyTＫxyｚ＝λＫyｚの固有値の上位ｒ個に対応する固有ベクトル｛ｚi｝ri=1を用いて、（１５）と表される。
【０１４２】
【数１８】
【０１４３】
ここで、Ｚ＝［ｚ1…ｚr］であり、ｚ1，…，ｚrは，＜ｚi｜Ｋyｚi＞＝１，（ｉ＝１，…，ｒ）となるように正規化されているとする。
【０１４４】
Ｒ（Ｔ）⊂Ｒ（Ｓ）ならば、ＳＫＰＣＡによるｒ次元への変換は、ＵS＝ＺTＴ*で与えられ、入力をｘとすると（１６）である。
【０１４５】
【数１９】
【０１４６】
以上の議論より、ＳＫＰＣＡは大きさｍの（一般）固有値問題に帰着する。評価時についてもｍ回のカーネル関数の計算で評価でき，ｍ個の標本を保存しておくだけでよい。これ以降、Ｒ（Ｔ）⊂Ｒ（Ｓ）である、すなわちｙ1，…，ｙmはｘ1，…，ｘnの部分集合で与えられているとする。
【０１４７】
ＩＫＰＣＡでは、変換作用素は、異なる考え方から導き出され、（１７）で与えられている。
【０１４８】
【数２０】
【０１４９】
非特許文献２ではＩＫＰＣＡは識別のための特徴抽出に（１８）の形で用いられている。
【０１５０】
【数２１】
【０１５１】
すなわち、ＳＫＰＣＡと比較すると各次元の要素が（１／√ξi）倍されている。この影響は後段の識別性能にはあまり影響しないと考えられるが、機械学習のための特徴抽出以外の応用においては大きな誤差が生じる。
【０１５２】
３．正則化
標本数が少ないとき、統計的推定量はしばしば標本に過適合（overfit）することがある。特に高次元中で問題を解くカーネルトリックは自由度が高いため過適合しやすい。線形逆問題や回帰問題（リッジ回帰）などにおいては、過適合を防ぐために（Tikhonov）正則化がよく用いられる。そこで、本手法にも過適合を防ぐために正則化を導入する。
【０１５３】
正則化は、パラメータベクトルあるいは行列のノルムを最小化することにより実現される。ＳＫＰＣＡの正則化を正則化ＳＫＰＣＡ（ＲＳＫＰＣＡ：ＲｅｇｕｌａｒｉｚｅｄＳＫＰＣＡ）と呼び、以下に定義を示す。
【０１５４】
［定義２］（正則化ＳＫＰＣＡ）
【０１５５】
【数２２】
【０１５６】
ここで、μ＞０は正則化パラメータ、‖・‖FはＦｒｏｂｅｎｉｕｓノルムを表す。目的関数の第２項は適合の度合いを表し、正則化パラメータμが適合度合いと近似誤差のトレードオフの調整を行う。
【０１５７】
正則化ＳＫＰＣＡの解を示す。
【０１５８】
［定理２］（正則化ＳＫＰＣＡの解）ＫxyTＫxy＋μＫyが正則であるとき、正則化ＳＫＰＣＡの解の１つは、（２０）で与えられる。
【０１５９】
【数２３】
【０１６０】
ここで、ｑiは固有値分解（２１）の固有ベクトルとして与えられる。
【０１６１】
【数２４】
【０１６２】
ここで、固有値σiは降順に並んでいるものとする。ＫxyTＫxy＋μＫyが特異の場合には、逆行列は一般逆行列に置き換わる。
【０１６３】
ＳＫＰＣＡと同様に、（ＫyT/2）+ｑiは、一般固有値問題ＫxyTＫxy（ＫxyTＫxy＋μＫy）-1ＫxyTＫxyｐ＝λＫyｐで計算することができる。
【０１６４】
［命題２］Ｒ（Ｔ）⊂Ｒ（Ｓ）であるとき、ＫxyTＫxy＋μＫyが正則となる必要十分条件は、Ｎ（Ｔ）＝｛０｝である。
【０１６５】
（証明）（ｉ）必要条件を示す。ＫxyTＫxyとＫyは非負定値であるため、どちらかが正定値であれば，ＫxyTＫxy＋μＫyは正則である。Ｎ（Ｔ）＝｛０｝であれば明らかにＫyは正則である。（ｉｉ）十分条件を示す。正則であれば明らかに正定値であり、任意のｘ∈Ｒm，ｘ≠０に対して＜ｘ｜（ＫxyTＫxy＋μＫy）ｘ＞＝‖Ｓ*Ｔｘ‖2＋μ‖Ｔｘ‖2が正であるため、Ｎ（Ｔ）＝｛０｝である（証明終）。
【０１６６】
［命題３］ＧａｕｓｓｉａｎカーネルＫ（ｘ，ｙ）＝ｅｘｐ（−ｃ‖ｘ−ｙ‖2）を用いた場合、ｙ1，…，ｙmが互いに異なる場合、Ｋyは正定値である。
【０１６７】
以上の命題より、ＫxyTＫxy＋μＫyはGaussianカーネルを用いた場合は、標本が重複しなければ正則であり、また、多項式カーネルｋ（ｘ，ｙ）＝＜ｘ，ｙ＋ｃ＞dの場合においても、ｄが充分に大きければＫyは正則となる。
【０１６８】
正則化パラメータμは交差検定（cross validation）などにより実験的に決定出来る。教師なし学習として用いる場合は、学習用標本とは別に検定用標本ｚ1，…，ｚlを用意し、その検定用標本の近似誤差Σli=1‖Φ（ｚi）−ＰRΦ（ｚi）‖2を最小化するようなμを探す。検定用標本と学習用標本を入れ替え、交差検定を行うことで最適な正則化パラメータが推定できる。
【０１６９】
〔数値実験による検証〕
次に、本発明の効果を検証するための数値実験について説明する。
【０１７０】
１．人工データによる実験
本実験においては、図７の（ａ）に示す２次元の１０００個のデータを利用した。図７の（ｂ）は、Gaussianカーネルｋ（ｘ，ｙ）＝ｅｘｐ（−０．１‖ｘ−ｙ‖2）を用いたＫＰＣＡで５次元の部分空間への射影ノルムの等高線を表す。図７の（ｃ）に、標本を無作為に５０個選び出し、ＫＰＣＡを行った結果を表す。簡便のため、これ以降、（ｂ）のように標本すべてを用いたＫＰＣＡを全体ＫＰＣＡと呼び、（ｃ）のように標本の一部を使って求めたＫＰＣＡを縮減ＫＰＣＡと呼ぶ。図７の（ｄ）は縮減ＫＰＣＡと同じ標本を基底に用いたＳＫＰＣＡ、図７の（ｅ）は縮減ＫＰＣＡと同じ標本を基底に用いたＲＳＫＰＣＡ（μ＝０．１）、図７の（ｆ）はＩＫＰＣＡを示す。
【０１７１】
定量的評価には２つの手法を用いた。１つ目の評価基準は、全標本に対する近似誤差（残差）である。
【０１７２】
【数２５】
【０１７３】
Ｘは、ＰKやＰSで置き換えられる。もう１つの評価基準は全体ＫＰＣＡからの距離である。ｍ次元中のｒ次元空間部分空間への正射影行列の全体の集合はGrassmann多様体（Grassmannian）と呼ばれ、その空間における距離の１つは射影行列の差のノルムで与えられる。ＫＰＣＡとＳＫＰＣＡも特徴空間Ｆ中のｒ次元部分空間への正射影作用素であるため、差のノルムで距離を定義する。ノルムにはFrobeniusノルムを用いた。全体ＫＰＣＡの正射影作用素をＰFKとおき、距離の２乗を全体ＫＰＣＡのノルムの２乗‖ＰFK‖F2で割ったものを正規化距離Ｄとして定義する。
【０１７４】
【数２６】
【０１７５】
図より、提案手法は全体ＫＰＣＡとほとんど変わらない等高線を描くことが分かる。一方、縮減ＫＰＣＡとＩＫＰＣＡの等高線はＫＰＣＡから大きくかけ離れたものとなっている。近似誤差Ｅと正規化距離Ｄにおいても、提案手法は全体ＫＰＣＡとほぼ変わらないのに対して、縮減ＫＰＣＡ、ＩＫＰＣＡはどちらも大きな値となっている。
【０１７６】
Ｉｎｔｅｌ Ｃｏｒｅ２ Ｄｕｏ（登録商標） ２．１６ＧＨｚを搭載したＭａｃＢｏｏｋ（登録商標）でＭＡＴＬＡＢ（登録商標）のＩＲＡＭを利用した固有値分解の関数ｅｉｇｓでは、１０００×１０００の対称行列の固有値分解で０．４秒、５０×５０の対称行列の固有値分解で０．０２秒であった。この実験の場合では、提案法と縮減ＫＰＣＡは、全体ＫＰＣＡの場合より２０倍高速に固有値分解を行うことができる。さらにカーネル関数の評価回数も２０ 分の１の回数で計算ができる。
【０１７７】
２．実データを用いた実験
UCI Machine learning repositoryにあるconcreteとhousingのデータセットを用いて実験を行った。concreteは９次元、１０３０個の標本があり、housingは１４次元、５０６個の標本がある。カーネル関数にはGaussianカーネルを用い、パラメータｃには、全要素の分散σ から、ｃ＝１／２σ2を求めて用いた。階数は入力次元と同じ次元で実験を行った。
【０１７８】
図８の（ａ）及び（ｂ）は全体ＫＰＣＡからの正規化距離を示す。ＩＫＰＣＡは、比較できない程度の大きな距離を示したため省略した。標本の部分集合は無作為抽出で行い、１００回の無作為抽出の平均、標準偏差、最小値を示す。図より、標本数が全体の２０％から１０％程度であっても提案手法は全体ＫＰＣＡから１％程度の誤差しか生じないことが確認できる。図９の（ａ）及び（ｂ）は式(２２)の誤差を全体ＫＰＣＡの誤差で割った正規化誤差を示す。標本数が２０％程度に減少した場合でも、正規化誤差の最小値は、全体ＫＰＣＡの誤差の１．０１倍程度であることが分かる。図１０の（ａ）及び（ｂ）は設計時の計算時間を示す。数値計算はIntel Math Kernel Library（登録商標）を用い、全体ＫＰＣＡ、縮減ＫＰＣＡの固有値分解にはdsyevx、ＳＫＰＣＡの一般固有値分解には、dsygvxを用いた。図より、基底に用いる標本数が全標本数の７０％以下の場合、提案手法の方が高速に設計を行うことができることがわかる。
【０１７９】
図１１はhousingにおける無作為抽出と特徴空間中Ｋ‐ｍｅａｎｓによる標本選択の比較である。Ｋ‐ｍｅａｎｓの方が平均的によい近似性能を示し、標準偏差も小さいことが確認できる。
【０１８０】
最後に、正則化の効果を示すためにhousingの部分集合として１００個（ｎ＝１００）、基底にｍ個の標本を用いて、この部分集合に対する近似を最小にするＳＫＰＣＡ及び、正則化ＳＫＰＣＡを求め、それぞれの全体ＫＰＣＡからの正規化距離を比較した。１００回の試行における正規化距離の平均値を表１に示す。実験結果より、正則化を用いた方がよい近似を与えることが確認できた。
【０１８１】
【表１】
【０１８２】
〔定理の証明〕
［補題２］制約条件Ｎ（Ｘ）⊃Ｒ（Ｔ）⊥は、ある作用素Ａ：Ｒm→Ｆを用いて、Ｘ＝ＡＴ*と表されることと同値である。
【０１８３】
【数２７】
【０１８４】
（証明）（ｉ）<＝はＲ（Ｔ）⊥＝Ｎ（Ｔ）⊂Ｎ（Ｘ）より明らかである。（ｉｉ）⇒の方向、すなわち、「Ｎ（Ｘ）⊃Ｒ（Ｔ）⊥のときＸ＝ＡＴ*となるＡが存在する。」を示す。Ａ＝Ｘ（Ｔ*）+とおけばＡＴ*＝Ｘ（Ｔ*）+Ｔ*＝ＸＰR(T)＝Ｘ（Ｉ−ＰR(T)⊥）である。ここで、Ｐχはχへの正射影を表す。Ｎ（Ｘ）⊃Ｒ（Ｔ）⊥より、ＸＰR(T)⊥＝０であり、ＡＴ*＝Ｘとなる。（ｉ），（ｉｉ）より補題を得る（証明終）。
【０１８５】
［補題３］制約条件Ｒ（Ｘ）⊂Ｒ（Ｔ）⊥は、ある作用素Ａ：Ｆ→Ｒmを用いて、Ｘ＝ＴＡと表されることと同値である。
【０１８６】
証明は、補題２と同様なので省略する。
【０１８７】
定理２でμ＝０とおけば定理１となるため、μ≧０として定理２のみを証明する。補題２，３よりＸ＝ＴＡＴ*とおける。作用素Ａに対するトレースを任意の正規直交基底ｕi，ｉ＝１，…を用いてＴｒａｃｅ（Ａ）＝Σi＜ｕi｜Ａｕi＞と定義する。最適化問題（１９）の目的関数をＪとおくと以下のとおりである。
【０１８８】
【数２８】
【０１８９】
Ｒ（ＫxyT）＝Ｒ（Ｔ*Ｓ）⊂Ｒ（Ｔ*）＝Ｒ（Ｋy）より、ＫxyT＝ＰR(Ky)ＫxyT＝ＫyＫy+ＫxyTであり、Ｋxy＝Ｋxy（ＫyT/2）+ＫyT/2である。また、Ｒ（ＫxyTＫxy＋μＫy）⊃Ｒ（ＫxyT）より、ＰR(KxyTKxy+μKy)ＫxyT＝（ＫxyTＫxy＋μＫy）1/2（（ＫxyTＫxy＋μＫy）1/2）+ＫxyT＝ＫxyTである。これより、以下のようになる。
【０１９０】
【数２９】
【０１９１】
ここでＣはＡに対する定数項を表す。式（２１） の固有値分解より、特異値分解を以下のように置く。
【０１９２】
【数３０】
【０１９３】
Schmidtの近似定理（Eckart−Youngの定理）より、以下を満たすＡが存在するときは、上式を満たすＡがＪを最小化する。
【０１９４】
【数３１】
【０１９５】
【数３２】
【０１９６】
以上より、式（Ａ．１）を満たすＡが存在し、その解の一つは、以下で与えられる。
【０１９７】
【数３３】
【０１９８】
証明終
【産業上の利用可能性】
【０１９９】
本発明は、カーネル主成分分析を用いたデータ解析一般に広く適用することができる。特徴抽出、雑音除去、パターン分類、パターン検出、欠損値補完、統計的推定、統計的予測などに好適に利用することができる。
【符号の説明】
【０２００】
１ ＳＫＰＣＡ装置（カーネル主成分分析装置）
１０ 選択部（選択手段）
２０ 設計部（設計手段）
２１ 自己カーネルグラム行列計算部
２２ 相互カーネルグラム行列計算部
２３ 一般化固有値問題計算部
３０ 評価部（評価手段）
３１ 経験カーネル写像ベクトル計算部
３２ 特徴ベクトル計算部
【特許請求の範囲】
【請求項１】
カーネル主成分分析装置を用いて、標本として与えられたｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝に基づき、評価対象として与えられたｄ次元ベクトルｘの特徴を表現するｒ次元ベクトルｙを算出するカーネル主成分分析方法であって、
上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝からｍ個（ｍはｎ未満の自然数）のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を選択する選択ステップ、又は、上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝とは別にｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を取得する取得ステップと、
カーネル関数値ｋ（ｙi，ｙj）を第ｉｊ成分とするｍ×ｍ行列をＫyとし、カーネル関数値ｋ（ｘi，ｙj）を第ｉｊ成分とするｎ×ｍ行列をＫxyとして、一般化固有値問題ＫxyTＫxyｚ＝λＫyz、又は、正則化一般化固有値問題ＫxyTＫxy（ＫxyTＫxy＋μＫy）ＫxyTＫxyｚ＝λＫyzを解くことによって、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出する設計ステップと、
上記ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を用いて、上記ｄ次元ベクトルｘを上記ｒ次元ベクトルｙに変換する評価ステップと、を含んでいる、
ことを特徴とするカーネル主成分分析方法。
【請求項２】
上記設計ステップは、
上記ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝から、上記ｍ×ｍ行列Ｋyを算出するステップと、
上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝、及び、上記ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝から、上記ｎ×ｍ行列Ｋxyを算出するステップと、
上記一般化固有値問題ＫxyTＫxyｚ＝λＫyz、又は、上記正則化一般化固有値問題ＫxyTＫxy（ＫxyTＫxy＋μＫy）ＫxyTＫxyｚ＝λＫyzを解くことによって、上記ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出するステップと、を含んでおり、
上記評価ステップは、
上記ｄ次元ベクトルｘ、及び、上記ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝から、カーネル関数値ｋ（ｘ，ｙj）を第ｊ成分とするｍ次元ベクトルｘΦを算出するステップと、
上記ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝、及び、上記ｍ次元ベクトルｘΦから、上記ｒ次元ベクトルｙ＝［ｚ1，ｚ2，…，ｚr］TｘΦを算出するステップと、を含んでおり、
当該カーネル主成分分析方法は、上記評価ステップにて算出された上記ｒ次元ベクトルｙを他の装置に出力する出力ステップを更に含んでいる、
ことを特徴とする請求項１に記載のカーネル主成分分析方法。
【請求項３】
上記選択ステップは、
上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝をｍ個のクラスタ｛Ｃ1，Ｃ2，…，Ｃm｝にクラスタリングするステップと、
上記ｍ個のクラスタ｛Ｃ1，Ｃ2，…，Ｃm｝の各々から上記ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝の各々を選択するステップと、を含んでいる、
ことを特徴とする請求項１または２に記載のカーネル主成分分析方法。
【請求項４】
上記選択ステップにおいて、上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝をｍ個のクラスタ｛Ｃ1，Ｃ2，…，Ｃm｝にクラスタリングする際に、特徴空間上のノルムを使用する、
ことを特徴とする請求項３に記載のカーネル主成分分析方法。
【請求項５】
上記ベクトルｙの次元ｒは、上記ベクトルｘの次元ｄよりも小さい、
ことを特徴とする請求項１から４までの何れか１項に記載のカーネル主成分分析方法。
【請求項６】
上記ベクトルｙの次元ｒは、上記ベクトルｘの次元ｄよりも大きい、
ことを特徴とする請求項１から４までの何れか１項に記載のカーネル主成分分析方法。
【請求項７】
標本として与えられたｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝に基づき、評価対象として与えられたｄ次元ベクトルｘの特徴を表現するｒ次元ベクトルｙを算出するカーネル主成分分析装置であって、
上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝からｍ個（ｍはｎ未満の自然数）のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を選択する選択手段、又は、上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝とは別にｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を取得する取得手段と、
カーネル関数値ｋ（ｙi，ｙj）を第ｉｊ成分とするｍ×ｍ行列をＫyとし、カーネル関数値ｋ（ｘi，ｙj）を第ｉｊ成分とするｎ×ｍ行列をＫxyとして、一般化固有値問題ＫxyTＫxyｚ＝λＫyz、又は、正則化一般化固有値問題ＫxyTＫxy（ＫxyTＫxy＋μＫy）ＫxyTＫxyｚ＝λＫyzを解くことによって、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出する設計手段と、
上記ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を用いて、上記ｄ次元ベクトルｘを上記ｒ次元ベクトルｙに変換する評価手段と、とを備えている、
ことを特徴とするカーネル主成分分析装置。
【請求項８】
コンピュータを請求項７に記載のカーネル主成分分析装置として動作させるためのためのカーネル主成分分析プログラムであって、上記コンピュータを上記カーネル主成分分析装置が備えている各手段として機能させるカーネル主成分分析プログラム。
【請求項１】
カーネル主成分分析装置を用いて、標本として与えられたｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝に基づき、評価対象として与えられたｄ次元ベクトルｘの特徴を表現するｒ次元ベクトルｙを算出するカーネル主成分分析方法であって、
上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝からｍ個（ｍはｎ未満の自然数）のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を選択する選択ステップ、又は、上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝とは別にｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を取得する取得ステップと、
カーネル関数値ｋ（ｙi，ｙj）を第ｉｊ成分とするｍ×ｍ行列をＫyとし、カーネル関数値ｋ（ｘi，ｙj）を第ｉｊ成分とするｎ×ｍ行列をＫxyとして、一般化固有値問題ＫxyTＫxyｚ＝λＫyz、又は、正則化一般化固有値問題ＫxyTＫxy（ＫxyTＫxy＋μＫy）ＫxyTＫxyｚ＝λＫyzを解くことによって、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出する設計ステップと、
上記ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を用いて、上記ｄ次元ベクトルｘを上記ｒ次元ベクトルｙに変換する評価ステップと、を含んでいる、
ことを特徴とするカーネル主成分分析方法。
【請求項２】
上記設計ステップは、
上記ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝から、上記ｍ×ｍ行列Ｋyを算出するステップと、
上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝、及び、上記ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝から、上記ｎ×ｍ行列Ｋxyを算出するステップと、
上記一般化固有値問題ＫxyTＫxyｚ＝λＫyz、又は、上記正則化一般化固有値問題ＫxyTＫxy（ＫxyTＫxy＋μＫy）ＫxyTＫxyｚ＝λＫyzを解くことによって、上記ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出するステップと、を含んでおり、
上記評価ステップは、
上記ｄ次元ベクトルｘ、及び、上記ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝から、カーネル関数値ｋ（ｘ，ｙj）を第ｊ成分とするｍ次元ベクトルｘΦを算出するステップと、
上記ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝、及び、上記ｍ次元ベクトルｘΦから、上記ｒ次元ベクトルｙ＝［ｚ1，ｚ2，…，ｚr］TｘΦを算出するステップと、を含んでおり、
当該カーネル主成分分析方法は、上記評価ステップにて算出された上記ｒ次元ベクトルｙを他の装置に出力する出力ステップを更に含んでいる、
ことを特徴とする請求項１に記載のカーネル主成分分析方法。
【請求項３】
上記選択ステップは、
上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝をｍ個のクラスタ｛Ｃ1，Ｃ2，…，Ｃm｝にクラスタリングするステップと、
上記ｍ個のクラスタ｛Ｃ1，Ｃ2，…，Ｃm｝の各々から上記ｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝の各々を選択するステップと、を含んでいる、
ことを特徴とする請求項１または２に記載のカーネル主成分分析方法。
【請求項４】
上記選択ステップにおいて、上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝をｍ個のクラスタ｛Ｃ1，Ｃ2，…，Ｃm｝にクラスタリングする際に、特徴空間上のノルムを使用する、
ことを特徴とする請求項３に記載のカーネル主成分分析方法。
【請求項５】
上記ベクトルｙの次元ｒは、上記ベクトルｘの次元ｄよりも小さい、
ことを特徴とする請求項１から４までの何れか１項に記載のカーネル主成分分析方法。
【請求項６】
上記ベクトルｙの次元ｒは、上記ベクトルｘの次元ｄよりも大きい、
ことを特徴とする請求項１から４までの何れか１項に記載のカーネル主成分分析方法。
【請求項７】
標本として与えられたｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝に基づき、評価対象として与えられたｄ次元ベクトルｘの特徴を表現するｒ次元ベクトルｙを算出するカーネル主成分分析装置であって、
上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝からｍ個（ｍはｎ未満の自然数）のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を選択する選択手段、又は、上記ｎ個のｄ次元ベクトル｛ｘ1，ｘ2，…，ｘn｝とは別にｍ個のｄ次元ベクトル｛ｙ1，ｙ2，…，ｙm｝を取得する取得手段と、
カーネル関数値ｋ（ｙi，ｙj）を第ｉｊ成分とするｍ×ｍ行列をＫyとし、カーネル関数値ｋ（ｘi，ｙj）を第ｉｊ成分とするｎ×ｍ行列をＫxyとして、一般化固有値問題ＫxyTＫxyｚ＝λＫyz、又は、正則化一般化固有値問題ＫxyTＫxy（ＫxyTＫxy＋μＫy）ＫxyTＫxyｚ＝λＫyzを解くことによって、ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を算出する設計手段と、
上記ｒ個の固有ベクトル｛ｚ1，ｚ2，…，ｚr｝を用いて、上記ｄ次元ベクトルｘを上記ｒ次元ベクトルｙに変換する評価手段と、とを備えている、
ことを特徴とするカーネル主成分分析装置。
【請求項８】
コンピュータを請求項７に記載のカーネル主成分分析装置として動作させるためのためのカーネル主成分分析プログラムであって、上記コンピュータを上記カーネル主成分分析装置が備えている各手段として機能させるカーネル主成分分析プログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【公開番号】特開２０１１−２２９１２（Ｐ２０１１−２２９１２Ａ）
【公開日】平成２３年２月３日（２０１１．２．３）
【国際特許分類】
【出願番号】特願２００９−１６９０７２（Ｐ２００９−１６９０７２）
【出願日】平成２１年７月１７日（２００９．７．１７）
【出願人】（５０３３５９８２１）独立行政法人理化学研究所 (1,056)
【Ｆターム（参考）】
【公開日】平成２３年２月３日（２０１１．２．３）
【国際特許分類】
【出願日】平成２１年７月１７日（２００９．７．１７）
【出願人】（５０３３５９８２１）独立行政法人理化学研究所 (1,056)
【Ｆターム（参考）】
[ Back to top ]