The Quadratic Cycle Cover Problem: special cases and efficient bounds

TL;DR

This paper investigates the quadratic cycle cover problem, providing new linearization conditions, bounds, and an iterative algorithm that improves lower bounds and outperforms existing methods.

Contribution

It introduces novel sufficient conditions for linearizability, an iterative bounding procedure, and demonstrates improved bounds over existing methods.

Findings

The iterative bounding algorithm computes the best linearizable matrix.

The classical Gilmore-Lawler bound is shown to be a linearization-based bound.

The new bounds outperform existing bounds in quality and efficiency.

Abstract

The quadratic cycle cover problem is the problem of finding a set of node-disjoint cycles visiting all the nodes such that the total sum of interaction costs between consecutive arcs is minimized. In this paper we study the linearization problem for the quadratic cycle cover problem and related lower bounds. In particular, we derive various sufficient conditions for the quadratic cost matrix to be linearizable, and use these conditions to compute bounds. We also show how to use a sufficient condition for linearizability within an iterative bounding procedure. In each step, our algorithm computes the best equivalent representation of the quadratic cost matrix and its optimal linearizable matrix with respect to the given sufficient condition for linearizability. Further, we show that the classical Gilmore-Lawler type bound belongs to the family of linearization based bounds, and therefore apply the above mentioned iterative reformulation technique. We also prove that the linearization vectors resulting from this iterative approach satisfy the constant value property. The best among here introduced bounds outperform existing lower bounds when taking both quality and efficiency into account.

Tables

Table 1. Table 1: Polyhedral descriptions of the sets τ 1 ​ ( Q ) , τ 2 ​ ( Q ) subscript 𝜏 1 𝑄 subscript 𝜏 2 𝑄 \tau_{1}(Q),\tau_{2}(Q) and τ 3 ​ ( Q ) subscript 𝜏 3 𝑄 \tau_{3}(Q) .

Set	Type of cost matrix	(In)equalities which describe the set
${\tau }_1\left(Q\right)$	Incident weak sum matrix	$b_e+c_f\le Q_{ef}$ ${\hat{p}}_e=b_e+c_e$ $b,c\in {ℝ}^m$	$\forall e\in A,f\in {\delta }^{+}\left(e^{-}\right)$ $\forall e\in A$
${\tau }_2\left(Q\right)$	Restricted generalized weak sum matrix	$b_{ij}+c_{ij,l}+d_{i,jl}+t_{jl}\le Q_{ij,jl}$ $c_{ij,l}+d_{i,kl}\le Q_{ij,kl}$ ${\hat{p}}_{ij}=b_{ij}+{\sum }_{k=1}^n{c}_{ij,k}+{\sum }_{k=1}^n{d}_{k,ij}+t_{ij}$ $b,t\in {ℝ}^m,C\in {ℝ}^{m\times n},D\in {ℝ}^{n\times m}$	$\forall (i,j),(j,l)\in A$ $\forall (i,j),(k,l)\in A,j\ne k$ $\forall (i,j)\in A$
${\tau }_3\left(Q\right)$	Generalized weak sum matrix	$b_{ij,k}+c_{ij,l}+d_{i,kl}+t_{j,kl}\le Q_{ij,kl}$ ${\hat{p}}_{ij}={\sum }_{k=1}^n{b}_{ij,k}+{\sum }_{k=1}^n{c}_{ij,k}+{\sum }_{k=1}^n{d}_{k,ij}+{\sum }_{k=1}^n{t}_{k,ij}$ $B,C\in {ℝ}^{m\times n},D,T\in {ℝ}^{n\times m}$	$\forall (i,j),(k,l)\in A$ $\forall (i,j)\in A$

Table 2. Table 2: Bounds and computation times in seconds of linearization based bounds on Erdős-Rényi instances.

				$LBB1$		$LBB2$		$LBB3$
$p$	$n$	$m$	$OPT$	bound	time	bound	time	bound	time
0.1	20	44	923	923	0.047	923	0.264	923	0.051
	25	76	1039	971	0.005	971	0.477	999	1.227
	30	100	1082	1066	0.013	1066	0.899	1082	0.787
0.3	15	61	485	478	0.010	478	0.635	485	0.714
	20	118	438	377	0.031	384	1.265	390	77.21
	25	172	382	291	0.050	295	2.531	295	869.0
0.5	15	116	226	215	0.034	215	1.407	215	90.39
	20	177	255	189	0.059	190	12.05	190	2105
	25	306	n.a.	172	0.516	173	353.6	n.a.	3600
0.7	15	149	173	127	0.017	128	0.986	128	580.5
	20	263	n.a.	116	0.094	116	7.778	n.a.	3600
	25	396	n.a.	129	0.194	129	18.26	n.a.	3600

Table 3. Table 3: Bounds and computation times in seconds of GL type bounds on Erdős-Rényi instances.

				$GL$		$RG{L}^{sym}$		$RGL$
$p$	$n$	$m$	$OPT$	bound	time	bound	time	bound	time
0.1	20	44	923	923	0.017	923	0.015	923	0.088
	25	76	1039	681	0.039	864	4.652	1018	95.52
	30	100	1082	781	0.053	899	2.412	1082	15.18
0.3	15	61	485	347	0.061	368	1.481	485	7.140
	20	118	438	223	0.068	263	3.482	418	3600
	25	172	382	176	0.136	190	5.105	276	3600
0.5	15	116	226	102	0.336	110	2.602	222	1835
	20	177	255	93	0.118	103	5.365	169	3600
	25	306	n.a.	66	0.296	75	13.80	105	3600
0.7	15	149	173	63	0.080	67	3.530	117	1236
	20	263	n.a.	52	0.181	54	9.624	63	269.3
	25	396	n.a.	56	0.365	62	18.77	79	1085

Table 4. Table 4: Bounds and computation times in seconds of R ​ L ​ T ​ 1 𝑅 𝐿 𝑇 1 RLT1 , C ​ G 𝐶 𝐺 CG , L ​ B ​ B ​ 1 𝐿 𝐵 𝐵 1 LBB1 and R ​ G ​ L 𝑅 𝐺 𝐿 RGL on small Erdős-Rényi instances.

			$OPT$		$RLT1$		$CG$		$LBB1$		$RGL$
$p$	$n$	$m$	value	time	value	time	value	time	value	time	value	time
0.3	20	119	319	10.28	301	4.825	289	102.3	260	0.020	285	3600
	25	177	386	19.04	331	20.09	331	928.9	305	0.040	280	3600
	30	280	n.a.	3600	284	70.62	n.a.	3600	274	0.134	185	3600
0.5	20	195	236	211.0	181	17.00	n.a.	3600	175	0.121	129	3600
	25	327	n.a.	3600	141	82.52	n.a.	3600	136	0.233	89	3600
	30	442	n.a.	3600	168	385.0	n.a.	3600	162	0.322	99	3600

Table 5. Table 5: Bounds and computation times in seconds of R ​ L ​ T ​ 1 𝑅 𝐿 𝑇 1 RLT1 , C ​ G 𝐶 𝐺 CG , L ​ B ​ B ​ 1 𝐿 𝐵 𝐵 1 LBB1 and R ​ G ​ L 𝑅 𝐺 𝐿 RGL on small Manhattan instances.

			$OPT$		$RLT1$		$CG$		$LBB1$		$RGL$
Instance	$n$	$m$	value	time	value	time	value	time	value	time	value	time
$(5,5)$	25	50	103	0.483	103	1.534	103	1.484	103	0.006	103	5.756
$(10,10)$	100	200	418	2.335	418	1.974	418	1645	418	0.022	371	16.06
$(4,4,4)$	64	192	199	6.312	193	9.415	196	691.4	193	0.081	175	3600
$(6,6,6)$	216	648	700	23.67	683	1152	n.a.	3600	683	0.568	551	3600
$(8,8,8)$	512	1536	1566	394.1	n.a.	3600	n.a.	3600	1530	1.213	n.a.	3600
$(10,10,10)$	1000	3000	n.a.	3600	n.a.	3600	n.a.	3600	3113	3.754	n.a.	3600

Table 6. Table 6: Bounds and computation times in seconds of R ​ L ​ T ​ 1 𝑅 𝐿 𝑇 1 RLT1 , C ​ G 𝐶 𝐺 CG , L ​ B ​ B ​ 1 𝐿 𝐵 𝐵 1 LBB1 and R ​ G ​ L 𝑅 𝐺 𝐿 RGL on small Angle-distance instances.

			$OPT$		$RLT1$		$CG$		$LBB1$		$RGL$
$p$	$n$	$m$	value	time	value	time	value	time	value	time	value	time
0.3	20	114	474	2.490	474	1.719	474	416.5	474	0.002	474	64.35
	25	180	553	323.9	553	4.119	n.a.	3600	552	0.004	553	1559
	30	261	512	2951.0	512	19.52	n.a.	3600	512	0.079	494	3600
0.5	20	190	276	177.8	276	6.732	n.a.	3600	276	0.053	274	1319
	25	300	342	2163.6	340	53.43	n.a.	3600	338	0.142	320	3600
	30	435	n.a.	3600	381	490.5	n.a.	3600	377	0.332	355	3600

Table 7. Table 7: Bounds and computation times in seconds of G ​ L 𝐺 𝐿 GL , M ​ I ​ L ​ P 𝑀 𝐼 𝐿 𝑃 MILP , L ​ B ​ B ​ 1 𝐿 𝐵 𝐵 1 LBB1 and L ​ B ​ B ​ 1 s ​ k ​ e ​ w 𝐿 𝐵 𝐵 superscript 1 𝑠 𝑘 𝑒 𝑤 LBB1^{skew} on large Erdős-Rényi instances.

			$GL$		$MILP$		$LBB1$		$LBB{1}^{skew}$
$p$	$n$	$m$	value	time	value	time	value	time	value	time
0.3	30	284	111	0.272	122	0.435	230	0.083	232	0.766
	40	468	117	0.645	131	1.006	265	0.179	278	1.711
	50	754	121	1.598	130	2.410	267	0.404	274	4.184
	60	1062	103	4.068	118	5.788	272	0.726	272	8.048
	70	1481	114	8.910	123	12.94	255	1.660	258	15.38
	80	1842	113	14.26	122	20.82	263	2.740	267	24.67
	90	2385	114	23.25	122	37.61	259	5.296	261	41.74
	100	2962	119	36.03	126	63.49	269	13.32	270	69.90
0.5	30	434	73	0.557	79	0.783	161	0.182	163	9.218
	40	793	69	1.607	74	2.364	166	0.554	169	10.38
	50	1197	72	4.323	77	6.682	165	1.185	167	15.74

Table 8. Table 8: Bounds and computation times in seconds of G ​ L 𝐺 𝐿 GL , M ​ I ​ L ​ P 𝑀 𝐼 𝐿 𝑃 MILP and L ​ B ​ B ​ 1 𝐿 𝐵 𝐵 1 LBB1 on large Manhattan instances.

			$GL$		$MILP$		$LBB1$
Instance	$n$	$m$	value	time	value	time	value	time
$(20,20)$	400	800	1237	5.31	1472	7.491	1537	0.100
$(30,30)$	900	1800	2813	56.24	3343	86.46	3517	0.410
$(40,40)$	1600	3200	5101	346.8	6028	553.5	6302	1.388
$(50,50)$	2500	5000	7983	1225.3	9424	1897.8	9828	2.838
$(17,17,17)$	4913	14739	n.a.	1800	n.a.	n.a.	15398	54.79

