0
Not a bug

randn() precision limited to 55 digits

dattorro 3 months ago updated 3 weeks ago 22

%randn ignores high precision

mp.Digits(1034)

t=randn(1,'mp')
t =0.1048747160164939795645722142580780200660228729248046875

t - 0.1048747160164939795645722142580780200660228729248046875
ans = 0

MATLAB Version 9.14 (R2023a)
Advanpix Multiprecision Computing Toolbox Version 5.2.8.15537

Answer

Answer
Not a bug

The rand(...,'mp') generates numbers in double precision and then converts them to mp. Hence digits are not filled up to the 1034th digit. We do this in order to produce the same random numbers as built-in MATLAB's rand function (for reproducibility).

If you want random numbers to the full precision (all 1034 digits are filled), please use mprand/mprandn.

>> mp.Digits(100);
>> randn(1,'mp')
ans =
0.8621733203681205548463140075909905135631561279296875

>> mprandn(1,'mp')
ans =
0.1868724737876430395187949415904261118721051078748075419355524734700877850288747581541290481709356331
Answer
Not a bug

The rand(...,'mp') generates numbers in double precision and then converts them to mp. Hence digits are not filled up to the 1034th digit. We do this in order to produce the same random numbers as built-in MATLAB's rand function (for reproducibility).

If you want random numbers to the full precision (all 1034 digits are filled), please use mprand/mprandn.

>> mp.Digits(100);
>> randn(1,'mp')
ans =
0.8621733203681205548463140075909905135631561279296875

>> mprandn(1,'mp')
ans =
0.1868724737876430395187949415904261118721051078748075419355524734700877850288747581541290481709356331

mprand() and mprandn() appear nowhere in the documentation nor on the Advanpix site;

Google search is empty: mprandn site:advanpix.com

Syntactically, function calls to mprandn(n) are expected equivalent to randn(n,'mp').

Sure looks like a bug to me.

Agree, this is lack of documentation. The mprand(n) was an experimental name, I just left it for the time being.

Please suggest other name, unambiguous enough. Thank you.

I just found that the mprand is not connected with rng, so there is no way how to control random number generator mprand(n)? How to setup user defined seed for mprand(n)?

From my point of view, the following two commands are now fully equivalent: mp(rand(1)) and rand(1,'mp') and final random number is effectively the same, and in double precision.

So I think, that rand(1,'mp') should represent random number in full precision defined by mp.Digits(100) by default, and be equivalent to your experimental function mprand(n).

But in the 'mp' (full precision) case should be possible to control random generator by rng command, too.


So, the suggested way is elimination of mprand(n) at all, and use the rand(1,'mp') at full precision
regime.

+1

Thank you very much for sharing your opinion. Yes, this is the way. This only difficulty is how to monitor when user calls "rng" to set the new seed. The easiest way would be to overload the 'rng' and make it part of toolbox. Or keep track of the 'rng' state internally in toolbox.... In any case, will include this into next toolbox version.

I have just released new version with this feature - 5.3.4.15817. Please update your toolbox.

The "rng" still not working with rand(1,'mp'):

>> mp.Digits(100)

ans =

300

>> rng(1)
rmp = rand(1,'mp')
rng(1)
rmp = rand(1,'mp')

rmp =

0.1019434353836572895054626158415646843010934355718609316510899011279499964199600542084573004435709827


rmp =

0.7686462842760086723268602992488197377426292102713405281965513370345202336516897460693092237330618525

Good catch!. It is not working only if the seed is the same. Now toolbox is monitoring the case when seed is changing. This is kind of light approach - without deep intrusions to MATLAB core.

To reset the sequence for the same seed would require full-blown hook over 'rng' and actually monitoring its calls. This will make it slow(er). 

Especially for me, it is important to have full control over random number generation.

Could you add some clear and simple example how "rng" works with rand(x,'mp') now? Or better, could you describe this specific behavior at function reference WWW page?

And finally, what about rng('default') and rng('shuffle'), these two options are extremely useful? 

