Linear Operators

Finite dimensional linear operators 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 can often be more convenient to derive.

[1]:
# 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

plt.style.use("../../probnum.mplstyle")

Sparsity and Operator Arithmetic

We begin by sampling a random sparse matrix, represent it as a LinearOperator and perform some simple arithmetic.

[2]:
import numpy as np
import scipy.sparse
from probnum import linops

# Linear operator from sparse matrix
n = 50
mat = scipy.sparse.rand(m=n, n=n, density=0.05, random_state=42)
A = linops.Matrix(A=mat)

# Linear operator arithmetic
Id = linops.Identity(shape=n)
B = A + 1.5 * Id
[3]:
# Plot
matdict = {"$A$": A.todense(), "$B$": B.todense()}
vmin = np.min([np.min(mat) for mat in list(matdict.values())])
vmax = np.max([np.max(mat) for mat in list(matdict.values())])

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(6, 3.5))
for i, (title, mat) in enumerate(matdict.items()):
    axes[i].imshow(mat, vmin=vmin, vmax=vmax)
    axes[i].set_axis_off()
    axes[i].title.set_text(title)
fig.tight_layout()
../../_images/tutorials_linalg_linear_operators_4_0.svg

The Kronecker Product

A frequently occuring linear operator is the Kronecker product of two matrices. Example application areas are matrix-variate distributions, signal processing, image processing, semi-definite programming and deep learning. We will illustrate some basic properties here.

[4]:
# Kronecker Product
A = np.array([[4, 1, 4], [2, 3, 2]])
B = np.array([[1], [4]])
W = linops.Kronecker(A, B)
V = linops.Kronecker(B, A)
[5]:
# Plot
matdict = {
    "$A$": A,
    "$B$": B,
    "$A \otimes B$": W.todense(),
    "$B \otimes A$": V.todense(),
}
vmin = np.min([np.min(mat) for mat in list(matdict.values())])
vmax = np.max([np.max(mat) for mat in list(matdict.values())])

fig, axes = plt.subplots(
    nrows=1, ncols=4, figsize=(6, 3.5), gridspec_kw={"width_ratios": [3, 1, 3, 3]}
)
for i, (title, mat) in enumerate(matdict.items()):
    axes[i].imshow(mat, vmin=vmin, vmax=vmax)
    axes[i].set_axis_off()
    axes[i].title.set_text(title)
fig.tight_layout()
../../_images/tutorials_linalg_linear_operators_7_0.svg

The elements of the Kronecker product \(A \otimes B\) consist of all possible products of an entry of \(B \in \mathbb{R}^{m_2 \times n_2}\) with an entry of \(A \in \mathbb{R}^{m_1 \times n_1}\). If this was represented as a full matrix instead of a linear operator, one would have to store and multiply with a matrix of dimensions \(m_1m_2 \times n_1n_2\). However, by just defining the matrix-vector product instead, it suffices to store the matrices \(A\) and \(B\).

A useful relative of the Kronecker product is the symmetric Kronecker product. As we will see later it can be used to define a distribution over symmetric matrices.

[6]:
# Symmetric Kronecker Product
C = np.array([[5, 1], [2, 10]])
D = np.array([[1, 2], [3, 4]])
Ws = linops.SymmetricKronecker(C, D)
Vs = linops.SymmetricKronecker(D, C)
[7]:
# Plot
matdict = {
    "$C$": C,
    "$D$": D,
    "$C \otimes_s D$": Ws.todense(),
    "$D \otimes_s C$": Vs.todense(),
}
vmin = np.min([np.min(mat) for mat in list(matdict.values())])
vmax = np.max([np.max(mat) for mat in list(matdict.values())])

fig, axes = plt.subplots(
    nrows=1, ncols=4, figsize=(6, 3.5), gridspec_kw={"width_ratios": [2, 2, 3, 3]}
)
for i, (title, mat) in enumerate(matdict.items()):
    axes[i].imshow(mat, vmin=vmin, vmax=vmax)
    axes[i].set_axis_off()
    axes[i].title.set_text(title)
fig.tight_layout()
../../_images/tutorials_linalg_linear_operators_10_0.svg

Note, that the symmetric Kronecker product generally does not have a symmetric matrix representation. If its arguments are symmetric, then so is their product. We also observe another property of the symmetric Kronecker product, namely that \(C \otimes_s D = D \otimes_s C\). This is not the case for the Kronecker product as we saw above.

Linear Operators and Random Variables

When defining distributions over \(m \times n\) matrices, the covariance has dimension \((mn)^2\), which quickly becomes prohibitively large. This can be somewhat remedied by using linear operators as parameters, since for most computations with matrix-variate random variables the full covariance matrix never needs to be formed explicitly. An example is the matrix-variate normal distribution, where one often assumes Kronecker-structured covariances. One can use the symmetric Kronecker product in this context to define a distribution over symmetric matrices.

[8]:
import scipy.sparse
from probnum import randvars

import probnum.problems.zoo.linalg

# Linear operators as random variable parameters
n = 15
rng = np.random.default_rng(seed=42)
mat = probnum.problems.zoo.linalg.random_spd_matrix(rng=rng, dim=n)

opmean = linops.Matrix(A=mat)
opcov = linops.SymmetricKronecker(2 * linops.Identity(n))

Y = randvars.Normal(mean=opmean, cov=opcov)
print(Y)

# Draw samples
Ysamples = Y.sample(rng=rng, size=3)
<Normal with shape=(15, 15), dtype=float64>
[9]:
# Plot
rvdict = {
    "$\mathbb{E}(\mathsf{Y})$": Y.mean.todense(),
    "$\mathsf{Y}_1$": Ysamples[0],
    "$\mathsf{Y}_2$": Ysamples[1],
    "$\mathsf{Y}_3$": Ysamples[2],
}
vmin = np.min([np.min(mat) for mat in list(rvdict.values())])
vmax = np.max([np.max(mat) for mat in list(rvdict.values())])

fig, axes = plt.subplots(nrows=1, ncols=4, figsize=(10, 5), sharey=True)
for i, (title, rv) in enumerate(rvdict.items()):
    axes[i].imshow(rv, vmin=vmin, vmax=vmax)
    axes[i].set_axis_off()
    axes[i].title.set_text(title)
plt.tight_layout()
../../_images/tutorials_linalg_linear_operators_14_0.svg
[ ]: