X = 3×2
8 9
9 6
1 0
特異値分解と主成分分析 (1)
ー MATLAB 版
ベータ版 3.0
早稲田大学教育・総合科学学術院 新井仁之
概要
今回は特異値分解を主成分分析に応用する．前回の「特異値分解のアルゴリズムと次元圧縮」（URLは本ノート末尾参照）の続編である。特異値分解による主成分分析についての数学の解説は
新井仁之，線形代数 基礎と応用，日本評論社，2006
の第17章「多変量解析と線形代数」を元にしている．
本ノートはMATLABのライブエディターを使用して作成した．
必要な用語解説
をデータからなる行列とする．ここで N は個体数（観測対象）の個数であり，p は変量（観測値）の個数である．たとえばある学校（あるいはクラス）の生徒の成績をイメージしてほしい．ここで N は生徒数 p は教科数，
は学生番号 j の生徒の科目 k の点数を表す．変量 k の平均値を 
とおき，
とおく．
とおく．
とおき， とおく． とおく．
を平均偏差行列という．
を分散共分散行列という．
とおくと
となっている．
をデータ 
の分散，
を 
と の共分散という． また，相関係数 
を用いて
を用いて
を相関係数行列という．データ X(:,j) と X(:,k) は相関係数 
が 1 に近いときに相関は強いといい，
のとき相関はないという．
が 1 に近いときに相関は強いといい， のとき相関はないという．
とおき，
を X の標準得点行列という．標準化したものは，各変量の平均が 0 であり，相関係数行列は変わらない．標準化したデータ行列の分散共分散行列は元のデータの標準化係数行列と一致している．
相関係数行列と分散共分散行列には次の関係式がある．
clear
平均偏差行列，分散共分散行列，相関係数行列，標準得点行列の計算
function A = MDmatrix(X)
% 平均偏差行列（mean deviation matrix）
A = X - mean(X,1);
end
function A = VCmatrix(X)
% 分散共分散行列（variance-covariance matrix）
M = MDmatrix(X);
A = (1/size(M,1))*(M')*M;
end
function A = CRmatrix(X)
% 相関係数行列（correleation matrix）
% MATLAB には corrcoef(X) がある．
S = VCmatrix(X);
p = size(S,1);
A = zeros(p,p);
for j =1:p
for k = 1:p
A(j,k) = S(j,k)/(sqrt(S(j,j))*sqrt(S(k,k)));
end
end
end
function D = MD_CR(X)
A = diag(VCmatrix(X));
A = sqrt(A);
A = 1./A;
D = diag(A);
end
function Y = stand(X)
% 標準得点行列
Y = MDmatrix(X)*MD_CR(X);
end
コメント：このうち，相関係数行列に関しては MATLAB の corrcoef があるが，それと上記の CRmatrix は同じであることを確認する．
たとえば次のデータ行列 X を考える．
X = [8 9;9 6;1 0]
disp(char('CRmatrix(X) = ', num2str(CRmatrix(X))));
disp(char('corrcoef(X) = ', num2str(corrcoef(X))));
データの標準化の具体例
データ行列としては前掲の X を考える．
XS = stand(X);
disp(char('X の標準化 = ', num2str(stand(X))));
disp(char('X の標準化の変量の平均 = ', num2str(mean(XS))));
標準化したデータ行列の分散共分散行列は元のデータの標準化係数行列と一致している
disp(char('X の標準化の分散共分散行列 = ', num2str(VCmatrix(XS))));
disp(char('X の相関係数行列 = ', num2str(CRmatrix(X))));
特異値分解による主成分ベクトル
X を 
行列で，N は個体数，p は変量の数とする．もしデータの単位などが不揃いの場合は，標準化したデータを考えていく．X の分散共分散行列（ 標準化したデータの場合の分散共分散行列は元データの相関係数行列）の特異値分解をする．
行列で，N は個体数，p は変量の数とする．もしデータの単位などが不揃いの場合は，標準化したデータを考えていく．X の分散共分散行列（ 標準化したデータの場合の分散共分散行列は元データの相関係数行列）の特異値分解をする．N 次ユニタリー行列 U と p 次ユニタリー行列 V により，その特異値標準形が 
により与えられる．
により与えられる．
の場合を考える．
を
の特異値とする．
，
とする．このとき特異値分解は
で与えられる．
を第 k 主成分（第 k 主成分ベクトル），
を固体 n の第k主成分得点という． 
を第 k 主成分の寄与率，
を第 k 主成分までの累積寄与率という．一般には累積寄与率が 80%を超えるような k までの主成分を考慮することが望ましいとされている．
を第 k 主成分までの累積寄与率という．一般には累積寄与率が 80%を超えるような k までの主成分を考慮することが望ましいとされている．具体例の計算
次の仮想的な成績のデータで試してみよう．
testkekka.xlsx （readtableで読み込んだものを下に表示する）
なおこれは数学的に作った架空の数値データで，実在の学校・個人とは一切関係ない．
データは，50人の架空の人たちの 架空の学生番号，数学，物理，アート，英語の点数からなっている．
filename = 'Score.xlsx'; % ここは自分の好きなデータを使う．
% その場合はこのセクションはデータの形態に合わせて変更
Score = readtable(filename,"TextType","string","VariableNamingRule","preserve")
各科目の平均点
mean(Score(:,["Mathematics" "Physics" "Art" "English"]))
各教科の最低点(Min)，中央値(Median)，最高点(Max)
summary(Score(:,2:end))
個人の平均点を計算して表に付加
Score.Average = round(mean(Score{:,2:5},2),1) %小数点2桁目を四捨五入
平均点のヒストグラム
histogram(Score.Average,10)
主成分分分析
X = table2array(Score(:,2:5));
X_org = X;
X = MDmatrix(X);
% データの標準化をしたいときは % を外す．
% X = stand(X);
% 特異値標準化
[Sigma, U, V]=SingNormal(X); % 前回の特異値分解のレクチャーで作った特異値標準形を使う．
% もちろんMATLAB の svd を使っても良い．
PCAdata = X*V; %主成分得点
主成分ベクトルの表示
for j = 1:size(V,1)
disp(char(['第',num2str(j),'主成分ベクトル'], num2str(V(:,j))));
end
寄与率と累積寄与率の計算
%寄与率と累積寄与率の計算
Kiyo = zeros(length(Sigma),1);
Ruiseki_Kiyo = zeros(length(Sigma),1);
RK = 0;
for j = 1:length(Sigma)
Kiyo(j)=Sigma(j)^2/sum(Sigma.^2);
RK = RK + Kiyo(j);
Ruiseki_Kiyo(j) = RK;
end
for j=1:length(Sigma)
disp(['第',num2str(j),'主成分寄与率 =',num2str(Kiyo(j))]);
end
for j=1:length(Sigma)
disp(['第',num2str(j),'累積寄与率 =',num2str(Ruiseki_Kiyo(j))]);
end
平均値による順位(order_avg)，第一主成分得点(fristPC)，第二主成分得点(secondPC)，第三主成分得点(thirdPC)を表に付加．
Score_avg = sortrows(Score,6,'descend');
a = table2array(Score_avg(:,1));
[~,p] = sort(a,'ascend');
Score.order_avg = p; %平均点による順位
Score.firstPC = PCAdata(:,1); %第1主成分得点
Score.secondPC = PCAdata(:,2); %第2主成分得点
Score.thirdPC = PCAdata(:,3); %第3主成分得点
Score_sub = Score(:,["Student_ID" "Average" "order_avg" "firstPC" "secondPC" "thirdPC"])
補足
以下は前回のレクチャー「特異値分解のアルゴリズムと次元圧縮」より必要なコマンド
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 = 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
特異値標準形を与えるアルゴリズム
以上の準備のもとに特異値標準系を与えるアルゴリズムを記す．このアルゴリズム自体は線形代数でよく知られたものである．（たとえば [1] 参照．）
function [Sigma, U, W]=SingNormal(A)
% 行列 A を与える．Aは MxN 行列
% Sigma は A の特異値 min(M,N)x1 配列
% U,W は U'*A*W が特異値標準形を与えるユニタリー行列
[S1,S2]=size(A);
% A_org = A;
if S1>S2
A=A';
end
[N1,N2]=size(A);
B=A'*A;
L = min(N1,N2);
%N = size(B,1);
[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
参考文献
[1] 新井仁之，線形代数 基礎と応用，日本評論社，2006. Kindle版（2022年版）があります．
履歴
ベータ版 1.1 2024/7/19 公開
ベータ版 2.0 2024/7/20 不十分な箇所をいくつか修正及び補足の追加
ベータ版 3.0 2024/8/2 minor changes