DEIM Forum 2023 2b-2-2

時系列テンソルデータのための将来予測
小幡 紘平†

川畑 光希†

松原 靖子†

櫻井 保志†† 大阪大学 産業科学研究所 産業科学 AI センター
E-mail: †obata88@sanken.osaka-u.ac.jp
あらまし 本論文では，大規模時系列テンソルデータのためのパターン検出手法である DDNF について述べる．DDNF
は，（sensor, location, time）の三つ組で構成される時系列テンソルデータに対し，sensor 間の直接的な依存関係を表す下層ネットワークと，location 間の直接的な依存関係を表す上層ネットワークから構成される階層ネットワークによってモデル化された重要なパターンを発見し，将来予測を実現する．具体的に，提案手法は，（a）時系列テンソルデータか
ら下層／上層ネットワークに基づいた解釈性の高いクラスタを発見する．
（b）計算量は入力データのサイズに対して線形である．そして，（c）ネットワークを用いた高精度な将来予測を実現する．人工データを用いたクラスタリング精度評価実験では最新の既存手法と比較して DDNF が大幅な精度向上を達成していることを明らかにした．また，実デー
タを用いた実験では提案手法が解釈性の高いクラスタを発見し，予測精度の向上を達成していることを確認した．
キーワード 時系列テンソル，特徴抽出，将来予測，グラフィカルラッソ
要であり，得られたクラスタの解釈性が高いことも望ましい．

1 まえがき車両走行センサ [1]，生体信号 [2], [3]，株価に代表される金融データ [4] など，様々なアプリケーションにおいて時系列デー
タが生成される．多くの場合，これらのデータは，タイムスタ
ンプが付与された複数のセンサ値が複数の地域（またはユーザー）から取得される．本論文では，二つの属性とタイムスタ
ンプから構成されるデータを時系列テンソルと呼ぶ．データ解析にとって，最も重要な課題の一つとして，将来予測に有用な
解釈可能なパターン（クラスタ）の発見がある．例えば，大気
成分濃度を正確に予測するためには，センサ間の依存関係と，センサの種類を考慮した地域間の依存関係によって特徴付けられた時間変化するパターンを発見することが重要である．本研
究では，時系列テンソルの依存関係を表現する要約情報を抽出することで，複数のセンサが互いにどのように相関するのか，
また，それぞれのセンサの振る舞いを考慮した地域間の相関関係を明らかにするパターンを発見する．さらに，パターンを用いた将来予測を達成する．
本論文で扱う問題は以下のとおりである．
問題：p 個のセンサ，d の地域，n 個のタイムスタンプから構成される時系列テンソル X ∈ Rp×d×n が与えられたとき，•

セグメントと分割点を発見する•

ネットワークに基づいたセグメントをグループ化し，類似

DTW [5] などの従来のクラスタリング手法は，値の距離由来の指標に依存するものが多く，クラスタの解釈が困難である．ク
ラスタリング手法の中には，変数間の関係性を考慮することで，解釈可能なクラスタを見つけることに重点を置いた手法もある [6], [7]．しかし，クラスタ数の自動決定，時系列テンソルの扱い，クラスタを利用した将来予測の達成には至っていない．
時系列テンソルは複雑なデータ構造を持つため，PARAFAC [8]は変数間の関係を犠牲にし，テンソルを低ランクのコアテン
ソルと行列の集合に分解する．近年のグラフニューラルネット
ワーク（GNN）[9] の進歩により，時系列データからより良い予測結果が得られるグラフ構造を構築することができるが，これらの多くはデータから一つの静的グラフを発見するため，クラスタを発見することはできない．
そこで本研究では，時系列テンソルから解釈性の高いクラス
タを発見し予測に用いる新しい手法として DDNF を提案する．
提案手法は時系列テンソルを二つの疎なネットワーク，下層／
上層ネットワーク，から構成される階層ネットワークによってモデル化する．ネットワークはグラフ構造をしており，変数がノードを，二つの変数間の直接的な依存関係がエッジを表す．
ここで，値に対して影響力の大きい次元を下位変数（センサ），
もう一方を上位変数（地域）と呼ぶ．下層ネットワークは下位変数間の関係を表し，上層ネットワークは下位変数の差を考慮した上位変数間の関係を表す．提案手法は時系列テンソルからパターン（クラスタ）を発見する•

時系列テンソルを予測する

ネットワークが異なるパターン／クラスタを発見し，それらのグラフを予測に利用する．
具体例. 図 1 は大気データと DDNF の出力結果例である．このデータは p = 6 個のセンサから取得された大気成分濃度が上記の問題を達成するには，隣接点が同じクラスタに属する傾向があるという，シーケンスの時間類似性を考慮する必要がある．また，データについて事前知識があることは稀である
ため，クラスタ数を指定せずにクラスタを発見することも必d = 12 の地域から取得された n = 8760 ポイントのデータである．図 1（a）は時系列テンソルのクラスタリング結果であり．同一のクラスタに含まれるセグメントは同一の色で表現されている．DDNF はテンソルを 4 つのセグメントと 3 つのクパターン発見. 時系列データの解析に関する研究はさまざまな分野で進められている [10] [11] [12]．中でも，時系列サブシー
ケンスのクラスタリングはデータを理解するために有用である．時系列データの教師なしクラスタリングの代表的な技術である，DTW（Dynamic Time Warping）[5] と K-menas は値の距
離由来の指標に基づいたクラスタリングを行い，データの構造(a) 3rd-order tensor

よりも実値を比べることに焦点を置いている．Li ら [13] が提案した，DynaMMo は線形動的システム（LDS: Linear Dynam-

ical System）に基づく手法で欠損を含む大規模時系列データ集
合から時系列のパターンを発見できる．Wang ら [14] によるpHMM（pattern-based hidden Markov model）は隠れマルコフモ
デル（HMM: Hidden Markov model）に基づく手法であり，時系列のセグメント化とクラスタリングのための動的モデルである．Mastubara ら [15] は多階層 HMM モデルを使ったパラメー
タフリーの手法として AutoPlait を提案している．これらの手法は，時系列の複雑な動的パターンを表現する能力はあるが，その一方で，ネットワーク構造を考慮していないため，クラス(b) Lower networks and upper networks

タの解釈には困難が伴う．

図 1 (a) 大気データにおける DDNF のクラスタリング結果．(b) 地図

