{'fraction_non_alphanumeric': 0.051417361451692005, 'fraction_numerical': 0.0182834722903384, 'mean_word_length': 4.038149831982605, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Derivative-Free Milstein Type Approximation Method for SPDEs covering the Non-Commutative Noise case

Higher order schemes for stochastic partial differential equations that do not possess commutative noise require the simulation of iterated stochastic integrals. In this work, we propose a derivative-free Milstein type scheme to approximate the mild solution of stochastic partial differential equations that need not to fulfill a commutativity condition for the noise term and which can flexibly be combined with some approximation method for the involved iterated integrals. Recently, the authors introduced two algorithms to simulate such iterated stochastic integrals; these clear the way for the implementation of the proposed higher order scheme. We prove the mean-square convergence of the introduced derivative-free Milstein type scheme which attains the same order as the original Milstein scheme. The original scheme, however, is definitely outperformed when the computational cost is taken into account additionally, that is, in terms of the effective order of convergence. We derive the effective order of convergence for the derivative-free Milstein type scheme analytically in the case that one of the recently proposed algorithms for the approximation of the iterated stochastic integrals is applied. Compared to the exponential Euler scheme and the original Milstein scheme, the proposed derivative-free Milstein type scheme possesses at least the same and in most cases even a higher effective order of convergence depending on the particular SPDE under consideration. These analytical results are illustrated and confirmed with numerical simulations.


Motivation
For the approximation of stochastic partial differential equations (SPDEs) with commutative noise, some higher order schemes such as the Milstein schemes in [1,2,10,15], the derivative-free versions [17] and [31], or the Wagner-Platen type scheme [3] were derived and implemented in the last years. Concerning equations that do not need to possess commutative noise, see [4,14,21,24,28] for some applications, it was, however, an open question how to implement a higher order scheme due to the iterated stochastic integrals that are involved and the numerical scheme of choice was so far some Euler scheme, for example, the exponential Euler or the linear implicit Euler, see [8,13,20]. Recently, the authors presented two algorithms to obtain an approximation of such stochastic integrals, see [19]. In [7], the Milstein scheme proposed by A. Jentzen and M. Röckner [10] has been analyzed for noncommutative equations in the case that it is combined with the algorithms proposed in [19]. However, as the main drawback the Milstein scheme requires the evaluation of the derivative of an operator in each time step. This is the reason that its computational complexity increases quadratically w.r.t. the dimension of the state space compared to the Euler scheme with linearly growing computational complexity. In the present paper, we propose a derivative-free numerical scheme to efficiently approximate the mild solution of SPDEs which do not need to have commutative noise, that is, the commutativity condition for all v ∈ H β , u,ũ ∈ U 0 has not to be fulfilled. Our goal is to approximate the mild solution to SPDEs of type dX t = AX t + F (X t ) dt + B(X t ) dW t , t ∈ (0, T ], X 0 = ξ (2) with a scheme that obtains the same temporal order of convergence as the Milstein scheme, however, without the need to evaluate any derivative and with significantly reduced computational complexity which is of the same order of magnitude as for the Euler scheme, i.e., which depends only linearly on the dimension of the state space. For details on the notation, we refer to Section 2.1. In general, the Milstein scheme proposed in [10] applied to (2) (4) for h > 0, t, t + h ∈ [0, T ] and K ∈ N in terms of increments of the approximated Q-Wiener process (W K t ) t∈[0,T ] like in the commutative case, see [10]. Therefore, methods such as the derivative-free Milstein type scheme presented in [18], which was developed based on this assumption, are not applicable to approximate the mild solution of these equations. In [19], we introduced two methods to approximate iterated stochastic integrals with t ≥ 0, h > 0 for some operators Ψ ∈ L(H, L(U, H) U 0 ), Φ ∈ L(U, H) U 0 , and a Q-Wiener process (W t ) t∈[0,T ] of trace class. Therewith, it is possible to implement the Milstein scheme (3) from [10], we refer to [7] for details. However, the evaluation of the derivative in the Milstein scheme is costly. Precisely, the computational cost needed to evaluate this term is of order O(N 2 K) in each time step, see [18,7]. This computational effort can be reduced by one order of magnitude if the derivative is replaced by some customized approximation -see also the detailed discussion of this issue in [18].
In this work, we design a derivative-free numerical scheme to approximate the mild solution of equation (2) which can be combined with any method to simulate the iterated stochastic integrals involved in the scheme, see also [16]. First, we introduce the setting in which we work and state results on the convergence of the proposed scheme -both, with and without an approximation of the iterated integrals. The same theoretical order of convergence as for the Milstein scheme can be obtained. Moreover, we illustrate the advantages of such a higher order derivative-free scheme with a concrete example in Section 3. We combine the scheme with Algorithm 1 presented in [19], which is based on a truncated Fourier series expansion, and derive the effective order of convergence for this schemea concept that combines the theoretical order of convergence with the computational effort based on a cost model introduced in [18]. In terms of this effective order of convergence, the original Milstein scheme (3) is outperformed by the proposed derivative-free Milstein type scheme. Compared to the exponential Euler scheme, the proposed scheme obtains a higher effective order of convergence for a large set of parameter values when combined with Algorithm 1 from [19]. In Section 4, we analyze the mean-square error and the computational cost for the derivative-free Milstein type scheme numerically. The presented simulations confirm a higher effective order of convergence in contrast to the original Milstein scheme and at least the same or even higher effective order of convergence in contrast to the Euler scheme for the examples considered in Section 4. Finally, in Sections 5 and 6, we give some concluding remarks and the proofs for the convergence results.

Approximation of Solutions for SPDEs
In this section, we present a derivative-free Milstein type scheme for SPDE (2) which does not need to have commutative noise. Precisely, we introduce a scheme which can be coupled with an arbitrary method for the approximation of the involved iterated stochastic integrals. For example, when combined with the algorithms introduced in [19] for the simulation of twice-iterated integrals, the theoretical order of convergence of the original Milstein scheme can be maintained.

Framework
Throughout this work, we assume the framework presented in the following. Let (H, ·, · H ) and (U, ·, · U ) denote some separable real-valued Hilbert spaces and let T ∈ (0, ∞) be some fixed time point. Further, let the operator Q ∈ L(U ) be non-negative, symmetric and have finite trace. Then, the subspace U 0 ⊂ U is defined as U 0 = Q 1 2 U . Moreover, we consider some complete probability space (Ω, F, P) and a U -valued Q-Wiener process (W t ) t∈[0,T ] with respect to the filtration (F t ) t∈[0,T ] which fulfills the usual conditions. In terms of the eigenvalues of Q, denoted as η j , with corresponding eigenvectorsẽ j for j ∈ J with some countable index set J forming an orthonormal basis {ẽ j : j ∈ J} of U (see [23]), we obtain the following series representation of the Q-Wiener process, see [23], In this representation, the stochastic processes (β j t ) t∈[0,T ] denote independent real-valued Brownian motions for all j ∈ J with η j = 0. Below, the following notation is used for different sets of linear operators. The space of linear and bounded operators mapping from U to H that are restricted to the subspace U 0 is called (L(U, H) U 0 , · L(U,H) ) with L(U, H) U 0 := {T : U 0 → H | T ∈ L(U, H)}, by L HS (U, H), we denote the set of Hilbert-Schmidt operators mapping from U to H and, finally, we denominate L (2) (U, H) = L(U, L(U, H)) and L For the existence and uniqueness of a mild solution of SPDE (2) and the validity of the proofs of convergence in Section 6, we assume the following conditions.
(A1) The linear operator A : D(A) ⊂ H → H is the generator of an analytic C 0 -semigroup S(t) = e At for all t ≥ 0. We denote the eigenvalues of −A by λ i ∈ (0, ∞) and the corresponding eigenvectors by e i for i ∈ I and some countable index set I, that is, −Ae i = λ i e i for all i ∈ I. Furthermore, let inf i∈I λ i > 0 and let the eigenfunctions {e i : i ∈ I} of −A form an orthonormal basis of H, see [29], and Av = i∈I −λ i v, e i H e i for all v ∈ D(A). We introduce the real Hilbert spaces H r : (A2) Let β ∈ [0, 1) and assume that F : (A3) The operator B : H β → L(U, H) U 0 is assumed to be twice continuously Fréchet differentiable such that sup v∈H β B ′ (v) L(H,L(U,H)) < ∞ and sup v∈H β B ′′ (v) L (2) (H,L(U,H)) < ∞. Further, let B(H δ ) ⊂ L(U, H δ ) for some δ ∈ (0, 1 2 ) and assume that , any projection operator P : H → span{e i : i ∈Ĩ} ⊂ H with finite index setĨ ⊂ I and the case that P is the identity.
(A5) Assume that at least one of the following conditions is fulfilled: and any projection operator P : H → span{e i : i ∈Ĩ} ⊂ H with finite index setĨ ⊂ I.
In this work, we do not make a difference between the operator B and its extensionB : The operatorB is globally Lipschitz continuous as H β ⊂ H is dense. With F , we deal analogously. Note that assumptions (A1)-(A4) are the same as for the scheme for SPDEs with commutative noise introduced in [18] and similar to the conditions imposed in [10] and [7] for the original Milstein scheme. However, the commutativity condition (1), which is essential in [10] and [18], needs not to be fulfilled in our setting. On the other hand, assumption (A5) is required. Assumptions (A1)-(A4) assure the existence of a unique mild solution X : [0, T ] × Ω → H γ for SPDE (2), see [9,10]. Moreover, it holds for every r ∈ [0, γ] and p ∈ [2,4], see [9].

The Derivative-Free Milstein Type Scheme
We derive a numerical scheme to approximate the mild solution of SPDE (2) in this section. At first, the infinite dimensional spaces have to be discretized. For the solution space H, we introduce the projection operator P N : H → H N that maps H to the finite dimensional subspace H N := span{e i : i ∈ I N } for some fixed N ∈ N with some index set I N ⊂ I and |I N | = N . We define this operator as Analogously, we define the projection operator P K : U → U K to approximate the Q-Wiener process for some fixed K ∈ N by with U K := span{ẽ j : j ∈ J K } for some index set J K ⊂ J , |J K | = K, and η j = 0 for j ∈ J K . In order to discretize the time interval, we work with an equidistant time step for legibility of the representation. Let h = T M for some M ∈ N and define t m = m · h for m ∈ {0, . . . , M }. The increments of the approximated Q-Wiener process are then denoted by where the increments of the real-valued Brownian motions are given by ∆β j m = β j t m+1 − β j tm for m ∈ {0, . . . , M − 1}, j ∈ J K . The Milstein scheme (3) is computationally expensive due to the derivative that has to be evaluated in each step, see also [18] and [7]. In order to compare numerical methods in this work, we consider the so-called effective order of convergence first introduced in [27]. This number combines the theoretical order of convergence with the computational cost involved in the calculation of an approximation by a particular scheme. As in [18], the goal is to raise the effective order of convergence by means of a customized approximation that is free of derivatives and, in addition, computationally less expensive. Here, however, we do not need to assume that the operator B ′ B fulfills a commutativity condition, that is, condition (1) is not required. In Theorem 2.1, we state that for the proposed derivative-free Milstein type scheme, the theoretical order of convergence is the same as for the Milstein scheme given in [10]. At the same time, compared to the Milstein scheme in [10], the computational effort is significantly reduced for the derivative-free Milstein type scheme. That means that the effective order of convergence is a priori larger for the proposed derivative-free method. Moreover, compared to the Euler schemes like in [8,13,20], the computational cost is of the same magnitude while the order of convergence w.r.t. step size h is at least the same or even significantly higher. Thus, the scheme that we derive in the following is more efficient in terms of the effective order of convergence than the Euler type schemes for most parameter sets determined by the SPDE under consideration if we combine it, for example, with the algorithms for the simulation of the iterated stochastic integrals introduced in [19], see Table 2. Precisely, compared to the Euler schemes, the increase in the computational cost that results from the approximation of the iterated stochastic integrals can be neglected and we get, in many cases, a significantly higher effective order of convergence due to the higher theoretical order of convergence in the time step that the derivative-free Milstein type scheme features.
The main idea for the derivative-free Milstein type scheme is alike to that in the commutative case, see [18], which in turn is based on the work for the finite dimensional setting in [25,26,27]. The operator B ′ B is approximated by a customized difference operator in such a way that the overall computational cost is decreased by one order of magnitude. However, the stage values have to be chosen differently compared to the commutative case. Compared to [18], the main distinction is that we employ one difference term only. The derivative-free Milstein type scheme yields a discrete process which we denote by At this point, we assume that the iterated stochastic integrals are given exactly in order to consider the error estimate independent of the approximation error for the iterated integrals. We consider the error resulting from the approximation of the mild solution of (2) without an approximation of the iterated stochastic integral since this is interchangeable. Then, in a second step, we conclude from Theorem 2.2 below that if an approximation of the iterated stochastic integral fulfills some specified conditions, this estimate remains valid.
Proof. The proof of Theorem 2.1 is stated in Section 6. This is the same estimate (apart from the constant) as for the Milstein scheme (3) proposed in [10] or the derivative-free Milstein type scheme for SPDEs with commutative noise in [18]. The computational effort, however, increases compared to the schemes for SPDEs with commutative noise as the iterated stochastic integrals have to be simulated. We discuss this issue below.

Approximation of Iterated Integrals
In Section 2.2, we implicitly assumed that the iterated stochastic integrals can be computed exactly. However, up to now there exists no algorithm for the exact simulation of the iterated stochastic integrals in a setting with non-commutative noise. Therefore, the iterated integrals have to be approximated appropriately. We prove the following general result.
Proof. The proof of this theorem is stated in Section 6.
Note that Theorem 2.2 applies to the Milstein scheme (3) as well, see also [7]. Now, we want to illustrate this statement with two exemplary choices -Algorithm 1 and Algorithm 2 as introduced in [19]. First, we consider Algorithm 1 which is based on a series representation of the iterated stochastic integral. This representation is truncated after D summands for some D ∈ N, see [12,19], which yields the approximation. The numerical scheme (7) is called DFM-A1 if the iterated integrals are approximated by Algorithm 1 -denoted asĪ Q (i,j),m =Ī Q,(D), (1) (i,j),m . For this method, there exists some constant C Q,T > 0 such that (10) is fulfilled with for all D, K, M ∈ N, see [19,Corollary 1]. If we approximate the integrals with Algorithm 2 instead, we denote the scheme (7) by DFM-A2 and the approximationĪ Q (i,j),m of I Q (i,j),m byĪ Q,(D), (2) (i,j),m . The series representation is not only truncated after D summands, but the remainder is approximated by a multivariate normally distributed random vector, for details, we refer to [19,32]. For this algorithm, (10) holds with for all D, K, M ∈ N and some constant C Q,T > 0, see [19,Corollary 2,Theorem 4]. This estimate shows that the error converges in D with a higher order, compared to the estimate for Algorithm 1.
Note that the error estimate also depends on the number K, which controls the accuracy of the approximation of the Q-Wiener process, and on the eigenvalues of the operator Q. For a proof of the error estimates (13) and (14), we refer to [19]. Moreover, conditions (11) and (12) are fulfilled for Algorithm 1 and 2, which can be easily seen from the definition of the algorithms in [19].
In order to determine which of the two algorithms obtains a higher order of convergence, one has to analyze the computational costs that are involved, see also [19,7] for a comparison. The goal is that the DFM scheme combined with Algorithm 1 or Algorithm 2 preserves the error estimate stated in Theorem 2.1. This requires a choice of D ≥ D 1 = ⌈M 2 min(2(γ−β),γ)−1 ⌉ for Algorithm 1, whereas for Algorithm 2 we need ) which results in the same orders of magnitude for the choice of D 1 and D 2 , respectively. These considerations show that the computational effort for the two schemes DFM-A1 and DFM-A2 is determined by the parameters which in turn are specified by the equation. Therefore, the choice of the optimal scheme depends on the SPDE that has to be solved. From now on, we assume that D ∈ N is chosen such that the temporal order of convergence is not decreased, i.e., such that D = D 1 for Algorithm 1 or D = D 2 for Algorithm 2, respectively.
Remark 2.1. Note that Algorithm 1 and Algorithm 2 proposed in [19] merely represent examples and that Theorem 2.2 is valid if the derivative-free Milstein type scheme DFM is combined with any approximation for the iterated stochastic integrals such that conditions (10) together with (11) or (12) are fulfilled.

