JP5017666B2

JP5017666B2 - 固有値分解装置、及び固有値分解方法

Info

Publication number: JP5017666B2
Application number: JP2008528725A
Authority: JP
Inventors: 佳正中村; 洋明坪井; 太朗譽田; 雅史岩▲崎▼; 雅美高田
Original assignee: Kyoto University
Current assignee: Kyoto University
Priority date: 2006-08-08
Filing date: 2007-01-31
Publication date: 2012-09-05
Anticipated expiration: 2027-01-31
Also published as: US8255447B2; WO2008018188A1; JPWO2008018188A1; US20100185716A1

Description

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

固有値分解は、例えば、分子軌道法による化学計算、統計計算、情報検索等の非常に広い分野で利用されている。また、固有値分解の方法も知られている（例えば、非特許文献１参照）。
Ｊ．Ｊ．Ｍ．Ｃｕｐｐｅｎ，「Ａｄｉｖｉｄｅａｎｄｃｏｎｑｕｅｒｍｅｔｈｏｄｆｏｒｔｈｅｓｙｍｍｅｔｒｉｃｔｒｉｄｉａｇｏｎａｌｅｉｇｅｎｐｒｏｂｌｅｍ」、ＮｕｍｅｒｉｓｃｈｅＭａｔｈｅｍａｔｉｋ，Ｖｏｌ．３６、ｐ．１７７−１９５、１９８１年

近年、これらの応用分野におけるデータ量の増大などに伴って、高速・高精度な固有値分解が求められている。また、固有値分解を並列処理することができる固有値分解装置等の開発が望まれていた。

本発明は、上記状況の下になされたものであり、並列性に優れた高速、高精度な固有値分解を実行可能な固有値分解装置等を提供することを目的とする。

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

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

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

また、本発明による固有値分解装置では、前記固有ベクトル算出部は、前記固有値記憶部から前記対称３重対角行列Ｔの固有値を読み出し、当該固有値のいずれかが負の値である場合に、前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、当該対称３重対角行列Ｔの固有ベクトルを変化させることなく、当該対称３重対角行列Ｔのすべての固有値が正の値となるように、当該対称３重対角行列Ｔを正定値化し、正定値化した後の固有値を前記固有値記憶部に蓄積し、正定値化した後の対称３重対角行列Ｔを前記対角行列記憶部に蓄積する正定値化部と、前記対角行列記憶部からすべての固有値が正の値である対称３重対角行列Ｔを読み出し、前記固有値記憶部から前記対称３重対角行列Ｔの正の値の固有値を読み出し、前記対称３重対角行列Ｔの各要素に関してｑｄ型変換を行うことによって、前記対称３重対角行列Ｔを上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）にコレスキー分解するコレスキー分解部と、前記上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）の各要素と、前記対称３重対角行列Ｔの正の値の固有値とを用いて固有ベクトルを算出して前記固有ベクトル記憶部に蓄積するベクトル算出部と、を備えてもよい。

このような構成により、負の値の固有値が含まれる場合には、対称３重対角行列Ｔを正定値化してからｑｄ型ツイスト分解法を実行することができる。また、ＬＶ型ツイスト分解法よりも高速に計算することができうる。

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

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

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

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

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

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

また、本発明による固有値分解装置では、対称行列Ａが記憶される行列記憶部と、前記対称行列Ａを前記行列記憶部から読み出し、前記対称行列Ａを３重対角化した前記対称３重対角行列Ｔを算出して前記対角行列記憶部に蓄積する対角化部と、をさらに備えてもよい。

このような構成により、任意の対称行列Ａについて固有値を算出することができる。また、対称３重対角行列Ｔの固有ベクトルを用いることにより、対称行列Ａの固有ベクトルを算出することもできうる。

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

また、本発明による固有値分解装置では、前記固有値算出部は、対称３重対角行列Ｔの全ての固有値を算出してもよい。

また、本発明による固有値分解装置では、前記固有ベクトル算出部は、対称３重対角行列Ｔの全ての固有ベクトルを算出してもよい。

また、本発明による固有値分解装置では、前記対称行列Ａは、顔画像データを示すベクトルの平均から算出された共分散行列であってもよい。
このような構成により、固有値分解装置を顔画像データの認識の処理に用いることができうる。

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

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

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

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

対角化部１２は、対称行列Ａを行列記憶部１１から読み出し、その読み出した対称行列Ａを３重対角化した対称３重対角行列Ｔを算出する。そして、対角化部１２は、その算出した対称３重対角行列Ｔを対角行列記憶部１３に蓄積する。対角化部１２は、例えば、ハウスホルダー（Ｈｏｕｓｅｈｏｌｄｅｒ）変換を必要なだけ繰り返し行う方法や、ランチョス（Ｌａｎｃｚｏｓ）法を用いた方法、その他の３重対角化法を用いて、対称行列Ａを３重対角化する。

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

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

固有値分解部１５は、あらかじめ決められた大きさ以下の各対称３重対角行列を対角行列記憶部１３から読み出し、その各対称３重対角行列に対して固有値分解を行い、各対称３重対角行列の固有値と、その各対称３重対角行列の固有ベクトルを算出する。固有値分解部１５は、例えば、２分法と逆反復法を組み合わせた方法、ＭＲ＾３法、ＱＲ法等を用いて固有値分解を行ってもよい。ここで、ＭＲ＾３法による固有値分解を行う場合には、ＬＡＰＡＣＫ３．０（ｈｔｔｐ：／／ｗｗｗ．ｎｅｔｌｉｂ．ｏｒｇ／ｌａｐａｃｋ／）で提供されているＤＳＴＥＧＲを用いてもよい。また、ＱＲ法による固有値分解を行う場合には、ＬＡＰＡＣＫ３．０で提供されているＤＳＴＥＱＲを用いてもよい。これらの固有値分解の方法については、すでに公知であり、その詳細な説明を省略する。固有値分解部１５は、算出した固有値を固有値分解記憶部１６に蓄積する。また、固有値分解部１５は、算出した固有ベクトルからなる直交行列の一部の要素である行列要素も固有値分解記憶部１６に蓄積する。固有ベクトルからなる直交行列とは、固有ベクトルを各列とする直交行列のことである。行列要素の詳細については後述する。

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

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

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

固有ベクトル算出部１９は、対角行列記憶部１３から対称３重対角行列Ｔを読み出し、固有値記憶部１８から対称３重対角行列Ｔの固有値を読み出す。そして、固有ベクトル算出部１９は、対称３重対角行列Ｔとその固有値とから、ツイスト分解法を用いて対称３重対角行列Ｔの固有ベクトルを算出して固有ベクトル記憶部２０に蓄積する。固有ベクトル算出部１９は、ＬＶ型ツイスト分解法により固有ベクトルを算出してもよく、あるいは、ｑｄ型ツイスト分解法により固有ベクトルを算出してもよい。なお、具体的なツイスト分解の手法については後述する。

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

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

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

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

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

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

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

次に、図２のフローチャートのステップＳ１０１，Ｓ１０２の処理、すなわち、対称３重対角行列Ｔの固有値を算出するまでの処理について、より詳細に説明する。まず、図２のフローチャートのステップＳ１０２の処理について、図３のフローチャートを用いて説明する。

（ステップＳ２０１）行列分割部１４は、対角行列記憶部１３から対称３重対角行列Ｔを読み出し、その対称３重対角行列Ｔを２個の対称３重対角行列に分割して対角行列記憶部１３に蓄積する。行列分割部１４は、その対称３重対角行列を２個の対称３重対角行列に分割する処理を、分割後の各対称３重対角行列があらかじめ決められた大きさ以下となるまで繰り返す。この処理の詳細については後述する。

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

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

［対称行列Ａの対角化］
対称行列Ａは、例えば、ハウスホルダー変換等を用いて、以下に示されるように３重対角化することができることが知られている。ここで、対称３重対角行列Ｔは、行の数と列の数とが一致する正方行列である。

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

［対称３重対角行列Ｔの分割］
対称３重対角行列Ｔを分割する処理について説明する。次式で示されるように、ｎ×ｎの対称３重対角行列Ｔが与えられた場合に、それらを２個の対称３重対角行列Ｔ_１，Ｔ_２と、その他の要素（下記のθとβ）とに分けることができる。

ここで、ｎ_１は、１<ｎ_１<ｎとなる整数である。例えば、８×８の対称３重対角行列Ｔを半分の大きさに分割する具体的な処理は、次のようになる。

