test: Add Valley-Bowl, optimizer and LR options to trajectories

tests/trajectories/_constants.py

@@ -4,6 +4,7 @@
 from torchjd.aggregation import (
     IMTLG,
     MGDA,
+    Aggregator,
     AlignedMTL,
     CAGrad,
     DualProj,
@@ -18,9 +19,10 @@
     ElementWiseQuadratic,
     HomogenousQuadraticFunction,
     QuadraticFunction,
+    ValleyBowl,
 )
 
-AGGREGATORS = {
+AGGREGATORS: dict[str, Aggregator] = {
     "upgrad": UPGrad(projector=QuadprogProjector(reg_eps=1e-7, norm_eps=1e-9)),
     "mgda": MGDA(),
     "cagrad": CAGrad(c=0.5),
@@ -33,6 +35,10 @@
     "random": Random(),
     "mean": Mean(),
 }
+OPTIMIZERS: dict[str, type[torch.optim.SGD] | type[torch.optim.Adam]] = {
+    "sgd": torch.optim.SGD,
+    "adam": torch.optim.Adam,
+}
 LR_MULTIPLIERS = {
     "upgrad": 1.0,
     "mgda": 2.0,
@@ -88,9 +94,31 @@
         "xlim": (-0.125, 2.625),
         "ylim": (-0.425, 8.925),
     },
+    "VB": {
+        "xlim": (-21.0, 11.0),
+        "ylim": (-21.0, 11.0),
+    },
+    "VB01": {
+        "xlim": (-2.1, 1.1),
+        "ylim": (-21.0, 11.0),
+    },
+}
+
+# Fixed xlim and ylim of the param plot, used instead of limits fitted to the content
+PLOT_PARAMS_LIMS = {
+    "VB": {
+        "xlim": (-10.5, 10.5),
+        "ylim": (-10.5, 10.5),
+    },
+    "VB01": {
+        "xlim": (-10.5, 10.5),
+        "ylim": (-10.5, 10.5),
+    },
 }
 
 OBJECTIVES = {
+    "VB": ValleyBowl(),
+    "VB01": ValleyBowl(scale=0.1),
     "EWQ": ElementWiseQuadratic(2),
     "CQF": QuadraticFunction(
         As=[torch.tensor([[1.0, 0.2], [0.2, 0.05]]), torch.tensor([[3.0, -0.6], [-0.6, 0.2]])],
@@ -103,11 +131,27 @@
     ),
 }
 BASE_LEARNING_RATES = {
+    "VB": 0.001,
+    "VB01": 0.001,
     "EWQ": 0.075,
     "CQF": 0.125,
     "HQF": 0.005,
 }
 INITIAL_POINTS = {
+    "VB": [
+        [-8.5, 7.5],
+        [0.0, 0.0],
+        [9.0, 9.0],
+        [-7.5, -0.5],
+        [9.0, -1.0],
+    ],
+    "VB01": [
+        [-8.5, 7.5],
+        [0.0, 0.0],
+        [9.0, 9.0],
+        [-7.5, -0.5],
+        [9.0, -1.0],
+    ],
     "EWQ": [
         [3.0, -2.0],
         [0.0, -3.0],
@@ -128,6 +172,8 @@
     ],
 }
 N_ITERS = {
+    "VB": 35000,
+    "VB01": 35000,
     "EWQ": 50,
     "CQF": 200,
     "HQF": 100,

tests/trajectories/_objectives.py

@@ -1,6 +1,7 @@
 from __future__ import annotations
 
 from abc import ABC, abstractmethod
+from typing import cast
 
 import torch
 from torch import Tensor
@@ -130,3 +131,59 @@ def __init__(self, n_values: int) -> None:
 
     def __call__(self, w: Tensor) -> Tensor:  # noqa: ARG002
         return torch.zeros(self.n_values)
+
+
+class ValleyBowl(Objective):
+    """Two-objective toy problem from CAGrad (Liu et al., 2021), also used by Nash-MTL
+    (Navon et al., 2022).
+
+    Above :math:`x_2 = 0`, each objective is a log-shaped valley; below it, a quadratic bowl.
+
+    Reference implementations:
+
+    - CAGrad: https://github.com/Cranial-XIX/CAGrad/blob/main/toy.py
+    - Nash-MTL: https://github.com/AvivNavon/nash-mtl/tree/main/experiments/toy
+
+    :param scale: Multiplier for the first objective.
+    :param scale_both_losses: Multiplier applied to both objectives.
+    """
+
+    # Value used for the Pareto set in CAGrad's plotting code (``toy.py``, ``plot_contour``).
+    PARETO_SET_X2 = -8.3552
+
+    def __init__(self, scale: float = 1.0, scale_both_losses: float = 1.0) -> None:
+        super().__init__(n_params=2, n_values=2)
+        self.scale = scale
+        self.scale_both_losses = scale_both_losses
+
+    def __call__(self, x: Tensor) -> Tensor:
+        if x.shape != (2,):
+            raise ValueError("x must have shape (2,).")
+        x1, x2 = x.unbind()
+        valley1 = (0.5 * (-x1 - 7) - torch.tanh(-x2)).abs().clamp_min(5e-6).log() + 6
+        valley2 = (0.5 * (-x1 + 3) + torch.tanh(-x2) + 2).abs().clamp_min(5e-6).log() + 6
+        above = torch.tanh(x2 * 0.5).clamp_min(0)
+        below = torch.tanh(-x2 * 0.5).clamp_min(0)
+        bowl1 = ((-x1 + 7).square() + 0.1 * (-x2 - 8).square()) / 10 - 20
+        bowl2 = ((-x1 - 7).square() + 0.1 * (-x2 - 8).square()) / 10 - 20
+        return (
+            torch.stack(
+                ((valley1 * above + bowl1 * below) * self.scale, valley2 * above + bowl2 * below)
+            )
+            * self.scale_both_losses
+        )
+
+    def jacobian(self, x: Tensor) -> Tensor:
+        return cast(Tensor, torch.func.jacrev(self)(x))
+
+    def __repr__(self) -> str:
+        return f"ValleyBowl(scale={self.scale}, scale_both_losses={self.scale_both_losses})"
+
+    def pareto_set(self, n: int = 1000) -> Tensor:
+        """Sample the reference Pareto set used by the Valley-Bowl example.
+
+        :param n: Number of points to sample along the reference set.
+        """
+        return torch.stack(
+            (torch.linspace(-7.0, 7.0, n), torch.full((n,), self.PARETO_SET_X2)), dim=1
+        )

tests/trajectories/_optimization.py

@@ -14,11 +14,21 @@ def optimize(
     aggregator: Aggregator,
     lr: float,
     n_iters: int,
+    optimizer_cls: type[torch.optim.SGD] | type[torch.optim.Adam] = torch.optim.SGD,
 ) -> tuple[list[Tensor], list[Tensor]]:
+    """Run an aggregator with the selected PyTorch optimizer.
+
+    :param objective: Vector-valued objective to optimize.
+    :param initial_x: Starting point in parameter space.
+    :param aggregator: Aggregator applied to the objective Jacobian.
+    :param lr: Optimizer learning rate.
+    :param n_iters: Number of optimization steps.
+    :param optimizer_cls: Optimizer class to instantiate for this trajectory.
+    """
     xs = []
     ys = []
     x = initial_x.clone().requires_grad_()
-    optimizer = torch.optim.SGD([x], lr=lr)
+    optimizer = optimizer_cls([x], lr=lr)
     for i in range(n_iters):
         print(f"{i + 1}/{n_iters}", end="\r")
         xs.append(x.detach().clone())

tests/trajectories/_plotters.py

@@ -312,6 +312,55 @@ def __call__(self, ax: plt.Axes) -> None:
         )
 
 
