extfieldsofstudy:
Computer Science
Mathematics

text:
Private optimization in the interpolation regime: faster rates and hardness results

In non-private stochastic convex optimization, stochastic gradient methods converge much faster on interpolation problems -- problems where there exists a solution that simultaneously minimizes all of the sample losses -- than on non-interpolating ones; we show that generally similar improvements are impossible in the private setting. However, when the functions exhibit quadratic growth around the optimum, we show (near) exponential improvements in the private sample complexity. In particular, we propose an adaptive algorithm that improves the sample complexity to achieve expected error $\alpha$ from $\frac{d}{\varepsilon \sqrt{\alpha}}$ to $\frac{1}{\alpha^\rho} + \frac{d}{\varepsilon} \log\left(\frac{1}{\alpha}\right)$ for any fixed $\rho>0$, while retaining the standard minimax-optimal sample complexity for non-interpolation problems. We prove a lower bound that shows the dimension-dependent term is tight. Furthermore, we provide a superefficiency result which demonstrates the necessity of the polynomial term for adaptive algorithms: any algorithm that has a polylogarithmic sample complexity for interpolation problems cannot achieve the minimax-optimal rates for the family of non-interpolation problems.

Introduction
We study differentially private stochastic convex optimization (DP-SCO), where given a dataset S = S n 1 iid ∼ P we wish to solve minimize f (x) = E P [F (x; S)] = Ω F (x; s)dP (s) subject to x ∈ X , (1) while guaranteeing differential privacy. In problem (1), X ⊂ R d is the parameter space, Ω is a sample space, and {F (·; s) : s ∈ S} is a collection of convex losses. We study the interpolation setting, where there exists a solution that simultaneously minimizes all of the sample losses.
Interpolation problems are ubiquitous in machine learning applications: for example, least squares problems with consistent solutions [SV09, NWS14], and problems with over-parametrized models where a perfect predictor exists [MBB18, BHM18, BRT19]. This has led to a great deal of work on the advantages and implications of interpolation [SST10, CSSS11, BHM18, BRT19].
For non-private SCO, interpolation problems allow significant improvements in convergence rates over generic problems [SST10, CSSS11, MBB18, VBS19, WS21]. For general convex functions, [SST10] develop algorithms that obtain O( 1 n ) sub-optimality, improving over the minimax-optimal rate O( 1 √ n ) for non-interpolation problems. Even more dramatic improvements are possible when the functions exhibit growth around the minimizer, as [VBS19] show that SGD achieves exponential rates in this setting compared to polynomial rates without interpolation. [AD19, ACCD20, CCD22] extend these fast convergence results to model-based optimization methods. Despite the recent progress and increased interest in interpolation problems, in the private setting they remain poorly understood. In spite of the substantial progress in characterizing the tight convergence guarantees for a variety of settings in DP optimization [BST14, BFTT19, FKT20, AFKT21, ALD21], we have little understanding of private optimization in the growing class of interpolation problems.
Given (i) the importance of differential privacy and interpolation problems in modern machine learning, (ii) the (often) paralyzingly slow rates of private optimization algorithms, and (iii) the faster rates possible for non-private interpolation problems, the interpolation setting provides a reasonable opportunity for significant speedups in the private setting. This motivates the following two questions: first, is it possible to improve the rates for DP-SCO in the interpolation regime? And, what are the optimal rates?

Our contributions
We answer both questions. In particular, we show that 1. No improvements in general (Section 3): our first result is a hardness result demonstrating that the rates cannot be improved for DP-SCO in the interpolation regime with general convex functions. More precisely, we prove a lower bound of Ω( d nε ) on the excess loss for pure differentially private algorithms. This shows that existing algorithms achieve optimal private rates for this setting.
[BST14] study the closely related problem of differentially private empirical risk minimization (DP-ERM) where the goal is to minimize the empirical loss, and obtain (minimax) optimal rates of d/nε for pure DP and d log(1/δ)/nε for (ε, δ)-DP. Recently, more papers have moved beyond DP-ERM to privately minimizing the population loss (DP-SCO) [BFTT19, FKT20, AFKT21, ADF + 21, BGN21,ALD21]. [BFTT19] was the first paper to obtain the optimal rate 1/ √ n + d log(1/δ)/nε for (ε, δ)-DP, and subsequent papers develop more efficient algorithms that achieve the same rates [FKT20, BFGT20]. Moreover, other papers study DP-SCO under different settings including non-Euclidean geometry [AFKT21, ADF + 21], heavy-tailed data [WXDX20], and functions with growth [ALD21]. However, to the best of our knowledge, there has not been any work in private optimization that studies the problem in the interpolation regime.
[SST10] propose algorithms that roughly achieve the rate 1/n + f ⋆ /n for smooth and convex functions where f ⋆ = min x∈X f (x). In the interpolation regime with f ⋆ = 0, this result obtains loss 1/n improving over the standard 1/ √ n rate for non-interpolation problems. Moreover, [VBS19] studied the interpolation regime for functions with growth and show that SGD enjoys linear convergence (exponential rates). More recently, several papers investigated and developed acceleration-based algorithms in the interpolation regime [LB20, WS21].

Preliminaries
We begin with notation that will be used throughout the paper and provide some standard definitions from convex analysis and differential privacy.
Notation We let n denote the sample size and d the dimension. We let x denote the optimization variable and X ⊂ R d the constraint set. s are samples from Ω, and S is an Ω-valued random variable. For each sample s ∈ Ω, F (·; s) : R d → R ∪ {+∞} is a closed convex function. Let ∂F (x; s) denote the subdifferential of F (·; s) at x. We let Ω n denote the collection of datasets S = (s 1 , . . . , s n ) with n data points from Ω. We let f S (x) := 1 n s∈S F (x, s) denote the empirical loss and f (x) := E[F (x; S)] denote the population loss. The distance of a point to a set is dist (x, Y ) = min y∈Y x − y 2 . We use Diam(X ) = sup x,y∈X x − y 2 to denote the diameter of parameter space X and use D as a bound on the diameter of our parameter space.
Definition 2.1. A randomized mechanism M is (ε, δ)-differentially private ((ε, δ)-DP) if for all datasets S, S ′ ∈ Ω n that differ in a single data point and for all events O in the output space of M , we have We define ε-differential privacy (ε-DP) to be (ε, 0)-differential privacy.
We now recall a couple of standard convex analysis definitions.

A function
Equivalently, a function is L-Lipschitz if ∇f (x) 2 ≤ L for all x ∈ X .

A function
h is H-smooth if it has H-Lipschitz gradient: for all x, y ∈ X ∇h(x) − ∇h(y) 2 ≤ H x − y 2 .

A function
We formally define interpolation problems: Definition 2.3 (Interpolation Problem). Let X ⋆ := argmin x∈X f (x). Then problem (1) is an interpolation problem if there exists x ⋆ ∈ X ⋆ such that for P -almost all s ∈ Ω, we have 0 ∈ ∂F (x ⋆ ; s).
Interpolation problems are common in modern machine learning, where models are overparameterized. One simple example is overparameterized linear regression: there exists a solution that minimizes each individual sample function. Classification problems with margin are another example.
Crucial to our results is the following quadratic growth assumption: Definition 2.4. We say that a function f satisfies the quadratic growth condition if for all x ∈ X This assumption is natural with interpolation and holds for many important applications including noiseless linear regression [SV09, NWS14]. Past work ([VBS19, WS21]) uses this assumption with interpolation to get faster rates of convergence for non-private optimization.
Finally, the adaptivity of our algorithms will crucially depend on an innovation leveraging Lipchitizian extensions, defined as follows. (2) The Lipschitzian extension (2) essentially transforms a general convex function into an L-Lipschitz convex function. We now present a few properties of the Lipschitzian extension that are relevant to our development.
Lemma 2.1. Let f : X → R be convex. Then its Lipschitzian extension satisfies the following: We use the Lipschitzian extension as a substitute for gradient clipping to ensure differential privacy. Unlike gradient clipping, which may alter the geometry of a convex problem to a non-convex one, the Lipschitzian extension of a function remains convex and thus retains other nice properties that we leverage in our algorithms in Section 4.

Hardness of private interpolation
In non-private stochastic convex optimization, for smooth functions it is well known that interpolation problems enjoy the fast rate O(1/n) [SST10] compared to the minimax-optimal O(1/ √ n) without interpolation [Duc18]. In this section, we show that such an improvement is not generally possible with privacy. The same lower bound of private non-interpolation problems, d/nε, holds for interpolation problems.
To state our lower bounds, we present some notation that we will use throughout of the paper. We let S denote the family of function F and dataset S pairs such that F : X × Ω → R is convex and H-smooth in its first argument, |S| = n, and f S (y) = 1 n s∈S F (y, s) is an interpolation problem (Definition 2.3). We define the constrained minimax risk to be where M (ε,δ) be the collection of (ε, δ)-differentially private mechanisms from Ω n to X . We use M (ε,0) to denote the collection of ε-DP mechanisms from Ω n to X . Here, the expectation is taken over the randomness of the mechanism, while the dataset S n is fixed.
We have the following lower bound for private interpolation problems; the proof is deferred to Appendix B.1.
Recall the optimal rate for pure DP optimization problems without interpolation is O( 1 √ n + d nε ). The first term is the non-private rate, as this is the rate one would get if ε = 0. The second term is the private rate, as this is the price algorithms have to pay for privacy. In modern machine learning, problems are often high dimensional, so we often think of the dimension d scaling with some function of the number of samples n. Thus, the private rate is often thought to dominate the non-private rate. For this reason, in this section, we focus on the private rate. The lower bounds of Theorem 1 show that it is not possible to improve the private rate for interpolation problems in general. Similarly, for approximate (ε, δ)-DP, the lower bound shows that improvements are not possible for d = 1. For completeness, as we alluded to earlier, we note that our results do not preclude the possibility of improving the non-private rate from O(1/ √ n) to O(1/n). We leave this as an open problem of independent interest for future work.
Despite this pessimistic result, in the next section we show that substantial improvements are possible for private interpolation problems with additional growth conditions.

