SHOGUN  v1.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
KernelPCA.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2011 Soeren Sonnenburg
8  * Copyright (C) 2011 Berlin Institute of Technology
9  */
10 
12 #ifdef HAVE_LAPACK
13 #include <shogun/lib/config.h>
15 
16 #include <string.h>
17 #include <stdlib.h>
18 
20 #include <shogun/lib/common.h>
21 #include <shogun/kernel/Kernel.h>
25 #include <shogun/io/SGIO.h>
26 
27 using namespace shogun;
28 
30 {
31  init();
32 }
33 
35 {
36  init();
37  set_kernel(k);
38 }
39 
41 {
42  m_initialized = false;
43  m_init_features = NULL;
44  m_transformation_matrix = SGMatrix<float64_t>(NULL, 0, 0, false);
45  m_bias_vector = SGVector<float64_t>(NULL, 0, false);
46 
47  m_parameters->add(&m_transformation_matrix, "transformation matrix",
48  "matrix used to transform data");
49  m_parameters->add(&m_bias_vector, "bias vector",
50  "bias vector used to transform data");
51 }
52 
54 {
58  m_initialized = false;
59 }
60 
62 {
66 }
67 
68 bool CKernelPCA::init(CFeatures* features)
69 {
70  if (!m_initialized && m_kernel)
71  {
72  SG_REF(features);
73  m_init_features = features;
74 
75  m_kernel->init(features,features);
77  m_kernel->cleanup();
78  int32_t n = kernel_matrix.num_cols;
79  int32_t m = kernel_matrix.num_rows;
80  ASSERT(n==m);
81 
82  float64_t* bias_tmp = CMath::get_column_sum(kernel_matrix.matrix, n,n);
83  CMath::scale_vector(-1.0/n, bias_tmp, n);
84  float64_t s = CMath::sum(bias_tmp, n)/n;
85  CMath::add_scalar(-s, bias_tmp, n);
86 
87  CMath::center_matrix(kernel_matrix.matrix, n, m);
88 
89  float64_t* eigenvalues=CMath::compute_eigenvectors(kernel_matrix.matrix, n, n);
90 
91  for (int32_t i=0; i<n; i++)
92  {
93  //normalize and trap divide by zero and negative eigenvalues
94  for (int32_t j=0; j<n; j++)
95  kernel_matrix.matrix[i*n+j]/=CMath::sqrt(CMath::max(1e-16,eigenvalues[i]));
96  }
97 
98  SG_FREE(eigenvalues);
99 
102 
105  CMath::fill_vector(m_bias_vector.vector, m_bias_vector.vlen, 0.0);
106 
108  CblasTrans, bias_tmp, 0.0, m_bias_vector.vector);
109 
111  CMath::scale_vector(1.0/n, rowsum, n);
112 
113  for (int32_t i=0; i<n; i++)
114  {
115  for (int32_t j=0; j<n; j++)
116  m_transformation_matrix.matrix[j+n*i] -= rowsum[i];
117  }
118  SG_FREE(rowsum);
119  SG_FREE(bias_tmp);
120 
121  m_initialized=true;
122  SG_INFO("Done\n");
123  return true;
124  }
125  return false;
126 }
127 
128 
130 {
132  CSimpleFeatures<float64_t>* simple_features = (CSimpleFeatures<float64_t>*)features;
133 
134  int32_t num_vectors = simple_features->get_num_vectors();
135  int32_t i,j,k;
136  int32_t n = m_transformation_matrix.num_cols;
137 
138  m_kernel->init(features,m_init_features);
139 
140  float64_t* new_feature_matrix = SG_MALLOC(float64_t, m_target_dim*num_vectors);
141 
142  for (i=0; i<num_vectors; i++)
143  {
144  for (j=0; j<m_target_dim; j++)
145  new_feature_matrix[i*m_target_dim+j] = m_bias_vector.vector[j];
146 
147  for (j=0; j<n; j++)
148  {
149  float64_t kij = m_kernel->kernel(i,j);
150 
151  for (k=0; k<m_target_dim; k++)
152  new_feature_matrix[k+i*m_target_dim] += kij*m_transformation_matrix.matrix[(n-k-1)*n+j];
153  }
154  }
155 
156  m_kernel->cleanup();
157  simple_features->set_feature_matrix(SGMatrix<float64_t>(new_feature_matrix,m_target_dim,num_vectors));
158  return ((CSimpleFeatures<float64_t>*)features)->get_feature_matrix();
159 }
160 
162 {
167 
168  int32_t j,k;
169  int32_t n = m_transformation_matrix.num_cols;
170 
171  for (j=0; j<m_target_dim; j++)
172  result.vector[j] = m_bias_vector.vector[j];
173 
174  for (j=0; j<n; j++)
175  {
176  float64_t kj = m_kernel->kernel(0,j);
177 
178  for (k=0; k<m_target_dim; k++)
179  result.vector[k] += kj*m_transformation_matrix.matrix[(n-k-1)*n+j];
180  }
181 
182  m_kernel->cleanup();
183  return result;
184 }
185 
187 {
189 
190  int32_t num_vectors = features->get_num_vectors();
191  int32_t i,j,k;
192  int32_t n = m_transformation_matrix.num_cols;
193 
194  m_kernel->init(features,m_init_features);
195 
196  float64_t* new_feature_matrix = SG_MALLOC(float64_t, m_target_dim*num_vectors);
197 
198  for (i=0; i<num_vectors; i++)
199  {
200  for (j=0; j<m_target_dim; j++)
201  new_feature_matrix[i*m_target_dim+j] = m_bias_vector.vector[j];
202 
203  for (j=0; j<n; j++)
204  {
205  float64_t kij = m_kernel->kernel(i,j);
206 
207  for (k=0; k<m_target_dim; k++)
208  new_feature_matrix[k+i*m_target_dim] += kij*m_transformation_matrix.matrix[(n-k-1)*n+j];
209  }
210  }
211 
212  m_kernel->cleanup();
213 
214  return new CSimpleFeatures<float64_t>(SGMatrix<float64_t>(new_feature_matrix,m_target_dim,num_vectors));
215 }
216 
217 #endif

SHOGUN Machine Learning Toolbox - Documentation