Current situation. If user sets new seed with 'rng' (using rng(seed), rng('default') or rng('shuffle')) - then toolbox resets the state of its random number generator to this new seed. The 'new seed' means different seed from what was previously set.

In your example, first call rng(1) resets the state of generator. The second call rng(1) doesn't trigger reset because the seed is the same. This is the bug and needs to be fixed (working on this). 

  

Here is generator reset when seeds are different:

>> rng(0);
>> randn(2,1,'mp')
ans = 
        0.1868724736939647133172833335569502
        0.8829593750178176873331742398432483
 
>> rng(1);
>> randn(2,1,'mp')
ans = 
        -0.06238805020702318459649989108952415
        -0.51254018310017212737352467081772230
 
>> rng(0);
>> randn(2,1,'mp')
ans = 
        0.1868724736939647133172833335569502    
        0.8829593750178176873331742398432483    
 
>> rng(1);
>> randn(2,1,'mp')
ans = 
        -0.06238805020702318459649989108952415
        -0.51254018310017212737352467081772230

In the next version, the toolbox will work with the rng the same as built-in rand/randn (call to rng with the same seed will trigger the generator reset).

Well, here is update after few hours of trying to monitor/overload 'rng' function calls in MATLAB.

As it turns out, there is no such a way :(


 If, say, 'rng' is overloaded, then it is impossible to call the native MATLAB's function 'rng' anymore (from the overloaded). Which basically prevents us from monitoring each call to 'rng' = impossible to track each call to 'rng' to handle the case when seed is the same.

So there are two ways possible (please decide what you like better):

1. Same as now. Toolbox relies on the seed from 'rng', but resets state only iff the seed is different (see previous comment for details).

2. New, special function is introduced to reset the state of random number generator in toolbox. Say mp.ResetRandom.
Then it can be used as: 

mp.ResetRandom(seed)
mp.ResetRandom('default')
mp.ResetRandom('shuffle')

So that user needs to use mp.ResetRandom if he/she wants to reset the full-precision rand/randn(...,'mp').

Let me know what you prefer.

I think mp.ResetRandom is the way to go, if we want a full control.

This is a hard question, because for a regular MCT user like me, it is not clear how difficult it is to create fully functional mp.ResetRandom function as an full-precision alternative to the "rng" function. Moreover, full-precision random generators (rand(n) + rng ... 'mp' functions) are so specific, that it is not clear how important these functions really are for other MCT users in general.


From my point of view, is the new "mp.ResetRandom" function optimal solution, I just prefer better name: like "mp.rng" (if is it possible, of course, due to the mentioned overloading problem), to be close to standard MATLAB as possible. I am working now on the full-precision random sampling, where is this topic crucial, including the 'mp' alternative of "rng" to get full control over random generation process.

+1

I have just released new toolbox version (5.3.5.15874) with new function mp.RandState. It allows precise control of the random generator state. Its usage syntax and semantic is the same as for built-in 'rng', e.g.: 


mp.RandState(seed)
mp.RandState('default')
mp.RandState('shuffle')
mp.RandState(seed,'twister')
mp.RandState(state)


The mp.RandState returns full state of the generator, therefore it is exactly reproducible.

Simple example can be found here: https://www.advanpix.com/documentation/version-history/  

Just a perfect, Pavel!!! Thanks...:)
Final request: Please add a description of the mp.RandState function to the Web function reference as a basic MCT functions info source.

Thanks again ...

clear all

clc

mp_angle = 2*rand(1,'mp') - 1

dp_angle = 2*rand(1) - 1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Run this program multiple times.

mp_angle is always the same.

dp_angle is always different.

Cause of bad mp_angle is `clear all' statement.

Despite Mathworks' recommendation against `clear all' usage, 

there are circumstances where it becomes necessary.

MCT must continue to work exactly the same as Matlab does.

