Linear Algebra¶
This chapter describes functions for solving linear systems. The
library provides linear algebra operations which operate directly on
the gsl_vector
and gsl_matrix
objects. These routines
use the standard algorithms from Golub & Van Loan’s Matrix
Computations with Level1 and Level2 BLAS calls for efficiency.
The functions described in this chapter are declared in the header file
gsl_linalg.h
.
LU Decomposition¶
A general by matrix has an decomposition
where is an by permutation matrix, is by and is by. For square matrices, is a lower unit triangular matrix and is upper triangular. For , is a unit lower trapezoidal matrix of size by. For , is upper trapezoidal of size by. For square matrices this decomposition can be used to convert the linear system into a pair of triangular systems (, ), which can be solved by forward and backsubstitution. Note that the decomposition is valid for singular matrices.

int
gsl_linalg_LU_decomp
(gsl_matrix * A, gsl_permutation * p, int * signum)¶ 
int
gsl_linalg_complex_LU_decomp
(gsl_matrix_complex * A, gsl_permutation * p, int * signum)¶ These functions factorize the matrix
A
into the decomposition . On output the diagonal and upper triangular (or trapezoidal) part of the input matrixA
contain the matrix . The lower triangular (or trapezoidal) part of the input matrix (excluding the diagonal) contains . The diagonal elements of are unity, and are not stored.The permutation matrix is encoded in the permutation
p
on output. The th column of the matrix is given by the th column of the identity matrix, where the th element of the permutation vector. The sign of the permutation is given bysignum
. It has the value , where is the number of interchanges in the permutation.The algorithm used in the decomposition is Gaussian Elimination with partial pivoting (Golub & Van Loan, Matrix Computations, Algorithm 3.4.1), combined with a recursive algorithm based on Level 3 BLAS (Peise and Bientinesi, 2016).

int
gsl_linalg_LU_solve
(const gsl_matrix * LU, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x)¶ 
int
gsl_linalg_complex_LU_solve
(const gsl_matrix_complex * LU, const gsl_permutation * p, const gsl_vector_complex * b, gsl_vector_complex * x)¶ These functions solve the square system using the decomposition of into (
LU
,p
) given bygsl_linalg_LU_decomp()
orgsl_linalg_complex_LU_decomp()
as input.

int
gsl_linalg_LU_svx
(const gsl_matrix * LU, const gsl_permutation * p, gsl_vector * x)¶ 
int
gsl_linalg_complex_LU_svx
(const gsl_matrix_complex * LU, const gsl_permutation * p, gsl_vector_complex * x)¶ These functions solve the square system inplace using the precomputed decomposition of into (
LU
,p
). On inputx
should contain the righthand side , which is replaced by the solution on output.

int
gsl_linalg_LU_refine
(const gsl_matrix * A, const gsl_matrix * LU, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x, gsl_vector * work)¶ 
int
gsl_linalg_complex_LU_refine
(const gsl_matrix_complex * A, const gsl_matrix_complex * LU, const gsl_permutation * p, const gsl_vector_complex * b, gsl_vector_complex * x, gsl_vector_complex * work)¶ These functions apply an iterative improvement to
x
, the solution of , from the precomputed decomposition of into (LU
,p
). Additional workspace of lengthN
is required inwork
.

int
gsl_linalg_LU_invert
(const gsl_matrix * LU, const gsl_permutation * p, gsl_matrix * inverse)¶ 
int
gsl_linalg_complex_LU_invert
(const gsl_matrix_complex * LU, const gsl_permutation * p, gsl_matrix_complex * inverse)¶ These functions compute the inverse of a matrix from its decomposition (
LU
,p
), storing the result in the matrixinverse
. The inverse is computed by computing the inverses , and finally forming the product . Each step is based on Level 3 BLAS calls.It is preferable to avoid direct use of the inverse whenever possible, as the linear solver functions can obtain the same result more efficiently and reliably (consult any introductory textbook on numerical linear algebra for details).

int
gsl_linalg_LU_invx
(gsl_matrix * LU, const gsl_permutation * p)¶ 
int
gsl_linalg_complex_LU_invx
(gsl_matrix_complex * LU, const gsl_permutation * p)¶ These functions compute the inverse of a matrix from its decomposition (
LU
,p
), storing the result inplace in the matrixLU
. The inverse is computed by computing the inverses , and finally forming the product . Each step is based on Level 3 BLAS calls.It is preferable to avoid direct use of the inverse whenever possible, as the linear solver functions can obtain the same result more efficiently and reliably (consult any introductory textbook on numerical linear algebra for details).

double
gsl_linalg_LU_det
(gsl_matrix * LU, int signum)¶ 
gsl_complex
gsl_linalg_complex_LU_det
(gsl_matrix_complex * LU, int signum)¶ These functions compute the determinant of a matrix from its decomposition,
LU
. The determinant is computed as the product of the diagonal elements of and the sign of the row permutationsignum
.

double
gsl_linalg_LU_lndet
(gsl_matrix * LU)¶ 
double
gsl_linalg_complex_LU_lndet
(gsl_matrix_complex * LU)¶ These functions compute the logarithm of the absolute value of the determinant of a matrix , , from its decomposition,
LU
. This function may be useful if the direct computation of the determinant would overflow or underflow.

int
gsl_linalg_LU_sgndet
(gsl_matrix * LU, int signum)¶ 
gsl_complex
gsl_linalg_complex_LU_sgndet
(gsl_matrix_complex * LU, int signum)¶ These functions compute the sign or phase factor of the determinant of a matrix , , from its decomposition,
LU
.
QR Decomposition¶
A general rectangular by matrix has a decomposition into the product of an orthogonal by square matrix (where ) and an by righttriangular matrix ,
This decomposition can be used to convert the square linear system into the triangular system , which can be solved by backsubstitution. Another use of the decomposition is to compute an orthonormal basis for a set of vectors. The first columns of form an orthonormal basis for the range of , , when has full column rank.
When , the bottom rows of are zero, and so can be naturally partioned as
where is by upper triangular, is by, and is by. is sometimes called the thin or reduced QR decomposition. The solution of the least squares problem when has full rank is:
where is the first elements of . If is rank deficient, see QR Decomposition with Column Pivoting and Complete Orthogonal Decomposition.
GSL offers two interfaces for the decomposition. The first proceeds by zeroing out columns below the diagonal of , one column at a time using Householder transforms. In this method, the factor is represented as a product of Householder reflectors:
where each for a scalar and column vector . In this method, functions which compute the full matrix or apply to a right hand side vector operate by applying the Householder matrices one at a time using Level 2 BLAS.
The second interface is based on a Level 3 BLAS block recursive algorithm developed by Elmroth and Gustavson. In this case, is written in block form as
where is an by matrix of the column vectors
and is an by upper triangular matrix, whose diagonal elements
are the . Computing the full , while requiring more flops than
the Level 2 approach, offers the advantage that all standard operations can take advantage
of cache efficient Level 3 BLAS operations, and so this method often performs faster
than the Level 2 approach. The functions for the recursive block algorithm have a
_r
suffix, and it is recommended to use this interface for performance
critical applications.

int
gsl_linalg_QR_decomp
(gsl_matrix * A, gsl_vector * tau)¶ This function factorizes the by matrix
A
into the decomposition . On output the diagonal and upper triangular part of the input matrix contain the matrix . The vectortau
and the columns of the lower triangular part of the matrixA
contain the Householder coefficients and Householder vectors which encode the orthogonal matrixQ
. The vectortau
must be of length . The matrix is related to these components by, where and is the Householder vector . This is the same storage scheme as used by LAPACK.The algorithm used to perform the decomposition is Householder QR (Golub & Van Loan, “Matrix Computations”, Algorithm 5.2.1). It is recommended to use this function only for small to modest sized matrices, or if .

int
gsl_linalg_QR_decomp_r
(gsl_matrix * A, gsl_matrix * T)¶ This function factors the by matrix
A
into the decomposition using the recursive Level 3 BLAS algorithm of Elmroth and Gustavson. On output the diagonal and upper triangular part ofA
contain the matrix . The by matrixT
stores the upper triangular factor appearing in . The matrix is given by , where the elements below the diagonal ofA
contain the columns of on output.This algorithm requires and performs best for “tallskinny” matrices, i.e. .

int
gsl_linalg_QR_solve
(const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x)¶ 
int
gsl_linalg_QR_solve_r
(const gsl_matrix * QR, const gsl_matrix * T, const gsl_vector * b, gsl_vector * x)¶ These functions solve the square system using the decomposition of held in (
QR
,tau
) or (QR
,T
) which must have been computed previously withgsl_linalg_QR_decomp()
orgsl_linalg_QR_decomp_r()
. The leastsquares solution for rectangular systems can be found usinggsl_linalg_QR_lssolve()
orgsl_linalg_QR_lssolve_r()
.

