linalg.texi

开放gsl矩阵运算

@cindex linear algebra
@cindex solution of linear systems, Ax=b
@cindex matrix factorization
@cindex factorization of matrices

This chapter describes functions for solving linear systems.  The
library provides simple linear algebra operations which operate directly
on the @code{gsl_vector} and @code{gsl_matrix} objects.  These are
intended for use with "small" systems where simple algorithms are
acceptable.

@cindex LAPACK, recommended for linear algebra
Anyone interested in large systems will want to use the sophisticated
routines found in @sc{lapack}. The Fortran version of @sc{lapack} is
recommended as the standard package for linear algebra.  It supports
blocked algorithms, specialized data representations and other
optimizations.

The functions described in this chapter are declared in the header file
@file{gsl_linalg.h}.


@menu
* LU Decomposition::            
* QR Decomposition::            
* QR Decomposition with Column Pivoting::  
* Singular Value Decomposition::  
* Cholesky Decomposition::      
* Tridiagonal Decomposition of Real Symmetric Matrices::  
* Tridiagonal Decomposition of Hermitian Matrices::  
* Bidiagonalization::           
* Householder solver for linear systems::  
* Tridiagonal Systems::         
* Linear Algebra Examples::     
* Linear Algebra References and Further Reading::  
@end menu

@node LU Decomposition
@section LU Decomposition
@cindex LU decomposition

A general square matrix @math{A} has an @math{LU} decomposition into
upper and lower triangular matrices,

@tex
\beforedisplay
$$
P A = L U
$$
\afterdisplay
@end tex
@ifinfo
@example
P A = L U
@end example
@end ifinfo
@noindent
where @math{P} is a permutation matrix, @math{L} is unit lower
triangular matrix and @math{U} is upper triangular matrix. For square
matrices this decomposition can be used to convert the linear system
@math{A x = b} into a pair of triangular systems (@math{L y = P b},
@math{U x = y}), which can be solved by forward and back-substitution.

@deftypefun int gsl_linalg_LU_decomp (gsl_matrix * @var{A}, gsl_permutation * @var{p}, int *@var{signum})
@deftypefunx int gsl_linalg_complex_LU_decomp (gsl_matrix_complex * @var{A}, gsl_permutation * @var{p}, int *@var{signum})
These functions factorize the square matrix @var{A} into the @math{LU}
decomposition @math{PA = LU}.  On output the diagonal and upper
triangular part of the input matrix @var{A} contain the matrix
@math{U}. The lower triangular part of the input matrix (excluding the
diagonal) contains @math{L}.  The diagonal elements of @math{L} are
unity, and are not stored.

The permutation matrix @math{P} is encoded in the permutation
@var{p}. The @math{j}-th column of the matrix @math{P} is given by the
@math{k}-th column of the identity matrix, where @math{k = p_j} the
@math{j}-th element of the permutation vector. The sign of the
permutation is given by @var{signum}. It has the value @math{(-1)^n},
where @math{n} is the number of interchanges in the permutation.

The algorithm used in the decomposition is Gaussian Elimination with
partial pivoting (Golub & Van Loan, @cite{Matrix Computations},
Algorithm 3.4.1).
@end deftypefun

@cindex linear systems, solution of
@deftypefun int gsl_linalg_LU_solve (const gsl_matrix * @var{LU}, const gsl_permutation * @var{p}, const gsl_vector * @var{b}, gsl_vector * @var{x})
@deftypefunx int gsl_linalg_complex_LU_solve (const gsl_matrix_complex * @var{LU}, const gsl_permutation * @var{p}, const gsl_vector_complex * @var{b}, gsl_vector_complex * @var{x})
These functions solve the system @math{A x = b} using the @math{LU}
decomposition of @math{A} into (@var{LU}, @var{p}) given by
@code{gsl_linalg_LU_decomp} or @code{gsl_linalg_complex_LU_decomp}.
@end deftypefun

