diff --git a/.github/workflows/portfolio-pgd.yml b/.github/workflows/portfolio-pgd.yml
new file mode 100644
index 00000000..9cebad76
--- /dev/null
+++ b/.github/workflows/portfolio-pgd.yml
@@ -0,0 +1,30 @@
+name: Portfolio PGD experiment
+
+on:
+  pull_request:
+    paths:
+      - "baseline/experiments/portfolio_pgd/**"
+      - ".github/workflows/portfolio-pgd.yml"
+  push:
+    branches: [main]
+    paths:
+      - "baseline/experiments/portfolio_pgd/**"
+      - ".github/workflows/portfolio-pgd.yml"
+
+jobs:
+  validate:
+    runs-on: ubuntu-latest
+    timeout-minutes: 15
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: "3.11"
+          cache: pip
+          cache-dependency-path: baseline/experiments/portfolio_pgd/pyproject.toml
+      - name: Install experiment
+        run: python -m pip install -e './baseline/experiments/portfolio_pgd[notebook]'
+      - name: Run tests
+        run: python baseline/experiments/portfolio_pgd/tests/run_tests.py
+      - name: Execute notebook cells
+        run: python baseline/experiments/portfolio_pgd/scripts/execute_notebooks.py
diff --git a/README.md b/README.md
index ea81c056..b9b46979 100644
--- a/README.md
+++ b/README.md
@@ -248,3 +248,20 @@ correspondence without overloading the Marchenko–Pastur aspect-ratio symbol
 `baseline/` contains no RG intervention. Each optimizer package remains an
 independent experiment that must be evaluated against the same frozen reference
 suite and its own README.
