I have just checked the source code of the nchoosek in the newest versions of MATLAB.
MathWorks changed it to do all the computations in double. Here is some code from the function: 

% Do the computation in doubles.  
nd = full(double(n)); 
kd = full(double(k));           

I am not sure, but the code might start work fine after replacing all the "double" to "mp".

Looks like MW has no plans to support OOP/custom classes properly at all.

Meanwhile I will try to implement the function in toolbox's core.

Jon, Happy New Year!

We have just released new version 5.3.6.15902: https://bit.ly/3mvxXqg
Now, its least-squares solver generates bitwise deterministic results every time.

Would appreciate it if you would update me on your test results.

Repeated runs of the same problem in the same environment will yield the same bitwise result unless something in the environment changes; e.g, CPU type, number of threads,...

I think with more advanced algorithms and many-core CPUs this will change in the future.

For example. Even though we use the same number of threads in a solver - it might split the overall work in different number  of jobs and distribute them differently among cores on each call. This might happen because some CPU cores might be busy with other processes on the computer. More flexible/adaptive algorithms take this into account (as we did in least-square solver before 5.3.6).

Thank you for all your support and contributions!

Best wishes for the New Year!

Pavel

Jon, thank you very much for the report!

Generally, some degree of variability of the solution is to be expected (due to matrix conditioning, multi-threading, flow of divide-and-conquer/bulge chasing computations, etc.).

The condition number of a matrix is log10(cond(Amat)) = ~21, roughly meaning that last ~21 digits are expected to be lost during computation of the solution.

However, even if we use extended precision, the last ~21 digits of the computed solution are not accurate (=can vary every time we solve the system). What counts is the relative accuracy of the system delivered by the solution. 

I added few changes to your script:

clear all

ndigits=34;
mp.Digits(ndigits);

load vL
vL = mp(vL,ndigits);
Amat=mp(Amat,ndigits);

%vL = double(vL);
%Amat = double(Amat);

n0 = norm(Amat,1)*norm(vL, 1); % PH: norms to compute the relative accuracy

xi = Amat\vL;  %backslash error is nonzero when differenced with same backslash later on

for i=1:10
    notabug = Amat\vL;
    
    acc = norm(Amat*notabug-vL,1) / n0;     % PH: system rel. accuracy w.r.t solution vector.
    xacc = norm(xi - notabug,1)/norm(xi,1); % PH: rel. difference from the first solution.
    [m,idx] = max(abs(xi - notabug));
    
    fprintf("rel.acc(sys) = %.5e, rel.acc(x) = %.5e, [xi(%d) = %.8f, notabug(%d) = %.8f, rel.acc = %.5e]\n",acc,xacc,idx,xi(idx),idx,notabug(idx),abs(xi(idx)-notabug(idx))/abs(xi(idx)));
    
    %sum(abs(xi - notabug))
endv

Here is its results for precisions: 34, 50, 100 and 500 decimal digits.

% ndigits = 34
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 2.54932e-19, [xi(40) = 33835.88114407, notabug(40) = 33835.88114407, rel.acc = 3.10975e-19]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.24177e-10, [xi(41) = 5382.66065303, notabug(41) = 5382.66064798, rel.acc = 9.37861e-10]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 8.72727e-11, [xi(41) = 5382.66065303, notabug(41) = 5382.66065016, rel.acc = 5.33988e-10]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.04572e-10, [xi(41) = 5382.66065303, notabug(41) = 5382.66064912, rel.acc = 7.26616e-10]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 8.19852e-11, [xi(41) = 5382.66065303, notabug(41) = 5382.66064979, rel.acc = 6.01897e-10]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.23477e-10, [xi(41) = 5382.66065303, notabug(41) = 5382.66064801, rel.acc = 9.32217e-10]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.23477e-10, [xi(41) = 5382.66065303, notabug(41) = 5382.66064801, rel.acc = 9.32217e-10]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.24177e-10, [xi(41) = 5382.66065303, notabug(41) = 5382.66064798, rel.acc = 9.37861e-10]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 7.31347e-11, [xi(41) = 5382.66065303, notabug(41) = 5382.66065016, rel.acc = 5.33995e-10]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.24177e-10, [xi(41) = 5382.66065303, notabug(41) = 5382.66064798, rel.acc = 9.37861e-10]

