Computing the center-outward quantile/superquantile function¶

The following codes are intended to help practicioners to test the entropic map as an empirical center-outward quantile function.

Without going any further, we state here the definitions we use. For more details, we refer to our paper defining center-outward superquantiles.

The entropic map can be seen as a regularized version of the center-outward quantile map. It is defined from $\pi_\epsilon$, the solution of Kantorovich regularized OT problem. $\pi_\epsilon$ is a joint probability measure for $(U,X)$ if $U\sim U_d$ and $X \sim \nu$. By barycentric projection, one can define, for any $u\in\mathcal{B}(0,1)$,

\begin{equation} Q_{\epsilon}(u) = \mathbb{E}_{\pi_\epsilon}\big( X \vert U = u \big), \end{equation} and \begin{equation} S_{\epsilon}(u) = \frac{1}{1-\Vert u \Vert} \int_{\Vert u \Vert}^1 Q_\epsilon\Big(t\frac{u}{\Vert u \Vert}\Big) dt. \end{equation}

These are our regularized versions of the center-outward quantile and superquantile functions respectively.

In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()   

def quantile(eps,x,Y,v):
    '''
    Returns the entropic map at point x with regularization parameter eps.
    v is the optimal potential obtained with the Robbins-Monro algorithm for semi-discrete EOT 
    '''
    arg = (v - cost(x, Y))/eps
    to_sum = np.exp(arg-np.max(arg))
    weights = to_sum/(np.sum(to_sum))
    return(np.sum(Y*weights.reshape(len(weights),1),axis=0))

def superquantile(eps,x,Y,v,N=100): 
    '''
    Returns the superquantile map at point x from the "quantile" function defined just before
    N = number of points in the Riemann sum
    '''
    SQx = 0
    tau = np.linalg.norm(x) 
    x = x/tau # norm x w.r.t. Ud
    t = np.linspace(tau,1,N) # subdivision of the interval [norm(x),1]
    for i in range(N):
        Qx = quantile(eps,t[i]*x,Y,v)
        SQx = SQx + Qx
    SQx = SQx / N
    return(SQx)

def cost(x,Y):
    """
    Squared euclidean distance between x and the J points of Y.
    if x is the ith point of the support of the distribution, 
    the result is the ith line of the cost matrix.
    """
    diff = Y-x
    return(0.5*np.linalg.norm(diff, axis = 1)**2)

def c_transform(v, x, Y, nu, eps=0):
    """Calculate the c_transform of g"""
    if eps > 0:
        arg = (v - cost(x, Y))/eps
        M = np.max(arg)
        to_sum = np.exp(arg-M)*nu
        return(-eps*np.log(np.sum(to_sum)) - eps*M )
    else:
        cost_x = cost(x,Y)
        return(np.min(cost_x-v))
    
def grad_heps(v, x, Y, nu, eps=0):
    """
    Calculate the gradient h_eps 
    """
    if eps > 0:
        arg = (v-cost(x,Y))/eps
        M = np.max(arg)
        pi = nu * np.exp(arg-M)
        pi = pi/pi.sum()
        return(nu-pi)
    else:
        cost_x = cost(x,Y)
        j_star = np.argmin(cost_x - v)
        to_return = nu.copy()
        to_return[j_star] = nu[j_star] - 1
        return(to_return)

def h_eps(v, x, Y, nu, eps=0):
    """
    Calculate the function h_eps whose expectation equals H_eps.
    """
    return(np.sum(v*nu)+c_transform(v,x,Y,nu,eps)-eps)