ただし、β＝ｔ_８であり、θは計算誤差が少なくなるように任意に設定することができる。また、上記各行列において、明記されている以外の行列の要素は０である。なお、上記各行列の明記されている各要素は０以外であってもよく、０であってもよい。ここで、並列処理を適切に実行するためには、ｎ_１を、ｎ／２を超えない最大の整数にとるか、あるいは、ｎ／２を下まわらない最小の整数にとることが好適である（このｎ_１の値を用いて行列を２個の行列に分割する場合を、「行列を略半分の２個の行列に分割する」と呼ぶことにする）。

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

したがって、行列分割部１４は、上述のようにして、対角行列記憶部１３が記憶している対称３重対角行列Ｔを読み出し、その対称３重対角行列Ｔを２個の対称３重対角行列Ｔ_１，Ｔ_２と、２個の値β、θに分割することができ、その分割の処理を繰り返すことができる。行列分割部１４は、分割後の２個の対称３重対角行列Ｔ_１，Ｔ_２と、２個の値β、θとを、対角行列記憶部１３に蓄積する。なお、２個の値β、θは、どの分割で発生した値であるのかがわかるように対角行列記憶部１３で記憶されることが好適である。図４は、図３のフローチャートのステップＳ２０１における行列分割部１４による行列を分割する処理を示すフローチャートである。

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

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

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

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

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

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

［分割された行列の固有値分解］
行列Ｔ_ｉが対称行列であるとすると、行列Ｔ_ｉの固有値分解は、次のようになる。

ここで、Ｄ_ｉは、行列Ｔ_ｉと同じ大きさの対角行列である。Ｄ_ｉの各対角成分はＴ_ｉの固有値である。また、Ｑ_ｉは、行列Ｔ_ｉの固有ベクトルからなる直交行列である。

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

なお、後述するように、固有値を算出するときに用いるために必要なのは、直交行列Ｑ_ｉの全体ではなく、その直交行列Ｑ_ｉの第１行と最終行である。したがって、この固有値分解において、直交行列Ｑ_ｉの全体を求めるのではなく、直交行列Ｑ_ｉの第１行と最終行を求めるようにしてもよい。そのためには、例えば、ＬＡＰＡＣＫ３．０で提供されているＱＲ法の固有値分解ルーチンであるＤＳＴＥＱＲにおいて、直交行列Ｑ_ｉの算出における初期値として、単位行列Ｉを用いる代わりに、次式で示される行列Ｙを用いてもよい。

このようにすることで、ＤＳＴＥＱＲで直交行列Ｑ_ｉの第１行と最終行を求めるための計算領域を削減することができる。直交行列Ｑ_ｉそのものを算出する場合には、Ｏ（ｎ^２）の作業領域が必要であるが、上記の行列Ｙを用いて直交行列Ｑ_ｉの第１行と最終行を算出する場合には、Ｏ（ｎ）の作業領域でよいことになる。その結果、計算量も少なくなり、固有値分解の処理速度も向上する。このように、固有値分解部１５は、分割後の対称３重対角行列に対して、直交行列Ｑ_ｉそのものを算出する固有値分解を行ってもよく、あるいは、直交行列の一部の成分（ここでは、第１行と最終行）を算出する処理を行ってもよい。この後者の処理も固有値分解と呼ぶことにする。このように、固有値分解部１５は、あらかじめ決められた大きさ以下の各対称３重対角行列に対して固有値分解を行い、各対称３重対角行列の固有値と、各対称３重対角行列の固有ベクトルからなる直交行列の行列要素（ここでは、第１行と最終行）とを少なくとも算出し、その固有値と、その行列要素とを固有値分解記憶部１６に蓄積してもよい。

［固有値の算出］
まず、図６で示されるように、分割された２個の対称３重対角行列Ｔ_１，Ｔ_２から、分割元の対称３重対角行列Ｔ_０の固有値等を算出する処理について説明する。対称３重対角行列Ｔ_１，Ｔ_２は、次のように固有値分解されていたとする。ここで、対称３重対角行列Ｔ_１，Ｔ_２は、両者共に正方行列であるとする。

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

上記の式において、β、θは、対称３重対角行列の分割の処理において説明した値である。ここで、行列（Ｄ_１２＋ρｚｚ^Ｔ）の固有値分解が得られれば、行列Ｔ_０の固有値分解が得られたことになる。なぜなら、行列（Ｄ_１２＋ρｚｚ^Ｔ）の固有値分解Ｄ_１２＋ρｚｚ^Ｔ＝ＱＤＱ^Ｔが求まると、次式のようになり、行列Ｔ_０の固有値分解が完了するからである。

次に、行列（Ｄ_１２＋ρｚｚ^Ｔ）の固有値を算出する方法について説明する。行列（Ｄ_１２＋ρｚｚ^Ｔ）の固有値をλ_ｉ、固有ベクトルをｑ_ｉとすると、

となる。

（１）ｚの成分に０が含まれず、Ｄ_１２の対角成分に同じ値が含まれない場合
まず、ｚの成分に０が含まれず、Ｄ_１２の対角成分に同じ値が含まれない場合について説明する。その場合には、上記の式を、次に示すように順次、変形することができる。

ここで、

とすると、次の固有方程式を得ることができる。固有値λ_ｉ（ｉ＝１，２，…，ｎ）は、その固有方程式ｎ個の解である。下記のｆ（λ）は、単調増加であり、かつｆ（λ）＝０の解の範囲が既知であるため、ニュートン法を用いて計算することができる。

なお、Ｄ＝ｄｉａｇ（λ_１，λ_２，…，λ_ｎ）となる。行列Ｑは、行列（Ｄ_１２＋ρｚｚ^Ｔ）の固有ベクトルからなる直交行列であるため、Ｑ＝（ｑ_１，ｑ_２，…，ｑ_ｎ）となる。したがって、固有ベクトルｑ_ｉを求めると、行列Ｑを求めることができる。なお、固有ベクトルｑ_ｉは、次式で与えられる。

（２）ｚの成分に０が含まれる場合
ｚの成分に０が含まれる場合には、行列の次数を下げることができ、その次数を下げた行列を用いて固有値、及び固有ベクトルを求めることができる。この行列の次数を下げる操作をデフレーションと呼ぶ。

ｚ_ｊ＝０の場合には、次式で示すように、（ｄ_ｊ，ｅ_ｊ）が一組の固有値と固有ベクトルとなる。ただし、ｅ_ｊは、ｊ番目の成分が１であり、その他の成分が０である単位ベクトルである。

また、行列（Ｄ_１２＋ρｚｚ^Ｔ）の固有方程式は、次式のようになる。

したがって、ｚ_ｊ＝０の場合の固有値は、ｄ_ｊと、行列（Ｄ_１２＋ρｚｚ^Ｔ）の第ｊ行、第ｊ列を除いた小行列について、上記式１の固有方程式を解いた値である。すなわち、固有値は、ｄ_ｊと、次の固有方程式を解いた値となる。なお、次の固有方程式を解いた固有値を、固有値λ_ｉ（ｉ＝１，２，…，ｊ−１，ｊ＋１，…，ｎ）とする。

固有ベクトルは、行列（Ｄ_１２＋ρｚｚ^Ｔ）の第ｊ行、第ｊ列を除いた小行列について上記式２で算出したｎ−１次元の固有ベクトルの第（ｊ−１）行と、第ｊ行との間に０を挿入した次式のｎ次元のベクトルと、上述のｅ_ｊとなる。

このように、ｚの成分に０が含まれる場合には、行列（Ｄ_１２＋ρｚｚ^Ｔ）よりも次元の小さい行列について固有値分解を行うことになるため、計算を高速化することが可能であり、結果として、計算時間が短縮されうる。

（３）Ｄ_１２の対角成分に同じ値が含まれる場合
ｄ_ｊ＋１＝ｄ_ｊ＋２＝…＝ｄ_ｋの場合にも、上記説明と同様のデフレーションを行うことによって、固有値と固有ベクトルとを算出することができる。

に対して、次の２個の直交行列を考える。

すると、ｄ_ｊ＋１＝ｄ_ｊ＋２＝…＝ｄ_ｋより、

となる。したがって、次式が得られる。

ここで、

の固有値と固有ベクトルとは、上記のｚの成分に０が含まれるデフレーションの場合と同様にして求めることができる。また、

は、符号の任意性を除いて一意に定まる。ここで、その算出方法について説明する。まず、次式を満たす直交行列について考える。

この直交行列を次式のようにとると、上記式を満たすことになる。

ここで、上記の式の＊の値は、次式から求まる。

したがって、

となる。その結果、

も、符号の任意性を除いて求めることができる。このようにして、ｄ_ｊ＋１＝ｄ_ｊ＋２＝…＝ｄ_ｋの場合にも行列（Ｄ_１２＋ρｚｚ^Ｔ）の固有値と固有ベクトルとを求めることができる。すなわち、固有値は、ｄ_ｊ＋２，ｄ_ｊ＋３，…，ｄ_ｋと、次の固有方程式を解いた値となる。なお、次の固有方程式を解いた固有値を、固有値λ_ｉ（ｉ＝１，２，…，ｊ＋１，ｋ＋１，…，ｎ）とする。

