Merge pull request #3532 from adowling2/pyomo-doe-fix-initialization

Pyomo.DoE: Corrected initialization when using only lower diagonal of FIM

Author: blnicho

pyomo/contrib/doe/doe.py

@@ -49,6 +49,19 @@
 
 from pyomo.opt import SolverStatus
 
+# This small and positive tolerance is used when checking
+# if the prior is negative definite or approximately
+# indefinite. It is defined as a tolerance here to ensure
+# consistency between the code below and the tests. The
+# user should not need to adjust it.
+_SMALL_TOLERANCE_DEFINITENESS = 1e-6
+
+# This small and positive tolerance is used to check
+# the FIM is approximately symmetric. It is defined as
+# a tolerance here to ensure consistency between the code
+# below and the tests. The user should not need to adjust it.
+_SMALL_TOLERANCE_SYMMETRY = 1e-6
+
 
 class ObjectiveLib(Enum):
     determinant = "determinant"
@@ -309,7 +322,31 @@ def run_doe(self, model=None, results_file=None):
                 (len(model.parameter_names), len(model.parameter_names))
             )
 
-            L_vals_sq = np.linalg.cholesky(fim_np)
+            # Need to compute the full FIM before initializing the Cholesky factorization
+            if self.only_compute_fim_lower:
+                fim_np = fim_np + fim_np.T - np.diag(np.diag(fim_np))
+
+            # Check if the FIM is positive definite
+            # If not, add jitter to the diagonal
+            # to ensure positive definiteness
+            min_eig = np.min(np.linalg.eigvals(fim_np))
+
+            if min_eig < _SMALL_TOLERANCE_DEFINITENESS:
+                # Raise the minimum eigenvalue to at least _SMALL_TOLERANCE_DEFINITENESS
+                jitter = np.min(
+                    [
+                        -min_eig + _SMALL_TOLERANCE_DEFINITENESS,
+                        _SMALL_TOLERANCE_DEFINITENESS,
+                    ]
+                )
+            else:
+                # No jitter needed
+                jitter = 0
+
+            # Add jitter to the diagonal to ensure positive definiteness
+            L_vals_sq = np.linalg.cholesky(
+                fim_np + jitter * np.eye(len(model.parameter_names))
+            )
             for i, c in enumerate(model.parameter_names):
                 for j, d in enumerate(model.parameter_names):
                     model.L[c, d].value = L_vals_sq[i, j]
@@ -1338,7 +1375,28 @@ def check_model_FIM(self, model=None, FIM=None):
                 )
             )
 
-        self.logger.info("FIM provided matches expected dimensions from model.")
+        # Compute the eigenvalues of the FIM
+        evals = np.linalg.eigvals(FIM)
+
+        # Check if the FIM is positive definite
+        if np.min(evals) < -_SMALL_TOLERANCE_DEFINITENESS:
+            raise ValueError(
+                "FIM provided is not positive definite. It has one or more negative eigenvalue(s) less than -{:.1e}".format(
+                    _SMALL_TOLERANCE_DEFINITENESS
+                )
+            )
+
+        # Check if the FIM is symmetric
+        if not np.allclose(FIM, FIM.T, atol=_SMALL_TOLERANCE_SYMMETRY):
+            raise ValueError(
+                "FIM provided is not symmetric using absolute tolerance {}".format(
+                    _SMALL_TOLERANCE_SYMMETRY
+                )
+            )
+
+        self.logger.info(
+            "FIM provided matches expected dimensions from model and is approximately positive (semi) definite."
+        )
 
     # Check the jacobian shape against what is expected from the model.
     def check_model_jac(self, jac=None):

pyomo/contrib/doe/tests/test_doe_errors.py

@@ -165,6 +165,55 @@ def test_reactor_check_bad_prior_size(self):
         ):
             doe_obj.create_doe_model()
 
