Skip to content

Solvers

All solver classes share the setup / solve / update workflow and use the same Settings. See Re-solving with new data for the fixed-structure update pattern and Differentiation for the backward() workflow. They differ only in the accepted storage format for P, A, G and the KKT factorization used. See Backends for guidance on choosing one.

DenseSolver

DenseSolver 

DenseSolver(
    dtype: Literal["float32", "float64"] = "float64",
)

Bases: SolverBase

GPU solver for general dense convex quadratic programs that solves a QP - or a whole batch of QPs - of the form

\[ \begin{aligned} \min_{x}\quad & \tfrac{1}{2}\, x^\top P x + c^\top x \\ \text{s.t.}\quad & A x = b, \\ & h_l \le G x \le h_u, \\ & x_l \le x \le x_u, \end{aligned} \]

using the proximal interior-point method, running entirely on the GPU.

Inputs. P, A, G and every vector (c, b, h_l, h_u, x_l, x_u) must be dense arrays that live on the GPU. Any object exposing the __cuda_array_interface__ protocol is accepted - a cupy.ndarray, a CUDA torch.Tensor, a CUDA JAX array, a Numba device array, and so on.

cuPIQP is GPU-only: CPU data (numpy.ndarray, CPU torch tensors, CPU JAX arrays) is rejected with a TypeError rather than copied to the device behind your back.

Batching. DenseSolver is natively batched: solve B independent QPs in a single GPU call by giving every array a leading batch dimension - P of shape (B, n, n), c of shape (B, n), and so on. A single problem is simply B = 1, and solver.result.x then has shape (B, n).

Parameters:

Name Type Description Default
dtype (float64, float32)

Floating-point precision used throughout the solve. "float32" is faster and uses less memory but converges to looser tolerances; the default convergence tolerances are chosen to match the dtype.

 "float64"

Examples:

A small inequality-constrained QP (one row is one-sided via -inf):

import cupy as cp
from cupiqp import DenseSolver

P = cp.eye(2)
c = cp.array([-1.0, -4.0])
G = cp.array([[1.0, 1.0]])      # constrain x1 + x2
h_l = cp.array([-cp.inf])       # no lower bound on the row
h_u = cp.array([1.0])           # x1 + x2 <= 1

solver = DenseSolver()
solver.setup(P=P, c=c, G=G, h_l=h_l, h_u=h_u)
solver.solve()

print(solver.result.info.status[0].name)    # CUPIQP_SOLVED
x = solver.result.x.get()[0]                 # bring the solution to the host
See Also

SparseSolver: solver for general sparse problems.

MultistageSolver: structure-exploiting solver for multistage optimization (e.g. optimal-control) problems.

Notes

The problem structure - array shapes, which constraint blocks are present, and which bounds are finite vs. +/-inf - is fixed by setup and can only be set up once per instance. To re-solve with new numerical values of the same structure (e.g. a moving target b in receding-horizon control), call update and then solve again, which reuses all GPU allocations; for a different structure, create a new DenseSolver. Solver behaviour (tolerances, verbosity, iteration cap, ...) is configured through solver.settings.

Source code in cupiqp/dense/dense_solver.py 
114
115
116
def __init__(self, dtype: Literal["float32", "float64"] = "float64"):
    super().__init__(dtype=dtype)
    self._settings.kkt_solver = "dense_cholesky"

setup 

setup(
    P: CudaArray,
    c: CudaArray,
    A: Optional[CudaArray] = None,
    b: Optional[CudaArray] = None,
    G: Optional[CudaArray] = None,
    h_u: Optional[CudaArray] = None,
    h_l: Optional[CudaArray] = None,
    x_u: Optional[CudaArray] = None,
    x_l: Optional[CudaArray] = None,
) -> None

Bind the problem data and prepare the solver for solve().

Fixes the problem structure - array shapes, which constraint blocks are present, the sparsity pattern (sparse backend), and the finite/infinite pattern of the bounds - and allocates all GPU buffers, the KKT system, and the preconditioner. Call this once per solver instance, then call solve().

Pass a single problem (2D P) or a batch (3D P with a leading batch axis, or a list of matrices); the batch size is inferred here and sets the shape of the result.

Parameters:

Name Type Description Default
P GPU array

Quadratic cost, shape (n, n) or batched (B, n, n). Must be symmetric positive semidefinite. Required.

 required
c GPU array

Linear cost, shape (n,) or (B, n). Required.

 required
A GPU array

Equality constraints A x = b; shapes (p, n) and (p,) (or batched). Omit for no equality constraints.

 None
b GPU array

Equality constraints A x = b; shapes (p, n) and (p,) (or batched). Omit for no equality constraints.

 None
G GPU array

Two-sided inequalities h_l <= G x <= h_u; G is (m, n) (or batched) and the bounds are (m,). Use -inf / +inf entries for one-sided rows.

 None
h_l GPU array

Two-sided inequalities h_l <= G x <= h_u; G is (m, n) (or batched) and the bounds are (m,). Use -inf / +inf entries for one-sided rows.

 None
h_u GPU array

Two-sided inequalities h_l <= G x <= h_u; G is (m, n) (or batched) and the bounds are (m,). Use -inf / +inf entries for one-sided rows.

 None
x_l GPU array

Element-wise box bounds x_l <= x <= x_u, shape (n,) (or batched). Use +/-inf for unbounded entries.

 None
x_u GPU array

Element-wise box bounds x_l <= x <= x_u, shape (n,) (or batched). Use +/-inf for unbounded entries.

 None

Raises:

Type Description
RuntimeError

If setup() has already been called on this instance. The structure is fixed after setup - create a new solver for a different structure, or use update() to change only the numerical values.
TypeError

If an input is not a GPU array of the kind this backend expects (e.g. a CPU numpy array, or a dense matrix passed to the sparse backend). See the backend's class docstring for the exact accepted types.
See Also

solve : run the solver after setup. update : change numerical data without a full re-setup.

Source code in cupiqp/dense/dense_solver.py 
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
def setup(
    self,
    P: CudaArray,
    c: CudaArray,
    A: Optional[CudaArray] = None,
    b: Optional[CudaArray] = None,
    G: Optional[CudaArray] = None,
    h_u: Optional[CudaArray] = None,
    h_l: Optional[CudaArray] = None,
    x_u: Optional[CudaArray] = None,
    x_l: Optional[CudaArray] = None,
) -> None:
    super().setup(P, c, A, b, G, h_u, h_l, x_u, x_l)
    if self.settings.enable_grad:
        d = self._data
        B = d.batch_size
        dtype = d.dtype
        self._dense_data_gradients_kernel = create_dense_data_gradients_kernel(
            d.n, d.p, d.m, d.num_hu, d.num_xu, dtype=dtype)
        self._grad_data = DenseData(dtype=dtype, device=self.settings.device)
        self._grad_data.init(
            P=cp.zeros((B, d.n, d.n), dtype=dtype),
            c=cp.zeros((B, d.n), dtype=dtype),
            A=cp.zeros((B, d.p, d.n), dtype=dtype) if d.p > 0 else None,
            b=cp.zeros((B, d.p), dtype=dtype) if d.p > 0 else None,
            G=cp.zeros((B, d.m, d.n), dtype=dtype) if d.m > 0 else None,
            h_u=cp.zeros((B, d.m), dtype=dtype) if d.num_hu > 0 else None,
            h_l=cp.zeros((B, d.m), dtype=dtype) if d.num_hl > 0 else None,
            x_u=cp.zeros((B, d.n), dtype=dtype) if d.num_xu > 0 else None,
            x_l=cp.zeros((B, d.n), dtype=dtype) if d.num_xl > 0 else None,
        )

solve 

solve() -> List[Status]

Solve the QP set up by setup() and return the solve status.

Runs the proximal interior-point iterations on the GPU. The full solution (primal x, dual, and slack variables) and per-problem diagnostics are written to solver.result; this method returns the status for convenience.

Returns:

Type Description
Status or list of Status

For a single problem, the Status. For a batched setup(), a list of B of them (one per problem). CUPIQP_SOLVED means the problem converged to tolerance. The list form is always available as solver.result.info.status.
Notes

Read the solution from solver.result after solving - e.g. solver.result.x (shape (B, n)) and solver.result.info.status. Set solver.settings.verbose = True to print a per-iteration log. After setup() you may solve() repeatedly, optionally calling update() in between to change the numerical data.

Source code in cupiqp/solver.py 
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
def solve(self) -> List[Status]:
    """Solve the QP set up by ``setup()`` and return the solve status.

    Runs the proximal interior-point iterations on the GPU. The full
    solution (primal ``x``, dual, and slack variables) and per-problem
    diagnostics are written to ``solver.result``; this method returns the
    status for convenience.

    Returns
    -------
    Status or list of Status
        For a single problem, the ``Status``. For a batched ``setup()``,
        a list of ``B`` of them (one per problem). ``CUPIQP_SOLVED``
        means the problem converged to tolerance. The list form is always
        available as ``solver.result.info.status``.

    Notes
    -----
    Read the solution from ``solver.result`` after solving - e.g.
    ``solver.result.x`` (shape ``(B, n)``) and
    ``solver.result.info.status``. Set ``solver.settings.verbose = True``
    to print a per-iteration log. After ``setup()`` you may ``solve()``
    repeatedly, optionally calling ``update()`` in between to change the
    numerical data.
    """
    if self.settings.verbose:
        try:
            from importlib.metadata import version
            _ver = version("cupiqp")
        except Exception:
            _ver = ""
        _w = 58
        print("-" * _w)
        print(f"cuPIQP v{_ver} - GPU-accelerated PIQP solver".strip().center(_w))
        print("(c) Fenglong Song".center(_w))
        print("Ecole Polytechnique Federale de Lausanne (EPFL) 2026".center(_w))
        print("-" * _w)
        if self.settings.kkt_solver == "dense_cholesky":
            print("dense backend:")
            print(f"batch size B = {self._data.batch_size}")
            print(f"variables n = {self._data.n}")
            print(f"equality constraints p = {self._data.p}")
            print(f"inequality constraints m = {self._data.m}")
        elif self.settings.kkt_solver == "sparse_ldlt":
            print("sparse backend:")
            print(f"batch size B = {self._data.batch_size}")
            print(f"variables n = {self._data.n}, nnz(P) = {self._data.P.nnz}")
            print(f"equality constraints p = {self._data.p}, nnz(A) = {self._data.A.nnz}")
            print(f"inequality constraints m = {self._data.m}, nnz(G) = {self._data.G.nnz}")
        elif self.settings.kkt_solver == "multistage_block_cholesky":
            print("multistage backend:")
            print(f"batch size B = {self._data.batch_size}")
            print(f"variables n = {self._data.n}, num_diag_blocks(P) = {self._data.P.num_diag_blocks}, block_size(P) = ({self._data.P.block_size}, {self._data.P.block_size})")
            print(f"equality constraints p = {self._data.p}, num_diag_blocks(A) = {self._data.A.N}, block_size(A) = ({self._data.A.rows_of_blocks}, {self._data.A.cols_of_blocks})")
            print(f"inequality constraints m = {self._data.m}, num_diag_blocks(G) = {self._data.G.N}, block_size(G) = ({self._data.G.rows_of_blocks}, {self._data.G.cols_of_blocks})")
        else:
            raise ValueError(f"Unsupported kkt_solver type: {self.settings.kkt_solver}")

        print(f"inequality lower bounds n_h_l = {self._data.num_hl}")
        print(f"inequality upper bounds n_h_u = {self._data.num_hu}")
        print(f"variable lower bounds n_x_l = {self._data.num_xl}")
        print(f"variable upper bounds n_x_u = {self._data.num_xu}")
        print("")
    return self._solve_impl()