時系列ネットワーク推定. 時系列情報を加味したネットワークにプロットされた上層ネットワークと下層ネットワーク．各測定

推定は経済データ，生体信号データの解析手法として研究され地点が上層ネットワークのノードを形成している．各測定地点で測定された物質が下層ネットワークのノードを形成している．変数間の関係値はエッジの太さと色で表されている．

ている [16]．グラフィカルラッソ [17] は静的なネットワーク推定手法であり，損失関数に ℓ1 正則化項を加味することで解釈が容易なスパースなネットワーク構造が推定できる [18]．Hallac

ラスタで表現している．DDNF は事前知識を必要とせず，大気ら [19] は文献 [17] に時系列情報を考慮したネットワーク推定

データを季節によって分類している．図 1（b）は各クラスタの手法である TVGL（Time Varying Graphical Lasso）を提案し，下層／上層ネットワークを示しており，ノードが変数，エッジ

Harutyunyan ら [20] は文献 [21] を改良した共分散行列推定手法の太さと色で変数間の偏相関係数を表している．提案手法が発として T-CorEx を提案した．Tomasi ら [22] は文献 [23] を時系見した 3 つクラスタはそれぞれ異なるネットワークを持つこと

列データに適応し，潜在状態を考慮した動的なネットワーク構がわかる．ノードが地図上にプロットされた上層ネットワーク

造を推定する手法である，LTGL（Latent variable Time-varying

では，近い場所同士のノードにしかエッジが存在しないことがGraphical Lasso）を提案した．これらの手法は，ネットワーク

わかる。7 章では提案手法がネットワークをテンソルの将来予構造の時系列変化をモデル化しており，前後のネットワーク構測に利用していることを示す．

造を比較することで変化点の検知は可能だが，クラスタリン

1. 1 本論文の貢献
DDNF は以下の特長を持つ．

グする能力はない．ネットワーク構造を基にしたクラスタリン
グ手法として Hallac ら [6] が提案した，TICC（Toeplitz Inverse

Covariance-based Clustering）と Tozzo ら [7] が提案した，TAGM(a) DDNF は類似したネットワークをもつ時系列パターン（ク（Time Adaptive Gaussian Model）がある．TICC はマルコフランラスタ）の個数と種類を把握し，データから解釈性の高い

ダムフィールド（MRF: Markov Random Field）とテプリッツ行

クラスタを発見する列を用い変数間に内在する関係を捉える手法であり，TAGM は

HMM と混合ガウスモデル（GMM: Gaussian Mixture Model）を(b) DDNF はクラスタリング精度と予測精度において，最新の既存手法と比較し性能向上を達成している(c) 提案したクラスタリングアルゴリズムは入力データの長さに対して線形である．つまり，長期間，高次元のデータに適応可能である融合した手法である．これらの手法は各サブシーケンスのネッ
トワーク構造に応じたクラスタを発見する．これにより，クラスタに解釈性を持たせ，従来のクラスタリング手法では発見できなかったパターンを発見することができる．しかしながら，これらの手法はクラスタ数を指定する必要があり，さらには時系列テンソルのクラスタリングに適していない．
時系列予測. 時系列データの将来予測は重要な課題であり，さ2 背景2. 1 関 連 研 究
関連研究は以下の 3 つに分類される．

まざまな研究がされている．近年，深層学習を用いた手法が非線形な特徴を捉える能力の高さから，よく研究されている．

LSTM [24] は時系列データから非線形性と時間依存性を捉える手法である．TCN [25] と WaveNet [26] は畳み込み層により，時系列データの時間的な相互作用を潜在状態にして捉える手法である．いくつかの深層学習由来の手法 [27] は変数間の依存性を捉えることができる．特に近年，GNN 由来の手法がグラフ畳み込み層により，大きな成功を収めている [28], [29]．しかしながら，これらの手法は変数間の関係を表すグラフを入力として必要とする．グラフを持たない時系列データに対応するため，予測に有用なグラフを生成する GNN が研究されている [9], [30], [31]．これらの手法はグラフがパターンに従って動的に変化するにも関わらず，一つの静的なグラフを使用する．

図2

DDNF は時系列テンソルから分割点とクラスタを発見する．各

結果として，動的なパターンと予測に有用なグラフを同時に発クラスタはグラフィカルラッソによって特徴づけられらスパース

見することは依然として困難な課題である．さらに，これまでな下層ネットワークと上層ネットワークから構成される．DDNFのところ，時系列テンソルに対するグラフに着目した研究はは下層ネットワークと上層ネットワークを階層ネットワークとして利用することにより，GNN を用いた時系列テンソルの予測をない．

可能にする．

3 問題定義され，それぞれが K 個の動的ネットワーク（クラスタ）のうち

ここでは本論文で取り組む問題について詳細に定義する．

一つに属するとする．cp を m 個のセグメントの開始点の集合

X ∈ Rp×d×n を p と d（センサと地域）と長さ n から構成とする，cp = {cp1 , cp2 , ..., cpm }．i 番目のセグメントの X はされる時系列テンソルとする．また，X のベクトル化として

Xcpi :cpi+1 （cpm+1 = n + 1）と表す．n ポイントを K に割り

X ∈ Rpd×n ，X = {X1 , X2 , . . . , Xn } を用いる．まず，X の静当てる，クラスタ割り当て集合を F = {f1 , f2 , . . . , fK } とする．

的なネットワーク（pd 変数間の関係）を推定する典型的な手法

fk ⊂ {1, 2, . . . , n} であり，クラスタ k に属する全てのデータはについて述べる，そして，X のより複雑なネットワークを表現

X[fk ] ⊂ X となる．また，各クラスタのパターンを表現するモデ

する．

ルパラメータ集合を Θ = {θ1 , θ2 , . . . , θK } とし，θk ∈ Rpd×pd はX[fk ] の変数間の関係を要約表現する疎なガウス逆共分散行列で

3. 1 グラフィカルラッソ
グラフィカルラッソは高い解釈性を持つグラフ表現を推定できる．多変量データ X が与えられたとき，グラフィカルラッ
ソは ℓ1 正則化項により，精度行列とも呼ばれるスパースな逆共分散行列 θ ∈ Rpd×pd , を推定する．グラフ表現から pd 変数

