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 * <https://shark-ml.github.io/Shark/>
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>
42namespace 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/// \ingroup kernels
52template<class InputType, class M, class Device>
53void calculateRegularizedKernelMatrix(
54 AbstractKernelFunction<InputType>const& kernel,
55 Data<InputType> const& dataset,
56 blas::matrix_expression<M, Device>& matrix,
57 double regularizer = 0
58){
59 SHARK_RUNTIME_CHECK(regularizer >= 0, "regularizer must be >=0");
60 std::size_t B = dataset.numberOfBatches();
61 //get start of all batches in the matrix
62 //also include the past the end position at the end
63 std::vector<std::size_t> batchStart(B+1,0);
64 for(std::size_t i = 1; i != B+1; ++i){
65 batchStart[i] = batchStart[i-1]+ batchSize(dataset.batch(i-1));
66 }
67 SIZE_CHECK(batchStart[B] == dataset.numberOfElements());
68 std::size_t N = batchStart[B];//number of elements
69 ensure_size(matrix,N,N);
70
71
72 for (std::size_t i=0; i<B; i++){
73 std::size_t startX = batchStart[i];
74 std::size_t endX = batchStart[i+1];
75 SHARK_PARALLEL_FOR(int j=0; j < (int)B; j++){
76 std::size_t startY = batchStart[j];
77 std::size_t endY = batchStart[j+1];
78 RealMatrix submatrix = kernel(dataset.batch(i), dataset.batch(j));
79 noalias(subrange(matrix(),startX,endX,startY,endY))=submatrix;
80 //~ if(i != j)
81 //~ noalias(subrange(matrix(),startY,endY,startX,endX))=trans(submatrix);
82 }
83 for(std::size_t k = startX; k != endX; ++k){
84 matrix()(k,k) += static_cast<typename M::value_type>(regularizer);
85 }
86 }
87}
88
89/// \brief Calculates the kernel gram matrix between two data sets.
90///
91/// \param kernel the kernel for which to calculate the kernel gram matrix
92/// \param dataset1 the set of points corresponding to rows of the Gram matrix
93/// \param dataset2 the set of points corresponding to columns of the Gram matrix
94/// \param matrix the target kernel matrix
95/// \ingroup kernels
96template<class InputType, class M, class Device>
97void calculateMixedKernelMatrix(
98 AbstractKernelFunction<InputType>const& kernel,
99 Data<InputType> const& dataset1,
100 Data<InputType> const& dataset2,
101 blas::matrix_expression<M, Device>& matrix
102){
103 std::size_t B1 = dataset1.numberOfBatches();
104 std::size_t B2 = dataset2.numberOfBatches();
105 //get start of all batches in the matrix
106 //also include the past the end position at the end
107 std::vector<std::size_t> batchStart1(B1+1,0);
108 for(std::size_t i = 1; i != B1+1; ++i){
109 batchStart1[i] = batchStart1[i-1]+ batchSize(dataset1.batch(i-1));
110 }
111 std::vector<std::size_t> batchStart2(B2+1,0);
112 for(std::size_t i = 1; i != B2+1; ++i){
113 batchStart2[i] = batchStart2[i-1]+ batchSize(dataset2.batch(i-1));
114 }
115 SIZE_CHECK(batchStart1[B1] == dataset1.numberOfElements());
116 SIZE_CHECK(batchStart2[B2] == dataset2.numberOfElements());
117 std::size_t N1 = batchStart1[B1];//number of elements
118 std::size_t N2 = batchStart2[B2];//number of elements
119 ensure_size(matrix,N1,N2);
120
121 for (std::size_t i=0; i<B1; i++){
122 std::size_t startX = batchStart1[i];
123 std::size_t endX = batchStart1[i+1];
124 SHARK_PARALLEL_FOR(int j=0; j < (int)B2; j++){
125 std::size_t startY = batchStart2[j];
126 std::size_t endY = batchStart2[j+1];
127 RealMatrix submatrix = kernel(dataset1.batch(i), dataset2.batch(j));
128 noalias(subrange(matrix(),startX,endX,startY,endY))=submatrix;
129 //~ if(i != j)
130 //~ noalias(subrange(matrix(),startY,endY,startX,endX))=trans(submatrix);
131 }
132 }
133}
134
135/// \brief Calculates the regularized kernel gram matrix of the points stored inside a dataset.
136///
137/// Regularization is applied by adding the regularizer on the diagonal
138/// \param kernel the kernel for which to calculate the kernel gram matrix
139/// \param dataset the set of points used in the gram matrix
140/// \param regularizer the regularizer of the matrix which is always >= 0. default is 0.
141/// \return the kernel gram matrix
142/// \ingroup kernels
143template<class InputType>
144RealMatrix calculateRegularizedKernelMatrix(
145 AbstractKernelFunction<InputType>const& kernel,
146 Data<InputType> const& dataset,
147 double regularizer = 0
148){
149 SHARK_RUNTIME_CHECK(regularizer >= 0, "regularizer must be >=0");
150 RealMatrix M;
151 calculateRegularizedKernelMatrix(kernel,dataset,M,regularizer);
152 return M;
153}
154
155/// \brief Calculates the kernel gram matrix between two data sets.
156///
157/// \param kernel the kernel for which to calculate the kernel gram matrix
158/// \param dataset1 the set of points corresponding to rows of the Gram matrix
159/// \param dataset2 the set of points corresponding to columns of the Gram matrix
160/// \return matrix the target kernel matrix
161/// \ingroup kernels
162template<class InputType>
163RealMatrix calculateMixedKernelMatrix(
164 AbstractKernelFunction<InputType>const& kernel,
165 Data<InputType> const& dataset1,
166 Data<InputType> const& dataset2
167){
168 RealMatrix M;
169 calculateMixedKernelMatrix(kernel,dataset1,dataset2,M);
170 return M;
171}
172
173
174/// \brief Efficiently calculates the weighted derivative of a Kernel Gram Matrix w.r.t the Kernel Parameters
175///
176/// 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
177/// the datapoints defining the gram matrix and k is the kernel. For efficiency it is assumd that w_ij = w_ji.
178///This method is only useful when the whole Kernel Gram Matrix neds to be computed to get the weights w_ij and
179///only computing smaller blocks is not sufficient.
180/// \param kernel the kernel for which to calculate the kernel gram matrix
181/// \param dataset the set of points used in the gram matrix
182/// \param weights the weights of the derivative, they must be symmetric!
183/// \return the weighted derivative w.r.t the parameters.
184/// \ingroup kernels
185template<class InputType,class WeightMatrix>
186RealVector calculateKernelMatrixParameterDerivative(
187 AbstractKernelFunction<InputType> const& kernel,
188 Data<InputType> const& dataset,
189 WeightMatrix const& weights
190){
191 std::size_t kp = kernel.numberOfParameters();
192 RealMatrix block;//stores the kernel results of the block which we need to compute to get the State :(
193 RealVector kernelGradient(kp);//weighted gradient summed over the whole kernel matrix
194 kernelGradient.clear();
195
196 //calculate the gradint blockwise taking symmetry into account.
197 RealVector blockGradient(kp);//weighted gradient summed over the whole block
198 boost::shared_ptr<State> state = kernel.createState();
199 std::size_t startX = 0;
200 for (std::size_t i=0; i<dataset.numberOfBatches(); i++){
201 std::size_t sizeX= batchSize(dataset.batch(i));
202 std::size_t startY = 0;//start of the current batch in y-direction
203 for (std::size_t j=0; j <= i; j++){
204 std::size_t sizeY= batchSize(dataset.batch(j));
205 kernel.eval(dataset.batch(i), dataset.batch(j),block,*state);
206 kernel.weightedParameterDerivative(
207 dataset.batch(i), dataset.batch(j),//points
208 subrange(weights,startX,startX+sizeX,startY,startY+sizeY),//weights
209 *state,
210 blockGradient
211 );
212 if(i != j)
213 kernelGradient+=2*blockGradient;//Symmetry!
214 else
215 kernelGradient+=blockGradient;//middle blocks are symmetric
216 startY+= sizeY;
217 }
218 startX+=sizeX;
219 }
220 return kernelGradient;
221}
222
223}
224#endif