JP5767576B2

JP5767576B2 - システム同定のための行列格納方法、プログラム及びシステム

Info

Publication number: JP5767576B2
Application number: JP2011277384A
Authority: JP
Inventors: 陽一初鳥; 尚宮下
Original assignee: International Business Machines Corp
Current assignee: International Business Machines Corp
Priority date: 2011-12-19
Filing date: 2011-12-19
Publication date: 2015-08-19
Anticipated expiration: 2031-12-19
Also published as: US20130159373A1; US9164960B2; JP2013127735A

Description

この発明は、パラメトリック・モデルをベースとしたシステム同定に関し、より詳細には、システムズ・モデルをベースとしたシステム同定のための行列格納技法に関する。

従来より、パラメトリック・モデルをベースとしたシステム同定において、予測誤差を含むパラメータ推定のための目的関数が構成され、この目的関数を最小とするようにパラメータを決定する処理が行われている。最も普通に使用される目的関数は二次関数であり、その場合、目的関数を最小とするようにパラメータを決定する処理は、最小二乗法推定と呼ばれる。

同定対象のパラメータを

プラント・モデルを

と定義する。最小二乗法によるパラメータ決定問題とは、観測値f_i (i = 1,...,n)との二乗誤差ε²が最小となるパラメータを決定する問題として下記の式で定義される。ここで、m,nは自然数とし、m < nとする。

そこで、二乗誤差ε²を各パラメータで偏微分した以下のm本の連立方程式からm個のパラメータを同定する問題となり、この連立方程式を正規方程式と定義する。

従って、最小二乗法で解くべき問題は以下のように定義される。

ここで、偏微分の項を行列Aであらわすと以下の式になる。

特に最小二乗法においては、下記の行列の項A^TAは、疎行列になることが知られている。ここで行列Aはm×mとし、a_ijは行列Aのij要素、nは測定回数であるとする。

モデル化する人と解析コードを書く人は一般的には別人であるため、従来は疎行列A^TAが与えられた際、いかに効率よく計算するかという研究が行われてきた。しかしながら、対象とする問題ごとに疎行列の形が大きく異なることや、局所的に異なる性質を持つこともあり、計算時にいくつかの格納形式を試験して格納形式を判断する処理や、事前に格納形式を判断するための実験が必要とされていた。また、疎行列A^TA全体に同じ格納形式を選択するため、効率よく疎行列を格納することは困難になっている。

久保田他, "GPUにおける格納形式自動選択による疎行列ベクトル積の高速化", 情報処理学会研究報告[ハイパフォーマンスコンピューティング], 2010-HPC-128(19), pp. 1-7, 2010は、非零要素率と行毎の非零要素率のばらつきに基づき疎行列における格納形式を決定することを開示する。

R. H. Bisseling et al., “Communication Balancing in Parallel Sparse Matrix-Vector Multiplication”, Electronic Transactions on Numerical Analysis, Volume 21, pp. 47-65, 2005は、マルチプロセッサ環境で、並列処理のために、プロセッサ間の通信コストを下げるように疎行列を分割する技法を開示する。

久保田他, "GPUにおける格納形式自動選択による疎行列ベクトル積の高速化", 情報処理学会研究報告[ハイパフォーマンスコンピューティング], 2010-HPC-128(19), pp. 1-7, 2010 R. H. Bisseling et al., "Communication Balancing in Parallel Sparse Matrix-Vector Multiplication", Electronic Transactions on Numerical Analysis, Volume 21, pp. 47-65, 2005

しかし、上記従来技術の技法は、疎行列に格納する際に、システムズ・モデルの情報を考慮に入れることができないので、結果の疎行列は、計算速度を向上させる点で十分ではなかった。

従って、この発明の目的は、システムズ・モデルの情報を考慮に入れて、システム同定のパラメータ計算のための最小二乗法で使用される疎行列に、計算速度を向上させるように効率よく値を格納する技法を提供することにある。

本発明は、先ず観測が複数の要素からなることを前提とする。観測の要素とは、それぞれが異なるパラメータの組と数からなる観測の部分集合である。そこで、本発明に従うシステムは、最小二乗法で使用される疎行列を、観測の要素の数に応じて小行列に分割する。

