Navigation

When SciPy is built using the optimized ATLAS LAPACK and BLAS
libraries, it has very fast linear algebra capabilities. If you dig
deep enough, all of the raw lapack and blas libraries are available
for your use for even more speed. In this section, some easier-to-use
interfaces to these routines are described.

All of these linear algebra routines expect an object that can be
converted into a 2-dimensional array. The output of these routines is
also a two-dimensional array. There is a matrix class defined in
Numpy, which you can initialize with an appropriate Numpy array in
order to get objects for which multiplication is matrix-multiplication
instead of the default, element-by-element multiplication.

The matrix class is initialized with the SciPy command mat
which is just convenient short-hand for matrix. If you are going to be doing a lot of matrix-math, it
is convenient to convert arrays into matrices using this command. One
advantage of using the mat command is that you can enter
two-dimensional matrices using MATLAB-like syntax with commas or
spaces separating columns and semicolons separting rows as long as the
matrix is placed in a string passed to mat .

The inverse of a matrix is the matrix
such that where
is the identity matrix consisting of ones down the
main diagonal. Usually is denoted
. In SciPy, the matrix inverse of
the Numpy array, A, is obtained using linalg.inv(A) , or
using A.I if A is a Matrix. For example, let

Solving linear systems of equations is straightforward using the scipy
command linalg.solve. This command expects an input matrix and
a right-hand-side vector. The solution vector is then computed. An
option for entering a symmetrix matrix is offered which can speed up
the processing when applicable. As an example, suppose it is desired
to solve the following simultaneous equations:

We could find the solution vector using a matrix inverse:

However, it is better to use the linalg.solve command which can be
faster and more numerically stable. In this case it however gives the
same answer as shown in the following example:

The determinant of a square matrix is often denoted
and is a quantity often used in linear
algebra. Suppose are the elements of the matrix
and let
be the determinant of the matrix left by removing the
row and column from
. Then for any row

This is a recursive way to define the determinant where the base case
is defined by accepting that the determinant of a matrix is the only matrix element. In SciPy the determinant can be
calculated with linalg.det . For example, the determinant of

Matrix and vector norms can also be computed with SciPy. A wide range
of norm definitions are available using different parameters to the
order argument of linalg.norm . This function takes a rank-1
(vectors) or a rank-2 (matrices) array and an optional order argument
(default is 2). Based on these inputs a vector or matrix norm of the
requested order is computed.

For vector x , the order parameter can be any real number including
inf or -inf. The computed norm is

For matrix the only valid values for norm are inf, and ‘fro’ (or ‘f’) Thus,

Linear least-squares problems occur in many branches of applied
mathematics. In this problem a set of linear scaling coefficients is
sought that allow a model to fit data. In particular it is assumed
that data is related to data
through a set of coefficients and model functions
via the model

where represents uncertainty in the data. The
strategy of least squares is to pick the coefficients to
minimize

Theoretically, a global minimum will occur when

or

where

When is invertible, then

where is called the pseudo-inverse of
Notice that using this definition of
the model can be written

The command linalg.lstsq will solve the linear least squares
problem for given and
. In addition linalg.pinv or
linalg.pinv2 (uses a different method based on singular value
decomposition) will find given

The following example and figure demonstrate the use of
linalg.lstsq and linalg.pinv for solving a data-fitting
problem. The data shown below were generated using the model:

where for , ,
and Noise is added to and the
coefficients and are estimated using
linear least squares.

The generalized inverse is calculated using the command
linalg.pinv or linalg.pinv2. These two commands differ
in how they compute the generalized inverse. The first uses the
linalg.lstsq algorithm while the second uses singular value
decomposition. Let be an matrix,
then if the generalized inverse is

The eigenvalue-eigenvector problem is one of the most commonly
employed linear algebra operations. In one popular form, the
eigenvalue-eigenvector problem is to find for some square matrix
scalars and corresponding vectors
such that

For an matrix, there are (not necessarily
distinct) eigenvalues — roots of the (characteristic) polynomial

The eigenvectors, , are also sometimes called right
eigenvectors to distinguish them from another set of left eigenvectors
that satisfy

