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'

[Jose Pujol]

% To understand the relation between GSVD(A,B) and GSVD(B,A) it is necessary to underestand the

% derivation of the GSVD, which relies on the CS decomposition (or CSD), which is discussed first.

% Consider the vertically partitioned matrix M = [A;B], where A is m1 × n, B is m2 × n, m1 ≥ n, and

% m2 ≥ n. These assumptions are not essential. The following is based on Algorithm 1 ofVan Loan

% (1985, Numer. Math., 479-491). The economy-size SVD of M will be written as

% M = [A;B]= Q*D*Z’ = [Q1; Q1]*D*Z’=[Q1*D*Z’;Q2*D*Z’], (1)

% where Q = [Q1;Q2], with Q1 and Q2 equal to the first m1 and last m2 rows of Q.

% Because the columns of Q constitute an orthonormal set

% Q’*Q = Q1’*Q1 + Q2’*Q2 = In=eye(n). (2)

% The singular values of Q are all equal to 1, while those of Q1 and Q2 are between 0 and 1.

% The CS decomposition of [Q1;Q2] is given by

% Q1 = R1*C*W’ (3)

% Q2 = R2*S*W’, (4)

% where R1 is m1 × m1, R2 is m2 × m2, and W is n × n. These three matrices are orthonormal. The

% matrices C are S are m1 × n and m2 × n and diagonal, with elements c_i and s_i, respectively. These% two matrices satisfy

% C’*C + S’*S = In (5)

% Eqs. (3) and (4) are the SVDs of Q1 and Q2. The first step in the derivation of the CSD is

% Eq. (4). The goal is to find expressions for R1 and C. We need the following results and steps

% (Q2*W)’*Q2*W = S’*S. (6)

% Introduce the m1 × n matrix X1 (assumed to be nonsingular)

% X1 = Q1*W. (7)

% X1’*X1 = W’*W − (Q2*W)’*Q2*W = I − S’*S. (8)

% X1’*X1 = diag(1 − s^2_i ) ≡ diag(c^2_i ), c^2_i ≡ 1 − s^2_i , c^2_i ≤ 1, i = 1, ..., n. (9)

% Introduce the n × n matrix C

% C = diag(c_i), i = 1, ..., n. (10)

% Introduce the m1 × n matrix R1

% R1 = X1*inv(C) = Q1*W*inv(C) ⇒ Q1 = R1*C*W’. (11)

% This completes the derivation of the CSD.

% Eq.(5) follows from Eqs. (8)-(10).

% The GSVD is derived in terms of the CSD of Q. The matrices C and S introduced above

% and the GSVD matrices C and S have the same properties. From Eq. (1)

% A = Q1*D*Z’ , B = Q2*D*Z’ . (12)

% Using Eqs. (3) and (4) the GSVD equations:

% A = R1*C*X’ , B = R2*S*X’, X = Z*D*W, (13)

% The following is an alternative derivation of the CSD, motivated by the previous one. It was used to

% establish the relation between GSVD(A,B) and GSVD(B,A). The starting point is the economy-size

% SVDsof the Q1 and Q2 introduced in Eq. (1), which now will be written as

% Q1= U1*Λ1*V1’, Q2= U2*Λ2*V2’. (14)

% Q1 and U1 are m1 × n, Q2 and U2 are m2 × n, and the other matrices are n × n. V2=W. We have

% Q’*Q = In = V1*Λ1^2*V1’ + V2*Λ2^2*V2’ . (15)

% This expression will be satisfied if V1 = V2 and Λ1^2 + Λ2^2 = In. Λ1 and Λ2 should

% have properties similar to those of S and C. This requires flipping the order of the elements of Λ1,

% and the columns of U1 and V1. Let Λ1r, U1r and V1r be the corresponding matrices. Then,

% U1r*Λ1r*V1r’=Q1. (16)

% This result does not mean that V1r will be computationally equal to V2. In fact, hundreds of

% thousands of computer runs involving different random matrices show that

% V1r = V2*Sg1, Λ1^2+ Λ2^2 = In. (17)

% The matrix Sg1 is diagonal with nonzero elements equal to either 1 or −1 and is needed because

% some of the columns of V2 and V1r have opposite signs. This ambiguity in the signs can be verified % using the Matlab functions svd and svds (Bro et al, Sandia Report SAND2007-6422).

% To satisfy the SVD of Q1 when the columns of V1r have a change of sign it is necessary to

% apply the same change to the columns of U1r. After these changes, Eq. (16) is recovered.

% In terms of V2 we have

% Q1’*Q1 = V2*Λ1r^2*V2’. (18)

