# Multivariate Jacobi ensemble¶

Important

For the details please refer to:

1. the extensive documentation of MultivariateJacobiOPE below

2. the associated Jupyter notebook which showcases MultivariateJacobiOPE

3. our NeurIPS‘19 paper [GBV19] On two ways to use determinantal point processes for Monte Carlo integration

The figures below display a sample of a $$d=2$$ dimensional Jacobi ensemble MultivariateJacobiOPE with $$N=200$$ points. The red and green dashed curves correspond to the normalized base densities proportional to $$(1-x)^{a_1} (1+x)^{b_1}$$ and $$(1-y)^{a_2} (1+y)^{b_2}$$, respectively.

import numpy as np
import matplotlib.pyplot as plt
from dppy.multivariate_jacobi_ope import MultivariateJacobiOPE

# The .plot() method outputs smtg only in dimension d=1 or 2

# Number of points / dimension
N, d = 200, 2
# Jacobi parameters in [-0.5, 0.5]^{d x 2}
jac_params = np.array([[0.5, 0.5],
[-0.3, 0.4]])

dpp = MultivariateJacobiOPE(N, jac_params)

# Get an exact sample
sampl = dpp.sample()

# Display
# the cloud of points
# the base probability densities
# the marginal empirical histograms
dpp.plot(sample=sampl, weighted=False)
plt.tight_layout()

# Attach a weight 1/K(x,x) to each of the points
dpp.plot(sample=sampl, weighted=True)
plt.tight_layout()


Fig. 42 (png, hires.png, pdf)

Fig. 43 (png, hires.png, pdf)

• In the first plot, observe that the empirical marginal density converges to the arcsine density $$\frac{1}{\pi\sqrt{1-x^2}}$$, displayed in orange.

• In the second plot, we take the same sample and attach a weight $$\frac{1}{K(x,x)}$$ to each of the points. This illustrates the choice of the weights defining the estimator of [BH16] as a proxy for the reference measure.

Implementation of the class MultivariateJacobiOPE used in [GBV19] for Monte Carlo with Determinantal Point Processes

It has 3 main methods:

class dppy.multivariate_jacobi_ope.MultivariateJacobiOPE(N, jacobi_params)[source]

Bases: object

Multivariate Jacobi Orthogonal Polynomial Ensemble used in [GBV19] for Monte Carlo with Determinantal Point Processes

This corresponds to a continuous multivariate projection DPP with state space $$[-1, 1]^d$$ with respect to

• reference measure $$\mu(dx) = w(x) dx$$ (see also eval_w()), where

$w(x) = \prod_{i=1}^{d} (1-x)^{a_i} (1+x)^{b_i}$
• kernel $$K$$ (see also K())

$K(x, y) = \sum_{\mathfrak{b}(k)=0}^{N-1} \frac{P_{k}(x)P_{k}(y)} {\|P_{k}\|^2}$

where

• $$k \in \mathbb{N}^d$$ is a multi-index ordered according to the ordering $$\mathfrak{b}$$ (see compute_ordering_BaHa16())

• $$P_{k}(x) = \prod_{i=1}^d P_{k_i}^{(a_i, b_i)}(x_i)$$ is the product of orthogonal Jacobi polynomials w.r.t. $$\mu$$, (see eval_poly_multiD())

$\int_{-1}^{1} P_{k}^{(a_i,b_i)}(x) P_{\ell}^{(a_i,b_i)}(x) (1-x)^{a_i} (1+x)^{b_i} d x \propto \delta_{k\ell}$
Parameters
• N (int) – Number of points $$N \geq 1$$

• jacobi_params (array_like) – Jacobi parameters $$[(a_i, b_i)]_{i=1}^d$$ The number of rows $$d$$ prescribes the ambient dimension of the points i.e. $$x_{1}, \dots, x_{N} \in [-1, 1]^d$$. - when $$d=1$$, $$a_1, b_1 > -1$$ - when $$d \geq 2$$, $$|a_i|, |b_i| \leq \frac{1}{2}$$

• Multivariate Jacobi ensemble

• when $$d=1$$, the univariate Jacobi ensemble is sampled by computing the eigenvalues of a properly randomized tridiagonal matrix of [KN04]

• [BH16] initiated the use of the multivariate Jacobi ensemble for Monte Carlo integration. In particular, they proved CLT with variance decay of order $$N^{-(1+1/d)}$$ which is faster that the $$N^{-1}$$ rate of vanilla Monte Carlo where the points are drawn i.i.d. from the base measure.

K(X, Y=None)[source]

Evaluate the orthogonal projection kernel $$K$$. It is based on the 3-terms recurrence relations satisfied by each univariate orthogonal Jacobi polynomial, see also SciPy eval_jacobi()