+class FilledContourPlotter(Plotter):
+    """
+    Plotter that draws filled contours of the given values on a regular grid extending between the
+    provided coordinates, with thin contour lines at the boundaries between levels.
+
+    :param values: Array of shape ``[n_x, n_y]``, where ``values[i, j]`` is the value at the i-th x
+        coordinate and the j-th y coordinate.
+    :param x_min: Smallest x coordinate of the grid.
+    :param x_max: Largest x coordinate of the grid.
+    :param y_min: Smallest y coordinate of the grid.
+    :param y_max: Largest y coordinate of the grid.
+    :param cmap: Name of the matplotlib colormap.
+    :param alpha: Opacity of the filled contours.
+    :param levels: Number of contour levels.
+    """
+
+    def __init__(
+        self,
+        values: np.ndarray,
+        x_min: float,
+        x_max: float,
+        y_min: float,
+        y_max: float,
+        cmap: str,
+        alpha: float,
+        levels: int = 35,
+    ) -> None:
+        self.values = values
+        self.x_min = x_min
+        self.x_max = x_max
+        self.y_min = y_min
+        self.y_max = y_max
+        self.cmap = cmap
+        self.alpha = alpha
+        self.levels = levels
+
+    def __call__(self, ax: plt.Axes) -> None:
+        xs = np.linspace(self.x_min, self.x_max, self.values.shape[0])
+        ys = np.linspace(self.y_min, self.y_max, self.values.shape[1])
+        ax.contourf(
+            xs, ys, self.values.T, levels=self.levels, cmap=self.cmap, alpha=self.alpha, zorder=0
+        )
+        # Draw the level boundaries explicitly, since vector outputs (e.g. PDF) render adjacent
+        # filled levels without any visible seam.
+        ax.contour(
+            xs, ys, self.values.T, levels=self.levels, cmap=self.cmap, linewidths=0.3, zorder=0
+        )
+
+
 class LimAdjuster(Plotter):
     """Plotter that adjusts the xlim and ylim of the plot to the specified xlim and ylim."""
 
