JP5011545B2 - 特異値分解装置、及び特異値分解方法 - Google Patents

Description

本発明は、特異値分解を行う特異値分解装置等に関する。

特異値分解（ＳＶＤ：Ｓｉｎｇｕｌａｒ Ｖａｌｕｅ Ｄｅｃｏｍｐｏｓｉｔｉｏｎ）はデータ処理の中心的な行列演算として、画像処理やデータ検索等の多くの分野に応用されている。
なお、特異値分解の従来例としては、例えば、下記の非特許文献１のものなどが知られている。
Ｍｉｎｇ Ｇｕ ａｎｄ Ｓｔａｎｌｅｙ Ｃ． Ｅｉｓｅｎｓｔａｔ、「Ａ Ｄｉｖｉｄｅ−ａｎｄ−Ｃｏｎｑｕｅｒ Ａｌｇｏｒｉｔｈｍ ｆｏｒ ｔｈｅ Ｂｉｄｉａｇｏｎａｌ ＳＶＤ」、ＳＩＡＭ Ｊｏｕｒｎａｌ ｏｎ Ｍａｔｒｉｘ Ａｎａｌｙｓｉｓ ａｎｄ Ａｐｐｌｉｃａｔｉｏｎｓ、Ｖｏｌ．１６，Ｎ０．１，ｐｐ．７９−９２（１９９５）

近年、これらの応用分野におけるデータ量の増大などに伴って、高速・高精度な特異値分解が求められている。また、特異値分解を並列処理することができる特異値分解装置等の開発が望まれていた。
本発明は、上記状況の下になされたものであり、並列性に優れた高速、高精度な特異値分解を実行可能な特異値分解装置等を提供することを目的とする。

上記目的を達成するため、本発明による特異値分解装置は、２重対角行列Ｂが記憶される対角行列記憶部と、前記２重対角行列Ｂが２個の２重対角行列に分割され、その２重対角行列が２個の２重対角行列に分割される処理が、分割後の各２重対角行列があらかじめ決められた大きさ以下となるまで繰り返され、当該あらかじめ決められた大きさ以下の各２重対角行列に対して特異値分解が行われた結果である、前記各２重対角行列の特異値と、前記各２重対角行列の特異ベクトルからなる左右直交行列の一部の要素である行列要素とが記憶される特異値分解記憶部と、前記２重対角行列Ｂの特異値が記憶される特異値記憶部と、各２重対角行列の特異値と、行列要素とを前記特異値分解記憶部から読み出し、前記特異値と、前記行列要素とから、分割元の２重対角行列の特異値と、分割元の２重対角行列の行列要素とを算出して前記特異値分解記憶部に蓄積し、その分割元の２重対角行列の特異値と、行列要素とを算出する処理を、２重対角行列Ｂの少なくとも１個の特異値を算出するまで繰り返し、前記２重対角行列Ｂの少なくとも１個の特異値を前記特異値記憶部に蓄積する特異値算出部と、前記２重対角行列Ｂの特異ベクトルが記憶される特異ベクトル記憶部と、前記対角行列記憶部から２重対角行列Ｂを読み出し、前記特異値記憶部から前記２重対角行列Ｂの特異値を読み出し、２重対角行列Ｂとその特異値とから、ツイスト分解法を用いて２重対角行列Ｂの少なくとも１個の特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する特異ベクトル算出部と、を備えたものである。

このような構成により、まず特異値を算出し、その特異値に基づいてツイスト分解法を用いて特異ベクトルを算出することによって、高速で高精度な特異値分解を実現することができうる。また、並列性にも優れている。

また、本発明による特異値分解装置では、前記対角行列記憶部には、前記あらかじめ決められた大きさ以下の各２重対角行列も記憶され、前記あらかじめ決められた大きさ以下の各２重対角行列を前記対角行列記憶部から読み出し、前記各２重対角行列に対して特異値分解を行い、前記各２重対角行列の特異値と、前記各２重対角行列の特異ベクトルとを算出し、当該特異値と、当該特異ベクトルからなる左右直交行列の一部の要素である行列要素とを前記特異値分解記憶部に蓄積する特異値分解部をさらに備えてもよい。

このような構成により、特異値を算出する際に、あらかじめ決められた大きさ以下の各２重対角行列に対する特異値分解を、特異値分解装置において行うことができる。

また、本発明による特異値分解装置は、２重対角行列Ｂが記憶される対角行列記憶部と、前記対角行列記憶部から前記２重対角行列Ｂを読み出し、当該２重対角行列Ｂを２個の２重対角行列に分割して前記対角行列記憶部に蓄積し、その２重対角行列を２個の２重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各２重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割部と、前記あらかじめ決められた大きさ以下の各２重対角行列を前記対角行列記憶部から読み出し、前記各２重対角行列に対して特異値分解を行い、前記各２重対角行列の特異値と、前記各２重対角行列の特異ベクトルとを算出する特異値分解部と、前記特異値分解部によって特異値分解された特異値と、特異ベクトルからなる左右直交行列の一部の要素である行列要素とが記憶される特異値分解記憶部と、前記２重対角行列Ｂの特異値が記憶される特異値記憶部と、各２重対角行列の特異値と、行列要素とを前記特異値分解記憶部から読み出し、前記特異値と、前記行列要素とから、分割元の２重対角行列の特異値と、分割元の２重対角行列の行列要素とを算出して前記特異値分解記憶部に蓄積し、その分割元の２重対角行列の特異値と、行列要素とを算出する処理を、２重対角行列Ｂの少なくとも１個の特異値を算出するまで繰り返し、前記２重対角行列Ｂの少なくとも１個の特異値を前記特異値記憶部に蓄積する特異値算出部と、前記２重対角行列Ｂの特異ベクトルが記憶される特異ベクトル記憶部と、前記対角行列記憶部から２重対角行列Ｂを読み出し、前記特異値記憶部から前記２重対角行列Ｂの特異値を読み出し、２重対角行列Ｂとその特異値とから、ツイスト分解法を用いて２重対角行列Ｂの少なくとも１個の特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する特異ベクトル算出部と、を備えたものである。

このような構成により、まず特異値を算出し、その特異値に基づいてツイスト分解法を用いて特異ベクトルを算出することによって、高速で高精度な特異値分解を実現することができうる。また、並列性にも優れている。さらに、全ての特異値や特異ベクトルを算出する必要がない場合には、必要な範囲で特異値や特異ベクトルを算出することができ、処理負荷を軽減することができる。

また、本発明による特異値分解装置では、前記特異ベクトル算出部は、ｑｄ型ツイスト分解法により特異ベクトルを算出してもよい。
このような構成により、高速で高精度な特異値分解を実現することができうる。また、並列性にも優れている。

また、本発明による特異値分解装置では、前記特異ベクトル算出部は、ＬＶ型ツイスト分解法により特異ベクトルを算出してもよい。
このような構成により、数値安定的に特異ベクトルの算出を行うことができうる。

また、本発明による特異値分解装置では、前記特異ベクトル算出部は、前記対角行列記憶部から２重対角行列Ｂを読み出し、前記特異値記憶部から前記２重対角行列Ｂの特異値を読み出し、前記２重対角行列Ｂの各要素に関してミウラ変換、ｄＬＶｖ型変換、逆ミウラ変換を行うことによって、前記２重対角行列Ｂを上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）にコレスキー分解するコレスキー分解部と、前記上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）の各要素と、前記２重対角行列Ｂの特異値とを用いて一方の左右直交行列を構成する特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する第１特異ベクトル算出部と、前記第１特異ベクトル算出部が算出した一方の左右直交行列を構成する特異ベクトルと、前記２重対角行列Ｂの特異値と、前記２重対角行列Ｂとを用いて他方の左右直交行列を構成する特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する第２特異ベクトル算出部と、をさらに備えてもよい。
このような構成により、コレスキー分解の処理において、複数の補助変数を用いることによって、数値安定的に特異ベクトルの算出を行うことができうる。

また、本発明による特異値分解装置では、前記コレスキー分解部は、複数のコレスキー分解手段を備え、前記複数のコレスキー分解手段が、前記２重対角行列Ｂをコレスキー分解する処理を並列実行してもよい。
このような構成により、コレスキー分解の処理を短時間で行うことができうる。

また、本発明による特異値分解装置では、前記第１特異ベクトル算出部は、複数の第１特異ベクトル算出手段を備え、前記複数の第１特異ベクトル算出手段が、特異ベクトルを算出する処理を並列実行してもよい。
このような構成により、特異ベクトルを算出する処理を短時間で行うことができうる。

また、本発明による特異値分解装置では、前記第２特異ベクトル算出部は、複数の第２特異ベクトル算出手段を備え、前記複数の第２特異ベクトル算出手段が、特異ベクトルを算出する処理を並列実行してもよい。
このような構成により、特異ベクトルを算出する処理を短時間で行うことができうる。

また、本発明による特異値分解装置では、前記特異値算出部は、複数の特異値算出手段を備え、前記複数の特異値算出手段が、分割元の２重対角行列の特異値と、行列要素とを算出する処理を並列実行してもよい。
このような構成により、特異値を算出する処理を短時間で行うことができうる。

また、本発明による特異値分解装置では、前記特異値分解部は、複数の特異値分解手段を備え、前記複数の特異値分解手段が、２重対角行列に対して特異値分解を行う処理を並列実行してもよい。
このような構成により、特異値分解する処理を短時間で行うことができうる。

また、本発明による特異値分解装置では、行列Ａが記憶される行列記憶部と、前記行列Ａを前記行列記憶部から読み出し、前記行列Ａを２重対角化した前記２重対角行列Ｂを算出して前記対角行列記憶部に蓄積する対角化部と、をさらに備えてもよい。
このような構成により、任意の行列Ａについて特異値を算出することができる。また、２重対角行列Ｂの特異ベクトルを用いることにより、行列Ａの特異ベクトルを算出することもできうる。

また、本発明による特異値分解装置では、前記行列分割部は、２重対角行列を略半分の２個の２重対角行列に分割してもよい。
このような構成により、並列処理を適切に行うことができうる。

本発明による特異値分解装置等によれば、高速で高精度な特異値分解を行うことができる。また、並列性にも優れている。

以下、本発明による特異値分解装置について、実施の形態を用いて説明する。なお、以下の実施の形態において、同じ符号を付した構成要素及びステップは同一または相当するものであり、再度の説明を省略することがある。

（実施の形態１）
本発明の実施の形態１による特異値分解装置について、図面を参照しながら説明する。
図１は、本実施の形態による特異値分解装置１の構成を示すブロック図である。図１において、本実施の形態による特異値分解装置１は、行列記憶部１１と、対角化部１２と、対角行列記憶部１３と、行列分割部１４と、特異値分解部１５と、特異値分解記憶部１６と、特異値算出部１７と、特異値記憶部１８と、特異ベクトル算出部１９と、特異ベクトル記憶部２０とを備える。

行列記憶部１１では、任意の行列Ａが記憶される。この行列Ａは、各要素が実数である実行列である。なお、行列Ａが記憶されているとは、行列Ａを示すデータが記憶されている、という意味である。後述する記憶部においても同様である。行列記憶部１１は、所定の記録媒体（例えば、半導体メモリや磁気ディスク、光ディスクなど）によって実現されうる。行列記憶部１１での記憶は、ＲＡＭ等における一時的な記憶でもよく、あるいは、長期的な記憶でもよい。行列記憶部１１に行列Ａが記憶される過程は問わない。例えば、記録媒体を介して行列Ａが行列記憶部１１で記憶されるようになってもよく、通信回線等を介して送信された行列Ａが行列記憶部１１で記憶されるようになってもよく、あるいは、キーボードやマウス等の入力デバイスを介して入力された行列Ａが行列記憶部１１で記憶されるようになってもよい。

対角化部１２は、行列Ａを行列記憶部１１から読み出し、その読み出した行列Ａを２重対角化した２重対角行列Ｂを算出する。そして、対角化部１２は、その算出した２重対角行列Ｂを対角行列記憶部１３に蓄積する。対角化部１２は、例えば、ハウスホルダー（Ｈｏｕｓｅｈｏｌｄｅｒ）変換を必要なだけ繰り返し行う方法や、その他の２重対角化法を用いて、行列Ａを２重対角化する。ここで、２重対角行列Ｂは、上２重対角行列であってもよく、下２重対角行列であってもよい。本実施の形態では、２重対角行列Ｂが上２重対角行列である場合について説明する。

対角行列記憶部１３では、２重対角行列Ｂが記憶される。対角行列記憶部１３は、所定の記録媒体（例えば、半導体メモリや磁気ディスク、光ディスクなど）によって実現されうる。対角行列記憶部１３での記憶は、ＲＡＭ等における一時的な記憶でもよく、あるいは、長期的な記憶でもよい。

行列分割部１４は、対角行列記憶部１３から２重対角行列Ｂを読み出し、その２重対角行列Ｂを２個の２重対角行列に分割して対角行列記憶部１３に蓄積する。行列分割部１４は、その２重対角行列を２個の２重対角行列に分割して対角行列記憶部１３に蓄積する処理を、分割後の各２重対角行列があらかじめ決められた大きさ以下となるまで再帰的に繰り返す。

特異値分解部１５は、あらかじめ決められた大きさ以下の各２重対角行列を対角行列記憶部１３から読み出し、その各２重対角行列に対して特異値分解を行い、各２重対角行列の特異値と、その各２重対角行列の特異ベクトルを算出する。特異値分解部１５は、例えば、２分法と逆反復法を組み合わせた方法、ＭＲ＾３法、ＱＲｓ法等を用いて特異値分解を行ってもよい。ここで、ＭＲ＾３法において、例えば、ｄｑｄｓやｐｑｄｓ等のｑｄｓ法によって特異値を算出してもよい。また、ＱＲｓ法による特異値分解を行う場合には、ＦＯＲＴＲＡＮにおいて提供されているＤＢＤＳＱＲを用いてもよい。また、ｑｄｓ法によって特異値を算出する場合には、ＦＯＲＴＲＡＮにおいて提供されているＤＬＡＳＱを用いてもよい。
これらの特異値分解の方法については、すでに公知であり、その詳細な説明を省略する。特異値分解部１５は、算出した特異値を特異値分解記憶部１６に蓄積する。また、特異値分解部１５は、算出した特異ベクトルからなる左右直交行列の一部の要素である行列要素も特異値分解記憶部１６に蓄積する。特異ベクトルからなる左右直交行列とは、特異ベクトルを各列とする左直交行列、及び特異ベクトルを各列とする右直交行列のことである。行列要素の詳細については後述する。

特異値分解記憶部１６では、特異値分解部によって特異値分解された特異値と、前述の行列要素とが記憶される。また、特異値分解記憶部１６では、特異値算出部１７によって算出された特異値や行列要素も記憶される。特異値分解記憶部１６は、所定の記録媒体（例えば、半導体メモリや磁気ディスク、光ディスクなど）によって実現されうる。特異値分解記憶部１６での記憶は、ＲＡＭ等における一時的な記憶でもよく、あるいは、長期的な記憶でもよい。

特異値算出部１７は、特異値分解部１５によって算出された特異値と、行列要素とを特異値分解記憶部１６から読み出し、その特異値と行列要素とから、分割元の２重対角行列の特異値と、分割元の２重対角行列の行列要素とを算出して特異値分解記憶部１６に蓄積する。特異値算出部１７は、その分割元の２重対角行列の特異値と、行列要素とを算出する処理を、２重対角行列Ｂの特異値を算出するまで再帰的に繰り返す。そして、特異値算出部１７は、算出した２重対角行列Ｂの特異値を特異値記憶部１８に蓄積する。