一方、観測の要素毎に観測ＩＤが付与され、各パラメータにはパラメータＩＤが付与され、ＩＤマッピングとして要素のパラメータと関連付けられている。そこで本発明に従うシステムは、ＩＤマッピングの有無、観測ＩＤとパラメータＩＤの対応関係、及び小行列の位置に応じて、非ゼロ要素の位置を判断し、それに基づき、各小行列の格納形式を選択する。

ここで、使用される格納形式には、CSR、ELL、DIA、BSR、COO及びDNSがある。本発明に従うシステムは、対角要素であるかどうか、ＩＤマッピングを介さない項が存在するかどうか、同じＩＤマッピングを参照しているかどうか、などの条件に応じて格納形式を選択する。

本発明に従うシステムは、こうして格納形式を選択して要素を行列に格納すると、ＳＣＧ法などの既知の方法で計算して、パラメータの値を求める。

この発明によれば、システムズ・モデルの情報を考慮に入れて、システム同定のパラメータ計算のための最小二乗法で使用される疎行列に要素を格納することにより、システム同定のパラメータ計算の計算速度を向上させるという効果が得られる。

本発明を実施するためのシステムのハードウェア構成を示すブロック図である。本発明を実施するためのシステムの機能構成を示すブロック図である。観測値とパラメータの値の対応を示す図である。ＩＤマッピングの対応を説明する図である。観測値の要素から小行列を生成することと、小行列に対する格納形式を示す図である。格納形式を選択して小行列に要素を格納する処理のフローチャートを示す図である。小行列の位置における定義を示す図である。ＩＤマッピングなしの場合の、観測とパラメータの値の対応を示す図である。ＩＤマッピングありの場合の、観測とパラメータの値の対応を示す図である。ＩＤマッピングなしとありの両方を含む、観測とパラメータの値の対応を示す図である。単一のＩＤマッピングが複数のパラメータの値を含む場合の図である。ＩＤマッピングを複数含む場合の図である。小行列に対する格納結果の例を示す図である。

以下、図面に基づき、この発明の実施例を説明する。特に断わらない限り、同一の参照番号は、図面を通して、同一の対象を指すものとする。尚、以下で説明するのは、本発明の一実施形態であり、この発明を、この実施例で説明する内容に限定する意図はないことを理解されたい。

図１を参照すると、本発明の一実施例に係るシステム構成及び処理を実現するためのコンピュータ・ハードウェアのブロック図が示されている。図１において、システム・パス１０２には、ＣＰＵ１０４と、主記憶（ＲＡＭ）１０６と、ハードディスク・ドライブ（ＨＤＤ）１０８と、キーボード１１０と、マウス１１２と、ディスプレイ１１４が接続されている。ＣＰＵ１０４は、好適には、３２ビットまたは６４ビットのアーキテクチャに基づくものであり、例えば、インテル社のPentium（商標）４、Core(商標)2 Duo、Xeon(商標)、AMD社のAthlon（商標）などを使用することができる。主記憶１０６は、好適には、４ＧＢ以上の容量をもつものである。ハードディスク・ドライブ１０８は、大量のデータを格納できるように、例えば、５００ＧＢ以上の容量をもつものであることが望ましい。

ハードディスク・ドライブ１０８には、個々に図示しないが、オペレーティング・システムが、予め格納されている。オペレーティング・システムは、Linux(商標)、マイクロソフト社のWindows XP(商標)、Windows(商標)7、アップルコンピュータのMac OS（商標）などの、ＣＰＵ１０４に適合する任意のものでよい。

ハードディスク・ドライブ１０８にはさらに、システム同定を行うための観測データ２０４、観測とパラメータの対応を示すＩＤマッピングのデータ２０６、行列格納ルーチン２０８、格納された行列に基づくパラメータ計算ルーチン２１０などが格納されている。行列格納ルーチン２０８及びパラメータ計算ルーチン２１０は、Java(R)、C、C++、C#などの既存の任意のプログラミング言語で作成することができる。これらのパラメータや処理ルーチンの詳細は、図２を参照して後で説明する。

