{'fraction_non_alphanumeric': 0.0417828736700754, 'fraction_numerical': 0.011104224770168371, 'mean_word_length': 4.653430656934306, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Algorithmic Linearly Constrained Gaussian Processes

We algorithmically construct multi-output Gaussian process priors which satisfy linear differential equations. Our approach attempts to parametrize all solutions of the equations using Gr\"obner bases. If successful, a push forward Gaussian process along the paramerization is the desired prior. We consider several examples from physics, geomathematics and control, among them the full inhomogeneous system of Maxwell's equations. By bringing together stochastic learning and computer algebra in a novel way, we combine noisy observations with precise algebraic computations.


Introduction
In recent years, Gaussian process regression has become a prime regression technique [37]. Roughly, a Gaussian process can be viewed as a suitable 1 probability distribution on a set of functions, which we can condition on observations using Bayes' rule. The resulting mean function is used for regression. The strength of Gaussian process regression lies in avoiding overfitting while still finding functions complex enough to describe any behavior present in given observations, even in noisy or unstructured data. Gaussian processes are usually applied when observations are rare or expensive to produce. Applications range, among many others, from robotics [9], biology [19], global optimization [33], astrophysics [13] to engineering [47].
Incorporating justified assumptions into the prior helps these applications: the full information content of the scarce observations can be utilized to create a more precise regression model. Examples of such assumptions are smooth or rough behavior, trends, homogeneous or heterogeneous noise, local or global behavior, and periodicity (cf. §4 in [37], [11]). Such assumptions are usually incorporated in the covariance structure of the Gaussian process.
Even certain physical laws, given by certain linear differential equations, could be incorporated into the covariance structures of Gaussian process priors. Thereby, despite their random nature, all realizations and the mean function of the posterior strictly adhere to these physical laws 2 . For example, [29,41] constructed covariance structures for divergence-free and curl-free vector fields, which [50,45] used to model electromagnetic phenomena.
A first step towards systematizing this construction was achieved in [24]. In certain cases, a map into the solution set for physical laws could be found by a computation that does not necessarily terminate. Having found such a map, one could assume a Gaussian process prior in its domain and push it forward. This results in a Gaussian process prior for the solutions of the physical laws.
In Section 2, we stress that the map from [24] into the solution set should be a parametrization, i.e., surjective. In Section 3, we combine this with an algorithm which computes this parametrization if it  exists or reports failure if it does not exist. This algorithm is a homological result in algebraic system theory (cf. §7. (25) in [32]). Using Gröbner bases, it is fully algorithmic and works for a wide variety of operator rings; among them the polynomial ring in the partial derivatives, which models linear systems of differential equations with constant coefficients; the (various) Weyl algebras which model linear systems of differential equations with variable coefficients (cf. Example 4.2); and similar rings for difference equations and combined delay differential equations. To demonstrate the use of Gröbner bases, Example 4.3 shows explicit computer algebra code.
Using the results of this paper, one can add information to Gaussian processes 3 not only by (i) conditioning on observations (Bayes' rule), but also by (ii) restricting to solutions of linear operator matrices by constructing a suitable prior.
Since these two constructions are compatible, we can combine strict, global information from equations with noisy, local information from observations. The author views this combination of techniques from homological algebra and machine learning as the main result of this paper, and the construction of covariance functions satisfying physical laws as a proof of concept.
Even though Gaussian processes are a highly precise interpolation tool, they lack in two regards: missing extrapolation capabilities and high computation time, cubically in the amount of observations. These problems have, to a certain degree, been addressed: more powerfull covariance structures [25,21,51,53,6] and several fast approximations to Gaussian process regression [48,18,52,20,10] have been proposed. This paper addresses these two problems from a complementary angle. The linear differential equations allow to extrapolate and reduce the needed amount of observations, which improves computation time.
The promises in this introduction are demontrated in Example 4.1. It constructs a Gaussian process such that all of its realizations satisfy the inhomogeneous Maxwell equations of electromagnetism. Conditioning this Gaussian process on a single observation of electric current yields, as expected, a magnetic field circling around this electric current. This shows how to combine data (the electric current) with differential equations for a global model, which extrapolates away from the data.

Differential Equations and Gaussian Processes
This section is mostly expository and summarizes Gaussian processes and how differential operators act on them. Subsection 2.1 summarizes Gaussian process regression. We then introduce differential (Subsection 2.2) and other operators (Subsection 2.3).

Gaussian processes
A Gaussian process g = GP(µ, k) is a distribution on the set of functions R d → R such that the function values g(x 1 ), . . . , g(x n ) at x 1 , . . . , x n ∈ R d have a joint Gaussian distribution. It is specified by a mean function µ : and a positive semidefinite covariance function Assume the regression model y i = g(x i ) and condition on n observations Denote by k(x, X) ∈ R × n resp. k(X, X) ∈ R n× n 0 the (covariance) matrices obtained by concatenating the matrices k(x, x j ) resp. the positive semidefinite block partitioned matrix with blocks k(x i , x j ). Write µ(X) ∈ R ×n for the matrix obtained by concatenating the vectors µ(x i ) and y ∈ R 1× n for the row vector obtained by concatenating the rows y i . The posterior GP x → µ(x) + (y − µ(X))k(X, X) −1 k(x, X) T , is again a Gaussian process and its mean function is used as regression model.

Differential equations
Roughly speaking, Gaussian processes are the linear objects among stochastic processes. Hence, we find a rich interplay of Gaussian processes and linear operators.
be the polynomial ring in the partial differential operators. For different or more general operator rings see Subsection 2.3. This ring models linear (partial) differential equations with constant coefficients, as it acts on the vector space F = C ∞ (R d , R) of smooth functions, where ∂ xi acts by partial derivative w.r.t. x i . The set of realizations 4 of a Gaussian process with squared exponential covariance function is dense in F (cf. Thm. 12, Prop. 42 in [43]).
The class of Gaussian processes is closed under matrices B ∈ R × of linear differential operators with constant coefficients. Let g = GP(µ, k) be a Gaussian process with realizations in a space F of vectors with functions in F as entries. Define the Gaussian process B * g as the Gaussian process induced by the pushforward measure under B of the Gaussian measure induced by g. It holds that where B denotes the operation of B on functions with argument The covariance function k for such Gaussian processes B * g as in (1) is often singular. This is to be expected, as B * g is rarely dense in F . For numerical stability, we tacitly assume the model y i = g(x i ) + ε for small Gaussian white noise term ε and adopt k by adding var(ε) to k(x i , x i ) for observations x i . Example 2.1. Let g = GP(0, k(x, x )) be a scalar univariate Gaussian process with differentiable realizations. Then, the Gaussian process of derivatives of functions is given by One can interpret this Gaussian process ∂ ∂x * g as taking derivatives as measurement data and producing a regression model of derivatives.
We say that a Gaussian process is in a function space, if its realizations are contained in said space.
Such solution sets and Gaussian processes are connected in an almost tautological way. Lemma 2.2. Let g = GP(µ, k) be a Gaussian process in F ×1 . Then g is also a Gaussian process in sol F (A) for A ∈ R × if and only if µ ∈ sol F (A) and A * (g − µ) is the constant zero process.
Proof. Assume that g is a Gaussian process in sol F (A). Then, the mean function is a realization, thus µ ∈ sol F (A). Furthermore, forg := (g − µ) = GP(0, k) all realizations are annihilated by A, and hence A * g is the constant zero process.
Conversely, assume that µ ∈ sol F (A) and A * (g − µ) is the constant zero process. This implies 0 = A * (g − µ) = A * g − A * µ = A * g, i.e. all realizations of g become zero after a pushforward by A. In particular, all realizations of g are contained in sol F (A).
This lemma implies another advantage of choosing a zero mean function: it is always a solution of the linear differential equations.
Our goal is to construct Gaussian processes with realizations dense in the solution set sol F (A) of an operator matrix A ∈ R × . The following remark, implicit in [24], is a first step towards an answer.
Parametrizations yield the denseness of the realizations of a Gaussian process B * g in sol F (A).
This proposition is a consequence of partial derivatives being bounded, and hence continuous, when F is equipped with the Fréchet topology generated by the family of seminorms [49]). Now, the continuous surjective map induced by B maps a dense set to a dense set.

Further operator rings
The theory presented for differential equations with constant coefficients also holds for other rings R For ordinary linear differential equations with rational 5 coefficients consider the Weyl algebra R = R(t) ∂ t , with the non-commutative relation ∂ t t = t∂ t + 1 representing the product rule of differentiation. We consider solutions in the set F of smooth functions defined on a co-finite set.
The polynomial ring R = R[σ x1 , . . . , σ x d ] models linear shift equations with constant coefficients when it acts on the set F = R Z d ≥0 of d-dimensional sequences by translation of the arguments.

Computing parametrizations
By the last section, constructing a parametrization B of sol F (A) yields a Gaussian process dense in the solution set sol F (A) of an operator matrix A ∈ R × . Subsection 3.1 gives necessary and sufficient conditions for a parametrization to exist and Subsection 3.2 describes their computation.

Existence of parametrizations
It turns out that we can decide whether a parametrization exists purely algebraically, only using operations over R that do not involve F.
By r-ker(A) we denote the right kernel of A ∈ R × , i.e. r-ker(A) = {m ∈ R ×1 | Am = 0}. By l-ker(A) we denote the left kernel of A, i.e. l-ker(A) = {m ∈ R 1× | mA = 0}. Abusing notation, denote any matrix as left resp. right kernel if its rows resp. columns generate the kernel as an R-module.    The approach assigns a prior to the parametrising functions and pushes this prior forward to a prior of the solution set sol F (A). The paramerization is not canonical, and hence different parametrizations might lead to different priors.

Algorithms
Summarizing Theorem 3.1 and Corollary 3.2 algorithmically, we need to compute right kernels (of A), compute left kernels (of B), and decide whether rows (of A ) are contained in a row module (generated by the rows of A). All these computations are an application of Gröbner basis algorithms.
In the recent decades, Gröbner bases algorithms have become one of the core algorithms of computer algebra, with manifold applications in geometry, system theory, natural sciences, automatic theorem proving, post-quantum cryptography, and many others. Reduced Gröbner bases generalize the reduced echelon form from linear systems to systems of polynomial (and hence linear operator) equations, by bringing them into a standard form 8 . They are computed by Buchberger's algorithm, which is a generalization of the Gaussian and Euclidean algorithm and a special case of the Knuth-Bendix completion algorithm.
Similar to the reduced echelon form, Gröbner bases allow to compute all solutions over R (not F) of the homogeneous system and compute, if it exists, a particular solution over R (not F) for an inhomogeneous system. Solving homogeneous systems is the same as computing its right resp. left kernel. Solving inhomogeneous equations decides whether an element is contained in a module. Alternatively, the uniqueness of reduced Gröbner bases also decides submodule equality.
A formal description of Gröbner bases would exceed the scope of this note. Instead, we refer to the excellent literature [46,12,1,17,14,5]. Gröbner basis algorithms exist for many rings R. They historically emerged from polynomial rings, and have since been generalized to the Weyl algebra, the shift algebra, and, more generally, G-algebras [26,27] and Ore-algebras [39,38]. They are implemented in various computer algebra systems, Singular [8] and Macaulay2 [16] are two well-known examples. Even though the complexity of Gröbner bases is in the vicinity of EXPSPACE completeness (cf. [30,31,3]), the "average interesting example" (e.g. every example in this paper) usually terminates instantaneously. This holds in particular since the Gröbner basis computations only involve the operator equations, but not the data in any way.

Hyperparameters
Many covariance functions 9 incorporate hyperparameters and advanced methods specifically add more hyperparameters to Gaussian processes, see e.g. [44,6,51], for additional flexibility. The approach in this paper is the opposite by restricting the Gaussian process prior to solutions of an operator matrix. Of course, the prior of the parametrizing functions can still contain hyperparameters, which can be determined by maximizing the likelihood. Many important applications contain unknown parameters in the equations. Such parameters can also be estimated by the likelihood.
Consider ordinary differential equations, with constant resp. variable coefficients. The solution set of an operator matrix is a direct sum of parametrizable functions and a finite dimensional set of functions, due to the Smith form resp. Jacobson form. In many cases, in particular the case of constant coefficients, the solution set of the finite dimensional summand can easily be computed. This paper also allows to compute with the parametrizable summand of the solution set and estimate parameters and hyperparameters of both summands together.

Examples
Example 4.1. Maxwell's equations of electromagnetism uses curl and divergence operators as building blocks. It is a well-known result that the solutions of the inhomogeneous Maxwell equations are parametrized by the electric and magnetic potentials. We verify this and use the parametrization to construct a Gaussian process, such that its realizations adhere to Maxwell's equations. In Figure 1, we condition this prior on a single observation of flowing electric current, which leads to the magnetic field circling around the current. This usage of differential equations shows an extrapolation away from the data point in space and into other components.
The inhomogenous Maxwell equations are given by the operator matrix applied to three components of the electric field, three components of the magnetic (pseudo) field, three components of electric current, and one component of electric flux. We have set all constants to 1.
Using Gröbner bases, one computes the right kernel of A and verifies that it parametrizes the set of solutions of the inhomogeneous Maxwell equations.
For the demonstration in Figure 1 we assume squared exponential covariance functions and a zero mean function for four uncorrelated parametrising functions (electric potential and magnetic potentials).  This system, given by the matrix A := ∂ t −t 3 , is parametrizable by For a parametrizing functions with squared exponential covariance functions k(t 1 , t 2 ) = exp(− 1 2 (t 1 − t 2 ) 2 ) and a zero mean function, the covariance function for (x, u) is For a demonstration of how to observe resp. control such a system see Figures 2 resp. 3. x(t) Figure 3: We control the system in Example 4.2 by specifying a desired behavior for the state x(t) and letting the Gaussian process construct a suitable input u(t), which is completely unspecified by us. Starting with x(1) = 1 we give u(t) one time step to control x(t) to zero, e.g., by setting x(t) = 0 for t ∈ { 20 10 , 21 10 , . . . , 5}.
be the polynomial ring in three indeterminates, which we can both interpret as the polynomial ring in the coordinates resp. in the differential operators.
Consider the matrix A = [x 1 x 2 x 3 ] representing the normals of circles centered around the origin resp. the divergence. The right kernel of A is given by the operator representing tangent spaces of circles centered around the origin resp. the curl, and these parametrize the solutions of A. A posterior mean field is demonstrated in Figure 4 when assuming equal covariance functions k for 3 uncorrelated parametrizing functions and the covariance function for the tangent field is We demonstrate how to compute B and A for this example using Macaulay2 [16]. We construct a prior for smooth tangent fields on the sphere without sources and sinks. We work in the third polynomial Weyl algebra R = R[x, y, z] ∂ x , ∂ y , ∂ z . I.e., we are interested in sol A (F) = {v ∈ C ∞ (S 2 , R 3 )|Av = 0} for The right kernel can be checked to yield a parametrization of sol F (A) Assuming a squared exponential covariance functions k for the parametrizing function, a demonstration can be found in Figure 4

Conclusion
The paper constructs multi-output Gaussian process priors, which adhere to linear operator equations.
With these priors few observations yield a precise regression model with strong extrapolation capabilities (cf. Examples 4.1, 4.3, and 4.4). This construction is fully algorithmic and rather general, as it allows linear systems of differential equations with constant or variable coefficients, shift equations, or multiplications with variables. It could be applied to settings from physics (cf. Examples 4.1), geometric settings with potential applications in geomathematics and weather prediction (cf. Examples 4.1, 4.3, and 4.4), or to observe and control systems (cf. Example 4.2). The main restriction is that the solutions of the system of equations must be parametrizable.
The author hopes that the results can be generalized from parametrizable solution sets to the general case using a Monge parametrization (computable via the purity filtration [36,35,2]) and right hand sides [15]. It would also be interesting to apply to parameter estimation (cf. Example 4.2), boundary conditions [15], and to clarify the connection between the algebra, functional analysis, topology, and measure theory used in this paper. Finally, experimental results would be interesting which covariance function for the parametrizing functions is most suitable.