SparseFourierApproximation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.linalg as scl\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sparse Fourier approximation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Discrete Fourier transform as a linear operator\n",
    "Let's now consider two discretized signals $x_1, x_2$ that are periodic of period of period N: $x_i[n+N]=x_i[n]$.\n",
    "\n",
    "The discrete fourier transform reads:\n",
    "\\begin{align*}\n",
    "  X[k]= \\sum_{n=0}^{N-1} x[n]e^{-2 \\pi j \\frac{kn}{N}}\n",
    "\\end{align*}\n",
    "And the inverse discrete fourier transform reads:\n",
    "\\begin{align*}\n",
    "  x[n]= \\frac{1}{N} \\sum_{k=0}^{N-1} X[k]e^{2 \\pi j \\frac{kn}{N}}\n",
    "\\end{align*}\n",
    "\n",
    "it is very interesting to see that each element of the direct space is the result of a dot product with Fourier space basis, and conversely, each element of the Fourier space is also the result of a dot product with a vector of the analysis Fourier basis.\n",
    "This is much more clear if we write the discrete Fourier transform as a matrix F:\n",
    "$X = F^{-1} x$ and $x = F X$ where we have:\n",
    "\n",
    "\\begin{align}\n",
    "    F = \\frac{1}{N} \\begin{pmatrix}\n",
    "    1 & 1 & \\cdots & 1\\\\\n",
    "    1 & e^{2\\pi j\\frac{1}{N}} & \\cdots & e^{2\\pi j\\frac{N-1}{N}}\\\\\n",
    "    \\vdots & \\vdots & \\vdots & \\vdots \\\\\n",
    "    1 & e^{2\\pi j\\frac{N-1}{N}} & \\cdots & e^{2\\pi j\\frac{(N-1)(N-1)}{N}}\\\\\n",
    "    \\end{pmatrix}\n",
    "\\end{align}\n",
    "where k indices are constant across a row, but are different along a column\n",
    "and\n",
    " \n",
    "\\begin{align}\n",
    "    F^{-1} = \\begin{pmatrix}\n",
    "    1 & 1 & \\cdots & 1\\\\\n",
    "    1 & e^{-2\\pi j\\frac{1}{N}} & \\cdots & e^{-2\\pi j\\frac{N-1}{N}}\\\\\n",
    "    \\vdots & \\vdots & \\vdots & \\vdots \\\\\n",
    "    1 & e^{-2\\pi j\\frac{N-1}{N}} & \\cdots & e^{-2\\pi j\\frac{(N-1)(N-1)}{N}}\\\\\n",
    "    \\end{pmatrix}\n",
    "\\end{align}\n",
    "where n indices are constant across row, but are different along a column\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simple numerical check of the definitions above\n",
    "N = 4\n",
    "F = np.zeros((N,N), dtype=np.complex64)\n",
    "Fm = np.zeros_like(F)\n",
    "ks = np.arange(N)\n",
    "ns = np.arange(N)\n",
    "\n",
    "# definition of F\n",
    "for n in ns:\n",
    "    F[n,:] = np.exp(2*np.pi*1j*ks*n/N)\n",
    "F = F/N\n",
    "\n",
    "# definition of F^(-1)\n",
    "for k in ks:\n",
    "    Fm[k,:] = np.exp(-2*np.pi*1j*k*ns/N)\n",
    "\n",
    "# Small random test with vector, notice there is no FFTSHIFT here !\n",
    "x=np.random.rand(N)+1j*np.random.rand(N)\n",
    "r=np.dot(F,x)\n",
    "assert np.allclose(r,np.fft.ifft(x)), 'F is not numerically close to ifft'\n",
    "r=np.dot(Fm,x)\n",
    "assert np.allclose(r,np.fft.fft(x)), 'Fm is not numerically close to fft'\n",
    "# Check if they are inverse of each other\n",
    "assert np.allclose(np.dot(F,Fm), np.identity(N, dtype=np.complex64), atol=1e-5), 'F and Fm are not numerically each other inv'\n",
    "assert np.allclose(np.dot(Fm,F), np.identity(N, dtype=np.complex64), atol=1e-5), 'Fm and F are not numerically each other inv'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alternative definition of Fourier transform as a unitary matrix, and its property\n",
    "\n",
    "### Definition of $G$\n",
    "It is interesting to notice that, depending on a scaling factor used for our definition, we can write that $G = \\frac{1}{\\sqrt{N}} F$ is a unitary matrix.\n",
    "Ie the following equalities holds:\n",
    "* $G.G^{\\star} = Id(N)$\n",
    "* $G^{\\star}.G = Id(N)$\n",
    "\n",
    "We recall that $G^{\\star}$ is called the adjoint matrix of $G$, it is the transpose of G where each coefficient has been conjugated.\n",
    "This construction is also sometimes called the Hermitian conjugate.\n",
    "Showing that $\\frac{1}{\\sqrt{N}} F$ is the actual adjoint of $\\sqrt{N} F^{-1}$ is rather easy if you notice that those matrices are both symmetric, and that $\\overline{e^{\\alpha j}} = e^{-\\alpha j}$\n",
    "\n",
    "The fact that the product makes them unitary is also pretty easy when dividing the output in two case:\n",
    "* on diagonal: each term in the sum has the form $e^{\\alpha j} \\times e^{-\\alpha j} = 1$ and the $\\frac{1}{\\sqrt{N}$ makes the overall norm to be 1\n",
    "* off diagonal: We can reorder all the terms to exhibit 1 on one side minus the formula for partition of unity on the other side which amounts to zero.\n",
    "\n",
    "### Orthonogonality for complex valued vectors\n",
    "The very nice property that arise when one matrix verify the unitary property is the orthogonality of its column vector. When working on the field of real numbers, unitary matrices and orthonormal matrices have the same definition because transpose (orthonormal) is the same as hermitian transpose (unitary).\n",
    "\n",
    "On the field of complex number, we talk about unitary matrices, but the semantic behind is close: you can see directly from $G.G^{\\star} = Id(N)$ that each column vector of $G$ either has its hermitian inner product with the other that amounts to 0, or hermitian inner product with itself that amounts to 1, not simply to a complex number of norm one.\n",
    "\n",
    "Working on a field equipped with an hermitian product helps define useful algebraic structure like Hilbert spaces.\n",
    "\n",
    "### Unitary matrices are stable\n",
    "The product of two unitary matrices is another unitary matrix. The inverse of a unitary matrix is another unitary matrix, and identity matrices are unitary. Hence the set of unitary matrices form a group, called the unitary group.\n",
    "\n",
    "This property is especially usefull when composing transformations.\n",
    "\n",
    "### Unitary does not mean hermitian\n",
    "Hermitian symmetry extends some of the properties of symmetric matrices to the complex field, like positive semi-definiteness:  $<z, Mz> = z^{\\star} M z \\geq 0 \\forall z$.\n",
    "\n",
    "In the case of complex number, hermitian matrices are defined as: $A \\; \\text{s.t} \\; a_{ij} = \\overline{a_{ji}}$\n",
    "\n",
    "All hermitian matrices are positive semi definite. But there is no reason for our Fourier matrix $G$ to be hermitian\n",
    "\n",
    "\n",
    "### Spectral theorem for normal matrices\n",
    "We also recall that unitary matrices are also commuting matrices by construction: $AB = BA$ when one matrix commutes with its adjoint we also say that this matrix is normal.\n",
    "\n",
    "One can show that A is normal if and only if it is unitarily diagonalizable:\n",
    "\\begin{align}\n",
    "    &G = V D V^{\\star}\n",
    "\\end{align}\n",
    "Where $V$ is also a unitary matrix. The proof uses Schurr decomposition.\n",
    "\n",
    "Interestingly, given the unitary property in addition, both $G^{\\star}$ and $G$  have all their eigenvalues lying on the unit circle. This can be shown using the fact that unitary matrices preserve the Hermitian inner product:\n",
    "\n",
    "\\begin{align}\n",
    "    &<\\vec{v},\\vec{w}>=\\overline{<\\vec{w},\\vec{v}>} \\; \\text{and we have} \\\\\n",
    "    &<\\vec{v},\\vec{w}>=<G\\vec{v},G\\vec{w}>\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lets define G and numerically check previous properties\n",
    "F = np.zeros((N,N), dtype=np.complex64)\n",
    "for n in ns:\n",
    "    F[n,:] = np.exp(2*np.pi*1j*ks*n/N)\n",
    "G = F/np.sqrt(N)\n",
    "\n",
    "# definition of G^(-1)\n",
    "Fm = np.zeros_like(F)\n",
    "for k in ks:\n",
    "    Fm[k,:] = np.exp(-2*np.pi*1j*k*ns/N)\n",
    "Gm = Fm/np.sqrt(N)\n",
    "\n",
    "def check_unitary_matrices(A, B, atol=1e-5):\n",
    "    # Check if they are adjoint of each other.\n",
    "    assert np.allclose(A.T.conj(), B, atol=atol), 'A adjoint is not numerically close to B'\n",
    "    assert np.allclose(B.T.conj(), A, atol=atol), 'B adjoint is not numerically close to A'\n",
    "\n",
    "    # Now check invert\n",
    "    assert np.allclose(np.dot(A.T.conj(), A), np.identity(N, dtype=np.complex64), atol=atol), 'A adjoint is not numerically close to A inverse'\n",
    "    assert np.allclose(np.dot(A, A.T.conj()), np.identity(N, dtype=np.complex64), atol=atol), 'A adjoint is not numerically close to A inverse'\n",
    "    \n",
    "check_unitary_matrices(G, Gm)\n",
    "\n",
    "# This was not working IN GENERAL with original definition of F and Fm\n",
    "assert not np.allclose(np.dot(F.T.conj(), F), np.identity(N, dtype=np.complex64), atol=1e-5), 'F adjoint is numerically close to Gm inverse in general !'\n",
    "assert not np.allclose(np.dot(Fm.T.conj(), Fm), np.identity(N, dtype=np.complex64), atol=1e-5), 'Fm adjoint is not numerically close to F inverse in general !'\n",
    "\n",
    "# Check for diagonalization with unitary matrices with SVD\n",
    "u,s,vh = np.linalg.svd(G, full_matrices=False)\n",
    "Ddiag = np.diag(s)\n",
    "assert np.allclose(np.abs(Ddiag), np.identity(N, dtype=np.complex64), atol=1e-5), 'Ddiag has not all unit singular values'\n",
    "assert np.allclose(np.dot(u, np.dot(Ddiag, vh)), G, atol=1e-5), 'SVD reconstruction is not ok'\n",
    "check_unitary_matrices(u, u.T.conj())\n",
    "check_unitary_matrices(vh.T.conj(), vh)\n",
    "# Surprisingly here, SVD does not give a numerically accurate estimate: u and vh are not adjoint of each other\n",
    "assert not np.allclose(u.T.conj(), vh, atol=1e-5), 'SVD does not provide numerically close result'\n",
    "\n",
    "# Check for diagonalization with unitary matrices with ED\n",
    "w,v = np.linalg.eig(G)\n",
    "Ddiag = np.diag(w)\n",
    "#print(np.dot(v, v.T.conj()))\n",
    "assert np.allclose(np.abs(Ddiag), np.identity(N, dtype=np.complex64), atol=1e-5), 'Ddiag has not all unit singular values'\n",
    "# This is the expected reconstruction formula for eigen decomposition\n",
    "assert np.allclose(np.dot(v, np.dot(Ddiag, np.linalg.inv(v))), G, atol=1e-5), 'QDQ-1 reconstruction is not ok'\n",
    "\n",
    "#It is not clear why those two properties related to unitary diagonalization do not hold up to a good precision ?\n",
    "check_unitary_matrices(v, v.T.conj(), atol=1e-1)\n",
    "assert np.allclose(np.dot(v, np.dot(Ddiag, v.T.conj())), G, atol=1e-1), 'DQt reconstruction is not ok'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fourier matrices diagonalizes circulant matrices\n",
    "\n",
    "Although this is not relevant to our current problem, we would like to recall this important property of the discrete Fourier transform, especially for \"numerical\" people use this powerful tool.\n",
    "\n",
    "### Circular convolution as a linear operator\n",
    "We also recall that discrete circular convolution reads:\n",
    "\\begin{align*}\n",
    "  (s_1 \\star s_2)[n] = \\sum_{k=0}^{N-1} s_1[k] s_2[(n-k)\\%N]\n",
    "\\end{align*}\n",
    "\n",
    "It is interesting to notice that this convolution operator can be expressed as a linear operation, $(s_1 \\star s_2) = S_2 s_1$ where we have:\n",
    "\\begin{align}\n",
    "    S_2 = \\begin{pmatrix}\n",
    "    s_2[0] & s_2[N-1] & \\cdots & s_2[1] \\\\\n",
    "    s_2[1] & s_2[0] & \\cdots & s_2[2] \\\\\n",
    "    \\vdots & \\vdots & \\vdots & \\vdots \\\\\n",
    "    s_2[N-1] & s_2[N-2] & \\cdots & s_2[0] \\\\\n",
    "    \\end{pmatrix}\n",
    "\\end{align}\n",
    "And one can easily see that this matrix is a circulant matrix\n",
    "\n",
    "### Diagonalization of circulant matrices\n",
    "From there, one can show that the $F, F^{-1}$ pair of matrices can be used to diagonalize any circulant matrices $C$. such that we can write $C = F D F^{-1}$ or, alternatively $D = F^{-1} C F$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fafbbd23978>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAACBCAYAAADpLPAWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAJ6ElEQVR4nO3d249dZRnH8d+vezqdAjPQCyV2SgKxYAI3rZ00GsQYjYhihAsTIdFEbxpjIGhMDPwF3nlIPKVBvBHtBQIhSjgk2KgJAjO0jdKCGRoSppxK0nZGmXMfL7oJY207692z1pq37/5+kqbde9693mfl1zxd3bPXM44IAQDytWG9CwAAXBiNGgAyR6MGgMzRqAEgcwN1H9D2AUkfkjRZ97GRZLuk4xGxs46DkWs2yLVMF8zVdX/qw/bU5SMe3XHDpsqvOfTsZbXWAGlGJ9XRgOZj1nUcr7dch+vYGit8kOt7teXqoaHRwdGtdRwOPVo49oY2bBrU0slT58y19itqSZM7btg0+szD2yq/4AtbdzRQRn8bj/11HzI919FaLvqwwvjp/XUfcnJwdOvoR+7+Tt3HXZvU68da/tlqQMXzePNnv7jg13mPGgAyV6lR277F9iu2J23f23RRaAe5lolcy7Nqo7bdkfRzSV+UdL2kO21f33RhaBa5lolcy1Tlinq3pMmIOBoRC5L2Sbqt2bLQAnItE7kWqEqjHpX0+orHU93n/oftPbbHJe164+3lmspDg8i1TMm5Lk9Pt1UbelTbNxMjYm9EjEma2Hplp67DYp2Ra5lW5toZGVnvcrCKKo36mKSrVjze1n0OFzdyLRO5FqhKo35B0rW2r7E9KOkOSY81WxZaQK5lItcCrXrDS0Qs2b5L0pOSOpIeiIiXGq8MjSLXMpFrmSrdmRgRj0t6vOFa0DJyLRO5lqeJW8h16NnLkm4LP31T+q3GG/56IPk1WJtDzw4n3RbugY3Je8TSYvJr0K4f3vr75Nfc96c7016Q6y3hqWo6D24hB4DM0agBIHM0agDIHI0aADJHowaAzNGoASBzNGoAyByNGgAyR6MGgMzRqAEgczRqAMgcjRoAMtfIUKZUvQxYmr19d9L6zY8+n7wH1qanAUsRaetdyvSei0fygCVJ27/396T1kz/+RPIeyRL/qvU0YCl1j/PgihoAMkejBoDMrdqobV9l+8+2D9t+yfY9bRSGZpFrmci1TFXeo16S9P2IeNH2sKQJ209HxOGGa0OzyLVM5FqgVa+oI+LNiHix++cZSUckjTZdGJpFrmUi1zIlvUdt+2pJOyU9d46v7bE9LmnXvOZqKQ7tqJ7rbMuVYS2q5ro8Pd1yZUhVuVHbvkzSHyR9NyL+L9mI2BsRY5ImNmmoxhLRpLRcN7deH3qTkmtnZKT1+pCmUqO2vVFnQn8wIh5utiS0hVzLRK7lqfKpD0v6taQjEfGj5ktCG8i1TORapipX1DdK+oakz9o+2P31pYbrQvPItUzkWqBVP54XEX9TbzdPImPkWiZyLVMWsz56kTq7o7NlS9L65RMnktajJomzOzrDw0nrl2dmktajHqmzOwZPpN00vbDldNJ6Se38c1bTHtxCDgCZo1EDQOZo1ACQORo1AGSORg0AmaNRA0DmaNQAkDkaNQBkjkYNAJmjUQNA5mjUAJC5i3bWR6rU2R2LN48l77HxqfHk12BtUmd3LNy8K3mPwacmkl9z0YmEtS3MyOhpdkfBuKIGgMzRqAEgczRqAMhcyg+37dg+YPuPTRaEdpFrmci1LClX1PdIOtJUIVg35Fomci1I1Z9Cvk3SrZLub7YctIlcy0Su5al6Rf0TST+QdN7PzNjeY3tc0q55zdVRG5qXmOtsa4VhTZJyXZ6ebq0w9GbVRm37y5LeiYgLfpg0IvZGxJikiU0aqqs+NKS3XDe3Uxx61kuunZGRdopDz6pcUd8o6Su2X5O0T2d+DP1vG60KbSDXMpFrgVZt1BFxX0Rsi4irJd0h6ZmI+HrjlaFR5Fomci0Tn6MGgMwlzfqIiP2S9jdSCdYNuZaJXMvRN0OZUvUyYKmzZUvS+tRBUVi7XgYsdYaHk9anDorqCylDn6RWBj+dvmIxaf2GkxvTN0k97/PtXc9hAABNoVEDQOZo1ACQORo1AGSORg0AmaNRA0DmaNQAkDkaNQBkjkYNAJmjUQNA5mjUAJA5Zn3UKHV2x+ztu5P32Pzo88mvwdokz+6IHgY8uIXhFuspw9NLnd3x6td+lbzHR/d9O/k158IVNQBkjkYNAJmr+lPIr7D9kO2XbR+x/cmmC0PzyLVM5Fqequ9R/1TSExHxVduDki5psCa0h1zLRK6FWbVR275c0qclfVOSImJB0kKzZaFp5Fomci1Tlbc+rpF0XNJvbB+wfb/tS89eZHuP7XFJu+Y1V3edqF8Puc62XiSSJee6PD3depFIU6VRD0j6uKRfRsROSf+RdO/ZiyJib0SMSZrYpKF6q0QTesh1c8slogfJuXZGRlouEamqNOopSVMR8Vz38UM68xcBFzdyLRO5FmjVRh0Rb0l63fbHuk99TtLhRqtC48i1TORapqqf+rhb0oPd7yAflfSt5kpCi8i1TORamEqNOiIOShpruBa0jFzLRK7l4c5EAMgcQ5nWUS8Dlk7ftLPawkMTycdGTXoYsOSBagOCvNjQdKOmhyalzqlKraeHOVipehmwdN0DJyutO/XW0gW/zhU1AGSORg0AmaNRA0DmaNQAkDkaNQBkjkYNAJmjUQNA5mjUAJA5GjUAZI5GDQCZc0S9917anupoYHRYV9R6XHRdPlxp2cy/31SnM6j5+VO13Bx8JteN5NoQV7ztfCZOqKMBzcV7teXqoaHRwdGtdRwOZ7lklVvD3zcz+7Y6GwY1tzh9zlybmPVxfFlLOql3J896fkf394MN7Jmr+s/51LtVV25fXp47Xtu+0vFlLeqkjpNrE+dc/Xpp+5IWa8015uY0/+pRcj2j1vOer750+9Lp+fPmWvsV9Xk3svdLUkR8ppUNM9AP59wP53i2fjjnfjjHc8n1vHmPGgAyR6MGgMy19tYHAKA3XFEDQOZo1ACQORo1AGSulUZt+xbbr9ietH1vG3vmwPZrtv9h+6Dt8fWup27kWmauUn9mm3OujX8z0XZH0r8kfV7SlKQXJN0ZEYcb3TgDtl+TNBYRle9SuViQa5m5Sv2bbc65tnFFvVvSZEQcjYgFSfsk3dbCvmgWuZaLbDPTRqMelfT6isdT3ef6QUh6yvaE7T3rXUzNyLXMXKX+zTbbXJuY9YEPfCoijtn+sKSnbb8cEX9Z76KwZuRapmxzbeOK+pikq1Y83tZ9rngRcaz7+zuSHtGZ/1KWglzLzFXq02xzzrWNRv2CpGttX2N7UNIdkh5rYd91ZftS28Pv/1nSzZL+ub5V1Ypcy8xV6sNsc8+18bc+ImLJ9l2SnpTUkfRARLzU9L4ZuFLSI905wwOSfhcRT6xvSfUh1zJzlfo226xzZdYHAGSOOxMBIHM0agDIHI0aADJHowaAzNGoASBzNGoAyByNGgAy91/nYj2hm2sgVAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\"\n",
    "  Check if fourier transform allows to diagonalize a circulant matrix\n",
    "\"\"\"\n",
    "N=8\n",
    "\n",
    "#now generate the Fourier matrix\n",
    "x = np.linspace(0, N-1, N) \n",
    "y = np.linspace(0, N-1, N)\n",
    "xg, yg = np.meshgrid(x, y)\n",
    "F=(1/N)*np.exp(2*np.pi*1j*xg*yg/N)\n",
    "Fm=np.exp(-2*np.pi*1j*xg*yg/N)\n",
    "\n",
    "#check that F and F_1 are inverse of each other\n",
    "assert np.allclose(np.dot(F,Fm), np.identity(N)), 'Fm is not an inverse of F'\n",
    "\n",
    "#Lets check in practice that Fourier matrix can be used to diagonalize circulant matrices\n",
    "circulant=scl.circulant(np.random.randint(0,100,N))\n",
    "\n",
    "#Now show the diagonalized version of the circulant matrix:\n",
    "diag = np.dot(Fm,np.dot(circulant,F))\n",
    "\n",
    "# Check reconstruction\n",
    "assert np.allclose(circulant, np.dot(F, np.dot(diag,Fm))), 'Cannot reconstruct '\n",
    "diag = np.dot(Fm,np.dot(circulant,F))\n",
    "\n",
    "fig,ax = plt.subplots(1,3)\n",
    "ax[0].imshow(np.absolute(diag))\n",
    "ax[1].imshow(np.real(diag))\n",
    "ax[2].imshow(np.imag(diag))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sparse approximation in Fourier domain\n",
    "\n",
    "### Introduction:\n",
    "\n",
    "In a quite famous paper, Bruckstein, Donoho and Elad have explored the link between the fundamental uncertainty principle (see our notebook on the Heisenberg–Pauli–Weyl inequality) and the resolution of systems of linear equations, the name of the paper was: From sparse solutions of systems of equations to sparse modeling of signals and images.\n",
    "\n",
    "### Original problem\n",
    "\n",
    "A natural way to formalize the sparsity problem in Fourier domain, given a linear operator $M$ in the context of inverse problem, is to use a $L_0$ norm penalty on the Fourier coefficient.\n",
    "It is interesting to notice that the problem can either be casted as an \"Analysis\" problem:\n",
    "\n",
    "\\begin{align}\n",
    " \\underset{x\\in \\mathbb{R}^N}{\\text{argmin}} \\qquad \\| Mx - y \\|_2^2 + \\lambda \\|G^{-1}x\\|_0\n",
    "\\end{align}\n",
    "\n",
    "or a \"Synthesis\" problem:\n",
    "\\begin{align*}\n",
    " \\underset{x\\in \\mathbb{C}^N \\; s.t \\; Gx \\in \\mathbb{R}^N}{\\text{argmin}} \\qquad & \\| MGx - y \\|_2^2 + \\lambda \\|x\\|_0 \\\\\n",
    " \\underset{x\\in \\mathbb{C}^N}{\\text{argmin}} \\qquad & \\| MGx - y \\|_2^2  + \\lambda \\|x\\|_0 + \\chi_{\\mathbb{R}^N}(Gx)\n",
    "\\end{align*}\n",
    "\n",
    "With $\\chi_{\\mathbb{R}^N}(x)$ the convex indicator function of the set $\\mathbb{R}^N$, ie:\n",
    "\n",
    "\\begin{align*}\n",
    "  \\chi_{\\mathbb{R}^{N}}(x) &= \n",
    "    \\begin{cases}0 \\qquad &\\text{if } x \\in \\mathbb{R}^{N} \\\\\n",
    "    +\\infty & \\text{otherwise} \\end{cases}\n",
    "\\end{align*}\n",
    "\n",
    "It is not even easy to see here if the analysis problem make sense. At first sight, it seems more simple as it allows to restricts the search space to $\\mathbb{R}^N$ instead of $\\mathbb{C}^N$ and has less constraint.\n",
    "\n",
    "On the other side, it is not sure that this formulation will still make sens when we will try to solve the problem: we may have to search inside $\\mathbb{C}^N$ and eventually we also might need to use the real signal set indicator: $\\chi_{\\mathbb{R}^N}(Gx)$.\n",
    "\n",
    "For people familiar to signal processing, it is well know that real signal have a specific property in Fourier domain:\n",
    "* Real part of Fourier transform is even\n",
    "* Imaginary part of Fourier transform is odd\n",
    "\n",
    "In practice, one can see that the underlying dimension of the problem given the realness of the function in direct domain is the same wether we are trying to solve it in the direct domain or in the Fourier domain."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $L_1$ relaxation\n",
    "\n",
    "### Introduction\n",
    "We won't recall the whole history of sparsity inducing norms, that dates back to the study of variable selection by Robert Tibshirani & Al with the Lasso and elastic net algorithms. Later on, Emmanuel Candes, Donoho, Elad, Terry Tao provided more in depth understanding on the role of sparsity inducing norms, and how they are linked to linear sampling matrices.\n",
    "\n",
    "Let's just assume here that $\\|L_0\\|$ norm based functional are generally non convex, and thus cannot be optimized easily. Instead, we will replace $\\|L_0\\|$ norm by $\\|L_1\\|$ norm on the vectors coefficient, that favors sparse solutions for the whole system.\n",
    "\n",
    "The two functionals now reads:\n",
    "\n",
    "For the analysis problem:\n",
    "\\begin{align}\n",
    " \\underset{x\\in \\mathbb{R}^N}{\\text{argmin}} \\qquad \\| Mx - y \\|_2^2 + \\lambda \\|G^{-1}x\\|_1\n",
    "\\end{align}\n",
    "\n",
    "and for the synthesis problem:\n",
    "\\begin{align*}\n",
    " \\underset{x\\in \\mathbb{C}^N}{\\text{argmin}} \\qquad & \\| MGx - y \\|_2^2  + \\lambda \\|x\\|_0 + \\chi_{\\mathbb{R}^N}(Gx)\n",
    "\\end{align*}\n",
    "\n",
    "### Convexity and differentiability, when the complex plane hits you hard.\n",
    "\n",
    "$L_1$ norm is usually always associated with real coefficients,\n",
    "\n",
    "We recall the definition of the $L_1$ norm:\n",
    "\\begin{align*}\n",
    " \\|x\\|_1 = \\sum_{k=0}^{N-1} x_k\n",
    "\\end{align*}\n",
    "It simply means the norm of the absolute values, and is valid for both real and complex valued vector as long as you use the proper inner product definition, which is, for the complex case:\n",
    "\\begin{align*}\n",
    " |x| &= \\left( <x,x> \\right)^{\\frac{1}{2}} \\\\\n",
    " &= \\left( x^{\\star}x \\right)^{\\frac{1}{2}} \\\\\n",
    " &= \\left( \\overline{x^T}x \\right)^{\\frac{1}{2}}\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $L_1$ norm for complex valued vectors and differentiability\n",
    "\n",
    "Unfortunately, the one norm of a complex vector (the sum of the absolute values of the entries of the vector) is not a differentiable function. In fact, the absolute value of a scalar complex number z=x+i*y is not a differentiable function.\n",
    "\n",
    "To see this, use the Cauchy-Riemann conditions. Write the absolute value as\n",
    "\n",
    "abs(z)=abs(x+i*y)=sqrt(x^2+y^2)+i*0\n",
    "\n",
    "let u(x,y)=sqrt(x^2+y^2) (the real part) and v(x,y)=0 (the imaginary part.)\n",
    "\n",
    "If the absolute value was differentiable then it would satisfy the Cauchy-Riemann conditions. In particular, the partial derivative of u with respect to x would have to equal the partial derivative of v with respect to y. Since this clearly doesn't hold, there's no need to check the other half of the CR conditions,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}