Skip to content

Step operators

Step operators are the building blocks of a first-order method, covering affine steps, ReLU, saturated linear, and soft thresholding.

Step

A base class for thresholding and activation steps. This class defines a common interface for step-like operators that transform inputs element-wise.

A step represents the relation: output = f(arg), where f is the step function (e.g., ReLU, soft threshold, saturated linear).

Attributes:

Name Type Description
arg LinearExpression

Right-hand side (input) linear expression.
output LinearExpression

Left-hand side (output) linear expression. Bounds on the output are stored in output.lb and output.ub.
relax_binary_vars bool

If True, binary decision variables associated with the operator are relaxed to take continuous values.
skip_obbt bool

If True, OBBT bound tightening is skipped for this step type. Override in subclasses that compute bounds analytically.
obbt_output bool

If true, OBBT is also run on the step's output leaves once the step's relaxation is in the OBBT LP. Enabled for Affine with non-identity output_mat, whose output bounds cannot otherwise see the tightened input box without M^{-1}.
Notes

Input vs output bounds: - step.lb, step.ub: bounds on the INPUT (arg) - step.output.lb, step.output.ub: bounds on the OUTPUT (output)

When update_bounds(lb, ub) is called with tightened input bounds, subclasses should update both self.lb/self.ub (input) and self.output.lb/self.output.ub (output, by applying the step function).

Source code in fomtools/steps/step.py 
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
class Step:
    """
    A base class for thresholding and activation steps. This class
    defines a common interface for step-like operators that transform
    inputs element-wise.

    A step represents the relation: output = f(arg), where f is the step
    function (e.g., ReLU, soft threshold, saturated linear).

    Attributes
    ----------
    arg : LinearExpression
        Right-hand side (input) linear expression.
    output : LinearExpression
        Left-hand side (output) linear expression. Bounds on the output
        are stored in ``output.lb`` and ``output.ub``.
    relax_binary_vars : bool
        If True, binary decision variables associated with the operator
        are relaxed to take continuous values.
    skip_obbt : bool
        If True, OBBT bound tightening is skipped for this step type.
        Override in subclasses that compute bounds analytically.
    obbt_output : bool
        If true, OBBT is also run on the step's output leaves once the step's
        relaxation is in the OBBT LP. Enabled for `Affine` with non-identity
        ``output_mat``, whose output bounds cannot otherwise see the tightened
        input box without ``M^{-1}``.

    Notes
    -----
    Input vs output bounds:
        - ``step.lb``, ``step.ub``: bounds on the INPUT (arg)
        - ``step.output.lb``, ``step.output.ub``: bounds on the OUTPUT (output)

    When ``update_bounds(lb, ub)`` is called with tightened input bounds,
    subclasses should update both ``self.lb``/``self.ub`` (input) and
    ``self.output.lb``/``self.output.ub`` (output, by applying the step function).
    """

    skip_obbt: bool = False
    obbt_output: bool = False

    def __init__(self, arg: LinearExpression, relax_binary_vars: bool = False):
        self.arg = arg
        self.relax_binary_vars = relax_binary_vars
        # Default: every output coordinate is "projected" (i.e., the step
        # transformation applies to all rows). Subclasses that support
        # partial projection (e.g. ``Relu`` with ``proj_ranges``) override
        # this in their own ``__init__``.
        self.proj_indices = set(range(arg.ndim))

    @cached_property
    def output(self) -> LinearExpression:
        """Left-hand side linear expression in the operator definition."""
        return LinearExpression.from_int(self.arg.ndim)

    @property
    def output_leaves(self) -> list[LinearExpression]:
        """Leaves this step introduces and that must be registered with the
        canonicalizer/step decomposer.

        Default: ``[self.output]`` (the output is itself a leaf). Steps that expose
        a non-leaf composite ``output`` (currently only partial-projection
        ``Relu``) override this to return the smaller fresh leaf actually
        introduced.
        """
        return [self.output]

    def update_bounds(self, lb, ub):
        raise NotImplementedError

    # TODO: figure out a better way to hash these objects
    def __hash__(self):
        return hash((self.arg, self.output))

output cached property 

output: LinearExpression

Left-hand side linear expression in the operator definition.

output_leaves property 

output_leaves: list[LinearExpression]

Leaves this step introduces and that must be registered with the canonicalizer/step decomposer.

Default: [self.output] (the output is itself a leaf). Steps that expose a non-leaf composite output (currently only partial-projection Relu) override this to return the smaller fresh leaf actually introduced.

Affine

Bases: Step

An affine operator of the following form:

\[ \phi(v,x) = \{ w \in \mathbf{R}^d \mid Mw = By,\ \ y = (v,x) \}, \]

where \(M \in \mathbf{R}^{ d \times d }\) is an invertible matrix and \(B \in \mathbf{R}^{d \times (p+d)}\).

It can also be constructed via weight and optional bias to represent the explicit transform \(w = W \cdot \text{arg} + b\), where arg is the input expression. In this form output_mat is set to identity internally.

Attributes:

Name Type Description
output LinearExpression

Left-hand side linear expression in the operator definition.
arg LinearExpression

Right-hand side linear expression in the operator definition.
lb ndarray

Lower bounds on the INPUT (arg), before the affine transform.
ub ndarray

Upper bounds on the INPUT (arg), before the affine transform.
relax_binary_vars bool

If True, binary decision variables associated with the operator are relaxed to take continuous values.
output_mat matrix

Coefficient matrix on the left-hand side of the operator equation.
obbt_output bool

