A Method for Convex Black-Box Integer Global Optimization

TL;DR

This paper introduces a new convex underestimator method for integer black-box optimization that does not require gradient information, enabling global minimization on finite integer domains.

Contribution

It proposes a novel secant-based underestimator and an algorithm that alternates between updating this underestimator and evaluating the objective, with proven convergence to the global minimum.

Findings

The method effectively underestimates the objective in disconnected regions.

The algorithm converges to the global minimum on the feasible set.

Comparative results show computational advantages over existing methods.

Abstract

We study the problem of minimizing a convex function on a nonempty, finite subset of the integer lattice when the function cannot be evaluated at noninteger points. We propose a new underestimator that does not require access to (sub)gradients of the objective but, rather, uses secant linear functions that interpolate the objective function at previously evaluated points. These linear mappings are shown to underestimate the objective in disconnected portions of the domain. Therefore, the union of these conditional cuts provides a nonconvex underestimator of the objective. We propose an algorithm that alternates between updating the underestimator and evaluating the objective function. We prove that the algorithm converges to a global minimum of the objective function on the feasible set. We present two approaches for representing the underestimator and compare their computational effectiveness. We also compare implementations of our algorithm with existing methods for minimizing functions on a subset of the integer lattice. We discuss the difficulty of this problem class and provide insights into why a computational proof of optimality is challenging even for moderate problem sizes.

Tables

Table 1. Table 1: Number of primitive directions, # = | 𝒩 ​ ( x c , 1 ) | # 𝒩 subscript 𝑥 𝑐 1 \#=\lvert\mathcal{N}(x_{c},1)\rvert , that emanate from the origin x c subscript 𝑥 𝑐 x_{c} of the domain Ω = [ − k , k ] n ∩ ℤ n Ω superscript 𝑘 𝑘 𝑛 superscript ℤ 𝑛 {\Omega=\left[-k,k\right]^{n}\cap\mathbb{Z}^{n}} and that correspond to points in Ω Ω \Omega .

	$n=2$		$n=3$		$n=4$		$n=5$
$k$	$\left|\Omega \right|$	$\#$	$\left|\Omega \right|$	$\#$	$\left|\Omega \right|$	$\#$	$\left|\Omega \right|$	$\#$
1	9	8	27	26	81	80	243	242
2	25	16	125	98	625	544	3,125	2,882
3	49	32	343	290	2,403	2,240	16,807	16,322
4	81	48	729	578	6,561	5,856	59,049	55,682

Table 2. Table 2: Description of how SUCIL versions choose X 𝑋 X in ( PLP ) and the next iterate x k + 1 superscript 𝑥 𝑘 1 x^{k+1} (breaking ties in arg ​ min arg min \operatorname*{arg\,min} arbitrarily). G k superscript 𝐺 𝑘 G^{k} is defined in Definition  3 , and X k superscript 𝑋 𝑘 X^{k} is all points evaluated before iteration k 𝑘 k .

Method	$X$ in (PLP)?	$x^{k+1}=$?
SUCIL	$G^k$	$\underset{x\in {\Omega }_k:{\Vert x-\hat{x}\Vert }_{\infty }\le {\Delta }_k}{\arg \min }{\eta }_k\left(x\right)$
SUCIL-noTR	$G^k$	$\underset{x\in {\Omega }_k}{\arg \min }{\eta }_k\left(x\right)$
SUCIL-ideal1	$X^k$	$\underset{x\in \Omega \setminus X^k}{\arg \min }f\left(x\right)$
SUCIL-ideal2	$G^k$	$\underset{x\in \Omega \setminus X^k}{\arg \min }f\left(x\right)$

Table 3. Table A1: Set of convex test problems.

Name	Expression	$f\left(x^{\ast }\right)$	$x^{\ast }$
CB3II [11]	$\max \{\sum _{i=1}^{n-1}{x}_i^4+x_{i+1}^2,\sum _{i=1}^{n-1}{\left(2-x_i\right)}^2+{\left(2-x_{i+1}\right)}^2,\sum _{i=1}^{n-1}2e^{-x_i+x_{i+1}}\}$	$2\left(n-1\right)$	$e$
CB3I [10]	$\sum _{i=1}^{n-1}\max \{x_i^4+x_{i+1}^2,{\left(2-x_i\right)}^2+{\left(2-x_{i+1}\right)}^2,2e^{-x_i+x_{i+1}}\}$	$2\left(n-1\right)$	$e$
KLT [29]	$\underset{i\in \{1,\dots ,n\}}{\max }\left\{{\Vert x-c_i-2e\Vert }^2\right\},c_i=2e_i-e$	$n$	$2e$
LQ [22]	$\sum _{i=1}^{n-1}\max \{-x_i-x_{i+1},-x_i-x_{i+1}+x_i^2+x_{i+1}^2-1\}$	$-\left(n-1\right)$	many
abhi [1]	$\sum _{i=1}^n\left[64{\left(c_1\left(x_i-2\right)-c_2\left(x_{i+1}-2\right)\right)}^2+{\left(c_2\left(x_i-2\right)-c_1\left(x_{i+1}-2\right)\right)}^2\right]$,                   $c_1=\cos \left(\frac{\pi }{8}\right),c_2=\sin \left(\frac{\pi }{8}\right)$	0	$2e$
maxq [22]	$\underset{i\in \{1,\dots ,n\}}{\max }\left\{x_i^2\right\}$	0	0
mxhilb [28]	$\underset{i\in \{1,\dots ,n\}}{\max }\left\{\sum _{j=1}^n\left|\frac{x_j}{i+j-1}\right|\right\}$	0	0
quad	$\sum _{i=1}^n{\left(x_i-2\right)}^2$	0	$2e$

Table 4. Table B1: Characteristics of the first 12 instances of ( CPF ) generated by Algorithm  1 minimizing the convex quadratic function abhi on Ω = [ − 2 , 2 ] 3 ∩ ℤ 3 Ω superscript 2 2 3 superscript ℤ 3 \Omega=[-2,2]^{3}\cap\mathbb{Z}^{3} . ( CPF ) instances are generated by AMPL and solved by CPLEX; times are the mean of five replications.

k	sHyp	LB	UB	time	simIter	nodes	bVars	cVars	cons	$\hat{x}$
1	20	-616.3	79.9	0.1	374	0	335	268	960	$[2;2;-2]$
2	52	-555.1	79.9	3.9	13,180	8,089	847	685	2,466	$[2;2;-1]$
3	100	-475.2	44.7	10.9	31,423	8,555	1,615	1,310	4,724	$[2;1;-2]$
4	172	-434.4	44.7	7.3	19,267	1,728	2,767	2,247	8,110	$[1;2;-2]$
5	276	-413.9	19.1	30.8	68,264	5,874	4,431	3,600	13,000	$[2;1;-1]$
6	418	-373.1	19.1	95.1	84,031	7,933	6,703	5,447	19,676	$[1;2;-1]$
7	611	-311.7	19.1	59.6	83,102	5,440	9,791	7,957	28,749	$[2;-2;-2]$
8	866	-293.2	19.1	99.6	86,318	3,933	13,871	11,273	40,736	$[1;1;-2]$
9	1,196	-232.0	19.1	154.2	84,440	4,568	19,151	15,564	56,248	$[1;1;-1]$
10	1,532	-199.5	19.1	452.3	235,473	6,400	24,527	19,933	72,042	$[2;-2;-1]$
11	2,038	-192.9	19.1	1,006	387,491	9,686	32,623	26,512	95,826	$[2;-1;-2]$
12	2,605	-140.9	19.1	964.3	455,939	29,279	41,695	33,884	122,477	$[1;-1;-2]$

Table 5. Table C1: Number of function evaluations before solvers terminate for n = 3 𝑛 3 n=3 test problems. Parentheses show number of evaluations until x ∗ superscript 𝑥 x^{*} is evaluated.

	SUCIL-ideal1	SUCIL	DFLINT	DFLINT-M	NOMAD	NOMAD-dm	MATSuMoTo
abhi	19 (8)	30 (17)	161 (57)	150 (10)	59 (20)	129 (56)	60 (31)
quad	25 (8)	39 (17)	157 (54)	150 (10)	53 (18)	119 (35)	45 (16)
KLT	22 (8)	28 (13)	152 (48)	146 (15)	51 (24)	116 (34)	46 (17)
maxq	14 (1)	14 (1)	181 (1)	181 (1)	34 (1)	107 (1)	50 (21)
mxhilb	16 (1)	21 (1)	181 (1)	181 (1)	35 (1)	106 (1)	48 (19)
LQ	18 (8)	36 (7)	182 (7)	181 (7)	48 (7)	107 (7)	41 (12)
CB3I	22 (8)	25 (10)	184 (22)	181 (7)	56 (12)	108 (12)	60 (31)
CB3II	22 (8)	34 (11)	184 (22)	181 (7)	44 (12)	108 (12)	60 (31)

Table 6. Table C2: Number of function evaluations before solvers terminate for n = 4 𝑛 4 n=4 test problems. Parentheses show number of evaluations until x ∗ superscript 𝑥 x^{*} is evaluated.

	SUCIL-ideal1	SUCIL	DFLINT	DFLINT-M	NOMAD	NOMAD-dm	MATSuMoTo
abhi	45 (10)	75 (41)	993 (137)	959 (13)	110 (16)	453 (88)	89 (60)
quad	51 (10)	95 (21)	982 (124)	959 (13)	111 (12)	460 (89)	56 (25)
KLT	51 (10)	67 (21)	954 (141)	953 (20)	116 (42)	457 (55)	54 (23)
maxq	30 (1)	33 (1)	1,001 (1)	1,001 (1)	91 (1)	451 (1)	89 (60)
mxhilb	32 (1)	65 (1)	1,001 (1)	1,001 (1)	90 (1)	450 (1)	63 (34)
LQ	39 (10)	109 (15)	1,000 (17)	1,001 (9)	145 (9)	453 (10)	47 (18)
CB3I	50 (10)	58 (14)	1,001 (58)	1,001 (9)	149 (42)	456 (23)	126 (81)
CB3II	48 (10)	91 (14)	1,001 (50)	1,001 (19)	125 (30)	460 (35)	131 (76)

Table 7. Table C3: Number of function evaluations before solvers terminate for n = 5 𝑛 5 n=5 test problems. Parentheses show number of evaluations until x ∗ superscript 𝑥 x^{*} is evaluated.

	SUCIL-ideal1	SUCIL	DFLINT	DFLINT-M	NOMAD	NOMAD-dm	MATSuMoTo
abhi	105 (12)	154 (113)	1,001 (167)	1,000 (16)	301 (41)	1,000 (156)	214 (157)
quad	108 (12)	146 (58)	1,001 (186)	1,000 (16)	286 (26)	1,000 (58)	113 (62)
KLT	108 (12)	121 (69)	1,001 (193)	1,001 (25)	296 (43)	1,000 (176)	108 (77)
maxq	75 (1)	80 (1)	1,001 (1)	1,001 (1)	257 (1)	1,000 (1)	284 (255)
mxhilb	96 (1)	154 (1)	1,001 (1)	1,001 (1)	257 (1)	1,000 (1)	131 (102)
LQ	83 (12)	126 (17)	1,001 (44)	1,001 (11)	425 (15)	1,000 (15)	56 (27)
CB3I	114 (12)	155 (68)	1,001 (103)	1,001 (11)	417 (18)	1,000 (18)	266 (237)
CB3II	100 (12)	135 (66)	1,001 (130)	1,001 (26)	465 (69)	1,000 (54)	281 (224)

Table 8. Table C4: Number of function evaluations taken by 20 replications of MATSuMoTo for each of the 8 8 8 convex test problems for n = 3 𝑛 3 n=3 .

	abhi	quad	KLT	maxq	mxhilb	LQ	CB3I	CB3II
