Optimal Graph Laplacian

This paper provides a construction method of the nearest graph Laplacian to a matrix identified from measurement data of graph Laplacian dynamics that include biochemical systems, synchronization systems, and multi-agent systems. We consider the case where the network structure, i.e., the connection relationship of edges of a given graph, is known. A problem of finding the nearest graph Laplacian is formulated as a convex optimization problem. Thus, our problem can be solved using interior point methods. However, the complexity of each iteration by interior point methods is $O(n^6)$, where $n$ is the number of nodes of the network. That is, if $n$ is large, interior point methods cannot solve our problem within a practical time. To resolve this issue, we propose a simple and efficient algorithm with the calculation complexity $O(n^2)$. Simulation experiments demonstrate that our method is useful to perform data-driven modeling of graph Laplacian dynamics.


Introduction
Many networked systems can be described as graph Laplacian dynamicṡ where x(t) ∈ R n and L ∈ R n×n denote the state and graph Laplacian, respectively; e.g., system (1) includes biochemical systems [1,2,3,4], synchronization systems [5,6,7], and multi-agent systems [8,9,10]. The graph Laplacian L is determined by the connection relationship of edges of a given graph (i.e., the network structure) and weights of the edges. In contrast to the identification of the network structure, it is difficult to identify the weights due to the lack of sensor measurements and sensor noises. That is, it is difficult to identify the graph Laplacian L using the existing identification methods such as subspace identification methods [11,12] and dynamic mode decomposition methods [13,14,15]. In other words, the existing identification methods may provide a matrix A ∈ R n×n that dose not have the graph Laplacian property. For this reason, this paper proposes a construction method of the nearest graph Laplacian L to a given matrix A under the assumption that the network structure is known. This assumption is motivated from the fact that it is possible to identify the network structure for many examples [16,17]. We can find related works in [18,19,20,21,22]. The authors of [19,21] studied the problem of enforcing the stability of a system a posteriori. That is, [19,21] proposed methods for finding the nearest stable matrix to a given unstable matrix. Reference [18] considered the problem of finding the nearest stable Metzler matrix to a given non-Metzler matrix. The authors of [20,22] proposed methods for finding the nearest symmetric positive semidefinite matrix with unit diagonal to a given symmetric matrix. However, to the best of our knowledge, there has been no previous work regarding construction of optimal graph Laplacian.
The contributions of this paper are summarized as follows. We show that the aforementioned optimal graph Laplacian construction problem can be formulated as a convex optimization problem with an entrywise 1-norm objective function in contrast to the previous works in [18,19,21]. In fact, the problems in [18,19,21] are non-convex optimization problems with an entrywise 2-norm objective function due to the stability constraint, while the constraint is not required for our problem. Thus, in contrast to [18,19,21], we can obtain a global optimal solution using interior point solvers such as CVX [23]. However, the complexity of each iteration by interior point methods is O(n 6 ) [24]. That is, if n is large, interior point methods cannot solve our convex optimization problem within a practical time. In order to solve our problem with n ≥ 1000 within a practical time, we develop a simple and effective algorithm with the calculation complexity O(n 2 ). Furthermore, we demonstrate that if we replace our objective function with an entrywise 2-norm function, we cannot develop such simple algorithm. Simulation experiments illustrate that our proposed method is useful to perform data-driven modeling of graph Laplacian dynamics (1).
Notation: The set of real numbers is denoted by R. The symbols 0 n ∈ R n and 1 n ∈ R n are column vectors with all zero entries and all one entries, respectively. For any real number a, |a| denotes the absolute value of a. Given a matrix A ∈ R n×n , ||A|| 1 and ||A|| 2 denote the entrywise 1-norm and 2-norm, respectively. That is,

