cvxpy port

README.md

@@ -25,17 +25,30 @@ conda install -n base -c conda-forge ipopt=3.11.1 pkg-config
      - If you have any installation issues, [see the docs](https://docs.astral.sh/uv/getting-started/installation/) for troubleshooting
 2. Install optimizer dependencies.
      - On MacOS: `brew install ipopt pkg-config`
-     - Linux: `sudo apt-get install coinor-libipopt-dev pkg-config`
+     - Linux: `sudo apt-get install coinor-libipopt-dev libblas-dev liblapack-dev pkg-config`
      - Windows: See note above if you have not previously installed Ipopt or JAX
 3. Synchronize the environment
 ```bash
   uv sync
 ```
-4. Finally, you will want to make sure you activate the python environment each time you use it.
+4. **Linux only — obtain the `ipopt` executable.** On Linux the `coinor-libipopt-dev` package provides the Ipopt *library* (needed to build `cyipopt` during `uv sync`) but **not** the `ipopt` command-line executable, which Pyomo's `SolverFactory("ipopt")` invokes to solve the models. The simplest way to get a prebuilt binary is via the IDAES extensions:
+```bash
+  uv run --with idaes-pse idaes get-extensions
+```
+   - This installs `ipopt` (and related solvers) into `~/.idaes/bin`. Add that directory to your `PATH` so Pyomo can find it, e.g. `export PATH="$HOME/.idaes/bin:$PATH"`.
+   - On MacOS (`brew install ipopt`) and Windows (`conda install ... ipopt`) this step is unnecessary because those installs already include the `ipopt` executable.
+5. Finally, you will want to make sure you activate the python environment each time you use it.
    - In VS Code you can activate the default environment with `>Python: Select Interpreter` to be the `.venv` local to the directory
    - If the debugger isn't working in that case, sometimes setting the vscode `terminal.integrated.shellIntegration.enabled: true` in the settings can help
    - Outside of vscode, a simple, platform specific CLI line will [activate .venv](https://docs.python.org/3/tutorial/venv.html#creating-virtual-environments) in your terminal.
 
+## Optimization solvers
+The models use two solver paths:
+
+- **Convex / DCP models** (e.g. asset pricing, `asset_pricing_matern_cvxpy.py`) are solved with [`cvxpy`](https://www.cvxpy.org/) using open-source solvers — **Clarabel, OSQP, SCS, and HiGHS**. These require **no extra setup**: `cvxpy` bundles all four as dependencies, so `uv sync` installs them automatically (HiGHS arrives via `highspy`).
+- **Nonconvex / NLP models** (neoclassical growth, neoclassical human capital, optimal advertising, concave-convex growth — the `*_matern_cvxpy.py` files) are solved through `cvxpy`'s DNLP interface (`prob.solve(nlp=True, ...)`), which hands the smooth nonlinear program to a nonlinear solver. The default is **UNO**, via [`unopy`](https://pypi.org/project/unopy/) — a self-contained binary wheel that statically bundles `libuno`, so `uv sync` installs it with no extra setup. This includes the concave-convex model: its exact complementarity (MPCC) reformulation uses a flat warm start with the marginal product seeded on the *active* production branch (consistent with the complementarity constraints), and UNO resolves the bistable threshold structure. **IPOPT** (via `cyipopt`, also in-process) stays selectable. No solver binary on `PATH` is needed for this path. The figure/table scripts take `--implementation cvxpy|pyomo` (default `cvxpy`), including the concave-convex **threshold** figure.
+- **Pyomo path (alternative).** The same nonconvex models also ship a `pyomo` implementation (e.g. `neoclassical_growth_matern.py`, selected with `--implementation pyomo`) that drives the `ipopt` *binary* on `PATH` — see step 4 above (the IDAES `ipopt` executable on Linux; `brew`/`conda` on MacOS/Windows). UNO is used only through `cvxpy` (the `unopy` wheel above), so no standalone solver binary needs to be installed for it.
+
 **Troubleshooting**:
 - If you receive JAX errors about DLL load failures, you may need to update [https://learn.microsoft.com/en-us/cpp/windows/latest-supported-vc-redist?view=msvc-170#visual-studio-2015-2017-2019-and-2022](https://learn.microsoft.com/en-us/cpp/windows/latest-supported-vc-redist?view=msvc-170#visual-studio-2015-2017-2019-and-2022)
 - See notes above on Windows Ipopt installation challenges

asset_pricing_matern_cvxpy.py

@@ -0,0 +1,114 @@
+import time
+import jax
+import jax.numpy as jnp
+import numpy as np
+import cvxpy as cp
+import jsonargparse
+from jax import config
+from kernels import integrated_matern_kernel_matrices
+from asset_pricing_benchmark import mu_f_array
+from typing import List, Optional
+
+config.update("jax_enable_x64", True)
+
+# Pyomo-free DCP port of asset_pricing_matern. The model is a convex QP:
+#   minimize  alpha' K alpha
+#   s.t.      K alpha == r (mu_0 + K_tilde alpha) - x,   mu_0 >= 0
+# solver_type selects an open-source cvxpy backend. OSQP (native QP) is the
+# default; CLARABEL is the more robust choice for ill-conditioned kernels
+# (large nu/rho/N), where OSQP's ADMM can report infeasibility.
+SOLVER_OPTIONS = {
+    "OSQP": (cp.OSQP, dict(eps_abs=1e-12, eps_rel=1e-12, max_iter=5000)),
+    "CLARABEL": (
+        cp.CLARABEL,
+        dict(tol_gap_abs=1e-12, tol_gap_rel=1e-12, tol_feas=1e-12),
+    ),
+    "SCS": (cp.SCS, dict(eps=1e-9, max_iters=20000)),
+    "HIGHS": (cp.HIGHS, dict(primal_feasibility_tolerance=1e-4)),
+}
+
+
+def asset_pricing_matern_cvxpy(
+    r: float = 0.1,
+    c: float = 0.02,
+    g: float = -0.2,
+    x_0: float = 0.01,
+    nu: float = 0.5,
+    sigma: float = 1.0,
+    rho: float = 10,
+    solver_type: str = "OSQP",
+    train_T: float = 40.0,
+    train_points: int = 41,
+    test_T: float = 50.0,
+    test_points: int = 100,
+    train_points_list: Optional[List[float]] = None,
+    verbose: bool = False,
+):
+    # if passing in `train_points_list` then doesn't use a grid.  Otherwise, uses linspace
+    if train_points_list is None:
+        train_data = jnp.linspace(0, train_T, train_points)
+    else:
+        train_data = jnp.array(train_points_list)
+    test_data = jnp.linspace(0, test_T, test_points)
+
+    # Construct kernel matrices.  cvxpy needs numpy (not jax) arrays.
+    N = len(train_data)
+    K, K_tilde = integrated_matern_kernel_matrices(
+        train_data, train_data, nu, sigma, rho
+    )
+    K = np.asarray((K + K.T) / 2)  # symmetrize -> exactly PSD for quad_form
+    K_tilde = np.asarray(K_tilde)
+    x = (x_0 + c / g) * np.exp(g * np.asarray(train_data)) - c / g
+
+    # Solve the QP.  psd_wrap asserts the (provably PSD) Gram matrix K so cvxpy
+    # skips its O(N^3) PSD certification, which both dominates canonicalization
+    # time and fails to converge for ill-conditioned K.
+    alpha_mu = cp.Variable(N)
+    mu_0 = cp.Variable(nonneg=True)
+    prob = cp.Problem(
+        cp.Minimize(cp.quad_form(alpha_mu, cp.psd_wrap(K))),
+        [K @ alpha_mu == r * (mu_0 + K_tilde @ alpha_mu) - x],
+    )
+    solver, options = SOLVER_OPTIONS[solver_type]
+    start = time.perf_counter()
+    prob.solve(solver=solver, verbose=verbose, **options)
+    print(f"elapsed solve(s) = {time.perf_counter() - start}")
+    if prob.status not in ("optimal", "optimal_inaccurate"):
+        print(f"solver status: {prob.status}")
+
+    alpha_mu = jnp.array(alpha_mu.value)
+    mu_0 = float(mu_0.value)
+
+    # Interpolator using training data
+    @jax.jit
+    def kernel_solution(test_data):
+        # pointwise comparison test_data to train_data
+        _, K_tilde_test = integrated_matern_kernel_matrices(
+            test_data, train_data, nu, sigma, rho
+        )
+        mu_test = mu_0 + K_tilde_test @ alpha_mu
+        return mu_test
+
+    # Generate test_data and compare to the benchmark
+    mu_benchmark = mu_f_array(test_data, c, g, r, x_0)
+    mu_test = kernel_solution(test_data)
+
+    mu_rel_error = jnp.abs(mu_benchmark - mu_test) / mu_benchmark
+    print(
+        f"solve_time(s) = {prob.solver_stats.solve_time}, E(|rel_error(p)|) = {mu_rel_error.mean()}"
+    )
+    return {
+        "t_train": train_data,
+        "t_test": test_data,
+        "p_test": mu_test,
+        "p_benchmark": mu_benchmark,
+        "p_rel_error": mu_rel_error,
+        "alpha": alpha_mu,
+        "p_0": mu_0,
+        "solve_time": prob.solver_stats.solve_time,
+        "kernel_solution": kernel_solution,  # interpolator
+    }
+
+
+if __name__ == "__main__":
+    jsonargparse.CLI(asset_pricing_matern_cvxpy)

figures_asset_pricing.py

@@ -1,7 +1,9 @@
 import jax.numpy as jnp
 import matplotlib.pyplot as plt
 import os
+import jsonargparse
 from asset_pricing_matern import asset_pricing_matern
+from asset_pricing_matern_cvxpy import asset_pricing_matern_cvxpy
 
 from mpl_toolkits.axes_grid1.inset_locator import (
     zoomed_inset_axes,
@@ -126,8 +128,16 @@ def plot_asset_pricing(
     plt.savefig(output_path, format="pdf")
 
 
-# Plots with various parameters
-sol_matern = asset_pricing_matern()
-plot_asset_pricing(
-    sol_matern, "figures/asset_pricing_contiguous.pdf"
-)
+# Plots with various parameters.  implementation selects the solve backend:
+# "cvxpy" (default, open-source DCP/QP solvers) or "pyomo" (Ipopt).
+def main(implementation: str = "cvxpy"):
+    solve = {
+        "cvxpy": asset_pricing_matern_cvxpy,
+        "pyomo": asset_pricing_matern,
+    }[implementation]
+    sol_matern = solve()
+    plot_asset_pricing(sol_matern, "figures/asset_pricing_contiguous.pdf")
+
+
+if __name__ == "__main__":
+    jsonargparse.CLI(main)