bayesmix/utils

Utils

Collection of miscellaneous auxiliary tools for the library.

Clustering utilities

The cluster_utils.h file includes some utilities for cluster estimation. These functions only use Eigen objects.

namespace bayesmix

Functions

Eigen::MatrixXd posterior_similarity(const Eigen::MatrixXd &alloc_chain)

Computes the posterior similarity matrix the data.

Eigen::VectorXi cluster_estimate(const Eigen::MatrixXi &alloc_chain)

Estimates the clustering structure of the data via LS minimization.

Distribution-related utilities

The distributions.h file includes several useful functions related to probability distributions, including categorical variables, popular multivariate distributions, and distribution distances. Some of these functions make use of OpenMP parallelism to achieve better efficiency.

namespace bayesmix

Functions

int categorical_rng(const Eigen::VectorXd &probas, std::mt19937 &rng, const int start = 0)

Returns a pseudorandom categorical random variable on the set {start, …, start + k} where k is the size of the given probability vector

Parameters:

• probas – Probabilities for each category

• rng – Random number generator

• start – (default = 0)

Returns:

categorical r.v. with values on {start, …, start + k}

double multi_normal_prec_lpdf(const Eigen::VectorXd &datum, const Eigen::VectorXd &mean, const Eigen::MatrixXd &prec_chol, const double prec_logdet)

Evaluates the log probability density function of a multivariate Gaussian distribution parametrized by mean and precision matrix on a single point

Parameters:

• datum – Point in which to evaluate the the lpdf

• mean – The mean of the Gaussian distribution

• prec_chol – The (lower) Cholesky factor of the precision matrix

• prec_logdet – The logarithm of the determinant of the precision matrix

Returns:

The evaluation of the lpdf

Eigen::VectorXd multi_normal_prec_lpdf_grid(const Eigen::MatrixXd &data, const Eigen::VectorXd &mean, const Eigen::MatrixXd &prec_chol, const double prec_logdet)

Evaluates the log probability density function of a multivariate Gaussian distribution parametrized by mean and precision matrix on multiple points

Parameters:

• data – Grid of points (by row) on which to evaluate the lpdf

• mean – The mean of the Gaussian distribution

• prec_chol – The (lower) Cholesky factor of the precision matrix

• prec_logdet – The logarithm of the determinant of the precision matrix

Returns:

The evaluation of the lpdf

Eigen::VectorXd multi_normal_diag_rng(const Eigen::VectorXd &mean, const Eigen::DiagonalMatrix<double, Eigen::Dynamic> &cov_diag, std::mt19937 &rng)

Returns a pseudorandom multivariate normal random variable with diagonal covariance matrix

Parameters:

• mean – The mean of the Gaussian r.v.

• cov_diag – The diagonal covariance matrix

• rng – Random number generator

Returns:

Multivariate normal r.v.

Eigen::VectorXd multi_normal_prec_chol_rng(const Eigen::VectorXd &mean, const Eigen::LLT<Eigen::MatrixXd> &prec_chol, std::mt19937 &rng)

Returns a pseudorandom multivariate normal random variable parametrized through mean and Cholesky decomposition of precision matrix

Parameters:

• mean – The mean of the Gaussian r.v.

• prec_chol – The (lower) Cholesky factor of the precision matrix

• rng – Random number generator

Returns:

Multivariate normal r.v.

double multi_normal_lpdf_woodbury_chol(const Eigen::RowVectorXd &datum, const Eigen::VectorXd &mean, const Eigen::DiagonalMatrix<double, Eigen::Dynamic> &sigma_diag_inverse, const Eigen::MatrixXd &wood_factor, const double &cov_logdet)

Evaluates the log probability density function of a multivariate Gaussian distribution with the following covariance structure: Sigma + Lambda * Lambda^T where Sigma is a diagonal matrix and Lambda a (p x d) one. Usually, p >> d. y^T*(Sigma + Lambda * Lambda^T)^{-1}*y = y^T*Sigma^{-1}*y - ||wood_factor*y||^2

Parameters:

• datum – Point on which to evaluate the lpdf

• mean – The mean of the Gaussian distribution

• sigma_diag_inverse – The inverse of the diagonal of Sigma matrix

• wood_factor – Computed as L^{-1} * Lambda^T * Sigma^{-1}, where L is the (lower) Cholesky factor of I + Lambda^T * Sigma^{-1} * Lambda

• cov_logdet – The logarithm of the determinant of the covariance matrix

Returns:

The evaluation of the lpdf

double multi_normal_lpdf_woodbury(const Eigen::VectorXd &datum, const Eigen::VectorXd &mean, const Eigen::VectorXd &sigma_diag, const Eigen::MatrixXd &lambda)

