Linear Operators Quickstart

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.

ProbNum’s linear operators also integrate with the probnum.randvars module and can be applied to random variable objects.

This tutorial shows

  • how to construct a ProbNum linear operator from a dense NumPy and a sparse SciPy matrix.

  • how to define a custom matrix-free linear operator.

  • some of ProbNum’s linear operator arithmetics.

  • how the operators are applied to ProbNum’s random variable objects.

import warnings

# Make inline plots vector graphics instead of raster graphics
%matplotlib inline
from IPython.display import set_matplotlib_formats

set_matplotlib_formats("pdf", "svg")

# Plotting
import matplotlib.pyplot as plt"../../probnum.mplstyle")

Linear Operators

A finite-dimensional linear operator is a map between two finite-dimensional vector spaces. 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\).

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.

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.

import numpy as np
from scipy.sparse import diags

n = 5  # size of the permulation
P = diags([np.ones(n-1)], [-1]).toarray()
P[0, -1] = 1

x = np.arange(0., n, 1)
b = P @ x  # apply the permutation operator

print("P:\n", P)
print("x:\n", x)
print("b=Px:\n", b)
 [[0. 0. 0. 0. 1.]
 [1. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0.]
 [0. 0. 1. 0. 0.]
 [0. 0. 0. 1. 0.]]
 [0. 1. 2. 3. 4.]
 [4. 0. 1. 2. 3.]

Creating a ProbNum LinearOperator from a Dense Matrix.

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.

from probnum.linops import Matrix

P_op = Matrix(P)
<Matrix with shape=(5, 5) and dtype=float64>
P_op @ x
array([4., 0., 1., 2., 3.])

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:

# apply transpose of P to x
P_op.T @ x
array([1., 2., 3., 4., 0.])
# apply sum of P and P transpose to x
(P_op + P_op.transpose()) @ x
array([5., 2., 4., 6., 3.])
# scalar multiplication
(2 * P_op) @ x
array([8., 0., 2., 4., 6.])
# applying P twice (this moved the elements 2 spots forward)
(P_op @ P_op) @ x
array([3., 4., 0., 1., 2.])

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.

P_op + P_op
<Matrix with shape=(5, 5) and dtype=float64>
P_op @ P_op
ProductLinearOperator [
        <Matrix with shape=(5, 5) and dtype=float64>,
        <Matrix with shape=(5, 5) and dtype=float64>,

Efficient LinearOperators in ProbNum

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.

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.

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.

Sparse LinearOperators

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.

import scipy.sparse

# Create a random sparse matrix using SciPy
n = 20
A_scipy = scipy.sparse.rand(m=n, n=n, density=0.05, random_state=42)

# create a ProbNum linear operator
A_op = Matrix(A=A_scipy)

# Some linear operator arithmetic
from probnum.linops import Identity
x = np.random.randn(n)
Id = Identity(shape=n)
(A_op + 1.5 * Id) @ x
array([ 1.49421769, -1.35451937,  1.05551543, -0.41823967,  0.42934955,
       -0.82155968, -1.93141743, -4.31860989, -1.70475714,  4.36385187,
        2.36850628, -2.94034717,  0.39821307, -1.08656905,  0.36490375,
       -0.86441656, -0.44778464, -0.44155178,  0.55687361,  0.17178464])

Matrix-Free LinearOperators

Now we create an efficient, matrix-free version of a linear operator.

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.

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. 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.

from probnum.linops import LinearOperator, LambdaLinearOperator

def mv(v):
    return np.roll(v, 1)

n = 5
P_op = LambdaLinearOperator(shape=(n, n), dtype=np.float_, matmul=mv)
x = np.arange(0., n, 1)

<LambdaLinearOperator with shape=(5, 5) and dtype=float64>
P_op @ x
array([4., 0., 1., 2., 3.])

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.

array([[0., 0., 0., 0., 1.],
       [1., 0., 0., 0., 0.],
       [0., 1., 0., 0., 0.],
       [0., 0., 1., 0., 0.],
       [0., 0., 0., 1., 0.]])

Linear Operators and Random Variables

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. Hence, here we only show a simple example that computes the marginal of a 3D normal random variable (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.

from probnum.randvars import Normal

# Create the 3D normal random variable
n = 3
rv = Normal(mean = np.arange(n, 0, -1), cov = np.diag(np.arange(1, n+1, 1)))
array([3., 2., 1.])
array([[1., 0., 0.],
       [0., 2., 0.],
       [0., 0., 3.]])
# Define the projection operator
def mv(v):
    return v[:n-1]

Pr = LambdaLinearOperator(shape=(n-1, n), dtype=np.float_, matmul=mv)

# Apply the operator to the 3D normal random variable
rv_projected = Pr @ rv
# the result is a 2D normal random variable
<Normal with shape=(2,), dtype=float64>
array([3., 2.])
array([[1., 0.],
       [0., 2.]])
# Dense 2x3 projection operator for completeness
array([[1., 0., 0.],
       [0., 1., 0.]])