何时不对重复线性变换进行特征分解

Question

假设我们有一个点 p，例如(1, 2, 3) 我们要对其应用 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

Answer 1

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

因为本例中的矩阵幂运算涉及 20 次乘法，我明白为什么您可能认为它比特征分解慢。但是，如果您查看 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 进行 2*log(p) 次以上的乘法运算。当 p 的所有位都为 1 时，即当 p 小于 2 的幂时达到此最大值。

结果是，虽然特征分解在理论上可能比重复矩阵乘法更快，但它带来了恒定的开销，导致效率降低，直到 p 变得非常大——可能比任何实际值都大。

我应该补充一点：直接乘以向量不会比提高矩阵的幂更快吗？二十次向量乘法仍然是 O(n^2)，不是吗？但也许您真正想要做的是对 10k 个向量执行此操作，在这种情况下，矩阵幂方法显然更优越。

Answer 2

my_code_1 和 my_code_2 都只包含一个 def 语句。因此，您对 timeit 的调用只是计算定义函数所需的时间；永远不会 调用函数 .

将函数 definitions 移动到设置代码中，并将要计时的语句替换为仅调用适当的函数，例如

mycode_1 = ''' 
run_sim_1()
'''

那么你应该降低（很多）传递给 timeit 的 number 的值。然后你必须修复 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 的回答。

何时不对重复线性变换进行特征分解

When not to do eigendecomposition for repeated linear transformation

python

numpy

linear-algebra