Index Terms—: Brain networks, functional magnetic resonance imaging, information theory, optimal control, optimal mass transport (OMT), Riemannian metric, system identification

A. On Optimal Mass Transport

Letp₀ (x) andp₁ (x) denote two probability density functions on${ℝ}^n$. The Wasserstein-2 metric between the two, denoted by w₂ (p₀,p₁), is defined by the square root of the optimal value as

wherem (x,y) represents a probability density function on the joint space${ℝ}^n\times {ℝ}^n$ with the marginals specified byp₀ andp₁ [17], [18]. A fluid-mechanics interpretation of w₂ (p₀, p₁)² was introduced in [21] and [22], which provided a Riemannian structure of the manifold of probability density functions. To introduce this formula, we consider the following continuity equation:

wherev_t (x) represents a time-varying velocity field andp_t (x) represents the time-varying density function. Then, w₂ (p₀, p₁)² is equal to (see [22])

The optimal solution ofp_t (x) is the geodesic on the manifold of probability density functions that connects the endpointsp₀ (x) andp₁ (x).

B. Hellinger–Bures Metric

In the special case whenp₀ (x) andp₁ (x) are zero-mean Gaussian probability density functions with

the corresponding geodesicp_t (x) at any fixed timet is also zero-mean Gaussian with the corresponding covariance matrix given by$P_t^{\text{omt}}$ [19], [20]. The geodesic distance is equal to the Wasserstein-2 metric w₂ (p₀, p₁), which also induces the following distance measure on the covariance matrices:

where ∥·∥_F denotes the Frobenius norm of a matrix.

The covariance path$P_t^{\text{omt}}$ is also equal to the geodesic induced by the Hellinger–Bures metric from quantum mechanics [23], [24] on the manifold of positive-definite matrices. In particular, let Δ ∈ Symⁿ, which represents a tangent vector at$P\in {\text{Sym}}_{++}^n$. Then, the Hellinger–Bures metric takes the form

whereM is the unique matrix in Symⁿ, which satisfies$\frac{1}{2}(PM+MP)=\Delta$ (see, e.g., [25] and [26]). The trajectory$P_t^{\text{omt}}$ in (2) is the shortest path connectingP₀ andP₁. Thus, it satisfies

with a given pair of endpoints$P_0,P_1\in {\text{Sym}}_{++}^n$ (see, e.g., [23]). The geodesic distance, also known as the Bures distance, is equal to${\text{d}}_{{\text{w}}_2}\left(P_0,P_1\right)$.

The Hellinger–Bures metric was originally proposed in quantum mechanics to compare density matrices, which are positive-definite matrices whose traces are equal to one. The density matrices are noncommutative analogues of probability vectors. In the commutative case whenP is restricted to be a diagonal matrix whose diagonal entries consist of a probability vector, then the Hellinger–Bures metric g_P,Bures(Δ) is equal the Fisher information metric, which will be discussed in Section III in detail. A generalization of the Hellinger–Bures metric was also introduced in [27] for multivariate spectral densities.

C. Hellinger–Bures-Based Linear Systems

Here, we present an alternative expression of (9) using the linear system in (3). For this purpose, we define that

which is equal to${\mathcal{E}}_{p_t}\left({\Vert {\dot{x}}_t\Vert }^2\right)$ ifẋ_t is given by (3). The following theorem relates the geodesics$P_t^{\text{omt}}$ to a linear system.

Theorem 1: Given$P_0,P_1\in {\text{Sym}}_{++}^n$. Let f_P be defined as above. Let$P_t^{\text{omt}}$,$A_t^{\text{omt}}$ be a pair of minimizer of

Then,$P_t^{\text{omt}}$ is equal to (2) and

where$Q=I-P_0^{-\frac{1}{2}}{(P_0^{\frac{1}{2}}{P}_1{P}_0^{\frac{1}{2}})}^{\frac{1}{2}}{P}_0^{-\frac{1}{2}}$. The optimal value of the objective function in (10) is equal to${\text{d}}_{{\text{w}}_2}{\left(P_0,P_1\right)}^2$.

Proof: The optimization problem (10) appears as a special case of (7) with the additional constraint that the velocity fieldv_t (x) =A_tx. Thus,${\text{d}}_{{\text{w}}_2}{\left(P_0,P_1\right)}^2$ is a lower bound of (10). Therefore, we need to show that$A_t^{\text{omt}}$ satisfies