% Then the second of Eqs. (17) show that the second equality in Eqs. (15) is satisfied.

% The following are expressions for A and B in terms of the previous results:

% A = U1r*Λ1r*V1r’*D*Z’, B = U2*Λ2*V2’* D*Z’ ≡ U2*Λ2*X’ . (19)

% Compare these results to the economy version of GSVD(A, B), written as

% A = Ug*C*Xg’, B = Vg*S*Xg’. (20)

% The numerical computations described in the preceding also produce these results

% Ug = U1r*Sg2, C = Λ1r, Xg*Sg2 = Z*D*V1r, (21)

% Vg = U2*Sg3, S = Λ2, Xg*Sg3 = Z*D*V2, (22)

% where Sg2 and Sg3 are similar toSg1 and Sg3=Sg1*Sg2 . Then

% A = Ug*Sg2*C*Sg2*Xg’ = Ug*C*Xg’, B = Vg*Sg3*S*Sg3*Xg’ = Vg*S*Xg’, (23)

% in agreement with Eqs. (20).

% The relation between GSVD(A,B) and GSVD(B,A) can be derived using Eqs. (19), (21) and

% (22).

% The final result for the case m1=m≥n, m2=p≥n, A and B nonsingular is as follows. Let

% [U1,V1,X1,C1,S1] = gsvd(A,B,0), alpha1=max(C1,[],1), beta1=max(S1,[],1),

% [U2,V2,X2,C2,S2] = gsvd(B,A,0), alpha2=max(C2,[],1), beta2=max(S2,[],1).

% U2flip=fliplr(U2), V2flip=fliplr(V2), X2flip=fliplr(X2),

% alpha2flip=fliplr(alpha2), beta2flip=fliplr(beta2).

% For Example 1 we have the following relations, whichapply to the n columns:

% U1=V2flip*Sg1, V1=U2flip*Sg1, X1=X2flip*Sg1, alpha1=beta2flip, beta1=alpha2flip (24)

% For Example 2, matrix A is singular. Eqs. (24) only apply to a subset of columns.

% For example 3, n≥ m, p. Matrices A and B are nonsingular. Eqs. (24) apply to different

% subsets of columns but with the same number of columns.

clear all

close all

time=clock;

DATE=datestr(datenum(time), 'mmm.dd,yyyy HH:MM'); %:SS')

fid=fopen('Matlab_compare_GSVD_out.txt','w');

fprintf(fid,'OUTPUT OF Matlab_compare_GSVD.m\n');

fprintf(fid,'DATE: %s\n\n',DATE);

niter=1 %number of iterations

%EXAMPLE 1. m>=n, p>=n. Full rank matrices.

fprintf('--------- EXAMPLE 1')

clearvars -except fid niter

m=48, p=47, n=40

sseed=11, rand('seed',sseed);

for ii=1:niter

A=rand(m,n);

B=rand(p,n);

[U1,V1,X1,C1,S1] = gsvd(A,B,0);

alpha1=max(C1,[],1);

beta1=max(S1,[],1);

[U2,V2,X2,C2,S2] = gsvd(B,A,0);

alpha2=max(C2,[],1);

beta2=max(S2,[],1);

U1flip=fliplr(U1);

V1flip=fliplr(V1);

X1flip=fliplr(X1);

U2flip=fliplr(U2);

V2flip=fliplr(V2);

X2flip=fliplr(X2);

alpha1flip=fliplr(alpha1);

beta1flip=fliplr(beta1);

sg1=sign(U1(1,1:n)./V2flip(1,1:n)); %vector of sign differences in the columns

Sg1=diag(sg1); %matrix of sign differences

Sg1r=diag(fliplr(sg1));

if max(max(U1(:,1:n)-V2flip(:,1:n)*Sg1))>1.0e-10

fprintf('Ex. 1; max(max(U1(:,1:n)-V2flip(:,1:n)*Sg1)) --- Should be close to the machine 0:\n')

max(max(U1(:,1:n)-V2flip(:,1:n)*Sg1))

ii, 1

return

end

if max(max(V1(:,1:n)-U2flip(:,1:n)*Sg1))>1.0e-10

fprintf('Ex. 1; max(max(V1(:,1:n)-U2flip(:,1:n)*Sg1)) --- Should be close to the machine 0:\n')

max(max(V1(:,1:n)-U2flip(:,1:n)*Sg1))

ii, 2

return

end

if max(max(X1(:,1:n)-X2flip(:,1:n)*Sg1))>1.0e-10

fprintf('Ex. 1; max(max(X1(:,1:n)-X2flip(:,1:n)*Sg1)) --- Should be close to the machine 0:\n')

