Monte Carlo construction of cubature on Wiener space

In this paper, we investigate application of mathematical optimization to construction of a cubature formula on Wiener space, which is a weak approximation method of stochastic differential equations introduced by Lyons and Victoir (Cubature on Wiener Space, Proc. R. Soc. Lond. A 460, 169--198). After giving a brief review of the cubature theory on Wiener space, we show that a cubature formula of general dimension and degree can be obtained through a Monte Carlo sampling and linear programming. This paper also includes an extension of stochastic Tchakaloff's theorem, which technically yields the proof of our primary result.


Introduction
Cubature on Wiener space [19] is a certain family of numerical formula for approximating the expectation of functionals of diffusion processes, which are important in mathematical finance and other related fields. The existence of cubature formula of general dimension and degree has been known, but the constructions given in the literature were based on algebraic structure among continuous paths and Brownian motion, and limited to ones of low degree and dimension. Our aim in this paper is to obtain a method of constructing general cubature formula on Wiener space through mathematical optimization.
More concretely, we are interested in approximating values in the form E [ f (X T )], where (X t ) 0≤t≤T is a solution of stochastic differential equation (SDE) driven by a multidimensional Brownian motion, and f is a function with some regularity (e.g., Lipschitz continuity). The standard approach consists of two steps: (1) approximate the distribution of X t+∆t from the information of (approximated) X t , where ∆t ≪ 1; (2) split the time interval [0, T ] by 0 = t 0 < t 1 < · · · < t k = T and sequentially apply (1) )] − E [ f (X T )] = O(1/k) with respect to k, the number of partitions [12].
To achieve higher precision of approximation such as O(1/k 2 ), Lyons and Victoir [19] introduced the cubature on Wiener space, which uses the higher order of stochastic Taylor expansion in step (1). The basic scheme for (1) in cubature on Wiener space is as follows: (c1) approximate the distribution of Brownian motion (B s ) t≤s≤t+∆t by a weighted discrete set of deterministic paths; (c2) locally (over [t,t +∆t]) solve ordinary differential equations (ODEs) driven by the deterministic paths in (c1) instead of the stochastic differential equations for approximating the distribution of X t+∆t ; Among the steps (c1), (c2), (2) in cubature on Wiener space, this paper primarily focuses on finding a good weighted set of deterministic paths described in (c1).
In the field of high-order approximation of SDEs, there is also a relative of cubature on Wiener space called Kusuoka approximation [15] using a random paths in (c1), which is followed by the papers [22,20] presenting concrete second-order schemes. The main challenge shared by these approaches (cubature on Wiener space, Kusuoka approximation) is that constructing such a (random) set of paths has been performed by actually solving Lie-algebraic equations and is limited to low dimension and degree of precision. Therefore, our objective is to find a way of generally constructing such a formula in arbitrary dimension and degree.
Contribution of this study Broadly speaking, the contribution of this study comprises the following two items: -(main contribution) We show that one can construct a cubature formula on Wiener space of general dimension and degree with a randomized algorithm.
-(technical contribution) To apply the technique of [10] to the problem of cubature on Wiener space, we characterize the affine hull of the distribution of iterated Stratonovich integrals and prove stochastic Tchakaloff's theorem in a stronger way.
The main result with a simple Monte Carlo construction is given in Proposition 19. It asserts that a certain random generation of piecewise linear paths combined with a linear programming yields a cubature formula on Wiener space.
As a technical contribution, we extend stochastic Tchakaloff's theorem [19], which assures the existence of cubature formulas on Wiener space. Although the original statement was just that there exists a cubature formula, we show that the expectation of iterated Stratonovich integrals of a Brownian motion (E [ϕ W (B)] in (4)) is contained in the relative interior of conv{ϕ W (w)} with a valid range of bounded variation (BV) paths w (Theorem 17). This stronger statement with "relative interior" follows our characterization (Proposition 15) of aff supp P ϕ W (B) in terms of (4) and is essential in exploiting the existing construction of general cubature formula [10].
We only treat the part (c1) in this study, but it is important to consider combining our construction with ODE solvers (c2) and some techniques reducing computational complexity (2) such as recombination [17], which needs further investigation and is deferred to future research.
Outline We give a brief overview of the following sections.
Section 2 describes the idea and background of this study in a more mathematical way. We briefly explain the concept of cubature on Wiener space in Section 2.1 and the overview of a preceding result recently given by one of the authors [10] in Section 2.2. Section 3 is devoted to a theoretical review of cubature on Wiener space by [19]. We introduce the facts around vector fields and relevant notations in Section 3.1 used throughout the paper. The precise definition and error estimate of cubature on Wiener space is given in Section 3.2, and Section 3.3 provides information about known constructions based on algebraic arguments.
In Section 4, we give the extended statement of stochastic Tchakaloff's theorem and its proof. Sections 4.1 and 4.2 provide algebraic background of cubature construction. Section 4.3 is devoted to the proof of our version of stochastic Tchakaloff's theorem, and it also includes the characterization of the distribution of iterated Stratonovich integrals from the viewpoint of the affine hull of the support.
We discuss a way to obtain piecewise linear cubature formula on Wiener space in Section 5. After giving some general properties of continuous BV paths in our context in Section 5.1, we prove our main result in Section 5.2 that we can generally construct cubature formulas on Wiener space with simple randomized algorithms. Section 5.3 represents the numerical verification of our result in a small range of parameters.
Finally we summarize our conclusion in Section 6.