True when output_mat is not the identity: OBBT is then also run on the output leaf after the affine equality is added to the OBBT LP (see StepDecomposer._process_primitive). For identity output_mat the output equals the arg, so sync_output_bounds already copies the OBBT-tightened input box and the extra LPs would be redundant.
Source code in fomtools/steps/affine.py 
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
class Affine(Step):
    r"""
    An affine operator of the following form:

    $$
      \phi(v,x) = \{ w \in \mathbf{R}^d \mid Mw = By,\ \ y = (v,x) \},
    $$

    where $M \in \mathbf{R}^{ d \times d }$ is an invertible matrix and
    $B \in \mathbf{R}^{d \times (p+d)}$.

    It can also be constructed via ``weight`` and optional ``bias`` to represent
    the explicit transform $w = W \cdot \text{arg} + b$, where ``arg`` is the
    input expression. In this form ``output_mat`` is set to identity internally.

    Attributes
    ----------
    output: LinearExpression
        Left-hand side linear expression in the operator definition.
    arg: LinearExpression
        Right-hand side linear expression in the operator definition.
    lb: np.ndarray
        Lower bounds on the INPUT (arg), before the affine transform.
    ub: np.ndarray
        Upper bounds on the INPUT (arg), before the affine transform.
    relax_binary_vars: bool
        If True, binary decision variables associated with the operator
        are relaxed to take continuous values.
    output_mat: np.matrix
        Coefficient matrix on the left-hand side of the operator equation.
    obbt_output: bool
        True when ``output_mat`` is not the identity: OBBT is then also run on the
        output leaf after the affine equality is added to the OBBT LP (see
        ``StepDecomposer._process_primitive``). For identity ``output_mat`` the
        output equals the arg, so ``sync_output_bounds`` already copies the
        OBBT-tightened input box and the extra LPs would be redundant.
    """

    def __init__(
        self,
        arg: LinearExpression,
        output_mat: np.matrix = None,
        weight: np.ndarray = None,
        bias: np.ndarray = None,
        relax_binary_vars: bool = False,
    ):
        if weight is not None:
            arg = weight @ arg
            if bias is not None:
                arg = arg + bias
            output_mat = np.eye(weight.shape[0])
        elif output_mat is None:
            raise ValueError("Either output_mat or weight must be provided.")
        super().__init__(arg, relax_binary_vars)
        self.output_mat = output_mat
        self.lb, self.ub = self.arg.bound_propagation()
        # from_int defaults lb/ub to scalars; force array form so subsequent
        # steps' bound_propagation() does not hit sparse @ scalar NaN errors.
        ndim = self.arg.ndim
        self.output.update_bounds(np.full(ndim, -np.inf), np.full(ndim, np.inf))

    @cached_property
    def obbt_output(self) -> bool:
        return not is_identity_matrix(self.output_mat)

    def update_bounds(self, lb: float | np.ndarray, ub: float | np.ndarray):
        """
        Update INPUT bounds (on arg) and propagate analytically to OUTPUT bounds on output.
        """
        self.lb, self.ub = lb, ub
        self.sync_output_bounds()

    def sync_output_bounds(self) -> None:
        """
        Tighten output bounds from the current input bounds.

        Interval-propagates ``output_mat_inverse`` (built by linear solves, never an
        explicit inverse). When ``output_mat`` is the identity, ``output == arg``, so
        the (possibly OBBT-tightened) input box ``self.lb``/``self.ub`` also bounds the
        output directly. ``output.update_bounds`` intersects with the existing output
        bounds, so tightening applied earlier is never discarded. For general
        ``output_mat`` the tightened input box is exploited by OBBT on the output leaf
        (see ``StepDecomposer._process_primitive``).
        """
        out_lb, out_ub = self.output_mat_inverse.bound_propagation()
        if is_identity_matrix(self.output_mat):
            out_lb = np.maximum(out_lb, self.lb)
            out_ub = np.minimum(out_ub, self.ub)
        self.output.update_bounds(lb=out_lb, ub=out_ub)
        # Repair the *combined* result (update_bounds intersects with the existing
        # output bounds, which is where noise crossings arise). Repairing before
        # update_bounds would be undone by the intersection, hence the direct write;
        # update_bounds is a plain elementwise max/min with no other side effects.
        self.output.lb, self.output.ub = repair_crossed_bounds(self.output.lb, self.output.ub)

    @cached_property
    def output_mat_inverse(self) -> LinearExpression:
        r"""
        Return a `LinearExpression` formed by applying $A^{-1}$, where $A$ is
        ``self.output_mat``, to each right-hand-side component $b_k$—i.e.,
        solving $A x_k = b_k$ for all $k$.
        Supports both dense and sparse ``output_mat``.
        """

        def _to_dense(b):
            return np.asarray(b.todense()) if spa.issparse(b) else np.asarray(b)

        if spa.issparse(self.output_mat):
            lu = spa.linalg.splu(self.output_mat.tocsc())

            def solve(b):
                return lu.solve(_to_dense(b))
        else:

            def solve(b):
                return np.asarray(np.linalg.solve(self.output_mat, _to_dense(b)))

        return LinearExpression(
            decomposition_dict={
                k: np.squeeze(solve(v)) for k, v in self.arg.decomposition_dict.items()
            },
            const=np.squeeze(solve(self.arg.const)),
        )

output_mat_inverse cached property 

output_mat_inverse: LinearExpression

Return a LinearExpression formed by applying \(A^{-1}\), where \(A\) is self.output_mat, to each right-hand-side component \(b_k\)—i.e., solving \(A x_k = b_k\) for all \(k\). Supports both dense and sparse output_mat.

sync_output_bounds 

sync_output_bounds() -> None

Tighten output bounds from the current input bounds.

Interval-propagates output_mat_inverse (built by linear solves, never an explicit inverse). When output_mat is the identity, output == arg, so the (possibly OBBT-tightened) input box self.lb/self.ub also bounds the output directly. output.update_bounds intersects with the existing output bounds, so tightening applied earlier is never discarded. For general output_mat the tightened input box is exploited by OBBT on the output leaf (see StepDecomposer._process_primitive).

Source code in fomtools/steps/affine.py 
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def sync_output_bounds(self) -> None:
    """
    Tighten output bounds from the current input bounds.

    Interval-propagates ``output_mat_inverse`` (built by linear solves, never an
    explicit inverse). When ``output_mat`` is the identity, ``output == arg``, so
    the (possibly OBBT-tightened) input box ``self.lb``/``self.ub`` also bounds the
    output directly. ``output.update_bounds`` intersects with the existing output
    bounds, so tightening applied earlier is never discarded. For general
    ``output_mat`` the tightened input box is exploited by OBBT on the output leaf
    (see ``StepDecomposer._process_primitive``).
    """
    out_lb, out_ub = self.output_mat_inverse.bound_propagation()
    if is_identity_matrix(self.output_mat):
        out_lb = np.maximum(out_lb, self.lb)
        out_ub = np.minimum(out_ub, self.ub)
    self.output.update_bounds(lb=out_lb, ub=out_ub)
    # Repair the *combined* result (update_bounds intersects with the existing
    # output bounds, which is where noise crossings arise). Repairing before
    # update_bounds would be undone by the intersection, hence the direct write;
    # update_bounds is a plain elementwise max/min with no other side effects.
    self.output.lb, self.output.ub = repair_crossed_bounds(self.output.lb, self.output.ub)

update_bounds 

update_bounds(lb: float | ndarray, ub: float | ndarray)

Update INPUT bounds (on arg) and propagate analytically to OUTPUT bounds on output.

Source code in fomtools/steps/affine.py 
77
78
79
80
81
82
def update_bounds(self, lb: float | np.ndarray, ub: float | np.ndarray):
    """
    Update INPUT bounds (on arg) and propagate analytically to OUTPUT bounds on output.
    """
    self.lb, self.ub = lb, ub
    self.sync_output_bounds()

Relu

Bases: Step

Rectified Linear Unit (ReLU) step operator, defined as:

\[ (\phi(v,x))_i = \mathbf{ReLU}(a^T y) = \max\{a^T y, 0\}. \]

It sets all negative values to zero while keeping the positive ones unchanged.

Attributes:

Name Type Description
arg LinearExpression

Right-hand side (input) linear expression.
lb ndarray

Lower bounds on the INPUT (arg), before ReLU is applied.
ub ndarray

Upper bounds on the INPUT (arg), before ReLU is applied.
relax_binary_vars bool

If True, binary decision variables associated with the operator are relaxed to take continuous values.
idx_with_binary_vars set[int]

Indices corresponding to binary decision variables.
proj_indices set[int]

Indices of output coordinates subjected to ReLU projection.
Notes

Partial projection (proj_ranges covers only a subset of the input coordinates) is encoded with a composite output: only len(proj_indices) fresh leaf variables are introduced (stored on self._output_proj_leaf), and the non-projected output coordinates share decomposition with arg directly. This keeps the SDP variable count, MIP variable count, and simulator leaf count at the minimum. For the full-projection default, output remains a single leaf of size arg.ndim (current behavior).

