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_zbdsqr (f08ms)

Purpose

nag_lapack_zbdsqr (f08ms) computes the singular value decomposition of a complex general matrix which has been reduced to bidiagonal form.

Syntax

[d, e, vt, u, c, info] = f08ms(uplo, d, e, vt, u, c, 'n', n, 'ncvt', ncvt, 'nru', nru, 'ncc', ncc)
[d, e, vt, u, c, info] = nag_lapack_zbdsqr(uplo, d, e, vt, u, c, 'n', n, 'ncvt', ncvt, 'nru', nru, 'ncc', ncc)

Description

nag_lapack_zbdsqr (f08ms) computes the singular values and, optionally, the left or right singular vectors of a real upper or lower bidiagonal matrix BB. In other words, it can compute the singular value decomposition (SVD) of BB as
B = U Σ VT .
B = U Σ VT .
Here ΣΣ is a diagonal matrix with real diagonal elements σiσi (the singular values of BB), such that
σ1 σ2 σn 0 ;
σ1 σ2 σn 0 ;
UU is an orthogonal matrix whose columns are the left singular vectors uiui; VV is an orthogonal matrix whose rows are the right singular vectors vivi. Thus
Bui = σi vi   and   BT vi = σi ui ,   i = 1,2,,n .
Bui = σi vi   and   BT vi = σi ui ,   i = 1,2,,n .
To compute UU and/or VTVT, the arrays u and/or vt must be initialized to the unit matrix before nag_lapack_zbdsqr (f08ms) is called.
The function stores the real orthogonal matrices UU and VTVT in complex arrays u and vt, so that it may also be used to compute the SVD of a complex general matrix AA which has been reduced to bidiagonal form by a unitary transformation: A = QBPHA=QBPH. If AA is mm by nn with mnmn, then QQ is mm by nn and PHPH is nn by nn; if AA is nn by pp with n < pn<p, then QQ is nn by nn and PHPH is nn by pp. In this case, the matrices QQ and/or PHPH must be formed explicitly by nag_lapack_zungbr (f08kt) and passed to nag_lapack_zbdsqr (f08ms) in the arrays u and/or vt respectively.
nag_lapack_zbdsqr (f08ms) also has the capability of forming UHCUHC, where CC is an arbitrary complex matrix; this is needed when using the SVD to solve linear least squares problems.
nag_lapack_zbdsqr (f08ms) uses two different algorithms. If any singular vectors are required (i.e., if ncvt > 0ncvt>0 or nru > 0nru>0 or ncc > 0ncc>0), the bidiagonal QRQR algorithm is used, switching between zero-shift and implicitly shifted forms to preserve the accuracy of small singular values, and switching between QRQR and QLQL variants in order to handle graded matrices effectively (see Demmel and Kahan (1990)). If only singular values are required (i.e., if ncvt = nru = ncc = 0ncvt=nru=ncc=0), they are computed by the differential qd algorithm (see Fernando and Parlett (1994)), which is faster and can achieve even greater accuracy.
The singular vectors are normalized so that ui = vi = 1ui=vi=1, but are determined only to within a complex factor of absolute value 11.

References

Demmel J W and Kahan W (1990) Accurate singular values of bidiagonal matrices SIAM J. Sci. Statist. Comput. 11 873–912
Fernando K V and Parlett B N (1994) Accurate singular values and differential qd algorithms Numer. Math. 67 191–229
Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore

Parameters

Compulsory Input Parameters

1:     uplo – string (length ≥ 1)
Indicates whether BB is an upper or lower bidiagonal matrix.
uplo = 'U'uplo='U'
BB is an upper bidiagonal matrix.
uplo = 'L'uplo='L'
BB is a lower bidiagonal matrix.
Constraint: uplo = 'U'uplo='U' or 'L''L'.
2:     d( : :) – double array
Note: the dimension of the array d must be at least max (1,n)max(1,n).
The diagonal elements of the bidiagonal matrix BB.
3:     e( : :) – double array
Note: the dimension of the array e must be at least max (1,n1)max(1,n-1).
The off-diagonal elements of the bidiagonal matrix BB.
4:     vt(ldvt, : :) – complex array
The first dimension, ldvt, of the array vt must satisfy
  • if ncvt > 0ncvt>0, ldvt max (1,n) ldvt max(1,n) ;
  • otherwise ldvt1ldvt1.