The Effective Order of Convergence -A Comparison
In the following, we compare the performance of the derivative-free Milstein type (DFM) scheme to the performance of the original Milstein (MIL) scheme (3), the exponential Euler (EXE) scheme and the linear implicit Euler (LIE) scheme. For example, one can combine the DFM scheme and the MIL scheme with Algorithm 1 or Algorithm 2 in order to approximate the solution of SPDEs that need not fulfill the commutativity condition (1). However, the analysis can be done similarly for any other approximation method for the iterated stochastic integrals as specified in Theorem 2.2. In the following, we restrict our analysis to Algorithm 1 as an example. The LIE scheme is considered in [13,30] and the EXE scheme is introduced in [20] which are combined with a Galerkin approximation. The proof of the following theorem is detailed in [16] and the main idea can be found in [10].
Compared to the DFM scheme, the EXE scheme requires less restrictive assumptions as we do not need (A5) or conditions on the second derivative of B and the estimate for B ′ (v)P B(v)−B ′ (w)P B(w) can be omitted. For the LIE scheme, similar results as in Proposition 3.1 can be obtained analogously. Below, q stands for the order of convergence w.r.t. the step size However, to compare the performance of the schemes we have to take into account their computational cost in combination with their error estimates as, e.g., iterated stochastic integrals have to be simulated for the higher order schemes DFM and MIL only.

