1+ function  [eig_vectors , eig_values , eig_ratio ] =  computePCA(X , num_sub_dim , varargin )
2+  %  computePCA: compute PCA of X
3+  %  Input:
4+  %  X: num_dim x num_samples data matrix
5+  %  num_sub_dim: number of sub-dimensions to keep (Principal Components). if num_sub_dim < 1
6+  %  then it is the ratio of variance to keep
7+  %  varargin: if not specified PCA is computed on covariance matrix, 
8+  %  if set to 'R' then PCA is computed on correlation matrix
9+  %  Output:
10+  %  Z: projected data to the principal components
11+  %  eig_vectors: eigenvectors of the covariance/correlation matrix
12+  %  eig_values: eigenvalues of the covariance/correlation matrix
13+  %  eig_ratio: ratio of variance explained by each principal component
14+  %  Last Update: 2023/10 by Santos Enoque
15+  %  Computer vison laboratory, University of Tsukuba
16+  %  http://www.cvlab.cs.tsukuba.ac.jp/ 
17+  X =  X(: , : );
18+  [num_dim , num_samples ] =  size(X );
19+ 
20+  use_covariance_matrix =  true ;
21+  if  nargin  >  2 
22+  if  strcmp(varargin{1 }, ' R' 
23+  use_covariance_matrix =  false ;
24+  end 
25+  end 
26+ 
27+  apply_matrix_normalization =  use_covariance_matrix ;
28+ 
29+  if  num_sub_dim  <  1 
30+  c_ratio =  num_sub_dim ;
31+  if  num_dim  <  num_samples 
32+  if  use_covariance_matrix 
33+  C =  cov(X ' , 1 );
34+  else 
35+  C =  X * X ' /(num_samples - 1 );
36+  end 
37+  [eig_vectors , temp_diagonal ] =  eig(C );
38+  [eig_values , ind ] =  sort(diag(temp_diagonal ), ' descend' 
39+  eig_vectors =  eig_vectors(: , ind );
40+  num_sub_dim =  find(cumsum(eig_values )/sum(eig_values ) >=  c_ratio , 1 );
41+  eig_vectors =  eig_vectors(: , 1 : num_sub_dim );
42+  eig_values =  eig_values(1 : num_sub_dim );
43+  eig_ratio =  sum(eig_values )/trace(C );
44+  else 
45+  if  apply_matrix_normalization 
46+  K =  X ' *X ;
47+  IN =  ones(num_samples , num_samples )/num_samples ; 
48+  K =  K  -  IN * K  -  K * IN  +  IN * K * IN ;
49+  else 
50+  K =  X ' *X ;
51+  end 
52+  [A , B ] =  eig(K );
53+  [B , ind ] =  sort(diag(B ), ' descend' 
54+  A =  A(: , ind );
55+  eig_values =  B / num_samples ;
56+  num_sub_dim =  find(cumsum(eig_values )/sum(eig_values ) >=  c_ratio , 1 );
57+  A =  A(: , 1 : num_sub_dim );
58+  B =  B(1 : num_sub_dim );
59+  A =  A / sqrt(diag(B ));
60+  eig_vectors =  X * A ;
61+  eig_values =  eig_values(1 : num_sub_dim );
62+  eig_ratio =  sum(eig_values );
63+  end 
64+  elseif  num_sub_dim  >=  1  
65+  num_sub_dim =  floor(num_sub_dim );
66+  if  num_dim  <  num_samples 
67+  if  use_covariance_matrix 
68+  C =  cov(X ' , 1 );
69+  else 
70+  C =  X * X ' /(num_samples - 1 );
71+  end 
72+ 
73+  OPTS.disp =  0 ;
74+  if  num_sub_dim  <  num_dim 
75+  [eig_vectors , temp_diagonal ] =  eigs(C , num_sub_dim , ' LM' OPTS );
76+  else 
77+  [eig_vectors , temp_diagonal ] =  eig(C );
78+  end 
79+  [eig_values , ind ] =  sort(diag(temp_diagonal ), ' descend' 
80+  eig_vectors =  eig_vectors(: , ind );
81+  eig_ratio =  sum(eig_values )/trace(C );
82+  else 
83+  if  apply_matrix_normalization 
84+  K =  X ' *X ;
85+  IN =  ones(num_samples , num_samples )/num_samples ; 
86+  K =  K  -  IN * K  -  K * IN  +  IN * K * IN ;
87+  else 
88+  K =  X ' *X ;
89+  end 
90+  OPTS.disp =  0 ;
91+  [A , B ] =  eigs(K , num_sub_dim , ' LM' OPTS );
92+  [B , ind ] =  sort(diag(B ), ' descend' 
93+  A =  A(: , ind );
94+  eig_values =  B / num_samples ;
95+  A =  A / sqrt(diag(B ));
96+  eig_vectors =  X * A ;
97+  eig_ratio =  sum(eig_values );
98+  end 
99+  end 
100+  
0 commit comments