"""Symmetric matrix and signature utilities.
Includes diagonalization, projective/affine normalization and signature
classes.
"""
import numpy as np
import ddg.math.projective as pmath
from ddg import nonexact
from ddg.math import linalg
# The following private functions are helper functions for affine_normalization.
def _decompose_matrix(Q):
"""
Decompose Q::
Q1 | q
Q = -----+---
q.T | r
Parameters
----------
Q : numpy.ndarray of shape (n + 1, n + 1)
Returns
-------
Q1 : numpy.ndarray of shape (n, n)
q : numpy.ndarray of shape (n,)
r : float
"""
Q1 = Q[:-1, :-1]
q = Q[:-1, -1]
r = Q[-1, -1]
return Q1, q, r
def _compute_translation(Q, rk):
"""Compute translation part of affine transformation.
Also computes the rank of Q1.
Satisfies Q1 @ b = -q in the non-parabolic case. Otherwise, computes least squares
solution of this equation.
Parameters
----------
Q : numpy.ndarray of shape (n + 1, n + 1)
Returns
-------
b : numpy.ndarray of shape (n,)
rk : int
"""
Q1, q, _ = _decompose_matrix(Q)
n = Q1.shape[0]
# Determine translation
if rk == n:
b = np.linalg.solve(Q1, -q)
else:
# This either means that Q is degenerate or parabolic
# We can still try to do our best to eliminate q by a translation
# Perform a low-rank approximation, since np.lstsq does not deal well
# with matrices that only have small entries / singular values.
u, s, vh = np.linalg.svd(Q1, hermitian=True)
s[rk:] = 0 # The singular values in s come sorted in descending order.
# Calculate pseudoinverse
s = np.array([0 if x == 0 else 1 / x for x in s])
# like b = np.linalg.lstsq(Q1, -q)[0]
b = vh.T @ np.diag(s) @ u.T @ -q
return b
def _check_parabolic(Q, b, atol, rtol):
"""Check whether we are in the parabolic case.
The case distinction is whether Q1 @ b + q = 0 has a solution or not.
Parameters
----------
Q : numpy.ndarray of shape (n + 1, n + 1)
b : numpy.ndarray of shape (n,)
Returns
-------
bool
"""
Q1, q, _ = _decompose_matrix(Q)
if not nonexact.isclose(np.linalg.norm(Q1 @ b + q), 0.0, atol=atol, rtol=rtol):
return True
ns = linalg.nullspace(Q1, atol=atol, rtol=rtol)
assert nonexact.isclose(np.linalg.norm(q.T @ ns), 0.0, atol=atol, rtol=rtol)
return False
def _parabolic_reduce_q(A1, Q, b, rk, atol, rtol):
"""Compute affine transformation that reduces q to e_n.
Parameters
----------
A1 : numpy.ndarray of shape (n, n)
Previously computed transformation that diagonalizes Q1.
Q : numpy.ndarray of shape (n + 1, n + 1)
Matrix to be normalized
b : numpy.ndarray of shape (n,)
Previously computed translation (least squares solution to Q1 @ b + q = 0).
rk : int
Previously computed rank of Q1
atol, rtol : float
Returns
-------
numpy.ndarray of shape (n + 1, n + 1)
"""
Q1, q, _ = _decompose_matrix(Q)
n = Q1.shape[0]
q2 = A1.T @ (Q1 @ b + q)
assert nonexact.isclose(np.linalg.norm(q2[:rk]), 0.0, atol=atol)
# compute affine transformation to make q have only one non-zero component
B1 = np.identity(n)
R = linalg.extend_to_basis(q2[rk:]).T
B1[rk:, rk:] = np.linalg.inv(R)
assert nonexact.isclose(
np.linalg.norm(B1.T @ q2 - linalg.e(-1, n)), 0.0, atol=atol, rtol=rtol
)
B = pmath.affine_transformation(B1, np.zeros(n))
return B
def _parabolic_eliminate_r(Q, b):
"""Compute affine transformation that eliminates last diagonal entry.
Parameters
----------
Q : numpy.ndarray of shape (n + 1, n + 1)
Matrix to be normalized
b : numpy.ndarray of shape (n,)
Returns
-------
numpy.ndarray of shape (n + 1, n + 1)
"""
# compute affine transformation (translation) to make r2 vanish
Q1, q, r = _decompose_matrix(Q)
n = Q1.shape[0]
r2 = np.dot(b, Q1 @ b + 2 * q) + r
b2 = -r2 / 2.0 * linalg.e(-1, n)
B = pmath.affine_transformation(np.identity(n), b2)
return B
def _normal_form_permutation(sgn, parabolic):
"""Return permutation matrix B such that A @ B produces a diagonal matrix of the
normal form specified in the AffineSignature docstring.
Parameters
----------
sgn : Sequence of 1, -1 and 0
This has to be sorted in the order 1, -1, 0 except for the last entry (or last
two entries, in the parabolic case).
parabolic : bool
Returns
-------
numpy.ndarray of shape (n + 1, n + 1)
Permutation matrix.
"""
B = np.eye(len(sgn))
mask = np.full(len(sgn), True)
if parabolic:
mask[[-2, -1]] = False
else:
mask[-1] = False
p = np.count_nonzero(sgn[mask] == 1)
q = np.count_nonzero(sgn[mask] == -1)
z = np.count_nonzero(sgn[mask] == 0)
# Determine whether to swap the first two blocks of diag[mask], i.e. go
# from (1,...,1, -1,...,-1, 0,...,0) to (-1,...,-1, 1,...,1, 0,...,0).
# The case sgn[-1] == 1 is the default since sgn is pre-sorted
# in the order 1, -1, 0.
if (sgn[-1] == 0 and p < q) or sgn[-1] == -1:
permutation = np.concatenate(
(np.arange(q) + p, np.arange(p), np.arange(z) + p + q)
)
B[:, mask] = B[:, mask][:, permutation]
if parabolic:
B[-2, -2] = -1 # Make the off-diagonal entries -1
return B
[docs]def signature_sort_key(x):
"""Sorting key function for signatures.
Desired order of values: positive, negative, zero with each category sorted
by absolute value (decreasing)
Parameters
----------
x : float
Returns
-------
tuple(int, float)
Examples
--------
>>> from ddg.math.symmetric_matrices import signature_sort_key
>>> signature_sort_key(1)
(0, -1)
>>> signature_sort_key(-1)
(1, -1)
>>> signature_sort_key(0)
(2, 0)
>>> signature_sort_key(15)
(0, -15)
>>> signature_sort_key(15.1)
(0, -15.1)
"""
sgn = 2
if np.sign(x) == +1:
sgn = 0
elif np.sign(x) == -1:
sgn = 1
return (sgn, -np.abs(x))
[docs]def diagonalize(A, sort_key=signature_sort_key, atol=0, rtol=0):
"""Diagonalize a symmetric matrix.
Diagonalization of symmetric matrix but with customizable order of
eigenvalues. Uses np.linalg.eigh, then eigenvalues and eigenvectors.
Makes eigenvalues close to zero exact.
Parameters
----------
A : numpy.ndarray of shape (n, n)
Symmetric matrix
sort_key : Callable (default=signature_sort_key)
Sorting key function. Eigenvalues are sorted according to the result of
this function for each eigenvalue.
atol, rtol : float (default=0)
Eigenvalues close to 0 according to these tolerances are set to 0.
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See :py:mod:`ddg.nonexact` for details.
Returns
-------
ev : list
List of `n` Eigenvalues, sorted according to `sort_key`.
S : numpy.ndarray of shape (n, n)
Columns of `S` are the `n` Eigenvectors in the order of `ev`.
Examples
--------
>>> import numpy as np
>>> from ddg.math.symmetric_matrices import diagonalize
>>> A = np.array([[10, 0, 0], [0, 5, 0], [0, 0, 3]])
>>> A
array([[10, 0, 0],
[ 0, 5, 0],
[ 0, 0, 3]])
>>> ev, S = diagonalize(A)
>>> ev
(10.0, 5.0, 3.0)
>>> S
array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]])
>>> B = np.array([[1, 0, 2], [0, 5, 0], [2, 0, 3]])
>>> B
array([[1, 0, 2],
[0, 5, 0],
[2, 0, 3]])
>>> ev, S = diagonalize(B)
>>> ev
(5.0, 4.23606797749979, -0.2360679774997897)
>>> S
array([[ 0. , -0.52573111, -0.85065081],
[-1. , 0. , 0. ],
[ 0. , -0.85065081, 0.52573111]])
"""
ev, S = np.linalg.eigh(A)
ev = [0.0 if nonexact.isclose(x, 0.0, atol=atol, rtol=rtol) else x for x in ev]
# enumerate eigenvalues, sort and retrieve permutation
ev_enumerated = [(ev[i], i) for i in range(len(ev))]
ev_enumerated.sort(key=lambda x: sort_key(x[0]))
ev, permutation = zip(*ev_enumerated)
# apply permutation to eigenvectors
S = np.array([S.T[i] for i in permutation]).T
return ev, S
################################
# Symmetry and diagonal matrices
################################
[docs]def symmetrize(A):
"""Symmetrize square matrix.
Parameters
----------
A : numpy.ndarary of shape (n, n)
Returns
-------
numpy.ndarray of shape (n, n)
``(A.T + A) / 2``
Examples
--------
>>> import numpy as np
>>> from ddg.math.symmetric_matrices import symmetrize
>>> A = np.array([[1, 0, 2], [0, 5, 0], [2, 0, 3]])
>>> A
array([[1, 0, 2],
[0, 5, 0],
[2, 0, 3]])
>>> symmetrize(A)
array([[1., 0., 2.],
[0., 5., 0.],
[2., 0., 3.]])
>>> B = np.array([[1, 1, 2], [0, 5, 0], [2, 0, 3]])
>>> B
array([[1, 1, 2],
[0, 5, 0],
[2, 0, 3]])
>>> symmetrize(B)
array([[1. , 0.5, 2. ],
[0.5, 5. , 0. ],
[2. , 0. , 3. ]])
"""
return (A.T + A) / 2.0
[docs]def is_symmetric(A, atol=1e-13, rtol=None):
"""Checks whether a square matrix is symmetric.
Parameters
----------
A : numpy.ndarray of shape (n, n)
atol : float (default=1e-13)
rtol : float (default=None)
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See :py:mod:`ddg.nonexact` for details.
Returns
-------
bool
Examples
--------
>>> import numpy as np
>>> from ddg.math.symmetric_matrices import is_symmetric
>>> A = np.array([[1, 0, 2], [0, 5, 0], [2, 0, 3]])
>>> A
array([[1, 0, 2],
[0, 5, 0],
[2, 0, 3]])
>>> is_symmetric(A)
True
>>> B = np.array([[1, 1, 2], [0, 5, 0], [2, 0, 3]])
>>> B
array([[1, 1, 2],
[0, 5, 0],
[2, 0, 3]])
>>> is_symmetric(B)
False
"""
return nonexact.isclose(np.linalg.norm(A - A.T), 0.0, atol=atol, rtol=rtol)
[docs]def symmetric_matrix_from_diagonal(sgn, parabolic=False):
"""
A normalized affine quadric can be reconstructed from its diagonal
and the additional knowledge whether it is parabolic.
For what "parabolic" means, see the docstring of AffineSignature.
Parameters
----------
sgn : array_like of shape (n,)
Diagonal of the matrix.
parabolic : bool (default=False)
If this is True, `sgn` must have 0 at positions -1 and -2.
Raises
------
ValueError
If `parabolic` is True and `sgn` does not have 0 at positions -1 and
-2.
See also
--------
AffineSignature
Examples
--------
>>> from ddg.math.symmetric_matrices import symmetric_matrix_from_diagonal
>>> symmetric_matrix_from_diagonal((1, 0, 1))
array([[1, 0, 0],
[0, 0, 0],
[0, 0, 1]])
>>> symmetric_matrix_from_diagonal([10, 2, 15, 1])
array([[10, 0, 0, 0],
[ 0, 2, 0, 0],
[ 0, 0, 15, 0],
[ 0, 0, 0, 1]])
>>> symmetric_matrix_from_diagonal((1, 2, 0, 0), parabolic=True)
array([[1, 0, 0, 0],
[0, 2, 0, 0],
[0, 0, 0, 1],
[0, 0, 1, 0]])
"""
Q = np.diag(sgn)
if parabolic:
if sgn[-1] != 0 or sgn[-2] != 0:
raise ValueError(
"This is not the diagonal of a normalized parabolic quadric: "
+ str(sgn)
)
Q[-1, -2] = Q[-2, -1] = 1.0
return Q