update 

update(
    P: Optional[Any] = None,
    c: Optional[Any] = None,
    A: Optional[Any] = None,
    b: Optional[Any] = None,
    G: Optional[Any] = None,
    h_u: Optional[Any] = None,
    h_l: Optional[Any] = None,
    x_u: Optional[Any] = None,
    x_l: Optional[Any] = None,
    check_validity: bool = False,
)

Change the numerical problem data, then solve() again.

The fast path for re-solving a problem of the same structure - for example a moving target b or a re-linearized P in receding-horizon control. It reuses every GPU allocation from setup(), so only the values change; shapes, sparsity patterns, and which blocks are present must stay the same (create a new solver for a structural change). Bound values may change freely, including which entries are +/-inf - a bound can flip between finite and infinite without re-setup().

Any argument left as None keeps its current value. After update(), call solve() to get the new solution.

Parameters:

Name Type Description Default
P GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
c GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
A GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
b GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
G GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
h_u GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
h_l GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
x_u GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
x_l GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
check_validity bool

If True, validate the dimensions and sparsity of the new data. Defaults to False for speed (validation forces device-to-host syncs in the sparse backend). When False, you must still keep the shapes and sparsity patterns of P/A/ G unchanged; bound values (including which entries are +/-inf) may change.

 False
Source code in cupiqp/solver.py 
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
def update(self,
           P: Optional[Any] = None,
           c: Optional[Any] = None,
           A: Optional[Any] = None,
           b: Optional[Any] = None,
           G: Optional[Any] = None,
           h_u: Optional[Any] = None,
           h_l: Optional[Any] = None,
           x_u: Optional[Any] = None,
           x_l: Optional[Any] = None,
           check_validity: bool = False,
           ):
    """Change the numerical problem data, then ``solve()`` again.

    The fast path for re-solving a problem of the **same structure** -
    for example a moving target ``b`` or a re-linearized ``P`` in
    receding-horizon control. It reuses every GPU allocation from
    ``setup()``, so only the values change; shapes, sparsity patterns,
    and which blocks are present must stay the same (create a new solver
    for a structural change). Bound *values* may change freely, including
    which entries are ``+/-inf`` - a bound can flip between finite and
    infinite without re-``setup()``.

    Any argument left as ``None`` keeps its current value. After
    ``update()``, call ``solve()`` to get the new solution.

    Parameters
    ----------
    P, c, A, b, G, h_u, h_l, x_u, x_l : GPU array, optional
        New values for the corresponding problem block. ``None`` (the
        default) leaves that block unchanged. Must match the original
        shapes / sparsity pattern set at ``setup()``.
    check_validity : bool, default: False
        If ``True``, validate the dimensions and sparsity of the new
        data. Defaults to ``False`` for speed (validation forces
        device-to-host syncs in the sparse backend). When ``False``, you
        must still keep the shapes and sparsity patterns of ``P``/``A``/
        ``G`` unchanged; bound values (including which entries are
        ``+/-inf``) may change.
    """
    if not self._setup_done:
        raise RuntimeError("Solver not setup yet. Call setup() first.")

    if self.settings.preconditioner_iter > 0:
        self._preconditioner.unscale_data(self._data)

    if P is not None:
        self._data.set_P(P, check=check_validity)
    if c is not None:
        self._data.set_c(c, check=check_validity)
    if A is not None:
        self._data.set_A(A, check=check_validity)
    if b is not None:
        self._data.set_b(b, check=check_validity)
    if G is not None:
        self._data.set_G(G, check=check_validity)
    if h_u is not None:
        self._data.set_h_u(h_u, check=check_validity)
    if h_l is not None:
        self._data.set_h_l(h_l, check=check_validity)
    if x_u is not None:
        self._data.set_x_u(x_u, check=check_validity)
    if x_l is not None:
        self._data.set_x_l(x_l, check=check_validity)

    matrix_changed = P is not None or A is not None or G is not None

    # NOTE: Since we allow changing h_l/h_u containing arbitrary +inf/-inf, 
    # an inequality row G[i] can switch between active (a finite
    # bound) and inactive (both bounds infinite) between updates.
    # If either of h_l or h_u are updated, we need to update G 
    # because for sparse kkt solver we need to set the inactive rows to 0
    ineq_bound_pattern_may_change = h_l is not None or h_u is not None

    # Apply preconditioner scaling to updated data.
    preconditioner_did_fresh_ruiz = False
    if self.settings.preconditioner_iter > 0:
        reuse = self.settings.preconditioner_reuse_on_update or not matrix_changed
        if reuse:
            self._preconditioner.reuse_scaling(self._data)
        else:
            self._preconditioner.reset()
            self._preconditioner.scale_data(
                self._data,
                self.settings.preconditioner_scale_cost,
                self.settings.preconditioner_iter,
            )
            preconditioner_did_fresh_ruiz = True

    self._preconditioner.compute_constraints_rhs_inf_norm_unscaled(
        self._data, self._constraints_rhs_inf_norm_unscaled,
    )
    # Fresh Ruiz produces new factors that re-scale ALL of P/A/G in place,
    # even matrices the user didn't pass. The KKT solver caches things
    # like A^T A keyed off those scaled values, so flag everything as
    # changed in that case.
    self._kkt_system.update_data(
        self._data,
        (P is not None) or preconditioner_did_fresh_ruiz,
        (A is not None) or preconditioner_did_fresh_ruiz,
        (G is not None) or preconditioner_did_fresh_ruiz or ineq_bound_pattern_may_change,
    )

backward 

backward(
    grad_x=None,
    grad_y=None,
    grad_z_u=None,
    grad_z_l=None,
    grad_z_bu=None,
    grad_z_bl=None,
    grad_s_u=None,
    grad_s_l=None,
    grad_s_bu=None,
    grad_s_bl=None,
)

Compute gradients of an outer scalar :math:L w.r.t. problem data, given upstream cotangents on the solution variables.

Orchestration (backend-agnostic):

  1. Pack the per-field cotangent kwargs into self._grad_in (a pre-allocated Variables); missing kwargs are treated as zeros.
  2. Solve the adjoint KKT system via :meth:_compute_adjoint, producing user-space adjoint vectors.
  3. Scatter the four active-size lambda groups (z_u, z_l, z_bu, z_bl) and the two active-size ineq result groups into full-m / full-n buffers (self._lam_z*_full, self._z*_full). Both the dG outer product and the dh_* / dx_* vector gradients consume these full-layout buffers.
  4. Delegate to :meth:_compute_data_gradients for backend- specific matrix-gradient assembly + Data subclass construction.

Returns the backend's Data subclass populated with the gradients in user space. Cotangents and returned gradients are interpreted in user (un-scaled) space throughout — the adjoint solve and scatter chain handles all preconditioner bookkeeping internally.

Raises:

Type Description
RuntimeError

If :meth:solve has not been called yet (no cached KKT factor).
Source code in cupiqp/solver.py 
1828
1829
1830
1831
1832
1833
1834
1835
1836
1837
1838
1839
1840
1841
1842
1843
1844
1845
1846
1847
1848
1849
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880
1881
1882
1883
1884
1885
1886
1887
1888
1889
1890
1891
1892
1893
1894
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
1907
1908
1909
1910
1911
1912
1913
1914
1915
1916
1917
1918
1919
1920
1921
1922
1923
1924
1925
1926
@nvtx.annotate("Solver::grad")
def backward(self,
         grad_x=None, grad_y=None,
         grad_z_u=None, grad_z_l=None, grad_z_bu=None, grad_z_bl=None,
         grad_s_u=None, grad_s_l=None, grad_s_bu=None, grad_s_bl=None):
    r"""Compute gradients of an outer scalar :math:`L` w.r.t. problem
    data, given upstream cotangents on the solution variables.

    Orchestration (backend-agnostic):

    1. Pack the per-field cotangent kwargs into ``self._grad_in``
       (a pre-allocated ``Variables``); missing kwargs are treated
       as zeros.
    2. Solve the adjoint KKT system via :meth:`_compute_adjoint`,
       producing user-space adjoint vectors.
    3. Scatter the four active-size lambda groups (``z_u, z_l,
       z_bu, z_bl``) and the two active-size ineq result groups
       into full-``m`` / full-``n`` buffers
       (``self._lam_z*_full``, ``self._z*_full``). Both the dG
       outer product and the ``dh_*`` / ``dx_*`` vector gradients
       consume these full-layout buffers.
    4. Delegate to :meth:`_compute_data_gradients` for backend-
       specific matrix-gradient assembly + ``Data`` subclass
       construction.

    Returns the backend's ``Data`` subclass populated with the
    gradients in user space. Cotangents and returned gradients are
    interpreted in user (un-scaled) space throughout — the
    adjoint solve and scatter chain handles all preconditioner
    bookkeeping internally.

    Raises
    ------
    RuntimeError
        If :meth:`solve` has not been called yet (no cached KKT factor).
    """
    if not self.settings.enable_grad:
        raise RuntimeError("Set enable_grad to True to enable gradient computation.")

    if not getattr(self, "_setup_done", False) or self._result is None:
        raise RuntimeError(
            f"{type(self).__name__}.grad() requires a prior solve(); "
            f"call setup() and solve() before grad()."
        )

    data = self._data
    B = data.batch_size
    wp_stream = wp.Stream(cuda_stream=cp.cuda.get_current_stream().ptr)

    # ---- Step 1
    zeros = self._zero_grad_in
    pack_total = data.n + data.p + 2 * zeros.num_ineq
    if pack_total > 0:
        wp.launch(
            kernel=self._backward_copy_kernel,
            dim=(B, pack_total),
            inputs=[
                grad_x    if grad_x    is not None else zeros.x,
                grad_y    if grad_y    is not None else zeros.y,
                grad_z_u  if grad_z_u  is not None else zeros.z_u,
                grad_z_l  if grad_z_l  is not None else zeros.z_l,
                grad_z_bu if grad_z_bu is not None else zeros.z_bu,
                grad_z_bl if grad_z_bl is not None else zeros.z_bl,
                grad_s_u  if grad_s_u  is not None else zeros.s_u,
                grad_s_l  if grad_s_l  is not None else zeros.s_l,
                grad_s_bu if grad_s_bu is not None else zeros.s_bu,
                grad_s_bl if grad_s_bl is not None else zeros.s_bl,
                self._grad_in.x,    self._grad_in.y,
                self._grad_in.z_u,  self._grad_in.z_l,
                self._grad_in.z_bu, self._grad_in.z_bl,
                self._grad_in.s_u,  self._grad_in.s_l,
                self._grad_in.s_bu, self._grad_in.s_bl,
            ],
            device="cuda",
            stream=wp_stream,
        )

    # ---- Step 2: adjoint KKT solve
    self._compute_adjoint(self._grad_in, self._backward_adjoint_vector)

    # ---- Step 3: write into full-layout buffers
    wp.launch(
        kernel=self._backward_pack_full_layout_kernel,
        dim=(B, 4 * data.m + data.num_xu + data.num_xl),
        inputs=[
            self._backward_adjoint_vector.z_u, self._backward_adjoint_vector.z_l,
            self._backward_adjoint_vector.z_bu, self._backward_adjoint_vector.z_bl,
            self._result.z_u, self._result.z_l,
            self._lam_zu_full, self._lam_zl_full,
            self._lam_zbu_full, self._lam_zbl_full,
            self._zu_full, self._zl_full,
        ],
        device="cuda",
        stream=wp_stream,
    )

    # ---- Step 4: backend-specific matrix and vector gradient
    # assembly + Data subclass construction.
    return self._compute_data_gradients(self._backward_adjoint_vector)