@deftypefun int gsl_linalg_LU_svx (const gsl_matrix * @var{LU}, const gsl_permutation * @var{p}, gsl_vector * @var{x})
@deftypefunx int gsl_linalg_complex_LU_svx (const gsl_matrix_complex * @var{LU}, const gsl_permutation * @var{p}, gsl_vector_complex * @var{x})
These functions solve the system @math{A x = b} in-place using the
@math{LU} decomposition of @math{A} into (@var{LU},@var{p}). On input
@var{x} should contain the right-hand side @math{b}, which is replaced
by the solution on output.
@end deftypefun

@cindex refinement of solutions in linear systems
@cindex iterative refinement of solutions in linear systems
@cindex linear systems, refinement of solutions
@deftypefun int gsl_linalg_LU_refine (const gsl_matrix * @var{A}, const gsl_matrix * @var{LU}, const gsl_permutation * @var{p}, const gsl_vector * @var{b}, gsl_vector * @var{x}, gsl_vector * @var{residual})
@deftypefunx int gsl_linalg_complex_LU_refine (const gsl_matrix_complex * @var{A}, const gsl_matrix_complex * @var{LU}, const gsl_permutation * @var{p}, const gsl_vector_complex * @var{b}, gsl_vector_complex * @var{x}, gsl_vector_complex * @var{residual})
These functions apply an iterative improvement to @var{x}, the solution
of @math{A x = b}, using the @math{LU} decomposition of @math{A} into
(@var{LU},@var{p}). The initial residual @math{r = A x - b} is also
computed and stored in @var{residual}.
@end deftypefun

@cindex inverse of a matrix, by LU decomposition
@cindex matrix inverse
@deftypefun int gsl_linalg_LU_invert (const gsl_matrix * @var{LU}, const gsl_permutation * @var{p}, gsl_matrix * @var{inverse})
@deftypefunx int gsl_complex_linalg_LU_invert (const gsl_matrix_complex * @var{LU}, const gsl_permutation * @var{p}, gsl_matrix_complex * @var{inverse})
These functions compute the inverse of a matrix @math{A} from its
@math{LU} decomposition (@var{LU},@var{p}), storing the result in the
matrix @var{inverse}. The inverse is computed by solving the system
@math{A x = b} for each column of the identity matrix.
@end deftypefun

@cindex determinant of a matrix, by LU decomposition
@cindex matrix determinant
@deftypefun double gsl_linalg_LU_det (gsl_matrix * @var{LU}, int @var{signum})
@deftypefunx gsl_complex gsl_linalg_complex_LU_det (gsl_matrix_complex * @var{LU}, int @var{signum})
These functions compute the determinant of a matrix @math{A} from its
@math{LU} decomposition, @var{LU}. The determinant is computed as the
product of the diagonal elements of @math{U} and the sign of the row
permutation @var{signum}.
@end deftypefun

@cindex logarithm of the determinant of a matrix
@deftypefun double gsl_linalg_LU_lndet (gsl_matrix * @var{LU})
@deftypefunx double gsl_linalg_complex_LU_lndet (gsl_matrix_complex * @var{LU})
These functions compute the logarithm of the absolute value of the
determinant of a matrix @math{A}, @math{\ln|det(A)|}, from its @math{LU}
decomposition, @var{LU}.  This function may be useful if the direct
computation of the determinant would overflow or underflow.
@end deftypefun

@cindex sign of the determinant of a matrix
@deftypefun int gsl_linalg_LU_sgndet (gsl_matrix * @var{LU}, int @var{signum})
@deftypefunx gsl_complex gsl_linalg_complex_LU_sgndet (gsl_matrix_complex * @var{LU}, int @var{signum})
These functions compute the sign or phase factor of the determinant of a
matrix @math{A}, @math{det(A)/|det(A)|}, from its @math{LU} decomposition,
@var{LU}.
@end deftypefun

@node QR Decomposition
@section QR Decomposition
@cindex QR decomposition

