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

明細書 :特異値分解装置、及び特異値分解方法

発行国 日本国特許庁(JP)
公報種別 特許公報(B2)
特許番号 特許第5011545号 (P5011545)
登録日 平成24年6月15日(2012.6.15)
発行日 平成24年8月29日(2012.8.29)
発明の名称または考案の名称 特異値分解装置、及び特異値分解方法
国際特許分類 G06F  17/16        (2006.01)
G06T   1/00        (2006.01)
FI G06F 17/16 K
G06T 1/00 A
請求項の数または発明の数 23
全頁数 54
出願番号 特願2007-549031 (P2007-549031)
出願日 平成18年9月21日(2006.9.21)
新規性喪失の例外の表示 特許法第30条第1項適用 2005年9月25日発行日本応用数理学会論文誌にて発表
国際出願番号 PCT/JP2006/318713
国際公開番号 WO2007/066445
国際公開日 平成19年6月14日(2007.6.14)
優先権出願番号 2005351089
優先日 平成17年12月5日(2005.12.5)
優先権主張国 日本国(JP)
審査請求日 平成21年8月7日(2009.8.7)
特許権者または実用新案権者 【識別番号】504132272
【氏名又は名称】国立大学法人京都大学
発明者または考案者 【氏名】中村 佳正
【氏名】譽田 太朗
【氏名】岩▲崎▼ 雅史
【氏名】阪野 真也
【氏名】高田 雅美
個別代理人の代理人 【識別番号】100115749、【弁理士】、【氏名又は名称】谷川 英和
審査官 【審査官】田中 幸雄
参考文献・文献 国際公開第2005/119507(WO,A1)
調査した分野 G06F 17/16
G06T 1/00
特許請求の範囲 【請求項1】
2重対角行列Bが記憶される対角行列記憶部と、
前記2重対角行列Bが2個の2重対角行列に分割され、その2重対角行列が2個の2重対角行列に分割される処理が、分割後の各2重対角行列があらかじめ決められた大きさ以下となるまで繰り返され、当該あらかじめ決められた大きさ以下の各2重対角行列に対して特異値分解が行われた結果である、前記各2重対角行列の特異値と、前記各2重対角行列の特異ベクトルからなる左右直交行列の一部の要素である行列要素とが記憶される特異値分解記憶部と、
前記2重対角行列Bの特異値が記憶される特異値記憶部と、
各2重対角行列の特異値と、行列要素とを前記特異値分解記憶部から読み出し、前記特異値と、前記行列要素とから、分割元の2重対角行列の特異値と、分割元の2重対角行列の行列要素とを算出して前記特異値分解記憶部に蓄積し、その分割元の2重対角行列の特異値と、行列要素とを算出する処理を、2重対角行列Bの少なくとも1個の特異値を算出するまで繰り返し、前記2重対角行列Bの少なくとも1個の特異値を前記特異値記憶部に蓄積する特異値算出部と、
前記2重対角行列Bの特異ベクトルが記憶される特異ベクトル記憶部と、
前記対角行列記憶部から2重対角行列Bを読み出し、前記特異値記憶部から前記2重対角行列Bの特異値を読み出し、2重対角行列Bとその特異値とから、ツイスト分解法を用いて2重対角行列Bの少なくとも1個の特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する特異ベクトル算出部と、を備えた特異値分解装置。
【請求項2】
前記対角行列記憶部には、前記あらかじめ決められた大きさ以下の各2重対角行列も記憶され、
前記あらかじめ決められた大きさ以下の各2重対角行列を前記対角行列記憶部から読み出し、前記各2重対角行列に対して特異値分解を行い、前記各2重対角行列の特異値と、前記各2重対角行列の特異ベクトルとを算出し、当該特異値と、当該特異ベクトルからなる左右直交行列の一部の要素である行列要素とを前記特異値分解記憶部に蓄積する特異値分解部をさらに備えた請求項1記載の特異値分解装置。
【請求項3】
2重対角行列Bが記憶される対角行列記憶部と、
前記対角行列記憶部から前記2重対角行列Bを読み出し、当該2重対角行列Bを2個の2重対角行列に分割して前記対角行列記憶部に蓄積し、その2重対角行列を2個の2重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各2重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割部と、
前記あらかじめ決められた大きさ以下の各2重対角行列を前記対角行列記憶部から読み出し、前記各2重対角行列に対して特異値分解を行い、前記各2重対角行列の特異値と、前記各2重対角行列の特異ベクトルとを算出する特異値分解部と、
前記特異値分解部によって特異値分解された特異値と、特異ベクトルからなる左右直交行列の一部の要素である行列要素とが記憶される特異値分解記憶部と、
前記2重対角行列Bの特異値が記憶される特異値記憶部と、
各2重対角行列の特異値と、行列要素とを前記特異値分解記憶部から読み出し、前記特異値と、前記行列要素とから、分割元の2重対角行列の特異値と、分割元の2重対角行列の行列要素とを算出して前記特異値分解記憶部に蓄積し、その分割元の2重対角行列の特異値と、行列要素とを算出する処理を、2重対角行列Bの少なくとも1個の特異値を算出するまで繰り返し、前記2重対角行列Bの少なくとも1個の特異値を前記特異値記憶部に蓄積する特異値算出部と、
前記2重対角行列Bの特異ベクトルが記憶される特異ベクトル記憶部と、
前記対角行列記憶部から2重対角行列Bを読み出し、前記特異値記憶部から前記2重対角行列Bの特異値を読み出し、2重対角行列Bとその特異値とから、ツイスト分解法を用いて2重対角行列Bの少なくとも1個の特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する特異ベクトル算出部と、を備えた特異値分解装置。
【請求項4】
前記特異ベクトル算出部は、qd型ツイスト分解法により特異ベクトルを算出する、請求項3記載の特異値分解装置。
【請求項5】
前記特異ベクトル算出部は、LV型ツイスト分解法により特異ベクトルを算出する、請求項3記載の特異値分解装置。
【請求項6】
前記特異ベクトル算出部は、
前記対角行列記憶部から2重対角行列Bを読み出し、前記特異値記憶部から前記2重対角行列Bの特異値を読み出し、前記2重対角行列Bの各要素に関してミウラ変換、dLVv型変換、逆ミウラ変換を行うことによって、前記2重対角行列Bを上2重対角行列B(+1)及び下2重対角行列B(-1)にコレスキー分解するコレスキー分解部と、
前記上2重対角行列B(+1)及び下2重対角行列B(-1)の各要素と、前記2重対角行列Bの特異値とを用いて一方の左右直交行列を構成する特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する第1特異ベクトル算出部と、
前記第1特異ベクトル算出部が算出した一方の左右直交行列を構成する特異ベクトルと、前記2重対角行列Bの特異値と、前記2重対角行列Bとを用いて他方の左右直交行列を構成する特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する第2特異ベクトル算出部と、をさらに備えた、請求項5記載の特異値分解装置。
【請求項7】
前記コレスキー分解部は、複数のコレスキー分解手段を備え、
前記複数のコレスキー分解手段が、前記2重対角行列Bをコレスキー分解する処理を並列実行する、請求項6記載の特異値分解装置。
【請求項8】
前記第1特異ベクトル算出部は、複数の第1特異ベクトル算出手段を備え、
前記複数の第1特異ベクトル算出手段が、特異ベクトルを算出する処理を並列実行する、請求項6または請求項7記載の特異値分解装置。
【請求項9】
前記第2特異ベクトル算出部は、複数の第2特異ベクトル算出手段を備え、
前記複数の第2特異ベクトル算出手段が、特異ベクトルを算出する処理を並列実行する、請求項6から請求項8のいずれか記載の特異値分解装置。
【請求項10】
前記特異値算出部は、複数の特異値算出手段を備え、
前記複数の特異値算出手段が、分割元の2重対角行列の特異値と、行列要素とを算出する処理を並列実行する、請求項3から請求項9のいずれか記載の特異値分解装置。
【請求項11】
前記特異値分解部は、複数の特異値分解手段を備え、
前記複数の特異値分解手段が、2重対角行列に対して特異値分解を行う処理を並列実行する、請求項3から請求項10のいずれか記載の特異値分解装置。
【請求項12】
行列Aが記憶される行列記憶部と、
前記行列Aを前記行列記憶部から読み出し、前記行列Aを2重対角化した前記2重対角行列Bを算出して前記対角行列記憶部に蓄積する対角化部と、をさらに備えた請求項3から請求項11のいずれか特異値分解装置。
【請求項13】
前記行列分割部は、2重対角行列を略半分の2個の2重対角行列に分割する、請求項3から請求項12のいずれか記載の特異値分解装置。
【請求項14】
前記特異値算出部は、2重対角行列Bの全ての特異値を算出する、請求項3から請求項13のいずれか記載の特異値分解装置。
【請求項15】
前記特異ベクトル算出部は、2重対角行列Bの全ての特異ベクトルを算出する、請求項14記載の特異値分解装置。
【請求項16】
前記行列Aは、2次元画像j(j=1,…,m、mは3以上の整数)から抽出された特徴点i(i=1,…,n、nは2以上の整数)の座標(x,y)から構成される行列である、請求項12記載の特異値分解装置。
【請求項17】
前記行列Aは、索引語の重みを要素とするベクトルであって、検索対象となる文書を示すベクトルであるd,…,d(nは2以上の整数)を各列に有する行列である索引語文書行列である、請求項12記載の特異値分解装置。
【請求項18】
2重対角行列Bが記憶される対角行列記憶部と、前記2重対角行列Bが2個の2重対角行列に分割され、その2重対角行列が2個の2重対角行列に分割される処理が、分割後の各2重対角行列があらかじめ決められた大きさ以下となるまで繰り返され、当該あらかじめ決められた大きさ以下の各2重対角行列に対して特異値分解が行われた結果である、前記各2重対角行列の特異値と、前記各2重対角行列の特異ベクトルからなる左右直交行列の一部の要素である行列要素とが記憶される特異値分解記憶部と、前記2重対角行列Bの特異値が記憶される特異値記憶部と、特異値算出部と、前記2重対角行列Bの特異ベクトルが記憶される特異ベクトル記憶部と、特異ベクトル算出部とを備えた特異値分解装置で用いられる特異値分解方法であって、
前記特異値算出部が、各2重対角行列の特異値と、行列要素とを前記特異値分解記憶部から読み出し、前記特異値と、前記行列要素とから、分割元の2重対角行列の特異値と、分割元の2重対角行列の行列要素とを算出して前記特異値分解記憶部に蓄積し、その分割元の2重対角行列の特異値と、行列要素とを算出する処理を、2重対角行列Bの少なくとも1個の特異値を算出するまで繰り返し、前記2重対角行列Bの少なくとも1個の特異値を前記特異値記憶部に蓄積する特異値算出ステップと、
前記特異ベクトル算出部が、前記対角行列記憶部から2重対角行列Bを読み出し、前記特異値記憶部から前記2重対角行列Bの特異値を読み出し、2重対角行列Bとその特異値とから、ツイスト分解法を用いて2重対角行列Bの少なくとも1個の特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する特異ベクトル算出ステップと、を備えた特異値分解方法。
【請求項19】
2重対角行列Bが記憶される対角行列記憶部と、行列分割部と、特異値分解部と、前記特異値分解部によって特異値分解された特異値と、特異ベクトルからなる左右直交行列の一部の要素である行列要素とが記憶される特異値分解記憶部と、前記2重対角行列Bの特異値が記憶される特異値記憶部と、特異値算出部と、前記2重対角行列Bの特異ベクトルが記憶される特異ベクトル記憶部と、特異ベクトル算出部と、を備えた特異値分解装置で用いられる特異値分解方法であって、
前記行列分割部が、前記対角行列記憶部から前記2重対角行列Bを読み出し、当該2重対角行列Bを2個の2重対角行列に分割して前記対角行列記憶部に蓄積し、その2重対角行列を2個の2重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各2重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割ステップと、
前記特異値分解部が、前記あらかじめ決められた大きさ以下の各2重対角行列を前記対角行列記憶部から読み出し、前記各2重対角行列に対して特異値分解を行い、前記各2重対角行列の特異値と、前記各2重対角行列の特異ベクトルとを算出する特異値分解ステップと、
前記特異値算出部が、各2重対角行列の特異値と、行列要素とを前記特異値分解記憶部から読み出し、前記特異値と、前記行列要素とから、分割元の2重対角行列の特異値と、分割元の2重対角行列の行列要素とを算出して前記特異値分解記憶部に蓄積し、その分割元の2重対角行列の特異値と、行列要素とを算出する処理を、2重対角行列Bの少なくとも1個の特異値を算出するまで繰り返し、前記2重対角行列Bの少なくとも1個の特異値を前記特異値記憶部に蓄積する特異値算出ステップと、
前記特異ベクトル算出部が、前記対角行列記憶部から2重対角行列Bを読み出し、前記特異値記憶部から前記2重対角行列Bの特異値を読み出し、2重対角行列Bとその特異値とから、ツイスト分解法を用いて2重対角行列Bの少なくとも1個の特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する特異ベクトル算出ステップと、を備えた特異値分解方法。
【請求項20】
前記特異ベクトル算出部は、コレスキー分解部と、第1特異ベクトル算出部と、第2特異ベクトル算出部と、をさらに備えており、
前記特異ベクトル算出ステップは、
前記コレスキー分解部が、前記対角行列記憶部から2重対角行列Bを読み出し、前記特異値記憶部から前記2重対角行列Bの特異値を読み出し、前記2重対角行列Bの各要素に関してミウラ変換、dLVv型変換、逆ミウラ変換を行うことによって、前記2重対角行列Bを上2重対角行列B(+1)及び下2重対角行列B(-1)にコレスキー分解するコレスキー分解ステップと、
前記第1特異ベクトル算出部が、前記上2重対角行列B(+1)及び下2重対角行列B(-1)の各要素と、前記2重対角行列Bの特異値とを用いて一方の左右直交行列を構成する特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する第1特異ベクトル算出ステップと、
前記第2特異ベクトル算出部が、前記第1特異ベクトル算出ステップで算出した一方の左右直交行列を構成する特異ベクトルと、前記2重対角行列Bの特異値と、前記2重対角行列Bとを用いて他方の左右直交行列を構成する特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する第2特異ベクトル算出ステップと、をさらに備えた、請求項19記載の特異値分解方法。
【請求項21】
コンピュータに、
2重対角行列Bが2個の2重対角行列に分割され、その2重対角行列が2個の2重対角行列に分割される処理が、分割後の各2重対角行列があらかじめ決められた大きさ以下となるまで繰り返され、当該あらかじめ決められた大きさ以下の各2重対角行列に対して特異値分解が行われた結果である、前記各2重対角行列の特異値と、前記各2重対角行列の特異ベクトルからなる左右直交行列の一部の要素である行列要素とが記憶される特異値分解記憶部から、各2重対角行列の特異値と、行列要素とを読み出し、前記特異値と、前記行列要素とから、分割元の2重対角行列の特異値と、分割元の2重対角行列の行列要素とを算出して前記特異値分解記憶部に蓄積し、その分割元の2重対角行列の特異値と、行列要素とを算出する処理を、2重対角行列Bの少なくとも1個の特異値を算出するまで繰り返し、前記2重対角行列Bの少なくとも1個の特異値を特異値記憶部に蓄積する特異値算出ステップと、
前記2重対角行列Bが記憶される対角行列記憶部から2重対角行列Bを読み出し、前記特異値記憶部から前記2重対角行列Bの特異値を読み出し、2重対角行列Bとその特異値とから、ツイスト分解法を用いて2重対角行列Bの少なくとも1個の特異ベクトルを算出して特異ベクトル記憶部に蓄積する特異ベクトル算出ステップと、を実現させるためのプログラム。
【請求項22】
コンピュータに、
2重対角行列Bが記憶される対角行列記憶部から前記2重対角行列Bを読み出し、当該2重対角行列Bを2個の2重対角行列に分割して前記対角行列記憶部に蓄積し、その2重対角行列を2個の2重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各2重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割ステップと、
前記あらかじめ決められた大きさ以下の各2重対角行列を前記対角行列記憶部から読み出し、前記各2重対角行列に対して特異値分解を行い、前記各2重対角行列の特異値と、前記各2重対角行列の特異ベクトルとを算出し、特異値分解された特異値と、特異ベクトルからなる左右直交行列の一部の要素である行列要素とを特異値分解記憶部に蓄積する特異値分解ステップと、
各2重対角行列の特異値と、行列要素とを前記特異値分解記憶部から読み出し、前記特異値と、前記行列要素とから、分割元の2重対角行列の特異値と、分割元の2重対角行列の行列要素とを算出して前記特異値分解記憶部に蓄積し、その分割元の2重対角行列の特異値と、行列要素とを算出する処理を、2重対角行列Bの少なくとも1個の特異値を算出するまで繰り返し、前記2重対角行列Bの少なくとも1個の特異値を特異値記憶部に蓄積する特異値算出ステップと、
前記対角行列記憶部から2重対角行列Bを読み出し、前記特異値記憶部から前記2重対角行列Bの特異値を読み出し、2重対角行列Bとその特異値とから、ツイスト分解法を用いて2重対角行列Bの少なくとも1個の特異ベクトルを算出して特異ベクトル記憶部に蓄積する特異ベクトル算出ステップと、を実行させるためのプログラム。
【請求項23】
前記特異ベクトル算出ステップは、
前記対角行列記憶部から2重対角行列Bを読み出し、前記特異値記憶部から前記2重対角行列Bの特異値を読み出し、前記2重対角行列Bの各要素に関してミウラ変換、dLVv型変換、逆ミウラ変換を行うことによって、前記2重対角行列Bを上2重対角行列B(+1)及び下2重対角行列B(-1)にコレスキー分解するコレスキー分解ステップと、
前記上2重対角行列B(+1)及び下2重対角行列B(-1)の各要素と、前記2重対角行列Bの特異値とを用いて一方の左右直交行列を構成する特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する第1特異ベクトル算出ステップと、
前記第1特異ベクトル算出ステップで算出した一方の左右直交行列を構成する特異ベクトルと、前記2重対角行列Bの特異値と、前記2重対角行列Bとを用いて他方の左右直交行列を構成する特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する第2特異ベクトル算出ステップと、をさらに備えた、請求項22記載のプログラム。
発明の詳細な説明 【技術分野】
【0001】
本発明は、特異値分解を行う特異値分解装置等に関する。
【背景技術】
【0002】
特異値分解(SVD:Singular Value Decomposition)はデータ処理の中心的な行列演算として、画像処理やデータ検索等の多くの分野に応用されている。
なお、特異値分解の従来例としては、例えば、下記の非特許文献1のものなどが知られている。

