|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: linear algebra
|
|
Packit |
67cb25 |
single: solution of linear systems, Ax=b
|
|
Packit |
67cb25 |
single: matrix factorization
|
|
Packit |
67cb25 |
single: factorization of matrices
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
**************
|
|
Packit |
67cb25 |
Linear Algebra
|
|
Packit |
67cb25 |
**************
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: include.rst
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This chapter describes functions for solving linear systems. The
|
|
Packit |
67cb25 |
library provides linear algebra operations which operate directly on
|
|
Packit |
67cb25 |
the :type:`gsl_vector` and :type:`gsl_matrix` objects. These routines
|
|
Packit |
67cb25 |
use the standard algorithms from Golub & Van Loan's *Matrix
|
|
Packit |
67cb25 |
Computations* with Level-1 and Level-2 BLAS calls for efficiency.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The functions described in this chapter are declared in the header file
|
|
Packit |
67cb25 |
:file:`gsl_linalg.h`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: LU decomposition
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
LU Decomposition
|
|
Packit |
67cb25 |
================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A general :math:`N`-by-:math:`N` square matrix :math:`A` has an :math:`LU` decomposition into
|
|
Packit |
67cb25 |
upper and lower triangular matrices,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: P A = L U
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`P` is a permutation matrix, :math:`L` is unit lower
|
|
Packit |
67cb25 |
triangular matrix and :math:`U` is upper triangular matrix. For square
|
|
Packit |
67cb25 |
matrices this decomposition can be used to convert the linear system
|
|
Packit |
67cb25 |
:math:`A x = b` into a pair of triangular systems (:math:`L y = P b`,
|
|
Packit |
67cb25 |
:math:`U x = y`), which can be solved by forward and back-substitution.
|
|
Packit |
67cb25 |
Note that the :math:`LU` decomposition is valid for singular matrices.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_LU_decomp (gsl_matrix * A, gsl_permutation * p, int * signum)
|
|
Packit |
67cb25 |
int gsl_linalg_complex_LU_decomp (gsl_matrix_complex * A, gsl_permutation * p, int * signum)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions factorize the square matrix :data:`A` into the :math:`LU`
|
|
Packit |
67cb25 |
decomposition :math:`PA = LU`. On output the diagonal and upper
|
|
Packit |
67cb25 |
triangular part of the input matrix :data:`A` contain the matrix
|
|
Packit |
67cb25 |
:math:`U`. The lower triangular part of the input matrix (excluding the
|
|
Packit |
67cb25 |
diagonal) contains :math:`L`. The diagonal elements of :math:`L` are
|
|
Packit |
67cb25 |
unity, and are not stored.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The permutation matrix :math:`P` is encoded in the permutation
|
|
Packit |
67cb25 |
:data:`p` on output. The :math:`j`-th column of the matrix :math:`P`
|
|
Packit |
67cb25 |
is given by the :math:`k`-th column of the identity matrix, where
|
|
Packit |
67cb25 |
:math:`k = p_j` the
|
|
Packit |
67cb25 |
:math:`j`-th element of the permutation vector. The sign of the
|
|
Packit |
67cb25 |
permutation is given by :data:`signum`. It has the value :math:`(-1)^n`,
|
|
Packit |
67cb25 |
where :math:`n` is the number of interchanges in the permutation.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The algorithm used in the decomposition is Gaussian Elimination with
|
|
Packit |
67cb25 |
partial pivoting (Golub & Van Loan, *Matrix Computations*,
|
|
Packit |
67cb25 |
Algorithm 3.4.1).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: linear systems, solution of
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_LU_solve (const gsl_matrix * LU, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
int gsl_linalg_complex_LU_solve (const gsl_matrix_complex * LU, const gsl_permutation * p, const gsl_vector_complex * b, gsl_vector_complex * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions solve the square system :math:`A x = b` using the :math:`LU`
|
|
Packit |
67cb25 |
decomposition of :math:`A` into (:data:`LU`, :data:`p`) given by
|
|
Packit |
67cb25 |
:func:`gsl_linalg_LU_decomp` or :func:`gsl_linalg_complex_LU_decomp` as input.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_LU_svx (const gsl_matrix * LU, const gsl_permutation * p, gsl_vector * x)
|
|
Packit |
67cb25 |
int gsl_linalg_complex_LU_svx (const gsl_matrix_complex * LU, const gsl_permutation * p, gsl_vector_complex * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions solve the square system :math:`A x = b` in-place using the
|
|
Packit |
67cb25 |
precomputed :math:`LU` decomposition of :math:`A` into (:data:`LU`, :data:`p`). On input
|
|
Packit |
67cb25 |
:data:`x` should contain the right-hand side :math:`b`, which is replaced
|
|
Packit |
67cb25 |
by the solution on output.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: refinement of solutions in linear systems
|
|
Packit |
67cb25 |
single: iterative refinement of solutions in linear systems
|
|
Packit |
67cb25 |
single: linear systems, refinement of solutions
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: 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)
|
|
Packit |
67cb25 |
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)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions apply an iterative improvement to :data:`x`, the solution
|
|
Packit |
67cb25 |
of :math:`A x = b`, from the precomputed :math:`LU` decomposition of :math:`A` into
|
|
Packit |
67cb25 |
(:data:`LU`, :data:`p`). Additional workspace of length :data:`N` is required in :data:`work`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: inverse of a matrix, by LU decomposition
|
|
Packit |
67cb25 |
single: matrix inverse
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_LU_invert (const gsl_matrix * LU, const gsl_permutation * p, gsl_matrix * inverse)
|
|
Packit |
67cb25 |
int gsl_linalg_complex_LU_invert (const gsl_matrix_complex * LU, const gsl_permutation * p, gsl_matrix_complex * inverse)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions compute the inverse of a matrix :math:`A` from its
|
|
Packit |
67cb25 |
:math:`LU` decomposition (:data:`LU`, :data:`p`), storing the result in the
|
|
Packit |
67cb25 |
matrix :data:`inverse`. The inverse is computed by solving the system
|
|
Packit |
67cb25 |
:math:`A x = b` for each column of the identity matrix. It is preferable
|
|
Packit |
67cb25 |
to avoid direct use of the inverse whenever possible, as the linear
|
|
Packit |
67cb25 |
solver functions can obtain the same result more efficiently and
|
|
Packit |
67cb25 |
reliably (consult any introductory textbook on numerical linear algebra
|
|
Packit |
67cb25 |
for details).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: determinant of a matrix, by LU decomposition
|
|
Packit |
67cb25 |
single: matrix determinant
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: double gsl_linalg_LU_det (gsl_matrix * LU, int signum)
|
|
Packit |
67cb25 |
gsl_complex gsl_linalg_complex_LU_det (gsl_matrix_complex * LU, int signum)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions compute the determinant of a matrix :math:`A` from its
|
|
Packit |
67cb25 |
:math:`LU` decomposition, :data:`LU`. The determinant is computed as the
|
|
Packit |
67cb25 |
product of the diagonal elements of :math:`U` and the sign of the row
|
|
Packit |
67cb25 |
permutation :data:`signum`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: logarithm of the determinant of a matrix
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: double gsl_linalg_LU_lndet (gsl_matrix * LU)
|
|
Packit |
67cb25 |
double gsl_linalg_complex_LU_lndet (gsl_matrix_complex * LU)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions compute the logarithm of the absolute value of the
|
|
Packit |
67cb25 |
determinant of a matrix :math:`A`, :math:`\ln|\det(A)|`, from its :math:`LU`
|
|
Packit |
67cb25 |
decomposition, :data:`LU`. This function may be useful if the direct
|
|
Packit |
67cb25 |
computation of the determinant would overflow or underflow.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: sign of the determinant of a matrix
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_LU_sgndet (gsl_matrix * LU, int signum)
|
|
Packit |
67cb25 |
gsl_complex gsl_linalg_complex_LU_sgndet (gsl_matrix_complex * LU, int signum)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions compute the sign or phase factor of the determinant of a
|
|
Packit |
67cb25 |
matrix :math:`A`, :math:`\det(A)/|\det(A)|`, from its :math:`LU` decomposition,
|
|
Packit |
67cb25 |
:data:`LU`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: QR decomposition
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
QR Decomposition
|
|
Packit |
67cb25 |
================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A general rectangular :math:`M`-by-:math:`N` matrix :math:`A` has a
|
|
Packit |
67cb25 |
:math:`QR` decomposition into the product of an orthogonal
|
|
Packit |
67cb25 |
:math:`M`-by-:math:`M` square matrix :math:`Q` (where :math:`Q^T Q = I`) and
|
|
Packit |
67cb25 |
an :math:`M`-by-:math:`N` right-triangular matrix :math:`R`,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: A = Q R
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This decomposition can be used to convert the linear system :math:`A x = b`
|
|
Packit |
67cb25 |
into the triangular system :math:`R x = Q^T b`, which can be solved by
|
|
Packit |
67cb25 |
back-substitution. Another use of the :math:`QR` decomposition is to
|
|
Packit |
67cb25 |
compute an orthonormal basis for a set of vectors. The first :math:`N`
|
|
Packit |
67cb25 |
columns of :math:`Q` form an orthonormal basis for the range of :math:`A`,
|
|
Packit |
67cb25 |
:math:`ran(A)`, when :math:`A` has full column rank.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_QR_decomp (gsl_matrix * A, gsl_vector * tau)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function factorizes the :math:`M`-by-:math:`N` matrix :data:`A` into
|
|
Packit |
67cb25 |
the :math:`QR` decomposition :math:`A = Q R`. On output the diagonal and
|
|
Packit |
67cb25 |
upper triangular part of the input matrix contain the matrix
|
|
Packit |
67cb25 |
:math:`R`. The vector :data:`tau` and the columns of the lower triangular
|
|
Packit |
67cb25 |
part of the matrix :data:`A` contain the Householder coefficients and
|
|
Packit |
67cb25 |
Householder vectors which encode the orthogonal matrix :data:`Q`. The
|
|
Packit |
67cb25 |
vector :data:`tau` must be of length :math:`k=\min(M,N)`. The matrix
|
|
Packit |
67cb25 |
:math:`Q` is related to these components by, :math:`Q = Q_k ... Q_2 Q_1`
|
|
Packit |
67cb25 |
where :math:`Q_i = I - \tau_i v_i v_i^T` and :math:`v_i` is the
|
|
Packit |
67cb25 |
Householder vector :math:`v_i = (0,...,1,A(i+1,i),A(i+2,i),...,A(m,i))`.
|
|
Packit |
67cb25 |
This is the same storage scheme as used by |lapack|.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The algorithm used to perform the decomposition is Householder QR (Golub
|
|
Packit |
67cb25 |
& Van Loan, "Matrix Computations", Algorithm 5.2.1).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_QR_solve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the square system :math:`A x = b` using the :math:`QR`
|
|
Packit |
67cb25 |
decomposition of :math:`A` held in (:data:`QR`, :data:`tau`) which must
|
|
Packit |
67cb25 |
have been computed previously with :func:`gsl_linalg_QR_decomp`.
|
|
Packit |
67cb25 |
The least-squares solution for
|
|
Packit |
67cb25 |
rectangular systems can be found using :func:`gsl_linalg_QR_lssolve`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_QR_svx (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the square system :math:`A x = b` in-place using
|
|
Packit |
67cb25 |
the :math:`QR` decomposition of :math:`A` held in (:data:`QR`, :data:`tau`)
|
|
Packit |
67cb25 |
which must have been computed previously by
|
|
Packit |
67cb25 |
:func:`gsl_linalg_QR_decomp`. On input :data:`x` should contain the
|
|
Packit |
67cb25 |
right-hand side :math:`b`, which is replaced by the solution on output.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_QR_lssolve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_vector * b, gsl_vector * x, gsl_vector * residual)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function finds the least squares solution to the overdetermined
|
|
Packit |
67cb25 |
system :math:`A x = b` where the matrix :data:`A` has more rows than
|
|
Packit |
67cb25 |
columns. The least squares solution minimizes the Euclidean norm of the
|
|
Packit |
67cb25 |
residual, :math:`||Ax - b||`.The routine requires as input
|
|
Packit |
67cb25 |
the :math:`QR` decomposition
|
|
Packit |
67cb25 |
of :math:`A` into (:data:`QR`, :data:`tau`) given by
|
|
Packit |
67cb25 |
:func:`gsl_linalg_QR_decomp`. The solution is returned in :data:`x`. The
|
|
Packit |
67cb25 |
residual is computed as a by-product and stored in :data:`residual`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_QR_QTvec (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * v)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function applies the matrix :math:`Q^T` encoded in the decomposition
|
|
Packit |
67cb25 |
(:data:`QR`, :data:`tau`) to the vector :data:`v`, storing the result :math:`Q^T v`
|
|
Packit |
67cb25 |
in :data:`v`. The matrix multiplication is carried out directly using
|
|
Packit |
67cb25 |
the encoding of the Householder vectors without needing to form the full
|
|
Packit |
67cb25 |
matrix :math:`Q^T`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_QR_Qvec (const gsl_matrix * QR, const gsl_vector * tau, gsl_vector * v)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function applies the matrix :math:`Q` encoded in the decomposition
|
|
Packit |
67cb25 |
(:data:`QR`, :data:`tau`) to the vector :data:`v`, storing the result :math:`Q v`
|
|
Packit |
67cb25 |
in :data:`v`. The matrix multiplication is carried out directly using
|
|
Packit |
67cb25 |
the encoding of the Householder vectors without needing to form the full
|
|
Packit |
67cb25 |
matrix :math:`Q`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_QR_QTmat (const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * A)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function applies the matrix :math:`Q^T` encoded in the decomposition
|
|
Packit |
67cb25 |
(:data:`QR`, :data:`tau`) to the matrix :data:`A`, storing the result :math:`Q^T A`
|
|
Packit |
67cb25 |
in :data:`A`. The matrix multiplication is carried out directly using
|
|
Packit |
67cb25 |
the encoding of the Householder vectors without needing to form the full
|
|
Packit |
67cb25 |
matrix :math:`Q^T`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_QR_Rsolve (const gsl_matrix * QR, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the triangular system :math:`R x = b` for
|
|
Packit |
67cb25 |
:data:`x`. It may be useful if the product :math:`b' = Q^T b` has already
|
|
Packit |
67cb25 |
been computed using :func:`gsl_linalg_QR_QTvec`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_QR_Rsvx (const gsl_matrix * QR, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the triangular system :math:`R x = b` for :data:`x`
|
|
Packit |
67cb25 |
in-place. On input :data:`x` should contain the right-hand side :math:`b`
|
|
Packit |
67cb25 |
and is replaced by the solution on output. This function may be useful if
|
|
Packit |
67cb25 |
the product :math:`b' = Q^T b` has already been computed using
|
|
Packit |
67cb25 |
:func:`gsl_linalg_QR_QTvec`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_QR_unpack (const gsl_matrix * QR, const gsl_vector * tau, gsl_matrix * Q, gsl_matrix * R)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function unpacks the encoded :math:`QR` decomposition
|
|
Packit |
67cb25 |
(:data:`QR`, :data:`tau`) into the matrices :data:`Q` and :data:`R`, where
|
|
Packit |
67cb25 |
:data:`Q` is :math:`M`-by-:math:`M` and :data:`R` is :math:`M`-by-:math:`N`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_QR_QRsolve (gsl_matrix * Q, gsl_matrix * R, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the system :math:`R x = Q^T b` for :data:`x`. It can
|
|
Packit |
67cb25 |
be used when the :math:`QR` decomposition of a matrix is available in
|
|
Packit |
67cb25 |
unpacked form as (:data:`Q`, :data:`R`).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_QR_update (gsl_matrix * Q, gsl_matrix * R, gsl_vector * w, const gsl_vector * v)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function performs a rank-1 update :math:`w v^T` of the :math:`QR`
|
|
Packit |
67cb25 |
decomposition (:data:`Q`, :data:`R`). The update is given by :math:`Q'R' = Q (R + w v^T)`
|
|
Packit |
67cb25 |
where the output matrices :math:`Q` and :math:`R` are also
|
|
Packit |
67cb25 |
orthogonal and right triangular. Note that :data:`w` is destroyed by the
|
|
Packit |
67cb25 |
update.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_R_solve (const gsl_matrix * R, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the triangular system :math:`R x = b` for the
|
|
Packit |
67cb25 |
:math:`N`-by-:math:`N` matrix :data:`R`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_R_svx (const gsl_matrix * R, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the triangular system :math:`R x = b` in-place. On
|
|
Packit |
67cb25 |
input :data:`x` should contain the right-hand side :math:`b`, which is
|
|
Packit |
67cb25 |
replaced by the solution on output.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: QR decomposition with column pivoting
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
QR Decomposition with Column Pivoting
|
|
Packit |
67cb25 |
=====================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The :math:`QR` decomposition of an :math:`M`-by-:math:`N` matrix :math:`A`
|
|
Packit |
67cb25 |
can be extended to the rank deficient case by introducing a column permutation :math:`P`,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: A P = Q R
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The first :math:`r` columns of :math:`Q` form an orthonormal basis
|
|
Packit |
67cb25 |
for the range of :math:`A` for a matrix with column rank :math:`r`. This
|
|
Packit |
67cb25 |
decomposition can also be used to convert the linear system :math:`A x = b`
|
|
Packit |
67cb25 |
into the triangular system :math:`R y = Q^T b, x = P y`, which can be
|
|
Packit |
67cb25 |
solved by back-substitution and permutation. We denote the :math:`QR`
|
|
Packit |
67cb25 |
decomposition with column pivoting by :math:`QRP^T` since :math:`A = Q R P^T`.
|
|
Packit |
67cb25 |
When :math:`A` is rank deficient with :math:`r = {\rm rank}(A)`, the matrix
|
|
Packit |
67cb25 |
:math:`R` can be partitioned as
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
R = \left(
|
|
Packit |
67cb25 |
\begin{matrix}
|
|
Packit |
67cb25 |
R_{11} & R_{12} \\
|
|
Packit |
67cb25 |
0 & R_{22}
|
|
Packit |
67cb25 |
\end{matrix}
|
|
Packit |
67cb25 |
\right) \approx
|
|
Packit |
67cb25 |
\left(
|
|
Packit |
67cb25 |
\begin{matrix}
|
|
Packit |
67cb25 |
R_{11} & R_{12} \\
|
|
Packit |
67cb25 |
0 & 0
|
|
Packit |
67cb25 |
\end{matrix}
|
|
Packit |
67cb25 |
\right)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
R = [ R11 R12 ] =~ [ R11 R12 ]
|
|
Packit |
67cb25 |
[ 0 R22 ] [ 0 0 ]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`R_{11}` is :math:`r`-by-:math:`r` and nonsingular. In this case,
|
|
Packit |
67cb25 |
a *basic* least squares solution for the overdetermined system :math:`A x = b`
|
|
Packit |
67cb25 |
can be obtained as
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
x = P \left(
|
|
Packit |
67cb25 |
\begin{matrix}
|
|
Packit |
67cb25 |
R_{11}^{-1} c_1 \\
|
|
Packit |
67cb25 |
0
|
|
Packit |
67cb25 |
\end{matrix}
|
|
Packit |
67cb25 |
\right)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
x = P [ R11^-1 c1 ]
|
|
Packit |
67cb25 |
[ 0 ]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`c_1` consists of the first :math:`r` elements of :math:`Q^T b`.
|
|
Packit |
67cb25 |
This basic solution is not guaranteed to be the minimum norm solution unless
|
|
Packit |
67cb25 |
:math:`R_{12} = 0` (see :ref:`Complete Orthogonal Decomposition <cod>`).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_QRPT_decomp (gsl_matrix * A, gsl_vector * tau, gsl_permutation * p, int * signum, gsl_vector * norm)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function factorizes the :math:`M`-by-:math:`N` matrix :data:`A` into
|
|
Packit |
67cb25 |
the :math:`QRP^T` decomposition :math:`A = Q R P^T`. On output the
|
|
Packit |
67cb25 |
diagonal and upper triangular part of the input matrix contain the
|
|
Packit |
67cb25 |
matrix :math:`R`. The permutation matrix :math:`P` is stored in the
|
|
Packit |
67cb25 |
permutation :data:`p`. The sign of the permutation is given by
|
|
Packit |
67cb25 |
:data:`signum`. It has the value :math:`(-1)^n`, where :math:`n` is the
|
|
Packit |
67cb25 |
number of interchanges in the permutation. The vector :data:`tau` and the
|
|
Packit |
67cb25 |
columns of the lower triangular part of the matrix :data:`A` contain the
|
|
Packit |
67cb25 |
Householder coefficients and vectors which encode the orthogonal matrix
|
|
Packit |
67cb25 |
:data:`Q`. The vector :data:`tau` must be of length :math:`k=\min(M,N)`. The
|
|
Packit |
67cb25 |
matrix :math:`Q` is related to these components by, :math:`Q = Q_k ... Q_2 Q_1`
|
|
Packit |
67cb25 |
where :math:`Q_i = I - \tau_i v_i v_i^T` and :math:`v_i` is the
|
|
Packit |
67cb25 |
Householder vector
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: v_i = (0,...,1,A(i+1,i),A(i+2,i),...,A(m,i))
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is the same storage scheme
|
|
Packit |
67cb25 |
as used by |lapack|. The vector :data:`norm` is a workspace of length
|
|
Packit |
67cb25 |
:data:`N` used for column pivoting.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The algorithm used to perform the decomposition is Householder QR with
|
|
Packit |
67cb25 |
column pivoting (Golub & Van Loan, "Matrix Computations", Algorithm
|
|
Packit |
67cb25 |
5.4.1).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: 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)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function factorizes the matrix :data:`A` into the decomposition
|
|
Packit |
67cb25 |
:math:`A = Q R P^T` without modifying :data:`A` itself and storing the
|
|
Packit |
67cb25 |
output in the separate matrices :data:`q` and :data:`r`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_QRPT_solve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the square system :math:`A x = b` using the :math:`QRP^T`
|
|
Packit |
67cb25 |
decomposition of :math:`A` held in (:data:`QR`, :data:`tau`, :data:`p`) which must
|
|
Packit |
67cb25 |
have been computed previously by :func:`gsl_linalg_QRPT_decomp`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_QRPT_svx (const gsl_matrix * QR, const gsl_vector * tau, const gsl_permutation * p, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the square system :math:`A x = b` in-place using the
|
|
Packit |
67cb25 |
:math:`QRP^T` decomposition of :math:`A` held in
|
|
Packit |
67cb25 |
(:data:`QR`, :data:`tau`, :data:`p`). On input :data:`x` should contain the
|
|
Packit |
67cb25 |
right-hand side :math:`b`, which is replaced by the solution on output.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: 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)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function finds the least squares solution to the overdetermined
|
|
Packit |
67cb25 |
system :math:`A x = b` where the matrix :data:`A` has more rows than
|
|
Packit |
67cb25 |
columns and is assumed to have full rank. The least squares solution minimizes
|
|
Packit |
67cb25 |
the Euclidean norm of the residual, :math:`||b - A x||`. The routine requires as input
|
|
Packit |
67cb25 |
the :math:`QR` decomposition of :math:`A` into (:data:`QR`, :data:`tau`, :data:`p`) given by
|
|
Packit |
67cb25 |
:func:`gsl_linalg_QRPT_decomp`. The solution is returned in :data:`x`. The
|
|
Packit |
67cb25 |
residual is computed as a by-product and stored in :data:`residual`. For rank
|
|
Packit |
67cb25 |
deficient matrices, :func:`gsl_linalg_QRPT_lssolve2` should be used instead.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: 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)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function finds the least squares solution to the overdetermined
|
|
Packit |
67cb25 |
system :math:`A x = b` where the matrix :data:`A` has more rows than
|
|
Packit |
67cb25 |
columns and has rank given by the input :data:`rank`. If the user does not
|
|
Packit |
67cb25 |
know the rank of :math:`A`, the routine :func:`gsl_linalg_QRPT_rank` can be
|
|
Packit |
67cb25 |
called to estimate it. The least squares solution is
|
|
Packit |
67cb25 |
the so-called "basic" solution discussed above and may not be the minimum
|
|
Packit |
67cb25 |
norm solution. The routine requires as input
|
|
Packit |
67cb25 |
the :math:`QR` decomposition of :math:`A` into (:data:`QR`, :data:`tau`, :data:`p`) given by
|
|
Packit |
67cb25 |
:func:`gsl_linalg_QRPT_decomp`. The solution is returned in :data:`x`. The
|
|
Packit |
67cb25 |
residual is computed as a by-product and stored in :data:`residual`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_QRPT_QRsolve (const gsl_matrix * Q, const gsl_matrix * R, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the square system :math:`R P^T x = Q^T b` for
|
|
Packit |
67cb25 |
:data:`x`. It can be used when the :math:`QR` decomposition of a matrix is
|
|
Packit |
67cb25 |
available in unpacked form as (:data:`Q`, :data:`R`).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_QRPT_update (gsl_matrix * Q, gsl_matrix * R, const gsl_permutation * p, gsl_vector * w, const gsl_vector * v)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function performs a rank-1 update :math:`w v^T` of the :math:`QRP^T`
|
|
Packit |
67cb25 |
decomposition (:data:`Q`, :data:`R`, :data:`p`). The update is given by
|
|
Packit |
67cb25 |
:math:`Q'R' = Q (R + w v^T P)` where the output matrices :math:`Q'` and
|
|
Packit |
67cb25 |
:math:`R'` are also orthogonal and right triangular. Note that :data:`w` is
|
|
Packit |
67cb25 |
destroyed by the update. The permutation :data:`p` is not changed.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_QRPT_Rsolve (const gsl_matrix * QR, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the triangular system :math:`R P^T x = b` for the
|
|
Packit |
67cb25 |
:math:`N`-by-:math:`N` matrix :math:`R` contained in :data:`QR`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_QRPT_Rsvx (const gsl_matrix * QR, const gsl_permutation * p, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the triangular system :math:`R P^T x = b` in-place
|
|
Packit |
67cb25 |
for the :math:`N`-by-:math:`N` matrix :math:`R` contained in :data:`QR`. On
|
|
Packit |
67cb25 |
input :data:`x` should contain the right-hand side :math:`b`, which is
|
|
Packit |
67cb25 |
replaced by the solution on output.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: size_t gsl_linalg_QRPT_rank (const gsl_matrix * QR, const double tol)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function estimates the rank of the triangular matrix :math:`R` contained in :data:`QR`.
|
|
Packit |
67cb25 |
The algorithm simply counts the number of diagonal elements of :math:`R` whose absolute value
|
|
Packit |
67cb25 |
is greater than the specified tolerance :data:`tol`. If the input :data:`tol` is negative,
|
|
Packit |
67cb25 |
a default value of :math:`20 (M + N) eps(max(|diag(R)|))` is used.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_QRPT_rcond (const gsl_matrix * QR, double * rcond, gsl_vector * work)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function estimates the reciprocal condition number (using the 1-norm) of the :math:`R` factor,
|
|
Packit |
67cb25 |
stored in the upper triangle of :data:`QR`. The reciprocal condition number estimate, defined as
|
|
Packit |
67cb25 |
:math:`1 / (||R||_1 \cdot ||R^{-1}||_1)`, is stored in :data:`rcond`.
|
|
Packit |
67cb25 |
Additional workspace of size :math:`3 N` is required in :data:`work`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: complete orthogonal decomposition
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. _cod:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Complete Orthogonal Decomposition
|
|
Packit |
67cb25 |
=================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The complete orthogonal decomposition of a :math:`M`-by-:math:`N` matrix
|
|
Packit |
67cb25 |
:math:`A` is a generalization of the QR decomposition with column pivoting, given by
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A P = Q
|
|
Packit |
67cb25 |
\left(
|
|
Packit |
67cb25 |
\begin{matrix}
|
|
Packit |
67cb25 |
R_{11} & 0 \\
|
|
Packit |
67cb25 |
0 & 0
|
|
Packit |
67cb25 |
\end{matrix}
|
|
Packit |
67cb25 |
\right) Z^T
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A P = Q [ R11 0 ] Z^T
|
|
Packit |
67cb25 |
[ 0 0 ]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`P` is a :math:`N`-by-:math:`N` permutation matrix,
|
|
Packit |
67cb25 |
:math:`Q` is :math:`M`-by-:math:`M` orthogonal, :math:`R_{11}` is
|
|
Packit |
67cb25 |
:math:`r`-by-:math:`r` upper triangular, with :math:`r = {\rm rank}(A)`,
|
|
Packit |
67cb25 |
and :math:`Z` is :math:`N`-by-:math:`N` orthogonal. If :math:`A`
|
|
Packit |
67cb25 |
has full rank, then :math:`R_{11} = R`, :math:`Z = I` and this reduces to the
|
|
Packit |
67cb25 |
QR decomposition with column pivoting.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
For a rank deficient least squares problem, :math:`\min_x{|| b - Ax||^2}`, the solution vector
|
|
Packit |
67cb25 |
:math:`x` is not unique. However if we further require that :math:`||x||^2` is minimized,
|
|
Packit |
67cb25 |
then the complete orthogonal decomposition gives us the ability to compute
|
|
Packit |
67cb25 |
the unique minimum norm solution, which is given by
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
x = P Z
|
|
Packit |
67cb25 |
\left(
|
|
Packit |
67cb25 |
\begin{matrix}
|
|
Packit |
67cb25 |
R_{11}^{-1} c_1 \\
|
|
Packit |
67cb25 |
0
|
|
Packit |
67cb25 |
\end{matrix}
|
|
Packit |
67cb25 |
\right)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
x = P Z [ R11^-1 c1 ]
|
|
Packit |
67cb25 |
[ 0 ]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
and the vector :math:`c_1` is the first :math:`r` elements of :math:`Q^T b`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The COD also enables a straightforward solution of regularized least squares problems
|
|
Packit |
67cb25 |
in Tikhonov standard form, written as
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: \min_x ||b - A x||^2 + \lambda^2 ||x||^2
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`\lambda > 0` is a regularization parameter which represents a tradeoff between
|
|
Packit |
67cb25 |
minimizing the residual norm :math:`||b - A x||` and the solution norm :math:`||x||`. For this system,
|
|
Packit |
67cb25 |
the solution is given by
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
x = P Z
|
|
Packit |
67cb25 |
\left(
|
|
Packit |
67cb25 |
\begin{matrix}
|
|
Packit |
67cb25 |
y_1 \\
|
|
Packit |
67cb25 |
0
|
|
Packit |
67cb25 |
\end{matrix}
|
|
Packit |
67cb25 |
\right)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
x = P Z [ y1 ]
|
|
Packit |
67cb25 |
[ 0 ]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`y_1` is a vector of length :math:`r` which is found by solving
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
\left(
|
|
Packit |
67cb25 |
\begin{matrix}
|
|
Packit |
67cb25 |
R_{11} \\
|
|
Packit |
67cb25 |
\lambda I_r
|
|
Packit |
67cb25 |
\end{matrix}
|
|
Packit |
67cb25 |
\right) y_1 =
|
|
Packit |
67cb25 |
\left(
|
|
Packit |
67cb25 |
\begin{matrix}
|
|
Packit |
67cb25 |
c_1 \\
|
|
Packit |
67cb25 |
0
|
|
Packit |
67cb25 |
\end{matrix}
|
|
Packit |
67cb25 |
\right)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
[ R11 ] y_1 = [ c_1 ]
|
|
Packit |
67cb25 |
[ \lambda I_r ] [ 0 ]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
and :math:`c_1` is defined above. The equation above can be solved efficiently for different
|
|
Packit |
67cb25 |
values of :math:`\lambda` using QR factorizations of the left hand side matrix.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: 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)
|
|
Packit |
67cb25 |
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)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions factor the :math:`M`-by-:math:`N` matrix :data:`A` into the decomposition :math:`A = Q R Z P^T`. The rank of :data:`A`
|
|
Packit |
67cb25 |
is computed as the number of diagonal elements of :math:`R` greater than the tolerance :data:`tol` and output in :data:`rank`.
|
|
Packit |
67cb25 |
If :data:`tol` is not specified, a default value is used (see :func:`gsl_linalg_QRPT_rank`). On output, the permutation
|
|
Packit |
67cb25 |
matrix :math:`P` is stored in :data:`p`. The matrix :math:`R_{11}` is stored in the upper :data:`rank`-by-:data:`rank` block of :data:`A`.
|
|
Packit |
67cb25 |
The matrices :math:`Q` and :math:`Z` are encoded in packed storage in :data:`A` on output. The vectors :data:`tau_Q` and :data:`tau_Z`
|
|
Packit |
67cb25 |
contain the Householder scalars corresponding to the matrices :math:`Q` and :math:`Z` respectively and must be
|
|
Packit |
67cb25 |
of length :math:`k = \min(M,N)`. The vector :data:`work` is additional workspace of length :math:`N`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: 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)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function finds the unique minimum norm least squares solution to the overdetermined
|
|
Packit |
67cb25 |
system :math:`A x = b` where the matrix :data:`A` has more rows than
|
|
Packit |
67cb25 |
columns. The least squares solution minimizes the Euclidean norm of the
|
|
Packit |
67cb25 |
residual, :math:`||b - A x||` as well as the norm of the solution :math:`||x||`. The routine requires as input
|
|
Packit |
67cb25 |
the :math:`QRZT` decomposition of :math:`A` into (:data:`QRZT`, :data:`tau_Q`, :data:`tau_Z`, :data:`p`, :data:`rank`)
|
|
Packit |
67cb25 |
given by :func:`gsl_linalg_COD_decomp`. The solution is returned in :data:`x`. The
|
|
Packit |
67cb25 |
residual, :math:`b - Ax`, is computed as a by-product and stored in :data:`residual`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: 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)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function finds the solution to the regularized least squares problem in Tikhonov
|
|
Packit |
67cb25 |
standard form, :math:`\min_x ||b - Ax||^2 + \lambda^2 ||x||^2`. The routine requires as input
|
|
Packit |
67cb25 |
the :math:`QRZT` decomposition of :math:`A` into (:data:`QRZT`, :data:`tau_Q`, :data:`tau_Z`, :data:`p`, :data:`rank`)
|
|
Packit |
67cb25 |
given by :func:`gsl_linalg_COD_decomp`. The parameter :math:`\lambda` is supplied in :data:`lambda`. The solution
|
|
Packit |
67cb25 |
is returned in :data:`x`. The residual, :math:`b - Ax`, is stored in :data:`residual` on output. :data:`S` is additional
|
|
Packit |
67cb25 |
workspace of size :data:`rank`-by-:data:`rank`. :data:`work` is additional workspace of length :data:`rank`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: 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)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function unpacks the encoded :math:`QRZT` decomposition
|
|
Packit |
67cb25 |
(:data:`QRZT`, :data:`tau_Q`, :data:`tau_Z`, :data:`rank`) into the matrices
|
|
Packit |
67cb25 |
:data:`Q`, :data:`R`, and :data:`Z`, where :data:`Q` is :math:`M`-by-:math:`M`,
|
|
Packit |
67cb25 |
:data:`R` is :math:`M`-by-:math:`N`, and :data:`Z` is :math:`N`-by-:math:`N`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_COD_matZ (const gsl_matrix * QRZT, const gsl_vector * tau_Z, const size_t rank, gsl_matrix * A, gsl_vector * work)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function multiplies the input matrix :data:`A` on the right by :data:`Z`,
|
|
Packit |
67cb25 |
:math:`A' = A Z` using the encoded :math:`QRZT` decomposition
|
|
Packit |
67cb25 |
(:data:`QRZT`, :data:`tau_Z`, :data:`rank`). :data:`A` must have :math:`N` columns but may
|
|
Packit |
67cb25 |
have any number of rows. Additional workspace of length :math:`M` is provided
|
|
Packit |
67cb25 |
in :data:`work`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: SVD, singular value decomposition
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Singular Value Decomposition
|
|
Packit |
67cb25 |
============================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A general rectangular :math:`M`-by-:math:`N` matrix :math:`A` has a
|
|
Packit |
67cb25 |
singular value decomposition (SVD) into the product of an
|
|
Packit |
67cb25 |
:math:`M`-by-:math:`N` orthogonal matrix :math:`U`, an :math:`N`-by-:math:`N`
|
|
Packit |
67cb25 |
diagonal matrix of singular values :math:`S` and the transpose of an
|
|
Packit |
67cb25 |
:math:`N`-by-:math:`N` orthogonal square matrix :math:`V`,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: A = U S V^T
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The singular values :math:`\sigma_i = S_{ii}`
|
|
Packit |
67cb25 |
are all non-negative and are
|
|
Packit |
67cb25 |
generally chosen to form a non-increasing sequence
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: \sigma_1 \ge \sigma_2 \ge ... \ge \sigma_N \ge 0
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: \sigma_1 >= \sigma_2 >= ... >= \sigma_N >= 0
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The singular value decomposition of a matrix has many practical uses.
|
|
Packit |
67cb25 |
The condition number of the matrix is given by the ratio of the largest
|
|
Packit |
67cb25 |
singular value to the smallest singular value. The presence of a zero
|
|
Packit |
67cb25 |
singular value indicates that the matrix is singular. The number of
|
|
Packit |
67cb25 |
non-zero singular values indicates the rank of the matrix. In practice
|
|
Packit |
67cb25 |
singular value decomposition of a rank-deficient matrix will not produce
|
|
Packit |
67cb25 |
exact zeroes for singular values, due to finite numerical
|
|
Packit |
67cb25 |
precision. Small singular values should be edited by choosing a suitable
|
|
Packit |
67cb25 |
tolerance.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
For a rank-deficient matrix, the null space of :math:`A` is given by
|
|
Packit |
67cb25 |
the columns of :math:`V` corresponding to the zero singular values.
|
|
Packit |
67cb25 |
Similarly, the range of :math:`A` is given by columns of :math:`U`
|
|
Packit |
67cb25 |
corresponding to the non-zero singular values.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Note that the routines here compute the "thin" version of the SVD
|
|
Packit |
67cb25 |
with :math:`U` as :math:`M`-by-:math:`N` orthogonal matrix. This allows
|
|
Packit |
67cb25 |
in-place computation and is the most commonly-used form in practice.
|
|
Packit |
67cb25 |
Mathematically, the "full" SVD is defined with :math:`U` as an
|
|
Packit |
67cb25 |
:math:`M`-by-:math:`M` orthogonal matrix and :math:`S` as an
|
|
Packit |
67cb25 |
:math:`M`-by-:math:`N` diagonal matrix (with additional rows of zeros).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_SV_decomp (gsl_matrix * A, gsl_matrix * V, gsl_vector * S, gsl_vector * work)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function factorizes the :math:`M`-by-:math:`N` matrix :data:`A` into
|
|
Packit |
67cb25 |
the singular value decomposition :math:`A = U S V^T` for :math:`M \ge N`.
|
|
Packit |
67cb25 |
On output the matrix :data:`A` is replaced by
|
|
Packit |
67cb25 |
:math:`U`. The diagonal elements of the singular value matrix :math:`S`
|
|
Packit |
67cb25 |
are stored in the vector :data:`S`. The singular values are non-negative
|
|
Packit |
67cb25 |
and form a non-increasing sequence from :math:`S_1` to :math:`S_N`. The
|
|
Packit |
67cb25 |
matrix :data:`V` contains the elements of :math:`V` in untransposed
|
|
Packit |
67cb25 |
form. To form the product :math:`U S V^T` it is necessary to take the
|
|
Packit |
67cb25 |
transpose of :data:`V`. A workspace of length :data:`N` is required in
|
|
Packit |
67cb25 |
:data:`work`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This routine uses the Golub-Reinsch SVD algorithm.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_SV_decomp_mod (gsl_matrix * A, gsl_matrix * X, gsl_matrix * V, gsl_vector * S, gsl_vector * work)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function computes the SVD using the modified Golub-Reinsch
|
|
Packit |
67cb25 |
algorithm, which is faster for :math:`M \gg N`.
|
|
Packit |
67cb25 |
It requires the vector :data:`work` of length :data:`N` and the
|
|
Packit |
67cb25 |
:math:`N`-by-:math:`N` matrix :data:`X` as additional working space.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: Jacobi orthogonalization
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_SV_decomp_jacobi (gsl_matrix * A, gsl_matrix * V, gsl_vector * S)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function computes the SVD of the :math:`M`-by-:math:`N` matrix :data:`A`
|
|
Packit |
67cb25 |
using one-sided Jacobi orthogonalization for :math:`M \ge N`.
|
|
Packit |
67cb25 |
The Jacobi method can compute singular values to higher
|
|
Packit |
67cb25 |
relative accuracy than Golub-Reinsch algorithms (see references for
|
|
Packit |
67cb25 |
details).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_SV_solve (const gsl_matrix * U, const gsl_matrix * V, const gsl_vector * S, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the system :math:`A x = b` using the singular value
|
|
Packit |
67cb25 |
decomposition (:data:`U`, :data:`S`, :data:`V`) of :math:`A` which must
|
|
Packit |
67cb25 |
have been computed previously with :func:`gsl_linalg_SV_decomp`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Only non-zero singular values are used in computing the solution. The
|
|
Packit |
67cb25 |
parts of the solution corresponding to singular values of zero are
|
|
Packit |
67cb25 |
ignored. Other singular values can be edited out by setting them to
|
|
Packit |
67cb25 |
zero before calling this function.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
In the over-determined case where :data:`A` has more rows than columns the
|
|
Packit |
67cb25 |
system is solved in the least squares sense, returning the solution
|
|
Packit |
67cb25 |
:data:`x` which minimizes :math:`||A x - b||_2`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_SV_leverage (const gsl_matrix * U, gsl_vector * h)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function computes the statistical leverage values :math:`h_i` of a matrix :math:`A`
|
|
Packit |
67cb25 |
using its singular value decomposition (:data:`U`, :data:`S`, :data:`V`) previously computed
|
|
Packit |
67cb25 |
with :func:`gsl_linalg_SV_decomp`. :math:`h_i` are the diagonal values of the matrix
|
|
Packit |
67cb25 |
:math:`A (A^T A)^{-1} A^T` and depend only on the matrix :data:`U` which is the input to
|
|
Packit |
67cb25 |
this function.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: Cholesky decomposition
|
|
Packit |
67cb25 |
single: square root of a matrix, Cholesky decomposition
|
|
Packit |
67cb25 |
single: matrix square root, Cholesky decomposition
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. _sec_cholesky-decomposition:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Cholesky Decomposition
|
|
Packit |
67cb25 |
======================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A symmetric, positive definite square matrix :math:`A` has a Cholesky
|
|
Packit |
67cb25 |
decomposition into a product of a lower triangular matrix :math:`L` and
|
|
Packit |
67cb25 |
its transpose :math:`L^T`,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: A = L L^T
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is sometimes referred to as taking the square-root of a matrix. The
|
|
Packit |
67cb25 |
Cholesky decomposition can only be carried out when all the eigenvalues
|
|
Packit |
67cb25 |
of the matrix are positive. This decomposition can be used to convert
|
|
Packit |
67cb25 |
the linear system :math:`A x = b` into a pair of triangular systems
|
|
Packit |
67cb25 |
(:math:`L y = b`, :math:`L^T x = y`), which can be solved by forward and
|
|
Packit |
67cb25 |
back-substitution.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
If the matrix :math:`A` is near singular, it is sometimes possible to reduce
|
|
Packit |
67cb25 |
the condition number and recover a more accurate solution vector :math:`x`
|
|
Packit |
67cb25 |
by scaling as
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: \left( S A S \right) \left( S^{-1} x \right) = S b
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: ( S A S ) ( S^(-1) x ) = S b
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`S` is a diagonal matrix whose elements are given by
|
|
Packit |
67cb25 |
:math:`S_{ii} = 1/\sqrt{A_{ii}}`. This scaling is also known as
|
|
Packit |
67cb25 |
Jacobi preconditioning. There are routines below to solve
|
|
Packit |
67cb25 |
both the scaled and unscaled systems.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_cholesky_decomp1 (gsl_matrix * A)
|
|
Packit |
67cb25 |
int gsl_linalg_complex_cholesky_decomp (gsl_matrix_complex * A)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions factorize the symmetric, positive-definite square matrix
|
|
Packit |
67cb25 |
:data:`A` into the Cholesky decomposition :math:`A = L L^T` (or
|
|
Packit |
67cb25 |
:math:`A = L L^{\dagger}`
|
|
Packit |
67cb25 |
for the complex case). On input, the values from the diagonal and lower-triangular
|
|
Packit |
67cb25 |
part of the matrix :data:`A` are used (the upper triangular part is ignored). On output
|
|
Packit |
67cb25 |
the diagonal and lower triangular part of the input matrix :data:`A` contain the matrix
|
|
Packit |
67cb25 |
:math:`L`, while the upper triangular part is unmodified. If the matrix is not
|
|
Packit |
67cb25 |
positive-definite then the decomposition will fail, returning the
|
|
Packit |
67cb25 |
error code :macro:`GSL_EDOM`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
When testing whether a matrix is positive-definite, disable the error
|
|
Packit |
67cb25 |
handler first to avoid triggering an error.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_cholesky_decomp (gsl_matrix * A)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function is now deprecated and is provided only for backward compatibility.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_cholesky_solve (const gsl_matrix * cholesky, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
int gsl_linalg_complex_cholesky_solve (const gsl_matrix_complex * cholesky, const gsl_vector_complex * b, gsl_vector_complex * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions solve the system :math:`A x = b` using the Cholesky
|
|
Packit |
67cb25 |
decomposition of :math:`A` held in the matrix :data:`cholesky` which must
|
|
Packit |
67cb25 |
have been previously computed by :func:`gsl_linalg_cholesky_decomp` or
|
|
Packit |
67cb25 |
:func:`gsl_linalg_complex_cholesky_decomp`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_cholesky_svx (const gsl_matrix * cholesky, gsl_vector * x)
|
|
Packit |
67cb25 |
int gsl_linalg_complex_cholesky_svx (const gsl_matrix_complex * cholesky, gsl_vector_complex * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions solve the system :math:`A x = b` in-place using the
|
|
Packit |
67cb25 |
Cholesky decomposition of :math:`A` held in the matrix :data:`cholesky`
|
|
Packit |
67cb25 |
which must have been previously computed by
|
|
Packit |
67cb25 |
:func:`gsl_linalg_cholesky_decomp` or
|
|
Packit |
67cb25 |
:func:`gsl_linalg_complex_cholesky_decomp`. On input :data:`x` should
|
|
Packit |
67cb25 |
contain the right-hand side :math:`b`, which is replaced by the
|
|
Packit |
67cb25 |
solution on output.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_cholesky_invert (gsl_matrix * cholesky)
|
|
Packit |
67cb25 |
int gsl_linalg_complex_cholesky_invert (gsl_matrix_complex * cholesky)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions compute the inverse of a matrix from its Cholesky
|
|
Packit |
67cb25 |
decomposition :data:`cholesky`, which must have been previously computed
|
|
Packit |
67cb25 |
by :func:`gsl_linalg_cholesky_decomp` or
|
|
Packit |
67cb25 |
:func:`gsl_linalg_complex_cholesky_decomp`. On output, the inverse is
|
|
Packit |
67cb25 |
stored in-place in :data:`cholesky`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_cholesky_decomp2 (gsl_matrix * A, gsl_vector * S)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function calculates a diagonal scaling transformation :math:`S` for
|
|
Packit |
67cb25 |
the symmetric, positive-definite square matrix :data:`A`, and then
|
|
Packit |
67cb25 |
computes the Cholesky decomposition :math:`S A S = L L^T`.
|
|
Packit |
67cb25 |
On input, the values from the diagonal and lower-triangular part of the matrix :data:`A` are
|
|
Packit |
67cb25 |
used (the upper triangular part is ignored). On output the diagonal and lower triangular part
|
|
Packit |
67cb25 |
of the input matrix :data:`A` contain the matrix :math:`L`, while the upper triangular part
|
|
Packit |
67cb25 |
of the input matrix is overwritten with :math:`L^T` (the diagonal terms being
|
|
Packit |
67cb25 |
identical for both :math:`L` and :math:`L^T`). If the matrix is not
|
|
Packit |
67cb25 |
positive-definite then the decomposition will fail, returning the
|
|
Packit |
67cb25 |
error code :macro:`GSL_EDOM`. The diagonal scale factors are stored in :data:`S`
|
|
Packit |
67cb25 |
on output.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
When testing whether a matrix is positive-definite, disable the error
|
|
Packit |
67cb25 |
handler first to avoid triggering an error.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_cholesky_solve2 (const gsl_matrix * cholesky, const gsl_vector * S, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the system :math:`(S A S) (S^{-1} x) = S b` using the Cholesky
|
|
Packit |
67cb25 |
decomposition of :math:`S A S` held in the matrix :data:`cholesky` which must
|
|
Packit |
67cb25 |
have been previously computed by :func:`gsl_linalg_cholesky_decomp2`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_cholesky_svx2 (const gsl_matrix * cholesky, const gsl_vector * S, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the system :math:`(S A S) (S^{-1} x) = S b` in-place using the
|
|
Packit |
67cb25 |
Cholesky decomposition of :math:`S A S` held in the matrix :data:`cholesky`
|
|
Packit |
67cb25 |
which must have been previously computed by
|
|
Packit |
67cb25 |
:func:`gsl_linalg_cholesky_decomp2`. On input :data:`x` should
|
|
Packit |
67cb25 |
contain the right-hand side :math:`b`, which is replaced by the
|
|
Packit |
67cb25 |
solution on output.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_cholesky_scale (const gsl_matrix * A, gsl_vector * S)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function calculates a diagonal scaling transformation of the
|
|
Packit |
67cb25 |
symmetric, positive definite matrix :data:`A`, such that
|
|
Packit |
67cb25 |
:math:`S A S` has a condition number within a factor of :math:`N`
|
|
Packit |
67cb25 |
of the matrix of smallest possible condition number over all
|
|
Packit |
67cb25 |
possible diagonal scalings. On output, :data:`S` contains the
|
|
Packit |
67cb25 |
scale factors, given by :math:`S_i = 1/\sqrt{A_{ii}}`.
|
|
Packit |
67cb25 |
For any :math:`A_{ii} \le 0`, the corresponding scale factor :math:`S_i`
|
|
Packit |
67cb25 |
is set to :math:`1`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_cholesky_scale_apply (gsl_matrix * A, const gsl_vector * S)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function applies the scaling transformation :data:`S` to the matrix :data:`A`. On output,
|
|
Packit |
67cb25 |
:data:`A` is replaced by :math:`S A S`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_cholesky_rcond (const gsl_matrix * cholesky, double * rcond, gsl_vector * work)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function estimates the reciprocal condition number (using the 1-norm) of the symmetric positive
|
|
Packit |
67cb25 |
definite matrix :math:`A`, using its Cholesky decomposition provided in :data:`cholesky`.
|
|
Packit |
67cb25 |
The reciprocal condition number estimate, defined as :math:`1 / (||A||_1 \cdot ||A^{-1}||_1)`, is stored
|
|
Packit |
67cb25 |
in :data:`rcond`. Additional workspace of size :math:`3 N` is required in :data:`work`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: Cholesky decomposition, pivoted
|
|
Packit |
67cb25 |
single: Pivoted Cholesky Decomposition
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Pivoted Cholesky Decomposition
|
|
Packit |
67cb25 |
==============================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A symmetric, positive definite square matrix :math:`A` has an alternate
|
|
Packit |
67cb25 |
Cholesky decomposition into a product of a lower unit triangular matrix :math:`L`,
|
|
Packit |
67cb25 |
a diagonal matrix :math:`D` and :math:`L^T`, given by :math:`L D L^T`. This is
|
|
Packit |
67cb25 |
equivalent to the Cholesky formulation discussed above, with
|
|
Packit |
67cb25 |
the standard Cholesky lower triangular factor given by :math:`L D^{1 \over 2}`.
|
|
Packit |
67cb25 |
For ill-conditioned matrices, it can help to use a pivoting strategy to
|
|
Packit |
67cb25 |
prevent the entries of :math:`D` and :math:`L` from growing too large, and also
|
|
Packit |
67cb25 |
ensure :math:`D_1 \ge D_2 \ge \cdots \ge D_n > 0`, where :math:`D_i` are
|
|
Packit |
67cb25 |
the diagonal entries of :math:`D`. The final decomposition is given by
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: P A P^T = L D L^T
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`P` is a permutation matrix.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_pcholesky_decomp (gsl_matrix * A, gsl_permutation * p)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function factors the symmetric, positive-definite square matrix
|
|
Packit |
67cb25 |
:data:`A` into the Pivoted Cholesky decomposition :math:`P A P^T = L D L^T`.
|
|
Packit |
67cb25 |
On input, the values from the diagonal and lower-triangular part of the matrix :data:`A` are
|
|
Packit |
67cb25 |
used to construct the factorization. On output the diagonal of the input matrix :data:`A` stores
|
|
Packit |
67cb25 |
the diagonal elements of :math:`D`, and the lower triangular portion of :data:`A`
|
|
Packit |
67cb25 |
contains the matrix :math:`L`. Since :math:`L` has ones on its diagonal these
|
|
Packit |
67cb25 |
do not need to be explicitely stored. The upper triangular portion of :data:`A` is
|
|
Packit |
67cb25 |
unmodified. The permutation matrix :math:`P` is stored in :data:`p` on output.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_pcholesky_solve (const gsl_matrix * LDLT, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the system :math:`A x = b` using the Pivoted Cholesky
|
|
Packit |
67cb25 |
decomposition of :math:`A` held in the matrix :data:`LDLT` and permutation
|
|
Packit |
67cb25 |
:data:`p` which must have been previously computed by :func:`gsl_linalg_pcholesky_decomp`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_pcholesky_svx (const gsl_matrix * LDLT, const gsl_permutation * p, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the system :math:`A x = b` in-place using the Pivoted Cholesky
|
|
Packit |
67cb25 |
decomposition of :math:`A` held in the matrix :data:`LDLT` and permutation
|
|
Packit |
67cb25 |
:data:`p` which must have been previously computed by :func:`gsl_linalg_pcholesky_decomp`.
|
|
Packit |
67cb25 |
On input, :data:`x` contains the right hand side vector :math:`b` which is
|
|
Packit |
67cb25 |
replaced by the solution vector on output.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_pcholesky_decomp2 (gsl_matrix * A, gsl_permutation * p, gsl_vector * S)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function computes the pivoted Cholesky factorization of the matrix
|
|
Packit |
67cb25 |
:math:`S A S`, where the input matrix :data:`A` is symmetric and positive
|
|
Packit |
67cb25 |
definite, and the diagonal scaling matrix :data:`S` is computed to reduce the
|
|
Packit |
67cb25 |
condition number of :data:`A` as much as possible. See
|
|
Packit |
67cb25 |
:ref:`Cholesky Decomposition <sec_cholesky-decomposition>` for more information on the matrix :data:`S`.
|
|
Packit |
67cb25 |
The Pivoted Cholesky decomposition satisfies :math:`P S A S P^T = L D L^T`.
|
|
Packit |
67cb25 |
On input, the values from the diagonal and lower-triangular part of the matrix :data:`A` are
|
|
Packit |
67cb25 |
used to construct the factorization. On output the diagonal of the input matrix :data:`A` stores the diagonal
|
|
Packit |
67cb25 |
elements of :math:`D`, and the lower triangular portion of :data:`A`
|
|
Packit |
67cb25 |
contains the matrix :math:`L`. Since :math:`L` has ones on its diagonal these
|
|
Packit |
67cb25 |
do not need to be explicitely stored. The upper triangular portion of :data:`A`
|
|
Packit |
67cb25 |
is unmodified. The permutation matrix :math:`P` is stored in :data:`p` on output.
|
|
Packit |
67cb25 |
The diagonal scaling transformation is stored in :data:`S` on output.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_pcholesky_solve2 (const gsl_matrix * LDLT, const gsl_permutation * p, const gsl_vector * S, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the system :math:`(S A S) (S^{-1} x) = S b` using the Pivoted Cholesky
|
|
Packit |
67cb25 |
decomposition of :math:`S A S` held in the matrix :data:`LDLT`, permutation
|
|
Packit |
67cb25 |
:data:`p`, and vector :data:`S`, which must have been previously computed by
|
|
Packit |
67cb25 |
:func:`gsl_linalg_pcholesky_decomp2`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_pcholesky_svx2 (const gsl_matrix * LDLT, const gsl_permutation * p, const gsl_vector * S, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the system :math:`(S A S) (S^{-1} x) = S b` in-place using the Pivoted Cholesky
|
|
Packit |
67cb25 |
decomposition of :math:`S A S` held in the matrix :data:`LDLT`, permutation
|
|
Packit |
67cb25 |
:data:`p` and vector :data:`S`, which must have been previously computed by
|
|
Packit |
67cb25 |
:func:`gsl_linalg_pcholesky_decomp2`.
|
|
Packit |
67cb25 |
On input, :data:`x` contains the right hand side vector :math:`b` which is
|
|
Packit |
67cb25 |
replaced by the solution vector on output.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_pcholesky_invert (const gsl_matrix * LDLT, const gsl_permutation * p, gsl_matrix * Ainv)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function computes the inverse of the matrix :math:`A`, using the Pivoted
|
|
Packit |
67cb25 |
Cholesky decomposition stored in :data:`LDLT` and :data:`p`. On output, the
|
|
Packit |
67cb25 |
matrix :data:`Ainv` contains :math:`A^{-1}`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_pcholesky_rcond (const gsl_matrix * LDLT, const gsl_permutation * p, double * rcond, gsl_vector * work)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function estimates the reciprocal condition number (using the 1-norm) of the symmetric positive
|
|
Packit |
67cb25 |
definite matrix :math:`A`, using its pivoted Cholesky decomposition provided in :data:`LDLT`.
|
|
Packit |
67cb25 |
The reciprocal condition number estimate, defined as :math:`1 / (||A||_1 \cdot ||A^{-1}||_1)`, is stored
|
|
Packit |
67cb25 |
in :data:`rcond`. Additional workspace of size :math:`3 N` is required in :data:`work`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: Cholesky decomposition, modified
|
|
Packit |
67cb25 |
single: Modified Cholesky Decomposition
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Modified Cholesky Decomposition
|
|
Packit |
67cb25 |
===============================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The modified Cholesky decomposition is suitable for solving systems
|
|
Packit |
67cb25 |
:math:`A x = b` where :math:`A` is a symmetric indefinite matrix. Such
|
|
Packit |
67cb25 |
matrices arise in nonlinear optimization algorithms. The standard
|
|
Packit |
67cb25 |
Cholesky decomposition requires a positive definite matrix and would
|
|
Packit |
67cb25 |
fail in this case. Instead of resorting to a method like QR or SVD,
|
|
Packit |
67cb25 |
which do not take into account the symmetry of the matrix, we can
|
|
Packit |
67cb25 |
instead introduce a small perturbation to the matrix :math:`A` to
|
|
Packit |
67cb25 |
make it positive definite, and then use a Cholesky decomposition on
|
|
Packit |
67cb25 |
the perturbed matrix. The resulting decomposition satisfies
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: P (A + E) P^T = L D L^T
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`P` is a permutation matrix, :math:`E` is a diagonal
|
|
Packit |
67cb25 |
perturbation matrix, :math:`L` is unit lower triangular, and
|
|
Packit |
67cb25 |
:math:`D` is diagonal. If :math:`A` is sufficiently positive
|
|
Packit |
67cb25 |
definite, then the perturbation matrix :math:`E` will be zero
|
|
Packit |
67cb25 |
and this method is equivalent to the pivoted Cholesky algorithm.
|
|
Packit |
67cb25 |
For indefinite matrices, the perturbation matrix :math:`E` is
|
|
Packit |
67cb25 |
computed to ensure that :math:`A + E` is positive definite and
|
|
Packit |
67cb25 |
well conditioned.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_mcholesky_decomp (gsl_matrix * A, gsl_permutation * p, gsl_vector * E)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function factors the symmetric, indefinite square matrix
|
|
Packit |
67cb25 |
:data:`A` into the Modified Cholesky decomposition :math:`P (A + E) P^T = L D L^T`.
|
|
Packit |
67cb25 |
On input, the values from the diagonal and lower-triangular part of the matrix :data:`A` are
|
|
Packit |
67cb25 |
used to construct the factorization. On output the diagonal of the input matrix :data:`A` stores the diagonal
|
|
Packit |
67cb25 |
elements of :math:`D`, and the lower triangular portion of :data:`A`
|
|
Packit |
67cb25 |
contains the matrix :math:`L`. Since :math:`L` has ones on its diagonal these
|
|
Packit |
67cb25 |
do not need to be explicitely stored. The upper triangular portion of :data:`A`
|
|
Packit |
67cb25 |
is unmodified. The permutation matrix :math:`P` is
|
|
Packit |
67cb25 |
stored in :data:`p` on output. The diagonal perturbation matrix is stored in
|
|
Packit |
67cb25 |
:data:`E` on output. The parameter :data:`E` may be set to NULL if it is not
|
|
Packit |
67cb25 |
required.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_mcholesky_solve (const gsl_matrix * LDLT, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the perturbed system :math:`(A + E) x = b` using the Cholesky
|
|
Packit |
67cb25 |
decomposition of :math:`A + E` held in the matrix :data:`LDLT` and permutation
|
|
Packit |
67cb25 |
:data:`p` which must have been previously computed by :func:`gsl_linalg_mcholesky_decomp`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_mcholesky_svx (const gsl_matrix * LDLT, const gsl_permutation * p, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the perturbed system :math:`(A + E) x = b` in-place using the Cholesky
|
|
Packit |
67cb25 |
decomposition of :math:`A + E` held in the matrix :data:`LDLT` and permutation
|
|
Packit |
67cb25 |
:data:`p` which must have been previously computed by :func:`gsl_linalg_mcholesky_decomp`.
|
|
Packit |
67cb25 |
On input, :data:`x` contains the right hand side vector :math:`b` which is
|
|
Packit |
67cb25 |
replaced by the solution vector on output.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_mcholesky_rcond (const gsl_matrix * LDLT, const gsl_permutation * p, double * rcond, gsl_vector * work)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function estimates the reciprocal condition number (using the 1-norm) of the perturbed matrix
|
|
Packit |
67cb25 |
:math:`A + E`, using its pivoted Cholesky decomposition provided in :data:`LDLT`.
|
|
Packit |
67cb25 |
The reciprocal condition number estimate, defined as :math:`1 / (||A + E||_1 \cdot ||(A + E)^{-1}||_1)`, is stored
|
|
Packit |
67cb25 |
in :data:`rcond`. Additional workspace of size :math:`3 N` is required in :data:`work`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: tridiagonal decomposition
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Tridiagonal Decomposition of Real Symmetric Matrices
|
|
Packit |
67cb25 |
====================================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A symmetric matrix :math:`A` can be factorized by similarity
|
|
Packit |
67cb25 |
transformations into the form,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: A = Q T Q^T
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`Q` is an orthogonal matrix and :math:`T` is a symmetric
|
|
Packit |
67cb25 |
tridiagonal matrix.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_symmtd_decomp (gsl_matrix * A, gsl_vector * tau)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function factorizes the symmetric square matrix :data:`A` into the
|
|
Packit |
67cb25 |
symmetric tridiagonal decomposition :math:`Q T Q^T`. On output the
|
|
Packit |
67cb25 |
diagonal and subdiagonal part of the input matrix :data:`A` contain the
|
|
Packit |
67cb25 |
tridiagonal matrix :math:`T`. The remaining lower triangular part of the
|
|
Packit |
67cb25 |
input matrix contains the Householder vectors which, together with the
|
|
Packit |
67cb25 |
Householder coefficients :data:`tau`, encode the orthogonal matrix
|
|
Packit |
67cb25 |
:math:`Q`. This storage scheme is the same as used by |lapack|. The
|
|
Packit |
67cb25 |
upper triangular part of :data:`A` is not referenced.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_symmtd_unpack (const gsl_matrix * A, const gsl_vector * tau, gsl_matrix * Q, gsl_vector * diag, gsl_vector * subdiag)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function unpacks the encoded symmetric tridiagonal decomposition
|
|
Packit |
67cb25 |
(:data:`A`, :data:`tau`) obtained from :func:`gsl_linalg_symmtd_decomp` into
|
|
Packit |
67cb25 |
the orthogonal matrix :data:`Q`, the vector of diagonal elements :data:`diag`
|
|
Packit |
67cb25 |
and the vector of subdiagonal elements :data:`subdiag`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_symmtd_unpack_T (const gsl_matrix * A, gsl_vector * diag, gsl_vector * subdiag)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function unpacks the diagonal and subdiagonal of the encoded
|
|
Packit |
67cb25 |
symmetric tridiagonal decomposition (:data:`A`, :data:`tau`) obtained from
|
|
Packit |
67cb25 |
:func:`gsl_linalg_symmtd_decomp` into the vectors :data:`diag` and :data:`subdiag`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: tridiagonal decomposition
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Tridiagonal Decomposition of Hermitian Matrices
|
|
Packit |
67cb25 |
===============================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A hermitian matrix :math:`A` can be factorized by similarity
|
|
Packit |
67cb25 |
transformations into the form,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: A = U T U^T
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`U` is a unitary matrix and :math:`T` is a real symmetric
|
|
Packit |
67cb25 |
tridiagonal matrix.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_hermtd_decomp (gsl_matrix_complex * A, gsl_vector_complex * tau)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function factorizes the hermitian matrix :data:`A` into the symmetric
|
|
Packit |
67cb25 |
tridiagonal decomposition :math:`U T U^T`. On output the real parts of
|
|
Packit |
67cb25 |
the diagonal and subdiagonal part of the input matrix :data:`A` contain
|
|
Packit |
67cb25 |
the tridiagonal matrix :math:`T`. The remaining lower triangular part of
|
|
Packit |
67cb25 |
the input matrix contains the Householder vectors which, together with
|
|
Packit |
67cb25 |
the Householder coefficients :data:`tau`, encode the unitary matrix
|
|
Packit |
67cb25 |
:math:`U`. This storage scheme is the same as used by |lapack|. The
|
|
Packit |
67cb25 |
upper triangular part of :data:`A` and imaginary parts of the diagonal are
|
|
Packit |
67cb25 |
not referenced.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: 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)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function unpacks the encoded tridiagonal decomposition (:data:`A`,
|
|
Packit |
67cb25 |
:data:`tau`) obtained from :func:`gsl_linalg_hermtd_decomp` into the
|
|
Packit |
67cb25 |
unitary matrix :data:`U`, the real vector of diagonal elements :data:`diag` and
|
|
Packit |
67cb25 |
the real vector of subdiagonal elements :data:`subdiag`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_hermtd_unpack_T (const gsl_matrix_complex * A, gsl_vector * diag, gsl_vector * subdiag)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function unpacks the diagonal and subdiagonal of the encoded
|
|
Packit |
67cb25 |
tridiagonal decomposition (:data:`A`, :data:`tau`) obtained from the
|
|
Packit |
67cb25 |
:func:`gsl_linalg_hermtd_decomp` into the real vectors
|
|
Packit |
67cb25 |
:data:`diag` and :data:`subdiag`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: Hessenberg decomposition
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Hessenberg Decomposition of Real Matrices
|
|
Packit |
67cb25 |
=========================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A general real matrix :math:`A` can be decomposed by orthogonal
|
|
Packit |
67cb25 |
similarity transformations into the form
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: A = U H U^T
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`U` is orthogonal and :math:`H` is an upper Hessenberg matrix,
|
|
Packit |
67cb25 |
meaning that it has zeros below the first subdiagonal. The
|
|
Packit |
67cb25 |
Hessenberg reduction is the first step in the Schur decomposition
|
|
Packit |
67cb25 |
for the nonsymmetric eigenvalue problem, but has applications in
|
|
Packit |
67cb25 |
other areas as well.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_hessenberg_decomp (gsl_matrix * A, gsl_vector * tau)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function computes the Hessenberg decomposition of the matrix
|
|
Packit |
67cb25 |
:data:`A` by applying the similarity transformation :math:`H = U^T A U`.
|
|
Packit |
67cb25 |
On output, :math:`H` is stored in the upper portion of :data:`A`. The
|
|
Packit |
67cb25 |
information required to construct the matrix :math:`U` is stored in
|
|
Packit |
67cb25 |
the lower triangular portion of :data:`A`. :math:`U` is a product
|
|
Packit |
67cb25 |
of :math:`N - 2` Householder matrices. The Householder vectors
|
|
Packit |
67cb25 |
are stored in the lower portion of :data:`A` (below the subdiagonal)
|
|
Packit |
67cb25 |
and the Householder coefficients are stored in the vector :data:`tau`.
|
|
Packit |
67cb25 |
:data:`tau` must be of length :data:`N`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_hessenberg_unpack (gsl_matrix * H, gsl_vector * tau, gsl_matrix * U)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function constructs the orthogonal matrix :math:`U` from the
|
|
Packit |
67cb25 |
information stored in the Hessenberg matrix :data:`H` along with the
|
|
Packit |
67cb25 |
vector :data:`tau`. :data:`H` and :data:`tau` are outputs from
|
|
Packit |
67cb25 |
:func:`gsl_linalg_hessenberg_decomp`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_hessenberg_unpack_accum (gsl_matrix * H, gsl_vector * tau, gsl_matrix * V)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function is similar to :func:`gsl_linalg_hessenberg_unpack`, except
|
|
Packit |
67cb25 |
it accumulates the matrix :data:`U` into :data:`V`, so that :math:`V' = VU`.
|
|
Packit |
67cb25 |
The matrix :data:`V` must be initialized prior to calling this function.
|
|
Packit |
67cb25 |
Setting :data:`V` to the identity matrix provides the same result as
|
|
Packit |
67cb25 |
:func:`gsl_linalg_hessenberg_unpack`. If :data:`H` is order :data:`N`, then
|
|
Packit |
67cb25 |
:data:`V` must have :data:`N` columns but may have any number of rows.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_hessenberg_set_zero (gsl_matrix * H)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function sets the lower triangular portion of :data:`H`, below
|
|
Packit |
67cb25 |
the subdiagonal, to zero. It is useful for clearing out the
|
|
Packit |
67cb25 |
Householder vectors after calling :func:`gsl_linalg_hessenberg_decomp`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: Hessenberg triangular decomposition
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Hessenberg-Triangular Decomposition of Real Matrices
|
|
Packit |
67cb25 |
====================================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A general real matrix pair (:math:`A`, :math:`B`) can be decomposed by
|
|
Packit |
67cb25 |
orthogonal similarity transformations into the form
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A &= U H V^T \\
|
|
Packit |
67cb25 |
B &= U R V^T
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A = U H V^T
|
|
Packit |
67cb25 |
B = U R V^T
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`U` and :math:`V` are orthogonal, :math:`H` is an upper
|
|
Packit |
67cb25 |
Hessenberg matrix, and :math:`R` is upper triangular. The
|
|
Packit |
67cb25 |
Hessenberg-Triangular reduction is the first step in the generalized
|
|
Packit |
67cb25 |
Schur decomposition for the generalized eigenvalue problem.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_hesstri_decomp (gsl_matrix * A, gsl_matrix * B, gsl_matrix * U, gsl_matrix * V, gsl_vector * work)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function computes the Hessenberg-Triangular decomposition of the
|
|
Packit |
67cb25 |
matrix pair (:data:`A`, :data:`B`). On output, :math:`H` is stored in :data:`A`,
|
|
Packit |
67cb25 |
and :math:`R` is stored in :data:`B`. If :data:`U` and :data:`V` are provided
|
|
Packit |
67cb25 |
(they may be null), the similarity transformations are stored in them.
|
|
Packit |
67cb25 |
Additional workspace of length :math:`N` is needed in :data:`work`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: bidiagonalization of real matrices
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Bidiagonalization
|
|
Packit |
67cb25 |
=================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A general matrix :math:`A` can be factorized by similarity
|
|
Packit |
67cb25 |
transformations into the form,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: A = U B V^T
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`U` and :math:`V` are orthogonal matrices and :math:`B` is a
|
|
Packit |
67cb25 |
:math:`N`-by-:math:`N` bidiagonal matrix with non-zero entries only on the
|
|
Packit |
67cb25 |
diagonal and superdiagonal. The size of :data:`U` is :math:`M`-by-:math:`N`
|
|
Packit |
67cb25 |
and the size of :data:`V` is :math:`N`-by-:math:`N`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_bidiag_decomp (gsl_matrix * A, gsl_vector * tau_U, gsl_vector * tau_V)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function factorizes the :math:`M`-by-:math:`N` matrix :data:`A` into
|
|
Packit |
67cb25 |
bidiagonal form :math:`U B V^T`. The diagonal and superdiagonal of the
|
|
Packit |
67cb25 |
matrix :math:`B` are stored in the diagonal and superdiagonal of :data:`A`.
|
|
Packit |
67cb25 |
The orthogonal matrices :math:`U` and :data:`V` are stored as compressed
|
|
Packit |
67cb25 |
Householder vectors in the remaining elements of :data:`A`. The
|
|
Packit |
67cb25 |
Householder coefficients are stored in the vectors :data:`tau_U` and
|
|
Packit |
67cb25 |
:data:`tau_V`. The length of :data:`tau_U` must equal the number of
|
|
Packit |
67cb25 |
elements in the diagonal of :data:`A` and the length of :data:`tau_V` should
|
|
Packit |
67cb25 |
be one element shorter.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: 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)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function unpacks the bidiagonal decomposition of :data:`A` produced by
|
|
Packit |
67cb25 |
:func:`gsl_linalg_bidiag_decomp`, (:data:`A`, :data:`tau_U`, :data:`tau_V`)
|
|
Packit |
67cb25 |
into the separate orthogonal matrices :data:`U`, :data:`V` and the diagonal
|
|
Packit |
67cb25 |
vector :data:`diag` and superdiagonal :data:`superdiag`. Note that :data:`U`
|
|
Packit |
67cb25 |
is stored as a compact :math:`M`-by-:math:`N` orthogonal matrix satisfying
|
|
Packit |
67cb25 |
:math:`U^T U = I` for efficiency.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_bidiag_unpack2 (gsl_matrix * A, gsl_vector * tau_U, gsl_vector * tau_V, gsl_matrix * V)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function unpacks the bidiagonal decomposition of :data:`A` produced by
|
|
Packit |
67cb25 |
:func:`gsl_linalg_bidiag_decomp`, (:data:`A`, :data:`tau_U`, :data:`tau_V`)
|
|
Packit |
67cb25 |
into the separate orthogonal matrices :data:`U`, :data:`V` and the diagonal
|
|
Packit |
67cb25 |
vector :data:`diag` and superdiagonal :data:`superdiag`. The matrix :data:`U`
|
|
Packit |
67cb25 |
is stored in-place in :data:`A`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_bidiag_unpack_B (const gsl_matrix * A, gsl_vector * diag, gsl_vector * superdiag)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function unpacks the diagonal and superdiagonal of the bidiagonal
|
|
Packit |
67cb25 |
decomposition of :data:`A` from :func:`gsl_linalg_bidiag_decomp`, into
|
|
Packit |
67cb25 |
the diagonal vector :data:`diag` and superdiagonal vector :data:`superdiag`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: Givens rotation
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Givens Rotations
|
|
Packit |
67cb25 |
================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A Givens rotation is a rotation in the plane acting on two elements
|
|
Packit |
67cb25 |
of a given vector. It can be represented in matrix form as
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
G(i,j,\theta) =
|
|
Packit |
67cb25 |
\left(
|
|
Packit |
67cb25 |
\begin{matrix}
|
|
Packit |
67cb25 |
1 & \ldots & 0 & \ldots & 0 & \ldots & 0 \\
|
|
Packit |
67cb25 |
\vdots & \ddots & \vdots & & \vdots & & \vdots \\
|
|
Packit |
67cb25 |
0 & \ldots & \cos{\theta} & \ldots & -\sin{\theta} & \ldots & 0 \\
|
|
Packit |
67cb25 |
\vdots & & \vdots & \ddots & \vdots & & \vdots \\
|
|
Packit |
67cb25 |
0 & \ldots & \sin{\theta} & \ldots & \cos{\theta} & \ldots & 0 \\
|
|
Packit |
67cb25 |
\vdots & & \vdots & & \vdots & \ddots & \vdots \\
|
|
Packit |
67cb25 |
0 & \ldots & 0 & \ldots & 0 & \ldots & 1
|
|
Packit |
67cb25 |
\end{matrix}
|
|
Packit |
67cb25 |
\right)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where the :math:`\cos{\theta}` and :math:`\sin{\theta}` appear at
|
|
Packit |
67cb25 |
the intersection of the :math:`i`-th and :math:`j`-th rows and columns.
|
|
Packit |
67cb25 |
When acting on a vector :math:`x`, :math:`G(i,j,\theta) x` performs
|
|
Packit |
67cb25 |
a rotation of the :math:`(i,j)` elements of :math:`x`. Givens
|
|
Packit |
67cb25 |
rotations are typically used to introduce zeros in
|
|
Packit |
67cb25 |
vectors, such as during the QR decomposition of a matrix. In this
|
|
Packit |
67cb25 |
case, it is typically desired to find :math:`c` and :math:`s` such that
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
\left(
|
|
Packit |
67cb25 |
\begin{matrix}
|
|
Packit |
67cb25 |
c & -s \\
|
|
Packit |
67cb25 |
s & c
|
|
Packit |
67cb25 |
\end{matrix}
|
|
Packit |
67cb25 |
\right)
|
|
Packit |
67cb25 |
\left(
|
|
Packit |
67cb25 |
\begin{matrix}
|
|
Packit |
67cb25 |
a \\
|
|
Packit |
67cb25 |
b
|
|
Packit |
67cb25 |
\end{matrix}
|
|
Packit |
67cb25 |
\right) =
|
|
Packit |
67cb25 |
\left(
|
|
Packit |
67cb25 |
\begin{matrix}
|
|
Packit |
67cb25 |
r \\
|
|
Packit |
67cb25 |
0
|
|
Packit |
67cb25 |
\end{matrix}
|
|
Packit |
67cb25 |
\right)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
[ c -s ] [ a ] = [ r ]
|
|
Packit |
67cb25 |
[ s c ] [ b ] [ 0 ]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
with :math:`r = \sqrt{a^2 + b^2}`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: void gsl_linalg_givens (const double a, const double b, double * c, double * s)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function computes :math:`c = \cos{\theta}` and :math:`s = \sin{\theta}`
|
|
Packit |
67cb25 |
so that the Givens matrix :math:`G(\theta)` acting on the
|
|
Packit |
67cb25 |
vector :math:`(a,b)` produces :math:`(r, 0)`, with :math:`r = \sqrt{a^2 + b^2}`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: void gsl_linalg_givens_gv (gsl_vector * v, const size_t i, const size_t j, const double c, const double s)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function applies the Givens rotation defined by
|
|
Packit |
67cb25 |
:math:`c = \cos{\theta}` and :math:`s = \sin{\theta}` to the :data:`i`
|
|
Packit |
67cb25 |
and :data:`j` elements of :data:`v`. On output,
|
|
Packit |
67cb25 |
:math:`(v(i),v(j)) \leftarrow G(\theta) (v(i),v(j))`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: Householder matrix
|
|
Packit |
67cb25 |
single: Householder transformation
|
|
Packit |
67cb25 |
single: transformation, Householder
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Householder Transformations
|
|
Packit |
67cb25 |
===========================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A Householder transformation is a rank-1 modification of the identity
|
|
Packit |
67cb25 |
matrix which can be used to zero out selected elements of a vector. A
|
|
Packit |
67cb25 |
Householder matrix :math:`P` takes the form,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: P = I - \tau v v^T
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`v` is a vector (called the *Householder vector*) and
|
|
Packit |
67cb25 |
:math:`\tau = 2/(v^T v)`. The functions described in this section use the
|
|
Packit |
67cb25 |
rank-1 structure of the Householder matrix to create and apply
|
|
Packit |
67cb25 |
Householder transformations efficiently.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: double gsl_linalg_householder_transform (gsl_vector * w)
|
|
Packit |
67cb25 |
gsl_complex gsl_linalg_complex_householder_transform (gsl_vector_complex * w)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function prepares a Householder transformation :math:`P = I - \tau v v^T`
|
|
Packit |
67cb25 |
which can be used to zero all the elements of the input vector :data:`w`
|
|
Packit |
67cb25 |
except the first. On output the Householder vector :data:`v` is stored in
|
|
Packit |
67cb25 |
:data:`w` and the scalar :math:`\tau` is returned. The householder vector
|
|
Packit |
67cb25 |
:data:`v` is normalized so that :code:`v[0] = 1`, however this 1 is not
|
|
Packit |
67cb25 |
stored in the output vector. Instead, :code:`w[0]` is set to
|
|
Packit |
67cb25 |
the first element of the transformed vector, so that if
|
|
Packit |
67cb25 |
:math:`u = P w`, :code:`w[0] = u[0]` on output and the remainder
|
|
Packit |
67cb25 |
of :math:`u` is zero.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_householder_hm (double tau, const gsl_vector * v, gsl_matrix * A)
|
|
Packit |
67cb25 |
int gsl_linalg_complex_householder_hm (gsl_complex tau, const gsl_vector_complex * v, gsl_matrix_complex * A)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function applies the Householder matrix :math:`P` defined by the
|
|
Packit |
67cb25 |
scalar :data:`tau` and the vector :data:`v` to the left-hand side of the
|
|
Packit |
67cb25 |
matrix :data:`A`. On output the result :math:`P A` is stored in :data:`A`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_householder_mh (double tau, const gsl_vector * v, gsl_matrix * A)
|
|
Packit |
67cb25 |
int gsl_linalg_complex_householder_mh (gsl_complex tau, const gsl_vector_complex * v, gsl_matrix_complex * A)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function applies the Householder matrix :math:`P` defined by the
|
|
Packit |
67cb25 |
scalar :data:`tau` and the vector :data:`v` to the right-hand side of the
|
|
Packit |
67cb25 |
matrix :data:`A`. On output the result :math:`A P` is stored in :data:`A`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_householder_hv (double tau, const gsl_vector * v, gsl_vector * w)
|
|
Packit |
67cb25 |
int gsl_linalg_complex_householder_hv (gsl_complex tau, const gsl_vector_complex * v, gsl_vector_complex * w)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function applies the Householder transformation :math:`P` defined by
|
|
Packit |
67cb25 |
the scalar :data:`tau` and the vector :data:`v` to the vector :data:`w`. On
|
|
Packit |
67cb25 |
output the result :math:`P w` is stored in :data:`w`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. @deftypefun int gsl_linalg_householder_hm1 (double tau, gsl_matrix * A)
|
|
Packit |
67cb25 |
.. This function applies the Householder transform, defined by the scalar
|
|
Packit |
67cb25 |
.. :data:`tau` and the vector :data:`v`, to a matrix being build up from the
|
|
Packit |
67cb25 |
.. identity matrix, using the first column of :data:`A` as a householder vector.
|
|
Packit |
67cb25 |
.. @end deftypefun
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: solution of linear system by Householder transformations
|
|
Packit |
67cb25 |
single: Householder linear solver
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Householder solver for linear systems
|
|
Packit |
67cb25 |
=====================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_HH_solve (gsl_matrix * A, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the system :math:`A x = b` directly using
|
|
Packit |
67cb25 |
Householder transformations. On output the solution is stored in :data:`x`
|
|
Packit |
67cb25 |
and :data:`b` is not modified. The matrix :data:`A` is destroyed by the
|
|
Packit |
67cb25 |
Householder transformations.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_HH_svx (gsl_matrix * A, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the system :math:`A x = b` in-place using
|
|
Packit |
67cb25 |
Householder transformations. On input :data:`x` should contain the
|
|
Packit |
67cb25 |
right-hand side :math:`b`, which is replaced by the solution on output. The
|
|
Packit |
67cb25 |
matrix :data:`A` is destroyed by the Householder transformations.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: tridiagonal systems
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Tridiagonal Systems
|
|
Packit |
67cb25 |
===================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The functions described in this section efficiently solve symmetric,
|
|
Packit |
67cb25 |
non-symmetric and cyclic tridiagonal systems with minimal storage.
|
|
Packit |
67cb25 |
Note that the current implementations of these functions use a variant
|
|
Packit |
67cb25 |
of Cholesky decomposition, so the tridiagonal matrix must be positive
|
|
Packit |
67cb25 |
definite. For non-positive definite matrices, the functions return
|
|
Packit |
67cb25 |
the error code :macro:`GSL_ESING`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_solve_tridiag (const gsl_vector * diag, const gsl_vector * e, const gsl_vector * f, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the general :math:`N`-by-:math:`N` system :math:`A x = b`
|
|
Packit |
67cb25 |
where :data:`A` is tridiagonal (:math:`N \geq 2`).
|
|
Packit |
67cb25 |
The super-diagonal and
|
|
Packit |
67cb25 |
sub-diagonal vectors :data:`e` and :data:`f` must be one element shorter
|
|
Packit |
67cb25 |
than the diagonal vector :data:`diag`. The form of :data:`A` for the 4-by-4
|
|
Packit |
67cb25 |
case is shown below,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A =
|
|
Packit |
67cb25 |
\left(
|
|
Packit |
67cb25 |
\begin{matrix}
|
|
Packit |
67cb25 |
d_0&e_0& 0& 0\\
|
|
Packit |
67cb25 |
f_0&d_1&e_1& 0\\
|
|
Packit |
67cb25 |
0 &f_1&d_2&e_2\\
|
|
Packit |
67cb25 |
0 &0 &f_2&d_3
|
|
Packit |
67cb25 |
\end{matrix}
|
|
Packit |
67cb25 |
\right)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A = ( d_0 e_0 0 0 )
|
|
Packit |
67cb25 |
( f_0 d_1 e_1 0 )
|
|
Packit |
67cb25 |
( 0 f_1 d_2 e_2 )
|
|
Packit |
67cb25 |
( 0 0 f_2 d_3 )
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_solve_symm_tridiag (const gsl_vector * diag, const gsl_vector * e, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the general :math:`N`-by-:math:`N` system :math:`A x = b`
|
|
Packit |
67cb25 |
where :data:`A` is symmetric tridiagonal (:math:`N \geq 2`).
|
|
Packit |
67cb25 |
The off-diagonal vector
|
|
Packit |
67cb25 |
:data:`e` must be one element shorter than the diagonal vector :data:`diag`.
|
|
Packit |
67cb25 |
The form of :data:`A` for the 4-by-4 case is shown below,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A =
|
|
Packit |
67cb25 |
\left(
|
|
Packit |
67cb25 |
\begin{matrix}
|
|
Packit |
67cb25 |
d_0&e_0& 0& 0\\
|
|
Packit |
67cb25 |
e_0&d_1&e_1& 0\\
|
|
Packit |
67cb25 |
0 &e_1&d_2&e_2\\
|
|
Packit |
67cb25 |
0 &0 &e_2&d_3
|
|
Packit |
67cb25 |
\end{matrix}
|
|
Packit |
67cb25 |
\right)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A = ( d_0 e_0 0 0 )
|
|
Packit |
67cb25 |
( e_0 d_1 e_1 0 )
|
|
Packit |
67cb25 |
( 0 e_1 d_2 e_2 )
|
|
Packit |
67cb25 |
( 0 0 e_2 d_3 )
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: 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)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the general :math:`N`-by-:math:`N` system :math:`A x = b`
|
|
Packit |
67cb25 |
where :data:`A` is cyclic tridiagonal (:math:`N \geq 3`).
|
|
Packit |
67cb25 |
The cyclic super-diagonal and
|
|
Packit |
67cb25 |
sub-diagonal vectors :data:`e` and :data:`f` must have the same number of
|
|
Packit |
67cb25 |
elements as the diagonal vector :data:`diag`. The form of :data:`A` for the
|
|
Packit |
67cb25 |
4-by-4 case is shown below,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A =
|
|
Packit |
67cb25 |
\left(
|
|
Packit |
67cb25 |
\begin{matrix}
|
|
Packit |
67cb25 |
d_0&e_0& 0 &f_3\\
|
|
Packit |
67cb25 |
f_0&d_1&e_1& 0 \\
|
|
Packit |
67cb25 |
0 &f_1&d_2&e_2\\
|
|
Packit |
67cb25 |
e_3& 0 &f_2&d_3
|
|
Packit |
67cb25 |
\end{matrix}
|
|
Packit |
67cb25 |
\right)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A = ( d_0 e_0 0 f_3 )
|
|
Packit |
67cb25 |
( f_0 d_1 e_1 0 )
|
|
Packit |
67cb25 |
( 0 f_1 d_2 e_2 )
|
|
Packit |
67cb25 |
( e_3 0 f_2 d_3 )
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_solve_symm_cyc_tridiag (const gsl_vector * diag, const gsl_vector * e, const gsl_vector * b, gsl_vector * x)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function solves the general :math:`N`-by-:math:`N` system :math:`A x = b`
|
|
Packit |
67cb25 |
where :data:`A` is symmetric cyclic tridiagonal (:math:`N \geq 3`).
|
|
Packit |
67cb25 |
The cyclic
|
|
Packit |
67cb25 |
off-diagonal vector :data:`e` must have the same number of elements as the
|
|
Packit |
67cb25 |
diagonal vector :data:`diag`. The form of :data:`A` for the 4-by-4 case is
|
|
Packit |
67cb25 |
shown below,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A =
|
|
Packit |
67cb25 |
\left(
|
|
Packit |
67cb25 |
\begin{matrix}
|
|
Packit |
67cb25 |
d_0&e_0& 0 &e_3\\
|
|
Packit |
67cb25 |
e_0&d_1&e_1& 0 \\
|
|
Packit |
67cb25 |
0 &e_1&d_2&e_2\\
|
|
Packit |
67cb25 |
e_3& 0 &e_2&d_3
|
|
Packit |
67cb25 |
\end{matrix}
|
|
Packit |
67cb25 |
\right)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A = ( d_0 e_0 0 e_3 )
|
|
Packit |
67cb25 |
( e_0 d_1 e_1 0 )
|
|
Packit |
67cb25 |
( 0 e_1 d_2 e_2 )
|
|
Packit |
67cb25 |
( e_3 0 e_2 d_3 )
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: triangular systems
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Triangular Systems
|
|
Packit |
67cb25 |
==================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_tri_upper_invert (gsl_matrix * T)
|
|
Packit |
67cb25 |
int gsl_linalg_tri_lower_invert (gsl_matrix * T)
|
|
Packit |
67cb25 |
int gsl_linalg_tri_upper_unit_invert (gsl_matrix * T)
|
|
Packit |
67cb25 |
int gsl_linalg_tri_lower_unit_invert (gsl_matrix * T)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions calculate the in-place inverse of the triangular matrix :data:`T`. When
|
|
Packit |
67cb25 |
the :code:`upper` prefix is specified, then the upper triangle of :data:`T` is used, and when
|
|
Packit |
67cb25 |
the :code:`lower` prefix is specified, the lower triangle is used. If the :code:`unit`
|
|
Packit |
67cb25 |
prefix is specified, then the diagonal elements of the matrix :data:`T` are taken as
|
|
Packit |
67cb25 |
unity and are not referenced. Otherwise the diagonal elements are used in the inversion.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_tri_upper_rcond (const gsl_matrix * T, double * rcond, gsl_vector * work)
|
|
Packit |
67cb25 |
int gsl_linalg_tri_lower_rcond (const gsl_matrix * T, double * rcond, gsl_vector * work)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions estimate the reciprocal condition number, in the 1-norm, of the upper or lower
|
|
Packit |
67cb25 |
:math:`N`-by-:math:`N` triangular matrix :data:`T`. The reciprocal condition number
|
|
Packit |
67cb25 |
is stored in :data:`rcond` on output, and is defined by :math:`1 / (||T||_1 \cdot ||T^{-1}||_1)`.
|
|
Packit |
67cb25 |
Additional workspace of size :math:`3 N` is required in :data:`work`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: balancing matrices
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. _balancing:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Balancing
|
|
Packit |
67cb25 |
=========
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The process of balancing a matrix applies similarity transformations
|
|
Packit |
67cb25 |
to make the rows and columns have comparable norms. This is
|
|
Packit |
67cb25 |
useful, for example, to reduce roundoff errors in the solution
|
|
Packit |
67cb25 |
of eigenvalue problems. Balancing a matrix :math:`A` consists
|
|
Packit |
67cb25 |
of replacing :math:`A` with a similar matrix
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: A' = D^{-1} A D
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where :math:`D` is a diagonal matrix whose entries are powers
|
|
Packit |
67cb25 |
of the floating point radix.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_linalg_balance_matrix (gsl_matrix * A, gsl_vector * D)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function replaces the matrix :data:`A` with its balanced counterpart
|
|
Packit |
67cb25 |
and stores the diagonal elements of the similarity transformation
|
|
Packit |
67cb25 |
into the vector :data:`D`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Examples
|
|
Packit |
67cb25 |
========
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The following program solves the linear system :math:`A x = b`. The
|
|
Packit |
67cb25 |
system to be solved is,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
\left(
|
|
Packit |
67cb25 |
\begin{matrix}
|
|
Packit |
67cb25 |
0.18& 0.60& 0.57& 0.96\\
|
|
Packit |
67cb25 |
0.41& 0.24& 0.99& 0.58\\
|
|
Packit |
67cb25 |
0.14& 0.30& 0.97& 0.66\\
|
|
Packit |
67cb25 |
0.51& 0.13& 0.19& 0.85
|
|
Packit |
67cb25 |
\end{matrix}
|
|
Packit |
67cb25 |
\right)
|
|
Packit |
67cb25 |
\left(
|
|
Packit |
67cb25 |
\begin{matrix}
|
|
Packit |
67cb25 |
x_0\\
|
|
Packit |
67cb25 |
x_1\\
|
|
Packit |
67cb25 |
x_2\\
|
|
Packit |
67cb25 |
x_3
|
|
Packit |
67cb25 |
\end{matrix}
|
|
Packit |
67cb25 |
\right)
|
|
Packit |
67cb25 |
=
|
|
Packit |
67cb25 |
\left(
|
|
Packit |
67cb25 |
\begin{matrix}
|
|
Packit |
67cb25 |
1.0\\
|
|
Packit |
67cb25 |
2.0\\
|
|
Packit |
67cb25 |
3.0\\
|
|
Packit |
67cb25 |
4.0
|
|
Packit |
67cb25 |
\end{matrix}
|
|
Packit |
67cb25 |
\right)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
[ 0.18 0.60 0.57 0.96 ] [x0] [1.0]
|
|
Packit |
67cb25 |
[ 0.41 0.24 0.99 0.58 ] [x1] = [2.0]
|
|
Packit |
67cb25 |
[ 0.14 0.30 0.97 0.66 ] [x2] [3.0]
|
|
Packit |
67cb25 |
[ 0.51 0.13 0.19 0.85 ] [x3] [4.0]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
and the solution is found using LU decomposition of the matrix :math:`A`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: examples/linalglu.c
|
|
Packit |
67cb25 |
:code:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Here is the output from the program,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: examples/linalglu.txt
|
|
Packit |
67cb25 |
:code:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This can be verified by multiplying the solution :math:`x` by the
|
|
Packit |
67cb25 |
original matrix :math:`A` using |octave|,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
octave> A = [ 0.18, 0.60, 0.57, 0.96;
|
|
Packit |
67cb25 |
0.41, 0.24, 0.99, 0.58;
|
|
Packit |
67cb25 |
0.14, 0.30, 0.97, 0.66;
|
|
Packit |
67cb25 |
0.51, 0.13, 0.19, 0.85 ];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
octave> x = [ -4.05205; -12.6056; 1.66091; 8.69377];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
octave> A * x
|
|
Packit |
67cb25 |
ans =
|
|
Packit |
67cb25 |
1.0000
|
|
Packit |
67cb25 |
2.0000
|
|
Packit |
67cb25 |
3.0000
|
|
Packit |
67cb25 |
4.0000
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This reproduces the original right-hand side vector, :math:`b`, in
|
|
Packit |
67cb25 |
accordance with the equation :math:`A x = b`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
References and Further Reading
|
|
Packit |
67cb25 |
==============================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Further information on the algorithms described in this section can be
|
|
Packit |
67cb25 |
found in the following book,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* G. H. Golub, C. F. Van Loan, "Matrix Computations" (3rd Ed, 1996),
|
|
Packit |
67cb25 |
Johns Hopkins University Press, ISBN 0-8018-5414-8.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The |lapack| library is described in the following manual,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* *LAPACK Users' Guide* (Third Edition, 1999), Published by SIAM,
|
|
Packit |
67cb25 |
ISBN 0-89871-447-8
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The |lapack| source code can be found at http://www.netlib.org/lapack,
|
|
Packit |
67cb25 |
along with an online copy of the users guide.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The Modified Golub-Reinsch algorithm is described in the following paper,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* T.F. Chan, "An Improved Algorithm for Computing the Singular Value
|
|
Packit |
67cb25 |
Decomposition", ACM Transactions on Mathematical Software, 8
|
|
Packit |
67cb25 |
(1982), pp 72--83.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The Jacobi algorithm for singular value decomposition is described in
|
|
Packit |
67cb25 |
the following papers,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* J.C. Nash, "A one-sided transformation method for the singular value
|
|
Packit |
67cb25 |
decomposition and algebraic eigenproblem", Computer Journal,
|
|
Packit |
67cb25 |
Volume 18, Number 1 (1975), p 74--76
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* J.C. Nash and S. Shlien "Simple algorithms for the partial singular
|
|
Packit |
67cb25 |
value decomposition", Computer Journal, Volume 30 (1987), p
|
|
Packit |
67cb25 |
268--275.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* J. Demmel, K. Veselic, "Jacobi's Method is more accurate than
|
|
Packit |
67cb25 |
QR", Lapack Working Note 15 (LAWN-15), October 1989. Available
|
|
Packit |
67cb25 |
from netlib, http://www.netlib.org/lapack/ in the :code:`lawns` or
|
|
Packit |
67cb25 |
:code:`lawnspdf` directories.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The algorithm for estimating a matrix condition number is described in
|
|
Packit |
67cb25 |
the following paper,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* N. J. Higham, "FORTRAN codes for estimating the one-norm of
|
|
Packit |
67cb25 |
a real or complex matrix, with applications to condition estimation",
|
|
Packit |
67cb25 |
ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
|