The Cost Model
In order to compare the efficiency of different approximation algorithms, one is usually interested in the dependency of the errors on their computational cost. Therefore, we consider a theoretical cost model proposed in [18]. It is assumed that any standard arithmetic operation or evaluation of sine, cosine or exponential function etc. produces unit cost 1. Further, the simulation of any realization of an N (0, 1)-distributed real valued random variable is assumed to produce cost one as well. However, the evaluation φ(v) of a functional φ ∈ V * with V = H or V = U is assumed to be usually more costly with cost(φ, v) = cost(φ) ≡ c for all v ∈ V and for some c ≥ 1 where typically c ≫ 1. Such functionals are needed for, e.g., the calculation of Fourier coefficients As a result, we obtain for any v, y ∈ H N and u ∈ U K the following computational costs due to |I N | = N and needs the evaluation of the functionals B(y)ẽ j , e i H for i ∈ I N and j ∈ J K which results in the functionals (B ′ (y)e k )ẽ j , e l H have to be evaluated for all k, l ∈ I N and j ∈ J K and it follows that cost( Considering the computational cost for one time step of the Milstein scheme (3), one evaluation of )ẽ j for some basis elementẽ j ∈ U K is for free because they are in the j-th column of the matrix representation )ẽ j , e i H which are already determined. The same applies to the operator  In contrast to the Milstein scheme, in each time step the proposed derivative-free Milstein type scheme DFM needs one evaluation of P N • F (·), one evaluation of P N • B(·)| U K , and the calculation of Observe that the calculation of each summand requires the computation of the functionals for k ∈ I N with cost({φ j k : k ∈ I N }) = cN for each j ∈ J K due to |I N | = N . Therefore, the evaluation of (16) can be done with cost O(N K). Here, the crucial point is that although the argument of B depends on the index j, the resulting operator is then applied to the basis functionẽ j only, which is responsible for the fundamental reduction of the computational complexity. In addition, the linear operators P N •e Ah | H N ∈ L(H, H) N,N and P N •B(Y N,K,M m )| U K ∈ L(U, H) K,N (note again that calculating, e.g., P N B(Y N,K,M m )ẽ j for a basisẽ j ∈ U K is for free) have to be applied to the corresponding elements of the Hilbert spaces. Thus, the total computational cost for M time steps of the DFM scheme is cost(DFM(N, K, M )) = O(N KM ) if the cost for the simulation of iterated stochastic integrals is not taken into account.
Analogously to the derivative-free Milstein type schemes DFM-A1 and DFM-A2, we denote the Milstein scheme by MIL-A1 and MIL-A2 if it is combined with either Algorithm 1 or Algorithm 2 proposed in [19] for the approximation of the iterated stochastic integrals, respectively. Then, the dominating computational cost due to necessary evaluations of real-valued functionals and the simulation of random numbers for each time step can be found in Table 1 for the linear implicit Euler scheme LIE as well as for the EXE, MIL-A1, MIL-A2, DFM-A1 and DFM-A2 schemes. It has to be pointed out that, in contrast to finite-dimensional stochastic differential equations (SDEs), the computational effort of each numerical scheme depends not only on the number of time steps M , but also on the dimensions N and K of the subspaces H N and U K which have to increase in order to decrease the approximation error, compare Theorem 2.1 and Proposition 3.1. However, it turns out that different schemes can attain the same error estimates like in (9) and (15), however with significantly different computational cost, see also the discussion in [18]. In order to compare the performance of different numerical schemes, one has to compare the accuracy of each scheme versus the needed computational cost instead of just comparing their error estimates w.r.t. N , K and M given in Theorem 2.1 and Proposition 3.1. For example, the Milstein scheme and the derivative-free Milstein type scheme both attain the same error estimate (9), however, the computational effort for MIL is cost(MIL(N, K, M )) = O(N 2 KM ) whereas for the DFM scheme it is only cost(DFM(N, K, M )) = O(N KM ) if the random numbers are assumed not to be the dominating cost. Therefore, in this case, the DFM scheme performs a priori with a higher order of convergence compared to the MIL scheme if errors versus costs are considered. Compared with the LIE scheme and the EXE scheme, the DFM scheme belongs to the same class O(N KM ) of computational complexity, which is in some sense optimal for one-step approximations for SPDEs of type (2). Although the LIE scheme as well as the EXE scheme have worse error bounds given in (15) compared to the one for the DFM and the MIL scheme in (9), it is not clear which scheme should be preferred because the computational cost for simulating the iterated stochastic integrals for the DFM and the MIL scheme have to be taken into account as well. Therefore, we derive the effective order of convergence for each scheme under consideration. This concept is also detailed in [18].

