東京都市大学 データ解析入門 3 行列分解 2

固有値問題の復習
主成分分析 (principal component analysis)
固有顔 (eigenface)
次元打ち切り（dimensional truncation）
ランダム化特異値分解（randomized SVD）
講師： 東京都市大学 田中宏和
講義ビデオ： https://www.youtube.com/playlist?list=PLXAfiwJfs0jGOvZnwUdAykZvSdRFd7K2p

1. 大規模データ解析応用事例
３. 行列分解２
情報工学部 知能情報工学科 田中宏和

2. 講義スケジュール
1.講義概要 ＆MATLAB入門
2.行列分解１：特異値分解、行列近似、最小二乗法、擬逆行列
3.行列分解２：主成分分析、固有顔、次元打ち切り、ランダム化SVD
4.スパース性と圧縮センシング１：フーリエ変換、圧縮センシング
5.スパース性と圧縮センシング２：スパース回帰、スパース分類、RPCA
6.回帰分析とモデル選択１：線形回帰、非線形回帰、数値最適化
7.回帰分析とモデル選択２：モデル選択、交差検証法、情報量基準
8.クラスタリングと分類分析１：特徴抽出、クラスタリング法
9.クラスタリングと分類分析２：教師あり学習、分類分析
10.ニューラルネットワーク1：パーセプトロン、誤差逆伝播法
11.ニューラルネットワーク2：確率勾配法、深層ネットワーク
12.発展学習：神経データ解析

3. 行列分解２：特異値分解（Singular value decomposition, SVD）
1.1 Overview
1.2 Matrix Approximation
1.3 Mathematical Properties and
Manipulations
1.4 Pseudo-Inverse, Least-Squares, and
Regression
1.5 Principal Component Analysis (PCA)
1.6 Eigenface Example
1.7 Truncation and Alignment
1.8 Randomized PCA
1.9 Tensor Decomposition and N- Way Data
Analysis

4. 行列分解２：特異値分解（Singular value decomposition, SVD）
% 1.2 SVD Example: Image Compression
CH01_SEC02_production.m
% 1.4 Psudo -Inverse, Least-Squares, and Regression
CH01_SEC04_1_Linear.m
CH01_SEC04_2_Cement.m
CH01_SEC04_3_Housing.m
% 1.5 Principal Component Analysis
CH01_SEC05_1_PCAGaussian.m
CH01_SEC05_2_OvarianCancer.m
% 1.6 Eigenface Example
CH01_SEC06_1.m
CH01_SEC06_2_3_4.m
CH01_SEC06_2_3_4_production.m
% 1.7 Truncation and Alignment
CH01_SEC07_1_Truncation.m
CH01_SEC07_1_production.m
CH01_SEC07_2_Truncation.m
CH01_SEC07_2_production.m
CH01_SEC07_3_Alignment.m
CH01_SEC07_3_production.m
% 1.8 Randomized PCA
CH01_SEC08_RSVD.m

5. 行列分解：特異値分解（Singular value decomposition, SVD）
1.特異値分解（Singular value decomposition ）
2.行列近似（Matrix approximation ）
3.応用例：画像圧縮（Image compression ）
4.擬逆行列、最小二乗、回帰問題
5.主成分分析、固有顔
6.次元打ち切り（dimensional truncation ）
7.ランダム化特異値分解（randomized SVD）
5/26
行列分解１
6/2
行列分解２

6. 【本日の内容】特異値分解（Singular value decomposition, SVD）
1.特異値分解、主成分分析、固有顔（ Eigenface）
-主成分分析と特異値分解の関係
-固有顔による顔分析
2.次元打ち切り（dimensional truncation ）
-打ち切り次元 rの決め方 二通り
-エネルギーによる打ち切り方
-ガウス雑音モデルに基づく打ち切り方
3.ランダム化特異値分解（randomized SVD）
-大規模行列の効率的な特異値分解法
-ランダム射影行列による低次元化

7. 【本日の内容】特異値分解（Singular value decomposition, SVD）
1.特異値分解、主成分分析、固有顔（ Eigenface）
-主成分分析と特異値分解の関係
-固有顔による顔分析
2.次元打ち切り（dimensional truncation ）
-打ち切り次元 rの決め方 二通り
-エネルギーによる打ち切り方
-ガウス雑音モデルに基づく打ち切り方
3.ランダム化特異値分解（randomized SVD）
-大規模行列の効率的な特異値分解法
-ランダム射影行列による低次元化