Preliminaries
In this section, we shall briefly explain the background and motivation of this study. Section 2.1 gives a brief explanation of the cubature theory on Wiener space [19] and identifies the problem to be solved. Section 2.2 introduces the recent Monte Carlo approach by [10] to general cubature construction, which our method in this paper is based on. We explain the relation between cubature on Wiener space and a generalized one, as well as the particular difficulty that arise on the Wiener space.

Cubature on Wiener space
be the space of infinitely differentiable R N -valued functions defined on R N whose every order of derivative is bounded. Let us consider the following N-dimensional Stratonovich SDE: As the process X t is dependent on the initial value x, we denote it by X t (x) if necessary. We may assume the solution X t (x) is continuous with respect to t and x. Our aim is to efficiently compute or approximate the expectation E [ f (X t )] with t > 0 and some smooth or Lipschitz f . This sort of approximation is called a weak approximation of SDE and well-studied in the literature [12,13,15,19,20,22,23].
We here focus on the approach introduced in [19] called cubature on Wiener space. Broadly speaking, a cubature formula on Wiener space (of the time interval [0, T ]) is the approximation where P is the Wiener measure on the Wiener space for j = 1, . . . , n, and λ 1 , . . ., λ n are positive real weights whose sum equals one. Instead of polynomials in conventional cubature formulas, we adopt iterated integrals as the test functionals, that is, we want to find paths w i satisfying over some set of multiindices (i 1 , . . . , i k ), where the iterated Stratonovich integral appears in the left-hand side. Precisely speaking, we formally set B 0 t = t for t ≥ 0 and assume w j is a path of BV for each j = 1, . . . , n. Although w 0 j (t) = t is also assumed in [19], here we may generalize and remove this condition.
The iterated integrals appearing in (2) have a rich algebraic structure (see Section 4.1), so algebraic approaches have been adopted in the literature [19,22,8,24,29,21]. However, solving complicated equations of Lie algebra is required in those approaches, and constructions of the formula are limited to a small range (see Section 3.3). Our objective is to give a construction method for a general setting where there are no limitations on the number of iterations of the integral (k in (2)). For this purpose, we adopt an optimization-based viewpoint instead of algebraic ones and extend to our situation the result of [10], which gives a randomized construction of general cubature formula.
We should mention the Kusuoka approximation [13,15,29], which is closely related to cubature on Wiener space. However, our objective in this paper is limited to the construction of the cubature on Wiener space, and the optimization-based approach to the general Kusuoka approximation is deferred for future work.

