quanttoolbox.optim

optim.proximal / optim.projection

Python alternatives

Keep — pyproximal covers some overlapping norm/constraint proximal operators, but is oriented toward signal-processing/inverse-problems use cases, not portfolio constraints (turnover, combined Dykstra-projected linear+box systems). No strong off-the-shelf equivalent for this exact operator set.

quanttoolbox.optim.proximal

Proximal operators for L1/L2/Linf norms, box/equality/inequality constraints, and combined ("Dykstra alternating projection") constraint sets.

Ported from QuantToolBox/optim/{proximal_L1,proximal_L2,proximal_Linfinity, proximal_max,proximal_bounds,proximal_equality,proximal_inequality, proximal_linear_constraints,proximal_turnover,soft_thresholding}.m

Translation notes:

• proximal_L1 and soft_thresholding (2-argument form) are algebraically identical in the original; only soft_thresholding is kept here (as the more common name), with proximal_l1 as an alias.
• proximal_bounds' original had two branches: a trivial closed-form clip (Proximal_Algorithm == 1) and a redundant quadprog call for the same box projection (Proximal_Algorithm == 2, which solves the exact same problem -- projecting onto a box has a closed-form solution, so a QP solver is never actually needed). Only the closed-form clip is ported.
• proximal_equality's Dykstra-loop branch was commented out in the original in favor of the closed-form pinv solution (which is exact for a single equality-constraint projection); that closed form is what's ported here.
• proximal_inequality/proximal_linear_constraints/ proximal_turnover's Dykstra alternating-projection loops (combining multiple constraint sets) are preserved, since combined constraint projection generally has no closed form.
• MATLAB's global Proximal_MaxIters is replaced by quanttoolbox.config.ProximalConfig.

proximal_bounds(v, lb, ub)

Projection onto a box [lb, ub] (closed-form clip).

Original: optim/proximal_bounds.m

Source code in src/quanttoolbox/optim/proximal.py

def proximal_bounds(v: np.ndarray, lb: np.ndarray | float, ub: np.ndarray | float) -> np.ndarray:
    """Projection onto a box [lb, ub] (closed-form clip).

    Original: optim/proximal_bounds.m
    """
    v = np.asarray(v, dtype=float)
    return np.clip(v, lb, ub)

proximal_equality(v, a_eq, b_eq)

Projection onto the affine subspace {x : A_eq @ x == b_eq}.

Original: optim/proximal_equality.m

Returns (x, retcode) with retcode always 0 (kept for API parity with proximal_inequality/proximal_linear_constraints, which can fail).

Source code in src/quanttoolbox/optim/proximal.py

def proximal_equality(v: np.ndarray, a_eq: np.ndarray, b_eq: np.ndarray) -> tuple[np.ndarray, int]:
    """Projection onto the affine subspace {x : A_eq @ x == b_eq}.

    Original: optim/proximal_equality.m

    Returns (x, retcode) with retcode always 0 (kept for API parity with
    proximal_inequality/proximal_linear_constraints, which can fail).
    """
    v = np.asarray(v, dtype=float)
    a_eq = np.asarray(a_eq, dtype=float)
    b_eq = np.asarray(b_eq, dtype=float).flatten()
    x = v - np.linalg.pinv(a_eq) @ (a_eq @ v - b_eq)
    return x, 0

proximal_inequality(v, c_ineq, d_ineq, config=None)

Projection onto a polyhedron {x : C_ineq @ x <= D_ineq} via Dykstra's alternating projection algorithm (cycling through each row/half-space).

Original: optim/proximal_inequality.m

Returns (x, retcode); retcode is -1 if max_iters was reached without convergence, 0 otherwise.

Source code in src/quanttoolbox/optim/proximal.py

def proximal_inequality(
    v: np.ndarray, c_ineq: np.ndarray, d_ineq: np.ndarray, config: ProximalConfig | None = None
) -> tuple[np.ndarray, int]:
    """Projection onto a polyhedron {x : C_ineq @ x <= D_ineq} via Dykstra's
    alternating projection algorithm (cycling through each row/half-space).

    Original: optim/proximal_inequality.m

    Returns (x, retcode); retcode is -1 if max_iters was reached without
    convergence, 0 otherwise.
    """
    if config is None:
        config = ProximalConfig()

    v = np.asarray(v, dtype=float)
    c_ineq = np.asarray(c_ineq, dtype=float)
    d_ineq = np.asarray(d_ineq, dtype=float).flatten()
    n = v.shape[0]
    n_ineq = c_ineq.shape[0]

    x = v.copy()
    z = np.zeros((n, n_ineq))
    retcode = 0

    for _it in range(config.max_iters):
        x_old = x.copy()
        for i in range(n_ineq):
            vi = x + z[:, i]
            c = c_ineq[i, :]
            d = d_ineq[i]
            u = np.dot(c, vi) - d
            x_new = vi - u * (u >= 0) * c / np.dot(c, c)
            z[:, i] = x + z[:, i] - x_new
            x = x_new
        if np.allclose(x_old, x):
            break
    else:
        retcode = -1

    return x, retcode

proximal_l2(v, lambda_)

Proximal operator of the (scaled) L2 norm: shrinks v toward the origin by at most lambda_ in Euclidean length.

Original: optim/proximal_L2.m

Source code in src/quanttoolbox/optim/proximal.py

def proximal_l2(v: np.ndarray, lambda_: float) -> np.ndarray:
    """Proximal operator of the (scaled) L2 norm: shrinks v toward the
    origin by at most lambda_ in Euclidean length.

    Original: optim/proximal_L2.m
    """
    v = np.asarray(v, dtype=float)
    norm_v = float(np.linalg.norm(v))
    return (1 - lambda_ / max(lambda_, norm_v)) * v

proximal_linear_constraints(v, a_eq=None, b_eq=None, c_ineq=None, d_ineq=None, lb=None, ub=None, config=None)

Projection onto the intersection of an affine subspace, a polyhedron, and a box, via Dykstra's alternating projection algorithm. Pass None for any constraint set to omit it.

Original: optim/proximal_linear_constraints.m

Note: like the original, the exact-equality convergence check (x1 == x4 in MATLAB, np.allclose here) can occasionally exit early during a temporary plateau in the iterate sequence, before reaching the true intersection point -- this is most likely for inputs sitting near a "corner" where the constraint sets meet tangentially. If a result looks suspicious, verify constraint satisfaction directly and re-run with a larger max_iters and a perturbed starting point if needed.

Source code in src/quanttoolbox/optim/proximal.py

