明細書 :行列の高速高精度特異値分解方法、プログラムおよび装置

発行国	日本国特許庁(JP)
公報種別	特許公報(B2)
特許番号	特許第4325877号 (P4325877)
登録日	平成21年6月19日(2009.6.19)
発行日	平成21年9月2日(2009.9.2)
発明の名称または考案の名称	行列の高速高精度特異値分解方法、プログラムおよび装置
国際特許分類	G06F  17/16        (2006.01) G06T   1/00        (2006.01) G06F  17/30        (2006.01)
ＦＩ	G06F 17/16 K G06T 1/00 315 G06F 17/30 220Z G06F 17/30 350C
請求項の数または発明の数	6
全頁数	39
出願番号	特願2006-514122 (P2006-514122)
出願日	平成17年6月1日(2005.6.1)
国際出願番号	PCT/JP2005/010084
国際公開番号	WO2005/119507
国際公開日	平成17年12月15日(2005.12.15)
優先権出願番号	2004166437
優先日	平成16年6月3日(2004.6.3)
優先権主張国	日本国(JP)
審査請求日	平成18年12月1日(2006.12.1)
特許権者または実用新案権者	【識別番号】503360115 【氏名又は名称】独立行政法人科学技術振興機構
発明者または考案者	【氏名】中村 佳正 【氏名】岩▲崎▼ 雅史 【氏名】阪野 真也
個別代理人の代理人	【識別番号】100078282、【弁理士】、【氏名又は名称】山本 秀策 【識別番号】100062409、【弁理士】、【氏名又は名称】安村 高明 【識別番号】100113413、【弁理士】、【氏名又は名称】森下 夏樹
審査官	【審査官】鳥居 稔
参考文献・文献	特開平１０－０３１７４７（ＪＰ，Ａ） 特開２００２－２０２９８３（ＪＰ，Ａ） Nakamura, Yoshimasa，A new approach to numerical algorithms in terms of integrable systems，Informatics Research for Development of Knowledge Society Infrastructure, 2004. ICKS 2004. International Conference on ，２００４年 ３月 １日，194 - 205 Masashi Iwasaki, Yoshimasa Nakamura，An application of the discrete Lotka?Volterra system with variable step-size to singular value computation，2004 Inverse Problems 20 553-563，２００４年 ２月２７日，553-563 誉田 太朗 Taro Konda, 中村 佳正 Yoshimasa Nakamura ，ロトカ・ボルテラ系による特異値計算アルゴリズムの並列化 A Parallelization of Singular Value Computation Algorithm by the Lotka-Volterra System，情報処理学会研究報告 Vol.2004 No.128 IPSJ SIG Technical Reports，日本，社団法人情報処理学会 Information Processing Society of Japan，２００４年１２月１７日，第2004巻，13-18 中村 佳正 Yoshimasa NAKAMURA，高機能化を目指した非線形システム制御 Nonlinear Engineering and Control，システム／制御／情報 第43巻 第11号 SYSTEMS,CONTROL AND INFORMATION，日本，システム制御情報学会 The Institute of Systems,Control and Information Engineers，１９９９年１１月１５日，第43巻，8-15 中村 佳正 Yoshimasa NAKAMURA，非線形現象の不思議 ソリトン理論と数値計算法，電子情報通信学会誌 第８０巻 第１１号 THE JOURNAL OF THE INSTITUTE OF ELECTRONICS,INFORMATION AND COMMUNICATION ENGINEERS，日本，社団法人電子情報通信学会 THE INSTITUTE OF ELECTRONICS,INFORMATION AND COMMUNICATION ENGINEERS，１９９７年１１月２５日，第80巻，1143-1146
調査した分野	G06F 17/16 G06F 17/30 G06T 1/00 IEEE Xplore JSTPlus(JDreamII)
特許請求の範囲	【請求項7】 物体の複数の２次元画像から３次元画像を復元する方法であって、 ２次元画像ｊ（ｊ＝１，・・・，ｍ、ｍは３以上の整数）における特徴点ｉ（ｉ＝１，・・・，ｎ、ｎは２以上の整数）の座標ｄｉｊ（ｘｉｊ，ｙｉｊ）を抽出するステップと、 該特徴点の２次元座標ｄｉｊ（ｘｉｊ，ｙｉｊ）より、該特徴点の３次元座標ｓｉ（Ｘｉ，Ｙｉ，Ｚｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算するステップと を包含し、 該特徴点の３次元座標ｓｉ（Ｘｉ，Ｙｉ，Ｚｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算するステップは、 行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求めるステップであって、行列Ｄは、 【数1】 と定義される、ステップと、 行列Ｄの特異値として行列Ｂの特異値σｊ（σ１≧σ２≧・・・≧σｒ＞０、ｒは、行列Ｄのランクに等しい）を求めるステップと、 σ１、σ２およびσ３に対する行列Ｄの特異ベクトルを求めるステップと、 Ｍ＝Ｍ’Ｃなる行列Ｃに対して、Ｅ＝ＣＣＴを満たす行列Ｅを計算するステップであって、Ｍ’＝Ｌ’（Σ’）１／２、Σ’はσ１、σ２およびσ３を対角要素に持ちそれ以外の要素が０である３×３行列、Ｌ’はσ１、σ２およびσ３に対応するＤの特異ベクトルを左から順に並べた行列である、ステップと、 行列Ｅから、行列Ｃを計算するステップと、 行列Ｃから、該３次元座標ｓｉ（Ｘｉ，Ｙｉ，Ｚｉ）および該変換を表す行列Ｍを計算するステップと を含み、 該σ１、σ２およびσ３に対する行列Ｄの特異ベクトルを求めるステップは、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列ＢＴＢ－σｊ２Ｉに対してツイスティッド分解を実行することにより、行列ＢＴＢを対角化するステップを含み、Ｉは単位行列である、方法。 【請求項8】 与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する方法であって、 該キーワードに対応する質問ベクトルｑを受け取るステップと、 該文書に対応する索引語文書行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求めるステップと、 行列Ｄの特異値として行列Ｂの特異値σｊ（σ１≧σ２≧・・・≧σｒ＞０、ｒは、行列Ｄのランクに等しい）を求めるステップと、 ｋ＜ｒなるｋを選択するステップと、 行列Ｄの擬似行列Ｄｋを計算するステップであって、行列Ｄｋは、σ１，σ２，・・・，σｋを対角要素としそれ以外の要素が０である行列Σｋ、σ１，σ２，・・・，σｋに対応する特異ベクトルのみを左から順に並べた左右直交行列Ｕｋ，Ｖｋを用いて、Ｄｋ＝ＵｋΣｋＶｋＴと定義される、ステップと、 行列Ｄｋと質問ベクトルｑとの類似度を計算するステップと、 該計算された類似度を基準に、検索結果を出力するステップと を包含し、 該行列Ｄｋの左右直交行列Ｕｋ，Ｖｋを求めるステップは、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列ＢＴＢ－σｊ２Ｉ（ｊ＝１，２，・・・，ｋ）に対してツイスティッド分解を実行することにより、行列ＢＴＢを対角化するステップを含み、Ｉは単位行列である、方法。 【請求項10】 物体の複数の２次元画像から３次元画像を復元する画像復元処理をコンピュータに実行させるプログラムであって、 該画像復元処理は、 ２次元画像ｊ（ｊ＝１，・・・，ｍ、ｍは３以上の整数）における特徴点ｉ（ｉ＝１，・・・，ｎ、ｎは２以上の整数）の座標ｄｉｊ（ｘｉｊ，ｙｉｊ）を抽出するステップと、 該特徴点の２次元座標ｄｉｊ（ｘｉｊ，ｙｉｊ）より、該特徴点の３次元座標ｓｉ（Ｘｉ，Ｙｉ，Ｚｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算するステップと を包含し、 該特徴点の３次元座標ｓｉ（Ｘｉ，Ｙｉ，Ｚｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算するステップは、 行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求めるステップであって、行列Ｄは、 【数2】 と定義される、ステップと、 行列Ｄの特異値として行列Ｂの特異値σｊ（σ１≧σ２≧・・・≧σｒ＞０、ｒは、行列Ｄのランクに等しい）を求めるステップと、 σ１、σ２およびσ３に対する行列Ｄの特異ベクトルを求めるステップと、 Ｍ＝Ｍ’Ｃなる行列Ｃに対して、Ｅ＝ＣＣＴを満たす行列Ｅを計算するステップであって、Ｍ’＝Ｌ’（Σ’）１／２、Σ’はσ１、σ２およびσ３を対角要素に持ちそれ以外の要素が０である３×３行列、Ｌ’はσ１、σ２およびσ３に対応するＤの特異ベクトルを左から順に並べた行列である、ステップと、 行列Ｅから、行列Ｃを計算するステップと、 行列Ｃから、該３次元座標ｓｉ（Ｘｉ，Ｙｉ，Ｚｉ）および該変換を表す行列Ｍを計算するステップと を含み、 該σ１、２およびσ３に対する行列Ｄの特異ベクトルを求めるステップは、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列ＢＴＢ－σｊ２Ｉに対してツイスティッド分解を実行することにより、行列ＢＴＢを対角化するステップを含み、Ｉは単位行列である、プログラム。 【請求項11】 与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する文書検索処理をコンピュータに実行させるプログラムであって、 該文書検索処理は、 該キーワードに対応する質問ベクトルｑを受け取るステップと、 該文書に対応する索引語文書行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求めるステップと、 行列Ｄの特異値として行列Ｂの特異値σｊ（σ１≧σ２≧・・・≧σｒ＞０、ｒは、行列Ｄのランクに等しい）を求めるステップと、 ｋ＜ｒなるｋを選択するステップと、 行列Ｄの擬似行列Ｄｋを計算するステップであって、行列Ｄｋは、σ１，σ２，・・・，σｋを対角要素としそれ以外の要素が０である行列Σｋ、σ１，σ２，・・・，σｋに対応する特異ベクトルのみを左から順に並べた左右直交行列Ｕｋ，Ｖｋを用いて、Ｄｋ＝ＵｋΣｋＶｋＴと定義される、ステップと、 行列Ｄｋと質問ベクトルｑとの類似度を計算するステップと、 該計算された類似度を基準に、検索結果を出力するステップと を包含し、 該行列Ｄｋの左右直交行列Ｕｋ，Ｖｋを求めるステップは、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列ＢＴＢ－σｊ２Ｉ（ｊ＝１，２，・・・，ｋ）に対してツイスティッド分解を実行することにより、行列ＢＴＢを対角化するステップを含み、Ｉは単位行列である、プログラム。 【請求項13】 物体の複数の２次元画像を３次元画像に復元する装置であって、 ｍ枚（ｍは３以上の整数）の２次元画像を受け取る手段と、 ２次元画像ｊ（ｊ＝１，・・・，ｍ）における特徴点ｉ（ｉ＝１，・・・，ｎ、ｎは２以上の整数）の座標ｄｉｊ（ｘｉｊ，ｙｉｊ）を抽出する手段と、 該特徴点の２次元座標ｄｉｊ（ｘｉｊ，ｙｉｊ）より、該特徴点の３次元座標ｓｉ（Ｘｉ，Ｙｉ，Ｚｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算する手段と を備え、 該特徴点の３次元座標ｓｉ（Ｘｉ，Ｙｉ，Ｚｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算する手段は、 行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求める手段であって、行列Ｄは、 【数3】 と定義される、手段と、 行列Ｄの特異値として行列Ｂの特異値σｊ（σ１≧σ２≧・・・≧σｒ＞０、ｒは、行列Ｄのランクに等しい）を求める手段と、 σ１、σ２およびσ３に対する行列Ｄの特異ベクトルを求める手段と、 Ｍ＝Ｍ’Ｃなる行列Ｃに対して、Ｅ＝ＣＣＴを満たす行列Ｅを計算する手段であって、Ｍ’＝Ｌ’（Σ’）１／２、Σ’はσ１、σ２およびσ３を対角要素に持ちそれ以外の要素が０である３×３行列、Ｌ’はσ１、σ２およびσ３に対応するＤの特異ベクトルを左から順に並べた行列である、手段と、 行列Ｅから、行列Ｃを計算する手段と、 行列Ｃから、該３次元座標ｓｉ（Ｘｉ，Ｙｉ，Ｚｉ）および該変換を表す行列Ｍを計算する手段と を含み、 該σ１、σ２およびσ３に対する行列Ｄの特異ベクトルを求める手段は、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列ＢＴＢ－σｊ２Ｉに対してツイスティッド分解を実行することにより、行列ＢＴＢを対角化する手段を含み、Ｉは単位行列である、装置。 【請求項14】 与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する装置であって、 該キーワードに対応する質問ベクトルｑを受け取る手段と、 該文書に対応する索引語文書行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求める手段と、 行列Ｄの特異値として行列Ｂの特異値σｊ（σ１≧σ２≧・・・≧σｒ＞０、ｒは、行列Ｄのランクに等しい）を求める手段と、 ｋ＜ｒなるｋを選択する手段と、 行列Ｄの擬似行列Ｄｋを計算する手段であって、行列Ｄｋは、σ１，σ２，・・・，σｋを対角要素としそれ以外の要素が０である行列Σｋ、σ１，σ２，・・・，σｋに対応する特異ベクトルのみを左から順に並べた左右直交行列Ｕｋ，Ｖｋを用いて、Ｄｋ＝ＵｋΣｋＶｋＴと定義される、手段と、 行列Ｄｋと質問ベクトルｑとの類似度を計算する手段と、 該計算された類似度を基準に、検索結果を出力する手段と を包含し、 該行列Ｄｋの左右直交行列Ｕｋ，Ｖｋを求める手段は、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列ＢＴＢ－σｊ２Ｉ（ｊ＝１，２，・・・，ｋ）に対してツイスティッド分解を実行することにより、行列ＢＴＢを対角化する手段を含み、Ｉは単位行列である、装置。
発明の詳細な説明	【技術分野】 【0001】 本発明は、任意の行列Ａに対して高速かつ高精度に特異値分解を実行する方法、プログラムおよび装置に関し、より詳細には、行列の特異値を計算し、計算された特異値に対する特異ベクトルを計算することによって、高速かつ高精度に特異値分解を実行する方法、プログラムおよび装置に関する。 【背景技術】 【0002】 現在、コンピュータ上で特異値分解を行うために用いられる標準的な方法は、国際的な数値計算ライブラリＬＡＰＡＣＫにより公開されているＤＢＤＳＱＲルーチンである。ＤＢＤＳＱＲルーチンは、ＱＲ法に基づいて作られたものであり、前処理、メインループ、後処理の３つの部分に分けることができる。前処理では、特異値の上限および下限を計算し、収束判定に用いる精度を計算する。メインループでは、ＱＲ法を繰り返しながら、除々に行列の分割および減次を行い、最終的に行列サイズが１となった時点でループを終了する。途中で２×２行列のブロックが現れたならば別処理を行う。後処理では、特異値として計算された値が負ならば正に直し、必要ならば特異ベクトルをスカラー倍する。特異値を降順に、それに対応するように特異ベクトルも並べ替える。ＤＢＤＳＱＲは、計算量が非常に多く、特に大規模問題では計算時間の増大が避けられない。ＤＢＤＳＱＲルーチンでは、特異値および特異ベクトルが同時に計算される。また、ＬＡＰＡＣＫには、特異値を計算するＤＬＡＳＱルーチン、および、計算された特異値を用いて特異ベクトルを計算する際に用いる対角化のためのＤＳＴＥＧＲルーチンがある。ＤＬＡＳＱルーチンは、高速高精度に特異値を求めることができるが、特異ベクトルを計算できない。ＤＳＴＥＧＲルーチンは、その数値的な欠点から、実際に特異値分解に使用することは難しい。 【0003】 ＬＡＰＡＣＫによるＤＢＤＳＱＲルーチンを例として、従来の特異値分解方法を説明する。ＤＢＤＳＱＲルーチンでは、ｌ１×ｌ２の一般行列Ａを特異値分解するために、まず、行列Ａに対してＨｏｕｓｅｈｏｌｄｅｒ（ハウスホルダー）変換を適用する。すなわち、行列Ａは、直交行列ＵＡ、ＶＡを用いて、 【0004】 【数4】 と表わすことができる。このとき得られるＢを、上２重対角行列という。ここで、Ａの特異値＝Ｂの特異値であることに注意する。このように、一般行列Ａに対する特異値分解問題を上２重対角行列Ｂに対する特異値分解問題 Ｂ＝ＵＢΣＶＢＴ 行列ＵＢ、ＶＢは、それぞれ左および右直交行列 Σ≡ｄｉａｇ（σ１，σ２，・・・，σｍ） σ１≧σ２≧・・・≧σｍ≧０ σｊはＢの特異値 に置き換える。 【0005】 次に、行列ＢＴＢを考える。この行列ＢＴＢの対角化 Λ＝ＶＴＢＴＢＶ Λ≡ｄｉａｇ（λ１，λ２，・・・，λｍ） λ１≧λ２≧・・・≧λｍ≧０ Ｖ≡（ｖ１，ｖ２，・・・，ｖｍ） λｊはＢＴＢの固有値 ｖｊは固有値λｊに対する固有ベクトル を行う。ここで、一般に、以下のことが成り立つ。（１）ＢＴＢは対称な３重対角行列である。（２）ＢＴＢの固有値は全て正であり、Ｂの特異値σｊ（σ１≧σ２≧・・・≧σｍ≧０）は、ＢＴＢの固有値λｊ（λ１≧λ２≧・・・≧λｍ≧０）の正の平方根に等しい。（３）ＶＢ＝Ｖ、すなわち、ＢＴＢの固有ベクトルｖｊは、Ｂの右特異ベクトルに等しい。従って、行列ＢＴＢの対角化が求まると、（２）よりΛ＝Σ２であることから、Σが求まり、さらに（３）より左直交行列ＵＢ＝ＢＶＢΣ－１＝ＢＶΣ－１が求まるので、Ｂが特異値分解される。すなわち、Ｂの特異値分解を、ＢＴＢの対角化の問題に置き換えることができる。この原理は、ｍ個全ての特異値および特異ベクトルを求める場合だけでなく、少なくとも１つの特異値および特異ベクトルを求める場合にも適用され得る。 【0006】 以上のように、一般行列Ａの特異値分解は、対称な３重対角行列ＢＴＢの対角化の問題を含む。 【発明の開示】 【発明が解決しようとする課題】 【0007】 ＤＢＤＳＱＲルーチンでは、計算量が非常に多いため、特に大規模問題では、高速に特異値分解を実行することが困難であった。一方、ＤＬＡＳＱルーチンは、高速かつ高精度に特異値を求めることができるが、ＤＳＴＥＧＲルーチンは、場合によっては特異ベクトルを低精度で計算するため、ＤＳＴＥＧＲルーチンでは、常に高精度に特異値分解を実行することは困難であった。 【0008】 本発明の目的の１つは、対称な３重対角行列ＢＴＢの対角化を高速かつ高精度で実行することにより、任意の行列Ａの特異値分解を高速かつ高精度に実行することを可能にする方法、プログラム、および装置を提供することにある。ここで、行列Ｂを、行列Ａをハウスホルダー変換することで得られる上２重対角行列とする。 【課題を解決するための手段】 【0009】 本発明では、ミウラ型逆変換、ｓｄＬＶｖｓ変換、ｒｄＬＶｖｓ変換およびミウラ型変換によるＴｗｉｓｔｅｄ（ツイスティッド）分解を実行することにより、行列ＢＴＢを対角化する。ここで、行列Ｂを、行列Ａと特異値が同じ上２重対角行列とする。行列Ａから行列Ｂへの変換は、例えばハウスホルダー法で実行できる。本発明による行列ＢＴＢの対角化では、固有値および固有ベクトルを同時に求めるＤＢＤＳＱＲルーチンとは異なり、ＤＳＴＥＧＲルーチンのように、まず固有値を計算し、次に計算された固有値を用いて固有ベクトルを求める。 【0010】 本発明により、コンピュータを用いて任意の行列Ａに対して特異値分解を実行する方法であって、行列Ａに対して上２重対角化を行い、行列Ａの上２重対角行列Ｂを求めるステップと、行列Ａの特異値として行列Ｂの少なくとも１つの特異値σを求めるステップと、σに対する行列Ａの特異ベクトルを求めるステップとを包含し、該行列Ａの特異ベクトルを求めるステップは、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列ＢＴＢ－σ２Ｉに対してツイスティッド分解を実行することにより、行列ＢＴＢを対角化するステップを含み、Ｉは単位行列である、方法が提供され、これにより上記目的が達成される。 【0011】 前記行列Ａの特異ベクトルを求めるステップは、前記行列ＢＴＢを対角化するステップの後に、逆反復法を実行するステップを含んでもよい。 【0012】 前記行列Ａの特異ベクトルを求めるステップは、前記逆反復法を実行するステップの後に、グラムシュミット法を実行するステップを含んでもよい。 【0013】 前記行列Ｂの特異値σを求めるステップは、ｄＬＶｓルーチンを実行するステップを含んでもよい。 【0014】 前記行列Ｂの特異値σを求めるステップは、要求される計算時間および精度に応じて、ｄＬＶｓルーチンもしくはＤＬＡＳＱルーチンのどちらを実行するのかを判定するステップを含んでもよい。 【0015】 前記行列Ａに対して上２重対角化を行い、上２重対角行列Ｂを求めるステップは、上２重対角化法（たとえばハウスホルダー法）を実行するステップを含んでもよい。 【0016】 本発明により、物体の複数の２次元画像から３次元画像を復元する方法であって、２次元画像ｊ（ｊ＝１，・・・，ｍ、ｍは３以上の整数）における特徴点ｉ（ｉ＝１，・・・，ｎ、ｎは２以上の整数）の座標ｄｉｊ（ｘｉｊ，ｙｉｊ）を抽出するステップと、該特徴点の２次元座標ｄｉｊ（ｘｉｊ，ｙｉｊ）より、該特徴点の３次元座標ｓｉ（Ｘｉ，Ｙｉ，Ｚｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算するステップとを包含し、該特徴点の３次元座標ｓｉ（Ｘｉ，Ｙｉ，Ｚｉ）および２次元座標から３次元 座標への変換を表す行列Ｍを計算するステップは、行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求めるステップであって、行列Ｄは、 【0017】 【数5】 と定義される、ステップと、行列Ｄの特異値として行列Ｂの特異値σｊ（σ１≧σ２≧・・・≧σｒ＞０、ｒは、行列Ｄのランクに等しい）を求めるステップと、σ１、σ２およびσ３に対する行列Ｄの特異ベクトルを求めるステップと、Ｍ＝Ｍ’Ｃなる行列Ｃに対して、Ｅ＝ＣＣＴを満たす行列Ｅを計算するステップであって、Ｍ’＝Ｌ’（Σ’）１／２、Σ’はσ１、σ２およびσ３を対角要素に持ちそれ以外の要素が０である３×３行列、Ｌ’はσ１、σ２およびσ３に対応するＤの特異ベクトルを左から順に並べた行列である、ステップと、行列Ｅから、行列Ｃを計算するステップと、行列Ｃから、該３次元座標ｓｉ（Ｘｉ，Ｙｉ，Ｚｉ）および該変換を表す行列Ｍを計算するステップとを含み、該σ１、σ２およびσ３に対する行列Ｄの特異ベクトルを求めるステップは、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列ＢＴＢ－σｊ２Ｉに対してツイスティッド分解を実行することにより、行列ＢＴＢを対角化するステップを含み、Ｉは単位行列である、方法が提供され、これにより上記目的が達成される。 【0018】 本発明により、与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する方法であって、該キーワードに対応する質問ベクトルｑを受け取るステップと、該文書に対応する索引語文書行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求めるステップと、行列Ｄの特異値として行列Ｂの特異値σｊ（σ１≧σ２≧・・・≧σｒ＞０、ｒは、行列Ｄのランクに等しい）を求めるステップと、ｋ＜ｒなるｋを選択するステップと、行列Ｄの擬似行列Ｄｋを計算するステップであって、行列Ｄｋは、σ１，σ２，・・・，σｋを対角要素としそれ以外の要素が０である行列Σｋ、σ１，σ２，・・・，σｋに対応する特異ベクトルのみを左から順に並べた左右直交行列Ｕｋ，Ｖｋを用いて、Ｄｋ＝ＵｋΣｋＶｋＴと定義される、ステップと、行列Ｄｋと質問ベクトルｑとの類似度を計算するステップと、該計算された類似度を基準に、検索結果を出力するステップとを包含し、該行列Ｄｋの左右直交行列Ｕｋ，Ｖｋを求めるステップは、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列ＢＴＢ－σｊ２Ｉ（ｊ＝１，２，・・・，ｋ）に対してツイスティッド分解を実行することにより、行列ＢＴＢを対角化するステップを含み、Ｉは単位行列である、方法が提供され、これにより上記目的が達成される。 【0019】 本発明により、任意の行列Ａに対して特異値分解を実行する特異値分解処理をコンピュータに実行させるプログラムであって、該特異値分解処理は、行列Ａに対して上２重対角化を行い、行列Ａの上２重対角行列Ｂを求めるステップと、行列Ａの特異値として行列Ｂの少なくとも１つの特異値σを求めるステップと、σに対する行列Ａの特異ベクトルを求めるステップとを包含し、該行列Ａの特異ベクトルを求めるステップは、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列ＢＴＢ－σ２Ｉに対してツイスティッド分解を実行することにより、行列ＢＴＢを対角化するステッ プを含み、Ｉは単位行列である、プログラムが提供され、これにより上記目的が達成される。 【0020】 本発明により、物体の複数の２次元画像から３次元画像を復元する画像復元処理をコンピュータに実行させるプログラムであって、該画像復元処理は、２次元画像ｊ（ｊ＝１，・・・，ｍ、ｍは３以上の整数）における特徴点ｉ（ｉ＝１，・・・，ｎ、ｎは２以上の整数）の座標ｄｉｊ（ｘｉｊ，ｙｉｊ）を抽出するステップと、該特徴点の２次元座標ｄｉｊ（ｘｉｊ，ｙｉｊ）より、該特徴点の３次元座標ｓｉ（Ｘｉ，Ｙｉ，Ｚｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算するステップとを包含し、該特徴点 の３次元座標ｓｉ（Ｘｉ，Ｙｉ，Ｚｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算するステップは、行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求めるステップであって、行列Ｄは、 【0021】 【数6】 と定義される、ステップと、行列Ｄの特異値として行列Ｂの特異値σｊ（σ１≧σ２≧・・・≧σｒ＞０、ｒは、行列Ｄのランクに等しい）を求めるステップと、σ１、σ２およびσ３に対する行列Ｄの特異ベクトルを求めるステップと、Ｍ＝Ｍ’Ｃなる行列Ｃに対して、Ｅ＝ＣＣＴを満たす行列Ｅを計算するステップであって、Ｍ’＝Ｌ’（Σ’）１／２、Σ’はσ１、σ２およびσ３を対角要素に持ちそれ以外の要素が０である３×３行列、Ｌ’はσ１、σ２およびσ３に対応するＤの特異ベクトルを左から順に並べた行列である、ステップと、行列Ｅから、行列Ｃを計算するステップと、行列Ｃから、該３次元座標ｓｉ（Ｘｉ，Ｙｉ，Ｚｉ）および該変換を表す行列Ｍを計算するステップとを含み、該σ１、σ２およびσ３に対する行列Ｄの特異ベクトルを求めるステップは、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列ＢＴＢ－σｊ２Ｉに対してツイスティッド分解を実行することにより、行列ＢＴＢを対角化するステップを含み、Ｉは単位行列である、プログラムが提供され、これにより上記目的が達成される。 【0022】 本発明により、与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する文書検索処理をコンピュータに実行させるプログラムであって、該文書検索処理は、該キーワードに対応する質問ベクトルｑを受け取るステップと、該文書に対応する索引語文書行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求めるステップと、行列Ｄの特異値として行列Ｂの特異値σｊ（σ１≧σ２≧・・・≧σｒ＞０、ｒは、行列Ｄのランクに等しい）を求めるステップと、ｋ＜ｒなるｋを選択するステップと、行列Ｄの擬似行列Ｄｋを計算するステップであって、行列Ｄｋは、σ１，σ２ ，・・・，σｋを対角要素としそれ以外の要素が０である行列Σｋ、σ１，σ２，・・・，σｋに対応する特異ベクトルのみを左から順に並べた左右直交行列Ｕｋ，Ｖｋを用いて、Ｄｋ＝ＵｋΣｋＶｋＴと定義される、ステップと、行列Ｄｋと質問ベクトルｑとの類似度を計算するステップと、該計算された類似度を基準に、検索結果を出力するステップとを包含し、該行列Ｄｋの左右直交行列Ｕｋ，Ｖｋを求めるステップは、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列ＢＴＢ－σｊ ２Ｉ（ｊ＝１，２，・・・，ｋ）に対してツイスティッド分解を実行することにより、行列ＢＴＢを対角化するステップを含み、Ｉは単位行列である、プログラムが提供され、これにより上記目的が達成される。 【0023】 本発明により、任意の行列Ａに対して特異値分解を実行する装置であって、行列Ａを入力として受け取る手段と、行列Ａに対して上２重対角化を行い、行列Ａの上２重対角行列Ｂを計算する手段と、行列Ａの特異値として行列Ｂの少なくとも１つの特異値σを計算する手段と、σに対する行列Ａの特異ベクトルを求める手段とを備え、該行列Ａの特異ベクトルを求める手段は、ミウラ型逆変換、ｓｄＬＶｖｓ変換およびｒｄＬＶｖｓ変換、ならびにミウラ型変換を用いて行列ＢＴＢ－σ２Ｉに対してツイスティッド分解を実行することにより、行列ＢＴＢを対角化する手段を含み、Ｉは単位行列である、装置が提供され、これにより上記目的が達成される。 【0024】 本発明により、物体の複数の２次元画像を３次元画像に復元する装置であって、ｍ枚（ｍは３以上の整数）の２次元画像を受け取る手段と、２次元画像ｊ（ｊ＝１，・・・，ｍ）における特徴点ｉ（ｉ＝１，・・・，ｎ、ｎは２以上の整数）の座標ｄｉｊ（ｘｉｊ，ｙｉｊ）を抽出する手段と、該特徴点の２次元座標ｄｉｊ（ｘｉｊ，ｙｉｊ）より、該特徴点の３次元座標ｓｉ（Ｘｉ，Ｙｉ，Ｚｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算する手段とを備え、該特徴点の３次元座標ｓｉ（Ｘｉ，Ｙｉ，Ｚｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算する手段は、行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求める手段であって、行列Ｄは、 【0025】 【数7】 と定義される、手段と、行列Ｄの特異値として行列Ｂの特異値σｊ（σ１≧σ２≧・・・≧σｒ＞０、ｒは、行列Ｄのランクに等しい）を求める手段と、σ１、σ２およびσ３に対する行列Ｄの特異ベクトルを求める手段と、Ｍ＝Ｍ’Ｃなる行列Ｃに対して、Ｅ＝ＣＣＴを満たす行列Ｅを計算する手段であって、Ｍ’＝Ｌ’（Σ’）１／２、Σ’はσ１、σ２およびσ３を対角要素に持ちそれ以外の要素が０である３×３行列、Ｌ’はσ１、σ２およびσ３に対応するＤの特異ベクトルを左から順に並べた行列である、手段と、行列Ｅから、行列Ｃを計算する手段と、行列Ｃから、該３次元座標ｓｉ（Ｘｉ，Ｙｉ，Ｚｉ）および該変換を表す行列Ｍを計算する手段とを含み、該σ１、σ２およびσ３に対する行列Ｄの特異ベクトルを求める手段は、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列ＢＴＢ－σｊ２Ｉに対してツイスティッド分解を実行することにより、行列ＢＴＢを対角化する手段を含み、Ｉは単位行列である、装置が提供され、これにより上記目的が達成される。 【0026】 本発明により、与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する装置であって、該キーワードに対応する質問ベクトルｑを受け取る手段と、該文書に対応する索引語文書行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求める手段と、行列Ｄの特異値として行列Ｂの特異値σｊ（σ１≧σ２≧・・・≧σｒ＞０、ｒは、行列Ｄのランクに等しい）を求める手段と、ｋ＜ｒなるｋを選択する手段と、行列Ｄの擬似行列Ｄｋを計算する手段であって、行列Ｄｋは、σ１，σ２，・・・，σｋを対角要素としそれ以外の要素が０である行列Σｋ、σ１，σ２，・・・，σｋに対応する特異ベクトルのみを左から順に並べた左右直交行列Ｕｋ，Ｖｋを用いて、Ｄｋ＝ＵｋΣｋＶｋＴと定義される、手段と、行列Ｄｋと質問ベクトルｑとの類似度を計算する手段と、該計算された類似度を基準に、検索結果を出力する手段とを包含し、該行列Ｄｋの左右直交行列Ｕｋ，Ｖｋを求める手段は、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列ＢＴＢ－σｊ２Ｉ（ｊ＝１，２，・・・，ｋ）に対してツイスティッド分解を実行することにより、行列ＢＴＢを対角化する手段を含み、Ｉは単位行列である、装置が提供され、これにより上記目的が達成される。 【発明の効果】 【0027】 本発明によれば、ミウラ型逆変換、ｓｄＬＶｖｓ変換、ｒｄＬＶｖｓ変換およびミウラ型変換によるツイスティッド分解を実行することにより、行列ＢＴＢが対角化される。それにより、ＤＳＴＥＧＲルーチンと比較した場合、精度の点ではるかに優れ、実用化に耐え得る固有値計算を実現し、高精度の特異値分解が可能になる。また、逆反復法、Ｇｒａｍ－Ｓｃｈｍｉｄｔ（グラムシュミット）法、本発明において提供するｄＬＶｓルーチンのうちの少なくとも１つを併用することにより、ＤＢＤＳＱＲと比較して精度の点では同等になることもあるが、速度の点では圧倒的に高速である特異値分解を実現することができる。さらに、ＤＢＤＳＱＲルーチンと比較すると、固有値および固有ベクトルを同時の求めるのか、または、計算された固有値を用いて固有ベクトルを求めるのかの計算フローの違いにより、以下で説明する２次元から３次元への画像復元問題および文書検索問題等においては、計算時間をさらに短縮することができる。 【図面の簡単な説明】 【0028】 【図1】図１は、本発明により提供される特異値分解のためのＩ－ＳＶＤルーチンの処理の手順を示すフローチャート 【図2】図２は、Ｐａｒｌｅｔｔらの手法によるｓｑｄｓ変換およびｒｐｑｄｓ変換の処理の手順を示すフローチャート 【図3】図３は、本発明によるミウラ型逆変換、ｓｄＬＶｖｓ変換、ｒｄＬＶｖｓ変換およびミウラ型変換の処理の手順を示すフローチャート 【図4】図４は、本発明によるミウラ型逆変換、ｓｄＬＶｖｓ変換、ｒｄＬＶｖｓ変換およびミウラ型変換の処理の手順を示すフローチャート 【図5】図５は、本発明によるミウラ型逆変換、ｓｄＬＶｖｓ変換、ｒｄＬＶｖｓ変換およびミウラ型変換の処理の手順を示すフローチャート。 【図6】図６は、本発明によるミウラ型逆変換、ｓｄＬＶｖｓ変換、ｒｄＬＶｖｓ変換およびミウラ型変換の処理の手順を示すフローチャート 【図7】図７は、本発明によるミウラ型逆変換、ｓｄＬＶｖｓ変換、ｒｄＬＶｖｓ変換およびミウラ型変換の処理の手順を示すフローチャート 【図8】図８は、本発明によるミウラ型逆変換、ｓｄＬＶｖｓ変換、ｒｄＬＶｖｓ変換およびミウラ型変換の処理の手順を示すフローチャート 【図9】図９は、本発明において誤差の低減をもたらす、Ｐａｒｌｅｔｔの手法と本発明による方法との計算経路の違いを示す図 【図10】図１０は、本発明による特異値分解装置の１つの実施形態であるコンピュータを示す図 【図11】図１１は、本発明による特異値分解方法を利用した物体の複数の２次元画像を３次元画像へ復元する画像処理の１つの実施形態を示すフローチャート 【図12】図１２は、本発明による特異値分解装置を利用した文書検索方法の１つの実施形態を示すフローチャート 【符号の説明】 【0029】 １０ コンピュータ １００１ ＣＰＵ １００２ メモリ １００３ 入力インターフェース部 １００４ 出力インターフェース部 １００５ バス 【発明を実施するための最良の形態】 【0030】 以下、図面を参照しながら、本発明の実施形態を説明する。 【0031】 １．特異値分解アルゴリズムＩ－ＳＶＤルーチン 図１を参照する。図１は、本発明により提供される特異値分解のためのＩ－ＳＶＤルーチンの処理の手順を示す。図１の処理は、コンピュータプログラムの形式で提供され得る。 【0032】 ステップ１０２において、一般行列Ａに対してハウスホルダー法を用いて上２重対角化を行う。このとき利用される上２重対角化方法は、ハウスホルダー法であってもよいし、他の上２重対角化方法であってもよい。これにより、上２重対角行列Ｂが得られる。次に、行列Ｂの特異値σｊ（σ１≧σ２≧・・・≧σｍ≧０、ｍは、行列ＢＴＢの行列サイズに等しい）を計算する。ここで、本発明により提供される特異値計算ルーチンｄＬＶｓおよびＬＡＰＡＣＫにより提供されるＤＬＡＳＱルーチンを利用し得る。なお、行列Ｂの特異値は、行列Ａの特異値に等しく、かつ、行列ＢＴＢの固有値の正の平方根λｊ１／２に等しい。ｄＬＶｓルーチンは、Ｂの特異値σｊを高精度で計算する。ｄＬＶｓルーチンの詳細は、後述される。ｄＬＶｓで求められる特異値は、従来のＬＡＰＣＫルーチン（ＤＢＤＳＱＲ、ＤＬＡＳＱ等）と比較して、最も高精度である。速度に関しては、ほとんどの場合ＤＬＡＳＱルーチンには劣るが、信頼性の高いＤＢＤＳＱＲルーチンの約２分の１の時間で計算できる(ＣＰＵがＩｔａｎｉｕｍ２の場合はｄＬＶｓルーチンが最速)。そこで、ステップ１０４において、高精度性が必要とされるかどうかを判定する。最も高精度な特異値分解を行いたい場合には、ステップ１０６においてｄＬＶｓを、高精度性および高速性を両立したい場合には、ステップ１０８において従来のＤＬＡＳＱルーチンを、特異値計算ルーチンとして用いるのがよい。さらに、ＤＬＡＳＱルーチンが有限回の演算で停止するかは定かでないが、ｄＬＶｓルーチンは収束することが保証されているので、信頼性を重視するならば、ｄＬＶｓルーチンを用いるのがよい。 【0033】 次に、ステップ１１０において、計算されたＢの特異値σｊ（すなわち、Ａの特異値に等しい）用いて、１．１において説明されるｓｄＬＶｖｓ変換およびｒｄＬＶｖｓ変換によるツイスティッド分解により、ＢＴＢの対角化を行う。本発明による行列ＢＴＢの対角化ルーチンの詳細は、１．１で説明する。さらに、ステップ１１２において、以下の１．１にて詳細が示される簡略化連立１次方程式 【0034】 【数8】 を求めると、行列Ｂの特異ベクトルｖｊが得られる（すなわち、ＢＴＢの対角化行列Ｖが得られる）。すなわち、ＢＴＢが対角化される。ステップ１１４において、Σ≡ｄｉａｇ（σ１，σ２，・・・，σｍ）（ただし、σ１≧σ２≧・・・≧σｍ≧０）とすると、行列Ｂの右直交行列ＶＢ＝Ｖであるので、行列Ｂの左直交行列ＵＢ＝ＢＶＢΣ－１＝ＢＶΣ－１からＵＢを得る。すなわち、行列Ｂが特異値分解される。ここでは、ｍ個全ての特異値および特異ベクトルを計算したが、この原理は、少なくとも１つの特異値および特異ベクトルを計算する場合にも適用され得る。 【0035】 このとき得られる左右直交行列、すなわち特異ベクトル計算の結果は、ＤＢＤＳＱＲルーチンによって得られる結果と比較して精度の点で劣る。そこで、ステップ１１６において、逆反復法を用いることにより、ＤＢＤＳＱＲルーチンと同程度の精度を得ることができる。さらなる高精度性が必要とされる場合は、ステップ１２０において、グラムシュミット法による再直交化を行えば、ＤＢＤＳＱＲよりも高精度での計算を実現し得る。Ｉ－ＳＶＤルーチンでは、速度に関しては、表２から、ＤＢＤＳＱＲと比較して大幅に短縮され得ることが理解される。また、表１から理解されるように、Ｉ－ＳＶＤルーチンは、ＤＳＴＥＧＲによる極度の精度悪化がみられる場合でも、高精度に特異値分解することができる。なお、表１、２の誤差および計算時間に関するデータは、逆反復法を用いているが、グラムシュミット法は用いていない場合の結果である。また、ステップ１１６および１２０を実行する順序は逆であってもよいし、あるいは、両ステップは、省略してもよい。 【0036】 このように、図１に示される各ステップの機能は、ソフトウェア（例えば、プログラム）によって実現される。しかし、本発明は、これに限定されない。図１に示される各ステップの機能をハードウェア（例えば、回路、ボード、半導体チップ）によって実現してもよいし、ソフトウェアとハードウェアとの組み合わせによって実現してもよい。 【0037】 以上により、Ｉ－ＳＶＤルーチンを用いると、従来技術と比較して高速かつ高精度の特異値分解が実行される。 【0038】 表１は、本発明による行列ＢＴＢ対角化ルーチンを実行した後に、逆反復法を実行したときの従来技術および本発明の誤差を比較した表である。 【0039】 【表1】 表２は、本発明による行列ＢＴＢ対角化ルーチンを実行した後に、逆反復法を実行したときの従来技術および本発明の計算時間を比較した表である。 【0040】 【表2】 また、１.３に、本発明の方法において、ＤＳＴＥＧＲルーチンと比較して、誤差が低減される理由を説明する。 【0041】 １．１．本発明による行列ＢＴＢの対角化ルーチン ＤＢＤＳＱＲはＢＴＢの固有値λｊとともに対応する固有ベクトルｖｊを求めることができるので、他のルーチンを併用しなくてもＢＴＢの対角化（Ｂの特異値分解）が得られる。しかしながら、ＤＬＡＳＱおよびｄＬＶｓは上述したように固有ベクトルｖｊを計算する機能はない。よって、別の固有ベクトル計算ルーチンが必要となる。ただし、ＢＴＢの固有値λｊ（Ｂの特異値σｊ）が求まっているものとする。 【0042】 まず、ＤＳＴＥＧＲによる固有ベクトルの計算方法を示す。Ｐａｒｌｅｔｔらは以下の（１）、（２）、（３）の手順で固有ベクトル計算ができることを示している。 （１）ｓｔａｔｉｏｎａｒｙ ｑｄ ｗｉｔｈ ｓｈｉｆｔ変換（ｓｔｑｄｓ変換）およびｒｅｖｅｒｓｅ－ｔｉｍｅ ｐｒｏｇｒｅｓｓｉｖｅ ｑｄ ｗｉｔｈ ｓｈｉｆｔ変換（ｒｐｑｄｓ変換）を用いた（Ｂ（０））ＴＢ（０）－λｊＩのＣｈｏｌｅｓｋｙ（コレスキー）分解（Ｂ（±１））ＴＢ（±１）＝（Ｂ（０））ＴＢ（０）－λｊＩを求める。ただし、Ｂ（０）＝Ｂであり、Ｂ（＋１）およびＢ（－１）はそれぞれ上２重対角および下２重対角行列とする。すなわち 【0043】 【数9】 となる。すなわち、｛ｑｋ（±１），ｅｋ（±１）｝を計算する。 （２）コレスキー分解より（Ｂ（０））ＴＢ（０）－λｋＩのツイスティッド分解 【0044】 【数10】 を形成する。ただし 【0045】 【数11】 とし、γｋ≡ｑｋ（１）＋ｑｋ（－１）－（ｅｋ－１（０）＋ｑｋ（０）－λｋ）が最小になるｋをｋ０とする。 （３）簡約化された連立１次方程式 【0046】 【数12】 を解く。 【0047】 図２は、Ｐａｒｌｅｔｔらの手法で最も核となる（１）の処理手順を示す。（１）で求めたいのは、｛ｑｋ（ｎ），ｅｋ（ｎ）｝，ｎ＝±１であるが、そのための変換としてｓｔａｔｉｏｎａｒｙ ｑｄ ｗｉｔｈ ｓｈｉｆｔ（ｓｔｑｄｓ）変換 ｑｋ（１）＋ｅｋ－１（１）＝ｑｋ（０）＋ｅｋ－１（０）－λｊ，ｋ＝１，２，・・・，ｍ， ｑｋ（１）ｅｋ（１）＝ｑｋ（０）ｅｋ（０），ｋ＝１，２，・・・，ｍ－１， ｅ０（０）≡０，ｅ０（１）≡０ およびｒｅｖｅｒｓｅ－ｔｉｍｅ ｐｒｏｇｒｅｓｓｉｖｅ ｑｄ ｗｉｔｈ ｓｈｉｆｔ（ｒｐｑｄｓ）変換 ｑｋ（－１）＋ｅｋ（－１）＝ｑｋ（０）＋ｅｋ－１（０）－λｊ，ｋ＝１，２，・・・，ｍ， ｑｋ＋１（－１）ｅｋ（－１）＝ｑｋ（０）ｅｋ（０），ｋ＝１，２，・・・，ｍ－１， ｅ０（０）≡０，ｅｍ（－１）≡０ が採用されている。 【0048】 固有値λｊが既知ならば反復的な計算は不要なため計算量がＱＲアルゴリズムに比べて圧倒的に少なく抑えられるが、常に数値安定かつ精度のよい方法とはいえない。ｓｔｑｄｓ変換とｒｐｑｄｓ変換ともに減算による桁落ちが発生する可能性がある。例えば、ｓｔｑｄｓ変換においてｑｋ（０）＋ｅｋ－１（０）－λｊ～ｅｋ－１（１）ならば、ｑｋ（１）（＝ｑｋ（０）＋ｅｋ－１（０）－λｊ－ｅｋ－１（１））を求める際、倍精度計算であってもｑｋ（１）の有効桁はわずか１桁になることもある。その場合、ｑｋ（０）ｅｋ（０）／ｑｋ（１）を計算すると誤差が生じる、つまり、ｅｋ（１）が精度よく計算できないことになる。また、ｑｋ＋１（１）を求めるのにｅｋ（１）が必要、ｅｋ（１）を求めるのにｑｋ（１）が必要と逐次的な計算が要求されるので、１個所で発生した桁落ちによる誤差が波及し、さらなる誤差増大の可能性も秘めている。その結果、理論上はｑｋ（１）≠０であるが誤差蓄積によりｑｋ（１）＝０となり、ｑｋ（０）ｅｋ（０）／ｑｋ（１）の計算においてオーバーフローが起こるといった数値不安定な状況も想定される。 Ｂ（０）＝Ｂの成分｛ｂ２ｋ－１（０），ｂ２ｋ（０）｝が与えられる、すなわち｛ｑｋ（０），ｅｋ（０）｝が与えられると、λｊおよびｅｋ－１（１）が一意的に決まるので、この状況は避けることはできない。ｒｐｑｄｓ変換も同様の性質を持つため、実用的なレベルにまで達したとはいいがたい。ＬＡＰＡＣＫにおいてＦＯＲＴＲＡＮルーチンＤＳＴＥＧＲとして改良版が公開されているものの欠点は完全に解決されてはいない。 【0049】 次に、本発明による対角化のルーチンを説明する。本発明によるルーチンでは、ｓｔｑｄｓ変換およびｒｐｑｄｓ変換にかわる新たな変換を採用している。本発明では、行列ＢＴＢを対角化するために、ｓｄＬＶｖｓ変換およびｒｄＬＶｖｓ変換を用いた行列ＢＴＢ－σｊ２Ｉに対するツイスティッド分解を行う。正確には、上記（１）を以下の（１ａ）、（１ｂ）、（１ｃ）の３つの工程にわけることで数値安定なコレスキー分解を実現している。 【0050】 図３～図８は、本発明によるミウラ型逆変換、ｓｄＬＶｖｓ変換、ｒｄＬＶｖｓ変換およびミウラ型変換を実行する処理手順を示す。 （１ａ）ミウラ型逆変換 【0051】 【数13】 を行う。 【0052】 ｕ２ｋ－１（０）＝η（０）ｔｋ（０），ｋ＝１，２，・・・，ｍ， ｕ２ｋ（０）＝ｅｋ（０）／ｔｋ（０），ｋ＝１，２，・・・，ｍ－１， ｔｋ（０）≡ｑｋ（０）／（η（０）＋ｕ２ｋ－２（０））－１，η（０）≡１／δ（０） ただし、δ（０）は任意に選ぶことができる。 （１ｂ）ｓｔａｔｉｏｎａｒｙ不等間隔差分Ｌｏｔｋａ－Ｖｏｌｔｅｒｒａ（ｓｄＬＶｖｓ）変換 ｕｋ（１）＝ｕｋ（０）＊（１＋δ（０）ｕｋ－１（０））／（１＋δ（１）ｕｋ－１（１））， ｋ＝１，２，・・・，２ｍ－１， によって 【0053】 【数14】 ｒｅｖｅｒｓｅ－ｔｉｍｅ不等間隔差分Ｌｏｔｋａ－Ｖｏｌｔｅｒｒａ（ｒｄＬＶｖｓ）変換 ｕｋ（－１）＝ｕｋ（０）＊（１＋δ（０）ｕｋ－１（０））／（１＋δ（－１）ｕ ｋ＋１（－１））， ｋ＝１，２，・・・，２ｍ－１， によって 【0054】 【数15】 を実行する。ただし、ｕｋ（ｎ）≡０，ｕ２ｍ（０）≡０とし、λｊ＝１／δ（０）－１／δ（ｎ），ｎ＝±１を満たすようにδ（ｎ），ｎ＝±１を設定する。 （１ｃ）ミウラ型変換 【0055】 【数16】 ｑｋ（ｎ）＝η（ｎ）＊（１＋δ（ｎ）ｕ２ｋ－２（ｎ））＊（１＋δ（ｎ）ｕ２ｋ－１（ｎ））， ｋ＝１，２，・・・，ｍ， ｅｋ（ｎ）＝δ（ｎ）＊ｕ２ｋ－１（ｎ）＊ｕ２ｋ（ｎ） ｋ＝１，２，・・・，ｍ－１， を行う。 【0056】 ｑｄ型変換では見られない離散Ｌｏｔｋａ－Ｖｏｌｔｅｒｒａ型変換の大きな特徴は、任意パラメータを持つことである。すなわち、λｊ＝１／δ（０）－１／δ（±１）を満たす範囲でδ（ｎ）の値を任意に設定できる。δ（ｎ）を変動させると補助変数ｕｋ（ｎ）の値も変化するが、桁落ちによる誤差や数値不安定が発生するかどうかは事前に判定できる。この判定は、ｉｆ文によって実装されてもよい。この場合は、δ（ｎ）を再設定後に再度計算すればよい。また、ｕｋ（±１）が求まればミウラ変換によってｑｋ（±１）およびｅｋ（±１）が独立に計算されるので、誤差が伝播しないという性質を持つ。なお、ミウラ型逆変換をミウラ変換、ミウラ型変換を逆ミウラ変換、ｓｄＬＶｖｓ変換をｓｔＬＶｖ変換、ｒｄＬＶｖ変換をｒＬＶｖ変換と呼び変えても問題はない。 【0057】 入力および出力の変数の対応は、Ｐａｒｌｅｔｔの手法と同じである。メモリ消費を抑えるために、補助変数ｕｋ（ｎ）のための配列は必ずしも用意する必要はない。一方、１＋δ（０）ｕ（０）のためのメモリ領域を確保し、（１ａ）～（１ｃ）のステップにまたがってこの値を利用することで、メモリ消費を抑え、計算量を低減することができる。これにより、誤差も低減される。このように、１．３で説明する方法を利用して誤差を低減すること等、さらなる工夫を行い得る。 【0058】 上記の方法でＢＴＢ－λｊＩのコレスキー分解を求めた後は、Ｐａｒｌｅｔｔらと同じ方法を用いて特異ベクトルｖｊを計算する。そうすると、ＢＴＢの固有値分解Λ＝ＶＢＴＢＴＢＶＢが得られる。ただしＶＢの列ベクトルがｖｊである。ＵＢ＝ＢＶＢΣ－１よりＵＢを求めると、Ｂの特異値分解Ｂ＝ＵＢΣＶＢＴも得られる。 【0059】 １．２．高精度特異値算出ｄＬＶｓルーチン 本発明によるｄＬＶｓルーチンは、基本的な枠組みではＤＬＡＳＱルーチンと同じであり、ＢＴＢの固有値（あるいはＢの特異値）のみが得られＢＴＢの固有ベクトル（あるいはＢの特異ベクトル）は求まらない。ただ、ルーチンの核となる部分に使用される計算法はＤＬＡＳＱと異なる。 【0060】 まず、ＤＬＡＳＱルーチンによる固有値の計算方法を説明する。ｑｄ法は固有ベクトルが計算できないが、１ステップは上２重対角行列の要素を更新するのみなので高速に特異値を求めることができる。かつ、計算量が少ないため丸め誤差の蓄積が抑えられ高精度計算が期待できる。ここで、Ｙ（ｎ）を対称な３重対角行列 【0061】 【数17】 とする。また、変換Ｙ（ｎ）→Ｙ（ｎ＋１）がｑｄ法の漸化式 ｑｋ（ｎ＋１）＝ｑｋ（ｎ）＋ｅｋ（ｎ）－ｅｋ－１（ｎ＋１），ｋ＝１，２，・・・，ｍ， ｅｋ（ｎ＋１）＝ｅｋ（ｎ）＊ｑｋ＋１（ｎ）／ｑｋ（ｎ＋１），ｋ＝１，２，・・・，ｍ－１， ｅ０（ｎ）≡０，ｅｍ（ｎ）≡０，ｎ＝０，１，・・・ によって実現されるとする。そのとき、適切なＲ（ｎ）が存在してＹ（ｎ＋１）＝Ｒ（ｎ）Ｙ（ｎ）（Ｒ（ｎ））－１となることをＲｕｔｉｓｈａｕｓｅｒは発見している。これは、Ｙ（ｎ＋１）とＹ（ｎ）の固有値が同じであることを意味し、すなわち、上記の漸化式によって固有値保存変形がなされる。この変形を繰り返せば非対角成分が０に近づき、Ｙ（ｎ）が対角化されることも証明されている。よってＹ（０）＝ＢＴＢとするとｌｉｍ ｎ→∞Ｙ（ｎ）＝ｄｉａｇ（λ１，λ２，・・・，λｍ）となる。ただし、λｋはＢＴＢの固有値である。さらに、λｋの正の平方根をとるとＢの特異値σｋが得られる。 【0062】 ＬＡＰＣＡＫの特異値計算ルーチンＤＬＡＳＱでは、ｑｄ法のｄｉｆｆｅｒｅｎｔｉａｌ型が採用されている。これはｄｑｄ法と呼ばれる。１ステップは 【0063】 【数18】 で与えられ、変換）Ｙ（ｎ）→Ｙ（ｎ＋１）が実現される。ｄｉｆｆｅｒｅｎｔｉａｌ型は減算を含まないため、ｑｄ法では起こりうる桁落ち誤差の心配もない。さらに高速版として原点シフト付きｄｑｄ（ｄｑｄｓ）法も組み込まれている。高速化に伴う計算量の減少で丸め誤差の発生も抑えることができる。ただし、収束するかどうかの保証はない。原点シフト版の１ステップは 【0064】 【数19】 で与えられる。ｓは原点シフト量であり、Ｇｅｒｓｇｏｒｉｎ型境界条件により見積もられたＹ（ｎ）の最小固有値とする。よって、ＤＬＡＳＱルーチンのメインループの核となる部分では最小固有値の見積もりを行い、ｓ＞ε（ε：非常に小さいとみなせる正の数）ならばｄｑｄｓ法が、そうでなければｄｑｄ法が用いられる。メインループの他の部分（行列の分割・減次）はほとんどＤＢＤＳＱＲルーチンと変わらない。すなわち、ＤＢＤＳＱＲルーチンに含まれるＱＲ法をｄｑｄ法に、原点シフトＱＲ法をｄｑｄｓ法に置き換えたものがＤＬＡＳＱルーチンとなる。（ＱＲ法とｄｑｄ法、原点シフト付きＱＲ法とｄｑｄｓ法で扱う変数が異なるので少しだけ判定式が異なる部分もあるが本質的に同じである。） 次に本発明によるｄＬＶｓルーチンを説明する。ｄＬＶｓルーチンではｄｑｄ法のかわりにｄＬＶ法を、ｄｑｄｓ法のかわりにｄＬＶｓ法を採用する。ここで、ｄＬＶｓルーチンはＤＬＡＳＱルーチンとは違って収束することが保証されている。その１ステップは以下で与えられる。 【0065】 【数20】 ＤＬＡＳＱルーチン同様、メインループの核の部分では最小特異値の見積もりによって原点シフト量ｓを決定し、ｄＬＶ法を使うかｄＬＶｓ法を使うかを決定する。ただし、補助変数ｖｋ（ｎ）を求めたあとでｓを決定する点はＤＬＡＳＱルーチンと異なる。 【0066】 １．３．本発明の方法による誤差の低減 人間が理想的な状況、すなわち、無限桁の計算をいくらでもできるとすると、Ｐａｒｌｅｔｔらの手法、本発明による特異ベクトル計算ルーチン（１．１参照）のどちらを使っても問題ない。しかしながら、コンピュータで計算を行う場合は注意が必要である。有限桁の計算しか行えないコンピュータ上では、数学的に正しい計算法を使ったとしても必ずしも正しい結果が得られる訳ではない。そればかりかいつまでたっても計算が終了しないといった思わぬ数値的な問題が発生する場合もある。 【0067】 コンピュータ計算による誤差には、丸め誤差および桁落ちによる誤差などが知られている。丸め誤差単独では、高々有効桁の最後の桁が真値と比べて異なる程度で大きな誤差にはならない。また、指数部が異なる２つの実数の加算、乗算、除算の少なくとも１つの演算を行えばやはり丸め誤差が生じるが、それ以上の誤差は発生しない。さらに、このような丸め誤差を発生するような操作が繰り返されても、丸めモードがｎｅａｒ（四捨五入）ならば、一方的に切り上げられたりあるいは切り捨てられたりして誤差が極端に蓄積することは少ない。よって、多くの数値計算法は加算、乗算、除算の少なくとも１つの演算によって発生する丸め誤差を特別注意することは少なく、ｄＬＶｓルーチンによる特異値計算でも結果的に丸め誤差は一様に増大しない。 【0068】 問題となるのは、同符号の実数の減算および異符号の実数の加算により生じる、桁落ちである。桁落ちによる誤差で値が０となった後、その値による除算を行うと、０が分母にくるような不定形となり計算不可能となる。こうなるといつまでたっても計算が終了しない。減算→除算と計算する部分がＰａｒｌｅｔｔの手法、本発明による方法ともに存在するので、桁落ち誤差には十分に注意する必要がある。 【0069】 本発明による計算方法では、上述の桁落ちによる誤差を含んでいるかどうかは減算によって得られた値が小さいかどうかで判断できる。Ｐａｒｌｅｔｔの手法の場合、図２のＤＯ文中のｑ１［ｋ］、ｑ２［ｋ］をチェックすればよい。ところが、Ｐａｒｌｅｔｔの手法は桁落ち誤差を含むことが分かったとしても、それを回避することはできない。なぜならば、初期値として｛ｑ０[ｋ]，ｅ０[ｋ]｝が与えられると、λは一意的に決定され、｛ｑ１[ｋ]，ｅ１[ｋ]｝および｛ｑ２[ｋ]，ｅ２[ｋ]｝も一意的に導出されるためである。すなわち、任意パラメータを持たない自由度のない計算法であるためである。 【0070】 それに対して本発明による方法は自由に設定できるパラメータη１（＝１／δ（０））をもつため、補助変数ｕｋ（ｎ）の値を様々に変化させることができる（図９参照）。すなわち、様々な経路で｛ｑ１[ｋ]，ｅ１[ｋ]｝および｛ｑ２[ｋ]，ｅ２[ｋ]｝を計算することができる。よって、桁落ちが発生する場合も回避できる。図６～図８の条件判定によって桁落ちの影響をチェックし、減算によって得られた値の絶対値が小さな数εより大きいという条件が満たされなければ、パラメータη１の設定に戻るというものである。この処理は、条件が満たされるまで繰り返される。なお、精度よりも高速性を重視する場合は、数回条件が満たさなければ（ｑ１[ｋ]＝０あるいはｑ２[ｋ]＝０ならば）、例外処理（文献ＰＤＦ参照）を行ってもよい。 【0071】 ２．本発明による特異値分解装置 特異値分解を実行する装置を説明する。特異値分解を実行する装置は、例えば、図１に示される各ステップの機能を有する装置である。特異値分解を実行する装置の１つの実施形態は、特異値分解プログラムを実行するコンピュータである。そのプログラムは、例えば、図１に示される各ステップをＣＰＵに実行させるプログラムである。 【0072】 図１０は、図１の処理をコンピュータプログラムとして実行するコンピュータ１０を示す。コンピュータ１０は、ＣＰＵ１００１と、メモリ１００２と、入力インターフェース部１００３と、出力インターフェース部１００４と、バス１００５とを含む。 【0073】 ＣＰＵ１００１は、プログラムを実行する。例えば、そのプログラムは、図１の処理を実行するプログラムである。そのプログラムおよびそのプログラムの実行に必要なデータは、例えば、メモリ１００２に格納される。そのプログラムは、任意の様態でメモリ１００２に含まれ得る。例えば、メモリ１００２が書き換え可能なメモリである場合には、コンピュータ１００２の外部からプログラムをローディングすることにより、そのプログラムをメモリ１００２に格納してもよい。あるいは、メモリ１００２が読み出し専用メモリである場合には、メモリ１００２に焼き付ける形式でそのプログラムをメモリ１００２に格納してもよい。 【0074】 さらに、入力インターフェース部１００３は、特異値分解を実行する対象である行列Ａを外部から受け取るためのインターフェースとして機能する。出力インターフェース部１００４は、計算結果を出力するためのインターフェースとして機能する。バス１００５は、コンピュータ１０内の構成要素１００１～１００４を相互に接続するために使用される。 【0075】 ３．本発明による特異値分解方法の応用 本発明による特異値分解方法は、様々な分野に応用可能である。以下に、２次元画像から３次元画像へ復元する画像処理への応用例、ならびに、文書検索への応用例を示す。しかし、これら２つの応用例は例示に過ぎず、本発明による特異値分解方法の応用は、これら２つの応用例に限定されない。本発明による特異値分解方法の応用は、特異値分解を利用するあらゆる分野に適用可能である。 【0076】 ３．１．２次元から３次元へ復元する画像処理への応用 図１１を参照する。図１１は、本発明による特異値分解方法を利用した物体の複数の２次元画像を３次元画像へ復元する画像処理の１つの実施形態を説明する。 【0077】 複数の２次元画像から３次元復元を行うために必要となるステップは、２次元画像から特徴点を抽出するステップと、特徴点データより形状（元の物体の特徴点の３次元座標データ）および回転（３次元データから特徴点データへの変換）に関するデータを計算するステップと、形状および回転のデータより可視化を行うステップとを包含する。 【0078】 ステップ１１０２において、２次元画像ｊ（ｊ＝１，・・・，ｍ、ｍは３以上の整数）から特徴点ｉ（ｉ＝１，・・・，ｎ、ｎは２以上の整数）の座標（ｘｉｊ，ｙｉｊ）を抽出する。取り扱う２次元画像は、弱中心射影画像であることが好ましい。このとき 【0079】 【数21】 が成り立つ。ここで、ｓｊは物体のスケールに相対するｊ番目の画像のスケール、ｒ１，ｊ，ｒ２，ｊはそれぞれ物体座標系に相対するｊ番目のカメラ座標系の回転行列の１番目と２番目の行ベクトル、［Ｘｉ，Ｙｉ，Ｚｉ］Ｔはｉ番目の点の３次元座標である。物体のスケールは１番目の画像のスケールと同じにし（ｓ１＝１）、物体の座標系の姿勢は１番目の画像のカメラ座標系と同じにする（ｒ１，１＝［１，０，０］Ｔ，ｒ２，１＝［０，１，０］Ｔ）。ｍ枚全ての画像におけるｎ個全ての点の座標を行列Ｄの要素として並べると、 Ｄ＝ＭＳ が得られる。 ただし、 【0080】 【数22】 である。ＭとＳの形から分かるように、Ｄのランクは３である。いま、ステップ１１０２よりＤが与えられている。以下、回転に関するデータＭおよび形状Ｓを求める。 【0081】 そこで、行列Ｄの特異値分解 Ｄ＝ＵΣＶＴ を考える。ここで、Σは特異値を大小順に対角線上に並べたもので、ＵおよびＶはそれぞれ左および右直交行列である。ステップ１１０４では、行列Ｄの特異値分解を計算するために、行列Ｄを上２重対角化して上２重対角化行列Ｂを得る。ステップ１１０６において、行列Ｂ、すなわちＤの特異値を計算する。ステップ１１０６では、本発明によるｄＬＶｓルーチン、ＤＬＡＳＱルーチン、その他の特異値計算ルーチン、あるいはその組み合わせを利用し得る。このとき、画像のデジタル誤差のため、ゼロでない特異値は３つ以上出てくる。しかし、４番目以降の特異値はノイズによるもので、最初の３つの特異値と比べて格段に小さい。 【0082】 そこで、ステップ１１０８において、最初の３つの特異値に対して特異ベクトルを計算する。ステップ１１０８は、図１のステップ１１０～１２０を利用し得る。採用する３個のベクトルをまとめると、次式となる。 【0083】 Ｄ’＝Ｌ’Σ’Ｒ’Ｔ＝Ｍ’Ｓ’ ただし、Ｍ’＝Ｌ’（Σ’）１／２、Ｓ’＝（Σ’）１／２Ｒ’Ｔ、Ｄ’は∥Ｄ－Ｄ’∥を最小にするランク３の行列である。 【0084】 次に、Ｄ’からＭおよびＳを求めたいが、その組合せは唯一ではない。なぜなら、任意の正則行列Ｃが Ｄ’＝（Ｍ’Ｃ）（Ｃ－１Ｓ’） を満たすからである。そこで、上式におけるＣをＭ＝Ｍ’Ｃを満たすように決める。Ｃは下記の式を満たす。 【0085】 【数23】 Ｅ＝ＣＣＴとすると、上式からＥの６つの要素に関する２ｍ＋１個の線形方程式が得られる。ｍ≧３であるので、Ｅの要素を一意に決めることができる。ステップ１１１０において、行列Ｅを求める。 【0086】 次に、ステップ１１１２において、ＥからＣを求める。Ｃの自由度（９）はＥの自由度（６）より多い。そこで、条件ｒ１ｊ＝［１，０，０］Ｔ，ｒ２ｊ＝［０，１，０］Ｔを加えれば、Ｃを決めることができる。このとき２つの解（Ｎｅｃｋｅｒ Ｒｅｖｅｒｓａｌ）が出る。 【0087】 次に、ステップ１１１４において、Ｍ＝Ｍ’ＣおよびＳ＝Ｃ－１Ｓ’より、回転に関するデータＭおよび形状Ｓが決まる。 【0088】 ＤＢＤＳＱＲルーチンでは上記のようにＵ、Ｖ（ｎ個のベクトル）を求めた後、３個のベクトル（ｎ－３個は不要）を採用する。ここでは、原則的にはｎ個のベクトルを求めるが、例外的にｎ個より少ない数のベクトルを求めてもよい。この実施形態では３個のベクトルだけ求めればよい。すなわち、計算コストを削減することができる。 【0089】 ３．２．本発明による特異値分解方法を利用した文書検索方法 文書中からその文書の内容に関連する索引語を抽出し、索引語の重みを計算する処理の後、ベクトル空間モデルでは、この索引語の重みを要素とするベクトルで文書を表現する。いま、検索対象となる文書をＤ１，Ｄ２，・・・，Ｄｎとし、これら文書集合全体を通して全部でｍ個の索引語ｗ１，ｗ２，・・・，ｗｍがあるとする。このとき、文書Ｄｊは、次のようなベクトルで表現されることになる。これを文書ベクトルと呼ぶ。 【0090】 【数24】 ここで、ｄｉｊは索引語ｗｉの文書Ｄｊにおける重みである。また、文書集合全体は、次のようなｍ×ｎ行列Ｄによって表現することができる。 【0091】 【数25】 行列Ｄを索引語文書行列と呼ぶ。索引語文書行列の各列は文書に関する情報を表している文書ベクトルであるが、同様に、索引語文書行列の各行は索引語に関する情報を表しているベクトルであり、これを索引語ベクトルと呼ぶ。検索質問も、文書と同様に、索引語の重みを要素とするベクトルで表現することができる。検索質問文に含まれる索引語ｗｉの重みをｑｉとすると、検索質問ベクトルｑは次のように表されることになる。 【0092】 【数26】 実際の文書検索においては、与えられた検索質問文と類似した文書を見つけ出す必要があるが、検索質問ベクトルｑと各文書ベクトルｄｊの間の類似度を計算することにより行う。ベクトル間の類似度の定義としてはさまざまなものが考えられるが、文書検索においてよく用いられているものはコサイン尺度（２つのベクトルのなす角度）または内積である。 【0093】 【数27】 なお、ベクトルの長さが１に正規化（コサイン正規化）されている場合には、コサイン尺度と内積とは一致する。 【0094】 図１２を参照する。図１２は、本発明による特異値分解装置を利用した、文書検索方法の１つの実施形態を説明する。 【0095】 ステップ１２０２において、質問ベクトルｑを受け取る。次に、ステップ１２０４では、行列Ｄの特異値分解を計算するために、行列Ｄを上２重対角化して上２重対角化行列Ｂを得る。ステップ１２０６では、行列Ｂ、すなわち行列Ｄの特異値を求める。ステップ１２０６では、本発明によるｄＬＶｓルーチン、ＤＬＡＳＱルーチン、その他の特異値計算ルーチン、あるいはその組み合わせを利用し得る。 【0096】 次に、Ｄの近似行列Ｄｋを使った検索を考える。ベクトル空間モデルでは、検索質問ベクトルｑと索引語文書行列Ｄ中の各文書ベクトルｄｊの間の類似度を計算することにより検索を行うが、ここではＤの代わりにＤｋを使う。ベクトル空間モデルでは、文書ベクトルの次元数は索引語の総数と等しい。したがって、検索対象となる文書の数が増えるに従い、文書ベクトルの次元数も増加する傾向にある。しかし、次元数が増加してくると、コンピュータのメモリによる制限や検索時間の増大などの問題が生じてくるばかりでなく、文書中に含まれる不必要な索引語がノイズ的な影響を及ぼし、検索精度を低下させてしまうという現象も起こってくる。潜在的意味インデキシング（ｌａｔｅｎｔ ｓｅｍａｎｔｉｃ ｉｎｄｅｘｉｎｇ；ＬＳＩ）は、高次元の空間にある文書ベクトルを低次元の空間へと射影することにより、検索精度の改善を図る技術である。高次元の空間では別々に扱われていた索引語が、低次元の空間では相互に関連を持ったものとして扱われる可能性もあるため、索引語の持つ意味や概念に基づく検索を行うことができる。たとえば、通常のベクトル空間モデルでは“ｃａｒ”という索引語と“ａｕｔｏｍｏｂｉｌｅ”という索引語はまったく別物であり、一方の索引語による質問ではもう片方の索引語を含んだ文書を検索することができない。しかし、低次元の空間ではこれらの意味的に関連した索引語は１つの次元に縮退することが期待できるため、“ｃａｒ”という検索質問によって“ｃａｒ”を含む文書ばかりでなく“ａｕｔｏｍｏｂｉｌｅ”を含む文書をも検索することが可能となる。潜在的意味インデキシングでは、特異債分解により高次元ベクトルの次元削減を行うが、これは基本的に多変量解析における主成分分析と等価である。 【0097】 行列ステップ１２０８では、ｋ＜ｒなるｋの値を選択する。ｋの値は、予め与えられていてもよいし、計算ごとに選択可能であってもよい。 【0098】 ステップ１２１０において、ステップ１２０６で計算した特異値のうち、大きい順に１番目からｋ番目のｋ個の特異値に対して、Ｄの特異ベクトルを計算する。すなわち、 Ｄｋ＝ＵｋΣＶｋＴ なるＵｋおよびＶｋを計算する。ここで、Ｕｋは、最初のｋ個の左特異ベクトルのみから構成されるｍ×ｋ行列であり、Ｖｋは、最初のｋ個の右特異ベクトルのみから構成されるｎ×ｋ行列であり、Σは、最初のｋ個の特異値のみから構成されるｋ×ｋ対角行列である。ステップ１２１０は、図１のステップ１１０～１２０を利用し得る。 【0099】 次に、ステップ１２１２において、行列Ｄｋと質問ベクトルｑとの類似度を計算する。いま、ベクトルｅｊをｎ次元の単位ベクトルとすると、Ｄｋのｊ番目の文書ベクトルはＤｋｅｊで表すことができる。文書ベクトルＤｋｅｊと検索質問ベクトルｑとの間の類似度計算は、 【0100】 【数28】 としてもよいが、別の定義を用いてもよい。上式では、ＤｋをＵｋ，Σｋ，Ｖｋから再構成する必要はなく特異値分解の結果から、直接、類似度を計算できることを示している。 数２９の中に現われるΣｋＶｋＴｅｊは、 ΣｋＶｋＴｅｊ＝ＵｋＴＤｋｅｊ と書き直すことができる。この式の右辺は、近似行列Ｄｋにおけるｊ番目の文書ベクトルの基底Ｕｋのもとでの座標（文書のｋ次元表現）を表している。同様に、数２９中のＵｋＴｑは、検索質問ベクトルｑの基底Ｕｋのもとでの座標（検索質問のｋ次元表現）である。 【0101】 ステップ１２１４において、ステップ１２１２において計算された類似度を基準に、検索結果を出力する。 【0102】 ＤＢＤＳＱＲルーチンではｒ個のベクトルを求めた後、ｋ個のベクトル（ｒ－ｋ個は不要）を採用するのに対して、この実施形態では、ｋ個のベクトルだけ求めればよい。すなわち、計算コストを削減することができる。 【0103】 以上のように、本発明の好ましい実施形態を用いて本発明を例示してきたが、本発明は、この実施形態に限定して解釈されるべきものではない。本発明は、特許請求の範囲によってのみその範囲が解釈されるべきであることが理解される。当業者は、本発明の具体的な好ましい実施形態の記載から、本発明の記載および技術常識に基づいて等価な範囲を実施することができることが理解される。本明細書において引用した特許、特許出願および文献は、その内容自体が具体的に本明細書に記載されているのと同様にその内容が本明細書に対する参考として援用されるべきであることが理解される。 【産業上の利用可能性】 【0104】 本発明は、任意の行列Ａの特異値分解を高速かつ高精度に実行することを可能にする方法、プログラム、および装置等として有用である。
図面	【図1】 【図2】 【図3】 【図4】 【図5】 【図6】 【図7】 【図8】 【図9】 【図10】 【図11】 【図12】