Re: Patch to residue.m

[Daniel J Sebald]

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...

> 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.


> > 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