naginterfaces.library.lapackeig.dggsvd3

naginterfaces.library.lapackeig.dggsvd3(jobu, jobv, jobq, a, b)

dggsvd3 computes the generalized singular value decomposition (GSVD) of an $m\times n$ real matrix $A$ and a $p\times n$ real matrix $B$.

For full information please refer to the NAG Library document for f08vc

https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/f08/f08vcf.html

Parameters

jobustr, length 1

If $jobu=\text{'U'}$, the orthogonal matrix $U$ is computed.

If $jobu=\text{'N'}$, $U$ is not computed.

jobvstr, length 1

If $jobv=\text{'V'}$, the orthogonal matrix $V$ is computed.

If $jobv=\text{'N'}$, $V$ is not computed.

jobqstr, length 1

If $jobq=\text{'Q'}$, the orthogonal matrix $Q$ is computed.

If $jobq=\text{'N'}$, $Q$ is not computed.

afloat, array-like, shape $(m,n)$

The $m\times n$ matrix $A$.

bfloat, array-like, shape $(p,n)$

The $p\times n$ matrix $B$.

Returns

kint

$k$ and $l$ specify the dimension of the subblocks $k$ and $l$ as described in Notes; $(k+l)$ is the effective numerical rank of (A;B).

lint

$k$ and $l$ specify the dimension of the subblocks $k$ and $l$ as described in Notes; $(k+l)$ is the effective numerical rank of (A;B).

afloat, ndarray, shape $(m,n)$

Contains the triangular matrix $R$, or part of $R$. See Notes for details.

bfloat, ndarray, shape $(p,n)$

Contains the triangular matrix $R$ if $m-k-l<0$. See Notes for details.

alphafloat, ndarray, shape $\left(n\right)$

See the description of $beta$.

betafloat, ndarray, shape $\left(n\right)$

$alpha$ and $beta$ contain the generalized singular value pairs of $A$ and $B$, ${\alpha }_i$ and ${\beta }_i$;

$alpha[0:k]=1$,

$beta[0:k]=0$,

and if $m-k-l\ge 0$,

$alpha[k:k+l]=C$,

$beta[k:k+l]=S$,

or if $m-k-l<0$,

$alpha[k:m]=C$,

$alpha[m:k+l]=0$,

$beta[k:m]=S$,

$beta[m:k+l]=1$, and

$alpha[k+l:n]=0$,

$beta[k+l:n]=0$.

The notation $alpha[k-1:n]$ above refers to consecutive elements $alpha[\text{i}-1]$, for $\text{i}=k,\dots ,n$.

ufloat, ndarray, shape $(:,:)$

If $jobu=\text{'U'}$, $u$ contains the $m\times m$ orthogonal matrix $U$.

If $jobu=\text{'N'}$, $u$ is not referenced.

vfloat, ndarray, shape $(:,:)$

If $jobv=\text{'V'}$, $v$ contains the $p\times p$ orthogonal matrix $V$.

If $jobv=\text{'N'}$, $v$ is not referenced.

qfloat, ndarray, shape $(:,:)$

If $jobq=\text{'Q'}$, $q$ contains the $n\times n$ orthogonal matrix $Q$.

If $jobq=\text{'N'}$, $q$ is not referenced.

iworkint, ndarray, shape $\left(n\right)$

Stores the sorting information. More precisely, if $I$ is the ordered set of indices of $alpha$ containing $C$ (denote as $alpha\left[I\right]$, see $beta$), then the corresponding elements $iwork\left[I\right]-1$ contain the swap pivots, $J$, that sorts $I$ such that $alpha\left[I\right]$ is in descending numerical order.

The following pseudocode sorts the set $I$:

$\begin{array}{c}\text{for}i\in I\\ j=J_i\\ \text{swap}{I}_i\text{and}{I}_j\\ \text{end}\end{array}$

Raises

NagValueError

(errno $-1$)

On entry, error in parameter $jobu$.

Constraint: $jobu=\text{'U'}$ or $\text{'N'}$.

(errno $-2$)

On entry, error in parameter $jobv$.

Constraint: $jobv=\text{'V'}$ or $\text{'N'}$.

(errno $-3$)

On entry, error in parameter $jobq$.

Constraint: $jobq=\text{'Q'}$ or $\text{'N'}$.

(errno $-4$)

On entry, error in parameter $\text{m}$.

Constraint: $m\ge 0$.

(errno $-5$)

On entry, error in parameter $\text{n}$.

Constraint: $n\ge 0$.

(errno $-6$)

On entry, error in parameter $\text{p}$.

Constraint: $p\ge 0$.

(errno $1$)

The Jacobi-type procedure failed to converge.

Notes

Given an $m\times n$ real matrix $A$ and a $p\times n$ real matrix $B$, the generalized singular value decomposition is given by

$$U^{T}AQ=D_1\left(\begin{array}{cc}0& R\end{array}\right)\text{,}{V}^{T}BQ=D_2\left(\begin{array}{cc}0& R\end{array}\right)\text{,}$$

