# Author : Fabian Jakobs
"""
Functions for solving linear systems.
The library provides simple linear algebra operations which operate
directly on the gsl_vector and gsl_matrix objects. These are intended
for use with "small" systems where simple algorithms are acceptable.
Anyone interested in large systems will want to use the sophisticated
routines found in LAPACK. The Fortran version of LAPACK is recommended
as the standard package for linear algebra. It supports blocked algorithms,
specialized data representations and other optimizations.
"""
from __future__ import print_function
import pygsl
from . import _gslwrap
from .permutation import Permutation
import pygsl._numobj as numx
zeros = numx.zeros
array = numx.array
get_typecode = pygsl.get_typecode
array_typed_copy = pygsl.array_typed_copy
Float = pygsl.Float
Complex = pygsl.Complex
#
# LU Decomposition
#
[docs]
def LU_decomp(A):
"""
returns (LU, P, signum)
This function factorizes the square matrix A into the LU decomposition
PA = LU. On output the diagonal and upper triangular part of the return
matrix contain the matrix U. The lower triangular part of the input matrix
(excluding the diagonal) contains L. The diagonal elements of L are unity,
and are not stored.
The permutation matrix P is encoded in the permutation p. The j-th column
of the matrix P is given by the k-th column of the identity matrix, where
k = p_j the j-th element of the permutation vector. The sign of the
permutation is given by signum. It has the value (-1)^n, where n is the
number of interchanges in the permutation.
The algorithm used in the decomposition is Gaussian Elimination with
partial pivoting (Golub & Van Loan, Matrix Computations, Algorithm 3.4.1).
Todo:
Check if creating a copy of array A is required
"""
p = Permutation(A.shape[1])
code = get_typecode(A)
tmp = array_typed_copy(A)
if code == Complex:
# Now all error flags are turned into python exceptions. So no
# unpack necessary any longer.
An, signum = _gslwrap.pygsl_linalg_complex_LU_decomp(tmp, p)
elif code == Float:
An, signum = _gslwrap.pygsl_linalg_LU_decomp(tmp, p)
else:
#print code, Float, Complex
raise TypeError("LU must be of type Float or Complex")
return (An, p, signum)
[docs]
def LU_unpack(LU):
"""
returns (L,U)
This function splits the matrix LU into the the upper matrix U and
the lower matrix L. The diagonal of L is the identity.
"""
code = get_typecode(LU)
u = zeros(LU.shape, code)
l = numx.identity(LU.shape[0], code)
for i in range(LU.shape[0]):
u[i, i: ] = LU[i, i:]
l[i, 0:i] = LU[i, :i]
return (l, u)
[docs]
def LU_solve(LU, p, b):
"""
This function solves the system A x = b using the LU decomposition of A
into (LU, p) given by LU_decomp.
"""
code = get_typecode(LU)
x = zeros(LU.shape[1], code)
if code == Complex:
_gslwrap.gsl_linalg_complex_LU_solve(LU, p, b, x)
elif code == Float:
_gslwrap.gsl_linalg_LU_solve(LU, p, b, x)
else:
raise TypeError( "LU must be of type Float or Complex")
return x
[docs]
def LU_refine(A, LU, p, b, x):
"""
returns (x, residual)
This functions applies an iterative improvement to x, the solution
of A x = b, using the LU decomposition of A into (LU,p). The initial
residual r = A x - b is also computed and stored in residual.
"""
code = get_typecode(LU)
raise NotImplementedError("This function is not (yet implemented)")
# if code == Numeric.Complex:
# _gslwrap.gsl_linalg_complex_LU_refine(A, LU, p, b, x, residual)
# elif code == Numeric.Float:
# _gslwrap.gsl_linalg_LU_refine
# else:
# raise TypeError, "LU must be of type Float or Complex"
[docs]
def LU_invert(LU, p):
"""
returns inverse
This function computes the inverse of a matrix A from its LU
decomposition (LU,p), storing the result in the matrix inverse.
The inverse is computed by solving the system A x = b for each column of
the identity matrix. It is preferable to avoid direct computation of the
inverse whenever possible.
"""
code = get_typecode(LU)
inverse = zeros(LU.shape, code)
if code == Complex:
_gslwrap.gls_linalg_complex_LU_invert(LU, p, inverse)
elif code == Float:
_gslwrap.gsl_linalg_LU_invert(LU, p, inverse)
else:
raise TypeError("LU must be of type Float or Complex")
return inverse
[docs]
def LU_det(LU, signum):
"""
returns determinant
This function computes the determinant of a matrix A from its LU
decomposition, LU. The determinant is computed as the product of
the diagonal elements of U and the sign of the row permutation signum.
"""
code = get_typecode(LU)
if code == Complex:
return _gslwrap.gls_linalg_complex_LU_det(LU, signum)
elif code == Float:
return _gslwrap.gsl_linalg_LU_det(LU, signum)
else:
raise TypeError("LU must be of type Float or Complex")
[docs]
def LU_lndet(LU):
r"""
This function computes the logarithm of the absolute value of the
determinant of a matrix A, \ln|det(A)|, from its LU decomposition, LU.
This function may be useful if the direct computation of the determinant
would overflow or underflow.
"""
code = get_typecode(LU)
if code == Complex:
return _gslwrap.gls_linalg_complex_LU_lndet(LU)
elif code == Float:
return _gslwrap.gsl_linalg_LU_lndet(LU)
else:
raise TypeError("LU must be of type Float or Complex")
[docs]
def LU_sgndet(LU, signum):
"""
This function computes the sign or phase factor of the determinant of a
matrix A, det(A)/|det(A)|, from its LU decomposition, LU.
"""
code = get_typecode(LU)
if code == Complex:
return _gslwrap.gls_linalg_complex_LU_sgndet(LU, signum)
elif code == Float:
return _gslwrap.gsl_linalg_LU_sgndet(LU, signum)
else:
raise TypeError("LU must be of type Float or Complex")
#
# QR Decomposition
#
[docs]
def QR_decomp(A):
r"""
returns (QR, tau)
Function: int gsl_linalg_QR_decomp (gsl_matrix * A, gsl_vector * tau)
This function factorizes the M-by-N matrix A into the QR
decomposition A = Q R. On output the diagonal and upper triangular
part of the input matrix contain the matrix R. The vector tau and the
columns of the lower triangular part of the matrix QR contain the
Householder coefficients and Householder vectors which encode the
orthogonal matrix Q. The vector tau must be of length k=\min(M,N).
The matrix Q is related to these components by, Q = Q_k ... Q_2 Q_1
where Q_i = I - \tau_i v_i v_i^T and v_i is the Householder vector
v_i = (0,...,1,A(i+1,i),A(i+2,i),...,A(m,i)).
This is the same storage scheme as used by LAPACK.
The algorithm used to perform the decomposition is Householder QR
(Golub & Van Loan, Matrix Computations, Algorithm 5.2.1).
"""
code = get_typecode(A)
qr = array_typed_copy(A, code)
tau = zeros((min(A.shape),), code)
_gslwrap.gsl_linalg_QR_decomp(qr, tau)
return (qr, tau)
[docs]
def QR_solve(QR, tau, b):
"""
returns x
This function solves the system A x = b using the QR decomposition of
A into (QR, tau) given by gsl_linalg_QR_decomp.
"""
code = get_typecode(b)
x = zeros(QR.shape[1], code)
_gslwrap.gsl_linalg_QR_solve(QR, tau, b, x)
return x
[docs]
def QR_lssolve(QR, tau, b):
"""
returns (x, residual)
This function finds the least squares solution to the overdetermined
system A x = b where the matrix A has more rows than columns. The least
squares solution minimizes the Euclidean norm of the residual, ||Ax - b||.
The routine uses the QR decomposition of A into (QR, tau) given by
gsl_linalg_QR_decomp. The solution is returned in x. The residual is
computed as a by-product and stored in residual.
"""
code = get_typecode(QR)
x = zeros(QR.shape[1], code)
residual = zeros(QR.shape[0], code)
_gslwrap.gsl_linalg_QR_lssolve(QR, tau, b, x, residual)
return (x, residual)
[docs]
def QR_QTvec(QR, tau, v):
"""
returns v'
This function applies the matrix Q^T encoded in the decomposition
(QR,tau) to the vector v, storing the result Q^T v in v'. The matrix
multiplication is carried out directly using the encoding of the
Householder vectors without needing to form the full matrix Q^T.
"""
vn = array_typed_copy(v)
_gslwrap.gsl_linalg_QR_QTvec(QR,tau,vn)
return vn
[docs]
def QR_Qvec(QR, tau, v):
"""
returns v'
This function applies the matrix Q encoded in the decomposition
(QR,tau) to the vector v, storing the result Q v in v'. The matrix
multiplication is carried out directly using the encoding of the
Householder vectors without needing to form the full matrix Q.
"""
vn = array_typed_copy(v)
_gslwrap.gsl_linalg_QR_Qvec(QR,tau,vn)
return vn
[docs]
def QR_Rsolve(QR, b):
"""
returns x
This function solves the triangular system R x = b for x.
It may be useful if the product b' = Q^T b has already been computed
using gsl.linalg.QR_QTvec.
"""
x = zeros(QR.shape[1], get_typecode(b))
_gslwrap.gsl_linalg_QR_Rsolve(QR, b, x)
return x
[docs]
def QR_unpack(QR, tau):
"""
returns (Q, R)
This function unpacks the encoded QR decomposition (QR,tau) into the
matrices Q and R, where Q is M-by-M and R is M-by-N.
"""
(m, n) = QR.shape
code = get_typecode(QR)
q = zeros([m,m], code)
r = zeros([m,n], code)
_gslwrap.gsl_linalg_QR_unpack(QR, tau, q, r)
return (q,r)
[docs]
def QR_QRsolve(Q, R, b):
"""
returns x
This function solves the system R x = Q^T b for x. It can be used when
the QR decomposition of a matrix is available in unpacked form as (Q,R).
"""
x = zeros(R.shape[1], get_typecode(R))
_gslwrap.gsl_linalg_QR_QRsolve(Q, R, b, x)
return x
[docs]
def QR_update(Q, R, w, v):
"""
This function performs a rank-1 update w v^T of the QR
decomposition (Q, R). The update is given by Q'R' = Q R + w v^T
where the output matrices Q' and R' are also orthogonal and right
triangular. Note that Q and R are overwritten with Q' and R'!
"""
raise NotImplementedError("Please verify the output of this function!")
wn = array_typed_copy(w)
_gslwrap.gsl_linalg_QR_update(Q, R, wn, v)
[docs]
def R_solve(R, b):
"""
returns x
This function solves the triangular system R x = b for the N-by-N
matrix R.
"""
x = zeros(R.shape[1], get_typecode(R))
_gslwrap.gsl_linalg_QR_QRsolve(R, b, x)
return x
#
# SVD Singular Value Decomposition
#
[docs]
def SV_decomp(A):
"""
returns (U, V, S)
This function factorizes the M-by-N matrix A into the singular value
decomposition A = U S V^T. The diagonal elements of the singular value
matrix S are stored in the vector S. The singular values are non-negative
and form a non-increasing sequence from S_1 to S_N. The matrix V
contains the elements of V in untransposed form. To form the product
U S V^T it is necessary to take the transpose of V. A workspace of
length N is required in work.
This routine uses the Golub-Reinsch SVD algorithm.
"""
code = get_typecode(A)
n = A.shape[1]
u = array_typed_copy(A, code)
s = zeros(n, code)
v = zeros((n, n), code)
work = zeros(A.shape[1], code)
_gslwrap.gsl_linalg_SV_decomp(u, v, s, work)
return (u, v, s)
[docs]
def SV_decomp_mod(A):
"""
returns (u, v, s)
This function computes the SVD using the modified Golub-Reinsch
algorithm, which is faster for M>>N. It requires the vector work
and the N-by-N matrix X as additional working space.
"""
code = get_typecode(A)
n = A.shape[1]
u = array_typed_copy(A, code)
s = zeros(n, code)
v = zeros((n, n), code)
x = zeros((n, n), code)
work = zeros(A.shape[1], code)
_gslwrap.gsl_linalg_SV_decomp_mod(u, x, v, s, work)
return (u, v, s)
[docs]
def SV_decomp_jacobi(A):
"""
returns (u, v, s)
This function computes the SVD using one-sided Jacobi orthogonalization
(see references for details). The Jacobi method can compute singular
values to higher relative accuracy than Golub-Reinsch algorithms.
"""
code = get_typecode(A)
n = A.shape[1]
u = array_typed_copy(A, code)
s = zeros(n, code)
v = zeros((n, n), code)
_gslwrap.gsl_linalg_SV_decomp_jacobi(u, v, s)
return (u, v, s)
[docs]
def SV_solve(U, V, S, b):
"""
returns x
This function solves the system A x = b using the singular value
decomposition (U, S, V) of A given by gsl_linalg_SV_decomp.
Only non-zero singular values are used in computing the solution.
The parts of the solution corresponding to singular values of zero
are ignored. Other singular values can be edited out by setting them
to zero before calling this function.
In the over-determined case where A has more rows than columns the
system is solved in the least squares sense, returning the solution x
which minimizes ||A x - b||_2.
"""
x = zeros(U.shape[1], get_typecode(b))
_gslwrap.gsl_linalg_SV_solve(U, V, S, b, x)
return x
#
# Cholesky Decomposition
#
[docs]
def cholesky_decomp(A):
"""
Function: int gsl_linalg_cholesky_decomp (gsl_matrix * A)
This function factorizes the positive-definite square matrix A into
the Cholesky decomposition A = L L^T. On output the diagonal and
lower triangular part of the input matrix A contain the matrix L.
The upper triangular part of the input matrix contains L^T, the diagonal
terms being identical for both L and L^T. If the matrix is not
positive-definite then the decomposition will fail, returning the
error code GSL_EDOM.
"""
An = array_typed_copy(A)
_gslwrap.gsl_linalg_cholesky_decomp(An)
return An
[docs]
def cholesky_unpack(L):
"""
returns (L, L^T)
This function splits the matrix L into the the upper matrix L^T and
the lower matrix L. The diagonal of L is the identical for both.
"""
lt = zeros(L.shape, get_typecode(L))
l = zeros(L.shape, get_typecode(L))
for i in range(L.shape[0]):
lt[i, i: ] = L[i, i:]
l[i, 0:i+1] = L[i, :i+1]
return (l, lt)
[docs]
def cholesky_solve(cholesky, b):
"""
returns x
This function solves the system A x = b using the Cholesky decomposition
of A into the matrix cholesky given by cholesky_decomp.
"""
x = zeros(b.shape, get_typecode(b))
_gslwrap.gsl_linalg_cholesky_solve(cholesky, b, x)
return x
#
# Tridiagonal Decomposition of Real Symmetric Matrices
#
# A symmetric matrix A can be factorized by similarity transformations
# into the form,
#
# A = Q T Q^T
#
# where Q is an orthogonal matrix and T is a symmetric tridiagonal matrix.
#
[docs]
def symmtd_decomp(A):
"""
returns (QT, tau)
This function factorizes the symmetric square matrix A into the
symmetric tridiagonal decomposition Q T Q^T. On output the diagonal
and subdiagonal part of the input matrix A contain the tridiagonal
matrix T. The remaining lower triangular part of the input matrix
contains the Householder vectors which, together with the Householder
coefficients tau, encode the orthogonal matrix Q. This storage scheme
is the same as used by LAPACK. The upper triangular part of A is not
referenced.
"""
code = get_typecode(A)
QT = array_typed_copy(A, code)
tau = zeros(A.shape[0]-1, code)
_gslwrap.gsl_linalg_symmtd_decomp(QT, tau)
return (QT, tau)
[docs]
def symmtd_unpack(A, tau):
"""
returns (Q, diag, subdiag)
This function unpacks the encoded symmetric tridiagonal decomposition
(A, tau) obtained from gsl_linalg_symmtd_decomp into the orthogonal
matrix Q, the vector of diagonal elements diag and the vector of
subdiagonal elements subdiag.
"""
n = A.shape[0]
code = get_typecode(A)
Q = zeros([n,n], code)
diag = zeros((n,), code)
subdiag = zeros((n-1,), code)
_gslwrap.gsl_linalg_symmtd_unpack(A, tau, Q, diag, subdiag)
return (Q, diag, subdiag)
def symmtd_unpack_T(A):
"""
returns (diag, subdiag)
This function unpacks the diagonal and subdiagonal of the encoded
symmetric tridiagonal decomposition (A, tau) obtained from
gsl_linalg_symmtd_decomp into the vectors diag and subdiag.
"""
n = A.shape[0]
code = get_typecode(A)
diag = zeros((n,), code)
subdiag = zeros((n-1,), code)
_gslwrap.gsl_linalg_symmtd_unpack_T(A, diag, subdiag)
return (diag, subdiag)
[docs]
def symmtd_unpack_diag(diag, subdiag):
"""
returns T
This functions unpacks the tridiagonal matrix T of the diagonal
and subdiagonal obtained from gsl_linalg_symmtd_unpack[_T]
"""
n = diag.shape[0]
T = numx.identity(n)*diag
sub = numx.identity(n-1)*subdiag
sub1 = numx.concatenate(
(zeros((1,n)), numx.concatenate((sub, zeros((n-1,1))), 1)))
sub2 = numx.concatenate(
(numx.concatenate((zeros([n-1,1]),sub), 1), zeros([1,n])))
T = sub1 + sub2 + T
return T
#
# Tridiagonal Decomposition of Hermitian Matrices
#
# A hermitian matrix A can be factorized by similarity transformations
# into the form,
#
# A = U T U^T
#
# where U is an unitary matrix and T is a real symmetric tridiagonal matrix.
#
[docs]
def hermtd_decomp(A):
"""
returns (QT, tau)
This function factorizes the hermitian matrix A into the symmetric
tridiagonal decomposition U T U^T. On output the real parts of the
diagonal and subdiagonal part of the input matrix A contain the
tridiagonal matrix T. The remaining lower triangular part of the input
matrix contains the Householder vectors which, together with the
Householder coefficients tau, encode the orthogonal matrix Q. This
storage scheme is the same as used by LAPACK. The upper triangular
part of A and imaginary parts of the diagonal are not referenced.
"""
code = get_typecode(A)
QT = array_typed_copy(A, code)
tau = zeros(A.shape[0]-1, code)
_gslwrap.gsl_linalg_hermtd_decomp(QT, tau)
return (QT, tau)
[docs]
def hermtd_unpack(A, tau):
"""
returns (Q, diag, subdiag)
This function unpacks the encoded tridiagonal decomposition (A, tau)
obtained from gsl_linalg_hermtd_decomp into the unitary matrix U, the
real vector of diagonal elements diag and the real vector of subdiagonal
elements subdiag.
"""
n = A.shape[0]
code = get_typecode(A)
Q = zeros([n,n], code)
diag = zeros((n,), Float)
subdiag = zeros((n-1,), Float)
_gslwrap.gsl_linalg_hermtd_unpack(A, tau, Q, diag, subdiag)
return (Q, diag, subdiag)
[docs]
def symmtd_unpack_T(A):
"""
returns (diag, subdiag)
This function unpacks the diagonal and subdiagonal of the encoded
tridiagonal decomposition (A, tau) obtained from
gsl_linalg_hermtd_decomp into the real vectors diag and subdiag.
"""
n = A.shape[0]
diag = zeros((n,), Float)
subdiag = zeros((n-1,), Float)
_gslwrap.gsl_linalg_hermtd_unpack_T(A, diag, subdiag)
return (diag, subdiag)
[docs]
def hermtd_unpack_diag(diag, subdiag):
"""
returns T
This functions unpacks the tridiagonal matrix T of the diagonal
and subdiagonal obtained from gsl_linalg_hermtd_unpack[_T]
"""
return symmtd_unpack_diag(diag, subdiag)
#
# Bidiagonalization
#
# A general matrix A can be factorized by similarity transformations into
# the form,
#
# A = U B V^T
#
# where U and V are orthogonal matrices and B is a N-by-N bidiagonal
# matrix with non-zero entries only on the diagonal and superdiagonal.
# The size of U is M-by-N and the size of V is N-by-N.
#
[docs]
def bidiag_decomp(A):
"""
returns (BUV, tau_U, tau_V)
This function factorizes the M-by-N matrix A into bidiagonal
form U B V^T. The diagonal and superdiagonal of the matrix B are
stored in the diagonal and superdiagonal of BUV. The orthogonal matrices
U and V are stored as compressed Householder vectors in the remaining
elements of BUV. The Householder coefficients are stored in the vectors
tau_U and tau_V. The length of tau_U must equal the number of elements
in the diagonal of A and the length of tau_V should be one element shorter.
"""
n = min(A.shape)
code = get_typecode(A)
BUV = array_typed_copy(A, code)
tau_U = zeros(n, code)
tau_V = zeros(n-1, code)
_gslwrap.gsl_linalg_bidiag_decomp(BUV, tau_U, tau_V)
return (BUV, tau_U, tau_V)
[docs]
def bidiag_unpack(A, tau_U, tau_V):
"""
returns (U, V, diag, superdiag)
This function unpacks the bidiagonal decomposition of A given by
gsl.linalg.bidiag_decomp, (A, tau_U, tau_V) into the separate
orthogonal matrices U, V and the diagonal vector diag and
superdiagonal superdiag.
"""
(m,n) = A.shape
code = get_typecode(A)
U = zeros([m,n], code)
V = zeros([n,n], code)
diag = zeros(n, code)
superdiag = zeros(n-1, code)
_gslwrap.gsl_linalg_bidiag_unpack(A, tau_U, U, tau_V, V, diag, superdiag)
return (U, V, diag, superdiag)
[docs]
def bidiag_unpack_B(A):
"""
returns (diag, superdiag)
This function unpacks the diagonal and superdiagonal of the bidiagonal
decomposition of A given by gsl_linalg_bidiag_decomp, into the diagonal
vector diag and superdiagonal vector superdiag.
"""
raise NotImplementedError("the GSL function for this is buggy!" )
n = n = A.shape[1]
code = get_typecode(A)
diag = zeros(n, code)
superdiag = zeros(n-1, code)
_gslwrap.gsl_linalg_bidiag_unpack_B(A, diag, superdiag)
return (diag, superdiag)
[docs]
def bidiag_unpack_diag(diag, superdiag):
"""
returns B
This functions unpacks the bidiagonal matrix T of the diagonal
and superdiagonal obtained from gsl_linalg_bidiag_unpack[_B]
"""
n = diag.shape[0]
B = numx.identity(n)*diag
sub = numx.identity(n-1)*superdiag
sub2 = numx.concatenate(
(numx.concatenate((zeros([n-1,1]),sub), 1), zeros([1,n])))
B = sub2 + B
return B
#
# Householder solver for linear systems
#
[docs]
def HH_solve(A,b):
"""
returns x
This function solves the system A x = b directly using Householder
transformations. On output the solution is stored in x and b is not
modified.
"""
code = get_typecode(A)
x = zeros(A.shape[1], code)
An = array_typed_copy(A, code)
_gslwrap.gsl_linalg_HH_solve(An, b, x)
return x
#
# Tridiagonal Systems
#
[docs]
def solve_tridiag(diag, e, f, b):
"""
returns x
This function solves the general N-by-N system A x = b where A is
tridiagonal. The form of A for the 4-by-4 case is shown below,
A = ( d_0 e_0 )
( f_0 d_1 e_1 )
( f_1 d_2 e_2 )
( f_2 d_3 )
"""
x = zeros(diag.shape, get_typecode(diag))
_gslwrap.gsl_linalg_solve_tridiag(diag, e, f, b, x)
return x
[docs]
def solve_symm_tridiag(diag, e, b):
"""
returns x
This function solves the general N-by-N system A x = b where A is
symmetric tridiagonal. The form of A for the 4-by-4 case is shown below,
A = ( d_0 e_0 )
( e_0 d_1 e_1 )
( e_1 d_2 e_2 )
( e_2 d_3 )
"""
x = zeros(diag.shape, get_typecode(diag))
_gslwrap.gsl_linalg_solve_symm_tridiag(diag, e, b, x)
return x
[docs]
def solve_symm_cyc_tridiag(diag, e, b):
"""
returns x
This function solves the general N-by-N system A x = b where A is
symmetric cyclic tridiagonal. The form of A for the 4-by-4 case is
shown below,
A = ( d_0 e_0 e_3 )
( e_0 d_1 e_1 )
( e_1 d_2 e_2 )
( e_3 e_2 d_3 )
"""
x = zeros(diag.shape, get_typecode(diag))
_gslwrap.gsl_linalg_solve_symm_cyc_tridiag(diag, e, b, x)
return x