行列計算を利用したデータ解析技術

九州数学教育会第 6 回算数数学教育研修会 (於. 福岡県立修猷館高等学校)

行列計算を利用したデータ解析技術
主成分分析・多次元尺度構成法

溝口佳寛

九州大学マス・フォア・インダストリ研究所

http://www.slideshare.net/yoshihiromizoguchi/

2012 年 4 月 15 日 (日)

溝口佳寛 (九大 IMI) 行列計算を利用したデータ解析技術 2012/04/15 1 / 44

目次

1 ..
自己紹介
2 はじめに
3 ..
最初の例
4 準備
5 ..
特異値分解 (SVD)
6 主成分分析 (PCA)
7 ..
多次元尺度構成法 (MDS)
8 まとめ
9 ..
参考文献

..

..

..

マス・フォア・インダストリ研究所 (正式名称)

http://www.imi.kyushu-u.ac.jp/ (産業数学研究所 (俗称))


研究者情報 (ResearchMap)

http://researchmap.jp/yoshihiro/
「Yoshihiro Mizoguchi」で検索.


プロフィール
計算について論理的・数学的に考察する研究を行っています. 1930 年代 Alan Turing は「チューリ
ング機械」と呼ばれる形式的な計算モデルを構成し, 計算可能性, 万能性に関する計算理論の礎を築
きました. 数の計算の実現に計算モデルのテープ上に記述された文字列が重要な役割を果たしました.
計算モデル「有限オートマトン」の研究では 1950 年代に入り文字列集合の認識機械としての考察が
始まり言語との重要な関係が導かれました. その後, 形式言語の研究は「自然言語の機械翻訳」「文献
データベース」「人工知能」などの研究へと発展しています.
数学 (代数学) は「群」
「環」
「体」などの一種の集合上の演算を持つ代数系を対象とします. 代数系
はモナドと呼ばれる圏の間の関手たちにより多種の集合間の関数を演算として持つ系に一般化され
ます. 多種代数系は「スタック」「木」
「グラフ」などのプログラミングに必要なデータ構造を含みま
すので, 多種代数系やモナドの研究成果は計算の性質を考察するための強力な武器になります.
近年の計算機の発展と普及により, 社会システムは計算機 (プログラム) にますます依存して来てい
ます. システムに要求される仕様をプログラムで実現し計算機で実行します. プログラムが仕様を正
しく実現していることは, 数理論理学において証明が定理を正しく導いていることに対応します. 論
理的に誤りがないプログラムのための仕様記述, 開発, 検証の技術は形式手法と呼ばれます. 証券取
引, 交通網, 航空宇宙工学, マイクロプロセッサなどプログラムの誤りにより多大な損失が生じるシス
テムに形式手法は活用されています.

最近は論理的・数学的な考察そのものにも計算機が利用され, 形式手法のための数学理論の構築と
実現がますます重要になっています. 数の取り扱いのための文字列の研究のように, 計算に必要な数
学概念の計算モデルの研究を続けています. 最近では, デジタル映像分野における計算のための新し
い計算モデルも模索しています.


デジタル映像数学の構築と表現技術の革新

CREST 研究プロジェクト (2010-2015): CG 固有の問題を数学的に解く
(MCG: Mathematics & Computer Graphics)
チーム HP は 5 月公開予定 http://mcg.imi.kyushu-u.ac.jp/


はじめに

FFT, RSA, DCT は知っている (聞いたことがある).

PCA, MDS, SVD って何だろう?


デジタルデータの話題

デジタル機器で利用される数学技術の例:

FFT (高速フーリエ変換)
J.W.Cooley and T.W.Turkey (1965).
信号処理 (通信や放送, 医療計測), 多倍精度数値計算などで利用
RSA (公開鍵暗号方式)
R.Rivest, A.Shamir, L.Adleman (1977).
インターネット・オンラインショッピング等で利用
DCT (離散コサイン変換)
N.Ahmed, T.Natarajan and K.R.Rao (1974).
デジタルカメラ等の画像ファイル圧縮 (JPEG) で利用


大量データの話題

「ビッグデータ」と呼ばれる大量データ処理で利用されている技術.
(cf. ギガ (106 ) バイト (G), テラ (1012 ) バイト (T), ペタ (1015 ) バイト (P))

PCA (主成分分析)
K.Pearson (1901).
高次元データをデータ間の特徴を表す低次元データに圧縮する方法
MDS (多次元尺度構成法)
W.S.Torgerson (1952)[3]
距離の関係がわかっているデータたちを表現する多次元データ座標
を構成する.
SVD (特異値分解)
PCA や MDS の基礎となる線形代数の定理 (行列の対角化の一般化).
C.Eckard and G.Young (1936)[1].
The approximation of one matrix by another of lower rank,
Psychometrika, 1(3). ⇐ 心理学の雑誌


