\section{Introduction}
The overarching goal of this paper is to find methods that help to make networks more robust against virus attacks (SIS type epidemics).  This is done by selecting nodes from the network to immunise or to remove such that the largest eigenvalue is reduced and thereby the epidemic threshold is improved \cite{chakrabarti2008epidemic}. In this paper, we derive a quadratic programming version of the problem of reducing the largest eigenvalue of a complex network (represented by the Netshield proxy-function \cite{chen2016node}) and compare the new exact method with heuristic methods proposed previously in \cite{maulana2017immunization}. This strategy is proposed for solving the  bi-objective cost-benefit optimization problem, and we discuss for the first time the exact Pareto front obtained by our new method for different example networks.

The problem that we discuss may for instance arise when combating or preventing the spread of viruses such as recently Ebola or SARS variants \cite{plaat}. To evaluate solutions an appropriate  vulnerability measure is necessary. For this paper the eigenvalue drop (abbreviated with 'eigen-drop') is used which is inversely proportional to the epidemic threshold. More precisely, the \marios{eigenvalue drop}  measures the decrease of \marios{the} maximum eigenvalue of the adjacency matrix
of a network
and the inverse relationship to the increase of the epidemic threshold \marios{,} under the SIS epidemic model. \cite{chakrabarti2008epidemic} \cite{li2013epidemic}.
The epidemic threshold is a critical value inherent to the given networks that determines the contagiousness a virus requires to infect the entire network (exponential growth) or to disappear (exponential decay). The SIS model, a is a special case of the so-called contact process where a virus can spread from an infected node to any of its susceptible nod neighbours and infect those neighbours.
At some later point in time the infected nodes recover and become susceptible to infection again.

A network or graph $G$ consists a pair $(V,E)$. Here $V$ is a set of nodes $V=(v_1, \dots, v_n)$. $E \subseteq V \times V$ is a set of edges representing connections between the nodes. A graph can also be represented as an adjacency matrix $A(V,E) \in \{0,1\}^{n \times n}$ with $a_{ij} = 1$ if $(v_i, v_j) \in E$ or $a_{ij} = 0$ if $(v_i, v_j) \notin E$. The first or maximum eigenvalue of this graph will be denoted with $\lambda$ and the corresponding eigenvector with $u$.

Note, that there that in the literature various strategies for reducing the vulnerability of a complex network to virus attacks have been suggested. Examples are, for instance,  random immunisation, acquaintance immunisation, and target immunisation.
Many of these strategies can be seen as heuristics and do not specifically relate to the spread dynamics of viruses in real world networks.
In Chakrabarti et al. \cite{chakrabarti2008epidemic} it is argued that it is \marios{of} paramount importance to take into account the dynamics of the contact process and to chose an appropriate measure with respect to the global structure of the complex networks.
Local measures, such as the degree of a vertex, do in general lead to mediocre performance and can be even misleading. For example targeted immunisation would choose the nodes with the highest degrees (hubs).  To focus the strategy, at least partially, on lower degree nodes, may seem counter-intuitive.  However, it is not always the case that the immunisation of the highest degree nodes will reduce the vulnerability of the network.
In contrast, it has been shown that focusing on the reduction of the largest eigenvalue of the adjacency matrix of the graph is a effective way to reduce the epidemic threshold which determines whether the number of infections grows exponentially or decreases exponentially.
In particular in the early stages of an outbreak there is a broad consensus of the effectiveness of this strategy.
\marios{In} \cite{chakrabarti2008epidemic} \marios{the authors} provide the example of the barbell graph \marios{(see Fig~\ref{fig:example})} that demonstrates this, where node 13 is of crucial importance although it has low degree.

\begin{figure}
\scalebox{0.8}{
\begin{tikzpicture}[<->,>=stealth,shorten >=1pt,auto,node distance=2cm,
                    thick,main node/.style={circle,draw,font=\sffamily\Large\bfseries}]
  \node[main node] (1) {$\ 1$};
  \node[main node] (2) [left of=1] {$\ 2$};
  \node[main node] (3) [below left of=2] {$\ 3$};
  \node[main node] (4) [below right of=3] {$\ 4$};
  \node[main node] (5) [right of=4] {$\ 5$};
  \node[main node] (6) [below right of=1] {$\ 6$};
  \node[main node] (13) [right of=6] {$13$};
  \node[main node] (7) [right of=13]{$\ 7$};
  \node[main node] (8) [above right of=7] {$\ 8$};
  \node[main node] (9) [right of=8] {$\ 9$};
  \node[main node] (10) [below right of=9] {$10$};
  \node[main node] (11) [below left of=10] {$11$};
  \node[main node] (12) [left of=11] {$12$};
\path[every node/.style={font=\sffamily\small}]
  (1) edge [right] (2)
  (1) edge [right] (3)
  (1) edge [right] (4)
  (1) edge [right] (5)
  (1) edge [right] (6)
  (2) edge [right] (3)
  (2) edge [right] (4)
  (2) edge [right] (5)
  (2) edge [right] (6)
  (3) edge [right] (4)
  (3) edge [right] (5)

  (4) edge [right] (5)
  (4) edge [right] (6)
  (5) edge [right] (6)

  (7) edge [right] (8)
  (7) edge [right] (9)

  (7) edge [right] (11)
  (7) edge [right] (12)
  (8) edge [right] (9)
  (8) edge [right] (10)
  (8) edge [right] (11)
  (8) edge [right] (12)
  (9) edge [right] (10)
  (9) edge [right] (11)
  (9) edge [right] (12)
  (10) edge [right] (11)
  (10) edge [right] (12)
  (11) edge [right] (12)
  (6) edge [right] (13)
  (13) edge [right] (7)
  ;
\end{tikzpicture}
}
\caption{\label{fig:example} A small version of the “bar-bell” graph. Two cliques of the same size connected with a bridge. Using this graph, it can be demonstrated that targeted immunization, that is focusing on the removal of high degree nodes, is not always the best strategy in immunization. Adapted from \cite{chakrabarti2008epidemic}.}
\end{figure}

\begin{definition}
Given a network $G$
and a network $G'$ where $G'$ is a sub-graph of $G$ with some its nodes and adjacent edges removed, $\Delta\lambda$ or eigen-drop
is defined as the difference between the maximum eigenvalue of the adjacency matrix of $G$ and the maximum eigenvalue of the adjacency matrix of $G'$.
\end{definition}
\begin{definition}
The $k$-node immunisation \marios{problem}: \marios{G}iven a graph $G=(V,E)$ and $k \marios{\in \mathbb{N}}$, \marios{the $k$-node immunisation problem aims in finding} a set of nodes $S \subseteq V$ with $|S| = k$, such that the removal of these nodes from $G$ maximises the eigen-drop $\Delta\lambda$.
\end{definition}

 \begin{figure}[t]
     \centering
     \includegraphics[width=\textwidth]{eigendropvsdegree}
     \caption{Table with eigen-drops vs. degrees of extended barbell graph}
     \label{fig:eigendroptable}
 \end{figure}
It has been shown in \cite{chen2016node} that this problem is NP-hard. Therefore heuristic methods have been suggested for solving this problem such as the NetShield and NetShield+ algorithms in \cite{chen2016node} and a problem specific genetic algorithm in \cite{maulana2017immunization}.
The drawback of these heuristics is that they are designed for the $k$-node immunisation problem which requires that a good value of $k$ is known in advance.
In addition, the heuristics treat\marios{this problem} as if every node requires the same effort to remove. Therefore in this paper, we reformulate the $k$-node immunisation problem to a multi-objective one.
This is similar to [7], where multi-objective optimisation using an evolutionary optimisation heuristic was used to take into account the cost of removal.
Indeed, the nodes with maximal eigen-drop do not coincide with the nodes with a high degree, as is seen in Figure \ref{fig:eigendroptable}. Node 13 in the earlier discussed barbell graph has only a degree of $2$ but it significantly reduces the Eigenvalue.
\begin{definition}
The multi-objective immunisation problem given a graph $G=(V,E)$, a cost denoted with $\mathrm{Cost}(v)$
for each $v \in V$ and $S \subseteq V$ reads
\begin{eqnarray}
    f_1(S) = \Delta\lambda &\rightarrow& \max\\
    f_2(S) = \sum_{v \in S} \mathrm{Cost}(v) &\rightarrow& \min
\end{eqnarray}
\end{definition}
This formulation requires no value for $k$ to be known a-priori and takes the cost of removal in account. As this is a multi-objective problem we are now interested in finding the efficient set and its corresponding Pareto front. To approximate this Pareto front we use and evaluate four different methods. The first two extend the NetShield and NetShield+ algorithms by substituting the first objective with the heuristic used by these methods and then applying the $\epsilon$-constraint method. The second two are two different genetic algorithms specifically designed for multi-objective optimisation problems.

\section{Multiobjective Quadratic Programming Formulation}

