(57)【特許請求の範囲】
【請求項１】
複数の第１オブジェクトの各々と複数の第２オブジェクトの各々との間の関係の強さを示す関係データの時系列を表わし、かつ、テンソルで表現される時間関係データを解析して、前記時間関係データを表わす未知のテンソルＸを求める時間関係データ解析装置であって、
前記第１オブジェクトと前記第２オブジェクトとの間の関係の強さのＭ個の観測値からなる観測ベクトルｙ、正の実数である正則化定数λ、及びテンソルＸの各モードｋ（ｋ＝1、2、3）に対する相対的な正則化の強さを表わす正則化バランス定数γkを記憶した記憶手段と、
前記記憶手段から、前記観測ベクトルｙ、前記正則化定数λ、前記正則化バランス定数γk（ｋ＝1、2、3）を取得し、以下の（Ｉ）式で表される目的関数を最小化する前記テンソルＸを求める解析手段と、
を含む時間関係データ解析装置であって、
前記記憶手段は、補助変数Ｚk（ｋ＝1、2、3）、前記テンソルＸの各要素の値ｘ、ラグランジュ乗数Ａk（ｋ＝1、2、3）、第１オブジェクト方向、第２オブジェクト方向、及び時間方向の各々のパターン基底Ｕk（ｋ＝1、2、3）、前記パターン基底Ｕk（ｋ＝1、2、3）に従って変化する前記第１オブジェクト、前記第２オブジェクト、及び時間の各ペアの間の関係の強さを表わす関係パターンＶk（ｋ＝1、2、3）を更に記憶し、
前記解析手段は、
各モードｋの補助係数Ｚｋ及びラグランジュ乗数Ａkに基づいて、以下の（ＩＩ）式に従ってテンソルの予測値^Ｘを算出し、算出されたテンソルの予測値^Ｘ及び前記観測ベクトルｙに基づいて、以下の（ＩＩＩ）式に従って前記テンソルＸの各要素の値ｘを更新して前記記憶手段に記憶する主変数更新手段と、
前記テンソルＸ及び各モードｋの前記ラグランジュ乗数Ａkに基づいて、各モードｋについて、以下の（ＩＶ）式に従って、補助変数の予測値^Ｚkを算出し、算出された前記補助変数の予測値^Ｚｋに基づいて、各モードｋについて、以下の（Ｖ）式に従って、特異値分解を行って、特異値が並んだ対角行列Ｓkを求めると共に、前記パターン基底Ｕk及び前記関係パターンＶkを更新して前記記憶手段に記憶し、
前記求められた前記対角行列Ｓk、前記パターン基底Ｕk、及び前記関係パターンＶkに基づいて、以下の（ＶＩ）式に従って、前記補助変数Ｚkを更新して、前記記憶手段に記憶する補助変数更新手段と、
前記テンソルＸ及び補助変数Ｚkに基づいて、各モードｋについて、以下の（ＶＩＩ）式に従って、前記ラグランジュ乗数Ａkを更新して前記記憶手段に記憶するラグランジュ乗数更新手段と、
前記主変数更新手段による更新、前記補助変数更新手段による更新、前記ラグランジュ乗数更新手段による更新を繰り返すことで、所定の解析終了条件を満たした場合に、前記記憶手段に記憶されている前記テンソルＸを出力する解析終了条件判定手段と、を含む
時間関係データ解析装置。
【数１】
0005684055.tif 000015
ただし、ｎ1，ｎ2，ｎ3は、前記時間関係データのサイズを表す自然数であり、ｘは、前記テンソルＸの要素である。Ω（Ｘ）は、テンソルＸから前記観測ベクトルの観測値に対応するＭ個の要素の取り出す操作である。ｍａｔ1は、n1×n2n3を返すモード１行列化であり、ｍａｔ2は、n2×n3n1を返すモード２行列化であり、ｍａｔ3は、n3×n1n2を返すモード３行列化である。‖‖trは、行列の特異値の線形和を表わす。
【数２】
0005684055.tif 000016
ただし、￣Ω（Ｘ）は、前記テンソルＸから前記観測値でない要素の取り出す操作であり、ηは定数である。
【請求項２】
複数の第１オブジェクトの各々と複数の第２オブジェクトの各々との間の関係の強さを示す関係データの時系列を表わし、かつ、テンソルで表現される時間関係データを解析して、前記時間関係データを表わす未知のテンソルＸを求める時間関係データ解析装置であって、
前記第１オブジェクトと前記第２オブジェクトとの間の関係の強さのＭ個の観測値からなる観測ベクトルｙ、正の実数である正則化定数λ、及びテンソルＸの各モードｋ（ｋ＝1、2、3）に対する相対的な正則化の強さを表わす正則化バランス定数γkを記憶した記憶手段と、
前記記憶手段から、前記観測ベクトルｙ、前記正則化定数λ、前記正則化バランス定数γk（ｋ＝1、2、3）を取得し、以下の（Ｉ）式で表される目的関数を最小化する前記テンソルＸを求める解析手段と、
前記解析手段によって求められたテンソルＸに基づいて、以下の（ＶＩＩＩ）式に従って、前記第１オブジェクトｉと前記第２オブジェクトｊとの間の関係における予測対象時刻ｔ’の関係データの予測値Ｘ’（ｉ，ｊ）を計算する関係データ予測手段と、
を含む時間関係データ解析装置。
【数３】
0005684055.tif 000017
ただし、ｎ1，ｎ2，ｎ3は、前記時間関係データのサイズを表す自然数であり、ｘは、前記テンソルＸの要素である。Ω（Ｘ）は、テンソルＸから前記観測ベクトルの観測値に対応するＭ個の要素の取り出す操作である。ｍａｔ1は、n1×n2n3を返すモード１行列化であり、ｍａｔ2は、n2×n3n1を返すモード２行列化であり、ｍａｔ3は、n3×n1n2を返すモード３行列化である。‖‖trは、行列の特異値の線形和を表わす。
【数４】
0005684055.tif 000018
ただし、予測対象時刻ｔ’の前後で観測値が観測された時刻を、Ｔ（ｋ）、Ｔ（ｋ＋１）とする。
【請求項３】
第１オブジェクトと第２オブジェクトとの間の関係の強さのＭ個の観測値からなる観測ベクトルｙ、正の実数である正則化定数λ、及びテンソルＸの各モードｋ（ｋ＝1、2、3）に対する相対的な正則化の強さを表わす正則化バランス定数γkを記憶した記憶手段を備えた時間関係データ解析装置において、複数の第１オブジェクトの各々と複数の第２オブジェクトの各々との間の関係の強さを示す関係データの時系列を表わし、かつ、テンソルで表現される時間関係データを解析して、前記時間関係データを表わす未知のテンソルＸを求める時間関係データ解析方法であって、
前記解析手段によって、前記記憶手段から、前記観測ベクトルｙ、前記正則化定数λ、前記正則化バランス定数γk（ｋ＝1、2、3）を取得し、以下の（Ｉ）式で表される目的関数を最小化する前記テンソルＸを求め、
前記記憶手段は、補助変数Ｚk（ｋ＝1、2、3）、前記テンソルＸの各要素の値ｘ、ラグランジュ乗数Ａk（ｋ＝1、2、3）、第１オブジェクト方向、第２オブジェクト方向、及び時間方向の各々のパターン基底Ｕk（ｋ＝1、2、3）、前記パターン基底Ｕk（ｋ＝1、2、3）に従って変化する前記第１オブジェクト、前記第２オブジェクト、及び時間の各ペアの間の関係の強さを表わす関係パターンＶk（ｋ＝1、2、3）を更に記憶し、
前記解析手段によって前記テンソルＸを求めることは、
主変数更新手段によって、各モードｋの補助係数Ｚｋ及びラグランジュ乗数Ａkに基づいて、以下の（ＩＩ）式に従ってテンソルの予測値^Ｘを算出し、算出されたテンソルの予測値^Ｘ及び前記観測ベクトルｙに基づいて、以下の（ＩＩＩ）式に従って前記テンソルＸの各要素の値ｘを更新して前記記憶手段に記憶するステップと、
前記補助変数更新手段によって、前記テンソルＸ及び各モードｋの前記ラグランジュ乗数Ａkに基づいて、各モードｋについて、以下の（ＩＶ）式に従って、補助変数の予測値^Ｚkを算出し、算出された前記補助変数の予測値^Ｚｋに基づいて、各モードｋについて、以下の（Ｖ）式に従って、特異値分解を行って、特異値が並んだ対角行列Ｓkを求めると共に、前記パターン基底Ｕk及び前記関係パターンＶkを更新して前記記憶手段に記憶し、
前記求められた前記対角行列Ｓk、前記パターン基底Ｕk、及び前記関係パターンＶkに基づいて、以下の（ＶＩ）式に従って、前記補助変数Ｚkを更新して、前記記憶手段に記憶するステップと、
ラグランジュ乗数更新手段によって、前記テンソルＸ及び補助変数Ｚkに基づいて、各モードｋについて、以下の（ＶＩＩ）式に従って、前記ラグランジュ乗数Ａkを更新して前記記憶手段に記憶するステップと、
前記解析終了条件判定手段によって、前記主変数更新手段による更新、前記補助変数更新手段による更新、前記ラグランジュ乗数更新手段による更新を繰り返すことで、所定の解析終了条件を満たした場合に、前記記憶手段に記憶されている前記テンソルＸを出力するステップと、を含む
時間関係データ解析方法。
【数５】
0005684055.tif 000019
ただし、ｎ1，ｎ2，ｎ3は、前記時間関係データのサイズを表す自然数であり、ｘは、前記テンソルＸの要素である。Ω（Ｘ）は、テンソルＸから前記観測ベクトルの観測値に対応するＭ個の要素の取り出す操作である。ｍａｔ1は、n1×n2n3を返すモード１行列化であり、ｍａｔ2は、n2×n3n1を返すモード２行列化であり、ｍａｔ3は、n3×n1n2を返すモード３行列化である。‖‖trは、行列の特異値の線形和を表わす。
【数６】
0005684055.tif 000020
ただし、￣Ω（Ｘ）は、前記テンソルＸから前記観測値でない要素の取り出す操作であり、ηは定数である。
【請求項４】
第１オブジェクトと第２オブジェクトとの間の関係の強さのＭ個の観測値からなる観測ベクトルｙ、正の実数である正則化定数λ、及びテンソルＸの各モードｋ（ｋ＝1、2、3）に対する相対的な正則化の強さを表わす正則化バランス定数γkを記憶した記憶手段を備えた時間関係データ解析装置において、複数の第１オブジェクトの各々と複数の第２オブジェクトの各々との間の関係の強さを示す関係データの時系列を表わし、かつ、テンソルで表現される時間関係データを解析して、前記時間関係データを表わす未知のテンソルＸを求める時間関係データ解析方法であって、
前記解析手段によって、前記記憶手段から、前記観測ベクトルｙ、前記正則化定数λ、前記正則化バランス定数γk（ｋ＝1、2、3）を取得し、以下の（Ｉ）式で表される目的関数を最小化する前記テンソルＸを求め、
関係データ予測手段によって、前記解析手段によって求められたテンソルＸに基づいて、以下の（ＶＩＩＩ）式に従って、前記第１オブジェクトｉと前記第２オブジェクトｊとの間の関係における予測対象時刻ｔ’の関係データの予測値Ｘ’（ｉ，ｊ）を計算することを更に含む
時間関係データ解析方法。
【数７】
0005684055.tif 000021
ただし、ｎ1，ｎ2，ｎ3は、前記時間関係データのサイズを表す自然数であり、ｘは、前記テンソルＸの要素である。Ω（Ｘ）は、テンソルＸから前記観測ベクトルの観測値に対応するＭ個の要素の取り出す操作である。ｍａｔ1は、n1×n2n3を返すモード１行列化であり、ｍａｔ2は、n2×n3n1を返すモード２行列化であり、ｍａｔ3は、n3×n1n2を返すモード３行列化である。‖‖trは、行列の特異値の線形和を表わす。
【数８】
0005684055.tif 000022
ただし、予測対象時刻ｔ’の前後で観測値が観測された時刻を、Ｔ（ｋ）、Ｔ（ｋ＋１）とする。
【請求項５】
請求項１又は２に記載の時間関係データ解析装置を構成する各手段として、コンピュータを機能させることを特徴とするプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、時間関係データ解析装置、方法、及びプログラムに係り、特に、オブジェクト間の関係データの時系列を解析する時間関係データ解析装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
大量のデータがインターネットなどを通じて入手可能な現在、そのような大規模データセットを用いた情報解析やビジネスが盛んに行われている。しかしながら、こういった大量のデータの内容を人間が把握するのは困難である。
【０００３】
特に最近、人間関係や企業間の取引関係といった、オブジェクトとオブジェクトの間の関係を示す「関係データ」に注目が集まっている。このような関係データを用いれば、例えば友人同士の関係からより精度の高い購買アイテムの推薦を実現したりすることができる。
【０００４】
オブジェクトの間の関係データの解析手法はこれまでにも多くの方法が提案されてきている。多くの手法は、２者の間の関係を行列の形で表現できる事に着目して、その行列の構造化を行う。非特許文献１、非特許文献２は、２つのオブジェクト間の関係の有無の情報が与えられたときに、似たような接続先をもつオブジェクトのクラスタリングを行う。また、これらの関係データでは、多くの似た傾向（パターン）を持つ人間やアイテムが存在するため、関係データを行列として見たときにランクが低いという特徴があることが知られている。非特許文献６はこの性質を仮定として用いて部分的に観測された行列の未観測部分を推定する方法である。
【０００５】
これらの手法で問題となるのは、例えば友人間の関係や企業間の取引は時間に応じて変化していく、ということである。したがって、さらに高精度な関係データ解析を実現するためには、一定のオブジェクト間の関係を複数の時間で観測し、その変化を含んだ時間関係データの解析手法が必要となる。そのためには、ただの行列データの解析ではなく、関係を表す行列が時間ごとに定義されるデータの解析手法が必要となる。時間関係データの解析手法はまだ多く提案されていない。非特許文献３、非特許文献４に記載の手法はそれぞれ非特許文献１、非特許文献２を時間変化する場合に適用できるように拡張した手法である。非特許文献３に記載の手法では、状態空間モデルを利用することで、関係が連続的に変化する状況を想定したモデル化を行っている。一方、非特許文献４はＨＭＭを利用することで急激な時間変化も対応可能な時間発展モデルを提案している。また非特許文献５に記載の手法は、時間関係データをテンソルとしてモデル化し、一種の自己回帰モデルを用いて時間変化をモデル化している。
【先行技術文献】
【非特許文献】
【０００６】
【非特許文献１】K. Nowicki and T. A. B. Snijders. Estimation and prediction for stochastic blockstructures. Journal of the American Statistical Association, 96(455):1077−1087, 2001.
【非特許文献２】C. Kemp, J. B. Tenenbaum, T. L. Griffiths, T. Yamada, and N. Ueda. Learning systems of concepts with an infinite relational model. In Proceedings of the 21st National Conference on Artificial Intelligence (AAAI), 2006.
【非特許文献３】T. Yang, Y. Chi, S. Zhu, Y. Gong, and R. Jin. A Bayesian approach toward finding communities and their evolutions in dynamic social networks. In Proceedings of SIAM International Conference on Data Mining (SDM), 2009.
【非特許文献４】K. Ishiguro, T. Iwata, N. Ueda and J. Tenenbaum. Dynamic Infinite Relational Model for Time-varying Relational Data Analysis. Advances in Neural Information Processing Systems 23, 2010.
【非特許文献５】L. Xiong, X. Chen, T.-K. Huang, J. Schneider, and J. G. Carbonell. Temporal Collaborative Filtering with Bayesian Probabilistic Tensor Factorization. In Proceedings of SIAM International Conference on Data Mining (SDM), 2010.
【非特許文献６】N. Srebro, J. D. M. Rennie, and T. S. Jaakkola. Maximum-Margin Matrix Factorization. In Advances in Neural Information Processing Systems 17, 2005.
【発明の概要】
【発明が解決しようとする課題】
【０００７】
以上の手法はそれぞれ以下のような問題点を持つ。
【０００８】
まず、非特許文献１、非特許文献２、及び非特許文献６に記載の手法では、時間関係データを解析することが不可能である、という問題がある。
【０００９】
非特許文献３及び非特許文献５に記載の手法は、時間関係データをパターンに分割することが可能であるが、事前にパターンの数を決定する必要がある、という問題がある。
【００１０】
非特許文献４に記載の手法は、パターンの数を自動的に決定することができる時間関係データ解析法である。しかしながら、オブジェクト間の関係の有無だけを解析の対象としている。したがって、例えば関係性の強弱の変化といった、より正確な関係の解析には不適切である。また、この手法では観測が行われてない時刻の推定は不可能である。たとえば５年毎に調査される産業連関表のような時間関係データが与えられた場合に、他の時期（３年目など）の関係の復元や、未来における関係の予測は原理上不可能である、という問題がある。
【００１１】
本発明は、上記の課題を解決するためになされたもので、時間変化するオブジェクト間の関係データを解析して、未観測のオブジェクト間の関係の強さを求めることができる時間関係データ解析装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【００１２】
上記の目的を達成するために時間関係データ解析装置は、複数の第１オブジェクトの各々と複数の第２オブジェクトの各々との間の関係の強さを示す関係データの時系列を表わし、かつ、テンソルで表現される時間関係データを解析して、前記時間関係データを表わす未知のテンソルＸを求める時間関係データ解析装置であって、前記第１オブジェクトと前記第２オブジェクトとの間の関係の強さのＭ個の観測値からなる観測ベクトルｙ、正の実数である正則化定数λ、及びテンソルＸの各モードｋ（ｋ＝1、2、3）に対する相対的な正則化の強さを表わす正則化バランス定数γkを記憶した記憶手段と、前記記憶手段から、前記観測ベクトルｙ、前記正則化定数λ、前記正則化バランス定数γk（ｋ＝1、2、3）を取得し、以下の（Ｉ）式で表される目的関数を最小化する前記テンソルＸを求める解析手段と、を含んで構成されている。
【００１３】
【数１】
0005684055.tif 000002
【００１４】
ただし、ｎ1，ｎ2，ｎ3は、前記時間関係データのサイズを表す自然数であり、ｘは、前記テンソルＸの要素である。Ω（Ｘ）は、テンソルＸから前記観測ベクトルの観測値に対応するＭ個の要素の取り出す操作である。ｍａｔ1は、n1×n2n3を返すモード１行列化であり、ｍａｔ2は、n2×n3n1を返すモード２行列化であり、ｍａｔ3は、n3×n1n2を返すモード３行列化である。‖‖trは、行列の特異値の線形和を表わす。
【００１５】
本発明に係る時間関係データ解析方法は、第１オブジェクトと第２オブジェクトとの間の関係の強さのＭ個の観測値からなる観測ベクトルｙ、正の実数である正則化定数λ、及びテンソルＸの各モードｋ（ｋ＝1、2、3）に対する相対的な正則化の強さを表わす正則化バランス定数γkを記憶した記憶手段を備えた時間関係データ解析装置において、複数の第１オブジェクトの各々と複数の第２オブジェクトの各々との間の関係の強さを示す関係データの時系列を表わし、かつ、テンソルで表現される時間関係データを解析して、前記時間関係データを表わす未知のテンソルＸを求める時間関係データ解析方法であって、前記解析手段によって、前記記憶手段から、前記観測ベクトルｙ、前記正則化定数λ、前記正則化バランス定数γk（ｋ＝1、2、3）を取得し、上記の（Ｉ）式で表される目的関数を最小化する前記テンソルＸを求めることを特徴としている。
【００１６】
本発明によれば、前記解析手段によって、前記記憶手段から、前記観測ベクトルｙ、前記正則化定数λ、前記正則化バランス定数γk（ｋ＝1、2、3）を取得し、上記の（Ｉ）式で表される目的関数を最小化する前記テンソルＸを求める。
【００１７】
このように、時間変化するオブジェクト間の関係データを解析して、目的関数を最小化するように、時間関係データを表わす未知のテンソルＸを求めることによって、未観測のオブジェクト間の関係の強さを求めることができる。
【００１８】
本発明に係る前記記憶手段は、補助変数Ｚk（ｋ＝1、2、3）、前記テンソルＸの各要素の値ｘ、ラグランジュ乗数Ａk（ｋ＝1、2、3）、第１オブジェクト方向、第２オブジェクト方向、及び時間方向の各々のパターン基底Ｕk（ｋ＝1、2、3）、前記パターン基底Ｕk（ｋ＝1、2、3）に従って変化する前記第１オブジェクト、前記第２オブジェクト、及び時間の各ペアの間の関係の強さを表わす関係パターンＶk（ｋ＝1、2、3）を更に記憶し、前記解析手段は、各モードｋの補助係数Ｚｋ及びラグランジュ乗数Ａkに基づいて、以下の（ＩＩ）式に従ってテンソルの予測値^Ｘを算出し、算出されたテンソルの予測値^Ｘ及び前記観測ベクトルｙに基づいて、以下の（ＩＩＩ）式に従って前記テンソルＸの各要素の値ｘを更新して前記記憶手段に記憶する主変数更新手段と、前記テンソルＸ及び各モードｋの前記ラグランジュ乗数Ａkに基づいて、各モードｋについて、以下の（ＩＶ）式に従って、補助変数の予測値^Ｚkを算出し、算出された前記補助変数の予測値^Ｚｋに基づいて、各モードｋについて、以下の（Ｖ）式に従って、特異値分解を行って、特異値が並んだ対角行列Ｓkを求めると共に、前記パターン基底Ｕk及び前記関係パターンＶkを更新して前記記憶手段に記憶し、前記求められた前記対角行列Ｓk、前記パターン基底Ｕk、及び前記関係パターンＶkに基づいて、以下の（ＶＩ）式に従って、前記補助変数Ｚkを更新して、前記記憶手段に記憶する補助変数更新手段と、前記テンソルＸ及び補助変数Ｚkに基づいて、各モードｋについて、以下の（ＶＩＩ）式に従って、前記ラグランジュ乗数Ａkを更新して前記記憶手段に記憶するラグランジュ乗数更新手段と、前記主変数更新手段による更新、前記補助変数更新手段による更新、前記ラグランジュ乗数更新手段による更新を繰り返すことで、所定の解析終了条件を満たした場合に、前記記憶手段に記憶されている前記テンソルＸを出力する解析終了条件判定手段と、を含むようにすることができる。
【００１９】
【数２】
0005684055.tif 000003
【００２０】
ただし、￣Ω（Ｘ）は、前記テンソルＸから前記観測値でない要素の取り出す操作であり、ηは定数である。
【００２１】
本発明に係る時間関係データ解析装置は、前記解析手段によって求められたテンソルＸに基づいて、以下の（ＶＩＩＩ）式に従って、前記第１オブジェクトｉと前記第２オブジェクトｊとの間の関係における予測対象時刻ｔ’の関係データの予測値Ｘ’（ｉ，ｊ）を計算する関係データ予測手段を更に含むようにすることができる。
【００２２】
【数３】
0005684055.tif 000004
【００２３】
ただし、予測対象時刻ｔ’の前後で観測値が観測された時刻を、Ｔ（ｋ）、Ｔ（ｋ＋１）とする。
【００２４】
本発明に係るプログラムは、上記の時間関係データ解析装置の各手段としてコンピュータを機能させるためのプログラムである。
【発明の効果】
【００２５】
以上説明したように、本発明の時間関係データ解析装置、方法、及びプログラムによれば、時間変化するオブジェクト間の関係データを解析して、目的関数を最小化するように、時間関係データを表わす未知のテンソルＸを求めることによって、未観測のオブジェクト間の関係の強さを求めることができる、という効果が得られる。
【図面の簡単な説明】
【００２６】
【図１】時間関係データの解析に関する模式図である。
【図２】本発明の実施の形態に係る時間関係データ解析装置の構成を示す概略図である。
【図３Ａ】モード１行列化及びモード１テンソル化を説明するための図である。
【図３Ｂ】モード２行列化及びモード２テンソル化を説明するための図である。
【図３Ｃ】モード３行列化及びモード３テンソル化を説明するための図である。
【図４】特異値分解によって得られた特異値を示すグラフである。
【図５】本発明の実施の形態に係る時間関係データ解析装置における時間関係データ解析処理ルーチンの内容を示すフローチャートである。
【図６】主変数を更新する処理の内容を示すフローチャートである。
【図７】補助変数、パターン基底、及び関係パターンを更新する処理の内容を示すフローチャートである。
【図８】ラグランジュ乗数を更新する処理の内容を示すフローチャートである。
【図９】ＫＫＴ条件の誤差を計算する処理の内容を示すフローチャートである。
【図１０】実験で得られた特異値を示すグラフである。
【図１１】実験で得られた真の1995年における逆行列係数と推定された逆行列係数を示す図である。
【図１２】実験で得られたモード３に関する特異値分解に付随する第１成分の特異ベクトル及び特異値を示す図である。
【図１３】実験で得られたモード３に関する特異値分解に付随する第２成分の特異ベクトル及び特異値を示す図である。
【発明を実施するための形態】
【００２７】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００２８】
＜概要＞
まず、本発明で提案する時間関係データ解析装置で用いる時間関係データについて説明する。図１は、提案する時間関係データ解析装置のコンセプトを模式的に表した図である。
【００２９】
図１において、左側にはユーザ（購買者）と商品の間の関係（買ったか否か）が時間軸に沿ってテンソルとして表現されている。黒い立方体は観測されている関係を表し、テンソルのそれ以外の要素は未知であるとする。図１の右側には左側の部分観測テンソルの低ランク分解を示している。中心を取り囲む３つの薄い直方体はそれぞれ商品、ユーザ、および時間変化を少数のパターンで表現したものであり、中心の小さい直方体はパターン間の相互作用を表現している。このように分解することにより、左側のテンソルの未知部分の推定も可能となる。
【００３０】
＜システム構成＞
次に、観測データを入力として、時間関係データを表わす未知のテンソルを求めて出力する時間関係データ解析装置に、本発明を適用した場合を例にして、本発明の実施の形態を説明する。
【００３１】
図２に示すように、本実施の形態に係る時間関係データ解析装置は、時間関係データの観測値の入力を受け付ける入力部１と、観測された時間関係データの解析を行う解析部２と、変数や定数を記憶する記憶部３と、予測値計算部４と、解析結果を出力する出力部５と、を備えている。
【００３２】
入力部１は、既知のビデオカメラ、マイクロフォン、記憶装置などの入力器により実現される。また、出力部５は、ディスプレイ、プリンタ、磁気ディスクなどで実装される。
【００３３】
解析部２、記憶部３、及び予測値計算部４は、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）と、ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）と、後述する時間関係データ解析処理ルーチンを実行するためのプログラムを記憶したＲＯＭ（Ｒｅａｄ  Ｏｎｌｙ Ｍｅｍｏｒｙ）とを備えたコンピュータで構成されている。
【００３４】
解析部２は、初期設定部２１、変数計算部２２、及び解析終了条件判定部２７を備えている。変数計算部２２は、主変数更新部２３、補助変数更新部２４、ラグランジュ乗数更新部２５、及び誤差計算部２６を備えている。
【００３５】
初期設定部２１は、入力された観測データから、記憶部３の各要素に値を設定する。
【００３６】
変数計算部２２は、記憶部３に保存された情報を用いて各変数の更新を行う。
【００３７】
主変数更新部２３は、記憶部３に保存された情報を用いて主変数の更新を行う。
【００３８】
補助変数更新部２４では、記憶部３に保存された情報を用いて、補助変数の更新とパターン基底の更新と関係パターンの更新とを行う。
【００３９】
ラグランジュ乗数更新部２５は、記憶部３に保存された情報を用いてラグランジュ乗数の更新を行う。
【００４０】
誤差計算部２６では、記憶部３に保存された情報を用いて誤差の計算を行う。
【００４１】
解析終了条件判定部２７では、変数計算部２２を用いて計算され記憶部３に記憶された数値をもちいて解析を終了してよいかどうかを判定する。
【００４２】
記憶部３は、定数３１、観測値３２、インデックス集合３３、主変数３４、補助変数３５、ラグランジュ乗数３６、誤差３７、パターン基底３８、関係パターン３９から構成される。
【００４３】
予測値計算部４では、指定された予測時刻t’における関係データの予測を行う。
【００４４】
＜データ構造＞
ｎ１、ｎ２、ｎ３は想定する時間関係データのサイズを表す自然数である。関係データは２種類のオブジェクトの間の接続関係を表現する。本実施の形態では、関係を構成する２種類のオブジェクトの一方を起点オブジェクト、他方を終点オブジェクトと呼ぶ。例えば、図１の例では商品が起点オブジェクト、ユーザが終点オブジェクトとなる。また、本実施の形態では、ｎ１、ｎ２はそれぞれ関係の起点と終点のオブジェクト数とし、ｎ３は観測が発生した時間の総ステップ数とする。また、可能なすべての要素数をＮとし、Ｎ＝ｎ１×ｎ２×ｎ３とする。また、起点となるオブジェクトの軸をテンソルの１番目の次元、終点となるオブジェクトの軸をテンソルの２番目の次元、時間の軸をテンソルの３番目の次元と呼ぶ。
【００４５】
Ｔはn3個の観測された時間点の時刻を格納するベクトルである。Ｔの要素は昇順に並んでいるものとする。
【００４６】
観測ベクトルｙは観測された関係データであり長さＭの縦ベクトルで表現される。Ｍはｙの要素数であり、ｙｍはｙのｍ番目（１≦ｍ≦Ｍ）の観測値を表す。なお、本実施の形態では観測値ｙｍはすべての実数を許容する。したがって、この値はある関係の強さを表すと理解できる。一方、上記の非特許文献３や非特許文献４では観測される関係データは１あるいは０の２値しか許容しないため、関係の有無しか表現できず、関係の強さを考慮した解析が不可能である。
【００４７】
Ｉ，Ｊ，Ｋはそれぞれ長さＭのインデックス集合３３で、それぞれのｍ番目の要素を並べた(im,jm,km)が未知テンソルのｍ番目の観測された要素を指し示す。すなわちｍ番目の観測値はim番目の起点オブジェクトとjm番目の終点オブジェクトの時刻kmでの関係である。
【００４８】
また同様にIc,Jc,Kcは、それぞれ長さ（Ｎ−Ｍ）のインデックス集合３３で、それぞれのｍ番目の要素の並べたもの(icm,jcm,kcm)は未知テンソルのｍ番目の未観測要素を指し示す。
【００４９】
λは正則化定数とよばれる正の実数である。分解・解析計算における正則化項の効果の大きさを表す。
【００５０】
Ｍ次元の観測ベクトルｙは、観測値３２に格納される。
【００５１】
主変数３４にはｎ１×ｎ２×ｎ３の要素からなる未知テンソルＸが格納される。初期値としてＸの各要素には適当な値を入力する。一般にはすべての値を０とすれば良い。未知テンソルＸは、上記図１に示すテンソルに相当するものであり、最終的には、要素(i,j,k)に、時刻kにおいて観測される起点オブジェクトiと終点オブジェクトjとの間の関係の強さを示す値が格納（上書き）されることになる。
【００５２】
以下に観測要素の取り出しΩ、未観測要素の取り出し￣Ω、観測要素の復元ΩT、モードk行列化matk、モードkテンソル化tensorkと呼ぶ4つの操作を定義する。これらの操作は以下の数式では行列演算として表現されるが、実際には行列として保持する必要はない。まず、観測要素の取り出しΩは、引数としてｎ１×ｎ２×ｎ３の要素からなるテンソルを取り、Ｍ次元ベクトルを返すもので、戻り値のｍ番目の要素は、引数であるテンソルの第(im,jm,km)要素である。同様に、未観測要素の取り出し￣Ωは、引数としてｎ１×ｎ２×ｎ３の要素からなるテンソルを取り、（Ｎ−Ｍ）次元ベクトルを返すもので、戻り値のｍ番目の要素は、引数であるテンソルの第(icm,jcm,kcm)要素である。観測要素の復元ΩTは観測要素の取り出しΩの逆の操作として定義され、Ｍ次元ベクトルを引数としてｎ１×ｎ２×ｎ３の要素からなるテンソルを返すもので、引数のｍ番目の要素は、戻り値であるテンソルの第(im,jm,km)要素に格納される。観測値に対応しない戻り値の要素はすべてゼロである。モードk行列化mat-kはｋ＝１,２,３の３種存在し、それぞれｎ１×ｎ２×ｎ３の要素からなるテンソルを引数として、n1×n2n3の行列、n2×n3n1の行列、n3×n1n2の行列を返すものである。具体的にはモード１行列化mat1は図３Ａに示すように与えられたテンソルの第1モード、すなわち起点オブジェクトが縦に連続するようにベクトルを取り出し、横にならべて行列にしたものを生成する操作である。また、第１モードテンソル化tensor-1はその逆の操作である。同様に、第２モード行列化mat-2は図３Ｂに示すように終点オブジェクトが縦に連続するようにベクトルを取り出し、横にならべて行列にしたものを生成する操作であり、第２モードテンソル化tensor2はその逆の操作である。第３モード行列化mat3は図３Ｃのように時間点が縦に連続するようにベクトルを取り出し、横に並べて行列にしたものを生成する操作であり、第３モードテンソル化tensor-3はその逆の操作である。
【００５３】
＜入力部１＞
入力部１では、解析にあたり事前に必要となる数値データの入力を受け付ける。入力されるデータは以下の通りである。
【００５４】
ｎ１、ｎ２、ｎ３、Ｎ、Ｍ、Ｉ，Ｊ，Ｋ，Ｉｃ、Ｊｃ，Ｋｃ、ｙ、λ、γ1、γ2、γ3、η、Ｌ，Ｅ、t’
【００５５】
ｎ１、ｎ２、ｎ３は想定する時間関係データのサイズを表す自然数である。可能なすべての要素数をＮとし、Ｎ＝ｎ１×ｎ２×ｎ３とする。
【００５６】
ｙは観測された関係データであり長さMの縦ベクトルで表現される観測ベクトルである。Ｍはｙの要素数であり、ｙｍはｍ番目の観測値を表す。
【００５７】
Ｉ，Ｊ，Ｋは観測値のインデックスであり、Ｍ個の観測値ｙと未知テンソルＸの関係を表す。同様に（Ｎ−Ｍ）個の未観測値のインデックスであるＩｃ、Ｊｃ，Ｋｃも入力される。
【００５８】
λは正則化定数とよばれる正の実数である。分解・解析計算における正則化項の効果の大きさを表す。γ1、γ2、γ3、は未知テンソルＸの各モードｋに対する相対的な正則化の強さを表現する正の定数である。ここではこれらを正則化バランス定数と呼ぶ。また、正則化定数λをある定数倍して正則化バランス定数γ1、γ2、γ3をその定数で割っても得られる解は同じであるため、規格化のためγ1＋γ2＋γ3＝１と仮定する。
【００５９】
ηはステップサイズであり、正の実数である。分解・解析計算においてラグランジュ乗数Ａ１、Ａ２、Ａ３の更新量を決める定数であり、主変数と補助変数の乖離に対するペナルティーの強さを表す。
【００６０】
Ｌは変数計算部２２における、最大繰り返し回数を表す自然数である。繰り返し回数はｌで表す。
【００６１】
Ｅは解析終了条件判定部２７における、許容される誤差の上限を表す定数である。
【００６２】
t’は予測したい時刻のインデックスを表す実数である。なお、予測が不要な場合にはt’は入力しなくても良い。
【００６３】
＜解析部２＞
解析部２は、入力された観測値や定数をもとにして、主変数である未知テンソルＸを求める。解析部２は、以下の（１）式で示される目的関数を、近似的に最小化することにより主変数Ｘを決定する。
【００６４】
【数４】
0005684055.tif 000005
【００６５】
ここで、上記（１）式の第１項は損失項であり、観測ベクトルyと、未知テンソルＸから観測値に対応する要素を取り出して得られるＭ次元ベクトルΩ(Ｘ)との差の二乗和として定義される。ただし、観測値に欠損値がある場合、第１項だけでは未知テンソルＸを一意に定めることができない。第２項は正則化項と呼ばれ、未知テンソルのランクが低いという仮定を表現する。ただし、ここで||Ｚk||trはトレースノルムと呼ばれ、行列Ｚkの特異値の線形和として定義される。ここで未知テンソルが上記図１のように少ない数のパターンとそれらの間の相互作用に分解されることは、行列化操作mat-kによって得られる行列が低ランクであることに対応することに注意する（例えば非特許文献８（T. G. Kolda and B. W. Bader. Tensor decompositions and applications. SIAM Review, 51(3):455−500, 2009.）を参照）。なぜトレースノルム正則化が低ランク性の仮定に対応しているかを説明するために、全要素が観測されている場合で簡単のために第１モードに対応する正則化項のみ考える（つまりγ1＝１,γ2＝０,γ3＝０）。この時、上記の最適化問題は、以下の（２）式のように書き直すことができる。
【００６６】
【数５】
0005684055.tif 000006
【００６７】
ここで、Ｘ、Ｙはそれぞれ、未知テンソルＸ、および観測ベクトルｙのモード１行列化によって得られる行列とする。上の最適化問題の解は、行列Ｙの特異値分解をＹ＝ＵＳＶTとして、解析的に以下の（３）式のように書けることが知られている（例えば、非特許文献９（R. Tomioka, T. Suzuki, and M. Sugiyama. Augmented lagrangian methods for learning, selecting, and combining features. In S. Sra, S. Nowozin, and S. J. Wright, editors, Optimization for Machine Learning. MIT Press, 2011.）を参照）。
【００６８】
【数６】
0005684055.tif 000007
【００６９】
ここでＳは行列Ｙの特異値が並んだ対角行列である。Ｕ、Ｖは各列が特異値に対応する特異ベクトルであるような行列である。また、ｍａｘの操作は行列の要素ごとに適用されるものとする。上記（３）式からこの最適化問題を解くと、図４に示すように、与えられた行列Ｙの特異値のうち正則化定数λ以下の特異値は雑音として切り捨てられ、雑音より十分大きい信号成分を取り出すことができることがわかる。また、上記の操作は観測値Ｙの両側から直交行列をかけても不変であり、これらの操作に対して頑健になっている。
【００７０】
上記最適化問題は、たとえば非特許文献７（S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein. Distributed optimization and statistical learning via the alternating direction method of multipliers. http://www.stanford.edu/~boyd/papers/distr_opt_stat_learning_admm.html）のようなAlternating Direction Method of Multipliers (ADMM)法を用いて解くことができる。さらに、上記最適化問題は凸最適化問題（凸最適化に関する一般論は非特許文献１０（ S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press. 2004.）を参照）なので、ADMM法が収束すれば、それが上記の最適化問題の大域最適解である。ADMM法の収束は下で詳述するようにKKT条件の誤差を用いて判定する。これにより、テンソルを幾つのパターンに分解するか（次元）を指定する必要がなく、最適化を通してこれらを自動的に決定することが可能となる。
【００７１】
＜初期設定部２１＞
初期設定部２１では、入力部１から入力された情報をもとに、記憶部３に記憶されている各要素の初期設定を行う。
【００７２】
まず、入力部１よりｎ１、ｎ２、ｎ３、Ｎ、Ｍ、ｙ、λ、γ1、γ2、γ3、η、Ｌ，Ｅ、t’を定数３１に記憶する。また繰り返し回数ｌを１に初期化する。
【００７３】
Ｍ次元の観測ベクトルｙを観測値３２に、観測値のインデックスＩ，Ｊ，Ｋおよび未観測値のインデックスＩｃ、Ｊｃ，Ｋｃをインデックス集合３３に保存する。
【００７４】
主変数３４にはｎ１×ｎ２×ｎ３の要素からなる未知テンソルＸを記憶する。初期値としてＸの各要素には適当な値を入力する。一般にはすべての値を０とすれば良い。
【００７５】
-  補助変数３５にはＺ１，Ｚ２，Ｚ３が記憶され、それぞれ大きさｎ１×ｎ２ｎ３、ｎ２×ｎ３ｎ１、ｎ３×ｎ１ｎ２の行列として表現される。そのベクトル化をそれぞれｚ１，ｚ２，ｚ３と書く。補助変数は変数計算部２２において主変数ｘに対する制約を緩和したものである。具体的にはＺkは最適化が進むに従ってテンソルＸのモードk行列化と一致していく量である。初期値として補助変数Ｚ１，Ｚ２，Ｚ３には適当な値を入力する。一般にはすべての値を０とすれば良い。
【００７６】
ラグランジュ乗数３６には、ラグランジュ乗数Ａ1、Ａ2、Ａ3が記憶され、ラグランジュ乗数Ａ1、Ａ2、Ａ3は補助変数Ｚ1、Ｚ２、Ｚ３と同じ次元を持つ行列である。ラグランジュ乗数は主変数Ｘと補助変数Ｚ１、Ｚ２、Ｚ３を一致させる役割を持つ。初期値としてラグランジュ乗数Ａ1、Ａ2、Ａ3には適当な値を入力する。一般にはすべての値を０とすれば良い。
【００７７】
誤差３７には、ε１、ε２、・・・、εＴが記憶され、ε１、ε２、・・・、εＴは誤差計算部２６において計算される誤差を表す。
【００７８】
パターン基底３８には、パターン基底Ｕ１、Ｕ２、Ｕ３が記憶される。パターン基底Ｕ１、Ｕ２、Ｕ３はそれぞれｎ１×ｒ１、ｎ２×ｒ２、ｎ３×ｒ３の行列であり、変数計算部２２で得られた起点オブジェクト方向、終点オブジェクト方向、および時間方向のパターンを表現する。関係パターン３９には、関係パターンＶ１、Ｖ２、Ｖ３が記憶される。関係パターンＶ１、Ｖ２、Ｖ３はそれぞれｎ２ｎ３×ｒ１、ｎ３ｎ１×ｒ２、ｎ１ｎ２×ｒ３の行列であり、パターン基底Ｕｋに従って変化する、終点オブジェクトおよび時間の間の関係の強さ、起点オブジェクトおよび時間の間の関係の強さ、及び起点オブジェクトおよび終点オブジェクトの間の関係の強さを表現する。ｒ１、ｒ２、ｒ３はそれぞれパターンの要素数を表す。この値はテンソルの分解解析が進むにつれて自動的に決定される。初期値として、パターン基底Ｕ１、Ｕ２、Ｕ３、関係パターンＶ１、Ｖ２、Ｖ３、パターンの要素数ｒ１、ｒ２、ｒ３には適当な値を入力する。一般にはパターン基底はすべての値を０とすれば良い。ｒ１、ｒ２、ｒ３は適当な自然数を選択する。
【００７９】
＜変数計算部２２＞
変数計算部２２は、記憶部３に記憶された各要素を利用して、ＡＤＭＭ法に基づいた主変数の推定を行う。実際の計算では、主変数更新部２３、補助変数更新部２４、ラグランジュ乗数更新部２５、及び誤差計算部２６の各処理をこの順番で行う解析処理を繰り返し行う。
【００８０】
繰り返しごとに解析終了条件判定部２７で解析の終了を判定し、解析終了とならない場合は繰り返し回数ｔを１増加させ、再度解析処理を行う。
【００８１】
＜主変数更新部２３＞
主変数更新部２３は、以下の（４）式および（５）を用いて主変数Ｘの更新を行う。
【００８２】
【数７】
0005684055.tif 000008
【００８３】
まず、主変数更新部２３は、上記の（４）式を用いて主変数の予測値^Ｘを構成する。予測値^Ｘは、主変数と同じサイズのテンソルであり、各モードｋについて補助変数Ｚkからラグランジュ乗数Ａkを差し引き、テンソル化したものの平均である。
【００８４】
次に、主変数更新部２３は、構成した主変数の予測値^Ｘを用いて、上記（５）式のように主変数Ｘを更新する。上記（５）式では、観測値の得られている要素については、新しい主変数の値が、その予測値と観測値を３λη：１の割合で内分したものに等しくなるように更新している。一方、観測値の得られていない要素については、新しい主変数の値がその予測値そのものとなるように更新している。得られた主変数Ｘを記憶部３の主変数３４に保存して、この主変数更新部２３の処理は終了である。
【００８５】
＜補助変数更新部２４＞
補助変数更新部２４は、以下の（６）式、（７）式、および（８）式を用いて、補助変数Ｚ１，Ｚ２，Ｚ３、パターン基底Ｕ１、Ｕ２、Ｕ３、及び関係パターンＶ１，Ｖ２，Ｖ３の更新を行う。
【００８６】
【数８】
0005684055.tif 000009
【００８７】
まず、上記（６）式のように、記憶部３に保存されている量を用いて補助変数の予測値＾Ｚkを構成する。補助変数の予測値＾Ｚkは各モードｋについて主変数Ｘを行列化したものにラグランジュ乗数Ａkを足したものである。
【００８８】
次に、このように構成した補助変数の予測値＾Ｚkに対し、上記（７）式に従って特異値分解を行う。特異値分解を行うと、ｋ番目の予測値＾Ｚkに対し、Ｕｋ、Ｓｋ、Ｖｋという行列が得られる。このうちＳｋは特異値が並んだ対角行列である。Ｕｋ、Ｖｋはそれぞれの列が特異値に対応する特異ベクトルである。ただし、特異値及び特異ベクトルとしては、次の（８）式の操作のために、γｋ／η以上の値を持つ特異値とそれに対応する特異ベクトルのみ計算すればよい。
【００８９】
最後に、上記（７）式で得られた特異値と特異ベクトルを用いて、上記（８）式に従って、補助変数Ｚkを更新する。また、新しいパターン基底Ｕｋは、max関数によって０となる列を除いたものである。すなわち、パターンの要素数ｒｋは特異値分解の結果によって繰り返しのたびに変動しうる。
【００９０】
得られた補助変数Ｚkをベクトル化して、記憶部３の補助変数３５に保存し、パターン基底Ｕkを記憶部３のパターン基底３８に保存し、さらに関係パターンＶkを記憶部３の関係パターン３９に保存して、この補助変数更新部２４の解析処理は終了となる。
【００９１】
＜ラグランジュ乗数更新部２５＞
ラグランジュ乗数更新部２５は、以下の（９）式を用いてラグランジュ乗数Ａ１、Ａ２、Ａ３の更新を行う。
【００９２】
【数９】
0005684055.tif 000010
【００９３】
上記（９）式では、ラグランジュ乗数Ａ１、Ａ２、Ａ３を、主変数Ｘを行列化したものとそれぞれの補助変数Ｚ１、Ｚ２、Ｚ３の乖離の量だけ更新する。得られたラグランジュ乗数を記憶部３のラグランジュ乗数３６に保存して、このラグランジュ乗数更新部２５の解析処理は終了である。
【００９４】
＜誤差計算部２６＞
誤差計算部２６では、ＫＫＴ条件誤差を計算する。ＫＫＴ条件については、例えば、上記の非特許文献１０に記載されているものと同様であるため、説明を省略する。
【００９５】
ＫＫＴ条件の誤差εtは、以下の（１０）式に示すように、２つの誤差項ε1およびε2の最大値として計算される。
【００９６】
【数１０】
0005684055.tif 000011
【００９７】
ただし、正則化定数λが正の値の場合、誤差項ε1は、以下の（１１）式に従って計算される。一方、正則化定数λがゼロの場合、誤差項ε1は、以下の（１２）式に従って計算される。
【００９８】
【数１１】
0005684055.tif 000012
【００９９】
また、誤差項ε2は、以下の（１３）式に従って計算される。
【０１００】
【数１２】
0005684055.tif 000013
【０１０１】
得られた誤差εを、記憶部３の誤差３７の対応する箇所に保存して、この誤差計算部２６の解析処理は終了となる。
【０１０２】
＜解析終了条件判定部２７＞
上記の主変数、補助変数、およびラグランジュ乗数の更新は反復して行われ、反復する毎に、主変数、補助変数、およびラグランジュ乗数の更新をした後、解析終了条件判定部２７は、ＫＫＴ条件の誤差および反復回数を用いて反復を終了するべきか判定する。具体的には、ＫＫＴ条件の誤差があらかじめ決められた定数ε以下であるか、又は反復回数が予め決められた回数Ｌに達した場合、解析終了条件を満たしたとして、反復を終了する。反復を継続する場合は反復回数tを１だけ増加させて主変数更新部２３による処理に戻る。
【０１０３】
解析終了条件判定部２７で反復を終了したときの未知テンソルＸの要素(i,j,k)には、時刻kにおいて観測される起点オブジェクトiと終点オブジェクトjとの間の関係の強さを表す値が格納されることになる。これにより、雑音の影響を受けたり、観測できなかったりしたオブジェクト間の関係データを復元したデータ（未知テンソルＸ）を獲得することができる。また、未知テンソルの各軸（次元）に注目することで、各パターン（上記図１の例だと商品パターン、ユーザパターン、時間パターン）の時間変化を知ることができる。
【０１０４】
さらに、任意の時刻の関係データの予測を行う場合には、解析終了条件判定部２７で反復を終了した後、以下の予測値計算部４を実行する。
【０１０５】
＜予測値計算部４＞
予測値計算部４では、まず入力部１で指定された予測時刻t’の前後の時刻Ｔ（ｋ）およびＴ（ｋ＋１）を、観測時刻を要素とするベクトルＴの中から検索する。これはＴの要素が昇順に並んでいるので２分探索を用いて容易に実行することができる。
【０１０６】
前後の時刻に対応するインデックスｋが見つかった場合、これを用いて時刻ｔ’における予測値Ｘ’の各要素（ｉ，ｊ）を以下の（１４）式に従って算出し、予測値Ｘ’を構成する。
【０１０７】
【数１３】
0005684055.tif 000014
【０１０８】
ここで、予測値Ｘ’はｎ1×ｎ2行列で、その要素(i,j)は、最適化の結果得られた主変数Ｘの要素(i,j,k)と要素(i,j,k+1)をＴ(k+1)−ｔ’：ｔ’−Ｔ(k)の割合で内分したものである。
【０１０９】
ｔ’＜Ｔ（１）の場合は、主変数Ｘの要素(i,j,1)を予測値Ｘ’の要素(i,j)とする。また、ｔ’＞Ｔ（ｎ3）の場合は主変数Ｘの要素(i,j,n3)を予測値Ｘ’の要素(i,j)とする。
【０１１０】
＜出力部５＞
出力部５は、最適化の結果得られた主変数であるテンソルＸ、誤差ε、予測値Ｘ’、およびパターン基底Ｕ1,Ｕ2,Ｕ3を出力する。
【０１１１】
＜時間関係データ解析装置の作用＞
次に、本実施の形態に係る時間関係データ解析装置の作用について説明する。まず、各種の値が入力部１を介して時間関係データ解析装置に入力されると、時間関係データ解析装置において、図５に示す時間関係データ解析処理ルーチンが実行される。
【０１１２】
まず、ステップ１００において、入力部１により入力された、データサイズｎ１、ｎ２、ｎ３、Ｎ、Ｍ、時刻ラベルＴ、観測ベクトルｙ、観測値のインデックスＩ，Ｊ，Ｋ、未観測値のインデックスＩｃ、Ｊｃ，Ｋｃ、正則化定数λ、正則化バランス定数γ1、γ2、γ3、ステップサイズη、最大繰り返し回数Ｌ，許容誤差の上限Ｅ、予測時刻t’を受け付けて、記憶部３に格納する。
【０１１３】
そして、ステップ１０２において、記憶部３に記憶されている主変数Ｘ、補助変数Ｚ１，Ｚ２，Ｚ３、パターン基底Ｕ１、Ｕ２、Ｕ３、及びラグランジュ乗数Ａ1、Ａ2、Ａ3をゼロに初期化する。ステップ１０４では、反復回数を示す変数ｌに初期値１を設定する。
【０１１４】
次のステップ１０６では、記憶部３に記憶されている主変数Ｘを更新する処理を行い、ステップ１０８では、記憶部３に記憶されている補助変数Ｚ１，Ｚ２，Ｚ３、パターン基底Ｕ１、Ｕ２、Ｕ３、及び関係パターンＶ1、Ｖ2、Ｖ3を更新する処理を行う。
【０１１５】
そして、ステップ１１０では、記憶部３に記憶されているラグランジュ乗数Ａ1、Ａ2、Ａ3を更新する処理を行い、ステップ１１２において、ＫＫＴ条件の誤差εを計算する。
【０１１６】
次のステップ１１４では、上記ステップ１１２で計算した誤差εが、許容誤差の上限Ｅより大きいか否かを判定する。誤差εが、許容誤差の上限Ｅ以下である場合には、解析終了条件を満たしたと判断し、ステップ１２０へ移行する。一方、誤差εが、許容誤差の上限Ｅより大きい場合には、ステップ１１６において、反復回数ｌが、最大繰り返し回数Ｌ未満であるか否かを判定する。反復回数ｌが、最大繰り返し回数Ｌに達した場合には、解析終了条件を満たしたと判断し、ステップ１２０へ移行する。一方、反復回数ｌが、最大繰り返し回数Ｌ未満である場合には、ステップ１１８へ移行し、反復回数ｌを１インクリメントして、上記ステップ１０６へ戻る。
【０１１７】
ステップ１２０では、記憶部３に記憶されている主変数のテンソルＸに基づいて、予測時刻ｔ’の関係データを予測する。そして、ステップ１２２において、最終的に得られた主変数Ｘ、誤差ε、予測値Ｘ’、およびパターン基底Ｕ1,Ｕ2,Ｕ3を出力して、時間関係データ解析処理ルーチンを終了する。
【０１１８】
上記ステップ１０６の処理は、図６に示す処理ルーチンによって実現される。
【０１１９】
まず、ステップ１０６０において、補助変数Ｚ１，Ｚ２，Ｚ３及びラグランジュ乗数Ａ1、Ａ2、Ａ3に基づいて、上記（４）式に従って、主変数の予測値＾Ｘを構成する。そして、ステップ１０６２において、上記ステップ１０６０で構成された主変数の予測値＾Ｘに基づいて、上記ステップ（５）式に従って主変数Ｘを更新して、記憶部３の主変数３４に格納して、処理ルーチンを終了する。
【０１２０】
上記ステップ１０８の処理は、図７に示す処理ルーチンによって実現される。
【０１２１】
まず、ステップ１０８０において、モードを示す変数ｋを初期値０に設定し、ステップ１０８２において、ｋを１インクリメントする。
【０１２２】
そして、ステップ１０８４において、主変数Ｘとラグランジュ乗数Ａ1、Ａ2、Ａ3に基づいて、上記（６）式に従って、補助変数の予測値＾Ｚkを構成する。次のステップ１０８６では、上記ステップ１０８４で構成された補助変数の予測値＾Ｚkに対して、上記（７）式に従って特異値分解を計算し、ステップ１０８８において、上記ステップ１０８６での特異値分解の結果に基づいて、パターン基底Ｕｋ、関係パターンＶkを更新して、記憶部３のパターン基底３８、関係パターン３９に格納する。
【０１２３】
そして、ステップ１０９０において、上記ステップ１０８６での特異値分解で得られる行列Ｓｋと上記ステップ１０８８で更新されたパターン基底Ｕｋ、関係パターンVｋとに基づいて、上記（８）式に従って、補助変数Ｚkを更新して、記憶部３の補助変数３５に格納する。
【０１２４】
そして、ステップ１０９２では、変数ｋが３未満であるか否かを判定し、変数ｋが３未満であれば、上記ステップ１０８２へ戻り、ｋを１インクメントする。一方、変数ｋが３に到達した場合には、各モードｋについて、上記ステップ１０８２〜１０９０の処理を実行したと判断し、処理ルーチンを終了する。
【０１２５】
上記ステップ１１０の処理は、図８に示す処理ルーチンによって実現される。
【０１２６】
まず、ステップ１１００において、モードを示す変数ｋを初期値０に設定し、ステップ１１０２において、ｋを１インクリメントする。
【０１２７】
そして、ステップ１１０４では、主変数Ｘと、補助係数Ｚkと、ラグランジュ乗数Ａｋとに基づいて、上記（９）式に従ってラグランジュ乗数Ａｋを更新して、記憶部３のラグランジュ乗数３６に格納する。
【０１２８】
そして、ステップ１１０６では、変数ｋが３未満であるか否かを判定し、変数ｋが３未満であれば、上記ステップ１１０２へ戻り、ｋを１インクメントする。一方、変数ｋが３に到達した場合には、各モードｋについて、上記ステップ１１０２〜１１０４の処理を実行したと判断し、処理ルーチンを終了する。
【０１２９】
上記ステップ１１２の処理は、図９に示す処理ルーチンによって実現される。
【０１３０】
まず、ステップ１１２０において、正則化定数λが０より大きいか否かを判定する。正則化定数λが０より大きい場合には、ステップ１１２２において、主変数Ｘ、観測ベクトルｙ、正則化定数λ、ステップサイズη、ラグランジュ乗数Ａ1、Ａ2、Ａ3に基づいて、上記（１１）式に従って、誤差項ε1を計算し、ステップ１１２６へ移行する。一方、正則化定数λが０である場合には、ステップ１１２４において、ラグランジュ乗数Ａ1、Ａ2、Ａ3に基づいて、上記（１２）式に従って、誤差項ε1を計算し、ステップ１１２６へ移行する。
【０１３１】
ステップ１１２６では、主変数Ｘ、及び補助変数Ｚ１，Ｚ２，Ｚ３に基づいて、上記（１３）式に従って、誤差項ε2を計算する。次のステップ１１２８では、上記ステップ１１２又は１１２４で計算した誤差項ε1と、上記ステップ１１２６で計算した誤差項ε2とに基づいて、上記（１０）式に従って、誤差εを計算し、処理ルーチンを終了する。
【０１３２】
＜実験例＞
次に、本実施の形態で提案する時間関係データの解析手法を用いた実験の結果について説明する。実験では、総務省統計局が公開している産業連関表データ（逆行列係数）で５年ごとに１９８５〜２００５年の間のデータを利用した。なお、もとのデータは32業種で集計されている期間と34業種で集計されている期間があるため、32業種に集計しなおしたものを使用した。ここでは1995年のデータを予測すべきデータ（すなわちＴ（ｔ’）＝１９９５）とした。さらに残りの1985, 1990, 2000, 2005年のデータのうちランダムに選んだ20%の要素を欠損データとした。この設定のもとでn1=32, n2=32, n3=4, N=4096, M=3277とした。また、正則化定数λ=0、正則化バランス定数γ1=0,γ2=0,γ3=1, ステップサイズη=1, 許容誤差Ｅ=0.001, 最大繰り返し回数L=1000とした。ただし、アルゴリズムは最大繰り返し回数に達する前に60回程度の反復で許容誤差Eを下回る誤差を達成した。さらに未観測要素に対する誤差err1を以下の（１５）式のように定義する。
【０１３３】
err1=‖true(Ic,Jc,Kc) − X(Ic,Jc,Kc) ‖／‖X(Ic,Jc,Kc) ‖      ・・・（１５）
【０１３４】
ここで、Xtrueは真の逆行列係数、X(Ic,Jc,Kc) はテンソルXの未観測要素を（N−M）次元のベクトルに並べたものを表し、‖ ‖は通常のユークリッドノルムを表す。また、予測点t'における誤差err2を以下の（１６）式のように定義する。
【０１３５】
err2=‖Xtrue(:,:,t')−X'(:)‖／‖Xtrue(:,:,t') ‖      ・・・（１６）
【０１３６】
ここで、Xtrue(:,:,t')は真の1995年における逆行列係数を1列目、2列目、・・・、32列目の順で縦に連結し、1024次元のベクトルにしたものであり、X'(:)は推定された1995年における逆行列係数を同様にベクトルにしたものである。
【０１３７】
実験結果は以下の通りである。まず、得られた解のランクはr1=32, r2=32, r3=4であった。モード１とモード２に関しては正則化バランス定数γ1=0, γ2=0であるからこれは妥当な結果である。モード３に関してはフルランクであるものの、得られた特異値は図１０に示すように第１特異値に集中していることがわかる。従って低ランクの仮定はあてはまっていると考えられる。観測要素（８０％）のランダムな選択を１００回繰り返して得られた誤差は、未観測要素に関してerr1=0.11、補間に関する誤差err2=0.050であった。図１１に、真の1995年における逆行列係数と推定された逆行列係数を示す。ここでは縦軸にのみ、業種の名前を示しているが、横軸も縦軸と同じ業種に対応する。この図からも視覚的に推定結果が真の逆行列係数とよく一致していることがわかる。図１２、１３に、モード３に関する特異値分解（上記図１０）に付随する特異ベクトルを示す。上記図１２は第１成分を表わし、上記図１３は、第２成分を表わしている。また、図１２、１３の右側にはパターン基底Ｕ3のそれぞれ１列目（第１特異ベクトル）および２列目（第２特異ベクトル）に対応する特異値を乗じたものを示す。なお、最適化の結果得られたＵ3のそれぞれの列ベクトルは4つの要素を持ち、時間点1985、1990、2000、2005年に対応し、1995年のデータ（丸印で示す）は1990年および2000年のデータの内分から得られる予測値である。上記図１２、１３の左側には、Ｚ3の第１、第２特異値に対応する特異ベクトル、すなわち関係パターンＶ３の第１列、第２列をそれぞれ32×32行列に変形したものを示しており、図の右側に示すパターン基底が表現するダイナミクスに従って時間変化する逆行列係数の成分に対応する。上記図１２は産業界全体での連関が2000年まで年々減少し、2000年を境にして増加に転じたことを示している。上記図１３は農林水産業や鉄鋼などの伝統的な産業から他業種への影響が減少し、逆に通信・放送や対事業所サービスのような非伝統的な産業の他業種へ与える影響が強まっていることを表している。
【０１３８】
以上説明したように、本発明の実施の形態に係る時間関係データ解析装置によれば、時間変化するオブジェクト間の関係の強さを示す時間関係データを解析して、目的関数を最小化するように、時間関係データを表わす未知のテンソルＸを求めることによって、未観測のオブジェクト間の関係の強さを求めることができる。
【０１３９】
また、時間変化する関係データを少数のパターンに分解して、内在するパターンを発見することができる。パターンの数は補助変数更新部２４の特異値分解において、γｋ／η以上の値を持つ特異値とそれに対応する特異ベクトルのみ計算する処理により、予めパターン数を設定せずとも自動的に最適な値が求まる。
【０１４０】
また、解析する関係データの時間発展を複数の時間的要因に分解することが可能である。また、解析する関係データを複数のパターンの重みづけとして表現するため、観測に失敗した未知数や雑音の重畳した観測データを高精度に復元・予測できる。
【０１４１】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０１４２】
例えば、上述の時間関係データ解析装置は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。
【０１４３】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。
【符号の説明】
【０１４４】
１     入力部
２     解析部
３     記憶部
４     予測値計算部
５     出力部
２１   初期設定部
２２   変数計算部
２３   主変数更新部
２４   補助変数更新部
２５   ラグランジュ乗数更新部
２６   誤差計算部
２７   解析終了条件判定部
【図１】
0005684055.tif 000023
【図２】
0005684055.tif 000024
【図３Ａ】
0005684055.tif 000025
【図３Ｂ】
0005684055.tif 000026
【図３Ｃ】
0005684055.tif 000027
【図４】
0005684055.tif 000028
【図５】
0005684055.tif 000029
【図６】
0005684055.tif 000030
【図７】
0005684055.tif 000031
【図８】
0005684055.tif 000032
【図９】
0005684055.tif 000033
【図１０】
0005684055.tif 000034
【図１１】
0005684055.tif 000035
【図１２】
0005684055.tif 000036
【図１３】
0005684055.tif 000037