最初の例 (成績処理)

国語数学英語 f (x1 , x2 , x3 )
x1 x2 x3 f (x)
= a1 x1 + a2 x2 + a3 x3
A 40 40 52 f ( A)
B 40 46 52 f (B) とするとき, a1 , a2 , a3 の値を
C 46 52 52 f (C) どのようにすれば f (x) の標準
D 46 58 52 f ( D) 偏差 (点数のばらつき) を最大
E 40 64 40 f (E) にすることが出来るか?
F 40 70 40 f (F) 但し, a2 + a2 + a2 = 1 の条件
G 46 76 40 f (G) 1 2 3
を加える.
H 46 82 40 f (H)
平均 43 61 46
標準偏差 3.21 14.7 6.41 最大化
a1 = a2 = a3 = 1
√
3
が単純な合計値に対応する.

注意: ばらつきを大きくすることと成績の善し悪しの評価とは無関係!
(例: ai の符号を全て反転しても, ばらつきは同じ.)

成績処理 (2)
国語数学英語
x1 x2 x3 f (x) g(x) h(x)
A 40 40 52 76.2 21.3 46.4
B 40 46 52 79.7 26.9 49.8
C 46 52 52 86.6 32.9 55.6
D 46 58 52 90.1 38.5 59.0
E 40 64 40 83.1 48.0 58.1
F 40 70 40 86.6 53.5 61.6
G 46 76 40 93.5 59.5 67.3
H 46 82 40 96.9 65.1 70.7
平均 43 61 46 86.6 43.2 58.6
標準偏差 3.21 14.7 6.41 6.93 15.8 8.16

1
f (x1 , x2 , x3 ) = √ (x1 + x2 + x3 )
3
g(x1 , x2 , x3 ) = 0.0795457x1 + 0.927542x2 − 0.365156x3
h(x1 , x2 , x3 ) = 0.395067x1 + 0.565842x3 + 0.154093x3

成績処理 (3)

A, B, · · · , H を座標 (g(x), h(x)) でプロットしたもの.
85

C D
80

G H
hx

75

A B
E F
70

65
20 30 40 50 60 70
gx

※ 3 科目の成績の特徴を 2 次元平面に上手にプロットする方法は?

3 次元データの特徴を維持しつつ, 2 次元データに圧縮する方法は?


記号と基礎統計量の復習
x = (x1 , x2 , · · · , x n), y = (y1 , y2 , · · · , y n), 1 = (1, 1, · · · , 1) とする.
 



y1 


 
∑
n ( ) 

 y2 



内積 s= xi yi = x, y = x · y t = x1 x2 ··· xn · 


 . 





 .
. 


i=1 
 

yn
∑
n
ノルム ||x|| = x, x = (xi )2
i=1
∑
n
合計 s= xi = 1, x
i=1
1
平均 x=
¯ 1, x
nn
∑
分散 σ2 = (xi − x)2 = x − x1, x − x1
¯ ¯ ¯
√i=1
標準偏差 σ= x − x1, x − x1
¯ ¯

データの分析, 数列 (Σ), 行列, ベクトル, · · · の勉強が必要.

データ行列 · 距離行列 · 内積行列

n 個のデータ oi = (xi1 , xi2 , · · · , xip) (i = 1, · · · , n) を考える. ここで
は, p ≤ n としておく. n × p 行列, X = (xi j ) (i = 1, · · · , n,
j = 1, · · · , p) に対して, n × n 行列 M と B を

M = (||oi − o j ||),
B = ( oi , o j ) = (oi · (o j ) t ) = X · X t

とする.
M を距離行列, B を内積行列, X をデータ行列と呼ぶことにする.
n × 1 の列ベクトルを xi = (x1i , x2i , · · · , x ni ) (i = 1, · · · , p) と書き,
X = (x1 , x2 , · · · , x p) と書くこともある.

X が n × p 行列のとき, B = X · X t は n × n 行列, X t · X は p × p 行列.


B と X の関係

問題
n × n 行列 B = (bi j ) が与えられ, X がわからないときに, 適当な次元 p を
.
定め, bi j = oi , o j (i = 1, · · · , n, j = 1, · · · , n), すなわち, B = X · X t とな
るような, n 個のベクトル oi = (xi1 , xi2 , · · · , xip) を求めよ.

B を直交対角化して UBU t = Λ と対角行列に出来たとする. このとき, B
.
は半正定値行列であり固有値は 0 または正であることが知られていて,
√
Λ (行列 Λ の対角成分に全て平方根を施した行列) を作ることが出来る.
√ √ √ √
そして, B = U t ( Λ)2 U = (U t Λ)(U t Λ) t であり, X c = U t Λ とすると
B = X c X c となる. このときの X c は n × n 行列である.
t

 b ···   u ···  λ ···  u ··· 
 11


