# Proof of relation between the generalized singular values of gsvd(A,B) and gsvd(B,A)

2 views (last 30 days)
Jose Pujol on 26 Jan 2023
Commented: Jose Pujol on 8 Feb 2023
First a comment. The gsvd documentation has the following sentence in the description of "sigma":
"When B is square and nonsingular, the generalized singular values, gsvd(A,B), correspond to the ordinary singular values, svd(A/B), but they are sorted in the opposite order. THEIR reciprocals are gsvd(B,A). "
I capitalized THEIR because I think it is ambiguous. It is unclear to me whether THEIR refers to "gsvd(A,B)" or to "svd(A/B)".
Now the question. Numerical comparison of gsvd(A,B) and gsvd(B,A) shows that their generalized singular values are reciprocal and in opposite order. This is a heuristic result that requires certain relations between the two pairs of matrices U,V,X,C,S corresponding to gsvd(A,B) and gsvd(B,A). Heuristically, I found those relations, but I wonder whether there is a theoretical analysis of my question.

Christine Tobler on 31 Jan 2023
The background for this is the 5-output form of the GSVD:
[U,V,X,C,S] = gsvd(A,B) returns unitary matrices U and V, a (usually) square matrix X, and nonnegative diagonal matrices C and S so that A = U*C*X', B = V*S*X', C'*C + S'*S = I.
In practice, here are the outputs for an example:
rng default; % Let's not have these change on every run
A = randn(3);
B = randn(3);
[U1, V1, X1, C1, S1] = gsvd(A, B)
U1 = 3×3
0.1977 -0.8555 -0.4785 0.9784 0.2027 0.0418 0.0612 -0.4764 0.8771
V1 = 3×3
0.5201 0.3942 -0.7577 -0.2706 -0.7654 -0.5839 0.8101 -0.5087 0.2914
X1 = 3×3
4.6140 1.1460 -2.2033 1.0532 -0.0580 -1.5759 1.2268 -1.4670 3.4249
C1 = 3×3
0.3819 0 0 0 0.8620 0 0 0 0.9811
S1 = 3×3
0.9242 0 0 0 0.5069 0 0 0 0.1933
Now say we want to swap the roles of A and B here - we can simply switch the roles of U and V, and those of C and S, and we'll have formulas that satisfy all three equations above.
U2 = V1;
V2 = U1;
C2 = S1;
S2 = C1;
X2 = X1;
norm(B - U2*C2*X2')
ans = 2.7914e-15
norm(A - V2*S2*X2')
ans = 2.1166e-15
norm(C2'*C2 + S2'*S2 - eye(3))
ans = 2.2204e-16
But we've ignored the fact that C1 had decreasing singular values and S1 had increasing values - so we'll want to flip the order of the rows and columns of C2 and S2, and make corresponding changes to U2, V2 and X2:
C2 = diag(flip(diag(C2)));
S2 = diag(flip(diag(S2)));
U2 = flip(U2, 2);
V2 = flip(V2, 2);
X2 = flip(X2, 2);
norm(B - U2*C2*X2')
ans = 2.8354e-15
norm(A - V2*S2*X2')
ans = 1.8216e-15
norm(C2'*C2 + S2'*S2 - eye(3))
ans = 2.2204e-16
Let's compare what we've constructed here to the output from gsvd with A and B swapped out:
[U3, V3, X3, C3, S3] = gsvd(B, A);
norm(C2 - C3)
ans = 4.4409e-16
norm(S2 - S3)
ans = 3.3307e-16
So we've shown that swapping the order of A and B will swap out the outputs C and S, and flip their order. In the one-output case (and when A and B are both square), this output is simply
gsvd(A, B)
ans = 3×1
0.4132 1.7005 5.0770
diag(C1)./diag(S1)
ans = 3×1
0.4132 1.7005 5.0770
It follows that if A and B are swapped, their inverses are returned, and the order of those inverses is flipped so that the numbers are again in increasing order.
Jose Pujol on 8 Feb 2023
% I was hoping for a theoretical answer because I am working on a book
% that will discuss the GSVD in detail and I am not comfortable stating
% a result without knowing whether there is a theoretical proof of it, or
% it is a heuristic result. The relation between the singular values of
% gsvd(A,B)and gsvd(B,A) is particularly important and I would think it
% could be derived by analysis of the gsvd algorithm. The relations
% between U1 and V2 and U2 and V1 are trickier and not as simple as
% indicated in the answer, as my code below shows.
% BACKGROUND
% From
% [U1,V1,X1,C1,S1] = gsvd(A,B)
% and
% [U2,V2,X2,C2,S2] = gsvd(B,A)
% we get
% A=U1*C1*X1=V2*S2*X2 (1)
% and
% B=V1*S1*X1=U2*C2*X2 (2)
% QUESTION: is it OK to write:
% V2=flip(U1) (3)
% S2=flip(C1) (4)
% X2=flip(X1) (5)
% A=flip(U1)*flip(C1)*flip(X1)' (6) (from eq. 1)
% and
% U2=flip(V1) (7)
% C2=flip(S1) (8)
% B=flip(V1)*flip(S1)*flip(X1)' (9) (from eq. 2)?
% ANSWER: in general, NO [except for equations (4) and (8), as expected].
% The code below has three examples. This is a summary of results:
% Example 1: A and B are square and nonsingular and the equations (3),
% (5), (7) are valid only in absolute value. Equations (6) and (9) apply
% as written.
% Example 2: A and B are rectangular and of full rank and the equations
% (3), (5), (7) only apply to subsets of columns. Equations (6) and (9)
% apply only when using the same subsets of columns.
% Example 3: A and B are square but A has two zero eigenvalues and the
% equations (3), (5), (7) only apply to subsets of columns in absolute
% value. Equations (6) and (9) apply % as written.
%EXAMPLE 1
fprintf('--------- EXAMPLE 1')
clear all
m=6, p=6, n=6
sseed=11, rand('seed',sseed);
A=rand(m,n);
B=rand(p,n);
[U1,V1,X1,C1,S1] = gsvd(A,B);
alpha1=max(C1,[],1);
beta1=max(S1,[],1);
[U2,V2,X2,C2,S2] = gsvd(B,A);
alpha2=max(C2,[],1);
beta2=max(S2,[],1);
U1flip=fliplr(U1);
V1flip=fliplr(V1);
U2flip=fliplr(U2);
V2flip=fliplr(V2);
X1flip=fliplr(X1);
X2flip=fliplr(X2);
%Norm applied to eqs. (3), (5), (7)
norm_V2_U1flip=norm(V2-U1flip) %this is equal to 2
norm_X2_X1flip=norm(X2-X1flip) %this is equal to 5.4
norm_U2_V1flip=norm(U2-V1flip) %this is equal to 2
%Norm applied to eqs. (3), (5), (7) in absolute value --- All ~= 1e-14
%Some columns have opposite signs
norm_V2_U1flipabs=norm(abs(V2)-abs(U1flip))
norm_X2_X1flipabs=norm(abs(X2)-abs(X1flip))
norm_U2_V1flipabs=norm(abs(U2)-abs(V1flip))
%Eqs. (6) and (9) and two variations are satisfied. All elements
%smaller than ~1e-14
C1flip=diag(fliplr(alpha1));
A-U1flip*C1flip*X1flip'
S1flip=diag(fliplr(beta1));
B-V1flip*S1flip*X1flip'
S2flip=diag(fliplr(beta2));
A-V2flip*S2flip*X2flip'
C2flip=diag(fliplr(alpha2));
B-U2flip*C2flip*X2flip'
%EXAMPLE 2
fprintf('--------- EXAMPLE 2')
clear all
m=6, p=5, n=4
sseed=11, rand('seed',sseed);
A=rand(m,n);
B=rand(p,n);
[U1,V1,X1,C1,S1] = gsvd(A,B);
alpha1=max(C1,[],1);
beta1=max(S1,[],1);
[U2,V2,X2,C2,S2] = gsvd(B,A);
alpha2=max(C2,[],1);
beta2=max(S2,[],1);
U1flip=fliplr(U1);
V1flip=fliplr(V1);
X1flip=fliplr(X1);
%Norm applied to eqs. (3), (5), (7)
norm_V2_U1flip=norm(V2-U1flip) %this is equal to 2
norm_X2_X1flip=norm(X2-X1flip) %this is equal to 5.1e-015
norm_U2_V1flip=norm(U2-V1flip) %this is equal to 1.9
%USE subsets of columns:
norm_V2_U1flipsub=norm(V2(:,1:4)-U1flip(:,3:6)) % ~= 1e-15
norm_U2_V1flipsub=norm(U2(:,1:4)-V1flip(:,2:5)) % ~= 1e-15
%Eqs. (6) and (9) and two variations are NOT SATISFIED ---
%All elements between about 0.02 and 2 in abs. %value
U1flip=fliplr(U1);
C1flip=C1;
C1flip(1:n,1:n)=diag(fliplr(alpha1));
fprintf('A-U1flip*C1flip*X1flip'':\n')
A-U1flip*C1flip*X1flip'
V1flip=fliplr(V1);
S1flip=S1;
S1flip(1:n,1:n)=diag(fliplr(beta1));
fprintf('B-V1flip*S1flip*X1flip'':\n')
B-V1flip*S1flip*X1flip'
S2flip=S2;
V2flip=fliplr(V2);
X2flip=fliplr(X2);
S2flip(1:n,1:n)=diag(fliplr(beta2));
fprintf('A-V2flip*S2flip*X2flip'':\n')
A-V2flip*S2flip*X2flip'
C2flip=C2;
U2flip=fliplr(U2);
C2flip(1:n,1:n)=diag(fliplr(alpha2));
fprintf('B-U2flip*C2flip*X2flip'':\n')
B-U2flip*C2flip*X2flip'
%Eqs. (6) and (9) and two variations are SATISFIED for subsets of columns
%All elements smalller than 1e-14 or 1e-15 in absolute value
fprintf('A-U1flip(:,3:6)*C1flip(1:4,1:4)*X1flip'':\n')
A-U1flip(:,3:6)*C1flip(1:4,1:4)*X1flip'
fprintf('B-V1flip(:,2:5)*S1flip(1:4,1:4)*X1flip'':\n')
B-V1flip(:,2:5)*S1flip(1:4,1:4)*X1flip'
fprintf('A-V2flip(:,3:6)*S2flip(1:4,1:4)*X2flip'':\n')
A-V2flip(:,3:6)*S2flip(1:4,1:4)*X2flip'
fprintf('B-U2flip(:,2:5)*C2flip(1:4,1:4)*X2flip'':\n')
B-U2flip(:,2:5)*C2flip(1:4,1:4)*X2flip'
%EXAMPLE 3
fprintf('--------- EXAMPLE 3')
clear all
m=6, p=6, n=6
sseed=11, rand('seed',sseed);
A1=rand(m,n);
B=rand(p,n);
% Create a matrix A with two singular values equal to zero
[UA1,LA1,VA1]=svd(A1);
LA1(1,1)=0;
LA1(n,n)=0;
A=UA1*LA1*VA1;
[U1,V1,X1,C1,S1] = gsvd(A,B);
alpha1=max(C1,[],1);
beta1=max(S1,[],1);
[U2,V2,X2,C2,S2] = gsvd(B,A);
alpha2=max(C2,[],1);
beta2=max(S2,[],1);
U1flip=fliplr(U1);
U2flip=fliplr(U2);
V1flip=fliplr(V1);
V2flip=fliplr(V2);
X1flip=fliplr(X1);
X2flip=fliplr(X2);
norm_V1_U2flip=norm(V1-U2flip) %this is equal to 2
norm_X2_X1flip=norm(X2-X1flip) %this is equal to 2.8
norm_V2_U1flip=norm(V2-U1flip) %this is equal to 2
%USE subsets of columns and absolute values.
%Some columns have opposite signs
norm_V1_U2flipsubabs=norm(abs(V1(:,3:6))-abs(U2flip(:,3:6))) % = 1.9e-14
norm_X1_X2flipsubabs=norm(abs(X1(:,3:6))-abs(X2flip(:,3:6))) % = 1.1e-14
norm_V2_U1flipsubabs=norm(abs(V2(:,1:4))-abs(U1flip(:,1:4))) % = 1.6e-14
%Eqs. (6) and (9) and two variations --- All elements smaller
%than 1e-15 or 1e-14
C1flip=diag(fliplr(alpha1));
A-fliplr(U1)*C1flip*X1flip'
S1flip=diag(fliplr(beta1));
B-fliplr(V1)*S1flip*X1flip'
S2flip=diag(fliplr(beta2));
A-V2flip*S2flip*X2flip'
C2flip=diag(fliplr(alpha2));
B-U2flip*C2flip*X2flip'