Table 9. Table 9: Bounds and computation times in seconds of G ​ L 𝐺 𝐿 GL , M ​ I ​ L ​ P 𝑀 𝐼 𝐿 𝑃 MILP , L ​ B ​ B ​ 1 𝐿 𝐵 𝐵 1 LBB1 and L ​ B ​ B ​ 1 s ​ k ​ e ​ w 𝐿 𝐵 𝐵 superscript 1 𝑠 𝑘 𝑒 𝑤 LBB1^{skew} on Angle-distance instances.

			$GL$		$MILP$		$LBB1$		$LBB{1}^{skew}$
$p$	$n$	$m$	value	time	value	time	value	time	value	time
0.3	30	261	456	0.238	467	0.410	525	0.054	525	6.957
	40	468	507	0.693	516	0.984	567	0.463	567	7.795
	50	735	622	1.567	631	2.223	709	0.317	709	9.982
	60	1062	609	3.684	618	5.401	684	0.694	684	13.97
	70	1449	656	7.436	666	12.37	746	1.331	747	21.63
	80	1896	749	13.03	756	23.77	867	2.613	867	31.96
	90	2403	815	20.33	826	39.99	933	4.838	933	48.81
	100	2970	810	30.35	823	66.04	951	12.50	952	78.62
0.5	30	435	339	0.516	343	0.876	373	0.168	373	16.59
	40	780	411	1.456	418	2.386	464	0.474	464	16.21
	50	1225	466	4.159	473	7.550	534	1.177	535	22.34

Equations

$$\text{OPT(Q)}:=\min$$

$$x\in X,$$

$$X:=\left\{x\in \{0,1{\}}^m|\sum _{e\in {\delta }^{+}(i)}{x}_e=\sum _{e\in {\delta }^{-}(i)}{x}_e=1,\forall i\in N\right\}.$$

$$Q_{e_{i{k}^{\prime }},f_{j{l}^{\prime }}}:=\{\begin{array}{ll}{d}_{k^{\prime }{l}^{\prime }}{w}_{ij}+d_{l^{\prime }{k}^{\prime }}{w}_{ji}& \text{if k′=l′}\\ \infty & \text{otherwise.}\end{array}$$

$$\underset{x\in \{0,1{\}}^m}{\min }\left\{p^{T}x|x\in X\right\},$$

$$\sum _{e\in A}\sum _{f\in {\delta }^{+}(e^{-})}{Q}_{ef}{x}_e{x}_f=\sum _{e\in A}{x}_e{b}_e\sum _{f\in {\delta }^{+}(e^{-})}{x}_f=\sum _{e\in A}{x}_e{b}_e=x^{T}b.$$

$$\sum _{e\in A}\sum _{f\in {\delta }^{+}(e^{-})}{Q}_{ef}{x}_e{x}_f$$

$$=\sum _{e\in A}{b}_e{x}_e\sum _{f\in {\delta }^{+}(e^{-})}{x}_f+\sum _{f\in A}{c}_f{x}_f\sum _{e\in {\delta }^{-}(f^{+})}{x}_e$$

$$=\sum _{e\in A}{b}_e{x}_e+\sum _{f\in A}{c}_f{x}_f=\sum _{e\in A}{p}_e{x}_e.$$

$$x^{T}Qx=x^{T}a{a}^{T}x=(a^{T}x{)}^T(a^{T}x)=(a^{T}x{)}^2.$$

$$x^{T}Qx$$

$$=\sum _{i\in N}\sum _{j\in N_i^{+}}\sum _{k\in N}\sum _{l\in N_k^{+}}{b}_{ij,k}{x}_{ij}{x}_{kl}+\sum _{i\in N}\sum _{j\in N_i^{+}}\sum _{l\in N}\sum _{k\in N_l^{-}}{c}_{ij,l}{x}_{ij}{x}_{kl}$$

$$+\sum _{k\in N}\sum _{l\in N_k^{+}}\sum _{i\in N}\sum _{j\in N_i^{+}}{d}_{i,kl}{x}_{ij}{x}_{kl}+\sum _{k\in N}\sum _{l\in N_k^{+}}\sum _{j\in N}\sum _{i\in N_j^{-}}{t}_{j,kl}{x}_{ij}{x}_{kl}$$

$$=\sum _{i\in N}\sum _{j\in N_i^{+}}{x}_{ij}\sum _{k\in N}{b}_{ij,k}\sum _{l\in N_k^{+}}{x}_{kl}+\sum _{i\in N}\sum _{j\in N_i^{+}}{x}_{ij}\sum _{l\in N}{c}_{ij,l}\sum _{k\in N_l^{-}}{x}_{kl}$$

$$+\sum _{k\in N}\sum _{l\in N_k^{+}}{x}_{kl}\sum _{i\in N}{d}_{i,kl}\sum _{j\in N_i^{+}}{x}_{ij}+\sum _{k\in N}\sum _{l\in N_k^{+}}{x}_{kl}\sum _{j\in N}{t}_{j,kl}\sum _{i\in N_j^{-}}{x}_{ij}$$

$$=\sum _{i\in N}\sum _{j\in N_i^{+}}({\bar{b}}_{ij}+{\bar{c}}_{ij}+{\bar{d}}_{ij}+{\bar{t}}_{ij})x_{ij}=\sum _{i\in N}\sum _{j\in N_i^{+}}{p}_{ij}{x}_{ij},$$

$$B_{ij,k}$$

$$T_{k,ij}$$

$$OPT(p)$$

$$=\underset{y\in {\mathbb{R}}^{2n}}{\max }\{\text{mathbbm}{1}_{2n}^{T}y|[U^T,V^T]y\le p\},$$

$$S_1(Q)$$

$$S_2(Q)$$

$$S_3(Q)$$

$${\tau }_i(Q):=\{\hat{p}\in {\mathbb{R}}^m|x^T\hat{p}=x^T\hat{Q}x\text{ for all }x\in X,\hat{Q}\in S_i(Q)\}.$$

$$OPT(Q)=\underset{x\in X}{\min }\{x^{T}Qx\}\ge \underset{x\in X}{\min }\{x^T\hat{Q}x\}=\underset{x\in X}{\min }\{x^T\hat{p}\}=OPT(\hat{p}).$$

$$v_{LBB}^i:=\underset{\hat{p}\in {\tau }_i(Q)}{\max }\{OPT(\hat{p})\}=\underset{\begin{array}{c}y\in {\mathbb{R}}^{2n}\\ \hat{p}\in {\mathbb{R}}^m\end{array}}{\max }\{\text{mathbbm}{1}_{2n}^{T}y|[U^T,V^T]y\le \hat{p},\hat{p}\in {\tau }_i(Q)\}.$$

$${\tau }_1^{skew}(Q):=$$

$$x^T{Q}_{0}x=x^T{\hat{Q}}_{0}x+x^T(Q_0-{\hat{Q}}_0)x=x^T{p}_1+x^T(Q_0-{\hat{Q}}_0)x$$

$$r_k:=\underset{\begin{array}{c}y\in {\mathbb{R}}^{2n}\\ p_k\in {\mathbb{R}}^m\\ \eta \in [0,1]\end{array}}{\max }\{\text{mathbbm}{1}_{2n}^{T}y|[U^T,V^T]y\le p_k,p_k\in \tau (Q_{k-1}^{\eta })\},$$

$$P_e$$

$$GL$$

$$(G{L}_{ILP})$$

$$\sum _{f\in {\delta }^{+}(i)}{y}_{ef}=\sum _{f\in {\delta }^{-}(i)}{y}_{ef}=x_e$$

$$y_{ee}=x_e$$

$$y_{ef}\in \{0,1\},x\in X$$

$$(DCG{L}_{ILP})\max$$

$${\lambda }_{e,f^{+}}+{\alpha }_{e,f^{-}}\le Q_{ef}$$

$${\lambda }_{e,e^{+}}+{\alpha }_{e,e^{-}}+{\theta }_e\le Q_{ee}$$

$$-\sum _{i\in N}{\lambda }_{e,i}-\sum _{i\in N}{\alpha }_{e,i}+{\gamma }_{e^{-}}+{\mu }_{e^{+}}-{\theta }_e\le 0$$

Full text

In particular, we derive various sufficient conditions for the quadratic cost matrix to be linearizable, and use these conditions to compute bounds. We also show how to use a sufficient condition for linearizability within an iterative bounding procedure. In each step, our algorithm computes the best equivalent representation of the quadratic cost matrix and its optimal linearizable matrix with respect to the given sufficient condition for linearizability. Further, we show that the classical Gilmore-Lawler type bound belongs to the family of linearization based bounds, and therefore apply the above mentioned iterative reformulation technique. We also prove that the linearization vectors resulting from this iterative approach satisfy the constant value property.

The best among here introduced bounds outperform existing lower bounds when taking both quality and efficiency into account.

Keywords: quadratic cycle cover problem, linearization problem, equivalent representations, Gilmore-Lawler bound

1 Introduction

A disjoint cycle cover in a directed graph is a set of node-disjoint cycles such that every node is on exactly one cycle. The quadratic cycle cover problem (QCCP) is the problem of finding a disjoint cycle cover in a graph such that the total sum of interaction costs between consecutive arcs is minimized. Since we assume that all cycle covers in this paper are disjoint, we use the term cycle cover to denote this concept throughout this work. The QCCP is proven to be $\mathcal{NP}$-hard [21]. The corresponding linear problem is called the cycle cover problem (CCP), in which one wants to find a minimum cycle cover with respect to linear arc costs. It is well known that the CCP is solvable in polynomial time.

In the literature several special cases with respect to the objective function of the QCCP are considered. In the angular metric cycle cover problem (Angle-CCP) the quadratic costs represent the change of the direction induced by two consecutive arcs. The goal of Angle-CCP is to find a cycle cover of the graph while minimizing the total angular costs. The Angle-CCP has applications in robotics [4]. In the same paper it is shown that Angle-CCP is $\mathcal{NP}$-hard. Only recently Galbiati et al. [14] introduced another special case of the QCCP: the minimum reload cost cycle cover problem (MinRC3). The MinRC3 problem asks for a minimum cycle cover in an arc-colored graph under the reload cost model. A reload cost is an interaction cost that is paid when two arcs of different colors are placed in succession on a cycle. The goal of the MinRC3 problem is to find a cycle cover such that the total reload cost is minimized. The problem is proven to be $\mathcal{NP}$-hard in the strong sense [14]. The notion of reload costs is introduced by Wirth and Steffan [37], and it has many applications in various fields, e.g. in cargo, energy and telecommunication networks [5, 37]. A detailed overview of the MinRC3 problem and its applications can be found in [7]. Several other combinatorial optimization problems including these reload costs have been investigated [5, 13, 15, 18, 37].

The QCCP is closely related to the quadratic traveling salesman problem (QTSP) which is introduced in [24]. The QTSP is the problem of finding a Hamiltonian cycle in a graph minimizing a quadratic cost function. It has applications in bioinformatics, robotics and telecommunication [20]. When we remove the subtour elimination constraints, the QTSP boils down to the QCCP. Therefore, the QCCP is often used to provide lower bounds for the QTSP [20, 24, 36]. For this reason, the quadratic cycle cover problem is an interesting optimization problem that has received more attention in the past few years.

Several papers have been written about solution methods for the QCCP or its related problems. Jäger and Molitor [24] introduced the QCCP in order to use the QCCP bounds as lower bounds in a branch-and-bound algorithm for the QTSP. Staněk et al. [36] use the QCCP in combination with a rounding procedure to construct heuristics for the QTSP. Aggarwal et al. [4] provide a $\mathcal{O}(\log n)$-approximation algorithm for the Angle-CCP. Fischer [19] studies the polyhedral properties of the QCCP by proving that some triangle inequalities are facet-defining. Galbiati et al. [14] derive various integer programming formulations for the MinRC3 problem. They exploit one of those formulations together with a column generation approach to compute lower bounds for the original problem. Moreover, in [14] a local search algorithm based on 2-exchange and 3-exchange neighbourhoods is constructed to compute upper bounds for the MinRC3 problem. Büyükçolak et al. [7] study the MinRC3 problem on complete graphs with an equitable or nearly equitable 2-edge coloring. For these types of graphs (except some special cases) a polynomial time algorithm is derived that constructs a monochromatic cycle cover.

We focus here on the linearization problem of the QCCP and its applications. An instance of the quadratic cycle cover problem is called linearizable if there exists an instance of the linear cycle cover problem such that the associated costs for both problems are equal for all feasible cycle covers. The linearization problem of the quadratic cycle cover problem asks whether a given instance of the QCCP is linearizable. To the best of our knowledge, this is the first paper about the linearization problem of the QCCP.

In the past few years linearization problems have become an active field of research for many combinatorial optimization problems. In [25, 30] the linearization problem of the quadratic assignment problem (QAP) is studied and polynomial time algorithms that solve it are provided. In particular, Kabadi and Punnen [25] (resp. Punnen and Kabadi [30]) present an ${\mathcal{O}}(n^{4})$ (resp. ${\mathcal{O}}(n^{2})$) algorithm for the general (resp. Koopmans-Beckmann) QAP linearization problem, where $n$ is the size of the problem. The linearization problem for the quadratic minimum spanning tree problem and the quadratic traveling salesman problem are studied by Ćustić and Punnen [11] and Punnen et al. [32], respectively. Hu and Sotirov [23] develop a polynomial time algorithm that solves the linearization problem of the quadratic shortest path problem on directed graphs.

Main results and outline.

In this paper, we first provide an elegant and compact proof that the quadratic cycle cover problem is strongly $\mathcal{NP}$-hard and not approximable within any constant factor unless $\mathcal{P}=\mathcal{NP}$. Then, we consider the linearization problem of the QCCP and derive various sufficient conditions for an instance of the QCCP to be linearizable. In particular, we provide three different types of weak sum conditions on the data matrix for which the corresponding instance can be solved in polynomial time. Further, we present a general framework in which each sufficient condition of linearizability can be used to construct a lower bound on the optimal objective value. Each of these bounds can be computed by a solving a linear programming problem, as long as the set of linearizable matrices is a polyhedron. These types of bounds are called linearization based bounds (LBB), and were recently introduced in [23] for general binary quadratic problems. However, our LBBs exploit sufficient conditions of linearizability suited for the QCCP.

Furthermore, we show how to use a sufficient condition of linearizability within an iterative bounding procedure. In each iteration, we search for the best equivalent representation of the objective and its optimal linearizable matrix that satisfies a particular sufficient condition of linearizability. We refer to the resulting bound as the reformulation based bound (RBB). Our iterative bounding procedure can be seen as a generalization of similar iterative procedures, see e.g., [8, 33, 34].

Finally, we consider the classical Gilmore-Lawler (GL) type bound [16, 28]. First, we show that the GL type bound for the QCCP can be obtained by solving a single linear programming problem instead of solving $m$ (integer) subproblems, where $m$ equals the number of arcs in the graph. Then, we prove that the GL type bound belongs to the family of linearization based bounds by providing the appropriate sufficient condition. We implement our iterative bounding procedure to compute the RBB using the GL type bound. In the literature, iterative approaches for various problems that are based on the GL type bounds use dual variables to obtain bounds, and do not search for equivalent reformulations that provide best bounds in each iteration. Clearly, our approach outperforms others in terms of strength of the bound. Another interesting result is that the linearization vectors resulting from this iterative procedure satisfy the constant value property. Yet another important property for linearizability.

Our numerical results show that the introduced bounding approaches are efficient and provide strong bounds compared to several methods from the literature. In particular, our most prominent bound can be computed within 60 seconds for instances up to 15000 arcs. Interestingly, the GL type bound that is known to be one of the computationally cheapest bounds for quadratic optimization problems cannot be computed on such large instances.

This paper is organized as follows. In Section 2, we formally introduce the QCCP and prove its $\mathcal{NP}$-hardness. In Section 3, the linearization problem for the QCCP is introduced and several sufficient conditions for linearizability are derived. The general framework for the computation of the linearization based bounds is discussed in Section 4. These bounds are used in Section 5 to construct an iterative bounding procedure for each sufficient condition. In Section 6, we consider the classical GL type bound and prove that it belongs to the family of linearization based bounds. We also show how the iterative procedure for this linearization based bound boils down to the computation of the strongest GL type bound in each step. In Section 7, we briefly discuss several other bounds from the literature. Numerical results are given in Section 8.

Notation

A directed graph $G=(N,A)$ is given by a node set $N$ and an arc set $A\subseteq N\times N$. For all nodes $i\in N$ we denote by $\delta^{+}(i)$ the set of arcs that are leaving $i$. Similarly, $\delta^{-}(i)$ denotes the set of arcs that are entering $i$. For all arcs $e\in A$ we let $e^{+}$ and $e^{-}$ denote the starting and ending node of $e$, respectively. To avoid confusion, the letters $e,f$ and $g$ are only used to denote arcs in this work.

For any square matrix $M$, we introduce the operator $\text{diag}:\mathbb{R}^{n\times n}\rightarrow\mathbb{R}^{n}$ that maps a matrix to a vector consisting of its diagonal elements. Moreover, we denote by $\text{Diag}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n\times n}$ its adjoint operator. That is, for any $v\in\mathbb{R}^{n}$ the matrix $\text{Diag}(v)$ equals a diagonal matrix with the entries of $v$ on its main diagonal.

