TimeSiam — Siameseネットワークで時系列の時間的相関を学習する【ICML 2024】

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# --- 周期的な合成テレメトリデータの生成 ---
def generate_telemetry(n_series=20, length=512, period=90):
    """軌道周期を持つ衛星テレメトリを模擬"""
    t = np.arange(length)
    series_list = []
    for i in range(n_series):
        base_freq = 2 * np.pi / period
        amplitude = 1.0 + 0.2 * np.random.randn()
        phase = np.random.uniform(0, 2 * np.pi)
        trend = 0.001 * np.random.randn() * t
        noise = np.random.randn(length) * 0.1
        s = amplitude * np.sin(base_freq * t + phase) + trend + noise
        # 一部の系列に異常を挿入
        if i % 5 == 0:
            anomaly_start = np.random.randint(200, 400)
            s[anomaly_start:anomaly_start+20] += 3.0
        series_list.append(s)
    return np.array(series_list)

telemetry = generate_telemetry(n_series=30, length=512, period=90)

# --- パッチの抽出 ---
def extract_patch_pairs(series, patch_size=32, delta_range=(10, 100)):
    """past-as-targetでパッチペアを生成"""
    pairs = []
    for s in series:
        for t in range(delta_range[1] + patch_size, len(s) - patch_size):
            delta = np.random.randint(delta_range[0], delta_range[1])
            current = s[t:t+patch_size]
            past = s[t-delta:t-delta+patch_size]
            pairs.append((current, past, delta))
    return pairs

pairs = extract_patch_pairs(telemetry, patch_size=32, delta_range=(10, 100))
print(f"生成されたパッチペア数: {len(pairs)}")

# サンプリングして学習データにする
np.random.shuffle(pairs)
pairs = pairs[:2000]

X_current = np.array([p[0] for p in pairs])
X_past = np.array([p[1] for p in pairs])
deltas = np.array([p[2] for p in pairs])

# --- 簡易Siameseエンコーダ ---
class SimpleSiamese:
    def __init__(self, input_dim=32, latent_dim=16, lr=0.01, ema_decay=0.99):
        self.latent_dim = latent_dim
        self.lr = lr
        self.ema_decay = ema_decay

        # オンラインエンコーダ
        self.W_online = np.random.randn(input_dim, latent_dim) * 0.05
        self.b_online = np.zeros(latent_dim)

        # ターゲットエンコーダ（EMA）
        self.W_target = self.W_online.copy()
        self.b_target = self.b_online.copy()

        # 射影ヘッド
        self.W_proj = np.random.randn(latent_dim, latent_dim) * 0.05
        self.b_proj = np.zeros(latent_dim)

        # Lineage embedding（時間ラグの埋め込み）
        self.lineage_scale = np.random.randn(latent_dim) * 0.01

    def encode_online(self, x):
        z = np.tanh(x @ self.W_online + self.b_online)
        return z

    def encode_target(self, x, delta):
        z = np.tanh(x @ self.W_target + self.b_target)
        # Lineage embedding: 正弦波ベース
        lineage = np.sin(delta[:, None] * np.arange(self.latent_dim)[None, :] * 0.01)
        z = z + 0.1 * lineage
        return z

    def project(self, z):
        return np.tanh(z @ self.W_proj + self.b_proj)

    def cosine_sim(self, a, b):
        a_norm = a / (np.linalg.norm(a, axis=1, keepdims=True) + 1e-8)
        b_norm = b / (np.linalg.norm(b, axis=1, keepdims=True) + 1e-8)
        return np.sum(a_norm * b_norm, axis=1)

    def train_step(self, x_current, x_past, deltas_batch):
        # 順伝播
        z_online = self.encode_online(x_current)
        z_target = self.encode_target(x_past, deltas_batch)
        pred = self.project(z_online)

        # 損失: -cosine_similarity
        sim = self.cosine_sim(pred, z_target)
        loss = -np.mean(sim)

        # 簡易的な勾配更新
        pred_norm = pred / (np.linalg.norm(pred, axis=1, keepdims=True) + 1e-8)
        tgt_norm = z_target / (np.linalg.norm(z_target, axis=1, keepdims=True) + 1e-8)

        # 射影ヘッドの勾配
        grad_pred = -tgt_norm / len(x_current)
        dtanh_proj = 1 - pred ** 2
        grad_proj_input = grad_pred * dtanh_proj
        grad_W_proj = z_online.T @ grad_proj_input
        self.W_proj -= self.lr * grad_W_proj
        self.b_proj -= self.lr * np.mean(grad_proj_input, axis=0)

        # オンラインエンコーダの勾配
        grad_z_online = grad_proj_input @ self.W_proj.T
        dtanh = 1 - z_online ** 2
        grad_z_online *= dtanh
        grad_W_online = x_current.T @ grad_z_online
        self.W_online -= self.lr * grad_W_online
        self.b_online -= self.lr * np.mean(grad_z_online, axis=0)

        # ターゲットエンコーダのEMA更新
        self.W_target = self.ema_decay * self.W_target + (1 - self.ema_decay) * self.W_online
        self.b_target = self.ema_decay * self.b_target + (1 - self.ema_decay) * self.b_online

        return loss

