scico.numpy.linalg

Linear algebra functions.

Functions

cholesky(a)	Cholesky decomposition.
cond(x[, p])	Compute the condition number of a matrix.
eig(a)	Compute the eigenvalues and right eigenvectors of a square array.
eigh(a[, UPLO, symmetrize_input])	Return the eigenvalues and eigenvectors of a complex Hermitian
eigvals(a)	Compute the eigenvalues of a general matrix.
eigvalsh(a[, UPLO])	Compute the eigenvalues of a complex Hermitian or real symmetric matrix.
inv(a)	Compute the (multiplicative) inverse of a matrix.
lstsq(a, b[, rcond, numpy_resid])	Return the least-squares solution to a linear matrix equation.
matrix_power(a, n)	Raise a square matrix to the (integer) power n.
matrix_rank(M[, tol])	Return matrix rank of array using SVD method
multi_dot(arrays, *[, precision])	Compute the dot product of two or more arrays in a single function call,
norm(x[, ord, axis, keepdims])	Matrix or vector norm.
qr(a[, mode])	Compute the qr factorization of a matrix.
slogdet(a, *[, method])	Compute the sign and (natural) logarithm of the determinant of an array.
solve(a, b)	Solve a linear matrix equation, or system of linear scalar equations.
svd(a[, full_matrices, compute_uv, hermitian])	Singular Value Decomposition.
tensorinv(a[, ind])	Compute the 'inverse' of an N-dimensional array.
tensorsolve(a, b[, axes])	Solve the tensor equation a x = b for x.

scico.numpy.linalg.cholesky(a)

Cholesky decomposition.

LAX-backend implementation of numpy.linalg.cholesky.

Original docstring below.

Return the Cholesky decomposition, L * L.H, of the square matrix a, where L is lower-triangular and .H is the conjugate transpose operator (which is the ordinary transpose if a is real-valued). a must be Hermitian (symmetric if real-valued) and positive-definite. No checking is performed to verify whether a is Hermitian or not. In addition, only the lower-triangular and diagonal elements of a are used. Only L is actually returned.

Parameters:

a ((..., M, M) array_like) – Hermitian (symmetric if all elements are real), positive-definite input matrix.

Return type:

Array

Returns:

L ((…, M, M) array_like) – Lower-triangular Cholesky factor of a. Returns a matrix object if a is a matrix object.

scico.numpy.linalg.cond(x, p=None)

Compute the condition number of a matrix.

LAX-backend implementation of numpy.linalg.cond.

Original docstring below.

This function is capable of returning the condition number using one of seven different norms, depending on the value of p (see Parameters below).

Parameters:

• x ((..., M, N) array_like) – The matrix whose condition number is sought.

• p ({None, 1, -1, 2, -2, inf, -inf, 'fro'}, optional) –

  Order of the norm used in the condition number computation:

  p	norm for matrices
  None	2-norm, computed directly using the SVD
  ’fro’	Frobenius norm
  inf	max(sum(abs(x), axis=1))
  -inf	min(sum(abs(x), axis=1))
  1	max(sum(abs(x), axis=0))
  -1	min(sum(abs(x), axis=0))
  2	2-norm (largest sing. value)
  -2	smallest singular value

  inf means the numpy.inf object, and the Frobenius norm is the root-of-sum-of-squares norm.

Returns:

c ({float, inf}) – The condition number of the matrix. May be infinite.

scico.numpy.linalg.eig(a)

Compute the eigenvalues and right eigenvectors of a square array.

LAX-backend implementation of numpy.linalg.eig.

This differs from numpy.linalg.eig in that the return type of jax.numpy.linalg.eig is always complex64 for 32-bit input, and complex128 for 64-bit input.

At present, non-symmetric eigendecomposition is only implemented on the CPU backend. However eigendecomposition for symmetric/Hermitian matrices is implemented more widely (see jax.numpy.linalg.eigh).

Original docstring below.

Parameters:

a ((..., M, M) array) – Matrices for which the eigenvalues and right eigenvectors will be computed

