bugGNU Octave - Bugs: bug #60273, Wrong result in GSVD

 
 

bug #60273: Wrong result in GSVD

Submitted by:  Leonardo <pocker>
Submitted on:  Tue 23 Mar 2021 11:15:02 AM UTC  
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Confirmed 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)

Tue 23 Mar 2021 08:42:31 PM UTC, comment #12: 

You can try tests in gdvd.cc-tst. Or Example 1 in Matlab's doc.
Essentially they are not 0-padded.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Tue 23 Mar 2021 08:35:51 PM UTC, comment #11: 

I have other pending changes for function handles that will affect the API and that should probably go in the next stable release (assuming there is one) so I don't object to also making the GSVD changes on stable.  It makes sense to me to have the same (or similar) interface for all the functions, so changing liboctave is probably best.  Can we also preserve (but deprecate) the previous interface?

John W. Eaton <jwe>
Project Administrator
Tue 23 Mar 2021 08:23:51 PM UTC, comment #10: 

Solution #3 : change on stable in libinterp to preserve liboctave API, and then change in liboctave on dev branch.

Rik <rik5>
Project Administrator
Tue 23 Mar 2021 08:21:46 PM UTC, comment #9: 

@jwe: How do you think this bug should be addressed?

Solution #1 : Change liboctave version of gsvd to return Matlab-compatible 5 outputs.

Solution #2: Change libinterp version of gsvd to calculate the 5 Matlab-compatible return outputs from the 6 outputs provided by liboctave.

It seems like the change should go on the stable branch because Octave is definitely not producing the output that it says it is, hence this is more than just Matlab-compatibility issue which could wait until next release.

Rik <rik5>
Project Administrator
Tue 23 Mar 2021 08:17:56 PM UTC, comment #8: 

@Dmitri: Do you have an example that shows the size of the matrices are different from Matlab?  Always helpful if we can see what is supposed to happen.  There is already some code to adjust the size of the output matrices so it might just require tweaking of that.

Rik <rik5>
Project Administrator
Tue 23 Mar 2021 06:36:33 PM UTC, comment #7: 

> B = V * S * R * Q' (and not V*C*R*Q')


Yes, just a typo.  I meant what you wrote.

Rik <rik5>
Project Administrator
Tue 23 Mar 2021 06:27:23 PM UTC, comment #6: 

I think the size of S and C matrices is not correct (i.e. not Matlab compatible) when A and B are not square.

Dmitri.
-- 

Dmitri A. Sergatskov <dasergatskov>
Tue 23 Mar 2021 06:08:28 PM UTC, comment #5: 

Indeed.

>> [U, V, Q, C, S, R] = gsvd (A, B);
>> U*C*R*Q'

ans =

   1.0000   1.0000   3.0000
   3.0000   1.0000   2.0000
   2.0000   1.0000   1.0000

>> V*S*R*Q'

ans =

   2.0000e+00   3.8858e-16   1.0000e+00
   2.0000e+00   1.0000e+00  -1.0000e+00
   1.0000e+00   7.6328e-17  -9.3675e-17

and finally, as you pointed out,

>> Q*R' - A'*U*inv(C)

ans =

   8.8818e-16   1.7764e-15  -3.1225e-16
   2.2204e-16   6.6613e-16  -1.1102e-16
  -8.8818e-16   1.7451e-15   4.4409e-16

which confirms that

X = Q*R'

And I do agree that the documentation of GSVD should be adjusted.

BTW, it is

B = V * S * R * Q' (and not V*C*R*Q')

Thanks on this!

Leonardo <pocker>
Tue 23 Mar 2021 05:43:32 PM UTC, comment #4: 

Octave is definitely taking the matrix Q from LAPACK and returning it as the matrix X from gsvd which is incorrect.  Documentation on the routine that Octave calls is at https://www.math.utah.edu/software/lapack/lapack-d/dggsvd.html.

The variable naming in gsvd.cc shows that the original coder knew that 'q' was being returned.  Maybe there wasn't a realization that Matlab was not exactly mapping the outputs of LAPACK to the outputs of the Matlab function.

In any case, there is a 6th undocumented output from Octave's gsvd function which is the matrix R from LAPACK.  The documentation should really read

[U, V, Q, C, S, R] = gsvd (A, B)

If you do this then the results conform to the LAPACK documentation and

A = U * C * R * Q'
B = V * C * R * Q'

To recover Matlab's behavior

X = Q*R'

Rik <rik5>
Project Administrator
Tue 23 Mar 2021 04:56:30 PM UTC, comment #3: 

From

https://www.netlib.org/lapack/lug/node36.html

it seems that matrix X is just the matrix Q because

>> inv(C)*U'*A*X

ans =

   2.6026e+00  -3.5966e+00  -2.3962e+00
  -5.3213e-16  -3.8168e+00   5.0481e-01
  -7.7716e-16  -8.8818e-16   1.6511e+00

is upper triangular.

Also, matrix X is orthogonal (in general, it should not be)

>> X'*X

ans =

   1.0000e+00  -1.8041e-16  -5.5511e-17
  -1.8041e-16   1.0000e+00   1.2490e-16
  -5.5511e-17   1.2490e-16   1.0000e+00

Leonardo <pocker>
Tue 23 Mar 2021 04:41:40 PM UTC, comment #2: 
Dmitri A. Sergatskov <dasergatskov>
Tue 23 Mar 2021 04:16:58 PM UTC, comment #1: 

Confirmed.

The problem is only with matrix X as stated.  The code in question is in liboctave/numeric/gsvd.cc.  The code there calls LAPACK to perform the actual gsvd calculation.

There are several possibilities.  1) Octave is not calling LAPACK correctly, 2) there is a bug in LAPACK, 3) or matrix X is meant to be calculated outside of LAPACK via some other method such as the suggested

X = A'*U*inv(C)

Rik <rik5>
Project Administrator
Tue 23 Mar 2021 11:15:02 AM UTC, original submission:  

Hi.

From "help gsvd" it says that

[U, V, X, C, S] = gsvd (A, B)
A = U*C*X'
B = V*S*X'

But,

>> A = [1 1 3; 3 1 2; 2 1 1];
>> B = [2 0 1; 2 1 -1;1 0 0];
>> [U, V, X, C, S] = gsvd (A, B);
>> U*C*X'

ans =

0.206864 0.733236 0.430549
-0.510638 0.269817 0.312529
-0.294408 0.404150 -0.080575

which is clearly not matrix A.

The problem is with matrix X

>> X

X =

-0.186817 -0.981925 0.030379
0.865745 -0.149940 0.477496
0.464311 -0.115506 -0.878108

while it should be

>> X = A'*U*inv(C)

X =

2.9726e+00 3.7631e+00 5.0158e-02
1.6483e+00 8.1334e-01 7.8837e-01
3.7280e+00 -2.4142e-03 -1.4498e+00

After this, I get

>> U*C*X'

ans =

1.0000 1.0000 3.0000
3.0000 1.0000 2.0000
2.0000 1.0000 1.0000

which is matrix A and

>> V*S*X'

ans =

2.0000 0 1.0000
2.0000 1.0000 -1.0000
1.0000 -0.0000 -0.0000

which is matrix B.

This affects, at least, v. 6.2.0 and 7.0.0 (4a516dce916c+). Both compiled from source, debian 10.

Kind regards,
Leonardo

Leonardo <pocker>

 

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

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by rik5
  • -email is unavailable- added by dasergatskov (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by pocker (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 2 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-03-23 rik5 Carbon-Copy- => Added jwe
    2021-03-23 rik5 StatusNone => Confirmed

    Back to the top


    Powered by Savane 3.6