Separable Convex Optimization with Nested Lower and Upper Constraints

Thibaut Vidal; Daniel Gribel; Patrick Jaillet

arXiv:1703.01484·math.OC·September 11, 2018·INFORMS J. Optim.

Separable Convex Optimization with Nested Lower and Upper Constraints

Thibaut Vidal, Daniel Gribel, Patrick Jaillet

PDF

Open Access

TL;DR

This paper introduces an efficient, gradient-free divide-and-conquer algorithm for solving convex resource allocation problems with nested bounds, achieving state-of-the-art complexity and practical performance for large-scale instances.

Contribution

The paper presents a novel divide-and-conquer algorithm that handles nested lower and upper constraints without requiring strict convexity or differentiability, improving computational efficiency.

Findings

01

Algorithm solves large problems with up to 1,000,000 variables in seconds.

02

Achieves optimal complexity bounds for the problem class.

03

Demonstrates promising applications to support vector ordinal regression.

Abstract

We study a convex resource allocation problem in which lower and upper bounds are imposed on partial sums of allocations. This model is linked to a large range of applications, including production planning, speed optimization, stratified sampling, support vector machines, portfolio management, and telecommunications. We propose an efficient gradient-free divide-and-conquer algorithm, which uses monotonicity arguments to generate valid bounds from the recursive calls, and eliminate linking constraints based on the information from sub-problems. This algorithm does not need strict convexity or differentiability. It produces an $ϵ$ -approximate solution for the continuous problem in $O (n lo g m lo g \frac{n B}{ϵ})$ time and an integer solution in $O (n lo g m lo g B)$ time, where $n$ is the number of decision variables, $m$ is the number of constraints, and…

Figures14

Click any figure to enlarge with its caption.

Tables3

Table 1. Table 1: Detailed average CPU times for experiments with a linear objective

Variable $𝐦$		CPU Time(s)		Fixed $𝐦$		CPU Time(s)
$n$	$m$	FLOW	MDA	$n$	$m$	FLOW	MDA
10	10	$0.000 002 753$	$0.000 004 784$	100	100	$0.000 050 934$	$0.000 059 476$
20	20	$0.000 006 262$	$0.000 010 242$	200	100	$0.000 135 978$	$0.000 126 415$
50	50	$0.000 021 454$	$0.000 028 513$	500	100	$0.000 394 375$	$0.000 285 976$
100	100	$0.000 050 576$	$0.000 058 915$	1000	100	$0.000 907 237$	$0.000 551 754$
200	200	$0.000 126 376$	$0.000 125 758$	2000	100	$0.002 065 989$	$0.001 139 486$
500	500	$0.000 371 994$	$0.000 336 22$	5000	100	$0.006 162 977$	$0.002 958 873$
1000	1000	$0.000 842 986$	$0.000 757 082$	10000	100	$0.014 429 2$	$0.006 262 521$
2000	2000	$0.001 865 427$	$0.001 740 006$	20000	100	$0.031 677 652$	$0.015 715 67$
5000	5000	$0.005 426 528$	$0.005 199 759$	50000	100	$0.092 668 418$	$0.052 614 627$
10000	10000	$0.012 305 958$	$0.011 198 173$	100000	100	$0.203 549 512$	$0.108 438 407$
20000	20000	$0.026 222 298$	$0.032 138 957$	200000	100	$0.441 196 771$	$0.235 623 163$
50000	50000	$0.079 430 967$	$0.105 266 726$	500000	100	$1.195 609 09$	$0.718 873 748$
100000	100000	$0.151 968 169$	$0.225 724 083$	1000000	100	$2.564 233 333$	$1.598 551 587$
200000	200000	$0.366 558 44$	$0.485 572 598$
500000	500000	$0.967 679 801$	$1.369 780 555$
1000000	1000000	$1.989 558 333$	$2.9783$

Table 2. Table 2: Detailed average CPU-time for experiments with a separable convex objective

		CPU Time(s) – MDA			CPU Time(s) – MOSEK
$n$	$m$	[F]	[Crash]	[Fuel]	[F]	[Crash]	[Fuel]
10	10	$0.000 052 804$	$0.000 032 711$	$0.000 061 129$	$0.007 686 943$	$0.007 832 116$	$0.008 061 434$
20	20	$0.000 114 081$	$0.000 073 175$	$0.000 133 301$	$0.008 274 424$	$0.008 604 278$	$0.008 640 814$
50	50	$0.000 380 328$	$0.000 263 063$	$0.000 445 094$	$0.009 952 053$	$0.010 301 065$	$0.010 400 973$
100	100	$0.000 803 668$	$0.000 538 618$	$0.000 930 292$	$0.017 299 139$	$0.017 526 892$	$0.017 409 55$
200	200	$0.001 934 188$	$0.001 231 371$	$0.002 160 542$	$0.063 120 002$	$0.062 191 123$	$0.062 980 577$
500	500	$0.005 448 312$	$0.003 550 539$	$0.006 211 667$	$0.778 735 352$	$0.756 403 467$	$0.785 640 721$
1000	1000	$0.012 697 424$	$0.008 608 136$	$0.014 253 454$	$6.314 833 334$	$6.293 333 333$	$6.368 833 333$
2000	2000	$0.028 767 787$	$0.018 728 026$	$0.031 944 613$	$85.749$	$93.761$	$90.474$
5000	5000	$0.092 704 885$	$0.060 521 584$	$0.098 564 133$	$1703.926$	$1605.518$	$1554.495$
10000	10000	$0.201 499 218$	$0.133 855 399$	$0.212 543 624$	—	—	—
20000	20000	$0.469 045 625$	$0.304 190 384$	$0.482 383 221$	—	—	—
50000	50000	$1.305 475 757$	$0.873 814 239$	$1.327 877 777$	—	—	—
100000	100000	$3.1194$	$2.023 380 952$	$3.068 65$	—	—	—
200000	200000	$6.683 166 667$	$4.584$	$6.613 666 667$	—	—	—
500000	500000	$19.838$	$13.538$	$19.101$	—	—	—
1000000	1000000	$45.403$	$31.004$	$42.991$	—	—	—

Table 3. Table 3: SVOREX resolution – impact of the working-set size and solution features

