hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_lapack_dggsvd (f08va)

Purpose

nag_lapack_dggsvd (f08va) computes the generalized singular value decomposition (GSVD) of an mm by nn real matrix AA and a pp by nn real matrix BB.

Syntax

[k, l, a, b, alpha, beta, u, v, q, iwork, info] = f08va(jobu, jobv, jobq, a, b, 'm', m, 'n', n, 'p', p)
[k, l, a, b, alpha, beta, u, v, q, iwork, info] = nag_lapack_dggsvd(jobu, jobv, jobq, a, b, 'm', m, 'n', n, 'p', p)
Note: the interface to this routine has changed since earlier releases of the toolbox:
Mark 23: iwork is now an output parameter
.

Description

The generalized singular value decomposition is given by
UT A Q = D1
(0R)
,   VT B Q = D2
(0R)
,
UT A Q = D1 0 R ,   VT B Q = D2 0 R ,
where UU, VV and QQ are orthogonal matrices. Let (k + l)(k+l) be the effective numerical rank of the matrix
(A)
B
A B , then RR is a (k + l)(k+l) by (k + l)(k+l) nonsingular upper triangular matrix, D1D1 and D2D2 are mm by (k + l)(k+l) and pp by (k + l)(k+l) ‘diagonal’ matrices structured as follows:
if mkl0m-k-l0,
D1 =
kl
k(I0)
l 0 C
mkl 0 0
D1= klk(I0) l 0 C m-k-l 0 0
D2 =
kl
l(0S)
pl 0 0
D2= kll(0S) p-l 0 0
(0R)
=
nklkl
k(0R11R12)
l 0 0 R22
0R = n-k-lklk(0R11R12) l 0 0 R22
where
C = diag(αk + 1,,αk + l) ,
C = diag(αk+1,,αk+l) ,
S = diag(βk + 1,,βk + l) ,
S = diag(βk+1,,βk+l) ,
and
C2 + S2 = I .
C2 + S2 = I .
RR is stored as a submatrix of AA with elements RijRij stored as Ai,nkl + jAi,n-k-l+j on exit.
If mkl < 0 m-k-l<0 ,
D1 =
kmkk + lm
k(I00)
mk 0 C 0
D1= km-kk+l-mk(I00) m-k 0 C 0
D2 =
kmkk + lm
mk(0S0)
k + lm 0 0 I
pl 0 0 0
D2= km-kk+l-mm-k(0S0) k+l-m 0 0 I p-l 0 0 0
(0R)
=
nklkmkk + lm
k(0R11R12R13)
mk 0 0 R22 R23
k + lm 0 0 0 R33
0R = n-k-lkm-kk+l-mk(0R11R12R13) m-k 0 0 R22 R23 k+l-m 0 0 0 R33
where
C = diag(αk + 1,,αm) ,
C = diag(αk+1,,αm) ,
S = diag(βk + 1,,βm) ,
S = diag(βk+1,,βm) ,
and
C2 + S2 = I .
C2 + S2 = I .
(R11R12R13)
0 R22 R23
R11 R12 R13 0 R22 R23  is stored as a submatrix of AA with RijRij stored as Ai,nkl + jAi,n-k-l+j, and R33 R33  is stored as a submatrix of BB with (R33)ij(R33)ij stored as Bmk + i,n + mkl + jBm-k+i,n+m-k-l+j.
The function computes CC, SS, RR and, optionally, the orthogonal transformation matrices UU, VV and QQ.
In particular, if BB is an nn by nn nonsingular matrix, then the GSVD of AA and BB implicitly gives the SVD of AB1AB-1:
A B1 = U (D1D21) VT .
A B-1 = U ( D1 D2-1 ) VT .
If
(A)
B
A B  has orthonormal columns, then the GSVD of AA and BB is also equal to the CS decomposition of AA and BB. Furthermore, the GSVD can be used to derive the solution of the eigenvalue problem:
AT Ax = λ BT Bx .
AT Ax=λ BT Bx .
In some literature, the GSVD of AA and BB is presented in the form
UT A X =
(0D1)
,   VT B X =
(0D2)
,
UT A X = 0D1 ,   VT B X = 0D2 ,
where UU and VV are orthogonal and XX is nonsingular, and D1D1 and D2D2 are ‘diagonal’. The former GSVD form can be converted to the latter form by taking the nonsingular matrix XX as
X = Q
(I0)
0 R1
.
X = Q I 0 0 R-1 .

References