Faster rates for interpolation with growth
Having established our hardness result for general interpolation problems, in this section we show that when the functions satisfy additional growth conditions, we get (nearly) exponential improvements in the rates of convergence for private interpolation.
Our algorithms use recent localization techniques that yield optimal algorithms for DP-SCO [FKT20, ALD21] where the algorithm iteratively shrinks the diameter of the domain. However, to obtain faster rates for interpolation, we crucially build on the observation that the norm of the gradients is decreasing as we approach the optimal solution, since ∇F (x; s) 2 ≤ H x − x ⋆ 2 . Hence, by carefully localizing the domain and shrinking the Lipschitz constant accordingly, our algorithms improve the rates for interpolating datasets. However, this technique alone yields an algorithm that may not be private for non-interpolation problems, violating that privacy must hold for all inputs: the reduction in the Lipschitz constant may not hold for non-interpolation problems, and thus, the amount of noise added may not be enough to ensure privacy. To solve this issue, we use the Lipschitzian extension (Definition 2.5) to transform our potentially non-Lipschitz sample functions into Lipschitz ones and guarantee privacy even for non-interpolation problems.
We begin in Section 4.1 by presenting our Lipschitzian extension based algorithm, which recovers the standard optimal rates for (non-interpolation) L-Lipschitz functions while still guaranteeing privacy when the function is not Lipschitz. Then in Section 4.2 we build on this algorithm to develop a localization-based algorithm that obtains faster rates for interpolation-with-growth problems. Finally, in Section 4.3 we present our final adaptive algorithm, which obtains fast rates for interpolation-withgrowth problems while achieving optimal rates for non-interpolation growth problems.

Lipschitzian-extension based algorithms
Existing algorithms for DP-SCO with L-Lipschitz functions may not be private if the input function is not L-Lipschitz [BFGT20, FKT20, ALD21]. Given any DP-SCO algorithm M L (ε,δ) , which is private for L-Lipschitz functions, we present a framework that transforms M L (ε,δ) to an algorithm which is (i) private for all functions, even ones which are not L-Lipschitz functions and (ii) has the same utility guarantees as M L (ε,δ) for L-Lipschitz functions. In simpler terms, our algorithm essentially feeds M L (ε,δ) the Lipschitzian-extension of the sample functions as inputs. Algorithm 1 describes our Lipschitzianextension based framework.

Algorithm 1 Lipschitzian-Extension Algorithm
Require: Dataset S = (s 1 , . . . , s n ) ∈ S n ; 1: Let F L (x; s i ) be the Lipschitzian extension of F (x; s i ) for all i.
For this paper we consider M L (ε,δ) to be Algorithm 2 of [ALD21] (reproduced in Appendix A.2 as Algorithm 5). The following proposition summarizes our guarantees for Algorithm 1.
Proof For the first item, note that Lemma 2.1 implies that For the second item, Lemma 2.1 implies that F L = F when F is L-Lipschitz. Thus, in Algorithm 1, we apply M L (ε,δ) over F itself.
While clipped DP-SGD does ensure privacy for input functions which are not L-Lipschitz, our algorithm has some advantages over clipped DP-SGD: first, clipping does not result in optimal rates for Algorithm 2 Domain and Lipschitz Localization algorithm Require: Dataset S = (s 1 , . . . , s n ) ∈ S n , Lipschitz constant L, domain X , probability parameter β, initial point x 0 1: Set L 1 = L, D 1 = Diam(X ) and X 1 = X 2: Partition the dataset into T partitions (denoted by {S k } T k=1 ) of size m each; S k = (s (k−1)m+1 , . . . , s km ) 3: for i = 1 to T do 4: Set L i+1 = HD i+1 8: end for 9: return the final iterate x T pure DP, and second, clipped DP-SGD results in time complexity O(n 3/2 ). In contrast, our Lipschitzian extension approach is amenable to existing linear time algorithms [FKT20] allowing for almost linear time complexity algorithms for interpolation problems. Finally, while clipping the gradients and using the Lipschitzian extension both alter the effective function being optimized, only the Lipschitzian extension is able to preserve the convexity of said effective function (see item 2 in Lemma 2.1). We make a note about the computational efficiency of Algorithm 1. Recall that when the objective is in fact L-Lipschitz, computing gradients for the Lipschitzian extension (say in the context of a first-order method) is only as expensive as computing the gradients for the original function. In particular, one can first compute the gradient of the original function and use item 4 of Lemma 2.1; when the problem is L-Lipschitz, ∇f (x) 2 is always less than or equal to L and thus the gradient of the Lipschitzian extension is just the gradient of the original function.