% ndigits = 50
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 3.92138e-30, [xi(40) = 33835.88114392, notabug(40) = 33835.88114392, rel.acc = 3.26593e-30]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 9.09178e-27, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 6.17458e-26]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 5.67413e-27, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 3.50551e-26]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 5.65223e-27, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 3.38773e-26]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 5.21849e-27, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 3.40105e-26]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 3.92138e-30, [xi(40) = 33835.88114392, notabug(40) = 33835.88114392, rel.acc = 3.26593e-30]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 6.09511e-27, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 4.11961e-26]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 5.21855e-27, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 3.40105e-26]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 6.12574e-27, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 4.11961e-26]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 4.69512e-27, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 2.96222e-26]

% ndigits = 100
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 5.74485e-81, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 4.02796e-80]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.42945e-79, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 9.44583e-79]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 5.76549e-81, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 3.97130e-80]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.42444e-77, [xi(40) = 33835.88114392, notabug(40) = 33835.88114392, rel.acc = 1.04581e-77]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.42410e-77, [xi(40) = 33835.88114392, notabug(40) = 33835.88114392, rel.acc = 1.14507e-77]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.57315e-79, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 1.02690e-78]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.42394e-77, [xi(40) = 33835.88114392, notabug(40) = 33835.88114392, rel.acc = 1.14487e-77]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.99963e-77, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 1.18875e-76]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.99972e-77, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 1.18875e-76]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 1.47611e-77, [xi(40) = 33835.88114392, notabug(40) = 33835.88114392, rel.acc = 1.15290e-77]

% ndigits = 500
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 2.88681e-477, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 1.15224e-476]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 3.48807e-477, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 2.16831e-476]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 2.88681e-477, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 1.15224e-476]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 0.00000e+00, [xi(1) = 0.00000000, notabug(1) = 0.00000000, rel.acc = 0.00000e+00]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 0.00000e+00, [xi(1) = 0.00000000, notabug(1) = 0.00000000, rel.acc = 0.00000e+00]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 2.88681e-477, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 1.15224e-476]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 2.88681e-477, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 1.15224e-476]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 0.00000e+00, [xi(1) = 0.00000000, notabug(1) = 0.00000000, rel.acc = 0.00000e+00]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 0.00000e+00, [xi(1) = 0.00000000, notabug(1) = 0.00000000, rel.acc = 0.00000e+00]
rel.acc(sys) = 1.55559e-15, rel.acc(x) = 2.88681e-477, [xi(41) = 5382.66065164, notabug(41) = 5382.66065164, rel.acc = 1.15224e-476]

First column is the system accuracy w.r.t current solution vector, second - current solution vector difference from the first one. Lastly in square brackets we have the current solution vector element with max difference from the first vector and its relative difference.

The relative accuracy of all the solutions are the same (1.55559e-15), this is clearly the maximum achievable accuracy for the system. 

The solution vectors differ from the first one by last ~21 decimal digits. This coincides with the matrix conditioning (as we expected). The solution vector elements with max. difference also differs by the same (very small) margin from the first vector.

Overall, this perfectly matches the expected behavior of matrix solution.

***

Could you please send me references to MATLAB's requirements for solutions of overdetermined system to be bit-by-bit deterministic? In this particular case the MATLAB doesn't use any advanced QR algorithms/multithreading and probably this leads to exactly the same solution each time. But, mathematically it is not guaranteed or expected/justified. 

In toolbox we use advanced algorithms for the case (multi-threadfinfg + D&C + bulge chasing with different execution paths each time). Nonetheless the solutions provide the same relative accuracy and differ from each other by expected margin (which is log10(cond(A)) = ~21 digits). This is how it is expected to behave. Let me know your comments.

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/  

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.

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

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

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

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.

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.