Monte Carlo approach to generalized cubature
A cubature formula is originally a numerical integration formula on some Euclidean space that exactly integrates polynomials up to a certain degree [30], the theory of which underlies cubature on Wiener space introduced in the previous section. Then, we shall explain it in a generalized setting and briefly explain the idea of [10] for constructing general cubature formulas.
In the above statement, supp P ϕ(X) is the support of the distribution of the vector-valued random variable ϕ(X ). Equivalently, this is the smallest closed set A ⊂ R D satisfying P(ϕ(X ) ∈ A) = 1. Generalized Tchakaloff's theorem can be understood as an immediate consequence of a discretegeometric argument. Indeed, E [ϕ(X )] is contained in the convex hull of supp P ϕ(X) (the convex hull of an A ⊂ R D is defined as conv . Therefore, the generalized Tchakaloff's theorem follows Carathéodory's theorem (note that we assume ϕ 1 ≡ 1): Theorem 2 (Carathéodory) For an arbitrary A ⊂ R D and x ∈ conv A, there exists D + 1 points x 1 , . . ., x D+1 ∈ A such that x ∈ conv{x 1 , . . ., x D+1 }.
Although the above argument cannot directly be used in the construction of cubature, we can use its nature by introducing the concept of relative interior. For a set A ⊂ R D , its affine hull is defined by Then, the relative interior of A is the interior of A regarding the subspace topology on aff A and denoted by ri A. In terms of this relative interior, the following generalization of Carathéodory's theorem holds: Then, for each x ∈ ri conv A, there exists some subset B ⊂ A composed of at most 2k points such that aff B = aff A and x ∈ ri conv B.
From this generalization and the fact that E [ϕ(X )] ∈ ri conv supp P ϕ(X) holds (see, e.g., [10] or essentially [2] for proof), we obtain the following randomized construction of cubature formulas from i.i.d. copies of X : Though the weights remain undetermined, for a sufficiently large n, it suffices to take a basic feasible solution of the linear programming problem Indeed, its basic feasible solution satisfies the bound of points used in a cubature given in Tchakaloff's theorem (Theorem 1). This sort of technique reducing the number of points in a discrete measure is called Carathéodory-Tchakaloff subsampling [25].
Here, if we formally write the iterated integral 0<t 1 <···<t k <T dw i 1 (t 1 ) · · ·dw i k (t k ) appearing in (2) as ϕ (i 1 ,...,i k ) (w) for a valid w (and the Brownian motion B), then the cubature on Wiener space is a set of paths w j and weights λ j formally satisfying where ϕ W denotes a vector of some functions of the form ϕ (i 1 ,...,i k ) . Therefore, if we could directly generate sample paths of the Brownian motion, then Theorem 4 should be applicable. In reality, it is impossible to generate a Brownian motion on a computer, and it is not even a BV path. However, we assume the following variant, supporting our arguments.

Remark 1
The assumption that X 1 , X 2 , . . . possess the same distribution as X in Theorem 4 can be relaxed; the same conclusion yields from the following condition for the i.i.d. sequence: aff supp P ϕ(X 1 ) = aff supp P ϕ(X) , supp P ϕ(X 1 ) ⊃ supp P ϕ(X) .
From this fact, it is sufficient to investigate the distribution of iterated integrals, and we indeed show the Wiener space counterpart of the condition (5) in Proposition 15.

Theoretical background of cubature on Wiener space
In this section, we provide a theoretical review on the cubature theory on Wiener space, first introduced by [19]. We quickly introduce basic notions concerning multidimensional stochastic flows, and give the error estimate of cubature formula. Moreover, we will demonstrate few examples of concrete construction of cubature formula on Wiener space in Section 3.3.

Vector fields
In this section, we define vector fields on R N and show the correspondence between vector fields and vector-valued functions. Let C ∞ (R N ) be the set of real-valued smooth functions over R N .
Due to this condition, a vector field on R N has to be a differential operator and ∂ i denotes the i-th partial derivative for i = 1, . . ., d. Therefore, a vector field corresponds to the vector-valued smooth function (V 1 , . . .,V N ) ⊤ : R N → R N . By abuse of notation, we also denote this vector-valued function by V .
If A and B are vector fields on R N , we define the Lie bracket is also a vector field because the second derivatives vanish. Note that [A, B] corresponds to the vector where A, B are regarded as functions and ∂C denotes the Jacobian matrix of C (see, e.g., [9]).
Reciprocally, the coefficients of the SDE (1) can be regarded as vector fields. These vector fields are closely related to the behavior of X t .
Let V 0 , . . .,V d be the vector fields (operators) induced by the coefficients of (1), and define the operator L : Let us consider the parabolic partial differential equation with a Lipschitz function f : [11], we can exploit the numerical schemes in PDE theory to get E [ f (X T (x))] and vice versa.
We also introduce several conditions on the vector fields. They are assumed to obtain the estimate given in the Proposition 6. Before we state those, we introduce some notations based on [13].
We also define, for each integer m ≥ 1, We can now state the uniformly finitely generated (UFG) condition [16,13]: (UFG) There exists a positive integer L ≥ 1 such that, for an arbitrary α ∈ A 1 , there exists This is equivalent to the statement that the C ∞ b (R N )-module generated by {V [α] | α ∈ A 1 } is finitely generated. Note that (UFG) is known to be strictly weaker than the uniform Hörmander condition (see, e.g., Example 2 in [14]), which is one of the typical assumptions on vector fields.
Although only the condition (UFG) was assumed in [19], it was pointed out in [5] that the following condition is also essential: .
The following estimate is essential.
We finally state the stochastic Taylor formula in terms of the vector-field notation introduced above. By Itô's formula, we obtain for any where we denote ds by • dB 0 s . Therefore, the repetition of Ito's formula yields This is the stochastic Taylor formula, which is rigorously stated as follows. For a multiindex α = (α 1 , . . . , α k ) ∈ A , we denote by V α the operator V α 1 · · ·V α k .
and m be a positive integer. Then, we have where the remainder term satisfies, for some constant C > 0,

