TOP > 国内特許検索 > 固有値分解装置、及び固有値分解方法 > 明細書

明細書 :固有値分解装置、及び固有値分解方法

発行国 日本国特許庁(JP)
公報種別 特許公報(B2)
特許番号 特許第5017666号 (P5017666)
登録日 平成24年6月22日(2012.6.22)
発行日 平成24年9月5日(2012.9.5)
発明の名称または考案の名称 固有値分解装置、及び固有値分解方法
国際特許分類 G06F  17/16        (2006.01)
G06T   1/00        (2006.01)
FI G06F 17/16 K
G06T 1/00 A
請求項の数または発明の数 16
全頁数 63
出願番号 特願2008-528725 (P2008-528725)
出願日 平成19年1月31日(2007.1.31)
新規性喪失の例外の表示 特許法第30条第1項適用 平成18年6月26日発行の情報処理学会研究報告等で発表
国際出願番号 PCT/JP2007/051575
国際公開番号 WO2008/018188
国際公開日 平成20年2月14日(2008.2.14)
優先権出願番号 2006215660
優先日 平成18年8月8日(2006.8.8)
優先権主張国 日本国(JP)
審査請求日 平成21年10月30日(2009.10.30)
特許権者または実用新案権者 【識別番号】504132272
【氏名又は名称】国立大学法人京都大学
発明者または考案者 【氏名】中村 佳正
【氏名】坪井 洋明
【氏名】譽田 太朗
【氏名】岩▲崎▼ 雅史
【氏名】高田 雅美
個別代理人の代理人 【識別番号】100115749、【弁理士】、【氏名又は名称】谷川 英和
審査官 【審査官】田中 幸雄
参考文献・文献 国際公開第2005/119507(WO,A1)
桑島豊ほか,実対称三重対角固有値問題の分割統治法の拡張,日本応用数理学会論文誌,日本,日本応用数理学会,2005年 6月25日,vol. 15, no. 2,89-115頁
調査した分野 G06F 17/16
G06T 1/00
特許請求の範囲 【請求項1】
対称3重対角行列Tが記憶される対角行列記憶部と、
前記対角行列記憶部から前記対称3重対角行列Tを読み出し、当該対称3重対角行列Tを2個の対称3重対角行列に分割して前記対角行列記憶部に蓄積し、その対称3重対角行列を2個の対称3重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各対称3重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割部と、
前記あらかじめ決められた大きさ以下の各対称3重対角行列の固有値と、当該各対称3重対角行列の固有ベクトルからなる直交行列の一部の要素である行列要素とが記憶される固有値分解記憶部と、
前記あらかじめ決められた大きさ以下の各対称3重対角行列を前記対角行列記憶部から読み出し、前記各対称3重対角行列に対して固有値分解を行い、前記各対称3重対角行列の固有値と、前記各対称3重対角行列の固有ベクトルからなる直交行列の行列要素とを少なくとも算出し、当該固有値と当該行列要素とを前記固有値分解記憶部に蓄積する固有値分解部と、
前記対称3重対角行列Tの固有値が記憶される固有値記憶部と、
各対称3重対角行列の固有値と、行列要素とを前記固有値分解記憶部から読み出し、前記固有値と、前記行列要素とから、分割元の対称3重対角行列の固有値と、分割元の対称3重対角行列の行列要素とを算出して前記固有値分解記憶部に蓄積し、その分割元の対称3重対角行列の固有値と、行列要素とを算出する処理を、対称3重対角行列Tの少なくとも1個の固有値を算出するまで繰り返し、前記対称3重対角行列Tの少なくとも1個の固有値を前記固有値記憶部に蓄積する固有値算出部と、
前記対称3重対角行列Tの固有ベクトルが記憶される固有ベクトル記憶部と、
前記対角行列記憶部から前記対称3重対角行列Tを読み出し、前記固有値記憶部から前記対称3重対角行列Tの固有値を読み出し、前記対称3重対角行列Tとその固有値とから、ツイスト分解法を用いて前記対称3重対角行列Tの少なくとも1個の固有ベクトルを算出して前記固有ベクトル記憶部に蓄積する固有ベクトル算出部と、を備えた固有値分解装置。
【請求項2】
前記固有ベクトル算出部は、qd型ツイスト分解法により固有ベクトルを算出する、請求項1記載の固有値分解装置。
【請求項3】
前記固有ベクトル算出部は、
前記固有値記憶部から前記対称3重対角行列Tの固有値を読み出し、当該固有値のいずれかが負の値である場合に、前記対角行列記憶部から前記対称3重対角行列Tを読み出し、当該対称3重対角行列Tの固有ベクトルを変化させることなく、当該対称3重対角行列Tのすべての固有値が正の値となるように、当該対称3重対角行列Tを正定値化し、正定値化した後の固有値を前記固有値記憶部に蓄積し、正定値化した後の対称3重対角行列Tを前記対角行列記憶部に蓄積する正定値化部と、
前記対角行列記憶部からすべての固有値が正の値である対称3重対角行列Tを読み出し、前記固有値記憶部から前記対称3重対角行列Tの正の値の固有値を読み出し、前記対称3重対角行列Tの各要素に関してqd型変換を行うことによって、前記対称3重対角行列Tを上2重対角行列B(+1)及び下2重対角行列B(-1)にコレスキー分解するコレスキー分解部と、
前記上2重対角行列B(+1)及び下2重対角行列B(-1)の各要素と、前記対称3重対角行列Tの正の値の固有値とを用いて固有ベクトルを算出して前記固有ベクトル記憶部に蓄積するベクトル算出部と、を備えた、請求項2記載の固有値分解装置。
【請求項4】
前記固有ベクトル算出部は、LV型ツイスト分解法により固有ベクトルを算出する、請求項1記載の固有値分解装置。
【請求項5】
前記固有ベクトル算出部は、
前記対角行列記憶部から前記対称3重対角行列Tを読み出し、前記固有値記憶部から前記対称3重対角行列Tの固有値を読み出し、前記対称3重対角行列Tの各要素に関してミウラ変換、dLVv型変換、逆ミウラ変換を行うことによって、前記対称3重対角行列Tを上2重対角行列B(+1)及び下2重対角行列B(-1)にコレスキー分解するコレスキー分解部と、
前記上2重対角行列B(+1)及び下2重対角行列B(-1)の各要素と、前記対称3重対角行列Tの固有値とを用いて固有ベクトルを算出して前記固有ベクトル記憶部に蓄積するベクトル算出部と、を備えた、請求項4記載の固有値分解装置。
【請求項6】
前記コレスキー分解部は、複数のコレスキー分解手段を備え、
前記複数のコレスキー分解手段が、前記対称3重対角行列Tをコレスキー分解する処理を並列実行する、請求項3または請求項5記載の固有値分解装置。
【請求項7】
前記ベクトル算出部は、複数のベクトル算出手段を備え、
前記複数のベクトル算出手段が、固有ベクトルを算出する処理を並列実行する、請求項3、請求項5または請求項6記載の固有値分解装置。
【請求項8】
前記固有値算出部は、複数の固有値算出手段を備え、
前記複数の固有値算出手段が、分割元の対称3重対角行列の固有値と、行列要素とを算出する処理を並列実行する、請求項1から請求項7のいずれか記載の固有値分解装置。
【請求項9】
前記固有値分解部は、複数の固有値分解手段を備え、
前記複数の固有値分解手段が、対称3重対角行列に対して固有値分解を行う処理を並列実行する、請求項1から請求項8のいずれか記載の固有値分解装置。
【請求項10】
対称行列Aが記憶される行列記憶部と、
前記対称行列Aを前記行列記憶部から読み出し、前記対称行列Aを3重対角化した前記対称3重対角行列Tを算出して前記対角行列記憶部に蓄積する対角化部と、をさらに備えた請求項1から請求項9のいずれか記載の固有値分解装置。
【請求項11】
前記行列分割部は、対称3重対角行列を略半分の2個の対称3重対角行列に分割する、請求項1から請求項10のいずれか記載の固有値分解装置。
【請求項12】
前記固有値算出部は、対称3重対角行列Tの全ての固有値を算出する、請求項1から請求項11のいずれか記載の固有値分解装置。
【請求項13】
前記固有ベクトル算出部は、対称3重対角行列Tの全ての固有ベクトルを算出する、請求項12記載の固有値分解装置。
【請求項14】
前記対称行列Aは、顔画像データを示すベクトルの平均から算出された共分散行列である、請求項10記載の固有値分解装置。
【請求項15】
対称3重対角行列Tが記憶される対角行列記憶部と、行列分割部と、固有値分解部と、前記固有値分解部によって固有値分解された固有値と固有ベクトルからなる直交行列の一部の要素である行列要素とが記憶される固有値分解記憶部と、前記対称3重対角行列Tの固有値が記憶される固有値記憶部と、固有値算出部と、前記対称3重対角行列Tの固有ベクトルが記憶される固有ベクトル記憶部と、固有ベクトル算出部とを備えた固有値分解装置で用いられる固有値分解方法であって、
前記行列分割部が、前記対角行列記憶部から前記対称3重対角行列Tを読み出し、当該対称3重対角行列Tを2個の対称3重対角行列に分割して前記対角行列記憶部に蓄積し、その対称3重対角行列を2個の対称3重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各対称3重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割ステップと、
前記固有値分解部が、前記あらかじめ決められた大きさ以下の各対称3重対角行列を前記対角行列記憶部から読み出し、前記各対称3重対角行列に対して固有値分解を行い、前記各対称3重対角行列の固有値と、前記各対称3重対角行列の固有ベクトルからなる直交行列の行列要素とを少なくとも算出し、当該固有値と当該行列要素とを前記固有値分解記憶部に蓄積する固有値分解ステップと、
前記固有値算出部が、各対称3重対角行列の固有値と、行列要素とを前記固有値分解記憶部から読み出し、前記固有値と、前記行列要素とから、分割元の対称3重対角行列の固有値と、分割元の対称3重対角行列の行列要素とを算出して前記固有値分解記憶部に蓄積し、その分割元の対称3重対角行列の固有値と、行列要素とを算出する処理を、対称3重対角行列Tの少なくとも1個の固有値を算出するまで繰り返し、前記対称3重対角行列Tの少なくとも1個の固有値を前記固有値記憶部に蓄積する固有値算出ステップと、
前記固有ベクトル算出部が、前記対角行列記憶部から前記対称3重対角行列Tを読み出し、前記固有値記憶部から前記対称3重対角行列Tの固有値を読み出し、前記対称3重対角行列Tとその固有値とから、ツイスト分解法を用いて前記対称3重対角行列Tの少なくとも1個の固有ベクトルを算出して前記固有ベクトル記憶部に蓄積する固有ベクトル算出ステップと、を備えた固有値分解方法。
【請求項16】
コンピュータに、
対称3重対角行列Tが記憶される対角行列記憶部から前記対称3重対角行列Tを読み出し、当該対称3重対角行列Tを2個の対称3重対角行列に分割して前記対角行列記憶部に蓄積し、その対称3重対角行列を2個の対称3重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各対称3重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割ステップと、
前記あらかじめ決められた大きさ以下の各対称3重対角行列を前記対角行列記憶部から読み出し、前記各対称3重対角行列に対して固有値分解を行い、前記各対称3重対角行列の固有値と、前記各対称3重対角行列の固有ベクトルからなる直交行列の行列要素とを少なくとも算出し、当該固有値と当該行列要素とを固有値分解記憶部に蓄積する固有値分解ステップと、
各対称3重対角行列の固有値と、行列要素とを前記固有値分解記憶部から読み出し、前記固有値と、前記行列要素とから、分割元の対称3重対角行列の固有値と、分割元の対称3重対角行列の行列要素とを算出して前記固有値分解記憶部に蓄積し、その分割元の対称3重対角行列の固有値と、行列要素とを算出する処理を、対称3重対角行列Tの少なくとも1個の固有値を算出するまで繰り返し、前記対称3重対角行列Tの少なくとも1個の固有値を固有値記憶部に蓄積する固有値算出ステップと、
前記対角行列記憶部から前記対称3重対角行列Tを読み出し、前記固有値記憶部から前記対称3重対角行列Tの固有値を読み出し、前記対称3重対角行列Tとその固有値とから、ツイスト分解法を用いて前記対称3重対角行列Tの少なくとも1個の固有ベクトルを算出して固有ベクトル記憶部に蓄積する固有ベクトル算出ステップと、を実行させるためのプログラム。
発明の詳細な説明 【技術分野】
【0001】
本発明は、固有値分解を行う固有値分解装置等に関する。
【背景技術】
【0002】
固有値分解は、例えば、分子軌道法による化学計算、統計計算、情報検索等の非常に広い分野で利用されている。また、固有値分解の方法も知られている(例えば、非特許文献1参照)。