特異値記憶部１８では、２重対角行列Ｂの特異値が記憶される。特異値記憶部１８は、所定の記録媒体（例えば、半導体メモリや磁気ディスク、光ディスクなど）によって実現されうる。特異値記憶部１８での記憶は、ＲＡＭ等における一時的な記憶でもよく、あるいは、長期的な記憶でもよい。

特異ベクトル算出部１９は、対角行列記憶部１３から２重対角行列Ｂを読み出し、特異値記憶部１８から２重対角行列Ｂの特異値を読み出す。そして、特異ベクトル算出部１９は、２重対角行列Ｂとその特異値とから、ツイスト分解法を用いて２重対角行列Ｂの特異ベクトルを算出して特異ベクトル記憶部２０に蓄積する。特異ベクトル算出部１９は、それらの処理を実行する、コレスキー分解部２１と、第１特異ベクトル算出部２２と、第２特異ベクトル算出部２３とを備える。特異ベクトル算出部１９は、ＬＶ型ツイスト分解法により特異ベクトルを算出してもよく、あるいは、ｑｄ型ツイスト分解法により特異ベクトルを算出してもよい。本実施の形態では、前者の場合について説明する。

コレスキー分解部２１は、対角行列記憶部１３から２重対角行列Ｂを読み出し、特異値記憶部１８から２重対角行列Ｂの特異値を読み出す。そして、コレスキー分解部２１は、２重対角行列Ｂの各要素に関してミウラ変換、ｄＬＶｖ型変換、逆ミウラ変換を行うことによって、２重対角行列Ｂを上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）にコレスキー分解する。この処理の詳細については後述する。

第１特異ベクトル算出部２２は、コレスキー分解部２１が算出した上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）の各要素と、２重対角行列Ｂの特異値とを用いて一方の左右直交行列を構成する特異ベクトルを算出して特異ベクトル記憶部２０に蓄積する。

第２特異ベクトル算出部２３は、特異値記憶部１８から２重対角行列Ｂの特異値を読み出す。そして、第２特異ベクトル算出部２３は、第１特異ベクトル算出部２２が算出した一方の左右直交行列を構成する特異ベクトルと、２重対角行列Ｂの特異値と、２重対角行列Ｂとを用いて他方の左右直交行列を構成する特異ベクトルを算出して特異ベクトル記憶部２０に蓄積する。このように、第１特異ベクトル算出部２２と、第２特異ベクトル算出部２３とによって、左右直交行列のそれぞれが算出されることになる。

特異ベクトル記憶部２０では、２重対角行列Ｂの特異ベクトルが記憶される。特異ベクトル記憶部２０は、所定の記録媒体（例えば、半導体メモリや磁気ディスク、光ディスクなど）によって実現されうる。特異ベクトル記憶部２０での記憶は、ＲＡＭ等における一時的な記憶でもよく、あるいは、長期的な記憶でもよい。

なお、行列記憶部１１、対角行列記憶部１３、特異値分解記憶部１６、特異値記憶部１８、特異ベクトル記憶部２０の任意の２以上の記憶部は、同一の記録媒体によって実現されてもよく、あるいは、別々の記録媒体によって実現されてもよい。前者の場合には、例えば、行列Ａの記憶されている領域が行列記憶部１１となり、２重対角行列Ｂ等の記憶されている領域が対角行列記憶部１３となる。

また、行列記憶部１１、対角行列記憶部１３、特異値分解記憶部１６、特異値記憶部１８、特異ベクトル記憶部２０の各記憶部は、２以上の記録媒体から構成されてもよい。

次に、本実施の形態による特異値分解装置１の動作について、図２のフローチャートを用いて説明する。
（ステップＳ１０１）対角化部１２は、行列記憶部１１で記憶されている行列Ａを読み出し、その行列Ａを２重対角化して２重対角行列Ｂを算出して対角行列記憶部１３に蓄積する。

（ステップＳ１０２）行列分割部１４、特異値分解部１５、特異値算出部１７によって２重対角行列Ｂの特異値が算出され、特異値記憶部１８に蓄積される。この処理の詳細については後述する。

（ステップＳ１０３）特異ベクトル算出部１９は、対角行列記憶部１３から２重対角行列Ｂを読み出し、特異値記憶部１８から２重対角行列Ｂの特異値を読み出し、２重対角行列Ｂの特異ベクトルを算出して特異ベクトル記憶部２０に蓄積する。この処理の詳細については後述する。

このようにして、２重対角行列Ｂの特異値分解が終了する。ここで、行列Ａの特異値は２重対角行列Ｂの特異値に等しいため、行列Ａの特異値も算出されたことになる。また、後述するように、所定の変換を行うことによって、行列Ａの特異ベクトルも２重対角行列Ｂの特異ベクトルから容易に算出することができる。

次に、図２のフローチャートのステップＳ１０２の処理について、図３のフローチャートを用いて説明する。
（ステップＳ２０１）行列分割部１４は、対角行列記憶部１３から２重対角行列Ｂを読み出し、その２重対角行列Ｂを２個の２重対角行列に分割して対角行列記憶部１３に蓄積する。行列分割部１４は、その２重対角行列を２個の２重対角行列に分割する処理を、分割後の各２重対角行列があらかじめ決められた大きさ以下となるまで繰り返す。

（ステップＳ２０２）特異値分解部１５は、対角行列記憶部１３で記憶されているあらかじめ決められた大きさ以下の各２重対角行列について、特異値分解を行う。特異値分解部１５は、特異値分解の結果である特異値と、特異ベクトルからなる左右直交行列の一部の要素である行列要素とを特異値分解記憶部１６に蓄積する。

（ステップＳ２０３）特異値算出部１７は、２重対角行列の特異値と、行列要素とを特異値分解記憶部１６から読み出し、分割元の２重対角行列の特異値と、分割元の２重対角行列の行列要素とを算出して特異値分解記憶部１６に蓄積する。特異値算出部１７は、その分割元の２重対角行列の特異値と、行列要素とを算出する処理を、２重対角行列Ｂの特異値を算出するまで繰り返し、２重対角行列Ｂの特異値を特異値記憶部１８に蓄積する。このようにして、特異値を算出する処理が終了する。

次に、図２のフローチャートのステップＳ１０３の処理について、図４のフローチャートを用いて説明する。
（ステップＳ３０１）コレスキー分解部２１は、対角行列記憶部１３から２重対角行列Ｂを読み出し、特異値記憶部１８から２重対角行列Ｂの特異値を読み出す。そして、コレスキー分解部２１は、２重対角行列Ｂの各要素に関してミウラ変換、ｄＬＶｖ型変換、逆ミウラ変換を行うことによって、２重対角行列Ｂを上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）にコレスキー分解する。

（ステップＳ３０２）第１特異ベクトル算出部２２は、上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）の各要素を用いて、一方の左右直交行列を構成する特異ベクトルを算出する。本実施の形態では、第１特異ベクトル算出部２２は、右直交行列を構成する右特異ベクトルを算出するものとする。

（ステップＳ３０３）第１特異ベクトル算出部２２は、算出した特異ベクトルを正規化する。すなわち、第１特異ベクトル算出部２２は、算出した特異ベクトルのノルムを算出し、特異ベクトルを算出したノルムで割ったものを最終的な特異ベクトルとして特異ベクトル記憶部２０に蓄積する。

（ステップＳ３０４）第２特異ベクトル算出部２３は、第１特異ベクトル算出部２２が算出した特異ベクトルと、２重対角行列Ｂの特異値と、２重対角行列Ｂとを用いて、第１特異ベクトル算出部２２が算出した特異ベクトルと異なる方の特異ベクトルを算出する。２重対角行列Ｂを特異値分解した結果が、第１特異ベクトル算出部２２が算出した特異ベクトルからなる直交行列と、第２特異ベクトル算出部２３が算出する特異ベクトルからなる直交行列と、特異値記憶部１８が記憶している特異値であるため、その性質を用いて、第２特異ベクトル算出部２３は、特異ベクトルを算出することができる。

（ステップＳ３０５）第２特異ベクトル算出部２３は、算出した特異ベクトルを正規化する。すなわち、第２特異ベクトル算出部２３は、算出した特異ベクトルのノルムを算出し、特異ベクトルを算出したノルムで割ったものを最終的な特異ベクトルとして特異ベクトル記憶部２０に蓄積する。このようにして、２重対角行列Ｂを特異値分解する処理は終了となる。
なお、この後に、行列記憶部１１が記憶している行列Ａの特異ベクトルを算出する処理が行われてもよいが、ここでは省略している。

また、特異ベクトル算出部１９が算出した特異ベクトルの精度を上げるために、図５で示されるように、特異ベクトルに関する処理を実行してもよく、あるいは、実行しなくてもよい。すなわち、図示しない逆反復法処理部は、特異ベクトル記憶部２０から特異ベクトル算出部１９が算出した特異ベクトルを読み出し、その特異ベクトルに対して逆反復法の処理を実行し、その結果の特異ベクトルを特異ベクトル記憶部２０に蓄積する（ステップＳ４０１）。次に、図示しないグラムシュミット処理部は、高精度が必要かどうか判断する（ステップＳ４０２）。この判断は、あらかじめ高精度が必要かどうか設定されている記録媒体等から、その設定を読み出すことによって判断してもよい。そして、図示しないグラムシュミット処理部は、高精度が必要な場合には、特異ベクトル記憶部２０から逆反復法の処理の実行された特異ベクトルを読み出し、その特異ベクトルに対してグラムシュミット法の処理を実行し、その結果の特異ベクトルを特異ベクトル記憶部２０に蓄積する（ステップＳ４０３）。なお、逆反復法や、グラムシュミット法については、すでに公知であり、その詳細な説明を省略する。

次に、本実施の形態による特異値分解装置１の動作について、以下、より詳細に説明する。
［行列Ａの対角化］
行列Ａは、例えば、ハウスホルダー変換等を用いて、以下に示されるように２重対角化することができることが知られている。ここでは、上２重対角化する場合について示すが、下２重対角化も同様にして行うことができる。ここで、２重対角行列Ｂは、行の数と列の数とが一致する正方行列である。

したがって、対角化部１２は、上記のようにして、行列記憶部１１から行列Ａを読み出し、２重対角行列Ｂを算出することができる。その算出された２重対角行列は、対角行列記憶部１３で記憶される（ステップＳ１０１）。

［２重対角行列Ｂの分割］
２重対角行列Ｂを分割する処理について説明する。まず、図６で示されるように、上２重対角行列であって、正方行列であるＢの右端に、全ての要素が０である列を加えたものを新たに上２重対角行列Ｂとする。なお、この新たな２重対角行列Ｂの特異値は、元の正方行列である２重対角行列Ｂの特異値と同じである。したがって、今後、この新たなｎ×（ｎ＋１）行列の特異値を求める処理について説明する。

図７で示されるように、ｎ×（ｎ＋１）の上２重対角行列Ｂが与えられた場合に、それらを２個の上２重対角行列と、２個の要素（図７では、ｂ_７，ｂ_８）とに分けることができる。したがって、上２重対角行列Ｂは、以下のように分割されることになる。

ただし、上２重対角行列Ｂがｎ×（ｎ＋１）行列であるため、上２重対角行列Ｂ_１は、（ｋ−１）×ｋ行列であり、上２重対角行列Ｂ_２は、（ｎ−ｋ）×（ｎ−ｋ＋１）行列である。ｋは、１<ｋ<ｎとなる整数である。ｅ_ｊは、適切な次元におけるｊ番目の単位ベクトルである。ここで、並列処理を適切に実行するためには、ｋを、ｎ／２を超えない最大の整数にとるか、あるいは、ｎ／２を下まわらない最小の整数にとることが好適である（このｋの値を用いて行列を２個の行列に分割する場合を、「行列を略半分の２個の行列に分割する」と呼ぶことにする）。

本実施の形態では、ｋを上式のようにとるものとする。なお、ｋの値は、前述のように、１<ｋ<ｎの範囲内で任意であることは言うまでもない。

したがって、行列分割部１４は、上述のようにして、対角行列記憶部１３が記憶している上２重対角行列Ｂを２個の上２重対角行列と、２個の要素に分割することができ、その分割の処理を繰り返すことができる。図８は、図３のフローチャートのステップＳ２０１における行列分割部１４による行列を分割する処理を示すフローチャートである。

（ステップＳ５０１）行列分割部１４は、カウンタＩを「１」に設定する。
（ステップＳ５０２）行列分割部１４は、対角行列記憶部１３からＩ番目の分割をしていない２重対角行列を読み出し、その２重対角行列を２個の２重対角行列と、２個の要素とに分割する。そして、行列分割部１４は、分割した２個の２重対角行列と、２個の要素とを対角行列記憶部１３に蓄積する。

（ステップＳ５０３）行列分割部１４は、Ｉ番目の分割を行っていない２重対角行列が、対角行列記憶部１３で記憶されているかどうか判断する。そして、Ｉ番目の分割を行っていない２重対角行列が、対角行列記憶部１３で記憶されている場合には、ステップＳ５０２に戻り、そうでない場合には、ステップＳ５０４に進む。

（ステップＳ５０４）行列分割部１４は、Ｉ番目の分割を行った２重対角行列の大きさがあらかじめ決められている大きさ以下かどうか判断する。行列分割部１４は、例えば、目的とする行列の大きさ（例えば、２５×２６など）の記憶されている図示しない記録媒体から、その目的とする行列の大きさを読み出し、対角行列記憶部１３で記憶されている、Ｉ番目の分割後の２重対角行列がその大きさ以下であるかどうかを判断してもよい。そして、Ｉ番目の分割を行った２重対角行列の大きさがあらかじめ決められている大きさ以下である場合には、２重対角行列を分割する処理は終了となり、そうでない場合には、ステップＳ５０５に進む。
（ステップＳ５０５）行列分割部１４は、カウンタＩを１だけインクリメントする。そして、ステップＳ５０２に戻る。

なお、この図８のフローチャートでは、上述のように、行列分割部１４が各行列を、略半分の２個の行列に分割するため、Ｉ番目の分割後の各２重対角行列の大きさがほとんど同じである場合について説明したが、行列分割部１４が各行列を、略半分の２個の行列に分割しない場合には、行列分割部１４は、分割後の各行列があらかじめ決められた大きさ以下となるように、分割を繰り返すものとする。

また、図８のフローチャートでは、ステップＳ５０４において、行列分割部１４が、対角行列記憶部１３で記憶されている分割後の行列の大きさをあらかじめ決められた大きさと比較する処理を実行する場合について説明したが、これは一例であって、行列分割部１４は、ステップＳ５０４において、それ以外の処理を行ってもよい。例えば、上述のように、行列分割部１４が各行列を、略半分の２個の行列に分割する場合には、元の２重対角行列Ｂの大きさを知ることができれば、何番目の分割で目的とする行列の大きさとなるのかを知ることができる。したがって、Ｎ番目（Ｎは１以上の整数）の分割で目的とする行列の大きさとなる場合には、ステップＳ５０４において、ＩがＮであるかどうかを比較し、ＮでなければステップＳ５０５に進み、Ｎであれば一連の処理を終了するようにしてもよい。

