bugGNU Octave - Bugs: bug #53897, poly of conjugate complex pairs...

 
 

bug #53897: poly of conjugate complex pairs should return a real polynomial

Submitter:  Marco Caliari <caliari>
Submitted:  Mon 14 May 2018 09:11:56 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Unexpected Error or Warning
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 06 Mar 2019 06:48:46 PM UTC, comment #17: 

I checked in the patch here (https://hg.savannah.gnu.org/hgweb/octave/rev/a018fca707ca).

I made a few small modifications.  In place of tmp1, tmp2, and tmp I used what I hope are more meaningful variable names pos_imag, neg_imag, and is_equal.

In place of using all() and size() together I just used the C++ function size_equal().


-    if (all (size (tmp1) == size (tmp2)))
+    if (size_equal (pos_imag, neg_imag))


For the BIST test, I used commas to clearly delineate where each number stops which can sometimes be ambiguous.  Is "-9 - 3i" a single complex number or is it two numbers, "-9+0i AND 0-3i".


-%! assert (y, [1 +  0i   -9 -  3i   25 + 24i  -17 - 57i  -12 + 36i])
+%! assert (y, [1 + 0i, -9 - 3i, 25 + 24i, -17 - 57i, -12 + 36i]);


Thanks Marco for coming up with a high-performance patch (no isequal call).

This will be a part of the bug fix release 5.2.0 in a few months.

Rik <rik5>
Group administrator
Wed 06 Mar 2019 05:57:27 PM UTC, comment #16: 

Right, got it. Patch works for me and includes test case, looks good.

Mike Miller <mtmiller>
Group Member
Wed 06 Mar 2019 02:38:16 PM UTC, comment #15: 

There is probably a bug in cheby2, ok. But poly(p1) is valid syntax and should work. The problem is the use of == instead of the slower isequal. I attach a patch without using isequal.

(file #46448)

Marco Caliari <caliari>
Group Member
Wed 06 Mar 2019 01:37:10 AM UTC, comment #14: 

Robert - I think you may want to have a look at bug #49996 and the patch posted there. This has been a long-standing bug in cheby2, it's just exposed further by this fix.

Mike Miller <mtmiller>
Group Member
Wed 06 Mar 2019 01:17:06 AM UTC, comment #13: 

This patch fails for zp2tf from the signal package with values that are "nearly" real. Copying from line 84 of share/octave/5.1.0/m/polynomial/poly.m:

octave:14> [z1,p1,k1]=cheby2(5,40,0.1);
octave:15> p1
p1 =
   9.345539872652292e-01 - 1.826618879469199e-01i
   8.372613608092301e-01 - 1.304057652412007e-01i
   7.781404349401341e-01 - 1.537554717468438e-17i
   8.372613608092301e-01 + 1.304057652412007e-01i
   9.345539872652292e-01 + 1.826618879469199e-01i

octave:16> tmp = sort (p1(imag (p1) > 0)) == sort (conj (p1(imag (p1) < 0)))
error: mx_el_eq: nonconformant arguments (op1 is 2x1, op2 is 3x1)
octave:16>


Robert Jenssen <morgawr>
Fri 26 Oct 2018 02:46:06 AM UTC, comment #12: 

isequal is relatively slow, about 2X slower than '=='.  For a 100-element vector the absolute difference is still small, about 400 microseconds, but I think it is worth testing using isempty before using all.  I made that change and another small tweak or two and checked in the changeset here (https://hg.savannah.gnu.org/hgweb/octave/rev/e02f500a0980).

Marking as fixed and closing report.


Rik <rik5>
Group administrator
Fri 19 Oct 2018 04:48:48 PM UTC, comment #11: 

@Markus: thanks.
A patch is included.

(file #45248)

Marco Caliari <caliari>
Group Member
Fri 19 Oct 2018 04:25:35 PM UTC, comment #10: 

Would this do what you need?

x = [1+1i, 2];
isequal (sort (x(imag (x) > 0)), sort (conj (x(imag (x) < 0))))


Markus Mützel <mmuetzel>
Group administrator
Fri 19 Oct 2018 11:40:42 AM UTC, comment #9: 

If you take


x = [1+1i,2];
all(sort(x(imag(x)>0))==sort(conj(x(imag(x)<0))))


you get true, which is not correct. How do I check whether sort(x(imag(x)>0)) is the same vector of sort(conj(x(imag(x)<0))))?

Marco Caliari <caliari>
Group Member
Thu 18 Oct 2018 08:06:51 AM UTC, comment #8: 

I will do asap.

Marco Caliari <caliari>
Group Member
Wed 17 Oct 2018 06:49:20 PM UTC, comment #7: 

Okay, then I think this is fine to work up in to a true changeset and commit it.

Rik <rik5>
Group administrator
Wed 17 Oct 2018 06:37:57 PM UTC, comment #6: 

As far as I know, eig (that is, the undelying LAPACK routines) gives exact complex conjugates for real matrices. This is confirmed by experiments like


octave:1> A = randn (100);
octave:2> x = eig (A);
octave:3> all(sort(x(imag(x)>0))==sort(conj(x(imag(x)<0))))
ans = 1


Therefore, I would not introduce the tolerance.

Marco Caliari <caliari>
Group Member
Wed 17 Oct 2018 04:49:59 AM UTC, comment #5: 

That seems to work.  Would the call to eig ever return values which are not exactly complex conjugates such that we need to use a tolerance on the equality test?  In other words,


if (all (abs (sort (x(imag(x)>0)) - sort (conj (x(imag(x)<0)))) < tol))



Rik <rik5>
Group administrator
Tue 16 Oct 2018 07:00:30 AM UTC, comment #4: 

What about the following?


if all(sort(x(imag(x)>0))==sort(conj(x(imag(x)<0))))
  y = real (y);
endif


Marco Caliari <caliari>
Group Member
Mon 15 Oct 2018 04:19:17 PM UTC, comment #3: 

Yes, but is it on the right track?  Using your example from comment #2 I find that the largest imaginary value is 10 eps which is quite small.  At least there are no horrendous errors like 256i as in the first example.

It seems that sort followed by a test to replace small values with 0 would work.  The algorithm in poly.m is


  for j = 1:n;
    y(2:(j+1)) = y(2:(j+1)) - v(j) .* y(1:j);
  endfor


which has a multiplication and a subtraction per round of the algorithm.  To me that says worst case of 2 eps per cycle or 2 x n x eps for the algorithm.  It's not perfect, but for eps one could use the largest value in v.


epsmax = 2* n * eps (max (abs (v(:))));
idx = abs (imag (y)) < epsmax;
y(idx) = real (y(idx));


Of course, this is starting to look like a fair amount of work so if detecting the conjugate pairs is easier then we could do that as you suggested in comment #0.

Rik <rik5>
Group administrator
Mon 15 Oct 2018 06:29:23 AM UTC, comment #2: 

@Rik: your idea does not work for the following example


x=[1,sqrt(2)/2+sqrt(2)/2*i,1i,-sqrt(2)/2+sqrt(2)/2*i,-1,-sqrt(2)/2-sqrt(2)/2*i,-1i,sqrt(2)/2-sqrt(2)/2*i]


i.e., elements of the vector with the same absolute value.

Marco Caliari <caliari>
Group Member
Sat 13 Oct 2018 04:37:00 PM UTC, comment #1: 

The algorithm in poly.m uses a for loop and repeated multiplications.  It could be something as simple as, with the inversion of matrices, using balance() to put the largest factors first, etc.

In this example, if I use sort on the input x then the conjugates are paired and the result of poly is entirely real without spurious imaginary terms.

See the session below


octave:25> n = 3;
octave:26> fr = 436;
octave:27> x = 1i * 2 * pi * fr * [-n:-1 1:n]
x =

     -0.0 - 8218.4i     -0.0 - 5478.9i     -0.0 - 2739.5i      0.0 + 2739.5i      0.0 + 5478.9i      0.0 + 8218.4i

octave:28> N = poly (x)
N =

 Columns 1 through 6:

   1.0000e+00 + 0.0000e+00i   0.0000e+00 - 1.8190e-12i   1.0507e+08 + 0.0000e+00i   0.0000e+00 + 1.8311e-04i   2.7597e+15 + 0.0000e+00i   0.0000e+00 + 2.5600e+02i

 Column 7:

   1.5216e+22 + 0.0000e+00i

octave:29> x2 = sort (x)
x2 =

     -0.0 - 2739.5i      0.0 + 2739.5i     -0.0 - 5478.9i      0.0 + 5478.9i     -0.0 - 8218.4i      0.0 + 8218.4i

octave:30> N2 = poly (x2)
N2 =

   1.0000e+00   0.0000e+00   1.0507e+08   0.0000e+00   2.7597e+15   0.0000e+00   1.5216e+22

octave:31> roots (N2)
ans =

     -0.0 + 8218.4i
     -0.0 - 8218.4i
     -0.0 + 5478.9i
     -0.0 - 5478.9i
     -0.0 + 2739.5i
     -0.0 - 2739.5i


Should we just sort the input vector 'v' after the input validation in poly.m?


Rik <rik5>
Group administrator
Mon 14 May 2018 09:11:56 AM UTC, original submission:  

If you try the following


n = 3;
fr = 436;
x = 1i * 2 * pi * fr * [-n:-1 1:n];
N = poly (x);


you see that N has spurious imaginary parts (quite large, in absolute value). This should not happen, since the roots x are complex conjugate pairs. We should first detect this fact (I thought to imag (sum (sort (x))) == 0, it is probably possible to do without computations, but only with sorting and comparison). Then, either we use the formula

+verbatim
(y-x(1))*(y-conj(x(1)) = y^2-2*real(x(1))*y+abs(x(1))^2
-verbatim-

for each complex conjugate pair or we simply remove the spurious imaginary parts after the computation of the coefficients.

Marco Caliari <caliari>
Group Member

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #46448:  poly2.diff added by caliari (1KiB - application/x-tex)
file #45248:  poly.diff added by caliari (1KiB - application/x-tex)

 

Depends on the following items: None found

Digest:
   bug dependencies.

 

Carbon-Copy List
  • -email is unavailable- added by morgawr (Posted a comment)
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by caliari (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 14 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2019-03-08 mtmiller Dependencies- bugs #55862 is dependent
    2019-03-06 rik5 StatusPatch Reviewed Fixed
        Open/ClosedOpen Closed
    2019-03-06 mtmiller StatusPatch Submitted Patch Reviewed
    2019-03-06 caliari Attached File- Added poly2.diff, #46448
        Item GroupInaccurate Result Unexpected Error or Warning
        StatusFixed Patch Submitted
        Open/ClosedClosed Open
    2018-10-26 rik5 StatusPatch Submitted Fixed
        Open/ClosedOpen Closed
    2018-10-19 caliari Attached File- Added poly.diff, #45248
        StatusIn Progress Patch Submitted
    2018-10-18 caliari StatusNeed Info In Progress
    2018-10-13 rik5 StatusNone Need Info

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code