JP4325877B2

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

Info

Publication number: JP4325877B2
Application number: JP2006514122A
Authority: JP
Inventors: 佳正中村; 雅史岩▲崎▼; 真也阪野
Original assignee: Japan Science and Technology Agency; National Institute of Japan Science and Technology Agency
Current assignee: Japan Science and Technology Agency; National Institute of Japan Science and Technology Agency
Priority date: 2004-06-03
Filing date: 2005-06-01
Publication date: 2009-09-02
Anticipated expiration: 2025-06-01
Also published as: EP1752884A4; US20090028455A1; EP1752884B1; WO2005119507A1; CA2568852A1; CA2568852C; JPWO2005119507A1; US8306361B2; EP1752884A1

Description

本発明は、任意の行列Ａに対して高速かつ高精度に特異値分解を実行する方法、プログラムおよび装置に関し、より詳細には、行列の特異値を計算し、計算された特異値に対する特異ベクトルを計算することによって、高速かつ高精度に特異値分解を実行する方法、プログラムおよび装置に関する。

現在、コンピュータ上で特異値分解を行うために用いられる標準的な方法は、国際的な数値計算ライブラリＬＡＰＡＣＫにより公開されているＤＢＤＳＱＲルーチンである。ＤＢＤＳＱＲルーチンは、ＱＲ法に基づいて作られたものであり、前処理、メインループ、後処理の３つの部分に分けることができる。前処理では、特異値の上限および下限を計算し、収束判定に用いる精度を計算する。メインループでは、ＱＲ法を繰り返しながら、除々に行列の分割および減次を行い、最終的に行列サイズが１となった時点でループを終了する。途中で２×２行列のブロックが現れたならば別処理を行う。後処理では、特異値として計算された値が負ならば正に直し、必要ならば特異ベクトルをスカラー倍する。特異値を降順に、それに対応するように特異ベクトルも並べ替える。ＤＢＤＳＱＲは、計算量が非常に多く、特に大規模問題では計算時間の増大が避けられない。ＤＢＤＳＱＲルーチンでは、特異値および特異ベクトルが同時に計算される。また、ＬＡＰＡＣＫには、特異値を計算するＤＬＡＳＱルーチン、および、計算された特異値を用いて特異ベクトルを計算する際に用いる対角化のためのＤＳＴＥＧＲルーチンがある。ＤＬＡＳＱルーチンは、高速高精度に特異値を求めることができるが、特異ベクトルを計算できない。ＤＳＴＥＧＲルーチンは、その数値的な欠点から、実際に特異値分解に使用することは難しい。

ＬＡＰＡＣＫによるＤＢＤＳＱＲルーチンを例として、従来の特異値分解方法を説明する。ＤＢＤＳＱＲルーチンでは、ｌ_１×ｌ_２の一般行列Ａを特異値分解するために、まず、行列Ａに対してＨｏｕｓｅｈｏｌｄｅｒ（ハウスホルダー）変換を適用する。すなわち、行列Ａは、直交行列Ｕ_Ａ、Ｖ_Ａを用いて、

と表わすことができる。このとき得られるＢを、上２重対角行列という。ここで、Ａの特異値＝Ｂの特異値であることに注意する。このように、一般行列Ａに対する特異値分解問題を上２重対角行列Ｂに対する特異値分解問題
Ｂ＝Ｕ_ＢΣＶ_Ｂ ^Ｔ
行列Ｕ_Ｂ、Ｖ_Ｂは、それぞれ左および右直交行列
Σ≡ｄｉａｇ（σ_１，σ_２，・・・，σ_ｍ） σ_１≧σ_２≧・・・≧σ_ｍ≧０
σ_ｊはＢの特異値
に置き換える。

次に、行列Ｂ^ＴＢを考える。この行列Ｂ^ＴＢの対角化
Λ＝Ｖ^ＴＢ^ＴＢＶ
Λ≡ｄｉａｇ（λ_１，λ_２，・・・，λ_ｍ） λ_１≧λ_２≧・・・≧λ_ｍ≧０
Ｖ≡（ｖ_１，ｖ_２，・・・，ｖ_ｍ）
λ_ｊはＢ^ＴＢの固有値
ｖ_ｊは固有値λ_ｊに対する固有ベクトル
を行う。ここで、一般に、以下のことが成り立つ。（１）Ｂ^ＴＢは対称な３重対角行列である。（２）Ｂ^ＴＢの固有値は全て正であり、Ｂの特異値σ_ｊ（σ_１≧σ_２≧・・・≧σ_ｍ≧０）は、Ｂ^ＴＢの固有値λ_ｊ（λ_１≧λ_２≧・・・≧λ_ｍ≧０）の正の平方根に等しい。（３）Ｖ_Ｂ＝Ｖ、すなわち、Ｂ^ＴＢの固有ベクトルｖ_ｊは、Ｂの右特異ベクトルに等しい。従って、行列Ｂ^ＴＢの対角化が求まると、（２）よりΛ＝Σ^２であることから、Σが求まり、さらに（３）より左直交行列Ｕ_Ｂ＝ＢＶ_ＢΣ^−１＝ＢＶΣ^−１が求まるので、Ｂが特異値分解される。すなわち、Ｂの特異値分解を、Ｂ^ＴＢの対角化の問題に置き換えることができる。この原理は、ｍ個全ての特異値および特異ベクトルを求める場合だけでなく、少なくとも１つの特異値および特異ベクトルを求める場合にも適用され得る。

以上のように、一般行列Ａの特異値分解は、対称な３重対角行列Ｂ^ＴＢの対角化の問題を含む。

ＤＢＤＳＱＲルーチンでは、計算量が非常に多いため、特に大規模問題では、高速に特異値分解を実行することが困難であった。一方、ＤＬＡＳＱルーチンは、高速かつ高精度に特異値を求めることができるが、ＤＳＴＥＧＲルーチンは、場合によっては特異ベクトルを低精度で計算するため、ＤＳＴＥＧＲルーチンでは、常に高精度に特異値分解を実行することは困難であった。

本発明の目的の１つは、対称な３重対角行列Ｂ^ＴＢの対角化を高速かつ高精度で実行することにより、任意の行列Ａの特異値分解を高速かつ高精度に実行することを可能にする方法、プログラム、および装置を提供することにある。ここで、行列Ｂを、行列Ａをハウスホルダー変換することで得られる上２重対角行列とする。

本発明では、ミウラ型逆変換、ｓｄＬＶｖｓ変換、ｒｄＬＶｖｓ変換およびミウラ型変換によるＴｗｉｓｔｅｄ（ツイスティッド）分解を実行することにより、行列Ｂ^ＴＢを対角化する。ここで、行列Ｂを、行列Ａと特異値が同じ上２重対角行列とする。行列Ａから行列Ｂへの変換は、例えばハウスホルダー法で実行できる。本発明による行列Ｂ^ＴＢの対角化では、固有値および固有ベクトルを同時に求めるＤＢＤＳＱＲルーチンとは異なり、ＤＳＴＥＧＲルーチンのように、まず固有値を計算し、次に計算された固有値を用いて固有ベクトルを求める。

本発明により、コンピュータを用いて任意の行列Ａに対して特異値分解を実行する方法であって、行列Ａに対して上２重対角化を行い、行列Ａの上２重対角行列Ｂを求めるステップと、行列Ａの特異値として行列Ｂの少なくとも１つの特異値σを求めるステップと、σに対する行列Ａの特異ベクトルを求めるステップとを包含し、該行列Ａの特異ベクトルを求めるステップは、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列Ｂ^ＴＢ−σ^２Ｉに対してツイスティッド分解を実行することにより、行列Ｂ^ＴＢを対角化するステップを含み、Ｉは単位行列である、方法が提供され、これにより上記目的が達成される。

前記行列Ａの特異ベクトルを求めるステップは、前記行列Ｂ^ＴＢを対角化するステップの後に、逆反復法を実行するステップを含んでもよい。

前記行列Ａの特異ベクトルを求めるステップは、前記逆反復法を実行するステップの後に、グラムシュミット法を実行するステップを含んでもよい。

前記行列Ｂの特異値σを求めるステップは、ｄＬＶｓルーチンを実行するステップを含んでもよい。

前記行列Ｂの特異値σを求めるステップは、要求される計算時間および精度に応じて、ｄＬＶｓルーチンもしくはＤＬＡＳＱルーチンのどちらを実行するのかを判定するステップを含んでもよい。

前記行列Ａに対して上２重対角化を行い、上２重対角行列Ｂを求めるステップは、上２重対角化法（たとえばハウスホルダー法）を実行するステップを含んでもよい。

本発明により、物体の複数の２次元画像から３次元画像を復元する方法であって、２次元画像ｊ（ｊ＝１，・・・，ｍ、ｍは３以上の整数）における特徴点ｉ（ｉ＝１，・・・，ｎ、ｎは２以上の整数）の座標ｄ_ｉｊ（ｘ_ｉｊ，ｙ_ｉｊ）を抽出するステップと、該特徴点の２次元座標ｄ_ｉｊ（ｘ_ｉｊ，ｙ_ｉｊ）より、該特徴点の３次元座標ｓ_ｉ（Ｘ_ｉ，Ｙ_ｉ，Ｚ_ｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算するステップとを包含し、該特徴点の３次元座標ｓ_ｉ（Ｘ_ｉ，Ｙ_ｉ，Ｚ_ｉ）および２次元座標から３次元
座標への変換を表す行列Ｍを計算するステップは、行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求めるステップであって、行列Ｄは、