Otherwise, every MCT program ever written with `clear all' is broken.

The 'clear all' wields full power over the MEX - as it just kills all the MEX modules (including mp-toolbox).
The mp.RandState is provided exactly for such cases when precise and explicit control over generator is needed. Please use it.

Matlab’s built-in random number generator never treats “clear all” as a command to reset or reseed its generator.  That behavior has also been integral to all prior MCT versions, which aligns with user expectations and ensures backward compatibility.  This sudden unexpected change without advance notice, making “clear all” reinitialize the MCT random seed, wields significant deleterious impact:

  1. Breaks Established Code: Existing scripts rely on “clear all”, merely to remove variables and free memory, with the understanding that random seeds remain untouched.  Altering that assumption invalidates all results that previously depended on dynamic seeds.

  2. Conflicts with Matlab Conventions: Users familiar with Matlab expect “clear all” to leave the random seed alone. Departure from that convention forces MCT users to incorporate extra & unexpected steps (like manually restoring the seed) whenever “clear all” appears in a script.

  3. Impacts Reproducibility: Research code often depends on a fixed seed for reproducible Monte Carlo or simulation results. Changing that seed unintentionally, simply because “clear all” appears for housekeeping, undermines reproducibility and can lead to subtle or undetectable failures, erroneous or repetitive results, and false conclusions..

  4. `Cause Unknown': Even experienced users will not know to remove "clear all" or to invoke the new mp.RandState function. Hours of needless debug time per user inevitably ensue.

While mp.RandState offers a manual workaround, it introduces overhead and potential for confusion.  Restoring the original MCT policy —so that MCT behaves like Matlab in not resetting the random generator after “clear all”— keeps MCT code compatible, predictable, and consistent with core Matlab behavior.

@"Matlab’s built-in random number generator never treats “clear all” as a command to reset or reseed its generator."

MATLAB doesn't notify the toolboxes (nor MEX modules) when user calls 'clear' function on it. So that toolbox has no way to react to this event in some special way (e.g. to store/pass the seed to the next call). The 'clear' just kills the toolbox from executing, unloads and clears its memory. Next call to toolbox starts it fresh from default state.

Please send MT request to allow flexible handling of 'clear' function in 3rd party toolboxes/MEX.

Now toolbox provides two ways of using random number generators:

(1) Using the standard MATLAB functions, e.g. mp(rand(..)). The seed functionality is controlled by MATLAB. But precision is limited to double precision numbers.

(2) The generation of a full-precision random numbers requires special syntax (rand(..,'mp')) and a separate seed control  function (mp.RandState). The separate seed control function is needed because of the MATLAB's restrictions on overloading of 'rng' and 'clear' functions (it is not possible to do or there is no notification mechanism provided to handle the events).

User can choose the (1) or the (2). In case of (2) manual control of a seed is required using the special function mp.RandState. 

Option (3): MCT reads/writes current state/seed from/to disk whenever  rand(n,m,'mp')  or  randn(n,m,'mp) is called.


Rationale: Spinning rust disk drives are obsolete.  

High-end Samsung PRO NVMe M.2 drives can reach on the order of 500,000 to 1,000,000 random write IOPS (4 KB) in burst conditions.  In contrast, even the fastest 15,000 rpm HGST enterprise hard disks typically peak around 300 to 400 random write IOPS.  In other words, the Samsung PRO M.2 SSD’s best-case random burst write performance is roughly three orders of magnitude higher than that of the best rotating magnetic media.

Advantage:  mp.RandState remains useful but no longer compulsory.

This is not an option either. Putting the slowness aside (applying the seed/state is a non-trivial process), there is another reason:

Random state must be reset to default on MATLAB startup (when user starts up the MATLAB/starts work with toolbox).

MATLAB doesn't allow differentiating startup vs. 'clear all' events. These events must be treated differently, but we can choose only one. We choose the most flexible: "startup" behavior and allow the user to take care of random state as one wishes.  User just need to call mp.RandState when needed (as it was designed).

This doesn't lead to any slowdowns, undocumented hacks/hooks to MATLAB kernel, etc. etc. etc.

Update: I edited the comment after posting it (to remove logically incorrect statement).

Perhaps all present and past users of MCT should be notified (by direct email) that intrinsic behavior of random functions has changed with regard to: 1) Matlab `clear all' command, 2) extension from 55 digits to arbitrary precision, and 3) whether MCT randomization function algorithms now differ from Matlab.