Python`exm`是一个`(N,M,M)矩阵

内容导读

互联网集市收集整理的这篇技术教程文章主要介绍了Python`exm`是一个`(N,M,M)矩阵，小编现在分享给大家，供广大互联网技能从业者学习和参考。文章包含2166字，纯文字阅读大概需要4分钟。

内容图文

设A是(N,M,M)矩阵(N非常大),我想为范围(N)中的每个n计算scipy.linalg.expm(A [n,：,]].我当然可以使用for循环,但我想知道是否有一些技巧以更好的方式执行此操作(类似于np.einsum).

我对其他操作有相同的问题,例如反转矩阵(反馈在注释中解决).

解决方法:

根据矩阵的大小和结构,您可以做得比循环更好.

假设您的矩阵可以对角化为A = V D V ^( – 1)(其中D在其对角线上具有特征值,V包含相应的特征向量作为列),您可以将矩阵指数计算为

exp(A) = V exp(D) V^(-1)

其中exp(D)只包含对角线中每个特征值lambda的exp(lambda).如果我们使用指数函数的幂级数定义,这很容易证明.如果矩阵A是正常的,则矩阵V是单一的,因此可以通过简单地取其伴随来计算其逆.

好消息是numpy.linalg.eig和numpy.linalg.inv都可以正常使用堆叠矩阵：

import numpy as np
import scipy.linalg

A = np.random.rand(1000,10,10)

def loopy_expm(A):
    expmA = np.zeros_like(A)
    for n in range(A.shape[0]):
        expmA[n,...] = scipy.linalg.expm(A[n,...])
    return expmA

def eigy_expm(A):
    vals,vects = np.linalg.eig(A)
    return np.einsum('...ik, ...k, ...kj -> ...ij',
                     vects,np.exp(vals),np.linalg.inv(vects))

请注意,在指定einsum调用中的操作顺序时,可能还有一些优化空间,但我没有对此进行调查.

测试上面的随机数组：

In [59]: np.allclose(loopy_expm(A),eigy_expm(A))
Out[59]: True

In [60]: %timeit loopy_expm(A)
824 ms ± 55.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [61]: %timeit eigy_expm(A)
138 ms ± 992 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

那已经很好了.如果你足够幸运,你的矩阵都是正常的(例如,因为它们是真正对称的)：

A = np.random.rand(1000,10,10)
A = (A + A.transpose(0,2,1))/2

def eigy_expm_normal(A):
    vals,vects = np.linalg.eig(A)
    return np.einsum('...ik, ...k, ...jk -> ...ij',
                     vects,np.exp(vals),vects.conj())

注意输入矩阵的对称定义和einsum模式内的转置.结果：

In [80]: np.allclose(loopy_expm(A),eigy_expm_normal(A))
Out[80]: True

In [79]: %timeit loopy_expm(A)
878 ms ± 89.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [80]: %timeit eigy_expm_normal(A)
55.8 ms ± 868 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

对于上面的示例形状,这是15倍的加速.

应该注意的是,scipy.linalg.eigm根据文档使用Padé近似.这可能意味着如果你的矩阵是病态的,那么特征值分解可能会产生与scipy.linalg.eigm不同的结果.我不熟悉这个功能如何工作,但我希望它对病理输入更安全.

内容总结

以上是互联网集市为您收集整理的Python`exm`是一个`(N,M,M)矩阵全部内容，希望文章能够帮你解决Python`exm`是一个`(N,M,M)矩阵所遇到的程序开发问题。 如果觉得互联网集市技术教程内容还不错，欢迎将互联网集市网站推荐给程序员好友。

内容备注

版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 gblab@vip.qq.com 举报，一经查实，本站将立刻删除。

内容手机端