6月 24日講義資料の2 - 知能システム学講座

知能システム論1　(10_2)
2015.6.24
情報システム学研究科
情報メディアシステム学専攻
知能システム学講座
末廣尚士
- （擬似）逆行列と特異値分解

行列の特異値分解
A :(m×n)
は、
A=U W V T
U :(m×r)
V :(n×r )
の形で表現することが出来る
U, Ｖは列直交行列, ｒは行列のrank
W=diag (σ 1 , σ 2 ,... , σ r ) (σ 1≥σ 2≥...≥σ r≥0)
列直交行列
U=(e 1, e 2, , , , e r )
2015年度前学期
1i= j
0i≠ j
つまり
U T U=I
e i⋅e j =
I
は r×r
電気通信大学大学院　情報システム学研究科情報メディアシステム学専攻　知能システム論１
2
- 逆行列と特異値分解
m = n=r　の場合（Aが正則行列の場合）: 逆行列
−1
−1
A =V W U
W−1 =diag (1/ σ 1 ,1/ σ 2 ,... ,1/ σ n )
T
A−1 A=V W−1 U T U W VT =V W W−1 V T =V V T=I
そうでない場合：　擬似逆行列
+
−1
A =V W U
T
Aが正則
W−1 =diag (1/ σ 1 ,1/ σ 2 ,... ,1/ σ r )
A A + A=U W VT V W−1 U T U W VT =U W V T =A
2015年度前学期
電気通信大学大学院　情報システム学研究科情報メディアシステム学専攻　知能システム論１
3
- 他の逆行列を求める方法との比較
よくある方法
m > n　の場合
A +=(AT A)−1 A T
A A + A=A (A T A)−1 AT A=A
問題点：　T
T
−1
(A A)
A A
・　が正則とは限らない.つまり　・　正則であったとしても、危ない場合がある
が求められない．
特異値分解を行うとこれらの問題点が明確になるし、対応を工夫することが出来る
ただし、計算コストはかかる。
(疑似)逆行列の計算には特異値分解を使おう！
2015年度前学期
電気通信大学大学院　情報システム学研究科情報メディアシステム学専攻　知能システム論１
4
- 特異値分解でみるヤコビアンの性質
T
は、 J=U W V
J :(6×6)
U :(6×6) (直交行列)
V T :(6×6) (直交行列)
の形で表現することが出来る
W=diag (σ 1 , σ 2 ,... , σ 6 ) (σ 1≥σ 2≥...≥σ 6≥0)
Δ p=J Δ q=U W V T Δ q
Δ ̃p =W Δ q̃
ここで
=>
U T Δ p=W VT Δ q
Δ ̃p =U T Δ p , Δ q̃ =VT Δ q
 p と　の座標軸(直交基底)をうまく選ぶと、軸方向に
q
σ i 倍するだけの変換になる。
逆行列は軸方向に　1/ σ i 倍すればよいことになる。
2015年度前学期
電気通信大学大学院　情報システム学研究科情報メディアシステム学専攻　知能システム論１
5
- ヤコビアンの逆行列と特異点
Δ p=J Δ q=U W V T Δ q
J−1 =V W−1 U T
Δ q=V W−1 U T Δ p
W−1 =diag (1/ σ 1 ,1/ σ 2 ,... ,1/ σ 6 )
つまり
Δ q̃ =W−1 Δ ̃p
σ i が大きい：　動かしやすい方向
σ i が小さい：　動かしにくい方向
σ i がゼロ：　動かせない方向がある
σ i の最大と最小の比：特性値
特異点ということ
一般に擬似逆行列では、 σ i =0　のところは計算しない。すなわち対応する　の部分を0としたことと等しい。
1/
σi
さらに、アームの逆運動学解を求める場合などでは、　の場合も0に
σ i≈0
することで、解の発散を防ぐ。
2015年度前学期
電気通信大学大学院　情報システム学研究科情報メディアシステム学専攻　知能システム論１
6
- 特異点
特異点
p
p= f (q)
解が複数ある
q
山や谷（特異点）
を越えられない
どちらへもいけな
い、またはどちらへ
行くかわからない
pを少し動かそう
とするだけでqが
大きく変化する
2015年度前学期
qをうごかしてもpが
ほとんど変化しない
電気通信大学大学院　情報システム学研究科情報メディアシステム学専攻　知能システム論１
7
- 注意点

