Получение комплексных коэффициентов в ближайших SPD-матрицах

Question

Я пишу python 3.7 программу, и мне нужно получить симметричные c положительно определенные матрицы.

Я использовал этот код, чтобы получить ближайшее SPD (все собственные значения должны быть> 0) :

Python: преобразовать матрицу в положительную полуопределенную

Дело в том, что мне нужно вычислить римановы экспоненты симметрии c матрицы: Определение римановой экспоненты .

Но я получаю матрицы с комплексными коэффициентами. Как я могу избавиться от этого?

Примечание: Даже с вышеупомянутой реализацией я получаю матрицы с комплексными числами.

Я также пытался исследовать пакет geomstats, но я не делаю знать, как его использовать:

https://geomstats.github.io/geometry.html#module -geomstats.geometry.spd_matrices

Большое спасибо

Редактировать 1: мой код и что я ожидаю:

Это моя функция:

def expNearestSPD(Y, Delta):

    """
    Implementation of riemannian exponential using the nearestPD function
    I try to always keep SPD matrices
    But at the end it does not work
    """

    Y2 = nearestPD(Y)

    mul = np.matmul(  np.matmul(inv(sqrtm(Y2)), Delta), inv(sqrtm(Y2))  ) 

    mul = nearestPD(mul)

    z1 = expm(   mul   )

    z1 = nearestPD(z1)

    z = np.matmul(   np.matmul(sqrtm(Y2), z1), sqrtm(Y2)   )

    return nearestPD(z)

Например, вот P, матрица SPD:

P

array([[0.1, 0. , 0. , ..., 0. , 0. , 0. ],
       [0. , 0.1, 0. , ..., 0. , 0. , 0. ],
       [0. , 0. , 0.1, ..., 0. , 0. , 0. ],
       ...,
       [0. , 0. , 0. , ..., 0.1, 0. , 0. ],
       [0. , 0. , 0. , ..., 0. , 0.1, 0. ],
       [0. , 0. , 0. , ..., 0. , 0. , 0.1]])

Можно проверить:

np.isrealobj(P)

True

Но когда я вычисляю, например, expNearestSPD(P, P), я получаю:

expNearestSPD(P, P)

array([[0.27182818+0.j, 0.        +0.j, 0.        +0.j, ...,
        0.        +0.j, 0.        +0.j, 0.        +0.j],
       [0.        +0.j, 0.27182818+0.j, 0.        +0.j, ...,
        0.        +0.j, 0.        +0.j, 0.        +0.j],
       [0.        +0.j, 0.        +0.j, 0.27182818+0.j, ...,
        0.        +0.j, 0.        +0.j, 0.        +0.j],
       ...,
       [0.        +0.j, 0.        +0.j, 0.        +0.j, ...,
        0.27182818+0.j, 0.        +0.j, 0.        +0.j],
       [0.        +0.j, 0.        +0.j, 0.        +0.j, ...,
        0.        +0.j, 0.27182818+0.j, 0.        +0.j],
       [0.        +0.j, 0.        +0.j, 0.        +0.j, ...,
        0.        +0.j, 0.        +0.j, 0.27182818+0.j]])

Я получаю сложные, но очень маленькие коэффициенты.

Answer 1

Фактически, используя пакет geomstats, я нашел решение.

Это исключительно для SPD-матриц.

Есть функция для непосредственного вычисления римановой экспоненты, и даже одна для вычислить A**t для A SPD и t real.

Я рекомендую использовать последнюю версию geomstats (из GitHub).

Здесь N - это порядок квадратные матрицы (так что их размеры N x N).

import geomstats.backend as gs
from geomstats.geometry.spd_matrices import SPDMatrices, SPDMetricAffine

gs.random.seed(0)
space = SPDMatrices(N)
metric = SPDMetricAffine(N)

def Exp(Y, Delta):
    return metric.exp(Delta, Y)

# for the power of a SPD matrix

gs.linalg.powerm(M, t)