図９は、行列の分割について説明するための図である。まず、行列分割部１４は、１番目の分割として、２重対角行列Ｂを、２重対角行列Ｂ_１と、２重対角行列Ｂ_２とに分割する（ステップＳ５０１，Ｓ５０２）。２重対角行列Ｂは、１個しかないため、行列分割部１４は、１番目の分割をしていない行列がないと判断する（ステップＳ５０３）。また、２重対角行列Ｂ_１等はあらかじめ決められた大きさ以下の行列でないとすると（ステップＳ５０４）、行列分割部１４は、２番目の分割として、２重対角行列Ｂ_１を、２重対角行列Ｂ_１１と、２重対角行列Ｂ_１２とに分割する（ステップＳ５０５，Ｓ５０２）。この場合には、２番目の分割を行っていない２重対角行列Ｂ_２が存在するため、行列分割部１４は、２重対角行列Ｂ_２も、２重対角行列Ｂ_２１と、２重対角行列Ｂ_２２とに分割する（ステップＳ５０３，Ｓ５０２）。このようにして、分割後の各２重対角行列が目的とする大きさ以下になるまで、行列を分割する処理が繰り返される。なお、図９では、２重対角行列以外の２個の要素については省略している。また、図９において、各２重対角行列は、列の数が行の数よりも１だけ大きい行列である。

［分割された行列の特異値分解］
行列Ｂ_ｉがｎ×（ｎ＋１）行列であるとすると、行列Ｂ_ｉの特異値分解は、次のようになる。

ここで、Ｄ_ｉは、ｎ×ｎの対角行列である。（Ｄ_ｉ０）は、行列Ｄ_ｉの右側に全ての要素が０の列が１つある行列である。Ｄ_ｉの各対角成分はＢ_ｉの特異値である。また、ｖ_ｉは、行列Ｂ_ｉの特異値分解における右直交行列の一番右側のベクトルである。Ｖ_ｉは、行列Ｂ_ｉの特異値分解における右直交行列のｖ_ｉをのぞいた行列である。（Ｖ_ｉ ｖ_ｉ）^Ｔは、全体として（ｎ＋１）×（ｎ＋１）の行列である。

したがって、特異値分解部１５は、対角行列記憶部１３から分割後の各上２重対角行列を読み出し、上記のようにして、特異値分解を行う（ステップＳ２０２）。特異値分解の方法として、例えば、２分法と逆反復法を組み合わせた方法、ＭＲ＾３法、ＱＲｓ法等を用いてもよいことは前述の通りである。図９で示されるように上２重対角行列の分割が行われた場合には、特異値分解部１５は、各２重対角行列Ｂ_１１１，Ｂ_１１２，Ｂ_１２１，Ｂ_１２２，……，Ｂ_２２２に対して特異値分解を行うことになる。なお、特異値分解部１５は、特異値分解の結果得られた各特異値と、特異値分解の結果得られた左右直交行列の一部の要素である行列要素とを特異値分解記憶部１６に蓄積する。ここで、行列要素が、左右直交行列のどの要素を含むのかについては後述する。

［特異値の算出］
まず、図１０で示されるように、分割された２個の２重対角行列Ｂ_１，Ｂ_２から、分割元の２重対角行列Ｂ_０の特異値等を算出する処理について説明する。２重対角行列Ｂ_１，Ｂ_２は、次のように特異値分解されていたとする。ここで、２重対角行列Ｂ_１，Ｂ_２は、両者共に列の数が行の数よりも１だけ大きい行列であるとする。

この場合に、分割元の２重対角行列Ｂ_０は、次のようになる。

また、上式のｂ_２ｋ−１等は２重対角行列の分割の処理において説明した行列の要素である。ここで、Ｇｉｖｅｎｓ変換を行い、ｂ_２ｋφ_２を０にすると、次のようになる。

ただし、

である。したがって、行列Ｂ_０は、

によって、（Ｍ ０）に直交変換される。これらより、行列Ｂ_０の特異値分解を行うためには、行列Ｍの特異値分解を行えばよいことになる。ここで、ｎ×ｎの行列Ｍを次のようにおく。

上記の行列Ｍにおいて、ｚ_ｉ、ｄ_ｊ以外の各要素は０である。ただし、ｉ＝１，２，…，ｎ、ｊ＝２，３，…，ｎである。この行列Ｍの特異値分解は、次の定理によって行うことができる。
（定理１）

ここで、計算機によって求めることができるのは、行列Ｍの真の特異値ではなく、誤差を含む近似値である。したがって、

は、真の値である

と大きく異なる可能性がありうる。すなわち、上記のように計算した場合には、特異ベクトルの計算が数値不安定となることがある。この欠点は、次の定理によって克服できる。

（定理２）

したがって、定理１を用いて行列Ｍの近似特異値を算出した後に、定理２を用いて、その近似特異値を真の特異値とする行列Ｍを再構成する。特異ベクトルは、上記の式５を用いて算出した

と、上記の式２を用いて算出した

とを上記の式３，式４に代入することによって数値安定的に求めることができる。このようにして、次のように行列Ｍの特異値分解が求められたとする。

すると、Ｂ_０の特異値分解は以下のようになる。

このように、対象となる行列を２個の副行列に分割し、その分割した後の各副行列について特異値分解を行う処理を再帰的に行うことによって特異値分解を行う方法は、分割統治法（Ｄｉｖｉｄｅ ａｎｄ Ｃｏｎｑｕｅｒ：Ｄ&Ｃ）と呼ばれている。

分割統治法では、特異ベクトルまで算出することになり、その特異ベクトルを算出する処理において行列の計算をしなければならないため、非常に負荷の大きい処理となる。分割統治法において特異値分解をする場合には、例えば、計算時間の９５％程度が特異ベクトルを算出するベクトル更新の処理（例えば、上述のように、行列Ｍの左右直交行列から、行列Ｂ_０の左右直交行列を算出するために行列を掛け合わせる処理）に費やされることもありうる。

一方、特異値のみを算出する場合には、上記処理を簡略化することができる。次に、その方法について説明する。上記の結果より、次のようになる。

これより、

となる。ここで、ｆとφとは、行列Ｂ_０を特異値分解した結果の右直交行列の最初の行の各要素である。また、ｌとψとは、行列Ｂ_０を特異値分解した結果の右直交行列の最後の行の各要素である。この右直交行列の最初の行の各要素と、最後の行の各要素とが行列要素となる。

上記の結果から、図１０で示されるように、行列Ｂ_１，Ｂ_２から分割元の行列Ｂ_０を構成する場合に、行列Ｂ_０について、

をそれぞれ算出することができる。このような処理を繰り返すことにより、最終的に、行列Ａを２重対角化した２重対角行列Ｂの特異値を算出することができる。

したがって、特異値算出部１７は、上述のようにして、特異値分解記憶部１６から特異値と行列要素とを読み出し、対角行列記憶部１３から行列の分割時に発生した２個の要素（ｂ_２ｋ−１、ｂ_２ｋ）を読み出し、それらを用いることによって、分割元の２重対角行列の特異値と行列要素とを算出する。そして、それらを算出する処理を繰り返すことによって、最終的に２重対角行列Ｂの特異値を算出する。図１１は、図３のフローチャートのステップＳ２０３における特異値算出部１７が特異値を算出する処理を示すフローチャートである。

（ステップＳ６０１）特異値算出部１７は、カウンタＪを「１」に設定する。
（ステップＳ６０２）特異値算出部１７は、Ｊ番目の特異値の算出は最後の特異値の算出であるかどうか判断する。ここで、最後の特異値の算出とは、２重対角行列Ｂの特異値を算出することである。そして、最後の特異値の算出である場合には、ステップＳ６０６に進み、そうでない場合には、ステップＳ６０３に進む。

（ステップＳ６０３）特異値算出部１７は、分割元の２重対角行列の特異値等を算出する。この処理の詳細については後述する。
（ステップＳ６０４）特異値算出部１７は、Ｊ番目の特異値の算出において、分割元の全ての２重対角行列の特異値等を算出したかどうか判断する。そして、Ｊ番目の特異値の算出において、分割元の全ての２重対角行列の特異値等を算出した場合には、ステップＳ６０５に進み、そうでない場合には、ステップＳ６０３に戻る。

（ステップＳ６０５）特異値算出部１７は、カウンタＪを１だけインクリメントする。そして、ステップＳ６０２に戻る。
（ステップＳ６０６）特異値算出部１７は、２重対角行列Ｂの特異値を算出し、その算出した特異値を特異値記憶部１８に蓄積する。このようにして、２重対角行列Ｂの特異値を算出する一連の処理は終了となる。

図１２は、図１１のフローチャートにおけるステップＳ６０３の詳細な処理を示すフローチャートである。
（ステップＳ７０１）特異値算出部１７は、分割元の行列の特異値を式２を用いて算出する。そして、特異値算出部１７は、算出した特異値を特異値分解記憶部１６に蓄積する。

（ステップＳ７０２）特異値算出部１７は、ステップＳ７０１で算出した特異値を用いて、式５のｚを算出する。
（ステップＳ７０３）特異値算出部１７は、ステップＳ７０１で算出した特異値と、ステップＳ７０２で算出したｚとを用いて、式４のｖ_ｉを算出する。このｖ_ｉを算出することにより、ｖ_ｉを列ベクトルに有するＶ_Ｍを算出したことになる。そして、特異値算出部１７は、算出したＶ_Ｍを特異値分解記憶部１６に蓄積する。

（ステップＳ７０４）特異値算出部１７は、式６から式９を用いて、分割元の行列に関する行列要素を算出する。そして、特異値算出部１７は、算出した分割元の行列要素を特異値分解記憶部１６に蓄積する。このようにしてステップＳ６０３の処理は終了となる。

図１３は、特異値を算出する処理について説明するための図である。特異値分解記憶部１６では、特異値分解部１５によって行列Ｂ_１１１，Ｂ_１１２，...，Ｂ_２２２が特異値分解された結果、すなわち、それらの特異値と、行列要素とが記憶されているものとする。まず、特異値算出部１７は、特異値等の１番目の算出を開始する（ステップＳ６０１，Ｓ６０２）。ここで、特異値等の１番目の算出とは、図１３の一番下の行の行列の特異値と、行列要素とから、下から２番目の行の行列の特異値と、行列要素とを算出することである。特異値算出部１７は、行列Ｂ_１１１と、行列Ｂ_１１２との各特異値、及び行列Ｂ_１１１と、行列Ｂ_１１２との各行列要素を特異値分解記憶部１６から読み出す。また、行列Ｂ_１１を行列Ｂ_１１１と、行列Ｂ_１１２とに分解したときに発生した２個の要素を対角行列記憶部１３から読み出す。そして、特異値算出部１７は、それらの値を用いて、行列Ｂ_１１の特異値を算出して特異値分解記憶部１６に蓄積する（ステップＳ７０１）。次に、特異値算出部１７は、行列Ｂ_１１の特異値を用いてｚの値を算出する（ステップＳ７０２）。特異値算出部１７は、行列Ｂ_１１の特異値と、ｚの値とを用いて直交行列Ｖ_Ｍを算出して特異値分解記憶部１６に蓄積する（ステップＳ７０３）。最後に、特異値算出部１７は、行列Ｂ_１１の行列要素を算出して特異値分解記憶部１６に蓄積する（ステップＳ７０４）。

特異値等の１番目の算出はまだ終わっていないため（ステップＳ６０４）、特異値算出部１７は、上記説明と同様にして、行列Ｂ_１２１、行列Ｂ_１２２の特異値、行列要素等に基づいて、行列Ｂ_１２の特異値、行列要素を算出する（ステップＳ６０２，Ｓ６０３）。このようにして、特異値等の１番目の算出が終了すると、特異値算出部１７は、特異値等の２番目の算出、すなわち、行列Ｂ_１１、行列Ｂ_１２の特異値、行列要素等に基づいて行列Ｂ_１の特異値、行列要素の算出を行う（ステップＳ６０４，Ｓ６０５，Ｓ６０２，Ｓ６０３）。その後、特異値算出部１７は、同様にして、行列Ｂ_２１、行列Ｂ_２２の特異値、行列要素等に基づいて行列Ｂ_２の特異値、行列要素の算出を行う（ステップＳ６０４，Ｓ６０３）。

特異値算出部１７が特異値等の２番目の算出を終了すると（ステップＳ６０４）、次の特異値等の算出は、２重対角行列Ｂの特異値を求める最後の処理となるため（ステップＳ６０５，Ｓ６０２）、特異値算出部１７は、特異値の算出のみを行い、その算出した特異値を特異値記憶部１８に蓄積する（ステップＳ６０６）。このようにして、特異値の算出が終了する。

［特異値からの特異ベクトルの算出］
まず、行列Ｂ^ＴＢを考える。ここで、行列Ｂは、前述のｎ×ｎの上２重対角行列である。この行列Ｂ^ＴＢを次のように対角化したとする。

ここで、一般に次のことが成り立つ。
（１）行列Ｂ^ＴＢは対称な３重対角行列である。
（２）行列Ｂ^ＴＢの固有値は全て正であり、行列Ｂの特異値σ_ｊ（σ_１≧σ_２≧…≧σ_ｍ≧０）は、行列Ｂ^ＴＢの固有値と次の関係を有する。

（３）Ｖ_Ｂ＝Ｖである。ただし、Ｖ_Ｂは行列Ｂの右直交行列である。したがって、行列Ｂ^ＴＢの固有ベクトルｘ_ｊは、行列Ｂの右特異ベクトルに等しい。

したがって、行列Ｂ^ＴＢの固有ベクトルを求めると、行列Ｂの右直交行列が求まることになる。さらに、行列Ｂの固有値分解を、

とする。すると、右直交行列Ｖ_Ｂが求まり、特異値を対角成分に有する行列Σが求まることによって、

から、左直交行列も求まることになり、行列Ｂが特異値分解される。したがって、行列Ｂの特異ベクトルを求めることは、行列Ｂ^ＴＢ＝Ｔ_Ｓの固有ベクトルを求めることに置き換えることができる。すなわち、次の方程式の固有ベクトルｘ_ｊを求めればよいことになる。

本来であれば、上記式の右辺は０になるはずであるが、行列Ｂの特異値を算出するときに、いくらかの誤差を含むため、特異ベクトルｘ_ｊが真値であれば、上記のように右辺に残差項が存在する。

ここで、以下のようにコレスキー分解できたとする。

すると、次式のように書ける。

そして、次のようになる。

ここで、行列Ｎ_ｋをツイスト行列と呼ぶ。また、

であるので、（Ｂ^ＴＢ−σ_ｊ ^２Ｉ）ｘ_ｊ＝γ_ｋｅ_ｋは、

となる。この簡単な式を解くことによって特異ベクトルを算出することができる。具体的には、

のようにわずかな演算で特異ベクトルを算出することができる。なお、あるρ０について、Ｄ_ρ０ ^＋＝０あるいはＤ_ρ０ ⁻＝０であったとしても、行列Ｂの（ρ，ρ）成分であるｂ_２ρ−１ ^（０）と、行列Ｂの（ρ，ρ＋１）成分であるｂ_２ρ ^（０）を用いて、

と特異ベクトルを算出することができる（このように特異ベクトルを算出する処理を例外処理と呼ぶことにする）。なお、残差項のパラメータγ_ｋ及びｋの値は、

の絶対値が最小となるように決定する。したがって、上述のコレスキー分解を求め、ツイスト行列Ｎ_ｋを求めることができれば、特異ベクトルを求めることができることになる。そこで、次にコレスキー分解について説明する。

図１４で示すように、Ｂ^（０）ＴＢ^（０）−σ_ｊ ^２Ｉをコレスキー分解することは、Ｂ^（０）に対応する｛ｑ_ｋ ^（０），ｅ_ｋ ^（０）｝から、Ｂ^（±１）に対応する｛ｑ_ｋ ^（±１），ｅ_ｋ ^（±１）｝を求めることである。

（ｑｄ型ツイスト分解法）
まず、従来から知られているｑｄ型ツイスト分解法で用いるｑｄ型変換について説明する（図１４参照）。

この変換は、さらにｓｔａｔｉｏｎａｒｙ ｑｄ ｗｉｔｈ ｓｈｉｆｔ（ｓｔｑｄｓ）変換

