Creating a fast random array generator

Question

I want to create N x N array of random numbers. (N is an integer). The constraints are, each element should be between 0 to l(integer) ( here 0<= l <= N) and the sum of elements in each row should be less than equal to N All elements are integers.

I have to call this random array generator function many times ( e.g. 10^6 ) and thus this has to be fast.

I wrote the code with procedural programming (very slow), but looking for a way to write it without loops.

For e.g. N=10, l=8, I want to generate 10x10 array with each element between 0 to 8 such that each row sums up to 10. I want to generate this array a million times.

My Code:

Nport = 10;
l = 8;
y = {};
While[True,
  x = RandomChoice[Range[0, l], {10^5, Nport}];
 AppendTo[y, Select[x, Total[#] <= Nport &]];
  If [Length[y] >= Nport, Break[]]
 ] // AbsoluteTiming
Take[Flatten[y, 1], Nport]

Please Help

I am using Mathematica 9

Answer 1

I'm going to just scrap my previous answer and hopefully this will actually work to include your parameter l, which I'm writing L to avoid confusion (whereas your N is n to avoid syntax errors).

MyNorm[v_] := n v/Total[v]

MakeRow[] := Block[{v},
 v = MyNorm[RandomReal[1, {n}]];
 While[Max[v] > L, v = MyNorm[RandomReal[1, {n}]]];
 v
];

MakeArray[] := Table[MakeRow[], {n}]

n = 50; L = 5;
Timing[For[i = 1, i <= 100, i++, MakeArray[]]][[1]]/100
(* Output: 0.00062400 *)

So for $n=50$ and $l=5$, we can get these cooked up relatively quickly, so that $10^6$ of them will take 10 minutes. Notice that $l$ being close to 1 is a serious problem, so making $l$ smaller will definitely make this worse (in terms of speed):

n = 50; L = 1.6;
Timing[For[i = 1, i <= 100, i++, MakeArray[]]][[1]]/100
(* Output: 0.21372137 *)

We can test this out to make sure it works:

Total /@ MakeArray[]
(* Output: {50., 50., .... , 50.}

Max /@ MakeArray[]
(* Output: {1.58351, 1.54753, 1.57146, .... , 1.59688}

However, if $l\sim N$, there shouldn't be much of a slowdown. For example:

n = 1000; L = 100;
Timing[For[i = 1, i <= 100, i++, MakeArray[]]][[1]]/100
(* Output: 0.02948419 *)

According to my tests, this is as fast as David's answer, but of course includes the parameter L. Generally if L is not too small, this should be as fast as David's, I think.

This means that, on my PC anyway, you'll get 10^6 of these (when $n=1000$ and $L$ is not too small) over the course of a couple days. Not great, but not infeasible if this is a serious project that's going to be using that large a timescale. On the other hand, $1000\times 1000$ is quite big.. no idea if you want them that big.

Any input for us on exactly how big $n$ is going to be? What is the timeframe you are thinking about? What else is taking up computing time in your algorithm? If it's really $n=10$ and $L=8$, then you can generate $10^6$ of these in a minute or two.

EDIT FOR INTEGER VALUES

If you want integers in these arrays, the question is substantially different. I don't think you're going to get much faster than this:

Needs["Combinatorica`"];
n = 32; L = 20;
MakeRow[] := Block[{v = RandomComposition[n, n]},
  While[Max[v] > L,
   v = RandomComposition[n, n]
   ];
  v
  ]
MakeArray[] := Table[MakeRow[], {n}]

10^3 Timing[For[i = 1, i <= 10^3, i++, MakeArray[]]][[1]]

The output there indicates that generating $10^6$ of these will take about 45 minutes. I don't think you can beat RandomComposition, and I doubt even more strongly that you could bring this down to the efficiency of the non-integer version. Maybe this 45 minutes isn't best, but it's not going to be remotely close to 5 minutes (or 1 minute) without a very clever implementation.

Answer 2

I have scrapped this answer -- please consider my other answer instead. I'm leaving this here for reference.

You should simply use Mathematica's built-in random number generator. For example:

n = 50;
f[] := (n /Plus @@ #) # & /@ & /@ Table[Random[], {i, 1, n}, {j, 1, n}];

This function (with no argument) f[] will return such a matrix. For n=50, we can test the timing.

Timing[For[i = 1, i <= 10^4, i++, f[]]][[1]]/10^4
(* Output: 0.0013197685 *)

The average instance of $f[]$ for me takes about 0.0013s. So at scale, if you had $10^6$ of these and $n=50$, it's somewhat feasible (about 21 minutes total to generate them all).

It might be useful to know what your n (or N) value is. If n=5, this is a much different question than n=1000, in terms of how practically we can produce $10^6$ of these $n\times n$ arrays.

I think both of the answers here have still misconstrued your use of l however. I will be updating this answer if I can incorporate that constraint into this.

(Edited due to misreading the question.)

Answer 3

A stylistic suggestion: Never use the variable $N$, which is a function in Mathematica that gives the real value of an expression.

n = 10^3;
ii = 174;
n (#/Total[#] & /@ Table[Min[RandomReal[],ii/n], {n}, {n}])

This should be fast enough:

n = 10^3;
Timing[n (#/Total[#] & /@ Table[Min[RandomReal[],ii/n], {n}, {n}]);]

(* {0.451378, Null} *)

Answer 4

I'm not sure what distribution is desired. Here is one that chooses uniformly among all distinction permutations of partitions of all integers m <= n into nonnegative parts no greater than k. This is done by transposing the Young tableaux for partitions of 2n into at most k + 1 parts, and subtracting 1 to get the parts to be between 0 and k. We then have to select those whose sums are at most n.

ClearAll[getparts, weights];
mem : getparts[n_, k_] := mem = getpartsC[
    PadRight[IntegerPartitions[2 n, k + 1], {Automatic, k + 1}],
    n];
getpartsC = Compile[{{partitions, _Integer, 2}, {n, _Integer}},
   Select[
    Function[p, Total@Transpose[UnitStep[p - #] & /@ Range[n]] - 1] /@ partitions,
    Last[#] >= 0 &]
   ];
mem : weights[n_, k_] := mem = Multinomial @@@ (Tally /@ getparts[n, k])[[All, All, 2]];
randsamp = Function[{n, howmany}, Ordering /@ RandomReal[1, {howmany, n}]];

ClearAll[nparts];
nparts[n_, k_, nsamp_] :=
  With[{parts = getparts[n, k]},
   With[{p = RandomChoice[weights[n, k] -> parts, nsamp]},
    Compile[{{pp, _Integer, 2}, {samp, _Integer, 2}},
      MapThread[#1[[#2]] &, {pp, samp}]][p, randsamp[n, nsamp]]
    ]];

A one-time cost: computing the basic partitions:

getparts[32, 20] // Dimensions // AbsoluteTiming
(*  {32.0588, {43202, 32}}  *)

Generating samples is fairly quick (it seems to be an order of magnitude faster than the Combinatorica` code of Kellen Meyers, although there is a claim in a comment, without supporting code, that suggests it might be about the same). A million samples should take a total of about six minutes.

Table[nparts[32, 20, 32], {10^3}] // Dimensions // AbsoluteTiming
(*  {0.35359, {1000, 32, 32}}  *)

Getting them in chunks is faster, if you have the memory. A thousand 32x32 arrays is not much, but a million takes more memory than I have RAM, and therefore would be much slower to do all at once.

ArrayReshape[
   nparts[32, 20, 32*10^3],
   {10^3, 32, 32}] // Dimensions // AbsoluteTiming
(*  {0.070147, {1000, 32, 32}}  *)

---

Comparison of this approach and the Combinatorica` approach (MakeArray):

10^3 AbsoluteTiming[Do[MakeArray[], {10^3}]]
10^3 AbsoluteTiming[Do[nparts[32, 20, 32], {10^3}]]
10^3 AbsoluteTiming[ArrayReshape[nparts[32, 20, 32*10^3], {10^3, 32, 32}];]
(*
  {3146.06, 1000 Null}
  {357.585, 1000 Null}
  {78.877, 1000 Null}
*)

Note that you need to add 30 sec. to the nparts results for the one-time computation of the basic partitions.