@@ -383,8 +432,8 @@ def __call__(self, ax: plt.Axes) -> None:
         ax.set_title(self.title)
 
 
-class SquareBoxAspectSetter(Plotter):
-    """Plotter that sets a square box aspect."""
+class BoxAspectSetter(Plotter):
+    """Plotter that sets a slightly wider than tall box aspect."""
 
     def __call__(self, ax: plt.Axes) -> None:
-        ax.set_box_aspect(1)
+        ax.set_box_aspect(5 / 6)

tests/trajectories/optimize.py

@@ -3,7 +3,7 @@
 spaces.
 
 Usage:
-    uv run python tests/trajectories/optimize.py <objective> <aggregator>...
+    uv run python tests/trajectories/optimize.py <objective> <aggregator>... [--optimizer sgd|adam] [--lr LR]
 
 Arguments:
     <objective>        The key of the objective function (e.g., EWQ, CQF, HQF).
@@ -28,6 +28,7 @@
     LR_MULTIPLIERS,
     N_ITERS,
     OBJECTIVES,
+    OPTIMIZERS,
 )
 from trajectories._optimization import optimize
 from trajectories._paths import get_params_dir, get_values_dir
@@ -52,7 +53,19 @@ def main() -> None:
         metavar="aggregator",
         help=f"Keys of the aggregators to use. Choices: {list(AGGREGATORS)}",
     )
+    parser.add_argument(
+        "--optimizer",
+        "--optim",
+        choices=list(OPTIMIZERS),
+        default="sgd",
+        help="Optimizer (default: sgd).",
+    )
+    parser.add_argument(
+        "--lr", type=float, default=None, help="Use this learning rate for every aggregator."
+    )
     args = parser.parse_args()
+    if args.lr is not None and (not np.isfinite(args.lr) or args.lr <= 0):
+        parser.error("--lr must be finite and positive")
 
     objective_key = args.objective
     if objective_key not in OBJECTIVES:
@@ -71,8 +84,9 @@ def main() -> None:
     base_lr = BASE_LEARNING_RATES[objective_key]
     for key in LR_MULTIPLIERS:
         mult = lr_multiplier_overrides.get(key, LR_MULTIPLIERS[key])
-        learning_rates[key] = mult * base_lr
+        learning_rates[key] = args.lr if args.lr is not None else mult * base_lr
     n_iters = N_ITERS[objective_key]
+    optimizer_cls = OPTIMIZERS[args.optimizer]
 
     torch.use_deterministic_algorithms(True)
 
@@ -88,6 +102,7 @@ def main() -> None:
         "aggregator_reprs": {key: repr(AGGREGATORS[key]) for key in aggregator_keys},
         "learning_rates": learning_rates,
         "initial_points": initial_points,
+        "optimizer": args.optimizer,
     }
     with open(TRAJECTORIES_RESULTS_DIR / objective_key / "metadata.json", "w") as f:
         json.dump(metadata, f)
@@ -107,7 +122,12 @@ def main() -> None:
 
             initial_x = torch.tensor(initial_point)
             xs, ys = optimize(
-                objective, initial_x=initial_x, aggregator=aggregator, lr=lr, n_iters=n_iters
+                objective,
+                initial_x=initial_x,
+                aggregator=aggregator,
+                lr=lr,
+                n_iters=n_iters,
+                optimizer_cls=optimizer_cls,
             )
 
             xs_list.append(torch.stack(xs))

tests/trajectories/plot_distance_to_pf.py

@@ -20,8 +20,8 @@
 from trajectories._pareto_utils import make_2d_pf_distance_fn
 from trajectories._paths import get_distance_to_pf_plots_dir, get_values_dir
 from trajectories._plotters import (
+    BoxAspectSetter,
     MultiEvolutionPlotter,
-    SquareBoxAspectSetter,
     TitleSetter,
     XAxisLabeller,
     XTicksClearer,
@@ -65,7 +65,7 @@ def main() -> None:
     objective_key = metadata["objective_key"]
     objective = OBJECTIVES[objective_key]
 
-    common_plotter = SquareBoxAspectSetter()
+    common_plotter = BoxAspectSetter()
 
     aggregator_keys = metadata["aggregator_keys"]
     aggregator_to_Y = {key: np.load(values_dir / f"{key}.npy") for key in aggregator_keys}