ある．よって，クラスタ情報集合は M = {M1 , M2 , . . . , MK }とする，Mk = {fk , θk }．最適な M は X における新たな特徴を表すため，M を明示的にテンソル予測に使用する．クラス
タ観点の回帰モデルを Ψ = {Ψ1 , Ψ2 , . . . , ΨK } とする．よって，
本研究の課題は以下のように定義される．

間の条件付き独立を理解できる．
（θi,j = 0 のとき，変数 i と j
は他の全ての変数値を与えられたとき，条件付き独立である．
）
具体的には，以下の式を最適化する：

minimizeθ∈S p λ||θ||od,1 − ll(X, θ),++ll(X, θ) =

n
∑
{

−

1
(Xi − µ)T θ(Xi − µ)
2

+

}
1
pd
log detθ −log(2π) ,
2
2

i=1

(1)

［問題 1］ 時系列テンソル X が与えられたとき，•クラスタ情報集合 M = {Mk }K
k=1•

回帰パタメータ集合 Ψ = {Ψk }K
k=1•

最適なセグメント数 m，クラスタ数 Kを発見する．

(2)p
ただし，θ は正定値対称行列（S++）である．|| · ||od,1 は対角成分を除いた ℓ1 ノルムである．ll(X, θ) は対数尤度関数であり µは X の平均である．正則化ハイパーパラメータ λ >= 0 により損失関数と ℓ1 正則化項のバランスを調整することで，スパース性を制御する．式 1 は凸最適化問題であり，本論文では交互

4 提案モデル
本章では解釈性の高いクラスタを発見し，予測に利用する

DDNF を提案する．図 2 は提案手法のモデル図である．まず，提案手法のモデルである階層ネットワーク θ について詳細に述べる．続いて，M を使用したテンソル予測について述べる．そして，クラスタ数を決定する基準を定義する．

方向乗数法（ADMM）[32] を用いることで高速に解く．ADMM
は凸最適化問題が大域解に収束することを保証する．

3. 2 DDNF 問題
実データの X はそれぞれが異なる関係／ネットワークを持
つ複数のパターンを含むため，一つの静的ネットワークでは表現しきれない．そこで本研究では，動的ネットワークに基づくテンソル予測問題に取り組む．X が m 個のセグメントに分割

4. 1 二階層グラフィカルラッソ
M を発見する問題は以下のように書ける．
arg min
Θ,FK∑

λ||θk ||od,1 − ll(X[fk ], θk ).

(3)

k=1

K と F が与えられたとき，θk をどのように定義すれば良いのだろうか．一般的なグラフィカルラッソで定義される θk は，ベクトル化したテンソルの全ての変数（pd）に対して関係を推定利用する，テンソル予測のフレームワークについて述べる．提する．しかしながら，この表現は二つの非時間次元の関係性を案手法は K 個の回帰モデル Ψ = {Ψk }K
k=1 を構築し，データそれぞれ表現するには無駄が多い．過剰に複雑なネットワークの性質に基づき最適なモデルを選択することで予測する．例えの形成を避けるため，本研究では，変数間の関係性を二つの側ば，データに季節性が存在する場合は，テストデータと同じ季面から捉えることで階層ネットワーク θk を表現する．
階層ネットワーク θ は，下位変数の依存関係を表す下層ネッ
トワーク A ∈ Rp×p （ai,j は変数 i と j の関係値を表す）と，上位変数の依存関係を表す上層ネットワーク B ∈ Rd×d から構成される．下位階層では p 個の変数（センサ）間に共通の関係性を捉えることを目的とするため，A はデータ全体を表現したものとなり共有される．一方，上位階層では p 個の変数の違いを加味した上で，d 個の変数（地域）間の関係性を捉える．そして，θ の表現を簡潔にするため，ある下位変数は下位変数と上位変数が共に異なる変数と関係しないという制約を加える．
よって，θ は pd × pd のテプリッツ行列となり，以下のように表される．

節の Ψk を使用する（後述の大気データ）．データの特徴が明


A


C
 2,1


C3,1

θ=
 ..
 .

 .
 .
 .


Cp,1

C1,2

C1,3

A

C2,3
..
.

C3,2
..
.

..

.



···
..
.
..

.

..

.

..

.

..

.

Cp−2,p−1..

.

Cp−1,p−2

A..

.

Cp,p−2Cp,p−1

C1,p

.. 
. 


.. 
. 

,

Cp−2,p 
Cp−1,p A

らかでない場合は，最新のクラスターの Ψk を使用する（株価データ）．Mk が与えられたとき，Ψk を構築する．グラフ構造 θk を利用するため，Ψk として GNN を使用する．θk は訓練とファインチューニングにおいて共通して入力する．訓練において X を入力データとして使用する．ファインチューニング
では X[fk ] を入力データとして使用することで，異なるネット
ワークを持つデータを取り除く．

4. 3 特徴抽出とデータ圧縮
ネットワークに基づいたクラスタリングがされた M を得る
ことは，精良な Ψ の構築につながる．最適な分割点とクラスタ
数 K の発見のために，最小記述長（MDL: minimum description

length）の概念を用いる．MDL は情報理論に基づくモデル選択基準の一つであり，直感的には，データをより圧縮できれば良いモデルとみなすことができる．本論文の目的を解決するために，新しい符号体形をグラフィカルラッソモデルに対して定義する．
ここでは，時系列テンソルを表現するための符号化スキーム

Ci,j = bi,j ∈ Rp×p は対角成分が bi,j ∈ B（i.e.,Ci,j = bi,j ·δi.j ）を導入する．簡潔に表すと，MDL を用いてデータを表現するの対角行列であり，δi.j はクロネッカーのデルタである．

ために必要な θ の最小数を求めることを目標とする．時系列テ

時系列テンソルから θ の要素である A と B を推定するために．グラフィカルラッソを拡張する．各クラスタにおいて上位

ンソル X が与えられたときのモデル M のよさは次の式で表現できる：

変数は同じ A を持つ．それゆえ，A は以下の式を最適化するこ⟨X ; M⟩ = α · ⟨M⟩ + ⟨X |M⟩,とで得られる：minimizeA∈S p λ||A||od,1 − llA (X , A),++llA (X , A) =

d ∑
n
∑
{

−

j=1 i=1(4)

ここで，⟨M⟩ は M を表現するためのコストを示し，⟨X |M⟩
は M が与えられたときの X の符号化のコストを示す．ハイ1(x:,j,i − µ)T A(x:,j,i − µ)2}1p
+ log detA − log(2π) ,22