【非特許文献1】J.J.M.Cuppen,「A divide and conquer method for the symmetric tridiagonal eigenproblem」、Numerische Mathematik,Vol.36、p.177-195、1981年
【発明の開示】
【発明が解決しようとする課題】
【0003】
近年、これらの応用分野におけるデータ量の増大などに伴って、高速・高精度な固有値分解が求められている。また、固有値分解を並列処理することができる固有値分解装置等の開発が望まれていた。
【0004】
本発明は、上記状況の下になされたものであり、並列性に優れた高速、高精度な固有値分解を実行可能な固有値分解装置等を提供することを目的とする。
【課題を解決するための手段】
【0005】
上記目的を達成するため、本発明による固有値分解装置は、対称3重対角行列Tが記憶される対角行列記憶部と、前記対角行列記憶部から前記対称3重対角行列Tを読み出し、当該対称3重対角行列Tを2個の対称3重対角行列に分割して前記対角行列記憶部に蓄積し、その対称3重対角行列を2個の対称3重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各対称3重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割部と、前記あらかじめ決められた大きさ以下の各対称3重対角行列の固有値と、当該各対称3重対角行列の固有ベクトルからなる直交行列の一部の要素である行列要素とが記憶される固有値分解記憶部と、前記あらかじめ決められた大きさ以下の各対称3重対角行列を前記対角行列記憶部から読み出し、前記各対称3重対角行列に対して固有値分解を行い、前記各対称3重対角行列の固有値と、前記各対称3重対角行列の固有ベクトルからなる直交行列の行列要素とを少なくとも算出し、当該固有値と当該行列要素とを前記固有値分解記憶部に蓄積する固有値分解部と、前記対称3重対角行列Tの固有値が記憶される固有値記憶部と、各対称3重対角行列の固有値と、行列要素とを前記固有値分解記憶部から読み出し、前記固有値と、前記行列要素とから、分割元の対称3重対角行列の固有値と、分割元の対称3重対角行列の行列要素とを算出して前記固有値分解記憶部に蓄積し、その分割元の対称3重対角行列の固有値と、行列要素とを算出する処理を、対称3重対角行列Tの少なくとも1個の固有値を算出するまで繰り返し、前記対称3重対角行列Tの少なくとも1個の固有値を前記固有値記憶部に蓄積する固有値算出部と、前記対称3重対角行列Tの固有ベクトルが記憶される固有ベクトル記憶部と、前記対角行列記憶部から前記対称3重対角行列Tを読み出し、前記固有値記憶部から前記対称3重対角行列Tの固有値を読み出し、前記対称3重対角行列Tとその固有値とから、ツイスト分解法を用いて前記対称3重対角行列Tの少なくとも1個の固有ベクトルを算出して前記固有ベクトル記憶部に蓄積する固有ベクトル算出部と、を備えたものである。
【0006】
このような構成により、まず固有値を算出し、その固有値に基づいてツイスト分解法を用いて固有ベクトルを算出することによって、高速で高精度な固有値分解を実現することができうる。また、並列性にも優れている。さらに、全ての固有値や固有ベクトルを算出する必要がない場合には、必要な範囲で固有値や固有ベクトルを算出することができ、処理負荷を軽減することができる。
【0007】
また、本発明による固有値分解装置では、前記固有ベクトル算出部は、qd型ツイスト分解法により固有ベクトルを算出してもよい。
このような構成により、高速で高精度な固有値分解を実現することができうる。また、並列性にも優れている。
【0008】
また、本発明による固有値分解装置では、前記固有ベクトル算出部は、前記固有値記憶部から前記対称3重対角行列Tの固有値を読み出し、当該固有値のいずれかが負の値である場合に、前記対角行列記憶部から前記対称3重対角行列Tを読み出し、当該対称3重対角行列Tの固有ベクトルを変化させることなく、当該対称3重対角行列Tのすべての固有値が正の値となるように、当該対称3重対角行列Tを正定値化し、正定値化した後の固有値を前記固有値記憶部に蓄積し、正定値化した後の対称3重対角行列Tを前記対角行列記憶部に蓄積する正定値化部と、前記対角行列記憶部からすべての固有値が正の値である対称3重対角行列Tを読み出し、前記固有値記憶部から前記対称3重対角行列Tの正の値の固有値を読み出し、前記対称3重対角行列Tの各要素に関してqd型変換を行うことによって、前記対称3重対角行列Tを上2重対角行列B(+1)及び下2重対角行列B(-1)にコレスキー分解するコレスキー分解部と、前記上2重対角行列B(+1)及び下2重対角行列B(-1)の各要素と、前記対称3重対角行列Tの正の値の固有値とを用いて固有ベクトルを算出して前記固有ベクトル記憶部に蓄積するベクトル算出部と、を備えてもよい。
【0009】
このような構成により、負の値の固有値が含まれる場合には、対称3重対角行列Tを正定値化してからqd型ツイスト分解法を実行することができる。また、LV型ツイスト分解法よりも高速に計算することができうる。
【0010】
また、本発明による固有値分解装置では、前記固有ベクトル算出部は、LV型ツイスト分解法により固有ベクトルを算出してもよい。
このような構成により、数値安定的に固有ベクトルの算出を行うことができうる。
【0011】
また、本発明による固有値分解装置では、前記固有ベクトル算出部は、前記対角行列記憶部から前記対称3重対角行列Tを読み出し、前記固有値記憶部から前記対称3重対角行列Tの固有値を読み出し、前記対称3重対角行列Tの各要素に関してミウラ変換、dLVv型変換、逆ミウラ変換を行うことによって、前記対称3重対角行列Tを上2重対角行列B(+1)及び下2重対角行列B(-1)にコレスキー分解するコレスキー分解部と、前記上2重対角行列B(+1)及び下2重対角行列B(-1)の各要素と、前記対称3重対角行列Tの固有値とを用いて固有ベクトルを算出して前記固有ベクトル記憶部に蓄積するベクトル算出部と、を備えてもよい。
このような構成により、コレスキー分解の処理において、複数の補助変数を用いることによって、数値安定的に固有ベクトルの算出を行うことができうる。
【0012】
また、本発明による固有値分解装置では、前記コレスキー分解部は、複数のコレスキー分解手段を備え、前記複数のコレスキー分解手段が、前記対称3重対角行列Tをコレスキー分解する処理を並列実行してもよい。
このような構成により、コレスキー分解の処理を短時間で行うことができうる。
【0013】
また、本発明による固有値分解装置では、前記ベクトル算出部は、複数のベクトル算出手段を備え、前記複数のベクトル算出手段が、固有ベクトルを算出する処理を並列実行してもよい。
このような構成により、固有ベクトルを算出する処理を短時間で行うことができうる。
【0014】
また、本発明による固有値分解装置では、前記固有値算出部は、複数の固有値算出手段を備え、前記複数の固有値算出手段が、分割元の対称3重対角行列の固有値と、行列要素とを算出する処理を並列実行してもよい。
このような構成により、固有値を算出する処理を短時間で行うことができうる。
【0015】
また、本発明による固有値分解装置では、前記固有値分解部は、複数の固有値分解手段を備え、前記複数の固有値分解手段が、対称3重対角行列に対して固有値分解を行う処理を並列実行してもよい。
このような構成により、固有値分解する処理を短時間で行うことができうる。
【0016】
また、本発明による固有値分解装置では、対称行列Aが記憶される行列記憶部と、前記対称行列Aを前記行列記憶部から読み出し、前記対称行列Aを3重対角化した前記対称3重対角行列Tを算出して前記対角行列記憶部に蓄積する対角化部と、をさらに備えてもよい。
【0017】
このような構成により、任意の対称行列Aについて固有値を算出することができる。また、対称3重対角行列Tの固有ベクトルを用いることにより、対称行列Aの固有ベクトルを算出することもできうる。
【0018】
また、本発明による固有値分解装置では、前記行列分割部は、対称3重対角行列を略半分の2個の対称3重対角行列に分割してもよい。
このような構成により、並列処理を適切に行うことができうる。
【0019】
また、本発明による固有値分解装置では、前記固有値算出部は、対称3重対角行列Tの全ての固有値を算出してもよい。
【0020】
また、本発明による固有値分解装置では、前記固有ベクトル算出部は、対称3重対角行列Tの全ての固有ベクトルを算出してもよい。
【0021】
また、本発明による固有値分解装置では、前記対称行列Aは、顔画像データを示すベクトルの平均から算出された共分散行列であってもよい。
このような構成により、固有値分解装置を顔画像データの認識の処理に用いることができうる。
【発明の効果】
【0022】
本発明による固有値分解装置等によれば、高速で高精度な固有値分解を行うことができる。また、並列性にも優れている。
【発明を実施するための最良の形態】
【0023】
以下、本発明による固有値分解装置について、実施の形態を用いて説明する。なお、以下の実施の形態において、同じ符号を付した構成要素及びステップは同一または相当するものであり、再度の説明を省略することがある。
【0024】
(実施の形態1)
本発明の実施の形態1による固有値分解装置について、図面を参照しながら説明する。
図1は、本実施の形態による固有値分解装置1の構成を示すブロック図である。図1において、本実施の形態による固有値分解装置1は、行列記憶部11と、対角化部12と、対角行列記憶部13と、行列分割部14と、固有値分解部15と、固有値分解記憶部16と、固有値算出部17と、固有値記憶部18と、固有ベクトル算出部19と、固有ベクトル記憶部20とを備える。
【0025】
行列記憶部11では、任意の対称行列Aが記憶される。この対称行列Aは、各要素が実数である実行列である。なお、対称行列Aが記憶されているとは、対称行列Aを示すデータが記憶されている、という意味である。後述する記憶部においても同様である。行列記憶部11は、所定の記録媒体(例えば、半導体メモリや磁気ディスク、光ディスクなど)によって実現されうる。行列記憶部11での記憶は、RAM等における一時的な記憶でもよく、あるいは、長期的な記憶でもよい。行列記憶部11に対称行列Aが記憶される過程は問わない。例えば、記録媒体を介して対称行列Aが行列記憶部11で記憶されるようになってもよく、通信回線等を介して送信された対称行列Aが行列記憶部11で記憶されるようになってもよく、あるいは、キーボードやマウス等の入力デバイスを介して入力された対称行列Aが行列記憶部11で記憶されるようになってもよい。
【0026】
対角化部12は、対称行列Aを行列記憶部11から読み出し、その読み出した対称行列Aを3重対角化した対称3重対角行列Tを算出する。そして、対角化部12は、その算出した対称3重対角行列Tを対角行列記憶部13に蓄積する。対角化部12は、例えば、ハウスホルダー(Householder)変換を必要なだけ繰り返し行う方法や、ランチョス(Lanczos)法を用いた方法、その他の3重対角化法を用いて、対称行列Aを3重対角化する。
【0027】
対角行列記憶部13では、対称3重対角行列Tが記憶される。対角行列記憶部13は、所定の記録媒体(例えば、半導体メモリや磁気ディスク、光ディスクなど)によって実現されうる。対角行列記憶部13での記憶は、RAM等における一時的な記憶でもよく、あるいは、長期的な記憶でもよい。
【0028】
行列分割部14は、対角行列記憶部13から対称3重対角行列Tを読み出し、その対称3重対角行列Tを2個の対称3重対角行列に分割して対角行列記憶部13に蓄積する。行列分割部14は、その対称3重対角行列を2個の対称3重対角行列に分割して対角行列記憶部13に蓄積する処理を、分割後の各対称3重対角行列があらかじめ決められた大きさ以下となるまで再帰的に繰り返す。
【0029】
固有値分解部15は、あらかじめ決められた大きさ以下の各対称3重対角行列を対角行列記憶部13から読み出し、その各対称3重対角行列に対して固有値分解を行い、各対称3重対角行列の固有値と、その各対称3重対角行列の固有ベクトルを算出する。固有値分解部15は、例えば、2分法と逆反復法を組み合わせた方法、MR^3法、QR法等を用いて固有値分解を行ってもよい。ここで、MR^3法による固有値分解を行う場合には、LAPACK3.0(http://www.netlib.org/lapack/)で提供されているDSTEGRを用いてもよい。また、QR法による固有値分解を行う場合には、LAPACK3.0で提供されているDSTEQRを用いてもよい。これらの固有値分解の方法については、すでに公知であり、その詳細な説明を省略する。固有値分解部15は、算出した固有値を固有値分解記憶部16に蓄積する。また、固有値分解部15は、算出した固有ベクトルからなる直交行列の一部の要素である行列要素も固有値分解記憶部16に蓄積する。固有ベクトルからなる直交行列とは、固有ベクトルを各列とする直交行列のことである。行列要素の詳細については後述する。
【0030】
固有値分解記憶部16では、固有値分解部15によって固有値分解された固有値と、前述の行列要素とが記憶される。また、固有値分解記憶部16では、固有値算出部17によって算出された固有値や行列要素も記憶される。固有値分解記憶部16は、所定の記録媒体(例えば、半導体メモリや磁気ディスク、光ディスクなど)によって実現されうる。固有値分解記憶部16での記憶は、RAM等における一時的な記憶でもよく、あるいは、長期的な記憶でもよい。
【0031】
固有値算出部17は、固有値分解部15によって算出された固有値と、行列要素とを固有値分解記憶部16から読み出し、その固有値と行列要素とから、分割元の対称3重対角行列の固有値と、分割元の対称3重対角行列の行列要素とを算出して固有値分解記憶部16に蓄積する。固有値算出部17は、その分割元の対称3重対角行列の固有値と、行列要素とを算出する処理を、対称3重対角行列Tの固有値を算出するまで再帰的に繰り返す。そして、固有値算出部17は、算出した対称3重対角行列Tの固有値を固有値記憶部18に蓄積する。
【0032】
固有値記憶部18では、対称3重対角行列Tの固有値が記憶される。固有値記憶部18は、所定の記録媒体(例えば、半導体メモリや磁気ディスク、光ディスクなど)によって実現されうる。固有値記憶部18での記憶は、RAM等における一時的な記憶でもよく、あるいは、長期的な記憶でもよい。
【0033】
固有ベクトル算出部19は、対角行列記憶部13から対称3重対角行列Tを読み出し、固有値記憶部18から対称3重対角行列Tの固有値を読み出す。そして、固有ベクトル算出部19は、対称3重対角行列Tとその固有値とから、ツイスト分解法を用いて対称3重対角行列Tの固有ベクトルを算出して固有ベクトル記憶部20に蓄積する。固有ベクトル算出部19は、LV型ツイスト分解法により固有ベクトルを算出してもよく、あるいは、qd型ツイスト分解法により固有ベクトルを算出してもよい。なお、具体的なツイスト分解の手法については後述する。
【0034】
固有ベクトル記憶部20では、対称3重対角行列Tの固有ベクトルが記憶される。固有ベクトル記憶部20は、所定の記録媒体(例えば、半導体メモリや磁気ディスク、光ディスクなど)によって実現されうる。固有ベクトル記憶部20での記憶は、RAM等における一時的な記憶でもよく、あるいは、長期的な記憶でもよい。
【0035】
なお、行列記憶部11、対角行列記憶部13、固有値分解記憶部16、固有値記憶部18、固有ベクトル記憶部20の任意の2以上の記憶部は、同一の記録媒体によって実現されてもよく、あるいは、別々の記録媒体によって実現されてもよい。前者の場合には、例えば、対称行列Aの記憶されている領域が行列記憶部11となり、対称3重対角行列T等の記憶されている領域が対角行列記憶部13となる。
【0036】
また、行列記憶部11、対角行列記憶部13、固有値分解記憶部16、固有値記憶部18、固有ベクトル記憶部20の各記憶部は、2以上の記録媒体から構成されてもよい。
【0037】
次に、本実施の形態による固有値分解装置1の動作について、図2のフローチャートを用いて説明する。
(ステップS101)対角化部12は、行列記憶部11で記憶されている対称行列Aを読み出し、その行列Aを3重対角化して対称3重対角行列Tを算出して対角行列記憶部13に蓄積する。
【0038】
(ステップS102)行列分割部14、固有値分解部15、固有値算出部17によって対称3重対角行列Tの固有値が算出され、固有値記憶部18に蓄積される。この処理の詳細については後述する。
【0039】
(ステップS103)固有ベクトル算出部19は、対角行列記憶部13から対称3重対角行列Tを読み出し、固有値記憶部18から対称3重対角行列Tの固有値を読み出し、対称3重対角行列Tの固有ベクトルを算出して固有ベクトル記憶部20に蓄積する。この処理の詳細については後述する。
【0040】
このようにして、対称3重対角行列Tの固有値分解が終了する。ここで、対称行列Aの固有値は対称3重対角行列Tの固有値に等しいため、対称行列Aの固有値も算出されたことになる。また、後述するように、所定の変換を行うことによって、対称行列Aの固有ベクトルも対称3重対角行列Tの固有ベクトルから容易に算出することができる。
【0041】
次に、図2のフローチャートのステップS101,S102の処理、すなわち、対称3重対角行列Tの固有値を算出するまでの処理について、より詳細に説明する。まず、図2のフローチャートのステップS102の処理について、図3のフローチャートを用いて説明する。
【0042】
(ステップS201)行列分割部14は、対角行列記憶部13から対称3重対角行列Tを読み出し、その対称3重対角行列Tを2個の対称3重対角行列に分割して対角行列記憶部13に蓄積する。行列分割部14は、その対称3重対角行列を2個の対称3重対角行列に分割する処理を、分割後の各対称3重対角行列があらかじめ決められた大きさ以下となるまで繰り返す。この処理の詳細については後述する。
【0043】
(ステップS202)固有値分解部15は、対角行列記憶部13で記憶されているあらかじめ決められた大きさ以下の各対称3重対角行列について、固有値分解を行う。固有値分解部15は、固有値分解の結果である固有値と、固有ベクトルからなる直交行列の一部の要素である行列要素とを固有値分解記憶部16に蓄積する。
【0044】
(ステップS203)固有値算出部17は、対称3重対角行列の固有値と、行列要素とを固有値分解記憶部16から読み出し、分割元の対称3重対角行列の固有値と、分割元の対称3重対角行列の行列要素とを算出して固有値分解記憶部16に蓄積する。固有値算出部17は、その分割元の対称3重対角行列の固有値と、行列要素とを算出する処理を、対称3重対角行列Tの固有値を算出するまで繰り返し、対称3重対角行列Tの固有値を固有値記憶部18に蓄積する。このようにして、固有値を算出する処理が終了する。
【0045】
[対称行列Aの対角化]
対称行列Aは、例えば、ハウスホルダー変換等を用いて、以下に示されるように3重対角化することができることが知られている。ここで、対称3重対角行列Tは、行の数と列の数とが一致する正方行列である。
【0046】
【数1】
JP0005017666B2_000002t.gif

【0047】
したがって、対角化部12は、上記のようにして、行列記憶部11から対称行列Aを読み出し、対称3重対角行列Tを算出することができる。その算出された対称3重対角行列Tは、対角行列記憶部13で記憶される(ステップS101)。
【0048】
[対称3重対角行列Tの分割]
対称3重対角行列Tを分割する処理について説明する。次式で示されるように、n×nの対称3重対角行列Tが与えられた場合に、それらを2個の対称3重対角行列T,Tと、その他の要素(下記のθとβ)とに分けることができる。
【0049】
【数2】
JP0005017666B2_000003t.gif
ここで、nは、1<n<nとなる整数である。例えば、8×8の対称3重対角行列Tを半分の大きさに分割する具体的な処理は、次のようになる。【数3】
JP0005017666B2_000004t.gif

【0051】
ただし、β=tであり、θは計算誤差が少なくなるように任意に設定することができる。また、上記各行列において、明記されている以外の行列の要素は0である。なお、上記各行列の明記されている各要素は0以外であってもよく、0であってもよい。ここで、並列処理を適切に実行するためには、nを、n/2を超えない最大の整数にとるか、あるいは、n/2を下まわらない最小の整数にとることが好適である(このnの値を用いて行列を2個の行列に分割する場合を、「行列を略半分の2個の行列に分割する」と呼ぶことにする)。
【0052】
【数4】
JP0005017666B2_000005t.gif
本実施の形態では、nを上式のようにとるものとする。なお、nの値は、前述のように、1<n<nの範囲内で任意であることは言うまでもない。したがって、行列分割部14は、上述のようにして、対角行列記憶部13が記憶している対称3重対角行列Tを読み出し、その対称3重対角行列Tを2個の対称3重対角行列T,Tと、2個の値β、θに分割することができ、その分割の処理を繰り返すことができる。行列分割部14は、分割後の2個の対称3重対角行列T,Tと、2個の値β、θとを、対角行列記憶部13に蓄積する。なお、2個の値β、θは、どの分割で発生した値であるのかがわかるように対角行列記憶部13で記憶されることが好適である。図4は、図3のフローチャートのステップS201における行列分割部14による行列を分割する処理を示すフローチャートである。
【0054】
(ステップS301)行列分割部14は、カウンタIを「1」に設定する。
(ステップS302)行列分割部14は、対角行列記憶部13からI番目の分割をしていない対称3重対角行列を読み出し、その対称3重対角行列を2個の対称3重対角行列と、その他の要素とに分割する。そして、行列分割部14は、分割した2個の対称3重対角行列と、その他の要素とを対角行列記憶部13に蓄積する。
【0055】
(ステップS303)行列分割部14は、I番目の分割を行っていない対称3重対角行列が、対角行列記憶部13で記憶されているかどうか判断する。そして、I番目の分割を行っていない対称3重対角行列が、対角行列記憶部13で記憶されている場合には、ステップS302に戻り、そうでない場合には、ステップS304に進む。
【0056】
(ステップS304)行列分割部14は、I番目の分割を行った対称3重対角行列の大きさがあらかじめ決められている大きさ以下かどうか判断する。行列分割部14は、例えば、目的とする行列の大きさ(例えば、25×25など)の記憶されている図示しない記録媒体から、その目的とする行列の大きさを読み出し、対角行列記憶部13で記憶されている、I番目の分割後の対称3重対角行列がその大きさ以下であるかどうかを判断してもよい。そして、I番目の分割を行った対称3重対角行列の大きさがあらかじめ決められている大きさ以下である場合には、対称3重対角行列を分割する処理は終了となり、そうでない場合には、ステップS305に進む。
(ステップS305)行列分割部14は、カウンタIを1だけインクリメントする。そして、ステップS302に戻る。
【0057】
なお、この図4のフローチャートでは、上述のように、行列分割部14が各行列を、略半分の2個の行列に分割するため、I番目の分割後の各対称3重対角行列の大きさがほとんど同じである場合について説明したが、行列分割部14が各行列を、略半分の2個の行列に分割しない場合には、行列分割部14は、分割後の各行列があらかじめ決められた大きさ以下となるように、分割を繰り返すものとする。
【0058】
また、図4のフローチャートでは、ステップS304において、行列分割部14が、対角行列記憶部13で記憶されている分割後の行列の大きさをあらかじめ決められた大きさと比較する処理を実行する場合について説明したが、これは一例であって、行列分割部14は、ステップS304において、それ以外の処理を行ってもよい。例えば、上述のように、行列分割部14が各行列を、略半分の2個の行列に分割する場合には、元の対称3重対角行列Tの大きさを知ることができれば、何番目の分割で目的とする行列の大きさとなるのかを知ることができる。したがって、N番目(Nは1以上の整数)の分割で目的とする行列の大きさとなる場合には、ステップS304において、IがNであるかどうかを比較し、NでなければステップS305に進み、Nであれば一連の処理を終了するようにしてもよい。
【0059】
図5は、行列の分割について説明するための図である。まず、行列分割部14は、1番目の分割として、対称3重対角行列Tを、対称3重対角行列Tと、対称3重対角行列Tとに分割する(ステップS301,S302)。対称3重対角行列Tは、1個しかないため、行列分割部14は、1番目の分割をしていない行列がないと判断する(ステップS303)。また、対称3重対角行列T等はあらかじめ決められた大きさ以下の行列でないとすると(ステップS304)、行列分割部14は、2番目の分割として、対称3重対角行列Tを、対称3重対角行列T11と、対称3重対角行列T12とに分割する(ステップS305,S302)。この場合には、2番目の分割を行っていない対称3重対角行列Tが存在するため、行列分割部14は、対称3重対角行列Tも、対称3重対角行列T21と、対称3重対角行列T22とに分割する(ステップS303,S302)。このようにして、分割後の各対称3重対角行列が目的とする大きさ以下になるまで、行列を分割する処理が繰り返される。なお、図5では、対称3重対角行列以外の要素(上述のβ、θ)については省略している。また、図5において、各対称3重対角行列は、正方行列である。
【0060】
[分割された行列の固有値分解]
行列Tが対称行列であるとすると、行列Tの固有値分解は、次のようになる。
【数5】
JP0005017666B2_000006t.gif

【0061】
ここで、Dは、行列Tと同じ大きさの対角行列である。Dの各対角成分はTの固有値である。また、Qは、行列Tの固有ベクトルからなる直交行列である。
【0062】
したがって、固有値分解部15は、対角行列記憶部13から分割後の各対称3重対角行列を読み出し、上記のようにして、固有値分解を行う(ステップS202)。固有値分解の方法として、例えば、2分法と逆反復法を組み合わせた方法、MR^3法、QR法等を用いてもよいことは前述の通りである。図5で示されるように対称3重対角行列の分割が行われた場合には、固有値分解部15は、各対称3重対角行列T111,T112,T121,T122,…,T222に対して固有値分解を行うことになる。なお、固有値分解部15は、固有値分解の結果得られた各固有値と、固有値分解の結果得られた直交行列の一部の要素である行列要素とを固有値分解記憶部16に蓄積する。ここで、行列要素が、直交行列のどの要素を含むのかについては後述する。なお、固有値分解部15が固有値と行列要素とを固有値分解記憶部16に蓄積するのではなく、固有値と固有ベクトルからなる直交行列とを固有値分解記憶部16に蓄積してもよい。その場合でも、固有ベクトルからなる直交行列を蓄積することによって、少なくとも、行列要素が蓄積されたことになる。
【0063】
なお、後述するように、固有値を算出するときに用いるために必要なのは、直交行列Qの全体ではなく、その直交行列Qの第1行と最終行である。したがって、この固有値分解において、直交行列Qの全体を求めるのではなく、直交行列Qの第1行と最終行を求めるようにしてもよい。そのためには、例えば、LAPACK3.0で提供されているQR法の固有値分解ルーチンであるDSTEQRにおいて、直交行列Qの算出における初期値として、単位行列Iを用いる代わりに、次式で示される行列Yを用いてもよい。
【数6】
JP0005017666B2_000007t.gif

【0064】
このようにすることで、DSTEQRで直交行列Qの第1行と最終行を求めるための計算領域を削減することができる。直交行列Qそのものを算出する場合には、O(n)の作業領域が必要であるが、上記の行列Yを用いて直交行列Qの第1行と最終行を算出する場合には、O(n)の作業領域でよいことになる。その結果、計算量も少なくなり、固有値分解の処理速度も向上する。このように、固有値分解部15は、分割後の対称3重対角行列に対して、直交行列Qそのものを算出する固有値分解を行ってもよく、あるいは、直交行列の一部の成分(ここでは、第1行と最終行)を算出する処理を行ってもよい。この後者の処理も固有値分解と呼ぶことにする。このように、固有値分解部15は、あらかじめ決められた大きさ以下の各対称3重対角行列に対して固有値分解を行い、各対称3重対角行列の固有値と、各対称3重対角行列の固有ベクトルからなる直交行列の行列要素(ここでは、第1行と最終行)とを少なくとも算出し、その固有値と、その行列要素とを固有値分解記憶部16に蓄積してもよい。
【0065】
[固有値の算出]
まず、図6で示されるように、分割された2個の対称3重対角行列T,Tから、分割元の対称3重対角行列Tの固有値等を算出する処理について説明する。対称3重対角行列T,Tは、次のように固有値分解されていたとする。ここで、対称3重対角行列T,Tは、両者共に正方行列であるとする。
【数7】
JP0005017666B2_000008t.gif

【0066】
この場合に、分割元の対称3重対角行列Tは、次のようになる。
【数8】
JP0005017666B2_000009t.gif

【0067】
上記の式において、β、θは、対称3重対角行列の分割の処理において説明した値である。ここで、行列(D12+ρzz)の固有値分解が得られれば、行列Tの固有値分解が得られたことになる。なぜなら、行列(D12+ρzz)の固有値分解D12+ρzz=QDQが求まると、次式のようになり、行列Tの固有値分解が完了するからである。
【数9】
JP0005017666B2_000010t.gif

【0068】
次に、行列(D12+ρzz)の固有値を算出する方法について説明する。行列(D12+ρzz)の固有値をλ、固有ベクトルをqとすると、
【数10】
JP0005017666B2_000011t.gif
となる。
【0069】
(1)zの成分に0が含まれず、D12の対角成分に同じ値が含まれない場合
まず、zの成分に0が含まれず、D12の対角成分に同じ値が含まれない場合について説明する。その場合には、上記の式を、次に示すように順次、変形することができる。
【数11】
JP0005017666B2_000012t.gif

【0070】
ここで、
【数12】
JP0005017666B2_000013t.gif
とすると、次の固有方程式を得ることができる。固有値λ(i=1,2,…,n)は、その固有方程式n個の解である。下記のf(λ)は、単調増加であり、かつf(λ)=0の解の範囲が既知であるため、ニュートン法を用いて計算することができる。
【数13】
JP0005017666B2_000014t.gif

【0071】
なお、D=diag(λ,λ,…,λ)となる。行列Qは、行列(D12+ρzz)の固有ベクトルからなる直交行列であるため、Q=(q,q,…,q)となる。したがって、固有ベクトルqを求めると、行列Qを求めることができる。なお、固有ベクトルqは、次式で与えられる。
【数14】
JP0005017666B2_000015t.gif

【0072】
(2)zの成分に0が含まれる場合
zの成分に0が含まれる場合には、行列の次数を下げることができ、その次数を下げた行列を用いて固有値、及び固有ベクトルを求めることができる。この行列の次数を下げる操作をデフレーションと呼ぶ。
【0073】
=0の場合には、次式で示すように、(d,e)が一組の固有値と固有ベクトルとなる。ただし、eは、j番目の成分が1であり、その他の成分が0である単位ベクトルである。
【数15】
JP0005017666B2_000016t.gif

【0074】
また、行列(D12+ρzz)の固有方程式は、次式のようになる。
【数16】
JP0005017666B2_000017t.gif

【0075】
したがって、z=0の場合の固有値は、dと、行列(D12+ρzz)の第j行、第j列を除いた小行列について、上記式1の固有方程式を解いた値である。すなわち、固有値は、dと、次の固有方程式を解いた値となる。なお、次の固有方程式を解いた固有値を、固有値λ(i=1,2,…,j-1,j+1,…,n)とする。
【数17】
JP0005017666B2_000018t.gif

【0076】
固有ベクトルは、行列(D12+ρzz)の第j行、第j列を除いた小行列について上記式2で算出したn-1次元の固有ベクトルの第(j-1)行と、第j行との間に0を挿入した次式のn次元のベクトルと、上述のeとなる。
【数18】
JP0005017666B2_000019t.gif

【0077】
このように、zの成分に0が含まれる場合には、行列(D12+ρzz)よりも次元の小さい行列について固有値分解を行うことになるため、計算を高速化することが可能であり、結果として、計算時間が短縮されうる。
【0078】
(3)D12の対角成分に同じ値が含まれる場合
j+1=dj+2=…=dの場合にも、上記説明と同様のデフレーションを行うことによって、固有値と固有ベクトルとを算出することができる。
【0079】
【数19】
JP0005017666B2_000020t.gif
に対して、次の2個の直交行列を考える。
【数20】
JP0005017666B2_000021t.gif

【0080】
すると、dj+1=dj+2=…=dより、
【数21】
JP0005017666B2_000022t.gif
となる。したがって、次式が得られる。
【数22】
JP0005017666B2_000023t.gif

【0081】
ここで、
【数23】
JP0005017666B2_000024t.gif
の固有値と固有ベクトルとは、上記のzの成分に0が含まれるデフレーションの場合と同様にして求めることができる。また、
【数24】
JP0005017666B2_000025t.gif
は、符号の任意性を除いて一意に定まる。ここで、その算出方法について説明する。まず、次式を満たす直交行列について考える。
【0082】
【数25】
JP0005017666B2_000026t.gif
この直交行列を次式のようにとると、上記式を満たすことになる。
【数26】
JP0005017666B2_000027t.gif

【0083】
ここで、上記の式の*の値は、次式から求まる。
【数27】
JP0005017666B2_000028t.gif

【0084】
したがって、
【数28】
JP0005017666B2_000029t.gif
となる。その結果、
【数29】
JP0005017666B2_000030t.gif
も、符号の任意性を除いて求めることができる。このようにして、dj+1=dj+2=…=dの場合にも行列(D12+ρzz)の固有値と固有ベクトルとを求めることができる。すなわち、固有値は、dj+2,dj+3,…,dと、次の固有方程式を解いた値となる。なお、次の固有方程式を解いた固有値を、固有値λ(i=1,2,…,j+1,k+1,…,n)とする。
【数30】
JP0005017666B2_000031t.gif

【0085】
固有ベクトルは、次のようになる。
【数31】
JP0005017666B2_000032t.gif

【0086】
なお、上記(2)zの成分に0が含まれる場合、(3)D12の対角成分に同じ値が含まれる場合において、zの成分が0であるかどうかの判断や、D12の対角成分に同じ値が含まれているかどうかの判断は、微少な誤差を考慮して行ってもよい。例えば、次式が満たされる場合には、zの成分が0であると判断してもよく、また、D12の対角成分が同じ値であると判断してもよい。
【数32】
JP0005017666B2_000033t.gif

【0087】
なお、εは、マシン・イプシロンに定数(例えば、2,4,8などの実数であり、計算誤差を少なくするための値を選択することが好適である)を掛けた値を用いてもよく、他の適切な値であってもよい。
【0088】
このようにして、(1)zの成分に0が含まれず、D12の対角成分に同じ値が含まれない場合、(2)zの成分に0が含まれる場合、(3)D12の対角成分に同じ値が含まれる場合のそれぞれについて、固有値分解D12+ρzz=QDQが求まることになる。したがって、行列Qと行列Q12から行列Qを求めることができ、行列Dも求まっているため、分割元の対称3重対角行列Tの固有値分解を行うことができる。このように、対象となる行列を2個の副行列に分割し、その分割した後の各副行列について固有値分解を行い、その後に各副行列の固有値分解を用いて分割元の行列の固有値分解を行う方法は、分割統治法(Divide and Conquer:D&C)と呼ばれている。
【0089】
分割統治法では、固有ベクトルまで算出することになり、その固有ベクトルを算出する処理において行列の計算をしなければならないため、非常に負荷の大きい処理となる。分割統治法において固有値分解をする場合には、例えば、計算時間の95%程度が固有ベクトルを算出するベクトル更新の処理(例えば、上述のように、行列Qと行列Q12から行列Qを算出するために行列を掛け合わせる処理)に費やされることもありうる。
【0090】
一方、分割元の行列Tの固有値のみを算出する場合には、上記処理を簡略化することができる。次に、その方法について説明する。上記の結果より、次のようになる。
【0091】
【数33】
JP0005017666B2_000034t.gif

【0092】
ここで、上記のように、fは行列Tを固有値分解した結果の直交行列の最初の行であり、lは行列Tを固有値分解した結果の直交行列の最後の行である。このfとl、すなわち、直交行列の最初の行と、最後の行とが行列要素となる。なお、直交行列Q12について、後述する列の入れ替えを行った場合には、f12、l12はそれぞれ、その列を入れ替えた後の直交行列Q12の第1行、最終行となる。
【0093】
上記の結果から、図6で示されるように、行列T,Tから分割元の行列Tを構成する場合に、行列Tについて、
【数34】
JP0005017666B2_000035t.gif
をそれぞれ算出することができる。このような処理を繰り返すことにより、最終的に、対称行列Aを3重対角化した対称3重対角行列Tの固有値を算出することができる。
【0094】
なお、上記の(3)D12の対角成分に同じ値が含まれる場合には、式3,式4が次のようになる。
【数35】
JP0005017666B2_000036t.gif

【0095】
ここで、上記式を計算する場合に、
【数36】
JP0005017666B2_000037t.gif
を先に計算して、その後に
【数37】
JP0005017666B2_000038t.gif
を掛けた方が、計算量が少なくなるため、その順番で計算するようにしてもよい。
【0096】
したがって、固有値算出部17は、上述のようにして、固有値分解記憶部16から固有値と行列要素とを読み出し、対角行列記憶部13から行列の分割時に発生したその他の要素に関するβ、θを読み出し、それらを用いることによって、分割元の対称3重対角行列の固有値と行列要素とを算出する。そして、それらを算出する処理を繰り返すことによって、最終的に対称3重対角行列Tの固有値を算出する。図7は、図3のフローチャートのステップS203における固有値算出部17が固有値を算出する処理を示すフローチャートである。
【0097】
(ステップS401)固有値算出部17は、カウンタJを「1」に設定する。
(ステップS402)固有値算出部17は、J番目の固有値の算出は最後の固有値の算出であるかどうか判断する。ここで、最後の固有値の算出とは、対称3重対角行列Tの固有値を算出することである。そして、最後の固有値の算出である場合には、ステップS406に進み、そうでない場合には、ステップS403に進む。
【0098】
(ステップS403)固有値算出部17は、分割元の対称3重対角行列の固有値等を算出する。この処理の詳細については後述する。
(ステップS404)固有値算出部17は、J番目の固有値の算出において、分割元の全ての対称3重対角行列の固有値等を算出したかどうか判断する。そして、J番目の固有値の算出において、分割元の全ての対称3重対角行列の固有値等を算出した場合には、ステップS405に進み、そうでない場合には、ステップS403に戻る。
【0099】
(ステップS405)固有値算出部17は、カウンタJを1だけインクリメントする。そして、ステップS402に戻る。
(ステップS406)固有値算出部17は、対称3重対角行列Tの固有値を算出し、その算出した固有値を固有値記憶部18に蓄積する。このようにして、対称3重対角行列Tの固有値を算出する一連の処理は終了となる。
【0100】
図8は、図7のフローチャートにおけるステップS403の詳細な処理を示すフローチャートである。なお、この図8のフローチャートの処理を実行する前に、対角行列D12の対角成分が昇順、あるいは降順となるように並べ替えておくことが好適である。なお、その場合には、直交行列Q12の列、及びzの成分も、対角行列D12の並び替えと同様に並び替えるものとする。具体的には、対角行列D12の第i成分と、第j成分とを入れ替えた場合には、直交行列Q12の第i列と、第j列とを入れ替え、zの第i行と、第j行とを入れ替えればよい。
【0101】
(ステップS501)固有値算出部17は、分割元の行列の固有値を、式1等を用いて算出する。そして、固有値算出部17は、算出した固有値を固有値分解記憶部16に蓄積する。
【0102】
(ステップS502)固有値算出部17は、ステップS501で算出した固有値を用いて、式2等からqを算出する。このqを算出することにより、qを列ベクトルに有する行列Qを算出したことになる。固有値算出部17は、算出した行列Qを図示しない記録媒体において一次記憶してもよい。
【0103】
(ステップS503)固有値算出部17は、ステップS502で算出した行列Qを用いて、分割元の行列に関する行列要素f、lを算出する。この行列要素は、式3,式4で示されるものである。そして、固有値算出部17は、算出した分割元の行列要素を固有値分解記憶部16に蓄積する。このようにしてステップS403の処理は終了となる。
【0104】
なお、図8のフローチャートにおいて、zの成分に0が含まれず、D12の対角成分に同じ値が含まれない場合には、上記(1)で説明したようにして、固有値、固有ベクトルからなる行列Qを算出することができる。また、zの成分に0が含まれ、D12の対角成分に同じ値が含まれない場合には、上記(2)で説明したようにして、固有値、固有ベクトルからなる行列Qを算出することができる。また、D12の対角成分に同じ値が含まれる場合には、zの成分に0が含まれるかどうかに関わらず、上記(3)で説明したようにして、固有値、固有ベクトルからなる行列Qを算出することができる。上記(3)で説明したようにして固有値等を算出する場合には、前述のように、行列Qを算出するのではなく、行列Qを構成する行列を直接f12等に作用させることによって、f、lを算出するようにしてもよい。
【0105】
また、図8のフローチャートでは、行列Qから、行列要素f、lを一度に算出する場合について示したが、行列要素f、lをその成分ごとに計算してもよい。上記の式3,式4から明らかなように、f12、l12に固有ベクトルqを掛けたものが、行列要素f、lのi番目の成分となる。このような方法で行列要素f、lを計算する場合には、固有ベクトル固有ベクトルqのみをiの値を順次、変えながら記憶することによってその計算を行うことができる。したがって、行列要素f、lの計算において行列Q全体を記憶する必要がなくなり、その計算で用いるメモリ上の作業領域を削減することが可能である。
【0106】
図9は、固有値を算出する処理について説明するための図である。固有値分解記憶部16では、固有値分解部15によって行列T111,T112,…,T222が固有値分解された結果、すなわち、それらの固有値と、行列要素とが記憶されているものとする。まず、固有値算出部17は、固有値等の1番目の算出を開始する(ステップS401,S402)。ここで、固有値等の1番目の算出とは、図9の一番下の行の行列の固有値と、行列要素とから、下から2番目の行の行列の固有値と、行列要素とを算出することである。固有値算出部17は、行列T111と、行列T112との各固有値、及び行列T111と、行列T112との各行列要素を固有値分解記憶部16から読み出す。また、行列T11を行列T111と、行列T112とに分解したときに発生した2個の値β、θを対角行列記憶部13から読み出す。そして、固有値算出部17は、それらの値を用いて、行列T11の固有値を算出して固有値分解記憶部16に蓄積する(ステップS501)。次に、固有値算出部17は、行列T11の固有値を用いて直交行列Qを算出し(ステップS502)、その直交行列Qを用いて行列T11の行列要素を算出して固有値分解記憶部16に蓄積する(ステップS503)。
【0107】
固有値等の1番目の算出はまだ終わっていないため(ステップS404)、固有値算出部17は、上記説明と同様にして、行列T121、行列T122の固有値、行列要素等に基づいて、行列T12の固有値、行列要素を算出する(ステップS402,S403)。このようにして、固有値等の1番目の算出が終了すると、固有値算出部17は、固有値等の2番目の算出、すなわち、行列T11、行列T12の固有値、行列要素等に基づいて行列Tの固有値、行列要素の算出を行う(ステップS404,S405,S402,S403)。その後、固有値算出部17は、同様にして、行列T21、行列T22の固有値、行列要素等に基づいて行列Tの固有値、行列要素の算出を行う(ステップS404,S403)。
【0108】
固有値算出部17が固有値等の2番目の算出を終了すると(ステップS404)、次の固有値等の算出は、対称3重対角行列Tの固有値を求める最後の処理となるため(ステップS405,S402)、固有値算出部17は、固有値の算出のみを行い、その算出した固有値を固有値記憶部18に蓄積する(ステップS406)。このようにして、固有値の算出が終了する。
【0109】
なお、図9の行列T、行列Tの固有値、行列要素等に基づいて行列Tの固有値を算出する場合には、行列要素f、l、あるいは行列要素f、lを用いるだけでよい。したがって、行列T、行列Tの行列要素を算出する場合に、行列要素f、l、あるいは行列要素f、lを計算しなくてもよい。このようにすることで、計算量を削減することができ、処理スピードを向上させることが可能である。
【0110】
次に、固有ベクトル算出部19について説明する。本実施の形態では、固有ベクトル算出部19が(1)qd型ツイスト分解法により固有ベクトルを算出する場合と、(2)LV型ツイスト分解法により固有ベクトルを算出する場合とについて説明する。
【0111】
(1)qd型ツイスト分解法により固有ベクトルを算出する場合
図10は、qd型ツイスト分解法により固有ベクトルを算出する場合の固有ベクトル算出部19の構成を示す図である。図10において、固有ベクトル算出部19は、正定値化部21と、コレスキー分解部22と、ベクトル算出部23とを備える。
【0112】
正定値化部21は、固有値記憶部18から対称3重対角行列Tの固有値を読み出し、その固有値のいずれかが負の値であるかどうか判断する。そして、その固有値のいずれかが負の値である場合に、正定値化部21は、対角行列記憶部13から対称3重対角行列Tを読み出し、その対称3重対角行列Tの固有ベクトルを変化させることなく、対称3重対角行列Tのすべての固有値が正の値となるように、その対称3重対角行列Tを正定値化し、正定値化した後の固有値を固有値記憶部18に蓄積し、正定値化した後の対称3重対角行列Tを対角行列記憶部13に蓄積する。この正定値化の処理を行うのは、qd型ツイスト分解法により固有ベクトルを算出する場合には、対称3重対角行列Tの固有値がすべて正の値でなければならないからである。なお、対称3重対角行列Tを正定値化した行列も、対称3重対角行列Tと呼ぶことにする。また、この正定値化の処理が行われた場合には、コレスキー分解や固有ベクトルを算出する処理は、正定値化後の対称3重対角行列T、固有値に対して行われることになる。この処理の詳細については後述する。ここで、正定値化の行われた場合であっても、固有値記憶部18には、正定値化の前の固有値が記憶されているものとする。その正定値化の前の固有値が、固有値分解装置1が最終的に求める固有値だからである。また、正定値化の行われた場合であって、qd型ツイスト分解法による固有ベクトルの算出のみを行う場合には、対角行列記憶部13で記憶されている対称3重対角行列Tを、正定値化の後の対称3重対角行列Tで上書きしてもよい。一方、LV型ツイスト分解法による固有ベクトルの算出も行う場合には、対角行列記憶部13には、正定値化の前の対称3重対角行列Tと、正定値化の後の対称3重対角行列Tとが記憶されることになる。
【0113】
コレスキー分解部22は、対角行列記憶部13から、すべての固有値が正の値である対称3重対角行列Tを読み出し、固有値記憶部18から、対称3重対角行列Tの正の値の固有値を読み出す。そして、コレスキー分解部22は、対称3重対角行列Tの各要素に関してqd型変換を行うことによって、対称3重対角行列Tを上2重対角行列B(+1)及び下2重対角行列B(-1)にコレスキー分解する。この処理の詳細については後述する。なお、「対称3重対角行列Tを上2重対角行列B(+1)及び下2重対角行列B(-1)にコレスキー分解する」ことには、対称3重対角行列Tの固有値に単位ベクトルを掛けたものを対称3重対角行列Tから引いた行列を、上2重対角行列B(+1)及び下2重対角行列B(-1)にコレスキー分解することが含まれてもよい。
【0114】
ベクトル算出部23は、コレスキー分解部22が算出した上2重対角行列B(+1)及び下2重対角行列B(-1)の各要素と、対称3重対角行列Tの正の値の固有値とを用いて固有ベクトルを算出して固有ベクトル記憶部20に蓄積する。この固有ベクトルを各列に有する行列が、固有値分解の直交行列となる。
【0115】
次に、図2のフローチャートのステップS103の処理、すなわち、対称3重対角行列Tの固有ベクトルを算出する処理について、より詳細に説明する。まず、図2のフローチャートのステップS103の処理について、図11のフローチャートを用いて説明する。
【0116】
(ステップS601)正定値化部21は、固有値記憶部18から対称3重対角行列Tの固有値を読み出し、その固有値のいずれかが負の値であるかどうか判断する。そして、負の値がある場合には、ステップS602に進み、そうでない場合には、ステップS603に進む。
【0117】
(ステップS602)正定値化部21は、対角行列記憶部13から対称3重対角行列Tを読み出し、その対称3重対角行列Tの固有ベクトルを変化させることなく、対称3重対角行列Tのすべての固有値が正の値となるように、その対称3重対角行列Tを正定値化し、正定値化した後の対称3重対角行列Tを対角行列記憶部13に蓄積し、正定値化した後の固有値を固有値記憶部18に蓄積する。
【0118】
(ステップS603)コレスキー分解部22は、対角行列記憶部13から、すべての固有値が正の値である対称3重対角行列Tを読み出し、固有値記憶部18から、対称3重対角行列Tの正の値の固有値を読み出す。なお、正定値化が行われた場合には、正定値化の後の対称3重対角行列T、及び固有値が読み出されることになる。そして、コレスキー分解部22は、対称3重対角行列Tの各要素に関してqd型変換を行うことによって、対称3重対角行列Tを上2重対角行列B(+1)及び下2重対角行列B(-1)にコレスキー分解する。
【0119】
(ステップS604)ベクトル算出部23は、コレスキー分解部22が算出した上2重対角行列B(+1)及び下2重対角行列B(-1)の各要素と、対称3重対角行列Tの正の値の固有値とを用いて固有ベクトルを算出する。
【0120】
(ステップS605)ベクトル算出部23は、算出した固有ベクトルを正規化する。すなわち、ベクトル算出部23は、算出した固有ベクトルのノルムを算出し、固有ベクトルを算出したノルムで割ったものを最終的な固有ベクトルとして固有ベクトル記憶部20に蓄積する。このようにして、対称3重対角行列Tを固有値分解する処理は終了となる。
なお、この後に、行列記憶部11が記憶している対称行列Aの固有ベクトルを算出する処理が行われてもよいが、ここでは省略している。
【0121】
また、固有ベクトル算出部19が算出した固有ベクトルの精度を上げるために、図12で示されるように、固有ベクトルに関する処理を実行してもよく、あるいは、実行しなくてもよい。すなわち、図示しない逆反復法処理部は、固有ベクトル記憶部20から固有ベクトル算出部19が算出した固有ベクトルを読み出し、その固有ベクトルに対して逆反復法の処理を実行し、その結果の固有ベクトルを固有ベクトル記憶部20に蓄積する(ステップS701)。このように、逆反復法の処理を実行することによって、固有ベクトルの精度を向上させることができうる。次に、図示しないグラムシュミット処理部は、高精度が必要かどうか判断する(ステップS702)。この判断は、あらかじめ高精度が必要かどうか設定されている記録媒体等から、その設定を読み出すことによって判断してもよい。そして、図示しないグラムシュミット処理部は、高精度が必要な場合には、固有ベクトル記憶部20から逆反復法の処理の実行された固有ベクトルを読み出し、その固有ベクトルに対してグラムシュミット法の処理を実行し、その結果の固有ベクトルを固有ベクトル記憶部20に蓄積する(ステップS703)。このように、グラムシュミット法の処理を実行することによって、固有ベクトルの直交性を高めることができうる。なお、逆反復法や、グラムシュミット法については、すでに公知であり、その詳細な説明を省略する。
【0122】
[対称3重対角行列Tの正定値化]
行列Tを次のように固有値分解できたとする。
【数38】
JP0005017666B2_000039t.gif

【0123】
その行列Tを次のようにa倍してsIだけシフトしたとする。ここで、Iは行列Tと同じ大きさの単位行列であり、a,sは実数である。
【数39】
JP0005017666B2_000040t.gif

【0124】
すると、シフト後の固有値分解は、次のようになる。
【数40】
JP0005017666B2_000041t.gif

【0125】
このように、行列Tをa倍してsIだけシフトした場合には、各固有値がa倍されてsだけシフトされるが、固有ベクトルは変化しない。したがって、この実数倍とシフトとによって、すべての固有値を正の値にしてから、ツイスト分解法を用いた固有ベクトルの算出を行えばよいことになる。
【0126】
なお、算出された固有値に負の値が含まれる場合には、例えば、(A)最小の固有値によりシフトを行ってもよく、(B)最大の固有値によりシフトを行ってもよい。ここで、固有値には、負の値と正の値とが含まれているものとする。固有値がすべて正の値であれば、シフトを行う必要はなく、固有値がすべて負の値であれば、単に行列Tを-Tとすれば、正定値化することができる。
【0127】
(A)最小の固有値によるシフト
最小の固有値をλmin<0とすると、a=1、s=-λ+λmax×εとする。ここで、εとしては、マシン・イプシロンに定数を掛けた値を用いてもよく、他の適切な値であってもよい。このように、行列Tをシフトすることによって、行列Tのすべての固有値を正の値にすることができる。
【0128】
(B)最大の固有値によるシフト
最大の固有値をλmax>0とすると、a=-1、s=λmax-λmin×εとする。この場合には、次のように行列Tが変換されることになる。
【数41】
JP0005017666B2_000042t.gif

【0129】
上式の右辺の括弧内では、シフトによってすべての固有値が負の値になるが、全体にマイナスを掛けることによって、正定値化することができている。ここで、εは、(A)の場合と同様にマシン・イプシロンに定数(例えば、2,3,4などの実数であり、計算誤差を少なくするための値を選択することが好適である)を掛けた値を用いてもよく、他の適切な値であってもよい。
【0130】
固有値の最小値の絶対値が、固有値の最大値の絶対値よりも小さい場合には、上記(A)のシフトを行い、固有値の最小値の絶対値が、固有値の最大値の絶対値よりも大きい場合には、上記(B)のシフトを行うようにしてもよい。このようにすることで、シフトの量を少なくすることができ、丸め誤差を削減することが可能だからである。
なお、正定値化の後に固有ベクトルを算出する処理の詳細については、次のLV型ツイスト分解法の場合と一緒に説明する。
【0131】
(2)LV型ツイスト分解法により固有ベクトルを算出する場合
図13は、LV型ツイスト分解法により固有ベクトルを算出する場合の固有ベクトル算出部19の構成を示す図である。図13において、固有ベクトル算出部19は、コレスキー分解部31と、ベクトル算出部32とを備える。
【0132】
コレスキー分解部31は、対角行列記憶部13から、対称3重対角行列Tを読み出し、固有値記憶部18から、対称3重対角行列Tの固有値を読み出す。そして、コレスキー分解部31は、対称3重対角行列Tの各要素に関してミウラ変換、dLVv型変換、逆ミウラ変換を行うことによって、対称3重対角行列Tを上2重対角行列B(+1)及び下2重対角行列B(-1)にコレスキー分解する。この処理の詳細については後述する。なお、「対称3重対角行列Tを上2重対角行列B(+1)及び下2重対角行列B(-1)にコレスキー分解する」ことには、対称3重対角行列Tの固有値に単位ベクトルを掛けたものを対称3重対角行列Tから引いた行列を、上2重対角行列B(+1)及び下2重対角行列B(-1)にコレスキー分解することが含まれてもよい。
【0133】
ベクトル算出部32は、コレスキー分解部31が算出した上2重対角行列B(+1)及び下2重対角行列B(-1)の各要素と、対称3重対角行列Tの固有値とを用いて固有ベクトルを算出して固有ベクトル記憶部20に蓄積する。この固有ベクトルを各列に有する行列が、固有値分解の直交行列となる。
【0134】
次に、図2のフローチャートのステップS103の処理、すなわち、対称3重対角行列Tの固有ベクトルを算出する処理について、より詳細に説明する。まず、図2のフローチャートのステップS103の処理について、図14のフローチャートを用いて説明する。
【0135】
(ステップS801)コレスキー分解部31は、対角行列記憶部13から対称3重対角行列Tを読み出し、固有値記憶部18から対称3重対角行列Tの固有値を読み出す。そして、コレスキー分解部31は、対称3重対角行列Tの各要素に関してミウラ変換、dLVv型変換、逆ミウラ変換を行うことによって、対称3重対角行列Tを上2重対角行列B(+1)及び下2重対角行列B(-1)にコレスキー分解する。
【0136】
(ステップS802)ベクトル算出部32は、上2重対角行列B(+1)及び下2重対角行列B(-1)の各要素と、対称3重対角行列Tの固有値とを用いて固有ベクトルを算出する。
【0137】
(ステップS803)ベクトル算出部32は、算出した固有ベクトルを正規化する。すなわち、ベクトル算出部32は、算出した固有ベクトルのノルムを算出し、固有ベクトルを算出したノルムで割ったものを最終的な固有ベクトルとして固有ベクトル記憶部20に蓄積する。このようにして、対称3重対角行列Tを固有値分解する処理は終了となる。
【0138】
なお、この後に、行列記憶部11が記憶している対称行列Aの固有ベクトルを算出する処理が行われてもよいが、ここでは省略している。また、固有ベクトル算出部19が算出した固有ベクトルの精度を上げるために、図12で示されるように、固有ベクトルに関する処理を実行してもよく、あるいは、実行しなくてもよいことは、qd型ツイスト分解法の場合と同様である。
【0139】
[固有値からの固有ベクトルの算出]
固有値から固有ベクトルを算出する処理については、qd型ツイスト分解法を用いる方法と、LV型ツイスト分解法を用いる方法とについて説明する。図11のフローチャートで示されるように、qd型ツイスト分解法を用いて固有ベクトルを算出する場合に、負の固有値があれば、対称3重対角行列Tを正定値化する処理をあらかじめ実行しておく(ステップS601,S602)。qd型ツイスト分解法を用いる場合には、この後の処理において、その正定値化の後の対称3重対角行列T、及び固有値を用いるものとする。
【0140】
対角行列記憶部13では、次のようなm×mの対称3重対角行列Tが記憶されているものとする。
【数42】
JP0005017666B2_000043t.gif

【0141】
また、固有値記憶部18では、固有値λ(1≦i≦m)が記憶されているものとする。このような場合に、対称3重対角行列Tの固有ベクトルを求めることは、次の方程式の固有ベクトルxを求めることになる。
【数43】
JP0005017666B2_000044t.gif

【0142】
本来であれば、上記式の右辺は0になるはずであるが、行列Tの固有値を算出するときに、いくらかの誤差を含むため、固有ベクトルxが真値であれば、上記のように右辺に残差項が存在する。
【0143】
ここで、以下のようにコレスキー分解できたとする。
【数44】
JP0005017666B2_000045t.gif

【0144】
ここで、行列Tの各要素と、行列B(0)の各要素との対応は、次のようになる。
【数45】
JP0005017666B2_000046t.gif

【0145】
なお、行列Tの各要素と、行列B(0)の各要素との対応を次のように書くこともできる。下記の式から明らかなように、行列B(0)の各要素を、行列Tの各要素を用いて順次、算出することができる。
【数46】
JP0005017666B2_000047t.gif

【0146】
上記のコレスキー分解した結果を、次式のように書くことができる。
【数47】
JP0005017666B2_000048t.gif

【0147】
そして、次のようになる。
【数48】
JP0005017666B2_000049t.gif

【0148】
ここで、行列Nをツイスト行列と呼ぶ。また、
【数49】
JP0005017666B2_000050t.gif
であるので、(T-λI)x=γは、
【数50】
JP0005017666B2_000051t.gif
となる。この簡単な式を解くことによって固有ベクトルを算出することができる。具体的には、
【数51】
JP0005017666B2_000052t.gif
のようにわずかな演算で固有ベクトルを算出することができる。なお、あるρ0について、Dρ0=0あるいはDρ0=0であったとしても、行列B(0)の(ρ,ρ)成分であるb2ρ-1(0)と、行列Bの(ρ,ρ+1)成分であるb2ρ(0)を用いて、あるいは、行列Tの成分を用いて、
【数52】
JP0005017666B2_000053t.gif
と固有ベクトルを算出することができる(このように固有ベクトルを算出する処理を例外処理と呼ぶことにする)。なお、残差項のパラメータγ及びkの値は、
【数53】
JP0005017666B2_000054t.gif
の絶対値が最小となるように決定する。したがって、上述のコレスキー分解を求め、ツイスト行列Nを求めることができれば、固有ベクトルを求めることができることになる。そこで、次にコレスキー分解について説明する。
【0149】
図15で示すように、T-λIをコレスキー分解することは、B(0)に対応する{q(0),e(0)}から、B(±1)に対応する{q(±1),e(±1)}を求めることである。
【0150】
(qd型ツイスト分解法)
まず、qd型ツイスト分解法で用いるqd型変換について説明する(図15参照)。
【数54】
JP0005017666B2_000055t.gif
この変換は、さらにstationary qd with shift(stqds)変換
【数55】
JP0005017666B2_000056t.gif
及びreverse-time progressive qd with shift(rpqds)変換
【数56】
JP0005017666B2_000057t.gif
に分けられる。固有値λが既知であれば反復的な計算が不要なため、計算量を少なく抑えることができるが、常に数値安定性と精度が高い方法ではない。それは、stqds変換、rpqds変換は、共に減算による桁落ちが発生する可能性があるからである。例えば、stqds変換において、q(0)+ek-1(0)-λ~ek-1(+1)であれば、q(+1)を求める際に、倍精度計算でもの有効数字がわずか1桁になることもある。その場合には、q(0)(0)/q(+1)を計算すると誤差が生じる。つまり、e(+1)が精度よく計算できないことになる。また、qk+1(+1)を求めるのにe(+1)が必要であり、e(+1)を求めるのにq(+1)が必要であるといったように逐次的な計算が要求されるので、1箇所で発生した桁落ちによる誤差が波及し、さらなる誤差増大の可能性も秘めている。その結果、理論上はq(+1)≠0であるが、誤差蓄積によりq(+1)=0となり、q(0)(0)/q(+1)の計算においてオーバーフローが起こるといった数値不安定な状況も想定される。B(0)の成分{b2k-1(0),b2k(0)}が与えられる、すなわち{q(0),e(0)}が与えられると、λ及びek-1(+1)が一意的に決まるので、この状況は避けることはできない。rpqds変換も同様の性質を持つため、実用的なレベルにまで達したとはいいがたい。LAPACKにおいてFORTRANルーチンDSTEGRとして改良版が公開されているものの欠点は完全に解決されてはいない。
【0151】
なお、このように対角行列Tの各要素を、行列B(0)の各要素に変換して、その行列B(0)の各要素についてqd型変換を行う場合であっても、間接的に対角行列Tの各要素に関してqd型変換を行うことになるため、「対角行列Tの各要素に関してqd型変換を行う」と言うことにする。
また、ここではスタンダードなqd型変換について説明したが、qd型変換はdifferential qd型変換であってもよい。
【0152】
(LV型ツイスト分解法)
次に、LV型ツイスト分解法で用いるミウラ変換、dLVv型変換、逆ミウラ変換について説明する(図15参照)。
【数57】
JP0005017666B2_000058t.gif

【0153】
まず、ミウラ変換について説明する。この変換は、次のように示される。
【数58】
JP0005017666B2_000059t.gif

【0154】
次に、dLVv型変換について説明する。この変換は、さらにstationary discrete Lotka-Volterra with variable step-size(stdLVv)変換
【数59】
JP0005017666B2_000060t.gif
及びreverse-time discrete Lotka-Volterra with variable step-size(rdLVv)変換
【数60】
JP0005017666B2_000061t.gif
に分けられる。ただし、
【数61】
JP0005017666B2_000062t.gif
を満たす範囲内でδ(±1)は任意である。
【0155】
最後に、逆ミウラ変換について説明する。この変換は、次のように示される。
【数62】
JP0005017666B2_000063t.gif

【0156】
このようにして、qd型変換と同様に、コレスキー分解を行うことができる。qd型変換では見られない離散Lotka-Volterra型変換の大きな特徴は、任意パラメータを持つことである。すなわち、λ=1/δ(0)-1/δ(±1)を満たす範囲でδ(n)の値を任意に設定できる。δ(n)を変動させると補助変数u(n)の値も変化するが、桁落ちによる誤差や数値不安定が発生するかどうかは事前に判定できる。この判定は、if文によって実装されてもよい。この場合は、δ(n)を再設定後に再度計算すればよい。また、u(±1)が求まれば逆ミウラ変換によってq(±1)及びe(±1)が独立に計算されるので、誤差が伝播しないという性質を持つ。なお、逆ミウラ変換をミウラ変換、ミウラ変換を逆ミウラ変換と呼んでもよく、stdLVv変換をstLVv変換と呼んでもよく、rdLVv変換をrLVv変換と呼んでもよい。
【0157】
なお、このように対角行列Tの各要素を、行列B(0)の各要素に変換して、その行列B(0)の各要素についてLV型変換を行う場合であっても、間接的に対角行列Tの各要素に関してLV型変換を行うことになるため、「対角行列Tの各要素に関してLV型変換を行う」と言うことにする。
【0158】
ここで、LV型ツイスト分解法による処理のより詳細な処理の一例について説明する。図16~図21は、LV型ツイスト分解法による処理の一例を示すフローチャートである。
【0159】
図16は、コレスキー分解の全体の処理の一例を示す図である。
(ステップS901)コレスキー分解部31は、ミウラ変換を行う。この処理の詳細については後述する。
【0160】
(ステップS902)コレスキー分解部31は、1/δ(±1)を1/δ(0)-λとする。
(ステップS903)コレスキー分解部31は、後述するProcedure1の処理を実行する。
【0161】
(ステップS904)コレスキー分解部31は、後述するProcedure2の処理を実行する。
(ステップS905)コレスキー分解部31は、q(+1),e(+1)がすでに算出されているかどうか判断する。そして、すでに算出されている場合には、コレスキー分解の一連の処理は終了となり、一方、算出されていない場合には、ステップS901に戻る。
【0162】
図17は、図16のフローチャートにおけるステップS903の処理の詳細を示すフローチャートである。
(ステップS1001)コレスキー分解部31は、q(+1),e(+1)がすでに算出されているかどうか判断する。そして、すでに算出されている場合には、終了となり、算出されていない場合には、ステップS1002に進む。
(ステップS1002)コレスキー分解部31は、stdLVv変換等の処理を行う。この処理の詳細については後述する。
【0163】
図18は、図16のフローチャートにおけるステップS904の処理の詳細を示すフローチャートである。
(ステップS1101)コレスキー分解部31は、q(-1),e(-1)がすでに算出されているかどうか判断する。そして、すでに算出されている場合には、終了となり、算出されていない場合には、ステップS1102に進む。
(ステップS1102)コレスキー分解部31は、rdLVv変換等の処理を行う。この処理の詳細については後述する。
【0164】
図19は、図16のフローチャートのステップS901の処理の詳細を示すフローチャートである。
(ステップS1201)コレスキー分解部31は、1/δ(0)の値を決定する。この値は、前述のように任意に決定することができるため、行列を正定値化する値に決定する。例えば、コレスキー分解部31は、1/δ(0)の値を次式のように決定する。
【数63】
JP0005017666B2_000064t.gif

【0165】
なお、なお、ε(0)は、マシン・イプシロンに定数(例えば、2,3,4などの実数であり、計算誤差を少なくするための値を選択することが好適である)を掛けた値を用いてもよく、他の適切な値であってもよい。その後、例えば、ステップS1203等で桁落ちの発生する可能性があると判断された場合に、コレスキー分解部31は、1/δ(0)の値を決めるε(0)の値を、例えば、マシン・イプシロンに掛ける定数を1ずつ大きくしながら設定していってもよい。なお、1/δ(0)の値を決定する方法は、これに限定されないことは言うまでもない。このように、LV型ツイスト分解法では、この1/δ(0)の値を自由に設定することができるため、qd型ツイスト分解法の場合のように、行列Tの正定値化を行わなくてもよいことになる。なお、qd型ツイスト分解法の場合と同様に、行列Tの正定値化の処理を行ってから、LV型ツイスト分解法を実行してもよいことは言うまでもない。
(ステップS1202)コレスキー分解部31は、β1をq(0)-1/δ(0)に設定する。
【0166】
(ステップS1203)コレスキー分解部31は、β1の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップS1204に進み、そうでない場合には、ステップS1201に戻って、1/δ(0)の値を決定する処理を再度実行する。なお、このεも任意に定めることができる。例えば、εとしてマシン・イプシロンを用いてもよい。εの値を大きくすると精度が向上し、εの値を小さくすると精度が低下する。なお、このステップS1203で行っている処理は、桁落ちが発生する可能性について判断する処理である。β1の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。
【0167】
(ステップS1204)コレスキー分解部31は、β2をq(0)/(1+δ(0)(0))に設定する。なお、前述のようにu(0)=0であるため、β2はq(0)に設定されたことになる。
【0168】
(ステップS1205)コレスキー分解部31は、u(0)(1+δ(0)(0))をβ1に設定する。ここで、β1は、ステップS1202においてq(0)-1/δ(0)に設定されているが、前述のようにu(0)=0であるため、β1は、q(0)/(1+δ(0)(0)-1/δ(0)と等しく、これにミウラ変換を行うとu(0)=u2k-1(0)k=1となるからである。なお、u(0)=0から(1+δ(0)(0))=1である。
【0169】
(ステップS1206)コレスキー分解部31は、カウンタkを「1」に設定する。
(ステップS1207)コレスキー分解部31は、α1をe(0)/β1に設定する。前述のように、β1はu2k-1(0)となるから、e(0)/β1にミウラ変換を行うと、α1はδ(0)2k(0)と等しいことになる。
【0170】
(ステップS1208)コレスキー分解部31は、α2をα1+1に設定する。ステップS1207の説明からわかるように、α2は1+δ(0)2k(0)と等しいことになる。
【0171】
(ステップS1209)コレスキー分解部31は、α2の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップS1210に進み、そうでない場合には、ステップS1201に戻って、1/δ(0)の値を決定する処理を再度実行する。なお、このステップS1209で行っている処理は、ステップS1203の処理と同様に、桁落ちが発生する可能性について判断する処理である。α2の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。
【0172】
(ステップS1210)コレスキー分解部31は、α1×β2を算出してu2k(0)(1+δ(0)2k-1(0))に設定する。前述のように、α1はδ(0)2k(0)と等しく、β2はq(0)/(1+δ(0)2k-2(0))であるため、α1×β2にミウラ変換を実行すると、u2k(0)(1+δ(0)2k-1(0))に等しくなるからである。
【0173】
(ステップS1211)コレスキー分解部31は、β2をqk+1(0)/α2に設定する。前述のように、α2は1+δ(0)2k(0)と等しいため、β2は、qk+1(0)/(1+δ(0)2k(0))となり、β2におけるkの値を1だけインクリメントしたことになる。
【0174】
(ステップS1212)コレスキー分解部31は、β1をβ2-1/δ(0)に設定する。前述のように、β2はqk+1(0)/(1+δ(0)2k(0))と等しいため、β1は、qk+1(0)/(1+δ(0)2k(0))-1/δ(0)となり、ミウラ変換を実行すると、β1は、u2k+1(0)となる。したがって、β1におけるkの値を1だけインクリメントしたことになる。
【0175】
(ステップS1213)コレスキー分解部31は、β1の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップS1214に進み、そうでない場合には、ステップS1201に戻って、1/δ(0)の値を決定する処理を再度実行する。なお、このステップS1213で行っている処理は、ステップS1203の処理と同様に、桁落ちが発生する可能性について判断する処理である。β1の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。
【0176】
(ステップS1214)コレスキー分解部31は、α2×β1を算出してu2k+1(0)(1+δ(0)2k(0))に設定する。前述のように、α2は1+δ(0)2k(0)と等しく、β1は、u2k+1(0)に等しいからである。
(ステップS1215)コレスキー分解部31は、カウンタkを1だけインクリメントする。
【0177】
(ステップS1216)コレスキー分解部31は、カウンタkがmに等しいかどうか判断する。そして、mに等しい場合には、一連の処理は終了となり、そうでない場合には、ステップS1207に戻る。
【0178】
このようにして、u2k+1(0)(1+δ(0)2k(0))と、u2k(0)(1+δ(0)2k-1(0))とが算出されることになる。これらの値は、コレスキー分解部31が有する図示しないメモリ等において一時的に記憶されてもよい。
【0179】
図20は、図17のフローチャートのステップS1002の処理の詳細を示すフローチャートである。
(ステップS1301)コレスキー分解部31は、α2を1+δ(+1)(+1)に設定する。なお、前述のようにu(+1)=0であるため、α2は1に設定されたことになる。
【0180】
(ステップS1302)コレスキー分解部31は、β1をu(0)(1+δ(0)(0))/(1+δ(+1)(+1))に設定する。ここで、u(0)(1+δ(0)(0))の値としては、ステップS1205で算出したものを用いる。なお、前述のようにu(+1)=0である。また、このβ1の式にstdLVv変換を実行すると、β1は、u(+1)に設定されたことになる。
【0181】
(ステップS1303)コレスキー分解部31は、カウンタkを「1」に設定する。
(ステップS1304)コレスキー分解部31は、α1をβ1+1/δ(+1)に設定する。
【0182】
(ステップS1305)コレスキー分解部31は、α1の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップS1306に進み、そうでない場合には、図16のフローチャートのProcedure2(ステップS904)に進む。なお、このステップS1305で行っている処理は、ステップS1203の処理と同様に、桁落ちが発生する可能性について判断する処理である。α1の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。
【0183】
(ステップS1306)コレスキー分解部31は、β2をu2k(0)(1+δ(0)2k-1(0))/α1に設定する。ここで、u2k(0)(1+δ(0)2k-1(0))の値としては、ステップS1210で算出したものを用いる。また、このβ2の式にstdLVv変換を実行すると、β2は、δ(+1)2k(+1)に設定されたことになる。
【0184】
(ステップS1307)コレスキー分解部31は、α1×α2を算出する。前述のように、α1はβ1+1/δ(+1)=u2k-1(+1)+1/δ(+1)に等しく、α2は1+δ(+1)2k-2(+1)に等しいため、α1×α2に逆ミウラ変換を実行するとq(+1)に等しくなる。したがって、コレスキー分解部31は、α1×α2を算出してq(+1)に設定する。
【0185】
(ステップS1308)コレスキー分解部31は、α2を1+β2に設定する。前述のように、β2はδ(+1)2k(+1)と等しいため、α2は1+δ(+1)2k(+1)となる。したがって、α2におけるkの値を1だけインクリメントしたことになる。
【0186】
(ステップS1309)コレスキー分解部31は、α2の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップS1310に進み、そうでない場合には、図16のフローチャートのProcedure2(ステップS904)に進む。なお、このステップS1309で行っている処理は、ステップS1203の処理と同様に、桁落ちが発生する可能性について判断する処理である。α2の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。
【0187】
(ステップS1310)コレスキー分解部31は、β1×β2を算出する。前述のように、β1はu2k-1(+1)に等しく、β2はδ(+1)2k(+1)と等しいため、β1×β2に逆ミウラ変換を実行すると、e(+1)に等しくなる。したがって、コレスキー分解部31は、β1×β2を算出してe(+1)に設定する。
【0188】
(ステップS1311)コレスキー分解部31は、β1をu2k+1(0)(1+δ(0)2k(0))/α2に設定する。ここで、u2k+1(0)(1+δ(0)2k(0))の値としては、ステップS1214で算出したものを用いる。また、このβ1の式にstdLVv変換を実行すると、β1は、u2k+1(+1)に設定されたことになる。したがって、β1におけるkの値を1だけインクリメントしたことになる。
【0189】
(ステップS1312)コレスキー分解部31は、カウンタkを1だけインクリメントする。
(ステップS1313)コレスキー分解部31は、カウンタkがmに等しいかどうか判断する。そして、mに等しい場合には、ステップS1314に進み、そうでない場合には、ステップS1304に戻る。
【0190】
(ステップS1314)コレスキー分解部31は、α2×(β1+1/δ(+1))を算出する。これは、ステップS1304でα1を更新した後に、ステップS1307でα1×α2を計算することと等しい。したがって、コレスキー分解部31は、α2×(β1+1/δ(+1))を算出してq(+1)に設定する。このようにして、ミウラ変換された結果にstdLVv変換と、逆ミウラ変換とを実行して、q(+1)、e(+1)を算出する処理が終了する。これらの値は、コレスキー分解部31が有する図示しないメモリ等において一時的に記憶されてもよい。
【0191】
図21は、図18のフローチャートのステップS1102の処理の詳細を示すフローチャートである。
(ステップS1401)コレスキー分解部31は、β1をu2m-1(0)(1+δ(0)2m-2(0))/(1+δ(-1)2m(-1))に設定する。ここで、u2m-1(0)(1+δ(0)2m-2(0))の値としては、ステップS1214で算出したものを用いる。なお、前述のようにu2m(-1)=0である。また、このβ1の式にrdLVv変換を実行すると、β1は、u2m-1(-1)に設定されたことになる。
【0192】
(ステップS1402)コレスキー分解部31は、カウンタkを「m-1」に設定する。
(ステップS1403)コレスキー分解部31は、α1をβ1+1/δ(-1)に設定する。
【0193】
(ステップS1404)コレスキー分解部31は、α1の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップS1405に進み、そうでない場合には、図16のフローチャートのミウラ変換(ステップS901)に戻る。なお、このステップS1404で行っている処理は、ステップS1203の処理と同様に、桁落ちが発生する可能性について判断する処理である。α1の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。
【0194】
(ステップS1405)コレスキー分解部31は、β2をu2k(0)(1+δ(0)2k-1(0))/α1に設定する。ここで、u2k(0)(1+δ(0)2k-1(0))の値としては、ステップS1210で算出したものを用いる。また、このβ2の式にrdLVv変換を実行すると、β2は、δ(-1)2k(-1)に設定されたことになる。
【0195】
(ステップS1406)コレスキー分解部31は、α2を1+β2に設定する。前述のように、β2はδ(-1)2k(-1)と等しいため、α2は1+δ(-1)2k(-1)となる。
【0196】
(ステップS1407)コレスキー分解部31は、α2の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップS1408に進み、そうでない場合には、図16のフローチャートのミウラ変換(ステップS901)に戻る。なお、このステップS1407で行っている処理は、ステップS1203の処理と同様に、桁落ちが発生する可能性について判断する処理である。α2の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。
【0197】
(ステップS1408)コレスキー分解部31は、α1×α2を算出する。前述のように、α1はβ1+1/δ(-1)=u2k+1(-1)+1/δ(-1)に等しく、α2は1+δ(-1)2k(-1)に等しいため、α1×α2に逆ミウラ変換を実行するとqk+1(-1)に等しくなる。したがって、コレスキー分解部31は、α1×α2を算出してqk+1(-1)に設定する。
【0198】
(ステップS1409)コレスキー分解部31は、β1をu2k-1(0)(1+δ(0)2k-2(0))/α2に設定する。ここで、u2k-1(0)(1+δ(0)2k-2(0))の値としては、ステップS1205,S1214で算出したものを用いる。また、このβ1の式にstdLVv変換を実行すると、β1は、u2k-1(-1)に設定されたことになる。したがって、β1におけるkの値を1だけデクリメントしたことになる。
【0199】
(ステップS1410)コレスキー分解部31は、β1×β2を算出する。前述のように、β1はu2k-1(-1)に等しく、β2はδ(-1)2k(-1)と等しいため、β1×β2に逆ミウラ変換を実行すると、e(-1)に等しくなる。したがって、コレスキー分解部31は、β1×β2を算出してe(-1)に設定する。
【0200】
(ステップS1411)コレスキー分解部31は、カウンタkを1だけデクリメントする。
(ステップS1412)コレスキー分解部31は、カウンタkが0に等しいかどうか判断する。そして、0に等しい場合には、ステップS1413に進み、そうでない場合には、ステップS1403に戻る。
【0201】
(ステップS1413)コレスキー分解部31は、β1+1/δ(-1)を算出する。これは、ステップS1403でα1を更新した後に、ステップS1408でα1×α2を計算することと等しい。この場合、α2は1だからである。したがって、コレスキー分解部31は、β1+1/δ(-1)を算出してq(-1)に設定する。このようにして、ミウラ変換された結果にrdLVv変換と、逆ミウラ変換とを実行して、q(-1)、e(-1)を算出する処理が終了する。これらの値は、コレスキー分解部31が有する図示しないメモリ等において一時的に記憶されてもよい。
【0202】
なお、図16のフローチャートの処理によって算出することができるのは1個の固有値に対応する固有ベクトルであるため、全ての固有値に対応する固有ベクトルを算出する場合には、コレスキー分解部31は、図16のフローチャートの処理を固有値の数だけ繰り返すことになる。
【0203】
メモリ消費を抑えるために、補助変数u(n)のための配列は必ずしも用意する必要はない。一方、1+δ(0)(0)のためのメモリ領域を確保し、ミウラ変換、dLVv型変換、逆ミウラ変換のステップにまたがってこの値を利用することで、メモリ消費を抑え、計算量を低減することができる。これにより、誤差も低減される。
【0204】
ここで、誤差について説明する。人間が理想的な状況、すなわち、無限桁の計算をいくらでもできるとすると、qd型ツイスト分解法であっても、LV型ツイスト分解法であっても問題ない。しかしながら、コンピュータで計算を行う場合は注意が必要である。有限桁の計算しか行えないコンピュータ上では、数学的に正しい計算法を使ったとしても必ずしも正しい結果が得られる訳ではない。そればかりか、いつまでたっても計算が終了しないといった思わぬ数値的な問題が発生する場合もある。
【0205】
コンピュータ計算による誤差には、丸め誤差及び桁落ちによる誤差などが知られている。丸め誤差単独では、高々有効桁の最後の桁が真値と比べて異なる程度で大きな誤差にはならない。また、指数部が異なる2つの実数の加算、乗算、除算の少なくとも1つの演算を行えばやはり丸め誤差が生じるが、それ以上の誤差は発生しない。さらに、このような丸め誤差を発生するような操作が繰り返されても、丸めモードがnear(四捨五入)ならば、一方的に切り上げられたり、あるいは切り捨てられたりして誤差が極端に蓄積することは少ない。よって、多くの数値計算法は加算、乗算、除算の少なくとも1つの演算によって発生する丸め誤差を特別注意することは少なく、適切な固有値計算でも結果的に丸め誤差は一様に増大しない。
【0206】
問題となるのは、同符号の実数の減算及び異符号の実数の加算により生じる、桁落ちである。桁落ちによる誤差で値が0となった後、その値による除算を行うと、0が分母にくるような不定形となり計算不可能となる。こうなるといつまでたっても計算が終了しない。減算→除算と計算する部分がqd型ツイスト分解法、LV型ツイスト分解法の両方に存在するので、桁落ち誤差には十分に注意する必要がある。
【0207】
LV型ツイスト分解法では、上述の桁落ちによる誤差を含んでいるかどうかは減算によって得られた値が小さいかどうかで判断できる。qd型ツイスト分解法の場合、桁落ち誤差を含むことが分かったとしても、それを回避することはできない。なぜならば、初期値として{q(0),e(0)}が与えられると、λは一意的に決定され、{q(±1),e(±1)}も一意的に導出されるためである。すなわち、任意パラメータを持たない自由度のない計算法であるためである。
【0208】
それに対して、LV型ツイスト分解法は、自由に設定できるパラメータδ(0)を持つため、補助変数u(n)の値を様々に変化させることができる(図22参照)。すなわち、様々な経路で{q(±1),e(±1)}を計算することができる。よって、桁落ちが発生する場合も回避できる。図19~図21の条件判定によって桁落ちの影響をチェックし、減算によって得られた値の絶対値が小さな数εより大きいという条件が満たされなければ、パラメータδ(0)の設定に戻るというものである。この処理は、条件が満たされるまで繰り返される。なお、精度よりも高速性を重視する場合は、数回条件が満たさなければ(q(+1)=0あるいはq(-1)=0ならば)、例外処理を行ってもよい。
【0209】
なお、前述の行列Tの各要素と、行列B(0)の各要素との対応を示す式を用いると、qd型変換を次のように書き換えることもできる。
【数64】
JP0005017666B2_000065t.gif

【0210】
この変換は、さらにstqds変換
【数65】
JP0005017666B2_000066t.gif
及びrpqds変換
【数66】
JP0005017666B2_000067t.gif
に分けられる。このように、{t2k-1,t2k}から{q(±1),e(±1)}を求める場合には、{t2k-1,t2k}から{q(0),e(0)}を求め、さらにq(0),e(0)}から{q(±1),e(±1)}を求める場合に比べて、計算量が少なくなることとなり、計算速度を向上させることも可能である。
【0211】
なお、この場合には、対角行列Tの各要素について直接qd型変換が行われることになる。したがって、この変換が、「対角行列Tの各要素に関してqd型変換を行う」場合に含まれることは言うまでもない。
【0212】
また、qd型変換と同様に、前述の行列Tの各要素と、行列B(0)の各要素との対応を示す式を用いると、LV型変換におけるミウラ変換を次のように書き換えることもできる。
【数67】
JP0005017666B2_000068t.gif

【0213】
この場合のミウラ変換は、次のように示される。
【数68】
JP0005017666B2_000069t.gif
このように、{t2k-1,t2k}から{u2k-1(0),u2k(0)}を求める場合には、{t2k-1,t2k}から{q(0),e(0)}を求め、さらにq(0),e(0)}から{u2k-1(0),u2k(0)}を求める場合に比べて、計算量が少なくなることとなり、計算速度を向上させることも可能である。
【0214】
なお、この場合には、対角行列Tの各要素について直接LV型変換が行われることになる。したがって、この変換が、「対角行列Tの各要素に関してLV型変換を行う」場合に含まれることは言うまでもない。
【0215】
qd型ツイスト分解法あるいはLV型ツイスト分解法によってコレスキー分解をすることができると、上述のようにツイストされた行列Nを算出することができ、その行列NのN=eを解くことによって、xを算出することができる。ここで、
【数69】
JP0005017666B2_000070t.gif
とxを置き換えることにより、このxを正規化する。このようにして、固有ベクトルxを求めることができる。この固有ベクトルxを用いて、Q=(x,x,...,x)とすることにより、直交行列Qを算出することができる。
【0216】
したがって、次式のように対称3重対角行列Tについて固有値分解が行われたことになる。
【数70】
JP0005017666B2_000071t.gif