max(max(X1(:,1:n)-X2flip(:,1:n)*Sg1))

ii, 3

return

end

sg2=sign(U2(1,1:n)./V1flip(1,1:n)); %vector of sign differences in the columns

Sg2=diag(sg2); %matrix of sign differences

if sg2-fliplr(sg1) ~=0

fprintf('Ex. 1; sg2-(fliplr(sg1) ~=0 ... RETURN\n')

ii, 4

end

if max(sg1-fliplr(sg2))>1.e-10

fprintf('Ex. 1; max(sg1-fliplr(sg2)) --- Should be close to the machine 0 \n')

max(sg1-fliplr(sg2))

ii, 5

return

end

if max(max(U2(:,1:n)-V1flip(:,1:n)*Sg2))>1.0e-10

fprintf('Ex. 1; max(max(U2(:,1:n)-V1flip(:,1:n)*Sg2)) --- Should be close to the machine 0 \n')

max(max(U2(:,1:n)-V1flip(:,1:n)*Sg2))

ii, 6

return

end

if max(max(V2(:,1:n)-U1flip(:,1:n)*Sg2))>1.0e-10

fprintf('Ex. 1; max(max(V2(:,1:n)-U1flip(:,1:n)*Sg2)) --- Should be close to the machine 0:\n')

max(max(V2(:,1:n)-U1flip(:,1:n)*Sg2))

ii, 7

return

end

if max(max(X2(:,1:n)-X1flip(:,1:n)*Sg2))>1.0e-10

fprintf('Ex. 1; max(max(X2(:,1:n)-X1flip(:,1:n)*Sg2)) --- Should be close to the machine 0: \n')

max(max(X2(:,1:n)-X1flip(:,1:n)*Sg2))

ii, 8

return

end

if max(alpha1flip-beta2)>1.e-10

fprintf('Ex. 1; max(alpha1flip-beta2) --- Should be close to the machine 0 \n')

max(beta1flip-alpha2)

max(alpha1flip-beta2)

ii, 9

return

end

if max(beta1flip-alpha2)>1.e-10

fprintf('Ex. 1; max(beta1flip-alpha2) --- Should be close to the machine 0 \n')

max(beta1flip-alpha2)

ii, 10

return

end

end

%=======================================================================

%EXAMPLE 2. m>=n, p>=n. Matrix A rank deficient.

fprintf('--------- EXAMPLE 2')

clearvars -except fid niter

m=40, p=49, n=38

sseed=119, rand('seed',sseed);

for ii=1:niter;

A1=rand(m,n);

B=rand(p,n);

% Create a matrix A with three singular values equal to zero

[UA1,LA1,VA1]=svd(A1);

LA1(1,1)=0;

fn2=floor(n/2);

LA1(fn2,fn2)=0;

LA1(n,n)=0;

A=UA1*LA1*VA1;

[U1,V1,X1,C1,S1] = gsvd(A,B,0);

alpha1=max(C1,[],1);

beta1=max(S1,[],1);

[U2,V2,X2,C2,S2] = gsvd(B,A,0);

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);

ind=alpha1<1.e-10;

mm=sum(ind); %number of leading zeros in C1 and trailing zeros in S2

%Some columns have opposite signs. Use columns without leading zeros

sg1=sign(U1(1,mm+1:n)./V2flip(1,mm+1:n)); %vector of sign differences in the columns

Sg1=diag(sg1); %matrix of sign differences

if max(max((U1(:,mm+1:n))-V2flip(:,mm+1:n)*Sg1))>1.e-10;

fprintf('Example 2 --- 1')

ii

return

end

if max(max((V1(:,mm+1:n))-U2flip(:,mm+1:n)*Sg1))>1.e-10;

fprintf('Example 2 --- 2')

ii

return

end

if max(max((X1(:,mm+1:n))-X2flip(:,mm+1:n)*Sg1))>1.e-10;

fprintf('Example 2 --- 3')

ii

return

end

sg2=sign(U2(1,1:n-mm)./V1flip(1,1:n-mm)); %vector of sign differences in the columns

Sg2=diag(sg2); %matrix of sign differences

if max(sg1-fliplr(sg2))>1.e-10

fprintf('Ex. 2; max(sg1-fliplr(sg2)) --- Should be close to the machine 0 \n')

max(sg1-fliplr(sg2))

ii, 4

return

end

if max(max((U2(:,1:n-mm))-V1flip(:,1:n-mm)*Sg2))>1.e-10;

fprintf('Example 2 --- 5')

ii

return

end

if max(max((V2(:,1:n-mm))-U1flip(:,1:n-mm)*Sg2))>1.e-10;