パーパラメータ α > 0 はモデル表現コストと符号化コストのバ
ランスを調整し，モデルの複雑さを制御する．

4. 3. 1 モデル表現コスト(5)

モデル M の表現コストは以下の要素の総和から構成される．

x:j: は X の上位変数 j におけるデータである. µ は下位変数の平均である．B における制約のため, 対数尤度 Eq. (2) を拡張し，•クラスタの総数 K ：log∗ (K) 1•∑∗
各クラスタの観測値数： K
k=1 log (|fk |)•∑各クラスタの平均値 p × 1 と d × 1： K
k=1 ((p + d) × cF )各下位変数が同じ関係値を持つが，異なる平均を持つようにする．具体的には，以下の式を最適化する：

minimizeB∈S p λ||B||od,1 − llB (X , B),++

llB (X , B) =

p
n
∑
∑
{
j=1 i=1

−

(6)

•

1
(xj,:,i − µj )T B(xj,:,i − µj )2}d1
+ log detB − log(2π) ,
2
2

(8)

(7)

xj:: は X の下位変数 j におけるデータであり，µj ∈ R1 は x:j:

各クラスタの下層ネットワーク p × p と上層ネットワー∑∗
ク d × d： K
| (2 log(p) + cF ) + log (|Ak |=0| ) +
k=1 |Ak |=0∗|Bk |=0| (2 log(d) + cF ) + log (|Bk |=0| )ここで， | · |=0| は行列の非 0 要素の数を，cF は浮動小数点のコストを示す．2の平均である．Eq. (4)，Eq. (6) は凸最適化問題であり，ADMMによって解かれる．

4. 2 予測フレームワーク
GNN とファインチューニングによりクラスタ情報集合 M を

1：ここで， log∗ は整数のユニバーサル符号長を表す．
2：本論文では 4 × 8 ビットとする．

4. 3. 2 データ記述長
先述のとおり，本論文では θ を用いて X のパターンを表現するが，ここで重要なのは，推定したモデルが X を正しく表現しているかを判断する指標の導入である．ハフマン符号 [33] を用いた情報圧縮では，M が与えられた際の X の符号化コストを
負の対数尤度を用いて表現する．本論文では精度行列を用いる
代わりに分散共分散行列を用い，モデルの差異を正確に捉える．

⟨X |M⟩ =K∑−1llA (X[fk ], A−1k ) + llB (X[fk ], Bk ).

k=1

本論文の次の目標は上記のコスト関数 ⟨X ; M⟩ を最小化する
ようなクラスタ情報集合 M を発見することである．

(c) Right merge

図 3 CutPointDetector で比較される 3 つの分割点の候補の概略図．これらの分割点の候補の MDL コスト式（8）を比較する．

を増加させることで，最適なクラスタ数 K を発見する．MDL
コストには，F と Θ を推定する必要があり，片方の結果にはもう一方の結果が影響を与えるため，ClusterDetector では EM
アルゴリズムにより最適化を行う．E ステップでは，データ記
述長 ⟨X |M⟩ を最小化する F を計算する．具体的には，i 番目

cp

i+1
arg min ⟨X |θk , {j}j=cp
i

−1

⟩,

(9)

k∈{1,...,K}

前章では各クラスタを表現する階層ネットワークについて記述した．ここで重要なことは，どのようにして予測につながる
最適な分割点とクラスタ割り当てを発見するのかである．本章
では DDNF を実現するスケーラブルなアルゴリズムを提案する．アルゴリズム 1 に全体の処理の流れを示す．時系列テンソル X を与えられたとき，以下の二つの部分アルゴリズムにより式（8）を最小化する．

CutPointDetector ：セグメント数 m とその分割点 cp を発見する．

•(b) Left mergeのセグメントについて以下の式を解く．

5 提案アルゴリズム•(a) Original cp

ClusterDetector ：クラスタ数 K とクラスタ情報集合 M を
発見する．

cp

i+1
i 番目のセグメントのインデックス {j}j=cpi−1

を最適な割り当て先である fk ∈ F に挿入する．M ステップでは，各クラスタについて割り当てられたデータ X[fk ] のモデルパラメータ θ の要素である A と B を式（4）と式（6）により求める．ここで，

Θ の初期化は適当に行う．結果として，Ψ を推定するのに適切な M を得ることができる．

5. 3 理論的な分析［補助定理 1］ 提案手法の計算コストは O((p + d)n) である．
［証明 1］ DDNF の計算コストの大部分は CutPointDetector のイテレーション回数と，全ての A と B を推定する計算コス
トによる．最終的に分割点の数が 0 個になる場合を考える．

そして，DDNF は最終的なクラスタ情報集合 M を利用しテン

全ての A と B を推測するための合計計算コストは，それぞソル予測を行う（4. 2 章）.

れの計算コストの和であるため，A の場合について述べる．

5. 1 CutPointDetectorn ≫ p のとき，全ての A を推定する計算コストは O(pn)最初の目標はデータについての事前情報なしで X を m 個のセグメント（パターン）に分割することである．組合せ爆発を防ぐため，分割統治法を基にした CutPointDetector により最適な分割点を探索する．

CutPointDetector は小サイズのセグメント集合に分割した X
を，隣接セグメントが同じクラスタに属する傾向を利用し，隣接セグメントを MDL コストが減少する方向にマージする問である．CutPointDetector のイテレーション回数は，セグメント数が一つずつ減少する場合は |w| 回であるが，これは起こりそうにない．イテレーションごとにセグメント数が半減していく場合は log2 |w| 回である．n ≫ log2 |w|
であるため，B の推定を除いた計算コストは O(pn) とな
る．B の推定は A と同等のため，DDNF の計算コストはO((p + d)n) となる．

題を再帰的に解くことで分割点を発見する．図 3 に例を示す．

w = {wi }mi=1 はハイパーパラメータで，初期セグメントのサイ

6 評価実験

ズの集合（ひと月の日数ごと，小さい定数，など）とする．i 番
cpi+1
目のセグメントの MDL コストは ⟨X ; θi:i+1 , {j}j=cp
i

−1

⟩ とな