8. 【線形代数の復習】固有値問題とは
固有値問題の定義
あるn×n行列Cが与えられたとき、以下の固有方程式
を満たすn次元ベクトル v
iを固有ベクトル（eigenvector ）、λ
iを固有値（eigenvalue ）と
と呼ぶ。
i ii
λ=Cv v
対称行列の場合
行列Cがn×n対称行列である場合、すなわち
の場合、行列Cの固有ベクトル v
iは直交系をなし、λ
iは非負である。
=CC

12
1
,
0
0
i j ij n
j
ij
i
δ λλ λ
=
= ≥≥≥≥
≠
=

vv 


9. 固有値問題の具体例
固有値問題の例１
20
01

=


C
固有値問題の例２
11
1
,2
0
λ

= =


v
22
0
,1
1
λ

= =


v
31
22
31
22


=


C
11
2
2
,2
2
2
λ


= =



v
22
2
2
,1
2
2
λ
−

= =



v
1
1
2
1
1
v
2
v
2
Cv
1
Cv
1
1
1
v
2
v
2
Cv
1
Cv
x
x x
x
yy
y y

10. 行列対角化としての固有値問題
2
11
22
1
n nn
λ
λ
λ
=
=
=
Cv v
Cv v
Cv v

[ ][ ]
1
2
122
1
00
00
00
nn
n
λ
λ
λ


=



Cvv v vv v





n個の固有方程式
V V
Λ
=CV VΛ
対称行列Cの固有方程式
=V ΛCV

対称行列Cの対角化
Uは直交行列

11. MATLABを用いた固有値問題の解法
固有値問題の例１
20
01

=


C
固有値問題の例２
11
1
,2
0
λ

= =


v
22
0
,1
1
λ

= =


v
31
22
31
22


=


C
11
2
2
,2
2
2
λ


= =



v
22
2
2
,1
2
2
λ
−

= =



v
>> C = [2 0; 0 1];
>> [V, D] = eig(C); % 固有ベクトルVと固有値D
>> V’*C*V
ここで演習
左の行列の固有ベクトルと固有値を、MATLAB を用
いて数値的に求めよ。

12. 特異値分解＝固有値分解である
ˆˆ
=XUΣV
 =CV VΛ
特異値分解 固有値分解
ˆˆ
=XUΣV

ˆˆ
=XVΣU

2ˆˆˆˆ ˆ
= =XX UUVΣ ΣV VΣV
  
()
2ˆ
=VX VX Σ

とみなすと、行列Xの右特異ベクトル行列Vと特異値行列Σ を用いて、
行列Cの固有ベクトルは行列V、固有値行列はΣ
2
と書ける。
=C XX

C
Λ

13. (特異値)
2
=(固有値)
ˆˆ
=XUΣV
 =CV VΛ
特異値分解 固有値分解
2ˆ
=ΣΛ
2
11
2
2
2
2
0000
0000
0000
mm
λσ
λσ
λσ 
 
 
=
 
 




  

特異値行列
固有値行列

14. 主成分分析と固有値問題：右固有ベクトル
2ˆ ˆˆ ˆˆˆ ˆ
= =X Σ UΣVV ΣXU U U
  
( )
2ˆˆˆ
=XX U UΣ

XX

1 m
uu=
2
1
2
m
σ
σ

1 m
uu

15. 主成分分析と特異値分解
2ˆˆˆˆ ˆ
= =XX UUVΣ ΣV VΣV
  
()
2ˆ
=V X VX Σ
 
XX

1

m
v v
1

m
v v
2
1
2
m
σ
σ

=

16. 特異値分解と主成分分析
()
2ˆ
=VX VX Σ

( )
2ˆˆ
=XX U UΣ
 ( )
2
i ii
σ=XX u u

()
2
i ii
σ=XXv v

( )
1
12
2 21 22 2
12
1 11 1
2
mmm
m
m
m m m mm
×
  
  
  
= = =
  
  
  
  
∈
x xx xx xx
x xx xx xx
C XX x x x
x xx xx xx



  

  
  

  

