KernelExpansion.h
Go to the documentation of this file.
1 //===========================================================================
2 /*!
3  *
4  *
5  * \brief Affine linear kernel function expansion
6  *
7  * \par
8  * Affine linear kernel expansions resulting from Support
9  * vector machine (SVM) training and other kernel methods.
10  *
11  *
12  *
13  *
14  * \author T. Glasmachers
15  * \date 2007-2011
16  *
17  *
18  * \par Copyright 1995-2015 Shark Development Team
19  *
20  * <BR><HR>
21  * This file is part of Shark.
22  * <http://image.diku.dk/shark/>
23  *
24  * Shark is free software: you can redistribute it and/or modify
25  * it under the terms of the GNU Lesser General Public License as published
26  * by the Free Software Foundation, either version 3 of the License, or
27  * (at your option) any later version.
28  *
29  * Shark is distributed in the hope that it will be useful,
30  * but WITHOUT ANY WARRANTY; without even the implied warranty of
31  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
32  * GNU Lesser General Public License for more details.
33  *
34  * You should have received a copy of the GNU Lesser General Public License
35  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
36  *
37  */
38 //===========================================================================
39 
40 
41 #ifndef SHARK_MODELS_KERNELEXPANSION_H
42 #define SHARK_MODELS_KERNELEXPANSION_H
43 
44 #include <shark/Models/Converter.h>
46 #include <shark/Data/Dataset.h>
47 #include <shark/Data/DataView.h>
48 
49 
50 namespace shark {
51 
52 
53 ///
54 /// \brief Linear model in a kernel feature space.
55 ///
56 /// An affine linear kernel expansion is a model of the type
57 /// \f[ x : \mathbb{R}^d \to \mathbb{R}^d \enspace , \enspace x \mapsto \sum_{n=1}^{\ell} \alpha_n k(x_n, x) + b \enspace ,\f]
58 /// with parameters \f$ \alpha_n \in \mathbb{R}^{d} \f$ for all
59 /// \f$ n \in \{1, \dots, \ell\} \f$ and \f$ b \in \mathbb{R}^d \f$.
60 ///
61 /// One way in which the possibility for vector-valued input and output of dimension \f$ d \f$ may be interpreted
62 /// is as allowing for a KernelExpansion model for \f$ d \f$ different classes/outputs in multi-class problems. Then,
63 /// the i-th column of the matrix #m_alpha is the KernelExpansion for class/output i.
64 ///
65 /// \tparam InputType Type of basis elements supplied to the kernel
66 ///
67 template<class InputType>
68 class KernelExpansion : public AbstractModel<InputType, RealVector>
69 {
70 public:
72  typedef AbstractModel<InputType, RealVector> base_type;
75 
76  // //////////////////////////////////////////////////////////
77  // //////////// CONSTRUCTORS /////////////////////
78  // //////////////////////////////////////////////////////////
79 
81 
82  KernelExpansion(KernelType* kernel):mep_kernel(kernel){
83  SHARK_ASSERT(kernel != NULL);
84  }
85 
86  KernelExpansion(KernelType* kernel, Data<InputType> const& basis,bool offset, std::size_t outputs = 1){
87  SHARK_ASSERT(kernel != NULL);
88  setStructure(kernel, basis,offset,outputs);
89  }
90 
91  void setStructure(KernelType* kernel, Data<InputType> const& basis,bool offset, std::size_t outputs = 1){
92  SHARK_ASSERT(kernel != NULL);
94  if(offset)
95  m_b.resize(outputs);
96  m_basis = basis;
97  m_alpha.resize(basis.numberOfElements(), outputs);
98  m_alpha.clear();
99  }
100 
101  /// \brief From INameable: return the class name.
102  std::string name() const
103  { return "KernelExpansion"; }
104 
105  /// dimensionality of the output RealVector
106  size_t outputSize() const{
107  return m_alpha.size2();
108  }
109 
110  // //////////////////////////////////////////////////////////
111  // /////////// ALL THINGS KERNEL //////////////////////
112  // //////////////////////////////////////////////////////////
113 
114  KernelType const* kernel() const{
115  return mep_kernel;
116  }
117  KernelType* kernel(){
118  return mep_kernel;
119  }
120  void setKernel(KernelType* kernel){
121  mep_kernel = kernel;
122  }
123 
124  // //////////////////////////////////////////////////////////
125  // /////// ALL THINGS ALPHA AND OFFSET ////////////////
126  // //////////////////////////////////////////////////////////
127 
128  bool hasOffset() const{
129  return m_b.size() != 0;
130  }
131  RealMatrix& alpha(){
132  return m_alpha;
133  }
134  RealMatrix const& alpha() const{
135  return m_alpha;
136  }
137  double& alpha(std::size_t example, std::size_t cls){
138  return m_alpha(example, cls);
139  }
140  double const& alpha(std::size_t example, std::size_t cls) const{
141  return m_alpha(example, cls);
142  }
143  RealVector& offset(){
144  SHARK_CHECK(hasOffset(), "[KernelExpansion::offset] invalid call for object without offset term");
145  return m_b;
146  }
147  RealVector const& offset() const{
148  SHARK_CHECK(hasOffset(), "[KernelExpansion::offset] invalid call for object without offset term");
149  return m_b;
150  }
151  double& offset(std::size_t cls){
152  SHARK_CHECK(hasOffset(), "[KernelExpansion::offset] invalid call for object without offset term");
153  return m_b(cls);
154  }
155  double const& offset(std::size_t cls) const{
156  SHARK_CHECK(hasOffset(), "[KernelExpansion::offset] invalid call for object without offset term");
157  return m_b(cls);
158  }
159 
160  // //////////////////////////////////////////////////////////
161  // //////// ALL THINGS UNDERLYING DATA ////////////////
162  // //////////////////////////////////////////////////////////
163 
164 
165  Data<InputType> const& basis() const {
166  return m_basis;
167  }
168 
170  return m_basis;
171  }
172 
173  /// The sparsify method removes non-support-vectors from
174  /// its set of basis vectors and the coefficient matrix.
175  void sparsify(){
176  std::size_t ic = m_basis.numberOfElements();
177  std::vector<std::size_t> svIndices;
178  for (std::size_t i=0; i != ic; ++i){
179  if (blas::norm_1(RealMatrixRow(m_alpha, i)) > 0.0){
180  svIndices.push_back(i);
181  }
182  }
183  //project basis on the support vectors
184  m_basis = toDataset(subset(toView(m_basis),svIndices));
185 
186  //reduce alpha to it's support vector variables
187  RealMatrix a(svIndices.size(), outputSize());
188  for (std::size_t i=0; i!= svIndices.size(); ++i){
189  noalias(row(a,i)) = row(m_alpha,svIndices[i]);
190  }
191  swap(m_alpha,a);
192  }
193 
194  // //////////////////////////////////////////////////////////
195  // //////// ALL THINGS KERNEL PARAMETERS //////////////
196  // //////////////////////////////////////////////////////////
197 
198  RealVector parameterVector() const{
199  RealVector ret(numberOfParameters());
200  if (hasOffset()){
201  init(ret) << toVector(m_alpha), m_b;
202  }
203  else{
204  init(ret) << toVector(m_alpha);
205  }
206  return ret;
207  }
208 
209  void setParameterVector(RealVector const& newParameters){
210  SHARK_CHECK(newParameters.size() == numberOfParameters(), "[KernelExpansion::setParameterVector] invalid size of the parameter vector");
211 
212  if (hasOffset())
213  init(newParameters) >> toVector(m_alpha), m_b;
214  else
215  init(newParameters) >> toVector(m_alpha);
216  }
217 
218  std::size_t numberOfParameters() const{
219  if (hasOffset())
220  return m_alpha.size1() * m_alpha.size2() + m_b.size();
221  else
222  return m_alpha.size1() * m_alpha.size2();
223  }
224 
225  // //////////////////////////////////////////////////////////
226  // //////// ALL THINGS EVALUATION //////////////
227  // //////////////////////////////////////////////////////////
228 
229  boost::shared_ptr<State> createState()const{
230  return boost::shared_ptr<State>(new EmptyState());
231  }
232 
234  void eval(BatchInputType const& patterns, BatchOutputType& output)const{
235  std::size_t numPatterns = boost::size(patterns);
236  SHARK_ASSERT(mep_kernel != NULL);
237 
238  output.resize(numPatterns,outputSize());
239  if (hasOffset())
240  output = repeat(m_b,numPatterns);
241  else
242  output.clear();
243 
244  std::size_t batchStart = 0;
245  for (std::size_t i=0; i != m_basis.numberOfBatches(); i++){
246  std::size_t batchEnd = batchStart+boost::size(m_basis.batch(i));
247  //evaluate kernels
248  //results in a matrix of the form where a column consists of the kernel evaluation of
249  //pattern i with respect to the batch of the basis,this gives a good memory alignment
250  //in the following matrix matrix product
251  RealMatrix kernelEvaluations = (*mep_kernel)(m_basis.batch(i),patterns);
252 
253  //get the part of the alpha matrix which is suitable for this batch
254  ConstRealSubMatrix batchAlpha = subrange(m_alpha,batchStart,batchEnd,0,outputSize());
255  noalias(output) += prod(trans(kernelEvaluations),batchAlpha);
256  batchStart = batchEnd;
257  }
258  }
259  void eval(BatchInputType const& patterns, BatchOutputType& outputs, State & state)const{
260  eval(patterns, outputs);
261  }
262 
263  // //////////////////////////////////////////////////////////
264  // //////// ALL THINGS SERIALIZATION //////////////
265  // //////////////////////////////////////////////////////////
266 
267  /// From ISerializable, reads a model from an archive
268  void read( InArchive & archive ){
269  SHARK_ASSERT(mep_kernel != NULL);
270 
271  archive >> m_alpha;
272  archive >> m_b;
273  archive >> m_basis;
274  archive >> (*mep_kernel);
275  }
276 
277  /// From ISerializable, writes a model to an archive
278  void write( OutArchive & archive ) const{
279  SHARK_ASSERT(mep_kernel != NULL);
280 
281  archive << m_alpha;
282  archive << m_b;
283  archive << m_basis;
284  archive << const_cast<KernelType const&>(*mep_kernel);//prevent compilation warning
285  }
286 
287 // //////////////////////////////////////////////////////////
288 // //////// MEMBERS //////////////
289 // //////////////////////////////////////////////////////////
290 
291 protected:
292  /// kernel function used in the expansion
293  KernelType* mep_kernel;
294 
295  /// "support" basis vectors
297 
298  /// kernel coefficients
299  RealMatrix m_alpha;
300 
301  /// offset or bias term
302  RealVector m_b;
303 };
304 
305 ///
306 /// \brief Linear classifier in a kernel feature space.
307 ///
308 /// This model is a simple wrapper for the KernelExpansion calculating the arg max
309 /// of the outputs of the model. This is the model used by kernel classifier models like SVMs.
310 ///
311 template<class InputType>
312 struct KernelClassifier: public ArgMaxConverter<KernelExpansion<InputType> >{
316 
318  { }
320  : base_type(KernelExpansionType(kernel))
321  { }
322  KernelClassifier(KernelExpansionType const& decisionFunction)
323  : base_type(decisionFunction)
324  { }
325 
326  std::string name() const
327  { return "KernelClassifier"; }
328 };
329 
330 
331 }
332 #endif