1	48 (19)	48 (19)	48 (19)	43 (14)	48 (19)	43 (14)	59 (30)	83 (54)
2	63 (34)	44 (15)	44 (15)	48 (19)	44 (15)	43 (14)	43 (14)	54 (25)
3	44 (15)	43 (14)	45 (16)	48 (19)	49 (20)	43 (14)	48 (19)	106(77)
4	58 (29)	43 (14)	48 (19)	53 (24)	64 (35)	39 (10)	48 (19)	63 (34)
5	79 (50)	48 (19)	49 (20)	54 (25)	44 (15)	38 (9)	53 (24)	53 (24)
6	73 (44)	37 (9)	48 (19)	39 (10)	53 (24)	39 (10)	64 (35)	53 (24)
7	88 (59)	43 (14)	44 (15)	58 (29)	49 (20)	44 (15)	54 (25)	59 (30)
8	60 (31)	48 (19)	43 (14)	38 (9)	43 (14)	44 (15)	73 (44)	73 (44)
9	79 (50)	43 (14)	48 (19)	53 (24)	43 (14)	39 (10)	68 (39)	58 (29)
10	65 (36)	53 (24)	49 (20)	53 (24)	49 (20)	44 (15)	49 (20)	43 (14)
11	63 (34)	48 (19)	48 (19)	48 (19)	58 (29)	43 (14)	63 (34)	48 (19)
12	53 (24)	44 (15)	43 (14)	48 (19)	48 (19)	38 (9)	58 (29)	53 (24)
13	48 (19)	49 (20)	48 (19)	54 (25)	49 (20)	53 (24)	68 (39)	37 (9)
14	43 (14)	48 (19)	48 (19)	43 (14)	49 (20)	43 (14)	78 (49)	49 (20)
15	44 (15)	48 (19)	48 (19)	48 (19)	48 (19)	43 (14)	39 (10)	73 (44)
16	58 (29)	48 (19)	43 (14)	58 (29)	37 (9)	38 (9)	58 (29)	73 (44)
17	48 (19)	44 (15)	43 (14)	53 (24)	49 (20)	37 (9)	68 (39)	48 (19)
18	64 (35)	43 (14)	48 (19)	63 (34)	53 (24)	38 (9)	93 (64)	69 (40)
19	74 (45)	43 (14)	48 (19)	58 (29)	53 (24)	37 (9)	68 (39)	63 (34)
20	49 (20)	43 (14)	43 (14)	58 (29)	48 (19)	44 (15)	54 (25)	58 (29)
$\lfloor \mathrm{mean}\rfloor$	60 (31)	45 (16)	46 (17)	50 (21)	48 (19)	41 (12)	60 (31)	60 (31)

Table 9. Table C5: Number of function evaluations taken by 20 replications of MATSuMoTo for each of the 8 8 8 convex test problems for n = 4 𝑛 4 n=4 .

	abhi	quad	KLT	maxq	mxhilb	LQ	CB3I	CB3II
1	110 (81)	56 (27)	51 (22)	66 (37)	60 (31)	41 (12)	230 (20)	85 (56)
2	70 (41)	55 (26)	60 (31)	120 (91)	60 (31)	46 (17)	76 (47)	85 (56)
3	50 (21)	50 (21)	50 (21)	94 (65)	50 (21)	41 (12)	150 (121)	196 (11)
4	104 (75)	50 (21)	50 (21)	65 (36)	80 (51)	45 (16)	110 (81)	186 (16)
5	91 (62)	60 (11)	55 (26)	112 (83)	80 (51)	41 (12)	100 (11)	161 (132)
6	66 (37)	55 (26)	56 (27)	90 (61)	56 (27)	55 (26)	101 (72)	196 (167)
7	86 (57)	65 (36)	45 (16)	65 (36)	50 (21)	46 (17)	90 (61)	141 (112)
8	90 (61)	55 (26)	65 (11)	81 (52)	65 (36)	41 (12)	105 (63)	346 (161)
9	166 (137)	56 (27)	56 (27)	55 (26)	50 (21)	51 (22)	96 (67)	100 (71)
10	85 (56)	65 (21)	55 (11)	60 (31)	65 (36)	55 (26)	55 (26)	106 (77)
11	147 (118)	70 (41)	46 (17)	70 (41)	65 (36)	46 (17)	95 (66)	80 (51)
12	96 (67)	45 (16)	46 (17)	112 (83)	71 (42)	50 (21)	112 (83)	231 (202)
13	39 (11)	60 (31)	71 (42)	119 (90)	71 (42)	41 (12)	50 (21)	81 (52)
14	85 (56)	61 (32)	50 (21)	85 (56)	55 (26)	55 (26)	205 (176)	55 (26)
15	101 (72)	60 (31)	60 (31)	94 (65)	55 (26)	41 (12)	160 (131)	114 (85)
16	66 (37)	51 (22)	45 (16)	50 (21)	75 (46)	41 (12)	65 (36)	65 (36)
17	76 (47)	65 (36)	56 (27)	55 (26)	60 (31)	60 (31)	85 (56)	131 (102)
18	70 (41)	56 (27)	51 (22)	110 (81)	75 (46)	65 (36)	314 (285)	70 (41)
19	81 (52)	56 (11)	60 (31)	221 (192)	71 (42)	41 (12)	117 (16)	141 (52)
20	105 (76)	46 (17)	60 (31)	66 (37)	50 (21)	56 (27)	210 (13)	50 (21)
$\lfloor \mathrm{mean}\rfloor$	89 (60)	56 (25)	54 (23)	89 (60)	63 (34)	47 (18)	126 (81)	131 (76)

Table 10. Table C6: Number of function evaluations taken by 20 replications of MATSuMoTo for each of the 8 8 8 convex test problems for n = 5 𝑛 5 n=5 .

	abhi	quad	KLT	maxq	mxhilb	LQ	CB3I	CB3II
1	472 (443)	82 (53)	83 (54)	137 (108)	149 (120)	52 (23)	199 (170)	83 (54)
2	358 (329)	83 (54)	317 (288)	229 (200)	208 (179)	63 (34)	267 (238)	518 (489)
3	193 (164)	83 (54)	172 (143)	115 (86)	175 (146)	62 (33)	377 (348)	1,626 (133)
4	194 (165)	73 (44)	62 (33)	171 (142)	72 (43)	57 (28)	223 (194)	325 (296)
5	160 (131)	72 (43)	73 (44)	504 (475)	87 (58)	53 (24)	835 (806)	487 (458)
6	396 (367)	187 (127)	78 (49)	92 (63)	152 (123)	53 (24)	127 (98)	196 (18)
7	82 (53)	82 (53)	120 (91)	135 (106)	247 (218)	62 (33)	184 (155)	603 (574)
8	62 (33)	106 (15)	58 (29)	969 (940)	186 (157)	52 (23)	164 (135)	67 (38)
9	107 (78)	83 (54)	68 (39)	276 (247)	103 (74)	57 (28)	399 (370)	671 (642)
10	186 (13)	78 (18)	73 (44)	352 (323)	82 (53)	73 (44)	364 (335)	205 (18)
11	181 (152)	130 (65)	62 (33)	255 (226)	53 (24)	58 (29)	414 (385)	240 (211)
12	239 (210)	115 (86)	125 (96)	403 (374)	181 (152)	47 (18)	83 (54)	462 (433)
13	295 (211)	87 (18)	72 (43)	167 (138)	151 (122)	48 (19)	128 (99)	324 (33)
14	354 (325)	78 (49)	72 (43)	236 (207)	82 (53)	62 (33)	200 (171)	154 (125)
15	113 (13)	73 (44)	52 (14)	512 (483)	88 (59)	53 (24)	271 (242)	87 (58)
16	227 (23)	82 (53)	120 (91)	189 (160)	93 (64)	68 (39)	326 (297)	63 (34)
17	102 (73)	135 (106)	270 (241)	195 (166)	57 (28)	58 (29)	180 (151)	599 (570)
18	139 (110)	268 (80)	150 (121)	173 (144)	112 (83)	58 (29)	228 (199)	112 (83)
19	169 (19)	67 (38)	72 (43)	108 (79)	117 (88)	42 (13)	204 (175)	93 (64)
20	268 (239)	301 (195)	62 (13)	465 (436)	237 (208)	53 (24)	162 (133)	180 (151)
$\lfloor \mathrm{mean}\rfloor$	214 (157)	113 (62)	108 (77)	284 (255)	131 (102)	56 (27)	266 (237)	281 (224)

Equations

$${\text{mbox}minimize}_{x}f(x)\text{mbox}subjecttox\in \Omega \subset {\mathbb{Z}}^n.$$

$${\text{mbox}minimize}_{x,y}F(x,y)\text{mbox}subjectto(x,y)\in \Omega \times \Psi \subset {\mathbb{Z}}^n\times {\mathbb{R}}^m$$

$$f(x)=\underset{y\in \Psi }{\min }F(x,y)$$

$$m^{\mbox{\boldmath$i$}}(x)=(c^{\mbox{\boldmath$i$}})^{\operatorname{T}}x+b^{\mbox{\boldmath$i$}},$$

$$\left[\bar{X}^{\mbox{\boldmath$i$}}\;e\right]\left[\begin{array}[]{l}c^{\mbox{\boldmath$i$}}\\ b^{\mbox{\boldmath$i$}}\end{array}\right]=f^{\mbox{\boldmath$i$}},\;\text{where}\;\bar{X}^{\mbox{\boldmath$i$}}=\left[\begin{array}[]{c}(x^{i_{1}})^{\operatorname{T}}\\ \vdots\\ (x^{i_{n+1}})^{\operatorname{T}}\end{array}\right],\;e=\left[\begin{array}[]{c}1\\ \vdots\\ 1\end{array}\right],\;\mbox{and}\;f^{\mbox{\boldmath$i$}}=\left[\begin{array}[]{c}f(x^{i_{1}})\\ \vdots\\ f(x^{i_{n+1}})\end{array}\right].$$

$$\mathcal{U}^{\mbox{\boldmath$i$}}=\bigcup\limits_{i_{j}\in\mbox{\boldmath$i$}}\mathrm{cone}\big{(}x^{i_{j}}-X^{\mbox{\boldmath$i$}}\big{)},$$

$$\mathrm{cone}\big{(}x^{i_{j}}-X^{\mbox{\boldmath$i$}}\big{)}=\{x^{i_{j}}+\sum_{l=1,l\neq j}^{n+1}\lambda_{l}(x^{i_{j}}-x^{i_{l}}):\,i_{j}\in\mbox{\boldmath$i$},i_{l}\in\mbox{\boldmath$i$},\lambda_{l}\geq 0\}.$$

$$x=x^{i_j}+\sum _{l=1,l\ne j}^{n+1}{\lambda }_l\left(x^{i_j}-x^{i_l}\right),$$

$$x^{i_j}=\frac{1}{1+{\sum }_{k=1,k\ne j}^{n+1}{\lambda }_k}x+\frac{1}{1+{\sum }_{k=1,k\ne j}^{n+1}{\lambda }_k}\sum _{l=1,l\ne j}^{n+1}{\lambda }_l{x}^{i_l},$$

$$f(x^{i_j})\le \frac{1}{1+{\sum }_{k=1,k\ne j}^{n+1}{\lambda }_k}f(x)+\frac{1}{1+{\sum }_{k=1,k\ne j}^{n+1}{\lambda }_k}\sum _{l=1,l\ne j}^{n+1}{\lambda }_{l}f\left(x^{i_l}\right).$$

$$f(x)$$

$$\displaystyle=\left(1+\sum_{k=1,k\neq j}^{n+1}\lambda_{k}\right)m^{\mbox{\boldmath$i$}}(x^{i_{j}})-\sum_{l=1,l\neq j}^{n+1}\lambda_{l}m^{\mbox{\boldmath$i$}}(x^{i_{l}})$$

$$\displaystyle=m^{\mbox{\boldmath$i$}}(x^{i_{j}})+\sum_{l=1,l\neq j}^{n+1}\lambda_{l}\left(m^{\mbox{\boldmath$i$}}(x^{i_{j}})-m^{\mbox{\boldmath$i$}}(x^{i_{l}})\right)$$

$$\displaystyle=m^{\mbox{\boldmath$i$}}(x),$$

$$0$$

$$=x^{i_1}-x^{i_2}+\sum _{l=2}^{n+1}{\lambda }_l{x}^{i_1}-\sum _{l=2}^{n+1}{\lambda }_l{x}^{i_2}+\sum _{l=2}^{n+1}{\lambda }_l{x}^{i_2}-\sum _{l=2}^{n+1}{\lambda }_l{x}^{i_l}-\sum _{l=1,l\ne 2}^{n+1}{\sigma }_l(x^{i_2}-x^{i_l})$$