int
gsl_linalg_QR_svx
(const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * x)¶ This function solves the square system inplace using the decomposition of held in (
QR
,tau
) which must have been computed previously bygsl_linalg_QR_decomp()
. On inputx
should contain the righthand side , which is replaced by the solution on output.

int
gsl_linalg_QR_lssolve
(const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x, gsl_vector * residual)¶ This function finds the least squares solution to the overdetermined system where the matrix
A
has more rows than columns. The least squares solution minimizes the Euclidean norm of the residual, .The routine requires as input the decomposition of into (QR
,tau
) given bygsl_linalg_QR_decomp()
. The solution is returned inx
. The residual is computed as a byproduct and stored inresidual
.

int
gsl_linalg_QR_lssolve_r
(const gsl_matrix * QR, const gsl_matrix * T, const gsl_vector * b, gsl_vector * x, gsl_vector * work)¶ This function finds the least squares solution to the overdetermined system where the matrix
A
has more rows than columns. The least squares solution minimizes the Euclidean norm of the residual, .The routine requires as input the decomposition of into (QR
,T
) given bygsl_linalg_QR_decomp_r()
. The parameterx
is of length . The solution is returned in the first rows ofx
, i.e.x[0], x[1], ..., x[N1]
. The last rows ofx
contain a vector whose norm is equal to the residual norm . This similar to the behavior of LAPACK DGELS. Additional workspace of length is required inwork
.

int
gsl_linalg_QR_QTvec
(const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * v)¶ 
int
gsl_linalg_QR_QTvec_r
(const gsl_matrix * QR, const gsl_matrix * T, gsl_vector * v, gsl_vector * work)¶ These functions apply the matrix encoded in the decomposition (
QR
,tau
) or (QR
,T
) to the vectorv
, storing the result inv
. The matrix multiplication is carried out directly using the encoding of the Householder vectors without needing to form the full matrix .The recursive variant
QTvec_r
requires additional workspace of size inwork
.

int
gsl_linalg_QR_Qvec
(const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * v)¶ This function applies the matrix encoded in the decomposition (
QR
,tau
) to the vectorv
, storing the result inv
. The matrix multiplication is carried out directly using the encoding of the Householder vectors without needing to form the full matrix .

int
gsl_linalg_QR_QTmat
(const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * B)¶ 
int
gsl_linalg_QR_QTmat_r
(const gsl_matrix * QR, const gsl_matrix * T, gsl_matrix * B, gsl_matrix * work)¶ This function applies the matrix encoded in the decomposition (
QR
,tau
) or (QR
,T
) to the by matrixB
, storing the result inB
. The matrix multiplication is carried out directly using the encoding of the Householder vectors without needing to form the full matrix .The recursive variant
QTmat_r
requires additional workspace of size by inwork
.

int
gsl_linalg_QR_Rsolve
(const gsl_matrix * QR, const gsl_vector * b, gsl_vector * x)¶ This function solves the triangular system for
x
. It may be useful if the product has already been computed usinggsl_linalg_QR_QTvec()
.

int
gsl_linalg_QR_Rsvx
(const gsl_matrix * QR, gsl_vector * x)¶ This function solves the triangular system for
x
inplace. On inputx
should contain the righthand side and is replaced by the solution on output. This function may be useful if the product has already been computed usinggsl_linalg_QR_QTvec()
.

int
gsl_linalg_QR_unpack
(const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * Q, gsl_matrix * R)¶ This function unpacks the encoded decomposition (
QR
,tau
) into the matricesQ
andR
, whereQ
is by andR
is by.

int
gsl_linalg_QR_unpack_r
(const gsl_matrix * QR, const gsl_matrix * T, gsl_matrix * Q, gsl_matrix * R)¶ This function unpacks the encoded decomposition (
QR
,T
) as output fromgsl_linalg_QR_decomp_r()
into the matricesQ
andR
, whereQ
is by andR
is by. Note that the full matrix is by, however the lower trapezoidal portion is zero, so only the upper triangular factor is stored.

int
gsl_linalg_QR_QRsolve
(gsl_matrix * Q, gsl_matrix * R, const gsl_vector * b, gsl_vector * x)¶ This function solves the system for
x
. It can be used when the decomposition of a matrix is available in unpacked form as (Q
,R
).

int
gsl_linalg_QR_update
(gsl_matrix * Q, gsl_matrix * R, gsl_vector * w, const gsl_vector * v)¶ This function performs a rank1 update of the decomposition (
Q
,R
). The update is given by where the output matrices and are also orthogonal and right triangular. Note thatw
is destroyed by the update.

int
gsl_linalg_R_solve
(const gsl_matrix * R, const gsl_vector * b, gsl_vector * x)¶ This function solves the triangular system for the by matrix
R
.

int
gsl_linalg_R_svx
(const gsl_matrix * R, gsl_vector * x)¶ This function solves the triangular system inplace. On input
x
should contain the righthand side , which is replaced by the solution on output.

int
gsl_linalg_QR_rcond
(const gsl_matrix * QR, double * rcond, gsl_vector * work)¶ This function estimates the reciprocal condition number (using the 1norm) of the factor, stored in the upper triangle of
QR
. The reciprocal condition number estimate, defined as , is stored inrcond
. Additional workspace of size is required inwork
.
Triangle on Top of Rectangle¶
This section provides routines for computing the decomposition of the specialized matrix
where is an by upper triangular matrix, and is an by dense matrix. This type of matrix arises, for example, in the sequential TSQR algorithm. The Elmroth and Gustavson algorithm is used to efficiently factor this matrix. Due to the upper triangular factor, the matrix takes the form
with
and is dense and of the same dimensions as .

int
gsl_linalg_QR_TR_decomp
(gsl_matrix * S, gsl_matrix * A, gsl_matrix * T)¶ This function computes the decomposition of the matrix , where is by upper triangular and is by dense. On output, is replaced by the factor, and is replaced by .
QR Decomposition with Column Pivoting¶
The decomposition of an by matrix can be extended to the rank deficient case by introducing a column permutation ,
The first columns of form an orthonormal basis for the range of for a matrix with column rank . This decomposition can also be used to convert the square linear system into the triangular system , which can be solved by backsubstitution and permutation. We denote the decomposition with column pivoting by since . When is rank deficient with , the matrix can be partitioned as
where is by and nonsingular. In this case, a basic least squares solution for the overdetermined system can be obtained as
where consists of the first elements of . This basic solution is not guaranteed to be the minimum norm solution unless (see Complete Orthogonal Decomposition).

int
gsl_linalg_QRPT_decomp
(gsl_matrix * A, gsl_vector * tau, gsl_permutation * p, int * signum, gsl_vector * norm)¶ This function factorizes the by matrix
A
into the decomposition . On output the diagonal and upper triangular part of the input matrix contain the matrix . The permutation matrix is stored in the permutationp
. The sign of the permutation is given bysignum
. It has the value , where is the number of interchanges in the permutation. The vectortau
and the columns of the lower triangular part of the matrixA
contain the Householder coefficients and vectors which encode the orthogonal matrixQ
. The vectortau
must be of length . The matrix is related to these components by, where and is the Householder vectorThis is the same storage scheme as used by LAPACK. The vector
norm
is a workspace of lengthN
used for column pivoting.The algorithm used to perform the decomposition is Householder QR with column pivoting (Golub & Van Loan, “Matrix Computations”, Algorithm 5.4.1).

int
gsl_linalg_QRPT_decomp2
(const gsl_matrix * A, gsl_matrix * q, gsl_matrix * r, gsl_vector * tau, gsl_permutation * p, int * signum, gsl_vector * norm)¶ This function factorizes the matrix
A
into the decomposition without modifyingA
itself and storing the output in the separate matricesq
andr
.

int
gsl_linalg_QRPT_solve
(const gsl_matrix * QR, const gsl_vector * tau, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x)¶ This function solves the square system using the decomposition of held in (
QR
,tau
,p
) which must have been computed previously bygsl_linalg_QRPT_decomp()
.

int
gsl_linalg_QRPT_svx
(const gsl_matrix * QR, const gsl_vector * tau, const gsl_permutation * p, gsl_vector * x)¶ This function solves the square system inplace using the decomposition of held in (
QR
,tau
,p
). On inputx
should contain the righthand side , which is replaced by the solution on output.

