bugGNU Octave - Bugs: bug #53869, Residue returning an inaccurate...


bug #53869: Residue returning an inaccurate result

Submitter:  None
Submitted:  Wed 09 May 2018 03:09:21 PM UTC
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Matlab Compatibility
Status:  Fixed Assigned to:  None
Originator Name:  Tony Richardson Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * 4.4.0
Operating System:  * Any Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup

Jump to the original submission

Tue 15 May 2018 06:51:59 PM UTC, comment #9: 

I pushed the change for the pre-conditioner to the stable branch here (http://hg.savannah.gnu.org/hgweb/octave/rev/8e30ef611002).

It really looks like there should be some additional code to either call polyreduce on the output, or possibly use a tolerance to eliminate leading values that are very small, like 1e-33.

Rik <rik5>
Group administrator
Tue 15 May 2018 06:28:28 AM UTC, comment #8: 

@Rik: Matlab produces a 4-element vector, with a leading 1.3e-33. Our result is now better, since in this test we should recover exactly the original b. Of course, we cannot guarantee this always happens.

Marco Caliari <caliari>
Group Member
Mon 14 May 2018 11:43:11 PM UTC, comment #7: 

I think this is the correct thing to do.  However, when I make the change and then run "test residue" I get a failure in the BIST tests.

octave:14> test residue
***** test
 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, [0,b], 1e-12);
 assert (ar, a, 1e-12);
!!!!! test failed
ASSERT errors for:  assert (br,[0, b],1e-12)

  Location  |  Observed  |  Expected  |  Reason
     .          O(1x3)       E(1x4)      Dimensions don't match

The issue is that the call

 [br, ar] = residue (r, p, k);

used produce a 4-element vector

[ 0, 1, 0, 1]

Now it produces just a 3-element vector

[ 1, 0, 1]

Given that this vector represents a polynomial, it seems better that the leading 0 should be removed.

@Tony: Could you check that Matlab produces a 3-element vector for the code above?

Rik <rik5>
Group administrator
Sat 12 May 2018 07:37:56 AM UTC, comment #6: 

What about a left preconditioning like

## Solve for the residues.
d = max (abs (A), [], 2);
r = (diag (d) \ A) \ (B ./ d);

It seems to fix the problem, for me. Of course, it would be better to construct A and B already in a normalized way, but I do not know the details of the algorithm.

Marco Caliari <caliari>
Group Member
Thu 10 May 2018 12:59:33 PM UTC, comment #5: 

In an interesting twist of fate, it looks as if I was one of the original authors of residue.m around 25 years ago (something that I had forgotten).  I believe I used a home brew algorithm based on partial fraction expansion.  It appears to have evolved significantly (thanks Ben Abbott) since then, but apparently is still not as accurate as whatever MATLAB is doing.

Tony Richardson <amrichardson>
Wed 09 May 2018 07:10:19 PM UTC, comment #4: 

Confirmed.  Calling residue_bug with n = 1,2 works fine.  However, with n = 3 there is already a warning printed and the condition number of the matrix is scarily high.

octave:6> debug_on_warning (1)
octave:7> residue_bug (3)
warning: matrix singular to machine precision, rcond = 5.73398e-20
warning: called from
    residue at line 265 column 5
    residue_bug at line 6 column 11
stopped in /home/rik/wip/Projects_Mine/octave-dev/scripts/polynomial/residue.m at line 265
265:   r = A \ B;

The residue algorithm is coded as an m-file so it should be easy to see why A is nearly singlular.  However, I don't understand the algorithm well enough to say how it should be modified or replaced.

Rik <rik5>
Group administrator
Wed 09 May 2018 06:47:26 PM UTC, comment #3: 

Sorry.  Trying again.

function residue_bug(n)
    alpha = 100;
    fr = 436;
    N=poly(1i*2*pi*fr*[-n:-1 1:n]);
    D=poly(-alpha+1i*2*pi*fr*[-n:-1 1:n]);
    [R, P, K]=residue(N,D)

Tony Richardson <amrichardson>
Wed 09 May 2018 06:45:13 PM UTC, comment #2: 

Can't figure out how to attach file to comment.  Trying verbatim method.

function residue_bug(n)
  alpha = 100;
  fr = 436;
  N=poly(1i2pifr[-n:-1 1:n]);
  D=poly(-alpha+1i2pifr[-n:-1 1:n]);
  [R, P, K]=residue(N,D)

Tony Richardson <amrichardson>
Wed 09 May 2018 06:33:49 PM UTC, comment #1: 

Could you upload the m-file residue_bug.m?  The bug tracker has interpreted your code using markup conventions where an asterisk means bold text.  The other possibility is to use the verbatim tag to mark a block of text that should not be interpreted.

+ verbatim +
- verbatim -

Remove the spaces between the '+' or '-' and the "verbatim" keyword to have it take effect.

Rik <rik5>
Group administrator
Wed 09 May 2018 03:09:21 PM UTC, original submission:  

Octave and MATLAB return different results for the following script with argument (n) values of 3 or larger:

function residue_bug(n)
    alpha = 100;
    fr = 436;
    N=poly(1i*2*pi*fr*[-n:-1 1:n]);
    D=poly(-alpha+1i*2*pi*fr*[-n:-1 1:n]);
    [R, P, K]=residue(N,D)

Octave displays a warning about a singular matrix and returns results that are different than MATLAB. MATLAB does not display a warning until n is much larger (around 30 or so).  I can normalize the problem so that Octave results are valid for values of n up to about 10, but MATLAB results are then valid for much larger values of n.

I thought that the MATLAB values might be incorrect (even though it does not display a warning), but I computed my ultimate result using a method that does not use residue and it is consistent with the result obtain using MATLAB's residue function.

The basic issue appears to be that whatever method MATLAB uses to compute residue results is much better behaved than the Octave method.



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

Attach Files:

No files currently attached


Depends on the following items: None found

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 amrichardson (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -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 4 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-05-15 rik5 StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2018-05-09 rik5 StatusNeed Info Confirmed
    2018-05-09 rik5 StatusNone Need Info

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code