b12 b1n   11
 
 
u21 u n1   11
 0 0   11
 u12 u1n 


 b21

 b22 b2n   u12
 
 
  u22 u n2

 0


 λ22 0

  u21


 u22 u2n






  
    
 .


 .  
. = .
 

.  .



.  .

 .






 . .   .
 
 
  .  .

 .  .

 . 





 .   .
 
 
  .  .



 .  .



 . 




b n1 b n2 ··· b nn u1n u2n ··· u nn 0 0 ··· λ nn u n1 u n2 ··· u nn


アルゴリズム 1

B を作ったもとのデータ行列 X を探す:
X = (xi j ) のランクが p であるとき, rank(B) は p であることが知られてお
り, Λ の対角成分は λ1 ≥ λ2 ≥ · · · ≥ λ p で, λ p+1 = · · · = λn = 0 と考えて
良い. Λ = (λi j ), u = (ui j ) とし, r ≤ p とするとき, r × r 行列 Λ r = (λi j )
(i, j = 1, · · · , r), n × r 行列 U r = (ui j ) (i = 1, · · · , n, j = 1, · · · , p) とする
√
と B = U pΛ pU p となる. X c = U p Λ p とすると X c は n × p 行列で,
t

B = X c X c となる.
t

 b ···   u ···  λ ···  u ··· 
 11


b12 b1n   11
 
 
u21 u p1   11
 0 0   11
 u12 u1n 


 b21

 b22 b2n   u12
 
 
  u22 u p2  0



 λ22 0   u21



 u22 u2n






  
    
 .


 .  
. = .
 

.  .



.  .

 .






 . .   .
 
 
  .  .

 .  .

 . 





 .   .
 
 
  .  .



 .  .



 . 




b n1 b n2 ··· b nn u1n u p2 ··· u pn 0 0 ··· λ pp u p1 u2 p ··· u pn

√
 u  
 11


u21 ··· u p1 



λ11
√
0 ··· 0 



 u
 12

 u22 u p2 


 0 λ22 0





 





XC =  .



 . 


 . .




 .
 . 

 . . 


 .

 . 

 . . 

 
 √ 


u n1 u n2 ··· u pn 0 0 ··· λ pp


特異値分解

定理 (Singular Value Decomposition[1])
n × p 行列 X に対して, n × n 直交行列 U, p × p 直交行列 V が存在し,
X = U · L · V t とすることが出来る. 但し, n × p 行列 L = (li j ) は対角行列
li j = 0 (i j, i = 1, · · · , n, j = 1, · · · , p) である.
そして, このとき次の条件が満たされる.
X · X t は, U で直交対角化され 0 でない固有値の集合は
{l2 | lii 0, i = 1, · · · , p}.
ii
UU t = E n で, U(XX t )U t が対角行列になる.
X t · X は, V で直交対角化され 0 でない固有値の集合は
{l2 | lii 0, i = 1, · · · , p}.
ii
VV t = E n で, V(X t X)V が対角行列になる.
rank(X) = #{lii | lii
0, i = 1, · · · , p}. (0 でない特異値の数)
.
∑ p
Tr(X · X ) = Tr(X · X) =
t t 2
lii
i=1


特異値分解 (例)

8 個の 3 次元ベクトルのデータ, すなわち, 8 × 3 行列 X を考える.
6
7
 



0 0 0 


5 


 0 1 0






 




 1 1 0 




 


2 8 
 1 0 0 

X=









 0 0 1 




 


3 




0 1 1 






 1 1 1 


 
1 1 0 1

4

行列 X の特異値分解

行列 X の特異値分解 X = U · L · V t を計算する.
 0 0 0 0 0 0 0 1 


 





1
√ 0 1
√
1 − √1 − √3 1
√ 0 




 2 





2 6 3 2 2 2 10 2 10 


   2 √2

1 −1 1 −1 1
− √ 1 − √1



√ √ √ 0 
 


 6 2 2 3 2 2 2 2 10 √ 10 

 0
√
0 




 

 





5




0 2 0
√




 1 1
−√ 1
−√



 1 2 
 
 √ 


 √ −1 1
− √ 0 0 0 0 
 0 0 2 
 √6 



 2 6 2 2 3 2 

 


3 2 



X=
 1 1 1 3 

 0 0 0 

 1 2 




 √ − √ −1 1
− √ − √ 1
√ 0 

 


√ 0 3 




 2 6 2 2 3 2 2 2 √ 10
2 2 10 

 0 0 0 


3 




 

 


1
√
1
√
1
−√ 



