extfieldsofstudy:
Mathematics
Computer Science

text:
Entropic regularization of Wasserstein distance between infinite-dimensional Gaussian measures and Gaussian processes

This work studies the entropic regularization formulation of the 2-Wasserstein distance on an infinite-dimensional Hilbert space, in particular for the Gaussian setting. We first present the Minimum Mutual Information property, namely the joint measures of two Gaussian measures on Hilbert space with the smallest mutual information are joint Gaussian measures. This is the infinite-dimensional generalization of the Maximum Entropy property of Gaussian densities on Euclidean space. We then give closed form formulas for the optimal entropic transport plan, entropic 2-Wasserstein distance, and Sinkhorn divergence between two Gaussian measures on a Hilbert space, along with the fixed point equations for the barycenter of a set of Gaussian measures. Our formulations fully exploit the regularization aspect of the entropic formulation and are valid both in singular and nonsingular settings. In the infinite-dimensional setting, both the entropic 2-Wasserstein distance and Sinkhorn divergence are Fr\'echet differentiable, in contrast to the exact 2-Wasserstein distance, which is not differentiable. Our Sinkhorn barycenter equation is new and always has a unique solution. In contrast, the finite-dimensional barycenter equation for the entropic 2-Wasserstein distance fails to generalize to the Hilbert space setting. In the setting of reproducing kernel Hilbert spaces (RKHS), our distance formulas are given explicitly in terms of the corresponding kernel Gram matrices, providing an interpolation between the kernel Maximum Mean Discrepancy (MMD) and the kernel 2-Wasserstein distance.

Introduction
In this work, we study the entropic regularization formulation of the 2-Wasserstein distance in the Hilbert space setting, with a particular focus on Gaussian measures and covariance operators on Hilbert space. This is the infinite-dimensional generalization of recent work on the entropic 2-Wasserstein distance between Gaussian measures on R n , as reported in [55,45,24]. Our work is along the direction of entropic regularization in optimal transport, which has recently attracted much attention in various fields, in particular machine learning and statistics [20,83,30,35,36,38,56,75,77], with applications in computer vision, density functional theory, and inverse problems (e.g. [36,37,53,69]). This direction of research is also closely connected with the Schrödinger bridge problem [81], which has been studied extensively [11,18,21,31,78,91,33,52,79,80].
Our focus in the Gaussian setting stems not only from its use in elucidating various aspects of the abstract theory, since many quantities of interest admit closed form formulas, but also from numerous applications utilizing Gaussian measures and covariance matrices/operators. These include brain imaging [3,27], computer vision [88,89,87], and brain computer interfaces [16]. Many distances/divergences have been studied and employed in practice, including the affine-invariant Riemannian metric [71], corresponding to the Fisher-Rao distance between centered Gaussians, the Alpha Log-Determinant divergences [13], corresponding to Rényi divergences between centered Gaussians, the Log-Euclidean metric [4], and recent work attempting to unify them [2,15,86,63].
Infinite-dimensional setting. The generalization of distances/divergences for Gaussian measures and covariance matrices on R n to the infinite-dimensional setting of Gaussian measures and covariance operators on Hilbert spaces has been carried out by various authors. In general, the infinite-dimensional formulations are substantially more complex than the finite-dimensional ones and regularization is often necessary. This is the case for the affine-invariant Riemannian distance [51], the Log-Hilbert-Schmidt metric [60], the Alpha and Alpha-Beta Log-Determinant divergences [59,61,58,64,65]. The settings for these distances/divergences are the sets of positive definite unitized trace class/ Hilbert-Schmidt operators, which are positive trace class/Hilbert-Schmidt operators plus a positive scalar multiple of the identity operator so that operations such as inversion, logarithm, and determinant, are well-defined. A particular advantage of the 2-Wasserstein distance compared to the above distances/divergences is that the finite and infinite-dimensional distance formulas [34,19] are the same and no regularization is necessary.
Reproducing kernel Hilbert space (RKHS) setting. From the computational and practical viewpoint, this setting is particularly interesting since many quantities of interest admit closed forms via kernel Gram matrices which can be efficiently computed. Examples include the kernel Maximum Mean Discrepancy (MMD) [41] and the RKHS covariance operators, the latter resulting in powerful nonlinear algorithms with substantial improvements over finite-dimensional covariance matrices, see e.g. [93,43,60,66,92] for examples of applications in computer vision.
Contributions of this work.
1. We generalize the Maximum Entropy property of Gaussian densities in R n to the Hilbert space setting, namely the Minimum Mutual Information of joint measures of two Gaussian measures on Hilbert space.
2. For two Gaussian measures on a Hilbert space H, we provide closed form formulas for the optimal entropic transport plan, the entropic 2-Wasserstein distance and the Sinkhorn divergence, generalizing results in [55,45,24]. 3. For a set of Gaussian measures, we show a new Sinkhorn barycenter equation, with always a unique non-trivial solution. In contrast, we show that the finite-dimensional barycenter equation for the entropic 2-Wasserstein distance fails to generalize to the Hilbert space setting. 4. In the RKHS setting, we present closed form formulas for the distances via the finite kernel Gram matrices, providing an interpolation between the kernel MMD [41] and kernel Wasserstein distance [92,62]. 5. Our proofs and results fully exploit the regularization aspect of the entropic formulation and are valid both in singular and non-singular settings. This is novel also in the finite-dimensional setting (compared to [55,45,24]). 6. As we discuss in detail, many properties of the Gaussian case have not been proved in the general theory, due to (i) dim(H) = ∞, (ii) c(x, y) = ||x−y|| 2 is unbounded on H, (iii) the support of Gaussian measures is unbounded .
The Banach space Tr(H) of trace class operators on H is defined by (see e.g. [76]) Some key differences between the finite and infinite-dimensional settings are 1. On R n , for two random variables X, Y with joint and marginal measures µ XY , µ X , µ Y , having densities f (x, y), f X (x), f Y (y), respectively, with respect to the Lebesgue measure, their mutual information is defined as [17] I(X; Y ) = where H(X) = − R n log[f X (x)]f X (x)dx is the differential entropy of X.
The classical Maximum Entropy of Gaussian densities property ( [43], Theorem 9.6.5) states that if X has mean zero and covariance matrix C, then H(X) ≤ 1 2 log(2πe) n det(C), with equality if and only if X ∼ N (0, C).
Thus if both X and Y have Gaussian densities, then I(X; Y ) is minimum if and only if their joint density is Gaussian, so that for c(x, y) = ||x − y|| 2 , and ν 0 , ν 1 being Gaussian, a minimizing γ in (4) is necessarily a joint Gaussian measure of ν 0 , ν 1 . When C is a covariance operator on H with dim(H) = ∞, the quantity det(C) is no longer well-defined. However, I(X; Y ) = KL(µ XY ||µ X ⊗ µ Y ) is well-defined and finite whenever µ XY is absolutely continuous with respect to µ X ⊗ µ Y . In the following, we show that the above Minimum Mutual Information property of joint Gaussian measures generalizes to the infinite-dimensional setting.

If
A is a strictly compact operator on H, then A −1 is unbounded when dim(H) = ∞. Thus finite-dimensional methods that utilize matrix inversion extensively, e.g. in [55,45,24], are not applicable when dim(H) = ∞. Instead, we fully exploit the regularization aspect of problem (4) and invert operators of the form γI + A > 0, thus our proofs fully resolve this issue and are valid in the general setting when A can be singular.
3. The identity operator I is not trace class when dim(H) = ∞. This leads to the breakdown in the entropic barycenter problem (Theorem 12) and has consequences for the analysis of the existence of solutions of the barycenter equations (detail given in Section 11).

