Consider a vector of variables \(\color{darkred}x_i\). We want to impose the constraint: \[\text{sum of $K$ largest $\color{darkred}x_i$} \le 0.5\sum_i \color{darkred}x_i\] I.e. the sum of the largest \(K\) values are less than half of the total. Sometimes this is expressed mathematically as \[\sum_{i=1}^K \color{darkred}x_{[i]} \le 0.5\sum_i \color{darkred}x_i\] where \(\color{darkred}x_{[i]}\) is an ordered version of \(\color{darkred}x_i\) with \(\color{darkred}x_{[1]}\ge\color{darkred}x_{[2]} \ge \color{darkred}x_{[3]} \dots \). There are several formulations for this constraint, all of them interesting, I think.
The first observation is that, if there are additional constraints, \(\color{darkred}x_i=0\) is feasible. If we exclude this solution, e.g. by a constraint like \[\sum_i \color{darkred}x_i=1\] we need the vector \(\color{darkred}x_i\) to be of a certain length. Obviously, we need at least \(K\) elements to form the \(K\) largest ones. But for those \(K\) largest not to dominate, we need at least \(2K+1\) elements. Only then the \(K\) largest can be less than half of the total.
In my examples, I will use this constraint in a small portfolio problem: \[\begin{align} \min\>&\color{darkred}x'\color{darkblue}Q\color{darkred}x\\ & \sum_i \color{darkred}x_i=1\\&\sum_i \color{darkblue}\mu_i \color{darkred}x_i\ge \color{darkblue}R\\ & \color{darkred}x_i\ge 0\end{align}\] As we know \(\sum_i \color{darkred}x_i=1\), the right-hand side of the sum of \(K\) largest constraint will be constant: \[\text{sum of $K$ largest $\color{darkred}x_i$} \le 0.5\] however, all the methods also work when the right-hand side contains a variable.
As a baseline, we have a random data set with \(N=12\) assets. The min-risk solution is:
---- 57 ***baseline***
VARIABLE z.L = 0.672 objective
---- 57 VARIABLE x.L fractions
i1 0.301, i2 0.138, i6 0.006, i7 0.081, i8 0.222, i9 0.077, i10 0.066, i12 0.109
The sum of the \(K=3\) largest allocations is 0.301+0.222+0.138=0.661. We are violating our special constraint.
All possible combinations
The first formulation is based on the observation: bound the \(K\) largest values by 0.5 is the same as bounding all possible \(K\) combinations of indices \(i\) by 0.5. For example, if we denote the length of \(\color{darkred}x\) by \(N=12\), and use \(K=3\), then we can form the constraints: \[\color{darkred}x_i+\color{darkred}x_j+\color{darkred}x_k \le 0.5\] for all 3-combinations \((i,j,k)\). There are \[{N \choose K} = {12 \choose 3} = 220\] of them:
---- 68 SET combinations3 all possible 3-combinations
i1 .i2 .i3 , i1 .i2 .i4 , i1 .i2 .i5 , i1 .i2 .i6 , i1 .i2 .i7 , i1 .i2 .i8
i1 .i2 .i9 , i1 .i2 .i10, i1 .i2 .i11, i1 .i2 .i12, i1 .i3 .i4 , i1 .i3 .i5
i1 .i3 .i6 , i1 .i3 .i7 , i1 .i3 .i8 , i1 .i3 .i9 , i1 .i3 .i10, i1 .i3 .i11
i1 .i3 .i12, i1 .i4 .i5 , i1 .i4 .i6 , i1 .i4 .i7 , i1 .i4 .i8 , i1 .i4 .i9
i1 .i4 .i10, i1 .i4 .i11, i1 .i4 .i12, i1 .i5 .i6 , i1 .i5 .i7 , i1 .i5 .i8
i1 .i5 .i9 , i1 .i5 .i10, i1 .i5 .i11, i1 .i5 .i12, i1 .i6 .i7 , i1 .i6 .i8
i1 .i6 .i9 , i1 .i6 .i10, i1 .i6 .i11, i1 .i6 .i12, i1 .i7 .i8 , i1 .i7 .i9
i1 .i7 .i10, i1 .i7 .i11, i1 .i7 .i12, i1 .i8 .i9 , i1 .i8 .i10, i1 .i8 .i11
i1 .i8 .i12, i1 .i9 .i10, i1 .i9 .i11, i1 .i9 .i12, i1 .i10.i11, i1 .i10.i12
i1 .i11.i12, i2 .i3 .i4 , i2 .i3 .i5 , i2 .i3 .i6 , i2 .i3 .i7 , i2 .i3 .i8
i2 .i3 .i9 , i2 .i3 .i10, i2 .i3 .i11, i2 .i3 .i12, i2 .i4 .i5 , i2 .i4 .i6
i2 .i4 .i7 , i2 .i4 .i8 , i2 .i4 .i9 , i2 .i4 .i10, i2 .i4 .i11, i2 .i4 .i12
i2 .i5 .i6 , i2 .i5 .i7 , i2 .i5 .i8 , i2 .i5 .i9 , i2 .i5 .i10, i2 .i5 .i11
i2 .i5 .i12, i2 .i6 .i7 , i2 .i6 .i8 , i2 .i6 .i9 , i2 .i6 .i10, i2 .i6 .i11
i2 .i6 .i12, i2 .i7 .i8 , i2 .i7 .i9 , i2 .i7 .i10, i2 .i7 .i11, i2 .i7 .i12
i2 .i8 .i9 , i2 .i8 .i10, i2 .i8 .i11, i2 .i8 .i12, i2 .i9 .i10, i2 .i9 .i11
i2 .i9 .i12, i2 .i10.i11, i2 .i10.i12, i2 .i11.i12, i3 .i4 .i5 , i3 .i4 .i6
i3 .i4 .i7 , i3 .i4 .i8 , i3 .i4 .i9 , i3 .i4 .i10, i3 .i4 .i11, i3 .i4 .i12
i3 .i5 .i6 , i3 .i5 .i7 , i3 .i5 .i8 , i3 .i5 .i9 , i3 .i5 .i10, i3 .i5 .i11
i3 .i5 .i12, i3 .i6 .i7 , i3 .i6 .i8 , i3 .i6 .i9 , i3 .i6 .i10, i3 .i6 .i11
i3 .i6 .i12, i3 .i7 .i8 , i3 .i7 .i9 , i3 .i7 .i10, i3 .i7 .i11, i3 .i7 .i12
i3 .i8 .i9 , i3 .i8 .i10, i3 .i8 .i11, i3 .i8 .i12, i3 .i9 .i10, i3 .i9 .i11
i3 .i9 .i12, i3 .i10.i11, i3 .i10.i12, i3 .i11.i12, i4 .i5 .i6 , i4 .i5 .i7
i4 .i5 .i8 , i4 .i5 .i9 , i4 .i5 .i10, i4 .i5 .i11, i4 .i5 .i12, i4 .i6 .i7
i4 .i6 .i8 , i4 .i6 .i9 , i4 .i6 .i10, i4 .i6 .i11, i4 .i6 .i12, i4 .i7 .i8
i4 .i7 .i9 , i4 .i7 .i10, i4 .i7 .i11, i4 .i7 .i12, i4 .i8 .i9 , i4 .i8 .i10
i4 .i8 .i11, i4 .i8 .i12, i4 .i9 .i10, i4 .i9 .i11, i4 .i9 .i12, i4 .i10.i11
i4 .i10.i12, i4 .i11.i12, i5 .i6 .i7 , i5 .i6 .i8 , i5 .i6 .i9 , i5 .i6 .i10
i5 .i6 .i11, i5 .i6 .i12, i5 .i7 .i8 , i5 .i7 .i9 , i5 .i7 .i10, i5 .i7 .i11
i5 .i7 .i12, i5 .i8 .i9 , i5 .i8 .i10, i5 .i8 .i11, i5 .i8 .i12, i5 .i9 .i10
i5 .i9 .i11, i5 .i9 .i12, i5 .i10.i11, i5 .i10.i12, i5 .i11.i12, i6 .i7 .i8
i6 .i7 .i9 , i6 .i7 .i10, i6 .i7 .i11, i6 .i7 .i12, i6 .i8 .i9 , i6 .i8 .i10
i6 .i8 .i11, i6 .i8 .i12, i6 .i9 .i10, i6 .i9 .i11, i6 .i9 .i12, i6 .i10.i11
i6 .i10.i12, i6 .i11.i12, i7 .i8 .i9 , i7 .i8 .i10, i7 .i8 .i11, i7 .i8 .i12
i7 .i9 .i10, i7 .i9 .i11, i7 .i9 .i12, i7 .i10.i11, i7 .i10.i12, i7 .i11.i12
i8 .i9 .i10, i8 .i9 .i11, i8 .i9 .i12, i8 .i10.i11, i8 .i10.i12, i8 .i11.i12
i9 .i10.i11, i9 .i10.i12, i9 .i11.i12, i10.i11.i12
A disadvantage of writing something like: \(\color{darkred}x_i+\color{darkred}x_j+\color{darkred}x_k \le 0.5\) is that we hard-coded \(K=3\). If we want to make it more adaptive, we need to form a 2d set \(\color{darkblue}{\mathit{comb}}_{n,i}\), with \(n\) the combination id, and \(i\) the indices to include in this combination. This will look like:
---- 113 SET comb combinations
i1 i2 i3 i4 i5 i6 i7 i8 i9
n1 YES YES YES
n2 YES YES YES
n3 YES YES YES
n4 YES YES YES
n5 YES YES YES
n6 YES YES YES
n7 YES YES YES
n8 YES YES
n9 YES YES
n10 YES YES
n11 YES YES YES
n12 YES YES YES
n13 YES YES YES
n14 YES YES YES
n15 YES YES YES
n16 YES YES YES
n17 YES YES
n18 YES YES
n19 YES YES
n20 YES YES YES
. . .
Note that this is a partial view. We can easily generate these combinations with a bit of Python code. itertools.combinations is our friend here. With this set we can write: \[\sum_{i|\color{darkblue}{\mathit{comb}}_{n,i}} \color{darkred}x_{i} \le 0.5 \>\>\forall n\] The structure does not change if we use a different value for \(K\).
When we run one of these models, we see:
---- 75 ***all combinations***
VARIABLE z.L = 0.760 objective
---- 75 VARIABLE x.L fractions
i1 0.160, i2 0.139, i3 0.001, i7 0.101, i8 0.201, i9 0.121, i10 0.139, i12 0.139
The objective (risk) deteriorated from 0.672 to 0.760. The sum of the \(K=3\) largest allocations is binding at 0.201+0.160+0.139=0.5.
Obviously, the number of constraints can become rather large. E.g. for \(N=100, K=3\), we have\[{N \choose K} = {100 \choose 3} = 161,700\]
An integer formulation
A bound for the largest \(\color{darkred}x_i\), can be easily found using: \[\color{darkred}y_1 \ge \color{darkred}x_i\>\forall i\] For the second largest, we allow one exception in the inequalities: \[\begin{align}& \color{darkred}y_2 \ge \color{darkred}x_i - M\cdot\color{darkred}\delta_{2,i} & \forall i \\ & \sum_i \color{darkred}\delta_{2,i} = 1 \\ & \color{darkred}\delta_{2,i} \in \{0,1\} \end{align}\] The exception can take place at the largest, so we have a bound on the second largest. This may require a bit of thinking before this is clear. The third largest is similar:\[\begin{align}& \color{darkred}y_3 \ge \color{darkred}x_i - M\cdot\color{darkred}\delta_{3,i} & \forall i \\ & \sum_i \color{darkred}\delta_{3,i} =2 \\ & \color{darkred}\delta_{3,i} \in \{0,1\} \end{align}\] Here, we allow two exceptions on \(\color{darkred}y_3\ge \color{darkred}x_i\). As our values \(\color{darkred}x_i\) are between 0 and 1, we can simply use \(\color{darkblue}M=1\).
When we run this model, we see:
---- 154 ***MIQP model***
VARIABLE z.L = 0.760 objective
---- 154 VARIABLE x.L fractions
i1 0.160, i2 0.139, i3 0.001, i7 0.101, i8 0.201, i9 0.121, i10 0.139, i12 0.139
Compact formulation
A compact linear formulation is from [1]: \[\begin{align}&\color{darkblue}K\cdot\color{darkred}t +\sum_i \color{darkred}u_i \le 0.5 \\ & \color{darkred}t +\color{darkred}u_i \ge \color{darkred}x_i \\ & \color{darkred}u_i \ge 0 \\ & \color{darkred}t\> \text{free}\end{align}\] This also delivers the same solution.
---- 177 ***compact formulation***
VARIABLE z.L = 0.760 objective
---- 177 VARIABLE x.L fractions
i1 0.160, i2 0.139, i3 0.001, i7 0.101, i8 0.201, i9 0.121, i10 0.139, i12 0.139
This is not very intuitive (the formulation is based on a dual model). But it is the most efficient (linear, and just a few extra variables and constraints). Some modeling tools have this built-in as a function (yalmip has
sumk, and cvxpy has
sum_largest).
Notes
- Here we used "sum of k largest" in a \(\le\) constraint. It is also applicable when it appears in an objective, as long as we are minimizing (when maximizing, things become non-convex).
- If you want to minimize (or constrain) just the k-th largest, use the MIP formulation. I have seen such a construct used in some financial models because of regulatory frameworks.
References
- Stephen Boyd, Lieven Vandenberghe, Convex Optimization, 2004, Exercise 5.19.
Appendix: GAMS code
|
$onText
Implement
the constraint
sum of K largest x(i) <= 50% of sum(x)
We do
this as part of a simple random portfolio problem.
$offText
*------------------------------------------------------------------------
* Size of the problem
*------------------------------------------------------------------------
Set i 'number of x(i)' /i1*i12/;
alias(i,j,k);
$set numk 3
* macro is used to dimension sets further down
scalar numk 'number of largest values to consider' /%numk%/;
*------------------------------------------------------------------------
* Random data
* using example from
https://www.cvxpy.org/examples/index.html
*------------------------------------------------------------------------
Parameters
mu(i) 'expected
returns'
sigma(i,j) 'variance/covariance matrix'
;
mu(i) = normal(0,1);
sigma(i,j) = normal(0,1);
sigma(i,j) = sum(k, sigma(k,i)*sigma(k,j));
display mu,sigma;
*------------------------------------------------------------------------
* standard portfolio problem
*------------------------------------------------------------------------
positive variable x(i) 'fractions';
variable z 'objective';
equations
obj 'quadratic
objective'
revenue 'minimum expected revenue'
budget 'sum to 1'
;
obj.. z =e= sum((i,j), x(i)*x(j)*sigma(i,j));
budget.. sum(i,x(i)) =e= 1;
revenue.. sum(i, mu(i)*x(i)) =g= 0.1;
model m1 /all/;
solve m1
minimizing z using qcp;
display "***baseline***",z.l,x.l;
*------------------------------------------------------------------------
* formulation 1: all 3-combinations
* this generates 120 constraints
*------------------------------------------------------------------------
set combinations3(i,j,k) 'all possible 3-combinations';
combinations3(i,j,k) = ord(j)>ord(i) and ord(k)>ord(j);
option
combinations3:0:0:6;
display
combinations3;
equation ksum(i,j,k) 'sum over possible combinations of three terms';
ksum(combinations3(i,j,k))..
x(i)+x(j)+x(k) =l= 0.5;
model m2 /all/;
solve m2
minimizing z using qcp;
display "***all combinations***",z.l,x.l;
*------------------------------------------------------------------------
* formulation 1a: all 3-combinations
* same thing but a bit more general
*------------------------------------------------------------------------
Sets
n 'key
for combinations'
comb(n<,i) 'combinations'
;
$onEmbeddedCode Python:
import itertools
# list of labels from set i
seti = list(gams.get("i"))
# integer K
k = int(list(gams.get("numk"))[0])
# list of combinations
c = list(itertools.combinations(seti,k))
comb = []
for i in range(len(c)):
nn = f'n{i+1}'
thiscomb = c[i]
for j in range(k):
comb.append((nn,thiscomb[j]))
gams.set('comb',comb)
$offEmbeddedCode comb
display n,comb;
equation ksum2(n);
ksum2(n).. sum(comb(n,i), x(i)) =l= 0.5;
model m2a /m1+ksum2/;
solve m2a
minimizing z using qcp;
display "***all combinations v2***",z.l,x.l;
*------------------------------------------------------------------------
* formulation 2
* binary variables are used to say "y[k] >= x[i]
with k-1 exceptions"
*------------------------------------------------------------------------
set kbest 'best k-th' /k1*k%numk%/ ;
variable y(kbest) 'bound on value of k-th best';
binary variable delta(kbest,i) 'exceptions';
scalar M 'big-M' /1/;
Equation
bound(kbest,i) 'bound on
y(kbest)'
exceptions(kbest) 'number of exceptions in yk >= x(i)'
ysum 'limit sum k largest'
;
bound(kbest,i).. y(kbest) =g= x(i) - M*delta(kbest,i);
exceptions(kbest).. sum(i, delta(kbest,i)) =e= ord(kbest)-1;
ysum.. sum(kbest,y(kbest)) =l= 0.5;
model m3 /m1,bound,exceptions,ysum/;
option miqcp=cplex;
solve m3
minimizing z using miqcp;
display "***MIQP model***",z.l,x.l,y.l,delta.l;
*------------------------------------------------------------------------
* formulation 3
* compact formulation
*------------------------------------------------------------------------
variable t;
positive variable u(i);
equations
compact1
compact2(i)
;
compact1.. numK*t + sum(i,u(i)) =l= 0.5;
compact2(i).. t+u(i) =g= x(i);
model m4 /m1,compact1,compact2/;
solve m4
minimizing z using qcp;
display "***compact formulation***",z.l,x.l,t.l,u.l;
|
A somewhat intuitive interpretation of the compact formulation is that t represents the Kth largest value, which gets counted K times, and u_i = max(x_i - t, 0) represents the amount by which x_i exceeds t, to account for the rest of the sum of the K largest values.
ReplyDelete