$$=\left(1+\sum _{l=2}^{n+1}{\lambda }_l\right)\left(x^{i_1}-x^{i_2}\right)-\sum _{l=2}^{n+1}{\lambda }_l\left(x^{i_l}-x^{i_2}\right)+\sum _{l=1,l\ne 2}^{n+1}{\sigma }_l\left(x^{i_l}-x^{i_2}\right)$$

$$=\left(1+{\sigma }_1+\sum _{l=2}^{n+1}{\lambda }_l\right)\left(x^{i_1}-x^{i_2}\right)+\sum _{l=3}^{n+1}\left({\sigma }_l-{\lambda }_l\right)\left(x^{i_l}-x^{i_2}\right).$$

$$1+{\sigma }_1+\sum _{l=2}^{n+1}{\lambda }_l=0,$$

$$W(X)=\left\{\mbox{\boldmath$i$}:X^{\mbox{\boldmath$i$}}\subseteq X,X^{\mbox{\boldmath$i$}}\mbox{ is poised}\right\}.$$

$${\text{mbox}minimize}_{x,\eta }$$

$$\displaystyle\eta\geq(c^{\mbox{\boldmath$i$}})^{\operatorname{T}}x+b^{\mbox{\boldmath$i$}},\text{ for all }\mbox{\boldmath$i$}\in W(X)\text{ with }x\in\mathcal{U}^{\mbox{\boldmath$i$}}$$

$$x\in \Omega ,$$

$$(c^{\mbox{\boldmath$i$}})^{\operatorname{T}}x+b^{\mbox{\boldmath$i$}}\leq f(x),$$

$$x^0=\sum _{j=1}^{n+1}{\alpha }_j{x}^j,\text{mbox}where\sum _{j=1}^{n+1}{\alpha }_j=1\text{mbox}and{\alpha }_j>0\text{mbox}forj=1,\dots ,n+1.$$

$$\sum _{j=0}^{n-1}{\beta }_j(x^j-x^n)=0.$$

$$0$$

$$\displaystyle=\sum_{j=1}^{n-1}(\beta_{0}\alpha_{j}+\beta_{j})x^{j}+\bigg{(}\beta_{0}(\alpha_{n}-1)-\sum_{j=1}^{n-1}\beta_{j}\bigg{)}x^{n}+\beta_{0}\alpha_{n+1}x^{n+1}.$$

$$\sum _{j=1}^{n-1}({\beta }_0{\alpha }_j+{\beta }_j)+{\beta }_0({\alpha }_n-1)-\sum _{j=1}^{n-1}{\beta }_j+{\beta }_0{\alpha }_{n+1}={\beta }_0\sum _{j=1}^{n+1}{\alpha }_j-{\beta }_0=0,$$

$$\bigcup\limits_{\mbox{\boldmath$i$}\in W(X\cup\{x^{0}\})}\mathcal{U}^{\mbox{\boldmath$i$}}=\mathbb{R}^{n}.$$

$$0=(\sum _{j=1}^{n+1}{\alpha }_j)(x^0-x^0)=(\sum _{j=1}^{n+1}{\alpha }_j)x^0-\sum _{j=1}^{n+1}{\alpha }_j{x}^j=\sum _{j=1}^{n+1}{\alpha }_j(x^0-x^j),$$

$$x=\sum _{j=1}^n{\alpha }_j(x^0-x^j),$$

$$x-x^0=\sum _{j=1,j\ne l}^{n+1}{\lambda }_j(x^0-x^j),$$

$${\lambda }_{j^{\prime }}(x^0-x^{j^{\prime }})=\frac{-{\lambda }_{j^{\prime }}}{{\alpha }_{j^{\prime }}}\left(\sum _{j=1,j\ne j^{\prime }}^{n+1}{\alpha }_j(x^0-x^j)\right).$$

$$x-x^0=\sum _{j=1,j\ne j^{\prime }}^{n+1}{\mu }_j(x^0-x^j),$$

Full text

1 Introduction

We study the problem of minimizing a convex function on a finite subset of the integer lattice. In particular, we consider problems of the form

[TABLE]

We first define what it means for $f$ to be convex on a set $\Omega\subset\mathbb{Z}^{n}$.

Definition 1** (Convexity on Integer Subsets).**

A function $f$ is convex on $\Omega$ if for $x\in\Omega$ and any $p$ points $y^{i}\in\Omega$ satisfying $x=\sum_{i=1}^{p}\lambda_{i}y^{i}$ with $\lambda_{i}\in[0,1]$ and $\sum_{i=1}^{p}\lambda_{i}=1$, then $f(x)\leq\sum_{i=1}^{p}\lambda_{i}f(y^{i})$.

We make the following assumption about problem (1).

Assumption 1**.**

$f$* is convex on $\Omega$, $\Omega\subset\mathbb{Z}^{n}$ is nonempty and bounded, and $f$ cannot be evaluated at $x\notin\Omega$.*

Because we assume that $f$ cannot be evaluated at noninteger points, problem (1) can be referred to as a convex optimization problem with unrelaxable integer constraints [30]. We note that $\Omega$ need not contain all integer points in its convex hull (i.e., our approach allows for situations where $\mathrm{conv}\left(\Omega\right)\cap\mathbb{Z}^{n}\not=\Omega$). Problem (1) can also be viewed as minimizing an integer-convex111Although [24] considers integer convexity only for polynomials, this definition can be applied to more general classes of functions over the sets considered here. function [24, Definition 15.2] over a nonempty finite subset of $\mathbb{Z}^{n}$.

Admittedly, it is rare to know that $f$ is convex when $f$ is not given in closed form (although one may be able to detect convexity [27] or estimate a probability that $f$ is convex on the finite domain $\Omega$ [26]). Nevertheless, studying the convex case is important because we are unaware of any method (besides complete enumeration) for obtaining exact solutions to (1) when $f$ (convex or otherwise) cannot be evaluated at noninteger points.

One example where an objective is not given in closed form but is known to be convex arises in the combinatorial optimal control of partial differential equations (PDEs). For example, Buchheim et al. [9, Lemma 2] show that the solution operator of certain semilinear elliptic PDEs is a convex function of the controls provided that the nonlinearities in the PDE and boundary conditions are concave and nondecreasing. Thus, any linear function of the states of the PDE (e.g., the $\max$-function) is a convex function of the controls when the states are eliminated. The authors of [9] propose using adjoint information to compute subgradients of the (continuous relaxation of the) objective, but an alternative would be to consider a derivative-free approach.

We consider only pure-integer problems of the form (1); however, our developments are equally applicable to the mixed-integer case

[TABLE]

provided $F$ is convex on $\Omega\times\Psi$. If we define the function

[TABLE]

and if $f$ is well defined over $\Omega$, then (2) can be solved by minimizing $f$ on $\Omega\subset\mathbb{Z}^{n}$, where each evaluation $f(x)$ requires an optimization of the continuous variables $y$ for a fixed $x$. (When there is no additional information on $F$ and $\Psi$, each continuous optimization problem may be difficult to solve.) Because many of the results below rely only on the convexity of $f$ and not the discrete nature of $\Omega$, much of the analysis below readily applies to the mixed-integer case.

We are especially interested in problems where the cost of evaluating $f$ is large. Problems of the form (1) or (2) where the objective is expensive to evaluate and some integer constraints are unrelaxable arise in a range of simulation-based optimization problems. For example, the optimal design of concentrating solar power plants gives rise to computationally expensive simulations for each set of design parameters [41]. Furthermore, some of the design parameters (e.g., the number of panels on the power plant receiver) cannot be relaxed to noninteger values. Similar problems arise when tuning codes to run on high-performance computers [6]. In this case, $f(x)$ may be the memory footprint of a code that is compiled with settings $x$, which can correspond to decisions such as loop unrolling or tiling factors that do not have meaningful noninteger values. Optimal material design problems may also constrain the choice of atoms to a finite set, resulting in unrelaxable integer constraints; see [21] for a derivative-free optimization algorithm designed explicitly for such a problem.

Motivated by such applications, we develop a method that will certifiably converge to the solution of (1) under Assumption 1 without access to $\partial f$. Using only evaluations of $f$, we construct secants, which are linear functions that interpolate $f$ at a set of $n+1$ points. These secant functions underestimate $f$ in certain parts of $\Omega$. We use these secants to define conditional cuts that are valid in disconnected portions of the domain. The complete set of secants and the conditions that describe when they are valid are used to construct an underestimator of a convex $f$. While access to $\partial f$ (i.e., a (sub)gradient of a continuous relaxation of $f$) could strengthen such an underestimator, we do not address such considerations in this paper.

Solving (1) under Assumption 1 without access to $\partial f$ poses a number of theoretical and computational challenges. Because the integer constraints are unrelaxable, one cannot apply traditional branch-and-bound approaches. In particular, model-based continuous derivative-free methods would require evaluating the objective at noninteger points to ensure convergence for the continuous relaxation of (1). In addition, other traditional techniques for mixed-integer optimization—such as Benders decomposition [20] or outer approximation [17, 18]—cannot be used to solve (1) when $\partial f$ is unavailable. Since we know of no method (other than complete enumeration) for obtaining global minimizers of (1) under Assumption 1, we know of no potential algorithm to address this problem when a (sub)gradient is unavailable.

We make three contributions in this paper: (1) we develop a new underestimator for convex functions on subsets of the integer lattice that is based solely on function evaluations; (2) we present an algorithm that alternates between updating this underestimator and evaluating the objective in order to identify a global solution of problem (1) under Assumption 1; and (3) we show empirically that certifying global optimality in such cases is a challenging problem. In our experiments, we are unable to prove optimality for many problems when $n\geq 5$, and we provide insights into why a proof of optimality remains computationally challenging.

Outline.

Section 2 surveys recent methods for addressing (1). Section 3 introduces valid conditional cuts using only the function values of a convex objective and discusses the theoretical properties of these cuts. Section 4 presents an algorithm for solving (1) and shows that this algorithm identifies a global minimizer of (1) under Assumption 1. Section 5 considers two approaches for formulating the underestimator and presents the method SUCIL—secant underestimator of convex functions on the integer lattice. Section 6 provides detailed numerical studies for implementations of SUCIL on a set of convex problems. Section 7 discusses many of the challenges in obtaining global solutions to (1).

2 Background

Developing methods to solve (1) without access to derivatives of $f$ is an active area of research. Most methods address general (i.e., nonconvex) functions $f$, and heuristic approaches are commonly adopted to handle integer decision variables for such derivative-free optimization problems. For example, the method in [37] rounds noninteger components of candidate points to the nearest feasible integer values. The method’s asymptotic convergence results are based on the inclusion of points drawn uniformly from the finite domain (and rounding noninteger values as necessary).

Integer-constrained pattern-search methods [2, 4] generalize their continuous counterparts. These modified pattern-search methods can be shown to converge to mesh-isolated minimizers: points with function values that are better than all neighboring points on the integer lattice. Unfortunately, such mesh-isolated minimizers can be arbitrarily far from a global minimizer, even when $f$ is convex; see [1, Fig. 2] for an example of such a function $f$. Other methods that converge to mesh-isolated minimizers include direct-search methods that update the integer variables via a local search [31, 32] and mesh adaptive direct-search methods adapted to address discrete and granular variables (i.e., those that have a controlled number of decimals) [3, 5]. The direct-search method in [19] accounts for integer constraints by constructing a set of directions that have a nonnegative span of $\mathbb{R}^{n}$ and that ensure that all evaluated points will be integer valued. This method is shown to converge to a type of stationary point that, even in the convex case, may not be a global minimizer. See [38] for various definitions of local minimizers of (1) and a discussion of associated properties. The BFO method in [39] has a recursive step that explores points near the current iterate by fixing each of the discrete variables to its value plus or minus a step-size parameter.

