python - 什么时候不做重复线性变换的特征分解
问题描述
假设我们有一个点 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
解决方案
这有点难以权威地回答。矩阵乘法和特征分解的标准实现都是 O(n^3),因此没有先验理由期望一个比另一个快。有趣的是,我的经验是特征分解通常比单个矩阵乘法慢得多,所以这个结果并不完全让我感到惊讶。
因为在这种情况下,矩阵幂运算涉及20次乘法,所以我可以理解为什么您可能期望它比特征分解要慢。但是,如果您查看源代码,就会发现这个有趣的花絮:
# 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 次乘法!它使用分而治之的方法来减少这个数字。在仔细考虑了这个非常优雅的算法之后,我相信它永远不会2*log(p)
对给定的幂做更多的乘法运算p
。当 的所有位p
都为 1 时,即p
1 小于 2 的幂时,达到此最大值。
结果是,尽管在理论上特征分解可能比重复矩阵乘法更快,但它会带来恒定的开销,这使得它在p
变得非常大之前效率较低——可能大于任何实际值。
我应该补充一点:直接将向量相乘不会比将矩阵乘以幂吗?二十个向量乘法仍然是O(n^2)
,不是吗?但也许您真正想做的是对 10k 个向量执行此操作,在这种情况下,矩阵幂方法显然更胜一筹。
推荐阅读
- swift - Swift Cocoa 通过 PHP 向应用发送通知
- ios - 选择搜索栏导致表格视图和搜索栏消失
- scala - Scala 测试错误:值 / 不是 sbt.Configuration 的成员
- sql - SQL 使用新列创建视图以保存来自 DateDiff 的值
- python - 使用 pip 安装包后无法导入子包
- python - 如何打印 Locust 响应 (JSONDecodeError)
- javascript - Bokeh CustomJS 传递字形数组
- ruby-on-rails - 在 Ruby 中使用单个命令更新多条记录
- blockchain - 无法在节点和运行时之间转换参数“tx”:错误解码字段调用
- wix - 我需要从我的 wix 安装程序中删除功能选择屏幕