def proximal_linear_constraints(
    v: np.ndarray,
    a_eq: np.ndarray | None = None,
    b_eq: np.ndarray | None = None,
    c_ineq: np.ndarray | None = None,
    d_ineq: np.ndarray | None = None,
    lb: np.ndarray | float | None = None,
    ub: np.ndarray | float | None = None,
    config: ProximalConfig | None = None,
) -> tuple[np.ndarray, int]:
    """Projection onto the intersection of an affine subspace, a polyhedron,
    and a box, via Dykstra's alternating projection algorithm. Pass None
    for any constraint set to omit it.

    Original: optim/proximal_linear_constraints.m

    Note: like the original, the exact-equality convergence check
    (``x1 == x4`` in MATLAB, ``np.allclose`` here) can occasionally exit
    early during a temporary plateau in the iterate sequence, before
    reaching the true intersection point -- this is most likely for
    inputs sitting near a "corner" where the constraint sets meet
    tangentially. If a result looks suspicious, verify constraint
    satisfaction directly and re-run with a larger ``max_iters`` and a
    perturbed starting point if needed.
    """
    if config is None:
        config = ProximalConfig()

    v = np.asarray(v, dtype=float)
    n = v.shape[0]

    x1 = v.copy()
    delta1 = np.zeros(n)
    delta2 = np.zeros(n)
    delta3 = np.zeros(n)
    retcode = 0

    for _it in range(config.max_iters):
        if a_eq is not None:
            x2, _ = proximal_equality(x1 + delta1, a_eq, b_eq)
            delta1 = x1 + delta1 - x2
        else:
            x2 = x1

        if c_ineq is not None:
            x3, _ = proximal_inequality(x2 + delta2, c_ineq, d_ineq, config)
            delta2 = x2 + delta2 - x3
        else:
            x3 = x2

        if lb is not None:
            x4 = proximal_bounds(x3 + delta3, lb, ub)
            delta3 = x3 + delta3 - x4
        else:
            x4 = x3

        if np.allclose(x1, x4):
            x1 = x4
            break
        x1 = x4
    else:
        retcode = -1

    return x1, retcode

proximal_linfinity(v, lambda_)

Proximal operator of lambda * ||v||_infinity.

Original: optim/proximal_Linfinity.m

Source code in src/quanttoolbox/optim/proximal.py

def proximal_linfinity(v: np.ndarray, lambda_: float) -> np.ndarray:
    """Proximal operator of lambda * ||v||_infinity.

    Original: optim/proximal_Linfinity.m
    """
    v = np.asarray(v, dtype=float)
    return np.sign(v) * proximal_max(np.abs(v), lambda_)

proximal_max(v, lambda_)

