marimo-team/marimo-learn

1# /// script2# requires-python = ">=3.13"3# dependencies = [4#     "clarabel>=0.11.1",5#     "cvxpy-base>=1.8.2",6#     "marimo",7#     "matplotlib==3.10.8",8#     "numpy==2.4.3",9#     "wigglystuff==0.2.37",10# ]11# ///12 13import marimo14 15__generated_with = "0.18.4"16app = marimo.App()17 18 19@app.cell20def _():21    import marimo as mo22    return (mo,)23 24 25@app.cell(hide_code=True)26def _(mo):27    mo.md(r"""28    # Quadratic Program29 30    A quadratic program is an optimization problem with a quadratic objective and31    affine equality and inequality constraints. A common standard form is the32    following:33 34    \[35        \begin{array}{ll}36        \text{minimize}   & (1/2)x^TPx + q^Tx\\37        \text{subject to} & Gx \leq h \\38                          & Ax = b.39        \end{array}40    \]41 42    Here $P \in \mathcal{S}^{n}_+$, $q \in \mathcal{R}^n$, $G \in \mathcal{R}^{m \times n}$, $h \in \mathcal{R}^m$, $A \in \mathcal{R}^{p \times n}$, and $b \in \mathcal{R}^p$ are problem data and $x \in \mathcal{R}^{n}$ is the optimization variable. The inequality constraint $Gx \leq h$ is elementwise.43 44    **Why quadratic programming?** Quadratic programs are convex optimization problems that generalize both least-squares and linear programming.They can be solved efficiently and reliably, even in real-time.45 46    **An example from finance.** A simple example of a quadratic program arises in finance. Suppose we have $n$ different stocks, an estimate $r \in \mathcal{R}^n$ of the expected return on each stock, and an estimate $\Sigma \in \mathcal{S}^{n}_+$ of the covariance of the returns. Then we solve the optimization problem47 48    \[49        \begin{array}{ll}50        \text{minimize}   & (1/2)x^T\Sigma x - r^Tx\\51        \text{subject to} & x \geq 0 \\52                          & \mathbf{1}^Tx = 1,53        \end{array}54    \]55 56    to find a nonnegative portfolio allocation $x \in \mathcal{R}^n_+$ that optimally balances expected return and variance of return.57 58    When we solve a quadratic program, in addition to a solution $x^\star$, we obtain a dual solution $\lambda^\star$ corresponding to the inequality constraints. A positive entry $\lambda^\star_i$ indicates that the constraint $g_i^Tx \leq h_i$ holds with equality for $x^\star$ and suggests that changing $h_i$ would change the optimal value.59    """)60    return61 62 63@app.cell(hide_code=True)64def _(mo):65    mo.md(r"""66    ## Example67 68    In this example, we use CVXPY to construct and solve a quadratic program.69    """)70    return71 72 73@app.cell74def _():75    import cvxpy as cp76    import numpy as np77    return cp, np78 79 80@app.cell(hide_code=True)81def _(mo):82    mo.md("""83    First we generate synthetic data. In this problem, we don't include equality constraints, only inequality.84    """)85    return86 87 88@app.cell89def _(np):90    m = 491    n = 292 93    np.random.seed(1)94    q = np.random.randn(n)95    G = np.random.randn(m, n)96    h = G @ np.random.randn(n)97    return G, h, n, q98 99 100@app.cell(hide_code=True)101def _(mo, np):102    import wigglystuff103 104    P_widget = mo.ui.anywidget(105        wigglystuff.Matrix(np.array([[4.0, -1.4], [-1.4, 4]]), step=0.1)106    )107 108    mo.md(109        f"""110        The quadratic form $P$ is equal to the symmetrized version of this111        matrix:112 113        {P_widget.center()}114        """115    )116    return (P_widget,)117 118 119@app.cell120def _(P_widget, np):121    P = 0.5 * (np.array(P_widget.matrix) + np.array(P_widget.matrix).T)122    return (P,)123 124 125@app.cell(hide_code=True)126def _(mo):127    mo.md(r"""128    Next, we specify the problem. Notice that we use the `quad_form` function from CVXPY to create the quadratic form $x^TPx$.129    """)130    return131 132 133@app.cell134def _(G, P, cp, h, n, q):135    x = cp.Variable(n)136 137    problem = cp.Problem(138        cp.Minimize((1 / 2) * cp.quad_form(x, P) + q.T @ x),139        [G @ x <= h],140    )141    _ = problem.solve()142    return problem, x143 144 145@app.cell(hide_code=True)146def _(mo, problem, x):147    mo.md(148        f"""149        The optimal value is {problem.value:.04f}.150 151        A solution $x$ is {mo.as_html(list(x.value))}152        A dual solution is is {mo.as_html(list(problem.constraints[0].dual_value))}153        """154    )155    return156 157 158@app.cell159def _(G, P, h, plot_contours, q, x):160    plot_contours(P, G, h, q, x.value)161    return162 163 164@app.cell(hide_code=True)165def _(mo):166    mo.md(r"""167    In this plot, the gray shaded region is the feasible region (points satisfying the inequality), and the ellipses are level curves of the quadratic form.168 169    **🌊 Try it!** Try changing the entries of $P$ above with your mouse. How do the170    level curves and the optimal value of $x$ change? Can you explain what you see?171    """)172    return173 174 175@app.cell(hide_code=True)176def _(P, mo):177    mo.md(178        rf"""179        The above contour lines were generated with180 181        \[182        P= \begin{{bmatrix}}183        {P[0, 0]:.01f} & {P[0, 1]:.01f} \\184        {P[1, 0]:.01f} & {P[1, 1]:.01f} \\185        \end{{bmatrix}}186        \]187        """188    )189    return190 191 192@app.cell(hide_code=True)193def _(np):194    def plot_contours(P, G, h, q, x_star):195        import matplotlib.pyplot as plt196 197        # Create a grid of x and y values.198        x = np.linspace(-5, 5, 400)199        y = np.linspace(-5, 5, 400)200        X, Y = np.meshgrid(x, y)201 202        # Compute the quadratic form Q(x, y) = a*x^2 + 2*b*x*y + c*y^2.203        # Here, a = P[0,0], b = P[0,1] (and P[1,0]), c = P[1,1]204        Z = (205            0.5 * (P[0, 0] * X**2 + 2 * P[0, 1] * X * Y + P[1, 1] * Y**2)206            + q[0] * X207            + q[1] * Y208        )209 210        # --- Evaluate the constraints on the grid ---211        # We stack X and Y to get a list of (x,y) points.212        points = np.vstack([X.ravel(), Y.ravel()]).T213 214        # Start with all points feasible215        feasible = np.ones(points.shape[0], dtype=bool)216 217        # Apply the inequality constraints Gx <= h.218        # Each row of G and corresponding h defines a condition.219        for i in range(G.shape[0]):220            # For a given point x, the condition is: G[i,0]*x + G[i,1]*y <= h[i]221            feasible &= points.dot(G[i]) <= h[i] + 1e-8  # small tolerance222        # Reshape the boolean mask back to grid shape.223        feasible_grid = feasible.reshape(X.shape)224 225        # --- Plot the feasible region and contour lines---226        plt.figure(figsize=(8, 6))227 228        # Use contourf to fill the region where feasible_grid is True.229        # We define two levels, so that points that are True (feasible) get one230        # color.231        plt.contourf(232            X,233            Y,234            feasible_grid,235            levels=[-0.5, 0.5, 1.5],236            colors=["white", "gray"],237            alpha=0.5,238        )239 240        contours = plt.contour(X, Y, Z, levels=10, cmap="viridis")241        plt.clabel(contours, inline=True, fontsize=8)242        plt.title("Feasible region and level curves")243        plt.xlabel("$x_1$")244        plt.ylabel("$y_2$")245        # plt.colorbar(contours, label='Q(x, y)')246 247        ax = plt.gca()248        # Optionally, mark and label the point x_star.249        ax.plot(x_star[0], x_star[1], "ko", markersize=5)250        ax.text(251            x_star[0],252            x_star[1],253            r"$\mathbf{x}^\star$",254            color="black",255            fontsize=12,256            verticalalignment="bottom",257            horizontalalignment="right",258        )259        return plt.gca()260    return (plot_contours,)261 262 263if __name__ == "__main__":264    app.run()265