The method in [33] uses line searches over a set of primitive directions, that is, a set of scaled directions $D$ where no vector $d_{i}\in D$ is a positive multiple of a different $d_{j}\in D$. This method explores a discrete set of directions around the current iterate until finding a local minimum $x_{c}$ in a $\beta$-neighborhood, defined as $\mathcal{N}(x_{c},\beta)=\{x_{c}+\alpha d\in\Omega:d\in D,\alpha\in\mathbb{N},\alpha\leq\beta\}$ for $\beta\in\mathbb{N}$. Although the authors of [33] target nonconvex objectives, their approach will converge to a global minimum $x_{c}$ of a convex objective $f$ if all points in $\mathcal{N}(x_{c},1)$ are evaluated. Figure 1 illustrates such a discrete 1-neighborhood. Unfortunately, $|\mathcal{N}(x_{c},1)|$ can be large; see Table 1.

Model-based methods approximate objective functions on the integer lattice by using surrogate models; see, for example, [25], [40], and [13]. The surrogate model is used to determine points where the objective should be evaluated; the model is typically refined after each objective evaluation. The methodology in [13] specifically uses radial basis function surrogate models and does automatic model selection at each iteration. Mixed-integer nonlinear optimization solvers are used to minimize the surrogate to obtain the next integer point for evaluation. The model-based methods of [34, 35, 36] modify the sampling strategies and local searches typically used to solve continuous objective versions. The approaches in [34, 35] restart when a suitably defined local minimizer is encountered, continuing to evaluate the objective until the available budget of function evaluations is exhausted. These model-based methods differ in the initial sampling method, the type of surrogate model, and the sampling strategy used to select the next points to be evaluated. See [7], for a survey and taxonomy of continuous and integer model-based optimization approaches.

In a different line of research, Davis and Ierapetritou [14] propose a branch-and-bound framework to address binary variables; a solution to the relaxed nonlinear subproblems is obtained via a combination of global kriging models and local surrogate models. Similarly, Hemker et al. [23] replace the black-box portions of the objective function (and constraints) by a stochastic surrogate; the resulting mixed-integer nonlinear programs are solved by branch and bound. Both approaches assume that the integer constraints are relaxable.

3 Underestimator of Convex Functions on the Integer Lattice

To construct an underestimator of a convex objective function $f$, we now discuss secant functions, which are linear mappings that interpolate $f$ at $n+1$ points. We provide conditions for where these cuts will underestimate $f$. We then discuss necessary conditions on the set of evaluated points so that if all possible secants are constructed, these conditional cuts underestimate $f$ on the domain $\Omega$. As we will see in Section 4, this underestimator is essential for obtaining a global minimizer of (1) under Assumption 1. Throughout this section, $X\subseteq\Omega$ denotes a set of at least $n+1$ points in $\Omega$ at which the objective function $f$ has been evaluated.

3.1 Secant Functions and Conditional Cuts

Constructing a secant function requires a set of $n+1$ interpolation points where $f$ has been evaluated. To define a secant function for $f$, we introduce a multi-index $i$ of $n+1$ distinct indices, $1\leq i_{1}<\ldots<i_{n+1}\leq\left|X\right|$, as \mbox{\boldmathi}=(i_{1},\ldots,i_{n+1}). With a slight abuse of notation, we will refer to elements i_{j}\in\mbox{\boldmathi}.

Given the set of points X^{\mbox{\boldmathi}}=\left\{x^{i_{j}}\in X\subseteq\Omega:i_{j}\in\mbox{\boldmathi}\right\}, we construct the secant function

[TABLE]

where the coefficients c^{\mbox{\boldmathi}}\in\mathbb{R}^{n} and b^{\mbox{\boldmathi}}\in\mathbb{R} are the solution to the linear system

[TABLE]

The secant function m^{\mbox{\boldmathi}} is unique provided that the set X^{\mbox{\boldmathi}} is poised, which we now define.

Definition 2**.**

The set of points X^{\mbox{\boldmathi}} is poised if the matrix \left[\bar{X}^{\mbox{\boldmathi}}\;e\right] is nonsingular.

Note that Definition 2 is equivalent to X^{\mbox{\boldmathi}} being affinely independent. We now show that the secant function m^{\mbox{\boldmathi}} underestimates $f$ in certain polyhedral cones, namely, the cones

[TABLE]

where

[TABLE]

Lemma 1** (Conditional Cuts).**

If $f$ is convex on $\Omega$ and X^{\mbox{\boldmathi}} $\subseteq\Omega$ is poised, then the unique linear mapping m^{\mbox{\boldmathi}} satisfying m^{\mbox{\boldmathi}}(x^{i_{j}})=f(x^{i_{j}}) for each i_{j}\in\mbox{\boldmathi} satisfies m^{\mbox{\boldmathi}}(x)\leq f(x) for all x\in\mathcal{U}^{\mbox{\boldmathi}}\bigcap\Omega.

Proof.

The uniqueness of the linear mapping follows directly from the affine independence guaranteed by Definition 2 for poised X^{\mbox{\boldmathi}}.

Let $x$ be a point in \mathrm{cone}\big{(}x^{i_{j}}-X^{\mbox{\boldmathi}}\big{)}\bigcap\Omega for arbitrary x^{i_{j}}\in X^{\mbox{\boldmathi}}. By (5),

[TABLE]

with $\lambda_{l}\geq 0$ (for $l=1,\ldots,n+1;l\neq j$). Rearranging (6) yields

[TABLE]

showing that $x^{i_{j}}$ can be expressed as a convex combination of $\{x\}\cup\{x^{i_{l}}:l=1,\ldots,n+1;l\neq j\}$, all of which are points in $\Omega$. Therefore, by convexity of $f$ on $\Omega$ (see Definition 1),

[TABLE]

Solving for $f(x)$ and using the fact that m^{\mbox{\boldmathi}} interpolates $f$ at points in X^{\mbox{\boldmathi}}, we obtain

[TABLE]

where the last equality holds by (6) and the linearity of m^{\mbox{\boldmathi}}. Because $x$ is an arbitrary point in \mathrm{cone}\big{(}x^{i_{j}}-X^{\mbox{\boldmathi}}\big{)} for arbitrary $x^{i_{j}}$, the result is shown. ∎

We now prove that the cones in \mathcal{U}^{\mbox{\boldmathi}} do not intersect when X^{\mbox{\boldmathi}} is poised. We note that the following result holds for any poised set X^{\mbox{\boldmathi}}\subset\mathbb{R}^{n}.

Lemma 2** (A Point Is In At Most One Cone).**

If X^{\mbox{\boldmathi}} is poised, no point $x\in\mathbb{R}^{n}$ satisfies x\in\mathrm{cone}\big{(}x^{i_{j}}-X^{\mbox{\boldmathi}}\big{)} and x\in\mathrm{cone}\big{(}x^{i_{k}}-X^{\mbox{\boldmathi}}\big{)} for x^{i_{j}},x^{i_{k}}\in X^{\mbox{\boldmathi}} and $x^{i_{j}}\neq x^{i_{k}}$.

Proof.

Let $x^{i_{1}}$ and $x^{i_{2}}$ be different, but otherwise arbitrary, points in X^{\mbox{\boldmathi}}. To arrive at a contradiction, suppose that there exists x\in\mathrm{cone}\big{(}x^{i_{1}}-X^{\mbox{\boldmathi}}\big{)}\cap\mathrm{cone}\big{(}x^{i_{2}}-X^{\mbox{\boldmathi}}\big{)}. That is, $x=x^{i_{1}}+\sum_{l=2}^{n+1}\lambda_{l}\left(x^{i_{1}}-x^{i_{l}}\right)$ and $x=x^{i_{2}}+\sum_{l=1,l\neq 2}^{n+1}\sigma_{l}\left(x^{i_{2}}-x^{i_{l}}\right)$ for $\lambda_{l}\geq 0$ ($l\in\left\{2,\ldots,n+1\right\}$) and $\sigma_{l}\geq 0$ ($l\in\left\{1,3,\ldots,n+1\right\}$). Subtracting these two expressions yields

[TABLE]

Since X^{\mbox{\boldmathi}} is a poised set, Definition 2 ensures that the vectors \left\{x^{i_{l}}-x^{i_{2}}:\,i_{l}\in\mbox{\boldmathi},i_{l}\neq i_{2}\right\} are linearly independent. Hence the dependence relation in (3.1) can be satisfied only if the coefficient on $(x^{i_{1}}-x^{i_{2}})$ vanishes. That is,

[TABLE]

which contradicts $\lambda_{l}\geq 0$ (for $l=2,\ldots,n+1$) and $\sigma_{1}\geq 0$. Since $x^{i_{1}}$ and $x^{i_{2}}$ were arbitrary points in X^{\mbox{\boldmathi}}, the result is shown. ∎

For each poised set X^{\mbox{\boldmathi}}, Lemma 1 ensures that the secant function m^{\mbox{\boldmathi}} underestimates $f$ on \Omega\cap\mathcal{U}^{\mbox{\boldmathi}}. We can therefore underestimate $f$ via a model that consists of the pointwise maximum of the underestimators for which the point is in \mathcal{U}^{\mbox{\boldmathi}} for some poised set X^{\mbox{\boldmathi}} of previously evaluated points. Minimizing this nonconvex model on the integer $\Omega$ can then provide a lower bound for the global minimum of $f$ on $\Omega$.

Figure 1 shows a two-dimensional example of points that produce such secant functions and the regions in which they will underestimate any convex function $f$. For the $n+2$ points (circles), we consider three poised sets indicated by triangles linking $n+1$ points. The left image shows a poised set (blue line triangle), and the three cones (shaded blue area) in which the secant through these points is a valid underestimator. The right image shows that conditional cuts using $n+2$ points can cover all of $\mathbb{R}^{n}$.

3.2 Lower Bound for $f$

We now describe an optimization problem whose solution provides a lower bound for $f$ on $\Omega$. Let $W(X)$ denote the set of all multi-indices corresponding to poised subsets of $X$:

[TABLE]

If $f$ has been evaluated at every point in $X$, we can construct a secant function (c^{\mbox{\boldmathi}})^{\operatorname{T}}x+b^{\mbox{\boldmathi}} interpolating $f$ on X^{\mbox{\boldmathi}} for every multi-index \mbox{\boldmathi}\in W(X). We then collect all such conditional cuts in the piecewise linear program

[TABLE]

where \mathcal{U}^{\mbox{\boldmathi}} is defined in (4). For the set of points $X$ and corresponding $W(X)$, let $\eta(\bar{x})$ denote the value of (PLP) when the constraint $x=\bar{x}$ is added to (PLP) for a particular $\bar{x}\in\Omega$. As we will see below, $\eta$ represents the largest lower bound for $f$ induced by the set $X$, and the solution to (PLP) provides a lower bound for $f$ on $\Omega$.

Lemma 3** (Underestimator of $f$).**

If $f$ is convex on $\Omega$, then the optimal value $\eta_{*}$ of (PLP) satisfies $\eta_{*}\leq f(x)$ for all $x\in\Omega$.

Proof.

If $W(X)$ is empty, the result holds trivially since $\eta$ is unconstrained. Otherwise, since (PLP) minimizes $\eta$, it suffices to show that $\eta(x)\leq f(x)$ for arbitrary $x\in\Omega$. Two cases can occur. First, if x\notin\mathcal{U}^{\mbox{\boldmathi}} for every \mbox{\boldmathi}\in W(X), then no conditional cut exists at $x$. Thus $\eta(x)=-\infty$ and $\eta(x)<f(x)$. Second, if x\in\mathcal{U}^{\mbox{\boldmathi}} for some \mbox{\boldmathi}\in W(X), then

[TABLE]

by Lemma 1, where the poisedness of X^{\mbox{\boldmathi}} follows from the definition of $W(X)$. Therefore, $\eta(x)\leq f(x)$ for all $x\in\Omega$. Since $\eta_{*}=\displaystyle\min_{x\in\Omega}\eta(x)$, the result is shown. ∎

If $W(X)$ in (PLP) is replaced by a proper subset $W^{\prime}(X)\subset W(X)$ of multi-indices, then Lemma 3 still holds. (This relaxation of (PLP) associated with removing constraints cannot increase $\eta_{*}$.) Such a replacement may be necessary if $W(X)$ becomes too large to allow considering every poised subset of $X$ when forming (PLP).

3.3 Covering $\mathbb{R}^{n}$ with Conditional Cuts