The second dimension of the array must be at least max (1,ncvt)max(1,ncvt)
If ncvt > 0ncvt>0, vt must contain an nn by ncvtncvt matrix. If the right singular vectors of BB are required, ncvt = nncvt=n and vt must contain the unit matrix; if the right singular vectors of AA are required, vt must contain the unitary matrix PHPH returned by nag_lapack_zungbr (f08kt) with vect = 'P'vect='P'.
5:     u(ldu, : :) – complex array
The first dimension of the array u must be at least max (1,nru)max(1,nru)
The second dimension of the array must be at least max (1,n)max(1,n)
If nru > 0nru>0, u must contain an nrunru by nn matrix. If the left singular vectors of BB are required, nru = nnru=n and u must contain the unit matrix; if the left singular vectors of AA are required, u must contain the unitary matrix QQ returned by nag_lapack_zungbr (f08kt) with vect = 'Q'vect='Q' .
6:     c(ldc, : :) – complex array
The first dimension, ldc, of the array c must satisfy
  • if ncc > 0ncc>0, ldc max (1,n) ldc max(1,n) ;
  • otherwise ldc1ldc1.
The second dimension of the array must be at least max (1,ncc)max(1,ncc)
The nn by nccncc matrix CC if ncc > 0ncc>0.

Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The dimension of the arrays d, u.
nn, the order of the matrix BB.
Constraint: n0n0.
2:     ncvt – int64int32nag_int scalar
Default: The second dimension of the array vt.
ncvtncvt, the number of columns of the matrix VHVH of right singular vectors. Set ncvt = 0ncvt=0 of right singular vectors. Set ncvt = 0ncvt=0 if no right singular vectors are required.
Constraint: ncvt0ncvt0.
3:     nru – int64int32nag_int scalar
Default: The first dimension of the array u.
nrunru, the number of rows of the matrix UU of left singular vectors. Set nru = 0nru=0 if no left singular vectors are required.
Constraint: nru0nru0.
4:     ncc – int64int32nag_int scalar
Default: The second dimension of the array c.
nccncc, the number of columns of the matrix CC. Set ncc = 0ncc=0 if no matrix CC is supplied.
Constraint: ncc0ncc0.

Input Parameters Omitted from the MATLAB Interface

ldvt ldu ldc work

Output Parameters

1:     d( : :) – double array
Note: the dimension of the array d must be at least max (1,n)max(1,n).
The singular values in decreasing order of magnitude, unless INFO > 0INFO>0 (in which case see Section [Error Indicators and Warnings]).
2:     e( : :) – double array
Note: the dimension of the array e must be at least max (1,n1)max(1,n-1).
e is overwritten, but if INFO > 0INFO>0 see Section [Error Indicators and Warnings].
3:     vt(ldvt, : :) – complex array
The first dimension, ldvt, of the array vt will be
  • if ncvt > 0ncvt>0, ldvt max (1,n) ldvt max(1,n) ;
  • otherwise ldvt1ldvt1.
The second dimension of the array will be max (1,ncvt)max(1,ncvt)
The nn by ncvtncvt matrix VH VH or VH VH of right singular vectors, stored by rows.
If ncvt = 0ncvt=0, vt is not referenced.
4:     u(ldu, : :) – complex array
The first dimension of the array u will be max (1,nru)max(1,nru)
The second dimension of the array will be max (1,n)max(1,n)
ldu max (1,nru) ldu max(1,nru) .
The nrunru by nn matrix UU or QUQU of left singular vectors, stored as columns of the matrix.
If nru = 0nru=0, u is not referenced.
5:     c(ldc, : :) – complex array
The first dimension, ldc, of the array c will be
  • if ncc > 0ncc>0, ldc max (1,n) ldc max(1,n) ;
  • otherwise ldc1ldc1.
The second dimension of the array will be max (1,ncc)max(1,ncc)
c stores the matrix UHCUHC. If ncc = 0ncc=0, c is not referenced.
6:     info – int64int32nag_int scalar
info = 0info=0 unless the function detects an error (see Section [Error Indicators and Warnings]).

Error Indicators and Warnings

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

  info = iinfo=-i
If info = iinfo=-i, parameter ii had an illegal value on entry. The parameters are numbered as follows:
1: uplo, 2: n, 3: ncvt, 4: nru, 5: ncc, 6: d, 7: e, 8: vt, 9: ldvt, 10: u, 11: ldu, 12: c, 13: ldc, 14: work, 15: 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.
W INFO > 0INFO>0
The algorithm failed to converge and info specifies how many off-diagonals did not converge. In this case, d and e contain on exit the diagonal and off-diagonal elements, respectively, of a bidiagonal matrix orthogonally equivalent to BB.

Accuracy