【非特許文献1】Ming Gu and Stanley C. Eisenstat、「A Divide-and-Conquer Algorithm for the Bidiagonal SVD」、SIAM Journal on Matrix Analysis and Applications、Vol.16,N0.1,pp.79-92(1995)
【発明の開示】
【発明が解決しようとする課題】
【0003】
近年、これらの応用分野におけるデータ量の増大などに伴って、高速・高精度な特異値分解が求められている。また、特異値分解を並列処理することができる特異値分解装置等の開発が望まれていた。
本発明は、上記状況の下になされたものであり、並列性に優れた高速、高精度な特異値分解を実行可能な特異値分解装置等を提供することを目的とする。
【課題を解決するための手段】
【0004】
上記目的を達成するため、本発明による特異値分解装置は、2重対角行列Bが記憶される対角行列記憶部と、前記2重対角行列Bが2個の2重対角行列に分割され、その2重対角行列が2個の2重対角行列に分割される処理が、分割後の各2重対角行列があらかじめ決められた大きさ以下となるまで繰り返され、当該あらかじめ決められた大きさ以下の各2重対角行列に対して特異値分解が行われた結果である、前記各2重対角行列の特異値と、前記各2重対角行列の特異ベクトルからなる左右直交行列の一部の要素である行列要素とが記憶される特異値分解記憶部と、前記2重対角行列Bの特異値が記憶される特異値記憶部と、各2重対角行列の特異値と、行列要素とを前記特異値分解記憶部から読み出し、前記特異値と、前記行列要素とから、分割元の2重対角行列の特異値と、分割元の2重対角行列の行列要素とを算出して前記特異値分解記憶部に蓄積し、その分割元の2重対角行列の特異値と、行列要素とを算出する処理を、2重対角行列Bの少なくとも1個の特異値を算出するまで繰り返し、前記2重対角行列Bの少なくとも1個の特異値を前記特異値記憶部に蓄積する特異値算出部と、前記2重対角行列Bの特異ベクトルが記憶される特異ベクトル記憶部と、前記対角行列記憶部から2重対角行列Bを読み出し、前記特異値記憶部から前記2重対角行列Bの特異値を読み出し、2重対角行列Bとその特異値とから、ツイスト分解法を用いて2重対角行列Bの少なくとも1個の特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する特異ベクトル算出部と、を備えたものである。
【0005】
このような構成により、まず特異値を算出し、その特異値に基づいてツイスト分解法を用いて特異ベクトルを算出することによって、高速で高精度な特異値分解を実現することができうる。また、並列性にも優れている。
【0006】
また、本発明による特異値分解装置では、前記対角行列記憶部には、前記あらかじめ決められた大きさ以下の各2重対角行列も記憶され、前記あらかじめ決められた大きさ以下の各2重対角行列を前記対角行列記憶部から読み出し、前記各2重対角行列に対して特異値分解を行い、前記各2重対角行列の特異値と、前記各2重対角行列の特異ベクトルとを算出し、当該特異値と、当該特異ベクトルからなる左右直交行列の一部の要素である行列要素とを前記特異値分解記憶部に蓄積する特異値分解部をさらに備えてもよい。
【0007】
このような構成により、特異値を算出する際に、あらかじめ決められた大きさ以下の各2重対角行列に対する特異値分解を、特異値分解装置において行うことができる。
【0008】
また、本発明による特異値分解装置は、2重対角行列Bが記憶される対角行列記憶部と、前記対角行列記憶部から前記2重対角行列Bを読み出し、当該2重対角行列Bを2個の2重対角行列に分割して前記対角行列記憶部に蓄積し、その2重対角行列を2個の2重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各2重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割部と、前記あらかじめ決められた大きさ以下の各2重対角行列を前記対角行列記憶部から読み出し、前記各2重対角行列に対して特異値分解を行い、前記各2重対角行列の特異値と、前記各2重対角行列の特異ベクトルとを算出する特異値分解部と、前記特異値分解部によって特異値分解された特異値と、特異ベクトルからなる左右直交行列の一部の要素である行列要素とが記憶される特異値分解記憶部と、前記2重対角行列Bの特異値が記憶される特異値記憶部と、各2重対角行列の特異値と、行列要素とを前記特異値分解記憶部から読み出し、前記特異値と、前記行列要素とから、分割元の2重対角行列の特異値と、分割元の2重対角行列の行列要素とを算出して前記特異値分解記憶部に蓄積し、その分割元の2重対角行列の特異値と、行列要素とを算出する処理を、2重対角行列Bの少なくとも1個の特異値を算出するまで繰り返し、前記2重対角行列Bの少なくとも1個の特異値を前記特異値記憶部に蓄積する特異値算出部と、前記2重対角行列Bの特異ベクトルが記憶される特異ベクトル記憶部と、前記対角行列記憶部から2重対角行列Bを読み出し、前記特異値記憶部から前記2重対角行列Bの特異値を読み出し、2重対角行列Bとその特異値とから、ツイスト分解法を用いて2重対角行列Bの少なくとも1個の特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する特異ベクトル算出部と、を備えたものである。
【0009】
このような構成により、まず特異値を算出し、その特異値に基づいてツイスト分解法を用いて特異ベクトルを算出することによって、高速で高精度な特異値分解を実現することができうる。また、並列性にも優れている。さらに、全ての特異値や特異ベクトルを算出する必要がない場合には、必要な範囲で特異値や特異ベクトルを算出することができ、処理負荷を軽減することができる。
【0010】
また、本発明による特異値分解装置では、前記特異ベクトル算出部は、qd型ツイスト分解法により特異ベクトルを算出してもよい。
このような構成により、高速で高精度な特異値分解を実現することができうる。また、並列性にも優れている。
【0011】
また、本発明による特異値分解装置では、前記特異ベクトル算出部は、LV型ツイスト分解法により特異ベクトルを算出してもよい。
このような構成により、数値安定的に特異ベクトルの算出を行うことができうる。
【0012】
また、本発明による特異値分解装置では、前記特異ベクトル算出部は、前記対角行列記憶部から2重対角行列Bを読み出し、前記特異値記憶部から前記2重対角行列Bの特異値を読み出し、前記2重対角行列Bの各要素に関してミウラ変換、dLVv型変換、逆ミウラ変換を行うことによって、前記2重対角行列Bを上2重対角行列B(+1)及び下2重対角行列B(-1)にコレスキー分解するコレスキー分解部と、前記上2重対角行列B(+1)及び下2重対角行列B(-1)の各要素と、前記2重対角行列Bの特異値とを用いて一方の左右直交行列を構成する特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する第1特異ベクトル算出部と、前記第1特異ベクトル算出部が算出した一方の左右直交行列を構成する特異ベクトルと、前記2重対角行列Bの特異値と、前記2重対角行列Bとを用いて他方の左右直交行列を構成する特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する第2特異ベクトル算出部と、をさらに備えてもよい。
このような構成により、コレスキー分解の処理において、複数の補助変数を用いることによって、数値安定的に特異ベクトルの算出を行うことができうる。
【0013】
また、本発明による特異値分解装置では、前記コレスキー分解部は、複数のコレスキー分解手段を備え、前記複数のコレスキー分解手段が、前記2重対角行列Bをコレスキー分解する処理を並列実行してもよい。
このような構成により、コレスキー分解の処理を短時間で行うことができうる。
【0014】
また、本発明による特異値分解装置では、前記第1特異ベクトル算出部は、複数の第1特異ベクトル算出手段を備え、前記複数の第1特異ベクトル算出手段が、特異ベクトルを算出する処理を並列実行してもよい。
このような構成により、特異ベクトルを算出する処理を短時間で行うことができうる。
【0015】
また、本発明による特異値分解装置では、前記第2特異ベクトル算出部は、複数の第2特異ベクトル算出手段を備え、前記複数の第2特異ベクトル算出手段が、特異ベクトルを算出する処理を並列実行してもよい。
このような構成により、特異ベクトルを算出する処理を短時間で行うことができうる。
【0016】
また、本発明による特異値分解装置では、前記特異値算出部は、複数の特異値算出手段を備え、前記複数の特異値算出手段が、分割元の2重対角行列の特異値と、行列要素とを算出する処理を並列実行してもよい。
このような構成により、特異値を算出する処理を短時間で行うことができうる。
【0017】
また、本発明による特異値分解装置では、前記特異値分解部は、複数の特異値分解手段を備え、前記複数の特異値分解手段が、2重対角行列に対して特異値分解を行う処理を並列実行してもよい。
このような構成により、特異値分解する処理を短時間で行うことができうる。
【0018】
また、本発明による特異値分解装置では、行列Aが記憶される行列記憶部と、前記行列Aを前記行列記憶部から読み出し、前記行列Aを2重対角化した前記2重対角行列Bを算出して前記対角行列記憶部に蓄積する対角化部と、をさらに備えてもよい。
このような構成により、任意の行列Aについて特異値を算出することができる。また、2重対角行列Bの特異ベクトルを用いることにより、行列Aの特異ベクトルを算出することもできうる。
【0019】
また、本発明による特異値分解装置では、前記行列分割部は、2重対角行列を略半分の2個の2重対角行列に分割してもよい。
このような構成により、並列処理を適切に行うことができうる。
【発明の効果】
【0020】
本発明による特異値分解装置等によれば、高速で高精度な特異値分解を行うことができる。また、並列性にも優れている。
【発明を実施するための最良の形態】
【0021】
以下、本発明による特異値分解装置について、実施の形態を用いて説明する。なお、以下の実施の形態において、同じ符号を付した構成要素及びステップは同一または相当するものであり、再度の説明を省略することがある。
【0022】
(実施の形態1)
本発明の実施の形態1による特異値分解装置について、図面を参照しながら説明する。
図1は、本実施の形態による特異値分解装置1の構成を示すブロック図である。図1において、本実施の形態による特異値分解装置1は、行列記憶部11と、対角化部12と、対角行列記憶部13と、行列分割部14と、特異値分解部15と、特異値分解記憶部16と、特異値算出部17と、特異値記憶部18と、特異ベクトル算出部19と、特異ベクトル記憶部20とを備える。
【0023】
行列記憶部11では、任意の行列Aが記憶される。この行列Aは、各要素が実数である実行列である。なお、行列Aが記憶されているとは、行列Aを示すデータが記憶されている、という意味である。後述する記憶部においても同様である。行列記憶部11は、所定の記録媒体(例えば、半導体メモリや磁気ディスク、光ディスクなど)によって実現されうる。行列記憶部11での記憶は、RAM等における一時的な記憶でもよく、あるいは、長期的な記憶でもよい。行列記憶部11に行列Aが記憶される過程は問わない。例えば、記録媒体を介して行列Aが行列記憶部11で記憶されるようになってもよく、通信回線等を介して送信された行列Aが行列記憶部11で記憶されるようになってもよく、あるいは、キーボードやマウス等の入力デバイスを介して入力された行列Aが行列記憶部11で記憶されるようになってもよい。
【0024】
対角化部12は、行列Aを行列記憶部11から読み出し、その読み出した行列Aを2重対角化した2重対角行列Bを算出する。そして、対角化部12は、その算出した2重対角行列Bを対角行列記憶部13に蓄積する。対角化部12は、例えば、ハウスホルダー(Householder)変換を必要なだけ繰り返し行う方法や、その他の2重対角化法を用いて、行列Aを2重対角化する。ここで、2重対角行列Bは、上2重対角行列であってもよく、下2重対角行列であってもよい。本実施の形態では、2重対角行列Bが上2重対角行列である場合について説明する。
【0025】
対角行列記憶部13では、2重対角行列Bが記憶される。対角行列記憶部13は、所定の記録媒体(例えば、半導体メモリや磁気ディスク、光ディスクなど)によって実現されうる。対角行列記憶部13での記憶は、RAM等における一時的な記憶でもよく、あるいは、長期的な記憶でもよい。
【0026】
行列分割部14は、対角行列記憶部13から2重対角行列Bを読み出し、その2重対角行列Bを2個の2重対角行列に分割して対角行列記憶部13に蓄積する。行列分割部14は、その2重対角行列を2個の2重対角行列に分割して対角行列記憶部13に蓄積する処理を、分割後の各2重対角行列があらかじめ決められた大きさ以下となるまで再帰的に繰り返す。
【0027】
特異値分解部15は、あらかじめ決められた大きさ以下の各2重対角行列を対角行列記憶部13から読み出し、その各2重対角行列に対して特異値分解を行い、各2重対角行列の特異値と、その各2重対角行列の特異ベクトルを算出する。特異値分解部15は、例えば、2分法と逆反復法を組み合わせた方法、MR^3法、QRs法等を用いて特異値分解を行ってもよい。ここで、MR^3法において、例えば、dqdsやpqds等のqds法によって特異値を算出してもよい。また、QRs法による特異値分解を行う場合には、FORTRANにおいて提供されているDBDSQRを用いてもよい。また、qds法によって特異値を算出する場合には、FORTRANにおいて提供されているDLASQを用いてもよい。
これらの特異値分解の方法については、すでに公知であり、その詳細な説明を省略する。特異値分解部15は、算出した特異値を特異値分解記憶部16に蓄積する。また、特異値分解部15は、算出した特異ベクトルからなる左右直交行列の一部の要素である行列要素も特異値分解記憶部16に蓄積する。特異ベクトルからなる左右直交行列とは、特異ベクトルを各列とする左直交行列、及び特異ベクトルを各列とする右直交行列のことである。行列要素の詳細については後述する。
【0028】
特異値分解記憶部16では、特異値分解部によって特異値分解された特異値と、前述の行列要素とが記憶される。また、特異値分解記憶部16では、特異値算出部17によって算出された特異値や行列要素も記憶される。特異値分解記憶部16は、所定の記録媒体(例えば、半導体メモリや磁気ディスク、光ディスクなど)によって実現されうる。特異値分解記憶部16での記憶は、RAM等における一時的な記憶でもよく、あるいは、長期的な記憶でもよい。
【0029】
特異値算出部17は、特異値分解部15によって算出された特異値と、行列要素とを特異値分解記憶部16から読み出し、その特異値と行列要素とから、分割元の2重対角行列の特異値と、分割元の2重対角行列の行列要素とを算出して特異値分解記憶部16に蓄積する。特異値算出部17は、その分割元の2重対角行列の特異値と、行列要素とを算出する処理を、2重対角行列Bの特異値を算出するまで再帰的に繰り返す。そして、特異値算出部17は、算出した2重対角行列Bの特異値を特異値記憶部18に蓄積する。
【0030】
特異値記憶部18では、2重対角行列Bの特異値が記憶される。特異値記憶部18は、所定の記録媒体(例えば、半導体メモリや磁気ディスク、光ディスクなど)によって実現されうる。特異値記憶部18での記憶は、RAM等における一時的な記憶でもよく、あるいは、長期的な記憶でもよい。
【0031】
特異ベクトル算出部19は、対角行列記憶部13から2重対角行列Bを読み出し、特異値記憶部18から2重対角行列Bの特異値を読み出す。そして、特異ベクトル算出部19は、2重対角行列Bとその特異値とから、ツイスト分解法を用いて2重対角行列Bの特異ベクトルを算出して特異ベクトル記憶部20に蓄積する。特異ベクトル算出部19は、それらの処理を実行する、コレスキー分解部21と、第1特異ベクトル算出部22と、第2特異ベクトル算出部23とを備える。特異ベクトル算出部19は、LV型ツイスト分解法により特異ベクトルを算出してもよく、あるいは、qd型ツイスト分解法により特異ベクトルを算出してもよい。本実施の形態では、前者の場合について説明する。
【0032】
コレスキー分解部21は、対角行列記憶部13から2重対角行列Bを読み出し、特異値記憶部18から2重対角行列Bの特異値を読み出す。そして、コレスキー分解部21は、2重対角行列Bの各要素に関してミウラ変換、dLVv型変換、逆ミウラ変換を行うことによって、2重対角行列Bを上2重対角行列B(+1)及び下2重対角行列B(-1)にコレスキー分解する。この処理の詳細については後述する。
【0033】
第1特異ベクトル算出部22は、コレスキー分解部21が算出した上2重対角行列B(+1)及び下2重対角行列B(-1)の各要素と、2重対角行列Bの特異値とを用いて一方の左右直交行列を構成する特異ベクトルを算出して特異ベクトル記憶部20に蓄積する。
【0034】
第2特異ベクトル算出部23は、特異値記憶部18から2重対角行列Bの特異値を読み出す。そして、第2特異ベクトル算出部23は、第1特異ベクトル算出部22が算出した一方の左右直交行列を構成する特異ベクトルと、2重対角行列Bの特異値と、2重対角行列Bとを用いて他方の左右直交行列を構成する特異ベクトルを算出して特異ベクトル記憶部20に蓄積する。このように、第1特異ベクトル算出部22と、第2特異ベクトル算出部23とによって、左右直交行列のそれぞれが算出されることになる。
【0035】
特異ベクトル記憶部20では、2重対角行列Bの特異ベクトルが記憶される。特異ベクトル記憶部20は、所定の記録媒体(例えば、半導体メモリや磁気ディスク、光ディスクなど)によって実現されうる。特異ベクトル記憶部20での記憶は、RAM等における一時的な記憶でもよく、あるいは、長期的な記憶でもよい。
【0036】
なお、行列記憶部11、対角行列記憶部13、特異値分解記憶部16、特異値記憶部18、特異ベクトル記憶部20の任意の2以上の記憶部は、同一の記録媒体によって実現されてもよく、あるいは、別々の記録媒体によって実現されてもよい。前者の場合には、例えば、行列Aの記憶されている領域が行列記憶部11となり、2重対角行列B等の記憶されている領域が対角行列記憶部13となる。
【0037】
また、行列記憶部11、対角行列記憶部13、特異値分解記憶部16、特異値記憶部18、特異ベクトル記憶部20の各記憶部は、2以上の記録媒体から構成されてもよい。
【0038】
次に、本実施の形態による特異値分解装置1の動作について、図2のフローチャートを用いて説明する。
(ステップS101)対角化部12は、行列記憶部11で記憶されている行列Aを読み出し、その行列Aを2重対角化して2重対角行列Bを算出して対角行列記憶部13に蓄積する。
【0039】
(ステップS102)行列分割部14、特異値分解部15、特異値算出部17によって2重対角行列Bの特異値が算出され、特異値記憶部18に蓄積される。この処理の詳細については後述する。
【0040】
(ステップS103)特異ベクトル算出部19は、対角行列記憶部13から2重対角行列Bを読み出し、特異値記憶部18から2重対角行列Bの特異値を読み出し、2重対角行列Bの特異ベクトルを算出して特異ベクトル記憶部20に蓄積する。この処理の詳細については後述する。
【0041】
このようにして、2重対角行列Bの特異値分解が終了する。ここで、行列Aの特異値は2重対角行列Bの特異値に等しいため、行列Aの特異値も算出されたことになる。また、後述するように、所定の変換を行うことによって、行列Aの特異ベクトルも2重対角行列Bの特異ベクトルから容易に算出することができる。
【0042】
次に、図2のフローチャートのステップS102の処理について、図3のフローチャートを用いて説明する。
(ステップS201)行列分割部14は、対角行列記憶部13から2重対角行列Bを読み出し、その2重対角行列Bを2個の2重対角行列に分割して対角行列記憶部13に蓄積する。行列分割部14は、その2重対角行列を2個の2重対角行列に分割する処理を、分割後の各2重対角行列があらかじめ決められた大きさ以下となるまで繰り返す。
【0043】
(ステップS202)特異値分解部15は、対角行列記憶部13で記憶されているあらかじめ決められた大きさ以下の各2重対角行列について、特異値分解を行う。特異値分解部15は、特異値分解の結果である特異値と、特異ベクトルからなる左右直交行列の一部の要素である行列要素とを特異値分解記憶部16に蓄積する。
【0044】
(ステップS203)特異値算出部17は、2重対角行列の特異値と、行列要素とを特異値分解記憶部16から読み出し、分割元の2重対角行列の特異値と、分割元の2重対角行列の行列要素とを算出して特異値分解記憶部16に蓄積する。特異値算出部17は、その分割元の2重対角行列の特異値と、行列要素とを算出する処理を、2重対角行列Bの特異値を算出するまで繰り返し、2重対角行列Bの特異値を特異値記憶部18に蓄積する。このようにして、特異値を算出する処理が終了する。
【0045】
次に、図2のフローチャートのステップS103の処理について、図4のフローチャートを用いて説明する。
(ステップS301)コレスキー分解部21は、対角行列記憶部13から2重対角行列Bを読み出し、特異値記憶部18から2重対角行列Bの特異値を読み出す。そして、コレスキー分解部21は、2重対角行列Bの各要素に関してミウラ変換、dLVv型変換、逆ミウラ変換を行うことによって、2重対角行列Bを上2重対角行列B(+1)及び下2重対角行列B(-1)にコレスキー分解する。
【0046】
(ステップS302)第1特異ベクトル算出部22は、上2重対角行列B(+1)及び下2重対角行列B(-1)の各要素を用いて、一方の左右直交行列を構成する特異ベクトルを算出する。本実施の形態では、第1特異ベクトル算出部22は、右直交行列を構成する右特異ベクトルを算出するものとする。
【0047】
(ステップS303)第1特異ベクトル算出部22は、算出した特異ベクトルを正規化する。すなわち、第1特異ベクトル算出部22は、算出した特異ベクトルのノルムを算出し、特異ベクトルを算出したノルムで割ったものを最終的な特異ベクトルとして特異ベクトル記憶部20に蓄積する。
【0048】
(ステップS304)第2特異ベクトル算出部23は、第1特異ベクトル算出部22が算出した特異ベクトルと、2重対角行列Bの特異値と、2重対角行列Bとを用いて、第1特異ベクトル算出部22が算出した特異ベクトルと異なる方の特異ベクトルを算出する。2重対角行列Bを特異値分解した結果が、第1特異ベクトル算出部22が算出した特異ベクトルからなる直交行列と、第2特異ベクトル算出部23が算出する特異ベクトルからなる直交行列と、特異値記憶部18が記憶している特異値であるため、その性質を用いて、第2特異ベクトル算出部23は、特異ベクトルを算出することができる。
【0049】
(ステップS305)第2特異ベクトル算出部23は、算出した特異ベクトルを正規化する。すなわち、第2特異ベクトル算出部23は、算出した特異ベクトルのノルムを算出し、特異ベクトルを算出したノルムで割ったものを最終的な特異ベクトルとして特異ベクトル記憶部20に蓄積する。このようにして、2重対角行列Bを特異値分解する処理は終了となる。
なお、この後に、行列記憶部11が記憶している行列Aの特異ベクトルを算出する処理が行われてもよいが、ここでは省略している。
【0050】
また、特異ベクトル算出部19が算出した特異ベクトルの精度を上げるために、図5で示されるように、特異ベクトルに関する処理を実行してもよく、あるいは、実行しなくてもよい。すなわち、図示しない逆反復法処理部は、特異ベクトル記憶部20から特異ベクトル算出部19が算出した特異ベクトルを読み出し、その特異ベクトルに対して逆反復法の処理を実行し、その結果の特異ベクトルを特異ベクトル記憶部20に蓄積する(ステップS401)。次に、図示しないグラムシュミット処理部は、高精度が必要かどうか判断する(ステップS402)。この判断は、あらかじめ高精度が必要かどうか設定されている記録媒体等から、その設定を読み出すことによって判断してもよい。そして、図示しないグラムシュミット処理部は、高精度が必要な場合には、特異ベクトル記憶部20から逆反復法の処理の実行された特異ベクトルを読み出し、その特異ベクトルに対してグラムシュミット法の処理を実行し、その結果の特異ベクトルを特異ベクトル記憶部20に蓄積する(ステップS403)。なお、逆反復法や、グラムシュミット法については、すでに公知であり、その詳細な説明を省略する。
【0051】
次に、本実施の形態による特異値分解装置1の動作について、以下、より詳細に説明する。
[行列Aの対角化]
行列Aは、例えば、ハウスホルダー変換等を用いて、以下に示されるように2重対角化することができることが知られている。ここでは、上2重対角化する場合について示すが、下2重対角化も同様にして行うことができる。ここで、2重対角行列Bは、行の数と列の数とが一致する正方行列である。
【数1】
JP0005011545B2_000002t.gif