Comparison of the Effective Orders of Convergence
In order to compare the performance of different numerical schemes, we consider the so-called effective order of convergence which was proposed in [27] and also considered in [18]. In the following, we restrict our comparison to the schemes DFM-A1 and MIL-A1, both using Algorithm 1, as well as the EXE scheme in order to keep the analysis concise. For a detailed analysis and comparison of the effective order of convergence for the schemes MIL-A1, MIL-A2 and EXE we refer to [7]. Since the LIE scheme and the EXE scheme have the same order of convergence and similar computational cost, we restrict our analysis to the EXE scheme in the following because one can get exactly the same results for the LIE scheme. We want to point out that the focus of this article lies on the introduction and analysis of the derivative-free Milstein type scheme and a complete comparison taking into account further algorithms next to Algorithm 1 for the simulation of the iterated stochastic integrals would go beyond the scope of this article and may be object of future research.  Table 1 depends on N , K and M only. In the following, we assume that Then, we obtain the following expression for all N, K, M ∈ N and some C > 0, see also [18], Note that the parameter q > 0 is determined by the scheme that is considered. Given some computational costc > 0, the goal is to minimize the error under the constraint that the computational cost is bounded byc. Solving this optimization problem yields the effective order of convergence, denoted by EOC(SCHEME), which is then given by an expression of the form Next, we analyze the effective order of convergence for the DFM-A1 scheme and the MIL-A1 scheme, which make use of Algorithm 1 for the approximation of the iterated stochastic integrals, and the EXE scheme. Therefore, let q := q DFM = q MIL = min(2(γ − β), γ) and let D = O(M 2q−1 ) for Algorithm 1 in the following.
First, we consider the scheme DFM-A1. The computational cost for the calculation of one trajectory amounts toc = O(M KN ) + O(KM 2q ), see Table 1 and the discussion in the last section. Now, two cases have to be distinguished: We solve the optimization problem and obtain Further, the effective order of convergence is given by which is the same result as for the derivative-free Milstein type scheme in the case of SPDEs with commutative noise, see the computations in [18]. On the other hand, if γρ A (2q − 1) ≥ q holds, then c = O(KM 2q ) and optimization yields In this case, we obtain the effective order of convergence from As a result of this, we obtain the effective order of convergence from which is also the same effective order of convergence as for the Milstein scheme if it is applied to some SPDE with commutative noise, see also [18]. However, in the case of γρ A (2q − 1) ≥ 2q the computational effort for the MIL-A1 scheme isc = O(KM 2q ) and we obtain the same choice for M , N and K as given in (19) and also the same effective order of convergence as given by (20), see also [7].
Finally, we consider the EXE scheme where the optimal choice for M , N and K is given by (17), however, with q = q EES = min( 1 2 , 2(γ − β), γ) and the effective order of convergence for the EXE scheme was computed in [18] and is given by Here, we note that the same holds for the LIE scheme.
In order to determine the scheme which is most efficient for the approximation of the solution for an SPDE of type (2) that needs not to fulfill a commutativity condition for the noise, we have to compare the effective orders of convergence according to the distinct parameter settings for the schemes DFM-A1, MIL-A1 and EXE.
If γρ A (2q − 1) ≤ q and q ≤ 1 2 , then, it follows that q = q EES . Thus, the EXE scheme and the DFM-A1 scheme have the same effective order of convergence given in (18) and (23), whereas the MIL-A1 scheme obviously has a lower effective order of convergence given in (22).