Source code in fomtools/steps/relu.py 
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
class Relu(Step):
    r"""
    Rectified Linear Unit (ReLU) step operator, defined as:

    $$
        (\phi(v,x))_i = \mathbf{ReLU}(a^T y) = \max\{a^T y, 0\}.
    $$

    It sets all negative values to zero while keeping the positive ones unchanged.

    Attributes
    ----------
    arg : LinearExpression
        Right-hand side (input) linear expression.
    lb : np.ndarray
        Lower bounds on the INPUT (arg), before ReLU is applied.
    ub : np.ndarray
        Upper bounds on the INPUT (arg), before ReLU is applied.
    relax_binary_vars : bool
        If True, binary decision variables associated with the operator
        are relaxed to take continuous values.
    idx_with_binary_vars : set[int]
        Indices corresponding to binary decision variables.
    proj_indices : set[int]
        Indices of output coordinates subjected to ReLU projection.

    Notes
    -----
    Partial projection (``proj_ranges`` covers only a subset of the input
    coordinates) is encoded with a *composite* ``output``: only
    ``len(proj_indices)`` fresh leaf variables are introduced (stored on
    ``self._output_proj_leaf``), and the non-projected output coordinates
    share decomposition with ``arg`` directly. This keeps the SDP variable
    count, MIP variable count, and simulator leaf count at the minimum.
    For the full-projection default, ``output`` remains a single leaf of size
    ``arg.ndim`` (current behavior).
    """

    def __init__(
        self,
        arg: LinearExpression,
        proj_ranges=None,
        relax_binary_vars=False,
    ):
        super().__init__(arg, relax_binary_vars)
        self.proj_indices = (
            set(range(0, self.arg.ndim)) if proj_ranges is None else self._get_ranges(proj_ranges)
        )
        self.idx_with_binary_vars = set()
        # When True, the SDP canonicalizer does not register this ReLU as an
        # orthant prox for the subgradient-monotonicity cut. Set by decomposers
        # that build a composite operator (e.g. soft-threshold) from internal
        # ReLUs and register the composite prox instead. Declared here as a formal
        # attribute (rather than a dynamic one) so it survives copying.
        self._skip_prox_registration = False

        arg_lb, arg_ub = self.arg.bound_propagation()

        if self.proj_indices != set(range(self.arg.ndim)):
            # Initialize the proj-leaf with array-valued bounds so the composite
            # output's bound_propagation has a consistent shape to work with.
            proj_list = sorted(self.proj_indices)
            proj_lb = np.maximum(0, np.asarray(arg_lb)[proj_list])
            proj_ub = np.maximum(0, np.asarray(arg_ub)[proj_list])
            self._output_proj_leaf = LinearExpression.from_int(
                len(self.proj_indices), lb=proj_lb, ub=proj_ub
            )
            # Override Step.output cached_property with a composite expression
            # whose non-proj rows share storage with arg and whose proj rows
            # come from _output_proj_leaf.
            self.__dict__["output"] = self._build_composite_output()
        else:
            self._output_proj_leaf = None

        self.update_bounds(arg_lb, arg_ub)

    def __hash__(self):
        return hash((self.arg, self.output))

    @staticmethod
    def _get_ranges(ranges: Tuple[int, int] | List[Tuple[int, int]]) -> set[int]:
        """
        Return the set of integers covered by one or more half-open intervals.

        Parameters
        ----------
        ranges: Tuple[int, int] or List[Tuple[int, int]]
            Either a single (start, end) pair or an iterable of such pairs.
            Each pair is interpreted as a half-open interval ``[start, end)``.

        Returns
        -------
        set[int]
            The unique integers contained in the given half-open range(s).

        Examples
        --------
        >>> Relu._get_ranges(((1, 2), (5, 9)))
        {1, 5, 6, 7, 8}

        >>> Relu._get_ranges((4, 7))
        {4, 5, 6}
        """
        # Try the iterable-of-pairs form first. If ``ranges`` is a single
        # ``(start, end)`` pair, unpacking its int elements raises TypeError
        # ("cannot unpack non-iterable int") — fall back to the single-range
        # form. This avoids the pre-existing ``ranges.append((0, 0))``
        # workaround which silently mutated caller-owned lists.
        try:
            return {i for (start, end) in ranges for i in range(start, end)}
        except TypeError:
            return set(range(ranges[0], ranges[1]))

    def _build_composite_output(self) -> LinearExpression:
        """Construct the composite output: proj rows from ``_output_proj_leaf``,
        non-proj rows passed through directly from ``arg``.

        The composite has dimension ``arg.ndim`` and its decomposition_dict
        references both ``_output_proj_leaf`` and the leaves of ``arg``.
        """
        n = self.arg.ndim
        proj_list = sorted(self.proj_indices)
        nonproj_list = sorted(set(range(n)) - self.proj_indices)
        k = len(proj_list)

        # Selection matrix placing _output_proj_leaf (size k) into proj rows of output.
        S_proj = spa.csc_matrix(
            (np.ones(k), (np.array(proj_list, dtype=int), np.arange(k))),
            shape=(n, k),
        )

        # Selection matrix masking out proj rows of arg; non-proj rows of arg
        # pass through to output unchanged.
        m = len(nonproj_list)
        if m > 0:
            nonproj_arr = np.array(nonproj_list, dtype=int)
            S_nonproj = spa.csc_matrix(
                (np.ones(m), (nonproj_arr, nonproj_arr)),
                shape=(n, n),
            )
        else:
            S_nonproj = spa.csc_matrix((n, n))

        decomposition_dict: dict = {self._output_proj_leaf: S_proj}
        for leaf, coeff in self.arg.decomposition_dict.items():
            composed = S_nonproj @ coeff
            if leaf in decomposition_dict:
                decomposition_dict[leaf] = decomposition_dict[leaf] + composed
            else:
                decomposition_dict[leaf] = composed

        const = S_nonproj @ self.arg.const

        # Drop entries whose composed coefficient is all-zero — happens when
        # an arg leaf contributes only to projected rows, since ``S_nonproj``
        # zeros out those rows. Pruning prevents downstream canonicalizers
        # from looking up an irrelevant leaf in ``var_idx_map``.
        result = LinearExpression(decomposition_dict=prune_dict(decomposition_dict), const=const)
        # Initialize array-valued bounds so later update_bounds calls have a
        # consistent shape to maximize/minimize against.
        result.lb, result.ub = result.bound_propagation()
        return result

    def relu(self, x: np.ndarray) -> np.ndarray:
        """
        Apply ReLU activation selectively to specified indices.

        Parameters
        ----------
        x: np.ndarray
            A 1-D array of values on which to apply the activation.

        Returns
        -------
        np.ndarray
            A new array where elements at positions specified in
            ``self.proj_indices`` are passed through ReLU, while
            all other elements are left unchanged.
        """
        proj_mask = np.isin(np.arange(len(x)), list(self.proj_indices))
        return np.where(proj_mask, np.maximum(0, x), x)

    @property
    def output_leaves(self) -> list[LinearExpression]:
        """For partial projection, only ``_output_proj_leaf`` is fresh; non-proj
        rows of the composite output reuse the arg's leaves. For full projection,
        the output is itself a single fresh leaf.
        """
        if self._output_proj_leaf is not None:
            return [self._output_proj_leaf]
        return [self.output]

    def update_bounds(self, lb: float | np.ndarray, ub: float | np.ndarray):
        """
        Update the bounds of the left-hand side after applying
        the ReLU activation.

        Parameters
        ----------
        lb: float or np.ndarray
            The new lower bound(s) before applying ReLU.

        ub: float or np.ndarray
            The new upper bound(s) before applying ReLU.

        """
        self.lb, self.ub = lb, ub
        self.output.update_bounds(lb=self.relu(lb), ub=self.relu(ub))
        if self._output_proj_leaf is not None:
            proj_list = sorted(self.proj_indices)
            proj_lb = np.maximum(0, np.asarray(lb)[proj_list])
            proj_ub = np.maximum(0, np.asarray(ub)[proj_list])
            self._output_proj_leaf.update_bounds(proj_lb, proj_ub)

