4
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ This seems like a duplicate of: (54448) $\endgroup$ – Mr.Wizard Feb 3 '15 at 8:05
  • $\begingroup$ @Mr.Wizard, at first I though so too, but this one seems to have an additional constraint on the values (0<=x<=l) which is not in 54448. $\endgroup$ – Simon Woods Feb 3 '15 at 9:49
5
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ This method is very fast. Thank you. But I am looking for l integer. Check the code attached in my edited question. $\endgroup$ – akhileshsk Feb 3 '15 at 3:28
  • $\begingroup$ How is that a problem? Just because I used $L=1.6$ (as well as several integer values of $L$ like $8$ and $100$) doesn't mean the code will somehow fail if you never use non-integers. This code works just fine for all integer values of $L$... it just happens to work for any other value of $L$ too. $\endgroup$ – Kellen Myers Feb 3 '15 at 3:30
  • $\begingroup$ I will be happy if I can generate 10^6 in a minute for this particular configuration : n = 32 and l= 20 to 30 . Each element of the matrix has to be an integer. $\endgroup$ – akhileshsk Feb 3 '15 at 3:31
  • $\begingroup$ "Each element of the matrix has to be an integer." You're kidding right? You didn't say that at all. That makes this a completely different question. $\endgroup$ – Kellen Myers Feb 3 '15 at 3:33
  • $\begingroup$ Sorry, Just realized that. I said l is an integer , but that doesn't imply that all elements will be integers. $\endgroup$ – akhileshsk Feb 3 '15 at 3:34
1
$\begingroup$

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.)

$\endgroup$
  • $\begingroup$ This will not guarantee that "elements in each row add to N." $\endgroup$ – David G. Stork Feb 3 '15 at 2:39
  • $\begingroup$ I clearly read the question too quickly! $\endgroup$ – Kellen Myers Feb 3 '15 at 2:42
  • $\begingroup$ And not to be rude, but neither does yours! $\endgroup$ – Kellen Myers Feb 3 '15 at 2:49
  • $\begingroup$ Yep. Fixed now. $\endgroup$ – David G. Stork Feb 3 '15 at 2:54
1
$\begingroup$

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

$\endgroup$
  • $\begingroup$ But you'd need to normalize with respect to the sum of elements, i.e. Normalize[#, Total[#]&/n]& or something like that. This still doesn't guarantee that all numbers are less than l, though that is weird constraint and I'm not sure the OP thought it through. Also, you could use RandomReal[1, {n,n}] for brevity and performance. $\endgroup$ – Szabolcs Feb 3 '15 at 2:46
  • $\begingroup$ Okay, +1 now :-) $\endgroup$ – Szabolcs Feb 3 '15 at 2:58
  • $\begingroup$ Any thoughts on the parameter l? $\endgroup$ – Kellen Myers Feb 3 '15 at 3:41
  • $\begingroup$ Incorporated now. $\endgroup$ – David G. Stork Feb 3 '15 at 4:22
1
$\begingroup$

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.

$\endgroup$

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.