$K(x, y) = \sum_{\mathfrak{b}(k)=0}^{N-1} \frac{P_{k}(x)P_{k}(y)} {\|P_{k}\|^2}$

where

• $$k \in \mathbb{N}^d$$ is a multi-index ordered according to the ordering $$\mathfrak{b}$$, compute_ordering_BaHa16()

• $$P_{k}(x) = \prod_{i=1}^d P_{k_i}^{(a_i, b_i)}(x_i)$$ is the product of orthogonal Jacobi polynomials w.r.t. $$\mu(dx) = \prod_{i=1}^{d} (1-x_i)^{a_i} (1+x_i)^{b_i} d x_i$$

$\int_{-1}^{1} P_{k}^{(a_i,b_i)}(x) P_{\ell}^{(a_i,b_i)}(x) (1-x)^{a_i} (1+x)^{b_i} d x \propto \delta_{k\ell}$
Parameters
• X (array_like) – Array of points $$\in [-1, 1]^d$$, with size $$N\times d$$ where $$N$$ is the number of points

• Y (array_like (default None)) – Array of points $$\in [-1, 1]^d$$, with size $$N\times d$$ where $$N$$ is the number of points

Returns

Pointwise evaluation of $$K$$ as depicted in the following pseudo code output

• if Y is None

• K(X, X) if X.size $$=d$$

• [K(x, x) for x in X] otherwise

• otherwise

• K(X, Y) if X.size=Y.size$$=d$$

• [K(X, y) for y in Y] if X.size $$=d$$

• [K(x, y) for x, y in zip(X, Y)] otherwise

Return type

• float if Y is None and X.size $$=d$$

• array_like otherwise

eval_poly_multiD(X, normalize='norm')[source]

Evaluate (and potentially normalize) multivariate Jacobi polynomials $$P_{k}(x) = \prod_{i=1}^d P_{k_i}^{(a_i, b_i)}(x_i)$$

Parameters
• X (array_like) – Array of points $$\in [-1, 1]^d$$, with size $$N\times d$$ where $$N$$ is the number of points

• normalize (str (default 'norm')) –

• ‘norm’

• ’square_norm’

Returns

• normalize='norm' $$P_k(X) / \left\| P_k \right\|$$

• normalize='square_norm' $$P_k(X) / \left\| P_k \right\|^2$$

Return type

array_like

eval_w(X, jac_params=None)[source]

Evaluate $$w(x) = \prod_{i=1}^{d} (1-x_i)^{a_i} (1+x_i)^{b_i}$$ which corresponds to the density of the base measure $$\mu$$ if jac_params is None.

Parameters
• X (array_like) – Array of points $$\in [-1, 1]^d$$, with size $$N\times d$$ where $$N$$ is the number of points

• jac_params (array_like (default None)) –

• if None, use attribute MultivariateJacobiOPE.jacobi_params i.e. $$[(a_i, b_i)]_{i=1}^d \in [-0.5, 0.5]^{d \times 2}$$.

• else Jacobi parameters $$[(a_i, b_i)]_{i=1}^d \in [-1, \infty)^{d \times 2}$$.

Returns

• if jac_params is None, evaluation of $$w(x)$$ the density of the base measure

• else jac_params $$= [(a_i, b_i)]_{i=1}^d$$ evaluation of $$w(x) = \prod_{i=1}^{d} (1-x_i)^{a_i} (1+x_i)^{b_i}$$

Return type

array_like

sample(nb_trials_max=10000, random_state=None, tridiag_1D=True)[source]

Use the chain rule [HKPVirag06] (Algorithm 18) to sample $$\left(x_{1}, \dots, x_{N} \right)$$ with density

$\begin{split}& \frac{1}{N!} \left[K(x_n,x_p)\right]_{n,p=1}^{N} \prod_{n=1}^{N} w(x_n)\\ &= \frac{1}{N} K(x_1,x_1) w(x_1) \prod_{n=2}^{N} \frac{ K(x_n,x_n) - K(x_n,x_{1:n-1}) \left[\left[K(x_k,x_l)\right]_{k,l=1}^{n-1}\right]^{-1} K(x_{1:n-1},x_n) }{N-(n-1)} w(x_n)\end{split}$

The order in which the points were sampled can be forgotten to obtain a valid sample of the corresponding DPP

Each conditional density is sampled using rejection sampling with proposal density $$\frac{1}{N} K(x,x) w(x)$$

• $$x_1 \sim \frac{1}{N} K(x,x) w(x)$$ using sample_proposal_lev_score()

• $$x_n | Y = x_{1}, \dots, x_{n-1}$$

