[Scipy-svn] r3952 - trunk/scipy/linalg
scipy-svn at scipy.org
scipy-svn at scipy.org
Wed Feb 20 20:57:49 EST 2008
Author: wnbell
Date: 2008-02-20 19:57:37 -0600 (Wed, 20 Feb 2008)
New Revision: 3952
Modified:
trunk/scipy/linalg/basic.py
trunk/scipy/linalg/blas.py
trunk/scipy/linalg/decomp.py
trunk/scipy/linalg/matfuncs.py
Log:
applied patch from user pv
Reformat and ReSTify scipy.linalg docstrings
resolves ticket #596
Modified: trunk/scipy/linalg/basic.py
===================================================================
--- trunk/scipy/linalg/basic.py 2008-02-20 05:39:03 UTC (rev 3951)
+++ trunk/scipy/linalg/basic.py 2008-02-21 01:57:37 UTC (rev 3952)
@@ -25,22 +25,34 @@
def lu_solve((lu, piv), b, trans=0, overwrite_b=0):
- """ lu_solve((lu, piv), b, trans=0, overwrite_b=0) -> x
+ """Solve an equation system, a x = b, given the LU factorization of a
+
+ Parameters
+ ----------
+ (lu, piv)
+ Factorization of the coefficient matrix a, as given by lu_factor
+ b : array
+ Right-hand side
+ trans : {0, 1, 2}
+ Type of system to solve:
- Solve a system of equations given a previously factored matrix
+ ===== =========
+ trans system
+ ===== =========
+ 0 a x = b
+ 1 a^T x = b
+ 2 a^H x = b
+ ===== =========
- Inputs:
+ Returns
+ -------
+ x : array
+ Solution to the system
- (lu,piv) -- The factored matrix, a (the output of lu_factor)
- b -- a set of right-hand sides
- trans -- type of system to solve:
- 0 : a * x = b (no transpose)
- 1 : a^T * x = b (transpose)
- 2 a^H * x = b (conjugate transpose)
-
- Outputs:
-
- x -- the solution to the system
+ See also
+ --------
+ lu_factor : LU factorize a matrix
+
"""
b1 = asarray_chkfinite(b)
overwrite_b = overwrite_b or (b1 is not b and not hasattr(b,'__array__'))
@@ -54,18 +66,24 @@
'illegal value in %-th argument of internal gesv|posv'%(-info)
def cho_solve((c, lower), b, overwrite_b=0):
- """ cho_solve((c, lower), b, overwrite_b=0) -> x
+ """Solve an equation system, a x = b, given the Cholesky factorization of a
- Solve a system of equations given a previously cholesky factored matrix
+ Parameters
+ ----------
+ (c, lower)
+ Cholesky factorization of a, as given by cho_factor
+ b : array
+ Right-hand side
- Inputs:
+ Returns
+ -------
+ x : array
+ The solution to the system a x = b
- (c,lower) -- The factored matrix, a (the output of cho_factor)
- b -- a set of right-hand sides
-
- Outputs:
-
- x -- the solution to the system a*x = b
+ See also
+ --------
+ cho_factor : Cholesky factorization of a matrix
+
"""
b1 = asarray_chkfinite(b)
overwrite_b = overwrite_b or (b1 is not b and not hasattr(b,'__array__'))
@@ -81,22 +99,29 @@
# Linear equations
def solve(a, b, sym_pos=0, lower=0, overwrite_a=0, overwrite_b=0,
debug = 0):
- """ solve(a, b, sym_pos=0, lower=0, overwrite_a=0, overwrite_b=0) -> x
+ """Solve the equation a x = b for x
- Solve a linear system of equations a * x = b for x.
+ Parameters
+ ----------
+ a : array, shape (M, M)
+ b : array, shape (M,) or (M, N)
+ sym_pos : boolean
+ Assume a is symmetric and positive definite
+ lower : boolean
+ Use only data contained in the lower triangle of a, if sym_pos is true.
+ Default is to use upper triangle.
+ overwrite_a : boolean
+ Allow overwriting data in a (may enhance performance)
+ overwrite_b : boolean
+ Allow overwriting data in b (may enhance performance)
+
+ Returns
+ -------
+ x : array, shape (M,) or (M, N) depending on b
+ Solution to the system a x = b
- Inputs:
-
- a -- An N x N matrix.
- b -- An N x nrhs matrix or N vector.
- sym_pos -- Assume a is symmetric and positive definite.
- lower -- Assume a is lower triangular, otherwise upper one.
- Only used if sym_pos is true.
- overwrite_y - Discard data in y, where y is a or b.
-
- Outputs:
-
- x -- The solution to the system a * x = b
+ Raises LinAlgError if a is singular
+
"""
a1, b1 = map(asarray_chkfinite,(a,b))
if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
@@ -129,29 +154,37 @@
def solve_banded((l,u), ab, b, overwrite_ab=0, overwrite_b=0,
debug = 0):
- """ solve_banded((l,u), ab, b, overwrite_ab=0, overwrite_b=0) -> x
+ """Solve the equation a x = b for x, assuming a is banded matrix.
- Solve a linear system of equations a * x = b for x where
- a is a banded matrix stored in diagonal orded form
+ The matrix a is stored in ab using the matrix diagonal orded form::
- * * a1u
+ ab[u + i - j, j] == a[i,j]
- * a12 a23 ...
- a11 a22 a33 ...
- a21 a32 a43 ...
- .
- al1 .. *
+ Example of ab (shape of a is (6,6), u=1, l=2)::
- Inputs:
+ * a01 a12 a23 a34 a45
+ a00 a11 a22 a33 a44 a55
+ a10 a21 a32 a43 a54 *
+ a20 a31 a42 a53 * *
- (l,u) -- number of non-zero lower and upper diagonals, respectively.
- a -- An N x (l+u+1) matrix.
- b -- An N x nrhs matrix or N vector.
- overwrite_y - Discard data in y, where y is ab or b.
+ Parameters
+ ----------
+ (l, u) : (integer, integer)
+ Number of non-zero lower and upper diagonals
+ ab : array, shape (l+u+1, M)
+ Banded matrix
+ b : array, shape (M,) or (M, K)
+ Right-hand side
+ overwrite_ab : boolean
+ Discard data in ab (may enhance performance)
+ overwrite_b : boolean
+ Discard data in b (may enhance performance)
- Outputs:
-
- x -- The solution to the system a * x = b
+ Returns
+ -------
+ x : array, shape (M,) or (M, K)
+ The solution to the system a x = b
+
"""
a1, b1 = map(asarray_chkfinite,(ab,b))
overwrite_b = overwrite_b or (b1 is not b and not hasattr(b,'__array__'))
@@ -171,34 +204,48 @@
def solveh_banded(ab, b, overwrite_ab=0, overwrite_b=0,
lower=0):
- """ solveh_banded(ab, b, overwrite_ab=0, overwrite_b=0) -> c, x
+ """Solve equation a x = b. a is Hermitian positive-definite banded matrix.
- Solve a linear system of equations a * x = b for x where
- a is a banded symmetric or Hermitian positive definite
- matrix stored in lower diagonal ordered form (lower=1)
+ The matrix a is stored in ab either in lower diagonal or upper
+ diagonal ordered form:
- a11 a22 a33 a44 a55 a66
- a21 a32 a43 a54 a65 *
- a31 a42 a53 a64 * *
+ ab[u + i - j, j] == a[i,j] (if upper form; i <= j)
+ ab[ i - j, j] == a[i,j] (if lower form; i >= j)
- or upper diagonal ordered form
+ Example of ab (shape of a is (6,6), u=2)::
- * * a31 a42 a53 a64
- * a21 a32 a43 a54 a65
- a11 a22 a33 a44 a55 a66
+ upper form:
+ * * a02 a13 a24 a35
+ * a01 a12 a23 a34 a45
+ a00 a11 a22 a33 a44 a55
+
+ lower form:
+ a00 a11 a22 a33 a44 a55
+ a10 a21 a32 a43 a54 *
+ a20 a31 a42 a53 * *
- Inputs:
+ Cells marked with * are not used.
- ab -- An N x l
- b -- An N x nrhs matrix or N vector.
- overwrite_y - Discard data in y, where y is ab or b.
- lower - is ab in lower or upper form?
+ Parameters
+ ----------
+ ab : array, shape (M, u + 1)
+ Banded matrix
+ b : array, shape (M,) or (M, K)
+ Right-hand side
+ overwrite_ab : boolean
+ Discard data in ab (may enhance performance)
+ overwrite_b : boolean
+ Discard data in b (may enhance performance)
+ lower : boolean
+ Is the matrix in the lower form. (Default is upper form)
- Outputs:
-
- c: the Cholesky factorization of ab
- x: the solution to ab * x = b
-
+ Returns
+ -------
+ c : array, shape (M, u+1)
+ Cholesky factorization of a, in the same banded format as ab
+ x : array, shape (M,) or (M, K)
+ The solution to the system a x = b
+
"""
ab, b = map(asarray_chkfinite,(ab,b))
@@ -215,32 +262,40 @@
'illegal value in %d-th argument of internal pbsv'%(-info)
def cholesky_banded(ab, overwrite_ab=0, lower=0):
- """ cholesky_banded(ab, overwrite_ab=0, lower=0) -> c
-
- Compute the Cholesky decomposition of a
- banded symmetric or Hermitian positive definite
- matrix stored in lower diagonal ordered form (lower=1)
-
- a11 a22 a33 a44 a55 a66
- a21 a32 a43 a54 a65 *
- a31 a42 a53 a64 * *
-
- or upper diagonal ordered form
-
- * * a31 a42 a53 a64
- * a21 a32 a43 a54 a65
- a11 a22 a33 a44 a55 a66
-
- Inputs:
-
- ab -- An N x l
- overwrite_ab - Discard data in ab
- lower - is ab in lower or upper form?
-
- Outputs:
-
- c: the Cholesky factorization of ab
-
+ """Cholesky decompose a banded Hermitian positive-definite matrix
+
+ The matrix a is stored in ab either in lower diagonal or upper
+ diagonal ordered form:
+
+ ab[u + i - j, j] == a[i,j] (if upper form; i <= j)
+ ab[ i - j, j] == a[i,j] (if lower form; i >= j)
+
+ Example of ab (shape of a is (6,6), u=2)::
+
+ upper form:
+ * * a02 a13 a24 a35
+ * a01 a12 a23 a34 a45
+ a00 a11 a22 a33 a44 a55
+
+ lower form:
+ a00 a11 a22 a33 a44 a55
+ a10 a21 a32 a43 a54 *
+ a20 a31 a42 a53 * *
+
+ Parameters
+ ----------
+ ab : array, shape (M, u + 1)
+ Banded matrix
+ overwrite_ab : boolean
+ Discard data in ab (may enhance performance)
+ lower : boolean
+ Is the matrix in the lower form. (Default is upper form)
+
+ Returns
+ -------
+ c : array, shape (M, u+1)
+ Cholesky factorization of a, in the same banded format as ab
+
"""
ab = asarray_chkfinite(ab)
@@ -259,9 +314,30 @@
# matrix inversion
def inv(a, overwrite_a=0):
- """ inv(a, overwrite_a=0) -> a_inv
+ """Compute the inverse of a matrix.
+
+ Parameters
+ ----------
+ a : array-like, shape (M, M)
+ Matrix to be inverted
+
+ Returns
+ -------
+ ainv : array-like, shape (M, M)
+ Inverse of the matrix a
- Return inverse of square matrix a.
+ Raises LinAlgError if a is singular
+
+ Examples
+ --------
+ >>> a = array([[1., 2.], [3., 4.]])
+ >>> inv(a)
+ array([[-2. , 1. ],
+ [ 1.5, -0.5]])
+ >>> dot(a, inv(a))
+ array([[ 1., 0.],
+ [ 0., 1.]])
+
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
@@ -312,35 +388,39 @@
## matrix and Vector norm
import decomp
def norm(x, ord=None):
- """ norm(x, ord=None) -> n
+ """Matrix or vector norm.
- Matrix or vector norm.
+ Parameters
+ ----------
+ x : array, shape (M,) or (M, N)
+ ord : number, or {None, 1, -1, 2, -2, inf, -inf, 'fro'}
+ Order of the norm:
- Inputs:
+ ===== ============================ ==========================
+ ord norm for matrices norm for vectors
+ ===== ============================ ==========================
+ None Frobenius norm 2-norm
+ 'fro' Frobenius norm -
+ inf max(sum(abs(x), axis=1)) max(abs(x))
+ -inf min(sum(abs(x), axis=1)) min(abs(x))
+ 1 max(sum(abs(x), axis=0)) as below
+ -1 min(sum(abs(x), axis=0)) as below
+ 2 2-norm (largest sing. value) as below
+ -2 smallest singular value as below
+ other - sum(abs(x)**ord)**(1./ord)
+ ===== ============================ ==========================
+
+ Returns
+ -------
+ n : float
+ Norm of the matrix or vector
- x -- a rank-1 (vector) or rank-2 (matrix) array
- ord -- the order of the norm.
-
- Comments:
- For arrays of any rank, if ord is None:
- calculate the square norm (Euclidean norm for vectors, Frobenius norm for matrices)
-
- For vectors ord can be any real number including Inf or -Inf.
- ord = Inf, computes the maximum of the magnitudes
- ord = -Inf, computes minimum of the magnitudes
- ord is finite, computes sum(abs(x)**ord,axis=0)**(1.0/ord)
-
- For matrices ord can only be one of the following values:
- ord = 2 computes the largest singular value
- ord = -2 computes the smallest singular value
- ord = 1 computes the largest column sum of absolute values
- ord = -1 computes the smallest column sum of absolute values
- ord = Inf computes the largest row sum of absolute values
- ord = -Inf computes the smallest row sum of absolute values
- ord = 'fro' computes the frobenius norm sqrt(sum(diag(X.H * X),axis=0))
-
- For values ord < 0, the result is, strictly speaking, not a
- mathematical 'norm', but it may still be useful for numerical purposes.
+ Notes
+ -----
+ For values ord < 0, the result is, strictly speaking, not a
+ mathematical 'norm', but it may still be useful for numerical
+ purposes.
+
"""
x = asarray_chkfinite(x)
if ord is None: # check the default case first and handle it immediately
@@ -382,9 +462,20 @@
### Determinant
def det(a, overwrite_a=0):
- """ det(a, overwrite_a=0) -> d
+ """Compute the determinant of a matrix
- Return determinant of a square matrix.
+ Parameters
+ ----------
+ a : array, shape (M, M)
+
+ Returns
+ -------
+ det : float or complex
+ Determinant of a
+
+ Notes
+ -----
+ The determinant is computed via LU factorization, LAPACK routine z/dgetrf.
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
@@ -399,25 +490,38 @@
### Linear Least Squares
def lstsq(a, b, cond=None, overwrite_a=0, overwrite_b=0):
- """ lstsq(a, b, cond=None, overwrite_a=0, overwrite_b=0) -> x,resids,rank,s
-
- Return least-squares solution of a * x = b.
-
- Inputs:
-
- a -- An M x N matrix.
- b -- An M x nrhs matrix or M vector.
- cond -- Used to determine effective rank of a.
-
- Outputs:
-
- x -- The solution (N x nrhs matrix) to the minimization problem:
- 2-norm(| b - a * x |) -> min
- resids -- The residual sum-of-squares for the solution matrix x
- (only if M>N and rank==N).
- rank -- The effective rank of a.
- s -- Singular values of a in decreasing order. The condition number
- of a is abs(s[0]/s[-1]).
+ """Compute least-squares solution to equation :m:`a x = b`
+
+ Compute a vector x such that the 2-norm :m:`|b - a x|` is minimised.
+
+ Parameters
+ ----------
+ a : array, shape (M, N)
+ b : array, shape (M,) or (M, K)
+ cond : float
+ Cutoff for 'small' singular values; used to determine effective
+ rank of a. Singular values smaller than rcond*largest_singular_value
+ are considered zero.
+ overwrite_a : boolean
+ Discard data in a (may enhance performance)
+ overwrite_b : boolean
+ Discard data in b (may enhance performance)
+
+ Returns
+ -------
+ x : array, shape (N,) or (N, K) depending on shape of b
+ Least-squares solution
+ residues : array, shape () or (1,) or (K,)
+ Sums of residues, squared 2-norm for each column in :m:`b - a x`
+ If rank of matrix a is < N or > M this is an empty array.
+ If b was 1-d, this is an (1,) shape array, otherwise the shape is (K,)
+ rank : integer
+ Effective rank of matrix a
+ s : array, shape (min(M,N),)
+ Singular values of a. The condition number of a is abs(s[0]/s[-1]).
+
+ Raises LinAlgError if computation does not converge
+
"""
a1, b1 = map(asarray_chkfinite,(a,b))
if len(a1.shape) != 2:
@@ -457,9 +561,36 @@
def pinv(a, cond=None, rcond=None):
- """ pinv(a, rcond=None) -> a_pinv
+ """Compute the (Moore-Penrose) pseudo-inverse of a matrix.
+
+ Calculate a generalized inverse of a matrix using a least-squares
+ solver.
+
+ Parameters
+ ----------
+ a : array, shape (M, N)
+ Matrix to be pseudo-inverted
+ cond, rcond : float
+ Cutoff for 'small' singular values in the least-squares solver.
+ Singular values smaller than rcond*largest_singular_value are
+ considered zero.
+
+ Returns
+ -------
+ B : array, shape (N, M)
+
+ Raises LinAlgError if computation does not converge
- Compute generalized inverse of A using least-squares solver.
+ Examples
+ --------
+ >>> from numpy import *
+ >>> a = random.randn(9, 6)
+ >>> B = linalg.pinv(a)
+ >>> allclose(a, dot(a, dot(B, a)))
+ True
+ >>> allclose(B, dot(B, dot(a, B)))
+ True
+
"""
a = asarray_chkfinite(a)
b = numpy.identity(a.shape[0], dtype=a.dtype)
@@ -473,9 +604,39 @@
_array_precision = {'f': 0, 'd': 1, 'F': 0, 'D': 1}
def pinv2(a, cond=None, rcond=None):
- """ pinv2(a, rcond=None) -> a_pinv
+ """Compute the (Moore-Penrose) pseudo-inverse of a matrix.
+
+ Calculate a generalized inverse of a matrix using its
+ singular-value decomposition and including all 'large' singular
+ values.
+
+ Parameters
+ ----------
+ a : array, shape (M, N)
+ Matrix to be pseudo-inverted
+ cond, rcond : float or None
+ Cutoff for 'small' singular values.
+ Singular values smaller than rcond*largest_singular_value are
+ considered zero.
- Compute the generalized inverse of A using svd.
+ If None or -1, suitable machine precision is used.
+
+ Returns
+ -------
+ B : array, shape (N, M)
+
+ Raises LinAlgError if SVD computation does not converge
+
+ Examples
+ --------
+ >>> from numpy import *
+ >>> a = random.randn(9, 6)
+ >>> B = linalg.pinv2(a)
+ >>> allclose(a, dot(a, dot(B, a)))
+ True
+ >>> allclose(B, dot(B, dot(a, B)))
+ True
+
"""
a = asarray_chkfinite(a)
u, s, vh = decomp.svd(a)
@@ -498,8 +659,37 @@
#-----------------------------------------------------------------------------
def tri(N, M=None, k=0, dtype=None):
- """ returns a N-by-M matrix where all the diagonals starting from
- lower left corner up to the k-th are all ones.
+ """Construct (N, M) matrix filled with ones at and below the k-th diagonal.
+
+ The matrix has A[i,j] == 1 for i <= j + k
+
+ Parameters
+ ----------
+ N : integer
+ M : integer
+ Size of the matrix. If M is None, M == N is assumed.
+ k : integer
+ Number of subdiagonal below which matrix is filled with ones.
+ k == 0 is the main diagonal, k < 0 subdiagonal and k > 0 superdiagonal.
+ dtype : dtype
+ Data type of the matrix.
+
+ Returns
+ -------
+ A : array, shape (N, M)
+
+ Examples
+ --------
+ >>> from scipy.linalg import tri
+ >>> tri(3, 5, 2, dtype=int)
+ array([[1, 1, 1, 0, 0],
+ [1, 1, 1, 1, 0],
+ [1, 1, 1, 1, 1]])
+ >>> tri(3, 5, -1, dtype=int)
+ array([[0, 0, 0, 0, 0],
+ [1, 0, 0, 0, 0],
+ [1, 1, 0, 0, 0]])
+
"""
if M is None: M = N
if type(M) == type('d'):
@@ -514,8 +704,29 @@
return m.astype(dtype)
def tril(m, k=0):
- """ returns the elements on and below the k-th diagonal of m. k=0 is the
- main diagonal, k > 0 is above and k < 0 is below the main diagonal.
+ """Construct a copy of a matrix with elements above the k-th diagonal zeroed.
+
+ Parameters
+ ----------
+ m : array
+ Matrix whose elements to return
+ k : integer
+ Diagonal above which to zero elements.
+ k == 0 is the main diagonal, k < 0 subdiagonal and k > 0 superdiagonal.
+
+ Returns
+ -------
+ A : array, shape m.shape, dtype m.dtype
+
+ Examples
+ --------
+ >>> from scipy.linalg import tril
+ >>> tril([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], -1)
+ array([[ 0, 0, 0],
+ [ 4, 0, 0],
+ [ 7, 8, 0],
+ [10, 11, 12]])
+
"""
svsp = getattr(m,'spacesaver',lambda:0)()
m = asarray(m)
@@ -524,8 +735,29 @@
return out
def triu(m, k=0):
- """ returns the elements on and above the k-th diagonal of m. k=0 is the
- main diagonal, k > 0 is above and k < 0 is below the main diagonal.
+ """Construct a copy of a matrix with elements below the k-th diagonal zeroed.
+
+ Parameters
+ ----------
+ m : array
+ Matrix whose elements to return
+ k : integer
+ Diagonal below which to zero elements.
+ k == 0 is the main diagonal, k < 0 subdiagonal and k > 0 superdiagonal.
+
+ Returns
+ -------
+ A : array, shape m.shape, dtype m.dtype
+
+ Examples
+ --------
+ >>> from scipy.linalg import tril
+ >>> triu([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], -1)
+ array([[ 1, 2, 3],
+ [ 4, 5, 6],
+ [ 0, 8, 9],
+ [ 0, 0, 12]])
+
"""
svsp = getattr(m,'spacesaver',lambda:0)()
m = asarray(m)
@@ -534,16 +766,36 @@
return out
def toeplitz(c,r=None):
- """ Construct a toeplitz matrix (i.e. a matrix with constant diagonals).
-
- Description:
-
- toeplitz(c,r) is a non-symmetric Toeplitz matrix with c as its first
- column and r as its first row.
-
- toeplitz(c) is a symmetric (Hermitian) Toeplitz matrix (r=c).
-
- See also: hankel
+ """Construct a Toeplitz matrix.
+
+ The Toepliz matrix has constant diagonals, c as its first column,
+ and r as its first row (if not given, r == c is assumed).
+
+ Parameters
+ ----------
+ c : array
+ First column of the matrix
+ r : array
+ First row of the matrix. If None, r == c is assumed.
+
+ Returns
+ -------
+ A : array, shape (len(c), len(r))
+ Constructed Toeplitz matrix.
+ dtype is the same as (c[0] + r[0]).dtype
+
+ Examples
+ --------
+ >>> from scipy.linalg import toeplitz
+ >>> toeplitz([1,2,3], [1,4,5,6])
+ array([[1, 4, 5, 6],
+ [2, 1, 4, 5],
+ [3, 2, 1, 4]])
+
+ See also
+ --------
+ hankel : Hankel matrix
+
"""
isscalar = numpy.isscalar
if isscalar(c) or isscalar(r):
@@ -566,17 +818,37 @@
def hankel(c,r=None):
- """ Construct a hankel matrix (i.e. matrix with constant anti-diagonals).
-
- Description:
-
- hankel(c,r) is a Hankel matrix whose first column is c and whose
- last row is r.
-
- hankel(c) is a square Hankel matrix whose first column is C.
- Elements below the first anti-diagonal are zero.
-
- See also: toeplitz
+ """Construct a Hankel matrix.
+
+ The Hankel matrix has constant anti-diagonals, c as its first column,
+ and r as its last row (if not given, r == 0 os assumed).
+
+ Parameters
+ ----------
+ c : array
+ First column of the matrix
+ r : array
+ Last row of the matrix. If None, r == 0 is assumed.
+
+ Returns
+ -------
+ A : array, shape (len(c), len(r))
+ Constructed Hankel matrix.
+ dtype is the same as (c[0] + r[0]).dtype
+
+ Examples
+ --------
+ >>> from scipy.linalg import hankel
+ >>> hankel([1,2,3,4], [4,7,7,8,9])
+ array([[1, 2, 3, 4, 7],
+ [2, 3, 4, 7, 7],
+ [3, 4, 7, 7, 8],
+ [4, 7, 7, 8, 9]])
+
+ See also
+ --------
+ toeplitz : Toeplitz matrix
+
"""
isscalar = numpy.isscalar
if isscalar(c) or isscalar(r):
@@ -599,12 +871,32 @@
return map(Matrix,args)
def kron(a,b):
- """kronecker product of a and b
+ """Kronecker product of a and b.
- Kronecker product of two matrices is block matrix
- [[ a[ 0 ,0]*b, a[ 0 ,1]*b, ... , a[ 0 ,n-1]*b ],
- [ ... ... ],
- [ a[m-1,0]*b, a[m-1,1]*b, ... , a[m-1,n-1]*b ]]
+ The result is the block matrix::
+
+ a[0,0]*b a[0,1]*b ... a[0,-1]*b
+ a[1,0]*b a[1,1]*b ... a[1,-1]*b
+ ...
+ a[-1,0]*b a[-1,1]*b ... a[-1,-1]*b
+
+ Parameters
+ ----------
+ a : array, shape (M, N)
+ b : array, shape (P, Q)
+
+ Returns
+ -------
+ A : array, shape (M*P, N*Q)
+ Kronecker product of a and b
+
+ Examples
+ --------
+ >>> from scipy import kron, array
+ >>> kron(array([[1,2],[3,4]]), array([[1,1,1]]))
+ array([[1, 1, 1, 2, 2, 2],
+ [3, 3, 3, 4, 4, 4]])
+
"""
if not a.flags['CONTIGUOUS']:
a = reshape(a, a.shape)
Modified: trunk/scipy/linalg/blas.py
===================================================================
--- trunk/scipy/linalg/blas.py 2008-02-20 05:39:03 UTC (rev 3951)
+++ trunk/scipy/linalg/blas.py 2008-02-21 01:57:37 UTC (rev 3952)
@@ -22,12 +22,15 @@
_inv_type_conv = {'s':'f','d':'d','c':'F','z':'D'}
def has_column_major_storage(arr):
+ """Is array stored in column-major format"""
return arr.flags['FORTRAN']
def get_blas_funcs(names,arrays=(),debug=0):
"""Return available BLAS function objects with names.
arrays are used to determine the optimal prefix of
- BLAS routines."""
+ BLAS routines.
+
+ """
ordering = []
for i in range(len(arrays)):
t = arrays[i].dtype.char
Modified: trunk/scipy/linalg/decomp.py
===================================================================
--- trunk/scipy/linalg/decomp.py 2008-02-20 05:39:03 UTC (rev 3951)
+++ trunk/scipy/linalg/decomp.py 2008-02-21 01:57:37 UTC (rev 3952)
@@ -98,32 +98,54 @@
return w, vr
def eig(a,b=None, left=False, right=True, overwrite_a=False, overwrite_b=False):
- """ Solve ordinary and generalized eigenvalue problem
- of a square matrix.
+ """Solve an ordinary or generalized eigenvalue problem of a square matrix.
- Inputs:
+ Find eigenvalues w and right or left eigenvectors of a general matrix::
+
+ a vr[:,i] = w[i] b vr[:,i]
+ a.H vl[:,i] = w[i].conj() b.H vl[:,i]
+
+ where .H is the Hermitean conjugation.
+
+ Parameters
+ ----------
+ a : array, shape (M, M)
+ A complex or real matrix whose eigenvalues and eigenvectors
+ will be computed.
+ b : array, shape (M, M)
+ Right-hand side matrix in a generalized eigenvalue problem.
+ If omitted, identity matrix is assumed.
+ left : boolean
+ Whether to calculate and return left eigenvectors
+ right : boolean
+ Whether to calculate and return right eigenvectors
+
+ overwrite_a : boolean
+ Whether to overwrite data in a (may improve performance)
+ overwrite_b : boolean
+ Whether to overwrite data in b (may improve performance)
+
+ Returns
+ -------
+ w : double or complex array, shape (M,)
+ The eigenvalues, each repeated according to its multiplicity.
- a -- An N x N matrix.
- b -- An N x N matrix [default is identity(N)].
- left -- Return left eigenvectors [disabled].
- right -- Return right eigenvectors [enabled].
- overwrite_a, overwrite_b -- save space by overwriting the a and/or
- b matrices (both False by default)
+ (if left == True)
+ vl : double or complex array, shape (M, M)
+ The normalized left eigenvector corresponding to the eigenvalue w[i]
+ is the column v[:,i].
+
+ (if right == True)
+ vr : double or complex array, shape (M, M)
+ The normalized right eigenvector corresponding to the eigenvalue w[i]
+ is the column vr[:,i].
+
+ Raises LinAlgError if eigenvalue computation does not converge
- Outputs:
-
- w -- eigenvalues [left==right==False].
- w,vr -- w and right eigenvectors [left==False,right=True].
- w,vl -- w and left eigenvectors [left==True,right==False].
- w,vl,vr -- [left==right==True].
-
- Definitions:
-
- a * vr[:,i] = w[i] * b * vr[:,i]
-
- a^H * vl[:,i] = conjugate(w[i]) * b^H * vl[:,i]
-
- where a^H denotes transpose(conjugate(a)).
+ See Also
+ --------
+ eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
+
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
@@ -182,28 +204,44 @@
return w, vr
def eigh(a, lower=True, eigvals_only=False, overwrite_a=False):
- """ Solve real symmetric or complex hermitian eigenvalue problem.
+ """Solve the eigenvalue problem for a Hermitian or real symmetric matrix.
- Inputs:
+ Find eigenvalues w and optionally right eigenvectors v of a::
+
+ a v[:,i] = w[i] v[:,i]
+ v.H v = identity
+
+ Parameters
+ ----------
+ a : array, shape (M, M)
+ A complex Hermitian or symmetric real matrix whose eigenvalues
+ and eigenvectors will be computed.
+ lower : boolean
+ Whether the pertinent array data is taken from the lower or upper
+ triangle of a. (Default: lower)
+ eigvals_only : boolean
+ Whether to calculate only eigenvalues and no eigenvectors.
+ (Default: both are calculated)
+ overwrite_a : boolean
+ Whether data in a is overwritten (may improve performance).
+
+ Returns
+ -------
+ w : double array, shape (M,)
+ The eigenvalues, in ascending order, each repeated according to its
+ multiplicity.
+
+ (if eigvals_only == False)
+ v : double or complex double array, shape (M, M)
+ The normalized eigenvector corresponding to the eigenvalue w[i] is
+ the column v[:,i].
+
+ Raises LinAlgError if eigenvalue computation does not converge
- a -- A hermitian N x N matrix.
- lower -- values in a are read from lower triangle
- [True: UPLO='L' (default) / False: UPLO='U']
- eigvals_only -- don't compute eigenvectors.
- overwrite_a -- content of a may be destroyed
-
- Outputs:
-
- For eigvals_only == False (the default),
- w,v -- w: eigenvalues, v: eigenvectors
- For eigvals_only == True,
- w -- eigenvalues
-
- Definitions:
-
- a * v[:,i] = w[i] * vr[:,i]
- v.H * v = identity
-
+ See Also
+ --------
+ eig : eigenvalues and right eigenvectors for non-symmetric arrays
+
"""
if eigvals_only or overwrite_a:
a1 = asarray_chkfinite(a)
@@ -258,43 +296,74 @@
def eig_banded(a_band, lower=0, eigvals_only=0, overwrite_a_band=0,
select='a', select_range=None, max_ev = 0):
- """ Solve real symmetric or complex hermetian band matrix problem.
+ """Solve real symmetric or complex hermetian band matrix eigenvalue problem.
- Inputs:
+ Find eigenvalues w and optionally right eigenvectors v of a::
+
+ a v[:,i] = w[i] v[:,i]
+ v.H v = identity
+
+ The matrix a is stored in ab either in lower diagonal or upper
+ diagonal ordered form:
- a_band -- A hermitian N x M matrix in 'packed storage'.
- Packed storage looks like this: ('upper form')
- [ ... (more off-diagonals) ...,
- [ * * a13 a24 a35 a46 ... a(n-2)(n)],
- [ * a12 a23 a34 a45 a56 ... a(n-1)(n)],
- [ a11 a22 a33 a44 a55 a66 ... a(n)(n) ]]
- The cells denoted with * may contain anything.
- lower -- a is in lower packed storage
- (default: upper packed form)
- eigvals_only -- if True, don't compute eigenvectors.
- overwrite_a_band -- content of a may be destroyed
- select -- 'a', 'all', 0 : return all eigenvalues/vectors
- 'v', 'value', 1 : eigenvalues in the interval (min, max]
- will be found
- 'i', 'index', 2 : eigenvalues with index [min...max]
- will be found
- select_range -- select == 'v': eigenvalue limits as tuple (min, max)
- select == 'i': index limits as tuple (min, max)
- select == 'a': meaningless
- max_ev -- select == 'v': set to max. number of eigenvalues that is
- expected. In doubt, leave untouched.
- select == 'i', 'a': meaningless
+ ab[u + i - j, j] == a[i,j] (if upper form; i <= j)
+ ab[ i - j, j] == a[i,j] (if lower form; i >= j)
- Outputs:
+ Example of ab (shape of a is (6,6), u=2)::
- w,v -- w: eigenvalues, v: eigenvectors [for eigvals_only == False]
- w -- eigenvalues [for eigvals_only == True].
+ upper form:
+ * * a02 a13 a24 a35
+ * a01 a12 a23 a34 a45
+ a00 a11 a22 a33 a44 a55
+
+ lower form:
+ a00 a11 a22 a33 a44 a55
+ a10 a21 a32 a43 a54 *
+ a20 a31 a42 a53 * *
- Definitions:
+ Cells marked with * are not used.
- a_full * v[:,i] = w[i] * v[:,i] (with full matrix corresponding to a)
- v.H * v = identity
+ Parameters
+ ----------
+ a_band : array, shape (M, u+1)
+ Banded matrix whose eigenvalues to calculate
+ lower : boolean
+ Is the matrix in the lower form. (Default is upper form)
+ eigvals_only : boolean
+ Compute only the eigenvalues and no eigenvectors.
+ (Default: calculate also eigenvectors)
+ overwrite_a_band:
+ Discard data in a_band (may enhance performance)
+ select: {'a', 'v', 'i'}
+ Which eigenvalues to calculate
+ ====== ========================================
+ select calculated
+ ====== ========================================
+ 'a' All eigenvalues
+ 'v' Eigenvalues in the interval (min, max]
+ 'i' Eigenvalues with indices min <= i <= max
+ ====== ========================================
+ select_range : (min, max)
+ Range of selected eigenvalues
+ max_ev : integer
+ For select=='v', maximum number of eigenvalues expected.
+ For other values of select, has no meaning.
+
+ In doubt, leave this parameter untouched.
+
+ Returns
+ -------
+ w : array, shape (M,)
+ The eigenvalues, in ascending order, each repeated according to its
+ multiplicity.
+
+ v : double or complex double array, shape (M, M)
+ The normalized eigenvector corresponding to the eigenvalue w[i] is
+ the column v[:,i].
+
+ Raises LinAlgError if eigenvalue computation does not converge
+
"""
if eigvals_only or overwrite_a_band:
a1 = asarray_chkfinite(a_band)
@@ -374,32 +443,180 @@
return w, v
def eigvals(a,b=None,overwrite_a=0):
- """Return eigenvalues of square matrix."""
+ """Compute eigenvalues from an ordinary or generalized eigenvalue problem.
+
+ Find eigenvalues of a general matrix::
+
+ a vr[:,i] = w[i] b vr[:,i]
+
+ Parameters
+ ----------
+ a : array, shape (M, M)
+ A complex or real matrix whose eigenvalues and eigenvectors
+ will be computed.
+ b : array, shape (M, M)
+ Right-hand side matrix in a generalized eigenvalue problem.
+ If omitted, identity matrix is assumed.
+ overwrite_a : boolean
+ Whether to overwrite data in a (may improve performance)
+
+ Returns
+ -------
+ w : double or complex array, shape (M,)
+ The eigenvalues, each repeated according to its multiplicity,
+ but not in any specific order.
+
+ Raises LinAlgError if eigenvalue computation does not converge
+
+ See Also
+ --------
+ eigvalsh : eigenvalues of symmetric or Hemitiean arrays
+ eig : eigenvalues and right eigenvectors of general arrays
+ eigh : eigenvalues and eigenvectors of symmetric/Hermitean arrays.
+
+ """
return eig(a,b=b,left=0,right=0,overwrite_a=overwrite_a)
def eigvalsh(a,lower=1,overwrite_a=0):
- """Return eigenvalues of hermitean or real symmetric matrix."""
+ """Solve the eigenvalue problem for a Hermitian or real symmetric matrix.
+
+ Find eigenvalues w of a::
+
+ a v[:,i] = w[i] v[:,i]
+ v.H v = identity
+
+ Parameters
+ ----------
+ a : array, shape (M, M)
+ A complex Hermitian or symmetric real matrix whose eigenvalues
+ and eigenvectors will be computed.
+ lower : boolean
+ Whether the pertinent array data is taken from the lower or upper
+ triangle of a. (Default: lower)
+ overwrite_a : boolean
+ Whether data in a is overwritten (may improve performance).
+
+ Returns
+ -------
+ w : double array, shape (M,)
+ The eigenvalues, in ascending order, each repeated according to its
+ multiplicity.
+
+ Raises LinAlgError if eigenvalue computation does not converge
+
+ See Also
+ --------
+ eigvals : eigenvalues of general arrays
+ eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
+ eig : eigenvalues and right eigenvectors for non-symmetric arrays
+
+ """
return eigh(a,lower=lower,eigvals_only=1,overwrite_a=overwrite_a)
def eigvals_banded(a_band,lower=0,overwrite_a_band=0,
select='a', select_range=None):
- """Return eigenvalues of hermitean or real symmetric matrix."""
+ """Solve real symmetric or complex hermetian band matrix eigenvalue problem.
+
+ Find eigenvalues w of a::
+
+ a v[:,i] = w[i] v[:,i]
+ v.H v = identity
+
+ The matrix a is stored in ab either in lower diagonal or upper
+ diagonal ordered form:
+
+ ab[u + i - j, j] == a[i,j] (if upper form; i <= j)
+ ab[ i - j, j] == a[i,j] (if lower form; i >= j)
+
+ Example of ab (shape of a is (6,6), u=2)::
+
+ upper form:
+ * * a02 a13 a24 a35
+ * a01 a12 a23 a34 a45
+ a00 a11 a22 a33 a44 a55
+
+ lower form:
+ a00 a11 a22 a33 a44 a55
+ a10 a21 a32 a43 a54 *
+ a20 a31 a42 a53 * *
+
+ Cells marked with * are not used.
+
+ Parameters
+ ----------
+ a_band : array, shape (M, u+1)
+ Banded matrix whose eigenvalues to calculate
+ lower : boolean
+ Is the matrix in the lower form. (Default is upper form)
+ overwrite_a_band:
+ Discard data in a_band (may enhance performance)
+ select: {'a', 'v', 'i'}
+ Which eigenvalues to calculate
+
+ ====== ========================================
+ select calculated
+ ====== ========================================
+ 'a' All eigenvalues
+ 'v' Eigenvalues in the interval (min, max]
+ 'i' Eigenvalues with indices min <= i <= max
+ ====== ========================================
+ select_range : (min, max)
+ Range of selected eigenvalues
+
+ Returns
+ -------
+ w : array, shape (M,)
+ The eigenvalues, in ascending order, each repeated according to its
+ multiplicity.
+
+ Raises LinAlgError if eigenvalue computation does not converge
+
+ See Also
+ --------
+ eig_banded : eigenvalues and right eigenvectors for symmetric/Hermitian band matrices
+ eigvals : eigenvalues of general arrays
+ eigh : eigenvalues and right eigenvectors for symmetric/Hermitian arrays
+ eig : eigenvalues and right eigenvectors for non-symmetric arrays
+
+ """
return eig_banded(a_band,lower=lower,eigvals_only=1,
overwrite_a_band=overwrite_a_band, select=select,
select_range=select_range)
def lu_factor(a, overwrite_a=0):
- """Return raw LU decomposition of a matrix and pivots, for use in solving
- a system of linear equations.
+ """Compute pivoted LU decomposition of a matrix.
+
+ The decomposition is::
- Inputs:
+ A = P L U
- a --- an NxN matrix
+ where P is a permutation matrix, L lower triangular with unit
+ diagonal elements, and U upper triangular.
+
+ Parameters
+ ----------
+ a : array, shape (M, M)
+ Matrix to decompose
+ overwrite_a : boolean
+ Whether to overwrite data in A (may increase performance)
- Outputs:
+ Returns
+ -------
+ lu : array, shape (N, N)
+ Matrix containing U in its upper triangle, and L in its lower triangle.
+ The unit diagonal elements of L are not stored.
+ piv : array, shape (N,)
+ Pivot indices representing the permutation matrix P:
+ row i of matrix was interchanged with row piv[i].
- lu --- the lu factorization matrix
- piv --- an array of pivots
+ See also
+ --------
+ lu_solve : solve an equation system using the LU factorization of a matrix
+
+ Notes
+ -----
+ This is a wrapper to the *GETRF routines from LAPACK.
+
"""
a1 = asarray(a)
if len(a1.shape) != 2 or (a1.shape[0] != a1.shape[1]):
@@ -414,8 +631,24 @@
return lu, piv
def lu_solve(a_lu_pivots,b):
- """Solve a previously factored system. First input is a tuple (lu, pivots)
- which is the output to lu_factor. Second input is the right hand side.
+ """Solve an equation system, a x = b, given the LU factorization of a
+
+ Parameters
+ ----------
+ (lu, piv)
+ Factorization of the coefficient matrix a, as given by lu_factor
+ b : array
+ Right-hand side
+
+ Returns
+ -------
+ x : array
+ Solution to the system
+
+ See also
+ --------
+ lu_factor : LU factorize a matrix
+
"""
a_lu, pivots = a_lu_pivots
a_lu = asarray_chkfinite(a_lu)
@@ -435,28 +668,46 @@
def lu(a,permute_l=0,overwrite_a=0):
- """Return LU decompostion of a matrix.
+ """Compute pivoted LU decompostion of a matrix.
- Inputs:
+ The decomposition is::
- a -- An M x N matrix.
- permute_l -- Perform matrix multiplication p * l [disabled].
+ A = P L U
- Outputs:
+ where P is a permutation matrix, L lower triangular with unit
+ diagonal elements, and U upper triangular.
+
+ Parameters
+ ----------
+ a : array, shape (M, N)
+ Array to decompose
+ permute_l : boolean
+ Perform the multiplication P*L (Default: do not permute)
+ overwrite_a : boolean
+ Whether to overwrite data in a (may improve performance)
- p,l,u -- LU decomposition matrices of a [permute_l=0]
- pl,u -- LU decomposition matrices of a [permute_l=1]
-
- Definitions:
-
- a = p * l * u [permute_l=0]
- a = pl * u [permute_l=1]
-
- p - An M x M permutation matrix
- l - An M x K lower triangular or trapezoidal matrix
- with unit-diagonal
- u - An K x N upper triangular or trapezoidal matrix
- K = min(M,N)
+ Returns
+ -------
+ (If permute_l == False)
+ p : array, shape (M, M)
+ Permutation matrix
+ l : array, shape (M, K)
+ Lower triangular or trapezoidal matrix with unit diagonal.
+ K = min(M, N)
+ u : array, shape (K, N)
+ Upper triangular or trapezoidal matrix
+
+ (If permute_l == True)
+ pl : array, shape (M, K)
+ Permuted L matrix.
+ K = min(M, N)
+ u : array, shape (K, N)
+ Upper triangular or trapezoidal matrix
+
+ Notes
+ -----
+ This is a LU factorization routine written for Scipy.
+
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2:
@@ -471,29 +722,60 @@
return p,l,u
def svd(a,full_matrices=1,compute_uv=1,overwrite_a=0):
- """Compute singular value decomposition (SVD) of matrix a.
+ """Singular Value Decomposition.
- Description:
+ Factorizes the matrix a into two unitary matrices U and Vh and
+ an 1d-array s of singular values (real, non-negative) such that
+ a == U S Vh if S is an suitably shaped matrix of zeros whose
+ main diagonal is s.
+
+ Parameters
+ ----------
+ a : array, shape (M, N)
+ Matrix to decompose
+ full_matrices : boolean
+ If true, U, Vh are shaped (M,M), (N,N)
+ If false, the shapes are (M,K), (K,N) where K = min(M,N)
+ compute_uv : boolean
+ Whether to compute also U, Vh in addition to s (Default: true)
+ overwrite_a : boolean
+ Whether data in a is overwritten (may improve performance)
+
+ Returns
+ -------
+ U: array, shape (M,M) or (M,K) depending on full_matrices
+ s: array, shape (K,)
+ The singular values, sorted so that s[i] >= s[i+1]. K = min(M, N)
+ Vh: array, shape (N,N) or (K,N) depending on full_matrices
- Singular value decomposition of a matrix a is
- a = u * sigma * v^H,
- where v^H denotes conjugate(transpose(v)), u,v are unitary
- matrices, sigma is zero matrix with a main diagonal containing
- real non-negative singular values of the matrix a.
+ For compute_uv = False, only s is returned.
- Inputs:
+ Raises LinAlgError if SVD computation does not converge
- a -- An M x N matrix.
- compute_uv -- If zero, then only the vector of singular values
- is returned.
+ Examples
+ --------
+ >>> from scipy import random, linalg, allclose, dot
+ >>> a = random.randn(9, 6) + 1j*random.randn(9, 6)
+ >>> U, s, Vh = linalg.svd(a)
+ >>> U.shape, Vh.shape, s.shape
+ ((9, 9), (6, 6), (6,))
+
+ >>> U, s, Vh = linalg.svd(a, full_matrices=False)
+ >>> U.shape, Vh.shape, s.shape
+ ((9, 6), (6, 6), (6,))
+ >>> S = linalg.diagsvd(s, 6, 6)
+ >>> allclose(a, dot(U, dot(S, Vh)))
+ True
+
+ >>> s2 = linalg.svd(a, compute_uv=False)
+ >>> allclose(s, s2)
+ True
- Outputs:
-
- u -- An M x M unitary matrix [compute_uv=1].
- s -- An min(M,N) vector of singular values in descending order,
- sigma = diagsvd(s).
- vh -- An N x N unitary matrix [compute_uv=1], vh = v^H.
-
+ See also
+ --------
+ svdvals : return singular values of a matrix
+ diagsvd : return the Sigma matrix, given the vector s
+
"""
# A hack until full_matrices == 0 support is fixed here.
if full_matrices == 0:
@@ -520,11 +802,47 @@
return s
def svdvals(a,overwrite_a=0):
- """Return singular values of a matrix."""
+ """Compute singular values of a matrix.
+
+ Parameters
+ ----------
+ a : array, shape (M, N)
+ Matrix to decompose
+ overwrite_a : boolean
+ Whether data in a is overwritten (may improve performance)
+
+ Returns
+ -------
+ s: array, shape (K,)
+ The singular values, sorted so that s[i] >= s[i+1]. K = min(M, N)
+
+ Raises LinAlgError if SVD computation does not converge
+
+ See also
+ --------
+ svd : return the full singular value decomposition of a matrix
+ diagsvd : return the Sigma matrix, given the vector s
+
+ """
return svd(a,compute_uv=0,overwrite_a=overwrite_a)
def diagsvd(s,M,N):
- """Return sigma from singular values and original size M,N."""
+ """Construct the sigma matrix in SVD from singular values and size M,N.
+
+ Parameters
+ ----------
+ s : array, shape (M,) or (N,)
+ Singular values
+ M : integer
+ N : integer
+ Size of the matrix whose singular values are s
+
+ Returns
+ -------
+ S : array, shape (M, N)
+ The S-matrix in the singular value decomposition
+
+ """
part = diag(s)
typ = part.dtype.char
MorN = len(s)
@@ -536,14 +854,40 @@
raise ValueError, "Length of s must be M or N."
def cholesky(a,lower=0,overwrite_a=0):
- """Compute Cholesky decomposition of matrix.
-
- Description:
-
- For a hermitian positive-definite matrix a return the
- upper-triangular (or lower-triangular if lower==1) matrix,
- u such that u^H * u = a (or l * l^H = a).
-
+ """Compute the Cholesky decomposition of a matrix.
+
+ Returns the Cholesky decomposition, :lm:`A = L L^*` or :lm:`A = U^* U`
+ of a Hermitian positive-definite matrix :lm:`A`.
+
+ Parameters
+ ----------
+ a : array, shape (M, M)
+ Matrix to be decomposed
+ lower : boolean
+ Whether to compute the upper or lower triangular Cholesky factorization
+ (Default: upper-triangular)
+ overwrite_a : boolean
+ Whether to overwrite data in a (may improve performance)
+
+ Returns
+ -------
+ B : array, shape (M, M)
+ Upper- or lower-triangular Cholesky factor of A
+
+ Raises LinAlgError if decomposition fails
+
+ Examples
+ --------
+ >>> from scipy import array, linalg, dot
+ >>> a = array([[1,-2j],[2j,5]])
+ >>> L = linalg.cholesky(a, lower=True)
+ >>> L
+ array([[ 1.+0.j, 0.+0.j],
+ [ 0.+2.j, 1.+0.j]])
+ >>> dot(L, L.T.conj())
+ array([[ 1.+0.j, 0.-2.j],
+ [ 0.+2.j, 5.+0.j]])
+
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
@@ -557,8 +901,32 @@
return c
def cho_factor(a, lower=0, overwrite_a=0):
- """ Compute Cholesky decomposition of matrix and return an object
- to be used for solving a linear system using cho_solve.
+ """Compute the Cholesky decomposition of a matrix, to use in cho_solve
+
+ Returns the Cholesky decomposition, :lm:`A = L L^*` or :lm:`A = U^* U`
+ of a Hermitian positive-definite matrix :lm:`A`.
+
+ The return value can be directly used as the first parameter to cho_solve.
+
+ Parameters
+ ----------
+ a : array, shape (M, M)
+ Matrix to be decomposed
+ lower : boolean
+ Whether to compute the upper or lower triangular Cholesky factorization
+ (Default: upper-triangular)
+ overwrite_a : boolean
+ Whether to overwrite data in a (may improve performance)
+
+ Returns
+ -------
+ c : array, shape (M, M)
+ Upper- or lower-triangular Cholesky factor of A
+ lower : array, shape (M, M)
+ Flag indicating whether the factor is lower or upper triangular
+
+ Raises LinAlgError if decomposition fails
+
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
@@ -573,8 +941,29 @@
def cho_solve(clow, b):
"""Solve a previously factored symmetric system of equations.
+
+ The equation system is
+
+ A x = b, A = U^H U = L L^H
+
+ and A is real symmetric or complex Hermitian.
+
+ Parameters
+ ----------
+ clow : tuple (c, lower)
+ Cholesky factor and a flag indicating whether it is lower triangular.
+ The return value from cho_factor can be used.
+ b : array
+ Right-hand side of the equation system
+
First input is a tuple (LorU, lower) which is the output to cho_factor.
Second input is the right-hand side.
+
+ Returns
+ -------
+ x : array
+ Solution to the equation system
+
"""
c, lower = clow
c = asarray_chkfinite(c)
@@ -591,29 +980,60 @@
return b
def qr(a,overwrite_a=0,lwork=None,econ=False,mode='qr'):
- """QR decomposition of an M x N matrix a.
+ """Compute QR decomposition of a matrix.
- Description:
+ Calculate the decomposition :lm:`A = Q R` where Q is unitary/orthogonal
+ and R upper triangular.
- Find a unitary (orthogonal) matrix, q, and an upper-triangular
- matrix r such that q * r = a
+ Parameters
+ ----------
+ a : array, shape (M, N)
+ Matrix to be decomposed
+ overwrite_a : boolean
+ Whether data in a is overwritten (may improve performance)
+ lwork : integer
+ Work array size, lwork >= a.shape[1]. If None or -1, an optimal size
+ is computed.
+ econ : boolean
+ Whether to compute the economy-size QR decomposition, making shapes
+ of Q and R (M, K) and (K, N) instead of (M,M) and (M,N). K=min(M,N)
+ mode : {'qr', 'r'}
+ Determines what information is to be returned: either both Q and R
+ or only R.
+
+ Returns
+ -------
+ (if mode == 'qr')
+ Q : double or complex array, shape (M, M) or (M, K) for econ==True
- Inputs:
+ (for any mode)
+ R : double or complex array, shape (M, N) or (K, N) for econ==True
+ Size K = min(M, N)
- a -- the matrix
- overwrite_a=0 -- if non-zero then discard the contents of a,
- i.e. a is used as a work array if possible.
+ Raises LinAlgError if decomposition fails
- lwork=None -- >= shape(a)[1]. If None (or -1) compute optimal
- work array size.
- econ=False -- computes the skinny or economy-size QR decomposition
- only useful when M>N
- mode='qr' -- if 'qr' then return both q and r; if 'r' then just return r
+ Notes
+ -----
+ This is an interface to the LAPACK routines dgeqrf, zgeqrf,
+ dorgqr, and zungqr.
- Outputs:
- q,r - if mode=='qr'
- r - if mode=='r'
+ Examples
+ --------
+ >>> from scipy import random, linalg, dot
+ >>> a = random.randn(9, 6)
+ >>> q, r = linalg.qr(a)
+ >>> allclose(a, dot(q, r))
+ True
+ >>> q.shape, r.shape
+ ((9, 9), (9, 6))
+ >>> r2 = linalg.qr(a, mode='r')
+ >>> allclose(r, r2)
+
+ >>> q3, r3 = linalg.qr(a, econ=True)
+ >>> q3.shape, r3.shape
+ ((9, 6), (6, 6))
+
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2:
@@ -674,26 +1094,29 @@
def qr_old(a,overwrite_a=0,lwork=None):
- """QR decomposition of an M x N matrix a.
+ """Compute QR decomposition of a matrix.
- Description:
+ Calculate the decomposition :lm:`A = Q R` where Q is unitary/orthogonal
+ and R upper triangular.
- Find a unitary (orthogonal) matrix, q, and an upper-triangular
- matrix r such that q * r = a
+ Parameters
+ ----------
+ a : array, shape (M, N)
+ Matrix to be decomposed
+ overwrite_a : boolean
+ Whether data in a is overwritten (may improve performance)
+ lwork : integer
+ Work array size, lwork >= a.shape[1]. If None or -1, an optimal size
+ is computed.
+
+ Returns
+ -------
+ Q : double or complex array, shape (M, M)
+ R : double or complex array, shape (M, N)
+ Size K = min(M, N)
- Inputs:
-
- a -- the matrix
- overwrite_a=0 -- if non-zero then discard the contents of a,
- i.e. a is used as a work array if possible.
-
- lwork=None -- >= shape(a)[1]. If None (or -1) compute optimal
- work array size.
-
- Outputs:
-
- q, r -- matrices such that q * r = a
-
+ Raises LinAlgError if decomposition fails
+
"""
a1 = asarray_chkfinite(a)
if len(a1.shape) != 2:
@@ -725,26 +1148,30 @@
def rq(a,overwrite_a=0,lwork=None):
- """RQ decomposition of an M x N matrix a.
+ """Compute RQ decomposition of a square real matrix.
- Description:
+ Calculate the decomposition :lm:`A = R Q` where Q is unitary/orthogonal
+ and R upper triangular.
- Find an upper-triangular matrix r and a unitary (orthogonal)
- matrix q such that r * q = a
-
- Inputs:
-
- a -- the matrix
- overwrite_a=0 -- if non-zero then discard the contents of a,
- i.e. a is used as a work array if possible.
-
- lwork=None -- >= shape(a)[1]. If None (or -1) compute optimal
- work array size.
-
- Outputs:
-
- r, q -- matrices such that r * q = a
-
+ Parameters
+ ----------
+ a : array, shape (M, M)
+ Square real matrix to be decomposed
+ overwrite_a : boolean
+ Whether data in a is overwritten (may improve performance)
+ lwork : integer
+ Work array size, lwork >= a.shape[1]. If None or -1, an optimal size
+ is computed.
+ econ : boolean
+
+ Returns
+ -------
+ R : double array, shape (M, N) or (K, N) for econ==True
+ Size K = min(M, N)
+ Q : double or complex array, shape (M, M) or (M, K) for econ==True
+
+ Raises LinAlgError if decomposition fails
+
"""
# TODO: implement support for non-square and complex arrays
a1 = asarray_chkfinite(a)
@@ -783,13 +1210,40 @@
_double_precision = ['i','l','d']
def schur(a,output='real',lwork=None,overwrite_a=0):
- """Compute Schur decomposition of matrix a.
+ """Compute Schur decomposition of a matrix.
+
+ The Schur decomposition is
+
+ A = Z T Z^H
+
+ where Z is unitary and T is either upper-triangular, or for real
+ Schur decomposition (output='real'), quasi-upper triangular. In
+ the quasi-triangular form, 2x2 blocks describing complex-valued
+ eigenvalue pairs may extrude from the diagonal.
+
+ Parameters
+ ----------
+ a : array, shape (M, M)
+ Matrix to decompose
+ output : {'real', 'complex'}
+ Construct the real or complex Schur decomposition (for real matrices).
+ lwork : integer
+ Work array size. If None or -1, it is automatically computed.
+ overwrite_a : boolean
+ Whether to overwrite data in a (may improve performance)
- Description:
+ Returns
+ -------
+ T : array, shape (M, M)
+ Schur form of A. It is real-valued for the real Schur decomposition.
+ Z : array, shape (M, M)
+ An unitary Schur transformation matrix for A.
+ It is real-valued for the real Schur decomposition.
- Return T, Z such that a = Z * T * (Z**H) where Z is a
- unitary matrix and T is either upper-triangular or quasi-upper
- triangular for output='real'
+ See also
+ --------
+ rsf2csf : Convert real Schur form to complex Schur form
+
"""
if not output in ['real','complex','r','c']:
raise ValueError, "argument must be 'real', or 'complex'"
@@ -850,16 +1304,29 @@
raise LinAlgError, 'Array must be square'
def rsf2csf(T, Z):
- """Convert real schur form to complex schur form.
-
- Description:
-
- If A is a real-valued matrix, then the real schur form is
- quasi-upper triangular. 2x2 blocks extrude from the main-diagonal
- corresponding to any complex-valued eigenvalues.
-
- This function converts this real schur form to a complex schur form
- which is upper triangular.
+ """Convert real Schur form to complex Schur form.
+
+ Convert a quasi-diagonal real-valued Schur form to the upper triangular
+ complex-valued Schur form.
+
+ Parameters
+ ----------
+ T : array, shape (M, M)
+ Real Schur form of the original matrix
+ Z : array, shape (M, M)
+ Schur transformation matrix
+
+ Returns
+ -------
+ T : array, shape (M, M)
+ Complex Schur form of the original matrix
+ Z : array, shape (M, M)
+ Schur transformation matrix corresponding to the complex form
+
+ See also
+ --------
+ schur : Schur decompose a matrix
+
"""
Z,T = map(asarray_chkfinite, (Z,T))
if len(Z.shape) !=2 or Z.shape[0] != Z.shape[1]:
@@ -898,7 +1365,23 @@
# Orthonormal decomposition
def orth(A):
- """Return an orthonormal basis for the range of A using svd"""
+ """Construct an orthonormal basis for the range of A using SVD
+
+ Parameters
+ ----------
+ A : array, shape (M, N)
+
+ Returns
+ -------
+ Q : array, shape (M, K)
+ Orthonormal basis for the range of A.
+ K = effective rank of A, as determined by automatic cutoff
+
+ See also
+ --------
+ svd : Singular value decomposition of a matrix
+
+ """
u,s,vh = svd(A)
M,N = A.shape
tol = max(M,N)*numpy.amax(s)*eps
@@ -907,21 +1390,33 @@
return Q
def hessenberg(a,calc_q=0,overwrite_a=0):
- """ Compute Hessenberg form of a matrix.
+ """Compute Hessenberg form of a matrix.
+
+ The Hessenberg decomposition is
+
+ A = Q H Q^H
+
+ where Q is unitary/orthogonal and H has only zero elements below the first
+ subdiagonal.
+
+ Parameters
+ ----------
+ a : array, shape (M,M)
+ Matrix to bring into Hessenberg form
+ calc_q : boolean
+ Whether to compute the transformation matrix
+ overwrite_a : boolean
+ Whether to ovewrite data in a (may improve performance)
- Inputs:
+ Returns
+ -------
+ H : array, shape (M,M)
+ Hessenberg form of A
- a -- the matrix
- calc_q -- if non-zero then calculate unitary similarity
- transformation matrix q.
- overwrite_a=0 -- if non-zero then discard the contents of a,
- i.e. a is used as a work array if possible.
+ (If calc_q == True)
+ Q : array, shape (M,M)
+ Unitary/orthogonal similarity transformation matrix s.t. A = Q H Q^H
- Outputs:
-
- h -- Hessenberg form of a [calc_q=0]
- h, q -- matrices such that a = q * h * q^T [calc_q=1]
-
"""
a1 = asarray(a)
if len(a1.shape) != 2 or (a1.shape[0] != a1.shape[1]):
Modified: trunk/scipy/linalg/matfuncs.py
===================================================================
--- trunk/scipy/linalg/matfuncs.py 2008-02-20 05:39:03 UTC (rev 3951)
+++ trunk/scipy/linalg/matfuncs.py 2008-02-21 01:57:37 UTC (rev 3952)
@@ -20,7 +20,20 @@
feps = sb.finfo(single).eps
def expm(A,q=7):
- """Compute the matrix exponential using Pade approximation of order q.
+ """Compute the matrix exponential using Pade approximation.
+
+ Parameters
+ ----------
+ A : array, shape(M,M)
+ Matrix to be exponentiated
+ q : integer
+ Order of the Pade approximation
+
+ Returns
+ -------
+ expA : array, shape(M,M)
+ Matrix exponential of A
+
"""
A = asarray(A)
ss = True
@@ -61,6 +74,17 @@
def expm2(A):
"""Compute the matrix exponential using eigenvalue decomposition.
+
+ Parameters
+ ----------
+ A : array, shape(M,M)
+ Matrix to be exponentiated
+
+ Returns
+ -------
+ expA : array, shape(M,M)
+ Matrix exponential of A
+
"""
A = asarray(A)
t = A.dtype.char
@@ -72,7 +96,20 @@
return dot(dot(vr,diag(exp(s))),vri).astype(t)
def expm3(A,q=20):
- """Compute the matrix exponential using a Taylor series.of order q.
+ """Compute the matrix exponential using Taylor series.
+
+ Parameters
+ ----------
+ A : array, shape(M,M)
+ Matrix to be exponentiated
+ q : integer
+ Order of the Taylor series
+
+ Returns
+ -------
+ expA : array, shape(M,M)
+ Matrix exponential of A
+
"""
A = asarray(A)
t = A.dtype.char
@@ -91,6 +128,16 @@
_array_precision = {'i': 1, 'l': 1, 'f': 0, 'd': 1, 'F': 0, 'D': 1}
def toreal(arr,tol=None):
"""Return as real array if imaginary part is small.
+
+ Parameters
+ ----------
+ arr : array
+ tol : float
+ Absolute tolerance
+
+ Returns
+ -------
+ arr : double or complex array
"""
if tol is None:
tol = {0:feps*1e3, 1:eps*1e6}[_array_precision[arr.dtype.char]]
@@ -100,7 +147,19 @@
return arr
def cosm(A):
- """matrix cosine.
+ """Compute the matrix cosine.
+
+ This routine uses expm to compute the matrix exponentials.
+
+ Parameters
+ ----------
+ A : array, shape(M,M)
+
+ Returns
+ -------
+ cosA : array, shape(M,M)
+ Matrix cosine of A
+
"""
A = asarray(A)
if A.dtype.char not in ['F','D','G']:
@@ -110,7 +169,19 @@
def sinm(A):
- """matrix sine.
+ """Compute the matrix sine.
+
+ This routine uses expm to compute the matrix exponentials.
+
+ Parameters
+ ----------
+ A : array, shape(M,M)
+
+ Returns
+ -------
+ sinA : array, shape(M,M)
+ Matrix cosine of A
+
"""
A = asarray(A)
if A.dtype.char not in ['F','D','G']:
@@ -119,7 +190,19 @@
return -0.5j*(expm(1j*A) - expm(-1j*A))
def tanm(A):
- """matrix tangent.
+ """Compute the matrix tangent.
+
+ This routine uses expm to compute the matrix exponentials.
+
+ Parameters
+ ----------
+ A : array, shape(M,M)
+
+ Returns
+ -------
+ tanA : array, shape(M,M)
+ Matrix tangent of A
+
"""
A = asarray(A)
if A.dtype.char not in ['F','D','G']:
@@ -128,7 +211,19 @@
return solve(cosm(A), sinm(A))
def coshm(A):
- """matrix hyperbolic cosine.
+ """Compute the hyperbolic matrix cosine.
+
+ This routine uses expm to compute the matrix exponentials.
+
+ Parameters
+ ----------
+ A : array, shape(M,M)
+
+ Returns
+ -------
+ coshA : array, shape(M,M)
+ Hyperbolic matrix cosine of A
+
"""
A = asarray(A)
if A.dtype.char not in ['F','D','G']:
@@ -137,7 +232,19 @@
return 0.5*(expm(A) + expm(-A))
def sinhm(A):
- """matrix hyperbolic sine.
+ """Compute the hyperbolic matrix sine.
+
+ This routine uses expm to compute the matrix exponentials.
+
+ Parameters
+ ----------
+ A : array, shape(M,M)
+
+ Returns
+ -------
+ sinhA : array, shape(M,M)
+ Hyperbolic matrix sine of A
+
"""
A = asarray(A)
if A.dtype.char not in ['F','D']:
@@ -146,7 +253,19 @@
return 0.5*(expm(A) - expm(-A))
def tanhm(A):
- """matrix hyperbolic tangent.
+ """Compute the hyperbolic matrix tangent.
+
+ This routine uses expm to compute the matrix exponentials.
+
+ Parameters
+ ----------
+ A : array, shape(M,M)
+
+ Returns
+ -------
+ tanhA : array, shape(M,M)
+ Hyperbolic matrix tangent of A
+
"""
A = asarray(A)
if A.dtype.char not in ['F','D']:
@@ -155,11 +274,32 @@
return solve(coshm(A), sinhm(A))
def funm(A,func,disp=1):
- """matrix function for arbitrary callable object func.
+ """Evaluate a matrix function specified by a callable.
+
+ Returns the value of matrix-valued function f at A. The function f
+ is an extension of the scalar-valued function func to matrices.
+
+ Parameters
+ ----------
+ A : array, shape(M,M)
+ Matrix at which to evaluate the function
+ func : callable
+ Callable object that evaluates a scalar function f.
+ Must be vectorized (eg. using vectorize).
+ disp : boolean
+ Print warning if error in the result is estimated large
+ instead of returning estimated error. (Default: True)
+
+ Returns
+ -------
+ fA : array, shape(M,M)
+ Value of the matrix function specified by func evaluated at A
+
+ (if disp == False)
+ errest : float
+ 1-norm of the estimated error, ||err||_1 / ||A||_1
+
"""
- # func should take a vector of arguments (see vectorize if
- # it needs wrapping.
-
# Perform Shur decomposition (lapack ?gees)
A = asarray(A)
if len(A.shape)!=2:
@@ -209,7 +349,28 @@
return F, err
def logm(A,disp=1):
- """Matrix logarithm, inverse of expm."""
+ """Compute matrix logarithm.
+
+ The matrix logarithm is the inverse of expm: expm(logm(A)) == A
+
+ Parameters
+ ----------
+ A : array, shape(M,M)
+ Matrix whose logarithm to evaluate
+ disp : boolean
+ Print warning if error in the result is estimated large
+ instead of returning estimated error. (Default: True)
+
+ Returns
+ -------
+ logA : array, shape(M,M)
+ Matrix logarithm of A
+
+ (if disp == False)
+ errest : float
+ 1-norm of the estimated error, ||err||_1 / ||A||_1
+
+ """
# Compute using general funm but then use better error estimator and
# make one step in improving estimate using a rotation matrix.
A = mat(asarray(A))
@@ -239,7 +400,37 @@
return F, errest
def signm(a,disp=1):
- """matrix sign"""
+ """Matrix sign function.
+
+ Extension of the scalar sign(x) to matrices.
+
+ Parameters
+ ----------
+ A : array, shape(M,M)
+ Matrix at which to evaluate the sign function
+ disp : boolean
+ Print warning if error in the result is estimated large
+ instead of returning estimated error. (Default: True)
+
+ Returns
+ -------
+ sgnA : array, shape(M,M)
+ Value of the sign function at A
+
+ (if disp == False)
+ errest : float
+ 1-norm of the estimated error, ||err||_1 / ||A||_1
+
+ Examples
+ --------
+ >>> from scipy.linalg import signm, eigvals
+ >>> a = [[1,2,3], [1,2,1], [1,1,1]]
+ >>> eigvals(a)
+ array([ 4.12488542+0.j, -0.76155718+0.j, 0.63667176+0.j])
+ >>> eigvals(signm(a))
+ array([-1.+0.j, 1.+0.j, 1.+0.j])
+
+ """
def rounded_sign(x):
rx = real(x)
if rx.dtype.char=='f':
@@ -286,12 +477,29 @@
return S0, errest
def sqrtm(A,disp=1):
- """Matrix square root
+ """Matrix square root.
+
+ Parameters
+ ----------
+ A : array, shape(M,M)
+ Matrix whose square root to evaluate
+ disp : boolean
+ Print warning if error in the result is estimated large
+ instead of returning estimated error. (Default: True)
+
+ Returns
+ -------
+ sgnA : array, shape(M,M)
+ Value of the sign function at A
- If disp is non-zero display warning if singular matrix.
- If disp is zero then return residual ||A-X*X||_F / ||A||_F
+ (if disp == False)
+ errest : float
+ Frobenius norm of the estimated error, ||err||_F / ||A||_F
+ Notes
+ -----
Uses algorithm by Nicholas J. Higham
+
"""
A = asarray(A)
if len(A.shape)!=2:
More information about the Scipy-svn
mailing list