Stochastic Processes 12: Gaussian Processes

A Gaussian process (GP) is often described as an infinite-dimensional Gaussian distribution. That definition alone does not explain how a kernel constrains functions, how observations produce a posterior, or why prediction has a closed form. The essential tool is the conditional distribution of a multivariate Gaussian: GP marginalization, regression, and uncertainty updates all follow from the same linear-algebra formulas. Readers who want to review finite-dimensional distributions, covariance, and correlation can begin with Stochastic Processes 2: Joint Distributions and Dependence.

Most Gaussian-process inference is a direct consequence of multivariate Gaussian conditioning.

1. Start with Multivariate Gaussian Conditioning

1.1 Geometry, marginalization, and correlation

The density of a $d$-dimensional random vector $\mathbf{x}\sim N(\boldsymbol\mu,\Sigma)$ is $p(\mathbf{x})=(2\pi)^{-d/2}|\Sigma|^{-1/2}\exp[-\tfrac12(\mathbf{x}-\boldsymbol\mu)^\top\Sigma^{-1}(\mathbf{x}-\boldsymbol\mu)]$. The quadratic form $(\mathbf{x}-\boldsymbol\mu)^\top\Sigma^{-1}(\mathbf{x}-\boldsymbol\mu)=c$ describes an ellipsoid. In the eigendecomposition $\Sigma=U\Lambda U^\top$, the columns of $U$ give its principal directions and $\sqrt{\lambda_i}$ gives the corresponding semiaxis lengths.

All three marginal distributions in the figure are $N(0,1)$; only the correlation coefficient changes. Stronger correlation squeezes the probability cloud into a narrower band. Once $x_1$ is known, the plausible range of $x_2$ also narrows. A GP uses precisely this covariance between input locations to propagate information from observations.

If a subvector $\mathbf{x}_A$ is selected from $\mathbf{x}\sim N(\boldsymbol\mu,\Sigma)$, its marginal distribution remains Gaussian: $\mathbf{x}_A\sim N(\boldsymbol\mu_A,\Sigma_{AA})$. Selecting components is equivalent to multiplication by a selection matrix $S$, so $\mathbf{x}_A=S\mathbf{x}\sim N(S\boldsymbol\mu,S\Sigma S^\top)$. Marginalization merely extracts the relevant mean entries and covariance submatrix; no new integration is required.

1.2 Deriving the conditional distribution

Partition a jointly Gaussian vector into two blocks:

$$\begin{pmatrix}\mathbf{x}_1\\ \mathbf{x}_2\end{pmatrix} \sim N\!\left( \begin{pmatrix}\boldsymbol\mu_1\\ \boldsymbol\mu_2\end{pmatrix}, \begin{pmatrix}\Sigma_{11}&\Sigma_{12}\\ \Sigma_{21}&\Sigma_{22}\end{pmatrix} \right).$$

Set $A=\Sigma_{12}\Sigma_{22}^{-1}$ and construct the residual $\mathbf z=\mathbf x_1-A\mathbf x_2$. Its covariance with the second block is $\operatorname{Cov}(\mathbf z,\mathbf x_2)=\Sigma_{12}-A\Sigma_{22}=0$. The pair $(\mathbf z,\mathbf x_2)$ remains jointly Gaussian, and uncorrelated jointly Gaussian variables are independent. Consequently, $p(\mathbf z\mid\mathbf x_2)=p(\mathbf z)$.

The residual has mean $E[\mathbf z]=\boldsymbol\mu_1-A\boldsymbol\mu_2$ and variance $\operatorname{Var}(\mathbf z)=\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}$. Given $\mathbf x_2$, the relation $\mathbf x_1=\mathbf z+A\mathbf x_2$ yields

$$\mathbf{x}_1\mid\mathbf{x}_2\sim N\!\left( \boldsymbol\mu_1+\Sigma_{12}\Sigma_{22}^{-1}(\mathbf{x}_2-\boldsymbol\mu_2), \Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21} \right).$$

The derivation needs only two facts: linear transformations preserve Gaussianity, and uncorrelated components of a jointly Gaussian vector are independent.

1.3 Reading predictive behavior from the formula

