ENH: use rotated companion matrix to reduce error #13202
Conversation
|
It seems to me that this would depend heavily on the linear algebra library used under the hood. |
|
That's interesting, I get the same general result but even better accuracy. What linear algebra library are you linked with? Be interesting to try with Chebyshev polynomials. I wonder if it is a Fortran artifact? |
|
I am using OpenBLAS. It could be related to a starting with an upper hessian vs lower hessian matrix. I think most eigensolvers start by converting the matrix to either upper or lower hessian. |
|
I've looked at the LAPACK routines and they reduce the matrix to upper Hessenberg form in Fortran, our matrices are lower Hessenberg from the Fortran point of view. If that were the major cause of the difference, then using the transpose should also improve things, and it does, but not to the same extent. So there seems to be is something special about the Hessenberg eigenvalue solver. |
|
Could you also do the remaining 5 polynomial types? The all seem to benefit from this change. |
|
It would also be good to add a comment that the rotated matrix gives more accurate results. |
c3a69d3 to
732d52a
Compare
|
Here are the results and the code for rotating with the other bases. Most results were improved by similar margins, except that the hermite and hermite_e polynomials showed little or no improvement. import numpy as np
import numpy.polynomial
from numpy.polynomial.polynomial import polyroots, polycompanion, polyfromroots
from numpy.polynomial.chebyshev import chebcompanion, chebroots, chebfromroots
from numpy.polynomial.hermite import hermcompanion, hermroots, hermfromroots
from numpy.polynomial.hermite_e import hermecompanion, hermeroots, hermefromroots
from numpy.polynomial.laguerre import lagcompanion, lagroots, lagfromroots
from numpy.polynomial.legendre import legcompanion, legroots, legfromroots
from numpy.polynomial import polyutils as pu
from numpy import linalg as la
def lagroots2(c):
[c] = pu.as_series([c])
if len(c) <= 1:
return np.array([], dtype=c.dtype)
if len(c) == 2:
return np.array([1 + c[0]/c[1]])
m = lagcompanion(c)
r = la.eigvals(m[::-1,::-1])
r.sort()
return r
def hermeroots2(c):
[c] = pu.as_series([c])
if len(c) <= 1:
return np.array([], dtype=c.dtype)
if len(c) == 2:
return np.array([-c[0]/c[1]])
# rotated companion matrix reduces error
m = hermecompanion(c)[::-1,::-1]
r = la.eigvals(m)
r.sort()
return r
def legroots2(c):
[c] = pu.as_series([c])
if len(c) < 2:
return np.array([], dtype=c.dtype)
if len(c) == 2:
return np.array([-c[0]/c[1]])
m = legcompanion(c)
r = la.eigvals(m[::-1,::-1])
r.sort()
return r
def hermroots2(c):
[c] = pu.as_series([c])
if len(c) <= 1:
return np.array([], dtype=c.dtype)
if len(c) == 2:
return np.array([-.5*c[0]/c[1]])
m = hermcompanion(c)
r = la.eigvals(m[::-1,::-1])
r.sort()
return r
def chebroots2(c):
[c] = pu.as_series([c])
if len(c) < 2:
return np.array([], dtype=c.dtype)
if len(c) == 2:
return np.array([-c[0]/c[1]])
m = chebcompanion(c)
r = la.eigvals(m[::-1,::-1])
r.sort()
return r
def polyroots2(c):
# c is a trimmed copy
[c] = pu.as_series([c])
if len(c) < 2:
return np.array([], dtype=c.dtype)
if len(c) == 2:
return np.array([-c[0]/c[1]])
m = polycompanion(c)
r = la.eigvals(m[::-1,::-1])
r.sort()
return r
def sort_diff(a,b):
return np.array(sorted(a)) - np.array(sorted(b))
k = 1
deg = 30
seed = 786
unif = np.random.uniform
def error(solver, maker):
err = []
np.random.seed(seed)
for _ in range(1000):
roots = unif(-k,k,size=deg) + 1j*unif(-k,k,size=deg)
c = maker(roots)
r = solver(c)
err.extend(sort_diff(roots, r).tolist())
return err
def test(solver_orig, solver2, maker):
err1 = error(solver_orig, maker)
err2 = error(solver2, maker)
print('\t\tAbsolute Error')
print("Original mean: {:1.3e}".format(np.abs(err1).mean()))
print("Rotated mean : {:1.3e}".format(np.abs(err2).mean()))
print()
print("Original max : {:1.3e}".format(np.abs(err1).max()))
print("Rotated max : {:1.3e}".format(np.abs(err2).max()))
print("Poly")
test(polyroots, polyroots2, polyfromroots)
print("Cheb")
test(chebroots, chebroots2, chebfromroots)
print("Lag")
test(lagroots, lagroots2, lagfromroots)
print("Herm")
test(hermroots, hermroots2, hermfromroots)
print("Herm_e")
test(hermeroots, hermeroots2, hermefromroots)
print("Leg")
test(legroots, legroots2, legfromroots) |
|
Thanks. Someday I need to look into why the generated polynomials for the other polynomial types are so bad. There has got to be a better way of generating the polynomials. ISTR getting good results for Chebyshev and equally spaced points in [-1, 1]. I expected bad results for the polynomials not confined to the interval [-1, 1], but still... |
|
I do wonder if some of the problem is a bad ordering of the resulting zeros, although 30 zeros really should have some decent separation and the sorting should work just fine. Of course, something like Chebyshev might be expected to produce a good uniform approximation, which doesn't mean the roots are well approximated. But I suspect something can be improved. |
I work with a research team on developing a higher dimensional root finding algorithm. We happened to learn that in the one-dimensional case, using the rotated companion matrix produces smaller errors.
Here is some sample code and the results.