Return type:

Tuple[Array, Array]

Returns:

• w ((…, 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

• v ((…, M, M) array) – The normalized (unit “length”) eigenvectors, such that the column v[:,i] is the eigenvector corresponding to the eigenvalue w[i].

scico.numpy.linalg.eigh(a, UPLO=None, symmetrize_input=True)

Return the eigenvalues and eigenvectors of a complex Hermitian

LAX-backend implementation of numpy.linalg.eigh.

Original docstring below.

(conjugate symmetric) or a real symmetric matrix.

Returns two objects, a 1-D array containing the eigenvalues of a, and a 2-D square array or matrix (depending on the input type) of the corresponding eigenvectors (in columns).

Parameters:

• a ((..., M, M) array) – Hermitian or real 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 a (‘L’, default) or the upper triangular part (‘U’). Irrespective of this value only the real parts of the diagonal will be considered in the computation to preserve the notion of a Hermitian matrix. It therefore follows that the imaginary part of the diagonal will always be treated as zero.

Return type:

Tuple[Array, Array]

Returns:

• w ((…, M) ndarray) – The eigenvalues in ascending order, each repeated according to its multiplicity.

• v ({(…, M, M) ndarray, (…, M, M) matrix}) – The column v[:, i] is the normalized eigenvector corresponding to the eigenvalue w[i]. Will return a matrix object if a is a matrix object.

scico.numpy.linalg.eigvals(a)

Compute the eigenvalues of a general matrix.

LAX-backend implementation of numpy.linalg.eigvals.

Original docstring below.

Main difference between eigvals and eig: the eigenvectors aren’t returned.

Parameters:

a ((..., M, M) array_like) – A complex- or real-valued matrix whose eigenvalues will be computed.

Return type:

Array

Returns:

w ((…, M,) ndarray) – The eigenvalues, each repeated according to its multiplicity. They are not necessarily ordered, nor are they necessarily real for real matrices.

scico.numpy.linalg.eigvalsh(a, UPLO='L')

Compute the eigenvalues of a complex Hermitian or real symmetric matrix.

LAX-backend implementation of numpy.linalg.eigvalsh.

Original docstring below.

Main difference from eigh: the eigenvectors are not computed.

Parameters:

• a ((..., M, M) array_like) – A complex- or real-valued matrix whose eigenvalues are to be computed.

• UPLO ({'L', 'U'}, optional) – Specifies whether the calculation is done with the lower triangular part of a (‘L’, default) or the upper triangular part (‘U’). Irrespective of this value only the real parts of the diagonal will be considered in the computation to preserve the notion of a Hermitian matrix. It therefore follows that the imaginary part of the diagonal will always be treated as zero.

Return type:

Array

Returns:

w ((…, M,) ndarray) – The eigenvalues in ascending order, each repeated according to its multiplicity.

scico.numpy.linalg.inv(a)

Compute the (multiplicative) inverse of a matrix.

LAX-backend implementation of numpy.linalg.inv.

Original docstring below.

Given a square matrix a, return the matrix ainv satisfying dot(a, ainv) = dot(ainv, a) = eye(a.shape[0]).

Parameters:

a ((..., M, M) array_like) – Matrix to be inverted.

Return type:

Array

Returns:

ainv ((…, M, M) ndarray or matrix) – (Multiplicative) inverse of the matrix a.

scico.numpy.linalg.lstsq(a, b, rcond=None, *, numpy_resid=False)

Return the least-squares solution to a linear matrix equation.

LAX-backend implementation of numpy.linalg.lstsq.

It has two important differences:

1. In numpy.linalg.lstsq, the default rcond is -1, and warns that in the future the default will be None. Here, the default rcond is None.

2. In np.linalg.lstsq the returned residuals are empty for low-rank or over-determined solutions. Here, the residuals are returned in all cases, to make the function compatible with jit. The non-jit compatible numpy behavior can be recovered by passing numpy_resid=True.

The lstsq function does not currently have a custom JVP rule, so the gradient is poorly behaved for some inputs, particularly for low-rank a.

Original docstring below.

Computes the vector x that approximately solves the equation a @ x = b. The equation may be under-, well-, or over-determined (i.e., the number of linearly independent rows of a can be less than, equal to, or greater than its number of linearly independent columns). If a is square and of full rank, then x (but for round-off error) is the “exact” solution of the equation. Else, x minimizes the Euclidean 2-norm \(||b - ax||\). If there are multiple minimizing solutions, the one with the smallest 2-norm \(||x||\) is returned.

Parameters:

• a ((M, N) array_like) – “Coefficient” matrix.

• b ({(M,), (M, K)} array_like) – Ordinate or “dependent variable” values. If b is two-dimensional, the least-squares solution is calculated for each of the K columns of b.

• rcond (float, optional) –

  Cut-off ratio for small singular values of a. For the purposes of rank determination, singular values are treated as zero if they are smaller than rcond times the largest singular value of a.

  Changed in version 1.14.0: If not set, a FutureWarning is given. The previous default of -1 will use the machine precision as rcond parameter, the new default will use the machine precision times max(M, N). To silence the warning and use the new default, use rcond=None, to keep using the old behavior, use rcond=-1.

Return type:

Tuple[Array, Array, Array, Array]

Returns:

• x ({(N,), (N, K)} ndarray) – Least-squares solution. If b is two-dimensional, the solutions are in the K columns of x.

• residuals ({(1,), (K,), (0,)} ndarray) – Sums of squared residuals: Squared Euclidean 2-norm for each column in b - a @ x. If the rank of a is < N or M <= N, this is an empty array. If b is 1-dimensional, this is a (1,) shape array. Otherwise the shape is (K,).

• rank (int) – Rank of matrix a.

• s ((min(M, N),) ndarray) – Singular values of a.

scico.numpy.linalg.matrix_power(a, n)

Raise a square matrix to the (integer) power n.

LAX-backend implementation of numpy.linalg.matrix_power.

Original docstring below.

For positive integers n, the power is computed by repeated matrix squarings and matrix multiplications. If n == 0, the identity matrix of the same shape as M is returned. If n < 0, the inverse is computed and then raised to the abs(n).

Note

Stacks of object matrices are not currently supported.

Parameters:

• a ((..., M, M) array_like) – Matrix to be “powered”.

• n (int) – The exponent can be any integer or long integer, positive, negative, or zero.

Return type:

Array

Returns:

a**n ((…, M, M) ndarray or matrix object) – The return value is the same shape and type as M; if the exponent is positive or zero then the type of the elements is the same as those of M. If the exponent is negative the elements are floating-point.

scico.numpy.linalg.matrix_rank(M, tol=None)

Return matrix rank of array using SVD method

LAX-backend implementation of numpy.linalg.matrix_rank.

Original docstring below.

Rank of the array is the number of singular values of the array that are greater than tol.

Changed in version 1.14: Can now operate on stacks of matrices

Parameters:

tol ((...) array_like, float, optional) –

Threshold below which SVD values are considered zero. If tol is None, and S is an array with singular values for M, and eps is the epsilon value for datatype of S, then tol is set to S.max() * max(M, N) * eps.

Changed in version 1.14: Broadcasted against the stack of matrices

Return type:

Array

Returns:

rank ((…) array_like) – Rank of A.

scico.numpy.linalg.multi_dot(arrays, *, precision=None)

Compute the dot product of two or more arrays in a single function call,

LAX-backend implementation of numpy.linalg.multi_dot.

Original docstring below.

while automatically selecting the fastest evaluation order.

multi_dot chains numpy.dot and uses optimal parenthesization of the matrices. Depending on the shapes of the matrices, this can speed up the multiplication a lot.

If the first argument is 1-D it is treated as a row vector. If the last argument is 1-D it is treated as a column vector. The other arguments must be 2-D.

Think of multi_dot as:

def multi_dot(arrays): return functools.reduce(np.dot, arrays)

Parameters:

arrays (sequence of array_like) – If the first argument is 1-D it is treated as row vector. If the last argument is 1-D it is treated as column vector. The other arguments must be 2-D.

Returns:

output (ndarray) – Returns the dot product of the supplied arrays.

scico.numpy.linalg.norm(x, ord=None, axis=None, keepdims=False)

Matrix or vector norm.

LAX-backend implementation of numpy.linalg.norm.

Original docstring below.

This function is able to return one of eight different matrix norms, or one of an infinite number of vector norms (described below), depending on the value of the ord parameter.

Parameters:

• x (array_like) – Input array. If axis is None, x must be 1-D or 2-D, unless ord is None. If both axis and ord are None, the 2-norm of x.ravel will be returned.

• ord ({non-zero int, inf, -inf, 'fro', 'nuc'}, optional) – Order of the norm (see table under Notes). inf means numpy’s inf object. The default is None.

• axis ({None, int, 2-tuple of ints}, optional.) – If axis is an integer, it specifies the axis of x along which to compute the vector norms. If axis is a 2-tuple, it specifies the axes that hold 2-D matrices, and the matrix norms of these matrices are computed. If axis is None then either a vector norm (when x is 1-D) or a matrix norm (when x is 2-D) is returned. The default is None.

• keepdims (bool, optional) – If this is set to True, the axes which are normed over are left in the result as dimensions with size one. With this option the result will broadcast correctly against the original x.

Return type:

Array

Returns:

n (float or ndarray) – Norm of the matrix or vector(s).

scico.numpy.linalg.qr(a, mode='reduced')

Compute the qr factorization of a matrix.

LAX-backend implementation of numpy.linalg.qr.

Original docstring below.

Factor the matrix a as qr, where q is orthonormal and r is upper-triangular.

Parameters:

• a (array_like, shape (..., M, N)) – An array-like object with the dimensionality of at least 2.

• mode ({'reduced', 'complete', 'r', 'raw'}, optional) –

  If K = min(M, N), then

  • ’reduced’returns q, r with dimensions

    (…, M, K), (…, K, N) (default)

  • ’complete’ : returns q, r with dimensions (…, M, M), (…, M, N)

  • ’r’ : returns r only with dimensions (…, K, N)

  • ’raw’ : returns h, tau with dimensions (…, N, M), (…, K,)

  The options ‘reduced’, ‘complete, and ‘raw’ are new in numpy 1.8, see the notes for more information. The default is ‘reduced’, and to maintain backward compatibility with earlier versions of numpy both it and the old default ‘full’ can be omitted. Note that array h returned in ‘raw’ mode is transposed for calling Fortran. The ‘economic’ mode is deprecated. The modes ‘full’ and ‘economic’ may be passed using only the first letter for backwards compatibility, but all others must be spelled out. See the Notes for more explanation.

Return type:

Union[Array, Tuple[Array, Array]]

Returns:

• q (ndarray of float or complex, optional) – A matrix with orthonormal columns. When mode = ‘complete’ the result is an orthogonal/unitary matrix depending on whether or not a is real/complex. The determinant may be either +/- 1 in that case. In case the number of dimensions in the input array is greater than 2 then a stack of the matrices with above properties is returned.

• r (ndarray of float or complex, optional) – The upper-triangular matrix or a stack of upper-triangular matrices if the number of dimensions in the input array is greater than 2.

• (h, tau) (ndarrays of np.double or np.cdouble, optional) – The array h contains the Householder reflectors that generate q along with r. The tau array contains scaling factors for the reflectors. In the deprecated ‘economic’ mode only h is returned.

scico.numpy.linalg.slogdet(a, *, method=None)

Compute the sign and (natural) logarithm of the determinant of an array.

LAX-backend implementation of numpy.linalg.slogdet.

Original docstring below.

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_like) – Input array, has to be a square 2-D array.

Return type:

Tuple[Array, Array]

Returns:

• sign ((…) array_like) – A number representing the sign of the determinant. For a real matrix, this is 1, 0, or -1. For a complex matrix, this is a complex number with absolute value 1 (i.e., it is on the unit circle), or else 0.

• logdet ((…) array_like) – The natural log of the absolute value of the determinant.

• 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).