17. 共分散と共分散行列
=C XX

() ()cov ,
ij ijij
xx= =C xx

データx の共分散行列
データx
iとx
jの共分散
（注）通常、共分散は平均値からのずれとして定義される。
平均値を差し引いた値 を新たに と定義し直せば、
と書ける。ここでは平均値を差し引いた値を用いている。
( )( )( )cov ,
jij i
−=−xxxx x x

i
−xx
i
x
( )cov ,
i j ij
=xx xx


18. 主成分分析（Principal Component Analysis, PCA ）
主成分分析 ：データ から変動の大きさに応じて成分を分析する解析法
各データ に係数 で重み付けした を考える。
x
1
i
m
i
i
v
=
= =∑y x Xv
{}
i
x {}
i
v y
の分散は以下で与えられる。y
()var= =y y y v Cv

主成分分析の計算法： ノルムが１であるようなベクトルv で、上記の分散
var(y)を最大にするものを求めよ。
1
ˆargmax
=
=
v
v v Cv


19. 主成分分析（Principal Component Analysis, PCA ）
主成分分析で解くべき問題
もしくは
1
ˆargmax
=
=
v
v v Cv

maximize subject to 1=v Cv v v

解法：固有値方程式
もしくは
λ=Cvv
=CV VΛ
( ) ( )
11
diag , ,,
mm
λλ==Vv v Λ


20. 主成分分析の図形的解釈
主成分分析で解くべき問題
maximize subject to 1=v Cv v v

具体例（二次元）
二次元ベクトル とし、行列を と取ると、問題(*) は、
10
4
01

=


C
(*)
1
2
v
v

=

v
2
21
2
2
v
v

+


22
12
1vv+=
の制約条件のもと、 を最大化せよ
1
v
2
v
1
v
2
v
1
v
2
v
1
v
2
v

21. 主成分分析の例：ガウス分布
% データの生成
xC= [2; 1;]; %平均
sig = [2; .5;]; %標準偏差
theta = pi/3; %
R = [cos(theta) -sin(theta); % pi/3回転行列
sin(theta) cos(theta)];
nPoints= 10000; % 10,000のデータ点
X = R*diag(sig)* randn(2,nPoints) + diag (xC)*ones(2,nPoints);
% 主成分分析
Xavg= mean(X,2); % 平均値
B = X - Xavg*ones(1,nPoints); % データから平均値を差し引く
[U,S,V] = svd(B/sqrt(nPoints), ‘econ’); % 主成分分析（svd）
% [U,V,S2] = pca(X’); % 主成分分析（pca）
% [U, S2]= eig(cov(X’)); % 主成分分析（eig）
CH01_SEC05_1_PCAGaussian.m
目的：二次元正規分布のサンプル点から、主成分分析を用いて主軸を見つけ
る。
データ点
主成分１ U(:,1)
主成分２ U(:,2)