int
gsl_linalg_QRPT_lssolve
(const gsl_matrix * QR, const gsl_vector * tau, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x, gsl_vector * residual)¶ This function finds the least squares solution to the overdetermined system where the matrix
A
has more rows than columns and is assumed to have full rank. The least squares solution minimizes the Euclidean norm of the residual, . The routine requires as input the decomposition of into (QR
,tau
,p
) given bygsl_linalg_QRPT_decomp()
. The solution is returned inx
. The residual is computed as a byproduct and stored inresidual
. For rank deficient matrices,gsl_linalg_QRPT_lssolve2()
should be used instead.

int
gsl_linalg_QRPT_lssolve2
(const gsl_matrix * QR, const gsl_vector * tau, const gsl_permutation * p, const gsl_vector * b, const size_t rank, gsl_vector * x, gsl_vector * residual)¶ This function finds the least squares solution to the overdetermined system where the matrix
A
has more rows than columns and has rank given by the inputrank
. If the user does not know the rank of , the routinegsl_linalg_QRPT_rank()
can be called to estimate it. The least squares solution is the socalled “basic” solution discussed above and may not be the minimum norm solution. The routine requires as input the decomposition of into (QR
,tau
,p
) given bygsl_linalg_QRPT_decomp()
. The solution is returned inx
. The residual is computed as a byproduct and stored inresidual
.

int
gsl_linalg_QRPT_QRsolve
(const gsl_matrix * Q, const gsl_matrix * R, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x)¶ This function solves the square system for
x
. It can be used when the decomposition of a matrix is available in unpacked form as (Q
,R
).

int
gsl_linalg_QRPT_update
(gsl_matrix * Q, gsl_matrix * R, const gsl_permutation * p, gsl_vector * w, const gsl_vector * v)¶ This function performs a rank1 update of the decomposition (
Q
,R
,p
). The update is given by where the output matrices and are also orthogonal and right triangular. Note thatw
is destroyed by the update. The permutationp
is not changed.

int
gsl_linalg_QRPT_Rsolve
(const gsl_matrix * QR, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x)¶ This function solves the triangular system for the by matrix contained in
QR
.

int
gsl_linalg_QRPT_Rsvx
(const gsl_matrix * QR, const gsl_permutation * p, gsl_vector * x)¶ This function solves the triangular system inplace for the by matrix contained in
QR
. On inputx
should contain the righthand side , which is replaced by the solution on output.

size_t
gsl_linalg_QRPT_rank
(const gsl_matrix * QR, const double tol)¶ This function estimates the rank of the triangular matrix contained in
QR
. The algorithm simply counts the number of diagonal elements of whose absolute value is greater than the specified tolerancetol
. If the inputtol
is negative, a default value of is used.

int
gsl_linalg_QRPT_rcond
(const gsl_matrix * QR, double * rcond, gsl_vector * work)¶ This function estimates the reciprocal condition number (using the 1norm) of the factor, stored in the upper triangle of
QR
. The reciprocal condition number estimate, defined as , is stored inrcond
. Additional workspace of size is required inwork
.
LQ Decomposition¶
A general rectangular by matrix has a decomposition into the product of a lower trapezoidal by matrix and an orthogonal by square matrix :
If , then can be written as where is by lower triangular, and
where consists of the first rows of , and contains the remaining rows. The factorization of is essentially the same as the QR factorization of .
The factorization may be used to find the minimum norm solution of an underdetermined system of equations , where is by and . The solution is

int
gsl_linalg_LQ_decomp
(gsl_matrix * A, gsl_vector * tau)¶ This function factorizes the by matrix
A
into the decomposition . On output the diagonal and lower trapezoidal part of the input matrix contain the matrix . The vectortau
and the elements above the diagonal of the matrixA
contain the Householder coefficients and Householder vectors which encode the orthogonal matrixQ
. The vectortau
must be of length . The matrix is related to these components by, where and is the Householder vector . This is the same storage scheme as used by LAPACK.

int
gsl_linalg_LQ_lssolve
(const gsl_matrix * LQ, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x, gsl_vector * residual)¶ This function finds the minimum norm least squares solution to the underdetermined system , where the by matrix
A
has . The routine requires as input the decomposition of into (LQ
,tau
) given bygsl_linalg_LQ_decomp()
. The solution is returned inx
. The residual, , is computed as a byproduct and stored inresidual
.

int
gsl_linalg_LQ_unpack
(const gsl_matrix * LQ, const gsl_vector * tau, gsl_matrix * Q, gsl_matrix * L)¶ This function unpacks the encoded decomposition (
LQ
,tau
) into the matricesQ
andL
, whereQ
is by andL
is by.

int
gsl_linalg_LQ_QTvec
(const gsl_matrix * LQ, const gsl_vector * tau, gsl_vector * v)¶ This function applies to the vector
v
, storing the result inv
on output.
Complete Orthogonal Decomposition¶
The complete orthogonal decomposition of a by matrix is a generalization of the QR decomposition with column pivoting, given by
where is a by permutation matrix, is by orthogonal, is by upper triangular, with , and is by orthogonal. If has full rank, then , and this reduces to the QR decomposition with column pivoting.
For a rank deficient least squares problem, , the solution vector is not unique. However if we further require that is minimized, then the complete orthogonal decomposition gives us the ability to compute the unique minimum norm solution, which is given by
and the vector is the first elements of .
The COD also enables a straightforward solution of regularized least squares problems in Tikhonov standard form, written as
where is a regularization parameter which represents a tradeoff between minimizing the residual norm and the solution norm . For this system, the solution is given by
where is a vector of length which is found by solving
and is defined above. The equation above can be solved efficiently for different values of using QR factorizations of the left hand side matrix.

int
gsl_linalg_COD_decomp
(gsl_matrix * A, gsl_vector * tau_Q, gsl_vector * tau_Z, gsl_permutation * p, size_t * rank, gsl_vector * work)¶ 
int
gsl_linalg_COD_decomp_e
(gsl_matrix * A, gsl_vector * tau_Q, gsl_vector * tau_Z, gsl_permutation * p, double tol, size_t * rank, gsl_vector * work)¶ These functions factor the by matrix
A
into the decomposition . The rank ofA
is computed as the number of diagonal elements of greater than the tolerancetol
and output inrank
. Iftol
is not specified, a default value is used (seegsl_linalg_QRPT_rank()
). On output, the permutation matrix is stored inp
. The matrix is stored in the upperrank
byrank
block ofA
. The matrices and are encoded in packed storage inA
on output. The vectorstau_Q
andtau_Z
contain the Householder scalars corresponding to the matrices and respectively and must be of length . The vectorwork
is additional workspace of length .

int
gsl_linalg_COD_lssolve
(const gsl_matrix * QRZT, const gsl_vector * tau_Q, const gsl_vector * tau_Z, const gsl_permutation * p, const size_t rank, const gsl_vector * b, gsl_vector * x, gsl_vector * residual)¶ This function finds the unique minimum norm least squares solution to the overdetermined system where the matrix
A
has more rows than columns. The least squares solution minimizes the Euclidean norm of the residual, as well as the norm of the solution . The routine requires as input the decomposition of into (QRZT
,tau_Q
,tau_Z
,p
,rank
) given bygsl_linalg_COD_decomp()
. The solution is returned inx
. The residual, , is computed as a byproduct and stored inresidual
.

int
gsl_linalg_COD_lssolve2
(const double lambda, const gsl_matrix * QRZT, const gsl_vector * tau_Q, const gsl_vector * tau_Z, const gsl_permutation * p, const size_t rank, const gsl_vector * b, gsl_vector * x, gsl_vector * residual, gsl_matrix * S, gsl_vector * work)¶ This function finds the solution to the regularized least squares problem in Tikhonov standard form, . The routine requires as input the decomposition of into (
QRZT
,tau_Q
,tau_Z
,p
,rank
) given bygsl_linalg_COD_decomp()
. The parameter is supplied inlambda
. The solution is returned inx
. The residual, , is stored inresidual
on output.S
is additional workspace of sizerank
byrank
.work
is additional workspace of lengthrank
.

int
gsl_linalg_COD_unpack
(const gsl_matrix * QRZT, const gsl_vector * tau_Q, const gsl_vector * tau_Z, const size_t rank, gsl_matrix * Q, gsl_matrix * R, gsl_matrix * Z)¶ This function unpacks the encoded decomposition (
QRZT
,tau_Q
,tau_Z
,rank
) into the matricesQ
,R
, andZ
, whereQ
is by,R
is by, andZ
is by.

