{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Matrix Multiplication\n\nThis example shows how to use the :py:class:`pylops.MatrixMult` operator\nto perform *Matrix inversion* of the following linear system.\n\n\\begin{align}\\mathbf{y}= \\mathbf{A} \\mathbf{x}\\end{align}\n\nYou will see that since this operator is a simple overloading to a\n:py:func:`numpy.ndarray` object, the solution of the linear system\ncan be obtained via both direct inversion (i.e., by means explicit\nsolver such as :py:func:`scipy.linalg.solve` or :py:func:`scipy.linalg.lstsq`)\nand iterative solver (i.e., :py:func:`from scipy.sparse.linalg.lsqr`).\n\nNote that in case of rectangular $\\mathbf{A}$, an exact inverse does\nnot exist and a least-square solution is computed instead.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import warnings\n\nimport matplotlib.gridspec as pltgs\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom scipy.sparse import rand\nfrom scipy.sparse.linalg import lsqr\n\nimport pylops\n\nplt.close(\"all\")\nwarnings.filterwarnings(\"ignore\")\n# sphinx_gallery_thumbnail_number = 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's define the sizes of the matrix $\\mathbf{A}$ (``N`` and ``M``) and\nfill the matrix with random numbers\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "N, M = 20, 20\nA = np.random.normal(0, 1, (N, M))\nAop = pylops.MatrixMult(A, dtype=\"float64\")\n\nx = np.ones(M)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now apply the forward operator to create the data vector $\\mathbf{y}$\nand use ``/`` to solve the system by means of an explicit solver.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "y = Aop * x\nxest = Aop / y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's visually plot the system of equations we just solved.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "gs = pltgs.GridSpec(1, 6)\nfig = plt.figure(figsize=(7, 3))\nax = plt.subplot(gs[0, 0])\nax.imshow(y[:, np.newaxis], cmap=\"rainbow\")\nax.set_title(\"y\", size=20, fontweight=\"bold\")\nax.set_xticks([])\nax.set_yticks(np.arange(N - 1) + 0.5)\nax.grid(linewidth=3, color=\"white\")\nax.xaxis.set_ticklabels([])\nax.yaxis.set_ticklabels([])\nax = plt.subplot(gs[0, 1])\nax.text(\n 0.35,\n 0.5,\n \"=\",\n horizontalalignment=\"center\",\n verticalalignment=\"center\",\n size=40,\n fontweight=\"bold\",\n)\nax.axis(\"off\")\nax = plt.subplot(gs[0, 2:5])\nax.imshow(Aop.A, cmap=\"rainbow\")\nax.set_title(\"A\", size=20, fontweight=\"bold\")\nax.set_xticks(np.arange(N - 1) + 0.5)\nax.set_yticks(np.arange(M - 1) + 0.5)\nax.grid(linewidth=3, color=\"white\")\nax.xaxis.set_ticklabels([])\nax.yaxis.set_ticklabels([])\nax = plt.subplot(gs[0, 5])\nax.imshow(x[:, np.newaxis], cmap=\"rainbow\")\nax.set_title(\"x\", size=20, fontweight=\"bold\")\nax.set_xticks([])\nax.set_yticks(np.arange(M - 1) + 0.5)\nax.grid(linewidth=3, color=\"white\")\nax.xaxis.set_ticklabels([])\nax.yaxis.set_ticklabels([])\nplt.tight_layout()\n\ngs = pltgs.GridSpec(1, 6)\nfig = plt.figure(figsize=(7, 3))\nax = plt.subplot(gs[0, 0])\nax.imshow(x[:, np.newaxis], cmap=\"rainbow\")\nax.set_title(\"xest\", size=20, fontweight=\"bold\")\nax.set_xticks([])\nax.set_yticks(np.arange(M - 1) + 0.5)\nax.grid(linewidth=3, color=\"white\")\nax.xaxis.set_ticklabels([])\nax.yaxis.set_ticklabels([])\nax = plt.subplot(gs[0, 1])\nax.text(\n 0.35,\n 0.5,\n \"=\",\n horizontalalignment=\"center\",\n verticalalignment=\"center\",\n size=40,\n fontweight=\"bold\",\n)\nax.axis(\"off\")\nax = plt.subplot(gs[0, 2:5])\nax.imshow(Aop.inv(), cmap=\"rainbow\")\nax.set_title(r\"A$^{-1}$\", size=20, fontweight=\"bold\")\nax.set_xticks(np.arange(N - 1) + 0.5)\nax.set_yticks(np.arange(M - 1) + 0.5)\nax.grid(linewidth=3, color=\"white\")\nax.xaxis.set_ticklabels([])\nax.yaxis.set_ticklabels([])\nax = plt.subplot(gs[0, 5])\nax.imshow(y[:, np.newaxis], cmap=\"rainbow\")\nax.set_title(\"y\", size=20, fontweight=\"bold\")\nax.set_xticks([])\nax.set_yticks(np.arange(N - 1) + 0.5)\nax.grid(linewidth=3, color=\"white\")\nax.xaxis.set_ticklabels([])\nax.yaxis.set_ticklabels([])\nplt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's also plot the matrix eigenvalues\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure(figsize=(8, 3))\nplt.plot(Aop.eigs(), \"k\", lw=2)\nplt.title(\"Eigenvalues\", size=16, fontweight=\"bold\")\nplt.xlabel(\"#eigenvalue\")\nplt.xlabel(\"intensity\")\nplt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also repeat the same exercise for a non-square matrix\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "N, M = 200, 50\nA = np.random.normal(0, 1, (N, M))\nx = np.ones(M)\n\nAop = pylops.MatrixMult(A, dtype=\"float64\")\ny = Aop * x\nyn = y + np.random.normal(0, 0.3, N)\n\nxest = Aop / y\nxnest = Aop / yn\n\nplt.figure(figsize=(8, 3))\nplt.plot(x, \"k\", lw=2, label=\"True\")\nplt.plot(xest, \"--r\", lw=2, label=\"Noise-free\")\nplt.plot(xnest, \"--g\", lw=2, label=\"Noisy\")\nplt.title(\"Matrix inversion\", size=16, fontweight=\"bold\")\nplt.legend()\nplt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And we can also use a sparse matrix from the :obj:`scipy.sparse`\nfamily of sparse matrices.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "N, M = 5, 5\nA = rand(N, M, density=0.75)\nx = np.ones(M)\n\nAop = pylops.MatrixMult(A, dtype=\"float64\")\ny = Aop * x\nxest = Aop / y\n\nprint(f\"A= {Aop.A.todense()}\")\nprint(f\"A^-1= {Aop.inv().todense()}\")\nprint(f\"eigs= {Aop.eigs()}\")\nprint(f\"x= {x}\")\nprint(f\"y= {y}\")\nprint(f\"lsqr solution xest= {xest}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, in several circumstances the input model $\\mathbf{x}$ may\nbe more naturally arranged as a matrix or a multi-dimensional array and\nit may be desirable to apply the same matrix to every columns of the model.\nThis can be mathematically expressed as:\n\n .. math::\n \\mathbf{y} =\n \\begin{bmatrix}\n \\mathbf{A} \\quad \\mathbf{0} \\quad \\mathbf{0} \\\\\n \\mathbf{0} \\quad \\mathbf{A} \\quad \\mathbf{0} \\\\\n \\mathbf{0} \\quad \\mathbf{0} \\quad \\mathbf{A}\n \\end{bmatrix}\n \\begin{bmatrix}\n \\mathbf{x_1} \\\\\n \\mathbf{x_2} \\\\\n \\mathbf{x_3}\n \\end{bmatrix}\n\nand apply it using the same implementation of the\n:py:class:`pylops.MatrixMult` operator by simply telling the operator how we\nwant the model to be organized through the ``otherdims`` input parameter.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "A = np.array([[1.0, 2.0], [4.0, 5.0]])\nx = np.array([[1.0, 1.0], [2.0, 2.0], [3.0, 3.0]])\n\nAop = pylops.MatrixMult(A, otherdims=(3,), dtype=\"float64\")\ny = Aop * x.ravel()\n\nxest, istop, itn, r1norm, r2norm = lsqr(Aop, y, damp=1e-10, iter_lim=10, show=0)[0:5]\nxest = xest.reshape(3, 2)\n\nprint(f\"A= {A}\")\nprint(f\"x= {x}\")\nprint(f\"y={y}\")\nprint(f\"lsqr solution xest= {xest}\")" ] } ], "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.9.15" } }, "nbformat": 4, "nbformat_minor": 0 }