# --- 学習ループ ---
model = SimpleSiamese(input_dim=32, latent_dim=16, lr=0.003)
n_epochs = 80
batch_size = 128
losses = []

for epoch in range(n_epochs):
    perm = np.random.permutation(len(X_current))
    epoch_loss = 0
    n_batches = 0
    for i in range(0, len(X_current), batch_size):
        idx = perm[i:i+batch_size]
        if len(idx) < 2:
            continue
        loss = model.train_step(X_current[idx], X_past[idx], deltas[idx])
        epoch_loss += loss
        n_batches += 1
    losses.append(epoch_loss / max(n_batches, 1))

# --- 学習曲線 ---
plt.figure(figsize=(8, 4))
plt.plot(losses, color='#00d4ff')
plt.xlabel('Epoch')
plt.ylabel('Loss (negative cosine similarity)')
plt.title('TimeSiam Training Loss')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('timesiam_loss.png', dpi=150, bbox_inches='tight')
plt.show()

# --- 時間ラグと埋め込み距離の関係 ---
# テスト用パッチペアを様々なΔで生成
test_series = telemetry[0]
test_deltas = list(range(5, 200, 5))
embeddings_current = []
embeddings_past = []

patch_size = 32
t_center = 300

current_patch = test_series[t_center:t_center+patch_size].reshape(1, -1)
z_current = model.encode_online(current_patch)

distances = []
for d in test_deltas:
    past_patch = test_series[t_center-d:t_center-d+patch_size].reshape(1, -1)
    z_past = model.encode_target(past_patch, np.array([d]))
    cos_sim = model.cosine_sim(z_current, z_past)[0]
    distances.append(cos_sim)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# (a) 時間ラグ vs コサイン類似度
axes[0].plot(test_deltas, distances, 'o-', color='#00d4ff', markersize=3)
axes[0].set_xlabel('Time lag Δ')
axes[0].set_ylabel('Cosine similarity')
axes[0].set_title('Embedding Similarity vs Time Lag')
axes[0].axvline(x=90, color='#ff6b6b', linestyle='--', label='Orbital period (90)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# (b) 元の時系列で対応区間を表示
t_plot = np.arange(len(test_series))
axes[1].plot(t_plot, test_series, color='gray', alpha=0.5, label='Full series')
axes[1].axvspan(t_center, t_center+patch_size, alpha=0.3, color='#00d4ff', label='Current patch')
for d in [30, 90, 180]:
    axes[1].axvspan(t_center-d, t_center-d+patch_size, alpha=0.2,
                     color='#ff6b6b' if d == 90 else '#ffd93d')
    axes[1].annotate(f'Δ={d}', xy=(t_center-d, test_series[t_center-d]),
                     fontsize=8, color='red')
axes[1].set_xlabel('Time step')
axes[1].set_ylabel('Value')
axes[1].set_title('Telemetry Series with Patch Locations')
axes[1].legend(fontsize=8)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('timesiam_lag_analysis.png', dpi=150, bbox_inches='tight')
plt.show()

# --- 全パッチの埋め込みを計算 ---
all_patches = []
all_positions = []
all_series_ids = []

for sid in range(min(10, len(telemetry))):
    s = telemetry[sid]
    for t in range(0, len(s) - patch_size, patch_size // 2):
        patch = s[t:t+patch_size]
        all_patches.append(patch)
        all_positions.append(t)
        all_series_ids.append(sid)

all_patches = np.array(all_patches)
all_positions = np.array(all_positions)
all_series_ids = np.array(all_series_ids)

Z = model.encode_online(all_patches)

# 簡易t-SNE（PCAで代用）
from numpy.linalg import svd
Z_centered = Z - Z.mean(axis=0)
U, S, Vt = svd(Z_centered, full_matrices=False)
Z_2d = Z_centered @ Vt[:2].T

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 系列IDで色分け
scatter = axes[0].scatter(Z_2d[:, 0], Z_2d[:, 1],
                          c=all_series_ids, cmap='tab10', alpha=0.5, s=10)
axes[0].set_title('Embeddings colored by Series ID')
axes[0].set_xlabel('PC1')
axes[0].set_ylabel('PC2')
plt.colorbar(scatter, ax=axes[0], label='Series ID')
axes[0].grid(True, alpha=0.3)

# 時刻で色分け（軌道位相）
orbital_phase = (all_positions % 90) / 90.0
scatter2 = axes[1].scatter(Z_2d[:, 0], Z_2d[:, 1],
                           c=orbital_phase, cmap='hsv', alpha=0.5, s=10)
axes[1].set_title('Embeddings colored by Orbital Phase')
axes[1].set_xlabel('PC1')
axes[1].set_ylabel('PC2')
plt.colorbar(scatter2, ax=axes[1], label='Orbital phase (0-1)')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('timesiam_embeddings.png', dpi=150, bbox_inches='tight')
plt.show()