[SciPy-dev] roots, poly, comparison operators, etc.
eric
eric at scipy.org
Sat Feb 16 03:39:25 EST 2002
Hey crew,
I spent some quality time with roots() and poly() tonight trying to find a solution to the rounding issue acceptable to all. In the process, I've rewritten both of them. As Heiko suggested, it is easy to test if the roots are complex congugates and reals in poly(), and, if so, ignore the imaginary part. Things are never as easy as they seem...
Here is what I have learned:
roots() calculates eigenvalues under the covers. This is done using eig which in turn calls some LAPACK (Atlas) routine xgeev, where x is the appropriate numeric type. This is where the trouble starts. The values returned by eig() that should be complex conjugates are actually off by as much as 2-3 digits of precision. I've compared the output of SciPy to Matlab and Octave below.
# MATLAB
# restricts display to 14 digits of precision, 15 digits if you print in long e format
>> format long
>> roots([-1,1,2,1,2,3,4,5,6,7,8,9])
ans =
2.37248275302093
0.98438469660693 + 0.73802739091801i
0.98438469660693 - 0.73802739091801i
0.45252913029012 + 1.04491513664358i
0.45252913029012 - 1.04491513664358i
-1.14001283346666 + 0.26706945671060i
-1.14001283346666 - 0.26706945671060i
-0.77808817947927 + 0.77104779864671i
-0.77808817947927 - 0.77104779864671i
-0.20505419046158 + 1.06435796274941i
-0.20505419046158 - 1.06435796274941i
# OCTAVE has slightly different answers than MatLab, but complex conjugates are always equal.
octave:15> format long
octave:16> roots([-1,1,2,1,2,3,4,5,6,7,8,9])
ans =
2.372482753020929 + 0.000000000000000i
0.984384696606924 + 0.738027390918009i
0.984384696606924 - 0.738027390918009i
0.452529130290118 + 1.044915136643580i
0.452529130290118 - 1.044915136643580i
-1.140012833466655 + 0.267069456710605i
-1.140012833466655 - 0.267069456710605i
-0.778088179479270 + 0.771047798646705i
-0.778088179479270 - 0.771047798646705i
-0.205054190461583 + 1.064357962749411i
-0.205054190461583 - 1.064357962749411i
# SciPy -- neither the real or imag part of the shown complex pair are equal.
>>> r=roots.roots([-1,1,2,1,2,3,4,5,6,7,8,9])
>>> r[0]
(-0.2050541904615848+1.0643579627494104j)
>>> r[1]
(-0.20505419046158302-1.0643579627494113j)
So, Matlab, Octave generally give same results out to 15 digits and their complex conjugate pairs are always equal out to 15 digits. Python all give the same results up to 14 digits, but complex conjugates are off for the 15th digit. This causes them to fail as equal during comparisons.
Both Matlab and Octave equivalency tests for conjugate pairs pass. I'm not sure how these tools treat comparisons under the covers -- whether they only compare the first 14 or 15 digits or not. Python comparisons will obviously fail unless we do some kind of rounding here.
The fact that octave has 15 digits of precision and SciPy only has 14 on LAPACK output bugs me -- and makes things more difficult. I wonder if Atlas is less precise than the standard LAPACK? This needs some more investigation and really should be fixed. It was save a lot of headaches.
Other things:
*
Array printing is far superior in Matlab and Octave -- they generally always look nice. We should clean up how arrays are output. Also, the "format long" and "format short", etc options for specifying how arrays are printed are pretty nice.
*
Numeric not allowing the comparison of complex bit me again. When trying to compare complex conjugates, I had to write about 5 lines of code plus an extra sort function. The same thing is possible in a single line of Matlab/Octave code. In Numeric, sort() will not work for complex arrays, so you have to write your own sort() routines. I think we need to override sort() to work with complex. I feel the same way about comparison operators -- otherwise you have to special case *every* comparison operator in your code in order for it to work with generic arrays.
On Matlab/Octave, sort() sorts by magnitude of complex and then by angle. On the other hand ==, >, <, etc. seem to only compare the real part of complex numbers. These seem fine to me. I know they aren't mathematically correct, but they seem to be pragmatically correct. I'd like comments on these conventions and what others think.
*
Comparison of IEEE floats is hairy. It looks to me like Matlab and Octave have chosen to limit precision to 15 digits. This seems like a reasonable thing to do, for SciPy also, but we'd currently have to limit to 14 digits to deal with problems of imprecise LAPACK routines. Pearu and Fernando are arguing against this, but maybe they wouldn't mind a 15 digit limit for double and 6 for float. We'd have to modify Numeric to do this internally on comparisons. There could be a flag that enables exact comparisons.
*
After playing with Matlab/Octave code a little today, I remember how easy it is to do some things in a line or two that take many lines of Python code (the conjugate comparison comes to mind). I hope we can tune SciPy such that this difference disappears.
Here is the code I'm playing with:
from scipy import diag,eig,r1array,limits
from Numeric import *
def find_non_zero(a):
""" Mimics Matlab's find behavior.
"""
flattened = a.flat
return compress(flattened != 0,arange(len(flattened)))
def roots(p):
""" return the roots of the polynomial coefficients in p.
The values in the rank-1 array p are coefficients of a polynomial.
If the length of p is n+1 then the polynomial is
p[0] * x**n + p[1] * x**(n-1) + ... + p[n-1]*x + p[n]
"""
# If input is scalar, this makes it an array
p = r1array(p)
if len(p.shape) != 1:
raise ValueError,"Input must be a rank-1 array."
# there will be N-1 roots -- start with all of them as zeros
roots = zeros(len(p)-1,Complex)
# find non-zero array entries
non_zero = find_non_zero(p)
print non_zero
# find the number of trailing zeros -- this is the number of roots at 0.
# strip leading and trailing zeros
p = p[non_zero[0]:non_zero[-1]+1]
print 'p:', p
N = len(p)
if N > 1:
A = diag(ones((N-2,),'D'),-1)
A[0,:] = -p[1:] / (p[0]+0.0)
# fill in the end of the roots array with the new values
print A
print 'eig:',eig(A)
roots[-(N-1):] = eig(A)[0]
# sort roots from smallest magnitude to largest.
ind = argsort(abs(roots))
roots = take(roots,ind)
return roots
def sort_complex(a):
""" Doesn't currently work for integer arrays -- only float or complex.
"""
a = asarray(a,typecode=a.typecode().upper())
def complex_cmp(x,y):
res = cmp(x.real,y.real)
if res == 0:
res = cmp(x.imag,y.imag)
return res
l = a.tolist()
l.sort(complex_cmp)
return array(l)
def poly(seq_of_zeros):
"""Return a sequence representing a polynomial given a sequence of roots.
If the input is a matrix, return the characteristic polynomial.
Example:
>>> b = roots([1,3,1,5,6])
>>> poly(b)
array([1., 3., 1., 5., 6.])
"""
seq_of_zeros = r1array(seq_of_zeros)
if len(seq_of_zeros) == 0:
return 1.0
if len(seq_of_zeros.shape) == 2:
seq_of_zeros, vecs = MLab.eig(seq_of_zeros)
a = [1]
for k in range(len(seq_of_zeros)):
a = convolve(a,[1, -seq_of_zeros[k]], mode=2)
try:
output_type = seq_of_zeros.typecode().upper()
zero_limit = limits.epsilon(output_type.lower()) * 50 # rather arbitrarily chosen.
roots = asarray(seq_of_zeros,output_type)
pos_roots = sort_complex(compress(roots.imag > 0,roots))
neg_roots = conjugate(sort_complex(compress(roots.imag < 0,roots)))
print neg_roots.real - pos_roots.real
print neg_roots.imag - pos_roots.imag
if (alltrue(abs(neg_roots.real - pos_roots.real) < eps * abs(neg_roots.real + pos_roots.real)) and
alltrue(abs(neg_roots.imag - pos_roots.imag) < eps * abs(neg_roots.imag + pos_roots.imag))):
a = a.real
except ValueError:
# neg_roots and pos_roots had different number of entries.
pass
return a
eric
--
Eric Jones <eric at enthought.com>
Enthought, Inc. [www.enthought.com and www.scipy.org]
(512) 536-1057
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20020216/7fdecfc5/attachment.html>
More information about the SciPy-Dev
mailing list