{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# CGLS and LSQR Solvers\n\nThis example shows how to use the :py:func:`pylops.optimization.leastsquares.cgls`\nand :py:func:`pylops.optimization.leastsquares.lsqr` PyLops solvers\nto minimize the following cost function:\n\n\\begin{align}J = \\| \\mathbf{y} -  \\mathbf{Ax} \\|_2^2 + \\epsilon \\| \\mathbf{x} \\|_2^2\\end{align}\n\nNote that the LSQR solver behaves in the same way as the scipy's\n:py:func:`scipy.sparse.linalg.lsqr` solver. However, our solver is also able\nto operate on cupy arrays and perform computations on a GPU.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import warnings\n\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nimport pylops\n\nplt.close(\"all\")\nwarnings.filterwarnings(\"ignore\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's define a matrix $\\mathbf{A}$ or size (``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, 10\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 use the cgls solver to invert this matrix\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "y = Aop * x\nxest, istop, nit, r1norm, r2norm, cost_cgls = pylops.optimization.basic.cgls(\n    Aop, y, x0=np.zeros_like(x), niter=10, tol=1e-10, show=True\n)\n\nprint(f\"x= {x}\")\nprint(f\"cgls solution xest= {xest}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And the lsqr solver to invert this matrix\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "y = Aop * x\n(\n    xest,\n    istop,\n    itn,\n    r1norm,\n    r2norm,\n    anorm,\n    acond,\n    arnorm,\n    xnorm,\n    var,\n    cost_lsqr,\n) = pylops.optimization.basic.lsqr(Aop, y, x0=np.zeros_like(x), niter=10, show=True)\n\nprint(f\"x= {x}\")\nprint(f\"lsqr solution xest= {xest}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally we show that the L2 norm of the residual of the two solvers decays\nin the same way, as LSQR is algebrically equivalent to CG on the normal\nequations and CGLS\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure(figsize=(12, 3))\nplt.plot(cost_cgls, \"k\", lw=2, label=\"CGLS\")\nplt.plot(cost_lsqr, \"--r\", lw=2, label=\"LSQR\")\nplt.title(\"Cost functions\")\nplt.legend()\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Note that while we used a dense matrix here, any other linear operator\ncan be fed to cgls and lsqr as is the case for any other PyLops solver.\n\n"
      ]
    }
  ],
  "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
}