と定義される、ステップと、行列Ｄの特異値として行列Ｂの特異値σ_ｊ（σ_１≧σ_２≧・・・≧σ_ｒ＞０、ｒは、行列Ｄのランクに等しい）を求めるステップと、σ_１、σ_２およびσ_３に対する行列Ｄの特異ベクトルを求めるステップと、Ｍ＝Ｍ’Ｃなる行列Ｃに対して、Ｅ＝ＣＣ^Ｔを満たす行列Ｅを計算するステップであって、Ｍ’＝Ｌ’（Σ’）^１／２、Σ’はσ_１、σ_２およびσ_３を対角要素に持ちそれ以外の要素が０である３×３行列、Ｌ’はσ_１、σ_２およびσ_３に対応するＤの特異ベクトルを左から順に並べた行列である、ステップと、行列Ｅから、行列Ｃを計算するステップと、行列Ｃから、該３次元座標ｓ_ｉ（Ｘ_ｉ，Ｙ_ｉ，Ｚ_ｉ）および該変換を表す行列Ｍを計算するステップとを含み、該σ_１、σ_２およびσ_３に対する行列Ｄの特異ベクトルを求めるステップは、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列Ｂ^ＴＢ−σ_ｊ ^２Ｉに対してツイスティッド分解を実行することにより、行列Ｂ^ＴＢを対角化するステップを含み、Ｉは単位行列である、方法が提供され、これにより上記目的が達成される。

本発明により、与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する方法であって、該キーワードに対応する質問ベクトルｑを受け取るステップと、該文書に対応する索引語文書行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求めるステップと、行列Ｄの特異値として行列Ｂの特異値σ_ｊ（σ_１≧σ_２≧・・・≧σ_ｒ＞０、ｒは、行列Ｄのランクに等しい）を求めるステップと、ｋ＜ｒなるｋを選択するステップと、行列Ｄの擬似行列Ｄ_ｋを計算するステップであって、行列Ｄ_ｋは、σ_１，σ_２，・・・，σ_ｋを対角要素としそれ以外の要素が０である行列Σ_ｋ、σ_１，σ_２，・・・，σ_ｋに対応する特異ベクトルのみを左から順に並べた左右直交行列Ｕ_ｋ，Ｖ_ｋを用いて、Ｄ_ｋ＝Ｕ_ｋΣ_ｋＶ_ｋ ^Ｔと定義される、ステップと、行列Ｄ_ｋと質問ベクトルｑとの類似度を計算するステップと、該計算された類似度を基準に、検索結果を出力するステップとを包含し、該行列Ｄ_ｋの左右直交行列Ｕ_ｋ，Ｖ_ｋを求めるステップは、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列Ｂ^ＴＢ−σ_ｊ ^２Ｉ（ｊ＝１，２，・・・，ｋ）に対してツイスティッド分解を実行することにより、行列Ｂ^ＴＢを対角化するステップを含み、Ｉは単位行列である、方法が提供され、これにより上記目的が達成される。

本発明により、任意の行列Ａに対して特異値分解を実行する特異値分解処理をコンピュータに実行させるプログラムであって、該特異値分解処理は、行列Ａに対して上２重対角化を行い、行列Ａの上２重対角行列Ｂを求めるステップと、行列Ａの特異値として行列Ｂの少なくとも１つの特異値σを求めるステップと、σに対する行列Ａの特異ベクトルを求めるステップとを包含し、該行列Ａの特異ベクトルを求めるステップは、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列Ｂ^ＴＢ−σ^２Ｉに対してツイスティッド分解を実行することにより、行列Ｂ^ＴＢを対角化するステッ
プを含み、Ｉは単位行列である、プログラムが提供され、これにより上記目的が達成される。

本発明により、物体の複数の２次元画像から３次元画像を復元する画像復元処理をコンピュータに実行させるプログラムであって、該画像復元処理は、２次元画像ｊ（ｊ＝１，・・・，ｍ、ｍは３以上の整数）における特徴点ｉ（ｉ＝１，・・・，ｎ、ｎは２以上の整数）の座標ｄ_ｉｊ（ｘ_ｉｊ，ｙ_ｉｊ）を抽出するステップと、該特徴点の２次元座標ｄ_ｉｊ（ｘ_ｉｊ，ｙ_ｉｊ）より、該特徴点の３次元座標ｓ_ｉ（Ｘ_ｉ，Ｙ_ｉ，Ｚ_ｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算するステップとを包含し、該特徴点
の３次元座標ｓ_ｉ（Ｘ_ｉ，Ｙ_ｉ，Ｚ_ｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算するステップは、行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求めるステップであって、行列Ｄは、

と定義される、ステップと、行列Ｄの特異値として行列Ｂの特異値σ_ｊ（σ_１≧σ_２≧・・・≧σ_ｒ＞０、ｒは、行列Ｄのランクに等しい）を求めるステップと、σ_１、σ_２およびσ_３に対する行列Ｄの特異ベクトルを求めるステップと、Ｍ＝Ｍ’Ｃなる行列Ｃに対して、Ｅ＝ＣＣ^Ｔを満たす行列Ｅを計算するステップであって、Ｍ’＝Ｌ’（Σ’）^１／２、Σ’はσ_１、σ_２およびσ_３を対角要素に持ちそれ以外の要素が０である３×３行列、Ｌ’はσ_１、σ_２およびσ_３に対応するＤの特異ベクトルを左から順に並べた行列である、ステップと、行列Ｅから、行列Ｃを計算するステップと、行列Ｃから、該３次元座標ｓ_ｉ（Ｘ_ｉ，Ｙ_ｉ，Ｚ_ｉ）および該変換を表す行列Ｍを計算するステップとを含み、該σ_１、σ_２およびσ_３に対する行列Ｄの特異ベクトルを求めるステップは、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列Ｂ^ＴＢ−σ_ｊ ^２Ｉに対してツイスティッド分解を実行することにより、行列Ｂ^ＴＢを対角化するステップを含み、Ｉは単位行列である、プログラムが提供され、これにより上記目的が達成される。

本発明により、与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する文書検索処理をコンピュータに実行させるプログラムであって、該文書検索処理は、該キーワードに対応する質問ベクトルｑを受け取るステップと、該文書に対応する索引語文書行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求めるステップと、行列Ｄの特異値として行列Ｂの特異値σ_ｊ（σ_１≧σ_２≧・・・≧σ_ｒ＞０、ｒは、行列Ｄのランクに等しい）を求めるステップと、ｋ＜ｒなるｋを選択するステップと、行列Ｄの擬似行列Ｄ_ｋを計算するステップであって、行列Ｄ_ｋは、σ_１，σ_２
，・・・，σ_ｋを対角要素としそれ以外の要素が０である行列Σ_ｋ、σ_１，σ_２，・・・，σ_ｋに対応する特異ベクトルのみを左から順に並べた左右直交行列Ｕ_ｋ，Ｖ_ｋを用いて、Ｄ_ｋ＝Ｕ_ｋΣ_ｋＶ_ｋ ^Ｔと定義される、ステップと、行列Ｄ_ｋと質問ベクトルｑとの類似度を計算するステップと、該計算された類似度を基準に、検索結果を出力するステップとを包含し、該行列Ｄ_ｋの左右直交行列Ｕ_ｋ，Ｖ_ｋを求めるステップは、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列Ｂ^ＴＢ−σ_ｊ
^２Ｉ（ｊ＝１，２，・・・，ｋ）に対してツイスティッド分解を実行することにより、行列Ｂ^ＴＢを対角化するステップを含み、Ｉは単位行列である、プログラムが提供され、これにより上記目的が達成される。

本発明により、任意の行列Ａに対して特異値分解を実行する装置であって、行列Ａを入力として受け取る手段と、行列Ａに対して上２重対角化を行い、行列Ａの上２重対角行列Ｂを計算する手段と、行列Ａの特異値として行列Ｂの少なくとも１つの特異値σを計算する手段と、σに対する行列Ａの特異ベクトルを求める手段とを備え、該行列Ａの特異ベクトルを求める手段は、ミウラ型逆変換、ｓｄＬＶｖｓ変換およびｒｄＬＶｖｓ変換、ならびにミウラ型変換を用いて行列Ｂ^ＴＢ−σ^２Ｉに対してツイスティッド分解を実行することにより、行列Ｂ^ＴＢを対角化する手段を含み、Ｉは単位行列である、装置が提供され、これにより上記目的が達成される。

本発明により、物体の複数の２次元画像を３次元画像に復元する装置であって、ｍ枚（ｍは３以上の整数）の２次元画像を受け取る手段と、２次元画像ｊ（ｊ＝１，・・・，ｍ）における特徴点ｉ（ｉ＝１，・・・，ｎ、ｎは２以上の整数）の座標ｄ_ｉｊ（ｘ_ｉｊ，ｙ_ｉｊ）を抽出する手段と、該特徴点の２次元座標ｄ_ｉｊ（ｘ_ｉｊ，ｙ_ｉｊ）より、該特徴点の３次元座標ｓ_ｉ（Ｘ_ｉ，Ｙ_ｉ，Ｚ_ｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算する手段とを備え、該特徴点の３次元座標ｓ_ｉ（Ｘ_ｉ，Ｙ_ｉ，Ｚ_ｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算する手段は、行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求める手段であって、行列Ｄは、