る．図 3（a）のように 3 つのセグメントが与えられたとき，中央のセグメントを片端のセグメントのどちらかとマージ（図 3（b）（c））するかを判断する．マージすることで MDL コストが3 つのセグメントの合計から減少する場合，不要な分割点を省
き，新しいセグメントに対して θ を推定する．この処理を全てのセグメントにおいて繰り返すことで，m は減少し収束する．

5. 2 ClusterDetector
DDNF はコスト関数 ⟨X ; M⟩ が減少する間，K = 1, 2, . . . , m

本章では，人工データに対する DDNF のクラスタリング精度と計算コストの検証を行う．クラスタリング精度比較に用いられる典型的な実データはネットワーク構造に基づいた正解ラベルが与えられていない．一方で，人工データでは明確なネットワーク構造のあるデータが生成可能で，ネットワーク構造に基づいたクラスタリング精度の比較が可能である．
人工データ生成. 人工データの生成，実験設計は文献 [6], [34] に従った．多変量正規分布 X ∼ N (0, θ −1 ) に従う時系列テンソル Rp×d×n をランダムグラフに基づき生成した．ネットワーク
構造に基づいたクラスタリング精度を評価するため，K 個の各表 1 4 つの異なるデータセットにおける DDNF と比較手法の macro-

Algorithm 1 DDNF(X , w)
1: Input: 3rd-order tensor X and initial segment sizes set w
2: Output: Cluster parameters Θ and cluster assignments F

F1 スコアによるクラスタリング精度（高いほど高精度）
Datasize

Pattern

DDNF

TICC TAGM AutoPlait

p=5, d=1

1,2,1
1,2,3,2,1
1,2,3,4,1,2,3,4
1,2,2,1,3,3,1

0.972
0.960
0.936
0.986

0.872
0.907
0.788
0.907

0.827
0.793
0.843
0.814

0.400
0.190
0.100
0.182

1,2,1
1,2,3,2,1
p=5, d=10 1,2,3,4,1,2,3,4
1,2,2,1,3,3,1

0.988
0.992
0.977
0.996

0.711
0.696
0.678
0.655

0.524
0.432
0.427
0.457

0.400
0.190
0.100
0.182

3: Initialize cp with w;
4: repeat5:

cp = C UT P OINT D ETECTOR(X , cp);

6: until cp is stable;
7: K = 1; Initialize Θ = {θ1 }; F = {{1, . . . , n}};8: Compute ⟨X ; Θ, F ⟩;9: repeat
10:

K = K + 1; Initialize Θ for K clusters;

パース性を制限するハイパーパラメータ λ = 0.1 とした．さ11:{Θ, F } = C LUSTER D ETECTOR(X , cp, K);
Compute ⟨X ; Θ, F ⟩;らに，これらの手法はクラスタ数の指定が必要であるため，正12:

しいクラスタ数を与え実験した．AutoPlait [15] は多階層 HMM

13: until ⟨X ; Θ, F ⟩ converges;
14: Estimate Ψk with X for each k = 1, . . . , K;
15: Fine-tune Ψk with X[fk ] for each k = 1, . . . , K;

ベースの自動クラスタリングアルゴリズムである．なお，これらの手法はテンソルを扱えないため，ベクトル化されたデータ

X が与えられた．DDNF には α と w のハイパーパラメータ16: return {M, Ψ}がある．実験では全ての初期セグメントサイズを一定 wi (s.t.,
クラスタの平均値は ⃗0 とした．以下の手順で，各クラスタの θを作成した [34]．
（ 1 ） 下層ネットワーク A ∈ R

i = 1, . . . , m) とした．訓練データを用いたパラメータチューニ
ングでは，wi = 4 とし，α を p と d に応じて変化させた．α が決定されると，w をコスト関数 Eq. (8) により決定できる．テ

p×p

を Erdős-Rényi モデルに従っ

ストデータを用いて wi = 4, 8, 16, 32 と変化させ，最小コストて作成する．全ノードペアについて，確率 20% でエッジ

Eq. (8) を示した結果を採用した．

を形成する．

多種類の人工データにおけるクラスタリング精度. クラスタリ（ 2 ） A の選ばれたエッジについて，aij ∼ U nif orm([−0.6, −0.3]∪[0.3, 0.6]) を設定する．また，A は対称行列 aij = aji とする．

ング精度を macro-F1 スコアで比較した結果を表 1 に示す．それぞれのセグメントについて p = 5, d = 1, 10, n = 100 とした（例 “1,2,1” では n = 300). 本手法が全てのデータセットにおいて最も高い平均精度を記録した．AutoPlait はネットワーク構（ 3 ） 上層ネットワーク B ∈ Rd×d を A と同じ手順で作成する．

造を考慮しないため，クラスタを発見できなかった．TICC とTAGM は正しいクラスタ数を与えられたにも関わらず，d = 1（ 4 ） テプリッツ行列 D ∈ Rpd×pd を A と B によって作成する．

の場合においても，平均精度について DDNF と比較して 9%（ 5 ） θ を正定値行列とするために，θ = D + (0.1 + |c|)I とす

以上低かった．これらの結果は提案手法が式（8）により，異なる．c = λmin (D) は D の最小固有値で，I は pd × pd の単位行列である．

るネットワークを持つセグメントを発見し，正しくクラスタに割り当てることができることを示している．
上位/下位変数を変化させたときのクラスタリング精度. 本手

評価指標. 次のような異なるセグメントの組み合わせの 4 つの法の符号化コスト関数はデータの次元数による影響を受ける．

データセットについて実験を行った [6] （“1,2,1”，“1,2,3,2,1”，“1,2,3,4,1,2,3,4” を例に取り p（p = 5 ∼ 50, d = 1, n = 800）と“1,2,3,4,1,2,3,4”，“1,2,2,1,3,3,3,1”）．それぞれのデータセットにd（p = 5, d = 5 ∼ 50, n = 800）を変化させ，精度への影響をつき 10 回実験を行い，macro-F1 スコアの平均と標準偏差を評価した. 図 4（a）（b）は次元数に対する macro-F1 スコアを記録した．macro-F1 スコアは，適合率（Precision）と再現率プロットした結果である. DDNF が全ての値において比較手法（Recall）の調和平均を各クラスタについて求め，平均したものを上回っていることがわかる．TICC と TAGM は p の増加によ

で，1 に近い値は高いクラスタリング精度を意味する. 3

り精度が低下したが，DDNF は α により高精度を維持してい比較手法. 与えられた人工データに対する提案手法のクラスタ