1
√
1 1
√ 0 0 2
√
1 0 


0 0 0 




 √6
2 2 3 5 10 

 0 0 0 


3 2 6


 

 


 3 

 0 0 0


 2 




 0 0 0 √
1 − √1 1
− √ 0 




 2 



 1 1 1
2
1 1
10 2 10 

√ 0 − √ 2
− √ √ − √1 0
6 3 2 2 2 10 10


X · X t の対角化
特異値分解の行列 U は X · X t を直交対角化する.
(※ 対称行列 X · X t は直交行列で対角化可能).
 0 0 0 0 0 0 0 1 


 





1
√ −1 1
√
1 − √1 − √3 1
√ 0 




 2 2 





2 6 2 3 2 2 2 10 2 10 





1
√ −1 − √ 1 −1 1
− √ 1
√ − √1 0 




 6 2 2 3 2 2 2 2 10 √ 10





 




 5 




 2 





1
√ 0 1
−√ 0 0 0 0 




 2 6 3 2 



U=
 1 1 1 3 




 √ √ −1 − √1 − √ 1
√ 0 




 2 6 2 2 3 2 2 2 √ 10
2 2 10 




 





1
√ 0 1
√ 0 0 2
√
1 0 




 √6 3 5 10 




 




 3 




 2 




 0 0 0 √
1 − √1 1
− √ 0 




 2 



 1 1 1 1
2
1 1
10 2 10 

√ 2
− √ 2
− √ √ − √1 0
6 2 3 2 2 2 10 10

 0 0 0 0 0 0 0 0   8 0 0 0 0 0 0 0 


 

 

 0 




 0 1 1 0 0 1 1 0 

 

 2 0 0 0 0 0 0 




 

 
 0
 




 0 1 2 1 0 1 2 1 

 

 0 2 0 0 0 0 0 



 
  

t (X · X t )U = U t 


0 0 1 1 0 0 1 1 U =  0







0 0 0 0 0 0 0 




U 

 0 0 0 0 1 1 1 1 

  0

 0 0 0 0 0 0 0 




 

 

 0 




 0 1 1 0 1 2 2 1 

 

 0 0 0 0 0 0 0 




 

 
 0
 



 0 1 2 1 1 2 3 2 
 
 0 0 0 0 0 0 0 

0 0 1 1 1 1 2 2 0 0 0 0 0 0 0 0

注: 0 以外の固有値は特異値の平方になる.

X t · X の対角化

特異値分解の行列 V は X t · X を直交対角化する.
(※ 対称行列 X t · X は直交行列で対角化可能).
 


1
√ −√
1
−√
1




 3 2 √6 



 

V=



1
√ 0 2 





 3 3 


 1
√
1
√ − 1
√

3 2 6

   
 4 2 2 

 
  8 0

 0 


   
V (X · X)V = V  2 4 2  V =  0 2
t t 


t 

 

 0 



   
2 2 4 0 0 2
注: 0 以外の固有値は特異値の平方になる.


中心化行列

定義 (中心化行列)
自然数 n に対して, E n を n × n 単位行列, 1 n を全成分が 1 の n 次元列ベク
トル, 1 n×n を全成分が 1 の n × n 行列とする. このとき, n × n 行列
J n = E n − 1 1 n×n を中心化行列という.
n .

1∑
n
1
n × p 行列 X に対して, J n X = (E n X − 1 n×n X) = (xi j − x k j ) であ
n n
k=1
∑
n .
り, J n X = (x ) とすると xi j = 0 となる. すなわち, J n X はデータ列 X
ij
i=1
の重心を原点に平行移動したデータ列になる.


主成分分析 (PCA)

定理 (Principal Component Analysis)
1
n × p 行列 X の共分散行列を p × p 行列 S = (J n X) t (J n X) とする. S を直
n
交対角化して VSV t = Λ となり, Λ = (λi j ) は, λi j = 0 (i j),
λ11 ≥ λ22 ≥ · · · ≥ λ pp ≥ 0 とする. n × p 行列 Y を Y = (J n X) · V t とする
と Y = (y1 , y2 , · · · , y p) = (yi j ) に対して, 次の性質が成り立つ.
∑n
1 yi j = 0 ( j = 1, · · · , p), (yi の平均は 0)
i=1
1 t
2 Y Y = Λ (Y の共分散行列は対角行列で対角成分 (yi の分散) は固有値),
n .
3 Tr(Y t Y) = Tr(S),
4 |Y t Y| = |S|.

変換されたデータ列 Y について, 次のことがわかる. (1) より各成分の平均は全て
0 である. (2) より各成分の分散は第 1 成分から順番に小さくなっている.


PCA (例)

8 個の 3 次元ベクトルのデータ, すなわち, 8 × 3 行列 X を考える.