及びｒｅｖｅｒｓｅ−ｔｉｍｅ ｐｒｏｇｒｅｓｓｉｖｅ ｑｄ ｗｉｔｈ ｓｈｉｆｔ（ｒｐｑｄｓ）変換

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

（ＬＶ型ツイスト分解法）
次に、ＬＶ型ツイスト分解法で用いるミウラ変換、ｄＬＶｖ型変換、逆ミウラ変換について説明する（図１４参照）。

まず、ミウラ変換について説明する。この変換は、次のように示される。

次に、ｄＬＶｖ型変換について説明する。この変換は、さらにｓｔａｔｉｏｎａｒｙ ｄｉｓｃｒｅｔｅ Ｌｏｔｋａ−Ｖｏｌｔｅｒｒａ ｗｉｔｈ ｖａｒｉａｂｌｅ ｓｔｅｐ−ｓｉｚｅ（ｓｔｄＬＶｖ）変換

及びｒｅｖｅｒｓｅ−ｔｉｍｅ ｄｉｓｃｒｅｔｅ Ｌｏｔｋａ−Ｖｏｌｔｅｒｒａ ｗｉｔｈ ｖａｒｉａｂｌｅ ｓｔｅｐ−ｓｉｚｅ（ｒｄＬＶｖ）変換

に分けられる。ただし、

を満たす範囲内でδ^（±１）は任意である。

最後に、逆ミウラ変換について説明する。この変換は、次のように示される。

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

ここで、ＬＶ型ツイスト分解法による処理のより詳細な処理の一例について説明する。図１５〜図２０は、ＬＶ型ツイスト分解法による処理の一例を示すフローチャートである。

図１５は、コレスキー分解の全体の処理の一例を示す図である。
（ステップＳ９０１）コレスキー分解部２１は、ミウラ変換を行う。この処理の詳細については後述する。
（ステップＳ９０２）コレスキー分解部２１は、１／δ^（±１）を１／δ^（０）−σ_ｊ ^２とする。

（ステップＳ９０３）コレスキー分解部２１は、後述するＰｒｏｃｅｄｕｒｅ１の処理を実行する。
（ステップＳ９０４）コレスキー分解部２１は、後述するＰｒｏｃｅｄｕｒｅ２の処理を実行する。

（ステップＳ９０５）コレスキー分解部２１は、ｑ_ｋ ^（＋１），ｅ_ｋ ^（＋１）がすでに算出されているかどうか判断する。そして、すでに算出されている場合には、コレスキー分解の一連の処理は終了となり、一方、算出されていない場合には、ステップＳ９０１に戻る。

図１６は、図１５のフローチャートにおけるステップＳ９０３の処理の詳細を示すフローチャートである。
（ステップＳ１００１）コレスキー分解部２１は、ｑ_ｋ ^（＋１），ｅ_ｋ ^（＋１）がすでに算出されているかどうか判断する。そして、すでに算出されている場合には、終了となり、算出されていない場合には、ステップＳ１００２に進む。
（ステップＳ１００２）コレスキー分解部２１は、ｓｔｄＬＶｖ変換等の処理を行う。この処理の詳細については後述する。

図１７は、図１５のフローチャートにおけるステップＳ９０４の処理の詳細を示すフローチャートである。
（ステップＳ１１０１）コレスキー分解部２１は、ｑ_ｋ ^（−１），ｅ_ｋ ^（−１）がすでに算出されているかどうか判断する。そして、すでに算出されている場合には、終了となり、算出されていない場合には、ステップＳ１１０２に進む。
（ステップＳ１１０２）コレスキー分解部２１は、ｒｄＬＶｖ変換等の処理を行う。この処理の詳細については後述する。

図１８は、図１５のフローチャートのステップＳ９０１の処理の詳細を示すフローチャートである。
（ステップＳ１２０１）コレスキー分解部２１は、１／δ^（０）の値を決定する。この値は、前述のように任意に決定することができうる。算出された特異値の小さい順にσ_１，σ_２，…とする場合に、例えば、コレスキー分解部２１は、１／δ^（０）の値をσ_１ ^２よりも小さい値（例えば、「１」など）に設定し、その後、ステップＳ１２０３等で桁落ちの発生する可能性があると判断された場合に、１／δ^（０）の値をσ_ｊ ^２とσ_ｊ＋１ ^２の間（ｊ＝１，２，…）に設定する、というように、ｊを１ずつ大きくしながら１／δ^（０）の値を設定していってもよい。なお、１／δ^（０）の値を決定する方法は、これに限定されないことは言うまでもない。
（ステップＳ１２０２）コレスキー分解部２１は、β１をｑ_１ ^（０）−１／δ^（０）に設定する。

（ステップＳ１２０３）コレスキー分解部２１は、β１の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップＳ１２０４に進み、そうでない場合には、ステップＳ１２０１に戻って、１／δ^（０）の値を決定する処理を再度実行する。なお、このεも任意に定めることができる。例えば、εとしてマシン・イプシロンを用いてもよい。εの値を大きくすると精度が向上し、εの値を小さくすると精度が低下する。なお、このステップＳ１２０３で行っている処理は、桁落ちが発生する可能性について判断する処理である。β１の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。

（ステップＳ１２０４）コレスキー分解部２１は、β２をｑ_１ ^（０）／（１＋δ^（０）ｕ_０ ^（０））に設定する。なお、前述のようにｕ_０ ^（０）＝０であるため、β２はｑ_１ ^（０）に設定されたことになる。

（ステップＳ１２０５）コレスキー分解部２１は、ｕ_１ ^（０）（１＋δ^（０）ｕ_０ ^（０））をβ１に設定する。ここで、β１は、ステップＳ１２０２においてｑ_１ ^（０）−１／δ^（０）に設定されているが、前述のようにｕ_０ ^（０）＝０であるため、β１は、ｑ_１ ^（０）／（１＋δ^（０）ｕ_０ ^（０）−１／δ^（０）と等しく、これにミウラ変換を行うとｕ_１ ^（０）＝ｕ_２ｋ−１ ^（０）｜_ｋ＝１となるからである。なお、ｕ_０ ^（０）＝０から（１＋δ^（０）ｕ_０ ^（０））＝１である。

（ステップＳ１２０６）コレスキー分解部２１は、カウンタｋを「１」に設定する。
（ステップＳ１２０７）コレスキー分解部２１は、α１をｅ_ｋ ^（０）／β１に設定する。前述のように、β１はｕ_２ｋ−１ ^（０）となるから、ｅ_ｋ ^（０）／β１にミウラ変換を行うと、α１はδ^（０）ｕ_２ｋ ^（０）と等しいことになる。

（ステップＳ１２０８）コレスキー分解部２１は、α２をα１＋１に設定する。ステップＳ１２０７の説明からわかるように、α２は１＋δ^（０）ｕ_２ｋ ^（０）と等しいことになる。

（ステップＳ１２０９）コレスキー分解部２１は、α２の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップＳ１２１０に進み、そうでない場合には、ステップＳ１２０１に戻って、１／δ^（０）の値を決定する処理を再度実行する。なお、このステップＳ１２０９で行っている処理は、ステップＳ１２０３の処理と同様に、桁落ちが発生する可能性について判断する処理である。α２の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。

（ステップＳ１２１０）コレスキー分解部２１は、α１×β２を算出してｕ_２ｋ ^（０）（１＋δ^（０）ｕ_２ｋ−１ ^（０））に設定する。前述のように、α１はδ^（０）ｕ_２ｋ ^（０）と等しく、β２はｑ_ｋ ^（０）／（１＋δ^（０）ｕ_２ｋ−２ ^（０））であるため、α１×β２にミウラ変換を実行すると、ｕ_２ｋ ^（０）（１＋δ^（０）ｕ_２ｋ−１ ^（０））に等しくなるからである。

（ステップＳ１２１１）コレスキー分解部２１は、β２をｑ_ｋ＋１ ^（０）／α２に設定する。前述のように、α２は１＋δ^（０）ｕ_２ｋ ^（０）と等しいため、β２は、ｑ_ｋ＋１ ^（０）／（１＋δ^（０）ｕ_２ｋ ^（０））となり、β２におけるｋの値を１だけインクリメントしたことになる。

（ステップＳ１２１２）コレスキー分解部２１は、β１をβ２−１／δ^（０）に設定する。前述のように、β２はｑ_ｋ＋１ ^（０）／（１＋δ^（０）ｕ_２ｋ ^（０））と等しいため、β１は、ｑ_ｋ＋１ ^（０）／（１＋δ^（０）ｕ_２ｋ ^（０））−１／δ^（０）となり、ミウラ変換を実行すると、β１は、ｕ_２ｋ＋１ ^（０）となる。したがって、β１におけるｋの値を１だけインクリメントしたことになる。

（ステップＳ１２１３）コレスキー分解部２１は、β１の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップＳ１２１４に進み、そうでない場合には、ステップＳ１２０１に戻って、１／δ^（０）の値を決定する処理を再度実行する。なお、このステップＳ１２１３で行っている処理は、ステップＳ１２０３の処理と同様に、桁落ちが発生する可能性について判断する処理である。β１の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。

（ステップＳ１２１４）コレスキー分解部２１は、α２×β１を算出してｕ_２ｋ＋１ ^（０）（１＋δ^（０）ｕ_２ｋ ^（０））に設定する。前述のように、α２は１＋δ^（０）ｕ_２ｋ ^（０）と等しく、β１は、ｕ_２ｋ＋１ ^（０）に等しいからである。

（ステップＳ１２１５）コレスキー分解部２１は、カウンタｋを１だけインクリメントする。
（ステップＳ１２１６）コレスキー分解部２１は、カウンタｋがｍに等しいかどうか判断する。そして、ｍに等しい場合には、一連の処理は終了となり、そうでない場合には、ステップＳ１２０７に戻る。

このようにして、ｕ_２ｋ＋１ ^（０）（１＋δ^（０）ｕ_２ｋ ^（０））と、ｕ_２ｋ ^（０）（１＋δ^（０）ｕ_２ｋ−１ ^（０））とが算出されることになる。これらの値は、コレスキー分解部２１が有する図示しないメモリ等において一時的に記憶されてもよい。

図１９は、図１６のフローチャートのステップＳ１００２の処理の詳細を示すフローチャートである。
（ステップＳ１３０１）コレスキー分解部２１は、α２を１＋δ^（＋１）ｕ_０ ^（＋１）に設定する。なお、前述のようにｕ_０ ^（＋１）＝０であるため、α２は１に設定されたことになる。

（ステップＳ１３０２）コレスキー分解部２１は、β１をｕ_１ ^（０）（１＋δ^（０）ｕ_０ ^（０））／（１＋δ^（＋１）ｕ_０ ^（＋１））に設定する。ここで、ｕ_１ ^（０）（１＋δ^（０）ｕ_０ ^（０））の値としては、ステップＳ１２０５で算出したものを用いる。なお、前述のようにｕ_０ ^（＋１）＝０である。また、このβ１の式にｓｔｄＬＶｖ変換を実行すると、β１は、ｕ_１ ^（＋１）に設定されたことになる。

（ステップＳ１３０３）コレスキー分解部２１は、カウンタｋを「１」に設定する。
（ステップＳ１３０４）コレスキー分解部２１は、α１をβ１＋１／δ^（＋１）に設定する。

（ステップＳ１３０５）コレスキー分解部２１は、α１の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップＳ１３０６に進み、そうでない場合には、図１５のフローチャートのＰｒｏｃｅｄｕｒｅ２（ステップＳ９０４）に進む。なお、このステップＳ１３０５で行っている処理は、ステップＳ１２０３の処理と同様に、桁落ちが発生する可能性について判断する処理である。α１の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。

（ステップＳ１３０６）コレスキー分解部２１は、β２をｕ_２ｋ ^（０）（１＋δ^（０）ｕ_２ｋ−１ ^（０））／α１に設定する。ここで、ｕ_２ｋ ^（０）（１＋δ^（０）ｕ_２ｋ−１ ^（０））の値としては、ステップＳ１２１０で算出したものを用いる。また、このβ２の式にｓｔｄＬＶｖ変換を実行すると、β２は、δ^（＋１）ｕ_２ｋ ^（＋１）に設定されたことになる。

（ステップＳ１３０７）コレスキー分解部２１は、α１×α２を算出する。前述のように、α１はβ１＋１／δ^（＋１）＝ｕ_２ｋ−１ ^（＋１）＋１／δ^（＋１）に等しく、α２は１＋δ^（＋１）ｕ_２ｋ−２ ^（＋１）に等しいため、α１×α２に逆ミウラ変換を実行するとｑ_ｋ ^（＋１）に等しくなる。したがって、コレスキー分解部２１は、α１×α２を算出してｑ_ｋ ^（＋１）に設定する。

（ステップＳ１３０８）コレスキー分解部２１は、α２を１＋β２に設定する。前述のように、β２はδ^（＋１）ｕ_２ｋ ^（＋１）と等しいため、α２は１＋δ^（＋１）ｕ_２ｋ ^（＋１）となる。したがって、α２におけるｋの値を１だけインクリメントしたことになる。

（ステップＳ１３０９）コレスキー分解部２１は、α２の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップＳ１３１０に進み、そうでない場合には、図１５のフローチャートのＰｒｏｃｅｄｕｒｅ２（ステップＳ９０４）に進む。なお、このステップＳ１３０９で行っている処理は、ステップＳ１２０３の処理と同様に、桁落ちが発生する可能性について判断する処理である。α２の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。

（ステップＳ１３１０）コレスキー分解部２１は、β１×β２を算出する。前述のように、β１はｕ_２ｋ−１ ^（＋１）に等しく、β２はδ^（＋１）ｕ_２ｋ ^（＋１）と等しいため、β１×β２に逆ミウラ変換を実行すると、ｅ_ｋ ^（＋１）に等しくなる。したがって、コレスキー分解部２１は、β１×β２を算出してｅ_ｋ ^（＋１）に設定する。

（ステップＳ１３１１）コレスキー分解部２１は、β１をｕ_２ｋ＋１ ^（０）（１＋δ^（０）ｕ_２ｋ ^（０））／α２に設定する。ここで、ｕ_２ｋ＋１ ^（０）（１＋δ^（０）ｕ_２ｋ ^（０））の値としては、ステップＳ１２１４で算出したものを用いる。また、このβ１の式にｓｔｄＬＶｖ変換を実行すると、β１は、ｕ_２ｋ＋１ ^（＋１）に設定されたことになる。したがって、β１におけるｋの値を１だけインクリメントしたことになる。

（ステップＳ１３１２）コレスキー分解部２１は、カウンタｋを１だけインクリメントする。
（ステップＳ１３１３）コレスキー分解部２１は、カウンタｋがｍに等しいかどうか判断する。そして、ｍに等しい場合には、ステップＳ１３１４に進み、そうでない場合には、ステップＳ１３０４に戻る。

（ステップＳ１３１４）コレスキー分解部２１は、α２×（β１＋１／δ^（＋１））を算出する。これは、ステップＳ１３０４でα１を更新した後に、ステップＳ１３０７でα１×α２を計算することと等しい。したがって、コレスキー分解部２１は、α２×（β１＋１／δ^（＋１））を算出してｑ_１ ^（＋１）に設定する。このようにして、ミウラ変換された結果にｓｔｄＬＶｖ変換と、逆ミウラ変換とを実行して、ｑ_ｋ ^（＋１）、ｅ_ｋ ^（＋１）を算出する処理が終了する。これらの値は、コレスキー分解部２１が有する図示しないメモリ等において一時的に記憶されてもよい。