SparseSolver

SparseSolver 

SparseSolver(
    dtype: Literal["float32", "float64"] = "float64",
)

Bases: SolverBase

GPU solver for general sparse convex quadratic programs that solves a QP - or a whole batch of QPs - of the form

\[ \begin{aligned} \min_{x}\quad & \tfrac{1}{2}\, x^\top P x + c^\top x \\ \text{s.t.}\quad & A x = b, \\ & h_l \le G x \le h_u, \\ & x_l \le x \le x_u, \end{aligned} \]

using the proximal interior-point method with a sparse LDL^T factorization, running entirely on the GPU. It is built for large, structurally sparse P / A / G.

Inputs - CSR matrices, dense vectors. P, A, G must be GPU sparse matrices in CSR layout; the vectors (c, b, h_l, h_u, x_l, x_u) are dense GPU arrays. Accepted matrix types are:

  • cupyx.scipy.sparse.csr_matrix (a GPU CSR matrix),
  • a CUDA torch.sparse_csr_tensor, or
  • a UniformBatchedCsrMatrix - cuPIQP's own container holding a batch of CSR matrices that share one sparsity pattern (the efficient way to pass a batch).

A list / tuple of cupyx.scipy.sparse.csr_matrix (one per batch element, all sharing the same sparsity pattern) is also accepted, but discouraged: separate matrix objects cannot be organized with the uniform stride that batched linear-algebra routines require, so cuPIQP must copy them into a single contiguous UniformBatchedCsrMatrix at setup. Build and pass one of those yourself to avoid the copy.

cuPIQP is GPU-only and CSR-only, and never converts formats behind your back:

  • CPU sparse matrices (scipy.sparse.*, a CPU torch.sparse_csr_tensor) are rejected - lift them onto the GPU first, e.g. cupyx.scipy.sparse.csr_matrix(P_scipy).
  • Non-CSR GPU layouts (CSC, BSR, BSC, COO) are rejected - convert with .tocsr() before passing.

Batching. Solve B independent QPs in a single GPU call by passing P, A, G as a UniformBatchedCsrMatrix (the preferred, fastest batched input), with the vectors stacked along a leading batch dimension (c of shape (B, n), etc.). A single problem is simply B = 1, and solver.result.x then has shape (B, n). (A list of per-element CSR matrices also works but is slower - see above.)

Parameters:

Name Type Description Default
dtype (float64, float32)

Floating-point precision used throughout the solve. "float32" is faster and uses less memory but converges to looser tolerances; the default convergence tolerances are chosen to match the dtype.

 "float64"

Examples:

A sparse problem assembled on the host with SciPy, then lifted to the GPU:

import scipy.sparse as sp
import cupy as cp
from cupyx.scipy.sparse import csr_matrix
from cupiqp import SparseSolver

P = csr_matrix(sp.eye(4, format="csr"))     # lift scipy -> GPU CSR
c = cp.zeros(4)

solver = SparseSolver()
solver.setup(P=P, c=c)
solver.solve()

print(solver.result.info.status[0].name)    # CUPIQP_SOLVED
See Also

DenseSolver: solver for general dense problems.

MultistageSolver: structure-exploiting solver for multistage optimization (e.g. optimal-control) problems.

Notes

The problem structure - array shapes, the sparsity pattern of each matrix, which constraint blocks are present, and which bounds are finite vs. +/-inf - is fixed by setup and can only be set up once per instance. To re-solve after changing only the numerical values (keeping the same sparsity pattern), call update and then solve again, which reuses all GPU allocations; for a different structure, create a new SparseSolver. Solver behaviour (tolerances, verbosity, iteration cap, ...) is configured through solver.settings.

Source code in cupiqp/sparse/sparse_solver.py 
200
201
202
def __init__(self, dtype: Literal["float32", "float64"] = "float64"):
    super().__init__(dtype=dtype)
    self._settings.kkt_solver = "sparse_ldlt"

setup 

setup(
    P: CsrMatrixInput,
    c: CudaArray,
    A: Optional[CsrMatrixInput] = None,
    b: Optional[CudaArray] = None,
    G: Optional[CsrMatrixInput] = None,
    h_u: Optional[CudaArray] = None,
    h_l: Optional[CudaArray] = None,
    x_u: Optional[CudaArray] = None,
    x_l: Optional[CudaArray] = None,
) -> None

Bind the problem data and prepare the solver for solve().

Fixes the problem structure - array shapes, which constraint blocks are present, the sparsity pattern (sparse backend), and the finite/infinite pattern of the bounds - and allocates all GPU buffers, the KKT system, and the preconditioner. Call this once per solver instance, then call solve().

Pass a single problem (2D P) or a batch (3D P with a leading batch axis, or a list of matrices); the batch size is inferred here and sets the shape of the result.

Parameters:

Name Type Description Default
P GPU array

Quadratic cost, shape (n, n) or batched (B, n, n). Must be symmetric positive semidefinite. Required.

 required
c GPU array

Linear cost, shape (n,) or (B, n). Required.

 required
A GPU array

Equality constraints A x = b; shapes (p, n) and (p,) (or batched). Omit for no equality constraints.

 None
b GPU array

Equality constraints A x = b; shapes (p, n) and (p,) (or batched). Omit for no equality constraints.

 None
G GPU array

Two-sided inequalities h_l <= G x <= h_u; G is (m, n) (or batched) and the bounds are (m,). Use -inf / +inf entries for one-sided rows.

 None
h_l GPU array

Two-sided inequalities h_l <= G x <= h_u; G is (m, n) (or batched) and the bounds are (m,). Use -inf / +inf entries for one-sided rows.

 None
h_u GPU array

Two-sided inequalities h_l <= G x <= h_u; G is (m, n) (or batched) and the bounds are (m,). Use -inf / +inf entries for one-sided rows.

 None
x_l GPU array

Element-wise box bounds x_l <= x <= x_u, shape (n,) (or batched). Use +/-inf for unbounded entries.

 None
x_u GPU array

Element-wise box bounds x_l <= x <= x_u, shape (n,) (or batched). Use +/-inf for unbounded entries.

 None

Raises:

Type Description
RuntimeError

If setup() has already been called on this instance. The structure is fixed after setup - create a new solver for a different structure, or use update() to change only the numerical values.
TypeError

If an input is not a GPU array of the kind this backend expects (e.g. a CPU numpy array, or a dense matrix passed to the sparse backend). See the backend's class docstring for the exact accepted types.
See Also

solve : run the solver after setup. update : change numerical data without a full re-setup.