def Robbins_Monro_Algo(Y, eps=1, n_iter=10000):
    """
    robbins-monro algorithm for solving EOT
    """    
    c = 0.51
    gamma = 1
    
    J = len(Y)
    nu = 1/J*np.ones(J)
    
    # Storage of recursive estimates.
    W_hat_storage = np.zeros(n_iter)
    h_eps_storage = np.zeros(n_iter)
    #bar_W_hat_storage = np.zeros(n_iter)
    #bar_h_eps_storage = np.zeros(n_iter)
    
    # Initialisation of dual vector v
    v = np.random.random(J)
    v = v - np.mean(v)
    
    # First sample to initiate
    Z = np.random.normal(0,1,2)
    Z = Z / np.sqrt(np.sum(Z**2))
    r = np.random.random()
    x_0 = r*Z
    
    W_hat_storage[0] = h_eps(v, x_0, Y, nu, eps)
    h_eps_storage[0] = h_eps(v, x_0, Y, nu, eps)

    # Robbins Monro algorithm
    for k in range(1,n_iter):
        #if k % 2000 == 0:
        #    print(k)
        # Sample from mu
        Z = np.random.normal(0,1,2)
        Z = Z / np.sqrt(np.sum(Z**2))
        r = np.random.random()
        x = r*Z
    
        # Update v
        v = v + gamma/((k+1)**c) * grad_heps(v, x, Y, nu, eps) 
        #bar_v = ((k-1)/k)*bar_v + v/k
        #bar_h_eps_storage[k] = h_eps(bar_v, x, Y, nu, eps)
        
        # Storage of h_eps at point (x,g).
        h_eps_storage[k] = h_eps(v, x, Y, nu, eps)
    
        # approximation of Sinkhorn divergence
        W_hat_storage[k] = k/(k+1) * W_hat_storage[k-1] + 1/(k+1) * h_eps_storage[k]
        # bar_W_hat_storage[k] = k/(k+1) * bar_W_hat_storage[k-1] + 1/(k+1) * bar_h_eps_storage[k]
        
    L = [v, W_hat_storage]
    return(L)

Simulate a sample Y in $\mathbb{R}^d$¶

When reusing, this is where you should change the codes, by inserting your data.¶

It must be of size (n_samples,d).

In [ ]:
from sklearn.datasets import make_moons
In [ ]:
X = make_moons(n_samples=1000, noise=0.1)[0]

mean = np.array([0.5,0.25])
mat = np.array([[0.8,0.2],[0.2,0.5]])/4
cov = np.dot(mat,mat)
Y = np.random.multivariate_normal(mean,cov,500)


plt.figure(figsize=(5,5))
Y = np.concatenate((X,Y))
Y = Y - np.mean(Y,axis=0)
plt.scatter(Y.T[0],Y.T[1])
plt.title("Very-non-convex data")
Out[ ]:
Text(0.5, 1.0, 'Very-non-convex data')
No description has been provided for this image

1. Computing $Q_\epsilon$ for any $u \in \mathcal{B}(0,1)$.¶

When using these codes, solving OT is a first step, computing the quantile/superquantile map is the second step.

STEP 1 : solving EOT¶

One can solve OT using any method from the POT library.

In this notebook, we used the stochastic gradient descent corresponding to Bercu, Bigot, (2021), Asymptotic distribution and convergence rates of stochastic algorithms for entropic optimal transportation between probability measures, The Annals of Statistics.

In [ ]:
eps = 0.001 # Regularization parameter
v1 = Robbins_Monro_Algo(Y,eps=eps, n_iter=30000)[0]

STEP 2 : compute quantiles and superquantiles¶

After that, for any $u \in \mathcal{B}(0,1)$, computing its image by the center-outward quantile map and/or superquantile map is immediate.

The quantile/superquantile function of some u depends on

  • the data data Y,
  • the second entropic Kantorovich potential v1,
  • and the parameter eps.

2. Descriptive plots associated with Y in dimension d=2¶

The next cell details each step of the computation of descriptive plots associated with the center-outward quantile and superquantile functions.

For the sake of simplicity, one can reuse these codes by using the function defined at the end of this notebook, just after.

In [ ]:
eps = 0.001
rayons = np.array([0.1,0.3,0.5,0.7,0.9])

##################################
##################################
# STEP 1 : solving EOT

v1 = Robbins_Monro_Algo(Y,eps=eps, n_iter=30000)[0]

##################################
##################################
# STEP 2 : Transporting Ud with the MK quantile map and superquantile map
# (with the barycentric projection as a regularized quantile map)

sphere = []
Ud = []
for i in range(2000):
    Z = np.random.normal(0,1,2)
    Z = Z / np.sqrt(np.sum(Z**2))
    r = np.random.random()
    Ud.append(r*Z)
    sphere.append(Z)
Ud = np.array(Ud)
sphere = np.array(sphere)

