Working on some stuff

Author: gnthibault

.gitignore

@@ -1,7 +1,10 @@
-.ipynb_checkpoints
-*.swp
 *.dat
-.map
 *.dot
 *.dot.pdf
+*.fits
+*.ipynb_checkpoints
+*.jpg
+*.map
+*.swp
+*.tar.gz
 

DifferentiatingPerceptron.ipynb

@@ -267,12 +267,14 @@
    "source": [
     "### Quadratic form\n",
     "\n",
-    "The matrix definition of an ellipsoid : $x^{\\intercal}Mx$ can obviously be developped, and give a polynomial of degree 2 of the form\n",
+    "The matrix definition of an ellipsoid : $x^{\\intercal}Mx$ can obviously be developped, and give a multivariate polynomial of degree 2 of the form\n",
     "$$\n",
     "\\sum_{i=0, j=0}^{n-1, n-1} M_{ij} x_i x_j\n",
     "$$\n",
     "\n",
-    "This sum has $n^2$ terms, but, as all $M_{ij}$ and $M_{ji}$ terms can be factorized, the underlying quadratic form actually has $\\frac{n \\times (n-1)}{2}$ terms. Unfortunately, we won't talk much about quadratic forms here."
+    "This sum has $n^2$ terms, but, as all $M_{ij}$ and $M_{ji}$ terms can be factorized, the underlying quadratic form actually has $\\frac{n \\times (n-1)}{2}$ terms. Unfortunately, we won't talk much about quadratic forms here.\n",
+    "\n",
+    "More on quadratic forms in the notebook QuadraticFormsForFun"
    ]
   },
   {
@@ -457,6 +459,132 @@
     "plt.axis('equal')"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Ellipsoid estimation from boundary points\n",
+    "\n",
+    "In the field of image  processing, it might sometimes be needed to estimate an ellipse, not from a a set of points pertaining to that geometric objects, but directly to a set of points that are supposed to lie on the edge of the ellipsoid.\n",
+    "\n",
+    "### Origin centered ellipsoid\n",
+    "\n",
+    "Lets recall that a origin centered ellipsoid E equation can be written as:\n",
+    "\n",
+    "\\begin{align*}\n",
+    "    E &= x \\quad \\text{s.t.} \\sqrt{x^t M x} = r \\\\\n",
+    "    &= x \\quad \\text{s.t.} \\sqrt{x^t M x} - r = 0 \\\\\n",
+    "\\end{align*}\n",
+    "\n",
+    "In 2d, given we can also write this more explicitly as:\n",
+    "\n",
+    "\\begin{align*}\n",
+    "    (x_0, x_1) \\quad \\text{s.t.} \\sqrt{\n",
+    "    \\begin{pmatrix}\n",
+    "        x_0\\\\\n",
+    "        x_1\n",
+    "    \\end{pmatrix} ^t\n",
+    "    \\begin{pmatrix}\n",
+    "        M_a & M_b\\\\\n",
+    "        M_c & M_d\\\\\n",
+    "    \\end{pmatrix}\n",
+    "    \\begin{pmatrix}\n",
+    "        x_0\\\\\n",
+    "        x_1\n",
+    "    \\end{pmatrix}} - r = 0 \\\\\n",
+    "    \\sqrt{\n",
+    "    \\begin{pmatrix}\n",
+    "        x_0\\\\\n",
+    "        x_1\n",
+    "    \\end{pmatrix} ^t\n",
+    "    \\begin{pmatrix}\n",
+    "        M_a x_0 + M_b x_1\\\\\n",
+    "        M_c x_0 + M_d x_1\\\\\n",
+    "    \\end{pmatrix}} - r = 0 \\\\\n",
+    "    \\sqrt{M_a x_0^2 + M_d x_1^2 + M_b M_c x_0 x_1} - r = 0\n",
+    "\\end{align*}\n",
+    "\n",
+    "Lets take a closer look at the resulting multivariate polynomial:\n",
+    "\\begin{align*}\n",
+    "    \\sqrt{M_a x_0^2 + M_d x_1^2 + M_b M_c x_0 x_1} - r = 0 \\\\\n",
+    "    \\Leftrightarrow M_a x_0^2 + M_d x_1^2 + M_b M_c x_0 x_1 - r^2 = 0\n",
+    "\\end{align*}\n",
+    "\n",
+    "Given multiple instances of $(x_{0,i}, x_{1,i})$ pairs, we can easily setup a least square problem, in order to find an optimally fitting ellipsoid, in the least square sense:\n",
+    "\n",
+    "\\begin{align*}\n",
+    "    \\underset{{M_a, M_b, M_c, M_d, r} \\in \\mathbb{R}^5}{argmin} &\\|\\begin{pmatrix}\n",
+    "        x_{0,0}^2 & x_{1,0}^2 & x_{0,0} x_{1,0} & -1\\\\\n",
+    "        x_{0,0}^2 & x_{1,0}^2 & x_{0,0} x_{1,0} & -1\\\\\n",
+    "        \\vdots    & \\vdots    & \\vdots          & -1\\\\\n",
+    "        x_{0,n-1}^2 & x_{1,n-1}^2 & x_{0,n-1} x_{1,n-1} & -1\\\\\n",
+    "    \\end{pmatrix}\n",
+    "    \\begin{pmatrix}\n",
+    "        M_a \\\\ M_d \\\\ M_b M_c \\\\ r^2\n",
+    "    \\end{pmatrix} - \\vec{0}\\|_2^2 \\\\\n",
+    "    \\Leftrightarrow\n",
+    "    \\underset{{M_a, M_b, M_c, M_d, r} \\in \\mathbb{R}^5}{argmin} &\\|\\begin{pmatrix}\n",
+    "        x_{0,0}^2 & x_{1,0}^2 & x_{0,0} x_{1,0} & -1\\\\\n",
+    "        x_{0,0}^2 & x_{1,0}^2 & x_{0,0} x_{1,0} & -1\\\\\n",
+    "        \\vdots    & \\vdots    & \\vdots          & -1\\\\\n",
+    "        x_{0,n-1}^2 & x_{1,n-1}^2 & x_{0,n-1} x_{1,n-1} & -1\\\\\n",
+    "    \\end{pmatrix}\n",
+    "    \\begin{pmatrix}\n",
+    "        a \\\\ b \\\\ c \\\\ d\n",
+    "    \\end{pmatrix}\\|_2^2 \\\\\n",
+    "    \\Leftrightarrow\n",
+    "    \\underset{{M_a, M_b, M_c, M_d, r} \\in \\mathbb{R}^5}{argmin} &M_2\n",
+    "    \\begin{pmatrix}\n",
+    "        a \\\\ b \\\\ c \\\\ d\n",
+    "    \\end{pmatrix}\\|_2^2 \n",
+    "\\end{align*}\n",
+    "\n",
+    "On can easily see that this problem is equivalent to finding the kernel of matrix $M_2$ or $M_2^t M_2$.\n",
+    "\n",
+    "Of course, the zero vector is a trivial solution. This is why it is usually better to formulate the problem with the builtin constraint that the zero solution is excluded (that project the problem from $\\mathbb{R}^5$ onto the lower dimensional $\\mathbb{R}^{4*}$ space):\n",
+    "\n",
+    "\\begin{align*}\n",
+    "    \\sqrt{M_a x_0^2 + M_d x_1^2 + M_b M_c x_0 x_1} &= r \\\\\n",
+    "    \\Leftrightarrow M_a x_0^2 + M_d x_1^2 + M_b M_c x_0 x_1 - r^2 &= 0\\\\\n",
+    "    \\Leftrightarrow \\frac{M_a}{r^2} x_0^2 + \\frac{M_d}{r^2} x_1^2 + \\frac{M_b M_c}{r^2} x_0 x_1 - 1 &= 0  \\quad \\text{with } r\\neq 0\\\\\n",
+    "\\end{align*}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Arbitrary ellipsoid\n",
+    "\n",
+    "Lets now assume that the unknown ellipsoid can be arbitrary centered at $$c=\\begin{pmatrix}c_0\\\\c_1end{pmatrix}$$\n",
+    "\n",
+    "\\begin{align*}\n",
+    "    (x_0, x_1-c_1) \\quad \\text{s.t.} \\sqrt{\n",
+    "    \\begin{pmatrix}\n",
+    "        x_0 - c_0\\\\\n",
+    "        x_1 - c_1\n",
+    "    \\end{pmatrix} ^t\n",
+    "    \\begin{pmatrix}\n",
+    "        M_a & M_b\\\\\n",
+    "        M_c & M_d\\\\\n",
+    "    \\end{pmatrix}\n",
+    "    \\begin{pmatrix}\n",
+    "        x_0 - c_0\\\\\n",
+    "        x_1 - c_1\n",
+    "    \\end{pmatrix}} - r = 0 \\\\\n",
+    "    \\sqrt{\n",
+    "    \\begin{pmatrix}\n",
+    "        x_0\\\\\n",
+    "        x_1\n",
+    "    \\end{pmatrix} ^t\n",
+    "    \\begin{pmatrix}\n",
+    "        M_a x_0 + M_b x_1\\\\\n",
+    "        M_c x_0 + M_d x_1\\\\\n",
+    "    \\end{pmatrix}} - r = 0 \\\\\n",
+    "    \\left(M_a \\right) x_0^2 + \\left(M_a M_b \\right) x_1^2 + \\left(M_a M_b \\right) x_0 + \\left(M_a M_b \\right) x_1 + \\left(M_a M_b \\right) x_0 x_1 + \\left(M_a M_b \\right)\n",
+    "\\end{align*}\n"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -963,7 +1091,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.6"
+   "version": "3.6.7"
   }
  },
  "nbformat": 4,

