"""
Utility functions for symmetric matrices.
Includes diagonalization, and projective/affine normalization.
"""
import numpy as np
import ddg.math.projective as putils
import ddg.math.linalg as linalg
from dataclasses import dataclass
from typing import Optional, Union
from ddg.abc import NonExact
[docs]@dataclass(frozen=True, eq=False)
class Signature:
"""Projective signature class.
Two projective signatures are equal if and only if quadrics with those
signatures are related by a projective transformation.
We want to identify:
* Signatures
* Quadrics up to projective transformation
* Normalized matrices up to nonzero scalar multiplication (-1 in
particular), which will be the result of e.g.
:py:func:`projective_normalization`.
To achieve this, we define a normalized matrix as a diagonal matrix ::
I_n | |
-----|------|---
| -I_m |
-----|------|---
| | 0
or ::
-I_n | |
------|-----|---
| I_m |
------|-----|---
| | 0
where I_n is the (n x n)-identity matrix, 0 means a 0 matrix and
n >= m.
Parameters
----------
plus : int
minus : int
zero : int (default=0)
Number of 1, -1, 0 on the diagonal.
"""
plus: int
minus: int
zero: int = 0
@property
def matrix(self):
"""Returns the normalized matrix corresponding to the signature.
See the class docstring for how this is defined.
Returns
-------
numpy.ndarray of shape (n, n)
Notes
-----
This method respects the sign the signature is initialized with, even
though the equality operator does not.
"""
p, q, z = self.plus, self.minus, self.zero
diag = max([p, q]) * [1] + min([p, q]) * [-1] + z * [0]
# respect the actual counts
if p < q:
diag = -np.array(diag)
return symmetric_matrix_from_diagonal(diag)
@property
def is_degenerate(self):
"""Whether this is the signature of a degenerate quadric.
Returns
-------
bool
"""
return self.zero > 0
@property
def rank(self):
"""Rank of associated matrix.
Returns
-------
int
"""
return self.plus + self.minus
@property
def is_positive_definite(self):
"""Whether the signature is positive definite.
Meaning, if q is a quadratic form with this signature, q(p,p) > 0 for
all p != 0.
Returns
-------
bool
"""
return self.plus > 0 and self.minus == 0 and self.zero == 0
@property
def is_negative_definite(self):
"""Whether the signature is negative definite.
Meaning, if q is a quadratic form with this signature, q(p,p) < 0 for
all p != 0.
Returns
-------
bool
"""
return self.minus > 0 and self.plus == 0 and self.zero == 0
@property
def is_positive_semi_definite(self):
"""Whether the signature is positive semi-definite.
Meaning, if q is a quadratic form with this signature, q(p,p) >= 0 for
all p.
Returns
-------
bool
"""
return self.minus == 0
@property
def is_negative_semi_definite(self):
"""Whether the signature is negative semi-definite.
Meaning, if q is a quadratic form with this signature, q(p,p) <= 0 for
all p.
Returns
-------
bool
"""
return self.plus == 0
@property
def is_indefinite(self):
"""Whether the signature is positive definite.
Meaning, if q is a quadratic form with this signature, q can have both
positive and negative values.
Returns
-------
bool
"""
return self.plus > 0 and self.minus > 0
@property
def is_definite(self):
"""Whether the signature is either positive or negative definite.
Returns
-------
bool
"""
return self.is_positive_definite or self.is_negative_definite
@property
def is_semi_definite(self):
"""Whether the signature is either positive or negative semi-definite.
Returns
-------
bool
"""
return self.is_positive_semi_definite or self.is_negative_semi_definite
@property
def _counts_up_to_sign(self):
"""Similar to _counts, but ignores order of first two entries.
Returns
-------
tuple(frozenset(int, int), int)
"""
return (frozenset([self.plus, self.minus]), self.zero)
def __eq__(self, other):
if not isinstance(other, Signature):
return NotImplemented
return self._counts_up_to_sign == other._counts_up_to_sign
def __str__(self):
result = (self.plus, self.minus)
if self.is_degenerate:
result += (self.zero,)
return str(result)
def __hash__(self):
return hash(self._counts_up_to_sign)
[docs]@dataclass(frozen=True, eq=False)
class AffineSignature(Signature):
"""Affine signature class.
Two affine signatures are equal if and only if quadrics with those
signatures are related by an affine transformation.
We want to identify:
* Affine signatures
* Quadrics up to affine transformation
* Normalized matrices up to nonzero scalar multiplication (-1 in
particular), which will be the result of e.g.
:py:func:`affine_normalization`.
To achieve this, we define what a normalized matrix is as follows:
In the non-parabolic case: A diagonal matrix with only 1, -1 and 0 on the
diagonal. The diagonal **without the entry at i** looks as follows:
* If diag[i] == 0: (s,...,s, -s,...,-s, 0,...,0), and the number of s's
is greater than or equal to the number of -s's.
* If diag[i] == s != 0: (s,...,s, -s,...,-s, 0,...,0).
In the parabolic case: a matrix of the form ::
D1 | |
----|-----|----
| 0 s |
| s 0 | < i
----|-----|----
| | D2
^
i
for i != 0 or, for i == 0, ::
0 | | s
---|---|---
| D |
---|---|---
s | | 0
where D1, D2 and D are diagonal matrices with only 1, -1 and 0 on the
diagonal. The diagonal is sorted in the order s, -s, 0 except at
positions i and i-1 and the number of s's is greater than or equal to
the number of -s's.
Note that for parabolic cases, flipping the sign of the diagonal is the
same as flipping the sign of the off-diagonal entries, which as an affine
transformation is a reflection about the plane x_{i-1} = 0.
Parameters
----------
plus : int
minus : int
zero : int (default=0)
Number of 1, -1 or 0 on the diagonal of the **diagonalized** matrix.
Note that for parabolic cases, these will not be the counts on the
diagonal: A 1 and a -1 will be replaced by a 0 and entries on the
off-diagonal, see above.
affine_component_entry : 1, -1, 0 or 'parabolic' (default=None)
The entry at `affine_component` on the diagonal. A value of 'parabolic'
means the matrix [[0, 1], [1, 0]], see above for details. This can be
omitted (left as None) if only one of `plus`, `minus` or `zero` is
nonzero.
affine_component : int (default=-1)
Raises
------
ValueError
* If `affine_component_entry` is 'parabolic' and `plus` or `minus` are
0
* If `affine_component_entry` was not given and could not be inferred
uniquely.
* If `affine_component_entry` was given and is not 1, -1, 0 or
'parabolic'.
"""
plus: int
minus: int
zero: int = 0
affine_component_entry: Optional[Union[int, str]] = None
affine_component: int = -1
def __post_init__(self):
# Just consistency checking
p, q, r = self.plus, self.minus, self.zero
acs = self.affine_component_entry
# We can only check if the entry can be inferred, but can not set
# self.affine_component_entry due to immutability.
if acs is None:
entries = set(p * [1] + q * [-1] + r * [0])
if len(entries) != 1:
raise ValueError(
"Sign of affine component entry could not be inferred "
"uniquely. Please specify 1, -1, 0 or 'parabolic'."
)
elif acs not in [1, -1, 0, 'parabolic']:
raise ValueError(
"Sign of affine component must be 1, -1, 0 or 'parabolic'."
)
elif acs == 'parabolic' and (p == 0 or q == 0):
raise ValueError(
"Parabolic signatures must have a 1 and a -1."
)
elif ((acs == 1 and p == 0) or
(acs == -1 and q == 0) or
(acs == 0 and r == 0)):
raise ValueError(
"The affine componenent sign is not contained in the diagonal."
)
@property
def matrix(self):
"""Return normalized matrix corresponding to the signature.
See the class docstring for how this is defined.
Returns
-------
numpy.ndarray of shape (n, n)
Notes
-----
This method respects the sign the signature is initialized with, even
though the equality operator does not.
"""
p, q, z = self.plus, self.minus, self.zero
# We already checked in __post_init__ that only one is nonzero.
if self.affine_component_entry is None:
for count, s in zip([p, q, z], [1, -1, 0]):
if count > 0:
break
else:
s = self.affine_component_entry
i = self.affine_component
parabolic = (s == 'parabolic')
# We need this to change the sign of the off-diagonal entries in the
# parabolic case
sign = 1
# build the diagonal without affine_component_entry
if parabolic:
diag_ex = ((max([p, q]) - 1) * [1]
+ (min([p, q]) - 1) * [-1]
+ z * [0])
# respect the actual counts
if p < q:
sign = -1
elif s == 0:
diag_ex = max([p, q]) * [1] + min([p, q]) * [-1] + (z-1) * [0]
# respect the actual counts
if p < q:
sign = -1
elif s == 1:
diag_ex = (p-1) * [1] + q * [-1] + z * [0]
elif s == -1:
diag_ex = (q-1) * [-1] + p * [1] + z * [0]
# build the actual diagonal
# dtype int for consistency with Signature.matrix
diag = np.empty(p + q + z, dtype=int)
mask = np.full(p + q + z, True)
if parabolic:
mask[[i, i-1]] = False
diag[[i, i-1]] = 0
else:
mask[i] = False
diag[i] = s
diag[mask] = diag_ex
return sign * symmetric_matrix_from_diagonal(
diag,
parabolic=parabolic,
affine_component=i
)
def __eq__(self, other):
i = self.affine_component
if not isinstance(other, AffineSignature):
return NotImplemented
# convert to non-negative index before comparing
r = self.plus + self.minus + self.zero
if i % r != other.affine_component % r:
return False
return np.array_equal(self.matrix, other.matrix)\
or np.array_equal(self.matrix, -other.matrix)
def __str__(self):
lines = [super().__str__()]
if self.affine_component_entry is not None:
lines.append(
f"affine component entry: {self.affine_component_entry}"
)
lines.append(f"affine component: {self.affine_component}")
return "\n".join(lines)
def __hash__(self):
"""Hash function.
Computed from self.matrix and self.affine_component (mod length).
"""
r = self.plus + self.minus + self.zero
M = self.matrix.flatten()
for m in M:
if m != 0:
M *= m
break
return hash(tuple(M) + (self.affine_component % r,))
[docs]def signature_from_diagonal(diag):
"""Create signature from diagonal.
Parameters
----------
diag : numpy.ndarray of shape (n,) or list
1D array containing only 1, -1 and 0
Returns
-------
Signature
Raises
------
ValueError
If `diag` contains entries other than 1, -1 or 0.
"""
diag = np.array(diag)
p = np.count_nonzero(diag == 1)
q = np.count_nonzero(diag == -1)
z = np.count_nonzero(diag == 0)
if p + q + z != len(diag):
raise ValueError("Diagonal can only contain 1, -1 and 0")
return Signature(p, q, z)
[docs]def affine_signature_from_diagonal(diag, parabolic=False, affine_component=-1):
"""Create affine signature from diagonal.
Parameters
----------
diag : numpy.ndarray of shape (n,) or list
1D array containing only 1, -1 and 0. If `parabolic` is True, this must
have 0 at entries i and i-1, where i is `affine_component`.
parabolic : bool (default=False)
affine_component : int (default=-1)
Returns
-------
AffineSignature
Raises
------
ValueError
* If `diag` contains entries other than 1, -1 or 0.
* If `parabolic` is True and `diag` does not have 0 at positions
i and i-1, where i is `affine_component`.
"""
diag = np.array(diag)
p = np.count_nonzero(diag == 1)
q = np.count_nonzero(diag == -1)
z = np.count_nonzero(diag == 0)
if p + q + z != len(diag):
raise ValueError("Diagonal can only contain 1, -1 and 0")
if parabolic and (
diag[affine_component] != 0 or diag[affine_component-1] != 0):
raise ValueError("parabolic diagonals must have 0 at positions "
"affine_component and affine_component-1.")
if parabolic:
affine_component_entry = 'parabolic'
p += 1
q += 1
z -= 2
else:
affine_component_entry = diag[affine_component]
return AffineSignature(p, q, z,
affine_component_entry=affine_component_entry,
affine_component=affine_component)
[docs]def symmetrize(A):
"""Symmetrize square matrix"""
return (A.T + A) / 2.
[docs]@NonExact.nonexact_function
def is_symmetric(A, atol=1e-13, rtol=None):
"""Checks whether A is symmetric
Notes
-----
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See ddg.abc.NonExact for details.
"""
return np.isclose( np.linalg.norm(A - A.T), 0., atol=atol, rtol=rtol )
[docs]def signature_sort_key(x):
"""
Order of signature values: positive, negative, zero
For eigenvalues each class is sorted by absolute value (decreasing)
"""
sgn = 2
if np.sign(x) == +1:
sgn = 0
elif np.sign(x) == -1:
sgn = 1
return (sgn, -np.abs(x))
[docs]@NonExact.nonexact_function
def diagonalize(A, sort_key=signature_sort_key, atol=0, rtol=0):
"""
Diagonalization of symmetric matrix but with customizable order of
eigenvalues. Uses np.linalg.eigh, then eigenvalues and eigenvectors.
Makes eigenvalues close to zero exact.
Notes
-----
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See ddg.abc.NonExact for details.
"""
ev, S = np.linalg.eigh(A)
ev = [0. if np.isclose(x,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
[docs]@NonExact.nonexact_function
def projective_normalization(Q, atol=None, rtol=None):
"""Bring quadric to normal form by projective transformation.
Returns signature sgn and transformation matrix A such that ::
A.T @ Q @ A == sgn.matrix
See :py:class:`Signature` for more information.
Note that this means that the point transformation defined by A
transforms ``sgn.matrix`` into Q.
Parameters
----------
Q : numpy.ndarray
Symmetric matrix.
atol : float
rtol : float
Returns
-------
sgn : Signature
A : numpy.ndarray
Transformation matrix.
Notes
-----
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See ddg.abc.NonExact for details.
"""
if not is_symmetric(Q, atol=atol, rtol=rtol):
raise ValueError('The matrix is not symmetric: ' + str(Q) )
ev, A = diagonalize(Q, atol=atol, rtol=rtol)
# Add affine scaling to the transformation
B = np.diag([1. if x==0 else 1/np.sqrt(np.abs(x)) for x in ev])
A = np.dot(A, B)
# compute signature
sgn = np.array([int(np.sign(x)) for x in ev])
# permutation step to bring the matrix to normal form as defined in the
# Signature docstring.
# We might have to swap the first two blocks of diag, i.e. go from
# (1,...,1, -1,...,-1, 0,...,0) to (-1,...,-1, 1,...,1, 0,...,0).
B = np.eye(len(sgn))
p = np.count_nonzero(sgn == 1)
q = np.count_nonzero(sgn == -1)
z = np.count_nonzero(sgn == 0)
if q > p:
permutation = np.concatenate((np.arange(q) + p,
np.arange(p),
np.arange(z) + p + q))
B = B[:,permutation]
A = A @ B
sgn = signature_from_diagonal(sgn)
return sgn, A
[docs]@NonExact.nonexact_function
def affine_normalization(Q, affine_component=-1, atol=None, rtol=None):
"""Bring quadric to normal form by affine transformation.
Up to parabolic cases this means diagonalizing Q by an affine
transformation.
Returns affine transformation A with affine component i such that ::
A.T @ Q @ A == sgn.matrix
See :py:class:`AffineSignature` for more information.
Parameters
----------
Q : numpy.ndarray of shape (n,n)
Symmetric matrix.
i : int
Affine component.
atol : float
rtol : float
Returns
-------
sgn : AffineSignature
A : numpy.ndarray of shape (n,n)
Affine transformation.
Notes
-----
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See ddg.abc.NonExact for details.
"""
if not is_symmetric(Q, atol=atol, rtol=rtol):
raise ValueError('The matrix is not symmetric: ' + str(Q) )
# Decompose Q:
#
# Q11 | q1 | Q12
# -----|----|-----
# Q = q1.T | r | q2.T <--- i
# -----|----|-----
# Q13 | q2 | Q14
#
# ^
# i
#
# Resulting in a scalar r, and a vector and matrix
#
# q1
# q = --
# q2
#
# Q11 | Q12
# Q1 = ----|----
# Q12 | Q14
N = Q.shape[0] - 1
# make indices non-negative because of np.insert
i = affine_component
if -N - 1 <= i <= -1:
i += N + 1
indices = np.delete(np.arange(N + 1), i)
Q1 = Q[np.ix_(indices, indices)]
q = Q[indices, i]
r = Q[i,i]
rk = linalg.rank(Q1, atol=atol, rtol=rtol)
parabolic = False
# 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)
# TODO: s[rk:] = 0 should suffice, but due to a bug that was fixed in
# numpy 1.19.0, s is not sorted correctly when using hermitian=True.
# Set the requirement to numpy >= 1.19.0?
s[np.argsort(s)[:N-rk]] = 0
# 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
if not np.isclose(np.linalg.norm(Q1 @ b + q), 0.,
atol=atol, rtol=rtol):
# parabolic case (impossible to normalize q to 0)
parabolic = True
else:
ns = linalg.nullspace(Q1, atol=atol, rtol=rtol)
assert np.isclose(np.linalg.norm(q.T @ ns), 0.,
atol=atol, rtol=rtol)
# Determine eigenvalues of Q1 and orthogonal part of transformation
ev, A1 = diagonalize(Q1, atol=atol, rtol=rtol)
# Write Euclidean trafo in homogeneous coordinates
A = putils.affine_transformation(A1, b, affine_component=i)
# Adjust A for parabolic case / compute new entries of Q (if necessary)
if parabolic:
# compute new r
r = np.dot(b, Q1 @ b + 2 * q) + r
# compute new q
q = A1.T @ (Q1 @ b + q)
assert np.isclose(np.linalg.norm(q[:rk]), 0., atol=atol)
# compute affine transformation to make q have only one non-zero
# component
B1 = np.identity(N)
B1[rk:,rk:] = np.linalg.inv(linalg.extend_to_basis(q[rk:], i=-1).T)
# permute columns so the correct component (i-1) is non-zero
B1 = np.hstack((B1[:,:i-1], B1[:,[-1]], B1[:,i-1:-1]))
# we also have to permute the eigenvalues
ev = np.hstack((ev[:i-1], ev[-1], ev[i-1:-1]))
assert np.isclose(np.linalg.norm(B1.T @ q - linalg.e(i-1, N)), 0.,
atol=atol)
B = putils.affine_transformation(B1, np.zeros(N), affine_component=i)
A = A @ B
# compute affine transformation (translation) to make r vanish
b2 = np.zeros(N)
b2[i-1] = -r/2.
B = putils.affine_transformation(np.identity(N), b2,
affine_component=i)
A = A @ B
r = 0.
else:
# compute new r
r = np.dot(b, q) + r
r = 0. if np.isclose(r, 0, atol=atol, rtol=rtol) else r
# Add last eigenvalue from Q
ev = np.insert(ev, i, r)
# Add affine scaling to the transformation
B = np.diag([1. if x==0 else 1/np.sqrt(np.abs(x)) for x in ev])
A = A @ B
# compute signature
# This is sorted in the order 1, -1, 0, except for entry i or, in the
# parabolic case, entries i and i-1.
sgn = np.array([int(np.sign(x)) for x in ev])
# Final permutation step to the bring the matrix to normal form as defined
# in the AffineSignature docstring.
B = np.eye(len(sgn))
mask = np.full(len(sgn), True)
if parabolic:
mask[[i,i-1]] = False
else:
mask[i] = 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[i] == 1 is the default since sgn is pre-sorted
# in the order 1, -1, 0.
if (sgn[i] == 0 and p < q) or sgn[i] == -1:
permutation = np.concatenate((np.arange(q) + p,
np.arange(p),
np.arange(z) + p + q))
B[:,mask] = B[:,mask][:,permutation]
if parabolic:
B[i-1,i-1] = -1 # Make the off-diagonal entries -1
A = A @ B
sgn = affine_signature_from_diagonal(
sgn,
parabolic=parabolic,
affine_component=affine_component
)
return sgn, A
[docs]def symmetric_matrix_from_diagonal(sgn, parabolic=False, affine_component=-1):
"""
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 : numpy.ndarray of shape (n,)
Diagonal of the matrix.
parabolic : bool (default=False)
If this is True, `sgn` must have 0 at positions i and i-1, where i is
`affine_component`.
affine_component : int (default=-1)
Only used when `parabolic` is True.
Raises
------
ddg.geometry.projective.exceptions.ValueError
When `parabolic` is True but there are no 0s at the relevant positions,
see above.
See also
--------
AffineSignature
"""
i = affine_component
j = affine_component-1
Q = np.diag(sgn)
if parabolic:
if sgn[i] != 0 or sgn[j] != 0:
raise ValueError(
'This is not the diagonal of a normalized parabolic quadric: '
+ str(sgn)
)
Q[i,j] = Q[j,i] = 1.
return Q