Equal segmentation of a series of numbers

Question

How can a series of real numbers be split into subsets in a way that the sum of the real numbers within each subset is as equal as possible?

Coming across from StackOverflow this is the first time, I'm posting a question in this forum. Basically I'm chemist (PhD) and IT manager, so somewhat familiar using mathematics, but not in depth. So I apologize for any lack of precision in my question. So, what is it about?

I have a series of n real numbers: $$r_1, r_2, r_3, ... r_n$$ This series should be split into m segments. So I need m-1 segment limits: $$k_1, k_2, k_3, ..., k_{m-1}$$ These build the segments of the series of real numbers: $$[r_1, ..., r_{k_1}], [r_{k_1+1}, ..., r_{k_2}], ..., [r_{k_{m-2}+1}, ..., r_{k_{m-1}}], [r_{k_{m-1}+1}, ..., r_n]$$ To ensure consistency the following conditions have to be met: $$0 < k_i < k_{i+1} < n \quad \text{for} \; i=1, ..., m-2$$ Each segment has a sum of r-values: $$S_j = \sum_{i=k_{j-1}+1}^{k_j} r_i\quad \text{for} \; j=1, ..., m \,\text{ and }k_0 = 0, \,k_m = n$$ So the basic optimization target is to minimize $$target = \sum_{j=1}^mS_j^2 - \frac1m\cdot(\sum_{j=1}^mS_j)^2$$ for all possible segmentations.

Currently I do this more or less brute-force, i.e. scanning all possible segmentations, calculation of the target and choosing the minimum at the end. This works, but is not satisfying from a performance point of view. As I plan to build this into a website with JavaScript longer runtimes have to be avoided.

So is there another approach for solving this optimization problem?

Thanks in advance for any idea or advice.

Answer 1

It turns out that the algorithm I was trying works for the minimization of the sum $\sum_{j=1}^m F(r_{k_{j-1}+1}+\dots+r_{k_j})$ (we put $k_0=0$) for an arbitrary strictly convex function $F$, so I'll describe it and you'll be able to choose your $F$ later. I personally like $F(t)=1/t$ (minimizing the total height), but for reasonable $r_i$, your $F(t)=t^2$ works pretty well too and the corresponding partitions are fairly close.

The idea is the same dynamic programming. We will compute the best partition of all initial segments $[1,k]$ ($k=0,\dots,n$ into $j=1,\dots,m$ pieces to minimize the sum of the corresponding $F$'s. Let $N(k,j)$ be the corresponding minimum. Put $x_i=r_1+\dots+r_i$ ($i=0,\dots,n$).

We trivially have $N(k,1)=F(r_1+\dots+r_k)$. Now, for $j\ge 2$, we have the recursion $$ N(k,j)=\min_{0\le i\le k}[N(i,j-1)+F(x[k]-x[i])]\,. $$ However, if we run it just like that, it will have the running time $n^2m$, which is a bit too much, though even that is sort of tolerable in your range. To accelerate, let us notice the following thing.

If for some fixed $j$ and $k$, the minimum in the recursion is achieved at $i=I$, then for the same $j$ and $k'>k$, the minimum is achieved at $i\ge I$.

Indeed, for $i<I$, we have $$ N(I,j-1)-N(i,j-1)\le F(x_k-x_i)-F(x_k-x_I)\le F(x_{k'}-x_i)-F(x_{k'}-x_I) $$ (here the first inequality is just a restatement of the condition that the minimum for $k,j$ is achieved at $i=I$, and the second inequality follows from the convexity of $F$).

Thus, we can run the recursion not in the natural order $k=1,2,3,\dots$, but in the order of the highest power of $2$ in $k$. So, once we know the best $i$ for the adjacent multiples of $2^q$, we can search for the best $i$ for the multiple of $2^{q-1}$ in the middle only between these two values. That ensures that we check at most $2n$ values to find $N(k,j)$ for all odd multiples of the same power of $2$ regardless of what that power of $2$ is, thus giving the running time $n\log n$ for each $j$.

The rest is absolutely standard for the dynamic programming approach, so I'll just post the code. I programmed in Asymptote (which allowed me to see the resulting geometric layout easily). If translation to JavaScript presents any difficulties for you, just let me know and I'll do it myself. If the time is still unsatisfactory, let me know and I'll think more.

int n=300, m=30;
srand(seconds());


real[] r,x;
r[0]=0;
for(int k=1;k<=n;++k) r[k]=0.3+2*unitrand();
for(int k=rand()%7+1;k<n;k+=rand()%7+1) r[k]=10*unitrand();

x[0]=0; for(int k=1;k<=n;++k) x[k]=x[k-1]+r[k];

real F(real t)
{
//return t^2;
return 1/(t+0.0001*x[n]/n);
}


real[] S,T;
int[][] v;
real t=0;
for(int k=0;k<=n;++k) {t+=r[k]; S[k]=F(t);}
int count=0;

for(int j=2; j<=m;++j)
{
int[] q; q[0]=0; 
int K=floor(log(n)/log(2)+1); 
for(int qq=2^K; qq>0; qq=quotient(qq,2))
for(int k=qq;k<=n;k+=2*qq)
{
int imin=q[k-qq];
int imax=k; if(k+qq<=n) imax=min(k,q[k+qq]);
real M=-1;
for(int i=imin;i<=imax ;++i) 
{
++count;
real MM=S[i]+F(x[k]-x[i]);
if(MM<M || M<0) {M=MM; q[k]=i;}
}
T[k]=M;
}
T[0]=F(0);
S=copy(T);
v[j]=copy(q);
}
int[] kk; kk[m]=n; kk[0]=0; for(int j=m; j>=2; --j) kk[j-1]=v[j][kk[j]];
write(kk);
write(count, S[n]); pause();

size(400);
real H=0;
for(int mm=1;mm<=m;++mm)
{
real h=1/(x[kk[mm]]-x[kk[mm-1]]); H-=h;
for(int j=kk[mm-1]; j<kk[mm]; ++j) 
draw(box((h*(x[j]-x[kk[mm-1]]),H),(h*(x[j+1]-x[kk[mm-1]]),H+h)));  
}

Answer 2

Excellent Algorithm

Hi @fedja, after streamlining the JavaScript program a little bit and some performance testing I'm really happy with the solution provided by you. Whereas the performance of the brute force recursion is catastrophic, the runtimes of your algorithm are nearly independent of n and m. The algorithm is outstanding!

Here are the results of tests with n = 20 ... 200 and m = 3 ... 10

If I'll use, share or distribute my application, I would like to add an acknowledgement for the algorithm in the docs. How should I show a reference to you?

At the end there are still 2 minor questions concerning the algorithm:

1. What is the reason for the line for(int k=rand()%7+1;k<n;k+=rand()%7+1) r[k]=10*unitrand();

2. Why do you use 1/(t+0.0001*x[n]/n) instead of 1/t for the F-function? as t>0 there is no need for avoiding division by zero

This brings me to one point I omitted in the definition of the problem: I did not state that $r_i$ > 0 perhaps it's a good idea to edit this, what do you think?

Conclusion: Everything works fine THANK YOU

This was the first time I posted a problem here and I'm really happy on the outcome