る．TICC と TAGM はテンソルを扱えないため，図 4（b）におリング精度を検証するために，最新の時系列クラスタリング手

ける精度が相対的に低い．

法と比較する．TICC [6]，および TAGM [7] はネットワーク構

サンプル数を変化させたときのクラスタリング精度. サンプ造に基づいたクラスタリングを行う手法である．TICC にはス

ル数の増加に対して精度を維持することは大規模時系列デー

パース性を制限するハイパーパラメータを λ = 0.1 とし，隣接タを扱うにあたり重要である．”1,2,3,4,1,2,3,4”を例に取り n

ポイントを同じクラスタに割り当てない際の罰則コストであ（p = 5, d = 1, n = 800 ∼ 80000）を変化させ，精度への影響るハイパーパラメータ β = 0, . . . 1000 と変化させ，訓練デーを評価した．図 4（c）は n に対する macro-F1 スコアをプロッ

タにおいて最良の結果を示した値を選択した．TAGM にはス

トした結果である．DDNF は精度において他手法を上回ってお3：macro-F1 =

1 ∑K
1
i
K
1/presicioni + 1/recalli

り，n が大きいほど高い精度を示している．これは wi が変化

回を訓練，4 周回を検証，8 周回をテストに使用した．
大気データ. 中国北京の d = 12 の地域から計測された p = 6 個の大気成分濃度センサについて 2013 年 3 月 1 日から 2016 年

2 月 29 日まで 1 時間間隔で取得された（n = 8760hours）大気(a) vs. lower mode
図4(b) vs. upper mode(c) vs. number of samples

データを使用する．それぞれ 1 年分のデータを使って，訓練，DDNF の（a）下位変数 p，（b）上位変数 d，（c）サンプル数 n に検証，テストを行う．

対するクラスタリング精度（macro-F1 スコア）.

評価指標. それぞれのデータセットにつき各手法で実験を 10 回行い，平均絶対誤差（MAE）と平均二乗誤差（MSE）の平均を記載した．
比較手法. 以下の比較手法を用いて実験を行った．

(a) vs. lower mode
図5(b) vs. upper mode(c) vs. number of samples

DDNF の（a）下位変数 p，（b）上位変数 d，（c）サンプル数 n に対する計算コスト.

•

Linear：全結合層．

•

LSTM [24]：エンコーダに LSTM，出力層に全結合層を使用．

•

seq2seq：LSTM をエンコーダとデコーダ共にに使用．

•

TCN [25]：Causal dilated convolution と残差ブロックによる手法．

•

WaveNet STN [35]：WaveNet にを時間空間正規化を追加した手法．

•

MTGNN [9]：多次元時系列データから予測に最適な隣接行列を生成する，GNN ベースの予測手法．

(a) DDNF result on map(b) 下層ネットワーク

各手法について “+C” と手法名の後に記載したものは DDNF に

図 6 車両走行センサデータコース A についてのクラスタリング結果：

よるクラスタリング結果 F を加えた手法である．提案手法は（a） クラスタ割り当てを地図上にプロットした結果．. 同一色がMTGNN を Ψ として使用し，Θ を事前グラフとして与えた．

同一クラスタを表す．
（i.e.，#紫 → 狭路， #青 → 高速度域，#赤→ 減速，#緑 → 中速度域）．
（b）各クラスタの下層ネットワーク．

コスト関数におけるハイパーパタメータ α はどれだけの数のパターンを発見するかを決定する．実データには正解ネット

させることで各セグメントのネットワークを正確に推定できるワークが存在しないため，検証データにおいて最良の結果を示からである．

した α を採用した．車両走行センサデータでは同じ座標にお提案手法の計算コスト. ここでは，提案手法の計算コストについけるデータは同じクラスタに属するとし，wi = 1（コース A：て検証する. 補助定理 1 では DDNF の計算コストが O((p+d)n)α = 0.7，コース B：α = 0.75，コース H：α = 1.5）とした．

であることを示した．図 5 に “1,2,3,4,1,2,3,4” の下位変数 p，上また，ドライバーは同時に走行しないため，上層ネットワーク

位変数 d，サンプル数 n を変化させたときの DDNF の計算コは使用しなかった．大気データでは，同じ日時におけるデータ

ストを示す. 補助定理 1 において示した通り，DDNF は提案ア

は同じクラスタに属するとし，wi = 24（1day）， α = 1.7 とルゴリズムにより，入力データサイズ p, d, n に対して線形で

した．

ある．

予測精度.DDNF による時系列テンソルの予測精度を検証した．
車両走行センサデータと大気データにおける結果は表 2 に示し

7 ケーススタディ
本章では，実データを対象とした実験により，DDNF の予測精度と DDNF が意味のあるネットワークを発見可能であることを示す．
車両走行センサデータ.3 つの異なる状況のコースから取得された車両走行センサデータを用いる．コース A は市街地を走行したデータで（n = 3241m），コース B は高架下の車通りコー
スのデータで（n = 2593m），コース H は高速道路を走行した
データ（n = 4000m）である．7 つのセンサーの値が 1m ごとに計測された：（ブレーキ，速度，前後加速度，左右加速度，ア
クセル開度，ハンドル角，燃費）それぞれ 40 周回あり，28 周た．最良の比較手法に対する DDNF の改善率が Improvementsに記されている．DDNF が全てにおいて最も良い結果を示した
ことがわかる．これは DDNF がクラスタ割り当てと各クラス
タの変数間の依存関係を示すグラフが与えられたからである．
次に，データを分割して学習する効果を “+C” との比較で検証した．結果から，変数間の依存関係が異なるデータセットに分割して学習する有効性がわかる．
解釈性.DDNF におけるクラスタリング結果が理にかなっていることを示す．大気データについての結果は 1 章（図 1）で述
べた．図 6 は車両走行センサデータコース A についてのクラスタリング結果である．図 6（a）は地図上にプロットされたク
ラスタリング結果であり，図 6（b）は各クラスタの下層ネット

表 2 車両走行センサデータと大気データにおける精度比較．“+C” は提案手法によるクラスタリ
ング結果 F が使用された手法を表す．
Course A
Models5m
MAE MSE

10 m
MAE MSE

Course B
5m
10 m
MAE MSE MAE MSE

Course H
5m
10 m
MAE MSE MAE MSE

