カルマンフィルタの導出

はじめに

wiki: 「カルマンフィルタ\((Kalman Filter, KF\)と略す\()\) は、誤差のある観測値を用いて、ある動的システムの状態を推定あるいは制御するための、無限インパルス応答フィルタの一種である」。観測値に観測雑音、状態予測値にシステム雑音があって場合により時間とともに誤差がドンドン蓄積してそのまま使えないので、観測と予測のガウス性を利用した線形センサデータフュージョンがカルマンフィルタの原点である。

状態空間モデル

時系列解析の中で、予測値と観測値の間何らかの因果関係を見つけて、何らかの方法でそれらのデータを絡んでモデル化して状態を推定していく。ここで汎用的な状態方程式観測方程式は以下の式にする。
$$\begin{eqnarray*}&& x_t = F_{t-1}(x_{t-1} ) + q_{t-1}\\&& y_t = H_t(x_t) + r_t\end{eqnarray*}$$
ただし、システム雑音、観測雑音\((q_t, r_t)\)は期待値が\(0\)、分散が\((Q_t, R_t)\)の独立正規分布(ガウス分布)にする。
$$\begin{eqnarray*}&& p(q_t) ~ N(0,Qt)\\&& p(r_t) ~ N(0,Rt)\\&& E[q_t r_t]=0 \end{eqnarray*}$$

基本原理

カルマンフィルタの原理は以下のイメージに示すように、最小解析誤差推定(最小二乗)から予測値、観測値をそれぞれの割合\(I-K_t H_t , K_t\)で線形合成した値を状態推定値にする。

kalman_filter_principle_n
kalman_filter_principle_n

式の導出

ベイズの定理やガウス分布からいくつかのの導出方法はあるが、本記事はシンプルな最小二乗法からカルマンフィルタの黄金\(5\)式\(\space *\)を誘導する。

線形システム

\(F(.),H(.)\)とも線形関数の場合、以下線形カルマンフィルタの誘導となる。例えば\( F(x)=Ax+B\)、\(A,B\)とも実数行列、\(x\)は変数行列とする。

まずは、予測値(数学期待)を求めるのは出発点なので、すべては以下の式から始まる。
$$ \gamma_t = F_{t-1}(\gamma_{t-1}) \space \space c $$
推定値\(x_t^{est}\)は、重み\(K_t\)の調整で予測値\(x_t\)と観測値\(y_t\)の間にあろうとして、以下線形表現がある。
$$ x_t^{est} = x_t + K_t(y_t-H_t x_t) \\
{\small = F_{t-1}(x_{t-1}^{true}+p_{t-1}) + q_{t-1}+ K_t\{(H_t x_t^{true}+r_t) }\\
{\small – H_t(F_{t-1}(x_t^{true}+p_{t-1})+q_{t-1})\} }\\
{\scriptsize = x_t^{true}+F_{t-1}p_{t-1}+q_{t-1}+K_t(r_t-H_t F_{t-1}p_{t-1}-H_t q_{t-1}) }$$

解析誤差(共分散\(P_t\))は以下の式より求める。
$$ P_t = E[(x_t^{est}-x_t^{true})^2] \space \space \\
{\scriptsize = E[(F_{t-1}p_{t-1} + q_{t-1} + K_tr_t-K_t H_t F_{t-1}p_{t-1}-K_t H_t q_{t-1})^2] }\\
{\scriptsize = F_{t-1}P_{t-1}F_{t-1}^T-2F_{t-1}P_t F_{t-1}^TH_t^T K_t^T+Q_{t-1}-2Q_{t-1}H_t^T K_t^T }\\
{\scriptsize + K_t R_t K_t^T+K_t H_t F_{t-1}P_{t-1}F_{t-1}^TH_t^T K_t^T+K_t H_t Q_{t-1}H_t^T K_t^T } $$