固有ベクトルは、次のようになる。

なお、上記（２）ｚの成分に０が含まれる場合、（３）Ｄ_１２の対角成分に同じ値が含まれる場合において、ｚの成分が０であるかどうかの判断や、Ｄ_１２の対角成分に同じ値が含まれているかどうかの判断は、微少な誤差を考慮して行ってもよい。例えば、次式が満たされる場合には、ｚの成分が０であると判断してもよく、また、Ｄ_１２の対角成分が同じ値であると判断してもよい。

なお、εは、マシン・イプシロンに定数（例えば、２，４，８などの実数であり、計算誤差を少なくするための値を選択することが好適である）を掛けた値を用いてもよく、他の適切な値であってもよい。

このようにして、（１）ｚの成分に０が含まれず、Ｄ_１２の対角成分に同じ値が含まれない場合、（２）ｚの成分に０が含まれる場合、（３）Ｄ_１２の対角成分に同じ値が含まれる場合のそれぞれについて、固有値分解Ｄ_１２＋ρｚｚ^Ｔ＝ＱＤＱ^Ｔが求まることになる。したがって、行列Ｑと行列Ｑ_１２から行列Ｑ_０を求めることができ、行列Ｄも求まっているため、分割元の対称３重対角行列Ｔ_０の固有値分解を行うことができる。このように、対象となる行列を２個の副行列に分割し、その分割した後の各副行列について固有値分解を行い、その後に各副行列の固有値分解を用いて分割元の行列の固有値分解を行う方法は、分割統治法（ＤｉｖｉｄｅａｎｄＣｏｎｑｕｅｒ：Ｄ&Ｃ）と呼ばれている。

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

一方、分割元の行列Ｔ_０の固有値のみを算出する場合には、上記処理を簡略化することができる。次に、その方法について説明する。上記の結果より、次のようになる。

ここで、上記のように、ｆは行列Ｔ_０を固有値分解した結果の直交行列の最初の行であり、ｌは行列Ｔ_０を固有値分解した結果の直交行列の最後の行である。このｆとｌ、すなわち、直交行列の最初の行と、最後の行とが行列要素となる。なお、直交行列Ｑ_１２について、後述する列の入れ替えを行った場合には、ｆ_１２、ｌ_１２はそれぞれ、その列を入れ替えた後の直交行列Ｑ_１２の第１行、最終行となる。

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

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

なお、上記の（３）Ｄ_１２の対角成分に同じ値が含まれる場合には、式３，式４が次のようになる。

ここで、上記式を計算する場合に、

を先に計算して、その後に

を掛けた方が、計算量が少なくなるため、その順番で計算するようにしてもよい。

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

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

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

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

図８は、図７のフローチャートにおけるステップＳ４０３の詳細な処理を示すフローチャートである。なお、この図８のフローチャートの処理を実行する前に、対角行列Ｄ_１２の対角成分が昇順、あるいは降順となるように並べ替えておくことが好適である。なお、その場合には、直交行列Ｑ_１２の列、及びｚの成分も、対角行列Ｄ_１２の並び替えと同様に並び替えるものとする。具体的には、対角行列Ｄ_１２の第ｉ成分と、第ｊ成分とを入れ替えた場合には、直交行列Ｑ_１２の第ｉ列と、第ｊ列とを入れ替え、ｚの第ｉ行と、第ｊ行とを入れ替えればよい。

（ステップＳ５０１）固有値算出部１７は、分割元の行列の固有値を、式１等を用いて算出する。そして、固有値算出部１７は、算出した固有値を固有値分解記憶部１６に蓄積する。

（ステップＳ５０２）固有値算出部１７は、ステップＳ５０１で算出した固有値を用いて、式２等からｑ_ｉを算出する。このｑ_ｉを算出することにより、ｑ_ｉを列ベクトルに有する行列Ｑを算出したことになる。固有値算出部１７は、算出した行列Ｑを図示しない記録媒体において一次記憶してもよい。

（ステップＳ５０３）固有値算出部１７は、ステップＳ５０２で算出した行列Ｑを用いて、分割元の行列に関する行列要素ｆ、ｌを算出する。この行列要素は、式３，式４で示されるものである。そして、固有値算出部１７は、算出した分割元の行列要素を固有値分解記憶部１６に蓄積する。このようにしてステップＳ４０３の処理は終了となる。

なお、図８のフローチャートにおいて、ｚの成分に０が含まれず、Ｄ_１２の対角成分に同じ値が含まれない場合には、上記（１）で説明したようにして、固有値、固有ベクトルからなる行列Ｑを算出することができる。また、ｚの成分に０が含まれ、Ｄ_１２の対角成分に同じ値が含まれない場合には、上記（２）で説明したようにして、固有値、固有ベクトルからなる行列Ｑを算出することができる。また、Ｄ_１２の対角成分に同じ値が含まれる場合には、ｚの成分に０が含まれるかどうかに関わらず、上記（３）で説明したようにして、固有値、固有ベクトルからなる行列Ｑを算出することができる。上記（３）で説明したようにして固有値等を算出する場合には、前述のように、行列Ｑを算出するのではなく、行列Ｑを構成する行列を直接ｆ_１２等に作用させることによって、ｆ、ｌを算出するようにしてもよい。

また、図８のフローチャートでは、行列Ｑから、行列要素ｆ、ｌを一度に算出する場合について示したが、行列要素ｆ、ｌをその成分ごとに計算してもよい。上記の式３，式４から明らかなように、ｆ_１２、ｌ_１２に固有ベクトルｑ_ｉを掛けたものが、行列要素ｆ、ｌのｉ番目の成分となる。このような方法で行列要素ｆ、ｌを計算する場合には、固有ベクトル固有ベクトルｑ_ｉのみをｉの値を順次、変えながら記憶することによってその計算を行うことができる。したがって、行列要素ｆ、ｌの計算において行列Ｑ全体を記憶する必要がなくなり、その計算で用いるメモリ上の作業領域を削減することが可能である。

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

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

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

なお、図９の行列Ｔ_１、行列Ｔ_２の固有値、行列要素等に基づいて行列Ｔの固有値を算出する場合には、行列要素ｆ_１、ｌ_２、あるいは行列要素ｆ_２、ｌ_１を用いるだけでよい。したがって、行列Ｔ_１、行列Ｔ_２の行列要素を算出する場合に、行列要素ｆ_１、ｌ_２、あるいは行列要素ｆ_２、ｌ_１を計算しなくてもよい。このようにすることで、計算量を削減することができ、処理スピードを向上させることが可能である。

次に、固有ベクトル算出部１９について説明する。本実施の形態では、固有ベクトル算出部１９が（１）ｑｄ型ツイスト分解法により固有ベクトルを算出する場合と、（２）ＬＶ型ツイスト分解法により固有ベクトルを算出する場合とについて説明する。

（１）ｑｄ型ツイスト分解法により固有ベクトルを算出する場合
図１０は、ｑｄ型ツイスト分解法により固有ベクトルを算出する場合の固有ベクトル算出部１９の構成を示す図である。図１０において、固有ベクトル算出部１９は、正定値化部２１と、コレスキー分解部２２と、ベクトル算出部２３とを備える。

正定値化部２１は、固有値記憶部１８から対称３重対角行列Ｔの固有値を読み出し、その固有値のいずれかが負の値であるかどうか判断する。そして、その固有値のいずれかが負の値である場合に、正定値化部２１は、対角行列記憶部１３から対称３重対角行列Ｔを読み出し、その対称３重対角行列Ｔの固有ベクトルを変化させることなく、対称３重対角行列Ｔのすべての固有値が正の値となるように、その対称３重対角行列Ｔを正定値化し、正定値化した後の固有値を固有値記憶部１８に蓄積し、正定値化した後の対称３重対角行列Ｔを対角行列記憶部１３に蓄積する。この正定値化の処理を行うのは、ｑｄ型ツイスト分解法により固有ベクトルを算出する場合には、対称３重対角行列Ｔの固有値がすべて正の値でなければならないからである。なお、対称３重対角行列Ｔを正定値化した行列も、対称３重対角行列Ｔと呼ぶことにする。また、この正定値化の処理が行われた場合には、コレスキー分解や固有ベクトルを算出する処理は、正定値化後の対称３重対角行列Ｔ、固有値に対して行われることになる。この処理の詳細については後述する。ここで、正定値化の行われた場合であっても、固有値記憶部１８には、正定値化の前の固有値が記憶されているものとする。その正定値化の前の固有値が、固有値分解装置１が最終的に求める固有値だからである。また、正定値化の行われた場合であって、ｑｄ型ツイスト分解法による固有ベクトルの算出のみを行う場合には、対角行列記憶部１３で記憶されている対称３重対角行列Ｔを、正定値化の後の対称３重対角行列Ｔで上書きしてもよい。一方、ＬＶ型ツイスト分解法による固有ベクトルの算出も行う場合には、対角行列記憶部１３には、正定値化の前の対称３重対角行列Ｔと、正定値化の後の対称３重対角行列Ｔとが記憶されることになる。

