Mixed Integer Programming for Searching Maximum Quasi-Bicliques

Let us introduce several basic notions.

The density of a bipartite graph $G=(V,U,E)$ is $\rho=\frac{|E|}{|V||U|}$.

Let us consider properties and searching algorithms of cliques in a graph $G=(V,E)$.

For a graph $G=(V,E)$ and a fixed $\gamma\in(0,1]$ we need to find a $V^{\prime}\subseteq V$ such that $G[V^{\prime}]$ is a $\gamma$-quasi-clique and $|V^{\prime}|$ is maximal.

Problem of searching for maximum quasi-clique as well as the problem of searching for maximum clique is NP-hard (Liu et al., 2008b),(Pattillo et al., 2013). In addition to that, the assumption of graph incompleteness leads to the loss of useful properties of a clique. For instance, inheritance property which is used in most maximum clique searching algorithms does not hold. Namely, if $G[V]$ is a clique, then $G[V^{\prime}]$ is a clique as well, where $V^{\prime}$ is a subset of $V$. This property does not hold for $\gamma$-quasi-cliques: i.e. a subset of a $\gamma$-quasi-clique is not necessarily a $\gamma$-quasi-clique.

However, for quasi-cliques we can define the property of quasi-inheritance: $\gamma$-quasi-clique with $|V|>1$ is a strict superset to a $\gamma$-quasi-clique with $|V|-1$ vertices (Pattillo et al., 2013).

The problem of maximum quasi-biclique in a bipartite graph $G=(U,~{}V,~{}E)$ with fixed $\gamma\in(0,1]$ is to find $U^{\prime}\subseteq U$ and $V^{\prime}\subseteq V$ such that vertex-induced subgraph $G[U^{\prime},V^{\prime}]$ is a $\gamma$-quasi-biclique of size $|U^{\prime}|+|V^{\prime}|$, maximum for this graph. Let us denote a maximum $\gamma$-quasi-biclique in the graph $G$ by $\omega_{\gamma}\left(G\right)$

Let us consider several commonly met definitions of quasi-biclique. In Liu et al. (2008b), we can find the following definition.

In order to consider the third definition of quasi-biclique Sim et al. (2006), let us introduce some useful notations. The neighbourhood of a vertex $v\in V$ in a graph $G=(V,E)$ is a set of vertices $\Gamma(v)=\{u\in V|(u,v)\in E\}.$

For a vertex set $V^{\prime}\subseteq V$ and a vertex $v\in V\setminus V^{\prime}$ let us denote a set of vertices from $V^{\prime}$ adjacent to $v$ as $\Gamma_{V^{\prime}}(v)=\{u|(u,v)\in E\&u\in V^{\prime}\}$. By a set of vertices $\Gamma(V^{\prime})={\cup_{v\in V^{\prime}}\Gamma(v)}$ we denote a loose neighbourhood of subset $V^{\prime}$

Most properties of quasi-cliques naturally fulfils for quasi-bicliques as well. However, since the density definition of quasi-biclique differs from the case of quasi-clique and the maximum number of edges is a function of two variables with no convex properties, most algorithms searching for maximum quasi-clique are not directly applicable to search for maximum quasi-biclique.

Pattillo et al. (2013) established inequalities for upper bounds for the size of maximum quasi-clique shown below.

In order to obtained similar bound for a quasi-biclique, we need to allow the following conditions on quasi-biclique.

Now let us discuss a few chosen algorithms that implement maximum quasi-biclique search.

A greedy algorithm for searching maximum quasi-bicliques according to Definition 7 is discussed in detail by Liu et al. (2008b). The algorithm uses two parameters: 1) $\delta$ to control the size of the quasi-biclique ($\delta=1-\gamma$) and 2) $\tau$ to control the smallest possible number of vertices that belong to one of the partitions of a quasi-biclique. Let us denote by $U^{\prime}$ and $V^{\prime}$ vertex sets of quasi-biclique of the graph $G(U,V,E)$. At the beginning of the algorithm we set $U^{\prime}=\emptyset$ and $V^{\prime}=V$. From the vertex set $U\setminus U^{\prime}$ we choose such vertex $u$ that its degree is maximum and delete from $V$’ all vertices for which $d(v,U^{\prime})<(1-\delta)\cdot|U^{\prime}|$. This process continues as long as the size of $U^{\prime}<\tau$. However, this algorithm can miss possible vertex candidates, thus authors introduce the second step of the algorithm: if there is a vertex $u$ outside of the current vertex set $U^{\prime}$ such that its degree is maximum in $U\setminus U^{\prime}$ and $U^{\prime}\cup\{u\}$ remains a quasi-biclique, then it can be added to $U^{\prime}$. The same applies to $V^{\prime}$ as long as there is a vertex to add.