EllipsoidEstimation.ipynb

@@ -649,7 +649,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.3"
+   "version": "3.6.7"
   }
  },
  "nbformat": 4,

ForwardBackwardDual.ipynb

@@ -878,7 +878,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.6"
+   "version": "3.6.7"
   }
  },
  "nbformat": 4,

ForwardBackwardSparseL1.ipynb

@@ -1329,7 +1329,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.3"
+   "version": "3.6.7"
   }
  },
  "nbformat": 4,

GradientDescentDifferentiable.ipynb

@@ -1135,7 +1135,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.6"
+   "version": "3.6.7"
   }
  },
  "nbformat": 4,

HilbertTransform.ipynb

@@ -652,7 +652,7 @@
    "source": [
     "### Recall on Strong duality\n",
     "\n",
-    "We recall that strong duality holds for non constrained quadratic programs like the on we are studying, such that the optimal primal and dual costs satisfy\n",
+    "We recall that strong duality holds for non constrained quadratic programs like the one we are studying, such that the optimal primal and dual costs satisfy\n",
     "\n",
     "\\begin{align*}\n",
     "E^\\star= F(x^\\star) + G(Lx^\\star) = -  F^*( -L^* u^\\star ) - G^*(u^\\star)\n",
@@ -707,24 +707,24 @@
     "\n",
     "\\begin{align*}\n",
     "    c(z) &= \\quad \\langle u,z \\rangle - \\frac{1}{2}\\|Az-y\\|_2^2 \\\\\n",
-    "    &= \\quad \\langle u,z \\rangle - \\frac{1}{2} \\left( \\| Az,Az \\|_2^2 +\n",
-    "    \\|y\\|_2^2 + \\langle Az,y rangle \\right)\n",
+    "    &= \\quad \\langle u,z \\rangle - \\frac{1}{2} \\left( \\| Az \\|_2^2 +\n",
+    "    \\|y\\|_2^2 - 2 \\langle Az,y \\rangle \\right)\n",
     "\\end{align*}\n",
     "\n",
     "This nice concave function can be derived in order to find its critical point:\n",
     "\n",
     "\\begin{align*}\n",
     "    \\frac{\\partial c}{\\partial z} &= 0 \\\\\n",
-    "    u - A^t A z - A^t y &= 0 \\\\\n",
-    "    A^t A z &= u - A^t y \\\\\n",
-    "    z &= (A^t A)^{-1} \\left( u - A^t y \\right)\n",
+    "    u - A^t A z + A^t y &= 0 \\\\\n",
+    "    A^t A z &= u + A^t y \\\\\n",
+    "    z &= (A^t A)^{-1} \\left( u + A^t y \\right)\n",
     "\\end{align*}\n",
     "\n",
     "Reinjecting this expression in $ F^*(u)$ gives:\n",
     "\n",
     "\\begin{align*}\n",
     "    F^*(u) &= \\langle u,(A^t A)^{-1} \\left( u - A^t y \\right) \\rangle -\n",
-    "    \\frac{1}{2}\\|A \\left(A^t A\\right)^{-1} \\left( u - A^t y \\right) -y\\|_2^2\n",
+    "    \\frac{1}{2}\\|A \\left(A^t A\\right)^{-1} \\left( u + A^t y \\right) -y\\|_2^2\n",
     "\\end{align*}"
    ]
   },
