SignalProcessingGraph.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "# Generic stuff\n",
    "from IPython.display import IFrame\n",
    "\n",
    "#Documentation : http://docs.scipy.org/doc/scipy/reference/sparse.html\n",
    "\"\"\"\n",
    "Sparse matrices can be used in arithmetic operations: they support addition, subtraction, multiplication, division, and matrix power.\n",
    "Advantages of the CSR format\n",
    "\n",
    "        efficient arithmetic operations CSR + CSR, CSR * CSR, etc.\n",
    "        efficient row slicing\n",
    "        fast matrix vector products\n",
    "Disadvantages of the CSR format\n",
    "\n",
    "        slow column slicing operations (consider CSC)\n",
    "        changes to the sparsity structure are expensive (consider LIL or DOK)\n",
    "\"\"\"\n",
    "import scipy.sparse as sparse\n",
    "import scipy.sparse.csgraph as csgraph\n",
    "from scipy.sparse import csr_matrix\n",
    "from scipy.sparse import lil_matrix\n",
    "from scipy.sparse import dia_matrix\n",
    "import scipy.sparse.linalg as spl\n",
    "\n",
    "import scipy as sc\n",
    "import scipy.misc\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "%matplotlib inline\n",
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Signal processing on graph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction\n",
    "\n",
    "First, we would like to thank John D Cook for its [short introduction](https://www.johndcook.com/blog/2016/02/09/fourier-transform-of-a-function-on-a-graph/) on Fourier series on graph, and of course [Pierre Vandergheynst](https://scholar.google.ch/citations?user=1p9NOFEAAAAJ&hl=fr), for its interesting courses and its contributions to this field.\n",
    "\n",
    "Also, you'll find interesting references in this book : Lectures on Spectral Graph Theory by FanR.K.Chung: http://www.math.ucsd.edu/~fan/cbms.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#IFrame(\"doc/SpectralGraphTheory/cbms.pdf\", width=1200, height=800)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are two very good references from Petar Veličković :\n",
    "* [Theoretical Foundations of Graph Neural Networks](https://www.youtube.com/watch?v=uF53xsT7mjc&ab_channel=PetarVeli%C4%8Dkovi%C4%87)\n",
    "* [Intro to graph neural networks (ML Tech Talks](https://www.youtube.com/watch?v=8owQBFAHw7E&ab_channel=TensorFlow)\n",
    "* [Recent set of references](https://twitter.com/PetarV_93/status/1306689702020382720)\n",
    "* [Some introduction to Graph attention network](https://petar-v.com/GAT/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fourier series\n",
    "\n",
    "Fourier series properties can be analyzed within various mathenatical framework. One of them is related to the Laplacian operator.\n",
    "\n",
    "We recall that the Laplacian operator, sometimes written $\\nabla \\cdot \\nabla$, $\\nabla^2$ or $\\Delta$, where $\\nabla$ can be written $\\left( \\frac{\\partial}{\\partial x_0}, \\frac{\\partial}{\\partial x_1}, \\dots \\frac{\\partial}{\\partial x_{n-1}} \\right)$ and the laplacian operator applied to a function $f$ reads $\\Delta f = \\sum_{i=0}^{n-1} \\frac{\\partial f}{\\partial x_i^2}$.\n",
    "\n",
    "It is interesting to notice that the trigonometric polynomial that defines the (separable) Fourier series basis elements are eigenfunctions for the laplacian operators over euclidean spaces:\n",
    "\n",
    "lets define $\\delta_{nF}(x) = e^{2 \\pi j nFx}$\n",
    "\n",
    "And\n",
    "\\begin{align}\n",
    "  \\frac{\\partial \\delta_{nF}}{\\partial x} (x) &= \\frac{\\partial  e^{2 \\pi j nFx}}{\\partial x} \\\\\n",
    "  &= 2 \\pi j nF e^{2 \\pi j nFx} \\\\\n",
    "\\end{align}\n",
    "\n",
    "And\n",
    "\\begin{align}\n",
    "  \\frac{\\partial \\delta_{nF}}{\\partial x^2} (x) &=  -(2\\pi nF)^2 e^{2 \\pi j nFx} \\\\\n",
    "  &= -(2\\pi nF)^2 \\delta_{nF}\\\\\n",
    "  &= \\lambda \\delta_{nF}\n",
    "\\end{align}\n",
    "\n",
    "Moreover, eigenfunctions of the laplacian can also be defined for functions on finite domains (Bessel functions for radially symmetric domain like disks and sphere)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Elements of calculus on graphs\n",
    "\n",
    "### Definitions, operators\n",
    "Graphs are not euclidean spaces, lets define some calculus elements before jumping on more advanced stuff:\n",
    "\n",
    "| Name | ........................ Definition ........................ | Comment |\n",
    "|:----:|:----------:|:-------:|\n",
    "| Graph | $\\mathcal{G} = (\\mathcal{V}, \\mathcal{E})$ | |\n",
    "| Vertices | $\\mathcal{V} = {1,\\dots, n}$ | n is the cardinality of the vertices |\n",
    "| Vertex weights | $b_i > 0 \\forall i \\in \\mathcal{V}$| |\n",
    "| Edges | $\\mathcal{E} \\subseteq \\mathcal{V} \\times \\mathcal{V}$ | |\n",
    "| Edge weights | $a_{ij} \\geq 0 \\quad \\forall \\quad (i,j) \\quad \\in \\mathcal{E}$ | |\n",
    "| Vertex field | $L^2(\\mathcal{V}) = {f : \\mathcal{V} \\rightarrow \\mathbb{R}^n}$ | Functional space of all $L^2$-integrable functions <br> sampled over the graph (1 real value per vertex) |\n",
    "| graph function | $f = (f_1, \\dots, f_n)$ | One element of $L^2(\\mathcal{V})$ presented above  |\n",
    "| Inner product | $\\langle f,g \\rangle_{L^2(\\mathcal{V})} = \\sum_{i \\in \\mathcal{V}} b_i f_i g_i$ | Allows to define a functional Hilbert <br> space of graph-based functions |\n",
    "| Gradient operator | $(\\nabla f)_{ij} = \\sqrt{a_{ij}}(f_i-f_j)$ | $\\nabla : L^2(\\mathcal{V}) \\rightarrow L^2(\\mathcal{E})$ |\n",
    "| Divergence operator | $(div F)_{i} = \\frac{1}{b_i} \\sum_{j:(i,j) \\in \\mathcal{E}} \\sqrt{a_{ij}}(F_{ji}-F_{ij})$ |$div : L^2(\\mathcal{E}) \\rightarrow L^2(\\mathcal{V})$ |\n",
    "| Gradient Adjoint |$\\nabla^{\\star}F = -div F$ | $\\langle F,\\nabla f \\rangle_{L^2(\\mathcal{E})} = \\langle \\nabla^{\\star}F, f \\rangle_{L^2(\\mathcal{V})} = \\langle -div F, f \\rangle_{L^2(\\mathcal{V})}$ |\n",
    "| Weight matrix | $A = (a_{ij})$ | | |\n",
    "| Degree matrix | $D = diag\\left(d_i\\right)$ where $d_i=\\sum_{j\\neq i}a_{ij})$ | | |\n",
    "| (Unnormalized) Laplacian operator | $(\\Delta F)_{i} = \\frac{1}{b_i} \\sum_{j:(i,j) \\in \\mathcal{E}} a_{ij}(f_i-f_j)$ | $\\Delta : L^2(\\mathcal{V}) \\rightarrow L^2(\\mathcal{V})$, this one is $\\Delta = D - A$ |\n",
    "| Normalized Laplacian | $(\\Delta F)_{i} = \\sum_{j:(i,j) \\in \\mathcal{E}} \\frac{1}{b_i\\sqrt{d_i d_j}} a_{ij}(f_i-f_j)$ | $\\Delta : L^2(\\mathcal{V}) \\rightarrow L^2(\\mathcal{V})$, this one is $\\Delta = I - D^{-\\frac{1}{2}} A D^{-\\frac{1}{2}}$ |\n",
    "| Random walk Laplacian | $\\Delta = I - D^{-1} A $ |  $\\Delta : L^2(\\mathcal{V}) \\rightarrow L^2(\\mathcal{V})$, this one is $\\Delta = I - D^{-1} A $ |\n",
    "\n",
    "Those definitions are coming from slides from Xavier Bresson."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Some precision about finite difference operators\n",
    "\n",
    "#### Gradient\n",
    "It must be noticed that gradient operator is a finitie difference operator of type forward difference in case we consider directed graphs.\n",
    "It is positively valued if the graph function has higher value at the terminal vertex than at the initial vertex.\n",
    "\n",
    "To get an idea of what the gradient linear operator looks like, we tried to draw it:\n",
    "\n",
    "\\begin{align*}\n",
    "    \\begin{pmatrix}\n",
    "        0 & 0 & 0 & \\cdots & 0 \\\\\n",
    "        \\sqrt{a_{0,1}} & -\\sqrt{a_{0,1}} & 0 & & 0 \\\\\n",
    "        \\sqrt{a_{0,2}} & 0 & -\\sqrt{a_{0,2}} & & 0 \\\\\n",
    "        \\vdots & & & & \\\\\n",
    "        \\sqrt{a_{0,I-1}} & 0 & 0 & \\cdots & -\\sqrt{a_{0,I-1}}\\\\\n",
    "        -\\sqrt{a_{1,0}} & \\sqrt{a_{1,0}} & 0 & \\cdots & 0 \\\\\n",
    "        0 & 0 & 0 & & 0 \\\\\n",
    "        0 & \\sqrt{a_{1,2}} & -\\sqrt{a_{1,2}} & & 0 \\\\\n",
    "        \\vdots & & & & \\\\\n",
    "        0 & \\sqrt{a_{1,I-1}} & 0 & & -\\sqrt{a_{1,I-1}} \\\\\n",
    "        -\\sqrt{a_{2,0}} & 0 & \\sqrt{a_{0,2}} & & 0 \\\\\n",
    "        0 & -\\sqrt{a_{2,1}} & \\sqrt{a_{2,1}} & & 0 \\\\\n",
    "        0 & 0 & 0 & & 0 \\\\\n",
    "        \\vdots & & & & \\\\\n",
    "        0 & 0 & \\cdots & -\\sqrt{a_{I-1,I-2}} & \\sqrt{a_{I-1,I-2}}\\\\\n",
    "        0 & 0 & 0 & & 0\n",
    "    \\end{pmatrix}\n",
    "\\end{align*}\n",
    "\n",
    "#### Divergence\n",
    "We must notice, that we prefered the following notation for the divergence operator (although it is exactly the same thing):\n",
    "\n",
    "| Name | ............................. Definition ............................. | ...............Comment............... |\n",
    "|:----:|:----------:|:-------:|\n",
    "| Divergence operator | $(div F)_{i} = \\sum_{(k,l) \\in \\mathcal{E}} \\sqrt{a_{kl}}(\\delta_{i}(l)-\\delta_{i}(k)) F_{kl}$ |$div : L^2(\\mathcal{E}) \\rightarrow L^2(\\mathcal{V})$ |\n",
    "\n",
    "Its matrix would look like that:\n",
    "\n",
    "\\begin{align*}\n",
    "    \\begin{pmatrix}\n",
    "        0 & -\\sqrt{a_{0,1}} & -\\sqrt{a_{0,2}} & \\cdots & -\\sqrt{a_{0,I-1}} & \\sqrt{a_{1,0}} & 0 &0 & \\cdots & 0 & \\sqrt{a_{2,0}} & 0 & 0 & \\cdots & 0 & 0\\\\\n",
    "        0 & \\sqrt{a_{0,1}} & 0 &  \\cdots & 0 & -\\sqrt{a_{1,0}} & 0 & -\\sqrt{a_{1,2}} & \\cdots & -\\sqrt{a_{1,I-1}} & 0 & \\sqrt{a_{2,1}} &  0 & & 0 & 0\\\\\n",
    "        0 & 0 & \\sqrt{a_{0,2}} & \\cdots & 0 & 0 & 0 & \\sqrt{a_{1,2}} & \\cdots & 0 & -\\sqrt{a_{0,2}} &  -\\sqrt{a_{2,1}} &  0 & & \\vdots &  0\\\\\n",
    "        \\vdots &  & \\cdots &  & \\vdots &  &  & \\cdots & & & & & & & \\sqrt{a_{I-1,I-2}} &  0\\\\\n",
    "        0 & 0 & 0 & \\cdots & \\sqrt{a_{0,I-1}} & 0 & 0 & 0 & \\cdots & \\sqrt{a_{1,I-1}} & 0 & 0 &  0 & \\cdots & -\\sqrt{a_{I-1,I-2}} & 0\\\\\n",
    "    \\end{pmatrix}\n",
    "\\end{align*}\n",
    "\n",
    "#### Laplacian\n",
    "We recall that laplacian operator $\\Delta$ is usually expressed as\n",
    "\\begin{align*}\n",
    "    \\Delta f &= \\nabla^2 f \\\\\n",
    "    &= \\nabla \\cdot \\nabla f \\\\\n",
    "    &= \\nabla^{T} \\nabla f \\\\\n",
    "    &= -div(grad(f))\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [],
   "source": [
    "#lets define some utilities for testing on random/symmetric sparse matrices\n",
    "\n",
    "# define graph adjacency matrix as a random csr matrix\n",
    "def gen_random_spm(size):\n",
    "    return sparse.rand(size[0], size[1], density=0.9, format='csr', dtype=np.float32, random_state=None)\n",
    "\n",
    "def gen_random_graph(size):\n",
    "    A = gen_random_spm((size,size))\n",
    "    # Now make it symmetric\n",
    "    return A+A.T-sparse.diags(A.diagonal())\n",
    "\n",
    "# define graph adjacency matrix as a small symmetric csr matrix\n",
    "def gen_small_graph():\n",
    "    weight=[4]\n",
    "    rowIdx=[0]\n",
    "    colIdx=[1]\n",
    "    A = csr_matrix((np.array(weight), (np.array(rowIdx), np.array(colIdx))),dtype=np.float32, shape=(2, 2))\n",
    "    #Build a symmetric matrix\n",
    "    return A+A.T-sparse.diags(A.diagonal())\n",
    "\n",
    "def gen_iota_graph():\n",
    "    weight=[1,2,3,4]\n",
    "    rowIdx=[0,0,1,1]\n",
    "    colIdx=[0,1,0,1]\n",
    "    return csr_matrix((np.array(weight), (np.array(rowIdx), np.array(colIdx))),dtype=np.float32, shape=(2, 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets define some operators, elementwise\n",
    "def grad(G,f):\n",
    "    \"\"\"-G is the graph adjacency matrix\n",
    "       -f is the function defined on this graph\n",
    "       Output has the dimension of the adjacency matrix\n",
    "       This operator can be considered a forward difference operation\n",
    "       Elementwise version\n",
    "    \"\"\"\n",
    "    g = lil_matrix(G.shape)\n",
    "    #iterate over all edges, row_index is i, col_idx is j\n",
    "    for row_idx, col_idx, weight in zip(*G.nonzero(), G.data):\n",
    "        g[row_idx, col_idx] = np.sqrt(weight)*(f[row_idx]-f[col_idx])\n",
    "    return g\n",
    "\n",
    "def div(G,F,b=None):\n",
    "    \"\"\"-G is the graph adjacency matrix\n",
    "       -F is the gradient of function f on this graph\n",
    "       -b is the optional vertex_weight\n",
    "       Output is 1D and has its size equal to the cardinality of the vertices\n",
    "       Elementwise version\n",
    "    \"\"\"\n",
    "    # If b is not provided, assume all vertices have weight 1\n",
    "    if b is None:\n",
    "        v_weight = lambda idx: 1\n",
    "    else:\n",
    "        v_weight = lambda idx: 1/b[idx]\n",
    "        \n",
    "    def delta_i(i):\n",
    "        return lambda x: 1 if x==i else 0\n",
    "    \n",
    "    div = np.zeros(G.shape[0])\n",
    "    for row_idx in range(div.size):\n",
    "        d = 0\n",
    "        for g_row_idx, g_col_idx, weight in zip(*G.nonzero(), G.data):\n",
    "            #print('row_idx is {}, col_idx is {}, weight is {}'.format(g_row_idx, g_col_idx, weight))\n",
    "            #print('F is {}, delta+ is {}, delta- is {}'.format(F[g_row_idx, g_col_idx],\n",
    "            #                                                   delta_i(row_idx)(g_row_idx),\n",
    "            #                                                   delta_i(row_idx)(g_col_idx)))\n",
    "            d += np.sqrt(weight)*F[g_row_idx, g_col_idx]*(delta_i(row_idx)(g_col_idx)-delta_i(row_idx)(g_row_idx))\n",
    "            #d += np.sqrt(weight)*(F[row_idx, col_idx]-F[col_idx,row_idx])\n",
    "        div[row_idx] = d#v_weight(row_idx)*d\n",
    "    return div\n",
    "\n",
    "def laplacian(G,f,b=None):\n",
    "    \"\"\"-G is the graph adjacency matrix\n",
    "       -f is a funcion defined on this graph\n",
    "       -b is the optional vertex_weight\n",
    "       Output is 1D and has its size equal to the cardinality of the vertices\n",
    "       Elementwise version\n",
    "    \"\"\"\n",
    "    lap = np.zeros(G.shape[0])\n",
    "    for row_idx in range(lap.size):\n",
    "        row = G.getrow(row_idx)\n",
    "        l = 0\n",
    "        for _, col_idx, weight in zip(*row.nonzero(), row.data):\n",
    "            l += weight*(f[row_idx]-f[col_idx])\n",
    "        lap[row_idx] = l#v_weight(row_idx)*d\n",
    "    return lap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we are going to define the operators as linear operators. To do so, we will need to defined the incidence matrix (actually a weighted version of it).\n",
    "\n",
    "We recall the definition of the incidence matrix $\\nabla$:\n",
    "\n",
    "\\begin{align*}\n",
    "    \\nabla_{e,v} &= \\begin{cases}\n",
    "        1 \\; \\text{if} \\; v \\; \\text{is initial vertex of edge} \\; e\\\\\n",
    "        -1 \\; \\text{if} \\; v \\; \\text{is terminal vertex of edge} \\; e\\\\\n",
    "        0 \\; \\text{if} \\; e \\; \\text{is not incident on vertex } \\; v\\\\\n",
    "    \\end{cases}\n",
    "\\end{align*}\n",
    "\n",
    "This operator can be considered as a forward difference operator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now define those operators as linear operators\n",
    "\n",
    "def get_incident_matrix(G):\n",
    "    \"\"\"-G is the graph adjacency matrix\n",
    "    Returns the matrix that takes a function f, valued of the the network and returns the gradient\n",
    "    (as a linear index (0,0) (0,1) (0,2)... (1,0) (1,1) (1.2)...\n",
    "    \"\"\"\n",
    "    grad_mat = lil_matrix((np.prod(G.shape), G.shape[0]), dtype=np.float32)\n",
    "    #print('G shape is {} gradmat shape is {}'.format(G.shape,grad_mat.shape))\n",
    "    for i, j, weight in zip(*G.nonzero(), G.data):\n",
    "        weight = np.sqrt(weight)\n",
    "        grad_mat[i*G.shape[1]+j,i]+=weight #+1 if non weighted\n",
    "        grad_mat[i*G.shape[1]+j,j]-=weight #-1 if non weighted\n",
    "    return grad_mat\n",
    "\n",
    "def get_degree_matrix(G):\n",
    "    \"\"\"-G is the graph adjacency matrix\n",
    "    Returns the degree matrix of the graph\n",
    "    \"\"\"\n",
    "    diag = np.squeeze(np.array(G.sum(axis=1)))\n",
    "    return sparse.diags(diag)\n",
    "\n",
    "def grad_with_matrix(G,f):\n",
    "    \"\"\"-G is the graph adjacency matrix\n",
    "       -f is the function defined on this graph\n",
    "       Output has the dimension of the adjacency matrix\n",
    "       This operator can be considered a forward difference operation\n",
    "       matrix version\n",
    "    \"\"\"\n",
    "    return get_incident_matrix(G).dot(lil_matrix(f.reshape(f.size,1))).tolil().reshape(G.shape)\n",
    "\n",
    "def div_with_matrix(G,F,b=None):\n",
    "    \"\"\"-G is the graph adjacency matrix\n",
    "       -F is the gradient of function f on this graph\n",
    "       -b is the optional vertex_weight\n",
    "       Output is 1D and has its size equal to the cardinality of the vertices\n",
    "       matrix version\n",
    "    \"\"\"\n",
    "    d = (-get_incident_matrix(G).T).dot(F.tolil().reshape((np.prod(G.shape),1)))\n",
    "    return np.squeeze(np.array(d.todense()))\n",
    "\n",
    "def laplacian_with_matrix(G,f,b=None):\n",
    "    \"\"\"-G is the graph adjacency matrix\n",
    "       -f is a funcion defined on this graph\n",
    "       -b is the optional vertex_weight\n",
    "       Output is 1D and has its size equal to the cardinality of the vertices\n",
    "       matrix version\n",
    "    \"\"\"\n",
    "    return (get_degree_matrix(G)-G).dot(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check if both grad give the same result:\n",
    "size = 20\n",
    "G1 = gen_random_graph(size)\n",
    "#G1 = gen_small_graph()\n",
    "f1 = np.random.rand(size)\n",
    "grad1 = grad(G1, f1)\n",
    "gradm1 = grad_with_matrix(G1, f1)\n",
    "assert(np.allclose(grad1.todense(), gradm1.todense()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check if both div give the same result:\n",
    "size = 20\n",
    "G1 = gen_random_graph(size)\n",
    "F1 = gen_random_graph(size)\n",
    "div1 = div(G1, F1)\n",
    "divm1 = div_with_matrix(G1, F1)\n",
    "assert(np.allclose(div1, divm1, atol=1e-4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check if both laplacian give the same result\n",
    "size = 20\n",
    "G1 = gen_random_graph(size)\n",
    "f1 = np.random.rand(size)\n",
    "lap1 = laplacian(G1, f1)\n",
    "lapm1 = laplacian_with_matrix(G1, f1)\n",
    "assert(np.allclose(lap1, lapm1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the following assumption, for a given linear operator A, and its transpose B, where we cannot explicitly obtain $A^t$ neither $B^t$:\n",
    "\n",
    "\\begin{align*}\n",
    "    xA^t A y &= x^tBAy \\\\\n",
    "    x^t \\cdot (A^tA y) &= x^t \\cdot (BAy) \\\\\n",
    "    (x^tA^t) \\cdot (A y) &= (x^tB) \\cdot (Ay) \\\\\n",
    "    (Ax)^t \\cdot (Ay) &= (B^tx) \\cdot (Ay)\n",
    "\\end{align*}\n",
    "\n",
    "By comparing $a = x^t \\cdot (BAy)$ and $b = (Ax)^t \\cdot (Ay)$ that both give a scalar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# check if matrix grad and matrix div are transpose of each other (by definition they should...)\n",
    "f1 = np.random.rand(size)\n",
    "f2 = np.random.rand(size)\n",
    "a = f1.dot(-div_with_matrix(G1, grad_with_matrix(G1, f2)))\n",
    "b = grad_with_matrix(G1, f1).multiply(grad_with_matrix(G1, f2)).sum()\n",
    "assert(np.allclose(a,b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check if grad and div are transpose of each other\n",
    "f1 = np.random.rand(size)\n",
    "f2 = np.random.rand(size)\n",
    "\n",
    "a = f1.dot(-div(G1, grad(G1, f2)))\n",
    "b = grad(G1, f1).multiply(grad(G1, f2)).sum()\n",
    "assert(np.allclose(a,b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-14.18685254 -11.46867769 -17.00218327  10.39833276 -18.7286261\n",
      "  -1.93604717  -8.69733653  13.4095348   -4.05484683  14.60066756\n",
      "   9.94081769  -3.71789777  12.25901727   9.55552775  14.60443636\n",
      "   7.4685563  -11.21909322  -2.27185562  -6.73175711   7.77828333]\n",
      "[-7.09342632 -5.73433903 -8.5010918   5.19916655 -9.36431314 -0.96802357\n",
      " -4.34866828  6.70476759 -2.02742345  7.30033387  4.970409   -1.85894895\n",
      "  6.1295084   4.77776394  7.30221818  3.7342782  -5.60954655 -1.13592778\n",
      " -3.36587855  3.88914171]\n",
      "[-7.09342603 -5.73433827 -8.50109178  5.19916642 -9.36431299 -0.96802239\n",
      " -4.34866833  6.70476968 -2.02742277  7.30033465  4.97040805 -1.85894912\n",
      "  6.12950938  4.77776382  7.30221649  3.73427867 -5.6095466  -1.13592795\n",
      " -3.36587824  3.88914024]\n"
     ]
    },
    {
     "ename": "AssertionError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-118-5036c60a8bc6>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0mc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcsgraph\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlaplacian\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mG1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 12\u001b[0;31m \u001b[0;32massert\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mallclose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mAssertionError\u001b[0m: "
     ]
    }
   ],
   "source": [
    "# Check if div(grad()) for matrix based implementation gives laplacian\n",
    "G1 = gen_random_graph(size)\n",
    "#G1 = gen_small_graph()\n",
    "f1 = np.random.rand(size)\n",
    "#f1 = np.random.rand(2)\n",
    "a = -div(G1, grad(G1, f1))\n",
    "print(a)\n",
    "b = laplacian(G1, f1)\n",
    "print(b)\n",
    "c = csgraph.laplacian(G1)\n",
    "print(c.dot(f1))\n",
    "assert(np.allclose(a,b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "dimension mismatch",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-66-3ebeebf277eb>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0mf1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrand\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msize\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0ma\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdiv_with_matrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mG1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgrad_with_matrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mG1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mf1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      5\u001b[0m \u001b[0mb\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlaplacian_with_matrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mG1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mf1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0;32massert\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mallclose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-59-03bd03289157>\u001b[0m in \u001b[0;36mgrad_with_matrix\u001b[0;34m(G, f)\u001b[0m\n\u001b[1;32m     30\u001b[0m        \u001b[0mmatrix\u001b[0m \u001b[0mversion\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     31\u001b[0m     \"\"\"\n\u001b[0;32m---> 32\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mget_edge_grad_matrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mG\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlil_matrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msize\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtolil\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mG\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     33\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     34\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mdiv_with_matrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mG\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mF\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/scipy/sparse/base.py\u001b[0m in \u001b[0;36mdot\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m    320\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    321\u001b[0m         \"\"\"\n\u001b[0;32m--> 322\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mother\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    323\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    324\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mpower\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/scipy/sparse/base.py\u001b[0m in \u001b[0;36m__mul__\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m    437\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0missparse\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    438\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mother\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 439\u001b[0;31m                 \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'dimension mismatch'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    440\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_mul_sparse_matrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    441\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: dimension mismatch"
     ]
    }
   ],
   "source": [
    "# Check if div(grad()) for matrix based implementation gives laplacian\n",
    "f1 = np.random.rand(size)\n",
    "\n",
    "a = div_with_matrix(G1, grad_with_matrix(G1, f1))\n",
    "b = laplacian_with_matrix(G1, f1)\n",
    "assert(np.allclose(a,b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Studying graph thanks to those operators\n",
    "\n",
    "It is interesting to notice that the concept of smoothness/regularity can be extended to graph structured data as well as:\n",
    "\n",
    "smoothness($f$) =\n",
    "\\begin{align}\n",
    "  & \\quad \\langle f, \\Delta f \\rangle_{L^2(\\mathcal{V})} \\\\\n",
    "  = &\\sum_{(i,j) \\in \\mathcal{E} \\; \\text{s.t} \\; i\\neq j} (f_i-f_j)^2 \\\\\n",
    "  \\geq & \\qquad 0\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fourier series on graph\n",
    "\n",
    "Following the definition seen in the previous section, we can see that on cane define eigenfunctions of laplacian operator on various domain including graph, for instance in our case, where we consider our laplacian operator $L$ as: $L=D^{-\\frac{1}{2}}AD^{-\\frac{1}{2}}$ where $A$ is the adjacency matrix and $D$ is the degree matrix of the graph.\n",
    "\n",
    "Hopefully, in an undirected graph, $L$ is a symmetric positive semi definite matrix, hence is diagonalizable in an orthonornmal basis. A function that is defined by a set of values, one per node can then be analyzed over this graph, using eigenvectors of the graph ordered by the norm of their corresponding eigenvalues.\n",
    "\n",
    "On can then consider that projection of the function over the lower eigenvectors corresponds to low \"frequencies\", and conversely, the high \"frequencies\" are obtained when projecting the function over the highest eigenvalues\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import an image, and downsize it in order to make the problem numerically tractable\n",
    "im = sc.misc.imresize(sc.misc.ascent(),size=(128,128))\n",
    "im = im.astype(np.float32)\n",
    "plt.imshow(im, interpolation='nearest', cmap=cm.gray, vmin=0, vmax=255)\n",
    "figsize(8,8)\n",
    "plt.axis(\"off\")\n",
    "colorbar()       # displays the color bar close to the image\n",
    "title('Initial image')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Threshold: valued used to forbid disconnected nodes in the graph of the image\n",
    "thresh= 1e-4\n",
    "\n",
    "#Sigma of the gaussian distance metric\n",
    "sigma = 2\n",
    "sigmaSquare = sigma**2\n",
    "\n",
    "#Size of the stencil\n",
    "stenSize = 4\n",
    "\n",
    "#will contain the value of the edges\n",
    "data=[]\n",
    "rowIdx=[]\n",
    "colIdx=[]\n",
    "\n",
    "\"\"\"\n",
    "Now we are going to fill the graph matrix, neglecting small terms.\n",
    "Distance are actually computed twice, this can be optimized\n",
    "\"\"\"\n",
    "\n",
    "for idxI in range(im.shape[0]):\n",
    "  for idxJ in range(im.shape[1]):\n",
    "    curPixIdx=idxI+im.shape[0]*idxJ\n",
    "    for stenX in range(-stenSize,stenSize+1):\n",
    "      for stenY in range(-stenSize,stenSize+1):\n",
    "        stenIdxX = idxI+stenX\n",
    "        stenIdxY = idxJ+stenY\n",
    "        destPixIdx = stenIdxX+im.shape[0]*stenIdxY\n",
    "        if (stenIdxX>=0)and(stenIdxX<im.shape[0])and\\\n",
    "          (stenIdxY>=0)and(stenIdxY<im.shape[1])and\\\n",
    "          ((stenIdxX!=idxI)and(stenIdxY!=idxJ)):\n",
    "          dist = np.exp(-np.linalg.norm(im[idxI,idxJ]-im[stenIdxX,stenIdxY])**2/(2*sigmaSquare))\n",
    "          dist=max(thresh,dist)\n",
    "          data.append(dist)\n",
    "          rowIdx.append(curPixIdx)\n",
    "          colIdx.append(destPixIdx)\n",
    "\n",
    "\"\"\"\n",
    "   Practical implementation using scipy support for sparse matrices, The compressed sparse row (CSR) format\n",
    "   represents a matrix M by three (one-dimensional) arrays, that respectively contain:\n",
    "      -nonzero values\n",
    "      -the extents of rows\n",
    "      -column indices\n",
    "    This format allows fast row access and matrix-vector multiplications (Mx)\n",
    "\"\"\"\n",
    "A = csr_matrix((np.array(data), (np.array(rowIdx), np.array(colIdx))), dtype=np.float32)\n",
    "\n",
    "#Now the diagonal matrix\n",
    "D = dia_matrix((np.array(A.sum(axis=1)),0),shape=A.shape)\n",
    "\n",
    "#We can define the laplacian: L = D^-0.5 A D^-0.5\n",
    "L = D.power(-0.5)*A*D.power(-0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Compute the eigendecomposition of the laplacian of the graph\n",
    "\"\"\"\n",
    "\n",
    "# Get the first 100 eigenvalues of L, a symmetric positive definite matrix\n",
    "# See http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigsh.html and\n",
    "# http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.linalg.svds.html\n",
    "[w,X] = spl.eigsh(L, k=100) #size is now [imsize,nbcomponent]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Show the spectrum of the image, i.e its expression in the basis made of the first lapacian eigenvectors\n",
    "spectrum = np.abs(np.dot(X.T,im.reshape(im.size)))\n",
    "plt.plot(spectrum)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#Norm of each eigenvector\n",
    "#print(\"Norm is \"+str(np.linalg.norm(X[:,0])))\n",
    "\n",
    "for i in np.arange(X.shape[1])[::-1]:\n",
    "  figure(i)\n",
    "  imshow(np.reshape(X[:,i],im.shape), interpolation='nearest', cmap=cm.gray)"
   ]
  },
  {
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}