Different results when multiplying the same matrix by its transpose

I multiply the same matrix, "a", by its transpose two different times and compare the outputs:

>> a_ = a.'; %store the transpose in an intermediate variable
>> isequal(a_*a, a.'*a) %returns a logical '0'

Why is the same operation returning different outputs for the same input? Round-off error should not be the cause of this, because the same function should introduce the same degree of round-off error for the same input.

Answer 1

[MathWorks Support Team]

In short, the discrepancy is due to the fact that MATLAB is using standard matrix multiplication for the first product and an optimized algorithm for the second product, introducing different degrees of round-off error. This is expected behavior. Use an epsilon threshold instead of exact comparison to check for equality.

Explanation:

First, you compute the transpose of "a" and assign it to "a_". There are no round-off errors associated with this operation. Then, you multiply "a_" and "a" to get the first product, introducing some expected round-off error. At that point, MATLAB did not know that "a_" stored a transposed version of "a". It just treated "a_" like a unique matrix. 

For the second product, you multiply the transpose of "a" with "a" itself all at once, rather than storing the transpose of "a" in an intermediate variable. Because of this, MATLAB knows that you are multiplying "a" and its transpose in this case. Rather than computing standard matrix multiplication like the first product, the second product is computed using an optimized algorithm for multiplying a matrix by its transpose. 

So, the products were computed with different algorithms even though their inputs were equivalent from a theoretical standpoint. The different algorithms introduce different sources of round-off errors under the hood, which is the source of the discrepancy. 

Workaround:

As with all round-off related issues, you should use an epsilon threshold to check for equality. In your case, it would look like this:

>> M1 = a_*a;
>> M2 = a.'*a;
>> all(abs(M1(:)-M2(:)) < eps(min(abs(M1(:)),abs(M2(:))))) %returns a logical '1'

Line 3 compares the difference between corresponding elements in M1 and M2 to an appropriate epsilon value. The "eps" function should return epsilon values that are larger than any discrepancies caused by round-off errors, yet small enough to detect any substantial inequality. Please read more about "eps" here:

https://www.mathworks.com/help/matlab/ref/double.eps.html

The entire comparison is wrapped in an 'all' command, which will only return true if all of the element-wise comparisons are true.