Hankel Non-negative Matrix Factorizationによる時系列データの前処理

HNMF

先日ブログで時系列データをHankel Singular Value Decompositionという低周波成分と高周波成分に分解するという前処理のアルゴリズムを紹介しました． 今回は特異値分解（SVD）の低ランク近似の部分を非負値行列因子分解（NMF）に変えても分解できるのかを試しました． とりあえずHankel Non-negative Matrix Factorization(HNMF)と呼びます．

詳しいアルゴリズムは先日のブログを見てください． 差分はDecompositionのSVDをNMFに変えただけです．

例

NMFを使用するので，時系列データは非負値であるという制約があります． 以下のようなsin波形にノイズを付加したデータを使用します．

データ

raw

分解

低周波

low

高周波

high

まとめ

summary

HSVDと同様に低周波成分と高周波成分に分解できたことが確認できました．

実装

import numpy as np
from sklearn.decomposition import NMF

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 hnmf(x, window, rank, random_state = 0):
    assert(np.min(x) > 0.0)
    m = moving_window_matrix(x, window)
    model = NMF(n_components=rank, init='random', random_state=random_state)
    w = model.fit_transform(m)
    h = model.components_
    hankel = w @ h
    c = hankel[0,:]
    c = np.append(c, hankel[1:,-1])
    return c, x-c

githubにあげておきました．

github.com

Usage

import numpy as np
from tfilter import hnmf

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

window = 100
rank = 3
low_freq, high_freq = hnmf(x, window, rank)

参考

• suzuzusu.hatenablog.com