キーボード１１０及びマウス１１２は、オペレーティング・システムまたは、ハードディスク・ドライブ１０８から主記憶１０６にロードされ、ディスプレイ１１４に表示された・メインプログラム２０２上で操作したり、文字を打ち込んだりするために使用される。

ディスプレイ１１４は、好適には液晶ディスプレイであり、例えば、ＸＧＡ（１０２４×７６８の解像度）、またはＵＸＧＡ（１６００×１２００の解像度）などの任意の解像度のものを使用することができる。ディスプレイ１１４は、図示しないが、本発明の処理を開始するための操作ウインドウや、パラメータの計算結果等を表示するために使用される。

次に、図２の機能ブロック図を参照して、本発明の処理の論理的な構成について説明する。

図２において、メイン・プログラム２０２は、全体の処理を統合する機能をもつプログラムであり、ユーザの操作に従い、観測データ２０４としてデータを格納したり、ＩＤマッピングのデータ２０６を作成したり、行列格納ルーチン２０８を起動したり、パラメータ計算ルーチン２１０を起動したり、結果をディスプレイ１１４に表示したりするために使用される。

観測データ２０４は、図３に示すような形式で観測データを保持する。すなわち、観測がm回行われるとすると、図３はi番目(i = 1,...,m)の観測値θ_iのデータ構造を示す。すなわち、各々の観測値は、一般的には各々パラメータが異なる複数の要素からなり、ここではそれを各々要素j1、要素j2、要素j3とする。

各々の要素には、パラメータのセットと、観測値を与える関数が関連付けられている。例えば、要素j1の場合、パラメータP_ik = {p₁₁,...,p_1N1}と、関数f(P_ik)が関連付けられている。

ここで多重度というのは、当該パラメータで測定される対象の数である。例えば、パラメータを(x,y)の座標だとすると、多重度Mというのは、(x₁,y₁),(x₂,y₂),...,(x_M,y_M)というM個のパラメータの組が存在することを意味する。

このような要素j1、要素j2、要素j3により、θ_iは以下の式で計算される。

ＩＤマッピングのデータ２０６とは、観測に使用されるパラメータ・セット毎にパラメータIDを付与し、上記したm回の観測毎にIDを付与し、図４に示すように、i番目の観測の観測ID(i) (i = 1,...,m)の観測データについて、その観測で使用されるパラメータIDをIDマッピングとして対応付けて保存するものである。ＩＤマッピングのデータ２０６からは、ある観測IDで使用されているパラメータのパラメータIDを知ることができる。一方、あるパラメータについて、それに関連付けられた観測IDがあるかどうかも知ることができる。

最小二乗法の計算で使用される疎行列A^TAは、観測値θ_iに含まれる要素に数に応じて、図５のように分割される。すなわち、行列Aが３分割されるなら、A^TAは、図示されているように、３×３＝９個の小行列に分割される。行列格納ルーチン２０８は、図６に関連して後で説明する処理に基づき、観測データ２０４と、ＩＤマッピング２０６を参照して、所定の格納形式で行列要素を格納して、行列データ２１２を用意する機能を有する。

格納形式には以下のような種類がある。
− CSR (Compressed Sparse Row)
疎行列を行方向に走査し、零要素を省いた格納形式であり、非零要素の値、格納された非零要素の列番号、各行の開始位置を格納する。CRS (Compressed Row Storage)と表記されることもある
− ELL (ELLPACK / ITPACK)
疎行列の一辺のサイズをｎ、疎行列１行あたりの非零要素数の最大値をｋとし、疎行列をｎｋの密行列に格納する。１行あたりの非零要素数がｋに満たない場合は0で埋める。
− DIA (Diagonal)
非零な対角要素の値および主対角から各対角へのオフセットを格納する。
− BSR (Block Sparse Row)
r×cの大きさの部分行列（ブロックと呼ぶ）に分割する。BSRはCSRと同様の手順で、少なくとも１つの非零要素が存在する非零ブロックを格納する。非零ブロックの全要素を格納し、非零ブロックのブロック列番号、ブロック行の開始位置を格納する。
− COO (Coordinate)
非零要素の値、行番号、列番号をそれぞれの非零要素ごとに格納する。
− DNS (Dense)
密行列の格納方法。疎行列を配列として格納する。