【0052】
したがって、対角化部12は、上記のようにして、行列記憶部11から行列Aを読み出し、2重対角行列Bを算出することができる。その算出された2重対角行列は、対角行列記憶部13で記憶される(ステップS101)。
【0053】
[2重対角行列Bの分割]
2重対角行列Bを分割する処理について説明する。まず、図6で示されるように、上2重対角行列であって、正方行列であるBの右端に、全ての要素が0である列を加えたものを新たに上2重対角行列Bとする。なお、この新たな2重対角行列Bの特異値は、元の正方行列である2重対角行列Bの特異値と同じである。したがって、今後、この新たなn×(n+1)行列の特異値を求める処理について説明する。
【0054】
図7で示されるように、n×(n+1)の上2重対角行列Bが与えられた場合に、それらを2個の上2重対角行列と、2個の要素(図7では、b,b)とに分けることができる。したがって、上2重対角行列Bは、以下のように分割されることになる。
【数2】
JP0005011545B2_000003t.gif

【0055】
ただし、上2重対角行列Bがn×(n+1)行列であるため、上2重対角行列Bは、(k-1)×k行列であり、上2重対角行列Bは、(n-k)×(n-k+1)行列である。kは、1<k<nとなる整数である。eは、適切な次元におけるj番目の単位ベクトルである。ここで、並列処理を適切に実行するためには、kを、n/2を超えない最大の整数にとるか、あるいは、n/2を下まわらない最小の整数にとることが好適である(このkの値を用いて行列を2個の行列に分割する場合を、「行列を略半分の2個の行列に分割する」と呼ぶことにする)。
【数3】
JP0005011545B2_000004t.gif
本実施の形態では、kを上式のようにとるものとする。なお、kの値は、前述のように、1<k<nの範囲内で任意であることは言うまでもない。したがって、行列分割部14は、上述のようにして、対角行列記憶部13が記憶している上2重対角行列Bを2個の上2重対角行列と、2個の要素に分割することができ、その分割の処理を繰り返すことができる。図8は、図3のフローチャートのステップS201における行列分割部14による行列を分割する処理を示すフローチャートである。
【0057】
(ステップS501)行列分割部14は、カウンタIを「1」に設定する。
(ステップS502)行列分割部14は、対角行列記憶部13からI番目の分割をしていない2重対角行列を読み出し、その2重対角行列を2個の2重対角行列と、2個の要素とに分割する。そして、行列分割部14は、分割した2個の2重対角行列と、2個の要素とを対角行列記憶部13に蓄積する。
【0058】
(ステップS503)行列分割部14は、I番目の分割を行っていない2重対角行列が、対角行列記憶部13で記憶されているかどうか判断する。そして、I番目の分割を行っていない2重対角行列が、対角行列記憶部13で記憶されている場合には、ステップS502に戻り、そうでない場合には、ステップS504に進む。
【0059】
(ステップS504)行列分割部14は、I番目の分割を行った2重対角行列の大きさがあらかじめ決められている大きさ以下かどうか判断する。行列分割部14は、例えば、目的とする行列の大きさ(例えば、25×26など)の記憶されている図示しない記録媒体から、その目的とする行列の大きさを読み出し、対角行列記憶部13で記憶されている、I番目の分割後の2重対角行列がその大きさ以下であるかどうかを判断してもよい。そして、I番目の分割を行った2重対角行列の大きさがあらかじめ決められている大きさ以下である場合には、2重対角行列を分割する処理は終了となり、そうでない場合には、ステップS505に進む。
(ステップS505)行列分割部14は、カウンタIを1だけインクリメントする。そして、ステップS502に戻る。
【0060】
なお、この図8のフローチャートでは、上述のように、行列分割部14が各行列を、略半分の2個の行列に分割するため、I番目の分割後の各2重対角行列の大きさがほとんど同じである場合について説明したが、行列分割部14が各行列を、略半分の2個の行列に分割しない場合には、行列分割部14は、分割後の各行列があらかじめ決められた大きさ以下となるように、分割を繰り返すものとする。
【0061】
また、図8のフローチャートでは、ステップS504において、行列分割部14が、対角行列記憶部13で記憶されている分割後の行列の大きさをあらかじめ決められた大きさと比較する処理を実行する場合について説明したが、これは一例であって、行列分割部14は、ステップS504において、それ以外の処理を行ってもよい。例えば、上述のように、行列分割部14が各行列を、略半分の2個の行列に分割する場合には、元の2重対角行列Bの大きさを知ることができれば、何番目の分割で目的とする行列の大きさとなるのかを知ることができる。したがって、N番目(Nは1以上の整数)の分割で目的とする行列の大きさとなる場合には、ステップS504において、IがNであるかどうかを比較し、NでなければステップS505に進み、Nであれば一連の処理を終了するようにしてもよい。
【0062】
図9は、行列の分割について説明するための図である。まず、行列分割部14は、1番目の分割として、2重対角行列Bを、2重対角行列Bと、2重対角行列Bとに分割する(ステップS501,S502)。2重対角行列Bは、1個しかないため、行列分割部14は、1番目の分割をしていない行列がないと判断する(ステップS503)。また、2重対角行列B等はあらかじめ決められた大きさ以下の行列でないとすると(ステップS504)、行列分割部14は、2番目の分割として、2重対角行列Bを、2重対角行列B11と、2重対角行列B12とに分割する(ステップS505,S502)。この場合には、2番目の分割を行っていない2重対角行列Bが存在するため、行列分割部14は、2重対角行列Bも、2重対角行列B21と、2重対角行列B22とに分割する(ステップS503,S502)。このようにして、分割後の各2重対角行列が目的とする大きさ以下になるまで、行列を分割する処理が繰り返される。なお、図9では、2重対角行列以外の2個の要素については省略している。また、図9において、各2重対角行列は、列の数が行の数よりも1だけ大きい行列である。
【0063】
[分割された行列の特異値分解]
行列Bがn×(n+1)行列であるとすると、行列Bの特異値分解は、次のようになる。
【数4】
JP0005011545B2_000005t.gif