+    def test_reactor_check_bad_prior_negative_eigenvalue(self):
+        from pyomo.contrib.doe.doe import _SMALL_TOLERANCE_DEFINITENESS
+
+        fd_method = "central"
+        obj_used = "trace"
+        flag_val = 0  # Value for faulty model build mode - 0: full model
+
+        prior_FIM = -np.ones((4, 4))
+
+        experiment = FullReactorExperiment(data_ex, 10, 3)
+
+        DoE_args = get_standard_args(experiment, fd_method, obj_used, flag_val)
+        DoE_args['prior_FIM'] = prior_FIM
+
+        doe_obj = DesignOfExperiments(**DoE_args)
+
+        with self.assertRaisesRegex(
+            ValueError,
+            r"FIM provided is not positive definite. It has one or more negative eigenvalue\(s\) less than -{:.1e}".format(
+                _SMALL_TOLERANCE_DEFINITENESS
+            ),
+        ):
+            doe_obj.create_doe_model()
+
+    def test_reactor_check_bad_prior_not_symmetric(self):
+        from pyomo.contrib.doe.doe import _SMALL_TOLERANCE_SYMMETRY
+
+        fd_method = "central"
+        obj_used = "trace"
+        flag_val = 0  # Value for faulty model build mode - 0: full model
+
+        prior_FIM = np.zeros((4, 4))
+        prior_FIM[0, 1] = 1e-3
+
+        experiment = FullReactorExperiment(data_ex, 10, 3)
+
+        DoE_args = get_standard_args(experiment, fd_method, obj_used, flag_val)
+        DoE_args['prior_FIM'] = prior_FIM
+
+        doe_obj = DesignOfExperiments(**DoE_args)
+
+        with self.assertRaisesRegex(
+            ValueError,
+            "FIM provided is not symmetric using absolute tolerance {}".format(
+                _SMALL_TOLERANCE_SYMMETRY
+            ),
+        ):
+            doe_obj.create_doe_model()
+
     def test_reactor_check_bad_jacobian_init_size(self):
         fd_method = "central"
         obj_used = "trace"

pyomo/contrib/doe/tests/test_doe_solve.py

@@ -230,6 +230,37 @@ def test_reactor_obj_cholesky_solve(self):
         # Make sure FIM and Q.T @ sigma_inv @ Q are close (alternate definition of FIM)
         self.assertTrue(np.all(np.isclose(FIM, Q.T @ sigma_inv @ Q)))
 
+    def test_reactor_obj_cholesky_solve_bad_prior(self):
+
+        from pyomo.contrib.doe.doe import _SMALL_TOLERANCE_DEFINITENESS
+
+        fd_method = "central"
+        obj_used = "determinant"
+
+        experiment = FullReactorExperiment(data_ex, 10, 3)
+
+        DoE_args = get_standard_args(experiment, fd_method, obj_used)
+
+        # Specify a prior that is slightly negative definite
+        # Because it is less than the tolerance, it should be adjusted to be positive definite
+        # No error should be thrown
+        DoE_args['prior_FIM'] = -(_SMALL_TOLERANCE_DEFINITENESS / 100) * np.eye(4)
+
+        doe_obj = DesignOfExperiments(**DoE_args)
+
+        doe_obj.run_doe()
+
+        self.assertEqual(doe_obj.results["Solver Status"], "ok")
+
+        # assert that Q, F, and L are the same.
+        FIM, Q, L, sigma_inv = get_FIM_Q_L(doe_obj=doe_obj)
+
+        # Since Cholesky is used, there is comparison for FIM and L.T @ L
+        self.assertTrue(np.all(np.isclose(FIM, L @ L.T)))
+
+        # Make sure FIM and Q.T @ sigma_inv @ Q are close (alternate definition of FIM)
+        self.assertTrue(np.all(np.isclose(FIM, Q.T @ sigma_inv @ Q)))
+
     # This test ensure that compute FIM runs without error using the
     # `sequential` option with central finite differences
     def test_compute_FIM_seq_centr(self):