cvLabMATLAB/orzPCA.m at main · ComputerVisionLaboratory/cvLabMATLAB · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
function [Z,U,D,E] = orzPCA(X,Y,varargin)
%function [Z,U,D,E] = orzPCA(X,Y,varargin)
% 主成分分析 ver.1.00 by ohkawa
% input
%  X: 入力マトリックス
%  Y: 部分空間の次元もしくは寄与率
%  Yが1.0以下ならば， Yは寄与率と判断
%  Yが1.0超過ならば， Yは部分空間の次元と判断
%  第三引数：デフォルトでは，共分散行列のPCA
%  第三引数が'R'の時，自己相関行列のPCA
%
%　nDim< nNumの時は通常のPCA
%  nDim>=nNumのときは線形カーネルPCA（双対問題）を適用する
% output
%  Z: 主成分空間に射影されたデータ
%  U: 主成分ベクトル
%  D: 固有値
%  E: 寄与率


X = X(:,:);
[nDim,nNum] =size(X);

flgM = true;
if nargin == 3
    if varargin{1} == 'R'
        flgM = false;
    end
end

if Y <= 1 % 寄与率
    cRate = Y;
    if nDim<nNum
        if flgM==true
            C = cov(X',1);
        else
            C = X*X'/nNum;
        end
        [U,tmpD]= eig(C);
        [D ind]= sort(diag(tmpD),'descend');
        U = U(:,ind);
        nSubDim = find(cumsum(D)/sum(D)>=cRate, 1 );
        U=U(:,1:nSubDim);
        D = D(1:nSubDim);
        E = sum(D)/trace(C);
    else
        if flgM==true
            K = X'*X;
            IN =  ones(nNum,nNum)/nNum;
            K = K - IN*K - K*IN + IN*K*IN;
        else
            K = X'*X;
        end
        [A B] = eig(K);
        [B ind] = sort(diag(B),'descend');
        A=A(:,ind);
        D = B/nNum;
        nSubDim = find(cumsum(D)/sum(D)>=cRate, 1 );
        A=A(:,1:nSubDim);
        B=B(1:nSubDim);
        A = A/sqrt(diag(B));
        U = X*A;
        D = D(1:nSubDim);
        E = sum(D);
    end
elseif Y > 1 % 主成分空間の次元
    nSubDim = floor(Y);
    if nDim<nNum
        if flgM==true % 共分散行列
            C = cov(X',1);
        else % 自己相関行列
            C = X*X'/nNum;
        end

        OPTS.disp = 0;
        if nSubDim<nDim
           [U tmpD] = eigs(C,nSubDim,'lm',OPTS); %[U,tmpD]= eigs(C,nSubDim);
        else
           [U tmpD] = eig(C);
        end
        [D ind]= sort(diag(tmpD),'descend');
        U = U(:,ind);
        E = sum(D)/trace(C);
    else %カーネル
        if flgM==true % 共分散行列
            K = X'*X;
            IN =  ones(nNum,nNum)/nNum;
            K = K - IN*K - K*IN + IN*K*IN;
        else  % 自己相関行列
            K = X'*X;
        end
        OPTS.disp = 0;
        [A B] = eigs(K,nSubDim,'lm',OPTS);
        [B ind] = sort(diag(B),'descend');
        A=A(:,ind);
        D = B/nNum;
        A = A/sqrt(diag(B));
        U = X*A;
        E = sum(D);
    end
end
Z = U'*X;