"""Linear algebra utils"""
import numpy as np
from numpy.linalg import svd
from ddg.abc import NonExact
# Computation of rank and nullspace of a matrix with given tolerances.
# source: http://scipy-cookbook.readthedocs.io/items/RankNullspace.html
[docs]@NonExact.nonexact_function
def rank(A, atol=1e-13, rtol=None):
"""Estimate the rank (i.e. the dimension of the nullspace) of a matrix.
The algorithm used by this function is based on the singular value
decomposition of `A`.
Parameters
----------
A : ndarray
A should be at most 2-D. A 1-D array with length n will be treated
as a 2-D with shape (1, n)
atol : float
The absolute tolerance for a zero singular value. Singular values
smaller than `atol` are considered to be zero.
rtol : float
The relative tolerance. Singular values less than rtol*smax are
considered to be zero, where smax is the largest singular value.
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.
Returns
-------
r : int
The estimated rank of the matrix.
See also
--------
numpy.linalg.matrix_rank
matrix_rank is basically the same as this function, but it does not
provide the option of the absolute tolerance.
Notes
-----
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See ddg.abc.NonExact for details.
"""
A = np.atleast_2d(A)
s = svd(A, compute_uv=False)
tol = max(atol, rtol * s[0])
rank = int((s >= tol).sum())
return rank
[docs]@NonExact.nonexact_function
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 : ndarray
A should be at most 2-D. A 1-D array with length k will be treated
as a 2-D with shape (1, k)
atol : float
The absolute tolerance for a zero singular value. Singular values
smaller than `atol` are considered to be zero.
rtol : float
The relative tolerance. Singular values less than rtol*smax are
considered to be zero, where smax is the largest singular value.
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.
Returns
-------
ns : ndarray
If `A` is an array with shape (m, k), then `ns` will be an array
with shape (k, n), where n is the estimated dimension of the
nullspace of `A`. The columns of `ns` are a basis for the
nullspace; each element in numpy.dot(A, ns) will be approximately
zero.
Notes
-----
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See ddg.abc.NonExact for details.
"""
A = np.atleast_2d(A)
u, s, vh = svd(A)
tol = max(atol, rtol * s[0])
nnz = (s >= tol).sum()
ns = vh[nnz:].conj().T
return ns
[docs]@NonExact.nonexact_function
def row_basis(A, atol=1e-13, rtol=None):
"""
Notes
-----
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See ddg.abc.NonExact for details.
"""
A = np.atleast_2d(A)
# Nullspace. Avoids index error
if A.size == 0:
return A
u, s, vh = svd(A)
tol = max(atol, rtol * s[0])
nnz = (s >= tol).sum()
ns = vh[:nnz]
return ns
[docs]@NonExact.nonexact_function
def col_basis(A, atol=1e-13, rtol=None):
"""
Notes
-----
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See ddg.abc.NonExact for details.
"""
return row_basis(A.T, atol, rtol).T
[docs]def arrays_to_matrix(arrays):
"""Converts a list of arrays to a matrix.
The arrays will be the columns (!) of the matrix.
Parameters
----------
arrays : list
List of lists, tuples or numpy.ndarray with equal
dimensions.
Returns
-------
matrix : numpy.ndarray
Matrix whose columns are the given arrays.
"""
return np.array(arrays).transpose()
[docs]def matrix_to_arrays(matrix):
"""
Converts the columns (!) of a (m x n)-matrix to a
list of arrays.
Parameters
----------
matrix : numpy.ndarray of shape (m, n)
Returns
-------
arrays : list
List of n numpy.ndarray, each of shape (m,)
"""
return list(matrix.transpose())
[docs]def e(i, N):
"""
i-th canonical basis vector of size N
"""
v = np.zeros(N)
v[i] = 1.
return v
[docs]def permutation_matrix(i,j, N):
"""
Projective transformation that permutes homogeneous coordinates i and j
"""
P = np.identity(N)
P[i,i] = P[j,j] = 0
P[i,j] = P[j,i] = 1
return P
[docs]def extend_to_basis(v, i=-1):
"""
Complete one vector into a basis.
Parameters
----------
v : ndarray
Vector you want to complete into a basis
i : int
Position you want v to be in the returned basis
Returns
-------
b : ndarray
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.
"""
N = len(v)
j = None
for k in range(N):
if v[k] != 0:
j = k
b = np.identity(N)
b[:, j] = v
b = np.dot(b, permutation_matrix(i, j, N))
return b
[docs]@NonExact.nonexact_function
def coordinates(vectors, basis, atol=1e-13, rtol=None):
"""Compute coordinates of vectors in a basis of whole space or subspace.
Parameters
----------
vectors : numpy.ndarray of shape (n,k) or (n,)
A matrix whose columns are vectors in the space
basis : numpy.ndarray of shape (n,m)
Matrix whose linearly independent columns are the basis of the space
Returns
-------
coordinates : numpy.ndarray of shape (m,k) or (m,)
coordinate matrix such that basis @ coordinates == vectors.
Raises
------
ValueError
If a vector does not lie in the subspace.
"""
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 np.allclose(residuals, 0, atol=atol, rtol=rtol):
raise ValueError("A vector does not lie in the subspace.")
return coordinates
[docs]@NonExact.nonexact_function
def linear_dependence(points, atol=1e-4, rtol=None):
"""
Computes the coefficients of a linear dependence of points.
Parameters
----------
points : list
List of vectors
Returns
-------
coefficients of linear dependence of points
Notes
-----
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See ddg.abc.NonExact for details.
"""
u, s, vh = np.linalg.svd(points.T)
tol = max(atol, rtol * s[0])
nnz = (s >= tol).sum()
ns = vh[nnz:].conj()[0]
return ns