1function [CC] = train_lda_sparse(X,G,par,tol) 2% Linear Discriminant Analysis for the Small Sample Size Problem as described in 3% Algorithm 1 of J. Duintjer Tebbens, P. Schlesinger: 'Improving 4% Implementation of Linear Discriminant Analysis for the High Dimension/Small Sample Size 5% Problem', Computational Statistics and Data Analysis, vol. 52, no. 1, pp. 423-437, 2007. 6% Input: 7% X ...... (sparse) training data matrix 8% G ...... group coding matrix of the training data 9% test ...... (sparse) test data matrix 10% Gtest ...... group coding matrix of the test data 11% par ...... if par = 0 then classification exploits sparsity too 12% tol ...... tolerance to distinguish zero eigenvalues 13% Output: 14% err ...... Wrong classification rate (in %) 15% trafo ...... LDA transformation vectors 16% 17% Reference(s): 18% J. Duintjer Tebbens, P. Schlesinger: 'Improving 19% Implementation of Linear Discriminant Analysis for the High Dimension/Small Sample Size 20% Problem', Computational Statistics and Data Analysis, vol. 52, no. 1, 21% pp. 423-437, 2007. 22% 23% Copyright (C) by J. Duintjer Tebbens, Institute of Computer Science of the Academy of Sciences of the Czech Republic, 24% Pod Vodarenskou vezi 2, 182 07 Praha 8 Liben, 18.July.2006. 25% This work was supported by the Program Information Society under project 26% 1ET400300415. 27% 28% 29% Modified for the use with Matlab6.5 by A. Schloegl, 22.Aug.2006 30% 31% $Id$ 32% This function is part of the NaN-toolbox 33% http://pub.ist.ac.at/~schloegl/matlab/NaN/ 34 35% This program is free software; you can redistribute it and/or 36% modify it under the terms of the GNU General Public License 37% as published by the Free Software Foundation; either version 3 38% of the License, or (at your option) any later version. 39% 40% This program is distributed in the hope that it will be useful, 41% but WITHOUT ANY WARRANTY; without even the implied warranty of 42% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 43% GNU General Public License for more details. 44% 45% You should have received a copy of the GNU General Public License 46% along with this program; if not, write to the Free Software 47% Foundation, Inc., 51 Franklin Street - Fifth Floor, Boston, MA 02110-1301, USA. 48 49%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Step (1) 50%p = length(X(1,:));n = length(X(:,1));g = length(G(1,:)); 51G = sparse(G); 52[n,p]=size(X); 53g = size(G,2); 54 55for j=1:g 56 nj(j) = norm(G(:,j))^2; 57end 58Dtild = spdiags(nj'.^(-1),0,g,g); 59Xtild = X*X'; 60Xtild1 = Xtild*ones(n,1); 61help = ones(n,1)*Xtild1'/n - (ones(1,n)*Xtild'*ones(n,1))/(n^2); 62matrix = Xtild - Xtild1*ones(1,n)/n - help; 63% eliminate non-symmetry of matrix due to rounding error: 64matrix = (matrix+matrix')/2; 65[V0,S] = eig(matrix); 66% [s,I] = sort(diag(S),'descend'); 67[s,I] = sort(-diag(S)); s = -s; 68 69cc = sum(s<tol); 70 71count = n-cc; 72V1 = V0(:,I(1:count)); 73D1inv = diag(s(1:count).^(-1.0)); 74Dhalfinv = diag(s(1:count).^(-0.5)); 75Dhalf = diag(s(1:count).^(0.5)); 76%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Step (2) 77help2 = V1*D1inv; 78M1 = Dtild*G'*Xtild; 79B1 = (G*(M1*(speye(n)-1/n))-help)*help2; 80%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Step (3) 81opts.issym = 1;opts.isreal = 1;opts.disp = 0; 82%if 0, 83try, 84 [V0,S,flag] = eigs(B1'*B1,g-1,'lm',opts); 85 EV = Dhalfinv*V0; 86 [s,I] = sort(-diag(S)); s = -s; 87 %else 88catch 89 % needed as long as eigs is not supported by Octave 90 [V0,S] = eig(B1'*B1); 91 flag = 0; 92 [s,I] = sort(-diag(S)); s = -s(I(1:g-1)); 93 EV = Dhalfinv * V0(:,I(1:g-1)); 94 I = 1:g-1; 95end; 96%EV = Dhalfinv*V0; 97%[s,I] = sort((diag(S)),'descend'); 98%[s,I] = sort(-diag(S)); s = -s; 99if flag ~= 0, 100 'eigs did not converge'; 101end 102%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Step (4) 103for j=1:g-1, 104 C(:,j) = EV(:,I(j))/norm(EV(:,I(j))); 105end 106cc = 0; 107for j=1:g-1, 108 if (1-s(j))<tol 109 cc = cc+1; 110 V2(:,j) = EV(:,I(j)); 111 else 112 break 113 end 114end 115if cc > 0 116 [Q,R] = qr(V2,0); 117 matrix = B1*Dhalf*Q; 118 [V0,S] = eig(matrix'*matrix); 119 %[s,I] = sort(diag(S),'descend'); 120 [s,I] = sort(-diag(S)); s = -s; 121 for j=1:cc 122 C(:,j) = Q*V0(:,I(j)); 123 end 124end 125 126%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Step (5) 127C1 = help2*Dhalf*C; 128trafo(:,1:g-1) = X'*C1 - (X'*ones(n,1))*(ones(1,n)*C1/n); 129for j=1:g-1 130 trafo(:,j) = trafo(:,j)/norm(trafo(:,j)); 131end 132CC.trafo = trafo; 133 134if par == 0 135% X2 = full(test*X'); 136% [pred] = classifs(C1,M1,X2); 137 CC.C1 = C1; 138 CC.M1 = M1; 139 CC.X = X; 140else 141% M = Dtild*G'*X; 142% [pred] = classifs(trafo,M,test); 143 CC.C1 = trafo; 144 CC.M1 = Dtild*G'*X; 145end 146