$\frac{1}{N-|Y|} [K(x,x) - K(x, Y) K_Y^{-1} K(Y, x)] w(x) \leq \frac{N}{N-|Y|} \frac{1}{N} K(x,x) w(x)$
sample_proposal_lev_score(nb_trials_max=10000, random_state=None)[source]

Use a rejection sampling mechanism to sample

$\frac{1}{N} K(x, x) w(x) dx = \frac{1}{N} \sum_{\mathfrak{b}(k)=0}^{N-1} \left( \frac{P_k(x)}{\left\| P_k \right\|} \right)^2 w(x)$

with proposal distribution

$w_{eq}(x) d x = \prod_{i=1}^{d} \frac{1}{\pi\sqrt{1-(x_i)^2}} d x_i$

Since the target density is a mixture, we can sample from it by

1. Select a multi-index $$k$$ uniformly at random in $$\left\{ \mathfrak{b}^{-1}(0), \dots, \mathfrak{b}^{-1}(N-1) \right\}$$

2. Sample from $$\left( \frac{P_k(x)}{\left\| P_k \right\|} \right)^2 w(x) dx$$ with proposal $$w_{eq}(x) d x$$.

The acceptance ratio writes

$\frac{ \left( \frac{P_k(x)}{\left\| P_k \right\|} \right)^2 w(x)} {w_{eq}(x)} = \prod_{i=1}^{d} \pi \left( \frac{P_{k_i}^{(a_i, b_i)}(x)} {\left\| P_{k_i}^{(a_i, b_i)} \right\|} \right)^2 (1-x_i)^{a_i+\frac{1}{2}} (1+x_i)^{b_i+\frac{1}{2}} \leq C_{k}$

which can be bounded using the result of [Gau09] on Jacobi polynomials.

Note

Each of the rejection constant $$C_{k}$$ is computed at initialization of the MultivariateJacobiOPE object using compute_Gautschi_bounds()

Returns

A sample $$x\in[-1,1]^d$$ with probability distribution $$\frac{1}{N} K(x,x) w(x)$$

Return type

array_like

dppy.multivariate_jacobi_ope.compute_Gautschi_bounds(jacobi_params, ordering, log_scale=True)[source]

Compute the rejection constants for the acceptance/rejection mechanism used in sample_proposal_lev_score() to sample

$\frac{1}{N} K(x, x) w(x) dx = \frac{1}{N} \sum_{\mathfrak{b}(k)=0}^{N-1} \left( \frac{P_k(x)}{\left\| P_k \right\|} \right)^2 w(x)$

with proposal distribution

$w_{eq}(x) d x = \prod_{i=1}^{d} \frac{1}{\pi\sqrt{1-(x_i)^2}} d x_i$

To get a sample:

1. Draw a multi-index $$k$$ uniformly at random in $$\left\{ \mathfrak{b}^{-1}(0), \dots, \mathfrak{b}^{-1}(N-1) \right\}$$

2. Sample from $$P_k(x)^2 w(x) dx$$ with proposal $$w_{eq}(x) d x$$.

The acceptance ratio writes

$\frac{\left( \frac{P_k(x)}{\left\| P_k \right\|} \right)^2 w(x)} {w_{eq}(x)} = \prod_{i=1}^{d} \pi \left( \frac{P_{k_i}^{(a_i, b_i)}(x)} {\left\| P_{k_i}^{(a_i, b_i)} \right\|} \right)^2 (1-x_i)^{a_i+\frac{1}{2}} (1+x_i)^{b_i+\frac{1}{2}} \leq C_k$
• For $$k_i>0$$ we use a result on Jacobi polynomials given by, e.g., [Gau09], for $$\quad|a|,|b| \leq \frac{1}{2}$$

$\begin{split}& \pi (1-x)^{a+\frac{1}{2}} (1+x)^{b+\frac{1}{2}} \left( \frac{P_{n}^{(a, b)}(x)} {\left\| P_{n}^{(a, b)} \right\|} \right)^2\\ &\leq \frac{2} {n!(n+(a+b+1) / 2)^{2 \max(a,b)}} \frac{\Gamma(n+a+b+1) \Gamma(n+\max(a,b)+1)} {\Gamma(n+\min(a,b)+1)}\end{split}$
• For $$k_i=0$$, we use less involved properties of the Jacobi polynomials:

• $$P_{0}^{(a, b)} = 1$$

• $$\|P_{0}^{(a, b)}\|^2 = 2^{a+b+1} \operatorname{B}(a+1,b+1)$$