図２０は、図１７のフローチャートのステップＳ１１０２の処理の詳細を示すフローチャートである。
（ステップＳ１４０１）コレスキー分解部２１は、β１をｕ_２ｍ−１ ^（０）（１＋δ^（０）ｕ_２ｍ−２ ^（０））／（１＋δ^（−１）ｕ_２ｍ ^（−１））に設定する。ここで、ｕ_２ｍ−１ ^（０）（１＋δ^（０）ｕ_２ｍ−２ ^（０））の値としては、ステップＳ１２１４で算出したものを用いる。なお、前述のようにｕ_２ｍ ^（−１）＝０である。また、このβ１の式にｒｄＬＶｖ変換を実行すると、β１は、ｕ_２ｍ−１ ^（−１）に設定されたことになる。

（ステップＳ１４０２）コレスキー分解部２１は、カウンタｋを「ｍ−１」に設定する。
（ステップＳ１４０３）コレスキー分解部２１は、α１をβ１＋１／δ^（−１）に設定する。

（ステップＳ１４０４）コレスキー分解部２１は、α１の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップＳ１４０５に進み、そうでない場合には、図１５のフローチャートのミウラ変換（ステップＳ９０１）に戻る。なお、このステップＳ１４０４で行っている処理は、ステップＳ１２０３の処理と同様に、桁落ちが発生する可能性について判断する処理である。α１の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。

（ステップＳ１４０５）コレスキー分解部２１は、β２をｕ_２ｋ ^（０）（１＋δ^（０）ｕ_２ｋ−１ ^（０））／α１に設定する。ここで、ｕ_２ｋ ^（０）（１＋δ^（０）ｕ_２ｋ−１ ^（０））の値としては、ステップＳ１２１０で算出したものを用いる。また、このβ２の式にｒｄＬＶｖ変換を実行すると、β２は、δ^（−１）ｕ_２ｋ ^（−１）に設定されたことになる。

（ステップＳ１４０６）コレスキー分解部２１は、α２を１＋β２に設定する。前述のように、β２はδ^（−１）ｕ_２ｋ ^（−１）と等しいため、α２は１＋δ^（−１）ｕ_２ｋ ^（−１）となる。

（ステップＳ１４０７）コレスキー分解部２１は、α２の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップＳ１４０８に進み、そうでない場合には、図１５のフローチャートのミウラ変換（ステップＳ９０１）に戻る。なお、このステップＳ１４０７で行っている処理は、ステップＳ１２０３の処理と同様に、桁落ちが発生する可能性について判断する処理である。α２の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。

（ステップＳ１４０８）コレスキー分解部２１は、α１×α２を算出する。前述のように、α１はβ１＋１／δ^（−１）＝ｕ_２ｋ＋１ ^（−１）＋１／δ^（−１）に等しく、α２は１＋δ^（−１）ｕ_２ｋ ^（−１）に等しいため、α１×α２に逆ミウラ変換を実行するとｑ_ｋ＋１ ^（−１）に等しくなる。したがって、コレスキー分解部２１は、α１×α２を算出してｑ_ｋ＋１ ^（−１）に設定する。

（ステップＳ１４０９）コレスキー分解部２１は、β１をｕ_２ｋ−１ ^（０）（１＋δ^（０）ｕ_２ｋ−２ ^（０））／α２に設定する。ここで、ｕ_２ｋ−１ ^（０）（１＋δ^（０）ｕ_２ｋ−２ ^（０））の値としては、ステップＳ１２０５，Ｓ１２１４で算出したものを用いる。また、このβ１の式にｓｔｄＬＶｖ変換を実行すると、β１は、ｕ_２ｋ−１ ^（−１）に設定されたことになる。したがって、β１におけるｋの値を１だけデクリメントしたことになる。

（ステップＳ１４１０）コレスキー分解部２１は、β１×β２を算出する。前述のように、β１はｕ_２ｋ−１ ^（−１）に等しく、β２はδ^（−１）ｕ_２ｋ ^（−１）と等しいため、β１×β２に逆ミウラ変換を実行すると、ｅ_ｋ ^（−１）に等しくなる。したがって、コレスキー分解部２１は、β１×β２を算出してｅ_ｋ ^（−１）に設定する。

（ステップＳ１４１１）コレスキー分解部２１は、カウンタｋを１だけデクリメントする。
（ステップＳ１４１２）コレスキー分解部２１は、カウンタｋが０に等しいかどうか判断する。そして、０に等しい場合には、ステップＳ１４１３に進み、そうでない場合には、ステップＳ１４０３に戻る。

（ステップＳ１４１３）コレスキー分解部２１は、β１＋１／δ^（−１）を算出する。これは、ステップＳ１４０３でα１を更新した後に、ステップＳ１４０８でα１×α２を計算することと等しい。この場合、α２は１だからである。したがって、コレスキー分解部２１は、β１＋１／δ^（−１）を算出してｑ_１ ^（−１）に設定する。このようにして、ミウラ変換された結果にｒｄＬＶｖ変換と、逆ミウラ変換とを実行して、ｑ_ｋ ^（−１）、ｅ_ｋ ^（−１）を算出する処理が終了する。これらの値は、コレスキー分解部２１が有する図示しないメモリ等において一時的に記憶されてもよい。

なお、図１５のフローチャートの処理によって算出することができるのは１個の特異値に対応する特異ベクトルであるため、全ての特異値に対応する特異ベクトルを算出する場合には、コレスキー分解部２１は、図１５のフローチャートの処理を特異値の数だけ繰り返すことになる。

メモリ消費を抑えるために、補助変数ｕ_ｋ ^（ｎ）のための配列は必ずしも用意する必要はない。一方、１＋δ^（０）ｕ^（０）のためのメモリ領域を確保し、ミウラ変換、ｄＬＶｖ型変換、逆ミウラ変換のステップにまたがってこの値を利用することで、メモリ消費を抑え、計算量を低減することができる。これにより、誤差も低減される。

ここで、誤差について説明する。人間が理想的な状況、すなわち、無限桁の計算をいくらでもできるとすると、ｑｄ型ツイスト分解法であっても、ＬＶ型ツイスト分解法であっても問題ない。しかしながら、コンピュータで計算を行う場合は注意が必要である。有限桁の計算しか行えないコンピュータ上では、数学的に正しい計算法を使ったとしても必ずしも正しい結果が得られる訳ではない。そればかりか、いつまでたっても計算が終了しないといった思わぬ数値的な問題が発生する場合もある。

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

問題となるのは、同符号の実数の減算及び異符号の実数の加算により生じる、桁落ちである。桁落ちによる誤差で値が０となった後、その値による除算を行うと、０が分母にくるような不定形となり計算不可能となる。こうなるといつまでたっても計算が終了しない。減算→除算と計算する部分がｑｄ型ツイスト分解法、ＬＶ型ツイスト分解法の両方に存在するので、桁落ち誤差には十分に注意する必要がある。

ＬＶ型ツイスト分解法では、上述の桁落ちによる誤差を含んでいるかどうかは減算によって得られた値が小さいかどうかで判断できる。ｑｄ型ツイスト分解法の場合、桁落ち誤差を含むことが分かったとしても、それを回避することはできない。なぜならば、初期値として｛ｑ_ｋ ^（０），ｅ_ｋ ^（０）｝が与えられると、σ_ｊは一意的に決定され、｛ｑ_ｋ ^（±１），ｅ_ｋ ^（±１）｝も一意的に導出されるためである。すなわち、任意パラメータを持たない自由度のない計算法であるためである。

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

ｑｄ型ツイスト分解法あるいはＬＶ型ツイスト分解法によってコレスキー分解をすることができると、上述のようにツイストされた行列Ｎ_ｋを算出することができ、その行列Ｎ_ｋのＮ_ｋ ^Ｔｘ_ｊ＝ｅ_ｋを解くことによって、ｘ_ｊを算出することができる。ここで、

とｘ_ｊを置き換えることにより、このｘ_ｊを正規化する。このようにして、右特異ベクトルｘ_ｊを求めることができる。この右特異ベクトルｘ_ｊを用いて、Ｖ_Ｂ＝（ｘ_１，ｘ_２，...，ｘ_ｍ）とすることにより、右直交行列Ｖ_Ｂを算出することができる。

また、前述のように、

となるため、Ｖ_Ｂが求まれば、２重対角行列Ｂ、特異値を対角成分に有する行列Σは既知であるため、左直交行列を算出することができる。より具体的には、特異値σ_ｊが０でない場合には、

となる。一方、特異値σ_ｊが０である場合には、

を解くことにより、ｙ_ｊを算出することができる。ここで、ｘ_ｊの場合と同様に、

とｙ_ｊを置き換えることにより、このｙ_ｊを正規化する。このようにして、左特異ベクトルｙ_ｊを求めることができる。この左特異ベクトルｙ_ｊを用いて、Ｕ_Ｂ＝（ｙ_１，ｙ_２，...，ｙ_ｍ）とすることにより、左直交行列Ｕ_Ｂを算出することができる。このようにして、２重対角行列Ｂの特異値分解がなされる。なお、この処理の後に、逆反復法や、グラムシュミット法の処理を行ってもよいことは前述の通りである。また、ここでは右直交行列Ｖ_Ｂを先に算出する場合について説明したが、ＢＢ^Ｔの固有ベクトル、すなわち左直交行列Ｕ_Ｂを先に算出してもよい。

また、２重対角行列Ｂの特異値分解を行うことができれば、次のように、行列Ａの特異値分解も行われる。

特異ベクトル算出部１９は、上述のようにして、特異値記憶部１８が記憶している上２重対角行列Ｂの特異値と、対角行列記憶部１３が記憶している上２重対角行列Ｂとを用いて、上２重対角行列Ｂの特異ベクトルを算出する。まず、コレスキー分解部２１は、対角行列記憶部１３から上２重対角行列Ｂを読み出し、特異値記憶部１８から上２重対角行列Ｂの特異値を読み出す。そして、コレスキー分解部２１は、図２２のフローチャートで示されるように各変換を行い、コレスキー分解の処理を行う。ここでは、ＬＶ型ツイスト分解法を用いる場合について説明する。

コレスキー分解部２１は、まず、上２重対角行列Ｂの各要素の値から、ｑ_ｋ ^（０），ｅ_ｋ ^（０）を求める。そして、コレスキー分解部２１は、ミウラ変換を実行することにより、ｕ_１ ^（０），ｕ_２ ^（０）等を順次求めていく（ステップＳ８０１）。

次に、コレスキー分解部２１は、ミウラ変換で得られたｕ_ｋ ^（０）を用いて、ｓｔｄＬＶｖ変換を実行することにより、ｕ_１ ^（＋１）等を順次求めていく（ステップＳ８０２）。また、コレスキー分解部２１は、ミウラ変換で得られたｕ_ｋ ^（０）を用いて、ｒｄＬＶｖ変換を実行することにより、ｕ_２ｍ−１ ^（−１）等を順次求めていく（ステップＳ８０３）。

最後に、コレスキー分解部２１は、ｓｔｄＬＶｖ変換で得られたｕ_ｋ ^（＋１）及びｒｄＬＶｖ変換で得られたｕ_ｋ ^（−１）を用いて、逆ミウラ変換を実行することにより、ｑ_１ ^（±１）、ｅ_１ ^（±１）等を順次求めていく（ステップＳ８０４）。ｑ_ｋ ^（±１）及びｅ_ｋ ^（±１）が求まると、すなわち、上２重対角行列Ｂ^（＋１）と、下２重対角行列Ｂ^（−１）が求まると、コレスキー分解部２１は、それらを第１特異ベクトル算出部２２に渡す。なお、コレスキー分解部２１は、各特異値について、それぞれコレスキー分解を行うものとする。

なお、ここでは、説明の便宜上、図２２のフローチャートを用いてコレスキー分解の処理を説明したが、図１５〜図２０のフローチャートで示されるように処理を行ってもよいことは言うまでもない。

第１特異ベクトル算出部２２は、特異値記憶部１８から特異値を読み出し、式１０を用いてｋの値を決定する。第１特異ベクトル参集部２２は、そのｋの値を用いて、ｑ_ｋ ^（±１）及びｅ_ｋ ^（±１）からツイスト行列Ｎ_ｋを求め、Ｎ_ｋ ^Ｔｘ_ｊ＝ｅ_ｋを解くことによって、ｘ_ｊを算出する（ステップＳ３０２）。

次に、第１特異ベクトル算出部２２は、算出したｘ_ｊを正規化して、右特異ベクトルｘ_ｊを求めて特異ベクトル記憶部２０に蓄積する（ステップＳ３０３）。そして、第１特異ベクトル算出部２２は、算出した右特異ベクトルｘ_ｊを第２特異ベクトル算出部２３に渡す。なお、第１特異ベクトル算出部２２は、各特異値について、ｘ_ｊを算出する処理と、正規化する処理とを行うことによって、全ての右特異ベクトルを算出することができる。

第２特異ベクトル算出部２３は、上２重対角行列Ｂを対角行列記憶部１３から読み出し、その上２重対角行列Ｂの特異値を特異値記憶部１８から読み出す。そして、第２特異ベクトル算出部２３は、第１特異ベクトル算出部２２から受け取った右特異ベクトルｘ_ｊと、上２重対角行列Ｂと、特異値とを用いて、ｙ_ｊ＝Ｂｘ_ｊ／σ_ｊ、あるいは、Ｂ^Ｔｙ_ｊ＝０を解くことによって、ｙ_ｊを算出する（ステップＳ３０４）。次に、第２特異ベクトル算出部２３は、右特異ベクトルｘ_ｊの場合と同様に、算出したｙ_ｊを正規化して、左特異ベクトルｙ_ｊを求めて特異ベクトル記憶部２０に蓄積する（ステップＳ３０５）。このようにして、特異ベクトルの算出が終了する。なお、第２特異ベクトル算出部２３は、各特異値について、ｙ_ｊを算出する処理と、正規化する処理とを行うことによって、全ての左特異ベクトルを算出することができる。
この後、行列記憶部１１が記憶している行列Ａの特異ベクトルを算出してもよい。その算出の方法は、前述の通りである。

また、特異値算出部１７が算出した特異値、特異ベクトル算出部１９が算出した特異ベクトルを出力する図示しない出力部を特異値分解装置１が備えてもよい。ここで、図示しない出力部による出力は、例えば、特異値等の表示デバイス（例えば、ＣＲＴや液晶ディスプレイなど）への表示でもよく、特異値等の所定の機器への通信回線を介した送信でもよく、特異値等のプリンタによる印刷でもよく、特異値等の記録媒体への蓄積でもよい。なお、その出力部は、出力を行うデバイス（例えば、表示デバイスやプリンタなど）を含んでもよく、あるいは含まなくてもよい。また、その出力部は、ハードウェアによって実現されてもよく、あるいは、それらのデバイスを駆動するドライバ等のソフトウェアによって実現されてもよい。

なお、ここでは、コレスキー分解部２１がＬＶ型ツイスト分解法によってコレスキー分解を行う場合について説明したが、コレスキー分解部２１は、ｑｄ型ツイスト分解法によってコレスキー分解を行ってもよい。例えば、コレスキー分解部２１は、算出された特異値の分布を見て、分布が密でない場合には、ｑｄ型ツイスト分解法を用いるようにしてもよい。