Formulation and evaluation of cubature on Wiener space
We can now precisely define the cubature formula [19].
Definition 8 Let T > 0, and let m be a positive integer. BV paths w 1 , . . . , w n ∈ C 0 0 ([0, T ]; R ⊕ R d ) and weights λ 1 , . . . , λ n ≥ 0 with ∑ n i=1 λ j = 1 to define a cubature formula on Wiener space of degree m at time T , if only holds for all (i 1 , . . . , i k ) ∈ A (m). Here, • dB 0 s represents ds. To construct a cubature formula, it suffices to find it over [0, 1]. Indeed, when w 1 , . . . , with the same weights define the cubature over [0, T ]. This is an immediate consequence of the scaling property of the Brownian motion. Once such paths are given, we can easily compute each evolution driven by w i as it is just an . Indeed, the estimate in Proposition 7 holds for t = T if we replace the Wiener measure by the discrete measure ∑ n j=1 λ j δ w j . Therefore, by applying Cauchy-Schwarz we obtain the evaluation: with the constant C > 0 depending only on w 1 , . . . , w n . The above formula does not work as a good approximation unless T is small. Therefore, we divide [0, T ] into smaller time intervals as 0 = t 0 < t 1 < · · · < t k = T . If we consider the repeated application of the cubature formula over each subinterval [t ℓ−1 ,t ℓ ], we can use where w * v denotes the concatenation of two paths and s ℓ := t ℓ − t ℓ−1 for each ℓ = 1, . . .k, as an approximation of the expectation E [ f (X T (x))]. If we define discrete Markov random variables Y 0 , . . . ,Y k independent of the Brownian motion as coincides with the approximation (9). Then, combining an estimate for with Proposition 6, we can prove the following assertion.
for some constant C > 0, which is dependent only on m and w 1 , . . ., w n .
The equally spaced partition t ℓ = ℓT /k for ℓ = 0, . . . , k is not the optimal one in terms of asymptotic error bound with k → ∞. Consider taking t ℓ = T 1 − 1 − ℓ k γ with a constant γ > 0 independent of k (γ = 1 corresponds to the equally spaced partition). By taking γ > m − 1, we have the following estimate [13]: Therefore, a cubature formula of degree m with an appropriate time partition achieves the error rate O(k −(m−1)/2 ) where k is the number of partitions.
Remark 2 If we have to compute all the k-times concatenation of a cubature formula composed of n sample paths, we have to solve n k+1 −1 n−1 ODEs in total [19]. When number of ODEs is too large, we reduce the computational complexity through some Monte Carlo simulation or subsampling method [17,32,25]. In this paper, we do not consider efficient implementation of concatenation of cubature formula, but only consider their constructions.

Known constructions of cubature on Wiener space
We should note that some concrete examples of cubature formulas on Wiener space are already known. The simplest case treated in [19] is m = 3, where we have a cubature formula composed of linear paths (i.e., with only one linear segment).
Constructions for the case d in general and m = 5 are also given in [19], where the authors give cubature formula using only O(d 3 ) paths. Moreover, [8] constructed higher order cubature formula up to m = 11, but the construction is limited to one dimensional space-time (d = 1). Ref. [21] represents other concrete examples when m = 5 with general d and the case (d, m) = (2, 7).
All the aforementioned examples are derived by solving equations in terms of Lie algebra (see Section 4.1), which can be directly written using the Campbell-Baker-Hausdorff formula. However, as a different approach, we address an optimization-based construction in the next section.

