low-rank.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Low-Rank Approximation\n",
    "\n",
    "---\n",
    "\n",
    "> Author: <font color='#f78c40'>Samuel Farrens</font>    \n",
    "> Year: 2017  \n",
    "> Email: [samuel.farrens@gmail.com](mailto:samuel.farrens@gmail.com)  \n",
    "> Website: <a href=\"https://sfarrens.github.io\" target=\"_blank\">https://sfarrens.github.io</a>\n",
    "\n",
    "---\n",
    "\n",
    "\n",
    "## Contents\n",
    "\n",
    "---\n",
    " \n",
    "1. [Set-Up](#Set-Up)\n",
    "1. [Introduction](#Introduction)\n",
    "1. [Rank of a Matrix](#Rank-of-a-Matrix)\n",
    " * [Full Rank Matrix](#Full-Rank-Matrix) \n",
    " * [Reduced Rank Matrix](#Reduced-Rank-Matrix) \n",
    " * [Uneven Matrix](#Uneven-Matrix)\n",
    "1. [Singular Value Decomposition](#Singular-Value-Decomposition) \n",
    " * [Reducing the Rank](#Reducing-the-Rank)\n",
    " * [Hard Thresholding](#Hard-Thresholding)\n",
    "1. [Low Rank Approximation](#Low-Rank-Approximation)\n",
    " * [Nuclear Norm](#Nuclear-Norm)\n",
    "1. [Deconvolution Exercise](#Deconvolution-Exercise)\n",
    " * [Gradient](#Gradient)\n",
    " * [Rank Reduction](#Rank-Reduction)\n",
    " * [Optimisation](#Optimisation)\n",
    " * [Deconvolution](#Deconvolution)\n",
    "1. [Hints](#Hints)\n",
    "1. [Exercise Solutions](#Exercise-Solutions)\n",
    " * [Deconvolution Exercise Solution](#Deconvolution-Exercise-Solution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Set-Up\n",
    "\n",
    "Here we will import a couple of packages that we will need throughout the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tell Jupyter to display plots in this notebook.\n",
    "%matplotlib inline\n",
    "\n",
    "# Import the numpy package with the alias np.\n",
    "import numpy as np        \n",
    "\n",
    "# Import the FFT convolution function from Scipy.\n",
    "from scipy.signal import fftconvolve\n",
    "\n",
    "# Import the Gaussian kernel function from Astropy.\n",
    "from astropy.convolution import Gaussian2DKernel\n",
    "\n",
    "# Import the pyplot package from matplotlib with the alias plt.\n",
    "import matplotlib.pyplot as plt  \n",
    "from matplotlib import pylab\n",
    "\n",
    "# Ignore warnings.\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Introduction\n",
    "\n",
    "The knowledge that a given matrix contains a certain amount of redundancy can be used a powerful prior for solving inverse problems. The notebook provides a brief introduction to the concepts of matrix rank and the low-rank approximation for regularisation. This notebook does not provide a comprehensive review on these topics, rather it is intended to give the user some intuitive hands-on experience."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Rank of a Matrix\n",
    "\n",
    "The rank of a matrix can be defined in the following ways:\n",
    "\n",
    "1. the maximum number of linearly independent column vectors in a given matrix\n",
    "1. the maximum number of linearly independent row vectors in a given matrix\n",
    "\n",
    "Both of these definitions are equivalent.\n",
    "\n",
    "### <font color='blue'>Full Rank Matrix</font>\n",
    "\n",
    "For the matrix\n",
    "\n",
    "$$ A = \\begin{bmatrix} 6 & 1 & 2 \\\\ 5 & 6 & 1 \\\\ 5 & 7 & 7 \\end{bmatrix} $$\n",
    "\n",
    "the rank is given by\n",
    "\n",
    "$$\\textrm{rank}(A) = 3$$\n",
    "\n",
    "since all three rows are linearly idependent. $A$ is then said to have **full rank**.\n",
    "\n",
    "This can be implemented in Python as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[6, 1, 2], [5, 6, 1], [5, 7, 7]])\n",
    "\n",
    "print('A =')\n",
    "print(A)\n",
    "print('')\n",
    "print('rank(A) =', np.linalg.matrix_rank(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color='blue'>Reduced Rank Matrix</font>\n",
    "\n",
    "For the matrix\n",
    "\n",
    "$$ A = \\begin{bmatrix} 1 & 0 & 1 \\\\ -2 & -3 & 1 \\\\ 3 & 3 & 0 \\end{bmatrix} $$\n",
    "\n",
    "the rank is given by\n",
    "\n",
    "$$\\textrm{rank}(A) = 2$$\n",
    "\n",
    "since the third row is simply the difference of the first two, which are linearly idependent. In this case $A$ is **rank deficient** or has **reduced rank**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[1, 0, 1], [-2, -3, 1], [3, 3, 0]])\n",
    "\n",
    "print('A =')\n",
    "print(A)\n",
    "print('')\n",
    "print('rank(A) =', np.linalg.matrix_rank(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color='blue'>Uneven Matrix</font>\n",
    "\n",
    "For a matrix $\\in\\mathbb{R}^{m\\times n}$, if $m>n$ then the maximum rank of the matrix is $n$ (*i.e.* if the matrix has more rows than columns then the maximum value of the rank is the number of columns). Conversely, if $m<n$ then the maximum rank is $m$.\n",
    "\n",
    "For example, we can calculate the rank of a $5\\times 3$ matrix $A$ and a $3\\times 2$ matrix $B$ both of which have linearly independent columns as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "A = np.random.randint(-9, 9, (5, 3))\n",
    "B = np.random.randint(-9, 9, (3, 2))\n",
    "\n",
    "print('A =')\n",
    "print(A)\n",
    "print('')\n",
    "print('rank(A) =', np.linalg.matrix_rank(A))\n",
    "print('')\n",
    "print('B =')\n",
    "print(B)\n",
    "print('')\n",
    "print('rank(B) =', np.linalg.matrix_rank(B))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Singular Value Decomposition\n",
    "\n",
    "The rank of a matrix can also be measured by counting the number of non-zero singular values.\n",
    "\n",
    "The <a href=\"https://en.wikipedia.org/wiki/Singular_value_decomposition\" target=\"_blank\">SVD</a> of a matrix $A$ is given by\n",
    "\n",
    "$$A = U\\Sigma V^T$$\n",
    "\n",
    "where $\\Sigma$ are the singular values. Thus, the rank of $A$ is\n",
    "\n",
    "$$\\textrm{rank}(A) = \\|\\Sigma\\|_0$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the matrix A\n",
    "A = np.array([[1, 0, 1], [-2, -3, 1], [3, 3, 0]])\n",
    "\n",
    "# Calculate the SVD of A\n",
    "U, S, V = np.linalg.svd(A)\n",
    "\n",
    "# Round the values of S to two decimal places. Comment out this operation to see why this is done.\n",
    "S = np.around(S, 2)\n",
    "\n",
    "print('A =')\n",
    "print(A)\n",
    "print('')\n",
    "print('Sigma =')\n",
    "print(np.diag(S))\n",
    "print('')\n",
    "print('rank(A) =', sum(S > 0))\n",
    "print('(or rank(A) =', np.linalg.norm(S, 0), ')')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color='blue'>Reducing the Rank</font>\n",
    "\n",
    "By definition setting one or more singular values to zero reduces the rank of a given matrix.\n",
    "\n",
    "For example, for a matrix with rank 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the matrix A\n",
    "A = np.array([[7, 3, -4], [-5, 1, 2], [2, 4, -2]])\n",
    "\n",
    "# Calculate the SVD of A\n",
    "U, S, V = np.linalg.svd(A)\n",
    "S = np.around(S, 2)\n",
    "\n",
    "print('A =')\n",
    "print(A)\n",
    "print('')\n",
    "print('Sigma =')\n",
    "print(np.diag(S))\n",
    "print('')\n",
    "print('rank(A) =', np.linalg.matrix_rank(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "the rank is reduced to 1 when the second singular value is set to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Set the second singular value to 0\n",
    "S[1] = 0\n",
    "\n",
    "# Create a new matrix with the updated singular values\n",
    "B = np.dot(U, np.dot(np.diag(S), V))\n",
    "\n",
    "print('Sigma =')\n",
    "print(np.diag(S))\n",
    "print('')\n",
    "print('B =')\n",
    "print(B)\n",
    "print('')\n",
    "print('rank(B) =', np.linalg.matrix_rank(B))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color='blue'>Hard Thresholding</font>\n",
    "\n",
    "Another way of looking at this would be apply a **Hard Threshold** to the singular values of the matrix.\n",
    "\n",
    "$$ \n",
    "HT_{\\lambda}(\\mathbf{x}) = \n",
    "\\begin{cases} \n",
    "\\mathbf{x} & \\text{if}\\ |\\mathbf{x}| \\geq \\lambda \\\\ \n",
    "0 & \\text{otherwise} \\\\\n",
    "\\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hard_thresh(value, threshold):\n",
    "    '''Hard Threshold\n",
    "    \n",
    "    Returns hard threshold of input\n",
    "    \n",
    "    '''\n",
    "    \n",
    "    return value * (np.abs(value) >= threshold)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, for a full rank matrix $A$ ($\\textrm{rank}(A)=3$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a matrix of random values with seed 1\n",
    "np.random.seed(1)\n",
    "A = np.random.randint(-9, 9, (3, 3))\n",
    "\n",
    "# Calculate the SVD of A\n",
    "U, S, V = np.linalg.svd(A)\n",
    "S = np.around(S, 2)\n",
    "\n",
    "print('A =')\n",
    "print(A)\n",
    "print('')\n",
    "print('Sigma =')\n",
    "print(np.diag(S))\n",
    "print('')\n",
    "print('rank(A) =', np.linalg.matrix_rank(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "a new matrix $B$ with a rank of 1 can be obtained by applying a hard threshold with $\\lambda=6$ to the singular values of $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hard threshold values of sigma\n",
    "S_ht = hard_thresh(S, 6)\n",
    "\n",
    "# Create a new matrix with the updated singular values\n",
    "B = np.dot(U, np.dot(np.diag(S_ht), V))\n",
    "\n",
    "print('HT(Sigma) =')\n",
    "print(np.diag(S_ht))\n",
    "print('')\n",
    "print('B =')\n",
    "print(B)\n",
    "print('')\n",
    "print('rank(B) =', np.linalg.matrix_rank(B))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, with an appropriate threshold, one can reduce the rank of a given matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Low Rank Approximation\n",
    "\n",
    "The <a href=\"https://en.wikipedia.org/wiki/Low-rank_approximation\" target=\"_blank\">low-rank approximation</a> is a minimisation constraint that is subject to the condition that the matrix one wishes to obtain has reduced rank. In other words, we can use the prior knowledge that the matrix one wants to recover is of reduced rank when solving an inverse problem.\n",
    "\n",
    "### <font color='blue'>Nuclear Norm</font>\n",
    "\n",
    "As with the $l_0$ norm in sparse regularisation, the rank of the matrix cannot be used directly for regularisation as it is a non-convex function. Instead the nuclear norm (or trace norm) is used to promote low-rank solutions. The nuclear norm is given by the sum of singular values,\n",
    "\n",
    "$$\\|X\\|_* = \\sum_{k=1} \\sigma_k(X)$$\n",
    "\n",
    "where $\\sigma_k(X)$ denotes the $k^{\\textrm{th}}$ largest singular value of $X$. This can be implemented in Python as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def nuclear_norm(data):\n",
    "    '''Nuclear Norm\n",
    "    \n",
    "    This function returns the nuclear norm of the input data.\n",
    "    \n",
    "    '''\n",
    "\n",
    "    # Get SVD of the data.\n",
    "    u, s, v = np.linalg.svd(data)\n",
    "\n",
    "    # Return nuclear norm.\n",
    "    return np.sum(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, given a full rank matrix $A$ and a rank 2 matrix $B$ one can see that the nuclear norm of $B$ is smaller than that of $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[6, 1, 2], [5, 6, 1], [5, 7, 7]])\n",
    "B = np.array([[1, 0, 1], [-2, -3, 1], [3, 3, 0]])\n",
    "\n",
    "print('rank(A)=', np.linalg.matrix_rank(A))\n",
    "print('||A||*=', nuclear_norm(A))\n",
    "print('')\n",
    "print('rank(B)=', np.linalg.matrix_rank(B))\n",
    "print('||B||*=', nuclear_norm(B))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## <font color='red'>Deconvolution Exercise</font>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This exercise deals with the concept of deconvolution.\n",
    "\n",
    "We begin by defining a matrix $X$ with rank 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#############################\n",
    "# NO NEED TO EDIT THIS CELL #\n",
    "# JUST EXECUTE IT           #\n",
    "#############################\n",
    "\n",
    "# Define a 11x11 matrix X\n",
    "vals = np.array(list(range(6)) + list(range(5))[::-1]) \n",
    "X = np.array([i * vals for i in range(6)] + [i * vals for i in range(5)[::-1]], dtype='float')\n",
    "X /= X.sum()\n",
    "\n",
    "# Calculate the rank of X\n",
    "print('rank(X)=', np.linalg.matrix_rank(X))\n",
    "\n",
    "# Display X\n",
    "plt.imshow(X, interpolation='none', cmap='magma')\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A slightly blurred version of $X$ with noise is given by\n",
    "\n",
    "$$Y = MX + N$$\n",
    "\n",
    "where $M$ represents the convolution of $X$ with a blurring kernel and $N$ is additive Gaussian noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#############################\n",
    "# NO NEED TO EDIT THIS CELL #\n",
    "# JUST EXECUTE IT           #\n",
    "#############################\n",
    "\n",
    "def blur(data, kernel):\n",
    "    '''Blurring operator\n",
    "\n",
    "    This function convolves the input data with the blurring kernel.\n",
    "\n",
    "    '''\n",
    "    \n",
    "    return fftconvolve(data, kernel, 'same')\n",
    "\n",
    "def blur_transpose(data, kernel):\n",
    "    '''Adjoint blurring operator\n",
    "    \n",
    "    This function convolves the input data with the rotated blurring kernel.\n",
    "    \n",
    "    '''\n",
    "    \n",
    "    return fftconvolve(data, np.rot90(kernel, 2), 'same')\n",
    "\n",
    "# Define a 3x3 Gaussian kernel\n",
    "gauss_kernel = np.array(Gaussian2DKernel(0.5, x_size=3, y_size=3))\n",
    "\n",
    "# Blur X and add noise\n",
    "Y = blur(X, gauss_kernel) + 0.003 * np.random.randn(*X.shape)\n",
    "\n",
    "# Calculate the rank of Y\n",
    "print('rank(Y)=', np.linalg.matrix_rank(Y))\n",
    "\n",
    "# Display X and Y\n",
    "plt.figure(figsize=(15, 12))\n",
    "plt.subplot(121)\n",
    "plt.imshow(X, interpolation='none', cmap='magma')\n",
    "plt.title('X', fontsize=20)\n",
    "plt.subplot(122)\n",
    "plt.imshow(Y, interpolation='none', cmap='magma')\n",
    "plt.title('Y', fontsize=20)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$Y$ is clearly distorted with respect to $X$. This sort of distortion is present in virtually all images obtained with optical instruments (*e.g.* telescopes, microscopes, *etc.*).\n",
    "\n",
    "In order to recover $X$ we will aim to solve the problem using the following minimsation:\n",
    "\n",
    "$$\\begin{aligned} & \\underset{X}{\\text{argmin}} & \\frac{1}{2}\\|Y-MX\\|_2^2 + \\lambda\\|X\\|_*\\end{aligned}$$\n",
    "\n",
    "which imposes a low-rank solution of $X$ via the nuclear norm term.\n",
    "\n",
    "### <font color='red'>Gradient</font>\n",
    "\n",
    "Your first task is to define a function that calculates the gradient of the minimisation problem, which is given by\n",
    "\n",
    "$$\\nabla F(X) = M^T(MX-Y)$$\n",
    "\n",
    "> See the [Hints](#Hints) if you get stuck.\n",
    "\n",
    "<br />\n",
    "<font color='red'>EDIT THE CELL BELOW</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##############################\n",
    "# YOU NEED TO EDIT THIS CELL #\n",
    "##############################\n",
    "\n",
    "def grad(y, x_rec, kernel):\n",
    "    '''Gradient\n",
    "    \n",
    "    This function calculates the gradient given the observation and an approximation of the true signal.\n",
    "    \n",
    "    '''\n",
    "\n",
    "    return None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color='red'>Rank Reduction</font>\n",
    "\n",
    "Your next task is to define a function for imposing that the recovered matrix has a reduced rank.\n",
    "\n",
    "To complete this task you will need to:\n",
    "1. calculate some measure of the rank of the input data\n",
    "1. reduce the this rank using the input threshold\n",
    "1. return the updated matrix\n",
    "\n",
    "> See the [Hints](#Hints) if you get stuck.\n",
    "\n",
    "<br />\n",
    "<font color='red'>EDIT THE CELL BELOW</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##############################\n",
    "# YOU NEED TO EDIT THIS CELL #\n",
    "##############################\n",
    "\n",
    "def reduce_rank(x, threshold):\n",
    "    '''Reduce Rank\n",
    "    \n",
    "    This function reduces the rank of the input value x using the specified threshold.\n",
    "    \n",
    "    '''\n",
    "    \n",
    "    return None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  <font color='blue'>Optimisation</font>\n",
    "\n",
    "For this exercise we will use the forward-backward splitting algorithm with a rank reduction operation to promote low-rank solutions. It has the following form for each iteration:\n",
    "\n",
    "$$X_{n+1} = X_n - \\gamma * \\nabla F(X_n)$$\n",
    "$$X_{n+1} = \\text{RR}_\\lambda(X_{n+1})$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#############################\n",
    "# NO NEED TO EDIT THIS CELL #\n",
    "# JUST EXECUTE IT           #\n",
    "#############################\n",
    "\n",
    "def cost_func(y, x_rec, threshold, kernel):\n",
    "    '''Cost Function\n",
    "    \n",
    "    Return the cost of the given estimate of x\n",
    "    \n",
    "    '''\n",
    "    \n",
    "    return 0.5 * np.linalg.norm(y - blur(x_rec, kernel)) ** 2 + threshold * nuclear_norm(x_rec)\n",
    "\n",
    "def forwardBackward(observation, first_guess, threshold, kernel):\n",
    "    '''Forward-Backward Algorithm\n",
    "    \n",
    "    Return the recovered signal using a rank reduction operator\n",
    "    \n",
    "    '''\n",
    "            \n",
    "    x_rec = first_guess\n",
    "    \n",
    "    n_iter = 15\n",
    "    gamma = 1.0\n",
    "        \n",
    "    for i in range(n_iter):\n",
    "        \n",
    "        x_temp = x_rec - gamma * grad(observation, x_rec, kernel)\n",
    "        x_rec = reduce_rank(x_temp, threshold)\n",
    "        cost = cost_func(observation, x_rec, threshold, kernel)\n",
    "        print('Iteration', i + 1, 'Cost =', cost)\n",
    "        \n",
    "    return x_rec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color='red'>Deconvolution</font>\n",
    "\n",
    "Your final task is to perform the deconvolution.\n",
    "\n",
    "1. Set the low-rank threshold value ($\\lambda$).\n",
    "1. Set an initial value for $X$.\n",
    "1. Implement the forward-backward algorithm with the following inputs:\n",
    " 1. The blurred noisy matrix.\n",
    " 1. The initial value for $X$.\n",
    " 1. The threshold.\n",
    " 1. The blurring kernel.\n",
    "\n",
    "> See the [Hints](#Hints) if you get stuck.\n",
    "\n",
    "<br />\n",
    "<font color='red'>EDIT THE CELL BELOW</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##############################\n",
    "# YOU NEED TO EDIT THIS CELL #\n",
    "##############################\n",
    "\n",
    "# Set the threshold value.\n",
    "threshold = None\n",
    "\n",
    "# Set initial guess for X.\n",
    "first_guess = None\n",
    "\n",
    "# Recover the matrix.\n",
    "X_rec = None\n",
    "\n",
    "# Calculate the rank of recoverd matrix.\n",
    "rank_x_rec = None\n",
    "\n",
    "if not isinstance(rank_x_rec, type(None)):\n",
    "    print('')\n",
    "    print('rank(X_rec)=', rank_x_rec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#############################\n",
    "# NO NEED TO EDIT THIS CELL #\n",
    "# JUST EXECUTE IT           #\n",
    "#############################\n",
    "\n",
    "# Display\n",
    "if not isinstance(X_rec, type(None)):\n",
    "    plt.figure(figsize=(17, 12))\n",
    "    plt.subplot(131)\n",
    "    plt.imshow(X, interpolation='none', cmap='magma')\n",
    "    plt.title('$X$', fontsize=20)\n",
    "    plt.subplot(132)\n",
    "    plt.imshow(Y, interpolation='none', cmap='magma')\n",
    "    plt.title('$Y$', fontsize=20)\n",
    "    plt.subplot(133)\n",
    "    plt.imshow(X_rec, interpolation='none', cmap='magma')\n",
    "    plt.title('$\\hat{X}$', fontsize=20)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Does your final $\\hat{X}$ look more like $X$ or $Y$? Try ajusting $\\lambda$ to see if that improves your results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## <font color='orange'>Hints</font>\n",
    "\n",
    "** Gradient **\n",
    "\n",
    "1. The blurring operators are already defined.\n",
    "\n",
    "** Rank Reduction **\n",
    "\n",
    "1. Make sure you are thresholding something *singular*.\n",
    "1. Thresholding matrices can be quite *hard*.\n",
    "\n",
    "** Deconvolution **\n",
    "\n",
    "1. Make sure your first guess of $X$ is the right shape.\n",
    "1. Try adjusting your value of $lambda$\n",
    "\n",
    "---\n",
    "\n",
    "## Exercise Solutions\n",
    "\n",
    "**<font color=\"red\">Warning!</font>** Try your best to solve all of the exercises on your own before checking the answers.\n",
    "\n",
    "<br />\n",
    "\n",
    "### Deconvolution Exercise Solution\n",
    "\n",
    "** Gradient **\n",
    "\n",
    "```Python\n",
    "\n",
    "def grad(y, x_rec, kernel):\n",
    "    '''Gradient\n",
    "    \n",
    "    This function calculates the gradient given the observation and an approximation of the true signal.\n",
    "    \n",
    "    '''\n",
    "\n",
    "    return blur_transpose(blur(x_rec, kernel) - y, kernel)\n",
    "    \n",
    "```\n",
    "\n",
    "** Rank Reduction **\n",
    "\n",
    "```Python\n",
    "\n",
    "def reduce_rank(x, threshold):\n",
    "    '''Reduce Rank\n",
    "    \n",
    "    This function reduces the rank of the input value x using the specified threshold.\n",
    "    \n",
    "    '''\n",
    "    \n",
    "    u, s, v = np.linalg.svd(x)\n",
    "        \n",
    "    s_new = hard_thresh(s, threshold)\n",
    "    \n",
    "    return np.dot(u, np.dot(np.diag(s_new), v))\n",
    "\n",
    "```\n",
    "\n",
    "** Deconvolution **\n",
    "\n",
    "```Python\n",
    "\n",
    "# Set the threshold value.\n",
    "threshold = 0.1\n",
    "\n",
    "# Set initial guess for X.\n",
    "first_guess = np.ones(Y.shape)\n",
    "\n",
    "# Recover the matrix.\n",
    "X_rec = forwardBackward(Y, first_guess, threshold, gauss_kernel)\n",
    "\n",
    "# Calculate the rank of recoverd matrix.\n",
    "rank_x_rec = np.linalg.matrix_rank(X_rec)\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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}