Source code in cupiqp/sparse/sparse_solver.py 
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
def setup(
    self,
    P: CsrMatrixInput,
    c: CudaArray,
    A: Optional[CsrMatrixInput] = None,
    b: Optional[CudaArray] = None,
    G: Optional[CsrMatrixInput] = None,
    h_u: Optional[CudaArray] = None,
    h_l: Optional[CudaArray] = None,
    x_u: Optional[CudaArray] = None,
    x_l: Optional[CudaArray] = None,
) -> None:
    super().setup(P, c, A, b, G, h_u, h_l, x_u, x_l)
    # Cache CSR row decompressions for the backward-pass gather.
    # Sparsity patterns are fixed at setup, so this is done once.
    if self.settings.enable_grad:
        d = self._data
        B = d.batch_size
        dtype = d.dtype

        # Row indices for each nnz position (CSR-to-COO). All row/col
        # index arrays are cast to int32 to match the warp kernel's
        # index type (consistent with idx_hu etc.).
        P_csr = d._P
        nnz_P = int(P_csr.nnz)
        self._p_rows = (cp.searchsorted(
            P_csr.indptr,
            cp.arange(nnz_P, dtype=P_csr.indptr.dtype),
            side="right",
        ) - 1).astype(cp.int32)
        self._p_indices_arr = P_csr.indices.astype(cp.int32)

        if d.p > 0:
            A_csr = d._A
            nnz_A = int(A_csr.nnz)
            self._a_rows = (cp.searchsorted(
                A_csr.indptr,
                cp.arange(nnz_A, dtype=A_csr.indptr.dtype),
                side="right",
            ) - 1).astype(cp.int32)
            self._a_indices_arr = A_csr.indices.astype(cp.int32)
        else:
            nnz_A = 0
            self._a_rows = cp.empty(0, dtype=cp.int32)
            self._a_indices_arr = cp.empty(0, dtype=cp.int32)

        if d.m > 0:
            G_csr = d._G
            nnz_G = int(G_csr.nnz)
            self._g_rows = (cp.searchsorted(
                G_csr.indptr,
                cp.arange(nnz_G, dtype=G_csr.indptr.dtype),
                side="right",
            ) - 1).astype(cp.int32)
            self._g_indices_arr = G_csr.indices.astype(cp.int32)
        else:
            nnz_G = 0
            self._g_rows = cp.empty(0, dtype=cp.int32)
            self._g_indices_arr = cp.empty(0, dtype=cp.int32)

        # Eager-compile the fused sparse data-gradients kernel.
        self._sparse_data_gradients_kernel = create_sparse_data_gradients_kernel(
            nnz_P, nnz_A, nnz_G, d.p, d.m, d.n, d.num_hu, d.num_xu, dtype=dtype)

        # Pre-allocate the gradient SparseData. The matrix
        # UniformBatchedCsrMatrix views share the forward sparsity (same
        # indices/indptr); their values buffers become the kernel-
        # output targets. Vector grads (c, h_l, x_l) are filled via
        # slice-assign in :meth:`_compute_data_gradients`.
        P_grad_csr = UniformBatchedCsrMatrix(
            B, P_csr.indices, P_csr.indptr, cp.zeros((B, nnz_P), dtype=dtype),
            shape=(P_csr.rows, P_csr.cols), dtype=dtype,
        )
        A_grad_csr = (UniformBatchedCsrMatrix(
            B, A_csr.indices, A_csr.indptr, cp.zeros((B, nnz_A), dtype=dtype),
            shape=(A_csr.rows, A_csr.cols), dtype=dtype,
        ) if d.p > 0 else None)
        G_grad_csr = (UniformBatchedCsrMatrix(
            B, G_csr.indices, G_csr.indptr, cp.zeros((B, nnz_G), dtype=dtype),
            shape=(G_csr.rows, G_csr.cols), dtype=dtype,
        ) if d.m > 0 else None)
        self._grad_data = SparseData(dtype=dtype, device=self.settings.device)
        self._grad_data.init(
            P=P_grad_csr,
            c=cp.zeros((B, d.n), dtype=dtype),
            A=A_grad_csr,
            b=cp.zeros((B, d.p), dtype=dtype) if d.p > 0 else None,
            G=G_grad_csr,
            h_u=cp.zeros((B, d.m), dtype=dtype) if d.num_hu > 0 else None,
            h_l=cp.zeros((B, d.m), dtype=dtype) if d.num_hl > 0 else None,
            x_u=cp.zeros((B, d.n), dtype=dtype) if d.num_xu > 0 else None,
            x_l=cp.zeros((B, d.n), dtype=dtype) if d.num_xl > 0 else None,
        )
        # Kernel value-buffer inputs. SparseData._A / _G are always
        # allocated (empty BatchedCsr placeholders when the
        # corresponding block is absent), so their ``.data`` is always
        # a (B, nnz_*) array matching the compiled kernel signature.
        self._grad_P_values = self._grad_data._P.data
        self._grad_A_values = self._grad_data._A.data
        self._grad_G_values = self._grad_data._G.data

solve 

solve() -> List[Status]

Solve the QP set up by setup() and return the solve status.

Runs the proximal interior-point iterations on the GPU. The full solution (primal x, dual, and slack variables) and per-problem diagnostics are written to solver.result; this method returns the status for convenience.

Returns:

Type Description
Status or list of Status

For a single problem, the Status. For a batched setup(), a list of B of them (one per problem). CUPIQP_SOLVED means the problem converged to tolerance. The list form is always available as solver.result.info.status.
Notes

Read the solution from solver.result after solving - e.g. solver.result.x (shape (B, n)) and solver.result.info.status. Set solver.settings.verbose = True to print a per-iteration log. After setup() you may solve() repeatedly, optionally calling update() in between to change the numerical data.

Source code in cupiqp/solver.py 
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
def solve(self) -> List[Status]:
    """Solve the QP set up by ``setup()`` and return the solve status.

    Runs the proximal interior-point iterations on the GPU. The full
    solution (primal ``x``, dual, and slack variables) and per-problem
    diagnostics are written to ``solver.result``; this method returns the
    status for convenience.

    Returns
    -------
    Status or list of Status
        For a single problem, the ``Status``. For a batched ``setup()``,
        a list of ``B`` of them (one per problem). ``CUPIQP_SOLVED``
        means the problem converged to tolerance. The list form is always
        available as ``solver.result.info.status``.

    Notes
    -----
    Read the solution from ``solver.result`` after solving - e.g.
    ``solver.result.x`` (shape ``(B, n)``) and
    ``solver.result.info.status``. Set ``solver.settings.verbose = True``
    to print a per-iteration log. After ``setup()`` you may ``solve()``
    repeatedly, optionally calling ``update()`` in between to change the
    numerical data.
    """
    if self.settings.verbose:
        try:
            from importlib.metadata import version
            _ver = version("cupiqp")
        except Exception:
            _ver = ""
        _w = 58
        print("-" * _w)
        print(f"cuPIQP v{_ver} - GPU-accelerated PIQP solver".strip().center(_w))
        print("(c) Fenglong Song".center(_w))
        print("Ecole Polytechnique Federale de Lausanne (EPFL) 2026".center(_w))
        print("-" * _w)
        if self.settings.kkt_solver == "dense_cholesky":
            print("dense backend:")
            print(f"batch size B = {self._data.batch_size}")
            print(f"variables n = {self._data.n}")
            print(f"equality constraints p = {self._data.p}")
            print(f"inequality constraints m = {self._data.m}")
        elif self.settings.kkt_solver == "sparse_ldlt":
            print("sparse backend:")
            print(f"batch size B = {self._data.batch_size}")
            print(f"variables n = {self._data.n}, nnz(P) = {self._data.P.nnz}")
            print(f"equality constraints p = {self._data.p}, nnz(A) = {self._data.A.nnz}")
            print(f"inequality constraints m = {self._data.m}, nnz(G) = {self._data.G.nnz}")
        elif self.settings.kkt_solver == "multistage_block_cholesky":
            print("multistage backend:")
            print(f"batch size B = {self._data.batch_size}")
            print(f"variables n = {self._data.n}, num_diag_blocks(P) = {self._data.P.num_diag_blocks}, block_size(P) = ({self._data.P.block_size}, {self._data.P.block_size})")
            print(f"equality constraints p = {self._data.p}, num_diag_blocks(A) = {self._data.A.N}, block_size(A) = ({self._data.A.rows_of_blocks}, {self._data.A.cols_of_blocks})")
            print(f"inequality constraints m = {self._data.m}, num_diag_blocks(G) = {self._data.G.N}, block_size(G) = ({self._data.G.rows_of_blocks}, {self._data.G.cols_of_blocks})")
        else:
            raise ValueError(f"Unsupported kkt_solver type: {self.settings.kkt_solver}")

        print(f"inequality lower bounds n_h_l = {self._data.num_hl}")
        print(f"inequality upper bounds n_h_u = {self._data.num_hu}")
        print(f"variable lower bounds n_x_l = {self._data.num_xl}")
        print(f"variable upper bounds n_x_u = {self._data.num_xu}")
        print("")
    return self._solve_impl()

update 

update(
    P: Optional[Any] = None,
    c: Optional[Any] = None,
    A: Optional[Any] = None,
    b: Optional[Any] = None,
    G: Optional[Any] = None,
    h_u: Optional[Any] = None,
    h_l: Optional[Any] = None,
    x_u: Optional[Any] = None,
    x_l: Optional[Any] = None,
    check_validity: bool = False,
)

Change the numerical problem data, then solve() again.

The fast path for re-solving a problem of the same structure - for example a moving target b or a re-linearized P in receding-horizon control. It reuses every GPU allocation from setup(), so only the values change; shapes, sparsity patterns, and which blocks are present must stay the same (create a new solver for a structural change). Bound values may change freely, including which entries are +/-inf - a bound can flip between finite and infinite without re-setup().

Any argument left as None keeps its current value. After update(), call solve() to get the new solution.

Parameters:

Name Type Description Default
P GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
c GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
A GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
b GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
G GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
h_u GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
h_l GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
x_u GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
x_l GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
check_validity bool

If True, validate the dimensions and sparsity of the new data. Defaults to False for speed (validation forces device-to-host syncs in the sparse backend). When False, you must still keep the shapes and sparsity patterns of P/A/ G unchanged; bound values (including which entries are +/-inf) may change.

 False
