Regularization Selection Fixes (#79)

* regression test: parametric ROM + regselect

* delay test case processing

* fix unit tests related to test case TypeErrors

* ROM.fit() calls model.fit(), not _fit_solver()

* BayesianROM regselect: catch LinAlg errors and warnings

* regselect: reset solver; fix regselect for interpolatory models

* literature: Gkimisis 2025

* v0.5.15 -> v0.5.16

Author: shanemcq18

docs/literature.bib

@@ -666,3 +666,14 @@ @article{kang2025semiconductor
     doi = {10.48550/arXiv.2504.03990},
     category = {application}
 }
+
+@article{gkimisis2025spatiallylocalized,
+    title = {Non-intrusive reduced-order modeling for dynamical systems with spatially localized features},
+    author = {Leonidas Gkimisis and Nicole Aretz and Marco Tezzele and Thomas Richter and Peter Benner and Karen E. Willcox},
+    journal = {Computer Methods in Applied Mechanics and Engineering},
+    volume = {444},
+    pages = {118115},
+    year = {2025},
+    doi = {10.1016/j.cma.2025.118115},
+    category = {method}
+}

docs/source/opinf/changelog.md

@@ -5,6 +5,17 @@
 New versions may introduce substantial new features or API adjustments.
 :::
 
+## Version 0.5.16
+
+Backend improvements to the regularization selection procedure.
+
+- `ROM.fit()` calls `ROM.model.fit()` instead of `ROM.model._fit_solver()`, which works better for inheritance.
+- `BayesianROM.fit_regselect_*()` catches errors and suppresses warnings from the least-squares solver.
+- Fixed a bug for interpolatory parametric models where `fit_regselect_*()` would only update the regularizer for the overall `model.solver`, not for each of the individual solvers for the models to be interpolated.
+- Fixed a bug related to test cases being processed before the basis was initialized.
+- Removed unnecessary abstract methods from `models.mono._base._OpInfModel`.
+- Small updates to the literature page.
+
 ## Version 0.5.15
 
 - Improvement to `fit_regselect_*()` so that the regularization does not have to be initialized before fitting the model. This fixes a longstanding chicken/egg problem and makes using `fit_regselect_*()` much less cumbersome.

src/opinf/__init__.py

@@ -7,7 +7,7 @@
     https://github.com/Willcox-Research-Group/rom-operator-inference-Python3
 """
 
-__version__ = "0.5.15"
+__version__ = "0.5.16"
 
 from . import (
     basis,

src/opinf/lstsq/_base.py

@@ -287,6 +287,10 @@ def _load_dict(hf, attr="options"):
             options[key] = None if value == "NULL" else value
         return options
 
+    def reset(self) -> None:
+        """Reset the solver by deleting data matrices."""
+        SolverTemplate.__init__(self)
+
     def copy(self):
         """Make a copy of the solver."""
         return copy.deepcopy(self)

src/opinf/lstsq/_tikhonov.py

@@ -113,6 +113,11 @@ def regresidual(self, Ohat: np.ndarray) -> np.ndarray:
         raise NotImplementedError  # pragma: no cover
 
     # Persistence -------------------------------------------------------------
+    def reset(self) -> None:
+        """Reset the solver by deleting data matrices and the regularizer."""
+        super().reset()
+        self.regularizer = None
+
     def _save(self, savefile, overwrite=False, extras=tuple()):
         """Serialize the solver, saving it in HDF5 format.
         The model can be recovered with the :meth:`_load()` class method.

src/opinf/models/mono/_base.py

@@ -391,31 +391,6 @@ def _check_is_trained(self):
 
     # Fitting -----------------------------------------------------------------
     @abc.abstractmethod
-    def _process_fit_arguments(self, *args, **kwargs):
-        """Prepare training data, validate and set dimensions, etc."""
-        pass
-
-    @abc.abstractmethod
-    def _assemble_data_matrix(self, *args, **kwargs):
-        """Construct the Operator Inference data matrix."""
-        pass
-
-    @abc.abstractmethod
-    def _extract_operators(self, *args, **kwargs):
-        """Unpack the Operator Inference solution, the operator matrix."""
-        pass
-
-    @abc.abstractmethod
-    def _fit_solver(self, *args, **kwargs):
-        """Initialize the regression solver."""
-        pass
-
-    @abc.abstractmethod
-    def refit(self, *args, **kwargs):
-        """Solve the regression and unpack the results."""
-        pass
-
-    @abc.abstractmethod
-    def fit(self, *args, **kwargs):
+    def fit(self, *args, **kwargs):  # pragma: no cover
         """Learn model operators from data."""
         pass

src/opinf/roms/_base.py

@@ -433,22 +433,27 @@ def _process_training_data(
 
         return states, lhs, inputs
 
-    def _process_test_cases(self, test_cases, TestCaseClass):
-        if test_cases is not None:
-            if isinstance(test_cases, TestCaseClass):
-                test_cases = [test_cases]
-            processed_test_cases = []
-            for tcase in test_cases:
-                if not isinstance(tcase, TestCaseClass):
-                    raise TypeError(
-                        "test cases must be "
-                        f"'utils.{TestCaseClass.__name__}' objects"
-                    )
-                processed_test_cases.append(
-                    tcase.copy(self.encode(tcase.initial_conditions))
+    def _process_test_cases(
+        self,
+        test_cases,
+        TestCaseClass: utils._gridsearch._RegTest,
+    ):
+        if test_cases is None:
+            return []
+
+        if isinstance(test_cases, TestCaseClass):
+            test_cases = [test_cases]
+        processed_test_cases = []
+        for tcase in test_cases:
+            if not isinstance(tcase, TestCaseClass):
+                raise TypeError(
+                    "test cases must be "
+                    f"'utils.{TestCaseClass.__name__}' objects"
                 )
-            return processed_test_cases
-        return []
+            processed_test_cases.append(
+                tcase.copy(self.encode(tcase.initial_conditions))
+            )
+        return processed_test_cases
 
     def _get_stability_limits(self, states, stability_margin):
         shifts = [np.mean(Q, axis=1).reshape((-1, 1)) for Q in states]
@@ -466,18 +471,25 @@ def _get_stability_limits(self, states, stability_margin):
                 limits[ell] = np.inf
         return shifts, limits
 
-    def _fit_solver(
+    def _fit_model(
         self,
         parameters,
         states,
         lhs,
         inputs,
-        fit_transformer,
-        fit_basis,
+        fit_transformer: bool,
+        fit_basis: bool,
+        solver_only: bool = False,
     ):
         """Process the training data and fit the model solver.
         Returns the processed training data.
 
+        Parameters
+        ----------
+        solver_only : bool
+            If ``True``, call ``self.model._fit_solver()`` instead of
+            ``self.model.fit()``. This is useful for regularization selection.
+            If ``False`` (default), call ``self.model.fit()``.
         """
         self._check_fit_args(lhs=lhs, inputs=inputs)
         if parameters is None:
@@ -495,11 +507,17 @@ def _fit_solver(
 
         if parameters is None:
             inputdata = None if inputs is None else np.hstack(inputs)
-            self.model._fit_solver(
-                np.hstack(states), np.hstack(lhs), inputdata
-            )
+            if solver_only:
+                self.model._fit_solver(
+                    np.hstack(states), np.hstack(lhs), inputdata
+                )
+            else:
+                self.model.fit(np.hstack(states), np.hstack(lhs), inputdata)
         else:
-            self.model._fit_solver(parameters, states, lhs, inputs)
+            if solver_only:
+                self.model._fit_solver(parameters, states, lhs, inputs)
+            else:
+                self.model.fit(parameters, states, lhs, inputs)
 
         return states
 
@@ -649,23 +667,29 @@ def fit_regselect_continuous(
             raise ValueError("argument 'test_time_length' must be nonnegative")
         if regularizer_factory is None:
             regularizer_factory = _identity
-        processed_test_cases = self._process_test_cases(
-            test_cases, utils.ContinuousRegTest
-        )
+
+        # Reset the solver (in case the basis dimension changed between calls).
+        interp = modutils.is_interpolatory(self.model)
+        if hasattr(self.model, "reset"):
+            self.model.solver.reset()
 
         # Fit the model for the first time.
-        self._fit_solver(
+        self._fit_model(
             parameters=parameters,
             states=states,
             lhs=ddts,
             inputs=inputs,
             fit_transformer=fit_transformer,
             fit_basis=fit_basis,
+            solver_only=True,
         )
 
         # Set up the regularization selection.
         states_ = [self.encode(Q) for Q in states]
         shifts, limits = self._get_stability_limits(states_, stability_margin)
+        processed_test_cases = self._process_test_cases(
+            test_cases, utils.ContinuousRegTest
+        )
 
         def unstable(_Q, ell, size):
             """Return ``True`` if the solution is unstable."""
@@ -695,7 +719,12 @@ def unstable(_Q, ell, size):
 
         def update_model(reg_params):
             """Reset the regularizer and refit the model operators."""
-            self.model.solver.regularizer = regularizer_factory(reg_params)
+            reg = regularizer_factory(reg_params)
+            if interp:
+                for solver in self.model.solvers:
+                    solver.regularizer = reg
+            else:
+                self.model.solver.regularizer = reg
             self.model.refit()
 
         def training_error(reg_params):
@@ -861,22 +890,28 @@ def fit_regselect_discrete(
                     )
         if regularizer_factory is None:
             regularizer_factory = _identity
-        processed_test_cases = self._process_test_cases(
-            test_cases, utils.DiscreteRegTest
-        )
+
+        # Reset the solver (in case the basis dimension changed between calls).
+        interp = modutils.is_interpolatory(self.model)
+        if hasattr(self.model, "reset"):
+            self.model.solver.reset()
 
         # Fit the model for the first time and get compressed training data.
-        states_ = self._fit_solver(
+        states_ = self._fit_model(
             parameters=parameters,
             states=states,
             lhs=None,
             inputs=inputs,
             fit_transformer=fit_transformer,
             fit_basis=fit_basis,
+            solver_only=True,
         )
 
         # Set up the regularization selection.
         shifts, limits = self._get_stability_limits(states_, stability_margin)
+        processed_test_cases = self._process_test_cases(
+            test_cases, utils.DiscreteRegTest
+        )
 
         def unstable(_Q, ell):
             """Return ``True`` if the solution is unstable."""
@@ -897,7 +932,12 @@ def unstable(_Q, ell):
 
         def update_model(reg_params):
             """Reset the regularizer and refit the model operators."""
-            self.model.solver.regularizer = regularizer_factory(reg_params)
+            reg = regularizer_factory(reg_params)
+            if interp:
+                for solver in self.model.solvers:
+                    solver.regularizer = reg
+            else:
+                self.model.solver.regularizer = reg
             self.model.refit()
 
         def training_error(reg_params):

src/opinf/roms/_bayes.py

@@ -8,8 +8,8 @@
 
 import warnings
 import numpy as np
-import scipy.linalg
 import scipy.stats
+import scipy.linalg
 
 from .. import errors, lstsq, post, utils
 from ..models import _utils as modutils
@@ -224,12 +224,11 @@ def posterior(self) -> OperatorPosterior:
 
     def _initialize_posterior(self):
         """Set the operator posterior if numerically possible."""
-        means, precisions = self.model.solver.posterior()
         try:
+            means, precisions = self.model.solver.posterior()
             self.__posterior = OperatorPosterior(means, precisions)
-        except np.linalg.LinAlgError as ex:
-            if ex.args[0] == "Matrix is not positive definite":
-                self.__posterior = None
+        except np.linalg.LinAlgError:
+            self.__posterior = None
 
     def draw_operators(self):
         """Set the :attr:`model` operators to a new random draw from the
@@ -292,23 +291,26 @@ def fit_regselect_continuous(
             raise ValueError("argument 'test_time_length' must be nonnegative")
         if regularizer_factory is None:
             regularizer_factory = _identity
-        processed_test_cases = self._process_test_cases(
-            test_cases, utils.ContinuousRegTest
-        )
 
         # Fit the model for the first time.
-        self._fit_solver(
+        if hasattr(self.model.solver, "reset"):
+            self.model.solver.reset()
+        self._fit_model(
             parameters=parameters,
             states=states,
             lhs=ddts,
             inputs=inputs,
             fit_transformer=fit_transformer,
             fit_basis=fit_basis,
+            solver_only=True,
         )
 
         # Set up the regularization selection.
         states_ = [self.encode(Q) for Q in states]
         shifts, limits = self._get_stability_limits(states_, stability_margin)
+        processed_test_cases = self._process_test_cases(
+            test_cases, utils.ContinuousRegTest
+        )
 
         def unstable(_Q, ell, size):
             """Return ``True`` if the solution is unstable."""
@@ -339,7 +341,9 @@ def unstable(_Q, ell, size):
         def update_model(reg_params):
             """Reset the regularizer and refit the model operators."""
             self.model.solver.regularizer = regularizer_factory(reg_params)
-            self._initialize_posterior()
+            with warnings.catch_warnings():
+                warnings.simplefilter("ignore", scipy.linalg.LinAlgWarning)
+                self._initialize_posterior()
 
         def training_error(reg_params):
             """Compute the training error for a single regularization
@@ -439,22 +443,25 @@ def fit_regselect_discrete(
                     )
         if regularizer_factory is None:
             regularizer_factory = _identity
-        processed_test_cases = self._process_test_cases(
-            test_cases, utils.DiscreteRegTest
-        )
 
         # Fit the model for the first time.
-        states_ = self._fit_solver(
+        if hasattr(self.model.solver, "reset"):
+            self.model.solver.reset()
+        states_ = self._fit_model(
             parameters=parameters,
             states=states,
             lhs=None,
             inputs=inputs,
             fit_transformer=fit_transformer,
             fit_basis=fit_basis,
+            solver_only=True,
         )
 
         # Set up the regularization selection.
         shifts, limits = self._get_stability_limits(states_, stability_margin)
+        processed_test_cases = self._process_test_cases(
+            test_cases, utils.DiscreteRegTest
+        )
 
         def unstable(_Q, ell):
             """Return ``True`` if the solution is unstable."""
@@ -476,7 +483,9 @@ def unstable(_Q, ell):
         def update_model(reg_params):
             """Reset the regularizer and refit the model operators."""
             self.model.solver.regularizer = regularizer_factory(reg_params)
-            self._initialize_posterior()
+            with warnings.catch_warnings():
+                warnings.simplefilter("ignore", scipy.linalg.LinAlgWarning)
+                self._initialize_posterior()
 
         def training_error(reg_params):
             """Compute the mean training error for a single regularization