fprintf('Example 2 --- 6')

ii

return

end

if max(max((X2(:,1:n-mm))-X1flip(:,1:n-mm)*Sg2))>1.e-10;

fprintf('Example 2 --- 7')

ii

return

end

end

%=======================================================================

%EXAMPLE 3. m and p<=n. Full rank matrices.

fprintf('------------ Example 3\n')

%clearvars -except fid niter

clearvars -except fid

niter=1

sseed=11; rand('seed',sseed);

for ii=1:niter

for jj=1:2

if jj==1;

%m=1410; p=1470; n=1500;

m=11; p=13; n=18;

else

%m=470; p=410; n=500;

m=13; p=11; n=18;

end

if niter==1

fprintf(fid,'ii=%d, m=%d, p=%d, n=%d, seed=%d\n \n',ii,m,p,n,sseed);

end

A=rand(m,n);

B=rand(p,n);

[U1,V1,X1,C1,S1] = gsvd(A,B);

alpha1=max(C1,[],1);

beta1=max(S1,[],1);

if niter==1

fprintf(fid,'Matrix A size: %d %d\n',size(A));

fprintf(fid,'Matrix B size: %d %d\n\n',size(B));

fprintf(fid,'Matrix U1 size: %d %d\n',size(U1));

fprintf(fid,'Matrix C1 size: %d %d\n',size(C1));