In this section we will show how to adapt the model F3 from Veremyev et al. (2016) for searching for maximum quasi-bicliques. Let us consider disjoint sets $U^{\prime}\cup V^{\prime}$, $U^{\prime}\cap V^{\prime}$ = $\emptyset$ that form a quasi-biclique of a bipartite graph $G=(U,V,E)$. Using similar techniques as in Veremyev et al. (2016), we introduce the following variables:

We can refine the sizes of vertex sets of a quasi-biclique using Proposition 2. Then we build Model 1:

Model 1

As in the model F3 we can bound $z_{k}^{(1)}$ and $z_{k}^{(2)}$ and recast them from binary into continuous variables.

Suppose, that there exists an optimal solution $\left(u^{*},~{}v^{*},~{}y^{*},~{}\overline{z^{(1)}},~{}\overline{z^{(2)}}\right)$ of Model 1, where vectors $\overline{z^{(1)}}$ and $\overline{z^{(2)}}$ are not binary $\left(\overline{z_{n}^{(1)}}\geq 0,\ \overline{z_{n}^{(2)}}\geq 0\right)$. Let $\widehat{z^{(1)}}$ be a binary vector with $\widehat{z^{(1)}_{k}}=1\Leftrightarrow|U^{\prime}|=k$ and $\widehat{z^{(1)}_{k}}=0$ otherwise, where $k\in\{\omega_{l}^{(1)},\ldots,\omega_{u}^{(1)}\}$; analogously, vector $\widehat{z^{(2)}}$: $\displaystyle\widehat{z^{(2)}_{k}}=1\Leftrightarrow|V^{\prime}|=k$ and $0$ otherwise. Hence, it is obvious, that vectors $\widehat{z^{(1)}}$ and $\widehat{z^{(2)}}$ satisfy constraint 6. Constraints 3 and 5 can be rewritten as follows:

This means that $\left(u^{*},v^{*},y^{*},\widehat{z^{(1)}},\widehat{z^{(2)}}\right)$ is also an optimal solution of the problem and usage of continuous variables $z^{(1)}_{n}$ and $z^{(2)}_{m}$ in Model 1 is proved.

In the worst case, when $\omega_{l}^{(1)}=\omega_{l}^{(2)}=1$, $\omega_{u}^{(1)}=|U|$, $\omega_{u}^{(2)}=|V|$, the model has $|U|+|V|$ binary variables and $|E|+|U|+|V|$ continuous.

Let us look at different maximizing criteria for related Mixed Integer Programming models. In papers (Ignatov et al., 2015; Mirkin and Kramarenko, 2011) dedicated to triclustering generation, $K=(G,M,B,I)$ is a triadic context with $G$, the set of objects, $M$, the set of attributes, $B$, the set of conditions, and $I\subseteq G\times M\times B$, the ternary relation. The proposed triclustering algorithm searches for clusters that maximize the following criteria:

By narrowing this criteria for binary contexts, it is possible to obtain another maximising criteria for Model 7 GF3($f$) from Veremyev et al. (2016), p.191.

For a bipartite graph $G=(U,V,E)$ and its induced subgraph $G[C_{1},C_{2}]$, function $f$ is maximized over the density and size of biclique.

Using variables definitions from the previous model we can rewrite function $f$:

Since function $f$ is multiplicative, the direct way to transform it to an additive function is logarithmisation:

As in Model 1,

Now we introduce a new variable: $w_{k}=1\Leftrightarrow\lvert\{(i,j):i\in C_{1},j\in C_{2},(i,j)\in E\}\lvert=k$, then $\sum_{(i,j)\in E}{y_{ij}}=\sum_{(i,j)\in E}{kw_{k}}$.

Obviously, that equality $\displaystyle\log{\left(\sum_{(i,j)\in E}{kw_{k}}\right)}=\sum_{(i,j)\in E}{\log(k)w_{k}}$ because $w_{k}$ is binary variable and $\displaystyle\sum_{(i,j)\in E}{w_{k}}=1$. Thus there exists a unique number $k^{*}$ such that $w_{k^{*}}=1$. It follows, that $\displaystyle\log{\left(\sum_{(i,j)\in E}{kw_{k}}\right)}=\log{(k^{*})}=w_{k^{*}}\cdot\log{(k^{*})}=\sum_{(i,j)\in E}{\log(k)w_{k}}$. A similar statement is true for $\displaystyle\log{\left(\sum_{n=\omega_{l}^{(1)}}^{\omega_{u}^{(1)}}{nz_{n}^{(1)}}\right)}$ and  $\displaystyle\log{\left(\sum_{m=\omega^{(2)}_{l}}^{\omega_{u}^{(2)}}{mz_{m}^{(2)}}\right)}$.

Without extra conditions on the sizes of vertex sets of a quasi-biclique and its minimum number of edges, the model has $2\cdot\left(|U|+|V|\right)+2\cdot|E|$ variables.

Model 2

