{ "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\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import warnings\nimport numpy as np\nfrom scipy.sparse import rand\nfrom scipy.sparse.linalg import lsqr\n\nimport matplotlib.pyplot as plt\nimport matplotlib.gridspec as pltgs\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(0.35, 0.5, '=', horizontalalignment='center',\n verticalalignment='center', size=40, fontweight='bold')\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([])\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(0.35, 0.5, '=', horizontalalignment='center',\n verticalalignment='center', size=40, fontweight='bold')\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([])" ] }, { "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()" ] }, { "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('A= %s' % Aop.A.todense())\nprint('A^-1=', Aop.inv().todense())\nprint('eigs=', Aop.eigs())\nprint('x= %s' % x)\nprint('y= %s' % y)\nprint('lsqr solution xest= %s' % 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 ``dims`` input parameter.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "A = np.array([[1., 2.], [4., 5.]])\nx = np.array([[1., 1.],\n [2., 2.],\n [3., 3.]])\n\nAop = pylops.MatrixMult(A, dims=(3,), dtype='float64')\ny = Aop*x.flatten()\n\nxest, istop, itn, r1norm, r2norm = \\\n lsqr(Aop, y, damp=1e-10, iter_lim=10, show=0)[0:5]\nxest = xest.reshape(3, 2)\n\nprint('A= %s' % A)\nprint('x= %s' % x)\nprint('y= %s' % y)\nprint('lsqr solution xest= %s' % 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.6.12" } }, "nbformat": 4, "nbformat_minor": 0 }