収束性


停留


local minimum に落ち込む
特異点


下手をすると大きく変化しすぎて違うところに飛び出
してしまう。
行列が正則でなくなる
複数解の存在

複数解のどれになるかの保障がない
2015年度前学期
電気通信大学大学院　情報システム学研究科情報メディアシステム学専攻　知能システム論１
8
- 特異値分解と最小二乗法



連立一次方程式の解
疑似逆行列の最小二乗，ノルム最小性
ニュートン法と（局所）最小二乗解
（ガウス・ニュートン法）
2015年度前学期
電気通信大学大学院　情報システム学研究科情報メディアシステム学専攻　知能システム論１
9
- 連立一次方程式の解(1)
b=A x
b の次元を m, x の次元を n ,
Aはm x n 行列，r=rank(A)とする．
(1) Aが正方でフルランク(m=n=r)の場合：
x=A−1 b
2015年度前学期
電気通信大学大学院　情報システム学研究科情報メディアシステム学専攻　知能システム論１
10
- 連立一次方程式の解(2)
(2) m>n=r の場合：
x=A + b
∥A x−b ∥
このとき，x は，　を最小とする解(最小二乗解)となる．
これは多数のデータをモデルにフィッティングする場合などに相当する．
これを特異値分解で考えてみる
A +=V W−1 UT
A x−b=U W VT V W−1 U T b−b=U b̃ −b
̃ =UT b=W V T x
ここで b
(m>nなので　U U T=I とならない．)
x
これが意味するところは，　は　で変えることのできる　の部分空間への
b̃
b
写像ということ．
x
の残りの部分はその直交補空間であり　で変えることが出来ない．　b
ということで最小二乗(近似)解になっていることが分かる．　2015年度前学期
電気通信大学大学院　情報システム学研究科情報メディアシステム学専攻　知能システム論１
11
- 連立一次方程式の解(3)
(3) m=r<n)の場合：
x=A + b
A x=b で　を最小とする解(ノルム最小解)となる
このとき，x は，　∥x∥
．
これは冗長自由度をもつアームの制御の場合などに相当する．
これを特異値分解で考えてみる
+
−1
A =V W U
T
(rank(A)=mなので　UU
T
=I となる．)
A x=U W V T V W−1 U T b=U U T b=b
x̃ =V T x=W−1 U T b
x̃
x
b
これが意味するところは，　は　に影響する　の部分空間への写像とい
ここで
うこと．
b
の残りの部分はその直交補空間であり　に影響しない．
x
その部分を　としたことでノルム最小解になっていることが分かる．　0
2015年度前学期
電気通信大学大学院　情報システム学研究科情報メディアシステム学専攻　知能システム論１
12
- 連立一次方程式の解(4)
(4) m>r, n>r の場合：
x=A + b
∥x∥
このとき，x は，　および　を最小とする解(最小二乗，
∥A x−b ∥
ノルム最小解)となる．
これは計測データでは決めることのできないパラメタを含むモデルを多数
のデータでフィッティングするような場合に相当する．
これも同様に特異値分解で考えることができる．
2015年度前学期
電気通信大学大学院　情報システム学研究科情報メディアシステム学専攻　知能システム論１
13
- 非線形最小二乗法



非線形最小化のニュートン法
二次微分を用いて極値を求める
非線形最小二乗（ガウス・ニュートン法）
誤差関数の二乗ノルムを最小化
誤差関数に対する一次のニュートン法
最小二乗性は疑似逆行列の最小二乗性による
Levenberg-Marquardt（レベンバーグ・マー
カート）法
ガウス・ニュートン法の収束性の改良版
極値から遠いときは最急降下法
2015年度前学期
電気通信大学大学院　情報システム学研究科情報メディアシステム学専攻　知能システム論１
14
- ガウス・ニュートン法