int
gsl_linalg_COD_matZ
(const gsl_matrix * QRZT, const gsl_vector * tau_Z, const size_t rank, gsl_matrix * A, gsl_vector * work)¶ This function multiplies the input matrix
A
on the right byZ
, using the encoded decomposition (QRZT
,tau_Z
,rank
).A
must have columns but may have any number of rows. Additional workspace of length is provided inwork
.
Singular Value Decomposition¶
A general rectangular by matrix has a singular value decomposition (SVD) into the product of an by orthogonal matrix , an by diagonal matrix of singular values and the transpose of an by orthogonal square matrix ,
The singular values are all nonnegative and are generally chosen to form a nonincreasing sequence
The singular value decomposition of a matrix has many practical uses. The condition number of the matrix is given by the ratio of the largest singular value to the smallest singular value. The presence of a zero singular value indicates that the matrix is singular. The number of nonzero singular values indicates the rank of the matrix. In practice singular value decomposition of a rankdeficient matrix will not produce exact zeroes for singular values, due to finite numerical precision. Small singular values should be edited by choosing a suitable tolerance.
For a rankdeficient matrix, the null space of is given by the columns of corresponding to the zero singular values. Similarly, the range of is given by columns of corresponding to the nonzero singular values.
Note that the routines here compute the “thin” version of the SVD with as by orthogonal matrix. This allows inplace computation and is the most commonlyused form in practice. Mathematically, the “full” SVD is defined with as an by orthogonal matrix and as an by diagonal matrix (with additional rows of zeros).

int
gsl_linalg_SV_decomp
(gsl_matrix * A, gsl_matrix * V, gsl_vector * S, gsl_vector * work)¶ This function factorizes the by matrix
A
into the singular value decomposition for . On output the matrixA
is replaced by . The diagonal elements of the singular value matrix are stored in the vectorS
. The singular values are nonnegative and form a nonincreasing sequence from to . The matrixV
contains the elements of in untransposed form. To form the product it is necessary to take the transpose ofV
. A workspace of lengthN
is required inwork
.This routine uses the GolubReinsch SVD algorithm.

int
gsl_linalg_SV_decomp_mod
(gsl_matrix * A, gsl_matrix * X, gsl_matrix * V, gsl_vector * S, gsl_vector * work)¶ This function computes the SVD using the modified GolubReinsch algorithm, which is faster for . It requires the vector
work
of lengthN
and the by matrixX
as additional working space.

int
gsl_linalg_SV_decomp_jacobi
(gsl_matrix * A, gsl_matrix * V, gsl_vector * S)¶ This function computes the SVD of the by matrix
A
using onesided Jacobi orthogonalization for . The Jacobi method can compute singular values to higher relative accuracy than GolubReinsch algorithms (see references for details).

int
gsl_linalg_SV_solve
(const gsl_matrix * U, const gsl_matrix * V, const gsl_vector * S, const gsl_vector * b, gsl_vector * x)¶ This function solves the system using the singular value decomposition (
U
,S
,V
) of which must have been computed previously withgsl_linalg_SV_decomp()
.Only nonzero singular values are used in computing the solution. The parts of the solution corresponding to singular values of zero are ignored. Other singular values can be edited out by setting them to zero before calling this function.
In the overdetermined case where
A
has more rows than columns the system is solved in the least squares sense, returning the solutionx
which minimizes .

int
gsl_linalg_SV_leverage
(const gsl_matrix * U, gsl_vector * h)¶ This function computes the statistical leverage values of a matrix using its singular value decomposition (
U
,S
,V
) previously computed withgsl_linalg_SV_decomp()
. are the diagonal values of the matrix and depend only on the matrixU
which is the input to this function.
Cholesky Decomposition¶
A symmetric, positive definite square matrix has a Cholesky decomposition into a product of a lower triangular matrix and its transpose ,
This is sometimes referred to as taking the squareroot of a matrix. The Cholesky decomposition can only be carried out when all the eigenvalues of the matrix are positive. This decomposition can be used to convert the linear system into a pair of triangular systems (, ), which can be solved by forward and backsubstitution.
If the matrix is near singular, it is sometimes possible to reduce the condition number and recover a more accurate solution vector by scaling as
where is a diagonal matrix whose elements are given by . This scaling is also known as Jacobi preconditioning. There are routines below to solve both the scaled and unscaled systems.

int
gsl_linalg_cholesky_decomp1
(gsl_matrix * A)¶ 
int
gsl_linalg_complex_cholesky_decomp
(gsl_matrix_complex * A)¶ These functions factorize the symmetric, positivedefinite square matrix
A
into the Cholesky decomposition (or for the complex case). On input, the values from the diagonal and lowertriangular part of the matrixA
are used (the upper triangular part is ignored). On output the diagonal and lower triangular part of the input matrixA
contain the matrix , while the upper triangular part contains the original matrix. If the matrix is not positivedefinite then the decomposition will fail, returning the error codeGSL_EDOM
.When testing whether a matrix is positivedefinite, disable the error handler first to avoid triggering an error. These functions use Level 3 BLAS to compute the Cholesky factorization (Peise and Bientinesi, 2016).

int
gsl_linalg_cholesky_decomp
(gsl_matrix * A)¶ This function is now deprecated and is provided only for backward compatibility.

int
gsl_linalg_cholesky_solve
(const gsl_matrix * cholesky, const gsl_vector * b, gsl_vector * x)¶ 
int
gsl_linalg_complex_cholesky_solve
(const gsl_matrix_complex * cholesky, const gsl_vector_complex * b, gsl_vector_complex * x)¶ These functions solve the system using the Cholesky decomposition of held in the matrix
cholesky
which must have been previously computed bygsl_linalg_cholesky_decomp()
orgsl_linalg_complex_cholesky_decomp()
.

int
gsl_linalg_cholesky_svx
(const gsl_matrix * cholesky, gsl_vector * x)¶ 
int
gsl_linalg_complex_cholesky_svx
(const gsl_matrix_complex * cholesky, gsl_vector_complex * x)¶ These functions solve the system inplace using the Cholesky decomposition of held in the matrix
cholesky
which must have been previously computed bygsl_linalg_cholesky_decomp()
orgsl_linalg_complex_cholesky_decomp()
. On inputx
should contain the righthand side , which is replaced by the solution on output.

int
gsl_linalg_cholesky_invert
(gsl_matrix * cholesky)¶ 
int
gsl_linalg_complex_cholesky_invert
(gsl_matrix_complex * cholesky)¶ These functions compute the inverse of a matrix from its Cholesky decomposition
cholesky
, which must have been previously computed bygsl_linalg_cholesky_decomp()
orgsl_linalg_complex_cholesky_decomp()
. On output, the inverse is stored inplace incholesky
.

int
gsl_linalg_cholesky_decomp2
(gsl_matrix * A, gsl_vector * S)¶ This function calculates a diagonal scaling transformation for the symmetric, positivedefinite square matrix
A
, and then computes the Cholesky decomposition . On input, the values from the diagonal and lowertriangular part of the matrixA
are used (the upper triangular part is ignored). On output the diagonal and lower triangular part of the input matrixA
contain the matrix , while the upper triangular part of the input matrix is overwritten with (the diagonal terms being identical for both and ). If the matrix is not positivedefinite then the decomposition will fail, returning the error codeGSL_EDOM
. The diagonal scale factors are stored inS
on output.When testing whether a matrix is positivedefinite, disable the error handler first to avoid triggering an error.

int
gsl_linalg_cholesky_solve2
(const gsl_matrix * cholesky, const gsl_vector * S, const gsl_vector * b, gsl_vector * x)¶ This function solves the system using the Cholesky decomposition of held in the matrix
cholesky
which must have been previously computed bygsl_linalg_cholesky_decomp2()
.

int
gsl_linalg_cholesky_svx2
(const gsl_matrix * cholesky, const gsl_vector * S, gsl_vector * x)¶ This function solves the system inplace using the Cholesky decomposition of held in the matrix
cholesky
which must have been previously computed bygsl_linalg_cholesky_decomp2()
. On inputx
should contain the righthand side , which is replaced by the solution on output.

int
gsl_linalg_cholesky_scale
(const gsl_matrix * A, gsl_vector * S)¶ This function calculates a diagonal scaling transformation of the symmetric, positive definite matrix
A
, such that has a condition number within a factor of of the matrix of smallest possible condition number over all possible diagonal scalings. On output,S
contains the scale factors, given by . For any , the corresponding scale factor is set to .

int
gsl_linalg_cholesky_scale_apply
(gsl_matrix * A, const gsl_vector * S)¶ This function applies the scaling transformation
S
to the matrixA
. On output,A
is replaced by .

