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

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.

Answer 1

[Christine Tobler]

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]

% 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'