Stochastic Tchakaloff's theorem
In the previous section, we have demonstrated the theoretical support of cubature on Wiener space. However, it is important to know if such formulas can actually be constructed. In this section, we shall state stochastic Tchakaloff's theorem, which assures the existence of cubature formula on Wiener space. Though the stochastic Tchakaloff's theorem is originally given in [19], we state it in a stronger way by using the concept of relative interior.
Before doing so, we shall introduce the rich algebraic structures behind the theory of cubature on Wiener space in the following two sections, which are also essential in our proof of stochastic Tchakaloff's theorem.

Tensor and Lie algebra
We introduce a tensor algebra which is suitable to our case [18,19,15] for each positive integer n. Here, the condition for (i 1 , . . ., i k ) means that A i 1 ⊗ · · · ⊗ A i k takes all the arrangement of R and R d such that 2(# of R) + (# of R d ) = n. Then, we consider the tensor algebra of formal series where the direct sum hereafter is regarded as a series, i.e., T is the set of all the infinite sequences (a n ) ∞ n=0 where a n ∈ U n for each n ≥ 0. Let T (n) (E) := n k=0 U k , and let π n : T ((E)) → T (n) (E) be the canonical projection for each n ≥ 0. As we are only interested in these projections in practice, we do not need to differentiate the usual tensor algebra from that of series treated here [27].
It might be easier to understand T as the ring of formal power series R[[Z 0 , Z 1 , . . ., Z d ]] with noncommutative variables Z 0 , . . ., Z d (see, e.g., [1]). In that case, we redefine the degree of some monomial Y by degY : denotes the number of Z i appearing in Y ) and regard U n (E) as the subspace spanned by monomials of degree n for each n ≥ 0.
For any elements a = (a n ) ∞ n=0 , b = (b n ) ∞ n=0 ∈ T ((E)), we define the sum and product as follows: where the latter two operations are limited for a ∈ T ((E)) with a 0 = 0. Note that these operations commute with each projection homomorphism π n . Let us introduce the space of Lie series. Define a (a ∈ A, b ∈ B). L((E)) is the so-called free Lie algebra generated by E [26]. The elements of L((E)) are called Lie series. We also define L (n) (E) := π n (L((E))), the elements of which are called Lie polynomials.

Signature of a path
We shall introduce the signature (or Chen series [4]) of a path, which summarizes the algebraic structure of iterated integrals. Let w = (w 0 , . . . , w d ) ∈ C 0 0 ([0, T ]; R ⊕ R d ) be a BV path and we define its signature. The integration by dw is a where the integration by dw means the Lebesgue-Stieltjes integration. In both presentations, we think of the 0-th (or constant) term of S(w) s,t as 1. We call S(w) s,t the signature of w over [s,t].
The following is Chen's theorem.

Theorem 11 ([4, 18]) The process S(w) satisfies S(w) s,t ⊗ S(w) t,u = S(w) s,u for arbitrary
It also holds that log S(w) s,t ∈ L((E)) and therefore π n (log S(w) s,t ) ∈ L (n) (E). Moreover, the inverse of this correspondence holds, i.e., for an arbitrary Lie polynomial L ∈ L (n) (E) ⊂ T (E) and arbitrary 0 ≤ s < t ≤ T , there exists a bounded-variation path w ∈ C 0 0 ([0, T ]; R⊕ R d ) such that π n (log S(w) s,t ) = L . Remark 3 Regarding the latter part, a stronger result is known [7, Theorem 7.28]. Every Lie polynomial can be exactly (not approximately) represented as a (truncated) logarithm of some continuous piecewise linear path with a finite number of linear intervals.
By virtue of these assertions, we see that the problem of finding paths constructing a cubature formula is equivalent to the problem of finding the corresponding Lie polynomials. The following Brownian-motion version of this result is also important.