@@ -740,14 +740,14 @@
     "Lets use these shorthand to express $F^*$:\n",
     "\n",
     "\\begin{align*}\n",
-    "    F^*(u) &= \\langle u,P ( u - A^t y ) \\rangle - \\frac{1}{2}\\|A P ( u - A^t y ) -y\\|_2^2 \\\\\n",
-    "    &= \\langle u,Pu \\rangle - \\langle u,\\lambda \\rangle - \\frac{1}{2} \\| APu - A\\lambda - y\\|_2^2 \\\\\n",
+    "    F^*(u) &= \\langle u,P ( u - A^t y ) \\rangle - \\frac{1}{2}\\|A P ( u + A^t y ) -y\\|_2^2 \\\\\n",
+    "    &= \\langle u,Pu \\rangle - \\langle u,\\lambda \\rangle - \\frac{1}{2} \\| APu + A\\lambda - y\\|_2^2 \\\\\n",
     "\\end{align*}\n",
     "\n",
     "So that the derivative gives:\n",
     "\n",
     "\\begin{align*}\n",
-    "    \\frac{\\partial F^*}{\\partial u} &= - \\lambda - \\left( P^t A^t APu - P^t A^t (A\\lambda+y) \\right) \\\\\n",
+    "    \\frac{\\partial F^*}{\\partial u} &= Pu - \\lambda - \\left( P^t A^t APu - P^t A^t (A\\lambda-y) \\right) \\\\\n",
     "\\end{align*}"
    ]
   },
