Разница между scipy.linalg.expm и ручной кодировкой

Question

Я пытался реализовать матричную экспоненциальную функцию, как в scipy.linalg.expm. Я черпал вдохновение из репозитория github kaityo256 . Таким образом, я записал следующее.

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?

Я видел следующие вопросы:

• Матричное возведение в степень со scipy: expm, expm2 и expm3

Answer 1

из математического подхода определение экспоненты матрицы производится с использованием ряда экспоненты Тейлора, поэтому:

пусть A будет диагональной квадратной матрицей:

проблема возникает, когда A является общей c квадратной матрицей, поэтому перед выполнением экспоненты вам нужно сделать диагонализацию, используя собственные значения и собственные векторы:

с U матрицей собственных векторов и лямбда-матрицей с собственными значениями на диагонали.

В этой точке мы близки к нахождению экспоненциальной матрицы:

теперь давайте реализуем этот результат в простом сценарии:

>>> 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 для возведения в степень матрицы.

Я создал репо с ноутбуком для дальнейшего тестирования.