Evaluates the log probability density function of a multivariate Gaussian distribution with the following covariance structure: Sigma + Lambda * Lambda^T where Sigma is a diagonal matrix and Lambda a (p x d) one. Usually, p >> d. The Woodbury matrix identity (https://en.wikipedia.org/wiki/Woodbury_matrix_identity) is used to turn computation from being O(p^3) to being O(d^3 p) which gives a substantial speedup when p >> d

Parameters:

• datum – Point on which to evaluate the lpdf

• mean – The mean of the Gaussian distribution

• sigma_diag – The diagonal of Sigma matrix

• lambda – Rectangular matrix in Woodbury Identity

Returns:

The evaluation of the lpdf

std::pair<Eigen::MatrixXd, double> compute_wood_chol_and_logdet(const Eigen::DiagonalMatrix<double, Eigen::Dynamic> &sigma_diag_inverse, const Eigen::MatrixXd &lambda)

Returns the log-determinant of the matrix and the ‘wood_factor’, i.e. L^{-1} * Lambda^T * Sigma^{-1}, where L is the (lower) Cholesky factor of I + Lambda^T * Sigma^{-1} * Lambda

Parameters:

• sigma_diag_inverse – The inverse of the diagonal matrix Sigma

• lambda – The matrix Lambda

double multi_student_t_invscale_lpdf(const Eigen::VectorXd &datum, const double df, const Eigen::VectorXd &mean, const Eigen::MatrixXd &invscale_chol, const double scale_logdet)

Evaluates the log probability density function of a multivariate Student’s t distribution on a single point

Parameters:

• datum – Point in which to evaluate the the lpdf

• df – The degrees of freedom of the Student’s t distribution

• mean – The mean of the Student’s t distribution

• invscale_chol – The (lower) Cholesky factor of the inverse scale matrix

• scale_logdet – The logarithm of the determinant of the inverse scale matrix

Returns:

The evaluation of the lpdf

Eigen::VectorXd multi_student_t_invscale_lpdf_grid(const Eigen::MatrixXd &data, const double df, const Eigen::VectorXd &mean, const Eigen::MatrixXd &invscale_chol, const double scale_logdet)

Evaluates the log probability density function of a multivariate Student’s t distribution on multiple points

Parameters:

• data – Grid of points (by row) on which to evaluate the lpdf

• df – The degrees of freedom of the Student’s t distribution

• mean – The mean of the Student’s t distribution

• invscale_chol – The (lower) Cholesky factor of the inverse scale matrix

• scale_logdet – The logarithm of the determinant of the inverse scale matrix

Returns:

The evaluation of the lpdf

double gaussian_mixture_dist(const Eigen::VectorXd &means1, const Eigen::VectorXd &vars1, const Eigen::VectorXd &weights1, const Eigen::VectorXd &means2, const Eigen::VectorXd &vars2, const Eigen::VectorXd &weights2)

Computes the L^2 distance between the univariate mixture of Gaussian densities p1(x) = sum_{h=1}^m1 w1[h] N(x | mean1[h], var1[h]) and p2(x) = sum_{h=1}^m2 w2[h] N(x | mean2[h], var2[h])

The L^2 distance amounts to d(p, q) = (int (p(x) - q(x)^2 dx))^{1/2}

double gaussian_mixture_dist(const std::vector<Eigen::VectorXd> &means1, const std::vector<Eigen::MatrixXd> &precs1, const Eigen::VectorXd &weights1, const std::vector<Eigen::VectorXd> &means2, const std::vector<Eigen::MatrixXd> &precs2, const Eigen::VectorXd &weights2)

Computes the L^2 distance between the multivariate mixture of Gaussian densities p1(x) = sum_{h=1}^m1 w1[h] N(x | mean1[h], Prec[1]^{-1}) and p2(x) = sum_{h=1}^m2 w2[h] N(x | mean2[h], Prec2[h]^{-1})

The L^2 distance amounts to d(p, q) = (int (p(x) - q(x)^2 dx))^{1/2}

double gaussian_mixture_dist(const std::vector<bayesmix::AlgorithmState::ClusterState> &clus1, const Eigen::VectorXd &weights1, const std::vector<bayesmix::AlgorithmState::ClusterState> &clus2, const Eigen::VectorXd &weights2)

Computes the L^2 distance between the mixture of Gaussian densities p(x) and q(x). These could be either univariate or multivariate. The L2 distance amounts to d(p, q) = (int (p(x) - q(x)^2 dx))^{1/2}

Parameters:

• clus1, clus2 – Cluster-specific parameters of the mix. densities

• weights1, weights2 – Weigths of the mixture densities

Returns:

The L^2 distance between p and q

double sample_truncated_beta(double a, double b, double l, double u, std::mt19937 &rng)

Samples from a truncated beta distribution over the interval (l, u) with parameters (a, b) using the inverse-cdf method

Parameters:

• a – the first shape parameter of the beta distribution

• b – the second shape parameter of the beta distribution

• l – lower bound of the truncation interval

• u – upper bound of the truncation interval

• rng – the random number generator

Returns:

a sample from the truncated distribution

Eigen matrix manipulation utilities

The eigen_utils.h file implements a few methods to manipulate groups of matrices, mainly by joining different objects, as well as additional utilities for SPD checking and grid creation.

namespace bayesmix

Functions

Eigen::MatrixXd vstack(const std::vector<Eigen::MatrixXd> &mats)

Concatenates a vector of Eigen matrices along the rows

Parameters:

mats – The matrices to be concatenated

Throws:

std::invalid – argument if sizes mismatch

Returns:

The resulting matrix

void append_by_row(Eigen::MatrixXd *const a, const Eigen::MatrixXd &b)

Concatenates two matrices by row, modifying the first matrix in-place

Throws:

std::invalid_argument – if sizes mismatch

Eigen::MatrixXd append_by_row(const Eigen::MatrixXd &a, const Eigen::MatrixXd &b)

Concatenates two matrices by row

Parameters:

a, b – The matrices to be concatenated

Throws:

std::invalid_argument – if sizes mismatch

Returns:

The resulting matrix

template<template<typename...> class Container>
Eigen::MatrixXd stack_vectors(const Container<Eigen::VectorXd> &rows)

Creates an Eigen matrix from a collection of rows

Template Parameters:

Container – An std-compatible container implementing operator[]

Parameters:

rows – The rows of the matrix

Returns:

The resulting matrix

void check_spd(const Eigen::MatrixXd &mat)

Checks whether the matrix is symmetric and semi-positive definite.

Eigen::MatrixXd get_2d_grid(const double x1, const double x2, const int nx, const double y1, const double y2, const int ny)

Creates a 2d grid over rectangle [x1, x2] x [y1, y2], with nx * ny points

Parameters:

• x1, x2, y1, y2 – Bounds for the rectangle

• nx, ny – Number of points created along the x, y directions

Returns:

The resulting grid

Eigen input-output utilities

The io_utils.h file implements basic input-output utilities for Eigen matrices from and to text files.

namespace bayesmix

Functions

bool check_file_is_writeable(const std::string &filename)

Checks whether the given file is available for writing.

Eigen::MatrixXd read_eigen_matrix(const std::string &filename, const char delim = ',')

Returns an Eigen matrix after reading it from a file.

void write_matrix_to_file(const Eigen::MatrixXd &mat, const std::string &filename, const char delim = ',')

Writes the given Eigen matrix to a text file.

protobuf input-output utilities

The proto_utils.h file implements a few useful functions to manipulate Protobuf objects. For instance, this library implements its own version of vectors and matrices, and the functions implemented here convert from these types to the Eigen ones and viceversa. One can also read a Protobuf from a text file. This is mostly useful for algorithm configuration files.

namespace bayesmix

Functions

void to_proto(const Eigen::VectorXd &vec, bayesmix::Vector *const out)

Writes an Eigen vector to a bayesmix::Vector Protobuf object by pointer.

void to_proto(const Eigen::MatrixXd &mat, bayesmix::Matrix *const out)

Writes an Eigen matrix to a bayesmix::Matrix Protobuf object by pointer.

Eigen::VectorXd to_eigen(const bayesmix::Vector &vec)

Converts a bayesmix::Vector Protobuf object into an Eigen vector.

Eigen::MatrixXd to_eigen(const bayesmix::Matrix &mat)

Converts a bayesmix::Matrix Protobuf object into an Eigen matrix.

void read_proto_from_file(const std::string &filename, google::protobuf::Message *const out)

Writes from a given file to a Protobuf object via pointer.

RNG wrapper

The rng.h file defines a simple Random Number Generation class wrapper. This class wraps the C++ standard RNG object and allows the use of any RNG seed. It is implemented as a singleton, so that every object used in the library has access to the same exact RNG engine. This is needed to ensure that the rng stream is well defined and that every random number generation causes an update in the rng state. The main drawback is that this design does not allow for efficient parallelization, as calls to the Rng::Instance() from different threads could cause data races. A preferred solution would be to define the Rng to be thread-local if omp-parallelism over several cores is desired, see: https://stackoverflow.com/q/64937761

namespace bayesmix

class Rng

#include <rng.h>

Public Functions

inline std::mt19937 &get()

Returns a reference to the underlying RNG object.

inline void seed(const int seed_val)

Sets the RNG seed.

Public Static Functions

static inline Rng &Instance()

Returns (and creates if nonexistent) the singleton of this class.

Private Functions

inline Rng(const int seed_val = 20201103)

inline ~Rng()

Rng(Rng const&) = delete

Rng &operator=(Rng const&) = delete

Private Members

std::mt19937 mt

C++ standard library RNG object.