Specification :(In Japanese)行列の高速高精度特異値分解方法、プログラムおよび装置

Country	(In Japanese)日本国特許庁(JP)
Gazette	(In Japanese)特許公報(B2)
Patent Number	P4325877
Date of registration	Jun 19, 2009
Date of issue	Sep 2, 2009
Title of the invention, or title of the device	(In Japanese)行列の高速高精度特異値分解方法、プログラムおよび装置
IPC (International Patent Classification)	G06F  17/16        (2006.01) G06T   1/00        (2006.01) G06F  17/30        (2006.01)
FI (File Index)	G06F 17/16 K G06T 1/00 315 G06F 17/30 220Z G06F 17/30 350C
Number of claims or invention	6
Total pages	39
Application Number	P2006-514122
Date of filing	Jun 1, 2005
International application number	PCT/JP2005/010084
International publication number	WO2005/119507
Date of international publication	Dec 15, 2005
Application number of the priority	2004166437
Priority date	Jun 3, 2004
Claim of priority (country)	(In Japanese)日本国(JP)
Date of request for substantive examination	Dec 1, 2006
Patentee, or owner of utility model right	(In Japanese)【識別番号】503360115 【氏名又は名称】独立行政法人科学技術振興機構
Inventor, or creator of device	(In Japanese)【氏名】中村 佳正 【氏名】岩▲崎▼ 雅史 【氏名】阪野 真也
Representative	(In Japanese)【識別番号】100078282、【弁理士】、【氏名又は名称】山本 秀策 【識別番号】100062409、【弁理士】、【氏名又は名称】安村 高明 【識別番号】100113413、【弁理士】、【氏名又は名称】森下 夏樹
Examiner	(In Japanese)【審査官】鳥居 稔
Document or reference	(In Japanese)特開平10－031747（JP，A） 特開2002－202983（JP，A） 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 ，2004年 3月 1日，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，2004年 2月27日，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年12月17日，第2004巻，13-18 中村 佳正 Yoshimasa NAKAMURA，高機能化を目指した非線形システム制御 Nonlinear Engineering and Control，システム／制御／情報 第43巻 第11号 SYSTEMS,CONTROL AND INFORMATION，日本，システム制御情報学会 The Institute of Systems,Control and Information Engineers，1999年11月15日，第43巻，8-15 中村 佳正 Yoshimasa NAKAMURA，非線形現象の不思議 ソリトン理論と数値計算法，電子情報通信学会誌 第80巻 第11号 THE JOURNAL OF THE INSTITUTE OF ELECTRONICS,INFORMATION AND COMMUNICATION ENGINEERS，日本，社団法人電子情報通信学会 THE INSTITUTE OF ELECTRONICS,INFORMATION AND COMMUNICATION ENGINEERS，1997年11月25日，第80巻，1143-1146
Field of search	G06F 17/16 G06F 17/30 G06T 1/00 IEEE Xplore JSTPlus(JDreamII)
Scope of claims	(In Japanese)【請求項7】 物体の複数の2次元画像から3次元画像を復元する方法であって、 2次元画像j（j＝1，・・・，m、mは3以上の整数）における特徴点i（i＝1，・・・，n、nは2以上の整数）の座標dij（xij，yij）を抽出するステップと、 該特徴点の2次元座標dij（xij，yij）より、該特徴点の3次元座標si（Xi，Yi，Zi）および2次元座標から3次元座標への変換を表す行列Mを計算するステップと を包含し、 該特徴点の3次元座標si（Xi，Yi，Zi）および2次元座標から3次元座標への変換を表す行列Mを計算するステップは、 行列Dに対して上2重対角化を行い、行列Dの上2重対角行列Bを求めるステップであって、行列Dは、 【数1】 と定義される、ステップと、 行列Dの特異値として行列Bの特異値σj（σ1≧σ2≧・・・≧σr＞0、rは、行列Dのランクに等しい）を求めるステップと、 σ1、σ2およびσ3に対する行列Dの特異ベクトルを求めるステップと、 M＝M’Cなる行列Cに対して、E＝CCTを満たす行列Eを計算するステップであって、M’＝L’（Σ’）1／2、Σ’はσ1、σ2およびσ3を対角要素に持ちそれ以外の要素が0である3×3行列、L’はσ1、σ2およびσ3に対応するDの特異ベクトルを左から順に並べた行列である、ステップと、 行列Eから、行列Cを計算するステップと、 行列Cから、該3次元座標si（Xi，Yi，Zi）および該変換を表す行列Mを計算するステップと を含み、 該σ1、σ2およびσ3に対する行列Dの特異ベクトルを求めるステップは、ミウラ型逆変換と、sdLVvs変換と、rdLVvs変換と、ミウラ型変換とを用いて行列BTB－σj2Iに対してツイスティッド分解を実行することにより、行列BTBを対角化するステップを含み、Iは単位行列である、方法。 【請求項8】 与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する方法であって、 該キーワードに対応する質問ベクトルqを受け取るステップと、 該文書に対応する索引語文書行列Dに対して上2重対角化を行い、行列Dの上2重対角行列Bを求めるステップと、 行列Dの特異値として行列Bの特異値σj（σ1≧σ2≧・・・≧σr＞0、rは、行列Dのランクに等しい）を求めるステップと、 k＜rなるkを選択するステップと、 行列Dの擬似行列Dkを計算するステップであって、行列Dkは、σ1，σ2，・・・，σkを対角要素としそれ以外の要素が0である行列Σk、σ1，σ2，・・・，σkに対応する特異ベクトルのみを左から順に並べた左右直交行列Uk，Vkを用いて、Dk＝UkΣkVkTと定義される、ステップと、 行列Dkと質問ベクトルqとの類似度を計算するステップと、 該計算された類似度を基準に、検索結果を出力するステップと を包含し、 該行列Dkの左右直交行列Uk，Vkを求めるステップは、ミウラ型逆変換と、sdLVvs変換と、rdLVvs変換と、ミウラ型変換とを用いて行列BTB－σj2I（j＝1，2，・・・，k）に対してツイスティッド分解を実行することにより、行列BTBを対角化するステップを含み、Iは単位行列である、方法。 【請求項10】 物体の複数の2次元画像から3次元画像を復元する画像復元処理をコンピュータに実行させるプログラムであって、 該画像復元処理は、 2次元画像j（j＝1，・・・，m、mは3以上の整数）における特徴点i（i＝1，・・・，n、nは2以上の整数）の座標dij（xij，yij）を抽出するステップと、 該特徴点の2次元座標dij（xij，yij）より、該特徴点の3次元座標si（Xi，Yi，Zi）および2次元座標から3次元座標への変換を表す行列Mを計算するステップと を包含し、 該特徴点の3次元座標si（Xi，Yi，Zi）および2次元座標から3次元座標への変換を表す行列Mを計算するステップは、 行列Dに対して上2重対角化を行い、行列Dの上2重対角行列Bを求めるステップであって、行列Dは、 【数2】 と定義される、ステップと、 行列Dの特異値として行列Bの特異値σj（σ1≧σ2≧・・・≧σr＞0、rは、行列Dのランクに等しい）を求めるステップと、 σ1、σ2およびσ3に対する行列Dの特異ベクトルを求めるステップと、 M＝M’Cなる行列Cに対して、E＝CCTを満たす行列Eを計算するステップであって、M’＝L’（Σ’）1／2、Σ’はσ1、σ2およびσ3を対角要素に持ちそれ以外の要素が0である3×3行列、L’はσ1、σ2およびσ3に対応するDの特異ベクトルを左から順に並べた行列である、ステップと、 行列Eから、行列Cを計算するステップと、 行列Cから、該3次元座標si（Xi，Yi，Zi）および該変換を表す行列Mを計算するステップと を含み、 該σ1、2およびσ3に対する行列Dの特異ベクトルを求めるステップは、ミウラ型逆変換と、sdLVvs変換と、rdLVvs変換と、ミウラ型変換とを用いて行列BTB－σj2Iに対してツイスティッド分解を実行することにより、行列BTBを対角化するステップを含み、Iは単位行列である、プログラム。 【請求項11】 与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する文書検索処理をコンピュータに実行させるプログラムであって、 該文書検索処理は、 該キーワードに対応する質問ベクトルqを受け取るステップと、 該文書に対応する索引語文書行列Dに対して上2重対角化を行い、行列Dの上2重対角行列Bを求めるステップと、 行列Dの特異値として行列Bの特異値σj（σ1≧σ2≧・・・≧σr＞0、rは、行列Dのランクに等しい）を求めるステップと、 k＜rなるkを選択するステップと、 行列Dの擬似行列Dkを計算するステップであって、行列Dkは、σ1，σ2，・・・，σkを対角要素としそれ以外の要素が0である行列Σk、σ1，σ2，・・・，σkに対応する特異ベクトルのみを左から順に並べた左右直交行列Uk，Vkを用いて、Dk＝UkΣkVkTと定義される、ステップと、 行列Dkと質問ベクトルqとの類似度を計算するステップと、 該計算された類似度を基準に、検索結果を出力するステップと を包含し、 該行列Dkの左右直交行列Uk，Vkを求めるステップは、ミウラ型逆変換と、sdLVvs変換と、rdLVvs変換と、ミウラ型変換とを用いて行列BTB－σj2I（j＝1，2，・・・，k）に対してツイスティッド分解を実行することにより、行列BTBを対角化するステップを含み、Iは単位行列である、プログラム。 【請求項13】 物体の複数の2次元画像を3次元画像に復元する装置であって、 m枚（mは3以上の整数）の2次元画像を受け取る手段と、 2次元画像j（j＝1，・・・，m）における特徴点i（i＝1，・・・，n、nは2以上の整数）の座標dij（xij，yij）を抽出する手段と、 該特徴点の2次元座標dij（xij，yij）より、該特徴点の3次元座標si（Xi，Yi，Zi）および2次元座標から3次元座標への変換を表す行列Mを計算する手段と を備え、 該特徴点の3次元座標si（Xi，Yi，Zi）および2次元座標から3次元座標への変換を表す行列Mを計算する手段は、 行列Dに対して上2重対角化を行い、行列Dの上2重対角行列Bを求める手段であって、行列Dは、 【数3】 と定義される、手段と、 行列Dの特異値として行列Bの特異値σj（σ1≧σ2≧・・・≧σr＞0、rは、行列Dのランクに等しい）を求める手段と、 σ1、σ2およびσ3に対する行列Dの特異ベクトルを求める手段と、 M＝M’Cなる行列Cに対して、E＝CCTを満たす行列Eを計算する手段であって、M’＝L’（Σ’）1／2、Σ’はσ1、σ2およびσ3を対角要素に持ちそれ以外の要素が0である3×3行列、L’はσ1、σ2およびσ3に対応するDの特異ベクトルを左から順に並べた行列である、手段と、 行列Eから、行列Cを計算する手段と、 行列Cから、該3次元座標si（Xi，Yi，Zi）および該変換を表す行列Mを計算する手段と を含み、 該σ1、σ2およびσ3に対する行列Dの特異ベクトルを求める手段は、ミウラ型逆変換と、sdLVvs変換と、rdLVvs変換と、ミウラ型変換とを用いて行列BTB－σj2Iに対してツイスティッド分解を実行することにより、行列BTBを対角化する手段を含み、Iは単位行列である、装置。 【請求項14】 与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する装置であって、 該キーワードに対応する質問ベクトルqを受け取る手段と、 該文書に対応する索引語文書行列Dに対して上2重対角化を行い、行列Dの上2重対角行列Bを求める手段と、 行列Dの特異値として行列Bの特異値σj（σ1≧σ2≧・・・≧σr＞0、rは、行列Dのランクに等しい）を求める手段と、 k＜rなるkを選択する手段と、 行列Dの擬似行列Dkを計算する手段であって、行列Dkは、σ1，σ2，・・・，σkを対角要素としそれ以外の要素が0である行列Σk、σ1，σ2，・・・，σkに対応する特異ベクトルのみを左から順に並べた左右直交行列Uk，Vkを用いて、Dk＝UkΣkVkTと定義される、手段と、 行列Dkと質問ベクトルqとの類似度を計算する手段と、 該計算された類似度を基準に、検索結果を出力する手段と を包含し、 該行列Dkの左右直交行列Uk，Vkを求める手段は、ミウラ型逆変換と、sdLVvs変換と、rdLVvs変換と、ミウラ型変換とを用いて行列BTB－σj2I（j＝1，2，・・・，k）に対してツイスティッド分解を実行することにより、行列BTBを対角化する手段を含み、Iは単位行列である、装置。
Detailed description of the invention	(In Japanese)【技術分野】 【0001】 本発明は、任意の行列Aに対して高速かつ高精度に特異値分解を実行する方法、プログラムおよび装置に関し、より詳細には、行列の特異値を計算し、計算された特異値に対する特異ベクトルを計算することによって、高速かつ高精度に特異値分解を実行する方法、プログラムおよび装置に関する。 【背景技術】 【0002】 現在、コンピュータ上で特異値分解を行うために用いられる標準的な方法は、国際的な数値計算ライブラリLAPACKにより公開されているDBDSQRルーチンである。DBDSQRルーチンは、QR法に基づいて作られたものであり、前処理、メインループ、後処理の3つの部分に分けることができる。前処理では、特異値の上限および下限を計算し、収束判定に用いる精度を計算する。メインループでは、QR法を繰り返しながら、除々に行列の分割および減次を行い、最終的に行列サイズが1となった時点でループを終了する。途中で2×2行列のブロックが現れたならば別処理を行う。後処理では、特異値として計算された値が負ならば正に直し、必要ならば特異ベクトルをスカラー倍する。特異値を降順に、それに対応するように特異ベクトルも並べ替える。DBDSQRは、計算量が非常に多く、特に大規模問題では計算時間の増大が避けられない。DBDSQRルーチンでは、特異値および特異ベクトルが同時に計算される。また、LAPACKには、特異値を計算するDLASQルーチン、および、計算された特異値を用いて特異ベクトルを計算する際に用いる対角化のためのDSTEGRルーチンがある。DLASQルーチンは、高速高精度に特異値を求めることができるが、特異ベクトルを計算できない。DSTEGRルーチンは、その数値的な欠点から、実際に特異値分解に使用することは難しい。 【0003】 LAPACKによるDBDSQRルーチンを例として、従来の特異値分解方法を説明する。DBDSQRルーチンでは、l1×l2の一般行列Aを特異値分解するために、まず、行列Aに対してHouseholder（ハウスホルダー）変換を適用する。すなわち、行列Aは、直交行列UA、VAを用いて、 【0004】 【数4】 と表わすことができる。このとき得られるBを、上2重対角行列という。ここで、Aの特異値＝Bの特異値であることに注意する。このように、一般行列Aに対する特異値分解問題を上2重対角行列Bに対する特異値分解問題 B＝UBΣVBT 行列UB、VBは、それぞれ左および右直交行列 Σ≡diag（σ1，σ2，・・・，σm） σ1≧σ2≧・・・≧σm≧0 σjはBの特異値 に置き換える。 【0005】 次に、行列BTBを考える。この行列BTBの対角化 Λ＝VTBTBV Λ≡diag（λ1，λ2，・・・，λm） λ1≧λ2≧・・・≧λm≧0 V≡（v1，v2，・・・，vm） λjはBTBの固有値 vjは固有値λjに対する固有ベクトル を行う。ここで、一般に、以下のことが成り立つ。（1）BTBは対称な3重対角行列である。（2）BTBの固有値は全て正であり、Bの特異値σj（σ1≧σ2≧・・・≧σm≧0）は、BTBの固有値λj（λ1≧λ2≧・・・≧λm≧0）の正の平方根に等しい。（3）VB＝V、すなわち、BTBの固有ベクトルvjは、Bの右特異ベクトルに等しい。従って、行列BTBの対角化が求まると、（2）よりΛ＝Σ2であることから、Σが求まり、さらに（3）より左直交行列UB＝BVBΣ－1＝BVΣ－1が求まるので、Bが特異値分解される。すなわち、Bの特異値分解を、BTBの対角化の問題に置き換えることができる。この原理は、m個全ての特異値および特異ベクトルを求める場合だけでなく、少なくとも1つの特異値および特異ベクトルを求める場合にも適用され得る。 【0006】 以上のように、一般行列Aの特異値分解は、対称な3重対角行列BTBの対角化の問題を含む。 【発明の開示】 【発明が解決しようとする課題】 【0007】 DBDSQRルーチンでは、計算量が非常に多いため、特に大規模問題では、高速に特異値分解を実行することが困難であった。一方、DLASQルーチンは、高速かつ高精度に特異値を求めることができるが、DSTEGRルーチンは、場合によっては特異ベクトルを低精度で計算するため、DSTEGRルーチンでは、常に高精度に特異値分解を実行することは困難であった。 【0008】 本発明の目的の1つは、対称な3重対角行列BTBの対角化を高速かつ高精度で実行することにより、任意の行列Aの特異値分解を高速かつ高精度に実行することを可能にする方法、プログラム、および装置を提供することにある。ここで、行列Bを、行列Aをハウスホルダー変換することで得られる上2重対角行列とする。 【課題を解決するための手段】 【0009】 本発明では、ミウラ型逆変換、sdLVvs変換、rdLVvs変換およびミウラ型変換によるTwisted（ツイスティッド）分解を実行することにより、行列BTBを対角化する。ここで、行列Bを、行列Aと特異値が同じ上2重対角行列とする。行列Aから行列Bへの変換は、例えばハウスホルダー法で実行できる。本発明による行列BTBの対角化では、固有値および固有ベクトルを同時に求めるDBDSQRルーチンとは異なり、DSTEGRルーチンのように、まず固有値を計算し、次に計算された固有値を用いて固有ベクトルを求める。 【0010】 本発明により、コンピュータを用いて任意の行列Aに対して特異値分解を実行する方法であって、行列Aに対して上2重対角化を行い、行列Aの上2重対角行列Bを求めるステップと、行列Aの特異値として行列Bの少なくとも1つの特異値σを求めるステップと、σに対する行列Aの特異ベクトルを求めるステップとを包含し、該行列Aの特異ベクトルを求めるステップは、ミウラ型逆変換と、sdLVvs変換と、rdLVvs変換と、ミウラ型変換とを用いて行列BTB－σ2Iに対してツイスティッド分解を実行することにより、行列BTBを対角化するステップを含み、Iは単位行列である、方法が提供され、これにより上記目的が達成される。 【0011】 前記行列Aの特異ベクトルを求めるステップは、前記行列BTBを対角化するステップの後に、逆反復法を実行するステップを含んでもよい。 【0012】 前記行列Aの特異ベクトルを求めるステップは、前記逆反復法を実行するステップの後に、グラムシュミット法を実行するステップを含んでもよい。 【0013】 前記行列Bの特異値σを求めるステップは、dLVsルーチンを実行するステップを含んでもよい。 【0014】 前記行列Bの特異値σを求めるステップは、要求される計算時間および精度に応じて、dLVsルーチンもしくはDLASQルーチンのどちらを実行するのかを判定するステップを含んでもよい。 【0015】 前記行列Aに対して上2重対角化を行い、上2重対角行列Bを求めるステップは、上2重対角化法（たとえばハウスホルダー法）を実行するステップを含んでもよい。 【0016】 本発明により、物体の複数の2次元画像から3次元画像を復元する方法であって、2次元画像j（j＝1，・・・，m、mは3以上の整数）における特徴点i（i＝1，・・・，n、nは2以上の整数）の座標dij（xij，yij）を抽出するステップと、該特徴点の2次元座標dij（xij，yij）より、該特徴点の3次元座標si（Xi，Yi，Zi）および2次元座標から3次元座標への変換を表す行列Mを計算するステップとを包含し、該特徴点の3次元座標si（Xi，Yi，Zi）および2次元座標から3次元 座標への変換を表す行列Mを計算するステップは、行列Dに対して上2重対角化を行い、行列Dの上2重対角行列Bを求めるステップであって、行列Dは、 【0017】 【数5】 と定義される、ステップと、行列Dの特異値として行列Bの特異値σj（σ1≧σ2≧・・・≧σr＞0、rは、行列Dのランクに等しい）を求めるステップと、σ1、σ2およびσ3に対する行列Dの特異ベクトルを求めるステップと、M＝M’Cなる行列Cに対して、E＝CCTを満たす行列Eを計算するステップであって、M’＝L’（Σ’）1／2、Σ’はσ1、σ2およびσ3を対角要素に持ちそれ以外の要素が0である3×3行列、L’はσ1、σ2およびσ3に対応するDの特異ベクトルを左から順に並べた行列である、ステップと、行列Eから、行列Cを計算するステップと、行列Cから、該3次元座標si（Xi，Yi，Zi）および該変換を表す行列Mを計算するステップとを含み、該σ1、σ2およびσ3に対する行列Dの特異ベクトルを求めるステップは、ミウラ型逆変換と、sdLVvs変換と、rdLVvs変換と、ミウラ型変換とを用いて行列BTB－σj2Iに対してツイスティッド分解を実行することにより、行列BTBを対角化するステップを含み、Iは単位行列である、方法が提供され、これにより上記目的が達成される。 【0018】 本発明により、与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する方法であって、該キーワードに対応する質問ベクトルqを受け取るステップと、該文書に対応する索引語文書行列Dに対して上2重対角化を行い、行列Dの上2重対角行列Bを求めるステップと、行列Dの特異値として行列Bの特異値σj（σ1≧σ2≧・・・≧σr＞0、rは、行列Dのランクに等しい）を求めるステップと、k＜rなるkを選択するステップと、行列Dの擬似行列Dkを計算するステップであって、行列Dkは、σ1，σ2，・・・，σkを対角要素としそれ以外の要素が0である行列Σk、σ1，σ2，・・・，σkに対応する特異ベクトルのみを左から順に並べた左右直交行列Uk，Vkを用いて、Dk＝UkΣkVkTと定義される、ステップと、行列Dkと質問ベクトルqとの類似度を計算するステップと、該計算された類似度を基準に、検索結果を出力するステップとを包含し、該行列Dkの左右直交行列Uk，Vkを求めるステップは、ミウラ型逆変換と、sdLVvs変換と、rdLVvs変換と、ミウラ型変換とを用いて行列BTB－σj2I（j＝1，2，・・・，k）に対してツイスティッド分解を実行することにより、行列BTBを対角化するステップを含み、Iは単位行列である、方法が提供され、これにより上記目的が達成される。 【0019】 本発明により、任意の行列Aに対して特異値分解を実行する特異値分解処理をコンピュータに実行させるプログラムであって、該特異値分解処理は、行列Aに対して上2重対角化を行い、行列Aの上2重対角行列Bを求めるステップと、行列Aの特異値として行列Bの少なくとも1つの特異値σを求めるステップと、σに対する行列Aの特異ベクトルを求めるステップとを包含し、該行列Aの特異ベクトルを求めるステップは、ミウラ型逆変換と、sdLVvs変換と、rdLVvs変換と、ミウラ型変換とを用いて行列BTB－σ2Iに対してツイスティッド分解を実行することにより、行列BTBを対角化するステッ プを含み、Iは単位行列である、プログラムが提供され、これにより上記目的が達成される。 【0020】 本発明により、物体の複数の2次元画像から3次元画像を復元する画像復元処理をコンピュータに実行させるプログラムであって、該画像復元処理は、2次元画像j（j＝1，・・・，m、mは3以上の整数）における特徴点i（i＝1，・・・，n、nは2以上の整数）の座標dij（xij，yij）を抽出するステップと、該特徴点の2次元座標dij（xij，yij）より、該特徴点の3次元座標si（Xi，Yi，Zi）および2次元座標から3次元座標への変換を表す行列Mを計算するステップとを包含し、該特徴点 の3次元座標si（Xi，Yi，Zi）および2次元座標から3次元座標への変換を表す行列Mを計算するステップは、行列Dに対して上2重対角化を行い、行列Dの上2重対角行列Bを求めるステップであって、行列Dは、 【0021】 【数6】 と定義される、ステップと、行列Dの特異値として行列Bの特異値σj（σ1≧σ2≧・・・≧σr＞0、rは、行列Dのランクに等しい）を求めるステップと、σ1、σ2およびσ3に対する行列Dの特異ベクトルを求めるステップと、M＝M’Cなる行列Cに対して、E＝CCTを満たす行列Eを計算するステップであって、M’＝L’（Σ’）1／2、Σ’はσ1、σ2およびσ3を対角要素に持ちそれ以外の要素が0である3×3行列、L’はσ1、σ2およびσ3に対応するDの特異ベクトルを左から順に並べた行列である、ステップと、行列Eから、行列Cを計算するステップと、行列Cから、該3次元座標si（Xi，Yi，Zi）および該変換を表す行列Mを計算するステップとを含み、該σ1、σ2およびσ3に対する行列Dの特異ベクトルを求めるステップは、ミウラ型逆変換と、sdLVvs変換と、rdLVvs変換と、ミウラ型変換とを用いて行列BTB－σj2Iに対してツイスティッド分解を実行することにより、行列BTBを対角化するステップを含み、Iは単位行列である、プログラムが提供され、これにより上記目的が達成される。 【0022】 本発明により、与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する文書検索処理をコンピュータに実行させるプログラムであって、該文書検索処理は、該キーワードに対応する質問ベクトルqを受け取るステップと、該文書に対応する索引語文書行列Dに対して上2重対角化を行い、行列Dの上2重対角行列Bを求めるステップと、行列Dの特異値として行列Bの特異値σj（σ1≧σ2≧・・・≧σr＞0、rは、行列Dのランクに等しい）を求めるステップと、k＜rなるkを選択するステップと、行列Dの擬似行列Dkを計算するステップであって、行列Dkは、σ1，σ2 ，・・・，σkを対角要素としそれ以外の要素が0である行列Σk、σ1，σ2，・・・，σkに対応する特異ベクトルのみを左から順に並べた左右直交行列Uk，Vkを用いて、Dk＝UkΣkVkTと定義される、ステップと、行列Dkと質問ベクトルqとの類似度を計算するステップと、該計算された類似度を基準に、検索結果を出力するステップとを包含し、該行列Dkの左右直交行列Uk，Vkを求めるステップは、ミウラ型逆変換と、sdLVvs変換と、rdLVvs変換と、ミウラ型変換とを用いて行列BTB－σj 2I（j＝1，2，・・・，k）に対してツイスティッド分解を実行することにより、行列BTBを対角化するステップを含み、Iは単位行列である、プログラムが提供され、これにより上記目的が達成される。 【0023】 本発明により、任意の行列Aに対して特異値分解を実行する装置であって、行列Aを入力として受け取る手段と、行列Aに対して上2重対角化を行い、行列Aの上2重対角行列Bを計算する手段と、行列Aの特異値として行列Bの少なくとも1つの特異値σを計算する手段と、σに対する行列Aの特異ベクトルを求める手段とを備え、該行列Aの特異ベクトルを求める手段は、ミウラ型逆変換、sdLVvs変換およびrdLVvs変換、ならびにミウラ型変換を用いて行列BTB－σ2Iに対してツイスティッド分解を実行することにより、行列BTBを対角化する手段を含み、Iは単位行列である、装置が提供され、これにより上記目的が達成される。 【0024】 本発明により、物体の複数の2次元画像を3次元画像に復元する装置であって、m枚（mは3以上の整数）の2次元画像を受け取る手段と、2次元画像j（j＝1，・・・，m）における特徴点i（i＝1，・・・，n、nは2以上の整数）の座標dij（xij，yij）を抽出する手段と、該特徴点の2次元座標dij（xij，yij）より、該特徴点の3次元座標si（Xi，Yi，Zi）および2次元座標から3次元座標への変換を表す行列Mを計算する手段とを備え、該特徴点の3次元座標si（Xi，Yi，Zi）および2次元座標から3次元座標への変換を表す行列Mを計算する手段は、行列Dに対して上2重対角化を行い、行列Dの上2重対角行列Bを求める手段であって、行列Dは、 【0025】 【数7】 と定義される、手段と、行列Dの特異値として行列Bの特異値σj（σ1≧σ2≧・・・≧σr＞0、rは、行列Dのランクに等しい）を求める手段と、σ1、σ2およびσ3に対する行列Dの特異ベクトルを求める手段と、M＝M’Cなる行列Cに対して、E＝CCTを満たす行列Eを計算する手段であって、M’＝L’（Σ’）1／2、Σ’はσ1、σ2およびσ3を対角要素に持ちそれ以外の要素が0である3×3行列、L’はσ1、σ2およびσ3に対応するDの特異ベクトルを左から順に並べた行列である、手段と、行列Eから、行列Cを計算する手段と、行列Cから、該3次元座標si（Xi，Yi，Zi）および該変換を表す行列Mを計算する手段とを含み、該σ1、σ2およびσ3に対する行列Dの特異ベクトルを求める手段は、ミウラ型逆変換と、sdLVvs変換と、rdLVvs変換と、ミウラ型変換とを用いて行列BTB－σj2Iに対してツイスティッド分解を実行することにより、行列BTBを対角化する手段を含み、Iは単位行列である、装置が提供され、これにより上記目的が達成される。 【0026】 本発明により、与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する装置であって、該キーワードに対応する質問ベクトルqを受け取る手段と、該文書に対応する索引語文書行列Dに対して上2重対角化を行い、行列Dの上2重対角行列Bを求める手段と、行列Dの特異値として行列Bの特異値σj（σ1≧σ2≧・・・≧σr＞0、rは、行列Dのランクに等しい）を求める手段と、k＜rなるkを選択する手段と、行列Dの擬似行列Dkを計算する手段であって、行列Dkは、σ1，σ2，・・・，σkを対角要素としそれ以外の要素が0である行列Σk、σ1，σ2，・・・，σkに対応する特異ベクトルのみを左から順に並べた左右直交行列Uk，Vkを用いて、Dk＝UkΣkVkTと定義される、手段と、行列Dkと質問ベクトルqとの類似度を計算する手段と、該計算された類似度を基準に、検索結果を出力する手段とを包含し、該行列Dkの左右直交行列Uk，Vkを求める手段は、ミウラ型逆変換と、sdLVvs変換と、rdLVvs変換と、ミウラ型変換とを用いて行列BTB－σj2I（j＝1，2，・・・，k）に対してツイスティッド分解を実行することにより、行列BTBを対角化する手段を含み、Iは単位行列である、装置が提供され、これにより上記目的が達成される。 【発明の効果】 【0027】 本発明によれば、ミウラ型逆変換、sdLVvs変換、rdLVvs変換およびミウラ型変換によるツイスティッド分解を実行することにより、行列BTBが対角化される。それにより、DSTEGRルーチンと比較した場合、精度の点ではるかに優れ、実用化に耐え得る固有値計算を実現し、高精度の特異値分解が可能になる。また、逆反復法、Gram－Schmidt（グラムシュミット）法、本発明において提供するdLVsルーチンのうちの少なくとも1つを併用することにより、DBDSQRと比較して精度の点では同等になることもあるが、速度の点では圧倒的に高速である特異値分解を実現することができる。さらに、DBDSQRルーチンと比較すると、固有値および固有ベクトルを同時の求めるのか、または、計算された固有値を用いて固有ベクトルを求めるのかの計算フローの違いにより、以下で説明する2次元から3次元への画像復元問題および文書検索問題等においては、計算時間をさらに短縮することができる。 【図面の簡単な説明】 【0028】 【図1】図1は、本発明により提供される特異値分解のためのI－SVDルーチンの処理の手順を示すフローチャート 【図2】図2は、Parlettらの手法によるsqds変換およびrpqds変換の処理の手順を示すフローチャート 【図3】図3は、本発明によるミウラ型逆変換、sdLVvs変換、rdLVvs変換およびミウラ型変換の処理の手順を示すフローチャート 【図4】図4は、本発明によるミウラ型逆変換、sdLVvs変換、rdLVvs変換およびミウラ型変換の処理の手順を示すフローチャート 【図5】図5は、本発明によるミウラ型逆変換、sdLVvs変換、rdLVvs変換およびミウラ型変換の処理の手順を示すフローチャート。 【図6】図6は、本発明によるミウラ型逆変換、sdLVvs変換、rdLVvs変換およびミウラ型変換の処理の手順を示すフローチャート 【図7】図7は、本発明によるミウラ型逆変換、sdLVvs変換、rdLVvs変換およびミウラ型変換の処理の手順を示すフローチャート 【図8】図8は、本発明によるミウラ型逆変換、sdLVvs変換、rdLVvs変換およびミウラ型変換の処理の手順を示すフローチャート 【図9】図9は、本発明において誤差の低減をもたらす、Parlettの手法と本発明による方法との計算経路の違いを示す図 【図10】図10は、本発明による特異値分解装置の1つの実施形態であるコンピュータを示す図 【図11】図11は、本発明による特異値分解方法を利用した物体の複数の2次元画像を3次元画像へ復元する画像処理の1つの実施形態を示すフローチャート 【図12】図12は、本発明による特異値分解装置を利用した文書検索方法の1つの実施形態を示すフローチャート 【符号の説明】 【0029】 10 コンピュータ 1001 CPU 1002 メモリ 1003 入力インターフェース部 1004 出力インターフェース部 1005 バス 【発明を実施するための最良の形態】 【0030】 以下、図面を参照しながら、本発明の実施形態を説明する。 【0031】 1．特異値分解アルゴリズムI－SVDルーチン 図1を参照する。図1は、本発明により提供される特異値分解のためのI－SVDルーチンの処理の手順を示す。図1の処理は、コンピュータプログラムの形式で提供され得る。 【0032】 ステップ102において、一般行列Aに対してハウスホルダー法を用いて上2重対角化を行う。このとき利用される上2重対角化方法は、ハウスホルダー法であってもよいし、他の上2重対角化方法であってもよい。これにより、上2重対角行列Bが得られる。次に、行列Bの特異値σj（σ1≧σ2≧・・・≧σm≧0、mは、行列BTBの行列サイズに等しい）を計算する。ここで、本発明により提供される特異値計算ルーチンdLVsおよびLAPACKにより提供されるDLASQルーチンを利用し得る。なお、行列Bの特異値は、行列Aの特異値に等しく、かつ、行列BTBの固有値の正の平方根λj1／2に等しい。dLVsルーチンは、Bの特異値σjを高精度で計算する。dLVsルーチンの詳細は、後述される。dLVsで求められる特異値は、従来のLAPCKルーチン（DBDSQR、DLASQ等）と比較して、最も高精度である。速度に関しては、ほとんどの場合DLASQルーチンには劣るが、信頼性の高いDBDSQRルーチンの約2分の1の時間で計算できる(CPUがItanium2の場合はdLVsルーチンが最速)。そこで、ステップ104において、高精度性が必要とされるかどうかを判定する。最も高精度な特異値分解を行いたい場合には、ステップ106においてdLVsを、高精度性および高速性を両立したい場合には、ステップ108において従来のDLASQルーチンを、特異値計算ルーチンとして用いるのがよい。さらに、DLASQルーチンが有限回の演算で停止するかは定かでないが、dLVsルーチンは収束することが保証されているので、信頼性を重視するならば、dLVsルーチンを用いるのがよい。 【0033】 次に、ステップ110において、計算されたBの特異値σj（すなわち、Aの特異値に等しい）用いて、1．1において説明されるsdLVvs変換およびrdLVvs変換によるツイスティッド分解により、BTBの対角化を行う。本発明による行列BTBの対角化ルーチンの詳細は、1．1で説明する。さらに、ステップ112において、以下の1．1にて詳細が示される簡略化連立1次方程式 【0034】 【数8】 を求めると、行列Bの特異ベクトルvjが得られる（すなわち、BTBの対角化行列Vが得られる）。すなわち、BTBが対角化される。ステップ114において、Σ≡diag（σ1，σ2，・・・，σm）（ただし、σ1≧σ2≧・・・≧σm≧0）とすると、行列Bの右直交行列VB＝Vであるので、行列Bの左直交行列UB＝BVBΣ－1＝BVΣ－1からUBを得る。すなわち、行列Bが特異値分解される。ここでは、m個全ての特異値および特異ベクトルを計算したが、この原理は、少なくとも1つの特異値および特異ベクトルを計算する場合にも適用され得る。 【0035】 このとき得られる左右直交行列、すなわち特異ベクトル計算の結果は、DBDSQRルーチンによって得られる結果と比較して精度の点で劣る。そこで、ステップ116において、逆反復法を用いることにより、DBDSQRルーチンと同程度の精度を得ることができる。さらなる高精度性が必要とされる場合は、ステップ120において、グラムシュミット法による再直交化を行えば、DBDSQRよりも高精度での計算を実現し得る。I－SVDルーチンでは、速度に関しては、表2から、DBDSQRと比較して大幅に短縮され得ることが理解される。また、表1から理解されるように、I－SVDルーチンは、DSTEGRによる極度の精度悪化がみられる場合でも、高精度に特異値分解することができる。なお、表1、2の誤差および計算時間に関するデータは、逆反復法を用いているが、グラムシュミット法は用いていない場合の結果である。また、ステップ116および120を実行する順序は逆であってもよいし、あるいは、両ステップは、省略してもよい。 【0036】 このように、図1に示される各ステップの機能は、ソフトウェア（例えば、プログラム）によって実現される。しかし、本発明は、これに限定されない。図1に示される各ステップの機能をハードウェア（例えば、回路、ボード、半導体チップ）によって実現してもよいし、ソフトウェアとハードウェアとの組み合わせによって実現してもよい。 【0037】 以上により、I－SVDルーチンを用いると、従来技術と比較して高速かつ高精度の特異値分解が実行される。 【0038】 表1は、本発明による行列BTB対角化ルーチンを実行した後に、逆反復法を実行したときの従来技術および本発明の誤差を比較した表である。 【0039】 【表1】 表2は、本発明による行列BTB対角化ルーチンを実行した後に、逆反復法を実行したときの従来技術および本発明の計算時間を比較した表である。 【0040】 【表2】 また、1.3に、本発明の方法において、DSTEGRルーチンと比較して、誤差が低減される理由を説明する。 【0041】 1．1．本発明による行列BTBの対角化ルーチン DBDSQRはBTBの固有値λjとともに対応する固有ベクトルvjを求めることができるので、他のルーチンを併用しなくてもBTBの対角化（Bの特異値分解）が得られる。しかしながら、DLASQおよびdLVsは上述したように固有ベクトルvjを計算する機能はない。よって、別の固有ベクトル計算ルーチンが必要となる。ただし、BTBの固有値λj（Bの特異値σj）が求まっているものとする。 【0042】 まず、DSTEGRによる固有ベクトルの計算方法を示す。Parlettらは以下の（1）、（2）、（3）の手順で固有ベクトル計算ができることを示している。 （1）stationary qd with shift変換（stqds変換）およびreverse－time progressive qd with shift変換（rpqds変換）を用いた（B（0））TB（0）－λjIのCholesky（コレスキー）分解（B（±1））TB（±1）＝（B（0））TB（0）－λjIを求める。ただし、B（0）＝Bであり、B（＋1）およびB（－1）はそれぞれ上2重対角および下2重対角行列とする。すなわち 【0043】 【数9】 となる。すなわち、｛qk（±1），ek（±1）｝を計算する。 （2）コレスキー分解より（B（0））TB（0）－λkIのツイスティッド分解 【0044】 【数10】 を形成する。ただし 【0045】 【数11】 とし、γk≡qk（1）＋qk（－1）－（ek－1（0）＋qk（0）－λk）が最小になるkをk0とする。 （3）簡約化された連立1次方程式 【0046】 【数12】 を解く。 【0047】 図2は、Parlettらの手法で最も核となる（1）の処理手順を示す。（1）で求めたいのは、｛qk（n），ek（n）｝，n＝±1であるが、そのための変換としてstationary qd with shift（stqds）変換 qk（1）＋ek－1（1）＝qk（0）＋ek－1（0）－λj，k＝1，2，・・・，m， qk（1）ek（1）＝qk（0）ek（0），k＝1，2，・・・，m－1， e0（0）≡0，e0（1）≡0 およびreverse－time progressive qd with shift（rpqds）変換 qk（－1）＋ek（－1）＝qk（0）＋ek－1（0）－λj，k＝1，2，・・・，m， qk＋1（－1）ek（－1）＝qk（0）ek（0），k＝1，2，・・・，m－1， e0（0）≡0，em（－1）≡0 が採用されている。 【0048】 固有値λjが既知ならば反復的な計算は不要なため計算量がQRアルゴリズムに比べて圧倒的に少なく抑えられるが、常に数値安定かつ精度のよい方法とはいえない。stqds変換とrpqds変換ともに減算による桁落ちが発生する可能性がある。例えば、stqds変換においてqk（0）＋ek－1（0）－λj～ek－1（1）ならば、qk（1）（＝qk（0）＋ek－1（0）－λj－ek－1（1））を求める際、倍精度計算であってもqk（1）の有効桁はわずか1桁になることもある。その場合、qk（0）ek（0）／qk（1）を計算すると誤差が生じる、つまり、ek（1）が精度よく計算できないことになる。また、qk＋1（1）を求めるのにek（1）が必要、ek（1）を求めるのにqk（1）が必要と逐次的な計算が要求されるので、1個所で発生した桁落ちによる誤差が波及し、さらなる誤差増大の可能性も秘めている。その結果、理論上はqk（1）≠0であるが誤差蓄積によりqk（1）＝0となり、qk（0）ek（0）／qk（1）の計算においてオーバーフローが起こるといった数値不安定な状況も想定される。 B（0）＝Bの成分｛b2k－1（0），b2k（0）｝が与えられる、すなわち｛qk（0），ek（0）｝が与えられると、λjおよびek－1（1）が一意的に決まるので、この状況は避けることはできない。rpqds変換も同様の性質を持つため、実用的なレベルにまで達したとはいいがたい。LAPACKにおいてFORTRANルーチンDSTEGRとして改良版が公開されているものの欠点は完全に解決されてはいない。 【0049】 次に、本発明による対角化のルーチンを説明する。本発明によるルーチンでは、stqds変換およびrpqds変換にかわる新たな変換を採用している。本発明では、行列BTBを対角化するために、sdLVvs変換およびrdLVvs変換を用いた行列BTB－σj2Iに対するツイスティッド分解を行う。正確には、上記（1）を以下の（1a）、（1b）、（1c）の3つの工程にわけることで数値安定なコレスキー分解を実現している。 【0050】 図3～図8は、本発明によるミウラ型逆変換、sdLVvs変換、rdLVvs変換およびミウラ型変換を実行する処理手順を示す。 （1a）ミウラ型逆変換 【0051】 【数13】 を行う。 【0052】 u2k－1（0）＝η（0）tk（0），k＝1，2，・・・，m， u2k（0）＝ek（0）／tk（0），k＝1，2，・・・，m－1， tk（0）≡qk（0）／（η（0）＋u2k－2（0））－1，η（0）≡1／δ（0） ただし、δ（0）は任意に選ぶことができる。 （1b）stationary不等間隔差分Lotka－Volterra（sdLVvs）変換 uk（1）＝uk（0）＊（1＋δ（0）uk－1（0））／（1＋δ（1）uk－1（1））， k＝1，2，・・・，2m－1， によって 【0053】 【数14】 reverse－time不等間隔差分Lotka－Volterra（rdLVvs）変換 uk（－1）＝uk（0）＊（1＋δ（0）uk－1（0））／（1＋δ（－1）u k＋1（－1））， k＝1，2，・・・，2m－1， によって 【0054】 【数15】 を実行する。ただし、uk（n）≡0，u2m（0）≡0とし、λj＝1／δ（0）－1／δ（n），n＝±1を満たすようにδ（n），n＝±1を設定する。 （1c）ミウラ型変換 【0055】 【数16】 qk（n）＝η（n）＊（1＋δ（n）u2k－2（n））＊（1＋δ（n）u2k－1（n））， k＝1，2，・・・，m， ek（n）＝δ（n）＊u2k－1（n）＊u2k（n） k＝1，2，・・・，m－1， を行う。 【0056】 qd型変換では見られない離散Lotka－Volterra型変換の大きな特徴は、任意パラメータを持つことである。すなわち、λj＝1／δ（0）－1／δ（±1）を満たす範囲でδ（n）の値を任意に設定できる。δ（n）を変動させると補助変数uk（n）の値も変化するが、桁落ちによる誤差や数値不安定が発生するかどうかは事前に判定できる。この判定は、if文によって実装されてもよい。この場合は、δ（n）を再設定後に再度計算すればよい。また、uk（±1）が求まればミウラ変換によってqk（±1）およびek（±1）が独立に計算されるので、誤差が伝播しないという性質を持つ。なお、ミウラ型逆変換をミウラ変換、ミウラ型変換を逆ミウラ変換、sdLVvs変換をstLVv変換、rdLVv変換をrLVv変換と呼び変えても問題はない。 【0057】 入力および出力の変数の対応は、Parlettの手法と同じである。メモリ消費を抑えるために、補助変数uk（n）のための配列は必ずしも用意する必要はない。一方、1＋δ（0）u（0）のためのメモリ領域を確保し、（1a）～（1c）のステップにまたがってこの値を利用することで、メモリ消費を抑え、計算量を低減することができる。これにより、誤差も低減される。このように、1．3で説明する方法を利用して誤差を低減すること等、さらなる工夫を行い得る。 【0058】 上記の方法でBTB－λjIのコレスキー分解を求めた後は、Parlettらと同じ方法を用いて特異ベクトルvjを計算する。そうすると、BTBの固有値分解Λ＝VBTBTBVBが得られる。ただしVBの列ベクトルがvjである。UB＝BVBΣ－1よりUBを求めると、Bの特異値分解B＝UBΣVBTも得られる。 【0059】 1．2．高精度特異値算出dLVsルーチン 本発明によるdLVsルーチンは、基本的な枠組みではDLASQルーチンと同じであり、BTBの固有値（あるいはBの特異値）のみが得られBTBの固有ベクトル（あるいはBの特異ベクトル）は求まらない。ただ、ルーチンの核となる部分に使用される計算法はDLASQと異なる。 【0060】 まず、DLASQルーチンによる固有値の計算方法を説明する。qd法は固有ベクトルが計算できないが、1ステップは上2重対角行列の要素を更新するのみなので高速に特異値を求めることができる。かつ、計算量が少ないため丸め誤差の蓄積が抑えられ高精度計算が期待できる。ここで、Y（n）を対称な3重対角行列 【0061】 【数17】 とする。また、変換Y（n）→Y（n＋1）がqd法の漸化式 qk（n＋1）＝qk（n）＋ek（n）－ek－1（n＋1），k＝1，2，・・・，m， ek（n＋1）＝ek（n）＊qk＋1（n）／qk（n＋1），k＝1，2，・・・，m－1， e0（n）≡0，em（n）≡0，n＝0，1，・・・ によって実現されるとする。そのとき、適切なR（n）が存在してY（n＋1）＝R（n）Y（n）（R（n））－1となることをRutishauserは発見している。これは、Y（n＋1）とY（n）の固有値が同じであることを意味し、すなわち、上記の漸化式によって固有値保存変形がなされる。この変形を繰り返せば非対角成分が0に近づき、Y（n）が対角化されることも証明されている。よってY（0）＝BTBとするとlim n→∞Y（n）＝diag（λ1，λ2，・・・，λm）となる。ただし、λkはBTBの固有値である。さらに、λkの正の平方根をとるとBの特異値σkが得られる。 【0062】 LAPCAKの特異値計算ルーチンDLASQでは、qd法のdifferential型が採用されている。これはdqd法と呼ばれる。1ステップは 【0063】 【数18】 で与えられ、変換）Y（n）→Y（n＋1）が実現される。differential型は減算を含まないため、qd法では起こりうる桁落ち誤差の心配もない。さらに高速版として原点シフト付きdqd（dqds）法も組み込まれている。高速化に伴う計算量の減少で丸め誤差の発生も抑えることができる。ただし、収束するかどうかの保証はない。原点シフト版の1ステップは 【0064】 【数19】 で与えられる。sは原点シフト量であり、Gersgorin型境界条件により見積もられたY（n）の最小固有値とする。よって、DLASQルーチンのメインループの核となる部分では最小固有値の見積もりを行い、s＞ε（ε：非常に小さいとみなせる正の数）ならばdqds法が、そうでなければdqd法が用いられる。メインループの他の部分（行列の分割・減次）はほとんどDBDSQRルーチンと変わらない。すなわち、DBDSQRルーチンに含まれるQR法をdqd法に、原点シフトQR法をdqds法に置き換えたものがDLASQルーチンとなる。（QR法とdqd法、原点シフト付きQR法とdqds法で扱う変数が異なるので少しだけ判定式が異なる部分もあるが本質的に同じである。） 次に本発明によるdLVsルーチンを説明する。dLVsルーチンではdqd法のかわりにdLV法を、dqds法のかわりにdLVs法を採用する。ここで、dLVsルーチンはDLASQルーチンとは違って収束することが保証されている。その1ステップは以下で与えられる。 【0065】 【数20】 DLASQルーチン同様、メインループの核の部分では最小特異値の見積もりによって原点シフト量sを決定し、dLV法を使うかdLVs法を使うかを決定する。ただし、補助変数vk（n）を求めたあとでsを決定する点はDLASQルーチンと異なる。 【0066】 1．3．本発明の方法による誤差の低減 人間が理想的な状況、すなわち、無限桁の計算をいくらでもできるとすると、Parlettらの手法、本発明による特異ベクトル計算ルーチン（1．1参照）のどちらを使っても問題ない。しかしながら、コンピュータで計算を行う場合は注意が必要である。有限桁の計算しか行えないコンピュータ上では、数学的に正しい計算法を使ったとしても必ずしも正しい結果が得られる訳ではない。そればかりかいつまでたっても計算が終了しないといった思わぬ数値的な問題が発生する場合もある。 【0067】 コンピュータ計算による誤差には、丸め誤差および桁落ちによる誤差などが知られている。丸め誤差単独では、高々有効桁の最後の桁が真値と比べて異なる程度で大きな誤差にはならない。また、指数部が異なる2つの実数の加算、乗算、除算の少なくとも1つの演算を行えばやはり丸め誤差が生じるが、それ以上の誤差は発生しない。さらに、このような丸め誤差を発生するような操作が繰り返されても、丸めモードがnear（四捨五入）ならば、一方的に切り上げられたりあるいは切り捨てられたりして誤差が極端に蓄積することは少ない。よって、多くの数値計算法は加算、乗算、除算の少なくとも1つの演算によって発生する丸め誤差を特別注意することは少なく、dLVsルーチンによる特異値計算でも結果的に丸め誤差は一様に増大しない。 【0068】 問題となるのは、同符号の実数の減算および異符号の実数の加算により生じる、桁落ちである。桁落ちによる誤差で値が0となった後、その値による除算を行うと、0が分母にくるような不定形となり計算不可能となる。こうなるといつまでたっても計算が終了しない。減算→除算と計算する部分がParlettの手法、本発明による方法ともに存在するので、桁落ち誤差には十分に注意する必要がある。 【0069】 本発明による計算方法では、上述の桁落ちによる誤差を含んでいるかどうかは減算によって得られた値が小さいかどうかで判断できる。Parlettの手法の場合、図2のDO文中のq1［k］、q2［k］をチェックすればよい。ところが、Parlettの手法は桁落ち誤差を含むことが分かったとしても、それを回避することはできない。なぜならば、初期値として｛q0[k]，e0[k]｝が与えられると、λは一意的に決定され、｛q1[k]，e1[k]｝および｛q2[k]，e2[k]｝も一意的に導出されるためである。すなわち、任意パラメータを持たない自由度のない計算法であるためである。 【0070】 それに対して本発明による方法は自由に設定できるパラメータη1（＝1／δ（0））をもつため、補助変数uk（n）の値を様々に変化させることができる（図9参照）。すなわち、様々な経路で｛q1[k]，e1[k]｝および｛q2[k]，e2[k]｝を計算することができる。よって、桁落ちが発生する場合も回避できる。図6～図8の条件判定によって桁落ちの影響をチェックし、減算によって得られた値の絶対値が小さな数εより大きいという条件が満たされなければ、パラメータη1の設定に戻るというものである。この処理は、条件が満たされるまで繰り返される。なお、精度よりも高速性を重視する場合は、数回条件が満たさなければ（q1[k]＝0あるいはq2[k]＝0ならば）、例外処理（文献PDF参照）を行ってもよい。 【0071】 2．本発明による特異値分解装置 特異値分解を実行する装置を説明する。特異値分解を実行する装置は、例えば、図1に示される各ステップの機能を有する装置である。特異値分解を実行する装置の1つの実施形態は、特異値分解プログラムを実行するコンピュータである。そのプログラムは、例えば、図1に示される各ステップをCPUに実行させるプログラムである。 【0072】 図10は、図1の処理をコンピュータプログラムとして実行するコンピュータ10を示す。コンピュータ10は、CPU1001と、メモリ1002と、入力インターフェース部1003と、出力インターフェース部1004と、バス1005とを含む。 【0073】 CPU1001は、プログラムを実行する。例えば、そのプログラムは、図1の処理を実行するプログラムである。そのプログラムおよびそのプログラムの実行に必要なデータは、例えば、メモリ1002に格納される。そのプログラムは、任意の様態でメモリ1002に含まれ得る。例えば、メモリ1002が書き換え可能なメモリである場合には、コンピュータ1002の外部からプログラムをローディングすることにより、そのプログラムをメモリ1002に格納してもよい。あるいは、メモリ1002が読み出し専用メモリである場合には、メモリ1002に焼き付ける形式でそのプログラムをメモリ1002に格納してもよい。 【0074】 さらに、入力インターフェース部1003は、特異値分解を実行する対象である行列Aを外部から受け取るためのインターフェースとして機能する。出力インターフェース部1004は、計算結果を出力するためのインターフェースとして機能する。バス1005は、コンピュータ10内の構成要素1001～1004を相互に接続するために使用される。 【0075】 3．本発明による特異値分解方法の応用 本発明による特異値分解方法は、様々な分野に応用可能である。以下に、2次元画像から3次元画像へ復元する画像処理への応用例、ならびに、文書検索への応用例を示す。しかし、これら2つの応用例は例示に過ぎず、本発明による特異値分解方法の応用は、これら2つの応用例に限定されない。本発明による特異値分解方法の応用は、特異値分解を利用するあらゆる分野に適用可能である。 【0076】 3．1．2次元から3次元へ復元する画像処理への応用 図11を参照する。図11は、本発明による特異値分解方法を利用した物体の複数の2次元画像を3次元画像へ復元する画像処理の1つの実施形態を説明する。 【0077】 複数の2次元画像から3次元復元を行うために必要となるステップは、2次元画像から特徴点を抽出するステップと、特徴点データより形状（元の物体の特徴点の3次元座標データ）および回転（3次元データから特徴点データへの変換）に関するデータを計算するステップと、形状および回転のデータより可視化を行うステップとを包含する。 【0078】 ステップ1102において、2次元画像j（j＝1，・・・，m、mは3以上の整数）から特徴点i（i＝1，・・・，n、nは2以上の整数）の座標（xij，yij）を抽出する。取り扱う2次元画像は、弱中心射影画像であることが好ましい。このとき 【0079】 【数21】 が成り立つ。ここで、sjは物体のスケールに相対するj番目の画像のスケール、r1，j，r2，jはそれぞれ物体座標系に相対するj番目のカメラ座標系の回転行列の1番目と2番目の行ベクトル、［Xi，Yi，Zi］Tはi番目の点の3次元座標である。物体のスケールは1番目の画像のスケールと同じにし（s1＝1）、物体の座標系の姿勢は1番目の画像のカメラ座標系と同じにする（r1，1＝［1，0，0］T，r2，1＝［0，1，0］T）。m枚全ての画像におけるn個全ての点の座標を行列Dの要素として並べると、 D＝MS が得られる。 ただし、 【0080】 【数22】 である。MとSの形から分かるように、Dのランクは3である。いま、ステップ1102よりDが与えられている。以下、回転に関するデータMおよび形状Sを求める。 【0081】 そこで、行列Dの特異値分解 D＝UΣVT を考える。ここで、Σは特異値を大小順に対角線上に並べたもので、UおよびVはそれぞれ左および右直交行列である。ステップ1104では、行列Dの特異値分解を計算するために、行列Dを上2重対角化して上2重対角化行列Bを得る。ステップ1106において、行列B、すなわちDの特異値を計算する。ステップ1106では、本発明によるdLVsルーチン、DLASQルーチン、その他の特異値計算ルーチン、あるいはその組み合わせを利用し得る。このとき、画像のデジタル誤差のため、ゼロでない特異値は3つ以上出てくる。しかし、4番目以降の特異値はノイズによるもので、最初の3つの特異値と比べて格段に小さい。 【0082】 そこで、ステップ1108において、最初の3つの特異値に対して特異ベクトルを計算する。ステップ1108は、図1のステップ110～120を利用し得る。採用する3個のベクトルをまとめると、次式となる。 【0083】 D’＝L’Σ’R’T＝M’S’ ただし、M’＝L’（Σ’）1／2、S’＝（Σ’）1／2R’T、D’は∥D－D’∥を最小にするランク3の行列である。 【0084】 次に、D’からMおよびSを求めたいが、その組合せは唯一ではない。なぜなら、任意の正則行列Cが D’＝（M’C）（C－1S’） を満たすからである。そこで、上式におけるCをM＝M’Cを満たすように決める。Cは下記の式を満たす。 【0085】 【数23】 E＝CCTとすると、上式からEの6つの要素に関する2m＋1個の線形方程式が得られる。m≧3であるので、Eの要素を一意に決めることができる。ステップ1110において、行列Eを求める。 【0086】 次に、ステップ1112において、EからCを求める。Cの自由度（9）はEの自由度（6）より多い。そこで、条件r1j＝［1，0，0］T，r2j＝［0，1，0］Tを加えれば、Cを決めることができる。このとき2つの解（Necker Reversal）が出る。 【0087】 次に、ステップ1114において、M＝M’CおよびS＝C－1S’より、回転に関するデータMおよび形状Sが決まる。 【0088】 DBDSQRルーチンでは上記のようにU、V（n個のベクトル）を求めた後、3個のベクトル（n－3個は不要）を採用する。ここでは、原則的にはn個のベクトルを求めるが、例外的にn個より少ない数のベクトルを求めてもよい。この実施形態では3個のベクトルだけ求めればよい。すなわち、計算コストを削減することができる。 【0089】 3．2．本発明による特異値分解方法を利用した文書検索方法 文書中からその文書の内容に関連する索引語を抽出し、索引語の重みを計算する処理の後、ベクトル空間モデルでは、この索引語の重みを要素とするベクトルで文書を表現する。いま、検索対象となる文書をD1，D2，・・・，Dnとし、これら文書集合全体を通して全部でm個の索引語w1，w2，・・・，wmがあるとする。このとき、文書Djは、次のようなベクトルで表現されることになる。これを文書ベクトルと呼ぶ。 【0090】 【数24】 ここで、dijは索引語wiの文書Djにおける重みである。また、文書集合全体は、次のようなm×n行列Dによって表現することができる。 【0091】 【数25】 行列Dを索引語文書行列と呼ぶ。索引語文書行列の各列は文書に関する情報を表している文書ベクトルであるが、同様に、索引語文書行列の各行は索引語に関する情報を表しているベクトルであり、これを索引語ベクトルと呼ぶ。検索質問も、文書と同様に、索引語の重みを要素とするベクトルで表現することができる。検索質問文に含まれる索引語wiの重みをqiとすると、検索質問ベクトルqは次のように表されることになる。 【0092】 【数26】 実際の文書検索においては、与えられた検索質問文と類似した文書を見つけ出す必要があるが、検索質問ベクトルqと各文書ベクトルdjの間の類似度を計算することにより行う。ベクトル間の類似度の定義としてはさまざまなものが考えられるが、文書検索においてよく用いられているものはコサイン尺度（2つのベクトルのなす角度）または内積である。 【0093】 【数27】 なお、ベクトルの長さが1に正規化（コサイン正規化）されている場合には、コサイン尺度と内積とは一致する。 【0094】 図12を参照する。図12は、本発明による特異値分解装置を利用した、文書検索方法の1つの実施形態を説明する。 【0095】 ステップ1202において、質問ベクトルqを受け取る。次に、ステップ1204では、行列Dの特異値分解を計算するために、行列Dを上2重対角化して上2重対角化行列Bを得る。ステップ1206では、行列B、すなわち行列Dの特異値を求める。ステップ1206では、本発明によるdLVsルーチン、DLASQルーチン、その他の特異値計算ルーチン、あるいはその組み合わせを利用し得る。 【0096】 次に、Dの近似行列Dkを使った検索を考える。ベクトル空間モデルでは、検索質問ベクトルqと索引語文書行列D中の各文書ベクトルdjの間の類似度を計算することにより検索を行うが、ここではDの代わりにDkを使う。ベクトル空間モデルでは、文書ベクトルの次元数は索引語の総数と等しい。したがって、検索対象となる文書の数が増えるに従い、文書ベクトルの次元数も増加する傾向にある。しかし、次元数が増加してくると、コンピュータのメモリによる制限や検索時間の増大などの問題が生じてくるばかりでなく、文書中に含まれる不必要な索引語がノイズ的な影響を及ぼし、検索精度を低下させてしまうという現象も起こってくる。潜在的意味インデキシング（latent semantic indexing；LSI）は、高次元の空間にある文書ベクトルを低次元の空間へと射影することにより、検索精度の改善を図る技術である。高次元の空間では別々に扱われていた索引語が、低次元の空間では相互に関連を持ったものとして扱われる可能性もあるため、索引語の持つ意味や概念に基づく検索を行うことができる。たとえば、通常のベクトル空間モデルでは“car”という索引語と“automobile”という索引語はまったく別物であり、一方の索引語による質問ではもう片方の索引語を含んだ文書を検索することができない。しかし、低次元の空間ではこれらの意味的に関連した索引語は1つの次元に縮退することが期待できるため、“car”という検索質問によって“car”を含む文書ばかりでなく“automobile”を含む文書をも検索することが可能となる。潜在的意味インデキシングでは、特異債分解により高次元ベクトルの次元削減を行うが、これは基本的に多変量解析における主成分分析と等価である。 【0097】 行列ステップ1208では、k＜rなるkの値を選択する。kの値は、予め与えられていてもよいし、計算ごとに選択可能であってもよい。 【0098】 ステップ1210において、ステップ1206で計算した特異値のうち、大きい順に1番目からk番目のk個の特異値に対して、Dの特異ベクトルを計算する。すなわち、 Dk＝UkΣVkT なるUkおよびVkを計算する。ここで、Ukは、最初のk個の左特異ベクトルのみから構成されるm×k行列であり、Vkは、最初のk個の右特異ベクトルのみから構成されるn×k行列であり、Σは、最初のk個の特異値のみから構成されるk×k対角行列である。ステップ1210は、図1のステップ110～120を利用し得る。 【0099】 次に、ステップ1212において、行列Dkと質問ベクトルqとの類似度を計算する。いま、ベクトルejをn次元の単位ベクトルとすると、Dkのj番目の文書ベクトルはDkejで表すことができる。文書ベクトルDkejと検索質問ベクトルqとの間の類似度計算は、 【0100】 【数28】 としてもよいが、別の定義を用いてもよい。上式では、DkをUk，Σk，Vkから再構成する必要はなく特異値分解の結果から、直接、類似度を計算できることを示している。 数29の中に現われるΣkVkTejは、 ΣkVkTej＝UkTDkej と書き直すことができる。この式の右辺は、近似行列Dkにおけるj番目の文書ベクトルの基底Ukのもとでの座標（文書のk次元表現）を表している。同様に、数29中のUkTqは、検索質問ベクトルqの基底Ukのもとでの座標（検索質問のk次元表現）である。 【0101】 ステップ1214において、ステップ1212において計算された類似度を基準に、検索結果を出力する。 【0102】 DBDSQRルーチンではr個のベクトルを求めた後、k個のベクトル（r－k個は不要）を採用するのに対して、この実施形態では、k個のベクトルだけ求めればよい。すなわち、計算コストを削減することができる。 【0103】 以上のように、本発明の好ましい実施形態を用いて本発明を例示してきたが、本発明は、この実施形態に限定して解釈されるべきものではない。本発明は、特許請求の範囲によってのみその範囲が解釈されるべきであることが理解される。当業者は、本発明の具体的な好ましい実施形態の記載から、本発明の記載および技術常識に基づいて等価な範囲を実施することができることが理解される。本明細書において引用した特許、特許出願および文献は、その内容自体が具体的に本明細書に記載されているのと同様にその内容が本明細書に対する参考として援用されるべきであることが理解される。 【産業上の利用可能性】 【0104】 本発明は、任意の行列Aの特異値分解を高速かつ高精度に実行することを可能にする方法、プログラム、および装置等として有用である。
Drawing	(In Japanese)【図1】 (In Japanese)【図2】 (In Japanese)【図3】 (In Japanese)【図4】 (In Japanese)【図5】 (In Japanese)【図6】 (In Japanese)【図7】 (In Japanese)【図8】 (In Japanese)【図9】 (In Japanese)【図10】 (In Japanese)【図11】 (In Japanese)【図12】