2019-11-08

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

HNMF

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

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

例

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

データ

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

分解

低周波

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

高周波

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

まとめ

f:id:suzuzusu:20191108033830p:plain — 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

2019-11-07

Hankel Singular Value Decompositionによる時系列データの前処理

機械学習系技術系

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)

参考

2019-11-06

しょぼん(´･ω･｀)基底をGPLVMで獲得する

機械学習系技術系

データの作成

以下の画像から，しょぼん(´･ω･｀)基底を作成します．

f:id:suzuzusu:20191106035235p:plain — オリジナル画像

まず，画像を読み込みます．

import numpy as np
from PIL import Image
image = np.asarray(Image.open('./data/shobon.png').convert('L'))

次に，2値化を行いしょぼん(´･ω･｀)基底を作成します．

# binarization
a = np.where(image < 240)
shobon = np.c_[a[1], np.flipud(a[0])]
# shobon.shape = (896, 2)

プロットすると以下のような，しょぼん(´･ω･｀)基底が得られたことが確認できます．

f:id:suzuzusu:20191106035727p:plain — しょぼん(´･ω･｀)基底

高次元変換

ガウス分布に従う重みを用いて，100次元の高次元にしょぼん(´･ω･｀)基底を変換します．そしてノイズを加えておきます．

dim = 100
w = np.random.normal(size=2*dim).reshape(2, dim)
high = shobon @ w
high = high + np.random.normal(size=high.shape[0]*high.shape[1]).reshape(high.shape[0], high.shape[1])
# high.shape = (896, 100)

GPLVMによる低次元埋め込み

Gaussian Process Latent Variable Models(GPLVM)とはガウス過程を利用して潜在変数空間に変換する教師なし学習のアルゴリズムです．応用例として以下の動画のように人間のポーズを低次元埋め込みするなどの試みがあります．

今回はGPLVMを使うためにGPyを使用します．

GPLVMを使って100次元のデータを2次元の低次元埋め込みをします．

# latent
gplvm = GPLVM(high, 2)
latent = gplvm.X
# latent.shape = (896, 2)

結果

得られた潜在変数をプロットすると以下のようになります．

f:id:suzuzusu:20191106040414p:plain — 潜在変数

しょぼん(´･ω･｀)基底と似たものが獲得できたので成功ですね！！

実装

github.com

参考

2019-11-04

Higher-Order-SVDで画像の低ランク近似

機械学習系技術系

HOSVD

SVDは行列の低ランク近似手法であり，Higher-Order-SVD(HOSVD)はSVDをテンソルに拡張した手法です． HOSVDではタッカー分解というコアテンソルとそれぞれのモードの行列に分解する方法を使用します．以下にその概要を図示しました．

f:id:suzuzusu:20191104153601p:plain — タッカー分解

方法

$m \times n \times l$ の3階テンソルを例にしてタッカー分解をする方法を説明します．まず最初に， $m \times n \cdot l$ ， $n \times l \cdot m$ ， $l \times m \cdot n$ のそれぞれの行列に分解しSVDをします．

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

上記の式から，分解されたユニタリ行列 $U$ からコアテンソル $S$ を以下のように定義します．

$S :=\{ U_{m}^{T}, U_{n}^{T}, U_{l}^{T} \} \cdot M$

元のテンソルに復元するためには，以下のように再度ユニタリ行列 $U$ を使います．

$M :=\{ U_{m}, U_{n}, U_{l} \} \cdot S$

上記のユニタリ行列 $U$ を任意のランクに切り取ることで，低ランク近似が可能となります．

例

3階テンソルであるカラー画像を使って，画像の低ランク近似をしてみます．使用する画像は画像処理でよく使われるレナ画像を使用します．

オリジナル画像

f:id:suzuzusu:20191104160723p:plain

Rank1の画像

f:id:suzuzusu:20191104160855p:plain

Rank5の画像

f:id:suzuzusu:20191104160906p:plain

Rank10の画像

f:id:suzuzusu:20191104160914p:plain

Rank20の画像

f:id:suzuzusu:20191104160924p:plain

Rank30の画像

f:id:suzuzusu:20191104160934p:plain

Rank50の画像

f:id:suzuzusu:20191104160942p:plain

Rank100の画像

f:id:suzuzusu:20191104160950p:plain

実装コード

gist.github.com

github.com

参考

2019-10-08

PyTorchで最急降下法による最適化

機械学習系技術系

余談

Differentiable Programming(微分可能プログラミング)という捉え方がとても素晴らしいと個人的に思います．詳しくは以下のブログを見てください．

bonotake.hatenablog.com

最急降下法による最適化

機械学習フレームワークのPyTorchを使って単純に最急降下法で関数の最適化をする．最適化する関数は10000000次元のsphere関数を使用する．

Sphere Function
出典：Sphere Function

コード

# -*- coding: utf-8 -*-

import torch
import torch.nn as nn
import torch.optim as optim

# 評価関数
def _sphere(x):
    return torch.pow(x, 2).sum()

class Model(nn.Module):
    def __init__(self, dim=100, func=None):
        super(Model, self).__init__()
        if func is None:
            func = _sphere
        self.func = func
        # 設計変数
        self.x = nn.Parameter(torch.rand(dim))

    def forward(self):
        return self.func(self.x)
    
    def vars(self):
        return self.x.detach().numpy()

    def objective(self):
        with torch.no_grad():
            return self.func(self.x).numpy()

# 次元は10000000次元
model = Model(dim=10000000)
optimizer = optim.SGD(model.parameters(), lr=0.1)

N = 100

print('初期の評価値', model.objective())

