WO2005119507A1

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

Info

Publication number: WO2005119507A1
Application number: PCT/JP2005/010084
Authority: WO
Inventors: Yoshimasa Nakamura; Masashi Iwasaki; Shinya Sakano
Original assignee: Japan Science and Technology Agency
Current assignee: Japan Science and Technology Agency
Priority date: 2004-06-03
Filing date: 2005-06-01
Publication date: 2005-12-15
Anticipated expiration: 2006-12-03
Also published as: CA2568852C; EP1752884A4; US20090028455A1; JP4325877B2; CA2568852A1; EP1752884A1; EP1752884B1; US8306361B2; JPWO2005119507A1

Abstract

　本発明による特異値分解方法は、コンピュータを用いて任意の行列Ａに対して特異値分解を実行する方法であって、行列Ａに対して上２重対角化を行い、行列Ａの上２重対角行列Ｂを求めるステップと、行列Ａの特異値として行列Ｂの特異値σを求めるステップと、σに対する行列Ａの特異ベクトルを求めるステップとを包含する。行列Ａの特異ベクトルを求めるステップは、ミウラ型逆変換と、ｓｄＬＶｖｓ変換と、ｒｄＬＶｖｓ変換と、ミウラ型変換とを用いて行列ＢＴＢ－σ２Ｉ（Ｉは単位行列）に対してツイスティッド分解を実行することにより、行列ＢＴＢを対角化するステップを含む。

Description

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

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

背景技術

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

[0003] LAPACKによる DBDSQR/レーチンを例として、従来の特異値分解方法を説明する。 DBDSQRルーチンでは、 1 X Iの一般行列 Aを特異値分解するために、まず、

1 2

行列 Aに対して Householder (ハウスホルダー）変換を適用する。すなわち、行列 A は、直交行列 U 、 Vを用いて、

[0004] [数 4]

UjAV_A =

U U_A = I, V_A = I, m = mm{£₁,£₂} と表わすことができる。このとき得られる Bを、上 2重対角行列という。ここで、 Aの特異値 =Bの特異値であることに注意する。このように、一般行列 Aに対する特異値分解問題を上 2重対角行列 Bに対する特異値分解問題

B = U ∑V ^T

σは Βの特異値

に置き換える。

[0005] 次に、行列 Β^ΤΒを考える。この行列 Β^ΤΒの対角化

Λ =V^TB^TBV

A≡diag(l , λ , ···, λ ) λ ≥λ ≥-"≥λ ≥0

1 2 m 1 2 m

V≡(v , v , ···, v )

λは Β^ΤΒの固有値

j

vは固有値えに対する固有ベクトル

を行う。ここで、一般に、以下のことが成り立つ。（1)B^TBは対称な 3重対角行列である。（2) B^TBの固有値は全て正であり、 Bの特異値 σ (σ ≥ σ ≥···≥ σ ≥0)は、 j 1 2 m