The conditional distribution remains Gaussian, so a GP with fixed hyperparameters and Gaussian observation noise has a closed-form posterior. The correction to the conditional mean is proportional to the cross-covariance $\Sigma_{12}$; if $\Sigma_{12}=0$, observing $\mathbf x_2$ does not change the mean of $\mathbf x_1$. The conditional covariance is the Schur complement $\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}$ and cannot exceed the prior covariance because the subtracted matrix can be written as $A\Sigma_{22}A^\top\succeq0$.

With fixed hyperparameters, conditional covariance depends on observation locations rather than observed values. Active learning can therefore evaluate how much a candidate measurement location would reduce uncertainty before collecting its value. If the lengthscale, signal variance, or other hyperparameters are learned from data, observed values affect posterior covariance indirectly through those hyperparameters.

Take $\Sigma=\begin{pmatrix}1&0.85\\0.85&1\end{pmatrix}$ and observe $x_1=1.6$. The conditional mean changes from 0 to $\rho x_1=1.36$, while the conditional variance falls from 1 to $1-\rho^2\approx0.28$.

A normalized vertical slice through the joint density at $x_1=1.6$ is still bell-shaped, but its center moves to 1.36 and its width contracts substantially.

2. From Finite Gaussian Vectors to Distributions over Functions

2.1 The intuition from denser index sets

Plot a Gaussian vector with component indices on the horizontal axis and component values on the vertical axis. Choose $d$ locations $t_1<\cdots<t_d$, form a covariance matrix with $\Sigma_{ij}=\exp[-(t_i-t_j)^2/(2\ell^2)]$, draw a sample from $N(\mathbf0,\Sigma)$, and connect adjacent components.

For $d=2$, a sample is only a line segment. For $d=200$, the polyline already resembles a smooth function. Every construction step still samples from an ordinary finite-dimensional Gaussian. Increasing the dimension suggests what a distribution over functions means, but a more careful construction must establish that the limit exists and explain where the kernel comes from.

2.2 Deriving a kernel from Bayesian linear regression