と定義される、手段と、行列Ｄの特異値として行列Ｂの特異値σ_ｊ（σ_１≧σ_２≧・・・≧σ_ｒ＞０、ｒは、行列Ｄのランクに等しい）を求める手段と、σ_１、σ_２およびσ_３に対する行列Ｄの特異ベクトルを求める手段と、Ｍ＝Ｍ’Ｃなる行列Ｃに対して、Ｅ＝ＣＣ^Ｔを満たす行列Ｅを計算する手段であって、Ｍ’＝Ｌ’（Σ’）^１／２、Σ’はσ_１、σ_２およびσ_３を対角要素に持ちそれ以外の要素が０である３×３行列、Ｌ’はσ_１、σ_２およびσ_３に対応するＤの特異ベクトルを左から順に並べた行列である、手段と、行列Ｅから、行列Ｃを計算する手段と、行列Ｃから、該３次元座標ｓ_ｉ（Ｘ_ｉ，Ｙ_ｉ，Ｚ_ｉ）および該変換を表す行列Ｍを計算する手段とを含み、該σ_１、σ_２およびσ_３に対する行列Ｄの特異ベクトルを求める手段は、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列Ｂ^ＴＢ−σ_ｊ ^２Ｉに対してツイスティッド分解を実行することにより、行列Ｂ^ＴＢを対角化する手段を含み、Ｉは単位行列である、装置が提供され、これにより上記目的が達成される。

本発明により、与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する装置であって、該キーワードに対応する質問ベクトルｑを受け取る手段と、該文書に対応する索引語文書行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求める手段と、行列Ｄの特異値として行列Ｂの特異値σ_ｊ（σ_１≧σ_２≧・・・≧σ_ｒ＞０、ｒは、行列Ｄのランクに等しい）を求める手段と、ｋ＜ｒなるｋを選択する手段と、行列Ｄの擬似行列Ｄ_ｋを計算する手段であって、行列Ｄ_ｋは、σ_１，σ_２，・・・，σ_ｋを対角要素としそれ以外の要素が０である行列Σ_ｋ、σ_１，σ_２，・・・，σ_ｋに対応する特異ベクトルのみを左から順に並べた左右直交行列Ｕ_ｋ，Ｖ_ｋを用いて、Ｄ_ｋ＝Ｕ_ｋΣ_ｋＶ_ｋ ^Ｔと定義される、手段と、行列Ｄ_ｋと質問ベクトルｑとの類似度を計算する手段と、該計算された類似度を基準に、検索結果を出力する手段とを包含し、該行列Ｄ_ｋの左右直交行列Ｕ_ｋ，Ｖ_ｋを求める手段は、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列Ｂ^ＴＢ−σ_ｊ ^２Ｉ（ｊ＝１，２，・・・，ｋ）に対してツイスティッド分解を実行することにより、行列Ｂ^ＴＢを対角化する手段を含み、Ｉは単位行列である、装置が提供され、これにより上記目的が達成される。

本発明によれば、ミウラ型逆変換、ｓｄＬＶｖｓ変換、ｒｄＬＶｖｓ変換およびミウラ型変換によるツイスティッド分解を実行することにより、行列Ｂ^ＴＢが対角化される。それにより、ＤＳＴＥＧＲルーチンと比較した場合、精度の点ではるかに優れ、実用化に耐え得る固有値計算を実現し、高精度の特異値分解が可能になる。また、逆反復法、Ｇｒａｍ−Ｓｃｈｍｉｄｔ（グラムシュミット）法、本発明において提供するｄＬＶｓルーチンのうちの少なくとも１つを併用することにより、ＤＢＤＳＱＲと比較して精度の点では同等になることもあるが、速度の点では圧倒的に高速である特異値分解を実現することができる。さらに、ＤＢＤＳＱＲルーチンと比較すると、固有値および固有ベクトルを同時の求めるのか、または、計算された固有値を用いて固有ベクトルを求めるのかの計算フローの違いにより、以下で説明する２次元から３次元への画像復元問題および文書検索問題等においては、計算時間をさらに短縮することができる。

図１は、本発明により提供される特異値分解のためのＩ−ＳＶＤルーチンの処理の手順を示すフローチャート図２は、Ｐａｒｌｅｔｔらの手法によるｓｑｄｓ変換およびｒｐｑｄｓ変換の処理の手順を示すフローチャート図３は、本発明によるミウラ型逆変換、ｓｄＬＶｖｓ変換、ｒｄＬＶｖｓ変換およびミウラ型変換の処理の手順を示すフローチャート図４は、本発明によるミウラ型逆変換、ｓｄＬＶｖｓ変換、ｒｄＬＶｖｓ変換およびミウラ型変換の処理の手順を示すフローチャート図５は、本発明によるミウラ型逆変換、ｓｄＬＶｖｓ変換、ｒｄＬＶｖｓ変換およびミウラ型変換の処理の手順を示すフローチャート。図６は、本発明によるミウラ型逆変換、ｓｄＬＶｖｓ変換、ｒｄＬＶｖｓ変換およびミウラ型変換の処理の手順を示すフローチャート図７は、本発明によるミウラ型逆変換、ｓｄＬＶｖｓ変換、ｒｄＬＶｖｓ変換およびミウラ型変換の処理の手順を示すフローチャート図８は、本発明によるミウラ型逆変換、ｓｄＬＶｖｓ変換、ｒｄＬＶｖｓ変換およびミウラ型変換の処理の手順を示すフローチャート図９は、本発明において誤差の低減をもたらす、Ｐａｒｌｅｔｔの手法と本発明による方法との計算経路の違いを示す図図１０は、本発明による特異値分解装置の１つの実施形態であるコンピュータを示す図図１１は、本発明による特異値分解方法を利用した物体の複数の２次元画像を３次元画像へ復元する画像処理の１つの実施形態を示すフローチャート図１２は、本発明による特異値分解装置を利用した文書検索方法の１つの実施形態を示すフローチャート

符号の説明

１０コンピュータ
１００１ＣＰＵ
１００２メモリ
１００３入力インターフェース部
１００４出力インターフェース部
１００５バス

以下、図面を参照しながら、本発明の実施形態を説明する。

１．特異値分解アルゴリズムＩ−ＳＶＤルーチン
図１を参照する。図１は、本発明により提供される特異値分解のためのＩ−ＳＶＤルーチンの処理の手順を示す。図１の処理は、コンピュータプログラムの形式で提供され得る。

ステップ１０２において、一般行列Ａに対してハウスホルダー法を用いて上２重対角化を行う。このとき利用される上２重対角化方法は、ハウスホルダー法であってもよいし、他の上２重対角化方法であってもよい。これにより、上２重対角行列Ｂが得られる。次に、行列Ｂの特異値σ_ｊ（σ_１≧σ_２≧・・・≧σ_ｍ≧０、ｍは、行列Ｂ^ＴＢの行列サイズに等しい）を計算する。ここで、本発明により提供される特異値計算ルーチンｄＬＶｓおよびＬＡＰＡＣＫにより提供されるＤＬＡＳＱルーチンを利用し得る。なお、行列Ｂの特異値は、行列Ａの特異値に等しく、かつ、行列Ｂ^ＴＢの固有値の正の平方根λ_ｊ ^１／２に等しい。ｄＬＶｓルーチンは、Ｂの特異値σ_ｊを高精度で計算する。ｄＬＶｓルーチンの詳細は、後述される。ｄＬＶｓで求められる特異値は、従来のＬＡＰＣＫルーチン（ＤＢＤＳＱＲ、ＤＬＡＳＱ等）と比較して、最も高精度である。速度に関しては、ほとんどの場合ＤＬＡＳＱルーチンには劣るが、信頼性の高いＤＢＤＳＱＲルーチンの約２分の１の時間で計算できる(ＣＰＵがＩｔａｎｉｕｍ２の場合はｄＬＶｓルーチンが最速)。そこで、ステップ１０４において、高精度性が必要とされるかどうかを判定する。最も高精度な特異値分解を行いたい場合には、ステップ１０６においてｄＬＶｓを、高精度性および高速性を両立したい場合には、ステップ１０８において従来のＤＬＡＳＱルーチンを、特異値計算ルーチンとして用いるのがよい。さらに、ＤＬＡＳＱルーチンが有限回の演算で停止するかは定かでないが、ｄＬＶｓルーチンは収束することが保証されているので、信頼性を重視するならば、ｄＬＶｓルーチンを用いるのがよい。

次に、ステップ１１０において、計算されたＢの特異値σ_ｊ（すなわち、Ａの特異値に等しい）用いて、１．１において説明されるｓｄＬＶｖｓ変換およびｒｄＬＶｖｓ変換によるツイスティッド分解により、Ｂ^ＴＢの対角化を行う。本発明による行列Ｂ^ＴＢの対角化ルーチンの詳細は、１．１で説明する。さらに、ステップ１１２において、以下の１．１にて詳細が示される簡略化連立１次方程式

を求めると、行列Ｂの特異ベクトルｖ_ｊが得られる（すなわち、Ｂ^ＴＢの対角化行列Ｖが得られる）。すなわち、Ｂ^ＴＢが対角化される。ステップ１１４において、Σ≡ｄｉａｇ（σ_１，σ_２，・・・，σ_ｍ）（ただし、σ_１≧σ_２≧・・・≧σ_ｍ≧０）とすると、行列Ｂの右直交行列Ｖ_Ｂ＝Ｖであるので、行列Ｂの左直交行列Ｕ_Ｂ＝ＢＶ_ＢΣ^−１＝ＢＶΣ^−１からＵ_Ｂを得る。すなわち、行列Ｂが特異値分解される。ここでは、ｍ個全ての特異値および特異ベクトルを計算したが、この原理は、少なくとも１つの特異値および特異ベクトルを計算する場合にも適用され得る。