# 1st plot
plt.figure(figsize=(7,7))
plt.scatter(Ud.T[0],Ud.T[1],label="$U_d$")
N = 15
N2 = 20
angles = 2*np.pi*np.arange(N)/N
for k in range(N):
    theta = angles[k]
    r = np.arange(N2)/N2
    us = np.array([[r*np.cos(theta)][0],[r*np.sin(theta)][0]])
    plt.plot(us[0],us[1],color='red')

for k in range(len(rayons)):
    tau = rayons[k]
    contour = tau*sphere
    plt.scatter(contour.T[0],contour.T[1],color="blue",s=1)

##################################
# we compute the regularized quantile contours associated with entropic optimal transport
##################################

n = 1000

contours_quantiles = []
contours_superquantiles = []
for tau in rayons:
    CQ1 = []
    CSQ1 = []
    normCQ1 = np.zeros(n)
    normCSQ1 = np.zeros(n)

    for i in range(n):
        # simulate uniformly on the sphere 
        Z = np.random.normal(0,1,2)
        Z = Z / np.sqrt(np.sum(Z**2)) 
        u = tau * Z
        # Contour quantile of order tau for Y1
        TZ = quantile(eps,u,Y,v1)
        CQ1.append( TZ )
        # Contour superquantile of order tau for Y1
        TZ = superquantile(eps,u,Y,v1)
        CSQ1.append( TZ  )

    CQ1 = np.array(CQ1)
    CSQ1 = np.array(CSQ1)
    contours_quantiles.append(CQ1)
    contours_superquantiles.append(CSQ1)
    
contours_superquantiles = np.array(contours_superquantiles)
contours_quantiles = np.array(contours_quantiles)

##################################
# after that, it's all about vizualisation stuff
##################################

# 2nd plot
plt.figure(figsize=(7,7))
plt.scatter(Y.T[0],Y.T[1],label = "Y")

for k in range(len(rayons)):
    contour =  contours_quantiles[k] 
    plt.scatter(contour.T[0],contour.T[1],s=1,color = "blue")

signcurve = np.zeros((N2,2))
for k in range(N):
    theta = angles[k]
    r = np.arange(N2)/N2
    for j in range(N2):
        u = np.array([r[j]*np.cos(theta),r[j]*np.sin(theta)])
        TZ = quantile(eps,u,Y,v1)
        signcurve[j,:] = TZ
    plt.plot(signcurve.T[0],signcurve.T[1],color='red')

# 3rd plot
plt.figure(figsize=(7,7))
plt.scatter(Y.T[0],Y.T[1],label = "Y")
for k in range(len(rayons)):
    contour = contours_superquantiles[k]
    plt.scatter(contour.T[0],contour.T[1],s=1,color = "blue")

signcurve = np.zeros((N2,2))
for k in range(N):
    theta = angles[k]
    r = np.arange(N2)/N2
    for j in range(N2):
        u = np.array([r[j]*np.cos(theta),r[j]*np.sin(theta)])
        TZ = quantile(eps,u,Y,v1)
        signcurve[j,:] = TZ
    plt.plot(signcurve.T[0],signcurve.T[1],color='red')
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

A single function for descriptive plots¶

In dimension $d=2$, the next function can be used to draw these descriptive plots for any scatterplot $Y$ of size (n_samples,2).