output_leaves property 

output_leaves: list[LinearExpression]

For partial projection, only _output_proj_leaf is fresh; non-proj rows of the composite output reuse the arg's leaves. For full projection, the output is itself a single fresh leaf.

relu 

relu(x: ndarray) -> ndarray

Apply ReLU activation selectively to specified indices.

Parameters:

Name Type Description Default
x ndarray

A 1-D array of values on which to apply the activation.

 required

Returns:

Type Description
ndarray

A new array where elements at positions specified in self.proj_indices are passed through ReLU, while all other elements are left unchanged.
Source code in fomtools/steps/relu.py 
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
def relu(self, x: np.ndarray) -> np.ndarray:
    """
    Apply ReLU activation selectively to specified indices.

    Parameters
    ----------
    x: np.ndarray
        A 1-D array of values on which to apply the activation.

    Returns
    -------
    np.ndarray
        A new array where elements at positions specified in
        ``self.proj_indices`` are passed through ReLU, while
        all other elements are left unchanged.
    """
    proj_mask = np.isin(np.arange(len(x)), list(self.proj_indices))
    return np.where(proj_mask, np.maximum(0, x), x)

update_bounds 

update_bounds(lb: float | ndarray, ub: float | ndarray)

Update the bounds of the left-hand side after applying the ReLU activation.

Parameters:

Name Type Description Default
lb float | ndarray

The new lower bound(s) before applying ReLU.

 required
ub float | ndarray

The new upper bound(s) before applying ReLU.

 required
Source code in fomtools/steps/relu.py 
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
def update_bounds(self, lb: float | np.ndarray, ub: float | np.ndarray):
    """
    Update the bounds of the left-hand side after applying
    the ReLU activation.

    Parameters
    ----------
    lb: float or np.ndarray
        The new lower bound(s) before applying ReLU.

    ub: float or np.ndarray
        The new upper bound(s) before applying ReLU.

    """
    self.lb, self.ub = lb, ub
    self.output.update_bounds(lb=self.relu(lb), ub=self.relu(ub))
    if self._output_proj_leaf is not None:
        proj_list = sorted(self.proj_indices)
        proj_lb = np.maximum(0, np.asarray(lb)[proj_list])
        proj_ub = np.maximum(0, np.asarray(ub)[proj_list])
        self._output_proj_leaf.update_bounds(proj_lb, proj_ub)

SaturatedLinear

Bases: Step

Saturated Linear Unit step operator, defined as:

\[ (\phi(v, x))_i = S_{[x, c]}(a^T y) = \begin{cases} x & a^T y < x, \\\\ a^T y & x \leq a^T y \leq c, \\\\ c & a^T y > c, \end{cases} \quad \text{with } y = (v, x), \]

where \(x, c \in \mathbb{R}\).

Attributes:

Name Type Description
arg LinearExpression

Right-hand side (input) linear expression.
output LinearExpression

Left-hand side (output) linear expression.
lb ndarray

Lower bounds on the INPUT (arg), before saturation is applied.
ub ndarray

Upper bounds on the INPUT (arg), before saturation is applied.
relax_binary_vars bool

If True, binary decision variables associated with the operator are relaxed to take continuous values.
l ndarray[float] | LinearExpression

The lower saturation bound of the saturated linear operator. Can be either a numeric value or a LinearExpression object whose lower bound acts as the lower saturation bound.
u ndarray[float] | LinearExpression

The upper saturation bound of the saturated linear operator. Can be either a numeric value or a LinearExpression object whose upper bound acts as the upper saturation bound.
Source code in fomtools/steps/saturated_linear.py 
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
class SaturatedLinear(Step):
    r"""
    Saturated Linear Unit step operator, defined as:

    $$
        (\phi(v, x))_i = S_{[x, c]}(a^T y) =
        \begin{cases}
            x & a^T y < x, \\\\
            a^T y & x \leq a^T y \leq c, \\\\
            c & a^T y > c,
        \end{cases}
        \quad \text{with } y = (v, x),
    $$

    where $x, c \in \mathbb{R}$.

    Attributes
    ----------
    arg : LinearExpression
        Right-hand side (input) linear expression.
    output : LinearExpression
        Left-hand side (output) linear expression.
    lb : np.ndarray
        Lower bounds on the INPUT (arg), before saturation is applied.
    ub : np.ndarray
        Upper bounds on the INPUT (arg), before saturation is applied.
    relax_binary_vars : bool
        If True, binary decision variables associated with the operator
        are relaxed to take continuous values.
    l : np.ndarray[float] | LinearExpression
        The lower saturation bound of the saturated linear operator.
        Can be either a numeric value or a LinearExpression object whose
        lower bound acts as the lower saturation bound.
    u : np.ndarray[float] | LinearExpression
        The upper saturation bound of the saturated linear operator.
        Can be either a numeric value or a LinearExpression object whose
        upper bound acts as the upper saturation bound.
    """

    def __init__(
        self,
        arg: LinearExpression,
        l: LinearExpression | np.ndarray[float],
        u: LinearExpression | np.ndarray[float],
        relax_binary_vars: bool = False,
    ):
        super().__init__(arg, relax_binary_vars)
        self.l = l
        self.u = u
        self.update_bounds(*self.arg.bound_propagation())

    def __hash__(self):
        return hash((self.arg, self.output))

    def satlin_prop(self, x: np.ndarray, bound: str) -> np.ndarray:
        """
        Propagate a bound through the saturated linear operator.

        Since clip(x, l, u) is monotonically non-decreasing in all three
        arguments, the correct output bounds are:
            lower: clip(x_lb, l_lb, u_lb)
            upper: clip(x_ub, l_ub, u_ub)

        Parameters
        ----------
        x: np.ndarray
            A 1-D array of input bound values to propagate.
        bound: str
            Either "lb" or "ub", indicating which bound is being propagated.

        Returns
        -------
        np.ndarray
            The propagated bound on the output of the saturated linear operator.
        """
        try:
            l_bound = self.l.lb if bound == "lb" else self.l.ub
        except AttributeError:
            l_bound = self.l

        try:
            u_bound = self.u.lb if bound == "lb" else self.u.ub
        except AttributeError:
            u_bound = self.u
        return np.minimum(np.maximum(x, l_bound), u_bound)

    def update_bounds(self, lb: float | np.ndarray, ub: float | np.ndarray):
        """
        Update the bounds of the left-hand side after applying
        the Saturated Linear Unit activation.

        Parameters
        ----------
        lb: float or np.ndarray
            The new lower bound(s) before applying the Saturated Linear
            Unit.

        ub: float or np.ndarray
            The new upper bound(s) before applying the Saturated Linear
            Unit.

        """
        self.lb, self.ub = lb, ub
        self.output.update_bounds(
            lb=self.satlin_prop(lb, "lb"),
            ub=self.satlin_prop(ub, "ub"),
        )

satlin_prop 

satlin_prop(x: ndarray, bound: str) -> ndarray

Propagate a bound through the saturated linear operator.

Since clip(x, l, u) is monotonically non-decreasing in all three arguments, the correct output bounds are: lower: clip(x_lb, l_lb, u_lb) upper: clip(x_ub, l_ub, u_ub)

Parameters:

Name Type Description Default
x ndarray

A 1-D array of input bound values to propagate.

 required
bound str

Either "lb" or "ub", indicating which bound is being propagated.

 required

Returns:

Type Description
ndarray