2 The Quadratic Cycle Cover Problem

In this section, we formally introduce the quadratic cycle cover problem.

An instance $\mathcal{I}$ of the QCCP is specified by the pair $\mathcal{I}=(G,Q)$, where $G=(N,A)$ is a directed graph with $n$ vertices and $m$ arcs and $Q=(Q_{ef})\in\mathbb{R}^{m\times m}$ is a quadratic cost matrix. The entries in $Q$ are such that $Q_{ef}=0$ if $f$ is not a successor of $e$. In other words, the quadratic cost of two arcs $e$ and $f$ can be nonzero only if the starting node of $f$ equals the ending node of $e$. In case we also consider linear arc costs, i.e. we have a cost function $p:A\rightarrow\mathbb{R}$, we can put these arc costs on the diagonal of the quadratic cost matrix. Therefore, we assume that the cost structure of an instance of the QCCP is fully determined by its quadratic cost matrix.

Now let $x\in\{0,1\}^{m}$ be a vector with $x_{e}=1$ if arc $e$ belongs to a cycle cover, and 0 otherwise. Then the QCCP can be formulated as

[TABLE]

where $X$ denotes the set consisting of all disjoint cycle covers in $G$, i.e.

[TABLE]

The above set equals the set of directed 2-factors in $G$. For the existence of such a directed 2-factor in a directed graph, see e.g. Chiba and Yamashita [9].

The quadratic cycle cover problem is $\mathcal{NP}$-hard [21]. Also, the related problems Angle-CCP and the MinRC3 problem are shown to be $\mathcal{NP}$-hard [4] and strongly $\mathcal{NP}$-hard [14], respectively. We now provide an alternative reduction that establishes strong $\mathcal{NP}$-hardness which is based on a reduction from the quadratic assignment problem. We consider the Koopmans-Beckmann form of the QAP introduced in [26]. Let $F$ and $P$ be a set of $n$ facilities and $n$ locations, respectively, $w:F\times F\rightarrow\mathbb{R}$ a weight function and $d:P\times P\rightarrow\mathbb{R}$ a distance function. Without loss of generality, we assume that $d_{ii}=w_{ii}=0$ for all $i\in\{1,...,n\}$. Then, we search for a bijection $\pi:F\rightarrow P$ such that $\sum_{i=1}^{n}\sum_{j=1}^{n}d_{\pi(i)\pi(j)}w_{ij}$ is minimized. The QAP is $\mathcal{NP}$-hard in the strong sense and not approximable within any constant factor [35].

Theorem 1**.**

The $QCCP$ is $\mathcal{NP}$-hard in the strong sense and cannot be approximated within a constant factor unless $\mathcal{P}=\mathcal{NP}$.

Proof.

Let an instance $\mathcal{I}$ of the QAP be given, i.e., we consider sets $F=\{1,...,n\}$ and $P=\{1^{\prime},...,n^{\prime}\}$ with $|P|=|F|=n$, functions $w:F\times F\rightarrow\mathbb{R}$ and $d:P\times P\rightarrow\mathbb{R}$ and a positive integer $K$. We create an instance $\mathcal{I}^{\prime}$ of the QCCP.

For the reduction we create a directed graph $G=(N,A)$ that consists of cells. A cell belongs to a single facility and consists of $n$ nodes, each of them corresponding to an assignment to one of the $n$ locations. For each facility $i\in F$, we define a set of $n-1$ identical cells, which we call a group. The nodes corresponding to the same assignment within a group are placed on a directed cycle. In this way, we obtain $n$ cycles per group, which we call inner cycles. We set the interaction cost between each of the successive arcs within a group to zero for all groups. In Figure 1 the group corresponding to facility 1 is given.

We now specify the connections between the groups. Each group is connected to any other group via one of its cells. Since we have $n$ groups and each group consists of $n-1$ cells, this results in ${n\choose 2}$ connections. Connecting the cells of two groups is done by introducing a connection node and a relink node. Starting from the first group, we draw an arc from every node of one of its cells to the connection node. Successively, we draw an arc from the connection node to all the nodes of one of the cells of the second group. The same is done for the relink node, now in the reverse direction. Figure 2 depicts an overview of the connection between the last cell of group $i$ and the first cell of group $j$. We denote the cycles between the groups by outer cycles. In Figure 2 solid arcs are used for the outer cycles, while the inner cycles are drawn using dashed arcs. A similar connection via connection and relink nodes exists for all other pairs of groups.

Observe that any arc in $G$ either belongs to an inner or an outer cycle. The quadratic cost of a pair of successive arcs $(e,f)$ where $e$ belongs to an inner cycle and $f$ to an outer cycle or vice versa, is set to $\infty$. It remains to specify the interaction cost between successive arcs on an outer cycle. We only specify the quadratic cost between the arcs entering and leaving the connection node, other costs are set to zero.

Let $i$ and $j$ be two distinct groups associated with facility $i$ and $j$, respectively. Let a node in group $i$ be given by $(ik^{\prime})$ with $k^{\prime}\in P$. Similarly, a node in group $j$ is given by $(jl^{\prime})$ with $l^{\prime}\in P$. Let $e_{ik^{\prime}}$ denote the arc between $(ik^{\prime})$ and the connection node and let $f_{jl^{\prime}}$ denote the arc between the connection node and $(jl^{\prime})$. Then the quadratic cost between $e_{ik^{\prime}}$ and $f_{jl^{\prime}}$ is defined as follows:

[TABLE]