【0064】
ここで、Dは、n×nの対角行列である。(D0)は、行列Dの右側に全ての要素が0の列が1つある行列である。Dの各対角成分はBの特異値である。また、vは、行列Bの特異値分解における右直交行列の一番右側のベクトルである。Vは、行列Bの特異値分解における右直交行列のvをのぞいた行列である。(Vは、全体として(n+1)×(n+1)の行列である。
【0065】
したがって、特異値分解部15は、対角行列記憶部13から分割後の各上2重対角行列を読み出し、上記のようにして、特異値分解を行う(ステップS202)。特異値分解の方法として、例えば、2分法と逆反復法を組み合わせた方法、MR^3法、QRs法等を用いてもよいことは前述の通りである。図9で示されるように上2重対角行列の分割が行われた場合には、特異値分解部15は、各2重対角行列B111,B112,B121,B122,……,B222に対して特異値分解を行うことになる。なお、特異値分解部15は、特異値分解の結果得られた各特異値と、特異値分解の結果得られた左右直交行列の一部の要素である行列要素とを特異値分解記憶部16に蓄積する。ここで、行列要素が、左右直交行列のどの要素を含むのかについては後述する。
【0066】
[特異値の算出]
まず、図10で示されるように、分割された2個の2重対角行列B,Bから、分割元の2重対角行列Bの特異値等を算出する処理について説明する。2重対角行列B,Bは、次のように特異値分解されていたとする。ここで、2重対角行列B,Bは、両者共に列の数が行の数よりも1だけ大きい行列であるとする。
【数5】
JP0005011545B2_000006t.gif

【0067】
この場合に、分割元の2重対角行列Bは、次のようになる。
【数6】
JP0005011545B2_000007t.gif

【0068】
また、上式のb2k-1等は2重対角行列の分割の処理において説明した行列の要素である。ここで、Givens変換を行い、b2kφを0にすると、次のようになる。
【数7】
JP0005011545B2_000008t.gif

【0069】
ただし、
【数8】
JP0005011545B2_000009t.gif
である。したがって、行列Bは、
【数9】
JP0005011545B2_000010t.gif
によって、(M 0)に直交変換される。これらより、行列Bの特異値分解を行うためには、行列Mの特異値分解を行えばよいことになる。ここで、n×nの行列Mを次のようにおく。
【数10】
JP0005011545B2_000011t.gif

【0070】
上記の行列Mにおいて、z、d以外の各要素は0である。ただし、i=1,2,…,n、j=2,3,…,nである。この行列Mの特異値分解は、次の定理によって行うことができる。
(定理1)
【数11】
JP0005011545B2_000012t.gif

【0071】
ここで、計算機によって求めることができるのは、行列Mの真の特異値ではなく、誤差を含む近似値である。したがって、
【数12】
JP0005011545B2_000013t.gif
は、真の値である
【数13】
JP0005011545B2_000014t.gif
と大きく異なる可能性がありうる。すなわち、上記のように計算した場合には、特異ベクトルの計算が数値不安定となることがある。この欠点は、次の定理によって克服できる。
【0072】
(定理2)
【数14】
JP0005011545B2_000015t.gif

【0073】
したがって、定理1を用いて行列Mの近似特異値を算出した後に、定理2を用いて、その近似特異値を真の特異値とする行列Mを再構成する。特異ベクトルは、上記の式5を用いて算出した
【数15】
JP0005011545B2_000016t.gif
と、上記の式2を用いて算出した
【数16】
JP0005011545B2_000017t.gif
とを上記の式3,式4に代入することによって数値安定的に求めることができる。このようにして、次のように行列Mの特異値分解が求められたとする。
【数17】
JP0005011545B2_000018t.gif

【0074】
すると、Bの特異値分解は以下のようになる。
【数18】
JP0005011545B2_000019t.gif

【0075】
このように、対象となる行列を2個の副行列に分割し、その分割した後の各副行列について特異値分解を行う処理を再帰的に行うことによって特異値分解を行う方法は、分割統治法(Divide and Conquer:D&C)と呼ばれている。
【0076】
分割統治法では、特異ベクトルまで算出することになり、その特異ベクトルを算出する処理において行列の計算をしなければならないため、非常に負荷の大きい処理となる。分割統治法において特異値分解をする場合には、例えば、計算時間の95%程度が特異ベクトルを算出するベクトル更新の処理(例えば、上述のように、行列Mの左右直交行列から、行列Bの左右直交行列を算出するために行列を掛け合わせる処理)に費やされることもありうる。
【0077】
一方、特異値のみを算出する場合には、上記処理を簡略化することができる。次に、その方法について説明する。上記の結果より、次のようになる。
【数19】
JP0005011545B2_000020t.gif

【0078】
これより、
【数20】
JP0005011545B2_000021t.gif
となる。ここで、fとφとは、行列Bを特異値分解した結果の右直交行列の最初の行の各要素である。また、lとψとは、行列Bを特異値分解した結果の右直交行列の最後の行の各要素である。この右直交行列の最初の行の各要素と、最後の行の各要素とが行列要素となる。
【0079】
上記の結果から、図10で示されるように、行列B,Bから分割元の行列Bを構成する場合に、行列Bについて、
【数21】
JP0005011545B2_000022t.gif
をそれぞれ算出することができる。このような処理を繰り返すことにより、最終的に、行列Aを2重対角化した2重対角行列Bの特異値を算出することができる。
【0080】
したがって、特異値算出部17は、上述のようにして、特異値分解記憶部16から特異値と行列要素とを読み出し、対角行列記憶部13から行列の分割時に発生した2個の要素(b2k-1、b2k)を読み出し、それらを用いることによって、分割元の2重対角行列の特異値と行列要素とを算出する。そして、それらを算出する処理を繰り返すことによって、最終的に2重対角行列Bの特異値を算出する。図11は、図3のフローチャートのステップS203における特異値算出部17が特異値を算出する処理を示すフローチャートである。
【0081】
(ステップS601)特異値算出部17は、カウンタJを「1」に設定する。
(ステップS602)特異値算出部17は、J番目の特異値の算出は最後の特異値の算出であるかどうか判断する。ここで、最後の特異値の算出とは、2重対角行列Bの特異値を算出することである。そして、最後の特異値の算出である場合には、ステップS606に進み、そうでない場合には、ステップS603に進む。
【0082】
(ステップS603)特異値算出部17は、分割元の2重対角行列の特異値等を算出する。この処理の詳細については後述する。
(ステップS604)特異値算出部17は、J番目の特異値の算出において、分割元の全ての2重対角行列の特異値等を算出したかどうか判断する。そして、J番目の特異値の算出において、分割元の全ての2重対角行列の特異値等を算出した場合には、ステップS605に進み、そうでない場合には、ステップS603に戻る。
【0083】
(ステップS605)特異値算出部17は、カウンタJを1だけインクリメントする。そして、ステップS602に戻る。
(ステップS606)特異値算出部17は、2重対角行列Bの特異値を算出し、その算出した特異値を特異値記憶部18に蓄積する。このようにして、2重対角行列Bの特異値を算出する一連の処理は終了となる。
【0084】
図12は、図11のフローチャートにおけるステップS603の詳細な処理を示すフローチャートである。
(ステップS701)特異値算出部17は、分割元の行列の特異値を式2を用いて算出する。そして、特異値算出部17は、算出した特異値を特異値分解記憶部16に蓄積する。
【0085】
(ステップS702)特異値算出部17は、ステップS701で算出した特異値を用いて、式5のzを算出する。
(ステップS703)特異値算出部17は、ステップS701で算出した特異値と、ステップS702で算出したzとを用いて、式4のvを算出する。このvを算出することにより、vを列ベクトルに有するVを算出したことになる。そして、特異値算出部17は、算出したVを特異値分解記憶部16に蓄積する。
【0086】
(ステップS704)特異値算出部17は、式6から式9を用いて、分割元の行列に関する行列要素を算出する。そして、特異値算出部17は、算出した分割元の行列要素を特異値分解記憶部16に蓄積する。このようにしてステップS603の処理は終了となる。
【0087】
図13は、特異値を算出する処理について説明するための図である。特異値分解記憶部16では、特異値分解部15によって行列B111,B112,...,B222が特異値分解された結果、すなわち、それらの特異値と、行列要素とが記憶されているものとする。まず、特異値算出部17は、特異値等の1番目の算出を開始する(ステップS601,S602)。ここで、特異値等の1番目の算出とは、図13の一番下の行の行列の特異値と、行列要素とから、下から2番目の行の行列の特異値と、行列要素とを算出することである。特異値算出部17は、行列B111と、行列B112との各特異値、及び行列B111と、行列B112との各行列要素を特異値分解記憶部16から読み出す。また、行列B11を行列B111と、行列B112とに分解したときに発生した2個の要素を対角行列記憶部13から読み出す。そして、特異値算出部17は、それらの値を用いて、行列B11の特異値を算出して特異値分解記憶部16に蓄積する(ステップS701)。次に、特異値算出部17は、行列B11の特異値を用いてzの値を算出する(ステップS702)。特異値算出部17は、行列B11の特異値と、zの値とを用いて直交行列Vを算出して特異値分解記憶部16に蓄積する(ステップS703)。最後に、特異値算出部17は、行列B11の行列要素を算出して特異値分解記憶部16に蓄積する(ステップS704)。
【0088】
特異値等の1番目の算出はまだ終わっていないため(ステップS604)、特異値算出部17は、上記説明と同様にして、行列B121、行列B122の特異値、行列要素等に基づいて、行列B12の特異値、行列要素を算出する(ステップS602,S603)。このようにして、特異値等の1番目の算出が終了すると、特異値算出部17は、特異値等の2番目の算出、すなわち、行列B11、行列B12の特異値、行列要素等に基づいて行列Bの特異値、行列要素の算出を行う(ステップS604,S605,S602,S603)。その後、特異値算出部17は、同様にして、行列B21、行列B22の特異値、行列要素等に基づいて行列Bの特異値、行列要素の算出を行う(ステップS604,S603)。
【0089】
特異値算出部17が特異値等の2番目の算出を終了すると(ステップS604)、次の特異値等の算出は、2重対角行列Bの特異値を求める最後の処理となるため(ステップS605,S602)、特異値算出部17は、特異値の算出のみを行い、その算出した特異値を特異値記憶部18に蓄積する(ステップS606)。このようにして、特異値の算出が終了する。
【0090】
[特異値からの特異ベクトルの算出]
まず、行列BBを考える。ここで、行列Bは、前述のn×nの上2重対角行列である。この行列BBを次のように対角化したとする。
【数22】
JP0005011545B2_000023t.gif

【0091】
ここで、一般に次のことが成り立つ。
(1)行列BBは対称な3重対角行列である。
(2)行列BBの固有値は全て正であり、行列Bの特異値σ(σ≧σ≧…≧σ≧0)は、行列BBの固有値と次の関係を有する。
【数23】
JP0005011545B2_000024t.gif
(3)V=Vである。ただし、Vは行列Bの右直交行列である。したがって、行列BBの固有ベクトルxは、行列Bの右特異ベクトルに等しい。
【0092】
したがって、行列BBの固有ベクトルを求めると、行列Bの右直交行列が求まることになる。さらに、行列Bの固有値分解を、
【数24】
JP0005011545B2_000025t.gif
とする。すると、右直交行列Vが求まり、特異値を対角成分に有する行列Σが求まることによって、
【数25】
JP0005011545B2_000026t.gif
から、左直交行列も求まることになり、行列Bが特異値分解される。したがって、行列Bの特異ベクトルを求めることは、行列BB=Tの固有ベクトルを求めることに置き換えることができる。すなわち、次の方程式の固有ベクトルxを求めればよいことになる。
【数26】
JP0005011545B2_000027t.gif

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

【0095】
すると、次式のように書ける。
【数28】
JP0005011545B2_000029t.gif

【0096】
そして、次のようになる。
【数29】
JP0005011545B2_000030t.gif

【0097】
ここで、行列Nをツイスト行列と呼ぶ。また、
【数30】
JP0005011545B2_000031t.gif
であるので、(BB-σI)x=γは、
【数31】
JP0005011545B2_000032t.gif
となる。この簡単な式を解くことによって特異ベクトルを算出することができる。具体的には、
【数32】
JP0005011545B2_000033t.gif
のようにわずかな演算で特異ベクトルを算出することができる。なお、あるρ0について、Dρ0=0あるいはDρ0=0であったとしても、行列Bの(ρ,ρ)成分であるb2ρ-1(0)と、行列Bの(ρ,ρ+1)成分であるb2ρ(0)を用いて、
【数33】
JP0005011545B2_000034t.gif
と特異ベクトルを算出することができる(このように特異ベクトルを算出する処理を例外処理と呼ぶことにする)。なお、残差項のパラメータγ及びkの値は、
【数34】
JP0005011545B2_000035t.gif
の絶対値が最小となるように決定する。したがって、上述のコレスキー分解を求め、ツイスト行列Nを求めることができれば、特異ベクトルを求めることができることになる。そこで、次にコレスキー分解について説明する。
【0098】
図14で示すように、B(0)T(0)-σIをコレスキー分解することは、B(0)に対応する{q(0),e(0)}から、B(±1)に対応する{q(±1),e(±1)}を求めることである。
【0099】
(qd型ツイスト分解法)
まず、従来から知られているqd型ツイスト分解法で用いるqd型変換について説明する(図14参照)。
【数35】
JP0005011545B2_000036t.gif

【0100】
この変換は、さらにstationary qd with shift(stqds)変換
【数36】
JP0005011545B2_000037t.gif
及びreverse-time progressive qd with shift(rpqds)変換
【数37】
JP0005011545B2_000038t.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)=Bの成分{b2k-1(0),b2k(0)}が与えられる、すなわち{q(0),e(0)}が与えられると、σ及びek-1(+1)が一意的に決まるので、この状況は避けることはできない。rpqds変換も同様の性質を持つため、実用的なレベルにまで達したとはいいがたい。LAPACKにおいてFORTRANルーチンDSTEGRとして改良版が公開されているものの欠点は完全に解決されてはいない。
【0101】
(LV型ツイスト分解法)
次に、LV型ツイスト分解法で用いるミウラ変換、dLVv型変換、逆ミウラ変換について説明する(図14参照)。
【数38】
JP0005011545B2_000039t.gif

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

