bugGNU Octave - Bugs: bug #42742, polygcd fails valid test


bug #42742: polygcd fails valid test

Submitted by:  Michael Godfrey <godfrey>
Submitted on:  Sat 12 Jul 2014 01:04:56 PM UTC  
Category: Octave FunctionSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Incorrect Result
Status: NoneAssigned to: None
Originator Name: GodfreyOpen/Closed: Open
Release: devOperating System: Any

Add a New Comment (Rich MarkupRich Markup):

You are not logged in

Please log in, so followups can be emailed to you.


Wed 13 Aug 2014 01:41:52 AM UTC, comment #5:

The following works for me in that it enlarges the tolerance based on the input and all tests pass even when run repeatedly.

There is an intuitive logic to doing this, but I haven't been able to convince myself rigorously why it should be so.

Rik <rik5>
Project Administrator
Wed 06 Aug 2014 10:52:29 PM UTC, comment #4:

Reiterating Rik's comment, I guess that is pretty much the same as what I pointed out. The

is the issue. The vector r can have a very wide range in terms of magnitude. Yes, perhaps applying norm in some way would do it.

An alternative test for the roots()-based approach might be to incorporate the magnitude of the root. That is, in the patch I submitted previously is this test:

which tests if the roots of the different polynomials are the same. That could be changed to something like:

I put the '1' in there to prevent decreasing the tolerance to something potentially smaller than eps. This works, but again, the tests fail because of something like the following:

The solution is good, but because the poly() routine multiplies all those roots, the many small differences in the roots all accumulate in a multiplicative way. So, again, the tolerance of the results may need refinement, but also in some sense I wish the roots function produced some more accurate results too.

(file #31857)

Dan Sebald <sebald>
Wed 06 Aug 2014 09:24:49 PM UTC, comment #3:

I've attached a patch that implements polygcd() in the way I described. That is, uses the roots() function and finds common roots. However, the first impression is that it too has the same sorts of tolerance problems.

I increased the polygcd tolerance (i.e., the input parameter) and then get the results expected, but the resulting polynomial coefficients then don't pass the tolerance test, off by some factor slightly greater than 10*sqrt(eps).

I'm not sure about this numerical stability comment in the documentation for plygcd(). It seems to me the bigger issue is that applying the tolerance in the way it is done is the issue. Taking a look at some of the tests in polygcd(), it seems to me that in one case a coefficient like 1 should be within 10sqrt(eps), but then also a coefficient like 20,000 should also be accurate with 10sqrt(eps). It's easily possible with polynomial orders of say 8 to have numbers in the range of 1 to 2e4. There's nothing that takes the magnitude of numbers into account. One would prefer a relative test, say 0.001%, rather than absolute in this case.

(file #31856)

Dan Sebald <sebald>
Wed 06 Aug 2014 08:44:53 AM UTC, comment #2:

In the polygcd() routine is this line:

Are the zeros of r always assured to be at the front of the array? Perhaps so.

Anyway, is there another way of solving this instead of using deconvolution? The documentation states:

"This is equivalent to the polynomial found by multiplying together all the common roots."

It is pretty straightforward to compute the roots of both of the polynomials using the roots() function. Then if one were to put the roots in order there could be a fast method of finding common roots (as opposed to a factorial permutation comparison of all combinations of roots). Here's an illustration.

All that is needed is the short simple loop inside a loop that picks the common roots within some tolerance. (There is the intersect() command, but that might require equality, I'm not sure.)

The question is whether the proposed approach is any more robust in the case of polynomial coefficients being large.

Dan Sebald <sebald>
Tue 05 Aug 2014 08:36:48 PM UTC, comment #1:

The problem can be partially alleviated by using a tolerance. The default is sqrt (eps) but I find that if I use 10 times that value the frequency of an error drops from an average of 16/1000 iterations to a mean of 1/1000.

Test code I used was:

One possibility is that the tolerance, which defaults to sqrt (eps), should be based on some characteristic of the input vectors to polygcd such as the norm.

The relevant lines in polygcd are

In cases where it fails the remainder vector r is very nearly zero, i.e., the remainder polynomial is very small. But there is one coefficient which just manages to exceed the tolerance so the algorithm continues for one more trip through the while loop and gets the wrong answer.

A sample failing set of variables p, p1, p2, obs, exp is attached to the report as polybad.var for those who want a test case.

(file #31842)

Rik <rik5>
Project Administrator
Sat 12 Jul 2014 01:04:56 PM UTC, original submission:

This was first reported on the Mainbtainters List, but it appears
to be a bug in polygcd.

On 07/07/2014 06:57 PM, Rik wrote:

> All,
> I'm getting occasional failures from the test code in polygcd. It seems to
> happen about 10% of the time according to the following code:
> for i = 1:100
> bm(i) = test ("polygcd");
> endfor
> sum (bm)
> The test that fails is
> %!test
> %! for ii=1:10
> %! p = (unique (randn (10, 1)) * 10).';
> %! p1 = p(3:end);
> %! p2 = p(1:end-2);
> %! assert (polygcd (poly (-p1), poly (-p2)), poly (- intersect (p1, p2)),
> sqrt (eps));
> %! endfor
> I tried a few different random seeds to see if I could fix the value to
> something that would always pass, but no luck.
> The simplest thing is to make this an %!xtest which can occasionally fail.
> But if someone understands polygcd and could suggest a way to modify the
> test that would be preferable.
> The error I get is a dimensional mismatch:
> !!!!! test failed
> ASSERT errors for: assert (polygcd (poly (-p1), poly (-p2)),poly
> (-intersect (p1, p2)),sqrt (eps))
> Location | Observed | Expected | Reason
> . O(1x1) E(1x7) Dimensions don't match
> --Rik

Here is where this goes wrong:

If the line [d, r] = deconv (b, a); in polygcd returns a r vector whose first element
is zero, then the line a = r / r(1); in polygcd fails. At this point a ends up
as a = 1 and x = 1, and the assert fails as above. I do not right now see why a = 1
due to divide by zero, but this is the cause of the failure.

The following script will, for me, reliably produce the failure as described above:

for ii=1:1000
p = (unique (randn (10, 1)) * 10).';
p1 = p(3:end);
p2 = p(1:end-2);
assert (polygcd (poly (-p1), poly (-p2)), poly (- intersect (p1, p2)), sqrt (eps));

Obviously, the line a = r / r(1); should not be executed if r(1) = 0. But, it is
not clear to me right now what should be done instead.

I hope this helps.

Michael Godfrey <godfrey>


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

Attach File(s):

Attached Files
file #31857:  octave-polygcd_2014aug06_2.patch added by sebald (1kB - application/octet-stream)
file #31856:  octave-polygcd_2014aug06.patch added by sebald (1kB - application/octet-stream)
file #31842:  polybad.var added by rik5 (663B - application/octet-stream)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -unavailable- added by sebald (Posted a comment)
  • -unavailable- added by rik5 (Updated the item)
  • -unavailable- added by godfrey (Submitted the item)

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.


    Please enter the title of George Orwell's famous dystopian book (it's a date):



    Follow 3 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Wed 06 Aug 2014 10:52:29 PM UTCsebaldAttached File-=>Added octave-polygcd_2014aug06_2.patch, #31857
    Wed 06 Aug 2014 09:24:49 PM UTCsebaldAttached File-=>Added octave-polygcd_2014aug06.patch, #31856
    Tue 05 Aug 2014 08:36:48 PM UTCrik5Attached File-=>Added polybad.var, #31842

    Back to the top

    Powered by Savane 3.1-cleanup