Session 3 - Tikhonov Regularization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Session 3 - Tikhonov Regularization\n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    " \n",
    "**By the end of this session you will be able to**\n",
    "\n",
    "- Evaluate how measurement noise influences the least-squares inversion scheme.\n",
    "- Apply regularization to smooth the least-squares inversion results.\n",
    "- Apply the L-curve to identify the optimal regularization parameter.\n",
    "\n",
    "</div>\n",
    "\n",
    "\n",
    "<div class=\"alert alert-danger\">\n",
    "\n",
    "**Instructions**\n",
    "    \n",
    "- Go to Menu \"Cell\" -> \"Run All Below\"\n",
    "- The notebook will initialize all cells \n",
    "- Navigate back to this point and work through the notebook sequentially.\n",
    "- When finished, move to Session 4.\n",
    " \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "include(\"Scripts/startup.jl\")\n",
    "play(\"Audio/Tikhonov1.ogg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The DMA fo this example\n",
    "qsa,qsh = 1.66e-5, 8.33e-5                       # Qsample [m3 s-1], Qsheath [m3 s-1]\n",
    "t,p = 295.15, 1e5                                # Temperature [K], Pressure [Pa]\n",
    "r₁,r₂,l = 9.37e-3,1.961e-2,0.44369               # DMA geometry [m]\n",
    "Λ = DMAconfig(t,p,qsa,qsh,r₁,r₂,l,0.0,:-,6,:cylindrical)  \n",
    "\n",
    "bins,z₁,z₂ = 60, vtoz(Λ,10000), vtoz(Λ,10)    # bins, upper, lower mobility limit\n",
    "δ = setupDMA(Λ, z₁, z₂, bins)                 # Setup DMA grid\n",
    "𝐀 = δ.𝐀;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Least-Squares Inversion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "play(\"Audio/Tikhonov2.ogg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "The response function is measured. The number concentration can be found by solving the equation for 𝕟 <br><br>\n",
    "    \n",
    "    \n",
    "<center> $ \\rm{𝕟} = \\bf{A}^{-1}\\rm{𝕣}$ </center><br>\n",
    "This simple inversion approach can be sufficient in some cases, However, even seemingly small noise  is amplified leading to unusable estimates of $\\rm{𝕟}$ (e.g., Kandlikar and Ramachandran, 1999). \n",
    "    \n",
    "</div>\n",
    "\n",
    "\n",
    "Another way to look at the inversion described above is to cast the problem as a general least-squares problem\n",
    "\n",
    "$L_1 = \\left\\lVert \\bf{A}\\rm{𝕟} - \\rm{𝕣}\\right\\rVert_2$ <br>\n",
    "\n",
    "where $L_1$ is the Euclidean norm of the residual. The optimal inverse solution is defined such that $L_1^2$ is minimized.\n",
    "\n",
    "$𝕟_{inv} = \\arg \\min\\{\\left\\rVert\\bf{A}\\rm{𝕟}-\\rm{𝕣}\\right\\rVert_2^2$}\n",
    "\n",
    "$ \\rm{𝕟} = \\bf{A}^{-1}\\rm{𝕣}$ is the analytical solution that minimizes $L_1$. \n",
    "\n",
    "If the matrix $\\bf{A}$ is rank deficient and not invertible the pseudo-inverse will yield the solution with the lowest $L_1$ norm\n",
    "\n",
    "$\\rm{𝕟} = \\bf{({A^TA})^{-1}A^T} \\rm{𝕣}$ \n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "Noisy least-square inverse solutions have a large residual error and thus a large $L_1$ norm.\n",
    "    \n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-danger\">\n",
    "    \n",
    "The solution is often dominated by contributions from data errors and rounding errors. Measurement noise is amplified in the inversion and the simple inverse produces unusable results. Measurement noise increases with increasing binsize and decreasing sample flow rate. The interactive applet demonstrates how these factors influence the inverted size distribution.\n",
    "    \n",
    "</div>\n",
    "\n",
    "**Measurement noise:** Typically a condensation particle counter is used to measure the response function. Noise in the condensation particle counter is mostly from counting statistics. The Poisson counting statistic depends on the number of counts per bin, and thus the sample flow rate and integration time per bin. In the app below the scan time for the entire distribution is fixed to 2 min. Both the number of bins and the sample flow rate influence the measurement noise. \"Noise free\" CPC means that there is no counting error in the CPC.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "play(\"Audio/Tikhonov3.ogg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "1. Predict how changing total number concentration, CPC flow rate, and number of bins will influence the noise in the inverted size distribution.\n",
    "    \n",
    "2. Use the applet above to verify or falsify your prediction.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "include(\"Scripts/noise_app.jl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Figure.** Influence of measurement noise on the least-squares solution.\n",
    "\n",
    "\n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    "   \n",
    "**You met the following learning objective**\n",
    "    \n",
    "- Evaluate how measurement noise influences the least-squares inversion scheme.\n",
    "    \n",
    "</div>\n",
    "\n",
    "# 2. Tikhonov Inversion\n",
    "\n",
    "## 2.1 Analytical Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "play(\"Audio/Tikhonov4.ogg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "This is also known as zeroth order Tikhonov regularization, or Phillips–Twomey regularization. To filter noise a regularization term is added\n",
    "\n",
    "$L_2 = \\left\\lVert\\bf{L}(\\rm{𝕟} - \\rm{𝕟_i})\\right\\rVert_2$\n",
    "\n",
    "where $\\bf{L}$ is a weights matrix and $\\rm{𝕟_i}$ is an initial guess. The $L_2$ norm describes the deviation from a smooth initial guess. The inverse is a solution that balance between the $L_1$ and $L_2$ norms.\n",
    "\n",
    "$𝕟_{inv} = \\arg \\min\\{\\left\\rVert\\bf{A}\\rm{𝕟}-\\rm{𝕣}\\right\\rVert_2^2 + \\lambda^2 \\left\\rVert\\bf{L}(\\rm{𝕟}-\\rm{𝕟_i})\\right\\rVert_2^2$\n",
    "\n",
    "where $\\lambda$ is the regularization parameter. \n",
    " \n",
    " <div class=\"alert alert-warning\">\n",
    "\n",
    "The regularization parameter $\\lambda$ \"interpolates\" between the noisy least-square inverse ($\\lambda = 0$) and the initial guess ($\\lambda >> 1$) by placing weights on either the $L_1$ or the $L_2$ norm. \n",
    "    \n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "play(\"Audio/Tikhonov5.ogg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Taking $\\bf{L} = \\bf{I}$ (weight matrix equals the identity matrix) and $\\rm{𝕟_i} = \\bf{S}^{-1}\\rm{𝕣}$ as initial guess, the regularized inverse is computed via:\n",
    "\n",
    "$𝕟_{inv} = (\\bf{A}^\\rm{T}\\bf{A} + \\lambda^\\rm{2} \\bf{I})^\\rm{-1}(\\bf{A}^\\rm{T} \\rm{𝕣} + \\lambda^2\\bf{S}^{-1} \\rm{𝕣})$ <br>\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<div class=\"alert alert-danger\">\n",
    "    \n",
    "**To prove that** $𝕟_{inv} = (\\bf{A}^\\rm{T}\\bf{A} + \\lambda^\\rm{2} \\bf{I})^\\rm{-1}(\\bf{A}^\\rm{T} \\rm{𝕣} + \\lambda^2\\bf{S}^{-1} \\rm{𝕣})$ **is the correct solution**\n",
    " \n",
    "- Set $\\bf{L} = \\bf{I}$\n",
    "- Evaluate $\\frac{\\partial}{\\partial \\rm{𝕟}}\\left(\\left\\rVert\\bf{A}\\rm{𝕟}-\\rm{𝕣}\\right\\rVert_2^2 + \\lambda^2 \\left\\rVert\\bf{I}(\\rm{𝕟}-\\rm{𝕟_i})\\right\\rVert_2^2\\right) = 0$\n",
    "- Use http://www.matrixcalculus.org/ to evaluate the derivative if you forgot your symbolic matrix calculus\n",
    "- Substitute $\\rm{𝕟_i} = \\bf{S}^{-1}\\rm{𝕣}$ \n",
    "- Solve for 𝕟\n",
    "   \n",
    "</div>\n",
    "\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "Show that for $\\lambda = 0$, the Tikhonov inverse reduces to the regular-least-squares inverse: $\\rm{𝕟} = \\bf{({A^TA})^{-1}A^T} \\rm{𝕣}$.\n",
    "    \n",
    "</div>\n",
    "\n",
    "# 2.2 Regularization Paramater\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "play(\"Audio/Tikhonov6.ogg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "include(\"Scripts/regularization_app1.jl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Figure**. Left: Least-squares solution, Middle: Tikhonov with λ selected by the slider, Right: initial guess.\n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "- Talukdar and Swihart (2003) introduced a clever initial guess: sum the rows of 𝐀 and place the results on the diagonal of 𝐒. The least-squares inverse $\\rm{𝕟_i} = \\bf{S}^{-1}\\rm{𝕣}$ gives excellent results where singly charged particles dominate and slight error where multiply charged particles become important.\n",
    "- The constraint of the initial condition results in a robust inversion.\n",
    "- The algorithm provides a reasonable inversion for a wide range of λ.\n",
    "       \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "play(\"Audio/Tikhonov7.ogg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "include(\"Scripts/regularization_app2.jl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "play(\"Audio/Tikhonov8.ogg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Regularization in code\n",
    "λ, bins = 0.3, length(δ.Dp)\n",
    "\n",
    "# some response distribtion\n",
    "𝕣 = DMALognormalDistribution([[400, 30, 1.2],[500, 110, 1.7]], δ)\n",
    "\n",
    "# Identity, Talukda and Swihart, and Convolution Matrices for DMA grid \n",
    "𝐈, 𝐒, 𝐀 =  δ.𝐈, δ.𝐒, δ.𝐀\n",
    "\n",
    "𝕟ⁱⁿᵛ = (𝐀'𝐀 + λ^2𝐈)^(-1) * (𝐀'𝕣  + λ^2 * 𝐒^(-1)*𝕣)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<div class=\"alert alert-danger\">\n",
    "\n",
    "- For convenience, the matrices 𝐈, 𝐒, 𝐀 are computed when the DMA grid is initialized. \n",
    "- They are stored in the δ structure together with the other grid variables.\n",
    "       \n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    "   \n",
    "**You met the following learning objective**\n",
    "    \n",
    "- Apply regularization to smooth the least-squares inversion results.\n",
    "    \n",
    "</div>\n",
    "\n",
    "## 2.3 The L-Curve Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "play(\"Audio/Tikhonov9.ogg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "- This is the L-curve for the case with strong constraints.\n",
    "       \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "include(\"Scripts/lcurve1.jl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Figure.** Left: Illustration of the L-curve. Right: Size distribution from inversion using the optimum regularization parameter. The raw distribution corresponds to N = 500 cm-3, Q = 1 L min-1 and 120 bins in the and is the same distribution as in the \"regularization_app1\". \n",
    "\n",
    "\n",
    "The optimal $\\lambda_{opt}$ is found using the L-curve method. The L-curve is defined as a plot of $\\log(L_1)$ vs.  $\\log(L_2)$, where $L_1$ and $L_2$ are obtained for a series of discrete $\\lambda$ values. The $\\lambda_{opt}$ is found at the corner of the L-curve, which is mathematically defined as the point where where the curvature of the L-curve is maximum. Here, the corner of the L-curve is found using the iterative algorithm described in Talukdar and Swihart (2003). The curvature is calculated using Eq. (14) in Hansen (2000), which requires the first and second derivatives of $d\\ln(L_i)^2/d\\lambda$. These derivatives of $d\\ln(L_i)^2/d\\lambda$ are estimated numerically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "play(\"Audio/Tikhonov10.ogg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "- This is the L-curve for if the initial guess is zero.\n",
    "- Tikhonov + L-curve still finds an acceptable solution, though the solution is noiser.\n",
    "       \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "include(\"Scripts/lcurve2.jl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Figure.** Left: Illustration of the L-curve. Right: Size distribution from inversion using the optimum regularization parameter. The raw distribution corresponds to N = 500 cm-3, Q = 1 L min-1 and 120 bins in the and is the same distribution as in the \"regularization_app2\". \n",
    "\n",
    "\n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    "   \n",
    "**You met the following learning objective**\n",
    "\n",
    "- Apply the L-curve method to identify the optimal regularization parameter.\n",
    "\n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-danger\">\n",
    "\n",
    "- The L-curve search is slow due to the need to compute derivatives and rerun the matrix inversion many times. For this specific problem, the initial guess is very good. For \"fast\" inversion and using the initial guess, a characteristic regularization parameter ~0.2 to 0.6 can be used to give adequate results. The value to pick will depend on the CPC and DMA.\n",
    "    \n",
    "- A variety of methods are available to identify the optimal $\\lambda$. These may be consulted if the L-curve minimum approach fails due to difficulty of estimating derivatives.\n",
    "   \n",
    "- Higher order Tikhonov solutions use estimates other then $\\bf{L} = \\bf{I}$.\n",
    "    \n",
    "- The following papers give excellent additional information about Tikhonov regularization: Twomey (1963), Kandlikar & Ramachandran (1999), Hansen (2000), Sipkens et al. (2020).\n",
    "</div>"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.4.2",
   "language": "julia",
   "name": "julia-1.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}