Main Results
We first state the following generalization of the Maximum Entropy of Gaussian densities in R n . To the best of our knowledge, this property has not been explicitly and rigorously presented in the literature in the infinite-dimensional setting. It states that among all joint measures, with the same covariance operators, of two Gaussian measures µ X , µ Y on two separable Hilbert spaces H 1 , H 2 , the ones with the minimum Mutual Information are precisely the joint Gaussian measures on H 1 × H 2 . In the following, Gauss(H) denotes the set of all Gaussian measures on H and Gauss(µ X , µ Y ) denotes the set of joint Gaussian measures having marginals µ X and µ Y .
Theorem 1 (Minimum Mutual Information of Joint Gaussian Measures) Let H 1 , H 2 be two separable Hilbert spaces. Let µ X = N (m X , C X ) ∈ Gauss(H 1 ), µ Y = N (m Y , C Y ) ∈ Gauss(H 2 ), ker(C X ) = ker(C Y ) = {0}. Let γ ∈ Joint(µ X , µ Y ), γ 0 ∈ Gauss(µ X , µ Y ), γ 0 is equivalent to µ X ⊗ µ Y . Assume that γ and γ 0 have the same covariance operator Γ and that µ X ⊗ µ Y has covariance operator Γ 0 . Then Equality happens if and only if γ = γ 0 . Here V is the unique bounded linear operator satisfying V ∈ HS(H 2 , H 1 ), ||V || < 1, such that Γ = Γ 1/2 0 The operator V in Theorem 1 is defined in Section 7. It links the covariance operators C X and C Y with the cross-covariance operator C XY via the relation Y . In Eq. (9), det refers to the Fredholm determinant (see e.g. [82]). Let A ∈ Tr(H), then the Fredholm determinant of I + A is given by det(I + A) = ∞ j=1 (1 + λ j ), where {λ j } j∈N are the eigenvalues of A. Theorem 1 in turn follows from the following more general result on the KL divergence on Hilbert space. It states in particular that if µ is a Gaussian measure on H, then among all probability measures with the same mean and covariance operator, KL(γ||µ) is minimum if and only if γ is Gaussian.
In Theorem 2, det 2 refers to the Hilbert-Carleman determinant (see e.g. [82] A direct consequence of Theorem 1 is that if C 0 , C 1 are nonsingular, then a minimizer of problem (13) is necessarily a joint Gaussian measure of µ 0 and µ 1 . We show that this holds in the general setting, i.e. in both nonsingular and singular cases. The following result gives the explicit formula for this minimizer, which is unique for any ǫ > 0. It is proved in Section 8, using two different methods: (i) by directly solving the optimization (13), and (ii) by solving the corresponding Schödinger system.
For each fixed ǫ > 0, problem (13) has a unique minimizer γ ǫ , which is the Gaussian measure where The Radon-Nikodym derivative of γ ǫ with respect to µ 0 ⊗ µ 1 is given by where the functions α ǫ : R → R and β ǫ : R → R take the form The constants a, b ∈ R and the operators A, B : H → H are given by Here det is the Fredholm determinant and M ǫ ij : H → H, are defined by Remark 1 The operator M ǫ ij as defined in Eq. (20) can be rewritten as from which it follows that M ǫ ij ∈ Sym + (H) ∩ Tr(H). Thus for the operator the Fredholm determinant det I + 1 2 M ǫ ij is well-defined and positive.
We remark that in the case the cost function c(x, y) is bounded, much more can be said about the duality formulation, see [25].
is convex in each argument. In particular, let C 0 be fixed, then the function Theorem 6 (Lower bound for Entropic 2-Wasserstein distances between general probability measures on Hilbert space) Let µ 0 , µ 1 be two probability measures on H, with means m 0 , m 1 and covariance operators C 0 , C 1 , respectively. Then Equality happens if and only if µ 0 , µ 1 are Gaussian measures.
In [30], the Sinkhorn divergence was proved to be convex in each variable for either a compact metric space X or for measures with bounded support on R n , with cost function c(x, y) = ||x − y|| p , p = 1, 2. In [44], this was shown for sub-Gaussian measures on R n , c(x, y) = ||x − y|| 2 . The following shows strict convexity for Gaussian measures on H with c(x, y) = ||x − y|| 2 . This property is crucial for guaranteeing the uniqueness of the barycenter problem below.
In [30], positivity of the Sinkhorn divergence was proved for either a compact metric space X or for measures with bounded support on R n , with cost function c(x, y) = ||x − y|| p , p = 1, 2. We now show that in the Gaussian case, with c(x, y) = ||x − y|| 2 , this holds in the much more general Hilbert space setting.
Furthermore, F S is strictly convex and F E is convex on Ω, with strict convexity if C 0 is strictly positive.
The expressions for OT ǫ d 2 and S ǫ d 2 in Theorems 4 and 7 are not intuitively close to the exact OT formula, which is the same as in the finite-dimensional setting [34,19]. Theorem 11 below gives equivalent formulas that better express the connections between the exact and regularized settings. .
Given N probability measures µ i ∈ P(H), i = 1, 2, .., N , the entropic barycenterμ with weights w i > 0, The entropic barycenter problem illustrates clearly the bias effect of the entropic regularization term, as analyzed in the finite-dimensional case, e.g. [44,45]. This bias effect is sharp when dim(H) = ∞, with the finite-dimensional barycenter equation failing to generalize to this case. In the following, we call a barycenter trivial if it is a Dirac delta measure in H.

Theorem 12 (Entropic Barycenter of Gaussians)
.., N be a set of Gaussian measures on H. Then a minimizer of problem (38) in P 2 (H) is necessarily Gaussian. Assume at least one of the C ′ i s is strictly positive. Then on Gauss(H), problem (38) is strictly convex and the first order minimality condition is 1. If ǫI ≥ 2 N i=1 w i C i , the unique barycenter is the Dirac delta measure centered atm = N i=1 w i m i . 2. If dim(H) < ∞, a necessary condition for the existence of a non-trivial barycenter is A sufficient condition for the existence of a non-trivial barycenter is ∃α ∈ R, α > 0 such that C i ≥ αI, 1 ≤ i ≤ N, and 0 < ǫ < 2α. (41) In this case, the barycenter is unique and is the Gaussian measure N (m,C), whereC > 0 is the unique solution of Eq. (39). Equivalently,C is the unique strictly positive solution of the following equation 3. If dim(H) = ∞, then Eq. (39) has no solution in Sym + (H). If ǫ 2 N i=1 w i C i , then the barycenter, if it exists, is not the Dirac measure centered atm.
Discussion of results. Eq.(39) is the first order optimality condition for the strictly convex problem (38). It has sharply different behavior when dim(H) = ∞ compared with the case dim(H) < ∞, as we stated.
Conditions (40) and (41) show that, even in the finite-dimensional setting, a non-trivial entropic barycenter of Gaussian measures exists if and only if ǫ is sufficiently small. Condition (40), namely 0 < ǫI < 2 Under condition (41), Eq.(39) has a unique solution, which is strictly positive. Then Eq.(42) also has a unique strictly positive solution, but it also has the trivial solution X 0 = 0 and uncountably infinitely many positive solutions, which are singular (see Theorem 14 and Proposition 7).

Remark 3
In general, for ǫ > 0, the entropic barycenter problem (38) does not make much sense. Consider the following one-dimensional scenario, where Thus if ǫ ≥ 2σ 2 , then Equation (39) has no positive solution. The above solutionC is also obtained for the case N = 1, w 1 = 1, C 1 = σ 2 > 0, in which case a sensible solution should beC = C 1 , i.e. the barycenter of a set of one point should the point itself. This seemingly pathological behavior is not necessarily surprising, since OT ǫ d 2 is neither a distance nor a divergence.
Sinkhorn barycenter of Gaussian measures. We now consider the barycenter problem with respect to the Sinkhorn divergence. For a set of probability In the current work, we consider barycenter of the N Gaussian measures in the set of all Gaussian measures on H µ = arg min In contrast to the entropic barycenter problem, the debiased Sinkhorn barycenter problem has a consistent generalization to the infinite-dimensional setting, with a unique solution that is valid in both singular and nonsingular cases.
Furthermore,C is strictly positive if and only if Under the additional hypothesis thatC > 0,C is equivalently the unique strictly positive solution of the following equation Define the following map F : Sym + (H) → Sym + (H) by Then the unique solution of Eq.(45) is the unique fixed point of F .
Limiting cases. When ǫ → 0, both Eqs. (42) and (47) become In the finite-dimensional setting, this is the barycenter equation for the exact 2-Wasserstein distance [1], assuming thatC > 0. As of the current writing, to the best of our knowledge, a rigorous proof for the infinite-dimensional case has not yet been established. We note that the proof given in [54], which uses the transport map in [19] to compute gradients, is only applicable in the case dim(H) < ∞, since the transport map is generally unbounded when dim(H) = ∞, see also the discussion in [57].

Exact
Without the condition N i=1 w i C i > 0, this equation always has at least one positive, nonzero singular solution.

Entropic Wasserstein barycenter
3. Sinkhorn barycenter (second version), Without the condition Eq.(47) has the same form as the finite-dimensional version reported in [55] and [45]. It is, however, only applicable for finding the barycenter in the case it is strictly positive, since it is derived under this explicit assumption. If the solution of Eq.(45) is strictly positive, then it is also the unique strictly positive solution of Eq. (47). Eq.(47), however, always has the trivial solution X = 0. Furthermore, if dim(H) ≥ 2 and at least one of the C i 's is positive, then it has uncountably infinitely many positive solutions, which are singular (Proposition 14). The same phenomenon happens for the barycenter equation (50) in the exact, unregularized setting [1], i.e. when ǫ = 0.
As we discuss in detail in Section 11, it is not straightforward to extend the approach in [1] for Eq.(50) and [45] for Eq. (47) in the finite-dimensional setting, which requires all C i 's to be strictly positive for the existence ofC > 0, to the infinite-dimensional setting. This is because it is no longer possible to uniformly lower bound the C i 's by αI for some α > 0 and it is not clear whether this lower bound can be replaced by another strictly positive operator.
We also remark on our condition N i=1 w i C i > 0 for the strict positivity ofC, which is more general than requiring all C i 's to be strictly positive (e.g. [45]).
In fact, we can have N i=1 w i C i > 0, guaranteeingC > 0, with all C i 's being singular (see Section 10 for an example).
The RKHS setting. We now apply the abstract Hilbert space setting above to the reproducing kernel Hilbert space (RKHS) setting. In this case, we obtain an interpolation between Kernel Maximum Mean Discrepancy (MMD) and Kernelized L 2 -Wasserstein Distance. The RKHS formulas are expressed explicitly in terms of the kernel Gram matrices, which are readily computable.
Let X be a complete separable metric space. Let K be a continuous positive definite kernel on X × X . Then the reproducing kernel Hilbert space (RKHS) H K induced by K is separable ( [84], Lemma 4.33). Let Φ : X → H K be the corresponding canonical feature map, so that K(x, y) = Φ(x), Φ(y) HK ∀(x, y) ∈ X × X . Let ρ be a Borel probability measure on X such that Then the RKHS mean vector µ Φ ∈ H K and covariance operator C Φ : H K → H K induced by the feature map Φ are both well-defined and are given by Here the rank-one operator u ⊗ v is defined by (u ⊗ v)w = v, w HK u, u, v, w ∈ H K . Then C Φ is a positive trace class operator on H K (see e.g. [66]).
Let X = [x 1 , . . . , x m ],m ∈ N , be a data matrix randomly sampled from X according to a Borel probability distribution ρ, where m ∈ N is the number of observations. The feature map Φ on X defines the bounded linear operator Φ(X) : The corresponding empirical mean vector and covariance operator for Φ(X) are defined to be where J m = I m − 1 m 1 m 1 T m , 1 m = (1, . . . , 1) T ∈ R m , is the centering matrix, with J 2 m = J m and AJ m is the matrix obtained from the (possibly infinite) matrix A by subtracting the mean column.
Remark 4 To keep our expressions simple, we have assumed that the number of data points in are the same, i.e. m = n. The extension to the case m = n is straightforward.

From Gaussian measures to Gaussian processes
Let us discuss the translation of the results for Gaussian measures on an abstract Hilbert space H into the setting of Gaussian processes, see also [68,32,73,54,57]. Consider the following correspondence between Gaussian measures and Gaussian processes with paths in a Hilbert space, as established in [74]. Let (Ω, F , P ) be a probability space. Let T be an index set. Let (T, A, ν) be a measurable space, ν nonnegative, σ-finite, such that L 2 (T, A, ν) = L 2 (T, ν) is separable (e.g. T ⊂ R n measurable, A = B(T ), ν is the Lebesgue measure). Let ξ = (ξ t ) t∈T = (ξ(t, ω)) t∈T be a real F /A-measurable Gaussian process on (Ω, F , P ), with mean m(t) and covariance function K(s, t), denoted by GP(m, K). The sample paths ξ(·, ω) ∈ H = L 2 (T, ν) almost P -surely, i.e. T ξ 2 (t, ω)dν(t) < ∞, almost P -surely, (63) if and only if ( [74], Theorem 2 and Corollary 1) Then ξ induces the following Gaussian measure P ξ on (H, B(H)) with mean m ∈ H and covariance operator C K : H → H, defined by Conversely, let µ be a Gaussian measure on (H = L 2 (T, A, ν), B(H)). Then there is an F /A-measurable Gaussian process ξ = (ξ t ) t∈T on (Ω, F , P ) with sample paths in H, such that the induced probability measure is P ξ = µ.
Correspondence between covariance function and covariance operator via Mercer Theorem. Covariance functions, being positive definite kernels, can be fully expressed via their induced covariance operators, as follows.
In the following, let T be a σ-compact metric space, that is T = ∪ ∞ i=1 T i , where T 1 ⊂ T 2 ⊂ · · · , with each T i compact. Let ν be a positive, non-degenerate Borel measure on T , i.e. ν(B) > 0 for any open U ⊂ T , with ν(T i ) < ∞∀i ∈ N.
Theorem 16 (Mercer Theorem -version in [85]) Let T be a σ-compact metric space and ν a positive, non-degenerate Borel measure on T . Let K : T × T → R be continuous, positive definite. Assume furthermore that T K(s, t) 2 dν(t) < ∞, ∀s ∈ T. (67) Then C K is Hilbert-Schmidt, self-adjoint, positive. Let {λ k } ∞ k=1 be the eigenvalues of C K , with corresponding orthonormal eigenvectors {φ k } ∞ k=1 . Then where the series converges absolutely for each pair (s, t) ∈ T ×T and uniformly on each compact subset of T × T .
Mercer Theorem thus describes the covariance function K(s, t) fully and explicitly via its covariance operator C K . Since K is positive definite, K(s, t) 2 ≤ K(s, s)K(t, t) ∀s, t ∈ T × T . Thus the condition T K(t, t)dν(t) < ∞ in (64) implies both conditions (67) and (68) in Mercer Theorem and from (69) We now generalize ideas in [68,32,73,54,57], using the fact that Gaussian processes are fully determined by their mean and covariance functions. Most importantly, the following incorporates Mercer Theorem to quantify the correspondence GP(m, K) ⇐⇒ N (m, C K ).
Definition 1 (Divergence between Gaussian processes) Let T be a σ-compact metric space, ν a positive, non-degenerate Borel measure on T . Let H = L 2 (T, B(T ), ν). Let ξ i = GP(m i , K i ), i = 1, 2, be two Gaussian processes with mean m i ∈ H, covariance function K i continuous, and Let D be a divergence function on Gauss(H)×Gauss(H). The corresponding divergence D GP between ξ 1 and ξ 2 is defined to be Mercer Theorem immediately implies the following.
Theorem 17 Assume the hypothesis in Definition 1. Then Remark 5 In our current context, we can immediately apply Definition 1 and Theorem 17 to the exact Wasserstein distances and Sinkhorn divergences. Definition 1 also can also be extended to cover the entropic OT ǫ c distances, however since they are not metrics/divergences, Theorem 17 no longer holds.

Kullback-Leibler divergence between Gaussian measures
In this section, we briefly review the KL divergence between Gaussian measures on Hilbert space and prove Theorem 2. For a separable Hilbert space H, consider the set P(H) of probability measures on (H, B(H)), where B(H) denotes the Borel σ-algebra on H. We focus on the subset P 2 (H) defined by For µ ∈ P 2 (H), its mean vector m ∈ H and covariance operator C : H → H are well-defined and are given by In particular C is a self-adjoint, positive, and trace class operator on H.
We recall that for two measures µ and ν on a measure space (Ω, F ), with µ σ-finite, ν is said to be absolutely continuous with respect to µ, denoted by ν << µ, if for any A ∈ F , µ(A) = 0 ⇒ ν(A) = 0. In this case, the Radon-Nikodym derivative dν dµ ∈ L 1 (µ) is well-defined. The Kullback-Leibler (KL) divergence between ν and µ is defined by If µ << ν and ν << µ, then we say that µ and ν are equivalent, denoted by µ ∼ ν. We say that µ and ν are mutually singular, Theorem [29], [42], which states that two Gaussian measures µ = N (m 1 , Q) and ν = N (m 2 , R) are either mutually singular or equivalent. The necessary and sufficient conditions for the equivalence of the two Gaussian measures ν and µ are given by the following.
2. There exists S ∈ Sym(H) ∩ HS(H), without the eigenvalue 1, such that We now recall results on Kullback-Leibler divergences between two Gaussian If µ ∼ ν, then we have the following result.
For two equivalent Gaussian measures µ, ν on H, the Radon-Nikodym derivative involves only the means and covariance operators ( [65], Theorem 11). This motivates Theorem 2, which seeks to extend the validity of Eq.(78) by generalizing the expression H log ( dν dµ )dν to H log ( dν dµ )dγ where γ ∈ P 2 (H) is any probability measure with the same mean and covariance operator as ν.
To prove Theorem 2, in the following we utilize the concept of white noise mapping, see e.g. [22,23].
For any pair z 1 , z 2 ∈ Q 1/2 (H), we have by definition of the covariance operator Thus the map W : Since ker(Q) = {0}, the subspace Q 1/2 (H) is dense in H and the map W can be uniquely extended to all of H, as follows. For any z ∈ H, let {z n } n∈N be a sequence in Q 1/2 (H) with lim n→∞ ||z n − z|| H = 0. Then {z n } n∈N is a Cauchy sequence in H, so that by isometry, {W zn } n∈N is also a Cauchy sequence in L 2 (H, µ), thus converging to a unique element in L 2 (H, µ). Thus we can define by the following unique limit in L 2 (H, µ) The map W : H → L 2 (H, µ) is called the white noise mapping associated with the measure µ = N (m, Q). W z can be expressed explicitly in terms of the finite-rank orthogonal projections where {e k } k∈N are the orthonormal eigenvectors of Q. For any z ∈ H, we have Thus Q −1/2 P N z is always well-defined ∀z ∈ H. Furthermore, for all x, y ∈ H, The operator The Radon-Nikodym derivative between two equivalent Gaussian measures on H is expressed explicitly via the white noise mapping, as follows.
Let {α k } k∈N be the eigenvalues of S, with corresponding orthonormal eigenvectors {φ k } k∈N . Let W be the white noise mapping induced by µ. The Radon-Nikodym derivative dν dµ is given by where for each k ∈ N, The series ∞ k=1 Φ k converges in L 1 (H, µ) and L 2 (H, µ).
In general, for z ∈ H, in the By the assumption ν << µ, we have supp(ν) ⊂ supp(µ). The sequence {W PN z } N ∈N is a Cauchy sequence in L 2 (H, µ) converging to W z , with a subsequence converging pointwise µ-almost everywhere, hence ν-almost everywhere. This subsequence is also a Cauchy sequence in L 2 (H, ν), converging to a unique limit, with a subsubsequence converging pointwise ν-almost everywhere. Thus this limit must be W z .
For the first expression, for any z ∈ H, N ∈ N, Since ν is a probability measure, we also have lim N →∞ ||W PN z −W z || L 1 (H,ν) = 0 by Hölder Inequality. Thus This completes the proof. ⊓ ⊔ We note that if Φ is an orthonormal basis for L 2 (X, µ), then the hypothesis of Proposition 1 becomes f ∈ L 2 (X, µ) and the conclusion is immediate.
Proof Let S Φ = span{φ k } k∈N be the closed Hilbert subspace of L 2 (X, µ) with orthonormal basis Φ. We show that the following linear functional is well-defined and bounded. First, by assumption, It follows that A f is well-defined on the dense subspace span{φ k } k∈N , with

then by linearity
Using the expression for W φj , W φ k L 2 (H,ν) from Lemma 1, we have Thus g ∈ L 2 (H, γ). Furthermore, for any N ∈ N and Since γ is a probability measure, by Hölder Inequality, we have g ∈ L 1 (H, γ) Letting gives the final answer. ⊓ ⊔ We are now ready to prove Theorem 2.

Mutual information of Gaussian measures on Hilbert space
In this section, we prove Theorem 1. For completeness, we give a new, shorter proof of the mutual information between two Gaussian measures [7]. We start by reviewing joint measures and cross-covariance operators on Hilbert space.
Joint measures. Following [6], let (H 1 , , 1 ) and (H 2 , , 2 ) be two real separable Hilbert spaces with corresponding Borel σ-algebras B(H 1 ) and B(H 2 ). Let H 1 ×H 2 be the Hilbert space with inner product defined by (u 1 , ). Let µ X and µ Y be its marginal probability measures defined on (H 1 , B(H 1 )) and (H 2 , B(H 2 )), respectively. Assume further that µ XY ∈ P 2 (H 1 × H 2 ). Under this assumption, the mean vector m XY ∈ H 1 × H 2 and covariance operator it follows that µ XY ∈ P 2 (H 1 ×H 2 ) if and only µ X ∈ P 2 (H 1 ) and µ Y ∈ P 2 (H 2 ). Subsequently, throughout the paper, we assume that Cross-covariance operators. Let (m X , C X ) and (m Y , C Y ) be the means and covariance operators of µ X and µ Y , respectively. Since µ XY ∈ P 2 (H 1 × H 2 ), the following linear functional G : By the Riesz Representation Theorem, there exist single-valued, bounded, linear maps C XY : Clearly, C Y X = C * XY . In operator tensor product notation, C XY is called the cross-covariance operator of µ XY . It is closely related to the covariance operators of the marginals µ X and µ Y via the following.
Theorem 21 ( [6]) Consider the projection operators P X : Then C XY admits the following representation where V : H 2 → H 1 is a unique bounded linear operator such that ||V || ≤ 1 and V = P X V P Y .
In column vector notation, the mean vector of the joint measure µ XY is given In operator-valued matrix notation, the covari- Mutual information of Gaussian measures. We now consider joint Gaus- is an orthonormal basis for H 1 ×H 2 and 0 V V * 0 Proof By Lemma 3, the operator 0 V V * 0 is self-adjoint and Hilbert-Schmidt on H 1 ×H 2 , thus admits a countable sequence of eigenvalues, with corresponding eigenvectors forming an orthonormal basis for The nonzero eigenvalues of the last operator are precisely the nonzero eigenvalues of V V * : Combining these two equations, we obtain V V * a = λ 2 a and V * V b = λ 2 b. Thus a and b are necessarily eigenvectors of V V * and V * V , respectively, under the same eigenvalue λ 2 . Due to the square factor, it follows that −λ is an eigenvalue By definition of the Hilbert-Carleman determinant, By applying Formula (78) in Theorem 19 and Lemma 4 to our setting, we obtain the expression for KL(µ XY ||µ X ⊗ µ Y ), first proved in [7] (Proposition 2) by a direct approach. The case H 1 = H 2 = H was proved in [5]. As we now show, by employing the more general result on the KL divergence between Gaussian measures in Theorem 19, this result is obtained immediately.
Proof (of Theorem 22) The condition for µ XY ∼ µ X ⊗ µ Y follows from Lemmas 3,4 and Theorem 18. Assume that µ XY ∼ µ X ⊗ µ Y . Since the mean vector for µ XY and µ X ⊗ µ Y is the same, namely m X m Y , the first term in Eq. (78) is equal to zero. From Eqs. (106) and (78), we have where the last equality follows from Lemma 4.

⊓ ⊔
We are now ready to prove Theorem 1, which we restate here for clarity.
Assume that γ and γ 0 have the same covariance operator Γ and that µ X ⊗ µ Y has covariance operator Γ 0 . Then Equality happens if and only if γ = γ 0 . Here V is the unique bounded linear op- we also have γ << γ 0 and the Radon-Nikodym derivatives dγ dγ0 and dγ d(µX ⊗µY ) are both well-defined. Since γ, γ 0 , and µ X ⊗µ Y all have the same mean, namely m X m Y , and γ, γ 0 have the same covariance operators, we have by Theorem 2, The nonsingular case. We first solve (110) for c(x, y) = ||x − y|| 2 under the nonsingular Gaussian setting. Let µ X , m Y ∈ P 2 (H), with means m X , m Y and covariance operators C X , C Y , respectively, then C X , C Y ∈ Tr(H) and Then C XY ∈ Tr(H) and for any A, B ∈ L(H), Proof By Theorem 21, are both Hilbert-Schmidt. Thus it follows that C XY is trace class. Recall that C XY : H → H is given by By Lebesgue Dominated Convergence Theorem, Ax, e j By, e j dµ XY (x, y) If γ ǫ is the minimizer of problem (13), then necessarily γ ǫ ∈ Gauss(µ 0 , µ 1 ) and γ ǫ ∼ µ 0 ⊗ µ 1 .

Lemma 6 Let A : H → H be a compact operator. Then
Proof This is a special case of Lemma 10 and Corollary 2 in [61].
The following result is then immediate.
Lemma 8 Let C ∈ Sym + (H), X ∈ Sym + (H) be fixed. For any constant a ∈ R, a = 0, Proof The desired equality is This last equality is valid as a consequence of Lemma 6, which gives Together with the left hand side of the previous expression, this gives the desired equality.

⊓ ⊔
We apply the following result on the log-concavity of the Fredholm determinant from [59], which is a generalization of Ky Fan's inequality for the log-concavity of the determinant on the set of symmetric positive definite matrices [28].

⊓ ⊔
Let us now compute the Radon-Nikodym density of the optimal entropic transport plan γ ǫ in Theorem 24 with respect to µ 0 ⊗ µ 1 . We note that the optimal V in Theorem 24 is trace class. Consequently, we can apply the following result (see Proposition 1.3.11 in [23] or Corollary 2 in [65]).
where, in the L 1 (H, µ) sense, The following result can be obtained by direct verification.
Eq.(134) is equivalent to Eq. (16) in Theorem 3 via the identity Proof For ||V || < 1, (I − V V * ) is invertible and by Lemma 10, For V = c ǫ (I +(I +c 2 ǫ C 1/2 ǫ , we have V ∈ Tr(H). Similar to the proof of Theorem 24, Thus for any x, y ∈ H and N ∈ N, Combining this with det(I − S) and Proposition 4 gives the desired result. ⊓ ⊔ The general case. We note that while Theorem 24 and Proposition 5 are proved under the hypothesis that ker(C 0 ) = ker(C 1 ) = {0}, the optimal solution obtained is clearly mathematically valid without this assumption. We now confirm that this is indeed the case via a different approach.
Let k(x, y) = exp − c(x,y) ǫ , ǫ > 0. It has been shown that, see e.g. [11,18,25,38,78], problem (110) has a unique minimizer γ ǫ if and only if there exist functions α ǫ , β ǫ satisfying the Schrödinger system In this case, the unique minimizer γ ǫ is the probability measure whose Radon-Nikodym derivative with respect to µ 0 ⊗ µ 1 is given by Motivated by Proposition 5, we now solve the Schrödinger system (136) when c(x, y) = ||x − y|| 2 on H, leading to another proof of Theorem 3, which is valid in the general setting, where C 0 and C 1 can be singular.
We make use of the following results on Gaussian integrals on Hilbert spaces.
The following result then follows immediately Proof (of Theorem 3: optimal entropic transport plan -the general case) We first have The Schrödinger system (136) then becomes Thus for the system of equations (140) to hold ∀x, y ∈ H, we must have We claim that the following A ǫ and B ǫ solve the system of equations (141) Let us verify the first equation in (141) (the second one is analogous). We have Thus in order to have B ǫ = − 2 Replacing the right hand side with the expression in Eq.(143), this is This is precisely Lemma 8 with a = 4 ǫ . Similarly, we have Finally, for the expression of exp(a + b), we note that it is clear that from Eq.(142) that C Here we have used the fact that the eigenvalues of C 1/2 are the same as those of C 1/2 ⊓ ⊔ Lemma 11 Let µ ∈ P 2 (H) with mean m and covariance operator C. Let A ∈ L(H). Then Proof It suffices to consider the case m = 0.
Proof For the optimal Recall that It follows that Here we have invoked Lemma 11. For the first trace term, we have Similarly, the second term is Combining all the previous expressions, we obtain Here we have used the fact that the nonzero eigenvalues of C It follows that D(ϕ ǫ , ψ ǫ ) = ǫE µ0 [log α ǫ ] + ǫE µ1 [log β ǫ ] = OT ǫ d 2 (µ 0 , µ 1 ), as in the proof of Corollary 4.

Corollary 5 (Sinkhorn divergence between Gaussian measures on Hilbert space)
Proof By definition of the Sinkhorn divergence and Theorem 4, This completes the proof.
Using these expressions, we see that the formulas for OT ǫ d 2 and S ǫ d 2 coincide with those in Theorems 4 and 7. Furthermore, it is immediately clear that By L'Hopital's rule, we have ∀x ≥ 0, lim ǫ→0 ǫ log( 1 2 + 1 2 (1 + 16 Combining these limits give the limiting behavior of OT ǫ d 2 and S ǫ d 2 . ⊓ ⊔

Proof of Theorem 15 -RKHS setting. Let λ(A) denote the set of eigenvalues of a compact operator A.
Then Combining these with the expressions for OT ǫ d 2 and S ǫ d 2 in Theorems 4 and 7 gives the desired results.

Remark 6
In Theorem 15, terms of the form AJ m and J m AJ m give the same answer, since they have the same eigenvalues, leading to the same trace and determinant. However, for implementation, the choice of which term to use is important. For computing the log det terms, we usually use the Cholesky decomposition, which requires a symmetric positive definite matrix and hence one needs to use the term J m AJ m , which is symmetric.
We now prove Theorem 6, which we restate as Theorem 26 below.

then we show in Theorems 3 and 4 that
Hence equality happens if and only if µ 0 , µ 1 are Gaussians. ⊓ ⊔

Entropic 2-Wasserstein barycenter
In this section, we prove Theorem 5 on the convexity of OT ǫ d 2 and Theorem 12 on the OT ǫ d 2 -based barycenter. We first recall the concept of the Fréchet derivative on Banach spaces (see e.g. [46]). Let V, W be Banach spaces and L(V, W ) be the Banach space of bounded linear maps between V and W . Assume that f : Then the map f is said to be Fréchet differentiable at x 0 ∈ Ω if there exists a bounded linear map Df (x 0 ) : The map Df (x 0 ) is called the Fréchet derivative of f at x 0 . Let now W = R. If x 0 is a local minimizer for f , then necessarily (see e.g. Theorem 1.33 in [72]) Furthermore, if f is Fréchet differentiable, then (Proposition 3.11 in [72]) for Ω ⊂ V open and convex, Under this identification, In If, furthermore, X 0 is invertible and X ∈ Sym(H) ∩ Tr(H), then In general, let f : L(H) → L(H) be such that f (X 0 ) and X 0 commute, then Proof i) For the function sq(X) = X 2 , we have Let f (X) = X = sq(sqrt(X)), then by the chain rule This gives us the first identity.
ii) For the second identity, we note that Taking the trace operation on both sides gives which gives us the desired expression.
As in the proof of Lemma 18,

Thus it follows from Lemma 17 that
Df (X 0 )(X) = c 2 2 Tr  ||f We remark that the condition that B and C commute is crucial in Lemma 23. It can be verified numerically that, without this condition, the stated inequality is generally false even if A is also positive.
Proof By the assumption that B and C commute, If in addition Y 0 , Z 0 are invertible, then equality happens if and only if X = 0.
Proof By Lemma 17, for any X ∈ Sym(H) ∩ Tr(H), Pre-and post-multiplying both sides by Dsqrt(Z 0 )(X)Y 0 and Y 0 , respectively, Taking trace on both sides and applying Lemma 23 gives The Fréchet derivative Df (X 0 ) : Sym(H) → L(L(H), R) is given by For c = 0, equality happens if and only C 1/2 XC 1/2 = 0. If C is strictly positive, then equality happens if and only if X = 0.
Proof Let h(X) = (I+c 2 C 1/2 XC 1/2 ) 1/2 , g(X) = (I+h(X)) −1 , and f (X)(Y ) = Tr[C 1/2 g(X)C 1/2 Y ]. By the chain rule and Lemma 21, In particular, for Y = X, by Lemma 24, Since Z 0 and (I+Z Then f is convex on Ω. Furthermore, f is strictly convex if C is strictly positive. Proof For any X 0 ∈ Ω, X ∈ Sym(H), by Lemmas 19 and 18, Thus we have the map Df : Ω → L(L(H), R), with Differentiating this map gives the second-order Fréchet derivative where Z 0 = I + c 2 C 1/2 X 0 C 1/2 , by Lemma 25. In particular, for Y = X, Equality happens if and only if C 1/2 XC 1/2 = 0. If C is strictly positive, then equality happens if and only if X = 0 by Lemma 26. Thus f is convex on Ω and furthermore, it is strictly convex if C is strictly positive. ⊓ ⊔

Proof (of Theorem 5 -Convexity of entropic Wasserstein distance)
By the strict convexity of the square Hilbert norm || || 2 , the function m → ||m − m 0 || 2 is strictly convex in m. For the covariance part, by Theorem 4, Since Tr(X) is linear in X and the remaining part is convex in Sym + (H) ∩ Tr(H) by Proposition 6, F is convex in X ∈ Sym + (H) ∩ Tr(H), with strict convexity if C 0 is strictly positive. ⊓ ⊔ Then f is strictly convex. For the minimum of f , there are two scenarios Equivalently, x * is the unique positive solution of the equation which also has the solution x 0 = 0. Proof is a strictly increasing function on [0, ∞).
We have the following three scenarios In this case f ′ (0) > 0 and thus f ′ (x) > 0 ∀x > 0 and thus the minimum In this case f ′ (0) < 0 and thus there exists a unique x * > 0 such that f ′ (x * ) = 0. This is the unique global minimizer of f and We note the last equation also has the solution x = 0. ⊓ ⊔ Proof (of Theorem 12 -Infinite-dimensional setting) As a direct consequence of Theorem 6, the entropic barycenter of a set of Gaussian measures in P 2 (H), if it exists, must be another Gaussian measure on H. By Theorem 4, OT ǫ d 2 (N (m 0 , C 0 ), N (m 1 , C 1 )) decomposes into the squared Euclidean distance ||m 0 − m 1 || 2 and the distance OT ǫ d 2 (N (0, C 0 ), N (0, C 1 )). It follows that we can compute the barycentric mean and covariance operator separately.

The barycentric mean is obviously that the Euclidean meanm =
On Ω, F is Fréchet differentiable. Since we assume that at least one of the C ′ i s is strictly positive, F is strictly convex. Thus a minimizer X 0 ∈ Ω of F must necessarily be unique and satisfy DF (X 0 ) = 0.
(ii) Assume now that ǫI 2 N i=1 w i C i . We show that X 0 = 0 is not a global minimum of F in Sym + (H) ∩ Tr(H). For any X 0 ∈ Sym + (H) ∩ Tr(H), The assumption ǫI 2 N i=1 w i C i means that ∃u ∈ H, ||u|| = 1, such that According to Lemma 27, the property 0 < ǫ < 2 N i=1 w i u, C i u implies that there exists a unique x * > 0 such that Thus X 0 = 0 is not a global minimum of F in Sym + ∩ Tr(H).

⊓ ⊔
Proof (of Theorem 12 -Finite-dimensional setting) When dim(H) < ∞, if we impose the additional condition that ker(X 0 ) = {0}, that is X 0 is invertible, then the identity (180) is equivalent to By Lemma 8, on the left hand side,

Substituting into the previous equation, we obtain
This gives Eq. (42). Define the following map F : Sym + (H) → Sym + (H) by then we have by the monotonicity of the square root function Thus it is clear that (39) has the solution X 0 = 0, which is necessarily unique due to the strict convexity of the entropic OT distance.
(iii) The condition 0 < ǫI < 2 N i=1 w i C i is thus necessary if Eq.(39) is to have a solution X 0 ≥ 0, X = 0. By Proposition 9, if C i ≥ αI, 1 ≤ i ≤ N , and 0 < ǫ < 2α, then Eq.(39) has a strictly positive solution, which is necessarily unique due to the strict convexity of the entropic OT distance.
2. Let u ∈ R n , ||u|| = 1. If 0 < ǫ < 2 N i=1 w i u, C i u , then X u = x u (u ⊗ u) = x u uu T is a fixed point of G, where x u is the unique positive solution of the one-dimensional fixed point equation Hence under these assumptions, the map G as defined in Eq.(183), has uncountably infinitely many fixed points, which are positive but singular.

Proof (of Proposition 7)
For the first part, let γ ∈ R, γ > 0 be such that C i ≤ γI, 1 ≤ i ≤ N . Let β ǫ = α − ǫ 2 > 0 and consider the following set This is a compact, convex set in Sym(n). By Lemma 41, the operator monotonicity of the square root function, and the inequality (1 + a 2 ) 1/2 ≤ 1 + a, a ∈ R, a ≥ 0, we have with c ǫ = 4 ǫ , Thus the continuous map G maps the compact convex set K ǫ into itself. By Brouwer Fixed Point Theorem, G has at least a fixed point in K ǫ .
For the second part, by Lemma 40, Therefore, for X u = x(u ⊗ u), Since (I − u ⊗ u) 2 = (I − u ⊗ u) and (I − u ⊗ u)(u ⊗ u) = 0, by Lemma 40, It follows that Thus the fixed point equation X u = G(X u ) is equivalent to By the assumption 0 < ǫ < 2 N i=1 w i u, C i u , this equation has a unique positive solution by Lemma 27.

Sinkhorn barycenter of Gaussian measures
In this section, we prove Theorem 8 on the strict convexity of the Sinkhorn divergence and Theorem 13 on the barycenter of a set of Gaussian measures on H under the Sinkhorn divergence.
Proof By assumption, lim t→0 This proves the formula for Df (X 0 ). Here we use the fact the Fréchet Differentiability implies continuity. The proof for Dh(X 0 ) is entirely similar.
Proof Let h(X) = (I + c 2 X 2 ) 1/2 and g(X) = (I + h(X)) −1 . By the chain rule, By Lemma 30, with f (X) = g(X)X, Since g(X) and X commute, we also have f (X) = Xg(X) and thus This gives the second, equivalent, expression for Df (X 0 )(X).
Proposition 9 Consider the Gaussian measures N (0, C), N (0, C i ), 1 ≤ i ≤ N . Define the following function F : Then F is well-defined on the larger, open set Ω = {X ∈ Sym(H) ∩ Tr(H) : F is Fréchet differentiable on Ω and the condition DF (X 0 ) = 0 is equivalent to A solution X 0 of Eq.(198) must necessarily satisfy X 0 ∈ Sym + (H) ∩ Tr(H).
Under the additional hypothesis that X 0 > 0, Eq.(198) is equivalent to Proof Let F i (C) = S ǫ d 2 (N (0, C), N (0, C i )). By Theorem 7, Summing over i, 1 ≤ i ≤ N , we obtain the Fréchet derivative for F , namely By Lemma 20, the first order optimality condition DF (X 0 )(X) = 0 ∀X ∈ Sym(H) ∩ Tr(H) is equivalent to DF (X 0 ) = 0, which in turn is equivalent to This gives the first equation. Clearly any solution X 0 of this equation must necessarily satisfy X 0 ∈ Sym + (H) ∩ Tr(H).
Assume now X 0 > 0. We rewrite DF (X 0 ) = 0 as Under the condition X 0 > 0, by Lemma 26, pre-and post-multiplying X 1/2 0 on both sides gives the equivalent expression where the right hand side follows from Lemma 8. This in turn is This completes the proof.

⊓ ⊔
Existence of the Fixed Point. It is clear from Eq.(45) that if it has a solution X 0 , then necessarily X 0 ≥ 0. We now prove that Eq.(45) has at least one solution X 0 , which is then necessarily unique by the strict convexity of the Sinkhorn divergence. This is done via the Schauder Fixed Point Theorem. Let E be a Banach space and M ⊂ E. We recall that a mapping f : M → E is said to be compact if it is continuous and maps bounded subsets into relatively compact subsets of E, that is subsets whose closures are compact. Consider the map F : Sym + (H) → Sym + (H) as defined in Eq. (48). The proof of the existence of a fixed point of F consists of two steps 1. We show that F is compact.
2. Let γ ∈ R, γ > 0 be such that C i ≤ γI, 1 ≤ i ≤ N , and consider the set We show that F maps K into itself. We can then apply Schauder Fixed Point Theorem to obtain the existence of a fixed point of F in K.
Proof If B ≥ 0, then by Lemma 41, we have A * BA ≥ 0. Assume now that A * BA ≥ 0, which means that x, A * BAx = Ax, BAx ≥ 0 ∀x ∈ H. In particular, for y = Ax ∈ Im(A), we have y, By ≥ 0.
It follows that y, By = lim n→∞ y n , By n ≥ 0. Since this holds for all y ∈ H, we have B ≥ 0.
The following result is then immediate.
Corollary 6 Let C p (H) denote the set of pth Schatten class operators on H, 1 ≤ p ≤ ∞. For two operators A, B ∈ Sym + (H) ∩ C p (H), Theorem 28 (Theorem 2.3 in [48]) Let A, B be two positive operators on H and f any operator monotone function with f (0) = 0. Then The following result is then immediate.
Corollary 7 Let A, B be two positive operators on H. Then Proposition 10 Let C ∈ Sym + (H) ∩ Tr(H). The following maps are compact Proof (i) By Corollary 7, we have Thus the map F 1 is continuous on Sym + (H).
Since C 1/2 is a compact operator on H, it maps bounded subsets of H into relatively compact subsets of H. Consider the set Y = {C 1/2 x : x ∈ H, ||x|| ≤ 1} ⊂ H, then Y being relatively compact means that every sequence {y n = C 1/2 x n } n∈N in Y contains a subsequence {y n k = C 1/2 x n k } k∈N that converges in H, that is ∃y ∈ H such that lim k→∞ ||C 1/2 x n k − y|| = 0.
Let {λ k } k∈N be the eigenvalues of C, arranged in decreasing order, with corresponding orthonormal eigenvectors {e k } k∈N . Let N ∈ N be fixed and consider the finite-rank operator C . Then for any ǫ > 0, there exists N (ǫ) ∈ N such that Consequently, for ||X|| ≤ 1, N > N (ǫ), For a fixed N ∈ N, consider the set Z N = {X 1/2 C N } is the desired convergent subsequence, with limit Y N ∈ HS(H). This shows that the set Z N = {X 1/2 C N : ||X|| ≤ 1} is relatively compact ∀N ∈ N in HS(H), so that ∀ǫ > 0, there is a finite ǫ-net in HS(H) such that Consequently, for N > N (ǫ), This shows that the set {X 1/2 C 1/2 : ||X|| ≤ 1} is relatively compact in HS(H), hence in L(H) and thus F 1 is a compact map on Sym + (H). Moreover, each sequence {X where Y ∈ HS(H) is as defined in Part (i). Thus F 2 is a compact map on Sym + (H). (iv) Entirely analogous to F 3 , the map F 4 is compact on Sym + (H).
Proof Let us show that the product map is compact. Let M B be any bounded, non-empty subset of M . By assumption of compactness, each sequence {f (X n ), X n ∈ M B } n∈N contains a convergent subsequence {f (X n1 )} with limit Y 1 in E. Next, the sequence {g(X n1 )} contains a convergent subsequence {g(X n2 )} with limit Y 2 in E. We then have as n 2 → ∞. Thus the set {h 2 (X) = f (X)g(X) : X ∈ M B } is relatively compact. An analogous argument shows that h 2 is continuous, hence h 2 is a compact map on M . ⊓ ⊔ Proof Let c ǫ = 4 ǫ and g(X) = I + I + c 2 ǫ X 2 1/2 1/2 and for 1 ≤ i ≤ N , Using the inequality (1 + a 2 ) 1/2 ≤ 1 + a ∀a ≥ 0, we obtain Applying Corollary 6 twice gives By Corollary 7, Combining Eqs.(214), (215), (217),(216), and (218) gives This shows that F is continuous in the operator norm || ||. In particular, for This completes the proof.
Assume now that X is compact and singular, X = 0. Let {λ k } k∈N be the set of eigenvalues of X, λ k ≥ 0, arranged in decreasing order, with corresponding orthonormal eigenvectors {e k } k∈N forming an orthonormal basis in H. Then λ k e k ⊗ e k and X ≤ γI ⇐⇒ λ 1 ≤ γ.
The barycentric mean is obviously the Euclidean meanm = Then F is strictly convex, since S ǫ d 2 is strictly convex, thus its minimum, if it exists, is unique. By Proposition 9, where F : Sym + (H) → Sym + (H) is the map defined by Eq. (48). By Proposition 12, F is a compact map on Sym + (H). Let γ ∈ R, γ > 0 be such that C i ≤ γI, 1 ≤ i ≤ N and consider the set Then K is a closed, bounded, convex subset of L(H). By Proposition 13, X ∈ K ⇒ F (X) ∈ K. Thus by Schauder Fixed Point Theorem, there exists X 0 ∈ K such that X 0 = F (X 0 ), which must be the unique global minimizer of the strictly convex function F . Clearly, with c ǫ = 4 ǫ , Since 0 ≤ X 0 ≤ γI, 0 ≤ C i ≤ γI, 1 ≤ i ≤ N , Thus it follows that X 0 > 0 ⇐⇒ N i=1 w i C i > 0. ⊓ ⊔

Comparison of barycenter fixed point equations
We now show that for dim(H) ≥ 2, N i=1 w i C i > 0, Eq.(47) has uncountably infinitely many positive, singular solutions.
Consider first the case H = R and C i = σ 2 i , i = 1, . . . , N , has a unique fixed point x * , which satisfies x * > 0. The function has two fixed points, namely x * and x 0 = 0.
Proof The fixed point equation x = F (x) is equivalent to Consider the left hand side, which is f (x) = 1 + 1 + c 2 Thus f (x) is strictly increasing on [0, ∞) and hence there must exist a unique x * > 0 at which f (x * ) = 0 ⇐⇒ x * = F (x * ). Since x * > 0, by Proposition 9, It is obvious that G has another fixed point x 0 = 0. ⊓ ⊔ Lemma 40 Let u ∈ H, ||u|| = 1. Consider the rank-one operator u ⊗ u : H → H defined by (u ⊗ u)x = u, x u. Then for any A ∈ L(H) and any k ∈ N, If A is self-adjoint, positive, then [(u ⊗ u)A(u ⊗ u)] 1/k = u, Au 1/k (u ⊗ u).
Proof Consider the cases k = 1 and k = 2. For any x ∈ H, (u ⊗ u)A(u ⊗ u)x = (u ⊗ u)A u, x u = u, x u, Au u = u, Au (u ⊗ u)x.
Since (u ⊗ u) 2 = (u ⊗ u), we have For the first expression, the general case then follows by induction.
If A is self-adjoint, positive, then (u ⊗ u)A(u ⊗ u) is self-adjoint, positive, and [(u ⊗ u)A(u ⊗ u)] 1/k is well-defined and unique [12]. We have from which the second identity follows. ⊓ ⊔ Proposition 14 Let C i ∈ Sym + (H), 1 ≤ i ≤ N be fixed. Consider the following map G : Sym + (H) → Sym + (H), defined by Then X 0 = 0 is a fixed point of G. Let u ∈ H, ||u|| = 1 be such that N i=1 w i u, C i u > 0. Then X u = x u (u ⊗ u) is a fixed point of G, where x u is the unique positive solution of the following one-dimensional fixed point equation The condition N i=1 w i u, C i u > 0 is satisfied for at least one u ∈ H, u = 0, since otherwise C 1 = · · · C N = 0. If N i=1 w i C i > 0, then N i=1 w i u, C i u > 0∀u ∈ H, u = 0. Thus under this assumption, for dim(H) ≥ 2, G has uncountably infinitely many fixed points of the form X u = x u (u ⊗ u).
Therefore, for X u = x(u ⊗ u), I + c 2 ǫ X 1/2 u C i X 1/2 u = (I − u ⊗ u) + 1 + c 2 ǫ x u, C i u (u ⊗ u). Since (I − u ⊗ u) 2 = (I − u ⊗ u) and (I − u ⊗ u)(u ⊗ u) = 0, by Lemma 40, Applying the same argument and using the fact that N i=1 w i = 1, we have Thus the fixed point equation X u = G(X u ) becomes As shown in Lemma 39, under the condition N i= w i u, C i u > 0, this onedimensional fixed point equation has a unique positive solution x * u . Thus X u = x * u (u ⊗ u) is a fixed point of G. ⊓ ⊔ When ǫ = 0, the fixed points X u of G in Proposition 14 admit a closed form.
Proposition 15 Let C i ∈ Sym + (H), 1 ≤ i ≤ N be fixed. Consider the following map G : Sym + (H) → Sym + (H), defined by Then X 0 = 0 is a fixed point of G. Let u ∈ H, ||u|| = 1 be such that N i=1 w i u, C i u 1/2 > 0, then the following is a fixed point of G The condition N i=1 w i u, C i u 1/2 > 0 is satisfied for at least one u ∈ H, u = 0, since otherwise C 1 = · · · C N = 0. If N i=1 w i C i > 0, then N i=1 w i u, C i u > 0∀u ∈ H, u = 0. Since w i > 0, i = 1, . . . , N , for each u = 0, there must be at least one i for which u, C i u > 0. This implies that under this assumption, N i=1 w i u, C i u 1/2 > 0 is satisfied for all u ∈ H, u = 0. In this case, when dim(H) ≥ 2, G has uncountably infinitely many fixed points of the form X u .
One can immediately see that some of the solutions of the above equation include X = 0, X = C, X u = u, Cu (u ⊗ u), for any u ∈ H, ||u|| = 1, including X k = λ k (e k ⊗ e k ), k ∈ N, where {λ k } k∈N are the eigenvalues of C, with corresponding orthonormal eigenvectors {e k } k∈N .
Thus the fixed point equation X u = G(X u ) becomes by the assumption N i=1 w i u, C i u 1/2 > 0. Hence Assume that there exist α, β ∈ R, α > 0, β > 0 such that αI ≤ C i ≤ βI, 1 ≤ i ≤ N . Consider the set K 2 = {X ∈ Sym ++ (n) : αI ≤ X ≤ βI}. Then by Lemma 41 α 2 I ≤ αX = X 1/2 αX 1/2 ≤ X 1/2 C i X 1/2 ≤ X 1/2 βX 1/2 = βX ≤ β 2 I. It is not clear, however, whether this proof strategy can be extended to the infinite-dimensional setting. This is because one can no longer assume that there is a uniform lower bound of the form αI for the C i 's and X as above.
One might assume instead that there is an operator C > 0 such that C i ≥ C, 1 ≤ i ≤ N and consider the set K 3 = {X ∈ Sym + (H) : C ≤ X ≤ βI}. This does not help, however, since the condition X ≥ C > 0 does not imply that X 1/2 CX 1/2 ≥ C 2 . The following is a counterexample (232) It can be verified numerically that X ≥ C but X 1/2 CX 1/2 C 2 .
We note also that the condition 0 < αI ≤ C i ≤ β I , 1 ≤ i ≤ N , above is more restrictive that the condition N i=1 w i C i > 0 stated in Theorem 13, which guarantees the existence of a strictly positive solution of Eq. (47). As an example, let H = R N and {e i } N i=1 be an orthonormal basis in H. Define guaranteeing thatC > 0, even though all the C i 's are singular.

Miscellaneous Technical Results
Operator monotone functions. We recall the concept of operator monotone functions. Let A, B be two self-adjoint bounded operators on H, then we say . Some of the well-known operator monotone functions are the following.
The following is the generalization of Proposition V.1.6 in [8] to the infinitedimensional setting, with the additional assumption that A, B be invertible, since if dim(H) = ∞, then A > 0 does not imply that A is invertible.


Evaluation of the Instantaneous Power Demand of an Electric Charging Station in an Urban Scenario

: The existence of an efficient and sufficiently extensive charging infrastructure network appears to be of vital importance for the widespread acceptance of electric mobility by users. The present work aims to develop a tool based on big data analysis that helps to deploy a network of charging stations which can efficiently serve the potential demand, both from the user side, improving the level of service for charging and to cover the territory in a satisfactory way, and from the business side, allowing an analysis of the potential power load. The paper integrates real world traffic data and the results of an experimental campaign on an electric vehicle to evaluate the instantaneous power demand of a fast charging station, based on a procedure for the evaluation and proper time allocation of each charge request.

Introduction
Battery electric vehicles (BEV) can contribute to the improvement of air quality in urban areas, especially as regards NOx and particulate matter (PM) pollution. However, electric vehicle deployment is hindered by many factors [1], and the diffusion of fast charge networks is fragmented, especially in emerging markets. One of the obstacles to the large-scale diffusion of purely BEV is the traction battery, which has so far presented limitations due to high costs, reliable determination of the state of charge, and density of energy [2]. In fact, despite the continuous technological progress of the batteries, the potential users of the BEVs still perceive their limited range as a problem [3,4]. One possible solution is the development of an efficient and sufficiently extensive charging network, which could help reduce charging anxiety [5,6]. A survey in California showed that fast charging is not exclusively used in long distance journeys. The results show that 50% of quick recharges made by Nissan Leaf owners take place within 10 km from homes and are carried out by 34% of users, while 37% of Bolt users recharge within a distance of 12 km from their home [7]. Among the users who use fast charging more often there are those who do not have a home charge. Similar results were obtained in a study in the UK and Ireland, which shows that the median distance traveled on journeys using rapid recharge is 50 km, while the average is 61 km [8,9]. A user survey in Norway has shown that BEVs get about 4-6 percent of their energy from fast chargers and that the average fast charge power was 30.5 kW in 2017, although chargers can deliver 50 kW [10]. The average time of a charge event was 20.3 min and the average energy was 9.6 kWh, and both these quantities follow a normal distribution around these values. From these studies, the complexity of charging and driving behavior of BEV users emerges, due also to differences among different urban environments and also rural areas [11].
In the present work, we use the results obtained in the demand-side analysis presented in [12,13] to determine the charge demand deduced from actual traffic data in an urban zone of Rome; moreover, experimental data coming from a specific series of on-road tests performed with an electric vehicle to evaluate typical and realistic charge profiles are used. Based on the association of these data, and through the integration of the assignment model presented in [13], a method was carried out for evaluating the instantaneous power load of a charging station over a defined time period and consequently its maximum, minimum, and average values. The novelty of our method resides in the fact that instantaneous (and not average) values are obtainable; the estimation of these operational aspects is of paramount importance for station investors/operators, because they generally draw long-term contracts with the power utility where average and peak power level is agreed in return for a lower price.
This information can also be used for an economic analysis of the investment viability and of the necessary considerations regarding station configuration, because the elements derived from our method are of central importance in these aspects.
In this work, only fast charging stations are considered, which, in our opinion, could have a boost in short and medium-term development, due to the similarity with the traditional refueling of internal combustion engine (ICE) vehicles [14], although, it should be borne in mind that the payback period for a fast charging station may be more than ten years [15]. However, further analyses based on real data could change these predictions.
This work will also illustrate the impact of widespread BEVs and fast charging stations on the electricity network [16], which may require a review of the configuration and performance of the actual urban electricity network.
The structure of the paper is as follows: in Section 2 we give a brief overview of the literature about the charging infrastructure allocation problem; in Section 3 we present the methodology and data used in the paper; in Section 4 we illustrate the results and Section 5 presents the conclusions.

Literature Review
Several design tools have been proposed in the literature in order to solve the problem of an effective charging network.
In [17], the feasibility of a BEV-enabled parking lot, that is a parking lot equipped with chargers that acts also as an energy trader, is evaluated using an agent-based simulation. They propose several key performance indicators for the investment strategy and investigate a number of possible scenarios. In [18], a methodology that uses an amalgam of a variety of business and geographical data is proposed, to foresee the charging station utilization when some condition in the charging infrastructure or around it changes and to optimally locate a new charging station. In [19], an approach is proposed for the optimal placing and sizing of fast charging stations that minimizes the total cost of station development, including electrification costs, and the BEV electric loss to reach the station. The stations are located along urban roads, and the resulting mixed-integer non-linear problem is solved using a genetic algorithm. Different scenarios are reported to validate the model.
An analytical method to determine the location of charging stations is proposed in [20], where the optimization aims to minimize the travel cost of each BEV, including its charging time, and to equalize the electric demand for each charging station. However, the approach does not investigate the nature of the input data for the model. An approach that includes the concept of quality of service, related to the waiting time for a charging operation, is proposed in [21], while in [22] a distance satisfaction function is proposed. The work focuses on long BEV trips and takes into account the charging reliability of the system by including in the optimization problem the Euclidean distance among the BEV actual positions at the moments of the charge request and the charging stations.
The impact of operational uncertainty and variability of the electric grid to supply the required power to an electrified transportation network was investigated in [23]. In order to evaluate the impact, the capacity of the charging stations is determined using origin-destination (OD) traffic flow, where the average arrival rate at the charge stations has been considered. However, using average values can lead to an incorrect estimation of some quantities, as we will see later. The impact of a large-scale diffusion of electric vehicles on the urban distribution network was investigated in [24] using a spatial-temporal model to obtain the EV charging load of each load busbar over time throughout a day. OD analysis from intelligent transportation research was performed to model the EV mobility. However, the work does not consider the charging infrastructure location problem, since the authors consider that all users start charging immediately after coming back from their daily trip, with the possibility of "dumb" charging, which starts as soon as the vehicle is plugged in, or "smart" charging, which has a management system that takes into account the distribution network parameters.
It is indeed important to analyze carefully BEV driving patterns in order to properly estimate the energy demand [25]. The adoption of an electric vehicle involves an adaptation phase that depends on several factors, such as the skills to drive the electric vehicle and take advantage of regenerative braking, recharging strategies and travel planning to take into account the limited range as well as the interaction with other road users [26]. A study conducted by Rolim et al. [27] shows that the adoption of an electric car has an impact on driving style and, to a lesser extent, on lifestyle. However, the study was conducted on a very small sample (11 participants in the survey). In this study, it is noted that EV cars are mainly used in commuting and to run small errands, and most of the journeys are made in the urban area. The influence of daily activities on BEV users and on their charging behavior has been highlighted in different studies [28,29] showing how to minimize the missed trips or when and where BEVs can be recharged at minimum cost for their owners. Other researchers focus on the development of software applications that help BEV owners in the optimization of the charging process. In [30], to determine an optimal charging strategy, the stochastic nature of driving patterns is described by an inhomogeneous Markov model, using data collected from the utilization of an electric vehicle. A mobile application to communicate between the BEV and the charging station server is proposed in [31]. The details of the electric vehicle are sent to the station server, which in turn indicates the available parking slots to the users. The users can then choose a charging station among the ones within their residual range, thus assuring timesaving and a higher service quality. Another application that manages the BEV charge scheduling has been proposed in [32]. This interactive user application combines estimated battery parameters and data on slot availability from the charging stations. The user can then book an available slot at the most convenient station in terms of money, distance, and waiting time.
The problem of electric power and energy sizing and management for charging infrastructures is highly complex, and different strategies can be implemented to interact optimally with the grid resources [33]. For this purpose, it becomes very important to estimate the instantaneous energy and power demand of the charging stations [34]. This depends on several parameters, such as the number of users requesting recharging, the battery state of charge (SOC) of each vehicle charging, and the duration of each charge [35], which can vary depending on the vehicle model and the type of charger [36].

Materials and Methods
The aim of the present work was to obtain the time profile for the electric power load of a charging infrastructure (CI) starting from its crowding time profile, properly considering each user power demand in the time series. We define the CI as a charging station in which one or more charging points (CPs), which represent the single pole where a car can be connected to recharge, can be present.
The salient steps of the procedure presented here are the following: 1. We made a preliminary choice of an area of Rome and selected a sample of drivers whose trips ended in the selected area and that potentially needed an electric charge. The sample was extracted from a big data set of private vehicular traffic in Italy derived from the Octo Telematics database [37]. 2. We assigned each charge request to a CI, included in a pre-set scenario of stations and CPs. The assignment was conducted through a fuzzy model based on vehicle-station distances and on the actual waiting time for the charge, following a criterion of users' convenience. 3. For each charge event, we calculated the required electric load, dependent on battery SOC as a function of time, by interpolating the experimental profiles. 4. For each CP, the instantaneous electric load profile was obtained from the sequence of the charge profiles of all the vehicles assigned to it, strictly respecting the arrival times of the vehicles at the CP.
5. By adding the instantaneous load of each and all CPs, we obtained the total power and energy demand of each station.
As for the traffic data in point 1), the GPS data were recorded by the Octo Telematics Company. The data represent journeys made by privately owned ICE vehicles equipped with an acquisition device, which records and sends a series of trip information collected by GPS to a central database.
The sample corresponds to about 6% of the Italian private vehicle fleet. To obtain our sample, we applied a filtering procedure to the traffic data [13], to obtain all the travel trajectories that ended in the chosen area, regardless of their starting point, in a given week of May 2013 (from 6 to 12 May). It was hypothesized that the sample considered consisted entirely of electric cars and, in particular, of the medium segment, such as the Nissan Leaf. Since the average effective range of the Leaf is around 160-180 km [38], we selected only the trajectories whose length was between 70 and 140 km, when presumably the user begins to consider a charge, due to range anxiety [39,40]. Although the sample was registered in 2013 and for ICE vehicles, the urban patterns can be considered reliable enough to be used to build a hypothetical electric mobility scenario, on the basis that urban patterns would not change too much for conventional and EV cars [41]. The sample served as a numerical validation for the proposed model. However, the method is independent of the data source and can be fed with more recent data or related to other areas.
The results of this filtering procedure, which identified 3073 trajectories, are shown in Figure 1  In two previous works [12,13], a method for determining the size and position of electric charging stations in an urban area was presented, which also considered the need of the user to find a CP close to his usual destinations. The method was based on the study of traffic flows and users' final destinations in the city of Rome. In particular, we introduced some hypothetical scenarios for the charging stations, with different numbers and distribution among stations of CPs. Figure 1 shows the position of the three charging stations (chosen from a list of traditional filling stations existing in the area and marked with red pins) used in the study. By varying the number of stations and number of charging points, it is possible to construct a multitude of scenarios.

Modeling Users Behaviour
After selecting one of the possible scenarios, we used a fuzzy model to mimic the human decisionmaking process of choosing a station by a user. The fuzzy method allows you to assign an ʺattractiveness scoreʺ to a specific station for a hypothetical BEV driver. In our work, we considered Mamdaniʹs fuzzy inference method, implemented in Matlab® (The MathWorks, Inc., Natick, MA, USA).
The fuzzy model used for the attractiveness score consists of the following elements: A. Input and output variables The input variables are (see Figure 2): 1. the waiting time for the charge event in a given station, which represents the time lapse between the parking time and the start of the charging event. This value is known to the user via some kind of communication, e.g., a mobile application. The variable values between 0 and 30 min are normalized with a linear function that varies between 0 and 10. If the waiting time is longer than 30 min, it is normalized to 10 (most disadvantageous score); 2. the distance of the CI from the actual parking position. The variable is normalized with a linear function that varies between 0 (corresponding to 0 distance) and 10 (corresponding to the diagonal of the rectangular area selected), which is the most disadvantageous score.
The output is the station attractiveness score, which varies between 0 (lowest score) and 10 (highest score). The model can in principle include other input variables that can determine the attractiveness score, but we will not examine any further parameter, for the sake of clarity of the model interpretation.

B. Set of rules and membership functions (MFs)
The logic structure of the model is represented by a set of rules, while the membership functions (MFs) are the weights for the input variables, which can be expressed by any function, but must vary between 0 and 1. The function is usually chosen by the model designers on the basis of their experience. The MFs for the inputs in our case are two sigmoidal functions, which attribute a higher membership degree to lower values of the input, while for the output a triangular form was chosen as MF.
We chose a single simple rule for the model: If (Distance is MFDis) and (Waiting time is MFWT) then (Score is MFSc). (1) The rule corresponds to the AND operator: the results are the minimum between the input variable values. Then, to obtain the score, the "smallest of local maximum" rule is applied, according to which the output score is identified with the abscissa of the smallest among the local maximum (see Figure 3). The fuzzy model is applied to each trip in the sample, taking into account the temporal distribution of arrivals for the entire week. The first request for a charge in the selected week coincides with the start time of the simulation. According to the fuzzy model, each charge request is assigned to a station, obtaining a list of charge requests for each CP of each CI, which is updated for each new request. The waiting time is calculated accordingly to the CP occupancy for each station. The result of the simulation is the list of charging events and the corresponding waiting times.
The results of this simulation give us an estimation of the crowding profile for the stations in each scenario, reflecting the time sequence of the charge events, and allow us to determine the level of service (LoS), which we define as being inversely proportional to the average waiting time. For a set of predefined scenarios, this method allows the best one to be defined based on the value of the ratio ʺ(LoS) / (total number of Cps)/(number of CIs)ʺ.
For example, let us consider the scenario consisting of three stations distributed in the studied zone, as shown in Figure 1, with the following allocation of CPs: (10,10,20) for the stations (#1, #2, #3) respectively. Let us fix the fast charge duration in 30 min. The engagement of the CPs for this scenario is shown in Figure 4 for the entire week duration (expressed in hours). The LoS for this scenario is 94.37%. Different distributions of the CPs can change the LoS. For example, for the scenario (14,12,14) of CPs the LoS reduces to 87.6%. This represents a worsening compared to the (10, 10, 20) scenario due to a different redistribution of the users. Moreover, in general, lowering the total number of CPs lowers the LoS.
However, the simulations show that the dynamic of the system is quite complex, and a lower number of CPs does not always reduce the LoS. Indeed, for the distribution (8,8,16) the LoS is 87.7%, thus comparable to what was obtained for the previous scenario with 40 CPs. This shows that enhancing the LoS is not always necessary to increase the number of CPs.
We also investigated a different configuration of CIs. As an example, in Figure 5 we show the engagement when only stations #2 and #3 are selected, with a CP distribution (0, 12,20). The LoS for this scenario is 89.33%, which is even higher than what was obtained for the same number of CPs distributed among three stations. Next, we analyzed the sensitivity of the model with respect to the charging duration. In Table 1 we compare the results for the maximum waiting time and LoS obtained using 30 and 23 min as charging session durations. The last value is the average of the mean values for Norway and Sweden, according to [10]. The influence of the charging session durations on the maximum waiting time is more important for the scenarios with high CP availability: As a further investigation, we considered distribution for the charging session duration. In particular, we compared four probability distributions for the average CP availability scenario (8,8,16): 1. Poisson distribution with mean value λ = 1380; 2. Weibull distribution with parameters k = 2, λ = 1; 3. Weibull distribution with parameters k = 3, λ = 1; 4. A simple distribution for which the charge duration is inversely proportional to the residual battery charge or depth of discharge (DoD = 1-SOC).
The Poisson distribution is a discrete probability distribution, representing the probability for the occurrence of stochastic independent events whose average time between events λ is known. For high values of λ, the Poisson distribution tends to the Gaussian distribution, which is symmetrically distributed about the mean value. The Weibull distribution is a continuous distribution that is positive only for positive values of x and is zero otherwise. It is often used in the failure analysis and is characterized by two parameters: the shape parameter k > 0, and the scale parameter λ > 0. For this distribution, the median is smaller than the mean, which is consistent with the data reported in [10]. The two distributions considered have different skewness. The last one is an empirical distribution that assumes that the higher the SOC the shorter the charge duration, which is reasonable although not always true in the real world.
All these are simple approximations of the complex charging behavior of EV users. The results for the maximum waiting time and the LoS for the proposed distributions are presented in Table 2. For options 1 to 3, the results are the mean values of 50 trials, while for option 4 we performed a single trial since it is a deterministic function. The Weibull distributions, which are the most left skewed among the ones considered here, show the smallest values for the maximum waiting time, as is intuitive. Surprisingly enough, the LoS, which is related to the average waiting time, does not show the highest value for the Weibull distributions but for the simple distribution proportional to the residual battery charge, which is derived from the travelled distance and under the hypothesis of a fully charged battery at the beginning of the journey. These data show the need for further investigation of the charging behavior of EV users in the CI modeling. Due to the lack of more detailed information, we will use a fixed charging duration scenario in the following, as was previously used.
The problem of charging station planning presents a contrast between the interests of investors, who aim to lower initial spending by minimizing charging points, and the users, who wish to receive the service in the shortest possible time. The overlapping of charge requests, due to a limited number of CPs, creates inconveniences for users due to long waiting times, thus reducing the LoS provided. To simulate a more realistic scenario, we add the possibility for the users to decide whether to recharge their cars or not, based on the estimated waiting time at the stations and the car battery SOC value. To this end, another fuzzy model is added to the simulation. The input values are the waiting time and the depth of discharge (DoD) of the car battery. The set of rules is represented by the single rule:

If (Waiting time is MFWT) or (DoD is MFDoD) then (Score is MFSc).
( The rule and the MFs for the waiting time (MFWT), the DoD (MFDoD), and the score for the decision (MFSc) are graphically represented in Figure 6. The output of the OR operator is the maximum between the input variable values, and the operator for the defuzzifying rule is the "smallest of local maximum" of the MFSc, as was for the previous model. Using this fuzzy decision model, designated as the time-out (TO) model in the following, we allow the possibility for the potential users to decide not to charge when stopping if their SOC is sufficiently high and the waiting time at the CI is too long. Conversely, if the SOC is low or there is no or little queue at the CI, the driver would charge. Tables 3 and 4 present some figures that compare the performance of the original model with the TO one. For the 30 min charge duration, in all the scenarios analyzed the loss of users is around 0.9% and 2% of the total for the scenarios with high and low CP availability, respectively, and concerns primarily the peak of arrival on Monday. On the other hand, the maximum waiting time for each scenario decreases and the LoS increases, although they remain unsatisfactory in the low CP availability scenarios. For the 23 min charge duration, the user loss is around 0.3% and 0.9% for high and low CPs availability, respectively, and the obtained level of service is always very high.
This sensitivity analysis shows that the mean duration of the charge strongly influences the model outcome. In addition, a suitable dislocation of lesser CPs can be more efficient than a larger number not satisfactorily allocated. All of these aspects should be taken into account when planning a charging infrastructure in a predefined scenario.
In the present work, the approach is completed with a detailed analysis of the power demand of a station, which is important information that influences the station performance and operation.

Fast Charging Profiles for Different SOCs
Currently, the international standards prescribe a direct current charge (Mode 4 IEC 61851-1 or SAE J1772 DC level 2) for fast (and ultra-fast) charge, because the onboard conversion would be very costly, heavy, and cumbersome.
The charging process starts from a punctual check of the conditions presented by the vehicle battery through the battery management system (BMS). This control is guaranteed by a vehiclestation communication protocol that enables charging operations after a correct recognition of the vehicle. Moreover, since the protocol structure is of the master (vehicle)-slave (station) type, it is always the vehicle BMS that determines the battery charge conditions, according to the battery state of charge (SOC).
The fast charging profile is generally defined as a constant current-constant voltage (CC-CV) profile: at the beginning of the charge process, the battery accepts a constant current at its maximum value, set by the manufacturer and by the standard of supply (CC phase). When the voltage of the battery reaches a certain value, or when one or more cells of the package have reached a limit condition for voltage or temperature, the charging current is gradually decreased, so as to maintain the constant battery voltage (CV phase). If the initial SOC is high, the CC fraction may be short or even absent, because the battery cells quickly reach the threshold value at which the transition to CV charging takes place. The total duration of the charge will depend on a number of factors such as initial charging status, battery conditions, and temperature. The power profile of the CC-CV charge is relatively constant in the initial part and falls during the constant voltage phase, when the value of the current supplied by the battery charger decreases. The average power accepted by the battery during fast charging depends on the battery SOC.
In order to better investigate the charge power profiles for different initial SOCs and the battery behaviors, a series of experimental tests on a commercial electric vehicle of medium size (a 2014 Nissan Leaf, with a 24 kWh nominal battery capacity) was carried out at ENEA laboratories using fast charging research equipment.
The vehicle was equipped with a data acquisition system to detect the most significant battery parameters (total voltage and current, cell temperature, SOC, individual cell voltage) and road tests were conducted. A Controller Area Network (CAN) interface extracted and registered the values of interest from the CAN bus messages between vehicle and charging equipment.
The charging station used is a 50 kW CC circuitor station that uses the CHAdeMO EQC-50 standard. A LAN interface allows access to the setting and control section (see Figure 7), which are stored in a PC. The Nissan Leaf mileage is around eight thousand kilometers on urban and extraurban roads. The relevant electrical parameters for fast charges have been recorded for several different initial SOCs, with a sampling frequency of 1 Hz. As an example, Figure 8a shows the battery voltage (red line) and absorbed current (blue line) for a 2000 s fast charge and initial SOC of 21%. The initial plateau for the current (constant current phase) and the subsequent rapidly decreasing current phase are clearly visible; the battery voltage shows an initial rapid growth, followed by a constant phase at the final value. As already mentioned, at high SOC values the BMS can decide to charge the battery without the initial constant current phase. Figure 8b presents three different charging current profiles at about 21%, 40%, and 55% initial SOC and shows a variety of trends: the constant current segment is present only in the lowest SOC profile (21% SOC). The absence of the constant current phase is due to the conditions of the battery cells that require a charge current lower than the maximum deliverable from the station.
In Table 5, some electric parameters pertaining to fast charge events for another set of initial SOCs are presented. Here, Pm and PMAX are the average power and the maximum power delivered during the charge, respectively. Looking at the data, we can deduce that the battery is not fully charged after 30 min. For example, for an initial 20% SOC, corresponding to a residual energy of 4.8 kWh for the battery and considering a nominal capacity of 24 kWh, the total absorbed energy is 12.5 kWh, so the final SOC is about 72% after a 30 min fast charge. From the measurements, we selected a set of fast charge time profiles (in terms of absorbed power per second) with initial SOCs very close to the four chosen values, 20%, 40%, 55%, 74%, and we grouped them accordingly. We then averaged the time profiles in each group in order to obtain the four average charge profiles relative to the chosen SOC values and sampled each second for a total duration of 1800 s. Next, we applied an interpolation procedure among these four reference profiles to determine the absorbed charge power relative to any SOC in the 35%-70% range of the sample. In the present work, the interpolation is done through cubic spline functions with the csapi routine in Matlab environment. Iterating the process for each timestep, we obtain the entire charging power profile over time, as shown in Figure 9, where the interpolation for SOC = 30% is reported. Dedicated software in Matlab was developed in order to create a complete charge profile for the Nissan vehicle, for any given initial SOC value included in the range (35%-70%).

Results
In this section, an evaluation is shown of the instantaneous power load at a given CP as an effect of the charge requests. In this process, a fundamental constraint has to be respected: the sequence of the charge requests at each CP must be properly presented respecting the arrival order and avoiding overlaps during the sequential charges.
After the drivers' charge request assignment to a CP of one of the stations in a defined scenario, according to the methodology previously illustrated, we developed a software procedure that calculates the actual SOC from the travelled distance (under the hypothesis that they started their journey with a fully charged battery). It then determines the individual charge power profiles through interpolation and, based on each parking time, properly collocates the request on the time axis of the week 6 to 12 May 2013.
We are then able to obtain the electric power load of a CP as a function of time.
We consider three stations distributed in the studied zone as shown in Figure 1, and the scenarios presented in Table 4 for a fixed charge duration of 23 min, with and without the time-out option. For each scenario illustrated in the previous section, our software procedure allows the total electric power load profile of the charging stations to be built. For the sake of simplicity, we illustrate the procedure for station #3 of each scenario. In Figure 10 the time histories of the power load at CPs number 1, 15, and 20 pertaining to station #3 are presented for the scenario (0, 12,20), as a result of our software procedure; the profiles are referred to Monday 6 May 2013, in the time interval 12:00 to 24:00. The individual charge power profiles are reported in chronological order, and the different form of the profiles is due to the fact that each vehicle has its own initial SOC which determines a different power absorption.
We can observe a situation of almost complete saturation for CP number 1 and a situation of low duty for CPs number 15 and 20. In fact, according to our assignment rules of drivers to a station, the procedure saturates first the CPs of lower numbering, then moves toward the higher ones when necessary. For comparison, Figure 11 shows the power charge profiles for CPs 1 and 15 for the scenario (0, 10,15). The crowding level is higher, and the LoS is lower, but the power profiles are quite similar. The power between two charges does not go to zero because we assume a null switch time between two consecutive users. In Table 6 the maximum power load peak values at station #3 are reported. It can be seen that a minor availability of CP redistributes the power load in time, with the peak values of the electric power lowering with the number of CPs, but with a worsening of the LoS, as expected. The total energy supplied for the electric charges is identical in both the situations, but with a different time distribution. Except for the highest availability scenarios, there is an increase in the peak power load for station #3 when the TO is allowed. This is due to a different redistribution of the charging demand among the available charging stations. Indeed, the total energy demand is lowered when the TO for users is allowed, as can be seen in Table 7. In the model, since the charge time is the same whatever the SOC is, the service provider does not apply any SOC-based priority to the request. However, this aspect could be further investigated in business model oriented research. The sum of all the time-dependent load profiles of the CPs at station #3 gives the second by second station power load profile for the entire week. Figures 12 and 13 present the two profiles, which refer to the situations of high and low availability of CPs, respectively. It is interesting to note, by comparison between the two figures, how the latter load profile decreases and flattens out as a result of the increased time delay in serving the charge request.  In Figures 14 and 15, the number of drivers served by each CP of station #3 is presented by means of two bar diagrams referred to as scenario a) and scenario b), respectively. We can note from Figure 14 how the CPs beyond number 15 for station #3 in scenario a) take a poor part in charge operations, and due to redundancy, an excellent LoS remains assured (see [13]); from Figure 15 we see that the 15 CPs of station #3 in scenario b) are more crowded. This is due to the reduction of the availability of CPs (which passes globally from 80 to 25), which does not assure an acceptable LoS, causing time overlapping between charge requests and consequent spooling (due to the assumed rule of 'first arrived first served'). A key strength of our methodology lies in the fact that it is easy to evaluate the ratio between the LoS (in terms of inverse of the average or maximum value of waiting time) and the level of investor profitability (depending among other things on the costs of the activated charging points); in fact, investment costs can be contained by limiting the number of CPs at the expense of the provided LoS. Another important aspect is that the output is very sensitive to the behavior model for the EV users, which needs to be thoroughly investigated.

Conclusions
The paper presents a methodology to evaluate the electric demand and the expected quality of service for charging infrastructure in a given scenario, based on different simulation models for users' behavior. Several scenarios with a different number of stations and charging points have been analyzed. In general, limited availability of charging points redistributes the power load in time, resulting in a lowering of the peak power and a worsening of the quality of the service for the charging stations. Indeed, when more CPs are available, more EVs can be served at the same time, resulting in a rise of the peak demand power. We also investigated the quality of service for different distributions of CIs among the stations, which turns out to be very sensitive to the distribution of CIs and the number of total charging stations. Charging duration influences the overall level of service for any scenario, and passing from a fixed charge duration of 30 to 23 min can enhance the level of service up by 18% in the low charging points availability scenario. Furthermore, if we allow for a time distribution of charge event durations, the results show a better level of service. However, further analysis is needed to find the most appropriate distribution function. As a next step, we included the possibility for the users to decide whether to recharge their cars or not, based on the estimated waiting time at the stations and the car battery SOC value. We found that the loss of users is around 0.9% and 2% of the total for the scenarios with high and low CP availability for a fixed 30 min charge duration, which goes down to 0.3% and 0.9%, respectively, when the charge lasts 23 min. The illustrated model can be used to generate useful information in the design of a charging infrastructure scenario for BEVs, as the results may have implications on the demand side and especially on the investor side, for which the economic analysis necessary in planning new charging stations or upgrading an existing one should not ignore elements such as those presented above. In fact, in order to configure the structure of the station and the consequent economic exposure, investors must consider elements such as the attractiveness of the station, the necessary levels of electricity, and the quality of the service provided. The simulations show a highly complex dynamic of the system, which depends on many factors, sometimes in a non-linear way.
In our analysis, we hypothesized a BEV penetration of approximately 6% of the current fleet of Italian private vehicles. In a scenario of a charging infrastructure that provides high LoS, this penetration translates into the need for stations of several hundreds of kW (estimated from 300 kW to 500 kW of electric power), which represents a medium-high electric load in the context of urban areas. The rapid development of electric mobility implies an increase in the load for the stations and consequent greater demand for the electric grid. The widespread use of BEVs will increase the demand for electric power inside and outside the urban context. Consequently, it is necessary for the grid operators to implement notable structural variations to reallocate their generation points and to modify their energy distribution policies, with an increase in electric generation, integration of renewable energy sources, use of stationary storage, and a reconsideration of the distribution grid in the cities. These findings depend however on the sample chosen. With data on actual electric mobility, more realistic results can be obtained to be used in investor or policymaker analysis. Moreover, to better estimate the load demand, one should take into account the present EV market distributions and the related batteries characteristics. These points will be addressed in future work.