possible issue dggsvd3 function

Post here if you want to report a bug to the LAPACK team

possible issue dggsvd3 function

Postby tiagofrepereira » Thu Jan 12, 2017 1:00 pm

I'm trying to use the the generalized SVD function dggsvd3 in C and I'm having the following issue.

In a very simplified way I'm creating the following function:
[U,V,Q,C,S] = my_dggsvd3_in_C(A,B),

where A and B are my inputs and they are:
A = [[ 1., 1., 2.],
[ 0., 1., 3.],
[ 0., 4., 3.]])

B = [[ 0.2, 1.4, 2. ],
[ 0.4, 1. , 3.2],
[ 0.3, 4.5, 3.3]]

With these inputs I have the following outputs:
U = [[-0.45904757, 0.22321234, -0.8599137 ],
[-0.88154988, -0.23451243, 0.40972397],
[-0.11020501, 0.94613961, 0.30442518]]

V = [[-0.38233328, 0.1538968 , 0.91111856],
[-0.92047714, -0.149751 , -0.36096602],
[-0.0808894 , 0.97667313, -0.19891328]]

C = [ 0.69954115, 0.64582169, 0.99967773]

S = [ 0.71459232, 0.76348828, 0.02538585]

Q = [[-0.09327252, -0.55839886, 0.82431241],
[ 0.70757717, 0.54528323, 0.44944494],
[ 0.70045328, -0.6251855 , -0.34425035]]

My problem is with the output of Q.
This output is completely different than the one from MATLAB.
In MATLAB, for the same inputs, Q has the following value.

Q =
0.3456 -0.6562 0.8602
5.8426 -2.5466 -0.7678
3.9969 -5.5656 -0.4228

Can someone please explain me what is going on?
I can attach my C code if needed.

Thanks in advance.

Tiago
tiagofrepereira
 
Posts: 3
Joined: Thu Jan 12, 2017 12:40 pm

Re: possible issue dggsvd3 function

Postby resultant » Thu Mar 30, 2017 12:16 pm

I show the result of Scilab
C=[0.99967773 0 0 ; 0 0.69954115 0 ; 0 0 0.64582169]
S=[0.02538585 0 0 ; 0 0.71459232 0 ; 0 0 0.76348828]
U=[- 0.8599137 - 0.4590476 0.2232123;
0.4097240 - 0.8815499 - 0.2345124 ;
0.3044252 - 0.1102050 0.9461396]
V=[0.9111186 - 0.3823333 0.1538968 ;
- 0.3609660 - 0.9204771 - 0.149751;
- 0.1989133 - 0.0808894 0.9766731 ]
R=[ [- 0.8839363 - 0.1972057 - 0.8293800] ,
[0. - 3.1809123 5.2700225 ],
[0. 0. - 7.0873518]]
and
Q^T=[[0.9908940 0.0310158 - 0.1310230],
[0.1255023 - 0.5652072 0.8153465 ]
[- 0.0487665 - 0.8243657 - 0.563953]],
where R and Q are defined in ww.netlib.org/lapack/lug/node36.html/.
A and B are square. I can compute svd(A*B^(-1))=U*Sigma*V^{T}.
Then I can check U,V,C,S. U,V,C,S are true.
Also, I can check C^(-1)*U^{T}*A=S^(-1)*V^{T}*B.
From C^(-1)*U^{T}*A, we can obtain R*Q^{T}=C^(-1)*U^{T}*A.
Of course, Q*R^{T}=(C^(-1)*U^{T}*A)^{T}.
We use F=[[0 0 1],[0 1 0],[1 0 0]].
F*Q*F*F*R^{T}*F=F*(C^(-1)*U^{T}*A)^{T}*F.
QR decomposition Q0*R0 <= F*(C^(-1)*U^{T}*A)^{T}*F
Thus, Q0=F*Q*F(F*Q0*F=Q) and R0=F*R^{T}*F (same).
Now, we have U,V,C,S,R,Q.

At least, Q of MATLAB is not orthogonal.
Last edited by resultant on Fri Mar 31, 2017 9:32 pm, edited 2 times in total.
resultant
 
Posts: 13
Joined: Fri Nov 13, 2015 9:35 am

Re: possible issue dggsvd3 function

Postby resultant » Fri Mar 31, 2017 2:22 am

As another approach, from QR decomposition Q*R^{-1}=B^{-1}*V*S,
we can obtain Q and R^{-1}. Moreover, we compute R.
Last edited by resultant on Fri Mar 31, 2017 9:33 pm, edited 1 time in total.
resultant
 
Posts: 13
Joined: Fri Nov 13, 2015 9:35 am

Re: possible issue dggsvd3 function

Postby resultant » Fri Mar 31, 2017 3:29 pm

C^(-1)*U^{T}*A=R*Q^{T}=
[- 0.8601909 0.7677584 0.4227563 ;
- 0.6562124 - 2.5465514 - 5.5655908 ;
0.3456253 5.8425699 3.9969333 ]
Thus, MATLAB returns Q*R^{T}.

If we set,
C=[ 0.69954115 0 0 ; 0 0.64582169 0; 0 0 0.99967773]
S=[0.71459232 0 0 ; 0 0.76348828 0 ; 0 0 0.02538585]
U=[-0.45904757, 0.22321234, -0.8599137 ; -0.88154988, -0.23451243, 0.40972397 ; -0.11020501, 0.94613961, 0.30442518 ]
V=[-0.38233328, 0.1538968 , 0.91111856 ; -0.92047714, -0.149751 , -0.36096602 ; -0.0808894 , 0.97667313, -0.19891328]
Q=[-0.09327252, -0.55839886, 0.82431241 ; 0.70757717, 0.54528323, 0.44944494 ; 0.70045328, -0.6251855 , -0.34425035],
C^(-1)*U^(-1)*A*Q=S^(-1)*V^(-1)*B*Q=
[ - 3.1045876 - 4.3543393 3.0483774
9.640D-09 5.2268135 - 4.7865413
3.085D-09 - 2.153D-08 - 1.228049 ].
Thus, dggsvd3.f returns correct values.
resultant
 
Posts: 13
Joined: Fri Nov 13, 2015 9:35 am


Return to Bug report

Who is online

Users browsing this forum: No registered users and 0 guests

cron