矩阵对角化 python_python - 对角化+特征分解比朴素矩阵乘法花费更多的时间 - SO中文参考 - www.soinside.com...

对角化+特征分解比朴素矩阵乘法花费更多的时间

问题描述 投票：2回答：2

假设我们正在运行一些粒子模拟。我们有一个点p，例如(1, 1, 2)我们想应用一次线性变换N次。如果转换用矩阵A表示，那么最终的变换将由A^N . p给出。矩阵乘法的成本很高，我假设特征分解和对角线化将加快整个过程。但是令我惊讶的是，这种据说改进的方法花费了更多时间。我在这里错了吗？import timeit

mysetup = '''

import numpy as np

from numpy import linalg as LA

from numpy.linalg import matrix_power

EXP = 5 # no. of time linear transformation is applied

LT = 10 # range from which numbers are picked at random for matrices and points.

N = 100 # dimension of the vector space

A_init = np.random.randint(LT, size=(N, N))

A = (A_init + A_init.T)/2

p = np.random.randint(LT, size=N)

def run_sim_1():

An = matrix_power(A, EXP)

return An @ p

def run_sim_2():

λ, V = LA.eig(A)

Λ = np.diag(λ)

Λ[np.diag_indices(N)] = λ ** EXP

An = V @ Λ @ V.T

return An @ p

'''

# code snippet whose execution time is to be measured

# naive implementation

mycode_1 = '''run_sim_1()'''

print(timeit.timeit(setup = mysetup, stmt = mycode_1, number = 1000))

# time taken = 0.14894760597962886

# improved code snippet whose execution time is to be measured

# expecting this to take much less time.

mycode_2 = '''run_sim_2()'''

# timeit statement

print(timeit.timeit(setup = mysetup, stmt = mycode_2, number = 1000))

# time taken = 8.035318267997354

python

numpy

linear-algebra

2个回答

2

投票

这很难权威地回答。 matrix multiplication and eigendecomposition的标准实现都是O(n ^ 3)，因此没有先验的理由期望一个比另一个要快。有趣的是，我的经验是本征分解通常比单矩阵乘法慢得多，所以这个结果并不完全让我感到惊讶。

因为在这种情况下，矩阵幂运算涉及二十乘法，所以我可以理解为什么您会期望它比本征分解慢。但是，如果您查看source code，则会显示此有趣的花絮：# Use binary decomposition to reduce the number of matrix multiplications.

# Here, we iterate over the bits of n, from LSB to MSB, raise `a` to

# increasing powers of 2, and multiply into the result as needed.

z = result = None

while n > 0:

z = a if z is None else fmatmul(z, z)

n, bit = divmod(n, 2)

if bit:

result = z if result is None else fmatmul(result, z)

因此，实际上并没有进行20次乘法！它使用分而治之的方法来减少数量。看起来在最坏的情况下，对于给定的幂p，这仍然会进行p乘法。当p全为1时，即比2的幂小1时，就会发生这种情况。但是在平均情况下，它将大约进行p / 2 + log(p / 2)乘法(因为平均而言，只有一半的位是1)。并且在最佳情况下，它将仅执行log(p)乘法(当p为2的幂时)。

我不确定这些数字是否100％；可能值得测试一下p从15到16时性能是否突然急剧提高。如果是这种情况，那么您会发现性能的差异[取决于怪异这些细节。我应该加上这个：将向量直接相乘会不会比将矩阵提升为幂更快？二十个矢量乘法仍然是O(n^2)，不是吗？但是也许您真正想做的是对10k向量执行此操作，在这种情况下，矩阵幂方法显然更好。

1

投票

my_code_1和my_code_2都只包含一个def语句。因此，您对timeit的调用仅是定时定义函数所需的时间。函数从不

被调用。将函数

definitions移至设置代码，并仅通过调用适当的函数来替换要计时的语句，例如mycode_1 = '''

run_sim_1()

'''然后您应该(大幅度)降低传递给number的timeit的值。然后，您必须修复run_sim_2()才能执行正确的计算：

def run_sim_2():

λ, V = LA.eig(A)

Λ = np.diag(λ)

Λ[np.diag_indices(N)] = λ ** 20

An = V @ Λ @ V.T

return An @ p一旦进行了这些更改，您仍然会发现run_sim_1()更快。有关可能的原因，请参见@senderle的答案。

热门问题