beetroots.sampler.utils package

Submodules

beetroots.sampler.utils.mml module

Marginal maximum likelihood update for scalar regularization parameters.

class beetroots.sampler.utils.mml.EBayesMMLE(scale: float, N0: int | float, N1: int | float, dim: int, vmin: float, vmax: float, homogeneity: float = 1.0, exponent: float = 0.8)

Bases: ABC

Abstract class implementing an empirical Bayes maximum marginal likelihood approach to estimate a scalar regularization parameter.

Empirical Bayes maximum marginal likelihood estimation of a scalar regularization parameter \(\theta\) from within an MCMC algorithm using a projected gradient ascent. The implementation follows the description provided in Vidal et al. [2020].

BayesMMLE constructor.

Parameters:

• scale (float) – Scale parameter involved in the definition of the projected gradient stepsize.

• N0 (Union[int, float]) – Number of iterations defining the initial update phase.

• N1 (Union[int, float]) – Number of iterations for the stabilization phase.

• dim (int) – Dimension of the space over which the distribution \(p(x \mid \theta)\) is defined.

• vmin (float) – Lower limit of the admissible interval.

• vmax (float) – Upper limit of the admissible interval.

• homogeneity (float, optional) – Positive homogeneity factor of the function \(- \log p(x \mid \theta)\) with respect to \(x\). By default 1.

• exponent (float, optional) – Exponent involved in the evolution of the stepsize of the projected gradient. By default 0.8

N0

N1

dim

exponent

homogeneity

mean_theta

mean_theta_old

scale

sum_weights

update(theta: float, iteration: int, gX: float) → float

Update parameter \(\theta\) from iteration \(k\) to \(k+1\).

Parameters:

• theta (float) – Current value of the parameter to be optimized.

• iteration (int) – Current iteration of the sampler.

• gX (float) – Value of the potential function \(g\) evaluated using the current iterate.

Returns:

Updated value theta^{(k+1)}.

Return type:

float

vmax

vmin

class beetroots.sampler.utils.mml.EBayesMMLELogRate(scale: float, N0: int | float, N1: int | float, dim: int, vmin: float = 1e-08, vmax: float = 100000000.0, homogeneity: float = 1.0, exponent: float = 0.8)

Bases: EBayesMMLE

Empirical Bayes maximum marginal likelihood approach to estimate a parameter \(\theta\) in log scale when

\[p(x, y \mid \theta) \propto p(y \mid x) p(x \mid \theta) p(x \mid \theta) = Z^{-1}(\theta) \exp(-\theta g(x))\]

with \(g\) a proper, closed, convex and \(\alpha\)-positively homogeneous function.

BayesMMLE constructor.

Parameters:

• scale (float) – Scale parameter involved in the definition of the projected gradient stepsize.

• N0 (Union[int, float]) – Number of iterations defining the initial update phase.

• N1 (Union[int, float]) – Number of iterations for the stabilization phase.

• dim (int) – Dimension of the space over which the distribution \(p(x \mid \theta)\) is defined.

• vmin (float) – Lower limit of the admissible interval.

• vmax (float) – Upper limit of the admissible interval.

• homogeneity (float, optional) – Positive homogeneity factor of the function \(- \log p(x \mid \theta)\) with respect to \(x\). By default 1.

• exponent (float, optional) – Exponent involved in the evolution of the stepsize of the projected gradient. By default 0.8

N0

N1

dim

exponent

homogeneity

mean_theta

mean_theta_old

scale

sum_weights

vmax

vmin

class beetroots.sampler.utils.mml.EBayesMMLELogScale(scale: float, N0: int | float, N1: int | float, dim: int, vmin: float, vmax: float, homogeneity: float = 1.0, exponent: float = 0.8)

Bases: EBayesMMLE

Empirical Bayes maximum marginal likelihood approach to estimate a parameter \(\theta\) in log scale when

\[p(x, y \mid \theta) \propto p(y \mid x) p(x \mid \theta) p(x \mid \theta) = Z^{-1}(\theta) \exp(-\theta^{-1} g(x))\]

with \(g\) a proper, closed, convex and \(\alpha\)-positively homogeneous function.

BayesMMLE constructor.

Parameters:

• scale (float) – Scale parameter involved in the definition of the projected gradient stepsize.

• N0 (Union[int, float]) – Number of iterations defining the initial update phase.

• N1 (Union[int, float]) – Number of iterations for the stabilization phase.

• dim (int) – Dimension of the space over which the distribution \(p(x \mid \theta)\) is defined.

• vmin (float) – Lower limit of the admissible interval.

• vmax (float) – Upper limit of the admissible interval.

• homogeneity (float, optional) – Positive homogeneity factor of the function \(- \log p(x \mid \theta)\) with respect to \(x\). By default 1.

• exponent (float, optional) – Exponent involved in the evolution of the stepsize of the projected gradient. By default 0.8

N0

N1

dim

exponent

homogeneity

mean_theta

mean_theta_old

scale

sum_weights

vmax

vmin

class beetroots.sampler.utils.mml.EBayesMMLERate(scale: float, N0: int | float, N1: int | float, dim: int, vmin: float, vmax: float, homogeneity: float = 1.0, exponent: float = 0.8)

Bases: EBayesMMLE

Empirical Bayes maximum marginal likelihood approach to estimate a parameter \(\theta\) in linear scale when

\[p(x, y \mid \theta) \propto p(y \mid x) p(x \mid \theta) p(x \mid \theta) = Z^{-1}(\theta) \exp(-\theta g(x))\]

with \(g\) a proper, closed, convex and \(\alpha\)-positively homogeneous function.

BayesMMLE constructor.

Parameters:

• scale (float) – Scale parameter involved in the definition of the projected gradient stepsize.

• N0 (Union[int, float]) – Number of iterations defining the initial update phase.

• N1 (Union[int, float]) – Number of iterations for the stabilization phase.

• dim (int) – Dimension of the space over which the distribution \(p(x \mid \theta)\) is defined.

• vmin (float) – Lower limit of the admissible interval.

• vmax (float) – Upper limit of the admissible interval.

• homogeneity (float, optional) – Positive homogeneity factor of the function \(- \log p(x \mid \theta)\) with respect to \(x\). By default 1.

• exponent (float, optional) – Exponent involved in the evolution of the stepsize of the projected gradient. By default 0.8

N0

N1

dim

exponent

homogeneity

mean_theta

mean_theta_old

scale

sum_weights

vmax

vmin

class beetroots.sampler.utils.mml.EBayesMMLEScale(scale: float, N0: int | float, N1: int | float, dim: int, vmin: float, vmax: float, homogeneity: float = 1.0, exponent: float = 0.8)

Bases: EBayesMMLE

Empirical Bayes maximum marginal likelihood approach to estimate a parameter \(\theta\) in linear scale when

\[p(x, y \mid \theta) \propto p(y \mid x) p(x \mid \theta) p(x \mid \theta) = Z^{-1}(\theta) \exp(-\theta^{-1} g(x))\]

with \(g\) a proper, closed, convex and \(\alpha\)-positively homogeneous function.

BayesMMLE constructor.

Parameters:

• scale (float) – Scale parameter involved in the definition of the projected gradient stepsize.

• N0 (Union[int, float]) – Number of iterations defining the initial update phase.

• N1 (Union[int, float]) – Number of iterations for the stabilization phase.

• dim (int) – Dimension of the space over which the distribution \(p(x \mid \theta)\) is defined.

• vmin (float) – Lower limit of the admissible interval.

• vmax (float) – Upper limit of the admissible interval.

• homogeneity (float, optional) – Positive homogeneity factor of the function \(- \log p(x \mid \theta)\) with respect to \(x\). By default 1.

• exponent (float, optional) – Exponent involved in the evolution of the stepsize of the projected gradient. By default 0.8

N0

N1

dim

exponent

homogeneity

mean_theta

mean_theta_old

scale

sum_weights

vmax

vmin

beetroots.sampler.utils.mml.pgd_rate(theta: float, stepsize: float, vmin: float, vmax: float, dim: int, gX: float, homogeneity: float) → float

Update parameter \(\theta\) from iteration \(k\) to \(k+1\) using projected gradient ascent, assuming \(p(x \mid \theta)\) is of the form

\[p(x \mid \theta) \propto \exp(- \theta g(x)),\]