Optimal scheme
Optimal M , N , K EOC αρ Q 2αρ Q +1 Table 2: For a given parameter set, the conditions in this table have to be checked in order to determine the optimal scheme among the schemes DFM-A1, MIL-A1 and EXE. Here, q = q DFM = q MIL .
If γρ A (2q − 1) ≤ q and q > 1 2 , then, it follows that q > q EES = 1 2 . Here, the DFM-A1 scheme has obviously a higher effective order of convergence compared to the one of the EXE scheme. Further, comparing the effective order of convergence of the EXE scheme and the MIL-A1 scheme results in Here, it follows that the EXE scheme has a higher order of convergence than the MIL-A1 scheme and thus the DFM-A1 scheme attains the highest effective order of convergence in this case.
If q ≤ γρ A (2q − 1) ≤ 2q, then, it follows that q > q EES = 1 2 . In this case, it holds for the effective orders of convergence of the EXE, the MIL-A1 and the DFM-A1 scheme that Thus, the DFM-A1 scheme is the one with the highest effective order of convergence in the present case.
If 2q ≤ γρ A (2q − 1), it holds that q > q EES = 1 2 . In this case, the DFM-A1 and the MIL-A1 scheme attain the same effective order of convergence given in (20). As a result of this, a comparison of the effective order of the EXE scheme with the one of the MIL-A1 scheme and DFM-A1 scheme results in In this case, the DFM-A1 scheme and the MIL-A1 scheme have the same effective order of convergence which is higher than the one of the EXE scheme.
The same holds true if the EXE scheme is replaced by the LIE scheme. We summarize the results of our comparison in Table 2 which shows that the DFM-A1 scheme always attains the highest possible effective order of convergence. However, in the case of q DFM = q EES ≤ 1 2 , although the EXE scheme and the DFM-A1 scheme have the same effective order of convergence, one may prefer the EXE scheme because it requires less computational effort compared to the DFM-A1 scheme, see Table 1. On the other hand, in the case of 2q ≤ γρ A (2q − 1), both the DFM-A1 and the MIL-A1 scheme have the same optimal effective order of convergence, which is higher than that of the EXE scheme. Here, one may prefer the DFM-A1 scheme because it needs less computational effort compared to the MIL-A1 scheme, see Table 1, and because it is derivative-free whereas one has to calculate the derivative of the operator B for the MIL-A1 scheme.
Finally, it has to be pointed out that the maximal effective order of convergence that can be attained is always bounded by 1/2 independent of the given parameters whenever Algorithm 1 is applied to simulate the iterated stochastic integrals.
For completeness, we want to note that assumption (A5) as well as parts of (A3) do not have to be fulfilled for the exponential Euler scheme. This means that there might be parameter sets that are valid for the EXE scheme but not for the DFM scheme and in these situations the exponential Euler scheme would be the method of choice. Moreover, it is not clear if the obtained upper error bounds are sharp and thus if the effective order of convergence may be further improved.

