stacie.dq module

Designed Quadrature (DQ).

This DQ implementation assumes fixed quadrature weights, by default equal weights as in Chebyshev quadrature, but the user can specify any positive weights. Unlike standard Gaussian quadrature, the integration grids are optimal for noisy integrands, for which the main source of error is the variance of the integrand rather than the limited polynomial degree of the quadrature. This is achieved by optimizing the grid point positions assuming equal (or fixed) weights. Equal weights minimize the variance of the numerical integral when all function values on the grid have independent and identically distributed noise. Unlike Monte Carlo integration, the integration grid is deterministic and guarantees a certain level of polynomial degree for smooth integrands. The main contribution of this module is to provide a practical and efficient implementation for the optimization of the grid points. It only supports 1D grids. This module also offers a few additional features that may come in handy for practical applications:

  • One can specify fixed non-equal weights for the quadrature.

  • One can impose symmetry on the weights, such that uneven integrands are integrated exactly.

class Equations(funcs, funcs_d, funcs_dd, targets, weights0, symmetry)[source]

Bases: object

Implementation of the equations to solve for the DQ optimization.

The system of equations contains obviously the constraints that the basis functions must be integrated correctly. Additionally, it contains additional equations that regularize the optimal points to be as uniformly spaced as possible. However, these extra equations are constructed such that their jacobian is always orthogonal to the jacobian of the basis function constraints. This ensures that the regularization does not affect the basis constraints.

Optionally, symmetry can be imposed on the points.

Parameters:
__call__(x, *, deriv=0)[source]

Evaluate the equations to be zerod and optionally their Jacobian.

apply_sym(x)[source]

Return points and weights with their symmetric counterparts if symmetry is imposed.

backprop_sym(jac)[source]

Backpropagate the Jacobian through the symmetry operation.

compute_low(x, w, *, deriv=0)[source]

Compute the equations and optionally their Jacobian for the given points and weights.

funcs: list[Callable]
funcs_d: list[Callable]
funcs_dd: list[Callable]
pen(x, w, *, deriv=0)[source]

Construct the regularization penalty and optionally its gradient and Hessian.

proj(x, w, *, deriv=0)[source]

Construct projection on the orthogonal complement of the derivatives of the basis.

Returns:

  • mat1 – Weighted derivative of basis functions w.r.t. x.

  • apply_proj – A function that applies the projection matrix to a vector: \(P_{i,j}\, v_j\).

  • apply_proj_d – A function that applies the derivative of the projection matrix to a vector: \(\frac{\partial\,P_{i,j}}{\partial\,x_k} v_j\).

solve_weights(points1, weights0, rcond=1e-14)[source]

Assign quadrature weights to exactly match the constraints for the given points.

Parameters:

  • points1 – Fixed (and optimized) grid point positions. It is assumed that the symmetry has already been applied to these points.

  • weights0 – The desired weights for the quadrature, which are used as a regularization to assign the weights.

Returns:

weights1 – The assigned optimal quadrature weights.

symmetry: Symmetry
targets: ndarray[tuple[Any, ...], dtype[_ScalarT]]
weights0: ndarray[tuple[Any, ...], dtype[_ScalarT]]
class Symmetry(*values)[source]

Bases: Enum

Symmetry of the quadrature grid points (and weights) around zero.

NONE = 0
NONZERO = 2
ZERO = 1
construct_dq_empirical(samples, points0, nmoment, symmetry=Symmetry.NONE, weights0=None, rmsdtol=1e-14, maxiter=1000, maxridge=100, verbose=False, do_extra=False)[source]

Construct a DQ grid for an empirical distribution.

