I would like to find a algorithm that will do the following:

Given two sets $A, B \subseteq \mathbb{R}$, where $|B| > |A|$, find the largest subset $C \subseteq B$, such that:

$\qquad |\operatorname{mean}(A) - \operatorname{mean}(C) |< \delta$ and

$\qquad |\operatorname{std}(A) - \operatorname{std}(C) | < \epsilon$

I think this problem is NP-hard, but I would like a good approximation that does reasonably well. An even harder version of this is to find the larget subset $C$ that matches not just the moments, but the histogram of $A$ to some quality metric. If you have a solution, or can point to same papers that would great!

I am a neuroscientist, and this is for my research.

  • $\begingroup$ Have you tried expressing your problem as LP? $\endgroup$ – Raphael Feb 12 '14 at 8:40
  • $\begingroup$ What are $\text{mean}(A)$ and $\text{std}(A)$ (I assume the mean value and the standard deviation)? Are $\delta$ and $\epsilon$ fixed constants? What happens when you restrict yourself to $A,B \subseteq \mathbb{N}$? $\endgroup$ – Pål GD Feb 12 '14 at 12:03
  • $\begingroup$ I suspect (but I may be wrong) that you just want to find a subset whose mean and std are the "closest" to $mean(A)$ and $std(A)$. Is this correct? If so, you should reformulate the problem accordingly. Also, I suggest that you start with just the mean. $\endgroup$ – Bitwise Feb 12 '14 at 21:33
  • $\begingroup$ @Pal GD - mean(A) and std(A) are mean and standard deviation. The $\delta$ and $\epsilon$ are not fixed, but I am interested in how the subset size $|C|$ varies with the constraint. Bitwise - I am not interested in the "closest" but rather, "given a tolerance, find the largest subset". Getting the mean i think, is easy; just remove points from the extrema until the means are within $\delta$. $\endgroup$ – ht959 Feb 12 '14 at 22:33

I tried formulating this as an integer linear program, but ran into problems, because the variance is not a linear function: it is quadratic.

Here was my attempt.

Introduce variables $x_b$ (one for each $b \in B$), such that $x_b=1$ if $b \in C$ and $x_b=0$ if $b \notin C$. Thus, we have $|B|$ new zero-or-one variables. We can force them to be zero or one by adding the constraints $0 \le x_b \le 1$ for each $b \in B$.

Next, we'll force the mean of $C$ to be in the accepted range. Introduce a new unknown $\mu$ and add the constraint

$$\mu = (\sum_{b \in B} x_b b)/|B|.$$

This forces $\mu$ to be the mean of $C$. Also, it is expressible in a linear equation, since everything except for $\mu$ and the $x$'s are constants. Finally, add the two linear inequalities

$$\text{mean}(A)-\delta \le \mu \le \text{mean}(A)+\delta.$$

But it's not clear how to force the standard deviation of $C$ to be in the accepted range. Let $\sigma = \text{std}(C)$, i.e., $\sigma$ denotes the standard deviation of $C$. We could try to use the fact that

$$\sigma^2 = \frac{1}{N} \sum_{c \in C} c^2 - \frac{1}{N^2} \left( \sum_{c \in C} c \right)^2,$$

but now when we substitute in with our definitions of $x_b$, we get something ugly:

$$\sigma^2 = \frac{1}{N} \sum_{b \in B} b^2 x_b - \frac{1}{N^2} \left( \sum_{b \in B} b x_b \right)^2,$$

where the right-hand term has something quadratic. Worse still, $N$ also depends upon the $x$'s, which takes this even further away from something linear:

$$S = \sum_{b \in B} x_b.$$

Alternatively, we could try to use the fact that

$$\sigma^2 = \frac{1}{N} \sum_{c \in C} (c - \mu)^2,$$

but when we substitute in, again we get something that's not linear:

$$\sigma^2 = \frac{1}{N} \sum_{b \in B} x_b (b - \mu)^2,$$

where the non-linearity comes from the fact that $N$ is itself a function of the $b$'s, as well as the fact that the right-hand side involves a product of two linear variables ($x_b$ and $\mu$).

We want

$$\text{std}(A) - \epsilon \le \sigma^2 \le \text{std}(A) + \epsilon.$$

Plugging in, we find that we need

$$(\text{std}(A) - \epsilon)N \le \sum_{b \in B} x_b (b - \mu)^2 \le (\text{std}(A) + \epsilon)N.$$

This is not linear, because $x_b (b - \mu)^2$ is not a linear function of $x_b,\mu$.

We could deal with the fact that $N$ is not a constant by enumerating all candidate values of $N$ and solving an ILP for each. However, that still doesn't deal with the problem that $x_b (b - \mu)^2$ is non-linear.

So, I'm stuck. I wasn't able to formulate this cleanly as an ILP.

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.