(57)【要約】      （修正有）【課題】ＬＬＥ（Locally  Linear  Embedding）の計算に要するコストを低減させる。【解決手段】ＬＬＥ計算装置１０は、多次元行列のＳＶＤを計算し、複数のデータポイントのそれぞれについて、計算対象のデータポイントとの間の近似の距離を、ＳＶＤに基づいて計算し、近似の距離が所定値以下であるデータポイントのそれぞれと、計算対象のデータポイントとの間のユークリッド距離の推定値を計算し、ユークリッド距離の推定値が所定値以下であるデータポイントのそれぞれと、計算対象のデータポイントとの間のユークリッド距離を計算し、当該計算したユークリッド距離が所定値以下であるデータポイントを、計算対象のデータポイントの近傍のデータポイントに確定し、計算対象のデータポイントと計算対象のデータポイントの近傍のデータポイントのそれぞれとの間のエッジの重みに基づいて、ＬＬＥカーネルの固有ベクトルを計算する。【選択図】図１
2018190251.tif 000002
【特許請求の範囲】
【請求項１】
複数のデータポイントを有する多次元行列の特異値分解を計算する特異値分解計算部と、
前記複数のデータポイントから計算対象のデータポイントを選択するデータポイント選択部と、
前記複数のデータポイントのそれぞれについて、前記計算対象のデータポイントとの間の近似の距離を、前記特異値分解に基づいて計算する近似距離計算部と、
前記近似の距離が所定値以下であるデータポイントのそれぞれと、前記計算対象のデータポイントとの間のユークリッド距離の推定値を計算する距離推定部と、
前記ユークリッド距離の推定値が所定値以下であるデータポイントのそれぞれと、前記計算対象のデータポイントとの間のユークリッド距離を計算し、当該計算したユークリッド距離が所定値以下であるデータポイントを、前記計算対象のデータポイントの近傍のデータポイントに確定する近傍確定部と、
前記計算対象のデータポイントと前記計算対象のデータポイントの近傍のデータポイントのそれぞれとの間のエッジの重みを計算するエッジ重み計算部と、
前記エッジの重みに基づいて、前記複数のデータポイントについてのＬＬＥカーネルの固有ベクトルを計算するベクトル計算部と、
を有することを特徴とするＬＬＥ計算装置。
【請求項２】
前記エッジ重み計算部は、前記エッジの重みが計算済みである計算済みデータポイントの近傍データポイントに、前記計算対象のデータポイントの近傍のデータポイントのうちの少なくとも一部が含まれる場合、前記計算済みデータポイントについての前記エッジの重みの計算結果に基づいて、前記計算対象のデータポイントと前記計算対象のデータポイントの近傍のデータポイントのそれぞれとの間のエッジの重みを計算することを特徴とする請求項１に記載のＬＬＥ計算装置。
【請求項３】
前記エッジ重み計算部は、Ｗｏｏｄｂｕｒｙの公式を用いて、前記計算済みデータポイントの近傍データポイントに関する逆行列から、前記計算対象のデータポイントの近傍のデータポイントに関する逆行列を計算することを特徴とする請求項２に記載のＬＬＥ計算装置。
【請求項４】
前記複数のデータポイントについての前記エッジの重みを表す行列を、次数の昇順に並び替える並び替え部と、
前記複数のデータポイントの隣接行列のＬＵ分解を計算するＬＵ分解計算部と、
をさらに有し、
前記ベクトル計算部は、前記ＬＵ分解計算部によって計算されたＬＵ分解に基づいて前記ＬＬＥカーネルの固有ベクトルを計算することを特徴とする請求項１又は２に記載のＬＬＥ計算装置。
【請求項５】
ＬＬＥ計算装置によって実行されるＬＬＥ計算方法であって、
複数のデータポイントを有する多次元行列の特異値分解を計算する特異値分解計算工程と、
前記複数のデータポイントから計算対象のデータポイントを選択するデータポイント選択工程と、
前記複数のデータポイントのそれぞれについて、前記計算対象のデータポイントとの間の近似の距離を、前記特異値分解に基づいて計算する近似距離計算工程と、
前記近似の距離が所定値以下であるデータポイントのそれぞれと、前記計算対象のデータポイントとの間のユークリッド距離の推定値を計算する距離推定工程と、
前記ユークリッド距離の推定値が所定値以下であるデータポイントのそれぞれと、前記計算対象のデータポイントとの間のユークリッド距離を計算し、当該計算したユークリッド距離が所定値以下であるデータポイントを、前記計算対象のデータポイントの近傍のデータポイントに確定する近傍確定工程と、
前記計算対象のデータポイントと前記計算対象のデータポイントの近傍のデータポイントのそれぞれとの間のエッジの重みを計算するエッジ重み計算工程と、
前記エッジの重みに基づいて、前記複数のデータポイントについてのＬＬＥカーネルの固有ベクトルを計算するベクトル計算工程と、
を含んだことを特徴とするＬＬＥ計算方法。
【請求項６】
コンピュータを、請求項１から４のいずれか１項に記載のＬＬＥ計算装置として機能させるためのＬＬＥ計算プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、ＬＬＥ計算装置、ＬＬＥ計算方法及びＬＬＥ計算プログラムに関する。
【背景技術】
【０００２】
ＬＬＥ（Locally  Linear  Embedding）は非線形次元削減の代表的な手法の一つである。ＬＬＥの基本的アイデアはそれぞれのデータポイントをそれらの近傍のデータポイントによる回帰によって近似し、各データポイントの低次元への埋め込みを計算するというものである。ＬＬＥでは再構築コストと埋め込みコストを最小化する二つの最適化問題を解く。これらの最適化において乱数による初期値の生成や学習率の設定等は必要でないため、ＬＬＥは効果的に次元削減を行うことができる。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】Sam  T.  Roweis  and  Lawrence  K.  Saul,  Nonlinear  Dimensionality  Reduction  by  Locally  Linear  Embedding,  Science,  290(5500):2323-2326,  2000.
【発明の概要】
【発明が解決しようとする課題】
【０００４】
しかしながら、従来のＬＬＥには計算に要するコストが高いという問題があった。例えば、Ｋを各データポイントへの近傍のノードからのエッジの数とすると、従来のＬＬＥのアルゴリズムは各データポイントに対してＫ×Ｋの大きさのグラム行列を計算し、ラグランジュの未定乗数法を用いてグラム行列の逆行列からエッジの重みを計算する（例えば、非特許文献１を参照）。逆行列を計算するコストは各データポイントへのエッジの数の３乗であるため、Ｎをデータポイントの数としたときにエッジの重みを計算するコストはＯ（ＮＫ３）となる。
【０００５】
また各データポイントの低次元への埋め込みを計算するために従来のＬＬＥのアルゴリズムはＬＬＥカーネルの固有値ベクトルを計算する。ＬＬＥカーネルの大きさはＮ×Ｎであるため、固有ベクトルを計算するために必要なコストはＯ（Ｎ３）となる。さらに固有値ベクトルを保持するためにＯ（Ｎ２）のメモリコストが必要となる。そのため大規模なデータに対してＬＬＥを適用することは多くの計算コストとメモリコストが必要となるため非現実的であった。
【課題を解決するための手段】
【０００６】
本発明のＬＬＥ計算装置は、複数のデータポイントを有する多次元行列の特異値分解を計算する特異値分解計算部と、前記複数のデータポイントから計算対象のデータポイントを選択するデータポイント選択部と、前記複数のデータポイントのそれぞれについて、前記計算対象のデータポイントとの間の近似の距離を、前記特異値分解に基づいて計算する近似距離計算部と、前記近似の距離が所定値以下であるデータポイントのそれぞれと、前記計算対象のデータポイントとの間のユークリッド距離の推定値を計算する距離推定部と、前記ユークリッド距離の推定値が所定値以下であるデータポイントのそれぞれと、前記計算対象のデータポイントとの間のユークリッド距離を計算し、当該計算したユークリッド距離が所定値以下であるデータポイントを、前記計算対象のデータポイントの近傍のデータポイントに確定する近傍確定部と、前記計算対象のデータポイントと前記計算対象のデータポイントの近傍のデータポイントのそれぞれとの間のエッジの重みを計算するエッジ重み計算部と、前記エッジの重みに基づいて、前記複数のデータポイントについてのＬＬＥカーネルの固有ベクトルを計算するベクトル計算部と、を有することを特徴とする。
【発明の効果】
【０００７】
本発明によれば、ＬＬＥの計算に要するコストを低減させることができる。
【図面の簡単な説明】
【０００８】
【図１】図１は、第１の実施形態に係るＬＬＥ計算装置の構成の一例を示す図である。
【図２】図２は、補題１を示す図である。
【図３】図３は、補題２を示す図である。
【図４】図４は、補題３を示す図である。
【図５】図５は、第１の実施形態に係るＬＬＥ計算処理のアルゴリズムの一例を示す図である。
【図６】図６は、定理１、定理２及び定理３を示す図である。
【図７】図７は、第１の実施形態に係るＬＬＥ計算装置の処理の流れを示すフローチャートである。
【図８】図８は、ＬＬＥ計算プログラムを実行するコンピュータの一例を示す図である。
【発明を実施するための形態】
【０００９】
以下に、本願に係るＬＬＥ計算装置、ＬＬＥ計算方法及びＬＬＥ計算プログラムの実施形態を図面に基づいて詳細に説明する。なお、本発明は、以下に説明する実施形態により限定されるものではない。
【００１０】
［従来のＬＬＥ］
まず従来のＬＬＥの説明を行う。Ｘ∈ＲＮ×ＭをＭ次元からなるＮ個のデータポイントの行列とする。ｘｉ＝（ｘｉ［１］，ｘｉ［２］，．．．，ｘｉ［Ｍ］）を行列Ｘのｉ番目の行ベクトルとすると、ｘｉはｉ番目のデータポイントに対応する。
【００１１】
ＬＬＥは高次元のデータセットＸから低次元のデータセットＹ∈ＲＮ×ｍを計算する。ｙｉ＝（ｙｉ［１］，ｙｉ［２］，．．．，ｙｉ［ｍ］）を長さがｍの行ベクトルとすると、ｙｉはデータポイントｘｉの低次元へのマッピングに対応する。ＬＬＥにおける基本的考えはデータポイントは大域的に非線形な構造を持っていたとしても、各データポイントとそれらの近傍のデータポイントは局所的に線形に近似できるというものである。近傍から各データポイントを表現することでＬＬＥは高次元のデータを低次元で表現することができる。
【００１２】
従来のＬＬＥのアルゴリズムは再構築コストと埋め込みコストを最小化する二つの最適化問題を解くパートから構成される。はじめのパートでは各データポイントに対してユークリッド距離から近傍のデータポイントを計算し、それらの回帰分析によりエッジの重みを計算する。ｘｐを回帰分析によって表現するデータポイントとすると、従来のＬＬＥのアルゴリズムは式（１）の再構築コストを最小化する。
【００１３】
【数１】
2018190251.tif 000003
【００１４】
式（１）において｜｜・｜｜はベクトルのＬ２ノルムであり、〜ｘｐは式（２）のように計算されるｘｐの回帰分析の結果である。ここで〜・は・の直上に〜が記された記号を示すものとする。
【００１５】
【数２】
2018190251.tif 000004
【００１６】
式（２）においてＮ［ｘｐ］はデータポイントｘｐのＫ個の近傍のデータポイントの集合であり、Ｗ［ｉ，ｊ］はデータポイントｘｊからｘｉのエッジの重みである。重みを計算するために、式（１）においてデータポイントｘｐに対するエッジの重みの和が１になるように再構築コストを最小化する。再構築コストを最小化するためにグラム行列にラグランジュの未定乗数法を適用する。具体的にはｘｉとｘｊをデータポイントｘｐに対するＫ近傍とすると、データポイントｘｐに対するＫ×Ｋのグラム行列Ｇｐの要素は、Ｇｐ［ｉ，ｊ］をグラム行列Ｇｐの［ｉ，ｊ］成分としたときに式（３）のように計算される。
【００１７】
【数３】
2018190251.tif 000005
【００１８】
ここで＜・，・＞は二つのベクトルの内積とする。Ｇｐ−１を行列Ｇｐの逆行列としたときに各データポイントのエッジの重みは式（４）のように計算される。
【００１９】
【数４】
2018190251.tif 000006
【００２０】
式（４）のようにエッジの重みを計算するには逆行列Ｇｐ−１を求める必要がある。ＬＬＥのアルゴリズムの二つ目のパートでは式（５）の埋め込みコストを最小化することでデータポイントｘｐに対するベクトルｙｐを求める。
【００２１】
【数５】
2018190251.tif 000007
【００２２】
Ｉを単位行列としたときにＫ＝（Ｉ−Ｗ）Ｔ（Ｉ−Ｗ）をＮ×Ｎの大きさのＬＬＥカーネルとし、ＷをＮ×Ｎのエッジの重みの行列とすると、埋め込みコストはＹＴＫＹと書き換えることができる。ここで・Ｔは行列・の転置行列とする。
【００２３】
埋め込みコストには二つの制約がある。はじめの制約はＹＴＹ＝Ｉというものである。これは解がランクｍということである。また０を零ベクトルとしたときにΣｐｙｐ＝０であるというものである。よってＬＬＥのアルゴリズムの二つ目のパートは式（６）の最適化問題を解くこととなる。
【００２４】
【数６】
2018190251.tif 000008
【００２５】
式（６）はカーネルＫの底のｍ＋１個の固有ベクトルが解となる固有値問題である。なおここで最小の固有値の固有ベクトルは全ての要素が同一になるため破棄される。
【００２６】
従来のＬＬＥのアルゴリズムは高い計算コストが必要になる。はじめのパートでは行列Ｇｐを計算し、その逆行列を計算する。行列Ｇｐの大きさがＫ×Ｋであり、その要素が長さＭのベクトルの内積から計算されるため、行列Ｇｐを全てのデータポイントに対して求める計算コストはＯ（ＮＫ２Ｍ）となる。また行列Ｇｐの逆行列を計算する必要があるため、エッジの重みを求めるためにＯ（ＮＫ３）の計算コストが必要になる。さらにすべてのノードに対して近傍を求めるのにＯ（Ｎ２Ｍ）の計算コストが必要になる。また二つ目のパートではＯ（Ｎ３）の計算コストが必要になる。これはＮ×ＮのカーネルＫの固有値分解を行う必要があるからである。そのため従来のＬＬＥの計算コストはＯ（Ｎ（Ｋ２Ｍ＋Ｋ３）＋Ｎ２Ｍ＋Ｎ３）となる。またメモリコストはカーネルＫの固有値分解を行うためＯ（Ｎ２）となる。結果として大規模なデータに対してＬＬＥを適用するのは現実的ではない。
【００２７】
［第１の実施形態の構成］
図１を用いて、第１の実施形態の構成について説明する。図１は、第１の実施形態に係るＬＬＥ計算装置の構成の一例を示す図である。図１に示すように、ＬＬＥ計算装置１０は、特異値分解計算部１０１、データポイント選択部１０２、近似距離計算部１０３、距離推定部１０４、近傍確定部１０５、エッジ重み計算部１０６、並び替え部１０７、ＬＵ分解計算部１０８及びベクトル計算部１０９を有する。
【００２８】
［第１の実施形態のＬＬＥ］
まず本実施形態におけるＬＬＥの計算方法及び数理的な背景について説明する。その後、ＬＬＥ計算装置１０の各機能部によって行われる処理について説明する。
【００２９】
本実施形態では多次元空間において二つのデータポイントが近くにあればそれらの近傍の大部分は共通するという知見に基づいてＬＬＥを高速に計算する。ここで、式（３）のグラム行列の定義に見られるように、従来のＬＬＥにおいては、データポイントｘｉとｘｊをデータポイントｘｐのＫ近傍とすると、データポイントｘｐのグラム行列の［ｉ，ｊ］はＧｐ［ｉ，ｊ］＝＜ｘｐ−ｘｉ，ｘｐ−ｘｊ＞となる。そのため、ｘｐ＝ｘｑであるようなデータポイントｘｐとｘｑに対して、ｘｐとｘｑがデータポイントｘｉとｘｊを近傍として共有したとしてもそれらのグラム行列の［ｉ，ｊ］成分は同じにならない。
【００３０】
そこで、本実施形態では、高速に処理を行うため、オリジナルのアルゴリズムとは異なる方法でラグランジュの未定乗数法をデータポイントｘｐに対して適用し、行列Ｃｐを計算する。本実施形態では、行列Ｃｐの［ｉ，ｊ］成分はＣｐ［ｉ，ｊ］＝＜ｘｉ，ｘｊ＞となるため、もし同じデータポイントを近傍として共有していれば、異なるデータポイントに対して行列Ｃｐは同じ成分を保持することとなる。そのため高速に行列Ｃｐを計算することが可能になる。
【００３１】
さらにこの手法を用いることで逆行列を高速に計算することが可能になる。具体的には、この手法により行列Ｃｐは行列Ｃｑと同じ成分を有するため、Ｗｏｏｄｂｕｒｙの公式を用いることにより、行列Ｃｑ−１から漸進的に行列Ｃｐ−１を計算することができる。
【００３２】
埋め込みコストを下げるために、本発明では逆ベキ乗法を用いてカーネルＫの固有値ベクトルを高速に計算する（例えば、参考文献１（Brian  Bradie,  A  Friendly  Introduction  to  Numerical  Analysis,  2007,  Pearson.）を参照）。しかし本実施形態ではナイーブにこの手法は用いない。それはこの手法を用いるのにカーネルＫの逆行列を計算する必要があり、それによりＯ（Ｎ３）の計算コストが発生するからである。またＯ（Ｎ２）  のメモリコストが発生するという問題もある。先に述べたとおりカーネルＫはＫ＝（Ｉ−Ｗ）Ｔ（Ｉ−Ｗ）と計算される。ここでＷは近傍のデータポイントからのエッジの重みから計算されるため、カーネルＫは疎なデータ構造となる。しかしカーネルＫが疎であってもその逆行列は密な構造となり、高い計算コストとメモリコストが必要になる。
【００３３】
これらの問題を解決するために、本実施形態では行列Ｉ−ＷのＬＵ分解を計算し、カーネルＫの固有ベクトルを計算する。ここでＬＵ分解とは行列を下三角行列と上三角行列に分解する手法である（例えば、参考文献２（William  H.  Press,  Saul  A.  Teukolsky,  William  T  Vetterling,  Brian  P.  Flannery,  Numerical  Recipes  3rd  Edition,  2007,  Cambridge  University  Press.）を参照）。下三角行列と上三角行列から固有ベクトルを計算できるため、カーネルＫの計算を避けることができる。また下三角行列と上三角行列は疎なデータ構造を持つため、計算コストとメモリコストを抑えることができる。
【００３４】
本実施形態は従来のＬＬＥのアルゴリズムと同じ計算結果となる。これは本実施形態が再構築コストと埋め込みコストを最小化するような計算を行うからである。すなわち本実施形態はオリジナルのアルゴリズムの計算結果を変えることなく、少ない計算コストとメモリコストで次元削減を行うことができる。
【００３５】
［第１の実施形態のエッジの重みの計算］
次にエッジの重みを高速計算するための手法について述べる。本発明ではラグランジュの未定乗数法を用いてデータポイントｘｐの行列Ｃｐを計算する。ラグランジュの未定乗数法は、複数の変数を持つ関数の一定の制約下における定常状態を求めることができる（例えば、参考文献３（Christopher  M.  Bishop,  Pattern  Recognition  and  Machine  Learning,  2010,  Springer.）を参照）。ＬＬＥの場合、変数は再構築コストを最小化するエッジの重みであり、制約は各データポイントへのエッジの重み和が１になるというものである。γをラグランジュの未定乗数とすると、各データポイントに対してラグランジュ関数Ｌは以下のようになる。
【００３６】
【数７】
2018190251.tif 000009
【００３７】
よってｘｉ∈Ｎ［ｘｐ］であるようなデータポイントｘｉのエッジの重みＷ［ｐ，ｉ］に対して関数ＬがＷ［ｐ，ｉ］に対して定常状態になる条件は∂Ｌ／∂Ｗ［ｐ，ｉ］＝０となる。同様にγに対して関数Ｌの条件は∂Ｌ／∂γ＝０となる。そのため関数Ｌが定常状態になる条件は式（１２）のようになる。
【００３８】
【数８】
2018190251.tif 000010
【００３９】
式（１２）においてＣｐは（Ｋ＋１）×（Ｋ＋１）の行列であり、ｗとｐは以下で与えられる長さがＫ＋１の列ベクトルである。
【００４０】
【数９】
2018190251.tif 000011
【００４１】
ここでｘｉ，ｘｊ∈Ｎ［ｘｐ］である。式（１３）から行列Ｃｐの要素はＣｐ［ｉ，ｊ］＝＜ｘｉ，ｘｊ＞と計算することができ、行列Ｃｐの要素はｘｐとは独立に計算することができることがわかる。このためもしデータポイントｘｐとｘｑが同じ近傍を共有するとき、データポイントｘｐの行列Ｃｐの要素はデータポイントｘｑの行列Ｃｑを参照することで高速に計算することができる。もしｄがデータポイントｘｐとｘｑで異なる近傍の数とすると（すなわちｄ＝Ｋ−｜Ｎ［ｘｐ］∩Ｎ［ｘｑ］｜のとき）、本実施形態では行列Ｃｐを求めるためにＯ（ｄＫＭ）の計算コストが必要となる。
【００４２】
［第１の実施形態の逆行列の計算］
次にＷｏｏｄｂｕｒｙの公式を用いて行列Ｃｐの逆行列Ｃｐ−１を高速に求める手法について述べる。式（１４）の通り、ベクトルｗのｉの要素はエッジの重みＷ［ｐ，ｉ］に対応する。そのためＣｐ−１が行列Ｃｐの逆行列である式（１２）からエッジの重みをｗ＝Ｃｐ−１と計算することができる。しかし直接行列Ｃｐから行列Ｃｐ−１を求めるにはＯ（Ｋ３）の計算コストがかかるため、この手法の計算コストは高くなる。
【００４３】
Ｗｏｏｄｂｕｒｙの公式を用いてエッジの重みを高速に計算する手法を述べるが、まずｄ＝１である場合、すなわちデータポイントｘｐがデータポイントｘｑに対して一つだけ異なる近傍を持つ場合（Ｋ−｜Ｎ［ｘｐ］∩Ｎ［ｘｑ］｜＝１）について述べる。ｘｋをデータポイントｘｐにおけるただ一つの異なる近傍のデータポイントとし、ｘｋをデータポイントｘｑのｘｋに対する近傍のデータポイントとする。そのためｘｋ＝Ｎ［ｘｐ］−Ｎ［ｘｐ］∩Ｎ［ｘｑ］でありｘｋ＝Ｎ［ｘｑ］−Ｎ［ｘｐ］∩Ｎ［ｘｑ］となる。
【００４４】
ここで、データポイントに対応する各ベクトルの分散は１であること、すなわち＜ｘｉ，ｘｉ＞＝１であることを仮定する。この仮定は後に外す。行列ＣｐとＣｑのｋ番目の行及び列にデータポイントｘｋとｘｋ´が対応する逆行列Ｃｐ−１を更新するときに、以下の行列ΔＣを用いる。
【００４５】
【数１０】
2018190251.tif 000012
【００４６】
行列ΔＣに対して図２に示す補題１（補助定理１）が成り立つ。図２は、補題１を示す図である。また、補題１は、行列ΔＣは非対称行列ＣｐとＣｑから計算されるが対称行列であり、行列ΔＣのｋの行及び列のみが非ゼロ要素を持つことを示している。この補助定理より行列ΔＣについての図３に示す補題２（補助定理２）が成り立つ。図３は、補題２を示す図である。
【００４７】
ＶＴＶ＝Ｉであり、ΔＣｖ１＝｜｜ｆ｜｜ｖ１かつΔＣｖ２＝−｜｜ｆ｜｜ｖ２であるため、補助定理２は補助定理１とともに、式（１９）は行列ΔＣのランク２の固有値分解に対応していることがわかる。この式を用いることにより行列ＶとＦをＯ（Ｋ）の計算コストで求めることができる。これはベクトルｆの長さがＫ＋１だからである。そのため補助定理２から行列ΔＣの固有値分解をＯ（Ｋ）の計算コストで求めることができる。
【００４８】
Ｗｏｏｄｂｕｒｙの公式を用いることでデータポイントｘｑの逆行列Ｃｑ−１からデータポイントｘｐの逆行列Ｃｐ−１を漸進的に更新することができる。式（１６）からＣｐ＝Ｃｑ＋ΔＣであるため、補助定理２からＣｐ＝Ｃｑ＋ΔＣとなる。そのためＷｏｏｄｂｕｒｙの公式を用いることで逆行列Ｃｐ−１を以下のように計算できる。
【００４９】
【数１１】
2018190251.tif 000013
【００５０】
式（２０）においてＦ−１とＶＴＣｑ−１Ｖは２×２の行列であるため、逆行列（Ｆ−１＋ＶＴＣｑ−１Ｖ）−１を求める計算コストは一定になる。また行列Ｃｑ−１とＶのサイズはそれぞれ（Ｋ＋１）×（Ｋ＋１）と（Ｋ＋１）×２になる。そのため式（２０）からＣｐ−１を求めるにはＯ（Ｋ２）の計算コストが必要となる。
【００５１】
式（２０）から逆行列Ｃｐ−１を高速に求めることができるが、この式は行列Ｃｐにおいて＜ｘｉ，ｘｉ＞＝１が成り立つことを仮定している。ここではこの仮定を用いてまず行列Ｃｐ−１を更新する。そしてＣ´ｐを行列Ｃｐに対応する行列としたとき、行列Ｃ´ｐのために（Ｋ＋１）×（Ｋ＋１）の行列ΔＣを式（２１）のように計算する。
【００５２】
【数１２】
2018190251.tif 000014
【００５３】
行列ΔＣ´から行列Ｃ´ｐとＣｐの違いは［ｋ，ｋ］要素に限られることがわかり、また行列ΔＣ´は＜ｘｉ，ｘｉ＞＝１  ｉｎ  ｍａｔｒｉｘ  Ｃｐという仮定を落とすことに対応することがわかる。そのため、もしその仮定を落とせばＣｐ＝Ｃｑ＋ΔＣ＋ΔＣ´が成り立つことがわかる。もしｅｋをｋ番目の要素以外は０である単位ベクトルだとすると、行列ΔＣ´を式（２２）のように求めることができる。
【００５４】
【数１３】
2018190251.tif 000015
【００５５】
そのためＷｏｏｄｂｕｒｙの公式を用いることで式（２０）の後で行列Ｃ´ｐの逆行列（Ｃ´ｐ）−１を式（２３）のように計算することができる。
【００５６】
【数１４】
2018190251.tif 000016
【００５７】
式（２０）と同様、式（２３）から逆行列をＯ（Ｋ２）の計算コストで高速に計算することができる。式（２０）と（２３）において、ｄ＝１が成り立つこと、つまりデータポイントｘｑに対してデータポイントｘｐは一つだけ異なる近傍のデータポイントを有することを仮定した。しかしｄ＞１の場合であっても、繰り返し上記処理をｘｉ∈｛ｘｊ｜ｘｊ∈Ｎ［ｘｐ］  ａｎｄ  ｘｊ／∈Ｎ［ｘｑ］｝であるようなデータポイント一つ一つに繰り返し適用することによって逆行列を高速に求めることができる。なお、ここでは、ｘ／∈Ｎは、ｘが集合Ｎに含まれないことを示すものとする。
【００５８】
［第１の実施形態の近傍のデータポイントの計算］
次に各データポイントに対して近傍のデータポイントを高速に計算する手法について述べる。本発明ではＳＶＤを用いて近似の距離を計算し、近傍を計算する。ＳＶＤは距離の下限値を求めることができるため、正確に近傍を求めることができる。もし〜ｘｐ＝（〜ｘｐ［１］，〜ｘｐ［２］，．．．，〜ｘｐ［ｓ］）がデータポイントｘｐ＝（ｘｐ［１］，ｘｐ［２］，．．．，ｘｐ［Ｍ］）に対するランクがｓのＳＶＤによる次元削減とし、Ｅ［ｘｐ，ｘｑ］をデータポイントｘｐとｘｑのユークリッド距離とすると、ＳＶＤが距離の下限値を与えるとはＥ［ｘｐ，ｘｑ］≧Ｅ［〜ｘｐ，〜ｘｑ］が成り立つと言うことである。
【００５９】
近傍を求めるのに既存の研究ではまず近似の距離を計算してから、もしデータポイントが近似の距離により近傍になる可能性があると判定されたらユークリッド距離を計算する（例えば、参考文献４（Dennis  Shasha,  High  Performance  Discovery  in  Time  Series:  Techniques  and  Case  Studies,  2004,  Springer-Verlag  New  York  Inc..）を参照）。なお各データポイントのＳＶＤは既存の手法を用いてＯ（ＮＭｌｏｇｓ）の計算コストで求めることができる（例えば、参考文献５（Nathan  Halko,  Per-Gunnar  Martinsson,  Joel  A.  Tropp,  Finding  Structure  with  Randomness:Probabilistic  Algorithms  for  Constructing  Approｘimate  Matriｘ  Decompositions,  2011,  SIAM  Review.）を参照）。
【００６０】
近似の距離の精度を向上するために、ユークリッド距離を近傍グラフの構造を用いて推定する。データポイントｘｑはデータポイントｘｒの近傍とする。このときＳＶＤを用いることによってＥ［〜ｘｑ，〜ｘｒ］とＥ［ｘｑ，ｘｒ］が成り立つ。またもしデータポイントｘｐの近傍を求める過程でデータポイントｘｒのユークリッド距離を計算すれば、Ｅ［〜ｘｐ，〜ｘｒ］  ａｎｄ  Ｅ［ｘｐ，ｘｒ］が成り立つ。Ｅ＿［ｘｐ，ｘｑ］をデータポイントｘｐからｘｑの推定距離とすると、本実施形態では推定距離を式（２４）のように計算する。ここで・＿は・の直下に＿が記された記号を示すものとする。
【００６１】
【数１５】
2018190251.tif 000017
【００６２】
ここでｕｒ［ｘｐ］はｕｒ［ｘｐ］＝（〜ｘｐ［ｓ＋１］−〜ｘｒ［ｓ＋１］，〜ｘｐ［ｓ＋２］−〜ｘｒ［ｓ＋２］，．．．，〜ｘｐ［Ｍ］−〜ｘｒ［Ｍ］）であるような長さがＭ−ｓのベクトルである。推定距離について図４に示す補題３（補助定理３）が成り立つ。図４は、補題３を示す図である。
【００６３】
補助定理３は推定距離Ｅ［ｘｐ，ｘｑ］のＥ［ｘｐ，ｘｑ］に対する誤差はＥ［〜ｘｐ，〜ｘｑ］より小さいため、Ｅ［ｘｐ，ｘｑ］も用いることでＥ［〜ｘｐ，〜ｘｑ］の近似の精度を向上できることを示している。また補助定理３のＥ［ｘｐ，ｘｑ］≧Ｅ［ｘｐ，ｘｑ］は下限値の性質を満たす、すなわち正確に近傍を求めることができることを示している。
【００６４】
この手法において式（２４）で用いられるベクトルｕｒ［ｘｐ］とｕｒ［ｘｑ］のＬ２ノルムを式（２５）及び式（２６）のように計算する。
【００６５】
【数１６】
2018190251.tif 000018
【００６６】
ここでもしデータポイントｘｐの近傍としてｘｒをチェックするするためにデータポイントｘｐとｘｒのユークリッド距離を計算し、データポイントｘｑがｘｒの近傍であれば、式（２５）と（２６）を用いることによりＥ［ｘｐ，ｘｑ］を求めるための計算コストはＯ（１）となる。これはこの場合、Ｅ［ｘｐ，ｘｒ］とＥ［〜ｘｐ，〜ｘｒ］とＥ［ｘｑ，ｘｒ］とＥ［〜ｘｑ，〜ｘｒ］をＯ（１）の計算コストで求めることができるからである。
【００６７】
次に固有ベクトルを高速に計算する方法について述べる。従来のＬＬＥのアルゴリズムはデータセットＸの低次元への埋め込みＹを埋め込みコストへの最適解により計算する。エッジの重み行列Ｗに対してＬＬＥカーネルはＫ＝（Ｉ−Ｗ）Ｔ（Ｉ−Ｗ）と計算されるが、この最適解はカーネルＫの底の固有ベクトルに対応する。そのためもしλｉが｜λ１｜≧｜λ２｜≧．．．≧｜λＮ｜であるようなカーネルＫのｉ番目の固有値であり、またｚｉをそれに対応する固有ベクトルとすると、データセットＸに対するｍ次元の埋め込みＹはｍ＋１個の底の固有ベクトルからＹ＝［ｚＮ−ｍ，ｚＮ−ｍ＋１，．．．，ｚＮ−１］として計算することができる。
【００６８】
高速に固有ベクトルを計算するために、本実施形態では逆ベキ乗法を用いる。逆ベキ乗法は最も底の固有ベクトルを、固有値を計算する行列の逆行列から計算する方法で、具体的には逆行列Ｋ−１から固有ベクトルｚＮを計算することができる。そのためａ０を長さがＮの任意の列ベクトルとすると、底の固有値λＮａとその固有ベクトルｚＮはベキ乗法をａτ＝Ｋ−１ａτ−１と適用して計算することができる。しかし逆行列Ｋ−１は密な構造となるため、この手法に必要な計算コストはＯ（Ｎ３）となりメモリコストはＯ（Ｎ２）となる。
【００６９】
直接カーネルＫの逆行列を避けるために、本発明ではＩ−ＷのＬＵ分解を計算する。具体的にはＬとＵを下三角行列及び上三角行列としたとき、行列ＬとＵをＬＵ＝Ｉ−Ｗとして計算する。その結果式（２７）のようになる。
【００７０】
【数１７】
2018190251.tif 000019
【００７１】
本実施形態では行列ＬとＵにベキ乗法を用いてａτ−１からａτを式（２８）のように計算する。
【００７２】
【数１８】
2018190251.tif 000020
【００７３】
式（２８）においてｂとｂ´とｂ´´とａτを計算するために下三角行列と上三角行列に前進代入と後進代入を適用する（例えば、参考文献２を参照）。式（２８）は下三角行列と上三角行列からベクトルａτをベクトルａτ−１からａτ−１＝ＵＴＬＴＬＵａτという形で求められることを示している。式（２７）からＫ＝ＵＴＬＴＬＵが成り立つため、カーネルＫの逆行列を直接用いることなくａτ＝Ｋ−１ａτ−１を計算することができる。式（２８）から式（２９）のようにベクトルａτが収束するまでベキ乗法を適用することにより、底の固有値λＮと固有ベクトルｚＮを計算することができる。
【００７４】
【数１９】
2018190251.tif 000021
【００７５】
カーネルＫの下三角行列と上三角行列は疎な構造を持つため、本発明は高速にベクトルａτを式（２８）から計算することができる。次にその他の底の固有値λＮ−ｉと固有ベクトルｚＮ−ｉ（ｉ＝１，２，．．．，ｍ）を計算する方法を述べる。このような固有値と固有ベクトルを求める方法としてホテリング法が有名である（例えば、参考文献１を参照）。ホテリング法の基本的アイデアは固有値を計算する行列を最大固有値が０としその他の固有値が変わらないようにシフトするというものである。その結果２番目に大きな固有値がシフトした行列において最大の固有値になる。具体的にはベキ乗法を用いて式（３０）の行列Ｈｉの固有値と固有ベクトルを計算することでλＮ−ｉとｚＮ−ｉを計算する。
【００７６】
【数２０】
2018190251.tif 000022
【００７７】
ここでｈｉ，τを長さがＮの列ベクトルとしたとき、固有値λＮ−ｉと固有ベクトルｚＮ−ｉはｈｉ，τ＝Ｈｉｈｉ，τ−１とベキ乗法を適することにより求めることができる。しかしこの方法は行列Ｈｉが密な構造を持つため、Ｏ（Ｎ３）の計算コストとＯ（Ｎ２）のメモリコストが必要になるという問題がある。
【００７８】
行列Ｈｉを直接計算することを避けるために、本発明ではベクトルｈｉ，τをベクトルｈｉ，τ−１から式（３１）のように計算する。
【００７９】
【数２１】
2018190251.tif 000023
【００８０】
ここで式（２８）と同様に以下のようにベクトルａτは前進代入と後進代入を用いてベクトルｈｉ，τ−１から式（３２）のように計算する。
【００８１】
【数２２】
2018190251.tif 000024
【００８２】
式（３１）において長さがＮのベクトルｚＮ−ｊｚＮ−ｊＴｈｉ，τ−１／λＮ−ｊはベクトルを右から左へ計算することでＯ（Ｎ）の計算コストで求める。Ｋ−１＝（ＵＴＬＴＬＵ）−１となるため、式（３０）、（３１）、（３２）から式（３３）が成り立つ。
【００８３】
【数２３】
2018190251.tif 000025
【００８４】
式（３３）は式（３１）と（３２）からベクトルｈｉ，τをホテリング法の行列Ｈｉを用いて計算できることを示している。本実施形態ではカーネルＫの固有値λＮ−ｉと固有ベクトルｚＮ−ｉ（１≦ｉ≦ｍ）を式（３４）に基づきベキ乗法によりベクトルｈｉ，τが収束するまで計算することにより求める。
【００８５】
【数２４】
2018190251.tif 000026
【００８６】
直接行列Ｈｉを計算しないため、本発明は高速に固有値と固有ベクトルを計算することができる。
【００８７】
次に逆ベキ乗法において用いる行列Ｉ−ＷのＬＵ分解を高速に計算する方法について述べる。行列ＬとＵの要素はもし行列Ｗ´がＷ´＝Ｉ−Ｗと与えられるときクラウトのアルゴリズム（例えば、参考文献２を参照）により式（３５）及び式（３６）のように計算できる。
【００８８】
【数２５】
2018190251.tif 000027
【００８９】
式（３５）と（３６）は行列ＬとＵの列は左から右へ、また各列において要素は上から下へ計算できることがわかる。そのため行列ＬとＵの要素は対応する行列Ｗ´、Ｌ、Ｕの上と左の要素から計算することができる。
【００９０】
本実施形態では行列ＬとＵの非零要素の数を減らしてＬＵ分解を高速に計算する。これはもしＬ［ｉ，ｋ］＝０又はＵ［ｋ，ｊ］＝０となれば式（３５）と（３６）から効果的にＬ［ｉ，ｋ］Ｕ［ｋ，ｊ］の計算を省くことができるからである。本実施形態では行列ＬとＵの上と左の要素はもし対応する行列Ｗ´の要素が零であれば零になることと、行列Ｗ´の上と左の要素はＷ´＝Ｉ−Ｗであるためもし対応する行列Ｗの左と上の要素が零になれば零になることを利用する。行列Ｗはエッジの重みから計算されるため、行列Ｗのデータポイントを並び変えることで上と左の要素が疎になればＬＵ分解を高速に行うことができる。
【００９１】
本実施形態ではデータポイントを次数の昇順で並び変えることでＬＵ分解を高速に行う。ここで次数とはデータポイントから出ているエッジの本数である。並び変えることで行列Ｗの上と左の要素を疎にできるため、行列Ｉ−ＷのＬＵ  分解を高速行うことができる。データポイントを並び変えるのには分布数え上げソートを用いることでＯ（Ｎ）の計算コストで並び替えを行うことができる（例えば、参考文献６（Thomas  H.  Cormen,  Charles  E.  Leiserson,  Ronald  L.  Rivest,  Clifford  Stein,  Introduction  to  Algorithms,  2009,  The  MIT  Press.）を参照）。ｎを三角行列における平均的な各行ベクトルの非零要素の数とすると、行列ＬとＵを保持するメモリコストはＯ（Ｎｎ）となる。また式（３５）と（３６）から行列Ｉ−ＷのＬＵ分解はＯ（Ｎｎ２）の計算コストで行うことができる。
【００９２】
［第１の実施形態のアルゴリズム］
本実施形態のアルゴリズムを、ＬＬＥ計算装置１０の各機能部の処理とともに、図５を用いて説明する。図５は、第１の実施形態に係るＬＬＥ計算処理のアルゴリズムの一例を示す図である。
【００９３】
アルゴリズム１５において、Ｓはエッジの重みを計算するために選択したデータポイントの集合であり、Ｄは選択されたデータポイントの近傍のデータポイントの集合であり、Ｃ［ｘｐ］は近傍グラフにおいてデータポイントｘｐに直接つながっているデータポイントの集合である。アルゴリズム１５において近傍を高速に計算するためにＳＶＤを計算するが、もしｓをｓ＝２ｌであるようなＳＶＤのランクとしたとき、本実施形態ではランクを２１，２２，．．．，２ｌと徐々に増やしていきながら近傍を求めていく。計算コストの観点では粗い近似であれば高速に近似距離を計算することができるが近似の精度は高くない。一方、細かい精度であれば高い精度で近似距離を計算できるが、高い計算コストがかかってしまう。そのため、本実施形態では近傍にならないデータポイントは粗い近似で枝刈りし、近傍になる可能性が高いデータポイントは細かい近似で近似距離を計算する。結果として本発明は各データポイントの距離に応じて適切に近似の精度を設定して近傍を求めることができる。
【００９４】
アルゴリズム１５に示すように、まず、行列Ｘ、低次元数ｍ、ＳＶＤのランクｓ（＝２ｌ）が入力される。そして、特異値分解計算部１０１は、複数のデータポイントを有する多次元行列の特異値分解（ＳＶＤ）を計算する。具体的には、特異値分解計算部１０１はＳを初期化し、行列ＸのＳＶＤを計算する（１−２行目）。
【００９５】
データポイント選択部１０２は、複数のデータポイントから計算対象のデータポイントを選択する。具体的には、データポイント選択部１０２は、他のデータポイントへ多くのエッジを持つデータポイントを選択する（４行目）。これは本実施形態では近傍グラフを用いて他のデータポイントへの距離を推定するからである。そして、データポイント選択部１０２は、Ｄとθを初期化し（５−６行目）、Ｋ個の距離が∞になるデータポイントを追加することでＮ［ｘｐ］を初期化する（７行目）。
【００９６】
近似距離計算部１０３は、複数のデータポイントのそれぞれについて、計算対象のデータポイントとの間の近似の距離を、特異値分解に基づいて計算する。このとき、近似距離計算部１０３は、近傍を計算するためにランクを２１から２ｌに徐々に増やしていくが、もし近似の距離がθより大きくなればデータポイントを枝刈りする（９−１３行目）。
【００９７】
また、距離推定部１０４は、近似の距離が所定値以下であるデータポイントのそれぞれと、計算対象のデータポイントとの間のユークリッド距離の推定値を計算する。つまり、距離推定部１０４は、ＳＶＤを用いてもデータポイントを枝刈りできなかった場合はグラフ構造を用いて距離を推定する（１４−１８行目）。
【００９８】
近傍確定部１０５は、ユークリッド距離の推定値が所定値以下であるデータポイントのそれぞれと、計算対象のデータポイントとの間のユークリッド距離を計算し、当該計算したユークリッド距離が所定値以下であるデータポイントを、計算対象のデータポイントの近傍のデータポイントに確定する。近傍確定部１０５は、データポイントを、選択されたデータポイントの近傍に確定した場合は、Ｎ［ｘｐ］をそのデータポイントの距離を用いて更新する（１９−２５行目）。
【００９９】
そして、近傍確定部１０５は、最も多くの近傍を持つデータポイントを見つける（２６行目）。これはエッジ重み計算部１０６が、エッジの重みを他のデータポイントの近傍を用いて計算するためである。
【０１００】
エッジ重み計算部１０６は、計算対象のデータポイントと計算対象のデータポイントの近傍のデータポイントのそれぞれとの間のエッジの重みを計算する。エッジ重み計算部１０６は、もし共有する近傍がある場合はＷｏｏｄｂｕｒｙの公式を用いてエッジの重みを計算し（２７−２８行目）、そうでない場合はグラム行列を用いてエッジの重みを計算する（２９−３０行目）。このように、エッジ重み計算部１０６は、Ｗｏｏｄｂｕｒｙの公式を用いて、計算済みデータポイントの近傍データポイントに関する逆行列から、計算対象のデータポイントの近傍のデータポイントに関する逆行列を計算することができる。
【０１０１】
つまり、エッジ重み計算部１０６は、エッジの重みが計算済みである計算済みデータポイントの近傍データポイントに、計算対象のデータポイントの近傍のデータポイントのうちの少なくとも一部が含まれる場合、計算済みデータポイントについてのエッジの重みの計算結果に基づいて、計算対象のデータポイントと計算対象のデータポイントの近傍のデータポイントのそれぞれとの間のエッジの重みを計算することができる。
【０１０２】
並び替え部１０７は、複数のデータポイントについてのエッジの重みを表す行列を、次数の昇順に並び替える（３２−３３行目）。また、ＬＵ分解計算部１０８は、複数のデータポイントの隣接行列のＬＵ分解を計算する。
【０１０３】
ベクトル計算部１０９は、エッジの重みに基づいて、複数のデータポイントについてのＬＬＥカーネルの固有ベクトルを計算する。例えば、ベクトル計算部１０９は、ベクトル計算部１０９は、ＬＵ分解計算部１０８によって計算されたＬＵ分解に基づいて、逆ベキ乗法を用いてＬＬＥカーネルの固有ベクトル及び固有値を計算する（３４−３６行目）。
【０１０４】
なお、２８−３０行目で、エッジ重み計算部１０６は、常に２９行目を実行するようにしてもよい。この場合、従来のＬＬＥと同様に、エッジ重み計算部１０６は、グラム行列を用いてエッジの重みの計算を行う。また、並び替え部１０７及びＬＵ分解計算部１０８による処理が実行されないようにしてもよい。この場合、ベクトル計算部１０９は、従来のＬＬＥと同様の方法で固有ベクトル及び固有値を計算する。また、また、図５のアルゴリズム１５に対しては、図６に示す性質（定理１〜３）が成り立つ。図６は、定理１、定理２及び定理３を示す図である。
【０１０５】
［第１の実施形態の処理］
図７を用いて、ＬＬＥ計算装置１０の処理の流れについて説明する。図７は、第１の実施形態に係るＬＬＥ計算装置の処理の流れを示すフローチャートである。まず、ＬＬＥ計算装置１０には、行列Ｘ、次元数ｍ及びランクｓが入力される（ステップＳ１０１）。次に、特異値分解計算部１０１はＳを空集合に初期化する（ステップＳ１０２）。そして、特異値分解計算部１０１は行列ＸのランクｓのＳＶＤを計算する（ステップＳ１０３）。
【０１０６】
データポイント選択部１０２は、ｉを１からＮまで増やしながら複数のデータポイントから計算対象のデータポイントを選択する（ステップＳ１０４、Ｓ１０５、Ｓ１３２）。そして、データポイント選択部１０２は、Ｄとθを初期化し（ステップＳ１０６、Ｓ１０７）、さらに、Ｋ個の距離が∞になるデータポイントを追加することでＮ［ｘｐ］を初期化する（ステップＳ１０８）。
【０１０７】
近似距離計算部１０３は、複数のデータポイントのそれぞれについて、計算対象のデータポイントとの間の近似の距離を、特異値分解に基づいて計算する（ステップＳ１０９、Ｓ１１０、Ｓ１１１）。このとき、近似距離計算部１０３は、近傍を計算するためにランクを２１から２ｌに徐々に増やしていくが、もし近似の距離がθより大きくなれば（ステップＳ１１２、Ｙｅｓ）データポイントを枝刈りする（ステップＳ１１３）。また、近似の距離がθより大きくない場合（ステップＳ１１２、Ｎｏ）近似距離計算部１０３は次の処理に進む（ステップＳ１１４）。
【０１０８】
ここで、枝刈りされていないデータポイントが残っていない場合（ステップＳ１１５、Ｎｏ）、ＬＬＥ計算装置１０は、ステップＳ１２８へ進む。一方、枝刈りされていないデータポイントが残っている場合（ステップＳ１１５、Ｙｅｓ）、距離推定部１０４は、近似の距離が所定値以下であるデータポイントのそれぞれと、計算対象のデータポイントとの間のユークリッド距離の推定値を計算する（ステップＳ１１６、Ｓ１１７）。距離の推定値がθより大きくなれば（ステップＳ１１８、Ｙｅｓ）、距離推定部１０４はデータポイントを枝刈りする（ステップＳ１１９）。また、近似の距離がθより大きくない場合（ステップＳ１１８、Ｎｏ）距離推定部１０４は次の処理に進む（ステップＳ１２０）。
【０１０９】
ここで、枝刈りされていないデータポイントが残っていない場合（ステップＳ１２１、Ｎｏ）、ＬＬＥ計算装置１０は、ステップＳ１２８へ進む。一方、枝刈りされていないデータポイントが残っている場合（ステップＳ１２１、Ｙｅｓ）、近傍確定部１０５は、ユークリッド距離の推定値が所定値以下であるデータポイントのそれぞれと、計算対象のデータポイントとの間のユークリッド距離を計算し（ステップＳ１２２）、当該計算したユークリッド距離が所定値以下であるデータポイント（ステップＳ１２３、Ｙｅｓ）を、計算対象のデータポイントの近傍のデータポイントに確定し、Ｎ［ｘｐ］をそのデータポイントの距離を用いて更新する（ステップＳ１２４、Ｓ１２５、Ｓ１２６）。また、近傍確定部１０５は、θを更新する（ステップＳ１２７）。そして、近傍確定部１０５は、最も多くの近傍を持つデータポイントを見つける（ステップＳ１２８）。
【０１１０】
エッジ重み計算部１０６は、共有する近傍がある場合（ステップＳ１２９、Ｙｅｓ）はＷｏｏｄｂｕｒｙの公式を用いてエッジの重みを計算し（ステップＳ１３０）、そうでない場合（ステップＳ１２９、Ｎｏ）はグラム行列を用いてエッジの重みを計算する（ステップＳ１３１）。
【０１１１】
並び替え部１０７は、複数のデータポイントについてのエッジの重みを表す行列を、次数の昇順に並び替える（ステップＳ１３３）。また、ＬＵ分解計算部１０８は、複数のデータポイントの隣接行列のＬＵ分解を計算する（ステップＳ１３４）。
【０１１２】
ベクトル計算部１０９は、エッジの重みに基づいて、複数のデータポイントについてのＬＬＥカーネルの固有ベクトルを計算する。例えば、ベクトル計算部１０９は、ＬＵ分解計算部１０８によって計算されたＬＵ分解に基づいて、逆ベキ乗法を用いてＬＬＥカーネルの固有ベクトル及び固有値を計算する（ステップＳ１３５、Ｓ１３６、Ｓ１３７、Ｓ１３８）。そして、ＬＬＥ計算装置１０は、次元削減された行列Ｙを出力する（ステップＳ１３９）。
【０１１３】
［第１の実施形態の効果］
特異値分解計算部１０１は、複数のデータポイントを有する多次元行列の特異値分解を計算する。また、データポイント選択部１０２は、複数のデータポイントから計算対象のデータポイントを選択する。また、近似距離計算部１０３は、複数のデータポイントのそれぞれについて、計算対象のデータポイントとの間の近似の距離を、特異値分解に基づいて計算する。また、距離推定部１０４は、近似の距離が所定値以下であるデータポイントのそれぞれと、計算対象のデータポイントとの間のユークリッド距離の推定値を計算する。また、近傍確定部１０５は、ユークリッド距離の推定値が所定値以下であるデータポイントのそれぞれと、計算対象のデータポイントとの間のユークリッド距離を計算し、当該計算したユークリッド距離が所定値以下であるデータポイントを、計算対象のデータポイントの近傍のデータポイントに確定する。また、エッジ重み計算部１０６は、計算対象のデータポイントと計算対象のデータポイントの近傍のデータポイントのそれぞれとの間のエッジの重みを計算する。また、ベクトル計算部１０９は、エッジの重みに基づいて、複数のデータポイントについてのＬＬＥカーネルの固有ベクトルを計算する。このように、本実施形態では、近傍確定部１０５によってユークリッドの距離が実際に計算される対象を、近似の距離及びユークリッド距離の推定値を用いてあらかじめ削減している。このため、本実施形態によれば、ＬＬＥによる次元削減の際の計算コスト及びメモリコストを低減させることができる。
【０１１４】
エッジ重み計算部１０６は、エッジの重みが計算済みである計算済みデータポイントの近傍データポイントに、計算対象のデータポイントの近傍のデータポイントのうちの少なくとも一部が含まれる場合、計算済みデータポイントについてのエッジの重みの計算結果に基づいて、計算対象のデータポイントと計算対象のデータポイントの近傍のデータポイントのそれぞれとの間のエッジの重みを計算することができる。このように、同じ近傍を持つデータポイントについては、計算済みの結果を利用してすることで、エッジの重みの計算に要するコストを低減させることができる。
【０１１５】
エッジ重み計算部１０６は、Ｗｏｏｄｂｕｒｙの公式を用いて、計算済みデータポイントの近傍データポイントに関する逆行列から、計算対象のデータポイントの近傍のデータポイントに関する逆行列を計算することができる。このように、Ｗｏｏｄｂｕｒｙの公式を用いることで、特に逆行列の計算に要するコストを低減させることができる。
【０１１６】
並び替え部１０７は、複数のデータポイントについてのエッジの重みを表す行列を、次数の昇順に並び替える。また、ＬＵ分解計算部１０８は、複数のデータポイントの隣接行列のＬＵ分解を計算する。このとき、ベクトル計算部１０９は、ＬＵ分解計算部１０８によって計算されたＬＵ分解に基づいてＬＬＥカーネルの固有ベクトルを計算する。このように、ＬＵ分解を行うことで、固有ベクトルの計算に要するコストを低減させることができる。
【０１１７】
［システム構成等］
また、図示した各装置の各構成要素は機能概念的なものであり、必ずしも物理的に図示のように構成されていることを要しない。すなわち、各装置の分散・統合の具体的形態は図示のものに限られず、その全部又は一部を、各種の負荷や使用状況等に応じて、任意の単位で機能的又は物理的に分散・統合して構成することができる。さらに、各装置にて行われる各処理機能は、その全部又は任意の一部が、ＣＰＵ（Central  Processing  Unit）及び当該ＣＰＵにて解析実行されるプログラムにて実現され、あるいは、ワイヤードロジックによるハードウェアとして実現され得る。
【０１１８】
また、本実施形態において説明した各処理のうち、自動的に行われるものとして説明した処理の全部又は一部を手動的に行うこともでき、あるいは、手動的に行われるものとして説明した処理の全部又は一部を公知の方法で自動的に行うこともできる。この他、上記文書中や図面中で示した処理手順、制御手順、具体的名称、各種のデータやパラメータを含む情報については、特記する場合を除いて任意に変更することができる。
【０１１９】
［プログラム］
一実施形態として、ＬＬＥ計算装置１０は、パッケージソフトウェアやオンラインソフトウェアとして上記のＬＬＥの計算を実行するＬＬＥ計算プログラムを所望のコンピュータにインストールさせることによって実装できる。例えば、上記のＬＬＥ計算プログラムを情報処理装置に実行させることにより、情報処理装置をＬＬＥ計算装置１０として機能させることができる。ここで言う情報処理装置には、デスクトップ型又はノート型のパーソナルコンピュータが含まれる。また、その他にも、情報処理装置にはスマートフォン、携帯電話機やＰＨＳ（Personal  Handyphone  System）等の移動体通信端末、さらには、ＰＤＡ（Personal  Digital  Assistant）等のスレート端末等がその範疇に含まれる。
【０１２０】
また、ＬＬＥ計算装置１０は、ユーザが使用する端末装置をクライアントとし、当該クライアントに上記のＬＬＥの計算に関するサービスを提供するＬＬＥ計算サーバ装置として実装することもできる。例えば、ＬＬＥ計算サーバ装置は、多次元行列を入力とし、次元削減した行列を出力とするＬＬＥ計算サービスを提供するサーバ装置として実装される。この場合、ＬＬＥ計算サーバ装置は、Ｗｅｂサーバとして実装することとしてもよいし、アウトソーシングによって上記のＬＬＥの計算に関するサービスを提供するクラウドとして実装することとしてもかまわない。
【０１２１】
図８は、ＬＬＥ計算プログラムを実行するコンピュータの一例を示す図である。コンピュータ１０００は、例えば、メモリ１０１０、ＣＰＵ１０２０を有する。また、コンピュータ１０００は、ハードディスクドライブインタフェース１０３０、ディスクドライブインタフェース１０４０、シリアルポートインタフェース１０５０、ビデオアダプタ１０６０、ネットワークインタフェース１０７０を有する。これらの各部は、バス１０８０によって接続される。
【０１２２】
メモリ１０１０は、ＲＯＭ（Read  Only  Memory）１０１１及びＲＡＭ（Random  access  memory）１０１２を含む。ＲＯＭ１０１１は、例えば、ＢＩＯＳ（Basic  Input  Output  System）等のブートプログラムを記憶する。ハードディスクドライブインタフェース１０３０は、ハードディスクドライブ１０９０に接続される。ディスクドライブインタフェース１０４０は、ディスクドライブ１１００に接続される。例えば磁気ディスクや光ディスク等の着脱可能な記憶媒体が、ディスクドライブ１１００に挿入される。シリアルポートインタフェース１０５０は、例えばマウス１１１０、キーボード１１２０に接続される。ビデオアダプタ１０６０は、例えばディスプレイ１１３０に接続される。
【０１２３】
ハードディスクドライブ１０９０は、例えば、ＯＳ（Operating  System）１０９１、アプリケーションプログラム１０９２、プログラムモジュール１０９３、プログラムデータ１０９４を記憶する。すなわち、ＬＬＥ計算装置１０の各処理を規定するプログラムは、コンピュータにより実行可能なコードが記述されたプログラムモジュール１０９３として実装される。プログラムモジュール１０９３は、例えばハードディスクドライブ１０９０に記憶される。例えば、ＬＬＥ計算装置１０における機能構成と同様の処理を実行するためのプログラムモジュール１０９３が、ハードディスクドライブ１０９０に記憶される。なお、ハードディスクドライブ１０９０は、ＳＳＤにより代替されてもよい。
【０１２４】
また、上述した実施形態の処理で用いられる設定データは、プログラムデータ１０９４として、例えばメモリ１０１０やハードディスクドライブ１０９０に記憶される。そして、ＣＰＵ１０２０が、メモリ１０１０やハードディスクドライブ１０９０に記憶されたプログラムモジュール１０９３やプログラムデータ１０９４を必要に応じてＲＡＭ１０１２に読み出して実行する。
【０１２５】
なお、プログラムモジュール１０９３やプログラムデータ１０９４は、ハードディスクドライブ１０９０に記憶される場合に限らず、例えば着脱可能な記憶媒体に記憶され、ディスクドライブ１１００等を介してＣＰＵ１０２０によって読み出されてもよい。あるいは、プログラムモジュール１０９３及びプログラムデータ１０９４は、ネットワーク（ＬＡＮ（Local  Area  Network）、ＷＡＮ（Wide  Area  Network）等）を介して接続された他のコンピュータに記憶されてもよい。そして、プログラムモジュール１０９３及びプログラムデータ１０９４は、他のコンピュータから、ネットワークインタフェース１０７０を介してＣＰＵ１０２０によって読み出されてもよい。
【符号の説明】
【０１２６】
１０  ＬＬＥ計算装置
１０１  特異値分解計算部
１０２  データポイント選択部
１０３  近似距離計算部
１０４  距離推定部
１０５  近傍確定部
１０６  エッジ重み計算部
１０７  並び替え部
１０８  ＬＵ分解計算部
１０９  ベクトル計算部
【図１】
2018190251.tif 000028
【図２】
2018190251.tif 000029
【図３】
2018190251.tif 000030
【図４】
2018190251.tif 000031
【図５】
2018190251.tif 000032
【図６】
2018190251.tif 000033
【図７】
2018190251.tif 000034
【図８】
2018190251.tif 000035