Because the cuts in (PLP) are valid only within \mathcal{U}^{\mbox{\boldmathi}}, the resulting model takes an optimal value of $\eta_{*}=-\infty$ if there is a point $x\in\Omega$ that is not in the union of \mathcal{U}^{\mbox{\boldmathi}} over all \mbox{\boldmathi}\in W(X). Thus, we find it beneficial to ensure that $X$ contains points that result in a finite objective value for the underestimator described by (PLP). We now establish a condition that ensures that the union of conditional cuts induced by $X$ covers $\mathbb{R}^{n}$ and therefore $\Omega$.

We say that a point $x^{0}$ belongs to the interior of the convex hull of a set of points $X=\{x^{1},\ldots,x^{n+1}\}$ if scalars $\alpha_{j}$ exist such that

[TABLE]

This is denoted by $x^{0}\in\mathrm{int}\left(\mathrm{conv}\left(X\right)\right)$.

Lemma 4** (Poisedness of Initial Points).**

If $X=\{x^{1},\ldots,x^{n+1}\}\subset\mathbb{R}^{n}$ is a poised set and if $x^{0}$ satisfies $x^{0}\in\mathrm{int}\left(\mathrm{conv}\left(X\right)\right)$, then all subsets of $n+1$ points in $\{x^{0}\}\cup X$ are poised.

Proof.

For contradiction, suppose that the set $\{x^{0}\}\cup X\setminus\{x^{n+1}\}$ is not poised and therefore is affinely dependent. Then, there must exist scalars $\beta_{j}$ not all zero, and (without loss of generality) $x^{n}\in X$ such that

[TABLE]

Replacing $x^{0}$ with (9) in the left-hand side above yields

[TABLE]

Since $X$ is poised, the vectors $\{x^{1},\ldots,x^{n+1}\}$ are affinely independent; by definition of affine independence, the only solution to $\sum_{j=1}^{n+1}\gamma_{j}x^{j}=0$ and $\sum_{j=1}^{n+1}\gamma_{j}=0$ is $\gamma_{j}=0$ for $j=1,\ldots,n+1$. The sum of the coefficients from (3.3) satisfies

[TABLE]

because $\sum_{j=1}^{n+1}\alpha_{j}=1$. This means that all coefficients of $x^{j}$ in (3.3) are also equal to zero. Since $\alpha_{n+1}>0$, the last term from (3.3) implies that $\beta_{0}=0$. Considering the remaining coefficients in (3.3), we conclude that $\beta_{0}\alpha_{j}+\beta_{j}=0$, which implies that $\beta_{j}=0$ for $j=1,\ldots,n-1$. This contradicts the assumption that not all $\beta_{j}=0$. Hence, the result is proved. ∎

We now establish a simple set of points that produces conditional cuts that cover $\mathbb{R}^{n}$ and, therefore, the domain $\Omega$.

Lemma 5** (Initial Points and Coverage of $\Omega$).**

Let $X$ be a poised set of $n+1$ points, let $x^{0}\in\mathrm{int}\left(\mathrm{conv}\left(X\right)\right)$, and let $W(X\cup\{x^{0}\})$ be defined as in (8). Then,

[TABLE]

Proof.

Since $x^{0}\in\mathrm{int}\left(\mathrm{conv}\left(X\right)\right)$, there exist $\alpha_{j}>0$ such that

[TABLE]

where the second equality follows from (9). The existence of $\alpha_{j}>0$ such that $\sum_{j=1}^{n+1}\alpha_{j}(x^{0}-x^{j})=0$ implies that the vectors $\left\{x^{0}-x^{j}:j\in\{1,\ldots,n+1\}\right\}$ are a positive spanning set (see, e.g., [12, Theorem 2.3 (iii)]). Therefore any $x\in\mathbb{R}^{n}$ can be expressed as

[TABLE]

with $\alpha_{j}\geq 0$ for all $j$.

We will show that any $x\in\mathbb{R}^{n}$ belongs to \mathcal{U}^{\mbox{\boldmathi}} for some multi-index $i$ containing $x^{0}$. By Lemma 4, every set of $n$ distinct vectors of the form $(x^{0}-x^{j})$ for $x^{j}\in X$ is a linearly independent set. Thus we can express

[TABLE]

for some $l\in\{1,\ldots n+1\}$. If $\lambda_{j}\geq 0$ for each $j$, then we are done, and x\in\mathrm{cone}\big{(}x^{0}-X\setminus\{x^{l}\}\big{)}.

Otherwise, choose an index $j^{\prime}$ such that $\lambda_{j^{\prime}}$ is the most negative coefficient on the right of (13) (breaking ties arbitrarily). Using (12), we can exchange the indices $l$ and $j^{\prime}$ in (13) by observing that

[TABLE]

Note that $\frac{-\lambda_{j^{\prime}}}{\alpha_{j^{\prime}}}\alpha_{j}>0$ by (9), and we can rewrite (13) as

[TABLE]

with new coefficients $\mu_{j}$ that are strictly larger than $\lambda_{j}$:

[TABLE]

Observe that (14) has the same form as (13) but with coefficients $\mu_{j}$ that are strictly greater than $\lambda_{j}$. We can now define $\lambda=\mu$ and repeat the process. If there is some $\lambda_{j^{\prime}}<0$, the process will strictly increase all $\lambda_{j}$. Because there are only a finite number of subsets of size $n$, we must eventually have all $\lambda_{j}\geq 0$. Once $\lambda_{j^{\prime}}$ has been pivoted out, it can reenter only with a positive value (like $\mu_{l}$ above), so eventually all $\lambda_{j}$ will be nonnegative. ∎

Lemma 5 ensures that any poised set of $n+1$ points with an additional point in their interior will produce conditional cuts that cover $\mathbb{R}^{n}$. Figure 1 illustrates this for $n=2$. An alternative set of $n+2$ points is

[TABLE]

where $e_{i}$ is the $i$th unit vector and $e$ is the vector of ones. Larger sets, such as those of the form

[TABLE]

will similarly guarantee coverage of $\mathbb{R}^{n}$.

We note that the results in this section do not rely on $X$ or $\Omega$ being a subset of $\mathbb{Z}^{n}$. Therefore, the results are readily applicable to the case when $f$ has continuous and integer variables.

4 Convergence Analysis

We now present Algorithm 1 to identify global solutions to (1) under Assumption 1. This algorithm constructs a sequence of underestimators of the form (PLP). Section 5.1 shows two approaches for modeling the underestimator and Section 5.2 highlights other details that are important for an efficient implementation of Algorithm 1. For example, the next point evaluated can be a solution of (PLP) but this need not be the case.

Note that (PLP) provides a valid lower bound for $f$ on $\Omega$. If $X\subseteq\Omega$ are points where $f$ has been evaluated, then $\min\left\{f(x):x\in X\right\}$ is an upper bound for the minimum of $f$ on $\Omega$. Algorithm 1 terminates when the upper bound is equal to the lower bound provided by (PLP). We observe that Algorithm 1 produces a nondecreasing sequence of lower bounds provided that conditional cuts are not removed from (PLP); we show in Theorem 1 that this sequence of lower bounds will converge to the global minimum of $f$ on $\Omega$.

Algorithm 1 resembles a traditional outer-approximation approach [8, 17, 18] in that it obtains a sequence of lower bounds of (1) using an underestimator that is updated after each function evaluation. These function evaluations provide a nonincreasing sequence of upper bounds on the objective; when the upper bound equals the lower bound provided by the underestimator, the method can terminate with a certificate of optimality.

Algorithm 1 leaves open a number of important decisions concerning how (PLP) is formulated and solved and how the next iterate is selected. While we will discuss more involved options for addressing these concerns, a simple choice would be to add all new possible cuts and let the next iterate be a minimizer of (PLP). If this minimizer is not chosen, a (possibly difficult) separation problem may have to be solved to obtain a new iterate in Line 1, for example, when $\Omega$ is sparse. Although such choices can result in computational difficulties, these choices are useful for showing the behavior of Algorithm 1, which we do now. In Figure 2 we see three iterations of Algorithm 1 solving the one-dimensional problem

[TABLE]

Black dots indicate interpolation points where $f$ has been previously evaluated, and green dots indicate the solution to (PLP) in each iteration. The solid red lines show the piecewise linear underestimator of the function. We observe that Lemma 1 can be strengthened for one-dimensional problems where conditional cuts underestimate convex $f$ at all points outside the convex hull of the points used to determine the corresponding secant function. (This is not true for $n>1$.)

We now prove that Algorithm 1 identifies a global minimizer of $f$.

Theorem 1** (Convergence of Algorithm 1).**

If Assumption 1 holds, Algorithm 1 terminates at an optimal solution $x^{*}$ of (1) in finitely many iterations.

Proof.

Algorithm 1 will terminate in a finite number of iterations because Assumption 1 ensures that $\Omega$ is finite and Line 1 ensures that $x^{k}$ is not a previously evaluated element of $\Omega$.

For contradiction, assume that Algorithm 1 terminates at iteration $k^{\prime}$ with $f(\hat{x})>f(x^{*})$ for some $x^{*}\displaystyle\in\operatorname*{arg\,min}_{x\in\Omega}f(x)$. It follows from Line 1 that $x^{*}\notin X^{k^{\prime}}$, because $f(x^{*})<f(\hat{x})$. Lemma 3 ensures that the value of each conditional cut at $x^{*}$ is not larger than $f(x^{*})$, which implies that $\eta(x^{*})\leq f(x^{*})$. Thus, the lower bound satisfies

[TABLE]

Since $l_{k^{\prime}}<u_{k^{\prime}}$, Algorithm 1 did not terminate at iteration $k^{\prime}$, giving a contradiction. Therefore, the result is shown. ∎

A special case of Theorem 1 ensures that Algorithm 1 terminates with a global solution of (1) when $x^{k}$ is an optimal solution of (PLP). As in most integer optimization algorithms, the termination condition $l_{k}=u_{k}$ may be met before $X^{k}=\Omega$; that is, Algorithm 1 need not evaluate all the points in $\Omega$. Termination before $X^{k}=\Omega$ occurs when (PLP) is refined in Step 4 of Algorithm 1 and the lower bound at all points (and on the optimal value of (1)) is tightened. When the next iterate is chosen, the upper bound typically improves and convergence occurs faster than enumeration. Yet, in the worst-case scenario, when cuts are unable to refine (PLP), the algorithm will evaluate all of $\Omega$. The lower bound at each point in $\Omega$ will be equal to the function value at that point, which will imply $l_{k}=u_{k}=\min_{x\in\Omega}f(x)$ for (PLP) (ensuring termination of the algorithm).

5 Implementation Details

Algorithm 1 relies critically on the underestimator described by (PLP). Section 5.1 develops two approaches for formulating (PLP), and Section 5.2 discusses important details for efficiently implementing Algorithm 1. Section 5.3 combines these details in a description of our preferred method for solving (1), SUCIL.

5.1 Formulating (PLP)

We present two methods for encoding (PLP) and thereby obtain lower bounds for (1). The first approach formulates (PLP) as a mixed-integer linear program (MILP) using binary variables to indicate when a point $x$ is in \mathcal{U}^{\mbox{\boldmathi}} for some multi-index $i$. Unfortunately, the resulting MILP is difficult to solve for even small problem instances. This motivates the development of the second approach, which directly builds an enumerative model of (PLP) in the space of the original variables only.

5.1.1 MILP Approach

Formulating (PLP) as an MILP requires forming the secant function m^{\mbox{\boldmathi}}(x)=(c^{\mbox{\boldmathi}})^{\operatorname{T}}x+b^{\mbox{\boldmathi}} corresponding to each multi-index \mbox{\boldmathi}\in W(X). Since m^{\mbox{\boldmathi}} is valid only in \mathcal{U}^{\mbox{\boldmathi}} (see Lemma 3), we use binary variables to encode when x\in\mathcal{U}^{\mbox{\boldmathi}}. Explicitly, for each \mbox{\boldmathi}\in W(X) and each i_{j}\in\mbox{\boldmathi}, our MILP model sets the binary variable $z^{{i}_{j}}$ to be 1 if and only if x\in\mathrm{cone}\big{(}x^{i_{j}}-X^{\mbox{\boldmathi}}\big{)}. Although the forward implication can be easily modeled by using continuous variables $\lambda^{i_{j}}$, we must introduce additional binary variables $w^{i_{j}}$ for the reverse implication.