where $U$, $V$ and $Q$ are orthogonal matrices. Let $l$ be the effective numerical rank of $B$ and $(k+l)$ be the effective numerical rank of the matrix (A;B), then the first $k$ generalized singular values are infinite and the remaining $l$ are finite. $R$ is a $(k+l)\times (k+l)$ nonsingular upper triangular matrix, $D_1$ and $D_2$ are $m\times (k+l)$ and $p\times (k+l)$ ‘diagonal’ matrices structured as follows:

if $m-k-l\ge 0$,

$$\begin{array}{c}{D}_1=\begin{array}{c}\\ k\\ l\\ m-k-l\end{array}\left(\begin{array}{cc}k& l\\ I& 0\\ 0& C\\ 0& 0\end{array}\right)\end{array}$$

$$\begin{array}{c}{D}_2=\begin{array}{c}\\ l\\ p-l\end{array}\left(\begin{array}{cc}k& l\\ 0& S\\ 0& 0\end{array}\right)\end{array}$$

$$\begin{array}{c}\left(\begin{array}{cc}0& R\end{array}\right)=\begin{array}{c}\\ k\\ l\end{array}\left(\begin{array}{ccc}n-k-l& k& l\\ 0& R_{11}& R_{12}\\ 0& 0& R_{22}\end{array}\right)\end{array}$$

where

$$C=diag({\alpha }_{k+1},\dots ,{\alpha }_{k+l})\text{,}$$

$$S=diag({\beta }_{k+1},\dots ,{\beta }_{k+l})\text{,}$$

and

$$C^2+S^2=I\text{.}$$

$R$ is stored as a submatrix of $A$ with elements $R_{ij}$ stored as $A_{i,n-k-l+j}$ on exit.

If $m-k-l<0$,

$$\begin{array}{c}{D}_1=\begin{array}{c}\\ k\\ m-k\end{array}\left(\begin{array}{ccc}k& m-k& k+l-m\\ I& 0& 0\\ 0& C& 0\end{array}\right)\end{array}$$

$$\begin{array}{c}{D}_2=\begin{array}{c}\\ m-k\\ k+l-m\\ p-l\end{array}\left(\begin{array}{ccc}k& m-k& k+l-m\\ 0& S& 0\\ 0& 0& I\\ 0& 0& 0\end{array}\right)\end{array}$$

$$\begin{array}{c}\left(\begin{array}{cc}0& R\end{array}\right)=\begin{array}{c}\\ k\\ m-k\\ k+l-m\end{array}\left(\begin{array}{cccc}n-k-l& k& m-k& k+l-m\\ 0& R_{11}& R_{12}& R_{13}\\ 0& 0& R_{22}& R_{23}\\ 0& 0& 0& R_{33}\end{array}\right)\end{array}$$

where

$$C=diag({\alpha }_{k+1},\dots ,{\alpha }_m)\text{,}$$

$$S=diag({\beta }_{k+1},\dots ,{\beta }_m)\text{,}$$

and

$$C^2+S^2=I\text{.}$$

$\left(\begin{array}{ccc}{R}_{11}& R_{12}& R_{13}\\ 0& R_{22}& R_{23}\end{array}\right)$ is stored as a submatrix of $A$ with $R_{ij}$ stored as $A_{i,n-k-l+j}$, and $R_{33}$ is stored as a submatrix of $B$ with ${\left(R_{33}\right)}_{ij}$ stored as $B_{m-k+i,n+m-k-l+j}$.

The function computes $C$, $S$, $R$ and, optionally, the orthogonal transformation matrices $U$, $V$ and $Q$.

In particular, if $B$ is an $n\times n$ nonsingular matrix, then the GSVD of $A$ and $B$ implicitly gives the SVD of $A{B}^{-1}$:

$$A{B}^{-1}=U\left(D_1{D}_2^{-1}\right)V^T\text{.}$$

If (A;B) has orthonormal columns, then the GSVD of $A$ and $B$ is also equal to the CS decomposition of $A$ and $B$. Furthermore, the GSVD can be used to derive the solution of the eigenvalue problem:

$$A^{T}Ax=\lambda B^{T}Bx\text{.}$$

In some literature, the GSVD of $A$ and $B$ is presented in the form

$$U^{T}AX=\left(\begin{array}{cc}0& D_1\end{array}\right)\text{,}{V}^{T}BX=\left(\begin{array}{cc}0& D_2\end{array}\right)\text{,}$$

where $U$ and $V$ are orthogonal and $X$ is nonsingular, and $D_1$ and $D_2$ are ‘diagonal’. The former GSVD form can be converted to the latter form by setting

$$\begin{array}{c}X=Q\left(\begin{array}{cc}I& 0\\ 0& R^{-1}\end{array}\right)\text{.}\end{array}$$

A two stage process is used to compute the GSVD of the matrix pair $(A,B)$. The pair is first reduced to upper triangular form by orthogonal transformations using dggsvp3(). The GSVD of the resulting upper triangular matrix pair is then performed by dtgsja() which uses a variant of the Kogbetliantz algorithm (a cyclic Jacobi method).

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, https://www.netlib.org/lapack/lug

Golub, G H and Van Loan, C F, 2012, Matrix Computations, (4th Edition), Johns Hopkins University Press, Baltimore