KernelHelpers.h
Go to the documentation of this file.
1 //===========================================================================
2 /*!
3  *
4  *
5  * \brief Collection of functions dealing with typical tasks of kernels
6 
7 
8  *
9  *
10  * \author O. Krause
11  * \date 2007-2012
12  *
13  *
14  * \par Copyright 1995-2017 Shark Development Team
15  *
16  * <BR><HR>
17  * This file is part of Shark.
18  * <http://shark-ml.org/>
19  *
20  * Shark is free software: you can redistribute it and/or modify
21  * it under the terms of the GNU Lesser General Public License as published
22  * by the Free Software Foundation, either version 3 of the License, or
23  * (at your option) any later version.
24  *
25  * Shark is distributed in the hope that it will be useful,
26  * but WITHOUT ANY WARRANTY; without even the implied warranty of
27  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
28  * GNU Lesser General Public License for more details.
29  *
30  * You should have received a copy of the GNU Lesser General Public License
31  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
32  *
33  */
34 //===========================================================================
35 
36 #ifndef SHARK_MODELS_KERNELS_KERNELHELPERS_H
37 #define SHARK_MODELS_KERNELS_KERNELHELPERS_H
38 
39 #include <shark/Models/Kernels/AbstractKernelFunction.h>
40 #include <shark/Data/Dataset.h>
41 #include <shark/Core/OpenMP.h>
42 namespace shark{
43 
44 /// \brief Calculates the regularized kernel gram matrix of the points stored inside a dataset.
45 ///
46 /// Regularization is applied by adding the regularizer on the diagonal
47 /// \param kernel the kernel for which to calculate the kernel gram matrix
48 /// \param dataset the set of points used in the gram matrix
49 /// \param matrix the target kernel matrix
50 /// \param regularizer the regularizer of the matrix which is always >= 0. default is 0.
51 template<class InputType, class M, class Device>
52 void calculateRegularizedKernelMatrix(
53  AbstractKernelFunction<InputType>const& kernel,
54  Data<InputType> const& dataset,
55  blas::matrix_expression<M, Device>& matrix,
56  double regularizer = 0
57 ){
58  SHARK_RUNTIME_CHECK(regularizer >= 0, "regularizer must be >=0");
59  std::size_t B = dataset.numberOfBatches();
60  //get start of all batches in the matrix
61  //also include the past the end position at the end
62  std::vector<std::size_t> batchStart(B+1,0);
63  for(std::size_t i = 1; i != B+1; ++i){
64  batchStart[i] = batchStart[i-1]+ batchSize(dataset.batch(i-1));
65  }
66  SIZE_CHECK(batchStart[B] == dataset.numberOfElements());
67  std::size_t N = batchStart[B];//number of elements
68  ensure_size(matrix,N,N);
69 
70 
71  for (std::size_t i=0; i<B; i++){
72  std::size_t startX = batchStart[i];
73  std::size_t endX = batchStart[i+1];
74  SHARK_PARALLEL_FOR(int j=0; j < (int)B; j++){
75  std::size_t startY = batchStart[j];
76  std::size_t endY = batchStart[j+1];
77  RealMatrix submatrix = kernel(dataset.batch(i), dataset.batch(j));
78  noalias(subrange(matrix(),startX,endX,startY,endY))=submatrix;
79  //~ if(i != j)
80  //~ noalias(subrange(matrix(),startY,endY,startX,endX))=trans(submatrix);
81  }
82  for(std::size_t k = startX; k != endX; ++k){
83  matrix()(k,k) += static_cast<typename M::value_type>(regularizer);
84  }
85  }
86 }
87 
88 /// \brief Calculates the kernel gram matrix between two data sets.
89 ///
90 /// \param kernel the kernel for which to calculate the kernel gram matrix
91 /// \param dataset1 the set of points corresponding to rows of the Gram matrix
92 /// \param dataset2 the set of points corresponding to columns of the Gram matrix
93 /// \param matrix the target kernel matrix
94 template<class InputType, class M, class Device>
95 void calculateMixedKernelMatrix(
96  AbstractKernelFunction<InputType>const& kernel,
97  Data<InputType> const& dataset1,
98  Data<InputType> const& dataset2,
99  blas::matrix_expression<M, Device>& matrix
100 ){
101  std::size_t B1 = dataset1.numberOfBatches();
102  std::size_t B2 = dataset2.numberOfBatches();
103  //get start of all batches in the matrix
104  //also include the past the end position at the end
105  std::vector<std::size_t> batchStart1(B1+1,0);
106  for(std::size_t i = 1; i != B1+1; ++i){
107  batchStart1[i] = batchStart1[i-1]+ batchSize(dataset1.batch(i-1));
108  }
109  std::vector<std::size_t> batchStart2(B2+1,0);
110  for(std::size_t i = 1; i != B2+1; ++i){
111  batchStart2[i] = batchStart2[i-1]+ batchSize(dataset2.batch(i-1));
112  }
113  SIZE_CHECK(batchStart1[B1] == dataset1.numberOfElements());
114  SIZE_CHECK(batchStart2[B2] == dataset2.numberOfElements());
115  std::size_t N1 = batchStart1[B1];//number of elements
116  std::size_t N2 = batchStart2[B2];//number of elements
117  ensure_size(matrix,N1,N2);
118 
119  for (std::size_t i=0; i<B1; i++){
120  std::size_t startX = batchStart1[i];
121  std::size_t endX = batchStart1[i+1];
122  SHARK_PARALLEL_FOR(int j=0; j < (int)B2; j++){
123  std::size_t startY = batchStart2[j];
124  std::size_t endY = batchStart2[j+1];
125  RealMatrix submatrix = kernel(dataset1.batch(i), dataset2.batch(j));
126  noalias(subrange(matrix(),startX,endX,startY,endY))=submatrix;
127  //~ if(i != j)
128  //~ noalias(subrange(matrix(),startY,endY,startX,endX))=trans(submatrix);
129  }
130  }
131 }
132 
133 /// \brief Calculates the regularized kernel gram matrix of the points stored inside a dataset.
134 ///
135 /// Regularization is applied by adding the regularizer on the diagonal
136 /// \param kernel the kernel for which to calculate the kernel gram matrix
137 /// \param dataset the set of points used in the gram matrix
138 /// \param regularizer the regularizer of the matrix which is always >= 0. default is 0.
139 /// \return the kernel gram matrix
140 template<class InputType>
141 RealMatrix calculateRegularizedKernelMatrix(
142  AbstractKernelFunction<InputType>const& kernel,
143  Data<InputType> const& dataset,
144  double regularizer = 0
145 ){
146  SHARK_RUNTIME_CHECK(regularizer >= 0, "regularizer must be >=0");
147  RealMatrix M;
148  calculateRegularizedKernelMatrix(kernel,dataset,M,regularizer);
149  return M;
150 }
151 
152 /// \brief Calculates the kernel gram matrix between two data sets.
153 ///
154 /// \param kernel the kernel for which to calculate the kernel gram matrix
155 /// \param dataset1 the set of points corresponding to rows of the Gram matrix
156 /// \param dataset2 the set of points corresponding to columns of the Gram matrix
157 /// \return matrix the target kernel matrix
158 template<class InputType>
159 RealMatrix calculateMixedKernelMatrix(
160  AbstractKernelFunction<InputType>const& kernel,
161  Data<InputType> const& dataset1,
162  Data<InputType> const& dataset2
163 ){
164  RealMatrix M;
165  calculateMixedKernelMatrix(kernel,dataset1,dataset2,M);
166  return M;
167 }
168 
169 
170 /// \brief Efficiently calculates the weighted derivative of a Kernel Gram Matrix w.r.t the Kernel Parameters
171 ///
172 /// The formula is \f$ \sum_i \sum_j w_{ij} k(x_i,x_j)\f$ where w_ij are the weights of the gradient and x_i x_j are
173 /// the datapoints defining the gram matrix and k is the kernel. For efficiency it is assumd that w_ij = w_ji.
174 ///This method is only useful when the whole Kernel Gram Matrix neds to be computed to get the weights w_ij and
175 ///only computing smaller blocks is not sufficient.
176 /// \param kernel the kernel for which to calculate the kernel gram matrix
177 /// \param dataset the set of points used in the gram matrix
178 /// \param weights the weights of the derivative, they must be symmetric!
179 /// \return the weighted derivative w.r.t the parameters.
180 template<class InputType,class WeightMatrix>
181 RealVector calculateKernelMatrixParameterDerivative(
182  AbstractKernelFunction<InputType> const& kernel,
183  Data<InputType> const& dataset,
184  WeightMatrix const& weights
185 ){
186  std::size_t kp = kernel.numberOfParameters();
187  RealMatrix block;//stores the kernel results of the block which we need to compute to get the State :(
188  RealVector kernelGradient(kp);//weighted gradient summed over the whole kernel matrix
189  kernelGradient.clear();
190 
191  //calculate the gradint blockwise taking symmetry into account.
192  RealVector blockGradient(kp);//weighted gradient summed over the whole block
193  boost::shared_ptr<State> state = kernel.createState();
194  std::size_t startX = 0;
195  for (std::size_t i=0; i<dataset.numberOfBatches(); i++){
196  std::size_t sizeX= batchSize(dataset.batch(i));
197  std::size_t startY = 0;//start of the current batch in y-direction
198  for (std::size_t j=0; j <= i; j++){
199  std::size_t sizeY= batchSize(dataset.batch(j));
200  kernel.eval(dataset.batch(i), dataset.batch(j),block,*state);
201  kernel.weightedParameterDerivative(
202  dataset.batch(i), dataset.batch(j),//points
203  subrange(weights,startX,startX+sizeX,startY,startY+sizeY),//weights
204  *state,
205  blockGradient
206  );
207  if(i != j)
208  kernelGradient+=2*blockGradient;//Symmetry!
209  else
210  kernelGradient+=blockGradient;//middle blocks are symmetric
211  startY+= sizeY;
212  }
213  startX+=sizeX;
214  }
215  return kernelGradient;
216 }
217 
218 }
219 #endif