このとき得られる左右直交行列、すなわち特異ベクトル計算の結果は、ＤＢＤＳＱＲルーチンによって得られる結果と比較して精度の点で劣る。そこで、ステップ１１６において、逆反復法を用いることにより、ＤＢＤＳＱＲルーチンと同程度の精度を得ることができる。さらなる高精度性が必要とされる場合は、ステップ１２０において、グラムシュミット法による再直交化を行えば、ＤＢＤＳＱＲよりも高精度での計算を実現し得る。Ｉ−ＳＶＤルーチンでは、速度に関しては、表２から、ＤＢＤＳＱＲと比較して大幅に短縮され得ることが理解される。また、表１から理解されるように、Ｉ−ＳＶＤルーチンは、ＤＳＴＥＧＲによる極度の精度悪化がみられる場合でも、高精度に特異値分解することができる。なお、表１、２の誤差および計算時間に関するデータは、逆反復法を用いているが、グラムシュミット法は用いていない場合の結果である。また、ステップ１１６および１２０を実行する順序は逆であってもよいし、あるいは、両ステップは、省略してもよい。

このように、図１に示される各ステップの機能は、ソフトウェア（例えば、プログラム）によって実現される。しかし、本発明は、これに限定されない。図１に示される各ステップの機能をハードウェア（例えば、回路、ボード、半導体チップ）によって実現してもよいし、ソフトウェアとハードウェアとの組み合わせによって実現してもよい。

以上により、Ｉ−ＳＶＤルーチンを用いると、従来技術と比較して高速かつ高精度の特異値分解が実行される。

表１は、本発明による行列Ｂ^ＴＢ対角化ルーチンを実行した後に、逆反復法を実行したときの従来技術および本発明の誤差を比較した表である。

表２は、本発明による行列Ｂ^ＴＢ対角化ルーチンを実行した後に、逆反復法を実行したときの従来技術および本発明の計算時間を比較した表である。

また、１.３に、本発明の方法において、ＤＳＴＥＧＲルーチンと比較して、誤差が低減される理由を説明する。

１．１．本発明による行列Ｂ ^ＴＢの対角化ルーチン
ＤＢＤＳＱＲはＢ^ＴＢの固有値λ_ｊとともに対応する固有ベクトルｖ_ｊを求めることができるので、他のルーチンを併用しなくてもＢ^ＴＢの対角化（Ｂの特異値分解）が得られる。しかしながら、ＤＬＡＳＱおよびｄＬＶｓは上述したように固有ベクトルｖ_ｊを計算する機能はない。よって、別の固有ベクトル計算ルーチンが必要となる。ただし、Ｂ^ＴＢの固有値λ_ｊ（Ｂの特異値σ_ｊ）が求まっているものとする。

まず、ＤＳＴＥＧＲによる固有ベクトルの計算方法を示す。Ｐａｒｌｅｔｔらは以下の（１）、（２）、（３）の手順で固有ベクトル計算ができることを示している。
（１）ｓｔａｔｉｏｎａｒｙｑｄｗｉｔｈｓｈｉｆｔ変換（ｓｔｑｄｓ変換）およびｒｅｖｅｒｓｅ−ｔｉｍｅｐｒｏｇｒｅｓｓｉｖｅｑｄｗｉｔｈｓｈｉｆｔ変換（ｒｐｑｄｓ変換）を用いた（Ｂ^（０））^ＴＢ^（０）−λ_ｊＩのＣｈｏｌｅｓｋｙ（コレスキー）分解（Ｂ^（±１））^ＴＢ^（±１）＝（Ｂ^（０））^ＴＢ^（０）−λ_ｊＩを求める。ただし、Ｂ^（０）＝Ｂであり、Ｂ^（＋１）およびＢ^（−１）はそれぞれ上２重対角および下２重対角行列とする。すなわち

となる。すなわち、｛ｑ_ｋ ^（±１），ｅ_ｋ ^（±１）｝を計算する。
（２）コレスキー分解より（Ｂ^（０））^ＴＢ^（０）−λ_ｋＩのツイスティッド分解

を形成する。ただし

とし、γ_ｋ≡ｑ_ｋ ^（１）＋ｑ_ｋ ^（−１）−（ｅ_ｋ−１ ^（０）＋ｑ_ｋ ^（０）−λ_ｋ）が最小になるｋをｋ_０とする。
（３）簡約化された連立１次方程式

を解く。

図２は、Ｐａｒｌｅｔｔらの手法で最も核となる（１）の処理手順を示す。（１）で求めたいのは、｛ｑ_ｋ ^（ｎ），ｅ_ｋ ^（ｎ）｝，ｎ＝±１であるが、そのための変換としてｓｔａｔｉｏｎａｒｙｑｄｗｉｔｈｓｈｉｆｔ（ｓｔｑｄｓ）変換
ｑ_ｋ ^（１）＋ｅ_ｋ−１ ^（１）＝ｑ_ｋ ^（０）＋ｅ_ｋ−１ ^（０）−λ_ｊ，ｋ＝１，２，・・・，ｍ，
ｑ_ｋ ^（１）ｅ_ｋ ^（１）＝ｑ_ｋ ^（０）ｅ_ｋ ^（０），ｋ＝１，２，・・・，ｍ−１，
ｅ_０ ^（０）≡０，ｅ_０ ^（１）≡０
およびｒｅｖｅｒｓｅ−ｔｉｍｅｐｒｏｇｒｅｓｓｉｖｅｑｄｗｉｔｈｓｈｉｆｔ（ｒｐｑｄｓ）変換
ｑ_ｋ ^（−１）＋ｅ_ｋ ^（−１）＝ｑ_ｋ ^（０）＋ｅ_ｋ−１ ^（０）−λ_ｊ，ｋ＝１，２，・・・，ｍ，
ｑ_ｋ＋１ ^（−１）ｅ_ｋ ^（−１）＝ｑ_ｋ ^（０）ｅ_ｋ ^（０），ｋ＝１，２，・・・，ｍ−１，
ｅ_０ ^（０）≡０，ｅ_ｍ ^（−１）≡０
が採用されている。

固有値λ_ｊが既知ならば反復的な計算は不要なため計算量がＱＲアルゴリズムに比べて圧倒的に少なく抑えられるが、常に数値安定かつ精度のよい方法とはいえない。ｓｔｑｄｓ変換とｒｐｑｄｓ変換ともに減算による桁落ちが発生する可能性がある。例えば、ｓｔｑｄｓ変換においてｑ_ｋ ^（０）＋ｅ_ｋ−１ ^（０）−λ_ｊ〜ｅ_ｋ−１ ^（１）ならば、ｑ_ｋ ^（１）（＝ｑ_ｋ ^（０）＋ｅ_ｋ−１ ^（０）−λ_ｊ−ｅ_ｋ−１ ^（１））を求める際、倍精度計算であってもｑ_ｋ ^（１）の有効桁はわずか１桁になることもある。その場合、ｑ_ｋ ^（０）ｅ_ｋ ^（０）／ｑ_ｋ ^（１）を計算すると誤差が生じる、つまり、ｅ_ｋ ^（１）が精度よく計算できないことになる。また、ｑ_ｋ＋１ ^（１）を求めるのにｅ_ｋ ^（１）が必要、ｅ_ｋ ^（１）を求めるのにｑ_ｋ ^（１）が必要と逐次的な計算が要求されるので、１個所で発生した桁落ちによる誤差が波及し、さらなる誤差増大の可能性も秘めている。その結果、理論上はｑ_ｋ ^（１）≠０であるが誤差蓄積によりｑ_ｋ ^（１）＝０となり、ｑ_ｋ ^（０）ｅ_ｋ ^（０）／ｑ_ｋ ^（１）の計算においてオーバーフローが起こるといった数値不安定な状況も想定される。
Ｂ^（０）＝Ｂの成分｛ｂ_２ｋ−１ ^（０），ｂ_２ｋ ^（０）｝が与えられる、すなわち｛ｑ_ｋ ^（０），ｅ_ｋ ^（０）｝が与えられると、λ_ｊおよびｅ_ｋ−１ ^（１）が一意的に決まるので、この状況は避けることはできない。ｒｐｑｄｓ変換も同様の性質を持つため、実用的なレベルにまで達したとはいいがたい。ＬＡＰＡＣＫにおいてＦＯＲＴＲＡＮルーチンＤＳＴＥＧＲとして改良版が公開されているものの欠点は完全に解決されてはいない。

次に、本発明による対角化のルーチンを説明する。本発明によるルーチンでは、ｓｔｑｄｓ変換およびｒｐｑｄｓ変換にかわる新たな変換を採用している。本発明では、行列Ｂ^ＴＢを対角化するために、ｓｄＬＶｖｓ変換およびｒｄＬＶｖｓ変換を用いた行列Ｂ^ＴＢ−σ_ｊ ^２Ｉに対するツイスティッド分解を行う。正確には、上記（１）を以下の（１ａ）、（１ｂ）、（１ｃ）の３つの工程にわけることで数値安定なコレスキー分解を実現している。

図３〜図８は、本発明によるミウラ型逆変換、ｓｄＬＶｖｓ変換、ｒｄＬＶｖｓ変換およびミウラ型変換を実行する処理手順を示す。
（１ａ）ミウラ型逆変換

を行う。

ｕ_２ｋ−１ ^（０）＝η^（０）ｔ_ｋ ^（０），ｋ＝１，２，・・・，ｍ，
ｕ_２ｋ ^（０）＝ｅ_ｋ ^（０）／ｔ_ｋ ^（０），ｋ＝１，２，・・・，ｍ−１，
ｔ_ｋ ^（０）≡ｑ_ｋ ^（０）／（η^（０）＋ｕ_２ｋ−２ ^（０））−１，η^（０）≡１／δ^（０）
ただし、δ^（０）は任意に選ぶことができる。
（１ｂ）ｓｔａｔｉｏｎａｒｙ不等間隔差分Ｌｏｔｋａ−Ｖｏｌｔｅｒｒａ（ｓｄＬＶｖｓ）変換
ｕ_ｋ ^（１）＝ｕ_ｋ ^（０）＊（１＋δ^（０）ｕ_ｋ−１ ^（０））／（１＋δ^（１）ｕ_ｋ−１ ^（１）），
ｋ＝１，２，・・・，２ｍ−１，
によって

