求助：核主成分KPCA的程序代码

kaze

小女子目前正在写一片关于绩效评价的论文，需要用到核主成分分析法，不知道有没有高手曾今做过这种程序，希望不吝指教。。。在下先谢了

[ 本帖最后由 lxq 于 2006-9-15 08:27 编辑 ]

happy

KPCA的提出者亲自写的程序

1. % Kernel PCA toy example for k(x,y)=exp(-||x-y||^2/rbf_var), cf. Fig. 4 in
   
2. % @article{SchSmoMue98,
   
3. %   author    = "B.~{Sch"olkopf} and A.~Smola and K.-R.~{M"uller}",
   
4. %   title     = "Nonlinear component analysis as a kernel Eigenvalue problem",
   
5. %   journal =         {Neural Computation},
   
6. %   volume    = 10,
   
7. %   issue     = 5,
   
8. %   pages     = "1299 -- 1319",
   
9. %   year      = 1998}
   
10. % This file can be downloaded from http://www.kernel-machines.org.
    
11. % Last modified: 4 July 2003
    
12. 
    
13. 
    
14. % parameters
    
15. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
16. rbf_var = 0.1;
    
17. xnum = 4;
    
18. ynum = 2;
    
19. max_ev = xnum*ynum;
    
20. % (extract features from the first <max_ev> Eigenvectors)
    
21. x_test_num = 15;
    
22. y_test_num = 15;
    
23. cluster_pos = [-0.5 -0.2; 0 0.6; 0.5 0];
    
24. cluster_size = 30;
    
25. 
    
26. % generate a toy data set
    
27. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
28. num_clusters = size(cluster_pos,1);
    
29. train_num = num_clusters*cluster_size;
    
30. patterns = zeros(train_num, 2);
    
31. range = 1;
    
32. randn('seed', 0);
    
33. for i=1:num_clusters,
    
34.   patterns((i-1)*cluster_size+1:i*cluster_size,1) = cluster_pos(i,1)+0.1*randn(cluster_size,1);
    
35.   patterns((i-1)*cluster_size+1:i*cluster_size,2) = cluster_pos(i,2)+0.1*randn(cluster_size,1);
    
36. end
    
37. test_num = x_test_num*y_test_num;
    
38. x_range = -range:(2*range/(x_test_num - 1)):range;
    
39. y_offset = 0.5;
    
40. y_range = -range+ y_offset:(2*range/(y_test_num - 1)):range+ y_offset;
    
41. [xs, ys] = meshgrid(x_range, y_range);
    
42. test_patterns(:, 1) = xs(:);
    
43. test_patterns(:, 2) = ys(:);
    
44. cov_size = train_num;  % use all patterns to compute the covariance matrix
    
45. 
    
46. % carry out Kernel PCA
    
47. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
48. for i=1:cov_size,
    
49.   for j=i:cov_size,
    
50.     K(i,j) = exp(-norm(patterns(i,:)-patterns(j,:))^2/rbf_var);
    
51.     K(j,i) = K(i,j);
    
52.   end
    
53. end
    
54. unit = ones(cov_size, cov_size)/cov_size;
    
55. % centering in feature space!
    
56. K_n = K - unit*K - K*unit + unit*K*unit;
    
57. 
    
58. [evecs,evals] = eig(K_n);
    
59. evals = real(diag(evals));
    
60. for i=1:cov_size,
    
61.   evecs(:,i) = evecs(:,i)/(sqrt(evals(i)));
    
62. end
    
63. 
    
64. % extract features
    
65. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
66. %  do not need the following here - only need test point features
    
67. %unit_train = ones(train_num,cov_size)/cov_size;
    
68. %for i=1:train_num,
    
69. %  for j=1:cov_size,
    
70. %    K_train(i,j) = exp(-norm(patterns(i,:)-patterns(j,:))^2/rbf_var);
    
71. %  end
    
72. %end
    
73. %K_train_n = K_train - unit_train*K - K_train*unit + unit_train*K*unit;
    
74. %features = zeros(train_num, max_ev);
    
75. %features = K_train_n * evecs(:,1:max_ev);
    
76. 
    
77. unit_test = ones(test_num,cov_size)/cov_size;
    
78. K_test = zeros(test_num,cov_size);
    
79. for i=1:test_num,
    
80.   for j=1:cov_size,
    
81.     K_test(i,j) = exp(-norm(test_patterns(i,:)-patterns(j,:))^2/rbf_var);
    
82.   end
    
83. end
    
84. K_test_n = K_test - unit_test*K - K_test*unit + unit_test*K*unit;
    
85. test_features = zeros(test_num, max_ev);
    
86. test_features = K_test_n * evecs(:,1:max_ev);
    
87. 
    
88. % plot it
    
89. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
90. figure(1); clf
    
91. for n = 1:max_ev,
    
92.   subplot(ynum, xnum, n);
    
93.   axis([-range range -range+y_offset range+y_offset]);
    
94.   imag = reshape(test_features(:,n), y_test_num, x_test_num);
    
95.   axis('xy')
    
96.   colormap(gray);
    
97.   hold on;
    
98.   pcolor(x_range, y_range, imag);
    
99.   shading interp
    
100.   contour(x_range, y_range, imag, 9, 'b');
     
101.   plot(patterns(:,1), patterns(:,2), 'r.')
     
102.   text(-1,1.65,sprintf('Eigenvalue=%4.3f', evals(n)));
     
103.   hold off;
     
104. end

复制代码

fishermanymc

这个代码太经典太经典 每次看都陶醉

studyias

请问kpca的提出者是谁？有网站没有？

dyhkxydfbb

studyias 发表于 2012-4-22 20:04

                               
登录/注册后可看大图

请问kpca的提出者是谁？有网站没有？

http://wenku.baidu.com/view/6ea8c81a6bd97f192279e9a6.html

有原始论文