# Efficiently summing sets of piecewise linear functions

I am doing a type of permutation test on a collection of continuous piecewise linear functions, which I'm representing in Mathematica using Piecewise. Specifically, this involves randomly selecting subsets the functions of a given size $n$, computing the pointwise averages over the subset and its complement, then integrating the difference between these. So, mathematically, we're looking at computing $$I(S)=\int_{-\infty}^\infty\left(\frac{1}{|S|}\sum_{f\in S}f(x)-\frac{1}{|X\setminus S|}\sum_{f\in X\setminus S}f(x)\right)dx$$ which naturally translates pretty naturally into Mathematica code.

You'd think this wouldn't be too hard, since the pieces are linear and all, but it is turning out to be prohibitively slow. I think that the problem is Mathematica is adding the functions inefficiently during the summation steps. The intervals are not the same on almost any of the functions, so the resulting summed function will end up having a very fine interval mesh.

When I look at the outcome of the summation, it is a sum of piecewise functions, rather than a single piecewise function. The summation itself happens very fast (I assume because, under the hood, this is done symbolically with pointers) but when I go to integrate, it takes 20 minutes.

I have tried using PiecewiseExpand before integrating, but this is even slower, even if I do it step by step like this:

summedFunction = S[[1]];
For[i=2, i <= Length[S], i++,
summedFunction = PiecewiseExpand[summedFunction + S[[i]]];
];


I'm starting to worry that this is simply a hard ceiling on Mathematica's computation ability, and that I will have to switch to a faster language. But maybe somebody else knows of a better way to speed this up.

Since I'm sure it will be asked, here is the exact code I'm using.

n1 = Length[Data1];
n2 = Length[Data2];
n = n1 + n2;
Results = {};
Do[
With[{X=Permute[Data1~Join~Data2,RandomPermutation[SymmetricGroup[n]]]},
m1 = (1/n1)Total[X[[1;;n1]]];
m2 = (1/n2)Total[X[[n1+1;;-1]]];
m = m1-m2;
value = Integrate[m, {x, -\[Infinity], \[Infinity]}];
AppendTo[Results, value];
];
,
{PermutationSampleSize}
];
Results = Sort@Results;
With[{X = Data1~Join~Data2},
m1 = (1/n1)Total[X[[1;;n1]]];
m2 = (1/n2)Total[X[[n1+1;;-1]]];
m = m1-m2;
value = Integrate[m, {x, -\[Infinity], \[Infinity]}];
ActualResult = value;
];
pValue = N[Count[Results, _?(# > ActualResult &)]/PermutationSampleSize];
{pValue, ActualResult, Results}


The piecewise functions themselves are built in an irrelevant and complicated way, but, they are just continuous piecewise linear functions, $a_ix+b_i$ with $a_i$ and $b_i$ at numerical precision.

• If you suspect that computing the sum of the piecewise linear functions is the bottleneck, have you tried taking the sum of the integrals instead of the integral of the sum? ...Actually, it seems to me you could just replace each function with its integral and do your permutation test directly on the resulting scalars. – Rahul Oct 21 '15 at 23:44
• The problem is it's not the same with respect to translation when we take the difference. – Alexander Gruber Oct 21 '15 at 23:47
• I do actually need to compute the average function, as I need to use it for other stuff. – Alexander Gruber Oct 21 '15 at 23:48
• Are these functions $C^0$, or step functions, or a mixture of the two? – J. M. is away Oct 22 '15 at 0:04
• They are all $C^0$. – Alexander Gruber Oct 22 '15 at 5:13