Problem formulation
This section formulates our problem. To this end, let G = (V, E, W) be a weighted graph, where V = {1, 2, . . . , n} is the node set, E ⊂ V × V is the edge set, and W is the adjacency matrix consisting of non-negative elements w i j called the weights. That is, for each edge (i, j) ∈ E, the i-th row and j-th column entry of W equals the weight w i j > 0, and all other entries of W are equal to zero. The degree matrix of G is a diagonal matrix denoted by D = diag(d 1 , d 2 , . . . , d n ), with d i = n j=1 w i j . The graph Laplacian L of G is defined as L := D − W. Thus, the graph Laplacian L has the following properties: • The sum of row elements is equal to zero; i.e., • The diagonal elements of L are non-negative and non-diagonal elements are non-positive; i.e., Here, the set S 1 is defined as Conversely, we call a matrix L ∈ R n×n (n ≥ 2) graph Laplacian if L satisfies (2) and (3) [8]. It follows from (2) that the graph Laplacian L has at least one zero eigenvalue. Furthermore, the eigenvalues of the graph Laplacian L different from zero have strictly-positive real parts [8]; i.e., −L in (1) is a stable matrix. We assume that the network structure, i.e., the edge set E, is known. This assumption comes from the fact that it is possible to identify the network structure for many examples [16,17]. The optimal solution to the following problem provides the nearest graph Laplacian to a given matrix A in the case where the graph G is a directed graph. Here, the set S 2 is defined as which indicates the connection relationship of the edges of the graph G.
Problem 1 is a convex optimization problem. This is because the set of all L satisfying constraint (2) is a vector space, i.e., a convex set, and the sets S 1 and S 2 are also convex, and ||A − L|| 1 is a convex function [25]. Thus, Problem 1 can be solved by interior point solvers such as CVX [23]. However, the complexity of each iteration by interior point methods is O(n 6 ) [24]. That is, if n is large, interior point methods cannot solve Problem 1 within a practical time, as shown in Section 5.
Although the objective functions of [18,19,21] are entrywise 2-norms, the function of Problem 1 is an entrywise 1-norm. This is because if we replace ||A − L|| 1 with ||A − L|| 2 2 , an algorithm for solving the modified problem becomes more complicated than the case of Problem 1, as explained in Section 4.

Main results
This section proves the following theorem.
Theorem 1. Algorithm 1 provides a global optimal solution to Problem 1.
Here, Π S 1 ∩S 2 : R n×n → S 1 ∩ S 2 in Algorithm 1 denotes the projection onto the closed convex set S 1 ∩ S 2 . That is, where i j. Note that Theorem 1 holds for any weighted graph G and any matrix A. Algorithm 1 produces L with where α i := n j=1L i j . Because the matrix L satisfies the constraint conditions of Problem 1, we prove that the matrix L minimizes ||A − L|| 1 . If α 1 = α 2 = · · · = α n = 0, we have L =L. That is, in this case,L is the global optimal solution to Problem 1. Suppose that there exists i such that α i 0, i.e., L L , and L is not a global optimal solution to Problem 1. That is, a global optimal solution L * to Problem 1 satisfies Then, there exists i such that n j=1 |d i j | < |α i | and |α i | > 0, where d i j := L * i j −L i j . This is because if n j=1 |d i j | ≥ |α i | for any i, (5) does not hold. From the definitions of α i and d i j , we obtain that If α i > 0, (6) and (7) imply that If α i < 0, (6) and (7) also imply that Thus, if (5) holds, L * 1 n 0 n . This is a contradiction that L * is a solution to Problem 1. Hence, L satisfying (4) and the constraint conditions of Problem 1 is a global optimal solution to Problem 1. This completes the proof. The calculation complexity of Algorithm 1 is O(n 2 ). Thus, Algorithm 1 is considerably more efficient than interior point methods. In Section 5, we demonstrate this fact.

Comments on another objective function
This section explains that if we replace the objective function ||A − L|| 1 of Problem 1 with ||A − L|| 2 2 , then we cannot obtain a simple algorithm such as Algorithm 1. To this end, we consider the case where the graph G is a complete graph and A ∈ S 1 ∩ S 2 . That is, we demonstrate it using the simplest case.
We consider the following optimization problem.
This is because under the assumption that G is a complete graph, the matrix L constructed by using the solutions (L i1 , L i2 , . . . , L in ) (i = 1, 2, . . . , n) to Problem 2 is a global optimal solution to minimize L∈R n×n subject to (2), (3), and L ∈ S 2 .
To explain that an algorithm for solving Problem 2 becomes more complicated than Algorithm 1, we need the following lemma.
Lemma 1. The unique global optimal solution to minimize where a 1 , a 2 , . . . , a n ∈ R.
Proof. The above optimization problem can be reduced to minimize Thus, we obtain that , ∂g ∂x 2 = ∂g ∂x 3 = · · · = ∂g ∂x n = 0 and n i=1 x i = 0 imply (9). Because original problem (8) is a convex optimization problem, (9) is a global optimal solution to (8). Furthermore, a point satisfying the Karush-Kuhn-Tucker (KKT) condition [25] is unique. Thus, (9) is the unique global optimal solution to (8).  From Lemma 1, if n j=1 A i j > 0, then the unique global optimal solution to Problem 2 is given by This is because (10) implies that L ii ≥ 0 and L i j ≤ 0 (i j). However, if n j=1 A i j < 0, then (10) does not guarantee L i j ≤ 0 (i j). That is, L defined by (10) is not a global optimal solution to the modified problem of Problem 1 in general. For this reason, we cannot develop a simple algorithm such as Algorithm 1 for solving the modified problem, even if the graph G is a complete graph and A ∈ S 1 ∩ S 2 .

