Namespace that contains solver and auxiliary functions for Best Latent Decomposition (BLD) problem. More...

Classes
struct	EMResult
	Structure holding the results of EM procedures. More...

Functions
template<typename T , typename Scalar >
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 \(S\), the solution of BLD problem [5], from matrix \(X\) using additive gradient ascent updates. More...

template<typename T , typename Scalar >
std::tuple< tensor_t< T, 3 >, matrix_t< T >, matrix_t< T > >	bld_appr (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 \(S\), the solution of BLD problem [5], and optimization matrices \(\mu\) and \(\nu\) from matrix \(X\) using the approximate multiplicative algorithm given in [5]. More...

template<typename T , typename Scalar >
std::tuple< matrix_t< T >, matrix_t< T >, vector_t< T > >	bld_fact (const tensor_t< T, 3 > &S, const alloc_model::Params< Scalar > &model_params, double eps=1e-50)
	Compute matrices \(W, H\) and vector \(L\) from tensor \(S\) according to the allocation model. More...

template<typename T , typename Scalar >
tensor_t< T, 3 >	bld_mult (const matrix_t< T > &X, const size_t z, const alloc_model::Params< Scalar > &model_params, size_t max_iter=1000, bool use_psi_appr=false, double eps=1e-50)
	Compute tensor \(S\), the solution of BLD problem [5], from matrix \(X\) using multiplicative update rules. More...

template<typename T , typename Scalar >
tensor_t< T, 3 >	bld_mult_cuda (const matrix_t< T > &X, const size_t z, const alloc_model::Params< Scalar > &model_params, size_t max_iter=1000, bool use_psi_appr=false, double eps=1e-50)
	Compute tensor \(S\), the solution of BLD problem [5], from matrix \(X\) using multiplicative update rules. More...

template<typename T , typename Scalar >
util::Generator< tensor_t< T, 3 >, details::CollapsedGibbsComputer< T, Scalar > >	collapsed_gibbs (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 a sequence of \(S\) tensors using Collapsed Gibbs Sampler method. More...

template<typename T , typename Scalar >
tensor_t< T, 3 >	collapsed_icm (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 \(S\), the solution of BLD problem [5], from matrix \(X\) using collapsed iterated conditional modes. More...

template<typename T , typename Scalar >
EMResult< T >	online_EM (const matrix_t< T > &X, const std::vector< alloc_model::Params< Scalar >> &param_vec, const size_t max_iter=1000, const bool use_psi_appr=false)
	Complete a matrix containing unobserved values given as NaN using an EM procedure according to the allocation model proposed by [5] . More...

template<typename T , typename Scalar >
tensor_t< T, 3 >	seq_greedy_bld (const matrix_t< T > &X, size_t z, const alloc_model::Params< Scalar > &model_params)
	Compute tensor \(S\), the solution of BLD problem [5], from matrix \(X\) using sequential greedy method. More...

Detailed Description

Namespace that contains solver and auxiliary functions for Best Latent Decomposition (BLD) problem.

Function Documentation

template<typename T , typename Scalar >

tensor_t<T, 3> bnmf_algs::bld::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 \(S\), the solution of BLD problem [5], from matrix \(X\) using additive gradient ascent updates.

According to Allocation Model [5],

\[ L_j \sim \mathcal{G}(a, b) \qquad W_{:k} \sim \mathcal{D}(\alpha) \qquad H_{:j} \sim \mathcal{D}(\beta) \]

Each entry \(S_{ijk} \sim \mathcal{PO}(W_{ik}H_{kj}L_j)\) and overall \(X = S_{ij+}\).

In this context, Best Latent Decomposition (BLD) problem is [5],

\[ S^* = \underset{S_{::+}=X}{\arg \max}\text{ }p(S). \]

Additive BLD algorithm solves a constrained optimization problem over the set of tensors whose sum over their third index is equal to \(X\). Let's denote this set with \(\Omega_X = \{S | S_{ij+} = X_{ij}\}\). To make the optimization easier the constraints on the found tensors are relaxed by searching over \(\mathcal{R}^{x \times y \times z}\) and defining a projection function \(\phi : \mathcal{R}^{x \times y \times z} \rightarrow \Omega_X\) such that \(\phi(Z) = S\) where \(Z\) is the solution of the unconstrained problem and \(S\) is the solution of the original constrained problem.

In this algorithm, \(\phi\) is chosen such that

\[ S_{ijk} = \phi(Z_{ijk}) = X_{ij}\frac{\exp(Z_{ijk})}{\sum_{k'}\exp(Z_{ijk'})} \]

This particular choice of \(\phi\) leads to additive update rules implemented in this algorithm.

Template Parameters

T	Type of the matrix/tensor entries.
Scalar	Type of the scalars used as model parameters.

Parameters

X	Nonnegative matrix of size \(x \times y\) to decompose.
z	Number of matrices into which matrix \(X\) will be decomposed. This is the depth of the output tensor \(S\).
model_params	Allocation model parameters. See bnmf_algs::alloc_model::Params<double>.
max_iter	Maximum number of iterations.
eps	Floating point epsilon value to be used to prevent division by 0 errors.

Returns: Tensor \(S\) of size \(x \times y \times z\) where \(X = S_{ij+}\).

Exceptions

std::invalid_argument if X contains negative entries, if number of rows of X is not equal to number of alpha parameters, if z is not equal to number of beta parameters.

template<typename T , typename Scalar >

std::tuple<tensor_t<T, 3>, matrix_t<T>, matrix_t<T> > bnmf_algs::bld::bld_appr	(	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 \(S\), the solution of BLD problem [5], and optimization matrices \(\mu\) and \(\nu\) from matrix \(X\) using the approximate multiplicative algorithm given in [5].

According to Allocation Model [5],

\[ L_j \sim \mathcal{G}(a, b) \qquad W_{:k} \sim \mathcal{D}(\alpha) \qquad H_{:j} \sim \mathcal{D}(\beta) \]

Each entry \(S_{ijk} \sim \mathcal{PO}(W_{ik}H_{kj}L_j)\) and overall \(X = S_{ij+}\).

In this context, Best Latent Decomposition (BLD) problem is [5],

\[ S^* = \underset{S_{::+}=X}{\arg \max}\text{ }p(S). \]

This algorithm finds the optimal \(S\) using the multiplicative algorithm employing Lagrange multipliers as given in [kurutmazbayesian:]

\[ S_{ijk} = X_{ij}\frac{\nu_{ik}\mu_{kj}}{\sum_c \nu_{ic}\mu{cj}}. \]

Template Parameters

T	Type of the matrix/tensor entries.
Scalar	Type of the scalars used as model parameters.

Parameters

X	Nonnegative matrix of size \(x \times y\) to decompose.
z	Number of matrices into which matrix \(X\) will be decomposed. This is the depth of the output tensor \(S\).
model_params	Allocation model parameters. See bnmf_algs::alloc_model::Params<double>.
max_iter	Maximum number of iterations.
eps	Floating point epsilon value to be used to prevent division by 0 errors.

Returns: std::tuple of tensor \(S\) of size \(x \times y \times z\) where \(X = S_{ij+}\), and matrices \(\nu\) and \(\mu\).

Exceptions

std::invalid_argument if X contains negative entries, if number of rows of X is not equal to number of alpha parameters, if z is not equal to number of beta parameters.

template<typename T , typename Scalar >

std::tuple<matrix_t<T>, matrix_t<T>, vector_t<T> > bnmf_algs::bld::bld_fact	(	const tensor_t< T, 3 > &	S,
		const alloc_model::Params< Scalar > &	model_params,
		double	eps = `1e-50`
	)

Compute matrices \(W, H\) and vector \(L\) from tensor \(S\) according to the allocation model.

According to Allocation Model [5],

\[ L_j \sim \mathcal{G}(a, b) \qquad W_{:k} \sim \mathcal{D}(\alpha) \qquad H_{:j} \sim \mathcal{D}(\beta) \]

Each entry \(S_{ijk} \sim \mathcal{PO}(W_{ik}H_{kj}L_j)\).

Template Parameters

T	Type of the matrix/tensor entries.
Scalar	Type of the scalars used as model parameters.

Parameters

S	Tensor \(S\) of shape \(x \times y \times z\).
model_params	Allocation model parameters. See bnmf_algs::alloc_model::Params<double>.
eps	Floating point epsilon value to be used to prevent division by 0 errors.

Returns: std::tuple of matrices \(W_{x \times z}, H_{z \times y}\) and vector \(L_y\).

template<typename T , typename Scalar >

tensor_t<T, 3> bnmf_algs::bld::bld_mult	(	const matrix_t< T > &	X,
		const size_t	z,
		const alloc_model::Params< Scalar > &	model_params,
		size_t	max_iter = `1000`,
		bool	use_psi_appr = `false`,
		double	eps = `1e-50`
	)

Compute tensor \(S\), the solution of BLD problem [5], from matrix \(X\) using multiplicative update rules.

According to Allocation Model [5],

\[ L_j \sim \mathcal{G}(a, b) \qquad W_{:k} \sim \mathcal{D}(\alpha) \qquad H_{:j} \sim \mathcal{D}(\beta) \]

Each entry \(S_{ijk} \sim \mathcal{PO}(W_{ik}H_{kj}L_j)\) and overall \(X = S_{ij+}\).

In this context, Best Latent Decomposition (BLD) problem is [5],

\[ S^* = \underset{S_{::+}=X}{\arg \max}\text{ }p(S). \]

Multiplicative BLD algorithm solves a constrained optimization problem over the set of tensors whose sum over their third index is equal to \(X\). Let's denote this set with \(\Omega_X = \{S | S_{ij+} = X_{ij}\}\). To make the optimization easier the constraints on the found tensors are relaxed by searching over \(\mathcal{R}^{x \times y \times z}\) and defining a projection function \(\phi : \mathcal{R}^{x \times y \times z} \rightarrow \Omega_X\) such that \(\phi(Z) = S\) where \(Z\) is the solution of the unconstrained problem and \(S\) is the solution of the original constrained problem.

In this algorithm, \(\phi\) is chosen such that

\[ S_{ijk} = \phi(Z_{ijk}) = \frac{Z_{ijk}X_{ij}}{\sum_{k'}Z_{ijk'}} \]

This particular choice of \(\phi\) leads to multiplicative update rules implemented in this algorithm.

Template Parameters

T	Type of the matrix/tensor entries.
Scalar	Type of the scalars used as model parameters.

Parameters

X	Nonnegative matrix of size \(x \times y\) to decompose.
z	Number of matrices into which matrix \(X\) will be decomposed. This is the depth of the output tensor \(S\).
model_params	Allocation model parameters. See bnmf_algs::alloc_model::Params<double>.
max_iter	Maximum number of iterations.
use_psi_appr	If true, use util::psi_appr function to compute the digamma function approximately. If false, compute the digamma function exactly.
eps	Floating point epsilon value to be used to prevent division by 0 errors.

Returns: Tensor \(S\) of size \(x \times y \times z\) where \(X = S_{ij+}\).

Exceptions

assertion error if X contains negative entries, if number of rows of X is not equal to number of alpha parameters, if z is not equal to number of beta parameters.

template<typename T , typename Scalar >

tensor_t<T, 3> bnmf_algs::bld::bld_mult_cuda	(	const matrix_t< T > &	X,
		const size_t	z,
		const alloc_model::Params< Scalar > &	model_params,
		size_t	max_iter = `1000`,
		bool	use_psi_appr = `false`,
		double	eps = `1e-50`
	)

Compute tensor \(S\), the solution of BLD problem [5], from matrix \(X\) using multiplicative update rules.

According to Allocation Model [5],

\[ L_j \sim \mathcal{G}(a, b) \qquad W_{:k} \sim \mathcal{D}(\alpha) \qquad H_{:j} \sim \mathcal{D}(\beta) \]

Each entry \(S_{ijk} \sim \mathcal{PO}(W_{ik}H_{kj}L_j)\) and overall \(X = S_{ij+}\).

In this context, Best Latent Decomposition (BLD) problem is [5],

\[ S^* = \underset{S_{::+}=X}{\arg \max}\text{ }p(S). \]

Multiplicative BLD algorithm solves a constrained optimization problem over the set of tensors whose sum over their third index is equal to \(X\). Let's denote this set with \(\Omega_X = \{S | S_{ij+} = X_{ij}\}\). To make the optimization easier the constraints on the found tensors are relaxed by searching over \(\mathcal{R}^{x \times y \times z}\) and defining a projection function \(\phi : \mathcal{R}^{x \times y \times z} \rightarrow \Omega_X\) such that \(\phi(Z) = S\) where \(Z\) is the solution of the unconstrained problem and \(S\) is the solution of the original constrained problem.

In this algorithm, \(\phi\) is chosen such that

\[ S_{ijk} = \phi(Z_{ijk}) = \frac{Z_{ijk}X_{ij}}{\sum_{k'}Z_{ijk'}} \]

This particular choice of \(\phi\) leads to multiplicative update rules implemented in this algorithm.

This function is provided only when the library is built with CUDA support. Large matrix and tensor updates employed in bld_mult algorithm are performed using CUDA. In particular, updating the tensor S, its gradient tensors and certain highly parallelized \(O(xyz)\) operations are performed on the GPU to speed up the computation. Tensors and matrices are copied to/from the GPU only when needed. Additionally, objects that are only used during GPU updates are never allocated on the main memory; they reside only on the GPU for higher efficiency.

Template Parameters

T	Type of the matrix/tensor entries.
Scalar	Type of the scalars used as model parameters.

Parameters

X	Nonnegative matrix of size \(x \times y\) to decompose.
z	Number of matrices into which matrix \(X\) will be decomposed. This is the depth of the output tensor \(S\).
model_params	Allocation model parameters. See bnmf_algs::alloc_model::Params<double>.
max_iter	Maximum number of iterations.
use_psi_appr	If true, use util::psi_appr function to compute the digamma function approximately. If false, compute the digamma function exactly. Note that updates on the GPU use bnmf_algs::details::kernel::psi_appr function since they can't use any third party library functions.
eps	Floating point epsilon value to be used to prevent division by 0 errors.

Returns: Tensor \(S\) of size \(x \times y \times z\) where \(X = S_{ij+}\).

Exceptions

assertion error if X contains negative entries, if number of rows of X is not equal to number of alpha parameters, if z is not equal to number of beta parameters.

template<typename T , typename Scalar >

util::Generator<tensor_t<T, 3>, details::CollapsedGibbsComputer<T, Scalar> > bnmf_algs::bld::collapsed_gibbs	(	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 a sequence of \(S\) tensors using Collapsed Gibbs Sampler method.

According to Allocation Model [5],

\[ L_j \sim \mathcal{G}(a, b) \qquad W_{:k} \sim \mathcal{D}(\alpha) \qquad H_{:j} \sim \mathcal{D}(\beta) \]

Each entry \(S_{ijk} \sim \mathcal{PO}(W_{ik}H_{kj}L_j)\) and overall \(X = S_{ij+}\).

In this context, Best Latent Decomposition (BLD) problem is [5],

\[ S^* = \underset{S_{::+}=X}{\arg \max}\text{ }p(S). \]

Todo:: Explain collapsed gibbs sampler algorithm.

Template Parameters

T	Type of the matrix/tensor entries.
Scalar	Type of the scalars used as model parameters.

Parameters

X	Nonnegative matrix of size \(x \times y\) to decompose.
z	Number of matrices into which matrix \(X\) will be decomposed. This is the depth of the output tensor \(S\).
model_params	Allocation model parameters. See bnmf_algs::alloc_model::Params<double>.
max_iter	Maximum number of iterations.
eps	Floating point epsilon value to be used to prevent division by 0 errors.

Returns: util::Generator object that will generate a sequence of \(S\) tensors using details::CollapsedGibbsComputer as its Computer type.

Exceptions

std::invalid_argument if X contains negative entries, if number of rows of X is not equal to number of alpha parameters, if z is not equal to number of beta parameters.

template<typename T , typename Scalar >

tensor_t<T, 3> bnmf_algs::bld::collapsed_icm	(	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 \(S\), the solution of BLD problem [5], from matrix \(X\) using collapsed iterated conditional modes.

According to Allocation Model [5],

\[ L_j \sim \mathcal{G}(a, b) \qquad W_{:k} \sim \mathcal{D}(\alpha) \qquad H_{:j} \sim \mathcal{D}(\beta) \]

Each entry \(S_{ijk} \sim \mathcal{PO}(W_{ik}H_{kj}L_j)\) and overall \(X = S_{ij+}\).

In this context, Best Latent Decomposition (BLD) problem is [5],

\[ S^* = \underset{S_{::+}=X}{\arg \max}\text{ }p(S). \]

Todo:: Explain collapsed icm algorithm.

Template Parameters

T	Type of the matrix/tensor entries.
Scalar	Type of the scalars used as model parameters.

Parameters

X	Nonnegative matrix of size \(x \times y\) to decompose.
z	Number of matrices into which matrix \(X\) will be decomposed. This is the depth of the output tensor \(S\).
model_params	Allocation model parameters. See bnmf_algs::alloc_model::Params<double>.
max_iter	Maximum number of iterations.
eps	Floating point epsilon value to be used to prevent division by 0 errors.

Returns: Tensor \(S\) of size \(x \times y \times z\) where \(X = S_{ij+}\).

Exceptions

std::invalid_argument if X contains negative entries, if number of rows of X is not equal to number of alpha parameters, if z is not equal to number of beta parameters.

template<typename T , typename Scalar >

EMResult< T > bnmf_algs::bld::online_EM	(	const matrix_t< T > &	X,
		const std::vector< alloc_model::Params< Scalar >> &	param_vec,
		const size_t	max_iter = `1000`,
		const bool	use_psi_appr = `false`
	)

Complete a matrix containing unobserved values given as NaN using an EM procedure according to the allocation model proposed by [5] .

This procedure completes a matrix \(X\) containing unobserved values as NaNs. The completion is performed using an EM procedure derived according to the allocation model. The hidden tensor described in [5] is never stored explicitly; only its sums along its each dimension is stored. Therefore, space complexity of the method is \(O(x \times y)\) for an input matrix \(X_{x \times y}\).

Template Parameters

T	Type of the matrix entries.
Scalar	Type of the parameters.

Parameters

X	Incomplete matrix containing NaN values for every unobserved entry.
param_vec	Parameter vector of size \(z\), the rank of the decomposition. Each entry of the parameter vector must be an alloc_model::Params object containing \(x\) many \(\alpha\) and \(y\) many \(\beta\) values for the input matrix \(X_{x \times y}\). Note that for this procedure, the rank of the decomposition is not given explicitly; it is inferred as the size of the parameter vector.
max_iter	Maximum number of iterations to run the EM procedure.
use_psi_appr	If true, use util::psi_appr function to compute the digamma function approximately. If false, compute the digamma function exactly.

Returns: Results of the EM procedure as EMResult. See EMResult documentation for explanations of each result.

Remarks: Throws assertion error if matrix X contains negative values, if the given parameter vector is not compatible with the given matrix.

template<typename T , typename Scalar >

tensor_t<T, 3> bnmf_algs::bld::seq_greedy_bld	(	const matrix_t< T > &	X,
		size_t	z,
		const alloc_model::Params< Scalar > &	model_params
	)

Compute tensor \(S\), the solution of BLD problem [5], from matrix \(X\) using sequential greedy method.

According to Allocation Model [5],

\[ L_j \sim \mathcal{G}(a, b) \qquad W_{:k} \sim \mathcal{D}(\alpha) \qquad H_{:j} \sim \mathcal{D}(\beta) \]

Each entry \(S_{ijk} \sim \mathcal{PO}(W_{ik}H_{kj}L_j)\) and overall \(X = S_{ij+}\).

In this context, Best Latent Decomposition (BLD) problem is [5],

\[ S^* = \underset{S_{::+}=X}{\arg \max}\text{ }p(S). \]

Sequential greedy BLD algorithm works by greedily allocating each entry of input matrix, \(X_{ij}\), one by one to its corresponding bin \(S_{ijk}\). \(k\) is chosen by considering the increase in log marginal value of tensor \(S\) when \(S_{ijk}\) is increased by \(1\). At a given time, \(i, j\) is chosen as the indices of the largest entry of matrix \(X\).

Template Parameters

T	Type of the matrix/tensor entries.
Scalar	Type of the scalars used as model parameters.

Parameters

X	Nonnegative matrix of size \(x \times y\) to decompose.
z	Number of matrices into which matrix \(X\) will be decomposed. This is the depth of the output tensor \(S\).
model_params	Allocation model parameters. See bnmf_algs::alloc_model::Params<double>.

Returns: Tensor \(S\) of size \(x \times y \times z\) where \(X = S_{ij+}\).

Exceptions

std::invalid_argument if matrix \(X\) is not nonnegative, if size of alpha is not \(x\), if size of beta is not \(z\).

Remarks: Worst-case time complexity is \(O(xy + z(\log{xy})\sum_{ij}X_{ij})\) for a matrix \(X_{x \times y}\) and output tensor \(S_{x \times y \times z}\).

Classes

Functions

Detailed Description

Function Documentation