パラメータ計算ルーチン２１０は、格納された行列データ２１２を使用して、システム同定のためのパラメータを計算し、パラメータ・データ２１４として、好適にはハードディスク・ドライブ１０８に保存する機能を有する。パラメータ計算ルーチン２１０は、好適には、ＳＣＧ(Scaled Conjugate Gradient)法、すなわち、スケーリング共役勾配法により、下記の式を解くことにより、パラメータを計算する。

これの具体的な計算アルゴリズムは後で説明する。

次に、図６のフローチャートを参照して、行列格納ルーチン２０８の処理について説明するが、その前に、小行列の呼称の定義を図７を参照して与える。すなわち、図７において、３×３に分割された小行列の集まりとしてみたとき、対角を占める小行列を、主対角ブロックと呼び、主対角に隣接する小行列をオフセット１のブロック、オフセット１のブロックに隣接して、主対角からさらに離れる小行列をオフセット２のブロックと呼ぶ。図７においてNは、パラメータの値で定義されるパラメータの数である。これは、図３では、P_ik = {p₁₁,...,p_1N1}のぱらメータの数に相当する。なおここでは便宜上、各小行列が、２×２の升目に区切られて示されている。ここでの議論は、３×３に分割された小行列の集まりだけではなく、一般に、より多くの任意の数の小行列の集まりにも適用可能であることを理解されたい。

図６に戻って、行列格納ルーチン２０８は、図６のフローチャートの処理を、分割された小行列に順次適用していくものであることを理解されたい。

ステップ６０２で、行列格納ルーチン２０８は、目下の小行列が対角要素、すなわち主対角かどうか判断する。もし対角要素であるなら、行列格納ルーチン２０８はステップ６０４で、ＩＤマッピングが存在するかどうか判断する。

ここで、ＩＤマッピングが存在しないとは、図８に示すように、観測IDに対してパラメータが定義され、観測IDによらず、常に同じパラメータ・インデックスが使用される場合（例えば、定数項など）である。

一方、ＩＤマッピングが存在するとは、図９に示すように、観測IDに対してパラメータが定義され、ＩＤマッピングを介してパラメータの値が決定される場合である。

ＩＤマッピングが存在するなら、行列格納ルーチン２０８は、ステップ６０６で、N行×N列のブロックを主対角からオフセットMの位置に生成して、ユーザーの設定で、BSRかDIAのどちらかを選択して、要素を格納する。ここでMとは、パラメータ・インデックスの値であり、図３では、M1, M2及びM3に相当する。

ＩＤマッピングが存在しないなら、行列格納ルーチン２０８は、ステップ６０８で、N行×N列の密行列(DNS)が存在すると判断して、要素を格納する。

ステップ６０２に戻って、対角要素でないと判断したなら、行列格納ルーチン２０８は、ステップ６１０で、N_L×N_Rのブロックを生成する。ここで、N_L、N_Rはそれぞれ、非対角要素において、行列積の左要素と右要素で定義されるパラメータの数である。

ステップ６１２では行列格納ルーチン２０８は、ＩＤマッピングを介さない値が存在するかどうか判断する。ここで、ＩＤマッピングを介さない値が存在するとは、図１０に示すような場合である。

ＩＤマッピングを介さない値が存在するなら、行列格納ルーチン２０８は、ステップ６１４で、密行列を格納する。

ＩＤマッピングを介さない値が存在しないなら、行列格納ルーチン２０８は、ステップ６１６で、同じＩＤマッピングを参照するかどうか判断する。

同じＩＤマッピングを参照するとは、例えば図１１に示すような場合である。より具体的に言うと、例えば、測定する対象の位置(x,y)を測定する場合、要素j1がxの測定値、要素j2がyの測定値、要素j3がTの測定値を含むなら、その両方は同じＩＤマッピングを参照することになる。