The propagated bound on the output of the saturated linear operator.
Source code in fomtools/steps/saturated_linear.py 
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
def satlin_prop(self, x: np.ndarray, bound: str) -> np.ndarray:
    """
    Propagate a bound through the saturated linear operator.

    Since clip(x, l, u) is monotonically non-decreasing in all three
    arguments, the correct output bounds are:
        lower: clip(x_lb, l_lb, u_lb)
        upper: clip(x_ub, l_ub, u_ub)

    Parameters
    ----------
    x: np.ndarray
        A 1-D array of input bound values to propagate.
    bound: str
        Either "lb" or "ub", indicating which bound is being propagated.

    Returns
    -------
    np.ndarray
        The propagated bound on the output of the saturated linear operator.
    """
    try:
        l_bound = self.l.lb if bound == "lb" else self.l.ub
    except AttributeError:
        l_bound = self.l

    try:
        u_bound = self.u.lb if bound == "lb" else self.u.ub
    except AttributeError:
        u_bound = self.u
    return np.minimum(np.maximum(x, l_bound), u_bound)

update_bounds 

update_bounds(lb: float | ndarray, ub: float | ndarray)

Update the bounds of the left-hand side after applying the Saturated Linear Unit activation.

Parameters:

Name Type Description Default
lb float | ndarray

The new lower bound(s) before applying the Saturated Linear Unit.

 required
ub float | ndarray

The new upper bound(s) before applying the Saturated Linear Unit.

 required
Source code in fomtools/steps/saturated_linear.py 
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
def update_bounds(self, lb: float | np.ndarray, ub: float | np.ndarray):
    """
    Update the bounds of the left-hand side after applying
    the Saturated Linear Unit activation.

    Parameters
    ----------
    lb: float or np.ndarray
        The new lower bound(s) before applying the Saturated Linear
        Unit.

    ub: float or np.ndarray
        The new upper bound(s) before applying the Saturated Linear
        Unit.

    """
    self.lb, self.ub = lb, ub
    self.output.update_bounds(
        lb=self.satlin_prop(lb, "lb"),
        ub=self.satlin_prop(ub, "ub"),
    )

SoftThreshold

Bases: Step

A soft-thresholding operator for shrinkage, defined as:

\[ (\phi(v,x))_i=\mathcal{T}_{\lambda}(a^\top y)= \begin{cases} a^\top y-\lambda,& a^\top y>\lambda,\\[6pt] 0,& |a^\top y|\le \lambda,\\[6pt] a^\top y+\lambda,& a^\top y<-\lambda, \end{cases} \quad \text{with } y=(v,x), \]

where \(\mathcal{T}_{\lambda}\) is the soft-thresholding operator with constant \(\lambda>0\), and \(a\in\mathbf{R}^d\).

Attributes:

Name Type Description
arg LinearExpression

Right-hand side (input) linear expression.
output LinearExpression

Left-hand side (output) linear expression.
lb ndarray

Lower bounds on the INPUT (arg), before soft-thresholding is applied.
ub ndarray

Upper bounds on the INPUT (arg), before soft-thresholding is applied.
relax_binary_vars bool

If True, binary decision variables associated with the operator are relaxed to take continuous values.
lambd float

The threshold (\(\lambda\)) controlling shrinkage. Must be non-negative.
idx_with_right_binary_vars set[int]

Indices tied to "right-branch" indicators in a MILP reformulation (e.g., \(t>\lambda\)).
idx_with_left_binary_vars set[int]

Indices tied to "left-branch" indicators (e.g., \(t<-\lambda\)).
Source code in fomtools/steps/soft_threshold.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
class SoftThreshold(Step):
    r"""
    A soft-thresholding operator for shrinkage, defined as:

    $$
        (\phi(v,x))_i=\mathcal{T}_{\lambda}(a^\top y)=
        \begin{cases}
        a^\top y-\lambda,& a^\top y>\lambda,\\[6pt]
        0,& |a^\top y|\le \lambda,\\[6pt]
        a^\top y+\lambda,& a^\top y<-\lambda,
        \end{cases}
        \quad \text{with } y=(v,x),
    $$

    where $\mathcal{T}_{\lambda}$ is the soft-thresholding operator with constant $\lambda>0$,
    and $a\in\mathbf{R}^d$.

    Attributes
    ----------
    arg : LinearExpression
        Right-hand side (input) linear expression.
    output : LinearExpression
        Left-hand side (output) linear expression.
    lb : np.ndarray
        Lower bounds on the INPUT (arg), before soft-thresholding is applied.
    ub : np.ndarray
        Upper bounds on the INPUT (arg), before soft-thresholding is applied.
    relax_binary_vars : bool
        If True, binary decision variables associated with the operator
        are relaxed to take continuous values.
    lambd : float
        The threshold ($\lambda$) controlling shrinkage. Must be non-negative.
    idx_with_right_binary_vars : set[int]
        Indices tied to "right-branch" indicators in a MILP reformulation (e.g., $t>\lambda$).
    idx_with_left_binary_vars : set[int]
        Indices tied to "left-branch" indicators (e.g., $t<-\lambda$).
    """

    def __init__(self, arg: LinearExpression, lambd: float, relax_binary_vars: bool = False):
        super().__init__(arg, relax_binary_vars)
        self.lambd = lambd
        self.idx_with_right_binary_vars = set()
        self.idx_with_left_binary_vars = set()
        self.update_bounds(*self.arg.bound_propagation())

    def __hash__(self):
        return hash((self.arg, self.output))

    def soft_threshold(self, x: np.ndarray) -> np.ndarray:
        """
        Apply the soft-thresholding operator to the input array.

        Parameters
        ----------
        x: np.ndarray
            A 1-D array of values on which to apply the activation.

        Returns
        -------
        np.ndarray
            A new array obtained by applying the soft-thresholding
            operator element-wise. Values with magnitude less than `self.lambd`
            are mapped to zero, while others are shrunk toward zero by
            `self.lambd`.

        """
        return np.sign(x) * np.maximum(np.abs(x) - self.lambd, 0)

    def update_bounds(self, lb: float | np.ndarray, ub: float | np.ndarray):
        """
        Update the bounds of the left-hand side after applying
        the soft-thresholding operator.

        Parameters
        ----------
        lb: float or np.ndarray
            The new lower bound(s) before applying soft-thresholding.

        ub: float or np.ndarray
            The new upper bound(s) before applying soft-thresholding.
        """
        self.lb, self.ub = lb, ub
        self.output.update_bounds(lb=self.soft_threshold(lb), ub=self.soft_threshold(ub))

soft_threshold 

soft_threshold(x: ndarray) -> ndarray

Apply the soft-thresholding operator to the input array.

Parameters:

Name Type Description Default
x ndarray

A 1-D array of values on which to apply the activation.

 required

Returns:

Type Description
ndarray

A new array obtained by applying the soft-thresholding operator element-wise. Values with magnitude less than self.lambd are mapped to zero, while others are shrunk toward zero by self.lambd.
Source code in fomtools/steps/soft_threshold.py 
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
def soft_threshold(self, x: np.ndarray) -> np.ndarray:
    """
    Apply the soft-thresholding operator to the input array.

    Parameters
    ----------
    x: np.ndarray
        A 1-D array of values on which to apply the activation.

    Returns
    -------
    np.ndarray
        A new array obtained by applying the soft-thresholding
        operator element-wise. Values with magnitude less than `self.lambd`
        are mapped to zero, while others are shrunk toward zero by
        `self.lambd`.

    """
    return np.sign(x) * np.maximum(np.abs(x) - self.lambd, 0)