Anderson E, Bai Z, Bischof C, Blackford S, Demmel J, Dongarra J J, Du Croz J J, Greenbaum A, Hammarling S, McKenney A and Sorensen D (1999) LAPACK Users' Guide (3rd Edition) SIAM, Philadelphia http://www.netlib.org/lapack/lug
Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore

Parameters

Compulsory Input Parameters

1:     jobu – string (length ≥ 1)
If jobu = 'U'jobu='U', the orthogonal matrix UU is computed.
If jobu = 'N'jobu='N', UU is not computed.
Constraint: jobu = 'U'jobu='U' or 'N''N'.
2:     jobv – string (length ≥ 1)
If jobv = 'V'jobv='V', the orthogonal matrix VV is computed.
If jobv = 'N'jobv='N', VV is not computed.
Constraint: jobv = 'V'jobv='V' or 'N''N'.
3:     jobq – string (length ≥ 1)
If jobq = 'Q'jobq='Q', the orthogonal matrix QQ is computed.
If jobq = 'N'jobq='N', QQ is not computed.
Constraint: jobq = 'Q'jobq='Q' or 'N''N'.
4:     a(lda, : :) – double array
The first dimension of the array a must be at least max (1,m)max(1,m)
The second dimension of the array must be at least max (1,n)max(1,n)
The mm by nn matrix AA.
5:     b(ldb, : :) – double array
The first dimension of the array b must be at least max (1,p)max(1,p)
The second dimension of the array must be at least max (1,n)max(1,n)
The pp by nn matrix BB.

Optional Input Parameters

1:     m – int64int32nag_int scalar
Default: The first dimension of the array a.
mm, the number of rows of the matrix AA.
Constraint: m0m0.
2:     n – int64int32nag_int scalar
Default: The second dimension of the array a.
nn, the number of columns of the matrices AA and BB.
Constraint: n0n0.
3:     p – int64int32nag_int scalar
Default: The first dimension of the array b.
pp, the number of rows of the matrix BB.
Constraint: p0p0.

Input Parameters Omitted from the MATLAB Interface

lda ldb ldu ldv ldq work

Output Parameters

1:     k – int64int32nag_int scalar
2:     l – int64int32nag_int scalar
k and l specify the dimension of the subblocks kk and ll as described in Section [Description]; (k + l)(k+l) is the effective numerical rank of
(A)
B
AB.
3:     a(lda, : :) – double array
The first dimension of the array a will be max (1,m)max(1,m)
The second dimension of the array will be max (1,n)max(1,n)
ldamax (1,m)ldamax(1,m).
Contains the triangular matrix RR, or part of RR. See Section [Description] for details.
4:     b(ldb, : :) – double array
The first dimension of the array b will be max (1,p)max(1,p)
The second dimension of the array will be max (1,n)max(1,n)
ldbmax (1,p)ldbmax(1,p).
Contains the triangular matrix RR if mkl < 0m-k-l<0. See Section [Description] for details.
5:     alpha(n) – double array
See the description of beta.
6:     beta(n) – double array
alpha and beta contain the generalized singular value pairs of AA and BB, αi αi and βi βi ;
  • alpha(1 : k) = 1 alpha1:k = 1 ,
  • beta(1 : k) = 0 beta1:k = 0 ,
and if mkl0 m-k-l0 ,
  • alpha(k + 1 : k + l) = C alphak+1:k+l = C ,
  • beta(k + 1 : k + l) = S betak+1:k+l = S ,
or if mkl < 0 m-k-l<0 ,
  • alpha(k + 1 : m) = C alphak+1:m = C ,
  • alpha(m + 1 : k + l) = 0 alpham+1:k+l = 0 ,
  • beta(k + 1 : m) = S betak+1:m = S ,
  • beta(m + 1 : k + l) = 1 betam+1:k+l = 1 , and
  • alpha(k + l + 1 : n) = 0 alphak+l+1:n = 0 ,
  • beta(k + l + 1 : n) = 0 betak+l+1:n = 0 .
The notation alpha(k : n)alphak:n above refers to consecutive elements alpha(i)alphai, for i = k,,ni=k,,n.
7:     u(ldu, : :) – double array
The first dimension, ldu, of the array u will be
  • if jobu = 'U'jobu='U', ldu max (1,m) ldu max(1,m) ;
  • otherwise ldu1ldu1.