【0103】
次に、dLVv型変換について説明する。この変換は、さらにstationary discrete Lotka-Volterra with variable step-size(stdLVv)変換
【数40】
JP0005011545B2_000041t.gif
及びreverse-time discrete Lotka-Volterra with variable step-size(rdLVv)変換
【数41】
JP0005011545B2_000042t.gif
に分けられる。ただし、
【数42】
JP0005011545B2_000043t.gif
を満たす範囲内でδ(±1)は任意である。
【0104】
最後に、逆ミウラ変換について説明する。この変換は、次のように示される。
【数43】
JP0005011545B2_000044t.gif

【0105】
このようにして、qd型変換と同様に、コレスキー分解を行うことができる。qd型変換では見られない離散Lotka-Volterra型変換の大きな特徴は、任意パラメータを持つことである。すなわち、σ=1/δ(0)-1/δ(±1)を満たす範囲でδ(n)の値を任意に設定できる。δ(n)を変動させると補助変数u(n)の値も変化するが、桁落ちによる誤差や数値不安定が発生するかどうかは事前に判定できる。この判定は、if文によって実装されてもよい。この場合は、δ(n)を再設定後に再度計算すればよい。また、u(±1)が求まれば逆ミウラ変換によってq(±1)及びe(±1)が独立に計算されるので、誤差が伝播しないという性質を持つ。なお、逆ミウラ変換をミウラ変換、ミウラ変換を逆ミウラ変換と呼んでもよく、stdLVv変換をstLVv変換と呼んでもよく、rdLVv変換をrLVv変換と呼んでもよい。
【0106】
ここで、LV型ツイスト分解法による処理のより詳細な処理の一例について説明する。図15~図20は、LV型ツイスト分解法による処理の一例を示すフローチャートである。
【0107】
図15は、コレスキー分解の全体の処理の一例を示す図である。
(ステップS901)コレスキー分解部21は、ミウラ変換を行う。この処理の詳細については後述する。
(ステップS902)コレスキー分解部21は、1/δ(±1)を1/δ(0)-σとする。
【0108】
(ステップS903)コレスキー分解部21は、後述するProcedure1の処理を実行する。
(ステップS904)コレスキー分解部21は、後述するProcedure2の処理を実行する。
【0109】
(ステップS905)コレスキー分解部21は、q(+1),e(+1)がすでに算出されているかどうか判断する。そして、すでに算出されている場合には、コレスキー分解の一連の処理は終了となり、一方、算出されていない場合には、ステップS901に戻る。
【0110】
図16は、図15のフローチャートにおけるステップS903の処理の詳細を示すフローチャートである。
(ステップS1001)コレスキー分解部21は、q(+1),e(+1)がすでに算出されているかどうか判断する。そして、すでに算出されている場合には、終了となり、算出されていない場合には、ステップS1002に進む。
(ステップS1002)コレスキー分解部21は、stdLVv変換等の処理を行う。この処理の詳細については後述する。
【0111】
図17は、図15のフローチャートにおけるステップS904の処理の詳細を示すフローチャートである。
(ステップS1101)コレスキー分解部21は、q(-1),e(-1)がすでに算出されているかどうか判断する。そして、すでに算出されている場合には、終了となり、算出されていない場合には、ステップS1102に進む。
(ステップS1102)コレスキー分解部21は、rdLVv変換等の処理を行う。この処理の詳細については後述する。
【0112】
図18は、図15のフローチャートのステップS901の処理の詳細を示すフローチャートである。
(ステップS1201)コレスキー分解部21は、1/δ(0)の値を決定する。この値は、前述のように任意に決定することができうる。算出された特異値の小さい順にσ,σ,…とする場合に、例えば、コレスキー分解部21は、1/δ(0)の値をσよりも小さい値(例えば、「1」など)に設定し、その後、ステップS1203等で桁落ちの発生する可能性があると判断された場合に、1/δ(0)の値をσとσj+1の間(j=1,2,…)に設定する、というように、jを1ずつ大きくしながら1/δ(0)の値を設定していってもよい。なお、1/δ(0)の値を決定する方法は、これに限定されないことは言うまでもない。
(ステップS1202)コレスキー分解部21は、β1をq(0)-1/δ(0)に設定する。
【0113】
(ステップS1203)コレスキー分解部21は、β1の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップS1204に進み、そうでない場合には、ステップS1201に戻って、1/δ(0)の値を決定する処理を再度実行する。なお、このεも任意に定めることができる。例えば、εとしてマシン・イプシロンを用いてもよい。εの値を大きくすると精度が向上し、εの値を小さくすると精度が低下する。なお、このステップS1203で行っている処理は、桁落ちが発生する可能性について判断する処理である。β1の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。
【0114】
(ステップS1204)コレスキー分解部21は、β2をq(0)/(1+δ(0)(0))に設定する。なお、前述のようにu(0)=0であるため、β2はq(0)に設定されたことになる。
【0115】
(ステップS1205)コレスキー分解部21は、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である。
【0116】
(ステップS1206)コレスキー分解部21は、カウンタkを「1」に設定する。
(ステップS1207)コレスキー分解部21は、α1をe(0)/β1に設定する。前述のように、β1はu2k-1(0)となるから、e(0)/β1にミウラ変換を行うと、α1はδ(0)2k(0)と等しいことになる。
【0117】
(ステップS1208)コレスキー分解部21は、α2をα1+1に設定する。ステップS1207の説明からわかるように、α2は1+δ(0)2k(0)と等しいことになる。
【0118】
(ステップS1209)コレスキー分解部21は、α2の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップS1210に進み、そうでない場合には、ステップS1201に戻って、1/δ(0)の値を決定する処理を再度実行する。なお、このステップS1209で行っている処理は、ステップS1203の処理と同様に、桁落ちが発生する可能性について判断する処理である。α2の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。
【0119】
(ステップS1210)コレスキー分解部21は、α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))に等しくなるからである。
【0120】
(ステップS1211)コレスキー分解部21は、β2をqk+1(0)/α2に設定する。前述のように、α2は1+δ(0)2k(0)と等しいため、β2は、qk+1(0)/(1+δ(0)2k(0))となり、β2におけるkの値を1だけインクリメントしたことになる。
【0121】
(ステップS1212)コレスキー分解部21は、β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だけインクリメントしたことになる。
【0122】
(ステップS1213)コレスキー分解部21は、β1の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップS1214に進み、そうでない場合には、ステップS1201に戻って、1/δ(0)の値を決定する処理を再度実行する。なお、このステップS1213で行っている処理は、ステップS1203の処理と同様に、桁落ちが発生する可能性について判断する処理である。β1の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。
【0123】
(ステップS1214)コレスキー分解部21は、α2×β1を算出してu2k+1(0)(1+δ(0)2k(0))に設定する。前述のように、α2は1+δ(0)2k(0)と等しく、β1は、u2k+1(0)に等しいからである。
【0124】
(ステップS1215)コレスキー分解部21は、カウンタkを1だけインクリメントする。
(ステップS1216)コレスキー分解部21は、カウンタkがmに等しいかどうか判断する。そして、mに等しい場合には、一連の処理は終了となり、そうでない場合には、ステップS1207に戻る。
【0125】
このようにして、u2k+1(0)(1+δ(0)2k(0))と、u2k(0)(1+δ(0)2k-1(0))とが算出されることになる。これらの値は、コレスキー分解部21が有する図示しないメモリ等において一時的に記憶されてもよい。
【0126】
図19は、図16のフローチャートのステップS1002の処理の詳細を示すフローチャートである。
(ステップS1301)コレスキー分解部21は、α2を1+δ(+1)(+1)に設定する。なお、前述のようにu(+1)=0であるため、α2は1に設定されたことになる。
【0127】
(ステップS1302)コレスキー分解部21は、β1をu(0)(1+δ(0)(0))/(1+δ(+1)(+1))に設定する。ここで、u(0)(1+δ(0)(0))の値としては、ステップS1205で算出したものを用いる。なお、前述のようにu(+1)=0である。また、このβ1の式にstdLVv変換を実行すると、β1は、u(+1)に設定されたことになる。
【0128】
(ステップS1303)コレスキー分解部21は、カウンタkを「1」に設定する。
(ステップS1304)コレスキー分解部21は、α1をβ1+1/δ(+1)に設定する。
【0129】
(ステップS1305)コレスキー分解部21は、α1の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップS1306に進み、そうでない場合には、図15のフローチャートのProcedure2(ステップS904)に進む。なお、このステップS1305で行っている処理は、ステップS1203の処理と同様に、桁落ちが発生する可能性について判断する処理である。α1の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。
【0130】
(ステップS1306)コレスキー分解部21は、β2をu2k(0)(1+δ(0)2k-1(0))/α1に設定する。ここで、u2k(0)(1+δ(0)2k-1(0))の値としては、ステップS1210で算出したものを用いる。また、このβ2の式にstdLVv変換を実行すると、β2は、δ(+1)2k(+1)に設定されたことになる。
【0131】
(ステップS1307)コレスキー分解部21は、α1×α2を算出する。前述のように、α1はβ1+1/δ(+1)=u2k-1(+1)+1/δ(+1)に等しく、α2は1+δ(+1)2k-2(+1)に等しいため、α1×α2に逆ミウラ変換を実行するとq(+1)に等しくなる。したがって、コレスキー分解部21は、α1×α2を算出してq(+1)に設定する。
【0132】
(ステップS1308)コレスキー分解部21は、α2を1+β2に設定する。前述のように、β2はδ(+1)2k(+1)と等しいため、α2は1+δ(+1)2k(+1)となる。したがって、α2におけるkの値を1だけインクリメントしたことになる。
【0133】
(ステップS1309)コレスキー分解部21は、α2の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップS1310に進み、そうでない場合には、図15のフローチャートのProcedure2(ステップS904)に進む。なお、このステップS1309で行っている処理は、ステップS1203の処理と同様に、桁落ちが発生する可能性について判断する処理である。α2の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。
【0134】
(ステップS1310)コレスキー分解部21は、β1×β2を算出する。前述のように、β1はu2k-1(+1)に等しく、β2はδ(+1)2k(+1)と等しいため、β1×β2に逆ミウラ変換を実行すると、e(+1)に等しくなる。したがって、コレスキー分解部21は、β1×β2を算出してe(+1)に設定する。
【0135】
(ステップS1311)コレスキー分解部21は、β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だけインクリメントしたことになる。
【0136】
(ステップS1312)コレスキー分解部21は、カウンタkを1だけインクリメントする。
(ステップS1313)コレスキー分解部21は、カウンタkがmに等しいかどうか判断する。そして、mに等しい場合には、ステップS1314に進み、そうでない場合には、ステップS1304に戻る。
【0137】
(ステップS1314)コレスキー分解部21は、α2×(β1+1/δ(+1))を算出する。これは、ステップS1304でα1を更新した後に、ステップS1307でα1×α2を計算することと等しい。したがって、コレスキー分解部21は、α2×(β1+1/δ(+1))を算出してq(+1)に設定する。このようにして、ミウラ変換された結果にstdLVv変換と、逆ミウラ変換とを実行して、q(+1)、e(+1)を算出する処理が終了する。これらの値は、コレスキー分解部21が有する図示しないメモリ等において一時的に記憶されてもよい。
【0138】
図20は、図17のフローチャートのステップS1102の処理の詳細を示すフローチャートである。
(ステップS1401)コレスキー分解部21は、β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)に設定されたことになる。
【0139】
(ステップS1402)コレスキー分解部21は、カウンタkを「m-1」に設定する。
(ステップS1403)コレスキー分解部21は、α1をβ1+1/δ(-1)に設定する。
【0140】
(ステップS1404)コレスキー分解部21は、α1の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップS1405に進み、そうでない場合には、図15のフローチャートのミウラ変換(ステップS901)に戻る。なお、このステップS1404で行っている処理は、ステップS1203の処理と同様に、桁落ちが発生する可能性について判断する処理である。α1の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。
【0141】
(ステップS1405)コレスキー分解部21は、β2をu2k(0)(1+δ(0)2k-1(0))/α1に設定する。ここで、u2k(0)(1+δ(0)2k-1(0))の値としては、ステップS1210で算出したものを用いる。また、このβ2の式にrdLVv変換を実行すると、β2は、δ(-1)2k(-1)に設定されたことになる。
【0142】
(ステップS1406)コレスキー分解部21は、α2を1+β2に設定する。前述のように、β2はδ(-1)2k(-1)と等しいため、α2は1+δ(-1)2k(-1)となる。
【0143】
(ステップS1407)コレスキー分解部21は、α2の絶対値がεよりも大きいかどうか判断する。そして、大きい場合には、ステップS1408に進み、そうでない場合には、図15のフローチャートのミウラ変換(ステップS901)に戻る。なお、このステップS1407で行っている処理は、ステップS1203の処理と同様に、桁落ちが発生する可能性について判断する処理である。α2の絶対値がεよりも大きくない場合には、桁落ちの発生する可能性があると判断されることになる。
【0144】
(ステップS1408)コレスキー分解部21は、α1×α2を算出する。前述のように、α1はβ1+1/δ(-1)=u2k+1(-1)+1/δ(-1)に等しく、α2は1+δ(-1)2k(-1)に等しいため、α1×α2に逆ミウラ変換を実行するとqk+1(-1)に等しくなる。したがって、コレスキー分解部21は、α1×α2を算出してqk+1(-1)に設定する。
【0145】
(ステップS1409)コレスキー分解部21は、β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だけデクリメントしたことになる。
【0146】
(ステップS1410)コレスキー分解部21は、β1×β2を算出する。前述のように、β1はu2k-1(-1)に等しく、β2はδ(-1)2k(-1)と等しいため、β1×β2に逆ミウラ変換を実行すると、e(-1)に等しくなる。したがって、コレスキー分解部21は、β1×β2を算出してe(-1)に設定する。
【0147】
(ステップS1411)コレスキー分解部21は、カウンタkを1だけデクリメントする。
(ステップS1412)コレスキー分解部21は、カウンタkが0に等しいかどうか判断する。そして、0に等しい場合には、ステップS1413に進み、そうでない場合には、ステップS1403に戻る。
【0148】
(ステップS1413)コレスキー分解部21は、β1+1/δ(-1)を算出する。これは、ステップS1403でα1を更新した後に、ステップS1408でα1×α2を計算することと等しい。この場合、α2は1だからである。したがって、コレスキー分解部21は、β1+1/δ(-1)を算出してq(-1)に設定する。このようにして、ミウラ変換された結果にrdLVv変換と、逆ミウラ変換とを実行して、q(-1)、e(-1)を算出する処理が終了する。これらの値は、コレスキー分解部21が有する図示しないメモリ等において一時的に記憶されてもよい。
【0149】
なお、図15のフローチャートの処理によって算出することができるのは1個の特異値に対応する特異ベクトルであるため、全ての特異値に対応する特異ベクトルを算出する場合には、コレスキー分解部21は、図15のフローチャートの処理を特異値の数だけ繰り返すことになる。
【0150】
メモリ消費を抑えるために、補助変数u(n)のための配列は必ずしも用意する必要はない。一方、1+δ(0)(0)のためのメモリ領域を確保し、ミウラ変換、dLVv型変換、逆ミウラ変換のステップにまたがってこの値を利用することで、メモリ消費を抑え、計算量を低減することができる。これにより、誤差も低減される。
【0151】
ここで、誤差について説明する。人間が理想的な状況、すなわち、無限桁の計算をいくらでもできるとすると、qd型ツイスト分解法であっても、LV型ツイスト分解法であっても問題ない。しかしながら、コンピュータで計算を行う場合は注意が必要である。有限桁の計算しか行えないコンピュータ上では、数学的に正しい計算法を使ったとしても必ずしも正しい結果が得られる訳ではない。そればかりか、いつまでたっても計算が終了しないといった思わぬ数値的な問題が発生する場合もある。
【0152】
コンピュータ計算による誤差には、丸め誤差及び桁落ちによる誤差などが知られている。丸め誤差単独では、高々有効桁の最後の桁が真値と比べて異なる程度で大きな誤差にはならない。また、指数部が異なる2つの実数の加算、乗算、除算の少なくとも1つの演算を行えばやはり丸め誤差が生じるが、それ以上の誤差は発生しない。さらに、このような丸め誤差を発生するような操作が繰り返されても、丸めモードがnear(四捨五入)ならば、一方的に切り上げられたり、あるいは切り捨てられたりして誤差が極端に蓄積することは少ない。よって、多くの数値計算法は加算、乗算、除算の少なくとも1つの演算によって発生する丸め誤差を特別注意することは少なく、dLVsルーチンによる特異値計算でも結果的に丸め誤差は一様に増大しない。
【0153】
問題となるのは、同符号の実数の減算及び異符号の実数の加算により生じる、桁落ちである。桁落ちによる誤差で値が0となった後、その値による除算を行うと、0が分母にくるような不定形となり計算不可能となる。こうなるといつまでたっても計算が終了しない。減算→除算と計算する部分がqd型ツイスト分解法、LV型ツイスト分解法の両方に存在するので、桁落ち誤差には十分に注意する必要がある。
【0154】
LV型ツイスト分解法では、上述の桁落ちによる誤差を含んでいるかどうかは減算によって得られた値が小さいかどうかで判断できる。qd型ツイスト分解法の場合、桁落ち誤差を含むことが分かったとしても、それを回避することはできない。なぜならば、初期値として{q(0),e(0)}が与えられると、σは一意的に決定され、{q(±1),e(±1)}も一意的に導出されるためである。すなわち、任意パラメータを持たない自由度のない計算法であるためである。
【0155】
それに対して、LV型ツイスト分解法は、自由に設定できるパラメータδ(0)を持つため、補助変数u(n)の値を様々に変化させることができる(図21A、図21B参照)。すなわち、様々な経路で{q(±1),e(±1)}を計算することができる。よって、桁落ちが発生する場合も回避できる。図18~図20の条件判定によって桁落ちの影響をチェックし、減算によって得られた値の絶対値が小さな数εより大きいという条件が満たされなければ、パラメータδ(0)の設定に戻るというものである。この処理は、条件が満たされるまで繰り返される。なお、精度よりも高速性を重視する場合は、数回条件が満たさなければ(q(+1)=0あるいはq(-1)=0ならば)、例外処理を行ってもよい。
【0156】
qd型ツイスト分解法あるいはLV型ツイスト分解法によってコレスキー分解をすることができると、上述のようにツイストされた行列Nを算出することができ、その行列NのN=eを解くことによって、xを算出することができる。ここで、
【数44】
JP0005011545B2_000045t.gif
とxを置き換えることにより、このxを正規化する。このようにして、右特異ベクトルxを求めることができる。この右特異ベクトルxを用いて、V=(x,x,...,x)とすることにより、右直交行列Vを算出することができる。
【0157】
また、前述のように、
【数45】
JP0005011545B2_000046t.gif
となるため、Vが求まれば、2重対角行列B、特異値を対角成分に有する行列Σは既知であるため、左直交行列を算出することができる。より具体的には、特異値σが0でない場合には、
【数46】
JP0005011545B2_000047t.gif
となる。一方、特異値σが0である場合には、
【数47】
JP0005011545B2_000048t.gif
を解くことにより、yを算出することができる。ここで、xの場合と同様に、
【数48】
JP0005011545B2_000049t.gif
とyを置き換えることにより、このyを正規化する。このようにして、左特異ベクトルyを求めることができる。この左特異ベクトルyを用いて、U=(y,y,...,y)とすることにより、左直交行列Uを算出することができる。このようにして、2重対角行列Bの特異値分解がなされる。なお、この処理の後に、逆反復法や、グラムシュミット法の処理を行ってもよいことは前述の通りである。また、ここでは右直交行列Vを先に算出する場合について説明したが、BBの固有ベクトル、すなわち左直交行列Uを先に算出してもよい。
【0158】
また、2重対角行列Bの特異値分解を行うことができれば、次のように、行列Aの特異値分解も行われる。
【数49】
JP0005011545B2_000050t.gif

