Generating 250 random points crashes the kernel, but not 249

Question

Fixed in 10.1.0.

Consider the following function, which generates uniformly random points on the surface of the 2-sphere:

randSphere[] := Block[{z = RandomReal[{-1, 1}, 3]},
    If[Total[z^2] > 1, randSphere[], Normalize[z]]]

I can use this function to generate a Table of 249 points:

Table[randSphere[], {249}] (* works fine *)

but mysteriously, changing 249 to 250 consistently crashes the kernel. I am running Mathematica 10.0.2 on Windows. What's going on here? It's worth noting that I can also generate 249 pairs of points with no problems:

Table[{randSphere[], randSphere[]}, {249}] (* also works fine *)

and I can even generate 249 Tables of 249 points:

Table[Table[randSphere[], {249}], {249}] (* still fine *)

but changing any instance of 249 to 250 in each of the above examples crashes the kernel again.

Reproduced in 10.0.2 OS X. Table will attempt to auto-compile its first argument starting at 250, which is certainly related to the crash. The CompileOptions setting from SystemOptions controls this. — Szabolcs, Dec 30 '14 at 18:13
Reproduced in v9.0.1 as well. This crash is almost certainly related to the recursive call in randSphere. The workaround is to get rid of that (use a While[... > 1, (* recompute *)] sort of thing). Compile doesn't support recursive calls. I'm tagging this as a bug, please report it to support at wolfram.com. — Szabolcs, Dec 30 '14 at 18:15
@Szabolcs Interesting. I've never heard of this compiling behavior before; is it documented anywhere? The docs for Table don't mention anything special about 250. — David Zhang, Dec 30 '14 at 18:16
It probably isn't documented, I'm not sure. It's often discussed on Mathematica forums though, so it's "common knowledge" in the community. Try SetSystemOptions["CompileOptions" -> "TableCompileLength" -> 300] to set it to 300. I don't mean this as a workaround though. — Szabolcs, Dec 30 '14 at 18:21
You can get all default compile lengths by evaluating SystemOptions["CompileOptions"]. — Karsten7, Dec 30 '14 at 18:25
For other methods of generating random points on a sphere, you may be interested in this Q&A: http://mathematica.stackexchange.com/questions/13038/vectors-with-a-certain-magnitude-in-mathematica — Michael E2, Dec 30 '14 at 18:56
Note also that Compile[{}, randSphere[]] itself crashes my kernel (V10.0.1 Mac OSX 10.9.5). I'm not sure of the relation to the OP's problem, since autocompile might do a different command. — Michael E2, Dec 30 '14 at 19:09
@MichaelE2 Your simpler example indeed runs afoul of the same underlying problem. — Daniel Lichtblau, Dec 30 '14 at 19:44
@MichaelE2 By the way, your formulation will have a separate but possibly related issue in that the (weak) type inferencing mechanism in Compile will not know what is getting returned. The variant below makes this explicit, and has the pleasant side effect of circumventing the crash.ff = Compile[{}, randSphere[], {{_randSphere,_Real,1}}]; Table[ff[],{250}] I do not understand how this manages to evade the crash. Which is weird, since I now know what causes said crash. One of those mysteries. — Daniel Lichtblau, Dec 30 '14 at 20:57

score 25 · Accepted Answer · edited Dec 31 '14 at 02:07

25

I can reproduce this on OS X in M10.0.2 and M9.0.1, so it looks like a bug. Please report it to Wolfram support.

Table will automatically try to compile its argument above a table length threshold. This threshold is 250 by default and can be set to a different value using SetSystemOptions["CompileOptions" -> "TableCompileLength" -> ...]. It seems the crash happens only when Table compiles its argument.

The randSphere function is recursive but Compile doesn't support recursion. My guess is that the crash is related to this.

I recommend eliminating the recursion as a workaround:

randSphere[] := Module[{z},
  While[True,
   z = RandomReal[{-1, 1}, 3];
   If[Total[z^2] <= 1, Return@Normalize[z]]
   ]
  ]

This version won't crash.

edited Dec 31 '14 at 02:07

halirutan

112,764
7
263
474

answered Dec 30 '14 at 18:26

Szabolcs

234,956
30
623
1,263

5

We are aware of this issue. I don't have an ETA for when it will be fixed, though. – Stefan R Dec 30 '14 at 19:43
1

@StefanR Out of curiosity, is the problem really related to the recursion? – Szabolcs Dec 30 '14 at 20:04
4

Yes, it's related to the recursion. Not the actual running of the recursion, but in the preprocessing. – Daniel Lichtblau Dec 31 '14 at 00:18

score 2 · Answer 2 · answered Jan 01 '15 at 23:22

As an alternative, you can use the following function to generate pseudo-random points on the sphere:

PointOnTwoSphere[] := Module[
  {z, phi, rho, x, y},
  z = RandomReal[{-1, 1}];
  phi = RandomReal[{0, 2*Pi}];
  rho = Sqrt[1 - z^2];
  x = rho * Cos[phi];
  y = rho * Sin[phi];
  {x, y, z}
]

This algorithm, due to Marsaglia, has the virtue that it does not throw away any points.

Chen Stats Yu · Answer 3 · 2015-01-01T23:36:28.650

0

Not a solution but worth noting.

enter image description here

ParallelTable[randSphere[], {250}]

Does seem to work.

I dont know about its limit, but eventually, it crashed (as expected for "big" data). enter image description here

edited Jan 01 '15 at 23:36

answered Jan 01 '15 at 23:30

Chen Stats Yu

4,986
2
24
50

Generating 250 random points crashes the kernel, but not 249

3 Answers3

Linked