with \(g \in \Gamma_0(\mathbb{R^N})\).

Parameters:

• theta (float) – Current value of the rate parameter \(\theta^{(k)}\) to be optimized.

• stepsize (float) – Stepsize.

• vmin (float) – Lower limit of the admissible interval.

• vmax (float) – Upper limit of the admissible interval.

• dim (float) –

  Dimension of the space over which the distribution

  \(p(x \mid \theta)\) is defined.

• gX (float) – Value of the potential function \(g\) evaluated using the current iterate.

• homogeneity (float) – Positive homogeneity factor of the function \(- \log p(x \mid \theta)\) with respect to \(x\). By default 1.

Returns:

Updated value \(\theta^{(k+1)}\).

Return type:

float

beetroots.sampler.utils.mml.pgd_rate_log(theta: float, stepsize: float, vmin: float, vmax: float, dim: int, gX: float, homogeneity: float) → float

Update parameter \(\log \theta^{(k)}\) from iteration \(k\) to \(k+1\) using projected gradient ascent (in log scale), assuming \(p(x \mid \theta)\) is of the form

\[p(x \mid \theta) \propto \exp(- \theta g(x)),\]

with \(g \in \Gamma_0(\mathbb{R^N})\).

Parameters:

• theta (float) – Current value of the rate parameter \(\theta^{(k)}\) to be optimized.

• stepsize (float) – Stepsize.

• vmin (float) – Lower limit of the admissible interval for \(\log \theta\).

• vmax (float) – Upper limit of the admissible interval for \(\log \theta\).

• dim (float) – Dimension of the space over which the distribution \(p(x \mid \theta)\) is defined.

• gX (float) – Value of the potential function \(g\) evaluated using the current iterate.

• homogeneity (float) – Positive homogeneity factor of the function \(- \log p(x \mid \theta)\) with respect to \(x\). By default 1.

Returns:

Updated value \(\theta^{(k+1)}\) (in linear scale).

Return type:

float

beetroots.sampler.utils.mml.pgd_scale(theta: float, stepsize: float, vmin: float, vmax: float, dim: int, gX: float, homogeneity: float) → float

Update parameter \(\theta\) from iteration \(k\) to \(k+1\) using projected gradient ascent, assuming \(p(x \mid \theta)\) is of the form

\[p(x \mid \theta) \propto \exp(- g(x) / \theta),\]

with \(g \in \Gamma_0(\mathbb{R^N})\).

Parameters:

• theta (float) – Current value of the scale parameter \(\theta^{(k)}\) to be optimized.

• stepsize (float) – Stepsize.

• vmin (float) – Lower limit of the admissible interval.

• vmax (float) – Upper limit of the admissible interval.

• dim (float) –

  Dimension of the space over which the distribution

  \(p(x \mid \theta)\) is defined.

• gX (float) – Value of the potential function \(g\) evaluated using the current iterate.

• homogeneity (float) – Positive homogeneity factor of the function \(- \log p(x \mid \theta)\) with respect to \(x\). By default 1.

Returns:

Updated value \(\theta^{(k+1)}\).

Return type:

float

beetroots.sampler.utils.mml.pgd_scale_log(theta: float, stepsize: float, vmin: float, vmax: float, dim: int, gX: float, homogeneity: float) → float

Update parameter \(\theta\) from iteration \(k\) to \(k+1\) using projected gradient ascent (in log scale), assuming \(p(x \mid \theta)\) is of the form

\[p(x \mid \theta) \propto \exp(- g(x) / \theta),\]

with \(g \in \Gamma_0(\mathbb{R^N})\).

Parameters:

• theta (float) – Current value of the scale parameter \(\theta^{(k)}\) to be optimized.

• stepsize (float) – Stepsize.

• vmin (float) – Lower limit of the admissible interval.

• vmax (float) – Upper limit of the admissible interval.

• dim (float) –

  Dimension of the space over which the distribution

  \(p(x \mid \theta)\) is defined.

• gX (float) – Value of the potential function \(g\) evaluated using the current iterate.

• homogeneity (float) – Positive homogeneity factor of the function \(- \log p(x \mid \theta)\) with respect to \(x\). By default 1.

Returns:

Updated value \(\theta^{(k+1)}\) (in linear scale).

Return type:

float

beetroots.sampler.utils.my_sampler_params module

Defines a Dataclass that stores the parameters of the augmented PSGLD sampler

class beetroots.sampler.utils.my_sampler_params.MySamplerParams(initial_step_size: float, extreme_grad: float, history_weight: float, selection_probas: ndarray | List[float], k_mtm: int, is_stochastic: bool = True, compute_correction_term: bool = True)

Bases: object

Dataclass that stores the parameters of the sampler proposed in Palud et al. [2023]

Parameters:

• initial_step_size (float) – step size used in the Position-dependent MALA transition kernel, denoted \(\epsilon\) in the article

• extreme_grad (float) – limit value that avoids division by zero when computing the RMSProp preconditioner, denoted \(\eta\) in the article

• history_weight (float) – weight of past values of \(v\) in the exponential decay (cf RMSProp preconditioner), denoted \(\alpha\) in the article

• selection_probas (np.ndarray of shape (2,)) – vector of selection probabilities for the MTM and PMALA kernels, respectively, i.e., \([p_{MTM}, 1 - p_{MTM}]\)

• k_mtm (int) – number of candidates in the MTM kernel, denoted \(K\) in the article

• is_stochastic (bool) – if True, the algorithm performs sampling, and optimization otherwise, by default True

• compute_correction_term (bool) – wether or not to use the correction term (denoted \(\gamma\) in the article) during the sampling (only used if is_stochastic=True), by default True

compute_correction_term

extreme_grad

history_weight

initial_step_size

is_stochastic

k_mtm

selection_probas

beetroots.sampler.utils.utils module

beetroots.sampler.utils.utils.compute_nlpdf_spatial_proposal(candidates: ndarray, spatial_list_edges: ndarray, spatial_weights: ndarray, idx_pix: ndarray) → ndarray

evaluates the negative log ratio of prior proposal in the MTM-chromatic Gibbs kernel in MySampler

Note

The function works in four steps :

• step 1: get the values of the neighbors of each pixel in idx_pix using get_neighboring_pixels

• step 2: compute the vector distances between the candidate and the pixel neighbors

• step 3: compute list of sums of subsets using compute_sum_subsets_norms

• step 4: evaluate LogSumExp to obtain the log pdf of the Gaussian mixture, using the stable numba implementation.

Parameters:

• candidates (np.ndarray) – iterate

• spatial_list_edges (np.ndarray) – list of edges, defined as pairs of pixel indices. Generated by the class:SpatialPrior.

• spatial_weights (np.ndarray) – regularization weights of the spatial prior

• idx_pix (np.ndarray of shape (n_pix,)) – array of indices of pixels for which candidates are to be proposed

Returns:

negative log ratio of prior proposal

Return type:

np.ndarray of shape (n_pix, k_mtm)

beetroots.sampler.utils.utils.compute_sum_subsets_mean(dists: ndarray) → Tuple[ndarray, ndarray]

Computes the sums of the elements of each subset of neighbors, including the empty set, ordered such that the last element is the complement of the first, and so on. It uses bitwise operation and take advantage of the tree structure of the problem.

Parameters:

dists (np.ndarray of shape (N_neighbors, k_mtm, D)) – Squared distances from a point to each of its neighbors.

Returns:

• np.ndarray of shape ((2**N_neighbors, k_mtm, D)) – Array containing the sum of distances for all subsets of neighbors, including the empty set.

• np.ndarray – Array containing the weight of each subset of neighbors

beetroots.sampler.utils.utils.get_neighboring_pixels(current_Theta: ndarray, list_edges: ndarray, idx_pix: int) → ndarray

returns the array of current values of the neighboring pixels of a specified pixel

Parameters:

• current_Theta (np.ndarray) – current iterate

• list_edges (np.ndarray) – list of edges, defined as pairs of pixel indices. Generated by the class:SpatialPrior.

• idx_pix (int) – index of a pixel

Returns:

array of current values of the neighboring pixels

Return type:

np.ndarray

beetroots.sampler.utils.utils.numba_logsumexp_stable(x: ndarray, weights: ndarray) → float

stable numba implementation of the logsumexp function. Used in step 4 of compute_nlratio_prior_proposal.