Each singular value and singular vector is computed to high relative accuracy. However, the reduction to bidiagonal form (prior to calling the function) may exclude the possibility of obtaining high relative accuracy in the small singular values of the original matrix if its singular values vary widely in magnitude.
If σiσi is an exact singular value of BB and σ̃iσ~i is the corresponding computed value, then
|σ̃iσi| p (m,n) ε σi
| σ~i - σi | p (m,n) ε σi
where p(m,n)p(m,n) is a modestly increasing function of mm and nn, and εε is the machine precision. If only singular values are computed, they are computed more accurately (i.e., the function p(m,n)p(m,n) is smaller), than when some singular vectors are also computed.
If uiui is an exact left singular vector of BB, and iu~i is the corresponding computed left singular vector, then the angle θ(i,ui)θ(u~i,ui) between them is bounded as follows:
θ (i,ui) ( p (m,n) ε )/(relgapi)
θ (u~i,ui) p (m,n) ε relgapi
where relgapirelgapi is the relative gap between σiσi and the other singular values, defined by
relgapi = min ( |σiσj| )/((σi + σj)).
ij
relgapi = min ij | σi - σj | ( σi + σj ) .
A similar error bound holds for the right singular vectors.

Further Comments

The total number of real floating point operations is roughly proportional to n2n2 if only the singular values are computed. About 12n2 × nru12n2×nru additional operations are required to compute the left singular vectors and about 12n2 × ncvt12n2×ncvt to compute the right singular vectors. The operations to compute the singular values must all be performed in scalar mode; the additional operations to compute the singular vectors can be vectorized and on some machines may be performed much faster.
The real analogue of this function is nag_lapack_dbdsqr (f08me).

Example

function nag_lapack_zbdsqr_example
uplo = 'Upper';
d = [-3.087005021051958;
     2.066039276679068;
     1.873128891125711;
     2.002182866206992];
e = [2.112571007455839;
     1.262810106655224;
     -1.612633872800393];
vt = [complex(1),  0 + 0i,  0 + 0i,  0 + 0i;
      0 + 0i,  -0.2312345316176021 - 0.5404263814542949i, ...
     0.1786240755181332 - 0.5887280243319772i,  -0.4047984350247399 - 0.3348147213441872i;
      0 + 0i,  0.496225984754381 + 0.2648076085564618i, ...
    -0.3556231052094441 - 0.6888167652038528i,  0.2756552444973658 - 0.0819424169862172i;
      0 + 0i,  0.2917790897284148 + 0.5029628047064867i, ...
    0.04850708544077981 + 0.1349202975418933i,  -0.7155819195589249 - 0.3595545469781386i];
u = [ -0.3109810296560065 + 0.2623902437722555i, ...
    -0.6521001566525803 - 0.5532329309191114i,  -0.04266522773914702 - 0.03605518759134559i, ...
      0.2634342917361948 + 0.07411348303298232i;
      0.3174598011071733 - 0.6413983736655133i, ...
    -0.3487944910656935 - 0.07205673398466937i,  -0.2287385684714657 - 0.006908841577660436i, ...
      -0.1101402299076372 + 0.03261803709640307i;
      -0.2008419149861708 + 0.1490117433768365i, ...
    0.310251489958586 - 0.02303091424485477i,  -0.1854601676256336 + 0.181728049653941i, ...
      0.295608864759608 - 0.5647819449422838i;
      0.1198572718465858 - 0.1230966575721692i, ...
    0.004591621691001381 + 0.0004641170819608892i,  0.3304662059847925 - 0.4820711498283887i, ...
      0.06749320884731419 - 0.346397021614317i;
      -0.2688690152234222 - 0.1652086720047534i, ...
    -0.1793802257421584 + 0.05860743927321595i,  0.5235372995239013 + 0.2579777845606151i, ...
      -0.3926950316237215 - 0.1449707406311504i;
      -0.3498536583630073 + 0.09070280031633525i, ...
    -0.08288497844277821 + 0.05058867982255799i,  -0.3202392526929317 - 0.3037968928492991i, ...
      -0.3173573238191287 - 0.3241353242370706i];
c = [];
[dOut, eOut, vtOut, uOut, cOut, info] = nag_lapack_zbdsqr(uplo, d, e, vt, u, c)
 

dOut =

    3.9994
    3.0003
    1.9944
    0.9995


eOut =

     0
     0
     0


vtOut =

  -0.6971 + 0.0000i  -0.0867 - 0.3548i   0.0560 - 0.5400i  -0.1878 - 0.2253i
   0.2403 + 0.0000i   0.0725 - 0.2336i  -0.2477 - 0.5291i   0.7026 + 0.2177i
  -0.5123 + 0.0000i  -0.3030 - 0.1735i   0.0678 + 0.5162i   0.4418 + 0.3864i
  -0.4403 + 0.0000i   0.5294 + 0.6361i  -0.3027 - 0.0346i   0.1667 + 0.0258i