Numerical experiments
This section numerically compares Algorithm 1 and CVX [23] which is a popular solver for solving convex optimization problems. Furthermore, we discuss eigenvalues of the graph Laplacian L produced by Algorithm 1. All computations were carried out using MATLAB R2017b on an Intel(R) Xeon(R) CPU E5-2637 v4 @ 3.50 GHz 3.50 GHz and 128 GB RAM.
We generated the matrix A in Problem 1 by the following steps.
1. Generate the graph using the Watts and Strogats model [26]. Here, the number of the nodes is n.
2. Replace all nonzero elements of the adjacency matrix of the graph with 10 × rand, where rand is a single uniformly distributed random number in the interval (0, 1). Set the modified adjacency matrix as X. 4. L * := Y − X. 5. A := L * + s × randn(n), where s > 0, and randn(n) denotes an n × n matrix of normally distributed random numbers.
We here note that the Watts and Strogats model has three parameters (n, K, β), where K and β denote the mean degree and rewiring probability, respectively. Table 1 shows the computational times (seconds) for different n when s = 5, K = 10, and β = 0.3. When n = 1000, 5000, and 10000, CVX could not solve Problem 1 due to the out of memory. According to Table 1, Algorithm 1 is considerably faster than CVX. In particular, Algorithm 1 could solve Problem 1 with n = 10000 within a practical time.

Comparison of Algorithm 1 and CVX
Real part I maginary part Figure 1: Eigenvalues of L * , A, and L.  Fig. 1 illustrates eigenvalues of the matrices L * , A, and L when n = 300, s = 5, K = 10, and β = 0.3. The eigenvalues of L * and L were more similar than those of L * and A. That is, we could construct the graph Laplacian L near the graph Laplacian L * from the matrix A in the sense of the eigenvalues. This is a preferable result if the matrix A can be regarded as a perturbed matrix of the graph Laplacian L * . If this is the case, it is important that the second smallest real parts of eigenvalues of L * and L are near. This is because those determine the consensus speed of multi-agent system (1) [27]. According to Fig. 1, those of L * and L are close, although those of L * and A are too different.
Moreover, Table 2 shows relations among the parameter s, Ave, and Var when n = 300, K = 10, and β = 0.3, where Here, Re(λ 2,i ) and Re(λ * 2,i ) denote the second smallest real parts of eigenvalues of L and L * , respectively, at the i-th trial. Ave and Var both monotonically increase as s increases. According to Table 2, if 0 < s ≤ 5, we can expect that the second smallest real parts of eigenvalues of L * and L are sufficiently close. That is, even if each component of A is relatively different from the corresponding component of L * , Algorithm 1 can generate L close to L * in the sense of the second smallest real parts of eigenvalues of L * and L.

Comments on an application to data-driven modeling
The above results conclude that Algorithm 1 is useful to perform data-driven modeling of graph Laplacian dynamics (1). In fact, to perform the data-driven modeling, it is desirable that • we can construct the graph Laplacian L from the matrix A within a very short time.
• the true graph Laplacian L * and the constructed graph Laplacian L are sufficiently close. Table 1 indicates that Algorithm 1 can produce the graph Laplacian within a very short time in contrast to CVX even if n ≈ 1000. Furthermore, according to Section 5.2, L * and L are sufficiently close in the sense of the second smallest real parts of the eigenvalues, even if the matrix A is relatively far from L * .

Conclusion
We have provided a simple and efficient algorithm with the calculation complexity O(n 2 ) for solving a convex optimization problem of constructing the nearest graph Laplacian to a given matrix. Simulation results have demonstrated that our proposed method is useful to perform data-driven modeling of graph Laplacian dynamics.