1 Introduction↩︎

Complex networks have become a fundamental tool for modeling and understanding the structure and dynamics of interconnected systems across various domains, such as biology [1], social sciences, and information technology [2]. In many cases, networks are considered static entities, where the set of nodes and edges remains unchanged over time. However, numerous real-world networks evolve dynamically, with both topology and edge weights varying over time. Such networks are referred to as time-evolving or temporal networks [3]. These dynamic interactions capture a wide range of phenomena, including communication patterns in social networks, where two users may follow each other for a certain period before one of them stops following the other. Similarly, in transportation systems, flight routes between airports may only be active during specific months of the year [4]. In brain networks, neural pathways may be activated only under particular conditions [5], [6]. Understanding how these networks evolve is crucial for uncovering patterns, predicting future behavior, and controlling processes occurring within the network.

The study of time-evolving networks has led to the development of an array of methods for analyzing how information propagates through such systems. Concepts such as dynamic centrality measures [7]–[10], temporal motifs [11], and communicability matrices [12] extend classical static graph metrics [13] to the temporal domain. These tools enable a more accurate and comprehensive understanding of the roles of individual nodes and overall network behavior in time-varying structures. In this context, we focus on how to control network connections by adding, modifying, or removing edges to drive measures of network communicability toward a desired state. This work addresses the centrality problem in dynamic networks, specifically dynamic centrality [7]. This measure computes the centrality of a node by solving a matrix differential equation when computing communicability indices between any two nodes or a vector differential equation when determining dynamic broadcast or receive centrality. Our goal is to determine a minimal modification of the time-dependent network that achieves a prescribed centrality at a given final time \(T\). This generalizes the analogous problem for centralities in static networks [14]–[16]. Since dynamic centrality is formulated as the solution of an ordinary differential equation (ODE), the problem can be posed as an optimal control problem for an ODE, either matrix- or vector-based, where the control acts on network edges to steer the centrality toward the desired values.

Control theory has been widely applied to networks in various contexts. One example is the introduction of a dynamical process defined on network edges [17], demonstrating that the controllability properties of such edge-based dynamics differ significantly from those of traditional node-based dynamics. In the latter case, it is possible to identify a set of driver nodes with time-dependent control inputs that can steer the entire system’s dynamics [18]. This edge-based control framework is particularly relevant for monitoring the evolution of dynamic phenomena on networks, including epidemic spread, opinion dynamics and other compartmental models [19]–[21]. Recent work has further addressed controllability in dynamic networks by proposing an adaptive control algorithm that minimizes driver node switching costs under stochastic topology changes without requiring prior knowledge of network evolution [22]; while in the same direction also a graphical criterion to evaluate the controllable subspace of temporal networks via maximum matching on aggregated static representations has been introduced [23], and another approach adopts a temporal ant colony optimization–based approach to efficiently approximate maximum matching in temporal networks to identify driver nodes for structural controllability [24].

For the problems we consider here, we have either a nonlinear or a linear ODE whose state variable represents the time-dependent centrality measure. Hence, the problem we face is that of steering this state towards a prescribed value by modifying the time-dependent edge weights. This is indeed an instance of a bilinear control problem in which the control, or more precisely one of its functions, acts multiplicatively on the state of the system. As it is often the case, the control problem is highly dimensional, and due to the underlying state equation [7] involves the computation of large-scale matrix-function [25], and matrix-function vector products [26], [27]. We discuss and adapt the related techniques for working in our regime of interest—large scale, repeated computation, possibly non-normal matrices. To formulate the solution of the optimal control problem we make use of the Pontryagin Maximum Principle (PMP) which requires computing the evolution of the adjoint variable, due to the presence of the matrix function in the state equation this requires the computation of large-scale Fréchet derivatives of a matrix function [28], and again their application against a vector [29] through the construction of suitable rational Krylov methods. Henceforth, the main contributions of this work include a pattern-constrained optimal-control formulation for steering dynamic receive centrality in temporal networks and the derivation of sufficient linear feasibility constraints to preserve the well-posedness of the logarithmic model. Furthermore, we establish PMP-based first-order optimality conditions for both logarithmic and linear dynamics and present a matrix-free projected-gradient algorithm that utilizes Krylov and rational Krylov approximations for matrix logarithm and Fréchet derivative actions. Finally, we provide empirical validation of our approach using both synthetic and real temporal networks. The manuscript is organized as follows. In Section 1.1, we introduce the notation and context necessary to describe time-dependent networks. Section 2 presents the construction of the dynamic centrality measure we focus on [7], which serves as the foundation for formulating the problem of steering centrality toward a desired state in Section 3. Once the problem is formulated, we discuss its solution using the PMP in Section 4, followed by a description of its discretization and efficient implementation in Section 5. Section 6 applies the proposed methodology to synthetic and real-world temporal network test cases. Finally, in Section 7, we conclude our analysis and discuss potential extensions.

1.1 Notation and basic definitions↩︎

The edge set \(E_t\) and weight function \(w_t\) are both time-dependent, reflecting the fact that the structure and properties of the graph may change as time progresses. To represent it in a compact form we define the related concept of a time-dependent (weighted) adjacency matrix.

Throughout the paper we assume that all \(G_t\) are such that \(A(t) \leq \overline{A}\) for all \(t \in \mathbb{T}\). This boundedness property guarantees that the spectral radius of all \(A(t)\) is bounded as well, and will permit formulating the differential models in Section 2. We also note that, in many applications, the network \(G_t\) is observed only at a finite collection of time points \(0 = t_1 < t_2 < \cdots < t_k = T\). In this setting, we still take the time domain to be \(\mathbb{T} = [0,T]\) and extend the observed network to all \(t \in \mathbb{T}\) by defining \(G_t\) to be piecewise constant on the intervals \([t_i, t_{i+1})\), that is, \(G_t = G_{t_i}\) for \(t \in [t_i, t_{i+1})\).

The other set of tools we need to write the problem formulation and solution are matrix functions and some related properties.

This can be computed for a \((2n-1)\)-times differentiable function \(f\) in an open set containing the eigenvalues of \(A\) as the matrix function of the \(2 \times 2\) block matrix [25] \[\label{eq:explicit-frechet} \mathrm{L}_f(A, E) = \left[ \begin{matrix} I & O \end{matrix} \right] f\left( \begin{bmatrix} A & E \\ O & A \end{bmatrix} \right) \left[ \begin{matrix} O \\ I \end{matrix} \right],\tag{1}\] where \(I\) and \(O\) are the \(n \times n\) identity and zero matrices, respectively. For some particular functions \(f\) it is possible to find other expressions of \(\mathrm{L}_f(A,E)\), e.g., in the case of the logarithm [28].

Finally, we recall the definition of Katz’s centrality index [30] for a static network, which we will need for using the collection of Katz centralities computed for the different time steps \(t \in \mathbb{T}\).

2 The dynamical system formulation of the network centrality↩︎

With the goal of quantifying the ability of node \(v_i\) to communicate with node \(v_j\) over time in a dynamic network, in [7] a dynamical system formulation is introduced for defining a dynamic communicability matrix for the network. The foundation of this approach is the time-dependent adjacency matrix \(A(t)\) (Definition 2), where edges appear, disappear and, generally, modify their weights as the network evolves. The dynamic communicability matrix \(S(t)\) is constructed by counting dynamic walks that involve edges present at time \(t\) or earlier, while also accounting for the temporal decay of information relevance.

To update \(S(t)\) over small time intervals \(\delta t\), the following formula is introduced in [7]: \[\label{eq:discrete-formulation} S(t + \delta t) = (I + e^{-b \delta t} S(t)) (I - a A(t + \delta t))^{-\delta t} - I,\tag{2}\] where \(S(0) = 0\), \(I\) is the identity matrix, \(a > 0\) is the edge attenuation parameter that controls the weight assigned to walks based on the number of edges, and \(b > 0\) is the temporal downweighting parameter that governs how older walks are penalized. This formulation allows the matrix \(S(t)\) to count walks between nodes while diminishing the importance of longer and older walks.