Proposition 12 ([15, 19]) Define the (Stratonovich) signature of the Brownian motion as an element of T ((E)) = R[[Z
Then, log S(B) s,t is almost surely a Lie series.
As we mainly deal with the signature over [0, 1], hereafter let S(w) and S(B) represent S(w) 0,1 and S(B) 0,1 , respectively. We also define for each α = (i 1 , . . . , i k ) ∈ A , Note that we set I α (w) = I α (B) = 1 if α = / 0. If to define L := π n (log S(B)), Indeed, we obtain the expression As L is a random Lie polynomial from the previous assertion, roughly speaking, the generalized Tchakaloff's theorem (Theorem 1) and the inverse statement in Theorem 11 yield the existence of a cubature formula on Wiener space. However, we should point out that the surjectivity stated in Theorem 11 fails if we require that w 0 (t) is monotone (so the original proof of stochastic Tchakaloff's theorem in [19] should be modified).
Proposition 13 Let n ≥ 4. Then, there exists a Lie polynomial L ∈ L (n) (R⊕R) such that π n (exp L ) cannot be expressed as π n (S(w) s,t ) for any BV path w ∈ C 0 0 ([0, T ]; R ⊕ R) with strictly monotone w 0 .
Proof Consider an R ⊕ R-valued continuous BV path w = (w 0 , w 1 ) on [0, T ] that starts at the origin with w 0 strictly increasing. We have Because w 0 is strictly increasing, there exists a differentiation is clearly a Lie polynomial, the proof is complete.

⊓ ⊔
Although the surjectivity fails, we can actually prove the existence of a cubature formula with w 0 (t) = t in the following section. We use the following well-known approximation statement for the Brownian motion.
Proof We only give a sketch as we use arguments based on rough paths. If there is no timeterm (the zero-th entry of the path) and n = 2, then the result yields from the well-known dyadic piecewise-linear approximation for the Brownian rough path (see [6,Proposition 3.6] or [7,Proposition 13.18]). Adding the time (zero-th entry) is not difficult as it is sufficiently smooth and does not affect the regularity of the rough path. To generalize n from n = 2, it suffices to observe the continuity of the "Lyons lift" (also see [7,Chapter 9]). ⊓ ⊔

Proof of stochastic Tchakaloff's theorem
Throughout the section, we fix a positive integer m and consider elements in T (m) (E). Note that T (m) (E) can naturally be regarded in the same light as F := R A (m) . Define a set G (as a subset of F) by is a BV path, w 0 (1) = 1} From Theorem 11, this coincides with the set of exp(L ), where L is a Lie polynomial such that the coefficient of Z 0 is 1.
We denote the distribution of S(B) over F by P S(B) . We shall argue the relation of G and supp P S(B) in the following.
Proposition 15 It holds that aff supp P S(B) = aff G.
Proof From Proposition 12, supp P S(B) ⊂ aff G holds (as aff G includes the closure of G). Therefore, it is sufficient to show G ⊂ aff supp P S(B) .
As aff supp P S(B) is the intersection of all the hyperplanes, which includes supp P S(B) , it can be represented as where H is the family of all (c, d) ∈ F × R such that c ⊤ S(B) = d holds almost surely. The problem is now reduced to the statement for every w appearing in the definition of G. This results from the following lemma as π 0 (S(B)) = π 0 (S(w)) = 1 always holds.

⊓ ⊔
The following is the key lemma in the above proof. We give its proof in the appendix as it is elementary.
Proof By virtue of Carathéodory's theorem, the former part follows from the latter part. We here show the latter part.
From Theorem 3, (3) and Proposition 15, we can find n Lie polynomials L 1 , . . ., L n such that each π m (exp L i ) is contained in G, and E [π m (S(B))] is contained in the relative interior of their convex hull. Here, n can actually be taken such that n ≤ 2 dim G (≤ 2|A (m)|) because of Theorem 3. From the correspondence stated in Chen's theorem (Theorem 11), we can find a desired set of paths in C 0 0 ([0, 1]; R ⊕ R d ). Note that the condition π m (exp L i ) ⊂ G implies that the corresponding path satisfies w 0 i (1) = 1. ⊓ ⊔ Remark 4 By exploiting Proposition 14, we can also prove the same result even if we require w 0 i (t) = t for each i = 1, . . . , n and 0 ≤ t ≤ 1. Although Proposition 14 only asserts an approximation result, "relative interior" argument fills the gap.

