coarsening_demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a demo of coarsening code from the preprint: [Graph reduction with spectral and cut guarantees](https://arxiv.org/abs/1808.10650). The final version of the paper will appear at JMLR.\n",
    "\n",
    "Depedencies: pygsp, matplotlib, numpy, scipy, networkx, sortedcontainers\n",
    "\n",
    "This work was kindly supported by the Swiss National Science Foundation (grant number PZ00P2 179981).\n",
    "\n",
    "15 March 2019\n",
    "\n",
    "[Andreas Loukas](https://andreasloukas.wordpress.com)\n",
    "\n",
    "Released under the MIT license "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.core.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:90% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from libraries.coarsening_utils import *\n",
    "import libraries.graph_utils\n",
    "\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pylab as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "import networkx as nx\n",
    "import pygsp as gsp\n",
    "from pygsp import graphs\n",
    "gsp.plotting.BACKEND = 'matplotlib'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Construct the graph ## "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 400\n",
    "\n",
    "G = graphs.BarabasiAlbert(N)\n",
    "if not hasattr(G, 'coords'): \n",
    "    graph = nx.from_scipy_sparse_matrix(G.W)\n",
    "    pos = nx.nx_agraph.graphviz_layout(graph, prog='neato')  \n",
    "    G.set_coordinates(np.array(list(pos.values())))  \n",
    "    G.compute_fourier_basis() # this is for convenience (not really needed by coarsening)\n",
    "    \n",
    "N = G.N\n",
    "L = G.L.toarray()\n",
    "S = get_S(G).T\n",
    "\n",
    "plt.spy(G.W, markersize=0.2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Do some coarsening\n",
    "\n",
    "* Possible methods supported are: 'variation_edges', 'variation_neighborhood', 'algebraic_JC', 'heavy_edge', 'affinity_GS'\n",
    "* $r = 1 - n/N$ is the dimensionality reduction ratio\n",
    "* $k$ is the size of the subspace we are interested in\n",
    "* $k_{max}$ should be set to be larger or equal to $k$. It allows us to measure the quality of coarsening over subspaces larger than $k$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "method = 'variation_neighborhood'  \n",
    "\n",
    "# Parameters\n",
    "r    = 0.6 # the extend of dimensionality reduction (r=0 means no reduction)\n",
    "k    = 5  \n",
    "kmax = int(3*k)\n",
    "        \n",
    "C, Gc, Call, Gall = coarsen(G, K=k, r=r, method=method) \n",
    "metrics = coarsening_quality(G, C, kmax=kmax)\n",
    "n = Gc.N\n",
    "\n",
    "print('{:16} | r: {:1.4}, nedges: {}, levels: {}, epsilon: {:1.4}'.format(method, metrics['r'], metrics['m'], len(Call), metrics['error_subspace'][k-1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize the sequence of coarsening levels\n",
    "\n",
    "* $c+1$ graphs are plotted in total. The leftmost is the original graph and the rightmost is the final coarse graph. \n",
    "* Colors are used to indicate the size of each contraction set $C$: \n",
    "    * green is for $|C|=2$  blue is for $|C|=3$, red is for $|C|=4$, and yellow for $|C|>4$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_coarsening(Gall, Call, size=5, alpha=0.6, title=method);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Various metrics for coarsening quality\n",
    "\n",
    "* $\\epsilon$ is the restricted similarity constant such that, for every $x \\in span(U_k)$ we have $$(1 - \\epsilon) x^\\top L x \\leq x_c^\\top L_c x_c \\leq (1+\\epsilon) x^\\top L x $$   \n",
    "* the eigenvalue error is defined (for every $i = 1, \\ldots, k, \\ldots, kmax$ as $\\frac{\\lambda_i - \\tilde{\\lambda}_i}{\\lambda_i}$\n",
    "* the angle matrix contains the angles between the eigenvectors of $L$ (y-axis) and the lifted eigenvectors of $L_c$. The closer to counter-diagonal it is, the better."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "size = 2.04; fig, axes = plt.subplots(1, 3, figsize=(4*size*3, 3*size)); lineWidth = 1\n",
    "\n",
    "axes[0].plot(np.arange(kmax), np.abs(metrics['error_subspace']), 'or-') \n",
    "axes[0].set_xlabel('$k$'); axes[0].set_ylabel('$\\epsilon$')\n",
    "axes[0].plot( [k, k], [0, max(metrics['error_subspace'])], ':k') \n",
    "\n",
    "axes[1].boxplot(np.abs(metrics['error_eigenvalue'])) \n",
    "axes[1].set_ylabel('relative eigenvalue error')\n",
    "\n",
    "axes[2].imshow(abs(metrics['angle_matrix'][:,0:kmax]) )\n",
    "axes[2].plot( [k, k], [0, kmax], ':w') \n",
    "axes[2].plot( [0, kmax], [k, k], ':w') \n",
    "axes[2].plot( [0, N], [n-1, n-1], ':w') \n",
    "axes[2].set_xlim([0, kmax-1])\n",
    "axes[2].set_ylim([0, kmax-1])\n",
    "axes[2].set_xlabel('Lc eigenvectors lifted'); axes[2].set_ylabel('L eigenvectors');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coarsen and lift a vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "size = 2.04; fig, axes = plt.subplots(1, 4, figsize=(4*size*4, 3*size)); lineWidth = 1\n",
    "\n",
    "# a random smooth signal \n",
    "x = G.U[:,:k] @ np.random.randn(k,1)\n",
    "x = x / np.linalg.norm(x)\n",
    "G.plot_signal(x, ax=axes[0], title='signal')\n",
    "\n",
    "# coarsen it \n",
    "xc = coarsen_vector(x, C)\n",
    "Gc.plot_signal(xc, ax=axes[1], title='coarsened signal')\n",
    "\n",
    "# lift it \n",
    "xp = lift_vector(xc, C)\n",
    "G.plot_signal(xp, ax=axes[2], title='lifted signal')\n",
    "\n",
    "# difference\n",
    "G.plot_signal(np.abs(x-xp), ax=axes[3], title='|x - xp|')\n",
    "\n",
    "print('signal error: {}'.format(np.linalg.norm(x - xp)))"
   ]
  },
  {
   "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}