B^TBの固有値え (λ ≥λ ≥-"≥λ ≥0)の正の平方根に等しい。 (3)V =V、す j 1 2 m B なわち、 B^TBの固有ベクトル vは、 Bの右特異ベクトルに等しい。従って、行列 B^TBの対角化が求まると、（2)より Λ =∑²であることから、∑が求まり、さらに（3)より左直交行列 U =BV ∑^_1 = BV∑^_1が求まるので、 Bが特異値分解される。すなわち、 Bの

B B

特異値分解を、 B^TBの対角化の問題に置き換えることができる。この原理は、 m個全ての特異値および特異ベクトルを求める場合だけでなぐ少なくとも 1つの特異値および特異ベクトルを求める場合にも適用され得る。

[0006] 以上のように、一般行列 Aの特異値分解は、対称な 3重対角行列 B^TBの対角化の問題を含む。

発明の開示

発明が解決しょうとする課題

[0007] DBDSQR/レーチンでは、計算量が非常に多いため、特に大規模問題では、高速に特異値分解を実行することが困難であった。一方、 DLASQルーチンは、高速かつ高精度に特異値を求めることができる力 DSTEGRルーチンは、場合によっては特異ベクトルを低精度で計算するため、 DSTEGR/レーチンでは、常に高精度に特異値分解を実行することは困難であった。

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

課題を解決するための手段

[0009] 本発明では、ミウラ型逆変換、 sdLVvs変換、 rdLVvs変換およびミウラ型変換による Twisted (ッイスティッド)分解を実行することにより、行列 B^TBを対角化する。ここで、行列 Bを、行列 Aと特異値が同じ上 2重対角行列とする。行列 Aから行列 Bへの変換は、例えばハウスホルダー法で実行できる。本発明による行列 B^TBの対角化では、固有値および固有ベクトルを同時に求める DBDSQR/レーチンとは異なり、 DSTEG Rルーチンのように、まず固有値を計算し、次に計算された固有値を用いて固有べクトルを求める。

[0010] 本発明により、コンピュータを用いて任意の行列 Aに対して特異値分解を実行する方法であって、行列 Aに対して上 2重対角化を行い、行列 Aの上 2重対角行列 Bを求めるステップと、行列 Aの特異値として行列 Bの少なくとも 1つの特異値 σを求めるステツプと、 σに対する行列 Αの特異ベクトルを求めるステップとを包含し、該行列 Aの特異ベクトルを求めるステップは、ミウラ型逆変換と、 sdLVvs変換と、 rdLVvs変換と、ミウラ型変換とを用いて行列 Β^ΤΒ— σ ²Ιに対してッイスティッド分解を実行することにより、行列 B^TBを対角化するステップを含み、 Iは単位行列である、方法が提供され、これにより上記目的が達成される。

[0011] 前記行列 Aの特異ベクトルを求めるステップは、前記行列 B^TBを対角化するステツプの後に、逆反復法を実行するステップを含んでもょヽ。

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

[0013] 前記行列 Bの特異値 σを求めるステップは、 dLVsルーチンを実行するステップを含んでもよい。

[0014] 前記行列 Bの特異値 σを求めるステップは、要求される計算時間および精度に応じて、 dLVsルーチンもしくは DLASQルーチンのどちらを実行するのかを判定するステツプを含んでもよい。

[0015] 前記行列 Aに対して上 2重対角化を行い、上 2重対角行列 Bを求めるステップは、上 2重対角化法 (たとえばノヽウスホルダー法)を実行するステップを含んでもょ、。

[0016] 本発明により、物体の複数の 2次元画像から 3次元画像を復元する方法であって、 2次元画像; j (j = l, · · · , m、 mは 3以上の整数）における特徴点 i(i= l, · ' · , η、 ηは 2以上の整数)の座標 d (X , y )を抽出するステップと、該特徴点の 2次元座標 d (X

, y )より、該特徴点の 3次元座標 s (X , Y , Z )および 2次元座標から 3次元座標への変換を表す行列 Mを計算するステップとを包含し、該特徴点の 3次元座標 s (X , Y , Z )および 2次元座標から 3次元

座標への変換を表す行列 Mを計算するステップは、行列 Dに対して上 2重対角化を行い、行列 Dの上 2重対角行列 Bを求めるステップであって、行列 Dは、

[0017] [数 5]

と定義される、ステップと、行列 Dの特異値として行列 Bの特異値 σ ( σ ≥ σ ≥· · · j 1 2

≥ σ >0、rは、行列 Dのランクに等しい）を求めるステップと、 σ 、 σ および σ に対 r 1 2 3 する行列 Dの特異ベクトルを求めるステップと、 M = M' Cなる行列 Cに対して、 E = C C^Tを満たす行列 Eを計算するステップであって、 M， =L，（∑，） ^1/2、∑，は σ 、 σ お

1 2 よび σ を対角要素に持ちそれ以外の要素が 0である 3 X 3行列、 L'は σ 、 σ およ

3 1 2 び σ に対応する Dの特異ベクトルを左から順に並べた行列である、ステップと、行列

3

Εから、行列 Cを計算するステップと、行列じから、該 3次元座標 s (X , Υ , Ζ )および該変換を表す行列 Μを計算するステップとを含み、該 σ 、 σ および σ に対する行

1 2 3

列 Dの特異ベクトルを求めるステップは、ミウラ型逆変換と、 sdLVvs変換と、 rdLVvs 変換と、ミウラ型変換とを用いて行列 Β^ΤΒ— σ ²1に対してッイスティッド分解を実行することにより、行列 Β^ΤΒを対角化するステップを含み、 Iは単位行列である、方法が提供され、これにより上記目的が達成される。

本発明により、与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する方法であって、該キーワードに対応する質問ベクトル qを受け取るステップと、該文書に対応する索引語文書行列 Dに対して上 2重対角化を行い、行列 Dの上 2重対角行列 Bを求めるステップと、行列 Dの特異値として行列 Bの特異値 σ ( σ ≥ σ ≥· · ·≥ σ >0、 rは、行列 Dのランクに等しい）を求めるステップと、 k j 1 2 r

<rなる kを選択するステップと、行列 Dの擬似行列 Dを計算するステップであって、 k

行列 Dは、 σ , σ , · · · , σ を対角要素としそれ以外の要素が 0である行列∑、 σ k 1 2 k k

, σ , · · · , σ に対応する特異ベクトルのみを左力順に並べた左右直交行列 U ,

1 2 k k

Vを用いて、 D =U ∑ V ^Tと定義される、ステップと、行列 Dと質問ベクトル qとの類 k k k k k k

似度を計算するステップと、該計算された類似度を基準に、検索結果を出力するステップとを包含し、該行列 Dの左右直交行列 U , Vを求めるステップは、ミウラ型逆変 k k k

換と、 sdLVvs変換と、 rdLVvs変換と、ミウラ型変換とを用いて行列 B¾— σ ²l (j = l , 2, · · · , k)に対してッイスティッド分解を実行することにより、行列 Β^ΤΒを対角化するステップを含み、 Iは単位行列である、方法が提供され、これにより上記目的が達成される。

[0019] 本発明により、任意の行列 Αに対して特異値分解を実行する特異値分解処理をコンピュータに実行させるプログラムであって、該特異値分解処理は、行列 Aに対して上 2重対角化を行い、行列 Aの上 2重対角行列 Bを求めるステップと、行列 Aの特異値として行列 Bの少なくとも 1つの特異値 σを求めるステップと、 σに対する行列 Αの特異ベクトルを求めるステップとを包含し、該行列 Aの特異ベクトルを求めるステップは、ミウラ型逆変換と、 sdLVvs変換と、 rdLVvs変換と、ミウラ型変換とを用いて行列 B^TB— σ ²1に対してッイスティッド分解を実行することにより、行列 Β^ΤΒを対角化するステツ

プを含み、 Iは単位行列である、プログラムが提供され、これにより上記目的が達成される。

[0020] 本発明により、物体の複数の 2次元画像から 3次元画像を復元する画像復元処理をコンピュータに実行させるプログラムであって、該画像復元処理は、 2次元画像； j (j = 1, · · · , m、mは 3以上の整数）における特徴点 i (i= l, · · · , n、 nは 2以上の整数 )の座標 d (X , y )を抽出するステップと、該特徴点の 2次元座標 d (X , y )より、該特徴点の 3次元座標 s (X , Υ , Z )および 2次元座標から 3次元座標への変換を表す行列 Mを計算するステップとを包含し、該特徴点の 3次元座標 _Si (X Y , )および 2次元座標から 3次元座標への変換を表す行列 Μ を計算するステップは、行列 Dに対して上 2重対角化を行い、行列 Dの上 2重対角行列 Βを求めるステップであって、行列 Dは、

[数 6]

3

Εから、行列 Cを計算するステップと、行列じから、該 3次元座標 s (X , Υ , および該変換を表す行列 Mを計算するステップとを含み、該 σ 、 σ および σ に対する行

1 2 3

列 Dの特異ベクトルを求めるステップは、ミウラ型逆変換と、 sdLVvs変換と、 rdLVvs 変換と、ミウラ型変換とを用いて行列 Β^ΤΒ— σ ²1に対してッイスティッド分解を実行することにより、行列 Β^ΤΒを対角化するステップを含み、 Iは単位行列である、プログラムが提供され、これにより上記目的が達成される。

[0022] 本発明により、与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する文書検索処理をコンピュータに実行させるプログラムであって、該文書検索処理は、該キーワードに対応する質問ベクトル qを受け取るステップと、該文書に対応する索引語文書行列 Dに対して上 2重対角化を行い、行列 Dの上 2重対角行列 Bを求めるステップと、行列 Dの特異値として行列 Bの特異値 σ ( σ ≥σ j 1 2

≥· · ·≥ σ >0, rは、行列 Dのランクに等しい）を求めるステップと、 k<rなる kを選択するステップと、行列 Dの擬似行列 Dを計算するステップであって、行列 Dは、 σ , k k 1 σ

2

, · · · , σ を対角要素としそれ以外の要素が 0である行列∑ 、 σ , σ , · · · , σ に k k 1 2 k 対応する特異ベクトルのみを左力順に並べた左右直交行列 U , Vを用いて、 D = k k k

U ∑ V ^Tと定義される、ステップと、行列 Dと質問べ外ル qとの類似度を計算するス k k k k

テツプと、該計算された類似度を基準に、検索結果を出力するステップとを包含し、該行列 Dの左右直交行列 U , Vを求めるステップは、ミウラ型逆変換と、 sdLVvs変 k k k

換と、 rdLVvs変換と、ミウラ型変換とを用いて行列 B^TB— σ

²l (j = l, 2, · · · , k)に対してッイスティッド分解を実行することにより、行列 B^TBを対角化するステップを含み、 Iは単位行列である、プログラムが提供され、これにより上記目的が達成される。

[0023] 本発明により、任意の行列 Aに対して特異値分解を実行する装置であって、行列 A を入力として受け取る手段と、行列 Aに対して上 2重対角化を行い、行列 Aの上 2重対角行列 Bを計算する手段と、行列 Aの特異値として行列 Bの少なくとも 1つの特異値 σを計算する手段と、 σに対する行列 Αの特異ベクトルを求める手段とを備え、該行列 Aの特異ベクトルを求める手段は、ミウラ型逆変換、 sdLVvs変換および rdLVvs 変換、ならびにミウラ型変換を用いて行列 Β^ΤΒ— σ ²Ιに対してッイスティッド分解を実行することにより、行列 Β^ΤΒを対角化する手段を含み、 Iは単位行列である、装置が提供され、これにより上記目的が達成される。

[0024] 本発明により、物体の複数の 2次元画像を 3次元画像に復元する装置であって、 m 枚 (mは 3以上の整数)の 2次元画像を受け取る手段と、 2次元画像； j (j = l, · · · , m) における特徴点 i(i=l, ···, n、 nは 2以上の整数)の座標 d (x , y )を抽出する手段と、該特徴点の 2次元座標 d (X , y )より、該特徴点の 3次元座標 s (X, Υ, Z)および 2次元座標から 3次元座標への変換を表す行列 Mを計算する手段とを備え、該特徴点の 3次元座標 s (X, Υ, Z)および 2次元座標から 3次元座標への変換を表す行列 Mを計算する手段は、行列 Dに対して上 2重対角化を行い、行列 Dの上 2重対角行列 Bを求める手段であって、行列 Dは、

[数 7]

と定義される、手段と、行列 Dの特異値として行列 Bの特異値 σ (σ ≥σ ≥···≥ j 1 2 σ >0、rは、行列 Dのランクに等しい）を求める手段と、 σ 、 σ および σ に対する r 1 2 3 行列 Dの特異ベクトルを求める手段と、 M = M'Cなる行列 Cに対して、 E = CC^Tを満たす行列 Eを計算する手段であって、 M， =L，（∑，） ^1/2、∑，は σ 、 σ および σ を

1 2 3 対角要素に持ちそれ以外の要素が 0である 3X3行列、 L'は σ 、 σ および σ に対

1 2 3 応する Dの特異ベクトルを左力も順に並べた行列である、手段と、行列 Εから、行列 C を計算する手段と、行列じから、該 3次元座標 s (X, Υ, Ζ)および該変換を表す行列 Μを計算する手段とを含み、該 σ 、 σ および σ に対する行列 Dの特異ベクトルを求める手段は、ミウラ型逆変換と、 sdLVvs変換と、 rdLVvs変換と、ミウラ型変換とを用いて行列 B^TB— σ ²1に対してッイスティッド分解を実行することにより、行列 Β^ΤΒ を対角化する手段を含み、 Iは単位行列である、装置が提供され、これにより上記目的が達成される。

[0026] 本発明により、与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する装置であって、該キーワードに対応する質問ベクトル qを受け取る手段と、該文書に対応する索引語文書行列 Dに対して上 2重対角化を行い、行列 Dの上 2重対角行列 Bを求める手段と、行列 Dの特異値として行列 Bの特異値 σ ( σ ≥ σ ≥· · ·≥ σ >0、 rは、行列 Dのランクに等しい）を求める手段と、 k<rなる k

1 2 r

を選択する手段と、行列 Dの擬似行列 Dを計算する手段であって、行列 Dは、 σ , k k 1 σ , · · · , σ を対角要素としそれ以外の要素が 0である行列∑、 σ , σ , · · · , σ

2 k k 1 2 k に対応する特異ベクトルのみを左力順に並べた左右直交行列 U , Vを用いて、 D k k k

=U ∑ V ^Tと定義される、手段と、行列 Dと質問べ外ル qとの類似度を計算する手 k k k k

段と、該計算された類似度を基準に、検索結果を出力する手段とを包含し、該行列 D の左右直交行列 U , Vを求める手段は、ミウラ型逆変換と、 sdLVvs変換と、 rdLVv k k k

s変換と、ミウラ型変換とを用いて行列 B^TB— σ ²l (j = l, 2, · · · , k)に対してツイステイツド分解を実行することにより、行列 B^TBを対角化する手段を含み、 Iは単位行列である、装置が提供され、これにより上記目的が達成される。

発明の効果

[0027] 本発明によれば、ミウラ型逆変換、 sdLVvs変換、 rdLVvs変換およびミウラ型変換によるッイスティッド分解を実行することにより、行列 B^TBが対角化される。それにより、 DSTEGR/レーチンと比較した場合、精度の点ではるかに優れ、実用化に耐え得る固有値計算を実現し、高精度の特異値分解が可能になる。また、逆反復法、 Gram Schmidt (グラムシュミット）法、本発明にお、て提供する dLVsルーチンのうちの少なくとも 1つを併用することにより、 DBDSQRと比較して精度の点では同等になることもあるが、速度の点では圧倒的に高速である特異値分解を実現することができる。さらに、 DBDSQRルーチンと比較すると、固有値および固有ベクトルを同時の求めるの力、または、計算された固有値を用いて固有ベクトルを求めるのかの計算フローの違いにより、以下で説明する 2次元から 3次元への画像復元問題および文書検索問題等においては、計算時間をさらに短縮することができる。

図面の簡単な説明

[図 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の特異値 σ ( σ ≥ σ ≥· · ·≥ σ ≥0、 mは、行列 Β^ΤΒの行列サイズ j 1 2 m

に等しい)を計算する。ここで、本発明により提供される特異値計算ルーチン dLVsおよび LAPACKにより提供される DLASQルーチンを利用し得る。なお、行列 Bの特異値は、行列 Aの特異値に等しぐかつ、行列 B^TBの固有値の正の平方根え ^1/2に等しい。 dLVsルーチンは、 Bの特異値 σを高精度で計算する。 dLVsルーチンの詳細は、後述される。 dLVsで求められる特異値は、従来の LAPCKルーチン（DBDS QR、 DLASQ等）と比較して、最も高精度である。速度に関しては、ほとんどの場合 DLASQルーチンには劣るが、信頼性の高!、DBDSQR/レーチンの約 2分の 1の時間で計算できる (CPUが Itanium2の場合は dLVsルーチンが最速)。そこで、ステツプ 104において、高精度性が必要とされるかどうかを判定する。最も高精度な特異値分解を行いたい場合には、ステップ 106において dLVsを、高精度性および高速性を両立したい場合には、ステップ 108において従来の DLASQルーチンを、特異値計算ルーチンとして用いるのがよい。さらに、 DLASQルーチンが有限回の演算で停止するかは定かでないが、 dLVsルーチンは収束することが保証されているので、信頼性を重視するならば、 dLVsルーチンを用いるのがよい。

[0033] 次に、ステップ 110において、計算された Bの特異値 σ (すなわち、 Αの特異値に j

等しい）用いて、 1. 1において説明される sdLVvs変換および rdLVvs変換によるツイスティッド分解により、 B^TBの対角化を行う。本発明による行列 B^TBの対角ィヒルーチンの詳細は、 1. 1で説明する。さらに、ステップ 112において、以下の 1. 1にて詳細が示される簡略化連立 1次方程式

[0034] [数 8]

≡ ( ， 1,0,- ，0)^τ を求めると、行列 Βの特異ベクトル Vが得られる（すなわち、 Β^ΤΒの対角化行列 Vが得られる）。すなわち、 Β^ΤΒが対角化される。ステップ 114において、∑≡diag ( a

1， σ 2

, · · · , σ ) (ただし、 σ ≥σ ≥· · ·≥σ ≥0)とすると、行列 Βの右直交行列 V = m 1 2 m B

Vであるので、行列 Bの左直交行列 U =BV ∑^_1 = BV∑^{_ 1}から Uを得る。すなわ

B B B

ち、行列 Bが特異値分解される。ここでは、 m個全ての特異値および特異ベクトルを計算したが、この原理は、少なくとも 1つの特異値および特異ベクトルを計算する場合にも適用され得る。

[0035] このとき得られる左右直交行列、すなわち特異ベクトル計算の結果は、 DBDSQR ルーチンによって得られる結果と比較して精度の点で劣る。そこで、ステップ 116において、逆反復法を用いることにより、 DBDSQRルーチンと同程度の精度を得ることができる。さらなる高精度性が必要とされる場合は、ステップ 120において、グラムシュミット法による再直交化を行えば、 DBDSQRよりも高精度での計算を実現し得る。 I SVDルーチンでは、速度に関しては、表 2から、 DBDSQRと比較して大幅に短縮され得ることが理解される。また、表 1から理解されるように、 I— SVDルーチンは、 D STEGRによる極度の精度悪ィヒがみられる場合でも、高精度に特異値分解することができる。なお、表 1、 2の誤差および計算時間に関するデータは、逆反復法を用いているが、グラムシュミット法は用いていない場合の結果である。また、ステップ 116および 120を実行する順序は逆であってもよいし、あるいは、両ステップは、省略してもよい。 [0036] このように、図 1に示される各ステップの機能は、ソフトウェア（例えば、プログラム）によって実現される。しかし、本発明は、これに限定されない。図 1に示される各ステツプの機能をノ、一ドウエア（例えば、回路、ボード、半導体チップ）によって実現してもよ V、し、ソフトウェアとハードウェアとの組み合わせによって実現してもよ、。

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

[0038] 表 1は、本発明による行列 B^TB対角ィ匕ルーチンを実行した後に、逆反復法を実行したときの従来技術および本発明の誤差を比較した表である。

[0039] [表 1]

I -SVD DBDSQR DSTEGR

(dLVsを含む）

II B-U∑VT L 1.69E-13 1.31E-13 4.40E+04

ΙΙ υτυ-HL· 2.09E-13 1.26E-13 ―

Ιΐντν-HL· 4.16E-13 1.21E-13 2.04E+04

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

[表 2]

表 2 :計算時間 (秒）

I -SVD

行列サイズ DBDSQR DSTEGR

(dLVsを含む）

1000 1.12 56.85 0.20

2000 4.88 1175.34 0.84

3000 11.84 4879.39 1.78

また、 1.3に、本発明の方法において、 DSTEGRルーチンと比較して、誤差が低減される理由を説明する。

[0041] 1. 1.本発明による行列 B¾の対角化ルーチン

DBDSQRは B^TBの固有値えとともに対応する固有ベクトル vを求めることができるので、他のルーチンを併用しなくても B^TBの対角化 (Bの特異値分解）が得られる。し力しながら、 DLASQおよび dLVsは上述したように固有ベクトル Vを計算する機能はない。よって、別の固有ベクトル計算ルーチンが必要となる。ただし、 B^TBの固有値え (Bの特異値 σ )が求まっているものとする。

まず、 DSTEGRによる固有ベクトルの計算方法を示す。 Parlettらは以下の（1)、 ( 2)、 (3)の手順で固有ベクトル計算ができることを示して、る。

(1) stationary qd with shift変換、 stqds変換)およひ reverse— time progre ssive qd with shift変換 (rpqds変換）を用いた（B(°)) ^TB(°)— λ Iの Cholesky (コ

j

レスキー）分解 (B^(±1)) ^TB^(±1) = (B⁽⁰⁾) ^TB⁽⁰⁾— λ Iを求める。ただし、 B⁽⁰⁾ =Bであり、

j

B⁽⁺¹⁾および B(^_1)はそれぞれ上 2重対角および下 2重対角行列とする。すなわち

[0043] [数 9]

Γ

となる。すなわち、 {q ^{(± 1)}, e ^(±1) } ¾

k k

(2)コレスキー分解より（B(°)) ^TB(⁰) - λ Iのッイスティッド分解

[0044] [数 10]

を形成する。ただし

[0045] [数 11]

D_fc 三 diag^¹),..., ^ , l,g ¹ )，·· 一¹))，

¾ , ひ ί— _≡

(e ( ) +q. ^{0)-λ )が最小になる kを kとする。

とし、 k-l "k

(3)簡約化された連立 1次方程式

[0046] [数 12]

vJ = ^e^' 三 ( ，^ 1，0，一，0)丁

(t„ 1)Θ

を解く。

[0047] 図 2は、 Parlettらの手法で最も核となる（1)の処理手順を示す。 (1)で求めたいのは、 {q (ⁿ), e (ⁿ)}， n=±lである力そのための変換として stationary qd with s k k

hift(stqds)変換

(1) , (1) (o) (o) . _Ί

q 十 e =q 十 e —λ， ¾:=1， 2， ···， πι，

k k-l k k-l i

(o) (o) .

q e ^:q e ， k= -， 2， …， m- k k k k

(0)_ (1)

=0, e

および re verse— time progressive qd with shift (rpqds)変換

(—1)1 (- -1) (0)_

[ +e =q 十 e λ , k=l, 2, ··■

k k k k- - 1 j

(-D (- 1) (0) (0)

% [ e k=l_; , 2， ···， m— 1,

k+l k ^=qk ^ek ，

(- e ： =0, e "ョ 0

0 c n 1

が採用されている。

[0048] 固有値； lが既知ならば反復的な計算は不要なため計算量力 SQRアルゴリズムに比ベて圧倒的に少なく抑えられる力常に数値安定かつ精度のよい方法とはいえない。 stqds変換と rpqds変換ともに減算による桁落ちが発生する可能性がある。例えば、 stqds変換において q ⁽⁰⁾ + e ⁽⁰⁾ - λ〜e ("ならば、 q ⁽¹) (=q ⁽⁰⁾ + e ⁽⁰⁾ - λ

k k-1 j k-1 k k k-1 j

-e を求める際、倍精度計算であっても q ⁽¹⁾の有効桁はわず力 1桁になること k-1 k

もある。そ IIIの場合、 q ⁽⁰⁾e ⁽⁰⁾/q ⁽¹⁾を計算すると誤差が生じる、つまり、 e ⁽¹⁾が精度よ k k k k

く計算できないことになる。また、 q ⁽¹⁾を求めるのに e ⁽¹⁾が必要、 e ⁽¹⁾を求めるのに

\ k+1 k k

q ⁽¹⁾が必要と逐次的な計算が要求されるので、 1個所で発生した桁落ちによる誤差 k

が波及し、さらなる誤差増大の可能性も秘めている。その結果、理論上は q ⁽¹⁾≠0で k あるが誤差蓄積により q ⁽¹)=0となり、 q (^G)e ⁽⁰⁾/q ⁽¹)の計算においてオーバーフロ k k k k

一が起こるとヽつた数値不安定な状況も想定される。

B⁽⁰⁾ = Bの成分 {b ⁽⁰⁾, b (⁰⁾}が与えられる、すなわち {q ⁽⁰⁾, e ⁽⁰⁾}が与えられると

2k- 1 2k k k

、 λおよび e ⁽¹⁾がー意的に決まるので、この状況は避けることはできない。 rpqds j k-1

変換も同様の性質を持っため、実用的なレベルにまで達したとはいいがたい。 LAP ACKにお!/、て FORTRANルーチン DSTEGRとして改良版が公開されて!、るものの欠点は完全に解決されては、な、。

[0049] 次に、本発明による対角化のルーチンを説明する。本発明によるルーチンでは、 st qds変換および rpqds変換にかわる新たな変換を採用している。本発明では、行列 B TBを対角化するために、 sdLVvs変換および rdLVvs変換を用いた行列 Β^ΤΒ— σ ²Ι に対するッイスティッド分解を行う。正確には、上記（1)を以下の（la)、 (lb), (lc) の 3つの工程にわけることで数値安定なコレスキー分解を実現している。

[0050] 図 3〜図 8は、本発明によるミウラ型逆変換、 sdLVvs変換、 rdLVvs変換およびミゥラ型変換を実行する処理手順を示す。

(la)ミウラ型逆変換

[0051] [数 13]

を行う。

(0) (0)

[0052] u ⁽⁰⁾= 7？ k=l ·, m,

2k- 1 k

(0) (0)

u =e /t ⁽⁰⁾, k=l, ·, m— 1,

2k k k

/ (0) (0)) _ -, (o) _ / _s (0)

( +u 1, =1/ 0

2k- 2

ただし、 δ ^(C))は任意に選ぶことができる。 (lb) stationary不等間隔差分 Lotka— Volterra (sdLVvs)変換

u ⁽¹⁾=u ⁽⁰⁾>f (l+ 6 ⁽⁰⁾u ⁽⁰⁾)/(l+ 6⁽¹⁾u

k k k-1 k一 1

k=l, 2, ···, 2m— 1,

によって

[0053] [数 14]

4。) - 4¹)

reverse— time不等間隔差分 Lotka— Volterra (rdLVvs)変換

(-1) (0) ,_Λ (0) (0ヽ / (_Λ j. (-1)

u =u * (1+ δ u )/ (1+ δ u

k k k-1

k+1 (

k=l, 2, ···, 2m— 1,

によって

[0054] [数 15] ) - 4 - を実行する。ただし、 u ⁽ⁿ⁾≡0, u (⁰⁾≡0とし、 L = 1Z δ ⁽⁰⁾— 1/ δ ^(η)， η=士

k 2m j 1を満たすように δ ^ω, η=±1を設定する。

(lc)ミウラ型変換

[0055] [数 16]

q ^ω = * + _δ (ⁿ)u ) * + _δ (n)_u ^ω)，

k 2k- 2 2k- 1

k=l, 2, ···, m,

(n) _s (n) . (n) . (n)

e = o 水 u 水 u

k 2k- 1 2k

k=l, 2, ···， m-1,

を行う。

[0056] qd型変換では見られない離散 Lotka— Volterra型変換の大きな特徴は、任意パラメータを持つことである。すなわち、え =1/ δ ()— 1Z δ ^(±1)を満たす範囲で δ ^(η) の値を任意に設定できる。 δ ^(η)を変動させると補助変数 u ^(η)の値も変化するが、桁 k

落ちによる誤差や数値不安定が発生するかどうかは事前に判定できる。この判定は、 if文によって実装されてもよい。この場合は、 δ ^(η)を再設定後に再度計算すればよい。また、 u ^(±1)が求まればミウラ変換によって q ^(±1)および e ^(±1)が独立に計算されるの

k k k

で、誤差が伝播しないという性質を持つ。なお、ミウラ型逆変換をミウラ変換、ミウラ型変換を逆ミウラ変換、 sdLVvs変換を stLVv変換、 rdLVv変換を rLVv変換と呼び変えても問題はない。

[0057] 入力および出力の変数の対応は、 Parlettの手法と同じである。メモリ消費を抑えるために、補助変数 u ⁽ⁿ⁾のための配列は必ずしも用意する必要はない。一方、 1 + δ ⁽

k

のためのメモリ領域を確保し、 (la)〜（lc)のステップにまたがつてこの値を利用することで、メモリ消費を抑え、計算量を低減することができる。これにより、誤差も低減される。このように、 1. 3で説明する方法を利用して誤差を低減すること等、さらなる工夫を行い得る。

[0058] 上記の方法で B^TB— λ Iのコレスキー分解を求めた後は、 Parlettらと同じ方法を用いて特異ベクトル Vを計算する。そうすると、 B^TBの固有値分解 A =V ^TB^TBVが得

j B B られる。ただし Vの列ベクトル力^である。 U =BV ∑^_1より Uを求めると、 Bの特異

B

値分解 B = U ∑V ¹も得られる。

[0059] 1· 2.高精度特虽値 dLVsルーチン

本発明による dLVsルーチンは、基本的な枠組みでは DLASQルーチンと同じであり、 B^TBの固有値（あるいは Bの特異値）のみが得られ B^TBの固有ベクトル（あるいは B の特異ベクトル）は求まらない。ただ、ルーチンの核となる部分に使用される計算法は DLASQと異なる。

[0060] まず、 DLASQルーチンによる固有値の計算方法を説明する。 qd法は固有べタトルが計算できな!/、が、 1ステップは上 2重対角行列の要素を更新するのみなので高速に特異値を求めることができる。かつ、計算量が少ないため丸め誤差の蓄積が抑えられ高精度計算が期待できる。ここで、 Y⁽ⁿ⁾を対称な 3重対角行列

[0061] [数 17]

とする。また、変換 Y⁽ⁿ⁾→Y⁽ⁿ⁺¹⁾が qd法の漸化式

q =q 十 e — e ， k= 1, 2, · · m,

e =e 氺 q / q ， k= l， 2, · · m—

e =0, e =0, n— 0, 1，

によって実現されるとする。そのとき、適切な R⁽ⁿ⁾が存在して Y⁽ⁿ⁺¹⁾ =R⁽ⁿ⁾Y⁽ⁿ⁾ (R⁽ⁿ⁾) _ 1となることを Rutishauserは発見している。これは、 Y⁽ⁿ⁺¹⁾と Y⁽ⁿ⁾の固有値が同じであることを意味し、すなわち、上記の漸化式によって固有値保存変形がなされる。この変形を繰り返せば非対角成分が 0に近づき、 Υ^ωが対角化されることも証明されている。よって Υ(°) = Β^ΤΒとすると lim

Y⁽ⁿ⁾ =diag ( l , λ , · · · , λ )となる。ただし、 λ は Β^ΤΒの固有値である。さらに、 λ の正の平方根をとると Βの特異値 σ が得られる。

[0062] LAPCAKの特異値計算ルーチン DLASQでは、 qd法の differential型が採用されている。これは dqd法と呼ばれる。 1ステップは

[0063] [数 18] d - for k：= l, m— 1 ei"⁺¹⁾：= ef * ( "）

= (Οί""))

end for で与えられ、変換 ^Y^"が実現される。 differential型は減算を含まないため、 qd法では起こりうる桁落ち誤差の心配もない。さらに高速版として原点シフト付き dq d (dqds)法も組み込まれて!/、る。高速ィ匕に伴う計算量の減少で丸め誤差の発生も抑えることができる。ただし、収束するかどうかの保証はない。原点シフト版の 1ステップは

[0064] [数 19] j (n)

a := q\ ⁷― s

end for

+¹) := d

で与えられる。 sは原点シフト量であり、 Gersgorin型境界条件により見積もられた Y⁽ⁿ )の最小固有値とする。よって、 DLASQルーチンのメインループの核となる部分では最小固有値の見積もりを行い、 s > ε ( ε：非常に小さいとみなせる正の数)ならば dq ds法が、そうでなければ dqd法が用いられる。メインループの他の部分 (行列の分割. 減次）はほとんど DBDSQR/レーチンと変わらない。すなわち、 DBDSQR/レーチンに含まれる QR法を dqd法に、原点シフト QR法を dqds法に置き換えたものが DLAS Qルーチンとなる。（QR法と dqd法、原点シフト付き QR法と dqds法で扱う変数が異なるので少しだけ判定式が異なる部分もあるが本質的に同じである。 )

次に本発明による dLVsルーチンを説明する。 dLVsルーチンでは dqd法のかわりに dLV法を、 dqds法のかわりに dLVs法を採用する。ここで、 dLVsルーチンは DLA SQルーチンとは違って収束することが保証されている。その 1ステップは以下で与えられる。

[数 20]

for = 2, 2m— 1

u ： = + δ^{η)* ^η ); νΆ = * (i + <5 ⁽ⁿ⁾ * ⁿ⁾);

end for

(n) 一 M

u2m-l '― ^u2m~l

(原点シフト sの決定）

if (s > e) w^^{( +1)}： = v[ⁿ⁾― s;

for — 2, m end for

end if

else

for k = 1, 2m— 1 end for

end else

DLASQルーチン同様、メインループの核の部分では最小特異値の見積もりによつて原点シフト量 sを決定し、 dLV法を使うか dLVs法を使うかを決定する。ただし、補助変数 V ^ωを求めたあとで sを決定する点は DLASQルーチンと異なる。

k

[0066] 1. 3.本明の方法による差の低減

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

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

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

[0069] 本発明による計算方法では、上述の桁落ちによる誤差を含んで、るかどうかは減算によって得られた値が小さいかどうかで判断できる。 Parlettの手法の場合、図 2の D O文中の ql [k]、 q2 [k]をチェックすればよい。ところが、 Parlettの手法は桁落ち誤差を含むことが分力たとしても、それを回避することはできない。なぜならば、初期値として {qO[k]， eO[k]}が与えられると、 λは一意的に決定され、 {ql[k]， el[k] }および {q2[k]， e2[k]}も一意的に導出されるためである。すなわち、任意パラメータを持たない自由度のない計算法であるためである。

[0070] それに対して本発明による方法は自由に設定できるパラメータ 7? 1 ( = ΐΖ δ ( ))をもっため、補助変数 u ^(η)の値を様々に変化させることができる（図 9参照)。すなわち

k

、様々な経路で {ql[k]， el[k]}および {q2[k]， e2[k]}を計算することができる。よつて、桁落ちが発生する場合も回避できる。図 6〜図 8の条件判定によって桁落ちの影響をチェックし、減算によって得られた値の絶対値力、さな数 εより大きいという条件が満たされなければ、パラメータ r? 1の設定に戻るというものである。この処理は、条件が満たされるまで繰り返される。なお、精度よりも高速性を重視する場合は、数回条件が満たさなければ (ql[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 = l, · · ·, m、 mは 3以上の整数）から特徴点 i(i=l, ···, n、nは 2以上の整数)の座標 (x^j, y^j)を抽出する。取り扱う 2次元画像は、弱中心射影画像であることが好ましい。このとき

[0079] [数 21]

が成り立つ。ここで、 sは物体のスケールに相対する j番目の画像のスケール、 r , r

j 1, j 2 はそれぞれ物体座標系に相対する j番目のカメラ座標系の回転行列の 1番目と 2番

， j

目の行ベクトル、 [X, Y, Z.]^T«i番目の点の 3次元座標である。物体のスケールは 1 番目の画像のスケールと同じにし（s =1)、物体の座標系の姿勢は 1番目の画像のカメラ座標系と同じにする (r =[1, 0, 0]^T, r =[0, 1, 0]^T)。 m枚全ての画像における n個全ての点の座標を行列 Dの要素として並べると、

D = MS

が得られる。

ただし、

[0080] [数 22]

である。 Mと Sの形から分かるように、 Dのランクは 3である。いま、ステップ 1102より D が与えられている。以下、回転に関するデータ Mおよび形状 Sを求める。

そこで、行列 Dの特異値分解

D = U∑V^T

を考える。ここで、∑は特異値を大小順に対角線上に並べたもので、 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，は II D-D' IIを最小にするランク 3の行列である。

[0084] 次に、 D'から Mおよび Sを求めたいが、その組合せは唯一ではない。なぜなら、任意の正則行列じが

D， = (M，C) (C^{_ 1}S，）

を満たす力もである。そこで、上式における Cを M = M' Cを満たすように決める。ま下記の式を満たす。

[0085] [数 23]

E = CC^Tとすると、上式力 Eの 6つの要素に関する 2m+ 1個の線形方程式が得られる。 m≥ 3であるので、 Eの要素を一意に決めることができる。ステップ 1110において、行列 Eを求める。

[0086] 次に、ステップ 1112において、 Eから Cを求める。 Cの自由度（9)は Eの自由度（6) より多い。そこで、条件 = [1, 0, 0]^T, r = [0, 1, 0 を加えれば、 Cを決めること lj ¾

ができる。このとき 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.本発明による特異値分解方法を刺用した文書枪索方法

文書中からその文書の内容に関連する索引語を抽出し、索引語の重みを計算する処理の後、ベクトル空間モデルでは、この索引語の重みを要素とするベクトルで文書を表現する。いま、検索対象となる文書を D , D , · · · , Dとし、これら文書集合全体

1 2 n

を通して全部で m個の索引語 w , w , · · · , w があるとする。このとき、文書 Dは、次

1 2 m j のようなベクトルで表現されることになる。これを文書ベクトルと呼ぶ。

[0090] [数 24]

ここで、 dは索引語 wの文書 Dにおける重みである。また、文書集合全体は、次のような m X n行列 Dによって表現することができる。

[数 25]

行列 Dを索引語文書行列と呼ぶ。索引語文書行列の各列は文書に関する情報を表している文書ベクトルである力同様に、索引語文書行列の各行は索引語に関する情報を表しているベクトルであり、これを索引語ベクトルと呼ぶ。検索質問も、文書と同様に、索引語の重みを要素とするベ外ルで表現することができる。検索質問文に含まれる索引語 wの重みを qとすると、検索質問ベクトル qは次のように表されることになる。 [0092] [数 26]

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

[0093] [数 27] 尺度

内積

dJ 9 = > , difli なお、ベクトルの長さが 1に正規ィ匕 (コサイン正規化）されている場合には、コサイン尺度と内積とは一致する。

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

[0095] ステップ 1202において、質問ベクトル qを受け取る。次に、ステップ 1204では、行列 Dの特異値分解を計算するために、行列 Dを上 2重対角化して上 2重対角化行列 Bを得る。ステップ 1206では、行列 B、すなわち行列 Dの特異値を求める。ステップ 1 206では、本発明による dLVsルーチン、 DLASQルーチン、その他の特異値計算ルーチン、あるいはその組み合わせを利用し得る。

[0096] 次に、 Dの近似行列 Dを使った検索を考える。ベクトル空間モデルでは、検索質問

k

ベクトル qと索引語文書行列 D中の各文書ベクトル dの間の類似度を計算することにより検索を行うが、ここでは Dの代わりに Dを使う。ベクトル空間モデルでは、文書べ

k タトルの次元数は索引語の総数と等しい。したがって、検索対象となる文書の数が増えるに従い、文書ベクトルの次元数も増加する傾向にある。しかし、次元数が増加してくると、コンピュータのメモリによる制限や検索時間の増大などの問題が生じてくるばかりでなぐ文書中に含まれる不必要な索引語がノイズ的な影響を及ぼし、検索精度を低下させてしまうという現象も起こってくる。潜在的意味インデキシング (latent semantic indexing ; LSI)は、高次元の空間にある文書ベクトルを低次元の空間へと射影することにより、検索精度の改善を図る技術である。高次元の空間では別々に扱われていた索引語力低次元の空間では相互に関連を持ったものとして扱われる可能性もあるため、索引語の持つ意味や概念に基づく検索を行うことができる。たとえば、通常のベクトル空間モデルでは" car"という索引語と" automobile"という索引語はまったく別物であり、一方の索引語による質問ではもう片方の索引語を含んだ文書を検索することができない。しかし、低次元の空間ではこれらの意味的に関連した索引語は 1つの次元に縮退することが期待できるため、 "car"という検索質問によつて" car"を含む文書ば力りでなぐ' automobile"を含む文書をも検索することが可能となる。潜在的意味インデキシングでは、特異債分解により高次元ベクトルの次元削減を行うが、これは基本的に多変量解析における主成分分析と等価である。

[0097] 行列ステップ 1208では、 kく rなる kの値を選択する。 kの値は、予め与えられて!/ヽてもよいし、計算ごとに選択可能であってもよい。

[0098] ステップ 1210において、ステップ 1206で計算した特異値のうち、大きい順に 1番目力も k番目の k個の特異値に対して、 Dの特異ベクトルを計算する。すなわち、 D =U∑V ^T

k k k

なる Uおよび Vを計算する。ここで、 Uは、最初の k個の左特異ベクトルのみから構 k k k

成される m X k行列であり、 Vは、最初の k個の右特異ベクトルのみカゝら構成される n k

X k行列であり、∑は、最初の k個の特異値のみ力構成される k X k対角行列である。ステップ 1210は、図 1のステップ 110〜120を利用し得る。

[0099] 次に、ステップ 1212において、行列 Dと質問ベクトル qとの類似度を計算する。い

k

ま、ベクトル eを n次元の単位ベクトルとすると、 Dの j番目の文書ベクトルは D eで表

i k k j すことができる。文書ベクトル D eと検索質問べ外ル qとの間の類似度計算は、

k j [0100] [数 28]

(D_ke_f) - q (D_ke_f)^Tq

cos(D_ke q) =

11 ll< ll ΙΙ^Α^/Ι

(U_k∑_kV_k ^Tej)^Tq e]V_k∑_k

IIひ ^ II l \\∑_kV e_;\

(∑_kV ejy(U q) k ii ιι ιι

としてもよいが、別の定義を用いてもよい。上式では、 Dを U , ∑ , Vから再構成す k k k k

る必要はなく特異値分解の結果から、直接、類似度を計算できることを示している。数 29の中に現われる∑ V^Teは、

k k j

∑ V ^Te =U ^TD e

k k j k k j

と書き直すことができる。この式の右辺は、近似行列 D

kにおける j番目の文書ベクトルの基底 Uのもとでの座標（文書の k次元表現）を表している。同様に、数 29中の U ^Tq k k は、検索質問ベクトル qの基底 Uのもとでの座標 (検索質問の k次元表現)である。

k

[0101] ステップ 1214において、ステップ 1212において計算された類似度を基準に、検索結果を出力する。

[0102] DBDSQR/レーチンでは r個のベクトルを求めた後、 k個のベクトル（r—k個は不要）を採用するのに対して、この実施形態では、 k個のベクトルだけ求めればよい。すなわち、計算コストを削減することができる。

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

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

Claims

請求の範囲

[1] コンピュータを用いて任意の行列 Aに対して特異値分解を実行する方法であって、行列 Aに対して上 2重対角化を行、、行列 Aの上 2重対角行列 Bを求めるステップと、

行列 Aの特異値として行列 Bの少なくとも 1つの特異値 σを求めるステップと、 σに対する行列 Αの特異ベクトルを求めるステップと

を包含し、

該行列 Aの特異ベクトルを求めるステップは、ミウラ型逆変換と、 sdLVvs変換と、 rd LVvs変換と、ミウラ型変換とを用いて行列 Β^ΤΒ— σ ²Ιに対してッイスティッド分解を実行することにより、行列 Β^ΤΒを対角化するステップを含み、 Iは単位行列である、方法

[2] 前記行列 Αの特異ベクトルを求めるステップは、前記行列 B^TBを対角化するステツプの後に、逆反復法を実行するステップを含む、請求項 1に記載の方法。

[3] 前記行列 Aの特異ベクトルを求めるステップは、前記逆反復法を実行するステップの後に、グラムシュミット法を実行するステップを含む、請求項 2に記載の方法。

[4] 前記行列 Bの特異値 σを求めるステップは、 dLVsルーチンを実行するステップを含む、請求項 1に記載の方法。

[5] 前記行列 Bの特異値 σを求めるステップは、要求される計算時間および精度に応じて、 dLVsルーチンもしくは DLASQルーチンのどちらを実行するのかを判定するステツプを含む、請求項 1に記載の方法。

[6] 前記行列 Aに対して上 2重対角化を行い、上 2重対角行列 Bを求めるステップは、ノ、ウスホルダー法を実行するステップを含む、請求項 1に記載の方法。

[7] 物体の複数の 2次元画像から 3次元画像を復元する方法であって、

2次元画像; j (j = l, · · · , m、mは 3以上の整数）における特徴点 i (i= l, · ' · , η、η は 2以上の整数)の座標 d (X , y )を抽出するステップと、

該特徴点の 2次元座標 d (X , y )より、該特徴点の 3次元座標 s (X , Υ , Z )および

2次元座標から 3次元座標への変換を表す行列 Mを計算するステップと

を包含し、該特徴点の 3次元座標 _Si ( ， Y, および 2次元座標から 3次元座標への変換を表す行列 Mを計算するステップは、

行列 Dに対して上 2重対角化を行、、行列 Dの上 2重対角行列 Bを求めるステップであって、行列 Dは、

[数 1]

と定義される、ステップと、

行列 Dの特異値として行列 Bの特異値 σ (σ ≥ σ ≥···≥ σ >0、 rは、行列 Dの j 1 2 r

ランクに等しい）を求めるステップと、

σ 、 σ および σ に対する行列 Dの特異ベクトルを求めるステップと、

1 2 3

M = M， Cなる行列 Cに対して、 E = CC^Tを満たす行列 Eを計算するステップであつて、 M'=L'(∑')^1/2、∑'は σ 、 σ および σ を対角要素に持ちそれ以外の要素

1 2 3

が 0である 3X3行列、 L，は σ 、 σ および σ に対応する Dの特異ベクトルを左から

1 2 3

順に並べた行列である、ステップと、

行列 Eから、行列 Cを計算するステップと、

行列じから、該 3次元座標 s (X, Υ, Z)および該変換を表す行列 Mを計算するステツプと

を含み、

該 σ 、 σ および σ に対する行列 Dの特異ベクトルを求めるステップは、ミウラ型逆

1 2 3

変換と、 sdLVvs変換と、 rdLVvs変換と、ミウラ型変換とを用いて行列 B¾— σ ²1に対してッイスティッド分解を実行することにより、行列 Β^ΤΒを対角化するステップを含み、 Iは単位行列である、方法。

[8] 与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する方法であって、

該キーワードに対応する質問ベクトル qを受け取るステップと、

該文書に対応する索引語文書行列 Dに対して上 2重対角化を行い、行列 Dの上 2 重対角行列 Bを求めるステップと、

行列 Dの特異値として行列 Bの特異値 σ ( σ ≥ σ ≥· · ·≥ σ >0、 rは、行列 Dの j 1 2 r

ランクに等しい）を求めるステップと、

k< rなる kを選択するステップと、

行列 Dの擬似行列 Dを計算するステップであって、行列 Dは、 σ , σ , · · · , σ k k 1 2 k を対角要素としそれ以外の要素が 0である行列∑ 、 σ , σ , · · · , σ に対応する特 k 1 2 k

異ベクトルのみを左から順に並べた左右直交行列 U ， Vを用いて、 D 二 U ∑ V τと k k k k k k 定義される、ステップと、

行列 Dと質問べ外ル qとの類似度を計算するステップと、

k

該計算された類似度を基準に、検索結果を出力するステップと

を包含し、

該行列 Dの左右直交行列 U , Vを求めるステップは、ミウラ型逆変換と、 sdLVvs k k k

変換と、 rdLVvs変換と、ミウラ型変換とを用いて行列 B^TB— σ ²l (j = l, 2, · · · , k) に対してッイスティッド分解を実行することにより、行列 B^TBを対角化するステップを含み、 Iは単位行列である、方法。

[9] 任意の行列 Aに対して特異値分解を実行する特異値分解処理をコンピュータに実行させるプログラムであって、

該特異値分解処理は、行列 Aに対して上 2重対角化を行、、行列 Aの上 2重対角行列 Bを求めるステップと、

を包含し、

該行列 Aの特異ベクトルを求めるステップは、ミウラ型逆変換と、 sdLVvs変換と、 rd LVvs変換と、ミウラ型変換とを用いて行列 Β^ΤΒ— σ²Ιに対してッイスティッド分解を実行することにより、行列 Β^ΤΒを対角化するステップを含み、 Iは単位行列である、プログラム。

物体の複数の 2次元画像から 3次元画像を復元する画像復元処理をコンピュータに実行させるプログラムであって、

該画像復元処理は、

2次元画像; j(j = l, ···, m、mは 3以上の整数）における特徴点 i(i=l, ·'·, η、η は 2以上の整数)の座標 d (X , y )を抽出するステップと、

該特徴点の 2次元座標 d (X , y )より、該特徴点の 3次元座標 s (X, Υ, Z)および

を包含し、

該特徴点の 3次元座標 _Si (X, Υ, Z)および 2次元座標から 3次元座標への変換を表す行列 Mを計算するステップは、

[数 2]

と定義される、ステップと、

ランクに等しい）を求めるステップと、

1 2 3

順に並べた行列である、ステップと、

行列 Eから、行列 Cを計算するステップと、

を含み、

1 2 3

変換と、 sdLVvs変換と、 rdLVvs変換と、ミウラ型変換とを用いて行列 B¾— σ ²1に対してッイスティッド分解を実行することにより、行列 Β^ΤΒを対角化するステップを含み、 Iは単位行列である、プログラム。

[11] 与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する文書検索処理をコンピュータに実行させるプログラムであって、

該文書検索処理は、

ランクに等しい）を求めるステップと、

k< rなる kを選択するステップと、

行列 Dと質問べ外ル qとの類似度を計算するステップと、

k

を包含し、

変換と、 rdLVvs変換と、ミウラ型変換とを用いて行列 B^TB— σ ²l (j = l, 2, · · · , k) に対してッイスティッド分解を実行することにより、行列 B^TBを対角化するステップを含み、 Iは単位行列である、プログラム。

[12] 任意の行列 Aに対して特異値分解を実行する装置であって、

行列 Aを入力として受け取る手段と、

行列 Aに対して上 2重対角化を行、、行列 Aの上 2重対角行列 Bを計算する手段と行列 Aの特異値として行列 Bの少なくとも 1つの特異値 σを計算する手段と、 σに対する行列 Αの特異ベクトルを求める手段と

を備え、該行列 Aの特異ベクトルを求める手段は、ミウラ型逆変換、 sdLVvs変換および rdL Vvs変換、ならびにミウラ型変換を用いて行列 B^TB— σ ²1に対してッイスティッド分解を実行することにより、行列 Β^ΤΒを対角化する手段を含み、 Iは単位行列である、装置物体の複数の 2次元画像を 3次元画像に復元する装置であって、

m枚 (mは 3以上の整数)の 2次元画像を受け取る手段と、

2次元画像; j(j = l, ···, m)における特徴点 i(i=l, ···, n、nは 2以上の整数)の座標 d (X , y )を抽出する手段と、

2次元座標から 3次元座標への変換を表す行列 Mを計算する手段と

を備え、

該特徴点の 3次元座標 _Si (X, Υ, Z)および 2次元座標から 3次元座標への変換を表す行列 Mを計算する手段は、

行列 Dに対して上 2重対角化を行、、行列 Dの上 2重対角行列 Bを求める手段であつて、行列 Dは、

[数 3]

と定義される、手段と、

ランクに等し、）を求める手段と、

σ 、 σ および σ に対する行列 Dの特異ベクトルを求める手段と、

1 2 3

M = M， Cなる行列 Cに対して、 E = CC^Tを満たす行列 Eを計算する手段であって、 M'=L' (∑')^1/2、∑'は σ 、 σ および σ を対角要素に持ちそれ以外の要素が 0

1 2 3

である 3X3行列、 L，は σ 、 σ および σ に対応する Dの特異ベクトルを左から順に

1 2 3

並べた行列である、手段と、

行列 Eから、行列 Cを計算する手段と、

行列じから、該 3次元座標 s (X, Υ, Z)および該変換を表す行列 Mを計算する手段と

を含み、

該 σ 、 σ および σ に対する行列 Dの特異ベクトルを求める手段は、ミウラ型逆変

1 2 3

換と、 sdLVvs変換と、 rdLVvs変換と、ミウラ型変換とを用いて行列 B¾— σ に対 j してツイスティッド分解を実行することにより、行列 B^TBを対角化する手段を含み、 Iは単位行列である、装置。

与えられた文書に含まれる情報であって、与えられたキーワードに関連する情報を検索する装置であって、

該キーワードに対応する質問ベクトル qを受け取る手段と、

該文書に対応する索引語文書行列 Dに対して上 2重対角化を行い、行列 Dの上 2 重対角行列 Bを求める手段と、

ランクに等し、）を求める手段と、

k<rなる kを選択する手段と、

行列 Dの擬似行列 Dを計算する手段であって、行列 Dは、 σ , σ , · · · , σ を対 k k 1 2 k 角要素としそれ以外の要素が 0である行列∑、 σ , σ , · · · , σ に対応する特異べ k 1 2 k

タトルのみを左から順に並べた左右直交行列 U , Vを用いて、 D =U ∑ V ^Tと定義 k k k k k k される、手段と、

行列 Dと質問べ外ル qとの類似度を計算する手段と、

該計算された類似度を基準に、検索結果を出力する手段と

を包含し、

該行列 Dの左右直交行列 U , Vを求める手段は、ミウラ型逆変換と、 sdLVvs変 k k k

換と、 rdLVvs変換と、ミウラ型変換とを用いて行列 B^TB— σ ²I (j = l, 2, · · · , k)に対してッイスティッド分解を実行することにより、行列 B^TBを対角化する手段を含み、 I は単位行列である、装置。