Parameters:

  • samples (TypeAliasType) – The samples from the empirical distribution, used to compute the target moments.

  • points0 (int | TypeAliasType) – The number of grid points to optimize or the initial grid points. If symmetry is imposed, only the positive points should be given in points0. If symmetry is ZERO, the zero point is implicitly added, and should not be included in points0.

  • nmoment (int) – The number of moments to match, must be strictly positive and strictly less than npoint.

  • symmetry (Symmetry) – The symmetry of the quadrature grid points (and weights) around zero. If symmetry is imposed, only the positive points are given in points0, and the negative ones are added implicitly by symmetry.

  • weights0 (TypeAliasType | None) – The desired weights for the quadrature, which will be normalized to sum to 1. If not given, all weights are set equal. If symmetry is ZERO, the weight of the zero point must be included.

  • rmsdtol (float) – The convergence threshold for the RMSD of the equations being solved iteratively. The default is very strict.

  • maxiter (int) – The maximum number of iterations for the optimization.

  • maxridge (int) – The maximum number of ridge adjustments for each iteration of the optimization.

  • verbose (bool) – If True, print the optimization progress at each iteration.

  • do_extra (bool) – if True, also return an extra dictionary with additional information about the optimization, which can be used for debugging or analysis.

Return type:

tuple[ndarray[tuple[Any, ...], dtype[TypeVar(_ScalarT, bound= generic)]], ndarray[tuple[Any, ...], dtype[TypeVar(_ScalarT, bound= generic)]]] | tuple[ndarray[tuple[Any, ...], dtype[TypeVar(_ScalarT, bound= generic)]], ndarray[tuple[Any, ...], dtype[TypeVar(_ScalarT, bound= generic)]], dict[str]]

Returns:

  • points – The optimized grid points.

  • weights – The optimized quadrature weights, typically nearly proportional to the Gaussian weights.

  • extra – If do_extra is True, a dictionary containing additional information.

construct_dq_low(points0, funcs, funcs_d, funcs_dd, targets, symmetry=Symmetry.NONE, weights0=None, rmsdtol=1e-14, maxiter=1000, maxridge=100, verbose=False, do_extra=False)[source]

Construct a DQ grid, low-level interface.

Parameters:

  • points0 (TypeAliasType) –

    Initial grid points. The optimization starts from this grid. If symmetry is imposed, only the positive points should be given in points0. If symmetry is ZERO, the zero point is implicitly added, and should not be included in points0.

    Note that the algorithm assumes that the optimal points are pre-conditioned such that they have mean zero and standard deviation one. See construct_dq_empirical for an example of how to do this.

  • funcs (list[Callable]) – The basis functions whose integrals are to be matched by the quadrature. These are functions of a single variable, and are evaluated at the grid points. They must be able to vectorize over the grid points. Note that a constant basis function is automatically added to the list of basis functions.

  • funcs_d (list[Callable]) – The derivatives of the basis functions.

  • funcs_dd (list[Callable]) – The second derivatives of the basis functions.

  • targets (TypeAliasType) – The target integrals of the basis functions, which the quadrature should match.

  • symmetry (Symmetry) – The symmetry of the quadrature grid points (and weights) around zero. If symmetry is imposed, only the positive points are given in points0, and the negative ones are added implicitly by symmetry.

  • weights0 (TypeAliasType | None) – The target weights for the quadrature points. If not given, all weights are set equal. If symmetry is ZERO, the weight of the zero point must be included.

  • rmsdtol (float) – The convergence threshold for the RMSD of the equations being solved iteratively. The default is very strict.

  • maxiter (int) – The maximum number of iterations for the optimization.

  • maxridge (int) – The maximum number of ridge adjustments for each iteration of the optimization.

  • verbose (bool) – If True, print the optimization progress at each iteration.

  • do_extra (bool) – if True, also return an extra dictionary with additional information about the optimization, which can be used for debugging or analysis.

Return type:

tuple[ndarray[tuple[Any, ...], dtype[TypeVar(_ScalarT, bound= generic)]], ndarray[tuple[Any, ...], dtype[TypeVar(_ScalarT, bound= generic)]]] | tuple[ndarray[tuple[Any, ...], dtype[TypeVar(_ScalarT, bound= generic)]], ndarray[tuple[Any, ...], dtype[TypeVar(_ScalarT, bound= generic)]], dict[str]]

Returns:

  • points1 – The optimized grid points. If symmetry is imposed, the full grid is returned, including the negative points.

  • weights1 – The optimized quadrature weights corresponding to the optimized grid points. These weights are typically nearly proportional to the given weights.

  • extra – If do_extra is True, a dictionary containing intermediate results:

    • points0: The initial grid points.

    • weights0: The initial quadrature weights.

    • errors0: The errors of the equations at the initial grid points and weights.

    • jacobian0: The Jacobian of the equations at the initial grid points and weights.

    • funcs0: The basis functions evaluated at the initial grid points.

    • points1: The optimized grid points, same as the first return value.

    • weights1: The optimized quadrature weights.

    • errors1: The errors of the equations at the optimized grid points and weights.

    • jacobian1: The Jacobian of the equations at the optimized grid points and weights.

    • funcs1: The basis functions evaluated at the optimized grid points.