The greedy algorithm of searching for maximal $\gamma$-quasi-biclique in a bipartite graph was implemented in Python 2.7. The MIP models were implemented with the optimization package CPLEX, created by IBM. All computations were carried out on a laptop with macOS operating system, 2.7 GHz Intel Core i5 processor, and RAM 8 GB 1867 MHz.

The search for solutions in the CPLEX package was performed by means of the branch-and-cut method, which is similar to the branch-and-bound algorithmic approach. The method uses a search tree, where each node represents a subproblem that needs to be solved and possibly analysed further.

The branch procedure creates two new nodes from the active parent node. Generally, at this point, the boundaries of one variable are applied and stored for the current node and all its child nodes. In its turn, the cut procedure adds a new constraint to the model. As a result of any cut, the solution space for the subproblems, which are presented in the nodes, is reduced, and the number of branches needed to process decreases. CPLEX processes active nodes in the tree until no more active nodes are available or a certain limit is reached ¹¹1CPLEX user manual: https://www.ibm.com/support/knowledgecenter/SSSA5P_ 12.8.0/ilog.odms.cplex.help/CPLEX/homepages/usrmancplex.html .

The standard solution with the CPLEX software package assumes only one of the optimal solutions as the answer. However, in CPLEX it is possible to obtain a set of optimal solutions using the solution pool method, which allows one to find and store several solutions of MIP models.

The generation of multiple solutions works in two steps. The first step is identical to the usual solution search using the CPLEX software package. At this step, the algorithm finds the only optimal solution of the integer programming problem. It also saves nodes in the search tree that could potentially be useful; for example, if not all the variable constraints are taken into account or if all the nodes contain a suitable value, but the target function is not optimal.

In the second step, using previously calculated and stored information in the first stage, several solutions are generated, and the tree is traversed again, in particular within the branches rooted from the additional nodes stored in the first stage.

On a toy example of a graph with 12 vertices, we consider the search results for maximal $\gamma$-quasi-bicliques, $\gamma=0.8$, using Models 1 and 2, respectively (Fig. 1).

The results for both models are the same (w.r.t. to the solutions output order). Even for this small-sized problem, the time is tangible: the computations with Model 1 took 2.16 s, and for Model 2, it was 2.94 s. A comparison of the executed models and the greedy algorithm in terms of computational time is given below for the selected bipartite graphs.

The algorithm of searching for the maximal quasi-biclique using the CPLEX software package was implemented for Models 1 and 2 (see Section 3) and compared with the greedy algorithm from (Liu et al., 2008b) (let us denote it as Greedy Algorithm).

There are no comparison results presented for the model F3 (Veremyev et al., 2016): despite its fast work, the algorithm based on this model chose quasi-bicliques of very small size and maximum density (i.e. bicliques). This phenomenon is rather expected since the model F3 implies a completely different function of the density of the subgraph. Therefore, the comparison, in this case, is irrelevant. The description of Complete QB in  (Sim et al., 2006) lacks of important implementation details.

The weakness of the constructed MIP models was identified during the finding solution. Since the problem of enumerating all maximal quasi-bicliques in practice requires considerable time, the software package can discard some solutions, if it has found quite a few optimal ones already. First of all, the search is carried out among unbalanced quasi-bicliques (no constraints on the approximately equal size of the quasi-clique partitions have been given). For large graphs, this means that the number of vertices in one of the parts of the found optimal solution may exceed the number of vertices in the second part by hundreds of times or more.

This issue can be addressed in two ways. Firstly, one can set roughly equal limits on the size of the partitions. Secondly, it is possible to adapt the model for finding an almost balanced quasi-bicliques, that means that sizes of partitions of a quasi-biclique differ by $\theta$. To do this, the following conditions should be added to Model 1 or Model 2:

Models with additional conditions 13 and 14 have not been tested.

It has also been noted that small-sized quasi-bicliques can be useless in practice, but their recalculation is costly. Therefore, for each data set, we can establish minimum bounds on the size of a quasi-biclique (of the order of the smallest vertex degree with respect to the graph partitions).

The results the algorithms execution are presented in Table 1 for $\gamma=0.6$, Table 2 for $\gamma=0.7$ and in Table 3 for $\gamma=0.8$²²2The size column in Table 3 shown as the result of summation $|U^{\prime}|$ and $|V^{\prime}|$. For each algorithm its main parameters are indicated: the algorithm running time (time), the number of found maximum quasi-bicliques (count) and the maximum size of the found solution.

Two found $\gamma$-quasi bicliques for the dataset divorce in the US are shown in Fig. 2.

Dashes (”-”) in the following tables mean that the algorithm worked 10 hours and did not find a solution. If one of the partitions of the maximal quasi-biclique has a unit size, this is marked in the table as $(U^{\prime},V^{\prime})$, where $U^{\prime}$ and $V^{\prime}$ are the sizes of the partitions.