uOut =

  -0.5634 + 0.0016i  -0.2687 - 0.2749i   0.2451 + 0.4657i   0.3787 + 0.2987i
   0.1205 - 0.6108i  -0.2909 + 0.1085i   0.4329 - 0.1758i  -0.0182 - 0.0437i
  -0.0816 + 0.1613i  -0.1660 + 0.3885i  -0.4667 + 0.3821i  -0.0800 - 0.2276i
   0.1441 - 0.1532i   0.1984 - 0.1737i  -0.0034 + 0.1555i   0.2608 - 0.5382i
  -0.2487 - 0.0926i   0.6253 + 0.3304i   0.2643 - 0.0194i   0.1002 + 0.0140i
  -0.3758 + 0.0793i  -0.0307 - 0.0816i   0.1266 + 0.1747i  -0.4175 - 0.4058i


cOut =

     []


info =

                    0


function f08ms_example
uplo = 'Upper';
d = [-3.087005021051958;
     2.066039276679068;
     1.873128891125711;
     2.002182866206992];
e = [2.112571007455839;
     1.262810106655224;
     -1.612633872800393];
vt = [complex(1),  0 + 0i,  0 + 0i,  0 + 0i;
      0 + 0i,  -0.2312345316176021 - 0.5404263814542949i, ...
     0.1786240755181332 - 0.5887280243319772i,  -0.4047984350247399 - 0.3348147213441872i;
      0 + 0i,  0.496225984754381 + 0.2648076085564618i, ...
    -0.3556231052094441 - 0.6888167652038528i,  0.2756552444973658 - 0.0819424169862172i;
      0 + 0i,  0.2917790897284148 + 0.5029628047064867i, ...
    0.04850708544077981 + 0.1349202975418933i,  -0.7155819195589249 - 0.3595545469781386i];
u = [ -0.3109810296560065 + 0.2623902437722555i, ...
    -0.6521001566525803 - 0.5532329309191114i,  -0.04266522773914702 - 0.03605518759134559i, ...
      0.2634342917361948 + 0.07411348303298232i;
      0.3174598011071733 - 0.6413983736655133i, ...
    -0.3487944910656935 - 0.07205673398466937i,  -0.2287385684714657 - 0.006908841577660436i, ...
      -0.1101402299076372 + 0.03261803709640307i;
      -0.2008419149861708 + 0.1490117433768365i, ...
    0.310251489958586 - 0.02303091424485477i,  -0.1854601676256336 + 0.181728049653941i, ...
      0.295608864759608 - 0.5647819449422838i;
      0.1198572718465858 - 0.1230966575721692i, ...
    0.004591621691001381 + 0.0004641170819608892i,  0.3304662059847925 - 0.4820711498283887i, ...
      0.06749320884731419 - 0.346397021614317i;
      -0.2688690152234222 - 0.1652086720047534i, ...
    -0.1793802257421584 + 0.05860743927321595i,  0.5235372995239013 + 0.2579777845606151i, ...
      -0.3926950316237215 - 0.1449707406311504i;
      -0.3498536583630073 + 0.09070280031633525i, ...
    -0.08288497844277821 + 0.05058867982255799i,  -0.3202392526929317 - 0.3037968928492991i, ...
      -0.3173573238191287 - 0.3241353242370706i];
c = [];
[dOut, eOut, vtOut, uOut, cOut, info] = f08ms(uplo, d, e, vt, u, c)
 

dOut =

    3.9994
    3.0003
    1.9944
    0.9995


eOut =

     0
     0
     0


vtOut =

  -0.6971 + 0.0000i  -0.0867 - 0.3548i   0.0560 - 0.5400i  -0.1878 - 0.2253i
   0.2403 + 0.0000i   0.0725 - 0.2336i  -0.2477 - 0.5291i   0.7026 + 0.2177i
  -0.5123 + 0.0000i  -0.3030 - 0.1735i   0.0678 + 0.5162i   0.4418 + 0.3864i
  -0.4403 + 0.0000i   0.5294 + 0.6361i  -0.3027 - 0.0346i   0.1667 + 0.0258i


uOut =

  -0.5634 + 0.0016i  -0.2687 - 0.2749i   0.2451 + 0.4657i   0.3787 + 0.2987i
   0.1205 - 0.6108i  -0.2909 + 0.1085i   0.4329 - 0.1758i  -0.0182 - 0.0437i
  -0.0816 + 0.1613i  -0.1660 + 0.3885i  -0.4667 + 0.3821i  -0.0800 - 0.2276i
   0.1441 - 0.1532i   0.1984 - 0.1737i  -0.0034 + 0.1555i   0.2608 - 0.5382i
  -0.2487 - 0.0926i   0.6253 + 0.3304i   0.2643 - 0.0194i   0.1002 + 0.0140i
  -0.3758 + 0.0793i  -0.0307 - 0.0816i   0.1266 + 0.1747i  -0.4175 - 0.4058i


cOut =

     []


info =

                    0



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