bnmf-algs
bld_add.hpp
Go to the documentation of this file.
1 #pragma once
2 
4 #include "defs.hpp"
5 #include "util_details.hpp"
6 #include "util/util.hpp"
7 #include <gsl/gsl_sf_psi.h>
8 
9 namespace bnmf_algs {
10 namespace bld {
68 template <typename T, typename Scalar>
69 tensor_t<T, 3> bld_add(const matrix_t<T>& X, size_t z,
70  const alloc_model::Params<Scalar>& model_params,
71  size_t max_iter = 1000, double eps = 1e-50) {
72  details::check_bld_params(X, z, model_params);
73 
74  long x = X.rows(), y = X.cols();
75  tensor_t<T, 3> S(x, y, z);
76  // initialize tensor S
77  {
78  util::gsl_rng_wrapper rand_gen(gsl_rng_alloc(gsl_rng_taus),
79  gsl_rng_free);
80  std::vector<double> dirichlet_params(z, 1);
81  std::vector<double> dirichlet_variates(z);
82  for (int i = 0; i < x; ++i) {
83  for (int j = 0; j < y; ++j) {
84  gsl_ran_dirichlet(rand_gen.get(), z, dirichlet_params.data(),
85  dirichlet_variates.data());
86  for (size_t k = 0; k < z; ++k) {
87  S(i, j, k) = X(i, j) * dirichlet_variates[k];
88  }
89  }
90  }
91  }
92 
93  tensor_t<T, 2> S_pjk(y, z); // S_{+jk}
94  tensor_t<T, 2> S_ipk(x, z); // S_{i+k}
95  matrix_t<T> alpha_eph(x, z);
96  matrix_t<T> beta_eph(y, z);
97  tensor_t<T, 3> grad_S(x, y, z);
98  tensor_t<T, 3> eps_tensor(x, y, z);
99  tensor_t<T, 3> Z(x, y, z);
100  tensor_t<T, 2> S_mult(x, y);
101  eps_tensor.setConstant(eps);
102 
103  const auto eta = [](const size_t step) -> double {
104  return 0.1 / std::pow(step + 1, 0.55);
105  };
106 
107 #ifdef USE_OPENMP
108  Eigen::ThreadPool tp(std::thread::hardware_concurrency());
109  Eigen::ThreadPoolDevice thread_dev(&tp,
111 #endif
112 
113  for (size_t eph = 0; eph < max_iter; ++eph) {
114  // update S_pjk, S_ipk, S_ijp
115 #ifdef USE_OPENMP
116  S_pjk.device(thread_dev) = S.sum(shape<1>({0}));
117  S_ipk.device(thread_dev) = S.sum(shape<1>({1}));
118 #else
119  S_pjk = S.sum(shape<1>({0}));
120  S_ipk = S.sum(shape<1>({1}));
121 #endif
122  // update alpha_eph, beta_eph
123  #pragma omp parallel for schedule(static)
124  for (size_t k = 0; k < z; ++k) {
125  for (int i = 0; i < x; ++i) {
126  alpha_eph(i, k) = model_params.alpha[i] + S_ipk(i, k);
127  }
128  for (int j = 0; j < y; ++j) {
129  beta_eph(j, k) = model_params.beta[k] + S_pjk(j, k);
130  }
131  }
132 
133  // update grad_S
134  vector_t<T> alpha_eph_sum = alpha_eph.colwise().sum();
135  #pragma omp parallel for schedule(static)
136  for (int i = 0; i < x; ++i) {
137  for (int j = 0; j < y; ++j) {
138  for (size_t k = 0; k < z; ++k) {
139  grad_S(i, j, k) = util::psi_appr(beta_eph(j, k)) -
140  util::psi_appr(S(i, j, k) + 1) -
141  util::psi_appr(alpha_eph_sum(k)) +
142  util::psi_appr(alpha_eph(i, k));
143  }
144  }
145  }
146  // construct Z
147 #ifdef USE_OPENMP
148  Z.device(thread_dev) = S.log() + eps_tensor;
149  S_mult.device(thread_dev) = (S * grad_S).sum(shape<1>({2}));
150 #else
151  Z = S.log() + eps_tensor;
152  S_mult = (S * grad_S).sum(shape<1>({2}));
153 #endif
154 
155  #pragma omp parallel for schedule(static)
156  for (int i = 0; i < x; ++i) {
157  for (int j = 0; j < y; ++j) {
158  for (size_t k = 0; k < z; ++k) {
159  Z(i, j, k) += eta(eph) * S(i, j, k) *
160  (X(i, j) * grad_S(i, j, k) - S_mult(i, j));
161  }
162  }
163  }
164 #ifdef USE_OPENMP
165  Z.device(thread_dev) = Z.exp();
166 #else
167  Z = Z.exp();
168 #endif
170  // update S
171 
172  #pragma omp parallel for schedule(static)
173  for (int i = 0; i < x; ++i) {
174  for (int j = 0; j < y; ++j) {
175  for (size_t k = 0; k < z; ++k) {
176  S(i, j, k) = X(i, j) * Z(i, j, k);
177  }
178  }
179  }
180  }
181  return S;
182 }
183 } // namespace bld
184 } // namespace bnmf_algs
Structure to hold the parameters for the Allocation Model .
Definition: alloc_model_params.hpp:25
Eigen::Tensor< Scalar, N, Eigen::RowMajor > tensor_t
Tensor type used in the computations.
Definition: defs.hpp:52
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > matrix_t
Matrix type used in the computations.
Definition: defs.hpp:41
T hardware_concurrency(T...args)
Eigen::array< size_t, N > shape
Shape of vectors, matrices, tensors, etc.
Definition: defs.hpp:66
void normalize(tensor_t< Scalar, N > &input, size_t axis, NormType type=NormType::L1, double eps=1e-50)
Normalize a tensor in-place along the given axis.
Definition: util.hpp:321
T data(T...args)
std::vector< Scalar > alpha
Parameter vector of Dirichlet prior for matrix of size .
Definition: alloc_model_params.hpp:37
Real psi_appr(Real x)
Calculate digamma function approximately using the first 8 terms of the asymptotic expansion of digam...
Definition: util.hpp:440
tensor_t< T, 3 > bld_add(const matrix_t< T > &X, size_t z, const alloc_model::Params< Scalar > &model_params, size_t max_iter=1000, double eps=1e-50)
Compute tensor , the solution of BLD problem , from matrix using additive gradient ascent updates...
Definition: bld_add.hpp:69
T get(T...args)
void check_bld_params(const matrix_t< T > &X, size_t z, const alloc_model::Params< Scalar > &model_params)
Do parameter checks on BLD computing function parameters and throw an assertion error if the paramete...
Definition: util_details.hpp:20
T pow(T...args)
std::vector< Scalar > beta
Parameter vector of Dirichlet prior for matrix of size .
Definition: alloc_model_params.hpp:42
Eigen::Matrix< Scalar, 1, Eigen::Dynamic, Eigen::RowMajor > vector_t
Vector type used in the computations.
Definition: defs.hpp:27
Main namespace for bnmf-algs library.
Definition: alloc_model_funcs.hpp:12