ViennaCL - The Vienna Computing Library  1.7.0
Free open-source GPU-accelerated linear algebra and solver library.
nmf_operations.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_HOST_BASED_NMF_OPERATIONS_HPP_
2 #define VIENNACL_LINALG_HOST_BASED_NMF_OPERATIONS_HPP_
3 
4 /* =========================================================================
5  Copyright (c) 2010-2015, Institute for Microelectronics,
6  Institute for Analysis and Scientific Computing,
7  TU Wien.
8  Portions of this software are copyright by UChicago Argonne, LLC.
9 
10  -----------------
11  ViennaCL - The Vienna Computing Library
12  -----------------
13 
14  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
15 
16  (A list of authors and contributors can be found in the manual)
17 
18  License: MIT (X11), see file LICENSE in the base directory
19  ============================================================================= */
20 
25 #include "viennacl/vector.hpp"
26 #include "viennacl/matrix.hpp"
27 #include "viennacl/linalg/prod.hpp"
30 
32 
33 namespace viennacl
34 {
35 namespace linalg
36 {
37 
40 {
41 public:
42  nmf_config(double val_epsilon = 1e-4, double val_epsilon_stagnation = 1e-5,
43  vcl_size_t num_max_iters = 10000, vcl_size_t num_check_iters = 100) :
44  eps_(val_epsilon), stagnation_eps_(val_epsilon_stagnation), max_iters_(num_max_iters), check_after_steps_(
45  (num_check_iters > 0) ? num_check_iters : 1), print_relative_error_(false), iters_(0)
46  {
47  }
48 
50  double tolerance() const
51  {
52  return eps_;
53  }
54 
56  void tolerance(double e)
57  {
58  eps_ = e;
59  }
60 
62  double stagnation_tolerance() const
63  {
64  return stagnation_eps_;
65  }
66 
68  void stagnation_tolerance(double e)
69  {
70  stagnation_eps_ = e;
71  }
72 
75  {
76  return max_iters_;
77  }
80  {
81  max_iters_ = m;
82  }
83 
85  vcl_size_t iters() const
86  {
87  return iters_;
88  }
89 
92  {
93  return check_after_steps_;
94  }
97  {
98  if (c > 0)
99  check_after_steps_ = c;
100  }
101 
103  bool print_relative_error() const
104  {
105  return print_relative_error_;
106  }
108  void print_relative_error(bool b)
109  {
110  print_relative_error_ = b;
111  }
112 
113  template<typename ScalarType>
114  friend void nmf(viennacl::matrix_base<ScalarType> const & V,
116  nmf_config const & conf);
117 
118 private:
119  double eps_;
120  double stagnation_eps_;
121  vcl_size_t max_iters_;
122  vcl_size_t check_after_steps_;
123  bool print_relative_error_;
124 public:
126 };
127 
128 namespace host_based
129 {
131  template<typename NumericT>
132  void el_wise_mul_div(NumericT * matrix1,
133  NumericT const * matrix2,
134  NumericT const * matrix3, vcl_size_t size)
135  {
136 #ifdef VIENNACL_WITH_OPENMP
137 #pragma omp parallel for
138 #endif
139  for (long i2 = 0; i2 < long(size); i2++)
140  {
141  vcl_size_t i = vcl_size_t(i2);
142  NumericT val = matrix1[i] * matrix2[i];
143  NumericT divisor = matrix3[i];
144  matrix1[i] = (divisor > (NumericT) 0.00001) ? (val / divisor) : (NumericT) 0;
145  }
146  }
147 
155  template<typename NumericT>
159  viennacl::linalg::nmf_config const & conf)
160  {
161  vcl_size_t k = W.size2();
162  conf.iters_ = 0;
163 
166 
169 
173 
177 
179 
180  NumericT last_diff = 0;
181  NumericT diff_init = 0;
182  bool stagnation_flag = false;
183 
184  for (vcl_size_t i = 0; i < conf.max_iterations(); i++)
185  {
186  conf.iters_ = i + 1;
187 
188  hn = viennacl::linalg::prod(trans(W), V);
189  htmp = viennacl::linalg::prod(trans(W), W);
190  hd = viennacl::linalg::prod(htmp, H);
191 
192  NumericT * data_H = detail::extract_raw_pointer<NumericT>(H);
193  NumericT * data_hn = detail::extract_raw_pointer<NumericT>(hn);
194  NumericT * data_hd = detail::extract_raw_pointer<NumericT>(hd);
195 
197 
198  wn = viennacl::linalg::prod(V, trans(H));
199  wtmp = viennacl::linalg::prod(W, H);
200  wd = viennacl::linalg::prod(wtmp, trans(H));
201 
202  NumericT * data_W = detail::extract_raw_pointer<NumericT>(W);
203  NumericT * data_wn = detail::extract_raw_pointer<NumericT>(wn);
204  NumericT * data_wd = detail::extract_raw_pointer<NumericT>(wd);
205 
207 
208  if (i % conf.check_after_steps() == 0) //check for convergence
209  {
210  appr = viennacl::linalg::prod(W, H);
211 
212  appr -= V;
214 
215  if (i == 0)
216  diff_init = diff_val;
217 
218  if (conf.print_relative_error())
219  std::cout << diff_val / diff_init << std::endl;
220 
221  // Approximation check
222  if (diff_val / diff_init < conf.tolerance())
223  break;
224 
225  // Stagnation check
226  if (std::fabs(diff_val - last_diff) / (diff_val * NumericT(conf.check_after_steps())) < conf.stagnation_tolerance()) //avoid situations where convergence stagnates
227  {
228  if (stagnation_flag) // iteration stagnates (two iterates with no notable progress)
229  break;
230  else
231  // record stagnation in this iteration
232  stagnation_flag = true;
233  } else
234  // good progress in this iteration, so unset stagnation flag
235  stagnation_flag = false;
236 
237  // prepare for next iterate:
238  last_diff = diff_val;
239  }
240  }
241  }
242 
243 } //namespace host_based
244 } //namespace linalg
245 } //namespace viennacl
246 
247 #endif /* VIENNACL_LINALG_HOST_BASED_NMF_OPERATIONS_HPP_ */
Generic interface for the l^2-norm. See viennacl/linalg/vector_operations.hpp for implementations...
friend void nmf(viennacl::matrix_base< ScalarType > const &V, viennacl::matrix_base< ScalarType > &W, viennacl::matrix_base< ScalarType > &H, nmf_config const &conf)
The nonnegative matrix factorization (approximation) algorithm as suggested by Lee and Seung...
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
Implementation of the dense matrix class.
vcl_size_t check_after_steps() const
Number of steps after which the convergence of NMF should be checked (again)
Configuration class for the nonnegative-matrix-factorization algorithm. Specify tolerances, maximum iteration counts, etc., here.
void nmf(viennacl::matrix_base< NumericT > const &V, viennacl::matrix_base< NumericT > &W, viennacl::matrix_base< NumericT > &H, viennacl::linalg::nmf_config const &conf)
The nonnegative matrix factorization (approximation) algorithm as suggested by Lee and Seung...
float NumericT
Definition: bisect.cpp:40
Main namespace in ViennaCL. Holds all the basic types such as vector, matrix, etc. and defines operations upon them.
Definition: cpu_ram.hpp:34
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:102
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:235
Generic interface for the Frobenius norm.
vcl_size_t max_iterations() const
Returns the maximum number of iterations for the NMF algorithm.
void check_after_steps(vcl_size_t c)
Set the number of steps after which the convergence of NMF should be checked (again) ...
Represents a vector consisting of scalars 's' only, i.e. v[i] = s for all i. To be used as an initial...
Definition: matrix_def.hpp:93
double tolerance() const
Returns the relative tolerance for convergence.
void print_relative_error(bool b)
Specify whether the relative error should be printed at each convergence check after 'num_check_iters...
std::size_t vcl_size_t
Definition: forwards.h:75
size_type size2() const
Returns the number of columns.
Definition: matrix_def.hpp:226
Common routines for single-threaded or OpenMP-enabled execution on CPU.
void tolerance(double e)
Sets the relative tolerance for convergence, i.e. norm(V - W * H) / norm(V - W_init * H_init) ...
void el_wise_mul_div(NumericT *matrix1, NumericT const *matrix2, NumericT const *matrix3, vcl_size_t size)
Missing OpenMP kernel for nonnegative matrix factorization of a dense matrices.
size_type size1() const
Returns the number of rows.
Definition: matrix_def.hpp:224
double stagnation_tolerance() const
Relative tolerance for the stagnation check.
nmf_config(double val_epsilon=1e-4, double val_epsilon_stagnation=1e-5, vcl_size_t num_max_iters=10000, vcl_size_t num_check_iters=100)
vcl_size_t iters() const
Returns the number of iterations of the last NMF run using this configuration object.
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
bool row_major() const
Definition: matrix_def.hpp:248
scalar_expression< const matrix_base< NumericT >, const matrix_base< NumericT >, op_norm_frobenius > norm_frobenius(const matrix_base< NumericT > &A)
void trans(const matrix_expression< const matrix_base< NumericT, SizeT, DistanceT >, const matrix_base< NumericT, SizeT, DistanceT >, op_trans > &proxy, matrix_base< NumericT > &temp_trans)
void max_iterations(vcl_size_t m)
Sets the maximum number of iterations for the NMF algorithm.
void stagnation_tolerance(double e)
Sets the tolerance for the stagnation check (i.e. the minimum required relative change of the residua...
size_type internal_size2() const
Returns the internal number of columns. Usually required for launching OpenCL kernels only...
Definition: matrix_def.hpp:240
size_type internal_size1() const
Returns the internal number of rows. Usually required for launching OpenCL kernels only...
Definition: matrix_def.hpp:238
bool print_relative_error() const
Returns the flag specifying whether the relative tolerance should be printed in each iteration...