ｒｅｖｅｒｓｅ−ｔｉｍｅ不等間隔差分Ｌｏｔｋａ−Ｖｏｌｔｅｒｒａ（ｒｄＬＶｖｓ）変換
ｕ_ｋ ^（−１）＝ｕ_ｋ ^（０）＊（１＋δ^（０）ｕ_ｋ−１ ^（０））／（１＋δ^（−１）ｕ
_ｋ＋１ ^（−１）），
ｋ＝１，２，・・・，２ｍ−１，
によって

を実行する。ただし、ｕ_ｋ ^（ｎ）≡０，ｕ_２ｍ ^（０）≡０とし、λ_ｊ＝１／δ^（０）−１／δ^（ｎ），ｎ＝±１を満たすようにδ^（ｎ），ｎ＝±１を設定する。
（１ｃ）ミウラ型変換

ｑ_ｋ ^（ｎ）＝η^（ｎ）＊（１＋δ^（ｎ）ｕ_２ｋ−２ ^（ｎ））＊（１＋δ^（ｎ）ｕ_２ｋ−１ ^（ｎ）），
ｋ＝１，２，・・・，ｍ，
ｅ_ｋ ^（ｎ）＝δ^（ｎ）＊ｕ_２ｋ−１ ^（ｎ）＊ｕ_２ｋ ^（ｎ）
ｋ＝１，２，・・・，ｍ−１，
を行う。

ｑｄ型変換では見られない離散Ｌｏｔｋａ−Ｖｏｌｔｅｒｒａ型変換の大きな特徴は、任意パラメータを持つことである。すなわち、λ_ｊ＝１／δ^（０）−１／δ^（±１）を満たす範囲でδ^（ｎ）の値を任意に設定できる。δ^（ｎ）を変動させると補助変数ｕ_ｋ ^（ｎ）の値も変化するが、桁落ちによる誤差や数値不安定が発生するかどうかは事前に判定できる。この判定は、ｉｆ文によって実装されてもよい。この場合は、δ^（ｎ）を再設定後に再度計算すればよい。また、ｕ_ｋ ^（±１）が求まればミウラ変換によってｑ_ｋ ^（±１）およびｅ_ｋ ^（±１）が独立に計算されるので、誤差が伝播しないという性質を持つ。なお、ミウラ型逆変換をミウラ変換、ミウラ型変換を逆ミウラ変換、ｓｄＬＶｖｓ変換をｓｔＬＶｖ変換、ｒｄＬＶｖ変換をｒＬＶｖ変換と呼び変えても問題はない。

入力および出力の変数の対応は、Ｐａｒｌｅｔｔの手法と同じである。メモリ消費を抑えるために、補助変数ｕ_ｋ ^（ｎ）のための配列は必ずしも用意する必要はない。一方、１＋δ^（０）ｕ^（０）のためのメモリ領域を確保し、（１ａ）〜（１ｃ）のステップにまたがってこの値を利用することで、メモリ消費を抑え、計算量を低減することができる。これにより、誤差も低減される。このように、１．３で説明する方法を利用して誤差を低減すること等、さらなる工夫を行い得る。

上記の方法でＢ^ＴＢ−λ_ｊＩのコレスキー分解を求めた後は、Ｐａｒｌｅｔｔらと同じ方法を用いて特異ベクトルｖ_ｊを計算する。そうすると、Ｂ^ＴＢの固有値分解Λ＝Ｖ_Ｂ ^ＴＢ^ＴＢＶ_Ｂが得られる。ただしＶ_Ｂの列ベクトルがｖ_ｊである。Ｕ_Ｂ＝ＢＶ_ＢΣ^−１よりＵ_Ｂを求めると、Ｂの特異値分解Ｂ＝Ｕ_ＢΣＶ_Ｂ ^Ｔも得られる。

１．２．高精度特異値算出ｄＬＶｓルーチン
本発明によるｄＬＶｓルーチンは、基本的な枠組みではＤＬＡＳＱルーチンと同じであり、Ｂ^ＴＢの固有値（あるいはＢの特異値）のみが得られＢ^ＴＢの固有ベクトル（あるいはＢの特異ベクトル）は求まらない。ただ、ルーチンの核となる部分に使用される計算法はＤＬＡＳＱと異なる。

まず、ＤＬＡＳＱルーチンによる固有値の計算方法を説明する。ｑｄ法は固有ベクトルが計算できないが、１ステップは上２重対角行列の要素を更新するのみなので高速に特異値を求めることができる。かつ、計算量が少ないため丸め誤差の蓄積が抑えられ高精度計算が期待できる。ここで、Ｙ^（ｎ）を対称な３重対角行列

とする。また、変換Ｙ^（ｎ）→Ｙ^{（ｎ＋１）}がｑｄ法の漸化式
ｑ_ｋ ^{（ｎ＋１）}＝ｑ_ｋ ^（ｎ）＋ｅ_ｋ ^（ｎ）−ｅ_ｋ−１ ^{（ｎ＋１）}，ｋ＝１，２，・・・，ｍ，
ｅ_ｋ ^{（ｎ＋１）}＝ｅ_ｋ ^（ｎ）＊ｑ_ｋ＋１ ^（ｎ）／ｑ_ｋ ^{（ｎ＋１）}，ｋ＝１，２，・・・，ｍ−１，
ｅ_０ ^（ｎ）≡０，ｅ_ｍ ^（ｎ）≡０，ｎ＝０，１，・・・
によって実現されるとする。そのとき、適切なＲ^（ｎ）が存在してＹ^{（ｎ＋１）}＝Ｒ^（ｎ）Ｙ^（ｎ）（Ｒ^（ｎ））^−１となることをＲｕｔｉｓｈａｕｓｅｒは発見している。これは、Ｙ^{（ｎ＋１）}とＹ^（ｎ）の固有値が同じであることを意味し、すなわち、上記の漸化式によって固有値保存変形がなされる。この変形を繰り返せば非対角成分が０に近づき、Ｙ^（ｎ）が対角化されることも証明されている。よってＹ^（０）＝Ｂ^ＴＢとするとｌｉｍ
_ｎ→∞Ｙ^（ｎ）＝ｄｉａｇ（λ_１，λ_２，・・・，λ_ｍ）となる。ただし、λ_ｋはＢ^ＴＢの固有値である。さらに、λ_ｋの正の平方根をとるとＢの特異値σ_ｋが得られる。

ＬＡＰＣＡＫの特異値計算ルーチンＤＬＡＳＱでは、ｑｄ法のｄｉｆｆｅｒｅｎｔｉａｌ型が採用されている。これはｄｑｄ法と呼ばれる。１ステップは

で与えられ、変換^）Ｙ^（ｎ）→Ｙ^{（ｎ＋１）}が実現される。ｄｉｆｆｅｒｅｎｔｉａｌ型は減算を含まないため、ｑｄ法では起こりうる桁落ち誤差の心配もない。さらに高速版として原点シフト付きｄｑｄ（ｄｑｄｓ）法も組み込まれている。高速化に伴う計算量の減少で丸め誤差の発生も抑えることができる。ただし、収束するかどうかの保証はない。原点シフト版の１ステップは

で与えられる。ｓは原点シフト量であり、Ｇｅｒｓｇｏｒｉｎ型境界条件により見積もられたＹ^（ｎ）の最小固有値とする。よって、ＤＬＡＳＱルーチンのメインループの核となる部分では最小固有値の見積もりを行い、ｓ＞ε（ε：非常に小さいとみなせる正の数）ならばｄｑｄｓ法が、そうでなければｄｑｄ法が用いられる。メインループの他の部分（行列の分割・減次）はほとんどＤＢＤＳＱＲルーチンと変わらない。すなわち、ＤＢＤＳＱＲルーチンに含まれるＱＲ法をｄｑｄ法に、原点シフトＱＲ法をｄｑｄｓ法に置き換えたものがＤＬＡＳＱルーチンとなる。（ＱＲ法とｄｑｄ法、原点シフト付きＱＲ法とｄｑｄｓ法で扱う変数が異なるので少しだけ判定式が異なる部分もあるが本質的に同じである。）
次に本発明によるｄＬＶｓルーチンを説明する。ｄＬＶｓルーチンではｄｑｄ法のかわりにｄＬＶ法を、ｄｑｄｓ法のかわりにｄＬＶｓ法を採用する。ここで、ｄＬＶｓルーチンはＤＬＡＳＱルーチンとは違って収束することが保証されている。その１ステップは以下で与えられる。

ＤＬＡＳＱルーチン同様、メインループの核の部分では最小特異値の見積もりによって原点シフト量ｓを決定し、ｄＬＶ法を使うかｄＬＶｓ法を使うかを決定する。ただし、補助変数ｖ_ｋ ^（ｎ）を求めたあとでｓを決定する点はＤＬＡＳＱルーチンと異なる。

１．３．本発明の方法による誤差の低減
人間が理想的な状況、すなわち、無限桁の計算をいくらでもできるとすると、Ｐａｒｌｅｔｔらの手法、本発明による特異ベクトル計算ルーチン（１．１参照）のどちらを使っても問題ない。しかしながら、コンピュータで計算を行う場合は注意が必要である。有限桁の計算しか行えないコンピュータ上では、数学的に正しい計算法を使ったとしても必ずしも正しい結果が得られる訳ではない。そればかりかいつまでたっても計算が終了しないといった思わぬ数値的な問題が発生する場合もある。