ここから二乗の最小化より解析誤差(共分散\(P_t\))を最小値に至らせる\(K_t\)を求める。
$$ \frac {\partial P_t } {\partial K_t } = 0 \\
{\small \frac {\partial P_t } {\partial K_t } = -2F_{t-1} P_{t-1} F_{t-1}^T H_t^T -2Q_{t-1}H_t^T } \\
{\small +2K_t R_t + 2K_t H_t F_{t-1} P_{t-1} F_{t-1}^T H_t^T +2K_tH_tQ_{t-1}H_t^T } \\
{\small K_t = (F_{t-1}P_{t-1}F_{t-1}^T + Q_{t-1}) H_t^T \{R_t } \\
{\small + H_t(F_{t-1}P_{t-1}F_{t-1}^T+Q_{t-1})H_t^T\}^{-1} } $$

これで、状態予測時\( K_t = 0 \)より、式\(\cdot\)を\({\small P_t = F_{t-1}P_{t-1}F_{t-1}^T + Q_{t-1} \space ** }\)に、\( K_t\)は、以下のように簡略化する。
$$ K_t = P_t H_t^T (R_t + H_t P_t H_t^T)^{-1} \space \space *** $$
これから観測値を入れて、\(K_t\) を使って推定値を求める。
$$ x_t^{est} = x_t + K_t(y_t – H_t x_t) \space \space **** $$
推定値を求めたら、続いて\(P_t\)を更新する。前述した\(#\)式の展開の項らより、最後の式は誘導できる。
$$ P^{\prime}_t = (I – K_t H_t)P_t \space \space ***** $$
\(P^{\prime}_t\)は、\((t+1)\)時刻の\(P_t\)として利用される。

非線形システム

線形カルマンフィルタの\(F(.),H(.)\)は非線形関数の場合、\(f(.),h(.)\)と記す。例えば\(f(x)=Ax^2+B\)、\(A,B\)とも実数行列、\(x\)は変数行列として、ティラー展開の最初の\(2\)項のみ使って、\( \hat{f}(x)=f(a)+\frac {\partial {f(x)}} {\partial {x}}|_{x=a}x \) のように\(Aa^2+B+2A(Aa^2+B)x\)に線形化することができる。\(3\)項目以降は捨てるので大いに誤差を招く可能性ある。状態予測の際に線形カルマンフィルタの\(F\)を\(f\)、観測値の取り入れる際に線形カルマンフィルタの\(H\)を\(h\)に、共分散行列の計算の際に\(F,H\)を\( \frac {\partial {f(x)}} {\partial {x}}|_{x=x_{t-1}}, \frac {\partial {h(x)}} {\partial {x}}|_{x=x_t}\)に置き換えて済む。

これまでカルマンフィルタの誘導で、状態値推定の手順を以下のとおり整理しておく。

推定の手順

時系列の状態推定に応用して、以下状態予測\(Predict\)→測定更新 \(Measurement \space Update\)→状態更新\(Correct, Update\)→時間更新\(Time \space Update\)→次の状態予測\(Predict\)→\(…\)の繰り返すことによって、時系列とともに状態の数学期待最小解析誤差(共分散)を求める手順となる。

Step0=パラメータ初期化(t=0)

状態推定(数学期待)、解析誤差共分散の初期化
$$\begin{eqnarray*}&& \hat{x}_0 = E[x_0]\\&& \hat{P_0} = E[(x_0-\hat{x_0})(x_0-\hat{x_0})^T]\end{eqnarray*}$$

Step1=状態値、共分散行列予測(t>0)

$$\begin{eqnarray*}&& \bar{x_t} = F_{t-1}(\hat{x}_{t-1})\\&& \bar{P_t} = F_{t-1}\hat{P}_{t-1}F_{t-1}^T + Q_{t-1}\end{eqnarray*}$$
ただし、\(\bar{x_t}\)は予測値、\(\bar{P_t}\)は事前共分散行列、\(Q_{t-1}=E[q_{t-1}q^T_{t-1}]\)

Step2=カルマンゲイン、状態値、共分散行列更新

$$\begin{eqnarray*}&& K_t = \bar{P_t}H_t^T(H_t\bar{P_t}H_t^T + R_t)^{-1}\\
&& \hat{x_t} = \bar{x_t} + K_t(y_t – H_t\bar{x_t})\\&&\hat{P_t} = (I-K_tH_t)\bar{P_t} \end{eqnarray*}$$
ただし、\(K_t\)はカルマンゲイン、\(\hat{x_t}\)は状態推定値、\(\hat{P_t}\)は事後共分散行列、\(R_{t}=E[r_{t}r^T_{t}]\)

KFフローチャート

Kalman_Filter_Flowchart_2
Kalman_Filter_Flowchart_2

拡張カルマンフィルタ

拡張カルマンフィルタ\((Extended Kalman Filter, EKF\)と略す\()\)は、非線形フィルタリングである。前述した状態方程式、観測方程式より、以下の状態空間モデルの\(f(⋅)\)または\(h(⋅)\)が非線形関数であり、拡張カルマンフィルタが適用される。テイラー展開より、2次微分以降の項目を省略して、非線形である\(f(⋅), h(⋅)\)の1次微分を線形化とし、前述したカルマンフィルタのアルゴリズムが適用可能となる。しかし、\(f(⋅), h(⋅)\)の微分では\(f(⋅), h(⋅)\)の一部しか表現できず、この線形化処理(1次微分)が誤差を大きく招く可能性がある。
$$\begin{eqnarray*}&& F_{t-1} =\frac{\partial f_{t-1}(x)}{\partial x}|_{x=\hat{x}_{t-1}}\\&& H_{t} =\frac{\partial h_{t}(x)}{\partial x}|_{x=\hat{x}_{t}}\end{eqnarray*}$$

状態空間モデル

$$\begin{eqnarray*}&& x_t = f_{t-1}(x_{t-1}) + q_{t-1}\\&& y_t = h_t(x_t) + r_t\end{eqnarray*}$$

Step0=パラメータ初期化(t=0)

$$\begin{eqnarray*}&& \hat{x}_0 = E[x_0]\\&& \hat{P_0} = E[(x_0-\hat{x_0})(x_0-\hat{x_0})^T]\end{eqnarray*}$$

Step1=状態値、共分散行列予測(t>0)

$$\begin{eqnarray*}&& \bar{x_t} = f_{t-1}(\hat{x}_{t-1})\\&& \bar{P_t} = F_{t-1}\hat{P}_{t-1}F_{t-1}^T + Q_{t-1}\end{eqnarray*}$$
ただし、\(\bar{x_t}\)は予測値、\(\bar{P_t}\)は事前共分散行列、\(Q_{t-1}=E[q_{t-1}q^T_{t-1}]\)

Step2=カルマンゲイン、状態値、共分散行列更新

$$\begin{eqnarray*}&& K_t = \bar{P_t}H_t^T(H_t\bar{P_t}H_t^T + R_t)^{-1}\\
&& \hat{x_t} = \bar{x_t} + K_t(y_t – h_t\bar{x_t})\\&&\hat{P_t} = (I-K_tH_t)\bar{P_t} \end{eqnarray*}$$
ただし、\(K_t\)はカルマンゲイン、\(\hat{x_t}\)は状態推定値、\(\hat{P_t}\)は事後共分散行列、\(R_t=E[r_t r^T_t]\)

EKFフローチャート

Extended_Kalman_Filter_Flowchart_2
Extended_Kalman_Filter_Flowchart_2

拡張カルマンフィルタを6軸IMUへの適用

6軸IMUへの適用例として、本サイトの記事 6軸IMU~拡張カルマンフィルタ に載ってある。

カルマンフィルタの再考

カルマンフィルタでは、状態推定値を予測結果\(x_t\)(実装例ではジャイロセンサーデータ)と観測データ\(y_t\)(実装例では加速度センサデータ)の線形結合で合成し,その誤差を最小にする推定法だと分かる。これはシステム誤差、観測誤差の数学期待が0の正規分布との前提条件から由来した推定法である。しかし、カルマンフィルタをかけることで、状態推定値は予測結果と観測データの間にあるのは、真値からかなり乖離してしまう場合にあるのか。これはシステム誤差と観測誤差が無相関かつ直交という前提から、勿論真値が推定値と観測値の間にある結論を結ぶ考えである。また、ジャイロセンサーデータと、加速度センサデータとも観測値\(y_t\)にする方法がある。それぞれのパーフォーマンスの検証は、比較的精確な実験環境(比較用の高精度ジャイロセンサ、加速度センサ、エンコーダ、モータ)がないと、実は容易ではない。というよりも、シミュレーションをかけてカルマンフィルタのアルゴリズムを検証するのが、確実に可能である。

参考文献

1-wikipedia: カルマンフィルタオイラー角
2-Greg Welch氏、Gary Bishop氏: An Introduction to the Kalman Filter
3-田島洋氏: マルチボディダイナミクスの基礎―3次元運動方程式の立て方

関連記事

オイラー角~ジンバルロック~クォータニオン
SLAM~拡張カルマンフィルタ
SLAM~Unscentedカルマンフィルタ
9軸IMUセンサ 6軸/9軸フュージョン 低遅延 USB出力 補正済み ROS対応
9軸IMU ICM-20948をロボットに組み込もう
YDLIDAR G4=16m 薄型 ROS対応SLAM LIDAR
研究開発用 台車型ロボット キット

4+

ロボット・ドローン部品お探しなら
ROBOT翔・電子部品ストア

Mask R-CNNを試す

はじめに

Mask R-CNNとはICCV 2017 Best Paper に選出された手法で、物体検出Object Dectectionやセマンティック・セグメンテーションSemantic Segmentationを実現するための手法である。COCOデータセットにより学習した、Matterport Mask_RCNNモデルを利用して、デモ画像より物体検出、セグメンテーションデモをGoogle colabで動かしてみよう。

デモを動かそう

Matterport Mask_RCNN、COCO API・Datasetのインストール、デフォルトのデモを動かす手順の以下の通り。

・Mask_RCNNのインストール、セットアップ

%cd /content/drive/My Drive
!git clone https://github.com/matterport/Mask_RCNN.git
%cd ./Mask_RCNN
!pip install -r requirements.txt
%run -i setup.py install

・COCO APIのインストール、セットアップ

%cd ..
!git clone https://github.com/waleedka/coco.git
%cd ./coco/PythonAPI
%run -i setup.py build_ext --inplace

・COCO Datasetで学習したMask_RCNNモデルのインストール

import os
import sys
import random
import math
import numpy as np
import skimage.io
import matplotlib
import matplotlib.pyplot as plt

# Root directory of the project
ROOT_DIR = os.path.abspath("/content/drive/My Drive/Mask_RCNN")

<pre class="brush: actionscript3; gutter: false">
# Import Mask RCNN
sys.path.append(ROOT_DIR)  # To find local version of the library
from mrcnn import utils
import mrcnn.model as modellib
from mrcnn import visualize
# Import COCO config
sys.path.append(os.path.join(ROOT_DIR, "samples/coco/"))  # To find local version
import coco

%matplotlib inline 

# Directory to save logs and trained model
MODEL_DIR = os.path.join(ROOT_DIR, "logs")

# Local path to trained weights file
COCO_MODEL_PATH = os.path.join(ROOT_DIR, "mask_rcnn_coco.h5")
# Download COCO trained weights from Releases if needed
if not os.path.exists(COCO_MODEL_PATH):
    utils.download_trained_weights(COCO_MODEL_PATH)

# Directory of images to run detection on
IMAGE_DIR = os.path.join(ROOT_DIR, "images")

class InferenceConfig(coco.CocoConfig):
    # Set batch size to 1 since we'll be running inference on
    # one image at a time. Batch size = GPU_COUNT * IMAGES_PER_GPU
    GPU_COUNT = 1
    IMAGES_PER_GPU = 1

config = InferenceConfig()
# config.display()

いよいよ認識しようと、以下のpythonコードを実行する。

# Create model object in inference mode.
model = modellib.MaskRCNN(mode="inference", model_dir=MODEL_DIR, config=config)

# Load weights trained on MS-COCO
model.load_weights(COCO_MODEL_PATH, by_name=True)

# COCO Class names
# Index of the class in the list is its ID. For example, to get ID of
# the teddy bear class, use: class_names.index('teddy bear')
class_names = ['BG', 'person', 'bicycle', 'car', 'motorcycle', 'airplane',
               'bus', 'train', 'truck', 'boat', 'traffic light',
               'fire hydrant', 'stop sign', 'parking meter', 'bench', 'bird',
               'cat', 'dog', 'horse', 'sheep', 'cow', 'elephant', 'bear',
               'zebra', 'giraffe', 'backpack', 'umbrella', 'handbag', 'tie',
               'suitcase', 'frisbee', 'skis', 'snowboard', 'sports ball',
               'kite', 'baseball bat', 'baseball glove', 'skateboard',
               'surfboard', 'tennis racket', 'bottle', 'wine glass', 'cup',
               'fork', 'knife', 'spoon', 'bowl', 'banana', 'apple',
               'sandwich', 'orange', 'broccoli', 'carrot', 'hot dog', 'pizza',
               'donut', 'cake', 'chair', 'couch', 'potted plant', 'bed',
               'dining table', 'toilet', 'tv', 'laptop', 'mouse', 'remote',
               'keyboard', 'cell phone', 'microwave', 'oven', 'toaster',
               'sink', 'refrigerator', 'book', 'clock', 'vase', 'scissors',
               'teddy bear', 'hair drier', 'toothbrush']
# Load a random image from the images folder
file_names = next(os.walk(IMAGE_DIR))[2]
for file_name in file_names:
  image = skimage.io.imread(os.path.join(IMAGE_DIR, file_name))

  # Run detection
  results = model.detect([image], verbose=1)

  # Visualize results
  r = results[0]
  visualize.display_instances(image, r['rois'], r['masks'], r['class_ids'], 
                              class_names, r['scores'])

うまくいけばSemantic Segmentationでマスクしたデモ画像が表示される。

Mask_RCNNデモ
Mask_RCNNデモ

※ Mask_RCNNの利用は、以下のとおりtensorflowバージョンを1.xに、

%tensorflow_version 1.x
import tensorflow
print(tensorflow.__version__)

ランタイムのタイプをGPUに設定して再起動して、またもう1回tensorflowバージョンを1.xにして確認する。

%tensorflow_version 1.x
import tensorflow
print(tensorflow.__version__)

デモの画像を入れ替えてみよう

images直下のデモ画像を入れ替えて上記pythonコードを実行してたら、以下画像のようにDining tableの一部が未検出であった。

Mask_RCNNテーブル一部検出
Mask_RCNNテーブル一部検出

Notebook ipynbファイルがGithubへ公開済み。

感想

デモ画像がよさそうに検出できたように見えますが、入れ替えたらそうでもない結果となった。やはり検出の正確性が学習モデルに大いに相関することで、専用学習データで学習モデルを作成しないと納得いく結果が得られず。画像データアノテーションImage Data Annotation業務が請負可能な業者さんがドンドン増えているらしい。

参考文献

Matterport Mask_RCNN on Github.

以上

1+

ロボット・ドローン部品お探しなら
ROBOT翔・電子部品ストア

struct2depth~単眼カメラ2D camera Visual SLAM

はじめに

Google がTensorflowのResearch Modelとしてstruct2depth、vid2depthを公開したので、struct2depthを利用して単眼カメラMonocular Cameraで撮った写真から深度Depthを推定してみよう。struct2depth、vid2depthは、KITTIまたは、CITYSCAPEの学習データを通してVisual Odometry、Depthの推定を習得するモデルである。また他の学習データを入れ替えてもあり得ると考えられる。SFM:Structure From Motionに基づく技術で、Depth深度まで推定できれば、3D Recontruction3次元復元まで使われる。

実測

雑居ビール内、ドラッグストア前および、ホールで写真を撮って完了とした。

推定

画像サイズを416×128に縮小して、推定の時間を短縮する。

環境

・ Google Colab, 18.04.3 LTS Bionic Beaver, GPU Tesla k80
・ Tensorflow 1.15.2
・ Research model struct2depth/KITTI

手順

学習せずKITTIモデルをそのまま利用したので、推定手順は以下のとおり。
・tensorflow_versionを1.xに合わせる。

・ランタイムを再起動。

%tensorflow_version 1.x
import tensorflow
print(tensorflow.__version__)

・以下確認できるまで、またランタイムを再起動する。

TensorFlow 1.x selected.
1.15.2

・インファレンス

!python inference.py --logtostderr --file_extension png --depth --egomotion true --input_dir image --output_dir output --model_ckpt model/KITTI/model-199160

結果

単眼カメラで撮ったRGB写真、レンダリングした深度推定イメージを結果として出力される。点群データの3Dイメージは別途プログラムを作成してレンダリングRenderingとする。

うまくいく例

完璧ではないが、扉、旗まで殆ど良く推定できている。

struct2depth_depth_ok_case
struct2depth_depth_ok_case

mayaviで点群Point Cloudデータの3D表現

Mayaviは、matplotlibよりパワーアップして、強力なエンジンVTKを利用した3Dツールである。

point_cloud_3d_plot
point_cloud_3d_plot

上図のように3Dで写真を細かく表現できた。点群データ(npyファイル)による3D表現のpythonソースは、Githubへ公開済み。

うまくいかない例

左下に推定が失敗と見られる。他の場所はなんとなく推定てきている。

struct2depth_depth_ng_case
struct2depth_depth_ng_case

原因を探る

・ KITTIモデルは屋外モデルでそのままでは屋内に向かない場合ある。測定環境にふさわしい学習データセット(モデル)が必要である。
・ 照明の強弱、特徴量に大きく関わること。
・ ついてはまだ実験が不十分だが、商用可能なVisual SLAMに道が長く感じさせられる。

参考文献

Depth Prediction Without the Sensors: Leveraging Structure for Unsupervised Learning from Monocular Videos, Auther: Vincent Casser etc
github google tensorflow model struct2depth

以上

1+

ロボット・ドローン部品お探しなら
ROBOT翔・電子部品ストア

屋内3D地図で可視化、空間情報をスマート管理

インドアマップが活躍の時代に

オフィスビル、デパート、工場など建物の内部は構造が複雑で、通常の2D地図では各場所の違いを十分に表現できない課題がある。地理情報システムとビッグデータ技術の進歩により、3D地図が室内外の空間情報の可視化技術で設備管理のスマート化やIoTのソリューションを提供することで、この課題を解決する日が近づいている。商業施設、産業IoTの他、交通監視、観光、旅行、セキュリティ、消防、会議、展示、娯楽、公共サービス等の分野で活躍する時代が訪れる。

屋内3D地図が活躍の時代に
屋内3D地図が活躍の時代に

プロダクト・サービス

お客様に対して地図可視化プラットフォームを提供して、クラウド上で各シーンに対応した情報システムを構築する。弊社とFengMap社と連携して、屋内外の3Dマップの作成サービスを提供する。また、開発者向けには専用のエンジンを提供してより簡単に各OS環境に対応したマップアプリの開発ができるようになる。

プロダクト・サービス
プロダクト・サービス

商業施設へ展開の例

商業施設では、CADデータによって室内のデータモデルを構築して各店舗の経営内容を組み込むことで空間データモデルを形成する。それにより、ショッピング案内、店舗管理、経営状況などの情報を共有し、可視化できる。スマート現場クラウド監視プラットフォームを提供する。

スマート現場クラウド監視プラットフォーム
スマート現場クラウド監視プラットフォーム

産業IoTへ展開の例

産業用として、Fengmapは可視化技術をIoTと融合し、設備の位置確認機能構内の設備、車両、人員の所在地と状態を把握して作業のモニタリング、消費エネルギー量の管理、データ統計などを行うことができる。スマート工場可視化管理システムを提案する。

スマート工場可視化管理システム
スマート工場可視化管理システム

フェングマップ社について

Beijing FengMap Technology Co.LTDは、2013年に設立された、北京に拠点を置く技術会社です。同社は、屋内および屋外の空間情報の可視化研究と開発に焦点を当て、地図データの作成、地図の編集、 ストレージ、マップ統合ソフトウェアアプリケーション開発。 空間情報の可視化技術に基づいて、資産管理、人事管理、施設および環境の監視、リモート制御、データ分析を含むさまざまな管理ソフトウェアシステムを多くの顧客に提供しました。創業以来、同社は商業用不動産、工業用IoT、工業団地から、家庭および幅広い公共サービスまで、多くの顧客を獲得した。500社、8000案件を開発した実績をもつという。
英文サイト→ https://www.fengmap.com/en/

日本総代理店

フェングマップ3Dマッピング作成サービスおよびSDK販売。
株式会社翔雲 令和2年3月1日から新住所↓
〒260-0026 千葉市中央区千葉みなと2-2-1502
代表取締役 柳建雄 電話 050-3598-8286
会社サイト https://soarcloud.com
技術情報サイト https://memo.soarcloud.com
販売サイト https://store.soarcloud.com

令和2年4月~5月特別キャンペーンお知らせ

上記時間限定、利益なし特別価格で3D地図作成サービスをご利用いただけます。どうぞお気軽にお問合せされるよう宜しくお願い申し上げます。

1+

ロボット・ドローン部品お探しなら
ROBOT翔・電子部品ストア

機械学習の13、SVD特異値分解

はじめに

本文には、大文字表現=行列/マトリクス、\(\boldsymbol{bold}\)小文字表現=ベクトル\(R^d\)、普通小文字表現=スカラー\(R\) と記す。

SVD(Singular Value Decomposition)は機械学習の分野で広く使用されているアルゴリズムで、次元削減アルゴリズムの特徴分解だけでなく、推薦システム(Recommender System)や自然言語処理(Nature Language Process)にも使用される。

原理

SVDは前述した特徴分解同じく行列を分解するが、SVDは分解する行列が正方行列にする必要ない。行列Aの形状が\(m×n\)であると仮定すると、行列AのSVDを次のように定義する。
$$ A = UΣV^T $$
ただし、\(U\)は\(m×m\)の行列\((u_1, u_2,…, u_m)\)、\(Σ\)は\(m×n\)の行列、主対角線上の要素を除く全ての要素が実数ゼロ\(0\)であり、主対角線上の各要素は特異値(Singular Value)という、\(V\)は\(n×n\)行列\((v_1, v_2,…, v_n)\)。

\(A\)の転置\(A^T\)と\(A\)を行列で乗算すると、\(n×n\)の正方行列\(A^T A\)が得られる。\(A^T A\)は正方行列であるため、特徴の分解を実行でき、得られた固有値と固有ベクトルは次の式を満たす。
$$ (A^T A)\boldsymbol{v_i}=λ_i\boldsymbol{v_i} $$
これで行列\(A^T A\)の\(n\)個の固有値と対応する\(n\)個の固有ベクトル\(v_i\)を取得できる。

\(A\)と\(A\)の転置\(A^T\)を行列で乗算すると、\(m×m\)の正方行列\(AA^T\)が得られる。\(AA^T\)は正方行列であるため、特徴の分解を実行でき、得られた固有値と固有ベクトルは次の式を満たす。
$$ (AA^T) \boldsymbol{u_i}=λ_i\boldsymbol{u_i} $$
これで行列\(AA^T\)の\(m\)個の固有値と対応する\(m\)個の固有ベクトル\(u_i\)を取得できる。

\(Σ\)は対角線上の特異値を除いて全て\(0\)で、各特異値\(σ_i\)を見つけるだけで\(Σ\)が求められる。
$$ \sigma_i = \sqrt{λ_i} $$

各特異値\(σ_i\)のうち、比較的大きいほう(主成分)とそれに対応する特異ベクトル\(u_i, v_i\)を\(k\)個\((k << n)\)残すとAの次元を削減する。
$$ A_{m×n}=U_{m×m}Σ_{m×n}V^T_{n×n} ≈ U_{m×k}Σ_{k×k}V^T_{k×n} $$

実装

以下行列dataSetに対して、SVDアルゴリズムで\((U, \Sigma, V^T)\)を求めて、5次元→3次元つまり2次元を削減してが新しい\(\Sigma\)で\((U* \Sigma*V^T)\)が元の行列dataSetに戻せるかを確かめる。

from numpy import *
def loadExData():
    return[[0, 0, 0, 2, 2],
           [0, 0, 0, 3, 3],
           [0, 0, 0, 1, 1],
           [1, 1, 1, 0, 0],
           [2, 2, 2, 0, 0],
           [5, 5, 5, 0, 0],
           [1, 1, 1, 0, 0]]
dataSet = loadExData()
U, Sigma, VT = linalg.svd(dataSet)
print(f'dataSet:\n{dataSet}')
print(f'U:\n{U}\nSigma:\n{Sigma}\nVT:\n{VT}')
// 小さいSigmaを0にする(削除)
Sig3 = mat([[Sigma[0], 0, 0], [0, Sigma[1], 0], [0, 0, Sigma[2]]])
print(f'U[:,:3] * Sig3 * VT[:3,:]:\n{U[:,:3] * Sig3 * VT[:3,:]}')

ソースコード
https://github.com/soarbear/Machine_Learning/tree/master/svd

結果

\((U* \Sigma*V^T)\)は元のdataSetとほぼ同じ行列だと分かる。

dataSet:
[[0, 0, 0, 2, 2], [0, 0, 0, 3, 3], [0, 0, 0, 1, 1], [1, 1, 1, 0, 0], [2, 2, 2, 0, 0], [5, 5, 5, 0, 0], [1, 1, 1, 0, 0]]
U:
[[-2.22044605e-16  5.34522484e-01  8.41641151e-01 -1.37443101e-02
  -7.57428665e-02 -1.11022302e-16  1.38777878e-17]
 [ 0.00000000e+00  8.01783726e-01 -4.92426901e-01 -2.47257115e-01
   2.31349353e-01  3.15719673e-16 -2.77555756e-17]
 [ 0.00000000e+00  2.67261242e-01 -2.06001597e-01  7.69259966e-01
  -5.42562325e-01 -7.55450741e-16  1.09551769e-16]
 [-1.79605302e-01  2.77555756e-17 -3.00520660e-02 -2.15935735e-01
  -2.94749442e-01  9.05439185e-01 -1.16246358e-01]
 [-3.59210604e-01  5.55111512e-17 -6.01041319e-02 -4.31871470e-01
  -5.89498885e-01 -4.19124526e-01 -3.97074256e-01]
 [-8.98026510e-01  0.00000000e+00  3.60624791e-02  2.59122882e-01
   3.53699331e-01  5.40010673e-16 -6.71525577e-17]
 [-1.79605302e-01  2.77555756e-17 -3.00520660e-02 -2.15935735e-01
  -2.94749442e-01 -6.71901321e-02  9.10394870e-01]]
Sigma:
[9.64365076e+00 5.29150262e+00 8.36478329e-16 6.91811207e-17
 1.11917251e-33]
VT:
[[-5.77350269e-01 -5.77350269e-01 -5.77350269e-01  0.00000000e+00
   0.00000000e+00]
 [-2.46566547e-16  1.23283273e-16  1.23283273e-16  7.07106781e-01
   7.07106781e-01]
 [-7.01908483e-01 -1.02844064e-02  7.12192890e-01 -2.22044605e-16
  -1.66533454e-16]
 [-4.17122461e-01  8.16431808e-01 -3.99309347e-01  0.00000000e+00
  -1.11022302e-16]
 [-0.00000000e+00 -1.96261557e-16  1.96261557e-16  7.07106781e-01
  -7.07106781e-01]]
U[:,:3] * Sig3 * VT[:3,:]:
[[ 4.47427211e-17  1.57774942e-15  2.08638397e-15  2.00000000e+00
   2.00000000e+00]
 [-7.56974048e-16  5.27282824e-16  2.29691224e-16  3.00000000e+00
   3.00000000e+00]
 [-2.27747782e-16  1.76121044e-16  5.16267387e-17  1.00000000e+00
   1.00000000e+00]
 [ 1.00000000e+00  1.00000000e+00  1.00000000e+00  1.03851855e-16
   1.03851855e-16]
 [ 2.00000000e+00  2.00000000e+00  2.00000000e+00  2.07703709e-16
   2.07703709e-16]
 [ 5.00000000e+00  5.00000000e+00  5.00000000e+00 -6.69808260e-33
  -5.02356195e-33]
 [ 1.00000000e+00  1.00000000e+00  1.00000000e+00  1.03851855e-16
   1.03851855e-16]]

参考文献

[1] PeterHarrington. Machine Learning in Action.

2+

ロボット・ドローン部品お探しなら
ROBOT翔・電子部品ストア