コレスキー分解部２２は、対角行列記憶部１３から、すべての固有値が正の値である対称３重対角行列Ｔを読み出し、固有値記憶部１８から、対称３重対角行列Ｔの正の値の固有値を読み出す。そして、コレスキー分解部２２は、対称３重対角行列Ｔの各要素に関してｑｄ型変換を行うことによって、対称３重対角行列Ｔを上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）にコレスキー分解する。この処理の詳細については後述する。なお、「対称３重対角行列Ｔを上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）にコレスキー分解する」ことには、対称３重対角行列Ｔの固有値に単位ベクトルを掛けたものを対称３重対角行列Ｔから引いた行列を、上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）にコレスキー分解することが含まれてもよい。

ベクトル算出部２３は、コレスキー分解部２２が算出した上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）の各要素と、対称３重対角行列Ｔの正の値の固有値とを用いて固有ベクトルを算出して固有ベクトル記憶部２０に蓄積する。この固有ベクトルを各列に有する行列が、固有値分解の直交行列となる。

次に、図２のフローチャートのステップＳ１０３の処理、すなわち、対称３重対角行列Ｔの固有ベクトルを算出する処理について、より詳細に説明する。まず、図２のフローチャートのステップＳ１０３の処理について、図１１のフローチャートを用いて説明する。

（ステップＳ６０１）正定値化部２１は、固有値記憶部１８から対称３重対角行列Ｔの固有値を読み出し、その固有値のいずれかが負の値であるかどうか判断する。そして、負の値がある場合には、ステップＳ６０２に進み、そうでない場合には、ステップＳ６０３に進む。

（ステップＳ６０２）正定値化部２１は、対角行列記憶部１３から対称３重対角行列Ｔを読み出し、その対称３重対角行列Ｔの固有ベクトルを変化させることなく、対称３重対角行列Ｔのすべての固有値が正の値となるように、その対称３重対角行列Ｔを正定値化し、正定値化した後の対称３重対角行列Ｔを対角行列記憶部１３に蓄積し、正定値化した後の固有値を固有値記憶部１８に蓄積する。

（ステップＳ６０３）コレスキー分解部２２は、対角行列記憶部１３から、すべての固有値が正の値である対称３重対角行列Ｔを読み出し、固有値記憶部１８から、対称３重対角行列Ｔの正の値の固有値を読み出す。なお、正定値化が行われた場合には、正定値化の後の対称３重対角行列Ｔ、及び固有値が読み出されることになる。そして、コレスキー分解部２２は、対称３重対角行列Ｔの各要素に関してｑｄ型変換を行うことによって、対称３重対角行列Ｔを上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）にコレスキー分解する。

（ステップＳ６０４）ベクトル算出部２３は、コレスキー分解部２２が算出した上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）の各要素と、対称３重対角行列Ｔの正の値の固有値とを用いて固有ベクトルを算出する。

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

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

［対称３重対角行列Ｔの正定値化］
行列Ｔを次のように固有値分解できたとする。

その行列Ｔを次のようにａ倍してｓＩだけシフトしたとする。ここで、Ｉは行列Ｔと同じ大きさの単位行列であり、ａ，ｓは実数である。

すると、シフト後の固有値分解は、次のようになる。

このように、行列Ｔをａ倍してｓＩだけシフトした場合には、各固有値がａ倍されてｓだけシフトされるが、固有ベクトルは変化しない。したがって、この実数倍とシフトとによって、すべての固有値を正の値にしてから、ツイスト分解法を用いた固有ベクトルの算出を行えばよいことになる。

なお、算出された固有値に負の値が含まれる場合には、例えば、（Ａ）最小の固有値によりシフトを行ってもよく、（Ｂ）最大の固有値によりシフトを行ってもよい。ここで、固有値には、負の値と正の値とが含まれているものとする。固有値がすべて正の値であれば、シフトを行う必要はなく、固有値がすべて負の値であれば、単に行列Ｔを−Ｔとすれば、正定値化することができる。

（Ａ）最小の固有値によるシフト
最小の固有値をλ_ｍｉｎ<０とすると、ａ＝１、ｓ＝−λ_ｍｉｎ＋λ_ｍａｘ×εとする。ここで、εとしては、マシン・イプシロンに定数を掛けた値を用いてもよく、他の適切な値であってもよい。このように、行列Ｔをシフトすることによって、行列Ｔのすべての固有値を正の値にすることができる。

（Ｂ）最大の固有値によるシフト
最大の固有値をλ_ｍａｘ>０とすると、ａ＝−１、ｓ＝λ_ｍａｘ−λ_ｍｉｎ×εとする。この場合には、次のように行列Ｔが変換されることになる。

上式の右辺の括弧内では、シフトによってすべての固有値が負の値になるが、全体にマイナスを掛けることによって、正定値化することができている。ここで、εは、（Ａ）の場合と同様にマシン・イプシロンに定数（例えば、２，３，４などの実数であり、計算誤差を少なくするための値を選択することが好適である）を掛けた値を用いてもよく、他の適切な値であってもよい。

固有値の最小値の絶対値が、固有値の最大値の絶対値よりも小さい場合には、上記（Ａ）のシフトを行い、固有値の最小値の絶対値が、固有値の最大値の絶対値よりも大きい場合には、上記（Ｂ）のシフトを行うようにしてもよい。このようにすることで、シフトの量を少なくすることができ、丸め誤差を削減することが可能だからである。
なお、正定値化の後に固有ベクトルを算出する処理の詳細については、次のＬＶ型ツイスト分解法の場合と一緒に説明する。

（２）ＬＶ型ツイスト分解法により固有ベクトルを算出する場合
図１３は、ＬＶ型ツイスト分解法により固有ベクトルを算出する場合の固有ベクトル算出部１９の構成を示す図である。図１３において、固有ベクトル算出部１９は、コレスキー分解部３１と、ベクトル算出部３２とを備える。

コレスキー分解部３１は、対角行列記憶部１３から、対称３重対角行列Ｔを読み出し、固有値記憶部１８から、対称３重対角行列Ｔの固有値を読み出す。そして、コレスキー分解部３１は、対称３重対角行列Ｔの各要素に関してミウラ変換、ｄＬＶｖ型変換、逆ミウラ変換を行うことによって、対称３重対角行列Ｔを上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）にコレスキー分解する。この処理の詳細については後述する。なお、「対称３重対角行列Ｔを上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）にコレスキー分解する」ことには、対称３重対角行列Ｔの固有値に単位ベクトルを掛けたものを対称３重対角行列Ｔから引いた行列を、上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）にコレスキー分解することが含まれてもよい。

ベクトル算出部３２は、コレスキー分解部３１が算出した上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）の各要素と、対称３重対角行列Ｔの固有値とを用いて固有ベクトルを算出して固有ベクトル記憶部２０に蓄積する。この固有ベクトルを各列に有する行列が、固有値分解の直交行列となる。

次に、図２のフローチャートのステップＳ１０３の処理、すなわち、対称３重対角行列Ｔの固有ベクトルを算出する処理について、より詳細に説明する。まず、図２のフローチャートのステップＳ１０３の処理について、図１４のフローチャートを用いて説明する。

（ステップＳ８０１）コレスキー分解部３１は、対角行列記憶部１３から対称３重対角行列Ｔを読み出し、固有値記憶部１８から対称３重対角行列Ｔの固有値を読み出す。そして、コレスキー分解部３１は、対称３重対角行列Ｔの各要素に関してミウラ変換、ｄＬＶｖ型変換、逆ミウラ変換を行うことによって、対称３重対角行列Ｔを上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）にコレスキー分解する。

（ステップＳ８０２）ベクトル算出部３２は、上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）の各要素と、対称３重対角行列Ｔの固有値とを用いて固有ベクトルを算出する。

（ステップＳ８０３）ベクトル算出部３２は、算出した固有ベクトルを正規化する。すなわち、ベクトル算出部３２は、算出した固有ベクトルのノルムを算出し、固有ベクトルを算出したノルムで割ったものを最終的な固有ベクトルとして固有ベクトル記憶部２０に蓄積する。このようにして、対称３重対角行列Ｔを固有値分解する処理は終了となる。