Parameters:

x (np.ndarray) – array of real values

Returns:

result of the logsumexp function

Return type:

float

beetroots.sampler.utils.utils.sample_conditional_spatial_and_indicator_prior(current_Theta: ndarray, spatial_list_edges: ndarray, spatial_weights: ndarray, indicator_lower_bounds: ndarray, indicator_upper_bounds: ndarray, indicator_indicator_margin_scale: float, idx_pix: ndarray, k_mtm: int, seed: int) → ndarray

draw k_mtm iid samples for each candidate in the ìdx_pix` array.

Parameters:

• current_Theta (np.ndarray) – current iterate

• spatial_list_edges (np.ndarray) – list of edges, defined as pairs of pixel indices. Generated by the class:SpatialPrior.

• spatial_weights (np.ndarray) – regularization weights of the spatial prior

• indicator_lower_bounds (np.ndarray) – vector of lower bounds for the D physical parameters

• indicator_upper_bounds (np.ndarray) – vector of upper bounds for the D physical parameters

• indicator_indicator_margin_scale (float) – scaling parameter, quantifies how much values out of the validity intervals are penalized

• idx_pix (np.ndarray) – array of indices of pixels for which candidates are to be proposed

• k_mtm (int) – number of candidates in the MTM kernel

• seed (int) – random seed, used for sampling reproducibility

Returns:

samples drawn from the conditional prior

Return type:

np.ndarray

beetroots.sampler.utils.utils.sample_conditional_spatial_prior(current_Theta: ndarray, spatial_list_edges: ndarray, spatial_weights: ndarray, idx_pix: ndarray, k_mtm: int, seed: int) → ndarray

draw k_mtm iid samples for each candidate in the ìdx_pix` array.

Parameters:

• current_Theta (np.ndarray) – current iterate

• spatial_list_edges (np.ndarray) – list of edges, defined as pairs of pixel indices. Generated by the class:SpatialPrior.

• spatial_weights (np.ndarray) – regularization weights of the spatial prior

• indicator_lower_bounds (np.ndarray) – vector of lower bounds for the D physical parameters

• indicator_upper_bounds (np.ndarray) – vector of upper bounds for the D physical parameters

• indicator_indicator_margin_scale (float) – scaling parameter, quantifies how much values out of the validity intervals are penalized

• idx_pix (np.ndarray) – array of indices of pixels for which candidates are to be proposed

• k_mtm (int) – number of candidates in the MTM kernel

• seed (int) – random seed, used for sampling reproducibility

Returns:

samples drawn from the conditional prior

Return type:

np.ndarray

beetroots.sampler.utils.utils.sample_generalized_gaussian(alpha: float, a_alpha: float, size: Tuple[int, int], seed: int) → ndarray

sample from a generalized Gaussian distribution of pdf:

\[p(\theta) = \frac{2}{\delta \Gamma(1/4)} \exp \left\{- \left( \frac{\theta}{\delta}\right)^4\right\}\]

with here \(\delta = 1/A(\alpha)\). See Nardon and Pianca [2009] for more details.

Parameters:

• alpha (float) – power used in the distribution. alpha=2 is equivalent to the Gaussian distribution.

• a_alpha (float) – scaling parameter

• size (Tuple[int, int]) – number of samples to draw

• seed (int) – random seed for sampling reproducibility

Returns:

array of samples

Return type:

np.ndarray

beetroots.sampler.utils.utils.sample_smooth_indicator(lower_bounds: ndarray, upper_bounds: ndarray, indicator_margin_scale: float, size: Tuple[int, int], seed: int) → ndarray

algorithm to draw iid samples from the smooth indicator prior, detailed in Appendix B in Palud et al. [2023]

Parameters:

• lower_bounds (np.ndarray) – vector of lower bounds for each of the D physical parameters

• upper_bounds (np.ndarray) – vector of upper bounds for each of the D physical parameters

• indicator_margin_scale (float) – scaling parameter, quantifies how much values out of the validity intervals are penalized

• size (Tuple[int, int]) – number of samples to draw

• seed (int) – random seed, for sampling reproducibility

Returns:

iid samples drawn from the smooth uniform distribution

Return type:

np.ndarray

Module contents