## Thursday, June 2, 2022

### A subset selection problem

Selecting $$k$$ out of $$n$$ can be a computational difficult task. Here we look at a special case and see how we can exploit the properties of the data. Suppose we have $$n$$ data points. We want to select $$k$$ points such that the average of these $$k$$ points is as close to a target value $$\color{darkblue}t$$ as possible [1]. I.e. we want $\min_{\color{darkred}S} \left| \frac{\sum_{s \in \color{darkred}S}\color{darkblue}p_s}{\color{darkblue}k} - \color{darkblue}t \right|$ where $$|\color{darkred}S|=k$$. The special property of the data is that it is composed of just a few different values. Let's say our data looks like:

----     21 PARAMETER p  data

i1   1,    i2   3,    i3   2,    i4   1,    i5   1,    i6   1,    i7   2,    i8   3
i9   1,    i10  2,    i11  3,    i12  2,    i13  3,    i14  3,    i15  1,    i16  2
i17  1,    i18  1,    i19  3,    i20  2,    i21  2,    i22  2,    i23  1,    i24  1
i25  2,    i26  3,    i27  1,    i28  2,    i29  3,    i30  1,    i31  1,    i32  2
i33  1,    i34  3,    i35  1,    i36  1,    i37  2,    i38  3,    i39  2,    i40  2
i41  2,    i42  1,    i43  1,    i44  1,    i45  2,    i46  1,    i47  2,    i48  2
i49  3,    i50  1,    i51  2,    i52  3,    i53  2,    i54  1,    i55  1,    i56  1
i57  2,    i58  2,    i59  1,    i60  3,    i61  1,    i62  1,    i63  2,    i64  3
i65  1,    i66  1,    i67  2,    i68  2,    i69  2,    i70  2,    i71  1,    i72  1
i73  1,    i74  3,    i75  2,    i76  3,    i77  1,    i78  1,    i79  3,    i80  1
i81  1,    i82  1,    i83  1,    i84  2,    i85  1,    i86  1,    i87  1,    i88  1
i89  1,    i90  3,    i91  3,    i92  2,    i93  2,    i94  3,    i95  2,    i96  3
i97  1,    i98  3,    i99  1,    i100 2

----     21 PARAMETER k                    =           50  number to select
PARAMETER t                    =      2.12345  target for mean of selected values

In this post, I will start with a generic model using $$n=100$$ binary variables. A second model will exploit that we have just a few unique values in our data. We will exploit that in a model with just 3 integer variables. These models are selecting exactly $$k$$ points.

When we consider the problem: allow up to $$kmax=50$$ points to be selected, things become more difficult. The obvious approach would be to use a loop, try each $$k$$ and pick the best solution. A more complex, non-linear model can be formed. With some effort, we can linearize this model.

### Model 1: binary variables

First, let's look at a generic model. Here we don't exploit the properties of the data. We introduce binary variables: $\color{darkred}\delta_i = \begin{cases} 1 & \text{if point i is selected} \\ 0 & \text{otherwise}\end{cases}$ The absolute value is easily modeled through variable splitting:

Generic Model
\begin{align} \min\> & \color{darkred}d^+ + \color{darkred}d^- \\ & \color{darkred}d^+ - \color{darkred}d^- = \sum_i \frac{\color{darkred}\delta_i \cdot \color{darkblue} p_i}{\color{darkblue} k} - \color{darkblue}t \\ & \sum_i \color{darkred}\delta_i = \color{darkblue}k \\ & \color{darkred}\delta_i \in \{0,1\} \\ & \color{darkred}d^+,\color{darkred}d^- \ge 0 \end{align}

The output is:

----     62 VARIABLE delta.L  selection of points

i2  1,    i5  1,    i6  1,    i7  1,    i8  1,    i9  1
i11 1,    i12 1,    i13 1,    i14 1,    i15 1,    i16 1
i17 1,    i18 1,    i19 1,    i22 1,    i26 1,    i27 1
i28 1,    i29 1,    i30 1,    i31 1,    i33 1,    i34 1
i35 1,    i36 1,    i38 1,    i41 1,    i49 1,    i52 1
i55 1,    i56 1,    i58 1,    i60 1,    i61 1,    i62 1
i64 1,    i67 1,    i68 1,    i74 1,    i76 1,    i79 1
i84 1,    i90 1,    i91 1,    i92 1,    i93 1,    i94 1
i96 1,    i98 1