The Case of a Finite-Dimensional Q-Wiener Process
If the Q-Wiener process W is finite-dimensional, i.e., if |{j ∈ J : η j = 0}| < ∞, the error estimate only depends on M and N provided we choose K = |{j ∈ J : η j = 0}|. Then, we obtain new solutions for M and N solving the optimization problem that minimizes the error under the constraint of a prescribed computational cost budgetc. Therefore, we compare once more the DFM-A1 scheme, the MIL-A1 scheme and the EXE scheme. Now, the computational cost required to approximate one trajectory of the solution of SPDE (2)  If γρ A (2q − 1) ≤ q, the computational cost for the DFM-A1 scheme isc = O(M N ) and solving the optimization problem yields Then, the effective order of convergence is given by If γρ A (2q − 1) ≥ q, then q > 1 2 andc = O(M 2q ) for the DFM-A1 scheme. Here, optimization results in and the effective order of convergence can be calculated as Considering the MIL-A1 scheme, again two cases have to be distinguished: If γρ A (2q − 1) ≤ 2q, the MIL-A1 has computational costc = O(M N 2 ) and optimization yields and the effective order of convergence is given by If 2q ≤ γρ A (2q − 1), it follows that q > 1 2 and that the MIL-A1 scheme attains the same computational costc = O(M 2q ) as the DFM-A1 scheme in the second case. Thus, we also get (26) for M and N , and also the same effective order of convergence as given by (27).
Clearly, for the EXE scheme it holds q EES ≤ 1 2 and thus we get the same results as for the DFM-A1 scheme given in (24) for M and N as well as by (25) for the effective order of convergence with q = q EES .
Finally, comparing the effective orders of convergence for the schemes under consideration, we easily derive the results presented in Table 3. Here, again the DFM-A1 scheme performs better or at least as good as one of the other schemes. Clearly, in the case of q DFM = q MIL = q EES ≤ 1 2 one may prefer the EXE scheme or the LIE scheme although they have the same effective order of convergence as the

Optimal scheme
Optimal M , N EOC  Table 3: In case of a finite-dimensional Q-Wiener process and K = |{j ∈ J : η j = 0}| < ∞, the conditions in this table have to be checked in order to determine the optimal scheme among the schemes DFM-A1, MIL-A1 and EXE. Here, q = q DFM = q MIL .
DFM-A1 scheme because they are easier to implement. However, in the case of 2q ≤ γρ A (2q − 1) where the DFM-A1 scheme and the MIL-A1 scheme attain the same effective order of convergence one may prefer the DFM-A1 scheme because it needs less computational effort and because no derivative of the operator B is needed by the DFM-A1 scheme. Again, the effective order of convergence is always bounded by 1/2 as for the infinite-dimensional noise case.

Numerical Analysis
In this section, we compare the DFM-A1 scheme to the MIL-A1 and the EXE schemes to demonstrate the theoretical results presented above, summarized in Tables 2 and 3. Here, we approximate the mild solution of SPDE (2), that is, For the numerical analysis, we consider the following setting. We fix H = U = L 2 ((0, 1), R), set T = 1, and I = J = N. Let A be the Laplace operator with Dirichlet boundary conditions. To be precise, A = ∆ 100 with eigenvalues λ i = π 2 i 2 100 of −A and eigenvectors e i = √ 2 sin(iπx) for i ∈ N, x ∈ (0, 1) and on the boundary, we have X t (0) = X t (1) = 0 for all t ∈ (0, T ]. The covariance operator Q is defined by the eigenvalues η j = j −ρ Q for some ρ Q > 1 which is given separately for each example below and e j = √ 2 sin(jπx) for j ∈ N, x ∈ (0, 1). For the operator B, we present the general setting introduced for the numerical analysis in [18]. Define the functionals µ ij : H β → R, φ k ij : H β → R for i, k ∈ I, j ∈ J such that φ k ij is the Fréchet derivative of µ ij in direction e k and let as well as for y ∈ H β and u, v ∈ U 0 , see [18, Sec. 5.3] for details.
We choose µ ij (y) = y,e j H i p +j 4 for all i ∈ I, j ∈ J , y ∈ H and some p > 1 that differs in the examples presented below, which leads to φ k ij (y) = 0, k = j 1 i p +j 4 , k = j for all i, k ∈ I, j ∈ J , y ∈ H. This is the setting considered in [7]. The assumptions (A1), (A2), and (A4) are obviously fulfilled. We only elaborate on (A3). By the definition of the L(U, H δ )-norm and the operator B, we obtain In the next steps, we employ the Parseval equality and the triangle inequality It holds by Parseval's equality that for all j ∈ J . As | u,ẽ j U | 2 ≤ 1 by Parseval, we obtain Then, for some ε ∈ (0, 2δ), some C 1 = C 1 (ε, δ) > 0 and with r = If for δ ∈ (0, 1 2 ) it holds that p > 2+8δ 3+2δ , then it follows that B(y) L(U,H δ ) ≤ C(1+ y H δ ) for all y ∈ H δ .
Next, we compute the term for all y ∈ H γ . We rewrite the expression above to obtain Here, we employed the definition of the operators A and Q. In the next step, we insert the definition of the operator B By Parseval's equality and calculations as in (30), we obtain for some C 2 > 0 that Then, for all ε ∈ (4ϑ − 1, 4ϑ − 1 + 2p) such that r = 2p 1+ε−4ϑ > 1, q = 2p 2p−1−ε+4ϑ > 1, Young's inequality yields that and ε ∈ (max(0, 4ϑ − 1), 4ϑ − 1 + 2p). In the following examples, p > max( 2+8δ 3+2δ , 1−4ϑ 2 ) and ρ Q are specified and we select γ and α to be maximal. We do not state any other condition given in (A3) as these do not pose a restriction on the parameters but note that these are fulfilled as well. Finally, we examine the commutativity condition (1). On the one hand, it holds that for all y ∈ H and all i ∈ I, m, n ∈ J . Obviously, these two expressions differ for some choice of m, n ∈ J . Thus, the considered example does not fulfill the commutativity condition (1).
With these parameters, we can identify the scheme that is superior. For this example, it holds that q < γρ A (2q − 1) ≤ 2q for sufficiently smallε > 0. Thus, the DFM-A1 scheme is optimal, i.e., it is the scheme with the highest effective order of convergence according to Table 2. In order to compare the DFM-A1 scheme to the other schemes under consideration, we calculate the effective orders of convergence for each of the schemes. We expect that the scheme DFM-A1 obtains the highest effective order of convergence in this setting with given by (20), i.e., EOC(DFM-A1) ≈ 27 58 . Moreover, we fix some arbitrary N ∈ N and compute the relation M = N 2 and K = N as given in (19) for the implementation of the DFM-A1 scheme.
Considering the scheme MIL-A1, the effective order of convergence for this scheme is given by (22) with i.e., EOC(MIL-A1) ≈ 189 460 . For this example, the relations between N , K and M for the MIL-A1 scheme given in (21) are exactly the same as for the DFM-A1 scheme.
For the EXE scheme, on the other hand, we obtain from (17)  as an optimal choice. The effective order of convergence for the EXE scheme is given as as stated in (23), i.e., it holds EOC(EXE) ≈ 189 514 .
As a result of this, for this example, it holds that EOC(EXE) < EOC(MIL-A1) < EOC(DFM-A1) and thus the DFM-A1 scheme performs better than the other schemes. For the numerical evaluation, we compare the schemes DFM-A1, MIL-A1 and EXE to an approximation computed with the linear implicit Euler scheme with N = 2 6 , K = ⌈2 14 9 ⌉ and M = ⌈2 35 2 ⌉ that serves as the reference solution. We simulate 500 paths with each scheme and each prescribed computational cost to compare the meansquare error versus computational cost, see Figure 1. Then, the slope indicates the effective order of convergence if log-log-scales are used. This confirms that for this example the DFM-A1 scheme performs significantly better than the MIL-A1 scheme and the EXE scheme. The results are also stated in Table 4.
Checking the conditions in Table 2, we are in the case of q > 1 2 and γρ A (2q − 1) ≤ q for sufficiently smallε > 0. In this case, the optimal effective order of convergence is obtained by the DFM-A1 scheme according to Table 2. For the DFM-A1 scheme, we get from (18)   i.e., it holds that EOC(MIL-A1) ≈ 1309 3900 . Given some N ∈ N, the optimal choice for M and K is given in (21) and yields the same results as for the DFM-A1 scheme in this example. Now, if we compare the effective orders of convergence for the schemes under consideration, then we have EOC(MIL-A1) < EOC(EXE) < EOC(DFM-A1). In this example, again the DFM-A1 scheme performs best with the highest effective order of convergence and the original Milstein scheme MIL-A1 has the lowest effective order of convergence, which is even less than that of the EXE scheme. As in Example 1, we substitute the exact solution with an approximation computed by the linear implicit Euler scheme. Precisely, we choose N = 2 6 , K = ⌈2 102 77 ⌉ and M = ⌈2 85 6 ⌉ for the computation of the reference solution and compute 500 paths. Table 5 and Figure 2 show the difference in the effective order of convergence between the derivative-free Milstein type scheme DFM-A1, the Milstein scheme MIL-A1 and the exponential Euler scheme EXE.

