Finding largest subset that matches moments

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.

• Have you tried expressing your problem as LP? – Raphael Feb 12 '14 at 8:40
• 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}$? – Pål GD Feb 12 '14 at 12:03
• 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. – Bitwise Feb 12 '14 at 21:33
• @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$. – 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.