Precision quandaries: why can I print 64 digits?

I've been testing a Modified Gram-Schmidt algorithm that is part of a larger pipeline:

function [Q, R] = mgs(X)
% [Q, R] = MGS(X) performs Modified Gram-Schmidt on the m x s matrix X.

%%
% Pre-allocate memory for Q and R
[m, s] = size(X);
R = zeros(s,s);
Q = zeros(m,s);

% Modified Gram-Schmidt
R(1,1) = norm(X(:,1));
Q(:,1) = X(:,1)/R(1,1);
for k = 1:s-1
    w = X(:,k+1);
    for j = 1:k
        R(j,k+1) = Q(:,j)'*w;
        w = w - Q(:,j)*R(j,k+1);
    end
    R(k+1,k+1) = norm(w);
    Q(:,k+1) = w/R(k+1,k+1);
end
end

Farther down the software pipeline, I measure loss of orthogonality and backward error. My colleague ran our software on her machine and found that for several key examples I have nearly full orthogonality and she has none. I am 100% certain we are using the same X and software (we checked; everything is pushed and pulled through Git).

We found something surprising with the following script that may explain the issue:

% mgs test
X = magic(10); X = X(:,1:2);
[Q, R] = mgs(X);

fprintf('||I - Q''*Q|| = %0.64e\n', norm(eye(2) - Q'*Q)); % loss of orthogonality
fprintf('||Q*R - X|| = %0.64e\n', norm(Q*R - X)); % backward error

% Results on my machine:
% mgs:
% ||I - Q'*Q|| = 4.1017407521273322775592436112613973088682572487422006712876054735e-16
% ||Q*R - X|| = 1.7763568394002504646778106689453125000000000000000000000000000000e-15

The same script on her machine prints trailing 0's after the 16th digit for both quantities.

I had assumed that Matlab computed in double precision regardless of the machine, but it would seem that Matlab is somehow utilizing higher precision on my machine than on my colleague's. Or is there another explanation for why I can print 64 digits of precision?

My specs: Matlab 2019a, Windows 10 64-bit, Intel Core i7-8850U CPU @ 1.80GHz, 16GB RAM, Thinkpad X1 Carbon Gen 6

Colleague's: Matlab 2017a, Windows 7 (don't know the rest but her Thinkpad is about 6 years older than mine, so likely lacking some features that mine has)

Answer 1

[James Tursa]

What you describe sounds like it is just a display issue. Earlier versions of MATLAB on PC machines used a library print function that only printed significant digits. Later versions of MATLAB on PC machines switched to a library function that prints all of the requested digits in an exact binary floating point to decimal sense ... but still only about 15-16 of those digits are significant. The underlying double precision arithmetic and storage is the same. I don't recall the exact MATLAB version when they made the switch, but I can look it up to confirm.

Results on PCWIN64 R2017a:

||I - Q'*Q|| = 4.1017407521273323000000000000000000000000000000000000000000000000e-16
||Q*R - X|| = 1.7763568394002505000000000000000000000000000000000000000000000000e-15
>> num2strexact(norm(eye(2) - Q'*Q))
ans =
    '4.1017407521273322775592436112613973088682572487422006712876054734806530177593231201171875e-16'
>> num2strexact(norm(Q*R - X))
ans =
    '1.7763568394002504646778106689453125e-15'

Results on PCWIN64 R2017b:

||I - Q'*Q|| = 4.1017407521273322775592436112613973088682572487422006712876054735e-16
||Q*R - X|| = 1.7763568394002504646778106689453125000000000000000000000000000000e-15
>> num2strexact(norm(eye(2) - Q'*Q))
ans =
    '4.1017407521273322775592436112613973088682572487422006712876054734806530177593231201171875e-16'
>> num2strexact(norm(Q*R - X))
ans =
    '1.7763568394002504646778106689453125e-15'

Results on PCWIN64 R2019a:

||I - Q'*Q|| = 4.1017407521273322775592436112613973088682572487422006712876054735e-16
||Q*R - X|| = 1.7763568394002504646778106689453125000000000000000000000000000000e-15
>> num2strexact(norm(eye(2) - Q'*Q))
ans =
    '4.1017407521273322775592436112613973088682572487422006712876054734806530177593231201171875e-16'
>> num2strexact(norm(Q*R - X))
ans =
    '1.7763568394002504646778106689453125e-15'

So you can see that it is just a display issue, and the display change came in R2017b. The underlying calculations actually produce the exact same result in all three versions. (num2strexact is a mex routine that prints the exact binary floating point to decimal conversion)

But, regardless of how many decimal digits actually get printed, the double floating point arithmetic in the background is the same on your machine as it is on hers.

[Erin Carson]

I am the colleague in question in the original post. Running the exact same mgs.m function as posted and the same script as posted, I actually don't get the same answer for the loss of orthogonality, not even the same leading digit. I am running 64-bit Windows 7 and Matlab R2017a.

||I - Q'*Q|| = 1.1796721942137590000000000000000000000000000000000000000000000000e-16
||Q*R - X|| = 1.7763568394002505000000000000000000000000000000000000000000000000e-15

[James Tursa]

You copy & pasted that exact code above and get different results? I wonder if you have some rounding modes set differently on your machine, since your differences still appear to be in the noise at 1e-16.

[Erin Carson]

Yes, copy-and-pasted that exact code. Any idea how I can check the rounding modes?

[Kathryn Lund]

@JamesTursa, thanks for pointing out the display issue. I should also mention, we have a third machine for comparison that's running Matlab 2019b, Update 2, on the login node of a Linux cluster (Ubuntu 16.04.6 LTS GNU/Linux 4.13.0-37-generic x64), and it's even less accurate than my or my Erin Carson's machine:

||I - Q'*Q|| = 6.5163573027050223246570951928891080315475899369798140625675841875e-16
||Q*R - X|| = 7.1054273576010018587112426757812500000000000000000000000000000000e-15

[James Tursa]

Well, 1e-16 is still in the noise and could be accounted for by slightly different order of operations in the background. The only way to remove that would be to remove all vectorization from your code.

As far as rounding mode, I don't know how to do that in MATLAB. Do you have a C compiler?

[Walter Roberson]

Loop code on vectors or arrays can end up automatically vectorized if the array is large enough and MATLAB recognizes the pattern of operations as ones that can be sent to LAPACK or MKL

[Kathryn Lund]

We both have C compilers.

And along similar lines: the only way to standardize vectorization across our machines would be to change the way Matlab compiles C-libraries, right? Is that even possible?