Finding eigenvalues of log-represented matrix

Question

I am currently trying to find the eigenvalues of very poorly conditioned real symmetric matrices, whose elements can be extremely small/large (from 1e-8 to 1e15) and possibly negative. The matrices are fairly small though (4x4), so speed of execution is not a major concern.

All operations before diagonalization are done using variants of the "logsumexp" trick (adapted to matrix-matrix / matrix-vector multiplication from this thread numerically stable way to multiply log probability matrices in numpy, and further to accomodate negative coefficients), so I end up with two matrices (sgnM, logM) containing respectively the sign and log of absolute value of the elements in M. This part works well.

However, I did not find any documentation for an eigvalsh equivalent that would take this as an input, and keep using these numerical tricks all the way until it returns the eigenvalues. For now, I simply use

scipy.linalg.eigvalsh(sgnM*np.exp(logM))

which still suffers from precision issues (the eigenvalue that I am interested in goes from 1e-4 to around 1e-9 when I change a parameter, and I can see that the ones around 1e-9 are much noisier than the first ones to the point that the results derived from there stop making sense).

Is there some function ressembling what I am looking for hidden in an existing linear algebra engine ? If not, can I still rely on some of Lapack/MKL/Blas routines to implement it or does it need to be done from the ground up?

score 0 · Accepted Answer · answered Nov 28 '19 at 12:30