int
gsl_linalg_cholesky_rcond
(const gsl_matrix * cholesky, double * rcond, gsl_vector * work)¶ This function estimates the reciprocal condition number (using the 1norm) of the symmetric positive definite matrix , using its Cholesky decomposition provided in
cholesky
. The reciprocal condition number estimate, defined as , is stored inrcond
. Additional workspace of size is required inwork
.
Pivoted Cholesky Decomposition¶
A symmetric positive semidefinite square matrix has an alternate Cholesky decomposition into a product of a lower unit triangular matrix , a diagonal matrix and , given by . For postive definite matrices, this is equivalent to the Cholesky formulation discussed above, with the standard Cholesky lower triangular factor given by . For illconditioned matrices, it can help to use a pivoting strategy to prevent the entries of and from growing too large, and also ensure , where are the diagonal entries of . The final decomposition is given by
where is a permutation matrix.

int
gsl_linalg_pcholesky_decomp
(gsl_matrix * A, gsl_permutation * p)¶ This function factors the symmetric, positivedefinite square matrix
A
into the Pivoted Cholesky decomposition . On input, the values from the diagonal and lowertriangular part of the matrixA
are used to construct the factorization. On output the diagonal of the input matrixA
stores the diagonal elements of , and the lower triangular portion ofA
contains the matrix . Since has ones on its diagonal these do not need to be explicitely stored. The upper triangular portion ofA
is unmodified. The permutation matrix is stored inp
on output.

int
gsl_linalg_pcholesky_solve
(const gsl_matrix * LDLT, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x)¶ This function solves the system using the Pivoted Cholesky decomposition of held in the matrix
LDLT
and permutationp
which must have been previously computed bygsl_linalg_pcholesky_decomp()
.

int
gsl_linalg_pcholesky_svx
(const gsl_matrix * LDLT, const gsl_permutation * p, gsl_vector * x)¶ This function solves the system inplace using the Pivoted Cholesky decomposition of held in the matrix
LDLT
and permutationp
which must have been previously computed bygsl_linalg_pcholesky_decomp()
. On input,x
contains the right hand side vector which is replaced by the solution vector on output.

int
gsl_linalg_pcholesky_decomp2
(gsl_matrix * A, gsl_permutation * p, gsl_vector * S)¶ This function computes the pivoted Cholesky factorization of the matrix , where the input matrix
A
is symmetric and positive definite, and the diagonal scaling matrixS
is computed to reduce the condition number ofA
as much as possible. See Cholesky Decomposition for more information on the matrixS
. The Pivoted Cholesky decomposition satisfies . On input, the values from the diagonal and lowertriangular part of the matrixA
are used to construct the factorization. On output the diagonal of the input matrixA
stores the diagonal elements of , and the lower triangular portion ofA
contains the matrix . Since has ones on its diagonal these do not need to be explicitely stored. The upper triangular portion ofA
is unmodified. The permutation matrix is stored inp
on output. The diagonal scaling transformation is stored inS
on output.

int
gsl_linalg_pcholesky_solve2
(const gsl_matrix * LDLT, const gsl_permutation * p, const gsl_vector * S, const gsl_vector * b, gsl_vector * x)¶ This function solves the system using the Pivoted Cholesky decomposition of held in the matrix
LDLT
, permutationp
, and vectorS
, which must have been previously computed bygsl_linalg_pcholesky_decomp2()
.

int
gsl_linalg_pcholesky_svx2
(const gsl_matrix * LDLT, const gsl_permutation * p, const gsl_vector * S, gsl_vector * x)¶ This function solves the system inplace using the Pivoted Cholesky decomposition of held in the matrix
LDLT
, permutationp
and vectorS
, which must have been previously computed bygsl_linalg_pcholesky_decomp2()
. On input,x
contains the right hand side vector which is replaced by the solution vector on output.

int
gsl_linalg_pcholesky_invert
(const gsl_matrix * LDLT, const gsl_permutation * p, gsl_matrix * Ainv)¶ This function computes the inverse of the matrix , using the Pivoted Cholesky decomposition stored in
LDLT
andp
. On output, the matrixAinv
contains .

int
gsl_linalg_pcholesky_rcond
(const gsl_matrix * LDLT, const gsl_permutation * p, double * rcond, gsl_vector * work)¶ This function estimates the reciprocal condition number (using the 1norm) of the symmetric positive definite matrix , using its pivoted Cholesky decomposition provided in
LDLT
. The reciprocal condition number estimate, defined as , is stored inrcond
. Additional workspace of size is required inwork
.
Modified Cholesky Decomposition¶
The modified Cholesky decomposition is suitable for solving systems where is a symmetric indefinite matrix. Such matrices arise in nonlinear optimization algorithms. The standard Cholesky decomposition requires a positive definite matrix and would fail in this case. Instead of resorting to a method like QR or SVD, which do not take into account the symmetry of the matrix, we can instead introduce a small perturbation to the matrix to make it positive definite, and then use a Cholesky decomposition on the perturbed matrix. The resulting decomposition satisfies
where is a permutation matrix, is a diagonal perturbation matrix, is unit lower triangular, and is diagonal. If is sufficiently positive definite, then the perturbation matrix will be zero and this method is equivalent to the pivoted Cholesky algorithm. For indefinite matrices, the perturbation matrix is computed to ensure that is positive definite and well conditioned.

int
gsl_linalg_mcholesky_decomp
(gsl_matrix * A, gsl_permutation * p, gsl_vector * E)¶ This function factors the symmetric, indefinite square matrix
A
into the Modified Cholesky decomposition . On input, the values from the diagonal and lowertriangular part of the matrixA
are used to construct the factorization. On output the diagonal of the input matrixA
stores the diagonal elements of , and the lower triangular portion ofA
contains the matrix . Since has ones on its diagonal these do not need to be explicitely stored. The upper triangular portion ofA
is unmodified. The permutation matrix is stored inp
on output. The diagonal perturbation matrix is stored inE
on output. The parameterE
may be set to NULL if it is not required.

int
gsl_linalg_mcholesky_solve
(const gsl_matrix * LDLT, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x)¶ This function solves the perturbed system using the Cholesky decomposition of held in the matrix
LDLT
and permutationp
which must have been previously computed bygsl_linalg_mcholesky_decomp()
.

int
gsl_linalg_mcholesky_svx
(const gsl_matrix * LDLT, const gsl_permutation * p, gsl_vector * x)¶ This function solves the perturbed system inplace using the Cholesky decomposition of held in the matrix
LDLT
and permutationp
which must have been previously computed bygsl_linalg_mcholesky_decomp()
. On input,x
contains the right hand side vector which is replaced by the solution vector on output.

int
gsl_linalg_mcholesky_rcond
(const gsl_matrix * LDLT, const gsl_permutation * p, double * rcond, gsl_vector * work)¶ This function estimates the reciprocal condition number (using the 1norm) of the perturbed matrix , using its pivoted Cholesky decomposition provided in
LDLT
. The reciprocal condition number estimate, defined as , is stored inrcond
. Additional workspace of size is required inwork
.
LDLT Decomposition¶
If is a symmetric, nonsingular square matrix, then it has a unique factorization of the form
where is a unit lower triangular matrix and is diagonal. If is positive definite, then this factorization is equivalent to the Cholesky factorization, where the lower triangular Cholesky factor is . Some indefinite matrices for which no Cholesky decomposition exists have an decomposition with negative entries in . The algorithm is sometimes referred to as the square root free Cholesky decomposition, as the algorithm does not require the computation of square roots. The algorithm is stable for positive definite matrices, but is not guaranteed to be stable for indefinite matrices.

int
gsl_linalg_ldlt_decomp
(gsl_matrix * A)¶ This function factorizes the symmetric, nonsingular square matrix
A
into the decomposition . On input, the values from the diagonal and lowertriangular part of the matrixA
are used. The upper triangle ofA
is used as temporary workspace. On output the diagonal ofA
contains the matrix and the lower triangle ofA
contains the unit lower triangular matrix . The matrix 1norm, is stored in the upper right corner on output, for later use bygsl_linalg_ldlt_rcond()
.If the matrix is detected to be singular, the function returns the error code
GSL_EDOM
.

int
gsl_linalg_ldlt_solve
(const gsl_matrix * LDLT, const gsl_vector * b, gsl_vector * x)¶ This function solves the system using the decomposition of held in the matrix
LDLT
which must have been previously computed bygsl_linalg_ldlt_decomp()
.

int
gsl_linalg_ldlt_svx
(const gsl_matrix * LDLT, gsl_vector * x)¶ This function solves the system inplace using the decomposition of held in the matrix
LDLT
which must have been previously computed bygsl_linalg_ldlt_decomp()
. On inputx
should contain the righthand side , which is replaced by the solution on output.