Source code in cupiqp/solver.py 
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
def update(self,
           P: Optional[Any] = None,
           c: Optional[Any] = None,
           A: Optional[Any] = None,
           b: Optional[Any] = None,
           G: Optional[Any] = None,
           h_u: Optional[Any] = None,
           h_l: Optional[Any] = None,
           x_u: Optional[Any] = None,
           x_l: Optional[Any] = None,
           check_validity: bool = False,
           ):
    """Change the numerical problem data, then ``solve()`` again.

    The fast path for re-solving a problem of the **same structure** -
    for example a moving target ``b`` or a re-linearized ``P`` in
    receding-horizon control. It reuses every GPU allocation from
    ``setup()``, so only the values change; shapes, sparsity patterns,
    and which blocks are present must stay the same (create a new solver
    for a structural change). Bound *values* may change freely, including
    which entries are ``+/-inf`` - a bound can flip between finite and
    infinite without re-``setup()``.

    Any argument left as ``None`` keeps its current value. After
    ``update()``, call ``solve()`` to get the new solution.

    Parameters
    ----------
    P, c, A, b, G, h_u, h_l, x_u, x_l : GPU array, optional
        New values for the corresponding problem block. ``None`` (the
        default) leaves that block unchanged. Must match the original
        shapes / sparsity pattern set at ``setup()``.
    check_validity : bool, default: False
        If ``True``, validate the dimensions and sparsity of the new
        data. Defaults to ``False`` for speed (validation forces
        device-to-host syncs in the sparse backend). When ``False``, you
        must still keep the shapes and sparsity patterns of ``P``/``A``/
        ``G`` unchanged; bound values (including which entries are
        ``+/-inf``) may change.
    """
    if not self._setup_done:
        raise RuntimeError("Solver not setup yet. Call setup() first.")

    if self.settings.preconditioner_iter > 0:
        self._preconditioner.unscale_data(self._data)

    if P is not None:
        self._data.set_P(P, check=check_validity)
    if c is not None:
        self._data.set_c(c, check=check_validity)
    if A is not None:
        self._data.set_A(A, check=check_validity)
    if b is not None:
        self._data.set_b(b, check=check_validity)
    if G is not None:
        self._data.set_G(G, check=check_validity)
    if h_u is not None:
        self._data.set_h_u(h_u, check=check_validity)
    if h_l is not None:
        self._data.set_h_l(h_l, check=check_validity)
    if x_u is not None:
        self._data.set_x_u(x_u, check=check_validity)
    if x_l is not None:
        self._data.set_x_l(x_l, check=check_validity)

    matrix_changed = P is not None or A is not None or G is not None

    # NOTE: Since we allow changing h_l/h_u containing arbitrary +inf/-inf, 
    # an inequality row G[i] can switch between active (a finite
    # bound) and inactive (both bounds infinite) between updates.
    # If either of h_l or h_u are updated, we need to update G 
    # because for sparse kkt solver we need to set the inactive rows to 0
    ineq_bound_pattern_may_change = h_l is not None or h_u is not None

    # Apply preconditioner scaling to updated data.
    preconditioner_did_fresh_ruiz = False
    if self.settings.preconditioner_iter > 0:
        reuse = self.settings.preconditioner_reuse_on_update or not matrix_changed
        if reuse:
            self._preconditioner.reuse_scaling(self._data)
        else:
            self._preconditioner.reset()
            self._preconditioner.scale_data(
                self._data,
                self.settings.preconditioner_scale_cost,
                self.settings.preconditioner_iter,
            )
            preconditioner_did_fresh_ruiz = True

    self._preconditioner.compute_constraints_rhs_inf_norm_unscaled(
        self._data, self._constraints_rhs_inf_norm_unscaled,
    )
    # Fresh Ruiz produces new factors that re-scale ALL of P/A/G in place,
    # even matrices the user didn't pass. The KKT solver caches things
    # like A^T A keyed off those scaled values, so flag everything as
    # changed in that case.
    self._kkt_system.update_data(
        self._data,
        (P is not None) or preconditioner_did_fresh_ruiz,
        (A is not None) or preconditioner_did_fresh_ruiz,
        (G is not None) or preconditioner_did_fresh_ruiz or ineq_bound_pattern_may_change,
    )

backward 

backward(
    grad_x=None,
    grad_y=None,
    grad_z_u=None,
    grad_z_l=None,
    grad_z_bu=None,
    grad_z_bl=None,
    grad_s_u=None,
    grad_s_l=None,
    grad_s_bu=None,
    grad_s_bl=None,
)

Compute gradients of an outer scalar :math:L w.r.t. problem data, given upstream cotangents on the solution variables.

Orchestration (backend-agnostic):

  1. Pack the per-field cotangent kwargs into self._grad_in (a pre-allocated Variables); missing kwargs are treated as zeros.
  2. Solve the adjoint KKT system via :meth:_compute_adjoint, producing user-space adjoint vectors.
  3. Scatter the four active-size lambda groups (z_u, z_l, z_bu, z_bl) and the two active-size ineq result groups into full-m / full-n buffers (self._lam_z*_full, self._z*_full). Both the dG outer product and the dh_* / dx_* vector gradients consume these full-layout buffers.
  4. Delegate to :meth:_compute_data_gradients for backend- specific matrix-gradient assembly + Data subclass construction.

Returns the backend's Data subclass populated with the gradients in user space. Cotangents and returned gradients are interpreted in user (un-scaled) space throughout — the adjoint solve and scatter chain handles all preconditioner bookkeeping internally.

Raises:

Type Description
RuntimeError

If :meth:solve has not been called yet (no cached KKT factor).
Source code in cupiqp/solver.py 
1828
1829
1830
1831
1832
1833
1834
1835
1836
1837
1838
1839
1840
1841
1842
1843
1844
1845
1846
1847
1848
1849
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880
1881
1882
1883
1884
1885
1886
1887
1888
1889
1890
1891
1892
1893
1894
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
1907
1908
1909
1910
1911
1912
1913
1914
1915
1916
1917
1918
1919
1920
1921
1922
1923
1924
1925
1926
@nvtx.annotate("Solver::grad")
def backward(self,
         grad_x=None, grad_y=None,
         grad_z_u=None, grad_z_l=None, grad_z_bu=None, grad_z_bl=None,
         grad_s_u=None, grad_s_l=None, grad_s_bu=None, grad_s_bl=None):
    r"""Compute gradients of an outer scalar :math:`L` w.r.t. problem
    data, given upstream cotangents on the solution variables.

    Orchestration (backend-agnostic):

    1. Pack the per-field cotangent kwargs into ``self._grad_in``
       (a pre-allocated ``Variables``); missing kwargs are treated
       as zeros.
    2. Solve the adjoint KKT system via :meth:`_compute_adjoint`,
       producing user-space adjoint vectors.
    3. Scatter the four active-size lambda groups (``z_u, z_l,
       z_bu, z_bl``) and the two active-size ineq result groups
       into full-``m`` / full-``n`` buffers
       (``self._lam_z*_full``, ``self._z*_full``). Both the dG
       outer product and the ``dh_*`` / ``dx_*`` vector gradients
       consume these full-layout buffers.
    4. Delegate to :meth:`_compute_data_gradients` for backend-
       specific matrix-gradient assembly + ``Data`` subclass
       construction.

    Returns the backend's ``Data`` subclass populated with the
    gradients in user space. Cotangents and returned gradients are
    interpreted in user (un-scaled) space throughout — the
    adjoint solve and scatter chain handles all preconditioner
    bookkeeping internally.

    Raises
    ------
    RuntimeError
        If :meth:`solve` has not been called yet (no cached KKT factor).
    """
    if not self.settings.enable_grad:
        raise RuntimeError("Set enable_grad to True to enable gradient computation.")

    if not getattr(self, "_setup_done", False) or self._result is None:
        raise RuntimeError(
            f"{type(self).__name__}.grad() requires a prior solve(); "
            f"call setup() and solve() before grad()."
        )

    data = self._data
    B = data.batch_size
    wp_stream = wp.Stream(cuda_stream=cp.cuda.get_current_stream().ptr)

    # ---- Step 1
    zeros = self._zero_grad_in
    pack_total = data.n + data.p + 2 * zeros.num_ineq
    if pack_total > 0:
        wp.launch(
            kernel=self._backward_copy_kernel,
            dim=(B, pack_total),
            inputs=[
                grad_x    if grad_x    is not None else zeros.x,
                grad_y    if grad_y    is not None else zeros.y,
                grad_z_u  if grad_z_u  is not None else zeros.z_u,
                grad_z_l  if grad_z_l  is not None else zeros.z_l,
                grad_z_bu if grad_z_bu is not None else zeros.z_bu,
                grad_z_bl if grad_z_bl is not None else zeros.z_bl,
                grad_s_u  if grad_s_u  is not None else zeros.s_u,
                grad_s_l  if grad_s_l  is not None else zeros.s_l,
                grad_s_bu if grad_s_bu is not None else zeros.s_bu,
                grad_s_bl if grad_s_bl is not None else zeros.s_bl,
                self._grad_in.x,    self._grad_in.y,
                self._grad_in.z_u,  self._grad_in.z_l,
                self._grad_in.z_bu, self._grad_in.z_bl,
                self._grad_in.s_u,  self._grad_in.s_l,
                self._grad_in.s_bu, self._grad_in.s_bl,
            ],
            device="cuda",
            stream=wp_stream,
        )

    # ---- Step 2: adjoint KKT solve
    self._compute_adjoint(self._grad_in, self._backward_adjoint_vector)

    # ---- Step 3: write into full-layout buffers
    wp.launch(
        kernel=self._backward_pack_full_layout_kernel,
        dim=(B, 4 * data.m + data.num_xu + data.num_xl),
        inputs=[
            self._backward_adjoint_vector.z_u, self._backward_adjoint_vector.z_l,
            self._backward_adjoint_vector.z_bu, self._backward_adjoint_vector.z_bl,
            self._result.z_u, self._result.z_l,
            self._lam_zu_full, self._lam_zl_full,
            self._lam_zbu_full, self._lam_zbl_full,
            self._zu_full, self._zl_full,
        ],
        device="cuda",
        stream=wp_stream,
    )

    # ---- Step 4: backend-specific matrix and vector gradient
    # assembly + Data subclass construction.
    return self._compute_data_gradients(self._backward_adjoint_vector)

SparseSolver takes a batch of sparse matrices as a single UniformBatchedCsrMatrix — cuPIQP's own batched CSR container (the preferred, fastest batched input; see setup above):

UniformBatchedCsrMatrix 

UniformBatchedCsrMatrix(
    batch_size: int,
    indices: Sequence[int],
    indptr: Sequence[int],
    data: ndarray,
    shape: Optional[Tuple[int, int]] = None,
    dtype=cp.float64,
)

A batch of CSR matrices that share one sparsity pattern.

A batched extension of cupy's cupyx.scipy.sparse.csr_matrix: where a single CSR matrix stores indptr, indices, and a 1-D data array of length nnz, this type stores one shared indptr / indices pair plus a 2-D data buffer of shape (batch_size, nnz) - the values of all matrices stacked along a leading batch axis. Every matrix in the batch therefore has the same nonzero structure and differs only in its values.

This is the storage cuPIQP uses internally for batched sparse problems, and the preferred input for solving a batch with SparseSolver: because the values are contiguous with a uniform per-matrix stride, batched sparse linear-algebra routines can sweep the whole batch with no copy (unlike a Python list of separate csr_matrix objects, which must be copied into this layout first).

Parameters:

Name Type Description Default
batch_size int

Number of matrices in the batch, B (must be positive).

 required
indices sequence of int

The shared CSR column-index and row-pointer arrays - the single sparsity pattern used by every matrix in the batch.

 required
indptr sequence of int

The shared CSR column-index and row-pointer arrays - the single sparsity pattern used by every matrix in the batch.

 required
data ndarray

Values of shape (batch_size, nnz): row i holds the nonzeros of the i-th matrix, laid out against the shared indices / indptr.

 required
shape tuple of (int, int)

Dense (rows, cols) shape of each matrix; inferred from the pattern when omitted.

 None
dtype data - type

Value dtype.

 ``cupy.float64``

Attributes:

Name Type Description
batch_size, nnz, rows, cols, shape

Batch size B, nonzeros per matrix, and the shared dense shape.
indices, indptr ndarray

The shared CSR sparsity pattern.
data ndarray

The (batch_size, nnz) values buffer.
Source code in cupiqp/sparse/batched_csr.py 
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
def __init__(
    self,
    batch_size: int,
    indices: Sequence[int],
    indptr: Sequence[int],
    data: cp.ndarray,
    shape: Optional[Tuple[int, int]] = None,
    dtype=cp.float64,
):
    if not batch_size > 0:
        raise ValueError("batch_size must be a positive integer.")

    self._dtype = dtype
    data_cp = cp.asarray(data, dtype=dtype)
    try:
        if shape is None:
            self._template_matrix = csr_matrix((data_cp[0], indices, indptr))
        else:
            self._template_matrix = csr_matrix(
                (data_cp[0], indices, indptr), shape=shape,
            )
    except Exception as e:
        raise ValueError(
            "Invalid indices, indptr, and data combination."
        ) from e

    self._batch_size = batch_size
    self._indices = self._template_matrix.indices
    self._indptr = self._template_matrix.indptr
    self._nnz = self._template_matrix.nnz

    if data_cp.shape != (batch_size, self._nnz):
        raise ValueError(
            f"data must have shape ({batch_size}, {self._nnz}), got {data_cp.shape}."
        )
    self.data = cp.empty((batch_size, self._nnz), dtype=dtype)
    self.data[:] = data_cp

MultistageSolver

MultistageSolver 

MultistageSolver(
    dtype: Literal["float32", "float64"] = "float64",
)

Bases: SolverBase

GPU solver for multistage (block-structured) convex quadratic programs that solves a QP - or a whole batch of QPs - of the form

\[ \begin{aligned} \min_{x}\quad & \tfrac{1}{2}\, x^\top P x + c^\top x \\ \text{s.t.}\quad & A x = b, \\ & h_l \le G x \le h_u, \\ & x_l \le x \le x_u, \end{aligned} \]

using the proximal interior-point method with a block-Cholesky factorization that exploits block-tridiagonal / block-tridiagonal-arrow KKT structure, running entirely on the GPU. It is built for multistage problems such as optimal control (OCPs / MPC), where that structure arises from the stage-by-stage dynamics and costs.

Inputs - block-structured, end to end. Unlike the dense and sparse solvers, MultistageSolver takes its problem data as pre-built block objects (from cupiqp.multistage.multistage_utils), not arrays or CSR matrices - generic CSR is not auto-promoted to block form, because the structure can only be exploited if you build it explicitly:

  • P : a BlockTridiagMat,
  • A, G : a BlockBidiagMat (or None),
  • c, b, h_l, h_u, x_l, x_u : a BlockVec (or None where the block is absent).

Only P and c are required. As with the other backends, +/-inf entries in the bound blocks mark one-sided or free bounds and are dropped at no numerical cost.

Batching. Build the block objects with batch_size=B and the solver processes all B problems in one GPU call; B = 1 is a single problem, and solver.result.x then has shape (B, n).

Parameters:

Name Type Description Default
dtype (float64, float32)

Floating-point precision used throughout the solve. "float32" is faster and uses less memory but converges to looser tolerances; the default convergence tolerances are chosen to match the dtype.

 "float64"

Examples:

from cupiqp import MultistageSolver
from cupiqp.multistage.multistage_utils import (
    BlockTridiagMat, BlockBidiagMat, BlockVec,
)

P = BlockTridiagMat(num_diag_blocks=N, block_size=d)
A = BlockBidiagMat(rows_of_blocks=d, cols_of_blocks=d, N=N)
c = BlockVec(num_blocks=N, rows=d)
b = BlockVec(num_blocks=N, rows=d)
# ... fill block data ...

s = MultistageSolver()
s.setup(P=P, c=c, A=A, b=b)
s.solve()
See Also

DenseSolver: solver for general dense problems.

SparseSolver: solver for general sparse problems.

Notes

Requires the socu block-solver package (install the multistage extra: pip install ".[cuda13,multistage]"). The problem structure - the block sizes and which blocks are present - is fixed by setup and can only be set up once per instance; to re-solve with new numerical values of the same structure, call update and then solve again, which reuses all GPU allocations. Solver behaviour (tolerances, verbosity, iteration cap, ...) is configured through solver.settings.

Source code in cupiqp/multistage/multistage_solver.py 
142
143
144
def __init__(self, dtype: Literal["float32", "float64"] = "float64"):
    super().__init__(dtype=dtype)
    self._settings.kkt_solver = "multistage_block_cholesky"

setup 

setup(
    P: BlockTridiagMat,
    c: BlockVec,
    A: Optional[BlockBidiagMat] = None,
    b: Optional[BlockVec] = None,
    G: Optional[BlockBidiagMat] = None,
    h_u: Optional[BlockVec] = None,
    h_l: Optional[BlockVec] = None,
    x_u: Optional[BlockVec] = None,
    x_l: Optional[BlockVec] = None,
) -> None

Bind the problem data and prepare the solver for solve().

Fixes the problem structure - array shapes, which constraint blocks are present, the sparsity pattern (sparse backend), and the finite/infinite pattern of the bounds - and allocates all GPU buffers, the KKT system, and the preconditioner. Call this once per solver instance, then call solve().

Pass a single problem (2D P) or a batch (3D P with a leading batch axis, or a list of matrices); the batch size is inferred here and sets the shape of the result.

Parameters:

Name Type Description Default
P GPU array

Quadratic cost, shape (n, n) or batched (B, n, n). Must be symmetric positive semidefinite. Required.

 required
c GPU array

Linear cost, shape (n,) or (B, n). Required.

 required
A GPU array

Equality constraints A x = b; shapes (p, n) and (p,) (or batched). Omit for no equality constraints.

 None
b GPU array

Equality constraints A x = b; shapes (p, n) and (p,) (or batched). Omit for no equality constraints.

 None
G GPU array

Two-sided inequalities h_l <= G x <= h_u; G is (m, n) (or batched) and the bounds are (m,). Use -inf / +inf entries for one-sided rows.

 None
h_l GPU array

Two-sided inequalities h_l <= G x <= h_u; G is (m, n) (or batched) and the bounds are (m,). Use -inf / +inf entries for one-sided rows.

 None
h_u GPU array

Two-sided inequalities h_l <= G x <= h_u; G is (m, n) (or batched) and the bounds are (m,). Use -inf / +inf entries for one-sided rows.

 None
x_l GPU array

Element-wise box bounds x_l <= x <= x_u, shape (n,) (or batched). Use +/-inf for unbounded entries.

 None
x_u GPU array

Element-wise box bounds x_l <= x <= x_u, shape (n,) (or batched). Use +/-inf for unbounded entries.

 None

Raises:

Type Description
RuntimeError

If setup() has already been called on this instance. The structure is fixed after setup - create a new solver for a different structure, or use update() to change only the numerical values.
TypeError

If an input is not a GPU array of the kind this backend expects (e.g. a CPU numpy array, or a dense matrix passed to the sparse backend). See the backend's class docstring for the exact accepted types.
See Also

solve : run the solver after setup. update : change numerical data without a full re-setup.

Source code in cupiqp/multistage/multistage_solver.py 
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
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
def setup(
    self,
    P: BlockTridiagMat,
    c: BlockVec,
    A: Optional[BlockBidiagMat] = None,
    b: Optional[BlockVec] = None,
    G: Optional[BlockBidiagMat] = None,
    h_u: Optional[BlockVec] = None,
    h_l: Optional[BlockVec] = None,
    x_u: Optional[BlockVec] = None,
    x_l: Optional[BlockVec] = None,
) -> None:
    super().setup(P, c, A, b, G, h_u, h_l, x_u, x_l)
    if self.settings.enable_grad:
        d = self._data
        B = d.batch_size
        N = d.num_blocks
        d_sz = d.block_size
        dtype = d.dtype
        wp_dtype = to_warp_dtype(dtype)

        r_a = d._A.rows_of_blocks if d.p > 0 else 0
        N_a = d._A.N              if d.p > 0 else 0
        r_g = d._G.rows_of_blocks if d.m > 0 else 0
        N_g = d._G.N              if d.m > 0 else 0

        placeholder_P = BlockTridiagMat(
            num_diag_blocks=N, block_size=d_sz, batch_size=B, dtype=wp_dtype,
        )
        placeholder_A = BlockBidiagMat(
            rows_of_blocks=r_a, cols_of_blocks=d_sz, N=N_a, batch_size=B, dtype=wp_dtype,
        ) if d.p > 0 else None
        placeholder_G = BlockBidiagMat(
            rows_of_blocks=r_g, cols_of_blocks=d_sz, N=N_g, batch_size=B, dtype=wp_dtype,
        ) if d.m > 0 else None
        # c: always (B, n). Box blocks x_u / x_l are optional -- their
        # gradient block exists only when that side was provided at
        # setup(), matching the full-length contract (grad_data.x_u is
        # (B, 0) / absent when the forward problem omitted it).
        placeholder_c = BlockVec(num_blocks=N, rows=d_sz, batch_size=B, dtype=wp_dtype)
        placeholder_xu = BlockVec(num_blocks=N, rows=d_sz, batch_size=B, dtype=wp_dtype) if d.has_x_u else None
        placeholder_xl = BlockVec(num_blocks=N, rows=d_sz, batch_size=B, dtype=wp_dtype) if d.has_x_l else None
        # b: presence tied to A existing. h_u / h_l are optional inequality
        # sides -- their gradient block exists only when that side was
        # provided at setup() (grad_data.h_u is (B, 0) / absent when the
        # forward problem omitted it), matching the box-block contract.
        placeholder_b  = BlockVec(num_blocks=N_a + 1, rows=r_a, batch_size=B, dtype=wp_dtype) if d.p > 0 else None
        placeholder_hu = BlockVec(num_blocks=N_g + 1, rows=r_g, batch_size=B, dtype=wp_dtype) if d.has_h_u else None
        placeholder_hl = BlockVec(num_blocks=N_g + 1, rows=r_g, batch_size=B, dtype=wp_dtype) if d.has_h_l else None

        self._grad_data = MultistageData(dtype=dtype, device=self.settings.device)
        self._grad_data.init(
            P=placeholder_P, c=placeholder_c,
            A=placeholder_A, b=placeholder_b,
            G=placeholder_G, h_u=placeholder_hu, h_l=placeholder_hl,
            x_u=placeholder_xu, x_l=placeholder_xl,
        )

        # Empty placeholder warp buffers for when A or G are absent —
        # the kernel still needs valid array arguments even though its
        # corresponding dispatch sub-range collapses to size 0.
        empty_blocks = wp.zeros((B, 0, 0, 0), dtype=wp_dtype, device="cuda")
        g = self._grad_data
        self._grad_dA_D = g._A.D if g._A is not None else empty_blocks
        self._grad_dA_E = g._A.E if g._A is not None else empty_blocks
        self._grad_dG_D = g._G.D if g._G is not None else empty_blocks
        self._grad_dG_E = g._G.E if g._G is not None else empty_blocks

        # Eager-compile the fused multistage data-gradients kernel.
        self._multistage_data_gradients_kernel = create_multistage_data_gradients_kernel(
            N, d_sz, N_a, r_a, N_g, r_g, d.p, d.m, d.n, d.num_hu, d.num_hl, d.num_xu, d.num_xl, dtype=dtype)

