BUG scipy.linalg.lapack: potrf/ptroi interpret their 'lower' argument differently

[untom]

The second argument to these functions should define if the cholesky decomposition works on the lower/upper part of the matrix, however something doesn't work out.

This is a simple test to try it out:

import numpy as np
import scipy.linalg as la
x=np.random.normal(size=(3, 3))
a = x.dot(x.T)

The normal inverse of this would be:

>>> la.inv(a)
array([[ 0.21546368,  0.01704477, -0.12058592],
       [ 0.01704477,  0.06994512, -0.01658643],
       [-0.12058592, -0.01658643,  0.11951102]])

Using dprotf/dproti:

>>> dpotrf, dpotri = la.lapack.get_lapack_funcs(("potrf", "potri"), (a, ))
>>> c, info = dpotrf(a, False, False)
>>> dpotri(c, False, True)
(array([[ 0.09379177, -0.01492376,  3.29255627],
       [ 0.        ,  0.06764316,  0.53362104],
       [-0.        , -0.        ,  0.11951102]]), 0)

Clearly, this is not the same as the answer given above. In fact, the answer is only exact when the 2nd argument to dpotrf and dpotri are different:

>>> c, info = dpotrf(a, True, True)
>>> np.abs(dpotri(c, True, True)[0] - np.triu(la.inv(a))).max()
3.2925562689500336

>>> c, info = dpotrf(a, False, True)
>>> np.abs(dpotri(c, True, True)[0] - np.triu(la.inv(a))).max()
5.5511151231257827e-17

>>> c, info = dpotrf(a, True, True)
>>> np.abs(dpotri(c, False, True)[0] - np.tril(la.inv(a))).max()
8.3266726846886741e-17

>>> c, info = dpotrf(a, False, True)
>>> np.abs(dpotri(c, False, True)[0] - np.tril(la.inv(a))).max()
3.2925562689500336

I tried this with SciPy 0.11 and 0.12, the underlying BLAS library was OpenBLAS.

[pv]

Sounds like get_lapack_funcs is getting one of them from clapack and the other from flapack. check what the __module__ attribute says (IIRC, check dir(dpotrf))

[ev-br]

c, info = dpotrf(a, True, True) is being interpreted as lower=True and clean=True
while in dp_inv, _ = dpotri(c, True, True) the 2nd argument is lower and the 3rd one is rowmajor.

So, setting rowmajor explicitly seems to work:

>>> c, info = dpotrf(a, True, True, rowmajor=True)
>>> print  np.abs(dpotri(c, True, True)[0] - np.tril(la.inv(a))).max()
2.41584530158e-13

[pv]

If you look into flapack.pyf.src, you see the if (clean) ... stuff in callstatement --- there's probably bug there, it seems the behavior is quite broken as it works only for potrf(lower=False, clean=True) & potri(lower=True)

This doesn't seem to be an issue with clapack/flapack conflict, but check what your dpotrf.module_name and dpotri.module_name say

The rowmajor stuff is clapack-only stuff; we should actually just remove the potr* routines from clapack.pyf

[ev-br]

>>>  dpotrf.module_name, dpotri.module_name
clapack clapack

Leaving the 3rd argument out completely works fine --- looks like a bug in clean stuff indeed. Does one actually need this clean thing?