コンピュータ計算による誤差には、丸め誤差および桁落ちによる誤差などが知られている。丸め誤差単独では、高々有効桁の最後の桁が真値と比べて異なる程度で大きな誤差にはならない。また、指数部が異なる２つの実数の加算、乗算、除算の少なくとも１つの演算を行えばやはり丸め誤差が生じるが、それ以上の誤差は発生しない。さらに、このような丸め誤差を発生するような操作が繰り返されても、丸めモードがｎｅａｒ（四捨五入）ならば、一方的に切り上げられたりあるいは切り捨てられたりして誤差が極端に蓄積することは少ない。よって、多くの数値計算法は加算、乗算、除算の少なくとも１つの演算によって発生する丸め誤差を特別注意することは少なく、ｄＬＶｓルーチンによる特異値計算でも結果的に丸め誤差は一様に増大しない。

問題となるのは、同符号の実数の減算および異符号の実数の加算により生じる、桁落ちである。桁落ちによる誤差で値が０となった後、その値による除算を行うと、０が分母にくるような不定形となり計算不可能となる。こうなるといつまでたっても計算が終了しない。減算→除算と計算する部分がＰａｒｌｅｔｔの手法、本発明による方法ともに存在するので、桁落ち誤差には十分に注意する必要がある。

本発明による計算方法では、上述の桁落ちによる誤差を含んでいるかどうかは減算によって得られた値が小さいかどうかで判断できる。Ｐａｒｌｅｔｔの手法の場合、図２のＤＯ文中のｑ１［ｋ］、ｑ２［ｋ］をチェックすればよい。ところが、Ｐａｒｌｅｔｔの手法は桁落ち誤差を含むことが分かったとしても、それを回避することはできない。なぜならば、初期値として｛ｑ０[ｋ]，ｅ０[ｋ]｝が与えられると、λは一意的に決定され、｛ｑ１[ｋ]，ｅ１[ｋ]｝および｛ｑ２[ｋ]，ｅ２[ｋ]｝も一意的に導出されるためである。すなわち、任意パラメータを持たない自由度のない計算法であるためである。

それに対して本発明による方法は自由に設定できるパラメータη１（＝１／δ^（０））をもつため、補助変数ｕ_ｋ ^（ｎ）の値を様々に変化させることができる（図９参照）。すなわち、様々な経路で｛ｑ１[ｋ]，ｅ１[ｋ]｝および｛ｑ２[ｋ]，ｅ２[ｋ]｝を計算することができる。よって、桁落ちが発生する場合も回避できる。図６〜図８の条件判定によって桁落ちの影響をチェックし、減算によって得られた値の絶対値が小さな数εより大きいという条件が満たされなければ、パラメータη１の設定に戻るというものである。この処理は、条件が満たされるまで繰り返される。なお、精度よりも高速性を重視する場合は、数回条件が満たさなければ（ｑ１[ｋ]＝０あるいはｑ２[ｋ]＝０ならば）、例外処理（文献ＰＤＦ参照）を行ってもよい。

２．本発明による特異値分解装置
特異値分解を実行する装置を説明する。特異値分解を実行する装置は、例えば、図１に示される各ステップの機能を有する装置である。特異値分解を実行する装置の１つの実施形態は、特異値分解プログラムを実行するコンピュータである。そのプログラムは、例えば、図１に示される各ステップをＣＰＵに実行させるプログラムである。

図１０は、図１の処理をコンピュータプログラムとして実行するコンピュータ１０を示す。コンピュータ１０は、ＣＰＵ１００１と、メモリ１００２と、入力インターフェース部１００３と、出力インターフェース部１００４と、バス１００５とを含む。

ＣＰＵ１００１は、プログラムを実行する。例えば、そのプログラムは、図１の処理を実行するプログラムである。そのプログラムおよびそのプログラムの実行に必要なデータは、例えば、メモリ１００２に格納される。そのプログラムは、任意の様態でメモリ１００２に含まれ得る。例えば、メモリ１００２が書き換え可能なメモリである場合には、コンピュータ１００２の外部からプログラムをローディングすることにより、そのプログラムをメモリ１００２に格納してもよい。あるいは、メモリ１００２が読み出し専用メモリである場合には、メモリ１００２に焼き付ける形式でそのプログラムをメモリ１００２に格納してもよい。

さらに、入力インターフェース部１００３は、特異値分解を実行する対象である行列Ａを外部から受け取るためのインターフェースとして機能する。出力インターフェース部１００４は、計算結果を出力するためのインターフェースとして機能する。バス１００５は、コンピュータ１０内の構成要素１００１〜１００４を相互に接続するために使用される。

３．本発明による特異値分解方法の応用
本発明による特異値分解方法は、様々な分野に応用可能である。以下に、２次元画像から３次元画像へ復元する画像処理への応用例、ならびに、文書検索への応用例を示す。しかし、これら２つの応用例は例示に過ぎず、本発明による特異値分解方法の応用は、これら２つの応用例に限定されない。本発明による特異値分解方法の応用は、特異値分解を利用するあらゆる分野に適用可能である。

３．１．２次元から３次元へ復元する画像処理への応用
図１１を参照する。図１１は、本発明による特異値分解方法を利用した物体の複数の２次元画像を３次元画像へ復元する画像処理の１つの実施形態を説明する。

複数の２次元画像から３次元復元を行うために必要となるステップは、２次元画像から特徴点を抽出するステップと、特徴点データより形状（元の物体の特徴点の３次元座標データ）および回転（３次元データから特徴点データへの変換）に関するデータを計算するステップと、形状および回転のデータより可視化を行うステップとを包含する。

ステップ１１０２において、２次元画像ｊ（ｊ＝１，・・・，ｍ、ｍは３以上の整数）から特徴点ｉ（ｉ＝１，・・・，ｎ、ｎは２以上の整数）の座標（ｘ_ｉ ^ｊ，ｙ_ｉ ^ｊ）を抽出する。取り扱う２次元画像は、弱中心射影画像であることが好ましい。このとき

が成り立つ。ここで、ｓ_ｊは物体のスケールに相対するｊ番目の画像のスケール、ｒ_１，ｊ，ｒ_２，ｊはそれぞれ物体座標系に相対するｊ番目のカメラ座標系の回転行列の１番目と２番目の行ベクトル、［Ｘ_ｉ，Ｙ_ｉ，Ｚ_ｉ］^Ｔはｉ番目の点の３次元座標である。物体のスケールは１番目の画像のスケールと同じにし（ｓ_１＝１）、物体の座標系の姿勢は１番目の画像のカメラ座標系と同じにする（ｒ_１，１＝［１，０，０］^Ｔ，ｒ_２，１＝［０，１，０］^Ｔ）。ｍ枚全ての画像におけるｎ個全ての点の座標を行列Ｄの要素として並べると、
Ｄ＝ＭＳ
が得られる。
ただし、

である。ＭとＳの形から分かるように、Ｄのランクは３である。いま、ステップ１１０２よりＤが与えられている。以下、回転に関するデータＭおよび形状Ｓを求める。

そこで、行列Ｄの特異値分解
Ｄ＝ＵΣＶ^Ｔ
を考える。ここで、Σは特異値を大小順に対角線上に並べたもので、ＵおよびＶはそれぞれ左および右直交行列である。ステップ１１０４では、行列Ｄの特異値分解を計算するために、行列Ｄを上２重対角化して上２重対角化行列Ｂを得る。ステップ１１０６において、行列Ｂ、すなわちＤの特異値を計算する。ステップ１１０６では、本発明によるｄＬＶｓルーチン、ＤＬＡＳＱルーチン、その他の特異値計算ルーチン、あるいはその組み合わせを利用し得る。このとき、画像のデジタル誤差のため、ゼロでない特異値は３つ以上出てくる。しかし、４番目以降の特異値はノイズによるもので、最初の３つの特異値と比べて格段に小さい。

そこで、ステップ１１０８において、最初の３つの特異値に対して特異ベクトルを計算する。ステップ１１０８は、図１のステップ１１０〜１２０を利用し得る。採用する３個のベクトルをまとめると、次式となる。

Ｄ’＝Ｌ’Σ’Ｒ’^Ｔ＝Ｍ’Ｓ’
ただし、Ｍ’＝Ｌ’（Σ’）^１／２、Ｓ’＝（Σ’）^１／２Ｒ’^Ｔ、Ｄ’は‖Ｄ−Ｄ’‖を最小にするランク３の行列である。

次に、Ｄ’からＭおよびＳを求めたいが、その組合せは唯一ではない。なぜなら、任意の正則行列Ｃが
Ｄ’＝（Ｍ’Ｃ）（Ｃ^−１Ｓ’）
を満たすからである。そこで、上式におけるＣをＭ＝Ｍ’Ｃを満たすように決める。Ｃは下記の式を満たす。

Ｅ＝ＣＣ^Ｔとすると、上式からＥの６つの要素に関する２ｍ＋１個の線形方程式が得られる。ｍ≧３であるので、Ｅの要素を一意に決めることができる。ステップ１１１０において、行列Ｅを求める。

次に、ステップ１１１２において、ＥからＣを求める。Ｃの自由度（９）はＥの自由度（６）より多い。そこで、条件ｒ_１ｊ＝［１，０，０］^Ｔ，ｒ_２ｊ＝［０，１，０］^Ｔを加えれば、Ｃを決めることができる。このとき２つの解（ＮｅｃｋｅｒＲｅｖｅｒｓａｌ）が出る。

次に、ステップ１１１４において、Ｍ＝Ｍ’ＣおよびＳ＝Ｃ^−１Ｓ’より、回転に関するデータＭおよび形状Ｓが決まる。