We repeat this construction for any two connected cells. Figure 3 gives a simplified overview of $G$ for $n=4$. The circles in the center denote the connections between the cells, where the connection and relink nodes are drawn using the symbols ‘$\bullet$’ and ‘$\ast$’, respectively. The graph $G$ has $n^{2}(n-1)+2{n\choose 2}=\mathcal{O}(n^{3})$ nodes and $n^{2}(n-1)+4{n\choose 2}=\mathcal{O}(n^{3})$ arcs.

It remains to show that there exists a cycle cover in $G$ with cost at most $K$ if and only if there exists a feasible assignment in $\mathcal{I}$ with cost at most $K$.

First, we verify that a cycle cover with finite cost in $G$ corresponds to a feasible assignment of facilities and locations. Note that the connection and relink nodes must be covered by an outer cycle, since any other cycle would induce a cost of $\infty$. Besides the connection and relink node, this cycle contains two nodes that each correspond to an assignment of a different facility. Moreover, these facilities must be assigned to different locations, otherwise this implies an infinite cost. The nodes in a cell that are not covered by an outer cycle must be covered by an inner cycle. Consequently, nodes on these inner cycles cannot belong to an outer cycle. Therefore, for each group exactly one location is ‘chosen’ to be on an outer cycle, i.e., each facility is assigned to some location. Moreover, no two facilities are assigned to the same location, since this would imply a cost of $\infty$ at the connection node connecting these groups. We conclude that a cycle cover with finite cost corresponds to a feasible assignment and vice versa.

Observe that the objective value of a feasible assignment in the QAP instance equals the total cost of the corresponding cycle cover in the QCCP instance. Namely, the latter cost equals the sum of quadratic costs incurred at the ${n\choose 2}$ connection nodes. If facility $i$ (resp. $j$) where $i\neq j$ is assigned to location $k^{\prime}$ (resp. $l^{\prime}$) where $k^{\prime}\neq l^{\prime}$, then this cost equals $d_{k^{\prime}l^{\prime}}w_{ij}+d_{l^{\prime}k^{\prime}}w_{ji}$. Taking the sum over all connections, the total cost of the cycle cover equals $\sum_{i=1}^{n}\sum_{j=1}^{n}d_{\pi(i)\pi(j)}w_{ij}$ where $\pi:F\rightarrow P$ is the bijection corresponding to the assignment.

We conclude that there exists an assignment for the QAP instance with cost at most $K$ if and only if there exists a feasible cycle cover in the corresponding QCCP instance of cost at most $K$. Since the QAP is strongly $\mathcal{NP}$-hard and the numbers defined in the reduction are polynomially bounded (infinite costs can be replaced by an appropriate value which is polynomially bounded in the largest number and the size of $\mathcal{I}$), we conclude that the QCCP is strongly $\mathcal{NP}$-hard.

Moreover, as the QAP cannot be approximated within any constant factor [35] and the reduction above is clearly gap preserving, the result follows.

∎

3 The QCCP Linearization Problem

In this section, we formally introduce the linearization problem for the QCCP and derive various sufficient conditions for an instance of the QCCP to be linearizable. Several of these conditions are used later on to construct lower bounds for the optimal value of the QCCP.

Let us consider the (linear) cycle cover problem. Given a cost function $p:A\rightarrow\mathbb{R}$, the CCP is the problem of finding a cycle cover of minimum cost. It can be written as follows:

[TABLE]

where $X$ is given in (2). Since the constraint set of $X$ is totally unimodular, it follows that the CCP is solvable in polynomial time. We call an instance $\mathcal{I}=(G,Q)$ of the QCCP linearizable if there exists a cost vector $p\in\mathbb{R}^{m}$ such that $x^{T}Qx=p^{T}x$ for all cycle covers $x\in X$. If such a vector $p$ exists, we call $p$ a linearization vector of $Q$ for the QCCP.

The QCCP linearization problem can be stated as follows: Given an instance $\mathcal{I}=(G,Q)$ of the QCCP, verify whether it is linearizable and, if yes, compute a linearization vector $p$ of $Q$.

In the remaining part of this section we provide sufficient conditions for the quadratic cost matrix $Q$ to be linearizable. The first type of sufficient conditions for linearizability are related to the constant value property (CVP) for cost vectors or cost matrices. The definition associated with the CCP is stated below.

Definition 1**.**

A cost matrix $p$ satisfies the constant value property if $p^{T}x=p^{T}\bar{x}$ for all cycle covers $x,\bar{x}\in X$.

A similar definition holds for the quadratic version of the problem.

Definition 2**.**

A cost matrix $Q$ satisfies the constant value property if $x^{T}Qx=\bar{x}^{T}Q\bar{x}$ for all cycle covers $x,\bar{x}\in X$.

When $Q$ satisfies the constant value property then $Q$ is linearizable, as stated by the following proposition.

Proposition 1**.**

Assume that $Q$ satisfies the constant value property, i.e. $x^{T}Qx=\xi$ where $\xi\in\mathbb{R}$ for all $x\in X$, then $Q$ is linearizable with cost vector $p$ defined as $p_{e}=\xi/n$ for all $e\in A$.

Proof.

For all $x\in X$ we have $x^{T}Qx=\xi=n\frac{\xi}{n}=x^{T}p$ since $x$ has exactly $n$ nonzero elements. ∎

A more restricted version of the CVP is obtained when the interaction cost of a single arc with its successor or predecessor is constant for all cycle covers $x\in X$. We refer to these properties as the row and column constant value property, respectively. These definitions are based on similar definitions by Punnen et al. [32] for the QTSP.

Definition 3**.**

*A cost matrix $Q$ satisfies the row CVP if there exists some $b\in\mathbb{R}^{m}$ such that for all arcs $e\in A$ we have $Q_{ef}=Q_{eg}=b_{e}$ for all $f,g\in\delta^{+}(e^{-})$ and $Q_{ef}=0$ otherwise.

A cost matrix $Q$ satisfies the column CVP if there exists some $c\in\mathbb{R}^{m}$ such that for all arcs $e\in A$ we have $Q_{fe}=Q_{ge}=c_{e}$ for all $f,g\in\delta^{-}(e^{+})$ and $Q_{fe}=0$ otherwise.*

It is not hard to verify that an instance of the QCCP is linearizable if the cost matrix $Q$ satisfies the row or column CVP.

Proposition 2**.**

If $Q$ satisfies the row CVP or the column CVP, then $Q$ is linearizable.

Proof.

We prove the case when $Q$ satisfies the row CVP. Assume that $b\in\mathbb{R}^{m}$ is such that for all arcs $e\in A$, $Q_{ef}=Q_{eg}=b_{e}$ for all $f,g\in\delta^{+}(e^{-})$ and $Q_{ef}=0$ otherwise. Since $Q_{ef}=0$ when $e$ and $f$ are not successors, we know that $x^{T}Qx=\sum_{e\in A}\sum_{f\in\delta^{+}(e^{-})}Q_{ef}x_{e}x_{f}$. We have

[TABLE]

The proof for the column CVP is similar. ∎

A matrix $Q\in\mathbb{R}^{m\times m}$ is called a sum matrix if there exist $b,c\in\mathbb{R}^{m}$ such that $Q_{ef}=b_{e}+c_{f}$ for all $e,f$. A weak sum matrix is a matrix for which this property holds except for the entries on the diagonal, i.e. $Q_{ef}=b_{e}+c_{f}$ for all $e\neq f$. The weak sum property is used as a condition for linearizability for several quadratic problems, see e.g., [22, 32]. Since in this work we only incur a cost when two arcs are successive, we use a different form of the weak sum condition in which we only put a restriction on successive arcs. We call this condition the incident weak sum property.

Definition 4**.**

A matrix $Q$ is called incident weak sum if there exist vectors $b,c\in\mathbb{R}^{m}$ such that $Q_{ef}=b_{e}+c_{f}$ for all $e\in A$, $f\in\delta^{+}(e^{-})$ and $Q_{ef}=0$ otherwise, i.e. $Q_{ef}=b_{e}+c_{f}$ for all pairs of arcs $e,f$ such that $f$ is a successor of $e$. If such vectors $b$ and $c$ exist, these vectors are called supporting vectors of $Q$.

If the quadratic cost matrix $Q$ is an incident weak sum matrix, then $Q$ is linearizable as stated in the following proposition.

Proposition 3**.**

Let $Q$ be an incident weak sum matrix with supporting vectors $b,c\in\mathbb{R}^{m}$. Then $Q$ is linearizable with cost vector $p=b+c$.

Proof.

We show that for all $x\in X$ we have $x^{T}Qx=p^{T}x$ where $p=b+c$. Note that since $Q_{ef}=0$ for all arcs that are not successors, we have $x^{T}Qx=\sum_{e\in A}\sum_{f\in\delta^{+}(e^{-})}Q_{ef}x_{e}x_{f}$. Now,

[TABLE]

Here we use the fact that $\sum_{f\in\delta^{+}(e^{-})}x_{f}=\sum_{e\in\delta^{-}(f^{+})}x_{e}=1$, since $x$ is a cycle cover. ∎

From Proposition 3 it follows that the incident weak sum property is a sufficient condition for $Q$ to be linearizable. By including linear arc costs, this result remains valid, since we only increase the entries on the diagonal of $Q$.

Moreover, note that when $Q$ satisfies the row or column CVP, then $Q$ is an incident weak sum matrix. Next, we provide a special type of instance for which the cost matrix is not by definition linearizable, but for which we can still obtain its optimal value by solving a linear cycle cover problem.

Definition 5**.**

A matrix $Q\in\mathbb{R}^{m\times m}$ is called a symmetric product matrix if $Q=aa^{T}$ for some vector $a\in\mathbb{R}^{m}$.

Equivalently, we can say that $Q$ is a symmetric product matrix if it is a symmetric positive semidefinite matrix of rank one. Instances with such a quadratic cost matrix can be solved in polynomial time, as stated in the following proposition.

Proposition 4**.**

If the quadratic cost matrix $Q$ is a symmetric product matrix, i.e. $Q=aa^{T}$ for some $a\in\mathbb{R}^{m}$, then the optimal cycle cover can be computed in polynomial time.

Proof.

Let $Q$ be such that $Q=aa^{T}$ for some $a\in\mathbb{R}^{m}$. Then,

[TABLE]

Minimizing $x^{T}Qx$ over all $x\in X$ is then equivalent to minimizing $a^{T}x$ over all cycle covers $x\in X$. ∎

Until now we considered instances for which $Q$ is of the desired type, i.e., $Q_{ef}=0$ when $f$ is not a successor of $e$. Below we derive two sufficient conditions for linearizability of a matrix $Q$ which can have nonzero interaction cost between non-consecutive arcs. Although these cost matrices do not meet the assumptions of the QCCP, we can still use them to derive strong bounds for the objective value of the original problem. This is addressed in Section 4.

Punnen et al. [32] introduce a generalized version of the weak sum property for the QTSP. Their approach can be applied to the QCCP. However, since Punnen et al. [32] prove the condition to hold for complete graphs, we provide a proof for general digraphs.

First, we define some new terminology. Instead of writing $Q_{ef}$ for $e,f\in A$ we can also write $Q_{ij,kl}$ with $(i,j),(k,l)\in A$. Let $N^{+}_{i}$ denote the set of nodes $j$ for which there exists an arc $(i,j)\in A$, i.e., $N^{+}_{i}:=\{j\in N\,|\,\,(i,j)\in A\}$. Similarly, let $N^{-}_{i}$ be the set of nodes $j$ for which an arc $(j,i)\in A$ exists, i.e., $N^{-}_{i}:=\{j\in N\,|\,\,(j,i)\in A\}$. Now we can introduce the notion of a generalized weak sum matrix and prove that it is linearizable.

