A single precision matrix multiplication problem
12 vues (au cours des 30 derniers jours)
plot(c) %=> It's not all 0
The algorithm says y and z should be the same,But c is not all 0,why is this?
Even multiplying individual rows and columns makes a difference.
John D'Errico le 18 Mai 2023
Modifié(e) : John D'Errico le 18 Mai 2023
Why are you confused? NEVER trust the least significant bits of a floating point number.
The numbers in a are singles. So the least significant bits in a are of magnitude
Now you add and subtract them, multiplying by numbers that are on the order of 1e7. Do you appreciate that when MATLAB does the computations you wrote, there is no presumption that it will add and subtract them in EXACTLY the same sequence as the explicit for loops you wrote? And that means again, NEVER TRUST THE LEAST SIGNIFICANT BITS OF A FLOATING POINT RESULT.
The difference is now on the order of 1, actually, 0.25. WHAT A SURPRISE! NOT.
Floating point arithmetic is NOT exact mathematically correct arithmetic. In fact, it is just an approximation to exact mathematically correct arithmetic. A very good approximation in most cases. (Slightly less good for singles than doubles, but the difference is irrelevant in this respect.) As an approximation, there are always tiny errors made. And then when you magnify those errors, you can find them, IF you look very closely.
So, just for kicks, suppose I asked you to compute the result
2/3 - 1/3 - 1/3
BUT, you must use 4 digit decimal arithmetic, first approximating each number as a 4 digit decimal number. You would now get:
0.6667 - 0.3333 - 0.3333 = 0.0001
So even though the answer is mathematically zero, when we use a finite number of digits in any form, the result need not be zero.
Does order make a difference?
x = single([0.1 0.2 0.3]);
x(3) - x(1) - x(2)
-(x(1) + x(2)) + x(3)
Again, NEVER TRUST THE LEAST SIGNIFICANT BITS OF A RESULT. In both cases, the answer is effectively zero, to within a tolerance of eps('single').
The answer is to learn to use tolerances. Learn to not trust those least significant bits. Never test for exact equality of floating point numbers, at least not until you know enough about them that you fully understand when you can.
Plus de réponses (1)
James Tursa le 19 Mai 2023
Modifié(e) : James Tursa le 22 Mai 2023
Somewhat related side comments:
In general, the BLAS/LAPACK libraries do not contain mixed type argument routines. E.g., there are no BLAS matrix multiply routines to multiply a single precision matrix with a double precision matrix. To accomplish this with BLAS routines one must first make a temporary copy of the single matrix as a double matrix and then call the appropriate BLAS routine. This of course will chew up extra time and memory. The only exceptions to this are some accumulation routines (e.g., dot product of single precision inputs but using a double precision accumulator). That being said, I don't know if MATLAB ever makes use of these extra precision accumulator routines that are available.
Related to that, there are also no mixed real/complex argument routines in the BLAS/LAPACK libraries. Either all arguments have to be real or all arguments have to be complex. And all of the complex variables are assumed to be interleaved memory model (this was true even back in the R2017b and earlier days).
MATLAB R2017b and prior used a separate complex memory model, meaning the real part of a variable was held in contiguous memory that was separate from the imaginary part which was also held in contiguous memory. For MATLAB R2018a and later this changed to an interleaved complex memory model, where the real part is interleaved with the imaginary part in memory. This now aligns better with the BLAS/LAPACK libraries and how complex variables are typically held in memory in other languages such as Fortran and C/C++, but has implications for how the BLAS/LAPACK routines can be used in the background. A couple of examples:
(Real Matrix) * (Complex Matrix)
R2017b and earlier: This can be accomplished with a series of Real BLAS routine calls on the individual real and imaginary parts. No need for temporary copies.
R2018a and later: This requires that a temporary copy of the real matrix be made as an interleaved complex matrix with imaginary part 0 so that the Complex BLAS routine can be used. Takes extra time and memory compared to R2017b and earlier, and also lots of unnecessary 0 multiplies will be done (which might give some NaN results that don't appear in R2017b if there are inf's involved).
R2017b and earlier: This involves calling complex LAPACK routines, so the separate real and imaginary parts of the variable must first be copied into an equivalent interleaved variable. Then the appropriate LAPACK routines are called. Then the interleaved result must be copied into a separate memory model variable. Takes extra time and memory compared to R2018b and later.
R2018b and later: This can be accomplished with direct complex LAPACK routine calls. No need for temporary copies.
Symmetrix Matrix Results:
A = rand(5000);
AT = A'; % Explicit transpose constructed in memory
tic; x = A' * A; toc % Symmetrix matrix multiply routine called, A' not explicitly formed
tic; y = AT * A; toc % Generic matrix multiply routine called using explicitly formed A'
isequal(x,x') % Strictly symmetric result
isequal(y,y') % Might not be strictly symmetric result
So, with this last example, you can influence which BLAS routines are called in the background by how you write your code. Care must be taken to maintain strict symmetry of results. In the A' * A example, MATLAB can recognize the symmetry of the inputs and call a symmetric routine in the background without explicitly forming the transpose A'. Takes less time and guaranteed results in strict symmetric result. In the AT * A example, the transpose A' is explicitly formed but MATLAB has no clue that AT is the exact transpose of A, so a generic matrix multiply routine is called which takes longer and doesn't guarantee a strict symmetric result. In that last case we happened to get strict symmetry, but this will highly depend on your MATLAB/BLAS version and can't be guaranteed.