Metrics.h
Go to the documentation of this file.
1 /*!
2  *
3  *
4  * \brief Helper functions to calculate several norms and distances.
5  *
6  *
7  *
8  * \author O.Krause M.Thuma
9  * \date 2010-2011
10  *
11  *
12  * \par Copyright 1995-2017 Shark Development Team
13  *
14  * <BR><HR>
15  * This file is part of Shark.
16  * <http://shark-ml.org/>
17  *
18  * Shark is free software: you can redistribute it and/or modify
19  * it under the terms of the GNU Lesser General Public License as published
20  * by the Free Software Foundation, either version 3 of the License, or
21  * (at your option) any later version.
22  *
23  * Shark is distributed in the hope that it will be useful,
24  * but WITHOUT ANY WARRANTY; without even the implied warranty of
25  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26  * GNU Lesser General Public License for more details.
27  *
28  * You should have received a copy of the GNU Lesser General Public License
29  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
30  *
31  */
32 #ifndef SHARK_LINALG_METRICS_H
33 #define SHARK_LINALG_METRICS_H
34 
35 #include <shark/LinAlg/BLAS/remora.hpp>
36 #include <shark/Core/Math.h>
37 namespace remora{
38 
39 ///////////////////////////////////////NORMS////////////////////////////////////////
40 
41 /**
42 * \brief Normalized squared norm_2 (diagonal Mahalanobis).
43 *
44 * Contrary to some conventions, dimension-wise weights are considered instead of std. deviations:
45 * \f$ n^2(v) = \sum_i w_i v_i^2 \f$
46 * nb: the weights themselves are not squared, but multiplied onto the squared components
47 */
48 template<class VectorT, class WeightT, class Device>
49 typename VectorT::value_type diagonalMahalanobisNormSqr(
50  vector_expression<VectorT, Device> const& vector,
51  vector_expression<WeightT, Device> const& weights
52 ) {
53  SIZE_CHECK( vector().size() == weights().size() );
54  return inner_prod(weights(),sqr(vector()));
55 }
56 
57 /**
58 * \brief Normalized norm_2 (diagonal Mahalanobis).
59 *
60 * Contrary to some conventions, dimension-wise weights are considered instead of std. deviations:
61 * \f$ n^2(v) = \sqrt{\sum_i w_i v_i^2} \f$
62 * nb: the weights themselves are not squared, but multiplied onto the squared components
63 */
64 template<class VectorT, class WeightT, class Device>
65 typename VectorT::value_type diagonalMahalanobisNorm(
66  vector_expression<VectorT, Device> const& vector,
67  vector_expression<WeightT, Device> const& weights
68 ) {
69  SIZE_CHECK( vector().size() == weights().size() );
70  return std::sqrt(diagonalMahalanobisNormSqr(vector,weights));
71 }
72 
73 ////////////////////////////////////////DISTANCES/////////////////////////////////////////////////
74 
75 namespace detail{
76  /**
77  * \brief Normalized Euclidian squared distance (squared diagonal Mahalanobis)
78  * between two vectors, optimized for two Compressed arguments.
79  *
80  * Contrary to some conventions, dimension-wise weights are considered instead of std. deviations:
81  * \f$ d^2(v) = \sum_i w_i (x_i-z_i)^2 \f$
82  * NOTE: The weights themselves are not squared, but multiplied onto the squared components.
83  */
84  template<class VectorT, class VectorU, class WeightT>
85  typename VectorT::value_type diagonalMahalanobisDistanceSqr(
86  VectorT const& op1,
87  VectorU const& op2,
88  WeightT const& weights,
89  sparse_tag,
90  sparse_tag
91  ){
92  using shark::sqr;
93  typename VectorT::value_type sum=0;
94  typename VectorT::const_iterator iter1=op1.begin();
95  typename VectorU::const_iterator iter2=op2.begin();
96  typename VectorT::const_iterator end1=op1.end();
97  typename VectorU::const_iterator end2=op2.end();
98  //be aware of empty vectors!
99  while(iter1 != end1 && iter2 != end2)
100  {
101  std::size_t index1=iter1.index();
102  std::size_t index2=iter2.index();
103  if(index1==index2){
104  sum += weights(index1) * sqr(*iter1-*iter2);
105  ++iter1;
106  ++iter2;
107  }
108  else if(index1<index2){
109  sum += weights(index1) * sqr(*iter1);
110  ++iter1;
111  }
112  else {
113  sum += weights(index2) * sqr(*iter2);
114  ++iter2;
115  }
116  }
117  while(iter1 != end1)
118  {
119  std::size_t index1=iter1.index();
120  sum += weights(index1) * sqr(*iter1);
121  ++iter1;
122  }
123  while(iter2 != end2)
124  {
125  std::size_t index2=iter2.index();
126  sum += weights(index2) * sqr(*iter2);
127  ++iter2;
128  }
129  return sum;
130  }
131 
132  /**
133  * \brief Normalized Euclidian squared distance (squared diagonal Mahalanobis)
134  * between two vectors, optimized for one dense and one sparse argument
135  */
136  template<class VectorT, class VectorU, class WeightT>
137  typename VectorT::value_type diagonalMahalanobisDistanceSqr(
138  VectorT const& op1,
139  VectorU const& op2,
140  WeightT const& weights,
141  sparse_tag,
142  dense_tag
143  ){
144  using shark::sqr;
145  typename VectorT::const_iterator iter=op1.begin();
146  typename VectorT::const_iterator end=op1.end();
147 
148  std::size_t index = 0;
149  std::size_t pos = 0;
150  typename VectorT::value_type sum=0;
151 
152  for(;iter != end;++iter,++pos){
153  index = iter.index();
154  for(;pos != index;++pos){
155  sum += weights(pos) * sqr(op2(pos));
156  }
157  sum += weights(index) * sqr(*iter-op2(pos));
158  }
159  for(;pos != op2.size();++pos){
160  sum += weights(pos) * sqr(op2(pos));
161  }
162  return sum;
163  }
164  template<class VectorT, class VectorU, class WeightT>
165  typename VectorT::value_type diagonalMahalanobisDistanceSqr(
166  VectorT const& op1,
167  VectorU const& op2,
168  WeightT const& weights,
169  dense_tag arg1tag,
170  sparse_tag arg2tag
171  ){
172  return diagonalMahalanobisDistanceSqr(op2,op1,weights,arg2tag,arg1tag);
173  }
174 
175  template<class VectorT, class VectorU, class WeightT>
176  typename VectorT::value_type diagonalMahalanobisDistanceSqr(
177  VectorT const& op1,
178  VectorU const& op2,
179  WeightT const& weights,
180  dense_tag,
181  dense_tag
182  ){
183  return inner_prod(op1-op2,(op1-op2)*weights);
184  }
185 
186 
187  template<class MatrixT,class VectorU, class Result>
188  void distanceSqrBlockVector(
189  MatrixT const& operands,
190  VectorU const& op2,
191  Result& result
192  ){
193  typedef typename Result::value_type value_type;
194  scalar_vector< value_type, cpu_tag > one(op2.size(),static_cast<value_type>(1.0));
195  for(std::size_t i = 0; i != operands.size1(); ++i){
196  result(i) = diagonalMahalanobisDistanceSqr(
197  row(operands,i),op2,one,
198  typename MatrixT::evaluation_category::tag(),
199  typename VectorU::evaluation_category::tag()
200  );
201  }
202  }
203 
204  ///\brief implementation for two input blocks where at least one matrix has only a few rows
205  template<class MatrixX,class MatrixY, class Result>
206  void distanceSqrBlockBlockRowWise(
207  MatrixX const& X,
208  MatrixY const& Y,
209  Result& distances
210  ){
211  std::size_t sizeX=X.size1();
212  std::size_t sizeY=Y.size1();
213  if(sizeX < sizeY){//iterate over the rows of the block with less rows
214  for(std::size_t i = 0; i != sizeX; ++i){
215  matrix_row<Result> distanceRow = row(distances,i);
216  distanceSqrBlockVector(
217  Y,row(X,i),distanceRow
218  );
219  }
220  }else{
221  for(std::size_t i = 0; i != sizeY; ++i){
222  auto distanceCol = column(distances,i);
223  distanceSqrBlockVector(
224  X,row(Y,i),distanceCol
225  );
226  }
227  }
228  }
229 
230  ///\brief implementation for two dense input blocks
231  template<class MatrixX,class MatrixY, class Result>
232  void distanceSqrBlockBlock(
233  MatrixX const& X,
234  MatrixY const& Y,
235  Result& distances,
236  dense_tag,
237  dense_tag
238  ){
239  typedef typename Result::value_type value_type;
240  std::size_t sizeX=X.size1();
241  std::size_t sizeY=Y.size1();
242  ensure_size(distances,X.size1(),Y.size1());
243  if(sizeX < 10 || sizeY<10){
244  distanceSqrBlockBlockRowWise(X,Y,distances);
245  return;
246  }
247  //fast blockwise iteration
248  //uses: (a-b)^2 = a^2 -2ab +b^2
249  noalias(distances) = -2*prod(X,trans(Y));
250  //first a^2+b^2
251  vector<value_type> ySqr(sizeY);
252  for(std::size_t i = 0; i != sizeY; ++i){
253  ySqr(i) = norm_sqr(row(Y,i));
254  }
255  //initialize d_ij=x_i^2+y_i^2
256  for(std::size_t i = 0; i != sizeX; ++i){
257  value_type xSqr = norm_sqr(row(X,i));
258  noalias(row(distances,i)) += repeat(xSqr,sizeY) + ySqr;
259  }
260  }
261  //\brief default implementation used, when one of the arguments is not dense
262  template<class MatrixX,class MatrixY,class Result>
263  void distanceSqrBlockBlock(
264  MatrixX const& X,
265  MatrixY const& Y,
266  Result& distances,
267  sparse_tag,
268  sparse_tag
269  ){
270  distanceSqrBlockBlockRowWise(X,Y,distances);
271  }
272 }
273 
274 /**
275 * \ingroup shark_globals
276 *
277 * @{
278 */
279 
280 /**
281 * \brief Normalized Euclidian squared distance (squared diagonal Mahalanobis)
282 * between two vectors.
283 *
284 * NOTE: The weights themselves are not squared, but multiplied onto the squared components.
285 */
286 template<class VectorT, class VectorU, class WeightT, class Device>
287 typename VectorT::value_type diagonalMahalanobisDistanceSqr(
288  vector_expression<VectorT, Device> const& op1,
289  vector_expression<VectorU, Device> const& op2,
290  vector_expression<WeightT, Device> const& weights
291 ){
292  SIZE_CHECK(op1().size()==op2().size());
293  SIZE_CHECK(op1().size()==weights().size());
294  //dispatch given the types of the argument
295  return detail::diagonalMahalanobisDistanceSqr(
296  op1(), op2(), weights(),
297  typename VectorT::evaluation_category::tag(),
298  typename VectorU::evaluation_category::tag()
299  );
300 }
301 
302 /**
303 * \brief Squared distance between two vectors.
304 */
305 template<class VectorT,class VectorU, class Device>
306 typename VectorT::value_type distanceSqr(
307  vector_expression<VectorT, Device> const& op1,
308  vector_expression<VectorU, Device> const& op2
309 ){
310  SIZE_CHECK(op1().size()==op2().size());
311  typedef typename VectorT::value_type value_type;
312  scalar_vector< value_type, cpu_tag > one(op1().size(),static_cast<value_type>(1.0));
313  return diagonalMahalanobisDistanceSqr(op1,op2,one);
314 }
315 
316 /**
317 * \brief Squared distance between a vector and a set of vectors and stores the result in the vector of distances
318 *
319 * The squared distance between the vector and every row-vector of the matrix is calculated.
320 * This can be implemented much more efficiently.
321 */
322 template<class MatrixT,class VectorU, class VectorR, class Device>
323 void distanceSqr(
324  matrix_expression<MatrixT, Device> const& operands,
325  vector_expression<VectorU, Device> const& op2,
326  vector_expression<VectorR, Device>& distances
327 ){
328  SIZE_CHECK(operands().size2()==op2().size());
329  ensure_size(distances,operands().size1());
330  detail::distanceSqrBlockVector(
331  operands(),op2(),distances()
332  );
333 }
334 
335 /**
336 * \brief Squared distance between a vector and a set of vectors
337 *
338 * The squared distance between the vector and every row-vector of the matrix is calculated.
339 * This can be implemented much more efficiently.
340 */
341 template<class MatrixT,class VectorU, class Device>
342 vector<typename MatrixT::value_type> distanceSqr(
343  matrix_expression<MatrixT, Device> const& operands,
344  vector_expression<VectorU, Device> const& op2
345 ){
346  SIZE_CHECK(operands().size2()==op2().size());
347  vector<typename MatrixT::value_type> distances(operands().size1());
348  distanceSqr(operands,op2,distances);
349  return distances;
350 }
351 
352 /**
353 * \brief Squared distance between a vector and a set of vectors
354 *
355 * The squared distance between the vector and every row-vector of the matrix is calculated.
356 * This can be implemented much more efficiently.
357 */
358 template<class MatrixT,class VectorU, class Device>
359 vector<typename MatrixT::value_type> distanceSqr(
360  vector_expression<VectorU, Device> const& op1,
361  matrix_expression<MatrixT, Device> const& operands
362 ){
363  SIZE_CHECK(operands().size2()==op1().size());
364  vector<typename MatrixT::value_type> distances(operands().size1());
365  distanceSqr(operands,op1,distances);
366  return distances;
367 }
368 
369 /**
370 * \brief Squared distance between the vectors of two sets of vectors
371 *
372 * The squared distance between every row-vector of the first matrix x
373 * and every row-vector of the second matrix y is calculated.
374 * This can be implemented much more efficiently.
375 * The results are returned as a matrix, where the element in the i-th
376 * row and the j-th column is distanceSqr(x_i,y_j).
377 */
378 template<class MatrixT,class MatrixU, class Device>
379 matrix<typename MatrixT::value_type> distanceSqr(
380  matrix_expression<MatrixT, Device> const& X,
381  matrix_expression<MatrixU, Device> const& Y
382 ){
383  typedef matrix<typename MatrixT::value_type> Matrix;
384  SIZE_CHECK(X().size2()==Y().size2());
385  std::size_t sizeX=X().size1();
386  std::size_t sizeY=Y().size1();
387  Matrix distances(sizeX, sizeY);
388  detail::distanceSqrBlockBlock(
389  X(),Y(),distances,
390  typename MatrixT::evaluation_category::tag(),
391  typename MatrixU::evaluation_category::tag()
392  );
393  return distances;
394 
395 }
396 
397 
398 /**
399 * \brief Calculates distance between two vectors.
400 */
401 template<class VectorT,class VectorU, class Device>
402 typename VectorT::value_type distance(
403  vector_expression<VectorT, Device> const& op1,
404  vector_expression<VectorU, Device> const& op2
405 ){
406  SIZE_CHECK(op1().size()==op2().size());
407  return std::sqrt(distanceSqr(op1,op2));
408 }
409 
410 /**
411 * \brief Normalized euclidian distance (diagonal Mahalanobis) between two vectors.
412 *
413 * Contrary to some conventions, dimension-wise weights are considered instead of std. deviations:
414 * \f$ d(v) = \left( \sum_i w_i (x_i-z_i)^2 \right)^{1/2} \f$
415 * nb: the weights themselves are not squared, but multiplied onto the squared components
416 */
417 template<class VectorT, class VectorU, class WeightT, class Device>
418 typename VectorT::value_type diagonalMahalanobisDistance(
419  vector_expression<VectorT, Device> const& op1,
420  vector_expression<VectorU, Device> const& op2,
421  vector_expression<WeightT, Device> const& weights
422 ){
423  SIZE_CHECK(op1().size()==op2().size());
424  SIZE_CHECK(op1().size()==weights().size());
425  return std::sqrt(diagonalMahalanobisDistanceSqr(op1(), op2(), weights));
426 }
427 /** @}*/
428 }
429 
430 #endif