"""Linear algebra utils"""
from itertools import combinations, product
import numpy as np
from numpy.linalg import svd
from ddg import nonexact
# Computation of rank and nullspace of a matrix with given tolerances.
# source: http://scipy-cookbook.readthedocs.io/items/RankNullspace.html
[docs]def rank(A, atol=1e-13, rtol=None):
"""Estimate the rank of a matrix, i.e. the dimension of its image.
The algorithm used by this function is based on the singular value
decomposition of `A`.
Parameters
----------
A : numpy.ndarray of shape (m, n)
A 1-D array with length n will be treated as a 2-D array with
shape (1, n).
atol : float, (default=1e-13)
The absolute tolerance for a zero singular value. Singular values
smaller than `atol` are considered to be zero.
rtol : float, (default=None)
The relative tolerance. Singular values less than ``rtol*smax`` are
considered to be zero, where ``smax`` is the largest singular value.
Returns
-------
r : int
The estimated rank of the matrix.
See also
--------
numpy.linalg.matrix_rank
NumPy's ``matrix_rank`` is basically the same as this function,
but it does not provide the option of the absolute tolerance.
Notes
-----
If both `atol` and `rtol` are positive, the combined tolerance is
the maximum of the two, that is::
tol = max(atol, rtol * smax)
Singular values smaller than `tol` are considered to be zero.
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See :py:mod:`ddg.nonexact` for details.
"""
A = np.atleast_2d(A)
s = svd(A, compute_uv=False)
rank_ = nonexact.svd_rank(s, atol=atol, rtol=rtol)
return rank_
[docs]def nullspace(A, atol=1e-13, rtol=None):
"""Compute an approximate basis for the nullspace of A.
The algorithm used by this function is based on the singular value
decomposition of `A`.
Parameters
----------
A : numpy.ndarray of shape (m, n)
A 1-D array with length n will be treated as a 2-D array with
shape (1, n).
atol : float, (default=1e-13)
The absolute tolerance for a zero singular value. Singular values
smaller than `atol` are considered to be zero.
rtol : float, (default=None)
The relative tolerance. Singular values less than ``rtol*smax`` are
considered to be zero, where ``smax`` is the largest singular value.
Returns
-------
ns : numpy.ndarray
If `A` is an array with shape (m, n), then `ns` will be an array
with shape (n, k), where k is the estimated dimension of the
nullspace of `A`. The columns of `ns` form a basis for the
nullspace; each element in ``numpy.dot(A, ns)`` will be approximately
zero.
Notes
-----
If both `atol` and `rtol` are positive, the combined tolerance is
the maximum of the two, that is::
tol = max(atol, rtol * smax)
Singular values smaller than `tol` are considered to be zero.
This function uses the global tolerance defaults if `atol` or `rtol`
are set to None. See :py:mod:`ddg.nonexact` for details.
"""
A = np.atleast_2d(A)
*_, s, vh = svd(A)
nnz = nonexact.svd_rank(s, atol=atol, rtol=rtol)
ns = vh[nnz:].conj().T
return ns
[docs]def minors(A, k):
"""Compute all k-minors of a given (m, n)-matrix.
A k-minor is the determinant of a (k, k)-submatrix.
Parameters
----------
A : numpy.ndarray of shape (m, n)
A 1-D array with length n will be treated as a 2-D array with
shape (1, n).
k : int
Order of minors.
Returns
-------
minors : list of length d
The resulting minors where length `d` is given by::
d = (m choose k)*(n choose k)
Raises
------
ValueError
If the order of minors `k` is greater than ``min(m, n)`` for a
(m, n)-matrix.
"""
def determinant(rows_and_columns):
rows, cols = rows_and_columns
return np.linalg.det(A[np.ix_(rows, cols)])
m = A.shape[0]
n = A.shape[1]
if k > min(m, n):
raise ValueError(
f"Can not return a minor of order {k} for matrix with shape {(m, n)}"
)
all_rows = range(m)
all_cols = range(n)
m = map(determinant, product(combinations(all_rows, k), combinations(all_cols, k)))
return list(m)
[docs]def row_basis(A, atol=1e-13, rtol=None):
"""
Computes an approximate basis for the span of rows of a given matrix.
The algorithm used by this function is based on the singular value
decomposition of `A`.
Parameters
----------
A : numpy.ndarray of shape (m, n)
A 1-D array with length n will be treated as a 2-D array with
shape (1, n).
atol : float, (default=1e-13)
The absolute tolerance for a zero singular value. Singular values
smaller than `atol` are considered to be zero.
rtol : float, (default=None)
The relative tolerance. Singular values less than ``rtol*smax`` are
considered to be zero, where ``smax`` is the largest singular value.
Returns
-------
ns : numpy.ndarray of shape (k, n)
If `A` is an array with shape (m, n), then `ns` will be an array
with shape (k, n), where k is the estimated dimension of the
row space of `A`. The rows of `ns` form its basis; they will be
pairwise orthogonal and have unit length.
Notes
-----
If both `atol` and `rtol` are positive, the combined tolerance is
the maximum of the two, that is::
tol = max(atol, rtol * smax)
Singular values smaller than `tol` are considered to be zero.
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See :py:mod:`ddg.nonexact` for details.
"""
A = np.atleast_2d(A)
# Nullspace. Avoids index error
if A.size == 0:
return A
*_, s, vh = svd(A)
nnz = nonexact.svd_rank(s, atol=atol, rtol=rtol)
ns = vh[:nnz]
return ns
[docs]def col_basis(A, atol=1e-13, rtol=None):
"""
Computes an approximate basis for the span of columns of a given matrix.
The algorithm used by this function is based on the singular value
decomposition of `A`.
Parameters
----------
A : numpy.ndarray of shape (m, n)
A 1-D array with length n will be treated as a 2-D array with
shape (1, n).
atol : float, (default=1e-13)
The absolute tolerance for a zero singular value. Singular values
smaller than `atol` are considered to be zero.
rtol : float, (default=None)
The relative tolerance. Singular values less than ``rtol*smax`` are
considered to be zero, where ``smax`` is the largest singular value.
Returns
-------
ns : numpy.ndarray of shape (m, k)
If `A` is an array with shape (m, n), then `ns` will be an array
with shape (m, k), where k is the estimated dimension of the
column space of `A`. The columns of `ns` form its basis; they will
be pairwise orthogonal and have unit length.
Notes
-----
If both `atol` and `rtol` are positive, the combined tolerance is
the maximum of the two, that is::
tol = max(atol, rtol * smax)
Singular values smaller than `tol` are considered to be zero.
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See :py:mod:`ddg.nonexact` for details.
"""
return row_basis(A.T, atol, rtol).T
[docs]def matrix_to_arrays(matrix):
"""
Converts the columns(!) of a (m, n)-matrix to a
list of arrays.
Parameters
----------
matrix : numpy.ndarray of shape (m, n)
Returns
-------
list of length n of numpy.ndarray, each of shape (m,)
"""
return list(matrix.transpose())
[docs]def e(i, n):
"""
Returns the ith basis vector in the canonical basis of K^n
over K.
Parameters
----------
i : int
Determines which canonical vector is returned,
0 <= i < n.
n : int
Determines the size of canonical vector.
Returns
-------
numpy.ndarray of shape (n,)
"""
v = np.zeros(n)
v[i] = 1.0
return v
[docs]def extend_to_basis(v, i=-1):
"""
Complete a vector in R^n to a basis for R^n.
Parameters
----------
v : numpy.ndarray of shape (n,)
Vector to complete into a basis.
i : int, (default=-1)
Determines position of v in the returned basis.
Raises
------
ValueError
If the given vector `v` is the zero vector.
Returns
-------
B : numpy.ndarray of shape (n, n)
If `v` is a vector of size n, then `B` will be a matrix
of shape (n, n) where the columns are the basis vectors
and ``B[:, i] == v``.
"""
if not np.any(v):
raise ValueError("Zero vector can not be extended to a basis.")
n = len(v)
j = None
for k in range(n):
if v[k] != 0:
j = k
B = np.identity(n)
B[:, j] = v
B[:, [i, j]] = B[:, [j, i]]
return B
[docs]def coordinates(vectors, basis, atol=1e-13, rtol=None):
"""Compute the coordinates of vectors in a given basis.
The basis can be of a certain subspace or of the whole space.
Parameters
----------
vectors : numpy.ndarray of shape (n, k) or (n,)
Matrix whose columns are the given vectors.
basis : numpy.ndarray of shape (n, m)
Matrix whose linearly independent columns form a basis.
atol : float, (default=1e-13)
The absolute tolerance parameter. See numpy.allclose.
rtol : float, (default=None)
The relative tolerance parameter. See numpy.allclose.
Returns
-------
coordinates : numpy.ndarray of shape (m, k) or (m,)
Coordinate matrix such that ``basis @ coordinates == vectors``.
Raises
------
ValueError
If at least one vector does not lie in the subspace defined by
given basis.
Notes
-----
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See :py:mod:`ddg.nonexact` for details.
"""
if np.shape(basis)[0] == np.shape(basis)[1]:
return np.linalg.solve(basis, vectors)
coordinates, residuals, _, _ = np.linalg.lstsq(basis, vectors, rcond=None)
if not nonexact.allclose(residuals, 0, atol=atol, rtol=rtol):
raise ValueError(
"Input vector does not lie in the subspace defined by given basis."
)
return coordinates
[docs]def linear_dependence(points, atol=1e-4, rtol=None):
"""
Computes coefficients of a linear dependence of points, i.e
real coefficients cᵢ such that::
∑ᵢ cᵢ * pᵢ = 0
The algorithm used by this function is based on the singular value
decomposition of given matrix of points.
Parameters
----------
points : numpy.ndarray of shape (n, m)
Matrix with the given n points as rows.
atol : float, (default=1e-4)
The absolute tolerance for a zero singular value. Singular values
smaller than `atol` are considered to be zero.
rtol : float, (default=None)
The relative tolerance. Singular values less than ``rtol*smax`` are
considered to be zero, where ``smax`` is the largest singular value.
Returns
-------
numpy.ndarray of shape (n,)
A 1-D array of coefficients.
Notes
-----
If both `atol` and `rtol` are positive, the combined tolerance is
the maximum of the two; that is::
tol = max(atol, rtol * smax)
Singular values smaller than `tol` are considered to be zero.
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See :py:mod:`ddg.nonexact` for details.
"""
n = points.shape[0]
*_, s, vh = np.linalg.svd(points.T)
nnz = nonexact.svd_rank(s, atol=atol, rtol=rtol)
if nnz == vh.shape[0]:
ns = np.zeros(n)
else:
ns = vh[nnz:].conj()[0]
return ns