Example 3
Compared to the first two examples, we choose a different nonlinearity for the third example in order to obtain restrictions for the parameter β. Therefore, we consider the mapping F : H β → H given by for v ∈ H β with some f i : H β → R for i ∈ I. In this example, we choose f i (v) = i −s sin(i r v, e i H ) for v ∈ H β , s > 1 2 , r ≤ min(s, 2β + s 2 ) and i ∈ I. Then, we get Further, F is twice continuously Fréchet differentiable and it holds because r ≤ s. Further, considering the second Fréchet derivative, we get 100 2 i∈I i 4β | z, e i H | 2 for any z ∈ H β and because r ≤ min(s, 2β + s 2 ). Thus, assumption (A2) is fulfilled.
Since we have q ≤ 1 2 , the optimal schemes are the EXE scheme and the DFM-A1 scheme, both attaining the same effective order of convergence for this example, see Table 2. Taking into account all parameters, we get from (18) and (23)  ⌉ as the optimal choice. In this case, the computation of the double integrals is not expensive as it holds that D ≥ D 1 = M − 1 2 −ε for some ε > 0 such that D = 1 can be fixed or it can even be neglected. On the other hand, for the MIL-A1 scheme, we compute from (22) that which gives us the effective order of convergence EOC(MIL-A1) ≈ 7 36 . Moreover, the optimal choice for M and K given some N ∈ N can be calculated from (21) to be exactly the same as for the EXE scheme and the DFM-A1 scheme. For this example, we have EOC(MIL-A1) < EOC(EXE) = EOC(DFM-A1). For this example , the computational effort involved in computing a convergence plot is very high due to the relation M ≈ N 8 . Therefore, we do not present a convergence plot for this setting.
The examples presented above confirm the theoretical analysis that we conducted in Section 3. The numerical experiments show that the derivative-free Milstein type scheme for equations with noncommutative noise defined in (7), in combination with Algorithm 1, has always at least the same and it in many cases an even higher effective order of convergence compared to the exponential Euler scheme and the original Milstein scheme.

