2

我正在编写一个考虑速度的数值算法。我在 scipy/numpy (scipy.linalg.expm2, scipy.linalg.expm) 中遇到了两个矩阵指数函数。但是，我有一个事先知道是对角线的矩阵。这些 scipy 函数在运行之前是否检查矩阵是否对角线？显然，对于对角矩阵，求幂算法可以更快，我只是想确保它们在做一些聪明的事情——如果不是，有没有简单的方法来做到这一点？

4

3 回答 3

4

如果一个矩阵是对角矩阵，那么它的指数可以通过对主对角线上的每个条目取幂来获得，所以你可以通过以下方式计算它：

np.diag(np.exp(np.diag(a)))
于 2013-04-18T12:16:51.167 回答
3

如果您知道 A 是对角线并且您想要 k 次方：

def dpow(a, k):
    return np.diag(np.diag(a) ** k)

检查矩阵是否对角线：

def isdiag(a):
    return np.all(a == np.diag(np.diag(a)))

所以 ：

def pow(a, k):
    if isdiag(a):
        return dpow(a, k)
    else:
        return np.asmatrix(a) ** k

类似地，对于指数（您可以从一组 pow 的扩展数学上得到），您可以执行以下操作：

def dexp(a, k):
    return np.diag(np.exp(np.diag(a)))

def exp(a, k):
    if isdiag(a):
        return dexp(a, k)
    else:
        #use scipy.linalg.expm2 or whatever
于 2013-04-18T12:20:58.683 回答
1

我开发了一个工具，可以帮助更快地做与 HYRY 相同的事情，但通过就地进行：

def diagonal(array):
    """ Return a **view** of the diagonal elements of 'array' """
    from numpy.lib.stride_tricks import as_strided
    return as_strided(array,shape=(min(array.shape),),strides=(sum(array.strides),))

# generate a random diagonal array
d  = np.diag(np.random.random(4000))

# in-place exponent of the diagonal elements
ddiag = diagonal(d)
ddiag[:] = np.exp(ddiag)

# timeit comparison with HYRY's method
%timeit -n10 np.diag(np.exp(np.diag(d)))   
    # out> 10 loops, best of 3: 52.1 ms per loop
%timeit -n10 ddiag = diagonal(d); ddiag[:] = np.exp(ddiag)
    # out> 10 loops, best of 3: 108 µs per loop

现在，

  • HYRY 的方法是对角线长度的二次方（可能是因为新的数组内存分配），因此如果您的矩阵维度很小，则差异可能不会那么大。

  • 你需要对就地计算没问题

  • 最后，非对角元素是 0，所以它们的指数应该是 1，不是吗？在我们的两种方法中，非对角线都是 0。

对于最后一部分，如果您希望所有非对角线元素为 1，那么您可以执行以下操作：

d2 = np.ones_like(d); 
diagonal(d2)[:] = np.exp(np.diag(d))

print (d2==np.exp(d)).all()  # True

但这对数组大小是线性的，所以对角线长度是二次的。timeit 为 4000x4000 阵列提供约 90 毫秒，为 2000x2000 提供 22.3 毫秒。

最后，您也可以就地执行以加快速度：

diag = np.diag(d)
d[:]=1
diagonal(d)[:] = np.exp(diag)

Timeit 为 4000^2 数组提供 66.1ms，为 2000^2 提供 16.8ms

于 2013-04-19T12:18:44.673 回答