We now describe the constraints in the MILP model. The first set of constraints ensures that $\eta$ is no smaller than any of the conditional cuts that underestimate $f$:

[TABLE]

where $M_{\eta}$ is a sufficiently large constant. By Lemma 2, we can add constraints to ensure that $x\in\Omega$ belongs to no more than one of the cones in \mathcal{U}^{\mbox{\boldmathi}} for a given $i$:

[TABLE]

The following constraints define each point $x\in\Omega$ as a linear combination of the extreme rays of each \mathrm{cone}\big{(}x^{i_{j}}-X^{\mbox{\boldmathi}}\big{)}:

[TABLE]

To indicate that x\in\mathrm{cone}\big{(}x^{i_{j}}-X^{\mbox{\boldmathi}}\big{)}, the following constraints enforce a lower bound of 0 on $\lambda$ when the corresponding $z^{i_{j}}=1$:

[TABLE]

where $M_{\lambda}$ is a sufficiently large constant. Next, we introduce the binary variables $w_{l}^{i_{j}}$ that are 1 when the corresponding variable $\lambda_{l}^{i_{j}}$ is nonnegative. The following constraints model the condition that $w_{l}^{i_{j}}=0$ implies that the corresponding $\lambda_{l}^{i_{j}}$ takes a negative value:

[TABLE]

where $\epsilon_{\lambda}$ is a sufficiently small positive constant. The last set of constraints force at least one of the $w$ variables to be 0 if the corresponding $z$ is 0:

[TABLE]

The full MILP model encoding of (PLP) is

[TABLE]

The constants $M_{\eta},M_{\lambda}$, and $\epsilon_{\lambda}$ must be chosen carefully in order to avoid numerical issues when solving (CPF); see Appendix D for more details. In early numerical results, we observed that taking large values for $M_{\eta}$ and $M_{\lambda}$ and small values for $\epsilon_{\lambda}$ resulted in numerical issues for the MILP solvers. In an attempt to remedy this situation, we derived cuts (also described in Appendix D) in which c^{\mbox{\boldmathi}} and b^{\mbox{\boldmathi}} are integer valued. One can then show, for example, that 1/\|c^{\mbox{\boldmathi}}\|_{2} is a valid lower bound for $\epsilon_{\lambda}$, and similar tight bounds can be derived for $M_{\lambda}$. With these tighter constants, some numerical issues were resolved. Yet, the growth in the number of constraints in (CPF) prevented its application to problems with $n\geq 3$.

Initial versions of the MILP model (CPF) resulted in large times to solution. Figure 3 shows the behavior of Algorithm 1—when adding all possible cuts when updating (PLP) and choosing the next iterate be a minimizer of (PLP)—when minimizing the convex quadratic function abhi (defined in Table A1 of Appendix A) on $\Omega=[-2,2]^{3}\cap\mathbb{Z}^{3}$. We note that the variations in CPU time are consistent over five repeated runs and vary by less than 2.4% for the last two iterations. The solution time of MILP solvers depends critically on implementation features, including presolve operations, node selection rules, and branching preferences. After the additional set of cuts (constraints) are introduced in iteration 12 of this problem instance, the MILP solver was able to solve the problem in slightly less time than the previous iteration. (Such occurrences are not rare in MILPs: time to solution is not strictly increasing in problem size.) Overall, we find that the growth in CPU time is due to the increasing number of conditional cuts and the associated explosion in the number of binary and continuous variables. This trend appears to limit the applicability of the MILP approach. Note that the global minimum of abhi on $[-2,2]^{3}\cap\mathbb{Z}^{3}$ has not yet been encountered when the MILPs become too large to solve. (The iteration 13 MILP was not solved in 30 minutes.)

5.1.2 Enumerative Approach

Whereas the MILP from Section 5.1.1 encodes information about every conditional cut in a single model, this section considers an alternative approach of updating the value of $\eta(x)$ for each $x\in\Omega$ as new conditional cuts are encountered. After the information from a new secant function is used to update $\eta(x)$, the secant is discarded.

Ordering the finite set of feasible integer points as $\left\{x^{1},x^{2},\ldots,x^{\left|\Omega\right|}\right\}$, our approach maintains and updates a vector of bounds

[TABLE]

where $\eta(x^{j})$ is the value of (PLP) when $x=x^{j}$. The value of $\eta(x^{j})$ is initialized to $-\infty$; and as each secant is constructed, $\eta(x^{j})$ is set to the maximum of its current value and the value of the conditional cut at $x^{j}$. This procedure is described in Algorithm 2. Since the important information about each conditional cut will be stored in $\eta(x)$, the secants defining each cut do not need to be stored. Furthermore, if $\eta_{k}(x)$ is the value of the underestimator (21) at iteration $k$, then solving each instance of (PLP) corresponds to looking up $\displaystyle\operatorname*{arg\,min}_{j\in\{1,\ldots,\left|\Omega\right|\}}\eta_{k}(x^{j})$ (breaking ties arbitrarily). Similarly, termination of Algorithm 1 requires testing only that $\displaystyle\min_{j\in\{1,\ldots,\left|\Omega\right|\}}\;\eta_{k}(x^{j})\geq u_{k}$.

Note that when solving (1), updating $\eta(x)$ for all $x\in\Omega$ is unnecessary. Rather, one needs to update $\eta(x)$ only at points that could possibly be a global minimum of $f$ on $\Omega$. When $f$ is evaluated at $x^{k+1}$ and a multi-index \mbox{\boldmathi}\in W(X^{k}\cup x^{k+1}) is encountered that is not in $W(X^{k})$, we update the lower bound only at points in \mathcal{U}^{\mbox{\boldmathi}} that are also in

[TABLE]

That is, we update $\eta_{k}(x)$ for points in \mathcal{U}_{k}^{\mbox{\boldmathi}}=\Omega_{k}\cap\mathcal{U}^{\mbox{\boldmathi}} for each newly encountered $i$.

5.2 Other Implementation Details

The enumerative approach of maintaining the value of the underestimator $\eta(x)$ described in Section 5.1.2 avoids many of the computational pitfalls of the MILP model discussed in Section 5.1.1. Below, we discuss additional computational enhancements that lead to an efficient implementation of Algorithm 1 in conjunction with Algorithm 2.

5.2.1 Checking Whether X^{\mbox{\boldmathi}} Is Poised and Whether x\in\mathcal{U}^{\mbox{\boldmathi}}

We now describe a numerically efficient representation of \mathrm{cone}\big{(}x^{i_{j}}-X^{\mbox{\boldmathi}}\big{)} for i_{j}\in\mbox{\boldmathi}. Given a poised set of $n+1$ points, X^{\mbox{\boldmathi}}, for each i_{j}\in\mbox{\boldmathi} we define a secant function satisfying

[TABLE]

Only one such secant exists for each i_{j}\in\mbox{\boldmathi}; however, the representation of this secant is not unique since $(c^{i_{j}},b^{i_{j}})$ are obtained by solving an underdetermined system of equations. Given $(c^{i_{j}},b^{i_{j}})$ satisfying (23) and (24), we define the corresponding halfspace,

[TABLE]

We now show that \mathrm{cone}\big{(}x^{j}-X^{\mbox{\boldmathi}}\big{)}, defined in (5), can be represented as the intersection of $n$ such halfspaces.

Lemma 6** (Set Equality).**

For a poised set X^{\mbox{\boldmathi}}, \mathrm{cone}\big{(}x^{i_{j}}-X^{\mbox{\boldmathi}}\big{)}=F^{i_{j}}=\bigcap\limits_{i_{l}\neq i_{j}}H^{i_{l}} for each i_{j}\in\mbox{\boldmathi}.

Proof.

Let $i$ be given and i_{j}\in\mbox{\boldmathi} fixed. We first show that \mathrm{cone}\big{(}x^{i_{j}}-X^{\mbox{\boldmathi}}\big{)}\subseteq F^{i_{j}} by showing that an arbitrary x\in\mathrm{cone}\big{(}x^{i_{j}}-X^{\mbox{\boldmathi}}\big{)} satisfies (25) for each i_{l}\in\mbox{\boldmathi},i_{l}\neq i_{j}. Given $(c^{i_{l}},b^{i_{l}})$ satisfying (23) and (24), then using the definition (5) yields

[TABLE]

where we have used (23) in the last three equations. The final inequality holds because $\lambda_{l}\geq 0$ by (5) and $(c^{i_{l}})^{\operatorname{T}}x^{i_{l}}+b^{i_{l}}>0$ by (24). Because $i_{l}$ is arbitrary, it follows that any $x$ in \mathrm{cone}\big{(}x^{i_{j}}-X\big{)} is also in $F^{i_{j}}$.

We now show that F^{i_{j}}\subseteq\mathrm{cone}\big{(}x^{j}-X^{\mbox{\boldmathi}}\big{)} by contradiction. If x\notin\mathrm{cone}\big{(}x^{i_{j}}-X^{\mbox{\boldmathi}}\big{)} for a set of $n+1$ poised points X^{\mbox{\boldmathi}}, then $x$ can be represented as $x^{i_{j}}+\displaystyle\sum_{l=1,l\neq j}^{n+1}\lambda_{l}\left(x^{i_{j}}-x^{i_{l}}\right)$ only with some $\lambda_{l}<0$. Thus, (25) is violated for some $l$, and hence $x\notin F^{i_{j}}$. ∎

Lemma 6 gives a representation of each \mathrm{cone}\big{(}x^{i_{j}}-X^{\mbox{\boldmathi}}\big{)} involving $n$ halfspaces that differs from \mathrm{cone}\big{(}x^{i_{l}}-X^{\mbox{\boldmathi}}\big{)} for i_{l}\in\mbox{\boldmathi}, $i_{l}\neq i_{j}$ in only one component. Therefore, we can represent \mathcal{U}^{\mbox{\boldmathi}} via only $n+1$ halfspaces. We efficiently calculate these halfspaces by utilizing the QR factorization \left[Q^{\mbox{\boldmathi}}\;R^{\mbox{\boldmathi}}\right]=\left[\bar{X}^{\mbox{\boldmathi}}\;e\right]^{\operatorname{T}}. If R^{\mbox{\boldmathi}} has positive diagonal entries, then the multi-index $i$ corresponds to a poised set X^{\mbox{\boldmathi}}. The coefficients in each $(c^{i_{j}},b^{i_{j}})$ can be obtained by updating Q^{\mbox{\boldmathi}}, R^{\mbox{\boldmathi}} by deleting the corresponding column from \left[\bar{X}^{\mbox{\boldmathi}}\;e\right]^{\operatorname{T}}. The sign of $(c^{i_{j}},b^{i_{j}})$ can be changed in order to ensure that (24) holds.

5.2.2 Approximating $W(X^{k}\cup\{x^{k+1}\})$

The use of $\eta_{k}(x)$ to store the lower bound at each $x\in\Omega_{k}$ allows us to avoid encoding all secants in $W(X^{k})$. After $f$ has been evaluated at a new point $x^{k+1}$, constructing the tightest possible underestimator in $\eta_{k}$ requires considering multi-indices $i$ in $W(X^{k}\cup\{x^{k+1}\})$ that contain $x^{k+1}$. (Combinations not containing $x^{k+1}$ have already been considered in previous iterations.) While not storing secants is significantly more computationally efficient than encoding and storing all secants in $W(X^{k})$, it still results in checking the poisedness of prohibitively many sets of $n+1$ points. For example, if $\left|X^{k}\right|=100$ and $n=5$, over 75 million QR factorizations must be performed, as discussed in Section 5.2.1.

Therefore, as an alternative, we seek a small, representative subset of multi-indices of $W(X^{k})$ by identifying a subset of points that will yield the best conditional cuts.

Definition 3**.**