To this end, we rewrite (2) as$P_t^{\text{omt}}=(I-tQ)P_0(I-tQ)$. Taking the derivative of$P_t^{\text{omt}}$ gives

which proves (12). Since all the eigenvalues ofQ are smaller than 1, the matrixtQ −I is invertible for allt ∈ [0, 1]. Therefore, (12) holds. Moreover, we have

which completes the proof. ■

We note that the matrix$A_t^{\text{omt}}$ is symmetric. Moreover, the matrices$A_{t_1}^{\text{omt}}$ and$A_{t_2}^{\text{omt}}$ commute for anyt₁, t₂. Therefore, the eigenspace ofA_t is fixed on the intervalt ∈ [0, 1].

The results from Theorem 1 can be further extended to obtain the optimal solutions corresponding to the following objective function:

whereẋ_t =A_tx,$W\in {\text{Sym}}_{++}^n$, and${\Vert {\dot{x}}_t\Vert }_W^2≔x^{\prime }Wx$. By applying change of variables, we define

Thus,${\text{f}}_{P_t}^W\left(A_t\right)=\text{tr}\left(A_{W,t}{P}_{W,t}{A}_{W,t}^{\prime }\right)={\text{f}}_{P_{W,t}}\left(A_{W,t}\right)$. Moreover, if (4) holds, then

Therefore, if$A_t^{\text{omt}}$ is the optimal system matrix that steersP_W,0 toP_W,1 with respect to f_P (A) given by Theorem 1, then$W^{-\frac{1}{2}}{A}_t^{\text{omt}}{W}^{\frac{1}{2}}$ is the optimal solution with respect to${\text{f}}_P^W(A)$. In the following section, we investigate a further extension of${\text{f}}_P^W(\cdot )$ by using a time-dependent weighting matrix, which provides a point of contact between OMT and information geometry.

A. Fisher–Rao Metric

For two probability density functionsp(x) and$\widehat{p}(x)$ on${ℝ}^n$, the Kullback–Leibler (KL) divergence

represents a well-established notion of distance between the two [28], [29]. If$\widehat{p}=p+\delta$ withδ representing a small perturbation, then the quadratic term of the Taylor’s expansion of d_KL(p∥p +δ) in terms ofδ is the Fisher information metric

For a probability distributionp(x,θ) parameterized by a vector θ, the corresponding metric is referred to as the Fisher–Rao metric and is given by

Ifp(x) is a zero-mean Gaussian probability density function parameterized by a covariance matrixP, then the metric becomes

Given$P_0,P_1\in {\text{Sym}}_{++}^n$, the geodesic in (1) is equal to the solution

which is the shortest path connectingP₀ andP₁. The corresponding path length is equal to (see, e.g., [30, Th. 6.1.6])

B. Fisher–Rao-Metric-Based Linear Systems

Following (5), we will introduce some positive quadratic forms f(A) so that the corresponding optimal covariance path is equal to$P_t^{\text{info}}$. One choice for the quadratic form would be given by

which satisfies that${\text{f}}_P^{\text{info},1}(A)\ge 0$,$\forall A\in {ℝ}^{n\times n}$ and$\text{P}\in {\text{Sym}}_{++}^n$.

Theorem 2: Given$P_0,P_1\in {\text{Sym}}_{++}^n$. Let${\text{f}}_P^{\text{info},1}(\cdot )$ be defined as above. Define$P_t^{\text{info}}$,A^info as a pair of minimizer of

Then, the optimal solution$P_t^{\text{info}}$ is equal to (1), and

Moreover, the optimal value of the objective function is equal to d_info(P₀, P₁)².

Proof: We rewrite (1) as

where the last equation is obtained using$e^{XY{X}^{-1}}=X{e}^Y{X}^{-1}$. Take the derivative of$P_t^{\text{info}}$ to give that

which completes the proof. ■

Note that the metric d_info(·) is invariant with respect to congruence transforms, i.e., d_info(P₀,P₁) = d_info(TP₀T′,TP₁T′) for any invertible matrixT. If$A_t^{\text{info}}$ is the optimal solution of (20), thenTA^infoT⁻¹ is the optimal solution corresponding to the pairTP₀T′,TP₁T′.