int
gsl_linalg_ldlt_rcond
(const gsl_matrix * LDLT, double * rcond, gsl_vector * work)¶ This function estimates the reciprocal condition number (using the 1norm) of the symmetric nonsingular matrix , using its decomposition provided in
LDLT
. The reciprocal condition number estimate, defined as , is stored inrcond
. Additional workspace of size is required inwork
.
Tridiagonal Decomposition of Real Symmetric Matrices¶
A symmetric matrix can be factorized by similarity transformations into the form,
where is an orthogonal matrix and is a symmetric tridiagonal matrix.

int
gsl_linalg_symmtd_decomp
(gsl_matrix * A, gsl_vector * tau)¶ This function factorizes the symmetric square matrix
A
into the symmetric tridiagonal decomposition . On output the diagonal and subdiagonal part of the input matrixA
contain the tridiagonal matrix . The remaining lower triangular part of the input matrix contains the Householder vectors which, together with the Householder coefficientstau
, encode the orthogonal matrix . This storage scheme is the same as used by LAPACK. The upper triangular part ofA
is not referenced.

int
gsl_linalg_symmtd_unpack
(const gsl_matrix * A, const gsl_vector * tau, gsl_matrix * Q, gsl_vector * diag, gsl_vector * subdiag)¶ This function unpacks the encoded symmetric tridiagonal decomposition (
A
,tau
) obtained fromgsl_linalg_symmtd_decomp()
into the orthogonal matrixQ
, the vector of diagonal elementsdiag
and the vector of subdiagonal elementssubdiag
.

int
gsl_linalg_symmtd_unpack_T
(const gsl_matrix * A, gsl_vector * diag, gsl_vector * subdiag)¶ This function unpacks the diagonal and subdiagonal of the encoded symmetric tridiagonal decomposition (
A
,tau
) obtained fromgsl_linalg_symmtd_decomp()
into the vectorsdiag
andsubdiag
.
Tridiagonal Decomposition of Hermitian Matrices¶
A hermitian matrix can be factorized by similarity transformations into the form,
where is a unitary matrix and is a real symmetric tridiagonal matrix.

int
gsl_linalg_hermtd_decomp
(gsl_matrix_complex * A, gsl_vector_complex * tau)¶ This function factorizes the hermitian matrix
A
into the symmetric tridiagonal decomposition . On output the real parts of the diagonal and subdiagonal part of the input matrixA
contain the tridiagonal matrix . The remaining lower triangular part of the input matrix contains the Householder vectors which, together with the Householder coefficientstau
, encode the unitary matrix . This storage scheme is the same as used by LAPACK. The upper triangular part ofA
and imaginary parts of the diagonal are not referenced.

int
gsl_linalg_hermtd_unpack
(const gsl_matrix_complex * A, const gsl_vector_complex * tau, gsl_matrix_complex * U, gsl_vector * diag, gsl_vector * subdiag)¶ This function unpacks the encoded tridiagonal decomposition (
A
,tau
) obtained fromgsl_linalg_hermtd_decomp()
into the unitary matrixU
, the real vector of diagonal elementsdiag
and the real vector of subdiagonal elementssubdiag
.

int
gsl_linalg_hermtd_unpack_T
(const gsl_matrix_complex * A, gsl_vector * diag, gsl_vector * subdiag)¶ This function unpacks the diagonal and subdiagonal of the encoded tridiagonal decomposition (
A
,tau
) obtained from thegsl_linalg_hermtd_decomp()
into the real vectorsdiag
andsubdiag
.
Hessenberg Decomposition of Real Matrices¶
A general real matrix can be decomposed by orthogonal similarity transformations into the form
where is orthogonal and is an upper Hessenberg matrix, meaning that it has zeros below the first subdiagonal. The Hessenberg reduction is the first step in the Schur decomposition for the nonsymmetric eigenvalue problem, but has applications in other areas as well.

int
gsl_linalg_hessenberg_decomp
(gsl_matrix * A, gsl_vector * tau)¶ This function computes the Hessenberg decomposition of the matrix
A
by applying the similarity transformation . On output, is stored in the upper portion ofA
. The information required to construct the matrix is stored in the lower triangular portion ofA
. is a product of Householder matrices. The Householder vectors are stored in the lower portion ofA
(below the subdiagonal) and the Householder coefficients are stored in the vectortau
.tau
must be of lengthN
.

int
gsl_linalg_hessenberg_unpack
(gsl_matrix * H, gsl_vector * tau, gsl_matrix * U)¶ This function constructs the orthogonal matrix from the information stored in the Hessenberg matrix
H
along with the vectortau
.H
andtau
are outputs fromgsl_linalg_hessenberg_decomp()
.

int
gsl_linalg_hessenberg_unpack_accum
(gsl_matrix * H, gsl_vector * tau, gsl_matrix * V)¶ This function is similar to
gsl_linalg_hessenberg_unpack()
, except it accumulates the matrixU
intoV
, so that . The matrixV
must be initialized prior to calling this function. SettingV
to the identity matrix provides the same result asgsl_linalg_hessenberg_unpack()
. IfH
is orderN
, thenV
must haveN
columns but may have any number of rows.

int
gsl_linalg_hessenberg_set_zero
(gsl_matrix * H)¶ This function sets the lower triangular portion of
H
, below the subdiagonal, to zero. It is useful for clearing out the Householder vectors after callinggsl_linalg_hessenberg_decomp()
.
HessenbergTriangular Decomposition of Real Matrices¶
A general real matrix pair (, ) can be decomposed by orthogonal similarity transformations into the form
where and are orthogonal, is an upper Hessenberg matrix, and is upper triangular. The HessenbergTriangular reduction is the first step in the generalized Schur decomposition for the generalized eigenvalue problem.

int
gsl_linalg_hesstri_decomp
(gsl_matrix * A, gsl_matrix * B, gsl_matrix * U, gsl_matrix * V, gsl_vector * work)¶ This function computes the HessenbergTriangular decomposition of the matrix pair (
A
,B
). On output, is stored inA
, and is stored inB
. IfU
andV
are provided (they may be null), the similarity transformations are stored in them. Additional workspace of length is needed inwork
.
Bidiagonalization¶
A general matrix can be factorized by similarity transformations into the form,
where and are orthogonal matrices and is a
by bidiagonal matrix with nonzero entries only on the
diagonal and superdiagonal. The size of U
is by
and the size of V
is by.

int
gsl_linalg_bidiag_decomp
(gsl_matrix * A, gsl_vector * tau_U, gsl_vector * tau_V)¶ This function factorizes the by matrix
A
into bidiagonal form . The diagonal and superdiagonal of the matrix are stored in the diagonal and superdiagonal ofA
. The orthogonal matrices andV
are stored as compressed Householder vectors in the remaining elements ofA
. The Householder coefficients are stored in the vectorstau_U
andtau_V
. The length oftau_U
must equal the number of elements in the diagonal ofA
and the length oftau_V
should be one element shorter.

int
gsl_linalg_bidiag_unpack
(const gsl_matrix * A, const gsl_vector * tau_U, gsl_matrix * U, const gsl_vector * tau_V, gsl_matrix * V, gsl_vector * diag, gsl_vector * superdiag)¶ This function unpacks the bidiagonal decomposition of
A
produced bygsl_linalg_bidiag_decomp()
, (A
,tau_U
,tau_V
) into the separate orthogonal matricesU
,V
and the diagonal vectordiag
and superdiagonalsuperdiag
. Note thatU
is stored as a compact by orthogonal matrix satisfying for efficiency.

int
gsl_linalg_bidiag_unpack2
(gsl_matrix * A, gsl_vector * tau_U, gsl_vector * tau_V, gsl_matrix * V)¶ This function unpacks the bidiagonal decomposition of
A
produced bygsl_linalg_bidiag_decomp()
, (A
,tau_U
,tau_V
) into the separate orthogonal matricesU
,V
and the diagonal vectordiag
and superdiagonalsuperdiag
. The matrixU
is stored inplace inA
.

int
gsl_linalg_bidiag_unpack_B
(const gsl_matrix * A, gsl_vector * diag, gsl_vector * superdiag)¶ This function unpacks the diagonal and superdiagonal of the bidiagonal decomposition of
A
fromgsl_linalg_bidiag_decomp()
, into the diagonal vectordiag
and superdiagonal vectorsuperdiag
.
Givens Rotations¶
A Givens rotation is a rotation in the plane acting on two elements of a given vector. It can be represented in matrix form as
where the and appear at the intersection of the th and th rows and columns. When acting on a vector , performs a rotation of the elements of . Givens rotations are typically used to introduce zeros in vectors, such as during the QR decomposition of a matrix. In this case, it is typically desired to find and such that
with .