Instance	$𝐍$	$𝐃$	Solution Variables s.t.			$𝐧_{ws}$	$𝐈_{𝐰𝐬}$	T(s)
Instance	$𝐍$	$𝐃$	$α = 0$	$α = C$	$α \in] 0, C [$	$𝐧_{ws}$	$𝐈_{𝐰𝐬}$	T(s)
Abalone	1000	8	39%	32%	29%	2	118233	13.46
						4	96673	21.51
						6	78433	26.34
						10	60605	35.46
Bank	3000	32	25%	0%	75%	2	139468	68.41
						4	52073	63.02
						6	31452	45.22
						10	21310	47.66
Boston	300	13	41%	0%	59%	2	7207	0.43
						4	3697	0.40
						6	2840	0.46
						10	2076	0.54
California	5000	8	51%	43%	6%	2	250720	124.46
						4	189289	185.79
						6	166879	245.08
						10	146170	360.52
Census	6000	16	38%	4%	59%	2	349894	242.11
						4	206951	301.74
						6	180608	393.28
						10	155731	574.28
Computer	4000	21	64%	32%	4%	2	290207	168.94
						4	140270	161.45
						6	98948	153.56
						10	68616	193.10
Machine CPU	150	6	49%	9%	41%	2	28856	1.24
						4	11534	0.86
						6	8144	0.91
						10	6363	1.24
Pyrimidines	50	27	21%	0%	79%	2	935	0.035
						4	367	0.021
						6	218	0.018
						10	144	0.023

Equations109

min f (x) =

min f (x) =

s.t. a_{i} \leq

k = 1 \sum n x_{k} = B

c_{i} \leq

min f (x) =

min f (x) =

s.t. a_{i} \leq

k = 1 \sum n x_{k} = B

c_{i} \leq

min f (x, I) =

min f (x, I) =

s.t. I_{i}

I_{0}

0

0

min f (x) =

min f (x) =

k = 1 \sum i d_{k} - K \leq k = 1 \sum i x_{k} \leq k = 1 \sum i d_{k} + I_{i}^{\textsc ma x} - K

0 \leq x_{i} \leq x_{i}^{\textsc ma x}

min

min

a_{i} \leq t_{i} \leq b_{i}

t_{i - 1} + \frac{δ _{i - 1, i}}{v _{i - 1, i}} \leq t_{i}

v_{min} \leq v_{i - 1, i} \leq v_{ma x}

min

min

a_{i} \leq k = 1 \sum i x_{k} \leq b_{i}

\frac{δ _{i - 1, i}}{v _{ma x}} \leq x_{i}

with g_{i} (v) = {f_{i} (v_{i}^{\textsc o pt}) f_{i} (v) if v \leq v_{i}^{\textsc o pt} otherwise.

with g_{i} (v) = {f_{i} (v_{i}^{\textsc o pt}) f_{i} (v) if v \leq v_{i}^{\textsc o pt} otherwise.

min

min

i = 1 \sum n \frac{N _{i}^{2} σ _{i}^{2}}{N ^{2}} (\frac{1}{x _{i}} - \frac{1}{N _{i}}) \leq V

0 \leq

i \in S_{i} \sum \frac{N _{i}^{2} σ _{i}^{2}}{N ^{2}} (\frac{1}{x _{i}} - \frac{1}{N _{i}}) \leq V_{i}, i \in {1, \dots, m},

i \in S_{i} \sum \frac{N _{i}^{2} σ _{i}^{2}}{N ^{2}} (\frac{1}{x _{i}} - \frac{1}{N _{i}}) \leq V_{i}, i \in {1, \dots, m},

α, α^{*}, μ max j = 1 \sum r i = 1 \sum n^{j} (α_{i}^{j} + α_{i}^{* j}) - \frac{1}{2} j = 1 \sum r i = 1 \sum n^{j} j^{'} = 1 \sum r i^{'} = 1 \sum n^{j^{'}} (α_{i}^{* j} - α_{i}^{j}) (α_{i^{'}}^{* j^{'}} - α_{i^{'}}^{j^{'}}) K (x_{i}^{j}, x_{i^{'}}^{j^{'}}) \vspace * - 0.45 c m

α, α^{*}, μ max j = 1 \sum r i = 1 \sum n^{j} (α_{i}^{j} + α_{i}^{* j}) - \frac{1}{2} j = 1 \sum r i = 1 \sum n^{j} j^{'} = 1 \sum r i^{'} = 1 \sum n^{j^{'}} (α_{i}^{* j} - α_{i}^{j}) (α_{i^{'}}^{* j^{'}} - α_{i^{'}}^{j^{'}}) K (x_{i}^{j}, x_{i^{'}}^{j^{'}}) \vspace * - 0.45 c m

0 \leq α_{i}^{j} \leq C

0 \leq α_{i}^{j} \leq C

0 \leq α_{i}^{* j} \leq C

i = 1 \sum n^{j} α_{i}^{j} + μ^{j} = i = 1 \sum n^{j + 1} α_{i}^{* j + 1} + μ^{j + 1}

μ^{j} \geq 0

k = 1 \sum j i = 1 \sum n^{k} α_{i}^{k} - i = 1 \sum n^{k + 1} α_{i}^{* k + 1} \geq 0 j \in {1, \dots, r - 2}

k = 1 \sum j i = 1 \sum n^{k} α_{i}^{k} - i = 1 \sum n^{k + 1} α_{i}^{* k + 1} \geq 0 j \in {1, \dots, r - 2}

k = 1 \sum r - 1 i = 1 \sum n^{k} α_{i}^{k} - i = 1 \sum n^{k + 1} α_{i}^{* k + 1} = 0.

{max i = 1 \sum n x_{i} μ_{i}; min i = 1 \sum n j = 1 \sum n x_{i} x_{j} σ_{ij}}

{max i = 1 \sum n x_{i} μ_{i}; min i = 1 \sum n j = 1 \sum n x_{i} x_{j} σ_{ij}}

i = 1 \sum n x_{i} = 1

0 \leq x_{i}

\displaystyle\text{RAP--NC}_{v,w}(L,R):\hskip 8.5359pt\text{min}\ \

\displaystyle\text{RAP--NC}_{v,w}(L,R):\hskip 8.5359pt\text{min}\ \

s.t. a_{i} - L \leq

i = σ [v - 1] + 1 \sum σ [w] x_{i} = R - L

with \overset{ˉ}{f}_{i} (x) = ⎩ ⎨ ⎧ f_{i} (c_{i}) + M (c_{i} - x) f_{i} (x_{i}) f_{i} (d_{i}) + M (x - d_{i}) if x < c_{i} if x \in [c_{i}, d_{i}] if x > d_{i}

k = σ [v - 1] + 1 \sum s x_{k}^{*} > k = σ [v - 1] + 1 \sum s x_{k} or k = s \sum σ [w] x_{k}^{*} > k = s \sum σ [w] x_{k} . \vspace * - 0.2 c m

k = σ [v - 1] + 1 \sum s x_{k}^{*} > k = σ [v - 1] + 1 \sum s x_{k} or k = s \sum σ [w] x_{k}^{*} > k = s \sum σ [w] x_{k} . \vspace * - 0.2 c m

Φ_{i} = k \in {v, \dots, w} ∣ σ [k] \geq i \sum (κ_{k} - λ_{k}) \in \partial \overset{ˉ}{f}_{i} (x_{i})

Φ_{i} = k \in {v, \dots, w} ∣ σ [k] \geq i \sum (κ_{k} - λ_{k}) \in \partial \overset{ˉ}{f}_{i} (x_{i})

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Taxonomy

TopicsSparse and Compressive Sensing Techniques · Machine Learning and Algorithms · Complexity and Algorithms in Graphs

Full text

Separable Convex Optimization with Nested Lower and Upper Constraints

Thibaut Vidal, Daniel Gribel, Patrick Jaillet

Separable Convex Optimization with Nested Lower and Upper Constraints

Thibaut Vidal, Daniel Gribel*

Departamento de Informática,

Pontifícia Universidade Católica do Rio de Janeiro (PUC-Rio)

[email protected]

Patrick Jaillet

Department of Electrical Engineering and Computer Science,

Laboratory for Information and Decision Systems,

Operations Research Center, Massachusetts Institute of Technology

[email protected]

August 2018

Author Accepted Manuscript (AAM)

Accepted for publication in INFORMS Journal on Optimization

Abstract. We study a convex resource allocation problem in which lower and upper bounds are imposed on partial sums of allocations. This model is linked to a large range of applications, including production planning, speed optimization, stratified sampling, support vector machines, portfolio management, and telecommunications. We propose an efficient gradient-free divide-and-conquer algorithm, which uses monotonicity arguments to generate valid bounds from the recursive calls, and eliminate linking constraints based on the information from sub-problems. This algorithm does not need strict convexity or differentiability. It produces an $\epsilon$ -approximate solution for the continuous problem in ${\mathcal{O}}(n\log m\log\frac{nB}{\epsilon})$ time and an integer solution in ${\mathcal{O}}(n\log m\log B)$ time, where $n$ is the number of decision variables, $m$ is the number of constraints, and $B$ is the resource bound. A complexity of ${\mathcal{O}}(n\log m)$ is also achieved for the linear and quadratic cases. These are the best complexities known to date for this important problem class. Our experimental analyses confirm the good performance of the method, which produces optimal solutions for problems with up to 1,000,000 variables in a few seconds. Promising applications to the support vector ordinal regression problem are also investigated.

Keywords. Convex optimization, resource allocation, nested constraints, speed optimization, lot sizing, stratified sampling, machine learning, support vector ordinal regression

Corresponding author

1 Introduction

Resource allocation problems involve the distribution of a fixed quantity of a resource (e.g., time, workforce, money, energy) over a number of tasks in order to optimize a value function. In its most fundamental form, the simple resource allocation problem (RAP) is formulated as the minimization of a separable objective subject to one linear constraint representing the total resource bound. Despite its apparent simplicity, this model has been the focus of a considerable research effort over the years, with more than a hundred articles, as underlined by the surveys of Patriksson (2008), Katoh et al. (2013), and Patriksson et Strömberg (2015). This level of interest arises from its applications in engineering, production and manufacturing, military operations, machine learning, financial economics, and telecommunications, among many other areas.

In several applications, a single global resource bound is not sufficient to model partial budget or investment limits, release dates and deadlines, or inventory or workforce limitations. In these situations, the problem must be generalized to include additional constraints over nested sums of the resource variables. This often leads to the model given in Equations (1)–(4), where the sets $J_{i}\subseteq\{1,\dots,n\}$ follow a total order such that $J_{i}\subset J_{i+1}$ for $i\in\{1,\dots,m-2\}$ :

[TABLE]

This problem involves a separable convex objective, subject to lower and upper bounds on nested subsets of the variables (Equation 2) and a global resource constraint (Equation 3). Despite being a special case of the former inequalities, the latter constraint is included in the model to emphasize the resource bound $B$ . Re-ordering the indices of $x_{i}$ , $c_{i}$ , $d_{i}$ and $f_{i}$ , we obtain the formulation given by Equations (5)–(8), where $(\sigma[1],\dots,\sigma[m])$ is a subsequence of $(1,\dots,n)$ :

[TABLE]

We assume that the functions $f_{i}:[c_{i},d_{i}]\rightarrow\Re$ are Lipschitz continuous but not necessarily differentiable or strictly convex, and the coefficients $a_{i}$ , $b_{i}$ , $c_{i}$ , and $d_{i}$ are integers. To ease the presentation, we define $a_{m}=b_{m}=B$ , $\sigma[0]=0$ and $\sigma[m]=n$ . We will study this continuous optimization problem as well as its restriction to integer solutions.

We refer to this problem as the RAP with nested lower and upper constraints (RAP–NC). As highlighted in Section 2, the applications of this model include production and capacity planning (Love 1973), vessel speed optimization (Psaraftis et Kontovas 2013, 2014), machine learning (Chu et Keerthi 2007), portfolio management (Bienstock 1996), telecommunications (D’Amico et al. 2014) and power management (Gerards et al. 2016). Some of these applications involve large data sets with millions of variables, and in other contexts multiple RAP–NC must be repeatedly solved (e.g., thousands or millions of times) to produce bounds in a tree-search-based algorithm, to optimize vessel speeds over candidate routes within a heuristic search for ship routing, or to perform projection steps in a subgradient procedure for a nonseparable objective. In these situations, complexity improvements are a determining factor between success and failure.

The literature contains a rich set of studies and algorithms for a closely related problem: the NESTED resource allocation problem, a special case of the RAP–NC in which $a_{i}=-\infty$ for all $i\in\{1,\dots,m-1\}$ (or $b_{i}=\infty$ for all $i\in\{1,\dots,m-1\}$ ). With integer variables, NESTED can be solved in ${\mathcal{O}}(n\log n\log\frac{B}{n})$ time using a scaling algorithm (Hochbaum 1994), and ${\mathcal{O}}(n\log m\log\frac{B}{n})$ time using divide-and-conquer principles (Vidal et al. 2016). These algorithms, however, are not applicable for joint lower and upper nested constraints, a case which is essential for a large variety of applications, e.g., to model time windows, time-dependent inventory bounds, or investment ranges.

Moreover, in the presence of continuous variables, the notion of computational complexity for convex problems must be carefully defined, since optimal solutions can be irrational and thus not representable in a bit-size computational model. We will use the same conventions as Hochbaum (1994), by measuring the computational complexity of achieving an $\epsilon$ -approximate solution, guaranteed to be located in the solution space no further than $\epsilon$ from an optimal solution. We also assume that an oracle is available to evaluate each function $f_{i}$ in ${\mathcal{O}}(1)$ time. When considering such a model of computation, controlling algorithmic approximations can be a hard task. To circumvent this issue, we will use, as in Hochbaum (1994), a proximity theorem to transform a continuous problem into an integer problem scaled by an appropriate factor, and to translate the integer solution back to a continuous solution with the desired precision.

The main contributions of this paper are the following:

•

We propose an efficient decomposition algorithm for the convex RAP–NC with integer variables. This algorithm is a non-trivial generalization of that of Vidal et al. (2016), and attains the same complexity of ${\mathcal{O}}(n\log m\log B)$ . Based on a proximity theorem from Moriguchi et al. (2011), we extend this algorithm to solve the continuous problem in ${\mathcal{O}}(n\log m\log\frac{nB}{\epsilon})$ time. These are the best known complexities, to date, for both problem variants. The complexity depends on the magnitude of $\log(B/\epsilon)$ , a dependency which is known to be unavoidable in the arithmetic model of computation for general forms of convex functions (Renegar 1987). Moreover, the algorithm calls only the oracle for the objective function, without need of gradient information, and does not rely on strict convexity or differentiability. Finally, Lipschitz continuity is assumed for convenience in the proofs but is not mandatory, as an alternative proof line based on submodular optimization could be adopted otherwise.

•

For the specific case of quadratic functions, with continuous or integer variables, the method runs in ${\mathcal{O}}(n\log m)$ time, hence extending the short list of quadratic problems known to be solvable in strongly polynomial time. This also resolves an open question from Moriguchi et al. (2011): “It is an open question whether there exist ${\mathcal{O}}(n\log n)$ algorithms for (Nest) with quadratic objective functions”.

•

We present computational experiments that demonstrate the good performance of the method. We compare it with a known algorithm for the linear case and a general-purpose separable convex optimization solver for the convex case, using benchmark instances derived from three families of applications.

•

We finally integrate the proposed algorithm as a projection step in a projected gradient algorithm for the support vector ordinal regression problem, highlighting promising connections with the machine learning literature.

The rest of the paper is organized as follows. In Section 2, we review related works in a wide variety of application domains. In Section 3, we describe the proposed algorithm and prove its correctness. In Section 4, we report our computational experiments, considering linear and convex objectives, as well as support vector ordinal regression problems. Finally, in Section 5, we provide some concluding remarks.

2 Related Literature and Applications

We now review the many applications of the RAP–NC, starting with classical operations research and management science applications and then moving to statistics, machine learning, and telecommunications.

Resource Allocation. The resource allocation problem (Equations 1, 3, and 4) has long been studied as a prototypical problem. The fastest known algorithms (Frederickson et Johnson 1982, Hochbaum 1994) reach a complexity of ${\mathcal{O}}(n\log\frac{B}{n})$ for the integer problem, and can be extended to find an $\epsilon$ -approximate solution of the continuous problem in ${\mathcal{O}}(n\log\frac{B}{\epsilon})$ operations. This complexity is known to be optimal in the algebraic tree model (Hochbaum 1994). Other algorithms have been developed for several generalizations of the RAP in which the constraint set forms a polymatroid. In this context, the greedy algorithm is optimal (Federgruen et Groenevelt 1986), albeit only pseudo-polynomial, and efficient scaling algorithms can be developed (Hochbaum 1994, Moriguchi et al. 2011). The special case of the integer RAP–NC with $a_{i}=-\infty$ , called NESTED problem, can be solved in ${\mathcal{O}}(n\log n\log\frac{B}{n})$ time using a scaling algorithm (Hochbaum 1994) or in ${\mathcal{O}}(n\log m\log\frac{B}{n})$ time using divide-and-conquer principles (Vidal et al. 2016).

Production Planning. The formulation given by Equations (5)–(8) is also encountered in early literature on production planning over time with inventory and production costs (Wagner et Whitin 1958). One of the models most closely related to our work is that of Love (1973), with time-dependent inventory bounds. The general problem with concave costs (economies of scale) and production capacities is known to be NP-hard. The linear or convex model remains polynomial but is more limited in terms of applicability, although convex production costs can occur in the presence of a limited workforce with possible overtime. Two relatively recent articles have proposed polynomial algorithms for the linear problem with time-dependent inventory bounds (Sedeño-Noda et al. 2004, Ahuja et Hochbaum 2008). With upper bounds $x^{\textsc{max}}_{i}$ on the production quantities, and time-dependent inventory capacities $I^{\textsc{max}}_{i}$ , the problem can be stated as:

[TABLE]

Then, expressing the inventory variables as a function of the production quantities, using $I_{i}=K+\sum_{k=1}^{i}(x_{k}-d_{k})$ , reduces this problem to an RAP–NC:

[TABLE]

The objective includes production costs and inventory costs, and the nested constraints model the time-dependent inventory limit. The algorithm of Ahuja et Hochbaum (2008) can solve Equations (9)–(13) in ${\mathcal{O}}(n\log n)$ time via a reduction to a minimum-cost network flow problem. The method was extended to deal with possible backorders. However, this good complexity comes at the price of an advanced dynamic tree data structure (Tarjan et Werneck 2009) that is used to keep track of the inventory capacities.

Workforce Planning. In contrast with the above studies, which involve the production quantities as decision variables, Bellman et al. (1954) study the balancing of workforce capacity (human or technical resources) over a time horizon under hard production constraints. The variable $x_{i}$ now represents the workforce variation at period $i$ , and the nested constraints impose bounds on the minimum and maximum workforce in certain periods, e.g., to satisfy forecast production demand. The overall objective, to be minimized, is a convex separable cost function representing positive costs for positive or negative variations of the workforce.

Vessel Speed Optimization. In an effort to reduce fuel consumption and emissions, shipping companies have adopted slow-steaming practices, which moderate ship speeds to reduce costs. This line of research has led to several recent contributions on ship speed optimization, aiming to optimize the vessel speed $v_{i-1,i}$ over each trip segment of length $\delta_{i-1,i}$ while respecting a time-window constraint $[a_{i},b_{i}]$ at each destination $i$ . Let $f_{i}(v_{i-1,i})$ be convex functions representing the fuel costs per mile, over $(i-1,i)$ , as a function of $v_{i-1,i}$ , and let $t_{i}$ be the arrival time at $i$ . The overall speed optimization problem can be formulated as:

[TABLE]

Recent work has considered a constant fuel-speed trade-off function on each leg, i.e., $f_{i}=f_{j}$ for all $(i,j)$ . An ${\mathcal{O}}(n^{2})$ recursive smoothing algorithm (RSA) was proposed by Norstad et al. (2011) and Hvattum et al. (2013) for this case. However, assuming constant fuel-speed over the complete trip is unrealistic, since fuel consumption depends on many varying factors, such as sea condition, weather, current, water depth, and ship load (Psaraftis et Kontovas 2013, 2014). The model can be improved by dividing the trip into smaller segments and considering different functions $f_{i}\neq f_{j}$ . This more general model falls outside the scope of applicability of RSA.

Let $v_{i}^{\textsc{opt}}$ be the minimum of each function $f_{i}$ . With the change of variables $x_{1}=t_{1}$ and $x_{i}=t_{i}-t_{i-1}$ for $i\geq 2$ , the model can then be reformulated as:

[TABLE]

This model is a RAP–NC with separable convex cost. An efficient algorithm for this problem is critical, since a speed-optimization algorithm is not often used as a stand-alone tool but rather as a subprocedure in a route-planning algorithm (Psaraftis et Kontovas 2014). This subprocedure can be called several million times when embedded in a local search, on subproblems counting a few hundred variables due to the division of each trip into smaller segments with different sea conditions. Finally, the RAP–NC is also appropriate for variants of vehicle routing problems with emission control (Bektas et Laporte 2011, Kramer et al. 2015a, b) as well as a special case of project crashing for a known critical path (Foldes et Soumis 1993).

Stratified Sampling. Consider a population of $N$ units divided into subpopulations (strata) of $N_{1},\dots,N_{n}$ units such that $N_{1}+\dots+N_{n}=N$ . An optimized stratified sampling method aims to determine the sample size $x_{i}\in[0,N_{i}]$ for each stratum, in order to estimate a characteristic of the population while ensuring a maximum variance level $V$ and minimizing the total sampling cost. Each subpopulation may have a different variance $\sigma_{i}$ , so a sampling plan that is proportional to the size of the subpopulations is frequently suboptimal. The following mathematical model for this sampling design problem was proposed by Neyman (1934) and extended by Srikantan (1963), Hartley (1965), Huddleston et al. (1970), Sanathanan (1971), and others:

[TABLE]

This is a classical RAP formulation. Two extensions of this model are noteworthy in our context. Hartley (1965) and Huddleston et al. (1970) considered multipurpose stratified sampling where more than one characteristic is evaluated while ensuring variance bounds. This leads to several constraints of type (26), and thus to a continuous multidimensional knapsack problem. Sanathanan (1971) considered a hierarchy of strata, with variance bounds for the estimates at each level. This situation occurs for example in survey sampling, when one seeks an estimate of a characteristic at both the national level (first-stage stratum) and the regional level (second-stage stratum). When two stages are considered, we obtain the additional constraints:

[TABLE]

where the $S_{i}$ are disjoint sets of strata, i.e., $\bigcup_{i=1}^{m}S_{i}=\{1,\dots,n\}$ and $S_{i}\cap S_{j}=\varnothing$ for all $i,j$ . The inequalities (28) lead to constraints on the disjoint subsets, giving a resource allocation problem with generalized upper bounds (GUB – Hochbaum 1994, Katoh et al. 2013).

Machine Learning. The support vector machine (SVM) is a supervised learning model which, in its most classical form, seeks to separate a set of samples into two classes according to their labels. This problem is modeled as the search for a separating hyperplane between the projection of the two sample classes into a kernel space of higher dimension, in such a way that the classes are divided by a gap that is as wide as possible, and a penalty for misclassified samples is minimized (Cortes et Vapnik 1995).

As a generalization of the SVM, the support vector ordinal regression (SVOR) aims to find $r-1$ parallel hyperplanes so as to separate $r$ ordered classes of samples. As reviewed in Gutierrez et al. (2016), various models and algorithms have been proposed in recent years to fulfill this task. In particular, the SVOR approach with “explicit constraints on thresholds” (SVOREX) of Chu et Keerthi (2007) obtains a good trade-off between training speed and generalization capability. A dual formulation of SVOREX is presented in Equations (29)–(33). $\mathcal{K}$ is the kernel function, corresponding to a dot product in the kernel space, and $n^{j}$ is the number of samples in a class $j\in\{1,\dots,r\}$ . Each dual variable $\alpha_{i}^{j}$ takes a non-null value only when the $i^{\text{th}}$ sample of the $j^{\text{th}}$ class is active in the definition of the $j^{\text{th}}$ hyperplane, for $j\in\{1,\dots,r-1\}$ . Similarly, each dual variable $\alpha_{i}^{*j}$ takes a non-null value only when the $i^{\text{th}}$ sample of the $j^{\text{th}}$ class is active in the definition of the $(j-1)^{\text{th}}$ hyperplane, for $j\in\{2,\dots,r\}$ . Additional constraints and variables $\mu^{j}$ impose an order on the hyperplanes. For the sake of simplicity, the dummy variables $\boldsymbol{\alpha}^{*1}$ , $\boldsymbol{\alpha}^{r}$ , $\mu^{1}$ , and $\mu^{r}$ are defined and should be fixed to zero.

[TABLE]

The last two constraints of Equations (32)–(33) can be reformulated to eliminate the $\boldsymbol{\mu}$ variables, leading to nested constraints on the variables $\boldsymbol{\alpha}$ and $-\boldsymbol{\alpha}^{*}$ :

[TABLE]

Overall, the problem of Equations (29)–(31) and (34)–(35) is a nonseparable convex problem over the same constraint polytope as the RAP–NC. Note that the number of nested constraints, corresponding to the number of classes, is usually much smaller than the number of variables, which is proportional to the total number of samples, and thus $m\ll n$ .

The solutions of this formulation are usually sparse, since only a fraction of the samples (support vectors) define the active constraints and separating hyperplanes. Given this structure and the size of practical applications, modern solution methods rely on decomposition steps, in which a working set of variables is iteratively re-optimized by a method of choice. Such an approach is referred to as block-coordinate descent in Bertsekas et al. (2003). The convergence of the algorithm can be guaranteed by including in the working set the variables that most severely violate the KKT conditions. Chu et Keerthi (2007), in line with the work of Platt (1998), consider a minimal working set with only two variables at each iteration. The advantage is that the subproblem can be solved analytically in this case, the disadvantage is that a large number of working set selections can be needed for convergence, and the KKT condition check and gradient update may become the bottleneck instead of the optimization itself. To better balance the computational effort and reduce the number of decomposition steps, larger working sets could be considered (e.g., as in SVM ${}^{\textsc{light}}$ of Joachims 1999). Still, to be successful, the algorithm must solve each subproblem, here a non-separable RAP–NC, very efficiently. Such an alternative optimization approach will be investigated in Section 4.3.

Portfolio Management. The mean-variance optimization (MVO) model of Markowitz (1952) has been refined over the years to integrate a large variety of constraints. In its most classical form, the model aims to maximize expected return while minimizing a risk measure such as the variance of the return. This problem can be formulated as:

[TABLE]

where the $x_{i}$ variables, $i\in\{1,\dots,n\}$ , represent investments in different assets, $\mu_{i}$ is the expected return of asset $i$ , and $\sigma_{ij}$ the covariance between asset $i$ and $j$ . In this model, Equation (37) is used to normalize the total investment and Equation (38) prevents short-selling. The literature on these models is vast, and we refer to the recent surveys of Kolm et al. (2014) and Mansini et al. (2014) for more thorough descriptions. Two additional constraint families, often used in practical portfolio models, are closely linked with the RAP–NC:

•

Class constraints limit the investment amounts for certain classes of assets or economic sectors. These can result from regulatory requirements, managerial insights, or customer guidelines (see, e.g., Chang et al. 2000 and Anagnostopoulos et Mamanis 2010). The assets may also be ranked into different categories, e.g., based on their risk or ecological impact. Imposing investment bounds at each level leads to the nested constraints of Equation (6).

•

Fixed transaction costs, minimum transaction levels, and cardinality constraints either impose a fixed price or threshold quantity for any investment in an asset, or limit the number of positions on different assets. These constraints usually require the introduction of additional integer variables $y_{i}$ , taking value one if and only if the asset $i$ is included in the portfolio. This leads to quadratic MIPs, for which metaheuristics (Chang et al. 2000, Crama et Schyns 2003) and branch-and-cut methods (Bienstock 1996, Jobst et al. 2001) form the current state-of-the-art. Bienstock (1996) branches on the $y_{i}$ variables and solves a quadratic resource allocation problem, with additional surrogate constraints in the form of Equation (6), at each node of the search tree. Improved algorithms for the RAP–NC can thus also prove helpful as a methodological building block for more complex portfolio optimization algorithms.

Telecommunications. Constrained resource allocation problems also have a variety of applications in telecommunications. Mobile signals, for example, can be emitted in different directions with different power levels, but interference between signals emitted in the same direction reduces the quality of the communication. In this context, a power and direction must be determined for each signal, while respecting service-quality constraints and minimizing transmission costs. As underlined by Viswanath et Anantharam (2002) and Padakandla et Sundaresan (2009), this problem can be formulated as an instance of the RAP–NC. Given the large size of typical applications, the efficiency of the algorithm is of foremost importance.

A similar model arises for power minimization in multiple-input and multiple-output communication systems, as well as in various other applications of optimization to telecommunications (D’Amico et al. 2014). Moreover, the RAP–NC generalizes a family of multilevel water-filling problems, which have been the focus of significant research (Palomar et Fonollosa 2005). Other applications include power management on multimedia devices, discussed by Huang et Wang (2009) and Gerards et al. (2016). As illustrated by these example applications, the RAP–NC is a prototypical model and an elementary building block for various problems. Therefore, a new algorithmic breakthrough can have considerable impact in many contexts.

3 Proposed Methodology

In this section, we first describe the proposed methodology for the case of continuous variables, and then move on to the case with integer variables. We assume that $a_{i}\leq b_{i}$ for $i\in\{1,\dots,n\}$ , otherwise the problem is trivially infeasible. We will use boldface notation for vectors and normal font for scalars. Let $\mathbf{e}^{s}$ be the unit vector such that $e_{s}=1$ and $e_{i}=0$ for $i\neq s$ .

3.1 Continuous RAP–NC

The proposed algorithm for the RAP–NC is a divide-and-conquer approach over the indices of the nested constraints. It can be seen as a generalization of the method of Vidal et al. (2016), with some fundamental differences related to the number and the nature of the subproblems. For each range of indices $(v,w)$ considered during the search, such that $1\leq v\leq w\leq m$ , it solves four subproblems, $\text{RAP--NC}_{v,w}(L,R)$ for $L\in\{a_{v-1},b_{v-1}\}$ and $R\in\{a_{w},b_{w}\}$ , expressed in Equations (39)–(42), obtained by fixing the $(v-1)^{\text{th}}$ and $w^{\text{th}}$ nested constraints to their lower or upper bounds. $M$ is a large number, defined to be larger than the Lipschitz constant of each function $f_{i}$ .

[TABLE]

To solve these problems when $v<w$ , the algorithm relies on known optimal solutions obtained deeper in a recursion over the range $(v,u)$ and $(u+1,w)$ , with $u=\lfloor(v+w)/2\rfloor$ . When $v=w$ (at the bottom of the recursion), the $\text{RAP--NC}_{v,v}(L,R)$ does not contain any nested constraints from Equation (40) and thus reduces to a simple RAP. We will refer to this approach as the monotonic decomposition algorithm, MDA $(v,w)$ . The original RAP–NC is solved by MDA $(1,m)$ , and the maximum depth of the recursion is $\lceil\log m\rceil$ since the binary decomposition is performed over the $m$ nested constraints.

In the formulation given by Equations (39)–(42), note that the bounds $c_{i}\leq x_{i}\leq d_{i}$ are transferred into the objective via an exact L1 penalty function. This is possible since the functions $f_{i}$ satisfy the Lipschitz condition (Theorem 1), and it helps simplify the exposition and proofs.

Theorem 1 (Relaxation–Penalization).

If there exists a solution $x$ of $\text{RAP--NC}_{v,w}(L,R)$ such that $\mathbf{c}\leq\mathbf{x}\leq\mathbf{d}$ , then all optimal solutions of $\smash{\text{RAP--NC}_{v,w}(L,R)}$ satisfy $\mathbf{c}\leq\mathbf{x}\leq\mathbf{d}$ .

Proof.

Assume the existence of an optimal solution $\mathbf{x^{*}}$ of $\text{RAP--NC}_{v,w}(L,R)$ with an index $s\in\{\sigma[v-1]+1,\dots,\sigma[w]\}$ such that $x^{*}_{s}>d_{s}$ , and a solution $\mathbf{x}$ such that $\mathbf{c}\leq\mathbf{x}\leq\mathbf{d}$ . Since $x^{*}_{s}>d_{s}\geq x_{s}$ and $\sum_{k=\sigma[v-1]+1}^{\sigma[w]}x^{*}_{k}=\sum_{k=\sigma[v-1]+1}^{\sigma[w]}x_{k}=R-L$ , either

[TABLE]

In the first case, define $t=\min\{i\ |\ i>s\text{ and }\sum_{k=\sigma[v-1]+1}^{i}x^{*}_{k}\leq\sum_{k=\sigma[v-1]+1}^{i}x_{k}\}$ . Observe that $x^{*}_{t}<x_{t}$ and thus $d_{t}-x^{*}_{t}>0$ . Moreover, there exists $\Delta>0$ such that, for each $j$ such that $\sigma[j]\in\{s,\dots,t-1\},$ $\smash{\sum_{k=\sigma[v-1]+1}^{\sigma[j]}x^{*}_{k}-\Delta}>\smash{\sum_{k=\sigma[v-1]+1}^{\sigma[j]}x_{k}\geq a_{i}}$ . Defining $\Delta^{\prime}=\min\{\Delta,d_{t}-x^{*}_{t},x^{*}_{s}-d_{s}\}$ , the solution $\mathbf{x^{\prime*}}=\mathbf{x^{*}}+\Delta^{\prime}(\mathbf{e}^{t}-\mathbf{e}^{s})$ is feasible and such that $\bar{f}(\mathbf{x^{\prime*}})=\bar{f}(\mathbf{x^{*}})+\bar{f}_{t}(x^{*}_{t}+\Delta^{\prime})-\bar{f}_{t}(x^{*}_{t})+\bar{f}_{s}(x^{*}_{s}-\Delta^{\prime})-\bar{f}_{s}(x^{*}_{s}).$ Due to the Lipschitz condition, we have $\bar{f}_{t}(x^{*}_{t}+\Delta^{\prime})-\bar{f}_{t}(x^{*}_{t})<M\Delta^{\prime}$ . Moreover, $M\Delta^{\prime}=\bar{f}_{s}(x^{*}_{s})-\bar{f}_{s}(x^{*}_{s}-\Delta^{\prime})$ and thus $\bar{f}(\mathbf{x^{\prime*}})<\bar{f}(\mathbf{x^{*}})$ , contradicting the optimality of $\mathbf{x^{*}}$ .

The second case of Equation (43) is analogous.∎

The main challenge of the MDA is now to exploit the information gathered at deeper steps of the recursion to solve each RAP–NC efficiently. For this purpose, we introduce Theorem 2, which expresses a monotonicity property of the optimal solutions as a function of the resource bound $R$ . As shown subsequently in Theorem 3, this result allows to generate tighter bounds on the variables, which supersede the nested constraints of the RAP–NC and allow to solve all subproblems (at all recursion levels) as simple RAPs.

Theorem 2 (Monotonicity).

Consider three bounds $R^{\downarrow}\leq R\leq R^{\uparrow}$ . If $\mathbf{x^{\downarrow}}$ is an optimal solution of $\text{RAP--NC}_{v,w}(L,R^{\downarrow})$ and $\mathbf{x^{\uparrow}}$ is an optimal solution of $\text{RAP--NC}_{v,w}(L,R^{\uparrow})$ such that $\mathbf{x^{\downarrow}}\leq\mathbf{x^{\uparrow}}$ , then there exists an optimal solution $\mathbf{x^{*}}$ of $\text{RAP--NC}_{v,w}(L,R)$ such that $\mathbf{x^{\downarrow}}\leq\mathbf{x^{*}}\leq\mathbf{x^{\uparrow}}$ .

Proof.

Define $\bar{a}_{i}=a_{i}-L$ and $\bar{b}_{i}=b_{i}-L$ for $i\in\{v,\dots,w-1\}$ as well as $\bar{a}_{w}=\bar{b}_{w}=R-L$ . By the KKT conditions (in the presence of a convex objective over a set of linear constraints), if $\mathbf{x}$ is an optimal solution of $\text{RAP--NC}_{v,w}(L,R)$ , then there exist dual multipliers $(\boldsymbol{\kappa},\boldsymbol{\lambda})$ such that:

[TABLE]

Note the appearance of the subgradients $\partial\bar{f}_{i}$ in Equation (44), since the functions $\bar{f}_{i}$ are not necessarily differentiable. Let $(\boldsymbol{\kappa}^{\uparrow},\boldsymbol{\lambda}^{\uparrow},\boldsymbol{\Phi}^{\uparrow})$ be a set of multipliers associated with the optimal solution $\mathbf{x^{\uparrow}}$ of $\text{RAP--NC}_{v,w}(L,R^{\uparrow})$ , and $\mathbf{x}$ be an optimal solution of $\text{RAP--NC}_{v,w}(L,R)$ . Define $\mathcal{S}^{+}_{\mathbf{x}}=\{i\ |\ x_{i}>x^{\uparrow}_{i}\}$ , $\mathcal{S}^{-}_{\mathbf{x}}=\{i\ |\ x_{i}<x^{\downarrow}_{i}\}$ , and $\mathcal{S}_{\mathbf{x}}=\{i\ |\ x^{\downarrow}_{i}\leq x_{i}\leq x^{\uparrow}_{i}\}$ . We will present a construct that generates a sequence of solutions $(\mathbf{x}^{\textsc{k}})$ , starting from $\mathbf{x}^{0}=\mathbf{x}$ , such that $\big{|}\mathcal{S}^{+}_{\mathbf{x}^{\textsc{k}+1}}\big{|}<\big{|}\mathcal{S}^{+}_{\mathbf{x}^{\textsc{k}}}\big{|}$ and $\big{|}\mathcal{S}^{-}_{\mathbf{x}^{\textsc{k}+1}}\big{|}\leq\big{|}\mathcal{S}^{-}_{\mathbf{x}^{\textsc{k}}}\big{|}$ as long as $\big{|}\mathcal{S}^{+}_{\mathbf{x}^{\textsc{k}}}\big{|}>0$ , leading by recurrence to a solution $\mathbf{\bar{x}}$ such that $\mathbf{\bar{x}}\leq\mathbf{x}^{\uparrow}$ .

If $\big{|}\mathcal{S}^{+}_{\mathbf{x}^{\textsc{k}}}\big{|}>0$ , then there exists $s\in\{\sigma[v-1]+1,\dots,\sigma[w]\}$ such that $x^{\uparrow}_{s}<x^{\textsc{k}}_{s}$ . Let $r$ be the greatest index in $\{\sigma[v-1]+1,\dots,s\}$ such that $\sum_{k=\sigma[v-1]+1}^{r-1}x^{\textsc{k}}_{k}\geq\sum_{k=\sigma[v-1]+1}^{r-1}x^{\uparrow}_{k}$ , and let $t$ be the smallest index in $\{s,\dots,\sigma[w]\}$ such that $\smash{\sum_{k=\sigma[v-1]+1}^{t}x^{\textsc{k}}_{k}\leq\sum_{k=\sigma[v-1]+1}^{t}x^{\uparrow}_{k}}$ .

Since $\smash{R^{\uparrow}-L=\sum_{i=\sigma[v-1]+1}^{\sigma[w]}x^{\uparrow}_{i}\geq\sum_{i=\sigma[v-1]+1}^{\sigma[w]}x^{\textsc{k}}_{i}=R-L}$ , and by the definition of $r$ and $t$ , it follows that $\smash{\sum_{i=r}^{t}x^{\uparrow}_{i}\geq\sum_{i=r}^{t}x^{\textsc{k}}_{i}}$ . Moreover, $r<s\Rightarrow x^{\textsc{k}}_{r}<x^{\uparrow}_{r}$ , and $s<t\Rightarrow x^{\textsc{k}}_{t}<x^{\uparrow}_{t}$ . Finally, note that $r=s=t$ (jointly) is impossible.

•

When $r<s$ , the following statements are valid:

For each $j$ such that $\sigma[j]\in\{r,\dots,s-1\}$ , $\smash{\bar{a}_{j}\leq\sum_{k=\sigma[v-1]+1}^{\sigma[j]}x^{\textsc{k}}_{k}<\sum_{k=\sigma[v-1]+1}^{\sigma[j]}x^{\uparrow}_{k}\leq\bar{b}_{j}}$ (by the definition of $r$ ) and thus $\kappa^{\uparrow}_{j}=\lambda^{\textsc{k}}_{j}=0$ . As a consequence, $\Phi^{\textsc{k}}_{i}\geq\Phi^{\textsc{k}}_{i+1}$ and $\Phi^{\uparrow}_{i}\leq\Phi^{\uparrow}_{i+1}$ for $i\in\{r,\dots,s-1\}$ . The functions $\bar{f}_{i}$ are convex, and thus their (Clarke) subgradients are monotone (Rockafellar 1970), i.e., $\{x^{\uparrow}_{s}<x^{\textsc{k}}_{s},\Phi^{\uparrow}_{s}\in\partial\bar{f}_{s}(x^{\uparrow}_{s}),\Phi^{\textsc{k}}_{s}\in\partial\bar{f}_{s}(x^{\textsc{k}}_{s})\}\Rightarrow\Phi^{\uparrow}_{s}\leq\Phi^{\textsc{k}}_{s}$ . Similarly, we have $\{x^{\uparrow}_{r}>x^{\textsc{k}}_{r},\Phi^{\uparrow}_{r}\in\partial\bar{f}_{r}(x^{\uparrow}_{r}),\Phi^{\textsc{k}}_{r}\in\partial\bar{f}_{r}(x^{\textsc{k}}_{r})\}\Rightarrow\Phi^{\uparrow}_{r}\geq\Phi^{\textsc{k}}_{r}$ . Combining these relations leads to

[TABLE]

and thus there exists $\Psi\in\Re$ such that $\Phi^{\uparrow}_{i}=\Phi^{\textsc{k}}_{i}=\Psi$ for $i\in\{r,\dots,s\}$ .

•

$\bullet$ When $s<t$ , the following statements are valid:

For each $j$ such that $\sigma[j]\in\{s,\dots,t-1\}$ , $\smash{\bar{a}_{j}\leq\sum_{k=\sigma[v-1]+1}^{\sigma[j]}x^{\uparrow}_{k}<\sum_{k=\sigma[v-1]+1}^{\sigma[j]}x^{\textsc{k}}_{k}\leq\bar{b}_{j}}$ (by the definition of $t$ ) and thus $\lambda^{\uparrow}_{j}=\kappa^{\textsc{k}}_{j}=0$ . As a consequence, $\Phi^{\textsc{k}}_{i}\leq\Phi^{\textsc{k}}_{i+1}$ and $\Phi^{\uparrow}_{i}\geq\Phi^{\uparrow}_{i+1}$ for $i\in\{s,\dots,t-1\}$ . Furthermore, as before, $x^{\uparrow}_{s}<x^{\textsc{k}}_{s}$ and $x^{\uparrow}_{t}>x^{\textsc{k}}_{t}$ , and thus $\Phi^{\uparrow}_{s}\leq\Phi^{\textsc{k}}_{s}$ and $\Phi^{\uparrow}_{t}\geq\Phi^{\textsc{k}}_{t}$ . Combining these relations leads to

[TABLE]

and thus there exists $\Psi\in\Re$ such that $\Phi^{\uparrow}_{i}=\Phi^{\textsc{k}}_{i}=\Psi$ for $i\in\{s,\dots,t\}$ .

Overall, $\Phi^{\uparrow}_{i}=\Phi^{\textsc{k}}_{i}=\Psi$ for $i\in\{r,\dots,t\}$ , and thus $\Psi\in\partial\bar{f}_{i}(x^{\textsc{k}}_{i})\cap\partial\bar{f}_{i}(x^{\uparrow}_{i})$ for $i\in\{r,\dots,t\}$ . Define $x_{i}^{\textsc{min}}=\min\{x^{\textsc{k}}_{i},x^{\uparrow}_{i}\}$ and $x_{i}^{\textsc{max}}=\max\{x^{\textsc{k}}_{i},x^{\uparrow}_{i}\}$ . This implies that $\partial\bar{f}_{i}(x)=\{\Psi\}$ for $x\in(x_{i}^{\textsc{min}},x_{i}^{\textsc{max}})$ and thus these functions are affine with identical slope: $\bar{f}_{i}(x)=\bar{f}_{i}(x^{\textsc{k}}_{i})+\Psi(x-x^{\textsc{k}}_{i})$ for $x\in[x_{i}^{\textsc{min}},x_{i}^{\textsc{max}}]$ . We can thus transfer value from the variables of the set $\mathcal{S}^{+}=\mathcal{S}^{+}_{\mathbf{x}^{\textsc{k}}}\cap\{r,\dots,t\}$ to those of the set $\overline{\mathcal{S}^{+}}=\{r,\dots,t\}-\mathcal{S}^{+}$ , via $\textsc{Adjust}([r,\dots,t],\mathbf{x}^{\textsc{k}},\mathbf{x^{\uparrow}})$ (Algorithm 1), leading to a feasible solution $\mathbf{x}^{\textsc{k}+1}$ with the same cost as $\mathbf{x^{\textsc{k}}}$ , hence optimal, such that

[TABLE]

We observe that $\big{|}\mathcal{S}^{+}_{\mathbf{x}^{\textsc{k}+1}}\big{|}<\big{|}\mathcal{S}^{+}_{\mathbf{x}^{\textsc{k}}}\big{|}$ , moreover $\big{|}\mathcal{S}^{-}_{\mathbf{x}^{\textsc{k}+1}}\big{|}\leq\big{|}\mathcal{S}^{-}_{\mathbf{x}^{\textsc{k}}}\big{|}$ . By recurrence, repeating the previous transformation leads to a solution $\mathbf{\bar{x}}$ such that $\mathcal{S}^{+}_{\mathbf{\bar{x}}}=\varnothing$ . A similar principle can then be applied to generate a sequence of solutions $(\mathbf{\bar{x}}^{\textsc{k}})$ , starting from $\mathbf{\bar{x}}^{0}=\mathbf{\bar{x}}$ , such that $\big{|}\mathcal{S}^{-}_{\mathbf{\bar{x}}^{\textsc{k}+1}}\big{|}<\big{|}\mathcal{S}^{-}_{\mathbf{\bar{x}}^{\textsc{k}}}\big{|}$ and $\mathcal{S}^{+}_{\mathbf{\bar{x}}^{\textsc{k}+1}}=\varnothing$ as long as $\big{|}\mathcal{S}^{-}_{\mathbf{\bar{x}}^{\textsc{k}}}\big{|}>0$ , leading to an optimal solution $\mathbf{x^{*}}$ such that $\mathbf{x}^{\downarrow}\leq\mathbf{x}^{*}\leq\mathbf{x}^{\uparrow}$ .∎

Theorem 3 (Variable Bounds).

Let $\mathbf{x^{La}}$ , $\mathbf{x^{Lb}}$ , $\mathbf{x^{aR}}$ , and $\mathbf{x^{bR}}$ be optimal solutions of $\text{RAP--NC}_{v,u}(L,a_{u})$ , $\text{RAP--NC}_{v,u}(L,b_{u})$ , $\text{RAP--NC}_{u+1,w}(a_{u},R)$ , and $\text{RAP--NC}_{u+1,w}(b_{u},R)$ , respectively. If $\mathbf{x^{La}}\leq\mathbf{x^{Lb}}$ and $\mathbf{x^{bR}}\leq\mathbf{x^{aR}}$ , then there exists an optimal solution $\mathbf{x^{*}}$ of $\text{RAP--NC}_{v,w}(L,R)$ such that:*

[TABLE]

Proof.

Let $\mathbf{x}$ be an optimal solution of $\text{RAP--NC}_{v,w}(L,R)$ . As such, $(x_{\sigma[v-1]+1},\dots,x_{\sigma[u]})$ and $(x_{\sigma[u]+1},\dots,x_{\sigma[w]})$ must be optimal solutions of $\text{RAP--NC}_{v,u}(L,X)$ and $\text{RAP--NC}_{u+1,w}(X,R)$ with $\smash{X=L+\sum_{i=\sigma[v-1]+1}^{\sigma[u]}x_{i}}$ . Since $a_{u}\leq X\leq b_{u}$ , there exists an optimal solution $\mathbf{x^{*}}$ of $\text{RAP--NC}_{v,u}(L,X)$ such that $x_{i}^{La}\leq x^{*}_{i}\leq x_{i}^{Lb}$ for $i\in\{\sigma[v-1]+1,\dots,\sigma[u]\}$ via Theorem 2. The other inequality is obtained for $i\in\{\sigma[u]+1,\dots,\sigma[w]\}$ with the same argument, after re-indexing the variables downwards from $\sigma[w]$ to $\sigma[u]+1$ .∎

As a consequence of Theorems 2 and 3, the inequalities of Equations (50)–(51) are valid and can be added to the RAP–NC formulation given by Equations (39)–(42). Moreover, we show that these inequalities dominate the nested constraints of Equation (40). Indeed,

[TABLE]

Moreover, Equations (50)–(51) imply that:

[TABLE]

and combining Equation (53) and (54) leads to:

[TABLE]

Therefore, the nested constraints are superseded at each level of the recursion by the variable bounds obtained from the subproblems. The immediate consequence is a problem simplification: without nested constraints, the formulation reduces to a simple RAP given in Equations (56)–(58), which can be efficiently solved by the algorithm of Frederickson et Johnson (1982) or Hochbaum (1994).

[TABLE]

The pseudocode of the overall decomposition approach is summarized in Algorithm 2.

Two final discussions follow.

•

First, observe the occurrence of Algorithm 1 (Adjust function, introduced in the proof of Theorem 2) before setting the RAP bounds. This ${\mathcal{O}}(n)$ time function can only occur when the functions $f_{i}$ are not strictly convex; in these cases, the solutions of the subproblems may not directly satisfy $\mathbf{x^{La}}\leq\mathbf{x^{Lb}}$ and $\mathbf{x^{bR}}\leq\mathbf{x^{aR}}$ because of possible ties between resource-allocation choices. Alternatively, one could also use a stable RAP solver that guarantees that the solution variables increase monotonically with the resource bound.

•

Second, note the occurrence of the L1 penalty function associated with the original variables’ bounds $c_{i}$ and $d_{i}$ in $\bar{f}_{i}(x_{i})$ while $\bar{c}_{i}$ and $\bar{d}_{i}$ are maintained as hard constraints. Indeed, some subproblems (e.g., $\text{RAP--NC}_{v,v+1}(b_{v},a_{v+1})$ when $b_{v}\geq a_{v+1}$ and $\mathbf{c}=0$ ) may not have a solution respecting the bounds $c_{i}$ and $d_{i}$ . On the other hand, the $\bar{c}_{i}$ and $\bar{d}_{i}$ constraints can always be fulfilled, otherwise the original problem would be infeasible, and their validity is essential to guarantee the correctness of the algorithm.

Nevertheless, since efficient RAP algorithms exist for some specific forms of the objective function, e.g., quadratic (Brucker 1984, Ibaraki et Katoh 1988), we wish to avoid explicit penalty terms in the objective. Therefore we note that an optimal solution $\mathbf{x^{*}}$ of $\text{RAP}_{v,w}(L,R,\mathbf{\bar{c}},\mathbf{\bar{d}})$ can be obtained as follows:

[TABLE]

where $c^{\prime}_{i}=\max\{\bar{c}_{i},\min\{c_{i},\bar{d}_{i}\}\}$ , $d^{\prime}_{i}=\min\{\bar{d}_{i},\max\{d_{i},\bar{c}_{i}\}\}$ , and $\mathbf{x}$ is the solution of the same RAP with the hard constraints of Equation (60):

[TABLE]

Thus, the penalty functions for $c_{i}$ and $d_{i}$ are taken into account by a ${\mathcal{O}}(n)$ test during each RAP resolution, and they never appear in the objective. Experimentally, we observe that the subproblems that fall in the first two cases of Equation (59) are solved notably faster, since they do not even require finding the minimum of a convex function.

3.2 Integer Optimization and Proximity

The previous section has considered continuous decision variables and proven the validity of the algorithm when all the subproblems are solved to optimality. Still, this proof is of limited practical utility for bit-complexity computational models, since the solutions of separable convex problems can involve irrational numbers (e.g., $\min f(x)=x^{3}-6x,x\geq 0$ ) which have no finite binary representation. Therefore, assuming that a subproblem is solved to optimality without any assumption on the shape of the functions is impracticable.

For this reason, most articles that present computational complexity results for convex resource allocation and network flow problems rely on the notion of $\epsilon$ -approximate solutions, located in the proximity of a truly optimal but not necessarily representable solution. In a decomposition algorithm such as MDA, proving that the method produces an $\epsilon$ -approximate solution for a given $\epsilon$ would require to control the precision of the algorithm at each level of the recursion, which could be cumbersome. Therefore, we adopt another approach, typically used in scaling algorithms (Hochbaum 1994, Moriguchi et al. 2011), which consists in proving the validity of the algorithm for integer variables, and using a proximity theorem between the integer and continuous solutions. By solving an integer problem scaled by an appropriate factor, and translating back the integer solution into a continuous solution, any desired $\epsilon$ precision can be achieved.

We define the functions $\bar{f}^{\textsc{pl}}_{i}(x)=\bar{f}_{i}(\lfloor x\rfloor)+(x-\lfloor x\rfloor)\times(\bar{f}_{i}(\lceil x\rceil)-\bar{f}_{i}(\lfloor x\rfloor))$ , which correspond to an inner linearization of the objective using as base the set of integer values. We call the linearized problem $\text{RAP--NC}^{\textsc{pl}}_{v,w}(L,R)$ ; it aims to find the minimum of $\bar{f}^{\textsc{pl}}(\mathbf{x})=\sum_{i=\sigma[v-1]+1}^{\sigma[w]}\bar{f}^{\textsc{pl}}_{i}(x_{i})$ subject to Equations (40)–(42). Since $\bar{f}_{i}$ and $\bar{f}^{\textsc{pl}}_{i}$ coincide on the integer domain, the integer $\text{RAP--NC}_{v,w}(L,R)$ and $\text{RAP--NC}^{\textsc{pl}}_{v,w}(L,R)$ have the same set of optimal solutions. Beyond this, there is a close relationship between the solutions of the integer $\text{RAP--NC}^{\textsc{pl}}_{v,w}(L,R)$ and those of its continuous counterpart, as formulated in Theorem 4, allowing us to prove the validity of Algorithm 2 for integer variables (Theorem 5).

Theorem 4 (Reformulation).

Any optimal solution $\mathbf{x^{*}}$ of the integer $\text{RAP--NC}^{\textsc{pl}}_{v,w}(L,R)$ is also an optimal solution of the continuous $\text{RAP--NC}^{\textsc{pl}}_{v,w}(L,R)$ .

Proof.

By contradiction. Suppose that $\mathbf{x^{*}}$ is not an optimal solution of the continuous $\text{RAP--NC}^{\textsc{pl}}_{v,w}(L,R)$ . Hence, there exists $\mathbf{x}$ such that $\bar{f}^{\textsc{pl}}(\mathbf{x})<\bar{f}^{\textsc{pl}}(\mathbf{x^{*}})$ , and the set $\{i\ |\ x_{i}-\lfloor x_{i}\rfloor>0\}$ contains at least two elements since $\smash{\sum_{i=\sigma[v-1]+1}^{\sigma[w]}x_{i}}=R-L\in\mathbb{Z}$ . Let $s$ and $t$ be, respectively, the first and second indices in this set. We know that the functions $f^{\textsc{pl}}_{s}$ and $f^{\textsc{pl}}_{t}$ are linear in $[\lfloor x_{s}\rfloor,\lceil x_{s}\rceil]$ and $[\lfloor x_{t}\rfloor,\lceil x_{t}\rceil]$ , respectively, with slope $\Phi_{s}$ and $\Phi_{t}$ . Observe that the solution

[TABLE]

is feasible and such that $\bar{f}^{\textsc{pl}}(\mathbf{x^{\prime}})\leq\bar{f}^{\textsc{pl}}(\mathbf{x})$ . Also, note that the number of non-integer values of $\mathbf{x^{\prime}}$ has been strictly decreased (by one or two). Repeating this process, we obtain an integer solution $\mathbf{x^{**}}$ such that $\bar{f}^{\textsc{pl}}(\mathbf{x^{**}})\leq\bar{f}^{\textsc{pl}}(\mathbf{x})<\bar{f}^{\textsc{pl}}(\mathbf{x^{*}})$ . This contradicts the original assumption that $\mathbf{x^{*}}$ is an optimal solution of the integer $\text{RAP--NC}^{\textsc{pl}}_{v,w}(L,R)$ .∎

Theorem 5 (Integer variables).

Theorems 1, 2, 3 and Algorithm 2 remain valid for RAP–NCs with integer variables.

Proof.

The mathematical arguments used in these proofs are independent of the continuous or integer nature of the variables. Moreover, the solution transformation of Algorithm 1 preserves the integrality of the variables. The only element that requires continuous variables is the use of the (necessary) KKT conditions in Equations (44)–(47). However, as we have demonstrated in Theorem 4, an optimal solution of the $\text{RAP--NC}^{\textsc{pl}}_{v,w}(L,R)$ with integer variables is also an optimal solution of the continuous $\text{RAP--NC}^{\textsc{pl}}_{v,w}(L,R)$ . Thus, the KKT conditions with functions $f_{i}^{\textsc{pl}}$ are necessary, hence completing the proof.∎

Finally, we exploit a proximity result for the solutions of the continuous and integer $\text{RAP--NC}^{\textsc{pl}}_{v,w}(L,R)$ :

Theorem 6 (Proximity).

For any integer optimal solution $\mathbf{x^{*}}$ of RAP–NC with $n\geq 2$ variables, there is a continuous optimal solution $\mathbf{x}$ such that*

[TABLE]

This theorem allows us to search for an $\epsilon$ -approximate solution of the continuous problem by defining an integer RAP–NC in which all parameters ( $a_{i}$ , $b_{i}$ , $c_{i}$ , $d_{i}$ ) have been scaled by a factor $\lceil n/\epsilon\rceil$ , solving this problem, and transforming back the solution. It constitutes a special case of Theorem 1.3 from Moriguchi et al. (2011), as the RAP–NC can be shown to be a special case of resource allocation problem under submodular constraints. Moreover, without even relying on submodular optimization arguments, this result can also be obtained directly via first-order (KKT) optimality conditions. This alternative proof shares many similarities with that of Theorem 3, and it is made available in Appendix A for the interested reader.

3.3 Computational Complexity

Convex objective. Each call to the main algorithm $\textsc{MDA}(v,w)$ involves a recursive call to $\textsc{{MDA}}(v,u)$ and $\textsc{{MDA}}(u+1,w)$ with $u=\lfloor\frac{v+w}{2}\rfloor$ , as well as

•

the solution of $\textsc{Rap}_{v,w}(L,R,\mathbf{\bar{c}},\mathbf{\bar{d}})$ for $L\in\{a_{v-1},b_{v-1}\}$ and $R\in\{a_{w},b_{w}\}$ ;

•

up to four calls to the Adjust function;

•

a linear number of operations to set the bounds $\bar{c}_{i}$ and $\bar{d}_{i}$ .

The function Adjust uses a number of elementary operations which grows linearly with the number of variables. Moreover, in the presence of integer variables, each RAP subproblem with $n$ variables and bound $B=R-L$ can be solved in ${\mathcal{O}}(n\log\frac{B}{n})$ time using the algorithm of Frederickson et Johnson (1982) or Hochbaum (1994). As a consequence, the number of operations $\Phi(n,m,B)$ of MDA, as a function of the number of variables $n$ and constraints $m$ , is bounded as

[TABLE]

where $K$ and $K^{\prime}$ are constants and $h=1+\lceil\log_{2}m\rceil$ . Thus, $\Phi(n,m,B)\in{\mathcal{O}}(n\log m\log B)$ in the integer case. For the continuous case, after scaling all problem parameters by $\lceil n/\epsilon\rceil$ , the complexity of the algorithm for the search for an $\epsilon-$ approximate solution becomes ${\mathcal{O}}(n\log m\log\frac{nB}{\epsilon})$ .

Quadratic and linear objectives. More efficient RAP solution methods are known for specific forms of objective functions. The quadratic RAP with continuous variables, in particular, can be solved in ${\mathcal{O}}(n)$ time (Brucker 1984). For the integer case, reviewed in Katoh et al. (2013), an ${\mathcal{O}}(n)$ algorithm can be derived from Section 4.6 of Ibaraki et Katoh (1988). Finally, in the linear case, each RAP subproblem can be solved in ${\mathcal{O}}(n)$ time as a weighted median problem (see, e.g., Korte et Vygen 2012). All these cases lead to ${\mathcal{O}}(n\log m)$ algorithms for the corresponding RAP–NC. Note that no transformation or proximity theorem is needed for the continuous quadratic RAP, since the solutions of quadratic problems are representable.

4 Computational Experiments

We perform computational experiments to evaluate the performance of the proposed algorithm in the presence of a linear objective, and for two convex objective functions arising in project crashing and speed optimization applications. For linear problems, we compare with the network flow algorithm of Ahuja et Hochbaum (2008), which achieved the previous best-known complexity of ${\mathcal{O}}(n\log n)$ for the problem; this complexity is slightly improved to ${\mathcal{O}}(n\log m)$ by the proposed MDA. For general convex objectives, no dedicated algorithm is available and we compare with the interior-point-based algorithm of MOSEK v7.1 for separable convex optimization. We finally report experimental analyses to evaluate the potential of this solver within a projected gradient method for the SVOREX problem (Section 2), for ordinal regression. The algorithms are implemented in C++ and executed on a single core of a Xeon 3.07 GHz CPU. For accurate time measurements, any algorithm with a CPU time smaller than one second was executed multiple times in a loop (up to a total time of 10 seconds) to determine the average time of a run.

We generated benchmark instances with a number of variables $n\in\{10,20,50,100,200,\dots,10^{6}\}$ . Overall, $10$ random benchmark instances were produced for each problem size, leading to a total of 16 $\times$ 10 instances with the same number of nested constraints as decision variables ( $n=m$ ). For fine-grained complexity analyses in the case of the linear objective, we also removed random nested constraints to produce an additional set of 13 $\times$ 10 instances with $m=100$ constraints and $n\in\{100,200,500,\dots,10^{6}\}$ . For each instance, we generated the parameters $c_{i}$ and $d_{i}$ for $i\in\{1,\dots,n\}$ from uniform distributions in the range $[0.1,0.5]$ and $[0.5,0.9]$ , respectively. Then, we defined two sequences of values $v_{i}$ and $w_{i}$ , such that $v_{0}=w_{0}=0$ , $v_{i}=v_{i-1}+X^{v}_{i}$ , and $w_{i}=w_{i-1}+X^{w}_{i}$ for $i\in\{1,\dots,n\}$ , where $X^{v}_{i}$ and $X^{w}_{i}$ are random variables drawn from a uniform distribution in the range $[c_{i},d_{i}]$ . Finally, we set $a_{i}=\min\{v_{i},w_{i}\}$ and $b_{i}=\max\{v_{i},w_{i}\}$ for all $i$ . We also selected a random parameter $p_{i}$ in $[0,1]$ to characterize the objective function. We conducted the experiments with four classes of objectives: a linear objective $\sum_{i=1}^{n}p_{i}x_{i}$ , and three convex objectives defined as:

[TABLE]

where the last two objectives are representative of applications in project crashing (Foldes et Soumis 1993) and ship speed optimization (Ronen 1982).

4.1 Linear Objective

We start the experimental analyses with the linear RAP–NC. We will refer to the network-flow-based approach of Ahuja et Hochbaum (2008) as “FLOW” in the text and tables. This method was precisely described in the original article, but no computational experiments or practical implementation were reported, so we had to implement it. The authors suggest the use of a red-black tree to locate the minimum-cost paths and a dynamic tree (Tarjan 1997, Tarjan et Werneck 2009) to manage the capacity constraints. This advanced data structure requires significant implementation effort and can result in high CPU time constants. We thus adopted a simpler structure, a segment tree (Bentley 1977) with lazy propagation, which allows evaluating and updating these capacities with the same complexity of ${\mathcal{O}}(\log n)$ per operation (and possibly a higher speed in practice). The proposed MDA was implemented as in Algorithm 2, solving each linear RAP subproblem in ${\mathcal{O}}(n)$ time as a variant of a weighted median problem (Korte et Vygen 2012).

We executed both algorithms on each instance. The results for the instances with as many nested constraints as decision variables ( $n=m$ ) are reported in Figure 1. To evaluate the growth of the computational effort of the algorithms as a function of problem size, we fitted the computational time as a power law $f(n)=\alpha\times n^{\beta}$ of the number of variables $n$ , via a least-squares regression of an affine function on the log-log graph (left figure). We also display as boxplots the ratio of the computational time of MDA and FLOW (right figure). The same conventions are used to display the results of the experimental analyses with a fixed number of constraints ( $m=100$ ) and increasing number of variables $n$ in Figure 2. Finally, the detailed average computational times for each group of 10 instances are reported in Table 1.

From these experiments, we observe that the computational times of the two methods are very similar in terms of magnitude and growth rate. When $m=n$ , the algorithms have the same theoretical complexity of ${\mathcal{O}}(n\log n)$ , as confirmed by the power law regression, with an observed growth that is close to linear (in $n^{1.16}$ and $n^{1.18}$ ). The FLOW algorithm is on average $1.1\times$ to $1.4\times$ faster than MDA for $n\in[10,100]\cup[10^{4},10^{6}]$ , for instances with the same number of variables and constraints. On the other hand, MDA is on average $2\times$ faster than FLOW when $m$ is fixed and $n$ grows beyond $1000$ . This is due to the difference in computational complexity: ${\mathcal{O}}(n\log m)$ for MDA instead of ${\mathcal{O}}(n\log n)$ . MDA and FLOW solve the largest instances, with up to $n=m=10^{6}$ constraints and variables, in three and two seconds on average, respectively.

Overall, the two algorithms have similar performance for linear objectives, and the CPU differences are small. Since these algorithms are based on drastically different principles, they lead the way to different methodological extensions. The computational complexity of FLOW is tied to its efficient use of a dynamic tree data structure, while the complexity of the MDA stems from its “monotonic” divide-and-conquer strategy. Because of this structure, MDA should be a good choice for re-optimization after a change of a few parameters, as well as for the iterative solution of multiple RAP–NC, e.g., for speed optimization within an algorithm enumerating a large number of similar visit sequences, since it can reuse the solutions of smaller subproblems (see, e.g., Norstad et al. 2011 and Vidal et al. 2014, 2015).

4.2 Separable Convex Objectives

In contrast with the case of linear objective functions, no specialized algorithm has been designed for the RAP–NC with separable convex costs to date. To illustrate the possible gain achieved by the use of a dedicated algorithm rather than a general-purpose solver, we do a simple comparison of the CPU time of MDA to that of the MOSEK v7.1 solver on two sets of instances with objective functions derived from project crashing and speed optimization applications. MOSEK is based on an interior-point method and is a good representative of the current generation of separable convex optimization solvers. We set a time limit of one hour. To simplify the execution of these experiments, we use a binary search over the single dual variable to solve each continuous RAP subproblem (Patriksson 2008) within a precision of $10^{-9}$ . As these approximations may stack up $\log m$ times in the recursion, we obtain an overall good accuracy but do not guarantee $\epsilon$ -proximity in these tests.

The results are reported in Figure 3 and Table 2. In the figure, the power-law regressions are presented only for MDA, since MOSEK does not exhibit polynomial behavior, likely due to the computational effort related to the initialization of the solver for small problems. In contrast, the computational time of MDA grows steadily in ${\mathcal{O}}(n^{1.19})$ at most. This observation is consistent with the theoretical ${\mathcal{O}}(n\log m)$ complexity of the method. Within one hour, MOSEK solved all the instances up to $n=5,000$ decision variables. In contrast, MDA solved all the available instances with up to one million variables. For the largest benchmark instances, the CPU time of the method did not exceed $50$ seconds. As illustrated in the bottom-right subfigure, the ratio of the CPU time of MOSEK and MDA ranges between 16 and 28,000. For all the instances, significant CPU time is saved when using the monotonic divide-and-conquer algorithm instead of a general-purpose solver.

4.3 Non-Separable Convex Objective – Support Vector Ordinal Regression

Our last experimental analysis is concerned with the SVOREX model, presented in Section 2. It is a non-separable convex optimization problem over a special case of the RAP–NC constraint polytope. The current state-of-the-art algorithm for this problem, proposed by Chu et Keerthi (2007), is based on a working-set decomposition. Iteratively, a set of variables is selected to be optimized over, while the others remain fixed. This approach leads to a (non-separable) restricted problem with fewer variables which can be solved to optimality. The authors rely on a minimal working set containing the two variables which most violate the KKT conditions (see Chu et Keerthi 2007, pp. 799–800, for all equations involved).

The advantage of a minimal working set comes from the availability of analytical solutions for the restricted problems. On the other hand, larger working sets can be beneficial in order to reduce the number of iterations until convergence (see, e.g., Joachims 1999). However this would require an efficient method for the resolution of the reduced problems. This is how the proposed RAP–NC solver can provide a meaningful option along this direction. In order to evaluate such a proof of concept, we conduct a simple experiment which consists of generating larger working sets within the approach of Chu et Keerthi (2007) and solving the resulting reduced problems with the help of the RAP–NC algorithm. As these reduced problems are non-separable convex, the RAP–NC algorithm is being used for the projection steps within a projected gradient descent procedure. The overall solution approach is summarized in Algorithm 3, in which $\mathcal{W}$ is the working set, $z$ is the objective function, and $\gamma$ is the fixed step size of the gradient descent.

To obtain a larger working set, we repeatedly select the most-violated sample pair until either reaching the desired size or not finding any remaining violation. In our experiments, we consider working sets of size $n_{\textsc{ws}}\in\{2,4,6,10\}$ , a step size of $\gamma=0.2$ and $n_{\textsc{grad}}=20$ iterations for the projected gradient descent. We use the eight problem instances introduced in Chu et Keerthi (2007), with the same Gaussian kernel, penalty parameter, and guidelines for data preparation (normalizing the input vectors to zero mean and unit variance, and using equal-frequency binning to discretize the target values into five ordinal scales).

Table 3 gives the results of these experiments. The columns report, in turn, the problem instance name, its number of samples $N$ , the dimension $D$ of its feature space, and characteristics of the optimal solutions: the number of variables set to [math] (correct classification), to $C$ (misclassified), and to intermediate values (support vectors). For each working-set size $n_{\textsc{ws}}$ , the total number of working set selections $I_{\textsc{ws}}$ done by the algorithm is also presented, as well as the CPU time in seconds. The fastest algorithm version is underlined for each instance.

As measured in these experiments, the CPU time of the algorithms ranges between 0.018 seconds for the smallest problem instances (with 50 samples and 27 dimensions) and 574.28 seconds for the largest case (6000 samples and 16 dimensions). The size of the working set has a significant impact on the number of iterations of the method and its CPU time.

In all cases, the number of iterations decreases significantly when the size of the working set grows. In terms of CPU time, the fastest results are either achieved with a working set of size two or six (with three instances in each case). We observe that the three instances for which a larger working set contributed to CPU-time reductions are those with higher-dimension feature spaces (dimension 21 to 32). For these instances, using a larger working set helped reduce the CPU time by a factor of $1.1$ to $1.9$ as compared to using a two-samples working set. In comparison, Joachims (1999) reported a speedup of $1.5$ to $2.0$ when using ten-samples working sets for SVM, which can be viewed as a special case of SVOREX with two classes.

To achieve a gain in CPU time, the number of iterations should decrease more than linearly as a function of the working set size. This is due to the effort spent updating the gradient, necessary for the verification of the KKT conditions and the working-set selection, which grows linearly with the product $I_{\textsc{ws}}\times N\times n_{\textsc{ws}}$ (using efficient incremental updates), and which remains a major bottleneck for SVOREX and SVM algorithms (see, e.g., the discussions in Joachims 1999 and Chang et Lin 2011). Usually, instances with a feature space of high dimension exhibit a fast decrease in the number of iterations as a function of the working set size, as their solutions include a larger proportion of variables taking values in $(0,C)$ (support vectors), values which are more quickly reached via simultaneous optimizations of several variables. As such, larger working sets are likely to be more useful in feature spaces of higher dimension.

Moreover, future research avenues concern possible improvements of the algorithm (e.g., using shrinking or a double-loop scheme – Keerthi et al. 2001), adaptive choices of working-set size based on analyses of the structure of the data set and solutions, as well as more advanced selection rules, e.g., based on Zoutendijk’s descent direction (Joachims 1999) or second order information (Fan et al. 2005). These options for improvement are now possible due to the availability of a fast algorithm for the resolution of the restricted problems.

5 Concluding Remarks

In this article, we have highlighted the importance of the RAP–NC, which is a problem connected with a wide range of applications in production and transportation optimization, portfolio management, sampling optimization, telecommunications and machine learning. To solve this problem, we proposed a decomposition algorithm, based on monotonicity principles coupled with divide-and-conquer, leading to new complexity breakthroughs for a variety of objectives (linear, quadratic, and convex), with continuous or integer variables, and to the first known strongly polynomial algorithm for the quadratic integer RAP–NC. In terms of practical performance, the algorithm matches the best dedicated (flow-based) algorithm for the linear case, outperforms general-purpose solvers by several orders of magnitude in the convex case, and opens interesting perspectives of algorithmic improvements for the SVOREX problem in machine learning.

The algorithm can be seen as a generalization of the method of Vidal et al. (2016), with some fundamental differences related to the number and the nature of the subproblems. It is not based on classical greedy steps and scaling, or on flow propagation techniques, often exploited for this problem family. It is an important research question to see how far this decomposition technique can be generalized. In particular, the approach can very likely be extended to the resource allocation problem with a TREE of lower and upper constraints (Hochbaum 1994). Other optimization problems related to PERT (Program Evaluation and Review Technique) may exhibit monotonicity properties as a function of time constraints or budget bounds, and we should investigate how to decompose efficiently their variables and constraints while maintaining a low computational complexity. Similarly, extended formulations involving the intersection of two or more RAP–NC type of constraint polytopes deserve a closer look. These are all open important research directions which can be explored in the near future.

Acknowledgments

This research was partially supported by the National Counsel of Technological and Scientific Development and Fundação de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ) in Brazil, grants 308498/2015-1 and E-26/203.310/2016, and by the Office of Naval Research (ONR), USA, grant N00014-15-1-2083.

Appendix A – Proof of Theorem 6, based on KKT conditions

Proof.

The proof shares many similarities with that of Theorem 2. It exploits the fact that an integer solution $\mathbf{x^{*}}$ of the $\text{RAP--NC}^{\textsc{pl}}$ is also an optimal solution of the continuous problem (Theorem 4) and thus satisfies the KKT conditions of Equations (44)–(47) based on the functions $f_{i}^{\textsc{pl}}$ . We first state two lemmas, that will be used later to link the values of the subderivatives of $f_{i}$ and $f_{i}^{\textsc{pl}}$ .

Lemma 1.

Consider $y\in\mathbb{R}$ and $x\in\mathbb{Z}$ such that $y+1\leq x$ , and a convex function $f$ . If $\phi_{y}\in\partial f(y)$ and $\phi_{x}\in\partial f^{\textsc{pl}}(x)$ , then $\phi_{y}\leq\phi_{x}$ .

Proof of Lemma 1. By definition, $f(x)-f(y)\geq\phi_{y}(x-y)$ and $f^{\textsc{pl}}(x-1)-f^{\textsc{pl}}(x)\geq\phi_{x}(x-1-x)=-\phi_{x}$ . Moreover, $y\leq x-1\leq x$ , $f$ is convex, and $f$ coincides with $f^{\textsc{pl}}$ at $x$ and $x-1$ , so

[TABLE]

Lemma 2.

Consider $y\in\mathbb{Z}$ and $x\in\mathbb{R}$ such that $y+1\leq x$ , and a convex function $f$ . If $\phi_{y}\in\partial f^{\textsc{pl}}(y)$ and $\phi_{x}\in\partial f(x)$ , then $\phi_{y}\leq\phi_{x}$ .

Proof of Lemma 2. By definition, $f^{\textsc{pl}}(y+1)-f^{\textsc{pl}}(y)\geq\phi_{y}$ and $f(y)-f(x)\geq\phi_{x}(y-x)$ . Moreover, $y\leq y+1\leq x$ , $f$ is convex and $f$ coincides with $f^{\textsc{pl}}$ at $y$ and $y+1$ , so

[TABLE]

The main proof follows. Let $\mathbf{x}$ be an optimal continuous solution of the RAP–NC. If Equation (62) is satisfied, then the proof is complete; otherwise there exists $s\in\{1,\dots,n\}$ such that $|x_{s}-x^{*}_{s}|\geq n-1$ . We consider here the case where $x_{s}\geq x^{*}_{s}+n-1$ , the other case being symmetric. Let $r$ be the greatest index in $\{1,\dots,s\}$ such that $\sum_{k=1}^{r-1}x_{k}\geq\sum_{k=1}^{r-1}x^{*}_{k}$ , and $t$ be the smallest index in $\{s,\dots,n\}$ such that $\smash{\sum_{k=1}^{t}x_{k}\leq\sum_{k=1}^{t}x^{*}_{k}}$ . By the definition of $r$ and $t$ , it follows that $\smash{\sum_{i=r}^{t}x^{*}_{i}\geq\sum_{i=r}^{t}x_{i}}$ , and thus $\smash{\sum_{i\in\{r,\dots,t\}-s}x^{*}_{i}\geq\sum_{i\in\{r,\dots,t\}-s}x_{i}+n-1}$ . Since $|\{r,\dots,t\}-s|\leq n-1$ , there exists $u\in\{r,\dots,t\}-s$ such that $x^{*}_{u}\geq x_{u}+1$ .

Two cases can arise:

If $u<s$ , for each $j$ such that $\sigma[j]\in\{u,\dots,s-1\}$ , $\bar{a}_{j}\leq\sum_{k=1}^{\sigma[j]}x_{k}<\sum_{k=1}^{\sigma[j]}x^{*}_{k}\leq\bar{b}_{j}$ and thus $\kappa^{*}_{j}=\lambda_{j}=0$ . As a consequence, $\Phi_{i}\geq\Phi_{i+1}$ and $\Phi^{*}_{i}\leq\Phi^{*}_{i+1}$ for $i\in\{u,\dots,s-1\}$ .

If $u>s$ , for each $j$ such that $\sigma[j]\in\{s,\dots,u-1\}$ , $\smash{\bar{a}_{j}\leq\sum_{k=1}^{\sigma[j]}x^{*}_{k}<\sum_{k=1}^{\sigma[j]}x_{k}\leq\bar{b}_{j}}$ and thus $\lambda^{*}_{j}=\kappa_{j}=0$ . As a consequence, $\Phi_{i}\leq\Phi_{i+1}$ and $\Phi^{*}_{i}\geq\Phi^{*}_{i+1}$ for $i\in\{s,\dots,u-1\}$ .

Moreover, $\{x^{*}_{s}+n-1\leq x_{s},\Phi^{*}_{s}\in\partial f_{s}(x^{*}_{s}),\Phi_{s}\in\partial f^{\textsc{pl}}_{s}(x_{s})\}\Rightarrow\Phi^{*}_{s}\leq\Phi_{s}$ (Lemma 1), and $\{x_{u}+1\leq x^{*}_{u},\Phi_{u}\in\partial f^{\textsc{pl}}_{u}(x_{u}),\Phi^{*}_{u}\in\partial f_{u}(x^{*}_{u})\}\Rightarrow\Phi_{u}\leq\Phi^{*}_{u}$ (Lemma 2). Combining all the relations leads to $\Phi_{s}\leq\Phi_{u}\leq\Phi^{*}_{u}\leq\Phi^{*}_{s}\leq\Phi_{s},$ and thus there exists $\Psi\in\Re$ such that $\Phi^{*}_{i}=\Phi_{i}=\Psi$ for $i\in\{u,\dots,s\}$ if $u<s$ (or $i\in\{s,\dots,u\}$ if $s<u$ ). As in the proof of Theorem 2, this implies that the functions $f_{s}$ and $f_{u}$ are affine with slope $\Psi$ over $[x_{s}^{\textsc{min}},x_{s}^{\textsc{max}}]$ and $[x_{u}^{\textsc{min}},x_{u}^{\textsc{max}}]$ , respectively, where $x_{i}^{\textsc{min}}=\min\{x_{i},x^{*}_{i}\}$ . Observe that the new solution $\mathbf{x^{\prime}}=\mathbf{x}-\mathbf{e}^{s}+\mathbf{e}^{u}$ is a feasible solution with the same cost as $\mathbf{x}$ , hence optimal. Moreover, we note that $\sum_{i=1}^{n}\max\{|x^{\prime}_{i}-x^{*}_{i}|-(n-1),0\}\leq\sum_{i=1}^{n}\max\{|x_{i}-x^{*}_{i}|-(n-1),0\}-1$ and/or $\sum_{i=1}^{n}\mathbbm{1}\{|x^{\prime}_{i}-x^{*}_{i}|\leq(n-1)\}\leq\sum_{i=1}^{n}\mathbbm{1}\{|x_{i}-x^{*}_{i}|\leq(n-1)\}-1$ , where $\mathbbm{1}(p)=1$ if and only if $p$ is true. Repeating this process leads, in a finite number of steps, to a solution $\mathbf{x^{\prime\prime}}$ such that $|x^{\prime\prime}_{i}-x^{*}_{i}|<n-1,\text{ for }i\in\{1,\dots,n\}$ .∎

Bibliography60

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Ahuja et Hochbaum (2008) Ahuja R, Hochbaum D (2008) Technical note – Solving linear cost dynamic lot-sizing problems in O(n log n) time. Operations Research 56(1):255–261.
2Anagnostopoulos et Mamanis (2010) Anagnostopoulos K, Mamanis G (2010) A portfolio optimization model with three objectives and discrete variables. Computers & Operations Research 37(7):1285–1297.
3Bektas et Laporte (2011) Bektas T, Laporte G (2011) The pollution-routing problem. Transportation Research Part B: Methodological 45(8):1232–1250.
4Bellman et al. (1954) Bellman R, Glicksberg I, Gross O (1954) The theory of dynamic programming as applied to a smoothing problem. Journal of the Society for Industrial and Applied Mathematics 2(2):82–88.
5Bentley (1977) Bentley J (1977) Solutions to Klee’s rectangle problems. Technical report, Carnegie-Mellon University, Pittsburgh PA.
6Bertsekas et al. (2003) Bertsekas D, Nedi A, Ozdaglar A (2003) Nonlinear Programming (Nashua, NH: Athena Scientific).
7Bienstock (1996) Bienstock D (1996) Computational study of a family of mixed-integer quadratic programming problems. Mathematical Programming 74(2):121–140.
8Brucker (1984) Brucker P (1984) An O(n) algorithm for quadratic knapsack problems. Operations Research Letters 3(3):163–166.

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

Taxonomy

Separable Convex Optimization with Nested Lower and Upper Constraints

1 Introduction

2 Related Literature and Applications

3 Proposed Methodology

3.1 Continuous RAP–NC

Theorem 1** (Relaxation–Penalization).**

Proof.

Theorem 2** (Monotonicity).**

Proof.

Theorem 3** (Variable Bounds).**

Proof.

3.2 Integer Optimization and Proximity

Theorem 4** (Reformulation).**

Proof.

Theorem 5** (Integer variables).**

Proof.

Theorem 6** (Proximity).**

3.3 Computational Complexity

4 Computational Experiments

4.1 Linear Objective

4.2 Separable Convex Objectives

4.3 Non-Separable Convex Objective – Support Vector Ordinal Regression

5 Concluding Remarks

Acknowledgments

Appendix A – Proof of Theorem 6, based on KKT conditions

Proof.

Lemma 1**.**

Lemma 2**.**

Theorem 1 (Relaxation–Penalization).

Theorem 2 (Monotonicity).

Theorem 3 (Variable Bounds).

Theorem 4 (Reformulation).

Theorem 5 (Integer variables).

Theorem 6 (Proximity).

Lemma 1.

Lemma 2.