C. Weighted-Mass-Transport View

Since the map fromA_t to${\dot{P}}_t$ in (4) is injective, the quadratic forms ofA_t that lead to the geodesic$P_t^{\text{info}}$ are not unique. Here, we provide an alternative quadratic form, which provides an interesting relation between OMT and the Fisher–Rao metric. For this purpose, we define

which is a special form of (13) with the weighting matrixW =P_t and scaled by the factor of 4. The following lemma draws the relation between${\text{f}}_P^{\text{info},2}(A)$ and${\text{f}}_P^{\text{info},1}(A)$.

Lemma 1: Consider${\text{f}}_P^{\text{info},1}(\cdot )$ and${\text{f}}_P^{\text{info},2}(\cdot )$ be defined in (19) and (23), respectively. Then, we have

Proof: Take the difference

which is nonnegative. ■

We note that if$P^{-\frac{1}{2}}A{P}^{\frac{1}{2}}$ is symmetric, then${\text{f}}_P^{\text{info},1}(A)$ is equal to${\text{f}}_P^{\text{info},2}(A)$. This gives rise to the following proposition in parallel to Theorem 2.

Proposition 1: Given$P_0,P_1\in {\text{Sym}}_{++}^n$. Then$P_t^{\text{info}}$ andA^info given by (1) and (21), respectively, are the unique pair of minimizer of

with${\text{f}}_{P_t}^{\text{info},2}\left(A_t\right)$ defined as in (23).

Proof: From Lemma 1, we have

for any feasible pairs ofP_t andA_t. It is straightforward to verify that the above inequalities become equalities with the given$P_t^{\text{info}}$ andA^info. Therefore, the proposition holds. ■

D. Fluid-Mechanics Interpretation

Note that${\text{f}}_{P_t}^{\text{info,}2}\left(A_t\right)$ is a special case of${\text{f}}_{P_t}^W(A)$ in (13) when$W=4P_t^{-1}$. It is equal to

with the velocity field given byv_t(x) =A_tx. Thus, if the initial distributionp₀(x) is Gaussian, so isp_t(x), ∀t ≥ 0. Proposition 1 implies that among all the trajectories that connect two Gaussian probability density functionsp₀ andp₁, the lowest weighted-mass-transport cost is obtained by Gaussian density functions whose covariance matrices are equal to$P_t^{\text{info}}$. But this optimal trajectory is obtained under the linear constraint of velocity fields. Next, we remove this constraint and show that this trajectory is still optimal.

Theorem 3: Given two zero-mean Gaussian probability density functionsp₀ (x),p₁ (x) in${ℝ}^n$ with covariance matrices$P_0,P_1\in {\text{Sym}}_{++}^n$, respectively, define$p_t^{\text{info}}(x)$,$v_t^{\text{info}}(x)$ as the minimizer of

Then,$p_t^{\text{info}}(x)$ is zero-mean Gaussian whose covariance matrix is equal to$P_t^{\text{info}}$ in (1) and$v_t^{\text{info}}(x)=A^{\text{info}}x$ almost surely withA^info given by (21). Moreover, the optimal value is equal to d_info(P₀, P₁)².

Proof: First, we define

Then, applying integral by parts, we obtain

Therefore, the following optimization problem

provides a lower bound of (26) because the higher order moments of the probability density functions are not considered. On the other hand, (24) provides an upper bound of (26) because the velocity field is constrained to satisfy the linear system. Then, we show that the two bounds coincide.

Note that$V_t-C_t{P}_t^{-1}{C}_t^{\prime }\in {\text{Sym}}_{+}^n$. Thus, the optimalV_t of (29) should satisfy that$V_t=C_t{P}_t^{-1}{C}_t^{\prime }$. Therefore, (29) is equal to

Note that the constraint$P_t\in {\text{Sym}}_{+}^n$ is automatically satisfied due to the inverse barrier objective function. Therefore, we drop the constraint that$P_t\in {\text{Sym}}_{+}^n$ in the following analysis.