----     62 PARAMETER value  selected values

i2  3,    i5  1,    i6  1,    i7  2,    i8  3,    i9  1
i11 3,    i12 2,    i13 3,    i14 3,    i15 1,    i16 2
i17 1,    i18 1,    i19 3,    i22 2,    i26 3,    i27 1
i28 2,    i29 3,    i30 1,    i31 1,    i33 1,    i34 3
i35 1,    i36 1,    i38 3,    i41 2,    i49 3,    i52 3
i55 1,    i56 1,    i58 2,    i60 3,    i61 1,    i62 1
i64 3,    i67 2,    i68 2,    i74 3,    i76 3,    i79 3
i84 2,    i90 3,    i91 3,    i92 2,    i93 2,    i94 3
i96 3,    i98 3

----     62 PARAMETER vcount  count of each value

1 16,    2 12,    3 22

----     62 VARIABLE dev.L  deviation

- 0.003450

----     62 VARIABLE z.L                   =     0.003450  objective


This model has $$n=100$$ binary variables.

Note that we can scale the objective by a constant value, and write: $\color{darkred}d^+ - \color{darkred}d^- = \sum_i \color{darkred}\delta_i \cdot \color{darkblue} p_i - \color{darkblue}t\cdot\color{darkblue} k$ As we see later on, this is a valid approach if $$\color{darkblue}k$$ is a variable.

### Model 2: Exploit special properties of the data

Instead of using binary variables, we can keep track of the number of occurrences of each unique value $$1,2,3$$. So we define the variable: $\color{darkred}n_j \in \{0,\dots,\color{darkblue}N_j\}$ indicating how many of unique value $$j$$ we select. This reduces the size of the model to just three integer variables.

The model looks like:

Specialized Model
\begin{align} \min\> & \color{darkred}d^+ + \color{darkred}d^- \\ & \color{darkred}d^+ - \color{darkred}d^- = \sum_j \frac{\color{darkred}n_j \cdot \color{darkblue} j}{\color{darkblue} k} - \color{darkblue}t \\ & \sum_j \color{darkred}n_j = \color{darkblue}k \\ & \color{darkred}n_j \in \{0,1,\dots,\color{darkblue}N_j\} \\ & \color{darkred}d^+,\color{darkred}d^- \ge 0 \end{align}

Here $$\color{darkblue}N_j$$ is the number of data points with $$\color{darkblue}p_i = j$$. For our data set these look like:

----     70 PARAMETER tab  tabulation

1 45,    2 33,    3 22


A solution is:

----     88 VARIABLE num.L  number selected from bin j

1  6,    2 32,    3 12

----     88 VARIABLE dev.L  deviation

- 0.003450

----     88 VARIABLE z.L                   =     0.003450  objective


We note that this is a different solution, but with the same optimal objective function value.

### Extension: select up to $$kmax$$ points

A simple extension of this model can look like: instead of exactly selecting $$k$$ points, select up to $$kmax$$ points. This potentially should reduce the objective. The simplest would be to loop over $$k=1,\dots,kmax$$, solve each model, and pick the best solution. This means solving 50 MIP models. When I do this, the results are:

----    119 **** Loop

----    119 PARAMETER traceObj  trace objective values

i1  0.123450,    i2  0.123450,    i3  0.123450,    i4  0.123450
i5  0.076550,    i6  0.043217,    i7  0.019407,    i8  0.001550
i9  0.012339,    i10 0.023450,    i11 0.032541,    i12 0.040117
i13 0.030396,    i14 0.019407,    i15 0.009883,    i16 0.001550
i17 0.005803,    i18 0.012339,    i19 0.018187,    i20 0.023450
i21 0.019407,    i22 0.012914,    i23 0.006985,    i24 0.001550
i25 0.003450,    i26 0.008065,    i27 0.012339,    i28 0.016307
i29 0.014481,    i30 0.009883,    i31 0.005582,    i32 0.001550
i33 0.002238,    i34 0.005803,    i35 0.009164,    i36 0.012339
i37 0.011685,    i38 0.008129,    i39 0.004755,    i40 0.001550
i41 0.001499,    i42 0.004402,    i43 0.007171,    i44 0.009814
i45 0.009883,    i46 0.006985,    i47 0.004210,    i48 0.001550
i49 0.001001,    i50 0.003450