• $$m = \frac{b-a}{a+b+1}$$ is the mode of $$(1-x)^{a+\frac{1}{2}} (1+x)^{b+\frac{1}{2}}$$ (valid since $$a+\frac{1}{2}, b+\frac{1}{2} > 0$$)

So that,

$\begin{split} \pi (1-x)^{a+\frac{1}{2}} (1+x)^{b+\frac{1}{2}} \left(\frac{P_{0}^{(a, b)}(x)} {\|P_{0}^{(a, b)}\|}\right)^{2} &= \frac {\pi (1-x)^{a+\frac{1}{2}} (1+x)^{b+\frac{1}{2}}} {\|P_{0}^{(a, b)}\|^2} \\ &\leq \frac {\pi (1-m)^{a+\frac{1}{2}} (1+m)^{b+\frac{1}{2}}} {2^{a+b+1} \operatorname{B}(a+1,b+1)}\end{split}$
Parameters
• jacobi_params (array_like) –

Jacobi parameters $$[(a_i, b_i)]_{i=1}^d \in [-\frac{1}{2}, \frac{1}{2}]^{d \times 2}$$.

The number of rows $$d$$ prescribes the ambient dimension of the points i.e. $$x_{1}, \dots, x_{N} \in [-1, 1]^d$$

• ordering (array_like) –

Ordering of the multi-indices $$\in\mathbb{N}^d$$ defining the order between the multivariate monomials (see also compute_ordering_BaHa16())

• the number of rows corresponds to the number $$N$$ of monomials considered.

• the number of columns $$=d$$

• log_scale (bool) – If True, the rejection bound is computed using the logarithmic versions betaln, gammaln of beta and gamma functions to avoid overflows

Returns

The rejection bounds $$C_{k}$$ for $$k = \mathfrak{b}^{-1}(0), \dots, \mathfrak{b}^{-1}(N-1)$$

Return type

array_like

dppy.multivariate_jacobi_ope.compute_ordering_BaHa16(N, d)[source]

Compute the ordering of the multi-indices $$\in\mathbb{N}^d$$ defining the order between the multivariate monomials as described in Section 2.1.3 of [BH16].

Parameters
• N (int) – Number of polynomials $$(P_k)$$ considered to build the kernel K() (number of points of the corresponding MultivariateJacobiOPE)

• d (int) – Size of the multi-indices $$k\in \mathbb{N}^d$$ characterizing the _degree_ of $$P_k$$ (ambient dimension of the points x_{1}, dots, x_{N} in [-1, 1]^d)

Returns

Array of size $$N\times d$$ containing the first $$N$$ multi-indices $$\in\mathbb{N}^d$$ in the order prescribed by the ordering $$\mathfrak{b}$$ [BH16] Section 2.1.3

Return type

array_like

For instance, for $$N=12, d=2$$

[(0, 0), (0, 1), (1, 0), (1, 1), (0, 2), (1, 2), (2, 0), (2, 1), (2, 2), (0, 3), (1, 3), (2, 3)]


dppy.multivariate_jacobi_ope.compute_poly1D_square_norms(jacobi_params, deg_max)[source]

Compute the square norms $$\|P_{k}^{(a_i,b_i)}\|^2$$ of each (univariate) Jacobi polynomial for $$k=0$$ to deg_max and $$a_i, b_i =$$ jacobi_params[i, :] Recall that the Jacobi polynomials $$\left( P_{k}^{(a_i,b_i)} \right)$$ are orthogonal w.r.t. $$(1-x)^{a_i} (1+x)^{b_i}$$.

$\|P_{k}^{(a_i,b_i)}\|^2 = \int_{-1}^{1} \left( P_{k}^{(a_i,b_i)}(x) \right)^2 (1-x)^{a_i} (1+x)^{b_i} d x = \frac{2^{a_i+b_i+1}} {2k+a_i+b_i+1} \frac{\Gamma(k+a_i+1)\Gamma(k+b_i+1)} {\Gamma(k+a_i+b_i+1)n!}$
Parameters
• jacobi_params (array_like) – Jacobi parameters $$[(a_i, b_i)]_{i=1}^d \in [-\frac{1}{2}, \frac{1}{2}]^{d \times 2}$$ The number of rows $$d$$ prescribes the ambient dimension of the points i.e. $$x_{1}, \dots, x_{N} \in [-1, 1]^d$$

• deg_max (int) – Maximal degree of 1D Jacobi polynomials

Returns

Array of size deg_max + 1 $$\times d$$ with entry $$k,i$$ given by $$\|P_{k}^{(a_i,b_i)}\|^2$$

Return type

array_like

dppy.multivariate_jacobi_ope.poly_1D_degrees(max_degrees)[source]

:return:poly_1D_degrees[i, j] = i if i <= max_degrees[j] else 0