1 #ifndef VIENNACL_LINALG_OPENCL_NMF_OPERATIONS_HPP_
2 #define VIENNACL_LINALG_OPENCL_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/opencl/vector_operations.hpp
22 @brief Implementations of NMF operations using OpenCL
23 */
24
25 #include "viennacl/linalg/opencl/kernels/nmf.hpp"
26 #include "viennacl/linalg/opencl/nmf_operations.hpp"
27
28 #include "viennacl/linalg/host_based/nmf_operations.hpp"
29
30 namespace viennacl
31 {
32 namespace linalg
33 {
34 namespace opencl
35 {
36
37 /** @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.
38 *
39 * @param V Input matrix
40 * @param W First factor
41 * @param H Second factor
42 * @param conf A configuration object holding tolerances and the like
43 */
44 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)45 void nmf(viennacl::matrix_base<NumericT> const & V,
46 viennacl::matrix_base<NumericT> & W,
47 viennacl::matrix_base<NumericT> & H,
48 viennacl::linalg::nmf_config const & conf)
49 {
50 viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(V).context());
51
52 const std::string NMF_MUL_DIV_KERNEL = "el_wise_mul_div";
53
54 viennacl::linalg::opencl::kernels::nmf<NumericT>::init(ctx);
55
56 vcl_size_t k = W.size2();
57 conf.iters_ = 0;
58
59 if (viennacl::linalg::norm_frobenius(W) <= 0)
60 W = viennacl::scalar_matrix<NumericT>(W.size1(), W.size2(), NumericT(1), ctx);
61
62 if (viennacl::linalg::norm_frobenius(H) <= 0)
63 H = viennacl::scalar_matrix<NumericT>(H.size1(), H.size2(), NumericT(1), ctx);
64
65 viennacl::matrix_base<NumericT> wn(V.size1(), k, W.row_major(), ctx);
66 viennacl::matrix_base<NumericT> wd(V.size1(), k, W.row_major(), ctx);
67 viennacl::matrix_base<NumericT> wtmp(V.size1(), V.size2(), W.row_major(), ctx);
68
69 viennacl::matrix_base<NumericT> hn(k, V.size2(), H.row_major(), ctx);
70 viennacl::matrix_base<NumericT> hd(k, V.size2(), H.row_major(), ctx);
71 viennacl::matrix_base<NumericT> htmp(k, k, H.row_major(), ctx);
72
73 viennacl::matrix_base<NumericT> appr(V.size1(), V.size2(), V.row_major(), ctx);
74
75 NumericT last_diff = 0;
76 NumericT diff_init = 0;
77 bool stagnation_flag = false;
78
79 for (vcl_size_t i = 0; i < conf.max_iterations(); i++)
80 {
81 conf.iters_ = i + 1;
82 {
83 hn = viennacl::linalg::prod(trans(W), V);
84 htmp = viennacl::linalg::prod(trans(W), W);
85 hd = viennacl::linalg::prod(htmp, H);
86
87 viennacl::ocl::kernel & mul_div_kernel = ctx.get_kernel(viennacl::linalg::opencl::kernels::nmf<NumericT>::program_name(), NMF_MUL_DIV_KERNEL);
88 viennacl::ocl::enqueue(mul_div_kernel(H, hn, hd, cl_uint(H.internal_size1() * H.internal_size2())));
89 }
90 {
91 wn = viennacl::linalg::prod(V, trans(H));
92 wtmp = viennacl::linalg::prod(W, H);
93 wd = viennacl::linalg::prod(wtmp, trans(H));
94
95 viennacl::ocl::kernel & mul_div_kernel = ctx.get_kernel(viennacl::linalg::opencl::kernels::nmf<NumericT>::program_name(), NMF_MUL_DIV_KERNEL);
96
97 viennacl::ocl::enqueue(mul_div_kernel(W, wn, wd, cl_uint(W.internal_size1() * W.internal_size2())));
98 }
99
100 if (i % conf.check_after_steps() == 0) //check for convergence
101 {
102 appr = viennacl::linalg::prod(W, H);
103
104 appr -= V;
105 NumericT diff_val = viennacl::linalg::norm_frobenius(appr);
106
107 if (i == 0)
108 diff_init = diff_val;
109
110 if (conf.print_relative_error())
111 std::cout << diff_val / diff_init << std::endl;
112
113 // Approximation check
114 if (diff_val / diff_init < conf.tolerance())
115 break;
116
117 // Stagnation check
118 if (std::fabs(diff_val - last_diff) / (diff_val * NumericT(conf.check_after_steps())) < conf.stagnation_tolerance()) //avoid situations where convergence stagnates
119 {
120 if (stagnation_flag) // iteration stagnates (two iterates with no notable progress)
121 break;
122 else
123 // record stagnation in this iteration
124 stagnation_flag = true;
125 } else
126 // good progress in this iteration, so unset stagnation flag
127 stagnation_flag = false;
128
129 // prepare for next iterate:
130 last_diff = diff_val;
131 }
132 }
133 }
134
135 } //namespace opencl
136 } //namespace linalg
137 } //namespace viennacl
138
139 #endif /* VIENNACL_LINALG_OPENCL_NMF_OPERATIONS_HPP_ */
140