solve 

solve() -> List[Status]

Solve the QP set up by setup() and return the solve status.

Runs the proximal interior-point iterations on the GPU. The full solution (primal x, dual, and slack variables) and per-problem diagnostics are written to solver.result; this method returns the status for convenience.

Returns:

Type Description
Status or list of Status

For a single problem, the Status. For a batched setup(), a list of B of them (one per problem). CUPIQP_SOLVED means the problem converged to tolerance. The list form is always available as solver.result.info.status.
Notes

Read the solution from solver.result after solving - e.g. solver.result.x (shape (B, n)) and solver.result.info.status. Set solver.settings.verbose = True to print a per-iteration log. After setup() you may solve() repeatedly, optionally calling update() in between to change the numerical data.

Source code in cupiqp/solver.py 
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
def solve(self) -> List[Status]:
    """Solve the QP set up by ``setup()`` and return the solve status.

    Runs the proximal interior-point iterations on the GPU. The full
    solution (primal ``x``, dual, and slack variables) and per-problem
    diagnostics are written to ``solver.result``; this method returns the
    status for convenience.

    Returns
    -------
    Status or list of Status
        For a single problem, the ``Status``. For a batched ``setup()``,
        a list of ``B`` of them (one per problem). ``CUPIQP_SOLVED``
        means the problem converged to tolerance. The list form is always
        available as ``solver.result.info.status``.

    Notes
    -----
    Read the solution from ``solver.result`` after solving - e.g.
    ``solver.result.x`` (shape ``(B, n)``) and
    ``solver.result.info.status``. Set ``solver.settings.verbose = True``
    to print a per-iteration log. After ``setup()`` you may ``solve()``
    repeatedly, optionally calling ``update()`` in between to change the
    numerical data.
    """
    if self.settings.verbose:
        try:
            from importlib.metadata import version
            _ver = version("cupiqp")
        except Exception:
            _ver = ""
        _w = 58
        print("-" * _w)
        print(f"cuPIQP v{_ver} - GPU-accelerated PIQP solver".strip().center(_w))
        print("(c) Fenglong Song".center(_w))
        print("Ecole Polytechnique Federale de Lausanne (EPFL) 2026".center(_w))
        print("-" * _w)
        if self.settings.kkt_solver == "dense_cholesky":
            print("dense backend:")
            print(f"batch size B = {self._data.batch_size}")
            print(f"variables n = {self._data.n}")
            print(f"equality constraints p = {self._data.p}")
            print(f"inequality constraints m = {self._data.m}")
        elif self.settings.kkt_solver == "sparse_ldlt":
            print("sparse backend:")
            print(f"batch size B = {self._data.batch_size}")
            print(f"variables n = {self._data.n}, nnz(P) = {self._data.P.nnz}")
            print(f"equality constraints p = {self._data.p}, nnz(A) = {self._data.A.nnz}")
            print(f"inequality constraints m = {self._data.m}, nnz(G) = {self._data.G.nnz}")
        elif self.settings.kkt_solver == "multistage_block_cholesky":
            print("multistage backend:")
            print(f"batch size B = {self._data.batch_size}")
            print(f"variables n = {self._data.n}, num_diag_blocks(P) = {self._data.P.num_diag_blocks}, block_size(P) = ({self._data.P.block_size}, {self._data.P.block_size})")
            print(f"equality constraints p = {self._data.p}, num_diag_blocks(A) = {self._data.A.N}, block_size(A) = ({self._data.A.rows_of_blocks}, {self._data.A.cols_of_blocks})")
            print(f"inequality constraints m = {self._data.m}, num_diag_blocks(G) = {self._data.G.N}, block_size(G) = ({self._data.G.rows_of_blocks}, {self._data.G.cols_of_blocks})")
        else:
            raise ValueError(f"Unsupported kkt_solver type: {self.settings.kkt_solver}")

        print(f"inequality lower bounds n_h_l = {self._data.num_hl}")
        print(f"inequality upper bounds n_h_u = {self._data.num_hu}")
        print(f"variable lower bounds n_x_l = {self._data.num_xl}")
        print(f"variable upper bounds n_x_u = {self._data.num_xu}")
        print("")
    return self._solve_impl()

update 

update(
    P: Optional[Any] = None,
    c: Optional[Any] = None,
    A: Optional[Any] = None,
    b: Optional[Any] = None,
    G: Optional[Any] = None,
    h_u: Optional[Any] = None,
    h_l: Optional[Any] = None,
    x_u: Optional[Any] = None,
    x_l: Optional[Any] = None,
    check_validity: bool = False,
)

Change the numerical problem data, then solve() again.

The fast path for re-solving a problem of the same structure - for example a moving target b or a re-linearized P in receding-horizon control. It reuses every GPU allocation from setup(), so only the values change; shapes, sparsity patterns, and which blocks are present must stay the same (create a new solver for a structural change). Bound values may change freely, including which entries are +/-inf - a bound can flip between finite and infinite without re-setup().

Any argument left as None keeps its current value. After update(), call solve() to get the new solution.

Parameters:

Name Type Description Default
P GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
c GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
A GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
b GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
G GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
h_u GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
h_l GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
x_u GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
x_l GPU array

New values for the corresponding problem block. None (the default) leaves that block unchanged. Must match the original shapes / sparsity pattern set at setup().

 None
check_validity bool

If True, validate the dimensions and sparsity of the new data. Defaults to False for speed (validation forces device-to-host syncs in the sparse backend). When False, you must still keep the shapes and sparsity patterns of P/A/ G unchanged; bound values (including which entries are +/-inf) may change.

 False
Source code in cupiqp/solver.py 
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
def update(self,
           P: Optional[Any] = None,
           c: Optional[Any] = None,
           A: Optional[Any] = None,
           b: Optional[Any] = None,
           G: Optional[Any] = None,
           h_u: Optional[Any] = None,
           h_l: Optional[Any] = None,
           x_u: Optional[Any] = None,
           x_l: Optional[Any] = None,
           check_validity: bool = False,
           ):
    """Change the numerical problem data, then ``solve()`` again.

    The fast path for re-solving a problem of the **same structure** -
    for example a moving target ``b`` or a re-linearized ``P`` in
    receding-horizon control. It reuses every GPU allocation from
    ``setup()``, so only the values change; shapes, sparsity patterns,
    and which blocks are present must stay the same (create a new solver
    for a structural change). Bound *values* may change freely, including
    which entries are ``+/-inf`` - a bound can flip between finite and
    infinite without re-``setup()``.

    Any argument left as ``None`` keeps its current value. After
    ``update()``, call ``solve()`` to get the new solution.

    Parameters
    ----------
    P, c, A, b, G, h_u, h_l, x_u, x_l : GPU array, optional
        New values for the corresponding problem block. ``None`` (the
        default) leaves that block unchanged. Must match the original
        shapes / sparsity pattern set at ``setup()``.
    check_validity : bool, default: False
        If ``True``, validate the dimensions and sparsity of the new
        data. Defaults to ``False`` for speed (validation forces
        device-to-host syncs in the sparse backend). When ``False``, you
        must still keep the shapes and sparsity patterns of ``P``/``A``/
        ``G`` unchanged; bound values (including which entries are
        ``+/-inf``) may change.
    """
    if not self._setup_done:
        raise RuntimeError("Solver not setup yet. Call setup() first.")

    if self.settings.preconditioner_iter > 0:
        self._preconditioner.unscale_data(self._data)

    if P is not None:
        self._data.set_P(P, check=check_validity)
    if c is not None:
        self._data.set_c(c, check=check_validity)
    if A is not None:
        self._data.set_A(A, check=check_validity)
    if b is not None:
        self._data.set_b(b, check=check_validity)
    if G is not None:
        self._data.set_G(G, check=check_validity)
    if h_u is not None:
        self._data.set_h_u(h_u, check=check_validity)
    if h_l is not None:
        self._data.set_h_l(h_l, check=check_validity)
    if x_u is not None:
        self._data.set_x_u(x_u, check=check_validity)
    if x_l is not None:
        self._data.set_x_l(x_l, check=check_validity)

    matrix_changed = P is not None or A is not None or G is not None

    # NOTE: Since we allow changing h_l/h_u containing arbitrary +inf/-inf, 
    # an inequality row G[i] can switch between active (a finite
    # bound) and inactive (both bounds infinite) between updates.
    # If either of h_l or h_u are updated, we need to update G 
    # because for sparse kkt solver we need to set the inactive rows to 0
    ineq_bound_pattern_may_change = h_l is not None or h_u is not None

    # Apply preconditioner scaling to updated data.
    preconditioner_did_fresh_ruiz = False
    if self.settings.preconditioner_iter > 0:
        reuse = self.settings.preconditioner_reuse_on_update or not matrix_changed
        if reuse:
            self._preconditioner.reuse_scaling(self._data)
        else:
            self._preconditioner.reset()
            self._preconditioner.scale_data(
                self._data,
                self.settings.preconditioner_scale_cost,
                self.settings.preconditioner_iter,
            )
            preconditioner_did_fresh_ruiz = True

    self._preconditioner.compute_constraints_rhs_inf_norm_unscaled(
        self._data, self._constraints_rhs_inf_norm_unscaled,
    )
    # Fresh Ruiz produces new factors that re-scale ALL of P/A/G in place,
    # even matrices the user didn't pass. The KKT solver caches things
    # like A^T A keyed off those scaled values, so flag everything as
    # changed in that case.
    self._kkt_system.update_data(
        self._data,
        (P is not None) or preconditioner_did_fresh_ruiz,
        (A is not None) or preconditioner_did_fresh_ruiz,
        (G is not None) or preconditioner_did_fresh_ruiz or ineq_bound_pattern_may_change,
    )

