"""Utility module for inner products.
An inner product is a real (or complex) valued bilinear function over some
vector space with the signature inner_product(v, w).
Note that we do not require inner products to be positive-definite.
Some of the utility functions require the given inner product to be able to
handle matrices for its arguments v and w, i.e. for two matrices V, W of shape
(n,m) and (n,l), and with columns v_i and w_j respectively, inner_product(V,W)
would return a matrix VW of shape (m,l) with VW[i][j] = inner_product(v_i, w_j).
We call inner products which satisfy this condition vectorized.
Others do not strictly require this, but will assume this property, if a matrix is
passed instead of a vector for one of its arguments. Refer for the Notes sections
of the utility functions for more.
"""
import numpy as np
import numpy.linalg as la
import ddg.math.linalg as linalg
from ddg.nonexact import get_tol_defaults
[docs]def col_complement(A, inner_product=None):
"""Returns the orthogonal complement of the columns of a given
array with respect to the given inner product.
Parameters
----------
A : numpy.ndarray of shape (n,m)
Array for which to return the column complement
inner_product : Callable (optional, default=None)
Inner product to be used. When the argument is None, the standard
inner product will be used.
Returns
-------
numpy.ndarray of shape (n, k)
Column complement of the given array
Notes
-----
If the columns of the given matrix A span the whole space, a vector of
shape (n, 0) is returned
See also
--------
col_complement_oriented, row_complement, row_complement_oriented
"""
if inner_product is None:
def inner_product(v, w):
return np.array(v).T @ np.array(w)
B = linalg.col_basis(A)
# Return whole space. Avoids index error
if (B == 0).all(): # 0 matrix or empty matrix
return np.eye(np.shape(B)[0])
AL = np.dot(get_matrix(inner_product, len(B[:, 0])), B)
Q, _ = la.qr(AL, mode="complete")
ACN = orthonormalize(Q[:, len(B[0]) :], inner_product)
return ACN # A complement normalized
[docs]def col_complement_oriented(A, inner_product=None):
"""Returns the oriented orthogonal complement of the columns of a given
array with respect to the given inner product.
Parameters
----------
A : numpy.ndarray of shape (n, m)
Array for which to return the column complement
inner_product : Callable (optional, default=None)
Inner product to be used. When the argument is None, the standard
inner product will be used.
Returns
-------
numpy.ndarray of shape (n, k)
Column complement of the given array
Notes
-----
If the columns of the given matrix A span the whole space, a vector of
shape (n, 0) is returned
See also
--------
col_complement, row_complement, row_complement_oriented
"""
ACN = col_complement(A, inner_product)
mat = np.concatenate((A, ACN), axis=1)
if (
ACN.shape[1] != 0
and mat.shape[0] == mat.shape[1]
and np.linalg.slogdet(mat)[0] < 0
):
ACN.T[-1] = -ACN.T[-1]
return ACN # A complement normalized
[docs]def row_complement(A, inner_product=None):
"""Returns the orthogonal complement of the rows of a given
array with respect to the given inner product.
Parameters
----------
A : numpy.ndarray of shape (m, n)
Array for which to return the column complement
inner_product : Callable (optional, default=None)
Inner product to be used. When the argument is None, the standard
inner product will be used.
Returns
-------
numpy.ndarray
Column complement of the given array
Notes
-----
If the rows of the given matrix A span the whole space, a vector of
shape (0, n) is returned
See also
--------
col_complement, col_complement_oriented, row_complement_oriented
"""
return col_complement(A.T, inner_product).T
[docs]def row_complement_oriented(A, inner_product=None):
"""Returns the oriented orthogonal complement of the rows of a given
array with respect to the given inner product.
Parameters
----------
A : numpy.ndarray of shape (m, n)
Array for which to return the column complement
inner_product : Callable (optional, default=None)
Inner product to be used. When the argument is None, the standard
inner product will be used.
Returns
-------
numpy.ndarray
Column complement of the given array
Notes
-----
If the rows of the given matrix A span the whole space, a vector of
shape (0, n) is returned
See also
--------
col_complement, col_complement_oriented, row_complement
"""
return col_complement_oriented(A.T, inner_product).T
[docs]def orthonormalize(U, inner_product, atol=None):
"""Orthonormalize columns of matrix.
Parameters
----------
U : numpy.ndarray of shape (n, m)
Matrix of vectors to orthonormalize
inner_product : Callable
Inner product to orthonormalize for
atol : float (default=None)
If None is given, the global defaults are used. See ddg.abc.NonExact
for details.
Returns
-------
numpy.ndarray of shape (n, m)
Matrix of orthonormal vectors
"""
atol_ = get_tol_defaults()["atol"] if atol is None else atol
UN = np.array(U)
if len(np.shape(UN)) == 1:
return normalize(UN, inner_product)
ncol = np.shape(UN)[1]
for i in range(ncol):
u = UN[:, i]
uu = inner_product(u, u)
# taking absolute values since we assume no definiteness
if np.fabs(uu) < atol_:
for j in range(i + 1, len(U[0])):
v = UN[:, j]
vv = inner_product(v, v)
if np.fabs(vv) >= atol_:
copyu = np.copy(u)
copyv = np.copy(v)
UN[:, i], UN[:, j] = copyv, copyu
break
if np.fabs(inner_product(u, v)) > atol_:
UN[:, i] = u + v
UN[:, j] = u - v
UN[:, i] = normalize(u, inner_product)
un = UN[:, i]
unun = inner_product(un, un)
if np.fabs(unun) < atol_:
continue
for j in range(i + 1, ncol):
UN[:, j] = UN[:, j] - inner_product(UN[:, j], un) / unun * un
return UN
[docs]def normalize(v, inner_product, atol=None):
"""Normalize vector with respect to inner_product.
Parameters
----------
v : array_like of shape (n,m) or (n,)
Matrix of m n-vectors to normalize
inner_product: Callable
Inner product to normalize for
atol : float (default=None)
Tolerance to determine vectors of zero length
If None is given, the global defaults are used. See ddg.abc.NonExact
for details.
Returns
-------
numpy.ndarray of shape (n,m) or (n,)
Normalized vector(s)
Notes
-----
This function requires the inner product to be vectorized to handle
2-dimensional inputs.
"""
atol_ = get_tol_defaults()["atol"] if atol is None else atol
vv = np.fabs(inner_product(v, v)) # we assume no definiteness
if len(v.shape) == 1:
if vv > atol_:
return v / np.sqrt(vv)
else:
return v
vv = np.fabs(inner_product(v, v))
vv = vv.diagonal().copy() # numpy.diagonal returns non-writeable array
mask = np.where(vv <= atol_)
vv[mask] = 1.0
return v / np.sqrt(vv)
[docs]def vectorize_inner_product(inner_product, n):
"""Vectorize a given inner product for dimension n.
Parameters
----------
inner_product : Callable
Inner product to vectorize.
The function should have the signature inner_product(v, w),
and support numpy.ndarrays of shape (n,)
n : int
Dimension
Returns
-------
Callable
vectorized inner product. See module description for its signature.
"""
basis = np.eye(n, dtype=float)
M = gram_matrix(inner_product, basis)
def vectorized_inner_product(v, w):
return np.array(w).T @ M @ v
return vectorized_inner_product
[docs]def gram_matrix(inner_product, M):
"""Calculate Gram-Matrix for column vectors of M with respect to inner_product.
Parameters
----------
M : array_like of shape (n,m)
Matrix of vectors.
inner_product : Callable
Inner product to use.
The function should have the signature inner_product(v, w),
and support numpy.ndarrays of shape (n,)
Returns
-------
numpy.ndarray of shape (m,m)
Gram-Matrix for the vectors.
See Also
--------
get_matrix, polarity_matrix
"""
m = M.shape[1]
g = np.zeros((m, m))
for i in range(m):
for j in range(m):
g[i][j] = inner_product(M.T[i], M.T[j])
return g
[docs]def signature(inner_product, M):
"""Signature of matrix with respect to inner_product.
Parameters
----------
inner_product : Callable
Inner product to use
M : array_like of shape (n,m)
Matrix to calculate the signature for
Returns
-------
numpy.ndarray
Signature of matrix with respect to inner product
Notes
-----
This function requires the given inner product to be vectorized and will
give unexpected results if it is not.
"""
g = inner_product(M, M)
s = la.eigvalsh(g)
signature_ = np.empty(len(s))
for i, l in enumerate(s):
if l == 0:
signature_[i] = 0
else:
signature_[i] = l / np.fabs(l)
return np.sort(signature_)
[docs]def reflect(normal, source, inner_product, atol=None):
"""Reflect source on hyperplane.
The reflection x_ref of a vector x in a hyperplane
given by its normal n with respect to an inner product g is given by
x_ref = x - 2 * g(x, n) * n.
Parameters
----------
normal : array_like of shape (n,)
Vector defining the hyperplane
source : array_like of shape (n, m) or (n,)
Vectors to reflect
inner_product : Callable
Inner product to use for reflection
atol : float (default=None)
Tolerance used to determine whether normal is of non-zero length
If None is given, the global defaults are used. See ddg.abc.NonExact
for details.
Returns
-------
numpy.ndarray of shape (n, m) or (n,)
Reflected vector(s)
Notes
-----
This function requires the inner product to be vectorized to handle
2-dimensional inputs for the argument source.
See Also
--------
project_onto_complement
"""
atol_ = get_tol_defaults()["atol"] if atol is None else atol
n2 = inner_product(normal, normal)
if np.fabs(n2) > atol_: # we assume no definiteness
normalcomponents = np.squeeze(
np.outer(inner_product(source, normal) / n2, normal).T
)
return source - 2 * normalcomponents
else:
raise ValueError(
"Normal is too close to zero. " "Use a nonzero vector or adjust tolerances"
)
[docs]def project_onto_complement(normal, source, inner_product, atol=None):
"""Project source onto hyperplane given by its normal.
The projection x_proj of a vector x onto a hyperplane
given by its normal n with respect to an inner product g is given by
x_proj = x - g(x, n) * n.
Parameters
----------
normal : array_like of shape (n,)
Vector defining the hyperplane
source : array_like of shape (n, m) or (n,)
Vectors to project
inner_product : Callable
inner_product to use
atol : float (default=None)
Tolerance used to determine whether normal is of non-zero length
If None is given, the global defaults are used. See ddg.abc.NonExact
for details.
Returns
-------
numpy.ndarray of shape (n, m) or (n,)
Projected vector(s)
Notes
-----
This function requires the inner product to be vectorized to handle
2-dimensional inputs for the argument source.
See Also
--------
reflect
"""
atol_ = get_tol_defaults()["atol"] if atol is None else atol
n2 = inner_product(normal, normal)
if np.fabs(n2) > atol_: # we assume no definiteness
normalcomponents = np.squeeze(
np.outer(inner_product(source, normal) / n2, normal).T
)
return source - normalcomponents
else:
raise ValueError(
"Normal is too close to zero. " "Use a nonzero vector or adjust tolerances"
)
[docs]def get_matrix(inner_product, n):
"""Get representative matrix of inner_product for dimension n.
Parameters
----------
inner_product : Callable
Inner product
n : int
Dimension
Returns
-------
numpy.ndarray of shape (n,n)
Matrix representation of inner product
See Also
--------
gram_matrix, polarity_matrix
"""
eye = np.eye(n)
return gram_matrix(inner_product, eye)
[docs]def polarity_matrix(inner_product, n):
"""Polarity matrix of inner_product of dimension n.
This corresponds to the transpose of the representative matrix of the
inner product.
Parameters
----------
inner_product : Callable
Inner Product to use
n : int
Dimension
Returns
-------
numpy.ndarray of shape (n,n)
Polarity matrix of the inner product
See Also
--------
get_matrix, gram_matrix
"""
return get_matrix(inner_product, n).T
[docs]def light_like_vectors_from_onb(onb, inner_product, atol=None):
"""
Compute possible points on the quadric from an orthonormal basis by suitable
linear combinations of space- and time-like vectors. Light-like vectors will
be appended to the list.
Parameters
----------
onb : numpy.ndarray of shape (n, m)
Rows of the matrix form an orthonormal basis
with respect to the given inner product
inner_product : Callable
Inner product to use
atol : float (default=None)
Tolerance used to determine whether vectors are space/time/light-like.
If None is given, the global defaults are used. See ddg.abc.NonExact
for details.
Returns
-------
numpy.array (k,m)
All possible sums and differences of the space- and time-like vectors of the
orthonormal basis, with the light-like vectors appended, as rows of a matrix.
This means that k = k_space * k_time * 2 + k_light, where k_space, k_time and
k_light are the numbers of space-like, time-like and light-like vectors of the
onb respectively.
"""
atol_ = get_tol_defaults()["atol"] if atol is None else atol
vals = np.array([inner_product(c, c) for c in onb])
spacelike = atol_ < vals
timelike = vals < -atol_
lightlike = np.logical_not(spacelike | timelike)
# We change the shapes (k_space, n) and (k_time, n) to (1, k_space, n)
# and (k_time, 1, n). These arrays are broadcastable, resulting in an
# array of shape (k_time, k_space, n). Reshaping yields arrays of
# shape (k_space * k_time, n).
sums = np.expand_dims(onb[spacelike], 0) + np.expand_dims(onb[timelike], 1)
differences = np.expand_dims(onb[spacelike], 0) - np.expand_dims(onb[timelike], 1)
n = np.shape(onb)[1]
sums_ = sums.reshape(-1, n)
differences_ = differences.reshape(-1, n)
return np.vstack((sums_, differences_, onb[lightlike]))