image.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Energy auto-encoder: learn from an image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import numpy.linalg as la\n",
    "import matplotlib.pyplot as plt\n",
    "from pyunlocbox import functions, solvers\n",
    "import time\n",
    "import matplotlib, pyunlocbox  # For versions only.\n",
    "print('Software versions:')\n",
    "for pkg in [np, matplotlib, pyunlocbox]:\n",
    "    print('  %s: %s' % (pkg.__name__, pkg.__version__))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hyper-parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The $\\lambda$ are the relative importance of each term in the composite objective function.\n",
    "* The sparse code dimensionality $m$ should be greater than $n$ for an overcomplete representation but much smaller than $N$ to avoid over-fitting.\n",
    "* $N_p$ denotes the size of the (square) patches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "l_d = 10\n",
    "l_e = 10\n",
    "l_g = 0\n",
    "m = 200\n",
    "Np = 12"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The set of data vectors $X \\in R^{n \\times N}$ is given by patches extracted from a grayscale image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "image = matplotlib.image.imread('data/barbara.png')\n",
    "(Nx, Ny) = np.shape(image)\n",
    "\n",
    "plt.figure(figsize=(8,5))\n",
    "plt.imshow(image, cmap='gray')\n",
    "plt.title('Source image')\n",
    "plt.show()\n",
    "\n",
    "X = np.zeros((Np**2, Nx*Ny/Np**2))\n",
    "for x in np.arange(Nx/Np):\n",
    "    for y in np.arange(Ny/Np):\n",
    "        X[:,x*y] = image[x:x+Np, y:y+Np].reshape((Np**2,))\n",
    "(n, N) = np.shape(X)\n",
    "\n",
    "print('N = %d samples with dimensionality n = %d (patches of %dx%d).' % (N, n, Np, Np))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given $X \\in R^{n \\times N}$, solve $\\min\\limits_{Z \\in R^{m \\times N}, D \\in R^{n \\times m}, E \\in R^{m \\times n}} \\frac{\\lambda_d}{2} \\|X - DZ\\|_F^2 + \\frac{\\lambda_e}{2} \\|Z - EX\\|_F^2 + \\|Z\\|_1$ s.t. $\\|d_i\\|_2 \\leq 1$, $\\|e_k\\|_2 \\leq 1$, $i = 1, \\ldots, m$, $k = 1, \\ldots, n$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solver numeric parameters.\n",
    "N_outer = 20\n",
    "rtol = 1e-3\n",
    "\n",
    "# Static loss function definitions.\n",
    "g_z = functions.norm_l1()\n",
    "g_de = functions.proj_b2(epsilon=1)  # L2-ball indicator function.\n",
    "\n",
    "# Initialization.\n",
    "Z = np.random.normal(size=(m, N))\n",
    "D = np.random.normal(size=(n, m))\n",
    "E = np.random.normal(size=(m, n))\n",
    "objective_z = []\n",
    "objective_d = []\n",
    "objective_e = []\n",
    "objective_g = []\n",
    "tstart = time.time()\n",
    "\n",
    "# Multi-variate non-convex optimization (outer loop).\n",
    "for k in np.arange(N_outer):\n",
    "    \n",
    "    # Convex minimization for Z.\n",
    "    f_zd = functions.norm_l2(lambda_=l_d/2., A=D, y=X, tight=False)\n",
    "    f_ze = functions.norm_l2(lambda_=l_e/2., y=np.dot(E,X))\n",
    "    f_z = functions.func()\n",
    "    f_z._eval = lambda Z: f_zd.eval(Z) + f_ze.eval(Z)\n",
    "    f_z._grad = lambda Z: f_zd.grad(Z) + f_ze.grad(Z)\n",
    "    L = l_e + l_d * la.norm(np.dot(D.T, D))  # Lipschitz continuous gradient.\n",
    "    solver = solvers.forward_backward(step=1./L, method='FISTA')\n",
    "    ret = solvers.solve([f_z, g_z], Z, solver, rtol=rtol, verbosity='NONE')\n",
    "    Z = ret['sol']\n",
    "    objective_z.extend(ret['objective'])\n",
    "    objective_d.extend(np.zeros(np.shape(ret['objective'])))\n",
    "    objective_e.extend(np.zeros(np.shape(ret['objective'])))\n",
    "    \n",
    "    # Convex minimization for D.\n",
    "    f_d = functions.norm_l2(lambda_=l_d/2., A=Z.T, y=X.T, tight=False)\n",
    "    L = l_d * la.norm(np.dot(Z, Z.T))  # Lipschitz continuous gradient.\n",
    "    solver = solvers.forward_backward(step=1./L, method='FISTA')\n",
    "    ret = solvers.solve([f_d, g_de], D.T, solver, rtol=rtol, verbosity='NONE')\n",
    "    D = ret['sol'].T\n",
    "    objective_d.extend(ret['objective'])\n",
    "    objective_z.extend(np.zeros(np.shape(ret['objective'])))\n",
    "    objective_e.extend(np.zeros(np.shape(ret['objective'])))\n",
    "    \n",
    "    E = D.T\n",
    "    \n",
    "    # Convex minimization for E.\n",
    "    f_e = functions.norm_l2(lambda_=l_e/2., A=X.T, y=Z.T, tight=False)\n",
    "    L = l_e * la.norm(np.dot(X, X.T))  # Lipschitz continuous gradient.\n",
    "    solver = solvers.forward_backward(step=1./L, method='FISTA')\n",
    "    ret = solvers.solve([f_e, g_de], E.T, solver, rtol=rtol, verbosity='NONE')\n",
    "    E = ret['sol'].T\n",
    "    objective_e.extend(ret['objective'])\n",
    "    objective_z.extend(np.zeros(np.shape(ret['objective'])))\n",
    "    objective_d.extend(np.zeros(np.shape(ret['objective'])))\n",
    "    \n",
    "    D = E.T\n",
    "    \n",
    "    # Global objective (the indicators are 0).\n",
    "    objective_g.append(g_z.eval(Z) + f_d.eval(D.T) + f_e.eval(E.T))\n",
    "\n",
    "print('Elapsed time: %d seconds' % (time.time() - tstart))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convergence analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Although the overall multi-variate problem is not convex, it seems to converge toward a solution. As we optimally solve each sub-problem (which are convex), we can guarantee that the global objective function will monotically decrease, which is indeed the case.\n",
    "* Good news: the encoder seems to very well approximate the sparse code (low L2 reconstruction error) !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print('g_z(Z) = %e' % g_z.eval(Z))\n",
    "print('f_z(Z,D) = %e' % f_z.eval(Z))\n",
    "print('f_d(D,Z) = %e' % f_d.eval(D.T))\n",
    "print('f_e(E,Z) = %e' % f_e.eval(E.T))\n",
    "print('g_z(Z) + f_d(D,Z) + f_e(E,Z) = %e' % objective_g[-1])\n",
    "\n",
    "plt.figure(figsize=(8,5))\n",
    "plt.semilogy(np.array(objective_z)[:, 0], label='Z: data term')\n",
    "plt.semilogy(np.array(objective_z)[:, 1], label='Z: prior term')\n",
    "#plt.semilogy(np.sum(objective[:,0:2], axis=1), label='Z: sum')\n",
    "plt.semilogy(np.array(objective_d)[:, 0], label='D: data term')\n",
    "plt.semilogy(np.array(objective_e)[:, 0], label='E: data term')\n",
    "N = np.shape(objective_z)[0]\n",
    "plt.xlim(0, N-1)\n",
    "plt.title('Sub-problems convergence')\n",
    "plt.xlabel('Iteration number (inner loops)')\n",
    "plt.ylabel('Objective function value')\n",
    "plt.grid(True); plt.legend(); plt.show()\n",
    "print('Inner loop: %d iterations' % N)\n",
    "\n",
    "plt.figure(figsize=(8,5))\n",
    "plt.plot(objective_g)\n",
    "N = np.shape(objective_g)[0]\n",
    "plt.xlim(0, N-1)\n",
    "plt.title('Global convergence')\n",
    "plt.xlabel('Iteration number (outer loop)')\n",
    "plt.ylabel('Objective function value')\n",
    "plt.grid(True); plt.show()\n",
    "print('Outer loop: %d iterations\\n' % N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sparse codes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* They can be arbitrary sparse by decreasing $\\lambda_d$ and $\\lambda_e$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nnz = np.count_nonzero(Z)\n",
    "#nnz = np.sum(np.abs(Z) < 1e-4)\n",
    "print('Sparsity of Z: %d non-zero entries out of %d entries, i.e. %.1f%%.' % (nnz, Z.size, 100.*nnz/Z.size))\n",
    "\n",
    "plt.figure(figsize=(8,5))\n",
    "plt.spy(Z, precision=0, aspect='auto')\n",
    "plt.xlabel('N = %d samples' % N)\n",
    "plt.ylabel('m = %d atoms' % m)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dictionary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* All constraints are indeed honored.\n",
    "* Only few atoms are actually used. This can already be seen via the sparse code sparsity pattern.\n",
    "* Learned atoms seem plausible but highly repetitive (keep in mind that the images are normalized)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d = np.sqrt(np.sum(D*D, axis=0))\n",
    "print('Constraints on D: %s' % np.alltrue(d <= 1))\n",
    "\n",
    "plt.figure(figsize=(8,5))\n",
    "plt.semilogy(d, '.')\n",
    "plt.title('Dictionary atom norms')\n",
    "plt.xlabel('Atom [1,m]')\n",
    "plt.ylabel('Norm [0,1]')\n",
    "plt.grid(True); plt.show()\n",
    "plt.show()\n",
    "\n",
    "plt.figure(figsize=(8,5))\n",
    "plt.spy(D, precision=1e-7)\n",
    "plt.xlabel('m = %d atoms' % m)\n",
    "plt.ylabel('data dimensionality of n = %d' % n)\n",
    "plt.show()\n",
    "\n",
    "#plt.scatter to show intensity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,8))\n",
    "Nx = np.ceil(np.sqrt(m))\n",
    "Ny = np.ceil(m / float(Nx))\n",
    "for k in np.arange(m):\n",
    "    plt.subplot(Ny, Nx, k)\n",
    "    img = D[:,k].reshape(Np,Np)\n",
    "    plt.imshow(img, cmap='gray')  # vmin=0, vmax=1 to disable normalization.\n",
    "    plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Encoder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* All constraints are indeed honored.\n",
    "* Depending on the conditions, the encoder ressembles the transpose of the dictionary, i.e. $E \\approx D^T$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "e = np.sqrt(np.sum(E*E, axis=0))\n",
    "print('Constraints on E: %s' % np.alltrue(e <= 1))\n",
    "\n",
    "plt.figure(figsize=(8,5))\n",
    "plt.semilogy(e)\n",
    "plt.title('Encoder norms')\n",
    "plt.xlabel('[1,n]')\n",
    "plt.ylabel('Norm [0,1]')\n",
    "plt.grid(True); plt.show()\n",
    "plt.show()\n",
    "\n",
    "plt.figure(figsize=(8,5))\n",
    "plt.spy(E, precision=1e-7)\n",
    "plt.xlabel('data dimensionality of n = %d' % n)\n",
    "plt.ylabel('m = %d atoms' % m)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}