# How to generalize and speed up this program?

There is a vector whose length is $n*m$, we can part it into $n$ sections in order. I want to define a function to compare the $i$-th and $i+1$-th section and find the minimum. The program for $n=m=3$ is as following:

f[x_] := Min[
Total[Total[Boole[Thread[x[[1 ;; 3]] < ##]]] & /@ x[[4 ;; 6]]],
Total[Total[Boole[Thread[x[[4 ;; 6]] < ##]]] & /@ x[[7 ;; 9]]],
Total[Total[Boole[Thread[x[[7 ;; 9]] < ##]]] & /@ x[[1 ;; 3]]]];


Could we generalize it to

f[x_,n_,m_] := Min[
Total[Total[Boole[Thread[x[[1 ;; m]] < ##]]] & /@ x[[m+1 ;; 2*m]]],
Total[Total[Boole[Thread[x[[m+1 ;; 2*m]] < ##]]] & /@ x[[2*m+1 ;; 3*m]]],
....
Total[Total[Boole[Thread[x[[(n-1)*m+1 ;; n*m]] < ##]]] & /@ x[[1 ;; m]]]];


By the way, could we speed up this code? Any help or suggestion will be appreciated!

-

## A new solution

I realized that comparing each and every value in the sections might be inefficient, especially in cases where the sections are long. Instead we need only the relative ordering of these elements from which we can compute the number of Less pairs.

Here is my solution:

Edit: Ray Koopman provided a greatly improved counting method (applied to the binarized orderings) which I have now incorporated. To see my original code look at the edit history of this post and the compiled count function.

f5[v_, n_] :=
Min[UnitStep[(n + 0.5) - Ordering /@ Partition[Reverse@v, 2 n, n, -n]].Range[2 n]] -
n (n + 1)/2


With long segments this is much faster than my first method (f2) and also a lot faster than Michael's code.

If you desire LessEqual behavior you can use:

f5LE[v_, n_] :=
Min[UnitStep[Ordering /@ Partition[v, 2 n, n, -n] - (n + 0.5)].Range[2 n]] - n(n+1)/2


### Explanation

This method works on the principle that a direct comparison as performed in the question or the other answers is an $n^2$ operation, whereas a sort is on average an $n\log(n)$ operation.

I start by splitting the list into lengths of $2 n$ with an offset of $n$ and an alignment of $-n$, which looks like this:

v = CharacterRange["a", "h"];
n = 2;

Partition[v, 2n , n, -n]


{{g,h,a,b}, {a,b,c,d}, {c,d,e,f}, {e,f,g,h}}

Each partition represents two adjacent segments, i.e. 4;;6 and 7;;9 in the question example.

I actually need the elements of each of these partitions in reverse order (see below), and it's more efficient to do that before partitioning, and for my purpose equivalent:

Sort[Reverse /@ Partition[v, 2 n, n, -n]] ===
Sort[Partition[Reverse@v, 2 n, n, -n]]


True

(The Sorts are only to show that the same partitions exist.)

The next step, and the core of the method as explained above, is Ordering, which gives the relative position of every element in a sorted list. For example:

Ordering @ {"c", "a", "f", "e", "b", "d"}


{2, 5, 1, 6, 4, 3}

Note that "a" would appear at position 1 in a sorted list but is at position 2 in the original - therefore the first element of the ordering is 2; likewise "b" would appear at position 2 in a sorted list but is at position 5 in the original, therefore the second element of the ordering is 5.

In this ordering list I then want to round every value up to n down to 0 and every value above that to 1, which I did like this:

UnitStep[{2, 5, 1, 6, 4, 3} - 3.5]


{0, 1, 0, 1, 1, 0}

Every zero represents a value from the first half of the partition, representing the first segment, and every one represents a value from the second half of the partition, representing the second segment. Equivalent but slower is a decorate-and-sort as follows:

Transpose @ {{"c", "a", "f", "e", "b", "d"}, {0, 0, 0, 1, 1, 1}}
SortBy[%, {First}]
Last /@ %


{{"c", 0}, {"a", 0}, {"f", 0}, {"e", 1}, {"b", 1}, {"d", 1}}

{{"a", 0}, {"b", 1}, {"c", 0}, {"d", 1}, {"e", 1}, {"f", 0}}

{0, 1, 0, 1, 1, 0}

We can then count the number of ones that are to the left of every zero to determine how many values from segment 2 are less than a value from segment one. This is what the compiled count function does. Edit: now done with Ray Koopman's dot product method.

The Reverse was needed to get the behavior requested for the case of duplicate elements. Note the difference in the ranking of the "c" elements in this example:

s = {"c", "a", "f", "e", "c", "d"};

Transpose @ {s, {1, 1, 1, 0, 0, 0}};
SortBy[%, {First}]

Transpose @ {Reverse[s], {0, 0, 0, 1, 1, 1}};
SortBy[%, {First}]


{{"a", 1}, {"c", 1}, {"c", 0}, {"d", 0}, {"e", 0}, {"f", 1}}

{{"a", 1}, {"c", 0}, {"c", 1}, {"d", 0}, {"e", 0}, {"f", 1}}

### Room for improvement

I believe that there is still considerable redundancy in this code. Specifically, as you are seeking only the minimum value we would ideally abort an Ordering operation as soon as it was apparent that the existing minimum was reached or exceeded. This would require a manual implementation of the Ordering function and incorporation of the counting. Doing this efficiently it out of my experience and beyond my reach, as version 7 which I use cannot compile to C. Others with the experience and tools may which to implement this idea to see if it has merit.

## Timings

With short segments Michael's f0 is fastest:

SetAttributes[timeAvg, HoldFirst]

timeAvg[func_] :=
Do[If[# > 0.3, Return[#/5^i]] & @@ Timing@Do[func, {5^i}], {i, 0, 15}]

a = RandomInteger[9999, 50000];
f2[a, 5] // timeAvg
f0[a, 5] // timeAvg
f5[a, 5] // timeAvg


0.1934

0.0018224

0.008488

However, with longer segments the picture changes considerably:

a = RandomInteger[9999, 50000];
f2[a, 500] // timeAvg
f0[a, 500] // timeAvg
f5[a, 500] // timeAvg


7.363

0.1154

0.006992

With even longer segments f2 is too slow to compare, and we find that f0 consumes a great deal of memory, while f3 does not:

a = RandomInteger[9999, 150000];
f0[a, 1500] // timeAvg
MaxMemoryUsed[]


0.905

1816949160

a = RandomInteger[9999, 150000];
f5[a, 1500] // timeAvg
MaxMemoryUsed[]


0.02308

19170792

In summary f5 is efficient with all segment lengths on long vectors:

a = RandomInteger[99999, 1*^6];
timeAvg @ f5[a, #] & /@ {2, 10, 100, 1000, 10000, 100000}


{0.374, 0.1434, 0.1248, 0.1498, 0.1872, 0.312}

-
Looks neat. I'll have to study it later when I have time. Fixed my answer. – Michael E2 Aug 24 '13 at 11:52
Dear Mr. Wizard, caould you explain your code in details? It is perfect but too hard for me to read. – Eden Harder Aug 24 '13 at 12:30
@Eden I added an explanation but I don't know if it is clear. Let me know if any part needs more or different explanation. – Mr.Wizard Aug 24 '13 at 18:07
@Mr.Wizard Perfect! How could you come up with this code which is not intuitive? – Eden Harder Aug 25 '13 at 0:40
This is just f3, modified to not need the compiled subroutine: f4[v_, n_] := With[{z = Join[Table[1,{n}],Table[0,{n}]]}, Min[ z[[Ordering@#]]& /@ Partition[Reverse@v,2n,n,-n].Range[2n] ] - n(n+1)/2 ] – Ray Koopman Aug 25 '13 at 3:57

Here is another way, which is pretty fast. The number n of sections does not matter as an argument, so it is omitted.

f0[x_, m_] := With[{v = Partition[x, m]},
m^2 - Max @ Total[UnitStep[Transpose @ ConstantArray[v, m] - RotateLeft[v]], {2, 3}]]


Example 1:

x = {1, 16, 3, 6, 20, 5, 13, 14, 2};
f0[x, 3]
(* 4 *)


Example 2:

SeedRandom[1];
x2 = RandomInteger[10, 30 * 5000];
f0[x2, 30]
(* 244 *)


Explanation

In part, I was primed by this answer.

Basically, each element in your list, such as,

Total[Total[Boole[Thread[x[[1 ;; 3]] < ##]]] & /@ x[[4 ;; 6]]]


is equivalent to counting how many of the following elements in an mxm matrix are negative, where m = 3 here:

With[{x = {a, b, c, d, e, f, g, h, i}},
Thread[x[[1 ;; 3]] - #] & /@ x[[4 ;; 6]]
]
(* {{a - d, b - d, c - d},
{a - e, b - e, c - e},
{a - f, b - f, c - f}} *)


If we reverse the subtractions, the number we seek is equivalent to m^2 minus the total number of nonnegative entries. Unitstep conveniently replaces the matrix entries by 1 if the entry is nonnegative and 0 if negative. It is an efficient way to deal with inequalities on packed arrays. Total[matrices, {2, 3}] totals each matrix in a list of matrices and returns a list of the results.

I'll now turn to the crux of the matter, which is how to make the list of matrices. This bit

Transpose @ ConstantArray[v, m] - RotateLeft[v]


creates the list of such mxm matrices of differences. I'll try to explain how. Arithmetic functions have the attribute Listable. This means that when one or more arguments are lists, Mathematica tries to thread the function, say Plus, through the lists until it matches up numbers or other symbols that be added. The trick is to get the lists in the right shapes so that the numbers match up the way you desire.

In the second argument, RotateLeft lines up adjacent segments in v as well as the first and last segments. (RotateRight might have been adapted, too.)

v = {{a, b, c}, {d, e, f}};
RotateLeft[v]
(* {{d, e, f}, {a, b, c}} *)


To understand the first argument, consider the two operations:

{a, b, c} - {d, e, f}
{a, b, c} - d
(* {a - d, b - e, c - f} *)
(* {a - d, b - d, c - d} *)


In the first line, the elements are matched up pairwise; in the second each of a, b, c is matched with d. Clearly we want something like this to execute by our code:

{{a, b, c} - d,
{a, b, c} - e,
{a, b, c} - f}


"Unthreading" this, we get the following, which can be checked to give the desired output:

{{a, b, c}, {a, b, c}, {a, b, c}} - {d, e, f}


From this we see we need m = 3 copies of first list, from which we should subtract the second list. ConstantArray is used to generate the copies.

One snag is that v is a list of segments and ConstantArray does not line up the copies right.

v = {{a, b, c}, {d, e, f}}
ConstantArray[v, 3]
(* {{{a, b, c}, {d, e, f}},
{{a, b, c}, {d, e, f}},
{{a, b, c}, {d, e, f}}} *)


The copies of {a, b, c} are not together. The list v has length 2 and the ConstantArray has length 3. The elements do not match up (and ConstantArray[v, 3] - RotateLeft[v] will generate error messages). But they are all there, and they need to be rearranged. Transpose is one of the functions that can rearrange lists (if they are regular arrays). Flatten, ArrayFlatten and ArrayReshape are some others. They all have a second argument to deal with levels and dimensions. It turns out a simple Transpose works here.

Transpose@ConstantArray[v, 3]
(* {{{a, b, c}, {a, b, c}, {a, b, c}},
{{d, e, f}, {d, e, f}, {d, e, f}}} *)


Now if we subtract RotateLeft[v], each matrix in the list of copies matches up with a segment of RotateLeft[v]. The first step in threading will give

{ {{a, b, c}, {a, b, c}, {a, b, c}} - {d, e, f}, ... }


and we end up with the desired result as above.

One can also see that these copies would use a lot of memory if the length of the segments m is large, as Mr.Wizard demonstrates in his answer.

-
Thanks! Why did you define f0 but use f? Could you explain your elegant code in details? I can not understand how it works. – Eden Harder Aug 24 '13 at 3:44
Really nice code. I knew there would be a faster way that what I gave but thinking in these terms is a bit harder for me and I was out of time. – Mr.Wizard Aug 24 '13 at 6:10
@Mr.Wizard I understand what the code says now but why this code will be faster than yours? – Eden Harder Aug 24 '13 at 9:51
@Eden Two reasons: (1) optimal use of Mathematica's array optimizations because the data is handled in large chunks that are never unpacked; (2) Numeric operations that are equivalent to Less (subtraction and UnitStep), which also do not require unpacking of the array. However, I am working on a new approach that if successful should be considerably faster than this method when the segment length is large. – Mr.Wizard Aug 24 '13 at 10:12
@EdenHarder f was a typo. It started out as f in my notebook, but I wanted to give a unique name so people could cut and paste from the different question/answers (like Mr.Wizard, I suspect). – Michael E2 Aug 24 '13 at 11:10

Here's a function that I think does what you want:

cmp[{sec1_, sec2_}] :=
Total@Flatten@Table[Boole[x < y], {x, sec1}, {y, sec2}]
findMin[vector_, n_] :=
Min[cmp /@ Partition[Partition[vector, n], 2, 1, {1, 1}]]


For example:

vector = RandomInteger[100, 9]


{65, 33, 90, 16, 48, 46, 90, 53, 37}

findMin[vector, 3] (* Three being the length of each segment *)


2

-

Very similar to the code that I see Anon just posted, but different enough I think to post as well:

f2[x_, n_] :=
Total[Boole@Outer[Less, ##], 2] & @@@ Partition[Partition[x, n], 2, 1, 1] // Min


Parameter n is the length of each segment rather than the number of segments; if that's a problem it's not hard to use Length[x] / segments to get the other number.

Example:

x = {1, 16, 3, 6, 20, 5, 13, 14, 2};

f2[x, 3]


4

-

### Comment -- Visualization of Mr.Wizard's answer

I don't know whether Mr.Wizard thought of it this way in [his answer], but if you think about the OP's problem (long enough), you might recognize it as a combinatorial about placing two sets on points on a line. In my case, it was after thinking about Mr.Wizard's second solution that I was reminded of this sort of combinatorial problem. Viewed this way, the connection with Ordering looks more natural. Here is an image corresponding to the two sets {{3, 1, 6}, {5, 2, 4}}. The exact coordinates do not matter, only their relative order.

Here is one corresponding to {{5, 2, 4}, {3, 1, 6}}:

The OP's function sums up how many times an element in the second set is greater than or equal to an element in the first set (or vice versa, ones from the first set are less than ones in the second). The OP's function on the list {3, 1, 6, 5, 2, 4} with segment length n = 3 is thus 4, the minimum of the two sums.

Note that the order in each segment does not matter. It is really just a set. Now Mr.Wizard's f5 (née f3) computes the ordering of pairs of segments of length 2n. Our observation suggests one might be able to figure only the ordering of each segment (or equivalently sort them). Since sorting/ordering has complexity n Log[n], there is the potential for a slight speed-up.

### A slightly different solution

I may be overlooking something, but I could only achieve a fast implementation of the idea of sorting each segment through Compile. Basically, preventing unpacking of packed arrays could be avoid in Compile, and I haven't thought of a vectorized approach.

The function basically does what is shown above. It counts how many of the second segment y are greater than or equal to each element in the first set.

count0 = Compile[{{x, _Real, 1}, {y, _Real, 1}},
Length[x]^2 + Length[x] - Total @
Module[{i = 1},
(While[i <= Length[y] && # >= y[[i]], i++]; i) & /@ x],
CompilationTarget -> "C"
];
f00 = Compile[{{list, _Real, 1}, {m, _Integer}},
Module[{segments = Sort /@ Partition[list, m, m]},
CompilationOptions -> {"InlineExternalDefinitions" -> True},
CompilationTarget -> "C"];


### Timings

Compiling to C helps make up for not using vectorized operations. I used Mr.Wizard's timeAvg and test data. We can see that Mr.Wizard's f5 is a little faster than my f00 on short segments.

a = RandomInteger[9999, 50000];
f00[a, 5] // timeAvg
f5[a, 5] // timeAvg


0.00683867

0.00626927

But f00 is faster on longer segments.

a = RandomInteger[9999, 50000];
f00[a, 500] // timeAvg
f5[a, 500] // timeAvg


0.00549772

0.00873428

### Code for the visualization

upTri[x_] :=
Translate[Scale[Polygon[{{0, 0}, {-1, -1}, {1, -1}}], 0.3, {0, 0}],
x];
downTri[x_] :=
Translate[Scale[Polygon[{{0, 0}, {1, 1}, {-1, 1}}], 0.3, {0, 0}], x];
Manipulate[
With[{counts =
Function[p,
Count[pts[[2 ;; 2 m ;; 2]], _?(First@p < First@# &)]] /@
pts[[1 ;; 2 m ;; 2]]},
Graphics[
{Line[{{0, 0}, {11, 0}}],
Text[#1, #2 - {0, 0.6}] &, {counts, pts[[1 ;; 2 m ;; 2]]}],
EdgeForm[Directive[Thin, Black]], Blue,
upTri[pts[[1 ;; 2 m ;; 2]]], Red, downTri[pts[[2 ;; 2 m ;; 2]]]},
PlotRange -> {{-0.5, 11.5}, 0.7 {-1, 1}},
PlotLabel ->
Graphics[{Blue, upTri[{0, 0}]}, ImageSize -> 16],
" the number of ",
Graphics[{Red, downTri[{0, 0}]}, ImageSize -> 16],
" to the right, the sum is ", Total@counts}
]
],
{{pts, Table[{i, 0}, {i, 10}]}, {0, 0}, {11, 0}, Locator,
Appearance -> None},
{{m, 3}, 2, 5, 1}
]

-
@Mr.Wizard This started just as a way to help Eden Harder visualize your algorithm, but it gave me an idea how to tweak your solution a little. Somehow I feel it might be improved, but I'm drawing a blank. – Michael E2 Aug 25 '13 at 17:08
Yes, this is exactly how I came to visualize the problem when I was working on my second answer. Thanks for the illustration. I haven't read your code yet and I've got to leave, but if you haven't yet you might try my idea under the "Room for improvement" heading in my answer. (It's a bit vague, but better formed in my head; we should chat later.) – Mr.Wizard Aug 25 '13 at 18:16
@MichaelE2 Why f00[a, 500] // timeAvg is faster than f00[a, 5] // timeAvg ? – Eden Harder Aug 26 '13 at 1:11
@EdenHarder I'm not entirely sure, but perhaps it's not that strange. The data in both cases has 50000 elements. In the case of f00[a, 500], it consists of 100 segments of 500 elements. In f00[a, 5], it consists o 10000 segments of 5 elements. It processes the segments differently than the list of them. For instance f00[a, 5000] is slower than f00[a, 500], so it doesn't keep getting faster as the length of the segments grows. The time depends mostly on the total number of elements and varies with how it's segmented. – Michael E2 Aug 26 '13 at 2:16