2

There appears to be a major difference between numpy's slogdet and the exact result when computing the log determinant of Vanermonde matrix.

I compare against the exact log determinant, see eg here for proof.

The minimal code to see this is:

A = np.power.outer(np.linspace(0,1,50),range(50))

print np.linalg.slogdet(A)[1]

s = 0
for v1 in np.linspace(0,1,50):
    for v2 in np.linspace(0,1,50):
        if v1>v2:
            s+= np.log(v1-v2)

print s

Which yeilds:

-1191.88408998
-1706.99560647

I was wondering if there was a more accurate log determinant implementation which I could use in this situation but also in non-Vandermonde matrix situation.

j__
  • 248
  • 1
  • 11
  • For N=10 the values match, as does `np.log(np.det(A))`. It could be that by N=30, the `det(A)` is so small that even `slogdet` can't cope. – hpaulj Jun 08 '17 at 17:35

1 Answers1

0

You can use sympy and mpmath like this:

import numpy as np
import sympy as smp
import mpmath as mp

mp.mp.dps = 50 

linspace1 = list(map(smp.mpmath.mpf,np.linspace(0,1,50)))
A = np.power.outer(list(map(float,linspace1)),range(50))

first_print = smp.mpmath.mpf(np.linalg.slogdet(A)[1])
print(first_print)

s = 0
linspace2 = list(map(smp.mpmath.mpf,np.linspace(0,1,50)))
linspace3 = list(map(smp.mpmath.mpf,np.linspace(0,1,50)))
for v1 in linspace1:
    for v2 in linspace2:
        if v1>v2:
            s+= mp.log(v1-v2)

print(s)

RESULTS

first_print = -1178.272517342130186079884879291057586669921875

s = -1706.9956064674289001970168329846189154212781094939
developer_hatch
  • 15,898
  • 3
  • 42
  • 75
  • Thanks but these two numbers are still off by roughly 50%, any ideas why? – j__ Jun 08 '17 at 17:00
  • @j__ what is the value expected aprox? – developer_hatch Jun 08 '17 at 17:04
  • 1706.995 is the actual value using the formula specific to vandermonde matrices, but slogdet is roughly 50% off with 1706.272 – j__ Jun 08 '17 at 17:07
  • @j__ I updated the answer again, let me know if those values are correct! – developer_hatch Jun 08 '17 at 17:42
  • 1
    I don't think the OP is complaining about the iterative solution; it's the `slogdet` that's off. It's faster but doesn't handle the larger array with very small determinant. And `mpmath` does nothing for the `slogdet`. That's calculated with a compiled LAPACK function using float64. – hpaulj Jun 08 '17 at 17:48
  • @hpaulj so there is no obvious solution? Thanks for the comments – j__ Jun 08 '17 at 17:53
  • You could try the [linear algebra routines of the mpmath library](http://mpmath.org/doc/current/matrices.html#the-eigenvalue-problem) instead the numpy ones. – user8153 Jun 08 '17 at 19:37