Proximal operator of lambda * max(v): caps the largest entries of v so that the total amount "shaved off" sums to lambda_ (used inside proximal_Linfinity / projection_L1's simplex-style water-filling step).

Original: optim/proximal_max.m

Source code in src/quanttoolbox/optim/proximal.py

def proximal_max(v: np.ndarray, lambda_: float) -> np.ndarray:
    """Proximal operator of lambda * max(v): caps the largest entries of v
    so that the total amount "shaved off" sums to lambda_ (used inside
    proximal_Linfinity / projection_L1's simplex-style water-filling step).

    Original: optim/proximal_max.m
    """
    v = np.asarray(v, dtype=float)
    n = v.shape[0]
    sorted_v = np.sort(v)[::-1]
    c = (np.cumsum(sorted_v) - lambda_) / np.arange(1, n + 1)
    mask = sorted_v > c
    rk = np.where(mask)[0][-1]  # last index where condition holds
    s = c[rk]
    return np.minimum(v, s)

proximal_turnover(v, a_eq, b_eq, c_ineq, d_ineq, lb, ub, x0, tau, config=None)

Projection onto the intersection of an affine subspace, a polyhedron, a box, and a turnover constraint (||x - x0||_1 <= tau), via Dykstra's alternating projection algorithm.

Original: optim/proximal_turnover.m

Source code in src/quanttoolbox/optim/proximal.py

def proximal_turnover(
    v: np.ndarray,
    a_eq: np.ndarray | None,
    b_eq: np.ndarray | None,
    c_ineq: np.ndarray | None,
    d_ineq: np.ndarray | None,
    lb: np.ndarray | float | None,
    ub: np.ndarray | float | None,
    x0: np.ndarray,
    tau: float | None,
    config: ProximalConfig | None = None,
) -> tuple[np.ndarray, int]:
    """Projection onto the intersection of an affine subspace, a polyhedron,
    a box, and a turnover constraint (||x - x0||_1 <= tau), via Dykstra's
    alternating projection algorithm.

    Original: optim/proximal_turnover.m
    """
    if config is None:
        config = ProximalConfig()

    from quanttoolbox.optim.projection import projection_l1

    v = np.asarray(v, dtype=float)
    x0 = np.asarray(x0, dtype=float)
    n = v.shape[0]

    x1 = v.copy()
    delta1 = np.zeros(n)
    delta2 = np.zeros(n)
    delta3 = np.zeros(n)
    delta4 = np.zeros(n)
    retcode = 0

    for _it in range(config.max_iters):
        if a_eq is not None:
            x2, _ = proximal_equality(x1 + delta1, a_eq, b_eq)
            delta1 = x1 + delta1 - x2
        else:
            x2 = x1

        if c_ineq is not None:
            x3, _ = proximal_inequality(x2 + delta2, c_ineq, d_ineq, config)
            delta2 = x2 + delta2 - x3
        else:
            x3 = x2

        if lb is not None:
            x4 = proximal_bounds(x3 + delta3, lb, ub)
            delta3 = x3 + delta3 - x4
        else:
            x4 = x3

        if tau is not None:
            x5 = projection_l1((x4 - x0) + delta4, tau) + x0
            delta4 = x4 + delta4 - x5
        else:
            x5 = x4

        if np.allclose(x1, x5):
            x1 = x5
            break
        x1 = x5
    else:
        retcode = -1

    return x1, retcode

soft_thresholding(v, lambda_minus, lambda_plus=None)

Soft-thresholding / proximal operator of the L1 norm.

One-argument-lambda form: symmetric soft threshold, sign(v) * max(|v|-lambda, 0). Two-argument form: asymmetric threshold with separate positive/negative shrinkage.

Original: optim/soft_thresholding.m (also optim/proximal_L1.m, identical)

Source code in src/quanttoolbox/optim/proximal.py

def soft_thresholding(
    v: np.ndarray, lambda_minus: float, lambda_plus: float | None = None
) -> np.ndarray:
    """Soft-thresholding / proximal operator of the L1 norm.

    One-argument-lambda form: symmetric soft threshold, sign(v) * max(|v|-lambda, 0).
    Two-argument form: asymmetric threshold with separate positive/negative shrinkage.

    Original: optim/soft_thresholding.m (also optim/proximal_L1.m, identical)
    """
    v = np.asarray(v, dtype=float)
    if lambda_plus is None:
        lam = lambda_minus
        return np.sign(v) * np.maximum(np.abs(v) - lam, 0.0)
    positive_part = np.maximum(v - lambda_plus, 0.0)
    negative_part = np.maximum(-(v + lambda_minus), 0.0)
    return positive_part - negative_part

quanttoolbox.optim.projection

Projection operators onto L1/L2/Linf norm balls and a box-intersect-L2-ball set.

Ported from QuantToolBox/optim/{projection_L1,projection_L2, projection_Linfinity,projection_box_L2}.m

Translation notes:

• projection_L2 is literally v - proximal_l2(v, lambda) in the original (projection = identity minus the proximal/shrinkage step); reproduced as-is here rather than re-derived, to stay faithful.
• projection_box_L2 (project onto a box intersected with an L2 ball centered at c) uses Dykstra's alternating projection algorithm, same pattern as proximal_linear_constraints.

projection_box_l2(v, x_minus, x_plus, c, lambda_, config=None)

Projection onto the intersection of a box [x_minus, x_plus] and an L2 ball of radius lambda_ centered at c, via Dykstra's alternating projection algorithm.

Original: optim/projection_box_L2.m

Source code in src/quanttoolbox/optim/projection.py

def projection_box_l2(
    v: np.ndarray,
    x_minus: np.ndarray | float,
    x_plus: np.ndarray | float,
    c: np.ndarray,
    lambda_: float,
    config: ProximalConfig | None = None,
) -> tuple[np.ndarray, int]:
    """Projection onto the intersection of a box [x_minus, x_plus] and an L2
    ball of radius lambda_ centered at c, via Dykstra's alternating
    projection algorithm.

    Original: optim/projection_box_L2.m
    """
    if config is None:
        config = ProximalConfig()

    v = np.asarray(v, dtype=float)
    c = np.asarray(c, dtype=float)
    n = v.shape[0]

    x2 = v.copy()
    delta1 = np.zeros(n)
    delta2 = np.zeros(n)
    retcode = 1

    for _it in range(config.max_iters):
        x1 = proximal_bounds(x2 + delta1, x_minus, x_plus)
        delta1 = x2 + delta1 - x1

        x2_new = c + projection_l2((x1 + delta2 - c), lambda_)
        delta2 = x1 + delta2 - x2_new

        if np.allclose(x1, x2_new):
            x2 = x2_new
            break
        x2 = x2_new
    else:
        retcode = -1

    return x1, retcode

projection_l1(v, radius, method=1)

Euclidean projection of v onto the L1 ball of the given radius.

method=1 (default): exact sorted-cumsum water-filling algorithm. method=2: equivalent computation via proximal_max.

Original: optim/projection_L1.m

Source code in src/quanttoolbox/optim/projection.py

def projection_l1(v: np.ndarray, radius: float, method: int = 1) -> np.ndarray:
    """Euclidean projection of v onto the L1 ball of the given radius.

    method=1 (default): exact sorted-cumsum water-filling algorithm.
    method=2: equivalent computation via proximal_max.

    Original: optim/projection_L1.m
    """
    v = np.asarray(v, dtype=float)
    if np.sum(np.abs(v)) <= radius:
        return v.copy()

    if method == 1:
        n = v.shape[0]
        abs_v = np.sort(np.abs(v))[::-1]
        mu = (np.cumsum(abs_v) - radius) / np.arange(1, n + 1)
        rk = np.where(abs_v > mu)[0][-1]
        lambda_star = mu[rk]
        from quanttoolbox.optim.proximal import soft_thresholding

        return soft_thresholding(v, lambda_star)
    else:
        return v - proximal_max(np.abs(v), radius) * np.sign(v)

projection_l2(v, lambda_)

Euclidean projection removing at most lambda_ of v's L2 length (identity minus the L2 proximal/shrinkage step).

Original: optim/projection_L2.m

Source code in src/quanttoolbox/optim/projection.py

def projection_l2(v: np.ndarray, lambda_: float) -> np.ndarray:
    """Euclidean projection removing at most lambda_ of v's L2 length
    (identity minus the L2 proximal/shrinkage step).

    Original: optim/projection_L2.m
    """
    v = np.asarray(v, dtype=float)
    return v - proximal_l2(v, lambda_)

projection_linfinity(v, radius)

Euclidean projection of v onto the L-infinity ball of the given radius (simple elementwise clip to [-radius, radius]).

Original: optim/projection_Linfinity.m

Source code in src/quanttoolbox/optim/projection.py

def projection_linfinity(v: np.ndarray, radius: float) -> np.ndarray:
    """Euclidean projection of v onto the L-infinity ball of the given radius
    (simple elementwise clip to [-radius, radius]).

    Original: optim/projection_Linfinity.m
    """
    v = np.asarray(v, dtype=float)
    return np.clip(v, -radius, radius)

Examples

Proximal projection under bounds, (in)equalities, and combined constraints — optim/proximal1.py

"""Translated from Examples/optim/proximal1.m -- proximal projections onto
boxes, single/multiple inequality constraints, single/multiple equality
constraints, and combined linear-constraint sets, on a fixed random
starting point.

The original compares `Proximal_Algorithm = 1` (closed-form, where one
exists) against `Proximal_Algorithm = 2` (an alternate/QP-based
computation of the *same* projection) side by side. During porting, the
redundant `Proximal_Algorithm == 2` branches for `proximal_bounds` and
`proximal_equality` were dropped (see `proximal.py`'s module docstring --
they solve the exact same problem the closed form already solves exactly,
so keeping both added no information), and `proximal_inequality`/
`proximal_linear_constraints` only ever had one (Dykstra) implementation.
So there is no second variant left to compare here -- each constraint
type is projected onto once below, rather than twice.

The original draws x from MATLAB's unseeded `rand`; a fixed seed
(`np.random.default_rng(0)`) is substituted here."""

import numpy as np

from quanttoolbox.optim.proximal import (
    proximal_bounds,
    proximal_equality,
    proximal_inequality,
    proximal_linear_constraints,
)

rng = np.random.default_rng(0)
n = 10
x = rng.random(n)

print("x:", np.round(x, 4))

# Lower & upper bounds
lb, ub = np.zeros(n), 0.5 * np.ones(n)
x1 = proximal_bounds(x, lb, ub)
print("\nBounds [0, 0.5]:", np.round(x1, 4))

# One inequality constraint
c_ineq = np.zeros((1, n))
c_ineq[0, :4] = 0.25
d_ineq = np.array([0.5])
x1, rc = proximal_inequality(x, c_ineq, d_ineq)
print("\nInequality (1 constraint):", np.round(x1, 4), "retcode:", rc)

# Two inequality constraints
c_ineq = np.zeros((2, n))
d_ineq = np.zeros(2)
c_ineq[0, :4] = 0.25
c_ineq[1, 3:6] = -np.array([1, 2, 2])
d_ineq[0] = 0.5
d_ineq[1] = -2.0
x1, rc = proximal_inequality(x, c_ineq, d_ineq)
print("\nInequality (2 constraints):", np.round(x1, 4), "retcode:", rc)

# Four inequality constraints
c_ineq = np.zeros((4, n))
d_ineq = np.zeros(4)
c_ineq[0, :4] = 0.25
c_ineq[1, 3:6] = -np.array([1, 2, 2])
c_ineq[2, 9] = 1
c_ineq[3, 9] = -1
d_ineq[0] = 0.5
d_ineq[1] = -2.0
d_ineq[2] = 0.1
d_ineq[3] = -0.1
x1, rc = proximal_inequality(x, c_ineq, d_ineq)
print("\nInequality (4 constraints):", np.round(x1, 4), "retcode:", rc)

# One equality constraint
a_eq = np.ones((1, n))
b_eq = np.array([4.0])
x1, rc = proximal_equality(x, a_eq, b_eq)
print("\nEquality (1 constraint):", np.round(x1, 4), "retcode:", rc)

# Two equality constraints
a_eq = np.zeros((2, n))
b_eq = np.zeros(2)
a_eq[0, :] = 1.0
a_eq[1, [0, 1]] = [1, -1]
b_eq[0] = 4.0
x1, rc = proximal_equality(x, a_eq, b_eq)
print("\nEquality (2 constraints):", np.round(x1, 4), "retcode:", rc)

# Linear constraints: 1 equality + 1 inequality + bounds
a_eq = np.ones((1, n))
b_eq = np.array([4.0])
c_ineq = np.zeros((1, n))
d_ineq = np.array([0.5])
c_ineq[0, :4] = 0.25
lb, ub = np.zeros(n), 0.5 * np.ones(n)
x1, rc = proximal_linear_constraints(
    x, a_eq=a_eq, b_eq=b_eq, c_ineq=c_ineq, d_ineq=d_ineq, lb=lb, ub=ub
)
print("\nLinear constraints (1 eq, 1 ineq, bounds):", np.round(x1, 4), "retcode:", rc)

# Linear constraints: 2 equality + 2 inequality + bounds
a_eq = np.zeros((2, n))
b_eq = np.zeros(2)
a_eq[0, :] = 1.0
a_eq[1, [0, 1]] = [1, -1]
b_eq[0] = 4.0
c_ineq = np.zeros((2, n))
d_ineq = np.zeros(2)
c_ineq[0, :4] = 0.25
c_ineq[1, 3:6] = -np.array([1, 2, 2])
d_ineq[0] = 0.5
d_ineq[1] = -2.0
x1, rc = proximal_linear_constraints(
    x, a_eq=a_eq, b_eq=b_eq, c_ineq=c_ineq, d_ineq=d_ineq, lb=lb, ub=ub
)
print("\nLinear constraints (2 eq, 2 ineq, bounds):", np.round(x1, 4), "retcode:", rc)

Proximal-L1 vs. two L1-ball-projection turnover algorithms — optim/prox_turnover1.py

"""Translated from Examples/optim/prox_turnover1.m -- compares the L1
proximal (soft-thresholding) operator against the L1-ball projection (both
its default sorted-cumsum algorithm and its `proximal_max`-based
alternative) at two different starting points and lambda values.

The original draws `v` from MATLAB's unseeded `rand`; a fixed seed
(`np.random.default_rng(0)`) is substituted here for reproducibility, same
convention used elsewhere in this port (see matrix1.py, building_blocks.md)."""

import numpy as np

from quanttoolbox.optim.projection import projection_l1
from quanttoolbox.optim.proximal import soft_thresholding

rng = np.random.default_rng(0)
v = 5 * (rng.random(10) - 0.5)

lambda_ = 1.20
x0 = np.zeros(10)

x1 = soft_thresholding(v - x0, lambda_) + x0
x2 = projection_l1(v - x0, lambda_) + x0
x3 = projection_l1(v - x0, lambda_, method=0) + x0

x = np.column_stack([v, x1, x2, x3])
print("lambda =", lambda_, " x0 = 0")
print("columns: v, proximal_L1, projection_L1 (method 1), projection_L1 (method 0)")
print(x)
print("column sums of |x - x0|:", np.sum(np.abs(x - x0[:, None]), axis=0))

lambda_ = 2.00
x0 = np.ones(10)
x0 = x0 / np.sum(x0)

x1 = soft_thresholding(v - x0, lambda_) + x0
x2 = projection_l1(v - x0, lambda_) + x0
x3 = projection_l1(v - x0, lambda_, method=0) + x0

x = np.column_stack([v, x1, x2, x3])
print("\nlambda =", lambda_, " x0 = equal-weight (sums to 1)")
print("columns: v, proximal_L1, projection_L1 (method 1), projection_L1 (method 0)")
print(x)
print("column sums of |x - x0|:", np.sum(np.abs(x - x0[:, None]), axis=0))

optim.quadprog (solve_qp)

Python alternatives

Hybrid: for hot loops calling solve_qp many times (e.g. risk-budgeting's inner ADMM loop), calling qpsolvers.solve_qp directly — bypassing cvxpy's DSL-parsing overhead — would likely be faster. Worth profiling if performance matters; keep the cvxpy-based version for its more expressive constraint composition (ridge/lasso penalties, arbitrary constraints).

quanttoolbox.optim.quadprog

Quadratic programming: a general QP solver plus ridge/lasso/turnover penalized variants, and the closed-form QP-on-a-hyperplane solution.

Ported from QuantToolBox/optim/{quadprog_bc_ccd,quadprog_lasso, quadprog_ridge,quadprog_turnover,quadprog_mixed_norm, quadprog_mixed2_norm,qp_hyperplane}.m

Translation notes -- this is the single biggest architectural simplification in the whole port:

MATLAB's Optimization Toolbox quadprog cannot express an L1 penalty term directly, so the original toolbox works around this with a "variable-splitting" trick: introduce two extra non-negative variable blocks per L1 term (x = x+ - x-), turning each lasso/mixed-norm/turnover problem into a bigger, purely-quadratic QP that quadprog can solve. This produces four near-duplicate ~100-line MATLAB functions (quadprog_lasso/ridge/mixed_norm/mixed2_norm), each hand-building the block matrices for a different combination of penalties.

cvxpy supports L1/L2 norm terms and turnover (L1-ball) constraints natively in its objective/constraint DSL, so none of that variable-splitting machinery is needed here. All four MATLAB functions above collapse into one function, solve_qp, parameterized by which penalty terms and constraints are supplied:

solve_qp(q, r)                                    # plain QP  (replaces quadprog_ridge with no penalty)
solve_qp(q, r, ridge_penalty=(gamma, x_target))    # replaces quadprog_ridge
solve_qp(q, r, lasso_penalty=(gamma, x_target))    # replaces quadprog_lasso
solve_qp(q, r, ridge_penalty=..., lasso_penalty=...)  # replaces quadprog_mixed_norm
solve_qp(q, r, turnover=(x0, tau))                 # replaces the turnover-constrained
                                                    # branch of quadprog_turnover

(quadprog_mixed2_norm's "two separate lasso/ridge penalties toward two different targets" case is representable by simply adding the two ridge penalties together algebraically, since ridge penalties toward different targets combine into a single quadratic form -- not separately exposed as a distinct code path here.)

The objective throughout is 0.5 * x'Qx - R'x (MATLAB's convention, note the minus sign on the linear term), matching every one of the original functions' calling convention.

quadprog_bc_ccd (box-constrained QP via cyclical coordinate descent) is preserved as a separate, dependency-free NumPy implementation since it's a genuinely different (iterative, matrix-free-friendly) algorithm worth keeping available alongside the cvxpy path.

qp_hyperplane (closed-form QP with a single equality constraint) has an exact closed-form Lagrangian solution and needs no solver at all.

qp_bc_ccd(q, r, x_minus=None, x_plus=None, x_init=None, n_iters=100)

Box-constrained QP (min 0.5*x'Qx - R'x s.t. x_minus <= x <= x_plus) via cyclical coordinate descent -- a lightweight, dependency-free alternative to solve_qp for this specific (box-only) case.

Original: optim/quadprog_bc_ccd.m

Source code in src/quanttoolbox/optim/quadprog.py

def qp_bc_ccd(
    q: np.ndarray,
    r: np.ndarray,
    x_minus: np.ndarray | float | None = None,
    x_plus: np.ndarray | float | None = None,
    x_init: np.ndarray | None = None,
    n_iters: int = 100,
) -> tuple[np.ndarray, np.ndarray]:
    """Box-constrained QP (min 0.5*x'Qx - R'x s.t. x_minus <= x <= x_plus)
    via cyclical coordinate descent -- a lightweight, dependency-free
    alternative to solve_qp for this specific (box-only) case.

    Original: optim/quadprog_bc_ccd.m
    """
    q = np.asarray(q, dtype=float)
    r = np.asarray(r, dtype=float).flatten()
    n = q.shape[0]

    x = (
        np.asarray(x_init, dtype=float).copy()
        if x_init is not None and np.asarray(x_init).shape[0] == n
        else np.linalg.solve(q, r)
    )

    truncate = x_minus is not None or x_plus is not None
    if x_minus is not None:
        x_minus = np.full(n, x_minus) if np.isscalar(x_minus) else np.asarray(x_minus, dtype=float)
    if x_plus is not None:
        x_plus = np.full(n, x_plus) if np.isscalar(x_plus) else np.asarray(x_plus, dtype=float)

    x_path = np.zeros((n_iters, n))
    for it in range(n_iters):
        x_path[it] = x
        for i in range(n):
            new_x = x.copy()
            new_x[i] = 0.0
            new_x[i] = (r[i] - 0.5 * new_x @ q[:, i] - 0.5 * q[i, :] @ new_x) / q[i, i]
            if truncate:
                new_x[i] = min(max(x_minus[i], new_x[i]), x_plus[i])
            x = new_x

    return x, x_path

qp_hyperplane(q, r, a, b)

Closed-form solution of min 0.5*x'Qx - R'x s.t. a'x = b (single equality constraint), via the Lagrangian stationary conditions.

Original: optim/qp_hyperplane.m

Returns (x, lagrange_multiplier).

Source code in src/quanttoolbox/optim/quadprog.py

def qp_hyperplane(
    q: np.ndarray, r: np.ndarray, a: np.ndarray, b: float
) -> tuple[np.ndarray, float]:
    """Closed-form solution of min 0.5*x'Qx - R'x s.t. a'x = b (single
    equality constraint), via the Lagrangian stationary conditions.

    Original: optim/qp_hyperplane.m

    Returns (x, lagrange_multiplier).
    """
    q = np.asarray(q, dtype=float)
    r = np.asarray(r, dtype=float).flatten()
    a = np.asarray(a, dtype=float).flatten()

    inv_q = np.linalg.inv(q)
    lagrange = (b - a @ inv_q @ r) / (a @ inv_q @ a)
    x = inv_q @ (r + lagrange * a)
    return x, lagrange

solve_qp(q, r, a_eq=None, b_eq=None, c_ineq=None, d_ineq=None, lb=None, ub=None, ridge_penalty=None, lasso_penalty=None, turnover=None, default_budget_constraint=False)

Solve min 0.5*x'Qx - R'x subject to the given constraints and optional ridge/lasso penalty terms and turnover budget.

Original: optim/{quadprog_lasso,quadprog_ridge,quadprog_turnover, quadprog_mixed_norm,quadprog_mixed2_norm}.m (consolidated -- see module docstring)

Parameters:

Name	Type	Description	Default
q	QP objective 0.5*x'Qx - R'x.		required
r	QP objective 0.5*x'Qx - R'x.		required
a_eq	equality constraints A_eq @ x == B_eq.		None
b_eq	equality constraints A_eq @ x == B_eq.		None
c_ineq	inequality constraints C_ineq @ x <= D_ineq.		None
d_ineq	inequality constraints C_ineq @ x <= D_ineq.		None
lb	box bounds.		None
ub	box bounds.		None
ridge_penalty	(gamma, x_target) adds gamma * ||x - x_target||_2^2 to	the objective. gamma may be a scalar, a vector (diagonal), or a full matrix.	None
lasso_penalty	(gamma, x_target) adds gamma' * |x - x_target| to the	objective (elementwise L1, gamma may be scalar or vector).	None
turnover	(x0, tau) constrains ||x - x0||_1 <= tau.		None
default_budget_constraint	if True and a_eq/b_eq are not given, adds	the default sum(x) == 1 budget constraint (matches the original's isscalar(A) && A==0 sentinel convention).	False

Source code in src/quanttoolbox/optim/quadprog.py

def solve_qp(
    q: np.ndarray,
    r: np.ndarray,
    a_eq: np.ndarray | None = None,
    b_eq: np.ndarray | None = None,
    c_ineq: np.ndarray | None = None,
    d_ineq: np.ndarray | None = None,
    lb: np.ndarray | float | None = None,
    ub: np.ndarray | float | None = None,
    ridge_penalty: tuple[np.ndarray | float, np.ndarray] | None = None,
    lasso_penalty: tuple[np.ndarray | float, np.ndarray] | None = None,
    turnover: tuple[np.ndarray, float] | None = None,
    default_budget_constraint: bool = False,
) -> np.ndarray:
    """Solve min 0.5*x'Qx - R'x subject to the given constraints and
    optional ridge/lasso penalty terms and turnover budget.

    Original: optim/{quadprog_lasso,quadprog_ridge,quadprog_turnover,
    quadprog_mixed_norm,quadprog_mixed2_norm}.m (consolidated -- see
    module docstring)

    Parameters
    ----------
    q, r : QP objective 0.5*x'Qx - R'x.
    a_eq, b_eq : equality constraints A_eq @ x == B_eq.
    c_ineq, d_ineq : inequality constraints C_ineq @ x <= D_ineq.
    lb, ub : box bounds.
    ridge_penalty : (gamma, x_target) adds gamma * ||x - x_target||_2^2 to
        the objective. gamma may be a scalar, a vector (diagonal), or a
        full matrix.
    lasso_penalty : (gamma, x_target) adds gamma' * |x - x_target| to the
        objective (elementwise L1, gamma may be scalar or vector).
    turnover : (x0, tau) constrains ||x - x0||_1 <= tau.
    default_budget_constraint : if True and a_eq/b_eq are not given, adds
        the default sum(x) == 1 budget constraint (matches the original's
        ``isscalar(A) && A==0`` sentinel convention).
    """
    q = np.asarray(q, dtype=float)
    r = np.asarray(r, dtype=float).flatten()
    n = q.shape[0]

    x = cp.Variable(n)
    objective = 0.5 * cp.quad_form(x, cp.psd_wrap(q)) - r @ x

    if ridge_penalty is not None:
        gamma, x_target = ridge_penalty
        gamma = np.asarray(gamma, dtype=float)
        x_target = np.asarray(x_target, dtype=float).flatten()
        if gamma.ndim == 0:
            objective = objective + gamma.item() * cp.sum_squares(x - x_target)
        elif gamma.ndim == 1:
            objective = objective + cp.sum(cp.multiply(gamma, cp.square(x - x_target)))
        else:
            objective = objective + cp.quad_form(x - x_target, cp.psd_wrap(gamma))

    if lasso_penalty is not None:
        gamma, x_target = lasso_penalty
        gamma = np.asarray(gamma, dtype=float)
        x_target = np.asarray(x_target, dtype=float).flatten()
        if gamma.ndim == 0:
            objective = objective + gamma.item() * cp.norm1(x - x_target)
        else:
            objective = objective + cp.sum(cp.multiply(gamma, cp.abs(x - x_target)))

    constraints = []
    if a_eq is not None:
        constraints.append(
            np.asarray(a_eq, dtype=float) @ x == np.asarray(b_eq, dtype=float).flatten()
        )
    elif default_budget_constraint:
        constraints.append(cp.sum(x) == 1)

    if c_ineq is not None:
        constraints.append(
            np.asarray(c_ineq, dtype=float) @ x <= np.asarray(d_ineq, dtype=float).flatten()
        )

    if lb is not None:
        constraints.append(x >= lb)
    if ub is not None:
        constraints.append(x <= ub)

    if turnover is not None:
        x0, tau = turnover
        constraints.append(cp.norm1(x - np.asarray(x0, dtype=float).flatten()) <= tau)

    problem = cp.Problem(cp.Minimize(objective), constraints)
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        problem.solve()

        if problem.status not in ("optimal", "optimal_inaccurate"):
            # the default solver (OSQP) can struggle with poorly-conditioned
            # or rank-deficient Q (e.g. block-structured QPs like
            # epsilon-insensitive SVM regression's stacked [[Q,-Q],[-Q,Q]]);
            # CLARABEL is a robust interior-point fallback for exactly this
            # kind of case.
            problem.solve(solver=cp.CLARABEL)

    if problem.status not in ("optimal", "optimal_inaccurate"):
        raise RuntimeError(f"solve_qp: solver did not converge (status={problem.status})")

    if x.value is None:
        raise RuntimeError("solve_qp: solver returned no solution")

    return x.value

Examples

Mean-variance optimization plus ridge/lasso-penalized portfolios — rpb/test_lasso1.py

"""Translated from Examples/rpb/test_lasso1.m -- Roncalli [2013],
"Introduction to Risk Parity and Budgeting", Example 1 (page 53):
compares an unconstrained gamma-problem MVO portfolio against the same
problem with a budget (sum-to-1) constraint, and ridge/lasso-penalized
variants (toward 0 and toward equal-weight) of the unconstrained
problem.

`compute_mvo_portfolio`/`quadprog_ridge`/`quadprog_lasso` map to
`mvo_portfolio`/`solve_qp(..., ridge_penalty=...)`/`solve_qp(...,
lasso_penalty=...)` as established in test_lasso3.py/test_lasso5.py; the
original's `0, 0` sentinel arguments for "no equality constraint" /
"no bounds" map to simply omitting those keyword arguments."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.optim.quadprog import solve_qp
from quanttoolbox.portfolio.mean_variance import mvo_portfolio
from quanttoolbox.stats.moments import corr_to_cov

mu = np.array([0.05, 0.06, 0.08, 0.06])
sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)
n = 4
x0 = np.full(n, 1 / n)

# Case gamma-problem
gamma_x = 0.5

x1 = mvo_portfolio(mu, cov_matrix, gamma=gamma_x).weights
x2 = mvo_portfolio(
    mu, cov_matrix, gamma=gamma_x, a_eq=np.ones((1, n)), b_eq=np.array([1.0])
).weights

lambda_ridge = 0.03
s_ridge = lambda_ridge * np.eye(n)
x3 = solve_qp(cov_matrix, gamma_x * mu, ridge_penalty=(s_ridge, np.zeros(n)))
x4 = solve_qp(cov_matrix, gamma_x * mu, ridge_penalty=(s_ridge, x0))

lambda_lasso = 0.03 / 2
s_lasso = lambda_lasso * np.ones(n)
x5 = solve_qp(cov_matrix, gamma_x * mu, lasso_penalty=(s_lasso, np.zeros(n)))
x6 = solve_qp(cov_matrix, gamma_x * mu, lasso_penalty=(s_lasso, x0))

results = 100 * np.column_stack([x1, x2, x3, x4, x5, x6])
print("            x1      x2      x3      x4      x5      x6")
print(np.round(results, 2))

Mixed ridge+lasso penalties toward two different target vectors — rpb/test_lasso5.py

"""Translated from Examples/rpb/test_lasso5.m -- ridge/lasso/mixed
portfolios with two different ridge/lasso target vectors (equal-weight
and a custom 20/20/30/30 target)."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.optim.quadprog import solve_qp
from quanttoolbox.portfolio.mean_variance import mvo_portfolio
from quanttoolbox.stats.moments import corr_to_cov

mu = np.array([0.05, 0.06, 0.08, 0.06])
sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)
n = 4
a_eq, b_eq = np.ones((1, n)), np.array([1.0])
gamma_x = 0.5

y1 = np.full(n, 1 / n)
y2 = np.array([0.20, 0.20, 0.30, 0.30])
S_ridge = np.diag(np.diag(cov_matrix))
lambda_lasso = 0.005

x1 = mvo_portfolio(mu, cov_matrix, gamma=gamma_x, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0).weights
x2 = solve_qp(
    cov_matrix, gamma_x * mu, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0, ridge_penalty=(S_ridge, y1)
)
x4 = solve_qp(
    cov_matrix, gamma_x * mu, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0, lasso_penalty=(lambda_lasso, y1)
)
x6 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, y1),
    lasso_penalty=(lambda_lasso, y1),
)
# mixed with DIFFERENT targets for ridge (toward y1) vs lasso (toward y2)
x8 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, y1),
    lasso_penalty=(lambda_lasso, y2),
)

