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

 
 

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

Submitted by:  Marco Caliari <caliari>
Submitted on:  Mon 14 May 2018 09:11:56 AM UTC  
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Patch Submitted Assigned to:  None
Originator Name:  Open/Closed:  Open
Release:  dev Operating System:  GNU/Linux

Add a New Comment(Rich Markup)
   

You are not logged in

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

 

( Jump to the original submission)

Fri 19 Oct 2018 04:48:48 PM UTC, comment #11:

@Markus: thanks.
A patch is included.

(file #45248)

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

Would this do what you need?

Markus M├╝tzel <mmuetzel>
Project Member
Fri 19 Oct 2018 11:40:42 AM UTC, comment #9:

If you take

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>
Project Member
Thu 18 Oct 2018 08:06:51 AM UTC, comment #8:

I will do asap.

Marco Caliari <caliari>
Project 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>
Project 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

Therefore, I would not introduce the tolerance.

Marco Caliari <caliari>
Project 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,

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

What about the following?

Marco Caliari <caliari>
Project 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

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.

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>
Project Administrator
Mon 15 Oct 2018 06:29:23 AM UTC, comment #2:

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

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

Marco Caliari <caliari>
Project 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

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

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

If you try the following

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>
Project Member

 

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

Attach Files:
   
   
Comment:
   

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

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -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)
  •  

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

    Only project members can vote.

     

     

     

    Follow 4 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    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.3