【0159】
特異ベクトル算出部19は、上述のようにして、特異値記憶部18が記憶している上2重対角行列Bの特異値と、対角行列記憶部13が記憶している上2重対角行列Bとを用いて、上2重対角行列Bの特異ベクトルを算出する。まず、コレスキー分解部21は、対角行列記憶部13から上2重対角行列Bを読み出し、特異値記憶部18から上2重対角行列Bの特異値を読み出す。そして、コレスキー分解部21は、図22のフローチャートで示されるように各変換を行い、コレスキー分解の処理を行う。ここでは、LV型ツイスト分解法を用いる場合について説明する。
【0160】
コレスキー分解部21は、まず、上2重対角行列Bの各要素の値から、q(0),e(0)を求める。そして、コレスキー分解部21は、ミウラ変換を実行することにより、u(0),u(0)等を順次求めていく(ステップS801)。
【0161】
次に、コレスキー分解部21は、ミウラ変換で得られたu(0)を用いて、stdLVv変換を実行することにより、u(+1)等を順次求めていく(ステップS802)。また、コレスキー分解部21は、ミウラ変換で得られたu(0)を用いて、rdLVv変換を実行することにより、u2m-1(-1)等を順次求めていく(ステップS803)。
【0162】
最後に、コレスキー分解部21は、stdLVv変換で得られたu(+1)及びrdLVv変換で得られたu(-1)を用いて、逆ミウラ変換を実行することにより、q(±1)、e(±1)等を順次求めていく(ステップS804)。q(±1)及びe(±1)が求まると、すなわち、上2重対角行列B(+1)と、下2重対角行列B(-1)が求まると、コレスキー分解部21は、それらを第1特異ベクトル算出部22に渡す。なお、コレスキー分解部21は、各特異値について、それぞれコレスキー分解を行うものとする。
【0163】
なお、ここでは、説明の便宜上、図22のフローチャートを用いてコレスキー分解の処理を説明したが、図15~図20のフローチャートで示されるように処理を行ってもよいことは言うまでもない。
【0164】
第1特異ベクトル算出部22は、特異値記憶部18から特異値を読み出し、式10を用いてkの値を決定する。第1特異ベクトル参集部22は、そのkの値を用いて、q(±1)及びe(±1)からツイスト行列Nを求め、N=eを解くことによって、xを算出する(ステップS302)。
【0165】
次に、第1特異ベクトル算出部22は、算出したxを正規化して、右特異ベクトルxを求めて特異ベクトル記憶部20に蓄積する(ステップS303)。そして、第1特異ベクトル算出部22は、算出した右特異ベクトルxを第2特異ベクトル算出部23に渡す。なお、第1特異ベクトル算出部22は、各特異値について、xを算出する処理と、正規化する処理とを行うことによって、全ての右特異ベクトルを算出することができる。
【0166】
第2特異ベクトル算出部23は、上2重対角行列Bを対角行列記憶部13から読み出し、その上2重対角行列Bの特異値を特異値記憶部18から読み出す。そして、第2特異ベクトル算出部23は、第1特異ベクトル算出部22から受け取った右特異ベクトルxと、上2重対角行列Bと、特異値とを用いて、y=Bx/σ、あるいは、B=0を解くことによって、yを算出する(ステップS304)。次に、第2特異ベクトル算出部23は、右特異ベクトルxの場合と同様に、算出したyを正規化して、左特異ベクトルyを求めて特異ベクトル記憶部20に蓄積する(ステップS305)。このようにして、特異ベクトルの算出が終了する。なお、第2特異ベクトル算出部23は、各特異値について、yを算出する処理と、正規化する処理とを行うことによって、全ての左特異ベクトルを算出することができる。
この後、行列記憶部11が記憶している行列Aの特異ベクトルを算出してもよい。その算出の方法は、前述の通りである。
【0167】
また、特異値算出部17が算出した特異値、特異ベクトル算出部19が算出した特異ベクトルを出力する図示しない出力部を特異値分解装置1が備えてもよい。ここで、図示しない出力部による出力は、例えば、特異値等の表示デバイス(例えば、CRTや液晶ディスプレイなど)への表示でもよく、特異値等の所定の機器への通信回線を介した送信でもよく、特異値等のプリンタによる印刷でもよく、特異値等の記録媒体への蓄積でもよい。なお、その出力部は、出力を行うデバイス(例えば、表示デバイスやプリンタなど)を含んでもよく、あるいは含まなくてもよい。また、その出力部は、ハードウェアによって実現されてもよく、あるいは、それらのデバイスを駆動するドライバ等のソフトウェアによって実現されてもよい。
【0168】
なお、ここでは、コレスキー分解部21がLV型ツイスト分解法によってコレスキー分解を行う場合について説明したが、コレスキー分解部21は、qd型ツイスト分解法によってコレスキー分解を行ってもよい。例えば、コレスキー分解部21は、算出された特異値の分布を見て、分布が密でない場合には、qd型ツイスト分解法を用いるようにしてもよい。
【0169】
また、上記説明では、行列Bが上2重対角行列である場合について説明したが、行列Bは下2重対角行列であっても、その行列Bを上2重対角行列に変換することによって特異値分解を同様に実行することができうる。例えば、行列Cが下2重対角行列であるとすると、上2重対角行列B=Cとすることができる。その上2重対角行列Bを、
B=UΣV
と特異値分解できたとすると、
C=B=(UΣV=VΣU
となる。したがって、下2重対角行列の特異値分解において、特異値は上2重対角行列Bの特異値と同じであり、特異ベクトルは、左右が逆になることがわかる。
【0170】
一方、下2重対角行列Cについても、行列分割部14等が下2重対角行列Cの分割を行い、その分割された下2重対角行列について特異値分解部15が特異値分解を行い、その特異値分解の結果を用いて特異値算出部17が下2重対角行列Cの特異値を算出するようにしてもよい。また、特異ベクトル算出部19において、下2重対角行列Cと、その特異値とを用いてコレスキー分解を行い、その特異値に対応する各特異ベクトルを算出するようにしてもよい。
【0171】
また、上記説明では、特異値算出部17が全ての特異値を算出する場合について説明したが、一部の特異値のみを算出するようにしてもよい。例えば、図13において、特異値算出部17は、行列B,Bまでは、全ての特異値を算出する必要があるが、行列B,Bから上2重対角行列Bの特異値を算出する際に、必要な範囲で特異値を算出してもよい。したがって、特異値算出部17は、2重対角行列Bの少なくとも1個の特異値を算出するものであってもよい。その場合には、不必要な特異値を算出することに伴う余分な処理を実行しなくてよいため、処理負荷を軽減することができうる。このように、一部の特異値のみを算出する場合には、例えば、上2重対角行列Bの特異値を算出する処理に費やす計算時間を約(1+k/m)/2倍にすることができうる。ここで、mは行列サイズであり、kは求める特異値の個数である。
【0172】
また、上記説明では、特異ベクトル算出部19が算出された特異値に対応する全ての特異ベクトルを算出する場合について説明したが、上記説明から明らかなように、特異ベクトルを算出する処理は、特異値ごとに行うことができる。したがって、特異ベクトル算出部19は、2重対角行列Bの少なくとも1個の特異ベクトルを算出するものであってもよい。このように、特異ベクトル算出部19は、必要な範囲で特異ベクトルを算出することができ、不必要な特異ベクトルを算出することに伴う余分な処理を実行しなくてよいため、処理負荷を軽減することができうる。
【0173】
次に、本実施の形態による特異値分解装置1における並列処理について説明する。
本実施の形態による特異値分解装置1では、特異値の計算及び特異ベクトルの計算において、並列的に処理を実行することができる。例えば、図23で示されるように、特異値分解装置1において、特異値分解部15、特異値算出部17、コレスキー分解部21、第1特異ベクトル算出部22、第2特異ベクトル算出部23において、それぞれ並列処理を行ってもよい。
【0174】
特異値分解部15は、複数の特異値分解手段15a、15bを備え、その複数の特異値分解手段15a、15bが、分割後の2重対角行列の特異値と、行列要素とを算出する処理を並列実行してもよい。例えば、図13において、特異値分解手段15aが行列B111,B112,B121,B121の特異値分解を実行し、特異値分解手段15bが行列B211,B212,B221,B222の特異値分解を実行してもよい。
【0175】
特異値算出部17は、複数の特異値算出手段17a、17bを備え、その複数の特異値算出手段17a、17bが、分割元の2重対角行列の特異値と、行列要素とを算出する処理を並列実行してもよい。例えば、図13において、特異値算出手段17aが行列B11,B12,B,Bの特異値等を算出する処理を実行し、特異値算出手段17bが行列B21,B22,Bの特異値等を算出する処理を実行してもよい。
【0176】
コレスキー分解部21は、複数のコレスキー分解手段21a、21bを備え、その複数のコレスキー分解手段21a、21bが、2重対角行列Bをコレスキー分解する処理を並列実行してもよい。例えば、コレスキー分解手段21aが半分の特異値についてコレスキー分解する処理を実行し、コレスキー分解手段21bが残りの半分の特異値についてコレスキー分解する処理を実行してもよい。
【0177】
第1特異ベクトル算出部22は、複数の第1特異ベクトル算出手段22a、22bを備え、その複数の第1特異ベクトル算出手段22a、22bが、特異ベクトルを算出する処理を並列実行してもよい。例えば、第1特異ベクトル算出手段22aがコレスキー分解手段21aによってコレスキー分解された特異値に関して特異ベクトルを算出する処理を実行し、第1特異ベクトル算出手段22bがコレスキー分解手段21bによってコレスキー分解された特異値に関して特異ベクトルを算出する処理を実行してもよい。
【0178】
第2特異ベクトル算出部23は、複数の第2特異ベクトル算出手段23a、23bを備え、その複数の第2特異ベクトル算出手段23a、23bが、特異ベクトルを算出する処理を並列実行してもよい。例えば、第2特異ベクトル算出手段23aが第1特異ベクトル算出手段22aによって算出された特異ベクトルに対応する特異ベクトルを算出する処理を実行し、第2特異ベクトル算出手段23bが第1特異ベクトル算出手段22bによって算出された特異ベクトルに対応する特異ベクトルを算出する処理を実行してもよい。
【0179】
なお、特異値算出部17の複数の特異値算出手段17a、17bは、図10で示される分割された2個の行列B,Bの特異値と行列要素とから、分割元の行列Bの特異値と行列要素とを算出する処理を、並列実行してもよい。以下、その処理について説明する。
【0180】
まず、式2から、各特異値を求める処理については、並列実行可能なことがわかる。ここで、式5を次のように書き換えることにする。
【数50】
JP0005011545B2_000051t.gif

【0181】
この式11から、式5を計算するのに、一部の特異値に対応する式11の部分を計算しておき、他の特異値に対応する式11の部分を後から掛け合わせることによって、最終的に
【数51】
JP0005011545B2_000052t.gif
を算出できることがわかる。したがって、特異値算出手段17a、17bのそれぞれは、まず、2個の行列B,Bの特異値と行列要素とを特異値分解記憶部16から読み出し、2個の行列B,Bの分割時に発生した2個の要素を対角行列記憶部13から読み出す。その読み出された行列要素と、2個の要素とを用いて、zの値を知ることができる。その後、各特異値算出手段17a、17bは、式2を用いて、それぞれが担当する特異値を算出する。この処理は、並列実行することができる。
【0182】
次に、特異値算出手段17aは、算出した特異値を用いて計算することができる式11の部分について計算する。また同様に、特異値算出手段17bも、算出した特異値を用いて計算することができる式11の部分について計算する。そして、特異値算出手段17a、17bは、その計算した値を交換し、あらかじめ計算していた値と掛け合わせることによって、最終的に式11の値を計算することができる。このように、特異値算出手段17a、17bは、式11を計算する処理も並列実行することができる。
【0183】
次に、特異値算出手段17a、17bは、式4を用いて、それぞれが算出した特異値に対応する特異ベクトルを計算する。このように、特異値算出手段17a、17bは、右直交行列Vを計算する処理も並列実行することができる。
【0184】
その後、特異値算出手段17a、あるいは特異値算出手段17bが、式6から式9を用いて行列Bの行列要素を計算することによって、分割元の行列Bの特異値と行列要素とを算出する処理が終了する。このように、分割された2個の行列B,Bの特異値と行列要素とから、分割元の行列Bの特異値と行列要素とを算出する処理を、特異値算出手段17a、17bによって並列実行することができる。
【0185】
上記のような各構成要素における並列処理において、各手段が特異値等の情報を格納するメモリを用いる場合に、複数の手段が同一のメモリ、すなわち共有メモリを使用してもよく、あるいは、各手段がそれぞれ別のメモリを使用してもよい。
【0186】
なお、図23を用いて、特異値分解部15や特異値算出部17等が2個の手段によって並列処理を実行する場合について説明したが、特異値分解部15や特異値算出部17等が3以上の手段を備え、並列処理を実行してもよい。また、特異値分解部15、特異値算出部17、コレスキー分解部21、第1特異ベクトル算出部22、第2特異ベクトル算出部23のそれぞれにおいて並列処理が実行される場合について説明したが、いずれかの任意の1以上の部において、並列処理が実行されなくてもよい。例えば、特異値分解部15において並列処理が行われなくてもよい。
【0187】
また、その並列処理は、1の装置において、2以上のCPU等を用いて実行されてもよく、2以上の装置において実行されてもよい。例えば、図24で示されるように、装置Aと、装置Bとがそれぞれ特異値分解手段15a、15bを備え、各装置において、特異値分解の処理が並列実行されてもよい。この場合には、装置Aの特異値分解手段15aを備える特異値分解部15-1と、装置Bの特異値分解手段15bを備える特異値分解部15-2とによって特異値分解部15が構成されることになる。したがって、特異値分解装置1は、装置Aと、装置Bとからなるシステムを構成することになる。ここでは、特異値分解部15の並列処理について説明したが、その他の特異値算出部17やコレスキー分解部21等についても、2以上の装置による並列処理を行ってもよい。
【0188】
[2次元画像から3次元へ復元する画像処理への応用]
次に、物体の2次元画像から3次元へ復元する画像処理に特異値分解を応用する場合について説明する。
【0189】
複数の2次元画像から3次元復元を行う処理は、以下の各ステップによって行われる。
(1)2次元画像から特徴点を抽出するステップ。
(2)特徴点データより形状(元の物体の特徴点の3次元座標データ)及び回転(3次元データから特徴点データへの変換)に関するデータを計算するステップ。
(3)形状及び回転のデータより可視化を行うステップ。
【0190】
以下、上記(1)、(2)の各ステップについて、図25のフローチャートを用いて説明する。ここで、2次元画像は、例えば、スキャナやデジタルスチルカメラ、デジタルビデオカメラ等の光学的読み取り機器によって読み取られた2次元画像であってもよい。
【0191】
ステップS5001において、2次元画像j(j=1,・・・,m、mは3以上の整数)から特徴点i(i=1,・・・,n、nは2以上の整数)の座標(x,y)を抽出する。取り扱う2次元画像は、弱中心射影画像であることが好ましい。このとき、次の式が成り立つ。
【数52】
JP0005011545B2_000053t.gif

【0192】
ここで、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の要素として並べると、次式のようになる。
【数53】
JP0005011545B2_000054t.gif

【0193】
MとSの形から分かるように、Dのランクは3である。ここで、ステップ5001において、行列Dが与えられている。以下、回転に関するデータM及び形状Sを求める。
【0194】
そこで、行列Dの特異値分解
D=UΣV
を考える。ここで、Σは特異値を大小順に対角線上に並べたもので、U及びVはそれぞれ左直交行列、及び右直交行列である。この特異値分解として、前述の方法を用いることができる。すなわち、特異値分解装置1において、行列記憶部11が記憶している行列Aを、このたびの行列Dとすることにより、前述の説明のようにして、行列Aの特異値分解がなされることになる。なお、特異値分解装置1において、特異ベクトル記憶部20に蓄積されるのは、行列Dを変換した上2重対角行列Bの特異ベクトルであるので、その特異ベクトルを、前述の説明のようにして、行列Dの特異ベクトルに変換する必要がある。
【0195】
ここで、画像のデジタル誤差のため、ゼロでない特異値は3つ以上出てくる。しかし、4番目以降の特異値はノイズによるもので、最初の3つの特異値と比べて格段に小さい。そこで、ステップS5002では、最初の3つの特異値に対して特異ベクトルを計算する。採用する3個のベクトルをまとめると、次式となる。
D'=L'Σ'R'=M'S'
ただし、M'=L'(Σ')1/2、S'=(Σ')1/2R'、D'は∥D-D'∥を最小にするランク3の行列である。
【0196】
次に、D'からM及びSを求めたいが、その組合せは唯一ではない。なぜなら、任意の正則行列Cが
D'=(M'C)(C-1S')
を満たすからである。そこで、上式におけるCをM=M'Cを満たすように決める。Cは下記の式を満たす。
【数54】
JP0005011545B2_000055t.gif

【0197】
E=CCとすると、上式からEの6つの要素に関する2m+1個の線形方程式が得られる。m≧3であるので、Eの要素を一意に決めることができる。ステップS5003において、行列Eを求める。
【0198】
次に、ステップ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が決まる。
【0199】
[文書検索への応用]
次に、文書検索処理に特異値分解を応用する場合について説明する。文書中からその文書の内容に関連する索引語を抽出し、索引語の重みを計算する処理の後、ベクトル空間モデルでは、この索引語の重みを要素とするベクトルで文書を表現する。ここで、検索対象となる文書をd,d,・・・,dとし、これら文書集合全体を通して全部でm個の索引語w,w,・・・,wがあるとする。このとき、文書dは、次のようなベクトルで表現されることになる。これを文書ベクトルと呼ぶ。
【数55】
JP0005011545B2_000056t.gif

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

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

【0202】
実際の文書検索においては、与えられた検索質問文と類似した文書を見つけ出す必要があるが、検索質問ベクトルqと各文書ベクトルdの間の類似度を計算することにより行う。ベクトル間の類似度の定義としてはさまざまなものが考えられるが、文書検索においてよく用いられているものはコサイン尺度(2つのベクトルのなす角度)または内積である。
【数58】
JP0005011545B2_000059t.gif

【0203】
なお、ベクトルの長さが1に正規化(コサイン正規化)されている場合には、コサイン尺度と内積とは一致する。
【0204】
図26は、本実施の形態による特異値分解装置1を利用した文書検索方法の一例を示すフローチャートである。
ステップ6001において、質問ベクトルqを受け取る。
【0205】
ここでは、Dの近似行列Dを使った検索を考える。ベクトル空間モデルでは、検索質問ベクトルqと索引語文書行列D中の各文書ベクトルdの間の類似度を計算することにより検索を行うが、ここではDの代わりにDを使う。ベクトル空間モデルでは、文書ベクトルの次元数は索引語の総数と等しい。したがって、検索対象となる文書の数が増えるに従い、文書ベクトルの次元数も増加する傾向にある。しかし、次元数が増加してくると、コンピュータのメモリによる制限や検索時間の増大などの問題が生じてくるばかりでなく、文書中に含まれる不必要な索引語がノイズ的な影響を及ぼし、検索精度を低下させてしまうという現象も起こってくる。潜在的意味インデキシング(latent semantic indexing;LSI)は、高次元の空間にある文書ベクトルを低次元の空間へと射影することにより、検索精度の改善を図る技術である。高次元の空間では別々に扱われていた索引語が、低次元の空間では相互に関連を持ったものとして扱われる可能性もあるため、索引語の持つ意味や概念に基づく検索を行うことができる。たとえば、通常のベクトル空間モデルでは"car"という索引語と"automobile"という索引語はまったく別物であり、一方の索引語による質問ではもう片方の索引語を含んだ文書を検索することができない。しかし、低次元の空間ではこれらの意味的に関連した索引語は1つの次元に縮退することが期待できるため、"car"という検索質問によって"car"を含む文書ばかりでなく"automobile"を含む文書をも検索することが可能となる。潜在的意味インデキシングでは、特異債分解により高次元ベクトルの次元削減を行うが、これは基本的に多変量解析における主成分分析と等価である。
【0206】
ステップS6002において、kを選択する。ここでは、k<rなるkの値を選択する。r=min(n,m)である。kの値は、予め与えられていてもよいし、計算ごとに選択可能であってもよい。次に、ステップS6003では、行列Dの特異値分解を行う。この特異値分解として、前述の方法を用いることができる。すなわち、特異値分解装置1において、行列記憶部11が記憶している行列Aを、このたびの行列Dとすることにより、前述の説明のようにして、行列Aの特異値分解がなされることになる。なお、特異値分解装置1において、特異ベクトル記憶部20に蓄積されるのは、行列Dを変換した上2重対角行列Bの特異ベクトルであるので、その特異ベクトルを、前述の説明のようにして、行列Dの特異ベクトルに変換する必要がある。また、この特異値分解では、計算された特異値のうち、大きい順に1番目からk番目までのk個の特異値に対してDの特異ベクトルを算出する。kは、ステップS6002で選択された値である。
【0208】
すなわち、
=UΣV
なるU及びVを計算する。ここで、Uは、最初のk個の左特異ベクトルのみから構成されるm×k行列であり、Vは、最初のk個の右特異ベクトルのみから構成されるn×k行列であり、Σは、最初のk個の特異値のみから構成されるk×k対角行列である。
【0209】
次に、ステップ6004において、行列Dと質問ベクトルqとの類似度を計算する。いま、ベクトルeをn次元の単位ベクトルとすると、Dのj番目の文書ベクトルはDで表すことができる。文書ベクトルDと検索質問ベクトルqとの間の類似度計算は、
【数59】
JP0005011545B2_000060t.gif
としてもよいが、別の定義を用いてもよい。上式では、DをU,Σ,Vから再構成する必要はなく特異値分解の結果から、直接、類似度を計算できることを示している。上式の中に現われるΣは、
Σ=U
と書き直すことができる。この式の右辺は、近似行列Dにおけるj番目の文書ベクトルの基底Uのもとでの座標(文書のk次元表現)を表している。同様に、上式の中のUqは、検索質問ベクトルqの基底Uのもとでの座標(検索質問のk次元表現)である。
【0210】
ステップS6005において、ステップS6004において計算された類似度を基準に、検索結果を出力する。
最後に、数値実験による性能の評価について説明する。ここでは、本実施の形態による特異値分解装置1の方法と、標準的な分割統治法と、シフト付きQR法とを比較した。分割統治法としては、LAPACKで提供されているDBDSDCを用いた。QR法としては、LAPACKで提供されているDBDSQRを用いた。実験で用いるLAPACKはATLAS(Automatically Tuned Linear Algebra Software)で最適化したBLAS(Basic Linear Algebra Subprograms)を呼び出す。また、本実施の形態による特異値分解装置1の特異値分解部15では、QR法を用いて特異値分解を行った。計算機は、キャッシュメモリL1D:64KB,L1I:64KB,L2:1024KBを搭載したOpteron 1.8GHzを用いた。
【0211】
図27は、計算時間を示す図である。特異値分解装置1による特異値分解が、他の特異値分解方法よりも非常に高速であることがわかる。特に、分割統治法(D&C)と比較すると、ベクトル更新を省略した効果がわかる。
【0212】
図28は、異なる行列サイズの計算時間の変化を示す図である。特異値分解装置1による特異値分解は、ほぼO(n)という計算量を持つことがわかる。
なお、図27及び図28で示される計算時間の測定では、対角成分が2.001,非対角成分が2.0であり、近接特異値を持たない行列について計算を行った。
【0213】
図29は、計算制度を示す図である。特異値分解装置1による特異値の算出は、分割統治法と同程度の高い精度を有することがわかる。図29では、100個のランダム行列を解き、精度を評価した。行列サイズは、1000である。
【0214】
以上のように、本実施の形態による特異値分解装置1では、分割統治法を用いて特異値のみを算出するため、標準的な分割統治法で実行されるベクトル更新が必要なくなり、標準的な分割統治法よりも非常に高速である。また、特異値から特異ベクトルを算出する処理も、高速に処理することができる。さらに、QR法では特異値を算出する段階の並列化が困難であるが、本実施の形態による特異値分解装置1では、特異値を算出する段階、及び特異値から特異ベクトルを算出する段階のそれぞれにおいて、本質的に高い並列性を持つ。また、本実施の形態による特異値分解装置1の精度は、分割統治法やQR法とほぼ同等の精度が得られることがわかる。
【0215】
なお、本実施の形態による特異値分解装置1は、スタンドアロンの装置であってもよく、サーバクライアントシステムにおけるサーバ装置であってもよく、前述のように、複数の装置から構成されるシステムであってもよい。特異値分解装置1がサーバクライアントシステムにおけるサーバ装置である場合には、行列記憶部11で記憶される行列Aや、特異値記憶部18で記憶される特異値や特異ベクトル記憶部20で記憶される特異ベクトル等は、インターネットやイントラネット等の通信回線を介して送受信されてもよい。
【0216】
また、特異値分解装置1における一部の処理を、特異値分解装置1以外において行ってもよい。例えば、行列分解部14による2重対角行列Bを分割する処理、あるいは、特異値分解部15による特異値分解の処理を、特異値分解装置1において実行しなくてもよい。その場合には、例えば、2重対角行列Bが分割された結果である、あらかじめ決められた大きさ以下の各2重対角行列が、入力デバイスや通信回線、記録媒体等を介して受け付けられ、対角行列記憶部13で記憶されていてもよい。また、例えば、そのあらかじめ決められた大きさ以下の各2重対角行列に対して特異値分解が行われた結果である、各2重対角行列の特異値と、各2重対角行列の特異ベクトルからなる左右直交行列の一部の要素である行列要素とが、入力デバイスや通信回線、記録媒体等を介して受け付けられ、特異値分解記憶部15で記憶されていてもよい。
【0217】
また、上記実施の形態において、各処理または各機能は、単一の装置または単一のシステムによって集中処理されることによって実現されてもよく、あるいは、複数の装置または複数のシステムによって分散処理されることによって実現されてもよい。
【0218】
また、上記実施の形態において、各構成要素は専用のハードウェアにより構成されてもよく、あるいは、ソフトウェアにより実現可能な構成要素については、プログラムを実行することによって実現されてもよい。例えば、ハードディスクや半導体メモリ等の記録媒体に記録されたソフトウェア・プログラムをCPU等のプログラム実行部が読み出して実行することによって、各構成要素が実現され得る。なお、上記実施の形態における特異値分解装置を実現するソフトウェアは、以下のようなプログラムである。つまり、このプログラムは、コンピュータに、2重対角行列Bが2個の2重対角行列に分割され、その2重対角行列が2個の2重対角行列に分割される処理が、分割後の各2重対角行列があらかじめ決められた大きさ以下となるまで繰り返され、当該あらかじめ決められた大きさ以下の各2重対角行列に対して特異値分解が行われた結果である、前記各2重対角行列の特異値と、前記各2重対角行列の特異ベクトルからなる左右直交行列の一部の要素である行列要素とが記憶される特異値分解記憶部から、各2重対角行列の特異値と、行列要素とを読み出し、前記特異値と、前記行列要素とから、分割元の2重対角行列の特異値と、分割元の2重対角行列の行列要素とを算出して前記特異値分解記憶部に蓄積し、その分割元の2重対角行列の特異値と、行列要素とを算出する処理を、2重対角行列Bの少なくとも1個の特異値を算出するまで繰り返し、前記2重対角行列Bの少なくとも1個の特異値を特異値記憶部に蓄積する特異値算出ステップと、前記2重対角行列Bが記憶される対角行列記憶部から2重対角行列Bを読み出し、前記特異値記憶部から前記2重対角行列Bの特異値を読み出し、2重対角行列Bとその特異値とから、ツイスト分解法を用いて2重対角行列Bの少なくとも1個の特異ベクトルを算出して特異ベクトル記憶部に蓄積する特異ベクトル算出ステップと、を実現させるためのものである。
【0219】
また、このプログラムでは、コンピュータに、前記あらかじめ決められた大きさ以下の各2重対角行列も記憶される前記対角行列記憶部から、前記あらかじめ決められた大きさ以下の各2重対角行列を読み出し、前記各2重対角行列に対して特異値分解を行い、前記各2重対角行列の特異値と、前記各2重対角行列の特異ベクトルとを算出し、当該特異値と、当該特異ベクトルからなる左右直交行列の一部の要素である行列要素とを前記特異値分解記憶部に蓄積する特異値分解ステップをさらに実行させてもよい。
【0220】
また、上記実施の形態における特異値分解装置を実現するソフトウェアは、以下のようなプログラムである。つまり、このプログラムは、コンピュータに、2重対角行列Bが記憶される対角行列記憶部から前記2重対角行列Bを読み出し、当該2重対角行列Bを2個の2重対角行列に分割して前記対角行列記憶部に蓄積し、その2重対角行列を2個の2重対角行列に分割して前記対角行列記憶部に蓄積する処理を、分割後の各2重対角行列があらかじめ決められた大きさ以下となるまで繰り返す行列分割ステップと、前記あらかじめ決められた大きさ以下の各2重対角行列を前記対角行列記憶部から読み出し、前記各2重対角行列に対して特異値分解を行い、前記各2重対角行列の特異値と、前記各2重対角行列の特異ベクトルとを算出し、特異値分解された特異値と、特異ベクトルからなる左右直交行列の一部の要素である行列要素とを特異値分解記憶部に蓄積する特異値分解ステップと、各2重対角行列の特異値と、行列要素とを前記特異値分解記憶部から読み出し、前記特異値と、前記行列要素とから、分割元の2重対角行列の特異値と、分割元の2重対角行列の行列要素とを算出して前記特異値分解記憶部に蓄積し、その分割元の2重対角行列の特異値と、行列要素とを算出する処理を、2重対角行列Bの少なくとも1個の特異値を算出するまで繰り返し、前記2重対角行列Bの少なくとも1個の特異値を特異値記憶部に蓄積する特異値算出ステップと、前記対角行列記憶部から2重対角行列Bを読み出し、前記特異値記憶部から前記2重対角行列Bの特異値を読み出し、2重対角行列Bとその特異値とから、ツイスト分解法を用いて2重対角行列Bの少なくとも1個の特異ベクトルを算出して特異ベクトル記憶部に蓄積する特異ベクトル算出ステップと、を実行させるためのものである。
【0221】
また、このプログラムでは、前記特異ベクトル算出ステップは、前記対角行列記憶部から2重対角行列Bを読み出し、前記特異値記憶部から前記2重対角行列Bの特異値を読み出し、前記2重対角行列Bの各要素に関してミウラ変換、dLVv型変換、逆ミウラ変換を行うことによって、前記2重対角行列Bを上2重対角行列B(+1)及び下2重対角行列B(-1)にコレスキー分解するコレスキー分解ステップと、前記上2重対角行列B(+1)及び下2重対角行列B(-1)の各要素と、前記2重対角行列Bの特異値とを用いて一方の左右直交行列を構成する特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する第1特異ベクトル算出ステップと、前記第1特異ベクトル算出ステップで算出した一方の左右直交行列を構成する特異ベクトルと、前記2重対角行列Bの特異値と、前記2重対角行列Bとを用いて他方の左右直交行列を構成する特異ベクトルを算出して前記特異ベクトル記憶部に蓄積する第2特異ベクトル算出ステップと、をさらに備えてもよい。
なお、上記プログラムにおいて、情報を蓄積する処理などでは、ハードウェアでしか行われない処理は少なくとも含まれない。
【0222】
また、このプログラムは、サーバなどからダウンロードされることによって実行されてもよく、所定の記録媒体(例えば、CD-ROMなどの光ディスクや磁気ディスク、半導体メモリなど)に記録されたプログラムが読み出されることによって実行されてもよい。
【0223】
また、このプログラムを実行するコンピュータは、単数であってもよく、複数であってもよい。すなわち、集中処理を行ってもよく、あるいは分散処理を行ってもよい。
【0224】
図30は、上記プログラムを実行して、上記実施の形態による特異値分解装置1を実現するコンピュータの外観の一例を示す模式図である。上記実施の形態は、コンピュータハードウェア及びその上で実行されるコンピュータプログラムによって実現されうる。
【0225】
図30において、コンピュータシステム100は、CD-ROM(Compact Disk Read Only Memory)ドライブ105、FD(Flexible Disk)ドライブ106を含むコンピュータ101と、キーボード102と、マウス103と、モニタ104とを備える。
【0226】
図31は、コンピュータシステムを示す図である。図31において、コンピュータ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への接続を提供する図示しないネットワークカードを含んでもよい。
【0227】
コンピュータシステム100に、上記実施の形態による特異値分解装置1の機能を実行させるプログラムは、CD-ROM121、またはFD122に記憶されて、CD-ROMドライブ105、またはFDドライブ106に挿入され、ハードディスク114に転送されてもよい。これに代えて、そのプログラムは、図示しないネットワークを介してコンピュータ101に送信され、ハードディスク114に記憶されてもよい。プログラムは実行の際にRAM113にロードされる。なお、プログラムは、CD-ROM121やFD122、またはネットワークから直接、ロードされてもよい。
【0228】
また、行列記憶部11、対角行列記憶部13、特異値分解記憶部16、特異値記憶部18,特異ベクトル記憶部20は、RAM113やハードディスク114によって実現されてもよい。
【0229】
プログラムは、コンピュータ101に、上記実施の形態による特異値分解装置1の機能を実行させるオペレーティングシステム(OS)、またはサードパーティプログラム等を必ずしも含まなくてもよい。プログラムは、制御された態様で適切な機能(モジュール)を呼び出し、所望の結果が得られるようにする命令の部分のみを含んでいてもよい。コンピュータシステム100がどのように動作するのかについては周知であり、詳細な説明を省略する。
【0230】
また、本発明は、以上の実施の形態に限定されることなく、種々の変更が可能であり、それらも本発明の範囲内に包含されるものであることは言うまでもない。
また、本発明のほんのいくつかの典型的な実施例について上で詳細に説明したが、その典型的な実施例において、発明の利益と新規な技術から実質的にはずれることなく多くの変更が可能であることを当業者は容易に認識することができるであろう。したがって、そのようなすべての変更は、本発明の範囲に含まれるものである。
【産業上の利用可能性】
【0231】
以上のように、本発明による特異値分解装置等によれば、特異値分解を高速に処理することができ、画像処理や検索処理、その他の特異値分解を用いる処理を実行する装置等において有用である。
【図面の簡単な説明】
【0232】
【図1】本発明の実施の形態1による特異値分解装置の構成を示すブロック図
【図2】同実施の形態による特異値分解装置の動作を示すフローチャート
【図3】同実施の形態による特異値分解装置の動作を示すフローチャート
【図4】同実施の形態による特異値分解装置の動作を示すフローチャート
【図5】同実施の形態による特異値分解装置の動作を示すフローチャート
【図6】同実施の形態における上2重対角行列Bの分割について説明するための図
【図7】同実施の形態における上2重対角行列Bの分割について説明するための図
【図8】同実施の形態による特異値分解装置の動作を示すフローチャート
【図9】同実施の形態における上2重対角行列Bの分割について説明するための図
【図10】同実施の形態における上2重対角行列Bの特異値の算出について説明するための図
【図11】同実施の形態による特異値分解装置の動作を示すフローチャート
【図12】同実施の形態による特異値分解装置の動作を示すフローチャート
【図13】同実施の形態における上2重対角行列Bの特異値の算出について説明するための図
【図14】同実施の形態におけるコレスキー分解について説明するための図
【図15】同実施の形態による特異値分解装置の動作を示すフローチャート
【図16】同実施の形態による特異値分解装置の動作を示すフローチャート
【図17】同実施の形態による特異値分解装置の動作を示すフローチャート
【図18】同実施の形態による特異値分解装置の動作を示すフローチャート
【図19】同実施の形態による特異値分解装置の動作を示すフローチャート
【図20】同実施の形態による特異値分解装置の動作を示すフローチャート
【図21A】同実施の形態におけるqd型ツイスト分解法について説明するための図
【図21B】同実施の形態におけるLV型ツイスト分解法について説明するための図
【図22】同実施の形態による特異値分解装置の動作を示すフローチャート
【図23】同実施の形態による特異値分解装置の構成の他の一例を示すブロック図
【図24】同実施の形態による特異値分解装置の構成の他の一例を示すブロック図
【図25】同実施の形態における画像処理の一例を示すフローチャート
【図26】同実施の形態における文書検索処理の一例を示すフローチャート
【図27】同実施の形態による特異値分解と従来の特異値分解との比較について説明するための図
【図28】同実施の形態による特異値分解と従来の特異値分解との比較について説明するための図
【図29】同実施の形態による特異値分解と従来の特異値分解との比較について説明するための図
【図30】同実施の形態におけるコンピュータシステムの外観一例を示す模式図
【図31】同実施の形態におけるコンピュータシステムの構成の一例を示す図
図面
【図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
【図21A】
20
【図21B】
21
【図22】
22
【図23】
23
【図24】
24
【図25】
25
【図26】
26
【図27】
27
【図28】
28
【図29】
29
【図30】
30
【図31】
31