bugGNU Octave - Bugs: bug #32502, cov function returns incorrect...

 
 

bug #32502: cov function returns incorrect results for complex matrices

Submitter:  Ivan <lisp_hacker>
Submitted:  Wed 16 Feb 2011 09:51:18 AM UTC
   
 
Category:  Libraries Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Wont Fix Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * 3.2.4 Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Thu 29 Sep 2011 05:14:49 AM UTC, comment #6: 

uff... mathematical notation looks ugly =)

I mean:

$$
c = \left (\frac{x^H x}{n-1} \right)^H
$$

Ivan <lisp_hacker>
Thu 29 Sep 2011 05:11:22 AM UTC, comment #5: 

Still  I'm not sure that Octave do the things right. Once again here is the code (from Octave):

c = conj (x' * x / (n - 1));

Here we see a double complex conjugation. First operator' do conjugation, and after it there is an explicit call to conj. In mathematical notation we do the following:

        H    H
    / X   X \  
C = |--------|
    \  N-1  /

So I think that the extra call to conj is erroneous, am I right?

Ivan <lisp_hacker>
Fri 18 Feb 2011 11:54:11 AM UTC, comment #4: 

I think Octave does the right thing here. See the discussion on Matlab Newsreader http://www.mathworks.com/matlabcentral/newsreader/view_thread/301757#815851 where this issue was discussed. The question is whether we should be compatible with Matlab or not.

Stefan

Anonymous
Wed 16 Feb 2011 07:35:33 PM UTC, comment #3: 

Sorry missed some lines in octave code. The correct example would be:

 A_i = load("-ascii", "imag_input.ascii");
 A_r = load("-ascii", "real_input.ascii");
 A = complex(A_r, A_i);
 imag(cov(A))


Ivan <lisp_hacker>
Wed 16 Feb 2011 07:19:41 PM UTC, comment #2: 

Here goes octave:

 A_r = load("-ascii", "real_input.ascii");
 A_i = load("-ascii", "imag_input.ascii");
 A = complex(A_r, A_i);
 imag(A)

ans =

  -0.00000  -0.25112   0.16022  -0.44749  -0.06819   0.23258  -0.05861  -0.06876
   0.25112  -0.00000  -0.07418  -0.05551   0.11349   0.01903   0.15368  -0.09032
  -0.16022   0.07418  -0.00000   0.28853  -0.07469  -0.16844   0.18365   0.02444
   0.44749   0.05551  -0.28853  -0.00000   0.12424   0.04581  -0.23836   0.10005
   0.06819  -0.11349   0.07469  -0.12424  -0.00000   0.03534  -0.20799   0.32262
  -0.23258  -0.01903   0.16844  -0.04581  -0.03534  -0.00000  -0.34949   0.11833
   0.05861  -0.15368  -0.18365   0.23836   0.20799   0.34949  -0.00000   0.15882
   0.06876   0.09032  -0.02444  -0.10005  -0.32262  -0.11833  -0.15882  -0.00000


Since matlab is prohibited (and actually I don't have access to it right now), a C++ code snippet, doing the same (it uses armadillo library):

#include <armadillo>
#include <iostream>

int main() {
    using namespace arma;
    mat real_input, imag_input;
    real_input.load("data/real_input.ascii", raw_ascii);
    imag_input.load("data/imag_input.ascii", raw_ascii);
    cx_mat input(real_input, imag_input);
    std::cout << imag(cov(input)) << "\n";
}

$ g++ test.cc -I../../../contrib -llapack -lblas
$ ./a.out
        0   0.2511  -0.1602   0.4475   0.0682  -0.2326   0.0586   0.0688
  -0.2511        0   0.0742   0.0555  -0.1135  -0.0190  -0.1537   0.0903
   0.1602  -0.0742        0  -0.2885   0.0747   0.1684  -0.1837  -0.0244
  -0.4475  -0.0555   0.2885        0  -0.1242  -0.0458   0.2384  -0.1000
  -0.0682   0.1135  -0.0747   0.1242        0  -0.0353   0.2080  -0.3226
   0.2326   0.0190  -0.1684   0.0458   0.0353        0   0.3495  -0.1183
  -0.0586   0.1537   0.1837  -0.2384  -0.2080  -0.3495        0  -0.1588
  -0.0688  -0.0903   0.0244   0.1000   0.3226   0.1183   0.1588        0

As you can see this matrix are complex conjugates of each other. Matlab yield the same result, btw.

(file #22704, file #22705)

Ivan <lisp_hacker>
Wed 16 Feb 2011 06:19:56 PM UTC, comment #1: 

Please do not quote code from Matlab here, not matter how small.

Please do give us an example of calling cov that produces an incorrect output, and show the output you are expecting instead.

John W. Eaton <jwe>
Group administrator
Wed 16 Feb 2011 09:51:18 AM UTC, original submission:  

For some reason function cov conjugates it's result:
Here is the code from octave 'm/statistics/base/cov.m':

   x = center (x, 1);
   c = conj (x' * x / (n - 1));


Matlab (and other toolkits) do not conjugate it's results
(Matlab code):

   xy = (xc' * xc) / (m-1);


Ivan <lisp_hacker>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #22704:  real_input.ascii added by lisp_hacker (2KiB - application/octet-stream)
file #22705:  imag_input.ascii added by lisp_hacker (2KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by rik5 (Updated the item)
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by lisp_hacker (Submitted the item)
  • -email is unavailable- added by lisp_hacker
  •  

    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 6 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2011-09-28 rik5 StatusNeed Info Wont Fix
        Open/ClosedOpen Closed
    2011-02-16 lisp_hacker Attached File- Added real_input.ascii, #22704
        Attached File- Added imag_input.ascii, #22705
    2011-02-16 jwe StatusNone Need Info
    2011-02-16 lisp_hacker Carbon-Copy- Added lisp_hacker

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code