V_LDATRACE Calculates an LDA transform to maximize trace discriminant [a,f,B,W]=(b,w,n,c) If a feature vector X can come from one of several class and W and B are respectively the within-class and between-class covariance matrices, then the generalized Fisher discriminant F=trace(W\B) is a measure of how well the feature vector discriminates between the classes. If we choose a rectangular (tall, skinny) transformation matrix, we can define a smaller feature vector Y=A'*X. The aim of this routine is to choose A to maximize the Fisher discriminant. We assume that W is positive definite and B is positive semi-definite. The input argument C allows the uset to pre-specify some of the columns of A. Inputs: w[m,m] = within class covariance matrix of x b[m,m] = between class covariance matrix of x [default = I] n is the number of columns in output matrix A [default = M] c[m,r] specifies the first few columns of A to be predefined values [default = null) Outputs: a[m,n] is the transformation matrix: y=a'*x f[n,1] gives the incremental gain in f value for successive columns of A B(n,n) gives the between-class covariance matrix of y W[n,n] gives the within-class covariance matrix of y
0001 function [a,f,B,W]=v_ldatrace(b,w,n,c) 0002 %V_LDATRACE Calculates an LDA transform to maximize trace discriminant [a,f,B,W]=(b,w,n,c) 0003 % If a feature vector X can come from one of several class and W and B are respectively 0004 % the within-class and between-class covariance matrices, then the generalized Fisher discriminant 0005 % F=trace(W\B) is a measure of how well the feature vector discriminates between the classes. 0006 % If we choose a rectangular (tall, skinny) transformation matrix, we can define a smaller 0007 % feature vector Y=A'*X. The aim of this routine is to choose A to maximize the Fisher 0008 % discriminant. We assume that W is positive definite and B is positive semi-definite. 0009 % The input argument C allows the uset to pre-specify some of the columns of A. 0010 % 0011 % Inputs: 0012 % w[m,m] = within class covariance matrix of x 0013 % b[m,m] = between class covariance matrix of x [default = I] 0014 % n is the number of columns in output matrix A [default = M] 0015 % c[m,r] specifies the first few columns of A to be predefined values [default = null) 0016 % 0017 % Outputs: 0018 % a[m,n] is the transformation matrix: y=a'*x 0019 % f[n,1] gives the incremental gain in f value for successive columns of A 0020 % B(n,n) gives the between-class covariance matrix of y 0021 % W[n,n] gives the within-class covariance matrix of y 0022 0023 % Copyright (C) Mike Brookes 1997 0024 % Version: $Id: v_ldatrace.m 10865 2018-09-21 17:22:45Z dmb $ 0025 % 0026 % VOICEBOX is a MATLAB toolbox for speech processing. 0027 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0028 % 0029 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0030 % This program is free software; you can redistribute it and/or modify 0031 % it under the terms of the GNU General Public License as published by 0032 % the Free Software Foundation; either version 2 of the License, or 0033 % (at your option) any later version. 0034 % 0035 % This program is distributed in the hope that it will be useful, 0036 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0037 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0038 % GNU General Public License for more details. 0039 % 0040 % You can obtain a copy of the GNU General Public License from 0041 % http://www.gnu.org/copyleft/gpl.html or by writing to 0042 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0043 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0044 0045 m=size(b,1); % dimension of data vectors 0046 if nargin<4 0047 r=0; 0048 if nargin<3 0049 n=m; 0050 if nargin<2 0051 w=eye(m); 0052 end 0053 end 0054 else 0055 r=size(c,2); % number of columns that are pre-specified 0056 end 0057 if r 0058 if n>r % need to find additional vectors 0059 g=chol(w); 0060 v=g\null(c'*g'); 0061 [p,l,q]=svd(v'*b*v); 0062 a(:,r+1:n)=v*p(:,1:n-r); 0063 a(:,1:r)=c; 0064 else 0065 a=c; % no new vectors to find 0066 end 0067 if nargout>1 0068 ari=a/triu(qr(chol(a'*w*a))); % matrix a must be of full rank 0069 f=diag(ari'*b*ari); 0070 end 0071 else 0072 [g,d]=eig(b,w,'qz'); 0073 [ds,is]=sort(-diag(d)); 0074 a=g(:,is(1:n)); 0075 if nargout>1 0076 f=-ds(1:n); 0077 end 0078 end 0079 if nargout > 2 0080 B=a'*b*a; 0081 if nargout > 3 0082 W=a'*w*a; 0083 end 0084 0085 end