The optimization problem (30) is viewed as an optimal control problem, withC_t being matrix-valued control. Then, we derive the optimal solution using Pontryagin’s minimum principle (see, e.g., [31, Ch. 2]). A necessary condition for the optimal solution is that it must annihilate the variation of the Hamiltonian

with respect to the controlC_t. Here, Π_t is a symmetric matrix representing the Lagrange multiplier, i.e., the costate. By setting the partial derivative ofh₁ (·) with respect toC_t to zero, we obtain that

which provides a necessary condition that the optimalC_t is symmetric. Therefore,$C_t=\frac{1}{2}{\dot{P}}_t$ and the objective function in (29) becomes$\text{tr}\left(P_t^{-1}{\dot{P}}_t{P}_t^{-1}{\dot{P}}_t\right)={\text{g}}_{P,\text{Rao}}\left({\dot{P}}_t\right)$. Thus, the theorem directly follows from (17). For completeness, we finish the proof based on the Hamiltonian in the following.

The optimal${\dot{\Pi }}_t$ must annihilate the partial derivative ofh₁ (·) with respect toP_t. This gives rise to

Then, substituting (31) into (28) and (32), we obtain

Note that${\dot{P}}_t{\Pi }_t+P_t{\dot{\Pi }}_t=0$ for allt. Hence,P_tΠ_t is constant. We set

Thus, (31) is equal toC_t =AP_t =P_tA′. Multiplying both sides by$P_t^{-\frac{1}{2}}$ gives that$P_t^{-\frac{1}{2}}{C}_t{P}_t^{-\frac{1}{2}}=P_t^{-\frac{1}{2}}A{P}_t^{\frac{1}{2}}$, which is symmetric for allt. Substituting (34) into (33) gives

Therefore, we have

Multiplying$P_0^{-\frac{1}{2}}$ to both sides gives

By settingt = 1, we solve that

which is equal toA^info. Furthermore, from (22), the correspondingP_t is equal to$P_t^{\text{info}}$. Then, the optimal covariance matrix in (27) is singular and has rankn. Thus, the optimal velocity fieldv_t(x) is equal toA^infox almost surely, implying that the correspondingp_t(x) is Gaussian. Therefore, the theorem is proved. ■

Note that the system matrixA^info is constant. Thus, bothA^info and$A_t^{\text{omt}}$ have fixed eigenspaces. Next, we introduce a different quadratic form ofA_t, which leads to system matrices with rotating eigenspaces.

A. Weighted-Least-Squares Cost Functions

Note that ifX is an antisymmetric matrix, i.e.,X′ = −X, thene^Xt is a rotation matrix. Consequently, if the system matrixA is antisymmetric, then the state covariance matrix has rotating eigenspaces. In this regard, we decompose

where$A_{\text{s}}≔\frac{1}{2}\left(A+A^{\prime }\right),A_{\text{a}}≔\frac{1}{2}\left(A-A^{\prime }\right)$ are the symmetric and antisymmetric parts ofA, respectively. Then, we define the following weighted-least-squares (WLS) function:

where the scalarϵ > 0 scalar weighs the relative significance of symmetric and asymmetric parts ofA. IfA satisfies$\dot{P}=AP+P{A}^{\prime }$ for a given pair$\dot{P}$ andP, then fϵ (A) is considered as a quadratic form of the noncommutative division of$\dot{P}$ byP, similar to the Fisher–Rao metric. Actually, for scalar-valued covariances, fϵ (A) is equal to the Fisher–Rao metric.

Following (5), we consider the optimal solution to

for a given pair of endpoints$P_0,P_1\in {\text{Sym}}_{++}^n$ and a scalarϵ > 0.

B. Optimal Covariance Paths

To introduce the solution to (36), we define

The next lemma shows thatT_ϵ,t(·) is equal to the state transition matrix of a linear time-varying system.

Lemma 2: Given$A\in {ℝ}^{n\times n}$ and a scalarϵ > 0, define

Then,${\dot{T}}_{ϵ,t}(A)=A_{ϵ,t}{T}_{ϵ,t}(A)$.

Proof:

■

The following corollary is a direct result of Lemma 2.

Corollary 1: Given$P_0\in {\text{Sym}}_{++}^n$,$A\in {ℝ}^{n\times n}$ and a scalarϵ > 0, define

Then, the following equation holds:

The solution to (36) is presented as follows.