Let $\bar{W}_{k}$ be the set of multi-indices in $W(X^{k})$ that define the largest lower bound at some point in $\Omega_{k}$ (defined in (22)). That is, \bar{W}_{k}=\{\mbox{\boldmathi}:\exists x\in\Omega_{k}\text{ such that }\eta_{k}(x)=m^{\mbox{\boldmathi}}(x)\}. We denote to the generator set of points as G^{k}=\{x^{j}:\exists\mbox{\boldmathi}\in\bar{W}_{k}\text{ such that }j\in\mbox{\boldmathi}\}.

Hence, $G^{k}$ contains points that define $\eta_{k}(x)$ for at least one $x\in\Omega_{k}$. Using $W(G^{k})$ in place of $W(X^{k})$ does relax (PLP), yet the lower bounding property of (PLP) still remains. We show below that this change does not affect the finite termination property of Algorithm 1 provided at least one cut is added for every new $x^{k+1}$.

Figure 4 compares the growth of the number of subsets of indices that must be considered when determining whether a multi-index $i$ is poised or not when using Algorithm 1 to minimize quad (Table A1 in Appendix A) on $\Omega=[-4,4]^{3}\cap\mathbb{Z}^{3}$. Preliminary numerical experiments showed that although a high percentage of all combinations in $W(X^{k}\cup x^{k+1})$, which involve the new iterate $x^{k+1}$ at an iteration $k$, are poised, only a small fraction of these actually update the lower bound at any point in $\Omega_{k}$ (we elaborate more on this in Section 7).

5.2.3 Selecting $x^{k+1}$

Early experiments with our algorithm showed that it spent many early iterations evaluating points at the boundary of $\Omega$. Although Section 3.3 provides a method for ensuring that all $x\in\Omega$ are bounded by at least one conditional cut, the solution to (PLP) is often at the boundary of $\Omega$. Rather than moving so far from a candidate solution, we consider a trust-region approach to keep iterates close to the current incumbent. As long as we maintain a lower bound for $f$ on $\Omega$, the convergence proof in Theorem 1 does not depend on $x^{k+1}$ being the global minimizer of our lower bound.

In practice, we use an infinity-norm trust region and set the minimum trust-region radius, $\Delta_{\min}$ to 1. At iteration $k$, the maximum radius that must be considered is $\max_{x,y\in\Omega_{k},x\neq y}\left\|x-y\right\|_{\infty}$.

5.3 The SUCIL Method

We now present the SUCIL method for obtaining global solutions to (1) under Assumption 1. The algorithm using the trust-region step is shown in Algorithm 3. We observe that Algorithm 3 maintains a valid lower bound $\eta_{k}(x)$ at every point, $x\in\Omega_{k}$, and that the trust-region mechanism ensures that the algorithm terminates only when the lower bound equals the best observed function value.

We note that $G^{k}$ may not be a subset of $G^{k+1}$, because $\Omega_{k}$ can contain fewer points as the upper and lower bounds on $f$ are improved. However, the following generalization of Theorem 1 ensures that Algorithm 3 still returns a global minimizer of (1).

Theorem 2** (Convergence of Algorithm 3).**

If Assumption 1 holds and if $W(G^{k})$ includes at least one cut for every $x\in\Omega_{k}$, then Algorithm 3 terminates at an optimal solution $x^{*}$ of (1) in finitely many iterations.

Proof.

Algorithm 3 will terminate in a finite number of iterations because $\Omega$ is bounded and Line 3 ensures that $x^{k}$ is not a previously evaluated element of $\Omega$. Because $W(G^{k})\subset W(X^{k})$, it follows that $\eta_{k}(x)$ is a valid lower bound for $f$ on $\Omega$, and the trust-region mechanism in Line 3 ensures that Algorithm 3 terminates only if $l_{k+1}=u_{k}$. Therefore, the result is shown. ∎

6 Numerical Experiments

We now describe numerical experiments performed on multiple versions of SUCIL; see Table 2. These methods differ in how $x^{k+1}$ is selected and in the set of points used within (PLP). The last two methods are idealized because they assume access to the true function value at every point in $\Omega_{k}$. They are included in order to provide a best-case performance for a SUCIL implementation. In the numerical experiments to follow, we set $\Delta_{\min}\leftarrow 1$ in Algorithm 3 and use an infinity-norm trust region. All SUCIL instances begin by evaluating the starting point $\bar{x}$ and $\left\{\bar{x}\pm e_{1},\ldots,\bar{x}\pm e_{n}\right\}$ ensuring a finite lower bound at every point in $\Omega$.

Below, we compare SUCIL implementations with a direct-search method, DFLINT [33]; a model-based method, MATSuMoTo [34]; and a hybrid method, NOMAD [3]. We tested the default nonmonotone DFLINT in MATLAB, as well as the monotone version, denoted DFLINT-M. We tested the default C++ version of NOMAD (v.3.9.0) as well as the same version with DISABLE MODELS set to true, denoted NOMAD-dm; the rest of the settings are default. MATSuMoTo is a surrogate-model toolbox explicitly designed for computationally expensive, black-box, global optimization problems. Since MATSuMoTo has a restarting mechanism that ensures that any budget of function evaluations will be exhausted, we input the optimal objective function value to MATSuMoTo and allowed it to run (and make as many restarts as required) until the global optimal value was identified. The default settings were used: surrogate models using cubic radial basis functions, sampling at the minimum of the surrogate, and using an initial symmetric Latin hypercube design. We performed 20 replications of MATSuMoTo for each problem instance; the details of each run are shown in Tables C4–C6 in Appendix C. We report the floor of the average number of function evaluations incurred in the last row of these tables and use this statistic for our comparisons. A common starting point is given to all methods; the starting point for the maxq and mxhilb problems is the global minimizer. A maximum function evaluation limit of 1,000 is set for all the methods when $n=4$ or $5$.

We perform numerical experiments minimizing the convex objectives in Table A1 in Appendix A on the domains $\left[-4,4\right]^{n}$ for $n\in\{3,4,5\}$ to yield 24 problem instances. (The last row of Table 1 shows $\left|\Omega\right|$ for these test problems.) Of note is the KLT function that generalizes the example function from [29] that shows how coordinate search methods can fail to find descent. The function from [29] is itself a modification of the Dennis-Woods function [15], is strongly convex, and points $x$ along the line $x_{1}=\cdots=x_{n}$ satisfy $f(x)<f(x\pm\epsilon e_{i})$ for all $i$ and for all $\epsilon>0$. The problems CB3II, CB3I, LQ, maxq, and mxhilb were introduced in [22] and also used in [33]. These five problems are either summation or maximization of generalizations of simple convex functions, constructed by extending or chaining nonsmooth convex functions or making smooth functions nonsmooth. The function LQ takes a global minimum at any $x\in[0,1]^{n}\cap\mathbb{Z}^{n}$ that does not have zeros in consecutive coordinates. For example, for $n=3$, the points $[0,1,0]^{\operatorname{T}},[0,1,1]^{\operatorname{T}},[1,0,1]^{\operatorname{T}},[1,1,0]^{\operatorname{T}}$, and $[1,1,1]^{\operatorname{T}}$ are optimal but $[0,0,0]^{\operatorname{T}},[0,0,1]^{\operatorname{T}}$, and $[1,0,0]^{\operatorname{T}}$ are not.

If there are relatively few points in $\Omega$ and the time required to evaluate $f$ is small (as for our test instances), one could argue that an enumerative procedure itself could solve the problem in a reasonable time. However, we use these instances to thoroughly examine the behavior that might be seen on expensive-to-evaluate black-box functions. Therefore, we compare methods using performance profiles [16] that are based on the number of function evaluations required to satisfy the respective convergence criterion. For each method $s$, $\rho_{s}(\alpha)=\frac{\lvert\{p\in P:r_{p,s}\leq\alpha\}\rvert}{\lvert P\rvert}$, for a scalar $\alpha\geq 1$, $P$ is the collection of problems, and $r_{p,s}=\frac{N_{p,s}}{\min_{s\in S}\{N_{p,s}\}}$ is the performance ratio. We consider two measures of $N_{p,s}$: (1) the number of function evaluations before a method $s$ terminates on a problem $p$ and (2) the number of function evaluations taken by method $s$ to evaluate a global minimizer on problem $p$.

Figure 2 compares the number of evaluations required for four implementations of SUCIL to terminate (with a certificate of optimality) on the set of test problems. While SUCIL-ideal1 is no slower than any other implementation on all the test problems, it is not a realistic method in that it evaluates points based on their known function values. SUCIL requires no more than three times the evaluations as SUCIL-ideal1 for the set of test problems. We do observe that using a trust region in SUCIL is a significant advantage. For many of the problems considered, SUCIL-noTR spent many function evaluations in the corners of $\Omega$.

As a point of comparison with the results in Figure 2, a different estimate of the number of function evaluations (or primitive directions explorations) required for the proof of optimality for our instances can be seen in Table 1, in columns corresponding to $n\in\{3,4,5\}$ and $k=4$. As evident from the results in Tables C1–C3, our method incurs a remarkably low number of function evaluations, which can be attributed to exploitation of convexity and subsequent formation of the underestimators, as explained in Section 2.

We now analyze the performance of SUCIL compared with the other methods. Figure 5 (left) shows the performance profiles for methods to terminate on the 24 test problems; Figure 5 (right) compares the number of function evaluations required before each method first evaluates a global minimizer. This comparison is nontrivial because each solver has its own design considerations and notions of local optimality. Also, the other solvers do not assume convexity of the problem or exploit it. Hence, our results merely demonstrate that SUCIL provably converges to a global optimum and uses fewer function evaluations because of its exploitation of convexity. Figure 5 shows that our algorithm requires the least number of function evaluations for more than 65% of the instances and provides a global optimality certificate, in addition. In reaching the global optimal solution quickly, however, DFLINT-M wins for more than 85% of the instances. Although SUCIL is not particularly designed to greedily descend to the global optimum, it is still competitive with the rest of the methods on this front.

7 Discussion

The order of results in this paper tells the story of how we arrived at the implementation of SUCIL. We first attempted to classify where linear interpolation models provide lower bounds for convex functions, yielding the results in Section 3; we then proved that such linear functions can underlie a convergent algorithm, as in Section 4. We initially modeled the secants and the conditions in which they are valid as an MILP, as in Section 5.1.1. After observing that the number of variables in the MILP model was larger than the number of points in the domain $\Omega$, we were motivated to develop the enumerative model in Section 5.1.2.

Our computational developments expose a number of fundamental challenges for integer derivative-free optimization. The complexity of our piecewise linear model (PLP) is made worse by the fact that each secant function is valid only in the union of $n+1$ cones \mathcal{U}^{\mbox{\boldmathi}}, resulting in conditional cuts. We note that it may not be possible to derive unconditional cuts, that is, cuts that are valid in the whole domain $\Omega$. For example, we might initially consider secants interpolating a convex $f$ at the $n+1$ points $x\in\mathbb{Z}^{n}$ and $x\pm e_{i}\in\mathbb{Z}^{n}$, where for every $i$ we can choose either $+$ or $-$. Such points form a unit simplex that has no integer points in its interior. Consequently, one might suspect that the resulting cut is valid everywhere in $\Omega$. However, the following example shows that the resulting cut is not unconditionally valid. Consider $f(x)=x_{1}^{2}-x_{1}x_{2}+x_{2}^{2}$ and the set of points $\{[1,1]^{\operatorname{T}},[0,1]^{\operatorname{T}},[1,0]^{\operatorname{T}}\}$. It follows that $f(x)=1$ at these points, and hence the unique interpolating secant function is the constant function, $m(x)=1$. Now consider the point $x=[0,0]^{\operatorname{T}}$ for which $f(x)=0$, which is not underestimated by $m(x)=1$. Another limitation of our method is that it will have to evaluate all feasible points when $\Omega=\{0,1\}^{n}$ (a pure binary domain) for convergence because no point in \Omega\setminus X^{\mbox{\boldmathi}} belongs to \mathcal{U}^{\mbox{\boldmathi}}, where X^{\mbox{\boldmathi}} is an arbitrary poised set of $n+1$ points in $\{0,1\}^{n}$.

In Figure 6 we show the number of function evaluations needed to first evaluate a global minimizer and the additional number of evaluations used to prove it is a global minimizer. As is common, the effort required to certify optimality can be significantly larger than the cost of finding the optimum. In terms of number of function evaluations required, the proof of optimality is even more time consuming. Because the iterations where $X^{k}$ or $G^{k}$ is large require checking many potential secant functions, in SUCIL the computational cost of iterations can differ by orders of magnitude as the algorithm progresses.

