All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
nmf_operations.hpp
Go to the documentation of this file.
#ifndef VIENNACL_LINALG_CUDA_NMF_OPERATIONS_HPP_
#define VIENNACL_LINALG_CUDA_NMF_OPERATIONS_HPP_
/* =========================================================================
Copyright (c) 2010-2016, Institute for Microelectronics,
Institute for Analysis and Scientific Computing,
TU Wien.
Portions of this software are copyright by UChicago Argonne, LLC.
-----------------
ViennaCL - The Vienna Computing Library
-----------------
Project Head: Karl Rupp rupp@iue.tuwien.ac.at
(A list of authors and contributors can be found in the manual)
License: MIT (X11), see file LICENSE in the base directory
============================================================================= */
#include "viennacl/linalg/host_based/nmf_operations.hpp"
#include "viennacl/linalg/cuda/common.hpp"
namespace viennacl
{
namespace linalg
{
namespace cuda
{
template<typename NumericT>
__global__ void el_wise_mul_div(NumericT * matrix1,
NumericT const * matrix2,
NumericT const * matrix3,
unsigned int size)
{
for (unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; i < size; i +=gridDim.x * blockDim.x)
{
NumericT val = matrix1[i] * matrix2[i];
NumericT divisor = matrix3[i];
matrix1[i] = (divisor > (NumericT) 0.00001) ? (val / divisor) : NumericT(0);
}
}
template<typename NumericT>
void nmf(viennacl::matrix_base<NumericT> const & V,
viennacl::matrix_base<NumericT> & W,
viennacl::matrix_base<NumericT> & H,
viennacl::linalg::nmf_config const & conf)
{
vcl_size_t k = W.size2();
conf.iters_ = 0;
if (!viennacl::linalg::norm_frobenius(W))
W = viennacl::scalar_matrix<NumericT>(W.size1(), W.size2(), NumericT(1.0));
if (!viennacl::linalg::norm_frobenius(H))
H = viennacl::scalar_matrix<NumericT>(H.size1(), H.size2(), NumericT(1.0));
viennacl::matrix_base<NumericT> wn(V.size1(), k, W.row_major());
viennacl::matrix_base<NumericT> wd(V.size1(), k, W.row_major());
viennacl::matrix_base<NumericT> wtmp(V.size1(), V.size2(), W.row_major());
viennacl::matrix_base<NumericT> hn(k, V.size2(), H.row_major());
viennacl::matrix_base<NumericT> hd(k, V.size2(), H.row_major());
viennacl::matrix_base<NumericT> htmp(k, k, H.row_major());
viennacl::matrix_base<NumericT> appr(V.size1(), V.size2(), V.row_major());
viennacl::vector<NumericT> diff(V.size1() * V.size2());
NumericT last_diff = 0;
NumericT diff_init = 0;
bool stagnation_flag = false;
for (vcl_size_t i = 0; i < conf.max_iterations(); i++)
{
conf.iters_ = i + 1;
hn = viennacl::linalg::prod(trans(W), V);
htmp = viennacl::linalg::prod(trans(W), W);
hd = viennacl::linalg::prod(htmp, H);
el_wise_mul_div<<<128, 128>>>(viennacl::cuda_arg<NumericT>(H),
viennacl::cuda_arg<NumericT>(hn),
viennacl::cuda_arg<NumericT>(hd),
static_cast<unsigned int>(H.internal_size1() * H.internal_size2()));
VIENNACL_CUDA_LAST_ERROR_CHECK("el_wise_mul_div");
wn = viennacl::linalg::prod(V, trans(H));
wtmp = viennacl::linalg::prod(W, H);
wd = viennacl::linalg::prod(wtmp, trans(H));
el_wise_mul_div<<<128, 128>>>(viennacl::cuda_arg<NumericT>(W),
viennacl::cuda_arg<NumericT>(wn),
viennacl::cuda_arg<NumericT>(wd),
static_cast<unsigned int>( W.internal_size1() * W.internal_size2()));
VIENNACL_CUDA_LAST_ERROR_CHECK("el_wise_mul_div");
if (i % conf.check_after_steps() == 0) //check for convergence
{
appr = viennacl::linalg::prod(W, H);
appr -= V;
NumericT diff_val = viennacl::linalg::norm_frobenius(appr);
if (i == 0)
diff_init = diff_val;
if (conf.print_relative_error())
std::cout << diff_val / diff_init << std::endl;
// Approximation check
if (diff_val / diff_init < conf.tolerance())
break;
// Stagnation check
if (std::fabs(diff_val - last_diff) / (diff_val * conf.check_after_steps()) < conf.stagnation_tolerance()) //avoid situations where convergence stagnates
{
if (stagnation_flag) // iteration stagnates (two iterates with no notable progress)
break;
else
// record stagnation in this iteration
stagnation_flag = true;
} else
// good progress in this iteration, so unset stagnation flag
stagnation_flag = false;
// prepare for next iterate:
last_diff = diff_val;
}
}
}
} //namespace cuda
} //namespace linalg
} //namespace viennacl
#endif /* VIENNACL_LINALG_CUDA_NMF_OPERATIONS_HPP_ */
viennacl::linalg::cuda::nmf
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...
Definition: nmf_operations.hpp:59
viennacl::linalg::cuda::trans
void trans(matrix_expression< const matrix_base< NumericT, SizeT, DistanceT >, const matrix_base< NumericT, SizeT, DistanceT >, op_trans > const &proxy, matrix_base< NumericT > &temp_trans)
Definition: matrix_operations.hpp:94
nmf_operations.hpp
viennacl::linalg::nmf_config::check_after_steps
vcl_size_t check_after_steps() const
Number of steps after which the convergence of NMF should be checked (again)
Definition: nmf_operations.hpp:91
viennacl::matrix_base< NumericT >
viennacl::linalg::nmf_config
Configuration class for the nonnegative-matrix-factorization algorithm. Specify tolerances, maximum iteration counts, etc., here.
Definition: nmf_operations.hpp:39
NumericT
float NumericT
Definition: bisect.cpp:40
viennacl::linalg::nmf_config::iters_
vcl_size_t iters_
Definition: nmf_operations.hpp:125
viennacl::linalg::prod
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:102
viennacl::traits::size
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:239
viennacl::linalg::cuda::el_wise_mul_div
__global__ void el_wise_mul_div(NumericT *matrix1, NumericT const *matrix2, NumericT const *matrix3, unsigned int size)
Main CUDA kernel for nonnegative matrix factorization of a dense matrices.
Definition: nmf_operations.hpp:38
viennacl::linalg::nmf_config::max_iterations
vcl_size_t max_iterations() const
Returns the maximum number of iterations for the NMF algorithm.
Definition: nmf_operations.hpp:74
viennacl::scalar_matrix
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
viennacl::linalg::nmf_config::tolerance
double tolerance() const
Returns the relative tolerance for convergence.
Definition: nmf_operations.hpp:50
viennacl::vcl_size_t
std::size_t vcl_size_t
Definition: forwards.h:75
viennacl::matrix_base::size2
size_type size2() const
Returns the number of columns.
Definition: matrix_def.hpp:226
viennacl::vector< NumericT >
viennacl::matrix_base::size1
size_type size1() const
Returns the number of rows.
Definition: matrix_def.hpp:224
viennacl::linalg::nmf_config::stagnation_tolerance
double stagnation_tolerance() const
Relative tolerance for the stagnation check.
Definition: nmf_operations.hpp:62
common.hpp
Common routines for CUDA execution.
viennacl::matrix_base::row_major
bool row_major() const
Definition: matrix_def.hpp:248
viennacl::linalg::norm_frobenius
scalar_expression< const matrix_base< NumericT >, const matrix_base< NumericT >, op_norm_frobenius > norm_frobenius(const matrix_base< NumericT > &A)
Definition: norm_frobenius.hpp:61
VIENNACL_CUDA_LAST_ERROR_CHECK
#define VIENNACL_CUDA_LAST_ERROR_CHECK(message)
Definition: common.hpp:30
viennacl::matrix_base::internal_size2
size_type internal_size2() const
Returns the internal number of columns. Usually required for launching OpenCL kernels only...
Definition: matrix_def.hpp:240
viennacl::matrix_base::internal_size1
size_type internal_size1() const
Returns the internal number of rows. Usually required for launching OpenCL kernels only...
Definition: matrix_def.hpp:238
diff
ScalarType diff(ScalarType &s1, viennacl::scalar< ScalarType > &s2)
Definition: blas3_solve.cpp:50
viennacl::linalg::nmf_config::print_relative_error
bool print_relative_error() const
Returns the flag specifying whether the relative tolerance should be printed in each iteration...
Definition: nmf_operations.hpp:103