Definition 6**.**

$Q$* is called a generalized weak sum matrix if there exist $B,C\in\mathbb{R}^{m\times n}$ and $D,T\in\mathbb{R}^{n\times m}$ such that $Q_{ij,kl}=b_{ij,k}+c_{ij,l}+d_{i,kl}+t_{j,kl}$ for all $i,j,k,l$ with $(i,j),(k,l)\in A$. If such $B,C,D$ and $T$ exist, these matrices are called supporting matrices of $Q$.*

Now we can prove the following proposition.

Proposition 5**.**

If $Q$ is a generalized weak sum matrix with supporting matrices $B,C\in\mathbb{R}^{m\times n}$ and $D,T\in\mathbb{R}^{n\times m}$, then $Q$ is linearizable with cost vector $p$ given by $p_{ij}=\sum_{k=1}^{n}b_{ij,k}+\sum_{k=1}^{n}c_{ij,k}+\sum_{k=1}^{n}d_{k,ij}+\sum_{k=1}^{n}t_{k,ij}$.

Proof.

Let $\bar{b}_{ij}:=\sum_{k=1}^{n}b_{ij,k}$, $\bar{c}_{ij}:=\sum_{k=1}^{n}c_{ij,k}$, $\bar{d}_{ij}:=\sum_{k=1}^{n}d_{k,ij}$, $\bar{t}_{ij}:=\sum_{k=1}^{n}t_{k,ij}$ and $p_{ij}=\bar{b}_{ij}+\bar{c}_{ij}+\bar{d}_{ij}+\bar{t}_{ij}$ for all $(i,j)\in A$. Then, for $x\in X$,

[TABLE]

where we use the fact that $\sum_{l\in N^{+}_{k}}x_{kl}=\sum_{k\in N^{-}_{l}}x_{kl}=\sum_{j\in N^{+}_{i}}x_{ij}=\sum_{i\in N^{-}_{j}}x_{ij}=1$ since $x$ is a cycle cover. ∎

Note that an incident weak sum matrix can be seen as a special case of a generalized weak sum matrix. That is, for all $(i,j)\in A$ we set $b_{ij,j}=b_{ij}$ and $b_{ij,k}=0$ for all $k\neq j$ and for all $(k,l)\in A$ we set $t_{k,kl}=t_{kl}$ and $t_{j,kl}=0$ for all $j\neq k$. Moreover, let $C$ and $D$ be the zero matrix. Then we have $Q_{ij,jl}=b_{ij,j}+c_{ij,l}+d_{i,jl}+t_{j,jl}=b_{ij}+t_{jl}$ for all $(i,j),(j,l)\in A$ and $Q_{ij,kl}=0$ otherwise.

When $Q$ is a generalized weak sum matrix, we need $4mn$ parameters to describe $Q$. This number can be reduced by considering a more restricted version of a generalized weak sum matrix.

Definition 7**.**

$Q$* is called a restricted generalized weak sum matrix if there exist $C\in\mathbb{R}^{m\times n}$, $D\in\mathbb{R}^{n\times m}$ and $b,t\in\mathbb{R}^{m}$ such that $Q_{ij,jl}=b_{ij}+c_{ij,l}+d_{i,jl}+t_{jl}$ for all $i,j,l$ with $(i,j),(j,l)\in A$ and $Q_{ij,kl}=c_{ij,l}+d_{i,kl}$ otherwise. If such $C,D$ and $b,t$ exist, these are called supporting matrices and vectors, respectively.*

We can show that restricted generalized weak sum matrices are linearizable.

Proposition 6**.**

If $Q$ is a restricted generalized weak sum matrix with supporting matrices $C,D$ and supporting vectors $b,t$, then $Q$ is linearizable with vector $p$ given by $p_{ij}=b_{ij}+\sum_{k=1}^{n}c_{ij,k}+\sum_{k=1}^{n}d_{k,ij}+t_{ij}$.

Proof.

Define $B\in\mathbb{R}^{m\times n}$ and $T\in\mathbb{R}^{n\times m}$ as follows:

[TABLE]

Now the matrices $B,C,D$ and $T$ are such that they satisfy the conditions of Proposition 5. This implies that $Q$ is linearizable with vector $p^{\prime}$ given by $p^{\prime}_{ij}=\sum_{k=1}^{n}b_{ij,k}+\sum_{k=1}^{n}c_{ij,k}+\sum_{k=1}^{n}d_{k,ij}+\sum_{k=1}^{n}t_{k,ij}$. Since $\sum_{k=1}^{n}b_{ij,k}=b_{ij}$ and $\sum_{k=1}^{n}t_{k,ij}=t_{ij}$ it follows that $Q$ is linearizable with vector $p:=b_{ij}+\sum_{k=1}^{n}b_{ij,k}+\sum_{k=1}^{n}d_{k,ij}+t_{ij}$. ∎

4 Linearization Based Bounds for QCCP

In this section we show how the sufficient conditions for linearizability can be used to derive bounds for the optimal value of the QCCP. The construction of these bounds is provided in Section 4.1. Section 4.2 shows some preliminary numerical results of these bounding procedures.

4.1 Construction of Linearization Based Bounds

When an instance of the QCCP is linearizable, we can solve the problem in polynomial time by solving the corresponding linear cycle cover problem. When a quadratic cost matrix $Q$ is not linearizable, we can still use the sufficient conditions for linearizability to find lower bounds for the optimal value of the problem. This approach is introduced by Hu and Sotirov [23] for general binary quadratic problems. We here use tailor made sufficient conditions for the QCCP, which leads to efficient lower bounds as we show later in the numerical results.

Before we proceed, let us recall the linear cycle cover problem. We introduce the matrix $U\in\mathbb{R}^{n\times m}$ with $U_{ie}=1$ if node $i$ is the starting node of arc $e$ and 0 otherwise. Similarly, we define $V\in\mathbb{R}^{n\times m}$ with $V_{ie}=1$ if node $i$ is the ending node of arc $e$ and 0 otherwise. Since the matrix $[U^{T},V^{T}]^{T}$ is totally unimodular, the optimal value of the CCP using cost vector $p$ equals

[TABLE]

where $\mathbbm{1}_{2n}\in\mathbb{R}^{2n}$ equals the vector of ones. Note that (3) and (4) are equivalent optimization problems. The corresponding dual problem is given in (5).

When $Q$ is linearizable with linearization vector $p$, we can find the optimal value for the QCCP by computing $OPT(p)$ using (4) or (5). If $Q$ is not linearizable, we can search for a linearizable matrix $\hat{Q}$ that is as close as possible to $Q$. To guarantee that $\hat{Q}$ is indeed linearizable, it should satisfy one of the sufficient conditions for linearizability derived in Section 3. We define the sets $S_{i}(Q)$, for $i=1,2,3$, consisting of cost matrices $\hat{Q}$ such that $\hat{Q}$ is linearizable w.r.t. a sufficient condition for linearizability and $Q-\hat{Q}$ is elementwise nonnegative. We have

[TABLE]

Remark 1**.**

We do not consider the sets of cost matrices $Q$ satisfying the row or column CVP, since these are special types of incident weak sum matrices. These type of matrices are contained in $S_{1}$.

$S_{i}(Q)$ can be seen as the set of all the linearizable cost matrices of a specific type that are suitable for obtaining lower bounds for the optimal value of the problem. For this purpose, we define for $i=1,2,3$ the set $\tau_{i}(Q)$ of cost vectors $\hat{p}\in\mathbb{R}^{m}$ as

[TABLE]

It is clear that for all $i$ and all $\hat{p}\in\tau_{i}(Q)$ we have

[TABLE]

So, indeed, $OPT(\hat{p})$ is a lower bound for the optimal objective value of the QCCP for all $i$ and $\hat{p}\in\tau_{i}(Q)$. By maximizing over all cost vectors in $\tau_{i}(Q)$, we obtain the strongest linearization based bound with respect to the set $S_{i}(Q)$, which we denote by $v_{LBB}^{i}$, see also [23]:

[TABLE]

The corresponding bounding approaches are denoted by $LBB1$, $LBB2$ and $LBB3$, respectively.

Remark 2**.**

Recall that the matrices in $S_{2}(Q)$ and $S_{3}(Q)$ can have nonzero interaction cost for non-consecutive arcs, so they do not satisfy the assumptions on the quadratic cost matrix of the QCCP. Nevertheless, they can still be used to derive lower bounds for the original problem. In other words, we search for the general quadratic cost matrix that is linearizable and as close as possible to $Q$ that gives us the best lower bound.

As long as the set $\tau_{i}(Q)$ is a polyhedron, the corresponding bound $v^{i}_{LBB}$ can be calculated by solving the linear programming problem (6). The sets $\tau_{i}(Q)$ for $i=1,2,3$ are indeed nonempty polyhedra, since they can be described by a finite number of linear equalities and inequalities. These polyhedral descriptions are provided in Table 1.

By construction, we have $\tau_{1}(Q)\subseteq\tau_{2}(Q)\subseteq\tau_{3}(Q)$ for all quadratic cost matrices $Q$. Consequently, we can establish the following result about the quality of the corresponding linearization based bounds.

Theorem 2**.**

For all instances of the QCCP, we have $v^{1}_{LBB}\leq v^{2}_{LBB}\leq v^{3}_{LBB}$.

Proof.

The proof follows immediately from the construction of the linearization based bounds and the definitions of the incident weak sum matrix, the restricted generalized weak sum matrix and the generalized weak sum matrix. ∎

Hu and Sotirov [23] argue that the linearization based bounds can be improved by extending the sets $\tau_{i}(Q)$ using a skew symmetric matrix $M$. That is, since each skew symmetric matrix is linearizable, a matrix $\hat{Q}$ is linearizable if and only if $\hat{Q}+M$ is linearizable for all $M$ with $M+M^{T}=0$. Using this, the set $\tau_{1}(Q)$ can be extended to:

[TABLE]

Note that in $\tau^{skew}_{1}(Q)$ we only include skew symmetric matrices whose support corresponds to the pairs of successive arcs in $G$, since adding dense skew symmetric matrices would increase computational complexity. Since $\tau_{1}(Q)\subseteq\tau_{1}^{skew}(Q)$, it follows that we can obtain a stronger bound by maximizing over $\tau_{1}^{skew}(Q)$, see Section 8. The same extension can be applied to any set $\tau_{i}(Q)$.

4.2 Preliminary Results

In order to check the quality of the bounds derived above, we perform a preliminary numerical study. We create instances according to the $G(n,p)$ Erdős-Rényi model [12]. Here $n$ equals the number of nodes and $p$ equals the probability that an arc is included. We create instances for various values of $n$ and $p$. The interaction cost between any two successive arcs is drawn uniformly at random as an integer from $\{1,...,100\}$. In Table 2 we present the bounds $v_{LBB}^{1},v_{LBB}^{2}$ and $v_{LBB}^{3}$ and the computational times required for computing them. Experiments are performed using a pc with an Intel(R) Core(TM) i5-6500 CPU, 3.20 GHz and 8 GB memory using CPLEX 12.6 as the solver. The maximum computation time is set to 3600 seconds and we put ‘n.a.’ in the table when this maximum is reached before a solution is obtained.