----    119 PARAMETER bestObj              =     0.001001  best objective

----    119 PARAMETER bestNum  best selection

1  5,    2 33,    3 11


We see the objective goes up and down for different values of $$k$$. The best is using $$k=49$$ points with an objective value of 0.001.

The obvious question is: Can we do this in one model instead?

The model for this could look like:

Incorrect Extension Model
\begin{align} \min\> & \color{darkred}d^+ + \color{darkred}d^- \\ & \color{darkred}d^+ - \color{darkred}d^- = \sum_j \color{darkred}n_j \cdot \color{darkblue} j - \color{darkred}k \cdot \color{darkblue}t \\ &\color{darkred}k = \sum_j \color{darkred}n_j \\ & \color{darkred}n_j \in \{0,1,\dots,\color{darkblue}N_j\} \\ & \color{darkred}d^+,\color{darkred}d^- \ge 0 \\ & \color{darkred}k \in \{1,\dots,\color{darkblue}{\mathit{kmax}}\} \end{align}

Here $$\color{darkred}{k}$$ is now a variable, so we prefer not to divide by $$\color{darkred}k$$. That would create an MINLP model. It looks like, we can multiply the objective by $$\color{darkred}k$$: this makes the problem linear again. But, but.... we have changed the problem a bit: we now prefer smaller $$\color{darkred}k$$ above larger ones. Multiplication by $$\color{darkred}k$$ this way is only allowed for constant $$k$$.

We need to be a bit more careful. In model 1, we measure the deviation as: $\color{darkred}\Delta = \sum_i \frac{\color{darkblue}p_i\cdot \color{darkred}\delta_i}{\color{darkred}k}-\color{darkblue}t$ We can rewrite this as: $\color{darkred}\Delta\cdot \color{darkred}k = \sum_i \color{darkblue}p_i\cdot \color{darkred}\delta_i -\color{darkblue}t\cdot\color{darkred}k$ or $\color{darkred}\Delta\cdot \sum_i \color{darkred}\delta_i = \sum_i \color{darkblue}p_i\cdot \color{darkred}\delta_i -\color{darkblue}t\cdot\color{darkred}k$ The first term is a series of multiplications of a continuous variable with a binary variable. This we know how to linearize. A slightly simpler linearization can be applied to non-negative variables. We can actually use non-negative variables by writing $\sum_i \color{darkred}\delta_i\cdot \color{darkred}d^+ - \sum_i \color{darkred}\delta_i\cdot \color{darkred}d^-$ Linearization the product of a binary variable and a continuous, non-negative variable can be accomplished as follows:  \begin{aligned} & \color{darkred}y=\color{darkred}x\cdot\color{darkred}\delta \\ &\color{darkred}x,\color{darkred}y \in [0,\color{darkblue}U] \\ & \color{darkred}\delta \in \{0,1\} \end{aligned} \Longleftrightarrow \begin{aligned} & \color{darkred}y \le \color{darkred}\delta\cdot\color{darkblue}U \\ & \color{darkred}y \le \color{darkred}x \\ & \color{darkred}y \ge \color{darkred}x - (1-\color{darkred}\delta) \cdot \color{darkblue}U\\ & \color{darkred}x,\color{darkred}y \in [0,\color{darkblue}U] \\ & \color{darkred}\delta \in \{0,1\}\end{aligned} When we assemble this, we have:

Linearized Extension Model
\begin{align}\min\> & \color{darkred}d^+ + \color{darkred}d^- \\ & \sum_i \left(\color{darkred}y_i^+ - \color{darkred}y_i^- \right) = \sum_i \color{darkred}\delta_i\cdot \color{darkblue}p_i - \color{darkred}k \cdot \color{darkblue}t \\ & \color{darkred}y_i^{pm} \le \color{darkred}\delta_i\cdot\color{darkblue}U && \forall i, {pm}\in\{+,-\} \\ & \color{darkred}y_i^{pm} \le \color{darkred}d^{pm} && \forall i, {pm}\in\{+,-\} \\ &\color{darkred}y_i^{pm} \le \color{darkred}d^{pm} -(1-\color{darkred}\delta_i)\cdot \color{darkblue}U &&\forall i, {pm}\in\{+,-\} \\ & \color{darkred} k \in \{1,\dots,\color{darkblue}{\mathit{kmax}}\} \\ & \color{darkred}y_i^{pm} \in [0,\color{darkblue}U] \\ & \color{darkred}\delta_i \in \{0,1\} \end{align}

Before we can run this model we need to establish a value for $$\color{darkblue}U$$. I used the safe value: $\color{darkblue}U := \color{darkblue}{\mathit{kmax}}\cdot \left(\max_i \color{darkblue}p_i-\min_i \color{darkblue}p_i\right)$

Note that in this model we explicitly need to forbid $$\color{darkred}k=0$$.

The solution is:

----    162 VARIABLE kvar.L                =           49  k as variable

----    162 PARAMETER vcount  count of each value

1  5,    2 33,    3 11

----    162 VARIABLE z.L                   =     0.001001  objective


This is the same solution as we found using the loop approach.

### Conclusions

The first two models deal with the "select exactly $$k$$ points" case. The specialized model 2 exploits the property of the data that there are just a few unique values. This can make the model easier to solve than the generic model 1. Especially when using less sophisticated solvers. High-end solvers do like binary variables: I found no real gain for those solvers. For less sophisticated solvers, however, this trick can help a lot.

The problem "allow up to $$kmax$$ points" turned out to be way more challenging to formulate in a linear fashion than I anticipated. But with some effort, it can be done. This is an example where a small change in the model specification has a lot of consequences.

### Appendix: GAMS Model

 *------------------------------------------------------------ * data *------------------------------------------------------------ set    i 'data points' /i1*i100/    j 'unique values' /1,2,3/    pm /'+','-'/  ; * * data drawn from set j * parameter p(i) 'data'; p(i) = uniformint(1,card(j)); p(i) = sum(j$(ord(j)=p(i)),j.val); option p:0:0:8; scalars k 'number to select' /50/ t 'target for mean of selected values' /2.12345/ ; option k:0,t:5; display p,k,t; *------------------------------------------------------------ * model 1 * binary variables delta(i) *------------------------------------------------------------ variables z 'objective' delta(i) 'selection of points' dev(pm) 'deviation' ; binary variable delta; positive variable dev; Equations obj 'objective' deviation 'deviation of mean from target' count 'number of selected points' ; obj.. z =e= sum(pm,dev(pm)); deviation.. dev('+')-dev('-') =e= sum(i,delta(i)*p(i))/k-t; count.. sum(i,delta(i)) =e= k; model m1 /obj,deviation,count/; solve m1 minimizing z using mip; parameter value(i) 'selected values'; value(i) = p(i)*delta.l(i); parameter vcount(j) 'count of each value'; vcount(j) = sum(i$(value(i)=j.val), 1); option delta:0:0:6, value:0:0:6, vcount:0; option dev:6,z:6; display delta.l,value,vcount,dev.l,z.l; *------------------------------------------------------------ * model 2 * integer variable num(j) *------------------------------------------------------------ parameter tab(j) 'tabulation'; tab(j) = sum(i$(p(i)=j.val),1); option tab:0; display tab; integer variable num(j) 'number selected from bin j'; equation idev, icount; idev.. dev('+')-dev('-') =e= sum(j,num(j)*j.val)/k-t; icount.. sum(j,num(j)) =e= k; num.up(j) = tab(j); model m2 /obj,idev,icount/; solve m2 minimizing z using mip; option num:0; display num.l,dev.l,z.l; *------------------------------------------------------------ * model 3 * best k=1...Kmax using a loop *------------------------------------------------------------ scalar kmax; kmax = k; * fast solve loop option solprint=silent, solvelink=5; parameter bestObj 'best objective' /INF/ bestNum 'best selection' traceObj 'trace objective values' ; for(k=1 to kmax, solve m2 minimizing z using mip; traceObj(i)$(ord(i)=k) = z.l;    if (z.l