Pyomo.DoE: fix trace/Cholesky initialization consistency by adowling2 · Pull Request #3867 · Pyomo/pyomo

adowling2 · 2026-03-04T01:00:59Z

Fixes # .

Summary/Motivation:

This PR ~~adds advanced debugging/inspection controls for pyomo.contrib.doe.DesignOfExperiments.run_doe() while preserving default behavior for typical users~~ fixes a trace-objective initialization inconsistency that can produce large residuals in Cholesky-related constraints at the start of the final NLP solve.

The motivating workflow showed:

large early infeasibilities in optimize mode,
~~inability to use max_iter=0 for final-model probing because scenario solves reused the same solver options,~~
and inconsistent initialization among fim, L, fim_inv, L_inv, and cov_trace.

Changes proposed in this PR:

~~Added an advanced run_config (ConfigBlock/dict) argument to run_doe(model=None, results_file=None, run_config=None):~~
- ~~scenario_solver_options~~
- ~~final_solver_options~~
- ~~final_solve~~
- ~~inspection.enabled~~
- ~~inspection.top_constraints~~
~~Implemented scoped solver option application:~~
- ~~scenario-generation and square-initialization solves can use different options from the final optimization solve.~~
~~Added optional structured residual reporting utility:~~
- ~~reports top violated active constraints with fields:~~
  - ~~constraint_name~~
  - ~~body~~
  - ~~lower_bound~~
  - ~~upper_bound~~
  - ~~violation~~
  - ~~constraint_type~~
~~Added assemble-and-inspect path (final_solve=False) to inspect assembled NLP state without running final optimization.~~
Fixed Cholesky jitter computation and added robust helper usage.
Added _initialize_cholesky_from_fim() to re-synchronize L, L_inv, fim_inv, and cov_trace from current FIM immediately before final solve / inspection.
~~Updated docstrings and added advanced usage documentation + example path in reactor_example.py.~~
Added/updated regression tests in test_doe_initialization.py for jitter behavior and FIM synchronization.

Backward compatibility:

~~Default usage (run_doe() with no run_config) remains unchanged.~~
~~Advanced functionality is available only when run_config is provided.~~

Testing:

I ran:

~~pytest -q pyomo/contrib/doe/tests/test_doe_debug.py~~
- ~~Result: 8 passed~~
~~pytest -q pyomo/contrib/doe/tests/test_doe_solve.py -k "rooney_biegler_fd_central_solve"~~
- ~~Result: 1 passed, 17 deselected~~
python3 -m unittest -q pyomo.contrib.doe.tests.test_doe_initialization
- Result: 3 skipped in the current interpreter because SciPy/numpy are unavailable there
python3 -m py_compile pyomo/contrib/doe/doe.py pyomo/contrib/doe/tests/test_doe_initialization.py
- passed

New/updated tests cover:

negative-jitter helper behavior,
zero-jitter helper behavior,
Cholesky/FIM inverse consistency after initialization in trace mode.

Outstanding Tasks (before removing WIP):

~~Finalize proposed interface for advanced initialization diagnostics~~
~~Decide the best way to document the proposed interface. Are examples sufficient or should we build out the Pyomo documentation?~~
~~Run black~~
~~@adowling2 carefully reviews proposed changes~~

Why the strikethought?