print('optimization...')
for i in range(N):
    output = model()
    optimizer.zero_grad()
    loss = output
    loss.backward()
    optimizer.step()
    print('\rloss:', loss.item(), end='')
print()

print('評価値', model.objective())

実行結果

初期の評価値 3331740.8
optimization...
loss: 2.1601908358880734e-13
評価値 1.382506e-13

gist

gist.github.com

参考

2019-10-07

GPyTorchで多クラス分類

機械学習系技術系

GPyTorchとは

PyTorch上で実装されたガウス過程(GP)のライブラリである．ハイパーパラメータをPyTorchのシステムを使用してGPUなどで最適化できるのが特徴である． GPは逆行列を計算する際に $\mathcal{O}(N^ 3)$ の計算量がかかり，スケーラビリティに難があるが，PyTorchのシステムを使ってそれを補おうとしている．

GPでは尤度を変更することで回帰や二値分類，多クラス分類が可能となる．今回はアヤメ(iris)データセットを用いて3クラスの多クラス分類をGPyTorchで分類してみる．

コード

# -*- coding: utf-8 -*-

import torch
import gpytorch
import numpy as np

from gpytorch.models import AbstractVariationalGP
from gpytorch.variational import CholeskyVariationalDistribution
from gpytorch.mlls.variational_elbo import VariationalELBO

from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split


class GPMultiClassificationModel(AbstractVariationalGP):
    def __init__(self, num_dim=2, grid_bounds=(-10., 10.), grid_size=64):
        variational_distribution = gpytorch.variational.CholeskyVariationalDistribution(
            num_inducing_points=grid_size, batch_size=num_dim
        )
        variational_strategy = gpytorch.variational.AdditiveGridInterpolationVariationalStrategy(
            self, grid_size=grid_size, grid_bounds=[grid_bounds], num_dim=num_dim,
            variational_distribution=variational_distribution, mixing_params=False, sum_output=False
        )
        super().__init__(variational_strategy)
        self.mean_module = gpytorch.means.ConstantMean()
        self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())

    def forward(self, x):
        mean_x = self.mean_module(x)
        covar_x = self.covar_module(x)
        latent_pred = gpytorch.distributions.MultivariateNormal(mean_x, covar_x)
        return latent_pred

device = 'cuda' if torch.cuda.is_available() else 'cpu'

iris = load_iris()
# 特徴量は2次元
# testデータは3割
(train_x, test_x, train_y, test_y) = train_test_split(iris.data[:, [0, 2]].astype(np.float32), iris.target, test_size=0.3, random_state=0,)
# 3クラス分類
n_labels = len(np.unique(train_y))

train_x = torch.from_numpy(train_x).to(device)
train_y = torch.from_numpy(train_y).to(device)
test_x = torch.from_numpy(test_x).to(device)
test_y = torch.from_numpy(test_y).to(device)

model = GPMultiClassificationModel(num_dim=2).to(device)

# 多クラス分類なのでSoftmaxを使用する
likelihood = gpytorch.likelihoods.SoftmaxLikelihood(num_classes=n_labels, num_features=2).to(device)

# 周辺尤度
mll = VariationalELBO(likelihood, model, train_y.numel())

optimizer = torch.optim.Adam(model.parameters(), lr=0.1)

# train
model.train()
likelihood.train()
training_iter = 10000
for i in range(training_iter):
    optimizer.zero_grad()
    output = model(train_x)
    loss = -mll(output, train_y)
    loss.backward()
    print('\rIter %d/%d - Loss: %.3f' % (i + 1, training_iter, loss.item()), end='')
    optimizer.step()
print('')

# test
model.eval()
likelihood.eval()
correct = 0.0
# sample数を64とする
with torch.no_grad(), gpytorch.settings.num_likelihood_samples(64):
    data = test_x
    target = test_y 
    output = likelihood(model(data))
    pred = output.probs.mean(0).argmax(-1)
    correct += pred.eq(target.view_as(pred)).cpu().sum()
    acc = correct.numpy() / float(len(test_x)) * 100.0
    print('Accuracy:', acc)

出力結果

Iter 10000/10000 - Loss: 0.397
Accuracy: 86.66666666666667

gist

gist.github.com

参考

2019-10-03

カーネル密度推定とエントロピー

機械学習系技術系

カーネル密度推定をしてエントロピーを計算する方法を忘備録として書いておく．

標準正規分布を例にとって計算してみる．ちなみに理論値は $\ln ({\sqrt{2 \pi e}})$

import numpy as np
import matplotlib.pyplot as plt

from scipy import integrate
from scipy.stats import norm
from scipy.stats import gaussian_kde

x = norm.rvs(size=10000)
plt.hist(x, bins=50)
plt.title('sampling')
plt.show()

k = gaussian_kde(x)
X = np.linspace(-5, 5, 100)
plt.plot(X, norm().pdf(X), label='norm')
plt.plot(X, k(X), label='kernel estimation')
plt.legend()
plt.title('pdf')
plt.show()

def entropy(pdf):
    def f(x):
        px = pdf(x)
        return - px * np.log(px)
    return f

entropy_k = entropy(k)

print('entropy:')
print('カーネル密度推定', integrate.quad(entropy_k, -5.0, 5.0))
print('理論値', np.log(np.sqrt(2*np.pi*np.e)))

出力結果

entropy:
カーネル密度推定 (1.4426170314750473, 1.2370907729071283e-08)
理論値 1.4189385332046727

f:id:suzuzusu:20191003215447p:plain

f:id:suzuzusu:20191003215500p:plain

gist

gist.github.com