void
gsl_linalg_givens
(const double a, const double b, double * c, double * s)¶ This function computes and so that the Givens matrix acting on the vector produces , with .

void
gsl_linalg_givens_gv
(gsl_vector * v, const size_t i, const size_t j, const double c, const double s)¶ This function applies the Givens rotation defined by and to the
i
andj
elements ofv
. On output, .
Householder Transformations¶
A Householder transformation is a rank1 modification of the identity matrix which can be used to zero out selected elements of a vector. A Householder matrix takes the form,
where is a vector (called the Householder vector) and . The functions described in this section use the rank1 structure of the Householder matrix to create and apply Householder transformations efficiently.

double
gsl_linalg_householder_transform
(gsl_vector * w)¶ 
gsl_complex
gsl_linalg_complex_householder_transform
(gsl_vector_complex * w)¶ This function prepares a Householder transformation which can be used to zero all the elements of the input vector
w
except the first. On output the Householder vectorv
is stored inw
and the scalar is returned. The householder vectorv
is normalized so thatv[0] = 1
, however this 1 is not stored in the output vector. Instead,w[0]
is set to the first element of the transformed vector, so that if ,w[0] = u[0]
on output and the remainder of is zero.

int
gsl_linalg_householder_hm
(double tau, const gsl_vector * v, gsl_matrix * A)¶ 
int
gsl_linalg_complex_householder_hm
(gsl_complex tau, const gsl_vector_complex * v, gsl_matrix_complex * A)¶ This function applies the Householder matrix defined by the scalar
tau
and the vectorv
to the lefthand side of the matrixA
. On output the result is stored inA
.

int
gsl_linalg_householder_mh
(double tau, const gsl_vector * v, gsl_matrix * A)¶ 
int
gsl_linalg_complex_householder_mh
(gsl_complex tau, const gsl_vector_complex * v, gsl_matrix_complex * A)¶ This function applies the Householder matrix defined by the scalar
tau
and the vectorv
to the righthand side of the matrixA
. On output the result is stored inA
.

int
gsl_linalg_householder_hv
(double tau, const gsl_vector * v, gsl_vector * w)¶ 
int
gsl_linalg_complex_householder_hv
(gsl_complex tau, const gsl_vector_complex * v, gsl_vector_complex * w)¶ This function applies the Householder transformation defined by the scalar
tau
and the vectorv
to the vectorw
. On output the result is stored inw
.
Householder solver for linear systems¶

int
gsl_linalg_HH_solve
(gsl_matrix * A, const gsl_vector * b, gsl_vector * x)¶ This function solves the system directly using Householder transformations. On output the solution is stored in
x
andb
is not modified. The matrixA
is destroyed by the Householder transformations.

int
gsl_linalg_HH_svx
(gsl_matrix * A, gsl_vector * x)¶ This function solves the system inplace using Householder transformations. On input
x
should contain the righthand side , which is replaced by the solution on output. The matrixA
is destroyed by the Householder transformations.
Tridiagonal Systems¶
The functions described in this section efficiently solve symmetric,
nonsymmetric and cyclic tridiagonal systems with minimal storage.
Note that the current implementations of these functions use a variant
of Cholesky decomposition, so the tridiagonal matrix must be positive
definite. For nonpositive definite matrices, the functions return
the error code GSL_ESING
.

int
gsl_linalg_solve_tridiag
(const gsl_vector * diag, const gsl_vector * e, const gsl_vector * f, const gsl_vector * b, gsl_vector * x)¶ This function solves the general by system where
A
is tridiagonal (). The superdiagonal and subdiagonal vectorse
andf
must be one element shorter than the diagonal vectordiag
. The form ofA
for the 4by4 case is shown below,

int
gsl_linalg_solve_symm_tridiag
(const gsl_vector * diag, const gsl_vector * e, const gsl_vector * b, gsl_vector * x)¶ This function solves the general by system where
A
is symmetric tridiagonal (). The offdiagonal vectore
must be one element shorter than the diagonal vectordiag
. The form ofA
for the 4by4 case is shown below,

int
gsl_linalg_solve_cyc_tridiag
(const gsl_vector * diag, const gsl_vector * e, const gsl_vector * f, const gsl_vector * b, gsl_vector * x)¶ This function solves the general by system where
A
is cyclic tridiagonal (). The cyclic superdiagonal and subdiagonal vectorse
andf
must have the same number of elements as the diagonal vectordiag
. The form ofA
for the 4by4 case is shown below,

int
gsl_linalg_solve_symm_cyc_tridiag
(const gsl_vector * diag, const gsl_vector * e, const gsl_vector * b, gsl_vector * x)¶ This function solves the general by system where
A
is symmetric cyclic tridiagonal (). The cyclic offdiagonal vectore
must have the same number of elements as the diagonal vectordiag
. The form ofA
for the 4by4 case is shown below,
Triangular Systems¶

int
gsl_linalg_tri_invert
(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix * T)¶ 
int
gsl_linalg_complex_tri_invert
(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix_complex * T)¶ These functions compute the inplace inverse of the triangular matrix
T
, stored in the lower triangle whenUplo
=CblasLower
and upper triangle whenUplo
=CblasUpper
. The parameterDiag
=CblasUnit
,CblasNonUnit
specifies whether the matrix is unit triangular.

int
gsl_linalg_tri_LTL
(gsl_matrix * L)¶ 
int
gsl_linalg_complex_tri_LHL
(gsl_matrix_complex * L)¶ These functions compute the product (or ) inplace and stores it in the lower triangle of
L
on output.

int
gsl_linalg_tri_UL
(gsl_matrix * LU)¶ 
int
gsl_linalg_complex_tri_UL
(gsl_matrix_complex * LU)¶ These functions compute the product where is upper triangular and is unit lower triangular, stored in
LU
, as computed bygsl_linalg_LU_decomp()
orgsl_linalg_complex_LU_decomp()
. The product is computed inplace using Level 3 BLAS.

int
gsl_linalg_tri_rcond
(CBLAS_UPLO_t Uplo, const gsl_matrix * A, double * rcond, gsl_vector * work)¶ This function estimates the 1norm reciprocal condition number of the triangular matrix
A
, using the lower triangle whenUplo
isCblasLower
and upper triangle whenUplo
isCblasUpper
. The reciprocal condition number is stored inrcond
on output. Additional workspace of size is required inwork
.
Banded Systems¶
Band matrices are sparse matrices whose nonzero entries are confined to a diagonal band. From a storage point of view, significant savings can be achieved by storing only the nonzero diagonals of a banded matrix. Algorithms such as LU and Cholesky factorizations preserve the band structure of these matrices. Computationally, working with compact banded matrices is preferable to working on the full dense matrix with many zero entries.
General Banded Format¶
An example of a general banded matrix is given below.
This matrix has a lower bandwidth of 1 and an upper bandwidth of 2. The lower bandwidth is the number of nonzero subdiagonals, and the upper bandwidth is the number of nonzero superdiagonals. A banded matrix has a lower bandwidth and upper bandwidth . For example, diagonal matrices are , tridiagonal matrices are , and upper triangular matrices are banded matrices.
The corresponding by packed banded matrix looks like
where the superdiagonals are stored in columns, followed by the diagonal, followed by the subdiagonals. The entries marked by are not referenced by the banded routines. With this format, each row of corresponds to the nonzero entries of the corresponding column of . For an by matrix , the dimension of will be by.
Symmetric Banded Format¶
Symmetric banded matrices allow for additional storage savings. As an example, consider the following symmetric banded matrix with lower bandwidth :
The packed symmetric banded matrix will look like:
The entries marked by are not referenced by the symmetric banded routines.
Warning
Note that this format is the transpose of the symmetric banded format used by LAPACK. In order to develop efficient routines for symmetric banded matrices, it helps to have the nonzero elements in each column in contiguous memory locations. Since C uses rowmajor order, GSL stores the columns in the rows of the packed banded format, while LAPACK, written in Fortran, uses the transposed format.
Banded Cholesky Decomposition¶
The routines in this section are designed to factor and solve by linear systems of the form where is a banded, symmetric, and positive definite matrix with lower bandwidth . See Cholesky Decomposition for more information on the factorization. The lower triangular factor of the Cholesky decomposition preserves the same banded structure as the matrix , enabling an efficient algorithm which overwrites the original matrix with the factor.

