Peter M. Groen / openbr

Commit b860556d56a97ba0c8d330e77cfcc603b9efd673

Authored by Brendan K
2 parents 2a45c457 01ccd1bb

Merge pull request #394 from biometrics/matlabLDA 

Matlab lda
Inline Side-by-side
Showing 3 changed files with 263 additions and 0 deletions
scripts/matlab/LDA.m 0 → 100644
  View file @b860556
  1 +function [subspaceData]=LDA(X,classNo,varargin)
  2 +% [subspaceData]=LDA(X,classNo)
  3 +%
  4 +% LDA method for learning a subspace that seeks to maximize the Fisher
  5 +% seperability measure. 'X' is d x n matrix, where d is the feature
  6 +% vector size, and n is the number of instaces. 'classNo' is a n x 1
  7 +% vector that indicates which class/subject each instance belongs to.
  8 +%
  9 +% Optional parameters:
  10 +
  11 +% 'EnergyRetain' - percentage of variance (0.0 to 1.0) to retain in the initial
  12 +% PCA step. (default = 0.98)
  13 +%
  14 +% 'do_direct' - Whether or not to perform Direct LDA (default = false)
  15 +%
  16 +% 'FixedRetain' - Keep a fixed number of eigenvectors in the initial PCA.
  17 +% If used, then specify number to keep (e.g. 100)
  18 +%
  19 +%
  20 +% 'DominantEig' - whether or not to use the dominant eigenvector method. If
  21 +% the d >> n, then this should be set to true.
  22 +% (default = true)
  23 +%
  24 +% 'SkipPCA' - whether of not to skip the PCA step. If Sw is believed to
  25 +% be non-singular then PCA step can be safely skipped.
  26 +% (default = false)
  27 +%
  28 +% 'ScaleW' - the factor by which to scale the within-class scatter matrix.
  29 +% This controls the importance of the between and within
  30 +% class scatter to each other in the fisher criterion
  31 +%
  32 +% The output 'subspaceData' is a struct that contains the following important fields:
  33 +% subspaceData.mean - d x 1 vector containing the mean of the training data
  34 +% subspaceData.W - d x k LDA projection matrix, where k is the number of
  35 +% subspace dimensions
  36 +%
  37 +%
  38 +% Algorithm based on Fukanaga's LDA method. This code was orginally written by Zhifeng Li,
  39 +% and has since been modified and improved by Brendan Klare.
  40 +%
  41 +%
  42 +
  43 + useFixedEnergy = false;
  44 + doDominantEig = true;
  45 + doPCA = true;
  46 + useFixedEig = false;
  47 + energyPercentage = .98;
  48 + ScaleW = 1;
  49 + doDLDA = false;
  50 + do_null = false;
  51 +
  52 + cnt = 1;
  53 + while cnt < length(varargin)
  54 + switch varargin{cnt}
  55 + case 'EnergyRetain'
  56 + useFixedEnergy = true;
  57 + energyPercentage = varargin{cnt+1};
  58 + case 'FixedRetain'
  59 + useFixedEig = true;
  60 + fixedEig = varargin{cnt+1};
  61 + case 'DominantEig'
  62 + doDominantEig = varargin{cnt+1};
  63 + case 'SkipPCA'
  64 + doPCA = ~varargin{cnt+1};
  65 + case 'ScaleW'
  66 + ScaleW = varargin{cnt+1};
  67 + case 'do_direct'
  68 + do_null = varargin{cnt+1};
  69 + otherwise
  70 + fprintf('Error, unknown argument %s\n',varargin{cnt});
  71 + return
  72 + end
  73 + cnt = cnt + 2;
  74 + end
  75 +
  76 +
  77 + [FeatureLength SampleNumber]=size(X);
  78 + % ClassNum=round(SampleNumber/2);
  79 +
  80 + [a1 a2 classNo] = unique(classNo);
  81 + ClassNum = max(classNo);
  82 +
  83 + %Calculate eigenspace from X
  84 + if doPCA
  85 + if ~doDominantEig
  86 + [eigenvalues, eigenvectors, Mean_Vector]=PCA2(X);
  87 + else
  88 + [eigenvalues, eigenvectors, Mean_Vector, V1]=PCA(X);
  89 + end
  90 +
  91 + if useFixedEnergy
  92 + d1 = cumsum(eigenvalues)./sum(eigenvalues);
  93 + [a nEigs] = max(d1 > energyPercentage);
  94 + elseif useFixedEig
  95 + nEigs = fixedEig;
  96 + else
  97 + nEigs = min(FeatureLength,SampleNumber - 1);
  98 + end
  99 +
  100 + %Select eigenvectors
  101 + Select_eigenvectors=eigenvectors(:,1:nEigs);
  102 + eigenvalues = eigenvalues(1:nEigs);
  103 +
  104 + %Project the sample data on to the eigenvectors
  105 + W=Select_eigenvectors'*(X-repmat(Mean_Vector,1,SampleNumber));
  106 + else
  107 + W = X;
  108 + Mean_Vector = zeros(size(X,1),1);
  109 + nEigs = size(X,1);
  110 + Select_eigenvectors = eye(size(X,1));
  111 + eigenvalues = ones(size(X,1),1);
  112 + end
  113 +
  114 + %Caculate the centers for each class
  115 + ClassCenters = zeros(nEigs,ClassNum);
  116 + for i = 1:ClassNum
  117 + ClassCenters(:,i) = mean(W(:,classNo == i),2);
  118 + end
  119 +
  120 + for i = 1:ClassNum,
  121 + W(:,classNo == i) = W(:,classNo == i) - repmat(ClassCenters(:,i),1,sum(classNo == i));
  122 + end
  123 +
  124 + if ScaleW ~= 1
  125 + W = ScaleW .* W;
  126 + end
  127 +
  128 + [W_val, W_vec, W_m]=PCA2(W);
  129 +
  130 + if ~do_null
  131 + nDim2 = min(nEigs,SampleNumber - ClassNum);
  132 + SW_val=W_val(1:nDim2);
  133 + SW_vec=W_vec(:,1:nDim2);
  134 + SW_vec=SW_vec./(repmat(SW_val',[size(SW_vec,1) 1]).^0.5);
  135 + else
  136 + nDim2 = nEigs;
  137 + SW_val = W_val;
  138 + SW_vec = W_vec;
  139 + if nEigs > SampleNumber - ClassNum
  140 + SW_val(SampleNumber-ClassNum+1:end) = SW_val(SampleNumber-ClassNum)/2;
  141 + end
  142 +
  143 + d1 = cumsum(W_val)/sum(W_val);
  144 + [d1 start_idx] = max(d1 > .1);
  145 +
  146 + SW_vec = SW_vec(:,start_idx:end);
  147 + SW_val = SW_val(start_idx:end);
  148 + nDim2 = size(SW_vec,2);
  149 + SW_vec=SW_vec./(repmat(SW_val',[size(SW_vec,1) 1]).^0.15);
  150 + end
  151 +
  152 +
  153 + m = mean(W,2);
  154 + M=repmat(m(:),1,ClassNum);
  155 + mean2 = m;
  156 + B=SW_vec'*(ClassCenters-M);
  157 + % Between_Class_Matrix=B*B';
  158 +
  159 +
  160 + [B_val,B_vec,B_m]=PCA2(B);
  161 +
  162 + nDim3 = min(ClassNum-1,nDim2);
  163 + SB_vec=B_vec(:,1:nDim3);
  164 +
  165 + subspaceData.mean = Mean_Vector(:);
  166 + subspaceData.mean2 = mean2(:);
  167 + subspaceData.W1 = Select_eigenvectors;
  168 + subspaceData.D1 = eigenvalues;
  169 + subspaceData.W2 = SW_vec;
  170 + subspaceData.W3 = SB_vec;
  171 + subspaceData.W = (subspaceData.W3' * subspaceData.W2' * subspaceData.W1')';
scripts/matlab/PCA.m 0 → 100644
  View file @b860556
  1 +function [eigenvalues, eigenvectors, meanVector, V]=PCA(X,varargin)
  2 +%function [eigenvalues, eigenvectors, meanVector, V]=PCA(X)
  3 + cnt = 1;
  4 + doVar = false;
  5 + doEigs = false;
  6 +
  7 + if nargin < 2,
  8 + end
  9 + while cnt < length(varargin)
  10 + switch varargin{cnt}
  11 + case 'VarEnergy'
  12 + doVar = true;
  13 + varPercent = varargin{cnt+1};
  14 + case 'nEigs'
  15 + doEigs = true;
  16 + eigKeep = varargin{cnt+1};
  17 + otherwise
  18 + fprintf('Error, unknown argument %s\n',varargin{cnt});
  19 + return
  20 + end
  21 + cnt = cnt + 2;
  22 + end
  23 +
  24 +
  25 + [Row Column]=size(X);
  26 +
  27 + %Mean center X
  28 + meanVector = mean(X,2); meanVector = meanVector(:);
  29 + M=repmat(meanVector,1,Column);
  30 + X=X-M;
  31 +
  32 + C=X'*X./Column;
  33 + [V,D]=eig(C);
  34 + eigenvalues=diag(D);
  35 +
  36 + %Ordered by eigenvalues%
  37 + [eigenvalues,Index]=sort(eigenvalues,'descend');
  38 +
  39 + V=V(:,Index) ; %V1 is the the eigenvectors got from X'X;
  40 + eigenvectors=X*V;%eigenvectors is the eigenvectors for XX';
  41 +
  42 + %normalize%
  43 + NV=sum(eigenvectors.^2);
  44 + NV=NV.^(1/2);
  45 +
  46 + %normalize eigenvectors
  47 + NM=repmat(NV,Row,1);
  48 + eigenvectors=eigenvectors./NM;
  49 +
  50 + if doVar
  51 + d = cumsum(eigenvalues)/ sum(eigenvalues);
  52 + [a1 a2] = max(d > varPercent);
  53 + eigenvalues = eigenvalues(1:a2);
  54 + eigenvectors = eigenvectors(:,1:a2);
  55 + end
  56 +
  57 + if doEigs
  58 + eigenvalues = eigenvalues(1:eigKeep);
  59 + eigenvectors = eigenvectors(:,1:eigKeep);
  60 + end
  61 +
  62 + %normalize V1;
  63 + NN=repmat(NV,[Column,1]);
  64 + V=V./NN;
  65 +
  66 +
  67 +
  68 +end
scripts/matlab/PCA2.m 0 → 100644
  View file @b860556
  1 +function [eigenvalues, eigenvectors, Mean_Vector]=PCA2(X,varRetain)
  2 +% [eigenvalues, eigenvectors, Mean_Vector]=PCA2(X)
  3 +%
  4 +%Compute the eienvectors of X when the sample number is larger than the feature lenght
  5 +
  6 +[Row Column]=size(X);
  7 +Mean_Vector=mean(X,2);
  8 +m=repmat(Mean_Vector(:),1,Column);
  9 +X=X-m;
  10 +C=X*X'./Column;
  11 +[V,D]=eig(C);
  12 +eigenvalues=diag(D);
  13 +%Ordered by eigenvalues%
  14 +[eigenvalues,Index]=sort(eigenvalues);
  15 +eigenvalues=eigenvalues(end:-1:1);
  16 +Index=Index(end:-1:1);
  17 +eigenvectors=V(:,Index);
  18 +
  19 +if nargin == 2,
  20 + d = cumsum(eigenvalues)/ sum(eigenvalues);
  21 + [a1 a2] = max(d > varRetain);
  22 + eigenvalues = eigenvalues(1:a2);
  23 + eigenvectors = eigenvectors(:,1:a2);
  24 +end