Finite dimensional inverse problems A geometrical perspective.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finite dimensional linear inverse problems: A geometrical perspective\n",
    "## Exercise Session 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Created by Adam Scholefield, Michalina Pacholska, Ivan Dokmanic and Gilles Baechler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we will play around with some finite dimensional linear inverse problems.\n",
    "Lets first import all the packages we will need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# %pylab inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LogNorm\n",
    "from scipy import signal\n",
    "from scipy.fftpack import dct\n",
    "from skimage import data_dir\n",
    "from skimage.transform import radon, rescale\n",
    "from skimage.io import imread\n",
    "from tqdm import tqdm_notebook as tqdm\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1-D example: ECG signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read and display the ECG signal\n",
    "data_ecg = np.genfromtxt('ecg.dat', delimiter=',')\n",
    "x = data_ecg[:-1,1] #subtract one element to make N even\n",
    "t = data_ecg[:-1,0]\n",
    "N = len(x)\n",
    "plt.plot(t, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a box function\n",
    "box = np.zeros(N)\n",
    "box[int(N/2)-20:int(N/2)+21] = 1/41\n",
    "plt.plot(box)\n",
    "plt.title('Box function')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Square convolution matrix\n",
    "We are going to convolve the ECG signal ($x$) with box function ($s$). We could just use a built in convolve funtion but we want to write this as the matrix vector product $Ax$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convolve the ECG signal with the box.\n",
    "# We're going to creeate a matrix to do this\n",
    "A = np.zeros((N, N))\n",
    "for i in range(N):\n",
    "    A[i:]=np.roll(box, int(i-N/2)) #assuming N is even\n",
    "\n",
    "print(A.shape)\n",
    "plt.figure(figsize=(16,4))\n",
    "plt.imshow(A+0.001, interpolation='None', norm=LogNorm(vmin=0.001, vmax=10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task. calculate the convolution as $y=Ax$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "y = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(y)\n",
    "plt.title('Blurred signal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task. Use the inverse of $A$ to estimate $x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# TODO\n",
    "x_hat = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets plot the result and check that it's close to $x$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,2,figsize = (15,5))\n",
    "ax[0].plot(t,x)\n",
    "ax[1].plot(t,x_hat)\n",
    "ax[0].set_title('Original (x)')\n",
    "ax[1].set_title('Esitmate (x_hat)')\n",
    "plt.show()\n",
    "print('Mean squared error (MSE) is', np.mean((x-x_hat)**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spectrum of the operator\n",
    "We just blindly calculated the inverse. Whilst this worked with this blurring kernel and no noise, in general it's not a good strategy. To see this lets add a small amount of noise to $y$ and try to reconstruct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add some noise to y\n",
    "M = len(y)\n",
    "y_noisy = y+2e-3*np.random.randn(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task. Try to estimate $x$ using the inverse of $A$. Plot the result and calculate the MSE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "x_hat = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that something has definetly gone wrong. \n",
    "To understand what has happened, lets plot the singular values of $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the signular values of A\n",
    "U, s, Vh = np.linalg.svd(A)\n",
    "plt.plot(s)\n",
    "plt.title('Singular values of A')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The singular values decay quickly but stay large enough to allow invertability in the noiseless case. However, with a small amount of noise, it breaks down.\n",
    "\n",
    "Trimming the SVD provides a more stable way to invert into the range space of the trimmed operator. \n",
    "We know that $\\boldsymbol{A}=\\boldsymbol{U\\Sigma V}^\\top$. However, since some of the singular values are too small, lets remove them from $\\boldsymbol{\\Sigma}$. This will produce a new trimmed matrix, which we will call $\\boldsymbol{\\Sigma_T}$. In addition, we will need to trim the corresponding columns of $\\boldsymbol{U}$ and $\\boldsymbol{V}$.\n",
    "### Task. Choose a sensible place to trim the singular values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#TODO\n",
    "trim_len = \n",
    "plt.title('Singular values of A and cut off point')\n",
    "plt.plot(s)\n",
    "plt.axvline(x=trim_len)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now trim the vector $\\boldsymbol{s}$ and $\\boldsymbol{U}$ and $\\boldsymbol{V}$ matrices, producing $\\boldsymbol{s_T}$ and $\\boldsymbol{U_T}$ and $\\boldsymbol{V_T}$. If we choose a good place to trim, we should have $\\boldsymbol{A}\\simeq\\boldsymbol{U_T\\Sigma_T V_T}^\\top$.\n",
    "Note that the SVD function of numpy returns the transpose of V, which we call Vh (for V Hermitian)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_trimmed = s[:trim_len]\n",
    "U_trimmed = U[:,:trim_len]\n",
    "Vh_trimmed = Vh[:trim_len,:]\n",
    "print(U_trimmed.shape)\n",
    "print(Vh_trimmed.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets verify that $\\boldsymbol{A}\\simeq\\boldsymbol{U_T\\Sigma_T V_T}^\\top$. We will print the maximum difference of an entry (feel free to compare them in a different way)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.max(np.abs(A-U_trimmed@np.diag(s_trimmed)@Vh_trimmed)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the SVD, the previous estimator $\\hat{x}=A^T(AA^T)^{-1}y$ can be written as $\\hat{x}=V\\Sigma^{-1}U^T y$.\n",
    "\n",
    "### Task. Estimate x by using the trimmed matrices to produce an approximation of this projection. Your solution should be in the range of $\\boldsymbol{V_T}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "x_hat = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets plot the result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,2,figsize = (15,5))\n",
    "ax[0].plot(t,x)\n",
    "ax[1].plot(t,x_hat)\n",
    "ax[0].set_title('Original (x)')\n",
    "ax[1].set_title('Esitmate (x_hat)')\n",
    "plt.show()\n",
    "print('Mean squared error (MSE) is', np.mean((x-x_hat)**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Much better.\n",
    "\n",
    "Often, we also have to do this in the noiseless case. After the session, you can experiment with a Gaussian blur forward operator, such as below. \n",
    "\n",
    "**For now, just run the cell below to see the kernel. If you complete everything else, you can come back and play around with this kernel**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gaussian = signal.gaussian(N+1, std=10)\n",
    "gaussian = gaussian[:N]\n",
    "A_gaussian = np.zeros((N, N))\n",
    "for i in range(N):\n",
    "    A_gaussian[i:]=np.roll(gaussian, int(i-N/2))\n",
    "\n",
    "U, s, Vh = np.linalg.svd(A_gaussian)\n",
    "fig, ax = plt.subplots(1,3,figsize = (15,5))\n",
    "ax[0].plot(gaussian)\n",
    "ax[1].imshow(A_gaussian+0.001, interpolation='None', norm=LogNorm(vmin=0.001, vmax=10))\n",
    "ax[2].plot(s)\n",
    "ax[0].set_title('Gaussian kernel')\n",
    "ax[1].set_title('Visualisation of Gaussian blur matrix')\n",
    "ax[2].set_title('Singular values of Gaussian blur matrix')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the rest of this section using the ECG signal, we will use the box function and no noise. Therefore, we do not need to worry about truncating the small singular values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fat convolution matrix (downsampling)\n",
    "We will now remove half the rows of A so that we have a downsampling matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove half the rows of A\n",
    "down_sample_factor = 2\n",
    "A2 = A[::down_sample_factor,:]\n",
    "t2 = t[::down_sample_factor]\n",
    "print(A2.shape)\n",
    "plt.figure(figsize=(16,4))\n",
    "plt.imshow(A2+0.001, interpolation='None', norm=LogNorm(vmin=0.001, vmax=10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the downsampled measurements\n",
    "y2 = A2@x\n",
    "plt.figure()\n",
    "plt.plot(t2,y2)\n",
    "plt.title('Blurred and downsampled signal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task. Find the consistent estimate with minimum $\\\\ell_2$-norm; i.e., find the consistent estimate such that the analysis and syntheis operators are ideally matched."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "x_hat = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets plot the result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,2,figsize = (15,5))\n",
    "ax[0].plot(t,x)\n",
    "ax[1].plot(t,x_hat)\n",
    "ax[0].set_title('Original')\n",
    "ax[1].set_title('Orthogonal projection onto R(A2^T) - Consistent and ideally matched')\n",
    "plt.show()\n",
    "print('MSE is', np.mean((x-x_hat)**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task. Verify that your estimator is consistent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Oblique projection\n",
    "Suppose we believe that an ECG signal does not contain very high frequencies. We can use this prior to project onto the range space of a different operator, whilst maintaining consistency. Let's use the discrete cosine transform (DCT) matrix and keep the $M$ lowest frequencies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# construct matrix B so R(B) is M lowest frequencies \n",
    "[M, N] = A2.shape\n",
    "B = dct(np.eye(N)) #create an NxN DCT matrix\n",
    "B = B[:,:M] #remove all cols after freq_to_keep to only keep the freq_to_keep lowest frequencies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task. Estimate $x$ by projecting onto the range of $B$ whilst maintaining consistency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "x_hat = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,2,figsize = (15,5))\n",
    "ax[0].plot(t,x)\n",
    "ax[1].plot(t,np.real(x_hat))\n",
    "ax[0].set_title('Original')\n",
    "ax[1].set_title('Oblique projection onto R(B) - Consistent')\n",
    "plt.show()\n",
    "print('MSE is', np.mean((x-x_hat)**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the prior has helped recover the shape but we are including too high frequencies. We chose to use the $M$ lowest frequency terms so that the matrix $AB$ was square and could be inverted. However, a better prior would be to use less frequency terms. While enforcing this prior is more complicated, it really shows the merit of the geometrical way of thinking! Since I couldn't resist including this, it appears below. It is beyond what is expected for this first session but, if you're interested, ask an instructor to explain the geometry. It's beautiful :-)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optional. Ask an instructor about this (you may want to complete the tomography part before looking at this)\n",
    "Now we'll decrease the number of frequency components we use so that we reduce the dimension of the range space of $B$. In this case, the affine subspace of consistent solutions and the range space of $B$ will not intersect, in general. Therefore, we can either find the solution in the range space of $B$ that is closest to the affine subspace of consistent solutions or find the solution in the affine subspace of consistent solutions that is closest to the range space of $B$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# we will find both the solution in R(B) that is closest to the affine subspace of consistent solutions\n",
    "# and the solution in the affine subspace of consistent solutions that is closest to R(B)\n",
    "\n",
    "freq_to_keep = 120 #half the number of frequency terms we will use\n",
    "B = dct(np.eye(N)) #create an NxN DCT matrix\n",
    "B = B[:,:freq_to_keep] #remove all cols after freq_to_keep to only keep the freq_to_keep lowest frequencies\n",
    "\n",
    "U, s, Vh = np.linalg.svd(A2) #take the SVD of A2 so that we can abstract a bases for its null space\n",
    "basesNullspaceA = Vh[len(s):,:].T #abstract the null space\n",
    "T = np.hstack([B,-basesNullspaceA]) #concatonate a bases for B with a bases for the null space of A\n",
    "\n",
    "coeffs = np.linalg.inv(T.T@T)@T.T@A2.T@np.linalg.inv(A2@A2.T)@y2 #solve the least squares problem (first 2*half_len coeffs are for B and the rest for the null space of A)\n",
    "x_hat = B@coeffs[:freq_to_keep] # point in R(B) that is closest to affine subspace of consistent solutions\n",
    "x_hat2 = basesNullspaceA@coeffs[freq_to_keep:] + A2.T@np.linalg.inv(A2@A2.T)@y2 #consistent solution closest to R(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets plot both estimators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(3,2,figsize = (20,15))\n",
    "ax[0][0].plot(t,x)\n",
    "ax[0][1].plot(t2,y2)\n",
    "ax[0][0].set_title('Original (x)')\n",
    "ax[0][1].set_title('Downsampled and blurred (y)')\n",
    "ax[1][0].plot(t,x_hat)\n",
    "ax[1][1].plot(t,np.real(x_hat2))\n",
    "ax[1][0].set_title('Point in R(B) that is closest to affine subspace of consistent solutions')\n",
    "ax[1][1].set_title('Point in affine subspace of consistent solutions that is closest to R(B)')\n",
    "ax[2][0].plot(t,np.real(x_hat)-x)\n",
    "ax[2][1].plot(t,np.real(x_hat2)-x)\n",
    "ax[2][0].set_title('Error for point in R(B) that is closest to affine subspace of consistent solutions')\n",
    "ax[2][1].set_title('Error for point in affine subspace of consistent solutions that is closest to R(B)')\n",
    "plt.show()\n",
    "print('MSE for point in R(B) that is closest to affine subspace of consistent solutions is', np.mean((x-x_hat)**2))\n",
    "print('MSE for point in affine subspace of consistent solutions that is closest to R(B) is', np.mean((x-x_hat2)**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want, you can verify if the estimators are consistent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tomography example (in 2D)\n",
    "Let's now look at a 2D version of tomography."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's first load an image\n",
    "N = 64\n",
    "image = imread(data_dir + \"/phantom.png\", as_gray=True)\n",
    "image = rescale(image, (N / image.shape[0], N / image.shape[1]), mode='constant', multichannel = False)\n",
    "\n",
    "plt.imshow(image, interpolation='None')\n",
    "plt.gray()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As explained in class, in X-ray tomography, x-rays are fired through the object at different angles and the transmission is measured at the other side. To simulate these measurements, we want to be able to compute integrals at different angles. For example, it is very easy to do this horizontally and vertically by just summing the pixels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets sum the columns to give the projection for x-rays fired vertically\n",
    "# and sum the rows to give the projection for x-rays fired horizontally\n",
    "fig, ax = plt.subplots(1,2,figsize = (20,5))\n",
    "ax[0].plot(np.sum(image,0))\n",
    "ax[1].plot(np.sum(image,1))\n",
    "ax[0].set_title('Sum of columns')\n",
    "ax[1].set_title('Sum of rows')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to do this at different angles and to concatonate the resulting 1-D signals. Later in the class, this will be formalised using the **Radon transform**. For now, let's just continue without this machinery."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lets vectorise the image into a vector x\n",
    "x = image.reshape(N*N)\n",
    "print(x.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task. Create a matrix $A$ which, when multiplied by the vectorised image, produces the sum of the columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let visualise a few rows from A (change the value of row and check things make sense)\n",
    "print('The dimensions of A are',A.shape)\n",
    "row = 10\n",
    "plt.imshow(A[row,:].reshape(N,N))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# And we can recalculate the sum of the columns using A (we should get the same as we did before)\n",
    "plt.plot(A@x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task. Add rows to the bottom of $A$ to sum the rows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can now visualise any of the rows of the larger A\n",
    "print('The dimensions of A are',A.shape)\n",
    "row = 70\n",
    "plt.imshow(A[row,:].reshape(N,N))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It should now be clear that we want to continue to add rows to $A$ but with lines taken at different angles. The following function calculates a matrix $A$ with a certain number of projection directions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calcTomographyForwardOperator(numAngles, N):\n",
    "    theta = np.linspace(0, 180, numAngles, endpoint = False)\n",
    "    A = []\n",
    "    E = np.zeros((N, N))\n",
    "    for i_y in tqdm(range(N)):\n",
    "        for i_x in range(N):\n",
    "            E[i_y, i_x] = 1\n",
    "            R_E = radon(E, theta=theta, circle=False)\n",
    "            E[i_y, i_x] = 0\n",
    "            A.append(R_E.flatten())\n",
    "            \n",
    "    return np.array(A).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the forward operator\n",
    "A = calcTomographyForwardOperator(20, N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualise a row\n",
    "print(A.shape)\n",
    "row = 505\n",
    "plt.imshow(A[row,:].reshape(N,N))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that $A$ is a fat matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets calculate our measurements\n",
    "y=A@x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task. Before we use a right inverse, estimate $x$ using $A^\\top$. Plot the resulting image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "x_hat = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's not great but we see something. This reconstruction is known as the backprojection solution. You'll learn about this more when you study tomography.\n",
    "\n",
    "### Task. Estimate $x$ using a right inverse (or an approximation of one): you may need to remove small singular values. Plot the resulting image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is much better than the estimate from back projection. Let's quickly check that it is consistent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mean squared error (MSE) between y and A@x_hat\n",
    "print(np.mean((y - A@x_hat)**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task. Repeat as above but with 50 angles. Verify that you get a tall matrix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "A = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets calculate the measurements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y=A@x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task. Calculate the backprojection estimate. Plot the resulting image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "x_hat = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task. Use a left-inverse (or an approximation of one) to estimate $x$. Plot the resulting image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's it for now. However, you can imagine how using priors for the tomography example could also help the reconstruction. In particular, with well chosen priors, we should be able to use less angles and still achieve a good reconstruction. In practice, this means less x-rays passing through the object (people in the medical scenario) and quicker acquisitions (cheaper).\n",
    "\n",
    "Later in class, we will come back to the tomography example. By then, we will know how to switch between general Hilbert spaces and so we will be able to use more elaborate tools."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.3"
  },
  "widgets": {
   "state": {
    "74390696920b446491ccdaed2c6bc9a6": {
     "views": [
      {
       "cell_index": 4
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}