bnmf-algs
nmf.hpp
Go to the documentation of this file.
1 #pragma once
2 
3 #include <cmath>
4 #include <limits>
5 #include <utility>
6 
7 #include "defs.hpp"
8 
12 namespace bnmf_algs {
17 namespace nmf {
48 template <typename T, typename Real>
50  Real beta, size_t max_iter = 1000) {
51  const auto x = static_cast<size_t>(X.rows());
52  const auto y = static_cast<size_t>(X.cols());
53  const auto z = r;
54 
55  BNMF_ASSERT((X.array() >= 0).all(), "X matrix must be nonnegative");
56  BNMF_ASSERT(r > 0, "r must be positive");
57 
58  if (X.isZero(0)) {
59  return std::make_pair(matrixd::Zero(x, z), matrixd::Zero(z, y));
60  }
61 
64 
65  if (max_iter == 0) {
66  return std::make_pair(W, H);
67  }
68 
69  const Real p = 2 - beta;
70  const double eps = std::numeric_limits<double>::epsilon();
71  const auto p_int = static_cast<int>(p);
72  const bool p_is_int = std::abs(p_int - p) <= eps;
73 
74  matrix_t<T> X_hat, nom_W(x, z), denom_W(x, z), nom_H(z, y), denom_H(z, y);
75  while (max_iter-- > 0) {
76  // update approximation
77  X_hat = W * H;
78 
79  // update W
80  if (p_is_int && p == 0) {
81  // Euclidean NMF
82  nom_W = X * H.transpose();
83  denom_W = X_hat * H.transpose();
84  } else if (p_is_int && p == 1) {
85  // KL NMF
86  nom_W = (X.array() / X_hat.array()).matrix() * H.transpose();
87  denom_W = matrixd::Ones(x, y) * H.transpose();
88  } else if (p_is_int && p == 2) {
89  // Itakuro-Saito NMF
90  nom_W =
91  (X.array() / X_hat.array().square()).matrix() * H.transpose();
92  denom_W = (X_hat.array().inverse()).matrix() * H.transpose();
93  } else {
94  // general case
95  nom_W = (X.array() / X_hat.array().pow(p)).matrix() * H.transpose();
96  denom_W = (X_hat.array().pow(1 - p)).matrix() * H.transpose();
97  }
98  W = W.array() * nom_W.array() / denom_W.array();
99 
100  // update approximation
101  X_hat = W * H;
102 
103  // update H
104  if (p_is_int && p == 0) {
105  // Euclidean NMF
106  nom_H = W.transpose() * X;
107  denom_H = W.transpose() * X_hat;
108  } else if (p_is_int && p == 1) {
109  // KL NMF
110  nom_H = W.transpose() * (X.array() / X_hat.array()).matrix();
111  denom_H = W.transpose() * matrixd::Ones(x, y);
112  } else if (p_is_int && p == 2) {
113  // Itakuro-Saito NMF
114  nom_H =
115  W.transpose() * (X.array() / X_hat.array().square()).matrix();
116  denom_H = W.transpose() * (X_hat.array().inverse()).matrix();
117  } else {
118  // general case
119  nom_H = W.transpose() * (X.array() / X_hat.array().pow(p)).matrix();
120  denom_H = W.transpose() * (X_hat.array().pow(1 - p)).matrix();
121  }
122  H = H.array() * nom_H.array() / denom_H.array();
123  }
124  return std::make_pair(W, H);
125 }
126 
157 template <typename Real>
158 Real beta_divergence(Real x, Real y, Real beta, double eps = 1e-50) {
159  if (std::abs(beta) <= std::numeric_limits<double>::epsilon()) {
160  // Itakuro Saito
161  if (std::abs(x) <= std::numeric_limits<double>::epsilon()) {
162  return -1;
163  }
164  return x / (y + eps) - std::log(x / (y + eps)) - 1;
165  } else if (std::abs(beta - 1) <= std::numeric_limits<double>::epsilon()) {
166  // KL
167  Real logpart;
168  if (std::abs(x) <= std::numeric_limits<double>::epsilon()) {
169  logpart = 0;
170  } else {
171  logpart = x * std::log(x / (y + eps));
172  }
173  return logpart - x + y;
174  }
175 
176  // general case
177  Real nom = std::pow(x, beta) + (beta - 1) * std::pow(y, beta) -
178  beta * x * std::pow(y, beta - 1);
179  // beta != 0 && beta != 1
180  return nom / (beta * (beta - 1));
181 }
182 
200 template <typename Tensor>
201 typename Tensor::value_type beta_divergence(const Tensor& X, const Tensor& Y,
202  typename Tensor::value_type beta,
203  double eps = 1e-50) {
204  typename Tensor::value_type result = 0;
205  #pragma omp parallel for schedule(static) reduction(+:result)
206  for (long i = 0; i < X.rows(); ++i) {
207  for (long j = 0; j < X.cols(); ++j) {
208  result += beta_divergence(X(i, j), Y(i, j), beta, eps);
209  }
210  }
211 
212  return result;
213 }
214 
234 template <typename InputIterator1, typename InputIterator2>
235 auto beta_divergence_seq(InputIterator1 first_begin, InputIterator1 first_end,
236  InputIterator2 second_begin,
237  decltype(*first_begin) beta, double eps = 1e-50) {
238  // todo: parallelize
239  return std::inner_product(first_begin, first_end, second_begin, 0.0,
240  std::plus<>(), [beta, eps](double x, double y) {
241  return beta_divergence(x, y, beta, eps);
242  });
243 }
244 } // namespace nmf
245 } // namespace bnmf_algs
auto beta_divergence_seq(InputIterator1 first_begin, InputIterator1 first_end, InputIterator2 second_begin, decltype(*first_begin) beta, double eps=1e-50)
Compute the -divergence as defined in .
Definition: nmf.hpp:235
std::pair< matrix_t< T >, matrix_t< T > > nmf(const matrix_t< T > &X, size_t r, Real beta, size_t max_iter=1000)
Compute nonnegative matrix factorization of a matrix.
Definition: nmf.hpp:49
T epsilon(T...args)
T log(T...args)
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > matrix_t
Matrix type used in the computations.
Definition: defs.hpp:41
T make_pair(T...args)
T pow(T...args)
Real beta_divergence(Real x, Real y, Real beta, double eps=1e-50)
Compute the -divergence as defined in .
Definition: nmf.hpp:158
Main namespace for bnmf-algs library.
Definition: alloc_model_funcs.hpp:12
T inner_product(T...args)