A general rectangular @math{M}-by-@math{N} matrix @math{A} has a
@math{QR} decomposition into the product of an orthogonal
@math{M}-by-@math{M} square matrix @math{Q} (where @math{Q^T Q = I}) and
an @math{M}-by-@math{N} right-triangular matrix @math{R},

@tex
\beforedisplay
$$
A = Q R
$$
\afterdisplay
@end tex
@ifinfo
@example
A = Q R
@end example
@end ifinfo
@noindent
This decomposition can be used to convert the linear system @math{A x =
b} into the triangular system @math{R x = Q^T b}, which can be solved by
back-substitution. Another use of the @math{QR} decomposition is to
compute an orthonormal basis for a set of vectors. The first @math{N}
columns of @math{Q} form an orthonormal basis for the range of @math{A},
@math{ran(A)}, when @math{A} has full column rank.

@deftypefun int gsl_linalg_QR_decomp (gsl_matrix * @var{A}, gsl_vector * @var{tau})
This function factorizes the @math{M}-by-@math{N} matrix @var{A} into
the @math{QR} decomposition @math{A = Q R}.  On output the diagonal and
upper triangular part of the input matrix contain the matrix
@math{R}. The vector @var{tau} and the columns of the lower triangular
part of the matrix @var{A} contain the Householder coefficients and
Householder vectors which encode the orthogonal matrix @var{Q}.  The
vector @var{tau} must be of length @math{k=\min(M,N)}. The matrix
@math{Q} is related to these components by, @math{Q = Q_k ... Q_2 Q_1}
where @math{Q_i = I - \tau_i v_i v_i^T} and @math{v_i} is the
Householder vector @math{v_i =
(0,...,1,A(i+1,i),A(i+2,i),...,A(m,i))}. This is the same storage scheme
as used by @sc{lapack}.

The algorithm used to perform the decomposition is Householder QR (Golub
& Van Loan, @cite{Matrix Computations}, Algorithm 5.2.1).
@end deftypefun

@deftypefun int gsl_linalg_QR_solve (const gsl_matrix * @var{QR}, const gsl_vector * @var{tau}, const gsl_vector * @var{b}, gsl_vector * @var{x})
This function solves the system @math{A x = b} using the @math{QR}
decomposition of @math{A} into (@var{QR}, @var{tau}) given by
@code{gsl_linalg_QR_decomp}.
@end deftypefun

@deftypefun int gsl_linalg_QR_svx (const gsl_matrix * @var{QR}, const gsl_vector * @var{tau}, gsl_vector * @var{x})
This function solves the system @math{A x = b} in-place using the
@math{QR} decomposition of @math{A} into (@var{QR},@var{tau}) given by
@code{gsl_linalg_QR_decomp}. On input @var{x} should contain the
right-hand side @math{b}, which is replaced by the solution on output.
@end deftypefun

@deftypefun int gsl_linalg_QR_lssolve (const gsl_matrix * @var{QR}, const gsl_vector * @var{tau}, const gsl_vector * @var{b}, gsl_vector * @var{x}, gsl_vector * @var{residual})
This function finds the least squares solution to the overdetermined
system @math{A x = b} where the matrix @var{A} has more rows than
columns.  The least squares solution minimizes the Euclidean norm of the
residual, @math{||Ax - b||}.The routine uses the @math{QR} decomposition
of @math{A} into (@var{QR}, @var{tau}) given by
@code{gsl_linalg_QR_decomp}.  The solution is returned in @var{x}.  The
residual is computed as a by-product and stored in @var{residual}.
@end deftypefun

@deftypefun int gsl_linalg_QR_QTvec (const gsl_matrix * @var{QR}, const gsl_vector * @var{tau}, gsl_vector * @var{v})
This function applies the matrix @math{Q^T} encoded in the decomposition
(@var{QR},@var{tau}) to the vector @var{v}, storing the result @math{Q^T
v} in @var{v}.  The matrix multiplication is carried out directly using
the encoding of the Householder vectors without needing to form the full
matrix @math{Q^T}.
@end deftypefun