update_bounds 

update_bounds(lb: float | ndarray, ub: float | ndarray)

Update the bounds of the left-hand side after applying the soft-thresholding operator.

Parameters:

Name Type Description Default
lb float | ndarray

The new lower bound(s) before applying soft-thresholding.

 required
ub float | ndarray

The new upper bound(s) before applying soft-thresholding.

 required
Source code in fomtools/steps/soft_threshold.py 
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
def update_bounds(self, lb: float | np.ndarray, ub: float | np.ndarray):
    """
    Update the bounds of the left-hand side after applying
    the soft-thresholding operator.

    Parameters
    ----------
    lb: float or np.ndarray
        The new lower bound(s) before applying soft-thresholding.

    ub: float or np.ndarray
        The new upper bound(s) before applying soft-thresholding.
    """
    self.lb, self.ub = lb, ub
    self.output.update_bounds(lb=self.soft_threshold(lb), ub=self.soft_threshold(ub))

QuadraticProgram

Bases: Step

Step representing the solution of an arbitrary convex quadratic program

\[ x^\star = \arg\min_x\; \tfrac12 x^\top P x + q^\top x \quad\text{s.t.}\quad A x = b,\; C x \le d, \]

where \(P \in \mathbf{R}^{n\times n}\) is a fixed symmetric PSD matrix (dense numpy or scipy.sparse), and q, b, d may be numeric data or LinearExpression parameters (carrying dependence on prior iterates). The verification-side lowering (decompose_quadraticprogram in the step decomposer) posts the KKT system

\[ P x + q + A^\top \nu + C^\top \lambda = 0,\quad A x = b,\quad C x \le d,\quad 0 \le \lambda \perp (d - C x) \ge 0, \]

with \(\nu\) a free equality multiplier and \(\lambda \ge 0\), using the complementarity primitive. The MIP big-M encoding of the complementarity needs a finite upper bound on \(\lambda\); supply it via dual_ub whenever inequalities are present (OBBT tightens it further when enabled).

Invoked via Verifier.quadratic_program rather than the single-arg run_step path: an arbitrary QP has up to three parametric inputs (q, b, d), and each LinearExpression already references the prior iterates it was built from, so no input is forced into an arg slot.

Parameters:

Name Type Description Default
P np.ndarray or scipy.sparse matrix

The (n, n) quadratic term; must be symmetric and PSD.

 required
q None, float, np.ndarray, or LinearExpression

Linear cost of dimension n. If None (projection convenience), q defaults to -arg in the decomposer, which requires arg to be a dim-n LinearExpression.

 None
A ndarray or sparse

Equality matrix of shape (p, n). None for no equality constraints.

 None
b float, np.ndarray, or LinearExpression

Equality right-hand side of dimension p.

 None
C ndarray or sparse

Inequality matrix of shape (m, n). None for no inequalities.

 None
d float, np.ndarray, or LinearExpression

Inequality right-hand side of dimension m.

 None
dual_ub None, float, or np.ndarray

Finite upper bound(s) on the inequality dual lambda >= 0. Required (raised eagerly in the decomposer) whenever C is present.

 None
primal_bounds tuple[ndarray, ndarray]

(lb, ub) on the primal x. Needed when P is singular / not positive definite (the stationarity row cannot pin x bounds); otherwise derived analytically and tightened by OBBT.

 None
nu_bounds tuple[ndarray, ndarray]

(lb, ub) on the free equality multiplier nu. Defaults to +/- inf (nu only appears in the linear stationarity row).

 None
arg LinearExpression

The projection point used when q is None (then q := -arg). Optional otherwise; retained only as the nominal Step input.

 None
relax_binary_vars bool

Forwarded to the Step interface (relax big-M binaries when set).

 False