1 hour
MAE MSEAir quality
6 hour
MAE MSE12 hour
MAE MSE

LinearLinear+C

0.339
0.324

0.419
0.395

0.472
0.457

0.649
0.617

0.344
0.337

0.384
0.377

0.479
0.471

0.629
0.616

0.147
0.146

0.078
0.077

0.242
0.240

0.182
0.179

0.194
0.177

0.097
0.087

0.448
0.436

0.445
0.427

0.536
0.529

0.617
0.607

LSTM
LSTM+C

0.274
0.269

0.299
0.288

0.405
0.397

0.517
0.494

0.270
0.268

0.272
0.271

0.388
0.386

0.464
0.461

0.135
0.133

0.067
0.066

0.224
0.221

0.164
0.162

0.201
0.200

0.122
0.122

0.444
0.445

0.488
0.491

0.547
0.541

0.683
0.661

seq2seqseq2seq+C

0.261
0.256

0.280
0.268

0.394
0.382

0.504
0.482

0.260
0.255

0.261
0.258

0.373
0.367

0.460
0.453

0.133
0.130

0.066
0.065

0.219
0.215

0.163
0.161

0.152
0.151

0.080
0.079

0.447
0.441

0.492
0.483

0.560
0.554

0.714
0.699

TCN
TCN+C

0.287
0.279

0.334
0.313

0.418
0.407

0.554
0.521

0.284
0.282

0.297
0.293

0.408
0.403

0.496
0.488

0.137
0.136

0.069
0.068

0.229
0.225

0.165
0.162

0.240
0.241

0.153
0.153

0.414
0.411

0.415
0.411

0.514
0.511

0.594
0.596

WaveNet STN
WaveNet STN+C

0.263
0.257

0.290
0.273

0.396
0.385

0.506
0.478

0.253
0.248

0.257
0.252

0.374
0.369

0.457
0.450

0.127
0.125

0.063
0.062

0.216
0.209

0.155
0.152

0.141
0.141

0.071
0.071

0.386
0.386

0.377
0.376

0.490
0.485

0.539
0.536

MTGNN
MTGNN+C

0.250
0.246

0.274
0.268

0.381
0.373

0.487
0.465

0.246
0.239

0.252
0.244

0.361
0.356

0.440
0.431

0.124
0.121

0.063
0.062

0.206
0.203

0.154
0.151

0.143
0.141

0.074
0.073

0.382
0.382

0.379
0.373

0.482
0.480

0.531
0.534

DDNF

0.243

0.254

0.369

0.455

0.237

0.240

0.353

0.425

0.120

0.061

0.202

0.149

0.136

0.069

0.379

0.369

0.475

0.522

Improvements+3.0% +7.4% +3.2% +6.7% +3.4% +4.8% +2.2% +3.4% +3.2% +3.7% +1.7% +3.0% +3.2% +3.2% +1.0% +2.1% +1.5% +1.8%

ワークである．同一色は同一クラスタを表している．図 6（a）のクラスタ 1 はアクセル，ブレーキシステムに関するエッジが他のクラスタよりも多い．これはこのパターンが繊細な速度調[5]整が必要な狭路で観測されることを示し，地図からもそれがわかる．一方で，図 6（a）のクラスタ 3 はブレーキと前後加速度に太いエッジが形成されている．地図を確認することで，この[6]パターンが信号の前に多く観測されていることがわかる．
[7]8 む す び
本論文では，下層ネットワークと上層ネットワークにより特[8]徴付けられた解釈性の高いクラスタを検出し予測に利用する手法として，DDNF を提案した．DDNF は時系列テンソルを二[9]

つのネットワークによる階層構造に基いてモデル化することで解釈性の高いクラスタを発見できる．人工データと実データを[10]用いた実験により，DDNF は最新の既存手法と比べてより高いクラスタリング精度と予測精度を持ち，解釈性の高いネット
ワーク構造を持つクラスタを発見することを示した．さらに，[11]

提案手法のアルゴリムは入力データに対して線形であり，長期間，高次元のデータに適応可能であることを示した．
謝辞 本研究の一部は JSPS 科研費，JP20H00585，JP21H03446，[12]

JP22K17896, 国立 研 究 開発法 人 情報通 信 研 究 機 構 委 託 研 究
NICT 03501, 総務省 SCOPE JP192107004，JSTAIP 加速課題 JPMJCR21U4，ERCA 環境研究総合推進費 JPMEERF20201R02,の助成を受けたものです．

文

[13]

献

[1] Miyajima, C., Nishiwaki, Y., Ozawa, K., Wakita, T., Itou, K., Takeda,
K. and Itakura, F.: Driver modeling based on driving behavior and its
evaluation in driver identiﬁcation, IEEE, Vol. 95, No. 2, pp. 427–437(2007).
[2] Hirano, S. and Tsumoto, S.: Cluster analysis of time-series medical
data based on the trajectory representation and multiscale comparison techniques, ICDM, IEEE, pp. 896–901 (2006).
[3] Monti, R. P., Hellyer, P., Sharp, D., Leech, R., Anagnostopoulos, C.
and Montana, G.: Estimating time-varying brain connectivity networks from functional MRI time series, NeuroImage, Vol. 103, pp.
427–443 (2014).
[4] Chiang, T. C., Jeon, B. N. and Li, H.: Dynamic correlation anal-

[14]
[15]
[16]

[17]

[18]