また、異なるＩＤマッピングを参照するとは、例えば図１２に示すような場合である。例えば、測定する対象の位置(x,y)と温度Tを別個に測定する場合、位置(x,y)は１つのＩＤマッピング１で測定データに関連付けられるが、温度Tは別のＩＤマッピング２で測定データに関連付けられる。

同じＩＤマッピングを参照すると判断すると、行列格納ルーチン２０８はステップ６１８で、ブロックが主対角に存在すると判断して、ユーザーの設定で、BSRまたはDIAを選択する。

もし同じＩＤマッピングを参照していない、すなわち、図１２のような場合なら、行列格納ルーチン２０８はステップ６２０で、N_L > N_Rかどうか判断し、もしそうなら、ステップ６２２でELLで格納し、そうでないなら、ステップ６２４でCSRで格納する。

図１３は、このようにして選択した格納形式の例を示す。

次に、パラメータ計算ルーチン２１０が実行するＳＣＧ法のアルゴリズムについて説明する。このアルゴリズムは例えば、速水、原田、”ベクトル計算機におけるScaled CG法の有効性について”, 情報処理学会研究報告ハイパフォーマンスコンピューティング（HPC）,Vol. 1986, No. 40 (1986-HPC-017), pp. 1-6, 1986に基づくものである。

ここで改めて、Aを疎行列、xを求めるパラメータ・ベクトル、bを観測値ベクトルとすると、
Ax = bという線形連立方程式で書かれる。
そこで、行列Aの対角項を用いてスケーリングを行い、共役勾配法（ＣＧ法）を適用する。そのために先ず、下記の式で初期値を与える。

次に、下記の式を反復することで、xを更新し、収束したら停止する。なお、下記の式で、(A,B)のように表記するのは、ベクトルAとベクトルBの内積を計算するものとして定義する。

パラメータ計算ルーチン２１０は、計算が収束すると、計算結果のパラメータを、パラメータのデータ２１４として、好適にはハードディスク・ドライブ１０８に書き出す。なお、パラメータ計算ルーチン２１０は、ＳＣＧ法に限定されず、ＩＣＣＧ法、ＭＩＣＣＧ法など、ある程度の収束速度が期待できる任意の連立一次方程式計算法を使用してもよい。

以上、特定の実施例に従い本発明を説明してきたが、本発明はこれには限定されず、本発明の技術思想の範囲で、様々な変形例を考慮することができる。例えば、行列のサイズ、パラメータの数は任意でよく、使用するコンピュータのハードウェアも、利用可能な任意のプラット・フォーム、アーキテクチャのものでよい。

１０４・・・ＣＰＵ
１０８・・・ハードディスク・ドライブ
２０４・・・観測データ
２０６・・・ＩＤマッピング
２０８・・・行列格納ルーチン
２１０・・・パラメータ計算ルーチン

Claims