for name, x in [
    ("MVO", x1),
    ("Ridge->y1", x2),
    ("Lasso->y1", x4),
    ("Mixed(ridge->y1,lasso->y1)", x6),
    ("Mixed(ridge->y1,lasso->y2)", x8),
]:
    print(f"{name}: {np.round(x, 4)}")

Ridge, lasso, and mixed-norm penalized portfolios — rpb/test_lasso3.py

"""Translated from Examples/rpb/test_lasso3.m (byte-identical to
test_lasso4.m) -- MVO/ridge/lasso/mixed-penalty portfolios via the
consolidated solve_qp (replacing the original's separate quadprog_ridge/
quadprog_lasso/quadprog_mixed calls -- see optim/quadprog.py docstring)."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.optim.quadprog import solve_qp
from quanttoolbox.portfolio.mean_variance import mvo_portfolio
from quanttoolbox.stats.moments import corr_to_cov

mu = np.array([0.05, 0.06, 0.08, 0.06])
sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)
n = 4
x0 = np.full(n, 1 / n)
gamma_x = 0.5
a_eq, b_eq = np.ones((1, n)), np.array([1.0])

# gamma-problem (plain MVO)
r1 = mvo_portfolio(mu, cov_matrix, gamma=gamma_x, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0)
x1 = r1.weights

# ridge-problem: penalize toward zero, scaled by each asset's own variance
S_ridge = np.diag(np.diag(cov_matrix))
x2 = solve_qp(
    cov_matrix, gamma_x * mu, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0, ridge_penalty=(S_ridge, x0)
)

# lasso-problem: L1 penalty toward equal weight
lambda_lasso = 0.005
x3 = solve_qp(
    cov_matrix, gamma_x * mu, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0, lasso_penalty=(lambda_lasso, x0)
)

# mixed-problem: both ridge and lasso, toward various targets
x4 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, np.zeros(n)),
    lasso_penalty=(lambda_lasso, np.zeros(n)),
)
x5 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, x0),
    lasso_penalty=(lambda_lasso, np.zeros(n)),
)
x6 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, np.zeros(n)),
    lasso_penalty=(lambda_lasso, x0),
)
x7 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, x0),
    lasso_penalty=(lambda_lasso, x0),
)

for name, x in [
    ("MVO", x1),
    ("Ridge", x2),
    ("Lasso", x3),
    ("Mixed(0,0)", x4),
    ("Mixed(EW,0)", x5),
    ("Mixed(0,EW)", x6),
    ("Mixed(EW,EW)", x7),
]:
    print(f"{name}: {np.round(x, 4)}")

optim.bisection

Python alternatives

Hybrid: scipy.optimize.brentq/bisect converge faster (superlinear) for scalar root-finding and are compiled. Keep our vectorized (array-broadcast, many-roots-at-once) version — scipy's root finders are scalar-only.

quanttoolbox.optim.bisection

Scalar (or elementwise-vectorized) bisection root-finding, and linear constraint explicit<->implicit (null-space) parametrization conversion.

Ported from QuantToolBox/optim/{bisection,bisection2,explicit2implicit, implicit2explicit}.m

Translation notes:

• Both bisection variants are vectorized in the original (operating on arrays a/b elementwise, not just scalars) -- preserved here via plain NumPy elementwise operations.
• bisection2 carries an auxiliary state variable z through the function evaluations (useful when fhandle also needs to return some side computation to warm-start the next evaluation); the Python signature keeps the same (y, z) = fhandle(x, z) calling convention.
• explicit2implicit/implicit2explicit convert between an explicit linear-constraint representation (C @ x = c) and an implicit null-space parametrization (x = R @ r for free parameter r) -- MATLAB's null(...) maps to scipy.linalg.null_space.
• MATLAB's global BISECTION_Tol is replaced by quanttoolbox.config.BisectionConfig.

bisection(fhandle, a, b, config=None)

Find the root of fhandle within bracket [a, b] via bisection (elementwise, if a/b are arrays -- each element is bracketed and solved independently).

Original: optim/bisection.m

Source code in src/quanttoolbox/optim/bisection.py

def bisection(
    fhandle: Callable[[np.ndarray], np.ndarray],
    a: np.ndarray | float,
    b: np.ndarray | float,
    config: BisectionConfig | None = None,
) -> np.ndarray | float:
    """Find the root of fhandle within bracket [a, b] via bisection
    (elementwise, if a/b are arrays -- each element is bracketed and
    solved independently).

    Original: optim/bisection.m
    """
    if config is None:
        config = BisectionConfig()

    a = np.asarray(a, dtype=float)
    b = np.asarray(b, dtype=float)
    scalar_input = a.ndim == 0

    ya = fhandle(a)
    yb = fhandle(b)

    if np.any(ya * yb > 0):
        result = np.full(a.shape, np.nan)
        return float(result) if scalar_input else result

    if np.all(ya == 0):
        return float(a) if scalar_input else a.copy()
    if np.all(yb == 0):
        return float(b) if scalar_input else b.copy()

    increasing = ya >= 0  # if ya < 0, function increases a->b; else decreases

    c = (a + b) / 2
    for _ in range(config.max_iters):
        if np.max(np.abs(a - b)) <= config.tol:
            break
        c = (a + b) / 2
        yc = fhandle(c)
        e = np.where(increasing, yc > 0, yc < 0)
        a = np.where(e, c, a)
        b = np.where(e, b, c)

    c = (a + b) / 2
    return float(c) if scalar_input else c

bisection2(fhandle, a, b, z0, config=None)

Bisection root-finding where fhandle also threads an auxiliary state z through each evaluation: (y, z) = fhandle(x, z).

Original: optim/bisection2.m

Source code in src/quanttoolbox/optim/bisection.py

def bisection2(
    fhandle: Callable[[np.ndarray, np.ndarray], tuple[np.ndarray, np.ndarray]],
    a: np.ndarray | float,
    b: np.ndarray | float,
    z0: np.ndarray,
    config: BisectionConfig | None = None,
) -> tuple[np.ndarray | float, np.ndarray]:
    """Bisection root-finding where fhandle also threads an auxiliary state
    z through each evaluation: (y, z) = fhandle(x, z).

    Original: optim/bisection2.m
    """
    if config is None:
        config = BisectionConfig()

    a = np.asarray(a, dtype=float)
    b = np.asarray(b, dtype=float)
    scalar_input = a.ndim == 0

    ya, za = fhandle(a, z0)
    yb, zb = fhandle(b, z0)
    zc = (za + zb) / 2

    if np.any(ya * yb > 0):
        result = np.full(a.shape, np.nan)
        return (float(result) if scalar_input else result), zc

    if np.all(ya == 0):
        return (float(a) if scalar_input else a.copy()), za
    if np.all(yb == 0):
        return (float(b) if scalar_input else b.copy()), zb

    increasing = ya >= 0

    c = (a + b) / 2
    for _ in range(config.max_iters):
        if np.max(np.abs(a - b)) <= config.tol:
            break
        c = (a + b) / 2
        yc, zc = fhandle(c, zc)
        e = np.where(increasing, yc > 0, yc < 0)
        a = np.where(e, c, a)
        b = np.where(e, b, c)

    c = (a + b) / 2
    return (float(c) if scalar_input else c), zc

explicit_to_implicit(cc, c)

Convert an explicit linear-constraint system (CC @ x = c) into an implicit null-space parametrization x = RR @ r + r0, returning (RR, r0).

Original: optim/explicit2implicit.m

Source code in src/quanttoolbox/optim/bisection.py

def explicit_to_implicit(cc: np.ndarray, c: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Convert an explicit linear-constraint system (CC @ x = c) into an
    implicit null-space parametrization x = RR @ r + r0, returning (RR, r0).

    Original: optim/explicit2implicit.m
    """
    cc = np.asarray(cc, dtype=float)
    c = np.asarray(c, dtype=float).flatten()

    if cc.shape[0] != c.shape[0] or cc.shape[0] >= cc.shape[1]:
        raise ValueError("explicit_to_implicit: wrong size of CC and c")

    rr = linalg.null_space(cc)
    mx = np.max(np.abs(rr), axis=0)
    mx = mx + (mx == 0)
    rr = rr / mx

    r = np.linalg.pinv(cc.T @ cc) @ cc.T @ c
    return rr, r