なお、この後に、行列記憶部１１が記憶している対称行列Ａの固有ベクトルを算出する処理が行われてもよいが、ここでは省略している。また、固有ベクトル算出部１９が算出した固有ベクトルの精度を上げるために、図１２で示されるように、固有ベクトルに関する処理を実行してもよく、あるいは、実行しなくてもよいことは、ｑｄ型ツイスト分解法の場合と同様である。

［固有値からの固有ベクトルの算出］
固有値から固有ベクトルを算出する処理については、ｑｄ型ツイスト分解法を用いる方法と、ＬＶ型ツイスト分解法を用いる方法とについて説明する。図１１のフローチャートで示されるように、ｑｄ型ツイスト分解法を用いて固有ベクトルを算出する場合に、負の固有値があれば、対称３重対角行列Ｔを正定値化する処理をあらかじめ実行しておく（ステップＳ６０１，Ｓ６０２）。ｑｄ型ツイスト分解法を用いる場合には、この後の処理において、その正定値化の後の対称３重対角行列Ｔ、及び固有値を用いるものとする。

対角行列記憶部１３では、次のようなｍ×ｍの対称３重対角行列Ｔが記憶されているものとする。

また、固有値記憶部１８では、固有値λ_ｉ（１≦ｉ≦ｍ）が記憶されているものとする。このような場合に、対称３重対角行列Ｔの固有ベクトルを求めることは、次の方程式の固有ベクトルｘ_ｊを求めることになる。

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

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

ここで、行列Ｔの各要素と、行列Ｂ^（０）の各要素との対応は、次のようになる。

なお、行列Ｔの各要素と、行列Ｂ^（０）の各要素との対応を次のように書くこともできる。下記の式から明らかなように、行列Ｂ^（０）の各要素を、行列Ｔの各要素を用いて順次、算出することができる。

上記のコレスキー分解した結果を、次式のように書くことができる。

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

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

であるので、（Ｔ−λ_ｊＩ）ｘ_ｊ＝γ_ｋｅ_ｋは、

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

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

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

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

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

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

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

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

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

なお、このように対角行列Ｔの各要素を、行列Ｂ^（０）の各要素に変換して、その行列Ｂ^（０）の各要素についてｑｄ型変換を行う場合であっても、間接的に対角行列Ｔの各要素に関してｑｄ型変換を行うことになるため、「対角行列Ｔの各要素に関してｑｄ型変換を行う」と言うことにする。
また、ここではスタンダードなｑｄ型変換について説明したが、ｑｄ型変換はｄｉｆｆｅｒｅｎｔｉａｌｑｄ型変換であってもよい。

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

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

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

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

に分けられる。ただし、

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

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

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

なお、このように対角行列Ｔの各要素を、行列Ｂ^（０）の各要素に変換して、その行列Ｂ^（０）の各要素についてＬＶ型変換を行う場合であっても、間接的に対角行列Ｔの各要素に関してＬＶ型変換を行うことになるため、「対角行列Ｔの各要素に関してＬＶ型変換を行う」と言うことにする。

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

図１６は、コレスキー分解の全体の処理の一例を示す図である。
（ステップＳ９０１）コレスキー分解部３１は、ミウラ変換を行う。この処理の詳細については後述する。

（ステップＳ９０２）コレスキー分解部３１は、１／δ^（±１）を１／δ^（０）−λ_ｊとする。
（ステップＳ９０３）コレスキー分解部３１は、後述するＰｒｏｃｅｄｕｒｅ１の処理を実行する。

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

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

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

図１９は、図１６のフローチャートのステップＳ９０１の処理の詳細を示すフローチャートである。
（ステップＳ１２０１）コレスキー分解部３１は、１／δ^（０）の値を決定する。この値は、前述のように任意に決定することができるため、行列を正定値化する値に決定する。例えば、コレスキー分解部３１は、１／δ^（０）の値を次式のように決定する。

なお、なお、ε^（０）は、マシン・イプシロンに定数（例えば、２，３，４などの実数であり、計算誤差を少なくするための値を選択することが好適である）を掛けた値を用いてもよく、他の適切な値であってもよい。その後、例えば、ステップＳ１２０３等で桁落ちの発生する可能性があると判断された場合に、コレスキー分解部３１は、１／δ^（０）の値を決めるε^（０）の値を、例えば、マシン・イプシロンに掛ける定数を１ずつ大きくしながら設定していってもよい。なお、１／δ^（０）の値を決定する方法は、これに限定されないことは言うまでもない。このように、ＬＶ型ツイスト分解法では、この１／δ^（０）の値を自由に設定することができるため、ｑｄ型ツイスト分解法の場合のように、行列Ｔの正定値化を行わなくてもよいことになる。なお、ｑｄ型ツイスト分解法の場合と同様に、行列Ｔの正定値化の処理を行ってから、ＬＶ型ツイスト分解法を実行してもよいことは言うまでもない。
（ステップＳ１２０２）コレスキー分解部３１は、β１をｑ_１ ^（０）−１／δ^（０）に設定する。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

なお、前述の行列Ｔの各要素と、行列Ｂ^（０）の各要素との対応を示す式を用いると、ｑｄ型変換を次のように書き換えることもできる。

この変換は、さらにｓｔｑｄｓ変換

及びｒｐｑｄｓ変換

に分けられる。このように、｛ｔ_２ｋ−１，ｔ_２ｋ｝から｛ｑ_ｋ ^（±１），ｅ_ｋ ^（±１）｝を求める場合には、｛ｔ_２ｋ−１，ｔ_２ｋ｝から｛ｑ_ｋ ^（０），ｅ_ｋ ^（０）｝を求め、さらにｑ_ｋ ^（０），ｅ_ｋ ^（０）｝から｛ｑ_ｋ ^（±１），ｅ_ｋ ^（±１）｝を求める場合に比べて、計算量が少なくなることとなり、計算速度を向上させることも可能である。

なお、この場合には、対角行列Ｔの各要素について直接ｑｄ型変換が行われることになる。したがって、この変換が、「対角行列Ｔの各要素に関してｑｄ型変換を行う」場合に含まれることは言うまでもない。

また、ｑｄ型変換と同様に、前述の行列Ｔの各要素と、行列Ｂ^（０）の各要素との対応を示す式を用いると、ＬＶ型変換におけるミウラ変換を次のように書き換えることもできる。

この場合のミウラ変換は、次のように示される。

このように、｛ｔ_２ｋ−１，ｔ_２ｋ｝から｛ｕ_２ｋ−１ ^（０），ｕ_２ｋ ^（０）｝を求める場合には、｛ｔ_２ｋ−１，ｔ_２ｋ｝から｛ｑ_ｋ ^（０），ｅ_ｋ ^（０）｝を求め、さらにｑ_ｋ ^（０），ｅ_ｋ ^（０）｝から｛ｕ_２ｋ−１ ^（０），ｕ_２ｋ ^（０）｝を求める場合に比べて、計算量が少なくなることとなり、計算速度を向上させることも可能である。

なお、この場合には、対角行列Ｔの各要素について直接ＬＶ型変換が行われることになる。したがって、この変換が、「対角行列Ｔの各要素に関してＬＶ型変換を行う」場合に含まれることは言うまでもない。

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

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

したがって、次式のように対称３重対角行列Ｔについて固有値分解が行われたことになる。

また、対称３重対角行列Ｔの固有値分解を行うことができれば、次のように、対称行列Ａの固有値分解も行われる。

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

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

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

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

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

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

次に、ベクトル算出部３２は、算出したｘ_ｊを正規化して、固有ベクトルｘ_ｊを求めて固有ベクトル記憶部２０に蓄積する（ステップＳ８０３）。なお、ベクトル算出部３２は、各固有値について、ｘ_ｊを算出する処理と、正規化する処理とを行うことによって、全ての固有ベクトルを算出することができる。このようにして、固有ベクトルの算出が終了する。
この後、行列記憶部１１が記憶している対称行列Ａの固有ベクトルを算出してもよい。その算出の方法は、前述の通りである。

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

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

また、上記説明では、行列Ｂ^（０）が上２重対角行列である場合について説明したが、行列Ｂ^（０）は下２重対角行列であっても固有値分解を同様に実行することができうる。ただし、行列Ｔの各要素と、行列Ｂ^（０）の各要素との対応を示す式が上述のものと異なることになる。

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

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

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

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

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

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

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

コレスキー分解部３１が複数のコレスキー分解手段３１ａ、３１ｂを備えて、対称３重対角行列Ｔに関してコレスキー分解する処理を並列実行してもよいことは、コレスキー分解部２２と同様である。