The second dimension of the array will be max (1,m)max(1,m) if jobu = 'U'jobu='U', and at least 11 otherwise
If jobu = 'U'jobu='U', u contains the mm by mm orthogonal matrix UU.
If jobu = 'N'jobu='N', u is not referenced.
8:     v(ldv, : :) – double array
The first dimension, ldv, of the array v will be
  • if jobv = 'V'jobv='V', ldv max (1,p) ldv max(1,p) ;
  • otherwise ldv1ldv1.
The second dimension of the array will be max (1,p)max(1,p) if jobv = 'V'jobv='V', and at least 11 otherwise
If jobv = 'V'jobv='V', v contains the pp by pp orthogonal matrix VV.
If jobv = 'N'jobv='N', v is not referenced.
9:     q(ldq, : :) – double array
The first dimension, ldq, of the array q will be
  • if jobq = 'Q'jobq='Q', ldq max (1,n) ldq max(1,n) ;
  • otherwise ldq1ldq1.
The second dimension of the array will be max (1,n)max(1,n) if jobq = 'Q'jobq='Q', and at least 11 otherwise
If jobq = 'Q'jobq='Q', q contains the nn by nn orthogonal matrix QQ.
If jobq = 'N'jobq='N', q is not referenced.
10:   iwork(n) – int64int32nag_int array
Stores the sorting information. More precisely, the following loop will sort alpha
 for i=k+1:min(m,k+l) % swap alpha(i) and alpha(iwork(i)) end 
such that alpha(1)alpha(2)alpha(n)alpha1alpha2alphan.
11:   info – int64int32nag_int scalar
info = 0info=0 unless the function detects an error (see Section [Error Indicators and Warnings]).

Error Indicators and Warnings

  info = iinfo=-i
If info = iinfo=-i, parameter ii had an illegal value on entry. The parameters are numbered as follows:
1: jobu, 2: jobv, 3: jobq, 4: m, 5: n, 6: p, 7: k, 8: l, 9: a, 10: lda, 11: b, 12: ldb, 13: alpha, 14: beta, 15: u, 16: ldu, 17: v, 18: ldv, 19: q, 20: ldq, 21: work, 22: iwork, 23: info.
It is possible that info refers to a parameter that is omitted from the MATLAB interface. This usually indicates that an error in one of the other input parameters has caused an incorrect value to be inferred.
  INFO = 1INFO=1
If info = 1info=1, the Jacobi-type procedure failed to converge.

Accuracy

The computed generalized singular value decomposition is nearly the exact generalized singular value decomposition for nearby matrices (A + E) (A+E)  and (B + F) (B+F) , where
E2 = O(ε) A2 ​ and ​ F2 = O(ε) B2 ,
E2 = O(ε) A2 ​ and ​ F2 = O(ε) B2 ,
and ε ε  is the machine precision. See Section 4.12 of Anderson et al. (1999) for further details.

Further Comments

The complex analogue of this function is nag_lapack_zggsvd (f08vn).

Example

function nag_lapack_dggsvd_example
jobu = 'U';
jobv = 'V';
jobq = 'Q';
a = [1, 2, 3;
     3, 2, 1;
     4, 5, 6;
     7, 8, 8];
b = [-2, -3, 3;
      4,  6, 5];
% Compute the generalized singular value decomposition of (a, b)
% (a = U*D1*(0 R)*(Q^T), b = V*D2*(0 R)*(Q^T), m>=n)
[k, l, a, b, alpha, beta, u, v, q, iwork, info] = ...
    nag_lapack_dggsvd(jobu, jobv, jobq, a, b);

if (info == 0)
  irank = double(k + l);
  kk = double(k);
  fprintf('\nNumber of infinite generalized singular values (k)\n%5d\n', k);
  fprintf('Number of finite generalized singular values (l)\n%5d\n', l);
  fprintf('Numerical rank of (a^T b^T)^T (k+l)\n%5d\n', irank);

  fprintf('\nFinite generalized singular values\n');
  disp(transpose((alpha(kk+1:irank)./beta(kk+1:irank))));

  fprintf('Orthogonal matrix u\n');
  disp(u);

  fprintf('Orthogonal matrix v\n');
  disp(v);

  fprintf('Orthogonal matrix q\n');
  disp(q);

  fprintf('Non singular upper triangular matrix R\n');
  disp(triu(a(1:3,3-irank+1:3)));

  % Estimate the reciprocal condition number of R
  [rcond, info] = nag_lapack_dtrcon('Infinity-norm','Upper','Non-unit', a(1:3,3-irank+1:3));

  fprintf('Estimate of reciprocal condition number for R\n%11.1e\n', rcond);

  % So long as irank = n, get the machine precision, eps, from nag_machine_precision and compute
  % the approximate error bound for the computed generalized singular values
  if (irank==3)
    serrbd = nag_machine_precision/rcond;
    fprintf('\nError estimate for the generalized singular values\n%11.1e\n', serrbd);
  else
    fprintf('\n(A^T B^T)^T is not of full rank\n');
  end