Theorem 4: Given$P_0,P_1\in {\text{Sym}}_{++}^n$ and a scalarϵ > 0, if there exists a matrix Π₀ ∈ Symⁿ such that the covariance path

withT_ϵ,t(·) given by (37) and

satisfies$P_{ϵ,1}^{\text{wls}}=P_1$ att = 1, then$P_{ϵ,t}^{\text{wls}}$ is a minimizer of (36). The corresponding optimalA_t is equal to

Proof: Consider (36) as an optimal control problem, withA_t being matrix-valued control. Then, the Hamiltonian is as follows:

It is necessary that${\dot{\Pi }}_t$ annihilates the partial derivative ofh₂ (·) with respect toP_t, which gives rise to

Moreover, the partial derivative ofh₂ (·) with respect to the controlA_t vanishes, which leads to

SolvingA_t from (43), we obtain

Then, substituting (44) into (4) and (42), respectively, we obtain

Next, it can be verified that$\left(\stackrel{\cdot }{{\Pi }_t{P}_t}\right)-\left(\stackrel{\cdot }{P_t{\Pi }_t}\right)=0$. Thus, the asymmetric part ofA_t, which is equal to

is constant and denoted byA_a. Taking the derivative of its symmetric part

gives that

SinceA_a is constant, the solution to the above equation is equal to

Therefore, the optimalA_t has the form

withA =A_a +A_s being the initial value ofA_t. Next, we define a new variable

whose derivative is equal to

Thus, the solution to${\widehat{P}}_t$ is equal to

Substituting this solution to (45), we obtain that the optimalP_t has the form

In a similar way, we define

Then, we have

whose solution is equal to

If (40) holds, then$A_{\text{s}}+ϵA_{\text{a}}^{\prime }=-P_0{\Pi }_0$. It is straightforward to show that (44) holds for allt > 0 for the provided expressions forA_t, P_t, Π_t. In this case,${\text{f}}_{ϵ}\left(A_t\right)={\Vert A_{\text{s}}\Vert }_{\text{F}}^2+ϵ{\Vert A_{\text{a}}\Vert }_{\text{F}}^2$ for allt. Therefore, ifA is equal to$\widehat{A}$ in (40), then the proposed trajectories$A_{ϵ,t}^{\text{wls}}$ and$P_{ϵ,t}^{\text{wls}}$ are local minimizer, which completes the proof. ■

It is interesting to note that the geodesics$P_t^{\text{info}}$ are special cases of the covariance paths in (39) whenϵ = −1. The existence and uniqueness of the covariance paths that satisfy the conditions in Theorem 4 forϵ > 0 will be discussed somewhere else [32].

Next, we provide an upper bound of the optimal value of (36) for allϵ > 0. For this purpose, we define

which is symmetric. The corresponding covariance path

is a feasible solution to (36). Therefore, the following proposition holds.

Proposition 2: Given$P_0,P_1\in {\text{Sym}}_{++}^n$, the optimal value of (36) for anyϵ > 0 is not larger than${\Vert \text{log}(P_0^{-\frac{1}{2}}{(P_0^{\frac{1}{2}}{P}_1{P}_0^{\frac{1}{2}})}^{\frac{1}{2}}{P}_0^{-\frac{1}{2}})\Vert }_{\text{F}}^2$.

Ifϵ → ∞, then the symmetric matrix$\widehat{A}$ becomes the optimal solution. Moreover, ifP₀ andP₁ commute, thenP_t in (46) is equal to$P_t^{\text{info}}$.

A. Interpolating Covariance Matrices

In this example, we highlight the difference between$P_{ϵ,t}^{\text{wls}}$ and the other two types of trajectories, i.e.,$P_t^{\text{omt}}$ and$P_t^{\text{info}}$, using the following two matrices as the endpoints:

Applying (1) and (2), we obtain

which are all diagonal. On the other hand, ifϵ = 0, there are infinitely many antisymmetric matrices$A_0^{\text{wls}}$ of the following form:

withk being an arbitrary integer. All these matrices are able to steerP₀ toP₁. Therefore, the corresponding objective function${\text{f}}_0\left(A_0^{\text{wls}}\right)$ is equal to zero. The covariance path$P_{0,t}^{\text{wls}}$ is equal to