Consider the following code in Python, where multiplying a pre-transposed matrix yields faster execution time compared to multiplying a non-transposed matrix:
import numpy as np
import time
# Generate random matrix
matrix_size = 1000
matrix = np.random.rand(matrix_size, matrix_size)
# Transpose the matrix
transposed_matrix = np.transpose(matrix)
# Multiply non-transposed matrix
start = time.time()
result1 = np.matmul(matrix, matrix)
end = time.time()
execution_time1 = end - start
# Multiply pre-transposed matrix
start = time.time()
result2 = np.matmul(transposed_matrix, transposed_matrix)
end = time.time()
execution_time2 = end - start
print("Execution time (non-transposed):", execution_time1)
print("Execution time (pre-transposed):", execution_time2)
Surprisingly, multiplying the pre-transposed matrix is faster. One might assume that the order of multiplication should not affect the performance significantly, but there seems to be a difference.
Why does processing a pre-transposed matrix result in faster execution time compared to a non-transposed matrix? Is there any underlying reason or optimization that explains this behavior?
UPDATE
I've taken the comments about the cache
into consideration and I'm generating new matrices on each loop:
import numpy as np
import time
import matplotlib.pyplot as plt
# Generate random matrices
matrix_size = 3000
# Variables to store execution times
execution_times1 = []
execution_times2 = []
# Perform matrix multiplication A @ B^T and measure execution time for 50 iterations
num_iterations = 50
for _ in range(num_iterations):
matrix_a = np.random.rand(matrix_size, matrix_size)
start = time.time()
result1 = np.matmul(matrix_a, matrix_a)
end = time.time()
execution_times1.append(end - start)
# Perform matrix multiplication A @ B and measure execution time for 50 iterations
for _ in range(num_iterations):
matrix_b = np.random.rand(matrix_size, matrix_size)
start = time.time()
result2 = np.matmul(matrix_b, matrix_b.T)
end = time.time()
execution_times2.append(end - start)
# Print average execution times
avg_execution_time1 = np.mean(execution_times1)
avg_execution_time2 = np.mean(execution_times2)
#print("Average execution time (A @ B^T):", avg_execution_time1)
#print("Average execution time (A @ B):", avg_execution_time2)
# Plot the execution times
plt.plot(range(num_iterations), execution_times1, label='A @ A')
plt.plot(range(num_iterations), execution_times2, label='B @ B.T')
plt.xlabel('Iteration')
plt.ylabel('Execution Time')
plt.title('Matrix Multiplication Execution Time Comparison')
plt.legend()
plt.show()
# Display BLAS configuration
np.show_config()
Results:
blas_mkl_info:
libraries = ['mkl_rt']
library_dirs = ['C:/Users/User/anaconda3\\Library\\lib']
define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
include_dirs = ['C:/Users/User/anaconda3\\Library\\include']
blas_opt_info:
libraries = ['mkl_rt']
library_dirs = ['C:/Users/User/anaconda3\\Library\\lib']
define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
include_dirs = ['C:/Users/User/anaconda3\\Library\\include']
lapack_mkl_info:
libraries = ['mkl_rt']
library_dirs = ['C:/Users/User/anaconda3\\Library\\lib']
define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
include_dirs = ['C:/Users/User/anaconda3\\Library\\include']
lapack_opt_info:
libraries = ['mkl_rt']
library_dirs = ['C:/Users/User/anaconda3\\Library\\lib']
define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
include_dirs = ['C:/Users/User/anaconda3\\Library\\include']
Supported SIMD extensions in this NumPy install:
baseline = SSE,SSE2,SSE3
found = SSSE3,SSE41,POPCNT,SSE42,AVX,F16C,FMA3,AVX2
not found = AVX512F,AVX512CD,AVX512_SKX,AVX512_CLX,AVX512_CNL
np.show_config()
.b.T @ b.T
) to rule out caching effects related to the fact that you already performeda @ a
, which can leave reusable data in cache for computinga.T @ a.T
.np.matmul(matrix_b, matrix_b.T)
is symmetrical, by exploiting this behaviour it can be almost twice as fast asnp.matmul(matrix_a, matrix_a)
.A@A
vs[email protected]
. Which, as I've shown, is not showing any time different, but for the fact that the fist computation takes more time that the 2nd (whether the 1st isA@A
and the second[email protected]
, or the 1st is[email protected]
and the 2ndA@A
, it doesn't matter: the 1st is slower). Because of cache.