end
 

Number of infinite generalized singular values (k)
    1
Number of finite generalized singular values (l)
    2
Numerical rank of (a^T b^T)^T (k+l)
    3

Finite generalized singular values
    1.3151    0.0802

Orthogonal matrix u
   -0.1348    0.5252   -0.2092    0.8137
    0.6742   -0.5221   -0.3889    0.3487
    0.2697    0.5276   -0.6578   -0.4650
    0.6742    0.4161    0.6101    0.0000

Orthogonal matrix v
    0.3554   -0.9347
    0.9347    0.3554

Orthogonal matrix q
   -0.8321   -0.0946   -0.5466
    0.5547   -0.1419   -0.8199
         0   -0.9853    0.1706

Non singular upper triangular matrix R
   -2.0569   -9.0121   -9.3705
         0  -10.8822   -7.2688
         0         0   -6.0405

Estimate of reciprocal condition number for R
    4.2e-02

Error estimate for the generalized singular values
    2.6e-15

function f08va_example
jobu = 'U';
jobv = 'V';
jobq = 'Q';
a = [1, 2, 3;
     3, 2, 1;
     4, 5, 6;
     7, 8, 8];
b = [-2, -3, 3;
      4,  6, 5];
% Compute the generalized singular value decomposition of (a, b)
% (a = U*D1*(0 R)*(Q^T), b = V*D2*(0 R)*(Q^T), m>=n)
[k, l, a, b, alpha, beta, u, v, q, iwork, info] = ...
    f08va(jobu, jobv, jobq, a, b);

if (info == 0)
  irank = double(k + l);
  kk = double(k);
  fprintf('\nNumber of infinite generalized singular values (k)\n%5d\n', k);
  fprintf('Number of finite generalized singular values (l)\n%5d\n', l);
  fprintf('Numerical rank of (a^T b^T)^T (k+l)\n%5d\n', irank);

  fprintf('\nFinite generalized singular values\n');
  disp(transpose((alpha(kk+1:irank)./beta(kk+1:irank))));

  fprintf('Orthogonal matrix u\n');
  disp(u);

  fprintf('Orthogonal matrix v\n');
  disp(v);

  fprintf('Orthogonal matrix q\n');
  disp(q);

  fprintf('Non singular upper triangular matrix R\n');
  disp(triu(a(1:3,3-irank+1:3)));

  % Estimate the reciprocal condition number of R
  [rcond, info] = f07tg('Infinity-norm','Upper','Non-unit', a(1:3,3-irank+1:3));

  fprintf('Estimate of reciprocal condition number for R\n%11.1e\n', rcond);

  % So long as irank = n, get the machine precision, eps, from x02aj and compute
  % the approximate error bound for the computed generalized singular values
  if (irank==3)
    serrbd = x02aj/rcond;
    fprintf('\nError estimate for the generalized singular values\n%11.1e\n', serrbd);
  else
    fprintf('\n(A^T B^T)^T is not of full rank\n');
  end
end
 

Number of infinite generalized singular values (k)
    1
Number of finite generalized singular values (l)
    2
Numerical rank of (a^T b^T)^T (k+l)
    3

Finite generalized singular values
    1.3151    0.0802

Orthogonal matrix u
   -0.1348    0.5252   -0.2092    0.8137
    0.6742   -0.5221   -0.3889    0.3487
    0.2697    0.5276   -0.6578   -0.4650
    0.6742    0.4161    0.6101    0.0000

Orthogonal matrix v
    0.3554   -0.9347
    0.9347    0.3554

Orthogonal matrix q
   -0.8321   -0.0946   -0.5466
    0.5547   -0.1419   -0.8199
         0   -0.9853    0.1706

Non singular upper triangular matrix R
   -2.0569   -9.0121   -9.3705
         0  -10.8822   -7.2688
         0         0   -6.0405

Estimate of reciprocal condition number for R
    4.2e-02

Error estimate for the generalized singular values
    2.6e-15


PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013