Loss of accuracy when computing a product between 2 tensors?

Question

I want to iteratively compute the product of a tensor $B$ by contracting its indexes in the following way, where $a,b,c$ are the dimensions of the tensor along those axis, so that the product can be carried away.

The tensor $B$ would have this shape and its axis would be labelled with this convention, being called by $B[1,2,3,4,5,6]$ and thus with shape $(a,b,c,c,b,a)$ according to the previous picture.

One possible way would be to reshape the tensor into indexes $\alpha$ and $\beta$ outlined in the figure, by doing B.reshape(a*b*c,a*b*c) and now the product would be matrix-like. After that we could reshape back into the original shape.

I tried both the possibility in which I link the axis as outlined in the 1st photo with np.tensordot and I compared it to the 3rd picture using np.matmul and np.linalg.matrix_power. Both ways yield similar tensors in the sense that their norm is near, but not equal.

I also tried computing this operation with other numpy functions and compared the results.

• Is this a loss of accuracy in the sense that I should be doing this product in a specific way?

This is an example of the code:

import numpy as np

seed = 2023
a,b,c = 2,3,4

# Generate the tensor B:
shape = (a,b,c,c,b,a)
rn_generator = np.random.default_rng(seed=seed)
# Generate a random complex B tensor
# See this question for details (https://stackoverflow.com/questions/71956506/how-to-generate-uniform-random-complex-numbers-in-python
):
B = np.sqrt(rn_generator.uniform(0, 1, shape)) * np.exp(1.j * rn_generator.uniform(0, 2 * np.pi, shape))



# Reshaping the indexes to that we can perform the product in a 'matrix-like' form
Bp = B.reshape(a*b*c,a*b*c)
Bv0 = np.matmul(Bp,Bp).reshape(shape)
Bv1 = np.linalg.matrix_power(Bp,2).reshape(shape)

# The previous should be equivalent to summing over the correct indexes:
Bv2 = np.tensordot(B,B,axes=([3,4,5],[2,1,0]))



#However we see that they are different
np.allclose(Bv0,Bv2)
np.allclose(Bv1,Bv2)

# We see how the norm (invariant) is different
np.linalg.norm(Bv0)
np.linalg.norm(Bv1)
np.linalg.norm(Bv2)




# We compare with other methods:
Bv3 = np.einsum('ijklmn,nmlpqr->ijkpqr',B,B)

from ncon import ncon
# Here negative indexes are free and repeated are being summed
Bv4 = ncon([B,B],[[-1,-2,-3,1,2,3],[3,2,1,-4,-5,-6]])

# We check their norms
np.linalg.norm(Bv3)
np.linalg.norm(Bv4)

# We also see that they are equal to v2
np.allclose(Bv2,Bv3)
np.allclose(Bv2,Bv4)

The code yields:

False

False

55.10565356464396

55.10565356464396

58.73821610499646

58.73821610499646

58.73821610499646

True

True

We can see that indeed the 2 operations that should be equivalent yield different results.