{'fraction_non_alphanumeric': 0.04836167837678114, 'fraction_numerical': 0.011097248670300085, 'mean_word_length': 4.0977405857740585, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the Bayesian Solution of Differential Equations

The interpretation of numerical methods, such as finite difference methods for differential equations, as point estimators allows for formal statistical quantification of the error due to discretisation in the numerical context. Competing statistical paradigms can be considered and Bayesian probabilistic numerical methods (PNMs) are obtained when Bayesian statistical principles are deployed. Bayesian PNM are closed under composition, such that uncertainty due to different sources of discretisation can be jointly modelled and rigorously propagated. However, we argue that no strictly Bayesian PNM for the numerical solution of ordinary differential equations (ODEs) have yet been developed. To address this gap, we work at a foundational level, where a novel Bayesian PNM is proposed as a proof-of-concept. Our proposal is a synthesis of classical Lie group methods, to exploit the underlying structure of the gradient field, and non-parametric regression in a transformed solution space for the ODE. The procedure is presented in detail for first order ODEs and relies on a certain technical condition -- existence of a solvable Lie algebra -- being satisfied. Numerical illustrations are provided.


Introduction
Numerical methods underpin almost all of scientific, engineering and industrial output. In the abstract, a numerical task can be formulated as the approximation of a quantity of interest Q : Y → Q, subject to a finite computational budget. The true underlying state y † ∈ Y is typically high-or infinite-dimensional, so that only limited information 1 0 y(x)dx and A(y) might consist of a finite number of evaluations [y(x 1 ), . . . , y(x n )] of the integrand. In this abstract framework a numerical method (e.g. the trapezoidal rule) corresponds to a map b : A → Q, as depicted in Figure 1a, where b(a) represents an approximation to the quantity of interest based on the information a ∈ A.
The increasing ambition and complexity of contemporary applications is such that the computational budget can be extremely small compared to the precision that is required at the level of the quantity of interest. As such, in many important problems it is not possible to reduce the numerical error to a negligible level. Fields acutely exposed to this challenge include climate forecasting [Wedi, 2014], computational cardiology [Chabiniok et al., 2016] and molecular dynamics [Perilla et al., 2015]. In the presence of non-negligible numerical error, it is unclear how scientific interpretation of the output of computation can proceed.

Probabilistic Numerical Methods
The field of probabilistic numerics dates back to Larkin [1972] and aims to provide formal statistical foundations for numerical methods. Indeed, under the abstract framework just described, numerical methods can be interpreted as point estimators in a statistical context, where a state y ∈ Y can be thought of as a latent variable in a statistical model, and the 'data' consist of information A(y) that does not fully determine the quantity of interest Q(y) but is indirectly related to it. Hennig et al. [2015] provide an accessible introduction and survey of the field.
Let the notation Σ Y denote a σ-algebra on the space Y, assumed fixed, and let P Y denote the set of probability measures on (Y, Σ Y ). A probabilistic numerical method (PNM) is a procedure which takes as input a 'belief' distribution µ ∈ P Y , representing epistemic uncertainty with respect to the true (but unknown) value y † , along with a finite amount of information, A(y † ) ∈ A. The output is a distribution B(µ, A(y † )) ∈ P Q on (Q, Σ Q ), representing epistemic uncertainty with respect to the quantity of interest Q(y † ) after the information A(y † ) has been processed. For example, a PNM for an ordinary differential equation (ODE) takes an initial belief distribution defined on the solution space of the differential equation, together with information arising from a finite number of evaluations of the gradient field, plus the initial condition of the ODE, to produce a distribution over either the solution space of the ODE, or perhaps some derived quantity of interest. In this paper, the measurability of A and Q will be assumed.
Despite computational advances in this emergent field, until recently there had not been an attempt to establish rigorous statistical foundations for PNM. In Cockayne et al. [2017] the authors argued that Bayesian principles can be adopted. In brief, this framework requires that the output of a PNM agrees with a rational Bayesian agent: Let Q # : P Y → P Q denote the push-forward map associated to Q. i.e. Q # (µ)(S) = µ(Q −1 (S)) for all S ∈ Σ Q . Let {µ a } a∈A ⊂ P Y denote the disintegration, assumed to exist 1 , of µ ∈ P Y along the map A. Figure 1: Diagrams for a numerical method. (a) The traditional viewpoint of a numerical method is equivalent to a map b from a finite-dimensional information space A to the space of the quantity of interest Q. (b) The probabilistic viewpoint treats approximation of Q(y † ) in a statistical context, described by a map B(µ, ·) from A to the space of probability distributions on Q. The probabilistic numerical method (A, B) is Bayesian if and only if (b) is a commutative diagram.
This definition is intuitive; the output of the PNM should coincide with the marginal distribution for Q(y † ) according to the disintegration element µ a ∈ P Y , based on the information a ∈ A that was provided. The definition is equivalent to the statement that Figure 1b is a commutative diagram. In Cockayne et al. [2017] the map A was termed an information operator and the map B was termed a belief update operator ; we adhere to these definitions in our work.
The Bayesian approach to PNM confers several important benefits: • The input µ and output B(µ, a) belief distributions can be interpreted, respectively, as a prior and (marginal) posterior. As such, they automatically inherit the stronger formal semantics and philosophical foundations that underpin the Bayesian framework and, in this sense, are well-understood [see e.g. Gelman and Shalizi, 2013].
• The definition of Bayesian PNM is operational. Thus, if we are presented with a prior µ and information a then there is a unique Bayesian PNM and it is constructively defined.
• The class of Bayesian PNM is closed under composition, such that uncertainty due to different sources of discretisation can be jointly modelled and rigorously propagated. This point will not be discussed further in this work, but we refer the interested reader to Section 5 of Cockayne et al. [2017].
Nevertheless, the strict definition of Bayesian PNM limits scope to design convenient computational algorithms and indeed several proposed PNM are not Bayesian [see of the 'conditional distribution of y given A(y) = a' in the context of this work. The disintegration theorem, Thm. 1 of Chang and Pollard [1997], guarantees existence and uniqueness of the disintegration up to a A # µ-null set under the weak requirements that Y is a metric space, Σ Y is the Borel σ-algebra, µ is Radon, Σ A is countable generated and Σ A contains all singletons {a} for a ∈ A. Table 1 in Cockayne et al., 2017]. The challenge is two-fold; for a Bayesian PNM, the elicitation of an appropriate prior distribution µ and the exact computation of its disintegration {µ a } a∈A must both be addressed. In the next section we argue that -perhaps as a consequence of these constraints -a strictly Bayesian PNM for the numerical solution of an ODE does not yet exist.

Early Work on PNM for ODEs
Consider a generic univariate first-order initial value problem The notational convention used in this paper is that y denotes a generic function, whereas the italicised y denotes a generic value taken by that function. Throughout this paper all ODEs that we consider will be assumed to be well-defined and admit a unique solution y † ∈ Y where Y is some pre-specified set. In this paper the quantity of interest Q(y † ) will be the solution curve y † itself. The first probabilistic perspective on ODEs, of which we are aware, was Skilling [1992]. Originally described as 'Bayesian' by the author, we will argue that, at least in the strict sense of Definition 1, it is not a Bayesian PNM. In general, a PNM attempts to reason about y † based on evaluations f (x i , y i ) of the gradient field for certain input The selection of the input pairs (x i , y i ) on which f is evaluated is not constrained and several possibilities, of increasing complexity, were discussed in Skilling [1992]. To fix ideas, the simplest such approach is to proceed iteratively as follows: (0.1) The first pair (x 0 , y 0 ) is fully determined by the initial condition of the ODE. 3) The belief µ is updated according to a 0 , leading to a new belief µ 0 which is just the disintegration element µ a 0 .
(1) A discrete time step n > 0, is performed and a particular point estimate y 1 = y(x 1 )dµ 0 (y) for the unknown true value y † (x 1 ) is obtained. This specifies the second pair (x 1 , y 1 ). (Here y 1 is the predictive mean for y(x 1 ) based on the current belief µ 0 .) The process continues similarly, such that at time step i−1 we have a belief distribution µ i−1 = B(µ, a i−1 ) ∈ P Y , where the general belief update operator B is yet to be defined, and the following step is performed: The final output is a probability distribution µ n = B(µ, a n ) ∈ P Y .
The method just described is not actually a PNM in the concrete sense that we have defined. Indeed, the final output µ n is a deterministic function of the values a n of the gradient field that were obtained. However, the values of the gradient field f (x, y) outside any open neighbourhood of the true solution curve C = {(x, y) : y = y † (x), x ∈ [x 0 , x T ]} do not determine the solution of the ODE and, conversely, the solution of the ODE provides no information about the values of the gradient field outside any open neighbourhood of the true solution curve C. Thus it is not possible, in general, to write down an information operator A : Y → A that reproduces the information a n when applied to the solution curve y † (·) of the ODE.
The approach taken in Skilling [1992], cast in our abstract framework, was therefore to posit an approximate information operator Of course,Â(y † ) = a n in general. To acknowledge the approximation error, Skilling [1992] proposed a particular belief update operator B that attempts to model the information with a Gaussian potential: This Gaussian potential was referred to in Skilling [1992] as a 'likelihood' and, together with µ 0 = µ a 0 , the output µ n is completely specified. Here σ is a fixed positive constant, however in principle a non-diagonal covariance matrix can also be considered. The negative consequences of basing inferences on an approximate information operatorÂ are potentially twofold. First, in the special case where the gradient field f does not depend the second argument, the quantities dy dx (x i ) and f (x i , y i ) are identical. From this perspective, µ n represents inference under a mis-specified likelihood, since information is treated as erroneous when it is in fact exact. Second, recall that values of the gradient field that are not contained on the true solution curve C of the ODE are ancillary with respect to the quantity of interest y † . It follows that dependence of the output µ n on these ancillary values violates the conditionality principle (CP) [Cox and Hinkley, 1974, Sec. 2.2]. Since the CP is an implication of the likelihood principle, violation of the CP implies violation of the Bayesian framework. This confirms, through a different argument, that the approach of Skilling [1992] cannot be Bayesian in the sense of Definition 1.

Recent Work on PNM for ODEs
Inspired by Skilling [1992], several variations on the method have been proposed: • The approach of Schober et al. [2014] considered Eq. (2) in the σ ↓ 0 limit. The authors proved that if the input pair (x 1 , y 1 ) is taken as y 1 = y(x 1 )dµ 0 (y), and a particular belief µ was used, then the smoothing estimateŷ 1 = y(x 1 )dµ 1 (y), i.e. the posterior mean for y(x 1 ) based on information a 1 , coincides with the deterministic approximation to y † (x 1 ) that would be provided by a k-th order Runge-Kutta method. Similar connections to multistep methods of Nordsieck and Adams form were identified, respectively, in Schober et al. [2016] and Teymur et al. [2016,2018].
• The work of Kersting and Hennig [2016] proposed an adaptive choice of covariance matrix for use in Eq.
(2), in order to encourage uncertainty estimates to be better calibrated.
Exact Bayesian PNM ? Figure 2: Schematic of our proposed approach. An nth order ODE that admits a solvable Lie algebra can be transformed into n integrals, to which exact Bayesian probabilistic numerical methods can be applied. The posterior measure on the transformed space is then pushed back through the inverse transformation onto the original domain of interest.
• The original work of Chkrebtii et al. [2016] is somewhat related to Kersting and Hennig [2016], however instead of using the mean of the current posterior as input to the gradient field, the input pair (x i , y i ) was selected by sampling y i from the marginal distribution for y(x i ) implied by µ i−1 .
• The approaches proposed in Conrad et al. [2017], Abdulle and Garegnani [2018] are not motivated in the Bayesian framework, but instead seek to introduce a stochastic perturbation into a classical numerical method.
The above methods provide elegant and practical uncertainty quantification for ODEs, but do not fulfil the requirements of a Bayesian PNM.

Our Contribution
This paper presents a proof-of-concept PNM for the numerical solution of ODEs that is strictly Bayesian in the sense of Definition 1. For an nth order ODE, we require the ODE to admit an n-dimensional solvable Lie algebra. Then, our proposal is to perform exact Bayesian inference on a particular Lie-transformed ODE, whose gradient field is a function of the independent state variable only, in effect reducing the ODE to an integral (see Fig. 2). Only the case of first order ODE is presented in this paper, due to page limitations, but the method itself is general 2 . In addition to the benefits conferred in the Bayesian framework, detailed in Section 1.1 and in Cockayne et al. [2017], the method being proposed can be computationally realised, but at a substantially increased run-time compared to existing, non-Bayesian approaches. As such, we consider this work to be a proof-of-concept rather than an applicable Bayesian PNM.

Methods
In this section our novel Bayesian PNM is presented for the case of First Order ODEs. This allows some of the more technical details associated to the general case to be omitted, due to the fact that any one-dimensional Lie algebra is trivially solvable. Due to page restrictions, an overview of basic Lie group methods is reserved for the Appendix, and we also refer the reader to chapters 1-3 of Bluman and Anco [2002] for complete detail. The main result from Lie group methods that will allow us to construct an exact Bayesian PNM is as follows: Theorem 1 (Reduction of a First Order ODE to an Integral). If a first order ODE admits a one parameter Lie group of transformations, then there exists canonical coordinates r(x, y), s(x, y) such that for some explicit function G(r).
Proof. See Appendix.
Note that the transformed ODE in Eq. (4) is nothing more than an integral, for which exact Bayesian PNM can be used [e.g. Briol et al., 2015]. In general, if the Lie algebra exhibited by the symmetries of an nth order ODE contains a n-dimensional solvable subalgebra, then that ODE can be reduced to n nested integrals and the method we are about to present can be employed. At a high level, as indicated in Fig. 2, our proposed PNM performs exact Bayesian inference for the solution s ∈ S of Eq. (4) and then transforms the resultant posterior in P S back into the original coordinate system to obtain an element in P Y . From Def. 1, our Bayesian PNM is uniquely determined by the combination of the prior ν ∈ P S and the information operator which corresponds indirectly to n + 1 evaluations of the original gradient field f at certain input pairs (x i , y i ). The selection of the inputs r i is in principle unrestricted. The transformation of a first order ODE is now briefly illustrated: Example 1. Consider the first order ODE This ODE has the canonical coordinates s = log y, r = y x . The transformed ODE is then Lie transform (Lie transform) −1 Figure 3: Illustration of the implicit prior principle: A prior elicited for the function s(r) in the transformed coordinate system (r, s) must be supported on functions s(r) that correspond to well-defined functions y(x) in the original coordinate system (x, y). Thus the situation depicted would not be allowed.
Thus an evaluation G(r) corresponds to an evaluation of f (x, y) at an input (x, y) such that r = y x . The approach just described cannot proceed unless the prior ν ∈ P S corresponds, implicitly, to a well-defined distribution µ ∈ P Y in the original coordinate system Y. This precludes standard (e.g. Gaussian process) priors in general, as such priors assign mass to functions in (r, s)-space that do not correspond to well-defined functions in (x, y)-space (see Fig. 3).
For such an implicit prior to be well-defined we need to understand when a function in (r, s) space maps to a well-defined function in the original (x, y) domain of interest.
Principle (Implicit Prior). A distribution ν ∈ P S on the transformed solution space S corresponds to a well-defined implicit prior µ ∈ P Y provided that x(r, s(r)) is strictly monotone as a function of r, or equivalently (without loss of generality) dx/dr > 0 for all r.
Thus dx/dr > 0 ⇔ s (r) > 1/r and the implicit prior principle requires that we respect the constraint The prior ν ∈ P S must therefore be supported on differentiable functions s on r ∈ (0, ∞) and satisfying Eq. (2.1).

Experimental Results
In this section the proposed Bayesian PNM is experimentally tested. Again the case of a first order ODE is considered. Scope is limited to verifying the correctness of the procedure, as well as indicating how implicit prior distributions can be constructed.
To limit scope, we consider our running example from Sec. 2: Note the associated canonical coordinates for this class of ODE have already been derived in Ex. 1. In constructing a prior µ ∈ P Y , we must respect the implicit prior principle in Sec. 2. Indeed, recall from Ex. 1 that the ODE in Eq. (2.2) can be transformed into an ODE of the form ds dr = G(r), r ∈ (0, ∞), s(y 0 ) = log(y 0 ).
Then our approach constructs a distribution ν ∈ P S where, from Ex. 2, S is the set of differentiable functions s defined on r ∈ (0, ∞) and satisfying Eq. (2.1). Note that the constraints in Eq. (2.1) preclude the direct use of standard prior models, such as Gaussian processes. However, it is nevertheless possible to design priors that are convenient for a given set of canonical coordinates. Indeed, for the canonical coordinates r, s in our example, we can consider a prior of the form s(r) = log(r) + log(x T )ζ(r) where the function ζ : (0, ∞) → R satisfies ζ(y 0 ) = 0, ζ(r) ≤ 1 and dζ/dr ≥ 0. For this experiment, the approach of López-Lopera et al. [2017] was used as a prior model for the monotone, bounded function ζ; for brevity we refer the reader to the paper for further detail.
To obtain empirical results we consider the ODE with F (r) = r −1 + r and y 0 = 0.1. The function F was evaluated on a regular grid of n points. The posterior distributions that were obtained as the number n of data points was increased were plotted in the (r, s) plane in Fig. 4 and in the (x, y) plane in Fig. 5. Observe that the implicit prior principle ensures that all curves in the (x, y) plane are well-defined functions (i.e. there is at most one y value for each x value). Observe also that the posterior mass contracts to the true solution of the ODE as the number of evaluations n of the gradient field is increased.

Conclusion
This paper presented a foundational perspective on the solution of ODEs by a PNM. It was first argued that there did not yet exist a Bayesian PNM in this context. Then, to address this gap, a novel Bayesian PNM was developed. The Bayesian perspective that we have put forward sheds light on foundational issues which will need to be addressed going forward: As explained in Section 1.2, existing PNM for ODEs each take the underlying state space Y to be the solution space of the ODE. This appears to be problematic, in the sense that a generic evaluation f (x i , y i ) of the gradient field cannot be cast as information A(y † ) about the solution y † of the ODE unless the point (x i , y i ) lies exactly on the solution curve {(x, y † (x)) : x ∈ [x 0 , x T ]}. As a consequence, all existing PNM of which we are aware operate outside the Bayesian framework. The assumption of a solvable Lie algebra, used in this work, can be seen as a mechanism to ensure the existence of an exact information operator A, so that this problem is avoided. The proposed method was intended as a proof-of-concept and it is therefore useful to highlight the aspects in which it is limited. First, the route to obtain transformations admitted by the ODE demands that some aspects of the gradient field f are known, in contrast to other work in which f is treated as a black-box. Second, the class of ODEs for which a solvable Lie algebra is admitted is relatively small, though references such as Bluman and Anco [2002] document important cases where our method could be applied. Third, the principles for prior construction that we identified do not entail a unique prior and, as such, the question of prior elicitation must still be addressed. Nevertheless, as our proof-of-concept has demonstrated, the development of an exact Bayesian method, or a principled approximation to such a method, is not an insurmountable task.

A Appendix
This appendix contains those definitions and concepts from the theory of Lie groups that are required to understand and implement our proposed Bayesian PNM. In particular, a constructive proof of the key enabling result, Theorem 1 in the main text, is provided.

A.1 One-Parameter Lie Groups of Transformations
Definition 2 (One-Parameter Group of Transformations). A one-parameter group of transformations on a domain D is a map defined on D × S for some set S ⊂ R, together with a bivariate map such that the following hold: (1) For each ∈ S, the transformation X(·, ) is a bijection on D.
(3) If 0 is the identity element in (S, φ), then X(·, 0 ) is the identity map on D.
In what follows we continue to use the shorthand notation x * = X(x, ). The notion of a Lie group additionally includes smoothness assumptions on the maps that constitute a group of transformations. Recall that a real-valued function is analytic if it can be locally expressed as a convergent power series.
Definition 3 (One-Parameter Lie Group of Transformations). Let X, together with φ, form a one-parameter group of transformations on a domain D. Then we say that X, together with φ, form a one-parameter Lie group of transformations on D if, in addition, the following hold: (5) S is a (possibly unbounded) interval in R.
(7) For each x ∈ D, X(x, ·) is an analytic function on S.
Without loss of generality it will be assumed, through re-parametrisation if required, that S contains the origin and = 0 is the identity element in (S, φ). The definition is now briefly illustrated: Example 3 (Rotation Group). The one-parameter transformation for ∈ R forms a Lie group of transformations on D = R 2 with group composition law φ( , δ) = + δ.
Definition 4 (Infinitesimal Transformation). Let X be a one-parameter Lie group of transformations. Then the transformation is called the infinitesimal transformation associated to X and the map ξ is called an infinitesimal.
Definition 5 (Infinitesimal Generator). The infinitesimal generator of a one-parameter Lie group of transformations X is defined to be the operator: where ξ is the infinitesimal associated to X and ∇ = ( ∂ ∂x 1 , ∂ ∂x 2 , . . . , ∂ ∂xn ) is the gradient. Example 4 (Ex. 3, continued). For Ex. 3, we have so the infinitesimal generator for the rotation group is: The first fundamental theorem of Lie provides a constructive route to obtain the infinitesimal generator from the transformation itself: Theorem 2 (First Fundamental Theorem of Lie). A one parameter Lie group of transformations X is characterised by the initial value problem: where τ ( ) is a parametrisation of which satisfies τ (0) = 0 and, for = 0, Here δ −1 denotes the group inverse element for δ.
The next result can be viewed as a converse to Theorem 2, as it shows how to obtain the transformation from the infinitesimal generator: Theorem 3. A one parameter Lie group of transformations with infinitesimal generator X is equivalent to Proof. From Taylor's theorem we have that For any differentiable function F we have that and similarly so that the stated result is recovered.
The following is immediate from the proof of Theorem 3:

A.2 Invariance Under Transformation
In this section we explain what it means for a curve or a surface to be invariant under a Lie group of transformations and how this notion relates to the infinitesimal generator.
Definition 6 (Invariant Function). A function F : D → R is said to be invariant under a one parameter Lie group of transformations Based on the results in Section A.1, one might expect that invariance to a transformation can be expressed in terms of the infinitesimal generator of the transformation. This is indeed the case: Proof. The result is established as follows: where the last line follows since the coefficient of the O( ) term in the Taylor expansion must vanish. This completes the proof.
Theorem 5. For a function F : D → R and a one parameter Lie group of transformations x * = X(x, ), the relation F (x * ) = F (x) + holds for all x ∈ D and ∈ S if and only if XF (x) = 1 for all x ∈ D.
Proof. From Cor. 1, we have that F (x * ) = F (x) + XF (x) + O( 2 ). The result follows from inspection of the coefficient.
The following definition is fundamental to the method proposed in Section 2 of the main text: Definition 7 (Canonical Coordinates). Consider a coordinate system r = (r 1 (x), . . . , r n (x)) on a domain D. Then any one parameter Lie group of transformations x * = X(x, ) induces a transformation of the coordinates r * i = r i (x * ). The coordinate system r is called canonical for the transformation if: r * n−1 = r n−1 r * n = r n + . Example 5 (Ex. 3, continued). For the rotation group in Ex. 3, we have canonical coordinates r 1 (x 1 , x 2 ) = x 2 1 + x 2 2 , r 2 (x 1 , x 2 ) = arctan(x 2 /x 1 ). In canonical coordinates, a one parameter Lie group of transformations can be viewed as a straight-forward translation in the r n -axis. The existence of canonical coordinates is established in Thm. 2.3.5-2 of Bluman and Anco [2002]. Note that Thms. 4 and 5 imply that Xr * i = 0 for i = 1, 2, ..., n − 1, Xr * n = 1. Definition 8 (Invariant Surface). For a function F : D → R, a surface defined by F (x) = 0 is said to be invariant under a one parameter Lie group of transformation x * = X(x, ) if and only if F (x * ) = 0 whenever F (x) = 0 for all x ∈ D and ∈ S.
The invariance of a surface, as for a function, can be cast in terms of an infinitesimal generator: Corollary 2. A surface F (x) = 0 is invariant under a one parameter Lie group of transformations with infinitesimal generator X if and only if XF (x) = 0 whenever F (x) = 0. This is sufficient background on groups of transformations; next we turn to the use of Lie group methods in the ODE context.

A.3 Symmetry Methods for ODEs
The aim of this section is to relate abstract Lie transformations to those ODEs for which these transformations are admitted. These techniques form the basis for our proposed method in Section 2.
For an ODE of the form in Eq.
(1), one can consider the action of a transformation on the coordinates (x, y); i.e. a special case of the above framework where the generic coordinates x 1 and x 2 are respectively the independent (x) and dependent (y) variables of the ODE. It is clear that such a transformation also implies some kind of transformation of the derivative dy/dx; this is made explicit next. Note that in this paper we focus on ODEs, as opposed to PDEs, and as such the independent variable x is scalar. The notation y m := d m y dx m will be used. Consider a one-parameter Lie group of transformations (x * , y * ) = (X(x, y; ), Y (x, y; )).
Then we have from the chain rule that is a function of x, y, y 1 , . . . , y m and we denote y * m = Y m (x, y, y 1 , . . . , y m ; ). As an explicit example: In this sense a transformation defined on (x, y) can be naturally extended to a transformation on (x, y, y 1 , y 2 , . . . ) as required.
Example 6. Clearly any ODE of the form dy dx = F (x) admits the transformation (x * , y * ) = (x, y + ). Our next task is to understand how the infinitesimal generator of a transformation can be extended to act on derivatives y m .
Corollary 3. A differentiable function F : D m → R where D m is the phase space containing elements of the form (x, y, y 1 , . . . , y m ), is invariant under a one parameter Lie group of transformations with an extended infinitesimal generator X (m) if and only if X (m) F (x, y, y 1 , . . . , y m ) = 0 for all (x, y, y 1 , . . . , y m ) ∈ D m .
The reader now has sufficient background to address the case of a first-order ODE discussed in this paper. In particular, we can now present a proof of the key enabling result stated in the main text.
Proof of Theorem 1. Let the infinitesimal generator associated with the Lie group of transformations be denoted X. From the remarks after Def. 7, we can obtain canonical coordinates by solving the pair of first order partial differential equations Xr = 0, Xs = 1. By the chain rule we have ds dr = s x + s y y r x + r y y =: G(r, s) From the definition of canonical coordinates, the Lie group of transformations is r * = r, s * = s + in the transformed coordinate system, so ds * dr * = G(r * , s * ) =⇒ ds dr = G(r, s + ) for all , which implies G(r, s) = G(r) and thus Eq. (3) becomes ds dr = G(r) as claimed.