Monte Carlo approach to cubature on Wiener space
In this section, we investigate a way to construct cubature formulas, which is based on mathematical optimization instead of Lie algebra. We limit the arguments to cubature formula composed of continuous piecewise linear paths, and propose a construction based on Monte Carlo sampling, which is the application of [10] to our case. We also carry out numerical experiments in concrete cases.
Although existing constructions of cubature formulas on Wiener space are based on Lie-algebraic equations, we can simply regard the cubature construction as an optimization problem. One such way is to consider an LP problem, which is analogous to ordinary cubature problems treated in Section 2.2. From this viewpoint, we can naively generate many sample paths and then reduce their number by using Carathéodory-Tchakaloff subsampling. We later see that this approach is applicable at least theoretically (Section 5.2).

Signature of continuous BV paths
In this section, we see the properties of BV paths and their signature. We also see that the truncated signature of continuous BV paths can be approximated with any accuracy by that of piecewise linear paths.
Let w = (w 0 , w 1 , . . ., w d ) ∈ C 0 0 ([0, 1]; R ⊕ R d ) be BV paths. We define the total variation of w as where ∆ is the partition of [0, 1] by 0 = t 0 < t 1 < · · · < t k = 1 and k varies in sup ∆ . We call w a BV path if w 1 < ∞ holds. Note that other norms are also equivalent as the space R ⊕ R d is finite-dimensional though we are using the sup-norm In this case, the total variation of w can be written as Signature can also be represented by the derivatives as for each multiindex α = (i 1 , . . . , i k ) ∈ A . As a special case of BV paths, we are interested in (continuous) piecewise linear paths, which are easy to implement on computers. Let 0 = s 0 < s 1 < · · · < s n = 1 be a partition of [0, 1]. Then, we can define a path w ∈ C 0 0 ([0, 1]; R ⊕ R d ) which is linear on each interval [s j−1 , s j ] ( j = 1, . . . , n) by determining the slope vector in R ⊕ R d at each interval. For the sake of computation, here we give the calculation of the signature explicitly. Let g j = (g 0 j , g 1 j , . . ., g d j ) be the slope of w over (s j−1 , s j ). Then, for α = (i 1 , . . . , i k ) ∈ A , holds. We can derive this by dividing the integral domain into disjoint segments, which are compatible with the partition 0 = s 0 < s 1 < · · · < s n = 1. If we adopt the notation g α j := g i 1 j · · · g i k j for each α ∈ A , I α (w) can also be written as The latter expression can also be easily derived from Chen's theorem (Theorem 11).

Piecewise linear cubature
The following theorem assures the existence of a cubature formula on Wiener space composed of continuous piecewise linear paths.

Theorem 18
For each positive integer m, there exist n paths w 1 , . . ., w n ∈ C 0 0 (R ⊕ R d ), which are piecewise linear and n positive weights λ 1 , . . . , λ n whose sum is 1 that satisfy n ≤ |A (m)| and The statement still holds even if we require w 0 i (t) = t for i = 1, . . ., n and 0 ≤ t ≤ 1.
Proof From Theorem 17 and Remark 3 (also from Proposition 14 if we want w 0 (t) = t), we can easily deduce that there exist a set of at most 2|A (m)| continuous piecewise linear paths whose truncated (by π m ) signatures convex hull contains E [π m (S(B))] (in its relative interior). Rigorously, we can apply the same argument as the proof of Theorem 4. Finally, by Carathéodory's theorem, n ≤ |A (m)| can actually be achieved.