or

With it’s default optional arguments, the command linalg.eig
returns and However, it can also
return and just by itself (
linalg.eigvals returns just as well).

In addtion, linalg.eig can also solve the more general eigenvalue problem

for square matrices and The
standard eigenvalue problem is an example of the general eigenvalue
problem for When a generalized
eigenvalue problem can be solved, then it provides a decomposition of
as

where is the collection of eigenvectors into
columns and is a diagonal matrix of
eigenvalues.

By definition, eigenvectors are only defined up to a constant scale
factor. In SciPy, the scaling factor for the eigenvectors is chosen so
that

As an example, consider finding the eigenvalues and eigenvectors of
the matrix

The characteristic polynomial is

The roots of this polynomial are the eigenvalues of :

The eigenvectors corresponding to each eigenvalue can be found using
the original equation. The eigenvectors associated with these
eigenvalues can then be found.

Singular Value Decompostion (SVD) can be thought of as an extension of
the eigenvalue problem to matrices that are not square. Let
be an matrix with and
arbitrary. The matrices and
are square hermitian matrices [1] of
size and respectively. It is known
that the eigenvalues of square hermitian matrices are real and
non-negative. In addtion, there are at most
identical non-zero eigenvalues of
and
Define these positive eigenvalues as The
square-root of these are called singular values of
The eigenvectors of are collected by
columns into an unitary [2] matrix
while the eigenvectors of
are collected by columns in the
unitary matrix , the singular values are collected
in an zero matrix
with main diagonal entries set to
the singular values. Then

is the singular-value decomposition of Every
matrix has a singular value decomposition. Sometimes, the singular
values are called the spectrum of The command
linalg.svd will return ,
, and as an array of the
singular values. To obtain the matrix use
linalg.diagsvd. The following example illustrates the use of
linalg.svd .

where is an permutation matrix (a
permutation of the rows of the identity matrix), is
in lower triangular or trapezoidal matrix (
) with unit-diagonal, and
is an upper triangular or trapezoidal matrix. The
SciPy command for this decomposition is linalg.lu .

Such a decomposition is often useful for solving many simultaneous
equations where the left-hand-side does not change but the right hand
side does. For example, suppose we are going to solve

for many different . The LU decomposition allows this to be written as

Because is lower-triangular, the equation can be
solved for and finally
very rapidly using forward- and
back-substitution. An initial time spent factoring
allows for very rapid solution of similar systems of equations in the
future. If the intent for performing LU decomposition is for solving
linear systems then the command linalg.lu_factor should be used
followed by repeated applications of the command
linalg.lu_solve to solve the system for each new
right-hand-side.

Cholesky decomposition is a special case of LU decomposition
applicable to Hermitian positive definite matrices. When
and
for all ,
then decompositions of can be found so that

where is lower-triangular and is
upper triangular. Notice that The
command linagl.cholesky computes the cholesky
factorization. For using cholesky factorization to solve systems of
equations there are also linalg.cho_factor and
linalg.cho_solve routines that work similarly to their LU
decomposition counterparts.

For a square matrix, , the Schur
decomposition finds (not-necessarily unique) matrices
and such that

where is a unitary matrix and is
either upper-triangular or quasi-upper triangular depending on whether
or not a real schur form or complex schur form is requested. For a
real schur form both and are
real-valued when is real-valued. When
is a real-valued matrix the real schur form is only
quasi-upper triangular because blocks extrude from
the main diagonal corresponding to any complex- valued
eigenvalues. The command linalg.schur finds the Schur
decomposition while the command linalg.rsf2csf converts
and from a real Schur form to a
complex Schur form. The Schur form is especially useful in calculating
functions of matrices.

Finally, any arbitrary function that takes one complex number and
returns a complex number can be called as a matrix function using the
command linalg.funm. This command takes the matrix and an
arbitrary Python function. It then implements an algorithm from Golub
and Van Loan’s book “Matrix Computations “to compute function applied
to the matrix using a Schur decomposition. Note that the function
needs to accept complex numbers as input in order to work with this
algorithm. For example the following code computes the zeroth-order
Bessel function applied to a matrix.