SHOGUN  6.1.3
KRRNystrom.cpp
Go to the documentation of this file.
1 
2 /*
3  * Copyright (c) The Shogun Machine Learning Toolbox
4  * Written (W) 2016 Fredrik Hallgren
5  * All rights reserved.
6  *
7  * Redistribution and use in source and binary forms, with or without
8  * modification, are permitted provided that the following conditions are met:
9  *
10  * 1. Redistributions of source code must retain the above copyright notice, this
11  * list of conditions and the following disclaimer.
12  * 2. Redistributions in binary form must reproduce the above copyright notice,
13  * this list of conditions and the following disclaimer in the documentation
14  * and/or other materials provided with the distribution.
15  *
16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
17  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
18  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
19  * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
20  * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
21  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
22  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
23  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26  *
27  * The views and conclusions contained in the software and documentation are those
28  * of the authors and should not be interpreted as representing official policies,
29  * either expressed or implied, of the Shogun Development Team.
30  */
31 
32 #include <limits>
37 
38 using namespace shogun;
39 using namespace Eigen;
40 
42 {
43  init();
44 }
45 
47 : CKernelRidgeRegression(tau, k, lab)
48 {
49  init();
50 
52 
53  int32_t n=kernel->get_num_vec_lhs();
54 
55  REQUIRE(m_num_rkhs_basis <= n, "Number of sampled rows (%d) must be \
56 less than number of data points (%d)\n", m_num_rkhs_basis, n);
57 }
58 
59 void CKRRNystrom::init()
60 {
62 }
63 
65 {
66  int32_t n=kernel->get_num_vec_lhs();
67  SGVector<int32_t> temp(n);
68  temp.range_fill();
69  CMath::permute(temp);
71  for (index_t i=0; i<m_num_rkhs_basis; ++i)
72  col[i]=temp[i];
73  CMath::qsort(col.vector, m_num_rkhs_basis);
74 
75  return col;
76 }
77 
79 {
80  int32_t n=kernel->get_num_vec_lhs();
81 
82  if (m_num_rkhs_basis == 0)
83  {
84  set_num_rkhs_basis((int32_t)CMath::ceil(n/2.0));
85  SG_SWARNING("Number of sampled rows not set, default is half (%d) "
86  "of the number of data points (%d)\n", m_num_rkhs_basis, n);
87  }
88 
90  if (y==NULL)
91  SG_ERROR("Labels not set.\n");
95  #pragma omp parallel for
96  for (index_t j=0; j<m_num_rkhs_basis; ++j)
97  {
98  for (index_t i=0; i<n; ++i)
99  K_nm(i,j)=kernel->kernel(i,col[j]);
100  }
101  #pragma omp parallel for
102  for (index_t i=0; i<m_num_rkhs_basis; ++i)
103  sg_memcpy(K_mm.matrix+i*m_num_rkhs_basis, K_nm.get_row_vector(col[i]), m_num_rkhs_basis*sizeof(float64_t));
104 
105  Map<MatrixXd> K_mm_eig(K_mm.matrix, m_num_rkhs_basis, m_num_rkhs_basis);
106  Map<MatrixXd> K_nm_eig(K_nm.matrix, n, m_num_rkhs_basis);
107  MatrixXd K_mn_eig = K_nm_eig.transpose();
108  Map<VectorXd> y_eig(y.vector, n);
109  VectorXd alphas_eig(m_num_rkhs_basis);
110 
111  /* Calculate the Moore-Penrose pseudoinverse */
112  MatrixXd Kplus=K_mn_eig*K_nm_eig+m_tau*K_mm_eig;
113  SelfAdjointEigenSolver<MatrixXd> solver(Kplus);
114  if (solver.info()!=Success)
115  {
116  SG_WARNING("Eigendecomposition failed.\n")
117  return false;
118  }
119 
120  /* Solve the system for alphas */
121  MatrixXd D=solver.eigenvalues().asDiagonal();
122  MatrixXd eigvec=solver.eigenvectors();
123  float64_t dbl_epsilon=std::numeric_limits<float64_t>::epsilon();
124  const float64_t tolerance=m_num_rkhs_basis*dbl_epsilon*D.maxCoeff();
125  for (index_t i=0; i<m_num_rkhs_basis; ++i)
126  {
127  if (D(i,i)<tolerance)
128  D(i,i)=0;
129  else
130  D(i,i)=1/D(i,i);
131  }
132  MatrixXd pseudoinv=eigvec*D*eigvec.transpose();
133  alphas_eig=pseudoinv*K_mn_eig*y_eig;
134 
135  /* Expand alpha with zeros to size n */
136  SGVector<float64_t> alpha_n(n);
137  alpha_n.zero();
138  for (index_t i=0; i<m_num_rkhs_basis; ++i)
139  alpha_n[col[i]]=alphas_eig[i];
140  m_alpha=alpha_n;
141 
142  return true;
143 }
int32_t m_num_rkhs_basis
Definition: KRRNystrom.h:119
static void permute(SGVector< T > v, CRandom *rand=NULL)
Definition: Math.h:962
Class KernelRidgeRegression implements Kernel Ridge Regression - a regularized least square method fo...
SGVector< T > get_row_vector(index_t row) const
Definition: SGMatrix.cpp:1211
Real Labels are real-valued labels.
int32_t index_t
Definition: common.h:72
static float64_t ceil(float64_t d)
Definition: Math.h:384
SGVector< int32_t > subsample_indices()
Definition: KRRNystrom.cpp:64
The class Labels models labels, i.e. class assignments of objects.
Definition: Labels.h:43
#define SG_SWARNING(...)
Definition: SGIO.h:163
Definition: SGMatrix.h:25
CLabels * m_labels
Definition: Machine.h:436
#define SG_ERROR(...)
Definition: SGIO.h:128
#define REQUIRE(x,...)
Definition: SGIO.h:181
float64_t kernel(int32_t idx_a, int32_t idx_b)
virtual int32_t get_num_vec_lhs()
SGVector< float64_t > m_alpha
virtual bool solve_krr_system()
Definition: KRRNystrom.cpp:78
double float64_t
Definition: common.h:60
void range_fill(T start=0)
Definition: SGVector.cpp:223
void set_num_rkhs_basis(int32_t m)
Definition: KRRNystrom.h:88
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
The Kernel base class.
#define SG_WARNING(...)
Definition: SGIO.h:127
static void qsort(T *output, int32_t size)
Definition: Math.h:1134

SHOGUN Machine Learning Toolbox - Documentation