HSVD

Hankel Singular Value Decomposition(HSVD)は時系列データを低周波と高周波の2つの成分に分解する手法です．HSVDを使用して時系列データの前処理をします．

Embedding

まず最初に， $N$ 個の時系列データ $X(t)$ を以下のような窓サイズ $L$ のハンケル行列H（正確にはハンケル行列は正方行列なのでちょっと違います...）に変換します．

$K = N - L + 1$

f:id:suzuzusu:20191107220308p:plain

Decomposition

得られたハンケル行列 $H$ を特異値分解して低ランク近似をして低周波成分を抽出します．

$H = U \Sigma V^{\ast}$

低ランク近似したハンケル行列を $H'$ とします．

f:id:suzuzusu:20191107223125p:plain

Unembedding

低周波成分 $C_L$ を以下のように抽出します．

$C_L = [ H'(1,1), H'(1, 2), \cdots, H'(1, K), H'(2, K), H'(3, K), \cdots, H'(L, K)$ ]

以下の画像のように赤い要素が低周波成分になります．

f:id:suzuzusu:20191107222527p:plain

高周波成分 $C_H$ は元のデータから低周波成分を引いたものになります．

$C_H(t) = X(t) - C_L(t)$

例

例としてsin波形にノイズを加えたデータをHSVDで前処理を施します．

データ

f:id:suzuzusu:20191107222404p:plain — raw

分解

低周波

f:id:suzuzusu:20191107222419p:plain — low

高周波

f:id:suzuzusu:20191107222435p:plain — high

まとめ

f:id:suzuzusu:20191107222453p:plain — summary

上の画像から元のデータを低周波と高周波に分解できていることが分かると思います．元のデータ = 低周波 + 高周波なのでそれぞれの成分を時系列予測アルゴリズムにかけるなどの応用が可能です．

実装

def moving_window_matrix(x,window_size):
    # Fork from https://qiita.com/bauer/items/48ef4a57ff77b45244b6
    n = x.shape[0]
    stride = x.strides[0]
    return np.lib.stride_tricks.as_strided(x, shape=(n-window_size+1, window_size), strides=(stride,stride) ).copy()

def hsvd(x, window, rank):
    m = moving_window_matrix(x, window)
    u, s, vh = np.linalg.svd(m)
    h = u[:,:rank] @ np.diag(s[:rank]) @ vh[:rank,:]
    c = h[0,:]
    c = np.append(c, h[1:,-1])
    return c, x-c

実装をgithubにあげました

github.com

Usage

import numpy as np
from tfilter import hsvd

N = 500
x = np.sin(np.arange(N) * np.pi/50.0)
x = x + np.random.normal(0, 0.3, size=N)

window = 100
rank = 2
low_freq, high_freq = hsvd(x, window, rank)