8
7  
6 


0 0 0 


5 

 





0 1 0 




 





1 2 0 




 1 3 0 


X=









 0 4 2 




 


4 

 0 5 2 




 


3 


1 6 2 


2 1 7 2
1


中心化行列 (例)

データ X の重心は ( 1 , 7 , 1) である.
2 2

 7 −1 1
−8 −1 −1 −1 −1 −1   −1 −7 −1 


 8 8 8 8 8 8 8 

  

 



 
   −1
2 2 
  
  0 0 0
−1   
7 −1 −1 −1 −1 −1 −1
  −5 −1




8 8 8 8 8 8 8 8 


 0 1 0  
  12
 


2 






 −1 −1 7 −1 −1 −1 −1 −1 

  
  2
  −3 −1 




 8 8 8 8 8 8 8 8 

 1 2 0   1
 
  2 




 −1 −1 1
−8 7 −1 −1 −1 −1 

  
 
  2 −1 −1 



J8 X = 
 8 8 8 8 8 8 8 


1 3 0 =
 
  2 




 −1 −1 1
−8 −1 7 −1 −1 −1 

 0 4 2   −1
 
  1 




 8 8 8 8 8 8 8 

  
 
  2 2
1 




 

 0 5 2  
  −1
  3 



 −1 −1 1
−8 −1 −1 7 −1 −1 
  
  1 



8 8 8 8 8 8 8 

 1 6 2  
  12
 
2 





 −1 −1 1
−8 −1 −1 −1 7 −1 

 
 2

5 1 



 8 8 8 8 8 8 8 
 1 7 2 
 1 2 

7
−1
8
−18
1
−8 −1
8
−1
8
−18
−1
8
7
8 2 2
1

原点を重心に平行移動したときの座標に変換することが行列
1
J 8 = E8 − 18×8
8
を掛けることに対応する.


共分散行列

共分散行列 S を計算する.
 −1 −7 −1 


 





2 2 



 −1 −5 −1 




2 2 


 
 1 3
−2 −1 
  1
 −1 −1 1 1 −1 −1 1 1 

 
 
  1 


2   4


 2 
 1 1
−2 −1  
  0 


1 1  2
 −7
2 2 2 2 2 2 

7  
= 1
 
2 

S = (J 8 X) t (J 8 X) = 

 −5 −3 −1 1 3 5 


2  
  2
21 2 


8 8 
 2 2 2 2 2 2 2 2 

 −1 1 1  

  0 4 

−1 −1 −1 −1 1 1 1 1 
 2 2 

 2 1


 −21 3 1 




 2 




 1 5 




 2 2
1 


1 7 1
2 2

共分散行列 S を直交対角化する.
   
   0.0795457
1 1
 0.0795457


0.927542 0.365156 


0 

 −0.710589 0.699097 


 
4 2
 −0.710589

 −0.204143 
 1 21   0.927542
 −0.204143 −0.313037 



 0.673342 

 2 

 

0.699097 −0.313037 0.642862  2 4  0.365156 0.673342 0.642862
0 2 1
 
 6.08024


0 0 


= 

 



 0 0.393643 0 

0 0 0.026113


共分散行列の固有ベクトル
※ 固有ベクトル: (0.0795457, 0.927542, 0.365156),
(−0.710589, −0.204143, 0.673342), (0.699097, −0.313037, 0.642862).
8
7
6
5 Y = (J 8 X) · V
 
 −3.65132


0.396452 0.103221 


 −2.72378

 0.192309 −0.209817 




 −1.7167 




 −0.722422 0.176243 




 −0.789154 




 −0.926565 −0.136795 


= 
 0.789154
 




 0.926565 0.136795 




 1.7167 



 0.722422 −0.176243 

4 
 2.72378
 




 −0.192309 0.209817 


3 3.65132 −0.396452 −0.103221
2
1

各列データの平均は全て 0, 分散は 6.94885, 0.449878, 0.0298434.


M と B の関係
問題
n 個のデータ oi = (xi1 , xi2 , · · · , xip) からなる n × p 行列 X = (xi j ) により
.
作られた n × n の内積行列 B = X · X t = ( oi , o j ) と n × n の距離行列
M = (||oi − o j ||) を考える. データ行列 X がわからなくても, 内積行列 B
から距離行列 M を直接計算することが出来るか?

距離行列 M = (mi j ) に対して, 全成分を 2 乗した行列 M(2) = (m2 ) を考え
ij
.
る. (注. M(2) は各成分が 2 乗であり, 行列の 2 乗 M2 とは異なる!)

m2j = ||oi − o j ||2 = oi − o j , oi − o j = oi , oi − 2 oi , o j + o j , o j
i