⊓ ⊔
Based on this theorem, it is sufficient for us to look for cubature formula within piecewise linear paths. Our approach to construction of a piecewise linear cubature is an application of "Monte Carlo cubature construction" [10]. Of course, we are not able to generate a Brownian motion and use it as a candidate for sample points of cubature formulas, because it is not a BV path and it cannot be implemented on computers anyway. However, the methods in [10] are still applicable here as we see in the following proposition.
are independent (also from the ones of zero-th coordinate) and have a density on R, which is positive almost everywhere.
Then, with probability one, there exists an N such that a subset of {w 1 , w 2 , . . . , w N } can construct a cubature on Wiener space of degree m.
Proof Take M as large as we can find (at most) 2|A (m)| piecewise linear paths (denoted byw ℓ ) with at most M linear segments such that conv{π m (S(w ℓ ))} ℓ contains E [π m (S(B))] in its relative interior. The existence of an M is assured by the proof of Theorem 18 (or see Remark 3, Theorem 17, and Remark 4).
From (11), the truncated signature of each w i is a polynomial of random variables w M)). Our assumption assures that these variables take values in every neighborhood of some point with a positive probability. In particular, it implies that where · is the Euclidean norm on T (m) (E) ≃ F, holds for each i, ℓ and ε > 0. Note that the left-hand side probability does not depend on i by i.i.d. assumption. Therefore, the argument in the proof of Theorem 4 holds here again and we obtain the desired assertion.
⊓ ⊔ Remark 5 For the scheme (b), the condition w 0 i (1) = 1 is necessary to assure that each π m (S(w i )) is contained in G defined in Section 4.3.
Note that the generation rule of sample paths in this proposition is just one of infinitely many possible examples. We may alternatively, for instance, directly generate w j i (k/M) independently. Proposition 19 assures only the existence of "some large" M and N, so it might be numerically hard to find cubature formulas from this approach. However, it is beneficial to know that the construction can be reduced at least to the stage of machine power.

Numerical experiments
Explaining our simple numerical method based on a Monte Carlo approach we first describe the algorithm for computing the signature of piecewise linear paths, and then present our Monte Carlo approach and its result in some pairs of (d, m).
Calculation of signature Note that hereafter d and m are regarded as already given parameters. For positive integers M and N, we generate N piecewise linear paths (denoted by w 1 , . . . , w N ) with M intervals of time (see Proposition 19). The time complexity of this paths generation is O(NMd). Then, we address each component of generated paths by In all the experiments, we generated PATH[i, j, k] (with k = 0) so that it follows the centered normal distribution of variance 1/M. We set PATH[i, j, 0] = 1/M for each i, j and denote w i by just writing PATH[i] for each i.
As we consider not so large M in this study, we calculate the signature of generated paths by a simple dynamic programming (Algorithm 1). In the algorithm, we calculate the signature of w i over [0, k/M] for k = 1, . . ., M, by using the expression (11). The time complexity of this algorithm is O(M|A (m)| 2 ), though a pruning of possible multiindices (α, β ) helps a bit.
Algorithm 1 Calculation of (i-th) signature Initialize: Monte Carlo approach In the approach based on Monte Carlo sampling, we simply generate many paths and determine by solving an LP problem whether or not we can construct a cubature formula of desired degree from generated paths.
The part of solving an LP problem was performed using IBM ILOG CPLEX Optimization Studio (https://www.ibm.com/analytics/cplex-optimizer, version 12.10; CPLEX hereafter). From Proposition 19, for a sufficiently large N we can construct a cubature formula using a subset of paths {w 1 , . . . , w N }.
From these results, one may expect that the change in m is more essential than d where we need larger number of partitions and ratio N/|A (m)| as m gets larger, though more experiments are necessary.  In this paper, we have demonstrated that piecewise linear cubature formula can be constructed on Wiener space through a Monte Carlo sampling and an LP problem. Our construction is supported by the technical contribution, which extends stochastic Tchakaloff's theorem using our characterization of the distribution of Stratonovich iterated integrals. We confirmed that for small pairs of (d, m) our algorithm actually works in numerical experiments.
Although we have shown that one can theoretically construct cubature formulas of any dimension and degree, the number of paths used in our construction only attains the Tchakaloff bound, and therefore it requires too much computational cost for large (d, m) in practice. Therefore, we may consider reducing the number of paths by using additional optimization techniques.
Indeed, in the case k-times zeros appear in the suffix of α, we can show inductively 0<t 1 <···<t k <t t 1 0 u r dB i r dt 1 · · · dt k = · · · = 0<t j <···<t k <t so, we can finally apply the following Ito isometry: for progressively measurable and second mean integrable stochastic process u and v, it holds that where i, j ∈ {1, . . ., d} and δ i j is Kronecker's delta.
Therefore, by an inductive argument, it only remains to prove that E I α t I β t = 0 in the case α or β only contain zeros and α + = β + , but this case is trivial as the one is a constant and the other's expectation becomes zero from (13). We now have completed the proof of (14). where all but finite c α equal to zero and A + denotes the set of all multiindices that contain no zeros.