By construction, the optimal solution has always an integer objective value. Therefore, we round up all bounds to the nearest integer. The results of Table 2 show that the linearization based bounds $LBB1$, $LBB2$ and $LBB3$ do not differ significantly, especially for the larger instances. At the same time, the computation times differ significantly. It turns out that $LBB1$ is most efficient. Therefore, this bound can be preferred when taking both quality and efficiency into account.

5 Reformulated LBB Approach

In this section we discuss how a reformulation of the quadratic cost matrix can be used to obtain a non-decreasing sequence of lower bounds that are based on the linearization based bound. It is important to note that one can construct such a bounding procedure using any sufficient condition for linearizability, not only the ones discussed in Section 4.

Suppose we are given a sufficient condition for linearizability. Let $S(Q)$ and $\tau(Q)$ be as in Section 4, but now for a general sufficient condition. That is, $S(Q)$ is the set consisting of all linearizable cost matrices $\hat{Q}$ of this type with $\hat{Q}\leq Q$ and $\tau(Q)$ consists of the corresponding linearization vectors. Moreover, we assume that the set $\tau(Q)$ is a polyhedral set. Let $Q_{0}$ be the initial quadratic cost matrix. If $\hat{Q}_{0}\in S(Q_{0})$, we know there exists some $p_{1}\in\tau(Q_{0})$ such that $x^{T}\hat{Q}_{0}x=x^{T}p_{1}$ for all $x\in X$. This leads to the following reformulation of the objective function

[TABLE]

for all cycle covers $x\in X$. By letting $p_{0}$ be the $m\times 1$ zero vector and $Q_{1}:=Q_{0}-\hat{Q}_{0}$, this relation can be written as $x^{T}p_{0}+x^{T}Q_{0}x=x^{T}p_{1}+x^{T}Q_{1}x$ for all $x\in X$. The vector $p_{1}$ is taken to be the largest linearization of the matrix $Q_{0}$ (see (6)) with $Q_{1}$ being the residual quadratic part. Recall that the bound $v_{LBB}$ is calculated by only taking the linear part of this new objective function into account. By construction we have $x^{T}p_{0}\leq x^{T}p_{1}$ for all $x\in X$ and $0\leq Q_{1}\leq Q_{0}$.

Now we can proceed in a similar way by considering the linearization problem of the residual cost matrix $Q_{1}$. In other words, we search for a linearizable matrix in $S(Q_{1})$ and its corresponding linearization vector $p_{2}$. If we let $Q_{2}$ be the new residual matrix, then the objective function can be reformulated as $x^{T}p_{0}+x^{T}Q_{0}x=x^{T}(p_{1}+p_{2})+x^{T}Q_{2}x$. Since $x^{T}p_{1}\leq x^{T}(p_{1}+p_{2})$, a stronger bound can be obtained. This procedure can be repeated to obtain a sequence of bounds. However, it is in general not possible to find a vector $p_{2}$ for which this bound has strictly improved. That is, since $p_{1}+p_{2}\in\tau(Q_{0})$ this would imply that $p_{1}$ is not the optimal solution to (6). Thus, applying this procedure iteratively, the resulting sequence of bounds remains constant after the first iteration. To overcome this issue, we need to reformulate the residual cost matrix in each step.

In the literature, various iterative bounding procedures have been proposed, see e.g., [6, 31, 33, 34]. In this paper we introduce a new approach that is different in two ways. First, the existing bounding procedures are mainly based on the classical Gilmore-Lawler type bound. Our approach is based on general sufficient conditions for linearizability and we can show that the Gilmore-Lawler type bounding procedure is a special case of this approach, see Section 6. Second, the existing bounding procedures mostly use a fixed reformulation of the cost matrix in each iteration. However, using a fixed reformulation is in general not the best one can do. Here, we search for the reformulation of the cost matrix that results in the strongest bound in the next iteration. For this purpose, we define the notion of an equivalent representation of a matrix, see e.g. [31].

Definition 8**.**

Let $(G,Q)$ be an instance of the QCCP. Then $(G,R)$ is an equivalent representation of $(G,Q)$ if $x^{T}Qx=x^{T}Rx$ for all $x\in X$.

If there is no confusion about the graph $G$ under consideration, we say that $R$ is an equivalent representation of $Q$. It is easy to verify that a matrix $R$ is an equivalent representation of $Q$ if for all $e,f\in A$ we have $R_{ef}+R_{fe}=Q_{ef}+Q_{fe}$. Here, we focus on a specific type of equivalent representation, which we call an $\eta$-equivalent representation of $Q$.

Definition 9**.**

Given $\eta\in[0,1]$, an equivalent representation $Q^{\eta}:=\eta Q+(1-\eta)Q^{T}$ of $Q$ is called an $\eta$-equivalent representation.

In other words, an $\eta$-equivalent representation is obtained by taking a convex combination of $Q$ and $Q^{T}$. Moreover, it follows that if $R$ and $Q$ are equivalent representations, then a linearization of $R$ is also a linearization of $Q$ and vice versa.

Instead of considering the linearization problem of the residual matrix $Q_{1}$, we can consider the linearization problem of $Q_{1}^{\eta}$ for some $\eta\in[0,1]$. Since $Q_{1}^{\eta}$ has a different structure than $Q_{1}$, it is in general possible to find a linearizable matrix $\hat{Q}_{1}\in S(Q_{1}^{\eta})$ and a corresponding linearization vector that result in a strictly stronger bound.

As already mentioned above, many approaches in the literature are based on taking a fixed value for $\eta$, e.g. $\eta=\frac{1}{2}$ which corresponds to the case of symmetrizing. This does not give the best bound in general. Instead, we search for $\eta\in[0,1]$ that results in the strongest bound in each iteration. Suppose we are in step $k$ of the algorithm in which we consider the linearization problem of the residual matrix $Q_{k-1}$. Then the optimal equivalent representation of $Q_{k-1}$ and its corresponding vector $p_{k}\in\tau(Q^{\eta}_{k-1})$ can be computed simultaneously by solving the following problem:

[TABLE]

which equals the additional amount of quadratic cost that can be linearized in iteration $k$. Note that if the set $\tau(Q_{k-1})$ is a polyhedron, then $\tau(Q_{k-1}^{\eta})$ is also a polyhedron and the corresponding problem (9) can be solved in polynomial time. For the sufficient conditions mentioned in Section 4 this is indeed the case.

Finally, we provide a new bounding procedure that is based on iteratively finding the best equivalent representation of the residual cost matrix and its optimal linearizable matrix. Starting with $Q_{0}=Q$, the goal is to find the best linearizable matrix $\hat{Q}_{k-1}$ of an equivalent representation of the residual matrix $Q_{k-1}$ and its corresponding linearization vector $p_{k}$. In each iteration we compute $r_{k}$ by (9) and let $d_{k}=d_{k-1}+p_{k}$ which equals the total linearization vector of $Q$. The final bound is given by the sum of all $r_{k}$’s, which we call the reformulation based bound. The procedure is given in Algorithm 1.

Remark 3**.**

Note that steps 3 and 4 of Algorithm 1 depend on the specific sufficient condition for linearizability. For instance, for the incident weak sum condition we construct in step 4 the linearizable matrix $\hat{Q}_{k-1}$ in the following way $(\hat{Q}_{k-1})_{ef}:=b_{e}+c_{f}$ for all $e\in A,f\in\delta^{+}(e^{-})$ and $(\hat{Q}_{k-1})_{ef}:=0$ otherwise, where $b,c\in\mathbb{R}^{m}$ are obtained from (9).

Hu and Sotirov [23] show that in the case that the linearizable matrix $\hat{Q}$ is in the form $\hat{Q}=[U^{T},V^{T}]Y+\text{Diag}(z)$ for some $Y\in\mathbb{R}^{2n\times m}$ and $z\in\mathbb{R}^{m}$, the bound $v_{RBB}$ is dominated by the solution of the first level RLT relaxation introduced by Adams and Sherali [2, 3]. Here RLT stands for reformulation linearization technique. In [23] it is moreover shown that the first level RLT bound, denoted by $v_{RLT1}$, can be obtained by searching for the optimal linearizable matrix $\hat{Q}$ of the form $\hat{Q}=[U^{T},V^{T}]Y+M+\text{Diag}(z)$ where $M$ is a skew symmetric matrix.

Our preliminary numerical results show that the above algorithm does not improve significantly the $LBB1$ bound. However, in the next section we show that our approach outperforms known iterative approaches related to the Gilmore-Lawler bounds.

6 The Gilmore-Lawler Type Bound

In this section we consider the classical Gilmore-Lawler type bound. The GL procedure is a well-known approach to construct lower bounds for quadratic 0-1 optimization problems, see e.g. [16, 28, 33, 34]. We provide a compact formulation of the GL type bound that can be used to compute the bound by a single LP-problem, instead of solving $m$ subproblems. Moreover, we show that this bound in fact belongs to the family of linearization based bounds. Therefore, based on the results of Section 5, we provide a bounding procedure that computes the best GL type bound in each step of the algorithm. We conclude this section by testing this new bounding procedure on some preliminary test instances.

6.1 The classical GL type bound

In the objective function of the QCCP, see (1), we have the quadratic term $x_{e}x_{f}$ for each two arcs $e,f\in A$ placed in succession on a cycle. To get rid of this quadratic term, for each given arc $e\in A$, potentially in the solution, we consider the cycle cover containing $e$ with minimum interaction cost with $e$. We denote this minimum contribution of arc $e$ to a solution by $z_{e}$. In particular, for all $e\in A$ we have

[TABLE]

where $Q_{e,:}$ denotes the $e$th row of the cost matrix $Q$. The feasible set of $(P_{e})$ equals the set of all node-disjoint cycle covers containing arc $e$. If this set is empty, then we set $z_{e}$ equal to 0 since arc $e$ cannot contribute to a cycle cover.

Let $z\in\mathbb{R}^{m}$ be the vector consisting of the elements $z_{e}$ for all $e\in A$. Then the classical GL type bound is obtained by solving the following CCP:

[TABLE]

Note that the constraint matrices of $(P_{e})$ and $(GL)$ are totally unimodular. For this reason, we can drop the integrality constraints and solve $(GL)$ and $(P_{e})$ for all $e\in A$ as linear programming problems.

Besides computing the GL type bound by solving $(GL)$ and $(P_{e})$ for all $e\in A$, we can also obtain its value by solving an integer linear programming (ILP) problem. The problem $(GL_{ILP})$ is defined as follows:

[TABLE]

It follows from the constraints that if $x_{e}=1$, then $y_{e,:}:=[y_{e1},...,y_{em}]$ is the characteristic vector of the cheapest cycle cover containing arc $e$ and if $x_{e}=0$, then $y_{e,:}$ equals the zero vector.

Let $(CGL_{ILP})$ be the continuous relaxation of $(GL_{ILP})$. In this continuous relaxation we can omit the upper bounds on $x_{e}$ and $y_{ef}$ for all $e,f\in A$, since it is never optimal to set the value of these variables larger than one. We can compute the GL type bound by solving $(CGL_{ILP})$ as stated by the following theorem. This theorem is based on a similar result for the QMST, see [33].

Theorem 3**.**

The optimal value of $(CGL_{ILP})$ equals $v_{GL}$.

Proof.

Let $\lambda_{e,i}$ and $\alpha_{e,i}$ be the dual variables corresponding to constraints (10), i.e. $\lambda_{e,i}$ corresponds to $\sum_{f\in\delta^{+}(i)}y_{ef}=x_{e}$ and $\alpha_{e,i}$ corresponds to $\sum_{f\in\delta^{-}(i)}y_{ef}=x_{e}$. Similarly, let $\mu_{i}$ and $\gamma_{i}$ be the dual variables corresponding to the first and second equalities of constraints (2), and $\theta_{e}$ the dual variable corresponding to constraints (11). The dual problem of $(CGL_{ILP})$ is as follows:

[TABLE]

Constraint (15) can be rewritten as $\gamma_{e^{-}}+\mu_{e^{+}}\leq\sum_{i\in N}\lambda_{e,i}+\sum_{i\in N}\alpha_{e,i}+\theta_{e}$ for all $e\in A$. In order to maximize the objective function of $(DCGL_{ILP})$, we maximize the right hand side of this inequality subject to constraints (13)-(14). This gives for each $e\in A$ the following subproblem:

[TABLE]

For each fixed $e\in A$ the subproblem given above equals the dual of the continuous relaxation of $(P_{e})$. By strong duality we know $z^{\prime}_{e}=z_{e}$ for all $e\in A$. Substitution of this term into constraint (15) gives a rewritten formulation for $(DCGL_{ILP})$:

[TABLE]

This problem exactly equals the dual of the continuous relaxation of $(GL)$. Because of strong duality, it follows that the optimal objective value of $(CGL_{ILP})$ equals $v_{GL}$. ∎

We can show that the Gilmore-Lawler type bound for the QCCP in fact belongs to the family of linearization based bounds introduced in Section 4. That is, we can obtain $v_{GL}$ by searching for a linearizable quadratic cost matrix $\hat{Q}$ of a specific type that is as close as possible to $Q$. The required linearizability condition on $\hat{Q}$ is given below, and it differs from the sufficient conditions presented in Section 3.

Proposition 7**.**

If there exists $B,C\in\mathbb{R}^{m\times n}$ and $t\in\mathbb{R}^{m}$ such that $Q_{ef}=B_{e,f^{+}}+C_{e,f^{-}}$ for $e\neq f$ and $Q_{ee}=B_{e,e^{+}}+C_{e,e^{-}}+t_{e}$ for all $e\in A$, then $Q$ is linearizable with vector $p$ given by $p_{e}=t_{e}+\sum_{i=1}^{n}B_{e,i}+\sum_{i=1}^{n}C_{e,i}$.

Proof.

Let $\tilde{Q}$ be defined as $\tilde{Q}_{ef}=B_{e,f^{+}}+C_{e,f^{-}}$ for all $e,f\in A$. Then $\tilde{Q}$ can be seen as a generalized weak sum matrix where $D$ and $T$ are equal to the zero matrix, see Definition 6. According to Proposition 5, $\tilde{Q}$ is linearizable with vector $\tilde{p}=\sum_{i=1}^{n}B_{e,i}+\sum_{i=1}^{n}C_{e,i}$. Since $Q=\tilde{Q}+\text{Diag}(t)$, it follows that $Q$ is linearizable with vector $p$ given by $p_{e}=t_{e}+\sum_{i=1}^{n}B_{e,i}+\sum_{i=1}^{n}C_{e,i}$. ∎

Similar to the notation used in Section 4, let $S_{GL}(Q)$ denote the set of all linearizable cost matrices $\hat{Q}\in\mathbb{R}^{m\times m}$ that satisfy the conditions of Proposition 7. Moreover, let $\tau_{GL}(Q)$ be the following polyhedron:

[TABLE]

and

[TABLE]

Now we prove the main result of this section which states that the classical Gilmore-Lawler type bound can be seen as a special case of linearization based bound.

Theorem 4**.**

We have $v_{LBB}^{GL}=v_{GL}$.

Proof.

By using the polyhedral description of $S_{GL}(Q)$ following from Proposition 7, the optimization problem in (17) can be written as follows:

[TABLE]

We show that this optimization problem is equivalent to $(DCGL_{ILP})$, the dual problem of the continuous relaxation of $(GL_{ILP})$. Take $B_{e,i}=\lambda_{e,i}$ and $C_{e,i}=\alpha_{e,i}$ for all $e\in A$ and $i\in N$, where $\lambda$ and $\alpha$ denote the dual vectors belonging to constraints (10). Similarly, let $t=\theta$ where $\theta$ equals the dual vector to constraints (11). Finally, let $y=[\mu^{T},\gamma^{T}]^{T}$, where $\mu$ and $\gamma$ are the dual variables belonging to constraints (2). By substitution of these variables and combining constraints (19) and (22), we obtain the problem $(DCGL_{ILP})$, i.e., the dual of $(CGL_{ILP})$. Thus, we have $v_{LBB}^{GL}=v_{GL}$. ∎

Theorem 4 shows that the GL type bound belongs to the family of linearization based bounds. This is also shown by Hu and Sotirov [23] and Rostami et al. [34], however our proof is very different as we exploit the fact that $v_{GL}$ can be obtained by solving an LP problem, i.e., $(CGL_{ILP})$. Additionally, we show here that the computation of the GL type bound is equivalent to the search for the optimal linearizable cost matrix $\hat{Q}$ satisfying the properties of Proposition 7.

6.2 The best Gilmore-Lawler type bound

Section 6.1 shows that the calculation of the classical GL type bound fits in the general framework discussed in Section 4. In this section we apply the reformulation procedure of Section 5 to the GL type bound. We also show that our approach outperforms several iterative approaches from the literature.

In order to apply Algorithm 1 to the sufficient condition for linearizability of Proposition 7, we need to define how to calculate $r_{k}$ for each iteration $k$, see (9). We rewrite the set $\tau_{GL}(Q)$, see (16), as follows:

[TABLE]

which is clearly a polyhedron. Then for all $k\geq 1$ we calculate the additional amount of quadratic cost that is linearized by solving:

[TABLE]

Note that as opposed to the constraints in (9), we have replaced the constraint $[U^{T},V^{T}]y\leq p_{k}$ by an equality constraint. This does not change the value of $r_{k}$. To verify this, suppose we solve (24) using the inequality constraint $[U^{T},V^{T}]y\leq p_{k}$ and let $\hat{y},\hat{p}_{k}$ and $\hat{t}$ be the corresponding optimal solutions. Let $e\in A$ be such that the inequality constraint is satisfied with strict inequality. Then without changing $\hat{y}$, we can reduce $\hat{t}_{e}$ (and thus $\hat{p}_{e}$) such that we get equality for $e\in A$. Although it changes the linearization vector $\hat{p}$, the resulting bound remains equal. To verify this, notice that only the left hand side of constraint (21) is decreased, so the solution is still feasible and the optimal value $r_{k}$ remains unchanged. From this, it follows that one may replace $[U^{T},V^{T}]y\leq p_{k}$ by an equality constraint and solve $r_{k}$ as in (24).

Algorithm 1 using (24) to compute $r_{k},p_{k}$ and $\eta$ gives a new bounding procedure for the QCCP. We call the resulting bound the reformulated GL type bound ($RGL$) and denote its value by $v_{RGL}$. By construction, it iteratively computes the best Gilmore-Lawler type bound among all equivalent representations of the quadratic cost matrix.

The algorithm proposed in this section satisfies another interesting property, namely the vectors $d_{k}$ satisfy the constant value property for all $k\geq 0$. This is an important property for linearizability because the set of linearizable cost matrices for combinatorial optimization problems with interaction costs can be characterized by the constant value property, under certain conditions, see [27].

Theorem 5**.**

All $d_{k}$ where $k\geq 0$ computed during the $RGL$ approach, satisfy the constant value property, i.e., we have $x^{T}d_{k}=\bar{x}^{T}d_{k}$ for all feasible cycle covers $x,\bar{x}\in X$.

Proof.

We apply a proof by induction on $k$. Note that the vector $d_{0}$ equals the $m\times 1$ vector of zeros which trivially satisfies the constant value property.

Now assume that the induction hypothesis is true for iteration $k-1$, i.e., $x^{T}d_{k-1}=\bar{x}^{T}d_{k-1}$ for all feasible cycle covers $x,\bar{x}\in X$. In iteration $k$ we solve (24). Let $\hat{y}\in\mathbb{R}^{2n}$ and $\hat{p}\in\mathbb{R}^{m}$ be the optimal variables for this problem and split $\hat{y}=[\mu^{T},\lambda^{T}]^{T}$ with $\mu,\lambda\in\mathbb{R}^{n}$. It follows that $[U^{T},V^{T}]\hat{y}=U^{T}\mu+V^{T}\lambda=\hat{p}$. Now let $x\in X$ be any feasible cycle cover in $G$. Then we can sum up the rows of this system of equalities for all arcs $e\in A$ in the cycle cover implied by $x$:

[TABLE]

where we use the fact that each vertex is visited exactly once on a cycle cover. So the quantity $x^{T}\hat{p}$ is equal for all $x\in X$. As a result, $\hat{p}$ satisfies the constant value property.

The vector $d_{k}$ is constructed as $d_{k-1}+p_{k}$ with $p_{k}=\hat{p}$. Since $d_{k-1}$ and $\hat{p}$ satisfy the constant value property, it follows that $d_{k}$ satisfies the constant value property. ∎

Remark 4**.**

Since the GL type bound can be computed both as a linearization based bound and by solving $(CGL_{ILP})$ (see Theorem 4), the iterative approach derived in this section can also be defined in terms of $(CGL_{ILP})$. In that case, we iteratively compute $v_{GL}$ and reformulate the quadratic cost matrix using the dual variables of $(CGL_{ILP})$. The details of this equivalent approach can be found in [29].

Since the linearizable matrix $\hat{Q}$ of Proposition 7 can be written in the form $\hat{Q}=[U^{T},V^{T}]Y+\text{Diag}(z)$ for some $Y\in\mathbb{R}^{2n\times m}$ and $z\in\mathbb{R}^{m}$, it follows from [23] that we have $v_{RGL}\leq v_{RLT1}$.

6.3 Preliminary Results

For the instances considered in the preliminary results of Section 4, we now test our Gilmore-Lawler type bounds. First, we compute the classical GL type bound, after symmetrizing the quadratic cost matrix $Q$. This bound is denoted by $GL$. Moreover, we consider the iterative GL type bounding approach where we symmetrize the quadratic cost matrix in each iteration. That is, we apply Algorithm 1 using (24) where instead of optimizing over $\eta$, we set $\eta=\frac{1}{2}$. We denote this bound by $RGL^{sym}$. Finally, we report the bound $RGL$ which is introduced in Section 6.2. The maximum computation time is set at 3600 seconds. The results are given in Table 3.

From Table 3 it follows that the iterative approaches significantly improve the classical GL type bound. Among these iterative approaches, $RGL$ provides much stronger bounds than $RGL^{sym}$. We conclude that this new approach of calculating the best GL type bound in each step provides better bounds than when setting $\eta=\frac{1}{2}$ in the reformulation. However, it turns out that this improvement in the quality comes at the cost of efficiency. Clearly, we can stop our algorithm after a pre-specified number of iterations and/or time.

7 Other bounds for the QCCP

In this section we present several known bounding approaches from the literature that can be applied to the QCCP. In the next section, we compare those bounds with the bounds introduced in this paper. We consider a column generation approach and a mixed integer linear programming (MILP) based bound.

Galbiati et al. [14] construct a column generation approach for the MinRC3. This approach can be extended to the QCCP. To the best of our knowledge, this is the only implemented lower bounding approach for the MinRC3 in the literature.

Let $\mathcal{C}$ be the set of all directed cycles in $G$. Moreover, let $\overline{\mathcal{C}}\subseteq\mathcal{C}$ be a subset of cycles such that it contains at least one cycle cover. Let $w_{c}$ be the cost of a cycle $c\in\mathcal{C}$. Then the restricted master problem $(RMP)$ is given by:

[TABLE]

Let $\pi\in\mathbb{R}^{n}$ be the vector of dual variables corresponding to constraint (25). Then the subproblem $(SP)$ searches for the cycle in $\mathcal{C}$ with the smallest (negative) reduced costs with respect to $\pi$, i.e.

[TABLE]

where $z_{i}=1$ if vertex $i$ is on the cycle and 0 otherwise. As stated in [14], the problem $(SP)$ is strongly $\mathcal{NP}$-hard itself. The quadratic objective function can be linearized by standard linearization techniques. A lower bound on the optimal value of the QCCP can be obtained by iteratively solving the master problem and its corresponding subproblem. If a cycle with negative reduced cost is found, we add it to the set $\overline{\mathcal{C}}$. This procedure is repeated until no more cycle with negative reduced cost is found or after some predefined stopping criteria. The obtained bound is denoted by $v_{CG}$.

Based on a procedure by [17, 1], we present the QCCP as an MILP problem. Let us first fix an equivalent representation of $(G,Q)$. Let $z_{e}$ be computed as in $(P_{e})$ for all $e\in A$, see Section 6.1. Moreover, we define for all $e\in A$

[TABLE]

Note that $q_{e}^{\max}$ can be obtained by solving a linear programming problem. Then the QCCP can be formulated as an MILP:

[TABLE]

If we relax the binary constraint on $x$, then we obtain a lower bound for the QCCP. We call this bound the MILP based bound and we denote its value by $v_{MILP}$. The next result shows that $v_{MILP}$ is at least as large as the Gilmore-Lawler type bound.

Theorem 6**.**

The MILP based bound dominates the Gilmore-Lawler type bound, i.e., $v_{GL}\leq v_{MILP}$.

Proof.

Let $\beta,\delta\in\mathbb{R}^{m}_{+}$ denote the dual variables of (27) and (28), respectively. Moreover, let $\mu,\gamma\in\mathbb{R}^{n}$ denote the dual variables of the cycle cover constraints $\sum_{e\in\delta^{+}(i)}x_{e}=1$ and $\sum_{e\in\delta^{-}(i)}x_{e}=1$ for all $i\in N$, respectively. Then the dual of the MILP based bound equals

[TABLE]

This problem equals the dual of the continuous relaxation of $(GL)$. Hence, it follows that $v_{GL}\leq v_{MILP}$. ∎

Note the MILP based bound and the Gilmore-Lawler type bound are comparable if the same equivalent reformulation of $(G,Q)$ is used in their computations.

8 Numerical Results

In this section we test our bounding approaches on a set of test instances and compare them with several approaches from the literature. We take into account the linearization based bound $LBB1$ from Section 4.1, the classical GL type bound $GL$ from Section 6.1, the reformulated GL type bound $RGL$ discussed in Section 6.2, the column generation approach $CG$ and the MILP based bound $MILP$ from Section 7, and the first level RLT bound $RLT1$, see [2, 3]. The GL bound and the MILP based bound are computed after symmetrizing $Q$. Note that we do not take into account $LBB2$ and $LBB3$, since our preliminary experiments from Section 4.2 show that $LBB1$ is preferred when taking both quality and efficiency into account.

All lower bounds are implemented in Matlab on a pc with an Intel(R) Core(TM) i5-6500 CPU, 3.20 GHz and 8 GB memory using CPLEX 12.6 as the solver.

We consider the following types of instances:

• •

Erdős-Rényi instances: These instances are created via the $G(n,p)$ Erdős-Rényi model [12]. The number of nodes is fixed to $n$ and each arc is included with probability $p$ independent of the other arcs. The quadratic cost between any pair of successive arcs is chosen discrete uniformly at random out of $\{0,...,100\}$.

• •

Manhattan instances: The Manhattan instances are introduced in [10] and resemble the street pattern of modern cities like Manhattan or Barcelona. Given a finite set of positive integers $(n_{1},n_{2},...,n_{k})$, the graph consists of a $n_{1}\times n_{2}\times...\times n_{k}$ directed grid. Each node in the interior is adjacent to its $2k$ neighbours. The nodes on the boundary are also incident to the corresponding nodes on the opposite boundary. For each dimension $k$, the arcs belonging to the same layer of the grid point in the same direction. However, the arcs of two consecutive layers point in the opposite direction. This results in a graph containing a large number of cycles. The quadratic cost between any pair of successive arcs is chosen discrete uniformly at random out of $\{0,...,10\}$.

• •

Angle-distance instances: The Angle-distance instances are originally constructed for the QTSP in [19]. The number of nodes $n$ and the graph density $p$ are given. The $(x,y)$-coordinates of each node is chosen discrete uniformly at random out of $\{0,...,500\}^{2}$. Exactly $\lceil pn(n-1)\rceil$ arcs are chosen uniformly at random from the total set of arcs. For each arc $e\in A$, let $d_{e}$ denote the Euclidean distance between the endpoints of $e$. Moreover, for each two successive arcs $e$ and $f$, let $\alpha_{ef}$ denote the turning angle (in radians). Given some constant $\rho\in\mathbb{R}_{+}$, the quadratic cost of two successive arcs $e$ and $f$ is calculated as:

[TABLE]

Similar as in [19], we take $\rho=40$.

For Erdős-Rényi and Angle-distance instances, preliminary experiments show that instances up to approximately 300 arcs can be solved to optimality within one hour. For the Manhattan instances the limit is around 2000 arcs, due to the small density of these types of graphs.

In total we consider two sets of experiments: experiments on small instances and experiments on large instances. Since the optimum, $RLT1$ and $CG$ cannot be calculated for large instances, we only test these approaches on the smaller instances. Moreover, we include the bounds introduced in this paper, namely $LBB1$ and $RGL$. The value and computation times (in seconds) on small Erdős-Rényi instances can be found in Table 4. This table contains 6 instances for $n=20,25,30$ and $p=0.3,0.5$. The results on Manhattan and Angle-distance instances are reported in Tables 5 and 6, respectively. For the Angle-distance instances we take the same values for $n$ and $p$ as for the Erdős-Rényi instances, while for the Manhattan instances we consider several two- and three-dimensional instances. The maximum computation time is set to 3600 seconds. When after this time no bound is computed, we report ‘n.a.’ in the tables. Since the optimal value is always integer, we round up all bounds.

For the smaller instances, we see that $RLT1$ performs best in quality. When it can be computed, it is often close to the optimal value and it dominates the other bounds. $LBB1$ is often very close to $RLT1$, but can be computed much more efficiently. Namely, for the Erdős-Rényi and the Angle-distance instances the computation time of $LBB1$ for all small instances is below 0.4 seconds, whereas $RLT1$ cannot be computed within one hour for some of these instances. The column generation approach provides strong bounds, but in most cases it is not able to compute a bound in a time span of one hour. The reformulated GL type bound performs well on the Manhattan and Angle-distance instances, see Tables 5 and 6. Although its total computation time is large, the advantage of this approach is that it provides a bound in a short time and then iteratively improves the value. This makes it possible to stop the procedure at any given time and still obtain a bound. The bounds computed by $RGL$ are in almost all cases dominated by $LBB1$.

When taking both efficiency and quality into account, we conclude that the linearization based bound $LBB1$ outperforms the other approaches. Based on Tables 4, 5 and 6, the value of $LBB1$ is at least 75% of the optimal value for the Erdős-Rényi instances. For the Angle-distance and Manhattan instances, this percentage equals 98% and 96%, respectively.

For the larger instances, we only compute the bounds that can be computed efficiently. That is, we do not consider the iterative approaches, but only the bounds $GL$, $MILP$ and $LBB1$. We also investigate the effect of a reformulation by adding an optimal incident skew symmetric matrix to the cost matrix, see Section 4.1. We apply this reformulation to $LBB1$, which implies that we optimize over the set $\tau_{1}^{skew}(Q)$, see (7), instead of $\tau_{1}(Q)$. The resulting bound is denoted by $LBB1^{skew}$. For the Manhattan instances this bound is omitted, since preliminary experiments showed that this reformulation does not improve the bounds for most Manhattan instances. This could be due to the sparsity of Manhattan instances. The bounds and computation times (in seconds) for the Erdős-Rényi, Manhattan and Angle-distance instances are reported in Tables 7, 8 and 9, respectively. For the Erdős-Rényi and Angle-distance instances we take for $n$ values between 30 and 100 nodes and consider $p=0.3$ and $p=0.5$. For the Manhattan instances we consider large two-dimensional instances and one large three-dimensional instance. The maximum computation time for these bounds is set to 1800 seconds. Again, we round up all bound values.

For the larger instances, we see that $LBB1$ in all cases dominates $GL$ and $MILP$ in both quality and efficiency. The difference in quality is most present for the Erdős-Rényi instances, see Table 7. For the Manhattan instances, we see that $GL$ and $MILP$ can be calculated efficiently for instances up to 3000 arcs. However, $LBB1$ remains efficient even for larger instances. In particular, bounds for Manhattan instances up to 15000 arcs can be computed within 60 seconds.

Moreover, we conclude from Tables 7 and 9 that the addition of an incidence skew symmetric matrix to the set $\tau_{1}(Q)$ only improves the bounds for some of the instances. In general, it turns out that the Erdős-Rényi instances can successfully be improved by this method, whereas for the Angle-distance instances only in a few cases there is an improvement. Although the computation times of $LBB1^{skew}$ are larger than those of $LBB1$, bounds can still be computed in a reasonable time span.

9 Conclusion

In this paper we consider the linearization problem for the QCCP and its applications. We provide several sufficient conditions for linearizability, and show how these conditions can be used to obtain strong lower bounds for the QCCP. The linearization based bound $LBB1$, resulting from the incident weak sum property, is the most efficient LBB in terms of complexity and quality, see Table 2. We show here that the GL type bound for the QCCP also belongs to the family of linearization based bounds, see Theorem 4, by providing the appropriate sufficient condition, see Proposition 7.

The first level RLT bounds and/or the GL type bounds are the only linearization based bounds for quadratic binary optimization problems that are implemented for various binary quadratic optimization problems up to date. This paper shows that besides these two well-known bounds, the linearization based bounds introduced here are worth considering.

Here, we also present how each sufficient condition can be used in an iterative bounding procedure. In particular, we introduce a new reformulation technique in which we search for the best equivalent representation of the residual cost matrix and its optimal linearizable matrix, see Algorithm 1. We show how the resulting iterative procedure computes the best GL type bound in each iteration. Our approach outperforms known iterative bounding procedures that use the GL type bounds, see Table 3. Moreover, we prove that the resulting linearization vectors in each step satisfy the constant value property, see Theorem 5.

Finally, our numerical results show that our approach outperforms several other bounds from the literature if we take into account both quality and efficiency. Although the linearization based bounds $LBB1$ are dominated by the well known first level RLT bounds, they can be computed extremely fast. For the Manhattan instances, $LBB1$ bounds for instances up to 15000 arcs can be computed within 60 seconds. However, other approaches fail to provide bounds for instances of this large sizes.

We expect that similar bounding procedures can be successfully applied for other quadratic optimization problems, such as the quadratic assignment problem, the quadratic minimum spanning tree, and the quadratic traveling salesman problem. However, this is a subject of our future research.

Acknowledgements

We would like to thank two anonymous reviewers for their insightful comments and suggestions to improve an earlier version of this work.