scico.numpy.linalg.solve(a, b)

Solve a linear matrix equation, or system of linear scalar equations.

LAX-backend implementation of numpy.linalg.solve.

Original docstring below.

Computes the “exact” solution, x, of the well-determined, i.e., full rank, linear matrix equation ax = b.

Parameters:

• a ((..., M, M) array_like) – Coefficient matrix.

• b ({(..., M,), (..., M, K)}, array_like) – Ordinate or “dependent variable” values.

Return type:

Array

Returns:

x ({(…, M,), (…, M, K)} ndarray) – Solution to the system a x = b. Returned shape is identical to b.

scico.numpy.linalg.svd(a, full_matrices=True, compute_uv=True, hermitian=False)

Singular Value Decomposition.

LAX-backend implementation of numpy.linalg.svd.

Original docstring below.

When a is a 2D array, and full_matrices=False, then it is factorized as u @ np.diag(s) @ vh = (u * s) @ vh, where u and the Hermitian transpose of vh are 2D arrays with orthonormal columns and s is a 1D array of a’s singular values. When a is higher-dimensional, SVD is applied in stacked mode as explained below.

Parameters:

• a ((..., M, N) array_like) – A real or complex array with a.ndim >= 2.

• full_matrices (bool, optional) – If True (default), u and vh 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 vh in addition to s. True by default.

