数百万点のセンサーログがあるとします。手元には「この波形に似た区間を探したい」というクエリ波形が1つ。たとえば製造ラインの振動・電流・温度をまとめて記録した多変量時系列から、「過去に一度だけ起きた、あの特徴的なうねり」に似た区間を全部引き出したい——これが部分列類似検索（subsequence similarity search）です。
素朴にやるなら、長い系列の上をクエリと同じ長さの窓で1点ずつスライドさせ、すべての位置でユークリッド距離を計算して、近い順に並べればよいだけです。答えは確実に合っています。しかし系列が数百万点・チャネルが数十本になると、この全探索（brute force）は途方もない回数の距離計算を要し、現実的な時間で終わりません。しかも、「クエリの長さ」は事前に決まっていないことが多い。短いパターンも長いパターンも同じデータから探したいのです。
この「多変量・可変長の部分列 k-NN 検索を、答えを一切妥協せず（厳密に）、しかも高速に」という難題に答えるのが、本記事で扱う MULISSE（B. Pelok & J. E. d’Hondt, IEEE ICDEW 2025）です。MULISSE は、単変量時系列で最先端だった索引 ULISSE（M. Linardi & T. Palpanas, VLDB 2018）を多変量へ拡張したもので、心臓部は iSAX（記号表現）による索引 と 下界（lower bound）による厳密な枝刈り にあります。
この技術が活きる場面は広いです。
- 多変量センサー監視：振動・電流・温度などを束ねたログから、既知の異常パターンや関心パターンに似た区間を引き当てる
- 多誘導の生体信号：12誘導の心電図のように複数チャネルが連動する波形から、特定の波形が現れた区間を検索する
本記事の内容
- なぜ全探索ではダメで、なぜ「索引＋下界」で厳密性を保てるのか
- PAA → SAX → iSAX という記号化のパイプライン（スクラッチ実装）
- SAX下界がユークリッド距離を常に下回ることの実測（厳密性の核心）
- 簡易 iSAX 木による部分列 k-NN 検索の実装と、全探索と同じ答えを出しつつ距離計算を減らせる枝刈りの実測
- ULISSE の master series / Envelope による「可変長を1表現で覆う」仕組み
- MULISSE の多変量化（チャネルごと SAX → 結合）とチャネル考慮のノード分割
前提・関連記事
この記事を読む前に、以下を押さえておくと理解が深まります。
まず、MULISSE が何をするのかを1枚の絵で俯瞰しましょう。
左には、チャネル数 $d$・全長 $N$ の長い多変量センサー時系列と、「これに似た区間はどこ？」という可変長 $\ell$ のクエリ波形があります。MULISSE はこれを直接全探索しません。あらかじめ系列を iSAX 索引木（中央）に整理しておき、クエリが来たら下界を使って大半のノードを枝刈りし、残った少数の候補だけを厳密なユークリッド距離で照合します（右）。鍵は「枝刈りしても全探索と同じ答えが必ず返る」こと。この厳密性がどこから来るのかを、これから順に解きほぐしていきます。
なぜ「索引＋下界」で厳密性を保てるのか
全探索は確実ですが遅い。では索引で速くしたとき、なぜ答えが壊れないのでしょうか。ここがすべての出発点なので、直感から丁寧に押さえます。
イメージは「安い見積もりで足切りする」です。ある候補とクエリの本当の距離を計算するのはコストが高い。そこで、本当の距離より必ず小さい（か等しい） 安い見積もり——これを下界（lower bound, LB） と呼びます——をまず計算します。もし、いまある $k$ 番目の近傍までの距離（暫定の合格ライン）よりも、この安い見積もりですら大きいなら、その候補は本当の距離も当然それより大きいので、計算するまでもなく不合格と分かります。
この論法が成り立つ条件はただ一つ、「下界が本当の距離を決して超えない」ことです。下界が一度でも本当の距離を上回ると（過大評価すると）、本来は近いはずの候補を「遠い」と誤判定して捨ててしまい、答えが壊れます。逆に、下界が常に本当の距離以下であることさえ保証されていれば、足切りで捨てるのは「本当に遠い候補だけ」になり、全探索とまったく同じ k-NN が、より少ない距離計算で得られます。これが「索引＋下界」で厳密性（exactness）を保てる仕組みです。
ポイント：索引が速いのは、安い下界で大半の候補を計算前に捨てるから。厳密なのは、その下界が真の距離を決して超えないから。
もう少し具体的に流れを描くと、検索はおおむね次の3ステップになります。まず、長い系列をあらかじめ記号化して木に整理しておきます（索引構築、クエリより前の前処理）。次に、クエリが来たらノードごとに下界を計算し、近そうなノードから順に降りていき、見込みのないノードは丸ごと飛ばします（枝刈り）。最後に、生き残ったごく少数の候補に対してだけ本当のユークリッド距離を計算し、近い順に $k$ 個を返します（厳密照合）。重い計算（本物の距離）を最後の少数に集中させ、軽い計算（下界）で大半を片付けるのがコツです。下界が真の距離以下である限り、この3ステップは全探索の答えを一字一句変えずに返します。
では、その「安い下界」を時系列でどう作るのか。ここで登場するのが PAA → SAX → iSAX という記号化のパイプラインです。次節から1段ずつ実装していきます。
PAA：系列を粗く要約する
下界を作る第一歩は、長い系列を少ない数で粗く要約することです。粗い要約どうしの距離は計算が軽く、しかも「本物の距離の下界」になるように設計できます。その最も基本的な道具が PAA（Piecewise Aggregate Approximation, 区分平均近似） です。
PAA の発想はとても素直です。長さ $n$ の系列を $w$ 個の等しい区間に分け、各区間を平均値1つで代表させる。長さ $n$ のデータが $w$ 個の数に圧縮されます。式で書くと、$j$ 番目の区間（$j=1,\dots,w$）の PAA 係数 $\bar{x}_j$ は次の通りです。
$$ \begin{equation} \bar{x}_j = \frac{w}{n}\sum_{t=\frac{n}{w}(j-1)+1}^{\frac{n}{w}j} x_t \end{equation} $$
各区間の点をすべて足して点数で割る、つまり区間内の平均を取っているだけです。区間の幅 $n/w$ が大きいほど要約は粗く、$w$ を大きくすると要約は細かくなります。実装も数行で済みます。
import numpy as np
def znorm(x):
"""zノルム化(平均0・分散1)。波形の『形』だけを比べるため。"""
s = x.std()
return (x - x.mean()) / (s if s > 1e-9 else 1.0)
def paa(x, w):
"""長さnの系列をw区間に区分平均近似。nはwで割り切れる前提。"""
n = len(x); seg = n // w
return x[:seg * w].reshape(w, seg).mean(axis=1)
実際に1本の系列を PAA で要約してみます。
import matplotlib.pyplot as plt
rng = np.random.default_rng(7); n = 64
x = znorm(np.cumsum(rng.standard_normal(n)) + 3*np.sin(np.linspace(0, 6, n)))
p = paa(x, w=8); seg = n // 8
plt.figure(figsize=(9.5, 4))
plt.plot(np.arange(n), x, color="#94a3b8", lw=1.4, label="元の系列(長さ64)")
for i in range(8):
plt.hlines(p[i], i*seg, (i+1)*seg, color="#6366f1", lw=3)
plt.legend(); plt.show()
灰色の元系列（64点）が、紫色の8本の水平線（8区間の平均）に置き換わっています。細かなギザギザは消えますが、系列全体の「上がって下がって」という大局的な形は保たれています。重要なのは、PAA どうしのユークリッド距離が元の系列の距離の下界になることです。平均を取ると区間内のばらつきが消えるぶん、PAA間の距離は元の距離より小さくなる——この「縮む」性質こそが、安全な足切りに使える理由です。
なお、波形の絶対的な高さやスケールではなく「形」を比べるため、PAA の前に zノルム化（平均0・分散1に正規化）しておくのが定石です。これで「同じ形だが上下にずれた／拡大された波形」も近いと判定できます。
PAA で系列を $w$ 個の実数に圧縮できました。次は、この実数をさらに記号（離散シンボル）に置き換えます。それが SAX です。
SAX：PAA係数を記号に変える
PAA は実数の列でした。これを有限種類の記号に変えると、索引（ハッシュや木）に載せやすくなり、距離テーブルを事前計算して足切りを高速化できます。これを担うのが SAX（Symbolic Aggregate approXimation） です。
SAX のアイデアは「値の高さで記号に振り分ける」こと。zノルム済みの系列は近似的に標準正規分布 $N(0,1)$ に従うので、その分布を等確率に$a$ 個の帯に分割します（$a$ は記号の種類数＝基数）。各帯に記号 a, b, c, …を割り当て、PAA 係数がどの帯に落ちるかでその区間を記号化します。帯の境目（ブレークポイント）$\beta_i$ は、$N(0,1)$ の累積分布が $i/a$ になる点、すなわち分位点で決めます。
$$ \begin{equation} \beta_i = \Phi^{-1}\!\left(\frac{i}{a}\right), \quad i = 1, \dots, a-1 \end{equation} $$
ここで $\Phi^{-1}$ は標準正規分布の逆累積分布関数（分位点関数）です。等確率に分けるのは、どの記号も同じくらいの頻度で現れ、情報を均等に使うためです。実装は分位点を計算して searchsorted で帯を引くだけです。
from scipy.stats import norm
def sax_breakpoints(a):
"""等確率ブレークポイント(N(0,1)をa等分する内部境界 a-1個)。"""
return norm.ppf(np.linspace(0, 1, a+1)[1:-1])
def sax(x, w, a):
"""zノルム→PAA→各係数を記号(0..a-1)に。"""
bp = sax_breakpoints(a)
return np.searchsorted(bp, paa(znorm(x), w))
先ほどの系列を基数 $a=4$（記号 a, b, c, d）で SAX に変換してみます。
a = 4; bp = sax_breakpoints(a); syms = sax(x, 8, a)
print("SAX語:", "".join("abcd"[s] for s in syms))
橙色の破線が等確率ブレークポイントで、平面を4つの帯に分けています。各 PAA 係数（紫の水平線）がどの帯にあるかで記号が決まり、系列全体が短いSAX語（記号列）になります。これで時系列が「単語」になりました。値が近い PAA どうしは同じ帯＝同じ記号になりやすいので、SAX語が一致・近接する系列は「形が似ている」候補だと、文字列の比較だけで素早く絞り込めます。
記号化のもう一つの利点は、索引との相性です。実数の列は「どこで区切ってグループ化するか」が自明ではありませんが、記号列なら「同じ記号で始まる系列」をまとめる、といった操作が文字列の前方一致だけで済みます。これが次節以降の iSAX 木（記号の接頭辞でノードを分ける）につながります。基数 $a$ を上げれば記号が増えて表現は精密になりますが、そのぶん同じ記号に落ちる系列が減って絞り込みは細かくなる——この粗さと精密さのトレードオフを、解像度として木の深さで調整するのが iSAX の発想です。
ただし SAX語を比べて「似ている／違う」と言うだけでは、足切りの数値（下界）になりません。次に、SAX語どうしからユークリッド距離の下界を計算する MINDIST を導入します。ここが厳密性の本丸です。
SAX下界（MINDIST）：枝刈りが厳密性を壊さない核心
足切りに使うには、SAX語どうしから「本当のユークリッド距離より必ず小さい数」を出さなければなりません。これを与えるのが MINDIST です。
直感はこうです。2つの系列の同じ区間が、たとえば記号 a と 記号 d のように離れた帯にあるなら、それぞれの値は少なくとも「a の帯の上端」と「d の帯の下端」だけ離れています。この帯と帯の隙間を、その区間が必ず持つ距離の下限として使います。逆に、記号が同じか隣り合う帯なら、値が境界付近でほぼ重なっているかもしれないので、下限は安全側に 0 とみなします。区間ごとのこの下限を二乗して足し、長さで補正したものが MINDIST です。
$$ \begin{equation} \mathrm{MINDIST}(\hat{Q}, \hat{C}) = \sqrt{\frac{n}{w}} \sqrt{\sum_{j=1}^{w} \big(\mathrm{dist}(\hat{q}_j, \hat{c}_j)\big)^2} \end{equation} $$
ここで $\hat{q}_j, \hat{c}_j$ は区間 $j$ の SAX 記号、$\mathrm{dist}$ は記号間のセル距離（隣接以下なら 0、離れていれば帯の隙間）です。各部分の意味を分解します。
- $\mathrm{dist}(\hat{q}_j, \hat{c}_j)$：区間 $j$ で「2つの値が最低でもこれだけ離れている」という下限。記号が隣接（$|q_j – c_j| \le 1$）なら 0、それ以外は離れた側の帯境界どうしの差。
- $\sqrt{\sum_j (\cdot)^2}$：各区間の下限を、ユークリッド距離と同じく二乗和の平方根でまとめる。
- $\sqrt{n/w}$：PAA で $n$ 点を $w$ 区間に圧縮したぶんを、元の長さ $n$ のスケールに戻す補正係数。
このセル距離テーブルと MINDIST を実装します。
def cell_dist_table(a):
"""記号間セル距離。隣接以下は0、離れた記号は帯境界の隙間。"""
bp = sax_breakpoints(a)
ext = np.concatenate([[-np.inf], bp, [np.inf]])
T = np.zeros((a, a))
for r in range(a):
for c in range(a):
if abs(r - c) <= 1:
T[r, c] = 0.0
else:
hi = ext[max(r, c)] # 大きい記号の下側境界
lo = ext[min(r, c) + 1] # 小さい記号の上側境界
T[r, c] = hi - lo
return T
def mindist(saxQ, saxC, n, w, a, T):
"""SAX下界(ユークリッド距離の下界)。"""
d2 = sum(T[saxQ[j], saxC[j]]**2 for j in range(w))
return np.sqrt(n / w) * np.sqrt(d2)
この MINDIST が「本当に下界か（真の距離を超えないか）」を、ランダムな系列対 2500 組で実測してみます。1組でも上回れば厳密性が崩れるので、ここは妥協できません。
rng = np.random.default_rng(11); n, w, a = 128, 16, 8
T = cell_dist_table(a); true_d = []; lb_d = []
for _ in range(2500):
q = znorm(np.cumsum(rng.standard_normal(n)))
c = znorm(np.cumsum(rng.standard_normal(n)))
true_d.append(np.sqrt(((q - c)**2).sum()))
sq = np.searchsorted(sax_breakpoints(a), paa(q, w))
sc = np.searchsorted(sax_breakpoints(a), paa(c, w))
lb_d.append(mindist(sq, sc, n, w, a, T))
true_d, lb_d = np.array(true_d), np.array(lb_d)
print("違反件数:", int((lb_d > true_d + 1e-9).sum()), "/", len(true_d))
print("平均 LB/真:", round((lb_d / true_d).mean(), 3))
横軸が真のユークリッド距離、縦軸が SAX下界です。すべての点が対角線 $y=x$ より下にあり、違反は 0/2500 件でした。下界の平均は真距離の約 0.61 倍で、「真の距離より小さい（安い見積もり）」かつ「ゼロに潰れず識別力がある（離れた対ほど下界も大きい）」という、足切りに理想的な性質を満たしています。これが「枝刈りしても真の最近傍を捨てない＝厳密」の数値的な根拠です。下界が一度も真距離を超えないので、合格ラインを超えた下界の候補は安心して計算前に捨てられます。
下界が手に入ったので、あとはこれを使って「どの候補から見て、どこで打ち切るか」を効率化する索引が要ります。ここで SAX を可変解像度に拡張した iSAX が登場します。
iSAX：可変解像度で木にする
SAX は基数 $a$ を固定して記号化しました。しかし索引木では、根に近いノードは粗く（少ない記号で大雑把に）、葉に近づくほど細かく（多くの記号で精密に）分けたい。同じ系列を解像度を変えて見られるようにしたのが iSAX（indexable SAX） です。
仕組みは「基数を 2 の冪で増やす」ことです。基数 $a=2$ なら各区間は 1 ビット（0 か 1）、$a=4$ なら 2 ビット、$a=8$ なら 3 ビット……と、ビット数を増やすほど帯が細かくなります。しかも下位ビットを足すだけなので、粗い記号は細かい記号の接頭辞になっていて、解像度の異なる SAX 表現を矛盾なく同じ木に載せられます。
x4 = x # 同じ系列を基数2,4,8で見る
for a in [2, 4, 8]:
bp = sax_breakpoints(a)
syms = np.searchsorted(bp, paa(x4, 4))
bits = int(np.log2(a))
print(f"基数{a}:", [format(s, f"0{bits}b") for s in syms])
同じ PAA 係数（4区間）を、基数 $a=2, 4, 8$ で見た様子です。基数を上げると帯（橙破線）が細かくなり、記号のビット数が増えます。木の浅いノードでは粗い基数で大雑把にグループ分けし、葉に下りるほど基数を上げて分解能を高める——これが iSAX による可変解像度索引です。あるノードに系列が溜まりすぎたら、ある区間の基数を 1 段上げて記号を細分化し、そのノードを2つに分割します。この「どの区間を細分化するか」が後で出てくるノード分割戦略です。
ここまでで「固定長」の部分列を索引する道具は揃いました。しかし、現実のクエリは長さが事前に分かりません。ULISSE が解いたのは、まさにこの「1つの索引で、ある長さの範囲のクエリすべてに答える」問題です。
ULISSE の核心：master series と Envelope で可変長を覆う
クエリ長が 64 かもしれないし 96 かもしれない。長さごとに索引を作り直すのは非現実的です。ULISSE は、ある長さ範囲 $[\ell_{\min}, \ell_{\max}]$ の部分列をまとめて1つの表現で覆うことで、これを解決しました。
鍵になるのが master series（マスター系列） と Envelope（包絡） です。発想はこうです。ある開始位置から、長さ $\ell_{\min}$ の部分列、$\ell_{\min}+1$ の部分列、……$\ell_{\max}$ の部分列を取り出すと、それぞれを zノルム→PAA した係数は少しずつ違います。これらをひとまとめにして、各 PAA 区間ごとに「取りうる最小値 $L$ と最大値 $U$」で挟んだ帯を作ります。これが Envelope（含有領域）で、その元になる長さ範囲ぶんの系列が master series です。
$$ \begin{equation} \mathrm{paaENV}[D, \ell_{\min}, \ell_{\max}] = [\,\mathbf{L}, \mathbf{U}\,], \quad L_j = \min_{\ell} \bar{x}_j^{(\ell)}, \;\; U_j = \max_{\ell} \bar{x}_j^{(\ell)} \end{equation} $$
つまり Envelope は、その範囲のどの長さの部分列の PAA も必ずこの帯の中に入ることを保証する「外箱」です。クエリがこの範囲のどんな長さでも、外箱との距離はクエリと中身の真の距離の下界になります。実際に作ってみましょう。
rng = np.random.default_rng(5); Lmax = 96
base = np.cumsum(0.15*rng.standard_normal(Lmax)) + np.sin(np.linspace(0, 5, Lmax))
lmin, lmax, w = 48, 96, 8
segs = []
for L in range(lmin, lmax+1, 2):
seg = L // w
segs.append(paa(znorm(base[:seg*w]), w))
segs = np.array(segs)
U, Lo = segs.max(axis=0), segs.min(axis=0) # 包絡の上端・下端
灰色の細い線が、長さ 48〜96 の各部分列を zノルム→PAA したものです。それらを各区間で min/max で挟んだ水色の帯が Envelope（含有領域）です。この帯1つが、48〜96 のどの長さのクエリに対しても下界を与えられるので、長さごとに索引を作る必要がなくなります。ULISSE はこの Envelope を iSAX で記号化して木に格納し、可変長クエリに単一索引で答えます（ULISSE は “ULtra compact index for variable-length similarity SEarch” の略で、この圧縮の効きが名前の由来です）。
ここまでが単変量の ULISSE です。いよいよ本題の MULISSE——これを多変量へ拡張するとき、何が課題で、どう解くのかを見ていきます。
MULISSE：多変量への拡張
チャネルが $d$ 本になると、何が難しくなるのでしょうか。素朴には「全チャネルを連結して1本の長いベクトルにすればいい」と思えますが、それでは次元が $d$ 倍に膨れ、記号表現が粗くなって下界が緩み（足切りが効かなくなり）、索引の利点が消えてしまいます。下界が緩むと、本来なら捨てられたはずの候補まで「念のため計算」せざるを得なくなり、枝刈り率が落ちて全探索に近づいてしまうのです。さらに、チャネルごとに値の暴れ方やスケールが違うのに、連結して一律に記号化すると、おとなしいチャネルに記号の解像度を取られて肝心のチャネルの差が潰れる、という問題も起きます。MULISSE の工夫は、チャネルごとに SAX で記号化し、それを結合して多変量の記号表現とする点にあります。チャネルごとに独立に zノルムして記号化すれば、各チャネルの「形」が公平に表現され、下界の識別力を保ったまま多変量に対応できます。
多変量SAX：チャネルごと記号化 → 結合
各チャネル $c$ の部分列を独立に zノルム→PAA→SAX し、得られた SAX語をチャネル方向に並べて1つの多変量記号表現にします。距離（下界）は、チャネルごとの MINDIST を二乗和して平方根を取ったものになります。
$$ \begin{equation} \mathrm{MINDIST}_{\text{MV}}(\hat{Q}, \hat{C}) = \sqrt{\sum_{c=1}^{d} \mathrm{MINDIST}(\hat{Q}^{(c)}, \hat{C}^{(c)})^2} \end{equation} $$
これは、多変量ユークリッド距離が「各チャネルの距離の二乗和の平方根」であることに対応します。各チャネルの MINDIST が真距離の下界なら、その二乗和の平方根も多変量の真距離の下界になります。チャネルごとに記号化しても厳密性が保たれるのか、$d=3$ で実測します。
rng = np.random.default_rng(21); n, w, a, d = 128, 16, 8, 3
T = cell_dist_table(a); true_d = []; lb_mv = []
for _ in range(2000):
dd2 = 0.0; lbsum2 = 0.0
for _c in range(d):
q = znorm(np.cumsum(rng.standard_normal(n)))
c = znorm(np.cumsum(rng.standard_normal(n)))
dd2 += ((q - c)**2).sum()
sq = np.searchsorted(sax_breakpoints(a), paa(q, w))
sc = np.searchsorted(sax_breakpoints(a), paa(c, w))
lbsum2 += mindist(sq, sc, n, w, a, T)**2
true_d.append(np.sqrt(dd2)); lb_mv.append(np.sqrt(lbsum2))
true_d, lb_mv = np.array(true_d), np.array(lb_mv)
print("多変量下界の違反:", int((lb_mv > true_d + 1e-9).sum()), "/", len(true_d))
3チャネルの場合でも、多変量SAX下界はすべて対角線より下にあり、違反は 0 件でした。チャネルごとに別々に記号化して結合しても、下界の性質（真距離を超えない）はチャネル方向の二乗和の中で保たれます。つまり MULISSE は、多変量でも全探索と同じ厳密な答えを保証できるのです。連結して粗くするのではなく、チャネルごとに精度よく記号化してから束ねる——これが次元増加に対する MULISSE の答えです。
チャネル考慮のノード分割
iSAX 木では、ノードが満杯になると「ある区間の基数を上げて2分割」します。単変量なら「どの PAA 区間を細分化するか」だけを選べばよかったのですが、多変量では「どのチャネルの・どの区間を細分化するか」という選択肢が $d$ 倍に増えます。ここで適当に選ぶと、下界がほとんど締まらない無駄な分割をしてしまいます。
MULISSE は、含有領域（その記号が覆う値の幅）が最も広いチャネル×区間を優先して分割します。広い帯を細分化すると下界が最も大きく締まり、足切りが効くようになるからです。模式的に示します。
ヒートマップは「チャネル×PAA区間」ごとの含有領域の広さです。MULISSE は最も広いセル（青枠＝チャネル1の区間3）の基数を上げて分割します。そこを細かくすると下界が最も締まり、以降の枝刈りが効くようになります。「下界が最も締まる方向に木を成長させる」 という指針が、多変量での索引効率を支えています。チャネルごとに重要度（暴れ方）が違う現実のデータで、これは効いてきます。
なぜ「広いところから分割する」のが得なのかも、下界の式に戻ると腑に落ちます。含有領域が広い記号は、その区間で「実は値がもっと近いかもしれない」という曖昧さを大きく抱えており、安全側に倒すぶん下界が小さく（緩く）出ます。そこを2つに割って帯を半分にすれば、曖昧さが減って下界が大きく（厳しく）なり、足切りラインを越えやすくなります。逆に、もともと狭い（値がほぼ確定している）区間をいくら割っても下界はほとんど変わらず、木だけ無駄に深くなります。限られた分割回数を、最も効果の大きいチャネル×区間に投資する——これがチャネル考慮の分割戦略の狙いです。
理論が揃ったので、最後に「実際に索引で k-NN を引くと、全探索と同じ答えが、より少ない距離計算で得られる」ことを通しで実証します。
実装：簡易 iSAX 索引で部分列 k-NN を厳密に解く
ここまでの部品（PAA・SAX・MINDIST）を組み合わせ、簡易 iSAX 索引による部分列 k-NN 検索を実装します。狙いは2つ。(1) 全探索と完全に同じ答えを出すこと（厳密性）、(2) 距離計算の回数を減らせること（枝刈り）。
方針はこうです。長い系列の全部分列を SAX語でバケットにまとめ（同じ SAX語＝同じ葉に相当）、クエリの SAX語との MINDIST が小さいバケットから順に訪問します。best-first で近そうな候補から厳密距離を計算し、暫定の $k$ 番目距離より「これから訪れるバケットの下界の方が大きく」なったら、残りのバケットを丸ごと枝刈りします（下界が合格ライン以上なので中身も全部不合格）。
import heapq
def knn_search(longx, query, w=8, a=8, k=5):
"""簡易iSAX索引による厳密 部分列k-NN。返り値:(上位idx, 厳密距離計算回数, 全候補数)"""
n = len(query); L = len(longx)
bp = sax_breakpoints(a); T = cell_dist_table(a)
qz = znorm(query); sq = np.searchsorted(bp, paa(qz, w))
# 全部分列をSAX語でバケット化(=葉に相当)
buckets = {}; subz = np.empty((L-n+1, n))
for i in range(L-n+1):
z = znorm(longx[i:i+n]); subz[i] = z
sw = tuple(np.searchsorted(bp, paa(z, w)))
buckets.setdefault(sw, []).append(i)
# 各バケットの下界で昇順に訪問(best-first)
blist = sorted(((mindist(sq, np.array(sw), n, w, a, T), idxs)
for sw, idxs in buckets.items()), key=lambda r: r[0])
knn = []; exact = 0
for lb, idxs in blist:
if len(knn) == k and lb > -knn[0][0]:
break # 残りは下界が現状kth以上→全部枝刈り
for i in idxs:
d = np.sqrt(((qz - subz[i])**2).sum()); exact += 1
if len(knn) < k: heapq.heappush(knn, (-d, i))
elif d < -knn[0][0]: heapq.heapreplace(knn, (-d, i))
res = sorted([(-nd, i) for nd, i in knn])
return [i for _, i in res], exact, (L-n+1)
これを全探索と突き合わせます。8000 点の系列から長さ 80 のクエリで 5-NN を引きます。
def build_long_series(N=8000, seed=1):
rng = np.random.default_rng(seed); t = np.arange(N)
return (np.sin(2*np.pi*t/50) + 0.5*np.sin(2*np.pi*t/17)
+ 0.25*rng.standard_normal(N))
def brute_force(longx, query, k=5):
n = len(query); qz = znorm(query); ds = []
for i in range(len(longx)-n+1):
ds.append((np.sqrt(((qz - znorm(longx[i:i+n]))**2).sum()), i))
return sorted(i for _, i in sorted(ds)[:k])
longx = build_long_series(8000, 1); n = 80
rng = np.random.default_rng(99)
query = longx[3000:3000+n] + 0.1*rng.standard_normal(n)
bf = brute_force(longx, query, 5)
idx, exact, total = knn_search(longx, query, k=5)
print("全探索 :", bf)
print("iSAX索引:", sorted(idx), " 一致:", set(idx) == set(bf))
print(f"厳密距離計算: {exact}/{total} 回 枝刈り率: {1-exact/total:.1%}")
結果は明快です。iSAX 索引の答えは全探索と完全一致（どちらも同じ5区間）でありながら、厳密距離計算は 7,921 回中 1,911 回——約76%の距離計算を枝刈りで省けました。上段の緑帯がヒット区間、赤帯がクエリ由来の位置で、似た波形がきちんと引けています。下段の棒が示す通り、答えの正しさを1ミリも犠牲にせず、計算量だけを削減できています。これが「索引＋下界」の厳密かつ高速な検索の威力です。
最後に、MULISSE の眼目である可変長で枝刈りがどう効くかを確認します。クエリ長 $\ell$ を変えながら枝刈り率を測ります。
longx = build_long_series(6000, 4); rng = np.random.default_rng(55)
for n in [40, 60, 80, 120, 160]:
q = longx[2500:2500+n] + 0.1*rng.standard_normal(n)
_, exact, total = knn_search(longx, q, k=5)
print(f"クエリ長{n:3d}: 枝刈り率 {1-exact/total:.0%}")
クエリ長を 40〜160 と変えても、枝刈り率は 60〜82% を保ちました。短いクエリほど SAX 表現の識別力が相対的に高く枝刈りが強く効き、長いクエリでも十分な削減が得られています。どの長さでも厳密性（全探索一致）を保ったまま距離計算を大幅に削減できる——これが ULISSE/MULISSE の「可変長を1つの索引で」という設計の実利です。実際の MULISSE は、この簡易バケットの代わりに Envelope 付き iSAX 木とチャネル考慮の分割を用い、より大規模・多変量で同じ厳密性と枝刈りを実現します。
ULISSE から MULISSE へ：何が変わったか
ここまでの流れを、系譜として1枚にまとめます。
ULISSE（VLDB 2018）は、PAA/SAX/iSAX に master series の Envelope を組み合わせ、単変量・可変長の部分列検索を厳密かつ高速に解く最先端でした。MULISSE（ICDEW 2025）はこれを多変量へ拡張し、(1) チャネルごとに記号化して結合する多変量SAX、(2) チャネル別下界の二乗和を全体下界とする厳密性の維持、(3) 下界が最も締まる方向に木を成長させるチャネル考慮のノード分割、を加えました。引き継いだのは「下界による厳密な枝刈り」という背骨、拡張したのは「複数チャネルをどう記号化し、どう分割するか」という多変量特有の設計です。
なお、部分列の「形」を比べる発想は、距離プロファイルを全位置で一括計算する Matrix Profile（slug: matrix-profile）とも近い親戚です。Matrix Profile が「すべての部分列ペアの最近傍距離を求める」のに対し、MULISSE は「与えたクエリに対する k-NN を索引で厳密かつ可変長で引く」点で目的が異なります。用途に応じて使い分けるとよいでしょう。
まとめ
本記事では、多変量時系列の可変長・部分列類似検索を厳密かつ高速に解く MULISSE を、その土台 ULISSE とともに、スクラッチ実装で核心まで追いました。
- 問題：長い多変量センサー時系列から「クエリ波形に似た区間」を、可変長で・厳密に・高速に引きたい
- 索引＋下界：本当の距離より必ず小さい下界（MINDIST）で大半の候補を計算前に枝刈り。下界が真距離を超えないから、全探索と同じ答えが返る（厳密性）
- PAA→SAX→iSAX：系列を粗く要約→記号化→可変解像度の木に。SAX下界が真のユークリッド距離を常に下回ることを 0/2500 違反で実測
- ULISSE：master series の Envelope で、ある長さ範囲のクエリすべてを1つの索引で覆う（可変長対応）
- MULISSE：チャネルごと記号化→結合の多変量SAX（多変量でも下界の違反 0 件）＋チャネル考慮のノード分割で、多変量へ厳密に拡張
- 実測：簡易 iSAX 索引で 5-NN を引くと、全探索と完全一致しつつ距離計算を約76%削減。可変長でも 60〜82% の枝刈りを維持
「安い見積もりで足切りし、残りだけ厳密に照合する」というシンプルな原理が、数百万点・多チャネルの時系列検索を現実的な速度に引き下げます。次のステップとして、ここで引いた部分列を異常検知やパターン発見にどうつなげるか、関連記事も参考にしてください。
参考文献
- B. Pelok & J. E. d’Hondt. “MULISSE: Variable-Length Similarity Search for Multivariate Time Series.” 2025 IEEE 41st International Conference on Data Engineering Workshops (ICDEW), 2025.
- M. Linardi & T. Palpanas. “Scalable, Variable-Length Similarity Search in Data Series: The ULISSE Approach.” PVLDB 11(13), 2236–2248, 2018.
- J. Lin, E. Keogh, L. Wei, S. Lonardi. “Experiencing SAX: a novel symbolic representation of time series.” Data Mining and Knowledge Discovery, 15(2), 2007（SAX/MINDIST の原典）.