【0217】
また、対称3重対角行列Tの固有値分解を行うことができれば、次のように、対称行列Aの固有値分解も行われる。
【数71】
JP0005017666B2_000072t.gif

【0218】
固有ベクトル算出部19は、上述のようにして、固有値記憶部18が記憶している対称3重対角行列Tの固有値と、対角行列記憶部13が記憶している対称3重対角行列Tとを用いて、対称3重対角行列Tの固有ベクトルを算出する。まず、コレスキー分解部31は、対角行列記憶部13から対称3重対角行列Tを読み出し、固有値記憶部18から対称3重対角行列Tの固有値を読み出す。そして、コレスキー分解部31は、図23のフローチャートで示されるように各変換を行い、コレスキー分解の処理を行う。ここでは、LV型ツイスト分解法を用いる場合について説明する。
【0219】
コレスキー分解部31は、まず、対称3重対角行列Tの各要素の値から、q(0),e(0)を求める。そして、コレスキー分解部31は、ミウラ変換を実行することにより、u(0),u(0)等を順次求めていく(ステップS1501)。
【0220】
次に、コレスキー分解部31は、ミウラ変換で得られたu(0)を用いて、stdLVv変換を実行することにより、u(+1)等を順次求めていく(ステップS1502)。また、コレスキー分解部31は、ミウラ変換で得られたu(0)を用いて、rdLVv変換を実行することにより、u2m-1(-1)等を順次求めていく(ステップS1503)。
【0221】
最後に、コレスキー分解部31は、stdLVv変換で得られたu(+1)及びrdLVv変換で得られたu(-1)を用いて、逆ミウラ変換を実行することにより、q(±1)、e(±1)等を順次求めていく(ステップS1504)。q(±1)及びe(±1)が求まると、すなわち、上2重対角行列B(+1)と、下2重対角行列B(-1)が求まると、コレスキー分解部31は、それらをベクトル算出部32に渡す。なお、コレスキー分解部31は、各固有値について、それぞれコレスキー分解を行うものとする。
【0222】
なお、ここでは、説明の便宜上、図23のフローチャートを用いてコレスキー分解の処理を説明したが、図16~図21のフローチャートで示されるように処理を行ってもよいことは言うまでもない。
【0223】
ベクトル算出部32は、固有値記憶部18から固有値を読み出し、式5を用いてkの値を決定する。ベクトル算出部32は、そのkの値を用いて、q(±1)及びe(±1)からツイスト行列Nを求め、N=eを解くことによって、xを算出する(ステップS802)。
【0224】
次に、ベクトル算出部32は、算出したxを正規化して、固有ベクトルxを求めて固有ベクトル記憶部20に蓄積する(ステップS803)。なお、ベクトル算出部32は、各固有値について、xを算出する処理と、正規化する処理とを行うことによって、全ての固有ベクトルを算出することができる。このようにして、固有ベクトルの算出が終了する。
この後、行列記憶部11が記憶している対称行列Aの固有ベクトルを算出してもよい。その算出の方法は、前述の通りである。
【0225】
また、固有値算出部17が算出した固有値、固有ベクトル算出部19が算出した固有ベクトルを出力する図示しない出力部を固有値分解装置1が備えてもよい。ここで、図示しない出力部による出力は、例えば、固有値等の表示デバイス(例えば、CRTや液晶ディスプレイなど)への表示でもよく、固有値等の所定の機器への通信回線を介した送信でもよく、固有値等のプリンタによる印刷でもよく、固有値等の記録媒体への蓄積でもよい。なお、その出力部は、出力を行うデバイス(例えば、表示デバイスやプリンタなど)を含んでもよく、あるいは含まなくてもよい。また、その出力部は、ハードウェアによって実現されてもよく、あるいは、それらのデバイスを駆動するドライバ等のソフトウェアによって実現されてもよい。
【0226】
なお、ここでは、コレスキー分解部31がLV型ツイスト分解法によってコレスキー分解を行う場合について説明したが、コレスキー分解部31は、qd型ツイスト分解法によってコレスキー分解を行ってもよい。例えば、コレスキー分解部31は、算出された固有値の分布を見て、分布が密でない場合には、qd型ツイスト分解法を用いるようにしてもよい。
【0227】
また、上記説明では、行列B(0)が上2重対角行列である場合について説明したが、行列B(0)は下2重対角行列であっても固有値分解を同様に実行することができうる。ただし、行列Tの各要素と、行列B(0)の各要素との対応を示す式が上述のものと異なることになる。
【0228】
また、上記説明では、固有値算出部17が全ての固有値を算出する場合について説明したが、一部の固有値のみを算出するようにしてもよい。例えば、図9において、固有値算出部17は、行列T,Tまでは、全ての固有値を算出する必要があるが、行列T,Tから対称3重対角行列Tの固有値を算出する際に、必要な範囲で固有値を算出してもよい。したがって、固有値算出部17は、対称3重対角行列Tの少なくとも1個の固有値を算出するものであってもよい。その場合には、不必要な固有値を算出することに伴う余分な処理を実行しなくてよいため、処理負荷を軽減することができうる。このように、一部の固有値のみを算出する場合には、例えば、対称3重対角行列Tの固有値を算出する処理に費やす計算時間を約(1+k/m)/2倍にすることができうる。ここで、mは行列サイズであり、kは求める固有値の個数である。
【0229】
また、上記説明では、固有ベクトル算出部19が算出された固有値に対応する全ての固有ベクトルを算出する場合について説明したが、上記説明から明らかなように、固有ベクトルを算出する処理は、固有値ごとに行うことができる。したがって、固有ベクトル算出部19は、対称3重対角行列Tの少なくとも1個の固有ベクトルを算出するものであってもよい。このように、固有ベクトル算出部19は、必要な範囲で固有ベクトルを算出することができ、不必要な固有ベクトルを算出することに伴う余分な処理を実行しなくてよいため、処理負荷を軽減することができうる。
【0230】
次に、本実施の形態による固有値分解装置1における並列処理について説明する。
本実施の形態による固有値分解装置1では、固有値の計算及び固有ベクトルの計算において、並列的に処理を実行することができる。例えば、図24~図26で示されるように、固有値分解装置1において、固有値分解部15、固有値算出部17、コレスキー分解部22、ベクトル算出部23,コレスキー分解部31、ベクトル算出部32において、それぞれ並列処理を行ってもよい。
【0231】
固有値分解部15は、複数の固有値分解手段15a、15bを備え、その複数の固有値分解手段15a、15bが、分割後の対称3重対角行列の固有値と、行列要素とを算出する処理を並列実行してもよい。例えば、図9において、固有値分解手段15aが行列T111,T112,T121,T121の固有値分解を実行し、固有値分解手段15bが行列T211,T212,T221,T222の固有値分解を実行してもよい。
【0232】
固有値算出部17は、複数の固有値算出手段17a、17bを備え、その複数の固有値算出手段17a、17bが、分割元の対称3重対角行列の固有値と、行列要素とを算出する処理を並列実行してもよい。例えば、図9において、固有値算出手段17aが行列T11,T12,T,Tの固有値等を算出する処理を実行し、固有値算出手段17bが行列T21,T22,Tの固有値等を算出する処理を実行してもよい。
【0233】
コレスキー分解部22は、複数のコレスキー分解手段22a、22bを備え、その複数のコレスキー分解手段22a、22bが、対称3重対角行列Tに関してコレスキー分解する処理を並列実行してもよい。例えば、コレスキー分解手段22aが半分の固有値についてコレスキー分解する処理を実行し、コレスキー分解手段22bが残りの半分の固有値についてコレスキー分解する処理を実行してもよい。
【0234】
ベクトル算出部23は、複数のベクトル算出手段23a、23bを備え、その複数のベクトル算出手段23a、23bが、固有ベクトルを算出する処理を並列実行してもよい。例えば、ベクトル算出手段23aがコレスキー分解手段22aによってコレスキー分解された固有値に関して固有ベクトルを算出する処理を実行し、ベクトル算出手段23bがコレスキー分解手段22bによってコレスキー分解された固有値に関して固有ベクトルを算出する処理を実行してもよい。
【0235】
コレスキー分解部31が複数のコレスキー分解手段31a、31bを備えて、対称3重対角行列Tに関してコレスキー分解する処理を並列実行してもよいことは、コレスキー分解部22と同様である。
【0236】
ベクトル算出部32が複数のベクトル算出手段32a、32bを備えて、固有ベクトルを算出する処理を並列実行してもよいことは、ベクトル算出部23と同様である。
【0237】
なお、固有値算出部17の複数の固有値算出手段17a、17bは、図6で示される分割された2個の行列T,Tの固有値と行列要素とから、分割元の行列Tの固有値と行列要素とを算出する処理を、並列実行してもよい。以下、その処理について説明する。
【0238】
まず、式1から、各固有値を求める処理については、並列実行可能なことがわかる。例えば、固有値算出手段17aが半分の固有値を算出し、固有値算出手段17bが残りの半分の固有値を算出してもよい。なお、固有値算出手段17a、17bのそれぞれは、行列D12と、ベクトルzとの各成分の値を有するものとする。
【0239】
また、式2から、各固有ベクトルqを求める処理についても、並列実行可能なことがわかる。例えば、固有値算出手段17aは、算出した固有値に対応する半分の固有ベクトルqを算出し、固有値算出手段17bは、算出した固有値に対応する残りの半分の固有ベクトルqを算出してもよい。
【0240】
また、式3,式4から、f、lを求める処理についても、並列実行可能なことがわかる。例えば、固有値算出手段17aは、算出した固有ベクトルqを用いて、f、lの一部の要素を算出し、固有値算出手段17bは、算出した固有ベクトルqを用いて、f、lの残りの要素を算出してもよい。その後、固有値算出手段17a、あるいは固有値算出手段17bが、固有値算出手段17a、17bの算出したf、lの各要素を統合することによって、式3,式4で示されるf、lを算出することができる。
【0241】
ここでは、zの成分に0が含まれず、D12の対角成分に同じ値が含まれない場合における固有値算出手段17a、17bの処理について説明したが、zの成分に0が含まれる場合、D12の対角成分に同じ値が含まれる場合であっても、固有方程式を解く処理、固有ベクトルqを算出する処理、固有ベクトルから分割元の行列の行列要素を算出する処理を同様に並列実行することができうる。
【0242】
上記のような各構成要素における並列処理において、各手段が固有値等の情報を格納するメモリを用いる場合に、複数の手段が同一のメモリ、すなわち共有メモリを使用してもよく、あるいは、各手段がそれぞれ別のメモリを使用してもよい。
【0243】
なお、図24~図26を用いて、固有値分解部15や固有値算出部17等が2個の手段によって並列処理を実行する場合について説明したが、固有値分解部15や固有値算出部17等が3以上の手段を備え、並列処理を実行してもよい。また、固有値分解部15、固有値算出部17、コレスキー分解部22、ベクトル算出部23,コレスキー分解部31、ベクトル算出部32のそれぞれにおいて並列処理が実行される場合について説明したが、いずれかの任意の1以上の部において、並列処理が実行されなくてもよい。例えば、固有値分解部15において並列処理が行われなくてもよい。
【0244】
また、その並列処理は、1の装置において、2以上のCPU等を用いて実行されてもよく、2以上の装置において実行されてもよい。例えば、図27で示されるように、装置Aと、装置Bとがそれぞれ固有値分解手段15a、15bを備え、各装置において、固有値分解の処理が並列実行されてもよい。この場合には、装置Aの固有値分解手段15aを備える固有値分解部15-1と、装置Bの固有値分解手段15bを備える固有値分解部15-2とによって固有値分解部15が構成されることになる。したがって、固有値分解装置1は、装置Aと、装置Bとからなるシステムを構成することになる。ここでは、固有値分解部15の並列処理について説明したが、その他の固有値算出部17やコレスキー分解部22等についても、2以上の装置による並列処理を行ってもよい。
【0245】
次に、数値実験による性能の評価について説明する。ここでは、本実施の形態による固有値分解装置1の方法(DDC)と、標準的な分割統治法(DC)と、QR法(QR)と、二分法と逆反復法とを組み合わせた方法(B&I)とを比較した。なお、本実施の形態による固有値分解装置1の方法では、コレスキー分解をqd型ツイスト分解法で行い、逆反復(図12のステップS701参照)を1回用いた。分割統治法としては、LAPACKで提供されているDSTEDCを用いた。QR法としては、LAPACKで提供されているDSTEQRを用いた。二分法と逆反復法とを組み合わせた方法としては、LAPACKで提供されているDSTEVRのソースコードを一部変更したものを用いた。実験で用いるLAPACKはBLAS(Basic Linear Algebra Subprograms)を呼び出す。また、本実施の形態による固有値分解装置1の固有値分解部15では、QR法を用いて固有値分解を行った。また、この数値実験では、ランダムに与えた対角行列を、ランチョス法を用いて対称3重対角化した行列を用いた。数値実験結果のグラフにおいて、横軸は、数値実験で用いた行列の大きさである。
【0246】
図28は、計算の実行時間を示す図である。実行時間が短いほど、高速に計算できることになる。図28の数値実験結果のグラフからわかるように、本実施の形態による固有値分解装置1の方法では、他の方法に比べて高速に固有値分解を行うことができることがわかる。
【0247】
図29は、算出された固有ベクトルの直交性を示す図である。値の小さい方が、直交性が高いことになる。図29の数値実験結果のグラフからわかるように、本実施の形態による固有値分解装置1の方法では、他の方法に比べて、固有ベクトルの直交性は高くない。本実施の形態による固有値分解装置1の方法では、二分法と逆反復法とを組み合わせた方法と同程度の直交性を有することがわかる。これは、固有値をまず算出し、その後に、固有ベクトルを算出していることに起因していると考えられる。
【0248】
図30は、全固有値の相対誤差の総和を示す図である。値の小さい方が、誤差が小さく、精度が高いことになる。図30の数値実験結果のグラフからわかるように、本実施の形態による固有値分解装置1の方法では、他の方法に比べて固有値の精度が高いことがわかる。分割統治法により算出した固有値と同程度の固有値の精度を有する。本実施の形態による固有値分解装置1の方法では、分割統治法によって固有値を算出しているため、当然の結果であると考えられる。
【0249】
図31は、前固有ベクトルの誤差の総和を示す図である。値の小さい方が、誤差が小さく、精度が高いことになる。図31の数値実験結果のグラフからわかるように、本実施の形態による固有値分解装置1の方法では、他の方法に比べて固有ベクトルの精度が高いことがわかる。
【0250】
このように、本実施の形態による固有値分解装置1の方法では、高速に固有値分解をすることができ、さらに、固有値、固有ベクトルの精度の高いことがわかる。また、前述のように、並列性にも優れているため、並列処理を行うことによってさらに高速化を図ることも可能である。
【0251】
[主成分分析による顔画像認識]
次に、主成分分析による顔画像認識に固有値分解を応用する場合について説明する。
【0252】
主成分分析による顔画像認識の処理は、以下の各ステップによって行われる。
(1)射影行列の作成
(2)個人の顔画像に関する情報の登録
(3)入力された顔画像と、登録されている個人の顔画像に関する情報とを比較することによる認識
【0253】
まず、上記(1)の射影行列の作成について説明する。
種々の人物の顔画像データを用意する。その顔画像データは、例えば、スキャナやデジタルスチルカメラ、デジタルビデオカメラ等の光学的読み取り機器によって読み取られた2次元画像データであってもよい。各顔画像データは、各画素値を成分とするベクトルで表示することができる。そのベクトルを、次のようにする。ここでは、M個の顔画像データを扱うものとする。
【数72】
JP0005017666B2_000073t.gif