@deftypefun int gsl_linalg_QR_Qvec (const gsl_matrix * @var{QR}, const gsl_vector * @var{tau}, gsl_vector * @var{v})
This function applies the matrix @math{Q} encoded in the decomposition
(@var{QR},@var{tau}) to the vector @var{v}, storing the result @math{Q
v} in @var{v}.  The matrix multiplication is carried out directly using
the encoding of the Householder vectors without needing to form the full
matrix @math{Q}.
@end deftypefun

@deftypefun int gsl_linalg_QR_Rsolve (const gsl_matrix * @var{QR}, const gsl_vector * @var{b}, gsl_vector * @var{x})
This function solves the triangular system @math{R x = b} for
@var{x}. It may be useful if the product @math{b' = Q^T b} has already
been computed using @code{gsl_linalg_QR_QTvec}.
@end deftypefun

@deftypefun int gsl_linalg_QR_Rsvx (const gsl_matrix * @var{QR}, gsl_vector * @var{x})
This function solves the triangular system @math{R x = b} for @var{x}
in-place. On input @var{x} should contain the right-hand side @math{b}
and is replaced by the solution on output. This function may be useful if
the product @math{b' = Q^T b} has already been computed using
@code{gsl_linalg_QR_QTvec}.
@end deftypefun

@deftypefun int gsl_linalg_QR_unpack (const gsl_matrix * @var{QR}, const gsl_vector * @var{tau}, gsl_matrix * @var{Q}, gsl_matrix * @var{R})
This function unpacks the encoded @math{QR} decomposition
(@var{QR},@var{tau}) into the matrices @var{Q} and @var{R}, where
@var{Q} is @math{M}-by-@math{M} and @var{R} is @math{M}-by-@math{N}. 
@end deftypefun

@deftypefun int gsl_linalg_QR_QRsolve (gsl_matrix * @var{Q}, gsl_matrix * @var{R}, const gsl_vector * @var{b}, gsl_vector * @var{x})
This function solves the system @math{R x = Q^T b} for @var{x}. It can
be used when the @math{QR} decomposition of a matrix is available in
unpacked form as (@var{Q},@var{R}).
@end deftypefun

@deftypefun int gsl_linalg_QR_update (gsl_matrix * @var{Q}, gsl_matrix * @var{R}, gsl_vector * @var{w}, const gsl_vector * @var{v})
This function performs a rank-1 update @math{w v^T} of the @math{QR}
decomposition (@var{Q}, @var{R}). The update is given by @math{Q'R' = Q
R + w v^T} where the output matrices @math{Q'} and @math{R'} are also
orthogonal and right triangular. Note that @var{w} is destroyed by the
update.
@end deftypefun

@deftypefun int gsl_linalg_R_solve (const gsl_matrix * @var{R}, const gsl_vector * @var{b}, gsl_vector * @var{x})
This function solves the triangular system @math{R x = b} for the
@math{N}-by-@math{N} matrix @var{R}.
@end deftypefun

@deftypefun int gsl_linalg_R_svx (const gsl_matrix * @var{R}, gsl_vector * @var{x})
This function solves the triangular system @math{R x = b} in-place. On
input @var{x} should contain the right-hand side @math{b}, which is
replaced by the solution on output.
@end deftypefun

@node QR Decomposition with Column Pivoting
@section QR Decomposition with Column Pivoting
@cindex QR decomposition with column pivoting

The @math{QR} decomposition can be extended to the rank deficient case
by introducing a column permutation @math{P},

@tex
\beforedisplay
$$
A P = Q R
$$
\afterdisplay
@end tex
@ifinfo
@example
A P = Q R
@end example
@end ifinfo
@noindent
The first @math{r} columns of this @math{Q} form an orthonormal basis
for the range of @math{A} for a matrix with column rank @math{r}.  This
decomposition can also be used to convert the linear system @math{A x =