construct_dq_stdnormal(points0, nmoment, symmetry=Symmetry.ZERO, weights0=None, rmsdtol=1e-14, maxiter=1000, maxridge=100, verbose=False, do_extra=False)[source]

Construct a DQ grid for the standard normal distribution.

Parameters:

  • points0 (TypeAliasType) – The optimization starts from this grid. If symmetry is imposed, only the positive points should be given in points0. If symmetry is ZERO, the zero point is implicitly added, and should not be included in points0.

  • nmoment (int) – The number of moments to match, must be strictly positive and even. Only even moments are considered, as the distribution is symmetric.

  • symmetry (Symmetry) – The symmetry of the quadrature grid points (and weights) around zero. For the standard normal distribution, symmetry must be imposed. One can choose whether to include the zero point in the grid or not.

  • weights0 (TypeAliasType | None) – The desired weights for the quadrature. If not given, all weights are set equal. If symmetry is ZERO, the weight of the zero point must be included.

  • rmsdtol (float) – The convergence threshold for the RMSD of the equations being solved iteratively. The default is very strict.

  • maxiter (int) – The maximum number of iterations for the optimization.

  • maxridge (int) – The maximum number of ridge adjustments for each iteration of the optimization.

  • verbose (bool) – If True, print the optimization progress at each iteration.

  • do_extra (bool) – if True, also return an extra dictionary with additional information about the optimization, which can be used for debugging or analysis.

Return type:

tuple[ndarray[tuple[Any, ...], dtype[TypeVar(_ScalarT, bound= generic)]], ndarray[tuple[Any, ...], dtype[TypeVar(_ScalarT, bound= generic)]]] | tuple[ndarray[tuple[Any, ...], dtype[TypeVar(_ScalarT, bound= generic)]], ndarray[tuple[Any, ...], dtype[TypeVar(_ScalarT, bound= generic)]], dict[str]]

Returns:

  • points1 – The optimized grid points, symmetric around zero.

  • weights1 – The optimized quadrature weights, typically matching the given weights.

  • extra – If do_extra is True, a dictionary containing additional information.

dq3(mean, std, skew)[source]

Construct a 3-point quadrature grid to integrate degree-3 polynomials exactly.

The results are exact for an integral of a degree-3 polynomial times a distribution with the given mean, standard deviation and skewness.

Parameters:

  • mean (float) – The mean of the distribution.

  • std (float) – The standard deviation of the distribution.

  • skew (float) – The skewness of the distribution. This value must be in the range \([-1/\sqrt{2}, 1/\sqrt{2}]\), otherwise the quadrature points become complex. A ValueError is raised if the skewness is out of this range.

Return type:

ndarray[tuple[Any, ...], dtype[TypeVar(_ScalarT, bound= generic)]]

Returns:

grid – The 3-point quadrature grid.

plot_dq(extra, plot_dist=None, figsize=None)[source]

Plot the results of the DQ grid optimization.

Parameters:
solve_modified_lm(x, equations, rmsdtol=1e-14, maxiter=1000, maxridge=100, verbose=False)[source]

Solve the equations using a modified Levenberg-Marquardt algorithm.

Parameters:

  • x (ndarray[tuple[Any, ...], dtype[TypeVar(_ScalarT, bound= generic)]]) – The initial guess of the solution.

  • equations (Equations) – The equations to solve, which must be an instance of the Equations class.

  • rmsdtol (float) – The convergence threshold for the RMSD of the equations being solved iteratively.

  • maxiter (int) – The maximum number of iterations for the optimization.

  • maxridge (int) – The maximum number of ridge adjustments for each iteration of the optimization.

  • verbose (bool) – If True, print the optimization progress at each iteration.

Return type:

ndarray[tuple[Any, ...], dtype[TypeVar(_ScalarT, bound= generic)]]

Returns:

x – The optimized solution.