【0254】
次に、その顔画像データのベクトルの平均値μを算出する。
【数73】
JP0005017666B2_000074t.gif

【0255】
その算出したμを用いて、共分散行列Cを次のように算出する。
【数74】
JP0005017666B2_000075t.gif

【0256】
その後、上述の固有値分解装置1によって、その行列Cに対して、固有値分解を行う。すなわち、行列記憶部11に行列Cを蓄積し、その行列Cに対して対角化や、行列の分割等の処理を行うことによって、固有値分解を行う。その固有値分解された結果である固有値λは、固有値記憶部18に蓄積されることになる。また、固有値分解された結果である固有ベクトルvは、固有ベクトル記憶部20で記憶されている固有ベクトルに対して、対角化で用いた直交行列を作用させることによって、前述のようにして求めることができる。なお、添字jは、固有値の降順になるように設定されているものとする。すなわち、固有値λが最大の固有値となる。
【0257】
その求めた固有ベクトルのうち、対応する固有値が上位からd個の固有ベクトルを用いて、射影行列ηを次のように算出する。
【数75】
JP0005017666B2_000076t.gif

【0258】
次に、上記(2)の個人の顔画像に関する情報の登録について説明する。
各個人について、m個の顔画像データ
【数76】
JP0005017666B2_000077t.gif
を用意する。ここで、添字kは、個人を識別するための添字である。そのベクトルに射影行列ηをそれぞれ掛けることにより、次式のようにしてm個のd次元ベクトルξを算出し、その平均値μを算出する。そして、その平均値μを登録しておく。
【数77】
JP0005017666B2_000078t.gif

