Finding column with maximum sum of by-element products of column subsets efficiently?

A couple of recent questions and associated answers reminded me of an old problem:

As a toy example, take:

SeedRandom[1211]
matsize = 5;
subsize = 3;
mat = RandomInteger[{-10, 10}, {matsize, matsize}];
mat // MatrixForm


$\left( \begin{array}{ccccc} -2 & 8 & 2 & 7 & -7 \\ 5 & -2 & -4 & 2 & -3 \\ 3 & -9 & -1 & -3 & 5 \\ 4 & -7 & 8 & 9 & -7 \\ 4 & 3 & -3 & -8 & 7 \\ \end{array} \right)$

where rows represent some subject/object under test, columns are different sets of conditions, and the elements are an integer score on [-10,10].

It is desired to find the column index whose elements, taken over all subsets of size subsize, have the maximum value when the elements of each subset are multiplied and those products are totaled.

For the toy example, this is column 2, with a result of 487.

A naive way might be

Tr /@ Apply[Times, Subsets[#, {subsize}] & /@ Transpose[mat], {2}] // Ordering // Last


but that will obviously blow RAM for cases with large subset cardinality.

Much more memory efficient, but with its own issues is:

SymmetricPolynomial[subsize, #] & /@ Transpose[mat] // Ordering // Last


How might you do this, with the proviso that I'd prefer pure Mathematica (compilation to WVM is fine, compilation to C means I'rather just write it directly in C or whatever), and the actual use-case is matsize=100 and subsize=30.

Edit: Since adding bounty, I came up with

(Coefficient[Times @@ (1 + #1 x) // Expand, x^#2] // Position[#, Max@#] &) &[mat, subsize]


Fastest way so far (on my box), reasonable RAM use.

• Can one assume that the entries are nonzero ? – A.G. Oct 5 '16 at 2:46
• @A.G. - no, though I'm curious if you think that leads to an optimization of some sort... – ciao Oct 5 '16 at 5:00

Algorithm

Since the bottleneck is calculating elementary symmetric polynomial, let's search for an efficient algorithm to do it. In answer to "Algorithm(s) for computing an elementary symmetric polynomial" question, Ben Kuhn shows simple recurrence relation between elementary symmetric polynomials: \begin{align} \forall_{n\in\mathbb N} \quad & s^0_n(x_1, ..., x_n) = 1 \\ \forall_{k,n\in\mathbb N;k>n} \quad & s^k_n(x_1, ..., x_n) = 0 \\ \forall_{k,n\in\mathbb Z_+} \quad & s^k_n(x_1, ..., x_n) = s^k_{n-1}(x_1, ..., x_{n-1}) + x_n s^{k-1}_{n-1}(x_1, ..., x_{n-1}) \end{align} which relates $k$-th polynomial of $n$ variables with $k$-th and $(k-1)$-th polynomial of $n-1$ variables. Using this recurrence we can create simple iterative algorithm starting with polynomials of one variable, using them to calculate polynomials of two variables and so on.

For example if we want to calculate $k=3$ polynomial of $n=7$ variables we need to calculate polynomials with following $n$ and $k$:

\begin{array}{l|cccc} n & 1 & 2 & 3 & 4 & 5 & 6 & 7 \\ \hline k=0 & 0 & 0 & 0 & 0 \\ k=1 & 1 & 1 & 1 & 1 & 1 \\ k=2 & & 2 & 2 & 2 & 2 & 2 \\ k=3 & & & 3 & 3 & 3 & 3 & 3 \\ \end{array}

for $k=5$ we need

\begin{array}{l|cccc} n & 1 & 2 & 3 & 4 & 5 & 6 & 7 \\ \hline k=0 & 0 & 0 \\ k=1 & 1 & 1 & 1 \\ k=2 & & 2 & 2 & 2 \\ k=3 & & & 3 & 3 & 3 \\ k=4 & & & & 4 & 4 & 4 \\ k=5 & & & & & 5 & 5 & 5 \\ \end{array}

In general we need a list of at most Min[k, n - k] + 1 elements to store already calculated values. Iteration, with i denoting number of variables from 1 to n, can be split to three stages:

1. for i <= n-k we increase number of required polynomials,
2. for n-k+1<= i <=k we need exactly Min[k, n - k] + 1 polynomials,
3. for Max[k,n-k]+1<=i<=n number of required polynomials decreases.

Implementation

Simple procedural implementation is as follows:

symPol // ClearAll
symPol[0, vars_List] = 1;
symPol[1, vars_List] := Total@vars
symPol[k_Integer?Positive, vars_List] /; k > Length@vars = 0;
symPol[k_Integer?Positive, vars_List] /; k === Length@vars := Times @@ vars
symPol[k_Integer?Positive, vars_List] := Module[{n, l, pols, prev, i},
n = Length@vars;
l = Max[k, n - k];
pols = ConstantArray[0, Min[k, n - k] + 1];
pols[[1]] = 1;
pols[[2]] = vars // First;
Do[
pols[[2 ;;]] += vars[[i]] pols[[;; -2]]
,
{i, 2, n - k}
];
Do[
prev = pols[[2 ;;]];
pols *= vars[[i]];
pols[[;; -2]] += prev
,
{i, n - k + 1, k}
];
Do[
prev = pols[[2 ;; l - i]];
pols[[;; l - 1 - i]] *= vars[[i]];
pols[[;; l - 1 - i]] += prev
,
{i, l + 1, n}
];
pols // First
]


symPol gives the same results as built-in SymmetricPolynomial:

Table[
With[{vars = Table[Unique[], n]},
Expand@symPol[k, vars] == SymmetricPolynomial[k, vars] // Simplify
],
{n, 0, 10}, {k, 0, n}
];
And @@ Flatten@%
(* True *)


but it's much faster and uses much less memory:

SymmetricPolynomial[12, Range[25]] // MaxMemoryUsed // Timing
symPol[12, Range[25]] // MaxMemoryUsed // Timing
(* {6.748, 915253280} *)
(* {0., 7984} *)


For actual size matrix from OP:

SeedRandom[1211]
matsize = 100;
subsize = 30;
mat = RandomInteger[{-10, 10}, {matsize, matsize}];


we get:

(result = Ordering[symPol[subsize, #] & /@ Transpose[mat], -1]) // MaxMemoryUsed // RepeatedTiming
(* {0.093, 180552} *)
result
(* {70} *)


in 0.093 seconds using 180552 bytes.

Compilation

Procedural nature of symPol function makes it easily compilable. If machine precision is enough, then we can compile function taking array of reals:

symPolR = Compile[{{k, _Integer}, {vars, _Real, 1}},
Module[{n, l, len, pols, prev, i, j},
If[k == 0, Return[1.]];
n = Length@vars;
If[k == n, Return[Times @@ vars]];
If[k > n, Return[0.]];
len = Min[k, n - k] + 1;
l = Max[k, n - k];
pols = Table[0., {len}];
pols[[1]] = 1.;
pols[[2]] = vars // First;
Do[
pols[[-j]] += vars[[i]] pols[[-j - 1]],
{i, 2, n - k}, {j, len - 1}
];
Do[
Do[
pols[[j]] = pols[[j + 1]] + vars[[i]] pols[[j]],
{j, len - 1}
];
pols[[len]] *= vars[[i]]
,
{i, n - k + 1, k}
];
Do[
pols[[j]] = pols[[j + 1]] + vars[[i]] pols[[j]],
{i, l + 1, n}, {j, len + l - i}
];
pols // First
],
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed",
CompilationTarget -> "WVM"
]


With "WVM" version we get ten times faster function:

(result = Ordering[symPolR[subsize, Transpose@mat], -1]) // MaxMemoryUsed // RepeatedTiming
(* {0.009, 172424} *)
result
(* {70} *)


With CompilationTarget -> "C" we get another factor of ten:

(result = Ordering[symPolR[subsize, Transpose@mat], -1]) // MaxMemoryUsed // RepeatedTiming
(* {0.0007, 164728} *)
result
(* {70} *)

• Quite nice, +1. Very good RAM utilization, nearly as fast as method in my OP edit. Well done. – ciao Oct 3 '16 at 6:14
• @ciao On my computer your method is a bit slower: (Coefficient[Times @@ (1 + #1 x) // Expand, x^#2] // Position[#, Max@#] &) &[mat, subsize]// MaxMemoryUsed // RepeatedTiming (* {0.12, 1 975 496} *), but if I map function on transposed matrix I get same RAM usage as with my method: Coefficient[Times @@ (1 + # x) // Expand, x^subsize] & /@ Transpose@mat // Position[#, Max@#] & // MaxMemoryUsed // RepeatedTiming (* {0.12, 181 184} *) – jkuczm Oct 3 '16 at 13:15
• Both were nice answers - awarded you the bounty for the nice procedure. Thanks! – ciao Oct 6 '16 at 22:02

The (edited) OP gives the following implementation of the time-consuming part of the calculation. (I will ignore the problem of finding the maximum as trivial).

refimplementation = Coefficient[Times @@ (1 + #1 x) // Expand, x^#2] &;


I offer the following slightly different implementation

serimplementation = SeriesCoefficient[Times @@ (1 + #1 x), {x, 0, #2}] &;


For the following test case

SeedRandom[1211]
matsize = 100;
subsize = 30;
mat = RandomInteger[{-10, 10}, {matsize, matsize}];


this implementation is better than 5 times faster on my PC, and gives the same answer

RepeatedTiming[v1 = refimplementation[mat, subsize];]
RepeatedTiming[v2 = serimplementation[mat, subsize];]
SameQ[v1, v2]
(* {0.0871, Null} *)
(* {0.016, Null} *)
(* True *)

• I implemented the edit in the OP several ways, including using SeriesCoefficient,, Series operations, Coefficient, etc. Overall, I found the simple Coefficient the best performer over various scenarios. Be aware, you need to be careful timing/comparing coefficient related functionality - MMA does some sneaky caching (that is not cleared by ClearSystemCache, etc., and the only way I've found to get reliable timing is to do single tests, one at a time, on a clean kernel. +1 in any case. – ciao Oct 4 '16 at 22:28
• I agree with your assessment that Coefficient is in fact faster when tested properly. However, I think that using a new random matrix should be sufficient to provide a fair test (if there were any matrix independent caching, we probably shouldn't count that anyway). – mikado Oct 5 '16 at 20:07

Here is an implementation based on ListConvolve

lcimplementation[mat_, n_] :=
With[{c = ConstantArray[1, Length[First[mat]]]},
Last[Fold[
Take[ListConvolve[{c, #2}, #1, {1, -1}, 0, Times, Plus, 1],
UpTo[n + 1]] &, {c}, mat]]]


Compare this with the reference implementation

refimplementation =
Coefficient[Times @@ (1 + #1 x) // Expand, x^#2] &;


For timing tests use a new random matrix (and see that the same result is given)

matsize = 100;
subsize = 30;
testsize = 10000;
With[{mat = RandomInteger[{-10, 10}, {matsize, testsize}]},
{AbsoluteTiming[v1 = refimplementation[mat, subsize];],
AbsoluteTiming[v2 = lcimplementation[mat, subsize];]}]
v1 === v2
(* {{8.90806, Null}, {7.37183, Null}} *)
(* True *)


We see a small improvement in speed. However, this approach can benefit from working at machine precision. We see (in this case at least) it locates the same maximum, though there is a risk (for two nearly equal peaks) that it will get it wrong.

matsize = 100;
subsize = 30;
testsize = 10000;
With[{mat = RandomInteger[{-10, 10}, {matsize, testsize}]},
{AbsoluteTiming[
v1 = refimplementation[mat, subsize] // Position[#, Max@#] &;],
AbsoluteTiming[
v2 = lcimplementation[N[mat], subsize] // Position[#, Max@#] &;]}]
v1 === v2
(* {{9.1465, Null}, {2.14473, Null}} *)
(* True *)


On my computer, we have a speed improvement of more than 4

• +1 - I'd not even considered using convolution directly. – ciao Oct 8 '16 at 22:25