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

Question

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]]. — prog9910, May 31, 2020 at 23:39

kglr · Accepted Answer · 2020-06-01 01:49:28Z

8

Update: A faster alternative:

foo = Boole @* LessThan[4] @* Length @* Union @* Transpose @* 
   Developer`ToPackedArray @* List;

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

c0 == c

True

Original answer:

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

edited Jun 1, 2020 at 1:49

answered Jun 1, 2020 at 0:05

kglr

395k18 gold badges479 silver badges902 bronze badges

$\begingroup$ 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? $\endgroup$
– Stuart Poss
Jun 1, 2020 at 1:21
$\begingroup$ @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. $\endgroup$
– kglr
Jun 1, 2020 at 2:05

Add a comment |

xzczd · Accepted Answer · 2020-06-02 02:08:00Z

5

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.

edited Jun 2, 2020 at 2:08

answered Jun 1, 2020 at 5:03

xzczd

66.2k9 gold badges163 silver badges469 bronze badges

$\begingroup$ (as a non-Mathematica user) if there's a thing called Compile, why doesn't it also fix his FORTRAN-style nested loop? $\endgroup$
– Ben Jackson
Jun 1, 2020 at 19:42
1

$\begingroup$ @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. $\endgroup$
– xzczd
Jun 2, 2020 at 1:46

Add a comment |

prog9910 · Accepted Answer · 2020-06-01 02:05:48Z

4

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} *)

edited Jun 1, 2020 at 2:05

answered Jun 1, 2020 at 1:04

prog9910

6304 silver badges10 bronze badges

1

$\begingroup$ 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. $\endgroup$
– MarcoB
Jun 1, 2020 at 1:21
$\begingroup$ I still like to see Map[ opposed to /@, but that is more. $\endgroup$
– prog9910
Jun 1, 2020 at 2:08

Add a comment |

Stack Exchange Network

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

3 Answers 3

Update

Your Answer

Not the answer you're looking for? Browse other questions tagged
performance-tuning
functional-style
procedural-programming
or ask your own question.

Hot Network Questions

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

3 Answers 3

Update

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged performance-tuningfunctional-styleprocedural-programming or ask your own question.

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
performance-tuning
functional-style
procedural-programming
or ask your own question.