source: distools/private/cuteigs.m @ 34

Last change on this file since 34 was 13, checked in by bduin, 14 years ago
File size: 2.4 KB
RevLine 
[13]1%CUTEIGS Select significant eigenvalues from a list
2%
3%   J = CUTEIGS(L,PREC)
4%
5% INPUT
6%   L     List of eigenvalues
7%   PREC  Precision parameter (optional; default: 0.0001)
8%
9% OUTPUT
10%   J   Index of selected eigenvalues
11%
12% DESCRIPTION
13% This is a low-level routine for SELEIGS, which serves PSEM, KPSEM, KPCA and
14% PEPCA. It makes use of PCHIP determing piecewise cubic Hermite interpolating
15% polynomial P, built on [1:length(LL) LL], where LL = -sort(-abs(L).^3) folloowed
16% by a normalization to [0,1]. A third power is used to emphasize the differences
17% better, as sometimes there is a very long tail of slowly dropping eigenvalues.
18% The number M of significant eigenvalues is determined such that the ratio
19% of the estimated derivative of P versus the largest derivative value is smaller
20% then the given precision, PREC * cumsum(LL).
21%
22% Note that M will differ between different samples and for different sizes of L.
23%
24% SEE ALSO
25% MAPPINGS, DATASETS, PSEM, KPSEM, KPCA, PEPCA
26%
27
28% Elzbieta Pekalska, e.pekalska@ewi.tudelft.nl
29% Faculty of Electrical Engineering, Mathematics and Computer Science,
30% Delft University of Technology, The Netherlands
31
32
33function [J,prec] = cuteigs(L,prec)
34
35if nargin < 2 | isempty(prec),
36  prec = 0.0001;
37end
38
39tol = 1e-14;
40
41n  = length(L);
42LL = L.^3;
43[lambda,J] = sort(-abs(LL)/sum(abs(LL)));
44lambda = -lambda;
45
46if n > 1,
47  M = findJ(lambda,prec);
48else
49  M = n;
50end
51JJ = J(1:M);
52
53I1 = find (L(JJ) >= 0);
54I2 = find (L(JJ) < 0);
55
56J  = [JJ(I1); JJ(I2)];
57
58
59
60
61
62
63
64function M = findJ (lambda,prec)
65
66n = length(lambda);
67if n >= 400,
68  Z = [1:80 81:4:n];
69elseif n >= 200,
70  Z = [1:50 51:2:n];
71else
72  Z =1:n;
73end
74
75if n > 3
76  % Hermite piecewise cubic interpolating polynomial
77  H = pchip(Z,lambda(Z));
78
79  % Derivative of the Hermite piecewise cubic polynomial
80  Hder = H;
81  Hder.order = 3;
82  Hder.coefs = zeros(H.pieces,3);
83  for i=1:Hder.pieces
84    pp = polyder(H.coefs(i,:));
85    Hder.coefs(i,4-length(pp):3) = pp;
86  end
87  pder= -ppval(1:n,Hder);   % Evaluated derivative
88else
89  pder= [lambda(1) -diff(lambda')];   % Approximated derivative
90end
91
92% pder(1) is an extrapolated value; pder(2) is used
93I  = find(pder/pder(2) < prec*cumsum(lambda'));
94II = diff(I);
95Q  = find(II > 1);
96if ~isempty(Q),
97  I = I(Q(1));
98end
99if ~isempty(I),
100  M = I(1)-1;
101else
102  if n <= 25,
103    M = n;
104  else
105    M = n-1;
106  end
107end
108if M < 1, M = 1; end
Note: See TracBrowser for help on using the repository browser.