0

I am writing a python 3.7 program and I need to get symmetric positive-definite matrices.

I used this code to get the nearest SPD (all eigenvalues have to be > 0) :

Python: convert matrix to positive semi-definite

The fact is that I have to compute riemannian exponentials of symmetric matrices: Definition of the riemannian exponential.

But I get matrices with complex coefficients. How could I get rid of that ?

Note: Even with the aforementioned implementation, I get matrices with complex numbers.

I also tried to explore the geomstats package, but I do not know how to use it :

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

Thanks a lot

Edit 1: my code and what I expect:

This is my function:

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)

For example, here is P, a SPD matrix:

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]])

One can check:

np.isrealobj(P)

True

But when I compute expNearestSPD(P, P) for example, I get:

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]])

I get complex, but very little, coefficients.

Deep Learner
  • 1
  • 2
  • Could you post a sample matrix and expected results of what you need to get out of it? I think it will make a difference if all or only some values have complex coefficients. – Ethan Apr 01 '20 at 19:08
  • Ok, I will post that. The thing is that, _in theory_, I should always keep SPD matrices. But I does not happen (don't know how python makes computations) and I get (quite close to 0) complex coefficients. I think the inversions are the cause of that. Thanks for your help. Most of the values are complex, but the imaginary parts are close to 0. – Deep Learner Apr 02 '20 at 13:09
  • I’m having trouble executing the code. Are expm() and sqrtm() custom functions? If so, can you post them? – Ethan Apr 03 '20 at 15:35
  • They come from scipy.linalg . They are respectively the exponential and square root for matrices. Sorry for not having precised this. – Deep Learner Apr 06 '20 at 08:29
  • Sorry, I confused a couple of different posts I was looking at! – Ethan Apr 06 '20 at 10:54

1 Answers1

0

In fact, using the package geomstats, I found a solution.

This is exclusively for SPD matrices.

There is a function to directly compute the riemannian exponential, and even one to compute A**t for A SPD and t real.

I recommend using the last version of geomstats (from GitHub).

Here, N is the order of the square matrices (so that their sizes are 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)
Deep Learner
  • 1
  • 2