また、上記説明では、行列Ｂが上２重対角行列である場合について説明したが、行列Ｂは下２重対角行列であっても、その行列Ｂを上２重対角行列に変換することによって特異値分解を同様に実行することができうる。例えば、行列Ｃが下２重対角行列であるとすると、上２重対角行列Ｂ＝Ｃ^Ｔとすることができる。その上２重対角行列Ｂを、
Ｂ＝ＵΣＶ^Ｔ
と特異値分解できたとすると、
Ｃ＝Ｂ^Ｔ＝（ＵΣＶ^Ｔ）^Ｔ＝ＶΣＵ^Ｔ
となる。したがって、下２重対角行列の特異値分解において、特異値は上２重対角行列Ｂの特異値と同じであり、特異ベクトルは、左右が逆になることがわかる。

一方、下２重対角行列Ｃについても、行列分割部１４等が下２重対角行列Ｃの分割を行い、その分割された下２重対角行列について特異値分解部１５が特異値分解を行い、その特異値分解の結果を用いて特異値算出部１７が下２重対角行列Ｃの特異値を算出するようにしてもよい。また、特異ベクトル算出部１９において、下２重対角行列Ｃと、その特異値とを用いてコレスキー分解を行い、その特異値に対応する各特異ベクトルを算出するようにしてもよい。

また、上記説明では、特異値算出部１７が全ての特異値を算出する場合について説明したが、一部の特異値のみを算出するようにしてもよい。例えば、図１３において、特異値算出部１７は、行列Ｂ_１，Ｂ_２までは、全ての特異値を算出する必要があるが、行列Ｂ_１，Ｂ_２から上２重対角行列Ｂの特異値を算出する際に、必要な範囲で特異値を算出してもよい。したがって、特異値算出部１７は、２重対角行列Ｂの少なくとも１個の特異値を算出するものであってもよい。その場合には、不必要な特異値を算出することに伴う余分な処理を実行しなくてよいため、処理負荷を軽減することができうる。このように、一部の特異値のみを算出する場合には、例えば、上２重対角行列Ｂの特異値を算出する処理に費やす計算時間を約（１＋ｋ／ｍ）／２倍にすることができうる。ここで、ｍは行列サイズであり、ｋは求める特異値の個数である。

また、上記説明では、特異ベクトル算出部１９が算出された特異値に対応する全ての特異ベクトルを算出する場合について説明したが、上記説明から明らかなように、特異ベクトルを算出する処理は、特異値ごとに行うことができる。したがって、特異ベクトル算出部１９は、２重対角行列Ｂの少なくとも１個の特異ベクトルを算出するものであってもよい。このように、特異ベクトル算出部１９は、必要な範囲で特異ベクトルを算出することができ、不必要な特異ベクトルを算出することに伴う余分な処理を実行しなくてよいため、処理負荷を軽減することができうる。

次に、本実施の形態による特異値分解装置１における並列処理について説明する。
本実施の形態による特異値分解装置１では、特異値の計算及び特異ベクトルの計算において、並列的に処理を実行することができる。例えば、図２３で示されるように、特異値分解装置１において、特異値分解部１５、特異値算出部１７、コレスキー分解部２１、第１特異ベクトル算出部２２、第２特異ベクトル算出部２３において、それぞれ並列処理を行ってもよい。

特異値分解部１５は、複数の特異値分解手段１５ａ、１５ｂを備え、その複数の特異値分解手段１５ａ、１５ｂが、分割後の２重対角行列の特異値と、行列要素とを算出する処理を並列実行してもよい。例えば、図１３において、特異値分解手段１５ａが行列Ｂ_１１１，Ｂ_１１２，Ｂ_１２１，Ｂ_１２１の特異値分解を実行し、特異値分解手段１５ｂが行列Ｂ_２１１，Ｂ_２１２，Ｂ_２２１，Ｂ_２２２の特異値分解を実行してもよい。

特異値算出部１７は、複数の特異値算出手段１７ａ、１７ｂを備え、その複数の特異値算出手段１７ａ、１７ｂが、分割元の２重対角行列の特異値と、行列要素とを算出する処理を並列実行してもよい。例えば、図１３において、特異値算出手段１７ａが行列Ｂ_１１，Ｂ_１２，Ｂ_１，Ｂの特異値等を算出する処理を実行し、特異値算出手段１７ｂが行列Ｂ_２１，Ｂ_２２，Ｂ_２の特異値等を算出する処理を実行してもよい。

コレスキー分解部２１は、複数のコレスキー分解手段２１ａ、２１ｂを備え、その複数のコレスキー分解手段２１ａ、２１ｂが、２重対角行列Ｂをコレスキー分解する処理を並列実行してもよい。例えば、コレスキー分解手段２１ａが半分の特異値についてコレスキー分解する処理を実行し、コレスキー分解手段２１ｂが残りの半分の特異値についてコレスキー分解する処理を実行してもよい。

第１特異ベクトル算出部２２は、複数の第１特異ベクトル算出手段２２ａ、２２ｂを備え、その複数の第１特異ベクトル算出手段２２ａ、２２ｂが、特異ベクトルを算出する処理を並列実行してもよい。例えば、第１特異ベクトル算出手段２２ａがコレスキー分解手段２１ａによってコレスキー分解された特異値に関して特異ベクトルを算出する処理を実行し、第１特異ベクトル算出手段２２ｂがコレスキー分解手段２１ｂによってコレスキー分解された特異値に関して特異ベクトルを算出する処理を実行してもよい。

第２特異ベクトル算出部２３は、複数の第２特異ベクトル算出手段２３ａ、２３ｂを備え、その複数の第２特異ベクトル算出手段２３ａ、２３ｂが、特異ベクトルを算出する処理を並列実行してもよい。例えば、第２特異ベクトル算出手段２３ａが第１特異ベクトル算出手段２２ａによって算出された特異ベクトルに対応する特異ベクトルを算出する処理を実行し、第２特異ベクトル算出手段２３ｂが第１特異ベクトル算出手段２２ｂによって算出された特異ベクトルに対応する特異ベクトルを算出する処理を実行してもよい。

なお、特異値算出部１７の複数の特異値算出手段１７ａ、１７ｂは、図１０で示される分割された２個の行列Ｂ_１，Ｂ_２の特異値と行列要素とから、分割元の行列Ｂ_０の特異値と行列要素とを算出する処理を、並列実行してもよい。以下、その処理について説明する。

まず、式２から、各特異値を求める処理については、並列実行可能なことがわかる。ここで、式５を次のように書き換えることにする。

この式１１から、式５を計算するのに、一部の特異値に対応する式１１の部分を計算しておき、他の特異値に対応する式１１の部分を後から掛け合わせることによって、最終的に

を算出できることがわかる。したがって、特異値算出手段１７ａ、１７ｂのそれぞれは、まず、２個の行列Ｂ_１，Ｂ_２の特異値と行列要素とを特異値分解記憶部１６から読み出し、２個の行列Ｂ_１，Ｂ_２の分割時に発生した２個の要素を対角行列記憶部１３から読み出す。その読み出された行列要素と、２個の要素とを用いて、ｚ_ｉの値を知ることができる。その後、各特異値算出手段１７ａ、１７ｂは、式２を用いて、それぞれが担当する特異値を算出する。この処理は、並列実行することができる。

次に、特異値算出手段１７ａは、算出した特異値を用いて計算することができる式１１の部分について計算する。また同様に、特異値算出手段１７ｂも、算出した特異値を用いて計算することができる式１１の部分について計算する。そして、特異値算出手段１７ａ、１７ｂは、その計算した値を交換し、あらかじめ計算していた値と掛け合わせることによって、最終的に式１１の値を計算することができる。このように、特異値算出手段１７ａ、１７ｂは、式１１を計算する処理も並列実行することができる。

次に、特異値算出手段１７ａ、１７ｂは、式４を用いて、それぞれが算出した特異値に対応する特異ベクトルを計算する。このように、特異値算出手段１７ａ、１７ｂは、右直交行列Ｖ_Ｍを計算する処理も並列実行することができる。

その後、特異値算出手段１７ａ、あるいは特異値算出手段１７ｂが、式６から式９を用いて行列Ｂ_０の行列要素を計算することによって、分割元の行列Ｂ_０の特異値と行列要素とを算出する処理が終了する。このように、分割された２個の行列Ｂ_１，Ｂ_２の特異値と行列要素とから、分割元の行列Ｂ_０の特異値と行列要素とを算出する処理を、特異値算出手段１７ａ、１７ｂによって並列実行することができる。

上記のような各構成要素における並列処理において、各手段が特異値等の情報を格納するメモリを用いる場合に、複数の手段が同一のメモリ、すなわち共有メモリを使用してもよく、あるいは、各手段がそれぞれ別のメモリを使用してもよい。

なお、図２３を用いて、特異値分解部１５や特異値算出部１７等が２個の手段によって並列処理を実行する場合について説明したが、特異値分解部１５や特異値算出部１７等が３以上の手段を備え、並列処理を実行してもよい。また、特異値分解部１５、特異値算出部１７、コレスキー分解部２１、第１特異ベクトル算出部２２、第２特異ベクトル算出部２３のそれぞれにおいて並列処理が実行される場合について説明したが、いずれかの任意の１以上の部において、並列処理が実行されなくてもよい。例えば、特異値分解部１５において並列処理が行われなくてもよい。

また、その並列処理は、１の装置において、２以上のＣＰＵ等を用いて実行されてもよく、２以上の装置において実行されてもよい。例えば、図２４で示されるように、装置Ａと、装置Ｂとがそれぞれ特異値分解手段１５ａ、１５ｂを備え、各装置において、特異値分解の処理が並列実行されてもよい。この場合には、装置Ａの特異値分解手段１５ａを備える特異値分解部１５−１と、装置Ｂの特異値分解手段１５ｂを備える特異値分解部１５−２とによって特異値分解部１５が構成されることになる。したがって、特異値分解装置１は、装置Ａと、装置Ｂとからなるシステムを構成することになる。ここでは、特異値分解部１５の並列処理について説明したが、その他の特異値算出部１７やコレスキー分解部２１等についても、２以上の装置による並列処理を行ってもよい。

［２次元画像から３次元へ復元する画像処理への応用］
次に、物体の２次元画像から３次元へ復元する画像処理に特異値分解を応用する場合について説明する。

複数の２次元画像から３次元復元を行う処理は、以下の各ステップによって行われる。
（１）２次元画像から特徴点を抽出するステップ。
（２）特徴点データより形状（元の物体の特徴点の３次元座標データ）及び回転（３次元データから特徴点データへの変換）に関するデータを計算するステップ。
（３）形状及び回転のデータより可視化を行うステップ。

以下、上記（１）、（２）の各ステップについて、図２５のフローチャートを用いて説明する。ここで、２次元画像は、例えば、スキャナやデジタルスチルカメラ、デジタルビデオカメラ等の光学的読み取り機器によって読み取られた２次元画像であってもよい。

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

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

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

そこで、行列Ｄの特異値分解
Ｄ＝ＵΣＶ^Ｔ
を考える。ここで、Σは特異値を大小順に対角線上に並べたもので、Ｕ及びＶはそれぞれ左直交行列、及び右直交行列である。この特異値分解として、前述の方法を用いることができる。すなわち、特異値分解装置１において、行列記憶部１１が記憶している行列Ａを、このたびの行列Ｄとすることにより、前述の説明のようにして、行列Ａの特異値分解がなされることになる。なお、特異値分解装置１において、特異ベクトル記憶部２０に蓄積されるのは、行列Ｄを変換した上２重対角行列Ｂの特異ベクトルであるので、その特異ベクトルを、前述の説明のようにして、行列Ｄの特異ベクトルに変換する必要がある。

ここで、画像のデジタル誤差のため、ゼロでない特異値は３つ以上出てくる。しかし、４番目以降の特異値はノイズによるもので、最初の３つの特異値と比べて格段に小さい。そこで、ステップＳ５００２では、最初の３つの特異値に対して特異ベクトルを計算する。採用する３個のベクトルをまとめると、次式となる。
Ｄ'＝Ｌ'Σ'Ｒ'^Ｔ＝Ｍ'Ｓ'
ただし、Ｍ'＝Ｌ'（Σ'）^１／２、Ｓ'＝（Σ'）^１／２Ｒ'^Ｔ、Ｄ'は‖Ｄ−Ｄ'‖を最小にするランク３の行列である。

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

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

次に、ステップＳ５００４において、ＥからＣを求める。Ｃの自由度（９）はＥの自由度（６）より多い。そこで、条件ｒ_１ｊ＝（１，０，０）^Ｔ，ｒ_２ｊ＝（０，１，０）^Ｔを加えれば、Ｃを決めることができる。このとき２つの解（Ｎｅｃｋｅｒ Ｒｅｖｅｒｓａｌ）が出る。
次に、ステップＳ５００５において、Ｍ＝Ｍ'Ｃ及びＳ＝Ｃ^−１Ｓ'より、回転に関するデータＭ及び形状Ｓが決まる。

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

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

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

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

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

図２６は、本実施の形態による特異値分解装置１を利用した文書検索方法の一例を示すフローチャートである。
ステップ６００１において、質問ベクトルｑを受け取る。

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

ステップＳ６００２において、ｋを選択する。ここでは、ｋ<ｒなるｋの値を選択する。ｒ＝ｍｉｎ（ｎ，ｍ）である。ｋの値は、予め与えられていてもよいし、計算ごとに選択可能であってもよい。

次に、ステップＳ６００３では、行列Ｄの特異値分解を行う。この特異値分解として、前述の方法を用いることができる。すなわち、特異値分解装置１において、行列記憶部１１が記憶している行列Ａを、このたびの行列Ｄとすることにより、前述の説明のようにして、行列Ａの特異値分解がなされることになる。なお、特異値分解装置１において、特異ベクトル記憶部２０に蓄積されるのは、行列Ｄを変換した上２重対角行列Ｂの特異ベクトルであるので、その特異ベクトルを、前述の説明のようにして、行列Ｄの特異ベクトルに変換する必要がある。また、この特異値分解では、計算された特異値のうち、大きい順に１番目からｋ番目までのｋ個の特異値に対してＤの特異ベクトルを算出する。ｋは、ステップＳ６００２で選択された値である。

すなわち、
Ｄ_ｋ＝Ｕ_ｋΣＶ_ｋ ^Ｔ
なるＵ_ｋ及びＶ_ｋを計算する。ここで、Ｕ_ｋは、最初のｋ個の左特異ベクトルのみから構成されるｍ×ｋ行列であり、Ｖ_ｋは、最初のｋ個の右特異ベクトルのみから構成されるｎ×ｋ行列であり、Σは、最初のｋ個の特異値のみから構成されるｋ×ｋ対角行列である。

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

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

ステップＳ６００５において、ステップＳ６００４において計算された類似度を基準に、検索結果を出力する。
最後に、数値実験による性能の評価について説明する。ここでは、本実施の形態による特異値分解装置１の方法と、標準的な分割統治法と、シフト付きＱＲ法とを比較した。分割統治法としては、ＬＡＰＡＣＫで提供されているＤＢＤＳＤＣを用いた。ＱＲ法としては、ＬＡＰＡＣＫで提供されているＤＢＤＳＱＲを用いた。実験で用いるＬＡＰＡＣＫはＡＴＬＡＳ（Ａｕｔｏｍａｔｉｃａｌｌｙ Ｔｕｎｅｄ Ｌｉｎｅａｒ Ａｌｇｅｂｒａ Ｓｏｆｔｗａｒｅ）で最適化したＢＬＡＳ（Ｂａｓｉｃ Ｌｉｎｅａｒ Ａｌｇｅｂｒａ Ｓｕｂｐｒｏｇｒａｍｓ）を呼び出す。また、本実施の形態による特異値分解装置１の特異値分解部１５では、ＱＲ法を用いて特異値分解を行った。計算機は、キャッシュメモリＬ１Ｄ：６４ＫＢ，Ｌ１Ｉ：６４ＫＢ，Ｌ２：１０２４ＫＢを搭載したＯｐｔｅｒｏｎ １．８ＧＨｚを用いた。