であるので, B = ( oi , o j ) がわかっていれば, M(2) を計算することが出来
る. 従って, B から M も直接計算することが出来る.

問題
M から B は直接計算出来るか?

多次元尺度構成法 (MDS)

n 個のデータ oi = (xi1 , xi2 , · · · , xip) からなる n × p 行列 X = (xi j ) により
作られた n × n の距離行列 M = (||oi − o j ||) を考える.

問題 (Multidimensional Scaling)
n × n 行列 M = (mi j ) が与えられ, X がわからないときに, 適当な次元 p を
定め, mi j = ||oi − o j || (i = 1, · · · , n, j = 1, · · · , n) となるような, n 個のベ
クトル oi = (xi1 , xi2 , · · · , xip) を求めよ.

.



M を作ったもとのデータ行列 X を探す:

1
B c = − J n M(2) J n
t
2
とすると
B c = (J n X)(J n X) t
であることが計算出来る [4].
距離行列 M から内積行列 B c を得る変換 − 1 J n M(2) J n を
2
t

Young-Householder 変換と呼ぶ.
この B c にアルゴリズム 1 を使えばもとのデータ X(を中心化したもの) を見つけ
ることが出来る.


次元縮小

アルゴリズム 1 は, n × p のデータ行列から作られた内積行列 B = X · X t
が与えられたときに, B を実現する X を求めるアルゴリズムであった.
丁度 B = X · X t である必要はないので, より小さな次元 r(< p) のデータ
行列 X で B を近似出来ないだろうか?

問題
準備で述べた n × n 行列 B = (bi j ) が与えられ, X がわからないとする. 自
然数 r ≤ rank(B) が与えられたとき, Tr((B − X c X c )2 ) を最小にする n × r
t

行列 X c を求めよ.

溝口佳寛 (九大 IMI) 行列計算を利用したデータ解析技術 . 2012/04/15 31 / 44


B を作ったもとのデータ行列 X を近似する:

アルゴリズム 1 で, p を r < p に置き換えて,
√
Xc = Ur Λr

とすると X c は n × r 行列となりもとの X よりは次元が小さなデータにな
る. 但し, B = X c X c とは限らない.
t

 b ···   u ···  λ ···  u ··· u1n 
 11


b12 b1n   11
 
 
u r1 u p1   11


0 0   11


u12 


 b
 21
 b22 b2n   u
  12
  u r2 u p2 
 .
 
 .
 



  
 
  
 . . 
 . 

 .  . . .  . . 


 . . = .
 
  . .  .

 . .

 .
 . 




 . .   .
 
  .


. 


 





 .   .
 
 
 
  .

.  0


 ··· λ rr 0  u
  r1


 u r2 ··· u rn 




b n1 b n2 ··· b nn u1n ··· u rn u pn 0 ··· 0 λ pp u p1 u p1 ··· u pn


理論背景

定理
n × p 行列 X から作られた距離行列 M に対して, n × r 行列 X c の中で Φ
を最小にするのは, アルゴリズム 3 で得られる X c である. そして, Φ の最
∑
n
小値は 2n(λ r+1 + · · · + λ p) となる. 但し, Φ = (d2j − di 2 ) で, di j は X を
i j
i, j=1
.
構成するデータ間の距離 di j = ||oi − o j || であり, d は X c を構成するデー
ij
タ間の距離とする.

定理
n × p 行列 X から作られた内積行列 B に対して, n × r 行列 X c の中で Ψ
を最小にするのは, アルゴリズム 3 で得られる X c である. そして, Ψ の最
小値は (λ2 + · · · + λ2 ) となる. 但し, Ψ = Tr((B − X c X c )2 ) とする.
p
t
r+1 .
定理の証明等は, 文献 [2] を参照のこと.


アルゴリズム 2 (例) (距離行列 M)