Taking the limit as \(\delta t \to 0\) in 2 , the following ordinary differential equation (ODE) for the communicability matrix is obtained: \[\label{eq:ODE-formulation} \dot{C}(t) = -b (C(t) - I) - C(t) \log (I - a A(t)),\tag{3}\] where \(C(t) = I + S(t)\) and \(\log\) denotes the matrix logarithm. In this ODE, the term \(-b (C(t) - I)\) ensures that older walks are downweighted over time, while the term \(-C(t) \log (I - a A(t))\) accounts for the dynamic walks utilizing edges present at time \(t\).

The matrix \(C(t)\) encodes the communication capability at time \(t\) between node \(i\) and node \(j\) through dynamic walks in the network, applying downweighting to both the length of walks and the age of information.

To derive a more efficient collective measure of the importance of the nodes in the time-evolving network from the communicability matrix \(C(t)\), as in [7], we define the dynamic broadcast centrality and dynamic receive centrality vectors as: \[\label{eq:broadcast95and95receiving} \mathbf{b}(t) = C(t) \mathbf{1} \quad \text{and} \quad \mathbf{r}(t) = C(t)^\top \mathbf{1},\tag{4}\] where \(\mathbf{1}\) is the vector of ones. These measures capture the propensity of each node to broadcast and receive information, respectively.

2.1 Extensions and generalizations↩︎

All these formulations can be generalized to the case \[\label{eq:generalized} \dot{C}(t) = - b( C(t) - I ) + C(t) \log(H(A(t))), \qquad C(0) = I,\tag{5}\] where \(H(A(t))\) is a given matrix function, i.e., we consider a given combination of walk lengths in the model. A choice that greatly simplifies the construction is \(H(x) = \exp(a x)\). With this choice the entire formulation of the problem is reduced to the linear case \[\label{eq:simplified} \dot{C}(t) = - b( C(t) - I ) + a C(t) A(t), \qquad C(0) = I,\tag{6}\] if \(\log(\exp(a A(t) )) = a A(t)\), which is true when \(|\Im(\lambda_i( a A(t)))| < \pi\), for all \(i = 1,\ldots,n\) and \(t \in \mathbb{T}\), or, for any consistent matrix norm \(\|a A(t)\| < \pi\) for all \(t \in \mathbb{T}\); see [25]. In the case in which \(A(t) = A(t)^\top\), i.e., the network is undirected this always holds. Otherwise, we need to select \(a\) in such a way that \(|\Im(\lambda_i(a A(t)))| < \pi\) for all \(t \in \mathbb{T}\), \(i=1,\ldots,n\).

3 Steering the dynamic centralities to prescribed values↩︎

Suppose we now have access to a network that evolves over time, \(G_t\), and that we can selectively intervene at each time \(t\) on the weights of all or some of the edges. Our objective is to steer the network’s communicability, represented in matrix form, or its broadcast and reception centrality indices, towards a desired state at a given time \(T\). We can formulate these problems as optimal control problems for the ODEs 3 or 6 , or for the derived broadcast and receiving centralities 4 .

For the communicability matrix, given a desired state \(C^* \in \mathbb{R}^{n\times n}\) at time \(T\) and a set \(\mathcal{P} \subseteq V \times V\), we can express the problem as: \[\label{eq:matrix-formulation} \begin{array}{rl} \displaystyle \min_{U} & \displaystyle \mathrm{J}(C,U) = \frac{1}{2}\| C(T) - C^* \|_F^2 + \frac{\alpha}{2}\int_{0}^{T} \| U(t) \|_F^2 \, \mathrm{d}t, \qquad \alpha > 0\\ \text{s.t.} & \begin{cases}\dot{C}(t) = -b (C(t) - I) - C(t) \log \left(I - a (A(t) + U(t)) \right), \\ C(0) = I, \end{cases} \\ & A(t) + U(t) \geq 0, \\ & \rho(A(t) + U(t)) < \frac{1}{a}, \\ & [U(t)]_{i,j} = 0 \text{ if } (i,j) \not\in \mathcal{P}, \end{array}\tag{7}\] where we require that the control \(U(t)\) be a matrix of the same size as \(A(t)\), maintaining a sparsity pattern contained in the predetermined set \(\mathcal{P}\), that the weights of the perturbed network remain nonnegative, and that the uniqueness condition for the matrix logarithm remains satisfied along the solution trajectory as described in Remark 2. The \(\alpha \in \mathbb{R}_{> 0}\) value in the cost functional is a Tikhonov regularization parameter of the control.

A similar formulation can be applied to the dynamic receive case, where we selected a desired receive centrality vector \(\mathbf{r}^* \geq \boldsymbol{1} \in \mathbb{R}^n\): \[\label{eq:dynamic-receive-original} \begin{array}{rl} \displaystyle \min_{U} & \displaystyle \mathrm{J}(\mathbf{r},U) = \frac{1}{2}\| \mathbf{r}(T) - \mathbf{r}^* \|_{2}^2 + \frac{\alpha}{2}\int_{0}^{T} \| U(t) \|_F^2 \, \mathrm{d}t, \qquad \alpha > 0 \\ \text{s.t.} & \begin{cases}\dot{\mathbf{r}}(t) = -b (\mathbf{r}(t) - \mathbf{1}) - \log \left(I - a (A(t) + U(t)) \right)^\top \mathbf{r}(t), \\ \mathbf{r}(0) = \mathbf{1}, \end{cases}\\ & A(t) + U(t) \geq 0, \\ & \rho(A(t) + U(t)) < \frac{1}{a}, \\ & [U(t)]_{i,j} = 0 \text{ if } (i,j) \not\in \mathcal{P}. \end{array}\tag{8}\] Observe that the broadcast vector \(\mathbf{b}(t)\) tracks the information that has flowed out of each node, thus it requires access to the entries of the \(C(t)\) matrix. Consequently, the receive centrality is fundamentally simpler to solve than the broadcast centrality.

If we consider the extended framework discussed in [7] with a generic matrix function \(H(A(t))\) and select \(H(A(t))=\exp(a A(t))\), we may employ the simplified linear model 6 , so as to lead us to the solution of the control problem \[\label{eq:dynamic-receive-simplified-original} \begin{array}{rl} \displaystyle \min_{U} & \displaystyle \mathrm{J}(\mathbf{r},U) = \frac{1}{2}\| \mathbf{r}(T) - \mathbf{r}^*\|_2^2 + \frac{\alpha}{2} \int_{0}^{T} \| U(t)\|_F^2 \mathrm{d}t, \qquad \alpha > 0 \\ \text{ s.t. }& \begin{cases}\dot{\mathbf{r}}(t) = - b (\mathbf{r}(t) - \mathbf{1}) + a (A(t) + U(t))^\top \mathbf{r}(t),\\ \mathbf{r}(0) = \mathbf{1}, \end{cases}\\ & A(t) + U(t) \geq 0, \\ & |\Im(\lambda_i(A(t) + U(t)))| < \frac{\pi}{a}, \quad \forall\, i = 1,\ldots,n \\ & [U(t)]_{i,j} = 0 \text{ if } (i,j) \not\in \mathcal{P}. \end{array}\tag{9}\] Observe that although we no longer compute any matrix logarithm in the state equation of 9 , to have the transformation between 5 and 6 we need the condition \(\log(\exp(A)) = A\), hence the condition on the imaginary part of the spectrum. Observe also that if both \(G_t\) is undirected and that \(U(t) = U(t)^\top\) for all \(t \in \mathbb{T}\), then such a condition can be removed; see the discussion in Section 2.1. In the case of directed graphs, a sufficient condition for the constraint on the imaginary part of the eigenvalues can be stated in terms of the outdegree, as in Proposition 2.

3.1 Relaxing the spectral constraints↩︎

