長いセンサー波形を眺めていて、「この区間と似た形が、過去のどこかにもあったはずだ」と思ったことはありませんか。あるいは逆に、「この一瞬だけ、どこにも似ていない奇妙な形が出ている」と気づきたいことはないでしょうか。心電図で同じ拍動パターンを探す、製造ラインのセンサ波形から繰り返し現れる動作を見つける、いつもと違う異常な部分列を拾う——これらはすべて、「時系列の部分列(サブシーケンス)どうしを総当たりで比べる」という同じ計算に行き着きます。
素朴にやると、長さ $n$ の系列から取れる窓は約 $n$ 個、その全ペアを比べると $O(n^2)$ 通り。各比較も窓長 $m$ に比例するので $O(n^2 m)$ になり、少し長い系列でもすぐ手に負えなくなります。Matrix Profile(Yeh et al., ICDM 2016)は、この「部分列の全ペア比較」を、ある一本のベクトルに圧縮して持つことで、モチーフ発見・ディスコード(異常)発見・類似検索を一気に・厳密に・高速にこなす仕組みです。
Matrix Profile が使える場面は驚くほど広く、たとえば次のようなものです。
- 異常検知：多変量センサー時系列から「どこにも似ていない最も孤立した部分列(ディスコード)」を異常として拾う。教師ラベルなしで動く。
- モチーフ発見：繰り返し現れる典型パターン(モチーフ)を見つけ、心拍・歩行・機械の動作サイクルなどの「定型動作」を自動抽出する。
本記事では、Matrix Profile の中心にある z正規化ユークリッド距離、それを高速に計算する MASS(FFT を使って距離プロファイルを $O(n\log n)$ にする)、全体を組み立てる STAMP/STOMP/SCRIMP++、そしてモチーフとディスコードの抽出までを、数式の導出と Python のスクラッチ実装で順に組み上げます。すべてのコードは実際に動かし、本文に書いた検出位置や距離の値と一致することを確認しています。
本記事の内容
- 部分列・距離プロファイル・Matrix Profile という3段重ねの考え方
- z正規化ユークリッド距離の定義と、Pearson相関との関係
- MASS：スライディング内積を FFT で計算して $O(n\log n)$ にする導出
- STAMP/STOMP/SCRIMP++ の違いと、モチーフ・ディスコード抽出のスクラッチ実装
前提・関連記事
Matrix Profileの全体像：3段重ねで理解する
いきなり定義に入る前に、Matrix Profile が何を計算しているのか、絵で押さえておきましょう。鍵になる用語は3つだけです。
- 部分列(サブシーケンス)：長さ $n$ の時系列 $T$ から、長さ $m$ の連続した切り出し $T_{i,m} = (t_i, t_{i+1}, \dots, t_{i+m-1})$ を取り出したもの。窓を1ステップずつスライドさせると $n-m+1$ 個の部分列が得られます。
- 距離プロファイル：ある1つの部分列(クエリ窓)を、系列の全位置に重ねていったときの距離の列。長さ $n-m+1$ のベクトルで、谷になっている位置が「クエリに似ている場所」です。
- Matrix Profile：各部分列について「自分以外で最も似た部分列までの距離」だけを集めたベクトル。これも長さ $n-m+1$ です。
イメージとしては、まず全ペアの距離を並べた巨大な「距離行列(distance matrix)」を考えます。その各行が「ある窓から見た距離プロファイル」、各行の最小値を縦に集めたものが Matrix Profile です。つまり Matrix Profile は、$n\times n$ の行列を見ずに、その「列ごとの最小値プロファイル」だけを賢く計算する、という発想です。
上の図では、規則的に繰り返す背景の波の中に、同じ形をした「二つこぶ」の部分列が2か所(クエリ窓と最も似た窓)埋め込まれています。クエリ窓から見ると、最も似ているのはもう一方のこぶ窓であり、そこまでの距離が Matrix Profile のその位置の値になります。中央には背景にもこぶにも似ない高周波のバーストがありますが、これは後でディスコード(異常)として浮かび上がります。
Matrix Profile の威力は、このベクトルさえ手に入れば、谷を見ればモチーフ、山を見ればディスコードと、ひと目で構造が読めることです。1本のベクトルが、繰り返しパターンの場所も、異常の場所も、似た区間の対応関係(Profile Index)も同時に教えてくれる——この「情報密度の高さ」が、Matrix Profile が時系列解析の万能ツールと呼ばれる理由です。
なお、Matrix Profile が解いているのは正確には自己類似結合(self similarity join)という問題です。ひとつの系列の中で、すべての部分列ペアの類似度を求め、各部分列に最近傍を割り当てる。これは2つの異なる系列の間でも定義でき(AB結合)、「系列Aの各窓に最も似た系列Bの窓」を探す結合類似検索になります。製造ラインで「正常な基準波形(B)に対して、いま流れている波形(A)のどこが似ていないか」を測る、といった使い方がまさにこれです。本記事では理解しやすい自己結合を中心に進めます。では、その値を決める「距離」をきちんと定義するところから始めましょう。
z正規化ユークリッド距離：形だけを比べる
部分列どうしの「似ている／似ていない」を測るには距離が要ります。素朴には、2つの長さ $m$ の窓 $\bm{a}, \bm{b}$ のユークリッド距離
$$ d(\bm{a},\bm{b}) = \sqrt{\sum_{k=1}^{m}(a_k-b_k)^2} $$
を使えばよさそうです。ところがこれには大きな落とし穴があります。同じ形なのに、上下にずれていたり、振幅が違っていたりするだけで距離が大きくなるのです。たとえば心電図の同じ拍動でも、ベースラインが少し上がっているだけで「別物」と判定されてしまいます。私たちが本当に比べたいのは、高さやスケールではなく「形」です。
そこで Matrix Profile では、各窓を比較の前に z正規化(標準化) します。窓の平均 $\mu$ を引いて標準偏差 $\sigma$ で割り、平均0・分散1にそろえてから距離を測ります。
$$ \hat{a}_k = \frac{a_k – \mu_a}{\sigma_a}, \qquad \hat{b}_k = \frac{b_k – \mu_b}{\sigma_b} $$
このz正規化した窓どうしのユークリッド距離を、z正規化ユークリッド距離と呼びます。
$$ \begin{equation} d_{\text{z}}(\bm{a},\bm{b}) = \sqrt{\sum_{k=1}^{m}\left(\hat{a}_k – \hat{b}_k\right)^2} \end{equation} $$
平均を引くことで上下のオフセットが消え、標準偏差で割ることで振幅のスケールが消えます。残るのは純粋な「形」だけ。これが Matrix Profile の比較がオフセット・スケール不変になる理由です。
実際に確かめてみましょう。同じ形の波 A と、それを2.2倍して+3だけ持ち上げた波 B を比べます。
import numpy as np
def znorm(s):
s = np.asarray(s, float)
mu, sd = s.mean(), s.std()
return (s - mu) / sd if sd > 1e-9 else s - mu
def zdist(a, b): # z正規化ユークリッド距離
return np.sqrt(np.sum((znorm(a) - znorm(b)) ** 2))
m = 80
base = np.linspace(0, 1, m)
shape = np.sin(2*np.pi*2*base) + 0.5*base
a = shape
b = 2.2*shape + 3.0 # スケール2.2倍・オフセット+3
print("生の距離 :", round(np.sqrt(((a-b)**2).sum()), 3))
print("z正規化距離 :", round(zdist(a, b), 6))
このコードを実行すると、生のユークリッド距離は約30.4と大きいのに対し、z正規化距離は 0.000 になります。形がまったく同じなら、オフセットやスケールがどれだけ違っても距離はゼロ——これがz正規化の効果です。
左の生のままでは、形が同じなのに高さと振幅が違うだけで2つの波は大きく離れています。右のz正規化後は、AとBがぴたりと重なり、距離がほぼゼロになります。Matrix Profile が「背景の高さが違っても同じ動作パターンを同一視できる」のは、この正規化のおかげです。
z正規化距離は、実は Pearson相関係数と表裏一体です。次に、この関係を導出しておきましょう。これが後の高速化(MASS)の鍵になります。
z正規化距離とPearson相関の関係
z正規化距離をそのまま計算するには、毎回 $\hat{a}_k – \hat{b}_k$ を作って二乗和を取る必要があります。これを展開すると、もっと計算に都合のよい形になります。ゴールは、z正規化距離を「窓どうしの内積(相関)」だけで書くことです。
まず二乗を展開します。
$$ d_{\text{z}}^2 = \sum_{k=1}^{m}\left(\hat{a}_k – \hat{b}_k\right)^2 = \sum_{k=1}^{m}\hat{a}_k^2 – 2\sum_{k=1}^{m}\hat{a}_k\hat{b}_k + \sum_{k=1}^{m}\hat{b}_k^2 $$
ここで $\hat{a}_k, \hat{b}_k$ は z正規化されているので、定義から分散が1、すなわち $\sum_k \hat{a}_k^2 = m$ と $\sum_k \hat{b}_k^2 = m$ が成り立ちます(平均0・分散1なので二乗和は $m$)。これを代入すると、第1項と第3項が両方 $m$ になります。
$$ d_{\text{z}}^2 = m – 2\sum_{k=1}^{m}\hat{a}_k\hat{b}_k + m = 2m – 2\sum_{k=1}^{m}\hat{a}_k\hat{b}_k $$
残った $\sum_k \hat{a}_k\hat{b}_k$ は、z正規化されたベクトルどうしの内積です。これはまさに、もとの窓 $\bm{a}, \bm{b}$ の Pearson相関係数 $\rho_{ab}$ に $m$ を掛けたものです。
$$ \sum_{k=1}^{m}\hat{a}_k\hat{b}_k = \sum_{k=1}^{m}\frac{(a_k-\mu_a)(b_k-\mu_b)}{\sigma_a\sigma_b} = m\,\rho_{ab} $$
これを代入すると、z正規化距離は相関係数だけできれいに書けます。
$$ \begin{equation} d_{\text{z}}(\bm{a},\bm{b}) = \sqrt{2m\left(1 – \rho_{ab}\right)} \end{equation} $$
この式の意味は明快です。相関 $\rho_{ab}=1$(完全に同じ形)なら距離は0、$\rho_{ab}=0$(無相関)なら $\sqrt{2m}$、$\rho_{ab}=-1$(上下反転)なら最大の $2\sqrt{m}$ になります。距離を測る問題が、相関(=内積)を測る問題に化けた——これが本質的な進歩です。なぜなら内積、特に「ある窓を全位置にスライドさせながら取る内積」は、畳み込み(相互相関)そのものであり、FFTで一気に計算できるからです。
次の節では、この事実を使って距離プロファイルを $O(n\log n)$ で求める MASS を組み立てます。
距離プロファイルとMASS：FFTで一気に距離を出す
1つのクエリ窓 $Q$(長さ $m$)を決めて、それを系列 $T$(長さ $n$)の全位置にスライドさせ、各位置との距離を並べたものが距離プロファイルでした。これをまともに計算すると、$n-m+1$ か所それぞれで長さ $m$ の距離を測るので $O(nm)$ かかります。
ここで先ほどの式(2)を思い出します。距離は相関、つまり「z正規化したクエリと、各位置の窓の内積」で決まりました。各位置 $j$ での(生の)スライディング内積を
$$ \text{QT}[j] = \sum_{k=0}^{m-1} \hat{Q}_k \, t_{j+k} $$
と書くと(ここで $\hat{Q}$ はz正規化したクエリ)、これは $\hat{Q}$ を反転させたものと $T$ の畳み込みにほかなりません。畳み込みは FFT を使えば $O(n\log n)$ で計算できます。さらに、各位置の窓の平均 $\mu_j$ と標準偏差 $\sigma_j$ は、累積和(prefix sum)を一度作っておけば $O(1)$ で取り出せます。
$$ \mu_j = \frac{1}{m}\sum_{k=0}^{m-1} t_{j+k}, \qquad \sigma_j = \sqrt{\frac{1}{m}\sum_{k=0}^{m-1} t_{j+k}^2 – \mu_j^2} $$
z正規化クエリ $\hat{Q}$(平均0・分散1)を使うと、位置 $j$ での相関は $\rho_j = \dfrac{\text{QT}[j]}{m\,\sigma_j}$ となり、これを式(2)に入れると距離プロファイルが一発で出ます。
$$ \begin{equation} D[j] = \sqrt{2m\left(1 – \frac{\text{QT}[j]}{m\,\sigma_j}\right)} \end{equation} $$
この手順全体が MASS(Mueen’s Algorithm for Similarity Search) です。支配項は FFT の $O(n\log n)$ なので、ナイーブな $O(nm)$ より圧倒的に速くなります。
距離プロファイルがどんな形をしているか、先にイメージをつかんでおきましょう。クエリ窓を1か所のモチーフに固定し、それを系列全体にスライドさせて距離を測ると、次のようなプロファイルが得られます。
上が時系列、下が距離プロファイルです。クエリ窓(モチーフ1)を全位置に重ねると、同じ形のモチーフ2の位置で距離が深い谷になるのが見えます。クエリ自身の位置(と隣接窓)は自明に距離ゼロですが、それを除けば、最も深い谷が「最も似ている相手」を指します。距離プロファイルは「1つの窓から見た世界地図」であり、その谷の場所が類似検索の答えになります。では、これを高速に計算する MASS を実装しましょう。
import numpy as np
def mass(Q, T):
"""クエリQ(長さm)とT(長さn)の全位置に対するz正規化距離プロファイル。O(n log n)。"""
n, m = len(T), len(Q)
Qz = znorm(Q) # クエリをz正規化(平均0・分散1)
# 各窓の平均・標準偏差を累積和でO(n)に
cs = np.cumsum(np.insert(T, 0, 0.0))
cs2 = np.cumsum(np.insert(T**2, 0, 0.0))
s = cs[m:] - cs[:n-m+1]
s2 = cs2[m:] - cs2[:n-m+1]
mu = s / m
sig = np.sqrt(np.maximum(s2/m - mu**2, 1e-12))
# スライディング内積 QT[j]=Σ Qz[k]·T[j+k] をFFTで一括計算
Qr = np.concatenate([Qz[::-1], np.zeros(n - m)]) # 反転+零詰め=畳み込みでスライド内積に
QT = np.fft.irfft(np.fft.rfft(T) * np.fft.rfft(Qr, n), n)[m-1:n]
d2 = 2*m*(1 - QT/(m*sig)) # 式(3)
return np.sqrt(np.maximum(d2, 0))
このMASSが、素朴な総当たり(各位置で zdist を計算)と本当に同じ値を返すかを確認します。
def gen_demo(): # 反復背景 + 同形モチーフ2回 + 異常バースト
rng = np.random.default_rng(7)
n, m = 900, 80
t = np.arange(n)
x = np.sin(2*np.pi*t/60) + 0.07*rng.standard_normal(n)
base = np.linspace(0, 1, m)
motif = (np.sin(2*np.pi*2*base)
+ 1.1*np.exp(-((base-0.30)**2)/0.01)
+ 1.1*np.exp(-((base-0.70)**2)/0.01))
common = 0.04*rng.standard_normal(m) # 2コピーで共有する微小ノイズ
for c0 in [150, 560]:
x[c0:c0+m] = motif + common + 0.02*rng.standard_normal(m)
x[380:380+m] = 3.0*np.sin(2*np.pi*9*np.linspace(0,1,m))*np.hanning(m) # 異常
return x.astype(float), m
T, m = gen_demo()
qi = 150
Q = T[qi:qi+m]
dp_mass = mass(Q, T)
dp_bf = np.array([zdist(Q, T[j:j+m]) for j in range(len(T)-m+1)])
print("MASS vs ブルートフォース 最大誤差 :", float(np.max(np.abs(dp_mass - dp_bf))))
実行すると、最大誤差は $10^{-6}$ 未満(浮動小数点の丸め程度)で、MASSは厳密に同じ距離プロファイルを返していることが確認できます。違いは速度だけです。
左の図で、灰色の太線(ブルートフォース $O(nm)$)と赤の破線(MASS $O(n\log n)$)は完全に重なっています。距離プロファイルには2つの深い谷があり、それが埋め込んだ2か所のモチーフ位置に対応します。右の図は計算量の比較で、系列が長くなるほど MASS の優位が指数的に開いていくことが分かります。長いセンサー波形ほど、この差は効いてきます。
距離プロファイルは「1つのクエリ窓」から見た景色でした。次は、これを全部の窓についてまとめあげて Matrix Profile を作ります。
距離プロファイルからMatrix Profileへ
距離プロファイルを1本のクエリについて求められるようになりました。Matrix Profile は、これをすべての窓 $i$ について計算し、各 $i$ で「自分以外の最小距離」だけを残したものです。式で書くと、$i$ 番目の Matrix Profile 値 $P[i]$ と、その最近傍がどこかを示す Profile Index $I[i]$ は次のようになります。
$$ \begin{equation} P[i] = \min_{\,|i-j|\ge w} D_i[j], \qquad I[i] = \arg\min_{\,|i-j|\ge w} D_i[j] \end{equation} $$
ここで $D_i$ は窓 $i$ をクエリにした距離プロファイルです。条件 $|i-j|\ge w$ は重要で、これは自明な一致(trivial match)を除外するためのものです。窓 $i$ のすぐ隣の窓 $i+1$ はほぼ同じ形なので距離が小さくなりますが、それは「似た部分列を見つけた」のではなく「自分自身とほぼ重なっている」だけ。これを除かないと Matrix Profile が全部ゼロ近くになってしまいます。そこで、窓長の $1/4$ 程度($w = \lceil m/4 \rceil$)を除外帯(exclusion zone)として無視します。
距離プロファイルを全 $i$ について順に MASS で計算し、最小を取れば Matrix Profile が得られます。これが STAMP(Scalable Time series Anytime Matrix Profile) です。
def stamp(T, m):
"""STAMP: 各窓のMASSを順に計算してMatrix ProfileとProfile Indexを得る。"""
n = len(T)
L = n - m + 1
mp = np.full(L, np.inf)
mpi = np.full(L, -1, dtype=int)
excl = int(np.ceil(m/4)) # 自明一致の除外帯
for i in range(L):
D = mass(T[i:i+m], T)
lo, hi = max(0, i-excl), min(L, i+excl+1)
D[lo:hi] = np.inf # 自分の近傍を除外
j = int(np.argmin(D))
if D[j] < mp[i]:
mp[i], mpi[i] = D[j], j
return mp, mpi
mp, mpi = stamp(T, m)
L = len(mp)
motif_i = int(np.argmin(mp)) # MP最小 = モチーフ
motif_j = int(mpi[motif_i])
discord_i = int(np.argmax(mp)) # MP最大 = ディスコード(異常)
print("モチーフ対 (i, j) =", motif_i, motif_j, " 距離 =", round(mp[motif_i], 3))
print("ディスコード i =", discord_i, " MP =", round(mp[discord_i], 3))
このコードを実行すると、モチーフ対は $(560, 150)$ で距離 0.296、ディスコードは $i=342$ で MP 9.832 と出ます。モチーフ対が、こちらが埋め込んだ2か所(150 と 560)をぴたりと言い当てているのが分かります。距離0.296という小ささは「ほぼ同一の形が2回現れた」ことを意味します。一方ディスコードは、異常バースト(380近辺)に重なる窓で、MP が背景の数倍に跳ね上がっています。
上が時系列、下が Matrix Profile です。MP の谷(緑の三角)が、2か所の同形モチーフの位置に対応し、ほぼゼロまで落ちています。MP の山(赤の三角、プラトーの立ち上がり)が、どこにも似ていない異常バーストの位置です。ベクトルを1本見るだけで、繰り返しパターンも異常も同時に読める——これが Matrix Profile の最大の魅力です。
ここで自然な疑問が湧きます。MASSで賢く計算したSTAMPは、本当に「全ペアを愚直に比べたナイーブ版」と同じ答えになっているのでしょうか。確かめておきましょう。
ナイーブ版との一致確認とSTOMP/SCRIMP++
STAMP が正しいことを保証するため、全ペアの z正規化距離を二重ループで愚直に計算するナイーブ版と突き合わせます。これは遅いですが、定義そのままなので「正解」とみなせます。
def naive_mp(T, m):
"""ナイーブ版: 全ペアのz正規化距離を二重ループで。STAMPとの一致確認用。"""
n = len(T)
L = n - m + 1
excl = int(np.ceil(m/4))
mp = np.full(L, np.inf)
for i in range(L):
for j in range(L):
if abs(i-j) < excl: # 自明一致を除外
continue
d = zdist(T[i:i+m], T[j:j+m])
if d < mp[i]:
mp[i] = d
return mp
mp_naive = naive_mp(T, m)
print("ナイーブ vs STAMP 最大誤差 :", float(np.max(np.abs(mp - mp_naive))))
実行すると最大誤差は 0.0。STAMP は MASS による高速化を行っても、ナイーブな全ペア比較と完全に同じ Matrix Profile を返します。近似ではなく厳密解である点が、Matrix Profile が実務で信頼される理由の一つです。
灰色の太線(ナイーブ)と赤の破線(STAMP)が完全に一致しています。高速化しても答えがぶれないことが、ひと目で確認できます。
ここまでが STAMP ですが、Matrix Profile には一連の改良アルゴリズムがあります。考え方の違いを押さえておきましょう。
- STAMP：各窓について独立に MASS を呼ぶ。全体で $O(n^2\log n)$。距離プロファイルを「どの順で」評価してもよいので、途中で打ち切っても近似解が得られる anytime 性が長所です。
- STOMP(Scalable Time series Ordered-search Matrix Profile)：隣り合う距離プロファイルの間でスライディング内積を漸化式で更新します。$\text{QT}_{i+1}[j+1] = \text{QT}_i[j] – t_i t_j + t_{i+m} t_{j+m}$ という $O(1)$ の更新で、FFT すら使わずに済みます。全体は $O(n^2)$ になり、$\log n$ の係数も消えて STAMP より高速です。
- SCRIMP++：STOMP の漸化式更新を、行列の対角線(diagonal)に沿って進めることで効率化し、さらに最初に粗い事前計算(preprocessing)で大きな谷を素早く見つける手法を組み合わせます。anytime 性と高速性を両立し、巨大データでも早い段階で良いモチーフ候補を返せます。
これらはすべて同じ Matrix Profile を計算するための高速化テクニックであり、出力の意味は変わりません。本記事では理解しやすい STAMP を実装の中心に据えていますが、実データで大規模に回す際は STOMP/SCRIMP++ が標準です。
漸化式更新の核を、STOMP の内積更新だけ抜き出して確認してみましょう。
def sliding_dot_update(T, m):
"""STOMPの核: 隣接窓のスライディング内積を漸化式でO(1)更新できることの確認。"""
n = len(T); L = n - m + 1
# 直接計算したQT(窓0と全窓の内積)を基準にする
QT0 = np.array([np.dot(T[0:m], T[j:j+m]) for j in range(L)])
# 漸化式: QT[i,j] -> QT[i+1,j+1] = QT[i,j] - T[i]T[j] + T[i+m]T[j+m]
QT_prev = QT0.copy()
for i in range(0, 3): # 数ステップだけ更新して検証
QT_next = np.full(L, np.nan)
for j in range(L-1):
QT_next[j+1] = QT_prev[j] - T[i]*T[j] + T[i+m]*T[j+m]
# 端(j=0)は直接計算で埋める
QT_next[0] = np.dot(T[i+1:i+1+m], T[0:m])
# 直接計算と一致するか
direct = np.array([np.dot(T[i+1:i+1+m], T[j:j+m]) for j in range(L)])
err = np.nanmax(np.abs(QT_next - direct))
QT_prev = direct # 次ステップの基準に直接値を使う
print("STOMP漸化式更新の最大誤差 :", float(err))
sliding_dot_update(T, m)
実行すると誤差はほぼ $0$(浮動小数点の丸め程度)で、隣接する距離プロファイルが1ステップ $O(1)$ の足し引きで更新できることが確認できます。これが、STOMP が FFT 無しで $O(n^2)$ を達成できる仕組みです。
Matrix Profile が手に入ったので、いよいよ本題の「モチーフ」と「ディスコード」を取り出します。
モチーフ発見：繰り返しパターンを見つける
モチーフ(motif)とは、時系列の中で繰り返し現れる典型的な部分列のことです。心拍の1拍、歩行の1歩、機械の1サイクル——「いつもの形」がモチーフです。Matrix Profile の上では、モチーフは話が驚くほど簡単になります。
Matrix Profile の最小値を取る位置 $i$ と、その Profile Index $I[i]=j$ が指す位置のペアが、最も強いモチーフ(top-1 motif)である。
理由は明快です。$P[i]$ が小さいということは「窓 $i$ には、自分以外にとてもよく似た窓 $j$ が存在する」ということ。その最小ペアこそ、系列全体で最も似た2つの部分列、すなわち最頻・最典型のパターンだからです。先ほどの実行結果では、モチーフ対は $(560, 150)$、距離 0.296 でした。この2つの窓を取り出して重ねてみます。
a = T[motif_i:motif_i+m]
b = T[motif_j:motif_j+m]
print("モチーフ対の生の距離 :", round(np.sqrt(((a-b)**2).sum()), 3))
print("モチーフ対のz正規化距離 :", round(zdist(a, b), 3)) # = mp[motif_i]
z正規化距離は0.296で、Matrix Profile の最小値と一致します。生の距離もごく小さく、2つの「二つこぶ」窓がほとんど同じ形であることを示しています。背景の高さやスケールが違う一般のケースでも、z正規化が形だけを取り出すので、同一パターンは確実に小さい距離になります。
左の生の部分列では、2つの窓は背景の高さの違いで少しずれて見えます。右のz正規化後では、2つの「二つこぶ」波形がぴたりと重なり、これらが同一のパターンであることが確認できます。Matrix Profile はこのペアを、ラベルも事前知識もなしに自動で見つけ出しました。
top-1 だけでなく、2番目・3番目のモチーフが欲しいときは、見つけたモチーフ近傍を除外して次の最小値を探す、という手順を繰り返します。多変量センサー時系列なら、各チャネルでモチーフを取り、よく一緒に出るパターンを束ねる、といった応用もできます。
モチーフが「最も似ているペア」なら、その正反対——最も似ていない、孤立した部分列は何を意味するでしょうか。それがディスコード、すなわち異常です。
ディスコード発見：異常部分列を見つける
ディスコード(discord)とは、系列全体を見渡しても、どこにも似た相手がいない最も孤立した部分列のことです。直感的に言えば「一度きりしか現れない、いつもと違う形」。これはまさに異常そのものです。Matrix Profile の上では、モチーフと鏡像の関係になります。
Matrix Profile の最大値を取る位置が、top-1 ディスコード(最も強い異常部分列)である。
$P[i]$ が大きいということは「窓 $i$ には、自分以外で似た窓が一つもない(最近傍までの距離が大きい)」ということ。つまり、その窓は系列の中で最も浮いた存在です。先ほどの結果では $i=342$、MP=9.832 でした。これは埋め込んだ高周波バースト(380近辺)に重なる窓で、背景のモチーフ(MP がほぼゼロ〜数程度)とは桁が違います。
print("ディスコード位置 i :", discord_i)
print("そのMP値 :", round(mp[discord_i], 3))
print("背景のMP中央値 :", round(float(np.median(mp)), 3))
ディスコードの MP は背景の中央値の何倍にもなり、突出していることが分かります。この窓を系列上で見てみましょう。
赤く塗った区間が、Matrix Profile が「最も孤立している」と判定した部分列です。背景の規則的な波ともモチーフのこぶとも似ない高周波バーストが、見事に異常として拾われています。しきい値も教師ラベルも使わず、ただ Matrix Profile の最大値を取るだけで異常部分列が見つかる——これが Matrix Profile による異常検知の手軽さです。
ディスコードによる異常検知が優れているのは、「正常がどんな形か」を一切モデル化しない点です。再構成ベースのオートエンコーダや予測ベースの手法は「正常を学んで、外れたものを異常とする」ので、訓練データに紛れたノイズや偏りに引きずられます。ディスコードは逆に「他のどの部分列とも似ていない」という相対的な孤立度だけを見るので、正常パターンが複数あっても、データの大半が正常でありさえすれば機能します。ラベルも訓練フェーズも要らず、Matrix Profile を1本計算するだけ。これが、現場で最初に試す異常検知としてディスコードが好まれる理由です。
実用では、Matrix Profile をそのまま点ごとの異常スコアに変換すると扱いやすくなります。各時刻について、それを含む窓の MP の最大値をスコアとし、しきい値を超えた区間を異常と判定します。
# 点ごとの異常スコア: その点を含む窓のMPの最大値
pt = np.zeros(len(T))
for i in range(L):
pt[i:i+m] = np.maximum(pt[i:i+m], mp[i])
med = np.median(pt)
thr = med + 0.6*(pt.max() - med) # 背景と異常の中間に閾値
det = np.where(pt > thr)[0]
print("検出区間 :", det.min(), "〜", det.max(), " 点数 =", len(det))
実行すると、検出区間は埋め込んだ異常(380〜460)を覆う形(窓幅 $m$ の分だけ前後に滲みます)で得られ、背景や正常なモチーフは一切引っかかりません。
下のスコアは、正常な区間では低く平らに保たれ、異常区間でだけ大きく跳ね上がります。閾値(橙の破線)を超えるのは異常部分列のところだけで、誤検出は出ていません。なお、検出区間が真の異常より少し広く出るのは、窓を使う手法に共通の性質です。異常を含む窓は、その窓が始まる位置すべてで MP が高くなるため、幅 $m$ だけ滲みます。
ここまで窓長 $m$ を固定して話してきましたが、この $m$ の選び方こそ Matrix Profile を使う上での最重要ハイパーパラメータです。最後にその効きを見ておきましょう。
窓長mの選び方
Matrix Profile で唯一きちんと決める必要があるのが、窓長 $m$ です。$m$ は「どのくらいの長さのパターンを探したいか」を表します。短すぎれば細かい揺れに反応し、長すぎれば大局的な形しか見なくなります。同じ系列でも $m$ を変えると、見える構造が変わります。
import matplotlib.pyplot as plt
for mm in [30, 80, 160]:
mpx, _ = stamp(T, mm)
di = int(np.argmax(mpx)) # ディスコード
mi = int(np.argmin(mpx)) # モチーフ
print(f"m={mm:3d}: モチーフ i={mi:3d} ディスコード i={di:3d} 最大MP={mpx[di]:.2f}")
この実行から、窓長によってモチーフ・ディスコードの位置や Matrix Profile の形が変わることが確認できます。ただし、本記事の合成データではどの $m$ でも異常バーストは最大値として拾われ、検出自体は安定しています。
上から $m=30, 80, 160$ の Matrix Profile です。短い窓($m=30$)では背景の細かい繰り返しまで谷として現れ、プロファイルがにぎやかになります。長い窓($m=160$)では大局的な形だけが残り、なめらかになります。中間の $m=80$ は、埋め込んだモチーフ(長さ80)とぴったり合っているため、モチーフの谷も異常の山も最もくっきり出ます。$m$ は「探したいパターンの典型的な長さ」に合わせるのが基本です。心拍なら1拍の長さ、機械なら1サイクルの長さ、というように、ドメイン知識から決めるのが定石です。
実務では複数の $m$ で Matrix Profile を計算し、どの解像度でも一貫して山になる区間を異常とみなす、といったマルチスケールの使い方もよく行われます。これで Matrix Profile の基本道具が一通りそろいました。
まとめ
本記事では、Matrix Profile を z正規化距離から MASS、STAMP、そしてモチーフ・ディスコード抽出まで、数式とスクラッチ実装で組み上げました。
- 3段重ね：部分列 → 距離プロファイル(1つの窓から全位置への距離) → Matrix Profile(各窓の最近傍距離だけを集めたベクトル)。1本のベクトルに全ペア比較を圧縮する。
- z正規化ユークリッド距離：平均と標準偏差をそろえ、オフセット・スケールに左右されず「形」だけを比べる。$d_{\text{z}}=\sqrt{2m(1-\rho)}$ と Pearson相関で書ける。
- MASS：距離が相関(=スライディング内積)で書けることを利用し、内積を FFT、平均・分散を累積和で求めて、距離プロファイルを $O(n\log n)$ に。STAMP/STOMP/SCRIMP++ はこれを全窓に広げる高速化で、いずれも厳密な同じ Matrix Profile を返す。
- モチーフとディスコード：MP の最小値が繰り返しパターン(モチーフ)、最大値が孤立した異常部分列(ディスコード)。しきい値も教師ラベルも要らない。
Matrix Profile の強みは、厳密・高速・パラメータが窓長 $m$ ほぼ一つという素直さです。教師なしでモチーフも異常も同時に読めるため、多変量センサー時系列の探索的な分析や、ラベルの付けにくい現場での異常検知の第一手として非常に有効です。
次のステップとして、Matrix Profile によるディスコード検出を、再構成ベースや予測ベースの手法と並べて位置づける以下の記事も参考にしてください。
参考文献
- C.-C. M. Yeh et al. “Matrix Profile I: All Pairs Similarity Joins for Time Series.” ICDM 2016.
- A. Mueen et al. “The Fastest Similarity Search Algorithm for Time Series Subsequences under Euclidean Distance”（MASS）.
- Y. Zhu et al. “Matrix Profile II”（STOMP）, “Matrix Profile XI: SCRIMP++”（SCRIMP++）.