ある 8 個のデータ oi = (xi1 , · · · , xip) (i = 1, · · · , 8, p < 8) により構成され
た, 距離行列 M(= (||oi − o j ||) が与えられているとする. データ行列
X = (xi j ) (i = 1, · · · , 8, p = 1, · · · , p) はわからないことにする.
 √ √ √ √ 


 0 1 2 1 1 2 3 2 




 √ √ √ √ 




 1 0 1 2 2 1 2 3 




 √ √ √ √ 




 2 1 0 1 3 2 1 2 




 √ √ √ √ 




 


M= 
 1 2 1 0 2 3 2 1 

 √ √ √ √ 



 1 2 3 2 0 1 2 1 




 √ √ √ √ 




 




 √
2 1
√
2
√
3 1
√
0 1 2 




 




 √
3
√
2 1
√
2 2 1
√
0 1 


 
2 3 2 1 1 2 1 0


アルゴリズム 2 (例) (内積行列 B)

距離行列 M に Young-Householder 変換を施し内積行列 B を得る.

1
B = − J n M(2) J n
t
2
1 1 1
= − (E n − 1 n×n) · M(2) · (E n − 1 n×n) t
2 n n
 

 0.75 0.25 −0.25 0.25 0.25 −0.25 −0.75 −0.25 


 0.25
 




 0.75 0.25 −0.25 −0.25 0.25 −0.25 −0.75 


 −0.25 0.25

 0.25 −0.75 −0.25 −0.25 





0.75 0.25 


 0.25 −0.25 0.25

 0.75 −0.25 −0.75 −0.25 0.25 


= 

 0.25 −0.25 −0.75 −0.25 0.75 




 0.25 −0.25 0.25 




 −0.25 0.25 −0.25 −0.75 0.25 




 0.75 0.25 −0.25 



 −0.75 −0.25 0.25 −0.25 −0.25
 




 0.25 0.75 0.25 


−0.25 −0.75 −0.25 0.25 0.25 −0.25 0.25 0.75


アルゴリズム 2 (例) ( B の対角化)

内積行列 B を直交対角化する.
 −0.45 −0.088 −0.54 



0 0.088 0.54 0.45 0 




 0.36 0 −0.57 −0.21 0.57 0.21 −0.36 0 




 




 0.2 0.61 0.2 −0.2 −0.2 0.2 −0.2 −0.61 




 


 −0.16 −0.16 −0.16 
G= 
 0.16 0.47 0.16 0.16 0.79 


 



 −0.63 −0.08 0.14 0.13 0.058 0.21 −0.71 0 



 




 −0.27 0.56 −0.12 0.56 0.44 0.0026 0.29 0 




 −0.099 −0.2 −0.098 



 0.13 0.73 0.63 0.027 0 

0.33 −0.26 0.19 0.48 −0.07 0.74 0.071 0

とすると
 



2. 0 0 0 0 0 0 0 




 0 2. 0 0 0 0 0 0 




 




 0 0 2. 0 0 0 0 0 




 


 
G BG =  
t  0 0 0 0 0 0 0 0 


 





0 0 0 0 0 0 0 0 




 





0 0 0 0 0 0 0 0 




 0 0 0 0 0 0 0 0 


 
0 0 0 0 0 0 0 0


アルゴリズム 2 (例) ( B の対角化 (つづき))

 



2. 0 0 0 0 0 0 0 




 0 2. 0 0 0 0 0 0 




 




 0 0 2. 0 0 0 0 0 




 


t  
G · ·G
 0 0 0 0 0 0 0 0 
B = 

 





0 0 0 0 0 0 0 0 




 





0 0 0 0 0 0 0 0 




 0 0 0 0 0 0 0 0 


 
0 0 0 0 0 0 0 0
 −0.45 0.2 



0.36 




 0 0 0.61 


 




 −0.088 −0.57 0.2   √2

 
 



  
−0.2  
   0 0 

 −0.54 −0.21  
· √ 


= 

   0
  



 0.088 0.57 −0.2  
 
 

2
√
0 



 



 0.54 0.21 0.2 



0 0 2


 −0.36 −0.2 


 0.45 

0 0 −0.61
 √   


 2 0 0  
  −0.45
  0 −0.088 −0.54 0.088 0.54 0.45 0 


 √  
 
 ·  0.36 

· 0


 2 0 
 
  0 −0.57 −0.21 0.57 0.21 −0.36 0 




 √  
 
0 0 2 0.2 0.61 0.2 −0.2 −0.2 0.2 −0.2 −0.61

(注意: G は最初の 3 行しか必要ない.)

アルゴリズム 2 (例) (データ行列 X)

 −0.45 



0.36 0.2 




 0 0 0.61 




 
  √


 −0.088 −0.57 0.2 
  


  
 
  2 0 0 


 −0.54 −0.21 −0.2  
  √ 


= 

 ·
  0
 

X 
 −0.2  
  2 0 





0.088 0.57  

 √ 


 

 0 0 2



0.54 0.21 0.2 




 0.45 −0.36 −0.2 


 
0 0 −0.61
 −0.64 



0.51 0.29 




 0. 0. 0.87 




 




 −0.12 −0.81 0.29 




 




 −0.76 −0.3 −0.29 


= 

 



 0.12 0.81 −0.29 



 




 0.76 0.3 0.29 




 −0.51 −0.29 



 0.64 

0. 0. −0.87

とすると, B = X · X t となる.
与えられた内積行列 B を実現するデータ X を求めることが出来た.


アルゴリズム 2 (例) ( M の確認)

   



x11 x12 x13   −0.64 0.51
 
 
0.29 




 x21 x22 x23   0.
 
  0. 0.87 




  
 
  



 x31 x32 x33   −0.12 −0.81
 
  0.29 





  
 
  −0.76 −0.3 



 x41 x42 x43  
  −0.29 

X=



=
 
 
  0.12






 x51 x52 x53  
 
  0.81 −0.29 




  
  0.76
  




x61 x62 x63  
 
 
0.3 0.29 





 x71 x72 x73   0.64 −0.51
 
  −0.29 



  
  

x81 x82 x83 0. 0. −0.87
に対して
√
||(xi1 , xi2 , xi3 ) − (x j1 , x j2 , x j3 )|| = (xi1 − x j1 )2 + (xi2 − x j2 )2 + (xi3 − x j3 )2

を計算すると最初の距離行列 M の i, j 成分に一致することが確認出来る.

与えられた距離行列 M を実現するデータ X を求めることが出来た.


アルゴリズム 2 (例) (距離行列 M(再び))

実は M は 1 辺の長さが 1 の立方体の 8 個の頂点間の距離行列である.

6
7
 √ √ √ √ 






0 1 2 1
√
1
√
2
√
3
√
2 





5


 1 0 1 2 2 1 2 3 



 √ √ √ √ 



 2 1 0 1 3 2 1 2 




 √ √ √ √ 




 1 2 1 0 2 3 2 1 

 2 8
M=


 √ √ √ √ 





 1 2 3 2 0 1 2 1 




 √ √ √ √ 



 







√
2 1
√
2
√
3 1
√
0 1 2 





3


 √
3
√
2 1
√
2 2 1
√
0 1 


2 3 2 1 1 2 1 0
1

4


アルゴリズム 3 (例)
G は 3 行まで減らすことが出来たが, さらに減らして 2 行にしてみる.
 −0.45   −0.64 



0.36 

 


0.51 




 0 0 

 

 0. 0. 




 

 

 




 −0.088 −0.57   √

  

 −0.12 −0.81 



  
    

 −0.54 −0.21  
 0  
=
 −0.76 −0.3 


X= · √   
   2


  

  
 









0.088 0.57 

 0 2 


0.12 0.81 




 

 

 





0.54 0.21 

 


0.76 0.3 




 0.45 −0.36 

 

 0.64 −0.51 


   
0 0 0. 0.

2 つの値を x- y 座標と考え接続情報を加えて 2 次元平面に描く.

5

1

6

8
2

4

7

3


PCA と MDS の関係
n × p 行列 J n X を特異値分解して,
J nX = U · L · V t
となったとする. n × n 行列 U は (J n X)(J n X) t を対角化し,
U(J n X)(J n X) t U t = Λ n
となり, p × p 行列 V は (J n X) t (J n X) を対角化し,
V(J n X) t (J n X)V t = Λ p
となる. V は,
1
S=
(J n X) t (J n X)
n
も直交対角化するので, Y = (J n X)V は, 主成分 (PCA) である.
このとき,
Y = (J n X)V = (U · LV t )V = U · L
となる. √
ここで, U · L = U p Λ p なので, Y はアルゴリズム 1 の X c と一致する.

まとめ

「数学」は「人」によって作られる.
「人」から「人」に伝えられる物語がある.
新しい「数学」が出来るまでの「物語」(誰がいつ)
新しい「数学」技術が確立するまでの「物語」(動機)
「数学」そのものだけでなく「物語」を伝えたい.

私は大学ではコンピュータ上にプログラムとして新しい数学技術を実現し
たり, 実現されたプログラムを利用して別の数学の問題を考えたりします.
一方で, 大学生に「数学」や「プログラミング」を教育しています. 大学
生には「物語」を語るようにしています.
注: 数学の計算を計算機を相手にプログラミングで実現するときは物語は
不要です. アルゴリズムと数式を忠実にプログラム言語で記述することが
肝要です.


参考文献
[1] C. Eckart and G. Young.
The approximation of one matrix by another of lower rank.
Psychometrika, 1(3):211–218, 1936.
[2] K.V. Mardia, J.T. Kent, and J.M. Bibby.
Multivariate Analysis.
Academic Press, 1994.
[3] W.S. Torgerson.
Multidimensional scaling: I. theory and method.
Psychometrika, 17(4):401–419, 1952.
[4] G. Young and A.S. Householder.
Discussion of a set of points in terms of their mutual distances.
Psychometrika, 3(1):19–22, 1938.
[5] 小西貞則.
多変量解析入門.
岩波書店, 2010.

行列計算を利用したデータ解析技術

Recommended

More Related Content

What's hot (20)

Similar to 行列計算を利用したデータ解析技術 (20)

More from Yoshihiro Mizoguchi (20)

Recently uploaded (7)

行列計算を利用したデータ解析技術