{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "*This notebook contains material from [cbe67701-uncertainty-quantification](https://ndcbe.github.io/cbe67701-uncertainty-quantification);\n", "content is available [on Github](https://github.com/ndcbe/cbe67701-uncertainty-quantification.git).*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "< [6.3 Sensitivity Analysis with Adjoint Operators](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.03-Sensitivity-Analysis-with-Adjoint-Operators.html) | [Contents](toc.html) | [7.0 Sampling-Based Uncertainty Quantification: Monte Carlo and Beyond](https://ndcbe.github.io/cbe67701-uncertainty-quantification/07.00-Sampling-Based-Uncertainty-Quantification.html)

\"Open

\"Download\"" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "NVXfmm073Zm6", "nbpages": { "level": 1, "link": "[6.4 Adjoint Sensitivity Notes on Numerical Computation](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.04-Adjoint-Sensitivity-Notes-on-Numerical-Computation.html#6.4-Adjoint-Sensitivity-Notes-on-Numerical-Computation)", "section": "6.4 Adjoint Sensitivity Notes on Numerical Computation" } }, "source": [ "# 6.4 Adjoint Sensitivity Notes on Numerical Computation" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "pRj0ROY4N4HT", "nbpages": { "level": 2, "link": "[6.4.1 Notes on Numerical Solution](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.04-Adjoint-Sensitivity-Notes-on-Numerical-Computation.html#6.4.1-Notes-on-Numerical-Solution)", "section": "6.4.1 Notes on Numerical Solution" } }, "source": [ "## 6.4.1 Notes on Numerical Solution\n", "\n", "For applications that can be reduced to the solution of linear equations, the essential computations for the adjoint sensitivity method is to compute solutions to two systems of linear equations\n", "\n", "\\begin{align*}\n", "Ax & = b \\\\\n", "A^Ty & = c\n", "\\end{align*}\n", "\n", "The most obvious why to proceed is to solve these two systems of equation independently of each other. The computational burden to compute that adjoint sensitivity is about the same as to compute the forward solution.\n", "\n", "Alternatively, the solution to these equations can be found by solving the symmetric system\n", "\n", "$$\\underbrace{\\begin{bmatrix} 0 & A \\\\ A^T & 0 \\end{bmatrix}}_P\n", "\\underbrace{\\begin{bmatrix} y \\\\ x \\end{bmatrix}}_z = \\underbrace{\\begin{bmatrix} b \\\\ c \\end{bmatrix}}_f$$\n", "\n", "Iterative algorithms have been developed that simulateously solve both systems of equations with potentially less work then solving them individually. Quoting Lu and Darmofal (2003)\n", "\n", "> * The primal and dual problems are solved simultaneously with essentially the\n", "same computational work as solving only one of the problems with the original\n", "QMR algorithm.\n", "\n", "The purpose of this notebook is to compare the solution time of various iterative algorithms for the solution of sparse systems of linear equations. This naive comparison looks at the relative speedup of solving the combined block matrix equations to the independent solution of the forward and adjoint equations.\n", "\n", "**Note that these computations do not exploit the special structure of the joint forward/adjoint problem as suggested in the following references.**\n" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "gQtMy9cp3l9Z", "nbpages": { "level": 2, "link": "[6.4.2 References](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.04-Adjoint-Sensitivity-Notes-on-Numerical-Computation.html#6.4.2-References)", "section": "6.4.2 References" } }, "source": [ "\n", "## 6.4.2 References\n", "\n", "* Golub, Gene H., Martin Stoll, and Andy Wathen. \"Approximation of the scattering amplitude and linear systems.\" Electron. Trans. Numer. Anal 31.2008 (2008): 178-203. [http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.439.5358&rep=rep1&type=pdf](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.439.5358&rep=rep1&type=pdf)\n", "\n", "* Lu, James, and David L. Darmofal. \"A quasi-minimal residual method for simultaneous primal-dual solutions and superconvergent functional estimates.\" SIAM Journal on Scientific Computing 24.5 (2003): 1693-1709. [http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.144.9176&rep=rep1&type=pdf](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.144.9176&rep=rep1&type=pdf)\n", "\n", "* Pierce, Niles A., and Michael B. Giles. \"Adjoint recovery of superconvergent functionals from PDE approximations.\" SIAM review 42.2 (2000): 247-264.ll [https://epubs.siam.org/doi/pdf/10.1137/S0036144598349423?casa_token=gU8GJSWSZnYAAAAA:jOFJizYd-tUiPVK3jLkOZNRlasVAQ_tasIVaNq8RGOgSJ2MW85CvTnxzW90KnLDin7HW21gnRCQL](https://epubs.siam.org/doi/pdf/10.1137/S0036144598349423?casa_token=gU8GJSWSZnYAAAAA:jOFJizYd-tUiPVK3jLkOZNRlasVAQ_tasIVaNq8RGOgSJ2MW85CvTnxzW90KnLDin7HW21gnRCQL)\n", "\n", "* Stoll, Martin. Solving linear systems using the adjoint. Diss. Oxford University, 2008. See the motivating examples in Chapter 1, and the extensive discussion in Chapter 5. [http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.709.8774&rep=rep1&type=pdf](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.709.8774&rep=rep1&type=pdf)\n", "\n", "* Stoll, Martin, and Andrew J. Wathen. \"All-at-once solution of time-dependent PDE-constrained optimization problems.\" (2010).\n" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "52Da5mk43p73", "nbpages": { "level": 2, "link": "[6.4.3 Comparison of interative solution algorithms from scipy.sparse.linalg](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.04-Adjoint-Sensitivity-Notes-on-Numerical-Computation.html#6.4.3-Comparison-of-interative-solution-algorithms-from-scipy.sparse.linalg)", "section": "6.4.3 Comparison of interative solution algorithms from scipy.sparse.linalg" } }, "source": [ "## 6.4.3 Comparison of interative solution algorithms from scipy.sparse.linalg" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "esAl4Fdy4A1Q", "nbpages": { "level": 3, "link": "[6.4.3.1 Simple test problems](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.04-Adjoint-Sensitivity-Notes-on-Numerical-Computation.html#6.4.3.1-Simple-test-problems)", "section": "6.4.3.1 Simple test problems" } }, "source": [ "### 6.4.3.1 Simple test problems" ] }, { "cell_type": "code", "execution_count": 46, "metadata": { "colab": {}, "colab_type": "code", "executionInfo": { "elapsed": 745, "status": "ok", "timestamp": 1593887491995, "user": { "displayName": "Jeffrey Kantor", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gg_n8V7bVINy02QRuRgOoMo11Ri7NKU3OUKdC1bkQ=s64", "userId": "09038942003589296665" }, "user_tz": 300 }, "id": "TJzRyAnfN48S", "nbpages": { "level": 3, "link": "[6.4.3.1 Simple test problems](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.04-Adjoint-Sensitivity-Notes-on-Numerical-Computation.html#6.4.3.1-Simple-test-problems)", "section": "6.4.3.1 Simple test problems" } }, "outputs": [], "source": [ "import numpy as np\n", "import scipy.sparse\n", "import scipy.sparse.linalg\n", "import matplotlib.pyplot as plt\n", "\n", "spmatrix = sparse.csr_matrix\n", "\n", "# toy problem\n", "def toy_problem():\n", " A = np.array([[1, 1, -1, 0], [2, 1, 0, -1], [0, 0, 1, 0 ], [0, 0, 0, 1]])\n", " b = np.array([0, 0, 80, 100]).reshape((4,1))\n", " c = np.array([170, 120, -50, -40]).reshape((4,1))\n", " return spmatrix(A), b, c\n", "\n", "# random problem\n", "def test_problem(N=25, density=0.2):\n", " A = scipy.sparse.random(N, N, density=density)\n", " b = np.random.rand(N, 1)\n", " c = np.random.rand(N, 1)\n", " return spmatrix(A), b, c" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "5vdoDORv36Tu", "nbpages": { "level": 3, "link": "[6.4.3.2 Create combined sparse block matrix](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.04-Adjoint-Sensitivity-Notes-on-Numerical-Computation.html#6.4.3.2-Create-combined-sparse-block-matrix)", "section": "6.4.3.2 Create combined sparse block matrix" } }, "source": [ "### 6.4.3.2 Create combined sparse block matrix" ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "colab": {}, "colab_type": "code", "executionInfo": { "elapsed": 918, "status": "ok", "timestamp": 1593887492178, "user": { "displayName": "Jeffrey Kantor", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gg_n8V7bVINy02QRuRgOoMo11Ri7NKU3OUKdC1bkQ=s64", "userId": "09038942003589296665" }, "user_tz": 300 }, "id": "rZ5Fp5xWiEPc", "nbpages": { "level": 3, "link": "[6.4.3.2 Create combined sparse block matrix](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.04-Adjoint-Sensitivity-Notes-on-Numerical-Computation.html#6.4.3.2-Create-combined-sparse-block-matrix)", "section": "6.4.3.2 Create combined sparse block matrix" } }, "outputs": [], "source": [ "def combined(A, b, c):\n", " Z = spmatrix(np.zeros(A.shape))\n", " P = scipy.sparse.vstack(\n", " (scipy.sparse.hstack((Z, spmatrix(A))),\n", " scipy.sparse.hstack((spmatrix(A.T), Z))\n", " )\n", " )\n", " P = spmatrix(P)\n", " xs, ys = A.shape\n", " #plt.figure(figsize=(xs/40, ys/40))\n", " #plt.spy(P, ms=1)\n", " f = scipy.sparse.vstack((spmatrix(b), spmatrix(c))).toarray()\n", " return P, f" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "KwoZyBJk4EMa", "nbpages": { "level": 3, "link": "[6.4.3.3 Solver comparisons](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.04-Adjoint-Sensitivity-Notes-on-Numerical-Computation.html#6.4.3.3-Solver-comparisons)", "section": "6.4.3.3 Solver comparisons" } }, "source": [ "### 6.4.3.3 Solver comparisons" ] }, { "cell_type": "code", "execution_count": 73, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 816 }, "colab_type": "code", "executionInfo": { "elapsed": 29629, "status": "ok", "timestamp": 1593887520916, "user": { "displayName": "Jeffrey Kantor", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gg_n8V7bVINy02QRuRgOoMo11Ri7NKU3OUKdC1bkQ=s64", "userId": "09038942003589296665" }, "user_tz": 300 }, "id": "pevh8emwmKie", "nbpages": { "level": 3, "link": "[6.4.3.3 Solver comparisons](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.04-Adjoint-Sensitivity-Notes-on-Numerical-Computation.html#6.4.3.3-Solver-comparisons)", "section": "6.4.3.3 Solver comparisons" }, "outputId": "5988d74e-075d-426e-e0af-7d7281cf522f" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Matrix size = (500, 500). Average of 10 runs.\n", "\n", "Solver: spsolve \n", "\tforward solution: 37.20 ms\n", "\tadjoint solution: 31.19 ms\n", "\tsum of forward and adjoint: 68.39 ms\n", "\tcombined solution: 85.01 ms\n", "\tspeedup: 0.80x\n", "\n", "Solver: cg \n", "\tforward solution: 643.98 ms\n", "\tadjoint solution: 546.12 ms\n", "\tsum of forward and adjoint: 1190.09 ms\n", "\tcombined solution: 288.91 ms\n", "\tspeedup: 4.12x\n", "\n", "Solver: minres \n", "\tforward solution: 94.84 ms\n", "\tadjoint solution: 94.48 ms\n", "\tsum of forward and adjoint: 189.32 ms\n", "\tcombined solution: 422.99 ms\n", "\tspeedup: 0.45x\n", "\n", "Solver: qmr \n", "\tforward solution: 974.82 ms\n", "\tadjoint solution: 967.60 ms\n", "\tsum of forward and adjoint: 1942.43 ms\n", "\tcombined solution: 577.50 ms\n", "\tspeedup: 3.36x\n" ] } ], "source": [ "import time\n", "import warnings\n", "import gc\n", "warnings.filterwarnings(\"ignore\", category=DeprecationWarning) \n", "\n", "solvers = [\n", " scipy.sparse.linalg.spsolve,\n", " scipy.sparse.linalg.cg,\n", " #scipy.sparse.linalg.gmres,\n", " #scipy.sparse.linalg.lgmres,\n", " scipy.sparse.linalg.minres,\n", " scipy.sparse.linalg.qmr,\n", "]\n", "\n", "def solve_time(A, b, solver, repeats=10):\n", " gc.disable()\n", " tic = time.time()\n", " for k in range(repeats):\n", " solver(A, b)\n", " toc = time.time()\n", " gc.enable()\n", " return (toc-tic)/repeats\n", "\n", "def compare_solvers(A, b, c, solvers=solvers, repeats=10):\n", " print(f\"Matrix size = {A.shape}. Average of {repeats} runs.\")\n", " P, f = combined(A, b, c)\n", " for solver in solvers:\n", " print(f\"\\nSolver: {solver.__name__:15s}\")\n", " tf = solve_time(A, b, solver, repeats)\n", " print(f\"\\tforward solution: {1000*tf:8.2f} ms\")\n", " ta = solve_time(A.T, c, solver, repeats)\n", " print(f\"\\tadjoint solution: {1000*ta:8.2f} ms\")\n", " print(f\"\\tsum of forward and adjoint: {1000*(tf + ta):8.2f} ms\")\n", " tc = solve_time(P, f, solver, repeats)\n", " print(f\"\\tcombined solution: {1000*tc:8.2f} ms\")\n", " print(f\"\\tspeedup: {(tf + ta)/tc:8.2f}x\")\n", " \n", "compare_solvers(*test_problem(500, density=0.05))" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 3, "link": "[6.4.3.4 Observations](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.04-Adjoint-Sensitivity-Notes-on-Numerical-Computation.html#6.4.3.4-Observations)", "section": "6.4.3.4 Observations" } }, "source": [ "### 6.4.3.4 Observations\n", "\n", "A speed up of more than 2x suggests that solving the combined problem is actually faster than solving the individual forward or adjoint problems. There are a number of possible explanatins:\n", "\n", "* The symmetry of the combined problem is better conditioned?\n", "* Enhanced convergence?\n", "\n", "A more careful analysis is required that would include factors such as initial conditioning, starting guess, and solution error tolerances." ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[6.4.4 PDE constrained test problem](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.04-Adjoint-Sensitivity-Notes-on-Numerical-Computation.html#6.4.4-PDE-constrained-test-problem)", "section": "6.4.4 PDE constrained test problem" } }, "source": [ "## 6.4.4 PDE constrained test problem\n", "\n", "(Under construction)" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[6.4.4 PDE constrained test problem](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.04-Adjoint-Sensitivity-Notes-on-Numerical-Computation.html#6.4.4-PDE-constrained-test-problem)", "section": "6.4.4 PDE constrained test problem" } }, "source": [ "\\begin{align*}\n", "\\text{PDE:}\\qquad & 0 = \\kappa \\frac{\\partial^2 u}{\\partial z^2} + f & 0 < z < L \\\\\n", "\\text{BC1:}\\qquad & u(0) = u_0 \\\\\n", "\\text{BC2:}\\qquad & u(L) = u_L \\\\\n", "\\end{align*}\n", "\n", "with " ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "colab": {}, "colab_type": "code", "executionInfo": { "elapsed": 29625, "status": "ok", "timestamp": 1593887520918, "user": { "displayName": "Jeffrey Kantor", "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gg_n8V7bVINy02QRuRgOoMo11Ri7NKU3OUKdC1bkQ=s64", "userId": "09038942003589296665" }, "user_tz": 300 }, "id": "IgmjVua-wixO", "nbpages": { "level": 2, "link": "[6.4.4 PDE constrained test problem](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.04-Adjoint-Sensitivity-Notes-on-Numerical-Computation.html#6.4.4-PDE-constrained-test-problem)", "section": "6.4.4 PDE constrained test problem" } }, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "import numpy as np\n", "import scipy.sparse\n", "import scipy.sparse.linalg\n", "import matplotlib.pyplot as plt\n", "\n", "def heat_problem():\n", "\n", " kappa = 10.0\n", " L = 10.0\n", " Nz = 100\n", "\n", " z = np.linspace(0, L, Nz+1)\n", " dz = z[1] - z[0]\n", "\n", " def f(z):\n", " return 1.0\n", "\n", " # diagonals\n", " main = np.zeros(Nz+1)\n", " upper = np.zeros(Nz)\n", " lower = np.zeros(Nz)\n", "\n", " main[1:Nz] = 2*kappa/dz/dz\n", " upper[1:] = -kappa/dz/dz\n", " lower[:-1] = -kappa/dz/dz\n", " main[0] = 1.0\n", " main[-1] = 1.0\n", "\n", " b = np.array([f(z) for z in z])\n", " b[0] = 1.0\n", " b[-1] = 0.0\n", "\n", " c = np.zeros(Nz+1)\n", " c[:] = 1.0\n", "\n", " A = scipy.sparse.diags(\n", " diagonals = [main, lower, upper],\n", " offsets=[0, -1, 1], shape=(Nz+1, Nz+1), format=\"csr\")\n", "\n", " u = scipy.sparse.linalg.spsolve(A, b)\n", " plt.plot(z, u)\n", " \n", " return A, b, c" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "nbpages": { "level": 2, "link": "[6.4.4 PDE constrained test problem](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.04-Adjoint-Sensitivity-Notes-on-Numerical-Computation.html#6.4.4-PDE-constrained-test-problem)", "section": "6.4.4 PDE constrained test problem" } }, "outputs": [ { "data": { "text/plain": [ "array([[0., 0., 0., ..., 0., 0., 0.],\n", " [0., 0., 0., ..., 0., 0., 0.],\n", " [0., 0., 0., ..., 0., 0., 0.],\n", " ...,\n", " [0., 0., 0., ..., 0., 0., 0.],\n", " [0., 0., 0., ..., 0., 0., 0.],\n", " [0., 0., 0., ..., 0., 0., 0.]])" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "nbpages": { "level": 2, "link": "[6.4.4 PDE constrained test problem](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.04-Adjoint-Sensitivity-Notes-on-Numerical-Computation.html#6.4.4-PDE-constrained-test-problem)", "section": "6.4.4 PDE constrained test problem" } }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "< [6.3 Sensitivity Analysis with Adjoint Operators](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.03-Sensitivity-Analysis-with-Adjoint-Operators.html) | [Contents](toc.html) | [7.0 Sampling-Based Uncertainty Quantification: Monte Carlo and Beyond](https://ndcbe.github.io/cbe67701-uncertainty-quantification/07.00-Sampling-Based-Uncertainty-Quantification.html)

\"Open

\"Download\"" ] } ], "metadata": { "colab": { "authorship_tag": "ABX9TyNZ/FAZkGSCLLHgULMr1Up2", "collapsed_sections": [], "name": "06.04-Adjoint-Sensitivity-Notes-on-Numerical-Computation.ipynb", "provenance": [] }, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.4" } }, "nbformat": 4, "nbformat_minor": 4 }