Labs/Lab1_Monday/inverseProblemPrototype.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ill-posedness: An illustrative example\n",
    "\n",
    "## 1D Inverse Heat Equation\n",
    "\n",
    "An illustrative example of ill-posedness is the inversion for the initial condition for a one-dimensional heat equation.\n",
    "\n",
    "Specifically, we consider a rod of length $L$, and let $u(x,t)$ denote the temperature of the rod at point $x$ and time $t$. \n",
    "We are interested in reconstructing the initial temperature profile $m(x) = u(x, 0)$ given some noisy measurements $d$ of the temperature profile at a later time $T$.\n",
    "\n",
    "### Forward problem\n",
    "\n",
    "Given\n",
    "- the initial temperature profile $u(x,0) = m(x)$,\n",
    "- the thermal diffusivity $k$,\n",
    "- a prescribed temperature $u(0,t) = u(L,t) = 0$ at the ends of the rod,\n",
    "\n",
    "solve the heat equation\n",
    "\n",
    "$$\n",
    "\\left\\{\n",
    "\\begin{array}{ll}\n",
    "\\frac{\\partial u}{\\partial t} - k \\frac{\\partial^2}{\\partial x^2} u = 0 & \\forall x\\,\\in\\,(0, L)\\; \\forall t \\in (0,T)\\\\\n",
    "u(x, 0) = m(x) & \\forall x \\in [0,L] \\\\\n",
    "u(0,t) = u(L,t) = 0 & \\forall t \\in (0, T],\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "and observe the temperature at the final time $T$:\n",
    "\n",
    "$$ \\mathcal{F}(m) = u(x, T). $$\n",
    "\n",
    "#### Analytical solution to the forward problem\n",
    "Verify that if\n",
    "\n",
    "$$ m(x) = \\sin\\left(i\\, \\frac{\\pi}{L} x \\right), \\quad i = 1,2,3, \\ldots ,$$\n",
    "\n",
    "then\n",
    "\n",
    "$$ u(x,t) = e^{ -k\\left(i\\, \\frac{\\pi}{L} \\right)^2 t} \\sin\\left(i\\,\\frac{\\pi}{L} x \\right) $$\n",
    "\n",
    "is the unique solution to the heat equation.\n",
    "\n",
    "### Inverse problem\n",
    "\n",
    "Given the forward model $\\mathcal{F}$ and a noisy measurement $d$ of the temperature profile at time $T$, find the initial temperature profile $m$ such that\n",
    "\n",
    "$$ \\mathcal{F}(m) = d. $$\n",
    "\n",
    "### Ill-posedness of the inverse problem\n",
    "\n",
    "Consider a perturbation\n",
    "\n",
    "$$ \\delta m(x) = \\varepsilon \\sin\\left(i \\, \\frac{\\pi}{L} x \\right), $$\n",
    "\n",
    "where $\\varepsilon > 0$ and $i = 1, 2, 3, \\ldots$.\n",
    "\n",
    "Then, by linearity of the forward model $\\mathcal{F}$, the corresponding perturbation $\\delta d(x) = \\mathcal{F}(m + \\delta m) - \\mathcal{F}(m)$ is\n",
    "\n",
    "$$ \\delta d(x) = \\varepsilon\\, e^{ -k\\left(i \\, \\frac{\\pi}{L}\\right)^2 T} \\sin\\left(i \\, \\frac{\\pi}{L} x \\right),$$\n",
    "\n",
    "which converges to zero as $i \\rightarrow +\\infty$.\n",
    "\n",
    "Hence the ratio between $\\delta m$ and $\\delta d$ can become arbitrary large, which shows that the stability requirement for well-posedness can not be satisfied."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Discretization\n",
    "\n",
    "To discretize the problem, we use finite differences in space and Implicit Euler in time.\n",
    "\n",
    "#### Semidiscretization in space\n",
    "We divide the $[0, L]$ interval in $n_x$ subintervals of the same lenght $h = \\frac{L}{n_x}$, and we denote with $u_j(t) := u( jh, t)$ the value of the temperature at point $x_j = jh$ and time $t$.\n",
    "\n",
    "We then use a centered finite difference approximation of the second derivative in space and write\n",
    "\n",
    "$$ \\frac{\\partial u_j(t)}{\\partial t} - k \\frac{u_{j-1}(t) - 2u_j(t) + u_{j+1}(t)}{h^2} \\quad \\text{for } j=1,2,\\ldots,n_x-1,$$\n",
    "\n",
    "with the boundary condition $u_0(t) = u_{n_x}(t) = 0$.\n",
    "\n",
    "We let $n = n_x-1$ be the number of discretization points in the interior of the interval $[0, L]$, and let \n",
    "\n",
    "$$ \\mathbf{u}(t) = \\begin{bmatrix}u_1(t)\\\\u_2(t)\\\\ \\ldots\\\\ u_{n_x-1}(t) \\end{bmatrix} \\in \\mathbb{R}^n $$\n",
    "\n",
    "be the vector collecting the values of the temperature $u(x,t)$ at the points $x_j = j\\,h$ with $j=1,\\ldots,n_x-1$.\n",
    "\n",
    "We then write the system of ordinary differential equations (ODEs):\n",
    "$$ \\frac{\\partial}{\\partial t} \\mathbf{u}(t) + K \\mathbf{u}(t) = 0,$$\n",
    "where $K \\in \\mathbb{R}^{n \\times n}$ is the tridiagonal matrix given by\n",
    "\n",
    "$$ K = \\frac{k}{h^2}\\begin{bmatrix}  2 & -1 &       &        &        &    \\\\\n",
    "                                    -1 &  2 & -1    &        &        &    \\\\\n",
    "                                       & -1 &  2    & -1     &        &    \\\\\n",
    "                                       &    &\\ldots & \\ldots & \\ldots &    \\\\\n",
    "                                       &    &       & -1     &     2  & -1 \\\\ \n",
    "                                       &    &       &        &     -1 & 2  \\\\\n",
    "                     \\end{bmatrix}.$$\n",
    "                     \n",
    "#### Time discretization\n",
    "We subdivide the time interval $(0, T]$ in $n_t$ time step of size $\\Delta t = \\frac{T}{n_t}$.\n",
    "By letting $\\mathbf{u}^{(i)} = \\mathbf{u}(i\\,\\Delta t)$ denote the discretized temperature profile at time $t_i = i\\,\\Delta t$, the Implicit Euler scheme reads\n",
    "\n",
    "$$ \\frac{\\mathbf{u}^{(i+1)} - \\mathbf{u}^{(i)}}{\\Delta t} + K\\mathbf{u}^{(i+1)} = 0, \\quad \\text{for } i=0,1,\\ldots, n_t-1.$$\n",
    "\n",
    "After simple algebraic manipulations and exploiting the initial condition $u(x,0) = m(x)$, we then obtain\n",
    "\n",
    "$$\n",
    "\\left\\{\n",
    "\\begin{array}\n",
    "\\mathbf{u}^{(0)} = \\mathbf{m} \\\\\n",
    "\\mathbf{u}^{(i+1)} = \\left( I + \\Delta t\\, K\\right)^{-1} \\mathbf{u}^{(i)},\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "or equivalently\n",
    "\n",
    "$$ \\mathbf{u}^{(i)} = \\left( I + \\Delta t\\, K\\right)^{-i} \\mathbf{m}.$$\n",
    "\n",
    "In the code below, the function `assembleMatrix` generates the finite difference matrix $\\left( I + \\Delta t\\, K \\right)$ and the function `solveFwd` evaluates the forward model\n",
    "\n",
    "$$ F\\, \\mathbf{m} = \\left( I + \\Delta t\\, K\\right)^{-n_t}\\, \\mathbf{m}. $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function, absolute_import, division\n",
    "\n",
    "import numpy as np\n",
    "import scipy.sparse as sp\n",
    "import scipy.sparse.linalg as la\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "def plot(f, style, **kwargs):\n",
    "    x = np.linspace(0., L, nx+1)\n",
    "    f_plot = np.zeros_like(x)\n",
    "    f_plot[1:-1] = f\n",
    "    plt.plot(x,f_plot, style, **kwargs)\n",
    "    \n",
    "def assembleMatrix(k, h, dt, n):\n",
    "    diagonals = np.zeros((3, n))   # 3 diagonals\n",
    "    diagonals[0,:] = -1.0/h**2\n",
    "    diagonals[1,:] =  2.0/h**2\n",
    "    diagonals[2,:] = -1.0/h**2\n",
    "    K = k*sp.spdiags(diagonals, [-1,0,1], n,n)\n",
    "    M = sp.spdiags(np.ones(n), 0, n,n)\n",
    "    \n",
    "    return M + dt*K\n",
    "    \n",
    "\n",
    "def solveFwd(m, k, h, dt, n, nt):\n",
    "    A = assembleMatrix(k, h, dt, n)\n",
    "    u_old = m.copy()\n",
    "    for i in np.arange(nt):\n",
    "        u = la.spsolve(A, u_old)\n",
    "        u_old[:] = u\n",
    "        \n",
    "    return u        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A naive solution to the inverse problem\n",
    "\n",
    "If $\\mathcal{F}$ is invertible a naive solution to the inverse problem $\\mathcal{F} m = d$ is simply to set\n",
    "\n",
    "$$ m = \\mathcal{F}^{-1} d. $$\n",
    "\n",
    "The function `naiveSolveInv` computes the solution of the discretized inverse problem $\\mathbf{m} = F^{-1} \\mathbf{d}$ as\n",
    "\n",
    "$$ \\mathbf{m} = \\left( I + \\Delta t\\,K\\right)^{n_t} \\mathbf{d}. $$\n",
    "\n",
    "The code below shows that:\n",
    "- for a very coarse mesh (`nx = 20`) and no measurement noise (`noise_std_dev = 0.0`) the naive solution is quite good\n",
    "- for a finer mesh (`nx = 100`) and/or even small measurement noise (`noise_std_dev = 1e-4`) the naive solution is garbage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def naiveSolveInv(d, k, h, dt, n, nt):\n",
    "    A = assembleMatrix(k, h, dt, n)\n",
    "    \n",
    "    p_i = d.copy()\n",
    "    for i in np.arange(nt):\n",
    "        p = A*p_i\n",
    "        p_i[:] = p\n",
    "        \n",
    "    return p\n",
    "\n",
    "T = 1.0\n",
    "L = 1.0\n",
    "k = 0.005"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Very coarse mesh and no measurement noise\")\n",
    "nx = 20\n",
    "nt = 100\n",
    "\n",
    "noise_std_dev = 0.\n",
    "\n",
    "h = L/float(nx)\n",
    "dt = T/float(nt)\n",
    "\n",
    "x = np.linspace(0.+h, L-h, nx-1) #place nx-1 equispace point in the interior of [0,L] interval\n",
    "m_true = 0.5 - np.abs(x-0.5)\n",
    "u_true = solveFwd(m_true, k, h, dt, nx-1, nt)\n",
    "\n",
    "d = u_true + noise_std_dev*np.random.randn(u_true.shape[0])\n",
    "\n",
    "m = naiveSolveInv(d, k, h, dt, nx-1, nt)\n",
    "\n",
    "plt.figure(figsize=(8,4))\n",
    "plt.subplot(1,2,1)\n",
    "plot(m_true, \"-r\", label = 'm_true')\n",
    "plot(m, \"-b\", label = 'm')\n",
    "plt.legend()\n",
    "plt.subplot(1,2,2)\n",
    "plot(u_true, \"-b\", label = 'u(T)')\n",
    "plot(d, \"og\", label = 'd')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Fine mesh and small measurement noise\")\n",
    "nx = 100\n",
    "nt = 100\n",
    "\n",
    "noise_std_dev = 1.e-4\n",
    "\n",
    "h = L/float(nx)\n",
    "dt = T/float(nt)\n",
    "\n",
    "x = np.linspace(0.+h, L-h, nx-1) #place nx-1 equispace point in the interior of [0,L] interval\n",
    "m_true = 0.5 - np.abs(x-0.5)\n",
    "u_true = solveFwd(m_true, k, h, dt, nx-1, nt)\n",
    "\n",
    "d = u_true + noise_std_dev*np.random.randn(u_true.shape[0])\n",
    "\n",
    "m = naiveSolveInv(d, k, h, dt, nx-1, nt)\n",
    "\n",
    "plt.figure(figsize=(8,4))\n",
    "plt.subplot(1,2,1)\n",
    "plot(m_true, \"-r\", label = 'm_true')\n",
    "plot(m, \"-b\", label = 'm')\n",
    "plt.legend()\n",
    "plt.subplot(1,2,2)\n",
    "plot(u_true, \"-b\", label = 'u(T)')\n",
    "plot(d, \"og\", label = 'd')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Why does the naive solution fail?\n",
    "\n",
    "#### Spectral property of the parameter to observable map\n",
    "\n",
    "Let $v_i = \\sqrt{\\frac{2}{L}} \\sin\\left( i \\, \\frac{\\pi}{L} x \\right)$ with $i=1,2,3, \\ldots$, then we have that\n",
    "\n",
    "$$ \\mathcal{F} v_i = \\lambda_i v_i, \\quad \\text{where the eigenvalues } \\lambda_i = e^{-kT\\left(\\frac{\\pi}{L} i \\right)^2}. $$\n",
    "\n",
    "**Note**:\n",
    "- Large eigenvalues $\\lambda_i$ corresponds to smooth eigenfunctions $v_i$;\n",
    "- Small eigenvalues $\\lambda_i$ corresponds to oscillatory eigenfuctions $v_i$.\n",
    "\n",
    "The figure below shows that the eigenvalues $\\lambda_i$ of the continuous parameter to obervable map $\\mathcal{F}$ decays extremely (exponentially) fast."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "T = 1.0\n",
    "L = 1.0\n",
    "k = 0.005\n",
    "\n",
    "i = np.arange(1,50)\n",
    "lambdas = np.exp(-k*T*np.power(np.pi/L*i,2))\n",
    "\n",
    "plt.semilogy(i, lambdas, 'ob')\n",
    "plt.xlabel('i')\n",
    "plt.ylabel('lambda_i')\n",
    "plt.title(\"Eigenvalues of continuous F\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a similar way, the figure below show the eigenvalues of the discrete parameter to observable map $F$: their fast decay means that $F$ is extremely ill conditioned.\n",
    "\n",
    "In the code below we assemble the matrix $F$ column-by-column, by computing its actions on the canonical vectors \n",
    "$$\\mathbf{m}_i = \\begin{bmatrix} 0 \\\\ \\ldots \\\\ 0\\\\ 1 \\\\ 0\\\\ \\ldots \\\\0 \\end{bmatrix}, \\quad i = 1,\\ldots,n,$$\n",
    "where the $i$th entry is the only non-zero component of $\\mathbf{m}_i$.\n",
    "\n",
    "> **Disclaimer**: $F$ is a large dense implicitly defined operator and should never be built explicitly for a real problem (since it would require $\\mathcal{O}(n)$ evaluations of the forward problem and $\\mathcal{O}( n^2)$ storage); instead --- as you will learn later this week --- scalable algorithms for the solution of the inverse problem only require the ability to compute the action of $F$ on a few given directions $\\mathbf{m}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def computeEigendecomposition(k, h, dt, n, nt):\n",
    "    ## Compute F as a dense matrix\n",
    "    F = np.zeros((n,n))\n",
    "    m_i = np.zeros(n)\n",
    "    \n",
    "    for i in np.arange(n):\n",
    "        m_i[i] = 1.0\n",
    "        F[:,i] = solveFwd(m_i, k, h, dt, n, nt)\n",
    "        m_i[i] = 0.0\n",
    "    \n",
    "    ## solve the eigenvalue problem\n",
    "    lmbda, U = np.linalg.eigh(F)\n",
    "    ## sort eigenpairs in decreasing order\n",
    "    lmbda[:] = lmbda[::-1]\n",
    "    lmbda[lmbda < 0.] = 0.0\n",
    "    U[:] = U[:,::-1]\n",
    "    \n",
    "    return lmbda, U \n",
    "\n",
    "## Compute eigenvector and eigenvalues of the discretized forward operator\n",
    "lmbda, U = computeEigendecomposition(k, h, dt, nx-1, nt)\n",
    "\n",
    "plt.semilogy(lmbda, 'ob')\n",
    "plt.title(\"Eigenvalues of discrete F\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Informed and uninformed modes\n",
    "\n",
    "The functions $v_i$ ($i=1,2,3, \\ldots$) form an orthonormal basis of $L^2([0,1])$. \n",
    "\n",
    "That is, every function $f \\in L^2([0,1])$ can be written as\n",
    "\n",
    "$$ f = \\sum_{i=1}^\\infty \\alpha_i v_i, \\text{ where } \\alpha_i = \\int_0^1 f v_i dx.$$\n",
    "\n",
    "Consider now the noisy problem\n",
    "\n",
    "$$ d = \\mathcal{F}\\,m_{\\rm true} + \\eta, $$\n",
    "\n",
    "where\n",
    "- $d$ is the data (noisy measurements)\n",
    "- $\\eta$ is the noise: $\\eta(x) = \\sum_{n=1}^\\infty \\eta_n v_n(x)$\n",
    "- $m_{\\rm true}$ is the true value of the parameter that generated the data\n",
    "- $\\mathcal{F}$ is the forward heat equation\n",
    "\n",
    "Then, the naive solution to the inverse problem $\\mathcal{F}m = d$ is\n",
    "\n",
    "$$ m = \\mathcal{F}^{-1}d = \\mathcal{F}^{-1}\\left( \\mathcal{F}\\,m_{\\rm true} + \\eta \\right) = m_{\\rm true} + \\mathcal{F}^{-1} \\eta = m_{\\rm true} + \\mathcal{F}^{-1} \\sum_{i=1}^{\\infty} \\eta_i v_i = m_{\\rm true} +  \\sum_{i=1}^{\\infty} \\frac{\\eta_i}{\\lambda_i} v_i. $$\n",
    "\n",
    "If the coefficients $\\eta_i = \\int_0^1 \\eta(x) \\, v_i(x) \\, dx$ do not decay sufficiently fast with respect to the eigenvalues $\\lambda_i$, then the naive solution is unstable.\n",
    "\n",
    "This implies that oscillatory components can not reliably be reconstructed from noisy data since they correspond to small eigenvalues.\n",
    "\n",
    "### Tikhonov regularization\n",
    "\n",
    "The remedy is to find a parameter $m$ that solves the inverse problem $\\mathcal{F}\\, m = d$ in a *least squares sense*.\n",
    "Specifically, we solve the minimization problem\n",
    "\n",
    "$$ \\min_m \\frac{1}{2} \\int_0^L (\\mathcal{F}\\, m - d )^2 dx + \\frac{\\alpha}{2} \\mathcal{R}(m), $$\n",
    "\n",
    "where the Tikhonov regularization $\\mathcal{R}(m)$ is a quadratic functional of $m$.\n",
    "In what follow, we will penalize the $L^2$-norm of the initial condition and let $\\mathcal{R}(m) = \\int_0^L m^2 dx$. However other choices are possible; for example by penalizing the $L^2$-norm of the gradient of $m$ ($\\mathcal{R}(m) = \\int_0^L m_x^2$), one will favor smoother solutions.\n",
    "\n",
    "The regularization parameter $\\alpha$ needs to be chosen appropriately. If $\\alpha$ is small, the computation of the initial condition $m$ is unstable as in the naive approach. On the other hand, if $\\alpha$ is too large, information is lost in the reconstructed $m$. Various criteria --- such as Morozov's discrepancy principle or L-curve criterion  --- can be used to find the optimal amount of regularization (see the **Hands on** section).\n",
    "\n",
    "#### Morozov's discrepancy principle\n",
    "\n",
    "The discrepancy principle, due to Morozov, chooses the regularization parameter to be the largest value of $\\alpha$ such that the norm of the misfit is bounded by the noise level in the data, i.e.,\n",
    "\n",
    "$$ \\| \\mathcal{F}\\,m_\\alpha - d \\| \\leq \\delta, $$ \n",
    "\n",
    "where $\\delta$ is the noise level. Here, $m_\\alpha$ denotes the parameter found minimizing the Tikhonov regularized minimization problem with parameter $\\alpha$. This choice aims to avoid overfitting of the data, i.e., fitting the noise.\n",
    "\n",
    "#### L-curve criterion\n",
    "\n",
    "Choosing $\\alpha$ parameter using the L-curve criterion requires the solution of inverse problems for a sequence of regularization parameters. Then, for each $\\alpha$, the norm of the data misfit (also called residual) $\\| \\mathcal{F}\\,m_\\alpha - d \\|$ is plotted against the norm of the regularization term $\\|m_\\alpha\\|$ in a log-log plot.\n",
    "This curve usually is found to be L-shaped and thus has an *elbow*, i.e. a point of greatest curvature. The L-curve criterion chooses the regularization parameter corresponding to that point. The idea behind the L-curve criterion is that this choice for the regularization parameter is a good compromise between fitting the data and controlling the stability of the parameters. A smaller $\\alpha$, which correspond to points to the left of the optimal value, only leads to a slightly better data fit while significantly increasing the norm of the parameters. Conversely, a larger $\\alpha$, corresponding to points to the right of the optimal value, slightly decrease the norm of the solution, but they increase the data misfit significantly.\n",
    "\n",
    "\n",
    "### Discretization\n",
    "In the discrete setting, the Tikhonov regularized solution $\\mathbf{m}_{\\alpha}$ solves the penalized least squares problem\n",
    "\n",
    "$$ \\min_{\\mathbf{m}} \\frac{1}{2} \\| F\\, \\mathbf{m} - \\mathbf{d} \\|^2 + \\frac{\\alpha}{2} \\| \\mathbf{m} \\|^2, $$\n",
    "\n",
    "where $\\| \\cdot \\|$ denotes the Euclidean vector norm in $\\mathbb{R}^n$.\n",
    "\n",
    "$\\mathbf{m}_{\\alpha}$ can then be computed by solving the normal equations\n",
    "\n",
    "$$ ( F^t F + \\alpha I) \\mathbf{m}_{\\alpha} = F^t \\mathbf{d}.$$\n",
    "\n",
    "The code below, find the Tikhonov regularized solution $\\mathbf{m}_{\\alpha}$ for $\\alpha = 10^{-3}$.\n",
    "\n",
    "> **Disclaimer**: In the code below, for simplicity, we explicitly construct and factorize the matrix $F^t F + \\alpha I$. This approach is not feasible and should **never** be used to solve real problems (the computational cost is  $\\mathcal{O}(n^3)$). Instead, as we will see tomorrow, one should solve the normal equations using the conjugate gradient algorithm, which only requires the ability to compute the action of $F^t F + \\alpha I$ on a few given directions $\\mathbf{m}$, and is guaranteed to converge in a number of iterations that is independent of $n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def assembleF(k, h, dt, n, nt):\n",
    "    F = np.zeros((n,n))\n",
    "    m_i = np.zeros(n)\n",
    "    \n",
    "    for i in np.arange(n):\n",
    "        m_i[i] = 1.0\n",
    "        F[:,i] = solveFwd(m_i, k, h, dt, n, nt)\n",
    "        m_i[i] = 0.0\n",
    "\n",
    "    return F\n",
    "\n",
    "def solveTikhonov(d, F, alpha):    \n",
    "    H = np.dot( F.transpose(), F) + alpha*np.identity(F.shape[1])\n",
    "    rhs = np.dot( F.transpose(), d)\n",
    "    return np.linalg.solve(H, rhs)\n",
    "\n",
    "## Setup the problem\n",
    "T = 1.0\n",
    "L = 1.0\n",
    "k = 0.005\n",
    "\n",
    "nx = 100\n",
    "nt = 100\n",
    "\n",
    "noise_std_dev = 1e-3\n",
    "\n",
    "h = L/float(nx)\n",
    "dt = T/float(nt)\n",
    "\n",
    "## Compute the data d by solving the forward model\n",
    "x = np.linspace(0.+h, L-h, nx-1)\n",
    "#m_true = np.power(.5,-36)*np.power(x,20)*np.power(1. - x, 16)\n",
    "m_true = 0.5 - np.abs(x-0.5)\n",
    "u_true = solveFwd(m_true, k, h, dt, nx-1, nt)\n",
    "d = u_true + noise_std_dev*np.random.randn(u_true.shape[0])\n",
    "\n",
    "alpha = 1e-3\n",
    "\n",
    "F = assembleF(k, h, dt, nx-1, nt)\n",
    "m_alpha = solveTikhonov(d, F, alpha)\n",
    "\n",
    "plot(m_true, \"-r\", label = 'm_true')\n",
    "plot(m_alpha, \"-g\", label = 'm_tikh')\n",
    "plt.title(\"Solution\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hands on\n",
    "\n",
    "## Problem 1: Spectrum of the parameter to observable map\n",
    "\n",
    "### Question a:\n",
    "Set $T=1$, $L=1$, $n_x = 200$, and $n_t = 100$, and plot the\n",
    "decay of the eigenvalues of both the continuous ($\\mathcal{F}$) and\n",
    "the discrete ($F$) parameter-to-observable maps as a function of the index $i$.\n",
    "\n",
    "Do this for the following values of $k$: 0.0001, 0.001, 0.01, 0.1, 1.0 (all on the same plot)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question b:\n",
    "Set $L=1$, $k = 0.01$ and $T=0.1$; plot the decay of the discrete eigenvalues as a function of $n_x, n_t$ for different resolutions in the space and time discretization. \n",
    "\n",
    "Use $(n_x, n_t) = (20, 20), (40, 40), (80,80), (160, 160)$.\n",
    "\n",
    "What do you observe as you increase the resolution? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 2: Tikhonov Regularization\n",
    "\n",
    "Consider the same inverse heat equation problem above\n",
    "with $L=1$, $T=0.1$, $k=0.01$.\n",
    "\n",
    "Discretize the problem using $n_x = 200$ intervals in space and $n_t = 100$ time steps.\n",
    "\n",
    "As initial condition use the (discrete version of) the true initial temperature profile \n",
    "\n",
    "$$\n",
    "m_{\\text{true} } = \\max(0, 1 - |1 - 4x|) + 100\\, x^{10}(1-x)^2.\n",
    "$$\n",
    "\n",
    "Use the code below to implement the above function in Python:\n",
    "```\n",
    "import numpy as np\n",
    "x = np.arange(1,n_x, dtype=np.float64)*h\n",
    "m_true = np.maximum( np.zeros_like(x), 1. - np.abs(1. - 4.*x)) \\\n",
    "         + 100.*np.power(x,10)*np.power(1.-x,2)\n",
    "```\n",
    "\n",
    "Add normally distributed noise $\\mathbf{n}$ with mean zero and\n",
    "variance $\\sigma^2 = 10^{-4}$. The resulting noisy observation of the\n",
    "final time temperature profile is $\\mathbf{d} = F \\mathbf{m} + \\mathbf{n}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question a\n",
    "Use Tikhonov regularization with $\\alpha = 0.0001,0.001,0.01,0.1,1$ to compute the regularized\n",
    "reconstructions $\\mathbf{m}_\\alpha$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question b\n",
    "Determine the (approximate) optimal value of the regularization parameter $\\alpha$ in the Tikhonov regularization using the L-curve criterion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question c\n",
    "Determine the (approximate) optimal value of the regularization parameter $\\alpha$ in the Tikhonov regularization \n",
    "using Morozov's discrepancy criterion, i.e., find the largest value of $\\alpha$ such that \n",
    "\n",
    "$$\n",
    "\\|F\\,\\mathbf{m}_\\alpha - \\mathbf{d}\\| \\le \\delta\n",
    "$$\n",
    "\n",
    "where $\\delta=\\| \\mathbf{n}\\|$ and $\\mathbf{m}_\\alpha$ is the solution of the Tikhonov-regularized inverse problem with regularization parameter $\\alpha$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question d\n",
    "Plot the  $L_2$ norm error in the reconstruction, $\\|\\mathbf{m}_{\\text{true}}-\\mathbf{m}_\\alpha\\|$, as a function of $\\alpha$, where $\\mathbf{m}_\\alpha$ is the Tikhonov regularized solution. Which value of $\\alpha$ (approximately) minimizes this error? Compare the *optimal* values of $\\alpha$ obtained using the L-curve and Morozov's discrepancy criterion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 3: Tikhonov regularization on the gradient of $m$\n",
    "\n",
    "Consider the Tikhonov regularized solution $m_\\alpha$ where we penalize the norm of the *gradient* of $m$.\n",
    "That is, $m_\\alpha$ solves minimization problem\n",
    "\n",
    "$$ \\min_m \\frac{1}{2} \\int_0^L (\\mathcal{F}\\, m - d )^2 dx + \\frac{\\alpha}{2} \\int_0^L m_x^2 \\, dx. $$\n",
    "\n",
    "### Question a\n",
    "\n",
    "Discretize and solve the inverse problem using the same parameters and true initial condition as in Problem 2.\n",
    "\n",
    "### Question b\n",
    "\n",
    "Determine the *optimal* amount of regularization using either the L-curve criterion or the discrepancy principle.\n",
    "How does this solution differs from the one computed in Problem 2?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Copyright &copy; 2018, The University of Texas at Austin & University of California, Merced. All Rights reserved. See file COPYRIGHT for details.\n",
    "\n",
    "This file is part of the hIPPYlib library. For more information and source code availability see https://hippylib.github.io.\n",
    "\n",
    "hIPPYlib is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License (as published by the Free Software Foundation) version 2.0 dated June 1991."
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}