ベクトル算出部３２が複数のベクトル算出手段３２ａ、３２ｂを備えて、固有ベクトルを算出する処理を並列実行してもよいことは、ベクトル算出部２３と同様である。

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

まず、式１から、各固有値を求める処理については、並列実行可能なことがわかる。例えば、固有値算出手段１７ａが半分の固有値を算出し、固有値算出手段１７ｂが残りの半分の固有値を算出してもよい。なお、固有値算出手段１７ａ、１７ｂのそれぞれは、行列Ｄ_１２と、ベクトルｚとの各成分の値を有するものとする。

また、式２から、各固有ベクトルｑ_ｉを求める処理についても、並列実行可能なことがわかる。例えば、固有値算出手段１７ａは、算出した固有値に対応する半分の固有ベクトルｑ_ｉを算出し、固有値算出手段１７ｂは、算出した固有値に対応する残りの半分の固有ベクトルｑ_ｉを算出してもよい。

また、式３，式４から、ｆ、ｌを求める処理についても、並列実行可能なことがわかる。例えば、固有値算出手段１７ａは、算出した固有ベクトルｑ_ｉを用いて、ｆ、ｌの一部の要素を算出し、固有値算出手段１７ｂは、算出した固有ベクトルｑ_ｉを用いて、ｆ、ｌの残りの要素を算出してもよい。その後、固有値算出手段１７ａ、あるいは固有値算出手段１７ｂが、固有値算出手段１７ａ、１７ｂの算出したｆ、ｌの各要素を統合することによって、式３，式４で示されるｆ、ｌを算出することができる。

ここでは、ｚの成分に０が含まれず、Ｄ_１２の対角成分に同じ値が含まれない場合における固有値算出手段１７ａ、１７ｂの処理について説明したが、ｚの成分に０が含まれる場合、Ｄ_１２の対角成分に同じ値が含まれる場合であっても、固有方程式を解く処理、固有ベクトルｑ_ｉを算出する処理、固有ベクトルから分割元の行列の行列要素を算出する処理を同様に並列実行することができうる。

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

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

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

次に、数値実験による性能の評価について説明する。ここでは、本実施の形態による固有値分解装置１の方法（ＤＤＣ）と、標準的な分割統治法（ＤＣ）と、ＱＲ法（ＱＲ）と、二分法と逆反復法とを組み合わせた方法（Ｂ&Ｉ）とを比較した。なお、本実施の形態による固有値分解装置１の方法では、コレスキー分解をｑｄ型ツイスト分解法で行い、逆反復（図１２のステップＳ７０１参照）を１回用いた。分割統治法としては、ＬＡＰＡＣＫで提供されているＤＳＴＥＤＣを用いた。ＱＲ法としては、ＬＡＰＡＣＫで提供されているＤＳＴＥＱＲを用いた。二分法と逆反復法とを組み合わせた方法としては、ＬＡＰＡＣＫで提供されているＤＳＴＥＶＲのソースコードを一部変更したものを用いた。実験で用いるＬＡＰＡＣＫはＢＬＡＳ（ＢａｓｉｃＬｉｎｅａｒＡｌｇｅｂｒａＳｕｂｐｒｏｇｒａｍｓ）を呼び出す。また、本実施の形態による固有値分解装置１の固有値分解部１５では、ＱＲ法を用いて固有値分解を行った。また、この数値実験では、ランダムに与えた対角行列を、ランチョス法を用いて対称３重対角化した行列を用いた。数値実験結果のグラフにおいて、横軸は、数値実験で用いた行列の大きさである。

図２８は、計算の実行時間を示す図である。実行時間が短いほど、高速に計算できることになる。図２８の数値実験結果のグラフからわかるように、本実施の形態による固有値分解装置１の方法では、他の方法に比べて高速に固有値分解を行うことができることがわかる。

図２９は、算出された固有ベクトルの直交性を示す図である。値の小さい方が、直交性が高いことになる。図２９の数値実験結果のグラフからわかるように、本実施の形態による固有値分解装置１の方法では、他の方法に比べて、固有ベクトルの直交性は高くない。本実施の形態による固有値分解装置１の方法では、二分法と逆反復法とを組み合わせた方法と同程度の直交性を有することがわかる。これは、固有値をまず算出し、その後に、固有ベクトルを算出していることに起因していると考えられる。

図３０は、全固有値の相対誤差の総和を示す図である。値の小さい方が、誤差が小さく、精度が高いことになる。図３０の数値実験結果のグラフからわかるように、本実施の形態による固有値分解装置１の方法では、他の方法に比べて固有値の精度が高いことがわかる。分割統治法により算出した固有値と同程度の固有値の精度を有する。本実施の形態による固有値分解装置１の方法では、分割統治法によって固有値を算出しているため、当然の結果であると考えられる。

図３１は、前固有ベクトルの誤差の総和を示す図である。値の小さい方が、誤差が小さく、精度が高いことになる。図３１の数値実験結果のグラフからわかるように、本実施の形態による固有値分解装置１の方法では、他の方法に比べて固有ベクトルの精度が高いことがわかる。

このように、本実施の形態による固有値分解装置１の方法では、高速に固有値分解をすることができ、さらに、固有値、固有ベクトルの精度の高いことがわかる。また、前述のように、並列性にも優れているため、並列処理を行うことによってさらに高速化を図ることも可能である。

［主成分分析による顔画像認識］
次に、主成分分析による顔画像認識に固有値分解を応用する場合について説明する。

主成分分析による顔画像認識の処理は、以下の各ステップによって行われる。
（１）射影行列の作成
（２）個人の顔画像に関する情報の登録
（３）入力された顔画像と、登録されている個人の顔画像に関する情報とを比較することによる認識

まず、上記（１）の射影行列の作成について説明する。
種々の人物の顔画像データを用意する。その顔画像データは、例えば、スキャナやデジタルスチルカメラ、デジタルビデオカメラ等の光学的読み取り機器によって読み取られた２次元画像データであってもよい。各顔画像データは、各画素値を成分とするベクトルで表示することができる。そのベクトルを、次のようにする。ここでは、Ｍ個の顔画像データを扱うものとする。

次に、その顔画像データのベクトルの平均値μを算出する。

その算出したμを用いて、共分散行列Ｃを次のように算出する。

その後、上述の固有値分解装置１によって、その行列Ｃに対して、固有値分解を行う。すなわち、行列記憶部１１に行列Ｃを蓄積し、その行列Ｃに対して対角化や、行列の分割等の処理を行うことによって、固有値分解を行う。その固有値分解された結果である固有値λ_ｊは、固有値記憶部１８に蓄積されることになる。また、固有値分解された結果である固有ベクトルｖ_ｊは、固有ベクトル記憶部２０で記憶されている固有ベクトルに対して、対角化で用いた直交行列を作用させることによって、前述のようにして求めることができる。なお、添字ｊは、固有値の降順になるように設定されているものとする。すなわち、固有値λ_１が最大の固有値となる。

その求めた固有ベクトルのうち、対応する固有値が上位からｄ個の固有ベクトルを用いて、射影行列ηを次のように算出する。

次に、上記（２）の個人の顔画像に関する情報の登録について説明する。
各個人について、ｍ個の顔画像データ

を用意する。ここで、添字ｋは、個人を識別するための添字である。そのベクトルに射影行列ηをそれぞれ掛けることにより、次式のようにしてｍ個のｄ次元ベクトルξ_ｉを算出し、その平均値μ^ｋを算出する。そして、その平均値μ^ｋを登録しておく。

次に、上記（３）の入力された顔画像と、登録されている個人の顔画像に関する情報とを比較することによる認識について説明する。
ある人物の顔画像データのベクトルｘが入力されたとする。すると、その顔画像データのベクトルｘに射影行列ηを掛けることにより、ｄ次元ベクトルξを算出する。

次に、そのｄ次元ベクトルξと、各個人のｄ次元ベクトルξ_ｉの平均値μ^ｋとの差のノルムを計算する。

そのノルムがあらかじめ設定されているしきい値より小さい場合に、入力された顔画像が、添字ｋで識別される個人の顔画像であると判断してもよい。あるいは、すべての添字ｋについて上記のノルムを計算し、入力された顔画像が、最小のノルムに対応する添字ｋで識別される個人の顔画像であると判断してもよい。このようにして、固有値分解装置１を用いて、主成分分析による顔画像認識を行うことができうる。

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

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

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

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

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

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