int
gsl_linalg_cholesky_band_decomp
(gsl_matrix * A)¶ This function factorizes the symmetric, positivedefinite square matrix
A
into the Cholesky decomposition . The input matrixA
is given in symmetric banded format, and has dimensions by, where is the lower bandwidth of the matrix. On output, the entries ofA
are replaced by the entries of the matrix in the same format. In addition, the lower right element ofA
is used to store the matrix 1norm, used later bygsl_linalg_cholesky_band_rcond()
to calculate the reciprocal condition number.If the matrix is not positivedefinite then the decomposition will fail, returning the error code
GSL_EDOM
. When testing whether a matrix is positivedefinite, disable the error handler first to avoid triggering an error.

int
gsl_linalg_cholesky_band_solve
(const gsl_matrix * LLT, const gsl_vector * b, gsl_vector * x)¶ This function solves the symmetric banded system using the Cholesky decomposition of held in the matrix
LLT
which must have been previously computed bygsl_linalg_cholesky_band_decomp()
.

int
gsl_linalg_cholesky_band_svx
(const gsl_matrix * LLT, gsl_vector * x)¶ This function solves the symmetric banded system inplace using the Cholesky decomposition of held in the matrix
LLT
which must have been previously computed bygsl_linalg_cholesky_band_decomp()
. On inputx
should contain the righthand side , which is replaced by the solution on output.

int
gsl_linalg_cholesky_band_invert
(const gsl_matrix * LLT, gsl_matrix * Ainv)¶ This function computes the inverse of a symmetric banded matrix from its Cholesky decomposition
LLT
, which must have been previously computed bygsl_linalg_cholesky_band_decomp()
. On output, the inverse is stored inAinv
, using both the lower and upper portions.

int
gsl_linalg_cholesky_band_unpack
(const gsl_matrix * LLT, gsl_matrix * L)¶ This function unpacks the lower triangular Cholesky factor from
LLT
and stores it in the lower triangular portion of the by matrixL
. The upper triangular portion ofL
is not referenced.

int
gsl_linalg_cholesky_band_rcond
(const gsl_matrix * LLT, double * rcond, gsl_vector * work)¶ This function estimates the reciprocal condition number (using the 1norm) of the symmetric banded positive definite matrix , using its Cholesky decomposition provided in
LLT
. The reciprocal condition number estimate, defined as , is stored inrcond
. Additional workspace of size is required inwork
.
Banded LDLT Decomposition¶
The routines in this section are designed to factor and solve by linear systems of the form where is a banded, symmetric, and nonsingular matrix with lower bandwidth . See LDLT Decomposition for more information on the factorization. The lower triangular factor of the decomposition preserves the same banded structure as the matrix , enabling an efficient algorithm which overwrites the original matrix with the and factors.

int
gsl_linalg_ldlt_band_decomp
(gsl_matrix * A)¶ This function factorizes the symmetric, nonsingular square matrix
A
into the decomposition . The input matrixA
is given in symmetric banded format, and has dimensions by, where is the lower bandwidth of the matrix. On output, the entries ofA
are replaced by the entries of the matrices and in the same format.If the matrix is singular then the decomposition will fail, returning the error code
GSL_EDOM
.

int
gsl_linalg_ldlt_band_solve
(const gsl_matrix * LDLT, const gsl_vector * b, gsl_vector * x)¶ This function solves the symmetric banded system using the decomposition of held in the matrix
LDLT
which must have been previously computed bygsl_linalg_ldlt_band_decomp()
.

int
gsl_linalg_ldlt_band_svx
(const gsl_matrix * LDLT, gsl_vector * x)¶ This function solves the symmetric banded system inplace using the decomposition of held in the matrix
LDLT
which must have been previously computed bygsl_linalg_ldlt_band_decomp()
. On inputx
should contain the righthand side , which is replaced by the solution on output.

gsl_linalg_ldlt_band_unpack
(const gsl_matrix * LDLT, gsl_matrix * L, gsl_vector * D)¶ This function unpacks the unit lower triangular factor from
LDLT
and stores it in the lower triangular portion of the by matrixL
. The upper triangular portion ofL
is not referenced. The diagonal matrix is stored in the vectorD
.

int
gsl_linalg_ldlt_band_rcond
(const gsl_matrix * LDLT, double * rcond, gsl_vector * work)¶ This function estimates the reciprocal condition number (using the 1norm) of the symmetric banded nonsingular matrix , using its decomposition provided in
LDLT
. The reciprocal condition number estimate, defined as , is stored inrcond
. Additional workspace of size is required inwork
.
Balancing¶
The process of balancing a matrix applies similarity transformations to make the rows and columns have comparable norms. This is useful, for example, to reduce roundoff errors in the solution of eigenvalue problems. Balancing a matrix consists of replacing with a similar matrix
where is a diagonal matrix whose entries are powers of the floating point radix.

int
gsl_linalg_balance_matrix
(gsl_matrix * A, gsl_vector * D)¶ This function replaces the matrix
A
with its balanced counterpart and stores the diagonal elements of the similarity transformation into the vectorD
.
Examples¶
The following program solves the linear system . The system to be solved is,
and the solution is found using LU decomposition of the matrix .
#include <stdio.h>
#include <gsl/gsl_linalg.h>
int
main (void)
{
double a_data[] = { 0.18, 0.60, 0.57, 0.96,
0.41, 0.24, 0.99, 0.58,
0.14, 0.30, 0.97, 0.66,
0.51, 0.13, 0.19, 0.85 };
double b_data[] = { 1.0, 2.0, 3.0, 4.0 };
gsl_matrix_view m
= gsl_matrix_view_array (a_data, 4, 4);
gsl_vector_view b
= gsl_vector_view_array (b_data, 4);
gsl_vector *x = gsl_vector_alloc (4);
int s;
gsl_permutation * p = gsl_permutation_alloc (4);
gsl_linalg_LU_decomp (&m.matrix, p, &s);
gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x);
printf ("x = \n");
gsl_vector_fprintf (stdout, x, "%g");
gsl_permutation_free (p);
gsl_vector_free (x);
return 0;
}
Here is the output from the program,
x =
4.05205
12.6056
1.66091
8.69377
This can be verified by multiplying the solution by the original matrix using GNU octave,
octave> A = [ 0.18, 0.60, 0.57, 0.96;
0.41, 0.24, 0.99, 0.58;
0.14, 0.30, 0.97, 0.66;
0.51, 0.13, 0.19, 0.85 ];
octave> x = [ 4.05205; 12.6056; 1.66091; 8.69377];
octave> A * x
ans =
1.0000
2.0000
3.0000
4.0000
This reproduces the original righthand side vector, , in accordance with the equation .
References and Further Reading¶
Further information on the algorithms described in this section can be found in the following book,
G. H. Golub, C. F. Van Loan, “Matrix Computations” (3rd Ed, 1996), Johns Hopkins University Press, ISBN 0801854148.
The LAPACK library is described in the following manual,
LAPACK Users’ Guide (Third Edition, 1999), Published by SIAM, ISBN 0898714478
The LAPACK source code can be found at http://www.netlib.org/lapack, along with an online copy of the users guide.
Further information on recursive Level 3 BLAS algorithms may be found in the following paper,
E. Peise and P. Bientinesi, “Recursive algorithms for dense linear algebra: the ReLAPACK collection”, http://arxiv.org/abs/1602.06763, 2016.
The recursive Level 3 BLAS QR decomposition is described in the following paper,
E. Elmroth and F. G. Gustavson, 2000. Applying recursion to serial and parallel QR factorization leads to better performance. IBM Journal of Research and Development, 44(4), pp.605624.
The Modified GolubReinsch algorithm is described in the following paper,
T.F. Chan, “An Improved Algorithm for Computing the Singular Value Decomposition”, ACM Transactions on Mathematical Software, 8 (1982), pp 72–83.
The Jacobi algorithm for singular value decomposition is described in the following papers,
J.C. Nash, “A onesided transformation method for the singular value decomposition and algebraic eigenproblem”, Computer Journal, Volume 18, Number 1 (1975), p 74–76
J.C. Nash and S. Shlien “Simple algorithms for the partial singular value decomposition”, Computer Journal, Volume 30 (1987), p 268–275.
J. Demmel, K. Veselic, “Jacobi’s Method is more accurate than QR”, Lapack Working Note 15 (LAWN15), October 1989. Available from netlib, http://www.netlib.org/lapack/ in the
lawns
orlawnspdf
directories.
The algorithm for estimating a matrix condition number is described in the following paper,
N. J. Higham, “FORTRAN codes for estimating the onenorm of a real or complex matrix, with applications to condition estimation”, ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381396, December 1988.