(57)【要約】【課題】高速・高精度であり、並列処理可能な固有値分解装置を提供する。【解決手段】対称３重対角行列Ｔを２個の対称３重対角行列に分割する処理を繰り返し行う行列分割部１４と、分割後の対称３重対角行列を固有値分解する固有値分解部１５と、固有値分解部１５によって固有値分解された固有値と、固有ベクトルからなる直交行列の一部の要素である行列要素とから、分割元の対称３重対角行列の固有値と、分割元の対称３重対角行列の行列要素とを算出する処理を、対称３重対角行列Ｔの固有値を算出するまで繰り返す固有値算出部１７と、対称３重対角行列Ｔとその固有値とから、ツイスト分解法を用いて対称３重対角行列Ｔの固有ベクトルを算出する固有ベクトル算出部１９と、を備える。【選択図】図１
2008018188.tif 000002
【特許請求の範囲】
【請求項１】
対称３重対角行列Ｔが記憶される対角行列記憶部と、
前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、当該対称３重対角行列Ｔを２個の対称３重対角行列に分割して前記対角行列記憶部に蓄積し、その対称３重対角行列を２個の対称３重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各対称３重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割部と、
前記あらかじめ決められた大きさ以下の各対称３重対角行列の固有値と、当該各対称３重対角行列の固有ベクトルからなる直交行列の一部の要素である行列要素とが記憶される固有値分解記憶部と、
前記あらかじめ決められた大きさ以下の各対称３重対角行列を前記対角行列記憶部から読み出し、前記各対称３重対角行列に対して固有値分解を行い、前記各対称３重対角行列の固有値と、前記各対称３重対角行列の固有ベクトルからなる直交行列の行列要素とを少なくとも算出し、当該固有値と当該行列要素とを前記固有値分解記憶部に蓄積する固有値分解部と、
前記対称３重対角行列Ｔの固有値が記憶される固有値記憶部と、
各対称３重対角行列の固有値と、行列要素とを前記固有値分解記憶部から読み出し、前記固有値と、前記行列要素とから、分割元の対称３重対角行列の固有値と、分割元の対称３重対角行列の行列要素とを算出して前記固有値分解記憶部に蓄積し、その分割元の対称３重対角行列の固有値と、行列要素とを算出する処理を、対称３重対角行列Ｔの少なくとも１個の固有値を算出するまで繰り返し、前記対称３重対角行列Ｔの少なくとも１個の固有値を前記固有値記憶部に蓄積する固有値算出部と、
前記対称３重対角行列Ｔの固有ベクトルが記憶される固有ベクトル記憶部と、
前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、前記固有値記憶部から前記対称３重対角行列Ｔの固有値を読み出し、前記対称３重対角行列Ｔとその固有値とから、ツイスト分解法を用いて前記対称３重対角行列Ｔの少なくとも１個の固有ベクトルを算出して前記固有ベクトル記憶部に蓄積する固有ベクトル算出部と、を備えた固有値分解装置。
【請求項２】
前記固有ベクトル算出部は、ｑｄ型ツイスト分解法により固有ベクトルを算出する、請求項１記載の固有値分解装置。
【請求項３】
前記固有ベクトル算出部は、
前記固有値記憶部から前記対称３重対角行列Ｔの固有値を読み出し、当該固有値のいずれかが負の値である場合に、前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、当該対称３重対角行列Ｔの固有ベクトルを変化させることなく、当該対称３重対角行列Ｔのすべての固有値が正の値となるように、当該対称３重対角行列Ｔを正定値化し、正定値化した後の固有値を前記固有値記憶部に蓄積し、正定値化した後の対称３重対角行列Ｔを前記対角行列記憶部に蓄積する正定値化部と、
前記対角行列記憶部からすべての固有値が正の値である対称３重対角行列Ｔを読み出し、前記固有値記憶部から前記対称３重対角行列Ｔの正の値の固有値を読み出し、前記対称３重対角行列Ｔの各要素に関してｑｄ型変換を行うことによって、前記対称３重対角行列Ｔを上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）にコレスキー分解するコレスキー分解部と、
前記上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）の各要素と、前記対称３重対角行列Ｔの正の値の固有値とを用いて固有ベクトルを算出して前記固有ベクトル記憶部に蓄積するベクトル算出部と、を備えた、請求項２記載の固有値分解装置。
【請求項４】
前記固有ベクトル算出部は、ＬＶ型ツイスト分解法により固有ベクトルを算出する、請求項１記載の固有値分解装置。
【請求項５】
前記固有ベクトル算出部は、
前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、前記固有値記憶部から前記対称３重対角行列Ｔの固有値を読み出し、前記対称３重対角行列Ｔの各要素に関してミウラ変換、ｄＬＶｖ型変換、逆ミウラ変換を行うことによって、前記対称３重対角行列Ｔを上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）にコレスキー分解するコレスキー分解部と、
前記上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）の各要素と、前記対称３重対角行列Ｔの固有値とを用いて固有ベクトルを算出して前記固有ベクトル記憶部に蓄積するベクトル算出部と、を備えた、請求項４記載の固有値分解装置。
【請求項６】
前記コレスキー分解部は、複数のコレスキー分解手段を備え、
前記複数のコレスキー分解手段が、前記対称３重対角行列Ｔをコレスキー分解する処理を並列実行する、請求項３または請求項５記載の固有値分解装置。
【請求項７】
前記ベクトル算出部は、複数のベクトル算出手段を備え、
前記複数のベクトル算出手段が、固有ベクトルを算出する処理を並列実行する、請求項３、請求項５または請求項６記載の固有値分解装置。
【請求項８】
前記固有値算出部は、複数の固有値算出手段を備え、
前記複数の固有値算出手段が、分割元の対称３重対角行列の固有値と、行列要素とを算出する処理を並列実行する、請求項１から請求項７のいずれか記載の固有値分解装置。
【請求項９】
前記固有値分解部は、複数の固有値分解手段を備え、
前記複数の固有値分解手段が、対称３重対角行列に対して固有値分解を行う処理を並列実行する、請求項１から請求項８のいずれか記載の固有値分解装置。
【請求項１０】
対称行列Ａが記憶される行列記憶部と、
前記対称行列Ａを前記行列記憶部から読み出し、前記対称行列Ａを３重対角化した前記対称３重対角行列Ｔを算出して前記対角行列記憶部に蓄積する対角化部と、をさらに備えた請求項１から請求項９のいずれか記載の固有値分解装置。
【請求項１１】
前記行列分割部は、対称３重対角行列を略半分の２個の対称３重対角行列に分割する、請求項１から請求項１０のいずれか記載の固有値分解装置。
【請求項１２】
前記固有値算出部は、対称３重対角行列Ｔの全ての固有値を算出する、請求項１から請求項１１のいずれか記載の固有値分解装置。
【請求項１３】
前記固有ベクトル算出部は、対称３重対角行列Ｔの全ての固有ベクトルを算出する、請求項１２記載の固有値分解装置。
【請求項１４】
前記対称行列Ａは、顔画像データを示すベクトルの平均から算出された共分散行列である、請求項１０記載の固有値分解装置。
【請求項１５】
対称３重対角行列Ｔが記憶される対角行列記憶部と、行列分割部と、固有値分解部と、前記固有値分解部によって固有値分解された固有値と固有ベクトルからなる直交行列の一部の要素である行列要素とが記憶される固有値分解記憶部と、前記対称３重対角行列Ｔの固有値が記憶される固有値記憶部と、固有値算出部と、前記対称３重対角行列Ｔの固有ベクトルが記憶される固有ベクトル記憶部と、固有ベクトル算出部とを備えた固有値分解装置で用いられる固有値分解方法であって、
前記行列分割部が、前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、当該対称３重対角行列Ｔを２個の対称３重対角行列に分割して前記対角行列記憶部に蓄積し、その対称３重対角行列を２個の対称３重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各対称３重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割ステップと、
前記固有値分解部が、前記あらかじめ決められた大きさ以下の各対称３重対角行列を前記対角行列記憶部から読み出し、前記各対称３重対角行列に対して固有値分解を行い、前記各対称３重対角行列の固有値と、前記各対称３重対角行列の固有ベクトルからなる直交行列の行列要素とを少なくとも算出し、当該固有値と当該行列要素とを前記固有値分解記憶部に蓄積する固有値分解ステップと、
前記固有値算出部が、各対称３重対角行列の固有値と、行列要素とを前記固有値分解記憶部から読み出し、前記固有値と、前記行列要素とから、分割元の対称３重対角行列の固有値と、分割元の対称３重対角行列の行列要素とを算出して前記固有値分解記憶部に蓄積し、その分割元の対称３重対角行列の固有値と、行列要素とを算出する処理を、対称３重対角行列Ｔの少なくとも１個の固有値を算出するまで繰り返し、前記対称３重対角行列Ｔの少なくとも１個の固有値を前記固有値記憶部に蓄積する固有値算出ステップと、
前記固有ベクトル算出部が、前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、前記固有値記憶部から前記対称３重対角行列Ｔの固有値を読み出し、前記対称３重対角行列Ｔとその固有値とから、ツイスト分解法を用いて前記対称３重対角行列Ｔの少なくとも１個の固有ベクトルを算出して前記固有ベクトル記憶部に蓄積する固有ベクトル算出ステップと、を備えた固有値分解方法。
【請求項１６】
コンピュータに、
対称３重対角行列Ｔが記憶される対角行列記憶部から前記対称３重対角行列Ｔを読み出し、当該対称３重対角行列Ｔを２個の対称３重対角行列に分割して前記対角行列記憶部に蓄積し、その対称３重対角行列を２個の対称３重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各対称３重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割ステップと、
前記あらかじめ決められた大きさ以下の各対称３重対角行列を前記対角行列記憶部から読み出し、前記各対称３重対角行列に対して固有値分解を行い、前記各対称３重対角行列の固有値と、前記各対称３重対角行列の固有ベクトルからなる直交行列の行列要素とを少なくとも算出し、当該固有値と当該行列要素とを固有値分解記憶部に蓄積する固有値分解ステップと、
各対称３重対角行列の固有値と、行列要素とを前記固有値分解記憶部から読み出し、前記固有値と、前記行列要素とから、分割元の対称３重対角行列の固有値と、分割元の対称３重対角行列の行列要素とを算出して前記固有値分解記憶部に蓄積し、その分割元の対称３重対角行列の固有値と、行列要素とを算出する処理を、対称３重対角行列Ｔの少なくとも１個の固有値を算出するまで繰り返し、前記対称３重対角行列Ｔの少なくとも１個の固有値を固有値記憶部に蓄積する固有値算出ステップと、
前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、前記固有値記憶部から前記対称３重対角行列Ｔの固有値を読み出し、前記対称３重対角行列Ｔとその固有値とから、ツイスト分解法を用いて前記対称３重対角行列Ｔの少なくとも１個の固有ベクトルを算出して固有ベクトル記憶部に蓄積する固有ベクトル算出ステップと、を実行させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、固有値分解を行う固有値分解装置等に関する。
【背景技術】
【０００２】
固有値分解は、例えば、分子軌道法による化学計算、統計計算、情報検索等の非常に広い分野で利用されている。また、固有値分解の方法も知られている（例えば、非特許文献１参照）。
【非特許文献１】Ｊ．Ｊ．Ｍ．Ｃｕｐｐｅｎ，「Ａ  ｄｉｖｉｄｅ  ａｎｄ  ｃｏｎｑｕｅｒ  ｍｅｔｈｏｄ  ｆｏｒ  ｔｈｅ  ｓｙｍｍｅｔｒｉｃ  ｔｒｉｄｉａｇｏｎａｌ  ｅｉｇｅｎｐｒｏｂｌｅｍ」、Ｎｕｍｅｒｉｓｃｈｅ  Ｍａｔｈｅｍａｔｉｋ，Ｖｏｌ．３６、ｐ．１７７−１９５、１９８１年
【発明の開示】
【発明が解決しようとする課題】
【０００３】
近年、これらの応用分野におけるデータ量の増大などに伴って、高速・高精度な固有値分解が求められている。また、固有値分解を並列処理することができる固有値分解装置等の開発が望まれていた。【０００４】
本発明は、上記状況の下になされたものであり、並列性に優れた高速、高精度な固有値分解を実行可能な固有値分解装置等を提供することを目的とする。
【課題を解決するための手段】
【０００５】
上記目的を達成するため、本発明による固有値分解装置は、対称３重対角行列Ｔが記憶される対角行列記憶部と、前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、当該対称３重対角行列Ｔを２個の対称３重対角行列に分割して前記対角行列記憶部に蓄積し、その対称３重対角行列を２個の対称３重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各対称３重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割部と、前記あらかじめ決められた大きさ以下の各対称３重対角行列の固有値と、当該各対称３重対角行列の固有ベクトルからなる直交行列の一部の要素である行列要素とが記憶される固有値分解記憶部と、前記あらかじめ決められた大きさ以下の各対称３重対角行列を前記対角行列記憶部から読み出し、前記各対称３重対角行列に対して固有値分解を行い、前記各対称３重対角行列の固有値と、前記各対称３重対角行列の固有ベクトルからなる直交行列の行列要素とを少なくとも算出し、当該固有値と当該行列要素とを前記固有値分解記憶部に蓄積する固有値分解部と、前記対称３重対角行列Ｔの固有値が記憶される固有値記憶部と、各対称３重対角行列の固有値と、行列要素とを前記固有値分解記憶部から読み出し、前記固有値と、前記行列要素とから、分割元の対称３重対角行列の固有値と、分割元の対称３重対角行列の行列要素とを算出して前記固有値分解記憶部に蓄積し、その分割元の対称３重対角行列の固有値と、行列要素とを算出する処理を、対称３重対角行列Ｔの少なくとも１個の固有値を算出するまで繰り返し、前記対称３重対角行列Ｔの少なくとも１個の固有値を前記固有値記憶部に蓄積する固有値算出部と、前記対称３重対角行列Ｔの固有ベクトルが記憶される固有ベクトル記憶部と、前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、前記固有値記憶部から前記対称３重対角行列Ｔの固有値を読み出し、前記対称３重対角行列Ｔとその固有値とから、ツイスト分解法を用いて前記対称３重対角行列Ｔの少なくとも１個の固有ベクトルを算出して前記固有ベクトル記憶部に蓄積する固有ベクトル算出部と、を備えたものである。【０００６】
このような構成により、まず固有値を算出し、その固有値に基づいてツイスト分解法を用いて固有ベクトルを算出することによって、高速で高精度な固有値分解を実現することができうる。また、並列性にも優れている。さらに、全ての固有値や固有ベクトルを算出する必要がない場合には、必要な範囲で固有値や固有ベクトルを算出することができ、処理負荷を軽減することができる。【０００７】
また、本発明による固有値分解装置では、前記固有ベクトル算出部は、ｑｄ型ツイスト分解法により固有ベクトルを算出してもよい。
このような構成により、高速で高精度な固有値分解を実現することができうる。また、並列性にも優れている。【０００８】
また、本発明による固有値分解装置では、前記固有ベクトル算出部は、前記固有値記憶部から前記対称３重対角行列Ｔの固有値を読み出し、当該固有値のいずれかが負の値である場合に、前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、当該対称３重対角行列Ｔの固有ベクトルを変化させることなく、当該対称３重対角行列Ｔのすべての固有値が正の値となるように、当該対称３重対角行列Ｔを正定値化し、正定値化した後の固有値を前記固有値記憶部に蓄積し、正定値化した後の対称３重対角行列Ｔを前記対角行列記憶部に蓄積する正定値化部と、前記対角行列記憶部からすべての固有値が正の値である対称３重対角行列Ｔを読み出し、前記固有値記憶部から前記対称３重対角行列Ｔの正の値の固有値を読み出し、前記対称３重対角行列Ｔの各要素に関してｑｄ型変換を行うことによって、前記対称３重対角行列Ｔを上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）にコレスキー分解するコレスキー分解部と、前記上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）の各要素と、前記対称３重対角行列Ｔの正の値の固有値とを用いて固有ベクトルを算出して前記固有ベクトル記憶部に蓄積するベクトル算出部と、を備えてもよい。【０００９】
このような構成により、負の値の固有値が含まれる場合には、対称３重対角行列Ｔを正定値化してからｑｄ型ツイスト分解法を実行することができる。また、ＬＶ型ツイスト分解法よりも高速に計算することができうる。【００１０】
また、本発明による固有値分解装置では、前記固有ベクトル算出部は、ＬＶ型ツイスト分解法により固有ベクトルを算出してもよい。
このような構成により、数値安定的に固有ベクトルの算出を行うことができうる。【００１１】
また、本発明による固有値分解装置では、前記固有ベクトル算出部は、前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、前記固有値記憶部から前記対称３重対角行列Ｔの固有値を読み出し、前記対称３重対角行列Ｔの各要素に関してミウラ変換、ｄＬＶｖ型変換、逆ミウラ変換を行うことによって、前記対称３重対角行列Ｔを上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）にコレスキー分解するコレスキー分解部と、前記上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）の各要素と、前記対称３重対角行列Ｔの固有値とを用いて固有ベクトルを算出して前記固有ベクトル記憶部に蓄積するベクトル算出部と、を備えてもよい。
このような構成により、コレスキー分解の処理において、複数の補助変数を用いることによって、数値安定的に固有ベクトルの算出を行うことができうる。【００１２】
また、本発明による固有値分解装置では、前記コレスキー分解部は、複数のコレスキー分解手段を備え、前記複数のコレスキー分解手段が、前記対称３重対角行列Ｔをコレスキー分解する処理を並列実行してもよい。
このような構成により、コレスキー分解の処理を短時間で行うことができうる。【００１３】
また、本発明による固有値分解装置では、前記ベクトル算出部は、複数のベクトル算出手段を備え、前記複数のベクトル算出手段が、固有ベクトルを算出する処理を並列実行してもよい。
このような構成により、固有ベクトルを算出する処理を短時間で行うことができうる。【００１４】
また、本発明による固有値分解装置では、前記固有値算出部は、複数の固有値算出手段を備え、前記複数の固有値算出手段が、分割元の対称３重対角行列の固有値と、行列要素とを算出する処理を並列実行してもよい。
このような構成により、固有値を算出する処理を短時間で行うことができうる。【００１５】
また、本発明による固有値分解装置では、前記固有値分解部は、複数の固有値分解手段を備え、前記複数の固有値分解手段が、対称３重対角行列に対して固有値分解を行う処理を並列実行してもよい。
このような構成により、固有値分解する処理を短時間で行うことができうる。【００１６】
また、本発明による固有値分解装置では、対称行列Ａが記憶される行列記憶部と、前記対称行列Ａを前記行列記憶部から読み出し、前記対称行列Ａを３重対角化した前記対称３重対角行列Ｔを算出して前記対角行列記憶部に蓄積する対角化部と、をさらに備えてもよい。【００１７】
このような構成により、任意の対称行列Ａについて固有値を算出することができる。また、対称３重対角行列Ｔの固有ベクトルを用いることにより、対称行列Ａの固有ベクトルを算出することもできうる。【００１８】
また、本発明による固有値分解装置では、前記行列分割部は、対称３重対角行列を略半分の２個の対称３重対角行列に分割してもよい。
このような構成により、並列処理を適切に行うことができうる。【００１９】
また、本発明による固有値分解装置では、前記固有値算出部は、対称３重対角行列Ｔの全ての固有値を算出してもよい。【００２０】
また、本発明による固有値分解装置では、前記固有ベクトル算出部は、対称３重対角行列Ｔの全ての固有ベクトルを算出してもよい。【００２１】
また、本発明による固有値分解装置では、前記対称行列Ａは、顔画像データを示すベクトルの平均から算出された共分散行列であってもよい。
このような構成により、固有値分解装置を顔画像データの認識の処理に用いることができうる。
【発明の効果】
【００２２】
本発明による固有値分解装置等によれば、高速で高精度な固有値分解を行うことができる。また、並列性にも優れている。
【発明を実施するための最良の形態】
【００２３】
以下、本発明による固有値分解装置について、実施の形態を用いて説明する。なお、以下の実施の形態において、同じ符号を付した構成要素及びステップは同一または相当するものであり、再度の説明を省略することがある。【００２４】
（実施の形態１）
本発明の実施の形態１による固有値分解装置について、図面を参照しながら説明する。
図１は、本実施の形態による固有値分解装置１の構成を示すブロック図である。図１において、本実施の形態による固有値分解装置１は、行列記憶部１１と、対角化部１２と、対角行列記憶部１３と、行列分割部１４と、固有値分解部１５と、固有値分解記憶部１６と、固有値算出部１７と、固有値記憶部１８と、固有ベクトル算出部１９と、固有ベクトル記憶部２０とを備える。【００２５】
行列記憶部１１では、任意の対称行列Ａが記憶される。この対称行列Ａは、各要素が実数である実行列である。なお、対称行列Ａが記憶されているとは、対称行列Ａを示すデータが記憶されている、という意味である。後述する記憶部においても同様である。行列記憶部１１は、所定の記録媒体（例えば、半導体メモリや磁気ディスク、光ディスクなど）によって実現されうる。行列記憶部１１での記憶は、ＲＡＭ等における一時的な記憶でもよく、あるいは、長期的な記憶でもよい。行列記憶部１１に対称行列Ａが記憶される過程は問わない。例えば、記録媒体を介して対称行列Ａが行列記憶部１１で記憶されるようになってもよく、通信回線等を介して送信された対称行列Ａが行列記憶部１１で記憶されるようになってもよく、あるいは、キーボードやマウス等の入力デバイスを介して入力された対称行列Ａが行列記憶部１１で記憶されるようになってもよい。【００２６】
対角化部１２は、対称行列Ａを行列記憶部１１から読み出し、その読み出した対称行列Ａを３重対角化した対称３重対角行列Ｔを算出する。そして、対角化部１２は、その算出した対称３重対角行列Ｔを対角行列記憶部１３に蓄積する。対角化部１２は、例えば、ハウスホルダー（Ｈｏｕｓｅｈｏｌｄｅｒ）変換を必要なだけ繰り返し行う方法や、ランチョス（Ｌａｎｃｚｏｓ）法を用いた方法、その他の３重対角化法を用いて、対称行列Ａを３重対角化する。【００２７】
対角行列記憶部１３では、対称３重対角行列Ｔが記憶される。対角行列記憶部１３は、所定の記録媒体（例えば、半導体メモリや磁気ディスク、光ディスクなど）によって実現されうる。対角行列記憶部１３での記憶は、ＲＡＭ等における一時的な記憶でもよく、あるいは、長期的な記憶でもよい。【００２８】
行列分割部１４は、対角行列記憶部１３から対称３重対角行列Ｔを読み出し、その対称３重対角行列Ｔを２個の対称３重対角行列に分割して対角行列記憶部１３に蓄積する。行列分割部１４は、その対称３重対角行列を２個の対称３重対角行列に分割して対角行列記憶部１３に蓄積する処理を、分割後の各対称３重対角行列があらかじめ決められた大きさ以下となるまで再帰的に繰り返す。【００２９】
固有値分解部１５は、あらかじめ決められた大きさ以下の各対称３重対角行列を対角行列記憶部１３から読み出し、その各対称３重対角行列に対して固有値分解を行い、各対称３重対角行列の固有値と、その各対称３重対角行列の固有ベクトルを算出する。固有値分解部１５は、例えば、２分法と逆反復法を組み合わせた方法、ＭＲ＾３法、ＱＲ法等を用いて固有値分解を行ってもよい。ここで、ＭＲ＾３法による固有値分解を行う場合には、ＬＡＰＡＣＫ３．０（ｈｔｔｐ：／／ｗｗｗ．ｎｅｔｌｉｂ．ｏｒｇ／ｌａｐａｃｋ／）で提供されているＤＳＴＥＧＲを用いてもよい。また、ＱＲ法による固有値分解を行う場合には、ＬＡＰＡＣＫ３．０で提供されているＤＳＴＥＱＲを用いてもよい。これらの固有値分解の方法については、すでに公知であり、その詳細な説明を省略する。固有値分解部１５は、算出した固有値を固有値分解記憶部１６に蓄積する。また、固有値分解部１５は、算出した固有ベクトルからなる直交行列の一部の要素である行列要素も固有値分解記憶部１６に蓄積する。固有ベクトルからなる直交行列とは、固有ベクトルを各列とする直交行列のことである。行列要素の詳細については後述する。【００３０】
固有値分解記憶部１６では、固有値分解部１５によって固有値分解された固有値と、前述の行列要素とが記憶される。また、固有値分解記憶部１６では、固有値算出部１７によって算出された固有値や行列要素も記憶される。固有値分解記憶部１６は、所定の記録媒体（例えば、半導体メモリや磁気ディスク、光ディスクなど）によって実現されうる。固有値分解記憶部１６での記憶は、ＲＡＭ等における一時的な記憶でもよく、あるいは、長期的な記憶でもよい。【００３１】
固有値算出部１７は、固有値分解部１５によって算出された固有値と、行列要素とを固有値分解記憶部１６から読み出し、その固有値と行列要素とから、分割元の対称３重対角行列の固有値と、分割元の対称３重対角行列の行列要素とを算出して固有値分解記憶部１６に蓄積する。固有値算出部１７は、その分割元の対称３重対角行列の固有値と、行列要素とを算出する処理を、対称３重対角行列Ｔの固有値を算出するまで再帰的に繰り返す。そして、固有値算出部１７は、算出した対称３重対角行列Ｔの固有値を固有値記憶部１８に蓄積する。【００３２】
固有値記憶部１８では、対称３重対角行列Ｔの固有値が記憶される。固有値記憶部１８は、所定の記録媒体（例えば、半導体メモリや磁気ディスク、光ディスクなど）によって実現されうる。固有値記憶部１８での記憶は、ＲＡＭ等における一時的な記憶でもよく、あるいは、長期的な記憶でもよい。【００３３】
固有ベクトル算出部１９は、対角行列記憶部１３から対称３重対角行列Ｔを読み出し、固有値記憶部１８から対称３重対角行列Ｔの固有値を読み出す。そして、固有ベクトル算出部１９は、対称３重対角行列Ｔとその固有値とから、ツイスト分解法を用いて対称３重対角行列Ｔの固有ベクトルを算出して固有ベクトル記憶部２０に蓄積する。固有ベクトル算出部１９は、ＬＶ型ツイスト分解法により固有ベクトルを算出してもよく、あるいは、ｑｄ型ツイスト分解法により固有ベクトルを算出してもよい。なお、具体的なツイスト分解の手法については後述する。【００３４】
固有ベクトル記憶部２０では、対称３重対角行列Ｔの固有ベクトルが記憶される。固有ベクトル記憶部２０は、所定の記録媒体（例えば、半導体メモリや磁気ディスク、光ディスクなど）によって実現されうる。固有ベクトル記憶部２０での記憶は、ＲＡＭ等における一時的な記憶でもよく、あるいは、長期的な記憶でもよい。【００３５】
なお、行列記憶部１１、対角行列記憶部１３、固有値分解記憶部１６、固有値記憶部１８、固有ベクトル記憶部２０の任意の２以上の記憶部は、同一の記録媒体によって実現されてもよく、あるいは、別々の記録媒体によって実現されてもよい。前者の場合には、例えば、対称行列Ａの記憶されている領域が行列記憶部１１となり、対称３重対角行列Ｔ等の記憶されている領域が対角行列記憶部１３となる。【００３６】
また、行列記憶部１１、対角行列記憶部１３、固有値分解記憶部１６、固有値記憶部１８、固有ベクトル記憶部２０の各記憶部は、２以上の記録媒体から構成されてもよい。【００３７】
次に、本実施の形態による固有値分解装置１の動作について、図２のフローチャートを用いて説明する。
（ステップＳ１０１）対角化部１２は、行列記憶部１１で記憶されている対称行列Ａを読み出し、その行列Ａを３重対角化して対称３重対角行列Ｔを算出して対角行列記憶部１３に蓄積する。【００３８】
（ステップＳ１０２）行列分割部１４、固有値分解部１５、固有値算出部１７によって対称３重対角行列Ｔの固有値が算出され、固有値記憶部１８に蓄積される。この処理の詳細については後述する。【００３９】
（ステップＳ１０３）固有ベクトル算出部１９は、対角行列記憶部１３から対称３重対角行列Ｔを読み出し、固有値記憶部１８から対称３重対角行列Ｔの固有値を読み出し、対称３重対角行列Ｔの固有ベクトルを算出して固有ベクトル記憶部２０に蓄積する。この処理の詳細については後述する。【００４０】
このようにして、対称３重対角行列Ｔの固有値分解が終了する。ここで、対称行列Ａの固有値は対称３重対角行列Ｔの固有値に等しいため、対称行列Ａの固有値も算出されたことになる。また、後述するように、所定の変換を行うことによって、対称行列Ａの固有ベクトルも対称３重対角行列Ｔの固有ベクトルから容易に算出することができる。【００４１】
次に、図２のフローチャートのステップＳ１０１，Ｓ１０２の処理、すなわち、対称３重対角行列Ｔの固有値を算出するまでの処理について、より詳細に説明する。まず、図２のフローチャートのステップＳ１０２の処理について、図３のフローチャートを用いて説明する。【００４２】
（ステップＳ２０１）行列分割部１４は、対角行列記憶部１３から対称３重対角行列Ｔを読み出し、その対称３重対角行列Ｔを２個の対称３重対角行列に分割して対角行列記憶部１３に蓄積する。行列分割部１４は、その対称３重対角行列を２個の対称３重対角行列に分割する処理を、分割後の各対称３重対角行列があらかじめ決められた大きさ以下となるまで繰り返す。この処理の詳細については後述する。【００４３】
（ステップＳ２０２）固有値分解部１５は、対角行列記憶部１３で記憶されているあらかじめ決められた大きさ以下の各対称３重対角行列について、固有値分解を行う。固有値分解部１５は、固有値分解の結果である固有値と、固有ベクトルからなる直交行列の一部の要素である行列要素とを固有値分解記憶部１６に蓄積する。【００４４】
（ステップＳ２０３）固有値算出部１７は、対称３重対角行列の固有値と、行列要素とを固有値分解記憶部１６から読み出し、分割元の対称３重対角行列の固有値と、分割元の対称３重対角行列の行列要素とを算出して固有値分解記憶部１６に蓄積する。固有値算出部１７は、その分割元の対称３重対角行列の固有値と、行列要素とを算出する処理を、対称３重対角行列Ｔの固有値を算出するまで繰り返し、対称３重対角行列Ｔの固有値を固有値記憶部１８に蓄積する。このようにして、固有値を算出する処理が終了する。【００４５】
［対称行列Ａの対角化］
対称行列Ａは、例えば、ハウスホルダー変換等を用いて、以下に示されるように３重対角化することができることが知られている。ここで、対称３重対角行列Ｔは、行の数と列の数とが一致する正方行列である。【００４６】
【数１】
2008018188.tif 000003
【００４７】
したがって、対角化部１２は、上記のようにして、行列記憶部１１から対称行列Ａを読み出し、対称３重対角行列Ｔを算出することができる。その算出された対称３重対角行列Ｔは、対角行列記憶部１３で記憶される（ステップＳ１０１）。【００４８】
［対称３重対角行列Ｔの分割］
対称３重対角行列Ｔを分割する処理について説明する。次式で示されるように、ｎ×ｎの対称３重対角行列Ｔが与えられた場合に、それらを２個の対称３重対角行列Ｔ１，Ｔ２と、その他の要素（下記のθとβ）とに分けることができる。【００４９】
【数２】
2008018188.tif 000004
ここで、ｎ１は、１<ｎ１<ｎとなる整数である。例えば、８×８の対称３重対角行列Ｔを半分の大きさに分割する具体的な処理は、次のようになる。【００５０】
【数３】
2008018188.tif 000005
【００５１】
ただし、β＝ｔ８であり、θは計算誤差が少なくなるように任意に設定することができる。また、上記各行列において、明記されている以外の行列の要素は０である。なお、上記各行列の明記されている各要素は０以外であってもよく、０であってもよい。ここで、並列処理を適切に実行するためには、ｎ１を、ｎ／２を超えない最大の整数にとるか、あるいは、ｎ／２を下まわらない最小の整数にとることが好適である（このｎ１の値を用いて行列を２個の行列に分割する場合を、「行列を略半分の２個の行列に分割する」と呼ぶことにする）。【００５２】
【数４】
2008018188.tif 000006
本実施の形態では、ｎ１を上式のようにとるものとする。なお、ｎ１の値は、前述のように、１<ｎ１<ｎの範囲内で任意であることは言うまでもない。【００５３】
したがって、行列分割部１４は、上述のようにして、対角行列記憶部１３が記憶している対称３重対角行列Ｔを読み出し、その対称３重対角行列Ｔを２個の対称３重対角行列Ｔ１，Ｔ２と、２個の値β、θに分割することができ、その分割の処理を繰り返すことができる。行列分割部１４は、分割後の２個の対称３重対角行列Ｔ１，Ｔ２と、２個の値β、θとを、対角行列記憶部１３に蓄積する。なお、２個の値β、θは、どの分割で発生した値であるのかがわかるように対角行列記憶部１３で記憶されることが好適である。図４は、図３のフローチャートのステップＳ２０１における行列分割部１４による行列を分割する処理を示すフローチャートである。【００５４】
（ステップＳ３０１）行列分割部１４は、カウンタＩを「１」に設定する。
（ステップＳ３０２）行列分割部１４は、対角行列記憶部１３からＩ番目の分割をしていない対称３重対角行列を読み出し、その対称３重対角行列を２個の対称３重対角行列と、その他の要素とに分割する。そして、行列分割部１４は、分割した２個の対称３重対角行列と、その他の要素とを対角行列記憶部１３に蓄積する。【００５５】
（ステップＳ３０３）行列分割部１４は、Ｉ番目の分割を行っていない対称３重対角行列が、対角行列記憶部１３で記憶されているかどうか判断する。そして、Ｉ番目の分割を行っていない対称３重対角行列が、対角行列記憶部１３で記憶されている場合には、ステップＳ３０２に戻り、そうでない場合には、ステップＳ３０４に進む。【００５６】
（ステップＳ３０４）行列分割部１４は、Ｉ番目の分割を行った対称３重対角行列の大きさがあらかじめ決められている大きさ以下かどうか判断する。行列分割部１４は、例えば、目的とする行列の大きさ（例えば、２５×２５など）の記憶されている図示しない記録媒体から、その目的とする行列の大きさを読み出し、対角行列記憶部１３で記憶されている、Ｉ番目の分割後の対称３重対角行列がその大きさ以下であるかどうかを判断してもよい。そして、Ｉ番目の分割を行った対称３重対角行列の大きさがあらかじめ決められている大きさ以下である場合には、対称３重対角行列を分割する処理は終了となり、そうでない場合には、ステップＳ３０５に進む。
（ステップＳ３０５）行列分割部１４は、カウンタＩを１だけインクリメントする。そして、ステップＳ３０２に戻る。【００５７】
なお、この図４のフローチャートでは、上述のように、行列分割部１４が各行列を、略半分の２個の行列に分割するため、Ｉ番目の分割後の各対称３重対角行列の大きさがほとんど同じである場合について説明したが、行列分割部１４が各行列を、略半分の２個の行列に分割しない場合には、行列分割部１４は、分割後の各行列があらかじめ決められた大きさ以下となるように、分割を繰り返すものとする。【００５８】
また、図４のフローチャートでは、ステップＳ３０４において、行列分割部１４が、対角行列記憶部１３で記憶されている分割後の行列の大きさをあらかじめ決められた大きさと比較する処理を実行する場合について説明したが、これは一例であって、行列分割部１４は、ステップＳ３０４において、それ以外の処理を行ってもよい。例えば、上述のように、行列分割部１４が各行列を、略半分の２個の行列に分割する場合には、元の対称３重対角行列Ｔの大きさを知ることができれば、何番目の分割で目的とする行列の大きさとなるのかを知ることができる。したがって、Ｎ番目（Ｎは１以上の整数）の分割で目的とする行列の大きさとなる場合には、ステップＳ３０４において、ＩがＮであるかどうかを比較し、ＮでなければステップＳ３０５に進み、Ｎであれば一連の処理を終了するようにしてもよい。【００５９】
図５は、行列の分割について説明するための図である。まず、行列分割部１４は、１番目の分割として、対称３重対角行列Ｔを、対称３重対角行列Ｔ１と、対称３重対角行列Ｔ２とに分割する（ステップＳ３０１，Ｓ３０２）。対称３重対角行列Ｔは、１個しかないため、行列分割部１４は、１番目の分割をしていない行列がないと判断する（ステップＳ３０３）。また、対称３重対角行列Ｔ１等はあらかじめ決められた大きさ以下の行列でないとすると（ステップＳ３０４）、行列分割部１４は、２番目の分割として、対称３重対角行列Ｔ１を、対称３重対角行列Ｔ１１と、対称３重対角行列Ｔ１２とに分割する（ステップＳ３０５，Ｓ３０２）。この場合には、２番目の分割を行っていない対称３重対角行列Ｔ２が存在するため、行列分割部１４は、対称３重対角行列Ｔ２も、対称３重対角行列Ｔ２１と、対称３重対角行列Ｔ２２とに分割する（ステップＳ３０３，Ｓ３０２）。このようにして、分割後の各対称３重対角行列が目的とする大きさ以下になるまで、行列を分割する処理が繰り返される。なお、図５では、対称３重対角行列以外の要素（上述のβ、θ）については省略している。また、図５において、各対称３重対角行列は、正方行列である。【００６０】
［分割された行列の固有値分解］
行列Ｔｉが対称行列であるとすると、行列Ｔｉの固有値分解は、次のようになる。
【数５】
2008018188.tif 000007
【００６１】
ここで、Ｄｉは、行列Ｔｉと同じ大きさの対角行列である。Ｄｉの各対角成分はＴｉの固有値である。また、Ｑｉは、行列Ｔｉの固有ベクトルからなる直交行列である。【００６２】
したがって、固有値分解部１５は、対角行列記憶部１３から分割後の各対称３重対角行列を読み出し、上記のようにして、固有値分解を行う（ステップＳ２０２）。固有値分解の方法として、例えば、２分法と逆反復法を組み合わせた方法、ＭＲ＾３法、ＱＲ法等を用いてもよいことは前述の通りである。図５で示されるように対称３重対角行列の分割が行われた場合には、固有値分解部１５は、各対称３重対角行列Ｔ１１１，Ｔ１１２，Ｔ１２１，Ｔ１２２，…，Ｔ２２２に対して固有値分解を行うことになる。なお、固有値分解部１５は、固有値分解の結果得られた各固有値と、固有値分解の結果得られた直交行列の一部の要素である行列要素とを固有値分解記憶部１６に蓄積する。ここで、行列要素が、直交行列のどの要素を含むのかについては後述する。なお、固有値分解部１５が固有値と行列要素とを固有値分解記憶部１６に蓄積するのではなく、固有値と固有ベクトルからなる直交行列とを固有値分解記憶部１６に蓄積してもよい。その場合でも、固有ベクトルからなる直交行列を蓄積することによって、少なくとも、行列要素が蓄積されたことになる。【００６３】
なお、後述するように、固有値を算出するときに用いるために必要なのは、直交行列Ｑｉの全体ではなく、その直交行列Ｑｉの第１行と最終行である。したがって、この固有値分解において、直交行列Ｑｉの全体を求めるのではなく、直交行列Ｑｉの第１行と最終行を求めるようにしてもよい。そのためには、例えば、ＬＡＰＡＣＫ３．０で提供されているＱＲ法の固有値分解ルーチンであるＤＳＴＥＱＲにおいて、直交行列Ｑｉの算出における初期値として、単位行列Ｉを用いる代わりに、次式で示される行列Ｙを用いてもよい。
【数６】
2008018188.tif 000008
【００６４】
このようにすることで、ＤＳＴＥＱＲで直交行列Ｑｉの第１行と最終行を求めるための計算領域を削減することができる。直交行列Ｑｉそのものを算出する場合には、Ｏ（ｎ２）の作業領域が必要であるが、上記の行列Ｙを用いて直交行列Ｑｉの第１行と最終行を算出する場合には、Ｏ（ｎ）の作業領域でよいことになる。その結果、計算量も少なくなり、固有値分解の処理速度も向上する。このように、固有値分解部１５は、分割後の対称３重対角行列に対して、直交行列Ｑｉそのものを算出する固有値分解を行ってもよく、あるいは、直交行列の一部の成分（ここでは、第１行と最終行）を算出する処理を行ってもよい。この後者の処理も固有値分解と呼ぶことにする。このように、固有値分解部１５は、あらかじめ決められた大きさ以下の各対称３重対角行列に対して固有値分解を行い、各対称３重対角行列の固有値と、各対称３重対角行列の固有ベクトルからなる直交行列の行列要素（ここでは、第１行と最終行）とを少なくとも算出し、その固有値と、その行列要素とを固有値分解記憶部１６に蓄積してもよい。【００６５】
［固有値の算出］
まず、図６で示されるように、分割された２個の対称３重対角行列Ｔ１，Ｔ２から、分割元の対称３重対角行列Ｔ０の固有値等を算出する処理について説明する。対称３重対角行列Ｔ１，Ｔ２は、次のように固有値分解されていたとする。ここで、対称３重対角行列Ｔ１，Ｔ２は、両者共に正方行列であるとする。
【数７】
2008018188.tif 000009
【００６６】
この場合に、分割元の対称３重対角行列Ｔ０は、次のようになる。
【数８】
2008018188.tif 000010
【００６７】
上記の式において、β、θは、対称３重対角行列の分割の処理において説明した値である。ここで、行列（Ｄ１２＋ρｚｚＴ）の固有値分解が得られれば、行列Ｔ０の固有値分解が得られたことになる。なぜなら、行列（Ｄ１２＋ρｚｚＴ）の固有値分解Ｄ１２＋ρｚｚＴ＝ＱＤＱＴが求まると、次式のようになり、行列Ｔ０の固有値分解が完了するからである。
【数９】
2008018188.tif 000011
【００６８】
次に、行列（Ｄ１２＋ρｚｚＴ）の固有値を算出する方法について説明する。行列（Ｄ１２＋ρｚｚＴ）の固有値をλｉ、固有ベクトルをｑｉとすると、
【数１０】
2008018188.tif 000012
となる。【００６９】
（１）ｚの成分に０が含まれず、Ｄ１２の対角成分に同じ値が含まれない場合
まず、ｚの成分に０が含まれず、Ｄ１２の対角成分に同じ値が含まれない場合について説明する。その場合には、上記の式を、次に示すように順次、変形することができる。
【数１１】
2008018188.tif 000013
【００７０】
ここで、
【数１２】
2008018188.tif 000014
とすると、次の固有方程式を得ることができる。固有値λｉ（ｉ＝１，２，…，ｎ）は、その固有方程式ｎ個の解である。下記のｆ（λ）は、単調増加であり、かつｆ（λ）＝０の解の範囲が既知であるため、ニュートン法を用いて計算することができる。
【数１３】
2008018188.tif 000015
【００７１】
なお、Ｄ＝ｄｉａｇ（λ１，λ２，…，λｎ）となる。行列Ｑは、行列（Ｄ１２＋ρｚｚＴ）の固有ベクトルからなる直交行列であるため、Ｑ＝（ｑ１，ｑ２，…，ｑｎ）となる。したがって、固有ベクトルｑｉを求めると、行列Ｑを求めることができる。なお、固有ベクトルｑｉは、次式で与えられる。
【数１４】
2008018188.tif 000016
【００７２】
（２）ｚの成分に０が含まれる場合
ｚの成分に０が含まれる場合には、行列の次数を下げることができ、その次数を下げた行列を用いて固有値、及び固有ベクトルを求めることができる。この行列の次数を下げる操作をデフレーションと呼ぶ。【００７３】
ｚｊ＝０の場合には、次式で示すように、（ｄｊ，ｅｊ）が一組の固有値と固有ベクトルとなる。ただし、ｅｊは、ｊ番目の成分が１であり、その他の成分が０である単位ベクトルである。
【数１５】
2008018188.tif 000017
【００７４】
また、行列（Ｄ１２＋ρｚｚＴ）の固有方程式は、次式のようになる。
【数１６】
2008018188.tif 000018
【００７５】
したがって、ｚｊ＝０の場合の固有値は、ｄｊと、行列（Ｄ１２＋ρｚｚＴ）の第ｊ行、第ｊ列を除いた小行列について、上記式１の固有方程式を解いた値である。すなわち、固有値は、ｄｊと、次の固有方程式を解いた値となる。なお、次の固有方程式を解いた固有値を、固有値λｉ（ｉ＝１，２，…，ｊ−１，ｊ＋１，…，ｎ）とする。
【数１７】
2008018188.tif 000019
【００７６】
固有ベクトルは、行列（Ｄ１２＋ρｚｚＴ）の第ｊ行、第ｊ列を除いた小行列について上記式２で算出したｎ−１次元の固有ベクトルの第（ｊ−１）行と、第ｊ行との間に０を挿入した次式のｎ次元のベクトルと、上述のｅｊとなる。
【数１８】
2008018188.tif 000020
【００７７】
このように、ｚの成分に０が含まれる場合には、行列（Ｄ１２＋ρｚｚＴ）よりも次元の小さい行列について固有値分解を行うことになるため、計算を高速化することが可能であり、結果として、計算時間が短縮されうる。【００７８】
（３）Ｄ１２の対角成分に同じ値が含まれる場合
ｄｊ＋１＝ｄｊ＋２＝…＝ｄｋの場合にも、上記説明と同様のデフレーションを行うことによって、固有値と固有ベクトルとを算出することができる。【００７９】
【数１９】
2008018188.tif 000021
に対して、次の２個の直交行列を考える。
【数２０】
2008018188.tif 000022
【００８０】
すると、ｄｊ＋１＝ｄｊ＋２＝…＝ｄｋより、
【数２１】
2008018188.tif 000023
となる。したがって、次式が得られる。
【数２２】
2008018188.tif 000024
【００８１】
ここで、
【数２３】
2008018188.tif 000025
の固有値と固有ベクトルとは、上記のｚの成分に０が含まれるデフレーションの場合と同様にして求めることができる。また、
【数２４】
2008018188.tif 000026
は、符号の任意性を除いて一意に定まる。ここで、その算出方法について説明する。まず、次式を満たす直交行列について考える。【００８２】
【数２５】
2008018188.tif 000027
この直交行列を次式のようにとると、上記式を満たすことになる。
【数２６】
2008018188.tif 000028
【００８３】
ここで、上記の式の＊の値は、次式から求まる。
【数２７】
2008018188.tif 000029
【００８４】
したがって、
【数２８】
2008018188.tif 000030
となる。その結果、
【数２９】
2008018188.tif 000031
も、符号の任意性を除いて求めることができる。このようにして、ｄｊ＋１＝ｄｊ＋２＝…＝ｄｋの場合にも行列（Ｄ１２＋ρｚｚＴ）の固有値と固有ベクトルとを求めることができる。すなわち、固有値は、ｄｊ＋２，ｄｊ＋３，…，ｄｋと、次の固有方程式を解いた値となる。なお、次の固有方程式を解いた固有値を、固有値λｉ（ｉ＝１，２，…，ｊ＋１，ｋ＋１，…，ｎ）とする。
【数３０】
2008018188.tif 000032
【００８５】
固有ベクトルは、次のようになる。
【数３１】
2008018188.tif 000033
【００８６】
なお、上記（２）ｚの成分に０が含まれる場合、（３）Ｄ１２の対角成分に同じ値が含まれる場合において、ｚの成分が０であるかどうかの判断や、Ｄ１２の対角成分に同じ値が含まれているかどうかの判断は、微少な誤差を考慮して行ってもよい。例えば、次式が満たされる場合には、ｚの成分が０であると判断してもよく、また、Ｄ１２の対角成分が同じ値であると判断してもよい。
【数３２】
2008018188.tif 000034
【００８７】
なお、εは、マシン・イプシロンに定数（例えば、２，４，８などの実数であり、計算誤差を少なくするための値を選択することが好適である）を掛けた値を用いてもよく、他の適切な値であってもよい。【００８８】
このようにして、（１）ｚの成分に０が含まれず、Ｄ１２の対角成分に同じ値が含まれない場合、（２）ｚの成分に０が含まれる場合、（３）Ｄ１２の対角成分に同じ値が含まれる場合のそれぞれについて、固有値分解Ｄ１２＋ρｚｚＴ＝ＱＤＱＴが求まることになる。したがって、行列Ｑと行列Ｑ１２から行列Ｑ０を求めることができ、行列Ｄも求まっているため、分割元の対称３重対角行列Ｔ０の固有値分解を行うことができる。このように、対象となる行列を２個の副行列に分割し、その分割した後の各副行列について固有値分解を行い、その後に各副行列の固有値分解を用いて分割元の行列の固有値分解を行う方法は、分割統治法（Ｄｉｖｉｄｅ  ａｎｄ  Ｃｏｎｑｕｅｒ：Ｄ&Ｃ）と呼ばれている。【００８９】
分割統治法では、固有ベクトルまで算出することになり、その固有ベクトルを算出する処理において行列の計算をしなければならないため、非常に負荷の大きい処理となる。分割統治法において固有値分解をする場合には、例えば、計算時間の９５％程度が固有ベクトルを算出するベクトル更新の処理（例えば、上述のように、行列Ｑと行列Ｑ１２から行列Ｑ０を算出するために行列を掛け合わせる処理）に費やされることもありうる。【００９０】
一方、分割元の行列Ｔ０の固有値のみを算出する場合には、上記処理を簡略化することができる。次に、その方法について説明する。上記の結果より、次のようになる。【００９１】
【数３３】
2008018188.tif 000035
【００９２】
ここで、上記のように、ｆは行列Ｔ０を固有値分解した結果の直交行列の最初の行であり、ｌは行列Ｔ０を固有値分解した結果の直交行列の最後の行である。このｆとｌ、すなわち、直交行列の最初の行と、最後の行とが行列要素となる。なお、直交行列Ｑ１２について、後述する列の入れ替えを行った場合には、ｆ１２、ｌ１２はそれぞれ、その列を入れ替えた後の直交行列Ｑ１２の第１行、最終行となる。【００９３】
上記の結果から、図６で示されるように、行列Ｔ１，Ｔ２から分割元の行列Ｔ０を構成する場合に、行列Ｔ０について、
【数３４】
2008018188.tif 000036
をそれぞれ算出することができる。このような処理を繰り返すことにより、最終的に、対称行列Ａを３重対角化した対称３重対角行列Ｔの固有値を算出することができる。【００９４】
なお、上記の（３）Ｄ１２の対角成分に同じ値が含まれる場合には、式３，式４が次のようになる。
【数３５】
2008018188.tif 000037
【００９５】
ここで、上記式を計算する場合に、
【数３６】
2008018188.tif 000038
を先に計算して、その後に
【数３７】
2008018188.tif 000039
を掛けた方が、計算量が少なくなるため、その順番で計算するようにしてもよい。【００９６】
したがって、固有値算出部１７は、上述のようにして、固有値分解記憶部１６から固有値と行列要素とを読み出し、対角行列記憶部１３から行列の分割時に発生したその他の要素に関するβ、θを読み出し、それらを用いることによって、分割元の対称３重対角行列の固有値と行列要素とを算出する。そして、それらを算出する処理を繰り返すことによって、最終的に対称３重対角行列Ｔの固有値を算出する。図７は、図３のフローチャートのステップＳ２０３における固有値算出部１７が固有値を算出する処理を示すフローチャートである。【００９７】
（ステップＳ４０１）固有値算出部１７は、カウンタＪを「１」に設定する。
（ステップＳ４０２）固有値算出部１７は、Ｊ番目の固有値の算出は最後の固有値の算出であるかどうか判断する。ここで、最後の固有値の算出とは、対称３重対角行列Ｔの固有値を算出することである。そして、最後の固有値の算出である場合には、ステップＳ４０６に進み、そうでない場合には、ステップＳ４０３に進む。【００９８】
（ステップＳ４０３）固有値算出部１７は、分割元の対称３重対角行列の固有値等を算出する。この処理の詳細については後述する。
（ステップＳ４０４）固有値算出部１７は、Ｊ番目の固有値の算出において、分割元の全ての対称３重対角行列の固有値等を算出したかどうか判断する。そして、Ｊ番目の固有値の算出において、分割元の全ての対称３重対角行列の固有値等を算出した場合には、ステップＳ４０５に進み、そうでない場合には、ステップＳ４０３に戻る。【００９９】
（ステップＳ４０５）固有値算出部１７は、カウンタＪを１だけインクリメントする。そして、ステップＳ４０２に戻る。
（ステップＳ４０６）固有値算出部１７は、対称３重対角行列Ｔの固有値を算出し、その算出した固有値を固有値記憶部１８に蓄積する。このようにして、対称３重対角行列Ｔの固有値を算出する一連の処理は終了となる。【０１００】
図８は、図７のフローチャートにおけるステップＳ４０３の詳細な処理を示すフローチャートである。なお、この図８のフローチャートの処理を実行する前に、対角行列Ｄ１２の対角成分が昇順、あるいは降順となるように並べ替えておくことが好適である。なお、その場合には、直交行列Ｑ１２の列、及びｚの成分も、対角行列Ｄ１２の並び替えと同様に並び替えるものとする。具体的には、対角行列Ｄ１２の第ｉ成分と、第ｊ成分とを入れ替えた場合には、直交行列Ｑ１２の第ｉ列と、第ｊ列とを入れ替え、ｚの第ｉ行と、第ｊ行とを入れ替えればよい。【０１０１】
（ステップＳ５０１）固有値算出部１７は、分割元の行列の固有値を、式１等を用いて算出する。そして、固有値算出部１７は、算出した固有値を固有値分解記憶部１６に蓄積する。【０１０２】
（ステップＳ５０２）固有値算出部１７は、ステップＳ５０１で算出した固有値を用いて、式２等からｑｉを算出する。このｑｉを算出することにより、ｑｉを列ベクトルに有する行列Ｑを算出したことになる。固有値算出部１７は、算出した行列Ｑを図示しない記録媒体において一次記憶してもよい。【０１０３】
（ステップＳ５０３）固有値算出部１７は、ステップＳ５０２で算出した行列Ｑを用いて、分割元の行列に関する行列要素ｆ、ｌを算出する。この行列要素は、式３，式４で示されるものである。そして、固有値算出部１７は、算出した分割元の行列要素を固有値分解記憶部１６に蓄積する。このようにしてステップＳ４０３の処理は終了となる。【０１０４】
なお、図８のフローチャートにおいて、ｚの成分に０が含まれず、Ｄ１２の対角成分に同じ値が含まれない場合には、上記（１）で説明したようにして、固有値、固有ベクトルからなる行列Ｑを算出することができる。また、ｚの成分に０が含まれ、Ｄ１２の対角成分に同じ値が含まれない場合には、上記（２）で説明したようにして、固有値、固有ベクトルからなる行列Ｑを算出することができる。また、Ｄ１２の対角成分に同じ値が含まれる場合には、ｚの成分に０が含まれるかどうかに関わらず、上記（３）で説明したようにして、固有値、固有ベクトルからなる行列Ｑを算出することができる。上記（３）で説明したようにして固有値等を算出する場合には、前述のように、行列Ｑを算出するのではなく、行列Ｑを構成する行列を直接ｆ１２等に作用させることによって、ｆ、ｌを算出するようにしてもよい。【０１０５】
また、図８のフローチャートでは、行列Ｑから、行列要素ｆ、ｌを一度に算出する場合について示したが、行列要素ｆ、ｌをその成分ごとに計算してもよい。上記の式３，式４から明らかなように、ｆ１２、ｌ１２に固有ベクトルｑｉを掛けたものが、行列要素ｆ、ｌのｉ番目の成分となる。このような方法で行列要素ｆ、ｌを計算する場合には、固有ベクトル固有ベクトルｑｉのみをｉの値を順次、変えながら記憶することによってその計算を行うことができる。したがって、行列要素ｆ、ｌの計算において行列Ｑ全体を記憶する必要がなくなり、その計算で用いるメモリ上の作業領域を削減することが可能である。【０１０６】
図９は、固有値を算出する処理について説明するための図である。固有値分解記憶部１６では、固有値分解部１５によって行列Ｔ１１１，Ｔ１１２，…，Ｔ２２２が固有値分解された結果、すなわち、それらの固有値と、行列要素とが記憶されているものとする。まず、固有値算出部１７は、固有値等の１番目の算出を開始する（ステップＳ４０１，Ｓ４０２）。ここで、固有値等の１番目の算出とは、図９の一番下の行の行列の固有値と、行列要素とから、下から２番目の行の行列の固有値と、行列要素とを算出することである。固有値算出部１７は、行列Ｔ１１１と、行列Ｔ１１２との各固有値、及び行列Ｔ１１１と、行列Ｔ１１２との各行列要素を固有値分解記憶部１６から読み出す。また、行列Ｔ１１を行列Ｔ１１１と、行列Ｔ１１２とに分解したときに発生した２個の値β、θを対角行列記憶部１３から読み出す。そして、固有値算出部１７は、それらの値を用いて、行列Ｔ１１の固有値を算出して固有値分解記憶部１６に蓄積する（ステップＳ５０１）。次に、固有値算出部１７は、行列Ｔ１１の固有値を用いて直交行列Ｑを算出し（ステップＳ５０２）、その直交行列Ｑを用いて行列Ｔ１１の行列要素を算出して固有値分解記憶部１６に蓄積する（ステップＳ５０３）。【０１０７】
固有値等の１番目の算出はまだ終わっていないため（ステップＳ４０４）、固有値算出部１７は、上記説明と同様にして、行列Ｔ１２１、行列Ｔ１２２の固有値、行列要素等に基づいて、行列Ｔ１２の固有値、行列要素を算出する（ステップＳ４０２，Ｓ４０３）。このようにして、固有値等の１番目の算出が終了すると、固有値算出部１７は、固有値等の２番目の算出、すなわち、行列Ｔ１１、行列Ｔ１２の固有値、行列要素等に基づいて行列Ｔ１の固有値、行列要素の算出を行う（ステップＳ４０４，Ｓ４０５，Ｓ４０２，Ｓ４０３）。その後、固有値算出部１７は、同様にして、行列Ｔ２１、行列Ｔ２２の固有値、行列要素等に基づいて行列Ｔ２の固有値、行列要素の算出を行う（ステップＳ４０４，Ｓ４０３）。【０１０８】
固有値算出部１７が固有値等の２番目の算出を終了すると（ステップＳ４０４）、次の固有値等の算出は、対称３重対角行列Ｔの固有値を求める最後の処理となるため（ステップＳ４０５，Ｓ４０２）、固有値算出部１７は、固有値の算出のみを行い、その算出した固有値を固有値記憶部１８に蓄積する（ステップＳ４０６）。このようにして、固有値の算出が終了する。【０１０９】
なお、図９の行列Ｔ１、行列Ｔ２の固有値、行列要素等に基づいて行列Ｔの固有値を算出する場合には、行列要素ｆ１、ｌ２、あるいは行列要素ｆ２、ｌ１を用いるだけでよい。したがって、行列Ｔ１、行列Ｔ２の行列要素を算出する場合に、行列要素ｆ１、ｌ２、あるいは行列要素ｆ２、ｌ１を計算しなくてもよい。このようにすることで、計算量を削減することができ、処理スピードを向上させることが可能である。【０１１０】
次に、固有ベクトル算出部１９について説明する。本実施の形態では、固有ベクトル算出部１９が（１）ｑｄ型ツイスト分解法により固有ベクトルを算出する場合と、（２）ＬＶ型ツイスト分解法により固有ベクトルを算出する場合とについて説明する。【０１１１】
（１）ｑｄ型ツイスト分解法により固有ベクトルを算出する場合
図１０は、ｑｄ型ツイスト分解法により固有ベクトルを算出する場合の固有ベクトル算出部１９の構成を示す図である。図１０において、固有ベクトル算出部１９は、正定値化部２１と、コレスキー分解部２２と、ベクトル算出部２３とを備える。【０１１２】
正定値化部２１は、固有値記憶部１８から対称３重対角行列Ｔの固有値を読み出し、その固有値のいずれかが負の値であるかどうか判断する。そして、その固有値のいずれかが負の値である場合に、正定値化部２１は、対角行列記憶部１３から対称３重対角行列Ｔを読み出し、その対称３重対角行列Ｔの固有ベクトルを変化させることなく、対称３重対角行列Ｔのすべての固有値が正の値となるように、その対称３重対角行列Ｔを正定値化し、正定値化した後の固有値を固有値記憶部１８に蓄積し、正定値化した後の対称３重対角行列Ｔを対角行列記憶部１３に蓄積する。この正定値化の処理を行うのは、ｑｄ型ツイスト分解法により固有ベクトルを算出する場合には、対称３重対角行列Ｔの固有値がすべて正の値でなければならないからである。なお、対称３重対角行列Ｔを正定値化した行列も、対称３重対角行列Ｔと呼ぶことにする。また、この正定値化の処理が行われた場合には、コレスキー分解や固有ベクトルを算出する処理は、正定値化後の対称３重対角行列Ｔ、固有値に対して行われることになる。この処理の詳細については後述する。ここで、正定値化の行われた場合であっても、固有値記憶部１８には、正定値化の前の固有値が記憶されているものとする。その正定値化の前の固有値が、固有値分解装置１が最終的に求める固有値だからである。また、正定値化の行われた場合であって、ｑｄ型ツイスト分解法による固有ベクトルの算出のみを行う場合には、対角行列記憶部１３で記憶されている対称３重対角行列Ｔを、正定値化の後の対称３重対角行列Ｔで上書きしてもよい。一方、ＬＶ型ツイスト分解法による固有ベクトルの算出も行う場合には、対角行列記憶部１３には、正定値化の前の対称３重対角行列Ｔと、正定値化の後の対称３重対角行列Ｔとが記憶されることになる。【０１１３】
コレスキー分解部２２は、対角行列記憶部１３から、すべての固有値が正の値である対称３重対角行列Ｔを読み出し、固有値記憶部１８から、対称３重対角行列Ｔの正の値の固有値を読み出す。そして、コレスキー分解部２２は、対称３重対角行列Ｔの各要素に関してｑｄ型変換を行うことによって、対称３重対角行列Ｔを上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）にコレスキー分解する。この処理の詳細については後述する。なお、「対称３重対角行列Ｔを上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）にコレスキー分解する」ことには、対称３重対角行列Ｔの固有値に単位ベクトルを掛けたものを対称３重対角行列Ｔから引いた行列を、上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）にコレスキー分解することが含まれてもよい。【０１１４】
ベクトル算出部２３は、コレスキー分解部２２が算出した上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）の各要素と、対称３重対角行列Ｔの正の値の固有値とを用いて固有ベクトルを算出して固有ベクトル記憶部２０に蓄積する。この固有ベクトルを各列に有する行列が、固有値分解の直交行列となる。【０１１５】
次に、図２のフローチャートのステップＳ１０３の処理、すなわち、対称３重対角行列Ｔの固有ベクトルを算出する処理について、より詳細に説明する。まず、図２のフローチャートのステップＳ１０３の処理について、図１１のフローチャートを用いて説明する。【０１１６】
（ステップＳ６０１）正定値化部２１は、固有値記憶部１８から対称３重対角行列Ｔの固有値を読み出し、その固有値のいずれかが負の値であるかどうか判断する。そして、負の値がある場合には、ステップＳ６０２に進み、そうでない場合には、ステップＳ６０３に進む。【０１１７】
（ステップＳ６０２）正定値化部２１は、対角行列記憶部１３から対称３重対角行列Ｔを読み出し、その対称３重対角行列Ｔの固有ベクトルを変化させることなく、対称３重対角行列Ｔのすべての固有値が正の値となるように、その対称３重対角行列Ｔを正定値化し、正定値化した後の対称３重対角行列Ｔを対角行列記憶部１３に蓄積し、正定値化した後の固有値を固有値記憶部１８に蓄積する。【０１１８】
（ステップＳ６０３）コレスキー分解部２２は、対角行列記憶部１３から、すべての固有値が正の値である対称３重対角行列Ｔを読み出し、固有値記憶部１８から、対称３重対角行列Ｔの正の値の固有値を読み出す。なお、正定値化が行われた場合には、正定値化の後の対称３重対角行列Ｔ、及び固有値が読み出されることになる。そして、コレスキー分解部２２は、対称３重対角行列Ｔの各要素に関してｑｄ型変換を行うことによって、対称３重対角行列Ｔを上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）にコレスキー分解する。【０１１９】
（ステップＳ６０４）ベクトル算出部２３は、コレスキー分解部２２が算出した上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）の各要素と、対称３重対角行列Ｔの正の値の固有値とを用いて固有ベクトルを算出する。【０１２０】
（ステップＳ６０５）ベクトル算出部２３は、算出した固有ベクトルを正規化する。すなわち、ベクトル算出部２３は、算出した固有ベクトルのノルムを算出し、固有ベクトルを算出したノルムで割ったものを最終的な固有ベクトルとして固有ベクトル記憶部２０に蓄積する。このようにして、対称３重対角行列Ｔを固有値分解する処理は終了となる。
なお、この後に、行列記憶部１１が記憶している対称行列Ａの固有ベクトルを算出する処理が行われてもよいが、ここでは省略している。【０１２１】
また、固有ベクトル算出部１９が算出した固有ベクトルの精度を上げるために、図１２で示されるように、固有ベクトルに関する処理を実行してもよく、あるいは、実行しなくてもよい。すなわち、図示しない逆反復法処理部は、固有ベクトル記憶部２０から固有ベクトル算出部１９が算出した固有ベクトルを読み出し、その固有ベクトルに対して逆反復法の処理を実行し、その結果の固有ベクトルを固有ベクトル記憶部２０に蓄積する（ステップＳ７０１）。このように、逆反復法の処理を実行することによって、固有ベクトルの精度を向上させることができうる。次に、図示しないグラムシュミット処理部は、高精度が必要かどうか判断する（ステップＳ７０２）。この判断は、あらかじめ高精度が必要かどうか設定されている記録媒体等から、その設定を読み出すことによって判断してもよい。そして、図示しないグラムシュミット処理部は、高精度が必要な場合には、固有ベクトル記憶部２０から逆反復法の処理の実行された固有ベクトルを読み出し、その固有ベクトルに対してグラムシュミット法の処理を実行し、その結果の固有ベクトルを固有ベクトル記憶部２０に蓄積する（ステップＳ７０３）。このように、グラムシュミット法の処理を実行することによって、固有ベクトルの直交性を高めることができうる。なお、逆反復法や、グラムシュミット法については、すでに公知であり、その詳細な説明を省略する。【０１２２】
［対称３重対角行列Ｔの正定値化］
行列Ｔを次のように固有値分解できたとする。
【数３８】
2008018188.tif 000040
【０１２３】
その行列Ｔを次のようにａ倍してｓＩだけシフトしたとする。ここで、Ｉは行列Ｔと同じ大きさの単位行列であり、ａ，ｓは実数である。
【数３９】
2008018188.tif 000041
【０１２４】
すると、シフト後の固有値分解は、次のようになる。
【数４０】
2008018188.tif 000042
【０１２５】
このように、行列Ｔをａ倍してｓＩだけシフトした場合には、各固有値がａ倍されてｓだけシフトされるが、固有ベクトルは変化しない。したがって、この実数倍とシフトとによって、すべての固有値を正の値にしてから、ツイスト分解法を用いた固有ベクトルの算出を行えばよいことになる。【０１２６】
なお、算出された固有値に負の値が含まれる場合には、例えば、（Ａ）最小の固有値によりシフトを行ってもよく、（Ｂ）最大の固有値によりシフトを行ってもよい。ここで、固有値には、負の値と正の値とが含まれているものとする。固有値がすべて正の値であれば、シフトを行う必要はなく、固有値がすべて負の値であれば、単に行列Ｔを−Ｔとすれば、正定値化することができる。【０１２７】
（Ａ）最小の固有値によるシフト
最小の固有値をλｍｉｎ<０とすると、ａ＝１、ｓ＝−λｍｉｎ＋λｍａｘ×εとする。ここで、εとしては、マシン・イプシロンに定数を掛けた値を用いてもよく、他の適切な値であってもよい。このように、行列Ｔをシフトすることによって、行列Ｔのすべての固有値を正の値にすることができる。【０１２８】
（Ｂ）最大の固有値によるシフト
最大の固有値をλｍａｘ>０とすると、ａ＝−１、ｓ＝λｍａｘ−λｍｉｎ×εとする。この場合には、次のように行列Ｔが変換されることになる。
【数４１】
2008018188.tif 000043
【０１２９】
上式の右辺の括弧内では、シフトによってすべての固有値が負の値になるが、全体にマイナスを掛けることによって、正定値化することができている。ここで、εは、（Ａ）の場合と同様にマシン・イプシロンに定数（例えば、２，３，４などの実数であり、計算誤差を少なくするための値を選択することが好適である）を掛けた値を用いてもよく、他の適切な値であってもよい。【０１３０】
固有値の最小値の絶対値が、固有値の最大値の絶対値よりも小さい場合には、上記（Ａ）のシフトを行い、固有値の最小値の絶対値が、固有値の最大値の絶対値よりも大きい場合には、上記（Ｂ）のシフトを行うようにしてもよい。このようにすることで、シフトの量を少なくすることができ、丸め誤差を削減することが可能だからである。
なお、正定値化の後に固有ベクトルを算出する処理の詳細については、次のＬＶ型ツイスト分解法の場合と一緒に説明する。【０１３１】
（２）ＬＶ型ツイスト分解法により固有ベクトルを算出する場合
図１３は、ＬＶ型ツイスト分解法により固有ベクトルを算出する場合の固有ベクトル算出部１９の構成を示す図である。図１３において、固有ベクトル算出部１９は、コレスキー分解部３１と、ベクトル算出部３２とを備える。【０１３２】
コレスキー分解部３１は、対角行列記憶部１３から、対称３重対角行列Ｔを読み出し、固有値記憶部１８から、対称３重対角行列Ｔの固有値を読み出す。そして、コレスキー分解部３１は、対称３重対角行列Ｔの各要素に関してミウラ変換、ｄＬＶｖ型変換、逆ミウラ変換を行うことによって、対称３重対角行列Ｔを上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）にコレスキー分解する。この処理の詳細については後述する。なお、「対称３重対角行列Ｔを上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）にコレスキー分解する」ことには、対称３重対角行列Ｔの固有値に単位ベクトルを掛けたものを対称３重対角行列Ｔから引いた行列を、上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）にコレスキー分解することが含まれてもよい。【０１３３】
ベクトル算出部３２は、コレスキー分解部３１が算出した上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）の各要素と、対称３重対角行列Ｔの固有値とを用いて固有ベクトルを算出して固有ベクトル記憶部２０に蓄積する。この固有ベクトルを各列に有する行列が、固有値分解の直交行列となる。【０１３４】
次に、図２のフローチャートのステップＳ１０３の処理、すなわち、対称３重対角行列Ｔの固有ベクトルを算出する処理について、より詳細に説明する。まず、図２のフローチャートのステップＳ１０３の処理について、図１４のフローチャートを用いて説明する。【０１３５】
（ステップＳ８０１）コレスキー分解部３１は、対角行列記憶部１３から対称３重対角行列Ｔを読み出し、固有値記憶部１８から対称３重対角行列Ｔの固有値を読み出す。そして、コレスキー分解部３１は、対称３重対角行列Ｔの各要素に関してミウラ変換、ｄＬＶｖ型変換、逆ミウラ変換を行うことによって、対称３重対角行列Ｔを上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）にコレスキー分解する。【０１３６】
（ステップＳ８０２）ベクトル算出部３２は、上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）の各要素と、対称３重対角行列Ｔの固有値とを用いて固有ベクトルを算出する。【０１３７】
（ステップＳ８０３）ベクトル算出部３２は、算出した固有ベクトルを正規化する。すなわち、ベクトル算出部３２は、算出した固有ベクトルのノルムを算出し、固有ベクトルを算出したノルムで割ったものを最終的な固有ベクトルとして固有ベクトル記憶部２０に蓄積する。このようにして、対称３重対角行列Ｔを固有値分解する処理は終了となる。【０１３８】
なお、この後に、行列記憶部１１が記憶している対称行列Ａの固有ベクトルを算出する処理が行われてもよいが、ここでは省略している。また、固有ベクトル算出部１９が算出した固有ベクトルの精度を上げるために、図１２で示されるように、固有ベクトルに関する処理を実行してもよく、あるいは、実行しなくてもよいことは、ｑｄ型ツイスト分解法の場合と同様である。【０１３９】
［固有値からの固有ベクトルの算出］
固有値から固有ベクトルを算出する処理については、ｑｄ型ツイスト分解法を用いる方法と、ＬＶ型ツイスト分解法を用いる方法とについて説明する。図１１のフローチャートで示されるように、ｑｄ型ツイスト分解法を用いて固有ベクトルを算出する場合に、負の固有値があれば、対称３重対角行列Ｔを正定値化する処理をあらかじめ実行しておく（ステップＳ６０１，Ｓ６０２）。ｑｄ型ツイスト分解法を用いる場合には、この後の処理において、その正定値化の後の対称３重対角行列Ｔ、及び固有値を用いるものとする。【０１４０】
対角行列記憶部１３では、次のようなｍ×ｍの対称３重対角行列Ｔが記憶されているものとする。
【数４２】
2008018188.tif 000044
【０１４１】
また、固有値記憶部１８では、固有値λｉ（１≦ｉ≦ｍ）が記憶されているものとする。このような場合に、対称３重対角行列Ｔの固有ベクトルを求めることは、次の方程式の固有ベクトルｘｊを求めることになる。
【数４３】
2008018188.tif 000045
【０１４２】
本来であれば、上記式の右辺は０になるはずであるが、行列Ｔの固有値を算出するときに、いくらかの誤差を含むため、固有ベクトルｘｊが真値であれば、上記のように右辺に残差項が存在する。【０１４３】
ここで、以下のようにコレスキー分解できたとする。
【数４４】
2008018188.tif 000046
【０１４４】
ここで、行列Ｔの各要素と、行列Ｂ（０）の各要素との対応は、次のようになる。
【数４５】
2008018188.tif 000047
【０１４５】
なお、行列Ｔの各要素と、行列Ｂ（０）の各要素との対応を次のように書くこともできる。下記の式から明らかなように、行列Ｂ（０）の各要素を、行列Ｔの各要素を用いて順次、算出することができる。
【数４６】
2008018188.tif 000048
【０１４６】
上記のコレスキー分解した結果を、次式のように書くことができる。
【数４７】
2008018188.tif 000049
【０１４７】
そして、次のようになる。
【数４８】
2008018188.tif 000050
【０１４８】
ここで、行列Ｎｋをツイスト行列と呼ぶ。また、
【数４９】
2008018188.tif 000051
であるので、（Ｔ−λｊＩ）ｘｊ＝γｋｅｋは、
【数５０】
2008018188.tif 000052
となる。この簡単な式を解くことによって固有ベクトルを算出することができる。具体的には、
【数５１】
2008018188.tif 000053
のようにわずかな演算で固有ベクトルを算出することができる。なお、あるρ０について、Ｄρ０＋＝０あるいはＤρ０−＝０であったとしても、行列Ｂ（０）の（ρ，ρ）成分であるｂ２ρ−１（０）と、行列Ｂの（ρ，ρ＋１）成分であるｂ２ρ（０）を用いて、あるいは、行列Ｔの成分を用いて、
【数５２】
2008018188.tif 000054
と固有ベクトルを算出することができる（このように固有ベクトルを算出する処理を例外処理と呼ぶことにする）。なお、残差項のパラメータγｋ及びｋの値は、
【数５３】
2008018188.tif 000055
の絶対値が最小となるように決定する。したがって、上述のコレスキー分解を求め、ツイスト行列Ｎｋを求めることができれば、固有ベクトルを求めることができることになる。そこで、次にコレスキー分解について説明する。【０１４９】
図１５で示すように、Ｔ−λｊＩをコレスキー分解することは、Ｂ（０）に対応する｛ｑｋ（０），ｅｋ（０）｝から、Ｂ（±１）に対応する｛ｑｋ（±１），ｅｋ（±１）｝を求めることである。【０１５０】
（ｑｄ型ツイスト分解法）
まず、ｑｄ型ツイスト分解法で用いるｑｄ型変換について説明する（図１５参照）。
【数５４】
2008018188.tif 000056
この変換は、さらにｓｔａｔｉｏｎａｒｙ  ｑｄ  ｗｉｔｈ  ｓｈｉｆｔ（ｓｔｑｄｓ）変換
【数５５】
2008018188.tif 000057
及びｒｅｖｅｒｓｅ−ｔｉｍｅ  ｐｒｏｇｒｅｓｓｉｖｅ  ｑｄ  ｗｉｔｈ  ｓｈｉｆｔ（ｒｐｑｄｓ）変換
【数５６】
2008018188.tif 000058
に分けられる。固有値λｊが既知であれば反復的な計算が不要なため、計算量を少なく抑えることができるが、常に数値安定性と精度が高い方法ではない。それは、ｓｔｑｄｓ変換、ｒｐｑｄｓ変換は、共に減算による桁落ちが発生する可能性があるからである。例えば、ｓｔｑｄｓ変換において、ｑｋ（０）＋ｅｋ−１（０）−λｊ〜ｅｋ−１（＋１）であれば、ｑｋ（＋１）を求める際に、倍精度計算でもの有効数字がわずか１桁になることもある。その場合には、ｑｋ（０）ｅｋ（０）／ｑｋ（＋１）を計算すると誤差が生じる。つまり、ｅｋ（＋１）が精度よく計算できないことになる。また、ｑｋ＋１（＋１）を求めるのにｅｋ（＋１）が必要であり、ｅｋ（＋１）を求めるのにｑｋ（＋１）が必要であるといったように逐次的な計算が要求されるので、１箇所で発生した桁落ちによる誤差が波及し、さらなる誤差増大の可能性も秘めている。その結果、理論上はｑｋ（＋１）≠０であるが、誤差蓄積によりｑｋ（＋１）＝０となり、ｑｋ（０）ｅｋ（０）／ｑｋ（＋１）の計算においてオーバーフローが起こるといった数値不安定な状況も想定される。Ｂ（０）の成分｛ｂ２ｋ−１（０），ｂ２ｋ（０）｝が与えられる、すなわち｛ｑｋ（０），ｅｋ（０）｝が与えられると、λｊ及びｅｋ−１（＋１）が一意的に決まるので、この状況は避けることはできない。ｒｐｑｄｓ変換も同様の性質を持つため、実用的なレベルにまで達したとはいいがたい。ＬＡＰＡＣＫにおいてＦＯＲＴＲＡＮルーチンＤＳＴＥＧＲとして改良版が公開されているものの欠点は完全に解決されてはいない。【０１５１】
なお、このように対角行列Ｔの各要素を、行列Ｂ（０）の各要素に変換して、その行列Ｂ（０）の各要素についてｑｄ型変換を行う場合であっても、間接的に対角行列Ｔの各要素に関してｑｄ型変換を行うことになるため、「対角行列Ｔの各要素に関してｑｄ型変換を行う」と言うことにする。
また、ここではスタンダードなｑｄ型変換について説明したが、ｑｄ型変換はｄｉｆｆｅｒｅｎｔｉａｌ  ｑｄ型変換であってもよい。【０１５２】
（ＬＶ型ツイスト分解法）
次に、ＬＶ型ツイスト分解法で用いるミウラ変換、ｄＬＶｖ型変換、逆ミウラ変換について説明する（図１５参照）。
【数５７】
2008018188.tif 000059
【０１５３】
まず、ミウラ変換について説明する。この変換は、次のように示される。
【数５８】
2008018188.tif 000060
【０１５４】
次に、ｄＬＶｖ型変換について説明する。この変換は、さらにｓｔａｔｉｏｎａｒｙ  ｄｉｓｃｒｅｔｅ  Ｌｏｔｋａ−Ｖｏｌｔｅｒｒａ  ｗｉｔｈ  ｖａｒｉａｂｌｅ  ｓｔｅｐ−ｓｉｚｅ（ｓｔｄＬＶｖ）変換
【数５９】
2008018188.tif 000061
及びｒｅｖｅｒｓｅ−ｔｉｍｅ  ｄｉｓｃｒｅｔｅ  Ｌｏｔｋａ−Ｖｏｌｔｅｒｒａ  ｗｉｔｈ  ｖａｒｉａｂｌｅ  ｓｔｅｐ−ｓｉｚｅ（ｒｄＬＶｖ）変換
【数６０】
2008018188.tif 000062
に分けられる。ただし、
【数６１】
2008018188.tif 000063
を満たす範囲内でδ（±１）は任意である。【０１５５】
最後に、逆ミウラ変換について説明する。この変換は、次のように示される。
【数６２】
2008018188.tif 000064
【０１５６】
このようにして、ｑｄ型変換と同様に、コレスキー分解を行うことができる。ｑｄ型変換では見られない離散Ｌｏｔｋａ−Ｖｏｌｔｅｒｒａ型変換の大きな特徴は、任意パラメータを持つことである。すなわち、λｊ＝１／δ（０）−１／δ（±１）を満たす範囲でδ（ｎ）の値を任意に設定できる。δ（ｎ）を変動させると補助変数ｕｋ（ｎ）の値も変化するが、桁落ちによる誤差や数値不安定が発生するかどうかは事前に判定できる。この判定は、ｉｆ文によって実装されてもよい。この場合は、δ（ｎ）を再設定後に再度計算すればよい。また、ｕｋ（±１）が求まれば逆ミウラ変換によってｑｋ（±１）及びｅｋ（±１）が独立に計算されるので、誤差が伝播しないという性質を持つ。なお、逆ミウラ変換をミウラ変換、ミウラ変換を逆ミウラ変換と呼んでもよく、ｓｔｄＬＶｖ変換をｓｔＬＶｖ変換と呼んでもよく、ｒｄＬＶｖ変換をｒＬＶｖ変換と呼んでもよい。【０１５７】
なお、このように対角行列Ｔの各要素を、行列Ｂ（０）の各要素に変換して、その行列Ｂ（０）の各要素についてＬＶ型変換を行う場合であっても、間接的に対角行列Ｔの各要素に関してＬＶ型変換を行うことになるため、「対角行列Ｔの各要素に関してＬＶ型変換を行う」と言うことにする。【０１５８】
ここで、ＬＶ型ツイスト分解法による処理のより詳細な処理の一例について説明する。図１６〜図２１は、ＬＶ型ツイスト分解法による処理の一例を示すフローチャートである。【０１５９】
図１６は、コレスキー分解の全体の処理の一例を示す図である。
（ステップＳ９０１）コレスキー分解部３１は、ミウラ変換を行う。この処理の詳細については後述する。【０１６０】
（ステップＳ９０２）コレスキー分解部３１は、１／δ（±１）を１／δ（０）−λｊとする。
（ステップＳ９０３）コレスキー分解部３１は、後述するＰｒｏｃｅｄｕｒｅ１の処理を実行する。【０１６１】
（ステップＳ９０４）コレスキー分解部３１は、後述するＰｒｏｃｅｄｕｒｅ２の処理を実行する。
（ステップＳ９０５）コレスキー分解部３１は、ｑｋ（＋１），ｅｋ（＋１）がすでに算出されているかどうか判断する。そして、すでに算出されている場合には、コレスキー分解の一連の処理は終了となり、一方、算出されていない場合には、ステップＳ９０１に戻る。【０１６２】
図１７は、図１６のフローチャートにおけるステップＳ９０３の処理の詳細を示すフローチャートである。
（ステップＳ１００１）コレスキー分解部３１は、ｑｋ（＋１），ｅｋ（＋１）がすでに算出されているかどうか判断する。そして、すでに算出されている場合には、終了となり、算出されていない場合には、ステップＳ１００２に進む。
（ステップＳ１００２）コレスキー分解部３１は、ｓｔｄＬＶｖ変換等の処理を行う。この処理の詳細については後述する。【０１６３】
図１８は、図１６のフローチャートにおけるステップＳ９０４の処理の詳細を示すフローチャートである。
（ステップＳ１１０１）コレスキー分解部３１は、ｑｋ（−１），ｅｋ（−１）がすでに算出されているかどうか判断する。そして、すでに算出されている場合には、終了となり、算出されていない場合には、ステップＳ１１０２に進む。
（ステップＳ１１０２）コレスキー分解部３１は、ｒｄＬＶｖ変換等の処理を行う。この処理の詳細については後述する。【０１６４】
図１９は、図１６のフローチャートのステップＳ９０１の処理の詳細を示すフローチャートである。
（ステップＳ１２０１）コレスキー分解部３１は、１／δ（０）の値を決定する。この値は、前述のように任意に決定することができるため、行列を正定値化する値に決定する。例えば、コレスキー分解部３１は、１／δ（０）の値を次式のように決定する。
【数６３】
2008018188.tif 000065
【０１６５】
なお、なお、ε（０）は、マシン・イプシロンに定数（例えば、２，３，４などの実数であり、計算誤差を少なくするための値を選択することが好適である）を掛けた値を用いてもよく、他の適切な値であってもよい。その後、例えば、ステップＳ１２０３等で桁落ちの発生する可能性があると判断された場合に、コレスキー分解部３１は、１／δ（０）の値を決めるε（０）の値を、例えば、マシン・イプシロンに掛ける定数を１ずつ大きくしながら設定していってもよい。なお、１／δ（０）の値を決定する方法は、これに限定されないことは言うまでもない。このように、ＬＶ型ツイスト分解法では、この１／δ（０）の値を自由に設定することができるため、ｑｄ型ツイスト分解法の場合のように、行列Ｔの正定値化を行わなくてもよいことになる。なお、ｑｄ型ツイスト分解法の場合と同様に、行列Ｔの正定値化の処理を行ってから、ＬＶ型ツイスト分解法を実行してもよいことは言うまでもない。
（ステップＳ１２０２）コレスキー分解部３１は、β１をｑ１（０）−１／δ（０）に設定する。【０１６６】
（ステップＳ１２０３）コレスキー分解部３１は、β１の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップＳ１２０４に進み、そうでない場合には、ステップＳ１２０１に戻って、１／δ（０）の値を決定する処理を再度実行する。なお、このεも任意に定めることができる。例えば、εとしてマシン・イプシロンを用いてもよい。εの値を大きくすると精度が向上し、εの値を小さくすると精度が低下する。なお、このステップＳ１２０３で行っている処理は、桁落ちが発生する可能性について判断する処理である。β１の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。【０１６７】
（ステップＳ１２０４）コレスキー分解部３１は、β２をｑ１（０）／（１＋δ（０）ｕ０（０））に設定する。なお、前述のようにｕ０（０）＝０であるため、β２はｑ１（０）に設定されたことになる。【０１６８】
（ステップＳ１２０５）コレスキー分解部３１は、ｕ１（０）（１＋δ（０）ｕ０（０））をβ１に設定する。ここで、β１は、ステップＳ１２０２においてｑ１（０）−１／δ（０）に設定されているが、前述のようにｕ０（０）＝０であるため、β１は、ｑ１（０）／（１＋δ（０）ｕ０（０）−１／δ（０）と等しく、これにミウラ変換を行うとｕ１（０）＝ｕ２ｋ−１（０）｜ｋ＝１となるからである。なお、ｕ０（０）＝０から（１＋δ（０）ｕ０（０））＝１である。【０１６９】
（ステップＳ１２０６）コレスキー分解部３１は、カウンタｋを「１」に設定する。
（ステップＳ１２０７）コレスキー分解部３１は、α１をｅｋ（０）／β１に設定する。前述のように、β１はｕ２ｋ−１（０）となるから、ｅｋ（０）／β１にミウラ変換を行うと、α１はδ（０）ｕ２ｋ（０）と等しいことになる。【０１７０】
（ステップＳ１２０８）コレスキー分解部３１は、α２をα１＋１に設定する。ステップＳ１２０７の説明からわかるように、α２は１＋δ（０）ｕ２ｋ（０）と等しいことになる。【０１７１】
（ステップＳ１２０９）コレスキー分解部３１は、α２の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップＳ１２１０に進み、そうでない場合には、ステップＳ１２０１に戻って、１／δ（０）の値を決定する処理を再度実行する。なお、このステップＳ１２０９で行っている処理は、ステップＳ１２０３の処理と同様に、桁落ちが発生する可能性について判断する処理である。α２の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。【０１７２】
（ステップＳ１２１０）コレスキー分解部３１は、α１×β２を算出してｕ２ｋ（０）（１＋δ（０）ｕ２ｋ−１（０））に設定する。前述のように、α１はδ（０）ｕ２ｋ（０）と等しく、β２はｑｋ（０）／（１＋δ（０）ｕ２ｋ−２（０））であるため、α１×β２にミウラ変換を実行すると、ｕ２ｋ（０）（１＋δ（０）ｕ２ｋ−１（０））に等しくなるからである。【０１７３】
（ステップＳ１２１１）コレスキー分解部３１は、β２をｑｋ＋１（０）／α２に設定する。前述のように、α２は１＋δ（０）ｕ２ｋ（０）と等しいため、β２は、ｑｋ＋１（０）／（１＋δ（０）ｕ２ｋ（０））となり、β２におけるｋの値を１だけインクリメントしたことになる。【０１７４】
（ステップＳ１２１２）コレスキー分解部３１は、β１をβ２−１／δ（０）に設定する。前述のように、β２はｑｋ＋１（０）／（１＋δ（０）ｕ２ｋ（０））と等しいため、β１は、ｑｋ＋１（０）／（１＋δ（０）ｕ２ｋ（０））−１／δ（０）となり、ミウラ変換を実行すると、β１は、ｕ２ｋ＋１（０）となる。したがって、β１におけるｋの値を１だけインクリメントしたことになる。【０１７５】
（ステップＳ１２１３）コレスキー分解部３１は、β１の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップＳ１２１４に進み、そうでない場合には、ステップＳ１２０１に戻って、１／δ（０）の値を決定する処理を再度実行する。なお、このステップＳ１２１３で行っている処理は、ステップＳ１２０３の処理と同様に、桁落ちが発生する可能性について判断する処理である。β１の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。【０１７６】
（ステップＳ１２１４）コレスキー分解部３１は、α２×β１を算出してｕ２ｋ＋１（０）（１＋δ（０）ｕ２ｋ（０））に設定する。前述のように、α２は１＋δ（０）ｕ２ｋ（０）と等しく、β１は、ｕ２ｋ＋１（０）に等しいからである。
（ステップＳ１２１５）コレスキー分解部３１は、カウンタｋを１だけインクリメントする。【０１７７】
（ステップＳ１２１６）コレスキー分解部３１は、カウンタｋがｍに等しいかどうか判断する。そして、ｍに等しい場合には、一連の処理は終了となり、そうでない場合には、ステップＳ１２０７に戻る。【０１７８】
このようにして、ｕ２ｋ＋１（０）（１＋δ（０）ｕ２ｋ（０））と、ｕ２ｋ（０）（１＋δ（０）ｕ２ｋ−１（０））とが算出されることになる。これらの値は、コレスキー分解部３１が有する図示しないメモリ等において一時的に記憶されてもよい。【０１７９】
図２０は、図１７のフローチャートのステップＳ１００２の処理の詳細を示すフローチャートである。
（ステップＳ１３０１）コレスキー分解部３１は、α２を１＋δ（＋１）ｕ０（＋１）に設定する。なお、前述のようにｕ０（＋１）＝０であるため、α２は１に設定されたことになる。【０１８０】
（ステップＳ１３０２）コレスキー分解部３１は、β１をｕ１（０）（１＋δ（０）ｕ０（０））／（１＋δ（＋１）ｕ０（＋１））に設定する。ここで、ｕ１（０）（１＋δ（０）ｕ０（０））の値としては、ステップＳ１２０５で算出したものを用いる。なお、前述のようにｕ０（＋１）＝０である。また、このβ１の式にｓｔｄＬＶｖ変換を実行すると、β１は、ｕ１（＋１）に設定されたことになる。【０１８１】
（ステップＳ１３０３）コレスキー分解部３１は、カウンタｋを「１」に設定する。
（ステップＳ１３０４）コレスキー分解部３１は、α１をβ１＋１／δ（＋１）に設定する。【０１８２】
（ステップＳ１３０５）コレスキー分解部３１は、α１の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップＳ１３０６に進み、そうでない場合には、図１６のフローチャートのＰｒｏｃｅｄｕｒｅ２（ステップＳ９０４）に進む。なお、このステップＳ１３０５で行っている処理は、ステップＳ１２０３の処理と同様に、桁落ちが発生する可能性について判断する処理である。α１の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。【０１８３】
（ステップＳ１３０６）コレスキー分解部３１は、β２をｕ２ｋ（０）（１＋δ（０）ｕ２ｋ−１（０））／α１に設定する。ここで、ｕ２ｋ（０）（１＋δ（０）ｕ２ｋ−１（０））の値としては、ステップＳ１２１０で算出したものを用いる。また、このβ２の式にｓｔｄＬＶｖ変換を実行すると、β２は、δ（＋１）ｕ２ｋ（＋１）に設定されたことになる。【０１８４】
（ステップＳ１３０７）コレスキー分解部３１は、α１×α２を算出する。前述のように、α１はβ１＋１／δ（＋１）＝ｕ２ｋ−１（＋１）＋１／δ（＋１）に等しく、α２は１＋δ（＋１）ｕ２ｋ−２（＋１）に等しいため、α１×α２に逆ミウラ変換を実行するとｑｋ（＋１）に等しくなる。したがって、コレスキー分解部３１は、α１×α２を算出してｑｋ（＋１）に設定する。【０１８５】
（ステップＳ１３０８）コレスキー分解部３１は、α２を１＋β２に設定する。前述のように、β２はδ（＋１）ｕ２ｋ（＋１）と等しいため、α２は１＋δ（＋１）ｕ２ｋ（＋１）となる。したがって、α２におけるｋの値を１だけインクリメントしたことになる。【０１８６】
（ステップＳ１３０９）コレスキー分解部３１は、α２の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップＳ１３１０に進み、そうでない場合には、図１６のフローチャートのＰｒｏｃｅｄｕｒｅ２（ステップＳ９０４）に進む。なお、このステップＳ１３０９で行っている処理は、ステップＳ１２０３の処理と同様に、桁落ちが発生する可能性について判断する処理である。α２の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。【０１８７】
（ステップＳ１３１０）コレスキー分解部３１は、β１×β２を算出する。前述のように、β１はｕ２ｋ−１（＋１）に等しく、β２はδ（＋１）ｕ２ｋ（＋１）と等しいため、β１×β２に逆ミウラ変換を実行すると、ｅｋ（＋１）に等しくなる。したがって、コレスキー分解部３１は、β１×β２を算出してｅｋ（＋１）に設定する。【０１８８】
（ステップＳ１３１１）コレスキー分解部３１は、β１をｕ２ｋ＋１（０）（１＋δ（０）ｕ２ｋ（０））／α２に設定する。ここで、ｕ２ｋ＋１（０）（１＋δ（０）ｕ２ｋ（０））の値としては、ステップＳ１２１４で算出したものを用いる。また、このβ１の式にｓｔｄＬＶｖ変換を実行すると、β１は、ｕ２ｋ＋１（＋１）に設定されたことになる。したがって、β１におけるｋの値を１だけインクリメントしたことになる。【０１８９】
（ステップＳ１３１２）コレスキー分解部３１は、カウンタｋを１だけインクリメントする。
（ステップＳ１３１３）コレスキー分解部３１は、カウンタｋがｍに等しいかどうか判断する。そして、ｍに等しい場合には、ステップＳ１３１４に進み、そうでない場合には、ステップＳ１３０４に戻る。【０１９０】
（ステップＳ１３１４）コレスキー分解部３１は、α２×（β１＋１／δ（＋１））を算出する。これは、ステップＳ１３０４でα１を更新した後に、ステップＳ１３０７でα１×α２を計算することと等しい。したがって、コレスキー分解部３１は、α２×（β１＋１／δ（＋１））を算出してｑ１（＋１）に設定する。このようにして、ミウラ変換された結果にｓｔｄＬＶｖ変換と、逆ミウラ変換とを実行して、ｑｋ（＋１）、ｅｋ（＋１）を算出する処理が終了する。これらの値は、コレスキー分解部３１が有する図示しないメモリ等において一時的に記憶されてもよい。【０１９１】
図２１は、図１８のフローチャートのステップＳ１１０２の処理の詳細を示すフローチャートである。
（ステップＳ１４０１）コレスキー分解部３１は、β１をｕ２ｍ−１（０）（１＋δ（０）ｕ２ｍ−２（０））／（１＋δ（−１）ｕ２ｍ（−１））に設定する。ここで、ｕ２ｍ−１（０）（１＋δ（０）ｕ２ｍ−２（０））の値としては、ステップＳ１２１４で算出したものを用いる。なお、前述のようにｕ２ｍ（−１）＝０である。また、このβ１の式にｒｄＬＶｖ変換を実行すると、β１は、ｕ２ｍ−１（−１）に設定されたことになる。【０１９２】
（ステップＳ１４０２）コレスキー分解部３１は、カウンタｋを「ｍ−１」に設定する。
（ステップＳ１４０３）コレスキー分解部３１は、α１をβ１＋１／δ（−１）に設定する。【０１９３】
（ステップＳ１４０４）コレスキー分解部３１は、α１の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップＳ１４０５に進み、そうでない場合には、図１６のフローチャートのミウラ変換（ステップＳ９０１）に戻る。なお、このステップＳ１４０４で行っている処理は、ステップＳ１２０３の処理と同様に、桁落ちが発生する可能性について判断する処理である。α１の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。【０１９４】
（ステップＳ１４０５）コレスキー分解部３１は、β２をｕ２ｋ（０）（１＋δ（０）ｕ２ｋ−１（０））／α１に設定する。ここで、ｕ２ｋ（０）（１＋δ（０）ｕ２ｋ−１（０））の値としては、ステップＳ１２１０で算出したものを用いる。また、このβ２の式にｒｄＬＶｖ変換を実行すると、β２は、δ（−１）ｕ２ｋ（−１）に設定されたことになる。【０１９５】
（ステップＳ１４０６）コレスキー分解部３１は、α２を１＋β２に設定する。前述のように、β２はδ（−１）ｕ２ｋ（−１）と等しいため、α２は１＋δ（−１）ｕ２ｋ（−１）となる。【０１９６】
（ステップＳ１４０７）コレスキー分解部３１は、α２の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップＳ１４０８に進み、そうでない場合には、図１６のフローチャートのミウラ変換（ステップＳ９０１）に戻る。なお、このステップＳ１４０７で行っている処理は、ステップＳ１２０３の処理と同様に、桁落ちが発生する可能性について判断する処理である。α２の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。【０１９７】
（ステップＳ１４０８）コレスキー分解部３１は、α１×α２を算出する。前述のように、α１はβ１＋１／δ（−１）＝ｕ２ｋ＋１（−１）＋１／δ（−１）に等しく、α２は１＋δ（−１）ｕ２ｋ（−１）に等しいため、α１×α２に逆ミウラ変換を実行するとｑｋ＋１（−１）に等しくなる。したがって、コレスキー分解部３１は、α１×α２を算出してｑｋ＋１（−１）に設定する。【０１９８】
（ステップＳ１４０９）コレスキー分解部３１は、β１をｕ２ｋ−１（０）（１＋δ（０）ｕ２ｋ−２（０））／α２に設定する。ここで、ｕ２ｋ−１（０）（１＋δ（０）ｕ２ｋ−２（０））の値としては、ステップＳ１２０５，Ｓ１２１４で算出したものを用いる。また、このβ１の式にｓｔｄＬＶｖ変換を実行すると、β１は、ｕ２ｋ−１（−１）に設定されたことになる。したがって、β１におけるｋの値を１だけデクリメントしたことになる。【０１９９】
（ステップＳ１４１０）コレスキー分解部３１は、β１×β２を算出する。前述のように、β１はｕ２ｋ−１（−１）に等しく、β２はδ（−１）ｕ２ｋ（−１）と等しいため、β１×β２に逆ミウラ変換を実行すると、ｅｋ（−１）に等しくなる。したがって、コレスキー分解部３１は、β１×β２を算出してｅｋ（−１）に設定する。【０２００】
（ステップＳ１４１１）コレスキー分解部３１は、カウンタｋを１だけデクリメントする。
（ステップＳ１４１２）コレスキー分解部３１は、カウンタｋが０に等しいかどうか判断する。そして、０に等しい場合には、ステップＳ１４１３に進み、そうでない場合には、ステップＳ１４０３に戻る。【０２０１】
（ステップＳ１４１３）コレスキー分解部３１は、β１＋１／δ（−１）を算出する。これは、ステップＳ１４０３でα１を更新した後に、ステップＳ１４０８でα１×α２を計算することと等しい。この場合、α２は１だからである。したがって、コレスキー分解部３１は、β１＋１／δ（−１）を算出してｑ１（−１）に設定する。このようにして、ミウラ変換された結果にｒｄＬＶｖ変換と、逆ミウラ変換とを実行して、ｑｋ（−１）、ｅｋ（−１）を算出する処理が終了する。これらの値は、コレスキー分解部３１が有する図示しないメモリ等において一時的に記憶されてもよい。【０２０２】
なお、図１６のフローチャートの処理によって算出することができるのは１個の固有値に対応する固有ベクトルであるため、全ての固有値に対応する固有ベクトルを算出する場合には、コレスキー分解部３１は、図１６のフローチャートの処理を固有値の数だけ繰り返すことになる。【０２０３】
メモリ消費を抑えるために、補助変数ｕｋ（ｎ）のための配列は必ずしも用意する必要はない。一方、１＋δ（０）ｕ（０）のためのメモリ領域を確保し、ミウラ変換、ｄＬＶｖ型変換、逆ミウラ変換のステップにまたがってこの値を利用することで、メモリ消費を抑え、計算量を低減することができる。これにより、誤差も低減される。【０２０４】
ここで、誤差について説明する。人間が理想的な状況、すなわち、無限桁の計算をいくらでもできるとすると、ｑｄ型ツイスト分解法であっても、ＬＶ型ツイスト分解法であっても問題ない。しかしながら、コンピュータで計算を行う場合は注意が必要である。有限桁の計算しか行えないコンピュータ上では、数学的に正しい計算法を使ったとしても必ずしも正しい結果が得られる訳ではない。そればかりか、いつまでたっても計算が終了しないといった思わぬ数値的な問題が発生する場合もある。【０２０５】
コンピュータ計算による誤差には、丸め誤差及び桁落ちによる誤差などが知られている。丸め誤差単独では、高々有効桁の最後の桁が真値と比べて異なる程度で大きな誤差にはならない。また、指数部が異なる２つの実数の加算、乗算、除算の少なくとも１つの演算を行えばやはり丸め誤差が生じるが、それ以上の誤差は発生しない。さらに、このような丸め誤差を発生するような操作が繰り返されても、丸めモードがｎｅａｒ（四捨五入）ならば、一方的に切り上げられたり、あるいは切り捨てられたりして誤差が極端に蓄積することは少ない。よって、多くの数値計算法は加算、乗算、除算の少なくとも１つの演算によって発生する丸め誤差を特別注意することは少なく、適切な固有値計算でも結果的に丸め誤差は一様に増大しない。【０２０６】
問題となるのは、同符号の実数の減算及び異符号の実数の加算により生じる、桁落ちである。桁落ちによる誤差で値が０となった後、その値による除算を行うと、０が分母にくるような不定形となり計算不可能となる。こうなるといつまでたっても計算が終了しない。減算→除算と計算する部分がｑｄ型ツイスト分解法、ＬＶ型ツイスト分解法の両方に存在するので、桁落ち誤差には十分に注意する必要がある。【０２０７】
ＬＶ型ツイスト分解法では、上述の桁落ちによる誤差を含んでいるかどうかは減算によって得られた値が小さいかどうかで判断できる。ｑｄ型ツイスト分解法の場合、桁落ち誤差を含むことが分かったとしても、それを回避することはできない。なぜならば、初期値として｛ｑｋ（０），ｅｋ（０）｝が与えられると、λｊは一意的に決定され、｛ｑｋ（±１），ｅｋ（±１）｝も一意的に導出されるためである。すなわち、任意パラメータを持たない自由度のない計算法であるためである。【０２０８】
それに対して、ＬＶ型ツイスト分解法は、自由に設定できるパラメータδ（０）を持つため、補助変数ｕｋ（ｎ）の値を様々に変化させることができる（図２２参照）。すなわち、様々な経路で｛ｑｋ（±１），ｅｋ（±１）｝を計算することができる。よって、桁落ちが発生する場合も回避できる。図１９〜図２１の条件判定によって桁落ちの影響をチェックし、減算によって得られた値の絶対値が小さな数εより大きいという条件が満たされなければ、パラメータδ（０）の設定に戻るというものである。この処理は、条件が満たされるまで繰り返される。なお、精度よりも高速性を重視する場合は、数回条件が満たさなければ（ｑｋ（＋１）＝０あるいはｑｋ（−１）＝０ならば）、例外処理を行ってもよい。【０２０９】
なお、前述の行列Ｔの各要素と、行列Ｂ（０）の各要素との対応を示す式を用いると、ｑｄ型変換を次のように書き換えることもできる。
【数６４】
2008018188.tif 000066
【０２１０】
この変換は、さらにｓｔｑｄｓ変換
【数６５】
2008018188.tif 000067
及びｒｐｑｄｓ変換
【数６６】
2008018188.tif 000068
に分けられる。このように、｛ｔ２ｋ−１，ｔ２ｋ｝から｛ｑｋ（±１），ｅｋ（±１）｝を求める場合には、｛ｔ２ｋ−１，ｔ２ｋ｝から｛ｑｋ（０），ｅｋ（０）｝を求め、さらにｑｋ（０），ｅｋ（０）｝から｛ｑｋ（±１），ｅｋ（±１）｝を求める場合に比べて、計算量が少なくなることとなり、計算速度を向上させることも可能である。【０２１１】
なお、この場合には、対角行列Ｔの各要素について直接ｑｄ型変換が行われることになる。したがって、この変換が、「対角行列Ｔの各要素に関してｑｄ型変換を行う」場合に含まれることは言うまでもない。【０２１２】
また、ｑｄ型変換と同様に、前述の行列Ｔの各要素と、行列Ｂ（０）の各要素との対応を示す式を用いると、ＬＶ型変換におけるミウラ変換を次のように書き換えることもできる。
【数６７】
2008018188.tif 000069
【０２１３】
この場合のミウラ変換は、次のように示される。
【数６８】
2008018188.tif 000070
このように、｛ｔ２ｋ−１，ｔ２ｋ｝から｛ｕ２ｋ−１（０），ｕ２ｋ（０）｝を求める場合には、｛ｔ２ｋ−１，ｔ２ｋ｝から｛ｑｋ（０），ｅｋ（０）｝を求め、さらにｑｋ（０），ｅｋ（０）｝から｛ｕ２ｋ−１（０），ｕ２ｋ（０）｝を求める場合に比べて、計算量が少なくなることとなり、計算速度を向上させることも可能である。【０２１４】
なお、この場合には、対角行列Ｔの各要素について直接ＬＶ型変換が行われることになる。したがって、この変換が、「対角行列Ｔの各要素に関してＬＶ型変換を行う」場合に含まれることは言うまでもない。【０２１５】
ｑｄ型ツイスト分解法あるいはＬＶ型ツイスト分解法によってコレスキー分解をすることができると、上述のようにツイストされた行列Ｎｋを算出することができ、その行列ＮｋのＮｋＴｘｊ＝ｅｋを解くことによって、ｘｊを算出することができる。ここで、
【数６９】
2008018188.tif 000071
とｘｊを置き換えることにより、このｘｊを正規化する。このようにして、固有ベクトルｘｊを求めることができる。この固有ベクトルｘｊを用いて、ＱＴ＝（ｘ１，ｘ２，...，ｘｍ）とすることにより、直交行列ＱＴを算出することができる。【０２１６】
したがって、次式のように対称３重対角行列Ｔについて固有値分解が行われたことになる。
【数７０】
2008018188.tif 000072
【０２１７】
また、対称３重対角行列Ｔの固有値分解を行うことができれば、次のように、対称行列Ａの固有値分解も行われる。
【数７１】
2008018188.tif 000073
【０２１８】
固有ベクトル算出部１９は、上述のようにして、固有値記憶部１８が記憶している対称３重対角行列Ｔの固有値と、対角行列記憶部１３が記憶している対称３重対角行列Ｔとを用いて、対称３重対角行列Ｔの固有ベクトルを算出する。まず、コレスキー分解部３１は、対角行列記憶部１３から対称３重対角行列Ｔを読み出し、固有値記憶部１８から対称３重対角行列Ｔの固有値を読み出す。そして、コレスキー分解部３１は、図２３のフローチャートで示されるように各変換を行い、コレスキー分解の処理を行う。ここでは、ＬＶ型ツイスト分解法を用いる場合について説明する。【０２１９】
コレスキー分解部３１は、まず、対称３重対角行列Ｔの各要素の値から、ｑｋ（０），ｅｋ（０）を求める。そして、コレスキー分解部３１は、ミウラ変換を実行することにより、ｕ１（０），ｕ２（０）等を順次求めていく（ステップＳ１５０１）。【０２２０】
次に、コレスキー分解部３１は、ミウラ変換で得られたｕｋ（０）を用いて、ｓｔｄＬＶｖ変換を実行することにより、ｕ１（＋１）等を順次求めていく（ステップＳ１５０２）。また、コレスキー分解部３１は、ミウラ変換で得られたｕｋ（０）を用いて、ｒｄＬＶｖ変換を実行することにより、ｕ２ｍ−１（−１）等を順次求めていく（ステップＳ１５０３）。【０２２１】
最後に、コレスキー分解部３１は、ｓｔｄＬＶｖ変換で得られたｕｋ（＋１）及びｒｄＬＶｖ変換で得られたｕｋ（−１）を用いて、逆ミウラ変換を実行することにより、ｑ１（±１）、ｅ１（±１）等を順次求めていく（ステップＳ１５０４）。ｑｋ（±１）及びｅｋ（±１）が求まると、すなわち、上２重対角行列Ｂ（＋１）と、下２重対角行列Ｂ（−１）が求まると、コレスキー分解部３１は、それらをベクトル算出部３２に渡す。なお、コレスキー分解部３１は、各固有値について、それぞれコレスキー分解を行うものとする。【０２２２】
なお、ここでは、説明の便宜上、図２３のフローチャートを用いてコレスキー分解の処理を説明したが、図１６〜図２１のフローチャートで示されるように処理を行ってもよいことは言うまでもない。【０２２３】
ベクトル算出部３２は、固有値記憶部１８から固有値を読み出し、式５を用いてｋの値を決定する。ベクトル算出部３２は、そのｋの値を用いて、ｑｋ（±１）及びｅｋ（±１）からツイスト行列Ｎｋを求め、ＮｋＴｘｊ＝ｅｋを解くことによって、ｘｊを算出する（ステップＳ８０２）。【０２２４】
次に、ベクトル算出部３２は、算出したｘｊを正規化して、固有ベクトルｘｊを求めて固有ベクトル記憶部２０に蓄積する（ステップＳ８０３）。なお、ベクトル算出部３２は、各固有値について、ｘｊを算出する処理と、正規化する処理とを行うことによって、全ての固有ベクトルを算出することができる。このようにして、固有ベクトルの算出が終了する。
この後、行列記憶部１１が記憶している対称行列Ａの固有ベクトルを算出してもよい。その算出の方法は、前述の通りである。【０２２５】
また、固有値算出部１７が算出した固有値、固有ベクトル算出部１９が算出した固有ベクトルを出力する図示しない出力部を固有値分解装置１が備えてもよい。ここで、図示しない出力部による出力は、例えば、固有値等の表示デバイス（例えば、ＣＲＴや液晶ディスプレイなど）への表示でもよく、固有値等の所定の機器への通信回線を介した送信でもよく、固有値等のプリンタによる印刷でもよく、固有値等の記録媒体への蓄積でもよい。なお、その出力部は、出力を行うデバイス（例えば、表示デバイスやプリンタなど）を含んでもよく、あるいは含まなくてもよい。また、その出力部は、ハードウェアによって実現されてもよく、あるいは、それらのデバイスを駆動するドライバ等のソフトウェアによって実現されてもよい。【０２２６】
なお、ここでは、コレスキー分解部３１がＬＶ型ツイスト分解法によってコレスキー分解を行う場合について説明したが、コレスキー分解部３１は、ｑｄ型ツイスト分解法によってコレスキー分解を行ってもよい。例えば、コレスキー分解部３１は、算出された固有値の分布を見て、分布が密でない場合には、ｑｄ型ツイスト分解法を用いるようにしてもよい。【０２２７】
また、上記説明では、行列Ｂ（０）が上２重対角行列である場合について説明したが、行列Ｂ（０）は下２重対角行列であっても固有値分解を同様に実行することができうる。ただし、行列Ｔの各要素と、行列Ｂ（０）の各要素との対応を示す式が上述のものと異なることになる。【０２２８】
また、上記説明では、固有値算出部１７が全ての固有値を算出する場合について説明したが、一部の固有値のみを算出するようにしてもよい。例えば、図９において、固有値算出部１７は、行列Ｔ１，Ｔ２までは、全ての固有値を算出する必要があるが、行列Ｔ１，Ｔ２から対称３重対角行列Ｔの固有値を算出する際に、必要な範囲で固有値を算出してもよい。したがって、固有値算出部１７は、対称３重対角行列Ｔの少なくとも１個の固有値を算出するものであってもよい。その場合には、不必要な固有値を算出することに伴う余分な処理を実行しなくてよいため、処理負荷を軽減することができうる。このように、一部の固有値のみを算出する場合には、例えば、対称３重対角行列Ｔの固有値を算出する処理に費やす計算時間を約（１＋ｋ／ｍ）／２倍にすることができうる。ここで、ｍは行列サイズであり、ｋは求める固有値の個数である。【０２２９】
また、上記説明では、固有ベクトル算出部１９が算出された固有値に対応する全ての固有ベクトルを算出する場合について説明したが、上記説明から明らかなように、固有ベクトルを算出する処理は、固有値ごとに行うことができる。したがって、固有ベクトル算出部１９は、対称３重対角行列Ｔの少なくとも１個の固有ベクトルを算出するものであってもよい。このように、固有ベクトル算出部１９は、必要な範囲で固有ベクトルを算出することができ、不必要な固有ベクトルを算出することに伴う余分な処理を実行しなくてよいため、処理負荷を軽減することができうる。【０２３０】
次に、本実施の形態による固有値分解装置１における並列処理について説明する。
本実施の形態による固有値分解装置１では、固有値の計算及び固有ベクトルの計算において、並列的に処理を実行することができる。例えば、図２４〜図２６で示されるように、固有値分解装置１において、固有値分解部１５、固有値算出部１７、コレスキー分解部２２、ベクトル算出部２３，コレスキー分解部３１、ベクトル算出部３２において、それぞれ並列処理を行ってもよい。【０２３１】
固有値分解部１５は、複数の固有値分解手段１５ａ、１５ｂを備え、その複数の固有値分解手段１５ａ、１５ｂが、分割後の対称３重対角行列の固有値と、行列要素とを算出する処理を並列実行してもよい。例えば、図９において、固有値分解手段１５ａが行列Ｔ１１１，Ｔ１１２，Ｔ１２１，Ｔ１２１の固有値分解を実行し、固有値分解手段１５ｂが行列Ｔ２１１，Ｔ２１２，Ｔ２２１，Ｔ２２２の固有値分解を実行してもよい。【０２３２】
固有値算出部１７は、複数の固有値算出手段１７ａ、１７ｂを備え、その複数の固有値算出手段１７ａ、１７ｂが、分割元の対称３重対角行列の固有値と、行列要素とを算出する処理を並列実行してもよい。例えば、図９において、固有値算出手段１７ａが行列Ｔ１１，Ｔ１２，Ｔ１，Ｔの固有値等を算出する処理を実行し、固有値算出手段１７ｂが行列Ｔ２１，Ｔ２２，Ｔ２の固有値等を算出する処理を実行してもよい。【０２３３】
コレスキー分解部２２は、複数のコレスキー分解手段２２ａ、２２ｂを備え、その複数のコレスキー分解手段２２ａ、２２ｂが、対称３重対角行列Ｔに関してコレスキー分解する処理を並列実行してもよい。例えば、コレスキー分解手段２２ａが半分の固有値についてコレスキー分解する処理を実行し、コレスキー分解手段２２ｂが残りの半分の固有値についてコレスキー分解する処理を実行してもよい。【０２３４】
ベクトル算出部２３は、複数のベクトル算出手段２３ａ、２３ｂを備え、その複数のベクトル算出手段２３ａ、２３ｂが、固有ベクトルを算出する処理を並列実行してもよい。例えば、ベクトル算出手段２３ａがコレスキー分解手段２２ａによってコレスキー分解された固有値に関して固有ベクトルを算出する処理を実行し、ベクトル算出手段２３ｂがコレスキー分解手段２２ｂによってコレスキー分解された固有値に関して固有ベクトルを算出する処理を実行してもよい。【０２３５】
コレスキー分解部３１が複数のコレスキー分解手段３１ａ、３１ｂを備えて、対称３重対角行列Ｔに関してコレスキー分解する処理を並列実行してもよいことは、コレスキー分解部２２と同様である。【０２３６】
ベクトル算出部３２が複数のベクトル算出手段３２ａ、３２ｂを備えて、固有ベクトルを算出する処理を並列実行してもよいことは、ベクトル算出部２３と同様である。【０２３７】
なお、固有値算出部１７の複数の固有値算出手段１７ａ、１７ｂは、図６で示される分割された２個の行列Ｔ１，Ｔ２の固有値と行列要素とから、分割元の行列Ｔ０の固有値と行列要素とを算出する処理を、並列実行してもよい。以下、その処理について説明する。【０２３８】
まず、式１から、各固有値を求める処理については、並列実行可能なことがわかる。例えば、固有値算出手段１７ａが半分の固有値を算出し、固有値算出手段１７ｂが残りの半分の固有値を算出してもよい。なお、固有値算出手段１７ａ、１７ｂのそれぞれは、行列Ｄ１２と、ベクトルｚとの各成分の値を有するものとする。【０２３９】
また、式２から、各固有ベクトルｑｉを求める処理についても、並列実行可能なことがわかる。例えば、固有値算出手段１７ａは、算出した固有値に対応する半分の固有ベクトルｑｉを算出し、固有値算出手段１７ｂは、算出した固有値に対応する残りの半分の固有ベクトルｑｉを算出してもよい。【０２４０】
また、式３，式４から、ｆ、ｌを求める処理についても、並列実行可能なことがわかる。例えば、固有値算出手段１７ａは、算出した固有ベクトルｑｉを用いて、ｆ、ｌの一部の要素を算出し、固有値算出手段１７ｂは、算出した固有ベクトルｑｉを用いて、ｆ、ｌの残りの要素を算出してもよい。その後、固有値算出手段１７ａ、あるいは固有値算出手段１７ｂが、固有値算出手段１７ａ、１７ｂの算出したｆ、ｌの各要素を統合することによって、式３，式４で示されるｆ、ｌを算出することができる。【０２４１】
ここでは、ｚの成分に０が含まれず、Ｄ１２の対角成分に同じ値が含まれない場合における固有値算出手段１７ａ、１７ｂの処理について説明したが、ｚの成分に０が含まれる場合、Ｄ１２の対角成分に同じ値が含まれる場合であっても、固有方程式を解く処理、固有ベクトルｑｉを算出する処理、固有ベクトルから分割元の行列の行列要素を算出する処理を同様に並列実行することができうる。【０２４２】
上記のような各構成要素における並列処理において、各手段が固有値等の情報を格納するメモリを用いる場合に、複数の手段が同一のメモリ、すなわち共有メモリを使用してもよく、あるいは、各手段がそれぞれ別のメモリを使用してもよい。【０２４３】
なお、図２４〜図２６を用いて、固有値分解部１５や固有値算出部１７等が２個の手段によって並列処理を実行する場合について説明したが、固有値分解部１５や固有値算出部１７等が３以上の手段を備え、並列処理を実行してもよい。また、固有値分解部１５、固有値算出部１７、コレスキー分解部２２、ベクトル算出部２３，コレスキー分解部３１、ベクトル算出部３２のそれぞれにおいて並列処理が実行される場合について説明したが、いずれかの任意の１以上の部において、並列処理が実行されなくてもよい。例えば、固有値分解部１５において並列処理が行われなくてもよい。【０２４４】
また、その並列処理は、１の装置において、２以上のＣＰＵ等を用いて実行されてもよく、２以上の装置において実行されてもよい。例えば、図２７で示されるように、装置Ａと、装置Ｂとがそれぞれ固有値分解手段１５ａ、１５ｂを備え、各装置において、固有値分解の処理が並列実行されてもよい。この場合には、装置Ａの固有値分解手段１５ａを備える固有値分解部１５−１と、装置Ｂの固有値分解手段１５ｂを備える固有値分解部１５−２とによって固有値分解部１５が構成されることになる。したがって、固有値分解装置１は、装置Ａと、装置Ｂとからなるシステムを構成することになる。ここでは、固有値分解部１５の並列処理について説明したが、その他の固有値算出部１７やコレスキー分解部２２等についても、２以上の装置による並列処理を行ってもよい。【０２４５】
次に、数値実験による性能の評価について説明する。ここでは、本実施の形態による固有値分解装置１の方法（ＤＤＣ）と、標準的な分割統治法（ＤＣ）と、ＱＲ法（ＱＲ）と、二分法と逆反復法とを組み合わせた方法（Ｂ&Ｉ）とを比較した。なお、本実施の形態による固有値分解装置１の方法では、コレスキー分解をｑｄ型ツイスト分解法で行い、逆反復（図１２のステップＳ７０１参照）を１回用いた。分割統治法としては、ＬＡＰＡＣＫで提供されているＤＳＴＥＤＣを用いた。ＱＲ法としては、ＬＡＰＡＣＫで提供されているＤＳＴＥＱＲを用いた。二分法と逆反復法とを組み合わせた方法としては、ＬＡＰＡＣＫで提供されているＤＳＴＥＶＲのソースコードを一部変更したものを用いた。実験で用いるＬＡＰＡＣＫはＢＬＡＳ（Ｂａｓｉｃ  Ｌｉｎｅａｒ  Ａｌｇｅｂｒａ  Ｓｕｂｐｒｏｇｒａｍｓ）を呼び出す。また、本実施の形態による固有値分解装置１の固有値分解部１５では、ＱＲ法を用いて固有値分解を行った。また、この数値実験では、ランダムに与えた対角行列を、ランチョス法を用いて対称３重対角化した行列を用いた。数値実験結果のグラフにおいて、横軸は、数値実験で用いた行列の大きさである。【０２４６】
図２８は、計算の実行時間を示す図である。実行時間が短いほど、高速に計算できることになる。図２８の数値実験結果のグラフからわかるように、本実施の形態による固有値分解装置１の方法では、他の方法に比べて高速に固有値分解を行うことができることがわかる。【０２４７】
図２９は、算出された固有ベクトルの直交性を示す図である。値の小さい方が、直交性が高いことになる。図２９の数値実験結果のグラフからわかるように、本実施の形態による固有値分解装置１の方法では、他の方法に比べて、固有ベクトルの直交性は高くない。本実施の形態による固有値分解装置１の方法では、二分法と逆反復法とを組み合わせた方法と同程度の直交性を有することがわかる。これは、固有値をまず算出し、その後に、固有ベクトルを算出していることに起因していると考えられる。【０２４８】
図３０は、全固有値の相対誤差の総和を示す図である。値の小さい方が、誤差が小さく、精度が高いことになる。図３０の数値実験結果のグラフからわかるように、本実施の形態による固有値分解装置１の方法では、他の方法に比べて固有値の精度が高いことがわかる。分割統治法により算出した固有値と同程度の固有値の精度を有する。本実施の形態による固有値分解装置１の方法では、分割統治法によって固有値を算出しているため、当然の結果であると考えられる。【０２４９】
図３１は、前固有ベクトルの誤差の総和を示す図である。値の小さい方が、誤差が小さく、精度が高いことになる。図３１の数値実験結果のグラフからわかるように、本実施の形態による固有値分解装置１の方法では、他の方法に比べて固有ベクトルの精度が高いことがわかる。【０２５０】
このように、本実施の形態による固有値分解装置１の方法では、高速に固有値分解をすることができ、さらに、固有値、固有ベクトルの精度の高いことがわかる。また、前述のように、並列性にも優れているため、並列処理を行うことによってさらに高速化を図ることも可能である。【０２５１】
［主成分分析による顔画像認識］
次に、主成分分析による顔画像認識に固有値分解を応用する場合について説明する。【０２５２】
主成分分析による顔画像認識の処理は、以下の各ステップによって行われる。
（１）射影行列の作成
（２）個人の顔画像に関する情報の登録
（３）入力された顔画像と、登録されている個人の顔画像に関する情報とを比較することによる認識【０２５３】
まず、上記（１）の射影行列の作成について説明する。
種々の人物の顔画像データを用意する。その顔画像データは、例えば、スキャナやデジタルスチルカメラ、デジタルビデオカメラ等の光学的読み取り機器によって読み取られた２次元画像データであってもよい。各顔画像データは、各画素値を成分とするベクトルで表示することができる。そのベクトルを、次のようにする。ここでは、Ｍ個の顔画像データを扱うものとする。
【数７２】
2008018188.tif 000074
【０２５４】
次に、その顔画像データのベクトルの平均値μを算出する。
【数７３】
2008018188.tif 000075
【０２５５】
その算出したμを用いて、共分散行列Ｃを次のように算出する。
【数７４】
2008018188.tif 000076
【０２５６】
その後、上述の固有値分解装置１によって、その行列Ｃに対して、固有値分解を行う。すなわち、行列記憶部１１に行列Ｃを蓄積し、その行列Ｃに対して対角化や、行列の分割等の処理を行うことによって、固有値分解を行う。その固有値分解された結果である固有値λｊは、固有値記憶部１８に蓄積されることになる。また、固有値分解された結果である固有ベクトルｖｊは、固有ベクトル記憶部２０で記憶されている固有ベクトルに対して、対角化で用いた直交行列を作用させることによって、前述のようにして求めることができる。なお、添字ｊは、固有値の降順になるように設定されているものとする。すなわち、固有値λ１が最大の固有値となる。【０２５７】
その求めた固有ベクトルのうち、対応する固有値が上位からｄ個の固有ベクトルを用いて、射影行列ηを次のように算出する。
【数７５】
2008018188.tif 000077
【０２５８】
次に、上記（２）の個人の顔画像に関する情報の登録について説明する。
各個人について、ｍ個の顔画像データ
【数７６】
2008018188.tif 000078
を用意する。ここで、添字ｋは、個人を識別するための添字である。そのベクトルに射影行列ηをそれぞれ掛けることにより、次式のようにしてｍ個のｄ次元ベクトルξｉを算出し、その平均値μｋを算出する。そして、その平均値μｋを登録しておく。
【数７７】
2008018188.tif 000079
【０２５９】
次に、上記（３）の入力された顔画像と、登録されている個人の顔画像に関する情報とを比較することによる認識について説明する。
ある人物の顔画像データのベクトルｘが入力されたとする。すると、その顔画像データのベクトルｘに射影行列ηを掛けることにより、ｄ次元ベクトルξを算出する。
【数７８】
2008018188.tif 000080
【０２６０】
次に、そのｄ次元ベクトルξと、各個人のｄ次元ベクトルξｉの平均値μｋとの差のノルムを計算する。
【数７９】
2008018188.tif 000081
【０２６１】
そのノルムがあらかじめ設定されているしきい値より小さい場合に、入力された顔画像が、添字ｋで識別される個人の顔画像であると判断してもよい。あるいは、すべての添字ｋについて上記のノルムを計算し、入力された顔画像が、最小のノルムに対応する添字ｋで識別される個人の顔画像であると判断してもよい。このようにして、固有値分解装置１を用いて、主成分分析による顔画像認識を行うことができうる。【０２６２】
［２次元画像から３次元へ復元する画像処理への応用］
次に、物体の２次元画像から３次元へ復元する画像処理に固有値分解を応用する場合について説明する。【０２６３】
複数の２次元画像から３次元復元を行う処理は、以下の各ステップによって行われる。
（１）２次元画像から特徴点を抽出するステップ。
（２）特徴点データより形状（元の物体の特徴点の３次元座標データ）及び回転（３次元データから特徴点データへの変換）に関するデータを計算するステップ。
（３）形状及び回転のデータより可視化を行うステップ。【０２６４】
以下、上記（１）、（２）の各ステップについて、図３２のフローチャートを用いて説明する。ここで、２次元画像は、例えば、スキャナやデジタルスチルカメラ、デジタルビデオカメラ等の光学的読み取り機器によって読み取られた２次元画像であってもよい。【０２６５】
ステップＳ５００１において、２次元画像ｊ（ｊ＝１，・・・，ｍ、ｍは３以上の整数）から特徴点ｉ（ｉ＝１，・・・，ｎ、ｎは２以上の整数）の座標（ｘｉｊ，ｙｉｊ）を抽出する。取り扱う２次元画像は、弱中心射影画像であることが好ましい。このとき、次の式が成り立つ。
【数８０】
2008018188.tif 000082
【０２６６】
ここで、ｓｊは物体のスケールに相対するｊ番目の画像のスケール、ｒ１，ｊ，ｒ２，ｊはそれぞれ物体座標系に相対するｊ番目のカメラ座標系の回転行列の１番目と２番目の行ベクトル、（Ｘｉ，Ｙｉ，Ｚｉ）Ｔはｉ番目の点の３次元座標である。物体のスケールは１番目の画像のスケールと同じにし（ｓ１＝１）、物体の座標系の姿勢は１番目の画像のカメラ座標系と同じにする（ｒ１，１＝（１，０，０）Ｔ，ｒ２，１＝（０，１，０）Ｔ）。ｍ枚全ての画像におけるｎ個全ての点の座標を行列Ｄの要素として並べると、次式のようになる。【０２６７】
【数８１】
2008018188.tif 000083
【０２６８】
ＭとＳの形から分かるように、Ｄのランクは３である。ここで、ステップＳ５００１において、行列Ｄが与えられている。以下、回転に関するデータＭ及び形状Ｓを求める。【０２６９】
そこで、行列Ｄの特異値分解
Ｄ＝ＵΣＶＴ
を考える。ここで、Σは特異値を大小順に対角線上に並べたもので、Ｕ及びＶはそれぞれ左直交行列、及び右直交行列である。この特異値分解において、前述の固有値分解を用いることができる。すなわち、固有値分解装置１において、行列記憶部１１が記憶している対称行列Ａを、ＤＴＤとすることにより、前述の説明のようにして、対称行列Ａの固有値分解がなされることになる。なお、固有値分解装置１において、固有ベクトル記憶部２０に蓄積されるのは、対称行列ＤＴＤを変換した対称３重対角行列Ｔの固有ベクトルであるので、その固有ベクトルを、前述の説明のようにして、行列ＤＴＤの固有ベクトルに変換する必要がある。また、その行列ＤＴＤの各固有値の１／２乗（平方根）が各特異値となる。すなわち、特異値σｊ＝（λｊ）１／２となる。また、固有ベクトルと、右直交行列Ｖの各列を構成する右特異ベクトルｘｊとは等しいため、固有ベクトルを各列に有する行列が、右直交行列Ｖとなる。なお、特異値σｊが０でない場合には、左直交行列Ｕの各列を構成する左特異ベクトルｙｊを次のようにして求めることができる。【０２７０】
【数８２】
2008018188.tif 000084
となる。一方、特異値σｊが０である場合には、
【数８３】
2008018188.tif 000085
を解くことにより、ｙｊを算出することができる。ここで、次の置き換えを行うことによって、左特異ベクトルｙｊを正規化する。
【数８４】
2008018188.tif 000086
【０２７１】
このようにして、左特異ベクトルｙｊを求めることができ、この左特異ベクトルｙｊを用いて、Ｕ＝（ｙ１，ｙ２，…，ｙｍ）とすることにより、左直交行列Ｕを算出することができる。【０２７２】
ここで、画像のデジタル誤差のため、ゼロでない特異値は３つ以上出てくる。しかし、４番目以降の特異値はノイズによるもので、最初の３つの特異値と比べて格段に小さい。そこで、ステップＳ５００２では、最初の３つの特異値に対して特異ベクトルを計算する。すなわち、本実施の形態による固有値分解装置１では、最初の３個の固有値に対して固有ベクトルを計算する。採用する３個のベクトルをまとめると、次式となる。【０２７３】
Ｄ'＝Ｌ'Σ'Ｒ'Ｔ＝Ｍ'Ｓ'
ただし、Ｍ'＝Ｌ'（Σ'）１／２、Ｓ'＝（Σ'）１／２Ｒ'Ｔ、Ｄ'は‖Ｄ−Ｄ'‖を最小にするランク３の行列である。【０２７４】
次に、Ｄ'からＭ及びＳを求めたいが、その組合せは唯一ではない。なぜなら、任意の正則行列Ｃが
Ｄ'＝（Ｍ'Ｃ）（Ｃ−１Ｓ'）
を満たすからである。そこで、上式におけるＣをＭ＝Ｍ'Ｃを満たすように決める。Ｃは下記の式を満たす。
【数８５】
2008018188.tif 000087
【０２７５】
Ｅ＝ＣＣＴとすると、上式からＥの６つの要素に関する２ｍ＋１個の線形方程式が得られる。ｍ≧３であるので、Ｅの要素を一意に決めることができる。ステップＳ５００３において、行列Ｅを求める。【０２７６】
次に、ステップＳ５００４において、ＥからＣを求める。Ｃの自由度（９）はＥの自由度（６）より多い。そこで、条件ｒ１ｊ＝（１，０，０）Ｔ，ｒ２ｊ＝（０，１，０）Ｔを加えれば、Ｃを決めることができる。このとき２つの解（Ｎｅｃｋｅｒ  Ｒｅｖｅｒｓａｌ）が出る。
次に、ステップＳ５００５において、Ｍ＝Ｍ'Ｃ及びＳ＝Ｃ−１Ｓ'より、回転に関するデータＭ及び形状Ｓが決まる。【０２７７】
［文書検索への応用］
次に、文書検索処理に固有値分解を応用する場合について説明する。文書中からその文書の内容に関連する索引語を抽出し、索引語の重みを計算する処理の後、ベクトル空間モデルでは、この索引語の重みを要素とするベクトルで文書を表現する。ここで、検索対象となる文書をｄ１，ｄ２，・・・，ｄｎとし、これら文書集合全体を通して全部でｍ個の索引語ｗ１，ｗ２，・・・，ｗｍがあるとする。このとき、文書ｄｊは、次のようなベクトルで表現されることになる。これを文書ベクトルと呼ぶ。
【数８６】
2008018188.tif 000088
【０２７８】
ここで、ｄｉｊは索引語ｗｉの文書ｄｊにおける重みである。また、文書集合全体は、次のようなｍ×ｎ行列Ｄによって表現することができる。
【数８７】
2008018188.tif 000089
【０２７９】
行列Ｄを索引語文書行列と呼ぶ。索引語文書行列の各列は文書に関する情報を表している文書ベクトルであるが、同様に、索引語文書行列の各行は索引語に関する情報を表しているベクトルであり、これを索引語ベクトルと呼ぶ。検索質問も、文書と同様に、索引語の重みを要素とするベクトルで表現することができる。検索質問文に含まれる索引語ｗｉの重みをｑｉとすると、検索質問ベクトルｑは次のように表されることになる。
【数８８】
2008018188.tif 000090
【０２８０】
実際の文書検索においては、与えられた検索質問文と類似した文書を見つけ出す必要があるが、検索質問ベクトルｑと各文書ベクトルｄｊの間の類似度を計算することにより行う。ベクトル間の類似度の定義としてはさまざまなものが考えられるが、文書検索においてよく用いられているものはコサイン尺度（２つのベクトルのなす角度）または内積である。【０２８１】
【数８９】
2008018188.tif 000091
なお、ベクトルの長さが１に正規化（コサイン正規化）されている場合には、コサイン尺度と内積とは一致する。【０２８２】
図３３は、本実施の形態による固有値分解装置１を利用した文書検索方法の一例を示すフローチャートである。
ステップＳ６００１において、質問ベクトルｑを受け取る。【０２８３】
ここでは、Ｄの近似行列Ｄｋを使った検索を考える。ベクトル空間モデルでは、検索質問ベクトルｑと索引語文書行列Ｄ中の各文書ベクトルｄｊの間の類似度を計算することにより検索を行うが、ここではＤの代わりにＤｋを使う。ベクトル空間モデルでは、文書ベクトルの次元数は索引語の総数と等しい。したがって、検索対象となる文書の数が増えるに従い、文書ベクトルの次元数も増加する傾向にある。しかし、次元数が増加してくると、コンピュータのメモリによる制限や検索時間の増大などの問題が生じてくるばかりでなく、文書中に含まれる不必要な索引語がノイズ的な影響を及ぼし、検索精度を低下させてしまうという現象も起こってくる。潜在的意味インデキシング（ｌａｔｅｎｔ  ｓｅｍａｎｔｉｃ  ｉｎｄｅｘｉｎｇ；ＬＳＩ）は、高次元の空間にある文書ベクトルを低次元の空間へと射影することにより、検索精度の改善を図る技術である。高次元の空間では別々に扱われていた索引語が、低次元の空間では相互に関連を持ったものとして扱われる可能性もあるため、索引語の持つ意味や概念に基づく検索を行うことができる。たとえば、通常のベクトル空間モデルでは"ｃａｒ"という索引語と"ａｕｔｏｍｏｂｉｌｅ"という索引語はまったく別物であり、一方の索引語による質問ではもう片方の索引語を含んだ文書を検索することができない。しかし、低次元の空間ではこれらの意味的に関連した索引語は１つの次元に縮退することが期待できるため、"ｃａｒ"という検索質問によって"ｃａｒ"を含む文書ばかりでなく"ａｕｔｏｍｏｂｉｌｅ"を含む文書をも検索することが可能となる。潜在的意味インデキシングでは、特異値分解により高次元ベクトルの次元削減を行うが、これは基本的に多変量解析における主成分分析と等価である。【０２８４】
ステップＳ６００２において、ｋを選択する。ここでは、ｋ<ｒなるｋの値を選択する。ｒ＝ｍｉｎ（ｎ，ｍ）である。ｋの値は、予め与えられていてもよいし、計算ごとに選択可能であってもよい。【０２８５】
次に、ステップＳ６００３では、行列Ｄの特異値分解を行う。この特異値分解として、前述の２次元画像から３次元へ復元する画像処理への応用で説明したように、固有値分解による方法を用いることができる。すなわち、固有値分解装置１において、行列記憶部１１が記憶している対称行列Ａを、このたびの行列ＤＴＤとすることにより、前述の説明のようにして、対称行列Ａの固有値分解がなされ、その結果として、行列Ｄの特異値分解を行うことができうる。なお、この特異値分解では、計算された特異値のうち、大きい順に１番目からｋ番目までのｋ個の特異値に対してＤの特異ベクトルを算出する。すなわち、固有値分解において、計算された固有値のうち、大きい順に１番目からｋ番目までのｋ個の固有値に対して固有ベクトルを算出し、この固有ベクトルを用いて特異ベクトルを算出すればよい。ｋは、ステップＳ６００２で選択された値である。【０２８６】
すなわち、
Ｄｋ＝ＵｋΣＶｋＴ
なるＵｋ及びＶｋを計算する。ここで、Ｕｋは、最初のｋ個の左特異ベクトルのみから構成されるｍ×ｋ行列であり、Ｖｋは、最初のｋ個の右特異ベクトルのみから構成されるｎ×ｋ行列であり、Σは、最初のｋ個の特異値のみから構成されるｋ×ｋ対角行列である。【０２８７】
次に、ステップ６００４において、行列Ｄｋと質問ベクトルｑとの類似度を計算する。いま、ベクトルｅｊをｎ次元の単位ベクトルとすると、Ｄｋのｊ番目の文書ベクトルはＤｋｅｊで表すことができる。文書ベクトルＤｋｅｊと検索質問ベクトルｑとの間の類似度計算は、
【数９０】
2008018188.tif 000092
としてもよいが、別の定義を用いてもよい。上式では、ＤｋをＵｋ，Σｋ，Ｖｋから再構成する必要はなく特異値分解の結果から、直接、類似度を計算できることを示している。上式の中に現われるΣｋＶｋＴｅｊは、
ΣｋＶｋＴｅｊ＝ＵｋＴＤｋｅｊ
と書き直すことができる。この式の右辺は、近似行列Ｄｋにおけるｊ番目の文書ベクトルの基底Ｕｋのもとでの座標（文書のｋ次元表現）を表している。同様に、上式の中のＵｋＴｑは、検索質問ベクトルｑの基底Ｕｋのもとでの座標（検索質問のｋ次元表現）である。【０２８８】
ステップＳ６００５において、ステップＳ６００４において計算された類似度を基準に、検索結果を出力する。
このように、本実施の形態による固有値分解装置１は、画像処理や検索処理等の種々の処理に対して用いることができうる。なお、上記以外の処理に用いてもよいことは言うまでもない。【０２８９】
以上のように、本実施の形態による固有値分解装置１では、分割統治法を用いて固有値のみを算出するため、標準的な分割統治法で実行されるベクトル更新が必要なくなり、標準的な分割統治法よりも非常に高速である。また、固有値から固有ベクトルを算出する処理も、高速に処理することができる。さらに、ＱＲ法では固有値を算出する段階の並列化が困難であるが、本実施の形態による固有値分解装置１では、固有値を算出する段階、及び固有値から固有ベクトルを算出する段階のそれぞれにおいて、本質的に高い並列性を持つ。また、本実施の形態による固有値分解装置１では、他の方法と比べて高い精度の得られることがわかる。【０２９０】
なお、本実施の形態による固有値分解装置１は、スタンドアロンの装置であってもよく、サーバクライアントシステムにおけるサーバ装置であってもよく、前述のように、複数の装置から構成されるシステムであってもよい。固有値分解装置１がサーバクライアントシステムにおけるサーバ装置である場合には、行列記憶部１１で記憶される対称行列Ａや、固有値記憶部１８で記憶される固有値や固有ベクトル記憶部２０で記憶される固有ベクトル等は、インターネットやイントラネット等の通信回線を介して送受信されてもよい。【０２９１】
また、本実施の形態による固有値分解装置１における固有ベクトルの算出において、ある固有値を用いた固有ベクトルの算出では、ｑｄ型ツイスト分解法を用い、他の固有値を用いた固有ベクトルの算出では、ＬＶ型ツイスト分解法を用いるようにしてもよい。例えば、値の近接している固有値については、ＬＶ型ツイスト分解法を用いて固有ベクトルの算出を行い、値が近接していない固有値については、ｑｄ型ツイスト分解法を用いて固有ベクトルの算出を行ってもよい。固有値の値が近接しているかどうかは、所定の記録媒体等においてあらかじめ設定されているしきい値と比較することによって判断してもよい。【０２９２】
また、上記実施の形態において、各処理または各機能は、単一の装置または単一のシステムによって集中処理されることによって実現されてもよく、あるいは、複数の装置または複数のシステムによって分散処理されることによって実現されてもよい。【０２９３】
また、上記実施の形態において、各構成要素は専用のハードウェアにより構成されてもよく、あるいは、ソフトウェアにより実現可能な構成要素については、プログラムを実行することによって実現されてもよい。例えば、ハードディスクや半導体メモリ等の記録媒体に記録されたソフトウェア・プログラムをＣＰＵ等のプログラム実行部が読み出して実行することによって、各構成要素が実現され得る。なお、上記実施の形態における固有値分解装置を実現するソフトウェアは、以下のようなプログラムである。つまり、このプログラムは、コンピュータに、対称３重対角行列Ｔが記憶される対角行列記憶部から前記対称３重対角行列Ｔを読み出し、当該対称３重対角行列Ｔを２個の対称３重対角行列に分割して前記対角行列記憶部に蓄積し、その対称３重対角行列を２個の対称３重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各対称３重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割ステップと、前記あらかじめ決められた大きさ以下の各対称３重対角行列を前記対角行列記憶部から読み出し、前記各対称３重対角行列に対して固有値分解を行い、前記各対称３重対角行列の固有値と、前記各対称３重対角行列の固有ベクトルからなる直交行列の行列要素とを少なくとも算出し、当該固有値と当該行列要素とを固有値分解記憶部に蓄積する固有値分解ステップと、各対称３重対角行列の固有値と、行列要素とを前記固有値分解記憶部から読み出し、前記固有値と、前記行列要素とから、分割元の対称３重対角行列の固有値と、分割元の対称３重対角行列の行列要素とを算出して前記固有値分解記憶部に蓄積し、その分割元の対称３重対角行列の固有値と、行列要素とを算出する処理を、対称３重対角行列Ｔの少なくとも１個の固有値を算出するまで繰り返し、前記対称３重対角行列Ｔの少なくとも１個の固有値を固有値記憶部に蓄積する固有値算出ステップと、前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、前記固有値記憶部から前記対称３重対角行列Ｔの固有値を読み出し、前記対称３重対角行列Ｔとその固有値とから、ツイスト分解法を用いて前記対称３重対角行列Ｔの少なくとも１個の固有ベクトルを算出して固有ベクトル記憶部に蓄積する固有ベクトル算出ステップと、を実行させるためのものである。
また、このプログラムでは、前記固有ベクトル算出ステップにおいて、ｑｄ型ツイスト分解法により固有ベクトルを算出してもよい。【０２９４】
また、このプログラムでは、前記固有ベクトル算出ステップが、前記固有値記憶部から前記対称３重対角行列Ｔの固有値を読み出し、当該固有値のいずれかが負の値である場合に、前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、当該対称３重対角行列Ｔの固有ベクトルを変化させることなく、当該対称３重対角行列Ｔのすべての固有値が正の値となるように、当該対称３重対角行列Ｔを正定値化し、正定値化した後の固有値を前記固有値記憶部に蓄積し、正定値化した後の対称３重対角行列Ｔを前記対角行列記憶部に蓄積する正定値化ステップと、前記対角行列記憶部からすべての固有値が正の値である対称３重対角行列Ｔを読み出し、前記固有値記憶部から前記対称３重対角行列Ｔの正の値の固有値を読み出し、前記対称３重対角行列Ｔの各要素に関してｑｄ型変換を行うことによって、前記対称３重対角行列Ｔを上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）にコレスキー分解するコレスキー分解ステップと、前記上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）の各要素と、前記対称３重対角行列Ｔの正の値の固有値とを用いて固有ベクトルを算出して前記固有ベクトル記憶部に蓄積するベクトル算出ステップと、を備えてもよい。
また、このプログラムでは、前記固有ベクトル算出ステップにおいて、ＬＶ型ツイスト分解法により固有ベクトルを算出してもよい。【０２９５】
また、このプログラムでは、前記固有ベクトル算出ステップが、前記対角行列記憶部から前記対称３重対角行列Ｔを読み出し、前記固有値記憶部から前記対称３重対角行列Ｔの固有値を読み出し、前記対称３重対角行列Ｔの各要素に関してミウラ変換、ｄＬＶｖ型変換、逆ミウラ変換を行うことによって、前記対称３重対角行列Ｔを上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）にコレスキー分解するコレスキー分解ステップと、前記上２重対角行列Ｂ（＋１）及び下２重対角行列Ｂ（−１）の各要素と、前記対称３重対角行列Ｔの固有値とを用いて固有ベクトルを算出して前記固有ベクトル記憶部に蓄積するベクトル算出ステップと、を備えてもよい。【０２９６】
また、このプログラムでは、コンピュータに、対称行列Ａが記憶される行列記憶部から前記対称行列Ａを読み出し、前記対称行列Ａを３重対角化した前記対称３重対角行列Ｔを算出して前記対角行列記憶部に蓄積する対角化ステップをさらに実行させてもよい。
なお、上記プログラムにおいて、情報を蓄積するステップなどでは、ハードウェアでしか行われない処理は少なくとも含まれない。【０２９７】
また、このプログラムは、サーバなどからダウンロードされることによって実行されてもよく、所定の記録媒体（例えば、ＣＤ−ＲＯＭなどの光ディスクや磁気ディスク、半導体メモリなど）に記録されたプログラムが読み出されることによって実行されてもよい。【０２９８】
また、このプログラムを実行するコンピュータは、単数であってもよく、複数であってもよい。すなわち、集中処理を行ってもよく、あるいは分散処理を行ってもよい。【０２９９】
図３４は、上記プログラムを実行して、上記実施の形態による固有値分解装置１を実現するコンピュータの外観の一例を示す模式図である。上記実施の形態は、コンピュータハードウェア及びその上で実行されるコンピュータプログラムによって実現されうる。【０３００】
図３４において、コンピュータシステム１００は、ＣＤ−ＲＯＭ（Ｃｏｍｐａｃｔ  Ｄｉｓｋ  Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）ドライブ１０５、ＦＤ（Ｆｌｅｘｉｂｌｅ  Ｄｉｓｋ）ドライブ１０６を含むコンピュータ１０１と、キーボード１０２と、マウス１０３と、モニタ１０４とを備える。【０３０１】
図３５は、コンピュータシステムを示す図である。図３５において、コンピュータ１０１は、ＣＤ−ＲＯＭドライブ１０５、ＦＤドライブ１０６に加えて、ＣＰＵ（Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）１１１と、ブートアッププログラム等のプログラムを記憶するためのＲＯＭ（Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）１１２と、ＣＰＵ１１１に接続され、アプリケーションプログラムの命令を一時的に記憶すると共に、一時記憶空間を提供するＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）１１３と、アプリケーションプログラム、システムプログラム、及びデータを記憶するハードディスク１１４と、ＣＰＵ１１１、ＲＯＭ１１２等を相互に接続するバス１１５とを備える。なお、コンピュータ１０１は、ＬＡＮへの接続を提供する図示しないネットワークカードを含んでもよい。【０３０２】
コンピュータシステム１００に、上記実施の形態による固有値分解装置１の機能を実行させるプログラムは、ＣＤ−ＲＯＭ１２１、またはＦＤ１２２に記憶されて、ＣＤ−ＲＯＭドライブ１０５、またはＦＤドライブ１０６に挿入され、ハードディスク１１４に転送されてもよい。これに代えて、そのプログラムは、図示しないネットワークを介してコンピュータ１０１に送信され、ハードディスク１１４に記憶されてもよい。プログラムは実行の際にＲＡＭ１１３にロードされる。なお、プログラムは、ＣＤ−ＲＯＭ１２１やＦＤ１２２、またはネットワークから直接、ロードされてもよい。【０３０３】
また、行列記憶部１１、対角行列記憶部１３、固有値分解記憶部１６、固有値記憶部１８，固有ベクトル記憶部２０は、ＲＡＭ１１３やハードディスク１１４によって実現されてもよい。【０３０４】
プログラムは、コンピュータ１０１に、上記実施の形態による固有値分解装置１の機能を実行させるオペレーティングシステム（ＯＳ）、またはサードパーティプログラム等を必ずしも含まなくてもよい。プログラムは、制御された態様で適切な機能（モジュール）を呼び出し、所望の結果が得られるようにする命令の部分のみを含んでいてもよい。コンピュータシステム１００がどのように動作するのかについては周知であり、詳細な説明を省略する。
また、本発明は、以上の実施の形態に限定されることなく、種々の変更が可能であり、それらも本発明の範囲内に包含されるものであることは言うまでもない。
また、本発明のほんのいくつかの典型的な実施例について上で詳細に説明したが、その典型的な実施例において、発明の利益と新規な技術から実質的にはずれることなく多くの変更が可能であることを当業者は容易に認識することができるであろう。したがって、そのようなすべての変更は、本発明の範囲に含まれるものである。
【産業上の利用可能性】
【０３０５】
以上のように、本発明による固有値分解装置等によれば、固有値分解を高速に処理することができ、分子軌道法による化学計算処理や統計計算処理、情報検索処理、その他の固有値分解を用いる処理を実行する装置等において有用である。
【図面の簡単な説明】
【０３０６】
【図１】本発明の実施の形態１による固有値分解装置の構成を示すブロック図
【図２】同実施の形態による固有値分解装置の動作を示すフローチャート
【図３】同実施の形態による固有値分解装置の動作を示すフローチャート
【図４】同実施の形態による固有値分解装置の動作を示すフローチャート
【図５】同実施の形態における対称３重対角行列Ｔの分割について説明するための図
【図６】同実施の形態における対称３重対角行列Ｔの固有値の算出について説明するための図
【図７】同実施の形態による固有値分解装置の動作を示すフローチャート
【図８】同実施の形態による固有値分解装置の動作を示すフローチャート
【図９】同実施の形態における対称３重対角行列Ｔの固有値の算出について説明するための図
【図１０】同実施の形態における固有ベクトル算出部の構成を示すブロック図
【図１１】同実施の形態による固有値分解装置の動作を示すフローチャート
【図１２】同実施の形態による固有値分解装置の動作を示すフローチャート
【図１３】同実施の形態における固有ベクトル算出部の構成を示すブロック図
【図１４】同実施の形態による固有値分解装置の動作を示すフローチャート
【図１５】同実施の形態におけるコレスキー分解について説明するための図
【図１６】同実施の形態による固有値分解装置の動作を示すフローチャート
【図１７】同実施の形態による固有値分解装置の動作を示すフローチャート
【図１８】同実施の形態による固有値分解装置の動作を示すフローチャート
【図１９】同実施の形態による固有値分解装置の動作を示すフローチャート
【図２０】同実施の形態による固有値分解装置の動作を示すフローチャート
【図２１】同実施の形態による固有値分解装置の動作を示すフローチャート
【図２２】同実施の形態におけるｑｄ型ツイスト分解法とＬＶ型ツイスト分解法とについて説明するための図
【図２３】同実施の形態による固有値分解装置の動作を示すフローチャート
【図２４】同実施の形態による固有値分解装置の構成の他の一例を示すブロック図
【図２５】同実施の形態における固有ベクトル算出部の構成の他の一例を示すブロック図
【図２６】同実施の形態における固有ベクトル算出部の構成の他の一例を示すブロック図
【図２７】同実施の形態による固有値分解装置の構成の他の一例を示すブロック図
【図２８】同実施の形態による固有値分解装置の数値実験の結果を示す図
【図２９】同実施の形態による固有値分解装置の数値実験の結果を示す図
【図３０】同実施の形態による固有値分解装置の数値実験の結果を示す図
【図３１】同実施の形態による固有値分解装置の数値実験の結果を示す図
【図３２】同実施の形態における画像処理の一例を示すフローチャート
【図３３】同実施の形態における文書検索処理の一例を示すフローチャート
【図３４】同実施の形態におけるコンピュータシステムの外観一例を示す模式図
【図３５】同実施の形態におけるコンピュータシステムの構成の一例を示す図
【図１】
2008018188.tif 000093
【図２】
2008018188.tif 000094
【図３】
2008018188.tif 000095
【図４】
2008018188.tif 000096
【図５】
2008018188.tif 000097
【図６】
2008018188.tif 000098
【図７】
2008018188.tif 000099
【図８】
2008018188.tif 000100
【図９】
2008018188.tif 000101
【図１０】
2008018188.tif 000102
【図１１】
2008018188.tif 000103
【図１２】
2008018188.tif 000104
【図１３】
2008018188.tif 000105
【図１４】
2008018188.tif 000106
【図１５】
2008018188.tif 000107
【図１６】
2008018188.tif 000108
【図１７】
2008018188.tif 000109
【図１８】
2008018188.tif 000110
【図１９】
2008018188.tif 000111
【図２０】
2008018188.tif 000112
【図２１】
2008018188.tif 000113
【図２２】
2008018188.tif 000114
【図２３】
2008018188.tif 000115
【図２４】
2008018188.tif 000116
【図２５】
2008018188.tif 000117
【図２６】
2008018188.tif 000118
【図２７】
2008018188.tif 000119
【図２８】
2008018188.tif 000120
【図２９】
2008018188.tif 000121
【図３０】
2008018188.tif 000122
【図３１】
2008018188.tif 000123
【図３２】
2008018188.tif 000124
【図３３】
2008018188.tif 000125
【図３４】
2008018188.tif 000126
【図３５】
2008018188.tif 000127
【国際調査報告】
2008018188.tif 000128
2008018188.tif 000129
2008018188.tif 000130
2008018188.tif 000131