IFT 6760A - Homework 1 (practical part).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# IFT 6760A - Homework 1 - Practical part"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 0: required packages\n",
    "\n",
    "This homework should be done in python 3 using this notebook. If you do not have python installed on your machine, then the [Anaconda distribution](https://www.anaconda.com/distribution/) is the way to go (make sure to choose the 3.7 version). This distribution contains most of the commonly used packages.  \n",
    "\n",
    "You will also need the `tensorly` extra packages (which you can install easily through pip: `pip install tensorly` from the command line).\n",
    "\n",
    "Run the cell below to make sure you have all the right packages installed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pylab as plt\n",
    "import scipy.sparse.linalg \n",
    "import PIL.Image\n",
    "import tensorly as tl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part I: The magical `einsum`function\n",
    "\n",
    "*Einstein notation* is a convenient way to lighten mathematical expressions involving summations and relies on the simple convention that indices appearing twice in an expression are summed over (this is actually a naive simplification where we omit the distinction between so-called _covariant_ and _contravariant_ indices). \n",
    "\n",
    "For example, the inner product between two vectors $u$ and $v$ is written as\n",
    "$$u_iv_i := \\sum_i u_iv_i.$$\n",
    "\n",
    "When an index appears only once, it corresponds to free indices of the resulting tensor. For example $A_{ik}B_{kj}$ has two free indices, $i$ and $j$, and thus corresponds to a second-tensor~(aka matrix) whose components $(i,j)$ is given by $\\sum_k A_{ik}B_{kj}$ (aka the matrix product of $A$ and $B$). Similarly, the outer product of two vectors $u$ and $v$ is noted $u_iv_j$.\n",
    "\n",
    "\n",
    "The `einsum` function from the numpy package is inspried from this convention and makes it easy to write contractions between tensors. The function uses *signature* strings to specify contractions between tensors. For example, `'ij,jk->ki'` denotes a contraction between two 2nd-order tensors~(i.e., matrices $A$ and $B$), the indices of the first one are $i$ and $j$ and the indices of the second one are $j$ and $k$. Since $j$ appears twice it is summed over and the indices appearing after the right arrow `->` specify the order of the indices in the resulting tensor. Thus, `'ij,jk->ki'` corresponds to $(AB)^\\top$.\n",
    "\n",
    "More precisely, a signature strings is interpreted using these two simple rules:\n",
    "\n",
    "- if an index appears in several input arrays, it means that the corresponding modes (are axes in numpy terminology) should be mutliplied together\n",
    "- if an index appears at the left of the arrow sign `->` but not at the right, it means that is should be summed over\n",
    "\n",
    "\n",
    "Here are a couple examples (`u`, `v` and `w` are numpy 1d arrays, `A` and `B` are two 2d arrays and `T` is a 3d array):\n",
    "\n",
    "- inner product: `np.einsum('i,i->',u,v)`\n",
    "- outer product: `np.einsum('i,j->ij',u,v)`\n",
    "- outer product between three vectors: `np.einsum('i,j,k->ijk',u,v,w)`\n",
    "- \"3-way contraction\" between three vectors: `np.einsum('i,i,i->',u,v,w)`\n",
    "- sum of elements of a vector: `np.einsum('i->')\n",
    "- transpose: `np.einsum('ij->ji',A)`\n",
    "- trace of a matrix: `np.einsum('ii->',A)`\n",
    "- diagonal of a matrix: `np.einsum('ii->i',A)`\n",
    "- element-wise multiplication of two matrices: `np.einsum('ij,ij->ij',A,B)`\n",
    "- sum of the product of the diagonal elements of two matrices: `np.einsum('ii,ii->',A,B)`\n",
    "- mode-2 product: `np.einsum('ijk,lj->ilk',T,A)`\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**\n",
    "\n",
    "Write the following functions with only one line using the `einsum` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "def inner_product(S,T):\n",
    "    \"\"\"\n",
    "    S,T: two tensors of order 3\n",
    "    \n",
    "    returns the inner product between the two\n",
    "    \"\"\"\n",
    "    pass\n",
    "\n",
    "def tucker_to_tensor(G,A,B,C):\n",
    "    \"\"\"\n",
    "    G: 3rd order core tensor \n",
    "    A,B,C: factor matrices\n",
    "    \n",
    "    returns the tensor with Tucker decomposition $$G \\times_1 A \\times_2 B \\times_3 C$$\n",
    "    \"\"\"\n",
    "    pass\n",
    "\n",
    "def cp_to_tensor(A,B,C):\n",
    "    \"\"\"\n",
    "    A,B,C: factor matrices (same number of columns in each)\n",
    "    \n",
    "    returns the tensor with CP decomposition $$[[A,B,C]]$$\n",
    "    \"\"\"\n",
    "    pass\n",
    "\n",
    "def tt_to_tensor(A,B,C,D):\n",
    "    \"\"\"\n",
    "    A,B,C,D: core tensors (A and D are order 2, and B and C are order 3)\n",
    "    \n",
    "    returns the tensor with tensor train decomposition $$<<A,B,C,D>>$$\n",
    "    \"\"\"\n",
    "    pass\n",
    "    \n",
    "def tr_to_tensor(A,B,C,D):\n",
    "    \"\"\"\n",
    "    A,B,C,D: core tensors (all of order 3)\n",
    "    \n",
    "    returns the tensor with tensor ring decomposition $$((A,B,C,D))$$. Recall that the tensor ring decomposition\n",
    "    is defined by\n",
    "    ((A,B,C,D))_{ijkl} = \\sum_{r_1,r_2,r_3,r_4} A_{r_1,i,r_2}B_{r_2,i,r_3}C_{r_3,i,r_4}D_{r_4,l,r_1}\n",
    "    \"\"\"\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the following cell to check whether your implementation of these functions is correct:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "G = tl.tensor(np.random.normal(0,1,[10,11,12]))\n",
    "T = tl.tensor(np.random.normal(0,1,[10,11,12]))\n",
    "A,B,C = [np.random.normal(0,1, [d,4]) for d in [10,11,12]]\n",
    "D,E,F = [np.random.normal(0,1, [d,r]) for d,r in zip([15,15,14],[10,11,12])]\n",
    "G1,G2,G3,G4 = [np.random.normal(0,1, [r1,d,r2]) for r1,d,r2 in zip([3,4,5,6],[9,8,7,8],[4,5,6,3])]\n",
    "\n",
    "\n",
    "try:\n",
    "    print(\"inner_product: ok\" if np.allclose(tl.tenalg.inner(G,T), inner_product(G,T)) else \"inner_product: failed\")\n",
    "except Exception as e:\n",
    "    print(\"inner_product: exception raised \\n\",e)\n",
    "print()                \n",
    "try:\n",
    "    print(\"tucker_to_tensor: ok\" if np.allclose(tl.tucker_to_tensor((G,[D,E,F])), tucker_to_tensor(G,D,E,F)) else \"tucker_to_tensor: failed\")\n",
    "except Exception as e:\n",
    "    print(\"tucker_to_tensor: exception raised \\n\",e)\n",
    "print() \n",
    "try:\n",
    "    print(\"cp_to_tensor: ok\" if np.allclose(tl.kruskal_to_tensor((np.ones(4),(A,B,C))), cp_to_tensor(A,B,C)) else \"cp_to_tensor: failed\")\n",
    "except Exception as e:\n",
    "    print(\"cp_to_tensor: exception raised \\n\",e)\n",
    "print() \n",
    "try:\n",
    "    print(\"tt_to_tensor: ok\" if np.allclose(tl.mps_to_tensor((G1[:1,:,:],G2,G3,G4[:,:,:1])), tt_to_tensor(G1[0,:,:],G2,G3,G4[:,:,0])) else \"tt_to_tensor: failed\")\n",
    "except Exception as e:\n",
    "    print(\"tt_to_tensor: exception raised \\n\",e)\n",
    "print() \n",
    "try:\n",
    "    print(\"tr_to_tensor: ok\" if np.allclose(sum([tl.mps_to_tensor([G1[i:i+1,:,:],G2,G3,G4[:,:,i:i+1]]) for i in range(3)]), tr_to_tensor(G1,G2,G3,G4)) else \"tr_to_tensor: failed\")\n",
    "except Exception as e:\n",
    "    print(\"tr_to_tensor: exception raised \\n\",e)    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part II: Alternating Least Squares and Gradient Descent for the CP Decomposition\n",
    "\n",
    "In this first part of the homework, you will compare two optimization algorithms to compute approximate solutions for the problem of approximating a tensor using the CP decomposition: Alternating Least Squares and Gradient Descent.\n",
    "\n",
    "For a given tensor $T$ of size $d_1\\times d_2 \\times d_3$ and a target rank $R$, the low rank CP decomposition problem can be formulated as follows:\n",
    "\n",
    "$${\\arg\\min}_{A\\in \\mathbb{R}^{d_1\\times R},B\\in \\mathbb{R}^{d_2\\times R}, C\\in \\mathbb{R}^{d_3\\times R}} \\| T - [\\!| A,B,C |\\!]\\|_{F}^2$$\n",
    "where $[\\!| A,B,C |\\!]$ denotes the CP decomposition with factor matrices $A,B$ and $C$.\n",
    "\n",
    "\n",
    "**Question** Start by implementing ALS by filling the code skeleton below. To reconstruct a CP tensor from the 3 factor matrices $A,B,C$, use the function cp_to_tensor you wrote in the previous question.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0. 0. 0. ... 0. 0. 0.]\n",
      " [0. 0. 0. ... 0. 0. 0.]\n",
      " [0. 0. 0. ... 0. 0. 0.]\n",
      " ...\n",
      " [0. 0. 0. ... 0. 0. 0.]\n",
      " [0. 0. 0. ... 0. 0. 0.]\n",
      " [0. 0. 0. ... 0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "\n",
    "\n",
    "def khatri_rao(A,B):\n",
    "    assert A.shape[1] == B.shape[1], \"A and B must have the same number of columns\"\n",
    "    return np.einsum('ir,jr->jir',A,B).reshape(A.shape[0]*B.shape[0],A.shape[1])\n",
    "\n",
    "\n",
    "\n",
    "A,B,C = [np.random.normal(0,1, [d,4]) for d in [10,11,12]]\n",
    "print(tl.unfold(tl.kruskal_to_tensor((np.ones(4),(A,B,C))),0) - A.dot(khatri_rao(C,B).T))\n",
    "\n",
    "def CP_ALS(T, rank, epsilon=1e-5,max_iters=100,verbose=True,errors_at_iterations=[]):\n",
    "    # Initialize A, B and C\n",
    "    # YOUR CODE HERE\n",
    "    \n",
    "    T_approx_old = tl.kruskal_to_tensor((np.ones(rank),[A,B,C]))\n",
    "    it = 1\n",
    "    while True: # repeat until convergence\n",
    "        # update the factor matrices A, B and C. Here are a couple of indications:\n",
    "        # - to efficiently solve a least square problem you can use numpy.lstsq\n",
    "        # - matricizations (unfolding) of T can be obtained using tl.unfold\n",
    "        # - to compute the Khatri-Rao product use the function khatri_rao defined above\n",
    "        # YOUR CODE HERE\n",
    "        \n",
    "        # check for convergence:\n",
    "        T_approx_new = tl.kruskal_to_tensor((np.ones(rank),[A,B,C]))\n",
    "        reconstruction_error =  tl.norm(T_approx_new-T)\n",
    "        errors_at_iterations.append(reconstruction_error)        \n",
    "        if tl.norm(T_approx_old - T_approx_new)/tl.norm(T_approx_old) < epsilon or it > max_iters:\n",
    "            return (A,B,C)\n",
    "        if verbose:\n",
    "            print(it, \" - reconstruction error: \",reconstruction_error, \"convergence: \", tl.norm(T_approx_old - T_approx_new)/tl.norm(T_approx_old))\n",
    "\n",
    "        \n",
    "        T_approx_old = T_approx_new\n",
    "        it += 1\n",
    "        \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Plot the reconstruction error as a function of the iterations for 10 different initialization for ALS to find an approximate CP decomposition of the tensor T below of rank 2. Plot all the curves in the same figure and use a maximum number of iterations of 100. Observations?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "T = tl.tensor(np.random.normal(0,1,[5,10,3]))\n",
    "\n",
    "# - to plot a curve, use the plt.plot function\n",
    "# - once all plt.plot functions have been run, use plt.show() to show the figure\n",
    "# - to recover the errors at each iteration, you can use the following code:\n",
    "# errors = []\n",
    "# CP_ALS(...[other arguments here]..., errors_at_iterations=errors)\n",
    "# - when the method CP_ALS returns, errors while be a list containing the error at each iteration\n",
    "\n",
    "# YOUR CODE HERE\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Write your observations in this cell*\n",
    "\n",
    "... ... ...\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Now take inspiration from the CP_ALS function you wrote above to implement gradient descent to solve the CP low rank approximation problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def CP_GD(T, rank, epsilon=1e-5,max_iters=100,lr=1e-2,verbose=True,errors_at_iterations=[]):\n",
    "    # Initialize A, B and C\n",
    "    # YOUR CODE HERE\n",
    "    \n",
    "    \n",
    "    T_approx_old = # YOUR CODE HERE #\n",
    "    it = 1\n",
    "    while True: # repeat until convergence\n",
    "        # compute gradients wrt A B and C and do a gradient step\n",
    "        # YOUR CODE HERE\n",
    "        \n",
    "        \n",
    "        # check for convergence:\n",
    "        T_approx_new = # YOUR CODE HERE #\n",
    "        reconstruction_error =  # YOUR CODE HERE #\n",
    "        errors_at_iterations.append(reconstruction_error)\n",
    "        \n",
    "                \n",
    "        if tl.norm(T_approx_old - T_approx_new)/tl.norm(T_approx_old) < epsilon or it > max_iters:\n",
    "            return (A,B,C)\n",
    "        if verbose:\n",
    "            print(it, \" - reconstruction error: \",reconstruction_error, \"convergence: \", tl.norm(T_approx_old - T_approx_new)/tl.norm(T_approx_old))\n",
    "        \n",
    "        T_approx_old = T_approx_new\n",
    "        it += 1\n",
    "        \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Plot the reconstruction error as a function of the iterations for 10 different initialization of ALS and 10 initializations of gradient descent to find an approximate CP decomposition of the tensor T below of rank 2. Plot all the curves in the same figure and use a maximum number of iterations of 100. Observations?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "T = tl.tensor(np.random.normal(0,1,[5,10,3]))\n",
    "\n",
    "### YOUR CODE HERE\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Write your observations in this cell*\n",
    "\n",
    "... ... ...\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part III: Image Compression with Matrix and Tensor Factorization\n",
    "\n",
    "In this part of the homework you will compare 5 different ways to compress an image using matrix and tensor factorization techniques:\n",
    "\n",
    "- SVD\n",
    "- CP decomposition\n",
    "- Tucker decomposition \n",
    "- Tensor Train decomposition\n",
    "- Tensor Train decomposition with reshaping to high-order tensor\n",
    "\n",
    "For the tensor decomposition methods, we will use the `tensorly` python package.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use two images as example for this part of the homework: an image of a building, and a of a racoon. Run the cell below to see the two images. Note that those images are color images, thus they are tensors of order 3 (`width`$\\times$ `height`$\\times 3$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget -nc https://www-labs.iro.umontreal.ca/~grabus/files/building2.jpg\n",
    "!wget -nc https://www-labs.iro.umontreal.ca/~grabus/files/racoon.jpg\n",
    "\n",
    "\n",
    "image1 = tl.tensor(np.array(PIL.Image.open('building2.jpg')), dtype='float64')/255\n",
    "image2 = tl.tensor(np.array(PIL.Image.open('racoon.jpg')), dtype='float64')/255\n",
    "\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.subplot(1,2,1)\n",
    "imshow(image1)\n",
    "plt.title(\"building.jpg\")\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "imshow(image2)\n",
    "plt.title(\"racoon.jpg\")\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** First you will use a simple low rank approximation with SVD. For this you will first unfold the image/tensor into a matrix of shape `width`$\\times 3$`height`) and then perform trunated SVD usimg the `svds` function from the `scipy.sparse.linalg` pakcage.\n",
    "\n",
    "Complete the code skeleton below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.sparse.linalg import svds\n",
    "\n",
    "def svd_compression(image, rank):\n",
    "    \"\"\"\n",
    "    image is a third order tensor of shape width x height x 3\n",
    "    rank is the rank of the truncated SVD\n",
    "    \n",
    "    This function should return an order 3 tensor of shape width x height x 3 corresponding to the initial \n",
    "    tensor compressed using truncated SVD on the appropriate matricization/unfolding \n",
    "    \"\"\"\n",
    "    # - Recall that unfolding a tensor can be done with the tl.unfold function\n",
    "    # - The inverse operation (i.e., reconstructing a tensor from an unfolding/matricization) can be done\n",
    "    # with the tl.fold function.\n",
    "    \n",
    "    ### YOUR CODE HERE\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Complete the code skeleton below to write a function that returns the compression ratio achieved by truncated SVD for a given input image and target rank. The compression ratio is given by the number of parameters of the low rank factorization divided by the number of parameters in the orginal image. *You should not count the size of the diagonal matrix of singular values in this computation since it can be incorporated in one of the matrices of singular vectors*.\n",
    "\n",
    "At the end of the cell, a test will let you know whether your implementation is correct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def svd_compress_ratio(image, rank):\n",
    "    \"\"\"\n",
    "    This function should return the compression ratio achieved by truncated SVD on the input image\n",
    "    using the input rank.\n",
    "    Note: this function *should not call the svd_compression function you wrote previously*, it should computes\n",
    "    the compression ratio using the size of the image, the rank, and your knowledge of the SVD.\n",
    "    \"\"\"\n",
    "    ### YOUR CODE HERE\n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "try:\n",
    "    image = np.zeros([200,300,3])\n",
    "    print(\"test passed\" if svd_compress_ratio(image,3)== 3300/(200*300*3) else \"test failed\")\n",
    "except Exception as e:\n",
    "    print(\"exception raised \\n\",e)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the cell below to see the compressed images using different values for the rank of the truncated SVD."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20,12))\n",
    "for i,rank in enumerate([2,10,20]):\n",
    "    plt.subplot(2,3,i+1)\n",
    "    imshow(svd_compression(image1,rank))\n",
    "    plt.title(\"rank:%i - compression: %.2f%%\" %(rank,svd_compress_ratio(image1,rank)*100))\n",
    "for i,rank in enumerate([2,10,20]):\n",
    "    plt.subplot(2,3,i+4)\n",
    "    imshow(svd_compression(image2,rank))\n",
    "    plt.title(\"rank:%i - compression: %.2f%%\" %(rank,svd_compress_ratio(image2,rank)*100))\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you will use a CP decomposition. For this, you will use the function `cp_decomposition` defined below. To reconstruct a tensor from the factor matrices, use the function `cp_to_tensor` you wrote in the previous question.\n",
    "\n",
    "Similarly to truncated SVD, you need to write a function which returns the number of parameters of a CP decomposition (a simple test is executed at the end of the celle to check that your implementation is correct).\n",
    "\n",
    "**Question** Complete the code skeleton below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from tensorly.decomposition import parafac\n",
    "\n",
    "def cp_decomposition(tensor, rank):\n",
    "    (weights,factors) = parafac(tensor, rank=rank, init='random', tol=10e-6)\n",
    "    L = [factors[0]@np.diag(weights)] + factors[1:]\n",
    "    return L\n",
    "\n",
    "\n",
    "def cp_compression(image, rank):\n",
    "    \"\"\"\n",
    "    image is a third order tensor of shape width x height x 3\n",
    "    rank is the rank of the CP decomposition\n",
    "    \n",
    "    This function should return an order 3 tensor of shape width x height x 3 corresponding to the initial \n",
    "    tensor compressed using CP decomposition\n",
    "    \"\"\"\n",
    "    ### YOUR CODE HERE\n",
    "    \n",
    "    \n",
    "    \n",
    "def cp_compress_ratio(image, rank):\n",
    "    \"\"\"\n",
    "    This function should return the compression ratio achieved by CP on the input image\n",
    "    using the input rank.\n",
    "    Note: this function *should not call the cp_compression function you wrote previously*, it should computes\n",
    "    the compression ratio using the size of the image, the rank, and your knowledge of the definition of the \n",
    "    CP decomposition.\n",
    "    \"\"\"\n",
    "    ### YOUR CODE HERE\n",
    "    \n",
    "\n",
    "\n",
    "try:\n",
    "    image = np.zeros([200,300,3])\n",
    "    print(\"test passed\" if cp_compress_ratio(image,4)== 2012/(200*300*3) else \"test failed\")\n",
    "except Exception as e:\n",
    "    print(\"exception raised \\n\",e)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the cell below to see the compressed images using different values for the rank of the CP decomposition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20,12))\n",
    "for i,rank in enumerate([2,10,20]):\n",
    "    plt.subplot(2,3,i+1)\n",
    "    imshow(cp_compression(image1,rank))\n",
    "    plt.title(\"rank:%i - compression: %.2f%%\" %(rank,cp_compress_ratio(image1,rank)*100))\n",
    "for i,rank in enumerate([2,10,20]):\n",
    "    plt.subplot(2,3,i+4)\n",
    "    imshow(cp_compression(image2,rank))\n",
    "    plt.title(\"rank:%i - compression: %.2f%%\" %(rank,cp_compress_ratio(image2,rank)*100))\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Now you will use a Tucker decomposition. Use the function `tucker` from the `tensorly` package  and the function `tucker_to_tensor` to fill in the code skeleton below. \n",
    "\n",
    "For the rank parameter, we will use a unform rank for the two first modes of the decomposition and use 3 for the rank of the 3rd mode of the decomposition. E.g., if the rank argument is 10, the Tucker decomposition should be of rank (10,10,3).\n",
    "\n",
    "As before, you need to write a function which returns the number of parameters of a Tucker decomposition (a simple test is executed at the end of the celle to check that your implementation is correct).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tensorly.decomposition import tucker\n",
    "\n",
    "def tucker_compression(image, rank):\n",
    "    \"\"\"\n",
    "    image is a third order tensor of shape width x height x 3\n",
    "    rank is the rank of the Tucker decomposition for the first two modes, the third rank should be 3 \n",
    "    (i.e., the Tucker rank of the tucker_decomposition is (rank,rank,3))\n",
    "    \n",
    "    This function should return an order 3 tensor of shape width x height x 3 corresponding to the initial \n",
    "    tensor compressed using Tucker decomposition\n",
    "    \"\"\"\n",
    "    ### YOUR CODE HERE\n",
    "    \n",
    "    \n",
    "    \n",
    "def tucker_compress_ratio(image, rank):\n",
    "    \"\"\"\n",
    "    This function should return the compression ratio achieved by Tucker on the input image\n",
    "    using the input rank.\n",
    "    Note 1: this function *should not call the tucker_compression function you wrote previously*, it should computes\n",
    "    the compression ratio using the size of the image, the rank, and your knowledge of the definition of the \n",
    "    Tucker decomposition.\n",
    "    \"\"\"\n",
    "    ### YOUR CODE HERE\n",
    "    \n",
    "\n",
    "\n",
    "try:\n",
    "    image = np.zeros([200,300,3])\n",
    "    print(\"test passed\" if tucker_compress_ratio(image,4)== 2057/(200*300*3) else \"test failed\")\n",
    "except Exception as e:\n",
    "    print(\"exception raised \\n\",e)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the cell below to see the compressed images using different values for the rank of the Tucker decomposition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20,12))\n",
    "for i,rank in enumerate([2,10,20]):\n",
    "    plt.subplot(2,3,i+1)\n",
    "    imshow(tucker_compression(image1,rank))\n",
    "    plt.title(\"rank:%i - compression: %.2f%%\" %(rank,tucker_compress_ratio(image1,rank)*100))\n",
    "for i,rank in enumerate([2,10,20]):\n",
    "    plt.subplot(2,3,i+4)\n",
    "    imshow(tucker_compression(image2,rank))\n",
    "    plt.title(\"rank:%i - compression: %.2f%%\" %(rank,tucker_compress_ratio(image2,rank)*100))\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Now you will use a TT decomposition. Use the function `matrix_product_states` and `mps_to_tensor` from the `tensorly` package to fill in the code skeleton below. As before, you need to write a function which returns the number of parameters of a TT decomposition (a simple test is executed at the end of the celle to check that your implementation is correct).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tensorly.decomposition import matrix_product_state\n",
    "from tensorly import mps_to_tensor\n",
    "\n",
    "def tt_compression(image, rank):\n",
    "    \"\"\"\n",
    "    image is a third order tensor of shape width x height x 3\n",
    "    rank is the rank of the TT decomposition (i.e., the TT rank of the TT_decomposition is (rank,rank,...,rank))\n",
    "    \n",
    "    This function should return an order 3 tensor of shape width x height x 3 corresponding to the initial \n",
    "    tensor compressed using TT decomposition\n",
    "    \"\"\"\n",
    "    ### YOUR CODE HERE\n",
    "    \n",
    "    \n",
    "    \n",
    "def tt_compress_ratio(image, rank):\n",
    "    \"\"\"\n",
    "    This function should return the compression ratio achieved by TT on the input image\n",
    "    using the input rank.\n",
    "    Note 1: this function *should not call the tt_compression function you wrote previously*, it should computes\n",
    "    the compression ratio using the size of the image, the rank, and your knowledge of the definition of the \n",
    "    TT decomposition.\n",
    "    Note 2: remember that the ranks of a TT decomposition are bounded by the dimension of some \n",
    "    unfolding/matricizations of the input tensor. You need to take this into account in your computations!\n",
    "    \"\"\"\n",
    "    ### YOUR CODE HERE\n",
    "    \n",
    "\n",
    "\n",
    "try:\n",
    "    image = np.zeros([200,300,3])\n",
    "    print(\"test passed\" if tt_compress_ratio(image,4)== 4409/(200*300*3) else \"test failed\")\n",
    "except Exception as e:\n",
    "    print(\"exception raised \\n\",e)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the cell below to see the compressed images using different values for the rank of the TT decomposition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20,12))\n",
    "for i,rank in enumerate([2,10,20]):\n",
    "    plt.subplot(2,3,i+1)\n",
    "    imshow(tt_compression(image1,rank))\n",
    "    plt.title(\"rank:%i - compression: %.2f%%\" %(rank,tt_compress_ratio(image1,rank)*100))\n",
    "for i,rank in enumerate([2,10,20]):\n",
    "    plt.subplot(2,3,i+4)\n",
    "    imshow(tt_compression(image2,rank))\n",
    "    plt.title(\"rank:%i - compression: %.2f%%\" %(rank,tt_compress_ratio(image2,rank)*100))\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To conclude this homework, we will see a last decomposition technique where we reshape the input image into a high order tensor before applying TT decomposition for the compression.\n",
    "\n",
    "For this we will implement a particular way of reshaping a tensor which takes into account the local correlation in images. The technique is borrowed from https://arxiv.org/abs/1807.01589:\n",
    "\n",
    "*The idea is simple, if the\n",
    "size of a RGB image is U ×V ×3, and U = u1×u2×· · ·×ul\n",
    "and V = v1 × v2 × · · · × vl are satisfied, then the image\n",
    "can be tensorized to a (l + 1)-dimension tensor of size\n",
    "u1v1 × u2v2 × · · · × ulvl × 3. The first dimension of the\n",
    "higher-dimension tensor stands for a pixel block of the image,\n",
    "and the following dimensions are the extension blocks of the\n",
    "image. The higher-dimension tensor generated by the above\n",
    "tensorization scheme is considered to be a better structure of\n",
    "the image.*\n",
    "\n",
    "**Question** Write two functions: \n",
    "- one that takes an image as input with a list of dimensions to factorize the rows and a list of dimensions to factorize the columns, and returns the tensor obtained by reshaping the image using the scheme described above.\n",
    "- a second that implements the reverse operation. \n",
    "\n",
    "Note that to properly implement this function, you need to first reshape the image into u1 × u2× ... × ul x v1 × v2 × ... × vl x 3, then transpose the tensor into shape u1 x v1 x u2 x v2 x ... x ul x vl x 3, before finally reshaping into shape u1v1 × u2v2 × · · · × ulvl × 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def reshape_image_to_tensor(img, dims_rows, dims_cols):\n",
    "    \"\"\"\n",
    "    img: a 3rd order tensor of shape width x height x 3\n",
    "    dims_rows: a list of dimension u1,u2,...,ul such that width=u1u2...ul\n",
    "    dims_cols: a list of dimension v1,v2,...,vl such that width=v1v2...vl\n",
    "    \n",
    "    returns a tensor of shape u1v1 x u2v2 x ... x ulvl x 3\n",
    "    \"\"\"\n",
    "    ### YOUR CODE HERE\n",
    "    \n",
    "\n",
    "def reshape_tensor_to_image(tensor, dims_rows, dims_cols):\n",
    "    \"\"\"\n",
    "    tensor: a tensor of shape u1v1 x u2v2 x ... x ulvl x 3\n",
    "    dims_rows: a list of dimension u1,u2,...,ul such that width=u1u2...ul\n",
    "    dims_cols: a list of dimension v1,v2,...,vl such that width=v1v2...vl\n",
    "    \n",
    "    returns a tensor T of shape width x height x 3 such that for an initial image img:\n",
    "    img = reshape_tensor_to_image(reshape_image_to_tensor(img, dims_rows, dims_cols), dims_rows, dims_cols)\n",
    "    \"\"\"\n",
    "    \n",
    "    ### YOUR CODE HERE\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question** Now you will use the TT decomposition on the reshaped tensor. Use the function `matrix_product_states` and `mps_to_tensor` from the `tensorly` package as well as the functions `reshape_image_to_tensor` and `reshape_tensor_to_image` you wrote previously to fill in the code skeleton below. As before, you need to write a function which returns the number of parameters of the TT decomposition, on the reshaping (a simple test is executed at the end of the celle to check that your implementation is correct).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tensorly.decomposition import matrix_product_state\n",
    "from tensorly import mps_to_tensor\n",
    "\n",
    "def tt_reshape_compression(image, rank, dims_rows, dims_cols):\n",
    "    \"\"\"\n",
    "    image is a third order tensor of shape width x height x 3\n",
    "    rank is the rank of the TT decomposition (i.e., the TT rank of the TT_decomposition is (rank,rank,...,rank))\n",
    "    dims_rows anddims_cols are the lists of dimensions to use for the reshaping operation discussed above.\n",
    "    \n",
    "    This function should return an order 3 tensor of shape width x height x 3 corresponding to the initial \n",
    "    tensor compressed using TT decomposition on the reshaping of the image\n",
    "    \"\"\"\n",
    "    ### YOUR CODE HERE\n",
    "    \n",
    "def tt_reshape_compress_ratio(image, rank,dims_rows, dims_cols):\n",
    "    \"\"\"\n",
    "    This function should return the compression ratio achieved by TT on the input image\n",
    "    using the input rank.\n",
    "    Note 1: this function *should not call the tt_compression function you wrote previously*, it should computes\n",
    "    the compression ratio using the size of the image, the rank, and your knowledge of the definition of the \n",
    "    TT decomposition.\n",
    "    Note 2: remember that the ranks of a TT decomposition are bounded by the dimension of some \n",
    "    unfolding/matricizations of the input tensor. You need to take this into account in your computations!\n",
    "    \"\"\"\n",
    "    ### YOUR CODE HERE\n",
    "\n",
    "try:\n",
    "    image = np.zeros([200,300,3])\n",
    "    print(\"test passed\" if tt_reshape_compress_ratio(image,4,[2,5,20],[3,20,5])== 2833/(200*300*3) else \"test failed\")\n",
    "except Exception as e:\n",
    "    print(\"exception raised \\n\",e)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the cell below to see the compressed images using different values for the rank of the TT decomposition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20,12))\n",
    "\n",
    "for i,rank in enumerate([2,10,20]):\n",
    "    plt.subplot(2,3,i+1)\n",
    "    dims_rows,dims_cols=[2,2,2,2,2,2,2,2],[2,2,2,2,2,2,2,2]\n",
    "    imshow(tt_reshape_compression(image1,rank,dims_rows,dims_cols))\n",
    "    plt.title(\"rank:%i - compression: %.2f%%\" %(rank,tt_reshape_compress_ratio(image1,rank,dims_rows,dims_cols)*100))\n",
    "\n",
    "for i,rank in enumerate([2,10,20]):\n",
    "    plt.subplot(2,3,i+4)\n",
    "    dims_rows,dims_cols=[2,2,2,2,2,2,2,2],[2,2,2,2,2,2,2,2]\n",
    "    imshow(tt_reshape_compression(image2,rank,dims_rows,dims_cols))\n",
    "    plt.title(\"rank:%i - compression: %.2f%%\" %(rank,tt_reshape_compress_ratio(image2,rank,dims_rows,dims_cols)*100))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To conclude, you can run the cell below to compare the various decomposition techniques comparable approximation ratio (change the ratio variable value to compare different ones)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ratio = 5/100\n",
    "\n",
    "\n",
    "def compress_ratio(image,method,rank,reshape=None):\n",
    "    if method=='SVD':\n",
    "        return svd_compress_ratio(image,rank)\n",
    "    if method=='CP':\n",
    "        return cp_compress_ratio(image,rank)\n",
    "    if method=='TT':\n",
    "        return tt_compress_ratio(image,rank)\n",
    "    if method=='TUCKER':\n",
    "        return tucker_compress_ratio(image,rank)\n",
    "    if method=='TT_reshape':\n",
    "        return tt_reshape_compress_ratio(image,rank,reshape[0],reshape[1])\n",
    "    \n",
    "def compress(image,method,rank,reshape=None):\n",
    "    if method=='SVD':\n",
    "        return svd_compression(image,rank)\n",
    "    if method=='CP':\n",
    "        return cp_compression(image,rank)\n",
    "    if method=='TT':\n",
    "        return tt_compression(image,rank)\n",
    "    if method=='TUCKER':\n",
    "        return tucker_compression(image,rank)\n",
    "    if method=='TT_reshape':\n",
    "        return tt_reshape_compression(image,rank,reshape[0],reshape[1])\n",
    "    \n",
    "def approx_ratio_to_rank(image,method,target_ratio,reshape=None):\n",
    "    rank = 1\n",
    "    while True:\n",
    "        ratio = compress_ratio(image,method,rank,reshape)\n",
    "        if ratio > target_ratio:\n",
    "            return rank-1\n",
    "        rank += 1\n",
    "        \n",
    "\n",
    "def plot_compressed(image,method,rank,reshape=None):\n",
    "    compressed_img = compress(image,method,rank,reshape)\n",
    "    plt.imshow(compressed_img)\n",
    "    plt.title(\"%s (rank: %i - compress ratio: %.4f)\" % (method,rank,compress_ratio(image,method,rank,reshape)))\n",
    "    return compressed_img\n",
    "              \n",
    "\n",
    "plt.figure(figsize=(15,30))\n",
    "\n",
    "reshape=([2,2,2,2,2,2,2,2],[2,2,2,2,2,2,2,2])\n",
    "for i,method in enumerate([\"SVD\",\"CP\",\"TUCKER\",\"TT\",\"TT_reshape\"]):\n",
    "    plt.subplot(5,2,2*i+1)\n",
    "    rank = approx_ratio_to_rank(image1,method,ratio,reshape)\n",
    "    plot_compressed(image1,method,rank,reshape)\n",
    "\n",
    "for i,method in enumerate([\"SVD\",\"CP\",\"TUCKER\",\"TT\",\"TT_reshape\"]):\n",
    "    plt.subplot(5,2,2*i+2)\n",
    "    rank = approx_ratio_to_rank(image2,method,ratio,reshape)\n",
    "    plot_compressed(image2,method,rank,reshape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Write your observations in this cell (if any)*\n",
    "\n",
    "... ... ..."
   ]
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}