fprintf(fid,'Matrix X1'' size: %d %d\n',size(X1'));

fprintf(fid,'Matrix V1 size: %d %d\n',size(V1));

fprintf(fid,'Matrix S1 size: %d %d\n',size(S1));

fprintf(fid,'Matrix X1'' size: %d %d\n\n',size(X1'));

end

[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);

if niter==1

fprintf(fid,'Matrix U2 size: %d %d\n',size(U2));

fprintf(fid,'Matrix C2 size: %d %d\n',size(C2));

fprintf(fid,'Matrix X2'' size: %d %d\n',size(X2'));

fprintf(fid,'Matrix V2 size: %d %d\n',size(V2));

fprintf(fid,'Matrix S2 size: %d %d\n',size(S2));

fprintf(fid,'Matrix X2'' size: %d %d\n\n',size(X2'));

end

ind=alpha1<1.e-10;

mm=sum(ind); %number of leading zeros in C1 and trailing zeros in S2, =n-m

k=0 %USE AS A TEST: change k=0 to k=1. This should not work

mm=mm-k

if niter==1

fprintf(fid,'mm= number of leading zeros in C1 and trailing zeros in S2: %d\n',mm);

end

%Some columns have opposite signs.

sg1=sign(U1(1,1:p-mm)./V2flip(1,1:p-mm)); %vector of sign differences in the columns

Sg1=diag(sg1); %matrix of sign differences

%fprintf(fid,'Matrix Sg1 size: %d %d\n\n',size(Sg1));

if niter==1

fprintf(fid,'Matrix U1(:,1:p-mm) size: %d %d\n\n',size(U1(:,1:p-mm)));

end

if max(max(U1(:,1:p-mm)-V2flip(:,1:p-mm)*Sg1))>1.0e-8

fprintf('max(max(U1(:,1:p-mm)-V2flip(:,1:p-mm)*Sg1)) --- Should be close to the machine 0 \n');

max(max(U1(:,1:p-mm)-V2flip(:,1:p-mm)*Sg1))

1

return

end

if max(alpha1-fliplr(beta2))>1.0e-8

fprintf('Ex. 3; max(alpha1-fliplr(beta2)) --- Should be close to the machine 0 \n')

max(alpha1-fliplr(beta2))

ii, 2

return

end

if max(alpha2-fliplr(beta1))>1.0e-8

fprintf('Ex. 3; max(alpha2-fliplr(beta1)) --- Should be close to the machine 0 \n')

max(alpha1-fliplr(beta2))

ii, 2

return

end

if niter==1

fprintf(fid,'Matrices V1 and X1 (:,mm+1:p) sizes: %d %d\n\n',size(V1(:,mm+1:p)));

end

if max(max(V1(:,mm+1:p)-U2flip(:,mm+1:p)*Sg1))>1.0e-8

fprintf('Ex. 3; max(max(V1(:,mm+1:p)-U2flip(:,mm+1:p)*Sg1)) --- Should be close to the machine 0 \n')

max(max(V1(:,mm+1:p)-U2flip(:,mm+1:p)*Sg1))

ii, 2

return

end

if niter==1

fprintf(fid,'Matrices V1 and X1 (:,mm+1:p) sizes: %d %d\n\n',size(V1(:,mm+1:p)));

end

if max(max(X1(:,mm+1:p)-X2flip(:,mm+1:p)*Sg1))>1.0e-8

fprintf('Ex. 3; max(max(X1(:,mm+1:m)-X2flip(:,mm+1:m)*Sg1)) --- Should be close to the machine 0 \n')

max(max(X1(:,mm+1:p)-X2flip(:,mm+1:p)*Sg1))

ii, 3

return

end

sg2=sign(U2(1,1:p-mm)./V1flip(1,1:1:p-mm)); %vector of sign differences in the columns

Sg2=diag(sg2); %matrix of sign differences

%fprintf(fid,'Matrix Sg2 size: %d %d\n',size(Sg2));

if max(sg1-fliplr(sg2))>1.e-8

fprintf('Ex. 3; max(sg1-fliplr(sg2)) --- Should be close to the machine 0 \n')

max(sg1-fliplr(sg2))

ii, 4

return

end

%fprintf(fid,'Matrix fliplr(Sg2) = Matrix Sg1\n\n');

if niter==1

fprintf(fid,'Matrix U2(:,1:p-mm) size: %d %d\n\n',size(U2(:,1:p-mm)));

end

if max(max(U2(:,1:p-mm)-V1flip(:,1:p-mm)*Sg2))>1.0e-8

fprintf('Ex. 3; max(max(U2(:,1:m-mm)-V1flip(:,1:m-mm)*Sg2)) --- Should be close to the machine 0 \n');

max(max(U2(:,1:m-mm)-V1flip(:,1:m-mm)*Sg2))

ii, 5

return

end

if m>=p

if niter==1

fprintf(fid,'--------- m>=p ---------\n');

end

ind2=alpha2<1.e-10;

mm2=sum(ind2); %number of leading zeros in C2 and trailing zeros in S1, =n-p

if niter==1

fprintf(fid,'mm2= number of leading zeros in C2 and trailing zeros in S1: %d\n\n',mm2);

fprintf(fid,'Matrices V2 and X2 (:,mm2+1:m) sizes: %d %d\n\n',size(V2(:,mm2+1:m)));

end

if max(max(V2(:,mm2+1:m)-U1flip(:,mm2+1:m)*Sg2))>1.0e-8

fprintf('Ex. 3; max(max(V2(:,mm2+1:m)-U1flip(:,mm2+1:m)*Sg2)) --- Should be close to the machine 0 \n');

max(max(V2(:,mm2+1:m)-U1flip(:,mm2+1:m)*Sg2))

ii, 6

return

end

if max(max(X2(:,mm2+1:m)-X1flip(:,mm2+1:m)*Sg2))>1.0e-8

fprintf('Ex. 3; max(max(X2(:,mm2+1:m)-X1flip(:,mm2+1:m)*Sg2)) --- Should be close to the machine 0 \n');

max(max(X2(:,mm2+1:m)-X1flip(:,mm2+1:m)*Sg2))

ii, 7

return

end

elseif m<p

if niter==1

fprintf(fid,'--------- m<p ---------\n');

end

np=n-p;

%fprintf(fid,'np=n-p: %d\n\n',np);

sg3=sign(V2(1,np+1:m)./U1flip(1,np+1:m)); %vector of sign differences in the columns =sg2

Sg3=diag(sg3); %matrix of sign differences

if max(sg3-sg2)>1.e-8

fprintf('Ex. 3; max(sg3-sg2) --- Should be close to the machine 0 \n')

max(sg3-sg2)

ii, 8

return

end

%fprintf(fid,'Matrix Sg3 = Matrix Sg2\n\n');

if niter==1

fprintf(fid,'Matrices V2 and X2 (:,np+1:m) sizes: %d %d\n\n',size(V2(:,np+1:m)));

end

if max(max(V2(:,np+1:m)-U1flip(:,np+1:m)*Sg3))>1.0e-8

fprintf('Ex. 3; max(max(V2(:,np+1:m)-U1flip(:,np+1:m)*Sg3)) --- Should be close to the machine 0 \n');

max(max(V2(:,np+1:m)-U1flip(:,np+1:m)*Sg3))

ii, 9

return

end

if max(max(X2(:,np+1:m)-X1flip(:,np+1:m)*Sg3))>1.0e-8

fprintf('Ex. 3; max(max(X2(:,np+1:m)-X1flip(:,np+1:m)*Sg3)) --- Should be close to the machine 0 \n');

max(max(X2(:,np+1:m)-X1flip(:,np+1:m)*Sg3))

ii, 10

return

end

end

%fprintf(fid,'================================================= \n');

end

end

fclose(fid);

fclose('all'); % The fclose above are not sufficient to open the files