# How can I replace a set of Nested For Loops with a Functional Code Alternative?

Despite having looked over various posts concerning optimizing nested For loop expressions into a functional programming equivalent I still can not figure out how to actually to do this. Such constructs are inefficient in Mathematica and hence, I seek to find a more efficient approach. Sadly, having started with programming with FORTRAN (prior to IV), some ideas die hard in an aging brain.

I have a very large array that I must evaluate. This array m, is 31x2754, but for purposes of example I have down-sampled so that it is only a 30 x 32 matrix. It produces a square matrix whose order is the size of the columns (here 32, when downsampled, but 2754x2754 when fully sampled). Although a doubly nested set of For loops provides the answer I seek and is relatively efficient when the matrix is only 32x32, it's efficiency degrades rapidly as n approaches 2754. I am at a loss to how to convert the following procedural code snippet into a functional programming equivalent that is more efficient.

The test input (ignoring how I get there, which itself is rather complicated) is:

m = {{0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1,
0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1}, {1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0,
0}, {1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0,
0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 1, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 1, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 1, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1,
0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0}, {0, 1, 1, 0, 1, 1, 1, 0, 0, 0,
0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0,
0}, {0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1,
1, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0}, {0, 1, 0, 0, 0, 1, 1, 0, 0, 0,
0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1,
0}, {0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1,
0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 1, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1,
0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 1, 0, 0, 0, 0, 1, 1,
1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1}}

Dimensions[m]

{30, 32}


Although not essential, to be sure that the output is properly created I create a constant array and fill the matrix c, which will be the eventual output matrix with a value that WILL NOT be present in the matrix upon completion, so that by quick observation I can be sure that the matrix has been properly populated.

This is done as follows and establishes the order of the square matrix, where the length of columnlabels for this snippet is 32 (the actual column labels are irrelevant here):

 c = ConstantArray[3, {Length[columnlabels], Length[columnlabels]}];

{characterrows, charactercolumns} = Dimensions[c]


The loop construct I wish to replace with a functional approach effectively compares two column of m, each of which has only two distinct values (either 0 or 1), and for which their product partial order potentially generates 4 distinct possibilities [vertices] ({0,0},{1,0},{0,1},{1,1}} and for which I want to assign to the each pair of columns a 1 if less than 4 of these possibilities are actually present in the comparison and a 0 if all 4 possibilities are present when comparing the two columns). This construct is as follows:

 t0 = AbsoluteTime[];

For[i = 1, i <= Length[columnlabels], i++,
For[j = 1, j <= Length[columnlabels], j++,
f = m[[All, i]];
g = m[[All, j]];
c[[i, j]] = If[Length[Tally[Table[{g[[n]], f[[n]]}, {n, 1, Length[g]}]]] < 4, 1, 0];
]
]

t1 = AbsoluteTime[];
timeelapsed = UnitConvert[Quantity[t1 - t0, "Seconds"], "Minutes"]


What is the most efficient functional programming expression needed to replace this nested set of For loops?

• What if? You could show what you are asking, with a 32-element vector. ie. Let's try to simplify the work. (+) Otherwise, your code looks like a double loop. Which functional code means Nest[ within Nest[. (+) Then, you are performing a conditional function in the long line of c[[i,j]]. May 31, 2020 at 23:39

Update: A faster alternative:

foo = Boole @* LessThan[4] @* Length @* Union @* Transpose @*
DeveloperToPackedArray @* List;

c0 = Outer[foo, mt, mt, 1];

c0 == c

True


