为什么在预转置矩阵上执行矩阵乘法比在非转置矩阵上执行矩阵乘法更快？
zaz*_*azz 10 python transpose numpy matrix
考虑 Python 中的以下代码，其中与乘以非转置矩阵相比，乘以预转置矩阵会产生更快的执行时间：
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)
令人惊讶的是，预转置矩阵的乘法速度更快。人们可能会认为乘法的顺序不会显着影响性能，但似乎存在差异。
为什么与非转置矩阵相比，处理预转置矩阵会导致更快的执行时间？是否有任何根本原因或优化可以解释这种行为？
更新
我已经考虑了有关的评论cache，并且在每个循环上生成新的矩阵：
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()
结果：
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
chr*_*slg 12
在我的机器上这似乎并不明显。
\n运行 1000 次。我得到这些计时（x=非转置，y=转置）。\n（在 y=x 轴下）红点多于蓝点。685/315 更准确。因此，从 p 值角度来看，毫无疑问，这不仅仅是随机效应。\n（抽出 1000 个硬币，其中 685 个头像是一个明显的异常现象）
\n
但从时间上来看，这一点并不明显。该簇主要以 y=x 轴为中心。
\n现在我开始回答这个问题，因为我很确定这是一个缓存问题。当我在工程师学校时（很久以前，当时这些考虑因素现在变得更加重要，并且由老师教授，他们自己可以追溯到更重要的时代），在 HPC 课程中，我们由于缓存效应，从 Fortran 切换到 C 时要非常小心：迭代数组时，按照其在内存中的顺序对其进行交互非常重要（在 numpy 中仍然称为“C”顺序与“ fortran”命令，证明对于比我更关心的人来说，它仍然是一个重要的考虑因素 - 我很少需要在日常工作中关心，因此我调用学校记忆而不是工作记忆的原因）。
\n因为当处理与内存中之前处理过的数字紧邻的数字时，该数字可能已经加载到高速缓存中。而如果您处理的下一个数字是 1 行（按 C 顺序，因此在内存中更远），那么它很可能不在缓存中。对于当今的缓存大小，它需要大矩阵，因此它会产生影响。
\n由于transpose不移动任何数据，而只是调整步幅，因此处理转置矩阵的效果是更改已处理数据在内存中的顺序。所以，如果你考虑朴素算法
for i in range(N):\n for j in range(N):\n res[i,j]=0\n for k in range(N):\n res[i,j] += A[i,k] * B[k,j]\n如果A和B是 C 顺序的，那么矩阵 A 的迭代是按照内存顺序完成的（我们沿着一行，一列又一列地迭代，因此内存中的相邻数字一个接一个），而B不是。
如果顺序颠倒了，例如，因为它们被调换了位置，那么情况就是相反的。B 按照不会造成缓存问题的顺序进行迭代，而 A 则不会。
\n好吧，不需要在这个问题上停留太久，因为我告诉所有这些都是为了解释为什么我想调查缓存问题的可能性（我的目的是将相同的乘法与转置矩阵的副本进行比较，因此它是相同的矩阵乘法，仅改变顺序。并且还要尝试查看矩阵大小是否存在阈值，在该阈值下该现象不可见，这也将验证缓存问题，因为就这一点而言，整个矩阵不得适合缓存）
\n但是，这样做的第一步也是要开始避免偏差，因为第一次计算使用尚未在缓存中的数据，而第二次计算使用已经在缓存中的数据（特别是在整个矩阵适合缓存的情况下）。
\n所以，这是我尝试的第一件事：只是颠倒计算顺序。首先在 transpose_matrix 上计算，然后在矩阵上计算。
\n
这次，移位有利于蓝点（当然，我只改变了计算顺序，没有改变轴的含义。所以x仍然是matrix@matrix计时，y仍然是transposed_matrix
这次的红点数量是 318 对 682。所以，几乎和之前完全相反。
\n所以，结论（至少对我的机器有效）：这确实是一个缓存问题。但是缓存问题仅由以下事实引起transposed_matrix：当您使用它进行计算时，它已经在缓存中（因为数据与矩阵数据相同）。
编辑：关于问题更新。
\n正如我在评论中所说（但由于问题已更新，并且这个答案已经得到了相当多的支持，我认为对于未来的读者来说，它也出现在该答案中很重要），更新是不同的。
\n你的第一个问题是关于A@Avs 的A.T@A.T。第二个似乎更快。但只进行了 1 次操作。因此，正如我所表明的，原因只是因为当第二个操作完成时，A它已经在缓存内存中（当第一个操作完成时，情况并非如此）。因为A.T与 是相同的数据A（不是副本。而是相同的数据，在相同的内存地址）。\n我之前的答案表明，如果您反转并计算A.T@A.T拳头，则A@A，那么它是，相反，A.T@A.T速度较慢，并且以完全相同的比例。
另一种显示方式是
\nimport numpy as np\nimport timeit \n\nA=np.random.normal(0,1,(1000,1000))\nB=A.copy()\n\nA@A\nprint(timeit.timeit(lambda: A@A, number=20))\nA.T@A.T\nprint(timeit.timeit(lambda: A.T@A.T, number=20))\nB@B\nprint(timeit.timeit(lambda: B@B, number=20))\n（事实上，A@A在 timeit 之前执行只是为了确保 20 次计算中的第一个计算不会因缓存考虑而变慢）
在我的计算机上，所有这些操作几乎完全需要 1 秒（number=20选择的因此需要 1 秒）
这次没有缓存效果，因为我们每次运行 21 次，不计算第一次运行的时间。并且没有影响.T
现在，对于您的问题更新，这是另一回事
\nfor i in range(N):\n for j in range(N):\n res[i,j]=0\n for k in range(N):\n res[i,j] += A[i,k] * B[k,j]\n这次，前两次操作只需要 650 毫秒。并不是因为缓存：无论这些操作的顺序如何，它都是相同的。
\n这是因为 numpy 能够检测到A.T和A是相同的矩阵，但需要一次转置操作。\n（它很容易检测到：它是相同的数据地址，但步幅和形状（这里的形状是无论如何，都是正方形；但更重要的是，步幅是反转的）是反转的：A.strides\xe2\x86\x92 (8000,8)、A.T.strides\xe2\x86\x92 (8, 8000)。
所以，numpy很容易意识到这是一种A@A.T情况。因此要应用一种计算速度更快的算法。正如评论中所说（并且之前我在评论中说过，而且几天前其他人也说过，他们误读了你的第一个问题……但这样做是正确的，因为他们提前回答了现在的更新） :A@A.T是对称的。因此，这里有一些简单的修正。
注意
\ntimeit.timeit(lambda: A@B, number=20)\ntimeit.timeit(lambda: A@B.T, number=20)\n均为 1 秒（如A@A、 和A.T@A.T）。所以，很容易理解A@B，A@A， ，A.T@A.T都只是使用一种标准的“矩阵乘法”算法。其中A@A.T，A.T@A使用更快的。
由于B是 的副本A，A@B.T因此具有与 相同的对称结果A@A.T。但这一次，因为它是一个副本，numpy无法意识到它是一种A@A.T情况，无法意识到它是一个对称结果（在计算结果之前）。因此A@B.T具有与 相同的标准“1 秒”计时A@A。虽然A@A.T还没有。
这证实了它确实依赖于这些same address, inverted strides标准。只要不是同一个地址，或者反转的步幅不一样，标准算法，1秒。如果两者是相同的地址，但步幅相反，则采用特殊算法，650 ms。
| 归档时间： |
|
| 查看次数： |
539 次 |
| 最近记录： |