非線形最小二乗（ガウス・ニュートン法）のモデ
ルフィッティングへの応用
x
q
がパラメタ　を持つモデル
に当てはまることを以下で表現する．
f ( x , q)=0
xk
このときモデル　に当てはまる複数の計測値　から，モデルパラメタ
f
q
を推定することをここではモデルフィッティングと呼ぶことにする．
すなわち，　2
argmin ∑ ∥ f k ( x k , q)∥
q
2015年度前学期
電気通信大学大学院　情報システム学研究科情報メディアシステム学専攻　知能システム論１
15
- ガウス・ニュートン法
2
argmin ∑ ∥ f k ( x k , q)∥
q
fk
xk
これを　に関して分解し，計測値　が固定であることを考慮して簡
略化して，改めてすべてを連立させてまとめると，
2
argmin ∥ f (q)∥
q
本来，モデルや計測値に誤差がなく，パラメタが正しく求められていれば
，
f (q)=0
q0 から始めて一次のニュートン法で解く．
初期パラメタ　e 0= f (q 0 )
2015年度前学期
電気通信大学大学院　情報システム学研究科情報メディアシステム学専攻　知能システム論１
16
- ガウス・ニュートン法
e 0= f (q 0 )
q0 のまわりでテーラー展開して
( )
∂f
f (q 0+ Δ q)≈ f (q 0 )+
Δ q=0
∂q
これを解くと
+
( )
∂f
Δ q=−
∂q
e0
これを繰り返しにより収束させる．この（局所）最小二乗性は疑似逆行列
が与える解の最小二乗性によって保証される．
2015年度前学期
電気通信大学大学院　情報システム学研究科情報メディアシステム学専攻　知能システム論１
17
- 特異値分解(numpy)

vpythonでも使っている numpyモジュールに
線形代数のルーチンがある。その中のsvd関
数が特異値分解
(1) モジュールの読み込み
from numpy import *
from numpy.linalg import *
(2) 分解される行列
a_mat=empty((M,N),dtype=float64)
(3) 特異値分解の実行
u_mat,w_vec,v_trn=svd(a_mat,compute_uv=1,full_matrices=0)
2015年度前学期
電気通信大学大学院　情報システム学研究科情報メディアシステム学専攻　知能システム論１
18
- 特異値分解の使い方

b=A x
連立一次方程式　を解く
(1) 行列の準備
a_mat=empty((M,N),dtype=float64)
b_vec=empty((M),dtype=float64)
w_size=min(M,N)
w_inv=empty((w_size),dtype=float64)
(2) a_mat, b_mat の作成
for i in range(M) :
for j in range(N):
a_mat[i][j]=....
b_vec[i]= ....
(3) 特異値分解の実行
u_mat,w_vec,v_trn=svd(a_mat,compute_uv=1,full_matrices=0)
2015年度前学期
電気通信大学大学院　情報システム学研究科情報メディアシステム学専攻　知能システム論１
19
- 特異値分解の使い方
(
x=A + b

x=V W−1 U T b
)の計算
(4) 小さい特異値の切り捨て
w_max=max(w_vec) #w_max=w_vec[0] could be ok
w_min=w_max*rlim
for j in range(w_size) :
if w_vec[j]<w_min :
w_inv[j]=0.0
else :
w_inv[j]=1.0/w_vec[j]
(5)
V W−1　UT b
の計算
u_trn=u_mat.transpose()
tmp=dot(u_trn,b_vec)
tmp=tmp*w_inv
v_mat=v_trn.transpose()
x_vec=dot(v_mat,tmp)
2015年度前学期
電気通信大学大学院　情報システム学研究科情報メディアシステム学専攻　知能システム論１
20
- 次回予告

アームのヤコビアン
2015年度前学期
電気通信大学大学院　情報システム学研究科情報メディアシステム学専攻　知能システム論１
21

Download Report