implicit_to_explicit(rr, r)

Convert an implicit null-space parametrization (x = RR @ r) back into an explicit linear-constraint system (CC @ x = c).

Original: optim/implicit2explicit.m

Source code in src/quanttoolbox/optim/bisection.py

def implicit_to_explicit(rr: np.ndarray, r: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Convert an implicit null-space parametrization (x = RR @ r) back into
    an explicit linear-constraint system (CC @ x = c).

    Original: optim/implicit2explicit.m
    """
    rr = np.asarray(rr, dtype=float)
    r = np.asarray(r, dtype=float).flatten()

    if rr.shape[0] != r.shape[0] or rr.shape[0] <= rr.shape[1]:
        raise ValueError("implicit_to_explicit: wrong size of RR and r")

    cc = linalg.null_space(rr.T).T
    mx = np.max(np.abs(cc), axis=1)
    mx = mx + (mx == 0)
    cc = cc / mx[:, None]

    c = cc @ r
    return cc, c

Examples

Explicit-to-implicit conversion for three simultaneous restrictions — optim/explicit3.py

"""Translated from Examples/optim/explicit3.m -- explicit-to-implicit
conversion for three simultaneous zero-restrictions on an 8-parameter
vector (as would arise e.g. from AR1_12=0, MA1_11=0, MA1_21=0 in a VAR
specification)."""

import numpy as np

from quanttoolbox.optim.bisection import explicit_to_implicit

CC = np.zeros((3, 8))
CC[0, 2] = 1  # beta[2] = 0  (AR1_12)
CC[1, 4] = 1  # beta[4] = 0  (MA1_11)
CC[2, 5] = 1  # beta[5] = 0  (MA1_21)
c = np.zeros(3)

RR, r = explicit_to_implicit(CC, c)

print("R:")
print(RR)
print("\nr:", r)

Explicit/implicit constraint round-trip, plus a design() demo — optim/explicit2.py

"""Translated from Examples/optim/explicit2.m -- explicit<->implicit
constraint conversion for a single equality constraint (x[0] = x[1]), plus
a `design` matrix built from a category-index vector."""

import numpy as np

from quanttoolbox.linalg.special_matrices import design
from quanttoolbox.optim.bisection import explicit_to_implicit, implicit_to_explicit

CC = np.array([[1.0, -1.0, 0, 0, 0, 0, 0, 0]])  # constraint: x[0] - x[1] = 0
c = np.array([0.0])

RR, r = explicit_to_implicit(CC, c)
CC2, c2 = implicit_to_explicit(RR, r)

print("R:")
print(RR)
print("\nr:", r)
print("\nC:")
print(CC2)
print("\nc:", c2)

w = np.concatenate([[1.0, 1.0], np.arange(2, 8, dtype=float)])  # seqa(2,1,6)
RR2 = design(w)
r2 = np.zeros(8)

print("\nR (design(w)):")
print(RR2)
print("\nr:", r2)