開かない瓶の蓋をバケットクラッシャーで破砕する

機械学習・バイオインフォマティクス関連の気になった内容を記事にしていきます．

特異値分解(singular value decomposition)でベクトルを直交化する

はじめに

特異値分解を使った直交化をする機会があったので書き留めておきます．数式で説明するか文章で説明するかはフィーリングです．

モチベーション

行列Mが与えられたとします．Mの列ベクトルを，変換後の空間のスケール(厳密な用語が分かりません...)を変えないように直交したベクトルに変換することを行いたいという話です．なぜかについては説明しません(面倒ですので・・・)．

説明

間違ってたらご一報頂けると幸いです．

特異値分解

行列分解の一種．行列MをM=U \Sigma Vのように分解します．ここで， U, Vは直交行列， \Sigmaは対角行列です．

直交行列の性質

行列が直交行列である，とは以下のことと同値です．

  1. 各行・列ベクトルが直交しておりノルムが1
  2. 逆行列は転置した行列
  3. 変換が等長写像(ノルムを変えない写像)になる

今回の直交化においてこの性質が大事になってきます．

特異値分解(SVD)による直交化

SVDによる直交化は M=U \Sigma V を利用して M'=MV ^T=U \Sigma のように直交化します．ここで，性質2より，Vの逆行列をかける代わりに転置行列をかけています．

特異値分解(SVD)による直交化の性質

SVDによる直交化では，被変換行列Mによる変換後の単位円を変えずに各ベクトルを直交にすることができます．以下，説明です．全てのノルムが1のベクトルを考えたとき，座標とみなしたベクトルの集合によって円ができます．これを単位円と言います．Mによってその各座標を変換すると単位円の形が変わります(楕円もしくは円に)．SVDによる直交化では，この楕円の形を変えないベクトルに直交化できるのです．もちろん，各ベクトルについて見てみれば直交化する前と後では値が変わってしまいます．ただ，ノルムが一定のベクトルによってできる，集合レベルでは変わらないようになっているということですね． この話の直感的な図がWikipediaにあります．

f:id:oaietnbvgq30:20191209170006p:plain
赤と黄色のベクトルは変換対象のベクトルです．Mによる変換(左上⇨右上)とU \Sigmaによる変換(左下⇨右上)では単位円の形を変えないことが分かります．
en.wikipedia.org

この直交化を行列Mに対して行うと，

MV ^T=U \Sigma

となり，SVDによって得られるUと\Sigmaで変換後の行列を得られることになります(もちろん元の行列にV ^Tをかけてもよい)． 更にU \Sigmaを見てみると直交行列と対角行列の行列積になっているので，性質1より，変換後の行列は列ベクトルが直交していることが分かります(行ベクトルは直交するとは限らない)．

要するに，V ^Tによる変換をしてあげると，変換後の空間のスケール(厳密な用語がry)を変えずに列ベクトルを直交化させることができる，というわけですね．

実験

実際にSVDを使って得られるベクトルを眺めてみます．

実験条件

pythonで動かしていきます． まずimportします．numpyと可視化関連です(個人的な好みでseabornを使ってますがimportしなくてもこの後のコードは実行できます，figureの見た目が良くなるだけです)．

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()

適当な行列を作ります．

M = np.array([[1,4],[3,4]])

SVDしてU \Sigmaを計算することでMを直交化します(numpyの関数だと \Sigmaは対角成分のみのベクトルで返ってくるためnp.diagで行列にしている)．

U, sigma, V = np.linalg.svd(M, full_matrices=True)
M_prime = np.dot(U, np.diag(sigma))

一応直交しているか確認．

np.dot(M_prime.T, M_prime)
# array([[  4.04164878e+01,  -1.21399145e-17],
# [ -1.21399145e-17,   1.58351216e+00]])

大丈夫そう．最後に可視化します．

fig, ax = plt.subplots(1, 1, figsize=(5, 5))
ax.set_xlim([-6,5])
ax.set_ylim([-6,5])
ax.plot([0,M.transpose()[0][0]], [0,M.transpose()[0][1]], color='b')
ax.plot([0,M.transpose()[1][0]], [0,M.transpose()[1][1]], color='b')
ax.plot([0,M_prime.transpose()[0][0]], [0,M_prime.transpose()[0][1]], color='r')
ax.plot([0,M_prime.transpose()[1][0]], [0,M_prime.transpose()[1][1]], color='r')
plt.show()

f:id:oaietnbvgq30:20191205012324p:plain

直交化後のM'の列ベクトルが描く単位円が変わってなさそうなことが分かります(適当)．