コンピュータの処理により、最小二乗法に基づきパラメトリック・システム同定を行う方法であって、
測定値と、観測モデルと、パラメータを用意するステップと、
前記パラメータに対して、測定値と、パラメータの値を関連付けるＩＤマッピングを作成するステップと、
最小二乗法によりパラメータを計算するために行列にパラメータを格納するステップであって、該行列を小行列に分割し、
(a) 該小行列の対角小行列においては、前記ＩＤマッピングが存在する場合に主対角からオフセットの位置にN×N（Nは、パラメータの値で定義されているパラメータの数）の疎行列の成分をオフセットM（Mはパラメータ・インデックスの値）の位置に生成してBSRまたはDIAを選択し、前記ＩＤマッピングが存在しない場合にN×Nの密行列を選択し、
(b) 前記小行列の対角小行列でない場合で、前記ＩＤマッピングを介さない項が存在するなら、密行列を格納し、
(c) 前記小行列の対角小行列でない場合で、前記ＩＤマッピングを介さない項が存在しないなら、
(c1) 同じＩＤマッピングを参照しているなら、ブロックが主対角に存在し、BSRまたはDIAを選択し、
(c2) 同じＩＤマッピングを参照しているのでないなら、左要素のパラメータ数が右要素のパラメータ数より大きいならELLで要素を格納し、そうでないならCSRで要素を格納するステップを有する、
方法。
前記小行列の対角小行列でない場合、前もって左要素のパラメータ数×右要素のパラメータ数のブロックを生成するステップをさらに有する、請求項１に記載の方法。
格納された上記行列を含む方程式を、スケーリング共役勾配法で解くことによりパラメータを決定するステップをさらに有する、請求項１に記載の方法。
コンピュータの処理により、最小二乗法に基づきパラメトリック・システム同定を行う方法であって、
前記コンピュータに、
測定値と、観測モデルと、パラメータを用意するステップと、
前記パラメータに対して、測定値と、パラメータの値を関連付けるＩＤマッピングを作成するステップと、
最小二乗法によりパラメータを計算するために行列にパラメータを格納するステップであって、該行列を小行列に分割し、
(a) 該小行列の対角小行列においては、前記ＩＤマッピングが存在する場合に主対角からオフセットの位置にN×N（Nは、パラメータの値で定義されているパラメータの数）の疎行列の成分をオフセットM（Mはパラメータ・インデックスの値）の位置に生成してBSRまたはDIAを選択し、前記ＩＤマッピングが存在しない場合にN×Nの密行列を選択し、
(b) 前記小行列の対角小行列でない場合で、前記ＩＤマッピングを介さない項が存在するなら、密行列を格納し、
(c) 前記小行列の対角小行列でない場合で、前記ＩＤマッピングを介さない項が存在しないなら、
(c1) 同じＩＤマッピングを参照しているなら、ブロックが主対角に存在し、BSRまたはDIAを選択し、
(c2) 同じＩＤマッピングを参照しているのでないなら、左要素のパラメータ数が右要素のパラメータ数より大きいならELLで要素を格納し、そうでないならCSRで要素を格納するステップを実行させる、
プログラム。
前記コンピュータに、前記小行列の対角小行列でない場合、前もって左要素のパラメータ数×右要素のパラメータ数のブロックを生成するステップをさらに実行させる、請求項４に記載のプログラム。
前記コンピュータに、格納された上記行列を含む方程式を、スケーリング共役勾配法で解くことによりパラメータを決定するステップをさらに実行させる、請求項４に記載のプログラム。
コンピュータの処理により、最小二乗法に基づきパラメトリック・システム同定を行うシステムであって、
記憶手段と、
前記記憶手段に保存された、測定値と、観測モデルと、パラメータのデータと、
前記記憶手段に保存された、前記パラメータに対して、測定値と、パラメータの値を関連付けるＩＤマッピングのデータと、
最小二乗法によりパラメータを計算するために行列にパラメータを格納するステップであって、該行列を小行列に分割し、
(a) 該小行列の対角小行列においては、前記ＩＤマッピングが存在する場合に主対角からオフセットの位置にN×N（Nは、パラメータの値で定義されているパラメータの数）の疎行列の成分をオフセットM（Mはパラメータ・インデックスの値）の位置に生成してBSRまたはDIAを選択し、前記ＩＤマッピングが存在しない場合にN×Nの密行列を選択し、
(b) 前記小行列の対角小行列でない場合で、前記ＩＤマッピングを介さない項が存在するなら、密行列を格納し、
(c) 前記小行列の対角小行列でない場合で、前記ＩＤマッピングを介さない項が存在しないなら、
(c1) 同じＩＤマッピングを参照しているなら、ブロックが主対角に存在し、BSRまたはDIAを選択し、
(c2) 同じＩＤマッピングを参照しているのでないなら、左要素のパラメータ数が右要素のパラメータ数より大きいならELLで要素を格納し、そうでないならCSRで要素を格納する手段を有する、
システム。
前記小行列の対角小行列でない場合、前もって左要素のパラメータ数×右要素のパラメータ数のブロックを生成する手段をさらに有する、請求項７に記載のシステム。
格納された上記行列を含む方程式を、スケーリング共役勾配法で解くことによりパラメータを決定する手段をさらに有する、請求項７に記載のシステム。