"""Symmetric matrix and signature utilities.
Includes diagonalization, projective/affine normalization and signature
classes.
"""
from dataclasses import dataclass
from typing import Literal, Optional
import numpy as np
import ddg.math.projective as pmath
from ddg import nonexact
from ddg.math import linalg
###################
# Signature classes
###################
[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)
These three arguments mean the number of 1, -1, 0 on the diagonal.
Notes
-----
This class implements the equals relation ``==``, with the meaning explained above.
It is also hashable.
"""
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.
This method respects the global sign the signature is initialized with,
even though the equality operator does not.
Returns
-------
numpy.ndarray of shape (n, n)
"""
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 indefinite.
Meaning, that a quadratic form with this signature, 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):
"""Returns counts of 1, -1 and 0, ignoring overall sign.
Ignoring overall sign means that we want the signatures ``(p, q, z)`` and
``(q, p, z)`` to be treated as equal. This method encodes the counts in
such a way that this is the case.
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 looks as follows (where ``s`` is either 1, -1 or 0):
* If the last entry on the diagonal is 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 the last entry on the diagonal is ``s != 0``:
``(s,...,s, -s,...,-s, 0,...,0, s)``.
In the parabolic case: a matrix of the form ::
D |
--+-----
| 0 s
| s 0
where ``D`` is a diagonal matrix with only 1, -1 and 0 on the diagonal. The diagonal
is sorted in the order ``s``, ``-s``, ``0`` 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[-2] = 0``.
Parameters
----------
plus : int
minus : int
zero : int (default=0)
These three arguments mean the 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.
last_entry : 1, -1, 0 or "parabolic" (default=None)
The last entry 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.
Raises
------
ValueError
* If `last_entry` is "parabolic" and `plus` or `minus` are 0.
* If `last_entry` was not given and could not be inferred uniquely.
* If `last_entry` was given and is not 1, -1, 0 or 'parabolic'.
Notes
-----
This class implements the equals relation ``==``, with the meaning explained above.
It is also hashable.
"""
plus: int
minus: int
zero: int = 0
last_entry: Optional[Literal[0, 1, -1, "parabolic"]] = None
def __post_init__(self):
# Just consistency checking
p, q, r = self.plus, self.minus, self.zero
s = self.last_entry
# We can only check if the entry can be inferred, but can not set
# self.last_entry due to immutability.
if s is None:
entries = set(p * [1] + q * [-1] + r * [0])
if len(entries) != 1:
raise ValueError(
"Sign of last entry could not be inferred uniquely. Please specify "
"1, -1, 0 or 'parabolic'."
)
elif s not in [1, -1, 0, "parabolic"]:
raise ValueError("Sign of last entry must be 1, -1, 0 or 'parabolic'.")
elif s == "parabolic" and (p == 0 or q == 0):
raise ValueError("Parabolic signatures must have a 1 and a -1.")
elif (s == 1 and p == 0) or (s == -1 and q == 0) or (s == 0 and r == 0):
raise ValueError("The last 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.last_entry is None:
for count, s in zip([p, q, z], [1, -1, 0]):
if count > 0:
break
else:
s = self.last_entry
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 last_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)
if parabolic:
diag[-2:] = 0
diag[:-2] = diag_ex
else:
diag[-1] = s
diag[:-1] = diag_ex
return sign * symmetric_matrix_from_diagonal(diag, parabolic=parabolic)
def __eq__(self, other):
if not isinstance(other, AffineSignature):
return NotImplemented
return np.array_equal(self.matrix, other.matrix) or np.array_equal(
self.matrix, -other.matrix
)
def __str__(self):
lines = [super().__str__()]
if self.last_entry is not None:
lines.append(f"last entry: {self.last_entry}")
return "\n".join(lines)
def __hash__(self):
"""Hash function computed from self.matrix()."""
M = self.matrix.flatten()
for m in M:
if m != 0:
M *= m
break
return hash(tuple(M))
[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):
"""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 -1 and -2.
parabolic : bool (default=False)
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 -1 and
-2.
"""
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[-1] != 0 or diag[-2] != 0):
raise ValueError("parabolic diagonals must have 0 at positions -1 and -2.")
if parabolic:
last_entry = "parabolic"
p += 1
q += 1
z -= 2
else:
last_entry = diag[-1]
return AffineSignature(p, q, z, last_entry=last_entry)
###################################
# Normalization and diagonalization
###################################
[docs]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, 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
-------
sgn : Signature
A : numpy.ndarray
Transformation matrix.
"""
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.0 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
# 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 affine_normalization(Q, atol=None, rtol=None):
"""Bring quadric to normal form by affine transformation.
Returns affine signature `sgn` and affine transformation `A` such that ::
A.T @ Q @ A == sgn.matrix
For non-parabolic cases, this means diagonalizing `Q` with an affine transformation
such that there is only 1, -1 or 0 on the diagonal. For parabolic cases, normalizes
to a matrix of the form ::
D1 |
(+/-) ---+-----
| 0 1
| 1 0
Where ``D1`` is a diagonal matrix with only 1, -1 or 0 on the diagonal. For more
information about the order of entries on the diagonal etc, see
:py:class:`.AffineSignature`.
Parameters
----------
Q : numpy.ndarray of shape (n + 1, n + 1)
Symmetric matrix. Can be interpreted as a quadric in n-dim. projective space.
atol, 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
-------
sgn : AffineSignature
A : numpy.ndarray of shape (n + 1, n + 1)
Affine transformation.
Raises
------
ValueError
If *Q* is not symmetric.
Notes
-----
For a detailed explanation of the implementation, see :ref:`supplemental material
<affine_normalization>`.
See also
--------
ddg.math.symmetric_matrices.AffineSignature
"""
if not is_symmetric(Q, atol=atol, rtol=rtol):
raise ValueError("The matrix is not symmetric: " + str(Q))
# We treat this as a special case because the decomposition will produce empty
# arrays Q1, q which many functions can't deal with.
if Q.shape[0] == 1:
r = Q[0, 0]
if nonexact.isclose(r, 0, atol=atol, rtol=rtol):
return AffineSignature(0, 0, 1), np.array([[1]])
else:
return AffineSignature(1, 0), np.array([[1 / np.sqrt(np.abs(r))]])
# List of affine transformations that must be applied to normalize Q.
trafos = []
Q1, q, r = _decompose_matrix(Q)
rk = linalg.rank(Q1, atol=atol, rtol=rtol)
b = _compute_translation(Q, rk)
parabolic = _check_parabolic(Q, b, atol, rtol)
ev_, A1 = diagonalize(Q1, atol=atol, rtol=rtol)
trafos.append(pmath.affine_transformation(A1, b))
if parabolic:
trafos.append(_parabolic_reduce_q(A1, Q, b, rk, atol, rtol))
trafos.append(_parabolic_eliminate_r(Q, b))
ev = np.append(ev_, 0.0)
else:
r2 = np.dot(b, q) + r
r2 = 0.0 if nonexact.isclose(r2, 0, atol=atol, rtol=rtol) else r2
ev = np.append(ev_, r2)
# We now normalize the diagonal entries to 1, -1 or 0.
# diagonalize already sets eigenvalues close to 0 according to tolerances to 0 and
# we did the same for the last eigenvalue manually, so an exact check is fine here.
trafos.append(np.diag([1.0 if x == 0 else 1 / np.sqrt(np.abs(x)) for x in ev]))
# compute signature
# This is sorted in the order 1, -1, 0, except for the last entry (or last two, in
# the parabolic case).
sgn_diag = np.sign(ev).astype(int)
# Final permutation step to achieve normal form as described in AffineSignature.
trafos.append(_normal_form_permutation(sgn_diag, parabolic))
return (
affine_signature_from_diagonal(sgn_diag, parabolic=parabolic),
np.linalg.multi_dot(trafos),
)
[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)
"""
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`.
"""
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``
"""
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
"""
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 : numpy.ndarray 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
"""
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