Let $\boldsymbol\phi(x)=(\phi_1(x),\ldots,\phi_m(x))^\top$ be a fixed vector of basis functions. Define $f(x)=\boldsymbol\phi(x)^\top\mathbf w$ and assign a Gaussian prior $\mathbf w\sim N(\mathbf0,\Sigma_p)$ to the weights. For any finite inputs $x_1,\ldots,x_n$, the function-value vector is $\mathbf f=\Phi\mathbf w$ and is therefore Gaussian. Its mean is zero, and its covariance is $\operatorname{Cov}[f(x),f(x')]=\boldsymbol\phi(x)^\top\Sigma_p\boldsymbol\phi(x')$. The basis functions and weight prior induce the kernel; a valid kernel is more than an arbitrary similarity score.

As a concrete example, take Gaussian basis functions $\phi_c(x)=\exp[-(x-c)^2/(2\ell^2)]$ with width $\ell$ and centers $c$. As the center spacing approaches zero and the weight variance is scaled appropriately, the discrete sum converges to

$$\begin{aligned} k(x,x') &=s^2\!\int_{-\infty}^{\infty} \exp\!\left[-\frac{(x-c)^2}{2\ell^2}\right] \exp\!\left[-\frac{(x'-c)^2}{2\ell^2}\right]dc\\ &=s^2\ell\sqrt\pi\, \exp\!\left[-\frac{(x-x')^2}{4\ell^2}\right]. \end{aligned}$$

The second line follows by completing the square with $(x-c)^2+(x'-c)^2=2[c-(x+x')/2]^2+(x-x')^2/2$. The result is a squared exponential (SE) kernel with lengthscale $\sqrt2\ell$ and variance $s^2\ell\sqrt\pi$.

With few basis functions, the implied kernel is uneven and the function is confined to an $m$-dimensional space. The numerical kernel approaches the analytic limit as $m$ grows. If $m<n$, then $K=\Phi\Sigma_p\Phi^\top$ has rank at most $m$ and must be singular. A nonparametric model does not lack parameters; instead, the function can use infinitely many parameters, and inference has integrated the weights out.

2.3 Definition and consistency

A stochastic process $\{X_t\}_{t\in T}$ is a Gaussian process if every finite vector $(X_{t_1},\ldots,X_{t_n})$ has a multivariate Gaussian distribution. The definition uses finite-dimensional distributions because any practical prediction also queries only finitely many input locations.

The distributions defined on different finite point sets must be consistent: marginalizing a distribution on a larger set must reproduce the distribution defined directly on its subset. Gaussian marginalization extracts the corresponding blocks of the mean vector and covariance matrix, so finite-dimensional distributions generated by a shared mean function and kernel are automatically consistent. The Kolmogorov extension theorem then guarantees a stochastic process with those finite-dimensional distributions.

Writing $m(t)=E[X_t]$ and $k(s,t)=\operatorname{Cov}(X_s,X_t)$ gives $X\sim\mathcal{GP}(m,k)$. For a general stochastic process, the first two moments do not determine the complete distribution. Every finite-dimensional distribution of a GP is Gaussian, so its mean and kernel determine its full probability law. In practice, data are often centered and $m\equiv0$ is used; the kernel and observations then generate the posterior mean function.

3. What Prior Assumptions Does a Kernel Express?

3.1 Positive semidefiniteness is mandatory

A kernel $k$ must make $K_{ij}=k(t_i,t_j)$ positive semidefinite for every finite point set, because the variance $\mathbf a^\top K\mathbf a$ of any linear combination cannot be negative. A function that decreases with distance is not necessarily a valid kernel. For example, use $k_p(s,t)=\exp(-|s-t|^p)$ on 12 equally spaced points in $[0,3]$:

$p$	Smallest eigenvalue of $K$
1	$+0.1378$
2	approximately $0$
3	$-0.2534$
4	$-0.4876$

The function $\exp(-|r|^p)$ is a valid kernel only for $0<p\le2$. Values $p=3$ and $p=4$ produce negative eigenvalues, implying a linear combination with a negative “variance”; Cholesky factorization also fails. The $p=2$ case is the SE kernel. Its smallest eigenvalues can approach machine precision, so a mathematically valid kernel can still be numerically ill-conditioned.

3.2 Quantifying smoothness through increment variance

For a zero-mean stationary GP, increments satisfy $E|f(t+h)-f(t)|^2=2[k(0)-k(h)]$. The expansion of the kernel near $h=0$ directly determines mean-square smoothness:

• The SE kernel has $k(h)=1-h^2/(2\ell^2)+O(h^4)$, so its first-increment variance is $O(h^2)$.
• The Matérn 3/2 kernel has $k(h)=1-3h^2/(2\ell^2)+O(|h|^3)$, also giving $O(h^2)$ first-increment variance.
• The Matérn 1/2 kernel has $k(h)=1-|h|/\ell+O(h^2)$, giving $O(|h|)$ first-increment variance.
• Brownian motion has $E|B_{t+h}-B_t|^2=|h|$; its paths are continuous but almost surely nowhere differentiable.

The SE and Matérn 3/2 kernels both have quadratic first-increment variance. Their higher-order smoothness can be distinguished with the second difference $E|f(t+h)-2f(t)+f(t-h)|^2=6k(0)-8k(h)+2k(2h)$.

Kernel	First-difference slope	Second-difference slope	Mean-square derivatives
SE	2	4	$\infty$
Matérn 3/2	2	3	1
Matérn 1/2	1	1	0
Brownian	1	1	0

Matérn-$\nu$ samples are generally mean-square differentiable $\lceil\nu\rceil-1$ times. The SE kernel assumes infinite differentiability, a very strong constraint. Matérn 3/2 or 5/2 often better represents data with rough changes or only finite smoothness.

3.3 Lengthscales, common kernels, and composition

The lengthscale $\ell$ controls the distance over which correlation decays. A small $\ell$ correlates only nearby inputs and produces rapidly varying sample functions; a large $\ell$ preserves correlation across longer distances and produces smoother variation. The signal variance $\sigma_f^2$ controls vertical amplitude. An SE sample on an interval of length $L$ has roughly $L/(2\pi\ell)$ zero-crossing-scale oscillations, so a useful lengthscale should match a physically meaningful variation scale.

| Name | Expression, $r=|s-t|$ | Main property | | --- | --- | --- | | Squared Exponential | $\sigma_f^2\exp[-r^2/(2\ell^2)]$ | Infinitely differentiable and extremely smooth | | Matérn 3/2 | $\sigma_f^2(1+\sqrt3r/\ell)e^{-\sqrt3r/\ell}$ | Once mean-square differentiable | | Matérn 1/2 (OU) | $\sigma_f^2e^{-r/\ell}$ | Continuous, nondifferentiable, and Markov | | Periodic | $\sigma_f^2\exp[-2\sin^2(\pi r/p)/\ell^2]$ | Exactly periodic with period $p$ | | Brownian | $\min(s,t)$ | Nonstationary, with variance growing over time |

A wider bright band in a covariance heatmap means stronger long-range correlation and smoother sample functions. A periodic kernel produces parallel bands at integer multiples of the period. The Brownian kernel depends on the earlier time rather than only the time difference, so its heatmap is angular and nonstationary.

Sums and products of positive-semidefinite kernels remain positive semidefinite. A sum $k_1+k_2$ superposes two independent structures, a product $k_1k_2$ enforces both structures, and $k(g(x),g(x'))$ first transforms the inputs. Adding trend, periodic, and short-term disturbance kernels yields a model with all three effects. The product $k_{\mathrm{per}}k_{\mathrm{SE}}$ lets seasonal patterns drift gradually and is often more realistic than exact periodic repetition.

The kernel contains the model’s principal prior assumptions about function space. “Nonparametric” does not mean “assumption-free”: choosing an SE, Matérn, or periodic kernel directly determines smoothness, extrapolation, and correlation structure.

4. GP Regression, Hyperparameters, and Stable Computation

4.1 Deriving the regression equations

Let the training inputs be $X=\{x_i\}_{i=1}^n$, with observations $y_i=f(x_i)+\varepsilon_i$, independent noise $\varepsilon_i\sim N(0,\sigma_n^2)$, and prior $f\sim\mathcal{GP}(0,k)$. Denote latent function values at test inputs $X_*$ by $\mathbf f_*$. The training outputs and test values have joint distribution

$$\begin{pmatrix}\mathbf y\\\mathbf f_*\end{pmatrix} \sim N\!\left(\mathbf0, \begin{pmatrix} K+\sigma_n^2I&K_*\\ K_*^\top&K_{**} \end{pmatrix}\right),$$

where $K=K(X,X)$, $K_*=K(X,X_*)$, and $K_{**}=K(X_*,X_*)$. Substitution into the conditional Gaussian formula from Section 1 gives

$$\begin{aligned} E[\mathbf f_*\mid\mathbf y]&=K_*^\top(K+\sigma_n^2I)^{-1}\mathbf y,\\ \operatorname{Cov}(\mathbf f_*\mid\mathbf y)&=K_{**}-K_*^\top(K+\sigma_n^2I)^{-1}K_*. \end{aligned}$$

Once the kernel and noise parameters are fixed, inference requires solving linear systems rather than gradient descent or iterative sampling. The posterior mean depends on $\mathbf y$, while posterior covariance depends only on input locations and hyperparameters.

4.2 A prediction that can be calculated by hand

Take $k(x,x')=\exp[-(x-x')^2/2]$ and $\sigma_n^2=0.1$. Observe $(x_1,y_1)=(-1,0.5)$ and $(x_2,y_2)=(1,-0.3)$, then predict at $x_*=0$. Since $k(-1,1)=e^{-2}=0.1353$ and $k(-1,0)=k(1,0)=e^{-0.5}=0.6065$, the required matrices are

$$K+\sigma_n^2I= \begin{pmatrix}1.1&0.1353\\0.1353&1.1\end{pmatrix},\qquad \mathbf k_*= \begin{pmatrix}0.6065\\0.6065\end{pmatrix}.$$

Solving $\boldsymbol\alpha=(K+\sigma_n^2I)^{-1}\mathbf y$ gives $\boldsymbol\alpha=(0.4956,-0.3337)^\top$. The posterior mean is $\bar f_*=\mathbf k_*^\top\boldsymbol\alpha=0.0982$, and the posterior variance is $1-\mathbf k_*^\top(K+\sigma_n^2I)^{-1}\mathbf k_*=0.4044$, corresponding to a standard deviation of about 0.636.

The predictive mean is close to the arithmetic mean 0.1 but is not a simple average. The kernel matrix accounts for the correlation of each observation with the test point and for correlation between observations. Because the test location was not directly observed, posterior variance falls from its prior value 1 to 0.4044 rather than to zero. The variance calculation does not use $\mathbf y$, exactly as the Gaussian conditioning formula predicts.

4.3 Two views of the posterior mean

Define $\boldsymbol\alpha=(K+\sigma_n^2I)^{-1}\mathbf y$. The posterior mean becomes $\bar f(x_*)=\sum_{i=1}^n\alpha_i k(x_*,x_i)$, a linear combination of kernel functions centered on the training inputs. A zero prior mean therefore does not force the posterior mean to remain zero.

The same equation can be written as $\bar f(x_*)=\mathbf w(x_*)^\top\mathbf y$, where $\mathbf w(x_*)^\top=\mathbf k_*^\top(K+\sigma_n^2I)^{-1}$. GP regression is consequently a linear smoother: predictions are weighted sums of observations, although weights may be negative and need not be restricted to nearest neighbors. The posterior mean is also identical to the kernel-ridge-regression solution, while the GP probability model additionally supplies a conditional covariance.

With no observations, the posterior equals the prior. One observation shifts the mean and lowers variance nearby. As observations cover the input interval, the posterior band contracts inside the observed region but returns toward the prior width when extrapolating. The noise variance $\sigma_n^2$ lets the posterior mean smooth the data; the model approaches exact interpolation as $\sigma_n^2\to0$.

4.4 Cholesky factorization and ill-conditioning

A GP implementation should never form the matrix inverse explicitly. Factor $A=K+\sigma_n^2I=LL^\top$ with Cholesky decomposition, then use triangular solves for the mean, variance, and log determinant:

L = cholesky(K + σn² I)
α = solve(Lᵀ, solve(L, y))
mean = k*ᵀ α
v = solve(L, k*)
variance = k** - vᵀv
log_determinant = 2 sum(log(diag(L)))

Cholesky factorization is faster and more stable than explicit inversion. The expression $k_{**}-v^\top v$ is also less likely to produce a negative variance from floating-point error. Exact GP computation still costs $O(n^3)$ time and $O(n^2)$ storage.

The eigenvalues of an SE kernel decay rapidly. On the same 12 points, increasing the lengthscale through 0.5, 1, 2, and 4 gives condition numbers of approximately $2.0\times10^5$, $3.9\times10^{11}$, $5.2\times10^{17}$, and $6.7\times10^{17}$. Double precision has a dynamic range of only about $10^{16}$, so the final two matrices are numerically singular. A standard remedy adds diagonal jitter, $K+\epsilon I$, often with $\epsilon\approx10^{-6}\sigma_f^2$. The observation-noise term $\sigma_n^2I$ also improves conditioning, which is why noiseless GP interpolation can be more difficult than noisy regression.

4.5 Learning hyperparameters with marginal likelihood

The lengthscale $\ell$, signal standard deviation $\sigma_f$, and noise standard deviation $\sigma_n$ can change predictions substantially.

A common approach maximizes the marginal likelihood after integrating out the latent function $f$. With $A=K_{\boldsymbol\theta}+\sigma_n^2I$, the log marginal likelihood is

$$\log p(\mathbf y\mid X,\boldsymbol\theta) =-\tfrac12\mathbf y^\top A^{-1}\mathbf y -\tfrac12\log|A| -\tfrac n2\log(2\pi).$$

The first term measures data fit, the second penalizes model capacity, and the third is constant. A larger $\ell$ makes functions more rigid and reduces the complexity penalty, but the fit term deteriorates rapidly once the model becomes too rigid to explain observed variation.

The total objective in the figure reaches its maximum near $\ell\approx1.40$. As $\ell\to0$, $K\to\sigma_f^2I$ and the fit term approaches the constant $-\|\mathbf y\|^2/[2(\sigma_f^2+\sigma_n^2)]$. A short-lengthscale model does not fit perfectly; it treats most structure as independent variation. Marginal likelihood is generally nonconvex and may have competing optima such as “long lengthscale with low noise” and “short lengthscale with high noise.” Multiple starting points are advisable, and small data sets require prior knowledge or cross-validation to check the optimized result.

5. Brownian Bridges, Stationarity, and Spectra

5.1 Constructing a Brownian bridge by conditioning

Brownian motion $\{B_t\}_{t\ge0}$ is itself a GP because every finite-dimensional distribution is jointly Gaussian. If $s<t$, independent increments give $\operatorname{Cov}(B_s,B_t)=\operatorname{Cov}[B_s,B_s+(B_t-B_s)]=\operatorname{Var}(B_s)=s$, so its kernel is $k(s,t)=\min(s,t)$.

Conditioning the entire process on $B_1=0$ produces a Brownian bridge. Take $\mathbf x_1=(B_s,B_t)^\top$ and $\mathbf x_2=B_1$. The cross-covariance is $\Sigma_{12}=(s,t)^\top$ and $\Sigma_{22}=1$. The conditional mean remains zero, while $\operatorname{Cov}(B_s,B_t\mid B_1=0)=\min(s,t)-st$. A Brownian bridge is therefore $\mathcal{GP}(0,\min(s,t)-st)$, with pointwise variance $t(1-t)$ attaining its maximum $1/4$ at $t=1/2$.

Brownian-bridge samples are fixed at zero at both endpoints and have their greatest uncertainty near the midpoint. The example turns a relationship between classical stochastic processes into covariance-matrix conditioning and also demonstrates that a GP need not be stationary or smooth.

5.2 Stationary kernels and Bochner’s theorem

With a constant mean, a kernel $k(s,t)$ that depends only on the lag $s-t$ makes the corresponding GP weakly stationary. If the kernel depends only on the distance $|s-t|$, it is also isotropic. Gaussian-process distributions are fully determined by their means and covariances, so a constant mean and translation-invariant covariance make weak stationarity imply strict stationarity. General stochastic processes do not share that equivalence; Stochastic Processes 3: Stationarity develops the distinction.

Bochner’s theorem characterizes stationary kernels in the frequency domain: a continuous function $k(\tau)$ is a stationary positive-semidefinite kernel if and only if it is the Fourier transform of a finite positive measure, $k(\tau)=\int e^{i\omega\tau}\,dS(\omega)$. The measure $S$ is the process’s power spectrum, so choosing a stationary kernel is equivalent to choosing a nonnegative spectrum.

The SE kernel has a Gaussian spectrum whose high frequencies decay faster than any polynomial, yielding infinitely many mean-square derivatives. A Matérn-$\nu$ spectrum decays as $(1+\omega^2)^{-(\nu+1/2)}$ and retains only finitely many spectral moments, limiting sample smoothness. Peaks at selected frequencies express periodic structure; Spectral Mixture kernels use mixtures of such peaks to approximate stationary kernels. For $p>2$, the Fourier transform of $\exp(-|\tau|^p)$ takes negative values, explaining in the frequency domain why the function is no longer a valid covariance kernel.

6. Practical Limits and a Modeling Checklist

Exact GP computation requires $O(n^3)$ time and $O(n^2)$ storage, becoming a bottleneck before sample sizes reach the tens of thousands. Inducing-point approximations such as FITC and SVGP, Kronecker structure, random Fourier features, and KISS-GP can reduce the cost. Every approximation trades among accuracy, memory, and implementation complexity.

Isotropic kernels also degrade in high-dimensional inputs. Pairwise distances concentrate as dimension grows, and an SE kernel may make most off-diagonal entries nearly equal, leaving little local structure for the model to distinguish. Automatic relevance determination assigns a separate lengthscale to each input dimension and can identify some irrelevant dimensions, but it does not replace meaningful dimensionality reduction, feature design, or a structured prior.

Posterior variance measures uncertainty under the selected kernel, noise model, and hyperparameters; it is not the true prediction error. A misspecified kernel can produce a narrow interval far from reality, and additional data may make the wrong model even more confident. Closed-form regression also assumes independent Gaussian observation noise. Heteroscedastic, heavy-tailed, or correlated noise must be modeled explicitly and often requires approximate inference.

Before applying a GP, check four questions: does the kernel express plausible smoothness and periodic structure; is the covariance matrix positive semidefinite and acceptably conditioned; does the observation-noise model match the data-generating mechanism; and can the sample size support exact factorization? Once those conditions are addressed, the core model remains compact: the mean and covariance functions define a joint Gaussian structure, observations update means and uncertainties through conditioning, and stable linear-system solves perform the actual computation.