LDA.m
上传用户:wo4099
上传日期:2008-12-29
资源大小:1k
文件大小:2k
源码类别:

生物技术

开发平台:

Matlab

  1. function [eigvector, eigvalue, Y] = LDA(X,gnd)
  2. % LDA: Linear discriminant analysis (Fisherfaces approach PCA+LDA)
  3. %
  4. %       [eigvector, eigvalue] = LDA(X, gnd)
  5. %
  6. %             Input:
  7. %               X     - Data matrix. Each row vector of fea is a data point.
  8. %               gnd   - Colunm vector of the label information for each
  9. %                       data point. 
  10. %
  11. %             Output:
  12. %               eigvector - Each column is an embedding function, for a new
  13. %                           data point (row vector) x,  y = x*eigvector
  14. %                           will be the embedding result of x.
  15. %               eigvalue  - The eigvalue of LDA eigen-problem. 
  16. %       [eigvector, eigvalue, Y] = LDA(X, gnd) 
  17. %               
  18. %               Y:  The embedding results, Each row vector is a data point.
  19. %                   Y = X*eigvector
  20. %
  21. % Reference:
  22. %   
  23. %         P. N. Belhumeur, J. P. Hespanha, and D. J. Kriegman, eigenfaces
  24. %         vs. fisherfaces: recognition using class specific linear
  25. %         projection,?IEEE Transactions on Pattern Analysis and Machine
  26. %         Intelligence, vol. 19, no. 7, pp. 711-720, July 1997.  
  27. %
  28. %    Written by Deng Cai (dengcai@gmail.com), April/2004, Feb/2006
  29. old_X = X;
  30. % ====== Initialization
  31. [nSmp,nFea] = size(X);
  32. classLabel = unique(gnd);
  33. nClass = length(classLabel);
  34. bPCA = 0;
  35. if nFea > (nSmp - nClass)
  36.     PCAoptions = [];
  37.     PCAoptions.ReducedDim = nSmp - nClass;
  38.     [eigvector_PCA, eigvalue_PCA, meanData, new_X] = PCA(X,PCAoptions);
  39.     X = new_X;
  40.     [nSmp,nFea] = size(X);
  41.     bPCA = 1;
  42. end
  43. sampleMean = mean(X);
  44. MMM = zeros(nFea, nFea);
  45. for i = 1:nClass
  46. index = find(gnd==classLabel(i));
  47. classMean = mean(X(index, :));
  48. MMM = MMM + length(index)*classMean'*classMean;
  49. end
  50. W = X'*X - MMM;
  51. B = MMM - nSmp*sampleMean'*sampleMean;
  52. W = (W + W')/2;
  53. B = (B + B')/2;
  54. option = struct('disp',0);
  55. [eigvector, eigvalue] = eigs(B,W,nClass-1,'la',option);
  56. eigvalue = diag(eigvalue);
  57. for i = 1:size(eigvector,2)
  58.     eigvector(:,i) = eigvector(:,i)./norm(eigvector(:,i));
  59. end
  60. if bPCA 
  61.     eigvector =eigvector_PCA*eigvector;
  62. end
  63. if nargout == 3
  64.    Y = old_X * eigvector;
  65. end