I managed to implement it, unfortunately it did not give much improvement in my case. I still leave this here in case someone wants to try it on another problem (this is based on code from https://github.com/mateuv/MetodosNumericos/blob/master/python/NumericalMethodsInEngineeringWithPython/jacobi.py:

import numpy as np
from scipy.special import logsumexp


def eigsh_jacobi_stable(log_a, sgn_a, tol = 1.0e-14): # Jacobi method

    def maxElem(log_a): # Find largest off-diag. element a[k,l]
        n = len(log_a)
        log_aMax = np.log(1e-36)
        for i in range(n-1):
            for j in range(i+1,n):
                if log_a[i,j] >= log_aMax:
                    log_aMax = log_a[i,j]
                    k = i; l = j
        return log_aMax,k,l

    def rotate(log_a, sgn_a, log_p, sgn_p, k, l): # Rotate to make a[k,l] = 0
        n = len(log_a)
        log_aDiff, sgn_aDiff = logsumexp(a=[log_a[l,l], log_a[k,k]],
                                         b=[sgn_a[l,l], -sgn_a[k,k]],
                                         return_sign=True)

        if log_a[k,l] - log_aDiff < np.log(1.0e-36):
            log_t = log_a[k,l] - log_aDiff
            sgn_t = sgn_a[k,l] * sgn_aDiff
        else:
            log_phi = log_aDiff - np.log(2) - log_a[k,l]
            sgn_phi = sgn_aDiff * sgn_a[k,l]
            log_t = -logsumexp(a=[log_phi, .5*logsumexp(a=[2*log_phi, 0])])
            sgn_t = 1.
            if sgn_phi < 0:
                sgn_t = -1

        log_c = - .5 * logsumexp(a=[2*log_t, 0])
        sgn_c = 1.
        log_s = log_t + log_c
        sgn_s = sgn_t * sgn_c
        log_tau = log_s - logsumexp(a=[log_c, 0])
        sgn_tau = sgn_s

        log_temp, sgn_temp = log_a[k,l], sgn_a[k,l]
        log_a[k,l] = np.log(1.0e-36)
        sgn_a[k, l] = 0
        log_a[k,k], sgn_a[k,k] = logsumexp(a=[log_a[k,k], log_t + log_temp],
                                           b=[sgn_a[k,k], -sgn_t * sgn_temp],
                                           return_sign=True)
        log_a[l,l], sgn_a[l,l] = logsumexp(a=[log_a[l,l], log_t + log_temp],
                                           b=[sgn_a[l,l], sgn_t * sgn_temp],
                                           return_sign=True)

        for i in range(k):      # Case of i < k
            log_temp, sgn_temp = log_a[i,k], sgn_a[i,k]
            log_plop, sgn_plop = logsumexp(a=[log_a[i, l], log_tau + log_temp],
                                           b=[sgn_a[i, l], sgn_tau * sgn_temp],
                                           return_sign=True)
            log_a[i,k], sgn_a[i,k] = logsumexp(a=[log_temp, log_s + log_plop],
                                               b=[sgn_temp, -sgn_plop * sgn_s],
                                               return_sign=True)
            log_plop, sgn_plop = logsumexp(a=[log_temp, log_tau + log_a[i, l]],
                                           b=[sgn_temp, -sgn_tau * sgn_a[i, l]],
                                           return_sign=True)
            log_a[i,l], sgn_a[i,l] = logsumexp(a=[log_a[i, l], log_s + log_plop],
                                               b=[sgn_a[i, l], sgn_plop * sgn_s],
                                               return_sign=True)

        for i in range(k+1,l):  # Case of k < i < l
            log_temp, sgn_temp = log_a[k, i], sgn_a[k,i]
            log_plop, sgn_plop = logsumexp(a=[log_a[i, l], log_tau + log_a[k,i]],
                                           b=[sgn_a[i, l], sgn_tau * sgn_a[k,i]],
                                           return_sign=True)
            log_a[k, i], sgn_a[k, i] = logsumexp(a=[log_temp, log_s + log_plop],
                                                 b=[sgn_temp, -sgn_plop * sgn_s],
                                                 return_sign=True)
            log_plop, sgn_plop = logsumexp(a=[log_temp, log_tau + log_a[i, l]],
                                           b=[sgn_temp, -sgn_tau * sgn_a[i, l]],
                                           return_sign=True)
            log_a[i,l], sgn_a[i,l] = logsumexp(a=[log_a[i, l], log_s + log_plop],
                                               b=[sgn_a[i, l], sgn_plop * sgn_s],
                                               return_sign=True)
        for i in range(l+1,n):  # Case of i > l
            log_temp, sgn_temp = log_a[k, i], sgn_a[k,i]
            log_plop, sgn_plop = logsumexp(a=[log_a[l,i], log_tau + log_temp],
                                           b=[sgn_a[l,i], sgn_tau * sgn_temp],
                                           return_sign=True)
            log_a[k, i], sgn_a[k, i] = logsumexp(a=[log_temp, log_s + log_plop],
                                                 b=[sgn_temp, -sgn_plop * sgn_s],
                                                 return_sign=True)
            log_plop, sgn_plop = logsumexp(a=[log_temp, log_tau + log_a[l,i]],
                                           b=[sgn_temp, -sgn_tau * sgn_a[l,i]],
                                           return_sign=True)
            log_a[l,i], sgn_a[l,i] = logsumexp(a=[log_a[l,i], log_s + log_plop],
                                               b=[sgn_a[l, i], sgn_plop * sgn_s],
                                               return_sign=True)
        for i in range(n):      # Update transformation matrix
            log_temp, sgn_temp = log_p[i,k], sgn_p[i,k]
            log_plop, sgn_plop = logsumexp(a=[log_p[i,l], log_tau + log_p[i,k]],
                                           b=[sgn_p[i,l], sgn_tau * sgn_p[i,k]],
                                           return_sign=True)
            log_p[k, i], sgn_p[k, i] = logsumexp(a=[log_temp, log_s + log_plop],
                                                 b=[sgn_temp, -sgn_plop * sgn_s],
                                                 return_sign=True)
            log_plop, sgn_plop = logsumexp(a=[log_temp, log_tau + log_p[i, l]],
                                           b=[sgn_temp, -sgn_tau * sgn_p[i, l]],
                                           return_sign=True)
            log_p[i,l], sgn_p[i,l] = logsumexp(a=[log_p[i, l], log_s + log_plop],
                                               b=[sgn_p[i, l], sgn_plop * sgn_s],
                                               return_sign=True)

    n = len(log_a)
    maxRot = 5*(n**2)       # Set limit on number of rotations
    log_p, sgn_p = np.log(1e-36) * np.ones_like(log_a), np.zeros_like(log_a)
    np.fill_diagonal(log_p, 0)
    np.fill_diagonal(sgn_p, 1)

    for i in range(maxRot): # Jacobi rotation loop
        log_aMax,k,l = maxElem(log_a)
        # print(log_aMax)
        if np.exp(log_aMax) < tol:
            return np.sort(np.exp(np.diag(log_a)) * np.diag(sgn_a))
        rotate(log_a, sgn_a, log_p, sgn_p ,k,l)
    raise RuntimeError('Jacobi method did not converge')



if __name__ == '__main__':
    from scipy.linalg import eigvalsh

    for _ in range(100):
        test = np.random.randn(10, 10)
        test += test.transpose()

        log_test, sgn_test = np.log(np.abs(test)), np.sign(test)
        eigs_jacobi = np.sort(eigsh_jacobi_stable(log_test, sgn_test)) # this sort is unnecessary, already done internally
        eigs_scipy = np.sort(eigvalsh(test))
        print(eigs_jacobi-eigs_scipy)

Finding eigenvalues of log-represented matrix

1 Answers1