From 26bda239c0ccf69df3b5c6a3774df55ee24c6058 Mon Sep 17 00:00:00 2001 From: Aaditya Chandrasekhar Date: Wed, 12 Nov 2025 21:00:35 -0600 Subject: [PATCH] more tests on lagrange multipliers --- .gitignore | 233 +++++++++++++++++++ mmapy/mma_test.py | 560 +++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 788 insertions(+), 5 deletions(-) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..01ab9f1 --- /dev/null +++ b/.gitignore @@ -0,0 +1,233 @@ +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[codz] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py.cover +.hypothesis/ +.pytest_cache/ +cover/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 +db.sqlite3-journal + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +.pybuilder/ +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py + +# pyenv +# For a library or package, you might want to ignore these files since the code is +# intended to run in multiple environments; otherwise, check them in: +# .python-version + +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock + +# UV +# Similar to Pipfile.lock, it is generally recommended to include uv.lock in version control. +# This is especially recommended for binary packages to ensure reproducibility, and is more +# commonly ignored for libraries. +#uv.lock + +# poetry +# Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control. +# This is especially recommended for binary packages to ensure reproducibility, and is more +# commonly ignored for libraries. +# https://python-poetry.org/docs/basic-usage/#commit-your-poetrylock-file-to-version-control +#poetry.lock +#poetry.toml + +# pdm +# Similar to Pipfile.lock, it is generally recommended to include pdm.lock in version control. +# pdm recommends including project-wide configuration in pdm.toml, but excluding .pdm-python. +# https://pdm-project.org/en/latest/usage/project/#working-with-version-control +#pdm.lock +#pdm.toml +.pdm-python +.pdm-build/ + +# pixi +# Similar to Pipfile.lock, it is generally recommended to include pixi.lock in version control. +#pixi.lock +# Pixi creates a virtual environment in the .pixi directory, just like venv module creates one +# in the .venv directory. It is recommended not to include this directory in version control. +.pixi + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow and github.com/pdm-project/pdm +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.envrc +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# Pyre type checker +.pyre/ + +# pytype static type analyzer +.pytype/ + +# Cython debug symbols +cython_debug/ + +# PyCharm +# JetBrains specific template is maintained in a separate JetBrains.gitignore that can +# be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore +# and can be added to the global gitignore or merged into this file. For a more nuclear +# option (not recommended) you can uncomment the following to ignore the entire idea folder. +#.idea/ + +# Abstra +# Abstra is an AI-powered process automation framework. +# Ignore directories containing user credentials, local state, and settings. +# Learn more at https://abstra.io/docs +.abstra/ + +# Visual Studio Code +# Visual Studio Code specific template is maintained in a separate VisualStudioCode.gitignore +# that can be found at https://github.com/github/gitignore/blob/main/Global/VisualStudioCode.gitignore +# and can be added to the global gitignore or merged into this file. However, if you prefer, +# you could uncomment the following to ignore the entire vscode folder +# .vscode/ + +# Ruff stuff: +.ruff_cache/ + +# PyPI configuration file +.pypirc + +# Cursor +# Cursor is an AI-powered code editor. `.cursorignore` specifies files/directories to +# exclude from AI features like autocomplete and code analysis. Recommended for sensitive data +# refer to https://docs.cursor.com/context/ignore-files +.cursorignore +.cursorindexingignore + +# Marimo +marimo/_static/ +marimo/_lsp/ +__marimo__/ + + +# Ignore image files +*.jpg +*.jpeg +*.gif +*.bmp +*.svg +*.png + +# Ignore video files +*.mp4 +*.mov +*.avi +*.mkv +*.wmv + +*.DS_Store + + +*.xlsx +*.ods +*.docx*.doc +*.pptx*.ppt +*.pkl +*.h5 \ No newline at end of file diff --git a/mmapy/mma_test.py b/mmapy/mma_test.py index 3234d43..e687e21 100644 --- a/mmapy/mma_test.py +++ b/mmapy/mma_test.py @@ -275,7 +275,6 @@ def objfn(xy): while not mma_state.is_converged: objective, grad_obj = objective_fn(mma_state.x) - print(f"{mma_state.epoch} obj, {objective}") constr, grad_cons = dummy_constraint(mma_state.x) mma_state = mma.update_mma( mma_state, mma_params, objective, grad_obj, constr, grad_cons @@ -388,8 +387,6 @@ def objfn(x): mma_state, mma_params, objective, grad_obj, constr, grad_cons ) - print(f"epoch {mma_state.epoch} , obj {objective[0]:.2E}") - np.testing.assert_allclose(objective, 0.0, atol=1e-2) np.testing.assert_allclose(mma_state.x[0, 0], 3.0, rtol=1e-1) np.testing.assert_allclose(mma_state.x[1, 0], 0.5, rtol=1e-1) @@ -442,8 +439,6 @@ def objfn(x): mma_state, mma_params, objective, grad_obj, constr, grad_cons ) - print(f"epoch {mma_state.epoch} , obj {objective[0]:.2E}") - np.testing.assert_allclose(objective, 0.0, atol=1e-2) np.testing.assert_allclose(mma_state.x[0, 0], 1.0, rtol=1e-1) np.testing.assert_allclose(mma_state.x[1, 0], 3, rtol=1e-1) @@ -707,6 +702,561 @@ def objfn(x): msg="Inactive upper bound should have ~0 multiplier", ) + def test_stationarity_with_active_general_constraint(self): + """Test stationarity when a general constraint is active. + + Theory: At optimum with active constraint: + ∇f(x*) + λ*∇g(x*) + bound terms = 0 + + Problem: + min x^2 + y^2 + s.t. x + y ≥ 1.5 (equivalently: -x - y + 1.5 ≤ 0) + -2 ≤ x, y ≤ 2 + + Expected: + x* = y* = 0.75 (on constraint boundary) + Constraint active: -0.75 - 0.75 + 1.5 = 0 + λ > 0 (positive multiplier) + Stationarity: ∇f + λ*∇g ≈ 0 + """ + num_design_var = 2 + num_cons = 1 + + def constraint_fn(x): + def confn(x): + # Constraint: x + y ≥ 1.5 → -x - y + 1.5 ≤ 0 + return -x[0, 0] - x[1, 0] + 1.5 + + c, dc = value_and_grad(confn)(x) + return ( + np.array(c).reshape((num_cons, 1)), + np.array(dc).reshape((num_cons, num_design_var)), + ) + + def objective_fn(xy): + def objfn(xy): + return xy[0, 0] ** 2 + xy[1, 0] ** 2 + + obj, grad_obj = value_and_grad(objfn)(xy) + return obj.reshape((-1)), grad_obj.reshape((-1, 1)) + + design_var = np.array([[0.5], [0.5]]) + lower_bound = -2.0 * np.ones((num_design_var, 1)) + upper_bound = 2.0 * np.ones((num_design_var, 1)) + + mma_params = mma.MMAParams( + max_iter=200, + kkt_tol=1e-5, + step_tol=1e-5, + move_limit=1e-2, + num_design_var=num_design_var, + num_cons=num_cons, + lower_bound=lower_bound, + upper_bound=upper_bound, + ) + mma_state = mma.init_mma(design_var, mma_params) + + while not mma_state.is_converged: + objective, grad_obj = objective_fn(mma_state.x) + constr, grad_cons = constraint_fn(mma_state.x) + mma_state, lm = mma.update_mma( + mma_state, + mma_params, + objective, + grad_obj, + constr, + grad_cons, + return_lagrange_multipliers=True, + ) + + x_opt = mma_state.x + + # Recompute gradients at optimum + _, grad_f = objective_fn(x_opt) + c_val, grad_c = constraint_fn(x_opt) + + # Check constraint is active + self.assertLess( + abs(c_val[0, 0]), 0.1, msg="Constraint should be approximately active" + ) + + # Check multiplier is positive + self.assertGreater( + lm.general_constraints[0, 0], + 0.1, + msg="Active constraint should have positive multiplier", + ) + + # Compute stationarity: ∇f + λ*∇g + bound terms + grad_lagrangian = grad_f.copy() + grad_lagrangian += lm.general_constraints[0, 0] * grad_c.T + grad_lagrangian -= lm.lower_bounds + grad_lagrangian += lm.upper_bounds + + stationarity_norm = np.linalg.norm(grad_lagrangian) + + # Verify stationarity + self.assertLess( + stationarity_norm, + 1e-3, + msg=f"Stationarity violated: ||∇L|| = {stationarity_norm}", + ) + + def test_stationarity_at_lower_bound(self): + """Test stationarity when solution is at lower bound. + + Theory: When x_j is at lower bound: + (∇f)_j - ξ_j ≈ 0, with ξ_j > 0 + + Problem: + min (x - 0.5)^2 + (y - 0.5)^2 + s.t. 1 ≤ x, y ≤ 10 + + Expected: + x* = y* = 1 (at lower bounds) + ∇f(x*) = [1.0, 1.0] + ξ ≈ ∇f for stationarity + """ + num_design_var = 2 + num_cons = 1 + + def dummy_constraint(x): + del x + return ( + np.array([0.0]).reshape((num_cons, 1)), + np.array([0.0, 0.0]).reshape((num_cons, num_design_var)), + ) + + def objective_fn(xy): + def objfn(xy): + return (xy[0, 0] - 0.5) ** 2 + (xy[1, 0] - 0.5) ** 2 + + obj, grad_obj = value_and_grad(objfn)(xy) + return obj.reshape((-1)), grad_obj.reshape((-1, 1)) + + design_var = np.array([[5.0], [5.0]]) + lower_bound = np.array([[1.0], [1.0]]) + upper_bound = np.array([[10.0], [10.0]]) + + mma_params = mma.MMAParams( + max_iter=200, + kkt_tol=1e-6, + step_tol=1e-6, + move_limit=1e-2, + num_design_var=num_design_var, + num_cons=num_cons, + lower_bound=lower_bound, + upper_bound=upper_bound, + ) + mma_state = mma.init_mma(design_var, mma_params) + + while not mma_state.is_converged: + objective, grad_obj = objective_fn(mma_state.x) + constr, grad_cons = dummy_constraint(mma_state.x) + mma_state, lm = mma.update_mma( + mma_state, + mma_params, + objective, + grad_obj, + constr, + grad_cons, + return_lagrange_multipliers=True, + ) + + x_opt = mma_state.x + + # Verify at lower bounds + self.assertAlmostEqual(x_opt[0, 0], 1.0, places=2) + self.assertAlmostEqual(x_opt[1, 0], 1.0, places=2) + + # Recompute gradient at optimum + _, grad_f = objective_fn(x_opt) + + # For bound constraints: ∇f - ξ ≈ 0 (when at lower bound) + residual = grad_f - lm.lower_bounds + residual_norm = np.linalg.norm(residual) + + # Both components should have positive multipliers + self.assertGreater(lm.lower_bounds[0, 0], 0.01) + self.assertGreater(lm.lower_bounds[1, 0], 0.01) + + # Stationarity: ∇f - ξ should be small + self.assertLess( + residual_norm, + 0.5, + msg=f"Stationarity at bounds violated: ||∇f - ξ|| = {residual_norm}", + ) + + def test_stationarity_at_upper_bound(self): + """Test stationarity when solution is at upper bound. + + Theory: When x_j is at upper bound: + (∇f)_j + η_j ≈ 0, with η_j > 0 + + Problem: + min (x - 10.5)^2 + (y - 10.5)^2 + s.t. 1 ≤ x, y ≤ 10 + + Expected: + x* = y* = 10 (at upper bounds) + ∇f(x*) = [-1.0, -1.0] + η ≈ -∇f for stationarity + """ + num_design_var = 2 + num_cons = 1 + + def dummy_constraint(x): + del x + return ( + np.array([0.0]).reshape((num_cons, 1)), + np.array([0.0, 0.0]).reshape((num_cons, num_design_var)), + ) + + def objective_fn(xy): + def objfn(xy): + return (xy[0, 0] - 10.5) ** 2 + (xy[1, 0] - 10.5) ** 2 + + obj, grad_obj = value_and_grad(objfn)(xy) + return obj.reshape((-1)), grad_obj.reshape((-1, 1)) + + design_var = np.array([[5.0], [5.0]]) + lower_bound = np.array([[1.0], [1.0]]) + upper_bound = np.array([[10.0], [10.0]]) + + mma_params = mma.MMAParams( + max_iter=200, + kkt_tol=1e-6, + step_tol=1e-6, + move_limit=1e-2, + num_design_var=num_design_var, + num_cons=num_cons, + lower_bound=lower_bound, + upper_bound=upper_bound, + ) + mma_state = mma.init_mma(design_var, mma_params) + + while not mma_state.is_converged: + objective, grad_obj = objective_fn(mma_state.x) + constr, grad_cons = dummy_constraint(mma_state.x) + mma_state, lm = mma.update_mma( + mma_state, + mma_params, + objective, + grad_obj, + constr, + grad_cons, + return_lagrange_multipliers=True, + ) + + x_opt = mma_state.x + + # Verify at upper bounds + self.assertAlmostEqual(x_opt[0, 0], 10.0, places=2) + self.assertAlmostEqual(x_opt[1, 0], 10.0, places=2) + + # Recompute gradient at optimum + _, grad_f = objective_fn(x_opt) + + # For upper bound: ∇f + η ≈ 0 + residual = grad_f + lm.upper_bounds + residual_norm = np.linalg.norm(residual) + + # Both components should have positive multipliers + self.assertGreater(lm.upper_bounds[0, 0], 0.01) + self.assertGreater(lm.upper_bounds[1, 0], 0.01) + + # Stationarity: ∇f + η should be small + self.assertLess( + residual_norm, + 0.5, + msg=f"Stationarity at upper bounds violated: ||∇f + η|| = {residual_norm}", + ) + + def test_stationarity_mixed_bounds(self): + """Test stationarity with variables at different bounds. + + Problem: + min (x - 0.5)^2 + (y - 10.5)^2 + (z - 5)^2 + s.t. 1 ≤ x, y, z ≤ 10 + + Expected: + x* = 1 (at lower), y* = 10 (at upper), z* = 5 (interior) + ξ_x > 0, η_y > 0, ξ_z ≈ 0, η_z ≈ 0 + """ + num_design_var = 3 + num_cons = 1 + + def dummy_constraint(x): + del x + return ( + np.array([0.0]).reshape((num_cons, 1)), + np.array([0.0, 0.0, 0.0]).reshape((num_cons, num_design_var)), + ) + + def objective_fn(xyz): + def objfn(xyz): + return (xyz[0, 0] - 0.5) ** 2 + (xyz[1, 0] - 10.5) ** 2 + (xyz[2, 0] - 5.0) ** 2 + + obj, grad_obj = value_and_grad(objfn)(xyz) + return obj.reshape((-1)), grad_obj.reshape((-1, 1)) + + design_var = np.array([[5.0], [5.0], [5.0]]) + lower_bound = np.array([[1.0], [1.0], [1.0]]) + upper_bound = np.array([[10.0], [10.0], [10.0]]) + + mma_params = mma.MMAParams( + max_iter=200, + kkt_tol=1e-6, + step_tol=1e-6, + move_limit=1e-2, + num_design_var=num_design_var, + num_cons=num_cons, + lower_bound=lower_bound, + upper_bound=upper_bound, + ) + mma_state = mma.init_mma(design_var, mma_params) + + while not mma_state.is_converged: + objective, grad_obj = objective_fn(mma_state.x) + constr, grad_cons = dummy_constraint(mma_state.x) + mma_state, lm = mma.update_mma( + mma_state, + mma_params, + objective, + grad_obj, + constr, + grad_cons, + return_lagrange_multipliers=True, + ) + + x_opt = mma_state.x + + # Verify expected bounds + self.assertAlmostEqual(x_opt[0, 0], 1.0, places=2, msg="x should be at lower") + self.assertAlmostEqual(x_opt[1, 0], 10.0, places=2, msg="y should be at upper") + self.assertAlmostEqual(x_opt[2, 0], 5.0, places=2, msg="z should be interior") + + # Check multipliers + self.assertGreater(lm.lower_bounds[0, 0], 0.01, msg="x at lower → ξ_x > 0") + self.assertGreater(lm.upper_bounds[1, 0], 0.01, msg="y at upper → η_y > 0") + self.assertAlmostEqual( + lm.lower_bounds[2, 0], 0.0, places=1, msg="z interior → ξ_z ≈ 0" + ) + self.assertAlmostEqual( + lm.upper_bounds[2, 0], 0.0, places=1, msg="z interior → η_z ≈ 0" + ) + + def test_multiplier_scales_with_objective(self): + """Test that multipliers scale linearly with objective magnitude. + + Theory: For scaled problem min α*f(x) s.t. g(x) ≤ 0 + Solution x* is the same, but λ_scaled = α * λ_original + + Problem: + min α * (x - 0.5)^2 + s.t. 1 ≤ x ≤ 10 + + Test with α = 1, 10, 100 + + Expected: + x* = 1 for all α (same solution) + λ_bound(α) = α * λ_bound(1) + """ + num_design_var = 1 + num_cons = 1 + + def dummy_constraint(x): + del x + return ( + np.array([0.0]).reshape((num_cons, 1)), + np.array([0.0]).reshape((num_cons, num_design_var)), + ) + + def create_objective_fn(scale_factor): + def objective_fn(x): + def objfn(x): + return scale_factor * (x[0, 0] - 0.5) ** 2 + + obj, grad_obj = value_and_grad(objfn)(x) + return obj.reshape((-1)), grad_obj.reshape((-1, 1)) + + return objective_fn + + lower_bound = np.array([[1.0]]) + upper_bound = np.array([[10.0]]) + + scales = [1.0, 10.0, 100.0] + results = [] + + for scale in scales: + design_var = np.array([[5.0]]) + + mma_params = mma.MMAParams( + max_iter=200, + kkt_tol=1e-6, + step_tol=1e-6, + move_limit=1e-2, + num_design_var=num_design_var, + num_cons=num_cons, + lower_bound=lower_bound, + upper_bound=upper_bound, + ) + mma_state = mma.init_mma(design_var, mma_params) + + objective_fn = create_objective_fn(scale) + + while not mma_state.is_converged: + objective, grad_obj = objective_fn(mma_state.x) + constr, grad_cons = dummy_constraint(mma_state.x) + mma_state, lm = mma.update_mma( + mma_state, + mma_params, + objective, + grad_obj, + constr, + grad_cons, + return_lagrange_multipliers=True, + ) + + results.append( + { + "scale": scale, + "x_opt": mma_state.x[0, 0], + "lambda_lower": lm.lower_bounds[0, 0], + } + ) + + # All should converge to same x* + for r in results: + self.assertAlmostEqual( + r["x_opt"], 1.0, places=1, msg="Solution should be same for all scales" + ) + + # Check scaling relationship + base_lambda = results[0]["lambda_lower"] + for i, r in enumerate(results[1:], start=1): + expected_lambda = base_lambda * (r["scale"] / results[0]["scale"]) + actual_lambda = r["lambda_lower"] + relative_error = abs(actual_lambda - expected_lambda) / max(expected_lambda, 1e-6) + + self.assertLess( + relative_error, + 0.2, # 20% tolerance + msg=f"Multiplier should scale with objective. " + f"Expected {expected_lambda:.4f}, got {actual_lambda:.4f}", + ) + + def test_multiplier_scales_inversely_with_constraint(self): + """Test that constraint multipliers scale inversely with constraint scaling. + + Theory: Scaling constraint α*g doesn't change feasible region, + but λ scales as 1/α to maintain stationarity balance. + + Problem: + min x^2 + y^2 + s.t. α*(-x - y + 1.5) ≤ 0 + -2 ≤ x, y ≤ 2 + + Test with α = 1, 2, 5 + + Expected: + x*, y* same for all α + λ(α) * α ≈ constant + """ + num_design_var = 2 + num_cons = 1 + + def create_constraint_fn(scale_factor): + def constraint_fn(x): + def confn(x): + return scale_factor * (-x[0, 0] - x[1, 0] + 1.5) + + c, dc = value_and_grad(confn)(x) + return ( + np.array(c).reshape((num_cons, 1)), + np.array(dc).reshape((num_cons, num_design_var)), + ) + + return constraint_fn + + def objective_fn(xy): + def objfn(xy): + return xy[0, 0] ** 2 + xy[1, 0] ** 2 + + obj, grad_obj = value_and_grad(objfn)(xy) + return obj.reshape((-1)), grad_obj.reshape((-1, 1)) + + lower_bound = -2.0 * np.ones((num_design_var, 1)) + upper_bound = 2.0 * np.ones((num_design_var, 1)) + + scales = [1.0, 2.0, 5.0] + results = [] + + for scale in scales: + design_var = np.array([[0.5], [0.5]]) + + mma_params = mma.MMAParams( + max_iter=200, + kkt_tol=1e-5, + step_tol=1e-5, + move_limit=1e-2, + num_design_var=num_design_var, + num_cons=num_cons, + lower_bound=lower_bound, + upper_bound=upper_bound, + ) + mma_state = mma.init_mma(design_var, mma_params) + + constraint_fn = create_constraint_fn(scale) + + while not mma_state.is_converged: + objective, grad_obj = objective_fn(mma_state.x) + constr, grad_cons = constraint_fn(mma_state.x) + mma_state, lm = mma.update_mma( + mma_state, + mma_params, + objective, + grad_obj, + constr, + grad_cons, + return_lagrange_multipliers=True, + ) + + results.append( + { + "scale": scale, + "x_opt": mma_state.x[0, 0], + "y_opt": mma_state.x[1, 0], + "lambda_constraint": lm.general_constraints[0, 0], + } + ) + + # Solutions should be approximately the same + for i in range(len(results) - 1): + self.assertAlmostEqual( + results[i]["x_opt"], + results[i + 1]["x_opt"], + places=1, + msg="x* should be same for all scalings", + ) + self.assertAlmostEqual( + results[i]["y_opt"], + results[i + 1]["y_opt"], + places=1, + msg="y* should be same for all scalings", + ) + + # Check inverse scaling: λ * scale should be approximately constant + products = [r["lambda_constraint"] * r["scale"] for r in results] + mean_product = np.mean(products) + for prod in products: + relative_deviation = abs(prod - mean_product) / mean_product + self.assertLess( + relative_deviation, + 0.3, # 30% tolerance + msg="λ * scale should be approximately constant", + ) + if __name__ == "__main__": absltest.main()