• hermitian (bool, optional) – If True, a is assumed to be Hermitian (symmetric if real-valued), enabling a more efficient method for finding singular values. Defaults to False.

Return type:

Union[Array, Tuple[Array, Array, Array]]

Returns:

• u ({ (…, M, M), (…, M, K) } array) – Unitary array(s). The first a.ndim - 2 dimensions have the same size as those of the input a. The size of the last two dimensions depends on the value of full_matrices. Only returned when compute_uv is True.

• s ((…, K) array) – Vector(s) with the singular values, within each vector sorted in descending order. The first a.ndim - 2 dimensions have the same size as those of the input a.

• vh ({ (…, N, N), (…, K, N) } array) – Unitary array(s). The first a.ndim - 2 dimensions have the same size as those of the input a. The size of the last two dimensions depends on the value of full_matrices. Only returned when compute_uv is True.

scico.numpy.linalg.tensorinv(a, ind=2)

Compute the ‘inverse’ of an N-dimensional array.

LAX-backend implementation of numpy.linalg.tensorinv.

Original docstring below.

The result is an inverse for a relative to the tensordot operation tensordot(a, b, ind), i. e., up to floating-point accuracy, tensordot(tensorinv(a), a, ind) is the “identity” tensor for the tensordot operation.

Parameters:

• a (array_like) – Tensor to ‘invert’. Its shape must be ‘square’, i. e., prod(a.shape[:ind]) == prod(a.shape[ind:]).

• ind (int, optional) – Number of first indices that are involved in the inverse sum. Must be a positive integer, default is 2.

Returns:

b (ndarray) – a’s tensordot inverse, shape a.shape[ind:] + a.shape[:ind].

scico.numpy.linalg.tensorsolve(a, b, axes=None)

Solve the tensor equation a x = b for x.

LAX-backend implementation of numpy.linalg.tensorsolve.

Original docstring below.

It is assumed that all indices of x are summed over in the product, together with the rightmost indices of a, as is done in, for example, tensordot(a, x, axes=x.ndim).

Parameters:

• a (array_like) – Coefficient tensor, of shape b.shape + Q. Q, a tuple, equals the shape of that sub-tensor of a consisting of the appropriate number of its rightmost indices, and must be such that prod(Q) == prod(b.shape) (in which sense a is said to be ‘square’).

• b (array_like) – Right-hand tensor, which can be of any shape.

• axes (tuple of ints, optional) – Axes in a to reorder to the right, before inversion. If None (default), no reordering is done.

Returns:

x (ndarray, shape Q)