Clear Filters
Clear Filters

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

6 views (last 30 days)
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.

Answers (1)

Christine Tobler
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.
  2 Comments
Jose Pujol
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'
Jose Pujol
Jose Pujol on 21 Oct 2023
% 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, nm, 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

Sign in to comment.

Categories

Find more on Linear Algebra in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!