[Octave-bug-tracker] [bug #34266] Four issues with "residue" and residuez"

[Dan Sebald]

Follow-up Comment #15, bug #34266 (project octave):

Below is some observations about this residue() routine.  I didn't do a deep
dive, but here's my summary:

1) The algorithm for computing 'small' of rresidue() doesn't seem as nuanced
as it should be.

2) The tolerance for this type of operation may be unnecessarily small because
there are singular matrices occurring along the way.  Do we want to do a
deeper dive on this?



>> The failure in this case doesn't seem to be precision related, at least
>> not directly.  The failure is for dimension mismatch:
>>
>> octave:23> assert (br, b, 1e-8);
>> error: ASSERT errors for:  assert (br,b,1e-8)
>>
>>    Location  |  Observed  |  Expected  |  Reason
>>       .          O(1x1)       E(1x2)      Dimensions don't match
>> octave:23> br
>> br =    7.0373e+06
>> octave:24> b
>> b =
>>
>>     1.0000e+00   7.0373e+06
>>
>> octave:25> assert (ar, a, 1e-8);
>> error: ASSERT errors for:  assert (ar,a,1e-8)
>>
>>    Location  |  Observed  |  Expected  |  Reason
>>       .          O(1x4)       E(1x5)      Dimensions don't match
>> octave:28> ar
>> ar =
>>
>>     3.6412e+09   1.5697e+18   0.0000e+00   0.0000e+00
>>
>> octave:29> a
>> a =
>>
>>     1.0000e+00   3.6412e+09   1.5697e+18   0.0000e+00   0.0000e+00
>>
>>
>> Is there an assumed coefficient value of 1 for the zeroeth order of the
>> polynomial (a common assumption) somewhere along the line that is lost?
>>   Or is that 1 being lost because of some tiny loss of precision
>> somewhere?
>>
>> I'd point out that the poles in this example are far into the left half
>> plane and not real near one another so my intuition is that numerical
>> issues shouldn't be problem.
>
> Actually, I see that there are two poles at zero placed into the input
> array, but that's not the issue.
>
> The problem with this example is that there is something inherently bad
> about the example.  There is this subroutine rresidue() inside residue.m
> that appears to have some kind of limitation with poles that far away
> from the real(s) = 0 line.  Given the values in this example for p1 and
> p2, the first format of residue gives a matrix warning (and I've printed
> out A and B from line 268):
>
> octave:153> [rrr ppp kkk eee] = residue(b,a);
> A =
>
>     0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00
>     3.6412e+09   0.0000e+00   3.1416e+09   4.9965e+08
>     1.5697e+18   3.6412e+09   0.0000e+00   5.1200e+02
>     0.0000e+00   1.5697e+18   0.0000e+00  -1.6085e+12
>
> B =
>
>     0.0000e+00
>     0.0000e+00
>     1.0000e+00
>     7.0373e+06
>
> warning: matrix singular to machine precision
> warning: called from
>      residue at line 268 column 5
>
> Something is causing that rresidue() subroutine to prepad a zero in each
> column--something that doesn't happen if I choose p1 and p2 to be more
> reasonable values.
>
> So, I'm wonder if the routine is being given data for which some
> operation (deconvolution?) is failing.  I've tried changing tolerance,
> but that seems to have no effect.  Should such a test even be run?

I can see where the loss of the 1 at the front of ar and br is coming from. 
In the rresidue sub-routine is the following test:

  ## Check for leading zeros and trim the polynomial coefficients.
  if (isa (r, "single") || isa (p, "single") || isa (k, "single"))
    small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps ("single");
  else
    small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps;
  endif

which is followed by polyreduce, a routine meant to strip leading "zeros". 
The point of this code must be to get rid of some type of pathological case, I
guess.  However, for this example, those distant poles lead to a significant
value for the variable 'small':

octave:175>  [br, ar] = residue (r, p, k, e);
small =  348.54
pnum =

  -9.6296e-35  -4.2894e-25   1.0000e+00   7.0373e+06

pden =

   1.0000e+00   3.6412e+09   1.5697e+18   0.0000e+00   0.0000e+00

Consequently, that 1.000 at the front of the denominator is getting stripped
when it shouldn't.  (The two leading pnum values though, should be stripped.) 
So, that formula for 'small' doesn't implement exactly what it is intended.

Now, if I force small to something that makes pden remain as is, THEN I get
the numerical problem:

error: ASSERT errors for:  assert (br,b,1e-8)

  Location  |  Observed  |  Expected  |  Reason
    (2)      7037297.6777 7037297.6777   Abs err 9.7789e-08 exceeds tol 1e-08



    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?34266>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/