Faster non-adaptive algorithm
Building on the Lipschitzian-extension framework of the previous section, in this section, we present our epoch based algorithm, which obtains faster rates in the interpolation-with-growth regime. It uses Algorithm 1 with M L (ε,δ) as Algorithm 5 (reproduced in Appendix A.2) as a subroutine in each epoch, to localize and shrink the domain as the iterates get closer to the true minimizer. Simultaneously, the algorithm also reduces the Lipschitz constant, as the interpolation assumption implies that the norm of the gradient decreases for iterates near the minimizer. The detailed algorithm is given in Algorithm 2 where D i denotes the effective diameter and L i denotes the effective Lipschitz constant in epoch i.
The following theorem provides our upper bounds for Algorithm 2, demonstrating near-exponential rates for interpolation problems; we present the proof in Appendix C.
Theorem 2. Assume each sample function F is L-Lipschitz and H-smooth, and let the population function f satisfy quadratic growth (Definition 2.4). Let Problem (1) be an interpolation problem.
log n , T = n/m and any µ > 0, Algorithm 2 returns x T such that For δ > 0, β = 1 n µ , m = 256 log 2 n H log(1/β) , T = n/m and any µ > 0, The exponential rates in Theorem 2 show a significant improvement in the interpolation regime over the minimax-optimal O(( 1 ALD21]. To get the linear convergence rates, we run roughly n/ log n epochs with log n samples each. Thus, each call of the subroutine runs the algorithm on only logarithmic number of samples compared to the number of epochs. Intuitively, growth conditions improves the performance of the sub-algorithm, while growth and interpolation conditions reduce the search space. This in tandem leads to faster rates.
To better illustrate the improvement in rates compared to the non-private setting, the next corollary states the private sample complexity required to achieve error α in the interpolation regime.
Corollary 4.1. Let the conditions of Theorem 2 hold. For δ = 0 , Algorithm 2 is ε-DP and requires As the sample complexity of DP-SCO to achieve expected error α on general quadratic growth problems is [ALD21] Corollary 4.1 shows that we are able to improve the polynomial dependence on 1/α in the sample complexity to (nearly) logarithmic for interpolation problems.
Remark 1. In contrast to Corollary 4.1, we can tune the failure probability parameter β to get the sample complexity d ε log 2 1 α . Even though this sample complexity does not have the polynomial factor, it may be worse than 1 α ρ + d ε log 1 α , because generally the dimension term is the dominant one. We end this section by considering growth conditions that are weaker than quadratic growth.
Remark 2. (interpolation with κ-growth) We can extend our algorithms to work for the weaker κ- We present the full details of these algorithms in Appendix C.1 (see Algorithm 6). In this setting, we obtain excess loss for interpolation problems, improving over the minimax-optimal loss for non-interpolation problems which is As an example, when κ = 3, this corresponds to an improvement from roughly (d/nε) 3/2 to (d/nε) 3 . Like our previous results, we are again able to show similar improvements for (ε, δ)-DP with better dependence on the dimension. Finally, we note that we have not provided lower bounds for the interpolation-with-κ-growth setting for κ > 2. We leave this question as a direction for future research.

Adaptive algorithm
Though Algorithm 2 is private and enjoys faster rates of convergence in the interpolation regime, it is not necessarily adaptive to interpolation, i.e. it may perform poorly given a non-interpolation problem. In fact, since the shrinkage of the diameter and Lipschitz constants at each iteration hinges squarely on the interpolation assumption, the new domain may not include the optimizing set X ⋆ in the noninterpolation setting, so our algorithm may not even converge. Since in general we do not know a priori whether a dataset is interpolating, it is important to have an algorithm which adapts to interpolation.
To that end, we present an adaptive algorithm that achieves faster rates for interpolation-withgrowth problems while simultaneously obtaining the standard optimal rates for general growth problems. The algorithm consists of two steps. In the first step, our algorithm privately minimizes the objective without assuming it is an interpolation problem. Next, we run our non-adaptive interpolation algorithm of Section 4.2 over the localized domain returned by the first step. If our problem was an interpolating one, the second step recovers the faster rates in Section 4.2. If our problem was not an interpolating one, the first localization step ensures that we at least recover the non-interpolating convergence rate. We stress that the privacy of Algorithm 3 requires that the call to Algorithm 2 remains private even if the problem is non-interpolating. This is ensured by using our Lipschitzian extension based algorithm with M L (ε,δ) as Algorithm 5. The Lipschitzian extension allows us to continue preserving privacy. We present the full details of this algorithm in Algorithm 3.
The following theorem (Theorem 3) states the convergence guarantees of our adaptive algorithm (Algorithm 3) in both the interpolation and non-interpolation regimes for the pure DP setting. The results for approximate DP are similar and can be obtained by replacing d with d log(1/δ); we give the full details in Appendix C.
Theorem 3. Let each sample function F be L-Lipschitz and H-smooth, and let the population function f satisfy quadratic growth (Definition 2.4) with coefficient λ. Let x adapt be the output of Algorithm 3. Then 1. Algorithm 3 is ε-DP.

Without any additional interpolation assumption, x adapt satisfies
3. Let problem (1) be an interpolation problem. Thenx adapt satisfies Proof The privacy of Algorithm 3 follows from the privacy of Algorithms 1 and 2 and post-processing.
To prove the convergence guarantees, we first need to show that the optimal set X ⋆ is in the shrinked domain X int . Using the high probability guarantees of Algorithm 1, we know that with probability 1 − β/2, we have Using the quadratic growth condition, we immediately have x ⋆ − x 1 2 ≤ D int /2 and hence X ⋆ ⊂ X int .
Using smoothness, we have that for any x ∈ X int , Since Algorithm 2 always outputs a point in its input domain (in this case X int ), even in the noninterpolation setting that In the interpolation setting, the guarantees of Algorithm 2 hold and result is immediate.

Optimality and Superefficiency
We conclude this paper by providing a lower bound and a super-efficiency result that demonstrate the tightness of our upper bounds. Recall that our upper bound from Section 4 is roughly (up to constants) 1 for any arbitrarily large c. We begin with an exponential lower bound showing that the second term in (5) is tight. We then prove a superefficiency result that demonstrates that any private algorithm which avoids the first term in (5) cannot be adaptive to interpolation, that is, it can not achieve the minimax optimal rate for the family of non-interpolation problems. Theorem 4 below presents our exponential lower bounds for private interpolation problems with growth. We use the notation and proof structure of Theorem 1. We let S λ ⊂ S be the subcollection of function, data set pairs which also have functions f S n that have λ-quadratic growth (Definition 2.4). The proof of Theorem 4 is found in Appendix D.1.
This lower bound addresses the second term of (5); we now turn to our superefficiency results to lower bound the first term of (5). We start with defining some notation and making some simplifying assumptions. For a fixed function F : X , Ω → R which is convex, H-smooth with respect to the first argument, let S L λ (F ) be the set of datasets S of n data points sampled from Ω such that f S (x) := 1 n s∈S n F (x, s) is L-Lipschitz and have λ-strongly convex objectives. For simplicity, we will assume that 1. inf x∈X F (x; s) = 0 for all s ∈ Ω, 2. X = [−D, D] ⊂ R, and 3. the codomain of F is R + . With this setup, we present the formal statement of our result; the proof of Theorem 5 is found in Appendix D.2.
Theorem 5. Suppose we have some S ∈ S L λ (F ) with L = 2HD such that (F, S) satisfy Definition 2.3. Suppose there is an ε-DP estimator M such that for some t > 0 and absolute constant c. Then, for sufficiently large n, there exists another dataset S ′ ∈ S L λ (F ), where (F, S ′ ) may not satisfy Definition 2.3, such that To better contextualize this result, suppose there exists an algorithm which atttains a exp(− Θ (nε/d)) convergence rate on interpolation problems; i.e., the algorithm is able to avoid the 1/n c term in (5). Then Theorem 5 states that there exists some strongly convex, non-interpolation problem on which the aforementioned algorithm will optimize very poorly; in particular, the algorithm will only be able to return a solution that attains, on average, constant error on this "hard" problem. More generally, recall that in the non-interpolation quadratic growth setting, the optimal error rate is on the order of 1/(nε) 2 [ALD21]. Theorem 5 shows that attaining better-than-polynomial error complexity on quadratic growth interpolation problems implies that the algorithm cannot be minimax optimal in the non-interpolation quadratic growth setting. Thus, the rates our adaptive algorithms attain are the best we can hope for if we want an algorithm to perform well on both interpolation and non-interpolation quadratic growth problems.

A.2 Algorithms from [ALD21]
A.3 Theoretical results from [ALD21] We first reproduce the high probability guarantees of Algorithm 4 as proved in [ALD21].

5:
Letx i be the output of the optimization algorithm.
Now, we reproduce the high probability convergence guarantees of Algorithm 5.
Theorem 6. Let β ≤ 1/(n + d), D 2 (X ) ≤ D and F (x; s) be convex, L-Lipschitz for all s ∈ Ω. Assume that f has κ-growth (Assumption 2.4) with κ ≥ κ > 1. Setting T = 2 log n κ−1 , Algorithm 5 is ε-DP and has with probability 1 − β Theorem 7. Let β ≤ 1/(n + d), D 2 (X ) ≤ D and F (x; s) be convex, L-Lipschitz for all s ∈ Ω. Assume that f has κ-growth (Assumption 2.4) with κ ≥ κ > 1. Setting T = 2 log n κ−1 and δ > 0, Algorithm 5 is (ε, δ)-DP and has with probability 1 − β . We select V to be a γ-packing (with respect to the ℓ 2 norm) of diameter D ball contained in X . Define the separation between v, v ′ ∈ V with respect to the loss f v and f v ′ by For the sake of contradiction, suppose that where inequality (ii) follows from the definition of the separation, and (iii) follows from privacy and the disjoint nature of the events in the union. Rearranging, we get that τ ≥ kHγ 2 8n(1 + e kε (|V| − 1) −1 ) , which is a contradiction. By standard packing inequalities [Wai19], we know that |V| ≥ (D/2γ) d . Setting k = d/ε and γ = D/2e and using the fact that x/(x − 1) is decreasing in x gives We now prove the (ε, δ)-DP lower bound. Consider the following sample risk function We define the datasets S n v := {0} n−k ∪ {v} k inducing the corresponding population risk We select two points v, v ′ contained within the diameter D ball contained in X such that |v − v ′ | = D. Define the separation between v, v ′ ∈ V with respect to the loss f v and f v ′ as For the sake of contradiction, suppose that where inequality (ii) follows from the definition of the separation, and (iii) follows from group privacy of (ε, δ)-privacy [DR14]. Rearranging, we get that which is a contradiction. Setting k = 1/ε and using the fact δ ≤ εe ε−1 /2 gives the first result.

C Proofs from Section 4
We first prove a lemma that each time we shrink the domain size, the set of interpolating solutions still lies in the new domain with high probability, and the new Lipschitz constant we define is a valid Lipschitz constant for the loss defined on the new domain. We prove it in generality for κ-growth.
Lemma C.1. Let X ⋆ denote the set of interpolating solutions of problem (1). Then X ⋆ ⊂ X i for all i ∈ [T ] with probability 1 − β, and ∇F (y; s) 2 ≤ L i for all y ∈ X i .
Proof We prove this lemma for the case when δ = 0, the case when δ > 0 follows similarly. For epoch i, using Theorem 2 of [ALD21], we have with probability 1 − β/T , Using the growth condition on f (·), we have Using c κ = 2(C κ κ) 1/κ , we get x i − x ⋆ 2 ≤ D i+1 /2 with probability 1 − β/T . Thus, for each epoch i, with probability 1 − β/T , each point in the set X ⋆ of optimizers lies in the domain X i . Using a union bound on all epochs, we have X ⋆ ⊂ X i for all i ∈ [T ] with probability 1 − β.
We now prove the second part of the lemma. Using the smoothness of F (·; s) and that ∇F (x ⋆ ; s) = 0 for all x ⋆ ∈ X ⋆ , we have We now restate and prove the convergence rate of Algorithm 2 Theorem 2. Assume each sample function F is L-Lipschitz and H-smooth, and let the population function f satisfy quadratic growth (Definition 2.4). Let Problem (1) be an interpolation problem.
Then Algorithm 2 is (ε, δ)-DP. For δ = 0, β = 1 n µ , m = 256 log 2 n H log(1/β) log n , T = n/m and any µ > 0, Algorithm 2 returns x T such that (3) , T = n/m and any µ > 0, Algorithm 2 returns x T such that Proof First we prove the privacy guarantee of the algorithm. Each sample impacts only one of the iteratesx i , thus Algorithm 2 satisfies the same privacy guarantee as Algorithm 5 by postprocessing. We divide the utility proof into 2 main parts; first is to check the validity of the assumptions while applying Algorithm 5 and second is using its high probability convergence guarantees to get the final rates. To check this, we ensure that the optimum set lies in the new domain X i at step i and that the Lipschitz constant L i defined with respect to the domain is a valid lipschitz constant. This follows from Lemma C.1.
Next, we use the high probability convergence guarantees of the subalgorithm Algorithm 5 to get convergence rates for Algorithm 2. We prove it for the case when δ = 0, the case when δ > 0 is similar. We know that Thus we have

Algorithm 6 Epoch based epoch based epoch based clipped-GD
Require: number of epochs: T , samples in each round: m = n/T , Diameter at the start: D 1 , lipschitz constant at the start L 1 , domain X 1 , initial pointx 0 1: for i = 1 to T do 2:x i ← Output of Algorithm 5 when run on domain X i (diameter D i ), with lipschitz constant L i using m samples.

3:
if δ = 0 then 4: Set L i+1 = HD i+1 10: end for 11: return the final iterate x T Proof The privacy guarantee follows from the proof of Theorem 2. We divide the utility proof into 2 main parts; first is to check the validity of the assumptions while applying Algorithm 5 and second is using its high probability convergence guarantees to get the final rates. To check this, we ensure that the optimum set lies in the new domain defined at every step and that the lipschitz constant defined with respect to the domain is a valid lipschitz constant. This follows from Lemma C.1.
Next, we use the high probability convergence guarantees of the subalgorithm Algorithm 5 to get convergence rates for Algorithm 2.
We prove it for the case when δ = 0, the case when δ > 0 is similar. We know that Thus, we have We note that for T ∼ log n, 1 (κ−1) T −1 ≈ 1 n log(κ−1) and thus for large n, we ignore the terms of the form a − 1 n log(κ−1) since they are ≈ 1. Ignoring these terms by including an additional constant C ′ we can write Using Theorem 2 of [ALD21] on the last epoch, we have with probability 1 − β that Choosing T = log n and m = n/ log n, we have log(log n/β) log 3/2 (n/ log n) n/ log n , d log(log n/β) log(n/ log n) εn/ log n κ κ−2 L κ n 1 .
Now we write results in terms of sample complexity required to achieve a particular error. The sufficient number of samples.
Proof First, we note that the privacy of Algorithm 3 follows from the privacy of Algorithm 2 and Algorithm 1 and post-processing.
To prove the convergence guarantees, we first need to show that the optimal set X ⋆ is included in the shrinked domain X int . Using the high probability guarantees of Algorithm 1, we know that with probability 1 − β/2, we have Using the quadratic growth condition, we immediately have x ⋆ − x 1 2 ≤ D int /2 and hence X ⋆ ⊂ X int . Using smoothness, we have that for any x ∈ X int , Since Algorithm 2 always outputs a point in its input domain (in this case X int ), even in the noninterpolation setting we have that In the interpolation setting, the guarantees of Algorithm 2 hold and the result is immediate.

D Proofs from Section 5 D.1 Proof of Theorem 4
The proof is exactly the same as Theorem 1, except we set k = λn H to ensure that f v (x) for any v ∈ X has λ-quadratic growth. Finally we set γ = D 2 exp( −λnε Hd ) and use the fact that e λnε H ≥ 2 and the fact that x/(x − 1) is decreasing in x to give the desired lower bound.

D.2 Proof of Theorem 5
The proof of this result hinges on the two following supporting propositions. We first copy Proposition 2.2 from [AD20] (listed as Proposition 4 below) in our notation for convenience. We then state Proposition 5 which gives upper and lower bounds on the modulus of continuity (defined in Proposition 4). We note that the lower bound presented in Proposition 5 is one of the novel contributions of this paper; the proof of Proposition 5 can be found in Appendix D.2.1. We will first assume this to be true and prove Theorem 5 before returning prove its correctness.
Proposition 4. For some fixed F : X , Ω → R which is convex and H-smooth with respect to its first argument, let S ∈ S L λ (F ) for L = 2HD. Let Define the corresponding modulus of continuity Assume the mechanism M is ε-DP and for some γ ≤ Proposition 5. Let F : X , Ω → R be convex and H-smooth in its first argument and satisfying inf x∈X F (x; s) = 0 for all s ∈ Ω. Suppose we have some S ∈ S L λ (F ) with L = 2HD which also induces an interpolation problem (a problem which satisfies Definition 2.3). With respect to the dataset S, the modulus of continuity ω(S, 1/ε) satsifies D nε ≤ ω(S, 1/ε) ≤ 8HD λnε With these two results, we can now prove Theorem 5. Restating the conditions of the theorem formally, suppose for some constants c 0 and c 1 there is an ε-DP estimator M such that If t > 1, set t = min(1, t), then the bound certainly still holds for large enough n. If we let x ⋆ S = argmin x∈X f S (x), using the definition of strong convexity, we have that there exists some c 2 and c 3 such that To satisfy the expression from Proposition 4, we select γ such that γω(S; 1/ε) 2 = c 2 De −c3(nε) t .
Using Proposition 5 we must have λnε 4H c 2 exp(−c 3 (nε) t ) ≤ γ ≤ 2nεc 2 exp(−c 3 (nε) t ). Using Proposition 4, we have that Before performing a further lower bound on this quantity, we first verify that log(1/2γ) 2ε does not exceed the total size of the dataset, n. Using our bounds on γ, we see that For any t ∈ (0, 1], for sufficiently large n, this quantity is less than n. We now lower bound the modulus of continuity by using the fact that it is a non-decreasing function in its second argument: This is the desired result; the last inequality comes from another application of Proposition 5 but with c3(nε) t −log c2−log(4nε) 2ε in place of 1/ε.
By smoothness and strong convexity we have However, this implies that γ > λ − H nε which is a contradiction.
The next lemma is a standard result on the closure under addition of strongly convex functions.
This lemma provides some growth conditions on the gradient under smoothness, strong convexity and quadratic growth.
Proof We note that by first order optimality conditions, for all x ⋆ ∈ X ⋆ , ∇g(x ⋆ ) = 0. To prove the first inequality, we have that for any x ⋆ ∈ X ⋆ , the following is true: In particular, minimizing over x ⋆ on the right hand side and rearranging gives the desired result. To prove the second result, we know that by strong convexity for any To prove the last result, we know that by smoothness for any Minimizing over x ⋆ on the right hand side gives the desired result.
This lemma controls how much the minimizers of a function can change if another function is added. This will directly be useful in lower bounding the modulus of continuity. h be the largest minimizer of h and x ⋆ g be the smallest minimizer of g, and assume that x ⋆ h ≤ x ⋆ g . Let x ⋆ be any minimizer of h + g. Assume that h(x ⋆ h ) = 0 and g(x ⋆ g ) = 0. If h has λ h -quadratic growth and g is H g -smooth, then


text:
High order asymptotic preserving finite difference WENO schemes with constrained transport for MHD equations in all sonic Mach numbers

In this paper, a high-order semi-implicit (SI) asymptotic preserving (AP) and divergence-free finite difference weighted essentially nonoscillatory (WENO) scheme is proposed for magnetohydrodynamic (MHD) equations. We consider the sonic Mach number $\varepsilon$ ranging from $0$ to $\mathcal{O}(1)$. High-order accuracy in time is obtained by SI implicit-explicit Runge-Kutta (IMEX-RK) time discretization. High-order accuracy in space is achieved by finite difference WENO schemes with characteristic-wise reconstructions. A constrained transport method is applied to maintain a discrete divergence-free condition. We formally prove that the scheme is AP. Asymptotic accuracy (AA) in the incompressible MHD limit is obtained if the implicit part of the SI IMEX-RK scheme is stiffly accurate. Numerical experiments are provided to validate the AP, AA, and divergence-free properties of our proposed approach. Besides, the scheme can well capture discontinuities such as shocks in an essentially non-oscillatory fashion in the compressible regime, while it is also a good incompressible solver with uniform large-time step conditions in the low sonic Mach limit.

For MHD equations, divergence free of the magnetic field is a very important property. If the initial magnetic field is divergence-free, this condition will be maintained for all later times. Numerically, if this divergence-free condition is violated seriously, an unphysical force will be created to parallel the magnetic field [16], which can produce some numerical instabilities or nonphysical features in simulation. The constrained transport (CT) methodology, which was originally introduced by Evans and Hawley [26], is an effective way to maintain the (discrete) divergence-free condition up to machine precision. Many improvements and extensions [1, 4, 15-18, 20, 21, 29] are followed.
Recently, the MHD system in the low sonic Mach limit has attracted a lot of interests [19,25,33,34,42,45,48]. Jiang, Ju, and Li [33] have shown that a weak solution of compressible MHD equations will converge to a strong solution of their corresponding incompressible MHD equations in this low sonic Mach limit. They have also investigated the low Mach limit of full MHD equations with a heat conductivity [34]. Cui, Ou, and Ren [19] have established a uniform convergence from full compressible MHD equations to isentropic incompressible MHD equations with well-prepared initial conditions in the three-dimensional case. Numerically it is also very attractive to develop schemes for MHD flows at any speed. For the above mentioned explicit shock capturing schemes, when applied to the MHD system in the low sonic Mach regime, the time step is subject to a very strict CFL condition which is proportional to the sonic Mach number, making it very undesirable for all speed flows. An implicit time discretization can release the small time step restrictions from an explicit one, but usually it results in a highly nonlinear system, which is not easy to be solved efficiently with an iterative method and it is hard to guarantee a fast convergence. Furthermore, numerical viscosities for such schemes are inversely proportional to the sonic Mach number, which introduce excessive numerical dissipations. Instead, in recent years some researchers are in pursuit of semi-implicit (SI) schemes with low dissipations. In [25], Dumbser et. al. developed a pressure-based SI finite volume all Mach number flow solver for compressible MHD equations. Minoshima and Miyoshi [48] proposed a multistate low-dissipation advection upstream splitting method for MHD with both high and low Mach numbers, by using a Harten-Lax-van Leer discontinuities (HLLD) approximate Riemann solver. Leidi et al. [42] presented an SI finite volume solver with a 5-wave HLLD and a well-balanced method to efficiently simulate MHD flows at low Mach numbers with a gravitational source. For the MHD system with all sonic Mach numbers, it is important to design schemes with asymptotic stability and consistency as the sonic Mach number goes to 0 in the incompressible limit, namely, asymptotic preserving (AP). AP schemes have been widely used in a variety of areas, such as hydrodynamic or diffusive limits of kinetic models, relaxation limits of hyperbolic models, and low Mach number limits of compressible fluid models. We refer readers to the papers [23,24,27,[35][36][37][38] for more information. To achieve high order AP schemes in time, SI implicit-explicit Runge-Kutta (IMEX-RK) methods are widely used, which apply an explicit time discretization for nonstiff terms, while an implicit time discretization for stiff terms [5-11, 13, 39, 40, 47, 50]. With a suitable choice of explicit and implicit discretizations, both high efficiency and uniform stability independent of the stiff parameter can be obtained. As related, in [57] Tavelli and Dumbser developed an SI space-time discontinuous Galerkin method for all Mach-number flow, which involves staggered meshes and a Picard iteration for solving a large nonlinear system. In [8,10], Boscarino, Qiu, Russo, and Xiong presented high-order SI IMEX weighted essentially nonoscillatory (WENO) schemes for all-Mach isentropic Euler equations and all-Mach full Euler system, respectively. In [11], Boscheri and Pareschi developed high order pressure-based semi-implicit IMEX schemes for the three dimensional Navier-Stokes equations at all Mach numbers. In [30], Huang, Xing, and Xiong designed a high-order well-balanced SI AP scheme for shallow water equations in all Fraude numbers. We refer to many other related works from the references therein.
In this paper, we would like to develop a high-order AP finite difference WENO scheme with CT for MHD equations in all sonic Mach numbers. To our best knowledge, there are very few AP schemes for MHD systems in literature. The new contributions and innovations of this work are outlined as follows: • In the low sonic Mach number regime, SI IMEX-RK methods are adopted for time discretization to obtain a uniform time step independent of the sonic Mach number. Besides, to avoid a nonlinearity from the equation of state (EOS), an SI approach is used, leading to a linearized elliptic equation for the pressure. With carefully designed explicit and implicit discretizations, only a linear system needs to be solved. The resulting scheme is more efficient than explicit shock capturing schemes, especially in the low sonic Mach regime.
• Numerical viscosities are carefully designed to avoid excessive numerical dissipations, and a CT method is applied to maintain a discrete divergence-free property. The scheme can well capture discontinuities such as shocks in an essentially non-oscillatory fashion in the compressible regime. Meanwhile, the scheme is also a good incompressible solver as ε → 0. Due to less dissipation, our scheme can be shown to perform better than explicit schemes for low sonic Mach problems. • Formal AP and asymptotically accurate (AA) properties in the stiff limit as the sonic Mach number ε → 0 are proved, by assuming the implicit part of an SI IMEX-RK scheme is stiffly accurate (SA). The rest of the paper is as follows. In Section 2, we will briefly review the ideal MHD system, corresponding to its low sonic Mach limit based on asymptotic expansions, and a CT methodology for the divergence-free condition. A first-order SI scheme in time will first be introduced in Section 3, and then a high-order SI scheme with IMEX-RK in time and finite difference WENO in space is followed. Formal proofs of AP and AA properties are given in Section 4. Numerical experiments are performed in Section 5. A brief conclusion will be drawn in the last section.

Equations of compressible ideal MHD
The ideal MHD equations in a conservative form can be written as where ρ is the fluid mass density, u = (u, v, w) is the velocity, B = (B x , B y , B z ) is the magnetic field, E is the total energy density, p is the gas pressure, and · is the Euclidean vector norm. The total energy is given by where γ is the specific heat ratio. The system is subject to a divergence-free condition of the magnetic field If initially, the divergence-free condition holds, from (2.1), it also holds for all later times [18]. For this reason, (2.3) is usually not regarded as a constraint (like the ∇·u = 0 constraint for the incompressible Navier-Stokes equations), but rather an involution [16,18,29].
2.1. Eigenstructure of the MHD system. For (2.1), the eigenvalues of the Jacobian for the flux function along a normal direction n are given as: where λ 1,8 , λ 2,7 , λ 3,6 , λ 4 , and λ 5 are the left/right fast magnetosonic waves, left/right Alfvén waves, left/right slow magnetosonic waves, entropy wave, and divergence wave, correspondingly. The sound speed a, Alfvén speed c a , fast magnetosonic speed c f , and slow magnetosonic speed c s are defined as A fixed ordering for the MHD eigenvalues is Spectral decomposition for this system has already been described in [12,14,32,51,52,55], and we will follow the one in [32] for all simulations in the current paper.
2.2. Non-dimensionalized MHD. To derive a low sonic Mach limit for the ideal MHD system, we start with rewriting (2.1), (2.2), and (2.3) into a dimensionless form. We choose some reference or characteristic values, such as a length x 0 , a time t 0 , a fluid mass density ρ 0 , a velocity u 0 , a magnetic field B 0 , and a gas pressure p 0 , with u 0 = x 0 /t 0 . For the sake of simplicity, we consider a regime where the characteristic Alfvén Mach number u 0 √ ρ 0 /B 0 = 1. Then we define the following dimensionless Inserting these new variables into (2.1), we obtain the following non-dimensionalised MHD equations [42]: where we have dropped the hats of the dimensionless variables for ease of presentation. The EOS, divergence-free condition and the characteristic sonic Mach number ε are now given by In this dimensionless MHD system (2.6) with (2.7), the expressions of those eigenvalues (2.4) remain unchanged, except the sound speed a, Alfvén speed c a , fast magnetosonic speed c f , and slow magnetosonic speed c s from (2.5) now become Following the eigenstructure of the dimensionless MHD system discussed above, it is natural to find that if the background flow velocity u 0 is slow, then the sonic Mach number ε = u 0 √ ρ 0 / √ p 0 1, and the fast magnetosonic speed c f becomes very large. We can see that large magnetosonic speed c f will make the fast magnetosonic waves λ 1,8 to be much faster as compared to other waves. In this situation, for explicit shock-capturing schemes, the time step subjects to ∆t = CFL ∆x max where ∆t is the time step size, ∆x is the mesh size, and CFL is the time stability CFL number. On one hand, highly inaccurate solutions will appear due to the excessive numerical viscosity (scales as ε −1 ) in standard upwind schemes. On the other hand, this restriction results in an increasingly large computational cost for low sonic Mach fluid flows. A straight approach is to apply an implicit time discretization to avoid such severe time step constraints. However, a fully implicit scheme will result in a complicated nonlinear system which is usually very hard to be solved efficiently. Besides, it may not be able to guarantee a correct asymptotic limit. Thus, it is important to design efficient schemes with both asymptotic stability and consistency in the low sonic Mach limit, namely, the AP property.
2.3. Low sonic Mach limit. Next, we recall a formal derivation of the incompressible MHD system from the non-dimensionalized compressible MHD equations (2.6) with (2.7). We consider an asymptotic expansion ansatz for the following variables [46]: Inserting (2.8) into (2.6) and (2.7), equating to zero for different orders of ε, corresponding to the first three leading orders, we have where we have assumed a similar expansion of E as p in (2.8), and take E 0 = p 0 /(γ − 1). Denoting the material derivative d/dt = ∂/∂t + u 0 · ∇, from the energy equation in (2.9), E 0 = p 0 /(γ − 1), with ∇p 0 = 0, we have Taking into account periodic or no-slip u · n = 0 boundary conditions on the spatial domain Ω, and integrating (2.10) over Ω, we get Ω ∇ · u 0 dx = ∂Ω u 0 · nds = 0, where n is the unit outward normal vector along ∂Ω. This implies p 0 is constant not only in space but also in time, that is p 0 = Const. A direct result ∇ · u 0 = 0 can be obtained through (2.10) with p 0 = Const. Inserting ∇ · u 0 = 0 into the mass-continuity equation in (2.9), we get dρ 0 /dt = 0 which means that the density is constant along the characteristic trajectories. For the order O(ε) in (2.6), we have: Due to E 1 = p 1 /(γ − 1) from the EOS in (2.7), similarly as above, we get p 1 = Const and ∇ · u 1 = 0. In this case, we include εu 1 into u 0 , and εp 1 into p 0 , the asymptotic expansions of u and p become (2.11) u = u 0 + ε 2 u 2 + · · · , and (2.12) p = p 0 + ε 2 p 2 + · · · .
Using such expansions for u and p, as ε → 0, the leading order of the system (2.6) yields Numerically, in order to have a correct asymptotic limit from (2.6) to (2.13), it is important to impose well-prepared initial conditions, which are consistent with the expansions given in (2.8), with u and p from (2.11) and (2.12). According to [22,41], the well-prepared initial conditions for (2.6) are given as where ρ 0 (x) is a strictly positive function, p 0 = Const, ∇ · u 0 = 0, and ∇ · B 0 = 0.
2.4. CT for ∇ · B = 0. Many CT methods have been proposed to satisfy the discrete divergence-free condition ∇ · B = 0 for the ideal MHD equations [15][16][17][18]. Here we consider the one described in [18], in which the magnetic field is written as the curl of a magnetic vector potential A: , the equation of the magnetic field B in the MHD system can be rewritten as Inserting (2.15) into (2.16), we further get which implies the existence of a scalar function ψ such that Stable solutions can be obtained by introducing a Weyl gauge, i.e. ψ ≡ 0 [29], and the evolution equation for the vector potential becomes It is easy to check that ∇ · B = ∇ · (∇ × A) = 0, i.e. the divergence-free condition will be satisfied for all times. In summary, we get the following MHD system: Corresponding to (2.17), with well-prepared initial conditions (2.14), the low sonic Mach limit incompressible MHD system becomes

Numerical schemes
In this section, we will construct and analyze a high-order finite difference scheme with AP and divergence-free properties for the non-dimensionalized MHD system (2.17). The features of our scheme include the following: we design an SI IMEX-RK time discretization so that the scheme is stable with a time step constraint independent of the sonic Mach number ε; the AP property is preserved in the zero sonic Mach number limit; the scheme can be implemented efficiently without solving a nonlinear system. We will adopt the high order finite difference WENO scheme in [8,10,30], and the CT method in [18]. The final scheme can well capture discontinuities such as shocks in the compressible regime, while keeping a discrete divergence-free condition in the low sonic Mach regime, with a uniform time step condition independent of ε. We start with a first-order SI scheme. The high-order IMEX-RK scheme in time and a finite difference WENO reconstruction in space will be followed. After that, an overall scheme will be summarized.

3.2.
High order SI IMEX-RK scheme. In this subsection, we will follow a similar procedure as in [8,10] to extend the first-order SI IMEX scheme (3.1) to high order. For simplicity, we only consider (2.6) with CT methods, namely, (2.17). First, (2.17) can be written as an autonomous system: where U = (ρ, q, A, E) T and H : R n × R n × R m → R n is a sufficiently regular mapping. We use subscript "E" to express an explicit treatment for the first argument U and the third one B, and subscript "I" to express an implicit treatment for the second U . For the system (2.17), the function To obtain high-order accuracy in time, we can apply a multi-stage IMEX-RK time discretization with a double Butcher tableau,cÃ is an s × s matrix for an explicit scheme, withã ij = 0 for j ≥ i, and A = (a ij ) is an s × s matrix for an implicit scheme. For the implicit part of the methods, we use a diagonally implicit scheme, i.e. a ij = 0, for j > i, to guarantee simplicity and efficiency in solving the algebraic equations corresponding to the implicit part of the discretization. The vectorsc = ( Later, for the consideration of AP property, we consider implicit schemes with SA property, i.e., the implicit part of the Butcher table satisfies the condition b T = e T s A, with e s = (0, · · · 0, 1) T . We will see that the SA property guarantees that the numerical solution is identical to the last internal stage value of the scheme. Now, we may update the solutions as follows. Starting from U • Update the known values for the implicit part U and then solve The first order scheme (3.1) is the same as applying the following Butcher table to (2.17): namely U E = U n and U I = U n+1 .
3.3. High order finite difference WENO scheme. Below we will describe our spatial discretization strategies which incorporate a WENO mechanism in order to capture shocks in the compressible regime and produce a high order incompressible solver for the flow in the zero Mach limit. In particular, we propose to apply a fifth order characteristic-wise WENO procedure to F 1 E of the ideal MHD system so that in the compressible regime numerical oscillations could be well controlled, and to apply a component-wise WENO procedure for the flux functions of F 2 E and F SI . As for G E , which is a Hamiltonian, a WENO procedure described in [31] will be applied.
• Component-wise WENO ∇ W for ∇ · F 2 E and ∇ · F SI : Let v = (ρ, q, B, E) T and f = f (v) be the unknown variables and flux of the ideal MHD system (2.1). For simplicity, we consider the 1D case as an example, and we take a uniform mesh with N + 1 grid points: a = x 0 < x 1 < · · · < x N = b. A conservative finite difference WENO scheme for system (2.1) can be written in the following flux-difference form: are the numerical fluxes, which can be reconstructed by the following procedures: (1) Compute the physical flux at each grid point: (2) Perform a Lax-Friedrichs flux vector splitting for each component of the physical variables to get (3) Perform a finite difference WENO reconstruction to obtain upwind and downwind numerical fluxes: Here Φ WENO5 denotes a 5th-order WENO reconstruction [31], while other orders can also be used. Then we setf • Characteristic-wise WENO ∇ CW for ∇ · F 1 E : The Jacobian matrix for the flux function ∂f /∂v, has a spectral decomposition of the form ∂f /∂v = RDL [18], where D is a diagonal matrix of real eigenvalues, R is the matrix of right eigenvectors, and L = R −1 is the matrix of left eigenvectors. Then, for a characteristic-wise WENO reconstruction ∇ CW , we have the following procedures: (1) Compute the physical flux at each grid point: (2) Compute the average state v i+ 1 2 by using an arithmetic mean: (3) Compute the right and left eigenvectors of the flux Jacobian matrix ∂f /∂v by taking ).
(4) Project the solution and physical flux into the characteristic space: for all j in the numerical stencil associated with x = x i+ 1 2 . In the case of a 5th-order finite difference WENO scheme: j = i − 2, i − 1, i, i + 1, i + 2, i + 3.
where α has the same expression in the component-wise WENO (3.11). (6) Perform a WENO reconstruction Φ WENO5 on each component of g ± j to obtain their corresponding numerical fluxes, and then set (7) Project the numerical fluxĝ i+ 1 2 back to the physical space: However, when the sonic Mach number ε is small, the fast eigenvalue c f will make such a spectral characteristic decomposition to be unstable. Since a characteristic decomposition is mainly required in the compressible regime where ε = O(1), while it does not make any significant differences for small Mach numbers; in this case we directly use the eigenvectors in the expression of (2.1) by taking ε = 1, which works well from our numerical results.
• A Hamiltonian WENO reconstruction WENO-HJ (·) HJ for G E : We take a 1D Hamilton-Jacobi equation as an example where HT is the Hamiltonian. As described in [31], a semi-discrete form for (3.13) is where HT is a numerical Hamiltonian which approximates (3.13), and v − x,i , v + x,i are upwind and downwind approximations of v x at x = x i . A Lax-Friedrichs Hamiltonian from [31] is given as follows (3.14) HT Such a Hamiltonian WENO reconstruction helps us to control unphysical oscillations both in v x and v. For G E in 2D, the Hamiltonian is and the numerical Hamiltonian (3.14) can be applied along each direction for u(t, x, y) ∂Az ∂x and v(t, x, y) ∂A 3 ∂y , respectively. Here A z is the third component of A.

3.4.
Flowchart of a high order SI IMEX-RK finite difference WENO scheme. With the above space and time discretizations, we now summarize our high order scheme as follows: • Step 1. Starting from U n at time t n , for an intermediate stage (2) We can obtain ρ /2). Notice that in the process of substitution to obtain (3.16) or (3.18), the gradient and the divergence operators are kept to be continuous, obtaining the second order operator ∇ · (H (i) ∇p (i) I,2 ). The second-order spatial derivative is then discretized by a compact discretization as proposed in [10].
and successively update E Step 4. Finally, we update the numerical solution U n+1 = U (s) I , and B n+1 = ∇ × A n+1 . 4. Asymptotic preserving (AP) and Asymptotically Accurate (AA) properties 4.1. AP property. In this section, we will prove the AP property for the first order SI IMEX scheme (3.1), and the AA property for the high order SI IMEX-RK scheme (3.7)-(3.10). We focus on time discretizations while keeping the space continuous when we discuss the AP or AA property. First, we have the following theorem: Theorem 4.1. The first order SI IMEX scheme (3.1) with space continuous is AP, in the sense that, with well-prepared initial conditions (2.14), at the leading order asymptotic expansions, the scheme (3.1) is a consistent approximation of the incompressible MHD equations (2.18).

AA property.
Similarly to the AP property, here we focus on the AA analysis on time discretizations, while keeping space continuous. Then we have the following theorem: We consider a high order SI IMEX-RK scheme (3.7)-(3.10) of order r applied to system (2.17) on a bounded domain Ω ⊂ R with periodic or compact support boundary conditions, and assume that the implicit part of the SI IMEX-RK scheme is SA and the initial conditions (ρ 0 (x), ρ 0 (x)u 0 (x), p 0 , A 0 (x), B 0 (x)) T are well-prepared, namely in the form of (4.1). If we denote (ρ 1 (x; ε), ρ 1 (x; ε)u 1 (x; ε), p 1 (x; ε), A 1 (x; ε), B 1 (x; ε)) T as the numerical solution after one-time step, then, with p * being a constant, we have: ) T be the exact solution of the incompressible MHD equations (2.18) with the same initial data, one has the following one-step error estimate i.e. the scheme is AA.
Proof. We consider the first step from t 0 = 0 to t 1 = ∆t for the SI IMEX-RK scheme (3.7)-(3.10) of order r applied to system (2.17) with well-prepared initial data (4.1): where ρ 0 inc := ρ inc (x, 0), u 0 inc := u inc (x, 0), A 0 inc := A inc (x, 0), B 0 inc := B inc (x, 0). By well-prepared assumptions we have: p inc (x, 0) := p * is a constant independent of time and space, ∇ · u 0 inc = 0 and ∇ · B 0 (x) = 0. Now we consider a formal ε-expansion of the quantities: For example, the expansions of the density and pressure are To prove the theorem, we use mathematical induction.
• AA property for internal stages i = 1, · · · , s. When i = 1, it leads to the same AP analysis for the scheme (3.1) with ∆t replaced by a 11 ∆t. To prove the result for i > 1, we make use of the induction hypothesis, assuming the property holds for j ≤ i − 1, and prove that it still holds for j = i. For j = 1, · · · , i − 1, we have Now we insert the expansions (4.7) into the explicit step in (3.7), and from the energy equation we get with E 0 inc = p * /(γ − 1) and for j = 1, · · · , i − 1, , in the order of O(1), the density and momentum equations yield (4.12) where from (3.17) and (4.8), it follows ∇p (j) E = 0 for j = 1, · · · , i − 1. Using (4.10), we get (4.14) E Therefore, from (4.8) and (4.14), we get Now from (4.14) and (3.9), it follows for the energy equation

