SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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>
24 #include <shogun/io/SGIO.h>
25 
26 using namespace shogun;
27 
29 {
30  init();
31 }
32 
34 {
35  init();
36  set_kernel(k);
37 }
38 
40 {
41  m_initialized = false;
42  m_init_features = NULL;
43  m_transformation_matrix = SGMatrix<float64_t>(NULL, 0, 0, false);
44  m_bias_vector = SGVector<float64_t>(NULL, 0, false);
45 
46  SG_ADD(&m_transformation_matrix, "transformation_matrix",
47  "matrix used to transform data", MS_NOT_AVAILABLE);
48  SG_ADD(&m_bias_vector, "bias_vector",
49  "bias vector used to transform data", MS_NOT_AVAILABLE);
50 }
51 
53 {
55 
56  if (m_init_features)
58 
59  m_initialized = false;
60 }
61 
63 {
64  if (m_init_features)
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 = SGMatrix<float64_t>::get_column_sum(kernel_matrix.matrix, n,n);
83  SGVector<float64_t>::scale_vector(-1.0/n, bias_tmp, n);
84  float64_t s = SGVector<float64_t>::sum(bias_tmp, n)/n;
85  SGVector<float64_t>::add_scalar(-s, bias_tmp, n);
86 
87  SGMatrix<float64_t>::center_matrix(kernel_matrix.matrix, n, m);
88 
89  float64_t* eigenvalues=SGMatrix<float64_t>::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 
101  kernel_matrix.matrix = NULL;
104 
105  cblas_dgemv(CblasColMajor, CblasTrans,
106  n, n, 1.0, m_transformation_matrix.matrix, n,
107  bias_tmp, 1, 0.0, m_bias_vector.vector, 1);
108 
110  SGVector<float64_t>::scale_vector(1.0/n, rowsum, n);
111 
112  for (int32_t i=0; i<n; i++)
113  {
114  for (int32_t j=0; j<n; j++)
115  m_transformation_matrix.matrix[j+n*i] -= rowsum[i];
116  }
117  SG_FREE(rowsum);
118  SG_FREE(bias_tmp);
119 
120  m_initialized=true;
121  SG_INFO("Done\n")
122  return true;
123  }
124  return false;
125 }
126 
127 
129 {
131  CDenseFeatures<float64_t>* simple_features = (CDenseFeatures<float64_t>*)features;
132 
133  int32_t num_vectors = simple_features->get_num_vectors();
134  int32_t i,j,k;
135  int32_t n = m_transformation_matrix.num_cols;
136 
137  m_kernel->init(features,m_init_features);
138 
139  float64_t* new_feature_matrix = SG_MALLOC(float64_t, m_target_dim*num_vectors);
140 
141  for (i=0; i<num_vectors; i++)
142  {
143  for (j=0; j<m_target_dim; j++)
144  new_feature_matrix[i*m_target_dim+j] = m_bias_vector.vector[j];
145 
146  for (j=0; j<n; j++)
147  {
148  float64_t kij = m_kernel->kernel(i,j);
149 
150  for (k=0; k<m_target_dim; k++)
151  new_feature_matrix[k+i*m_target_dim] += kij*m_transformation_matrix.matrix[(n-k-1)*n+j];
152  }
153  }
154 
155  m_kernel->cleanup();
156  simple_features->set_feature_matrix(SGMatrix<float64_t>(new_feature_matrix,m_target_dim,num_vectors));
157  return ((CDenseFeatures<float64_t>*)features)->get_feature_matrix();
158 }
159 
161 {
166 
167  int32_t j,k;
168  int32_t n = m_transformation_matrix.num_cols;
169 
170  for (j=0; j<m_target_dim; j++)
171  result.vector[j] = m_bias_vector.vector[j];
172 
173  for (j=0; j<n; j++)
174  {
175  float64_t kj = m_kernel->kernel(0,j);
176 
177  for (k=0; k<m_target_dim; k++)
178  result.vector[k] += kj*m_transformation_matrix.matrix[(n-k-1)*n+j];
179  }
180 
181  m_kernel->cleanup();
182  return result;
183 }
184 
186 {
188 
189  int32_t num_vectors = features->get_num_vectors();
190  int32_t i,j,k;
191  int32_t n = m_transformation_matrix.num_cols;
192 
193  m_kernel->init(features,m_init_features);
194 
195  float64_t* new_feature_matrix = SG_MALLOC(float64_t, m_target_dim*num_vectors);
196 
197  for (i=0; i<num_vectors; i++)
198  {
199  for (j=0; j<m_target_dim; j++)
200  new_feature_matrix[i*m_target_dim+j] = m_bias_vector.vector[j];
201 
202  for (j=0; j<n; j++)
203  {
204  float64_t kij = m_kernel->kernel(i,j);
205 
206  for (k=0; k<m_target_dim; k++)
207  new_feature_matrix[k+i*m_target_dim] += kij*m_transformation_matrix.matrix[(n-k-1)*n+j];
208  }
209  }
210 
211  m_kernel->cleanup();
212 
213  return new CDenseFeatures<float64_t>(SGMatrix<float64_t>(new_feature_matrix,m_target_dim,num_vectors));
214 }
215 
216 #endif

SHOGUN Machine Learning Toolbox - Documentation