Although our method provides a practical iterative way to check sufficiency of a set of points (optimality conditions) for a given convex instance, each iteration involves construction and evaluation of a large number of combinations of different $n+1$ points, which limits the scalability of Algorithm 3 in solving instances of higher dimensions. Yet, in our numerical experiments, we observe that only a small fraction of the total cuts evaluated are useful. We call (c^{\mbox{\boldmathi}},b^{\mbox{\boldmathi}}) an updating cut at an iteration $k$ if there exists an $x\in\Omega_{k}$ such that m^{\mbox{\boldmathi}}(x)>\eta_{k}(x), that is, a cut that improves the lower bound at at least one $x\in\Omega_{k}$. In addition, if m^{\mbox{\boldmathi}}(x)\geq u_{k}, we call it a pruning cut. A pruning cut helps eliminate points to be considered in the next iteration ($\Omega_{k+1}$). Figure 7 shows the number of updating and pruning cuts generated per iteration of SUCIL when minimizing quad on $[-4,4]^{3}\cap\mathbb{Z}^{3}$. The fact that few cuts prune a point or update the lower bound at any point where the minimum could be suggests that there may be some way to exclude a large set of multi-indices from consideration, possibly yielding dramatic computational savings.

Ideally, we would like to evaluate only the combinations that yield updating or pruning cuts. However, this approach requires the solution of a separate problem that we believe is especially hard to solve. Even the following simpler problem of finding a pruning cut at a given candidate point seems difficult.

Problem 1**.**

Given a point $\bar{x}\in\mathbb{Z}^{n}$, a set of (integer) points $X$ where $f$ has been evaluated, and scalar $u$, find a cut that prunes $\bar{x}$. That is, find a multi-index $i$ such that \bar{x}\in\mathcal{U}^{\mbox{\boldmathi}} and {(c^{\mbox{\boldmathi}})}^{\operatorname{T}}\bar{x}+b^{\mbox{\boldmathi}}\geq u, and (c^{\mbox{\boldmathi}},b^{\mbox{\boldmathi}}) solves (3), or show that no such multi-index $i$ exists.

If we choose a small subset $\bar{X^{k}}$ of $X^{k}$ to form $W(\bar{X}^{k})$, the SUCIL algorithm can end up using a large number of function evaluations to obtain a certificate of optimality. The reason is that points are evaluated that would be ruled worse than optimal if secants were built by using all combinations of points in $X^{k}$. This situation occurred when setting $\bar{X^{k}}$ to be a random subset of $X^{k}$, a subset of the points closest to $\hat{x}$, or a subset of points with best function values. Using $G^{k}$ avoids discarding too many points from $X^{k}$; but we observe a significant increase in $\left|W(G^{k})\right|$, and thus we incur heavy computational costs during some iterations. The wall-clock time required per iteration for solving instances of dimension less than 5 in our setup is not significant, but we present the same for 5-dimensional instances using SUCIL on a 96-core Intel Xeon computer with 1.5 TB of RAM. The complexity of our approach is better quantified by counting the number of combinations of points (or potential secants) considered at iteration $k$. Using $G^{k}$, we typically produce a strict subset of all possible combinations in such a way that the size of $W(G^{k})$ decreases during the later iterations. This is shown in Figure 8: the number of secants added per iteration for all 5-dimensional test instances using SUCIL. Once $\Omega_{k}$, the number of points with $\eta(x)$ less than $f(\hat{x})$, starts decreasing, so do $G^{k}$ and $\left|W(G^{k})\right|$. In general, it is difficult to predict when the number of combinations (or the wall-clock time curve) would be at the peak, but we suspect this peak will be worse as $n$ increases, by both the size and the iteration number where it occurs. This limits the applicability of the current implementation of SUCIL on higher-dimensional problems.

Again, since nearly all cuts in $W(G^{k})$ do not update $\eta(x)$ at any point in $\Omega_{k}$ (see Figure 7), we believe there may be some approach for intelligently selecting points from $X^{k}$ using their geometry, their function values, and distance from $\hat{x}$ that will rule some multi-indices $i$ as unnecessary to consider. We did attempt to identify minimal sets of points that were necessary for SUCIL to certify optimality for a variety of $n=2$ test cases, but no general rule was apparent.

We note that the storage requirements for the enumerative model may be prohibitive, even for moderate problem sizes. For example, an array storing the value of $\eta(x)$ as an 8-byte scalar for all $x\in\Omega=[-10,10]^{10}\cap\mathbb{Z}^{10}$ would require over 200 GB of storage.

Ultimately, we believe further insights are yet to be discovered that will facilitate better algorithms for minimizing convex functions on integer domains.

Acknowledgements

We are grateful to Eric Ni for his insights on derivative-free algorithms for unrelaxable integer variables. Sven Leyffer also wishes to acknowledge the insightful discussions on an early draft of this work during the Dagstuhl seminar 18081. This material is based upon work supported by the applied mathematics and SciDAC activities of the Office of Advanced Scientific Computing Research, Office of Science, U.S. Department of Energy, under Contract DE-AC02-06CH11357.

Appendix A Test Problems

Appendix B Performance of MILP Model

Table B1 shows the size of the MILP model at each iteration and the computational effort required for solving it. The column k refers to the iteration of Algorithm 1, sHyp denotes the number of secants (i.e., $\left|W(X)\right|$), and LB and UB give the lower and upper bound on $f$ on $\Omega$, respectively. We show the computational effort needed to solve each MILP via time, the mean solution time (in seconds) for $5$ replications; simIter, the number of simplex iterations; and nodes, the number of branch-and-bound nodes explored by the MILP solver. The size of each MILP (after presolve) is shown in terms of bVars, the number of binary variables; cVars, the number of continuous variables; and cons, the number of constraints. Table B1 also shows the optimal solution $\hat{x}$ of each MILP. These experiments were performed by using CPLEX (v.12.6.1.0) on a 2.20 GHz, 12-core Intel Xeon computer with 64 GB of RAM. For this small problem we see that the size of the MILP grows exponentially as the iterations proceed, which results in an exponential growth in solution time as illustrated in Figure 3. The iteration 13 MILP was not solved after $30$ minutes.

Appendix C Detailed Numerical Results

Tables C1–C6 contain detailed numerical results for the interested reader. Note that some solvers do not respect the given budget of function evaluations. We have used a different stopping criterion for MATSuMoTo: it is set to stop only when a point with the optimal value has been identified. Also, although the global minimum is the starting point for maxq and mxhilb, MATSuMoTo instead uses its initial symmetric Latin hypercube design. The last row of Table 1 in Section 2 shows $\left|\Omega\right|$ for these problems.

Appendix D Derivation of Tolerances

In this section, we show how one can compute sufficient values of $M_{\eta}$, $\epsilon_{\lambda}$, and $M_{\lambda}$ for the MILP model (CPF). We first show that if we choose these parameters incorrectly, then the resulting MILP model no longer provides a valid lower bound.

Effect of an Insufficient $\epsilon_{\lambda}$ Value.

A large $\epsilon_{\lambda}$ or small $M_{\lambda}$ (or both) could result in an incorrect value of $z^{i_{j}}=1$ for a point x\notin\mathrm{cone}\big{(}x^{i_{j}}-X^{\mbox{\boldmathi}}\big{)}, violating the implication of $z^{i_{j}}=1$ and yielding an invalid lower bound on $f$. This is illustrated using a one-dimensional example in Figure D1. Similarly, one can encounter an invalid lower bound at an iteration if $M_{\lambda}$ is chosen to be smaller than required.

Bound on $M_{\eta}$

Let $l_{f}$ be a valid lower bound of $f$ on $\Omega$. With this lower bound, scalars M_{\mbox{\boldmathi}} can be defined as

[TABLE]

If $\Omega=\left\{x:l_{x}\leq x\leq u_{x}\right\}$, where $l_{x},u_{x}\in\mathbb{R}^{n}$ are known, then we can set

[TABLE]

Then, we can either set individual values of M_{\mbox{\boldmathi}} within each constraint (15), which would yield a tighter model, or use a single parameter,

[TABLE]

as shown in (CPF).

Bounds on $\epsilon_{\lambda}$ and $M_{\lambda}$

Sufficient values of $M_{\lambda}$ and $\epsilon_{\lambda}$ are not easy to calculate as they depend on the rays generated at $x^{i_{j}}$, and the domain $\Omega$. We show next, that bounds on these values can be computed by solving a set of optimization problems. A sufficient value for $\epsilon_{\lambda}$ can be computed based on the representations of the $n+1$ hyperplanes formed using different combinations of $n$ points, \forall\mbox{\boldmathi}\in W(X). We can use one of the representations of the hyperplane passing through points $X^{i}\setminus\{x^{i_{j}}\},j=1\ldots n+1$ to obtain bounds on $\epsilon_{\lambda}$, by solving an optimization problem for each i_{j}\in\mbox{\boldmathi}. Let c^{\mbox{\boldmathi}_{i_{j}-}} and b^{\mbox{\boldmathi}_{i_{j}-}} denote the solution of the following problem.

[TABLE]

Problem (P$-hyp$) can be easily cast as an integer program by replacing $\left\|c\right\|_{1}$ by $e^{T}y$, where $e=(1,\ldots,1)^{T}$ is the vector or all ones, and the variables $y$ satisfy the constraints $y\geq c,y\geq-c$. Because, the hyperplane is generated using $n$ integer points, it can be shown that there exists a solution c^{\mbox{\boldmathi}_{i_{j}-}} and b^{\mbox{\boldmathi}_{i_{j}-}} that is nonzero and integral, as stated in the following proposition.

Proposition D.1**.**

Consider a hyperplane $S=\{x:c^{T}x+b=0\}$ such that $c$ and $b$ are integral. The Euclidean distance between an arbitrary point $\hat{x}\in\mathbb{Z}^{n}\setminus S$, and $S$, is greater than or equal to $\frac{1}{\lvert\lvert c\rvert\rvert_{2}}$.

Proof.

The Euclidean distance between a point $\hat{x}$ and $S$ is the 2-norm of the projection of the line segment joining an arbitrary point $w\in S$ and $\hat{x}$, on the normal passing through $\hat{x}$, and can be expressed as

[TABLE]

By integrality of $c,b$ and $\hat{x}$, and since $\hat{x}\notin S$, $\lvert c^{T}\hat{x}+b\rvert\geq 1$, and the result follows. ∎∎

Proposition D.1 can be used directly to get a sufficient value of $\epsilon_{\lambda}$ as follows.

[TABLE]

The bound can be tightened using the following optimization problem.

[TABLE]

If we denote by \epsilon^{\mbox{\boldmathi}_{i_{j}-}}_{\lambda}, the optimal value of (P-$\epsilon$), then the following is a sufficient value of $\epsilon_{\lambda}$.

[TABLE]

Similarly, if we maximize the objective in (P-$\epsilon$), and denote the optimal value by M^{\mbox{\boldmathi}_{i_{j}-}}_{\lambda} we get a sufficient value for $M_{\lambda}$ as follows.

[TABLE]

No-Good Cuts and Stronger Objective Bounds.

For our initial experiments, we tried setting arbitrarily small values for $\epsilon_{\lambda}$ instead of obtaining the best possible values by solving a set of optimization problems as elaborated above. The problem with having such values of $\epsilon_{\lambda}>0$ too small is that it allows us to ignore the cuts at the interpolation points, $x^{(k)}$, causing a violation of the bound, $\eta\geq f^{(k)}$, for $x=x^{(k)}$ in the MILP, resulting in repetition of iterates in the algorithm. We fixed this by adding the following valid inequality to the MILP.

[TABLE]

We can write this inequality as a linear constraint, by introducing a binary representation of the variables, $x$, requiring $nU$ binary variables, $\xi_{ij},\;i=1,\ldots,n,\;j=1,\ldots,U$, where $n$ is the dimension of our problem, and $0\leq x\leq U$. With this representation, we obtain

[TABLE]

and write the $\eta$-constraint equivalently as

[TABLE]