Linear algebra ( numpy.linalg )#
The NumPy linear algebra functions rely on BLAS and LAPACK to provide efficient low level implementations of standard linear algebra algorithms. Those libraries may be provided by NumPy itself using C versions of a subset of their reference implementations but, when possible, highly optimized libraries that take advantage of specialized processor functionality are preferred. Examples of such libraries are OpenBLAS, MKL (TM), and ATLAS. Because those libraries are multithreaded and processor dependent, environmental variables and external packages such as threadpoolctl may be needed to control the number of threads or specify the processor architecture.
The SciPy library also contains a linalg submodule, and there is overlap in the functionality provided by the SciPy and NumPy submodules. SciPy contains functions not found in numpy.linalg , such as functions related to LU decomposition and the Schur decomposition, multiple ways of calculating the pseudoinverse, and matrix transcendentals such as the matrix logarithm. Some functions that exist in both have augmented functionality in scipy.linalg . For example, scipy.linalg.eig can take a second matrix argument for solving generalized eigenvalue problems. Some functions in NumPy, however, have more flexible broadcasting options. For example, numpy.linalg.solve can handle “stacked” arrays, while scipy.linalg.solve accepts only a single square array as its first argument.
The term matrix as it is used on this page indicates a 2d numpy.array object, and not a numpy.matrix object. The latter is no longer recommended, even for linear algebra. See the matrix object documentation for more information.
The @ operator#
Introduced in NumPy 1.10.0, the @ operator is preferable to other methods when computing the matrix product between 2d arrays. The numpy.matmul function implements the @ operator.
Matrix and vector products#
Dot product of two arrays.
Compute the dot product of two or more arrays in a single function call, while automatically selecting the fastest evaluation order.
Return the dot product of two vectors.
Inner product of two arrays.
Compute the outer product of two vectors.
Matrix product of two arrays.
Compute tensor dot product along specified axes.
einsum (subscripts, *operands[, out, dtype, . ])
Evaluates the Einstein summation convention on the operands.
einsum_path (subscripts, *operands[, optimize])
Evaluates the lowest cost contraction order for an einsum expression by considering the creation of intermediate arrays.
Raise a square matrix to the (integer) power n.
Kronecker product of two arrays.
Decompositions#
Compute the qr factorization of a matrix.
linalg.svd (a[, full_matrices, compute_uv, . ])
Singular Value Decomposition.
Matrix eigenvalues#
Compute the eigenvalues and right eigenvectors of a square array.
Return the eigenvalues and eigenvectors of a complex Hermitian (conjugate symmetric) or a real symmetric matrix.
Compute the eigenvalues of a general matrix.
Compute the eigenvalues of a complex Hermitian or real symmetric matrix.
Norms and other numbers#
Compute the condition number of a matrix.
Compute the determinant of an array.
Return matrix rank of array using SVD method
Compute the sign and (natural) logarithm of the determinant of an array.
trace (a[, offset, axis1, axis2, dtype, out])
Return the sum along diagonals of the array.
Solving equations and inverting matrices#
Solve a linear matrix equation, or system of linear scalar equations.
Solve the tensor equation a x = b for x.
Return the least-squares solution to a linear matrix equation.
Compute the (multiplicative) inverse of a matrix.
Compute the (Moore-Penrose) pseudo-inverse of a matrix.
Compute the ‘inverse’ of an N-dimensional array.
Exceptions#
Generic Python-exception-derived object raised by linalg functions.
Linear algebra on several matrices at once#
Several of the linear algebra routines listed above are able to compute results for several matrices at once, if they are stacked into the same array.
This is indicated in the documentation via input parameter specifications such as a : (. M, M) array_like . This means that if for instance given an input array a.shape == (N, M, M) , it is interpreted as a “stack” of N matrices, each of size M-by-M. Similar specification applies to return values, for instance the determinant has det : (. ) and will in this case return an array of shape det(a).shape == (N,) . This generalizes to linear algebra operations on higher-dimensional arrays: the last 1 or 2 dimensions of a multidimensional array are interpreted as vectors or matrices, as appropriate for each operation.
numpy.linalg.eig#
Compute the eigenvalues and right eigenvectors of a square array.
Parameters : a (…, M, M) array
Matrices for which the eigenvalues and right eigenvectors will be computed
Returns : A namedtuple with the following attributes: eigenvalues (…, M) array
The eigenvalues, each repeated according to its multiplicity. The eigenvalues are not necessarily ordered. The resulting array will be of complex type, unless the imaginary part is zero in which case it will be cast to a real type. When a is real the resulting eigenvalues will be real (0 imaginary part) or occur in conjugate pairs
eigenvectors (…, M, M) array
The normalized (unit “length”) eigenvectors, such that the column eigenvectors[:,i] is the eigenvector corresponding to the eigenvalue eigenvalues[i] .
If the eigenvalue computation does not converge.
eigenvalues of a non-symmetric array.
eigenvalues and eigenvectors of a real symmetric or complex Hermitian (conjugate symmetric) array.
eigenvalues of a real symmetric or complex Hermitian (conjugate symmetric) array.
Similar function in SciPy that also solves the generalized eigenvalue problem.
Best choice for unitary and other non-Hermitian normal matrices.
Broadcasting rules apply, see the numpy.linalg documentation for details.
This is implemented using the _geev LAPACK routines which compute the eigenvalues and eigenvectors of general square arrays.
The number w is an eigenvalue of a if there exists a vector v such that a @ v = w * v . Thus, the arrays a, eigenvalues, and eigenvectors satisfy the equations a @ eigenvectors[:,i] = eigenvalues[i] * eigenvalues[:,i] for \(i \in \\) .
The array eigenvectors may not be of maximum rank, that is, some of the columns may be linearly dependent, although round-off error may obscure that fact. If the eigenvalues are all different, then theoretically the eigenvectors are linearly independent and a can be diagonalized by a similarity transformation using eigenvectors, i.e, inv(eigenvectors) @ a @ eigenvectors is diagonal.
For non-Hermitian normal matrices the SciPy function scipy.linalg.schur is preferred because the matrix eigenvectors is guaranteed to be unitary, which is not the case when using eig . The Schur factorization produces an upper triangular matrix rather than a diagonal matrix, but for normal matrices only the diagonal of the upper triangular matrix is needed, the rest is roundoff error.
Finally, it is emphasized that eigenvectors consists of the right (as in right-hand side) eigenvectors of a. A vector y satisfying y.T @ a = z * y.T for some number z is called a left eigenvector of a, and, in general, the left and right eigenvectors of a matrix are not necessarily the (perhaps conjugate) transposes of each other.
G. Strang, Linear Algebra and Its Applications, 2nd Ed., Orlando, FL, Academic Press, Inc., 1980, Various pp.
>>> from numpy import linalg as LA
(Almost) trivial example with real eigenvalues and eigenvectors.
>>> eigenvalues, eigenvectors = LA.eig(np.diag((1, 2, 3))) >>> eigenvalues array([1., 2., 3.]) >>> eigenvectors array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]])
Real matrix possessing complex eigenvalues and eigenvectors; note that the eigenvalues are complex conjugates of each other.
>>> eigenvalues, eigenvectors = LA.eig(np.array([[1, -1], [1, 1]])) >>> eigenvalues array([1.+1.j, 1.-1.j]) >>> eigenvectors array([[0.70710678+0.j , 0.70710678-0.j ], [0. -0.70710678j, 0. +0.70710678j]])
Complex-valued matrix with real eigenvalues (but complex-valued eigenvectors); note that a.conj().T == a , i.e., a is Hermitian.
>>> a = np.array([[1, 1j], [-1j, 1]]) >>> eigenvalues, eigenvectors = LA.eig(a) >>> eigenvalues array([2.+0.j, 0.+0.j]) >>> eigenvectors array([[ 0. +0.70710678j, 0.70710678+0.j ], # may vary [ 0.70710678+0.j , -0. +0.70710678j]])
Be careful about round-off error!
>>> a = np.array([[1 + 1e-9, 0], [0, 1 - 1e-9]]) >>> # Theor. eigenvalues are 1 +/- 1e-9 >>> eigenvalues, eigenvectors = LA.eig(a) >>> eigenvalues array([1., 1.]) >>> eigenvectors array([[1., 0.], [0., 1.]])
numpy.linalg.eig#
Compute the eigenvalues and right eigenvectors of a square array.
Parameters : a (…, M, M) array
Matrices for which the eigenvalues and right eigenvectors will be computed
Returns : A namedtuple with the following attributes: eigenvalues (…, M) array
The eigenvalues, each repeated according to its multiplicity. The eigenvalues are not necessarily ordered. The resulting array will be of complex type, unless the imaginary part is zero in which case it will be cast to a real type. When a is real the resulting eigenvalues will be real (0 imaginary part) or occur in conjugate pairs
eigenvectors (…, M, M) array
The normalized (unit “length”) eigenvectors, such that the column eigenvectors[:,i] is the eigenvector corresponding to the eigenvalue eigenvalues[i] .
If the eigenvalue computation does not converge.
eigenvalues of a non-symmetric array.
eigenvalues and eigenvectors of a real symmetric or complex Hermitian (conjugate symmetric) array.
eigenvalues of a real symmetric or complex Hermitian (conjugate symmetric) array.
Similar function in SciPy that also solves the generalized eigenvalue problem.
Best choice for unitary and other non-Hermitian normal matrices.
Broadcasting rules apply, see the numpy.linalg documentation for details.
This is implemented using the _geev LAPACK routines which compute the eigenvalues and eigenvectors of general square arrays.
The number w is an eigenvalue of a if there exists a vector v such that a @ v = w * v . Thus, the arrays a, eigenvalues, and eigenvectors satisfy the equations a @ eigenvectors[:,i] = eigenvalues[i] * eigenvalues[:,i] for \(i \in \\) .
The array eigenvectors may not be of maximum rank, that is, some of the columns may be linearly dependent, although round-off error may obscure that fact. If the eigenvalues are all different, then theoretically the eigenvectors are linearly independent and a can be diagonalized by a similarity transformation using eigenvectors, i.e, inv(eigenvectors) @ a @ eigenvectors is diagonal.
For non-Hermitian normal matrices the SciPy function scipy.linalg.schur is preferred because the matrix eigenvectors is guaranteed to be unitary, which is not the case when using eig . The Schur factorization produces an upper triangular matrix rather than a diagonal matrix, but for normal matrices only the diagonal of the upper triangular matrix is needed, the rest is roundoff error.
Finally, it is emphasized that eigenvectors consists of the right (as in right-hand side) eigenvectors of a. A vector y satisfying y.T @ a = z * y.T for some number z is called a left eigenvector of a, and, in general, the left and right eigenvectors of a matrix are not necessarily the (perhaps conjugate) transposes of each other.
G. Strang, Linear Algebra and Its Applications, 2nd Ed., Orlando, FL, Academic Press, Inc., 1980, Various pp.
>>> from numpy import linalg as LA
(Almost) trivial example with real eigenvalues and eigenvectors.
>>> eigenvalues, eigenvectors = LA.eig(np.diag((1, 2, 3))) >>> eigenvalues array([1., 2., 3.]) >>> eigenvectors array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]])
Real matrix possessing complex eigenvalues and eigenvectors; note that the eigenvalues are complex conjugates of each other.
>>> eigenvalues, eigenvectors = LA.eig(np.array([[1, -1], [1, 1]])) >>> eigenvalues array([1.+1.j, 1.-1.j]) >>> eigenvectors array([[0.70710678+0.j , 0.70710678-0.j ], [0. -0.70710678j, 0. +0.70710678j]])
Complex-valued matrix with real eigenvalues (but complex-valued eigenvectors); note that a.conj().T == a , i.e., a is Hermitian.
>>> a = np.array([[1, 1j], [-1j, 1]]) >>> eigenvalues, eigenvectors = LA.eig(a) >>> eigenvalues array([2.+0.j, 0.+0.j]) >>> eigenvectors array([[ 0. +0.70710678j, 0.70710678+0.j ], # may vary [ 0.70710678+0.j , -0. +0.70710678j]])
Be careful about round-off error!
>>> a = np.array([[1 + 1e-9, 0], [0, 1 - 1e-9]]) >>> # Theor. eigenvalues are 1 +/- 1e-9 >>> eigenvalues, eigenvectors = LA.eig(a) >>> eigenvalues array([1., 1.]) >>> eigenvectors array([[1., 0.], [0., 1.]])