【0259】
次に、上記(3)の入力された顔画像と、登録されている個人の顔画像に関する情報とを比較することによる認識について説明する。
ある人物の顔画像データのベクトルxが入力されたとする。すると、その顔画像データのベクトルxに射影行列ηを掛けることにより、d次元ベクトルξを算出する。
【数78】
JP0005017666B2_000079t.gif

【0260】
次に、そのd次元ベクトルξと、各個人のd次元ベクトルξの平均値μとの差のノルムを計算する。
【数79】
JP0005017666B2_000080t.gif

【0261】
そのノルムがあらかじめ設定されているしきい値より小さい場合に、入力された顔画像が、添字kで識別される個人の顔画像であると判断してもよい。あるいは、すべての添字kについて上記のノルムを計算し、入力された顔画像が、最小のノルムに対応する添字kで識別される個人の顔画像であると判断してもよい。このようにして、固有値分解装置1を用いて、主成分分析による顔画像認識を行うことができうる。
【0262】
[2次元画像から3次元へ復元する画像処理への応用]
次に、物体の2次元画像から3次元へ復元する画像処理に固有値分解を応用する場合について説明する。
【0263】
複数の2次元画像から3次元復元を行う処理は、以下の各ステップによって行われる。
(1)2次元画像から特徴点を抽出するステップ。
(2)特徴点データより形状(元の物体の特徴点の3次元座標データ)及び回転(3次元データから特徴点データへの変換)に関するデータを計算するステップ。
(3)形状及び回転のデータより可視化を行うステップ。
【0264】
以下、上記(1)、(2)の各ステップについて、図32のフローチャートを用いて説明する。ここで、2次元画像は、例えば、スキャナやデジタルスチルカメラ、デジタルビデオカメラ等の光学的読み取り機器によって読み取られた2次元画像であってもよい。
【0265】
ステップS5001において、2次元画像j(j=1,・・・,m、mは3以上の整数)から特徴点i(i=1,・・・,n、nは2以上の整数)の座標(x,y)を抽出する。取り扱う2次元画像は、弱中心射影画像であることが好ましい。このとき、次の式が成り立つ。
【数80】
JP0005017666B2_000081t.gif