@@ -757,7 +757,7 @@
    "source": [
     "### Convex conjugate of $g$\n",
     "\n",
-    "As seen previously, $g$ is the trivial function $g(x)=0$, its proximity operator is the identity: $\\mathrm{prox}_{\\gamma g}(x) = x$. Its convex conjugate reads $G^*(u)= \\underset{z}{max} \\quad \\langle u, z \\rangle_{\\mathbb{R}^n}$ which reduces to the constraint $u=0$ that translate into the indicator function of 0 : $\\delta_0(u)$ \\\\\n",
+    "As seen previously, $g$ is the trivial function $g(x)=0$, its proximity operator is the identity: $\\mathrm{prox}_{\\gamma g}(x) = x$. Its convex conjugate reads $G^*(u)= \\underset{z}{max} \\quad \\langle u, z \\rangle_{\\mathbb{R}^n}$ which reduces to the constraint $u=0$ that translate into the indicator function of 0 : $\\delta_0(u)$\n",
     "\n",
     "Thanks to Moreau identity, we can esily derive the proximity operator of $G^*$:\n",
     "\n",
@@ -787,7 +787,7 @@
     "\\begin{align*}\n",
     "    & F^*( -L^* u ) + G^*(u)\\\\\n",
     "    =& F^*( -u ) + \\delta_0(u)\\\\\n",
-    "    =& \\langle u,Pu \\rangle + \\langle u,\\lambda \\rangle - \\frac{1}{2} \\| APu + A\\lambda + y\\|_2^2 + \\delta_0(u)\n",
+    "    =& \\langle u,Pu \\rangle + \\langle u,\\lambda \\rangle - \\frac{1}{2} \\| APu + A\\lambda - y\\|_2^2 + \\delta_0(u)\n",
     "\\end{align*}\n",
     "\n",
     "With $P = \\left(A^t A\\right)^{-1}$ and $\\lambda = \\left(A^t A\\right)^{-1}A^t y$\n",
@@ -796,7 +796,7 @@
     "\n",
     "\\begin{align*}\n",
     "    x^{k} &= \\nabla f^*( -L^* u^{k} )\\\\\n",
-    "    x^{k} &= - \\lambda - \\left( -P^t A^t APu^{k} - P^t A^t (A\\lambda+y) \\right)\n",
+    "    x^{k} &= Pu^{k} - \\lambda - \\left( -P^t A^t APu^{k} - P^t A^t (A\\lambda-y) \\right)\n",
     "\\end{align*}\n",
     "\n",
     "And the dual update\n",
@@ -809,14 +809,14 @@
     "\n",
     "The dual objective now reads\n",
     "\\begin{align*}\n",
-    "    &= \\langle u,Pu \\rangle + \\langle u,\\lambda \\rangle - \\frac{1}{2} \\| APu + A\\lambda + y\\|_2^2 + \\delta_0(u) \\\\\n",
-    "    &= - \\frac{1}{2} \\| A \\lambda + y\\|_2^2 \\\\\n",
-    "    &= - \\frac{1}{2} \\| A \\left(A^t A\\right)^{-1}A^t y + y\\|_2^2\n",
+    "    &= \\langle u,Pu \\rangle + \\langle u,\\lambda \\rangle - \\frac{1}{2} \\| APu + A\\lambda - y\\|_2^2 + \\delta_0(u) \\\\\n",
+    "    &= - \\frac{1}{2} \\| A \\lambda - y\\|_2^2 \\\\\n",
+    "    &= - \\frac{1}{2} \\| A \\left(A^t A\\right)^{-1}A^t y - y\\|_2^2\n",
     "\\end{align*}\n",
     "\n",
     "and the forward backward algorithm reads\n",
     "\\begin{align*}\n",
-    "    x^{k} &= - \\lambda + P^t A^t (A\\lambda+y)\\\\\n",
+    "    x^{k} &= - \\lambda + P^t A^t (A\\lambda-y)\\\\\n",
     "    u^{k+1} &= 0\n",
     "\\end{align*}"
    ]
@@ -847,7 +847,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.3"
+   "version": "3.6.7"
   }
  },
  "nbformat": 4,

LeastSquareWithOptimalityCertificate.ipynb

@@ -0,0 +1,45 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Numerical stuff\n",
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Quadratic forms for fun\n",
+    "\n",
+    "## Introduction\n",
+    "This notebook is a kind of follow up on the EllipsoidEstimation notebook."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

QuadraticFormsForFun.ipynb

@@ -6023,7 +6023,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.6"
+   "version": "3.6.7"
   }
  },
  "nbformat": 4,