"""Linear Algebra functions implemented as gufuncs, so they broadcast. .. warning:: This module is only for testing, the functionality will be integrated into numpy.linalg proper. ======================= gufuncs_linalg module ======================= gufuncs_linalg implements a series of linear algebra functions as gufuncs. Most of these functions are already present in numpy.linalg, but as they are implemented using gufunc kernels they can be broadcasting. Some parts that are python in numpy.linalg are implemented inside C functions, as well as the iteration used when used on vectors. This can result in faster execution as well. In addition, there are some ufuncs thrown in that implement fused operations over numpy vectors that can result in faster execution on large vector compared to non-fused versions (for example: multiply_add, multiply3). In fact, gufuncs_linalg is a very thin wrapper of python code that wraps the actual kernels (gufuncs). This wrapper was needed in order to provide a sane interface for some functions. Mostly working around limitations on what can be described in a gufunc signature. Things like having one dimension of a result depending on the minimum of two dimensions of the sources (like in svd) or passing an uniform keyword parameter to the whole operation (like UPLO on functions over symmetric/hermitian matrices). The gufunc kernels are in a c module named _umath_linalg, that is imported privately in gufuncs_linalg. ========== Contents ========== Here is an enumeration of the functions. These are the functions exported by the module and should appear in its __all__ attribute. All the functions contain a docstring explaining them in detail. General ======= - inner1d - innerwt - matrix_multiply - quadratic_form Lineal Algebra ============== - det - slogdet - cholesky - eig - eigvals - eigh - eigvalsh - solve - svd - chosolve - inv - poinv Fused Operations ================ - add3 - multiply3 - multiply3_add - multiply_add - multiply_add2 - multiply4 - multiply4_add ================ Error Handling ================ Unlike the numpy.linalg module, this module does not use exceptions to notify errors in the execution of the kernels. As these functions are thougth to be used in a vector way it didn't seem appropriate to raise exceptions on failure of an element. So instead, when an error computing an element occurs its associated result will be set to an invalid value (all NaNs). Exceptions can occur if the arguments fail to map properly to the underlying gufunc (due to signature mismatch, for example). ================================ Notes about the implementation ================================ Where possible, the wrapper functions map directly into a gufunc implementing the computation. That's not always the case, as due to limitations of the gufunc interface some functions cannot be mapped straight into a kernel. Two cases come to mind: - An uniform parameter is needed to configure the way the computation is performed (like UPLO in the functions working on symmetric/hermitian matrices) - svd, where it was impossible to map the function to a gufunc signature. In the case of uniform parameters like UPLO, there are two separate entry points in the C module that imply either 'U' or 'L'. The wrapper just selects the kernel to use by checking the appropriate keyword parameter. This way a function interface similar to numpy.linalg can be kept. In the case of SVD not only there were problems with the support of keyword arguments. There was the added problem of the signature system not being able to cope with the needs of this functions. Just for the singular values a a signature like (m,n)->(min(m,n)) was needed. This has been worked around by implementing different kernels for the cases where min(m,n) == m and where min(m,n) == n. The wrapper code automatically selects the appropriate one. """ from __future__ import division, absolute_import, print_function __all__ = ['inner1d', 'dotc1d', 'innerwt', 'matrix_multiply', 'det', 'slogdet', 'inv', 'cholesky', 'quadratic_form', 'add3', 'multiply3', 'multiply3_add', 'multiply_add', 'multiply_add2', 'multiply4', 'multiply4_add', 'eig', 'eigvals', 'eigh', 'eigvalsh', 'solve', 'svd', 'chosolve', 'poinv'] import numpy as np from . import _umath_linalg as _impl def inner1d(a, b, **kwargs): """ Compute the dot product of vectors over the inner dimension, with broadcasting. Parameters ---------- a : (..., N) array Input array b : (..., N) array Input array Returns ------- inner : (...) array dot product over the inner dimension. Notes ----- Numpy broadcasting rules apply when matching dimensions. Implemented for types single, double, csingle and cdouble. Numpy conversion rules apply. For single and double types this is equivalent to dotc1d. Maps to Blas functions sdot, ddot, cdotu and zdotu. See Also -------- dotc1d : dot product conjugating first vector. innerwt : weighted (i.e. triple) inner product. Examples -------- >>> a = np.arange(1,5).reshape(2,2) >>> b = np.arange(1,8,2).reshape(2,2) >>> res = inner1d(a,b) >>> res.shape (2,) >>> print res [ 7. 43.] """ return _impl.inner1d(a, b, **kwargs) def dotc1d(a, b, **kwargs): """ Compute the dot product of vectors over the inner dimension, conjugating the first vector, with broadcasting Parameters ---------- a : (..., N) array Input array b : (..., N) array Input array Returns ------- dotc : (...) array dot product conjugating the first vector over the inner dimension. Notes ----- Numpy broadcasting rules apply when matching dimensions. Implemented for types single, double, csingle and cdouble. Numpy conversion rules apply. For single and double types this is equivalent to inner1d. Maps to Blas functions sdot, ddot, cdotc and zdotc. See Also -------- inner1d : dot product innerwt : weighted (i.e. triple) inner product. Examples -------- >>> a = np.arange(1,5).reshape(2,2) >>> b = np.arange(1,8,2).reshape(2,2) >>> res = inner1d(a,b) >>> res.shape (2,) >>> print res [ 7. 43.] """ return _impl.dotc1d(a, b, **kwargs) def innerwt(a, b, c, **kwargs): """ Compute the weighted (i.e. triple) inner product, with broadcasting. Parameters ---------- a, b, c : (..., N) array Input arrays Returns ------- inner : (...) array The weighted (i.e. triple) inner product. Notes ----- Numpy broadcasting rules apply. Implemented for types single, double, csingle and cdouble. Numpy conversion rules apply. See Also -------- inner1d : inner product. dotc1d : dot product conjugating first vector. Examples -------- >>> a = np.arange(1,5).reshape(2,2) >>> b = np.arange(1,8,2).reshape(2,2) >>> c = np.arange(0.25,1.20,0.25).reshape(2,2) >>> res = innerwt(a,b,c) >>> res.shape (2,) >>> res array([ 3.25, 39.25]) """ return _impl.innerwt(a, b, c, **kwargs) def matrix_multiply(a,b,**kwargs): """ Compute matrix multiplication, with broadcasting Parameters ---------- a : (..., M, N) array Input array. b : (..., N, P) array Input array. Returns ------- r : (..., M, P) array matrix multiplication of a and b over any number of outer dimensions Notes ----- Numpy broadcasting rules apply. Matrix multiplication is computed using BLAS _gemm functions. Implemented for single, double, csingle and cdouble. Numpy conversion rules apply. Examples -------- >>> a = np.arange(1,17).reshape(2,2,4) >>> b = np.arange(1,25).reshape(2,4,3) >>> res = matrix_multiply(a,b) >>> res.shape (2, 2, 3) >>> res array([[[ 70., 80., 90.], [ 158., 184., 210.]], [[ 750., 792., 834.], [ 1030., 1088., 1146.]]]) """ return _impl.matrix_multiply(a,b,**kwargs) def det(a, **kwargs): """ Compute the determinant of arrays, with broadcasting. Parameters ---------- a : (NDIMS, M, M) array Input array. Its inner dimensions must be those of a square 2-D array. Returns ------- det : (NDIMS) array Determinants of `a` See Also -------- slogdet : Another representation for the determinant, more suitable for large matrices where underflow/overflow may occur Notes ----- Numpy broadcasting rules apply. The determinants are computed via LU factorization using the LAPACK routine _getrf. Implemented for single, double, csingle and cdouble. Numpy conversion rules apply. Examples -------- The determinant of a 2-D array [[a, b], [c, d]] is ad - bc: >>> a = np.array([[1, 2], [3, 4]]) >>> np.allclose(-2.0, det(a)) True >>> a = np.array([[[1, 2], [3, 4]], [[5, 6], [7, 8]] ]) >>> np.allclose(-2.0, det(a)) True """ return _impl.det(a, **kwargs) def slogdet(a, **kwargs): """ Compute the sign and (natural) logarithm of the determinant of an array, with broadcasting. If an array has a very small or very large determinant, then a call to `det` may overflow or underflow. This routine is more robust against such issues, because it computes the logarithm of the determinant rather than the determinant itself Parameters ---------- a : (..., M, M) array Input array. Its inner dimensions must be those of a square 2-D array. Returns ------- sign : (...) array An array of numbers representing the sign of the determinants. For real matrices, this is 1, 0, or -1. For complex matrices, this is a complex number with absolute value 1 (i.e., it is on the unit circle), or else 0. logdet : (...) array The natural log of the absolute value of the determinant. This is always a real type. If the determinant is zero, then `sign` will be 0 and `logdet` will be -Inf. In all cases, the determinant is equal to ``sign * np.exp(logdet)``. See Also -------- det Notes ----- Numpy broadcasting rules apply. The determinants are computed via LU factorization using the LAPACK routine _getrf. Implemented for types single, double, csingle and cdouble. Numpy conversion rules apply. For complex versions `logdet` will be of the associated real type (single for csingle, double for cdouble). Examples -------- The determinant of a 2-D array [[a, b], [c, d]] is ad - bc: >>> a = np.array([[1, 2], [3, 4]]) >>> (sign, logdet) = slogdet(a) >>> sign.shape () >>> logdet.shape () >>> np.allclose(-2.0, sign * np.exp(logdet)) True >>> a = np.array([[[1, 2], [3, 4]], [[5, 6], [7, 8]] ]) >>> (sign, logdet) = slogdet(a) >>> sign.shape (2,) >>> logdet.shape (2,) >>> np.allclose(-2.0, sign * np.exp(logdet)) True """ return _impl.slogdet(a, **kwargs) def inv(a, **kwargs): """ Compute the (multiplicative) inverse of matrices, with broadcasting. Given a square matrix `a`, return the matrix `ainv` satisfying ``matrix_multiply(a, ainv) = matrix_multiply(ainv, a) = Identity matrix`` Parameters ---------- a : (..., M, M) array Matrices to be inverted Returns ------- ainv : (..., M, M) array (Multiplicative) inverse of the `a` matrices. Notes ----- Numpy broadcasting rules apply. Implemented for types single, double, csingle and cdouble. Numpy conversion rules apply. Singular matrices and thus, not invertible, result in an array of NaNs. See Also -------- poinv : compute the multiplicative inverse of hermitian/symmetric matrices, using cholesky decomposition. Examples -------- >>> a = np.array([[1, 2], [3, 4]]) >>> ainv = inv(a) >>> np.allclose(matrix_multiply(a, ainv), np.eye(2)) True >>> np.allclose(matrix_multiply(ainv, a), np.eye(2)) True """ return _impl.inv(a, **kwargs) def cholesky(a, UPLO='L', **kwargs): """ Compute the cholesky decomposition of `a`, with broadcasting The Cholesky decomposition (or Cholesky triangle) is a decomposition of a Hermitian, positive-definite matrix into the product of a lower triangular matrix and its conjugate transpose. A = LL* where L* is the positive-definite matrix. Parameters ---------- a : (..., M, M) array Matrices for which compute the cholesky decomposition Returns ------- l : (..., M, M) array Matrices for each element where each entry is the lower triangular matrix with strictly positive diagonal entries such that a = ll* for all outer dimensions See Also -------- chosolve : solve a system using cholesky decomposition poinv : compute the inverse of a matrix using cholesky decomposition Notes ----- Numpy broadcasting rules apply. Implemented for types single, double, csingle and cdouble. Numpy conversion rules apply. Decomposition is performed using LAPACK routine _potrf. For elements where the LAPACK routine fails, the result will be set to NaNs. If an element of the source array is not a positive-definite matrix the result for that element is undefined. Examples -------- >>> A = np.array([[1,-2j],[2j,5]]) >>> A array([[ 1.+0.j, 0.-2.j], [ 0.+2.j, 5.+0.j]]) >>> L = cholesky(A) >>> L array([[ 1.+0.j, 0.+0.j], [ 0.+2.j, 1.+0.j]]) """ if 'L' == UPLO: gufunc = _impl.cholesky_lo else: gufunc = _impl.cholesky_up return gufunc(a, **kwargs) def eig(a, **kwargs): """ Compute the eigenvalues and right eigenvectors of square arrays, with broadcasting Parameters ---------- a : (..., M, M) array Matrices for which the eigenvalues and right eigenvectors will be computed Returns ------- w : (..., M) array The eigenvalues, each repeated according to its multiplicity. The eigenvalues are not necessarily ordered. The resulting array will be always be of complex type. When `a` is real the resulting eigenvalues will be real (0 imaginary part) or occur in conjugate pairs v : (..., M, M) array The normalized (unit "length") eigenvectors, such that the column ``v[:,i]`` is the eigenvector corresponding to the eigenvalue ``w[i]``. See Also -------- eigvals : eigenvalues of general arrays. eigh : eigenvalues and right eigenvectors of symmetric/hermitian arrays. eigvalsh : eigenvalues of symmetric/hermitian arrays. Notes ----- Numpy broadcasting rules apply. Implemented for types single, double, csingle and cdouble. Numpy conversion rules apply. Eigenvalues and eigenvectors for single and double versions will always be typed csingle and cdouble, even if all the results are real (imaginary part will be 0). This is implemented using the _geev LAPACK routines which compute the eigenvalues and eigenvectors of general square arrays. For elements where the LAPACK routine fails, the result will be set to NaNs. Examples -------- First, a utility function to check if eigvals/eigvectors are correct. This checks the definition of eigenvectors. For each eigenvector v with associated eigenvalue w of a matrix M the following equality must hold: Mv == wv >>> def check_eigen(M, w, v): ... '''vectorial check of Mv==wv''' ... lhs = matrix_multiply(M, v) ... rhs = w*v ... return np.allclose(lhs, rhs) (Almost) Trivial example with real e-values and e-vectors. Note the complex types of the results >>> M = np.diag((1,2,3)).astype(float) >>> w, v = eig(M) >>> check_eigen(M, w, v) True Real matrix possessing complex e-values and e-vectors; note that the e-values are complex conjugates of each other. >>> M = np.array([[1, -1], [1, 1]]) >>> w, v = eig(M) >>> check_eigen(M, w, v) True Complex-valued matrix with real e-values (but complex-valued e-vectors); note that a.conj().T = a, i.e., a is Hermitian. >>> M = np.array([[1, 1j], [-1j, 1]]) >>> w, v = eig(M) >>> check_eigen(M, w, v) True """ return _impl.eig(a, **kwargs) def eigvals(a, **kwargs): """ Compute the eigenvalues of general matrices, with broadcasting. Main difference between `eigvals` and `eig`: the eigenvectors aren't returned. Parameters ---------- a : (..., M, M) array Matrices whose eigenvalues will be computed Returns ------- w : (..., M) array The eigenvalues, each repeated according to its multiplicity. The eigenvalues are not necessarily ordered. The resulting array will be always be of complex type. When `a` is real the resulting eigenvalues will be real (0 imaginary part) or occur in conjugate pairs See Also -------- eig : eigenvalues and right eigenvectors of general arrays. eigh : eigenvalues and right eigenvectors of symmetric/hermitian arrays. eigvalsh : eigenvalues of symmetric/hermitian arrays. Notes ----- Numpy broadcasting rules apply. Implemented for types single, double, csingle and cdouble. Numpy conversion rules apply. Eigenvalues for single and double versions will always be typed csingle and cdouble, even if all the results are real (imaginary part will be 0). This is implemented using the _geev LAPACK routines which compute the eigenvalues and eigenvectors of general square arrays. For elements where the LAPACK routine fails, the result will be set to NaNs. Examples -------- Eigenvalues for a diagonal matrix are its diagonal elements >>> D = np.diag((-1,1)) >>> eigvals(D) array([-1.+0.j, 1.+0.j]) Multiplying on the left by an orthogonal matrix, `Q`, and on the right by `Q.T` (the transpose of `Q` preserves the eigenvalues of the original matrix >>> x = np.random.random() >>> Q = np.array([[np.cos(x), -np.sin(x)], [np.sin(x), np.cos(x)]]) >>> A = matrix_multiply(Q, D) >>> A = matrix_multiply(A, Q.T) >>> eigvals(A) array([-1.+0.j, 1.+0.j]) """ return _impl.eigvals(a, **kwargs) def quadratic_form(u,Q,v, **kwargs): """ Compute the quadratic form uQv, with broadcasting Parameters ---------- u : (..., M) array The u vectors of the quadratic form uQv Q : (..., M, N) array The Q matrices of the quadratic form uQv v : (..., N) array The v vectors of the quadratic form uQv Returns ------- qf : (...) array The result of the quadratic forms Notes ----- Numpy broadcasting rules apply. Implemented for types single, double, csingle and cdouble. Numpy conversion rules apply. This is similar to PDL inner2 Examples -------- The result in absence of broadcasting is just as np.dot(np.dot(u,Q),v) or np.dot(u, np.dot(Q,v)) >>> u = np.array([2., 3.]) >>> Q = np.array([[1.,1.], [0.,1.]]) >>> v = np.array([1.,2.]) >>> quadratic_form(u,Q,v) 12.0 >>> np.dot(np.dot(u,Q),v) 12.0 >>> np.dot(u, np.dot(Q,v)) 12.0 """ return _impl.quadratic_form(u, Q, v, **kwargs) def add3(a, b, c, **kwargs): """ Element-wise addition of 3 arrays: a + b + c. Parameters ---------- a, b, c : (...) array arrays with the addends Returns ------- add3 : (...) array resulting element-wise addition. Notes ----- Numpy broadcasting rules apply. Implemented for types single, double, csingle and cdouble. Numpy conversion rules apply. See Also -------- multiply3 : element-wise three-way multiplication. multiply3_add : element-wise three-way multiplication and addition. multiply_add : element-wise multiply-add. multiply_add2 : element-wise multiplication with two additions. multiply4 : element-wise four-way multiplication multiply4_add : element-wise four-way multiplication and addition, Examples -------- >>> a = np.linspace(1.0, 30.0, 30) >>> add3(a[0::3], a[1::3], a[2::3]) array([ 6., 15., 24., 33., 42., 51., 60., 69., 78., 87.]) """ return _impl.add3(a, b, c, **kwargs) def multiply3(a, b, c, **kwargs): """ Element-wise multiplication of 3 arrays: a*b*c. Parameters ---------- a, b, c : (...) array arrays with the factors Returns ------- m3 : (...) array resulting element-wise product Notes ----- Numpy broadcasting rules apply. Implemented for types single, double, csingle and cdouble. Numpy conversion rules apply. See Also -------- add3 : element-wise three-was addition multiply3_add : element-wise three-way multiplication and addition. multiply_add : element-wise multiply-add. multiply_add2 : element-wise multiplication with two additions. multiply4 : element-wise four-way multiplication multiply4_add : element-wise four-way multiplication and addition, Examples -------- >>> a = np.linspace(1.0, 10.0, 10) >>> multiply3(a, 1.01, a) array([ 1.01, 4.04, 9.09, 16.16, 25.25, 36.36, 49.49, 64.64, 81.81, 101. ]) """ return _impl.multiply3(a, b, c, **kwargs) def multiply3_add(a, b, c, d, **kwargs): """ Element-wise multiplication of 3 arrays adding an element of the a 4th array to the result: a*b*c + d Parameters ---------- a, b, c : (...) array arrays with the factors d : (...) array array with the addend Returns ------- m3a : (...) array resulting element-wise addition Notes ----- Numpy broadcasting rules apply. Implemented for types single, double, csingle and cdouble. Numpy conversion rules apply. See Also -------- add3 : element-wise three-was addition multiply3 : element-wise three-way multiplication. multiply3_add : element-wise three-way multiplication and addition. multiply_add : element-wise multiply-add. multiply_add2 : element-wise multiplication with two additions. multiply4 : element-wise four-way multiplication multiply4_add : element-wise four-way multiplication and addition, Examples -------- >>> a = np.linspace(1.0, 10.0, 10) >>> multiply3_add(a, 1.01, a, 42e-4) array([ 1.0142, 4.0442, 9.0942, 16.1642, 25.2542, 36.3642, 49.4942, 64.6442, 81.8142, 101.0042]) """ return _impl.multiply3_add(a, b, c, d, **kwargs) def multiply_add(a, b, c, **kwargs): """ Element-wise addition of 3 arrays Parameters ---------- a, b, c : (...) array arrays with the addends Returns ------- add3 : (...) array resulting element-wise addition Notes ----- Numpy broadcasting rules apply. Implemented for types single, double, csingle and cdouble. Numpy conversion rules apply. See Also -------- add3 : element-wise three-was addition multiply3 : element-wise three-way multiplication. multiply3_add : element-wise three-way multiplication and addition. multiply_add : element-wise multiply-add. multiply_add2 : element-wise multiplication with two additions. multiply4 : element-wise four-way multiplication multiply4_add : element-wise four-way multiplication and addition, Examples -------- >>> a = np.linspace(1.0, 10.0, 10) >>> multiply_add(a, a, 42e-4) array([ 1.0042, 4.0042, 9.0042, 16.0042, 25.0042, 36.0042, 49.0042, 64.0042, 81.0042, 100.0042]) """ return _impl.multiply_add(a, b, c, **kwargs) def multiply_add2(a, b, c, d, **kwargs): """ Element-wise addition of 3 arrays Parameters ---------- a, b, c : (...) array arrays with the addends Returns ------- add3 : (...) array resulting element-wise addition Notes ----- Numpy broadcasting rules apply. Implemented for types single, double, csingle and cdouble. Numpy conversion rules apply. See Also -------- add3 : element-wise three-was addition multiply3 : element-wise three-way multiplication. multiply3_add : element-wise three-way multiplication and addition. multiply_add : element-wise multiply-add. multiply_add2 : element-wise multiplication with two additions. multiply4 : element-wise four-way multiplication multiply4_add : element-wise four-way multiplication and addition, Examples -------- >>> a = np.linspace(1.0, 10.0, 10) >>> multiply_add2(a, a, a, 42e-4) array([ 2.0042, 6.0042, 12.0042, 20.0042, 30.0042, 42.0042, 56.0042, 72.0042, 90.0042, 110.0042]) """ return _impl.multiply_add2(a, b, c, d, **kwargs) def multiply4(a, b, c, d, **kwargs): """ Element-wise addition of 3 arrays Parameters ---------- a, b, c : (...) array arrays with the addends Returns ------- add3 : (...) array resulting element-wise addition Notes ----- Numpy broadcasting rules apply. Implemented for types single, double, csingle and cdouble. Numpy conversion rules apply. See Also -------- add3 : element-wise three-was addition multiply3 : element-wise three-way multiplication. multiply3_add : element-wise three-way multiplication and addition. multiply_add : element-wise multiply-add. multiply_add2 : element-wise multiplication with two additions. multiply4 : element-wise four-way multiplication multiply4_add : element-wise four-way multiplication and addition, Examples -------- >>> a = np.linspace(1.0, 10.0, 10) >>> multiply4(a, a, a[::-1], 1.0001) array([ 10.001 , 36.0036, 72.0072, 112.0112, 150.015 , 180.018 , 196.0196, 192.0192, 162.0162, 100.01 ]) """ return _impl.multiply4(a, b, c, d, **kwargs) def multiply4_add(a, b, c, d, e, **kwargs): """ Element-wise addition of 3 arrays Parameters ---------- a, b, c : (...) array arrays with the addends Returns ------- add3 : (...) array resulting element-wise addition Notes ----- Numpy broadcasting rules apply. Implemented for types single, double, csingle and cdouble. Numpy conversion rules apply. See Also -------- add3 : element-wise three-was addition multiply3 : element-wise three-way multiplication. multiply3_add : element-wise three-way multiplication and addition. multiply_add : element-wise multiply-add. multiply_add2 : element-wise multiplication with two additions. multiply4 : element-wise four-way multiplication multiply4_add : element-wise four-way multiplication and addition, Examples -------- >>> a = np.linspace(1.0, 10.0, 10) >>> multiply4_add(a, a, a[::-1], 1.01, 42e-4) array([ 10.1042, 36.3642, 72.7242, 113.1242, 151.5042, 181.8042, 197.9642, 193.9242, 163.6242, 101.0042]) """ return _impl.multiply4_add(a, b, c, d, e,**kwargs) def eigh(A, UPLO='L', **kw_args): """ Computes the eigenvalues and eigenvectors for the square matrices in the inner dimensions of A, being those matrices symmetric/hermitian. Parameters ---------- A : (..., M, M) array Hermitian/Symmetric matrices whose eigenvalues and eigenvectors are to be computed. UPLO : {'L', 'U'}, optional Specifies whether the calculation is done with the lower triangular part of the elements in `A` ('L', default) or the upper triangular part ('U'). Returns ------- w : (..., M) array The eigenvalues, not necessarily ordered. v : (..., M, M) array The inner dimensions contain matrices with the normalized eigenvectors as columns. The column-numbers are coherent with the associated eigenvalue. Notes ----- Numpy broadcasting rules apply. The eigenvalues/eigenvectors are computed using LAPACK routines _ssyevd, _heevd For elements where the LAPACK routine fails, the result will be set to NaNs. Implemented for single, double, csingle and cdouble. Numpy conversion rules apply. Unlike eig, the results for single and double will always be single and doubles. It is not possible for symmetrical real matrices to result in complex eigenvalues/eigenvectors See Also -------- eigvalsh : eigenvalues of symmetric/hermitian arrays. eig : eigenvalues and right eigenvectors for general matrices. eigvals : eigenvalues for general matrices. Examples -------- First, a utility function to check if eigvals/eigvectors are correct. This checks the definition of eigenvectors. For each eigenvector v with associated eigenvalue w of a matrix M the following equality must hold: Mv == wv >>> def check_eigen(M, w, v): ... '''vectorial check of Mv==wv''' ... lhs = matrix_multiply(M, v) ... rhs = w*v ... return np.allclose(lhs, rhs) A simple example that computes eigenvectors and eigenvalues of a hermitian matrix and checks that A*v = v*w for both pairs of eignvalues(w) and eigenvectors(v) >>> M = np.array([[1, -2j], [2j, 1]]) >>> w, v = eigh(M) >>> check_eigen(M, w, v) True """ if 'L' == UPLO: gufunc = _impl.eigh_lo else: gufunc = _impl.eigh_up return gufunc(A, **kw_args) def eigvalsh(A, UPLO='L', **kw_args): """ Computes the eigenvalues for the square matrices in the inner dimensions of A, being those matrices symmetric/hermitian. Parameters ---------- A : (..., M, M) array Hermitian/Symmetric matrices whose eigenvalues and eigenvectors are to be computed. UPLO : {'L', 'U'}, optional Specifies whether the calculation is done with the lower triangular part of the elements in `A` ('L', default) or the upper triangular part ('U'). Returns ------- w : (..., M) array The eigenvalues, not necessarily ordered. Notes ----- Numpy broadcasting rules apply. The eigenvalues are computed using LAPACK routines _ssyevd, _heevd For elements where the LAPACK routine fails, the result will be set to NaNs. Implemented for single, double, csingle and cdouble. Numpy conversion rules apply. Unlike eigvals, the results for single and double will always be single and doubles. It is not possible for symmetrical real matrices to result in complex eigenvalues. See Also -------- eigh : eigenvalues of symmetric/hermitian arrays. eig : eigenvalues and right eigenvectors for general matrices. eigvals : eigenvalues for general matrices. Examples -------- eigvalsh results should be the same as the eigenvalues returned by eigh >>> a = np.array([[1, -2j], [2j, 5]]) >>> w, v = eigh(a) >>> np.allclose(w, eigvalsh(a)) True eigvalsh on an identity matrix is all ones >>> eigvalsh(np.eye(6)) array([ 1., 1., 1., 1., 1., 1.]) """ if ('L' == UPLO): gufunc = _impl.eigvalsh_lo else: gufunc = _impl.eigvalsh_up return gufunc(A,**kw_args) def solve(A,B,**kw_args): """ Solve the linear matrix equations on the inner dimensions. Computes the "exact" solution, `x`. of the well-determined, i.e., full rank, linear matrix equations `ax = b`. Parameters ---------- A : (..., M, M) array Coefficient matrices. B : (..., M, N) array Ordinate or "dependent variable" values. Returns ------- X : (..., M, N) array Solutions to the system A X = B for all the outer dimensions Notes ----- Numpy broadcasting rules apply. The solutions are computed using LAPACK routine _gesv For elements where the LAPACK routine fails, the result will be set to NaNs. Implemented for single, double, csingle and cdouble. Numpy conversion rules apply. See Also -------- chosolve : solve a system using cholesky decomposition (for equations having symmetric/hermitian coefficient matrices) Examples -------- Solve the system of equations ``3 * x0 + x1 = 9`` and ``x0 + 2 * x1 = 8``: >>> a = np.array([[3,1], [1,2]]) >>> b = np.array([9,8]) >>> x = solve(a, b) >>> x array([ 2., 3.]) Check that the solution is correct: >>> np.allclose(np.dot(a, x), b) True """ if len(B.shape) == (len(A.shape) - 1): gufunc = _impl.solve1 else: gufunc = _impl.solve return gufunc(A,B,**kw_args) def svd(a, full_matrices=1, compute_uv=1 ,**kw_args): """ Singular Value Decomposition on the inner dimensions. Factors the matrices in `a` as ``u * np.diag(s) * v``, where `u` and `v` are unitary and `s` is a 1-d array of `a`'s singular values. Parameters ---------- a : (..., M, N) array The array of matrices to decompose. full_matrices : bool, optional If True (default), `u` and `v` have the shapes (`M`, `M`) and (`N`, `N`), respectively. Otherwise, the shapes are (`M`, `K`) and (`K`, `N`), respectively, where `K` = min(`M`, `N`). compute_uv : bool, optional Whether or not to compute `u` and `v` in addition to `s`. True by default. Returns ------- u : { (..., M, M), (..., M, K) } array Unitary matrices. The actual shape depends on the value of ``full_matrices``. Only returned when ``compute_uv`` is True. s : (..., K) array The singular values for every matrix, sorted in descending order. v : { (..., N, N), (..., K, N) } array Unitary matrices. The actual shape depends on the value of ``full_matrices``. Only returned when ``compute_uv`` is True. Notes ----- Numpy broadcasting rules apply. Singular Value Decomposition is performed using LAPACK routine _gesdd For elements where the LAPACK routine fails, the result will be set to NaNs. Implemented for types single, double, csingle and cdouble. Numpy conversion rules apply. Examples -------- >>> a = np.random.randn(9, 6) + 1j*np.random.randn(9, 6) Reconstruction based on full SVD: >>> U, s, V = svd(a, full_matrices=True) >>> U.shape, V.shape, s.shape ((9, 9), (6, 6), (6,)) >>> S = np.zeros((9, 6), dtype=complex) >>> S[:6, :6] = np.diag(s) >>> np.allclose(a, np.dot(U, np.dot(S, V))) True Reconstruction based on reduced SVD: >>> U, s, V = svd(a, full_matrices=False) >>> U.shape, V.shape, s.shape ((9, 6), (6, 6), (6,)) >>> S = np.diag(s) >>> np.allclose(a, np.dot(U, np.dot(S, V))) True """ m = a.shape[-2] n = a.shape[-1] if 1 == compute_uv: if 1 == full_matrices: if m < n: gufunc = _impl.svd_m_f else: gufunc = _impl.svd_n_f else: if m < n: gufunc = _impl.svd_m_s else: gufunc = _impl.svd_n_s else: if m < n: gufunc = _impl.svd_m else: gufunc = _impl.svd_n return gufunc(a, **kw_args) def chosolve(A, B, UPLO='L', **kw_args): """ Solve the linear matrix equations on the inner dimensions, using cholesky decomposition. Computes the "exact" solution, `x`. of the well-determined, i.e., full rank, linear matrix equations `ax = b`, where a is a symmetric/hermitian positive-definite matrix. Parameters ---------- A : (..., M, M) array Coefficient symmetric/hermitian positive-definite matrices. B : (..., M, N) array Ordinate or "dependent variable" values. UPLO : {'L', 'U'}, optional Specifies whether the calculation is done with the lower triangular part of the elements in `A` ('L', default) or the upper triangular part ('U'). Returns ------- X : (..., M, N) array Solutions to the system A X = B for all elements in the outer dimensions Notes ----- Numpy broadcasting rules apply. The solutions are computed using LAPACK routines _potrf, _potrs For elements where the LAPACK routine fails, the result will be set to NaNs. Implemented for single, double, csingle and cdouble. Numpy conversion rules apply. See Also -------- solve : solve a system using cholesky decomposition (for equations having symmetric/hermitian coefficient matrices) Examples -------- Solve the system of equations ``3 * x0 + x1 = 9`` and ``x0 + 2 * x1 = 8``: (note the matrix is symmetric in this system) >>> a = np.array([[3,1], [1,2]]) >>> b = np.array([9,8]) >>> x = solve(a, b) >>> x array([ 2., 3.]) Check that the solution is correct: >>> np.allclose(np.dot(a, x), b) True """ if len(B.shape) == (len(A.shape) - 1): if 'L' == UPLO: gufunc = _impl.chosolve1_lo else: gufunc = _impl.chosolve1_up else: if 'L' == UPLO: gufunc = _impl.chosolve_lo else: gufunc = _impl.chosolve_up return gufunc(A, B, **kw_args) def poinv(A, UPLO='L', **kw_args): """ Compute the (multiplicative) inverse of symmetric/hermitian positive definite matrices, with broadcasting. Given a square symmetic/hermitian positive-definite matrix `a`, return the matrix `ainv` satisfying ``matrix_multiply(a, ainv) = matrix_multiply(ainv, a) = Identity matrix``. Parameters ---------- a : (..., M, M) array Symmetric/hermitian postive definite matrices to be inverted. Returns ------- ainv : (..., M, M) array (Multiplicative) inverse of the `a` matrices. Notes ----- Numpy broadcasting rules apply. The inverse is computed using LAPACK routines _potrf, _potri For elements where the LAPACK routine fails, the result will be set to NaNs. Implemented for types single, double, csingle and cdouble. Numpy conversion rules apply. See Also -------- inv : compute the multiplicative inverse of general matrices. Examples -------- >>> a = np.array([[5, 3], [3, 5]]) >>> ainv = poinv(a) >>> np.allclose(matrix_multiply(a, ainv), np.eye(2)) True >>> np.allclose(matrix_multiply(ainv, a), np.eye(2)) True """ if 'L' == UPLO: gufunc = _impl.poinv_lo else: gufunc = _impl.poinv_up return gufunc(A, **kw_args); if __name__ == "__main__": import doctest doctest.testmod()