{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "3c094389",
            "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": "a4d75de8",
            "cell_type": "markdown",
            "source": "Construct a matrix `A` from singular values and vectors.\n",
            "metadata": {}
        },
        {
            "id": "c91670f6",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import numpy as np\nfrom scipy import sparse, linalg, stats\nfrom scipy.sparse.linalg import svds, aslinearoperator, LinearOperator",
            "outputs": []
        },
        {
            "id": "50eefed5",
            "cell_type": "markdown",
            "source": "Construct a dense matrix `A` from singular values and vectors.\n",
            "metadata": {}
        },
        {
            "id": "f3524db4",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rng = np.random.default_rng()\northogonal = stats.ortho_group.rvs(10, random_state=rng)\ns = [1e-3, 1, 2, 3, 4]  # non-zero singular values\nu = orthogonal[:, :5]         # left singular vectors\nvT = orthogonal[:, 5:].T      # right singular vectors\nA = u @ np.diag(s) @ vT",
            "outputs": []
        },
        {
            "id": "53c8b5a4",
            "cell_type": "markdown",
            "source": "With only four singular values/vectors, the SVD approximates the original\nmatrix.\n",
            "metadata": {}
        },
        {
            "id": "e6bd5287",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "u4, s4, vT4 = svds(A, k=4)\nA4 = u4 @ np.diag(s4) @ vT4\nnp.allclose(A4, A, atol=1e-3)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "1c51892d",
            "cell_type": "markdown",
            "source": "With all five non-zero singular values/vectors, we can reproduce\nthe original matrix more accurately.\n",
            "metadata": {}
        },
        {
            "id": "91cdeb9d",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "u5, s5, vT5 = svds(A, k=5)\nA5 = u5 @ np.diag(s5) @ vT5\nnp.allclose(A5, A)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "91b3601d",
            "cell_type": "markdown",
            "source": "The singular values match the expected singular values.\n",
            "metadata": {}
        },
        {
            "id": "06616e3f",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "np.allclose(s5, s)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "9d19bcbd",
            "cell_type": "markdown",
            "source": "Since the singular values are not close to each other in this example,\nevery singular vector matches as expected up to a difference in sign.\n",
            "metadata": {}
        },
        {
            "id": "9d3ff8a7",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "(np.allclose(np.abs(u5), np.abs(u)) and\n np.allclose(np.abs(vT5), np.abs(vT)))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "302b124e",
            "cell_type": "markdown",
            "source": "The singular vectors are also orthogonal.\n",
            "metadata": {}
        },
        {
            "id": "4ecf949c",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "(np.allclose(u5.T @ u5, np.eye(5)) and\n np.allclose(vT5 @ vT5.T, np.eye(5)))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "a673b974",
            "cell_type": "markdown",
            "source": "If there are (nearly) multiple singular values, the corresponding\nindividual singular vectors may be unstable, but the whole invariant\nsubspace containing all such singular vectors is computed accurately\nas can be measured by angles between subspaces via 'subspace_angles'.\n",
            "metadata": {}
        },
        {
            "id": "3d7942d8",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rng = np.random.default_rng()\ns = [1, 1 + 1e-6]  # non-zero singular values\nu, _ = np.linalg.qr(rng.standard_normal((99, 2)))\nv, _ = np.linalg.qr(rng.standard_normal((99, 2)))\nvT = v.T\nA = u @ np.diag(s) @ vT\nA = A.astype(np.float32)\nu2, s2, vT2 = svds(A, k=2, rng=rng)\nnp.allclose(s2, s)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "2028cf0d",
            "cell_type": "markdown",
            "source": "The angles between the individual exact and computed singular vectors\nmay not be so small. To check use:\n",
            "metadata": {}
        },
        {
            "id": "c07bca6d",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "(linalg.subspace_angles(u2[:, :1], u[:, :1]) +\n linalg.subspace_angles(u2[:, 1:], u[:, 1:]))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([0.06562513])  # may vary"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "dfcc0a3f",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "(linalg.subspace_angles(vT2[:1, :].T, vT[:1, :].T) +\n linalg.subspace_angles(vT2[1:, :].T, vT[1:, :].T))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([0.06562507])  # may vary"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "56e37cf6",
            "cell_type": "markdown",
            "source": "As opposed to the angles between the 2-dimensional invariant subspaces\nthat these vectors span, which are small for rights singular vectors\n",
            "metadata": {}
        },
        {
            "id": "f5b1dab4",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "linalg.subspace_angles(u2, u).sum() < 1e-6",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "bc931335",
            "cell_type": "markdown",
            "source": "as well as for left singular vectors.\n",
            "metadata": {}
        },
        {
            "id": "96af1d7e",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "linalg.subspace_angles(vT2.T, vT.T).sum() < 1e-6",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "1293d7ab",
            "cell_type": "markdown",
            "source": "The next example follows that of 'sklearn.decomposition.TruncatedSVD'.\n",
            "metadata": {}
        },
        {
            "id": "2858c64f",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rng = np.random.default_rng()\nX_dense = rng.random(size=(100, 100))\nX_dense[:, 2 * np.arange(50)] = 0\nX = sparse.csr_array(X_dense)\n_, singular_values, _ = svds(X, k=5, rng=rng)\nprint(singular_values)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[ 4.3221...  4.4043...  4.4907...  4.5858... 35.4549...]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "9303bab1",
            "cell_type": "markdown",
            "source": "The function can be called without the transpose of the input matrix\never explicitly constructed.\n",
            "metadata": {}
        },
        {
            "id": "3eb25c69",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rng = np.random.default_rng()\nG = sparse.random_array((8, 9), density=0.5, rng=rng)\nGlo = aslinearoperator(G)\n_, singular_values_svds, _ = svds(Glo, k=5, rng=rng)\n_, singular_values_svd, _ = linalg.svd(G.toarray())\nnp.allclose(singular_values_svds, singular_values_svd[-4::-1])",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "6ea2da53",
            "cell_type": "markdown",
            "source": "The most memory efficient scenario is where neither\nthe original matrix, nor its transpose, is explicitly constructed.\nOur example computes the smallest singular values and vectors\nof 'LinearOperator' constructed from the numpy function 'np.diff' used\ncolumn-wise to be consistent with 'LinearOperator' operating on columns.\n",
            "metadata": {}
        },
        {
            "id": "b017f78f",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "diff0 = lambda a: np.diff(a, axis=0)",
            "outputs": []
        },
        {
            "id": "bc1cd40e",
            "cell_type": "markdown",
            "source": "Let us create the matrix from 'diff0' to be used for validation only.\n",
            "metadata": {}
        },
        {
            "id": "f5a4b2ac",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "n = 5  # The dimension of the space.\nM_from_diff0 = diff0(np.eye(n))\nprint(M_from_diff0.astype(int))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[[-1  1  0  0  0]\n [ 0 -1  1  0  0]\n [ 0  0 -1  1  0]\n [ 0  0  0 -1  1]]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "2739be3b",
            "cell_type": "markdown",
            "source": "The matrix 'M_from_diff0' is bi-diagonal and could be alternatively\ncreated directly by\n",
            "metadata": {}
        },
        {
            "id": "e6141690",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "M = - np.eye(n - 1, n, dtype=int)\nnp.fill_diagonal(M[:,1:], 1)\nnp.allclose(M, M_from_diff0)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "76e14367",
            "cell_type": "markdown",
            "source": "Its transpose\n",
            "metadata": {}
        },
        {
            "id": "bcb5b670",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "print(M.T)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[[-1  0  0  0]\n [ 1 -1  0  0]\n [ 0  1 -1  0]\n [ 0  0  1 -1]\n [ 0  0  0  1]]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "c1fa0b66",
            "cell_type": "markdown",
            "source": "can be viewed as the incidence matrix; see\nIncidence matrix, (2022, Nov. 19), Wikipedia, https://w.wiki/5YXU,\nof a linear graph with 5 vertices and 4 edges. The 5x5 normal matrix\n``M.T @ M`` thus is\n",
            "metadata": {}
        },
        {
            "id": "b32404ff",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "print(M.T @ M)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[[ 1 -1  0  0  0]\n [-1  2 -1  0  0]\n [ 0 -1  2 -1  0]\n [ 0  0 -1  2 -1]\n [ 0  0  0 -1  1]]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "f1fdc569",
            "cell_type": "markdown",
            "source": "the graph Laplacian, while the actually used in 'svds' smaller size\n4x4 normal matrix ``M @ M.T``\n",
            "metadata": {}
        },
        {
            "id": "44cb3f0c",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "print(M @ M.T)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[[ 2 -1  0  0]\n [-1  2 -1  0]\n [ 0 -1  2 -1]\n [ 0  0 -1  2]]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "0367f0b3",
            "cell_type": "markdown",
            "source": "is the so-called edge-based Laplacian; see\nSymmetric Laplacian via the incidence matrix, in Laplacian matrix,\n(2022, Nov. 19), Wikipedia, https://w.wiki/5YXW.\n\nThe 'LinearOperator' setup needs the options 'rmatvec' and 'rmatmat'\nof multiplication by the matrix transpose ``M.T``, but we want to be\nmatrix-free to save memory, so knowing how ``M.T`` looks like, we\nmanually construct the following function to be\nused in ``rmatmat=diff0t``.\n",
            "metadata": {}
        },
        {
            "id": "8534e137",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def diff0t(a):\n    if a.ndim == 1:\n        a = a[:,np.newaxis]  # Turn 1D into 2D array\n    d = np.zeros((a.shape[0] + 1, a.shape[1]), dtype=a.dtype)\n    d[0, :] = - a[0, :]\n    d[1:-1, :] = a[0:-1, :] - a[1:, :]\n    d[-1, :] = a[-1, :]\n    return d",
            "outputs": []
        },
        {
            "id": "4fd9850d",
            "cell_type": "markdown",
            "source": "We check that our function 'diff0t' for the matrix transpose is valid.\n",
            "metadata": {}
        },
        {
            "id": "9e3d1ab3",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "np.allclose(M.T, diff0t(np.eye(n-1)))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "87bbe696",
            "cell_type": "markdown",
            "source": "Now we setup our matrix-free 'LinearOperator' called 'diff0_func_aslo'\nand for validation the matrix-based 'diff0_matrix_aslo'.\n",
            "metadata": {}
        },
        {
            "id": "26892ed4",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def diff0_func_aslo_def(n):\n    return LinearOperator(matvec=diff0,\n                          matmat=diff0,\n                          rmatvec=diff0t,\n                          rmatmat=diff0t,\n                          shape=(n - 1, n))\ndiff0_func_aslo = diff0_func_aslo_def(n)\ndiff0_matrix_aslo = aslinearoperator(M_from_diff0)",
            "outputs": []
        },
        {
            "id": "85445927",
            "cell_type": "markdown",
            "source": "And validate both the matrix and its transpose in 'LinearOperator'.\n",
            "metadata": {}
        },
        {
            "id": "0efa015c",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "np.allclose(diff0_func_aslo(np.eye(n)),\n            diff0_matrix_aslo(np.eye(n)))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "7767712f",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "np.allclose(diff0_func_aslo.T(np.eye(n-1)),\n            diff0_matrix_aslo.T(np.eye(n-1)))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "4ee01465",
            "cell_type": "markdown",
            "source": "Having the 'LinearOperator' setup validated, we run the solver.\n",
            "metadata": {}
        },
        {
            "id": "b27bef95",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "n = 100\ndiff0_func_aslo = diff0_func_aslo_def(n)\nu, s, vT = svds(diff0_func_aslo, k=3, which='SM')",
            "outputs": []
        },
        {
            "id": "179351d3",
            "cell_type": "markdown",
            "source": "The singular values squared and the singular vectors are known\nexplicitly; see\nPure Dirichlet boundary conditions, in\nEigenvalues and eigenvectors of the second derivative,\n(2022, Nov. 19), Wikipedia, https://w.wiki/5YX6,\nsince 'diff' corresponds to first\nderivative, and its smaller size n-1 x n-1 normal matrix\n``M @ M.T`` represent the discrete second derivative with the Dirichlet\nboundary conditions. We use these analytic expressions for validation.\n",
            "metadata": {}
        },
        {
            "id": "074d009b",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "se = 2. * np.sin(np.pi * np.arange(1, 4) / (2. * n))\nue = np.sqrt(2 / n) * np.sin(np.pi * np.outer(np.arange(1, n),\n                             np.arange(1, 4)) / n)\nnp.allclose(s, se, atol=1e-3)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "7ecd39cd",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "np.allclose(np.abs(u), np.abs(ue), atol=1e-6)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        }
    ]
}