Identification of critical nodes in large-scale spatial networks
IIdentification of critical nodes in large-scale spatial networks ∗ Vishaal Krishnan Sonia Mart´ınez † Abstract
The notion of network connectivity is used to characterize the robustness and failuretolerance of networks, with high connectivity being a desirable feature. In this paper, wedevelop a novel dynamical approach to the problem of identifying critical nodes in large-scale networks, with algebraic connectivity (the second smallest eigenvalue of the graphLaplacian) as the chosen metric. Employing a graph-embedding technique, we reduce theclass of considered weight-balanced graphs to spatial networks with uniformly distributednodes and nearest-neighbors communication topologies. Through a continuum approxi-mation, we consider the Laplace operator on a manifold (with the Neumann boundarycondition) as the limiting case of the graph Laplacian. We then reduce the critical nodeset identification problem to that of finding a ball of fixed radius, whose removal mini-mizes the second (Neumann) eigenvalue of the Laplace operator on the residual domain.This leads us to consider two functional and nested optimization problems. Resortingto the Min-max theorem, we first treat the problem of determining the second smallesteigenvalue for a fixed domain by minimizing an energy functional. We then obtain aclosed-form expression for a projected gradient flow that converges to the set of pointssatisfying the KKT conditions and provide a novel proof that the only locally asymptoti-cally stable critical point is the second eigenfunction of the Laplace operator. Building onthese results, we consider the critical ball identification problem and define novel dynam-ics to converge asymptotically to these points. Finally, we provide a characterization ofthe location of critical nodes (for infinitesimally-small balls) as those points which belongto the nodal set of the second eigenfunction of the Laplacian operator.
The identification of critical nodes in a network is motivated by the question of network ro-bustness and is crucial to improving its resilience to attacks and failures. The notion of criticalnodes refers to the subset of nodes in the network whose removal results in the maximumdeterioration of a given performance metric. In the context of robustness of networks/graphs,a widely studied metric [1, 2] is the second smallest eigenvalue of the graph Laplacian matrix(also called the algebraic connectivity of the graph). In addition to being an indicator ofhow well connected the graph is, it is typically of significance in the context of agreementdynamics on networks (such as consensus and synchronization), as it governs the convergencerate of the dynamics. ∗ This work has been partially supported by grant FA9550-18-1-0158. † The authors are with the Department of Mechanical and Aerospace Engineering, University of Californiaat San Diego, La Jolla CA 92093 USA (email: [email protected]; [email protected]). a r X i v : . [ m a t h . O C ] O c t he problem of identifying critical nodes in a network graph leads to combinatorial op-timization problems. Thus, for large-scale networks any algorithm that solves the problemexactly is of high complexity. Motivated by this, we study a relaxation of the problemthrough a continuum approximation of the network to the spatial domain where the nodesare distributed. Literature review.
We first cite some works that present combinatorial approaches tothe problem of critical node identification. In [3–6], the authors investigate the problem ofidentifying nodes whose deletion minimizes some network connectivity metric. An alternativeapproach to improving network robustness involves incorporating redundancy in the networkby adding nodes and links, also called network augmentation [7]. In [8], the authors studythe problem of network design as a function of the comparative costs of augmentation anddefense against attack/failure.The approximation of large networks by weighted graphs over a continuum set of infinitecardinality appears in previous literature. In this way, in [9] large networks are approxi-mated by the so-called graphons, which result from the limit of convergent sequences of largedense graphs. Extending this idea to spatial networks, where the nodes are embedded ina domain Ω ∈ R N , the nodes can be thought to be indexed by their positions x ∈ Ω, andinteractions restricted between the nearest spatial neighbors. Combining these notions in thecontext of network consensus dynamics, the object of interest is the continuum counterpartof the graph Laplacian, the Laplace operator on the domain. Theoretical results concerningthe convergence of the graph Laplacian to the Laplace operator can be found in [10] and [11],which motivates the approach adopted in this paper.There have been severals attempts to investigate problems linking the shape of a domainwith the sequence of eigenvalues of the Laplace operator, for various boundary conditions,although those related to the critical subset identification are fewer in number. The work [12]contains an overview of the literature on extremum problems for eigenvalues of elliptic (e.g.Laplace) operators. In [13], the authors consider the problem of placing small holes in adomain to optimize the smallest Neumann eigenvalue of the Laplace operator (but withDirichlet boundary condition on the hole).
Statement of Contributions.
In this paper, we aim to study a critical node set identificationproblem for large-scale spatial networks with an associated weight-balanced Laplacian matrix.By considering a graph embedding technique, we reduce the problem to spatial networkswith uniformly distributed nodes and nearest-neighbors communication topologies. Then weconsider a special case of a hole-placement problem, which consists of identifying the optimallocation of the center of a ball in the domain that minimizes the smallest positive eigenvalueof the Laplace operator for the residual domain. With the help of the Min-max theorem,we formulate our objective as an infinite-dimensional, non-convex and nested optimizationproblem. This limits our goal at the outset to achieving convergence to a local optimum.Since the solution is hard to obtain analytically, we develop an algorithmic approach to suchproblem. First, we consider the inner optimization or eigenvalue problem, whose KKT pointsinclude the eigenvalues of the Laplace operator. We then provide a closed-form expressionfor the projected gradient flow in a Banach space for this problem that converges to theset of KKT points. Exploiting further the special properties of these dynamics, we provethat the only locally asymptotically stable equilibrium point for the dynamics is the secondeigenfunction of the Laplace operator. Moreover, since the other KKT points are saddleoints that are non-degenerate, we infer almost global asymptotic stability of the secondeigenfunction. Building on these results, we then design a novel hole-placement dynamicsfor the nested-optimization problem, and prove its local asymptotic stability to strict localminima. Finally, we provide a characterization of critical balls in the interior of the domain,and study the limiting case when its radius approaches zero. We conclude that the locationof such critical nodes is at the nodal set of the second eigenfunction of the Laplace operator,which has an intuitive geometric interpretation in some cases. A partial account of the resultsof this paper were presented without technical proofs in [14]. In addition to presenting the fulltechnical proofs, we present further analysis on the limiting case of hole-placement problemand include additional simulation examples.
Organization.
This paper is organized as follows. Sections 2 and 3 introduce some notationand preliminaries respectively. This is followed by the problem formulation in Section 4 andmain analysis in Section 5. We present some simulation results in Section 6 and concludewith the summary and future directions in Section 7.
We now introduce some basic notation used in the sequel. First, we denote by n the vectorof ones (cid:62) n = (1 , . . . , (cid:62) ∈ R n , for some n ∈ N . For a graph G , we denote by L ( G ) thegraph Laplacian and by λ ( L ( G )) the algebraic connectivity of the graph. The correspondingeigenvector, also called the Fiedler eigenvector, is denoted by v F . The open ball of radius r > x ∈ R N is represented by B r ( x ), and | Ω | denotes the Lebesgue measureof the set Ω ⊂ R N . The set of square-integrable functions on Ω is denoted by L (Ω). Inother words, L (Ω) = { f : Ω → R | (cid:82) Ω | f | dν < ∞} , where dv is the standard Lebesguemeasure. When clear from the context, we will denote (cid:82) Ω f dν simply as (cid:82) Ω f , for some f ∈ L (Ω), with a slight abuse of notation. For f, g ∈ L (Ω), we let (cid:104) f, g (cid:105) = (cid:82) Ω f gdν denote the inner product and (cid:107) f (cid:107) = (cid:104) f, f (cid:105) denote the corresponding induced norm. Wedenote by H (Ω) = { f ∈ L (Ω) | (cid:82) Ω |∇ f | dν < ∞} . For a bounded domain Ω, we denoteby ∂ Ω the boundary of Ω and by n the outward normal to the boundary. We also let S denote the Lebesgue measure on the boundary (where the integral of f on the boundaryis written as (cid:82) ∂ Ω f dS ). Let ∂ denote the partial differential operator. For a differentiablefunction F : Ω × Ω → R , we denote by ∂ F ( x , y ) (resp. ∂ F ( x , y )) the partial derivativeof F w.r.t. the first argument (resp. the second argument), evaluated at ( x , y ). Finally,given Ω ⊂ R N , ∆(Ω) represents the Laplace operator on the domain Ω (we omit Ω in ∆(Ω)when it is clear from context). In this section, we present the necessary background for setting up the critical node identi-fication problem addressed in this paper. We begin by explaining how we employ a graphembedding along with a continuum approximation to go from the graph Laplacian to theLaplace operator on the domain. Using the Min-max theorem, we are then able to character-ize the second eigenvalue of the Laplace operator corresponding to the algebraic connectivityof the graph. We finally point out to a connection to agreement algorithms in networkedsystems.et G = ( V, E ) be a weight-balanced directed graph such that | V | = n , and w ij bethe edge weight corresponding to ( i, j ) ∈ E . A map x : V → Ω ⊂ R N , is called a graphembedding ( N (cid:28) n and Ω bounded), if x i = x ( i ) ∈ R N is the (spatial) position assigned tonode i ∈ V , and the map x preserves some proximity measure on the graph G . There existsa vast literature on graph embeddings [15, 16], of which we adopt the notion of the structure-preserving embedding. Starting with the unweighted, undirected graph corresponding to G (where the weighted directed edges in G are replaced by unweighted undirected edges), astructure preserving embedding can be constructed such that any node j which is a neighborof i in the graph G is within a ball of radius h centered at at x i in the embedding. Oncethe graph is embedded in Ω ⊂ R N , we view the nodes V as having been sampled from anunderlying distribution µ ∈ P (Ω) (with density function ρ , such that dµ = ρ dvol). It isalways possible to obtain the weighted adjacency matrix W = [ w ij ] of the digraph G as thediscretization of a smooth weight function W : Ω × Ω → R ≥ , such that w ij = W ( x i , x j ).The weight function W encodes the weights and directionality of the edges, and since thenumber of nodes V is finite, such a smooth weight function always exists. Let ϕ : Ω → R be a real-valued function on Ω and φ d : V → R such that φ di = φ d ( i ) = ϕ ( x i ). We definethe W -weighted average variation in ϕ around a point x ∈ Ω, averaged over a ball B h ( x ) ofradius h > x as follows:1 µ ( B h ( x )) (cid:90) B h ( x ) W ( x, y )( ϕ ( y ) − ϕ ( x )) dµ ( y ) . We see next that the weighted Laplace operator on Ω can be obtained as the limit of a W -weighted average variation as h →
0. We first let w ( x ) = W ( x, x ) and ∇ w ( x ) = ( ∂ W + ∂ W )( x, x ), and we obtain the following by means of a Taylor expansion:lim h → ch µ ( B h ( x )) (cid:90) B h ( x ) W ( x, y )( ϕ ( y ) − ϕ ( x )) dµ ( y )= 1 ρ ∇ · ( wρ ∇ ϕ ) , where c is a constant. The graph Laplacian matrix L ( G ) corresponding to G can now beviewed as the discretization of the (negative) w -weighted Laplace operator − ρ ∇ · ( wρ ∇ ).Alternatively, the w -weighted Laplace operator can be viewed as an approximation of L ( G ),with closer approximations obtained as n = | V | → ∞ and h → L ( G ) by the Laplace operator on Ωrequires the specification of a boundary condition. This condition is obtained by observingthat n ∈ Null( L (cid:62) ( G )), that is, (cid:10) , L ( G ) φ d (cid:11) = (cid:62) n L ( G ) φ d = 0 for any φ d . In the continuoussetting, this translates into the Neumann boundary condition ∇ ϕ · n = 0 on ∂ Ω. Thiscan be seen from an application of the Divergence theorem, that is, (cid:68) , ρ ∇ · ( wρ ∇ ϕ ) (cid:69) = (cid:82) Ω 1 ρ ∇ · ( wρ ∇ ϕ ) dµ = (cid:82) ∂ Ω wρ ∇ ϕ · n dS = 0 (if ∇ ϕ · n = 0). Thus, the Neumann boundarycondition is imposed as the natural boundary condition here. Remark 1. (Problem reduction to uniformly spatially embedded graphs).
Basedon the previous considerations, and without loss of generality, in the following we focus onnetworks that are spatially embedded in an open bounded domain Ω according to a uniformdistribution (the distribution µ is uniform above) and such that the underlying graph is undi-rected and unweighted. Note that the following derivations are analogous for the case of aon-uniform µ and weight-balanced directed graph: all results carry through by keeping theweights w and ρ in the weighted Laplace operator. The Laplace operator ∆ with the Neumann boundary condition, has an infinite sequenceof eigenvalues 0 = µ ≤ µ ≤ . . . ≤ µ m ≤ . . . , whose corresponding eigenfunctions { ψ i } ∞ i =1 form an orthonormal basis for L (Ω), [17]. Using the Min-max theorem [17] for the operator∆, one can determine: µ (Ω) = inf ψ ∈{ ψ } ⊥ (cid:104) ψ, ∆ ψ (cid:105) L (Ω) (cid:104) ψ, ψ (cid:105) L (Ω) , (1)where { ψ } ⊥ = { ψ ∈ H (Ω) | ψ (cid:54) = 0 , (cid:82) Ω ψ ψ dν = 0 } , and ψ is constant, the eigenfunctioncorresponding to µ = 0. This implies { ψ } ⊥ = { ψ ∈ H (Ω) | (cid:82) Ω ψ dν = 0 } . Thus, usingthe Divergence theorem, applying the Neumann boundary condition, and normalizing thefunctions, we obtain an equivalent reformulation of (1) as: µ (Ω) = inf ψ ∈ H (Ω) , (cid:82) Ω ψdν =0 , (cid:82) Ω | ψ | dν =1 (cid:90) Ω |∇ ψ | dν. (2) Remark 2. (Connection to agreement algorithms).
The second eigenvalue is alsoof relevance to Laplacian-based agreement/consensus algorithms in networked systems, as itgoverns the convergence rate of these algorithms.
We define in this section the notion of criticality adopted in this manuscript. We define criticalnodes as those nodes in the graph whose removal results in the maximum deterioration inalgebraic connectivity for the residual network, making them the most crucial nodes to beprotected.More precisely, this amounts to identifying a set K ∗ ⊂ Ω of given measure | K ∗ | = c > µ (Ω \ K ∗ ) is an infimum. The problem of identifying the critical nodes, K ∗ , canbe formulated as: K ∗ ∈ arg inf K ⊂ Ω , | K | = c inf ψ ∈ H (Ω \ K ) , (cid:82) Ω \ K ψdν =0 , (cid:82) Ω \ K | ψ | dν =1 (cid:90) Ω \ K |∇ ψ | dν. We restrict the search to a class of subsets K = B r ( x ) = { y ∈ Ω | | y − x | < r } ⊂ Ω, openballs of radius r (such that | B r ( x ) | = c ). This reduces the search space to ˜Ω r = { x ∈ Ω | dist( x, ∂ Ω) > r } , and the problem is reformulated as: x ∗ ∈ arg inf x ∈ ˜Ω r inf ψ ∈ H (Ω \ B r ( x )) , (cid:82) Ω \ Br ( x ) ψdν =0 , (cid:82) Ω \ Br ( x ) | ψ | dν =1 (cid:90) Ω \ B r ( x ) |∇ ψ | dν. (3)which we refer to as the hole-placement problem in the sequel. emark 3. (Generalization using multiple balls). We note that any compact sub-set K ⊂ Ω can be covered by a finite number m of open balls of a given radius r , and witharbitrary precision (as r → and m → ∞ ). Given a finite collection { B r ( x i ) } mi =1 of openballs, we can then formulate the above optimization w.r.t. ( x , . . . , x m ) , the positions of the m open balls. For simplicity, we just focus on the one-ball case. Here, we present our main results and algorithms to determine the most critical nodes in thenetwork, in a functional optimization framework. To do this, we begin with the eigenvalueproblem (2) (which is the inner optimization problem in (3)) for D , a fixed domain, and designa projected gradient flow to converge to a local minimizer of the problem. This algorithmwill help us build subsequently the dynamics that can be employed to solve the full holeplacement problem (3) in an algorithmic manner. The analysis of the projected gradient flowwill also be instrumental in evaluating the properties of the second dynamics. µ (Ω) In what follows, we study the eigenvalue problem (2), characterize its critical points, con-struct and analyze a novel projected gradient flow to converge to the infimum. We writethe optimization problem (for the smallest positive eigenvalue of the Laplace operator on adomain D with a C , Lipschitz boundary) as:inf ψ ∈ H ( D ) (cid:90) D |∇ ψ | , s.t (cid:90) D | ψ | = 1 , (cid:90) D ψ = 0 , ∇ ψ · n = 0 on ∂D. Let S D = { ψ ∈ H ( D ) | (cid:82) D | ψ | = 1 , (cid:82) D ψ = 0 , ∇ ψ · n = 0 on ∂D } and J ( ψ ) = (cid:82) D |∇ ψ | . Wecan now express the above problem as inf ψ ∈S D J ( ψ ). Lemma 1. (Minimizer of J ( ψ ) ). The eigenfunctions of ∆( D ) are the critical points ofthe functional J ( ψ ) , and the second eigenfunction ψ of ∆( D ) is the only minimizer of thefunctional J ( ψ ) in S D . Moreover, the critical points of J ( ψ ) are non-degenerate, i.e., theHessian of J ( ψ ) is non-singular at the critical points. The content of this Lemma follows from the Min-max theorem [17]. We refer the reader tothe Appendix for an alternative proof of this lemma, as well as for the proofs of other resultscontained in this paper. We explicitly compute the analytical expression for the Hessian of theobjective function J ( ψ ) in the proof of Lemma 1, which allows us to infer the non-degeneracyof the saddle points of J ( ψ ) which is useful in establishing almost-global convergence of theprojected gradient flow we present below.We now provide a novel closed-form expression for a projected gradient flow to con-verge to the minimum value of J ( ψ ) in S D . For smooth one-parameter families of functions ψ ( t, x ) } t ∈ R ≥ (with x ∈ D ), the derivative of the objective functional J is given by: ddt [ J ( ψ ( t ))] = 2 (cid:90) D ∇ ψ · ∇ ( ∂ t ψ ) = − (cid:90) D ∂ t ψ (∆ ψ ) . We obtain a gradient flow by setting ∂ t ψ = ∆ ψ . We project this flow onto the tangent spaceof the set S D . For ψ ∈ S D , we require that (cid:104) ψ, ∂ t ψ (cid:105) = 0 and (cid:82) D ∂ t ψ = 0, which are satisfiedif (this will be shown in Proposition 1): ∂ t ψ = ∆ ψ − (cid:104) ∆ ψ, ψ (cid:105)(cid:107) ψ (cid:107) ψ = ∆ ψ − (cid:104) ∆ ψ, ψ (cid:105) ψ, since (cid:107) ψ (cid:107) = 1 for ψ ∈ S D . Further, using J ( ψ ) = − (cid:104) ∆ ψ, ψ (cid:105) , we get the projected gradientflow: ∂ t ψ = ∆ ψ + J ( ψ ) ψ. (4)The equilibria ψ ∗ of (4) satisfy ∆ ψ ∗ + J ( ψ ∗ ) ψ ∗ = 0 and the Neumann boundary condition ∇ ψ ∗ = 0 on ∂D . Clearly, J ( ψ ∗ ) is an eigenvalue, and so let µ ∗ = J ( ψ ∗ ). It is also clear thatthe equilibria of the projected gradient flow are also the critical points of the functional J over the set S D . Proposition 1. (Convergence of gradient flow).
The set S D is invariant with respectto the flow (4) , and the solutions to (4) in S D converge in an L sense to the set of equilibriaof (4) . Moreover, the only locally asymptotically stable equilibrium in S D for (4) is the secondeigenfunction ψ . Remark 4. (Implication of Proposition 1).
Proposition 1 states that we have globalconvergence to the set of isolated equilibria of the gradient flow (4) and that only the secondeigenfunction ψ is locally asymptotically stable among the set of isolated equilibria. Moreover,as seen in the proof of Lemma 1, we have that the other equilibria are saddle points of J ( ψ ) and are non-degenerate (the Hessian of J at these saddle points are non-singular). From thiswe deduce almost global asymptotic stability of the second eigenfuction ψ for the flow (4) ,and we therefore have convergence from almost all initial conditions, see [18] for an overviewof this property. We now consider the full optimization problem (3), which can be expressed as: x ∗ ∈ arg inf x ∈ ˜Ω r µ (Ω \ B r ( x )) Assumption 1. (Simplicity of the second eigenvalue).
We assume that the secondeigenvalue µ (Ω \ B r ( x )) is simple for any x ∈ ˜Ω . Remark 5. (Relaxing Assumption 1).
The assumption that the eigenvalue µ is simpleis ensures differentiability of µ (Ω \ B r ( x )) w.r.t. x . The eigenvalues of ∆(Ω \ B r ( x )) existas branches x (cid:55)→ µ (Ω \ B r ( x )) , which can then be ordered as µ ≤ µ ≤ . . . for any given x .The branches x (cid:55)→ µ (Ω \ B r ( x )) of eigenvalues are differentiable w.r.t. x (more generally.r.t. the perturbation of domains with Lipschitz boundaries [12]). The case of a non-simpleeigenvalue µ occurs when multiple branches intersect, for some x , at which point the or-dering of the branches may change and we lose differentiability of µ . This situation canhowever be mitigated by considering the subdifferential of µ in place of the gradient of µ .The dynamics presented later in the paper can be modified in this sense, and the analysiswould require further investigation on the regularity/lower-semicontinuity properties of thesesubdifferentials. We nevertheless avoid this problem through Assumption 1, which we leaveas future work. The following lemma allows for a characterization of the critical points of the functional µ in the interior of the domain. Lemma 2. (Characterization of critical ball).
The first-order condition for a criticalpoint x ∗ of the functional µ in the interior of the domain is given by: µ ∗ (cid:32)(cid:90) ∂B r ( x ∗ ) | ψ ∗ | n (cid:33) = (cid:90) ∂B r ( x ∗ ) |∇ ψ ∗ | n , (5) where ( µ ∗ , ψ ∗ ) is the second eigenpair such that µ ∗ (cid:52) = µ (Ω \ B r ( x ∗ )) . We now construct the gradient dynamics to converge to a critical point of µ in theinterior of the domain. Note that the function µ (Ω \ B r ( x )) is not known explicitly for ageneral domain Ω \ B r ( x ). We reformulate the optimization problem (3) as: x ∗ = arg inf ( x,ψ ) ∈ ˜Ω × Ψ( x ) (cid:90) Ω \ B r ( x ) |∇ ψ | dν, (6)where the set Ψ( x ) is defined as: Ψ( x ) = (cid:40) ψ ∈ H (Ω \ B r ( x )) (cid:12)(cid:12)(cid:12)(cid:12) (cid:90) Ω \ B r ( x ) ψ = 0 , (cid:90) Ω \ B r ( x ) | ψ | = 1 (cid:41) , (7) where arg indicates the first argument x in ( x, ψ ). We also define the set Ψ = ∪ x ∈ ˜Ω r Ψ( x ).We recall that ˜Ω r = { x ∈ Ω | dist( x, ∂ Ω) > r } . Now let { x ( t ) } t ∈ R ≥ be a smooth curvein ˜Ω r and { ψ ( t, y ) } t ∈ R ≥ (with y ∈ Ω \ B r ( x ( t )),) a smooth one-parameter family of functionson Ω \ B r ( x ( t )). Also, let ˜ n ( x ) be the normal to the boundary ∂ ˜Ω r at x ∈ ∂ ˜Ω r . We nowconsider the following hole-placement dynamics for our nested optimization problem: dxdt = (cid:40) v int , x ∈ int ˜Ω r v int − ( v int · ˜ n )˜ n , x ∈ ∂ ˜Ω r v int = − (cid:90) ∂B r ( x ) |∇ ψ | n + J ( ψ ) (cid:90) ∂B r ( x ) | ψ | n ,∂ t ψ = ∆ ψ + J ( ψ ) ψ + aψ + b, ∇ ψ · n = 0 , on ∂ Ω ∪ ∂B r ( x ) , (8)where a = − v · (cid:16)(cid:82) ∂B r ( x ) | ψ | n (cid:17) and b = − | Ω |− c v · (cid:16)(cid:82) ∂B r ( x ) ψ n (cid:17) , with c = | B r ( x ) | , forall x ∈ ˜Ω r . heorem 1. (Convergence of the hole placement dynamics). The set Ψ in (7) is invariant with respect to the dynamics (8) . The solutions to the dynamics (8) convergeto a critical point of the objective functional µ in (6) . A critical point of µ is locallyasymptotically stable with respect to the dynamics (8) only if it is a strict local minimum. Remark 6. (Implication of Theorem 1).
Theorem 1 states that we have convergence tothe equilibria of the hole-placement dynamics which are also critical points of µ (Ω \ B r ( x )) .In addition, we have that among the critical points of µ (Ω \ B r ( x )) , only the strict localminima are locally asymptotically stable. For almost global convergence to these strict localminima, we additionally require non-degeneracy of the saddle points of µ (Ω \ B r ( x )) (i.e.,that the Hessian is non-singular at the critical point), but this additional characterization isnot contained in our result. We now consider the following question: if an initial failure happens with the removal of anode, what is the most critical node? This is appropriately posed in the continuum setting asthe hole placement problem where the size of the hole is very small, i.e., as the radius r → f ( x ) = lim r → | ∂B r ( x ) | ∂∂r µ (Ω \ B r ( x )),which quantifies as a function of the hole position, the rate of deterioration of the metric asfailure begins to occur. Theorem 2. (Connection to the nodal set of eigenfunction).
In the limit r → forthe radius of the hole, the hole-placement problem reduces to finding the minima x ∗ ∈ Ω ofthe function: f ( x ) = µ Ω2 | ψ Ω2 ( x ) | − |∇ ψ Ω2 ( x ) | , where µ Ω2 , ψ Ω2 ( x ) is the second eigenpair of the domain Ω . Moreover, a point x ∗ ∈ Ω is alocal minimizer of f , if ψ Ω2 ( x ∗ ) = 0 and the family of level sets of ψ Ω2 is locally flat at x ∗ .In other words, the nodal points of ψ Ω2 where the family of level sets of ψ Ω2 is locally flat arelocal minimizers of f . Remark 7. (Geometry of nodal sets).
The nodal sets of Neumann eigenfunctions havebeen extensively investigated [19]. It is known that if the domain is symmetric about a subset,then it contains the nodal set of ψ . The nodal set for the second eigenfunction ψ Ω2 dividesthe domain Ω into no more than two regions Ω a and Ω b . Now, µ Ω2 is the first eigenvalue λ ofthe Laplacian for Ω a and Ω b , with Neumann boundary condition on ∂ Ω ∩ ∂ Ω a and Dirichletboundary condition on ∂ Ω a ∩ ∂ Ω b . Remark 8. (Implication for networks).
Theorem 2 can be used to provide new insighton where the most critical nodes in a network with a finite number of nodes are located, viaa continuum approximation. This is based on the fact that the entries v Fi of the Fiedlereigenvector v F of the finite graph embedded in Ω can be approximated by the value of theeigenfunction ψ Ω2 at the location x i of the node i . That is, v Fi ≈ ψ Ω2 ( x i ) . Then the mostcritical nodes in the network correspond to the zero entries of the Fiedler eigenvector. TheFiedler eigenvector, however, does not necessarily contain zero entries for general finite graphs(this situation improves with the size of the graph), in which case we may expect the criticalnodes to be concentrated at the entries of lowest magnitude. This is a heuristic obtained fromthe fact that ψ Ω2 is smooth and that ψ Ω2 more closely approximates v F as n → ∞ . Simulation results
In this section, we present some numerical simulation results that can illustrate the conceptsand algorithms of the previous sections.First, we consider a disk-shaped domain Ω of unit radius, and the placement of a hole B of radius of 0 . µ for the residual domain Ω \ B as afunction of h (distance between the center of the disk and the center of the hole). Sincethe hole is of radius 0 . h ∈ [0 , . Figure 1: µ as a function of h for a disk-shaped domain.from Figure 1 that the second (also the smallest positive) eigenvalue of the Laplace operatorfor a disk-shaped domain with a hole increases with the distance between the centers of thedomain and the hole, but also appears to decrease as the hole approaches close to the domainboundary (around h = 0 .
85 units). Moreover, µ as a function of h appears to be a convexin the interval h ∈ [0 , .
85] and concave for h ∈ (0 . , . x (the center of the hole) as theslow-scale variable and ψ the fast-scale variable. We first consider the case of the disk-shapeddomain, that is, the dynamics (8) corresponds to hole placement for the disk-shaped domainto minimize µ of the residual domain.Figure 2 is a plot of x ( t ), the path of the center of the hole, on the spatial domain, fortwo different initial conditions x (0) = (0 . , .
5) and x (0) = ( − . , − . x ( t ), the path of the center of the hole (from the dynamics (8)) for aconvex polygonal spatial domain. The final location of the hole is also indicated in the figure.Figure 4 contains the results for a non-convex polygonal domain. The outer polygon isthe spatial domain Ω, while the inner polygon is the domain ˜Ω (the set of allowed positionsfor the center of the hole). The heatmap shows the value of µ of the residual domain (whichwas obtained by first sampling the domain uniformly at random at the points indicated bythe tiny circles, placing the hole at those points, computing µ of the residual domain, andthen interpolating to obtain the plot). The paths of the center of the hole x ( t ) (from the Figure 2: Path of the center of the hole, x ( t ) from two different initial conditions x (0) =(0 . , .
5) and x (0) = ( − . , − . -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1-1-0.8-0.6-0.4-0.200.20.40.60.81 Figure 3: Path of the center of the hole, x ( t ) from an initial condition x (0) = (0 . , − .
5) fora convex polygonal domain.dynamics (8)) from different initial conditions are also plotted. The paths do not all convergeto the same point in this case, but to a broader region (the darker region in the heatmap),which possibly contains more than one local minimum x ∗ . -1-0.8-0.6-0.4-0.200.20.40.60.81-1-0.8-0.6-0.4-0.200.20.40.60.81 Figure 4: Paths of the center of the hole, x ( t ) from different initial conditions.In Figure 5, we present a numerical validation of the discussion in Remark 8. Wefirst generated a random connected graph G with 50 nodes. We then computed the al-gebraic connectivities of the residual graphs obtained by the removal of one node from theraph λ ( L ( G \ { i } )), for each node, plotting it against the corresponding entry of the Fiedlereigenvector v Fi (the eigenvector corresponding to the second eigenvalue of the Laplacian, oralgebraic connectivity) of the original graph G . From the discussion in Remark 8, we expectthat the local minima of λ ( L ( G \ { i } )) are concentrated around nodes corresponding to theentries of the Fiedler eigenvector of lowest magnitude, which is illustrated in the figure. Wenote that in the corresponding hole-placement problem, the nodal sets of the second eigen-function ψ Ω2 are only the local minimizers of f ( x ) = µ Ω2 | ψ Ω2 ( x ) | − |∇ ψ Ω2 ( x ) | . We thereby donot expect all the zero entries of the Fiedler eigenvector to correspond necessarily to globalminimizers. However, the figure shows that the global minimum is indeed concentratedaround nodes corresponding to the entries of the Fiedler eigenvector of lowest magnitude.Figure 5: Plot of algebraic connectivity of residual network with the removal of one node vs.its corresponding entry in the Fiedler eigenvector, for a network with 50 nodes. In this paper, we studied the problem of identifying the critical nodes for consensus in large-scale spatial networks. We began by making a functional approximation of the Laplacianmatrix of the graph to the Laplace operator on the domain. In addition to being a naturalstep in the large- N limit, the real advantage of the approximation is that it does not concealthe geometry of the problem, which is important for spatial networks such as swarms andsensor networks. As a starting point, we analyzed the removal of balls of given measure fromthe domain. In future work, we would like to generalize the results to arbitrary sets overdomains with a non uniform distribution of nodes. Further generalization of the analysisrelaxing Assumption 1, as outlined in Remark 5, is also left for future work. We note thatthe proposed gradient dynamics were centralized in nature, the problem of distributed criticalnode set identification is also of interest and left for future work. Appendix
Proof. (Proof of Lemma 1).
The first variation of the Lagrangian L ( ψ, µ, λ ) = J ( ψ ) + µ (cid:0) − (cid:82) D | ψ | (cid:1) + λ (cid:82) D ψ , at a critical point ψ ∗ is zero (where (cid:82) D | ψ | = 1 and (cid:82) D ψ = 0 arethe constraints, as ψ ∈ S D and the Neumann boundary condition is assumed implicitly.)Thus, for any δψ ∈ T ψ ∗ S D the tangent space of S D at ψ ∗ , we have (cid:68) δLδψ , δψ (cid:69) ( ψ ∗ , µ ∗ , λ ∗ ) =2 (cid:82) D ∇ ψ ∗ ·∇ ( δψ ) − µ ∗ (cid:82) D ψ ∗ δψ + λ ∗ (cid:82) D δψ = − (cid:82) D (∆ ψ ∗ + µ ∗ ψ ∗ − λ ∗ ) δψ = 0, for any δψ (notethat the Neumann boundary condition was used in obtaining the equation.) Additionally,we also have (cid:68) ∂L∂µ , δµ (cid:69) ( ψ ∗ , µ ∗ , λ ∗ ) = 1 − (cid:82) D | ψ ∗ | = 0, and (cid:10) ∂L∂λ , δλ (cid:11) ( ψ ∗ , µ ∗ , λ ∗ ) = (cid:82) D ψ ∗ = 0.Thus, the critical points of the objective functional ψ ∗ ∈ S D are characterized by:∆ ψ ∗ + µ ∗ ψ ∗ − λ ∗ = 0 . Integrating the previous equation over D and using the Neumann boundary condition, weobtain λ ∗ = 0. Therefore, the critical points ψ ∗ satisfy:∆ ψ ∗ + µ ∗ ψ ∗ = 0 . (9)Let ψ ( x, (cid:15), η ), x ∈ D , be a smooth two-parameter family of functions in S D with (cid:82) D ψ ( x, (cid:15), η ) =0 for all (cid:15) and η . The first variation of J at (cid:15) = 0, η = 0 is given by: δJδ(cid:15) (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 ,η =0 ( ψ ) = 2 (cid:90) D ∇ ψ · ∂ (cid:15) ∇ ψ = 2 (cid:90) D ∇ ψ · ∇ ( ∂ (cid:15) ψ ) . We let ∂ (cid:15) ψ | (cid:15) =0 ,η =0 = X and ∂ η ψ | (cid:15) =0 ,η =0 = Y . The second variation of J at (cid:15) = 0, η = 0 isgiven by: δ Jδηδ(cid:15) ( X, Y ) = 2 (cid:90) D ∇ ( ∂ η ψ ) · ∇ ( ∂ (cid:15) ψ ) + 2 (cid:90) D ∇ ψ · ∇ ( ∂ η(cid:15) ψ )= 2 (cid:90) D ∇ ( ∂ η ψ ) · ∇ ( ∂ (cid:15) ψ ) − (cid:90) D ∆ ψ ( ∂ η(cid:15) ψ )= 2 (cid:90) D ∇ X · ∇ Y − (cid:90) D ∆ ψ ( ∂ η(cid:15) ψ ) . Evaluating the second variation at a critical point ψ ( x, ,
0) = ψ ∗ , and from (9), we obtain: δ Jδηδ(cid:15) ( X, Y ) = 2 (cid:90) D ∇ X · ∇ Y + 2 µ ∗ (cid:90) D ψ ∗ ( ∂ η(cid:15) ψ ∗ ) . (10)Since ψ ( x, (cid:15), η ) is a smooth two-parameter family of functions in S D , we have (cid:82) D | ψ ( x, (cid:15), η ) | =1 for all (cid:15), η , which implies that (cid:82) D ψ ( ∂ (cid:15) ψ ) = 0 and (cid:82) D ∂ η ψ∂ (cid:15) ψ + (cid:82) D ψ ( ∂ η(cid:15) ψ ) = (cid:82) D XY + (cid:82) D ψ ( ∂ η(cid:15) ψ ) = 0. Substituting in (10), we obtain: δ Jδηδ(cid:15) ( X, Y ) = 2 (cid:90) D ∇ X · ∇ Y − µ ∗ (cid:90) D XY. n particular, for X (cid:54) = 0, this implies: δ Jδηδ(cid:15) ( X, X ) = 2 (cid:90) D |∇ X | − µ ∗ (cid:90) D | X | = 2 (cid:18)(cid:90) D | X | (cid:19) (cid:18) (cid:82) D |∇ X | (cid:82) D | X | − µ ∗ (cid:19) . (11)We also have that (cid:82) D ψ ( x, (cid:15), η ) = 0, which leads to (cid:82) D ∂ (cid:15) ψ = (cid:82) D X = 0. From (2), we havethat inf (cid:82) D X =0 (cid:82) D |∇ X | (cid:82) D | X | = µ , which implies that if µ ∗ > µ in (11), by the definition ofinfimum, there exists an X such that δ Jδηδ(cid:15) (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 ,η =0 ( X, X ) <
0. Therefore, the only criticalpoint for which δ Jδηδ(cid:15) (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 ,η =0 ( X, X ) ≥ ψ ∗ = ψ . Note that, forthis case, δ Jδηδ(cid:15) (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 ,η =0 ( X, X ) = 0 if and only if X = kψ . Since (cid:82) D ψ X = 0, it must bethat k = 0, and therefore X = 0. Thus, for all X (cid:54) = 0, δ Jδηδ(cid:15) (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 ,η =0 ( X, X ) > ψ ∗ = ψ .Therefore, the second eigenfunction ψ is the only minimizer of the functional J ( ψ ) in S D .It further follows from the above argument that the Hessian δ Jδηδ(cid:15) (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 ,η =0 is non-degenerate(or non-singular) at the critical points of J ( ψ ), that is, δ Jδηδ(cid:15) (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 ,η =0 ( X, X ) = 0 at the criticalpoints of J ( ψ ) if and only if X = 0. Proof. (Proof of Proposition 1).
Recall that S D = { ψ ∈ H ( D ) | (cid:82) D | ψ | = 1 , (cid:82) D ψ = 0 } .Therefore, for a smooth one-parameter family { ψ ( t, x ) } t ∈ R ≥ , (with x ∈ D ) to be in S D , weneed to prove that (cid:82) D ψ ∂ t ψ = 0 and (cid:82) D ∂ t ψ = 0, assuming that the initial condition is in S D . (Note that it will later be shown that ddt (cid:107)∇ ψ (cid:107) ≤
0, thus ψ ( t, · ) ∈ H ( D ) for all t ≥ ψ (0 , · ) ∈ S D ).From Equation (4), we have (cid:82) D ψ ∂ t ψ = (cid:82) D ψ (∆ ψ + J ( ψ ) ψ ). Using the Divergence theoremand the Neumann boundary condition on ∂ Ω, we get (cid:82) D ψ ∂ t ψ = − (cid:82) D |∇ ψ | + J ( ψ ) (cid:82) D | ψ | =0 (since J ( ψ ) = (cid:82) D |∇ ψ | and (cid:82) D | ψ | = 1).We also have (cid:82) D ∂ t ψ = (cid:82) D ∆ ψ + J ( ψ ) (cid:82) D ψ = (cid:82) D ∇ ψ · n + J ( ψ ) (cid:82) D ψ = 0 because of theNeumann boundary condition, ∇ ψ · n = 0 on ∂D , and (cid:82) D ψ = 0.Let ψ ( t, x ) be a solution of (4) in S D , with t ∈ R ≥ , x ∈ D , such that ψ (0 , x ) ∈ S D . Wealso have (cid:82) D | ψ | = 1, for all t ≥
0. Thus, J ( ψ ) = (cid:82) D |∇ ψ | = (cid:82) D |∇ ψ | (cid:82) D | ψ | . The time derivativef J is given by: ddt J = 2 (cid:82) D | ψ | (cid:90) D ∇ ψ · ∇ ∂ t ψ − (cid:82) D |∇ ψ | (cid:0)(cid:82) D | ψ | (cid:1) (cid:90) D ψ∂ t ψ = − (cid:90) D ∆ ψ ∂ t ψ − J ( ψ ) (cid:90) D ψ∂ t ψ = − (cid:90) D (∆ ψ + J ( ψ ) ψ ) ∂ t ψ = − (cid:90) D | ∆ ψ + J ( ψ ) ψ | ≤ . We have that J ≥ ddt J ≤
0. We also have S D ⊂ H ( D ), D a bounded, open subsetof R N with ∂D being C . Thus by the Rellich-Kondrachov Compactness Theorem [17], weget that the orbit ψ is precompact in L ( D ). Therefore, by the LaSalle invariance principlefor infinite dimensional spaces [20], the solutions converge in an L sense to largest invariantset contained in { ψ ∗ ∈ S D | ∆ ψ ∗ + J ( ψ ∗ ) ψ ∗ = 0 } , the set of equilibria of (4).In what follows we use the shorthand ∂ t ψ = F ( ψ ), where F ( ψ ∗ ) = 0, for the dynamics (4).We consider perturbations δψ ∈ T D along the tangent space of S D at ψ ∗ (also note that ψ ∗ is an eigenfunction). Thus (cid:82) D δψ = 0 and (cid:82) D ψ ∗ δψ = 0. We have: F ( ψ ∗ + δψ ) = ∆( ψ ∗ + δψ ) + J ( ψ ∗ + δψ )( ψ ∗ + δψ ) . Since ψ ∗ is a critical point of J ( ψ ) it holds that J ( ψ ∗ + δψ ) = J ( ψ ∗ ) + O ( (cid:107) δψ (cid:107) ) = µ ∗ + O ( (cid:107) δψ (cid:107) ). Thus, up to first-order we have that: F ( ψ ∗ + δψ ) = ∆( ψ ∗ + δψ ) + J ( ψ ∗ + δψ )( ψ ∗ + δψ )= − µ ∗ ψ ∗ + ∆( δψ ) + µ ∗ ψ ∗ + µ ∗ δψ = ∆( δψ ) + µ ∗ δψ. Therefore, we have ∂ t ( δψ ) = ∆( δψ ) + µ ∗ δψ . Expressing δψ ( t ) = (cid:80) ∞ i =2 α i ( t ) ψ i , where ψ i arethe eigenfunctions which form an orthonormal basis for T D , we have that: ∂ t ( δψ ) = ∞ (cid:88) i =2 ddt α i ( t ) ψ i = ∆( δψ ) + µ ∗ δψ = ∞ (cid:88) i =2 α i ( t )( − µ i + µ ∗ ) ψ i , which implies that δψ ( t ) = (cid:80) ∞ i =2 e ( µ ∗ − µ i ) t α i (0) ψ i . (Note that, from orthogonality, the previ-ous equality leads to ddt α i ( t ) = α i ( t )( − µ i + µ ∗ ), for each i .) We claim that the latter convergesto δψ = 0 for all initial conditions δψ (0) ∈ T D at ψ ∗ if and only if µ ∗ = µ (correspond-ingly, ψ ∗ = ψ ). To see this, first observe that, if µ ∗ = µ (correspondingly, ψ ∗ = ψ ), wehave (cid:82) D ψ δψ (0) = 0 (since δψ ∈ T D at ψ ∗ = ψ ), which implies that α (0) = α ( t ) = 0.Hence δψ ( t ) = (cid:80) ∞ i =3 e µ − µ i α i (0) ψ i and the exponent µ − µ i < i ≥
3. Conversely,if δψ ( t ) = (cid:80) ∞ i =2 e ( µ ∗ − µ i ) t α i (0) ψ i converges to δψ = 0 for all initial conditions δψ (0) ∈ T D at ψ ∗ , and ψ ∗ = ψ i for some i ∈ { , , . . . } . We have that α i (0) = α i ( t ) = 0 (from orthogo-nality), and that δψ ( t ) = (cid:80) ∞ j =2 ,j (cid:54) = i e ( µ i − µ j ) t α j (0) ψ j , which converges to δψ = 0 only if i = 2.Therefore, the second eigenfunction ψ is the only locally asymptotically stable equilibriumin S D for the projected gradient flow. roof. (Proof of Lemma 2). Let x ( (cid:15) ) for (cid:15) ∈ R be a smooth curve contained in ˜Ω r . Let ψ (cid:15) be the second eigenfunction of the Laplace operator with Neumann boundary condition in thedomain Ω \ B r ( x ( (cid:15) )). Thus, we have µ (cid:15) = (cid:82) Ω (cid:15) |∇ ψ (cid:15) | , where Ω (cid:15) = Ω \ B r ( x ( (cid:15) )) and (cid:107) ψ (cid:15) (cid:107) Ω (cid:15) =1 .The derivative dµ (cid:15) d(cid:15) is given by: dµ (cid:15) d(cid:15) = dd(cid:15) (cid:90) Ω (cid:15) |∇ ψ (cid:15) | = 2 (cid:90) Ω (cid:15) ∇ ψ (cid:15) · ∇ (cid:18) ∂ψ (cid:15) ∂(cid:15) (cid:19) + (cid:90) ∂ Ω (cid:15) |∇ ψ (cid:15) | v · n , (12) where v = dx ( (cid:15) ) d(cid:15) , is constant on ∂B r ( x (cid:15) ). Equation (12) becomes: dµ (cid:15) d(cid:15) = 2 (cid:90) Ω (cid:15) ∇ ψ (cid:15) · ∇ (cid:18) ∂ψ (cid:15) ∂(cid:15) (cid:19) + v · (cid:32)(cid:90) ∂B r ( x (cid:15) ) |∇ ψ (cid:15) | n (cid:33) = − (cid:90) Ω (cid:15) ∂ψ (cid:15) ∂(cid:15) ∆ ψ (cid:15) + v · (cid:32)(cid:90) ∂B r ( x (cid:15) ) |∇ ψ (cid:15) | n (cid:33) = 2 (cid:90) Ω (cid:15) µ (cid:15) ψ (cid:15) ∂ψ (cid:15) ∂(cid:15) + v · (cid:32)(cid:90) ∂B r ( x (cid:15) ) |∇ ψ (cid:15) | n (cid:33) (13)= µ (cid:15) dd(cid:15) (cid:18)(cid:90) Ω (cid:15) | ψ (cid:15) | (cid:19) − µ (cid:15) v · (cid:32)(cid:90) ∂B r ( x (cid:15) ) | ψ (cid:15) | n (cid:33) + v · (cid:32)(cid:90) ∂B r ( x (cid:15) ) |∇ ψ (cid:15) | n (cid:33) = − µ (cid:15) v · (cid:32)(cid:90) ∂B r ( x (cid:15) ) | ψ (cid:15) | n (cid:33) + v · (cid:32)(cid:90) ∂B r ( x (cid:15) ) |∇ ψ (cid:15) | n (cid:33) , since (cid:82) Ω (cid:15) | ψ (cid:15) | = 1 for all (cid:15) ∈ R , which implies that dd(cid:15) (cid:16)(cid:82) Ω (cid:15) | ψ (cid:15) | (cid:17) = 0. Let x (0) = x ∗ ∈ ˜Ωbe a critical point of µ ( x ), such that µ ( x ∗ ) = µ ∗ , with ψ ∗ being the second eigenfunction.Thus we have dµ (cid:15) d(cid:15) (cid:12)(cid:12) (cid:15) =0 = 0 for all v , which implies that: µ ∗ (cid:32)(cid:90) ∂B r ( x ∗ ) | ψ ∗ | n (cid:33) = (cid:90) ∂B r ( x ∗ ) |∇ ψ ∗ | n . This is the first-order condition for critical points of µ in the interior of the domain. Proof. (Proof of Theorem 1).
Let { x ( t ) , ψ ( t, y ) } t ∈ R ≥ (with y ∈ Ω \ B r ( x ( t )),) be aone-parameter family of functions that is a solution to the dynamics (8), and let ψ (0 , · ) ∈ Ψ( x (0)). To prove the invariance of Ψ( x ( t )), we need to show that ddt (cid:16)(cid:82) Ω \ B r ( x ( t )) | ψ | (cid:17) = 0and ddt (cid:16)(cid:82) Ω \ B r ( x ( t )) ψ (cid:17) = 0 (Note that it will later be shown that ddt (cid:107)∇ ψ (cid:107) ≤
0, thus ψ ( t, · ) ∈ (Ω) for all t ≥ ψ (0 , · ) ∈ Ψ). From (8), we have (with Ω( t ) = Ω \ B r ( x ( t ))): ddt (cid:32)(cid:90) Ω( t ) | ψ | (cid:33) = 2 (cid:90) Ω( t ) ψ ∂ t ψ + v · (cid:90) ∂B r ( x ( t )) | ψ | n = 2 (cid:90) Ω( t ) ψ ∆ ψ + 2 J ( ψ ) (cid:90) Ω( t ) | ψ | + 2 a ( t ) × (cid:90) Ω( t ) | ψ | + 2 b ( t ) (cid:90) Ω( t ) ψ + v · (cid:90) ∂B r ( x ( t )) | ψ | n = − (cid:90) Ω( t ) |∇ ψ | + 2 J ( ψ ) + 2 a ( t )+ v · (cid:90) ∂B r ( x ( t )) | ψ | n = 0 , because J ( ψ ) = (cid:82) Ω( t ) |∇ ψ | , (cid:82) Ω( t ) | ψ | = 1 and (cid:82) Ω( t ) ψ = 0 (since ψ ( t, · ) ∈ Ψ( x ( t )).) We alsohave: ddt (cid:32)(cid:90) Ω( t ) ψ (cid:33) = (cid:90) Ω( t ) ∂ t ψ + v · (cid:90) ∂B r ( x ( t )) ψ n = (cid:90) Ω( t ) ∆ ψ + J ( ψ ) (cid:90) Ω( t ) ψ + a ( t ) (cid:90) Ω( t ) ψ + b ( | Ω | − c )+ v · (cid:90) ∂B r ( x ( t )) ψ n = 0 . Since we also have that ψ (0 , · ) ∈ Ψ( x (0)), we conclude that the set Ψ is invariant with respectto the dynamics (8).Let { x ( t ) , ψ ( t, y ) } t ∈ R ≥ (with y ∈ Ω \ B r ( x ( t ))), be a one-parameter family of functionsthat is a solution to the dynamics (8), and let ψ ( t, · ) ∈ Ψ( x ( t )) for all t ∈ R ≥ (this assumptionis justified by the invariance of Ψ). We have J ( ψ ) = (cid:82) Ω( t ) |∇ ψ | = (cid:82) Ω( t ) |∇ ψ | (cid:82) Ω( t ) | ψ | ≥ ψ ( t, · ) ∈ Ψ( x ( t )) (since (cid:82) Ω( t ) | ψ | = 1). Now: ddt J = 2 (cid:90) Ω( t ) ∇ ψ · ∇ ∂ t ψ + v · (cid:90) ∂B r ( x ( t )) |∇ ψ | n − J ( ψ ) (cid:90) Ω( t ) ψ ∂ t ψ − J ( ψ ) v · (cid:90) ∂B r ( x ( t )) | ψ | n = − (cid:90) Ω( t ) | ∆ ψ + J ( ψ ) ψ | − v · v int ≤ , where we have used (8) to obtain the second equality. By the Rellich-Kondrachov Compact-ness Theorem [17], we see that the orbit ψ is precompact in L (Ω). Thus, by the invariancerinciple [20], the solutions { x ( t ) , ψ ( t, y ) } t ∈ R ≥ (with y ∈ Ω \ B r ( x ( t ))), converge to x ∗ , ψ ∗ (theconvergence ψ ( t, · ) → ψ ∗ , is in the sense of L ) such that v = 0 and ∆ ψ ∗ + J ( ψ ∗ ) ψ ∗ = 0. Wealready have that the only asymptotically stable case is when ψ ∗ = ψ ∗ (the second eigenfunc-tion corresponding to Ω \ B r ( x ∗ )), which implies that J ( ψ ∗ ) = J ( ψ ∗ ) = µ ∗ . And v = 0 impliesthat (cid:82) ∂B r ( x ∗ ) |∇ ψ ∗ | n = µ ∗ (cid:82) ∂B r ( x ∗ ) | ψ ∗ | n , the critical point of the functional µ from (5).Consider perturbations δx and δψ , about an equilibrium ( x ∗ , ψ ∗ ) such that x ∗ + δx ∈ ˜Ωand ˜ ψ = ψ ∗ + δψ ∈ Ψ( x ∗ + δx ) is the second eigenfunction of the domain Ω \ B r ( x ∗ + δx ).In other words, we consider perturbations purely in x to investigate the local asymptoticstability of the critical points of µ ( x ). The dynamics in x in this case, referring to (8), aregiven by: ddt ( x ∗ + δx ) = − (cid:90) ∂B r ( x ∗ + δx ) (cid:16) |∇ ˜ ψ | − ˜ µ | ˜ ψ | (cid:17) n . This can be reduced to: ddt ( δx ) = − ∂∂x (cid:12)(cid:12)(cid:12)(cid:12) x = x ∗ (cid:32)(cid:90) ∂B r ( x ) (cid:16) |∇ ˜ ψ | − ˜ µ | ˜ ψ | (cid:17) n (cid:33) δx. (14)From Equation (13), we recognize that (cid:82) ∂B r ( x ) (cid:16) |∇ ˜ ψ | − ˜ µ | ˜ ψ | (cid:17) n = ∂µ ∂x . Therefore, thelinearized dynamics reduces to: ddt ( δx ) = − ∂ µ ∂x (cid:12)(cid:12)(cid:12)(cid:12) x = x ∗ δx, where ∂ µ ∂x (cid:12)(cid:12)(cid:12)(cid:12) x = x ∗ is the Hessian of µ at x = x ∗ . Therefore, we have that the linearizeddynamics is asymptotically stable if and only if the Hessian of µ is positive definite, in otherwords, if and only if x ∗ is a strict local minimum of µ . Therefore, the necessary conditionfor the local asymptotic stability of the primal-dual dynamics at a critical point of µ is thatit is a strict local minimum. Proof. (Proof of Theorem 2).
Let r : R → R ≥ with r (0) = 0 be a smooth non-negativefunction. Let Ω( t ) = Ω \ B r ( t ) ( x ) for some x ∈ Ω ⊂ R N , be a one parameter family ofspatial domains such that Ω(0) = Ω. Let µ ( t ) be the second eigenvalue of the domain Ω( t )and ψ ( t, · ) the corresponding normalized eigenfunction (we assume that the family of spatialdomains Ω( t ) have simple eigenvalues). Thus, we have µ ( t ) = (cid:82) x ∈ Ω( t ) |∇ ψ ( t, x ) | . From [21],we have that µ ( t ) and ψ are real-analytic locally at t = 0. Thus, for small τ >
0, we have: µ ( τ ) = µ (0) + ddt µ (cid:12)(cid:12)(cid:12)(cid:12) t =0 τ + . . .ψ ( τ, x ) = ψ (0 , x ) + ∂ t ψ ( t, x ) (cid:12)(cid:12)(cid:12)(cid:12) t =0 τ + . . . (15)We note that µ (0) and ψ (0 , · ) are the second eigenpair corresponding to Ω. At a given t >
0, let the deformation of the domain be characterized by v = − (cid:15) n , the velocity of pointsn the boundary of the hole, B r ( t ) ( x ), where n is the normal to the domain Ω( t ) on theboundary of B r ( t ) ( x ), and (cid:15) > ddt µ = ddt (cid:90) x ∈ Ω( t ) |∇ ψ ( t, x ) | = 2 (cid:90) Ω( t ) ∇ ψ ∇ ∂ t ψ + (cid:90) ∂B r ( t ) ( x ) |∇ ψ | v · n = − (cid:90) Ω( t ) ∆ ψ ∂ t ψ + (cid:90) ∂B r ( t ) ( x ) |∇ ψ | v · n = 2 µ ( t ) (cid:90) Ω( t ) ψ ∂ t ψ + (cid:90) ∂B r ( t ) ( x ) |∇ ψ | v · n = µ ( t ) (cid:32) ddt (cid:90) Ω( t ) | ψ | − (cid:90) ∂B r ( t ) ( x ) | ψ | v · n (cid:33) + (cid:90) ∂B r ( t ) ( x ) |∇ ψ | v · n = µ ( t ) (cid:15) (cid:90) ∂B r ( t ) ( x ) | ψ | − (cid:15) (cid:90) ∂B r ( t ) ( x ) |∇ ψ | , since (cid:82) Ω( t ) | ψ | = 1, for all t . For small τ >
0, we then substitute from (15) in the aboveequation, to obtain: ddt µ = (cid:15) (cid:18) µ (0) + ddt µ (cid:12)(cid:12)(cid:12)(cid:12) t =0 τ + . . . (cid:19) × (cid:90) y ∈ ∂B r ( τ ) ( x ) | ψ (0 , y ) + ∂ t ψ ( t, y ) (cid:12)(cid:12)(cid:12)(cid:12) t =0 τ + . . . | − (cid:15) (cid:90) y ∈ ∂B r ( τ ) ( x ) |∇ ( ψ (0 , y ) + ∂ t ψ ( t, y ) (cid:12)(cid:12)(cid:12)(cid:12) t =0 τ + . . . ) | = µ (0) (cid:15) (cid:90) y ∈ ∂B r ( τ ) ( x ) | ψ (0 , y ) | − (cid:15) (cid:90) y ∈ ∂B r ( τ ) ( x ) |∇ ψ (0 , y ) | + O ( τ )= µ (0) S N − r ( τ ) N − (cid:15) | ψ (0 , x ) | − S N − r ( τ ) N − (cid:15) |∇ ψ (0 , x ) | + O ( r ( τ ) N − τ ) , where S N is the surface area of th unit N -sphere. Now, given that v = − (cid:15) n , we have r ( τ ) = (cid:15)τ , and therefore: ddt µ = µ (0) S N − (cid:15) N τ N − | ψ (0 , x ) | − S N − (cid:15) N τ N − |∇ ψ (0 , x ) | + O ( τ N ) . Substituting for ddt µ from the above equation into µ ( τ ) = µ (0)+ ddt µ (cid:12)(cid:12) ¯ τ τ (where ¯ τ ∈ [0 , τ ]),e get: µ ( τ ) = µ (0) + S N − (cid:15) N ¯ τ N − τ (cid:0) µ (0) | ψ (0 , x ) | − |∇ ψ (0 , x ) | (cid:1) + O (¯ τ N τ ) ≤ µ (0) + S N − (cid:15) N τ N (cid:0) µ (0) | ψ (0 , x ) | − |∇ ψ (0 , x ) | (cid:1) + O ( τ N +1 ) ≈ µ (0) + c ( τ ) (cid:0) µ (0) | ψ (0 , x ) | − |∇ ψ (0 , x ) | (cid:1) , where we have ignored the O ( τ N +1 ) term in the final expression. We also have r ( τ ) = (cid:15)τ , andtherefore the above can also be written as µ ( r ) ≈ µ (0)+ c ( r ) (cid:0) µ (0) | ψ (0 , x ) | − |∇ ψ (0 , x ) | (cid:1) as a function of the radius of the hole. We also note that the function (cid:0) µ (0) | ψ (0 , x ) | − |∇ ψ (0 , x ) | (cid:1) =lim r → | ∂B r ( x ) | ∂∂r µ (Ω \ B r ( x )).We now show that the local minima of f ( x ) = µ Ω2 | ψ Ω2 | − |∇ ψ Ω2 | occur along the nodal setof ψ Ω2 , that is, in the set { x ∈ Ω | ψ Ω2 ( x ) = 0 } , when some conditions on the curvature of thelevel sets are satisfied. Let { r , t , . . . , t N − } be an orthonormal basis at x ∈ Ω, where r isthe unit normal to the level set of ψ Ω2 at x and { t , . . . , t N − } the unit tangents. We canexpress the gradient operator in this coordinate system as ∇ = r ∂∂r + (cid:80) N − i =1 t i ∂∂t i . We nowhave ∇ ψ Ω2 = ∂ψ Ω2 ∂r r (since the derivative of ψ Ω2 vanishes along the tangent space of its levelset). Moreover, the eigenvalue equation ∆ ψ Ω2 + µ Ω2 ψ Ω2 = 0 expressed in this coordinate systemis given by ∂ ψ Ω2 ∂r + ( N − H ∂ψ∂r + µ Ω2 ψ Ω2 = 0, where H ( x ) is the mean curvature at x ∈ Ωof the level set of ψ . Following some computation, we get that the gradient of f is givenby ∇ f = 4 µ Ω2 ψ Ω2 ∂ψ Ω2 ∂r r + ( N − H (cid:12)(cid:12)(cid:12) ∂ψ Ω2 ∂r (cid:12)(cid:12)(cid:12) r . Clearly, for any point x ∗ in the nodal set of ψ Ω2 satisfying H ( x ∗ ) = 0 (which is certainly the case when the level set of ψ Ω2 is locally flat at x ∗ ),we have ∇ f = 0, which implies that they are critical points of f . Moreover, in computingthe entries of the Hessian of f in this coordinate frame, we first have: ∂ f∂r = 4 µ Ω2 (cid:32)(cid:12)(cid:12)(cid:12)(cid:12) ∂ψ Ω2 ∂r (cid:12)(cid:12)(cid:12)(cid:12) − µ Ω2 | ψ Ω2 | (cid:33) − N − Hµ Ω2 ψ Ω2 ∂ψ Ω2 ∂r + (cid:18) ( N − ∂H∂r − N − H (cid:19) (cid:12)(cid:12)(cid:12)(cid:12) ∂ψ Ω2 ∂r (cid:12)(cid:12)(cid:12)(cid:12) . Furthermore, when the family of level sets of ψ Ω2 is locally flat at x ∗ , we have in particularthat ∂H∂r ( x ∗ ) = 0, and it follows that ∂ f∂r ( x ∗ ) ≥
0. Also, under local flatness of the family oflevel sets, the off-diagonal entries ∂ f∂r∂t i and ∂ f∂t i ∂t j vanish for all i ∈ { , . . . , N − } , and so dothe rest of the diagonal entries of the Hessian, i.e. ∂ f∂t i ( x ∗ ) = 0 for i ∈ { , . . . , N − } . Thisyields positive semidefiniteness of the Hessian at x ∗ . Therefore, points in the nodal set of ψ Ω2 where the family of level sets of ψ Ω2 is locally flat, are local minima of f . References [1] A. Jamakovic and P. V. Mieghem, “On the robustness of complex networks by using thealgebraic connectivity,” in
International conference on research in networking . Springer,2008, pp. 183–194.2] M. Fiedler, “Algebraic connectivity of graphs,”
Czechoslovak Mathematical Journal ,vol. 23, no. 98, pp. 298–305, 1973.[3] M. Ventresca and D. Aleman, “Efficiently identifying critical nodes in large complexnetworks,”
Computational Social Networks , vol. 2:6, no. 1, 2015. [Online]. Available:https://doi.org/10.1186/s40649-015-0010-y[4] A. Arulselvan, C. Commander, L. Elefteriadou, and P. Pardalos, “Detecting criticalnodes in sparse graphs,”
Computers and Operations Research , vol. 36, no. 7, pp. 2193–2200, 2009.[5] W. Abbas, A. Laszka, Y. Vorobeychik, and X. Koutsoukos, “Improving network connec-tivity using trusted nodes and edges,” in
American Control Conference , Seattle, USA,2017, pp. 328–333.[6] M. Sheng, J. Li, and Y. Shi, “Critical nodes detection in mobile ad hoc network,” in
Int.Conf. on Advanced Information Networking and Applications , vol. 2. IEEE, 2006, pp.336–340.[7] K. P. Eswaran and R. E. Tarjan, “Augmentation problems,”
SIAM Journal on Comput-ing , vol. 5, no. 4, pp. 653–665, 1976.[8] M. Dziubi´nski and S. Goyal, “Network design and defence,”
Games and Economic Be-havior , vol. 79, pp. 30–43, 2013.[9] L. Lov´asz,
Large networks and graph limits , ser. Coloquium Publications. AmericanMathematical Society, 2012, vol. 60.[10] M. Belkin and P. Niyogi, “Towards a theoretical foundation for Laplacian-based manifoldmethods,” in
Int. Conf. on Computational Learning Theory . Springer, 2005, pp. 486–500.[11] ——, “Convergence of Laplacian eigenmaps,”
Advances in Neural Information Process-ing Systems , vol. 19, pp. 129–136, 2007.[12] A. Henrot,
Extremum problems for eigenvalues of elliptic operators . Springer Science& Business Media, 2006.[13] T. Kolokolnikov, M. S. Titcombe, and M. J. Ward, “Optimizing the fundamental Neu-mann eigenvalue for the Laplacian in a domain with small traps,”
European Journal ofApplied Mathematics , vol. 16, no. 2, pp. 161– 200, 2005.[14] V. Krishnan and S. Mart´ınez, “Identification of critical nodes for consensus in large-scalespatial networks,” in
IFAC World Congress , Toulouse, France, July 2017, pp. 14 721–14 726.[15] P. Goyal and E. Ferrara, “Graph embedding techniques, applications, and performance:A survey,” arXiv preprint arXiv:1705.02801 , 2017.[16] B. Shaw and T. Jebara, “Structure preserving embedding,” in
Proceedings of the 26thAnnual International Conference on Machine Learning . ACM, 2009, pp. 937–944.17] L. C. Evans,
Partial differential equations , ser. Graduate studies in mathematics. Prov-idence, RI: American Mathematical Society, 1998.[18] R. Murray, B. Swenson, and S. Kar, “Revisiting normalized gradient descent: Evasionof saddle points,” arXiv preprint arXiv:1711.05224 , 2017.[19] R. Atar and K. Burdzy, “On nodal lines of Neumann eigenfunctions,”
Electronic Com-munications in Probability , vol. 7, pp. 129–139, 2002.[20] D. Henry,
Geometric theory of semilinear parabolic equations . Springer, 1981.[21] M. L. de Cristoforis, “Simple Neumann eigenvalues for the Laplace operator in a domainwith a small hole,”