scipy.linalg.expm 与手动编码之间的区别

Question

我试图在scipy.linalg.expm 中实现矩阵指数函数。我从kaityo256's github repository 获得灵感。因此，我写下了以下内容。

from scipy.linalg import expm
from scipy.linalg import eigh
from scipy.linalg import inv
from math import exp as math_exp
from numpy import array, zeros
from numpy.random import random_sample
from numpy.testing import assert_allclose

def diag2sqr(x):
    '''Makes an square matrix from a diagonal one.

    Takes a 1d matrix. Determines its data type.
    Finds out the shape of the 1d matrix.
    Makes an empty  square matrix with  both
    dimensions equal to largest (nonzero) dimension of
    the 1d matrix. It then fills the elements of the
    1d matrix into diagonal slots of the empty
    square one.


    Parameters
    ----------
    x : ndarray
        ndarray of be coverted to a square ndarray

    Returns
    -------
    xsqr : ndarray
        ndarray with diagonals sameas those of x
        all other elements are zero
        dtype same as that of x

    '''

    x_flat = x.ravel()
    xsqr = zeros((x_flat.shape[0], x_flat.shape[0]), dtype=x.dtype)
    # Making the empty matrix
    for i in range(x_flat.shape[0]):
        xsqr[i, i] = x_flat[i]
        # filling up the ith element

    print('xsqr', xsqr)
    return xsqr


def kaityo_expm(x, ):
    '''Exponentiates an ndarray (kaityo).

    Exponentiates a ndarray in the most naive way.

    Parameters
    ----------
    x : ndarray
        The ndarray to be exponentiated

    Returns
    -------
    kexpm : ndarray
        x after exponentiating

    '''

    rx, ux = eigh(x)
    # Find eigenvalues and eigenvectors
    # eigenvectors composed to form a unitary
    ux_inv = inv(ux)
    # Inverse of the unitary
    # tx = diag([array([math_exp(i) for i in rx]).ravel()])
    # tx = array([math_exp(i) for i in rx])
    tx = diag2sqr(array([math_exp(i) for i in rx]))
    # Constructing the diagonal matrix
    kexpm1 = tx@ux_inv
    kexpm = ux@kexpm1

    return kexpm

之后，我尝试测试上面的代码和scipy.linalg.expm。

x = random_sample((10, 10))
assert_allclose(expm(x), kaityo_expm(x))

这导致以下输出。

AssertionError: 
Not equal to tolerance rtol=1e-07, atol=0

Mismatch: 100%
Max absolute difference: 7.04655733
Max relative difference: 0.59875635
 x: array([[18.032424, 16.224408, 12.432163, 16.614248, 12.85653 , 13.705387,
        15.096966, 10.577946, 18.399573, 17.938062],
       [16.352809, 17.525898, 12.79079 , 16.295562, 13.512996, 14.407979,...
 y: array([[18.649103, 13.157682, 11.264763, 16.099163, 15.2293  , 17.854499,
        11.691586, 13.412066, 15.023189, 15.598455],
       [13.157682, 13.612502,  9.628261, 12.659313, 13.559437, 13.382417,..

显然，两种实现方式都不同。 问题如下：

• 他们可以接受不同吗？
• 我的实现有错吗？
• 如果我的实现有误，我该如何解决？
• 如果我的实现是正确的，什么时候可以安全使用scipy.linalg.expm？

我看到了以下问题：

• Matrix exponentiation with scipy: expm, expm2 and expm3

Answer 1

从数学方法上，矩阵的指数定义是使用指数的泰勒级数来定义的，所以：

设 A 为对角方阵：

当 A 是通用方阵时会出现问题，因此在进行指数运算之前，您需要使用特征值和特征向量对其进行对角化：

U 是特征向量矩阵，Lambda 是特征值在对角线上的矩阵。

在这一点上，我们接近于找到矩阵的指数：

现在让我们在一个简单的脚本中实现这个结果：

>>> import numpy as np
>>> import scipy.linalg as ln

>>> A = [[2/3, -4/3, 2],
         [5/6, 4/3, -2],
         [5/6, -2/3, 0]]
>>> A = np.matrix(A)
>>> print(A)
[[ 0.66666667 -1.33333333  2.        ]
 [ 0.83333333  1.33333333 -2.        ]
 [ 0.83333333 -0.66666667  0.        ]]

>>> eigvalue, eigvectors = np.linalg.eig(A)
>>> print("eigvalue: ", eigvalue)
>>> print("eigvectors:")
>>> print(eigvectors)
eigvalue:  [ 1. -1.  2.]
eigvectors:
[[ 0.81649658  0.27216553  0.87287156]
 [ 0.40824829 -0.68041382 -0.21821789]
 [ 0.40824829 -0.68041382  0.43643578]]

>>> e_Lambda = np.eye(np.size(A, 0))*(np.exp(eigvalue))
>>> print(e_Lambda)
[[2.71828183 0.         0.        ]
 [0.         0.36787944 0.        ]
 [0.         0.         7.3890561 ]]

>>> e_A = eigvectors*e_Lambda*eigvectors.I
>>> print(e_A)
[[ 2.3265481  -6.22769903  7.01116649]
 [ 0.97933433  4.27520659 -3.51559341]
 [ 0.97933433 -3.11384951  3.87346269]]

>>> e_A2 = ln.expm(A)
>>> print(e_A2)
[[ 2.3265481  -6.22769903  7.01116649]
 [ 0.97933433  4.27520659 -3.51559341]
 [ 0.97933433 -3.11384951  3.87346269]]

>>> np.testing.assert_allclose(e_A, e_A2)
>>> print(e_A - e_A2)
[[-1.77635684e-15  1.77635684e-15 -8.88178420e-16]
 [ 4.44089210e-16 -1.77635684e-15  8.88178420e-16]
 [-2.22044605e-16  0.00000000e+00  4.44089210e-16]]

我们看到结果基本相同，所以我认为使用scipy.linalg.expm 进行矩阵求幂是安全的。

我用笔记本创建了一个repo 以供进一步测试。