ViennaCL - The Vienna Computing Library  1.5.2
nmf.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_NMF_HPP
2 #define VIENNACL_LINALG_NMF_HPP
3 
4 /* =========================================================================
5  Copyright (c) 2010-2014, 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 PDF manual)
17 
18  License: MIT (X11), see file LICENSE in the base directory
19 ============================================================================= */
20 
28 #include "viennacl/vector.hpp"
29 #include "viennacl/matrix.hpp"
30 #include "viennacl/linalg/prod.hpp"
34 
35 namespace viennacl
36 {
37  namespace linalg
38  {
40  class nmf_config
41  {
42  public:
43  nmf_config(double val_epsilon = 1e-4,
44  double val_epsilon_stagnation = 1e-5,
45  vcl_size_t num_max_iters = 10000,
46  vcl_size_t num_check_iters = 100)
47  : eps_(val_epsilon), stagnation_eps_(val_epsilon_stagnation),
48  max_iters_(num_max_iters),
49  check_after_steps_( (num_check_iters > 0) ? num_check_iters : 1),
50  print_relative_error_(false),
51  iters_(0) {}
52 
54  double tolerance() const { return eps_; }
55 
57  void tolerance(double e) { eps_ = e; }
58 
60  double stagnation_tolerance() const { return stagnation_eps_; }
61 
63  void stagnation_tolerance(double e) { stagnation_eps_ = e; }
64 
66  vcl_size_t max_iterations() const { return max_iters_; }
68  void max_iterations(vcl_size_t m) { max_iters_ = m; }
69 
71  vcl_size_t iters() const { return iters_; }
72 
73 
75  vcl_size_t check_after_steps() const { return check_after_steps_; }
77  void check_after_steps(vcl_size_t c) { if (c > 0) check_after_steps_ = c; }
78 
80  bool print_relative_error() const { return print_relative_error_; }
82  void print_relative_error(bool b) { print_relative_error_ = b; }
83 
84  template <typename ScalarType>
85  friend void nmf(viennacl::matrix<ScalarType> const & V,
88  nmf_config const & conf);
89 
90  private:
91  double eps_;
92  double stagnation_eps_;
93  vcl_size_t max_iters_;
94  vcl_size_t check_after_steps_;
95  bool print_relative_error_;
96  mutable vcl_size_t iters_;
97  };
98 
99 
107  template <typename ScalarType>
111  nmf_config const & conf)
112  {
113  viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(V).context());
114 
115  const std::string NMF_MUL_DIV_KERNEL = "el_wise_mul_div";
116 
118 
119  assert(V.size1() == W.size1() && V.size2() == H.size2() && bool("Dimensions of W and H don't allow for V = W * H"));
120  assert(W.size2() == H.size1() && bool("Dimensions of W and H don't match, prod(W, H) impossible"));
121 
122  vcl_size_t k = W.size2();
123  conf.iters_ = 0;
124 
128 
131  viennacl::matrix<ScalarType> htmp(k, k);
132 
134  viennacl::vector<ScalarType> diff(V.size1() * V.size2());
135 
136  ScalarType last_diff = 0;
137  ScalarType diff_init = 0;
138  bool stagnation_flag = false;
139 
140 
141  for (vcl_size_t i = 0; i < conf.max_iterations(); i++)
142  {
143  conf.iters_ = i + 1;
144  {
145  hn = viennacl::linalg::prod(trans(W), V);
146  htmp = viennacl::linalg::prod(trans(W), W);
147  hd = viennacl::linalg::prod(htmp, H);
148 
150  viennacl::ocl::enqueue(mul_div_kernel(H, hn, hd, cl_uint(H.internal_size1() * H.internal_size2())));
151  }
152  {
153  wn = viennacl::linalg::prod(V, trans(H));
154  wtmp = viennacl::linalg::prod(W, H);
155  wd = viennacl::linalg::prod(wtmp, trans(H));
156 
158 
159  viennacl::ocl::enqueue(mul_div_kernel(W, wn, wd, cl_uint(W.internal_size1() * W.internal_size2())));
160  }
161 
162  if (i % conf.check_after_steps() == 0) //check for convergence
163  {
164  appr = viennacl::linalg::prod(W, H);
165 
166  appr -= V;
167  ScalarType diff_val = viennacl::linalg::norm_frobenius(appr);
168 
169  if (i == 0)
170  diff_init = diff_val;
171 
172  if (conf.print_relative_error())
173  std::cout << diff_val / diff_init << std::endl;
174 
175  // Approximation check
176  if (diff_val / diff_init < conf.tolerance())
177  break;
178 
179  // Stagnation check
180  if (std::fabs(diff_val - last_diff) / (diff_val * conf.check_after_steps()) < conf.stagnation_tolerance()) //avoid situations where convergence stagnates
181  {
182  if (stagnation_flag) // iteration stagnates (two iterates with no notable progress)
183  break;
184  else // record stagnation in this iteration
185  stagnation_flag = true;
186  }
187  else // good progress in this iteration, so unset stagnation flag
188  stagnation_flag = false;
189 
190  // prepare for next iterate:
191  last_diff = diff_val;
192  }
193  }
194 
195 
196  }
197  }
198 }
199 
200 #endif
viennacl::ocl::kernel & get_kernel(std::string const &program_name, std::string const &kernel_name)
Convenience function for retrieving the kernel of a program directly from the context.
Definition: context.hpp:470
std::size_t vcl_size_t
Definition: forwards.h:58
Generic interface for the l^2-norm. See viennacl/linalg/vector_operations.hpp for implementations...
viennacl::enable_if< viennacl::is_any_sparse_matrix< M1 >::value, matrix_expression< const M1, const M1, op_trans > >::type trans(const M1 &mat)
Returns an expression template class representing a transposed matrix.
Definition: sparse_matrix_operations.hpp:330
double stagnation_tolerance() const
Relative tolerance for the stagnation check.
Definition: nmf.hpp:60
vcl_size_t check_after_steps() const
Number of steps after which the convergence of NMF should be checked (again)
Definition: nmf.hpp:75
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
size_type size2() const
Returns the number of columns.
Definition: matrix.hpp:627
Represents an OpenCL kernel within ViennaCL.
Definition: kernel.hpp:59
Implementation of the dense matrix class.
Manages an OpenCL context and provides the respective convenience functions for creating buffers...
Definition: context.hpp:51
double tolerance() const
Returns the relative tolerance for convergence.
Definition: nmf.hpp:54
bool print_relative_error() const
Returns the flag specifying whether the relative tolerance should be printed in each iteration...
Definition: nmf.hpp:80
Configuration class for the nonnegative-matrix-factorization algorithm. Specify tolerances, maximum iteration counts, etc., here.
Definition: nmf.hpp:40
A dense matrix class.
Definition: forwards.h:293
void nmf(viennacl::matrix< ScalarType > const &V, viennacl::matrix< ScalarType > &W, viennacl::matrix< ScalarType > &H, nmf_config const &conf)
The nonnegative matrix factorization (approximation) algorithm as suggested by Lee and Seung...
Definition: nmf.hpp:108
void tolerance(double e)
Sets the relative tolerance for convergence, i.e. norm(V - W * H) / norm(V - W_init * H_init) ...
Definition: nmf.hpp:57
size_type internal_size2() const
Returns the internal number of columns. Usually required for launching OpenCL kernels only...
Definition: matrix.hpp:649
Main namespace in ViennaCL. Holds all the basic types such as vector, matrix, etc. and defines operations upon them.
Definition: cpu_ram.hpp:29
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:91
void enqueue(KernelType &k, viennacl::ocl::command_queue const &queue)
Enqueues a kernel in the provided queue.
Definition: enqueue.hpp:48
void print_relative_error(bool b)
Specify whether the relative error should be printed at each convergence check after 'num_check_iters...
Definition: nmf.hpp:82
void max_iterations(vcl_size_t m)
Sets the maximum number of iterations for the NMF algorithm.
Definition: nmf.hpp:68
Generic interface for the Frobenius norm.
vcl_size_t max_iterations() const
Returns the maximum number of iterations for the NMF algorithm.
Definition: nmf.hpp:66
Main kernel class for generating OpenCL kernels for nonnegative matrix factorization of a dense matri...
Definition: nmf.hpp:41
vcl_size_t iters() const
Returns the number of iterations of the last NMF run using this configuration object.
Definition: nmf.hpp:71
friend void nmf(viennacl::matrix< ScalarType > const &V, viennacl::matrix< ScalarType > &W, viennacl::matrix< ScalarType > &H, nmf_config const &conf)
The nonnegative matrix factorization (approximation) algorithm as suggested by Lee and Seung...
Definition: nmf.hpp:108
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)
Definition: nmf.hpp:43
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
void stagnation_tolerance(double e)
Sets the tolerance for the stagnation check (i.e. the minimum required relative change of the residua...
Definition: nmf.hpp:63
scalar_expression< const matrix_base< NumericT, F >, const matrix_base< NumericT, F >, op_norm_frobenius > norm_frobenius(const matrix< NumericT, F > &A)
Definition: norm_frobenius.hpp:61
OpenCL kernel file for nonnegative matrix factorization.
void check_after_steps(vcl_size_t c)
Set the number of steps after which the convergence of NMF should be checked (again) ...
Definition: nmf.hpp:77
static void init(viennacl::ocl::context &ctx)
Definition: nmf.hpp:48
size_type internal_size1() const
Returns the internal number of rows. Usually required for launching OpenCL kernels only...
Definition: matrix.hpp:647
size_type size1() const
Returns the number of rows.
Definition: matrix.hpp:625