Source code in fomtools/steps/quad_prog.py 
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
class QuadraticProgram(Step):
    r"""
    Step representing the solution of an arbitrary convex quadratic program

    $$
        x^\star = \arg\min_x\; \tfrac12 x^\top P x + q^\top x
        \quad\text{s.t.}\quad A x = b,\; C x \le d,
    $$

    where $P \in \mathbf{R}^{n\times n}$ is a fixed **symmetric PSD** matrix
    (dense ``numpy`` or ``scipy.sparse``), and ``q``, ``b``, ``d`` may be numeric
    data or `LinearExpression` parameters (carrying dependence on prior
    iterates). The verification-side lowering (``decompose_quadraticprogram`` in
    the step decomposer) posts the KKT system

    $$
        P x + q + A^\top \nu + C^\top \lambda = 0,\quad
        A x = b,\quad C x \le d,\quad
        0 \le \lambda \perp (d - C x) \ge 0,
    $$

    with $\nu$ a free equality multiplier and $\lambda \ge 0$, using
    the complementarity primitive. The MIP big-M encoding of the complementarity
    needs a finite upper bound on $\lambda$; supply it via ``dual_ub``
    whenever inequalities are present (OBBT tightens it further when enabled).

    Invoked via `Verifier.quadratic_program` rather than the single-``arg``
    ``run_step`` path: an arbitrary QP has up to three parametric inputs
    (``q``, ``b``, ``d``), and each ``LinearExpression`` already references the
    prior iterates it was built from, so no input is forced into an ``arg`` slot.

    Parameters
    ----------
    P : np.ndarray or scipy.sparse matrix
        The ``(n, n)`` quadratic term; must be symmetric and PSD.
    q : None, float, np.ndarray, or LinearExpression, optional
        Linear cost of dimension ``n``. If ``None`` (projection convenience),
        ``q`` defaults to ``-arg`` in the decomposer, which requires ``arg`` to be
        a dim-``n`` ``LinearExpression``.
    A : np.ndarray or scipy.sparse, optional
        Equality matrix of shape ``(p, n)``. ``None`` for no equality constraints.
    b : float, np.ndarray, or LinearExpression, optional
        Equality right-hand side of dimension ``p``.
    C : np.ndarray or scipy.sparse, optional
        Inequality matrix of shape ``(m, n)``. ``None`` for no inequalities.
    d : float, np.ndarray, or LinearExpression, optional
        Inequality right-hand side of dimension ``m``.
    dual_ub : None, float, or np.ndarray, optional
        Finite upper bound(s) on the inequality dual ``lambda >= 0``. Required
        (raised eagerly in the decomposer) whenever ``C`` is present.
    primal_bounds : tuple[np.ndarray, np.ndarray], optional
        ``(lb, ub)`` on the primal ``x``. Needed when ``P`` is singular / not
        positive definite (the stationarity row cannot pin ``x`` bounds);
        otherwise derived analytically and tightened by OBBT.
    nu_bounds : tuple[np.ndarray, np.ndarray], optional
        ``(lb, ub)`` on the free equality multiplier ``nu``. Defaults to
        ``+/- inf`` (``nu`` only appears in the linear stationarity row).
    arg : LinearExpression, optional
        The projection point used when ``q is None`` (then ``q := -arg``).
        Optional otherwise; retained only as the nominal Step input.
    relax_binary_vars : bool, default False
        Forwarded to the Step interface (relax big-M binaries when set).
    """

    skip_obbt: bool = True

    def __init__(
        self,
        P,
        q=None,
        A=None,
        b=None,
        C=None,
        d=None,
        dual_ub=None,
        primal_bounds=None,
        nu_bounds=None,
        arg: LinearExpression = None,
        relax_binary_vars: bool = False,
    ):
        # --- validate P: square, symmetric, PSD (dense or sparse) ---
        if P.ndim != 2 or P.shape[0] != P.shape[1]:
            raise ValueError(f"P must be a square 2-D matrix; got shape {P.shape}")
        if spa.issparse(P):
            P = P.tocsr()  # CSR supports T, arithmetic, and @ (LinearExpression.__rmatmul__)
        n = P.shape[0]
        sym_err = _sym_error(P)
        if sym_err > _PSD_TOL:
            raise ValueError(f"P must be symmetric; max|P - P^T| = {sym_err:.3e}")
        min_eig = _min_eig(P)
        if min_eig < -_PSD_TOL:
            raise ValueError(
                f"P must be positive semidefinite; smallest eigenvalue {min_eig:.3e} < -{_PSD_TOL:g}"
            )

        # --- validate parametric data dimensions ---
        if q is None:
            if arg is None or not _is_linexpr(arg) or arg.ndim != n:
                raise ValueError(
                    "q=None defaults to -arg, which requires arg to be a "
                    f"LinearExpression of dimension n={n}."
                )
        elif _vec_dim(q, n) != n:
            raise ValueError(f"q dim {_vec_dim(q, n)} must match P dimension n={n}")

        if A is not None:
            if A.shape[1] != n:
                raise ValueError(f"A must have {n} columns to match P; got shape {A.shape}")
            p = A.shape[0]
            if b is None:
                raise ValueError("b must be provided when A is given")
            if _vec_dim(b, p) != p:
                raise ValueError(f"b dim {_vec_dim(b, p)} must match A rows {p}")
        if C is not None:
            if C.shape[1] != n:
                raise ValueError(f"C must have {n} columns to match P; got shape {C.shape}")
            m = C.shape[0]
            if d is None:
                raise ValueError("d must be provided when C is given")
            if _vec_dim(d, m) != m:
                raise ValueError(f"d dim {_vec_dim(d, m)} must match C rows {m}")

        # Bypass Step.__init__ (it dereferences arg.ndim): arg is optional here.
        self.arg = arg
        self.relax_binary_vars = relax_binary_vars
        self.n = n
        self.proj_indices = set(range(n))

        self.P = P
        self.q = q
        self.A = A
        self.b = b
        self.C = C
        self.d = d
        self.dual_ub = dual_ub
        self.primal_bounds = primal_bounds
        self.nu_bounds = nu_bounds

    @cached_property
    def output(self) -> LinearExpression:
        """Output is the QP primal solution ``x`` of dimension ``n`` (``!= arg.ndim``)."""
        return LinearExpression.from_int(self.n)

    def update_bounds(self, lb, ub):
        """Store primal bounds (interface compliance; QP is a composite step that
        bypasses ``_process_primitive``, so this is rarely invoked)."""
        self.lb, self.ub = lb, ub
        self.output.update_bounds(lb=lb, ub=ub)

    def __hash__(self):
        return hash((self.arg, self.output))

    def solve(self, q=None, b=None, d=None) -> np.ndarray:
        r"""
        Numerically solve the QP with cvxpy (reference forward-evaluation).

        Numeric overrides for ``q``/``b``/``d`` may be supplied for the case where
        the stored values are ``LinearExpression`` parameters; otherwise the stored
        numeric data is used. ``P`` is densified for cvxpy.

        Returns
        -------
        np.ndarray
            The optimal primal ``x`` of dimension ``n``.
        """
        import cvxpy as cp

        def _resolve(override, stored, name):
            value = override if override is not None else stored
            if value is None:
                return None
            if _is_linexpr(value):
                raise ValueError(
                    f"{name} is a LinearExpression; pass a numeric override to solve()"
                )
            return np.asarray(value, dtype=float)

        P_dense = P.toarray() if spa.issparse(P := self.P) else np.asarray(P, dtype=float)
        q_val = _resolve(q, self.q, "q")
        if q_val is None:
            raise ValueError("q is None (defaults to -arg); pass a numeric q to solve()")
        q_val = np.broadcast_to(q_val, (self.n,)).astype(float)

        x = cp.Variable(self.n)
        objective = q_val @ x
        if np.any(P_dense):
            objective = 0.5 * cp.quad_form(x, cp.psd_wrap(P_dense)) + objective
        constraints = []
        if self.A is not None:
            A = self.A.toarray() if spa.issparse(self.A) else np.asarray(self.A, dtype=float)
            b_val = _resolve(b, self.b, "b")
            constraints.append(A @ x == np.broadcast_to(b_val, (A.shape[0],)))
        if self.C is not None:
            C = self.C.toarray() if spa.issparse(self.C) else np.asarray(self.C, dtype=float)
            d_val = _resolve(d, self.d, "d")
            constraints.append(C @ x <= np.broadcast_to(d_val, (C.shape[0],)))
        # Pin the conic solver (cvxpy's default may be an unavailable/expired backend).
        cp.Problem(cp.Minimize(objective), constraints).solve(solver=cp.CLARABEL)
        return np.asarray(x.value, dtype=float)

output cached property 

output: LinearExpression

Output is the QP primal solution x of dimension n (!= arg.ndim).

solve 

solve(q=None, b=None, d=None) -> ndarray

Numerically solve the QP with cvxpy (reference forward-evaluation).

Numeric overrides for q/b/d may be supplied for the case where the stored values are LinearExpression parameters; otherwise the stored numeric data is used. P is densified for cvxpy.

Returns:

Type Description
ndarray

The optimal primal x of dimension n.
Source code in fomtools/steps/quad_prog.py 
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
def solve(self, q=None, b=None, d=None) -> np.ndarray:
    r"""
    Numerically solve the QP with cvxpy (reference forward-evaluation).

    Numeric overrides for ``q``/``b``/``d`` may be supplied for the case where
    the stored values are ``LinearExpression`` parameters; otherwise the stored
    numeric data is used. ``P`` is densified for cvxpy.

    Returns
    -------
    np.ndarray
        The optimal primal ``x`` of dimension ``n``.
    """
    import cvxpy as cp

    def _resolve(override, stored, name):
        value = override if override is not None else stored
        if value is None:
            return None
        if _is_linexpr(value):
            raise ValueError(
                f"{name} is a LinearExpression; pass a numeric override to solve()"
            )
        return np.asarray(value, dtype=float)

    P_dense = P.toarray() if spa.issparse(P := self.P) else np.asarray(P, dtype=float)
    q_val = _resolve(q, self.q, "q")
    if q_val is None:
        raise ValueError("q is None (defaults to -arg); pass a numeric q to solve()")
    q_val = np.broadcast_to(q_val, (self.n,)).astype(float)

    x = cp.Variable(self.n)
    objective = q_val @ x
    if np.any(P_dense):
        objective = 0.5 * cp.quad_form(x, cp.psd_wrap(P_dense)) + objective
    constraints = []
    if self.A is not None:
        A = self.A.toarray() if spa.issparse(self.A) else np.asarray(self.A, dtype=float)
        b_val = _resolve(b, self.b, "b")
        constraints.append(A @ x == np.broadcast_to(b_val, (A.shape[0],)))
    if self.C is not None:
        C = self.C.toarray() if spa.issparse(self.C) else np.asarray(self.C, dtype=float)
        d_val = _resolve(d, self.d, "d")
        constraints.append(C @ x <= np.broadcast_to(d_val, (C.shape[0],)))
    # Pin the conic solver (cvxpy's default may be an unavailable/expired backend).
    cp.Problem(cp.Minimize(objective), constraints).solve(solver=cp.CLARABEL)
    return np.asarray(x.value, dtype=float)

update_bounds 

update_bounds(lb, ub)