そこで、行列Ｄの特異値分解
Ｄ＝ＵΣＶ^Ｔ
を考える。ここで、Σは特異値を大小順に対角線上に並べたもので、Ｕ及びＶはそれぞれ左直交行列、及び右直交行列である。この特異値分解において、前述の固有値分解を用いることができる。すなわち、固有値分解装置１において、行列記憶部１１が記憶している対称行列Ａを、Ｄ^ＴＤとすることにより、前述の説明のようにして、対称行列Ａの固有値分解がなされることになる。なお、固有値分解装置１において、固有ベクトル記憶部２０に蓄積されるのは、対称行列Ｄ^ＴＤを変換した対称３重対角行列Ｔの固有ベクトルであるので、その固有ベクトルを、前述の説明のようにして、行列Ｄ^ＴＤの固有ベクトルに変換する必要がある。また、その行列Ｄ^ＴＤの各固有値の１／２乗（平方根）が各特異値となる。すなわち、特異値σ_ｊ＝（λ_ｊ）^１／２となる。また、固有ベクトルと、右直交行列Ｖの各列を構成する右特異ベクトルｘ_ｊとは等しいため、固有ベクトルを各列に有する行列が、右直交行列Ｖとなる。なお、特異値σ_ｊが０でない場合には、左直交行列Ｕの各列を構成する左特異ベクトルｙ_ｊを次のようにして求めることができる。

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

を解くことにより、ｙ_ｊを算出することができる。ここで、次の置き換えを行うことによって、左特異ベクトルｙ_ｊを正規化する。

このようにして、左特異ベクトルｙ_ｊを求めることができ、この左特異ベクトルｙ_ｊを用いて、Ｕ＝（ｙ_１，ｙ_２，…，ｙ_ｍ）とすることにより、左直交行列Ｕを算出することができる。

ここで、画像のデジタル誤差のため、ゼロでない特異値は３つ以上出てくる。しかし、４番目以降の特異値はノイズによるもので、最初の３つの特異値と比べて格段に小さい。そこで、ステップＳ５００２では、最初の３つの特異値に対して特異ベクトルを計算する。すなわち、本実施の形態による固有値分解装置１では、最初の３個の固有値に対して固有ベクトルを計算する。採用する３個のベクトルをまとめると、次式となる。

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

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

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

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

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

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

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

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

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

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

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

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

次に、ステップＳ６００３では、行列Ｄの特異値分解を行う。この特異値分解として、前述の２次元画像から３次元へ復元する画像処理への応用で説明したように、固有値分解による方法を用いることができる。すなわち、固有値分解装置１において、行列記憶部１１が記憶している対称行列Ａを、このたびの行列Ｄ^ＴＤとすることにより、前述の説明のようにして、対称行列Ａの固有値分解がなされ、その結果として、行列Ｄの特異値分解を行うことができうる。なお、この特異値分解では、計算された特異値のうち、大きい順に１番目からｋ番目までのｋ個の特異値に対してＤの特異ベクトルを算出する。すなわち、固有値分解において、計算された固有値のうち、大きい順に１番目からｋ番目までのｋ個の固有値に対して固有ベクトルを算出し、この固有ベクトルを用いて特異ベクトルを算出すればよい。ｋは、ステップＳ６００２で選択された値である。

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

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

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

ステップＳ６００５において、ステップＳ６００４において計算された類似度を基準に、検索結果を出力する。
このように、本実施の形態による固有値分解装置１は、画像処理や検索処理等の種々の処理に対して用いることができうる。なお、上記以外の処理に用いてもよいことは言うまでもない。

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

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

また、本実施の形態による固有値分解装置１における固有ベクトルの算出において、ある固有値を用いた固有ベクトルの算出では、ｑｄ型ツイスト分解法を用い、他の固有値を用いた固有ベクトルの算出では、ＬＶ型ツイスト分解法を用いるようにしてもよい。例えば、値の近接している固有値については、ＬＶ型ツイスト分解法を用いて固有ベクトルの算出を行い、値が近接していない固有値については、ｑｄ型ツイスト分解法を用いて固有ベクトルの算出を行ってもよい。固有値の値が近接しているかどうかは、所定の記録媒体等においてあらかじめ設定されているしきい値と比較することによって判断してもよい。

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

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

また、このプログラムでは、前記固有ベクトル算出ステップが、前記固有値記憶部から前記対称３重対角行列Ｔの固有値を読み出し、当該固有値のいずれかが負の値である場合に、前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、当該対称３重対角行列Ｔの固有ベクトルを変化させることなく、当該対称３重対角行列Ｔのすべての固有値が正の値となるように、当該対称３重対角行列Ｔを正定値化し、正定値化した後の固有値を前記固有値記憶部に蓄積し、正定値化した後の対称３重対角行列Ｔを前記対角行列記憶部に蓄積する正定値化ステップと、前記対角行列記憶部からすべての固有値が正の値である対称３重対角行列Ｔを読み出し、前記固有値記憶部から前記対称３重対角行列Ｔの正の値の固有値を読み出し、前記対称３重対角行列Ｔの各要素に関してｑｄ型変換を行うことによって、前記対称３重対角行列Ｔを上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）にコレスキー分解するコレスキー分解ステップと、前記上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）の各要素と、前記対称３重対角行列Ｔの正の値の固有値とを用いて固有ベクトルを算出して前記固有ベクトル記憶部に蓄積するベクトル算出ステップと、を備えてもよい。
また、このプログラムでは、前記固有ベクトル算出ステップにおいて、ＬＶ型ツイスト分解法により固有ベクトルを算出してもよい。

また、このプログラムでは、前記固有ベクトル算出ステップが、前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、前記固有値記憶部から前記対称３重対角行列Ｔの固有値を読み出し、前記対称３重対角行列Ｔの各要素に関してミウラ変換、ｄＬＶｖ型変換、逆ミウラ変換を行うことによって、前記対称３重対角行列Ｔを上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）にコレスキー分解するコレスキー分解ステップと、前記上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）の各要素と、前記対称３重対角行列Ｔの固有値とを用いて固有ベクトルを算出して前記固有ベクトル記憶部に蓄積するベクトル算出ステップと、を備えてもよい。

また、このプログラムでは、コンピュータに、対称行列Ａが記憶される行列記憶部から前記対称行列Ａを読み出し、前記対称行列Ａを３重対角化した前記対称３重対角行列Ｔを算出して前記対角行列記憶部に蓄積する対角化ステップをさらに実行させてもよい。
なお、上記プログラムにおいて、情報を蓄積するステップなどでは、ハードウェアでしか行われない処理は少なくとも含まれない。

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

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

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

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

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

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

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

プログラムは、コンピュータ１０１に、上記実施の形態による固有値分解装置１の機能を実行させるオペレーティングシステム（ＯＳ）、またはサードパーティプログラム等を必ずしも含まなくてもよい。プログラムは、制御された態様で適切な機能（モジュール）を呼び出し、所望の結果が得られるようにする命令の部分のみを含んでいてもよい。コンピュータシステム１００がどのように動作するのかについては周知であり、詳細な説明を省略する。
また、本発明は、以上の実施の形態に限定されることなく、種々の変更が可能であり、それらも本発明の範囲内に包含されるものであることは言うまでもない。
また、本発明のほんのいくつかの典型的な実施例について上で詳細に説明したが、その典型的な実施例において、発明の利益と新規な技術から実質的にはずれることなく多くの変更が可能であることを当業者は容易に認識することができるであろう。したがって、そのようなすべての変更は、本発明の範囲に含まれるものである。

以上のように、本発明による固有値分解装置等によれば、固有値分解を高速に処理することができ、分子軌道法による化学計算処理や統計計算処理、情報検索処理、その他の固有値分解を用いる処理を実行する装置等において有用である。

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

Claims