22. 主成分分析の例：ガンに関わる遺伝子データ
load ovariancancer; % データの読み込み
[U,S,V] = svd(obs,‘econ‘);% 特異値分解
% 主成分への射影（三次元プロット）
figure, hold on
fori=1:size(obs,1)
% 被験者i の遺伝子データを3 次元に射影
x = V(:,1)'*obs(i,:)';
y = V(:,2)'*obs (i,:)';
z = V(:,3)'*obs(i,:)’;
if(grp{i}=='Cancer’)
plot3(x,y,z,'rx','LineWidth',2);
else
plot3(x,y,z,'bo','LineWidth',2);
end
end
遺伝子(4000)
被験者
(216)
CH01_SEC05_2_OvarianCancer.m
目的：被験者一人当たり4000次元の遺伝子データ
→ 主成分分析で3次元に次元縮約
患者
121名
健常者
95名

23. 主成分分析の例：ガンに関わる遺伝子データ
(216)
被験者
(216)
CH01_SEC05_2_OvarianCancer.m
目的：被験者一人当たり4000次元の遺伝子データ
→ 主成分分析で3次元に次元縮約
遺伝子(4000)
被験者
(216)
(3)
被験者
(216)
216 4000×
∈=XΣVU 
 216 216×
=∈X ΣVU  ()
216 3
:,1:3
×
∈UΣ 
x = V(:,1)'*obs(i,:)'; y = V(:,2)'*obs(i ,:)';
z = V(:,3)'*obs(i,:)’;

24. ここで演習
主成分分析のコードを走らせて、コードの中身を理解してください。
CH01_SEC05_1_PCAGaussian.m
CH01_SEC05_2_OvarianCancer.m

25. 主成分分析の例：固有顔（Eigenface ）
192
168
32256
ベクトル化

26. 主成分分析の例：固有顔（Eigenface ）
load ../DATA/allFaces.mat
% We use the first 36 people for training data
trainingFaces= faces(:,1:sum(nfaces (1:36)));
avgFace= mean(trainingFaces,2); % size n*m by 1; 36名の顔の平均を計算
% compute eigenfaces on mean- subtracted training data
X = trainingFaces -avgFace*ones(1,size(trainingFaces,2)); % データから平均顔を引く
[U,S,V] = svd(X,‘econ’ ); % size n*m by 1; エコノミーSVDを計算
figure, axes('position',[0 0 1 1]), axis off
imagesc(reshape(avgFace,n,m )), colormap gray
% 固有顔の最初50成分を表示
fori=1:50 % plot the first 50 eigenfaces
pause(0.1); % wait for 0.1 seconds
imagesc(reshape(U(:,i),n,m)); colormap gray ;
end
CH01_SEC06_2_3_4_production.m 固有顔の計算

27. 主成分分析の例：固有顔（Eigenface ）
第1主成分 u
1 第2主成分 u
2 第3主成分 u
3 第4主成分 u
4 第5主成分 u
5
ˆˆ
=XUΣV
データ行列の特異値分解 ˆ
=U 1
u
2
u
m
u左特異ベクトル行列
左特異ベクトルは、顔を特徴づける基底ベクトルである。

28. 主成分分析の例：固有顔（Eigenface ）
%% Now show eigenface reconstruction of image that was omitted from test set
testFace= faces(:,1+sum(nfaces(1:36))); % first face of person 37
axes('position',[0 0 1 1]), axis off
imagesc(reshape(testFace,n,m)), colormap gray
testFaceMS= testFace-avgFace;
forr=25:25:2275
reconFace = avgFace + (U(:,1:r)*(U(:,1:r)'*testFaceMS));
imagesc(reshape(reconFace,n,m)), colormap gray
title(['r=',num2str(r,'%d')]);
pause(0.1)
end
CH01_SEC06_2_3_4_production.m 固有顔による被験者37の顔再構成
approx average rr
= +x x UUx


29. 主成分分析の例：固有顔によるターゲット顔の再構成
平均顔 r=60 r=120 r=180 r=240 r=300 ターゲット顔
() ()
approx average
1
average 1
average 1 1


r
rr
rr
r
= +


= +


 
+ += +
x x UUx
u
xuu x
u
x uxu uxu







30. ここで演習
固有顔のコードを走らせて、コードの中身を理解してください。
CH01_SEC06_2_3_4_production.m

31. 【本日の内容】特異値分解（Singular value decomposition, SVD）
1.主成分分析、固有顔
-主成分分析と特異値分解の関係
-固有顔による顔分析
2.次元打ち切り（dimensional truncation ）
-打ち切り次元 rの決め方 二通り
3.ランダム化特異値分解（randomized SVD）
-大規模行列の効率的な特異値分解法

32. 特異値分解による行列近似（復習）
X= 1
u
1
v

1
σ
r
u
r
v

r
σ++
1r+
u
1r+
v

1r
σ
+ m
u
m
v

m
σ+++
=
1
u
1
v

1
σ
r
u
r
v

r
σ++
X

特異値分解における行列近似とは・・・
寄与の小さい成分 (r+1:m) を打ち切る
（truncate）こと
ではどれくらいの成分を残すべきか？
（＝打ち切り次元 rをどう選ぶべきか？）

33. 特異値分解による行列近似（復習）
2
2
F
1
m
j
j
σ
=
=∑X
1
2
2
F
j
r
j
σ
=
=∑X

データ行列 の変動
近似行列 の変動X

X
（データ行列 の変動）X
（近似行列 の変動）X

1
2
1
2
r
m
j
j
j
j
σ
σ
=
=
=
∑
∑

34. 特異値分解による行列近似
問題：打ち切り次元をどのように決めるべきか？
近似精度と圧縮率のトレードオフ
大きいr → 近似精度は高いが、圧縮率は低い
小さいr → 圧縮率は高いが、近似精度は低い
解法１：元の行列を十分よく近似できる程度に大きいr をとる
解法２：元の行列に含まれる雑音成分の大きさを推定し、雑音成分だけ
取り除く程度に大きいr をとる

35. 計算法：次の不等式
を満たすような最小の r を見つける。
打ちきり次元（r ）の決め方１
方針：近似行列が、元のデータ行列の変動の大部分（例え ば90%）を説明
できるように、打ち切り次元 r を決める。
（データ行列 の変動）X
（近似行列 の変動）X

1
1
0.9
j
j
r
j
m
j
σ
σ
=
=
= >
∑
∑

36. 方針：観測データ X
dataが、（ノイズのない）きれいな成分X
cleanとノイズ
成分X
noiseの線形和であると仮定する。
ノイズ成分から期待される特異値の最大値より十分大きいものを
取り出せば、残された成分が X
cleanであろう。
図を用いた説明 ⇒次スライド
打ちきり次元（r ）の決め方２：Gavish -Donohoの方法
data clean noise
γ= + XXX

37. 打ちきり次元（r ）の決め方２：Gavish -Donohoの方法
data
X
clean
X
noise
γX= +
= +

38. 打ちきり次元（r ）の決め方２：Gavish -Donohoの方法
この論文のやっていること
・与えられたデータ行列X
dataから、ノイズ行列X
noiseの標準偏差γ を推定
・その標準偏差γ から、期待されるノイズ部分の特異値の最大値を推定
・ノイズ部分の特異値の最大値を閾値として、その閾値より大きい特異成分を残す

39. 打ちきり次元（r ）の決め方２：Gavish -Donohoの方法
閾値τによる打ち切り次元r
・特異値を降順に並べて、ある閾値τ より大きい最小の特異値を見つけ、その固有
値をr番目とする
・r番目までの特異値成分を残し、それ以降の特異値成分を切り捨てる
τ
r
残す部分 切り捨てる部分

40. 打ちきり次元（r ）の決め方２：Gavish -Donohoの方法
ノイズの標準偏差 γが既知の場合
1.行列Xが正方行列であるとき、 として、
2.行列Xが矩形行列（正方行列ではない）であるとき、 として、行
列の縦横比 とおくと、
3
4
nτγ=
nn×
∈X
( )
nm
mn
×
<∈X
/mnβ=
()nτ λβ γ=
()()
()
( )
2
1/2
1/2
8
21
1 14 1β
λβ β
β ββ


= ++

++ + +


41. 打ちきり次元（r ）の決め方２：Gavish -Donohoの方法
ノイズの標準偏差 γが未知の場合
3. 行列Xが矩形行列（正方行列ではない）であるとき、 として、行列
の縦横比 とおく。行列Xの特異値の中央値をσ
medianとすると、閾値は
ここで であり、μ
βは以下の方程式で決まる。
( )
nm
mn
×
<∈X
/mnβ=
()
median
τωβσ=
()()/
β
ωβ λβ µ=
( )( )( )( )
()
2
1
1
2
/2
2
11
1
22t
tt
dt
β
µ
β
ββ
π
−

+ − −−


=
∫

42. 打ちきり次元（r ）の決め方２：Gavish -Donohoの方法
ここまで：大変複雑な式でした・・・
しかし・・・
実際のデータ解析の場合には、閾値τ をMatlab関数一つで導くことができます！
tau =
optimal_SVHT_coef(size(X,2)/size(X,1),0)*median(s);
以下、打ち切り次元の二つの方法を実際にやってみましょう。

43. 打ちきり次元（r ）の決め方２：Gavish -Donohoの方法
CH01_SEC07_1_Truncation.m 目的：打ち切り次元の決め方の方法の比較
% クリーンデータ（ランク２）の生成
t = (- 3:.01:3)';
Utrue= [cos(17*t).*exp (-t.^2) sin(11*t)];
Strue= [2 0; 0 .5];
Vtrue= [sin(5*t).*exp( -t.^2) cos(13*t)];
X = Utrue*Strue* Vtrue';
%観測データの生成（ガウスノイズ付加）
sigma = 1;
Xnoisy= X+sigma*randn(size(X));
%特異値分解
[U,S,V] = svd(Xnoisy);
%エネルギーに基づく打ち切り次元の決定
cdS = cumsum(diag(S))./sum(diag(S)); %累積エネルギーの計算
r90 = min(find(cdS>0.90)); %90%のエネルギーを説明する最小の次元を見つける
X90 = U(:,1:r90)*S(1:r90,1:r90)*V(:,1:r90)'; %データの再構成
%Gavish-Donoho法による打ち切り次元の決定
% 簡単のためノイズ分散は既知とする
N = size(Xnoisy,1);
cutoff = (4/sqrt(3))*sqrt(N)*sigma; % 公式
r = max(find(diag(S)>cutoff)); %
Xclean = U(:,1:r)*S(1:r,1:r)*V(:,1:r)'; %データの再構成

44. 打ちきり次元の決め方の比較
Gavish-Donohoの方法 r=2 エネルギーに基づく打ち切り法 r=403
data
X
clean
X

45. 打ちきり次元の決め方の比較：固有顔の場合
% 固有顔の計算
load ../DATA/allFaces.mat
trainingFaces= faces(:,1:sum(nfaces (1:36)));
avgFace= mean(trainingFaces,2); % size n*m by 1;
X = trainingFaces -avgFace*ones(1,size(trainingFaces,2));
[U,S,V] = svd(X,'econ' );
%Gavish- Donohoの方法
s = diag (S);
threshold = ...
optimal_SVHT_coef(size(X,2)/size(X,1),0)*median(s);
r = find(s>threshold,1,’last’ );
% r = 766
%エネルギーの方法
cumS= cumsum(s)/sum(s);
r90 = find(cumS <0.9,1,'last');
% r90 = 1122
Gavish-Donohoの方法
エネルギーの方法
r=1122
r=766

46. ここで演習
以下のコードを走らせて、コードの中身を理解してください。
CH01_SEC07_1_production.m

47. 【本日の内容】特異値分解（Singular value decomposition, SVD）
1.主成分分析、固有顔
-主成分分析と特異値分解の関係
-固有顔による顔分析
2.次元打ち切り（dimensional truncation ）
-打ち切り次元 rの決め方 二通り
3.ランダム化特異値分解（randomized SVD）
-大規模行列の効率的な特異値分解法

48. 大規模行列の特異値分解：ランダマイズSVD （Randomized SVD）
問題：大規模行列の特異値分解には大量の計算コストがかかる。
Goulb- ReinschSVD の計算コスト
（Golub & Van Loan, “Matrix Computation,” Figure 8.6.1, P.493）
考え方： データ量が増えても、データの本質であるランクは増えない。
（4K → 8Kになっても、ランクは増えない）
解法：（１）データを低次元空間にランダムに射影する。
（２）その低次元空間で特異値分解を行う。
（３）低次元空間で得られたベクトル空間を高次元に戻す。
ランダム化線形代数（randomized linear algebra ）という最新分野
( )( )
2 23
max , ,m n mnn

49. 大規模行列の特異値分解：ランダマイズSVD （Randomized SVD）
問題：大規模行列の特異値分解には大量の計算コストがかかる。
SVDにどれくらい時間が掛かるか、数値的に確かめてみよう。
n = 2.^[6:14]; %2^6から2^14まで
cputimes= zeros(size(n));
fork = 1:length(n)
disp(['Processing SVD for n=' num2str(n(k))
'...']);
X = randn(n(k)); % n(k)*n(k)のランダム行列
t0 = cputime(); % 計算開始の時刻
[U,S,V] = svd(X); %特異値分解
cputimes(k) = cputime() -t0; % 計算時間
end
n=2
14
=16384
•n=2
14
=16384では、計算時間はほぼ一時間
•計算時間はn
2.71
でスケール（ n
3
より速いーMatlab はより速いアルゴリズムを使って
いるのかもー）

50. 大規模行列の特異値分解：ランダマイズSVD （Randomized SVD）
データ行列の規則性 ＝低ランク構造
低ランク行列では列が独立ではない ランダムに列をサンプルした行列
アイディア：元々の大きい行列の低ランク構造は、ランダムに列をサンプルした小さ
い行列に含まれているはず
→ ランダムサンプリングで次元を落とした行列で行列分解すればよい

51. ランダマイズSVD （R-SVD）：アルゴリズム
入力：データ行列 およびランク r
ステップ１：ランダム射影行列 を用いて、行列X の列空間を射影したものを Zとする
ステップ２：射影した行列Zを QR分解し、データ行列の張る直交基底を見つける。
ステップ３：ステップ２の直交基底の空間に行列 Xを射影し、行列 Yを得る。
ステップ４：行列Y を特異値分解する。
ステップ５：行列Y の張る空間 U
Yを行列Xの張る空間 Uに引き戻す。
出力：データ行列の特異値分解
nm×
∈X
mr×
∈P
nr×
∈=Z XP 
,,
nr rr××
= ∈∈Z QR Q R
mr×
∈=Y QX 

Y
Y=UΣV

Y
U=QU
X=UΣV


52. ランダマイズSVD （R-SVD）：アルゴリズム
入力：データ行列 およびランク r
出力：データ行列の特異値分解
nm×
∈X
X=UΣV

① ②
③ ④
⑤

53. ランダマイズSVD （R-SVD）：MATLAB コード
function[U,S,V] = rsvd(X,r,q,p)
% ステップ①：データ行列 Xをランダム行列 Pで射影
ny= size(X,2);
P = randn(ny,r+p); % projection matrix ny*(r+p)
Z = X*P;
fork=1:q %power iteration
Z = X*(X'*Z);
end
% ステップ②： 行列Z のQR分解で行列X の列空間を見つける
[Q, R] = qr(Z,0); % QR decomposition of Z
% ステップ③： 行列X をQの張る空間に射影
Y = Q'*X;
%ステップ④： 行列Y を特異値分解
[UY,S,V] = svd(Y, 'econ' );
%ステップ⑤： 行列UYを引き戻す
U = Q*UY;
① ②
③ ④
⑤

54. ランダマイズSVD （R-SVD）による画像圧縮
MATLABコード：CH01_SEC08_RSVD
（注）”Jupiter.jpg”がないので、”dog.jpg”に置き換えています。
clear all, close all , clc
% A=imread('../data/jupiter.jpg');
A=imread('dog.jpg');
X=double(rgb2gray(A));
t0=cputime();
[U,S,V] = svd(X,'econ' ); % Deterministic
SVD
cputime() -t0
% 5.0938 seconds
r = 400; % Target rank
q = 1; % Power iterations
p = 5; % Oversampling parameter
t0 = cputime();
[rU,rS,rV] = rsvd( X,r,q,p); % Randomized SVD
cputime() -t0
% 0.4063 seconds
12倍以上の計算時間差！

55. ランダマイズSVD （R-SVD）による画像圧縮
%% Reconstruction
XSVD = U(:,1:r)*S(1:r,1:r)*V(:,1:r)'; % SVD approx.
errSVD = norm(X- XSVD,2)/norm(X,2);
% 9.1644e- 04 (0.09%)
XrSVD = rU(:,1:r)*rS(1:r,1:r)*rV(:,1:r)'; % rSVD approx.
errrSVD = norm(X-XrSVD,2)/norm(X,2);
% 0.0012 (0.12%)
原画像 SVD再構成画像 R-SVD再構成画像
ほぼ同等の再構成誤差

56. ランダマイズSVD （R-SVD）による画像圧縮
原画像 SVD再構成画像 R-SVD再構成画像
計算時間 21.23秒
残差 3.9×10
-4
計算時間 1.11秒
残差 5.1×10
-4

57. ここで演習
以下のコードを走らせて、コードの中身を理解してください。
CH01_SEC08_RSVD.m

58. 【本日のまとめ】特異値分解（Singular value decomposition, SVD）
1.主成分分析、固有顔
-主成分分析と特異値分解の関係
-固有顔による顔分析
2.次元打ち切り（dimensional truncation ）
-打ち切り次元 rの決め方 二通り
3.ランダム化特異値分解（randomized SVD）
-大規模行列の効率的な特異値分解法