8
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ 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]]. $\endgroup$ – prog9910 May 31 at 23:39
7
$\begingroup$

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
| improve this answer | |
$\endgroup$
  • $\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 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 at 2:05
4
$\begingroup$

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} *)
| improve this answer | |
$\endgroup$
  • 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 at 1:21
  • $\begingroup$ I still like to see Map[ opposed to /@, but that is more. $\endgroup$ – prog9910 Jun 1 at 2:08
4
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\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 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 at 1:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.