Numerical tests
In this section, we perform some numerical experiments to validate the high order accuracy, AP, AA, divergence-free, and good performances of our proposed scheme for the MHD system with all sonic Mach numbers. The third order SI IMEX-RK with an SA property from [8,30] is used for time discretizations. The fifth order finite difference WENO reconstruction [8,30,31,53] for the first order spatial derivatives, and a fourth order compact central difference discretization [8,10] for the second order derivatives will be applied for spatial discretizations. Our scheme is referred as the "IMEX" scheme. For some of the following examples, we will compare our results to reference solutions, which are produced by a fifth-order finite difference WENO scheme with a third-order strong stability preserving RK method. It is referred as an "ERK" scheme if without a CT method [32], while as an "ERKC" scheme if with a CT method [18]. To avoid excessive numerical dissipations for low sonic Mach numbers, for the explicit schemes "ERK" and "ERKC", here we also take a Lax-Friedrichs flux the same as in (3.11) and (3.12), and numerically we find that it leads to much better results as we can see in the following. Unless otherwise specified, the CFL number is taken as 0.25, and the gas constant is γ = 5/3. For 1D problems, the time step is with N + 1 computational grid points, and the definition ofĉ has been given in (3.12). For 2D problems, the time step is with (N x + 1) × (N y + 1) computational grid points. Here, the fast speeds in the x-direciton (ĉ (k,l) f,x ) and in the y-direciton (ĉ 5.1. One dimensional case. Since the divergence-free condition is satisfied automatically in 1D, therefore, the CT method will not need to be applied in the following 1D numerical simulations. with a periodic boundary condition. The exact solution to (5.1) propagates with a unit Alfvén speed (that is U (t, x) = U (0, x + t).). We set the time step to be ∆t = CFL∆x The sonic Mach number is ε = 1 and we compute up to a final time T = 1. The numerical L 1 , L 2 , and L ∞ errors and orders are shown in Table 5.1. We can observe that a fifth-order spatial accuracy is obtained.

5.2.
Two dimensional case. A CT method described in [18] will be used to maintain the divergencefree condition. In 2D, B = ∇ × A can be written as B x = ∂A z /∂y and B y = −∂A z /∂x. Therefore, we only need to update the third component A z of A, and ignore A x and A y . All solutions of the following tests maintain ∇ · B = 0 up to machine round off errors [16]. from Table 5.3 we can find around 5th-order accuracy for the momentum ρu, and 4th-order accuracy for the magnetic field component B x where a fourth-order central difference discretization is used for B = ∇ × A.
In this test, a loop of the magnetic field is advected. In Fig. 5.6, we show both the deviation of the density ρ away from 1 and the magnetic pressure B 2 /2 for the case of ε = 0.1 at T = √ 5. For the magnetic pressure, both schemes preserve the initial symmetric loop structure very well, and the IMEX scheme has almost the same results as the ERKC scheme. The one-dimensional cutting plots for the magnetic pressure along x = 0.2 show that both schemes slightly distort the profile at y = 0.15. For the deviation of the density, we can find the IMEX scheme is better than the ERKC scheme for preserving a symmetric structure, as compared with the results of the ERKC scheme on a uniform mesh of 384 × 192 cells. For this problem in the intermediate regime with ε = 0.1, we also compare the CPU cost of two schemes in Table 5.5. We can observe that the ERKC scheme costs less when ε is of O(1). However, the IMEX scheme is much more efficient in the intermediate and low sonic Mach regimes. The computational cost for the ERKC scheme increases very rapidly, as the sonic Mach number ε becomes smaller. This example demonstrates the main advantage of an AP scheme for such multiscale problems.  Example 5.9. (Magnetized Kelvin-Helmholtz instability) Finally, we consider a simulation of the magnetized Kelvin-Helmholtz instability problem [8,42]. The initial conditions are  Periodic boundary conditions are used for the conserved quantities, while a first-order extrapolation is utilized when the magnetic potential is outside of the computational domain. We run the solution up to T = 0.8 with γ = 1.4 and ε = 10 −6 on a mesh grid N x × N y = 256 × 128. We define the local sonic Mach number as M a = √ u 2 + v 2 / γp/ρ, and the ratio of the local sonic Mach number M ratio to the maximum value of M a as M ratio = M a / max(M a ), where the maximum is taken over all computational grid points. The vorticity ω = v x − u y , where v x and u y are discretized by the 4th order central difference, and M ratio are shown in Fig. 5.7. We observe that they are comparable to the results in [8,42], namely, our IMEX scheme can capture the incompressible MHD system in the low sonic Mach number regime very well, which verifies the AP property of our scheme.

Conclusion
In this paper, a high-order SI AP scheme with a discrete divergence-free property for the MHD equations has been developed. High-order accuracy in time has been obtained by SI temporal integrator based on an IMEX-RK framework. High-order accuracy in space has been achieved by finite difference WENO reconstructions. We have formally proved that the scheme is AP. We have also proved the AA property in the stiff limit as the sonic Mach number ε → 0, if the SI IMEX-RK scheme is SA. Numerical experiments in 1D and 2D have demonstrated the divergence-free for all ranges of sonic Mach numbers, AP and AA properties in the low sonic Mach limit. Compared to the explicit schemes as reference solutions, such as the ERK [32] and ERKC [18] schemes in 1D and 2D respectively, our IMEX scheme has almost the same performances as the ERK or ERKC scheme, but is much more efficient in the intermediate and low sonic Mach regimes.


An Empirical Evaluation of Convolutional Networks for Malaria Diagnosis

Malaria is a globally widespread disease caused by parasitic protozoa transmitted to humans by infected female mosquitoes of Anopheles. It is caused in humans only by the parasite Plasmodium, further classified into four different species. Identifying malaria parasites is possible by analysing digital microscopic blood smears, which is tedious, time-consuming and error prone. So, automation of the process has assumed great importance as it helps the laborious manual process of review and diagnosis. This work focuses on deep learning-based models, by comparing off-the-shelf architectures for classifying healthy and parasite-affected cells, by investigating the four-class classification on the Plasmodium falciparum stages of life and, finally, by evaluating the robustness of the models with cross-dataset experiments on two different datasets. The main contributions to the research in this field can be resumed as follows: (i) comparing off-the-shelf architectures in the task of classifying healthy and parasite-affected cells, (ii) investigating the four-class classification on the P. falciparum stages of life and (iii) evaluating the robustness of the models with cross-dataset experiments. Eleven well-known convolutional neural networks on two public datasets have been exploited. The results show that the networks have great accuracy in binary classification, even though they lack few samples per class. Moreover, the cross-dataset experiments exhibit the need for some further regulations. In particular, ResNet-18 achieved up to 97.68% accuracy in the binary classification, while DenseNet-201 reached 99.40% accuracy on the multiclass classification. The cross-dataset experiments exhibit the limitations of deep learning approaches in such a scenario, even though combining the two datasets permitted DenseNet-201 to reach 97.45% accuracy. Naturally, this needs further investigation to improve the robustness. In general, DenseNet-201 seems to offer the most stable and robust performance, offering as a crucial candidate to further developments and modifications. Moreover, the mobile-oriented architectures showed promising and satisfactory performance in the classification of malaria parasites. The obtained results enable extensive improvements, specifically oriented to the application of object detectors for type and stage of life recognition, even in mobile environments.

Introduction
Malaria is a globally widespread disease caused by parasitic protozoa transmitted to humans by infected female mosquitoes of Anopheles. In 2019, there were an estimated 229 million malaria cases worldwide, with an estimated 409,000 deaths due to malaria. Of them, 94% of malaria cases and deaths occurred in Africa [1]. In this context, children under five years of age are the most vulnerable group accounting for 67% (274,000) of all malaria deaths worldwide. Parasites of the genus Plasmodium (P.) cause malaria in humans by attacking red blood cells (RBCs). They spread to people through the bites of infected female Anopheles mosquitoes, called "malaria vectors". Five species of parasites cause malaria in humans: P. falciparum, P. vivax, P. ovale, P. malariae and P. knowlesi. P. falciparum and P. vivax are the two posing the most significant threat [1,2]. The former is most prevalent in Africa, while P. vivax is predominant in the Americas. Malaria plasmids within the human host have the following life stages: ring, trophozoite, schizont and gametocyte. The World Health Organization (WHO) defines human malaria as a preventable and treatable disease if diagnosed promptly. Still, the diagnosis must be made promptly, as the worsening illness can lead to disseminated intravascular thrombosis, tissue necrosis and spleen hypertrophy [1, [3][4][5].
Blood cell analysis using peripheral blood slides under a light microscope is considered the gold standard for the detection of leukaemia [6][7][8][9], blood cell counting [10][11][12][13][14] or the diagnosis of malaria [15][16][17]. Manual microscopic examination of peripheral blood smears (PBS) for malaria diagnosis has advantages such as high sensitivity and specificity compared to other methods. However, it requires about 15 minutes for microscopic examination of a single blood sample [18], and the quality of the diagnosis depends solely on the experience and knowledge of the microscopist. It is common for the microscopist to work in isolation without a rigorous system to ensure the quality of the diagnosis. In addition, the images analysed may be subject to variations in illumination and staining that can affect the results. In general, the manual process is tedious and time-consuming, and decisions dictated by misdiagnosis lead to unnecessary use of drugs and exposure to their side-effects or severe disease progression [19,20].
This work investigates the classification of malaria parasites using transfer learning (TL) to distinguish healthy and parasite-affected cells and classify the four P. falciparum stages of life. Moreover, the robustness of the models has been evaluated with cross-dataset experiments on two very different public datasets.
In this paper, transfer learning will be introduced by explaining how it works and discussing the pretrained networks selected to perform the comparative tests. The experiments are divided into (i) binary, (ii) multiclass and (iii) cross-domain classification. In the latter, networks trained on datasets from different domains were used to see if this improves accuracy over results obtained in a single domain.
The rest of the manuscript is organised as follows. Section 2 presents the literature on computer-aided diagnostic (CAD) systems for malaria analysis. Section 3 illustrates the datasets, methods and experimental setup. The results are presented and discussed in Section 4 and, finally, in Section 5, the findings and directions for future works are drawn.

Related Work
Several solutions for the automatic detection of malaria parasites have been developed in recent years. They aim to reduce the problems of manual analysis depicted in Section 1 and provide a more robust and standardised interpretation of blood samples while reducing the costs of diagnosis [15,21,22], mainly represented by CAD systems. They can be based on the combination of image processing and traditional machine learning techniques [23][24][25], and also deep learning approaches [16,[26][27][28], especially after the proposal of AlexNet's convolutional neural network (CNN) [29].
Since malaria parasites always affect the RBCs, any automatic malaria detection needs to analyse the erythrocytes to discover if they are infected or not by the parasite and, further, to find the stage of life or the type.
Among the more recent and classical solutions not employing CNNs, Somasekar et al. [23] and Rode et al. [25] proposed two malaria parasite segmentation methods. The first one used fuzzy clustering and connected component labelling followed by minimum perimeter polygon to segment parasite-infected erythrocytes and detect malaria, while the second one is based on image filtering and saturation separation followed by triangles thresholding.
Regarding the CNN-based approaches, Liang et al. [26] proposed a novel model for the classification of single cells as infected or uninfected, while Rajaraman et al. [27] studied the accuracy of CNN models, starting from pretrained networks, and proposed a novel architecture trained on a dataset available from the National Institutes of Health (NIH). They found that some pre-existing networks, by means of TL, can be more efficient than networks designed ad hoc. In particular, ResNet-50 obtained the best performance. Subsequently, they further improved through an ensemble of CNNs [28]. Rahman et al. [30] also exploited TL strategies using both natural and medical images and performed an extensive test of some off-the-shelf CNNs to realise a binary classification.
Some other techniques not explored in this work are based on the combination of CNNextracted features and handcrafted ones [31][32][33] or the direct use of object detectors [34]. For example, Kudisthalert et al. [33] proposed a malaria parasite detection system, based on the combination of handcrafted and deep features, extracted from pretrained AlexNet. Abdurahman et al. [34] realised a modified version of the YOLOV4 detector. Moreover, they generated new anchor box sizes with a K-means clustering algorithm to exploit the model on small objects.
Finally, a recent focus has been posed on mobile devices, which enable a cheaper and quicker diagnosis in the underdeveloped areas of the world, where more expensive laboratories do not exist. As an example, Bias et al. [24] realised an edge detection technique based on a novel histogram-based analysis, coupled with easily accessible hardware, focused on malaria-infected thin smear images.
The work in [30] is the most similar to the approach here proposed. In particular, they compared different off-the-shelf networks for a binary classification using two datasets, one is the Malaria Parasite Image Database for Image Processing and Analysis (MP-IDB) [35] and another is composed of synthetic and medical images. The task faced, however, is a binary classification. On the entire MP-IDB, they reported 85.18% accuracy with a fine-tuned version of VGG-19.
In summary, the main difference between our work and the state-of-art is that here an extended set of off-the-shelf CNNs on two very different public datasets have been exploited with a dual purpose: detect healthy and unhealthy RBCs and distinguish the various stages of life. Finally, it is the first baseline provided for the stages of life classification on the MP-IDB.

Materials and Methods
In this section, the datasets, the techniques and the employed experimental setup are described.

Datasets
Two well-known benchmark datasets were used: the National Institutes of Health (NIH) [27], proposed for malaria detection, and MP-IDB [35], a recently proposed dataset for malaria parasite types and stages of life classification.

NIH
The NIH is a public PBS images dataset from healthy individuals and malaria-affected patients. Image acquisition was performed at the Lister Hill National Center for Biomedical Communications from Giemsa-stained blood samples obtained from 150 patients infected with P. falciparum and 50 healthy patients. The resulting dataset consists of 27,558 images, uniformly subdivided between infected and healthy cells. The infected cells contain a ring-shaped parasite. Figure 1 shows a healthy and sick RBC extracted from NIH.

MP-IDB
The MP-IDB consists of four malaria parasite species, P. falciparum, P. malariae, P. ovale and P. vivax, represented by 122, 37, 29 and 46 images, respectively, for a total amount of 229. The images have been acquired with 2592 × 1944 resolution and 24-bit colour depth. Moreover, every species contains four distinct life stages: ring, trophozoite, schizont and gametocyte. Figure 2 shows four examples of the types of malaria parasites included in MP-IDB.

Deep Learning
In this work, deep learning approaches have been used as classifiers. In particular, eleven different off-the-shelf CNN architectures have been evaluated. All of them have been pretrained on the well-known natural image dataset ImageNet [36], and then adapted to this medical image task, as proposed in [37], following a transfer learning and finetuning procedure.
AlexNet [29] and VGG-16 [38] are pretty simple and similar architectures, composed of 8 and 16 layers, respectively. Nevertheless, they are widely used for transfer learning and fine-tuning [37], as they have gained popularity for their excellent performance in many classification tasks [29]. GoogLeNet [39] and Inceptionv3 [40] are both based on the inception layer; in fact, Inceptionv3 is a variant of GoogLeNet, using 140 levels, 40 more than GoogLeNet. The 3 ResNet architectures have 18, 50, 101 layers for ResNet-18, ResNet-50 and ResNet-101, respectively, based on residual learning. They are easier to optimise even when the depth increases considerably [41]. The building block of ResNet inspired DenseNet-201. It deals with the vanishing-gradient problem by introducing a connectivity pattern between layers. It is comprised of multiple densely connected layers, and their outputs are connected to all the successors in a dense block [42]. ShuffleNet [43], SqueezeNet [44] and MobileNetV2 [45] are lighter networks. In particular, the last two are oriented to real-time executions and mobile device usage.
Regarding the transfer learning strategy, the approach used in [37] was followed. All CNN layers were retained except for the last fully connected one. It was replaced with a new layer, initialised and set up to accommodate new categories according to the classification strategy exploited (two classes in the binary one and four in the multiclass one).

Image Preprocessing
All the images have a uniform background, even if different according to the dataset. Indeed, NIH contains only single-cell images surrounded by a black background, while the images of MP-IDB have several blood components with their actual plasma background. Both datasets are organised into classes, healthy or sick, for NIH and the four malaria types for MP-IDB. As it can be seen from Figures 1 and 2, there are many colour variations between the images of both datasets. This condition is undoubtedly due to the different acquisition conditions and from the status of the blood smear [35]. Therefore, it was considered appropriate to apply a preprocessing step to realise a colour balance. Furthermore, for MP-IDB, we created a single image for each parasite from the full-size images. The preprocessing step was designed mainly to adjust the colour components of the images and applied to all the RGB channels, using a colour-balancing technique, through Equation (1), where C out is the processed component, C in is the component to be processed, m im is the average of the average intensities of the three channels and, finally, m c is the average intensity of the component to be processed. This procedure was carried out on all three channels of the RGB image.
As far as MP-IDB is concerned, the work was also oriented to generate single-cell images from the full-size and, more specifically, the Falciparum class. From now on, we refer to this dataset as MP-IDB-Falciparum-Crops (MP-IDB-FC). MP-IDB contains 104 patientlevel images with the corresponding ground truths of the infected cells. As the desired classification was at the cellular level, individual cells had to be extracted from the images and labelled using the class references in the file names. The Falciparum class presents the following image distribution per stage: 1230 rings, 18 schizonts, 42 trophozoites and 7 gametocytes. In Figures 3 and 4, the preprocessing step results are shown.
Finally, to overcome the class imbalance issue in the dataset produced, a further augmentation step was applied and described below.
Sick RBC before.
Sick RBC after. Healthy RBC before. Healthy RBC after.

Data Augmentation
Different data augmentation techniques, such as flipping, shifting and rotating, are used to overcome the problem of a limited dataset and to produce further copies of the original image to give the algorithm more generalisation capability and reduce the error rate [29,[46][47][48].
MP-IDB-FC presented an unbalanced distribution of images per class; therefore, we proposed an offline data augmentation to oversample the underrepresented classes. The applied geometric transformations are random rotations between −90 • and 90 • , random translations between −10 and 10 pixels on the X-axis and Y-axis, and 50% reflection around the X-axis and Y-axis.

Experimental Setup
As previously introduced, the images to be classified are microscopic blood smear images. More specifically, NIH contains single-cell images representing healthy or ring stage malaria parasite-affected RBC, while in the case of MP-IDB-FC, the images represent a single cell containing one of the four stages of P. falciparum.
In the case of NIH, the dataset was divided into three parts, one for the training (80%), one for the validation (10%) and 10% for testing. A stratified sampling procedure to keep the splits balanced was used, as NIH classes are well balanced.
As for MP-IDB, in order to consider the class imbalance and to have a sufficient number of samples for the training process while preserving a sufficient number of samples for performance evaluation, the dataset was split first into two parts, namely training and testing set, with 80 and 20% of images, respectively. A total of 10% of the training set was used for validation due to the small number of images. Naturally, the subdivision of the sets was carried out before oversampling to avoid that train or validation images fell into the test set, thus compromising the results. The subdivision was made to ensure that 20% of each class made up the test set, thus avoiding the circumstance where no class images with few elements were in the test set. Oversampling and undersampling were then adopted to increase the number of images to 100, 200 or 300 per class in the training set. Moreover, the splits were not created randomly but by taking the images in lexicographic order from each class to further ease reproducibility. All the experiments have been conducted on a single machine with the following configuration: Intel(R) Core(TM) i9-8950HK @ 2.90 GHz CPU with 32 GB RAM and NVIDIA GTX1050 Ti 4GB GPU. Finally, no randomisation strategy has been applied to make all the experiments on both datasets reproducible.
This work represents a baseline for further investigation and searching for the best architecture for our purpose. The selected CNNs have been employed without any modification to their architecture. In particular, after empirical evaluation, Adam algorithm was adopted, which performed better than the other solvers. In addition, the maximum number of epochs was set to 10 due to the number of images.
As mentioned in Section 3.2.1, a fine-tuning process was applied on all the CNN architectures exploited on both datasets. The hyperparameters defined in Table 1 were used for all networks to evaluate potential performance variations. Furthermore, the regularisation factor L 2 was set to avoid a possible overfitting during the training phase.

Experimental Results
Three different experiments were conducted, according to the classification purpose: The results obtained in the analysis of each experiment were performed using the confusion matrix. The confusion matrix metric used in this study is Accuracy. The formula of this metric is given in Equation (2). The variables used in the equation are True Positive (TP), False Positive (FP), True Negative (TN), and False Negative (FN), parameters of the confusion matrix used to calculate the metrics [49,50].

Binary Classification Performance on NIH
To determine the training options to use, test trials were carried out. From them, it was mainly verified that: • Extending the training phase beyond ten epochs did not improve accuracy, as the network stored individual image features rather than class features, and overfitting compromised the results; • The ideal learning rate was 1 × 10 −4 . The accuracy increased too slowly for smaller values, and for larger ones, it did not converge to a specific value; • Empirically, Adam was found as the best solver. Table 2 shows that almost all the networks have an accuracy value close to the average. The standard deviation of the collected data is solely 0.16%. This aspect could be because the dataset used has many valuable images for training the network. In particular, ResNet-18 recorded the highest accuracy value, confirming the high performance expressed in [27]. MobileNetV2, SqueezeNet and ShuffleNet recorded average values, which is an important result since they are networks designed for mobile use.

Multiclass Classification Performance on MP-IDB-FC
The multiclass classification on MP-IDB-FC was designed to determine the life stage of the parasite: ring phase, adult trophozoite, schizont and gametocyte.
Like the binary classification, comparative tests with the same training set, validation set and test set were carried out to determine which networks perform best and allow comparison between them. We created three datasets with 100, 200 or 300 images per class in the training set. We refer to these sets as D1, D2 and D3, respectively. They were constructed by oversampling with augmentation of the remaining classes. Table 3 shows the performance in this experiment. Each test was cross-validated five times; then, we reported the mean accuracy and standard deviation considering each of the five folds. The most notable result is that the average performance from D2 to D3 worsens. However, DenseNet-201 and GoogLeNet are the only networks to benefit from increasing the dimensionality of the training set.

Cross-Dataset Classification Evaluation
The cross-domain classification was carried out to evaluate the CNNs robustness.

MP-IDB-FC Classification with NIH Models
Firstly, two different multiclass classifications were realised on MP-IDB-FC by: • Training on NIH and testing on MP-IDB-FC (Exp1); • Training on NIH + fine-tuning on MP-IDB and testing on MP-IDB-FC (Exp2).
When used for training or fine-tuning, the split on MP-IDB-FC was 50% for training (with 10% for validation) and 50% for testing. Every test was evaluated with five-fold crossvalidation. Therefore, we report the average accuracy of the five folds and the standard deviation. This cross-domain experiment tested whether the networks trained on the NIH dataset could be employed on the MP-IDB-FC dataset. It is helpful to point out that the significant difference is that, on the one hand, MP-IDB-FC has parasite crops and not healthy RBCs and, on the other hand, NIH contains only ring-stage parasites. For this reason, the objective of Exp1 was to discriminate between rings and the remaining stages of life, while Exp2 aimed to expand the knowledge of the NIH pretrained models with new information on the stages of life.
The results depicted in Table 4 show that when the target domain differs excessively from the source domain, it is hard to directly apply the models trained with NIH to MP-IDB-FC (Exp1), even if the task seemed feasible, as ring-stage parasites were contained in both datasets. Conversely, Exp2 shows that using the CNNs first trained on NIH and then fine-tuned on MP-IDB-FC led to an improvement in average accuracy. The information about healthy RBCs provided with NIH training does not affect the overall result. In addition, the standard deviation is under 4% for all the networks, leading to satisfactory performance stability.  Trophozoites, schizonts and gametocytes greatly vary between the two types, while the ring stages are pretty similar.
As it can be seen from Table 5, the classification of P. falciparum stages of life employing models trained P. vivax produced dreadfully low results due to the differences between all the stages except rings. On the other hand, using same-domain models (Exp4) had satisfactory results. Exp5 demonstrates that the fine-tuning strategy on the models pretrained on P. falciparum improved the accuracy, as already happened in Section 4.3.1 Exp2. In this task, DenseNet-201 provided the best performance, being the only CNN to overcome 85% and outperforming the average of 10% in both cases.

Conclusions
The results obtained in this work support the importance of deep learning in haematology. This work aimed to demonstrate that pretrained off-the-shelf networks can offer high accuracy for diagnosing malaria utilising transfer learning however showing several limitations of this approach. Several comparative tests were developed using a selection of pretrained networks differentiated by size, depth and the number of parameters. In particular, using the NIH dataset, it is possible to distinguish a healthy from an infected erythrocyte with an accuracy of over 97%. Small networks such as SqueezeNet and ShuffleNet performed well, consolidating a possible development of software for malaria diagnosis in small devices such as smartphones. On the other hand, MP-IDB has highlighted some critical issues: deep learning is not very effective when the dataset used for training is unbalanced. Some classes of parasites in the dataset have a small number of images. Nevertheless, the oversampling, augmentation and preprocessing methods still allowed us to exceed 90% accuracy on the test set for distinguishing the four life stages of the P. falciparum parasite. Finally, the cross-domain experiments have highlighted some critical points in classifying data from heterogeneous domains. It was counterproductive to apply the models trained with NIH to MP-IDB-FC, but the use of the CNNs firstly trained on NIH and fine-tuned on MP-IDB-FC led to an improvement in average accuracy. This aspect also applies to using the P. vivax dataset as the target domain, as most of the classes deviate too much from the corresponding P. falciparum classes. However, using both types of parasites as source domains produced better results than training on P. vivax only. In general, the extensive experimentation has highlighted how DenseNet-201 offers the most stable and robust performance, offering itself as a crucial candidate for further developments and modifications.
Among the possible developments of this work, we aim to propose a framework able to detect malaria parasites from blood smear images and classify different species of parasites and different stages of life, mainly focusing on high variation data. We also plan to use domain adaptation algorithms to improve cross-domain performance. Data Availability Statement: All the material used and developed for this work is available at the following GitHub repository. All the models trained for the different experimentations can be found at this repository.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: