docs/source/tutorials/linops/linear_operators_quickstart.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Operators Quickstart\n",
    "\n",
    "Finite-dimensional linear operators in ProbNum allow matrix algebra without explicitly constructing a full matrix representation. Instead it suffices to define a matrix-vector product and a shape attribute. This avoids unnecessary memory usage and generally yields a more efficient computation of matrix-vector products. \n",
    "\n",
    "ProbNum's linear operators also integrate with the `probnum.randvars` module and can be applied to random variable objects.\n",
    "\n",
    "This tutorial shows\n",
    "\n",
    "- how to construct a ProbNum linear operator from a dense NumPy and a sparse SciPy matrix.\n",
    "- how to define a custom matrix-free linear operator.\n",
    "- some of ProbNum's linear operator arithmetics.\n",
    "- how the operators are applied to ProbNum's random variable objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "# Make inline plots vector graphics instead of raster graphics\n",
    "%matplotlib inline\n",
    "from IPython.display import set_matplotlib_formats\n",
    "\n",
    "set_matplotlib_formats(\"pdf\", \"svg\")\n",
    "\n",
    "# Plotting\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.style.use(\"../../probnum.mplstyle\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Operators \n",
    "\n",
    "\n",
    "A finite-dimensional [linear operator](https://en.wikipedia.org/wiki/Linear_map) is a map between two finite-dimensional [vector spaces](https://en.wikipedia.org/wiki/Vector_space). Elements of the vector spaces can be represented as real-valued vectors $x\\in\\mathbb{R}^n$ and $b\\in\\mathbb{R}^m$ respectively, and the operator as a matrix $A\\in\\mathbb{R}^{m\\times n}$. Applying the operator to $x$ is then equal to performing a matrix-vector product $Ax = b$. \n",
    "\n",
    "For illustrative purposes we will use a simple example of a linear operator here: A permutation $P=A$ that moves each entry of a vector one spot forward and the last entry to the empty first spot. The matrix associated with that operator is square, has zeros everywhere and 1s on the subdiagonal and in the top right corner. \n",
    "\n",
    "We will create $P$ now as a NumPy `ndarray` for $n=m=5$ and apply it to a vector $x$. We observe that the permutation works as intended when we compute the matrix-vector product."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P:\n",
      " [[0. 0. 0. 0. 1.]\n",
      " [1. 0. 0. 0. 0.]\n",
      " [0. 1. 0. 0. 0.]\n",
      " [0. 0. 1. 0. 0.]\n",
      " [0. 0. 0. 1. 0.]]\n",
      "x:\n",
      " [0. 1. 2. 3. 4.]\n",
      "b=Px:\n",
      " [4. 0. 1. 2. 3.]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.sparse import diags\n",
    "\n",
    "n = 5  # size of the permulation\n",
    "P = diags([np.ones(n-1)], [-1]).toarray()\n",
    "P[0, -1] = 1\n",
    "\n",
    "x = np.arange(0., n, 1)\n",
    "b = P @ x  # apply the permutation operator\n",
    "\n",
    "print(\"P:\\n\", P)\n",
    "print(\"x:\\n\", x)\n",
    "print(\"b=Px:\\n\", b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating a ProbNum LinearOperator from a Dense Matrix.\n",
    "\n",
    "We will now create a ProbNum linear operator `P_op` from the dense NumPy array `P`. This is a naive but often useful way to create a linear operator in ProbNum. The functionality we are looking for is provided by the `Matrix` class. Even though `P_op` is an instance of `Matrix`, we can use the same notation (`@`) as above to apply the operator to the vector `x`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Matrix with shape=(5, 5) and dtype=float64>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from probnum.linops import Matrix\n",
    "\n",
    "P_op = Matrix(P)  \n",
    "P_op"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4., 0., 1., 2., 3.])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P_op @ x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use other arithmetic operations generally associated with matrices, such as adding them (`+`), transposing them (`.T`), or multiplying (`@`) them also with ProbNum's linear operator instances. Here are some examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1., 2., 3., 4., 0.])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# apply transpose of P to x\n",
    "P_op.T @ x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([5., 2., 4., 6., 3.])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# apply sum of P and P transpose to x\n",
    "(P_op + P_op.transpose()) @ x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([8., 0., 2., 4., 6.])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# scalar multiplication\n",
    "(2 * P_op) @ x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3., 4., 0., 1., 2.])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# applying P twice (this moved the elements 2 spots forward)\n",
    "(P_op @ P_op) @ x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When summing or multiplying linear operators, we get new linear operators, called `SumLinearOperator` and `ProductLinearOperator` that are built from their constituents. There are more operations that can be performed on linear operators not listed here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Matrix with shape=(5, 5) and dtype=float64>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P_op + P_op"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ProductLinearOperator [\n",
       "\t<Matrix with shape=(5, 5) and dtype=float64>, \n",
       "\t<Matrix with shape=(5, 5) and dtype=float64>, \n",
       "]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P_op @ P_op"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Efficient LinearOperators in ProbNum \n",
    "\n",
    "The `Matrix` class is quite useful as it allows for a quick creation of linear operators that have other benefits as well such as being compatible with ProbNum's random variable objects (see below). However, storing a dense matrix is not always feasible or desired, especially when the dimensionality of the matrix is large. \n",
    "\n",
    "ProbNum thus currently supports two alternative ways to create linear maps: i) The `Matrix` class is also compatible with SciPy's sparse matrices, and ii) custom implementations of matrix-free linear operations.  \n",
    "\n",
    "Both alternatives have two benefits: 1) they are memory-efficient, and ii) they generally yield efficient matrix-vector computations, reducing the trivial $\\mathcal{O}(nm)$ ($\\mathcal{O}(n^2)$ for square matrices) complexity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sparse LinearOperators \n",
    "\n",
    "Instead of dense matrices, `Matrix` can be used with SciPy's `sparse` matrices as well. The interface is analogous to above, we simply hand the sparse matrix to `Matrix` instead of the 2D NumPy array. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.49421769, -1.35451937,  1.05551543, -0.41823967,  0.42934955,\n",
       "       -0.82155968, -1.93141743, -4.31860989, -1.70475714,  4.36385187,\n",
       "        2.36850628, -2.94034717,  0.39821307, -1.08656905,  0.36490375,\n",
       "       -0.86441656, -0.44778464, -0.44155178,  0.55687361,  0.17178464])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import scipy.sparse\n",
    "\n",
    "# Create a random sparse matrix using SciPy\n",
    "n = 20\n",
    "A_scipy = scipy.sparse.rand(m=n, n=n, density=0.05, random_state=42)\n",
    "\n",
    "# create a ProbNum linear operator\n",
    "A_op = Matrix(A=A_scipy)\n",
    "\n",
    "# Some linear operator arithmetic\n",
    "from probnum.linops import Identity\n",
    "x = np.random.randn(n)\n",
    "Id = Identity(shape=n)\n",
    "(A_op + 1.5 * Id) @ x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matrix-Free LinearOperators \n",
    "\n",
    "Now we create an efficient, matrix-free version of a linear operator. \n",
    "\n",
    "In practice, it is often only necessary to know the result of a linear operator applied to an arbitrary vector (a function $b(x) = Ax$), but not the explicit matrix-representation of the linear operator itself. This is an implicit definition that does not require the storage of the $nm$ ($n^2$ if $n=m$) entries of the matrix form of $A$ in contrast to the explicit definition above.\n",
    "\n",
    "In our example (simple permutation), applying `P` to `x` is conveniently performed by NumPy's `roll` functionality; but any custom map could be used instead.\n",
    "First, we define the function handle $b(x)=Ax$ that represents the desired matrix-vector product (called `mv` below). Then, we hand it to ProbNum's general `LinearOperator` class to create an instance of the implicitly defined linear operator. The computation of the matrix-vector product is now as efficient as NumPy's `roll` function, and the matrix $A$ needs not be stored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<LambdaLinearOperator with shape=(5, 5) and dtype=float64>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from probnum.linops import LinearOperator, LambdaLinearOperator\n",
    "\n",
    "@LinearOperator.broadcast_matvec\n",
    "def mv(v):\n",
    "    return np.roll(v, 1)\n",
    "\n",
    "n = 5\n",
    "P_op = LambdaLinearOperator(shape=(n, n), dtype=np.float_, matmul=mv)\n",
    "x = np.arange(0., n, 1)\n",
    "\n",
    "P_op"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4., 0., 1., 2., 3.])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P_op @ x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can still create a dense matrix, even from the implicit definition if we really need to. However, most of the time this is not advised as the matrix may be very large. Other functionality, such as computing the determinant, rank or eigenvalues of a linear operator currently may also require to construct a dense matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0., 0., 0., 0., 1.],\n",
       "       [1., 0., 0., 0., 0.],\n",
       "       [0., 1., 0., 0., 0.],\n",
       "       [0., 0., 1., 0., 0.],\n",
       "       [0., 0., 0., 1., 0.]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P_op.todense()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Operators and Random Variables\n",
    "\n",
    "ProbNum's linear operators integrate with `probnum.randvars` which means that linear operators can be applied to random variable objects to obtain a transformed random variable object. This is already shown the the [Random Variables Quickstart](https://probnum.readthedocs.io/en/latest/tutorials/prob/random_variables_quickstart.html). Hence, here we only show a simple example that computes the marginal of a 3D [normal random variable](https://en.wikipedia.org/wiki/Normal_distribution) (marginalizing out the 3rd, and keeping the first 2 dimensions). For normal random variables this is equivalent to applying a projection operator `Pr` that projects on the the first two axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "from probnum.randvars import Normal\n",
    "\n",
    "# Create the 3D normal random variable\n",
    "n = 3  \n",
    "rv = Normal(mean = np.arange(n, 0, -1), cov = np.diag(np.arange(1, n+1, 1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3., 2., 1.])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rv.mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0., 0.],\n",
       "       [0., 2., 0.],\n",
       "       [0., 0., 3.]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rv.cov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the projection operator\n",
    "@LinearOperator.broadcast_matvec\n",
    "def mv(v):\n",
    "    return v[:n-1]\n",
    "\n",
    "Pr = LambdaLinearOperator(shape=(n-1, n), dtype=np.float_, matmul=mv)\n",
    "\n",
    "# Apply the operator to the 3D normal random variable\n",
    "rv_projected = Pr @ rv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Normal with shape=(2,), dtype=float64>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# the result is a 2D normal random variable\n",
    "rv_projected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3., 2.])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rv_projected.mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0.],\n",
       "       [0., 2.]])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rv_projected.cov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0., 0.],\n",
       "       [0., 1., 0.]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Dense 2x3 projection operator for completeness\n",
    "Pr.todense()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.6 (conda)",
   "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.10.6"
  },
  "vscode": {
   "interpreter": {
    "hash": "0457b12441837086dec1b475e0008c28e5fc37f4ffe0e5ee9f2b481cc28bc3c9"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}