Store primal bounds (interface compliance; QP is a composite step that bypasses _process_primitive, so this is rarely invoked).

Source code in fomtools/steps/quad_prog.py 
169
170
171
172
173
def update_bounds(self, lb, ub):
    """Store primal bounds (interface compliance; QP is a composite step that
    bypasses ``_process_primitive``, so this is rarely invoked)."""
    self.lb, self.ub = lb, ub
    self.output.update_bounds(lb=lb, ub=ub)

LinearProgram

Bases: QuadraticProgram

Step representing the solution of a linear program

\[ x^\star = \arg\min_x\; q^\top x \quad\text{s.t.}\quad A x = b,\; C x \le d. \]

A QuadraticProgram specialization with P = 0 (stored as a sparse zero matrix). The KKT system drops the quadratic term: stationarity becomes the dual-only equality \(q + A^\top \nu + C^\top \lambda = 0\), so the primal x is pinned by primal feasibility + complementarity rather than by stationarity. Because stationarity no longer bounds x, finite primal bounds must come from primal_bounds or OBBT over {Ax = b, Cx <= d}.

Parameters:

Name Type Description Default
q ndarray or LinearExpression

Linear cost; its dimension defines n. Required (no -arg default makes sense for an LP cost). A scalar q is rejected (cannot size the LP).

 required
A

See QuadraticProgram.

 None
b

See QuadraticProgram.

 None
C

See QuadraticProgram.

 None
d

See QuadraticProgram.

 None
dual_ub

See QuadraticProgram.

 None
primal_bounds

See QuadraticProgram.

 None
nu_bounds

See QuadraticProgram.

 None
arg

See QuadraticProgram.

 None
relax_binary_vars

See QuadraticProgram.

 None
Source code in fomtools/steps/lin_prog.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
class LinearProgram(QuadraticProgram):
    r"""
    Step representing the solution of a linear program

    $$
        x^\star = \arg\min_x\; q^\top x
        \quad\text{s.t.}\quad A x = b,\; C x \le d.
    $$

    A `QuadraticProgram` specialization with ``P = 0`` (stored as a sparse
    zero matrix). The KKT system drops the quadratic term: stationarity becomes
    the dual-only equality $q + A^\top \nu + C^\top \lambda = 0$, so the
    primal ``x`` is pinned by primal feasibility + complementarity rather than by
    stationarity. Because stationarity no longer bounds ``x``, finite primal
    bounds must come from ``primal_bounds`` or OBBT over ``{Ax = b, Cx <= d}``.

    Parameters
    ----------
    q : np.ndarray or LinearExpression
        Linear cost; its dimension defines ``n``. Required (no ``-arg`` default
        makes sense for an LP cost). A scalar ``q`` is rejected (cannot size the LP).
    A, b, C, d, dual_ub, primal_bounds, nu_bounds, arg, relax_binary_vars
        See `QuadraticProgram`.
    """

    def __init__(
        self,
        q,
        A=None,
        b=None,
        C=None,
        d=None,
        dual_ub=None,
        primal_bounds=None,
        nu_bounds=None,
        arg: LinearExpression = None,
        relax_binary_vars: bool = False,
    ):
        if _is_linexpr(q):
            n = q.ndim
        else:
            q_arr = np.asarray(q, dtype=float)
            if q_arr.ndim == 0:
                raise ValueError("LinearProgram requires a vector q (a scalar cannot size the LP)")
            n = q_arr.shape[0]

        P = spa.csr_matrix((n, n))  # P = 0
        super().__init__(
            P=P,
            q=q,
            A=A,
            b=b,
            C=C,
            d=d,
            dual_ub=dual_ub,
            primal_bounds=primal_bounds,
            nu_bounds=nu_bounds,
            arg=arg,
            relax_binary_vars=relax_binary_vars,
        )
        self.is_lp = True

StandardFormLinearProgram

Bases: LinearProgram

Step representing the solution of a standard-form linear program

\[ x^\star = \arg\min_x\; c^\top x \quad\text{s.t.}\quad A x = b,\; x \ge 0. \]

A LinearProgram with q = c, equality (A, b), and the non-negativity constraint x >= 0 written as the inequality -I x <= 0 (C = -I, d = 0). No new lowering logic — it reuses the entire LinearProgram / _qp_kkt machinery (the step decomposer resolves it to the LinearProgram decomposer via the class MRO).

Parameters:

Name Type Description Default
c ndarray or LinearExpression

Linear cost; its dimension defines n. A scalar is rejected (cannot size the LP).

 required
A optional

Equality matrix (p, n) and right-hand side (dim p).

 None
b optional

Equality matrix (p, n) and right-hand side (dim p).

 None
dual_ub

See QuadraticProgram. dual_ub is required (the x >= 0 inequality is always present, so the decomposer needs a finite bound on its multiplier).

 None
primal_bounds

See QuadraticProgram. dual_ub is required (the x >= 0 inequality is always present, so the decomposer needs a finite bound on its multiplier).

 None
nu_bounds

See QuadraticProgram. dual_ub is required (the x >= 0 inequality is always present, so the decomposer needs a finite bound on its multiplier).

 None
arg

See QuadraticProgram. dual_ub is required (the x >= 0 inequality is always present, so the decomposer needs a finite bound on its multiplier).

 None
relax_binary_vars

See QuadraticProgram. dual_ub is required (the x >= 0 inequality is always present, so the decomposer needs a finite bound on its multiplier).

 None
Source code in fomtools/steps/lin_prog.py 
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
class StandardFormLinearProgram(LinearProgram):
    r"""
    Step representing the solution of a **standard-form** linear program

    $$
        x^\star = \arg\min_x\; c^\top x
        \quad\text{s.t.}\quad A x = b,\; x \ge 0.
    $$

    A `LinearProgram` with ``q = c``, equality ``(A, b)``, and the non-negativity
    constraint ``x >= 0`` written as the inequality ``-I x <= 0`` (``C = -I``, ``d = 0``).
    No new lowering logic — it reuses the entire ``LinearProgram`` / ``_qp_kkt`` machinery
    (the step decomposer resolves it to the ``LinearProgram`` decomposer via the class MRO).

    Parameters
    ----------
    c : np.ndarray or LinearExpression
        Linear cost; its dimension defines ``n``. A scalar is rejected (cannot size the LP).
    A, b : optional
        Equality matrix ``(p, n)`` and right-hand side (dim ``p``).
    dual_ub, primal_bounds, nu_bounds, arg, relax_binary_vars
        See `QuadraticProgram`. ``dual_ub`` is required (the ``x >= 0`` inequality is
        always present, so the decomposer needs a finite bound on its multiplier).
    """

    def __init__(
        self,
        c,
        A=None,
        b=None,
        dual_ub=None,
        primal_bounds=None,
        nu_bounds=None,
        arg: LinearExpression = None,
        relax_binary_vars: bool = False,
    ):
        if _is_linexpr(c):
            n = c.ndim
        else:
            c_arr = np.asarray(c, dtype=float)
            if c_arr.ndim == 0:
                raise ValueError(
                    "StandardFormLinearProgram requires a vector c (a scalar cannot size the LP)"
                )
            n = c_arr.shape[0]

        super().__init__(
            q=c,
            A=A,
            b=b,
            C=-np.eye(n),  # x >= 0  <=>  -I x <= 0
            d=np.zeros(n),
            dual_ub=dual_ub,
            primal_bounds=primal_bounds,
            nu_bounds=nu_bounds,
            arg=arg,
            relax_binary_vars=relax_binary_vars,
        )