<jupyter_start><jupyter_text>ProbNum QuickstartProbNum implements probabilistic numerical methods in Python. Such methods quantify _uncertainty arising from finite computation_ or from _stochastic input_.Below we explain how to get started with ProbNum and its basic functionality.<jupyter_code># 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.rcParams['font.size'] = 18 
plt.rcParams['text.usetex'] = True
plt.rcParams['text.latex.preamble'] = [r'\usepackage{amsfonts}', 
                                       r'\usepackage{amsmath}', 
                                       r'\usepackage{bm}']<jupyter_output><empty_output><jupyter_text>InstallationYou can install ProbNum using `pip` (or `pip3`).```bashpip install probnum```Alternatively, you can install the package from source.```bashpip install git+https://github.com/probabilistic-numerics/probnum.git```Then in a python session simply import ProbNum.<jupyter_code>import probnum<jupyter_output><empty_output><jupyter_text>Basic ConceptsThe main objects of interest in ProbNum are random variables. `RandomVariable`s have a `Distribution` which models the (numerical) uncertainty on the variable in question.<jupyter_code>import numpy as np
from probnum.prob import RandomVariable, Normal

np.random.seed(1)
x = RandomVariable(distribution=Normal(0, 1))
print(x.sample())<jupyter_output>1.6243453636632417<jupyter_text>`RandomVariable`s behave similarly to NumPy arrays and support basic arithmetic, indexing and slicing.<jupyter_code>y = 2 * x + 1
print(f"Mean: {y.mean()} and covariance: {y.cov()} of y.")<jupyter_output>Mean: 1.0 and covariance: 4.0 of y.<jupyter_text>Probabilistic Numerical MethodsPN methods solve numerical problems (e.g. solution of linear systems, quadrature, differential equations, ...) by treating them as _statistical inference problems_ instead.At a basic level they can serve as drop-in replacements for classic numerical routines.<jupyter_code># Linear System Ax=b
A = np.array([[7.5, 2.0, 1.0],
              [2.0, 2.0, 0.5],
              [1.0, 0.5, 5.5]])
b = np.array([1, 2, -3])

# Solve using NumPy
x = np.linalg.solve(A, b)
print(x)

# Solve using ProbNum
x_rv, _, _, info = probnum.linalg.problinsolve(A, b)
print(x_rv.mean())<jupyter_output>[-0.12366738  1.28358209 -0.63965885]
[-0.12366738  1.28358209 -0.63965885]<jupyter_text>However, probabilistic numerical methods return random variables instead of just numbers. Their distribution models the uncertainty arising from finite computation or stochastic input.<jupyter_code># Solve with limited computational budget
x_rv, _, _, _ = probnum.linalg.problinsolve(A, b, maxiter=2)
# Covariance of solution representing uncertainty
print(f"Covariance matrix: \n{x_rv.cov().todense()}")
# Sample from output random variable
n_samples = 3
x_samples = x_rv.sample(n_samples)
# Plot of true solution, mean and samples
rvdict = {"$x_*$" : x, "$\mathbb{E}(\mathsf{x})$" : x_rv.mean(), 
          "$\mathsf{x}_1$" : x_samples[0], "$\mathsf{x}_2$" : x_samples[1], "$\mathsf{x}_3$" : x_samples[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=2 + n_samples, figsize=(8, 2.5), sharey=True)
for i, (title, rv) in enumerate(rvdict.items()):
    axes[i].imshow(rv[:, np.newaxis], vmin=vmin, vmax=vmax, cmap="bwr")
    axes[i].set_axis_off()
    axes[i].title.set_text(title)
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Here, the probabilistic linear solver has identified one component of the solution already with a high degree of confidence, while there is still some uncertainty about the others left due to early termination.  Encoding Prior KnowledgeIf we have prior knowledge about the problem setting, we can encode this into a PN method by specifying a prior distribution on the input. For this problem we observe that the matrix $A$ is symmetric. Additionally, suppose we are given an approximate inverse of the system matrix.<jupyter_code># Approximate inverse of A
Ainv_approx = np.array([[ 0.2  , -0.18, -0.015],
                        [-0.18 ,  0.7 , -0.03 ],
                        [-0.015, -0.03,  0.20 ]])
print(A @ Ainv_approx)<jupyter_output>[[1.125  0.02   0.0275]
 [0.0325 1.025  0.01  ]
 [0.0275 0.005  1.07  ]]<jupyter_text>For this problem we encode the symmetry of the system matrix and its inverse by using symmetric matrix-variate normal distributions as priors. Realizations of these priors are symmetric matrices. Since there is no stochasticity in the problem definition we choose a highly concentrated symmetric prior over $A$. This codifies our certainty about (multiplications with) $A$. For the prior over the inverse we choose the approximate inverse from above as a mean.<jupyter_code>from probnum.linalg.linops import SymmetricKronecker, Identity

# Prior distribution(s)
A0 = RandomVariable(distribution=Normal(mean=A, 
                                        cov=SymmetricKronecker(10 ** -6 * Identity(A.shape[0]))))
Ainv0 = RandomVariable(distribution=Normal(mean=Ainv_approx, 
                                           cov=SymmetricKronecker(0.1 * Identity(A.shape[0]))))<jupyter_output><empty_output><jupyter_text>We now solve the problem again, but this time we also pass the prior distributions to the solver.<jupyter_code># Solve linear system with limited computational budget and prior knowledge
x_rv, _, _, info = probnum.linalg.problinsolve(A, b, A0=A0, Ainv0=Ainv0, maxiter=2)
print(info)<jupyter_output>{'iter': 2, 'maxiter': 2, 'resid_l2norm': 0.00014783150371253874, 'trace_sol_cov': 0.007808118710123991, 'conv_crit': 'maxiter', 'rel_cond': None}<jupyter_text>Examining the PosteriorWe now take another look at the uncertainty in the solution, given prior knowledge about the inverse.<jupyter_code># Covariance of solution representing uncertainty
print(f"Covariance matrix: \n{x_rv.cov().todense()}")
# Sample from output random variable
n_samples = 3
x_samples = x_rv.sample(n_samples)
# Plot of true solution, mean and samples
rvdict = {"$x_*$" : x, "$\mathbb{E}(\mathsf{x})$" : x_rv.mean(), 
          "$\mathsf{x}_1$" : x_samples[0], "$\mathsf{x}_2$" : x_samples[1], "$\mathsf{x}_3$" : x_samples[2]}

fig, axes = plt.subplots(nrows=1, ncols=2 + n_samples, figsize=(8, 2.5), sharey=True)
for i, (title, rv) in enumerate(rvdict.items()):
    axes[i].imshow(rv[:, np.newaxis], vmin=vmin, vmax=vmax, cmap="bwr")
    axes[i].set_axis_off()
    axes[i].title.set_text(title)
plt.tight_layout()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11962505579526857, 'fraction_numerical': 0.03704805832465407, 'mean_word_length': 4.6821639898562974, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 36, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '70321', 'license': 'MIT', 'n_tokens_mistral': 2364, 'n_tokens_neox': 2143, 'n_words': 639, 'path': 'docs/source/introduction/quickstart.ipynb', 'repo_name': 'coldfix/probnum', 'size': '6721'}