The constraint on the spectral radius of \(A(t)+U(t)\) in 7 is nonlinear and nonconvex, hence its use complicates any optimization procedure one can use to solve the control problem. To obtain a more efficient procedure with a linear constraint, we consider the following relaxation which makes use of the collection of Katz centrality measures (Definition 5) for all \(t \in \mathbb{T}\), and establishes first the existence of a choice of the \(a\) parameter which makes Katz centrality well defined for a general time evolving graph (Definition 1).

According to Proposition 1, the linear inequality \((I - a(A(t)+U(t)))\boldsymbol{\mu}(t) \geq \mathbf{0}\), which is easier to handle computationally than the nonlinear spectral radius constraint, guarantees that \(\rho((A(t)+U(t))\le \frac{1}{a}\). To provide a sufficient condition guaranteeing that the strict inequality \(\rho(A(t)+U(t)) < \frac{1}{a}\) is satisfied, we then need to just add a safety guarantee of a small \(\frac{1}{a} \gg \epsilon>0\), asking that \(U(t) \boldsymbol{\mu}(t) \leq \left( \frac{1}{a} - \epsilon \right)\boldsymbol{1}\), where we used the fact that \((I-aA(t))\boldsymbol{\mu}(t) = \mathbf{1}\).

We need a relaxation also for the constraint of the imaginary part of the spectrum in problem 9 in the case of directed graphs. We can employ the observation on the consistent matrix norms made in Section 2.1.

As for the other relaxation, we need to add a safety guarantee of a small \(\frac{\pi}{a} \gg \varepsilon > 0\), asking that \[(A(t) + U(t)) \mathbf{1} \leq \left( \frac{\pi}{a} - \varepsilon \right) \mathbf{1},\] observe also that in network terms, this amounts to a constraint on the out-degree of the network. An immediate corollary of Proposition 2 is then obtained by changing from out-degree to in-degree and from the compatible norm \(\|\cdot\|_\infty\) to the compatible norm \(\|\cdot\|_1\).

We can therefore rewrite the two problems 8 , 9 using the simplified bounds, i.e., for the case of 8 , we rewrite it as \[\label{eq:dynamic-receive} \begin{array}{rl} \displaystyle \min_{U} & \displaystyle \mathrm{J}(\mathbf{r},U) = \frac{1}{2}\| \mathbf{r}(T) - \mathbf{r}^* \|_{2}^2 + \frac{\alpha}{2}\int_{0}^{T} \| U(t) \|_F^2 \, \mathrm{d}t, \qquad \alpha > 0, \; \epsilon > 0 \\ \text{s.t.} & \begin{cases}\dot{\mathbf{r}}(t) = -b (\mathbf{r}(t) - \mathbf{1}) - \log \left(I - a (A(t) + U(t)) \right)^\top \mathbf{r}(t), \\ \mathbf{r}(0) = \mathbf{1}, \end{cases}\\ & A(t) + U(t) \geq 0, \\ & U(t) \boldsymbol{\mu}(t) \leq \left( \frac{1}{a} - \epsilon \right)\boldsymbol{1}, \\ & [U(t)]_{i,j} = 0 \text{ if } (i,j) \not\in \mathcal{P}, \end{array}\tag{10}\] while we simplify the constraints in 9 for a directed graph as \[\label{eq:dynamic-receive-simplified} \begin{array}{rl} \displaystyle \min_{U} & \displaystyle \mathrm{J}(\mathbf{r},U) = \frac{1}{2}\| \mathbf{r}(T) - \mathbf{r}^*\|_2^2 + \frac{\alpha}{2} \int_{0}^{T} \| U(t)\|_F^2 \mathrm{d}t, \qquad \alpha > 0, \; \epsilon > 0 \\ \text{ s.t. }& \begin{cases}\dot{\mathbf{r}}(t) = - b (\mathbf{r}(t) - \mathbf{1}) + a (A(t) + U(t))^\top \mathbf{r}(t),\\ \mathbf{r}(0) = \mathbf{1}, \end{cases}\\ & A(t) + U(t) \geq 0, \\ & (A(t) + U(t)) \mathbf{1} \leq \left( \frac{\pi}{a} -\epsilon\right) \mathbf{1}, \\ & [U(t)]_{i,j} = 0 \text{ if } (i,j) \not\in \mathcal{P}, \end{array}\tag{11}\] while for an undirected graph, where \(A(t) = A^\top(t)\) for all times, we have \[\label{eq:dynamic-receive-simplified-undirected} \begin{array}{rl} \displaystyle \min_{U} & \displaystyle \mathrm{J}(\mathbf{r},U) = \frac{1}{2}\| \mathbf{r}(T) - \mathbf{r}^*\|_2^2 + \frac{\alpha}{2} \int_{0}^{T} \| U(t)\|_F^2 \mathrm{d}t, \qquad \alpha > 0 \\ \text{ s.t. }& \begin{cases}\dot{\mathbf{r}}(t) = - b (\mathbf{r}(t) - \mathbf{1}) + a (A(t) + U(t))^\top \mathbf{r}(t),\\ \mathbf{r}(0) = \mathbf{1}, \end{cases}\\ & A(t) + U(t) \geq 0, \\ & [U(t)]_{i,j} = 0 \text{ if } (i,j) \not\in \mathcal{P}, \\ & U(t) = U^\top(t). \end{array}\tag{12}\]

We can prove that the admissible set of controls in the three cases are not empty under reasonable choices of the modification pattern \(\mathcal{P}\).

We introduce the admissible control set for the linear model for the general case of directed graphs 11 : \[\label{eq:adm95set95directed} \mathcal{U}_{ad}(t) = \left\{ U(t) \in \mathbb{R}^{n \times n} \;\middle|\; \begin{align} & U(t) \ge -A(t), \\ & U(t) \mathbf{1} \leq \left( \frac{\pi}{a} - \varepsilon \right) \mathbf{1} - A(t)\mathbf{1}, \\ & [U(t)]_{i,j} = 0 \text{ if } (i,j) \notin \mathcal{P} \end{align} \right\},\tag{13}\] while for undirected graphs case 12 we have \[\label{eq:adm95set} \mathcal{U}_{ad}(t) = \left\{ U(t) \in \mathbb{R}^{n \times n} \;\middle|\; \begin{align} & U(t) \ge -A(t), \\ & U(t) = U^\top(t) , \\ & [U(t)]_{i,j} = 0 \text{ if } (i,j) \notin \mathcal{P} \end{align} \right\},\tag{14}\] which are not empty under a stricter assumption on the pattern of the admissible modifications.

Propositions 3 and 4 suggest that a reasonable choice for the modification pattern \(\mathcal{P}\) is therefore the collection of the patterns of the \(A(t)\) for all \(t \in \mathbb{T}\): from a modeling perspective, this choice means that we want to modify edges that existed at some point in the network’s temporal evolution.

3.2 Existence of the optimal control↩︎

To guarantee that the optimal control problems formulated in the previous sections are well-posed, we must establish the existence of an optimal control \(U^*(t)\) and its corresponding state trajectory. We state the existence result focusing on the simplified linear model 11 , utilizing the classic Filippov-Cesari existence theorem for Bolza problems [33].

4 Formulation of the solution via the PMP↩︎

The optimization problems described in 10 , 11 and 12 are instances of a general optimal control problem with constraints imposed on the control variable. To outline our solution approach, we first recall the general formulation of an optimal control problem.

Consider a system characterized by a state variable \(\mathbf{x}(t) \in \mathbb{R}^n\), which evolves according to the following system of differential equations: \[\begin{cases} \dot{\mathbf{x}}(t) = \mathbf{f}(\mathbf{x}(t), \mathbf{u}(t), t), \\ \mathbf{x}(0) = \mathbf{x}_0, \end{cases}\] where \(\mathbf{u}(t) \in \mathcal{U}_{ad} \subset \mathbb{R}^m\) is the control variable which is taken into a subset \(\mathcal{U}_{ad}\) of \(\mathbb{R}^m\) of admissible controls, and \(\mathbf{f} : \mathbb{R}^n \times \mathbb{R}^m \times \mathbb{R} \to \mathbb{R}^n\) is a continuous function defining the system dynamics. Our objective is to determine an optimal control function \(\mathbf{u}^*(t)\) that minimizes a given cost functional \(\mathrm{J}\), which typically takes the form: \[\mathrm{J}(\mathbf{x},\mathbf{u}) = \int_{0}^{T} L(\mathbf{x}(t), \mathbf{u}(t), t) \, \mathrm{d}t + \Phi(\mathbf{x}(T)),\] where \(L : \mathbb{R}^n \times \mathbb{R}^m \times \mathbb{R} \to \mathbb{R}\) represents the running cost, and \(\Phi : \mathbb{R}^n \to \mathbb{R}\) is the terminal cost. The optimal control function \(\mathbf{u}^*(t)\) must ensure that the state trajectory \(\mathbf{x}(t)\) satisfies the given dynamic constraints and initial conditions while minimizing the cost functional.

For a.e. \(t\in[0,T]\), we work with effective admissible control sets that are nonempty, closed, convex, and compact in all the cases under consideration. In addition, the corresponding state dynamics are measurable with respect to \(t\) and continuously differentiable with respect to the variables \((\mathbf{r},U)\) on the admissible set. Hence, the PMP applies and provides first-order necessary optimality conditions; see, e.g., [34]. It provides necessary conditions that an optimal control \(\mathbf{u}^*(t)\) and corresponding state \(\mathbf{x}^*(t)\) must satisfy. Define the Hamiltonian function \(\mathrm{H} : \mathbb{R}^n \times \mathbb{R}^m \times \mathbb{R}^n \times \mathbb{R} \to \mathbb{R}\) as: \[\label{eq:generic95hamiltonian} \mathrm{H}(\mathbf{x}, \mathbf{u}, \boldsymbol{\lambda}, t) = \boldsymbol{\lambda}(t)^\top \mathbf{f}(\mathbf{x}, \mathbf{u}, t) + L(\mathbf{x}, \mathbf{u}, t),\tag{15}\] where \(\boldsymbol{\lambda}(t) \in \mathbb{R}^n\) is the costate (or adjoint) variable.

The PMP asserts that for an optimal control \(\mathbf{u}^*(t)\) and state trajectory \(\mathbf{x}^*(t)\), there exists a nontrivial costate vector \(\boldsymbol{\lambda}(t)\) such that:

1. The state and costate dynamics satisfy: \[\begin{align} & \dot{\mathbf{x}}^*(t) = & \frac{\partial \mathrm{H}}{\partial \boldsymbol{\lambda}}(\mathbf{x}^*(t), \mathbf{u}^*(t), \boldsymbol{\lambda}(t), t), \tag{16}\\ & \dot{\boldsymbol{\lambda}}(t) = & -\frac{\partial \mathrm{H}}{\partial \mathbf{x}}(\mathbf{x}^*(t), \mathbf{u}^*(t), \boldsymbol{\lambda}(t), t). \tag{17} \end{align}\]

2. The optimal control \(\mathbf{u}^*(t)\) minimizes the Hamiltonian: \[\label{eq:hamiltonian95minimization} \mathbf{u}^*(t) = \arg \min_{\mathbf{u} \in \mathcal{U}_{ad}} \mathrm{H}(\mathbf{x}^*(t), \mathbf{u}, \boldsymbol{\lambda}(t), t),\tag{18}\] for the set of admissible controls \(\mathcal{U}_{ad} \subset \mathbb{R}^m\), for which the unconstrained gradient is given by \[\label{eq:hamiltonian95gradient} \mathbf{g}(t) = \frac{\partial \mathrm{H}}{\partial \mathbf{u}}.\tag{19}\]

3. The boundary conditions include the initial condition \(\mathbf{x}(0) = \mathbf{x}_0\), and if the terminal state \(\mathbf{x}(T)\) is free, the transversality condition for the costate: \[\label{eq:terminal95condition} \boldsymbol{\lambda}(T) = \frac{\partial \Phi}{\partial \mathbf{x}}(\mathbf{x}(T)).\tag{20}\]

In the following, we reformulate the problems 10 , 11 and 12 in a manner that allows us to leverage the PMP to derive a solution procedure. Specifically, we express these conditions as a constrained optimization problem: \[\min_{\mathbf{u} \in \mathcal{U}_{ad}} \mathrm{J}(\mathbf{u}) = \min_{\mathbf{u} \in \mathcal{U}_{ad}} \mathrm{J}(\mathbf{x}(\mathbf{u}), \mathbf{u}) = \min_{\mathbf{x}, \mathbf{u}} \left\{ \mathrm{J}(\mathbf{x}, \mathbf{u}) \mid e(\mathbf{x}, \mathbf{u}) = 0 \right\},\] where the constraint function \(e(\mathbf{x}, \mathbf{u})\) is given by \[e(\mathbf{x}, \mathbf{u}) = \begin{bmatrix} \dot{\mathbf{x}}(t) - \mathbf{f}(\mathbf{x}, \mathbf{u}, t) \\ \mathbf{x}(0) - \mathbf{x}_0 \end{bmatrix}.\] The unconstrained gradient of the objective function is denoted by \(\mathbf{g}(\mathbf{u})\). Suppose we have a projection operator \(\Pi_{\mathcal{U}_{ad}}\) that projects a vector from \(\mathbb{R}^m\) onto the feasible set \(\mathcal{U}_{ad}\), i.e., \[\Pi_{\mathcal{U}_{ad}}(\mathbf{y}) = \arg\min_{\mathbf{x} \in \mathcal{U}_{ad}} \| \mathbf{x} - \mathbf{y}\|_2.\] We can then construct a projected gradient method, incorporating Nesterov’s acceleration to optimize this problem (Algorithm [alg:nesterov95pg]).

The key idea is to maintain, alongside the iterates \(\mathbf{u}^{(k)}\), an extrapolated auxiliary sequence \(\mathbf{y}^{(k)}\) that incorporates a momentum term from the previous step: \[\mathbf{y}^{(k)} = \mathbf{u}^{(k)} + \beta_k \left( \mathbf{u}^{(k)} - \mathbf{u}^{(k-1)} \right), \qquad \beta_k = \frac{k-1}{k+2}.\] The gradient of the Hamiltonian and the state/costate equations are then evaluated at \(\mathbf{y}^{(k)}\) rather than at \(\mathbf{u}^{(k)}\), and the next iterate is obtained by a single projected gradient step from \(\mathbf{y}^{(k)}\): \[\label{eq:nesterov95update} \mathbf{u}^{(k+1)} = \Pi_{\mathcal{U}_{ad}}\!\left[ \mathbf{y}^{(k)} - \eta\,\hat{\mathbf{g}}^{(k)} \right],\tag{21}\] where \(\hat{\mathbf{g}}^{(k)}\) is the unconstrained gradient 19 evaluated at \((\mathbf{x}(\mathbf{y}^{(k)}), \mathbf{y}^{(k)}, \boldsymbol{\lambda}(\mathbf{y}^{(k)}))\).

Since the optimization problem associated with the controlled dynamic network is not convex in general, the extrapolated point \(\mathbf{y}^{(k)}\) may occasionally produce an accelerated step that increases the objective functional [35]. For this reason, we adopt a monotone restart strategy: if the candidate iterate obtained from 21 satisfies \[J\bigl(\mathbf{u}^{(k+1)}\bigr) > J\bigl(\mathbf{y}^{(k)}\bigr),\] then the momentum contribution is discarded for that iteration by resetting \(\mathbf{y}^{(k)} \gets \mathbf{u}^{(k)}\), recomputing the corresponding gradient at \(\mathbf{u}^{(k)}\), and replacing the accelerated trial with the standard projected-gradient step \[\mathbf{u}^{(k+1)} = \Pi_{\mathcal{U}_{ad}}\!\left[ \mathbf{u}^{(k)} - \eta\,\hat{\mathbf{g}}^{(k)} \right].\] This safeguard preserves the simplicity of the fixed-step iteration while preventing the momentum term from generating non-descent updates.

The convergence of Algorithm [alg:nesterov95pg] is monitored via the projected-gradient mapping norm, which is the standard optimality measure for projected-gradient methods [36]. At iteration \(k\), given step size \(\eta> 0\), the projected-gradient mapping is defined as \[\mathbf{g}^{(k)} = \dfrac{1}{\eta}\!\left(\mathbf{y}^{(k)} - \mathbf{u}^{(k+1)}\right),\] defined in line [line:project] of the algorithm. Note that \(\mathbf{g}^{(k)} = \mathbf{0}\) is a necessary optimality condition for the projected problem, and \(\| \mathbf{g}^{(k)} \|_2\) measures how far the current iterate is from satisfying this condition. The algorithm is declared converged when \[\frac{\|\mathbf{g}^{(k)}\|_2}{\max(1,\|\mathbf{y}^{(k)}\|_2)} \leq \varepsilon_g,\] where \(\varepsilon_g > 0\) is a user-specified tolerance. As a secondary safeguard against stagnation, a simultaneous check on both the relative change in the iterates, \[\frac{\|\mathbf{u}^{(k+1)} - \mathbf{u}^{(k)}\|_2}{\max(1,\|\mathbf{u}^{(k)}\|_2)} \leq \varepsilon_x,\] and the relative change in the objective value is also employed.

In the following subsection, we derive the expression for the state, costate, and gradient equations for the problems we formulated in Section 3.

4.1 The dynamic receive case↩︎

We start by applying the general theory to dynamic receive case from 10 . For this problem the Hamiltonian 15 can be written as \[\label{eq:our-hamiltonian} \begin{align} \mathrm{H}(\mathbf{r}, U, \boldsymbol{\lambda}, t) = \frac{\alpha}{2} \| U(t) \|^2_F + \boldsymbol{\lambda}(t)^\top \left( -b \left( \mathbf{r}(t) - \mathbf{1} \right) - \log \left( I - a (A(t) + U(t))^\top \right) \mathbf{r}(t) \right). \end{align}\tag{22}\] The next step is to derive the optimality conditions for the state 16 , costate 17 and for the gradient of the minimization problem in 18 . To obtain them we need to prove some preliminary results regarding the Fréchet derivative of the matrix logarithm function.

We can now express the entire derivative in matrix form using an expression for computing a bilinear form of the Fréchet derivative of a matrix function in a rank-1 direction [29].

With these preliminary results, we can then formulate the set of optimality conditions for the problem 10 .

In the linear case 11 , obtaining the optimality conditions is then an immediate corollary of Proposition 7.

At this stage, we have derived all the analytical expressions required to implement Algorithm [alg:nesterov95pg]. It remains to address their discretization and to develop an efficient computational strategy, which we discuss in the next Section 5.

5 Discretizing and solving the problem↩︎

To solve the problems in 10 , 11 and 12 using the PMP framework introduced in Section 4 and summarized in Algorithm [alg:nesterov95pg], we first need to discretize the expressions for the gradient, as well as the state and costate equations.

Once this is achieved, we apply a projected gradient method enhanced with the Nesterov’s acceleration to compute an approximation of the optimal control \(U(t)\) at the selected discretization time steps. In our implementation, all norms are induced by the Frobenius inner product on control matrices \(U(t)\), i.e., aggregated across all time steps. We use the tolerances \(\varepsilon_g\) and \(\varepsilon_x\) for monitoring the convergence of the projected gradient Algorithm [alg:nesterov95pg].

In the following sections, we detail the components of the resulting algorithm. Specifically, we begin by deriving an appropriate discretization scheme for the state and costate equations. We then address the efficient computation of matrix-function–vector products (Section 5.2), and finally discuss the evaluation of the projection operators (Section 5.3).

5.1 Selecting a time-stepping method for the state and costate equations↩︎

To evaluate the conditions in Proposition 7, i.e., update the state and costate in lines [line:solve_state] and [line:solve_costate] of Algorithm [alg:nesterov95pg], we need to march in time two non-autonomous ODEs \[\label{eq:state95for95r} \begin{cases} \displaystyle \dot{\mathbf{r}}(t) = \frac{\partial \mathrm{H}}{\partial \boldsymbol{\lambda}} = -b \left( \mathbf{r}(t) - \mathbf{1} \right) - \log \left( I - a (A(t) + U(t))^\top \right) \mathbf{r}(t), \\ \mathbf{r}(0) = \mathbf{1}, \end{cases}\tag{23}\] and \[\label{eq:costate95for95r} \begin{cases}\displaystyle \dot{\boldsymbol{\lambda}}(t) =-\frac{\partial \mathrm{H}}{\partial \mathbf{r}} = b \boldsymbol{\lambda}(t) + \log \left( I - a (A(t) + U(t)) \right) \boldsymbol{\lambda}(t), \\ \boldsymbol{\lambda}(T) = \mathbf{r}(T) - \mathbf{r}^*. \end{cases}\tag{24}\] In general, we consider the adjacency matrix \(A(t)\) at a finite set of time instants, \(\{A(t_k)\}_k\). However, these discrete observations are often insufficient to define a difference equation with reasonable accuracy, especially when the time intervals \(h_k = t_{k+1} - t_k\) are large. To address this, it is typically assumed that \(A(t)\) is a piecewise constant function over the intervals \([t_k, t_{k+1})\). While this assumption simplifies the representation, it results in the dynamics of the corresponding ODE being only piecewise continuous.

To effectively integrate the ODE systems 23 –24 , it is necessary to account for the discontinuities in the time-dependent adjacency matrix \(A(t)\). Under the piecewise constant assumption, the right-hand side of both equations is only piecewise continuous, with possible jumps at the time instants \(\{t_k\}_k\).

To handle this structure, we adopt an explicit Euler scheme applied separately on each interval \([t_k, t_{k+1})\) with a fixed step size \(h\). More precisely, within each interval, the dynamics are integrated assuming \(A(t)\) is constant and equal to \(A(t_k)\). The numerical solution is then advanced step by step until reaching \(t_{k+1}\).

At each switching time \(t_{k+1}\), the value of the numerical solution obtained at the end of the interval \([t_k, t_{k+1})\) is used as the initial condition for the subsequent interval \([t_{k+1}, t_{k+2})\). In this way, the integration procedure consistently propagates the solution across intervals while properly accounting for the discontinuities in \(A(t)\). While we selected the Euler method to control the computational cost, higher-order time-stepping methods can in principle improve the accuracy of the approximation.

5.2 Evaluating the matrix logarithm and its Fréchet derivative↩︎

To apply the procedure to a large time-dependent network we need a reliable way to evaluate for a large, sparse matrices \(M\) the matrix-vector product \(\log(M)\mathbf{v}\). Instead of explicitly forming \(\log(M)\), we employ Krylov-type subspace methods to provide an efficient way to approximate \(\log(M)\mathbf{v}\).

The key idea behind Krylov methods is to approximate the action of \(\log(M)\) on \(\mathbf{v}\) by working on the Krylov subspace of dimension \(k\) \[\mathcal{K}_k(M,\mathbf{v}) = \{\mathbf{v}, M\mathbf{v}, M^2\mathbf{v}, \ldots, M^{k-1}\mathbf{v}\}.\] The approximation in this subspace becomes more manageable by working on a smaller, projected, version of \(M\). The process begins by constructing an orthonormal basis \(Q_k\) for the Krylov subspace via the Lanczos method [37] for a symmetric \(M\) or via the Arnoldi method [37] for a general matrix. The result is a set of basis vectors that are used to form a small tridiagonal or Hessenberg matrix \(H_k\), respectively, which represents \(M\) in the reduced subspace as \(H_k = Q_k^\top M Q_k\). Once the projection is complete, we compute the matrix logarithm of \(H_k\) explicitly [28], [38] and obtain the approximation \[\log(M)\mathbf{v} \approx Q_k \log(H_k) (Q_k^\top \mathbf{v}).\] For the large number of evaluations that we need to march in time the solution to the state, costate equation and Hamiltonian gradient in Proposition 7 polynomial Krylov methods are particularly appealing. To march in time the solution of the differential equation for the state and the costate we can resort to the Lanczos method in the case where the underlying network is undirected at every instant of time and we seek a symmetric control \(U(t)\) to preserve the overall structure.

To evaluate the gradient of the Hamiltonian we need to compute the Fréchet derivative of the matrix logarithm. By 1 , the Fréchet derivative \(\mathrm{L}_{\log}(M,E)\) can be extracted from the matrix logarithm of a \(2n \times 2n\) block matrix. For each direction vector \(\mathbf{v}\), this amounts to computing the \(n\)-dimensional vector \[\mathrm{L}_{\log}(M,E)\mathbf{v} = \begin{bmatrix} I & O \end{bmatrix} \log \begin{pmatrix} M & E \\ O & M \end{pmatrix} \begin{bmatrix} \mathbf{0} \\ \mathbf{v} \end{bmatrix}, \quad M = I - a(A(t_k) + U_k),\] for \(U_k\) the approximation of \(U(t_k)\) induced by the choice of the time-stepping procedure. Although polynomial Krylov methods perform well for the state and costate equations, where the matrix \(M\) has spectrum in the right half-plane and the Arnoldi projection preserves this property reliably, the situation is different for the \(2n \times 2n\) block matrix above. For this doubled system, the projected Hessenberg matrix \(H_k\) can develop eigenvalues with negative real part, even when the spectrum of the original block matrix lies entirely in the right half-plane. Since the principal logarithm has a branch cut on \((-\infty,0]\), such spurious eigenvalues render the evaluation of \(\log(H_k)\) unreliable. To address this, we employ rational Krylov methods [27], [39], which replace the polynomial Krylov subspace by the rational Krylov subspace \[\mathcal{Q}_k(M,\mathbf{b},\Xi) = \operatorname{span}\left\{ \prod_{i=1}^{k}(M - \xi_i I)^{-1} \, p(M)\mathbf{b} \;:\; p \in \mathcal{P}_{k-1} \right\},\] where \(\Xi = \{\xi_1,\ldots,\xi_k\}\) is a prescribed set of poles. Poles at infinity recover the standard polynomial Krylov step. The rational Arnoldi decomposition [26] produces an orthonormal basis \(V_k\) and upper-Hessenberg pencil \((K_k,H_k)\) such that \[M \, V_k K_k = V_k H_k.\] The small projected matrix \(\hat{M}_k = K_k^{-1} H_k\) captures the spectral information of \(M\) relevant to the starting vector, and the approximation reads \[\label{eq:rat95krylov} \log(M)\mathbf{b} \approx V_k\, K_k \log(\hat{M}_k)\, (K_k^{-1} e_1 \|\mathbf{b}\|_2).\tag{25}\] In practice, the quality of the rational approximation depends crucially on the choice of poles. We select them automatically by computing an AAA rational approximant [40] of \(\log(x)\) on a set of sample points covering the interval \([1 - a \max_t \rho(A(t)), 1 + a \max_t \rho(A(t))]\), where \(a > 0\) is the edge attenuation parameter of the network model. The poles of the resulting rational function provide near-optimal shifts for the rational Krylov iteration. In this way, the number of poles (and therefore the Krylov subspace dimension) is determined adaptively by the AAA algorithm and is typically small (on the order of \(10\)–\(20\)) even when polynomial Krylov methods would require significantly more iterations.

We apply this rational Krylov strategy to the \(2n \times 2n\) block matrix from 1 . The block upper-triangular structure of this matrix is exploited in the operator-vector products and shifted system solves required by the rational Arnoldi iteration: each shifted solve reduces to two triangular back-substitutions with the \(n\times n\) shifted matrix \(\nu M - \mu I\), preserving the sparsity of the original problem.

5.3 Computing the projections↩︎

The projected gradient method in Algorithm [alg:nesterov95pg] requires the computation of the projection operator \(\Pi_{\mathcal{U}_{ad}}:\mathbb{R}^{m}\rightarrow\mathcal{U}_{ad}\) onto the admissible control set. The dimension \(m\) corresponds to the number of independent control variables after taking into account the sparsity pattern \(\mathcal{P}\). In the general directed case, \(m\) is equal to the number of admissible entries of the pattern. In the symmetric case 14 , the number of variables is further reduced by considering only one triangular part of the control matrix, since the remaining entries are imposed by the condition \(U(t)=U^\top(t)\).

The sparsity constraint \([U(t)]_{i,j}=0\) if \((i,j)\notin\mathcal{P}\) is therefore incorporated directly into the choice of optimization variables and does not need to be enforced explicitly during the projection step.

For the admissible sets 13 and 14 , the constraint \(U(t)\geq -A(t)\) is an elementwise lower bound. After discretizing the control variables in time and introducing the vectorization operator \(\operatorname{vec}(\cdot)\), the projection onto this constraint reduces to a componentwise clipping operation.

Let \(\mathbf{U}^{(k)}\) denote the vectorized control at iteration \(k\) and let \(G^{(k)}\) denote the gradient of the objective functional. The projected gradient update is \[\begin{align} \mathbf{U}^{(k+1)} &= \Pi_{\mathcal{U}_{ad}} \left( \mathbf{U}^{(k)} -\eta G^{(k)} \right) \nonumber \\ &= \max \left\{ \mathbf{U}^{(k)} -\eta \left( \alpha\mathbf{U}^{(k)} + a\,\operatorname{vec} (\boldsymbol{\lambda}\mathbf{r}^{\top}) \right), -\operatorname{vec}(\mathbf{A}) \right\}, \end{align}\] where the maximum is taken componentwise. Here \(\mathbf{A}\) denotes the vectorized adjacency matrix corresponding to the considered time step.

For the admissible set ?? , the constraint involving the Katz vector is given by \[U(t)\boldsymbol{\mu}(t) \leq \left(\frac{1}{a}-\epsilon\right)\mathbf{1}.\] Using the identity \(\operatorname{vec}(U\boldsymbol{\mu}) = (\boldsymbol{\mu}^{\top}\otimes I_n)\operatorname{vec}(U)\), the constraint can be written as a linear inequality in the vectorized control variables: \[\begin{align} (\boldsymbol{\mu}(t)^\top\otimes I_n)\mathbf{U}(t) \leq \left(\frac{1}{a}-\epsilon\right)\mathbf{1}. \end{align}\] The unconstrained gradient step is first computed independently at each time instant. E.g., for the logarithmic dynamics, this gives \[\begin{align} \widetilde{\mathbf{U}}^{(k)} = \mathbf{U}^{(k)} - \eta \left[ \alpha\mathbf{U}^{(k)} + a\,\operatorname{vec} \left( \mathrm{L}_{\log} \left( I-a(A(t_k)+U^{(k)})^\top, \boldsymbol{\lambda}(t_k)\mathbf{r}(t_k)^\top \right) \right) \right]. \end{align}\] The projection is then obtained by solving the quadratic program \[\begin{align} \min_{\mathbf{U}\in\mathbb{R}^{m}} \quad& \frac{1}{2} \left\| \mathbf{U}-\widetilde{\mathbf{U}}^{(k)} \right\|_2^2 \\ \text{s.t.}\quad& (\boldsymbol{\mu}(t_k)^\top\otimes I_n)\mathbf{U} \leq \left(\frac{1}{a}-\epsilon\right)\mathbf{1}, \\ & \mathbf{U}\geq-\mathbf{A}(t_k). \end{align}\] This projection is a convex quadratic program because both the objective and the constraints are convex. The case in 13 can be handled in an entirely analogous manner, leading to a similar quadratic program where the first linear constraint is replaced by \((\mathbf{1}^{\top} \otimes I_n)\mathbf{U} \leq \left(\frac{\pi}{a}-\epsilon\right)\mathbf{1} - A(t_k)\mathbf{1}\). Finally, for the symmetric admissible set 14 , the equality constraint \(U(t)=U^\top(t)\) is enforced by a reduced parametrization. Let \(\mathcal{I}_{\mathcal{P}}\) denote the set of admissible entries in the upper triangular part of the sparsity pattern. The independent control variables are collected in \(\mathbf{u}(t) = \{U_{ij}(t):(i,j)\in\mathcal{I}_{\mathcal{P}},\;i\leq j\}\). A linear expansion operator \(S\) reconstructs the complete symmetric matrix according to \(\operatorname{vec}(U(t)) = S\mathbf{u}(t)\). The projected gradient step is therefore performed in the reduced variable space: \(\mathbf{u}^{(k+1)} = \Pi_{\mathcal{U}_{ad}^{s}} \left( \mathbf{u}^{(k)} -\eta g^{(k)} \right)\), where the reduced gradient is obtained from the full matrix gradient through \(g^{(k)} = S^\top \operatorname{vec} \left( \nabla_U J(U^{(k)}) \right)\). For the elementwise admissible constraint, the projection remains a componentwise clipping operation, but it is applied only to the independent triangular variables \[\begin{align} u_{ij}^{(k+1)} = \max \left\{ u_{ij}^{(k)} -\eta g_{ij}^{(k)}, -A_{ij}(t_k) \right\}, \qquad (i,j)\in\mathcal{I}_{\mathcal{P}}. \end{align}\] The remaining entries are recovered automatically by symmetry, namely \(U_{ji}=U_{ij}\). Consequently, the equality constraint \(U=U^\top\) is satisfied by construction and no additional projection step is required.

6 Numerical examples↩︎

In this section we consider numerical experiments to validate the theory discussed in Section 4. In Section 6.1 we discuss a few small-size cases for which we can easily visualize the obtained controls and states to verify that the results are the expected ones. Then, in Section 6.2, we investigate a larger case scenario in which we employ the techniques discussed in Section 5 to reduce the computational costs. The implementation is done in Matlab (v9.10.0.1602886 (R2021a)), and the experiments are run on a single node (AMD EPYC 7763, 2 Threads per core, 64 cores per socket, 2 sockets, 2   of RAM, Ubuntu 22.04.4 LTS) of the Toeplitz cluster at the Green Data Center (GDC) of the University of Pisa. For the AAA algorithm [40] we employ the implementation available in the Chebfun [41] package (v5.7.0). For the generation of the Rational Krylov subspace we use instead the Rational Krylov Toolbox for Matlab (RKToolbox) [26] (v2.9).

6.1 Small descriptive examples↩︎

The focus of this section is on small-scale examples where the results can be effectively visualized. We begin with the case of the temporal network illustrated in Figure 3.

This network, originally presented in [7], is a synthetic telephone-call example with \(n = 17\) nodes evolving through seven snapshots. The nodes represent callers and are visibly organized in the form of a tree, at time step \(0\) two callers start the diffusion of an information that triggers a cascade of calls between the callers. We first compute the uncontrolled dynamic receive centrality \(\mathbf{r}(t)\) over \(t \in [0,7]\), which measures the importance as receiver at time \(t\) of the different callers, then we define a perturbed target by setting \(\mathbf{r}^* \equiv \mathbf{r}(7)\), \(\mathbf{r}_3^* = \mathbf{r}_4^* = 1.2\), while keeping all other components equal to the uncontrolled final state (Figure [fig:switching-tree-shortcut:states]). Namely, we want to just increase the importance as receivers of two of the callers.

The admissible pattern is \(\mathcal{P} = \{ (i,j) \in V\times V\, :\, (i,j) \in E_t \text{ for some } t\in[0,T]\}\), that is, every edge that appears at least once in the temporal sequence can be controlled. Time integration is performed with the forward Euler scheme discussed in Section 5.1, using step size \(h=10^{-2}\). We then solve problem 10 by Algorithm [alg:nesterov95pg] for three values of the Tikhonov weight \(\alpha\in\{1,0.5,0.05\}\), and \(\eta = 10^{-1}\).

The convergence curves in Figure 4 (a) show the expected trade-off between control effort and iteration count. For \(\alpha=1\) the method reaches \(\|\mathbf{g}^{(k)}\|\approx 8.43\times 10^{-7}\) in \(24\) iterations, while for \(\alpha=0.5\) it reaches \(\|\mathbf{g}^{(k)}\|\approx 8.80\times 10^{-7}\) in \(31\) iterations; the low-regularization case \(\alpha=0.05\) requires substantially more iterations (\(604\)) and reaches \(\|\mathbf{g}^{(k)}\|\approx 1.40\times 10^{-6}\).

The objective values in Figure 4 (b) are consistent with this behavior: reducing \(\alpha\) yields a lower final cost, from \(J\approx 4.80\times 10^{-3}\) (\(\alpha=1\)) to \(J\approx 3.87\times 10^{-3}\) (\(\alpha=0.5\)) and \(J\approx 9.31\times 10^{-4}\) (\(\alpha=0.05\)). Panel [fig:dynamic-receive-panel-c] confirms that the controlled terminal state moves the two modified components toward the prescribed target, and that lower regularization improves the terminal matching while preserving the qualitative structure of the uncontrolled profile on the remaining components.

We repeat the same test for the linear dynamics in 11 , using exactly the same target-design procedure as above. Namely, we first compute the uncontrolled trajectory up to \(t=7\), then define the desired terminal state by replacing only the components at nodes \(3\) and \(4\) with the value \(1.2\), while leaving all other entries unchanged. Figure [fig:switching-tree-shortcut-linear:results] reports the corresponding convergence histories and final objective values.

The linear case exhibits the same qualitative dependence on the Tikhonov weight \(\alpha\) observed in the logarithmic model. From Figure 5, for \(\alpha=1\) the method converges in \(54\) iterations with \(\|\mathbf{g}^{(k)}\|\approx 8.33\times 10^{-7}\), for \(\alpha=0.5\) in \(102\) iterations with \(\|\mathbf{g}^{(k)}\|\approx 9.93\times 10^{-7}\), and for \(\alpha=0.05\) in \(615\) iterations with \(\|\mathbf{g}^{(k)}\|\approx 1.75\times 10^{-6}\). As expected, weaker regularization reduces the final value of the cost functional: Figure [fig:switching-tree-shortcut-linear:results-b] gives \(J\approx 7.35\times 10^{-3}\) (\(\alpha=1\)), \(J\approx 6.12\times 10^{-3}\) (\(\alpha=0.5\)), and \(J\approx 1.92\times 10^{-3}\) (\(\alpha=0.05\)). Overall, this confirms the same accuracy-versus-effort trade-off: smaller \(\alpha\) yields better terminal matching at the price of more iterations.

6.2 Larger scale example↩︎

As a larger-scale example, we consider the CollegeMsg temporal network [42], a dataset of private messages exchanged on an online social network at the University of California, Irvine. The dataset records directed interactions: an edge \((u, v, t)\) indicates that user \(u\) sent a private message to user \(v\) at Unix timestamp¹ \(t\), yielding a total of \(n = 1899\) users and over \(59\,000\) timestamped communications.

To construct the time-varying adjacency matrices, we partition the full observation period uniformly into \(n_t = 30\) time windows and normalize the resulting breakpoints to the interval \([0, 10]\); see Figure [fig:messagenet-window]. For each window \(i = 1, \ldots, n_t - 1\), we collect all edges with timestamps in \([t_i, t_{i+1})\) and assemble the sparse adjacency matrix \(A_i \in \mathbb{R}^{n \times n}\), where each nonzero entry records an observed message. The communication activity varies considerably across windows, as reflected by the spread in the number of nonzeros \(\operatorname{nnz}(A_i)\) shown in Figure [fig:messagenet-nnz]. The parameter \(a\) in the centrality model is set to \(a = 1/(\rho_{\max} + 1)\), where \(\rho_{\max} = \max_i \rho(A_i)\) is the largest spectral radius across all snapshots, ensuring that \(I - aA_i\) admits a matrix logarithm throughout the horizon; see again Remark 1.

Before running the full optimization, we validate the Krylov method for the matrix logarithm action described in Section 5. For each snapshot \(A_i^\top\), we apply logmv with tolerance \(10^{-6}\) and a maximum of \(40\) Arnoldi iterations to compute \(\log(I - a A_i^\top)\mathbf{e}\), where \(\mathbf{e}\) is a fixed random Gaussian vector, and compare the result against the dense reference [28] \(\texttt{logm}(I - aA_i^\top)\mathbf{e}\). The absolute and relative errors, together with the number of Krylov iterations \(k\) required to meet the stopping criterion, are reported in Table 1. Across all \(29\) snapshots the method converges in at most \(18\) iterations, and both absolute and relative errors remain well below \(10^{-7}\), confirming that the Krylov-based procedure from Section 5 is an effective method for the matrices arising in this dataset.

We perform a similar validation for the Fréchet-derivative action used in the gradient evaluation, \[\mathrm{L}_{\log}(A_i^\top,\boldsymbol{\lambda}\mathbf{r}^\top)\mathbf{v}.\] For each snapshot \(A_i^\top\), we set \(M_i = I-aA_i\), generate fixed random Gaussian vectors \(\boldsymbol{\lambda},\mathbf{r},\mathbf{v}\in\mathbb{R}^n\), and compare the Krylov approximation produced by the rational Krylov approximation 25 against the dense reference \(\mathbf{y}_{\mathrm{ref}} = \mathrm{L}_{\log}(M_i^\top,\boldsymbol{\lambda}\mathbf{r}^\top)\mathbf{v}\), computed with logm_frechet_pade [28]. For the Fréchet action we now use the rational Krylov routine with a fixed set of \(11\) poles (obtained by AAA on the \([1 - a \max_t \rho(A(t)), 1 + a \max_t \rho(A(t))]\) interval). The absolute and relative errors are reported in Table 2. Across all \(29\) snapshots, the relative errors range from \(1.07\times10^{-14}\) to \(1.60\times10^{-7}\), confirming that the rational Krylov Fréchet action is highly accurate for gradient evaluations in this dataset.

We now consider the general problem using the CollegeMsg temporal network with approximately 1,899 nodes and 29 snapshots. The adjacency matrices are interpolated at a uniform time step \(\Delta t = 0.1\) across the normalized time domain, yielding \(116\) discrete time levels for optimization. Following Remark 2, the edge attenuation parameter is set to \(a = 1/(2\rho_{\max})\), where \(\rho_{\max} = \max_k \rho(A_k)\) is the maximum spectral radius. The temporal decay rate is \(b = 0.85\), and the control regularization weight is \(\alpha = 1\). For the Krylov approximations, we use tolerance \(10^{-6}\) and maximum 40 iterations for polynomial actions on state/costate, and a fixed set of 11 poles (chosen via AAA) for rational Krylov Fréchet actions. To decide the set of edges of the temporal network we can modify, we first build the sparsity pattern \(P\) as the union of all temporal adjacency patterns, then we restrict it to only the edges incident to 50 randomly selected nodes, creating \(814\) controllable entries per time step. Feasibility is enforced via the Katz-based linear constraint from Proposition 1, ensuring \(\rho(A(t) + U(t)) < 1/a\) at all times through convex quadratic projections with numerical margin \(10^{-6}\). To decide the target state, we first simulate the uncontrolled receive centrality \(\mathbf{r}^{\mathrm{nc}}(T)\) to terminal time \(T = 10\). The target state \(\mathbf{r}^*\) is a perturbed version: for each of the 50 controllable nodes, we set \(r^*_i = \max(2 \cdot r^{\mathrm{nc}}_i \cdot \xi, 1)\) where \(\xi \sim \mathrm{U}(0,1)\) is uniform random, creating a non-trivial steering objective.

We employ the Nesterov accelerated projected-gradient method with monotone restart from Algorithm [alg:nesterov95pg]. The algorithm runs for up to \(2000\) iterations with fixed step size \(\eta = 0.10\), momentum coefficient and the three stopping criteria: \(\|\nabla_U J\| \le 10^{-6}\), relative step size \(\le 10^{-8}\), or relative objective change \(\le 10^{-8}\). Monotone restart rejects momentum if an extrapolated gradient step increases the objective, enhancing robustness.

At each iteration, the algorithm computes the state trajectory (forward Euler), costate trajectory (backward adjoint, reverse time), and gradient via matrix-free Krylov evaluations. The control gradient is computed employing the Fréchet action via rational Krylov with AAA poles on \([1-a \max_{t}\rho(A(t)),1+a \max_{t}\rho(A(t))] \approx [a, 2]\). Feasibility projections via quadratic programming are applied to each time-step control after the gradient step. The accuracy of the Krylov approximations (as validated in Table 2) ensures that repeated state/costate/gradient evaluations remain stable throughout the optimization.

Figure 7 shows the convergence history of the Nesterov accelerated projected-gradient solver on the full temporal network problem.

We note that even in the case of a larger-scale problem the behavior with respect to the regularization parameter is analogous to what was seen in the cases in Section 6.1, the number of necessary iterations increases as \(\alpha\) decreases, and the method declares convergence at iteration 862 with a projected gradient norm \(\approx 3.663 \times 10^{-6}\), and a total convergence cost of \(J \approx 3.758922\). The gradient tolerance is not satisfied, but the difference between two successive iterates has dropped below the threshold of \(10^{-7}\) reaching \(9.9 \times 10^{-8}\) and causing convergence to be declared.

7 Conclusion and future perspectives↩︎

In this work we formulated and analyzed a pattern-constrained optimal control problem for dynamic receive centrality on temporal networks. Using the Pontryagin Maximum Principle, we derived state, costate, and gradient conditions for both the logarithmic model and its linearized counterpart, and we introduced a practical projection strategy based on a sufficient linear constraint that preserves well-posedness of the logarithmic dynamics. On the algorithmic side, we adopted a Nesterov accelerated projected-gradient scheme with monotone restart, and we combined it with matrix-free Krylov approximations for both matrix-logarithm actions and Fréchet-derivative actions. The numerical results show that this combination is effective on both synthetic and real temporal data: the optimization procedure remains stable under pattern constraints, and Krylov-based evaluations provide accurate approximations at moderate subspace dimensions, making repeated state/costate/gradient evaluations computationally viable.

Several extensions can follow from the present framework. A first direction is a sharper theoretical analysis of convergence for the projected accelerated method in the nonconvex setting induced by the logarithmic dynamics, including quantitative guarantees for restart strategies. A second direction is to further reduce the cost of the Fréchet term by exploiting block, low-rank, or recycling Krylov spaces across consecutive time steps and optimization iterations. A third direction is to investigate adaptive and higher-order time integration tailored to piecewise-smooth temporal networks, balancing accuracy and runtime in large horizons. From a modeling perspective, it is also of interest to consider multi-objective formulations that jointly steer receive centrality and other dynamic-network functionals, as well as robust or stochastic variants that account for uncertainty in temporal edges. Finally, extending the framework to distributed and parallel implementations would be essential for very large networks, where sparse linear algebra and matrix-function actions become the dominant computational bottleneck.

All the authors are members of the INdAM GNCS and acknowledge the MUR Excellence Department Project awarded to the Department of Mathematics, University of Pisa, CUP I57G22000700001.

This research has been partially funded by the Italian Ministry of University and Research (MUR) through the PRIN 2022 “MOLE: Manifold constrained Optimization and LEarning”, code: 2022ZK5ME7 MUR D.D. financing decree n. 20428 of November 6th, 2024 (CUP B53C24006410006) and Fondo Italiano per la Scienza (FIS2023-01334) advanced grant “ADvanced numerical Approaches for MUltiscale Systems with uncertainties"- ADAMUS.

The code to reproduce the numerical examples is available on the GitHub repository Cirdans-Home/optimal-network-control.

Declarations↩︎

Not applicable The authors declare no competing interests.

References↩︎

---

1. The Unix timestamp is the number of seconds between a particular date and the Unix Epoch (January 1^st, 1970 UTC).↩︎