V_DUALDIAG Simultaneous diagonalisation of two hermitian matrices [A,D,E]=(W,B) Inputs: W,B Two square hermitian matrices Outputs: A Diagonalizing matrix (not normally unitary or hermitian) D Real diagonal matrix elements: A'*B*A=diag(D) (see note below) E Real diagonal matrix elements: A'*W*A=diag(E) Note: At least one of W and B must be either positive definite or negative definite. If this is not the case, then D=A'*B*A and E=A'*W*A will be complex hermitian matrices rather than being vectors of real diagonal elements. The columns of A will be ordered so that abs(D./E) is in descending order. If two output arguments are given then A will be scaled to make diag(E)=I but, if W is singular, this will cause some elements of A to be infinite. If A'*B*A=diag(D) and A'*W*A=diag(E) then A'*W*A*diag(1./E)=I so A'*B*A=A'*W*A*diag(D./E) and hence B*A=W*A*diag(D./E) so the columns of A are the eigenvectors of W\A or equivalently the generalized eigenvectors of (B,W). Suppose we have several N-dimensional data row-vectors arising from each of C different classes of data. for each class, c, we can form the mean data vector m(c) and the within-class covariance matrix W(c) We can then form the between class covariance matrix B by taking the covariance of the mean vectors m(1), m(2), ... and also the averaged within-class covariance matrix W by averaging W(1), W(2), ... If we then take A=v_dualdiag(W,B) and postmultiply all our original data vectors by A, we obtain new data vectors for which the average within-class covariance matrix is the identity and for which the first few components contain most of the information that is useful in discriminating between classes.
0001 function [a,d,e]=v_dualdiag(w,b) 0002 %V_DUALDIAG Simultaneous diagonalisation of two hermitian matrices [A,D,E]=(W,B) 0003 % Inputs: W,B Two square hermitian matrices 0004 % 0005 % Outputs: A Diagonalizing matrix (not normally unitary or hermitian) 0006 % D Real diagonal matrix elements: A'*B*A=diag(D) (see note below) 0007 % E Real diagonal matrix elements: A'*W*A=diag(E) 0008 % 0009 % Note: At least one of W and B must be either positive definite or negative 0010 % definite. If this is not the case, then D=A'*B*A and E=A'*W*A will be 0011 % complex hermitian matrices rather than being vectors of real diagonal elements. 0012 % 0013 % The columns of A will be ordered so that abs(D./E) is in descending order. 0014 % If two output arguments are given then A will be scaled to make diag(E)=I 0015 % but, if W is singular, this will cause some elements of A to be infinite. 0016 % 0017 % If A'*B*A=diag(D) and A'*W*A=diag(E) then A'*W*A*diag(1./E)=I so A'*B*A=A'*W*A*diag(D./E) 0018 % and hence B*A=W*A*diag(D./E) so the columns of A are the eigenvectors of W\A or 0019 % equivalently the generalized eigenvectors of (B,W). 0020 % 0021 % Suppose we have several N-dimensional data row-vectors arising from each of C different classes of data. 0022 % for each class, c, we can form the mean data vector m(c) and the within-class covariance matrix W(c) 0023 % We can then form the between class covariance matrix B by taking the covariance of the mean vectors m(1), m(2), ... 0024 % and also the averaged within-class covariance matrix W by averaging W(1), W(2), ... 0025 % If we then take A=v_dualdiag(W,B) and postmultiply all our original data vectors by A, we obtain new 0026 % data vectors for which the average within-class covariance matrix is the identity and for which 0027 % the first few components contain most of the information that is useful in discriminating between classes. 0028 0029 % Copyright (C) Mike Brookes 1997-2013 0030 % Version: $Id: v_dualdiag.m 10865 2018-09-21 17:22:45Z dmb $ 0031 % 0032 % VOICEBOX is a MATLAB toolbox for speech processing. 0033 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0034 % 0035 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0036 % This program is free software; you can redistribute it and/or modify 0037 % it under the terms of the GNU General Public License as published by 0038 % the Free Software Foundation; either version 2 of the License, or 0039 % (at your option) any later version. 0040 % 0041 % This program is distributed in the hope that it will be useful, 0042 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0043 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0044 % GNU General Public License for more details. 0045 % 0046 % You can obtain a copy of the GNU General Public License from 0047 % http://www.gnu.org/copyleft/gpl.html or by writing to 0048 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0049 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0050 [a,l]=eig(b+b',w+w'); % generalized eigendecomposition 0051 if isreal(l) 0052 [d,i]=sort(abs(diag(l)),'descend'); % sort by absolute value 0053 if nargout==2 0054 q=sqrt(diag(a'*w*a))'.^(-1); % scale to make e=1 0055 a=a(:,i).*q(ones(size(w,1),1),i); % reorder and scale columns of a 0056 else 0057 a=a(:,i); % reorder columns of a to match eigenvalues 0058 e=real(diag(a'*w*a)); 0059 end 0060 d=real(diag(a'*b*a)); 0061 else 0062 d=a'*b*a; 0063 e=a'*w*a; 0064 end