The first method for approximating the Pareto front of the multi-objective immunisation problem is based on the NetShield and NetShield+ algorithms designed by Chen Chen et al. \cite{chen2016node}. These methods are designed for the $k$-node immunisation problem and are briefly discussed here. Core to the design of these algorithms is a function called Shield-value.

\begin{definition}
Given the adjacency matrix $A$ of a graph $G(V,E)$, its first eigenvalue $\lambda$, the corresponding eigenvector $u$ and an input set of nodes $S \subseteq V$, the Shield-value function is defined as:
\begin{equation}
    Sv(S) = \sum_{i \in S} 2 \lambda u_{i}^{2} - \sum_{i,j \in S} 2 u_{i}u_{j}A_{ij}
\end{equation}
\end{definition}

The Shield-value function gives an approximation of $\Delta\lambda$ if all nodes in the set $S$ were to be removed from the graph. The NetShield algorithm then finds a set of $k$ nodes that approximates the maximization of this function via greedy selection.

As the cardinality of $S$ grows, the Shield-value function becomes less accurate. The NetShield+ algorithm therefore introduces an extra parameter called the batch size. Instead of finding $k$ nodes at once, a set of $b$ nodes is found and added to the solution. Then these nodes are removed from the network. The new network is used to compute a new Shield-value function, that is again maximized for a set of $b$ nodes. These are then again added to the solution set and this process continues until $k$ nodes have been removed from the network.

To extend these algorithms to work with the multi-objective immunisation problem, we substitute the eigen-drop objective with the Shield-value function. Furthermore, we define the problem as a quadratic multi-objective program by representing the solution with a binary vector $x$. If the node $i$ is in the solution set, $x_i$ will be 1. Otherwise $x_i$ will be 0:

\begin{definition}
Given the adjacency matrix $A$ of a graph $G(V,E)$, its first eigenvalue $\lambda$, and the corresponding eigenvector $u$, the Shield-value function with cost objective is:
\begin{equation}
    f_{1}(x) = \sum_{i=1}^{m} 2 \lambda  u_{i}^{2} x_{i} - \sum_{i=1}^{m} \sum_{j=i+1}^{m} 2 u_{i} u_{j} A_{ij} x_{i} x_{j} \to \max \\
\end{equation}
\begin{equation}
    f_{2}(x) = \sum_{i=1}^{m} x_{i} \mathrm{Cost(i)} \to \min\\
\end{equation}
Subject to:
\begin{equation}
x \in \{0,1\}^{m}
\end{equation}
\end{definition}
To approximate the Pareto front the $\epsilon$-constraint method can be applied \cite{KaisaNMO}. For this method, one of the objectives is transformed into a constraint smaller or equal than $\epsilon$. In this case it will be the second cost objective.
\begin{definition}
Given the adjacency matrix $A$ of a graph $G(V,E)$, its first eigenvalue $\lambda$, the corresponding eigenvector $u$, and some value of $\epsilon$, the Shield-value function with cost constraint is:
\begin{equation}
    f_{1}(x) = \sum_{i=1}^{m} 2 \lambda  u_{i}^{2} x_{i} - \sum_{i=1}^{m} \sum_{j=i+1}^{m} 2 u_{i} u_{j} A_{ij} x_{i} x_{j} \to \max \\
\end{equation}
Subject to:
\begin{equation}
    f_{2}(x) = \sum_{i=1}^{m} x_{i}  \mathrm{Cost}(i) \leq \epsilon \\
\end{equation}
\begin{equation}
    x \in \{0,1\}^{m}
\end{equation}
\end{definition}

By choosing a concrete value of $\epsilon$, the problem is transformed into a quadratic program with a linear constraint. Problems such as these can be solved with a quadratic problem solver via branch-and-bound based methods. By solving multiple programs with different values of $\epsilon$, the Pareto front of the multi-objective Shield-value problem can be found. As the Shield-value is an approximation of $\Delta\lambda$, this Pareto front should therefore also be an approximation of the original problem.

This method can be extended analogously to how the original NetShield algorithm can be extended to NetShield+. Instead of finding a set of nodes that maximises the Shield-value objective at once, an extra batch size parameter $b$ can be introduced. Then a solution that maximises the Shield-value function with only $b$ nodes can be found. These nodes are then added to the complete solution and removed from the network. A new quadratic program can be created with a new Shield-value function computed from the new network. This process continues, adding $b$ nodes to the solution set at every step. This process stops when no more nodes can be added. This occurs when either all nodes have already been added to the solution set, or if any of the nodes not yet added would violate the cost constraint.

\section{Genetic Algorithms}
The advantage of using genetic algorithms over the NetShield based methods described in the previous section, is that they can be made to work directly on the eigen-drop. This sidesteps the need of using a possible inaccurate approximation of the eigen-drop. In addition, by sampling the search space in an efficient manner, genetic algorithms can also consider more candidate solutions that the NetShield methods will. Therefore, it is possible that better Pareto front approximations can be found by these meta-heuristics.
The GAs used in this paper are specifically designed for multi-objective problems. They use specialised selection operators that aim for both convergence to the Pareto front and spread over the Pareto front. The GAs used are NSGA-II \cite{deb2002fast} and SMS-EMOA \cite{emmerich2005emo}.
In addition to this, it is also possible to hybridise the GAs with the NetShield methods. This is done by initialising the GAs with the solutions found by the NetShield methods. This can cut out the potentially large search effort by the GAs to converge on the Pareto front by starting them from what already is a good approximation. Then the GAs may further refine the solutions using their advantages over the NetShield methods.

\section{Experiments and Results}

A specific cost function is required to define $f_2$.
This cost function should be a good local measure for the effort required for the removal of a node. The cost function we used is the degree of each node, as a highly connected node is likely to be more difficult to remove from the network than a node with less incoming and outgoing edges.

All of the GAs were run 5 times under each configuration, both when the population was initialised at random and when it was initialised with the NetShield solutions. The populations were set to a size of $100$ for the random initialisation. The mutation probability $p_m$ was set to $1/n$, with $n$ being the number of vertices of the graph. Crossover probability $p_c$ was set to $0.75$. All GAs were run with $10000$ iterations of the main loop.
All results of the GAs are plotted as the first attainment curve \cite{Fonseca96a}. All points on these curves are weakly dominated by only 1 run out of the 5 and are therefore a best case scenario of the GAs

Both the NetShield and NetShield+ methods with the $\epsilon$-constraint method were tested. For the NetShield+ method, the batch size was set to 1. The resolution of the Pareto front approximation depends on how many different values of $\epsilon$ are sampled. As the cost function chosen uses only non-negative integers, it is possible to get the best possible resolution by sampling only a finite amount of points: from 0 to the sum of all degrees increasing $\epsilon$ by 1 every step. The quadratic program solver used is Gurobi\cite{bixby2011gurobi}.

All results shown are for the following set of four graphs:

\begin{enumerate}
    \item  \textbf{Pandemic}: Based on the Pandemic board game in which a global virus outbreak is fought. The graph connects 27 cities in the world to each other with 93 edges. \cite{maulana2017immunization}. See Figure \ref{fig:conf1pand} (right).
     \item \textbf{Conference Day 1}: Interactions between members of a conference on the first day. Only the largest connected components has been selected from this graph. The graph consists of 190 nodes and 703 edges. See Figure \ref{fig:conf1pand} (left). Taken from \url{www.sociopatterns.org/datasets/infectioussociopatterns}.
     \item \textbf{Erd\H{o}s-R\'enyi graph}: Graph sampled from the Erd\H{o}s-R\'enyi random graph model. This graph has 100 nodes and 294 edges.
     \item \textbf{Barab\'asi-Albert graph}: Graph sampled from the Barab\'asi-Albert graph model.  This graph also has 100 nodes and 294 edges. See Fig. \ref{fig:bara_j1}
\end{enumerate}
\begin{figure}[t]
\includegraphics[width=\textwidth]{conf1-and-pandemic}
\caption{\label{fig:conf1pand}Conference Day 1 and Pandemic Network.}
\end{figure}
\subsection{NetShield with $\epsilon$-constraints and GAs}
The results of both the NetShield methods and the randomly initialised GAs are plotted in Figure \ref{fig:res_panconf} for the Pandemic and Conference day graph and in Figure \ref{fig:res_erbara} for the Erd\H{o}s-R\'enyi and Barab\'asi-Albert graph. When comparing the NSGA-II algorithm to the SMS-EMOA algorithm, no clear differences show. It changes from graph to graph which algorithm finds the better Pareto front approximation and they consistently lie very closely together.

Difference do show when comparing the NetShield with the NetShield+ method. Sometimes the difference are large, such as for the Conference day 1 graph and Pandemic graph. For the Barab\'asi-Albert and Erd\H{o}s-R\'enyi graphs the differences are smaller, but the NetShield+ method still tends to give the better results. This is likely due to the Shield-Value losing accuracy when the number of nodes removed increases. Initially the performance of NetShield is very similar to NetShield+. When the allowed cost increases and consequently more nodes can be selected, the NetShield+ method can find solutions that are significantly better.

At the rightmost extremes however, the NetShield method sometimes finds some solutions that dominate those found by the NetShield+ method. See, for example, the Erd\H{o}s-R\'enyi graph and the Barab\'asi-Albert graph. This may be because the NetShield+ method with a batch size of 1 is more greedy than the NetShield method. With a batch size of 1, the NetShield+ method selects at every step the node with the highest eigenscore that would not violate the $\epsilon$-constraint when added to the solution. It then recomputes a new Shield-value function with the node removed. The NetShield method however, only computes the Shield-value function at the beginning and selects multiple nodes at once to optimise this function. In this way it can take a more global view of the problem. If the Shield-Value then happens to still be a good approximation of the
eigen-drop, better solution may be found. A possible approach is to repeat the NetShield+ method several times with different batch sizes if time allows. Then all results can be combined for the most accurate Pareto front approximation.

When comparing both the NetShield methods with the GAs, the GAs give very competitive performance when the networks have relatively few nodes. This means the search space is smaller and the GAs have enough time to converge on the Pareto front. This results in some parts of the Pareto front being approximated better by the GAs, because they can work directly on the eigen-drop. This is most notably the case for the Pandemic graph. Here both NetShield method and the NetShield+ method to a lesser degree have difficulty approximating the Pareto front. This is likely caused by this graph having a low maximum degree. This means that the Shield-value approximation loses accuracy quickly \cite{chen2016node}.

The results also show that the Pareto fronts do not form a single distinctive shape. The Pareto front for the Erd\H{o}s-R\'enyi graph is mostly linear. At the rightmost extreme however, two solutions are found where large gains in eigen-drop can be made with a comparatively small cost increase. This results in the Pareto front having a concave section at the end. This is the opposite for the Pareto front for the Conference day 1 graph. For this graph there is a clear case of diminishing returns: it costs increasingly more to get the same improvement in terms of eigen-drop the further the eigen-drop increases.

The Barab\'asi-Albert graph Pareto front consists of several sections that are mostly linear, but with gaps in between this sections. At these gaps, large increases in eigen-drop are suddenly gained for low costs. This is the result of the preferential attachment model used to generate this graph. At those points the value of $\epsilon$ allows replacing a larger selection of smaller cost nodes with one of the highly connected hub nodes. While these nodes have high cost, their impact on eigen-drop is still disproportionate to their cost. Two solutions with this graph are visualised in figure \ref{fig:bara_j1}: one at the left side of a gap and one at the right side.

\begin{figure}
  \centering
    \includegraphics[width=\textwidth]{Images/panconf_attaintment_netshield.png}
  \caption{Results GAs and NetShield(+) with $\epsilon$-constraint method}
  \label{fig:res_panconf}
\end{figure}

\begin{figure}
  \centering
    \includegraphics[width=\textwidth]{Images/Erdos_Bara_V100_attaintment_netshield.png}
  \caption{Results GAs and NetShield(+) with $\epsilon$-constraint method}
  \label{fig:res_erbara}
\end{figure}

\begin{figure}[t]
  \centering
    \includegraphics[width=\textwidth]{Images/j1}
    \caption{Selected nodes are red and denoted with $\times$. Nodes have been scaled with degree. Left: 6 selected nodes, cost of 36, $\Delta\lambda$ of 0.975.
            Right: 1 selected node, cost of 37, $\Delta\lambda$ of 1.455}
  \label{fig:bara_j1}
\end{figure}

\subsection{Hybrid GA approach}
The results with the GAs initialised from the results from the NetShield methods for the Pandemic and Barab\'asi-Albert graph are shown in Figure \ref{fig:res_atins}. They are shown together with the initialisation sets. The most notable improvements found by the GAs are for the Pandemic graph. Here the results of the inaccuracies of the Shield-value have been corrected. It appears that in these cases, the GAs have the ability to repair such issues.
The improvements for the Barab\'asi-Albert graph are more minor. Either the initialisation sets are already close to the Pareto fronts or there may not be enough diversity in the initial populations for the GAs to find better solutions.
\begin{figure}
  \centering
    \includegraphics[width=\textwidth]{Images/pandbara_attainment_nsinit}
  \caption{Results hybrid GAs}
  \label{fig:res_atins}
\end{figure}
\section{Conclusion}
In this paper, it is shown that the NetShield and NetShield+ algorithms can be extended with an $\epsilon$-constraint method, using an exact quadratic programming solver (here: Gurobi). In this manner, a multi-objective variant of the node immunisation problem can be solved with a cost function added which is proportional to the effort of the node removal. The performance is mostly equivalent or better than two multi-objective genetic algorithms specifically designed for multi-objective optimisation, except for cases where the Shield-value function is inaccurate due to characteristics of the network. In general the NetShield+ method is more robust than the NetShield method, but there are exceptions. Combining the GAs with the NetShield algorithm as initial population only provided small further improvements. Therefore, if time permits and the most accurate Pareto front approximation is required, a valid approach would be to use all methods and combine the end results.
The results also show that there does not appear to be a typical shape to the Pareto fronts resulting from this problem. They are dependent on both the topology of the network and on the cost function.

A main contribution of this paper is the insight that, if time permits, exact and problem-specific methods approximation should be used, which are often far better than Pareto front approximations obtained by general meta-heuristics.
Based on these insight, it will be more effective to develop strategies for controlling real-world networks when the goal is to prevent epidemic outbreaks.
It should be noted, however, that we focused solely on the eigen-value drop. While this is an effective measure under the SIS infection model, the properties of real world epidemics may not be fully captured by this model. Such additional aspects are mentioned, for instance, in \cite{plaat}. In addition, we also assume that any nodes in the network can be immunised.  This also may not translate well to real world scenarios. Therefore, future work should also take a broader view of the problem than further improving the eigen-drop via the process of node removal.

A note of precaution shall be provided here, when applying the model to epidemics in a late stage. A crucial assumption is that the epidemic is in an earlier stage, which will imply that it is likely many neighboring nodes of an infected node are not yet infected. In a late stage of an epidemic this can no-longer be assumed. Moreover, the dynamics are \marios{of the} susceptible-infected-suceptible (SIS) type, and not of the susceptible-infected-recovered (SIR) type, which would be also a common real-world model, e.g., because it would be applicable to the recent SARS-CoV2 pandemic.
It is conjectured, that also in the SIR dynamics the largest eigenvalue is an important value to focus on when dampening or preventing a virus outbreak. When it comes to more realistic and fine grained models of network dynamics, it is probably not straightforward anymore to use a single indicator, such as the eigen-drop, and more complex simulation models are required such as the event-based simulation model in \cite{michalak2020evolutionary}. For an excellent overview various  epidemic models and dynamics on complex networks the reader is referred to \cite{pastor2015epidemic}.

Whereas in our study exact QP solvers were applicable due to the quadratic equations in the objective function, in the more general setting this will no longer be the case and black-box optimization, such as SMS-EMOA will be very useful.
\paragraph{Software*}: All source code (Python) of algorithm implementations and network data of this study is made free available under:\\
\url{https://github.com/joostnibbeling/node-immunisation}

\addcontentsline{toc}{section}{References}

\bibliographystyle{abbrv}

\section{Introduction}
The overarching goal of this paper is to find methods that help to make networks more robust against virus attacks (SIS type epidemics).  This is done by selecting nodes from the network to immunise or to remove such that the largest eigenvalue is reduced and thereby the epidemic threshold is improved \cite{chakrabarti2008epidemic}. In this paper, we derive a quadratic programming version of the problem of reducing the largest eigenvalue of a complex network (represented by the Netshield proxy-function \cite{chen2016node}) and compare the new exact method with heuristic methods proposed previously in \cite{maulana2017immunization}. This strategy is proposed for solving the  bi-objective cost-benefit optimization problem, and we discuss for the first time the exact Pareto front obtained by our new method for different example networks.

The problem that we discuss may for instance arise when combating or preventing the spread of viruses such as recently Ebola or SARS variants \cite{plaat}. To evaluate solutions an appropriate  vulnerability measure is necessary. For this paper the eigenvalue drop (abbreviated with 'eigen-drop') is used which is inversely proportional to the epidemic threshold. More precisely, the \marios{eigenvalue drop}  measures the decrease of \marios{the} maximum eigenvalue of the adjacency matrix
of a network
and the inverse relationship to the increase of the epidemic threshold \marios{,} under the SIS epidemic model. \cite{chakrabarti2008epidemic} \cite{li2013epidemic}.
The epidemic threshold is a critical value inherent to the given networks that determines the contagiousness a virus requires to infect the entire network (exponential growth) or to disappear (exponential decay). The SIS model, a is a special case of the so-called contact process where a virus can spread from an infected node to any of its susceptible nod neighbours and infect those neighbours.
At some later point in time the infected nodes recover and become susceptible to infection again.

A network or graph $G$ consists a pair $(V,E)$. Here $V$ is a set of nodes $V=(v_1, \dots, v_n)$. $E \subseteq V \times V$ is a set of edges representing connections between the nodes. A graph can also be represented as an adjacency matrix $A(V,E) \in \{0,1\}^{n \times n}$ with $a_{ij} = 1$ if $(v_i, v_j) \in E$ or $a_{ij} = 0$ if $(v_i, v_j) \notin E$. The first or maximum eigenvalue of this graph will be denoted with $\lambda$ and the corresponding eigenvector with $u$.

Note, that there that in the literature various strategies for reducing the vulnerability of a complex network to virus attacks have been suggested. Examples are, for instance,  random immunisation, acquaintance immunisation, and target immunisation.
Many of these strategies can be seen as heuristics and do not specifically relate to the spread dynamics of viruses in real world networks.
In Chakrabarti et al. \cite{chakrabarti2008epidemic} it is argued that it is \marios{of} paramount importance to take into account the dynamics of the contact process and to chose an appropriate measure with respect to the global structure of the complex networks.
Local measures, such as the degree of a vertex, do in general lead to mediocre performance and can be even misleading. For example targeted immunisation would choose the nodes with the highest degrees (hubs).  To focus the strategy, at least partially, on lower degree nodes, may seem counter-intuitive.  However, it is not always the case that the immunisation of the highest degree nodes will reduce the vulnerability of the network.
In contrast, it has been shown that focusing on the reduction of the largest eigenvalue of the adjacency matrix of the graph is a effective way to reduce the epidemic threshold which determines whether the number of infections grows exponentially or decreases exponentially.
In particular in the early stages of an outbreak there is a broad consensus of the effectiveness of this strategy.
\marios{In} \cite{chakrabarti2008epidemic} \marios{the authors} provide the example of the barbell graph \marios{(see Fig~\ref{fig:example})} that demonstrates this, where node 13 is of crucial importance although it has low degree.

\begin{figure}
\scalebox{0.8}{
\begin{tikzpicture}[<->,>=stealth,shorten >=1pt,auto,node distance=2cm,
                    thick,main node/.style={circle,draw,font=\sffamily\Large\bfseries}]
  \node[main node] (1) {$\ 1$};
  \node[main node] (2) [left of=1] {$\ 2$};
  \node[main node] (3) [below left of=2] {$\ 3$};
  \node[main node] (4) [below right of=3] {$\ 4$};
  \node[main node] (5) [right of=4] {$\ 5$};
  \node[main node] (6) [below right of=1] {$\ 6$};
  \node[main node] (13) [right of=6] {$13$};
  \node[main node] (7) [right of=13]{$\ 7$};
  \node[main node] (8) [above right of=7] {$\ 8$};
  \node[main node] (9) [right of=8] {$\ 9$};
  \node[main node] (10) [below right of=9] {$10$};
  \node[main node] (11) [below left of=10] {$11$};
  \node[main node] (12) [left of=11] {$12$};
\path[every node/.style={font=\sffamily\small}]
  (1) edge [right] (2)
  (1) edge [right] (3)
  (1) edge [right] (4)
  (1) edge [right] (5)
  (1) edge [right] (6)
  (2) edge [right] (3)
  (2) edge [right] (4)
  (2) edge [right] (5)
  (2) edge [right] (6)
  (3) edge [right] (4)
  (3) edge [right] (5)

  (4) edge [right] (5)
  (4) edge [right] (6)
  (5) edge [right] (6)

  (7) edge [right] (8)
  (7) edge [right] (9)

  (7) edge [right] (11)
  (7) edge [right] (12)
  (8) edge [right] (9)
  (8) edge [right] (10)
  (8) edge [right] (11)
  (8) edge [right] (12)
  (9) edge [right] (10)
  (9) edge [right] (11)
  (9) edge [right] (12)
  (10) edge [right] (11)
  (10) edge [right] (12)
  (11) edge [right] (12)
  (6) edge [right] (13)
  (13) edge [right] (7)
  ;
\end{tikzpicture}
}
\caption{\label{fig:example} A small version of the “bar-bell” graph. Two cliques of the same size connected with a bridge. Using this graph, it can be demonstrated that targeted immunization, that is focusing on the removal of high degree nodes, is not always the best strategy in immunization. Adapted from \cite{chakrabarti2008epidemic}.}
\end{figure}

\begin{definition}
Given a network $G$
and a network $G'$ where $G'$ is a sub-graph of $G$ with some its nodes and adjacent edges removed, $\Delta\lambda$ or eigen-drop
is defined as the difference between the maximum eigenvalue of the adjacency matrix of $G$ and the maximum eigenvalue of the adjacency matrix of $G'$.
\end{definition}
\begin{definition}
The $k$-node immunisation \marios{problem}: \marios{G}iven a graph $G=(V,E)$ and $k \marios{\in \mathbb{N}}$, \marios{the $k$-node immunisation problem aims in finding} a set of nodes $S \subseteq V$ with $|S| = k$, such that the removal of these nodes from $G$ maximises the eigen-drop $\Delta\lambda$.
\end{definition}

 \begin{figure}[t]
     \centering
     \includegraphics[width=\textwidth]{eigendropvsdegree}
     \caption{Table with eigen-drops vs. degrees of extended barbell graph}
     \label{fig:eigendroptable}
 \end{figure}
It has been shown in \cite{chen2016node} that this problem is NP-hard. Therefore heuristic methods have been suggested for solving this problem such as the NetShield and NetShield+ algorithms in \cite{chen2016node} and a problem specific genetic algorithm in \cite{maulana2017immunization}.
The drawback of these heuristics is that they are designed for the $k$-node immunisation problem which requires that a good value of $k$ is known in advance.
In addition, the heuristics treat\marios{this problem} as if every node requires the same effort to remove. Therefore in this paper, we reformulate the $k$-node immunisation problem to a multi-objective one.
This is similar to [7], where multi-objective optimisation using an evolutionary optimisation heuristic was used to take into account the cost of removal.
Indeed, the nodes with maximal eigen-drop do not coincide with the nodes with a high degree, as is seen in Figure \ref{fig:eigendroptable}. Node 13 in the earlier discussed barbell graph has only a degree of $2$ but it significantly reduces the Eigenvalue.
\begin{definition}
The multi-objective immunisation problem given a graph $G=(V,E)$, a cost denoted with $\mathrm{Cost}(v)$
for each $v \in V$ and $S \subseteq V$ reads
\begin{eqnarray}
    f_1(S) = \Delta\lambda &\rightarrow& \max\\
    f_2(S) = \sum_{v \in S} \mathrm{Cost}(v) &\rightarrow& \min
\end{eqnarray}
\end{definition}
This formulation requires no value for $k$ to be known a-priori and takes the cost of removal in account. As this is a multi-objective problem we are now interested in finding the efficient set and its corresponding Pareto front. To approximate this Pareto front we use and evaluate four different methods. The first two extend the NetShield and NetShield+ algorithms by substituting the first objective with the heuristic used by these methods and then applying the $\epsilon$-constraint method. The second two are two different genetic algorithms specifically designed for multi-objective optimisation problems.

\section{Multiobjective Quadratic Programming Formulation}

The first method for approximating the Pareto front of the multi-objective immunisation problem is based on the NetShield and NetShield+ algorithms designed by Chen Chen et al. \cite{chen2016node}. These methods are designed for the $k$-node immunisation problem and are briefly discussed here. Core to the design of these algorithms is a function called Shield-value.

\begin{definition}
Given the adjacency matrix $A$ of a graph $G(V,E)$, its first eigenvalue $\lambda$, the corresponding eigenvector $u$ and an input set of nodes $S \subseteq V$, the Shield-value function is defined as:
\begin{equation}
    Sv(S) = \sum_{i \in S} 2 \lambda u_{i}^{2} - \sum_{i,j \in S} 2 u_{i}u_{j}A_{ij}
\end{equation}
\end{definition}

The Shield-value function gives an approximation of $\Delta\lambda$ if all nodes in the set $S$ were to be removed from the graph. The NetShield algorithm then finds a set of $k$ nodes that approximates the maximization of this function via greedy selection.

As the cardinality of $S$ grows, the Shield-value function becomes less accurate. The NetShield+ algorithm therefore introduces an extra parameter called the batch size. Instead of finding $k$ nodes at once, a set of $b$ nodes is found and added to the solution. Then these nodes are removed from the network. The new network is used to compute a new Shield-value function, that is again maximized for a set of $b$ nodes. These are then again added to the solution set and this process continues until $k$ nodes have been removed from the network.

To extend these algorithms to work with the multi-objective immunisation problem, we substitute the eigen-drop objective with the Shield-value function. Furthermore, we define the problem as a quadratic multi-objective program by representing the solution with a binary vector $x$. If the node $i$ is in the solution set, $x_i$ will be 1. Otherwise $x_i$ will be 0:

\begin{definition}
Given the adjacency matrix $A$ of a graph $G(V,E)$, its first eigenvalue $\lambda$, and the corresponding eigenvector $u$, the Shield-value function with cost objective is:
\begin{equation}
    f_{1}(x) = \sum_{i=1}^{m} 2 \lambda  u_{i}^{2} x_{i} - \sum_{i=1}^{m} \sum_{j=i+1}^{m} 2 u_{i} u_{j} A_{ij} x_{i} x_{j} \to \max \\
\end{equation}
\begin{equation}
    f_{2}(x) = \sum_{i=1}^{m} x_{i} \mathrm{Cost(i)} \to \min\\
\end{equation}
Subject to:
\begin{equation}
x \in \{0,1\}^{m}
\end{equation}
\end{definition}
To approximate the Pareto front the $\epsilon$-constraint method can be applied \cite{KaisaNMO}. For this method, one of the objectives is transformed into a constraint smaller or equal than $\epsilon$. In this case it will be the second cost objective.
\begin{definition}
Given the adjacency matrix $A$ of a graph $G(V,E)$, its first eigenvalue $\lambda$, the corresponding eigenvector $u$, and some value of $\epsilon$, the Shield-value function with cost constraint is:
\begin{equation}
    f_{1}(x) = \sum_{i=1}^{m} 2 \lambda  u_{i}^{2} x_{i} - \sum_{i=1}^{m} \sum_{j=i+1}^{m} 2 u_{i} u_{j} A_{ij} x_{i} x_{j} \to \max \\
\end{equation}
Subject to:
\begin{equation}
    f_{2}(x) = \sum_{i=1}^{m} x_{i}  \mathrm{Cost}(i) \leq \epsilon \\
\end{equation}
\begin{equation}
    x \in \{0,1\}^{m}
\end{equation}
\end{definition}

By choosing a concrete value of $\epsilon$, the problem is transformed into a quadratic program with a linear constraint. Problems such as these can be solved with a quadratic problem solver via branch-and-bound based methods. By solving multiple programs with different values of $\epsilon$, the Pareto front of the multi-objective Shield-value problem can be found. As the Shield-value is an approximation of $\Delta\lambda$, this Pareto front should therefore also be an approximation of the original problem.

This method can be extended analogously to how the original NetShield algorithm can be extended to NetShield+. Instead of finding a set of nodes that maximises the Shield-value objective at once, an extra batch size parameter $b$ can be introduced. Then a solution that maximises the Shield-value function with only $b$ nodes can be found. These nodes are then added to the complete solution and removed from the network. A new quadratic program can be created with a new Shield-value function computed from the new network. This process continues, adding $b$ nodes to the solution set at every step. This process stops when no more nodes can be added. This occurs when either all nodes have already been added to the solution set, or if any of the nodes not yet added would violate the cost constraint.

\section{Genetic Algorithms}
The advantage of using genetic algorithms over the NetShield based methods described in the previous section, is that they can be made to work directly on the eigen-drop. This sidesteps the need of using a possible inaccurate approximation of the eigen-drop. In addition, by sampling the search space in an efficient manner, genetic algorithms can also consider more candidate solutions that the NetShield methods will. Therefore, it is possible that better Pareto front approximations can be found by these meta-heuristics.
The GAs used in this paper are specifically designed for multi-objective problems. They use specialised selection operators that aim for both convergence to the Pareto front and spread over the Pareto front. The GAs used are NSGA-II \cite{deb2002fast} and SMS-EMOA \cite{emmerich2005emo}.
In addition to this, it is also possible to hybridise the GAs with the NetShield methods. This is done by initialising the GAs with the solutions found by the NetShield methods. This can cut out the potentially large search effort by the GAs to converge on the Pareto front by starting them from what already is a good approximation. Then the GAs may further refine the solutions using their advantages over the NetShield methods.

\section{Experiments and Results}

A specific cost function is required to define $f_2$.
This cost function should be a good local measure for the effort required for the removal of a node. The cost function we used is the degree of each node, as a highly connected node is likely to be more difficult to remove from the network than a node with less incoming and outgoing edges.

All of the GAs were run 5 times under each configuration, both when the population was initialised at random and when it was initialised with the NetShield solutions. The populations were set to a size of $100$ for the random initialisation. The mutation probability $p_m$ was set to $1/n$, with $n$ being the number of vertices of the graph. Crossover probability $p_c$ was set to $0.75$. All GAs were run with $10000$ iterations of the main loop.
All results of the GAs are plotted as the first attainment curve \cite{Fonseca96a}. All points on these curves are weakly dominated by only 1 run out of the 5 and are therefore a best case scenario of the GAs

Both the NetShield and NetShield+ methods with the $\epsilon$-constraint method were tested. For the NetShield+ method, the batch size was set to 1. The resolution of the Pareto front approximation depends on how many different values of $\epsilon$ are sampled. As the cost function chosen uses only non-negative integers, it is possible to get the best possible resolution by sampling only a finite amount of points: from 0 to the sum of all degrees increasing $\epsilon$ by 1 every step. The quadratic program solver used is Gurobi\cite{bixby2011gurobi}.

All results shown are for the following set of four graphs:

\begin{enumerate}
    \item  \textbf{Pandemic}: Based on the Pandemic board game in which a global virus outbreak is fought. The graph connects 27 cities in the world to each other with 93 edges. \cite{maulana2017immunization}. See Figure \ref{fig:conf1pand} (right).
     \item \textbf{Conference Day 1}: Interactions between members of a conference on the first day. Only the largest connected components has been selected from this graph. The graph consists of 190 nodes and 703 edges. See Figure \ref{fig:conf1pand} (left). Taken from \url{www.sociopatterns.org/datasets/infectioussociopatterns}.
     \item \textbf{Erd\H{o}s-R\'enyi graph}: Graph sampled from the Erd\H{o}s-R\'enyi random graph model. This graph has 100 nodes and 294 edges.
     \item \textbf{Barab\'asi-Albert graph}: Graph sampled from the Barab\'asi-Albert graph model.  This graph also has 100 nodes and 294 edges. See Fig. \ref{fig:bara_j1}
\end{enumerate}
\begin{figure}[t]
\includegraphics[width=\textwidth]{conf1-and-pandemic}
\caption{\label{fig:conf1pand}Conference Day 1 and Pandemic Network.}
\end{figure}
\subsection{NetShield with $\epsilon$-constraints and GAs}
The results of both the NetShield methods and the randomly initialised GAs are plotted in Figure \ref{fig:res_panconf} for the Pandemic and Conference day graph and in Figure \ref{fig:res_erbara} for the Erd\H{o}s-R\'enyi and Barab\'asi-Albert graph. When comparing the NSGA-II algorithm to the SMS-EMOA algorithm, no clear differences show. It changes from graph to graph which algorithm finds the better Pareto front approximation and they consistently lie very closely together.

Difference do show when comparing the NetShield with the NetShield+ method. Sometimes the difference are large, such as for the Conference day 1 graph and Pandemic graph. For the Barab\'asi-Albert and Erd\H{o}s-R\'enyi graphs the differences are smaller, but the NetShield+ method still tends to give the better results. This is likely due to the Shield-Value losing accuracy when the number of nodes removed increases. Initially the performance of NetShield is very similar to NetShield+. When the allowed cost increases and consequently more nodes can be selected, the NetShield+ method can find solutions that are significantly better.

At the rightmost extremes however, the NetShield method sometimes finds some solutions that dominate those found by the NetShield+ method. See, for example, the Erd\H{o}s-R\'enyi graph and the Barab\'asi-Albert graph. This may be because the NetShield+ method with a batch size of 1 is more greedy than the NetShield method. With a batch size of 1, the NetShield+ method selects at every step the node with the highest eigenscore that would not violate the $\epsilon$-constraint when added to the solution. It then recomputes a new Shield-value function with the node removed. The NetShield method however, only computes the Shield-value function at the beginning and selects multiple nodes at once to optimise this function. In this way it can take a more global view of the problem. If the Shield-Value then happens to still be a good approximation of the
eigen-drop, better solution may be found. A possible approach is to repeat the NetShield+ method several times with different batch sizes if time allows. Then all results can be combined for the most accurate Pareto front approximation.

When comparing both the NetShield methods with the GAs, the GAs give very competitive performance when the networks have relatively few nodes. This means the search space is smaller and the GAs have enough time to converge on the Pareto front. This results in some parts of the Pareto front being approximated better by the GAs, because they can work directly on the eigen-drop. This is most notably the case for the Pandemic graph. Here both NetShield method and the NetShield+ method to a lesser degree have difficulty approximating the Pareto front. This is likely caused by this graph having a low maximum degree. This means that the Shield-value approximation loses accuracy quickly \cite{chen2016node}.

The results also show that the Pareto fronts do not form a single distinctive shape. The Pareto front for the Erd\H{o}s-R\'enyi graph is mostly linear. At the rightmost extreme however, two solutions are found where large gains in eigen-drop can be made with a comparatively small cost increase. This results in the Pareto front having a concave section at the end. This is the opposite for the Pareto front for the Conference day 1 graph. For this graph there is a clear case of diminishing returns: it costs increasingly more to get the same improvement in terms of eigen-drop the further the eigen-drop increases.

The Barab\'asi-Albert graph Pareto front consists of several sections that are mostly linear, but with gaps in between this sections. At these gaps, large increases in eigen-drop are suddenly gained for low costs. This is the result of the preferential attachment model used to generate this graph. At those points the value of $\epsilon$ allows replacing a larger selection of smaller cost nodes with one of the highly connected hub nodes. While these nodes have high cost, their impact on eigen-drop is still disproportionate to their cost. Two solutions with this graph are visualised in figure \ref{fig:bara_j1}: one at the left side of a gap and one at the right side.

\begin{figure}
  \centering
    \includegraphics[width=\textwidth]{Images/panconf_attaintment_netshield.png}
  \caption{Results GAs and NetShield(+) with $\epsilon$-constraint method}
  \label{fig:res_panconf}
\end{figure}

\begin{figure}
  \centering
    \includegraphics[width=\textwidth]{Images/Erdos_Bara_V100_attaintment_netshield.png}
  \caption{Results GAs and NetShield(+) with $\epsilon$-constraint method}
  \label{fig:res_erbara}
\end{figure}

\begin{figure}[t]
  \centering
    \includegraphics[width=\textwidth]{Images/j1}
    \caption{Selected nodes are red and denoted with $\times$. Nodes have been scaled with degree. Left: 6 selected nodes, cost of 36, $\Delta\lambda$ of 0.975.
            Right: 1 selected node, cost of 37, $\Delta\lambda$ of 1.455}
  \label{fig:bara_j1}
\end{figure}

\subsection{Hybrid GA approach}
The results with the GAs initialised from the results from the NetShield methods for the Pandemic and Barab\'asi-Albert graph are shown in Figure \ref{fig:res_atins}. They are shown together with the initialisation sets. The most notable improvements found by the GAs are for the Pandemic graph. Here the results of the inaccuracies of the Shield-value have been corrected. It appears that in these cases, the GAs have the ability to repair such issues.
The improvements for the Barab\'asi-Albert graph are more minor. Either the initialisation sets are already close to the Pareto fronts or there may not be enough diversity in the initial populations for the GAs to find better solutions.
\begin{figure}
  \centering
    \includegraphics[width=\textwidth]{Images/pandbara_attainment_nsinit}
  \caption{Results hybrid GAs}
  \label{fig:res_atins}
\end{figure}
\section{Conclusion}
In this paper, it is shown that the NetShield and NetShield+ algorithms can be extended with an $\epsilon$-constraint method, using an exact quadratic programming solver (here: Gurobi). In this manner, a multi-objective variant of the node immunisation problem can be solved with a cost function added which is proportional to the effort of the node removal. The performance is mostly equivalent or better than two multi-objective genetic algorithms specifically designed for multi-objective optimisation, except for cases where the Shield-value function is inaccurate due to characteristics of the network. In general the NetShield+ method is more robust than the NetShield method, but there are exceptions. Combining the GAs with the NetShield algorithm as initial population only provided small further improvements. Therefore, if time permits and the most accurate Pareto front approximation is required, a valid approach would be to use all methods and combine the end results.
The results also show that there does not appear to be a typical shape to the Pareto fronts resulting from this problem. They are dependent on both the topology of the network and on the cost function.

A main contribution of this paper is the insight that, if time permits, exact and problem-specific methods approximation should be used, which are often far better than Pareto front approximations obtained by general meta-heuristics.
Based on these insight, it will be more effective to develop strategies for controlling real-world networks when the goal is to prevent epidemic outbreaks.
It should be noted, however, that we focused solely on the eigen-value drop. While this is an effective measure under the SIS infection model, the properties of real world epidemics may not be fully captured by this model. Such additional aspects are mentioned, for instance, in \cite{plaat}. In addition, we also assume that any nodes in the network can be immunised.  This also may not translate well to real world scenarios. Therefore, future work should also take a broader view of the problem than further improving the eigen-drop via the process of node removal.

A note of precaution shall be provided here, when applying the model to epidemics in a late stage. A crucial assumption is that the epidemic is in an earlier stage, which will imply that it is likely many neighboring nodes of an infected node are not yet infected. In a late stage of an epidemic this can no-longer be assumed. Moreover, the dynamics are \marios{of the} susceptible-infected-suceptible (SIS) type, and not of the susceptible-infected-recovered (SIR) type, which would be also a common real-world model, e.g., because it would be applicable to the recent SARS-CoV2 pandemic.
It is conjectured, that also in the SIR dynamics the largest eigenvalue is an important value to focus on when dampening or preventing a virus outbreak. When it comes to more realistic and fine grained models of network dynamics, it is probably not straightforward anymore to use a single indicator, such as the eigen-drop, and more complex simulation models are required such as the event-based simulation model in \cite{michalak2020evolutionary}. For an excellent overview various  epidemic models and dynamics on complex networks the reader is referred to \cite{pastor2015epidemic}.

Whereas in our study exact QP solvers were applicable due to the quadratic equations in the objective function, in the more general setting this will no longer be the case and black-box optimization, such as SMS-EMOA will be very useful.
\paragraph{Software*}: All source code (Python) of algorithm implementations and network data of this study is made free available under:\\
\url{https://github.com/joostnibbeling/node-immunisation}

\addcontentsline{toc}{section}{References}

\bibliographystyle{abbrv}


\section{Introduction}

 Let $f(z_1, \ldots, z_n)$ be a holomorphic function defined near $0 \in \mathbf{C}^n$ such that the hypersurface $f= 0$ has an isolated singularity at $0$. In the context of a series of extensive research on such singularities (cf. \cite{T73}, \cite{T77}, \cite{T80}), Teissier considered the \emph{polar invariant} $\theta(f)$ which measures the rate of vanishing of the Jacobian ideal $J_f$ of $f$ with respect to that of the maximal ideal $\mathfrak{m}$ of $0 \in \mathbf{C}^n$. Moreover, he considered the (well-defined) polar invariants $\theta(f_1), \ldots, \theta(f_{n-1})$ where $f_j$ denotes the restriction of $f$ to general $j$-codimensional planes containing $0 \in \mathbf{C}^n$. In this paper, we give the following upper bound for the particular combination of these polar invariants from one of the questions in \cite{T80}  (cf. \cite{T12}, \cite{Li83}, \cite{L84}).

 \begin{thm}\label{main}

  Let $\lct(\mathfrak{m} \cdot J_f)$ be the log canonical threshold at $0 \in \mathbf{C}^n$ of the product ideal $\mathfrak{m} \cdot J_f$. We have

 \begin{equation}\label{main1}
 \frac{1}{1+ \theta (f)} + \frac{1}{1+ \theta (f_1)} + \ldots + \frac{1}{1+ \theta (f_{n-1})} \le  \lct (\mathfrak{m} \cdot J_f).
 \end{equation}

 \end{thm}

In fact, in his question~\cite[p.7]{T80}, Teissier conjectured that one can put the \emph{Arnold exponent} $\sigma(f)$ of $f$ at $0$, in the place of $ \lct (\mathfrak{m} \cdot J_f)$ in \eqref{main1}. The Arnold exponent $\sigma(f)$ is an invariant which is defined from certain natural asymptotic expansions of integrals (cf. \cite{M74}, \cite{AGV84}).     It is related to log canonical thresholds by $\lct (f) = \min \{ \sigma(f), 1 \}$, cf. \cite[\S 9]{Ko97}.  It is also known as (i.e. equal to) the \emph{minimal exponent} which is defined in terms of Bernstein-Sato polynomials, the equality being due to \cite{M75}, \cite{M74}.

  Since $f$ belongs to the integral closure of $\mathfrak{m} \cdot J_f$ in the local ring at $0$ by \cite{BS74} (cf. \cite[Cor. 11.19]{D11}), we have
\begin{equation}\label{lctm}
\lct (\mathfrak{m} \cdot J_f)  \ge \lct (f)
\end{equation} (note that strict inequality holds here for many examples, e.g. $f(x,y) = y^2 + x^3$). In view of these, Theorem~\ref{main} provides an upper bound, weaker than $\sigma(f)$, for the LHS of \eqref{main1} when the hypersurface $(f=0)$ fails to have rational singularities (which is equivalent to $\sigma(f) \ge 1$ by \cite{S93}). When $(f=0)$ does have rational singularities, it will be interesting to further compare $\lct (\mathfrak{m} \cdot J_f)$ and $\sigma(f)$. We also remark that, at least in a certain sense, $\lct (\mathfrak{m} \cdot J_f)$ can be regarded as `not too distant' from the singularity of $f$ itself, in view of  \cite[Prop. 3.8]{ELSV} which says that $J_f$ is included in the multiplier ideal $\mathcal{J}((f)^{1-\epsilon})$ for all $\epsilon > 0$.

 On the other hand, very recently the conjectured bound was shown by B. Dirks and M. Musta\c{t}\v{a} \cite{DM20}, building on the approaches of \cite{EM20}, \cite{L84} which, in particular, reduces the problem to the `one codimension' version. \footnote{This paper was written after we got interested from \cite{EM20} and  before the appearance of \cite{DM20}.}  The method of \cite{DM20} uses, among other things, Saito's theory of mixed Hodge modules and the theory of Hodge ideals (cf. \cite{MP19}, \cite{MP20}). In contrast, our method for the weaker result is completely different and does not use these theories. \footnote{But our method does not deal with the Arnold exponent directly. It will be interesting if the method of this paper can be further combined with the complex analytic aspect of the Arnold exponent = the minimal exponent.}    It uses singularities of plurisubharmonic  functions and  especially the following consequence of a fundamental result of Demailly and Pham~\cite{DH}.

\begin{theorem}

 Let $\varphi$ be a psh function defined near $0 \in \mathbf{C}^n$. Suppose that there exists $C \ge 0$ such that   $ \varphi \ge C \log \abs{z} + O(1)$ as in \eqref{bb}.
 Let $\lct$ be the log-canonical threshold at $0$ and $e_k$ be the $k$-th Lelong number at $0$. Let $L(\varphi)$ be the  Łojasiewicz exponent of $\varphi$ at $0$ (see Definition~\ref{loja}). For each $j= 1, \ldots, n-1$, let $\Lambda_{j} = H_1 \cap \ldots \cap H_{j}$ where $H_1, \ldots, H_j$ are general hyperplanes containing $0 \in \mathbf{C}^n$.
Then we have

\begin{align*}
 \lct (\varphi) &\ge  \frac{e_{n-1} (\varphi)}{e_n (\varphi)} +  \frac{e_{n-2} (\varphi)}{e_{n-1} (\varphi)} + \ldots +  \frac{1}{e_1 (\varphi)} \\
  &\ge  \frac{1}{ L( \varphi ) } + \frac{1}{ L( \varphi|_{\Lambda_1} ) }  + \ldots +  \frac{1}{ L( \varphi|_{\Lambda_{n-1}} ) }.
  \end{align*}

\end{theorem}

 \noindent Here the first inequality is from  \cite{DH} and the second due to Proposition~\ref{eel}. Our arguments for Theorem~\ref{main} use only the `algebraic' case of  plurisubharmonic functions, i.e. those arising from ideals : thus it may be possible that they can be translated in the algebraic language (cf. \cite{BF16}, \cite{B17} for some closely related developments in this regard). We remark that the above main result of \cite{DH} is proved for general psh functions (with isolated singularities) by reducing to the algebraic case (which was further elaborated  in \cite{B17}).
 It is also interesting that our method for Theorem~\ref{main} does not use the `one codimension' version of the main result of \cite{DH}.  In fact, such a statement is an open question :

\begin{question}[Hoang Hiep Pham~\cite{H19}]\label{Pham}
Let $\varphi$ be a psh function with isolated singularities at $0 \in \mathbf{C}^n$. Do we have

$$ \lct(\varphi) \ge \lct_{n-1} (\varphi) + \frac{e_{n-1} (\varphi)}{e_n (\varphi)} $$

\noindent where $\lct_{n-1} (\varphi)$ denotes the supremum of $\lct (\varphi|_H)$ when $H$ ranges over hyperplanes in $\mathbf{C}^n$ containing $0 \in \mathbf{C}^n$ ?

\end{question}

 \noindent We  hope that the interface of analytic and algebraic ideas as in this paper  will help shed light on Question~\ref{Pham}.
 \\

\noindent \textbf{Acknowledgements.}
The author is grateful to Mircea  Musta\c{t}\v{a} and Alexander Rashkovskii for  helpful comments.
This research was supported by Basic Science Research Program through NRF Korea funded by the Ministry of Education (2018R1D1A1B07049683).
\\

\section{Some preliminaries and the proof}

\subsection{Plurisubharmonic singularities}

 We refer to \cite{D11} (and its predecessors) for the definition and basic properties of plurisubharmonic (i.e. psh) functions.
 Let $u$ be a psh function defined near $0 \in \mathbf{C}^n$. We will say that $u$ is singular at $0$ if $u(0) = -\infty$.

Following \cite{D11},  we will say that two psh functions $u$ and $v$ have equivalent singularities if their difference is locally bounded, i.e. if $u = v + O(1)$ where $O(1)$ refers to a function which is locally bounded near every point.   Also we will say that $u$ is less singular (resp. more singular) than $v$ when $u \ge v + O(1)$ (resp. $u \le v + O(1)$).

 In this paper, we will be mostly concerned with psh singularities arising from ideals. When $\mathfrak{a}$ is an ideal (or an ideal sheaf) locally generated by holomorphic functions $f_1, \ldots, f_m$, we will use the notation

 \begin{equation}\label{ideal}
 u = c \log \abs{\mathfrak{a}}
 \end{equation}
 \noindent to refer to the psh function (or more precisely the equivalence class of psh functions under the above relation $u = v + O(1)$) defined by $u = c \log (\sum^m_{j =1} \abs{f_j}^{}) $.  It is well known from \cite{D11} that the singularity equivalence class of $u$ is well-defined by the ideal (or the ideal sheaf), independent of choices of generators $f_1, \ldots, f_m$ (see also \cite[Prop. 3.1]{K14} for some basic exposition). The following is convenient when dealing with integral closure of ideals.

\begin{proposition}\label{closure}

 Suppose that $\mathfrak{a}$ and $\mathfrak{b}$ are two coherent ideal sheaves whose integral closures coincide. Then we have $\log \abs{\mathfrak{a}} = \log \abs{\mathfrak{b}} + O(1)$ for the psh functions defined above as in \eqref{ideal} (for $c=1$).

\end{proposition}

\begin{proof}

 This is well-known, we just recall it explicitly from the literature. It suffices to consider the case when $\mathfrak{b} = \overline{\mathfrak{a}}$. Since $\mathfrak{a} \subset \overline{\mathfrak{a}}$,  we first have $\log \abs{\mathfrak{a}} \le \log \abs{\overline{\mathfrak{a}}} + O(1)$.

 On the other hand, suppose that $\mathfrak{a}$ is generated by $h_1, \ldots, h_m$. From \cite[p.799, Thm. 2.1]{LT08} (cf. \cite[(9.6.10)]{L}), we have $\log \abs{h_j} \le \log \abs{\mathfrak{a}}$. Thus we have  $\log \abs{\overline{\mathfrak{a}}} = \log \sum \abs{h_j} + O(1) = \log \max \abs{h_j} + O(1) \le \log \abs{\mathfrak{a}} + O(1) $.
\end{proof}

 For psh functions, Lelong numbers are  important and fundamental measure of their singularities  (cf. \cite{D11}). When a psh function $u$ has isolated singularities (i.e. locally bounded outside the point $0 \in \mathbf{C}^n$), not only the first (usual) Lelong number $e_1 (u)$, but higher Lelong numbers $e_2 (u), \ldots, e_n(u)$ are defined thanks to work of Demailly (cf. \cite{D93}), which can be expressed as follows:

\begin{equation}\label{kth}
  e_k (u) := e_k (u, 0) = \int_{\{0\}} (dd^c u)^k \wedge (dd^c \log \abs{z} )^{n-k}
\end{equation}

\noindent for $k= 1, \ldots, n$. When $u$ arises from an ideal $\mathfrak{a}$ (which then should be zero-dimensional)  in the sense of \eqref{ideal}, the higher Lelong numbers correspond exactly to mixed multiplicities defined from commutative algebra (cf. \cite{T73}, \cite{L}, \cite{P15}, \cite{B17}). In fact, this holds in the generality of mixed Monge-Ampère masses which include higher Lelong numbers as special cases. For psh functions $u_1, \ldots, u_n$ having isolated singularities at $0 \in \mathbf{C}^n$, we define and denote their (residual) mixed Monge-Ampère mass at $0 \in \mathbf{C}^n$ by

$$ e(u_1, \ldots, u_n) = \int_{\{0\}}  (dd^c u_1) \wedge \ldots \wedge (dd^c u_n) .$$

\noindent We have the following relation between mixed Monge-Ampère mass and mixed multiplicities.

\begin{proposition}

 Let $\mathfrak{a}_1, \ldots, \mathfrak{a}_n \subset \mathcal{O}_{0, \mathbf{C}^n}$ be zero-dimensional ideals at $0 \in \mathbf{C}^n$. For the psh functions $u_k = \log \abs{\mathfrak{a}_k}, k = 1, \ldots, n$, we have  the equality  $$ e(u_1, \ldots, u_n) = \mu (\mathfrak{a}_1, \ldots, \mathfrak{a}_n) $$ where the RHS is the mixed multiplicity of $\mathfrak{a}_1, \ldots, \mathfrak{a}_n$.

\end{proposition}

\begin{proof}

The case when $\mathfrak{a}_1 = \ldots = \mathfrak{a}_n$ is due to  the fundamental result \cite[Lem. 2.1]{D09}. The general case follows from  the polarization (cf. \cite{R11}, \cite{KR18}) as was noted in \cite[Cor. 4.2]{R11}.

\end{proof}

\subsection{Łojasiewicz exponents and polar invariants}

Now we define the Łojasiewicz exponent of a psh function.

\begin{definition}\label{loja}

  Let $\varphi$ be a psh function germ at $0 \in \mathbf{C}^n$ with isolated singularities, i.e. $\varphi$ is locally bounded outside $0$. Suppose that there exists $C \ge 0$ such that

 \begin{equation}\label{bb}
  \varphi \ge C \log \abs{z} + O(1)
  \end{equation}
where  $\abs{z}^2 = \abs{z_1}^2 + \ldots + \abs{z_n}^2$.
   We define the Łojasiewicz exponent of $\varphi$ to be the infimum of such $C$'s and denote it by $L(\varphi)$.

\end{definition}

 This agrees with the more usual Łojasiewicz exponents when $\varphi$ arises from an ideal (cf. \cite[(1.7)]{T77}, \cite{BF16}). An example of a psh function $\varphi$ with isolated singularities but not satisfying the condition in this definition can be found e.g. in \cite[Ex. 2.7]{KR20} (cf. \cite[p.354]{R10}). Note that if $\varphi$ satisfies \eqref{bb}, then so does its restriction to $\Lambda_j$ and thus we can define $L(\varphi|_{\Lambda_j})$.

 \begin{proposition}\label{eel}
 Let $\varphi$ be as in Definition~\ref{loja}.

\begin{enumerate}
\item
 We have $\displaystyle \frac{e_{n-1} (\varphi)}{e_n (\varphi)} \ge \frac{1}{L (\varphi)} $ where $L(\varphi)$ is the Łojasiewicz exponent of $\varphi$.

\item
 For each $j= 1, \ldots, n-1$, let $\Lambda_{j} = H_1 \cap \ldots \cap H_{j}$ where $H_1, \ldots, H_j$ are general hyperplanes containing $0 \in \mathbf{C}^n$.  Then we have
 $$ \frac{e_{n-j-1} (\varphi) }{e_{n-j} (\varphi)} \ge \frac{1}{ L( \varphi|_{\Lambda_j} ) } . $$

 \end{enumerate}

\end{proposition}

\noindent In (2), `general' is in the sense of `almost all' as in \cite{S74} which is used in the following proof.

\begin{proof}

 Let $u := \log \abs{\mathfrak{m}} = \log \abs{z}$.
 In terms of mixed Monge-Ampère masses, we have $e_{n-1} (\varphi) = e (\varphi, \ldots, \varphi, u)$ and $e_n (\varphi) = e(\varphi, \ldots, \varphi, \varphi)$. Whenever $\varphi \ge C u + O(1)$, we have $  e (\varphi, \ldots, \varphi, u) \ge \frac{1}{C} e(\varphi, \ldots, \varphi, \varphi)$ by Demailly's comparison theorem~\cite{D93}. Since $L(\varphi)$ is the infimum of such $C$'s, (1) follows.

 Now when we restrict (1) to $\Lambda_j$, we obtain (2) from \cite{S74}. More precisely, let $T$ and $S$ be the closed positive currents of bidegree $(n-j, n-j)$ given by

 $$ T = (dd^c \varphi) \wedge \ldots \wedge (dd^c \varphi) \wedge (dd^c u), \; \; S = (dd^c \varphi) \wedge \ldots \wedge (dd^c \varphi) \wedge (dd^c \varphi)   .$$

 From the fundamental slicing theory of closed positive currents of \cite[\S 11, p.136]{S74}, the Lelong number of $T$ (resp. of $S$) at $0$ is equal to the Lelong number of its slice by almost every $\Lambda_j$ (i.e. restriction to $\Lambda_j$), a plane of codimension $j$. Here the Lelong number of $T$ (resp. $S$) is nothing but $e(\varphi, \ldots, \varphi, u, u, \ldots, u)$
 (resp. $e(\varphi, \ldots, \varphi, \varphi, u, \ldots, u)$) where $u$ is repeated $j+1$ times (resp. $j$ times).  Hence they are equal to

 $$ e(\varphi|_{\Lambda_j}, \ldots, \varphi|_{\Lambda_j}, u|_{\Lambda_j} )   \quad \text{          resp. }    e(\varphi|_{\Lambda_j}, \ldots, \varphi|_{\Lambda_j}, \varphi|_{\Lambda_j} )  $$

\noindent which are mixed Monge-Ampère masses taken on $\Lambda_j$ (thus taking $n-j$ arguments, of course).
Therefore (1) applied on $\Lambda_j$ yields (2).

\end{proof}

\begin{corollary}\label{polar}

 Let $f$ be as in the beginning of the introduction. For Teissier's polar invariants, we have the following relation when we take $\varphi = \log \abs{\mathfrak{m} \cdot J_f}$, for every $j = 0, 1, \ldots, n-1$,

 $$ \frac{e_{n-j-1} (\varphi) }{e_{n-j} (\varphi)} \ge \frac{1}{1+ \theta (f_{j}) }  $$

\noindent where $ f_0 := f$ and for $j \ge 1$,  $f_j := f|_{\Lambda_j}$ and $\Lambda_{j} = H_1 \cap \ldots \cap H_{j}$ for general hyperplanes $H_1, \ldots, H_j$ containing $0 \in \mathbf{C}^n$.

\end{corollary}

\begin{proof}

 Let us first recall from  \cite[p.270, (1.7) Corollaire 2]{T77}  that, in general, the polar invariant $ \theta (f)$ is equal to the Łojasiewicz exponent of the Jacobian ideal $J_f$, which is equal to $L(\log \abs{J_f})$ in our notation. \footnote{The same argument can be also checked using \cite[Thm. 10.10]{B20} and the definition of $\theta(f)$.}  We will apply this to $f_j$, for which we also need the fundamental fact that the Jacobian ideal of the restriction  $f_j$ and the restriction of the Jacobian ideal $J_f$  coincide up to integral closure, taken from \cite{T73} (cf. \cite[p.21]{P15}). In our notation of psh functions, on $\Lambda_j$, we have (using Proposition~\ref{closure})

 \begin{equation}\label{jfj}
 \log \abs{J_{f_j}} = \psi|_{\Lambda_j} + O(1)
 \end{equation} where $\psi = \log \abs{J_f}$.
  Now since $\varphi =  \log \abs{\mathfrak{m} J_f} = \psi + \log \abs{\mathfrak{m}} + O(1)$, we have $L (\varphi |_{\Lambda_j} ) = L (\psi |_{\Lambda_j})  + L (\log \abs{\mathfrak{m}} |_{\Lambda_j}) = L (\psi |_{\Lambda_j})+ 1$ which is equal to $L(\log \abs{J_{f_j}}) + 1$ by \eqref{jfj}.  We then have $L (\varphi |_{\Lambda_j} ) = 1+ \theta (f_j)$ by the first sentence above : hence  the assertion follows from Proposition~\ref{eel}, (2).
\end{proof}

  Now we complete the proof of Theorem~\ref{main}.

\begin{proof}[Proof of Theorem~\ref{main}]

 Let $J_f  = ( \frac{\partial f}{\partial z_1}, \ldots,   \frac{\partial f}{\partial z_n} )$ be the Jacobian ideal of $f$.  Briançon and Skoda~\cite{BS74} showed (cf. \cite{D11}) that $f$ belongs to the integral closure of the ideal $$(z_1  \frac{\partial f}{\partial z_1}, \ldots, z_n  \frac{\partial f}{\partial z_n} ) \subset (z_1, \ldots, z_n) \cdot J_f =: \mathfrak{a} .$$

Letting $\varphi = \log \abs{\mathfrak{a}}$, we apply the main result of \cite{DH}, Theorem 1.2, to $\varphi$ which has isolated singularities. Together with Corollary~\ref{polar}, it gives

\begin{equation}\label{lct_polar}
 \lct (\varphi) \ge \frac{1}{e_1 (\varphi)} + \frac{e_1(\varphi)}{e_2(\varphi)} + \ldots + \frac{e_{n-1} (\varphi)}{e_n (\varphi)} \ge \sum^{n-1}_{j=0} \frac{1}{1+ \theta (f_{j}) } .
\end{equation}

\noindent Here we note that certainly $e_1 (\varphi) >0 $ holds, which then implies  $e_j (\varphi) >0$ by \cite[Cor. 2.2]{DH} for $j = 2, \ldots, n$.
\end{proof}

\begin{example}

 Let $f(z_1, \ldots, z_n) = z_1^d + \ldots + z_n^d$. Then equality holds in all of \eqref{main1}, \eqref{lctm}, \eqref{lct_polar}. Note that $(z_1^d, \ldots, z_n^d) \subset  \mathfrak{m} \cdot J_f  \subset \mathfrak{m}^d$ and thus the integral closure of $ \mathfrak{m} \cdot J_f $ equals $\mathfrak{m}^d$.

\end{example}

 We remark that further investigations on the equality cases of the above inequality of \cite{DH} (cf. \cite{R15}, \cite{B17}) and of the inequalities in Proposition~\ref{eel} would be of great interest.

\footnotesize

\bibliographystyle{amsplain}