図２７は、計算時間を示す図である。特異値分解装置１による特異値分解が、他の特異値分解方法よりも非常に高速であることがわかる。特に、分割統治法（Ｄ&Ｃ）と比較すると、ベクトル更新を省略した効果がわかる。

図２８は、異なる行列サイズの計算時間の変化を示す図である。特異値分解装置１による特異値分解は、ほぼＯ（ｎ^２）という計算量を持つことがわかる。
なお、図２７及び図２８で示される計算時間の測定では、対角成分が２．００１，非対角成分が２．０であり、近接特異値を持たない行列について計算を行った。

図２９は、計算制度を示す図である。特異値分解装置１による特異値の算出は、分割統治法と同程度の高い精度を有することがわかる。図２９では、１００個のランダム行列を解き、精度を評価した。行列サイズは、１０００である。

以上のように、本実施の形態による特異値分解装置１では、分割統治法を用いて特異値のみを算出するため、標準的な分割統治法で実行されるベクトル更新が必要なくなり、標準的な分割統治法よりも非常に高速である。また、特異値から特異ベクトルを算出する処理も、高速に処理することができる。さらに、ＱＲ法では特異値を算出する段階の並列化が困難であるが、本実施の形態による特異値分解装置１では、特異値を算出する段階、及び特異値から特異ベクトルを算出する段階のそれぞれにおいて、本質的に高い並列性を持つ。また、本実施の形態による特異値分解装置１の精度は、分割統治法やＱＲ法とほぼ同等の精度が得られることがわかる。

なお、本実施の形態による特異値分解装置１は、スタンドアロンの装置であってもよく、サーバクライアントシステムにおけるサーバ装置であってもよく、前述のように、複数の装置から構成されるシステムであってもよい。特異値分解装置１がサーバクライアントシステムにおけるサーバ装置である場合には、行列記憶部１１で記憶される行列Ａや、特異値記憶部１８で記憶される特異値や特異ベクトル記憶部２０で記憶される特異ベクトル等は、インターネットやイントラネット等の通信回線を介して送受信されてもよい。

また、特異値分解装置１における一部の処理を、特異値分解装置１以外において行ってもよい。例えば、行列分解部１４による２重対角行列Ｂを分割する処理、あるいは、特異値分解部１５による特異値分解の処理を、特異値分解装置１において実行しなくてもよい。その場合には、例えば、２重対角行列Ｂが分割された結果である、あらかじめ決められた大きさ以下の各２重対角行列が、入力デバイスや通信回線、記録媒体等を介して受け付けられ、対角行列記憶部１３で記憶されていてもよい。また、例えば、そのあらかじめ決められた大きさ以下の各２重対角行列に対して特異値分解が行われた結果である、各２重対角行列の特異値と、各２重対角行列の特異ベクトルからなる左右直交行列の一部の要素である行列要素とが、入力デバイスや通信回線、記録媒体等を介して受け付けられ、特異値分解記憶部１５で記憶されていてもよい。

また、上記実施の形態において、各処理または各機能は、単一の装置または単一のシステムによって集中処理されることによって実現されてもよく、あるいは、複数の装置または複数のシステムによって分散処理されることによって実現されてもよい。

また、上記実施の形態において、各構成要素は専用のハードウェアにより構成されてもよく、あるいは、ソフトウェアにより実現可能な構成要素については、プログラムを実行することによって実現されてもよい。例えば、ハードディスクや半導体メモリ等の記録媒体に記録されたソフトウェア・プログラムをＣＰＵ等のプログラム実行部が読み出して実行することによって、各構成要素が実現され得る。なお、上記実施の形態における特異値分解装置を実現するソフトウェアは、以下のようなプログラムである。つまり、このプログラムは、コンピュータに、２重対角行列Ｂが２個の２重対角行列に分割され、その２重対角行列が２個の２重対角行列に分割される処理が、分割後の各２重対角行列があらかじめ決められた大きさ以下となるまで繰り返され、当該あらかじめ決められた大きさ以下の各２重対角行列に対して特異値分解が行われた結果である、前記各２重対角行列の特異値と、前記各２重対角行列の特異ベクトルからなる左右直交行列の一部の要素である行列要素とが記憶される特異値分解記憶部から、各２重対角行列の特異値と、行列要素とを読み出し、前記特異値と、前記行列要素とから、分割元の２重対角行列の特異値と、分割元の２重対角行列の行列要素とを算出して前記特異値分解記憶部に蓄積し、その分割元の２重対角行列の特異値と、行列要素とを算出する処理を、２重対角行列Ｂの少なくとも１個の特異値を算出するまで繰り返し、前記２重対角行列Ｂの少なくとも１個の特異値を特異値記憶部に蓄積する特異値算出ステップと、前記２重対角行列Ｂが記憶される対角行列記憶部から２重対角行列Ｂを読み出し、前記特異値記憶部から前記２重対角行列Ｂの特異値を読み出し、２重対角行列Ｂとその特異値とから、ツイスト分解法を用いて２重対角行列Ｂの少なくとも１個の特異ベクトルを算出して特異ベクトル記憶部に蓄積する特異ベクトル算出ステップと、を実現させるためのものである。

また、このプログラムでは、コンピュータに、前記あらかじめ決められた大きさ以下の各２重対角行列も記憶される前記対角行列記憶部から、前記あらかじめ決められた大きさ以下の各２重対角行列を読み出し、前記各２重対角行列に対して特異値分解を行い、前記各２重対角行列の特異値と、前記各２重対角行列の特異ベクトルとを算出し、当該特異値と、当該特異ベクトルからなる左右直交行列の一部の要素である行列要素とを前記特異値分解記憶部に蓄積する特異値分解ステップをさらに実行させてもよい。

また、上記実施の形態における特異値分解装置を実現するソフトウェアは、以下のようなプログラムである。つまり、このプログラムは、コンピュータに、２重対角行列Ｂが記憶される対角行列記憶部から前記２重対角行列Ｂを読み出し、当該２重対角行列Ｂを２個の２重対角行列に分割して前記対角行列記憶部に蓄積し、その２重対角行列を２個の２重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各２重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割ステップと、前記あらかじめ決められた大きさ以下の各２重対角行列を前記対角行列記憶部から読み出し、前記各２重対角行列に対して特異値分解を行い、前記各２重対角行列の特異値と、前記各２重対角行列の特異ベクトルとを算出し、特異値分解された特異値と、特異ベクトルからなる左右直交行列の一部の要素である行列要素とを特異値分解記憶部に蓄積する特異値分解ステップと、各２重対角行列の特異値と、行列要素とを前記特異値分解記憶部から読み出し、前記特異値と、前記行列要素とから、分割元の２重対角行列の特異値と、分割元の２重対角行列の行列要素とを算出して前記特異値分解記憶部に蓄積し、その分割元の２重対角行列の特異値と、行列要素とを算出する処理を、２重対角行列Ｂの少なくとも１個の特異値を算出するまで繰り返し、前記２重対角行列Ｂの少なくとも１個の特異値を特異値記憶部に蓄積する特異値算出ステップと、前記対角行列記憶部から２重対角行列Ｂを読み出し、前記特異値記憶部から前記２重対角行列Ｂの特異値を読み出し、２重対角行列Ｂとその特異値とから、ツイスト分解法を用いて２重対角行列Ｂの少なくとも１個の特異ベクトルを算出して特異ベクトル記憶部に蓄積する特異ベクトル算出ステップと、を実行させるためのものである。

また、このプログラムでは、前記特異ベクトル算出ステップは、前記対角行列記憶部から２重対角行列Ｂを読み出し、前記特異値記憶部から前記２重対角行列Ｂの特異値を読み出し、前記２重対角行列Ｂの各要素に関してミウラ変換、ｄＬＶｖ型変換、逆ミウラ変換を行うことによって、前記２重対角行列Ｂを上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）にコレスキー分解するコレスキー分解ステップと、前記上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）の各要素と、前記２重対角行列Ｂの特異値とを用いて一方の左右直交行列を構成する特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する第１特異ベクトル算出ステップと、前記第１特異ベクトル算出ステップで算出した一方の左右直交行列を構成する特異ベクトルと、前記２重対角行列Ｂの特異値と、前記２重対角行列Ｂとを用いて他方の左右直交行列を構成する特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する第２特異ベクトル算出ステップと、をさらに備えてもよい。
なお、上記プログラムにおいて、情報を蓄積する処理などでは、ハードウェアでしか行われない処理は少なくとも含まれない。

また、このプログラムは、サーバなどからダウンロードされることによって実行されてもよく、所定の記録媒体（例えば、ＣＤ−ＲＯＭなどの光ディスクや磁気ディスク、半導体メモリなど）に記録されたプログラムが読み出されることによって実行されてもよい。

また、このプログラムを実行するコンピュータは、単数であってもよく、複数であってもよい。すなわち、集中処理を行ってもよく、あるいは分散処理を行ってもよい。

図３０は、上記プログラムを実行して、上記実施の形態による特異値分解装置１を実現するコンピュータの外観の一例を示す模式図である。上記実施の形態は、コンピュータハードウェア及びその上で実行されるコンピュータプログラムによって実現されうる。

図３０において、コンピュータシステム１００は、ＣＤ−ＲＯＭ（Ｃｏｍｐａｃｔ Ｄｉｓｋ Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）ドライブ１０５、ＦＤ（Ｆｌｅｘｉｂｌｅ Ｄｉｓｋ）ドライブ１０６を含むコンピュータ１０１と、キーボード１０２と、マウス１０３と、モニタ１０４とを備える。

図３１は、コンピュータシステムを示す図である。図３１において、コンピュータ１０１は、ＣＤ−ＲＯＭドライブ１０５、ＦＤドライブ１０６に加えて、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）１１１と、ブートアッププログラム等のプログラムを記憶するためのＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）１１２と、ＣＰＵ１１１に接続され、アプリケーションプログラムの命令を一時的に記憶すると共に、一時記憶空間を提供するＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）１１３と、アプリケーションプログラム、システムプログラム、及びデータを記憶するハードディスク１１４と、ＣＰＵ１１１、ＲＯＭ１１２等を相互に接続するバス１１５とを備える。なお、コンピュータ１０１は、ＬＡＮへの接続を提供する図示しないネットワークカードを含んでもよい。

コンピュータシステム１００に、上記実施の形態による特異値分解装置１の機能を実行させるプログラムは、ＣＤ−ＲＯＭ１２１、またはＦＤ１２２に記憶されて、ＣＤ−ＲＯＭドライブ１０５、またはＦＤドライブ１０６に挿入され、ハードディスク１１４に転送されてもよい。これに代えて、そのプログラムは、図示しないネットワークを介してコンピュータ１０１に送信され、ハードディスク１１４に記憶されてもよい。プログラムは実行の際にＲＡＭ１１３にロードされる。なお、プログラムは、ＣＤ−ＲＯＭ１２１やＦＤ１２２、またはネットワークから直接、ロードされてもよい。

また、行列記憶部１１、対角行列記憶部１３、特異値分解記憶部１６、特異値記憶部１８，特異ベクトル記憶部２０は、ＲＡＭ１１３やハードディスク１１４によって実現されてもよい。

プログラムは、コンピュータ１０１に、上記実施の形態による特異値分解装置１の機能を実行させるオペレーティングシステム（ＯＳ）、またはサードパーティプログラム等を必ずしも含まなくてもよい。プログラムは、制御された態様で適切な機能（モジュール）を呼び出し、所望の結果が得られるようにする命令の部分のみを含んでいてもよい。コンピュータシステム１００がどのように動作するのかについては周知であり、詳細な説明を省略する。

また、本発明は、以上の実施の形態に限定されることなく、種々の変更が可能であり、それらも本発明の範囲内に包含されるものであることは言うまでもない。
また、本発明のほんのいくつかの典型的な実施例について上で詳細に説明したが、その典型的な実施例において、発明の利益と新規な技術から実質的にはずれることなく多くの変更が可能であることを当業者は容易に認識することができるであろう。したがって、そのようなすべての変更は、本発明の範囲に含まれるものである。

以上のように、本発明による特異値分解装置等によれば、特異値分解を高速に処理することができ、画像処理や検索処理、その他の特異値分解を用いる処理を実行する装置等において有用である。

本発明の実施の形態１による特異値分解装置の構成を示すブロック図 同実施の形態による特異値分解装置の動作を示すフローチャート 同実施の形態による特異値分解装置の動作を示すフローチャート 同実施の形態による特異値分解装置の動作を示すフローチャート 同実施の形態による特異値分解装置の動作を示すフローチャート 同実施の形態における上２重対角行列Ｂの分割について説明するための図 同実施の形態における上２重対角行列Ｂの分割について説明するための図 同実施の形態による特異値分解装置の動作を示すフローチャート 同実施の形態における上２重対角行列Ｂの分割について説明するための図 同実施の形態における上２重対角行列Ｂの特異値の算出について説明するための図 同実施の形態による特異値分解装置の動作を示すフローチャート 同実施の形態による特異値分解装置の動作を示すフローチャート 同実施の形態における上２重対角行列Ｂの特異値の算出について説明するための図 同実施の形態におけるコレスキー分解について説明するための図 同実施の形態による特異値分解装置の動作を示すフローチャート 同実施の形態による特異値分解装置の動作を示すフローチャート 同実施の形態による特異値分解装置の動作を示すフローチャート 同実施の形態による特異値分解装置の動作を示すフローチャート 同実施の形態による特異値分解装置の動作を示すフローチャート 同実施の形態による特異値分解装置の動作を示すフローチャート 同実施の形態におけるｑｄ型ツイスト分解法について説明するための図 同実施の形態におけるＬＶ型ツイスト分解法について説明するための図 同実施の形態による特異値分解装置の動作を示すフローチャート 同実施の形態による特異値分解装置の構成の他の一例を示すブロック図 同実施の形態による特異値分解装置の構成の他の一例を示すブロック図 同実施の形態における画像処理の一例を示すフローチャート 同実施の形態における文書検索処理の一例を示すフローチャート 同実施の形態による特異値分解と従来の特異値分解との比較について説明するための図 同実施の形態による特異値分解と従来の特異値分解との比較について説明するための図 同実施の形態による特異値分解と従来の特異値分解との比較について説明するための図 同実施の形態におけるコンピュータシステムの外観一例を示す模式図 同実施の形態におけるコンピュータシステムの構成の一例を示す図

Claims (23)

1. ２重対角行列Ｂが記憶される対角行列記憶部と、
   前記２重対角行列Ｂが２個の２重対角行列に分割され、その２重対角行列が２個の２重対角行列に分割される処理が、分割後の各２重対角行列があらかじめ決められた大きさ以下となるまで繰り返され、当該あらかじめ決められた大きさ以下の各２重対角行列に対して特異値分解が行われた結果である、前記各２重対角行列の特異値と、前記各２重対角行列の特異ベクトルからなる左右直交行列の一部の要素である行列要素とが記憶される特異値分解記憶部と、
   前記２重対角行列Ｂの特異値が記憶される特異値記憶部と、
   各２重対角行列の特異値と、行列要素とを前記特異値分解記憶部から読み出し、前記特異値と、前記行列要素とから、分割元の２重対角行列の特異値と、分割元の２重対角行列の行列要素とを算出して前記特異値分解記憶部に蓄積し、その分割元の２重対角行列の特異値と、行列要素とを算出する処理を、２重対角行列Ｂの少なくとも１個の特異値を算出するまで繰り返し、前記２重対角行列Ｂの少なくとも１個の特異値を前記特異値記憶部に蓄積する特異値算出部と、
   前記２重対角行列Ｂの特異ベクトルが記憶される特異ベクトル記憶部と、
   前記対角行列記憶部から２重対角行列Ｂを読み出し、前記特異値記憶部から前記２重対角行列Ｂの特異値を読み出し、２重対角行列Ｂとその特異値とから、ツイスト分解法を用いて２重対角行列Ｂの少なくとも１個の特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する特異ベクトル算出部と、を備えた特異値分解装置。
