Remove cvxpy to compute the first point (#296)

* introduce the cla function

* obsolete empty cell?

* remove cvxpy

* remove the cvxpy construct

* notebook updated

* linting

Author: tschm

book/marimo/cla.py

@@ -3,14 +3,14 @@
 # dependencies = [
 #     "marimo==0.14.13",
 #     "numpy==2.3.0",
-#     "cvxcla==1.3.1"
+#     "cvxcla==1.3.2"
 # ]
 # ///
 """Little demo for the Critical Line Algorithm."""
 
 import marimo
 
-__generated_with = "0.13.15"
+__generated_with = "0.14.13"
 app = marimo.App()
 
 with app.setup:
@@ -35,13 +35,23 @@ def _():
 @app.cell
 def _():
     slider = mo.ui.slider(4, 100, step=1, value=10, label="Size of the problem")
+    # display the slider
     slider
     return (slider,)
 
 
-@app.cell(hide_code=True)
-def _(slider):
-    n = slider.value
+@app.function(hide_code=True)
+def cla(n):
+    """Compute using the Critical Line Algorithm (CLA) an efficient frontier.
+
+    Args:
+        n (int): The dimension size of the mean vector, lower and upper bounds
+            arrays, and covariance matrix used in the computation.
+
+    Returns:
+        numpy.ndarray: The efficient frontier generated by the CLA based on the
+            provided parameters.
+    """
     mean = np.random.randn(n)
     lower_bounds = np.zeros(n)
     upper_bounds = np.ones(n)
@@ -57,18 +67,14 @@ def _(slider):
         a=np.ones((1, len(mean))),
         b=np.ones(1),
     ).frontier
-    return (f1,)
-
-
-@app.cell
-def _(f1):
-    f1.interpolate(2).plot(volatility=True, markers=True)
-    f1.plot()
-    return
+    return f1
 
 
 @app.cell
-def _():
+def _(slider):
+    frontier = cla(slider.value)
+    frontier.interpolate(2).plot(volatility=True, markers=True)
+    frontier.plot()
     return
 
 

experiments/init_algo.py

@@ -15,7 +15,7 @@
 import numpy as np
 from loguru import logger
 
-from src.cvxcla.first import init_algo, init_algo_lp
+from cvxcla.first import init_algo
 
 if __name__ == "__main__":
     # Define a large-scale portfolio problem with 10,000 assets
@@ -32,9 +32,3 @@
     # Print the indices of free variables in the solution
     print("Free variables from init_algo:")
     print(np.where(tp.free)[0])
-
-    # Compute the first turning point using the linear programming implementation
-    tp = init_algo_lp(mean=mean, lower_bounds=np.zeros_like(upper_bound), upper_bounds=upper_bound)
-    # Print the indices of free variables in the solution
-    print("Free variables from init_algo_lp:")
-    print(np.where(tp.free)[0])

pyproject.toml

@@ -5,10 +5,8 @@ description = "Critical line algorithm for the efficient frontier"
 readme = "README.md"
 requires-python = ">=3.10"
 dependencies = [
-    "cvxpy-base>=1.5.1",
     "numpy>=2.0.0",
     "scipy>=1.14.1",
-    "clarabel>=0.9.0",
     "plotly>=6.0.1",
     "kaleido==1.0.0"
 ]
@@ -42,12 +40,9 @@ build-backend = "hatchling.build"
 packages = ["src/cvxcla"]
 
 [tool.deptry.per_rule_ignores]
-#DEP001 = ["cvxpy"]
-DEP002 = ["clarabel", "kaleido"]
+DEP002 = ["kaleido"]
 
 [tool.deptry]
 # see https://deptry.com/usage/#pep-621-dev-dependency-groups
 pep621_dev_dependency_groups = ["dev"]
 
-[tool.deptry.package_module_name_map]
-cvxpy-base = ["cvxpy"]

src/cvxcla/first.py

@@ -20,7 +20,6 @@
 
 from __future__ import annotations
 
-import cvxpy as cp
 import numpy as np
 from numpy.typing import NDArray
 