ysis of ﬁnancial contagion: Evidence from Asian markets, Journalof International Money and Finance, Vol. 26, No. 7, pp. 1206–1228(2007).
Berndt, D. J. and Clifford, J.: Using Dynamic Time Warping to Find
Patterns in Time Series, Knowledge Discovery in Databases: Papers
from the 1994 AAAI Workshop, Seattle, Washington, USA, July 1994.
Technical Report WS-94-03, pp. 359–370 (1994).
Hallac, D., Vare, S., Boyd, S. P. and Leskovec, J.: Toeplitz Inverse Covariance-Based Clustering of Multivariate Time Series Data,
KDD, pp. 215–223 (2017).
Tozzo, V., Ciech, F., Garbarino, D. and Verri, A.: Statistical Models
Coupling Allows for Complex Local Multivariate Time Series Analysis, KDD, pp. 1593–1603 (2021).
Kolda, T. G. and Bader, B. W.: Tensor decompositions and applications, SIAM review, Vol. 51, No. 3, pp. 455–500 (2009).
Wu, Z., Pan, S., Long, G., Jiang, J., Chang, X. and Zhang, C.: Connecting the dots: Multivariate time series forecasting with graph neural networks, KDD, pp. 753–763 (2020).
Kumar, S., Zhang, X. and Leskovec, J.: Predicting Dynamic Embedding Trajectory in Temporal Interaction Networks, Proceedingsof the 25th ACM SIGKDD international conference on Knowledge
discovery and data mining, ACM (2019).
Wen, Q., Gao, J., Song, X., Sun, L., Xu, H. and Zhu, S.: RobustSTL:
A robust seasonal-trend decomposition algorithm for long time series, Proceedings of the AAAI Conference on Artiﬁcial Intelligence,Vol. 33, pp. 5409–5416 (2019).
Hallac, D., Bhooshan, S., Chen, M. H., Abida, K., Sosic, R. andLeskovec, J.: Drive2Vec: Multiscale State-Space Embedding of Vehicular Sensor Data, 21st International Conference on Intelligent
Transportation Systems, ITSC 2018, Maui, HI, USA, November 47, 2018 (Zhang, W., Bayen, A. M., Medina, J. J. S. and Barth, M. J.,eds.), IEEE, pp. 3233–3238 (2018).
Li, L., McCann, J., Pollard, N. S. and Faloutsos, C.: DynaMMo:mining and summarization of coevolving sequences with missing
values, KDD, pp. 507–516 (2009).
Wang, P., Wang, H. and Wang, W.: Finding semantics in time series,
SIGMOD, pp. 385–396 (2011).
Matsubara, Y., Sakurai, Y. and Faloutsos, C.: AutoPlait: Automatic
Mining of Co-evolving Time Sequences, SIGMOD (2014).
Mohan, K., Chung, M., Han, S., Witten, D., Lee, S.-I. and Fazel,
M.: Structured learning of Gaussian graphical models, Advances inneural information processing systems, Vol. 25 (2012).
Friedman, J., Hastie, T. and Tibshirani, R.: Sparse inverse covarianceestimation with the graphical lasso, Biostatistics, Vol. 9, No. 3, pp.
432–441 (2008).
Tomasi, F., Tozzo, V., Verri, A. and Salzo, S.: Forward-Backward
Splitting for Time-Varying Graphical Models, Proceedings of theNinth International Conference on Probabilistic Graphical Models

[19]
[20]

[21]

[22]
[23]

[24]
[25]

[26]

[27]

[28]

[29]

[30]

[31]

[32]

[33]

[34]

[35]

(Kratochvı́l, V. and Studený, M., eds.), Proceedings of Machine
Learning Research, Vol. 72, PMLR, pp. 475–486 (2018).
Hallac, D., Park, Y., Boyd, S. P. and Leskovec, J.: Network Inference
via the Time-Varying Graphical Lasso, KDD, pp. 205–213 (2017).
Harutyunyan, H., Moyer, D., Khachatrian, H., Steeg, G. V. and
Galstyan, A.: Efﬁcient Covariance Estimation from Temporal Data,arXiv preprint arXiv:1905.13276 (2019).
Steeg, G. V., Harutyunyan, H., Moyer, D. and Galstyan, A.: Fast
Structure Learning with Modular Regularization, Curran Associates
Inc. (2019).
Tomasi, F., Tozzo, V., Salzo, S. and Verri, A.: Latent Variable Timevarying Network Inference, KDD, pp. 2338–2346 (2018).
Chandrasekaran, V., Parrilo, P. A. and Willsky, A. S.: Latent variable graphical model selection via convex optimization, 2010 48th
Annual Allerton Conference on Communication, Control, and Computing (Allerton), pp. 1610–1613 (2010).
Hochreiter, S. and Schmidhuber, J.: Long short-term memory, Neural computation, Vol. 9, No. 8, pp. 1735–1780 (1997).
Bai, S., Kolter, J. Z. and Koltun, V.: An Empirical Evaluation ofGeneric Convolutional and Recurrent Networks for Sequence Modeling, arXiv:1803.01271 (2018).
Oord, A. v. d., Dieleman, S., Zen, H., Simonyan, K., Vinyals,
O., Graves, A., Kalchbrenner, N., Senior, A. and Kavukcuoglu,
K.: Wavenet: A generative model for raw audio, arXiv preprint
arXiv:1609.03499 (2016).
Shih, S.-Y., Sun, F.-K. and Lee, H.-y.: Temporal pattern attention for
multivariate time series forecasting, Machine Learning, pp. 1421–
1441 (2019).
Yu, B., Yin, H. and Zhu, Z.: Spatio-temporal graph convolutional
networks: A deep learning framework for trafﬁc forecasting, arXiv
preprint arXiv:1709.04875 (2017).
Zheng, C., Fan, X., Wang, C. and Qi, J.: Gman: A graph multiattention network for trafﬁc prediction, Proceedings of the AAAI conference on artiﬁcial intelligence, pp. 1234–1241 (2020).
Wu, Z., Pan, S., Long, G., Jiang, J. and Zhang, C.: Graph Wavenet
for Deep Spatial-Temporal Graph Modeling, IJCAI, p. 1907–1913(2019).
Shang, C., Chen, J. and Bi, J.: Discrete graph structure learning forforecasting multiple time series, arXiv preprint arXiv:2101.06861(2021).
Boyd, S. P., Parikh, N., Chu, E., Peleato, B. and Eckstein, J.: Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers, Found. Trends Mach. Learn., Vol. 3,
No. 1, pp. 1–122 (2011).
Böhm, C., Faloutsos, C., Pan, J.-Y. and Plant, C.: Ric: Parameter-freenoise-robust clustering, ACM Transactions on Knowledge Discovery
from Data (TKDD), Vol. 1, No. 3, pp. 10–es (2007).
Mohan, K., London, P., Fazel, M., Witten, D. and Lee, S.-I.: NodeBased Learning of Multiple Gaussian Graphical Models, J. Mach.
Learn. Res., Vol. 15, No. 1, p. 445–488 (2014).
Deng, J., Chen, X., Jiang, R., Song, X. and Tsang, I. W.: ST-Norm:
Spatial and Temporal Normalization for Multi-variate Time Series
Forecasting, KDD, pp. 269–278 (2021).