特異値分解のアルゴリズムと次元圧縮
ー MATLAB 版
早稲田大学教育・総合科学学術院 新井仁之
今回は特異値分解を，MATLABを使いながら学んでいく．このノートではエルミート行列の対角化は既習として，特異値標準化・特異値分解の証明（= アルゴリズム）の一つをプログラミングしながら学んでいく．特異値分解の応用として画像データの次元圧縮についても述べる．
１．準備
基本的な記法を設定しておく．
ℝ を実数全体のなす集合，ℂ を複素数全体からなる集合とし，𝕂 は ℝ または ℂ を表す．𝕂 はスカラー体，𝕂 の要素はスカラーと呼ばれる．
N次元縦ベクトル 
（ただし 
）全体からなる集合を 
で表す．これにはスカラー積と和が定義され，N次元線形空間をなしている（説明は既知として省くが，たとえば [1] 参照）．
（ただし ）全体からなる集合を で表す．これにはスカラー積と和が定義され，N次元線形空間をなしている（説明は既知として省くが，たとえば [1] 参照）． には次のエルミート内積（あるいはユークリッド内積）が定義される．
（ただし 
は 
の複素共役）．また 
を 
のユークリッド・ノルムあるいは長さという．以下，特に混乱のない場合は
を のユークリッド・ノルムあるいは長さという．以下，特に混乱のない場合は
，
と略記する．
M行N列の行列とは
（ただし 
） である．このような行列全体のなす集合を本ノートでは
で表すことにする．各成分の複素共役をとり転置したものを 
で表し．A の共役行列という．
で表し．A の共役行列という．
である． は N行M列の行列になっている．
は N行M列の行列になっている．MATLABでは A' で実行できる．たとえば
A = randi([1 10],2,3)+1i*randi([1 10],2,3);
disp(char('A =', num2str(A)));
disp(char('A^* =', num2str(A')));
２．エルミート行列の対角化
特異値分解で使うので，エルミート行列の対角化について概説する．
A が N行N列の行列（N次正方行列）で，特に
をみたすものをエルミート行列という．また N次正方行列で
（ただし 
は単位行列）をみたすものをユニタリー行列という．
N 次正方行列 A を考える．複素数 λ に対する方程式 
が非自明解を持つとする．すなわち，
が非自明解を持つとする．すなわち，ある 
で， をみたすものが存在するとする．このとき λ を A の固有値， を λ に属する固有ベクトルという．λ に属する k 個の線形独立な固有ベクトルが存在するとき，λ の重複度は k であるという．このとき，この固有値の重複度を込めた個数は k 個であるともいう．
の場合，代数学の基本定理から，一般のN 次正方行列 A は重複度も込めて N 個の固有値を持つことが知られている（[1]参照）．
で， をみたすものが存在するとする．このとき λ を A の固有値， を λ に属する固有ベクトルという．λ に属する k 個の線形独立な固有ベクトルが存在するとき，λ の重複度は k であるという．このとき，この固有値の重複度を込めた個数は k 個であるともいう． の場合，代数学の基本定理から，一般のN 次正方行列 A は重複度も込めて N 個の固有値を持つことが知られている（[1]参照）．また，エルミート行列の固有値はすべて実数であることも知られている（[1]参照）．
特異値分解では次の定理を使っていく．
定理1 （エルミート行列の対角化） A を N 次正方行列で，エルミート行列であるとする．A の重複度も込めた N 個の固有値を 
とする．このとき，ある N 次 ユニタリー行列 U を
とする．このとき，ある N 次 ユニタリー行列 U を
をみたすようにとれる．
本ノートではエルミート行列の対角化の証明はしない．これについては通常の大学１年の線形代数のコースで学ぶ．たとえば [1] に詳しい証明が記されているので，興味のある読者は参照してほしい．
また，MATLAB では次のように eig コマンドを使ってエルミート行列の対角化ができるので，このコマンドを使っていくことにする．
eig を使って具体例でエルミート行列の対角化を行う．
まずエルミート行列をランダムに作る．
N = 3;
B = randi([0 9],N,N)+1i*randi([0 9],N,N);
A = B+B'; %これは作り方から明らかにエルミート行列
disp(char('A =', num2str(A)));
[U,D] = eig(A);
disp(char('ユニタリー行列 U =', num2str(U)));
disp(char('対角行列 D =', num2str(D)));
検算：
まず U がユニタリー行列であることの確認．
以下，行列 に対して，
とする．
に対して， とする．C1 = sum(abs(U*U'-eye(3)),'all');
disp(char('|UU*-I|=', num2str(C1)));
次に対角化の確認
C2 = sum(abs(U'*A*U-D),'all');
disp(char('|U*AU-D|=', num2str(C1)));
３．特異値標準化
とする．これに対する特異値標準化について説明していく．特異値標準化のアルゴリズム
とし，
とする．明らかに 
である．
は N 次正方行列になっている．さらに 
よりエルミート行列である． の固有値は実数であるが，さらにこのタイプの行列の固有値は非負であることが示せる．そこでその固有値を
と大きい順に並べる．
（
）とおき，これを A の特異値という．
であり，
の場合は 
である．さて， は N 次のエルミート行列であるから，定理1より N 次ユニタリー行列 V が存在し，
をみたすものがとれる．
とおき，
は N 次のエルミート行列であるから，定理1より N 次ユニタリー行列 V が存在し， をみたすものがとれる． とおき，
（
）と定める．このとき 
が 
の正規直交系であることが証明できる．
であるが，特に 
の場合には，適当に 
を追加して，
が の正規直交基底になるようにできる．
が の正規直交系であることが証明できる． であるが，特に の場合には，適当に を追加して， が の正規直交基底になるようにできる．具体的にどのようなベクトルを追加すればよいか？
そのプログラムは次のように書ける．
Lemma. 足りない正規直交系を補完して正規直交基底にするアルゴリズム
function A = HokanONB(Q)
% Q は Nxp 配列，ただし N>=p （行数 >= 列数）である必要がある.
% Q(:,1),...,Q(:,p) は正規直交系になるようにしておく．
% A は NxN 配列で，A(:,1),...,A(:,N)は正規直交基底であり，かつ A(:,k)=Q(:,k), 1<=k<=p
% Q は正規直交系でない場合，Q の部分も合わせて直交化される．
N=size(Q,1);
n=size(Q,2);
if N<n
disp('入力配列の列数が行数を超えてます．')
return
end
if n == N
A = GramSchmidt(Q);
return
end
A=zeros(N,N);
A(:,1:n) = Q;
NB = eye(N);
%m=0;
M=0;
for k = 1:N
l = k+M;
A(:,n+l)=NB(:,k);
rank_A = rank(A);
if rank_A < n+1
A(:,n+l)=zeros(N,1);
M=M-1;
end
if rank_A == N
A=GramSchmidt(A);
return
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q = GramSchmidt(A)
% グラム・シュミットの正規直交化
% 正規直交基底を保存するための行列
Q = zeros(size(A));
n = size(A, 2);
for i = 1:n
v = A(:, i);
for j = 1:i-1
v = v - (Q(:, j)' * A(:, i)) * Q(:, j);
end
% 正規化
Q(:, i) = v / norm(v);
end
end
試用例
で試してみる．まず二つの線形独立なベクトルを適当に作るv1 = [1 2 3 4]; v2 = [2 3 4 5];
グラム-シュミットでこれを正規直交系にする．
A = zeros(4,2);
A(:,1)=v1';
A(:,2)=v2';
Q = GramSchmidt(A);
disp(char('与えた正規直交系', num2str(Q)));
Q の正規直交性を確認
Q'*Q
Q を補完して， の正規直交基底を作る．
の正規直交基底を作る．X = HokanONB(Q);
disp(char('Xの最初の二つは与えられた正規直交系と一致してることの確認', num2str(X(:,1:2))));
X の正規直交性を確認
X'*X
disp(char('よって求める残りのベクトルは', num2str(X(:,3:4))));
特異値標準形を与えるアルゴリズム
以上の Lemma を使って特異値標準系を与えるアルゴリズムを記す．このアルゴリズム自体は線形代数でよく知られたものである．（たとえば [1] 参照．）
function [Sigma, U, W]=SingNormal(A)
% Version 1.1 (July 6, 2024) by Hitoshi Arai
% 行列 A を与える．Aは MxN 行列
% Sigma は A の特異値を並べた min(M,N)x1 配列
% U,W は U'*A*W が特異値標準形を与える MxM 及び NxN ユニタリー行列
% エルミート行列の対角化の部分でMATLABの eig は使用している．
[S1,S2]=size(A);
if S1>S2
A=A';
end
[N1,N2]=size(A);
B=A'*A;
L = min(N1,N2);
[V,D]=eig(B);
r = rank(B);
Sigma = sqrt(diag(D));
[Sigma,p] = sort(Sigma,'descend');
Sigma = Sigma(1:L);
VV = zeros(size(V));
for k = 1:size(V,1)
VV(:,k) = V(:,p(k));
end
U1 = zeros(N1,r);
for k=1:r
U1(:,k)=Sigma(k)^(-1)*A*VV(:,k);
end
MinSize = min(N1,N2);
if r == MinSize
U=U1;
W=VV;
else
U = HokanONB(U1);
W=VV;
end
if S1>S2
U_org = U;
U=W;
W=U_org;
end
end
特異値分解の検算
まず適当に 4x5行列を与える．ただし rank≤3．
A = randi([0 9],4,3);
A(:,4)=A(:,3);
A(:,5)=A(:,2)
rank(A)
ans = 3
[Sigma,U,V] = SingNormal(A)
Sigma %特異値
U'*A*V %特異値標準形
5x4 行列の場合の特異値分解の検証，
B=A' % 5x4行列
[Sigma2,U2,V2] = SingNormal(B)
U2'*B*V2
MATLABの svd との比較
MATLABに装備されている svd との比較をしておく．なお，前述の SingNormal は特異値標準化定理の証明（[1]）を再現した（特異値標準化定理の証明を学ぶための）プログラムで，高速化は考慮していないことに注意．
[U_matlab,S_matlab,V_matlab]=svd(A)
sum(abs(U_matlab'*A*V_matlab - U'*A*V),'all')
ans = 4.2387e-14
計算時間の比較：1000x2000行列でテスト．
% TestMatrix = rand(1000,2000);
% save('TestMatrix','TestMatrix');
load TestMatrix.mat
tic;
[S,U,V] = SingNormal(TestMatrix);
toc
tic;
[U_matlab,S_matlab,V_matlab]=svd(TestMatrix);
toc
コメント：SingNormal は特異値標準化定理の証明を理解するためのもので，高速化していないから，MATLABの svd の方が断然早い．実際に使う場合は MATLAB の svd がよい．
４．画像の特異値分解と次元圧縮
本節では，特異値分解を使った画像の次元圧縮について述べる．まず特異値分解から始める．
行列 
の特異値標準化 
を考える． とおく．U は M次ユニタリー行列であるから，
の特異値標準化 を考える． とおく．U は M次ユニタリー行列であるから，
は 
の正規直交基底，また
の正規直交基底，また
は 
の正規直交基底になっている．また 
とし，
を正の特異値とする．このとき
の正規直交基底になっている．また とし， を正の特異値とする．このとき
が成り立つ．これを A の特異値分解という．
特異値分解を用いた次元圧縮は次のようなものである．
に対して
に対して
このようにすると，
である．これをプログラムする．
である．これをプログラムする．function A_k = SingDecomp(A,k)
[s1,s2]=size(A);
A_k = zeros(s1,s2);
[S,U,V] = SingNormal(A);
for j=1:k
A_k = A_k + S(j)*U(:,j)*conj(V(:,j)');
end
end
画像の読み込み
A = imread('RIMG0030_AI_Gray.tif');
A= double(A);
検算 r=rank(A) のとき，A=A_r であることの確認
r = rank(A);
disp(['原画像 A の階数は',num2str(r),'です']);
A_r = SingDecomp(A,r);
sum(abs(A - A_r),'all')
ans = 6.6423e-05
特に 
と原画像 A との比較をする．
と原画像 A との比較をする．a = [5 10 120];
[s1,s2]=size(A);
AK = zeros(s1,s2,length(a));
j=1;
for k = 1:a(end)
if ismember(k,a)
AK(:,:,j) = SingDecomp(A,k);
j=j+1;
disp(['近似データ A_{',num2str(k),'} の階数は',num2str(k),' です。'])
end
end
figure
subplot(2,2,1)
imshow(uint8(A))
title('原画像')
for j = 1:length(a)
subplot(2,2,j+1)
imshow(AK(:,:,j),[])
title(['A_{',num2str(a(j)),'}'])
end
%print('CompressDim', '-dtiff', '-r300'); % 画像ファイルを保存する場合
参考文献
[1] 新井仁之，線形代数 基礎と応用，日本評論社，2006. Kindle版（2022年版）があります．