backward 

backward(
    grad_x=None,
    grad_y=None,
    grad_z_u=None,
    grad_z_l=None,
    grad_z_bu=None,
    grad_z_bl=None,
    grad_s_u=None,
    grad_s_l=None,
    grad_s_bu=None,
    grad_s_bl=None,
)

Compute gradients of an outer scalar :math:L w.r.t. problem data, given upstream cotangents on the solution variables.

Orchestration (backend-agnostic):

  1. Pack the per-field cotangent kwargs into self._grad_in (a pre-allocated Variables); missing kwargs are treated as zeros.
  2. Solve the adjoint KKT system via :meth:_compute_adjoint, producing user-space adjoint vectors.
  3. Scatter the four active-size lambda groups (z_u, z_l, z_bu, z_bl) and the two active-size ineq result groups into full-m / full-n buffers (self._lam_z*_full, self._z*_full). Both the dG outer product and the dh_* / dx_* vector gradients consume these full-layout buffers.
  4. Delegate to :meth:_compute_data_gradients for backend- specific matrix-gradient assembly + Data subclass construction.

Returns the backend's Data subclass populated with the gradients in user space. Cotangents and returned gradients are interpreted in user (un-scaled) space throughout — the adjoint solve and scatter chain handles all preconditioner bookkeeping internally.

Raises:

Type Description
RuntimeError

If :meth:solve has not been called yet (no cached KKT factor).
Source code in cupiqp/solver.py 
1828
1829
1830
1831
1832
1833
1834
1835
1836
1837
1838
1839
1840
1841
1842
1843
1844
1845
1846
1847
1848
1849
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880
1881
1882
1883
1884
1885
1886
1887
1888
1889
1890
1891
1892
1893
1894
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
1907
1908
1909
1910
1911
1912
1913
1914
1915
1916
1917
1918
1919
1920
1921
1922
1923
1924
1925
1926
@nvtx.annotate("Solver::grad")
def backward(self,
         grad_x=None, grad_y=None,
         grad_z_u=None, grad_z_l=None, grad_z_bu=None, grad_z_bl=None,
         grad_s_u=None, grad_s_l=None, grad_s_bu=None, grad_s_bl=None):
    r"""Compute gradients of an outer scalar :math:`L` w.r.t. problem
    data, given upstream cotangents on the solution variables.

    Orchestration (backend-agnostic):

    1. Pack the per-field cotangent kwargs into ``self._grad_in``
       (a pre-allocated ``Variables``); missing kwargs are treated
       as zeros.
    2. Solve the adjoint KKT system via :meth:`_compute_adjoint`,
       producing user-space adjoint vectors.
    3. Scatter the four active-size lambda groups (``z_u, z_l,
       z_bu, z_bl``) and the two active-size ineq result groups
       into full-``m`` / full-``n`` buffers
       (``self._lam_z*_full``, ``self._z*_full``). Both the dG
       outer product and the ``dh_*`` / ``dx_*`` vector gradients
       consume these full-layout buffers.
    4. Delegate to :meth:`_compute_data_gradients` for backend-
       specific matrix-gradient assembly + ``Data`` subclass
       construction.

    Returns the backend's ``Data`` subclass populated with the
    gradients in user space. Cotangents and returned gradients are
    interpreted in user (un-scaled) space throughout — the
    adjoint solve and scatter chain handles all preconditioner
    bookkeeping internally.

    Raises
    ------
    RuntimeError
        If :meth:`solve` has not been called yet (no cached KKT factor).
    """
    if not self.settings.enable_grad:
        raise RuntimeError("Set enable_grad to True to enable gradient computation.")

    if not getattr(self, "_setup_done", False) or self._result is None:
        raise RuntimeError(
            f"{type(self).__name__}.grad() requires a prior solve(); "
            f"call setup() and solve() before grad()."
        )

    data = self._data
    B = data.batch_size
    wp_stream = wp.Stream(cuda_stream=cp.cuda.get_current_stream().ptr)

    # ---- Step 1
    zeros = self._zero_grad_in
    pack_total = data.n + data.p + 2 * zeros.num_ineq
    if pack_total > 0:
        wp.launch(
            kernel=self._backward_copy_kernel,
            dim=(B, pack_total),
            inputs=[
                grad_x    if grad_x    is not None else zeros.x,
                grad_y    if grad_y    is not None else zeros.y,
                grad_z_u  if grad_z_u  is not None else zeros.z_u,
                grad_z_l  if grad_z_l  is not None else zeros.z_l,
                grad_z_bu if grad_z_bu is not None else zeros.z_bu,
                grad_z_bl if grad_z_bl is not None else zeros.z_bl,
                grad_s_u  if grad_s_u  is not None else zeros.s_u,
                grad_s_l  if grad_s_l  is not None else zeros.s_l,
                grad_s_bu if grad_s_bu is not None else zeros.s_bu,
                grad_s_bl if grad_s_bl is not None else zeros.s_bl,
                self._grad_in.x,    self._grad_in.y,
                self._grad_in.z_u,  self._grad_in.z_l,
                self._grad_in.z_bu, self._grad_in.z_bl,
                self._grad_in.s_u,  self._grad_in.s_l,
                self._grad_in.s_bu, self._grad_in.s_bl,
            ],
            device="cuda",
            stream=wp_stream,
        )

    # ---- Step 2: adjoint KKT solve
    self._compute_adjoint(self._grad_in, self._backward_adjoint_vector)

    # ---- Step 3: write into full-layout buffers
    wp.launch(
        kernel=self._backward_pack_full_layout_kernel,
        dim=(B, 4 * data.m + data.num_xu + data.num_xl),
        inputs=[
            self._backward_adjoint_vector.z_u, self._backward_adjoint_vector.z_l,
            self._backward_adjoint_vector.z_bu, self._backward_adjoint_vector.z_bl,
            self._result.z_u, self._result.z_l,
            self._lam_zu_full, self._lam_zl_full,
            self._lam_zbu_full, self._lam_zbl_full,
            self._zu_full, self._zl_full,
        ],
        device="cuda",
        stream=wp_stream,
    )

    # ---- Step 4: backend-specific matrix and vector gradient
    # assembly + Data subclass construction.
    return self._compute_data_gradients(self._backward_adjoint_vector)

MultistageSolver takes its problem data as the block-structured objects below — build them, fill in their data, and pass them to setup (see above for which argument expects which type):

BlockTridiagMat 

BlockTridiagMat(
    num_diag_blocks: int,
    block_size: int,
    dtype=wp.float64,
    device="cuda",
    batch_size: int = 1,
)

Symmetric block-tridiagonal matrix, batched.

Layout::

D   shape (B, N, d, d)     - diagonal blocks
E   shape (B, N-1, d, d)   - lower-diagonal blocks

The matrix is symmetric; the upper off-diagonal is implicitly the transpose of the lower one. batch_size defaults to 1.

diag_blocks and off_diag_blocks_lower may be assigned a Warp array or any CUDA array implementing cuda_array_interface (CuPy, dense CUDA PyTorch, JAX, Numba, ...). Non-Warp arrays are wrapped zero-copy; shape and dtype must match the construction arguments.

Source code in cupiqp/multistage/multistage_utils.py 
24
25
26
27
28
29
30
31
32
def __init__(self, num_diag_blocks: int, block_size: int,
             dtype=wp.float64, device="cuda", batch_size: int = 1):
    self.block_size = block_size
    self._device = device
    self._dtype = to_warp_dtype(dtype)
    self._D_shape = (batch_size, num_diag_blocks, block_size, block_size)
    self._E_shape = (batch_size, num_diag_blocks - 1, block_size, block_size)
    self.D = wp.zeros(self._D_shape, dtype=self._dtype, device=device)
    self.E = wp.zeros(self._E_shape, dtype=self._dtype, device=device)

BlockBidiagMat 

BlockBidiagMat(
    rows_of_blocks: int,
    cols_of_blocks: int,
    N: int,
    dtype=wp.float64,
    device="cuda",
    batch_size: int = 1,
)

Block lower-bidiagonal matrix, batched.

Stores A or G of the multistage problem with shape::

D shape (B, N, rows_of_blocks, cols_of_blocks)
E shape (B, N, rows_of_blocks, cols_of_blocks)

Logical structure (per batch)::

A =
[ D0
  E0  D1
      E1  D2
          ...
          E_{N-2} D_{N-1}
                  E_{N-1} ]

D and E may be assigned a warp array or any CUDA array implementing __cuda_array_interface__ (cupy, dense CUDA torch, jax, numba, ...). Non-warp arrays are wrapped zero-copy, so the container aliases the assigned buffer; shape and dtype must match the construction arguments.

Source code in cupiqp/multistage/multistage_utils.py 
117
118
119
120
121
122
123
124
125
126
def __init__(self, rows_of_blocks: int, cols_of_blocks: int, N: int,
             dtype=wp.float64, device="cuda", batch_size: int = 1):
    self.N = N
    self.cols_of_blocks = cols_of_blocks
    self.rows_of_blocks = rows_of_blocks
    self._device = device
    self._dtype = to_warp_dtype(dtype)
    self._shape = (batch_size, N, rows_of_blocks, cols_of_blocks)
    self._D = wp.zeros(self._shape, dtype=self._dtype, device=device)
    self._E = wp.zeros(self._shape, dtype=self._dtype, device=device)

BlockVec 

BlockVec(
    num_blocks: int,
    rows: int,
    dtype=wp.float64,
    device="cuda",
    batch_size: int = 1,
)

Block vector, batched. data shape (B, num_blocks, rows).

data may be assigned a warp array or any CUDA array implementing __cuda_array_interface__ (cupy, dense CUDA torch, jax, numba, ...). Non-warp arrays are wrapped zero-copy, so the container aliases the assigned buffer; shape and dtype must match the construction arguments.

Source code in cupiqp/multistage/multistage_utils.py 
178
179
180
181
182
183
184
185
def __init__(self, num_blocks: int, rows: int,
             dtype=wp.float64, device="cuda", batch_size: int = 1):
    self.num_blocks = num_blocks
    self.rows = rows
    self._device = device
    self._dtype = to_warp_dtype(dtype)
    self._shape = (batch_size, num_blocks, rows)
    self._data = wp.zeros(self._shape, dtype=self._dtype, device=device)