bugGNU Octave - Bugs: bug #50377, [octave forge] (statistics) pdist2...

 
 

bug #50377: [octave forge] (statistics) pdist2 function returns complex results

Submitter:  Georg Wiora <gwiora>
Submitted:  Wed 22 Feb 2017 08:09:21 AM UTC
   
 
Category:  Octave Package Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Fixed Assigned to:  None
Originator Name:  gwiora Open/Closed:  * Closed
Release:  * 4.2.0 Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Thu 25 Nov 2021 11:02:10 AM UTC, comment #9: 

Status should be changed to closed?

Georg Wiora <gwiora>
Wed 24 Nov 2021 11:51:23 PM UTC, comment #8: 
Nir Krakauer <nir_krakauer>
Mon 07 Jan 2019 02:25:29 PM UTC, comment #7: 

I was asking as the file should probally have a test added to prove the fix

John Donoghue <lostbard>
Group Member
Mon 07 Jan 2019 02:05:09 PM UTC, comment #6: 

Sorry, it is so long ago I posted this sample that I forgot where I found the problem. But the problem is a standard case in numeric. Therefore I recommend the implementation of the suggested prevention for complex results.

An expression like sqrt(x^2+y^2 - 2 x y) can always become complex if the value of x^2+y^2 and 2 x y are equal in a range of eps. In this case the argument of sqrt() may become negative and the result type changed from double to complex.

This description may be used as a guide to constructing examples. Any numeric pair of [x,y] is a good test case where the condition x^2+y^2 == 2 x y analytically is true and x~=0 and y~=0 and I would suggest x~=y.


Georg Wiora <gwiora>
Mon 07 Jan 2019 01:22:12 PM UTC, comment #5: 

Do have any examples that show the issue in the unmodified code ?

John Donoghue <lostbard>
Group Member
Mon 07 Jan 2019 07:33:34 AM UTC, comment #4: 

Anyone willing to assign a status to this bug?

Georg Wiora <gwiora>
Wed 01 Mar 2017 03:13:13 PM UTC, comment #3: 

I reviewed your "max()" suggestion and now aggree that this is the better solution. The code segement below should the be changed to the following code:


function dists = distEucSq (x, y)
  xx = sumsq (x, 2);
  yy = sumsq (y, 2)';
  dists = max(0, bsxfun (@plus, xx, yy) - 2 * x * (y'));
endfunction


Georg Wiora <gwiora>
Wed 22 Feb 2017 09:55:15 AM UTC, comment #2: 

I thought about that too. I think the numeric result is in both cases equal to a precision of eps. But the abs function is in terms of memory use and complexity maybe more efficient than the solution using max().

Please correct me if my understanding is wrong. The max solution is definitely "more" precise in its result.

Georg Wiora <gwiora>
Wed 22 Feb 2017 09:47:51 AM UTC, comment #1: 

What about using max (., 0) instead of abs (.) ?

Guillaume <gyom>
Wed 22 Feb 2017 08:09:21 AM UTC, original submission:  

For certain cases the pdist2 function from statistics package may return complex results leading to unexpected behaviour in following calculations. The complex results are caused by a negative result for "dists" in the distEucSq() local function:


function dists = distEucSq (x, y)
  xx = sumsq (x, 2);
  yy = sumsq (y, 2)';
  dists = bsxfun (@plus, xx, yy) - 2 * x * (y');
endfunction


The negative results are caused by rounding errors and are in the range of -1e-16. A single negative value in dists array causes the whole result array in case of euclidian distance calculation (default for pdist2) to become complex data type.

Suggested solution:
Add an abs() function for squared distances_


function dists = distEucSq (x, y)
  xx = sumsq (x, 2);
  yy = sumsq (y, 2)';
  dists = abs(bsxfun (@plus, xx, yy) - 2 * x * (y'));
endfunction


I have appended a patched version of pdist2.m.

Georg Wiora <gwiora>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #39810:  pdist2.m added by gwiora (5KiB - application/octet-stream - Patched version of pdist2 function.)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by nrjank (Updated the item)
  • -email is unavailable- added by nir_krakauer (Posted a comment)
  • -email is unavailable- added by lostbard (Posted a comment)
  • -email is unavailable- added by jwe (Updated the item)
  • -email is unavailable- added by gyom (Posted a comment)
  • -email is unavailable- added by gwiora (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
    2021-12-16 nrjank StatusNone Fixed
        Open/ClosedOpen Closed
    2017-08-12 jwe Summarypdist2 function returns complex results [octave forge] (statistics) pdist2 function returns complex results
    2017-02-22 gwiora Attached File- Added pdist2.m, #39810

    Back to the top

    Powered by Savane 3.13-3230.
    Corresponding source code