"複数の第1オブジェクトの各々と複数の第2オブジェクトの各々との間の関係の強さを示す関係データの時系列を表わし、かつ、テンソルで表現される時間関係データを解析して、前記時間関係データを表わす未知のテンソルXを求める時間関係データ解析装置であって、前記第1オブジェクトと前記第2オブジェクトとの間の関係の強さのM個の観測値からなる観測ベクトルy、正の実数である正則化定数λ、及びテンソルXの各モードkk=1、2、3に対する相対的な正則化の強さを表わす正則化バランス定数γkを記憶した記憶手段と、前記記憶手段から、前記観測ベクトルy、前記正則化定数λ、前記正則化バランス定数γkk=1、2、3を取得し、以下のI式で表される目的関数を最小化する前記テンソルXを求める解析手段と、を含む時間関係データ解析装置。ただし、n1,n2,n3は、前記時間関係データのサイズを表す自然数であり、xは、前記テンソルXの要素である。ΩXは、テンソルXから前記観測ベクトルの観測値に対応するM個の要素の取り出す操作である。mat1は、n1×n2n3を返すモード1行列化であり、mat2は、n2×n3n1を返すモード2行列化であり、mat3は、n3×n1n2を返すモード3行列化である。∥∥trは、行列の特異値の線形和を表わす。前記記憶手段は、補助変数Zkk=1、2、3、前記テンソルXの各要素の値x、ラグランジュ乗数Akk=1、2、3、第1オブジェクト方向、第2オブジェクト方向、及び時間方向の各々のパターン基底Ukk=1、2、3、前記パターン基底Ukk=1、2、3に従って変化する前記第1オブジェクト、前記第2オブジェクト、及び時間の各ペアの間の関係の強さを表わす関係パターンVkk=1、2、3を更に記憶し、前記解析手段は、各モードkの補助係数Zk及びラグランジュ乗数Akに基づいて、以下のII式に従ってテンソルの予測値Xを算出し、算出されたテンソルの予測値X及び前記観測ベクトルyに基づいて、以下のIII式に従って前記テンソルXの各要素の値xを更新して前記記憶手段に記憶する主変数更新手段と、前記テンソルX及び各モードkの前記ラグランジュ乗数Akに基づいて、各モードkについて、以下のIV式に従って、補助変数の予測値Zkを算出し、算出された前記補助変数の予測値Zkに基づいて、各モードkについて、以下のV式に従って、特異値分解を行って、特異値が並んだ対角行列Skを求めると共に、前記パターン基底Uk及び前記関係パターンVkを更新して前記記憶手段に記憶し、前記求められた前記対角行列Sk、前記パターン基底Uk、及び前記関係パターンVkに基づいて、以下のVI式に従って、前記補助変数Zkを更新して、前記記憶手段に記憶する補助変数更新手段と、前記テンソルX及び補助変数Zkに基づいて、各モードkについて、以下のVII式に従って、前記ラグランジュ乗数Akを更新して前記記憶手段に記憶するラグランジュ乗数更新手段と、前記主変数更新手段による更新、前記補助変数更新手段による更新、前記ラグランジュ乗数更新手段による更新を繰り返すことで、所定の解析終了条件を満たした場合に、前記記憶手段に記憶されている前記テンソルXを出力する解析終了条件判定手段と、を含む請求項1記載の時間関係データ解析装置。ただし、 ̄ΩXは、前記テンソルXから前記観測値でない要素の取り出す操作であり、ηは定数である。前記解析手段によって求められたテンソルXに基づいて、以下のVIII式に従って、前記第1オブジェクトiと前記第2オブジェクトjとの間の関係における予測対象時刻t’の関係データの予測値X’i,jを計算する関係データ予測手段を更に含む請求項1又は2記載の時間関係データ解析装置。ただし、予測対象時刻t’の前後で観測値が観測された時刻を、Tk、Tk1とする。第1オブジェクトと第2オブジェクトとの間の関係の強さのM個の観測値からなる観測ベクトルy、正の実数である正則化定数λ、及びテンソルXの各モードkk=1、2、3に対する相対的な正則化の強さを表わす正則化バランス定数γkを記憶した記憶手段を備えた時間関係データ解析装置において、複数の第1オブジェクトの各々と複数の第2オブジェクトの各々との間の関係の強さを示す関係データの時系列を表わし、かつ、テンソルで表現される時間関係データを解析して、前記時間関係データを表わす未知のテンソルXを求める時間関係データ解析方法であって、前記解析手段によって、前記記憶手段から、前記観測ベクトルy、前記正則化定数λ、前記正則化バランス定数γkk=1、2、3を取得し、以下のI式で表される目的関数を最小化する前記テンソルXを求めることを特徴とする時間関係データ解析方法。ただし、n1,n2,n3は、前記時間関係データのサイズを表す自然数であり、xは、前記テンソルXの要素である。ΩXは、テンソルXから前記観測ベクトルの観測値に対応するM個の要素の取り出す操作である。mat1は、n1×n2n3を返すモード1行列化であり、mat2は、n2×n3n1を返すモード2行列化であり、mat3は、n3×n1n2を返すモード3行列化である。∥∥trは、行列の特異値の線形和を表わす。前記記憶手段は、補助変数Zkk=1、2、3、前記テンソルXの各要素の値x、ラグランジュ乗数Akk=1、2、3、第1オブジェクト方向、第2オブジェクト方向、及び時間方向の各々のパターン基底Ukk=1、2、3、前記パターン基底Ukk=1、2、3に従って変化する前記第1オブジェクト、前記第2オブジェクト、及び時間の各ペアの間の関係の強さを表わす関係パターンVkk=1、2、3を更に記憶し、前記解析手段によって前記テンソルXを求めることは、主変数更新手段によって、各モードkの補助係数Zk及びラグランジュ乗数Akに基づいて、以下のII式に従ってテンソルの予測値Xを算出し、算出されたテンソルの予測値X及び前記観測ベクトルyに基づいて、以下のIII式に従って前記テンソルXの各要素の値xを更新して前記記憶手段に記憶するステップと、前記補助変数更新手段によって、前記テンソルX及び各モードkの前記ラグランジュ乗数Akに基づいて、各モードkについて、以下のIV式に従って、補助変数の予測値Zkを算出し、算出された前記補助変数の予測値Zkに基づいて、各モードkについて、以下のV式に従って、特異値分解を行って、特異値が並んだ対角行列Skを求めると共に、前記パターン基底Uk及び前記関係パターンVkを更新して前記記憶手段に記憶し、前記求められた前記対角行列Sk、前記パターン基底Uk、及び前記関係パターンVkに基づいて、以下のVI式に従って、前記補助変数Zkを更新して、前記記憶手段に記憶するステップと、ラグランジュ乗数更新手段によって、前記テンソルX及び補助変数Zkに基づいて、各モードkについて、以下のVII式に従って、前記ラグランジュ乗数Akを更新して前記記憶手段に記憶するステップと、前記解析終了条件判定手段によって、前記主変数更新手段による更新、前記補助変数更新手段による更新、前記ラグランジュ乗数更新手段による更新を繰り返すことで、所定の解析終了条件を満たした場合に、前記記憶手段に記憶されている前記テンソルXを出力するステップと、を含む請求項4記載の時間関係データ解析方法。ただし、 ̄ΩXは、前記テンソルXから前記観測値でない要素の取り出す操作であり、ηは定数である。関係データ予測手段によって、前記解析手段によって求められたテンソルXに基づいて、以下のVIII式に従って、前記第1オブジェクトiと前記第2オブジェクトjとの間の関係における予測対象時刻t’の関係データの予測値X’i,jを計算することを更に含む請求項4又は5記載の時間関係データ解析方法。ただし、予測対象時刻t’の前後で観測値が観測された時刻を、Tk、Tk1とする。請求項1~請求項3の何れか1項に記載の時間関係データ解析装置を構成する各手段として、コンピュータを機能させることを特徴とするプログラム。以上の手法はそれぞれ以下のような問題点を持つ。まず、非特許文献1、非特許文献2、及び非特許文献6に記載の手法では、時間関係データを解析することが不可能である、という問題がある。非特許文献3及び非特許文献5に記載の手法は、時間関係データをパターンに分割することが可能であるが、事前にパターンの数を決定する必要がある、という問題がある。非特許文献4に記載の手法は、パターンの数を自動的に決定することができる時間関係データ解析法である。しかしながら、オブジェクト間の関係の有無だけを解析の対象としている。したがって、例えば関係性の強弱の変化といった、より正確な関係の解析には不適切である。また、この手法では観測が行われてない時刻の推定は不可能である。たとえば5年毎に調査される産業連関表のような時間関係データが与えられた場合に、他の時期3年目などの関係の復元や、未来における関係の予測は原理上不可能である、という問題がある。本発明は、上記の課題を解決するためになされたもので、時間変化するオブジェクト間の関係データを解析して、未観測のオブジェクト間の関係の強さを求めることができる時間関係データ解析装置、方法、及びプログラムを提供することを目的とする。上記の目的を達成するために時間関係データ解析装置は、複数の第1オブジェクトの各々と複数の第2オブジェクトの各々との間の関係の強さを示す関係データの時系列を表わし、かつ、テンソルで表現される時間関係データを解析して、前記時間関係データを表わす未知のテンソルXを求める時間関係データ解析装置であって、前記第1オブジェクトと前記第2オブジェクトとの間の関係の強さのM個の観測値からなる観測ベクトルy、正の実数である正則化定数λ、及びテンソルXの各モードkk=1、2、3に対する相対的な正則化の強さを表わす正則化バランス定数γkを記憶した記憶手段と、前記記憶手段から、前記観測ベクトルy、前記正則化定数λ、前記正則化バランス定数γkk=1、2、3を取得し、以下のI式で表される目的関数を最小化する前記テンソルXを求める解析手段と、を含んで構成されている。ただし、n1,n2,n3は、前記時間関係データのサイズを表す自然数であり、xは、前記テンソルXの要素である。ΩXは、テンソルXから前記観測ベクトルの観測値に対応するM個の要素の取り出す操作である。mat1は、n1×n2n3を返すモード1行列化であり、mat2は、n2×n3n1を返すモード2行列化であり、mat3は、n3×n1n2を返すモード3行列化である。∥∥trは、行列の特異値の線形和を表わす。本発明に係る時間関係データ解析方法は、第1オブジェクトと第2オブジェクトとの間の関係の強さのM個の観測値からなる観測ベクトルy、正の実数である正則化定数λ、及びテンソルXの各モードkk=1、2、3に対する相対的な正則化の強さを表わす正則化バランス定数γkを記憶した記憶手段を備えた時間関係データ解析装置において、複数の第1オブジェクトの各々と複数の第2オブジェクトの各々との間の関係の強さを示す関係データの時系列を表わし、かつ、テンソルで表現される時間関係データを解析して、前記時間関係データを表わす未知のテンソルXを求める時間関係データ解析方法であって、前記解析手段によって、前記記憶手段から、前記観測ベクトルy、前記正則化定数λ、前記正則化バランス定数γkk=1、2、3を取得し、上記のI式で表される目的関数を最小化する前記テンソルXを求めることを特徴としている。本発明によれば、前記解析手段によって、前記記憶手段から、前記観測ベクトルy、前記正則化定数λ、前記正則化バランス定数γkk=1、2、3を取得し、上記のI式で表される目的関数を最小化する前記テンソルXを求める。このように、時間変化するオブジェクト間の関係データを解析して、目的関数を最小化するように、時間関係データを表わす未知のテンソルXを求めることによって、未観測のオブジェクト間の関係の強さを求めることができる。本発明に係る前記記憶手段は、補助変数Zkk=1、2、3、前記テンソルXの各要素の値x、ラグランジュ乗数Akk=1、2、3、第1オブジェクト方向、第2オブジェクト方向、及び時間方向の各々のパターン基底Ukk=1、2、3、前記パターン基底Ukk=1、2、3に従って変化する前記第1オブジェクト、前記第2オブジェクト、及び時間の各ペアの間の関係の強さを表わす関係パターンVkk=1、2、3を更に記憶し、前記解析手段は、各モードkの補助係数Zk及びラグランジュ乗数Akに基づいて、以下のII式に従ってテンソルの予測値Xを算出し、算出されたテンソルの予測値X及び前記観測ベクトルyに基づいて、以下のIII式に従って前記テンソルXの各要素の値xを更新して前記記憶手段に記憶する主変数更新手段と、前記テンソルX及び各モードkの前記ラグランジュ乗数Akに基づいて、各モードkについて、以下のIV式に従って、補助変数の予測値Zkを算出し、算出された前記補助変数の予測値Zkに基づいて、各モードkについて、以下のV式に従って、特異値分解を行って、特異値が並んだ対角行列Skを求めると共に、前記パターン基底Uk及び前記関係パターンVkを更新して前記記憶手段に記憶し、前記求められた前記対角行列Sk、前記パターン基底Uk、及び前記関係パターンVkに基づいて、以下のVI式に従って、前記補助変数Zkを更新して、前記記憶手段に記憶する補助変数更新手段と、前記テンソルX及び補助変数Zkに基づいて、各モードkについて、以下のVII式に従って、前記ラグランジュ乗数Akを更新して前記記憶手段に記憶するラグランジュ乗数更新手段と、前記主変数更新手段による更新、前記補助変数更新手段による更新、前記ラグランジュ乗数更新手段による更新を繰り返すことで、所定の解析終了条件を満たした場合に、前記記憶手段に記憶されている前記テンソルXを出力する解析終了条件判定手段と、を含むようにすることができる。ただし、 ̄ΩXは、前記テンソルXから前記観測値でない要素の取り出す操作であり、ηは定数である。本発明に係る時間関係データ解析装置は、前記解析手段によって求められたテンソルXに基づいて、以下のVIII式に従って、前記第1オブジェクトiと前記第2オブジェクトjとの間の関係における予測対象時刻t’の関係データの予測値X’i,jを計算する関係データ予測手段を更に含むようにすることができる。ただし、予測対象時刻t’の前後で観測値が観測された時刻を、Tk、Tk1とする。本発明に係るプログラムは、上記の時間関係データ解析装置の各手段としてコンピュータを機能させるためのプログラムである。以上説明したように、本発明の時間関係データ解析装置、方法、及びプログラムによれば、時間変化するオブジェクト間の関係データを解析して、目的関数を最小化するように、時間関係データを表わす未知のテンソルXを求めることによって、未観測のオブジェクト間の関係の強さを求めることができる、という効果が得られる。以下、図面を参照して本発明の実施の形態を詳細に説明する。<概要>まず、本発明で提案する時間関係データ解析装置で用いる時間関係データについて説明する。図1は、提案する時間関係データ解析装置のコンセプトを模式的に表した図である。図1において、左側にはユーザ購買者と商品の間の関係買ったか否かが時間軸に沿ってテンソルとして表現されている。黒い立方体は観測されている関係を表し、テンソルのそれ以外の要素は未知であるとする。図1の右側には左側の部分観測テンソルの低ランク分解を示している。中心を取り囲む3つの薄い直方体はそれぞれ商品、ユーザ、および時間変化を少数のパターンで表現したものであり、中心の小さい直方体はパターン間の相互作用を表現している。このように分解することにより、左側のテンソルの未知部分の推定も可能となる。<システム構成>次に、観測データを入力として、時間関係データを表わす未知のテンソルを求めて出力する時間関係データ解析装置に、本発明を適用した場合を例にして、本発明の実施の形態を説明する。図2に示すように、本実施の形態に係る時間関係データ解析装置は、時間関係データの観測値の入力を受け付ける入力部1と、観測された時間関係データの解析を行う解析部2と、変数や定数を記憶する記憶部3と、予測値計算部4と、解析結果を出力する出力部5と、を備えている。入力部1は、既知のビデオカメラ、マイクロフォン、記憶装置などの入力器により実現される。また、出力部5は、ディスプレイ、プリンタ、磁気ディスクなどで実装される。解析部2、記憶部3、及び予測値計算部4は、CPUCentralProcessingUnitと、RAMRandomAccessMemoryと、後述する時間関係データ解析処理ルーチンを実行するためのプログラムを記憶したROMReadOnlyMemoryとを備えたコンピュータで構成されている。解析部2は、初期設定部21、変数計算部22、及び解析終了条件判定部27を備えている。変数計算部22は、主変数更新部23、補助変数更新部24、ラグランジュ乗数更新部25、及び誤差計算部26を備えている。初期設定部21は、入力された観測データから、記憶部3の各要素に値を設定する。変数計算部22は、記憶部3に保存された情報を用いて各変数の更新を行う。主変数更新部23は、記憶部3に保存された情報を用いて主変数の更新を行う。補助変数更新部24では、記憶部3に保存された情報を用いて、補助変数の更新とパターン基底の更新と関係パターンの更新とを行う。ラグランジュ乗数更新部25は、記憶部3に保存された情報を用いてラグランジュ乗数の更新を行う。誤差計算部26では、記憶部3に保存された情報を用いて誤差の計算を行う。解析終了条件判定部27では、変数計算部22を用いて計算され記憶部3に記憶された数値をもちいて解析を終了してよいかどうかを判定する。記憶部3は、定数31、観測値32、インデックス集合33、主変数34、補助変数35、ラグランジュ乗数36、誤差37、パターン基底38、関係パターン39から構成される。予測値計算部4では、指定された予測時刻t’における関係データの予測を行う。<データ構造>n1、n2、n3は想定する時間関係データのサイズを表す自然数である。関係データは2種類のオブジェクトの間の接続関係を表現する。本実施の形態では、関係を構成する2種類のオブジェクトの一方を起点オブジェクト、他方を終点オブジェクトと呼ぶ。例えば、図1の例では商品が起点オブジェクト、ユーザが終点オブジェクトとなる。また、本実施の形態では、n1、n2はそれぞれ関係の起点と終点のオブジェクト数とし、n3は観測が発生した時間の総ステップ数とする。また、可能なすべての要素数をNとし、N=n1×n2×n3とする。また、起点となるオブジェクトの軸をテンソルの1番目の次元、終点となるオブジェクトの軸をテンソルの2番目の次元、時間の軸をテンソルの3番目の次元と呼ぶ。Tはn3個の観測された時間点の時刻を格納するベクトルである。Tの要素は昇順に並んでいるものとする。観測ベクトルyは観測された関係データであり長さMの縦ベクトルで表現される。Mはyの要素数であり、ymはyのm番目1≦m≦Mの観測値を表す。なお、本実施の形態では観測値ymはすべての実数を許容する。したがって、この値はある関係の強さを表すと理解できる。一方、上記の非特許文献3や非特許文献4では観測される関係データは1あるいは0の2値しか許容しないため、関係の有無しか表現できず、関係の強さを考慮した解析が不可能である。I,J,Kはそれぞれ長さMのインデックス集合33で、それぞれのm番目の要素を並べたが未知テンソルのm番目の観測された要素を指し示す。すなわちm番目の観測値はim番目の起点オブジェクトとjm番目の終点オブジェクトの時刻kmでの関係である。また同様にIcJcKcは、それぞれ長さN-Mのインデックス集合33で、それぞれのm番目の要素の並べたものは未知テンソルのm番目の未観測要素を指し示す。λは正則化定数とよばれる正の実数である。分解解析計算における正則化項の効果の大きさを表す。M次元の観測ベクトルyは、観測値32に格納される。主変数34にはn1×n2×n3の要素からなる未知テンソルXが格納される。初期値としてXの各要素には適当な値を入力する。一般にはすべての値を0とすれば良い。未知テンソルXは、上記図1に示すテンソルに相当するものであり、最終的には、要素に、時刻kにおいて観測される起点オブジェクトiと終点オブジェクトjとの間の関係の強さを示す値が格納上書きされることになる。以下に観測要素の取り出しΩ、未観測要素の取り出し ̄Ω、観測要素の復元ΩT、モードk行列化matk、モードkテンソル化tensorkと呼ぶ4つの操作を定義する。これらの操作は以下の数式では行列演算として表現されるが、実際には行列として保持する必要はない。まず、観測要素の取り出しΩは、引数としてn1×n2×n3の要素からなるテンソルを取り、M次元ベクトルを返すもので、戻り値のm番目の要素は、引数であるテンソルの第要素である。同様に、未観測要素の取り出し ̄Ωは、引数としてn1×n2×n3の要素からなるテンソルを取り、N-M次元ベクトルを返すもので、戻り値のm番目の要素は、引数であるテンソルの第要素である。観測要素の復元ΩTは観測要素の取り出しΩの逆の操作として定義され、M次元ベクトルを引数としてn1×n2×n3の要素からなるテンソルを返すもので、引数のm番目の要素は、戻り値であるテンソルの第要素に格納される。観測値に対応しない戻り値の要素はすべてゼロである。モードk行列化matkはk=123の3種存在し、それぞれn1×n2×n3の要素からなるテンソルを引数として、n1×n2n3の行列、n2×n3n1の行列、n3×n1n2の行列を返すものである。具体的にはモード1行列化mat1は図3Aに示すように与えられたテンソルの第1モード、すなわち起点オブジェクトが縦に連続するようにベクトルを取り出し、横にならべて行列にしたものを生成する操作である。また、第1モードテンソル化tensor1はその逆の操作である。同様に、第2モード行列化mat2は図3Bに示すように終点オブジェクトが縦に連続するようにベクトルを取り出し、横にならべて行列にしたものを生成する操作であり、第2モードテンソル化tensor2はその逆の操作である。第3モード行列化mat3は図3Cのように時間点が縦に連続するようにベクトルを取り出し、横に並べて行列にしたものを生成する操作であり、第3モードテンソル化tensor3はその逆の操作である。<入力部1>入力部1では、解析にあたり事前に必要となる数値データの入力を受け付ける。入力されるデータは以下の通りである。n1、n2、n3、N、M、I,J,K,Ic、Jc,Kc、y、λ、γ1、γ2、γ3、η、L,E、t’n1、n2、n3は想定する時間関係データのサイズを表す自然数である。可能なすべての要素数をNとし、N=n1×n2×n3とする。yは観測された関係データであり長さMの縦ベクトルで表現される観測ベクトルである。Mはyの要素数であり、ymはm番目の観測値を表す。I,J,Kは観測値のインデックスであり、M個の観測値yと未知テンソルXの関係を表す。同様にN-M個の未観測値のインデックスであるIc、Jc,Kcも入力される。λは正則化定数とよばれる正の実数である。分解解析計算における正則化項の効果の大きさを表す。γ1、γ2、γ3、は未知テンソルXの各モードkに対する相対的な正則化の強さを表現する正の定数である。ここではこれらを正則化バランス定数と呼ぶ。また、正則化定数λをある定数倍して正則化バランス定数γ1、γ2、γ3をその定数で割っても得られる解は同じであるため、規格化のためγ1γ2γ3=1と仮定する。ηはステップサイズであり、正の実数である。分解解析計算においてラグランジュ乗数A1、A2、A3の更新量を決める定数であり、主変数と補助変数の乖離に対するペナルティーの強さを表す。Lは変数計算部22における、最大繰り返し回数を表す自然数である。繰り返し回数はlで表す。Eは解析終了条件判定部27における、許容される誤差の上限を表す定数である。t’は予測したい時刻のインデックスを表す実数である。なお、予測が不要な場合にはt’は入力しなくても良い。<解析部2>解析部2は、入力された観測値や定数をもとにして、主変数である未知テンソルXを求める。解析部2は、以下の1式で示される目的関数を、近似的に最小化することにより主変数Xを決定する。ここで、上記1式の第1項は損失項であり、観測ベクトルyと、未知テンソルXから観測値に対応する要素を取り出して得られるM次元ベクトルΩとの差の二乗和として定義される。ただし、観測値に欠損値がある場合、第1項だけでは未知テンソルXを一意に定めることができない。第2項は正則化項と呼ばれ、未知テンソルのランクが低いという仮定を表現する。ただし、ここでZktrはトレースノルムと呼ばれ、行列Zkの特異値の線形和として定義される。ここで未知テンソルが上記図1のように少ない数のパターンとそれらの間の相互作用に分解されることは、行列化操作matkによって得られる行列が低ランクであることに対応することに注意する例えば非特許文献8TGKoldaandBWBaderTensordecompositionsandapplicationsSIAMReview51455-5002009を参照。なぜトレースノルム正則化が低ランク性の仮定に対応しているかを説明するために、全要素が観測されている場合で簡単のために第1モードに対応する正則化項のみ考えるつまりγ1=1γ2=0γ3=0。この時、上記の最適化問題は、以下の2式のように書き直すことができる。ここで、X、Yはそれぞれ、未知テンソルX、および観測ベクトルyのモード1行列化によって得られる行列とする。上の最適化問題の解は、行列Yの特異値分解をY=USVTとして、解析的に以下の3式のように書けることが知られている例えば、非特許文献9RTomiokaTSuzukiandMSugiyamaAugmentedlagrangianmethodsforlearningselectingandcombiningfeaturesInSSraSNowozinandSJWrighteditorsOptimizationforMachineLearningMITPress2011を参照。ここでSは行列Yの特異値が並んだ対角行列である。U、Vは各列が特異値に対応する特異ベクトルであるような行列である。また、maxの操作は行列の要素ごとに適用されるものとする。上記3式からこの最適化問題を解くと、図4に示すように、与えられた行列Yの特異値のうち正則化定数λ以下の特異値は雑音として切り捨てられ、雑音より十分大きい信号成分を取り出すことができることがわかる。また、上記の操作は観測値Yの両側から直交行列をかけても不変であり、これらの操作に対して頑健になっている。上記最適化問題は、たとえば非特許文献7SBoydNParikhEChuBPeleatoandJEcksteinDistributedoptimizationandstatisticallearningviathealternatingdirectionmethodofmultipliershttpwwwstanfordeduboydpapersdistroptstatlearningadmmhtmlのようなAlternatingDirectionMethodofMultipliers法を用いて解くことができる。さらに、上記最適化問題は凸最適化問題凸最適化に関する一般論は非特許文献10SBoydandLVandenbergheConvexOptimizationCambridgeUniversityPress2004を参照なので、ADMM法が収束すれば、それが上記の最適化問題の大域最適解である。ADMM法の収束は下で詳述するようにKKT条件の誤差を用いて判定する。これにより、テンソルを幾つのパターンに分解するか次元を指定する必要がなく、最適化を通してこれらを自動的に決定することが可能となる。<初期設定部21>初期設定部21では、入力部1から入力された情報をもとに、記憶部3に記憶されている各要素の初期設定を行う。まず、入力部1よりn1、n2、n3、N、M、y、λ、γ1、γ2、γ3、η、L,E、t’を定数31に記憶する。また繰り返し回数lを1に初期化する。M次元の観測ベクトルyを観測値32に、観測値のインデックスI,J,Kおよび未観測値のインデックスIc、Jc,Kcをインデックス集合33に保存する。主変数34にはn1×n2×n3の要素からなる未知テンソルXを記憶する。初期値としてXの各要素には適当な値を入力する。一般にはすべての値を0とすれば良い。補助変数35にはZ1,Z2,Z3が記憶され、それぞれ大きさn1×n2n3、n2×n3n1、n3×n1n2の行列として表現される。そのベクトル化をそれぞれz1,z2,z3と書く。補助変数は変数計算部22において主変数xに対する制約を緩和したものである。具体的にはZkは最適化が進むに従ってテンソルXのモードk行列化と一致していく量である。初期値として補助変数Z1,Z2,Z3には適当な値を入力する。一般にはすべての値を0とすれば良い。ラグランジュ乗数36には、ラグランジュ乗数A1、A2、A3が記憶され、ラグランジュ乗数A1、A2、A3は補助変数Z1、Z2、Z3と同じ次元を持つ行列である。ラグランジュ乗数は主変数Xと補助変数Z1、Z2、Z3を一致させる役割を持つ。初期値としてラグランジュ乗数A1、A2、A3には適当な値を入力する。一般にはすべての値を0とすれば良い。誤差37には、ε1、ε2、、εTが記憶され、ε1、ε2、、εTは誤差計算部26において計算される誤差を表す。パターン基底38には、パターン基底U1、U2、U3が記憶される。パターン基底U1、U2、U3はそれぞれn1×r1、n2×r2、n3×r3の行列であり、変数計算部22で得られた起点オブジェクト方向、終点オブジェクト方向、および時間方向のパターンを表現する。関係パターン39には、関係パターンV1、V2、V3が記憶される。関係パターンV1、V2、V3はそれぞれn2n3×r1、n3n1×r2、n1n2×r3の行列であり、パターン基底Ukに従って変化する、終点オブジェクトおよび時間の間の関係の強さ、起点オブジェクトおよび時間の間の関係の強さ、及び起点オブジェクトおよび終点オブジェクトの間の関係の強さを表現する。r1、r2、r3はそれぞれパターンの要素数を表す。この値はテンソルの分解解析が進むにつれて自動的に決定される。初期値として、パターン基底U1、U2、U3、関係パターンV1、V2、V3、パターンの要素数r1、r2、r3には適当な値を入力する。一般にはパターン基底はすべての値を0とすれば良い。r1、r2、r3は適当な自然数を選択する。<変数計算部22>変数計算部22は、記憶部3に記憶された各要素を利用して、ADMM法に基づいた主変数の推定を行う。実際の計算では、主変数更新部23、補助変数更新部24、ラグランジュ乗数更新部25、及び誤差計算部26の各処理をこの順番で行う解析処理を繰り返し行う。繰り返しごとに解析終了条件判定部27で解析の終了を判定し、解析終了とならない場合は繰り返し回数tを1増加させ、再度解析処理を行う。<主変数更新部23>主変数更新部23は、以下の4式および5を用いて主変数Xの更新を行う。まず、主変数更新部23は、上記の4式を用いて主変数の予測値Xを構成する。予測値Xは、主変数と同じサイズのテンソルであり、各モードkについて補助変数Zkからラグランジュ乗数Akを差し引き、テンソル化したものの平均である。次に、主変数更新部23は、構成した主変数の予測値Xを用いて、上記5式のように主変数Xを更新する。上記5式では、観測値の得られている要素については、新しい主変数の値が、その予測値と観測値を3λη:1の割合で内分したものに等しくなるように更新している。一方、観測値の得られていない要素については、新しい主変数の値がその予測値そのものとなるように更新している。得られた主変数Xを記憶部3の主変数34に保存して、この主変数更新部23の処理は終了である。<補助変数更新部24>補助変数更新部24は、以下の6式、7式、および8式を用いて、補助変数Z1,Z2,Z3、パターン基底U1、U2、U3、及び関係パターンV1,V2,V3の更新を行う。まず、上記6式のように、記憶部3に保存されている量を用いて補助変数の予測値^Zkを構成する。補助変数の予測値^Zkは各モードkについて主変数Xを行列化したものにラグランジュ乗数Akを足したものである。次に、このように構成した補助変数の予測値^Zkに対し、上記7式に従って特異値分解を行う。特異値分解を行うと、k番目の予測値^Zkに対し、Uk、Sk、Vkという行列が得られる。このうちSkは特異値が並んだ対角行列である。Uk、Vkはそれぞれの列が特異値に対応する特異ベクトルである。ただし、特異値及び特異ベクトルとしては、次の8式の操作のために、γk/η以上の値を持つ特異値とそれに対応する特異ベクトルのみ計算すればよい。最後に、上記7式で得られた特異値と特異ベクトルを用いて、上記8式に従って、補助変数Zkを更新する。また、新しいパターン基底Ukは、max関数によって0となる列を除いたものである。すなわち、パターンの要素数rkは特異値分解の結果によって繰り返しのたびに変動しうる。得られた補助変数Zkをベクトル化して、記憶部3の補助変数35に保存し、パターン基底Ukを記憶部3のパターン基底38に保存し、さらに関係パターンVkを記憶部3の関係パターン39に保存して、この補助変数更新部24の解析処理は終了となる。<ラグランジュ乗数更新部25>ラグランジュ乗数更新部25は、以下の9式を用いてラグランジュ乗数A1、A2、A3の更新を行う。上記9式では、ラグランジュ乗数A1、A2、A3を、主変数Xを行列化したものとそれぞれの補助変数Z1、Z2、Z3の乖離の量だけ更新する。得られたラグランジュ乗数を記憶部3のラグランジュ乗数36に保存して、このラグランジュ乗数更新部25の解析処理は終了である。<誤差計算部26>誤差計算部26では、KKT条件誤差を計算する。KKT条件については、例えば、上記の非特許文献10に記載されているものと同様であるため、説明を省略する。KKT条件の誤差εtは、以下の10式に示すように、2つの誤差項ε1およびε2の最大値として計算される。ただし、正則化定数λが正の値の場合、誤差項ε1は、以下の11式に従って計算される。一方、正則化定数λがゼロの場合、誤差項ε1は、以下の12式に従って計算される。また、誤差項ε2は、以下の13式に従って計算される。得られた誤差εを、記憶部3の誤差37の対応する箇所に保存して、この誤差計算部26の解析処理は終了となる。<解析終了条件判定部27>上記の主変数、補助変数、およびラグランジュ乗数の更新は反復して行われ、反復する毎に、主変数、補助変数、およびラグランジュ乗数の更新をした後、解析終了条件判定部27は、KKT条件の誤差および反復回数を用いて反復を終了するべきか判定する。具体的には、KKT条件の誤差があらかじめ決められた定数ε以下であるか、又は反復回数が予め決められた回数Lに達した場合、解析終了条件を満たしたとして、反復を終了する。反復を継続する場合は反復回数tを1だけ増加させて主変数更新部23による処理に戻る。解析終了条件判定部27で反復を終了したときの未知テンソルXの要素には、時刻kにおいて観測される起点オブジェクトiと終点オブジェクトjとの間の関係の強さを表す値が格納されることになる。これにより、雑音の影響を受けたり、観測できなかったりしたオブジェクト間の関係データを復元したデータ未知テンソルXを獲得することができる。また、未知テンソルの各軸次元に注目することで、各パターン上記図1の例だと商品パターン、ユーザパターン、時間パターンの時間変化を知ることができる。さらに、任意の時刻の関係データの予測を行う場合には、解析終了条件判定部27で反復を終了した後、以下の予測値計算部4を実行する。<予測値計算部4>予測値計算部4では、まず入力部1で指定された予測時刻t’の前後の時刻TkおよびTk1を、観測時刻を要素とするベクトルTの中から検索する。これはTの要素が昇順に並んでいるので2分探索を用いて容易に実行することができる。前後の時刻に対応するインデックスkが見つかった場合、これを用いて時刻t’における予測値X’の各要素i,jを以下の14式に従って算出し、予測値X’を構成する。ここで、予測値X’はn1×n2行列で、その要素は、最適化の結果得られた主変数Xの要素と要素をT-t’:t’-Tの割合で内分したものである。t’<T1の場合は、主変数Xの要素を予測値X’の要素とする。また、t’>Tn3の場合は主変数Xの要素を予測値X’の要素とする。<出力部5>出力部5は、最適化の結果得られた主変数であるテンソルX、誤差ε、予測値X’、およびパターン基底U1U2U3を出力する。<時間関係データ解析装置の作用>次に、本実施の形態に係る時間関係データ解析装置の作用について説明する。まず、各種の値が入力部1を介して時間関係データ解析装置に入力されると、時間関係データ解析装置において、図5に示す時間関係データ解析処理ルーチンが実行される。まず、ステップ100において、入力部1により入力された、データサイズn1、n2、n3、N、M、時刻ラベルT、観測ベクトルy、観測値のインデックスI,J,K、未観測値のインデックスIc、Jc,Kc、正則化定数λ、正則化バランス定数γ1、γ2、γ3、ステップサイズη、最大繰り返し回数L,許容誤差の上限E、予測時刻t’を受け付けて、記憶部3に格納する。そして、ステップ102において、記憶部3に記憶されている主変数X、補助変数Z1,Z2,Z3、パターン基底U1、U2、U3、及びラグランジュ乗数A1、A2、A3をゼロに初期化する。ステップ104では、反復回数を示す変数lに初期値1を設定する。次のステップ106では、記憶部3に記憶されている主変数Xを更新する処理を行い、ステップ108では、記憶部3に記憶されている補助変数Z1,Z2,Z3、パターン基底U1、U2、U3、及び関係パターンV1、V2、V3を更新する処理を行う。そして、ステップ110では、記憶部3に記憶されているラグランジュ乗数A1、A2、A3を更新する処理を行い、ステップ112において、KKT条件の誤差εを計算する。次のステップ114では、上記ステップ112で計算した誤差εが、許容誤差の上限Eより大きいか否かを判定する。誤差εが、許容誤差の上限E以下である場合には、解析終了条件を満たしたと判断し、ステップ120へ移行する。一方、誤差εが、許容誤差の上限Eより大きい場合には、ステップ116において、反復回数lが、最大繰り返し回数L未満であるか否かを判定する。反復回数lが、最大繰り返し回数Lに達した場合には、解析終了条件を満たしたと判断し、ステップ120へ移行する。一方、反復回数lが、最大繰り返し回数L未満である場合には、ステップ118へ移行し、反復回数lを1インクリメントして、上記ステップ106へ戻る。ステップ120では、記憶部3に記憶されている主変数のテンソルXに基づいて、予測時刻t’の関係データを予測する。そして、ステップ122において、最終的に得られた主変数X、誤差ε、予測値X’、およびパターン基底U1U2U3を出力して、時間関係データ解析処理ルーチンを終了する。上記ステップ106の処理は、図6に示す処理ルーチンによって実現される。まず、ステップ1060において、補助変数Z1,Z2,Z3及びラグランジュ乗数A1、A2、A3に基づいて、上記4式に従って、主変数の予測値^Xを構成する。そして、ステップ1062において、上記ステップ1060で構成された主変数の予測値^Xに基づいて、上記ステップ5式に従って主変数Xを更新して、記憶部3の主変数34に格納して、処理ルーチンを終了する。上記ステップ108の処理は、図7に示す処理ルーチンによって実現される。まず、ステップ1080において、モードを示す変数kを初期値0に設定し、ステップ1082において、kを1インクリメントする。そして、ステップ1084において、主変数Xとラグランジュ乗数A1、A2、A3に基づいて、上記6式に従って、補助変数の予測値^Zkを構成する。次のステップ1086では、上記ステップ1084で構成された補助変数の予測値^Zkに対して、上記7式に従って特異値分解を計算し、ステップ1088において、上記ステップ1086での特異値分解の結果に基づいて、パターン基底Uk、関係パターンVkを更新して、記憶部3のパターン基底38、関係パターン39に格納する。そして、ステップ1090において、上記ステップ1086での特異値分解で得られる行列Skと上記ステップ1088で更新されたパターン基底Uk、関係パターンVkとに基づいて、上記8式に従って、補助変数Zkを更新して、記憶部3の補助変数35に格納する。そして、ステップ1092では、変数kが3未満であるか否かを判定し、変数kが3未満であれば、上記ステップ1082へ戻り、kを1インクメントする。一方、変数kが3に到達した場合には、各モードkについて、上記ステップ1082~1090の処理を実行したと判断し、処理ルーチンを終了する。上記ステップ110の処理は、図8に示す処理ルーチンによって実現される。まず、ステップ1100において、モードを示す変数kを初期値0に設定し、ステップ1102において、kを1インクリメントする。そして、ステップ1104では、主変数Xと、補助係数Zkと、ラグランジュ乗数Akとに基づいて、上記9式に従ってラグランジュ乗数Akを更新して、記憶部3のラグランジュ乗数36に格納する。そして、ステップ1106では、変数kが3未満であるか否かを判定し、変数kが3未満であれば、上記ステップ1102へ戻り、kを1インクメントする。一方、変数kが3に到達した場合には、各モードkについて、上記ステップ1102~1104の処理を実行したと判断し、処理ルーチンを終了する。上記ステップ112の処理は、図9に示す処理ルーチンによって実現される。まず、ステップ1120において、正則化定数λが0より大きいか否かを判定する。正則化定数λが0より大きい場合には、ステップ1122において、主変数X、観測ベクトルy、正則化定数λ、ステップサイズη、ラグランジュ乗数A1、A2、A3に基づいて、上記11式に従って、誤差項ε1を計算し、ステップ1126へ移行する。一方、正則化定数λが0である場合には、ステップ1124において、ラグランジュ乗数A1、A2、A3に基づいて、上記12式に従って、誤差項ε1を計算し、ステップ1126へ移行する。ステップ1126では、主変数X、及び補助変数Z1,Z2,Z3に基づいて、上記13式に従って、誤差項ε2を計算する。次のステップ1128では、上記ステップ112又は1124で計算した誤差項ε1と、上記ステップ1126で計算した誤差項ε2とに基づいて、上記10式に従って、誤差εを計算し、処理ルーチンを終了する。<実験例>次に、本実施の形態で提案する時間関係データの解析手法を用いた実験の結果について説明する。実験では、総務省統計局が公開している産業連関表データ逆行列係数で5年ごとに1985~2005年の間のデータを利用した。なお、もとのデータは32業種で集計されている期間と34業種で集計されている期間があるため、32業種に集計しなおしたものを使用した。ここでは1995年のデータを予測すべきデータすなわちTt’=1995とした。さらに残りの1985199020002005年のデータのうちランダムに選んだ20の要素を欠損データとした。この設定のもとでn132n232n34N4096M3277とした。また、正則化定数λ0、正則化バランス定数γ10γ20γ31ステップサイズη1許容誤差E0001最大繰り返し回数L1000とした。ただし、アルゴリズムは最大繰り返し回数に達する前に60回程度の反復で許容誤差Eを下回る誤差を達成した。さらに未観測要素に対する誤差err1を以下の15式のように定義する。err1∥true-X∥/∥X∥15ここで、Xtrueは真の逆行列係数、XはテンソルXの未観測要素をN-M次元のベクトルに並べたものを表し、∥∥は通常のユークリッドノルムを表す。また、予測点tにおける誤差err2を以下の16式のように定義する。err2∥Xtrue-X∥/∥Xtrue∥16ここで、Xtrueは真の1995年における逆行列係数を1列目、2列目、、32列目の順で縦に連結し、1024次元のベクトルにしたものであり、Xは推定された1995年における逆行列係数を同様にベクトルにしたものである。実験結果は以下の通りである。まず、得られた解のランクはr132r232r34であった。モード1とモード2に関しては正則化バランス定数γ10γ20であるからこれは妥当な結果である。モード3に関してはフルランクであるものの、得られた特異値は図10に示すように第1特異値に集中していることがわかる。従って低ランクの仮定はあてはまっていると考えられる。観測要素80のランダムな選択を100回繰り返して得られた誤差は、未観測要素に関してerr1011、補間に関する誤差err20050であった。図11に、真の1995年における逆行列係数と推定された逆行列係数を示す。ここでは縦軸にのみ、業種の名前を示しているが、横軸も縦軸と同じ業種に対応する。この図からも視覚的に推定結果が真の逆行列係数とよく一致していることがわかる。図12、13に、モード3に関する特異値分解上記図10に付随する特異ベクトルを示す。上記図12は第1成分を表わし、上記図13は、第2成分を表わしている。また、図12、13の右側にはパターン基底U3のそれぞれ1列目第1特異ベクトルおよび2列目第2特異ベクトルに対応する特異値を乗じたものを示す。なお、最適化の結果得られたU3のそれぞれの列ベクトルは4つの要素を持ち、時間点1985、1990、2000、2005年に対応し、1995年のデータ丸印で示すは1990年および2000年のデータの内分から得られる予測値である。上記図12、13の左側には、Z3の第1、第2特異値に対応する特異ベクトル、すなわち関係パターンV3の第1列、第2列をそれぞれ32×32行列に変形したものを示しており、図の右側に示すパターン基底が表現するダイナミクスに従って時間変化する逆行列係数の成分に対応する。上記図12は産業界全体での連関が2000年まで年々減少し、2000年を境にして増加に転じたことを示している。上記図13は農林水産業や鉄鋼などの伝統的な産業から他業種への影響が減少し、逆に通信放送や対事業所サービスのような非伝統的な産業の他業種へ与える影響が強まっていることを表している。以上説明したように、本発明の実施の形態に係る時間関係データ解析装置によれば、時間変化するオブジェクト間の関係の強さを示す時間関係データを解析して、目的関数を最小化するように、時間関係データを表わす未知のテンソルXを求めることによって、未観測のオブジェクト間の関係の強さを求めることができる。また、時間変化する関係データを少数のパターンに分解して、内在するパターンを発見することができる。パターンの数は補助変数更新部24の特異値分解において、γk/η以上の値を持つ特異値とそれに対応する特異ベクトルのみ計算する処理により、予めパターン数を設定せずとも自動的に最適な値が求まる。また、解析する関係データの時間発展を複数の時間的要因に分解することが可能である。また、解析する関係データを複数のパターンの重みづけとして表現するため、観測に失敗した未知数や雑音の重畳した観測データを高精度に復元予測できる。なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。例えば、上述の時間関係データ解析装置は、内部にコンピュータシステムを有しているが、コンピュータシステムは、WWWシステムを利用している場合であれば、ホームページ提供環境あるいは表示環境も含むものとする。また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。"