ＤＢＤＳＱＲルーチンでは上記のようにＵ、Ｖ（ｎ個のベクトル）を求めた後、３個のベクトル（ｎ−３個は不要）を採用する。ここでは、原則的にはｎ個のベクトルを求めるが、例外的にｎ個より少ない数のベクトルを求めてもよい。この実施形態では３個のベクトルだけ求めればよい。すなわち、計算コストを削減することができる。

３．２．本発明による特異値分解方法を利用した文書検索方法
文書中からその文書の内容に関連する索引語を抽出し、索引語の重みを計算する処理の後、ベクトル空間モデルでは、この索引語の重みを要素とするベクトルで文書を表現する。いま、検索対象となる文書をＤ_１，Ｄ_２，・・・，Ｄ_ｎとし、これら文書集合全体を通して全部でｍ個の索引語ｗ_１，ｗ_２，・・・，ｗ_ｍがあるとする。このとき、文書Ｄ_ｊは、次のようなベクトルで表現されることになる。これを文書ベクトルと呼ぶ。

ここで、ｄ_ｉｊは索引語ｗ_ｉの文書Ｄ_ｊにおける重みである。また、文書集合全体は、次のようなｍ×ｎ行列Ｄによって表現することができる。

行列Ｄを索引語文書行列と呼ぶ。索引語文書行列の各列は文書に関する情報を表している文書ベクトルであるが、同様に、索引語文書行列の各行は索引語に関する情報を表しているベクトルであり、これを索引語ベクトルと呼ぶ。検索質問も、文書と同様に、索引語の重みを要素とするベクトルで表現することができる。検索質問文に含まれる索引語ｗ_ｉの重みをｑ_ｉとすると、検索質問ベクトルｑは次のように表されることになる。

実際の文書検索においては、与えられた検索質問文と類似した文書を見つけ出す必要があるが、検索質問ベクトルｑと各文書ベクトルｄ_ｊの間の類似度を計算することにより行う。ベクトル間の類似度の定義としてはさまざまなものが考えられるが、文書検索においてよく用いられているものはコサイン尺度（２つのベクトルのなす角度）または内積である。

なお、ベクトルの長さが１に正規化（コサイン正規化）されている場合には、コサイン尺度と内積とは一致する。

図１２を参照する。図１２は、本発明による特異値分解装置を利用した、文書検索方法の１つの実施形態を説明する。

ステップ１２０２において、質問ベクトルｑを受け取る。次に、ステップ１２０４では、行列Ｄの特異値分解を計算するために、行列Ｄを上２重対角化して上２重対角化行列Ｂを得る。ステップ１２０６では、行列Ｂ、すなわち行列Ｄの特異値を求める。ステップ１２０６では、本発明によるｄＬＶｓルーチン、ＤＬＡＳＱルーチン、その他の特異値計算ルーチン、あるいはその組み合わせを利用し得る。

次に、Ｄの近似行列Ｄ_ｋを使った検索を考える。ベクトル空間モデルでは、検索質問ベクトルｑと索引語文書行列Ｄ中の各文書ベクトルｄ_ｊの間の類似度を計算することにより検索を行うが、ここではＤの代わりにＤ_ｋを使う。ベクトル空間モデルでは、文書ベクトルの次元数は索引語の総数と等しい。したがって、検索対象となる文書の数が増えるに従い、文書ベクトルの次元数も増加する傾向にある。しかし、次元数が増加してくると、コンピュータのメモリによる制限や検索時間の増大などの問題が生じてくるばかりでなく、文書中に含まれる不必要な索引語がノイズ的な影響を及ぼし、検索精度を低下させてしまうという現象も起こってくる。潜在的意味インデキシング（ｌａｔｅｎｔｓｅｍａｎｔｉｃｉｎｄｅｘｉｎｇ；ＬＳＩ）は、高次元の空間にある文書ベクトルを低次元の空間へと射影することにより、検索精度の改善を図る技術である。高次元の空間では別々に扱われていた索引語が、低次元の空間では相互に関連を持ったものとして扱われる可能性もあるため、索引語の持つ意味や概念に基づく検索を行うことができる。たとえば、通常のベクトル空間モデルでは“ｃａｒ”という索引語と“ａｕｔｏｍｏｂｉｌｅ”という索引語はまったく別物であり、一方の索引語による質問ではもう片方の索引語を含んだ文書を検索することができない。しかし、低次元の空間ではこれらの意味的に関連した索引語は１つの次元に縮退することが期待できるため、“ｃａｒ”という検索質問によって“ｃａｒ”を含む文書ばかりでなく“ａｕｔｏｍｏｂｉｌｅ”を含む文書をも検索することが可能となる。潜在的意味インデキシングでは、特異債分解により高次元ベクトルの次元削減を行うが、これは基本的に多変量解析における主成分分析と等価である。

行列ステップ１２０８では、ｋ＜ｒなるｋの値を選択する。ｋの値は、予め与えられていてもよいし、計算ごとに選択可能であってもよい。

ステップ１２１０において、ステップ１２０６で計算した特異値のうち、大きい順に１番目からｋ番目のｋ個の特異値に対して、Ｄの特異ベクトルを計算する。すなわち、
Ｄ_ｋ＝Ｕ_ｋΣＶ_ｋ ^Ｔ
なるＵ_ｋおよびＶ_ｋを計算する。ここで、Ｕ_ｋは、最初のｋ個の左特異ベクトルのみから構成されるｍ×ｋ行列であり、Ｖ_ｋは、最初のｋ個の右特異ベクトルのみから構成されるｎ×ｋ行列であり、Σは、最初のｋ個の特異値のみから構成されるｋ×ｋ対角行列である。ステップ１２１０は、図１のステップ１１０〜１２０を利用し得る。

次に、ステップ１２１２において、行列Ｄ_ｋと質問ベクトルｑとの類似度を計算する。いま、ベクトルｅ_ｊをｎ次元の単位ベクトルとすると、Ｄ_ｋのｊ番目の文書ベクトルはＤ_ｋｅ_ｊで表すことができる。文書ベクトルＤ_ｋｅ_ｊと検索質問ベクトルｑとの間の類似度計算は、

としてもよいが、別の定義を用いてもよい。上式では、Ｄ_ｋをＵ_ｋ，Σ_ｋ，Ｖ_ｋから再構成する必要はなく特異値分解の結果から、直接、類似度を計算できることを示している。
数２９の中に現われるΣ_ｋＶ_ｋ ^Ｔｅ_ｊは、
Σ_ｋＶ_ｋ ^Ｔｅ_ｊ＝Ｕ_ｋ ^ＴＤ_ｋｅ_ｊ
と書き直すことができる。この式の右辺は、近似行列Ｄ_ｋにおけるｊ番目の文書ベクトルの基底Ｕ_ｋのもとでの座標（文書のｋ次元表現）を表している。同様に、数２９中のＵ_ｋ ^Ｔｑは、検索質問ベクトルｑの基底Ｕ_ｋのもとでの座標（検索質問のｋ次元表現）である。

ステップ１２１４において、ステップ１２１２において計算された類似度を基準に、検索結果を出力する。

ＤＢＤＳＱＲルーチンではｒ個のベクトルを求めた後、ｋ個のベクトル（ｒ−ｋ個は不要）を採用するのに対して、この実施形態では、ｋ個のベクトルだけ求めればよい。すなわち、計算コストを削減することができる。

以上のように、本発明の好ましい実施形態を用いて本発明を例示してきたが、本発明は、この実施形態に限定して解釈されるべきものではない。本発明は、特許請求の範囲によってのみその範囲が解釈されるべきであることが理解される。当業者は、本発明の具体的な好ましい実施形態の記載から、本発明の記載および技術常識に基づいて等価な範囲を実施することができることが理解される。本明細書において引用した特許、特許出願および文献は、その内容自体が具体的に本明細書に記載されているのと同様にその内容が本明細書に対する参考として援用されるべきであることが理解される。

本発明は、任意の行列Ａの特異値分解を高速かつ高精度に実行することを可能にする方法、プログラム、および装置等として有用である。

Claims

