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

7 visualizaciones (últimos 30 días)
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.

Respuestas (1)

Christine Tobler
Christine Tobler el 31 de En. de 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 comentarios
Jose Pujol
Jose Pujol el 8 de Feb. de 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 el 21 de Oct. de 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

Iniciar sesión para comentar.

Categorías

Más información sobre Linear Algebra en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by