{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "0c4f2eeb",
            "cell_type": "markdown",
            "source": "<div class='alert alert-warning'>\n\nSciPy's interactive examples with Jupyterlite are experimental and may not always work as expected. Execution of cells containing imports may result in large downloads (up to 60MB of content for the first import from SciPy). Load times when importing from SciPy may take roughly 10-20 seconds. If you notice any problems, feel free to open an [issue](https://github.com/scipy/scipy/issues/new/choose).\n\n</div>",
            "metadata": {}
        },
        {
            "id": "3f9005b2",
            "cell_type": "markdown",
            "source": "Our first example is minimalistic - find the largest eigenvalue of\na diagonal matrix by solving the non-generalized eigenvalue problem\n``A x = lambda x`` without constraints or preconditioning.\n",
            "metadata": {}
        },
        {
            "id": "d1bd24ff",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import numpy as np\nfrom scipy.sparse import diags_array\nfrom scipy.sparse.linalg import LinearOperator, aslinearoperator\nfrom scipy.sparse.linalg import lobpcg",
            "outputs": []
        },
        {
            "id": "ba50ef5e",
            "cell_type": "markdown",
            "source": "The square matrix size is\n",
            "metadata": {}
        },
        {
            "id": "1fe171e8",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "n = 100",
            "outputs": []
        },
        {
            "id": "fa1a4895",
            "cell_type": "markdown",
            "source": "and its diagonal entries are 1, ..., 100 defined by\n",
            "metadata": {}
        },
        {
            "id": "c7d43a54",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "vals = np.arange(1, n + 1).astype(np.int16)",
            "outputs": []
        },
        {
            "id": "00ad4355",
            "cell_type": "markdown",
            "source": "The first mandatory input parameter in this test is\nthe sparse diagonal matrix `A`\nof the eigenvalue problem ``A x = lambda x`` to solve.\n",
            "metadata": {}
        },
        {
            "id": "bf8fcfe8",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "A = diags_array(vals, offsets=0, shape=(n, n))\nA = A.astype(np.int16)\nA.toarray()",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[  1,   0,   0, ...,   0,   0,   0],\n       [  0,   2,   0, ...,   0,   0,   0],\n       [  0,   0,   3, ...,   0,   0,   0],\n       ...,\n       [  0,   0,   0, ...,  98,   0,   0],\n       [  0,   0,   0, ...,   0,  99,   0],\n       [  0,   0,   0, ...,   0,   0, 100]], shape=(100, 100), dtype=int16)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "9233b848",
            "cell_type": "markdown",
            "source": "The second mandatory input parameter `X` is a 2D array with the\nrow dimension determining the number of requested eigenvalues.\n`X` is an initial guess for targeted eigenvectors.\n`X` must have linearly independent columns.\nIf no initial approximations available, randomly oriented vectors\ncommonly work best, e.g., with components normally distributed\naround zero or uniformly distributed on the interval [-1 1].\nSetting the initial approximations to dtype ``np.float32``\nforces all iterative values to dtype ``np.float32`` speeding up\nthe run while still allowing accurate eigenvalue computations.\n",
            "metadata": {}
        },
        {
            "id": "17c3b306",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "k = 1\nrng = np.random.default_rng()\nX = rng.normal(size=(n, k))\nX = X.astype(np.float32)",
            "outputs": []
        },
        {
            "id": "1d94223b",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "eigenvalues, _ = lobpcg(A, X, maxiter=60)\neigenvalues",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([100.], dtype=float32)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "db8f5e44",
            "cell_type": "markdown",
            "source": "`lobpcg` needs only access the matrix product with `A` rather\nthen the matrix itself. Since the matrix `A` is diagonal in\nthis example, one can write a function of the matrix product\n``A @ X`` using the diagonal values ``vals`` only, e.g., by\nelement-wise multiplication with broadcasting in the lambda-function\n",
            "metadata": {}
        },
        {
            "id": "9475e536",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "A_lambda = lambda X: vals[:, np.newaxis] * X",
            "outputs": []
        },
        {
            "id": "4199598e",
            "cell_type": "markdown",
            "source": "or the regular function\n",
            "metadata": {}
        },
        {
            "id": "878d9c57",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def A_matmat(X):\n    return vals[:, np.newaxis] * X",
            "outputs": []
        },
        {
            "id": "1d7afbb8",
            "cell_type": "markdown",
            "source": "and use the handle to one of these callables as an input\n",
            "metadata": {}
        },
        {
            "id": "482903d9",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "eigenvalues, _ = lobpcg(A_lambda, X, maxiter=60)\neigenvalues",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([100.], dtype=float32)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "cbdad926",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "eigenvalues, _ = lobpcg(A_matmat, X, maxiter=60)\neigenvalues",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([100.], dtype=float32)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "debb9bca",
            "cell_type": "markdown",
            "source": "The traditional callable `LinearOperator` is no longer\nnecessary but still supported as the input to `lobpcg`.\nSpecifying ``matmat=A_matmat`` explicitly improves performance. \n",
            "metadata": {}
        },
        {
            "id": "2a5ac157",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "A_lo = LinearOperator((n, n), matvec=A_matmat, matmat=A_matmat, dtype=np.int16)\neigenvalues, _ = lobpcg(A_lo, X, maxiter=80)\neigenvalues",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([100.], dtype=float32)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "c7582801",
            "cell_type": "markdown",
            "source": "The least efficient callable option is `aslinearoperator`:\n",
            "metadata": {}
        },
        {
            "id": "9e75473a",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "eigenvalues, _ = lobpcg(aslinearoperator(A), X, maxiter=80)\neigenvalues",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([100.], dtype=float32)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "5990a5f2",
            "cell_type": "markdown",
            "source": "We now switch to computing the three smallest eigenvalues specifying\n",
            "metadata": {}
        },
        {
            "id": "b35cf17a",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "k = 3\nX = np.random.default_rng().normal(size=(n, k))",
            "outputs": []
        },
        {
            "id": "c17af5d0",
            "cell_type": "markdown",
            "source": "and ``largest=False`` parameter\n",
            "metadata": {}
        },
        {
            "id": "6ebe28ef",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "eigenvalues, _ = lobpcg(A, X, largest=False, maxiter=90)\nprint(eigenvalues)  ",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[1. 2. 3.]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "875011d9",
            "cell_type": "markdown",
            "source": "The next example illustrates computing 3 smallest eigenvalues of\nthe same matrix `A` given by the function handle ``A_matmat`` but\nwith constraints and preconditioning.\n\nConstraints - an optional input parameter is a 2D array comprising\nof column vectors that the eigenvectors must be orthogonal to\n",
            "metadata": {}
        },
        {
            "id": "69510d0c",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "Y = np.eye(n, 3)",
            "outputs": []
        },
        {
            "id": "4484704f",
            "cell_type": "markdown",
            "source": "The preconditioner acts as the inverse of `A` in this example, but\nin the reduced precision ``np.float32`` even though the initial `X`\nand thus all iterates and the output are in full ``np.float64``.\n",
            "metadata": {}
        },
        {
            "id": "cad0e552",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "inv_vals = 1./vals\ninv_vals = inv_vals.astype(np.float32)\nM = lambda X: inv_vals[:, np.newaxis] * X",
            "outputs": []
        },
        {
            "id": "cf1dbac6",
            "cell_type": "markdown",
            "source": "Let us now solve the eigenvalue problem for the matrix `A` first\nwithout preconditioning requesting 80 iterations\n",
            "metadata": {}
        },
        {
            "id": "0893bdce",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "eigenvalues, _ = lobpcg(A_matmat, X, Y=Y, largest=False, maxiter=80)\neigenvalues",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([4., 5., 6.])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "24d3a872",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "eigenvalues.dtype",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "dtype('float64')"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "5d2b5b34",
            "cell_type": "markdown",
            "source": "With preconditioning we need only 20 iterations from the same `X`\n",
            "metadata": {}
        },
        {
            "id": "51220525",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "eigenvalues, _ = lobpcg(A_matmat, X, Y=Y, M=M, largest=False, maxiter=20)\neigenvalues",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([4., 5., 6.])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "ae9c8e6b",
            "cell_type": "markdown",
            "source": "Note that the vectors passed in `Y` are the eigenvectors of the 3\nsmallest eigenvalues. The results returned above are orthogonal to those.\n\nThe primary matrix `A` may be indefinite, e.g., after shifting\n``vals`` by 50 from 1, ..., 100 to -49, ..., 50, we still can compute\nthe 3 smallest or largest eigenvalues.\n",
            "metadata": {}
        },
        {
            "id": "6eb5a4a5",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "vals = vals - 50\nX = rng.normal(size=(n, k))\neigenvalues, _ = lobpcg(A_matmat, X, largest=False, maxiter=99)\neigenvalues",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([-49., -48., -47.])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "16dfac56",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "eigenvalues, _ = lobpcg(A_matmat, X, largest=True, maxiter=99)\neigenvalues",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([50., 49., 48.])"
                    },
                    "execution_count": null
                }
            ]
        }
    ]
}