物体の複数の２次元画像から３次元画像を復元する方法であって、
２次元画像ｊ（ｊ＝１，・・・，ｍ、ｍは３以上の整数）における特徴点ｉ（ｉ＝１，・・・，ｎ、ｎは２以上の整数）の座標ｄ_ｉｊ（ｘ_ｉｊ，ｙ_ｉｊ）を抽出するステップと、
該特徴点の２次元座標ｄ_ｉｊ（ｘ_ｉｊ，ｙ_ｉｊ）より、該特徴点の３次元座標ｓ_ｉ（Ｘ_ｉ，Ｙ_ｉ，Ｚ_ｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算するステップと
を包含し、
該特徴点の３次元座標ｓ_ｉ（Ｘ_ｉ，Ｙ_ｉ，Ｚ_ｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算するステップは、
行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求めるステップであって、行列Ｄは、
と定義される、ステップと、
行列Ｄの特異値として行列Ｂの特異値σ_ｊ（σ_１≧σ_２≧・・・≧σ_ｒ＞０、ｒは、行列Ｄのランクに等しい）を求めるステップと、
σ_１、σ_２およびσ_３に対する行列Ｄの特異ベクトルを求めるステップと、
Ｍ＝Ｍ’Ｃなる行列Ｃに対して、Ｅ＝ＣＣ^Ｔを満たす行列Ｅを計算するステップであって、Ｍ’＝Ｌ’（Σ’）^１／２、Σ’はσ_１、σ_２およびσ_３を対角要素に持ちそれ以外の要素が０である３×３行列、Ｌ’はσ_１、σ_２およびσ_３に対応するＤの特異ベクトルを左から順に並べた行列である、ステップと、
行列Ｅから、行列Ｃを計算するステップと、
行列Ｃから、該３次元座標ｓ_ｉ（Ｘ_ｉ，Ｙ_ｉ，Ｚ_ｉ）および該変換を表す行列Ｍを計算するステップと
を含み、
該σ_１、σ_２およびσ_３に対する行列Ｄの特異ベクトルを求めるステップは、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列Ｂ^ＴＢ−σ_ｊ ^２Ｉに対してツイスティッド分解を実行することにより、行列Ｂ^ＴＢを対角化するステップを含み、Ｉは単位行列である、方法。
与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する方法であって、
該キーワードに対応する質問ベクトルｑを受け取るステップと、
該文書に対応する索引語文書行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求めるステップと、
行列Ｄの特異値として行列Ｂの特異値σ_ｊ（σ_１≧σ_２≧・・・≧σ_ｒ＞０、ｒは、行列Ｄのランクに等しい）を求めるステップと、
ｋ＜ｒなるｋを選択するステップと、
行列Ｄの擬似行列Ｄ_ｋを計算するステップであって、行列Ｄ_ｋは、σ_１，σ_２，・・・，σ_ｋを対角要素としそれ以外の要素が０である行列Σ_ｋ、σ_１，σ_２，・・・，σ_ｋに対応する特異ベクトルのみを左から順に並べた左右直交行列Ｕ_ｋ，Ｖ_ｋを用いて、Ｄ_ｋ＝Ｕ_ｋΣ_ｋＶ_ｋ ^Ｔと定義される、ステップと、
行列Ｄ_ｋと質問ベクトルｑとの類似度を計算するステップと、
該計算された類似度を基準に、検索結果を出力するステップと
を包含し、
該行列Ｄ_ｋの左右直交行列Ｕ_ｋ，Ｖ_ｋを求めるステップは、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列Ｂ^ＴＢ−σ_ｊ ^２Ｉ（ｊ＝１，２，・・・，ｋ）に対してツイスティッド分解を実行することにより、行列Ｂ^ＴＢを対角化するステップを含み、Ｉは単位行列である、方法。
物体の複数の２次元画像から３次元画像を復元する画像復元処理をコンピュータに実行させるプログラムであって、
該画像復元処理は、
２次元画像ｊ（ｊ＝１，・・・，ｍ、ｍは３以上の整数）における特徴点ｉ（ｉ＝１，・・・，ｎ、ｎは２以上の整数）の座標ｄ_ｉｊ（ｘ_ｉｊ，ｙ_ｉｊ）を抽出するステップと、
該特徴点の２次元座標ｄ_ｉｊ（ｘ_ｉｊ，ｙ_ｉｊ）より、該特徴点の３次元座標ｓ_ｉ（Ｘ_ｉ，Ｙ_ｉ，Ｚ_ｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算するステップと
を包含し、
該特徴点の３次元座標ｓ_ｉ（Ｘ_ｉ，Ｙ_ｉ，Ｚ_ｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算するステップは、
行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求めるステップであって、行列Ｄは、
と定義される、ステップと、
行列Ｄの特異値として行列Ｂの特異値σ_ｊ（σ_１≧σ_２≧・・・≧σ_ｒ＞０、ｒは、行列Ｄのランクに等しい）を求めるステップと、
σ_１、σ_２およびσ_３に対する行列Ｄの特異ベクトルを求めるステップと、
Ｍ＝Ｍ’Ｃなる行列Ｃに対して、Ｅ＝ＣＣ^Ｔを満たす行列Ｅを計算するステップであって、Ｍ’＝Ｌ’（Σ’）^１／２、Σ’はσ_１、σ_２およびσ_３を対角要素に持ちそれ以外の要素が０である３×３行列、Ｌ’はσ_１、σ_２およびσ_３に対応するＤの特異ベクトルを左から順に並べた行列である、ステップと、
行列Ｅから、行列Ｃを計算するステップと、
行列Ｃから、該３次元座標ｓ_ｉ（Ｘ_ｉ，Ｙ_ｉ，Ｚ_ｉ）および該変換を表す行列Ｍを計算するステップと
を含み、
該σ_１、_２およびσ_３に対する行列Ｄの特異ベクトルを求めるステップは、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列Ｂ^ＴＢ−σ_ｊ ^２Ｉに対してツイスティッド分解を実行することにより、行列Ｂ^ＴＢを対角化するステップを含み、Ｉは単位行列である、プログラム。
与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する文書検索処理をコンピュータに実行させるプログラムであって、
該文書検索処理は、
該キーワードに対応する質問ベクトルｑを受け取るステップと、
該文書に対応する索引語文書行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求めるステップと、
行列Ｄの特異値として行列Ｂの特異値σ_ｊ（σ_１≧σ_２≧・・・≧σ_ｒ＞０、ｒは、行列Ｄのランクに等しい）を求めるステップと、
ｋ＜ｒなるｋを選択するステップと、
行列Ｄの擬似行列Ｄ_ｋを計算するステップであって、行列Ｄ_ｋは、σ_１，σ_２，・・・，σ_ｋを対角要素としそれ以外の要素が０である行列Σ_ｋ、σ_１，σ_２，・・・，σ_ｋに対応する特異ベクトルのみを左から順に並べた左右直交行列Ｕ_ｋ，Ｖ_ｋを用いて、Ｄ_ｋ＝Ｕ_ｋΣ_ｋＶ_ｋ ^Ｔと定義される、ステップと、
行列Ｄ_ｋと質問ベクトルｑとの類似度を計算するステップと、
該計算された類似度を基準に、検索結果を出力するステップと
を包含し、
該行列Ｄ_ｋの左右直交行列Ｕ_ｋ，Ｖ_ｋを求めるステップは、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列Ｂ^ＴＢ−σ_ｊ ^２Ｉ（ｊ＝１，２，・・・，ｋ）に対してツイスティッド分解を実行することにより、行列Ｂ^ＴＢを対角化するステップを含み、Ｉは単位行列である、プログラム。
物体の複数の２次元画像を３次元画像に復元する装置であって、
ｍ枚（ｍは３以上の整数）の２次元画像を受け取る手段と、
２次元画像ｊ（ｊ＝１，・・・，ｍ）における特徴点ｉ（ｉ＝１，・・・，ｎ、ｎは２以上の整数）の座標ｄ_ｉｊ（ｘ_ｉｊ，ｙ_ｉｊ）を抽出する手段と、
該特徴点の２次元座標ｄ_ｉｊ（ｘ_ｉｊ，ｙ_ｉｊ）より、該特徴点の３次元座標ｓ_ｉ（Ｘ_ｉ，Ｙ_ｉ，Ｚ_ｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算する手段と
を備え、
該特徴点の３次元座標ｓ_ｉ（Ｘ_ｉ，Ｙ_ｉ，Ｚ_ｉ）および２次元座標から３次元座標への変換を表す行列Ｍを計算する手段は、
行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求める手段であって、行列Ｄは、
と定義される、手段と、
行列Ｄの特異値として行列Ｂの特異値σ_ｊ（σ_１≧σ_２≧・・・≧σ_ｒ＞０、ｒは、行列Ｄのランクに等しい）を求める手段と、
σ_１、σ_２およびσ_３に対する行列Ｄの特異ベクトルを求める手段と、
Ｍ＝Ｍ’Ｃなる行列Ｃに対して、Ｅ＝ＣＣ^Ｔを満たす行列Ｅを計算する手段であって、Ｍ’＝Ｌ’（Σ’）^１／２、Σ’はσ_１、σ_２およびσ_３を対角要素に持ちそれ以外の要素が０である３×３行列、Ｌ’はσ_１、σ_２およびσ_３に対応するＤの特異ベクトルを左から順に並べた行列である、手段と、
行列Ｅから、行列Ｃを計算する手段と、
行列Ｃから、該３次元座標ｓ_ｉ（Ｘ_ｉ，Ｙ_ｉ，Ｚ_ｉ）および該変換を表す行列Ｍを計算する手段と
を含み、
該σ_１、σ_２およびσ_３に対する行列Ｄの特異ベクトルを求める手段は、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列Ｂ^ＴＢ−σ_ｊ ^２Ｉに対してツイスティッド分解を実行することにより、行列Ｂ^ＴＢを対角化する手段を含み、Ｉは単位行列である、装置。
与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する装置であって、
該キーワードに対応する質問ベクトルｑを受け取る手段と、
該文書に対応する索引語文書行列Ｄに対して上２重対角化を行い、行列Ｄの上２重対角行列Ｂを求める手段と、
行列Ｄの特異値として行列Ｂの特異値σ_ｊ（σ_１≧σ_２≧・・・≧σ_ｒ＞０、ｒは、行列Ｄのランクに等しい）を求める手段と、
ｋ＜ｒなるｋを選択する手段と、
行列Ｄの擬似行列Ｄ_ｋを計算する手段であって、行列Ｄ_ｋは、σ_１，σ_２，・・・，σ_ｋを対角要素としそれ以外の要素が０である行列Σ_ｋ、σ_１，σ_２，・・・，σ_ｋに対応する特異ベクトルのみを左から順に並べた左右直交行列Ｕ_ｋ，Ｖ_ｋを用いて、Ｄ_ｋ＝Ｕ_ｋΣ_ｋＶ_ｋ ^Ｔと定義される、手段と、
行列Ｄ_ｋと質問ベクトルｑとの類似度を計算する手段と、
該計算された類似度を基準に、検索結果を出力する手段と
を包含し、
該行列Ｄ_ｋの左右直交行列Ｕ_ｋ，Ｖ_ｋを求める手段は、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列Ｂ^ＴＢ−σ_ｊ ^２Ｉ（ｊ＝１，２，・・・，ｋ）に対してツイスティッド分解を実行することにより、行列Ｂ^ＴＢを対角化する手段を含み、Ｉは単位行列である、装置。