bugGNU Octave - Bugs: bug #49291, residue function gives...

 
 

bug #49291: residue function gives incompatible result

Submitter:  None
Submitted:  Fri 07 Oct 2016 12:12:48 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Matlab Compatibility
Status:  Fixed Assigned to:  None
Originator Name:  Frank Wiedmann Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * 4.0.3
Operating System:  * Any Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Mon 10 Oct 2016 08:03:31 PM UTC, comment #11: 

I deleted the code that was removing small coefficients.  See http://hg.savannah.gnu.org/hgweb/octave/rev/dc872d5d74c5.  This also may have fixed bug #34266.  The fix was made on the development branch as we are very close to the 4.2.0 release.

Rik <rik5>
Group administrator
Mon 10 Oct 2016 06:46:06 PM UTC, comment #10: 

br has size 1x4, with br(1) around -9e-35 and br(2) around -4e-25.

Marco Caliari <caliari>
Group Member
Mon 10 Oct 2016 06:16:52 PM UTC, comment #9: 

@Marco: Thanks for testing.  Can you try another one from the BIST tests?  This is from https://savannah.gnu.org/bugs/?34266.


 z1 =  7.0372976777e6;
 p1 = -3.1415926536e9;
 p2 = -4.9964813512e8;
 r1 = -(1 + z1/p1)/(1 - p1/p2)/p2/p1;
 r2 = -(1 + z1/p2)/(1 - p2/p1)/p2/p1;
 r3 = (1 + (p2 + p1)/p2/p1*z1)/p2/p1;
 r4 = z1/p2/p1;
 r = [r1; r2; r3; r4];
 p = [p1; p2; 0; 0];
 k = [];
 e = [1; 1; 1; 2];
 b = [1, z1];
 a = [1, -(p1 + p2), p1*p2, 0, 0];
 [br, ar] = residue (r, p, k)


If we eliminate all scrubbing code for small values then Octave returns a value for br which is 1x4 where the first two entries are nearly zero.  If Matlab does the same thing then we can be reasonably certain they are not applying the filtering present in polyreduce.

There is still a possibility that they are replacing values less than eps(1) with 0.  But first I want to see if Matlab returns a 1x4 vector for br.


Rik <rik5>
Group administrator
Mon 10 Oct 2016 05:43:50 PM UTC, comment #8: 

Uhm... matlab's residue does not accept 4 outputs. All the asserts (except that on e) are correct, length(br)=4 with br(1) exactly 0.

Marco Caliari <caliari>
Group Member
Mon 10 Oct 2016 04:41:04 PM UTC, comment #7: 

@Marco: Could you try the following code in Matlab?  It is from one of the BIST tests for residue.  If Matlab is really not applying any filtering for small results than br will be a 4-element vector where the first entry is nearly 0.


b = [1, 0, 1];
a = [1, 0, 18, 0, 81];
[r, p, k, e] = residue (b, a);
r1 = [-5i; 12; +5i; 12]/54;
p1 = [+3i; +3i; -3i; -3i];
assert (r, r1, 1e-12);
assert (p, p1, 1e-12);
assert (isempty (k));
assert (e, [1; 2; 1; 2]);
[br, ar] = residue (r, p, k);
assert (br, b, 1e-12);
assert (ar, a, 1e-12);



Rik <rik5>
Group administrator
Mon 10 Oct 2016 11:52:47 AM UTC, comment #6: 

To comment #1

I think the code which tries to remove small coefficients should be removed, since too dangerous. I tried the following in matlab


a=10*[1, -5, 8, -4];
b=[1,0,1];
[r,p,k]=residue(b,a);
[b,a]=residue(r,p,k);


b(2) is about 6.66e-16 and less than max([abs(a),abs(b)])*eps. But it was not converted to 0 by matlab.

Marco Caliari <caliari>
Group Member
Fri 07 Oct 2016 10:15:18 PM UTC, comment #5: 

Opps. Never mind. I just checked. The no common roots.

Ben Abbott <bpabbott>
Group Member
Fri 07 Oct 2016 10:11:36 PM UTC, comment #4: 

@RIk,

The polynomials "a" and "b" are each 3rd order. I haven't checked, but perhaps the r, p, & k you proposed produced a numerator which shares a root with the denominator?

Ben Abbott <bpabbott>
Group Member
Fri 07 Oct 2016 07:26:00 PM UTC, comment #3: 

More thoughts, since residue is invertible it should be able to apply residue twice and arrive back at the original inputs.  This works for more ordinary inputs, but not for the special case presented here.  I made a script, tst_res.m, which is attached.


rf = [1e3 2e3 1e3 2e3];
cf = [316.2e-9 50e-9 31.6e-9 5e-9];
r = 1./cf
p = -1./(rf.*cf)
k = 0
[b, a] = residue (r, p, k);
[r2, p2, k2] = residue (b, a)