In [ ]:
def descriptiv_plots(Y,superplot=True,Udplot=True):
    d = 2
    eps = 0.001
    rayons = np.array([0.1,0.3,0.5,0.7,0.9])

    ##################################
    ##################################
    # STEP 1 : solving EOT
    v1 = Robbins_Monro_Algo(Y,eps=eps, n_iter=30000)[0]

    ##################################
    ##################################
    # STEP 2 : Transporting Ud with the MK quantile map and superquantile map
    # (with the barycentric projection as a regularized quantile map)

    sphere = []
    Ud = []
    for i in range(2000):
        Z = np.random.normal(0,1,d)
        Z = Z / np.sqrt(np.sum(Z**2))
        r = np.random.random()
        Ud.append(r*Z)
        sphere.append(Z)
    Ud = np.array(Ud)
    sphere = np.array(sphere)
    N = 15
    N2 = 20
    angles = 2*np.pi*np.arange(N)/N
        
    if (Udplot==True):
        plt.figure(figsize=(7,7))
        plt.scatter(Ud.T[0],Ud.T[1],label="$U_d$")
        for k in range(N):
            theta = angles[k]
            r = np.arange(N2)/N2
            us = np.array([[r*np.cos(theta)][0],[r*np.sin(theta)][0]])
            plt.plot(us[0],us[1],color='red')

        for k in range(len(rayons)):
            tau = rayons[k]
            contour = tau*sphere
            plt.scatter(contour.T[0],contour.T[1],color="blue",s=1)

    ##################################
    # we compute the regularized quantile contours associated with entropic optimal transport
    ##################################

    n = 1000

    contours_quantiles = []
    contours_superquantiles = []
    for tau in rayons:
        CQ1 = []
        CSQ1 = []

        for i in range(n):
            # simulate uniformly on the sphere 
            Z = np.random.normal(0,1,2)
            Z = Z / np.sqrt(np.sum(Z**2)) 
            u = tau * Z
            # Contour quantile of order tau for Y1
            TZ = quantile(eps,u,Y,v1)
            CQ1.append( TZ )
            # Contour superquantile of order tau for Y1
            TZ = superquantile(eps,u,Y,v1)
            CSQ1.append( TZ  )

        CQ1 = np.array(CQ1)
        CSQ1 = np.array(CSQ1)
        contours_quantiles.append(CQ1)
        contours_superquantiles.append(CSQ1)

    contours_superquantiles = np.array(contours_superquantiles)
    contours_quantiles = np.array(contours_quantiles)

    ##################################
    # after that, it's all about vizualisation stuff
    ##################################

    plt.figure(figsize=(7,7))
    plt.scatter((Y).T[0],(Y).T[1],label = "Y")

    for k in range(len(rayons)):
        contour = ( contours_quantiles[k] )
        plt.scatter(contour.T[0],contour.T[1],s=1,color = "blue")

    signcurve = np.zeros((N2,2))
    for k in range(N):
        theta = angles[k]
        r = np.arange(N2)/N2
        for j in range(N2):
            u = np.array([r[j]*np.cos(theta),r[j]*np.sin(theta)])
            TZ = quantile(eps,u,Y,v1)
            signcurve[j,:] = TZ
        plt.plot(signcurve.T[0],signcurve.T[1],color='red')

    if (superplot==True):
        plt.figure(figsize=(7,7))
        plt.scatter((Y).T[0],(Y).T[1],label = "Y")
        for k in range(len(rayons)):
            contour = contours_superquantiles[k] 
            plt.scatter(contour.T[0],contour.T[1],s=1,color = "blue")


        signcurve = np.zeros((N2,2))
        for k in range(N):
            theta = angles[k]
            r = np.arange(N2)/N2
            for j in range(N2):
                u = np.array([r[j]*np.cos(theta),r[j]*np.sin(theta)])
                TZ = quantile(eps,u,Y,v1)
                signcurve[j,:] = TZ
            plt.plot(signcurve.T[0],signcurve.T[1],color='red')

By default, the function draws the three descriptive plots.

If superplot=False, the one with superquantiles contours is hidden.

If Udplot=False, the one with the reference distribution is hidden.

By using the boxplot function from matplotlib, and by selecting the positions argument accordingly to your data, one can combine univariate and bivariate boxplots through a single vizualisation.

In [ ]:
descriptiv_plots(Y,superplot=False,Udplot=False)
plt.boxplot(Y.T[0],vert=False,positions=[-1])
plt.boxplot(Y.T[1],positions=[-2])
Out[ ]:
{'whiskers': [<matplotlib.lines.Line2D at 0x175539060>,
  <matplotlib.lines.Line2D at 0x175539300>],
 'caps': [<matplotlib.lines.Line2D at 0x1755395a0>,
  <matplotlib.lines.Line2D at 0x175539840>],
 'boxes': [<matplotlib.lines.Line2D at 0x175538dc0>],
 'medians': [<matplotlib.lines.Line2D at 0x175539ae0>],
 'fliers': [<matplotlib.lines.Line2D at 0x175539d80>],
 'means': []}
No description has been provided for this image