Transpose property of matrix multiplication does not equal exactly when last dimension of the array is 1

user_10381614 0
Viewed 95

I have a weird problem that I ran into the other day and I wondered if someone knows the reason for this weirdness. I'm sorry if it's a duplicate of some other posts, but I couldn't find similar posts on this.

I was testing the transpose property on matrix multiplications where (A @ B)T == BT @ AT holds true.

I wrote a simple code to test it out and to my surprise when the last dimension for the array is 1, the array is "close" but not "equal" in some cases (ignore the ugly code...).

def transpose_property(last_dim: int = 1):
    A = np.random.rand(4, 3)
    B = np.random.rand(3, last_dim)

    AT = A.transpose((1, 0)).copy()
    BT = B.transpose((1, 0)).copy()

    M_left = (A @ B).transpose(1, 0).copy()
    M_right = BT @ AT

    equal = np.array_equal(M_left, M_right)
    close = np.allclose(M_left, M_right)
    return [equal, close]

I ran some tests where I change the last dimension (last_dim) of B from 1 to 3 and counted how many times (A@B)T and BT @ AT were "equal" (np.array_equal) and "close" (np.allclose):

np.random.seed(0)

num_tests = 10000

for dim in range(1, 4):
    results = [transpose_property(last_dim=dim) for _ in range(num_tests)]
    equals = [r[0] for r in results if r[0]]
    closes = [r[1] for r in results if r[1]]
    print(f"dim={dim}:")
    print(f"\tequals: {len(equals)}/{num_tests}")
    print(f"\tcloses: {len(closes)}/{num_tests}")

Results of this script are shown below:

dim=1:
        equals: 3452/10000
        closes: 10000/10000
dim=2:
        equals: 10000/10000
        closes: 10000/10000
dim=3:
        equals: 10000/10000
        closes: 10000/10000

I'm puzzled that when the last dimension is 1, the number of equals is low, but when the last dimension is greater than 1, all of them match.

I think it might be due to floating-point precision rounding, but I don't understand why it's only when the last dimension is 1. How could this be?

I would like to note that:

  • this happens when you add more dimensions (e.g., A.shape: (3, 3, 3) and B.shape: (3, 3, 1) and transpose the last two dimensions).
  • this will not happen when the array is an integer (e.g., A = np.arange(12).reshape(6, 2) and B = np.arange(2).reshape(2, 1)).
  • This problem is not detrimental since the end arrays are "close". It's just that this has been on my mind and wanted to figure out why.
transpose matrix numpy arrays python
1 Answers

Here is a simple test to show you what is going on:

from itertools import count
import numpy as np

last_dim = 1

for i in count(1):
    A = np.random.rand(4, 3)
    B = np.random.rand(3, last_dim)

    AT = A.transpose((1, 0)).copy()
    BT = B.transpose((1, 0)).copy()

    M_left = (A @ B).transpose(1, 0).copy()
    M_right = BT @ AT

    equal = np.array_equal(M_left, M_right)
    if not equal:
        break

print(i)
print(M_left - M_right)

On my machine, I get a discrepancy on the first iteration:

array([[ 1.11022302e-16,  0.00000000e+00,  0.00000000e+00, -1.11022302e-16]])

Floats are finite-precision integers with a scale factor. The order of multiplication and addition can cause an error in the last bit to creep in, which is what you are seeing here:

>>> 2**-53
1.1102230246251565e-16

Remember that adding four numbers of order 1 changes the scale of the number to order 4, so you may end up losing up to two bits, depending on how the rounding works out. You will lose more bits as the size of the matrix increases.

user_29887300
Related