"""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.
Examples
--------
>>> import numpy as np
>>> from ddg.math.linalg import rank
>>> A = np.array([[1, 2], [2, 5]])
>>> rank(A)
2
>>> B = np.array([[1, 2], [2, 4]])
>>> rank(B)
1
"""
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.
Examples
--------
>>> import numpy as np
>>> from ddg.math.linalg import nullspace
>>> A = np.array([[1, 2], [2, 5]])
>>> nullspace(A)
array([], shape=(2, 0), dtype=float64)
>>> B = np.array([[1, 2], [2, 4]])
>>> nullspace(B)
array([[-0.89442719],
[ 0.4472136 ]])
"""
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.
Examples
--------
>>> import numpy as np
>>> from ddg.math.linalg import minors
>>> B = np.array([[1, 2], [2, 4]])
>>> minors(B, 1)
[1.0, 2.0, 2.0, 4.0]
>>> minors(B, 2)
[0.0]
>>>
"""
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.
Examples
--------
>>> import numpy as np
>>> from ddg.math.linalg import row_basis
>>> A = np.array([[1, 2], [2, 5]])
>>> row_basis(A)
array([[-0.38268343, -0.92387953],
[-0.92387953, 0.38268343]])
>>> B = np.array([[1, 2], [2, 4]])
>>> row_basis(B)
array([[-0.4472136 , -0.89442719]])
"""
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.
Examples
--------
>>> import numpy as np
>>> from ddg.math.linalg import col_basis
>>> A = np.array([[1, 2], [2, 5]])
>>> col_basis(A)
array([[-0.38268343, -0.92387953],
[-0.92387953, 0.38268343]])
>>> B = np.array([[1, 2], [2, 4]])
>>> col_basis(B)
array([[-0.4472136 ],
[-0.89442719]])
"""
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,)
Examples
--------
>>> import numpy as np
>>> from ddg.math.linalg import matrix_to_arrays
>>> A = np.array([[1, 2], [2, 5]])
>>> matrix_to_arrays(A)
[array([1, 2]), array([2, 5])]
>>> type(matrix_to_arrays(A))
<class 'list'>
"""
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,)
Examples
--------
>>> import numpy as np
>>> from ddg.math.linalg import e
>>> e(0, 3)
array([1., 0., 0.])
>>> e(2, 3)
array([0., 0., 1.])
>>> e(2, 5)
array([0., 0., 1., 0., 0.])
"""
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 : array_like 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``.
Examples
--------
>>> import numpy as np
>>> from ddg.math.linalg import extend_to_basis
>>> v = np.array([1, 0, 0])
>>> extend_to_basis(v)
array([[0., 0., 1.],
[0., 1., 0.],
[1., 0., 0.]])
>>> w = (1, 2, 3)
>>> extend_to_basis(w)
array([[1., 0., 1.],
[0., 1., 2.],
[0., 0., 3.]])
"""
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.
Examples
--------
>>> import numpy as np
>>> from ddg.math.linalg import coordinates
>>> B = np.array([[-1, 0, 0], [0, 2, 0], [0, 0, -2]])
>>> V = np.array([[1, 2, 3], [5, 1, 4], [0, 0, 1]])
>>> V
array([[1, 2, 3],
[5, 1, 4],
[0, 0, 1]])
>>> coordinates(V, B)
array([[-1. , -2. , -3. ],
[ 2.5, 0.5, 2. ],
[-0. , -0. , -0.5]])
"""
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.
Examples
--------
>>> import numpy as np
>>> from ddg.math.linalg import linear_dependence
>>> B = np.array([[1, 2], [2, 4]])
>>> linear_dependence(B)
array([-0.89442719, 0.4472136 ])
>>> A = np.array([[1, 2], [2, 5]])
>>> linear_dependence(A)
array([0., 0.])
"""
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
[docs]def laplace_determinant(M):
"""
Calculates the determinant for a 4x4 matrix M via Laplace expansion
along the first row.
The determinant is calculated with the formula::
4
det(M)= Σ m_1j det(M_1j)
j=1
where M_1j is the matrix M with first row and j-th column deleted.
Then det(M_1j) can be calculated as signed sum of the product in its
diagonals::
det((a_ij)_i,j=1,2,3)= a_00 a_11 a_22
+ a_01 a_12 a_20
+ a_02 a_10 a_21
- a_02 a_11 a_20
- a_01 a_10 a_22
- a_00 a_12 a_21
Then the result will be a polynomial function instead of a rational
function as opposed to the Gauss algorithm where there might be
division through entries
Parameters
----------
M : numpy.ndarray
of size 4 x 4
Returns
-------
deter : determinant of M as a float
"""
if not isinstance(M, np.ndarray) or M.shape != (4, 4):
raise NotImplementedError(
"Laplace determinant only implemented for a 4x4 np.ndarray"
)
# start with determinant 0, then add to it
deter = 0
# choose pivot in first row, i-th place
for i in range(4):
# make a copy of the original matrix
A = M.copy()
# delete the first row of the copy
A = np.delete(A, (0), axis=0)
# delete the corresponding column of the copy too
A = np.delete(A, (i), axis=1)
# we now have a 3x3 Matrix and can calculate the according summand
# by the 3x3 determinant formula as follows:
deter += (
(-1) ** i
* M[0][i]
* (
A[0][0] * A[1][1] * A[2][2]
+ A[0][1] * A[1][2] * A[2][0]
+ A[0][2] * A[1][0] * A[2][1]
- (
A[0][2] * A[1][1] * A[2][0]
+ A[0][1] * A[1][0] * A[2][2]
+ A[0][0] * A[1][2] * A[2][1]
)
)
)
return deter