Conclusion
We proposed the derivative-free Milstein type scheme DFM for the approximation of the mild solution of SPDEs that need not fulfill a commutativity condition for the noise and we proved an upper bound for the mean-square error. As the main novelty, the introduced DFM scheme is derivative-free and has computational cost O(N KM ) which is of the same magnitude as for the Euler schemes EXE and LIE. This is a significant reduction of the computational complexity compared to the original Milstein scheme MIL that is not derivative-free and which has computational cost O(N 2 KM ). In addition, the convergence of the DFM method is proved if it is combined with any suitable simulation method for the iterated stochastic integrals. As an example, the effective order of convergence of the DFM scheme combined with Algorithm 1 in [19] for the simulation of the iterated stochastic integrals is analyzed in detail. For Algorithm 1, the effective order of convergence of the DFM scheme is at least that for the Euler schemes or the Milstein scheme MIL and turns out to be even significantly higher for many parameter settings depending on the specific SPDE to be approximated. Thus, in many cases the proposed DFM scheme outperforms the Euler schemes as well as the original Milstein scheme.
The maximal possible effective order of convergence that can be attained by the DFM scheme combined with Algorithm 1 is bounded by 1/2, which is in accordance with the upper bound for the order of strong convergence in case of finite-dimensional SDEs if Algorithm 1 is applied, see also [5]. However, in contrast to the finite-dimensional SDEs setting, for SPDEs the Euler schemes often attain some effective order of convergence less than 1/2. This gap in the order of convergence for the Euler schemes is the reason why the use of higher order approximation methods can be reasonable and which is in strong contrast to the finite-dimensional SDE setting. To the best of the authors knowledge, this is the first attempt to give a rigorous analysis of the error versus computational cost for higher order approximation methods applied to SPDEs without any commutativity condition where the computational cost for the approximation of iterated stochastic integrals is incorporated within the framework of a cost model. It remains an open question whether the application of higher order numerical methods that incorporate further iterated stochastic integrals from the stochastic Taylor expansion may close the gap for the order of convergence to the upper bound of 1/2 if, e.g., naive approximations like Algorithm 1 are applied for the approximation of these iterated stochastic integrals. As a result of this, higher order approximation methods may be of strong interest, especially in the case of SPDEs. On the other hand, it may be possible to overcome the upper bound of 1/2 for the order of convergence if some more sophisticated algorithm for the simulation of the iterated stochastic integrals is combined with the DFM scheme, see e.g., Algorithm 2 in [19].

Proofs
Here, we give the proof of the convergence result for the derivative-free Milstein scheme (7) as stated in Theorem 2.1. Moreover, we prove the estimate given in Theorem 2.2 which incorporates the approximation of the stochastic double integrals additionally. In the following, we always denote Y m = Y N,K,M m for simplicity and let I (i,j),l = (η i η j ) − 1 2 I Q (i,j),l . Attention should be paid to the fact that for ease of notation the constants in our proofs may differ from line to line even though their denomination is not changed. We need the following estimate on the moments of the approximation process (Y m ) m∈{0,...,M } for the proof of Theorem 2.1. Note that, without loss of generality, we present the proofs with an equidistant time step h = h m for all m ∈ {0, . . . , M − 1}.
Case 1: Assume that assumption (A5a) is fulfilled. We estimate the individual terms by a Burkholder-Davis-Gundy type inequality [6,Theorem 4.37], and a Taylor expansion of the difference approximation. Precisely, we use The estimates on the analytic semigroup, see Lemma 6.3 and 6.13 in [22,Ch.2], and assumptions (A2), (A3), yield the difference approximation. Thus, we use 1]. With estimates on the analytic semigroup, see Lemma 6.3 and 6.13 in [22,Ch.2], we get that Making use of p = 2, assumptions (A2), (A3) and (A5b) yield j∈J K for all m ∈ {0, . . . , M }, M ∈ N to align the components with the corresponding terms in the approximation below. We define the following auxiliary processes for m ∈ {0, . . . , M }, M, N, K ∈ N The discrete process (Y m ) m∈{0,...,M } denotes the approximation obtained by the DFM scheme in (7). The auxiliary processes are introduced in order to split the approximation error such that we can employ some known prior estimates. We analyze the following terms separately for all m ∈ {0, . . . , M }, M ∈ N. The first term is similar to the error that results from the approximation of (2) with the Milstein scheme by Jentzen and Röckner presented in [10]. A slight difference arises as we introduce the projection operator P N in the definition of Y MIL m , see the computations in [17,18]. The main reasoning, however, is the same. In the error analysis in [10], the commutativity condition is not needed -it is only employed to facilitate implementation -whereas all conditions required in the proof in [10] are fulfilled due to assumptions (A1)-(A4). Therefore, the estimate for arbitrary N, M, K ∈ N is valid. For details, we refer to [10].
The error estimate of the second term in (33), H , m ∈ {0, . . . , M }, M ∈ N, can be obtained by the same means as in the proof of convergence of the Milstein scheme in [10] and mainly relies on the Lipschitz properties of the operators. We transfer this reasoning from [10, Section 6.3], which yields for all m ∈ {0, . . . , M }, M, N, K ∈ N.
Next, we analyze the third term in (33) which represents the error that results from the approximation of the derivative. We can show that the theoretical order of convergence that the Milstein scheme obtains is not reduced by this approximation. We rewrite the expression in (33) as We employ a Taylor approximation of first order for the second term, see (32), such that the first order derivatives cancel. Moreover, the triangle inequality and assumption (A3) imply Case 1: Assume that assumption (A5a) is fulfilled, i.e., Lemma 6.1 is valid for any p ≥ 2. Thus, it follows from (36) that Then, we obtain with Lemma 6.1 in the case that condition (A5a) is valid that by the distributional properties of I (i,j),l , l ∈ {0, . . . , m − 1}, i, j ∈ J K for all m ∈ {1, . . . , M }, M, K ∈ N, see [11].
Case 2: If assumption (A5b) is fulfilled, then Lemma 6.1 is valid for p = 2. Therefore, we need a customized proof to proceed and we get for p = 2 from (32) and (36) that Due to E I Q (i 1 ,j 1 ),l I Q (i 2 ,j 2 ),l = 1 2 η i 1 η j 1 h 2 l if i 1 = i 2 and j 1 = j 2 and E I Q (i 1 ,j 1 ),l I Q (i 2 ,j 2 ),l = 0 otherwise, we get For the terms inside the two integrals, we employ Taylor approximations of first order of the difference operators as in (32) where ξ(Y l , j, r) = Y l + r i∈J K B(Y l )ẽ i I Q (i,j),l for all j ∈ J K , l ∈ {0, . . . , m − 1} and r ∈ [0, 1]; belowξ(Ȳ l , j, r) is defined analogously. This yields . (43) The first term is the error that results from the approximation of the iterated stochastic integral.
Depending on the choice of the scheme, this error estimate may differ. Assumption (10) states that