refactor: Improve PSD typing (#522)

* Add PSDTensor
* Remove GeneralizedMatrix
* Use classes for Matrix, PSDMatrix and PSDTensor instead of type annotations
* Use casting in compute_gramian, normalize and regularize
* Add typeguard functions is_matrix, is_psd_tensor and is_psd_matrix
* Move normalize and regularize to _linalg
* Move _check_is_matrix to `Aggregator.__call__`
* Improve internal type hints
* Rename reshape_gramian to reshape and movedim_gramian to movedim
* Add _gramian_utils.flatten
* Rename a few tests
* Add some parametrizations to test_gramian_is_psd
* Add test_reshape_yields_psd, test_flatten_yields_matrix, test_flatten_yields_psd, test_movedim_yields_psd, test_normalize_yields_psd and test_regularize_yields_psd
* Add assert_is_psd_tensor

---------

Co-authored-by: Valérian Rey <valerian.rey@gmail.com>

Author: PierreQuinton

pyproject.toml

@@ -107,3 +107,9 @@ full = [
 
 [tool.pytest.ini_options]
 xfail_strict = true
+
+[tool.coverage.report]
+exclude_lines = [
+    "pragma: not covered",
+    "@overload",
+]

src/torchjd/_linalg/__init__.py

@@ -1,4 +1,14 @@
-from .gramian import compute_gramian
-from .matrix import Matrix, PSDMatrix
+from ._gramian import compute_gramian, normalize, regularize
+from ._matrix import Matrix, PSDMatrix, PSDTensor, is_matrix, is_psd_matrix, is_psd_tensor
 
-__all__ = ["compute_gramian", "Matrix", "PSDMatrix"]
+__all__ = [
+    "compute_gramian",
+    "normalize",
+    "regularize",
+    "Matrix",
+    "PSDMatrix",
+    "PSDTensor",
+    "is_matrix",
+    "is_psd_matrix",
+    "is_psd_tensor",
+]

src/torchjd/_linalg/_gramian.py

@@ -0,0 +1,70 @@
+from typing import Literal, cast, overload
+
+import torch
+from torch import Tensor
+
+from ._matrix import Matrix, PSDMatrix, PSDTensor
+
+
+@overload
+def compute_gramian(t: Tensor) -> PSDMatrix:
+    pass
+
+
+@overload
+def compute_gramian(t: Tensor, contracted_dims: Literal[-1]) -> PSDMatrix:
+    pass
+
+
+@overload
+def compute_gramian(t: Matrix, contracted_dims: Literal[1]) -> PSDMatrix:
+    pass
+
+
+def compute_gramian(t: Tensor, contracted_dims: int = -1) -> PSDTensor:
+    """
+    Computes the `Gramian matrix <https://en.wikipedia.org/wiki/Gram_matrix>`_ of the input.
+
+    `contracted_dims` specifies the number of trailing dimensions to contract. If negative,
+    it indicates the number of leading dimensions to preserve (e.g., ``-1`` preserves the
+    first dimension).
+    """
+
+    contracted_dims = contracted_dims if 0 <= contracted_dims else contracted_dims + t.ndim
+    indices_source = list(range(t.ndim - contracted_dims))
+    indices_dest = list(range(t.ndim - 1, contracted_dims - 1, -1))
+    transposed = t.movedim(indices_source, indices_dest)
+    gramian = torch.tensordot(t, transposed, dims=contracted_dims)
+    return cast(PSDTensor, gramian)
+
+
+def normalize(gramian: PSDMatrix, eps: float) -> PSDMatrix:
+    """
+    Normalizes the gramian `G=AA^T` with respect to the Frobenius norm of `A`.
+
+    If `G=A A^T`, then the Frobenius norm of `A` is the square root of the trace of `G`, i.e., the
+    sqrt of the sum of the diagonal elements. The gramian of the (Frobenius) normalization of `A` is
+    therefore `G` divided by the sum of its diagonal elements.
+    """
+    squared_frobenius_norm = gramian.diagonal().sum()
+    if squared_frobenius_norm < eps:
+        output = torch.zeros_like(gramian)
+    else:
+        output = gramian / squared_frobenius_norm
+    return cast(PSDMatrix, output)
+
+
+def regularize(gramian: PSDMatrix, eps: float) -> PSDMatrix:
+    """
+    Adds a regularization term to the gramian to enforce positive definiteness.
+
+    Because of numerical errors, `gramian` might have slightly negative eigenvalue(s). Adding a
+    regularization term which is a small proportion of the identity matrix ensures that the gramian
+    is positive definite.
+    """
+
+    regularization_matrix = eps * torch.eye(
+        gramian.shape[0], dtype=gramian.dtype, device=gramian.device
+    )
+    output = gramian + regularization_matrix
+    return cast(PSDMatrix, output)

src/torchjd/_linalg/_matrix.py

@@ -0,0 +1,40 @@
+from typing import TypeGuard
+
+from torch import Tensor
+
+# Note: we're using classes and inherittance instead of NewType because it's possible to have
+# multiple inherittance but there is no type intersection. However, these classes should never be
+# instantiated: they're only used for static type checking.
+
+
+class Matrix(Tensor):
+    """Tensor with exactly 2 dimensions."""
+
+
+class PSDTensor(Tensor):
+    """
+    Tensor representing a quadratic form. The first half of its dimensions matches the reversed
+    second half of its dimensions (e.g. shape=[4, 3, 3, 4]), and its reshaping into a matrix should
+    be positive semi-definite.
+    """
+
+
+class PSDMatrix(PSDTensor, Matrix):
+    """Positive semi-definite matrix."""
+
+
+def is_matrix(t: Tensor) -> TypeGuard[Matrix]:
+    return t.ndim == 2
+
+
+def is_psd_tensor(t: Tensor) -> TypeGuard[PSDTensor]:
+    half_dim = t.ndim // 2
+    return t.ndim % 2 == 0 and t.shape[:half_dim] == t.shape[: half_dim - 1 : -1]
+    # We do not check that t is PSD as it is expensive, but this must be checked in the tests of
+    # every function that uses this TypeGuard by using `assert_is_psd_tensor`.
+
+
+def is_psd_matrix(t: Tensor) -> TypeGuard[PSDMatrix]:
+    return t.ndim == 2 and t.shape[0] == t.shape[1]
+    # We do not check that t is PSD as it is expensive, but this must be checked in the tests of
+    # every function that uses this TypeGuard, by using `assert_is_psd_matrix`.

src/torchjd/_linalg/gramian.py

@@ -2,7 +2,7 @@
 
 from torch import Tensor, nn
 
-from torchjd._linalg import Matrix, PSDMatrix, compute_gramian
+from torchjd._linalg import Matrix, PSDMatrix, compute_gramian, is_matrix
 
 from ._weighting_bases import Weighting
 
@@ -18,20 +18,19 @@ def __init__(self):
 
     @staticmethod
     def _check_is_matrix(matrix: Tensor) -> None:
-        if len(matrix.shape) != 2:
+        if not is_matrix(matrix):
             raise ValueError(
                 "Parameter `matrix` should be a tensor of dimension 2. Found `matrix.shape = "
                 f"{matrix.shape}`."
             )
 
     @abstractmethod
-    def forward(self, matrix: Tensor) -> Tensor:
+    def forward(self, matrix: Matrix) -> Tensor:
         """Computes the aggregation from the input matrix."""
 
-    # Override to make type hints and documentation more specific
     def __call__(self, matrix: Tensor) -> Tensor:
         """Computes the aggregation from the input matrix and applies all registered hooks."""
-
+        Aggregator._check_is_matrix(matrix)
         return super().__call__(matrix)
 
     def __repr__(self) -> str:
@@ -54,7 +53,7 @@ def __init__(self, weighting: Weighting[Matrix]):
         self.weighting = weighting
 
     @staticmethod
-    def combine(matrix: Tensor, weights: Tensor) -> Tensor:
+    def combine(matrix: Matrix, weights: Tensor) -> Tensor:
         """
         Aggregates a matrix by making a linear combination of its rows, using the provided vector of
         weights.
@@ -63,8 +62,7 @@ def combine(matrix: Tensor, weights: Tensor) -> Tensor:
         vector = weights @ matrix
         return vector
 
-    def forward(self, matrix: Tensor) -> Tensor:
-        self._check_is_matrix(matrix)
+    def forward(self, matrix: Matrix) -> Tensor:
         weights = self.weighting(matrix)
         vector = self.combine(matrix, weights)
         return vector

src/torchjd/_linalg/matrix.py

@@ -12,8 +12,9 @@
 import torch
 from torch import Tensor
 
+from torchjd._linalg import normalize
+
 from ._aggregator_bases import GramianWeightedAggregator
-from ._utils.gramian import normalize
 from ._utils.non_differentiable import raise_non_differentiable_error
 
 

src/torchjd/aggregation/_aggregator_bases.py

@@ -28,6 +28,8 @@
 import torch
 from torch import Tensor
 
+from torchjd._linalg import Matrix
+
 from ._aggregator_bases import Aggregator
 from ._sum import SumWeighting
 from ._utils.non_differentiable import raise_non_differentiable_error
@@ -56,7 +58,7 @@ def __init__(self, pref_vector: Tensor | None = None):
         # This prevents computing gradients that can be very wrong.
         self.register_full_backward_pre_hook(raise_non_differentiable_error)
 
-    def forward(self, matrix: Tensor) -> Tensor:
+    def forward(self, matrix: Matrix) -> Tensor:
         weights = self.weighting(matrix)
         units = torch.nan_to_num((matrix / (matrix.norm(dim=1)).unsqueeze(1)), 0.0)
         best_direction = torch.linalg.pinv(units) @ weights

src/torchjd/aggregation/_cagrad.py

@@ -2,12 +2,11 @@
 
 from torch import Tensor
 
-from torchjd._linalg import PSDMatrix
+from torchjd._linalg import PSDMatrix, normalize, regularize
 
 from ._aggregator_bases import GramianWeightedAggregator
 from ._mean import MeanWeighting
 from ._utils.dual_cone import project_weights
-from ._utils.gramian import normalize, regularize
 from ._utils.non_differentiable import raise_non_differentiable_error
 from ._utils.pref_vector import pref_vector_to_str_suffix, pref_vector_to_weighting
 from ._weighting_bases import Weighting