{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"tags": []
},
"source": [
"[](https://colab.research.google.com/github/bgoujaud/PEPit/blob/master/ressources/demo/PEPit_demo_extracting_a_proof.ipynb)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# PEPit : numerical examples of worst-case analyses"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This notebook provides:\n",
"- two simple examples illustrating how to extract numerical values for the dual variables to analyse **gradient descent** in the smooth and strongly convex setting.\n",
"- We show how to search for simplified versions of the proofs, if possible, by explicitly removing certain constraints from the performance estimation problem and by observing they play no role in certain performance guarantees.\n",
"- Finally, we show how to leverage those numerical findings to find a convergence proof with SymPy. More details about such strategies can be found in the repository [learning performance estimation](https://github.com/PerformanceEstimation/Learning-Performance-Estimation), from which those examples were taken.\n",
"\n",
"The examples below are taken from [this blog post](https://francisbach.com/computer-aided-analyses/) and in the [habilitation thesis](https://hal.science/tel-04794552v2/file/HDR_ATaylor_V_HAL2.pdf) that contain a more mathematical introduction to performance estimation problems."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This notebook is organized as follows:\n",
"* [1. Gradient descent, distance to a solution](#example1) : numerical study of **gradient descent** for smooth strongly convex minimization, convergence in terms of $\\frac{\\|x_{k+1}-x_\\star\\|^2}{\\|x_k-x_\\star\\|^2}$ (basic example).\n",
"* [2. Gradient descent, function value accuracy](#example2) : numerical study of **gradient descent** for smooth strongly convex minimization, convergence in terms of $\\frac{f(x_{k+1})-f_\\star}{f(x_{k})-f_\\star}$ (more complicated example).\n",
"* [3. Function values: numerical proof simplification](#example3) : a base (greedy) approach to search for which inequalities are needed for a proof of convergence in terms of $\\frac{f(x_{k+1})-f_\\star}{f(x_{k})-f_\\star}$.\n",
"* [4. Using SymPy for the distance analysis](#example4) : leverage the SymPy package to get a closed form solution (and proof) for the analysis of gradient descent in terms of $\\frac{\\|x_{k+1}-x_\\star\\|^2}{\\|x_k-x_\\star\\|^2}$ (basic example).\n",
"* [5. Using SymPy for the function value accuracy analysis](#example5) : leverage the SymPy package to get a closed form solution (and proof) for the analysis of gradient descent in terms of $\\frac{f(x_{k+1})-f_\\star}{f(x_{k})-f_\\star}$ (more technical and illustrates why step [3.](#example3) is useful/needed.\n",
"* [6. Using symbolic regression (PySR package) for guessing closed-forms](#example6)"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# import PEPit and the required function class\n",
"from PEPit import PEP\n",
"from PEPit.functions import SmoothStronglyConvexFunction\n",
"# import numpy and matplotlib\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# set up the smoothness and strong convexity constant for the whole notebook\n",
"L, mu = 1, .1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## Table of Contents\n",
"- [1. Gradient descent, distance to a solution](#example1)\n",
"- [2. Gradient descent, function value accuracy](#example2)\n",
"- [3. Function values: numerical proof simplification](#example3)\n",
"- [4. Using SymPy for the distance analysis](#example4)\n",
"- [5. Using SymPy for the function value accuracy analysis](#example5)\n",
"- [6. Using symbolic regression (PySR)](#example6)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. Gradient descent, distance to a solution "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We start by a direct attempt: fix class parameters $L,\\mu$ (chosen above) and experiment with different step size values $\\gamma$. Verify that the obtained convergence rate is indeed smaller than one only when $\\gamma\\in\\left(0,\\frac{2}{L}\\right)$. Further verify that it matches the very well-known $\\max\\{(1-\\gamma L)^2,(1-\\gamma\\mu)^2\\}$."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"def wc_gradient_descent_distance(mu,L,gamma, verbose):\n",
" # Instantiate PEP\n",
" problem = PEP()\n",
"\n",
" # Declare a smooth convex function\n",
" f = problem.declare_function(SmoothStronglyConvexFunction, L=L, mu=mu)\n",
" \n",
" # Start by defining its unique optimal point xs = x_* and corresponding function value fs = f_*\n",
" xs = f.stationary_point()\n",
" fs = f(xs)\n",
" \n",
" # Then define the starting point x0 of the algorithm\n",
" x0 = problem.set_initial_point()\n",
" \n",
" # Set the initial constraint that is the distance between x0 and x^*\n",
" problem.set_initial_condition((x0 - xs) ** 2 <= 1)\n",
" \n",
" # Run n steps of the GD method\n",
" x1 = x0 - gamma * f.gradient(x0)\n",
" \n",
" # Set the performance metric to the function values accuracy\n",
" problem.set_performance_metric( (x1 - xs) ** 2 )\n",
" \n",
" # Solve the PEP\n",
" pepit_verbose = max(verbose, 0)\n",
" pepit_tau = problem.solve(verbose=pepit_verbose)\n",
" \n",
" return pepit_tau, problem._list_of_prepared_constraints # outputs optimal value and list of constraints (with their associated dual variables)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following pieces of code compute and plot the numerical worst-case guarantees as a function of the tested step size value $\\gamma$, for fixed $\\mu,L$."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"verbose = 0\n",
"\n",
"gamma_min, gamma_max = -1, 3\n",
"nb_gammas = 50\n",
"gamma_list = np.linspace(gamma_min,gamma_max,nb_gammas)\n",
"\n",
"pepit_worst_case_value = list()\n",
"known_worst_case_value = list()\n",
"\n",
"for gamma in gamma_list:\n",
" pepit_tau, _ = wc_gradient_descent_distance(mu,L,gamma, verbose)\n",
" pepit_worst_case_value.append(pepit_tau)\n",
" known_worst_case_value.append(max((1-gamma*L)**2,(1-gamma*mu)**2))"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"tags": []
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEICAYAAACeSMncAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA5nUlEQVR4nO3dd3gU1dvG8e+TkBB6r6GFLiA1otKlqBSpERCkKIIIImDDDtafL1YQRXqRFnoRVEBBQBEJCEgLJKEFgoTQW+p5/9gFYkxIIbuzyT6f65qLmd0zM/cuhCfTzhFjDEoppVRGeFgdQCmlVNalRUQppVSGaRFRSimVYVpElFJKZZgWEaWUUhmWw+oAzlS0aFFToUIFq2MopVSWsmPHjrPGmGLJvedWRaRChQoEBQVZHUMppbIUETmW0nt6OksppVSGaRFRSimVYVpElFJKZZgWEaWUUhmmRUQppVSGuXQRERFPEflLRL5P5r2cIhIoIiEisk1EKlgQUSml3JpLFxFgOHAghfcGAOeNMZWBL4D/c1iKCxfgu+9g9WqH7UIppRxizx4YPx7Cwx2yeZctIiJSBmgPTE2hSSdgln1+MdBKRCTTg6xeDcWLQ9++JHw8NtM3r5RSDjVrFgwfDmXLwocfZvrmXbaIAF8CrwIJKbzvC5wAMMbEAReBIkkbicggEQkSkaDIyMj0p/D3Z2NcE+7nD17e0hlOn07/NpRSygrGwJIlt5cbNMj0XbhkERGRDsAZY8yOu92WMWayMcbfGONfrFiyT+3fWYkS5KxTnT+5n8V0wyxddreRlFLKOXbsgGP2h80LFICWLTN9Fy5ZRIDGQEcROQosAFqKyJwkbU4CZQFEJAdQAIhyRJj7n6qBL+GcoBzbZ/ztiF0opVSmiwlcxjVy2RY6dQJv70zfh0sWEWPM68aYMsaYCkBP4BdjzJNJmq0E+tnnA+xtHDLWr0dAV7qyFIDFOypCRk6LKaWUMxnDyjmXKEYkoxkDAQEO2Y1LFpGUiMh7ItLRvjgNKCIiIcCLwGsO23Hp0gTUPAjAEtMFs3yFw3allFKZYs8eFp9uzDXykN87Gtq0cchuXL6IGGM2GmM62OffMcastM/fMMY8boypbIxpaIwJc2SOxk9VpQSnCaMSu2bucuSulFLqrl1fsILVtAeg2yNXwMfHIftx+SLiKjwDutAF20X1xVt94fx5ixMppVTK1s45wxXy4c92KvRr7rD9aBFJq/LlebbmFubSi1Hmf7BypdWJlFIqefv3szj8fgC6ea2Ctm0dtistIulQt28dejGf/Fz+973XSinlQqIDl7OKxwDo1vI85M7tsH1pEUmPbt1uz//0E1y6ZF0WpZRKwdY5oVykIHXYRZX+jR26Ly0i6VGpEmdrtWAAU3ko5kftS0sp5XoOH6ZF2HQOU5lvvEZA+/YO3Z0WkXTK360NS+jGRh4ieMbvVsdRSql/s59qr0wojdoWgHz5HLo7LSLp5N2zK52wPSeyZENhuHrV4kRKKXVbzMLltxcc9IBhYlpE0qt6dbqV3Q7AkriO8MMPFgdSSim7o0cZ9tdT3MseNni2hscec/gutYhkwMN9SpCXy+ykAWEzN1kdRymlAIhftJRldGEv91L4gapQsKDD96lFJAN8enamA7bBFpesyw/Xr1ucSCmlYPPMUCIpTmUOU/upzO/2PTlaRDKiVi0CSv4GwJKYDrB2rcWBlFJuLzycxfvvAaCbLEM6d3LKbrWIZIQIbZ8swpt8wESeg8WLrU6klHJzCYuXspSuAATcdwyK/GeMPofI4ZS9ZEO5e3bkg0/9bQsrwyA6GnLmtDaUUsptbZ0ZTASlKc9RGvS/12n71SORjKpfHypUsM1fuqSntJRS1jl1iqW7KwHQjaVIl85O27UWkYwSge7dmcQgGvEbx6evtzqRUspdLVrEh7zBUrow6P7dULKk03atReRu9OjBz7RiK41YuCav3qWllLJGYCA+RNOF5VR7qpFTd61F5G7Uq0ePkrbnRAJjOuuDh0op5zt+nIStf9jmPT3/3VGsE2gRuRsitOtXjLxcJoj7CJ220epESik3Ezd/EVU4zBPM40qLDlC0qFP3r0XkLuV6shsdsQ1QtXBdIe1LSynlVBunhxFGJXbQgDy9nPNsSGJaRO5WrVr0LGvrzXdBbFdYtcriQEoptxESwoJD9QDo6bHIqXdl3aRFJBM8/JQvBbjAHupwcOoWq+MopdxEzLzFtx4w7NHsFBQq5PQMLllERMRHRP4Ukd0isk9E3k2mTX8RiRSRXfbpGSuyAuTsHcB7vMNs+lBm83wd8VAp5RTrZ5zgPIWpyV5qPuPcu7JuctUn1qOBlsaYKyLiBWwRkR+MMX8kaRdojHnegnz/VrUqL9TdDLt2QQywYgX06WN1KqVUdnbwIIFHGwLQI8cS6PiiJTFc8kjE2FyxL3rZJ2NhpNT16HF7PjDQuhxKKbcQM3cRy+kMQI+HIh0+gmFKXLKIAIiIp4jsAs4A64wx25Jp1k1E9ojIYhEpm8J2BolIkIgERUZGOi5wjx7s5x6e5Vs++aEmnDvnuH0ppdybMXgvnsefNORrhlB1YHPLoogxrv0LvogUBJYBw4wxexO9XgS4YoyJFpFngR7GmJZ32pa/v78JCgpyWNZf7hlCq4PfUIkQDk/5FXlmgMP2pZRyY3v2QJ06tvk8eeDMGcid22G7E5Edxhj/5N5z2SORm4wxF4ANwKNJXo8yxkTbF6cCzhmB5Q6aD6hCCU4TSmV2Tt1pdRylVDZlFiQ6Zf7YYw4tIKlxySIiIsXsRyCISC6gDXAwSZtSiRY7AgecFjAFnj0CCMA2tkjgn37gyNNnSin3ZAyLp1+iLn8xh97/vh5rAZcsIkApYIOI7AG2Y7sm8r2IvCciHe1tXrDf/rsbeAHob1HW28qWpUet/QAEmscxi5dYHEgple3s3MmCf1qwm7pE5iwLjz6a+joO5JK3+Bpj9gD1knn9nUTzrwOvOzNXWjQeWAPf4eEcpzx/TPmbB5+zOpFSKju5/N1y1vAGAI93uA4+PpbmcdUjkSzLo3sA3VkEQOBfVeHUKYsTKaWyDWNYOecSN8hFEzZT5umHrU6kRSTTlSxJH/8DvMwnPMV0WLjQ6kRKqezi99+ZH9UGgB65VkHr1hYH0iLiEPUG3ccnvEod9sDcuVbHUUplE2enLucnHsGTOLp3iwdvb6sjaRFxiICA23+5QUFw6JC1eZRSWV9MDIsWQRxetGEdxQc6v9v35GgRcYRChTDt2vM1Q2jDWq7P1G5QlFJ36aefeOrqVywigNeKTYcmTaxOBGgRcRh5sjcz6c962vD9tH/AxXsGUEq5uLlz8SGaAJbQ/OlK4OEa/327RorsqH17evssBWDumdawLbmuv5RSKg0uX8asWHl7uXdv67IkoUXEUXx86Nn5Bh7Es4Z2nJu2zOpESqksyixdxv03NtKH2Zy9pynce6/VkW7RIuJAJQc+RmvWE4s3ixfEQWys1ZGUUlnQXxP/YDsN+ZFHKdCnY+orOJEWEUdq3pzehdYAMOdKJ1i3zuJASqks5/Rp5m6rDEBPAvF60tq+spLSIuJInp506ZOPXFxjM804NvknqxMppbKY+HmBzKcnAL3r7YeyyQ6dZBktIg6W76kA3uE9ZtKPImvnw5Urqa+klFJ2GycFE0FpKhHC/YP/06Wg5bSIOFqdOrxWYxX9mE3e65G28deVUiotgoOZe8g2FlQvj0Dk8QCLA/2XFhFHE/n37Xhz5liXRSmVpUTPWsASugHQu2UEFCpkcaL/0iLiDL16cYgqPMc3fPRTA9tQlkopdSfGkDNwNn/wABMYSrXBD1mdKFlaRJyhQgUia7fmW55jonmWhPnaDYpSKhXbtkFYGPdwkKEF5kL79lYnSpYWESdpNLg2FThCOGXZNMnykXyVUi4u4btEPYB362b54FMp0SLiJNL9cXrJAgDmHqgHISEWJ1JKuazYWKbO9qIeO1lKF5fq5iQpLSLOUqQIvZuHA7CYAG7MXGBxIKWUy1q7ljlXOrOLelwuWA6aN7c6UYq0iDhRjeeaU4+dXKAQq6ZEQEKC1ZGUUi4o7Osf2EwzcnGNLn3ygqen1ZFSpEXEmTp2pH8u20X1mWfawZYtFgdSSrmc8+eZ/VMJALqxhPyDeloc6M5csoiIiI+I/Ckiu0Vkn4i8m0ybnCISKCIhIrJNRCpYEDV9fHzo1T2e4XzJB7wFs2ZZnUgp5WIS5i1gVsKTAPSrshVq1bI40Z25ZBEBooGWxpg6QF3gURF5IEmbAcB5Y0xl4Avg/5wbMWOKDg7gS0ZSj12wcCFcvWp1JKWUC9n89R6O4kdZjvPQ8zWtjpMqlywixuZmJ1Ne9inp0ICdgJu/yi8GWomIOClixt1/P1SrZpu/cgWW6TgjSim7gweZd6AuAH095uLZ27VPZYGLFhEAEfEUkV3AGWCdMSbp0IC+wAkAY0wccBEoksx2BolIkIgERUZGOjh1GohAv37MpyeN2ULQ+N+tTqSUchWzZvElI5hPTwa0OQ5F/vNfmstx2SJijIk3xtQFygANRSRDJwaNMZONMf7GGP9ixYplasYM69OHbTzA7zRm1vYacPy41YmUUlaLj4fZs8nFDXoSiN+QtlYnShOXLSI3GWMuABuAR5O8dRIoCyAiOYACQJRTw2VUmTL0u/8gAPN4gugZ8ywOpJSy3M8/E33qrG2+WDFoq0Ukw0SkmIgUtM/nAtoAB5M0Wwn0s88HAL8YY5JeN3FZdYc1pTa7OUcRVk8Kh6wTXSnlAMe+/p7inGEoE2xPqHt5WR0pTVyyiAClgA0isgfYju2ayPci8p6I3BxgeBpQRERCgBeB1yzKmiHSpTP9ctqeWp8V0Qa2brU4kVLKMhcv8t2awlyiAOcoDP37W50ozXJYHSA5xpg9wH+G8DLGvJNo/gbwuDNzZarcuendLZpX58WxhnacmfgaxRs1sjqVUsoCJnAhs+Js/WP199sEdZ6wOFHaueqRiFso8VxX2rGGOLyYt8gLrl+3OpJSygK/T9hJCFUozUlaD61mdZx00SJipcaNean0AqbxNE9HfwPLl1udSCnlbIcPM/Pv+gD0kbl49ullcaD00SJiJRGaD76Hp5lBfi5rNyhKuaFr0+azkO4A9Gt+FIoXtzZQOmkRsVrfvrdmzdp1cPKkhWGUUk6VkMCG6Ue4RAEaso17hrW2OlG6aRGxWvny3Gj2MMMYT03zN7Gz9JkRpdzGxo20j5zJPmrwZb53XHYI3DvRIuICcj7dm19oyQFqsPrrI/rMiFLuYvp0AGpwgAf7VYWcOS0OlH5aRFyABHTj6Zy28ZSnnWqr44wo5Q7Onydq0S+3l596yrosd0GLiCvIk4e+PWLwIoY1tCN83BKrEymlHCz+u3nUj9nKffxJxL0PQ/36VkfKEC0iLqLYsJ50ZjkJeDJzRSG4cMHqSEopRzGG9V/u5TjliaIIJZ7tbHWiDNMi4ioaNOCZihsAmBbXl4Q5eoFdqWwrKIgpR1oBMCDHbDyezFrPhiSmRcRViNB6RC3Kc5Sj+PHb+B1WJ1JKOciZrwJZQSc8iKd/5wtQoIDVkTJMi4gL8XiyF197jeRP7qPJ4emwc6fVkZRSme3KFWYH5iQOL9qzGt/hAVYnuitaRFxJoUK075GX+whCAKZMsTqRUiqTmcCFTI3pA8Azvj9C48YWJ7o7WkRczcCBt2avzl0OV69al0UplemOfL2GU5SmFKdo90Jl25DZWZgWEVfTtCnhFZrQgg00vLwes2ix1YmUUpll3z4q/rWECErxfY7O5Oj/pNWJ7poWEVcjQolBnThIdfZTk61f/GF1IqVUZpk6FYA8XKN+5/JZrrPF5GgRcUFeT/ehv8wGYOqe++DAAYsTKaXuWnQ0h2b8xjVy2ZYTnbrOyjKliIhIfGZsR9mVKMGANscBCKQHl76ZY3EgpdRdW7aMnhe/pTSnCCrZAVpnvR57k5NZRyJZ+8qQC6oyoj0t2MA18jB/xg2IibE6klLqLuz8fCN/UR9P4qn1zAPgkT1OBN3xU4hIbxF5VUTyi0ibOzQ16WyvUvPwwzxTeBkAU6/2hJUrLQ6klMqwsDCmbK8DQB/m4DOobyorZB2plcKKwARgJPBoGrZXKZ3tkyUiZUVkg4jsF5F9IjI8mTYtROSiiOyyT+9kdH8uydOTroOLU5Dz7KMm4V8tszqRUiqDrk6czTxsXZs80zQYypa1OFHmSa2I7DDGXAPeA86kYXtB6WyfkjjgJWNMDeABYKiI1Eim3WZjTF379N5d7M8l5RrUh5V04hSlKbN5PoSGWh1JKZVeMTEETr7IJQpwP39Qa2T2OklzxyJijFlj/9MYY/4vtY2lt/0dthNhjNlpn78MHAB8M7q9LKt8eZq2y0dBLtoGqpo0yepESqn0Wr6ciZdsRyHPFZgPHTpYHChz5UitgYiMxn7NI6mUfvvPyDp32H8FoB6wLZm3HxSR3cAp4GVjzL5k1h8EDAIoV65cenbtGoYMgTVriMabE1M2UPm9G+DjY3UqpVQaXRg/GxhNIc7RfWgx8PKyOlKmEpPKUKwiUj6l94wxx+xt4o0xnulZJ03hRPICvwIfGmOWJnkvP5BgjLkiIu2AccaYKnfanr+/vwkKCkrr7l1DfDx7y7alVcR3FOcMe2btQvr2sTqVUiot9u+HmjUBOOVRhtLH/wDfrHdSRUR2GGP8k3sv1XvMjDHHUpoyc51kQnsBS4C5SQuIfR+XjDFX7PNrAC8RKZrW7WcZnp5UHdIawbCXe9nyf79ZnUgplVbffntrtnTnhlmygKQmQzcqi0hhkbT3GpaB9gJMAw4YYz5PoU3Jm9sUkYbYPktUWveRlXgP6s8zHjMAmLi/Gfz1l8WJlFKpunKFDdPCOEB12/Jzz1mbx0HSXUREpCSwFm4+u5+57e0aA32Alolu4W0nIoNFZLC9TQCw135NZDzQ06R2bi6rKl6cQY9F4EE8iwngn8/0CXalXF3CnHk8e+1zanCAX8v0hpYtrY7kEKleWE/KGHNaRFrZb+XN9Pb2dbaQylPwxpgJ2J5JcQvlXulBhxXfs5JOTF+Yl9e/vpilR0NTKlszhp8/3clhBlGGEzQecV+2eUI9qTR9KhHxTLxsjLmYnnXS0l6lolEjhlT4AYBvY58mfuZ3FgdSSqXojz+YGPowAM/mmE6Op7PPE+pJpbU0ThaR3AAi0syB66iUiNDmlbpUIoS8XOHkhGW2Z0eUUi4n/NMFrKATOYjlmYALUKiQ1ZEcJq1F5B1gmoh8B9znwHXUHXj06c2m3G3ZSy3KhfwCv/5qdSSlVFJnzzJleTES8KQrSyn5Um+rEzlUWovI+0AwtgcIFzpwHXUn+fJRuv/Dty8WffONlWmUUsmInTqLKQlPAzCk2i/gn+zjFdlGWovIq8aYMcBzwGgHrqNSY79N8AgV+GHJNYiIsDiQUuqWhAQufzuXtvyAP9tpNupBqxM5XJqKiDHmrP3Pq8CzjlpHpUGtWhxs0JtKhNInYSY3Js6wOpFS6qa1ayl87C+m8QzbCj6K9OxhdSKHS/M9ZzefBjfGpHkUw4yso1JX7aUO1GUXURRl8YTTEBdndSSlFPzrFLPH0/0hV3oej8ua0nPj8vQMbD8j66hUSLeuPJdvLgBfne8Ny5dbG0gpBWFhjF/lRyDdiSUHDB6c+jrZQHqKSEaGwNVhcx3B25veQwpQiHP8yf388f46qxMp5fYufDKFN/iQngQS3OhpqHLH/mCzjfQUkYw8lKAPMjhI7heeYZDHNADG7WkBO3ZYG0gpd3b5MtNnCFfJSyvWU+vtLlYncho9EsmqSpdmaMcTeBLHIh4n/KPZVidSym3FT5/FV9EDARjuuwQeftjiRM6TniLyega2n5F1VBqVfbMvPVnA4ywiZuWPeruvUlZISGDlx/s5ih+VCKH967WzbT9ZyUnzJzXG7E3vxjOyjkoHf3++e3Ai8+lFxbhD/xq7QCnlJGvW8OVp2628L+Sagkf/7NtPVnLcp1xmUzJi+O2FiRPhxg3rwijlhna9v4pNNCcfl+j/jBfkyWN1JKdKd1fwN4lIa2PM+lTaPAmUti+eNMbMzej+VAq6diWhTDnWhN/LysiOTJq/AHmqv9WplHIPe/dS6c95jMebG5KL/C8NtTqR06XrSERE5ovIqyIyChiWhlV8jTFjjTFjgTIZSqjuLEcO4gY/zyAmM4VB/PzBVu3dVylnGTeOfFxhGBN4pWsYlC9vdSKnS+94It/Yi8L/AW+mYZ3f7UXnZeD3u4uqUuL93ACG5JgCwLiwDrBpk8WJlHIDZ89ivks0yuiIEZZFsVK6xhMxxmy+OTZIGi6aTwZ22I9CthtjNt9NUHUHhQvzbK/L5OQG3/MYhz/UTpOVcrToidO5N3o7b/IB0XXvh8aNrY5kCWeNJ5K9+0J2AcVGPU1v7F2hrKsOR45YnEipbCw2loVfnGQftfieDniPHArino/F6Xgi2UWNGgx/cDsAM+jPxc+mWhxIqezLLF7Cl+dtt/IOzz8T6dHd4kTW0fFEspHab3XkIX7hCvmYPjUBLl+2OpJS2dJvH2xgJw0oSiS9hhWBnDmtjmSZVG/xFZHRgJHbh2oRIvIOgDHmvcxaJ8n6ZYHZQAlsRzKTjTHjkrQRYBzQDrgG9DfG7Ext29nao4/yuu8AHj65lr7RU2B6SRg+PPX1lFJpt3UrX+y3dWsy2GMKPsMGWhzIWml5TmRmGtokPRmYlnXuJA54yRizU0TyATtEZJ0xZn+iNm2BKvbpfmCi/U/35eFBmzfuo81Q+73qn38OQ4aAl5e1uZTKRg6/8x3LmIA30QwJOAMlSlgdyVKpns4yxhxLaUrUxiO966Syz4ibRxXGmMvAAcA3SbNOwGxj8wdQUERKpWX72Vr//lC0KADRx09jFi6yNo9S2cnBgyxcXxiDB32ZTal33PsoBLJAtyciUgGoB2xL8pYvcCLRcjj/LTSIyCARCRKRoMjISIfldBm5c8OwYUxiEH4c4Zd3NurDh0plls8+4w0+ZAMteP2hbVCzptWJLOfSRURE8gJLgBHGmEsZ2YYxZrIxxt8Y41+sWLHMDeiqhg7lrFcpIijN2LBusE4HrVLqrkVEwOzZCNCCX6n4bj+rE7kEly0iIuKFrYDMNcYsTabJSaBsouUy9tdUkSI81/8GebjCWh5h19tLrE6kVJZ35dNvORjjZ1t44AFo0sTaQC7CJYuI/c6racABY8znKTRbCfQVmweAi8YYHVDDrvAbgxkotpEPP/mzmY58qNTduHSJKd/Ecg8HeZv3YNQot324MCmXLCJAY6AP0FJEdtmndiIyWEQG29usAcKAEGAKMMSirK6pQgVGdgzFkzgC6cHRMTOtTqRUlhU7cSpf3LD919Ow9Eno2NHiRK4jw13BO5IxZgupDK1rjDGA+/W7nA7lxjzNEyvmM4c+fPF9FcaFhUHFilbHUipriYkh8OMjnKAc97Cf9qP93WrkwtToN5Gd1a3LKw/+BsDvPEjCZ19YHEiprMfMncfYC7ZbeV/JPxmPvk9anMi1aBHJ5mq//zi/0Yht3I/HjGngDrc5K5VZEhL4acxW/qY2pTlJr1d8wcfH6lQuRYtIdteyJY3qR+OBgevX4euvrU6kVNaxZg1jj9vGTx/hPZGcz+vDhUlpEcnuRODVVwE4Tlm2fLEdrl61OJRSWUPMx59TgaMUJZJBg4CCBa2O5HK0iLiDbt3YVqozFQmj36XxxE2ZYXUipVzf1q14/7aB6QzgmGclCowanPo6bkiLiDvIkQP/19tQgaOEUYnA94MhOtrqVEq5tg8+uDWbu3cXKFPGwjCuS4uIm/B85ilez2+7HvLhueeInzHb4kRKubAdO3hvTQOW0ZkEPOD1161O5LK0iLiLXLno83pZynGMA9RgyTu7ITbW6lRKuaTDo6byLqPpzkLC2z8L1atbHcllaRFxI97PD+L1PF8B8EHkIBJmz7E4kVIuaM8e/vfzfSTgSV9mU+5/z1mdyKVpEXEnefPy1Kji+BLO39Rm5Vt/Qlyc1amUcilHXp/Ed/TBg3heb7MD7r3X6kguTYuIm8k5fDCv5RpPV5ZQ+fRmCAy0OpJSrmP/fj5eU4c4vOjNXCp//IzViVyeFhF3kz8/Q0flZQkB1GKf7Q6U+HirUynlEk68MZEZ9EdI4I0WW6F+fasjuTwtIm5Ihr8A+fPbFg4ehCU63ohSHDrEuBUViMWb7iyk+v89ZXWiLEGLiDsqWBCGDWMLjWnJz6x9dT0kJFidSilrffQRYxjNWF7hrQd/gYYNrU6UJbhkV/DKCUaOZMtYDzbEtiT2mBdtlq9AunaxOpVS1ggLgzlzyEs8r/ApfLLF6kRZhh6JuKsiRRg6xFCIc2yhKb+OWgPGWJ1KKUtcfG8cV+LtvfM+9BA0bmxtoCxEi4gby/fGMEbmmADA+yE9YfVqixMpZYHjx/nou7JU4ChL6Apvv211oixFi4g7K16cYQNvkJ+L/EIrNr+8Qo9GlNs5+854vkl4liiKUq5OYWjRwupIWYoWETdX8O1hjPC0HY28FfwkZtlyawMp5UyhoXw8uzRXyEdb1nDf2MdtwyeoNNMi4u5KleLF565TiHPspD5HRn2rz40ot3Hy1XFMMEMA+KDeUmjTxuJEWY8WEUWB0SNY4vMkYVSkYshamDfP6khKOd6+fXywtAbR+BDAIuqP769HIRngkkVERKaLyBkR2ZvC+y1E5KKI7LJP7zg7Y7ZStCgPvXofxThrWx4zBmJiLI2klKOFjpzAVAbgQTzvNV0PTZpYHSlLcskiAswEHk2lzWZjTF379J4TMmVvL74IhQsTgxffhTUifqqOfqiyse3bubhuG/fyN32ZzT1fPmt1oizLJYuIMWYTcM7qHG6lQAF47TXa8gN9+Y55b+yF69etTqWUY7z1FvX5iyD8+arLL9pH1l1wySKSRg+KyG4R+UFEaqbUSEQGiUiQiARFRkY6M1/WM3QofQusBGD0xZHEjJtocSClHGDjRli7FgAPDyHvR29amyeLy6pFZCdQ3hhTB/gKWJ5SQ2PMZGOMvzHGv1ixYs7KlzXlzs2TH1TnHvZzhIpMe/8UXLpkdSqlMo8xbH/hO/ozg6OUh759ddTCu5Qli4gx5pIx5op9fg3gJSJFLY6VLXgOGsB7xexPsV97ketjv7I4kVKZaM0a3vy7B7Poz7ceQ2H0aKsTZXlZsoiISEkR2714ItIQ2+eIsjZVNuHtTdexD1CfHURQmq8/vQZnz1qdSqm7l5DAxuHLWMfD5Ocirz59FipUsDpVlueSRURE5gNbgWoiEi4iA0RksIgMtjcJAPaKyG5gPNDTGO2vI7N49OnNh2UnAfC/6Be59N6X1gZSKhOYRYt5M9Q2RsjLOcZR+P2RFifKHlyyK3hjzBOpvD8BmOCkOO7H05NHPn+E3o/PoTXryT1pEbw8CMqVszqZUhkTE8P3I3/mdyZRlEhGPB8HJUtanSpbcMkjEWU96daVOf7j6M8scsRcgzfesDqSUhkWO2ESr0TYjjze8vmMfG+PsDZQNqJFRCVPBD777Nbi+bmrYft2CwMplUHnzrF59HoOU4XKHOa5MSWgcGGrU2UbWkRUypo1g86dGcNofDnJX4MmalfxKuv54ANaXlnJLuoyo/SbeI8YYnWibEWLiLqzsWO5LAW4Tm5e2vUkZukyqxMplXYhITDBdvn0XvbSZFx3yJnT4lDZixYRdWdVqvDWoDMUJooNtGTVsJ+0c0aVZRwd9hmrY9tgwDbkbbduVkfKdrSIqFQV+ugVRucaC8ArES8S+9W3FidSKg22bGHUjy3owGrG8qrtGp929Z7ptIio1BUuzOB3S1OFQxyiGpPeDodz2j+mcmEJCWwdNIOF9MCH6/TqeBXuv9/qVNmSFhGVJt7Dn+OTkp8DMOb6q1x48xOLEymVMjN/AS8eeAaAlzzHUXb8KxYnyr60iKi08fam41dtaM5G/Ani0tSFtouWSrma69dZOPw3/uBBSnCaUcOuQfnyVqfKtrSIqDSTbl35/sGP+JG2lIsLg1GjrI6k1H/c+HQCr0W9DMB7eceSb8xLFifK3rSIqLQTIe+XH9xaNEuXwoYNFgZSKolTp5j04VmO4kct/ubp/1W1DbimHEaLiEqfhg2hd2/2UpOH2MC6fnMgNtbqVErZvPIKz0RP4HU+4vOyX5Jj8DNWJ8r2tIio9Bs7ltXeXfmVFjx/4lWiP9UxR5QL2LAB5s0jD9f4iDdpM7M35HDJPmazFS0iKv1Kl2bk+4WpzgEOUY3PR1+E8HCrUyl3FhvL/mc+5xL5bMs9e0LLltZmchNaRFSGeI8cyoTynwLwfuwojg3+n8WJlDuL/vQruoR9SnUOsj93A/j0U6sjuQ23P9aLjY0lPDycGzduWB3FZfj4+FCmTBm8vLxSbuTlRavZ/ejRfAGB9GTk6lYsXb8eWrd2XlClAMLD+Xz0RQ5RjWocpPLoJ8HX1+pUbkPcaUBAf39/ExQU9K/Xjhw5Qr58+ShSpAiiXSJgjCEqKorLly/j5+eXavuT3V6g+tIPuUI+1vgOpG3oBO3gTjnV8Q5DqL76U66Tm/XlB9Dq8Ldwp1+AVLqJyA5jjH9y77n96awbN25oAUlERChSpEiaj8x8v36DMTk/JhfXCD8JfPGFYwMqldj69Yxc3Yrr5KYHC2g1q68WECdz+yICaAFJIl3fR8mSvPC/UhykOgOZCu+/D8ePOy6cUjdFR/Nj/wUspRt5uMJnXX+H5s2tTuV2tIiou+Y1bDDlaheyLVy7Bi++aG0g5RZiPx3HsJO2XhPG5PwY3691CGcraBFxAZ6entStW5datWrx+OOPc+3atX+9fnP6+OOPAWjRogXVqlWjTp06NG7cmODgYAAaNWoEwNGjR5k3b57zPkCOHPD11yQgTGUAHZb0J/Th5+Dvv52XQbmXo0fx+uhdxvMCHVnB8I9KQMmSVqdyT8YYl5uA6cAZYG8K7wswHggB9gD107LdBg0amKT279//n9ecLU+ePLfme/XqZT777LP/vJ5Y8+bNzfbt240xxkyaNMk89thj/3p/w4YNpn379neVKSPfyz/dnzeFOWvAGE9izQCmmqMdhhpz4MBdZVHqXxISjGnd2hjbYM3G1K5tTGys1amyNSDIpPD/qqseicwEHr3D+22BKvZpEDAxU/Yq4rgpjZo2bUpIOnrHbdas2a32efPmBeC1115j8+bN1K1bly+ceKG7+LT/saPDGPozA4MwjQFU+f5zht7zCycDhsPhw07LorKv+Gkz2bk+yrbg4QGTJ+uT6RZyySJijNkE3GnUo07AbHuR/AMoKCKlnJPOceLi4vjhhx+49957Abh+/fq/TmcFBgb+Z51Vq1bdan/Txx9/TNOmTdm1axcjR450SnYA8ualwqqvmPFXPQ48NJTezCGOHHzDECot+T/Cqj4CxYpB06acfXIE5pNPYdUqW3GJi3NeTpV1nTrFV8MO4U8QH/AmjBihg01ZLKuWb1/gRKLlcPtrEUkbisggbEcrlCtXzinh0utmsQDbkciAAQMAyJUrF7t27Up2nd69e5MrVy4qVKjAV1+5WN9VdetS9ZdvmfPnn7zx4mDG/NaacxSmIkfgLJgtW/DbsgaZa6hGMNX5g2oec6hW8iLVqkHV+nnxqVkJqlWzTUWKWP2JlCswhtC+7/LGjS8weFCn5Bl4/0urU7m9rFpE0swYMxmYDLaHDVNp7IxI/3GnYpGSuXPn4u+f7LM/rqNhQ2psacjC334j5r2PYXMuuH6dcxTGmxjOUYQg7iOI+yABOGWbpm4YwAA+BGAbDdmRpxnVykdTrZYXvvWKI9XtxaVSJfD2tvQjKucxCwIZ+HMPrpObJ5jHY/N7Qe7cVsdye1m1iJwEyiZaLmN/TQH58uXj8uXLVse4rXFjvH9aBQkJEB5OkYMHiTo4l7O7wjm4O5qDoV4EXyxBMNUIpho12H9r1aV0ZezVUbAf2A95Fl6hKoeoxi7qymxGVVoMVavePmq5OV+qVLquRSkXFxnJlEF/soHPKUok4/ruhBa9rE6lyLpFZCXwvIgsAO4HLhpj/nMqK6tLfJoL4NFHH711m++d1K5dG09PT+rUqUP//v2de13kTjw8oFw52/TwwxQFmtgnLl+GQ4cgeAccbAPB5SA4mAcP7OTpmGm3CsxZivEX9fmL+gSbaowK+RhCQkhY8wPVCKYM4VRlFVW9j1KtzDWqVvfAr15BvO6pbCswVXWQoqwofOC7vHzlIwAmFB5Nsa9S/zlQzuGSfWeJyHygBVAU+AcYDXgBGGO+Fdsj1ROw3cF1DXjKGBOU/NZuS67vrAMHDnDPPfdkav7swGW+F/vRC8HBEBzMud0nCN59g0NhOcgZdYqeLADgBGUo96/LZLd5EscCehLAEgAOFG5MROkGVK3phW/dYki1qlCliu30WK5cTvtoKo1WraJ3x0vMozedWM6y1TmRdm2tTuVW7tR3lkseiRhjnkjlfQMMdVIcZaXERy9t2lAYeNA+cf06HH4dDh3C90AwYTtf5NCBeA4dz8mh62UJphqHqMpxyuGb6Gzn9HOd+PTcK7AXcgdepQqHqcIBqrKSOkXC6V7vsK2oVLUXlypVwM9P+2SywoULMHgwnxOHNzF82O0vpN14q1OpRFyyiCiVJrlyQe3aULs2HoCffXrEGIiKsh29HNrA9X1heIX4QmgtOHyYctHHacwWDlGVSIqzm7rspi4A/lHb6b6+IaxfTzweNOdXKvAnVWUeVYpfpEqlBKrUzkWBWmVvF5hy5cDT08IvIhsbOhROnaIEMKP4azBpf6qrKOfSIqKyHxEoWtQ2NW7Mv05QxcczLDycYYcOwaFFnN99nMN/206PHT5TgGKcudX0BGX5jSb8RhMw2E6s/gP8DsU4w3f04RHWgrc3oWWac6lMDSrXzk2+GokKTJkyWmAyKOG7ucyY50NfcuBFHHzzjd7u7YK0iCj34ukJ5cvbpjZtKAQ0tE/ExEBYGBx+BIKDKXngCJv+Gs7ho14cPl+EQ1TlEFUJoTKRFKcQ523bjIlhQlhbvgwbCZugJBFUJoQq/Exlz6PULhVJh/qnoHLl21OVKlC2rBaYlBw5wucDD/AK01hOZ1b1XwrdulmdSiVDi4hSN3l7Q/Xqtumxx/ABmtonrl2DkBA4FExC8Pec2nOWYuF5IbQE/PMPRYiiBvsIpRKnKcVpSrGFphAPD4RvpUO4rXPMODxpxc9UZAOVPY5QucRlqlRKoFKtXBSoWeZ2kSlf3n2vwcTF8Vfnd3kjejIAg0uugPF6HcRVaRFRKi1y5/7X9Zcyid+7dIm3QkJ4K2QfCcHLCd8dxeGD8YQc9+Lw5ZKUIfxW02OUZxPN2URz2wOWEfZpi+0U2Vx604b14OlJcOmHOF+6JpVreFOkVimkSmXbHWQVK4KPj1M/vjNdGzOWXntGEYs3Q+Vr2i8fCPnyWR1LpUCLiAvImzcvV65cAWDNmjWMGDGCdevWUb58eYuTqTTJnx/q14f69fEAytmnVgAXL9qOYEIehJAQSh08zs9/vUDIcW8OXy5BCJUJoTKhVCKS4hS+2WVcfDxfn3iMr068ANugABeoRCiV+ZtKLKdu4RN0v/eArahUtheXm1PBglZ9E3dv61Ze+qgwB7mHGuzjk3euaN9YLk6LiAv5+eefeeGFF/jpp5+0gGQXBQpAgwa2CcgNtLRPXL4MoaEQEkzCodVE/H2WYicLQWhpOHWKUkRQj52EUJmLFGQnDdiJbTuNz22h+69N4ddficGLBuzAjyNUYhOVcp+mcpnrVKriSfl78+NdtcLtI5jSpW23TbuiS5dY2Xk635opeBPNvHqfkuvtqVanUqnQIpLEnXrKmDQJBg2yzU+eDM8+m3Lb9D7DuWnTJgYOHMiaNWuoVKkSAP379yd//vwEBQVx+vRpxo4dS0BAAMYYXn31VX744QdEhLfeeosePXowdOhQHnnkETp27EiXLl0oVKgQ06dPZ/r06YSGhjJw4EDatm1LkyZN+P333/H19WXFihXk0gfsrJEvH9StC3Xr4oGtB9Fbrl3j9bAwXg8NxRzeQOTefwjZF03oUU9CzxagdKLnXo5Sgb3cy17svTlfAw7ZJo/V8aziMdrxAwB/ejfhWDF/KvklULFmLgreU+p2gfHzs/Zhy2HDWHimDQAf+7xLnWVj9MaDLECLiAuIjo6mc+fObNy4kerVq//rvYiICLZs2cLBgwfp2LEjAQEBLF26lF27drF7927Onj3LfffdR7NmzWjatCmbN2+mY8eOnDx5kogIW08wmzdvpmfPngAcPnyY+fPnM2XKFLp3786SJUt48sknnf6ZVSpy54ZataBWLQQobp8aAcTGwrFjENIVQkOpEHyEnbtHEBomhJ7OQ2hcOUKpRCiVOEFZSnPq1mZnxjzBxJNDbD3NbYHCRFGRMCqxB3+m83Lp+beKSnyFSnhW9rMVmIoVoUQJx/RHdvEivP8+zJ7NbL6jEyvoNq2b7eYC5fK0iCSR1iOIQYNuH5XcLS8vLxo1asS0adMYN27cv97r3LkzHh4e1KhRg3/++QeALVu28MQTT+Dp6UmJEiVo3rw527dvp2nTpnz55Zfs37+fGjVqcP78eSIiIti6dSvjx48nKioKPz+/W/1xNWjQgKNHj2bOh1DO4+V1+y4uwBuoZ59ISICICNtpstANxBw6imdYDTjiDaGh1D23i46sIJRKHMGPcxS51ZvyP5Tg5VOfwalT3Nj8JwW4SDmOU5EwKrKMijlOULHkNfz8oHqtHOSuVtZ29HLzKCZPnvR9jthY4r6dyvQ3DlPjyp80ATwwPN43N/TqmbnfmXIYLSIuwMPDg4ULF9KqVSs++ugj3njjjVvv5cyZ89Z8av2c+fr6cuHCBX788UeaNWvGuXPnWLhwIXnz5iVfvnxERUX9a3uenp5cv3498z+Qso6HB/j62qZmzUjaUf6gCxcYFBYGoQcwoav5Z99Zwg7GEHbMg/xRR2x3jAHHKUcMOQmhCiFUsb0Yh23knnBYvrkTnfgMgBV05Hca4ZfvLH6lY6hY2YPytfLhXbmcrbj4+dmeibl5y7IxmO9X8+OQlbwS/gL7eI767GA79+HRtAm42vg46o60iLiI3Llzs3r1apo2bUqJEiVuDUyVnKZNmzJp0iT69evHuXPn2LRpE5988gkADzzwAF9++SW//PILUVFRBAQEEBAQ4KyPoVxdwYK37iQToKR9unWa7PhxOHKEqmFhXD34Fkf3XSUs1BB60ocjN0pyBD/CqEhlbg/hvIZ2TOZZuAwE2yZZnYAvJ2nKZubR2nZto2xZfi/cAWKiGbM3gHW2YX7wI4xRRaYh4+bAEz1d98K/SpYWERdSuHDhW0cRxYoVS7Fdly5d2Lp1K3Xq1EFEGDt2LCVLlgRsBWbt2rVUrlyZ8uXLc+7cOZo2beqsj6CyMi+v27cJY7uTrIZ9AuD8edsT/UeCIawfHDkCoaH03LeZcqdPEpZQniP4cQQ/jlOOcMryDyVs68bHc/XoGRofvX2UUYALvJ3zE55/uxA5X/o8Wz/7kp25ZFfwjqJdwaedfi8qXeLj4eRJW2EJCyM25Bgn9l4k+vg/3PPPRoiI4CSl6cYSzlCcx+R73nkqnCIfvwJ3+IVJuYYs1xW8UiqL8fS83WV/8+Z4ARUTv3/9Or7HjvHHkSMQFQIPtr91xKOyNi0iSinHy5Xrdr9kKlvRK1ikfteTu9HvQymVVm5fRHx8fIiKitL/OO2MMURFReGjFzmVUmng9qezypQpQ3h4OJGRkVZHcRk+Pj6UKVMm9YZKKbfn9kXEy8sLPz8/q2MopVSW5Pans5RSSmWcFhGllFIZpkVEKaVUhrnVE+siEgkcy+DqRYGzmRgns7hqLnDdbJorfTRX+mTHXOWNMcl2LeBWReRuiEhQSo/9W8lVc4HrZtNc6aO50sfdcunpLKWUUhmmRUQppVSGaRFJu8lWB0iBq+YC182mudJHc6WPW+XSayJKKaUyTI9ElFJKZZgWEaWUUhmmRSQFIvK4iOwTkQQRSfG2OBF5VESCRSRERF5zQq7CIrJORA7b/yyUQrt4Edlln1Y6MM8dP7+I5BSRQPv720SkgqOypDNXfxGJTPQdPeOkXNNF5IyI7E3hfRGR8fbce0SkvovkaiEiFxN9X+84KVdZEdkgIvvtP4/Dk2nj1O8sjZms+r58RORPEdltz/ZuMm0y92fSGKNTMhNwD1AN2Aj4p9DGEwjFNoibN7AbqOHgXGOB1+zzrwH/l0K7K074jlL9/MAQ4Fv7fE8g0EVy9QcmWPDvqhlQH9ibwvvtgB8AAR4AtrlIrhbA9xZ8X6WA+vb5fMChZP4unfqdpTGTVd+XAHnt817ANuCBJG0y9WdSj0RSYIw5YIwJTqVZQyDEGBNmjIkBFgCdHBytEzDLPj8L6Ozg/d1JWj5/4ryLgVYiIi6QyxLGmE3AuTs06QTMNjZ/AAVFpJQL5LKEMSbCGLPTPn8ZOAD4Jmnm1O8sjZksYf8OrtgXvexT0runMvVnUovI3fEFTiRaDsfx/5hKGGMi7POngRIptPMRkSAR+UNEOjsoS1o+/602xpg44CJQxEF50pMLoJv99MdiESnr4ExpZcW/qbR60H6a5AcRqensndtPu9TD9tt1YpZ9Z3fIBBZ9XyLiKSK7gDPAOmNMit9XZvxMuvV4IiKyHiiZzFtvGmNWODvPTXfKlXjBGGNEJKV7tMsbY06KSEXgFxH52xgTmtlZs7BVwHxjTLSIPIvtN7OWFmdyZTux/Zu6IiLtgOVAFWftXETyAkuAEcaYS87a752kksmy78sYEw/UFZGCwDIRqWWMSfZaV2Zw6yJijGl9l5s4CST+DbaM/bW7cqdcIvKPiJQyxkTYD9nPpLCNk/Y/w0RkI7bfljK7iKTl899sEy4iOYACQFQm50h3LmNM4gxTsV1rcgUO+Td1txL/J2mMWSMi34hIUWOMwzsaFBEvbP9ZzzXGLE2midO/s9QyWfl9JdrvBRHZADwKJC4imfozqaez7s52oIqI+ImIN7aLVA67E8puJdDPPt8P+M8Rk4gUEpGc9vmiQGNgvwOypOXzJ84bAPxi7Ff0HCjVXEnOmXfEdl7bFawE+trvOHoAuJjo9KVlRKTkzfPmItIQ2/8djv5lAPs+pwEHjDGfp9DMqd9ZWjJZ+H0Vsx+BICK5gDbAwSTNMvdn0tl3D2SVCeiC7dxqNPAP8JP99dLAmkTt2mG7OyMU22kwR+cqAvwMHAbWA4Xtr/sDU+3zjYC/sd2V9DcwwIF5/vP5gfeAjvZ5H2AREAL8CVR00t9farn+B+yzf0cbgOpOyjUfiABi7f++BgCDgcH29wX42p77b1K4M9CCXM8n+r7+ABo5KVcTbBeG9wC77FM7K7+zNGay6vuqDfxlz7YXeMf+usN+JrXbE6WUUhmmp7OUUkplmBYRpZRSGaZFRCmlVIZpEVFKKZVhWkSUUkplmBYRpZRSGaZFRCmlVIZpEVHKIiJSRUSOikhl+7KXfewJV+kMUqlUaRFRyiLGmMPAZOAR+0vPAyuNMSdSXksp1+LWHTAq5QL2Aq1FpDC2rkbutziPUumiRyJKWesQthE0xwCfGmOuWhtHqfTRvrOUspC9S/FT2DoPbGSMSbA4klLpokciSlnIGBMLXAJe0wKisiItIkpZzwv41eoQSmWEFhGlLGQfo/uY0fPKKovSayJKKaUyTI9ElFJKZZgWEaWUUhmmRUQppVSGaRFRSimVYVpElFJKZZgWEaWUUhmmRUQppVSG/T8ldwAh1ssqwgAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(gamma_list, pepit_worst_case_value, color='red', linestyle='-', linewidth=3, label='PEPit')\n",
"\n",
"plt.plot(gamma_list, known_worst_case_value, color='blue', linestyle='--', linewidth=2, label='Known')\n",
"\n",
"plt.legend()\n",
"plt.xlabel(r'$\\gamma$')\n",
"plt.ylabel(r'$\\frac{\\|x_1 - x_\\star\\|^2}{\\|x_0 - x_\\star\\|^2}$')\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"One can admit that the match between PEPit's results and the known convergence rate is pretty good. For completeness, let us also extract the dual variables. The next lines examplify how to do this extraction for given values of the parameters, and we integrate it in a loop for plotting below."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Constraint \"Performance metric 1\" value: 0.9999999999866084\n",
"Constraint \"Initial condition\" value: 0.8100001805827125\n",
"Constraint \"IC_Function_0_smoothness_strong_convexity(Point_0, Point_1)\" value: 1.7999998617137993\n",
"Constraint \"IC_Function_0_smoothness_strong_convexity(Point_1, Point_0)\" value: 1.7999998617138078\n"
]
}
],
"source": [
"verbose = 0\n",
"gamma = 1/L\n",
"\n",
"pepit_tau, list_of_constraints = wc_gradient_descent_distance(mu,L,gamma, verbose)\n",
"\n",
"nb_cons = len(list_of_constraints)\n",
"\n",
"for i in range(nb_cons):\n",
" print('Constraint \\\"{}\\\" value: {}'.format(list_of_constraints[i].get_name(),\n",
" list_of_constraints[i]._dual_variable_value))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"One can observe that the performance estimation problem is actually involving 4 constraints, 2 of which are of interest for us: the 3rd and 4th ones (that encode the fact the function is smooth and strongly convex). We can therefore extract their values in a loop, plot them, and expect their closed-forms to be nice. Perhaps one could even guess their values."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"verbose = 0\n",
"\n",
"gamma_min, gamma_max = -1, 3\n",
"nb_gammas = 50\n",
"gamma_list = np.linspace(gamma_min,gamma_max,nb_gammas)\n",
"\n",
"pepit_worst_case_value = list()\n",
"pepit_dual_value1 = list()\n",
"pepit_dual_value2 = list()\n",
"known_worst_case_value = list()\n",
"\n",
"for gamma in gamma_list:\n",
" pepit_tau, list_of_constraints = wc_gradient_descent_distance(mu,L,gamma, verbose)\n",
" pepit_worst_case_value.append(pepit_tau)\n",
" known_worst_case_value.append(max((1-gamma*L)**2,(1-gamma*mu)**2))\n",
" pepit_dual_value1.append(list_of_constraints[2]._dual_variable_value)\n",
" pepit_dual_value2.append(list_of_constraints[3]._dual_variable_value)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Playing a bit with candidate expressions, one can guess a closed-form for the multipliers (which happen to be always equal): $\\lambda=2|\\gamma| \\rho(\\gamma)$ with $\\rho(\\gamma)=\\max\\{|1-\\gamma L|,|1-\\gamma\\mu|\\}$."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"tags": []
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEICAYAAABS0fM3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAywUlEQVR4nO3dd3gU1dvG8e+TRoDQQaSIgCIoSDOAIk2kCUhHehPpICJSQpOm9CY9FOm9d4EfzQYICEpvAtJDJxBKkvP+sWteRAIpm51N9vlc117ZnZ3s3Jkk++w5c+aMGGNQSinlvjysDqCUUspaWgiUUsrNaSFQSik3p4VAKaXcnBYCpZRyc15WB4iJtGnTmqxZs1odQyml4pW9e/deM8ake3p5vCwEWbNmZc+ePVbHUEqpeEVEzj5ruXYNKaWUm9NCoJRSbk4LgVJKubl4eYzgWR4/fsz58+d58OCB1VFULPj6+pI5c2a8vb2tjqKU20gwheD8+fMkS5aMrFmzIiJWx1ExYIzh+vXrnD9/nmzZslkdRym34bSuIRGZLiJXReTgE8uGichREflDRJaLSMqYvv6DBw9IkyaNFoF4TERIkyaNtuqUcjJnHiOYAVR4atkmII8xJi9wHAiIzQa0CMR/+jtUKhKbNsF338G9ew5/aacVAmPMDuDGU8s2GmNC7Q93ApmdlUcppeINY6B3b+jYEV59FTZvdujLu9KooU+B9ZE9KSItRWSPiOwJCgpyYqyo8/T0JH/+/OTJk4fatWtz//79fy3/5zZ48GAASpUqRc6cOcmXLx/vv/8+x44dA6Bo0aIAnDlzhnnz5lnzwyilXMeOHbBrl+3+3buQO7dDX94lCoGI9ARCgbmRrWOMCTTG+Btj/NOl+88Z0i4hceLE7N+/n4MHD+Lj48OkSZP+tfyfW/fu3SO+Z+7cuRw4cIAmTZrQpUsXAH755RdAC4FSyuZ836lcxf6+17QpZMjg0Ne3vBCISFOgMtDAJKDLpRUvXpyTJ09Gef0SJUpErO/n5wdA9+7d+fHHH8mfPz+jRo2Kk5xKKRe3fz89t5XhVc6ykDrw1VcO34SlhUBEKgBdgSrGmPsOfOG4u0VBaGgo69ev5+233wYgJCTkX11DCxcu/M/3rF69OmL9fwwePJjixYuzf/9+OnXqFPv9opSKd872mcY86vMYbwpXTAM5cjh8G047j0BE5gOlgLQich74GtsooUTAJvtokZ3GmNbOyuRo/7zhg61F0Lx5c+D/u4aepUGDBiROnJisWbMyduxYJyVVSsULp04xcnUOQvGmPnPJ1v/TONmM0wqBMabeMxZPc9b2neF5b/iRmTt3Lv7+/nETSCkVr10bOIkp9AOg27s74J0GcbKdBHNm8b8kkEMNyZIl4+7du1bHUEpZ4coVxs5OSQhJqMha8n5TJ842ZfnBYnfw9DGCJ0cNPU/evHnx9PQkX758erBYKTcTPGwi48JsPeXdc66ADz6Is20lzBaBRYKDg5+5PCws7JnLt23b9tzX8fb2ZsuWLQ7JppSKR27f5vrkJRTiXe6SjGIDK0R5sEpMaCFQSilXM3kyrwYfYgMfcf/1vEj1fXG6Oe0aUkopV/LgATzRFZyk++fg6Rmnm9RCoJRSLiR85mzaX+7JTopAxozQsGGcb1O7hpRSylWEhbG6717GM4lVVOF0xyV4JUoU55vVFoFSSrkIs3gJ315uBkDnxBPwatPCKdvVQqCUUq4gPJzNAf9jN0VIx1U+65AEkiVzyqa1EDhQZNNNP2nbtm1UrlzZodvdtm1bxIyl0ZE1a1auXbvm0CxKqRhavZqBZ2xnDn/pPY6kXds5bdN6jMCBYjLFhCNs27YNPz+/iOsYKKXiGWPY0W0tOwgkJTdp28ZAmjRO27y2CJxgw4YN5MqVi4IFC7Js2bKI5X379mX48OERj/PkycOZM2cAmDVrFnnz5iVfvnw0atQIsM1QWqRIEQoUKECZMmW4cuUKZ86cYdKkSYwaNYr8+fPz448/EhQURM2aNSlUqBCFChXi559/BuD69euUK1eO3Llz89lnnxHZrN/Tpk3jjTfeoHDhwrRo0YL27dsD0LRpU5YsWRKx3j/TZQMMGzaMQoUKkTdvXr7++msA7t27R6VKlciXLx958uSJmHW1e/fuvPXWW+TNm5ev4mBKXaXinY0bGXqsCgAdPceTPMB5rQFIwC2C552EN3kytGxpux8YCK1aRb5udKYtenL2UYCAgACqVq1KixYt2LJlC6+//jp16rx4vpBDhw4xcOBAfvnlF9KmTcuNG7YrfBYrVoydO3ciIkydOpWhQ4cyYsQIWrdujZ+fX8Sbav369enUqRPFihXj3LlzlC9fniNHjtCvXz+KFStGnz59WLt2LdOm/XfOv4sXLzJgwAD27dtHsmTJKF26NPny5Xtu3o0bN3LixAl2796NMYYqVaqwY8cOgoKCyJgxI2vXrgXg9u3bXL9+neXLl3P06FFEhFu3bkVx7yqVQBkDAwYwmb8YzRd83uwuvPyyUyMk2EJghWd1De3fv59s2bKRwz6HeMOGDQkMDHzu62zZsoXatWuTNm1aAFKnTg3A+fPnqVOnDpcuXeLRo0dky5btmd+/efNmDh8+HPH4zp07BAcHs2PHjogWSaVKlUiVKtV/vnf37t2ULFkyYpu1a9fm+PHjz827ceNGNm7cSIECBQDbFBknTpygePHidO7cmW7dulG5cmWKFy9OaGgovr6+NG/enMqVKzv8eIlS8c6OHfDzz2QChnn3hN5Rv6CVoyTYriFjIr/90xoA2/3nrRuXvLy8CA8Pj3j84MGD567foUMH2rdvz59//snkyZMjXT88PJydO3dGXBrzwoUL/+rGcUTe8PBwHj16BIAxhoCAgIjtnTx5kubNm/PGG2+wb98+3n77bXr16kX//v3x8vJi9+7d1KpVizVr1lChQoVY51IqPrv59WjC/nkrbtIEsmRxeoYEWwhcRa5cuThz5gynTp0CYP78+RHPZc2alX37bHOI7Nu3j7/++guA0qVLs3jxYq5fvw4Q0TV0+/ZtMmXKBMDMmTMjXufp6arLlSv3r4vc/NNKKVGiRMQ1kNevX8/Nmzf/k7dQoUJs376dmzdvEhoaytKlS/+Vd+/evQCsWrWKx48fA1C+fHmmT58eMVnehQsXuHr1KhcvXiRJkiQ0bNiQLl26sG/fPoKDg7l9+zYVK1Zk1KhRHDhwIHo7VKmEZOdO2myvw1scZqe8B1GcmdjRtGvIgZ4+RlChQgUGDx5MYGAglSpVIkmSJBQvXjziTbtmzZrMmjWL3LlzU6RIEd544w0AcufOTc+ePSlZsiSenp4UKFCAGTNm0LdvX2rXrk2qVKkoXbp0ROH4+OOPqVWrFitXrmTs2LF89913tGvXjrx58xIaGkqJEiWYNGkSX3/9NfXq1SN37twULVqULM/45JEpUyZ69OhB4cKFSZ06Nbly5SJFihQAtGjRgqpVq5IvXz4qVKhA0qRJAVvhOXLkCO+99x5gO4g8Z84cTp48SZcuXfDw8MDb25uJEydy9+5dqlatyoMHDzDGMHLkyDj7fSjl6o4FzGARE/AilMzV/OG11yzJIfHxevH+/v5mz549/1p25MgR3nzzTYsSJSzBwcH4+fkRGhpK9erV+fTTT6levbrTtq+/S+UW9u+naYH9zKQpLQlk8qFi8NZbcbpJEdlrjPnPJRG1a0j9R9++fcmfPz958uQhW7ZsVKtWzepISiU4f/Wcwhwa4kko3SociPMi8DzaNaT+48lzG5RSceDIEYasy0sYXjRiFtkHOWdOocgkqBZBfOzmUv+mv0PlDv7uPp7pNEMIJ6Dkr/DEsUUrJJhC4Ovry/Xr1/WNJB4zxnD9+nV8fX2tjqJU3DlyhNOrDvISV6nDQt4c2szqRAmnayhz5sycP3+eoKAgq6OoWPD19SVz5sxWx1Aq7gwYQEm2c4rXuPNhDSg8/8XfE8cSTCHw9vaO9ExbpZRyCYcPw4IFACTiEem+7WRxIBundQ2JyHQRuSoiB59YllpENonICfvX/855oJRSCcTZ7hPpY/pyg1RQqRIULmx1JMC5xwhmAE/PJ9Ad+J8xJgfwP/tjpZRKeA4fZtDqPAygD18xHOyz9LoCpxUCY8wO4MZTi6sC/8yVMBOo5qw8SinlTGe7TWA6zfAgjK4ld0OhQlZHimD1qKH0xphL9vuXgfSRrSgiLUVkj4js0QPCSql45dAhvl2Tl8f4UI/55BrW3OpE/2J1IYhgbOM+Ix37aYwJNMb4G2P806VL58RkSikVO2e7T4xoDfQq9bNLtQbA+kJwRUQyANi/XrU4j1JKOZa9NRCKt0u2BsD6QrAKaGK/3wRYaWEWpZRyuIsBYyNaA71L/QT+/5nzzXJOO49AROYDpYC0InIe+BoYDCwSkebAWeATZ+VRSqk4d/AgGVYHsp6T/EYhcg77zOpEz+S0QmCMqRfJUx86K4NSSjlV374IhjL8jzIfJwH/QVYneiaru4aUUiph2rePC0t//f/HLnTewNO0ECilVBw48sVksnKGxszEVK8B77xjdaRIaSFQSilH+/ln+vxYhlC8SUIIMqC/1YmeSwuBUko5kjH8/vn3LKE2iXhA7xqHIHduq1M9lxYCpZRypM2b6bXPdo3vdjKRTEM7WhzoxbQQKKWUoxjDL58vYB2V8OMu3Rueh9deszrVC2khUEopBzErV9HzaEMAvvAc5zLXG3gRLQRKKeUI4eE86tWf7JwmLUF0bnEH4snV9hLMFcqUUspSCxeS6NA+pvEZd5OkJ1m/P6xOFGXaIlBKqdgKDf3XCWPJvmgOL71kYaDo0UKglFKxFP79TOqf6MtqKmOSp4CvvrI6UrRo15BSSsVGSAgLu//OfMbxI8Up9+VcEqWKX5df1xaBUkrFwqPRE+h540sAvk42ikSd21ucKPq0ECilVEzdvMnkAVf4i+y8yWGafpMD/PysThVtWgiUUiqG7vQbRf+QLgAMevk7vFq75vUGXkQLgVJKxcTffzNiXCKukY73+Ykqo0uDt7fVqWJEC4FSSsVAaJ/+zAmrC8CQnN8jtWtZnCjmtBAopVR0HTqE16zpHCAf86jH+xMagEf8fTuNv8mVUsoqPXpAeDh+3KNe+ZtQurTViWJFC4FSSkXHTz+xeJUPwSQFERg82OpEsaaFQCmlosoYdrf5nk9YTF7+4FHdxpA/v9WpYk0LgVJKRZFZsZJuB23TTNfxWILPt32tDeQgWgiUUioqQkPZ8Pk6tvEBqbhBt5Y3IWtWq1M5hBYCpZSKgtCJU+hy/nMAeiQaQcoBnS1O5DguUQhEpJOIHBKRgyIyX0R8rc6klFIRbt9mWsBJDpGHbJymfc+UkDat1akcxvJCICKZgM8Bf2NMHsATqGttKqWU+n/3+g6jz72uAAxJOxzfLh0sTuRYrjINtReQWEQeA0mAixbnUUopm7/+Isn4YUzndxZTm1pjioNvwuq0sLwQGGMuiMhw4BwQAmw0xmx8ej0RaQm0BMiSJYtzQyql3FdAAPL4EZVYR6Ui16Her1YncjhX6BpKBVQFsgEZgaQi0vDp9YwxgcYYf2OMf7p06ZwdUynljn79lfMLf/r/xyNH2k4iS2AsLwRAGeAvY0yQMeYxsAwoanEmpZS7M4bdLaeSlTN0ZDTUrg1FE+ZbkysUgnPAuyKSREQE+BA4YnEmpZSbMwsX8eXBZoThRRKPhwliKonIWF4IjDG7gCXAPuBPbJkCLQ2llHJvDx6wtMM2fqYY6bhKQLs7kD271anijOUHiwGMMV8DX1udQymlAB6OGEfXa7Yrjw1IOoTk/XtbnChuWd4iUEopl3L1KmMH3OIvsvMWh2j+7WuQMqXVqeKUFgKllHrClU6DGfDQ1hoYnmk0Xm1aWJwo7rlE15BSSrmEPXtg3jwq8w63ScFHgdXj7XWIo0MLgVJKAYSHw+efk54rzKUhjz6qChVXWJ3KKbRrSCmlgPDZcwn7dZftgbc3PmOGWRvIibQQKKXU3bvM7bgbf/bwC+/Bl19CjhxWp3IaLQRKKbd3t88wut7uwX4KcDKFP/TsaXUkp9JCoJRybydOMPC75FwmA0XYScMxhSBZMqtTOZUWAqWUWzvRchijwm1XHhubJxCPRg0sTuR8WgiUUu5r/Xo6bavCY3z4lOkUmt4GPNzvbdH9fmKllAJ49Ih1ny1jLZVJzm2+rX8QChWyOpUl9DwCpZR7Gj2aOxfvkoob9Eo0nPSjuludyDJaCJRS7ufcOejXj7rcpyybSDagD7z0ktWpLKOFQCnldkzHL5D79wFI83Ym+KKtxYmspccIlFJuxaxZS50VdRlCVx7hDRMmuMV8Qs+jLQKllPsICWFF89UsZhI/UJ4mdR7ycrFiVqeynLYIlFJuI7jfCDpe7QHAN4m/4eWx7nUGcWS0ECil3MPx4wwY5svfZKEA+2gzMgekS2d1KpegXUNKqYTPGA41GcrI8IkI4Ux8axyeLadancplRLtFICJJRcQzLsIopVRcMAsX0XZnI0LxpiVTKDK7vVueQRyZF+4JEfEQkfoislZErgJHgUsiclhEhonI63EfUymlYujOHW590ZcwPElLEN+2+AsKFrQ6lUuJStfQVmAzEAAcNMaEA4hIauADYIiILDfGzIm7mEopFUO9e5PqylF2UILTaYuQetgGqxO5nKgUgjLGmMcikumfIgBgjLkBLAWWioh7D8JVSrmmXbtg7FgAPDC8PqYDpEhhcSjX88KuIWPMY/vdNSLSV0QSP2cdpZRyDY8fs63eZCqaNZzhVShbFurVszqVS4rO0ZJCwG1gl4g0dmQIEUkpIktE5KiIHBGR9xz5+kop9xPy7Sha/tWd9VRkjnczmDwZRKyO5ZKiXAiMMaHGmFFASeAdEflFRIo7KMcYYIMxJheQDzjioNdVSrmj48cZOMBwgjd4i0N0HZAcsmWzOpXLivJ5BCKSHSgP5LTfXge+tx8fOGOMKRmTACKSAigBNAUwxjwCHsXktZRSCmP4o8EQhoZNQghnSs4R+HQOtDqVS4tO19D/gBT2rx2BjMaY140xrwKx6SrKBgRhKyq/i8hUEUn69Eoi0lJE9ojInqCgoFhsTimVkIVN/Z4We1oSijdtZBJF53cALz139nmiUwjKGmMGG2NWG2OOGWNC/3nCGHM2Fhm8gILARGNMAeAe8J8rRBhjAo0x/sYY/3R6WrhS6lkuX2Zcx+PspgiZOM+gDpegQAGrU7m8qJxQJgDGmJMvWieGzgPnjTG77I+XYCsMSikVPR07EhIC3jxiQvr+JB8UYHWieCEqLYKtItJBRLI8uVBEfESktIjMBJrENIAx5jLwt4jktC/6EDgc09dTSrmp1ath0SK6M4STvE6V2bUhSRKrU8ULUek4qwB8CswXkWzALcAX8AQ2AqONMb/HMkcHYK6I+ACngWaxfD2llDu5dYvw1m0jPtlmafyB7bwBFSUvLATGmAfABGCCfYRQWiDEGHPLUSGMMfsBf0e9nlLKvVxu3ZcPLm5kAL2plW4HjBhhdaR4JVrT7xljHhtjLjmyCCilVGyYtetos7AkR3mTqXyGGT8B0qa1Ola8Eq0xVSLyCpAbyAO8DeQ2xugneaWUNW7dYkHDNaxgAsm5zZTKq5DaE6xOFe9EZdRQK/tZxLeA48BngB+wCqgft/GUUipyl1v3pf2tAQCM9PuaV6b3szhR/BSVFkEAUAe4BgwGEgPTjTHn4jKYUko9zz9dQjdIQ3k28On0YnrpyRiKyjGCysaYXcaYU8aY2sB4YLWIdBIRvcSPUsr5bt1ifsO1rKD6E11CtaxOFW9FZRrqg089Xg8UBlIDP8dRLqWUilynTrxy6w+yc0q7hBwgRhNwGGMeAr1FZLaD8yil1POtWwczZlAc+JO3Sfz9bO0SiqVYde0YY447KohSSr3QjRtc/fT/pyJLUqcKUqumhYESBu3jV0rFD8ZwrnEvcl7ZTjvG8ThtBhg3zupUCYIWAqVUvBA+ey5N1tbmFqk4T2a8pkzUE8ccRAuBUsr1nT3LyJZH2cYHvMQVpjTYjlSranWqBEOv1qCUcm1hYfxRoy89H04CYPrLPXlp0mhLIyU02iJQSrm0B4NH02DflzwiEa1lEpWWfwZ+flbHSlC0ECilXNfvvzO4z30O8jY5OM7w7tfh3XetTpXgaNeQUso1hYRAw4Z0DL/ISbLx+Vv/I2k/vQh9XNBCoJRyTQEBcPgwqYA5SVrBiv3g7W11qgRJu4aUUi7HrN/A7DHXeYiPbcGoUZAjh7WhEjAtBEop13LxItM++YHGzKYMmzGVKkOLFlanStC0a0gp5TrCwjhUvRefB9vOGG6TcgHy/XQQsThYwqYtAqWUy7jfdyif7O5MCEloxvfUX15bJ5RzAi0ESinXsGMHHQem4zC5ycURxgZchFKlrE7lFrRrSCllvWvXWFBtAVOZQCIesLDgUJL2n2J1KrehLQKllLWMgWbN+OFmIQBGJ+1F3pUDwEs/pzqLFgKllLVGj4Y1a5jOp6ylIq3ml4LMma1O5VZcphCIiKeI/C4ia6zOopRykt9+I6xrAAACVOyUC/m4srWZ3JDLFAKgI3AkTrdw5w4sWBCnm1BKRdG1ayytOI1Cob9wmmzwzjsweLDVqdySSxQCEckMVAKmxtlGDh+GQoV4XK8RZsXKONuMUioKwsI4Vq0bza4N5XcKsta3FixcCD4+VidzSy5RCIDRQFcgPLIVRKSliOwRkT1BQUHR30Lv3lw+fpvSbGFE3d1w7FiMwyqlYie45yBq/Pwld0lOHRbQfmFxeO01q2O5LcsLgYhUBq4aY/Y+bz1jTKAxxt8Y458uJieYBAby20uV+YnidHvYny3lh8DduzFMrZSKKbN6DS2GvMZhcvMmh5n61TGkysdWx3JrlhcC4H2gioicARYApUVkjsO3kiYNH29oRw/PIYTjSZ2zQzhXt6tt6JpSyjlOn2ZsnR9ZQD38uMuy94bjN7iX1ancnuWFwBgTYIzJbIzJCtQFthhjGsbJxgoUoP+0TJTjB66RjprrPuXB4NFxsiml1FNCQjhaqTOdQwYC8H2aLuRaOQQ8PS0OpiwvBM7m2aQh85pvISt/sYdCdOiRDLZssTqWUgmbMdC2LTmPrmAUnejmMYxaa5vpPEIuwqUKgTFmmzEmzgcRp5kwgGV5++FLCFP5jHXVp8C5c3G9WaXcV2AgzJiBAO0Zz+BxflCkiNWplJ1LFQKn8fGhwPpvCUzRlf70psKdhVCzJjx4YHUypRKeHTsY0vYMJ3jd9rhxY2jd2tpM6l/csxAAZMxIozV16O01GA8M7NkDbdrowWOlHOnMGaZUWk738EEU50fuvf0uTJyo1xdwMe5bCACKFbNdAg+4QEZazyjCgyFjLA6lVAJx9y47PuxH2+ChAAxO9i1J1yyEJEksDqae5t6FAKBdO0yTptRiCZNpTcuA1JjVOt2RUrESHs6Zmp2peXoooXjzpccomm6oC1myWJ1MPYMWAhFk8iQm5Z9MUoKZTWOG1PoN/vzT6mRKxVvBAd9QdVM7rpGO8mxgyORUULSo1bFUJLQQACRKRL6Nw5iT7kuEcAIe9WP5h+Pg6lWrkykV75gFC2k8NDd/kI+cHGVBmx14fdbU6ljqObQQ/CNdOqpt+ZxBPn0BaBg0kt/LdYOHD63NpVR8smcP8mkzPmI96bjKqmLDSPldf6tTqRfQQvCkPHnouqQwjZnJfZJS5UB/bjb5QkcSKfUixmAmToISJSAkhBZM5dTrFXhj9Qi90lg8oIXgKfJxZQIH36QYP9KecaRcOAmGDrU6llKu68oVTpRuxftt87I1xH6SWIoUJFszH1KmtDSaihotBM+QqGtHtjadRTeGIgDdu9vmSldK/YtZsZIprw8m/7ZR/EpRejMA81Zu2LEDcua0Op6KIi0EzyKC1+TxtmYucJpsDGuwH7ZtszSWUi7j7l2uNvySatWhZfAo7pOUesxjdet1yN49kDev1QlVNGjnXWR8fGD5ch4ULU2pY6v4OywLXh/1oNPuNPD221anU8r5Tp6E9esJW7+RYZvyMyq0G1dJTwpuMSFVL+ovqgZlvrU6pYoBMfHwQKi/v7/Zs2ePczZ29izz8g+lwa3xACxM1ZpP/ugFmTM7Z/tKWSUkBLZt4+ayrfy89haVL02JeCoXRzhGLkqxlZmVFpNl1kBIndrCsCoqRGSvMcb/P8u1EETBgQMMLbyEbo8G4MNDNmZrTcl9o/RAmEp4zpzBrF3HHwsOs3ZnataFluNX3iMcT86ShSz8DcBCPsEn00tUHfQuHg3r69xB8URkhUC7hqIiXz66rAni7/LjGWfaUfWvUfxUthN5fpoEiRJZnU6pmAsNhV9+gdWrubxqN32ON2AdVbhA24hVvHhMCXZw0zcjWcrmhwoVqFOhAmTPbl1u5VDaIoiGsFlzqd0kMcupQWb+5mD1PqRYMg089Ji7ikdu3oQNGzgz/1eObblA+XvLAAgmKWm4ziMSkYGLVGQdFTP9QZlqfiSvVhqKF9cPPvGctggcwLNxA+aeG0H53juozWJSLJ8B7XxhwgRtGivXduIEYStWs3PuKVb/8SprTEUOUY+U3CSIlXgRhh/3mO7dmrcKJSV/3VxIpYqQ/TOrkysn0EIQTYl7fsnWSx3xnDDOtmDSJEiWDIYM0WKgXEdYGOzaBatWcWTxQYacrsVaGnGN/780ZDLuUJZN3MyQm3RVi0KlSjQoXVqniXZDWgiiSwTPsaPh9g2YO5ej5KTvsIJMTzSYJAMCrE6n3FlICGzezKV5W7m8YT8Fbm0F4DFvM5OmAGTnFB+zho9znaB4nYz4VKsI+fbrhxg3p4UgJjw8YMYMTPA9Gq7syV78uTnwB1alGEuirzpYnU65kxs3YM0ajs7Zw8ptKVjxuCI7GUkRdrKT9wB4mz8Z7f0VZYs94M0GBZHKdSF9eouDK1eihSCmvLyQhQuYW7oNJX4ZxEbKU6/LPRYlm4ZXq+ZWp1MJ2YULsGIFR2f/xozdb7HCVOEYjSOe9iWEl7lMaPpMeFWpiFSpQscPP4TEiS0MrVyZjhqKrXv3OFCsHaX2j+IWqWjELGbM9cGjfl2rk6mE5MQJQhcv586yzaTeuwmAedSjAfMASMUNPmY11TLtoVydVCStXREKF9YRbepf9ISyuHT7NjuLdKTMsXHcw4/WMpkJi9MhNWtYnUzFV8bAoUM8XLiCTbMusezcO6ykKtVZzlRaAHCb5PShP9XePE7xhq/iVaMK5MplcXDlyrQQxLVr19jq34WPzk7kIb5s9/iAEgvaQu3aVidT8YUxsH8/9+evZMPc6yy9+C6r+Zi7JI9Y5X1+4ifPUvDBB1CjBlStChkzWpdZxSsuex6BiLwCzALSAwYINMaMsTZVDKRNywe/fsvKQu04fcGHEuHboN6PtmF8dbWbSEXCGNi3DxYtgiVL4PRpvmUA3zA2YpX8/E5Nr5XULHGNN5sUhspXdV4f5VCWFwIgFOhsjNknIsmAvSKyyRhz2Opg0ZYhA+V/GwilS8NRICyMs/UDyPw4HM9G9a1Op1yFMbB3L8FzV7J27i0WB5WkLLdoxWkAarCMTZSlpvdqapa9w2tNi8NHX4Gfn8XBVUJleSEwxlwCLtnv3xWRI0AmIP4VAoAMGWzXLShdmiOHwylltlGu8UZmhM7Cs1njF367SqDs3T73Zi9jzRzbm/86AgjBdvLWZV6mFYGQPDkFq+RmV81LUL6njvRRTmF5IXiSiGQFCgC7LI4SO+nTw9atXHuvK/dOJ2UOjQj7dB6zQr/Hq0Uzq9MpZzp0yHZ1u4ULGXq8Kv34mvskjXj6PX7hk0SrqFn5ITRdDWXL6nw+yulcphCIiB+wFPjCGHPnGc+3BFoCZMmSxcnpYuCllyi+azg/vNeej05+x3zqE9ZyIXMeTMK7Q2ur06m4dOoUj+YsYuP082Q9t508HALgZS5zn6S8y6/USbSSmpUf8krTD6FsP33zV5ZyiVFDIuINrAF+MMaMfNH6LjlqKDI3brCzaCfKH/uOO6SgAutZ0uN3kg4M0NP6E5KLFwmbv4jtU44z/1gBllKTm6SmFZOYRBsAgpOm53qZOrzavAyUK6dv/srpXHnUkADTgCNRKQLxTurUvPvraLaU+JwKB4ezgY8o921ydtzogOe4MeDpaXVCFVM3bsDSpfw+eTcz9+ZhEZ9wif8fyvk2f/C29zGoVhvq1MGvYkX8tM9fuSDLCwHwPtAI+FNE9tuX9TDGrLMukoOlSsU7v47n5486Uu6n3jRmFp6TAuFGEMyapZ8M45OQEFizhvA58/BYvxYeP2YuwxhDR8A2qVs9j0XUK3Ge3C3ehyoDdLSPcnmWFwJjzE9Awu8j8fPjjf9N5GCDVvgtmWFbtmgRoddu4bViiW0qa+WawsJg61YuTlnLgpW+zH1Yi+ZkoC2PAWjMLMLxpK7/CQq1LIjUbKXj/FW8YnkhcCs+PvgtnAadksN333GQ3FTfMo45hTtQZMcwSJfuxa+hnOfAAe5OX8yyWcHMvlWZLYzAYJu7Jw3XactEKFSIvPXrM7JOHdvQYaXiIS0EzubhAaNHQ/r0jOmZjpPkoPTR8SzM+wWVt3bWuWKsduECzJ0Lc+Yw+M+K9KdPxFh/Hx5SibU0eHkLlZq/DI2PwRtvWBxYqdjTQmAFEejRgwmppvK47Uxm0oSqlycxvEBvvlhVGilbxuqE7uXePcyy5fw+7mdS7N7Ea5wC4GUKEkISSrCdhslWUaueN6ma14BCY3XEl0pQXGL4aHTFq+GjL2BWrGTAJ3/w9ePeALSUQMaNBe92LS1OlsCFh8P27VyctIq5K5Iy81FdDpGHNkxgAu0AuOebhqCy9cna5iPbiV5e+rlJxW8uO3zU3Um1qvTZlYWcH7ah6c2RBJqW/N1+HWuPdURGjtA3H0c7dYr7U+exYmoQM69VYjPDCcc2hDctQaThBnz4ITRqRNIaNUiqB/GVG9B3GVdQoAB1DvYm64dtqHZ0EPWZh4ydCydPwIIFkDz5i19DRe7uXVi8GGbMgB9/pBcjGIWtBebNI6qxgiaZ/keF1lnxadIEXullbV6lnEy7hlzJ/fvcqd+a5CtnRywKylGUdKum6UHk6LJ3/Zwft4LZq1OS8/Gf1GA5APsoQBsm0jTpEj6p60GaVrXA31/7/VWCp11D8UGSJCRfNgP6ZIFvvmE/+Sh24gd65RtG15m58aj7idUJXd+5czycNodVky4w/WplNjKScDx5n59shcDTk4KVXmFX04tQ6Rvw8bE6sVKW00Lgajw8YOBAyJWLH5vt416oHwGP+vFLvVXM3BZAqu/66ZvX0x4+hJUrOTJmIxN/ycdcWnGDNIBtyGcVltHs1a3w+Qho0MA2O6xSKoJ2DbmyAwdYW2EsjS4P5SapycZplubtT4G1AyFzZqvTWe/PP2HaNJg9G27cYCrNacFUAPKxn+aJ51G/HqRpWwcKFtSuH+X29JrF8dWtW5z5pCu1NrVkL/4k4gHj/brTfMXHttEt7ubuXcyChewYtZepR4qSiQsMJgCAOySjNwNpWuQIBb4oCdWqga+vtXmVciFaCOIzY3gwaBQde/kRaFqSlGCOk5OM3RpB//4Jv6vIGNi9m8tjFjJjiR/THjfiJDkASMM1LpEB71czQbNm0LQpvPqqtXmVclFaCBKC7duZWWUpvneuUIdFtmX589umRHjrLUujxYlbt2DuXP4Ys5W+J+qziiqE2Q9rZeI8zTxm8WmFi2T7oqqtdeThYW1epVycjhpKCEqWpMnRN6BxY9hsWzRx/7ucz7uar4dtxadjm/j/ZmgM/PoroZOm4rVkAYSE8IBCLKcGnoRSjeW0yLKR8h1z4dm4JaRNa3VipeI9bRHER+HhMHYsN7oOJsujE9zDj3fYw5z3J5FrUX/ImPHFr+Fqbt0ibNZcfhhxkMnnKnCfJGyiHAAGmOzdgao1PMnQ8RN491098KtUDGjXUEJ08CA/VhtBo1Nfc5asJOY+Q5P0o+343Hg0aeT6b5bGwG+/cWnUAqYvSc6U0KacJStgO+P3HFl4OX8GaNkS6teHFCmszatUPKeFIKF6+JA73b7h8zHZmUlTAIryM5OLfE+e2d0gRw5r8z1LcDDMm8eZMSv56nAzVlKVULwByMZpWnrPoFntYNJ3qg/vvOP6BU2peEILQUK3fTvLas+jXVBfLpOBUmxla6KPoE8f+Oor1xhZdPgwoeMn4zVnBty5w01SkpGLPMabKqyiVbZNlO2cF49GDXR+JaXigBYCdxAczK1ug+gxITPtGEduDgPw+K18eE+ZAEWLOj/To0eYZcv5bcgWJux/jx2U4Ci58LFf5nGldy0KVc9Mxk51oEgR/fSvVBzSQuBO9u2z9avv3YsBPmY1abjO4Nr7yDCqK2TKFPcZzp/n/rjpzJ94kwl3GrKPdyKe2kopSuW8DK1bQ5MmkCpV3OdRSmkhcDuhoTB2LMd6ziJfyK88xJekBNPVaxSdO0PSXp3Az8+x2zQGtm0j5Lsp9F7pz3TTlJvYLuKemut8KjNoVf4Mr3epDh98oJ/+lXIyLQTu6tw5Tjb7hi5bKrCC6gBk5ALfJB9Ko+H58Py0CXh6xm4bd+9iZs1GJoyHw4cxwJsc4Ri5KMwu2qWcxyftX8K3TbP4ObRVqQRCC4G727KF7a3m0flka/Zi+zt4j1/4KU8bPAZ/CxUrRv8T+okT3Bw+jekzPZn8sAmbKUMW/gZgMx+SsuBr+PcoB1WqgLe3o38ipVQ0RVYI4vlpqCrKSpem5NHJ7J52kNkpO5CZvynFNjwO/gGVKxNe0B+zdJntZLXnCQ+Hdes4WKw1rd7YQubA3nz18BtO8AZzaGjrbmrbljIHx+C/dzLUrKlFQCkX5xItAhGpAIwBPIGpxpjBz1tfWwSxdO8eIYPHED5iFElDrgEQSAtm04jer86m7KDSyCe1bV1GDx7AgQOwaxfs3s3ajd4MD2rMNj6IeLkybKJDpmVU6poHz6aNdOinUi7KZbuGRMQTOA6UBc4DvwH1jDGHI/seLQQOcukSDB0KkyfzTsiPESN7CrOLHhlm4JvMm90nU/NZ+GQycBmA5kxlOs3x4y5NmEW7kgd5s0d1KFMm/s9zpFQC58qF4D2grzGmvP1xAIAxZlBk36OFwMGuXOHOoPFMnGAY8bgDQbz0r6cXU4taLAVgA+U5kTgvjZv7kKLTp5A9uxWJlVIx4Mqzj2YC+xFGm/NAEYuyuKf06Uk+uj/del+n/dDxBI4JYerDhqTkFkXYxWuvPIYSDaBIESoULkyF/PkhUSKrUyulHMQVCkGUiEhLoCVAlixZLE6TQKVJQ9IhfegUcItOW7ZAkiRQuCmk/tLqZEqpOOQKheAC8MoTjzPbl/2LMSYQCARb15BzormplCmhRg2rUyilnMQVju79BuQQkWwi4gPUBVZZnEkppdyG5S0CY0yoiLQHfsA2fHS6MeaQxbGUUsptWF4IAIwx64B1VudQSil35ApdQ0oppSykhUAppdycFgKllHJzWgiUUsrNWT7FREyISBBwNobfnha45sA4jqK5okdzRY/mih5XzQWxy/aqMSbd0wvjZSGIDRHZ86y5NqymuaJHc0WP5ooeV80FcZNNu4aUUsrNaSFQSik3546FINDqAJHQXNGjuaJHc0WPq+aCOMjmdscIlFJK/Zs7tgiUUko9QQuBUkq5uQRfCESktogcEpFwEYl0yJWIVBCRYyJyUkS6OyFXahHZJCIn7F9TRbJemIjst9/ibHruF/38IpJIRBban98lIlnjKks0czUVkaAn9tFnTso1XUSuisjBSJ4XEfnOnvsPESnoIrlKicjtJ/ZXHydkekVEtorIYfv/YsdnrOP0/RXFXFbsL18R2S0iB+y5+j1jHcf+PxpjEvQNeBPICWwD/CNZxxM4BWQHfIADwFtxnGso0N1+vzswJJL1gp2wj1748wNtgUn2+3WBhS6SqykwzoK/qxJAQeBgJM9XBNYDArwL7HKRXKWANU7eVxmAgvb7yYDjz/g9On1/RTGXFftLAD/7fW9gF/DuU+s49P8xwbcIjDFHjDHHXrBaYeCkMea0MeYRsACoGsfRqgIz7fdnAtXieHvPE5Wf/8m8S4APRURcIJcljDE7gBvPWaUqMMvY7ARSikgGF8jldMaYS8aYffb7d4Ej2K5V/iSn768o5nI6+z4Itj/0tt+eHtXj0P/HBF8IoigT8PcTj88T938Q6Y0xl+z3LwPpI1nPV0T2iMhOEakWR1mi8vNHrGOMCQVuA2niKE90cgHUtHcnLBGRV57xvBWs+JuKqvfs3Q7rRSS3Mzds78IogO1T7pMs3V/PyQUW7C8R8RSR/cBVYJMxJtL95Yj/R5e4ME1sichm4OVnPNXTGLPS2Xn+8bxcTz4wxhgRiWwc76vGmAsikh3YIiJ/GmNOOTprPLYamG+MeSgirbB9SiptcSZXtg/b31SwiFQEVgA5nLFhEfEDlgJfGGPuOGObUfGCXJbsL2NMGJBfRFICy0UkjzHmmcd9HCFBFAJjTJlYvsQF4MlPkpnty2LleblE5IqIZDDGXLI3ga9G8hoX7F9Pi8g2bJ9aHF0IovLz/7POeRHxAlIA1x2cI9q5jDFPZpiK7diLK4iTv6nYevKNzhizTkQmiEhaY0ycTrAmIt7Y3mznGmOWPWMVS/bXi3JZtb+e2OYtEdkKVACeLAQO/X/UriGb34AcIpJNRHywHXyJsxE6dquAJvb7TYD/tFxEJJWIJLLfTwu8DxyOgyxR+fmfzFsL2GLsR6ri0AtzPdWPXAVbP68rWAU0to+GeRe4/URXoGVE5OV/+pJFpDC294A4Lej27U0DjhhjRkaymtP3V1RyWbS/0tlbAohIYqAscPSp1Rz7/+jMo+FW3IDq2PobHwJXgB/syzMC655YryK2UQOnsHUpxXWuNMD/gBPAZiC1fbk/MNV+vyjwJ7bRMn8CzeMwz39+fqA/UMV+3xdYDJwEdgPZnfT7e1GuQcAh+z7aCuRyUq75wCXgsf3vqznQGmhtf16A8fbcfxLJiDULcrV/Yn/tBIo6IVMxbAc7/wD2228Vrd5fUcxlxf7KC/xuz3UQ6POMv3uH/j/qFBNKKeXmtGtIKaXcnBYCpZRyc1oIlFLKzWkhUEopN6eFQCml3JwWAqWUcnNaCJRSys1pIVAqlkQkh4icEZHX7Y+97XPXu8oEeEo9lxYCpWLJGHMC2wXFy9sXtQdWGWP+jvy7lHIdCWLSOaVcwEGgjIikxjatQxGL8ygVZdoiUMoxjmO7El5fYLgx5p61cZSKOp1rSCkHsE9nfBHbpGlFjTHhFkdSKsq0RaCUAxhjHgN3sF2HWouAile0ECjlON7AdqtDKBVdWgiUcgD7NW/PGu1rVfGQHiNQSik3py0CpZRyc1oIlFLKzWkhUEopN6eFQCml3JwWAqWUcnNaCJRSys1pIVBKKTf3f6ztrrh9DDTWAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"candidate_lambda = [ max(abs(1-gamma*L),abs(1-gamma*mu))*2*abs(gamma) for gamma in gamma_list ]\n",
"\n",
"plt.plot(gamma_list, pepit_dual_value1, color='red', linestyle='-', linewidth=3, label='PEPit')\n",
"plt.plot(gamma_list, candidate_lambda, color='blue', linestyle='--', linewidth=2, label='Educated guess')\n",
"\n",
"plt.legend()\n",
"plt.xlabel(r'$\\gamma$')\n",
"plt.ylabel(r'$\\lambda(\\gamma)$')\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Section [4](#example4) below explores how to transform that into a formal mathematical proof (that include replacing the actual guessing game by formal mathematical steps). In the meantime, we explore a case where things are directly less convenient, that of function values. In fact, this slightly more complicated problem already show a typical feature appearing in computer-aided analyses: non-uniqueness of the proofs. This might seem like a good thing (in fact, it is, in many aspects), but it actually renders the formal translation and search more painful, as the numerical solvers actually generally do not provide the simplest proofs possible, but rather combinations of them and fail providing the *sparse* (and likely simpler) ones."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 6. Using symbolic regression (PySR) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Up to this point, two complementary approaches have been used to obtain the Lagrange multipliers needed to reconstruct a proof:\n",
"\n",
"1. **By inspection** — guess the functional form from numerical data \n",
"2. **Symbolic exploitation of the LMI** — exploit the LMI to solve for the multipliers symbolically \n",
"\n",
"This section focuses on automating the “inspection” step in (1).\n",
"\n",
"The task can be framed as **symbolic regression**: learning a function from measurements without restricting it to a prescribed functional form.\n",
"\n",
"Symbolic regression is hard in general, but effective heuristics are often available in practice. This demo uses [PySR](https://github.com/MilesCranmer/PySR). \n",
"Install it first (note: PySR’s backend is implemented in the Julia programming language)---the installation might take some time; so the first run might seem unnecessarily long."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, imagine that you do not know the closed form of $\\lambda_1(\\gamma)$ given at the end of Section 3. The next cell automatically learns the closed form of $\\lambda_1(\\gamma)$ using measurements from PEPit:"
]
},
{
"cell_type": "code",
"execution_count": 64,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"import numpy as np\n",
"from pysr import PySRRegressor\n",
"\n",
"# seed for reproducibility\n",
"np.random.seed(42)\n",
"\n",
"X = []\n",
"y = []\n",
"mu = 0.1\n",
"L = 1\n",
"gamma_list = np.linspace(0.01, 1.99, 20)\n",
"\n",
"for gamma in gamma_list:\n",
" pepit_tau, list_of_constraints = wc_gradient_descent_function_values_sparse_proof(mu, L, gamma, verbose=0)\n",
" l1 = list_of_constraints[2]._dual_variable_value\n",
" \n",
" X.append([np.sqrt(pepit_tau)])\n",
" y.append(l1)\n",
"\n",
"model = PySRRegressor(\n",
" maxsize=10,\n",
" niterations=10,\n",
" binary_operators=[\"+\", \"-\", \"*\"],\n",
" verbosity=0,\n",
" progress=False,\n",
" deterministic=True, # Just to maintain consistency in outputs\n",
" parallelism='serial',\n",
")\n",
" \n",
"result = model.fit(np.array(X), np.array(y), variable_names=[\"rho\"])"
]
},
{
"cell_type": "code",
"execution_count": 65,
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(1.0000063 - rho) * rho\n"
]
}
],
"source": [
"print(model.get_best().equation)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This recovers the expression $\\lambda_1(\\gamma)= (1-\\rho(\\gamma))\\rho(\\gamma)$. This was a simple, univariate function, which is where symbolic regression works best.\n",
"\n",
"The next cell uses the flexibility of the PySR API in order to learn the convergence rate of gradient descent."
]
},
{
"cell_type": "code",
"execution_count": 62,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"import itertools\n",
"np.random.seed(42)\n",
"\n",
"X = []\n",
"y = []\n",
"\n",
"# Use different values for L and mu to generate dataset\n",
"L_values = [1, 2, 10]\n",
"mu_values = [0.1, 0.9]\n",
"\n",
"for L_val in L_values:\n",
" for mu_val in mu_values:\n",
" \n",
" # Generate uniform points between 0 and 2/L (at which point we do not have convergence)\n",
" limit = 2.0 / L_val\n",
" gammas = np.linspace(0.01, limit - 0.01, 5)\n",
" \n",
" for g_val in gammas:\n",
" pepit_tau, _ = wc_gradient_descent_function_values_sparse_proof(mu_val, L_val, g_val, verbose=0)\n",
" X.append([mu_val, L_val, g_val])\n",
" y.append(pepit_tau)\n",
"\n",
"# Increased maxsize to allow for the expression complexity\n",
"model = PySRRegressor(\n",
" niterations=200,\n",
" binary_operators=[\"+\", \"-\", \"*\", \"max\"], # PySR supports many operators, such as max, min, abs, etc\n",
" unary_operators=[\"square\"],\n",
" maxsize=15,\n",
" verbosity=0,\n",
" progress=False,\n",
" deterministic=True,\n",
" parallelism='serial',\n",
")\n",
"\n",
"result = model.fit(np.array(X), np.array(y), variable_names=[\"mu\", \"L\", \"g\"])"
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"square(max(1.0 - (g * mu), (L * g) - 1.0))\n"
]
}
],
"source": [
"print(model.get_best().equation)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"PySR was able to effectively learn this convergence rate, which is a function of 3 different variables. This did, however, require guiding it to use the max and square operators. In other problems, a larger number of operators should be used, which does decrease the speed of convergence of the heuristic. This is why this type of approach works best for problems with relatively simple closed forms. When it works, it can save a lot of time."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"