対称３重対角行列Ｔが記憶される対角行列記憶部と、
前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、当該対称３重対角行列Ｔを２個の対称３重対角行列に分割して前記対角行列記憶部に蓄積し、その対称３重対角行列を２個の対称３重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各対称３重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割部と、
前記あらかじめ決められた大きさ以下の各対称３重対角行列の固有値と、当該各対称３重対角行列の固有ベクトルからなる直交行列の一部の要素である行列要素とが記憶される固有値分解記憶部と、
前記あらかじめ決められた大きさ以下の各対称３重対角行列を前記対角行列記憶部から読み出し、前記各対称３重対角行列に対して固有値分解を行い、前記各対称３重対角行列の固有値と、前記各対称３重対角行列の固有ベクトルからなる直交行列の行列要素とを少なくとも算出し、当該固有値と当該行列要素とを前記固有値分解記憶部に蓄積する固有値分解部と、
前記対称３重対角行列Ｔの固有値が記憶される固有値記憶部と、
各対称３重対角行列の固有値と、行列要素とを前記固有値分解記憶部から読み出し、前記固有値と、前記行列要素とから、分割元の対称３重対角行列の固有値と、分割元の対称３重対角行列の行列要素とを算出して前記固有値分解記憶部に蓄積し、その分割元の対称３重対角行列の固有値と、行列要素とを算出する処理を、対称３重対角行列Ｔの少なくとも１個の固有値を算出するまで繰り返し、前記対称３重対角行列Ｔの少なくとも１個の固有値を前記固有値記憶部に蓄積する固有値算出部と、
前記対称３重対角行列Ｔの固有ベクトルが記憶される固有ベクトル記憶部と、
前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、前記固有値記憶部から前記対称３重対角行列Ｔの固有値を読み出し、前記対称３重対角行列Ｔとその固有値とから、ツイスト分解法を用いて前記対称３重対角行列Ｔの少なくとも１個の固有ベクトルを算出して前記固有ベクトル記憶部に蓄積する固有ベクトル算出部と、を備えた固有値分解装置。
前記固有ベクトル算出部は、ｑｄ型ツイスト分解法により固有ベクトルを算出する、請求項１記載の固有値分解装置。
前記固有ベクトル算出部は、
前記固有値記憶部から前記対称３重対角行列Ｔの固有値を読み出し、当該固有値のいずれかが負の値である場合に、前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、当該対称３重対角行列Ｔの固有ベクトルを変化させることなく、当該対称３重対角行列Ｔのすべての固有値が正の値となるように、当該対称３重対角行列Ｔを正定値化し、正定値化した後の固有値を前記固有値記憶部に蓄積し、正定値化した後の対称３重対角行列Ｔを前記対角行列記憶部に蓄積する正定値化部と、
前記対角行列記憶部からすべての固有値が正の値である対称３重対角行列Ｔを読み出し、前記固有値記憶部から前記対称３重対角行列Ｔの正の値の固有値を読み出し、前記対称３重対角行列Ｔの各要素に関してｑｄ型変換を行うことによって、前記対称３重対角行列Ｔを上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）にコレスキー分解するコレスキー分解部と、
前記上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）の各要素と、前記対称３重対角行列Ｔの正の値の固有値とを用いて固有ベクトルを算出して前記固有ベクトル記憶部に蓄積するベクトル算出部と、を備えた、請求項２記載の固有値分解装置。
前記固有ベクトル算出部は、ＬＶ型ツイスト分解法により固有ベクトルを算出する、請求項１記載の固有値分解装置。
前記固有ベクトル算出部は、
前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、前記固有値記憶部から前記対称３重対角行列Ｔの固有値を読み出し、前記対称３重対角行列Ｔの各要素に関してミウラ変換、ｄＬＶｖ型変換、逆ミウラ変換を行うことによって、前記対称３重対角行列Ｔを上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）にコレスキー分解するコレスキー分解部と、
前記上２重対角行列Ｂ^（＋１）及び下２重対角行列Ｂ^（−１）の各要素と、前記対称３重対角行列Ｔの固有値とを用いて固有ベクトルを算出して前記固有ベクトル記憶部に蓄積するベクトル算出部と、を備えた、請求項４記載の固有値分解装置。
前記コレスキー分解部は、複数のコレスキー分解手段を備え、
前記複数のコレスキー分解手段が、前記対称３重対角行列Ｔをコレスキー分解する処理を並列実行する、請求項３または請求項５記載の固有値分解装置。
前記ベクトル算出部は、複数のベクトル算出手段を備え、
前記複数のベクトル算出手段が、固有ベクトルを算出する処理を並列実行する、請求項３、請求項５または請求項６記載の固有値分解装置。
前記固有値算出部は、複数の固有値算出手段を備え、
前記複数の固有値算出手段が、分割元の対称３重対角行列の固有値と、行列要素とを算出する処理を並列実行する、請求項１から請求項７のいずれか記載の固有値分解装置。
前記固有値分解部は、複数の固有値分解手段を備え、
前記複数の固有値分解手段が、対称３重対角行列に対して固有値分解を行う処理を並列実行する、請求項１から請求項８のいずれか記載の固有値分解装置。
対称行列Ａが記憶される行列記憶部と、
前記対称行列Ａを前記行列記憶部から読み出し、前記対称行列Ａを３重対角化した前記対称３重対角行列Ｔを算出して前記対角行列記憶部に蓄積する対角化部と、をさらに備えた請求項１から請求項９のいずれか記載の固有値分解装置。
前記行列分割部は、対称３重対角行列を略半分の２個の対称３重対角行列に分割する、請求項１から請求項１０のいずれか記載の固有値分解装置。
前記固有値算出部は、対称３重対角行列Ｔの全ての固有値を算出する、請求項１から請求項１１のいずれか記載の固有値分解装置。
前記固有ベクトル算出部は、対称３重対角行列Ｔの全ての固有ベクトルを算出する、請求項１２記載の固有値分解装置。
前記対称行列Ａは、顔画像データを示すベクトルの平均から算出された共分散行列である、請求項１０記載の固有値分解装置。
対称３重対角行列Ｔが記憶される対角行列記憶部と、行列分割部と、固有値分解部と、前記固有値分解部によって固有値分解された固有値と固有ベクトルからなる直交行列の一部の要素である行列要素とが記憶される固有値分解記憶部と、前記対称３重対角行列Ｔの固有値が記憶される固有値記憶部と、固有値算出部と、前記対称３重対角行列Ｔの固有ベクトルが記憶される固有ベクトル記憶部と、固有ベクトル算出部とを備えた固有値分解装置で用いられる固有値分解方法であって、
前記行列分割部が、前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、当該対称３重対角行列Ｔを２個の対称３重対角行列に分割して前記対角行列記憶部に蓄積し、その対称３重対角行列を２個の対称３重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各対称３重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割ステップと、
前記固有値分解部が、前記あらかじめ決められた大きさ以下の各対称３重対角行列を前記対角行列記憶部から読み出し、前記各対称３重対角行列に対して固有値分解を行い、前記各対称３重対角行列の固有値と、前記各対称３重対角行列の固有ベクトルからなる直交行列の行列要素とを少なくとも算出し、当該固有値と当該行列要素とを前記固有値分解記憶部に蓄積する固有値分解ステップと、
前記固有値算出部が、各対称３重対角行列の固有値と、行列要素とを前記固有値分解記憶部から読み出し、前記固有値と、前記行列要素とから、分割元の対称３重対角行列の固有値と、分割元の対称３重対角行列の行列要素とを算出して前記固有値分解記憶部に蓄積し、その分割元の対称３重対角行列の固有値と、行列要素とを算出する処理を、対称３重対角行列Ｔの少なくとも１個の固有値を算出するまで繰り返し、前記対称３重対角行列Ｔの少なくとも１個の固有値を前記固有値記憶部に蓄積する固有値算出ステップと、
前記固有ベクトル算出部が、前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、前記固有値記憶部から前記対称３重対角行列Ｔの固有値を読み出し、前記対称３重対角行列Ｔとその固有値とから、ツイスト分解法を用いて前記対称３重対角行列Ｔの少なくとも１個の固有ベクトルを算出して前記固有ベクトル記憶部に蓄積する固有ベクトル算出ステップと、を備えた固有値分解方法。
コンピュータに、
対称３重対角行列Ｔが記憶される対角行列記憶部から前記対称３重対角行列Ｔを読み出し、当該対称３重対角行列Ｔを２個の対称３重対角行列に分割して前記対角行列記憶部に蓄積し、その対称３重対角行列を２個の対称３重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各対称３重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割ステップと、
前記あらかじめ決められた大きさ以下の各対称３重対角行列を前記対角行列記憶部から読み出し、前記各対称３重対角行列に対して固有値分解を行い、前記各対称３重対角行列の固有値と、前記各対称３重対角行列の固有ベクトルからなる直交行列の行列要素とを少なくとも算出し、当該固有値と当該行列要素とを固有値分解記憶部に蓄積する固有値分解ステップと、
各対称３重対角行列の固有値と、行列要素とを前記固有値分解記憶部から読み出し、前記固有値と、前記行列要素とから、分割元の対称３重対角行列の固有値と、分割元の対称３重対角行列の行列要素とを算出して前記固有値分解記憶部に蓄積し、その分割元の対称３重対角行列の固有値と、行列要素とを算出する処理を、対称３重対角行列Ｔの少なくとも１個の固有値を算出するまで繰り返し、前記対称３重対角行列Ｔの少なくとも１個の固有値を固有値記憶部に蓄積する固有値算出ステップと、
前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、前記固有値記憶部から前記対称３重対角行列Ｔの固有値を読み出し、前記対称３重対角行列Ｔとその固有値とから、ツイスト分解法を用いて前記対称３重対角行列Ｔの少なくとも１個の固有ベクトルを算出して固有ベクトル記憶部に蓄積する固有ベクトル算出ステップと、を実行させるためのプログラム。