{ "cells": [ { "cell_type": "markdown", "id": "1e9639fd", "metadata": {}, "source": [ "# Distancias de Wasserstein\n", "\n", "Este cuaderno acompaña el **Capítulo 10** de las notas del curso (*Distancias de Wasserstein*). Exploramos computacionalmente:\n", "\n", "- El cálculo de $W_p$ en **dimensión uno** mediante las pseudoinversas $F^{[-1]}$, comparado con la solución del programa lineal.\n", "- La **fórmula cerrada para gaussianas** y su verificación por Monte Carlo.\n", "- La comparación entre **distintos exponentes**: $W_q \\le W_p$ para $q \\le p$ y la desigualdad recíproca en soportes acotados.\n", "- La relación entre **convergencia en $W_p$, convergencia débil y momentos**, con el ejemplo de la masa que se escapa.\n", "- La **interpolación por desplazamiento** (geodésicas de McCann) frente a la interpolación lineal.\n", "- La **dualidad de Kantorovich–Rubinstein** para $W_1$.\n", "\n", "Utilizamos la biblioteca [POT (Python Optimal Transport)](https://pythonot.github.io/) para resolver el problema de Kantorovich exacto, y `scipy` para los programas lineales de la última sección." ] }, { "cell_type": "code", "execution_count": null, "id": "19562288", "metadata": {}, "outputs": [], "source": [ "# @title\n", "pip install POT" ] }, { "cell_type": "code", "execution_count": null, "id": "ebf89e63", "metadata": {}, "outputs": [], "source": [ "# @title\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import ot\n", "from scipy.stats import norm\n", "from scipy.optimize import linprog\n", "from scipy.linalg import sqrtm\n", "\n", "rng = np.random.default_rng(0)" ] }, { "cell_type": "markdown", "id": "ca7702dd", "metadata": {}, "source": [ "## 1. $W_p$ en dimensión uno\n", "\n", "Recordemos la definición: para $\\mu,\\nu\\in\\mathcal P_p(\\mathbb R^d)$,\n", "\n", "$$\n", "W_p(\\mu,\\nu)=\\Bigl(\\min_{\\pi\\in\\Pi(\\mu,\\nu)}\\int |x-y|^p\\,d\\pi\\Bigr)^{1/p}.\n", "$$\n", "\n", "En dimensión uno el plan comonótono es óptimo para todo costo $h(x-y)$ con $h$ convexa, y en particular para $|x-y|^p$ con $p\\ge1$. Esto da la fórmula\n", "\n", "$$\n", "W_p(\\mu,\\nu)^p=\\int_0^1\\bigl|F_\\mu^{[-1]}(t)-F_\\nu^{[-1]}(t)\\bigr|^p\\,dt ,\n", "$$\n", "\n", "que reduce el cálculo de $W_p$ a una integral en $[0,1]$. Para medidas empíricas con $n$ átomos de igual peso, la pseudoinversa es constante a trozos y la integral es una suma sobre los datos **ordenados**:\n", "\n", "$$\n", "W_p(\\mu_n,\\nu_n)^p=\\frac1n\\sum_{k=1}^n |x_{(k)}-y_{(k)}|^p .\n", "$$\n", "\n", "Verificamos esto contra la solución del programa lineal que calcula POT (`ot.emd2`), que no sabe nada de la estructura unidimensional." ] }, { "cell_type": "code", "execution_count": null, "id": "a510136c", "metadata": {}, "outputs": [], "source": [ "def Wp_1d(x, y, p):\n", " \"\"\"W_p entre las medidas empíricas uniformes de las muestras x e y (mismo tamaño), vía cuantiles.\"\"\"\n", " xs, ys = np.sort(x), np.sort(y)\n", " return np.mean(np.abs(xs - ys)**p)**(1/p)\n", "\n", "def Wp_lp(x, y, p, a=None, b=None):\n", " \"\"\"W_p vía el programa lineal (POT). Funciona en cualquier dimensión.\"\"\"\n", " x = np.atleast_2d(x.T).T if x.ndim == 1 else x\n", " y = np.atleast_2d(y.T).T if y.ndim == 1 else y\n", " a = np.ones(len(x))/len(x) if a is None else a\n", " b = np.ones(len(y))/len(y) if b is None else b\n", " M = ot.dist(x, y, metric='euclidean')**p\n", " return ot.emd2(a, b, M, numItermax=10_000_000)**(1/p)\n", "\n", "n = 300\n", "x = rng.normal(0, 1, n) # muestra de N(0,1)\n", "y = rng.exponential(1.5, n) - 1 # muestra de Exp(1/1.5) desplazada\n", "\n", "print(f\"{'p':>3} {'cuantiles':>12} {'LP (POT)':>12}\")\n", "for p in [1, 2, 3, 5]:\n", " print(f\"{p:>3} {Wp_1d(x, y, p):>12.6f} {Wp_lp(x, y, p):>12.6f}\")" ] }, { "cell_type": "markdown", "id": "228e060e", "metadata": {}, "source": [ "Las dos columnas coinciden hasta la precisión del solver. El cálculo por cuantiles tiene costo $O(n\\log n)$ (ordenar); el programa lineal, en el peor caso, $O(n^3\\log n)$. Ilustramos la fórmula gráficamente: $W_1$ es el área entre las dos pseudoinversas, que coincide con el área entre las dos funciones de distribución." ] }, { "cell_type": "code", "execution_count": null, "id": "4100caa1", "metadata": {}, "outputs": [], "source": [ "t = (np.arange(n) + 0.5)/n\n", "xs, ys = np.sort(x), np.sort(y)\n", "\n", "fig, ax = plt.subplots(1, 2, figsize=(11, 4))\n", "ax[0].step(t, xs, where='mid', label=r'$F_\\mu^{[-1]}$')\n", "ax[0].step(t, ys, where='mid', label=r'$F_\\nu^{[-1]}$')\n", "ax[0].fill_between(t, xs, ys, alpha=0.25, step='mid')\n", "ax[0].set_xlabel('t'); ax[0].set_title(r'$W_1$ = área entre las pseudoinversas'); ax[0].legend()\n", "\n", "grid = np.linspace(min(xs.min(), ys.min()), max(xs.max(), ys.max()), 800)\n", "Fx = np.searchsorted(xs, grid, side='right')/n\n", "Fy = np.searchsorted(ys, grid, side='right')/n\n", "ax[1].plot(grid, Fx, label=r'$F_\\mu$'); ax[1].plot(grid, Fy, label=r'$F_\\nu$')\n", "ax[1].fill_between(grid, Fx, Fy, alpha=0.25)\n", "ax[1].set_xlabel('x'); ax[1].set_title(r'... = área entre las distribuciones'); ax[1].legend()\n", "plt.tight_layout(); plt.show()\n", "\n", "print(\"W_1 por cuantiles :\", Wp_1d(x, y, 1))\n", "print(\"W_1 = ∫|F_mu - F_nu| dx:\", np.trapezoid(np.abs(Fx - Fy), grid))" ] }, { "cell_type": "markdown", "id": "7f7dd771", "metadata": {}, "source": [ "**Ejercicio.** La igualdad $W_1(\\mu,\\nu)=\\int_{\\mathbb R}|F_\\mu(x)-F_\\nu(x)|\\,dx$ vale sólo para $p=1$. Verificar numéricamente que $\\int|F_\\mu-F_\\nu|^2\\,dx$ **no** es $W_2^2$, y explicar geométricamente por qué (las áreas entre las curvas coinciden, pero las integrales de las potencias no: en un caso se integra en $t$ y en el otro en $x$)." ] }, { "cell_type": "markdown", "id": "7ff9c103", "metadata": {}, "source": [ "## 2. Gaussianas\n", "\n", "Para $\\mu=N(m_0,\\Sigma_0)$ y $\\nu=N(m_1,\\Sigma_1)$ en $\\mathbb R^d$ vale\n", "\n", "$$\n", "W_2(\\mu,\\nu)^2=|m_0-m_1|^2+\\operatorname{tr}\\Bigl(\\Sigma_0+\\Sigma_1-2\\bigl(\\Sigma_0^{1/2}\\Sigma_1\\Sigma_0^{1/2}\\bigr)^{1/2}\\Bigr),\n", "$$\n", "\n", "y el mapa óptimo es afín, $T(x)=m_1+A(x-m_0)$ con $A=\\Sigma_0^{-1/2}(\\Sigma_0^{1/2}\\Sigma_1\\Sigma_0^{1/2})^{1/2}\\Sigma_0^{-1/2}$. Implementamos la fórmula y la comparamos con:\n", "\n", "1. la función `ot.gaussian.bures_wasserstein_distance` de POT (misma fórmula, otra implementación);\n", "2. una estimación **Monte Carlo**: $W_2$ entre medidas empíricas de $n$ muestras de cada gaussiana, resuelta por programación lineal.\n", "\n", "La estimación Monte Carlo converge a $W_2(\\mu,\\nu)$ cuando $n\\to\\infty$, pero lentamente: la distancia entre una medida y su medida empírica en $\\mathbb R^d$ decae como $n^{-1/d}$ para $d\\ge3$ (y como $n^{-1/2}$, con correcciones logarítmicas, en dimensión baja)." ] }, { "cell_type": "code", "execution_count": null, "id": "ff6949e3", "metadata": {}, "outputs": [], "source": [ "def W2_gauss(m0, S0, m1, S1):\n", " S0h = np.real(sqrtm(S0))\n", " C = np.real(sqrtm(S0h @ S1 @ S0h))\n", " return np.sqrt(np.sum((m0 - m1)**2) + np.trace(S0 + S1 - 2*C))\n", "\n", "def mapa_gauss(m0, S0, m1, S1):\n", " \"\"\"Matriz A del mapa óptimo T(x) = m1 + A (x - m0).\"\"\"\n", " S0h = np.real(sqrtm(S0)); S0hi = np.linalg.inv(S0h)\n", " return S0hi @ np.real(sqrtm(S0h @ S1 @ S0h)) @ S0hi\n", "\n", "m0, S0 = np.array([0., 0.]), np.array([[1.0, 0.3], [0.3, 0.5]])\n", "m1, S1 = np.array([3., 1.]), np.array([[0.6, -0.4], [-0.4, 1.5]])\n", "\n", "w_formula = W2_gauss(m0, S0, m1, S1)\n", "w_pot = ot.gaussian.bures_wasserstein_distance(m0, m1, S0, S1)\n", "print(\"W_2 fórmula de la traza :\", w_formula)\n", "print(\"W_2 POT (Bures) :\", w_pot)\n", "\n", "A = mapa_gauss(m0, S0, m1, S1)\n", "print(\"\\nVerificación A Σ0 A = Σ1 :\\n\", A @ S0 @ A)\n", "print(\"A simétrica:\", np.allclose(A, A.T), \" definida positiva:\", np.all(np.linalg.eigvalsh(A) > 0))" ] }, { "cell_type": "code", "execution_count": null, "id": "f1c98ab1", "metadata": {}, "outputs": [], "source": [ "ns = [50, 100, 200, 500, 1000, 2000]\n", "reps = 5\n", "est = np.zeros((len(ns), reps))\n", "for i, n in enumerate(ns):\n", " for r in range(reps):\n", " X = rng.multivariate_normal(m0, S0, n)\n", " Y = rng.multivariate_normal(m1, S1, n)\n", " est[i, r] = Wp_lp(X, Y, 2)\n", "\n", "plt.figure(figsize=(6, 4))\n", "plt.errorbar(ns, est.mean(1), yerr=est.std(1), fmt='o-', capsize=3, label=r'$W_2(\\mu_n,\\nu_n)$ Monte Carlo')\n", "plt.axhline(w_formula, color='k', ls='--', label='fórmula cerrada')\n", "plt.xscale('log'); plt.xlabel('n muestras'); plt.ylabel(r'$W_2$'); plt.legend(); plt.title('Convergencia de la estimación empírica')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "731ec24d", "metadata": {}, "source": [ "Notar que la estimación empírica **sobreestima** sistemáticamente la distancia para $n$ pequeño: la medida empírica está lejos de la gaussiana que la genera, y ese error se suma. Esto no es un defecto del solver sino un fenómeno estadístico, que en las aplicaciones motiva las regularizaciones entrópicas del último capítulo.\n", "\n", "Visualizamos el mapa óptimo afín $T$ sobre una muestra: los puntos de $\\mu$ y sus imágenes $T(x)$, que se distribuyen según $\\nu$." ] }, { "cell_type": "code", "execution_count": null, "id": "fb9c1ac1", "metadata": {}, "outputs": [], "source": [ "X = rng.multivariate_normal(m0, S0, 400)\n", "TX = m1 + (X - m0) @ A.T\n", "Y = rng.multivariate_normal(m1, S1, 400)\n", "\n", "plt.figure(figsize=(6, 5))\n", "plt.scatter(*X.T, s=8, alpha=0.6, label=r'$x\\sim\\mu$')\n", "plt.scatter(*Y.T, s=8, alpha=0.3, color='gray', label=r'$y\\sim\\nu$ (independiente)')\n", "plt.scatter(*TX.T, s=8, alpha=0.6, color='C3', label=r'$T(x)$')\n", "for k in range(0, 400, 20):\n", " plt.plot([X[k, 0], TX[k, 0]], [X[k, 1], TX[k, 1]], color='C3', lw=0.5, alpha=0.6)\n", "plt.axis('equal'); plt.legend(); plt.title(r'Mapa óptimo afín $T(x)=m_1+A(x-m_0)$'); plt.show()" ] }, { "cell_type": "markdown", "id": "000bfee4", "metadata": {}, "source": [ "## 3. Comparación entre exponentes\n", "\n", "Por la desigualdad de Jensen (o de Hölder), si $1\\le q\\le p$ entonces\n", "\n", "$$\n", "W_q(\\mu,\\nu)\\le W_p(\\mu,\\nu).\n", "$$\n", "\n", "En la dirección opuesta no hay una desigualdad general (piénsese en las medidas con masa que se escapa de la sección siguiente), pero si ambas medidas tienen soporte en un conjunto de diámetro $D$, entonces $|x-y|^p\\le D^{p-q}|x-y|^q$ sobre el soporte de cualquier plan, y por lo tanto\n", "\n", "$$\n", "W_p(\\mu,\\nu)^p\\le D^{\\,p-q}\\,W_q(\\mu,\\nu)^q .\n", "$$\n", "\n", "Verificamos ambas desigualdades con medidas discretas aleatorias en $[0,1]^2$ (diámetro $D=\\sqrt2$)." ] }, { "cell_type": "code", "execution_count": null, "id": "387632bd", "metadata": {}, "outputs": [], "source": [ "def W_discreta(X, a, Y, b, p):\n", " M = ot.dist(X, Y, metric='euclidean')**p\n", " return ot.emd2(a, b, M, numItermax=10_000_000)**(1/p)\n", "\n", "n, m = 40, 60\n", "X, Y = rng.random((n, 2)), rng.random((m, 2))\n", "a = rng.dirichlet(np.ones(n)); b = rng.dirichlet(np.ones(m))\n", "D = np.sqrt(2)\n", "\n", "ps = [1, 1.5, 2, 3, 4, 6]\n", "W = {p: W_discreta(X, a, Y, b, p) for p in ps}\n", "print(\"p -> W_p\"); [print(f\"{p:<4} {W[p]:.5f}\") for p in ps]\n", "\n", "print(\"\\nMonotonía W_q <= W_p (q\", W[p]**p <= D**(p-1)*W[1] + 1e-12)" ] }, { "cell_type": "markdown", "id": "c64545c0", "metadata": {}, "source": [ "## 4. Convergencia en $W_p$, convergencia débil y momentos\n", "\n", "El teorema central del capítulo dice que, en $\\mathcal P_p(\\mathbb R^d)$,\n", "\n", "$$\n", "W_p(\\mu_n,\\mu)\\to0\n", "\\quad\\Longleftrightarrow\\quad\n", "\\mu_n\\rightharpoonup\\mu\\ \\text{ débilmente}\\ \\ \\text{y}\\ \\ \\int|x|^p\\,d\\mu_n\\to\\int|x|^p\\,d\\mu .\n", "$$\n", "\n", "La convergencia débil sola **no** alcanza. El ejemplo canónico es una masa pequeña que se escapa al infinito:\n", "\n", "$$\n", "\\mu_n=\\Bigl(1-\\frac1n\\Bigr)\\delta_0+\\frac1n\\,\\delta_{a_n},\\qquad a_n\\to\\infty .\n", "$$\n", "\n", "Siempre $\\mu_n\\rightharpoonup\\delta_0$ (la masa en $a_n$ tiende a cero). Pero un cálculo directo da $W_p(\\mu_n,\\delta_0)^p=\\frac1n a_n^p$, de modo que la convergencia en $W_p$ depende de la velocidad de escape **y del exponente**:\n", "\n", "- con $a_n=n$: $W_1(\\mu_n,\\delta_0)=1$ no tiende a cero, y $W_2^2=n\\to\\infty$;\n", "- con $a_n=\\sqrt n$: $W_1=n^{-1/2}\\to0$ pero $W_2=1$ no tiende a cero;\n", "- con $a_n=n^{1/4}$: $W_1\\to0$ y $W_2\\to0$ pero $W_4=1$.\n", "\n", "En cada caso lo que falla es exactamente la convergencia del momento de orden $p$: $\\int|x|^p\\,d\\mu_n=\\frac1na_n^p$, que debería tender a $\\int|x|^p\\,d\\delta_0=0$. Verificamos con POT que el programa lineal reproduce estos valores." ] }, { "cell_type": "code", "execution_count": null, "id": "9a98123e", "metadata": {}, "outputs": [], "source": [ "def W_escape(n, a_n, p):\n", " \"\"\"W_p entre mu_n = (1-1/n) d_0 + (1/n) d_{a_n} y delta_0, vía POT.\"\"\"\n", " X = np.array([[0.], [a_n]]); a = np.array([1 - 1/n, 1/n])\n", " Y = np.array([[0.]]); b = np.array([1.])\n", " M = ot.dist(X, Y, metric='euclidean')**p\n", " return ot.emd2(a, b, M)**(1/p)\n", "\n", "ns = np.array([2, 5, 10, 20, 50, 100, 200, 500, 1000])\n", "fig, ax = plt.subplots(1, 3, figsize=(13, 3.8), sharey=True)\n", "for k, (nombre, f) in enumerate([(\"a_n = n\", lambda n: n), (\"a_n = √n\", np.sqrt), (\"a_n = n^{1/4}\", lambda n: n**0.25)]):\n", " for p in [1, 2, 4]:\n", " ax[k].plot(ns, [W_escape(n, f(n), p) for n in ns], 'o-', label=f'$W_{p}$')\n", " ax[k].set_xscale('log'); ax[k].set_yscale('log'); ax[k].set_title(f'${nombre}$'); ax[k].set_xlabel('n')\n", " ax[k].axhline(1, color='k', ls=':', lw=0.8)\n", "ax[0].legend(); ax[0].set_ylabel(r'$W_p(\\mu_n,\\delta_0)$')\n", "plt.suptitle(r'$\\mu_n\\rightharpoonup\\delta_0$ siempre; la convergencia en $W_p$ depende de $p$ y de la velocidad de escape')\n", "plt.tight_layout(); plt.show()" ] }, { "cell_type": "markdown", "id": "58f0223b", "metadata": {}, "source": [ "En el otro sentido, cuando **sí** hay convergencia débil con soportes uniformemente acotados (una de las hipótesis bajo las que demostramos la recíproca), la convergencia en $W_p$ es automática para todo $p$. Un ejemplo: la medida uniforme discreta sobre los puntos medios de la grilla, $\\mu_n=\\frac1n\\sum_{k=1}^n\\delta_{(k-1/2)/n}$, converge a la uniforme en $[0,1]$, y en dimensión uno podemos calcular $W_p$ exactamente con los cuantiles: $|F_n^{[-1]}(t)-t|\\le\\frac1{2n}$ para todo $t$, luego $W_p(\\mu_n,U[0,1])\\le\\frac1{2n}$ para todo $p$, y de hecho $W_p^p=\\int_0^1|F_n^{[-1]}(t)-t|^p\\,dt=\\frac{1}{(p+1)(2n)^p}$." ] }, { "cell_type": "code", "execution_count": null, "id": "c6b40e65", "metadata": {}, "outputs": [], "source": [ "def Wp_grilla_vs_uniforme(n, p, K=200_001):\n", " t = np.linspace(0, 1, K)[1:-1]\n", " Finv_n = (np.floor(t*n) + 0.5)/n # pseudoinversa de la uniforme en {(k-1/2)/n}: escalera\n", " Finv_U = t # pseudoinversa de U[0,1]\n", " return np.mean(np.abs(Finv_n - Finv_U)**p)**(1/p)\n", "\n", "print(f\"{'n':>6} {'W_1 numérico':>14} {'exacto':>10} {'W_2 numérico':>14} {'exacto':>10}\")\n", "for n in [5, 10, 50, 100, 500]:\n", " e1 = 1/(2*2*n); e2 = (1/(3*(2*n)**2))**0.5\n", " print(f\"{n:>6} {Wp_grilla_vs_uniforme(n,1):>14.6f} {e1:>10.6f} {Wp_grilla_vs_uniforme(n,2):>14.6f} {e2:>10.6f}\")" ] }, { "cell_type": "markdown", "id": "18cc9efd", "metadata": {}, "source": [ "## 5. Interpolación por desplazamiento\n", "\n", "Si $\\pi$ es un plan óptimo para $W_2$ entre $\\mu$ y $\\nu$, la curva\n", "\n", "$$\n", "\\mu_t=(e_t)_\\#\\pi,\\qquad e_t(x,y)=(1-t)x+ty,\\qquad t\\in[0,1],\n", "$$\n", "\n", "es una **geodésica** en $(\\mathcal P_2,W_2)$: $W_2(\\mu_s,\\mu_t)=|t-s|\\,W_2(\\mu,\\nu)$. Cuando el plan es inducido por un mapa $T$ (Brenier), $\\mu_t=\\bigl((1-t)\\,\\mathrm{id}+tT\\bigr)_\\#\\mu$: **cada partícula se mueve en línea recta** de $x$ a $T(x)$ a velocidad constante. Esto es lo que McCann llamó interpolación por desplazamiento, y contrasta con la interpolación lineal $(1-t)\\mu+t\\nu$, en la que la masa no se mueve: simplemente se desvanece en un lado y aparece en el otro.\n", "\n", "### 5.1 En dimensión uno\n", "\n", "Comparamos las dos interpolaciones entre dos mezclas de gaussianas. En dimensión uno $T=T_{\\mathrm{mon}}$, y para muestras de igual tamaño el mapa monótono simplemente empareja los datos ordenados." ] }, { "cell_type": "code", "execution_count": null, "id": "2d110f98", "metadata": {}, "outputs": [], "source": [ "n = 4000\n", "mu_s = np.concatenate([rng.normal(-3, 0.5, n//2), rng.normal(-1, 0.4, n//2)])\n", "nu_s = np.concatenate([rng.normal(2, 0.7, n//4), rng.normal(4, 0.3, 3*n//4)])\n", "xs, ys = np.sort(mu_s), np.sort(nu_s) # T_mon: x_(k) -> y_(k)\n", "\n", "ts = [0, 0.25, 0.5, 0.75, 1]\n", "bins = np.linspace(-5, 6, 120)\n", "fig, ax = plt.subplots(2, len(ts), figsize=(15, 5), sharex=True, sharey=True)\n", "for k, t in enumerate(ts):\n", " # desplazamiento: partícula k está en (1-t) x_(k) + t y_(k)\n", " ax[0, k].hist((1-t)*xs + t*ys, bins=bins, density=True, color='C3', alpha=0.8)\n", " ax[0, k].set_title(f't = {t}')\n", " # lineal: mezcla de mu y nu con pesos (1-t), t\n", " mezcla = np.concatenate([rng.choice(mu_s, int((1-t)*n)), rng.choice(nu_s, n - int((1-t)*n))])\n", " ax[1, k].hist(mezcla, bins=bins, density=True, color='C0', alpha=0.8)\n", "ax[0, 0].set_ylabel('desplazamiento\\n' + r'$((1-t)\\,\\mathrm{id}+tT)_\\#\\mu$'); ax[1, 0].set_ylabel('lineal\\n' + r'$(1-t)\\mu+t\\nu$')\n", "plt.tight_layout(); plt.show()" ] }, { "cell_type": "markdown", "id": "e91148a1", "metadata": {}, "source": [ "En la fila superior la masa **viaja**; en la inferior se desvanece y reaparece. Verificamos ahora la propiedad geodésica: $W_2(\\mu_s,\\mu_t)=|t-s|\\,W_2(\\mu,\\nu)$, mientras que la interpolación lineal **no** es geodésica: la distancia $W_2\\bigl((1-s)\\mu+s\\nu,(1-t)\\mu+t\\nu\\bigr)$ se comporta aquí como $\\sqrt{|t-s|}$ (hay que mover una fracción $|t-s|$ de la masa desde el soporte de $\\mu$ hasta el de $\\nu$, que están separados, y eso cuesta del orden de $|t-s|$ en costo cuadrático), que para $|t-s|$ pequeño es mucho **mayor** que $|t-s|\\,W_2(\\mu,\\nu)$. Para medidas discretas con soportes disjuntos, como estas muestras, la curva lineal tiene por eso longitud infinita (Ejercicio 3); para densidades suaves y positivas la longitud es finita, aunque la curva sigue sin ser geodésica." ] }, { "cell_type": "code", "execution_count": null, "id": "5bbd18d8", "metadata": {}, "outputs": [], "source": [ "W_total = Wp_1d(xs, ys, 2)\n", "print(f\"W_2(mu, nu) = {W_total:.4f}\\n\")\n", "print(f\"{'(s,t)':>12} {'|t-s| W_2':>10} {'desplaz.':>10} {'lineal':>10}\")\n", "sub = rng.choice(n, 800, replace=False) # submuestra para el LP de la interpolación lineal\n", "for s, t in [(0, 0.5), (0.25, 0.75), (0.5, 1), (0.2, 0.3)]:\n", " d_desp = Wp_1d((1-s)*xs + s*ys, (1-t)*xs + t*ys, 2)\n", " # lineal: medidas con pesos (no equiponderadas) sobre la unión de soportes -> LP\n", " Xu = np.concatenate([mu_s[sub], nu_s[sub]])\n", " ws = np.concatenate([np.full(800, (1-s)/800), np.full(800, s/800)])\n", " wt = np.concatenate([np.full(800, (1-t)/800), np.full(800, t/800)])\n", " M = ot.dist(Xu[:, None], Xu[:, None]) # métrica por defecto: euclídea al cuadrado\n", " d_lin = np.sqrt(ot.emd2(ws, wt, M, numItermax=10_000_000))\n", " print(f\"({s:.2f},{t:.2f}) {abs(t-s)*W_total:>10.4f} {d_desp:>10.4f} {d_lin:>10.4f}\")" ] }, { "cell_type": "markdown", "id": "30420009", "metadata": {}, "source": [ "### 5.2 Gaussianas en el plano\n", "\n", "Entre gaussianas la interpolación por desplazamiento permanece gaussiana: como $T(x)=m_1+A(x-m_0)$ es afín, $\\mu_t=N(m_t,\\Sigma_t)$ con\n", "\n", "$$\n", "m_t=(1-t)m_0+tm_1,\\qquad \\Sigma_t=\\bigl((1-t)I+tA\\bigr)\\,\\Sigma_0\\,\\bigl((1-t)I+tA\\bigr).\n", "$$\n", "\n", "Dibujamos las elipses de confianza de $\\mu_t$ y verificamos la propiedad geodésica con la fórmula cerrada de la Sección 2." ] }, { "cell_type": "code", "execution_count": null, "id": "1c009675", "metadata": {}, "outputs": [], "source": [ "def elipse(m, S, ax, **kw):\n", " w, V = np.linalg.eigh(S)\n", " th = np.linspace(0, 2*np.pi, 200)\n", " pts = (V * np.sqrt(w)) @ np.vstack([np.cos(th), np.sin(th)]) * 2 # 2 desvíos\n", " ax.plot(m[0] + pts[0], m[1] + pts[1], **kw)\n", "\n", "I = np.eye(2)\n", "fig, ax = plt.subplots(figsize=(7, 5))\n", "for t in np.linspace(0, 1, 6):\n", " Bt = (1-t)*I + t*A\n", " mt, St = (1-t)*m0 + t*m1, Bt @ S0 @ Bt\n", " elipse(mt, St, ax, color=plt.cm.viridis(t), lw=2, label=f't={t:.1f}')\n", " ax.plot(*mt, 'o', color=plt.cm.viridis(t))\n", "ax.axis('equal'); ax.legend(); ax.set_title(r'Geodésica de $W_2$ entre dos gaussianas'); plt.show()\n", "\n", "print(\"Verificación de W_2(mu_s, mu_t) = |t-s| W_2(mu_0, mu_1):\")\n", "for s, t in [(0, 0.5), (0.3, 0.8), (0.5, 1)]:\n", " Bs, Bt = (1-s)*I + s*A, (1-t)*I + t*A\n", " d = W2_gauss((1-s)*m0 + s*m1, Bs @ S0 @ Bs, (1-t)*m0 + t*m1, Bt @ S0 @ Bt)\n", " print(f\" (s,t)=({s},{t}): {d:.6f} vs {abs(t-s)*w_formula:.6f}\")" ] }, { "cell_type": "markdown", "id": "9877a1cb", "metadata": {}, "source": [ "## 6. Dualidad de Kantorovich–Rubinstein\n", "\n", "Para $p=1$ el teorema de dualidad toma una forma particularmente limpia:\n", "\n", "$$\n", "W_1(\\mu,\\nu)=\\sup\\Bigl\\{\\int\\varphi\\,d(\\mu-\\nu):\\ \\varphi\\ \\text{1-Lipschitz}\\Bigr\\}.\n", "$$\n", "\n", "Para medidas discretas sobre un conjunto finito de puntos $\\{z_k\\}$ (la unión de los dos soportes), el supremo se puede calcular como un programa lineal en las variables $\\varphi_k=\\varphi(z_k)$: maximizar $\\sum_k\\varphi_k(\\mu_k-\\nu_k)$ sujeto a $|\\varphi_k-\\varphi_l|\\le|z_k-z_l|$ para todo par $k,l$. (Toda función 1-Lipschitz sobre un subconjunto de $\\mathbb R^d$ se extiende a una 1-Lipschitz en $\\mathbb R^d$ —extensión de McShane—, así que restringirse a los valores en los puntos no pierde nada.)\n", "\n", "Resolvemos ese programa lineal con `scipy` y comparamos con el valor primal que calcula POT: el teorema afirma que coinciden." ] }, { "cell_type": "code", "execution_count": null, "id": "43c044df", "metadata": {}, "outputs": [], "source": [ "# medidas discretas en el plano con soportes distintos\n", "n, m = 12, 15\n", "X, Y = rng.random((n, 2))*3, rng.random((m, 2))*3 + np.array([1.0, 0.5])\n", "a, b = rng.dirichlet(np.ones(n)), rng.dirichlet(np.ones(m))\n", "\n", "# primal\n", "M = ot.dist(X, Y, metric='euclidean')\n", "W1_primal = ot.emd2(a, b, M)\n", "\n", "# dual K-R: variables phi_k en Z = X ∪ Y ; maximizar sum phi_k (mu_k - nu_k)\n", "Z = np.vstack([X, Y]); N = len(Z)\n", "diff = np.concatenate([a, -b]) # mu - nu como vector sobre Z\n", "Dz = ot.dist(Z, Z, metric='euclidean')\n", "rows, rhs = [], []\n", "for k in range(N):\n", " for l in range(N):\n", " if k != l:\n", " r = np.zeros(N); r[k], r[l] = 1, -1 # phi_k - phi_l <= |z_k - z_l|\n", " rows.append(r); rhs.append(Dz[k, l])\n", "res = linprog(-diff, A_ub=np.array(rows), b_ub=np.array(rhs), bounds=[(None, None)]*N, method='highs')\n", "phi = res.x\n", "W1_dual = diff @ phi\n", "\n", "print(\"W_1 primal (POT) :\", W1_primal)\n", "print(\"W_1 dual Kantorovich-Rubinstein:\", W1_dual)\n", "print(\"phi es 1-Lipschitz sobre Z:\", np.all(np.abs(phi[:, None] - phi[None, :]) <= Dz + 1e-9))" ] }, { "cell_type": "code", "execution_count": null, "id": "645d80d1", "metadata": {}, "outputs": [], "source": [ "# dibujamos el potencial 1-Lipschitz óptimo y el plan óptimo\n", "P = ot.emd(a, b, M)\n", "fig, ax = plt.subplots(figsize=(7, 5.5))\n", "sc = ax.scatter(*Z.T, c=phi, cmap='coolwarm', s=400*np.concatenate([a, b]), edgecolor='k', zorder=3)\n", "for i in range(n):\n", " for j in range(m):\n", " if P[i, j] > 1e-10:\n", " ax.plot([X[i, 0], Y[j, 0]], [X[i, 1], Y[j, 1]], color='gray', lw=8*P[i, j]/P.max(), alpha=0.6)\n", "ax.scatter(*X.T, marker='o', facecolor='none', edgecolor='C0', s=200, lw=2, label=r'sop $\\mu$')\n", "ax.scatter(*Y.T, marker='s', facecolor='none', edgecolor='C1', s=200, lw=2, label=r'sop $\\nu$')\n", "plt.colorbar(sc, label=r'$\\varphi$ (1-Lipschitz óptimo)')\n", "ax.axis('equal'); ax.legend(); ax.set_title(r'Plan óptimo (grises) y potencial de Kantorovich–Rubinstein (color)')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "c4dcb6b5", "metadata": {}, "source": [ "Observar que $\\varphi$ **decrece a lo largo de cada segmento activo del plan**: si $\\pi_{ij}>0$ entonces $\\varphi(x_i)-\\varphi(y_j)=|x_i-y_j|$, es decir, sobre el soporte del plan la desigualdad de Lipschitz se satura. Es la condición de holgura complementaria del programa lineal, y es la versión $p=1$ del criterio de optimalidad $\\varphi(x)+\\psi(y)=c(x,y)$ en el soporte, con $\\psi=\\varphi^c=-\\varphi$." ] }, { "cell_type": "code", "execution_count": null, "id": "9d42bee4", "metadata": {}, "outputs": [], "source": [ "sat = [(i, j, phi[i] - phi[n + j], M[i, j]) for i in range(n) for j in range(m) if P[i, j] > 1e-10]\n", "print(f\"{'i':>3} {'j':>3} {'phi(x_i)-phi(y_j)':>18} {'|x_i-y_j|':>10}\")\n", "for i, j, d, c in sat[:10]:\n", " print(f\"{i:>3} {j:>3} {d:>18.6f} {c:>10.6f}\")\n", "print(\"...\\nSaturación en todo el soporte del plan:\", all(abs(d - c) < 1e-7 for *_, d, c in sat))" ] }, { "cell_type": "markdown", "id": "d33e31f5", "metadata": {}, "source": [ "## Ejercicios computacionales\n", "\n", "Los enunciados siguientes figuran también en la sección de ejercicios del capítulo correspondiente de las notas.\n", "\n", "1. **Escala.** Probar numéricamente (y luego a mano) que $W_p(\\lambda_\\#\\mu,\\lambda_\\#\\nu)=\\lambda\\,W_p(\\mu,\\nu)$ para la dilatación $x\\mapsto\\lambda x$, y que $W_p$ es invariante por traslaciones simultáneas.\n", "\n", "2. **Traslaciones.** Para $\\nu=\\tau_{v\\,\\#}\\mu$ (trasladar $\\mu$ por el vector $v$) se tiene $W_p(\\mu,\\nu)\\le|v|$. Probar que para $p=2$ hay siempre igualdad, y verificarlo numéricamente con medidas discretas. (Sugerencia: probar primero la descomposición $W_2^2(\\mu,\\nu)=|m_\\mu-m_\\nu|^2+W_2^2(\\tilde\\mu,\\tilde\\nu)$, donde $m$ denota la media y $\\tilde\\mu,\\tilde\\nu$ las medidas centradas.)\n", "\n", "3. **Geodésicas lineales.** Calcular la longitud de la curva $t\\mapsto(1-t)\\mu+t\\nu$ en $(\\mathcal P_2,W_2)$ para $\\mu=\\delta_0$, $\\nu=\\delta_1$ en $\\mathbb R$: mostrar que $W_2\\bigl((1-s)\\mu+s\\nu,(1-t)\\mu+t\\nu\\bigr)=\\sqrt{|t-s|}$ y deducir que la longitud es infinita. Comparar con la Sección 5.1.\n", "\n", "4. **Un tercer exponente.** Extender la Sección 4 para exhibir, para cada par $1\\le q