【0266】
ここで、sは物体のスケールに相対するj番目の画像のスケール、r1,j,r2,jはそれぞれ物体座標系に相対するj番目のカメラ座標系の回転行列の1番目と2番目の行ベクトル、(Xi,Yi,Zi)はi番目の点の3次元座標である。物体のスケールは1番目の画像のスケールと同じにし(s=1)、物体の座標系の姿勢は1番目の画像のカメラ座標系と同じにする(r1,1=(1,0,0),r2,1=(0,1,0))。m枚全ての画像におけるn個全ての点の座標を行列Dの要素として並べると、次式のようになる。
【0267】
【数81】
JP0005017666B2_000082t.gif

【0268】
MとSの形から分かるように、Dのランクは3である。ここで、ステップS5001において、行列Dが与えられている。以下、回転に関するデータM及び形状Sを求める。
【0269】
そこで、行列Dの特異値分解
D=UΣV
を考える。ここで、Σは特異値を大小順に対角線上に並べたもので、U及びVはそれぞれ左直交行列、及び右直交行列である。この特異値分解において、前述の固有値分解を用いることができる。すなわち、固有値分解装置1において、行列記憶部11が記憶している対称行列Aを、DDとすることにより、前述の説明のようにして、対称行列Aの固有値分解がなされることになる。なお、固有値分解装置1において、固有ベクトル記憶部20に蓄積されるのは、対称行列DDを変換した対称3重対角行列Tの固有ベクトルであるので、その固有ベクトルを、前述の説明のようにして、行列DDの固有ベクトルに変換する必要がある。また、その行列DDの各固有値の1/2乗(平方根)が各特異値となる。すなわち、特異値σ=(λ1/2となる。また、固有ベクトルと、右直交行列Vの各列を構成する右特異ベクトルxとは等しいため、固有ベクトルを各列に有する行列が、右直交行列Vとなる。なお、特異値σが0でない場合には、左直交行列Uの各列を構成する左特異ベクトルyを次のようにして求めることができる。
【0270】
【数82】
JP0005017666B2_000083t.gif
となる。一方、特異値σが0である場合には、
【数83】
JP0005017666B2_000084t.gif
を解くことにより、yを算出することができる。ここで、次の置き換えを行うことによって、左特異ベクトルyを正規化する。
【数84】
JP0005017666B2_000085t.gif

【0271】
このようにして、左特異ベクトルyを求めることができ、この左特異ベクトルyを用いて、U=(y,y,…,y)とすることにより、左直交行列Uを算出することができる。
【0272】
ここで、画像のデジタル誤差のため、ゼロでない特異値は3つ以上出てくる。しかし、4番目以降の特異値はノイズによるもので、最初の3つの特異値と比べて格段に小さい。そこで、ステップS5002では、最初の3つの特異値に対して特異ベクトルを計算する。すなわち、本実施の形態による固有値分解装置1では、最初の3個の固有値に対して固有ベクトルを計算する。採用する3個のベクトルをまとめると、次式となる。
【0273】
D'=L'Σ'R'=M'S'
ただし、M'=L'(Σ')1/2、S'=(Σ')1/2R'、D'は∥D-D'∥を最小にするランク3の行列である。
【0274】
次に、D'からM及びSを求めたいが、その組合せは唯一ではない。なぜなら、任意の正則行列Cが
D'=(M'C)(C-1S')
を満たすからである。そこで、上式におけるCをM=M'Cを満たすように決める。Cは下記の式を満たす。
【数85】
JP0005017666B2_000086t.gif

【0275】
E=CCとすると、上式からEの6つの要素に関する2m+1個の線形方程式が得られる。m≧3であるので、Eの要素を一意に決めることができる。ステップS5003において、行列Eを求める。
【0276】
次に、ステップS5004において、EからCを求める。Cの自由度(9)はEの自由度(6)より多い。そこで、条件r1j=(1,0,0),r2j=(0,1,0)を加えれば、Cを決めることができる。このとき2つの解(Necker Reversal)が出る。
次に、ステップS5005において、M=M'C及びS=C-1S'より、回転に関するデータM及び形状Sが決まる。
【0277】
[文書検索への応用]
次に、文書検索処理に固有値分解を応用する場合について説明する。文書中からその文書の内容に関連する索引語を抽出し、索引語の重みを計算する処理の後、ベクトル空間モデルでは、この索引語の重みを要素とするベクトルで文書を表現する。ここで、検索対象となる文書をd,d,・・・,dとし、これら文書集合全体を通して全部でm個の索引語w,w,・・・,wがあるとする。このとき、文書dは、次のようなベクトルで表現されることになる。これを文書ベクトルと呼ぶ。
【数86】
JP0005017666B2_000087t.gif

【0278】
ここで、dijは索引語wの文書dにおける重みである。また、文書集合全体は、次のようなm×n行列Dによって表現することができる。
【数87】
JP0005017666B2_000088t.gif

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

【0280】
実際の文書検索においては、与えられた検索質問文と類似した文書を見つけ出す必要があるが、検索質問ベクトルqと各文書ベクトルdの間の類似度を計算することにより行う。ベクトル間の類似度の定義としてはさまざまなものが考えられるが、文書検索においてよく用いられているものはコサイン尺度(2つのベクトルのなす角度)または内積である。
【0281】
【数89】
JP0005017666B2_000090t.gif
なお、ベクトルの長さが1に正規化(コサイン正規化)されている場合には、コサイン尺度と内積とは一致する。
【0282】
図33は、本実施の形態による固有値分解装置1を利用した文書検索方法の一例を示すフローチャートである。
ステップS6001において、質問ベクトルqを受け取る。
【0283】
ここでは、Dの近似行列Dを使った検索を考える。ベクトル空間モデルでは、検索質問ベクトルqと索引語文書行列D中の各文書ベクトルdの間の類似度を計算することにより検索を行うが、ここではDの代わりにDを使う。ベクトル空間モデルでは、文書ベクトルの次元数は索引語の総数と等しい。したがって、検索対象となる文書の数が増えるに従い、文書ベクトルの次元数も増加する傾向にある。しかし、次元数が増加してくると、コンピュータのメモリによる制限や検索時間の増大などの問題が生じてくるばかりでなく、文書中に含まれる不必要な索引語がノイズ的な影響を及ぼし、検索精度を低下させてしまうという現象も起こってくる。潜在的意味インデキシング(latent semantic indexing;LSI)は、高次元の空間にある文書ベクトルを低次元の空間へと射影することにより、検索精度の改善を図る技術である。高次元の空間では別々に扱われていた索引語が、低次元の空間では相互に関連を持ったものとして扱われる可能性もあるため、索引語の持つ意味や概念に基づく検索を行うことができる。たとえば、通常のベクトル空間モデルでは"car"という索引語と"automobile"という索引語はまったく別物であり、一方の索引語による質問ではもう片方の索引語を含んだ文書を検索することができない。しかし、低次元の空間ではこれらの意味的に関連した索引語は1つの次元に縮退することが期待できるため、"car"という検索質問によって"car"を含む文書ばかりでなく"automobile"を含む文書をも検索することが可能となる。潜在的意味インデキシングでは、特異値分解により高次元ベクトルの次元削減を行うが、これは基本的に多変量解析における主成分分析と等価である。
【0284】
ステップS6002において、kを選択する。ここでは、k<rなるkの値を選択する。r=min(n,m)である。kの値は、予め与えられていてもよいし、計算ごとに選択可能であってもよい。次に、ステップS6003では、行列Dの特異値分解を行う。この特異値分解として、前述の2次元画像から3次元へ復元する画像処理への応用で説明したように、固有値分解による方法を用いることができる。すなわち、固有値分解装置1において、行列記憶部11が記憶している対称行列Aを、このたびの行列DDとすることにより、前述の説明のようにして、対称行列Aの固有値分解がなされ、その結果として、行列Dの特異値分解を行うことができうる。なお、この特異値分解では、計算された特異値のうち、大きい順に1番目からk番目までのk個の特異値に対してDの特異ベクトルを算出する。すなわち、固有値分解において、計算された固有値のうち、大きい順に1番目からk番目までのk個の固有値に対して固有ベクトルを算出し、この固有ベクトルを用いて特異ベクトルを算出すればよい。kは、ステップS6002で選択された値である。
【0286】
すなわち、
=UΣV
なるU及びVを計算する。ここで、Uは、最初のk個の左特異ベクトルのみから構成されるm×k行列であり、Vは、最初のk個の右特異ベクトルのみから構成されるn×k行列であり、Σは、最初のk個の特異値のみから構成されるk×k対角行列である。
【0287】
次に、ステップ6004において、行列Dと質問ベクトルqとの類似度を計算する。いま、ベクトルeをn次元の単位ベクトルとすると、Dのj番目の文書ベクトルはDで表すことができる。文書ベクトルDと検索質問ベクトルqとの間の類似度計算は、
【数90】
JP0005017666B2_000091t.gif
としてもよいが、別の定義を用いてもよい。上式では、DをU,Σ,Vから再構成する必要はなく特異値分解の結果から、直接、類似度を計算できることを示している。上式の中に現われるΣは、
Σ=U
と書き直すことができる。この式の右辺は、近似行列Dにおけるj番目の文書ベクトルの基底Uのもとでの座標(文書のk次元表現)を表している。同様に、上式の中のUqは、検索質問ベクトルqの基底Uのもとでの座標(検索質問のk次元表現)である。
【0288】
ステップS6005において、ステップS6004において計算された類似度を基準に、検索結果を出力する。
このように、本実施の形態による固有値分解装置1は、画像処理や検索処理等の種々の処理に対して用いることができうる。なお、上記以外の処理に用いてもよいことは言うまでもない。
【0289】
以上のように、本実施の形態による固有値分解装置1では、分割統治法を用いて固有値のみを算出するため、標準的な分割統治法で実行されるベクトル更新が必要なくなり、標準的な分割統治法よりも非常に高速である。また、固有値から固有ベクトルを算出する処理も、高速に処理することができる。さらに、QR法では固有値を算出する段階の並列化が困難であるが、本実施の形態による固有値分解装置1では、固有値を算出する段階、及び固有値から固有ベクトルを算出する段階のそれぞれにおいて、本質的に高い並列性を持つ。また、本実施の形態による固有値分解装置1では、他の方法と比べて高い精度の得られることがわかる。
【0290】
なお、本実施の形態による固有値分解装置1は、スタンドアロンの装置であってもよく、サーバクライアントシステムにおけるサーバ装置であってもよく、前述のように、複数の装置から構成されるシステムであってもよい。固有値分解装置1がサーバクライアントシステムにおけるサーバ装置である場合には、行列記憶部11で記憶される対称行列Aや、固有値記憶部18で記憶される固有値や固有ベクトル記憶部20で記憶される固有ベクトル等は、インターネットやイントラネット等の通信回線を介して送受信されてもよい。
【0291】
また、本実施の形態による固有値分解装置1における固有ベクトルの算出において、ある固有値を用いた固有ベクトルの算出では、qd型ツイスト分解法を用い、他の固有値を用いた固有ベクトルの算出では、LV型ツイスト分解法を用いるようにしてもよい。例えば、値の近接している固有値については、LV型ツイスト分解法を用いて固有ベクトルの算出を行い、値が近接していない固有値については、qd型ツイスト分解法を用いて固有ベクトルの算出を行ってもよい。固有値の値が近接しているかどうかは、所定の記録媒体等においてあらかじめ設定されているしきい値と比較することによって判断してもよい。
【0292】
また、上記実施の形態において、各処理または各機能は、単一の装置または単一のシステムによって集中処理されることによって実現されてもよく、あるいは、複数の装置または複数のシステムによって分散処理されることによって実現されてもよい。
【0293】
また、上記実施の形態において、各構成要素は専用のハードウェアにより構成されてもよく、あるいは、ソフトウェアにより実現可能な構成要素については、プログラムを実行することによって実現されてもよい。例えば、ハードディスクや半導体メモリ等の記録媒体に記録されたソフトウェア・プログラムをCPU等のプログラム実行部が読み出して実行することによって、各構成要素が実現され得る。なお、上記実施の形態における固有値分解装置を実現するソフトウェアは、以下のようなプログラムである。つまり、このプログラムは、コンピュータに、対称3重対角行列Tが記憶される対角行列記憶部から前記対称3重対角行列Tを読み出し、当該対称3重対角行列Tを2個の対称3重対角行列に分割して前記対角行列記憶部に蓄積し、その対称3重対角行列を2個の対称3重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各対称3重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割ステップと、前記あらかじめ決められた大きさ以下の各対称3重対角行列を前記対角行列記憶部から読み出し、前記各対称3重対角行列に対して固有値分解を行い、前記各対称3重対角行列の固有値と、前記各対称3重対角行列の固有ベクトルからなる直交行列の行列要素とを少なくとも算出し、当該固有値と当該行列要素とを固有値分解記憶部に蓄積する固有値分解ステップと、各対称3重対角行列の固有値と、行列要素とを前記固有値分解記憶部から読み出し、前記固有値と、前記行列要素とから、分割元の対称3重対角行列の固有値と、分割元の対称3重対角行列の行列要素とを算出して前記固有値分解記憶部に蓄積し、その分割元の対称3重対角行列の固有値と、行列要素とを算出する処理を、対称3重対角行列Tの少なくとも1個の固有値を算出するまで繰り返し、前記対称3重対角行列Tの少なくとも1個の固有値を固有値記憶部に蓄積する固有値算出ステップと、前記対角行列記憶部から前記対称3重対角行列Tを読み出し、前記固有値記憶部から前記対称3重対角行列Tの固有値を読み出し、前記対称3重対角行列Tとその固有値とから、ツイスト分解法を用いて前記対称3重対角行列Tの少なくとも1個の固有ベクトルを算出して固有ベクトル記憶部に蓄積する固有ベクトル算出ステップと、を実行させるためのものである。
また、このプログラムでは、前記固有ベクトル算出ステップにおいて、qd型ツイスト分解法により固有ベクトルを算出してもよい。
【0294】
また、このプログラムでは、前記固有ベクトル算出ステップが、前記固有値記憶部から前記対称3重対角行列Tの固有値を読み出し、当該固有値のいずれかが負の値である場合に、前記対角行列記憶部から前記対称3重対角行列Tを読み出し、当該対称3重対角行列Tの固有ベクトルを変化させることなく、当該対称3重対角行列Tのすべての固有値が正の値となるように、当該対称3重対角行列Tを正定値化し、正定値化した後の固有値を前記固有値記憶部に蓄積し、正定値化した後の対称3重対角行列Tを前記対角行列記憶部に蓄積する正定値化ステップと、前記対角行列記憶部からすべての固有値が正の値である対称3重対角行列Tを読み出し、前記固有値記憶部から前記対称3重対角行列Tの正の値の固有値を読み出し、前記対称3重対角行列Tの各要素に関してqd型変換を行うことによって、前記対称3重対角行列Tを上2重対角行列B(+1)及び下2重対角行列B(-1)にコレスキー分解するコレスキー分解ステップと、前記上2重対角行列B(+1)及び下2重対角行列B(-1)の各要素と、前記対称3重対角行列Tの正の値の固有値とを用いて固有ベクトルを算出して前記固有ベクトル記憶部に蓄積するベクトル算出ステップと、を備えてもよい。
また、このプログラムでは、前記固有ベクトル算出ステップにおいて、LV型ツイスト分解法により固有ベクトルを算出してもよい。
【0295】
また、このプログラムでは、前記固有ベクトル算出ステップが、前記対角行列記憶部から前記対称3重対角行列Tを読み出し、前記固有値記憶部から前記対称3重対角行列Tの固有値を読み出し、前記対称3重対角行列Tの各要素に関してミウラ変換、dLVv型変換、逆ミウラ変換を行うことによって、前記対称3重対角行列Tを上2重対角行列B(+1)及び下2重対角行列B(-1)にコレスキー分解するコレスキー分解ステップと、前記上2重対角行列B(+1)及び下2重対角行列B(-1)の各要素と、前記対称3重対角行列Tの固有値とを用いて固有ベクトルを算出して前記固有ベクトル記憶部に蓄積するベクトル算出ステップと、を備えてもよい。
【0296】
また、このプログラムでは、コンピュータに、対称行列Aが記憶される行列記憶部から前記対称行列Aを読み出し、前記対称行列Aを3重対角化した前記対称3重対角行列Tを算出して前記対角行列記憶部に蓄積する対角化ステップをさらに実行させてもよい。
なお、上記プログラムにおいて、情報を蓄積するステップなどでは、ハードウェアでしか行われない処理は少なくとも含まれない。
【0297】
また、このプログラムは、サーバなどからダウンロードされることによって実行されてもよく、所定の記録媒体(例えば、CD-ROMなどの光ディスクや磁気ディスク、半導体メモリなど)に記録されたプログラムが読み出されることによって実行されてもよい。
【0298】
また、このプログラムを実行するコンピュータは、単数であってもよく、複数であってもよい。すなわち、集中処理を行ってもよく、あるいは分散処理を行ってもよい。
【0299】
図34は、上記プログラムを実行して、上記実施の形態による固有値分解装置1を実現するコンピュータの外観の一例を示す模式図である。上記実施の形態は、コンピュータハードウェア及びその上で実行されるコンピュータプログラムによって実現されうる。
【0300】
図34において、コンピュータシステム100は、CD-ROM(Compact Disk Read Only Memory)ドライブ105、FD(Flexible Disk)ドライブ106を含むコンピュータ101と、キーボード102と、マウス103と、モニタ104とを備える。
【0301】
図35は、コンピュータシステムを示す図である。図35において、コンピュータ101は、CD-ROMドライブ105、FDドライブ106に加えて、CPU(Central Processing Unit)111と、ブートアッププログラム等のプログラムを記憶するためのROM(Read Only Memory)112と、CPU111に接続され、アプリケーションプログラムの命令を一時的に記憶すると共に、一時記憶空間を提供するRAM(Random Access Memory)113と、アプリケーションプログラム、システムプログラム、及びデータを記憶するハードディスク114と、CPU111、ROM112等を相互に接続するバス115とを備える。なお、コンピュータ101は、LANへの接続を提供する図示しないネットワークカードを含んでもよい。
【0302】
コンピュータシステム100に、上記実施の形態による固有値分解装置1の機能を実行させるプログラムは、CD-ROM121、またはFD122に記憶されて、CD-ROMドライブ105、またはFDドライブ106に挿入され、ハードディスク114に転送されてもよい。これに代えて、そのプログラムは、図示しないネットワークを介してコンピュータ101に送信され、ハードディスク114に記憶されてもよい。プログラムは実行の際にRAM113にロードされる。なお、プログラムは、CD-ROM121やFD122、またはネットワークから直接、ロードされてもよい。
【0303】
また、行列記憶部11、対角行列記憶部13、固有値分解記憶部16、固有値記憶部18,固有ベクトル記憶部20は、RAM113やハードディスク114によって実現されてもよい。
【0304】
プログラムは、コンピュータ101に、上記実施の形態による固有値分解装置1の機能を実行させるオペレーティングシステム(OS)、またはサードパーティプログラム等を必ずしも含まなくてもよい。プログラムは、制御された態様で適切な機能(モジュール)を呼び出し、所望の結果が得られるようにする命令の部分のみを含んでいてもよい。コンピュータシステム100がどのように動作するのかについては周知であり、詳細な説明を省略する。
また、本発明は、以上の実施の形態に限定されることなく、種々の変更が可能であり、それらも本発明の範囲内に包含されるものであることは言うまでもない。
また、本発明のほんのいくつかの典型的な実施例について上で詳細に説明したが、その典型的な実施例において、発明の利益と新規な技術から実質的にはずれることなく多くの変更が可能であることを当業者は容易に認識することができるであろう。したがって、そのようなすべての変更は、本発明の範囲に含まれるものである。
【産業上の利用可能性】
【0305】
以上のように、本発明による固有値分解装置等によれば、固有値分解を高速に処理することができ、分子軌道法による化学計算処理や統計計算処理、情報検索処理、その他の固有値分解を用いる処理を実行する装置等において有用である。
【図面の簡単な説明】
【0306】
【図1】本発明の実施の形態1による固有値分解装置の構成を示すブロック図
【図2】同実施の形態による固有値分解装置の動作を示すフローチャート
【図3】同実施の形態による固有値分解装置の動作を示すフローチャート
【図4】同実施の形態による固有値分解装置の動作を示すフローチャート
【図5】同実施の形態における対称3重対角行列Tの分割について説明するための図
【図6】同実施の形態における対称3重対角行列Tの固有値の算出について説明するための図
【図7】同実施の形態による固有値分解装置の動作を示すフローチャート
【図8】同実施の形態による固有値分解装置の動作を示すフローチャート
【図9】同実施の形態における対称3重対角行列Tの固有値の算出について説明するための図
【図10】同実施の形態における固有ベクトル算出部の構成を示すブロック図
【図11】同実施の形態による固有値分解装置の動作を示すフローチャート
【図12】同実施の形態による固有値分解装置の動作を示すフローチャート
【図13】同実施の形態における固有ベクトル算出部の構成を示すブロック図
【図14】同実施の形態による固有値分解装置の動作を示すフローチャート
【図15】同実施の形態におけるコレスキー分解について説明するための図
【図16】同実施の形態による固有値分解装置の動作を示すフローチャート
【図17】同実施の形態による固有値分解装置の動作を示すフローチャート
【図18】同実施の形態による固有値分解装置の動作を示すフローチャート
【図19】同実施の形態による固有値分解装置の動作を示すフローチャート
【図20】同実施の形態による固有値分解装置の動作を示すフローチャート
【図21】同実施の形態による固有値分解装置の動作を示すフローチャート
【図22】同実施の形態におけるqd型ツイスト分解法とLV型ツイスト分解法とについて説明するための図
【図23】同実施の形態による固有値分解装置の動作を示すフローチャート
【図24】同実施の形態による固有値分解装置の構成の他の一例を示すブロック図
【図25】同実施の形態における固有ベクトル算出部の構成の他の一例を示すブロック図
【図26】同実施の形態における固有ベクトル算出部の構成の他の一例を示すブロック図
【図27】同実施の形態による固有値分解装置の構成の他の一例を示すブロック図
【図28】同実施の形態による固有値分解装置の数値実験の結果を示す図
【図29】同実施の形態による固有値分解装置の数値実験の結果を示す図
【図30】同実施の形態による固有値分解装置の数値実験の結果を示す図
【図31】同実施の形態による固有値分解装置の数値実験の結果を示す図
【図32】同実施の形態における画像処理の一例を示すフローチャート
【図33】同実施の形態における文書検索処理の一例を示すフローチャート
【図34】同実施の形態におけるコンピュータシステムの外観一例を示す模式図
【図35】同実施の形態におけるコンピュータシステムの構成の一例を示す図
図面
【図1】
0
【図2】
1
【図3】
2
【図4】
3
【図5】
4
【図6】
5
【図7】
6
【図8】
7
【図9】
8
【図10】
9
【図11】
10
【図12】
11
【図13】
12
【図14】
13
【図15】
14
【図16】
15
【図17】
16
【図18】
17
【図19】
18
【図20】
19
【図21】
20
【図22】
21
【図23】
22
【図24】
23
【図25】
24
【図26】
25
【図27】
26
【図28】
27
【図29】
28
【図30】
29
【図31】
30
【図32】
31
【図33】
32
【図34】
33
【図35】
34