When running it in Octave I get,


tst_res
r =

   3.1626e+06   2.0000e+07   3.1646e+07   2.0000e+08

p =

  -3.1626e+03  -1.0000e+04  -3.1646e+04  -1.0000e+05

k = 0
r2 =

   3.1703e+06
   2.3710e+07
   2.1554e+07

p2 =

  -3.1573e+03
  -1.1024e+04
  -1.9857e+04

k2 =  1759.6


As you can see, r2 and p2 don't even have the correct number of elements, and k2 is wildly off from 0.


(file #38683)

Rik <rik5>
Group administrator
Fri 07 Oct 2016 04:49:42 PM UTC, comment #2: 

As a second example,


b = [2 1 0 0];
a = [1 0 1 1];
[r,p,k] = residue(b,a);
[b2, a2] = residue (r,p,k)
b2 =

 Columns 1 through 3:

   2.00000 + 0.00000i   1.00000 + 0.00000i   0.00000 - 0.00000i

 Column 4:

   0.00000 + 0.00000i

a2 =

   1.00000   0.00000   1.00000   1.00000

imag (b2)

ans =

   0.0000e+00   4.9865e-17  -3.4024e-17   7.3080e-17


As you can see, b2 is nearly equivalent to the starting variable b.  However, if the code to trim polynomial coefficients was changed to act separately on the real part and on the imaginary part then these small numbers would be removed and there would be complete equality.  Perhaps that is an addition to the code that is worth considering separately from the first issue reported.

Rik <rik5>
Group administrator
Fri 07 Oct 2016 04:34:35 PM UTC, comment #1: 

This seems to be entirely due to a cleanup routine in residue.m.

I used the debugger to step through the code while executing your example, and I have the Matlab-equivalent result until the end of the file.


pden =

   1.0000e+00   1.4481e+05   4.9290e+09   4.5817e+13   1.0008e+17


At that point, the m-file tries to remove elements which it thinks are likely due to round-off errors in the calculation.  In this case, anything less than the machine precision of the largest value.  The code starting at line 341 is


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

  pnum(abs (pnum) < small) = 0;
  pden(abs (pden) < small) = 0;

  pnum = polyreduce (pnum);
  pden = polyreduce (pden);

endfunction


The obvious question is whether this code is necessary.  If it is necessary, should it apply both to the numerator (pnum) and the denominator (pden) or just the numerator?  And if it is still necessary, and it should apply to both pnum and pden, then how should the code be changed to recover the original (and Matlab) result.

This code is really old and some form of this truncation has been in place since 2007.  I'm adding the authors listed in the file to this bug report in case they remember why this should be the way it is.

The Matlab documentation (http://www.mathworks.com/help/matlab/ref/residue.html) says that the denominator polynomial should be the size of 'r' or 'p' which in this case is 4.  So it is possible that this is also a bug/inconsistency in Matlab.

Rik <rik5>
Group administrator
Fri 07 Oct 2016 12:12:48 PM UTC, original submission:  

The Octave result


>> clear all
>> rf = [1e3 2e3 1e3 2e3];
>> cf = [316.2e-9 50e-9 31.6e-9 5e-9];
>> [num,den] = residue(1./cf,-1./(rf.*cf),0)
num =

  2.5481e+008  1.5687e+013  2.1808e+017  6.0049e+020

den =

  1.4481e+005  4.9290e+009  4.5817e+013  1.0008e+017

>>

seems to be incorrect (the first element of den is missing). With MATLAB, I get what seems to be the correct result:


>> format shorte
>> clear all
>> rf = [1e3 2e3 1e3 2e3];
>> cf = [316.2e-9 50e-9 31.6e-9 5e-9];
>> [num,den] = residue(1./cf,-1./(rf.*cf),0)

num =

   2.5481e+08   1.5687e+13   2.1808e+17   6.0049e+20


den =

   1.0000e+00   1.4481e+05   4.9290e+09   4.5817e+13   1.0008e+17

>>


Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #38683:  tst_res.m added by rik5 (152B - d2l/unknowntype)

 

Digest:
   bug dependencies.

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by bpabbott (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by rik5
  • -email is unavailable- added by rik5
  • -email is unavailable- added by None (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 10 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2016-10-10 rik5 StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2016-10-10 rik5 Operating SystemMicrosoft Windows Any
    2016-10-10 rik5 Dependencies- Depends on bugs #34266
    2016-10-07 rik5 Attached File- Added tst_res.m, #38683
    2016-10-07 rik5 Item GroupNone Matlab Compatibility
        StatusNone Confirmed
        Summaryresidue function gives incorrect result residue function gives incompatible result
    2016-10-07 rik5 Carbon-Copy- Added bpabbott
        Carbon-Copy- Added jwe

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code