Patch to residue.m

[Ben Abbott]

On Dec 31, 2007, at 12:04 PM, Daniel J Sebald wrote:

> Ben Abbott wrote:
>> Please consider this patch to residue.m which improves its  
>> accuracy  when pole multiplicity exists.
>> To capture the improvement, the test scripts were also modified to  
>> be  more stringent.
>
> Ben,
>
> This looks like not too much of a change.  Still, a few questions.   
> Why the change?  What motivated this?  Just improving your code?  Or  
> did you find an example that was failing?  (If so, give us the  
> example.)
>
> Observations:
>
> 1) Extra space in "group of pole  multiplicity,"
>
> 2)  p_group = cumsum (e == 1);
>   for ng = 1:max(p_group)
>
> Since (e == 1) is strictly non-negative, rather than max(p_group) it  
> could be p_group(end).  That is more efficient, but of course the  
> polynomial order is quite small so it doesn't really matter.

I'd not come across the "end" functionality before. Very cool! ... I  
was unhappy with p_group(numel(p_group)).

>
> 3)  Could you explain the reasoning of the various tolerance  
> parameters in the code?  For example, I see this code for  
> determining if a pole is zero or not:
>
> small = max (abs (p));
> small = max ([small, 1] ) * 1e-12 * (1 + numel (p))^2;
> p(abs (p) < small) = 0;
>
> I'm don't understand what the (1 + numel (p))^2 does.  The more  
> poles there are, the larger the tolerance becomes?  Also, might you  
> want to define 1e-12 with respect to the machine precision, say
>
> small = max ([small, 1]) * eps*1e4 * (1 + numel (p))^2;

Good point.

> (On my machine,
>
> octave:14> eps
> ans =  2.2204e-16
>
> so that eps*1e4 is about 1e-12.
> Why is it important to cast poles to 0?  Is there something later in  
> the code that depends on that?

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.

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

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.

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

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.

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.

Ben

p.s. I'll incorporate the changes and post an update to the patch.





-------------- next part --------------
An HTML attachment was scrubbed...
URL: https://www.cae.wisc.edu/pipermail/octave-maintainers/attachments/20071231/71394f8d/attachment-0001.html