Re: Patch to residue.m

[Doug Stewart]

Daniel J Sebald wrote:
> Ben Abbott wrote:
>
> > I don't find the part concerning zero, pure real, pure imaginary 
> > poles  well thought out. However, since it was there before, and 
> > others have  expressed a desire for such checks, I've left it alone. 
> > In any event,  the tolerance is small enough that I don't think it 
> > makes any  numerical difference.
> >
> > > Now, here is the part I'm not so comfortable with.  There is the  
> > > code above that says a pole is zero if it is within ~1e-12 of 
> > > zero.   But shortly thereafter the code does this:
> > >
> > > [e, indx] = mpoles (p, toler, 1);
> > >
> > > where poles within "toler" (in the code it is set to 0.001) are  
> > > equal.  One case is a more strict requirement to be declared 0 than  
> > > the other case of declaring two poles equal.  It's just slightly  
> > > strange.  Doesn't that seem a contradiction somehow?
> >
> >
> > The problem is that as multiplicity increases, the location of the  
> > poles become numerically difficult to determine. When the 
> > multiplicity  is unity, a relative tolerance of 1e-12 is possible. 
> > However, as the  multiplicity increases, numerical errors increase as 
> > well. The values  of "toler" and "small" were determined by empirical 
> > means (experience  and testing). I tried to make them as small as was 
> > comfortable.
> >
> > However, I must admit your questions have made me re-examine this  
> > subject, and it does appear that the checks for zero valued, real  
> > valued, and imaginary valued poles  should be done *after* the  
> > multiplicity has been determined ... if they are to be done at all.
>
> Perhaps that is an issue.
>
>
> > > I'm not sure I'm comfortable with averaging the poles either:
> > >
> > >   p(m) = mean (p(m));
> > >
> > > From a practical standpoint it may be a good idea, but from a  
> > > numerical standpoint it seems a little unorthodox.  I'd say let the  
> > > numerical effects remain and let the end user decide what to do.   
> > > Why is this needed?
> >
> >
> > I'll give a very simple example.
> >
> > octave:20> p = [3;3;3];
> > octave:21> a = poly(p)
> > a =
> >
> >     1   -9   27  -27
> >
> > octave:22> roots(a)
> > ans =
> >
> >    3.0000 + 0.0000i
> >    3.0000 + 0.0000i
> >    3.0000 - 0.0000i
> >
> > octave:23> roots(a)-p
> > ans =
> >
> >   3.3231e-05 + 0.0000e+00i
> >   -1.6615e-05 + 2.8778e-05i
> >   -1.6615e-05 - 2.8778e-05i
> >
> > octave:24> mean(roots(a))-p
> > ans =
> >
> >    1.7764e-15
> >    1.7764e-15
> >    1.7764e-15
> >
> > Thus the error from "roots" is about 1e-5 for each root. However, 
> > the  error for the average is about 10 orders of magnitude less!
>
> Well, the improvement from averaging probably isn't a surprise.  (And 
> one might be able to look at the rooting formula and analytically show 
> that the average is a good idea.)  The idea of then casting the poles 
> is what worries me because maybe the user is dealing with a situation 
> where he or she wouldn't want that... but...

This is  "numerical  method"  math and some times there is no exact 
answer. see below.


> > To be honest, I began from your position and was sure that taking 
> > the  mean would not ensure an improvement in the result. Doug was 
> > certain  otherwise. So I decided to justify my skepticism by 
> > designing an  example, where I stacked the cards in my favor.
> >
> > To my surprise, I failed. All tests confirmed Doug's conclusion. If  
> > you can devise an example that succeeds where I failed, I'd love to  
> > see it!
> >
> > >  I suppose you are dealing with finding higher order residues,
> >
> >
> > Yes.
> >
> > > but could you compute the mean value for use in the rest of the  
> > > routine but leave the values non-averaged when returned from the  
> > > routine?  What is the user expecting?
> >
> >
> > No, at least I don't think so. The residues and poles are uniquely  
> > paired. The returned residues should be paired with the returned  
> > poles. Regarding what the user expects, a better question would be  
> > what does the documentation reflect, and is it consistent with 
> > Matlab.  Regarding the returned variable, they are consistent with 
> > the  documentation and with Matlab.
>
>
> ...but, now that I think of it, if you are saying there is a root of 
> greater than one multiplicity I guess it makes sense that those poles 
> should come back as being exactly the same value.
Yes this is the point


> > > I guess my main concerns are 1) is there a standard way of dealing  
> > > with "nearly zero" in Octave?  Should "toler" be replaced by "small"?
> >
> >
> > Regarding "near zero" ... I don't know how to properly determine 
> > what  is "nearly zero", but would appreciate anyone telling us/me how.
>
> I'm trying to remember back to whether there was ever a discussion 
> here about this topic, i.e., when is something considered zero.
>
> > Regarding, "toler" and "small" these two have two different 
> > purposes.  "toler" is used to determine the pole multiplicity and 
> > "small" is used  to round off to zero, pure real, or pure imaginary 
> > values. If "toler"  is replace by "small" there will be no 
> > multiplicity of poles. So I  don't think that is practical.
>
> Well, that doesn't mean that small = toler = 1e-12.  Perhaps some 
> other value.  I'm just trying to justify in my mind the need for two 
> tolerance parameters vs. one.
>
>
> > What we really need is numerically more accurate method for  
> > determining the roots. This entire discussion is the consequence of  
> > inaccurate roots. So if we are to address the ills of residue.m, a  
> > better roots.m will be needed, in my opinion anyway.
>
> I think you are right.  The roots being off by o(1e-5) seems rather 
> poor.  Averaging gets to o(1e-15) and machine precision is o(1e-16)?  
> One would think maybe o(1e-10) or better is achievable.  Maybe you are 
> trying to improve on something that should be addressed elsewhere, and 
> if that where fixed your tolerance equations could be simplified.
>
> And considering the cast to zero after multiplicity is determined 
> seems worth investigating.
>
> Dan
For a quadratic there is an exact answer, and for cubics and for 
quartics:    but there is no exact method for 5th and higher polynomials.
Using residue one should recognize that it is a numerical approach to 
find approximate roots and residues. And if you recognize that they are 
approximate then the values used for tol and small seem appropriate. 
Yes, if you are looking for roots in the 10e-12 range  then you you will 
have problems  - but this was true before the improvements.

If you want exact answers try Maxima, I did and found that it didn't 
even know how to do the exact 3 and 4th order answers..

If someone would write a better roots.m file we will all benefit. But  
it is not an easy task. ( I have been thinking about it since 
"Numerical  method" course in university in 1982) Using VPA (Variable 
Precision Arithmetic) yes you can get the roots to the nth decimal point 
but at a cost of cpu time-- so the conclusion that I have come to is:    
that this is a step in the right direction and is an improvement over 
what was here before.




Doug