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