mt = Transpose[m];
c1 = Boole @ Outer[Length@Union@Transpose[{##}] < 4 &, mt, mt, 1];

c1 == c

True


Alternatively,

mt = Transpose[m];
c2 = ConstantArray[0, {32, 32}];
Do[c2[[i, j]] = c2[[j, i]] = Boole[Length@Union[Transpose[{mt[[i]], mt[[j]]}]] < 4],
{i, 1, Length @ columnlabels}, {j, 1, i}];

c2 == c

True


and

mt = Transpose[m];
c3 = SparseArray[{i_, j_} :>
Boole[Length@Union[Transpose[{mt[[ i]], mt[[j]]}]] < 4],
{1, 1} Length[columnlabels]]

Normal[c3] == c

True


and

mt = Transpose[m];
c4 = SymmetrizedArray[{i_, j_} :>
Boole[Length@Union@Transpose[{mt[[ i]], mt[[j]]}] < 4],
{1, 1} Length[columnlabels], Symmetric[{1, 2}]]

Normal[c4] == c

 True

• Marvelously useful, as always. Thanks. c1 seems to run the fastest @ .1984701 on my system, with c2, c3, c4 coming in at 0.2423850, 1.2566402, and 0.3739999, including transposing m to get mt in each.. I find c3 and c4, albeit a bit slower, have the most intriguing code, although the syntax is more difficult to digest as the expression segment defining the Held SymetrizedArray {1,1} Length[columnlabels] remains mysterious to me. Is it that Outer is more efficiently compiled that gives it the edge? Jun 1, 2020 at 1:21
• @StuartPoss, thank you for the accept. I don't know why Outer is faster. I guess the method with SymmetrizedArray is faster than SparseArray because (somehow taking advantage of symmetry) some calculations are avoided.
– kglr
Jun 1, 2020 at 2:05

An improved version of kglr's answer, makes use of the fact that m only consists of 0 and 1:

m = RandomInteger[{0, 1}, {31, 2754}];

mt = Transpose[m];

func = Composition[Length, Union, Plus];

c2 = 1 - (Outer[func, mt, 2 mt, 1] - 4 // UnitStep); // AbsoluteTiming

(* {22.1601, Null} *)


kglr's solution takes about 53 seconds. Tested on v12.1, Wolfram Cloud.

Remark

My solution is slower in v9.0.1. (72 seconds v.s. 39 seconds. ) Not sure about the reason.

# Update

A solution with Compile (fastest one so far):

help = Compile[{{mat, _Integer, 2}},
Table[If[4 > (lsti + 2 lstj // Union // Length), 1, 0], {lsti, mat}, {lstj, mat}](* ,
CompilationTarget -> C *)]

test = help@mt; // AbsoluteTiming
(* {9.29816, Null} *)


If you have a C compiler installed, add the CompilationTarget -> C option and the code will be faster.

P.S.

I didn't expect ContainsAll/SubsetQ is so slow.

• (as a non-Mathematica user) if there's a thing called Compile, why doesn't it also fix his FORTRAN-style nested loop? Jun 1, 2020 at 19:42
• @Ben Er… Are you asking whether OP's original approach can be Compiled or not? If so, the answer is yes, with a little adjustion: opfunc=Compile[{{m, _Integer, 2}}, Module[{col = Length@m[[1]],i,j,f,g,c = Table[0,{col},{col}]}, For[i=1,i<=col,i++,For[j = 1,j <= col,j++,f=m[[All, i]]; g=m[[All,j]];c[[i,j]]=If[Length[Union[Table[{g[[n]],f[[n]]},{n,1,Length[g]}]]] < 4, 1, 0];]]; c]]; optest = opfunc@m; // AbsoluteTiming (Notice Union has replaced Tally, and variables are localized. These are necessary for full compilation. ) The timing is 49.5249 s on Wolfram Cloud. Jun 2, 2020 at 1:46

Using Table instead of For.

SeedRandom[99]
dat1 = Table[Table[RandomInteger[8], {6}], {6}];
dat2 = Table[Tally[dat1[[i]][[All]]], {i, 1, 6}];

wrapFn[x_List, yLimit_Integer] :=
If[Length[x] < yLimit, 100, 0]

datOut = Table[wrapFn[dat2[[i]], 5], {i, 1, 6}]
(* Out: {100, 0, 100, 100, 100, 0} *)


Just if it helps the OP, with a pre-Fortran IV mind like myself? The O[?] is probably worse than the checked answer.

New Way:

SeedRandom[99]
dat3 = RandomInteger[8, {6, 6}]
dat4 = Map[Tally, dat3]
dat5Out = Map[If[Length[#] < 5, 100, 0] &, dat4]
(* Out: {100,0,100,100,100,0} *)

• Your first two expressions really should be dat1 = RandomInteger[8, {6, 6}]; dat2 = Tally /@ dat1. If you find yourself using nested Tables, remember that you can use a single one, with multiple iterators. Then, if you find yourself writing Table[function[ something[[i]] ], {i, 1, Length[something]}], i.e. iterating the same operation on each element of something, then consider replacing that with a Map` operation: that's exactly what they are for. Jun 1, 2020 at 1:21
• I still like to see Map[ opposed to /@, but that is more. Jun 1, 2020 at 2:08