1 #ifndef VIENNACL_LINALG_HOST_BASED_NMF_OPERATIONS_HPP_ 2 #define VIENNACL_LINALG_HOST_BASED_NMF_OPERATIONS_HPP_ 3 4 /* ========================================================================= 5 Copyright (c) 2010-2016, 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 21 /** @file viennacl/linalg/host_based/vector_operations.hpp 22 @brief Implementations of NMF operations using a plain single-threaded or OpenMP-enabled execution on CPU 23 */ 24 25 #include "viennacl/vector.hpp" 26 #include "viennacl/matrix.hpp" 27 #include "viennacl/linalg/prod.hpp" 28 #include "viennacl/linalg/norm_2.hpp" 29 #include "viennacl/linalg/norm_frobenius.hpp" 30 31 #include "viennacl/linalg/host_based/common.hpp" 32 33 namespace viennacl 34 { 35 namespace linalg 36 { 37 38 /** @brief Configuration class for the nonnegative-matrix-factorization algorithm. Specify tolerances, maximum iteration counts, etc., here. */ 39 class nmf_config 40 { 41 public: 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)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 49 /** @brief Returns the relative tolerance for convergence */ tolerance() const50 double tolerance() const 51 { 52 return eps_; 53 } 54 55 /** @brief Sets the relative tolerance for convergence, i.e. norm(V - W * H) / norm(V - W_init * H_init) */ tolerance(double e)56 void tolerance(double e) 57 { 58 eps_ = e; 59 } 60 61 /** @brief Relative tolerance for the stagnation check */ stagnation_tolerance() const62 double stagnation_tolerance() const 63 { 64 return stagnation_eps_; 65 } 66 67 /** @brief Sets the tolerance for the stagnation check (i.e. the minimum required relative change of the residual between two iterations) */ stagnation_tolerance(double e)68 void stagnation_tolerance(double e) 69 { 70 stagnation_eps_ = e; 71 } 72 73 /** @brief Returns the maximum number of iterations for the NMF algorithm */ max_iterations() const74 vcl_size_t max_iterations() const 75 { 76 return max_iters_; 77 } 78 /** @brief Sets the maximum number of iterations for the NMF algorithm */ max_iterations(vcl_size_t m)79 void max_iterations(vcl_size_t m) 80 { 81 max_iters_ = m; 82 } 83 84 /** @brief Returns the number of iterations of the last NMF run using this configuration object */ iters() const85 vcl_size_t iters() const 86 { 87 return iters_; 88 } 89 90 /** @brief Number of steps after which the convergence of NMF should be checked (again) */ check_after_steps() const91 vcl_size_t check_after_steps() const 92 { 93 return check_after_steps_; 94 } 95 /** @brief Set the number of steps after which the convergence of NMF should be checked (again) */ check_after_steps(vcl_size_t c)96 void check_after_steps(vcl_size_t c) 97 { 98 if (c > 0) 99 check_after_steps_ = c; 100 } 101 102 /** @brief Returns the flag specifying whether the relative tolerance should be printed in each iteration */ print_relative_error() const103 bool print_relative_error() const 104 { 105 return print_relative_error_; 106 } 107 /** @brief Specify whether the relative error should be printed at each convergence check after 'num_check_iters' steps */ print_relative_error(bool b)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, 115 viennacl::matrix_base<ScalarType> & W, viennacl::matrix_base<ScalarType> & H, 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: 125 mutable vcl_size_t iters_; 126 }; 127 128 namespace host_based 129 { 130 /** @brief Missing OpenMP kernel for nonnegative matrix factorization of a dense matrices. */ 131 template<typename NumericT> el_wise_mul_div(NumericT * matrix1,NumericT const * matrix2,NumericT const * matrix3,vcl_size_t size)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 148 /** @brief The nonnegative matrix factorization (approximation) algorithm as suggested by Lee and Seung. Factorizes a matrix V with nonnegative entries into matrices W and H such that ||V - W*H|| is minimized. 149 * 150 * @param V Input matrix 151 * @param W First factor 152 * @param H Second factor 153 * @param conf A configuration object holding tolerances and the like 154 */ 155 template<typename NumericT> nmf(viennacl::matrix_base<NumericT> const & V,viennacl::matrix_base<NumericT> & W,viennacl::matrix_base<NumericT> & H,viennacl::linalg::nmf_config const & conf)156 void nmf(viennacl::matrix_base<NumericT> const & V, 157 viennacl::matrix_base<NumericT> & W, 158 viennacl::matrix_base<NumericT> & H, 159 viennacl::linalg::nmf_config const & conf) 160 { 161 vcl_size_t k = W.size2(); 162 conf.iters_ = 0; 163 164 if (viennacl::linalg::norm_frobenius(W) <= 0) 165 W = viennacl::scalar_matrix<NumericT>(W.size1(), W.size2(), NumericT(1.0)); 166 167 if (viennacl::linalg::norm_frobenius(H) <= 0) 168 H = viennacl::scalar_matrix<NumericT>(H.size1(), H.size2(), NumericT(1.0)); 169 170 viennacl::matrix_base<NumericT> wn(V.size1(), k, W.row_major()); 171 viennacl::matrix_base<NumericT> wd(V.size1(), k, W.row_major()); 172 viennacl::matrix_base<NumericT> wtmp(V.size1(), V.size2(), W.row_major()); 173 174 viennacl::matrix_base<NumericT> hn(k, V.size2(), H.row_major()); 175 viennacl::matrix_base<NumericT> hd(k, V.size2(), H.row_major()); 176 viennacl::matrix_base<NumericT> htmp(k, k, H.row_major()); 177 178 viennacl::matrix_base<NumericT> appr(V.size1(), V.size2(), V.row_major()); 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 196 viennacl::linalg::host_based::el_wise_mul_div(data_H, data_hn, data_hd, H.internal_size1() * H.internal_size2()); 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 206 viennacl::linalg::host_based::el_wise_mul_div(data_W, data_wn, data_wd, W.internal_size1() * W.internal_size2()); 207 208 if (i % conf.check_after_steps() == 0) //check for convergence 209 { 210 appr = viennacl::linalg::prod(W, H); 211 212 appr -= V; 213 NumericT diff_val = viennacl::linalg::norm_frobenius(appr); 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_ */ 248