+
+## Portfolio projected-gradient experiment
+
+The independent portfolio-construction experiment is in
+[`baseline/experiments/portfolio_pgd`](baseline/experiments/portfolio_pgd). It validates projected
+gradient descent against exact KKT and SciPy reference solutions, then exercises nonlinear impact,
+factor and sector exposure controls, box constraints, turnover, and gross exposure.
+
+Mac setup and the complete foreground campaign use literal `/tmp`:
+
+```bash
+bash baseline/experiments/portfolio_pgd/scripts/setup_mac.sh
+bash baseline/experiments/portfolio_pgd/scripts/run_portfolio_pgd_experiment.sh
+```
+
+The campaign streams timestamped optimization progress to the terminal and saves the same log plus
+CSV/JSON metrics under `/tmp/portfolio-pgd-runs`.
diff --git a/baseline/experiments/portfolio_pgd/.gitignore b/baseline/experiments/portfolio_pgd/.gitignore
new file mode 100644
index 00000000..e2742770
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/.gitignore
@@ -0,0 +1,9 @@
+.venv/
+__pycache__/
+*.py[cod]
+.pytest_cache/
+.ruff_cache/
+.ipynb_checkpoints/
+build/
+dist/
+*.egg-info/
diff --git a/baseline/experiments/portfolio_pgd/LICENSE b/baseline/experiments/portfolio_pgd/LICENSE
new file mode 100644
index 00000000..f8007416
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/LICENSE
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2026 Calculation Consulting
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/baseline/experiments/portfolio_pgd/README.md b/baseline/experiments/portfolio_pgd/README.md
new file mode 100644
index 00000000..cff881ba
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/README.md
@@ -0,0 +1,178 @@
+# Portfolio PGD Optimizer Experiment
+
+A tested projected-gradient implementation for the portfolio-construction models in
+*Portfolio Construction and Information Flow: Transaction Costs, Constraints, and Forecast
+Persistence*.
+
+This experiment lives inside `CalculatedContent/rg_optimizers` at
+`baseline/experiments/portfolio_pgd`. Runtime environments, logs, metrics, and generated outputs
+are kept under literal `/tmp` on the Mac.
+
+## Complete Mac run from a fresh checkout
+
+```bash
+cd /tmp
+git clone https://github.com/CalculatedContent/rg_optimizers.git
+cd /tmp/rg_optimizers
+
+bash baseline/experiments/portfolio_pgd/scripts/setup_mac.sh
+bash baseline/experiments/portfolio_pgd/scripts/run_portfolio_pgd_experiment.sh
+```
+
+The foreground campaign runner prints every stage and command, streams PGD progress records with
+objective, projected-gradient norm, step size, constraint violation, and elapsed time, and saves the
+same terminal log under:
+
+```text
+/tmp/portfolio-pgd-runs/<run-id>/run.log
+```
+
+Machine-readable output is written beside it under `metrics/`: per-scenario histories, holdings,
+trades, `summary.csv`, and `summary.json`. The runner never backgrounds, detaches, replaces the
+calling shell, or sends process-killing signals.
+
+The package solves
+
+\[
+\min_h\;
+\frac{\lambda}{2}h^\top Vh-\alpha^\top h
++\frac{\theta}{2}(h-h_-)^\top Q(h-h_-)
++c(h-h_-)
+\]
+
+over an intersection of convex portfolio constraints. Here, \(h_-\) is the pre-trade portfolio and
+\(c\) may be a smooth square-root/power-law or smoothed bid-ask cost.
+
+## What is included
+
+- Projected gradient descent with a majorization line search and convergence diagnostics.
+- Exact KKT solver for the equality-constrained quadratic case.
+- Independent SciPy SLSQP benchmark solver. Turnover and gross exposure use exact lifted linear
+  formulations rather than nonsmooth finite differences.
+- Dykstra projection over analytic projectors for:
+  - full-investment, factor-neutrality, and other affine equalities;
+  - sector, industry, and other linear exposure bounds;
+  - long-only/short and per-name bounds;
+  - hard \(L_1\) turnover;
+  - hard gross exposure.
+- Convex power-law impact and smoothed bid-ask transaction costs.
+- Ten automated tests covering gradients, progress callbacks, exact quadratic solutions, nonlinear costs, and
+  realistic institutional constraint intersections.
+- Three detailed Jupyter notebooks.
+
+## Standalone installation
+
+Python 3.10 or later is required.
+
+```bash
+python3 -m venv .venv
+source .venv/bin/activate          # Windows: .venv\Scripts\activate
+python -m pip install --upgrade pip
+python -m pip install -e ".[notebook]"
+```
+
+The optimizer itself requires only NumPy and SciPy. The `notebook` extra adds JupyterLab,
+Matplotlib, and pandas.
+
+## Run the tests
+
+The tests use the standard-library `unittest` runner, so pytest is not required:
+
+```bash
+python tests/run_tests.py
+```
+
+or:
+
+```bash
+python -m unittest discover -s tests -p "test_*.py" -v
+```
+
+## Open the notebooks
+
+```bash
+jupyter notebook notebooks
+```
+
+Run them in order:
+
+1. `01_quadratic_pgd_vs_standard_solvers.ipynb`
+2. `02_nonlinear_transaction_costs.ipynb`
+3. `03_realistic_constraints.ipynb`
+
+Each notebook contains numerical assertions and can be run from a fresh kernel.
+
+## Minimal example
+
+```python
+import numpy as np
+from portfolio_pgd import ConstraintSet, PGDOptions, PortfolioProblem, solve_pgd
+
+n = 20
+rng = np.random.default_rng(7)
+raw = rng.normal(size=(n, n))
+covariance = raw @ raw.T / n + 0.1 * np.eye(n)
+
+problem = PortfolioProblem(
+    alpha=rng.normal(scale=0.02, size=n),
+    covariance=covariance,
+    previous_holdings=np.full(n, 1.0 / n),
+    risk_aversion=2.0,
+    quadratic_cost_matrix=np.ones(n),
+    quadratic_cost_aversion=0.5,
+)
+
+constraints = ConstraintSet(
+    n,
+    equality_matrix=np.ones((1, n)),
+    equality_target=np.array([1.0]),
+    lower_bounds=0.0,
+    upper_bounds=0.10,
+    turnover_limit=0.25,
+    turnover_center=problem.previous_holdings,
+)
+
+result = solve_pgd(problem, constraints, options=PGDOptions(tolerance=1e-8))
+print(result.status, result.utility, result.max_constraint_violation)
+```
+
+## Solver conventions
+
+- The implementation **minimizes negative utility**. `result.utility` is the economically familiar
+  maximized quantity; `result.objective` is its negative.
+- All matrices use the holdings convention \(A_{eq}h=b_{eq}\) and
+  \(A_{ub}h\le b_{ub}\).
+- Turnover is two-way turnover \(\lVert h-h_-\rVert_1\). Divide by two externally if your reporting
+  convention defines one-way turnover.
+- Power-law costs with \(1<p<2\) have singular curvature at zero. For numerical work, set a small
+  positive `epsilon`; the notebooks use `1e-3`.
+- Cardinality, minimum lots, fixed ticket costs, and integer positions are nonconvex and are outside
+  this package's global-convergence guarantees.
+
+## Numerical method
+
+For the quadratic case, the Hessian is
+
+\[
+H=\lambda V+\theta Q\succ0,
+\]
+
+so the solution is unique. The package validates PGD against both the exact KKT system and SciPy's
+SLSQP implementation. For nonlinear smooth convex costs, PGD uses the projected-majorization test
+
+\[
+F(h^+)\le F(h)+\nabla F(h)^\top(h^+-h)+\frac{\lVert h^+-h\rVert^2}{2\eta}
+\]
+
+to backtrack to a safe step.
+
+Projection onto a constraint intersection is performed by Dykstra's algorithm. Every constituent
+set has an analytic Euclidean projector; failure to reach the requested feasibility tolerance raises
+`ProjectionError` instead of silently returning an infeasible portfolio.
+
+## Scope
+
+This is a research-quality reference implementation intended for reproducible experiments and as a
+clear foundation for a production service. A production deployment should additionally integrate
+data validation, sparse matrix representations, warm starts across rebalance dates, monitoring,
+and organization-specific trading/risk controls.
diff --git a/baseline/experiments/portfolio_pgd/VALIDATION.md b/baseline/experiments/portfolio_pgd/VALIDATION.md
new file mode 100644
index 00000000..25c73694
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/VALIDATION.md
@@ -0,0 +1,40 @@
+# Validation record
+
+Validated on 2026-08-21 with Python 3 and the declared NumPy/SciPy runtime dependencies.
+
+## Automated tests
+
+Command:
+
+```bash
+python tests/run_tests.py
+```
+
+Result: **10 tests passed**.
+
+Coverage includes:
+
+- analytic objective and nonlinear-cost gradients versus centered finite differences;
+- live progress-callback initial/final state and final-result consistency;
+- PGD holdings and objective versus the exact quadratic KKT solution;
+- exact KKT versus SciPy SLSQP;
+- nonlinear power-law PGD versus SciPy SLSQP;
+- Dykstra projection feasibility and idempotence;
+- long-only, caps, affine equalities, sector inequalities, hard turnover;
+- long-short dollar/beta neutrality, box bounds, and gross exposure.
+
+## Notebook execution
+
+Command:
+
+```bash
+python scripts/execute_notebooks.py
+```
+
+Result: every code cell and every embedded numerical assertion passed in all three notebooks.
+
+## Packaging smoke test
+
+The project was built as a PEP 517 wheel with local build dependencies, installed into a clean
+temporary target, imported from that installed target, and used to run `examples/quick_start.py`.
+The example converged with constraint violation below `1e-9`.
diff --git a/baseline/experiments/portfolio_pgd/examples/quick_start.py b/baseline/experiments/portfolio_pgd/examples/quick_start.py
new file mode 100644
index 00000000..cf34bad6
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/examples/quick_start.py
@@ -0,0 +1,51 @@
+"""Run a constrained portfolio optimization and print audit diagnostics."""
+
+from __future__ import annotations
+
+import numpy as np
+
+from portfolio_pgd import (
+    ConstraintSet,
+    PGDOptions,
+    PortfolioProblem,
+    PowerLawCost,
+    capped_long_only_portfolio,
+    factor_covariance,
+    solve_pgd,
+)
+
+
+def main() -> None:
+    n = 30
+    covariance, _ = factor_covariance(n, 4, seed=700)
+    rng = np.random.default_rng(701)
+    previous = capped_long_only_portfolio(n, cap=0.06, seed=702)
+    problem = PortfolioProblem(
+        alpha=rng.normal(scale=0.025, size=n),
+        covariance=covariance,
+        previous_holdings=previous,
+        risk_aversion=2.0,
+        quadratic_cost_matrix=0.3 + rng.random(n),
+        quadratic_cost_aversion=0.4,
+        nonlinear_cost=PowerLawCost(eta=0.01, p=1.5, epsilon=1.0e-3),
+    )
+    constraints = ConstraintSet(
+        n,
+        equality_matrix=np.ones((1, n)),
+        equality_target=np.array([1.0]),
+        lower_bounds=0.0,
+        upper_bounds=0.075,
+        turnover_limit=0.20,
+        turnover_center=previous,
+    )
+    result = solve_pgd(problem, constraints, options=PGDOptions(tolerance=1.0e-8))
+    print(f"status:                  {result.status}")
+    print(f"iterations:              {result.iterations}")
+    print(f"utility:                 {result.utility:.10f}")
+    print(f"two-way turnover:        {np.sum(np.abs(result.trades)):.10f}")
+    print(f"projected-gradient norm: {result.projected_gradient_norm:.3e}")
+    print(f"constraint violation:    {result.max_constraint_violation:.3e}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/baseline/experiments/portfolio_pgd/notebooks/01_quadratic_pgd_vs_standard_solvers.ipynb b/baseline/experiments/portfolio_pgd/notebooks/01_quadratic_pgd_vs_standard_solvers.ipynb
new file mode 100644
index 00000000..aa6ce697
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/notebooks/01_quadratic_pgd_vs_standard_solvers.ipynb
@@ -0,0 +1,275 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Quadratic portfolio construction: PGD versus exact and standard solvers\n",
+    "\n",
+    "This notebook implements the quadratic model\n",
+    "\n",
+    "$$\n",
+    "\\min_h\\;\\frac{\\lambda}{2}h^\\top Vh-\\alpha^\\top h\n",
+    "+\\frac{\\theta}{2}(h-h_-)^\\top Q(h-h_-)\n",
+    "\\quad\\text{subject to}\\quad C^\\top h=c.\n",
+    "$$\n",
+    "\n",
+    "We compare three independent solution routes:\n",
+    "\n",
+    "1. projected gradient descent with exact affine projection;\n",
+    "2. the exact equality-constrained KKT system;\n",
+    "3. SciPy's standard SLSQP solver.\n",
+    "\n",
+    "The numerical assertions at the end turn the notebook into an executable validation document."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## PGD iteration\n",
+    "\n",
+    "With negative utility denoted by $F(h)$,\n",
+    "\n",
+    "$$\n",
+    "\\nabla F(h)=\\lambda Vh+\\theta Q(h-h_-)-\\alpha.\n",
+    "$$\n",
+    "\n",
+    "The Euclidean affine projection is\n",
+    "\n",
+    "$$\n",
+    "\\Pi(z)=z-C(C^\\top C)^{\\dagger}(C^\\top z-c).\n",
+    "$$\n",
+    "\n",
+    "The implementation uses a projected-majorization line search, records the projected-gradient norm,\n",
+    "and refuses to return a portfolio whose constraints exceed the requested tolerance."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pathlib import Path\n",
+    "import sys\n",
+    "\n",
+    "ROOT = Path.cwd().resolve()\n",
+    "if ROOT.name == \"notebooks\":\n",
+    "    ROOT = ROOT.parent\n",
+    "sys.path.insert(0, str(ROOT / \"src\"))\n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "\n",
+    "plt.style.use(\"seaborn-v0_8-whitegrid\")\n",
+    "pd.set_option(\"display.float_format\", lambda value: f\"{value:,.8g}\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from portfolio_pgd import (\n",
+    "    ConstraintSet,\n",
+    "    PGDOptions,\n",
+    "    PortfolioProblem,\n",
+    "    factor_covariance,\n",
+    "    solve_pgd,\n",
+    "    solve_quadratic_kkt,\n",
+    "    solve_scipy_slsqp,\n",
+    ")\n",
+    "\n",
+    "n_assets = 36\n",
+    "rng = np.random.default_rng(1201)\n",
+    "covariance, loadings = factor_covariance(n_assets, 4, seed=1202, specific_risk=0.12)\n",
+    "alpha = rng.normal(scale=0.025, size=n_assets)\n",
+    "previous = rng.normal(scale=0.01, size=n_assets)\n",
+    "q_diagonal = 0.4 + rng.random(n_assets)\n",
+    "\n",
+    "problem = PortfolioProblem(\n",
+    "    alpha=alpha,\n",
+    "    covariance=covariance,\n",
+    "    previous_holdings=previous,\n",
+    "    risk_aversion=2.25,\n",
+    "    quadratic_cost_matrix=q_diagonal,\n",
+    "    quadratic_cost_aversion=0.75,\n",
+    ")\n",
+    "\n",
+    "# Full investment and one factor-exposure target.\n",
+    "factor_direction = loadings[:, 0] - np.mean(loadings[:, 0])\n",
+    "C_transpose = np.vstack([np.ones(n_assets), factor_direction])\n",
+    "targets = np.array([1.0, 0.0])\n",
+    "constraints = ConstraintSet(\n",
+    "    n_assets,\n",
+    "    equality_matrix=C_transpose,\n",
+    "    equality_target=targets,\n",
+    ")\n",
+    "print(\"Minimum covariance eigenvalue:\", np.linalg.eigvalsh(covariance).min())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pgd = solve_pgd(\n",
+    "    problem,\n",
+    "    constraints,\n",
+    "    options=PGDOptions(max_iterations=25_000, tolerance=2e-9),\n",
+    ")\n",
+    "kkt = solve_quadratic_kkt(problem, constraints)\n",
+    "slsqp = solve_scipy_slsqp(problem, constraints)\n",
+    "\n",
+    "comparison = pd.DataFrame(\n",
+    "    {\n",
+    "        \"objective\": [pgd.objective, kkt.objective, slsqp.objective],\n",
+    "        \"utility\": [pgd.utility, -kkt.objective, -slsqp.objective],\n",
+    "        \"distance_to_KKT\": [\n",
+    "            np.linalg.norm(pgd.holdings - kkt.holdings),\n",
+    "            0.0,\n",
+    "            np.linalg.norm(slsqp.holdings - kkt.holdings),\n",
+    "        ],\n",
+    "        \"max_constraint_violation\": [\n",
+    "            constraints.max_violation(pgd.holdings),\n",
+    "            constraints.max_violation(kkt.holdings),\n",
+    "            constraints.max_violation(slsqp.holdings),\n",
+    "        ],\n",
+    "    },\n",
+    "    index=[\"PGD\", \"Exact KKT\", \"SciPy SLSQP\"],\n",
+    ")\n",
+    "print(comparison.to_string())\n",
+    "print(f\"\\nPGD status={pgd.status}; iterations={pgd.iterations}; SLSQP={slsqp.message}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Convergence audit\n",
+    "\n",
+    "For a convex quadratic, the KKT objective is the global minimum of the estimated problem. The first\n",
+    "panel plots the PGD objective gap; the second shows the norm of the projected-gradient mapping."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "history = pd.DataFrame(pgd.history)\n",
+    "objective_gap = np.maximum(history[\"objective\"] - kkt.objective, 1e-18)\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n",
+    "axes[0].semilogy(history[\"iteration\"], objective_gap)\n",
+    "axes[0].set(title=\"Objective gap to exact KKT\", xlabel=\"Iteration\", ylabel=\"F(h) - F(h*)\")\n",
+    "valid = history[\"projected_gradient_norm\"].notna()\n",
+    "axes[1].semilogy(\n",
+    "    history.loc[valid, \"iteration\"],\n",
+    "    np.maximum(history.loc[valid, \"projected_gradient_norm\"], 1e-18),\n",
+    ")\n",
+    "axes[1].set(title=\"First-order residual\", xlabel=\"Iteration\", ylabel=\"Projected-gradient norm\")\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Alpha, persistence, and constraint decomposition\n",
+    "\n",
+    "Let $H=\\lambda V+\\theta Q$ and $b=\\alpha+\\theta Qh_-$. The unconstrained target is $H^{-1}b$.\n",
+    "The exact constrained solution decomposes as\n",
+    "\n",
+    "$$\n",
+    "h^*=\\underbrace{H^{-1}\\alpha}_{\\text{alpha}}\n",
+    "+\\underbrace{H^{-1}\\theta Qh_-}_{\\text{persistence}}\n",
+    "-\\underbrace{H^{-1}C\\mu}_{\\text{constraint correction}}.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "H = problem.quadratic_hessian\n",
+    "Q = np.asarray(problem.quadratic_cost_matrix)\n",
+    "A = np.asarray(constraints.equality_matrix)\n",
+    "c = np.asarray(constraints.equality_target)\n",
+    "\n",
+    "alpha_component = np.linalg.solve(H, alpha)\n",
+    "persistence_component = np.linalg.solve(\n",
+    "    H, problem.quadratic_cost_aversion * Q @ previous\n",
+    ")\n",
+    "unconstrained = alpha_component + persistence_component\n",
+    "H_inv_A_T = np.linalg.solve(H, A.T)\n",
+    "mu = np.linalg.solve(A @ H_inv_A_T, A @ unconstrained - c)\n",
+    "constraint_component = -H_inv_A_T @ mu\n",
+    "reconstructed = alpha_component + persistence_component + constraint_component\n",
+    "\n",
+    "decomposition = pd.DataFrame(\n",
+    "    {\n",
+    "        \"L2 norm\": [\n",
+    "            np.linalg.norm(alpha_component),\n",
+    "            np.linalg.norm(persistence_component),\n",
+    "            np.linalg.norm(constraint_component),\n",
+    "            np.linalg.norm(reconstructed),\n",
+    "        ],\n",
+    "        \"net exposure\": [\n",
+    "            np.sum(alpha_component),\n",
+    "            np.sum(persistence_component),\n",
+    "            np.sum(constraint_component),\n",
+    "            np.sum(reconstructed),\n",
+    "        ],\n",
+    "    },\n",
+    "    index=[\"alpha\", \"persistence\", \"constraint correction\", \"total\"],\n",
+    ")\n",
+    "print(decomposition.to_string())\n",
+    "print(\"Reconstruction error:\", np.linalg.norm(reconstructed - kkt.holdings))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Executable acceptance tests"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "assert pgd.converged\n",
+    "assert slsqp.success\n",
+    "assert constraints.max_violation(pgd.holdings) < 1e-8\n",
+    "assert np.linalg.norm(pgd.holdings - kkt.holdings) < 5e-7\n",
+    "assert np.linalg.norm(slsqp.holdings - kkt.holdings) < 5e-6\n",
+    "assert np.linalg.norm(reconstructed - kkt.holdings) < 1e-10\n",
+    "print(\"All quadratic validation checks passed.\")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/baseline/experiments/portfolio_pgd/notebooks/02_nonlinear_transaction_costs.ipynb b/baseline/experiments/portfolio_pgd/notebooks/02_nonlinear_transaction_costs.ipynb
new file mode 100644
index 00000000..35279181
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/notebooks/02_nonlinear_transaction_costs.ipynb
@@ -0,0 +1,268 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Nonlinear convex transaction costs\n",
+    "\n",
+    "This notebook replaces the quadratic-only trading model with a separable power-law impact cost\n",
+    "\n",
+    "$$\n",
+    "c(t)=\\sum_i\\eta_i\\left[(t_i^2+\\epsilon^2)^{p/2}-\\epsilon^p\\right],\n",
+    "\\qquad t=h-h_-.\n",
+    "$$\n",
+    "\n",
+    "For $p=3/2$ this is a smooth approximation to square-root-impact total cost. The smoothing parameter\n",
+    "$\\epsilon>0$ avoids the divergent curvature of $|t|^{3/2}$ at zero without changing convexity."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pathlib import Path\n",
+    "import sys\n",
+    "\n",
+    "ROOT = Path.cwd().resolve()\n",
+    "if ROOT.name == \"notebooks\":\n",
+    "    ROOT = ROOT.parent\n",
+    "sys.path.insert(0, str(ROOT / \"src\"))\n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "\n",
+    "plt.style.use(\"seaborn-v0_8-whitegrid\")\n",
+    "pd.set_option(\"display.float_format\", lambda value: f\"{value:,.8g}\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from portfolio_pgd import (\n",
+    "    ConstraintSet,\n",
+    "    PGDOptions,\n",
+    "    PortfolioProblem,\n",
+    "    PowerLawCost,\n",
+    "    capped_long_only_portfolio,\n",
+    "    factor_covariance,\n",
+    "    solve_pgd,\n",
+    "    solve_scipy_slsqp,\n",
+    ")\n",
+    "\n",
+    "n_assets = 30\n",
+    "rng = np.random.default_rng(2201)\n",
+    "covariance, _ = factor_covariance(n_assets, 4, seed=2202, specific_risk=0.15)\n",
+    "previous = capped_long_only_portfolio(n_assets, cap=0.055, seed=2203)\n",
+    "eta = 0.008 + 0.012 * rng.random(n_assets)\n",
+    "\n",
+    "cost = PowerLawCost(eta=eta, p=1.5, epsilon=1e-3)\n",
+    "problem = PortfolioProblem(\n",
+    "    alpha=rng.normal(scale=0.03, size=n_assets),\n",
+    "    covariance=covariance,\n",
+    "    previous_holdings=previous,\n",
+    "    risk_aversion=1.8,\n",
+    "    quadratic_cost_matrix=0.2 + rng.random(n_assets),\n",
+    "    quadratic_cost_aversion=0.25,\n",
+    "    nonlinear_cost=cost,\n",
+    ")\n",
+    "constraints = ConstraintSet(\n",
+    "    n_assets,\n",
+    "    equality_matrix=np.ones((1, n_assets)),\n",
+    "    equality_target=np.array([1.0]),\n",
+    "    lower_bounds=0.0,\n",
+    "    upper_bounds=0.075,\n",
+    ")\n",
+    "\n",
+    "pgd = solve_pgd(\n",
+    "    problem,\n",
+    "    constraints,\n",
+    "    options=PGDOptions(max_iterations=25_000, tolerance=5e-8),\n",
+    ")\n",
+    "slsqp = solve_scipy_slsqp(problem, constraints)\n",
+    "print(f\"PGD: {pgd.status} in {pgd.iterations} iterations\")\n",
+    "print(f\"SLSQP: success={slsqp.success}; {slsqp.message}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Independent solver comparison\n",
+    "\n",
+    "The nonlinear objective is convex but no longer quadratic. Therefore the notebook compares PGD to\n",
+    "SciPy SLSQP rather than to a linear KKT solve."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "comparison = pd.DataFrame(\n",
+    "    {\n",
+    "        \"objective\": [pgd.objective, slsqp.objective],\n",
+    "        \"utility\": [pgd.utility, -slsqp.objective],\n",
+    "        \"turnover\": [np.sum(np.abs(pgd.trades)), np.sum(np.abs(slsqp.holdings - previous))],\n",
+    "        \"distance_to_SLSQP\": [np.linalg.norm(pgd.holdings - slsqp.holdings), 0.0],\n",
+    "        \"constraint_violation\": [\n",
+    "            constraints.max_violation(pgd.holdings),\n",
+    "            constraints.max_violation(slsqp.holdings),\n",
+    "        ],\n",
+    "    },\n",
+    "    index=[\"PGD\", \"SciPy SLSQP\"],\n",
+    ")\n",
+    "print(comparison.to_string())"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Cost, marginal cost, and curvature\n",
+    "\n",
+    "The gradient enters PGD directly. The Hessian diagonal is used only to initialize a conservative\n",
+    "step; the majorization line search supplies the actual global safeguard."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "trade_grid = np.linspace(-0.08, 0.08, 401)\n",
+    "unit_cost = PowerLawCost(eta=1.0, p=1.5, epsilon=1e-3)\n",
+    "cost_values = np.array([unit_cost.value(np.array([trade])) for trade in trade_grid])\n",
+    "marginal = np.array([unit_cost.gradient(np.array([trade]))[0] for trade in trade_grid])\n",
+    "curvature = np.array([unit_cost.hessian_diag(np.array([trade]))[0] for trade in trade_grid])\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n",
+    "axes[0].plot(trade_grid, cost_values)\n",
+    "axes[0].set(title=\"Power-law cost\", xlabel=\"Trade\", ylabel=\"Cost / eta\")\n",
+    "axes[1].plot(trade_grid, marginal)\n",
+    "axes[1].set(title=\"Marginal cost\", xlabel=\"Trade\", ylabel=\"dc/dt / eta\")\n",
+    "axes[2].plot(trade_grid, curvature)\n",
+    "axes[2].set(title=\"Local curvature\", xlabel=\"Trade\", ylabel=\"d\u00b2c/dt\u00b2 / eta\")\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Convergence and the realized trade distribution"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "history = pd.DataFrame(pgd.history)\n",
+    "valid = history[\"projected_gradient_norm\"].notna()\n",
+    "fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n",
+    "axes[0].semilogy(\n",
+    "    history.loc[valid, \"iteration\"],\n",
+    "    np.maximum(history.loc[valid, \"projected_gradient_norm\"], 1e-18),\n",
+    ")\n",
+    "axes[0].set(title=\"Projected-gradient residual\", xlabel=\"Iteration\", ylabel=\"Norm\")\n",
+    "axes[1].bar(np.arange(n_assets), pgd.trades)\n",
+    "axes[1].axhline(0.0, color=\"black\", linewidth=0.8)\n",
+    "axes[1].set(title=\"Optimal trades\", xlabel=\"Asset\", ylabel=\"Weight change\")\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Sensitivity to the impact exponent\n",
+    "\n",
+    "Holding every other input fixed, we resolve the portfolio for several convex exponents. This is an\n",
+    "algorithmic comparison\u2014not a claim that any one exponent is universally appropriate."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rows = []\n",
+    "for exponent in [1.25, 1.5, 2.0, 3.0]:\n",
+    "    exponent_problem = PortfolioProblem(\n",
+    "        alpha=problem.alpha,\n",
+    "        covariance=problem.covariance,\n",
+    "        previous_holdings=previous,\n",
+    "        risk_aversion=problem.risk_aversion,\n",
+    "        quadratic_cost_matrix=problem.quadratic_cost_matrix,\n",
+    "        quadratic_cost_aversion=problem.quadratic_cost_aversion,\n",
+    "        nonlinear_cost=PowerLawCost(eta=eta, p=exponent, epsilon=1e-3),\n",
+    "    )\n",
+    "    solved = solve_pgd(\n",
+    "        exponent_problem,\n",
+    "        constraints,\n",
+    "        options=PGDOptions(max_iterations=25_000, tolerance=2e-7),\n",
+    "    )\n",
+    "    rows.append(\n",
+    "        {\n",
+    "            \"p\": exponent,\n",
+    "            \"converged\": solved.converged,\n",
+    "            \"utility\": solved.utility,\n",
+    "            \"turnover\": np.sum(np.abs(solved.trades)),\n",
+    "            \"max_abs_trade\": np.max(np.abs(solved.trades)),\n",
+    "            \"iterations\": solved.iterations,\n",
+    "        }\n",
+    "    )\n",
+    "sensitivity = pd.DataFrame(rows).set_index(\"p\")\n",
+    "print(sensitivity.to_string())"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Executable acceptance tests"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "assert pgd.converged\n",
+    "assert slsqp.success\n",
+    "assert constraints.max_violation(pgd.holdings) < 2e-8\n",
+    "assert abs(pgd.objective - slsqp.objective) < 2e-7\n",
+    "assert np.linalg.norm(pgd.holdings - slsqp.holdings) < 8e-4\n",
+    "assert sensitivity[\"converged\"].all()\n",
+    "print(\"All nonlinear-cost validation checks passed.\")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/baseline/experiments/portfolio_pgd/notebooks/03_realistic_constraints.ipynb b/baseline/experiments/portfolio_pgd/notebooks/03_realistic_constraints.ipynb
new file mode 100644
index 00000000..3452199d
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/notebooks/03_realistic_constraints.ipynb
@@ -0,0 +1,317 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Realistic institutional constraints\n",
+    "\n",
+    "This notebook combines the smooth nonlinear objective with a realistic long-only mandate:\n",
+    "\n",
+    "- fully invested;\n",
+    "- one fixed factor exposure;\n",
+    "- long-only and per-name caps;\n",
+    "- sector lower and upper bounds;\n",
+    "- hard two-way turnover.\n",
+    "\n",
+    "The PGD projection is the Euclidean projection onto the **joint intersection**. Dykstra's algorithm\n",
+    "cycles over analytic projectors while retaining correction terms, so this is not naive sequential\n",
+    "clipping."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pathlib import Path\n",
+    "import sys\n",
+    "\n",
+    "ROOT = Path.cwd().resolve()\n",
+    "if ROOT.name == \"notebooks\":\n",
+    "    ROOT = ROOT.parent\n",
+    "sys.path.insert(0, str(ROOT / \"src\"))\n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "\n",
+    "plt.style.use(\"seaborn-v0_8-whitegrid\")\n",
+    "pd.set_option(\"display.float_format\", lambda value: f\"{value:,.8g}\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from portfolio_pgd import (\n",
+    "    ConstraintSet,\n",
+    "    PGDOptions,\n",
+    "    PortfolioProblem,\n",
+    "    PowerLawCost,\n",
+    "    capped_long_only_portfolio,\n",
+    "    factor_covariance,\n",
+    "    sector_membership,\n",
+    "    solve_pgd,\n",
+    "    solve_scipy_slsqp,\n",
+    ")\n",
+    "\n",
+    "n_assets = 40\n",
+    "n_sectors = 5\n",
+    "rng = np.random.default_rng(3301)\n",
+    "covariance, loadings = factor_covariance(n_assets, 5, seed=3302, specific_risk=0.18)\n",
+    "previous = capped_long_only_portfolio(n_assets, cap=0.045, seed=3303)\n",
+    "sectors = sector_membership(n_assets, n_sectors)\n",
+    "previous_sector = sectors @ previous\n",
+    "\n",
+    "# Sector bands are centered on the existing portfolio, making feasibility explicit.\n",
+    "sector_lower = np.maximum(0.12, previous_sector - 0.035)\n",
+    "sector_upper = np.minimum(0.30, previous_sector + 0.035)\n",
+    "A_ub = np.vstack([sectors, -sectors])\n",
+    "b_ub = np.concatenate([sector_upper, -sector_lower])\n",
+    "\n",
+    "factor_direction = loadings[:, 0] - np.mean(loadings[:, 0])\n",
+    "A_eq = np.vstack([np.ones(n_assets), factor_direction])\n",
+    "b_eq = np.array([1.0, float(factor_direction @ previous)])\n",
+    "\n",
+    "constraints = ConstraintSet(\n",
+    "    n_assets,\n",
+    "    equality_matrix=A_eq,\n",
+    "    equality_target=b_eq,\n",
+    "    inequality_matrix=A_ub,\n",
+    "    inequality_upper=b_ub,\n",
+    "    lower_bounds=0.0,\n",
+    "    upper_bounds=0.06,\n",
+    "    turnover_limit=0.22,\n",
+    "    turnover_center=previous,\n",
+    ")\n",
+    "\n",
+    "problem = PortfolioProblem(\n",
+    "    alpha=rng.normal(scale=0.035, size=n_assets),\n",
+    "    covariance=covariance,\n",
+    "    previous_holdings=previous,\n",
+    "    risk_aversion=1.6,\n",
+    "    quadratic_cost_matrix=0.2 + rng.random(n_assets),\n",
+    "    quadratic_cost_aversion=0.25,\n",
+    "    nonlinear_cost=PowerLawCost(eta=0.006 + 0.006 * rng.random(n_assets), p=1.5, epsilon=1e-3),\n",
+    ")\n",
+    "print(\"Starting portfolio violations:\", constraints.violations(previous))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pgd = solve_pgd(\n",
+    "    problem,\n",
+    "    constraints,\n",
+    "    options=PGDOptions(\n",
+    "        max_iterations=30_000,\n",
+    "        tolerance=1e-7,\n",
+    "        projection_tolerance=2e-10,\n",
+    "    ),\n",
+    ")\n",
+    "slsqp = solve_scipy_slsqp(problem, constraints, tolerance=1e-10)\n",
+    "\n",
+    "comparison = pd.DataFrame(\n",
+    "    {\n",
+    "        \"objective\": [pgd.objective, slsqp.objective],\n",
+    "        \"utility\": [pgd.utility, -slsqp.objective],\n",
+    "        \"turnover\": [np.sum(np.abs(pgd.trades)), np.sum(np.abs(slsqp.holdings - previous))],\n",
+    "        \"distance_to_SLSQP\": [np.linalg.norm(pgd.holdings - slsqp.holdings), 0.0],\n",
+    "        \"constraint_violation\": [\n",
+    "            constraints.max_violation(pgd.holdings),\n",
+    "            constraints.max_violation(slsqp.holdings),\n",
+    "        ],\n",
+    "    },\n",
+    "    index=[\"PGD\", \"SciPy SLSQP\"],\n",
+    ")\n",
+    "print(comparison.to_string())\n",
+    "print(f\"\\nPGD status={pgd.status}; iterations={pgd.iterations}; SLSQP={slsqp.message}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Constraint audit"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "audit = pd.DataFrame(\n",
+    "    {\n",
+    "        \"PGD violation\": constraints.violations(pgd.holdings),\n",
+    "        \"SLSQP violation\": constraints.violations(slsqp.holdings),\n",
+    "    }\n",
+    ")\n",
+    "print(audit.to_string())\n",
+    "\n",
+    "sector_exposure = sectors @ pgd.holdings\n",
+    "sector_table = pd.DataFrame(\n",
+    "    {\n",
+    "        \"lower\": sector_lower,\n",
+    "        \"previous\": previous_sector,\n",
+    "        \"optimized\": sector_exposure,\n",
+    "        \"upper\": sector_upper,\n",
+    "    },\n",
+    "    index=[f\"Sector {index}\" for index in range(n_sectors)],\n",
+    ")\n",
+    "print(\"\\nSector exposures:\\n\", sector_table.to_string())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axes = plt.subplots(1, 2, figsize=(13, 4))\n",
+    "asset_index = np.arange(n_assets)\n",
+    "axes[0].plot(asset_index, previous, \"o-\", label=\"Previous\", markersize=3)\n",
+    "axes[0].plot(asset_index, pgd.holdings, \"o-\", label=\"Optimized\", markersize=3)\n",
+    "axes[0].axhline(0.06, color=\"black\", linestyle=\"--\", linewidth=1, label=\"Per-name cap\")\n",
+    "axes[0].set(title=\"Holdings\", xlabel=\"Asset\", ylabel=\"Weight\")\n",
+    "axes[0].legend()\n",
+    "axes[1].bar(asset_index, pgd.trades)\n",
+    "axes[1].axhline(0.0, color=\"black\", linewidth=0.8)\n",
+    "axes[1].set(title=f\"Trades (L1={np.sum(np.abs(pgd.trades)):.4f})\", xlabel=\"Asset\", ylabel=\"Weight change\")\n",
+    "plt.tight_layout()\n",
+    "plt.show()\n",
+    "\n",
+    "sector_table.plot(kind=\"bar\", figsize=(11, 4))\n",
+    "plt.title(\"Sector exposure audit\")\n",
+    "plt.ylabel(\"Portfolio weight\")\n",
+    "plt.xticks(rotation=0)\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Optimization convergence"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "history = pd.DataFrame(pgd.history)\n",
+    "valid = history[\"projected_gradient_norm\"].notna()\n",
+    "fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n",
+    "axes[0].plot(history[\"iteration\"], history[\"objective\"])\n",
+    "axes[0].axhline(slsqp.objective, color=\"black\", linestyle=\"--\", label=\"SLSQP\")\n",
+    "axes[0].set(title=\"Objective\", xlabel=\"Iteration\", ylabel=\"Negative utility\")\n",
+    "axes[0].legend()\n",
+    "axes[1].semilogy(\n",
+    "    history.loc[valid, \"iteration\"],\n",
+    "    np.maximum(history.loc[valid, \"projected_gradient_norm\"], 1e-18),\n",
+    ")\n",
+    "axes[1].set(title=\"Projected-gradient residual\", xlabel=\"Iteration\", ylabel=\"Norm\")\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Long-short extension\n",
+    "\n",
+    "The same projection engine can construct a dollar-neutral, beta-neutral portfolio with per-name and\n",
+    "gross-exposure limits. This remains convex; gross exposure is an $L_1$ ball centered at zero."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "long_short_constraints = ConstraintSet(\n",
+    "    n_assets,\n",
+    "    equality_matrix=np.vstack([np.ones(n_assets), loadings[:, 1]]),\n",
+    "    equality_target=np.array([0.0, 0.0]),\n",
+    "    lower_bounds=-0.08,\n",
+    "    upper_bounds=0.08,\n",
+    "    gross_exposure_limit=1.0,\n",
+    ")\n",
+    "long_short_problem = PortfolioProblem(\n",
+    "    alpha=2.5 * problem.alpha,\n",
+    "    covariance=problem.covariance,\n",
+    "    previous_holdings=np.zeros(n_assets),\n",
+    "    risk_aversion=0.9,\n",
+    "    quadratic_cost_matrix=np.ones(n_assets),\n",
+    "    quadratic_cost_aversion=0.08,\n",
+    ")\n",
+    "long_short = solve_pgd(\n",
+    "    long_short_problem,\n",
+    "    long_short_constraints,\n",
+    "    options=PGDOptions(max_iterations=25_000, tolerance=1e-8),\n",
+    ")\n",
+    "print(\n",
+    "    pd.Series(\n",
+    "        {\n",
+    "            \"status\": long_short.status,\n",
+    "            \"net exposure\": np.sum(long_short.holdings),\n",
+    "            \"gross exposure\": np.sum(np.abs(long_short.holdings)),\n",
+    "            \"beta exposure\": loadings[:, 1] @ long_short.holdings,\n",
+    "            \"maximum position\": np.max(np.abs(long_short.holdings)),\n",
+    "            \"constraint violation\": long_short.max_constraint_violation,\n",
+    "        }\n",
+    "    ).to_string()\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Executable acceptance tests"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "assert pgd.converged\n",
+    "assert slsqp.success\n",
+    "assert constraints.max_violation(pgd.holdings) < 2e-7\n",
+    "assert constraints.max_violation(slsqp.holdings) < 2e-6\n",
+    "assert abs(pgd.objective - slsqp.objective) < 2e-5\n",
+    "assert np.sum(np.abs(pgd.trades)) <= 0.22 + 2e-7\n",
+    "assert long_short.converged\n",
+    "assert long_short_constraints.max_violation(long_short.holdings) < 2e-7\n",
+    "print(\"All realistic-constraint validation checks passed.\")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/baseline/experiments/portfolio_pgd/pyproject.toml b/baseline/experiments/portfolio_pgd/pyproject.toml
new file mode 100644
index 00000000..7af23ef8
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/pyproject.toml
@@ -0,0 +1,39 @@
+[build-system]
+requires = ["setuptools>=68", "wheel"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "portfolio-pgd"
+version = "1.0.0"
+description = "Projected-gradient solvers for constrained portfolio construction"
+readme = "README.md"
+requires-python = ">=3.10"
+license = {text = "MIT"}
+authors = [{name = "Calculation Consulting"}]
+dependencies = [
+  "numpy>=1.24",
+  "scipy>=1.10",
+]
+
+[project.optional-dependencies]
+notebook = [
+  "jupyterlab>=4.0",
+  "notebook>=7.0",
+  "matplotlib>=3.7",
+  "pandas>=2.0",
+]
+dev = [
+  "pytest>=7.4",
+  "ruff>=0.6",
+]
+
+[tool.setuptools.packages.find]
+where = ["src"]
+
+[tool.pytest.ini_options]
+testpaths = ["tests"]
+addopts = "-ra"
+
+[tool.ruff]
+line-length = 100
+target-version = "py310"
diff --git a/baseline/experiments/portfolio_pgd/requirements.txt b/baseline/experiments/portfolio_pgd/requirements.txt
new file mode 100644
index 00000000..06c0b330
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/requirements.txt
@@ -0,0 +1,6 @@
+numpy>=1.24
+scipy>=1.10
+matplotlib>=3.7
+pandas>=2.0
+jupyterlab>=4.0
+notebook>=7.0
diff --git a/baseline/experiments/portfolio_pgd/scripts/build_notebooks.py b/baseline/experiments/portfolio_pgd/scripts/build_notebooks.py
new file mode 100644
index 00000000..d7db68ac
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/scripts/build_notebooks.py
@@ -0,0 +1,588 @@
+"""Generate the checked-in notebooks using only the Python standard library."""
+
+from __future__ import annotations
+
+import json
+from pathlib import Path
+
+ROOT = Path(__file__).resolve().parents[1]
+NOTEBOOK_DIR = ROOT / "notebooks"
+
+
+def markdown(source: str) -> dict:
+    return {"cell_type": "markdown", "metadata": {}, "source": source.splitlines(keepends=True)}
+
+
+def code(source: str) -> dict:
+    return {
+        "cell_type": "code",
+        "execution_count": None,
+        "metadata": {},
+        "outputs": [],
+        "source": source.splitlines(keepends=True),
+    }
+
+
+def write_notebook(filename: str, cells: list[dict]) -> None:
+    payload = {
+        "cells": cells,
+        "metadata": {
+            "kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"},
+            "language_info": {"name": "python", "version": "3.10"},
+        },
+        "nbformat": 4,
+        "nbformat_minor": 5,
+    }
+    NOTEBOOK_DIR.mkdir(parents=True, exist_ok=True)
+    (NOTEBOOK_DIR / filename).write_text(json.dumps(payload, indent=1) + "\n", encoding="utf-8")
+
+
+COMMON_SETUP = r'''from pathlib import Path
+import sys
+
+ROOT = Path.cwd().resolve()
+if ROOT.name == "notebooks":
+    ROOT = ROOT.parent
+sys.path.insert(0, str(ROOT / "src"))
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+
+plt.style.use("seaborn-v0_8-whitegrid")
+pd.set_option("display.float_format", lambda value: f"{value:,.8g}")'''
+
+
+def build_quadratic() -> None:
+    cells = [
+        markdown(r'''# Quadratic portfolio construction: PGD versus exact and standard solvers
+
+This notebook implements the quadratic model
+
+$$
+\min_h\;\frac{\lambda}{2}h^\top Vh-\alpha^\top h
++\frac{\theta}{2}(h-h_-)^\top Q(h-h_-)
+\quad\text{subject to}\quad C^\top h=c.
+$$
+
+We compare three independent solution routes:
+
+1. projected gradient descent with exact affine projection;
+2. the exact equality-constrained KKT system;
+3. SciPy's standard SLSQP solver.
+
+The numerical assertions at the end turn the notebook into an executable validation document.'''),
+        markdown(r'''## PGD iteration
+
+With negative utility denoted by $F(h)$,
+
+$$
+\nabla F(h)=\lambda Vh+\theta Q(h-h_-)-\alpha.
+$$
+
+The Euclidean affine projection is
+
+$$
+\Pi(z)=z-C(C^\top C)^{\dagger}(C^\top z-c).
+$$
+
+The implementation uses a projected-majorization line search, records the projected-gradient norm,
+and refuses to return a portfolio whose constraints exceed the requested tolerance.'''),
+        code(COMMON_SETUP),
+        code(r'''from portfolio_pgd import (
+    ConstraintSet,
+    PGDOptions,
+    PortfolioProblem,
+    factor_covariance,
+    solve_pgd,
+    solve_quadratic_kkt,
+    solve_scipy_slsqp,
+)
+
+n_assets = 36
+rng = np.random.default_rng(1201)
+covariance, loadings = factor_covariance(n_assets, 4, seed=1202, specific_risk=0.12)
+alpha = rng.normal(scale=0.025, size=n_assets)
+previous = rng.normal(scale=0.01, size=n_assets)
+q_diagonal = 0.4 + rng.random(n_assets)
+
+problem = PortfolioProblem(
+    alpha=alpha,
+    covariance=covariance,
+    previous_holdings=previous,
+    risk_aversion=2.25,
+    quadratic_cost_matrix=q_diagonal,
+    quadratic_cost_aversion=0.75,
+)
+
+# Full investment and one factor-exposure target.
+factor_direction = loadings[:, 0] - np.mean(loadings[:, 0])
+C_transpose = np.vstack([np.ones(n_assets), factor_direction])
+targets = np.array([1.0, 0.0])
+constraints = ConstraintSet(
+    n_assets,
+    equality_matrix=C_transpose,
+    equality_target=targets,
+)
+print("Minimum covariance eigenvalue:", np.linalg.eigvalsh(covariance).min())'''),
+        code(r'''pgd = solve_pgd(
+    problem,
+    constraints,
+    options=PGDOptions(max_iterations=25_000, tolerance=2e-9),
+)
+kkt = solve_quadratic_kkt(problem, constraints)
+slsqp = solve_scipy_slsqp(problem, constraints)
+
+comparison = pd.DataFrame(
+    {
+        "objective": [pgd.objective, kkt.objective, slsqp.objective],
+        "utility": [pgd.utility, -kkt.objective, -slsqp.objective],
+        "distance_to_KKT": [
+            np.linalg.norm(pgd.holdings - kkt.holdings),
+            0.0,
+            np.linalg.norm(slsqp.holdings - kkt.holdings),
+        ],
+        "max_constraint_violation": [
+            constraints.max_violation(pgd.holdings),
+            constraints.max_violation(kkt.holdings),
+            constraints.max_violation(slsqp.holdings),
+        ],
+    },
+    index=["PGD", "Exact KKT", "SciPy SLSQP"],
+)
+print(comparison.to_string())
+print(f"\nPGD status={pgd.status}; iterations={pgd.iterations}; SLSQP={slsqp.message}")'''),
+        markdown(r'''## Convergence audit
+
+For a convex quadratic, the KKT objective is the global minimum of the estimated problem. The first
+panel plots the PGD objective gap; the second shows the norm of the projected-gradient mapping.'''),
+        code(r'''history = pd.DataFrame(pgd.history)
+objective_gap = np.maximum(history["objective"] - kkt.objective, 1e-18)
+
+fig, axes = plt.subplots(1, 2, figsize=(12, 4))
+axes[0].semilogy(history["iteration"], objective_gap)
+axes[0].set(title="Objective gap to exact KKT", xlabel="Iteration", ylabel="F(h) - F(h*)")
+valid = history["projected_gradient_norm"].notna()
+axes[1].semilogy(
+    history.loc[valid, "iteration"],
+    np.maximum(history.loc[valid, "projected_gradient_norm"], 1e-18),
+)
+axes[1].set(title="First-order residual", xlabel="Iteration", ylabel="Projected-gradient norm")
+plt.tight_layout()
+plt.show()'''),
+        markdown(r'''## Alpha, persistence, and constraint decomposition
+
+Let $H=\lambda V+\theta Q$ and $b=\alpha+\theta Qh_-$. The unconstrained target is $H^{-1}b$.
+The exact constrained solution decomposes as
+
+$$
+h^*=\underbrace{H^{-1}\alpha}_{\text{alpha}}
++\underbrace{H^{-1}\theta Qh_-}_{\text{persistence}}
+-\underbrace{H^{-1}C\mu}_{\text{constraint correction}}.
+$$'''),
+        code(r'''H = problem.quadratic_hessian
+Q = np.asarray(problem.quadratic_cost_matrix)
+A = np.asarray(constraints.equality_matrix)
+c = np.asarray(constraints.equality_target)
+
+alpha_component = np.linalg.solve(H, alpha)
+persistence_component = np.linalg.solve(
+    H, problem.quadratic_cost_aversion * Q @ previous
+)
+unconstrained = alpha_component + persistence_component
+H_inv_A_T = np.linalg.solve(H, A.T)
+mu = np.linalg.solve(A @ H_inv_A_T, A @ unconstrained - c)
+constraint_component = -H_inv_A_T @ mu
+reconstructed = alpha_component + persistence_component + constraint_component
+
+decomposition = pd.DataFrame(
+    {
+        "L2 norm": [
+            np.linalg.norm(alpha_component),
+            np.linalg.norm(persistence_component),
+            np.linalg.norm(constraint_component),
+            np.linalg.norm(reconstructed),
+        ],
+        "net exposure": [
+            np.sum(alpha_component),
+            np.sum(persistence_component),
+            np.sum(constraint_component),
+            np.sum(reconstructed),
+        ],
+    },
+    index=["alpha", "persistence", "constraint correction", "total"],
+)
+print(decomposition.to_string())
+print("Reconstruction error:", np.linalg.norm(reconstructed - kkt.holdings))'''),
+        markdown("## Executable acceptance tests"),
+        code(r'''assert pgd.converged
+assert slsqp.success
+assert constraints.max_violation(pgd.holdings) < 1e-8
+assert np.linalg.norm(pgd.holdings - kkt.holdings) < 5e-7
+assert np.linalg.norm(slsqp.holdings - kkt.holdings) < 5e-6
+assert np.linalg.norm(reconstructed - kkt.holdings) < 1e-10
+print("All quadratic validation checks passed.")'''),
+    ]
+    write_notebook("01_quadratic_pgd_vs_standard_solvers.ipynb", cells)
+
+
+def build_nonlinear() -> None:
+    cells = [
+        markdown(r'''# Nonlinear convex transaction costs
+
+This notebook replaces the quadratic-only trading model with a separable power-law impact cost
+
+$$
+c(t)=\sum_i\eta_i\left[(t_i^2+\epsilon^2)^{p/2}-\epsilon^p\right],
+\qquad t=h-h_-.
+$$
+
+For $p=3/2$ this is a smooth approximation to square-root-impact total cost. The smoothing parameter
+$\epsilon>0$ avoids the divergent curvature of $|t|^{3/2}$ at zero without changing convexity.'''),
+        code(COMMON_SETUP),
+        code(r'''from portfolio_pgd import (
+    ConstraintSet,
+    PGDOptions,
+    PortfolioProblem,
+    PowerLawCost,
+    capped_long_only_portfolio,
+    factor_covariance,
+    solve_pgd,
+    solve_scipy_slsqp,
+)
+
+n_assets = 30
+rng = np.random.default_rng(2201)
+covariance, _ = factor_covariance(n_assets, 4, seed=2202, specific_risk=0.15)
+previous = capped_long_only_portfolio(n_assets, cap=0.055, seed=2203)
+eta = 0.008 + 0.012 * rng.random(n_assets)
+
+cost = PowerLawCost(eta=eta, p=1.5, epsilon=1e-3)
+problem = PortfolioProblem(
+    alpha=rng.normal(scale=0.03, size=n_assets),
+    covariance=covariance,
+    previous_holdings=previous,
+    risk_aversion=1.8,
+    quadratic_cost_matrix=0.2 + rng.random(n_assets),
+    quadratic_cost_aversion=0.25,
+    nonlinear_cost=cost,
+)
+constraints = ConstraintSet(
+    n_assets,
+    equality_matrix=np.ones((1, n_assets)),
+    equality_target=np.array([1.0]),
+    lower_bounds=0.0,
+    upper_bounds=0.075,
+)
+
+pgd = solve_pgd(
+    problem,
+    constraints,
+    options=PGDOptions(max_iterations=25_000, tolerance=5e-8),
+)
+slsqp = solve_scipy_slsqp(problem, constraints)
+print(f"PGD: {pgd.status} in {pgd.iterations} iterations")
+print(f"SLSQP: success={slsqp.success}; {slsqp.message}")'''),
+        markdown(r'''## Independent solver comparison
+
+The nonlinear objective is convex but no longer quadratic. Therefore the notebook compares PGD to
+SciPy SLSQP rather than to a linear KKT solve.'''),
+        code(r'''comparison = pd.DataFrame(
+    {
+        "objective": [pgd.objective, slsqp.objective],
+        "utility": [pgd.utility, -slsqp.objective],
+        "turnover": [np.sum(np.abs(pgd.trades)), np.sum(np.abs(slsqp.holdings - previous))],
+        "distance_to_SLSQP": [np.linalg.norm(pgd.holdings - slsqp.holdings), 0.0],
+        "constraint_violation": [
+            constraints.max_violation(pgd.holdings),
+            constraints.max_violation(slsqp.holdings),
+        ],
+    },
+    index=["PGD", "SciPy SLSQP"],
+)
+print(comparison.to_string())'''),
+        markdown(r'''## Cost, marginal cost, and curvature
+
+The gradient enters PGD directly. The Hessian diagonal is used only to initialize a conservative
+step; the majorization line search supplies the actual global safeguard.'''),
+        code(r'''trade_grid = np.linspace(-0.08, 0.08, 401)
+unit_cost = PowerLawCost(eta=1.0, p=1.5, epsilon=1e-3)
+cost_values = np.array([unit_cost.value(np.array([trade])) for trade in trade_grid])
+marginal = np.array([unit_cost.gradient(np.array([trade]))[0] for trade in trade_grid])
+curvature = np.array([unit_cost.hessian_diag(np.array([trade]))[0] for trade in trade_grid])
+
+fig, axes = plt.subplots(1, 3, figsize=(15, 4))
+axes[0].plot(trade_grid, cost_values)
+axes[0].set(title="Power-law cost", xlabel="Trade", ylabel="Cost / eta")
+axes[1].plot(trade_grid, marginal)
+axes[1].set(title="Marginal cost", xlabel="Trade", ylabel="dc/dt / eta")
+axes[2].plot(trade_grid, curvature)
+axes[2].set(title="Local curvature", xlabel="Trade", ylabel="d²c/dt² / eta")
+plt.tight_layout()
+plt.show()'''),
+        markdown("## Convergence and the realized trade distribution"),
+        code(r'''history = pd.DataFrame(pgd.history)
+valid = history["projected_gradient_norm"].notna()
+fig, axes = plt.subplots(1, 2, figsize=(12, 4))
+axes[0].semilogy(
+    history.loc[valid, "iteration"],
+    np.maximum(history.loc[valid, "projected_gradient_norm"], 1e-18),
+)
+axes[0].set(title="Projected-gradient residual", xlabel="Iteration", ylabel="Norm")
+axes[1].bar(np.arange(n_assets), pgd.trades)
+axes[1].axhline(0.0, color="black", linewidth=0.8)
+axes[1].set(title="Optimal trades", xlabel="Asset", ylabel="Weight change")
+plt.tight_layout()
+plt.show()'''),
+        markdown(r'''## Sensitivity to the impact exponent
+
+Holding every other input fixed, we resolve the portfolio for several convex exponents. This is an
+algorithmic comparison—not a claim that any one exponent is universally appropriate.'''),
+        code(r'''rows = []
+for exponent in [1.25, 1.5, 2.0, 3.0]:
+    exponent_problem = PortfolioProblem(
+        alpha=problem.alpha,
+        covariance=problem.covariance,
+        previous_holdings=previous,
+        risk_aversion=problem.risk_aversion,
+        quadratic_cost_matrix=problem.quadratic_cost_matrix,
+        quadratic_cost_aversion=problem.quadratic_cost_aversion,
+        nonlinear_cost=PowerLawCost(eta=eta, p=exponent, epsilon=1e-3),
+    )
+    solved = solve_pgd(
+        exponent_problem,
+        constraints,
+        options=PGDOptions(max_iterations=25_000, tolerance=2e-7),
+    )
+    rows.append(
+        {
+            "p": exponent,
+            "converged": solved.converged,
+            "utility": solved.utility,
+            "turnover": np.sum(np.abs(solved.trades)),
+            "max_abs_trade": np.max(np.abs(solved.trades)),
+            "iterations": solved.iterations,
+        }
+    )
+sensitivity = pd.DataFrame(rows).set_index("p")
+print(sensitivity.to_string())'''),
+        markdown("## Executable acceptance tests"),
+        code(r'''assert pgd.converged
+assert slsqp.success
+assert constraints.max_violation(pgd.holdings) < 2e-8
+assert abs(pgd.objective - slsqp.objective) < 2e-7
+assert np.linalg.norm(pgd.holdings - slsqp.holdings) < 8e-4
+assert sensitivity["converged"].all()
+print("All nonlinear-cost validation checks passed.")'''),
+    ]
+    write_notebook("02_nonlinear_transaction_costs.ipynb", cells)
+
+
+def build_realistic() -> None:
+    cells = [
+        markdown(r'''# Realistic institutional constraints
+
+This notebook combines the smooth nonlinear objective with a realistic long-only mandate:
+
+- fully invested;
+- one fixed factor exposure;
+- long-only and per-name caps;
+- sector lower and upper bounds;
+- hard two-way turnover.
+
+The PGD projection is the Euclidean projection onto the **joint intersection**. Dykstra's algorithm
+cycles over analytic projectors while retaining correction terms, so this is not naive sequential
+clipping.'''),
+        code(COMMON_SETUP),
+        code(r'''from portfolio_pgd import (
+    ConstraintSet,
+    PGDOptions,
+    PortfolioProblem,
+    PowerLawCost,
+    capped_long_only_portfolio,
+    factor_covariance,
+    sector_membership,
+    solve_pgd,
+    solve_scipy_slsqp,
+)
+
+n_assets = 40
+n_sectors = 5
+rng = np.random.default_rng(3301)
+covariance, loadings = factor_covariance(n_assets, 5, seed=3302, specific_risk=0.18)
+previous = capped_long_only_portfolio(n_assets, cap=0.045, seed=3303)
+sectors = sector_membership(n_assets, n_sectors)
+previous_sector = sectors @ previous
+
+# Sector bands are centered on the existing portfolio, making feasibility explicit.
+sector_lower = np.maximum(0.12, previous_sector - 0.035)
+sector_upper = np.minimum(0.30, previous_sector + 0.035)
+A_ub = np.vstack([sectors, -sectors])
+b_ub = np.concatenate([sector_upper, -sector_lower])
+
+factor_direction = loadings[:, 0] - np.mean(loadings[:, 0])
+A_eq = np.vstack([np.ones(n_assets), factor_direction])
+b_eq = np.array([1.0, float(factor_direction @ previous)])
+
+constraints = ConstraintSet(
+    n_assets,
+    equality_matrix=A_eq,
+    equality_target=b_eq,
+    inequality_matrix=A_ub,
+    inequality_upper=b_ub,
+    lower_bounds=0.0,
+    upper_bounds=0.06,
+    turnover_limit=0.22,
+    turnover_center=previous,
+)
+
+problem = PortfolioProblem(
+    alpha=rng.normal(scale=0.035, size=n_assets),
+    covariance=covariance,
+    previous_holdings=previous,
+    risk_aversion=1.6,
+    quadratic_cost_matrix=0.2 + rng.random(n_assets),
+    quadratic_cost_aversion=0.25,
+    nonlinear_cost=PowerLawCost(eta=0.006 + 0.006 * rng.random(n_assets), p=1.5, epsilon=1e-3),
+)
+print("Starting portfolio violations:", constraints.violations(previous))'''),
+        code(r'''pgd = solve_pgd(
+    problem,
+    constraints,
+    options=PGDOptions(
+        max_iterations=30_000,
+        tolerance=1e-7,
+        projection_tolerance=2e-10,
+    ),
+)
+slsqp = solve_scipy_slsqp(problem, constraints, tolerance=1e-10)
+
+comparison = pd.DataFrame(
+    {
+        "objective": [pgd.objective, slsqp.objective],
+        "utility": [pgd.utility, -slsqp.objective],
+        "turnover": [np.sum(np.abs(pgd.trades)), np.sum(np.abs(slsqp.holdings - previous))],
+        "distance_to_SLSQP": [np.linalg.norm(pgd.holdings - slsqp.holdings), 0.0],
+        "constraint_violation": [
+            constraints.max_violation(pgd.holdings),
+            constraints.max_violation(slsqp.holdings),
+        ],
+    },
+    index=["PGD", "SciPy SLSQP"],
+)
+print(comparison.to_string())
+print(f"\nPGD status={pgd.status}; iterations={pgd.iterations}; SLSQP={slsqp.message}")'''),
+        markdown("## Constraint audit"),
+        code(r'''audit = pd.DataFrame(
+    {
+        "PGD violation": constraints.violations(pgd.holdings),
+        "SLSQP violation": constraints.violations(slsqp.holdings),
+    }
+)
+print(audit.to_string())
+
+sector_exposure = sectors @ pgd.holdings
+sector_table = pd.DataFrame(
+    {
+        "lower": sector_lower,
+        "previous": previous_sector,
+        "optimized": sector_exposure,
+        "upper": sector_upper,
+    },
+    index=[f"Sector {index}" for index in range(n_sectors)],
+)
+print("\nSector exposures:\n", sector_table.to_string())'''),
+        code(r'''fig, axes = plt.subplots(1, 2, figsize=(13, 4))
+asset_index = np.arange(n_assets)
+axes[0].plot(asset_index, previous, "o-", label="Previous", markersize=3)
+axes[0].plot(asset_index, pgd.holdings, "o-", label="Optimized", markersize=3)
+axes[0].axhline(0.06, color="black", linestyle="--", linewidth=1, label="Per-name cap")
+axes[0].set(title="Holdings", xlabel="Asset", ylabel="Weight")
+axes[0].legend()
+axes[1].bar(asset_index, pgd.trades)
+axes[1].axhline(0.0, color="black", linewidth=0.8)
+axes[1].set(title=f"Trades (L1={np.sum(np.abs(pgd.trades)):.4f})", xlabel="Asset", ylabel="Weight change")
+plt.tight_layout()
+plt.show()
+
+sector_table.plot(kind="bar", figsize=(11, 4))
+plt.title("Sector exposure audit")
+plt.ylabel("Portfolio weight")
+plt.xticks(rotation=0)
+plt.tight_layout()
+plt.show()'''),
+        markdown("## Optimization convergence"),
+        code(r'''history = pd.DataFrame(pgd.history)
+valid = history["projected_gradient_norm"].notna()
+fig, axes = plt.subplots(1, 2, figsize=(12, 4))
+axes[0].plot(history["iteration"], history["objective"])
+axes[0].axhline(slsqp.objective, color="black", linestyle="--", label="SLSQP")
+axes[0].set(title="Objective", xlabel="Iteration", ylabel="Negative utility")
+axes[0].legend()
+axes[1].semilogy(
+    history.loc[valid, "iteration"],
+    np.maximum(history.loc[valid, "projected_gradient_norm"], 1e-18),
+)
+axes[1].set(title="Projected-gradient residual", xlabel="Iteration", ylabel="Norm")
+plt.tight_layout()
+plt.show()'''),
+        markdown(r'''## Long-short extension
+
+The same projection engine can construct a dollar-neutral, beta-neutral portfolio with per-name and
+gross-exposure limits. This remains convex; gross exposure is an $L_1$ ball centered at zero.'''),
+        code(r'''long_short_constraints = ConstraintSet(
+    n_assets,
+    equality_matrix=np.vstack([np.ones(n_assets), loadings[:, 1]]),
+    equality_target=np.array([0.0, 0.0]),
+    lower_bounds=-0.08,
+    upper_bounds=0.08,
+    gross_exposure_limit=1.0,
+)
+long_short_problem = PortfolioProblem(
+    alpha=2.5 * problem.alpha,
+    covariance=problem.covariance,
+    previous_holdings=np.zeros(n_assets),
+    risk_aversion=0.9,
+    quadratic_cost_matrix=np.ones(n_assets),
+    quadratic_cost_aversion=0.08,
+)
+long_short = solve_pgd(
+    long_short_problem,
+    long_short_constraints,
+    options=PGDOptions(max_iterations=25_000, tolerance=1e-8),
+)
+print(
+    pd.Series(
+        {
+            "status": long_short.status,
+            "net exposure": np.sum(long_short.holdings),
+            "gross exposure": np.sum(np.abs(long_short.holdings)),
+            "beta exposure": loadings[:, 1] @ long_short.holdings,
+            "maximum position": np.max(np.abs(long_short.holdings)),
+            "constraint violation": long_short.max_constraint_violation,
+        }
+    ).to_string()
+)'''),
+        markdown("## Executable acceptance tests"),
+        code(r'''assert pgd.converged
+assert slsqp.success
+assert constraints.max_violation(pgd.holdings) < 2e-7
+assert constraints.max_violation(slsqp.holdings) < 2e-6
+assert abs(pgd.objective - slsqp.objective) < 2e-5
+assert np.sum(np.abs(pgd.trades)) <= 0.22 + 2e-7
+assert long_short.converged
+assert long_short_constraints.max_violation(long_short.holdings) < 2e-7
+print("All realistic-constraint validation checks passed.")'''),
+    ]
+    write_notebook("03_realistic_constraints.ipynb", cells)
+
+
+def main() -> None:
+    build_quadratic()
+    build_nonlinear()
+    build_realistic()
+    print(f"Wrote notebooks to {NOTEBOOK_DIR}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/baseline/experiments/portfolio_pgd/scripts/execute_notebooks.py b/baseline/experiments/portfolio_pgd/scripts/execute_notebooks.py
new file mode 100644
index 00000000..d473c8c9
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/scripts/execute_notebooks.py
@@ -0,0 +1,37 @@
+"""Execute every code cell without Jupyter to provide a lightweight CI smoke test."""
+
+from __future__ import annotations
+
+import json
+import os
+from pathlib import Path
+
+os.environ.setdefault("MPLBACKEND", "Agg")
+os.environ.setdefault("MPLCONFIGDIR", "/tmp/portfolio-pgd-matplotlib")
+
+ROOT = Path(__file__).resolve().parents[1]
+
+
+def execute(path: Path) -> None:
+    namespace: dict[str, object] = {"__name__": "__notebook__"}
+    payload = json.loads(path.read_text(encoding="utf-8"))
+    previous = Path.cwd()
+    os.chdir(ROOT)
+    try:
+        for index, cell in enumerate(payload["cells"]):
+            if cell["cell_type"] != "code":
+                continue
+            source = "".join(cell["source"])
+            exec(compile(source, f"{path.name}:cell-{index}", "exec"), namespace)
+    finally:
+        os.chdir(previous)
+    print(f"PASS {path.name}")
+
+
+def main() -> None:
+    for path in sorted((ROOT / "notebooks").glob("*.ipynb")):
+        execute(path)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/baseline/experiments/portfolio_pgd/scripts/run_portfolio_pgd_experiment.py b/baseline/experiments/portfolio_pgd/scripts/run_portfolio_pgd_experiment.py
new file mode 100755
index 00000000..b6e1d073
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/scripts/run_portfolio_pgd_experiment.py
@@ -0,0 +1,354 @@
+#!/usr/bin/env python3
+"""Run the complete portfolio-PGD validation campaign with visible progress logging."""
+
+from __future__ import annotations
+
+import argparse
+import csv
+import json
+import logging
+import platform
+import sys
+from dataclasses import asdict
+from datetime import datetime, timezone
+from pathlib import Path
+from time import perf_counter
+
+import numpy as np
+import scipy
+
+EXPERIMENT_DIR = Path(__file__).resolve().parents[1]
+SRC_DIR = EXPERIMENT_DIR / "src"
+if str(SRC_DIR) not in sys.path:
+    sys.path.insert(0, str(SRC_DIR))
+
+from portfolio_pgd import (  # noqa: E402
+    ConstraintSet,
+    PGDOptions,
+    PortfolioProblem,
+    PowerLawCost,
+    ProgressState,
+    capped_long_only_portfolio,
+    factor_covariance,
+    sector_membership,
+    solve_pgd,
+    solve_quadratic_kkt,
+    solve_scipy_slsqp,
+)
+
+LOGGER = logging.getLogger("portfolio_pgd.experiment")
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(
+        description="Run quadratic, nonlinear-cost, and realistic-constraint PGD benchmarks."
+    )
+    parser.add_argument("--output-dir", type=Path, required=True)
+    parser.add_argument("--seed", type=int, default=20260822)
+    parser.add_argument("--log-every", type=int, default=25)
+    parser.add_argument("--max-iterations", type=int, default=30_000)
+    parser.add_argument("--tolerance", type=float, default=1.0e-7)
+    return parser.parse_args()
+
+
+def configure_logging() -> None:
+    handler = logging.StreamHandler(sys.stdout)
+    handler.setFormatter(
+        logging.Formatter("%(asctime)s %(levelname)-8s %(message)s", datefmt="%Y-%m-%d %H:%M:%S")
+    )
+    LOGGER.handlers.clear()
+    LOGGER.addHandler(handler)
+    LOGGER.setLevel(logging.INFO)
+    LOGGER.propagate = False
+
+
+def stage(number: int, total: int, title: str) -> None:
+    LOGGER.info("=" * 72)
+    LOGGER.info("STAGE %d/%d %s", number, total, title)
+    LOGGER.info("=" * 72)
+
+
+def progress_logger(name: str):
+    def report(state: ProgressState) -> None:
+        residual = (
+            "initial"
+            if not np.isfinite(state.projected_gradient_norm)
+            else f"{state.projected_gradient_norm:.3e}"
+        )
+        LOGGER.info(
+            "PROGRESS scenario=%-11s iter=%6d objective=% .10e utility=% .10e "
+            "pg_norm=%s step=%.3e violation=%.3e elapsed=%.2fs",
+            name,
+            state.iteration,
+            state.objective,
+            state.utility,
+            residual,
+            state.step_size,
+            state.max_constraint_violation,
+            state.elapsed_seconds,
+        )
+
+    return report
+
+
+def write_history(path: Path, history: dict[str, list[float]]) -> None:
+    columns = list(history)
+    with path.open("w", newline="", encoding="utf-8") as stream:
+        writer = csv.writer(stream)
+        writer.writerow(columns)
+        writer.writerows(zip(*(history[column] for column in columns)))
+
+
+def solve_and_compare(
+    name: str,
+    problem: PortfolioProblem,
+    constraints: ConstraintSet,
+    options: PGDOptions,
+    output_dir: Path,
+    *,
+    exact_quadratic: bool,
+    objective_gap_tolerance: float,
+    holdings_distance_tolerance: float,
+) -> dict[str, object]:
+    LOGGER.info(
+        "START scenario=%s assets=%d initial_objective=% .10e initial_violation=%.3e",
+        name,
+        problem.dimension,
+        problem.value(constraints.project(problem.previous_holdings)),
+        constraints.max_violation(constraints.project(problem.previous_holdings)),
+    )
+    started = perf_counter()
+    pgd = solve_pgd(
+        problem,
+        constraints,
+        options=options,
+        progress_callback=progress_logger(name),
+    )
+    LOGGER.info("REFERENCE scenario=%s solver=scipy_slsqp status=running", name)
+    slsqp = solve_scipy_slsqp(problem, constraints, tolerance=1.0e-10)
+    if not slsqp.success:
+        raise RuntimeError(f"SciPy SLSQP failed for {name}: {slsqp.message}")
+
+    kkt = solve_quadratic_kkt(problem, constraints) if exact_quadratic else None
+    reference = kkt if kkt is not None else slsqp
+    holdings_distance = float(np.linalg.norm(pgd.holdings - reference.holdings))
+    objective_gap = float(pgd.objective - reference.objective)
+    turnover = float(np.sum(np.abs(pgd.trades)))
+    elapsed = perf_counter() - started
+    LOGGER.info(
+        "RESULT scenario=%s status=%s iterations=%d objective=% .10e reference=% .10e "
+        "objective_gap=% .3e holdings_distance=%.3e turnover=%.6f violation=%.3e elapsed=%.2fs",
+        name,
+        pgd.status,
+        pgd.iterations,
+        pgd.objective,
+        reference.objective,
+        objective_gap,
+        holdings_distance,
+        turnover,
+        pgd.max_constraint_violation,
+        elapsed,
+    )
+    if not pgd.converged:
+        raise RuntimeError(f"PGD failed to converge for {name}: {pgd.status}")
+    if pgd.max_constraint_violation > 2.0e-7:
+        raise RuntimeError(f"constraint violation too large for {name}")
+    if abs(objective_gap) > objective_gap_tolerance:
+        raise RuntimeError(
+            f"objective gap too large for {name}: {objective_gap:.3e} > "
+            f"{objective_gap_tolerance:.3e}"
+        )
+    if holdings_distance > holdings_distance_tolerance:
+        raise RuntimeError(
+            f"holdings distance too large for {name}: {holdings_distance:.3e} > "
+            f"{holdings_distance_tolerance:.3e}"
+        )
+
+    write_history(output_dir / f"{name}_history.csv", pgd.history)
+    np.savetxt(output_dir / f"{name}_holdings.csv", pgd.holdings, delimiter=",")
+    np.savetxt(output_dir / f"{name}_trades.csv", pgd.trades, delimiter=",")
+    return {
+        "scenario": name,
+        "converged": pgd.converged,
+        "status": pgd.status,
+        "iterations": pgd.iterations,
+        "objective": pgd.objective,
+        "utility": pgd.utility,
+        "reference_objective": reference.objective,
+        "objective_gap": objective_gap,
+        "holdings_distance_to_reference": holdings_distance,
+        "turnover": turnover,
+        "projected_gradient_norm": pgd.projected_gradient_norm,
+        "max_constraint_violation": pgd.max_constraint_violation,
+        "elapsed_seconds": elapsed,
+        "reference_solver": "exact_kkt" if kkt is not None else "scipy_slsqp",
+        "slsqp_message": slsqp.message,
+    }
+
+
+def quadratic_case(seed: int):
+    n = 36
+    covariance, loadings = factor_covariance(n, 4, seed=seed + 1, specific_risk=0.12)
+    rng = np.random.default_rng(seed + 2)
+    previous = rng.normal(scale=0.01, size=n)
+    problem = PortfolioProblem(
+        alpha=rng.normal(scale=0.025, size=n),
+        covariance=covariance,
+        previous_holdings=previous,
+        risk_aversion=2.25,
+        quadratic_cost_matrix=0.4 + rng.random(n),
+        quadratic_cost_aversion=0.75,
+    )
+    factor = loadings[:, 0] - np.mean(loadings[:, 0])
+    constraints = ConstraintSet(
+        n,
+        equality_matrix=np.vstack([np.ones(n), factor]),
+        equality_target=np.array([1.0, 0.0]),
+    )
+    return problem, constraints
+
+
+def nonlinear_case(seed: int):
+    n = 30
+    covariance, _ = factor_covariance(n, 4, seed=seed + 11, specific_risk=0.15)
+    rng = np.random.default_rng(seed + 12)
+    previous = capped_long_only_portfolio(n, cap=0.055, seed=seed + 13)
+    problem = PortfolioProblem(
+        alpha=rng.normal(scale=0.03, size=n),
+        covariance=covariance,
+        previous_holdings=previous,
+        risk_aversion=1.8,
+        quadratic_cost_matrix=0.2 + rng.random(n),
+        quadratic_cost_aversion=0.25,
+        nonlinear_cost=PowerLawCost(
+            eta=0.008 + 0.012 * rng.random(n), p=1.5, epsilon=1.0e-3
+        ),
+    )
+    constraints = ConstraintSet(
+        n,
+        equality_matrix=np.ones((1, n)),
+        equality_target=np.array([1.0]),
+        lower_bounds=0.0,
+        upper_bounds=0.075,
+    )
+    return problem, constraints
+
+
+def realistic_case(seed: int):
+    n = 40
+    covariance, loadings = factor_covariance(n, 5, seed=seed + 21, specific_risk=0.18)
+    rng = np.random.default_rng(seed + 22)
+    previous = capped_long_only_portfolio(n, cap=0.045, seed=seed + 23)
+    sectors = sector_membership(n, 5)
+    previous_sector = sectors @ previous
+    sector_lower = np.maximum(0.12, previous_sector - 0.035)
+    sector_upper = np.minimum(0.30, previous_sector + 0.035)
+    factor = loadings[:, 0] - np.mean(loadings[:, 0])
+    constraints = ConstraintSet(
+        n,
+        equality_matrix=np.vstack([np.ones(n), factor]),
+        equality_target=np.array([1.0, float(factor @ previous)]),
+        inequality_matrix=np.vstack([sectors, -sectors]),
+        inequality_upper=np.concatenate([sector_upper, -sector_lower]),
+        lower_bounds=0.0,
+        upper_bounds=0.06,
+        turnover_limit=0.22,
+        turnover_center=previous,
+    )
+    problem = PortfolioProblem(
+        alpha=rng.normal(scale=0.035, size=n),
+        covariance=covariance,
+        previous_holdings=previous,
+        risk_aversion=1.6,
+        quadratic_cost_matrix=0.2 + rng.random(n),
+        quadratic_cost_aversion=0.25,
+        nonlinear_cost=PowerLawCost(
+            eta=0.006 + 0.006 * rng.random(n), p=1.5, epsilon=1.0e-3
+        ),
+    )
+    return problem, constraints
+
+
+def main() -> int:
+    args = parse_args()
+    configure_logging()
+    args.output_dir.mkdir(parents=True, exist_ok=True)
+    total_stages = 5
+    started = perf_counter()
+
+    stage(1, total_stages, "runtime and configuration preflight")
+    LOGGER.info("experiment_dir=%s", EXPERIMENT_DIR)
+    LOGGER.info("output_dir=%s", args.output_dir)
+    LOGGER.info("python=%s", sys.version.replace("\n", " "))
+    LOGGER.info("platform=%s", platform.platform())
+    LOGGER.info("numpy=%s scipy=%s", np.__version__, scipy.__version__)
+    LOGGER.info(
+        "seed=%d tolerance=%.3e max_iterations=%d log_every=%d",
+        args.seed,
+        args.tolerance,
+        args.max_iterations,
+        args.log_every,
+    )
+
+    options = PGDOptions(
+        max_iterations=args.max_iterations,
+        tolerance=args.tolerance,
+        projection_tolerance=2.0e-10,
+        progress_interval=args.log_every,
+    )
+    summaries = []
+    cases = [
+        (2, "quadratic", quadratic_case(args.seed), True, 1.0e-9, 1.0e-6),
+        (3, "nonlinear", nonlinear_case(args.seed), False, 1.0e-6, 1.0e-3),
+        (4, "realistic", realistic_case(args.seed), False, 2.0e-5, 2.0e-3),
+    ]
+    for stage_number, name, (problem, constraints), exact, objective_tol, distance_tol in cases:
+        stage(stage_number, total_stages, f"{name} portfolio solve")
+        summaries.append(
+            solve_and_compare(
+                name,
+                problem,
+                constraints,
+                options,
+                args.output_dir,
+                exact_quadratic=exact,
+                objective_gap_tolerance=objective_tol,
+                holdings_distance_tolerance=distance_tol,
+            )
+        )
+
+    stage(5, total_stages, "write manifest and final acceptance summary")
+    manifest = {
+        "completed_utc": datetime.now(timezone.utc).isoformat(),
+        "experiment": "portfolio_pgd",
+        "seed": args.seed,
+        "options": asdict(options),
+        "python": sys.version,
+        "platform": platform.platform(),
+        "numpy": np.__version__,
+        "scipy": scipy.__version__,
+        "scenarios": summaries,
+        "elapsed_seconds": perf_counter() - started,
+    }
+    (args.output_dir / "summary.json").write_text(
+        json.dumps(manifest, indent=2, sort_keys=True) + "\n", encoding="utf-8"
+    )
+    with (args.output_dir / "summary.csv").open("w", newline="", encoding="utf-8") as stream:
+        writer = csv.DictWriter(stream, fieldnames=list(summaries[0]))
+        writer.writeheader()
+        writer.writerows(summaries)
+    for summary in summaries:
+        LOGGER.info(
+            "ACCEPT scenario=%s converged=%s objective_gap=% .3e distance=%.3e violation=%.3e",
+            summary["scenario"],
+            summary["converged"],
+            summary["objective_gap"],
+            summary["holdings_distance_to_reference"],
+            summary["max_constraint_violation"],
+        )
+    LOGGER.info("COMPLETE output_dir=%s elapsed=%.2fs", args.output_dir, manifest["elapsed_seconds"])
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/baseline/experiments/portfolio_pgd/scripts/run_portfolio_pgd_experiment.sh b/baseline/experiments/portfolio_pgd/scripts/run_portfolio_pgd_experiment.sh
new file mode 100755
index 00000000..9c9f181f
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/scripts/run_portfolio_pgd_experiment.sh
@@ -0,0 +1,74 @@
+#!/usr/bin/env bash
+
+# Foreground, terminal-safe campaign runner. It never backgrounds, detaches,
+# kills, or replaces the caller's shell. All output is streamed and persisted.
+
+set -uo pipefail
+
+SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
+EXPERIMENT_DIR="$(cd "$SCRIPT_DIR/.." && pwd)"
+ENV_PREFIX="${PORTFOLIO_PGD_ENV:-/tmp/portfolio-pgd-env}"
+RUN_ROOT="${PORTFOLIO_PGD_RUN_ROOT:-/tmp/portfolio-pgd-runs}"
+RUN_ID="${PORTFOLIO_PGD_RUN_ID:-$(date '+%Y%m%d_%H%M%S')_$$}"
+RUN_DIR="$RUN_ROOT/$RUN_ID"
+LOG_FILE="$RUN_DIR/run.log"
+
+mkdir -p "$RUN_DIR"
+
+timestamp() { date '+%Y-%m-%d %H:%M:%S'; }
+log() {
+  printf '%s %-8s %s\n' "$(timestamp)" "$1" "$2" | tee -a "$LOG_FILE"
+}
+
+run_stage() {
+  local stage_number="$1"
+  local title="$2"
+  shift 2
+  log INFO "========================================================================"
+  log INFO "CAMPAIGN STAGE $stage_number/3 $title"
+  log INFO "COMMAND $*"
+  log INFO "========================================================================"
+  "$@" 2>&1 | tee -a "$LOG_FILE"
+  local command_status="${PIPESTATUS[0]}"
+  if [[ "$command_status" -ne 0 ]]; then
+    log ERROR "FAILED stage=$stage_number status=$command_status"
+    return "$command_status"
+  fi
+  log INFO "PASSED stage=$stage_number"
+  return 0
+}
+
+log INFO "START portfolio-PGD complete campaign"
+log INFO "experiment_dir=$EXPERIMENT_DIR"
+log INFO "environment=$ENV_PREFIX"
+log INFO "run_dir=$RUN_DIR"
+log INFO "log_file=$LOG_FILE"
+
+if [[ ! -x "$ENV_PREFIX/bin/python" ]]; then
+  log ERROR "Missing environment. Run: bash $SCRIPT_DIR/setup_mac.sh"
+  exit 1
+fi
+
+PYTHON=(conda run --no-capture-output -p "$ENV_PREFIX" python -u)
+
+if ! run_stage 1 "automated unit and reference-solver tests" \
+  "${PYTHON[@]}" "$EXPERIMENT_DIR/tests/run_tests.py"; then
+  exit 1
+fi
+
+if ! run_stage 2 "logged quadratic, nonlinear, and realistic PGD solves" \
+  "${PYTHON[@]}" "$SCRIPT_DIR/run_portfolio_pgd_experiment.py" \
+  --output-dir "$RUN_DIR/metrics" \
+  --seed "${PORTFOLIO_PGD_SEED:-20260822}" \
+  --log-every "${PORTFOLIO_PGD_LOG_EVERY:-25}"; then
+  exit 1
+fi
+
+if ! run_stage 3 "execute all demonstration notebook cells" \
+  "${PYTHON[@]}" "$SCRIPT_DIR/execute_notebooks.py"; then
+  exit 1
+fi
+
+log INFO "COMPLETE all stages passed"
+log INFO "summary=$RUN_DIR/metrics/summary.json"
+log INFO "terminal remains active; this script did not detach or kill any process"
diff --git a/baseline/experiments/portfolio_pgd/scripts/setup_mac.sh b/baseline/experiments/portfolio_pgd/scripts/setup_mac.sh
new file mode 100755
index 00000000..995adeef
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/scripts/setup_mac.sh
@@ -0,0 +1,54 @@
+#!/usr/bin/env bash
+
+# Create the complete disposable Mac environment under literal /tmp.
+# Run as: bash baseline/experiments/portfolio_pgd/scripts/setup_mac.sh
+
+set -uo pipefail
+
+SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
+EXPERIMENT_DIR="$(cd "$SCRIPT_DIR/.." && pwd)"
+ENV_PREFIX="${PORTFOLIO_PGD_ENV:-/tmp/portfolio-pgd-env}"
+
+timestamp() { date '+%Y-%m-%d %H:%M:%S'; }
+log() { printf '%s %-8s %s\n' "$(timestamp)" "$1" "$2"; }
+
+log INFO "START portfolio-PGD Mac setup"
+log INFO "experiment_dir=$EXPERIMENT_DIR"
+log INFO "environment=$ENV_PREFIX"
+
+if ! command -v conda >/dev/null 2>&1; then
+  log ERROR "conda was not found on PATH"
+  exit 1
+fi
+
+if [[ -x "$ENV_PREFIX/bin/python" ]]; then
+  log INFO "REUSE existing environment"
+else
+  log INFO "CREATE conda environment with Python 3.11"
+  if ! conda create -p "$ENV_PREFIX" python=3.11 -y; then
+    log ERROR "conda environment creation failed"
+    exit 1
+  fi
+fi
+
+log INFO "INSTALL editable experiment and notebook dependencies"
+if ! conda run --no-capture-output -p "$ENV_PREFIX" \
+  python -m pip install -e "$EXPERIMENT_DIR[notebook]"; then
+  log ERROR "dependency installation failed"
+  exit 1
+fi
+
+log INFO "VERIFY imports"
+if ! conda run --no-capture-output -p "$ENV_PREFIX" python - <<'PY'
+import numpy
+import scipy
+import portfolio_pgd
+print(f"portfolio_pgd={portfolio_pgd.__version__} numpy={numpy.__version__} scipy={scipy.__version__}")
+PY
+then
+  log ERROR "import verification failed"
+  exit 1
+fi
+
+log INFO "COMPLETE setup succeeded"
+log INFO "Next: bash baseline/experiments/portfolio_pgd/scripts/run_portfolio_pgd_experiment.sh"
diff --git a/baseline/experiments/portfolio_pgd/src/portfolio_pgd/__init__.py b/baseline/experiments/portfolio_pgd/src/portfolio_pgd/__init__.py
new file mode 100644
index 00000000..eac698cf
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/src/portfolio_pgd/__init__.py
@@ -0,0 +1,29 @@
+"""Projected-gradient portfolio optimization."""
+
+from .constraints import ConstraintSet, ProjectionError
+from .costs import PowerLawCost, SmoothAbsoluteCost, TransactionCost
+from .problem import PortfolioProblem
+from .reference import ReferenceResult, solve_quadratic_kkt, solve_scipy_slsqp
+from .solver import PGDOptions, ProgressState, SolverResult, solve_pgd
+from .synthetic import capped_long_only_portfolio, factor_covariance, sector_membership
+
+__all__ = [
+    "ConstraintSet",
+    "PGDOptions",
+    "PortfolioProblem",
+    "PowerLawCost",
+    "ProjectionError",
+    "ProgressState",
+    "ReferenceResult",
+    "SmoothAbsoluteCost",
+    "SolverResult",
+    "TransactionCost",
+    "capped_long_only_portfolio",
+    "factor_covariance",
+    "sector_membership",
+    "solve_pgd",
+    "solve_quadratic_kkt",
+    "solve_scipy_slsqp",
+]
+
+__version__ = "1.0.0"
diff --git a/baseline/experiments/portfolio_pgd/src/portfolio_pgd/constraints.py b/baseline/experiments/portfolio_pgd/src/portfolio_pgd/constraints.py
new file mode 100644
index 00000000..1948b9c7
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/src/portfolio_pgd/constraints.py
@@ -0,0 +1,259 @@
+"""Convex portfolio constraints and Euclidean projection by Dykstra's algorithm."""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+import numpy as np
+from numpy.typing import ArrayLike, NDArray
+
+FloatArray = NDArray[np.float64]
+
+
+class ProjectionError(RuntimeError):
+    """Raised when the supplied convex constraints appear infeasible or fail to project."""
+
+
+def _vector(value: ArrayLike | None, size: int, name: str, default: float) -> FloatArray:
+    if value is None:
+        return np.full(size, default, dtype=float)
+    array = np.asarray(value, dtype=float)
+    if array.ndim == 0:
+        array = np.full(size, float(array), dtype=float)
+    if array.shape != (size,):
+        raise ValueError(f"{name} must be scalar or have shape ({size},)")
+    return array
+
+
+def project_l1_ball(vector: ArrayLike, radius: float, center: ArrayLike | None = None) -> FloatArray:
+    """Project a vector onto ``{x: ||x-center||_1 <= radius}``."""
+    x = np.asarray(vector, dtype=float)
+    if radius < 0.0:
+        raise ValueError("radius must be nonnegative")
+    c = np.zeros_like(x) if center is None else np.asarray(center, dtype=float)
+    if c.shape != x.shape:
+        raise ValueError("center and vector must have the same shape")
+    shifted = x - c
+    absolute = np.abs(shifted)
+    if float(np.sum(absolute)) <= radius:
+        return x.copy()
+    if radius == 0.0:
+        return c.copy()
+    ordered = np.sort(absolute)[::-1]
+    cumulative = np.cumsum(ordered)
+    indices = np.arange(1, ordered.size + 1)
+    active = np.nonzero(ordered - (cumulative - radius) / indices > 0.0)[0]
+    if active.size == 0:
+        return c.copy()
+    rho = int(active[-1])
+    threshold = (cumulative[rho] - radius) / float(rho + 1)
+    return c + np.sign(shifted) * np.maximum(absolute - threshold, 0.0)
+
+
+@dataclass(frozen=True)
+class ConstraintSet:
+    """Intersection of common convex institutional portfolio constraints.
+
+    Supported constraints are affine equalities, linear inequalities, per-name
+    lower/upper bounds, turnover around a reference portfolio, and gross
+    exposure.  Projection onto their intersection uses Dykstra's algorithm;
+    individual projections are analytic.
+    """
+
+    dimension: int
+    equality_matrix: ArrayLike | None = None
+    equality_target: ArrayLike | None = None
+    inequality_matrix: ArrayLike | None = None
+    inequality_upper: ArrayLike | None = None
+    lower_bounds: ArrayLike | None = None
+    upper_bounds: ArrayLike | None = None
+    turnover_limit: float | None = None
+    turnover_center: ArrayLike | None = None
+    gross_exposure_limit: float | None = None
+
+    def __post_init__(self) -> None:
+        n = int(self.dimension)
+        if n <= 0:
+            raise ValueError("dimension must be positive")
+
+        if self.equality_matrix is None:
+            a_eq = np.zeros((0, n), dtype=float)
+            b_eq = np.zeros(0, dtype=float)
+        else:
+            a_eq = np.atleast_2d(np.asarray(self.equality_matrix, dtype=float))
+            if a_eq.shape[1] != n:
+                raise ValueError(f"equality_matrix must have {n} columns")
+            if self.equality_target is None:
+                raise ValueError("equality_target is required with equality_matrix")
+            b_eq = np.atleast_1d(np.asarray(self.equality_target, dtype=float))
+            if b_eq.shape != (a_eq.shape[0],):
+                raise ValueError("equality_target has incompatible shape")
+
+        if self.inequality_matrix is None:
+            a_ub = np.zeros((0, n), dtype=float)
+            b_ub = np.zeros(0, dtype=float)
+        else:
+            a_ub = np.atleast_2d(np.asarray(self.inequality_matrix, dtype=float))
+            if a_ub.shape[1] != n:
+                raise ValueError(f"inequality_matrix must have {n} columns")
+            if self.inequality_upper is None:
+                raise ValueError("inequality_upper is required with inequality_matrix")
+            b_ub = np.atleast_1d(np.asarray(self.inequality_upper, dtype=float))
+            if b_ub.shape != (a_ub.shape[0],):
+                raise ValueError("inequality_upper has incompatible shape")
+
+        lower = _vector(self.lower_bounds, n, "lower_bounds", -np.inf)
+        upper = _vector(self.upper_bounds, n, "upper_bounds", np.inf)
+        if np.any(lower > upper):
+            raise ValueError("lower_bounds must not exceed upper_bounds")
+
+        if self.turnover_limit is not None and self.turnover_limit < 0.0:
+            raise ValueError("turnover_limit must be nonnegative")
+        if self.turnover_limit is not None:
+            if self.turnover_center is None:
+                raise ValueError("turnover_center is required with turnover_limit")
+            center = np.asarray(self.turnover_center, dtype=float)
+            if center.shape != (n,):
+                raise ValueError(f"turnover_center must have shape ({n},)")
+        else:
+            center = np.zeros(n, dtype=float)
+        if self.gross_exposure_limit is not None and self.gross_exposure_limit < 0.0:
+            raise ValueError("gross_exposure_limit must be nonnegative")
+
+        for name, array in (("equality_matrix", a_eq), ("equality_target", b_eq),
+                            ("inequality_matrix", a_ub), ("inequality_upper", b_ub)):
+            if np.any(~np.isfinite(array)):
+                raise ValueError(f"{name} must contain finite values")
+
+        object.__setattr__(self, "dimension", n)
+        object.__setattr__(self, "equality_matrix", a_eq)
+        object.__setattr__(self, "equality_target", b_eq)
+        object.__setattr__(self, "inequality_matrix", a_ub)
+        object.__setattr__(self, "inequality_upper", b_ub)
+        object.__setattr__(self, "lower_bounds", lower)
+        object.__setattr__(self, "upper_bounds", upper)
+        object.__setattr__(self, "turnover_center", center)
+
+    @property
+    def has_only_equalities(self) -> bool:
+        return bool(
+            np.asarray(self.inequality_matrix).shape[0] == 0
+            and np.all(np.isneginf(np.asarray(self.lower_bounds)))
+            and np.all(np.isposinf(np.asarray(self.upper_bounds)))
+            and self.turnover_limit is None
+            and self.gross_exposure_limit is None
+        )
+
+    def _projectors(self):
+        projectors = []
+        a_eq = np.asarray(self.equality_matrix)
+        b_eq = np.asarray(self.equality_target)
+        if a_eq.shape[0]:
+            gram_pinv = np.linalg.pinv(a_eq @ a_eq.T, rcond=1.0e-13)
+            correction = a_eq.T @ gram_pinv
+
+            def affine(x, a=a_eq, b=b_eq, k=correction):
+                return x - k @ (a @ x - b)
+
+            projectors.append(affine)
+
+        lower = np.asarray(self.lower_bounds)
+        upper = np.asarray(self.upper_bounds)
+        if np.any(np.isfinite(lower)) or np.any(np.isfinite(upper)):
+            projectors.append(lambda x, lo=lower, hi=upper: np.clip(x, lo, hi))
+
+        a_ub = np.asarray(self.inequality_matrix)
+        b_ub = np.asarray(self.inequality_upper)
+        for row, bound in zip(a_ub, b_ub):
+            norm_squared = float(row @ row)
+            if norm_squared == 0.0:
+                if bound < 0.0:
+                    raise ProjectionError("infeasible zero-row inequality")
+                continue
+
+            def halfspace(x, a=row.copy(), b=float(bound), denom=norm_squared):
+                excess = float(a @ x - b)
+                return x if excess <= 0.0 else x - (excess / denom) * a
+
+            projectors.append(halfspace)
+
+        if self.turnover_limit is not None:
+            center = np.asarray(self.turnover_center)
+            radius = float(self.turnover_limit)
+            projectors.append(lambda x, r=radius, c=center: project_l1_ball(x, r, c))
+
+        if self.gross_exposure_limit is not None:
+            radius = float(self.gross_exposure_limit)
+            projectors.append(lambda x, r=radius: project_l1_ball(x, r))
+        return projectors
+
+    def project(
+        self,
+        vector: ArrayLike,
+        *,
+        tolerance: float = 1.0e-10,
+        max_cycles: int = 5_000,
+        check_feasible: bool = True,
+    ) -> FloatArray:
+        """Return the Euclidean projection onto the full constraint intersection."""
+        z = np.asarray(vector, dtype=float)
+        if z.shape != (self.dimension,):
+            raise ValueError(f"vector must have shape ({self.dimension},)")
+        projectors = self._projectors()
+        if not projectors:
+            return z.copy()
+        x = z.copy()
+        corrections = [np.zeros_like(x) for _ in projectors]
+        converged = False
+        for cycle in range(max_cycles):
+            before = x.copy()
+            for index, projector in enumerate(projectors):
+                shifted = x + corrections[index]
+                projected = np.asarray(projector(shifted), dtype=float)
+                corrections[index] = shifted - projected
+                x = projected
+            scale = max(1.0, float(np.linalg.norm(x, ord=np.inf)))
+            delta = float(np.linalg.norm(x - before, ord=np.inf))
+            if delta <= tolerance * scale and self.max_violation(x) <= max(10.0 * tolerance, 1.0e-11):
+                converged = True
+                break
+        violation = self.max_violation(x)
+        if check_feasible and (not converged or violation > max(100.0 * tolerance, 1.0e-8)):
+            raise ProjectionError(
+                f"projection failed after {max_cycles} cycles; maximum violation={violation:.3e}. "
+                "The constraints may be infeasible."
+            )
+        return x
+
+    def violations(self, holdings: ArrayLike) -> dict[str, float]:
+        h = np.asarray(holdings, dtype=float)
+        if h.shape != (self.dimension,):
+            raise ValueError(f"holdings must have shape ({self.dimension},)")
+        a_eq = np.asarray(self.equality_matrix)
+        a_ub = np.asarray(self.inequality_matrix)
+        values: dict[str, float] = {
+            "equality": float(np.max(np.abs(a_eq @ h - np.asarray(self.equality_target))))
+            if a_eq.shape[0]
+            else 0.0,
+            "linear_inequality": float(
+                max(0.0, np.max(a_ub @ h - np.asarray(self.inequality_upper)))
+            )
+            if a_ub.shape[0]
+            else 0.0,
+            "lower_bound": float(max(0.0, np.max(np.asarray(self.lower_bounds) - h))),
+            "upper_bound": float(max(0.0, np.max(h - np.asarray(self.upper_bounds)))),
+            "turnover": 0.0,
+            "gross_exposure": 0.0,
+        }
+        if self.turnover_limit is not None:
+            values["turnover"] = float(
+                max(0.0, np.sum(np.abs(h - np.asarray(self.turnover_center))) - self.turnover_limit)
+            )
+        if self.gross_exposure_limit is not None:
+            values["gross_exposure"] = float(
+                max(0.0, np.sum(np.abs(h)) - self.gross_exposure_limit)
+            )
+        return values
+
+    def max_violation(self, holdings: ArrayLike) -> float:
+        return max(self.violations(holdings).values())
diff --git a/baseline/experiments/portfolio_pgd/src/portfolio_pgd/costs.py b/baseline/experiments/portfolio_pgd/src/portfolio_pgd/costs.py
new file mode 100644
index 00000000..94c216b6
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/src/portfolio_pgd/costs.py
@@ -0,0 +1,126 @@
+"""Differentiable convex transaction-cost models."""
+
+from __future__ import annotations
+
+from abc import ABC, abstractmethod
+from dataclasses import dataclass
+
+import numpy as np
+from numpy.typing import ArrayLike, NDArray
+
+FloatArray = NDArray[np.float64]
+
+
+class TransactionCost(ABC):
+    """Interface for an additive convex cost applied to the trade vector."""
+
+    @abstractmethod
+    def value(self, trades: FloatArray) -> float:
+        """Return the total transaction cost."""
+
+    @abstractmethod
+    def gradient(self, trades: FloatArray) -> FloatArray:
+        """Return the gradient with respect to trades."""
+
+    def hessian_diag(self, trades: FloatArray) -> FloatArray:
+        """Return a diagonal Hessian approximation, used only for step initialization."""
+        return np.zeros_like(trades, dtype=float)
+
+
+def _broadcast_parameter(value: ArrayLike, size: int, name: str) -> FloatArray:
+    array = np.asarray(value, dtype=float)
+    if array.ndim == 0:
+        array = np.full(size, float(array))
+    if array.shape != (size,):
+        raise ValueError(f"{name} must be scalar or have shape ({size},)")
+    if np.any(~np.isfinite(array)):
+        raise ValueError(f"{name} must contain finite values")
+    return array
+
+
+@dataclass(frozen=True)
+class PowerLawCost(TransactionCost):
+    r"""Separable cost :math:`\sum_i \eta_i |t_i|^p` with optional smoothing.
+
+    With ``epsilon > 0`` the implementation uses
+    ``eta * ((t**2 + epsilon**2)**(p/2) - epsilon**p)``.  This keeps the
+    objective smooth and convex for every ``p > 1`` while preserving zero cost
+    at zero trade.
+    """
+
+    eta: ArrayLike
+    p: float = 1.5
+    epsilon: float = 0.0
+
+    def _eta(self, size: int) -> FloatArray:
+        eta = _broadcast_parameter(self.eta, size, "eta")
+        if np.any(eta < 0.0):
+            raise ValueError("eta must be nonnegative")
+        if self.p <= 1.0:
+            raise ValueError("p must be greater than one for a differentiable convex cost")
+        if self.epsilon < 0.0:
+            raise ValueError("epsilon must be nonnegative")
+        return eta
+
+    def value(self, trades: FloatArray) -> float:
+        trades = np.asarray(trades, dtype=float)
+        eta = self._eta(trades.size)
+        if self.epsilon == 0.0:
+            return float(np.sum(eta * np.abs(trades) ** self.p))
+        radius2 = trades * trades + self.epsilon * self.epsilon
+        return float(np.sum(eta * (radius2 ** (0.5 * self.p) - self.epsilon**self.p)))
+
+    def gradient(self, trades: FloatArray) -> FloatArray:
+        trades = np.asarray(trades, dtype=float)
+        eta = self._eta(trades.size)
+        if self.epsilon == 0.0:
+            return self.p * eta * np.abs(trades) ** (self.p - 1.0) * np.sign(trades)
+        radius2 = trades * trades + self.epsilon * self.epsilon
+        return self.p * eta * trades * radius2 ** (0.5 * self.p - 1.0)
+
+    def hessian_diag(self, trades: FloatArray) -> FloatArray:
+        trades = np.asarray(trades, dtype=float)
+        eta = self._eta(trades.size)
+        if self.epsilon == 0.0:
+            magnitude = np.abs(trades)
+            with np.errstate(divide="ignore", invalid="ignore"):
+                diagonal = self.p * (self.p - 1.0) * eta * magnitude ** (self.p - 2.0)
+            return np.nan_to_num(diagonal, nan=0.0, posinf=np.finfo(float).max ** 0.25)
+        radius2 = trades * trades + self.epsilon * self.epsilon
+        return (
+            self.p
+            * eta
+            * radius2 ** (0.5 * self.p - 2.0)
+            * (self.epsilon * self.epsilon + (self.p - 1.0) * trades * trades)
+        )
+
+
+@dataclass(frozen=True)
+class SmoothAbsoluteCost(TransactionCost):
+    r"""Smooth bid-ask cost ``rate * (sqrt(t**2 + epsilon**2) - epsilon)``."""
+
+    rate: ArrayLike
+    epsilon: float = 1.0e-4
+
+    def _rate(self, size: int) -> FloatArray:
+        rate = _broadcast_parameter(self.rate, size, "rate")
+        if np.any(rate < 0.0):
+            raise ValueError("rate must be nonnegative")
+        if self.epsilon <= 0.0:
+            raise ValueError("epsilon must be strictly positive")
+        return rate
+
+    def value(self, trades: FloatArray) -> float:
+        trades = np.asarray(trades, dtype=float)
+        rate = self._rate(trades.size)
+        return float(np.sum(rate * (np.sqrt(trades * trades + self.epsilon**2) - self.epsilon)))
+
+    def gradient(self, trades: FloatArray) -> FloatArray:
+        trades = np.asarray(trades, dtype=float)
+        rate = self._rate(trades.size)
+        return rate * trades / np.sqrt(trades * trades + self.epsilon**2)
+
+    def hessian_diag(self, trades: FloatArray) -> FloatArray:
+        trades = np.asarray(trades, dtype=float)
+        rate = self._rate(trades.size)
+        return rate * self.epsilon**2 / (trades * trades + self.epsilon**2) ** 1.5
diff --git a/baseline/experiments/portfolio_pgd/src/portfolio_pgd/problem.py b/baseline/experiments/portfolio_pgd/src/portfolio_pgd/problem.py
new file mode 100644
index 00000000..32524097
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/src/portfolio_pgd/problem.py
@@ -0,0 +1,143 @@
+"""Portfolio objective and validation."""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+import numpy as np
+from numpy.typing import ArrayLike, NDArray
+
+from .costs import TransactionCost
+
+FloatArray = NDArray[np.float64]
+
+
+@dataclass(frozen=True)
+class PortfolioProblem:
+    r"""Single-period convex portfolio-construction problem.
+
+    The minimized objective is
+
+    .. math::
+
+        \tfrac{\lambda}{2}h^T Vh - \alpha^T h
+        + \tfrac{\theta}{2}(h-h_-)^TQ(h-h_-)
+        + c(h-h_-).
+    """
+
+    alpha: ArrayLike
+    covariance: ArrayLike
+    previous_holdings: ArrayLike
+    risk_aversion: float = 1.0
+    quadratic_cost_matrix: ArrayLike | None = None
+    quadratic_cost_aversion: float = 0.0
+    nonlinear_cost: TransactionCost | None = None
+
+    def __post_init__(self) -> None:
+        alpha = np.asarray(self.alpha, dtype=float)
+        covariance = np.asarray(self.covariance, dtype=float)
+        previous = np.asarray(self.previous_holdings, dtype=float)
+        if alpha.ndim != 1:
+            raise ValueError("alpha must be one-dimensional")
+        n = alpha.size
+        if covariance.shape != (n, n):
+            raise ValueError(f"covariance must have shape ({n}, {n})")
+        if previous.shape != (n,):
+            raise ValueError(f"previous_holdings must have shape ({n},)")
+        if self.risk_aversion <= 0.0:
+            raise ValueError("risk_aversion must be strictly positive")
+        if self.quadratic_cost_aversion < 0.0:
+            raise ValueError("quadratic_cost_aversion must be nonnegative")
+        if np.any(~np.isfinite(alpha)) or np.any(~np.isfinite(covariance)):
+            raise ValueError("problem data must be finite")
+        if not np.allclose(covariance, covariance.T, atol=1.0e-12):
+            raise ValueError("covariance must be symmetric")
+        try:
+            np.linalg.cholesky(covariance)
+        except np.linalg.LinAlgError as exc:
+            raise ValueError("covariance must be positive definite") from exc
+
+        if self.quadratic_cost_matrix is None:
+            q = np.zeros((n, n), dtype=float)
+        else:
+            raw_q = np.asarray(self.quadratic_cost_matrix, dtype=float)
+            if raw_q.ndim == 1:
+                if raw_q.shape != (n,):
+                    raise ValueError(f"quadratic cost diagonal must have shape ({n},)")
+                q = np.diag(raw_q)
+            else:
+                q = raw_q
+            if q.shape != (n, n):
+                raise ValueError(f"quadratic_cost_matrix must have shape ({n}, {n})")
+            if not np.allclose(q, q.T, atol=1.0e-12):
+                raise ValueError("quadratic_cost_matrix must be symmetric")
+            if np.min(np.linalg.eigvalsh(q)) < -1.0e-12:
+                raise ValueError("quadratic_cost_matrix must be positive semidefinite")
+
+        object.__setattr__(self, "alpha", alpha)
+        object.__setattr__(self, "covariance", covariance)
+        object.__setattr__(self, "previous_holdings", previous)
+        object.__setattr__(self, "quadratic_cost_matrix", q)
+
+    @property
+    def dimension(self) -> int:
+        return int(np.asarray(self.alpha).size)
+
+    @property
+    def quadratic_hessian(self) -> FloatArray:
+        return (
+            self.risk_aversion * np.asarray(self.covariance)
+            + self.quadratic_cost_aversion * np.asarray(self.quadratic_cost_matrix)
+        )
+
+    @property
+    def quadratic_linear_term(self) -> FloatArray:
+        return np.asarray(self.alpha) + self.quadratic_cost_aversion * (
+            np.asarray(self.quadratic_cost_matrix) @ np.asarray(self.previous_holdings)
+        )
+
+    def value(self, holdings: ArrayLike) -> float:
+        h = np.asarray(holdings, dtype=float)
+        if h.shape != (self.dimension,):
+            raise ValueError(f"holdings must have shape ({self.dimension},)")
+        trades = h - np.asarray(self.previous_holdings)
+        value = (
+            0.5 * self.risk_aversion * float(h @ np.asarray(self.covariance) @ h)
+            - float(np.asarray(self.alpha) @ h)
+            + 0.5
+            * self.quadratic_cost_aversion
+            * float(trades @ np.asarray(self.quadratic_cost_matrix) @ trades)
+        )
+        if self.nonlinear_cost is not None:
+            value += self.nonlinear_cost.value(trades)
+        return float(value)
+
+    def utility(self, holdings: ArrayLike) -> float:
+        """Return utility, the negative of the minimized objective."""
+        return -self.value(holdings)
+
+    def gradient(self, holdings: ArrayLike) -> FloatArray:
+        h = np.asarray(holdings, dtype=float)
+        trades = h - np.asarray(self.previous_holdings)
+        gradient = (
+            self.risk_aversion * (np.asarray(self.covariance) @ h)
+            - np.asarray(self.alpha)
+            + self.quadratic_cost_aversion
+            * (np.asarray(self.quadratic_cost_matrix) @ trades)
+        )
+        if self.nonlinear_cost is not None:
+            gradient = gradient + self.nonlinear_cost.gradient(trades)
+        return np.asarray(gradient, dtype=float)
+
+    def local_lipschitz_bound(self, holdings: ArrayLike) -> float:
+        """Return a conservative infinity-norm bound on local Hessian size."""
+        h = np.asarray(holdings, dtype=float)
+        base = self.quadratic_hessian
+        bound = float(np.max(np.sum(np.abs(base), axis=1)))
+        if self.nonlinear_cost is not None:
+            trades = h - np.asarray(self.previous_holdings)
+            diagonal = self.nonlinear_cost.hessian_diag(trades)
+            finite = diagonal[np.isfinite(diagonal)]
+            if finite.size:
+                bound += float(max(0.0, np.max(finite)))
+        return max(bound, np.finfo(float).eps)
diff --git a/baseline/experiments/portfolio_pgd/src/portfolio_pgd/reference.py b/baseline/experiments/portfolio_pgd/src/portfolio_pgd/reference.py
new file mode 100644
index 00000000..fe9a98df
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/src/portfolio_pgd/reference.py
@@ -0,0 +1,201 @@
+"""Exact quadratic and SciPy reference solvers used for validation."""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import Any
+
+import numpy as np
+from numpy.typing import ArrayLike, NDArray
+from scipy.optimize import Bounds, minimize
+
+from .constraints import ConstraintSet
+from .problem import PortfolioProblem
+
+FloatArray = NDArray[np.float64]
+
+
+@dataclass
+class ReferenceResult:
+    holdings: FloatArray
+    objective: float
+    success: bool
+    message: str
+    multipliers: FloatArray | None = None
+    raw_result: Any | None = None
+
+
+def solve_quadratic_kkt(
+    problem: PortfolioProblem,
+    constraints: ConstraintSet | None = None,
+) -> ReferenceResult:
+    """Solve the quadratic, equality-constrained problem through its KKT system."""
+    constraints = constraints or ConstraintSet(problem.dimension)
+    if problem.nonlinear_cost is not None:
+        raise ValueError("KKT solver supports quadratic costs only")
+    if not constraints.has_only_equalities:
+        raise ValueError("KKT solver supports affine equality constraints only")
+    hessian = problem.quadratic_hessian
+    linear = problem.quadratic_linear_term
+    a_eq = np.asarray(constraints.equality_matrix)
+    b_eq = np.asarray(constraints.equality_target)
+    if a_eq.shape[0] == 0:
+        holdings = np.linalg.solve(hessian, linear)
+        multipliers = np.zeros(0, dtype=float)
+    else:
+        kkt = np.block(
+            [
+                [hessian, a_eq.T],
+                [a_eq, np.zeros((a_eq.shape[0], a_eq.shape[0]))],
+            ]
+        )
+        rhs = np.concatenate([linear, b_eq])
+        try:
+            solution = np.linalg.solve(kkt, rhs)
+        except np.linalg.LinAlgError:
+            solution, residuals, rank, _ = np.linalg.lstsq(kkt, rhs, rcond=1.0e-13)
+            if rank < kkt.shape[0] and residuals.size and float(np.max(residuals)) > 1.0e-14:
+                raise ValueError("singular or inconsistent KKT system")
+        holdings = solution[: problem.dimension]
+        multipliers = solution[problem.dimension :]
+    return ReferenceResult(
+        holdings=np.asarray(holdings, dtype=float),
+        objective=problem.value(holdings),
+        success=True,
+        message="exact KKT solution",
+        multipliers=np.asarray(multipliers, dtype=float),
+    )
+
+
+def solve_scipy_slsqp(
+    problem: PortfolioProblem,
+    constraints: ConstraintSet | None = None,
+    *,
+    initial_holdings: ArrayLike | None = None,
+    max_iterations: int = 5_000,
+    tolerance: float = 1.0e-11,
+) -> ReferenceResult:
+    """Solve the same problem with SciPy SLSQP as an independent benchmark."""
+    constraints = constraints or ConstraintSet(problem.dimension)
+    initial = problem.previous_holdings if initial_holdings is None else np.asarray(initial_holdings, dtype=float)
+    initial = constraints.project(initial, tolerance=1.0e-10, max_cycles=10_000)
+    # Lift L1 constraints with auxiliary variables.  This makes the reference
+    # problem smooth with purely linear constraints instead of asking SLSQP to
+    # finite-difference an absolute value at zero.
+    n = problem.dimension
+    turnover_slice: slice | None = None
+    gross_slice: slice | None = None
+    total_dimension = n
+    if constraints.turnover_limit is not None:
+        turnover_slice = slice(total_dimension, total_dimension + n)
+        total_dimension += n
+    if constraints.gross_exposure_limit is not None:
+        gross_slice = slice(total_dimension, total_dimension + n)
+        total_dimension += n
+
+    lifted_initial = np.zeros(total_dimension, dtype=float)
+    lifted_initial[:n] = initial
+    if turnover_slice is not None:
+        lifted_initial[turnover_slice] = np.abs(initial - np.asarray(constraints.turnover_center))
+    if gross_slice is not None:
+        lifted_initial[gross_slice] = np.abs(initial)
+
+    def lifted_value(vector: FloatArray) -> float:
+        return problem.value(vector[:n])
+
+    def lifted_gradient(vector: FloatArray) -> FloatArray:
+        gradient = np.zeros(total_dimension, dtype=float)
+        gradient[:n] = problem.gradient(vector[:n])
+        return gradient
+
+    equality_rows: list[FloatArray] = []
+    equality_targets: list[FloatArray] = []
+    inequality_rows: list[FloatArray] = []
+    inequality_targets: list[FloatArray] = []
+    a_eq = np.asarray(constraints.equality_matrix)
+    b_eq = np.asarray(constraints.equality_target)
+    if a_eq.shape[0]:
+        lifted = np.zeros((a_eq.shape[0], total_dimension), dtype=float)
+        lifted[:, :n] = a_eq
+        equality_rows.append(lifted)
+        equality_targets.append(b_eq)
+    a_ub = np.asarray(constraints.inequality_matrix)
+    b_ub = np.asarray(constraints.inequality_upper)
+    if a_ub.shape[0]:
+        lifted = np.zeros((a_ub.shape[0], total_dimension), dtype=float)
+        lifted[:, :n] = a_ub
+        inequality_rows.append(lifted)
+        inequality_targets.append(b_ub)
+    if turnover_slice is not None:
+        center = np.asarray(constraints.turnover_center)
+        radius = float(constraints.turnover_limit)
+        positive = np.zeros((n, total_dimension), dtype=float)
+        positive[:, :n] = np.eye(n)
+        positive[:, turnover_slice] = -np.eye(n)
+        negative = np.zeros((n, total_dimension), dtype=float)
+        negative[:, :n] = -np.eye(n)
+        negative[:, turnover_slice] = -np.eye(n)
+        total = np.zeros((1, total_dimension), dtype=float)
+        total[:, turnover_slice] = 1.0
+        inequality_rows.extend([positive, negative, total])
+        inequality_targets.extend([center, -center, np.array([radius])])
+    if gross_slice is not None:
+        radius = float(constraints.gross_exposure_limit)
+        positive = np.zeros((n, total_dimension), dtype=float)
+        positive[:, :n] = np.eye(n)
+        positive[:, gross_slice] = -np.eye(n)
+        negative = np.zeros((n, total_dimension), dtype=float)
+        negative[:, :n] = -np.eye(n)
+        negative[:, gross_slice] = -np.eye(n)
+        total = np.zeros((1, total_dimension), dtype=float)
+        total[:, gross_slice] = 1.0
+        inequality_rows.extend([positive, negative, total])
+        inequality_targets.extend([np.zeros(n), np.zeros(n), np.array([radius])])
+
+    scipy_constraints: list[dict[str, Any]] = []
+    if equality_rows:
+        lifted_a_eq = np.vstack(equality_rows)
+        lifted_b_eq = np.concatenate(equality_targets)
+        scipy_constraints.append(
+            {
+                "type": "eq",
+                "fun": lambda y, a=lifted_a_eq, b=lifted_b_eq: a @ y - b,
+                "jac": lambda y, a=lifted_a_eq: a,
+            }
+        )
+    if inequality_rows:
+        lifted_a_ub = np.vstack(inequality_rows)
+        lifted_b_ub = np.concatenate(inequality_targets)
+        scipy_constraints.append(
+            {
+                "type": "ineq",
+                "fun": lambda y, a=lifted_a_ub, b=lifted_b_ub: b - a @ y,
+                "jac": lambda y, a=lifted_a_ub: -a,
+            }
+        )
+
+    lower = np.full(total_dimension, -np.inf, dtype=float)
+    upper = np.full(total_dimension, np.inf, dtype=float)
+    lower[:n] = np.asarray(constraints.lower_bounds)
+    upper[:n] = np.asarray(constraints.upper_bounds)
+    if turnover_slice is not None:
+        lower[turnover_slice] = 0.0
+    if gross_slice is not None:
+        lower[gross_slice] = 0.0
+    result = minimize(
+        lifted_value,
+        lifted_initial,
+        jac=lifted_gradient,
+        method="SLSQP",
+        bounds=Bounds(lower, upper),
+        constraints=scipy_constraints,
+        options={"maxiter": max_iterations, "ftol": tolerance, "disp": False},
+    )
+    holdings = np.asarray(result.x[:n], dtype=float)
+    return ReferenceResult(
+        holdings=holdings,
+        objective=problem.value(holdings),
+        success=bool(result.success),
+        message=str(result.message),
+        raw_result=result,
+    )
diff --git a/baseline/experiments/portfolio_pgd/src/portfolio_pgd/solver.py b/baseline/experiments/portfolio_pgd/src/portfolio_pgd/solver.py
new file mode 100644
index 00000000..f2177bea
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/src/portfolio_pgd/solver.py
@@ -0,0 +1,243 @@
+"""Projected-gradient solver with majorization backtracking and diagnostics."""
+
+from __future__ import annotations
+
+from dataclasses import dataclass, field
+from time import perf_counter
+from typing import Callable
+
+import numpy as np
+from numpy.typing import ArrayLike, NDArray
+
+from .constraints import ConstraintSet
+from .problem import PortfolioProblem
+
+FloatArray = NDArray[np.float64]
+
+
+@dataclass(frozen=True)
+class PGDOptions:
+    max_iterations: int = 10_000
+    tolerance: float = 1.0e-8
+    step_size: float | None = None
+    use_backtracking: bool = True
+    backtracking_factor: float = 0.5
+    step_growth: float = 1.15
+    minimum_step: float = 1.0e-16
+    maximum_step: float = 1.0e6
+    max_backtracking_steps: int = 60
+    projection_tolerance: float = 1.0e-10
+    projection_max_cycles: int = 5_000
+    record_every: int = 1
+    progress_interval: int = 50
+
+    def __post_init__(self) -> None:
+        if self.max_iterations <= 0:
+            raise ValueError("max_iterations must be positive")
+        if self.tolerance <= 0.0:
+            raise ValueError("tolerance must be positive")
+        if self.step_size is not None and self.step_size <= 0.0:
+            raise ValueError("step_size must be positive")
+        if not 0.0 < self.backtracking_factor < 1.0:
+            raise ValueError("backtracking_factor must lie in (0, 1)")
+        if self.step_growth < 1.0:
+            raise ValueError("step_growth must be at least one")
+        if self.record_every <= 0:
+            raise ValueError("record_every must be positive")
+        if self.progress_interval <= 0:
+            raise ValueError("progress_interval must be positive")
+
+
+@dataclass(frozen=True)
+class ProgressState:
+    """Immutable progress record emitted by a running PGD solve."""
+
+    iteration: int
+    objective: float
+    utility: float
+    projected_gradient_norm: float
+    step_size: float
+    max_constraint_violation: float
+    elapsed_seconds: float
+
+
+@dataclass
+class SolverResult:
+    holdings: FloatArray
+    trades: FloatArray
+    objective: float
+    utility: float
+    converged: bool
+    status: str
+    iterations: int
+    projected_gradient_norm: float
+    max_constraint_violation: float
+    history: dict[str, list[float]] = field(default_factory=dict)
+
+
+def solve_pgd(
+    problem: PortfolioProblem,
+    constraints: ConstraintSet | None = None,
+    *,
+    initial_holdings: ArrayLike | None = None,
+    options: PGDOptions | None = None,
+    progress_callback: Callable[[ProgressState], None] | None = None,
+) -> SolverResult:
+    """Solve a smooth convex portfolio problem by projected gradient descent."""
+    started = perf_counter()
+    options = options or PGDOptions()
+    constraints = constraints or ConstraintSet(problem.dimension)
+    if constraints.dimension != problem.dimension:
+        raise ValueError("problem and constraints have different dimensions")
+    x0 = problem.previous_holdings if initial_holdings is None else np.asarray(initial_holdings, dtype=float)
+    x = constraints.project(
+        x0,
+        tolerance=options.projection_tolerance,
+        max_cycles=options.projection_max_cycles,
+    )
+    objective = problem.value(x)
+    if not np.isfinite(objective):
+        raise FloatingPointError("initial objective is not finite")
+
+    if options.step_size is None:
+        step = 1.0 / problem.local_lipschitz_bound(x)
+    else:
+        step = float(options.step_size)
+    step = float(np.clip(step, options.minimum_step, options.maximum_step))
+
+    history: dict[str, list[float]] = {
+        "iteration": [0.0],
+        "objective": [objective],
+        "utility": [-objective],
+        "projected_gradient_norm": [np.nan],
+        "step_size": [step],
+        "max_constraint_violation": [constraints.max_violation(x)],
+    }
+    converged = False
+    status = "maximum_iterations_reached"
+    projected_gradient_norm = np.inf
+    last_callback_iteration = -1
+
+    if progress_callback is not None:
+        progress_callback(
+            ProgressState(
+                iteration=0,
+                objective=float(objective),
+                utility=float(-objective),
+                projected_gradient_norm=float("nan"),
+                step_size=step,
+                max_constraint_violation=constraints.max_violation(x),
+                elapsed_seconds=perf_counter() - started,
+            )
+        )
+        last_callback_iteration = 0
+
+    for iteration in range(1, options.max_iterations + 1):
+        gradient = problem.gradient(x)
+        if np.any(~np.isfinite(gradient)):
+            raise FloatingPointError("objective gradient is not finite")
+
+        trial_step = (
+            min(step * options.step_growth, options.maximum_step)
+            if options.use_backtracking
+            else step
+        )
+        candidate = x
+        candidate_objective = objective
+        accepted = False
+        for _ in range(options.max_backtracking_steps):
+            candidate = constraints.project(
+                x - trial_step * gradient,
+                tolerance=options.projection_tolerance,
+                max_cycles=options.projection_max_cycles,
+            )
+            displacement = candidate - x
+            candidate_objective = problem.value(candidate)
+            majorizer = (
+                objective
+                + float(gradient @ displacement)
+                + 0.5 * float(displacement @ displacement) / trial_step
+            )
+            slack = 1.0e-13 * max(1.0, abs(objective), abs(candidate_objective))
+            if (
+                not options.use_backtracking
+                or candidate_objective <= majorizer + slack
+            ):
+                accepted = True
+                break
+            trial_step *= options.backtracking_factor
+            if trial_step < options.minimum_step:
+                break
+        if not accepted:
+            status = "line_search_failed"
+            break
+
+        projected_gradient_norm = float(np.linalg.norm(x - candidate) / trial_step)
+        x = candidate
+        objective = candidate_objective
+        step = trial_step
+        violation = constraints.max_violation(x)
+
+        if iteration % options.record_every == 0 or iteration == options.max_iterations:
+            history["iteration"].append(float(iteration))
+            history["objective"].append(float(objective))
+            history["utility"].append(float(-objective))
+            history["projected_gradient_norm"].append(projected_gradient_norm)
+            history["step_size"].append(step)
+            history["max_constraint_violation"].append(violation)
+
+        should_report = iteration % options.progress_interval == 0
+
+        scale = max(1.0, float(np.linalg.norm(x)))
+        if projected_gradient_norm <= options.tolerance * scale:
+            converged = True
+            status = "converged"
+            should_report = True
+            if history["iteration"][-1] != float(iteration):
+                history["iteration"].append(float(iteration))
+                history["objective"].append(float(objective))
+                history["utility"].append(float(-objective))
+                history["projected_gradient_norm"].append(projected_gradient_norm)
+                history["step_size"].append(step)
+                history["max_constraint_violation"].append(violation)
+        if progress_callback is not None and should_report:
+            progress_callback(
+                ProgressState(
+                    iteration=iteration,
+                    objective=float(objective),
+                    utility=float(-objective),
+                    projected_gradient_norm=projected_gradient_norm,
+                    step_size=step,
+                    max_constraint_violation=violation,
+                    elapsed_seconds=perf_counter() - started,
+                )
+            )
+            last_callback_iteration = iteration
+        if converged:
+            break
+
+    iterations = iteration
+    if progress_callback is not None and last_callback_iteration != iterations:
+        progress_callback(
+            ProgressState(
+                iteration=iterations,
+                objective=float(objective),
+                utility=float(-objective),
+                projected_gradient_norm=projected_gradient_norm,
+                step_size=step,
+                max_constraint_violation=constraints.max_violation(x),
+                elapsed_seconds=perf_counter() - started,
+            )
+        )
+    return SolverResult(
+        holdings=np.asarray(x, dtype=float),
+        trades=np.asarray(x - problem.previous_holdings, dtype=float),
+        objective=float(objective),
+        utility=float(-objective),
+        converged=converged,
+        status=status,
+        iterations=iterations,
+        projected_gradient_norm=projected_gradient_norm,
+        max_constraint_violation=constraints.max_violation(x),
+        history=history,
+    )
diff --git a/baseline/experiments/portfolio_pgd/src/portfolio_pgd/synthetic.py b/baseline/experiments/portfolio_pgd/src/portfolio_pgd/synthetic.py
new file mode 100644
index 00000000..8be9e812
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/src/portfolio_pgd/synthetic.py
@@ -0,0 +1,58 @@
+"""Deterministic synthetic data helpers for examples and tests."""
+
+from __future__ import annotations
+
+import numpy as np
+from numpy.typing import NDArray
+
+FloatArray = NDArray[np.float64]
+
+
+def factor_covariance(
+    n_assets: int,
+    n_factors: int = 4,
+    *,
+    seed: int = 7,
+    specific_risk: float = 0.08,
+) -> tuple[FloatArray, FloatArray]:
+    """Create a positive-definite covariance and its asset-factor loadings."""
+    rng = np.random.default_rng(seed)
+    loadings = rng.normal(scale=0.35, size=(n_assets, n_factors))
+    raw = rng.normal(size=(n_factors, n_factors))
+    factor_cov = raw @ raw.T / n_factors + 0.2 * np.eye(n_factors)
+    idiosyncratic = specific_risk * (0.5 + rng.random(n_assets))
+    covariance = loadings @ factor_cov @ loadings.T + np.diag(idiosyncratic)
+    covariance = 0.5 * (covariance + covariance.T)
+    return covariance, loadings
+
+
+def capped_long_only_portfolio(n_assets: int, *, cap: float, seed: int = 11) -> FloatArray:
+    """Generate a fully invested long-only portfolio below a per-name cap."""
+    if cap * n_assets < 1.0:
+        raise ValueError("cap is too small for a fully invested portfolio")
+    rng = np.random.default_rng(seed)
+    weights = np.full(n_assets, 1.0 / n_assets)
+    perturbation = rng.normal(scale=0.15 / n_assets, size=n_assets)
+    perturbation -= np.mean(perturbation)
+    weights = np.clip(weights + perturbation, 0.0, cap)
+    for _ in range(100):
+        deficit = 1.0 - float(np.sum(weights))
+        if abs(deficit) < 1.0e-14:
+            break
+        if deficit > 0.0:
+            room = cap - weights
+            active = room > 1.0e-14
+            weights[active] += deficit * room[active] / np.sum(room[active])
+        else:
+            active = weights > 1.0e-14
+            weights[active] += deficit * weights[active] / np.sum(weights[active])
+        weights = np.clip(weights, 0.0, cap)
+    return weights / np.sum(weights)
+
+
+def sector_membership(n_assets: int, n_sectors: int) -> FloatArray:
+    """Return a sector-by-asset binary membership matrix."""
+    sectors = np.zeros((n_sectors, n_assets), dtype=float)
+    for asset in range(n_assets):
+        sectors[asset % n_sectors, asset] = 1.0
+    return sectors
diff --git a/baseline/experiments/portfolio_pgd/tests/__init__.py b/baseline/experiments/portfolio_pgd/tests/__init__.py
new file mode 100644
index 00000000..38bb211b
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/tests/__init__.py
@@ -0,0 +1 @@
+"""Test package."""
diff --git a/baseline/experiments/portfolio_pgd/tests/run_tests.py b/baseline/experiments/portfolio_pgd/tests/run_tests.py
new file mode 100644
index 00000000..b0c41208
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/tests/run_tests.py
@@ -0,0 +1,16 @@
+"""Dependency-free test entry point: python tests/run_tests.py"""
+
+from __future__ import annotations
+
+import pathlib
+import sys
+import unittest
+
+ROOT = pathlib.Path(__file__).resolve().parents[1]
+sys.path.insert(0, str(ROOT / "src"))
+
+
+if __name__ == "__main__":
+    suite = unittest.defaultTestLoader.discover(str(ROOT / "tests"), pattern="test_*.py")
+    result = unittest.TextTestRunner(verbosity=2).run(suite)
+    raise SystemExit(0 if result.wasSuccessful() else 1)
diff --git a/baseline/experiments/portfolio_pgd/tests/test_constraints.py b/baseline/experiments/portfolio_pgd/tests/test_constraints.py
new file mode 100644
index 00000000..bec7c88f
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/tests/test_constraints.py
@@ -0,0 +1,101 @@
+from __future__ import annotations
+
+import unittest
+
+import numpy as np
+
+from portfolio_pgd import (
+    ConstraintSet,
+    PGDOptions,
+    PortfolioProblem,
+    PowerLawCost,
+    capped_long_only_portfolio,
+    factor_covariance,
+    sector_membership,
+    solve_pgd,
+    solve_scipy_slsqp,
+)
+
+
+def realistic_case(n: int = 16):
+    rng = np.random.default_rng(301)
+    covariance, loadings = factor_covariance(n, 3, seed=302, specific_risk=0.20)
+    previous = capped_long_only_portfolio(n, cap=0.12, seed=303)
+    sectors = sector_membership(n, 4)
+    previous_sector = sectors @ previous
+    lower_sector = np.maximum(0.10, previous_sector - 0.05)
+    upper_sector = np.minimum(0.40, previous_sector + 0.05)
+    a_ub = np.vstack([sectors, -sectors])
+    b_ub = np.concatenate([upper_sector, -lower_sector])
+    factor = loadings[:, 0] - np.mean(loadings[:, 0])
+    equality = np.vstack([np.ones(n), factor])
+    target = np.array([1.0, float(factor @ previous)])
+    constraints = ConstraintSet(
+        n,
+        equality_matrix=equality,
+        equality_target=target,
+        inequality_matrix=a_ub,
+        inequality_upper=b_ub,
+        lower_bounds=0.0,
+        upper_bounds=0.14,
+        turnover_limit=0.24,
+        turnover_center=previous,
+    )
+    problem = PortfolioProblem(
+        alpha=rng.normal(scale=0.04, size=n),
+        covariance=covariance,
+        previous_holdings=previous,
+        risk_aversion=1.6,
+        quadratic_cost_matrix=0.2 + rng.random(n),
+        quadratic_cost_aversion=0.2,
+        nonlinear_cost=PowerLawCost(eta=0.008, p=1.5, epsilon=1.0e-3),
+    )
+    return problem, constraints, sectors
+
+
+class ConstraintProjectionTests(unittest.TestCase):
+    def test_dykstra_projection_satisfies_realistic_intersection(self) -> None:
+        problem, constraints, _ = realistic_case()
+        rng = np.random.default_rng(304)
+        projected = constraints.project(rng.normal(scale=0.5, size=problem.dimension))
+        self.assertLess(constraints.max_violation(projected), 2.0e-8)
+        reprojection = constraints.project(projected)
+        np.testing.assert_allclose(projected, reprojection, atol=2.0e-8, rtol=0.0)
+
+    def test_realistic_constraints_pgd_matches_standard_solver(self) -> None:
+        problem, constraints, _ = realistic_case()
+        result = solve_pgd(
+            problem,
+            constraints,
+            options=PGDOptions(
+                max_iterations=20_000,
+                tolerance=1.0e-7,
+                projection_tolerance=2.0e-10,
+            ),
+        )
+        standard = solve_scipy_slsqp(problem, constraints, tolerance=1.0e-10)
+        self.assertTrue(result.converged, result.status)
+        self.assertTrue(standard.success, standard.message)
+        self.assertLess(result.max_constraint_violation, 2.0e-7)
+        self.assertLess(constraints.max_violation(standard.holdings), 2.0e-6)
+        self.assertLess(abs(result.objective - standard.objective), 1.0e-5)
+
+    def test_long_short_gross_exposure_projection(self) -> None:
+        n = 12
+        rng = np.random.default_rng(305)
+        beta = rng.normal(size=n)
+        constraints = ConstraintSet(
+            n,
+            equality_matrix=np.vstack([np.ones(n), beta]),
+            equality_target=np.array([0.0, 0.0]),
+            lower_bounds=-0.20,
+            upper_bounds=0.20,
+            gross_exposure_limit=1.0,
+        )
+        projected = constraints.project(rng.normal(size=n))
+        self.assertLess(constraints.max_violation(projected), 2.0e-8)
+        self.assertLessEqual(np.sum(np.abs(projected)), 1.0 + 2.0e-8)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/baseline/experiments/portfolio_pgd/tests/test_nonlinear.py b/baseline/experiments/portfolio_pgd/tests/test_nonlinear.py
new file mode 100644
index 00000000..2caa1b1d
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/tests/test_nonlinear.py
@@ -0,0 +1,86 @@
+from __future__ import annotations
+
+import unittest
+
+import numpy as np
+
+from portfolio_pgd import (
+    ConstraintSet,
+    PGDOptions,
+    PortfolioProblem,
+    PowerLawCost,
+    SmoothAbsoluteCost,
+    capped_long_only_portfolio,
+    factor_covariance,
+    solve_pgd,
+    solve_scipy_slsqp,
+)
+
+
+class NonlinearCostTests(unittest.TestCase):
+    def test_power_law_gradient_matches_finite_difference(self) -> None:
+        rng = np.random.default_rng(200)
+        trades = rng.normal(scale=0.1, size=15)
+        cost = PowerLawCost(eta=np.linspace(0.01, 0.04, 15), p=1.5, epsilon=2.0e-4)
+        analytic = cost.gradient(trades)
+        epsilon = 1.0e-7
+        finite_difference = np.empty_like(trades)
+        for index in range(trades.size):
+            perturbation = np.zeros_like(trades)
+            perturbation[index] = epsilon
+            finite_difference[index] = (
+                cost.value(trades + perturbation) - cost.value(trades - perturbation)
+            ) / (2.0 * epsilon)
+        np.testing.assert_allclose(analytic, finite_difference, rtol=2.0e-6, atol=2.0e-8)
+
+    def test_smoothed_absolute_gradient_matches_finite_difference(self) -> None:
+        rng = np.random.default_rng(201)
+        trades = rng.normal(scale=0.05, size=10)
+        cost = SmoothAbsoluteCost(rate=0.003, epsilon=1.0e-3)
+        direction = rng.normal(size=10)
+        direction /= np.linalg.norm(direction)
+        epsilon = 1.0e-7
+        finite_difference = (
+            cost.value(trades + epsilon * direction)
+            - cost.value(trades - epsilon * direction)
+        ) / (2.0 * epsilon)
+        self.assertAlmostEqual(finite_difference, float(cost.gradient(trades) @ direction), places=8)
+
+    def test_nonlinear_pgd_matches_scipy(self) -> None:
+        n = 14
+        covariance, _ = factor_covariance(n, 3, seed=210, specific_risk=0.18)
+        rng = np.random.default_rng(211)
+        previous = capped_long_only_portfolio(n, cap=0.14, seed=212)
+        problem = PortfolioProblem(
+            alpha=rng.normal(scale=0.035, size=n),
+            covariance=covariance,
+            previous_holdings=previous,
+            risk_aversion=1.8,
+            quadratic_cost_matrix=0.15 + rng.random(n),
+            quadratic_cost_aversion=0.25,
+            nonlinear_cost=PowerLawCost(
+                eta=0.01 + 0.015 * rng.random(n), p=1.5, epsilon=1.0e-3
+            ),
+        )
+        constraints = ConstraintSet(
+            n,
+            equality_matrix=np.ones((1, n)),
+            equality_target=np.array([1.0]),
+            lower_bounds=0.0,
+            upper_bounds=0.16,
+        )
+        result = solve_pgd(
+            problem,
+            constraints,
+            options=PGDOptions(max_iterations=20_000, tolerance=5.0e-8),
+        )
+        standard = solve_scipy_slsqp(problem, constraints)
+        self.assertTrue(result.converged, result.status)
+        self.assertTrue(standard.success, standard.message)
+        self.assertLess(result.max_constraint_violation, 2.0e-8)
+        self.assertLess(abs(result.objective - standard.objective), 5.0e-8)
+        self.assertLess(np.linalg.norm(result.holdings - standard.holdings), 3.0e-4)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/baseline/experiments/portfolio_pgd/tests/test_quadratic.py b/baseline/experiments/portfolio_pgd/tests/test_quadratic.py
new file mode 100644
index 00000000..8c0ef9fc
--- /dev/null
+++ b/baseline/experiments/portfolio_pgd/tests/test_quadratic.py
@@ -0,0 +1,95 @@
+from __future__ import annotations
+
+import unittest
+
+import numpy as np
+
+from portfolio_pgd import (
+    ConstraintSet,
+    PGDOptions,
+    PortfolioProblem,
+    factor_covariance,
+    solve_pgd,
+    solve_quadratic_kkt,
+    solve_scipy_slsqp,
+)
+
+
+class QuadraticPortfolioTests(unittest.TestCase):
+    def setUp(self) -> None:
+        n = 18
+        covariance, loadings = factor_covariance(n, 3, seed=101, specific_risk=0.12)
+        rng = np.random.default_rng(102)
+        alpha = rng.normal(scale=0.025, size=n)
+        previous = rng.normal(scale=0.01, size=n)
+        q_diag = 0.5 + rng.random(n)
+        self.problem = PortfolioProblem(
+            alpha=alpha,
+            covariance=covariance,
+            previous_holdings=previous,
+            risk_aversion=2.5,
+            quadratic_cost_matrix=q_diag,
+            quadratic_cost_aversion=0.8,
+        )
+        factor_direction = loadings[:, 0]
+        factor_direction = factor_direction - np.mean(factor_direction)
+        self.constraints = ConstraintSet(
+            n,
+            equality_matrix=np.vstack([np.ones(n), factor_direction]),
+            equality_target=np.array([1.0, 0.0]),
+        )
+
+    def test_pgd_matches_exact_kkt_holdings_and_objective(self) -> None:
+        exact = solve_quadratic_kkt(self.problem, self.constraints)
+        result = solve_pgd(
+            self.problem,
+            self.constraints,
+            options=PGDOptions(max_iterations=20_000, tolerance=2.0e-9),
+        )
+        self.assertTrue(result.converged, result.status)
+        self.assertLess(result.max_constraint_violation, 1.0e-8)
+        self.assertLess(np.linalg.norm(result.holdings - exact.holdings), 2.0e-7)
+        self.assertAlmostEqual(result.objective, exact.objective, places=10)
+
+    def test_pgd_and_kkt_match_scipy_slsqp(self) -> None:
+        exact = solve_quadratic_kkt(self.problem, self.constraints)
+        standard = solve_scipy_slsqp(self.problem, self.constraints)
+        self.assertTrue(standard.success, standard.message)
+        self.assertLess(np.linalg.norm(standard.holdings - exact.holdings), 2.0e-6)
+        self.assertAlmostEqual(standard.objective, exact.objective, places=9)
+
+    def test_objective_gradient(self) -> None:
+        rng = np.random.default_rng(18)
+        h = rng.normal(size=self.problem.dimension)
+        direction = rng.normal(size=self.problem.dimension)
+        direction /= np.linalg.norm(direction)
+        epsilon = 1.0e-6
+        finite_difference = (
+            self.problem.value(h + epsilon * direction)
+            - self.problem.value(h - epsilon * direction)
+        ) / (2.0 * epsilon)
+        analytic = float(self.problem.gradient(h) @ direction)
+        self.assertAlmostEqual(finite_difference, analytic, places=7)
+
+    def test_progress_callback_reports_initial_and_final_state(self) -> None:
+        states = []
+        result = solve_pgd(
+            self.problem,
+            self.constraints,
+            options=PGDOptions(
+                max_iterations=20_000,
+                tolerance=2.0e-9,
+                progress_interval=25,
+            ),
+            progress_callback=states.append,
+        )
+        self.assertTrue(result.converged)
+        self.assertGreaterEqual(len(states), 2)
+        self.assertEqual(states[0].iteration, 0)
+        self.assertEqual(states[-1].iteration, result.iterations)
+        self.assertAlmostEqual(states[-1].objective, result.objective, places=14)
+        self.assertLess(states[-1].max_constraint_violation, 1.0e-8)
+
+
+if __name__ == "__main__":
+    unittest.main()