The original version of this PR included a new interface to run_doe() that makes debugging the initialization and, more broadly, solver failures easier. We decided that the new run_doe() interface added too much complexity and exposed some inconsistencies in other parts of the Pyomo solver interface. (This led to a real-time Codex experiment documented in #3871.) To expedite this current PR, I decided to remove these contributions and strike through the above parts of the PR text. These contributions are still available in the git history. Thus, this current PR only contributes a targeted patch to the Pyomo FIM initialization.

Legal Acknowledgement

By contributing to this software project, I have read the contribution guide and agree to the following terms and conditions for my contribution:

I agree my contributions are submitted under the BSD license.
I represent I am authorized to make the contributions and grant the license. If my employer has rights to intellectual property that includes these contributions, I represent that I have received permission to make contributions and grant the required license on behalf of that employer.

adowling2 · 2026-03-04T01:05:15Z

@smondal13 @sscini @slilonfe5 @snarasi2 @djlaky I found a small bug in how we initialize the Cholesky factorization when the FIM at the initial point needs "inertia correction". This PR fixes the initialization inconsistency. This PR also proposes an interface to allow for more advanced solver options when diagnosing initialization issues for the Pyomo.DoE (dynamic) optimization formulation.

Edit: This is a work in progress right now. I am not 100% happy with the proposed changes.

adowling2 · 2026-05-12T19:46:28Z

@sscini @blnicho @mrmundt @smondal13 @slilonfe5 This is ready for review. This PR is now narrowly targeted at fixing a Pyomo.DoE FIM initialization issue.

adowling2 · 2026-05-13T11:08:34Z

Test failures appear to be unrelated to this PR: GAMS solvers (Windows) and NEOS solver (all)

slilonfe5 · 2026-05-13T13:23:29Z

Test failures appear to be unrelated to this PR: GAMS solvers (Windows) and NEOS solver (all)

These failures are global. The Pyomo team mentioned that they are working on resolving the issues.

smondal13 · 2026-05-13T13:51:42Z

+                        model.fim_inv[c, d].value = fim_inv_vals[i, j]
+
+        if hasattr(model, "cov_trace"):
+            fim_inv_np = np.linalg.inv(fim_pd)


Since you are making the fim PD by adding the jitter, I believe the fim_pd will probably be non-singular. But if it is singular for some edge cases, then isn't it better to use the pseudo-inverse here just to be safe?

I like this suggestion. I am going to switch to pinv for most of the initialization code.

smondal13 · 2026-05-13T14:45:01Z

+
+
+@unittest.skipIf(
+    not (numpy_available and scipy_available),


I think scipy is not required for the following tests. So, it may unnecessarily skip tests. Only greybox and multi-experiment initialization with LHS require scipy as far as I can remember.

…/pyomo into fix-doe-initialization

mrmundt · 2026-05-20T15:47:15Z

+    solver = _SolverData()
+
+
+class _MutatingRecordingSolver:


Are there any tests in any of the other files that exercise initialization using a real solver? (I do see changes in other tests in other files, but it's unclear to me if any of them are actually exercising the initialization code.)

If not, can you please add a test that uses a real solver / real result as well?

blnicho · 2026-05-19T22:13:37Z

+        if (
+            hasattr(model, "dummy_obj")
+            and model.dummy_obj.active
+            and hasattr(model, "fim")
+            and len(list(model.parameter_names)) == 2
+        ):


You know the attributes of the model being tested with this dummy solver so is this if-statement really necessary?

blnicho · 2026-05-19T22:19:39Z

+        return _FakeResult()
+
+
+class _TwoParamExperiment:


Should this inherit from the Experiment class?

blnicho · 2026-05-19T22:29:56Z

+        jitter = doe_obj._compute_cholesky_jitter(min_eig)
+        self.assertAlmostEqual(
+            jitter, _SMALL_TOLERANCE_DEFINITENESS - min_eig, places=14
+        )


I know this is a simple method but in general we don't want to repeat the calculation implementation in the test and instead check against a hard-coded value. Checking up to 14 places also seems excessive. I would also recommend adding a test verifying that the return value is non-negative which is the intention of this test mentioned in the docstring.

blnicho · 2026-05-20T13:47:53Z

+        Finite-difference mode passed to DoE (for example, ``"central"``).
+    obj_used : str
+        Objective option name (for example, ``"trace"`` or ``"determinant"``).


Suggested change

Finite-difference mode passed to DoE (for example, ``"central"``).

obj_used : str

Objective option name (for example, ``"trace"`` or ``"determinant"``).

Finite-difference mode passed to DoE (for example, "central").

obj_used : str

Objective option name (for example, "trace" or "determinant").

I would recommend removing the extra formatting around the quotations

blnicho · 2026-05-20T13:52:56Z

    # Test that we can properly
    # set the inputs for the
    # Grey Box object
    def test_set_inputs(self):
+        """Verify packed upper-triangular inputs reconstruct the full FIM."""


This docstring doesn't really match the comment above the method describing the intended purpose. Maybe the comment should just be used as the docstring instead?

blnicho · 2026-05-20T16:57:03Z

+    def test_D_opt_greybox_build_with_triangular_fim(self):
+        """Determinant grey-box build accepts triangular FIM storage."""
+        objective_option = "determinant"
+        doe_obj, grey_box_object = make_greybox_and_doe_objects(
+            objective_option=objective_option
+        )
+
+        self.assertTrue(doe_obj.only_compute_fim_lower)
+        # Grey-box objective construction should not depend on the Cholesky
+        # build path. This fixture may still carry the shared default flag, so
+        # we do not treat it as a grey-box invariant here.
+        self.assertTrue(doe_obj.use_grey_box)
+
+        doe_obj.create_grey_box_objective_function()
+
+        self.assertTrue(hasattr(doe_obj.model.obj_cons, "egb_fim_block"))
+        self.assertIn("log-D-opt", list(grey_box_object.output_names()))
+        self.assertIsInstance(doe_obj.model.objective, pyo.Objective)
+        self.assertEqual(doe_obj.model.objective.sense, pyo.maximize)


I don't really understand the need for this test. It's checking for certain attributes and types but not checking any values or calculations. Aren't there separate unit tests verifying that objective function is being constructed with the right type and sense? Could these assertions be added to the previous test, test_D_opt_greybox_build?

The comment about "Grey-box objective construction not depending on the Cholesky build path" seems like it's trying to highlight something subtle in the implementation but I'm not following the subtlety based on the surrounding assertions.

blnicho · 2026-05-20T17:28:13Z

+    def _compute_cholesky_jitter(self, min_eig):
+        """
+        Compute diagonal regularization for Cholesky initialization.
+
+        Parameters
+        ----------
+        min_eig: float
+            Minimum eigenvalue of the current FIM estimate.
+
+        Returns
+        -------
+        float
+            Nonnegative diagonal shift needed so the minimum eigenvalue
+            is at least ``_SMALL_TOLERANCE_DEFINITENESS``.
+        """
+        return max(0.0, _SMALL_TOLERANCE_DEFINITENESS - float(min_eig))


Is there a reason to have this as a method on the DesignOfExperiments class instead of a standalone utility function? It isn't using anything from the class instance.

Also, it's a 1-line method that is only being used once on L555. Why did you even abstract this out as a separate method?

blnicho · 2026-05-20T17:47:13Z

+        if not hasattr(model, "L"):
+            return


Could you add a comment explaining why this is a valid check to skip the rest of this method?

blnicho · 2026-05-20T17:56:21Z

+        jitter = self._compute_cholesky_jitter(min_eig)
+        fim_pd = fim_np + jitter * np.eye(len(model.parameter_names))
+
+        L_vals = np.linalg.cholesky(fim_pd)
+        for i, c in enumerate(model.parameter_names):
+            for j, d in enumerate(model.parameter_names):
+                if i >= j:
+                    model.L[c, d].value = L_vals[i, j]
+                else:
+                    model.L[c, d].value = 0.0
+
+        if hasattr(model, "L_inv"):
+            L_inv_vals = np.linalg.inv(L_vals)
+            for i, c in enumerate(model.parameter_names):
+                for j, d in enumerate(model.parameter_names):
+                    if i >= j:
+                        model.L_inv[c, d].value = L_inv_vals[i, j]
+                    else:
+                        model.L_inv[c, d].value = 0.0


This implementation is slightly different from the implementation where these values are originally calculated on L1443-1465. Is there a reason for the difference? Can these be combined to avoid duplication?

blnicho · 2026-05-20T18:01:01Z

+        # Inputs store one triangular half
+        # of a symmetric FIM. Reconstruct
+        # the full symmetric matrix here,
+        # consistent with manuscript S5.


What is "manuscript S5"?

adowling2 added 5 commits March 2, 2026 14:39

Added prototype capabilities for debugging Pyomo.DoE initial points

5e99540

Updated jitter calculation for FIM

191ef45

Updated FIM initialization

5470b1d

Improved API for run_doe

d871e08

Added missing doc strings

d588397

adowling2 marked this pull request as draft March 4, 2026 01:01

adowling2 added 6 commits March 3, 2026 20:13

Add doc strings to the tests

6797a90

Adds comments

6f20842

Dedicated example for Pyomo.DoE initial point debugging

d889cc6

Ran Black

64a22ee

Updated tests

35e4b85

Improved test documentation

2ab8041

adowling2 marked this pull request as ready for review March 4, 2026 23:07

adowling2 changed the title ~~WIP: Pyomo.DoE: add advanced run_doe debug config and fix trace/Cholesky initialization consistency~~ Pyomo.DoE: add advanced run_doe debug config and fix trace/Cholesky initialization consistency Mar 4, 2026

adowling2 added 2 commits March 16, 2026 12:17

Additional edits to for GreyBox in Pyomo.DoE

5f5d0cc

Continued to debug GreyBox

7497df6

blnicho self-requested a review March 17, 2026 19:15

blnicho added the pyomo.DoE label Mar 18, 2026

github-project-automation Bot added this to ParmEst & Pyomo.DoE Development Mar 18, 2026

github-project-automation Bot moved this to Todo in ParmEst & Pyomo.DoE Development Mar 18, 2026

adowling2 and others added 5 commits April 2, 2026 21:31

Update to avoid NaN

a0005ec

Merge branch 'main' into fix-doe-initialization

b9ca381

Merge branch 'main' into fix-doe-initialization

a1770cb

Removed new run_doe() interface

f0fa351

Changed test file name

dc40a41

Ran black and typos

ec3d9c4

adowling2 changed the title ~~Pyomo.DoE: add advanced run_doe debug config and fix trace/Cholesky initialization consistency~~ Pyomo.DoE: fix trace/Cholesky initialization consistency May 13, 2026

Fixed two tests

eabb32f

smondal13 reviewed May 13, 2026

View reviewed changes

adowling2 added 6 commits May 14, 2026 20:45

Merge branch 'main' into fix-doe-initialization

3c9e291

Addressed review feedback

8b2f487

Merge branch 'fix-doe-initialization' of https://github.com/adowling2…

ff56322

…/pyomo into fix-doe-initialization

Added edge case / error tests

16345be

Added pseudo-A test for GreyBox

96973f3

Merge branch 'main' into fix-doe-initialization

9cc3111

mrmundt reviewed May 20, 2026

View reviewed changes

blnicho requested changes May 20, 2026

View reviewed changes

Conversation

adowling2 commented Mar 4, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Fixes # .

Summary/Motivation:

Changes proposed in this PR:

Backward compatibility:

Testing:

Outstanding Tasks (before removing WIP):

Why the strikethought?

Legal Acknowledgement

Uh oh!

adowling2 commented Mar 4, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

adowling2 commented May 12, 2026

Uh oh!

adowling2 commented May 13, 2026

Uh oh!

slilonfe5 commented May 13, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

6 participants

adowling2 commented Mar 4, 2026 •

edited

Loading

adowling2 commented Mar 4, 2026 •

edited

Loading