@@ -75,95 +74,6 @@ def init_algo(
     return TurningPoint(free=free, weights=weights)
 
 
-def init_algo_lp(
-    mean: NDArray[np.float64],
-    lower_bounds: NDArray[np.float64],
-    upper_bounds: NDArray[np.float64],
-    a_eq: NDArray[np.float64] | None = None,
-    b_eq: NDArray[np.float64] | None = None,
-    solver=cp.CLARABEL,
-    **kwargs,
-    # A_ub: NDArray[np.float64] | None = None,
-    # b_ub: NDArray[np.float64] | None = None,
-) -> TurningPoint:
-    """Compute the first turning point using linear programming.
-
-    This function formulates the problem of finding the first turning point as a linear
-    programming problem and solves it using a convex optimization solver. The objective
-    is to maximize the expected return subject to the constraints that the weights sum
-    to 1 and are within their bounds.
-
-    Args:
-        mean: Vector of expected returns for each asset.
-        lower_bounds: Vector of lower bounds for asset weights.
-        upper_bounds: Vector of upper bounds for asset weights.
-        a_eq: Matrix for additional linear equality constraints (Ax = b).
-            If None, only the fully invested constraint (sum(weights) = 1) is used.
-        b_eq: Vector for additional linear equality constraints (Ax = b).
-            If None, only the fully invested constraint (sum(weights) = 1) is used.
-        solver: The CVXPY solver to use for the optimization.
-        **kwargs: Additional keyword arguments to pass to the solver.
-
-    Returns:
-        A TurningPoint object representing the first point on the efficient frontier.
-
-    Raises:
-        ValueError: If the problem is infeasible or if lower bounds exceed upper bounds.
-
-    """
-    if a_eq is None:
-        a_eq = np.atleast_2d(np.ones_like(mean))
-
-    if b_eq is None:
-        b_eq = np.array([1.0])
-
-    # if A_ub is None:
-    #    A_ub = np.atleast_2d(np.zeros_like(mean))
-
-    # if b_ub is None:
-    #    b_ub = np.array([0.0])
-
-    w = cp.Variable(mean.shape[0], "weights")
-
-    objective = cp.Maximize(mean.T @ w)
-    constraints = [
-        a_eq @ w == b_eq,
-        # A_ub @ w <= b_ub,
-        lower_bounds <= w,
-        w <= upper_bounds,
-        cp.sum(w) == 1.0,
-    ]
-
-    problem = cp.Problem(objective, constraints)
-    problem.solve(solver=solver, **kwargs)
-    # check status of problem is optimal
-    if problem.status != cp.OPTIMAL:
-        raise ValueError("Could not construct a fully invested portfolio")
-
-    # assert problem.status == cp.OPTIMAL
-    # print(problem.status)
-    # print(status)
-
-    w = w.value
-
-    # compute the distance from the closest bound
-    # distance = np.min(
-    #    np.array([np.abs(w - lower_bounds), np.abs(upper_bounds - w)]), axis=0
-    # )
-
-    # which element has the largest distance to any bound?
-    # Even if all assets are at their bounds,
-    # we get a (somewhat random) free asset.
-    # index = np.argmax(distance)
-
-    # free = np.full_like(mean, False, dtype=np.bool_)
-    # free[index] = True
-
-    free = _free(w, lower_bounds, upper_bounds)
-
-    return TurningPoint(free=free, weights=w)
-
-
 def _free(
     w: NDArray[np.float64], lower_bounds: NDArray[np.float64], upper_bounds: NDArray[np.float64]
 ) -> NDArray[np.bool_]:

src/tests/test_cla.py

@@ -13,7 +13,7 @@
 from pandas import DataFrame
 
 from cvxcla import CLA
-from cvxcla.first import init_algo_lp
+from cvxcla.first import init_algo
 
 
 def test_solver(input_data: Any, results: Any) -> None:
@@ -71,7 +71,7 @@ def test_example(example: DataFrame, example_solution: DataFrame) -> None:
     upper_bounds = 0.5 * np.ones(ns)
 
     # covariance = example.cov().values
-    tp = init_algo_lp(mean=means.values, lower_bounds=lower_bounds, upper_bounds=upper_bounds)
+    tp = init_algo(mean=means.values, lower_bounds=lower_bounds, upper_bounds=upper_bounds)
     assert np.allclose(tp.weights, np.array([0.1, 0.5, 0.4]), atol=1e-9)
     assert np.allclose(tp.free, np.array([False, False, True]))
 

src/tests/test_first.py

@@ -10,7 +10,7 @@
 import numpy as np
 import pytest
 
-from cvxcla.first import init_algo, init_algo_lp
+from cvxcla.first import init_algo
 from cvxcla.types import TurningPoint
 
 
@@ -29,12 +29,6 @@ def test_init_algo(n: int) -> None:
     lower_bounds = np.zeros(n)
     upper_bounds = np.random.rand(n)
 
-    first = init_algo_lp(mean=mean, lower_bounds=lower_bounds, upper_bounds=upper_bounds)
-
-    assert np.sum(first.free) == 1
-    assert np.sum(first.weights) == pytest.approx(1.0)
-    assert isinstance(first, TurningPoint)
-
     first = init_algo(mean=mean, lower_bounds=lower_bounds, upper_bounds=upper_bounds)
 
     assert np.sum(first.free) == 1
@@ -51,11 +45,6 @@ def test_small() -> None:
     mean = np.array([1.0, 2.0, 3.0])
     lower_bound = np.array([0.0, 0.0, 0.0])
     upper_bound = np.array([0.4, 0.4, 0.4])
-    tp = init_algo_lp(mean=mean, lower_bounds=lower_bound, upper_bounds=upper_bound)
-
-    assert np.allclose(tp.weights, [0.2, 0.4, 0.4])
-    assert tp.lamb == np.inf
-    assert np.allclose(tp.free, [True, False, False])
 
     tp = init_algo(mean=mean, lower_bounds=lower_bound, upper_bounds=upper_bound)
 
@@ -74,9 +63,6 @@ def test_no_fully_invested_portfolio() -> None:
     lower_bounds = np.array([0.0, 0.0, 0.0])
     upper_bounds = np.array([0.2, 0.2, 0.2])
 
-    with pytest.raises(ValueError, match="Could not construct a fully invested portfolio"):
-        init_algo_lp(mean=mean, lower_bounds=lower_bounds, upper_bounds=upper_bounds)
-
     with pytest.raises(ValueError, match="Could not construct a fully invested portfolio"):
         init_algo(mean=mean, lower_bounds=lower_bounds, upper_bounds=upper_bounds)
 
@@ -91,8 +77,5 @@ def test_lb_ub_mixed() -> None:
     lower_bounds = np.ones(3)
     mean = np.ones(3)
 
-    with pytest.raises(ValueError):
-        init_algo_lp(mean=mean, lower_bounds=lower_bounds, upper_bounds=upper_bounds)
-
     with pytest.raises(ValueError):
         init_algo(mean=mean, lower_bounds=lower_bounds, upper_bounds=upper_bounds)