2. 前記対角行列記憶部には、前記あらかじめ決められた大きさ以下の各２重対角行列も記憶され、
   前記あらかじめ決められた大きさ以下の各２重対角行列を前記対角行列記憶部から読み出し、前記各２重対角行列に対して特異値分解を行い、前記各２重対角行列の特異値と、前記各２重対角行列の特異ベクトルとを算出し、当該特異値と、当該特異ベクトルからなる左右直交行列の一部の要素である行列要素とを前記特異値分解記憶部に蓄積する特異値分解部をさらに備えた請求項１記載の特異値分解装置。
3. ２重対角行列Ｂが記憶される対角行列記憶部と、
   前記対角行列記憶部から前記２重対角行列Ｂを読み出し、当該２重対角行列Ｂを２個の２重対角行列に分割して前記対角行列記憶部に蓄積し、その２重対角行列を２個の２重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各２重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割部と、
   前記あらかじめ決められた大きさ以下の各２重対角行列を前記対角行列記憶部から読み出し、前記各２重対角行列に対して特異値分解を行い、前記各２重対角行列の特異値と、前記各２重対角行列の特異ベクトルとを算出する特異値分解部と、
   前記特異値分解部によって特異値分解された特異値と、特異ベクトルからなる左右直交行列の一部の要素である行列要素とが記憶される特異値分解記憶部と、
   前記２重対角行列Ｂの特異値が記憶される特異値記憶部と、
   各２重対角行列の特異値と、行列要素とを前記特異値分解記憶部から読み出し、前記特異値と、前記行列要素とから、分割元の２重対角行列の特異値と、分割元の２重対角行列の行列要素とを算出して前記特異値分解記憶部に蓄積し、その分割元の２重対角行列の特異値と、行列要素とを算出する処理を、２重対角行列Ｂの少なくとも１個の特異値を算出するまで繰り返し、前記２重対角行列Ｂの少なくとも１個の特異値を前記特異値記憶部に蓄積する特異値算出部と、
   前記２重対角行列Ｂの特異ベクトルが記憶される特異ベクトル記憶部と、
   前記対角行列記憶部から２重対角行列Ｂを読み出し、前記特異値記憶部から前記２重対角行列Ｂの特異値を読み出し、２重対角行列Ｂとその特異値とから、ツイスト分解法を用いて２重対角行列Ｂの少なくとも１個の特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する特異ベクトル算出部と、を備えた特異値分解装置。
4. 前記特異ベクトル算出部は、ｑｄ型ツイスト分解法により特異ベクトルを算出する、請求項３記載の特異値分解装置。
5. 前記特異ベクトル算出部は、ＬＶ型ツイスト分解法により特異ベクトルを算出する、請求項３記載の特異値分解装置。
6. 前記特異ベクトル算出部は、
   前記対角行列記憶部から２重対角行列Ｂを読み出し、前記特異値記憶部から前記２重対角行列Ｂの特異値を読み出し、前記２重対角行列Ｂの各要素に関してミウラ変換、ｄＬＶｖ型変換、逆ミウラ変換を行うことによって、前記２重対角行列Ｂを上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）にコレスキー分解するコレスキー分解部と、
   前記上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）の各要素と、前記２重対角行列Ｂの特異値とを用いて一方の左右直交行列を構成する特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する第１特異ベクトル算出部と、
   前記第１特異ベクトル算出部が算出した一方の左右直交行列を構成する特異ベクトルと、前記２重対角行列Ｂの特異値と、前記２重対角行列Ｂとを用いて他方の左右直交行列を構成する特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する第２特異ベクトル算出部と、をさらに備えた、請求項５記載の特異値分解装置。
7. 前記コレスキー分解部は、複数のコレスキー分解手段を備え、
   前記複数のコレスキー分解手段が、前記２重対角行列Ｂをコレスキー分解する処理を並列実行する、請求項６記載の特異値分解装置。
8. 前記第１特異ベクトル算出部は、複数の第１特異ベクトル算出手段を備え、
   前記複数の第１特異ベクトル算出手段が、特異ベクトルを算出する処理を並列実行する、請求項６または請求項７記載の特異値分解装置。
9. 前記第２特異ベクトル算出部は、複数の第２特異ベクトル算出手段を備え、
   前記複数の第２特異ベクトル算出手段が、特異ベクトルを算出する処理を並列実行する、請求項６から請求項８のいずれか記載の特異値分解装置。
10. 前記特異値算出部は、複数の特異値算出手段を備え、
    前記複数の特異値算出手段が、分割元の２重対角行列の特異値と、行列要素とを算出する処理を並列実行する、請求項３から請求項９のいずれか記載の特異値分解装置。
11. 前記特異値分解部は、複数の特異値分解手段を備え、
    前記複数の特異値分解手段が、２重対角行列に対して特異値分解を行う処理を並列実行する、請求項３から請求項１０のいずれか記載の特異値分解装置。
12. 行列Ａが記憶される行列記憶部と、
    前記行列Ａを前記行列記憶部から読み出し、前記行列Ａを２重対角化した前記２重対角行列Ｂを算出して前記対角行列記憶部に蓄積する対角化部と、をさらに備えた請求項３から請求項１１のいずれか特異値分解装置。
13. 前記行列分割部は、２重対角行列を略半分の２個の２重対角行列に分割する、請求項３から請求項１２のいずれか記載の特異値分解装置。
14. 前記特異値算出部は、２重対角行列Ｂの全ての特異値を算出する、請求項３から請求項１３のいずれか記載の特異値分解装置。
15. 前記特異ベクトル算出部は、２重対角行列Ｂの全ての特異ベクトルを算出する、請求項１４記載の特異値分解装置。
16. 前記行列Ａは、２次元画像ｊ（ｊ＝１，…，ｍ、ｍは３以上の整数）から抽出された特徴点ｉ（ｉ＝１，…，ｎ、ｎは２以上の整数）の座標（ｘ_ｉ ^ｊ，ｙ_ｉ ^ｊ）から構成される行列である、請求項１２記載の特異値分解装置。
17. 前記行列Ａは、索引語の重みを要素とするベクトルであって、検索対象となる文書を示すベクトルであるｄ_１，…，ｄ_ｎ（ｎは２以上の整数）を各列に有する行列である索引語文書行列である、請求項１２記載の特異値分解装置。
18. ２重対角行列Ｂが記憶される対角行列記憶部と、前記２重対角行列Ｂが２個の２重対角行列に分割され、その２重対角行列が２個の２重対角行列に分割される処理が、分割後の各２重対角行列があらかじめ決められた大きさ以下となるまで繰り返され、当該あらかじめ決められた大きさ以下の各２重対角行列に対して特異値分解が行われた結果である、前記各２重対角行列の特異値と、前記各２重対角行列の特異ベクトルからなる左右直交行列の一部の要素である行列要素とが記憶される特異値分解記憶部と、前記２重対角行列Ｂの特異値が記憶される特異値記憶部と、特異値算出部と、前記２重対角行列Ｂの特異ベクトルが記憶される特異ベクトル記憶部と、特異ベクトル算出部とを備えた特異値分解装置で用いられる特異値分解方法であって、
    前記特異値算出部が、各２重対角行列の特異値と、行列要素とを前記特異値分解記憶部から読み出し、前記特異値と、前記行列要素とから、分割元の２重対角行列の特異値と、分割元の２重対角行列の行列要素とを算出して前記特異値分解記憶部に蓄積し、その分割元の２重対角行列の特異値と、行列要素とを算出する処理を、２重対角行列Ｂの少なくとも１個の特異値を算出するまで繰り返し、前記２重対角行列Ｂの少なくとも１個の特異値を前記特異値記憶部に蓄積する特異値算出ステップと、
    前記特異ベクトル算出部が、前記対角行列記憶部から２重対角行列Ｂを読み出し、前記特異値記憶部から前記２重対角行列Ｂの特異値を読み出し、２重対角行列Ｂとその特異値とから、ツイスト分解法を用いて２重対角行列Ｂの少なくとも１個の特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する特異ベクトル算出ステップと、を備えた特異値分解方法。
19. ２重対角行列Ｂが記憶される対角行列記憶部と、行列分割部と、特異値分解部と、前記特異値分解部によって特異値分解された特異値と、特異ベクトルからなる左右直交行列の一部の要素である行列要素とが記憶される特異値分解記憶部と、前記２重対角行列Ｂの特異値が記憶される特異値記憶部と、特異値算出部と、前記２重対角行列Ｂの特異ベクトルが記憶される特異ベクトル記憶部と、特異ベクトル算出部と、を備えた特異値分解装置で用いられる特異値分解方法であって、
    前記行列分割部が、前記対角行列記憶部から前記２重対角行列Ｂを読み出し、当該２重対角行列Ｂを２個の２重対角行列に分割して前記対角行列記憶部に蓄積し、その２重対角行列を２個の２重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各２重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割ステップと、
    前記特異値分解部が、前記あらかじめ決められた大きさ以下の各２重対角行列を前記対角行列記憶部から読み出し、前記各２重対角行列に対して特異値分解を行い、前記各２重対角行列の特異値と、前記各２重対角行列の特異ベクトルとを算出する特異値分解ステップと、
    前記特異値算出部が、各２重対角行列の特異値と、行列要素とを前記特異値分解記憶部から読み出し、前記特異値と、前記行列要素とから、分割元の２重対角行列の特異値と、分割元の２重対角行列の行列要素とを算出して前記特異値分解記憶部に蓄積し、その分割元の２重対角行列の特異値と、行列要素とを算出する処理を、２重対角行列Ｂの少なくとも１個の特異値を算出するまで繰り返し、前記２重対角行列Ｂの少なくとも１個の特異値を前記特異値記憶部に蓄積する特異値算出ステップと、
    前記特異ベクトル算出部が、前記対角行列記憶部から２重対角行列Ｂを読み出し、前記特異値記憶部から前記２重対角行列Ｂの特異値を読み出し、２重対角行列Ｂとその特異値とから、ツイスト分解法を用いて２重対角行列Ｂの少なくとも１個の特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する特異ベクトル算出ステップと、を備えた特異値分解方法。
20. 前記特異ベクトル算出部は、コレスキー分解部と、第１特異ベクトル算出部と、第２特異ベクトル算出部と、をさらに備えており、
    前記特異ベクトル算出ステップは、
    前記コレスキー分解部が、前記対角行列記憶部から２重対角行列Ｂを読み出し、前記特異値記憶部から前記２重対角行列Ｂの特異値を読み出し、前記２重対角行列Ｂの各要素に関してミウラ変換、ｄＬＶｖ型変換、逆ミウラ変換を行うことによって、前記２重対角行列Ｂを上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）にコレスキー分解するコレスキー分解ステップと、
    前記第１特異ベクトル算出部が、前記上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）の各要素と、前記２重対角行列Ｂの特異値とを用いて一方の左右直交行列を構成する特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する第１特異ベクトル算出ステップと、
    前記第２特異ベクトル算出部が、前記第１特異ベクトル算出ステップで算出した一方の左右直交行列を構成する特異ベクトルと、前記２重対角行列Ｂの特異値と、前記２重対角行列Ｂとを用いて他方の左右直交行列を構成する特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する第２特異ベクトル算出ステップと、をさらに備えた、請求項１９記載の特異値分解方法。
21. コンピュータに、
    ２重対角行列Ｂが２個の２重対角行列に分割され、その２重対角行列が２個の２重対角行列に分割される処理が、分割後の各２重対角行列があらかじめ決められた大きさ以下となるまで繰り返され、当該あらかじめ決められた大きさ以下の各２重対角行列に対して特異値分解が行われた結果である、前記各２重対角行列の特異値と、前記各２重対角行列の特異ベクトルからなる左右直交行列の一部の要素である行列要素とが記憶される特異値分解記憶部から、各２重対角行列の特異値と、行列要素とを読み出し、前記特異値と、前記行列要素とから、分割元の２重対角行列の特異値と、分割元の２重対角行列の行列要素とを算出して前記特異値分解記憶部に蓄積し、その分割元の２重対角行列の特異値と、行列要素とを算出する処理を、２重対角行列Ｂの少なくとも１個の特異値を算出するまで繰り返し、前記２重対角行列Ｂの少なくとも１個の特異値を特異値記憶部に蓄積する特異値算出ステップと、
    前記２重対角行列Ｂが記憶される対角行列記憶部から２重対角行列Ｂを読み出し、前記特異値記憶部から前記２重対角行列Ｂの特異値を読み出し、２重対角行列Ｂとその特異値とから、ツイスト分解法を用いて２重対角行列Ｂの少なくとも１個の特異ベクトルを算出して特異ベクトル記憶部に蓄積する特異ベクトル算出ステップと、を実現させるためのプログラム。
22. コンピュータに、
    ２重対角行列Ｂが記憶される対角行列記憶部から前記２重対角行列Ｂを読み出し、当該２重対角行列Ｂを２個の２重対角行列に分割して前記対角行列記憶部に蓄積し、その２重対角行列を２個の２重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各２重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割ステップと、
    前記あらかじめ決められた大きさ以下の各２重対角行列を前記対角行列記憶部から読み出し、前記各２重対角行列に対して特異値分解を行い、前記各２重対角行列の特異値と、前記各２重対角行列の特異ベクトルとを算出し、特異値分解された特異値と、特異ベクトルからなる左右直交行列の一部の要素である行列要素とを特異値分解記憶部に蓄積する特異値分解ステップと、
    各２重対角行列の特異値と、行列要素とを前記特異値分解記憶部から読み出し、前記特異値と、前記行列要素とから、分割元の２重対角行列の特異値と、分割元の２重対角行列の行列要素とを算出して前記特異値分解記憶部に蓄積し、その分割元の２重対角行列の特異値と、行列要素とを算出する処理を、２重対角行列Ｂの少なくとも１個の特異値を算出するまで繰り返し、前記２重対角行列Ｂの少なくとも１個の特異値を特異値記憶部に蓄積する特異値算出ステップと、
    前記対角行列記憶部から２重対角行列Ｂを読み出し、前記特異値記憶部から前記２重対角行列Ｂの特異値を読み出し、２重対角行列Ｂとその特異値とから、ツイスト分解法を用いて２重対角行列Ｂの少なくとも１個の特異ベクトルを算出して特異ベクトル記憶部に蓄積する特異ベクトル算出ステップと、を実行させるためのプログラム。
23. 前記特異ベクトル算出ステップは、
    前記対角行列記憶部から２重対角行列Ｂを読み出し、前記特異値記憶部から前記２重対角行列Ｂの特異値を読み出し、前記２重対角行列Ｂの各要素に関してミウラ変換、ｄＬＶｖ型変換、逆ミウラ変換を行うことによって、前記２重対角行列Ｂを上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）にコレスキー分解するコレスキー分解ステップと、
    前記上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）の各要素と、前記２重対角行列Ｂの特異値とを用いて一方の左右直交行列を構成する特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する第１特異ベクトル算出ステップと、
    前記第１特異ベクトル算出ステップで算出した一方の左右直交行列を構成する特異ベクトルと、前記２重対角行列Ｂの特異値と、前記２重対角行列Ｂとを用いて他方の左右直交行列を構成する特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する第２特異ベクトル算出ステップと、をさらに備えた、請求項２２記載のプログラム。

Images