Accept any shape in geo_mean, and fix row vectors along the way - #3493

@PTNobel on cvxpy#3493: one PR, one implementation. Three things came out of
the review and they all point at AxisAtom, so geo_mean now inherits from
it rather than stacking scalars.

axis and keepdims follow cp.mean's spelling, and tuples of axes come
with them: normalize_axis handles tuples, negative indices, and reducing
over every axis. Shapes and values were checked against
exp(mean(log(x), axis)) for 2-D and 3-D inputs, every single axis and
every pair, with and without keepdims.

The stacking is gone. gm was already vectorised — it takes t.size and
emits one SOC over all entries — and PowConeND takes one cone per column,
so both canonicalizers now cover every output with a single set of
constraints. Measured through geo_mean_approx_canon: 3 SOC constraints
for 2 outputs, and 3 for 200. A test pins that down.

Zero-weight entries used to be indexed out of the argument itself, which
put a stray [True, True, True] into every geo_mean's printed form and
would not have survived an axis. They are now dropped in the aligned view
the canonicalizers work on, so the argument keeps the shape the axis
bookkeeping is based on, and the two label expectations lose the stray
index.

This also fixes a row vector reaching the solver: geo_mean(Variable((1, 4)))
built and evaluated fine on master and then raised IndexError during
canonicalization, because the canonicalizer indexed along axis 0.

Closes the geo_mean half of cvxpy#3489

@PTNobel on #3493: one PR, one implementation. Three things came out of
the review and they all point at AxisAtom, so geo_mean now inherits from
it rather than stacking scalars.

axis and keepdims follow cp.mean's spelling, and tuples of axes come
with them: normalize_axis handles tuples, negative indices, and reducing
over every axis. Shapes and values were checked against
exp(mean(log(x), axis)) for 2-D and 3-D inputs, every single axis and
every pair, with and without keepdims.

The stacking is gone. gm was already vectorised — it takes t.size and
emits one SOC over all entries — and PowConeND takes one cone per column,
so both canonicalizers now cover every output with a single set of
constraints. Measured through geo_mean_approx_canon: 3 SOC constraints
for 2 outputs, and 3 for 200. A test pins that down.

Zero-weight entries used to be indexed out of the argument itself, which
put a stray [True, True, True] into every geo_mean's printed form and
would not have survived an axis. They are now dropped in the aligned view
the canonicalizers work on, so the argument keeps the shape the axis
bookkeeping is based on, and the two label expectations lose the stray
index.

This also fixes a row vector reaching the solver: geo_mean(Variable((1, 4)))
built and evaluated fine on master and then raised IndexError during
canonicalization, because the canonicalizer indexed along axis 0.

Closes the geo_mean half of #3489

(cherry picked from commit 1f673ee)

The pyright workflow ran `uvx pyright` with no environment configured, so
numpy/scipy imports resolved to Unknown. With #3493's geo_mean rewrite this
tripped reportOptionalIterable: `p = p[self._keep]` assigns an Unknown,
which resets p's narrowed type to its declared Optional parameter type.

Port the workflow half of #3339 from master: add the dev dependency group
(pyright, scipy-stubs) and run pyright via `uv run --no-sync` so it
analyzes real NumPy/SciPy types. The full run is at 0 errors on this branch.

Co-Authored-By: Claude Fable 5 <[email protected]>

* Add dense and parametric quad_form bindings to the NLP diffengine (#3365)

* Add dense and parametric quad_form bindings to the NLP diff engine

convert_quad_form now uses the format-string make_quad_form binding with
three paths: sparse constant P (CSR), dense constant P (permuted_dense),
and parametric P (fed to the dense path as a matrix-valued child that is
re-evaluated each solve). Previously only constant sparse P was supported
and a parametric P raised NotImplementedError.

Adds TestQuadFormDiffEngine (DerivativeChecker-based, no NLP solver
required) covering dense/sparse/parametric P over a leaf variable and
dense/sparse P over a composition (x[:k]); removes the obsolete
test_parameterized_quad_form_errors_clearly.

Bumps the sparsediffpy pin to >=0.4.0,<0.5.0, which provides the
format-string make_quad_form binding the new paths require.

Co-Authored-By: Claude Opus 4.8 (1M context) <[email protected]>

* Add quad_form diff-engine tests from review; drop permuted_dense comments

Address dance858's review on the dense/parametric quad_form binding:
- nonlinear composition: quad_form(sin(x .* x), P) with dense P
- diag(p) parametric matrix: quad_form(x, cp.diag(p)) with a vector parameter
- parametric P with an updated value (solve, update, re-solve, re-check) in the
  NLP parameter tests file

Also remove the diff-engine-internal "permuted_dense" name from convert_quad_form's
docstring/comment and the quad_form test comment.

Co-Authored-By: Claude Opus 4.8 (1M context) <[email protected]>

* Route a dense-but-sparse quad_form P to the sparse binding

convert_quad_form density-checks a dense constant P (nnz/size) and routes a mostly-zero
one (e.g. a diagonal written as a dense ndarray, H*np.eye) to the sparse make_quad_form
binding instead of building a dense Hessian block. Shares SPARSE_DENSITY_THRESHOLD with
the matmul sparse-dispatch. Adds TestQuadFormDenseSparseDispatch.

Co-Authored-By: Claude Opus 4.8 (1M context) <[email protected]>

* bump sparsediffpy version

* Clarify convert_quad_form docstring; drop misleading local-Hessian note

Address review (@PTNobel): the parenthetical claiming a parametric P keeps
"the Hessian ... 2P" referred to the local node Hessian of x'Px, not the
problem/parametric handling, and was misleading for block / known-sparsity P.
Rewrite the docstring to plainly describe the constant (sparse/dense, with
dense-mostly-zero routed to sparse) and parametric (matrix-valued child
re-evaluated each solve) routing. Comment-only; no behavior change.

Co-Authored-By: Claude Opus 4.8 (1M context) <[email protected]>

---------

Co-authored-by: Claude Opus 4.8 (1M context) <[email protected]>

* Route super-sparse constant left-matmul to the sparse diffengine binding (#3366)

* Route super-sparse constant left-matmul to the sparse diffengine binding

A constant dense matmul operand stored as a numpy array but mostly zeros
(e.g. the murray benchmark, ~0.2% nonzero) was always sent down the dense
permuted_dense path, which materializes a fully dense Jacobian/Hessian and
is orders of magnitude slower and heavier than the sparse CSR binding.

convert_matmul now checks the nonzero density of a constant dense left
operand and, when it falls below SPARSE_MATMUL_DENSITY_THRESHOLD (5%),
converts it to CSR and dispatches to the existing sparse left-matmul path.
The check is gated on pure constants (param_node is None): sparsifying a
parametric matrix would freeze its sparsity pattern to the current value.

The density count runs once per matmul node at C-graph construction time,
not per solver iteration.

Co-Authored-By: Claude Opus 4.8 (1M context) <[email protected]>

* Address review: move sparse-matmul threshold to settings, couple tests

Move SPARSE_MATMUL_DENSITY_THRESHOLD out of converters.py into cvxpy.settings
(alongside DPP_PARAM_THRESHOLD) so it is not a magic constant floating in the
module, and have convert_matmul reference s.SPARSE_MATMUL_DENSITY_THRESHOLD.

Tests now derive their densities from the constant instead of hardcoding 0.01 /
0.10, so they stay coupled to it: below-threshold cases use THRESHOLD / 5 and the
block-diagonal case asserts blk / n > THRESHOLD.

Co-Authored-By: Claude Opus 4.8 (1M context) <[email protected]>

* Clean up comments: drop diff-engine-internal permuted_dense references

permuted_dense is an internal name of the diff engine; refer to it neutrally as
"the dense path" in cvxpy-side comments and test assertions. Also trim the
settings.py threshold comment to a concise one-liner matching its neighbors.
Comment-only changes; no behavior change.

Co-Authored-By: Claude Opus 4.8 (1M context) <[email protected]>

* Tighten TestSparseMatmulDispatch comments and test names

Trim the verbose class/per-test comments to the non-obvious essentials and rename
the tests to short, symmetric names (below/above threshold, parametric). No
behavior change.

Co-Authored-By: Claude Opus 4.8 (1M context) <[email protected]>

* Address review: sparray + shared density constant; skip eager constant matmul operand

- Rename SPARSE_MATMUL_DENSITY_THRESHOLD -> SPARSE_DENSITY_THRESHOLD and reuse it in
  Leaf._sparse_high_fill_in (leaf.py) instead of the hardcoded 0.25 (Parth review).
- Migrate the sparse-dispatch path to scipy sparray (csr_array) over spmatrix
  (csr_matrix) in convert_matmul and make_sparse_{left,right}_matmul (Parth review).
- convert_expr no longer eagerly converts a pure-constant matmul operand: convert_matmul
  consumes it as raw arg.value, so building (and densifying a sparse constant into) a node
  it discards is wasted. Avoids a dense Jacobian/parameter blow-up (e.g. sqrt(D) in a
  factor-covariance model). Adds test_constant_matmul.py.

Co-Authored-By: Claude Opus 4.8 (1M context) <[email protected]>

* edits to the PR

---------

Co-authored-by: Claude Opus 4.8 (1M context) <[email protected]>

* Fix scalar boolean indexing to use NumPy mask semantics (#3393)

* Fix scalar boolean indexing to use NumPy mask semantics

x[True] previously meant three different things: is_special_slice
classified Python/NumPy bool scalars as ordinary numbers (bool passes
both numbers.Number and np.isscalar), so validate_key/format_slice
converted True to integer index 1 for shape inference and the
solve/linop path, while index.numeric applied the original key with
NumPy bool-mask semantics (new leading axis). Shape, .value, and solve
results all disagreed.

Fix: route keys containing scalar booleans to special_index, which
already implements full NumPy semantics self-consistently (it builds
its selection matrix by applying the key with NumPy itself). Now
x[True] == x[None] with shape (1, *shape), x[False] yields a 0-size
leading axis, and bools inside tuple keys work, with shape, value, and
solve all agreeing. Also make special_key_to_str tuple-wrap the key so
name() works for non-iterable keys such as a bare bool.

Adds test_scalar_bool_indices covering True/False/np.True_, tuple keys,
and shape/value/solve agreement via a CLARABEL solve; existing boolean
array-mask and integer/slice indexing tests are unchanged and pass.

Review follow-up: assert the solve status is OPTIMAL before checking the
solution value in the new test.

Co-Authored-By: Claude Fable 5 <[email protected]>

* Preserve elementwise label formatting for array index keys

special_key_to_str only needs to wrap scalar keys (the new scalar
boolean case); wrapping every non-tuple key changed the printed label
of boolean-array indices from x[True, True, True] to x[[ True  True
True]], breaking test_labels.

Co-Authored-By: Claude Fable 5 <[email protected]>

---------

Co-authored-by: Claude Fable 5 <[email protected]>

* Fix CUOPT crash on time limit with no incumbent solution (#3443)

When a MILP or LP solve hits a time/iteration limit (TimeLimit ->
USER_LIMIT) before cuOpt finds any feasible incumbent, cuOpt returns an
empty primal solution array. The CUOPT interface still reported the
status as USER_LIMIT (a SOLUTION_PRESENT status), so invert() tried to
reshape the 0-length array into the problem's variable shapes and crashed
with "ValueError: cannot reshape array of size 0 into shape (...)".

Downgrade the status to INFEASIBLE_INACCURATE when the primal solution is
empty, matching how gurobi_conif.py and copt_conif.py already handle the
same "USER_LIMIT with no incumbent" case. solve() now returns cleanly
with prob.value == None and a "solution may be inaccurate" warning instead
of raising an opaque reshape error.

Reported at NVIDIA-AI-Blueprints/portfolio-optimization#53.

Signed-off-by: Ramakrishna Prabhu <[email protected]>
Co-authored-by: Claude Sonnet 5 <[email protected]>

* Add regression tests for duplicate-index gather multiply (gh-3442) (#3446)

The gh-3442 shape — cp.multiply of gathers whose index arrays contain
duplicates, so each gather Jacobian has more nonzeros than its source
variable — corrupted the heap during derivative-structure
initialization. The fix (SparseDiffEngine#105) ships in sparsediffpy
0.6.0: bump the requirement to >= 0.6.0, < 0.7.0 and exercise the fixed
path directly with a solver-free DerivativeChecker test plus an IPOPT
solve test that must agree with the issue's scalar-indexed loop
formulation.

Co-authored-by: Transurgeon <[email protected]>
Co-authored-by: Claude Fable 5 <[email protected]>

* Fix FloatingPointError in Sum shape inference via $O(1)$ calculation (#3238)

* perf: replace memory-allocating shape inference in Sum with O(1) logic

* perf: replace memory-allocating shape inference in Sum with O(1) logic

* style: automated ruff cleanup of imports and formatting

* perf: revert to numpy shape inference with int8 dtype to fix benchmark slowdown

* Replace Sum shape inference with exact tuple formula (#3460)

* Updated the link to the Uno installation instructions (#3461)

* Make bmat promote scalars and 1-D blocks like numpy.block (#3466)

* Make bmat promote scalars and 1-D blocks like numpy.block

cp.bmat raised a ValueError when a block row mixed a scalar or vector
with a matrix (e.g. building an LMI / Schur complement), because the
blocks were stacked without dimension promotion. numpy.block handles
this by promoting every block to a common dimensionality.

Promote scalars to (1, 1) and 1-D blocks to (1, n) before stacking so
bmat matches numpy.block. Adds a regression test.

Fixes #2328

* Test bmat promotion of a scalar Variable

* Add cp.block wishlist entry to contributing docs

* Fix SCIP interface misclassifying single-SOC problems as continuous LPs (#3467)

* fix: SCIP interface misclassifying single-SOC problems as continuous LPs

* Rename SCIP SOCP regression test and tighten socp_1 tolerance to places=3

* Draft: Pass attributes to real variables in Complex2Real canonicalizer (#3342)

* Draft: Pass attributes to real variables in Complex2Real canonicalizer

* Handle complex PSD attributes in Complex2Real

* Fix Complex2Real attribute handling

* Add Complex2Real sparsity attribute tests

* Move Complex2Real test imports to module scope

* Move Complex2Real test imports to module scope

* Clarify Hermitian imaginary attribute handling

* Use allow-list for Complex2Real split attributes

* Use set intersection for complex/imag invalid attr check

* Add Complex2Real Hermitian attribute regression tests

* Vectorize Gurobi conic interface model building, drop pre-9.5 code paths (#3488)

* Begin development of CVXPY 1.10

* Vectorize Gurobi conic interface model building

The conic interface built the Gurobi model one API call at a time, which
made it far slower than the QP interface on the same problem sizes. Three
things dominated:

- Variable types were assigned with `if i in data[s.BOOL_IDX]`, a list
  membership test per column. That is quadratic in the number of integer
  variables: a 6000-boolean MISOCP spent ~0.5s here alone.
- Variables were added one `addVar` call at a time.
- Each SOC cone added its auxiliary variables one at a time and built one
  `LinExpr` per row, addressing the standing "TODO: add all SOC constrs at
  once!".

Variable types are now assigned by index, variables are created with a
single `addMVar`, and all cones are built together: one `addMVar` for the
auxiliary variables, one `addMConstr` for the coupling rows written as
[A | I] @ [x; aux] == b, and one `QuadExpr` per cone. Warm starts and
primal extraction use Gurobi's batch attribute API instead of per-variable
access, and the primal is now trimmed to the n problem variables rather
than carrying the internal SOC auxiliaries.

The scalar code paths are kept for Gurobi older than v9.5, which lacks the
matrix API. Variable naming (x_i, soc_t_i, soc_x_i) is unchanged; the
existing test_gurobi_socp_bound_attr depends on it.

Model construction time, best of 3 (Gurobi 13.0.2):

  SOCP 5000 cones x dim 10     244.0 ms -> 86.2 ms    2.8x
  SOCP 1 cone x dim 50000      154.7 ms -> 59.0 ms    2.6x
  MISOCP 1000 booleans          18.3 ms ->  4.8 ms    3.8x
  MISOCP 3000 booleans         131.0 ms -> 13.9 ms    9.4x
  MISOCP 6000 booleans         498.0 ms -> 28.5 ms   17.5x
  LP 50000 rows                 18.9 ms -> 19.1 ms    1.0x

MISOCP construction is now linear rather than quadratic in the number of
integer variables. LP problems are unchanged, as that path already used
addMConstr.

Co-Authored-By: Claude Opus 5 <[email protected]>

* Fix aliasing in cone test, drop getAttr/setAttr, guard the fast path

Review follow-ups.

The multiple-cones test built its CLARABEL reference from the same
Variable and constraint objects as the Gurobi problem. Solving the
reference overwrote x.value and every dual_value, so the primal
comparison compared x.value against itself and the dual checks inspected
CLARABEL's output rather than Gurobi's. Gurobi's results are now
snapshotted before the reference solve, and the cone duals are compared by
value instead of only by shape, which is what actually pins down their
ordering. Reversing the per-cone iteration order now fails the test; it
passed before. Cross-solver tolerance is 1e-3.

Replaced the getAttr/setAttr calls introduced by the previous commit with
direct attribute access: MVar.fromlist(...).X and MVar.Start on the
vectorized path, plain attribute loops on the scalar fallback. This costs
nothing, since the vectorized path already reads MVar attributes in bulk.
The two remaining getAttr calls extract constraint duals and are unchanged
from master.

Added test_gurobi_conic_uses_vectorized_api, which patches Model.addVar
and Model.addLConstr to raise so a silent regression to the scalar API
fails the suite rather than merely running slowly. addMVar/addMConstr do
not call either internally, so the guard is exact.

Co-Authored-By: Claude Opus 5 <[email protected]>

* Remove Gurobi pre-9.5 code paths from the conic interface

Finishes the job started in #2541, "Remove old Gurobi code paths", which
removed the same three-way addMConstr / addMConstrs / scalar-loop pattern
from the QP interface but only touched that file.

These hasattr checks were the last Gurobi version gates in the codebase,
and they were already unreachable. gurobi_qpif.py calls addMConstr (v9.5+)
and setMObjective (v9+) with no gate at all, so anyone on an older Gurobi
has a broken install for every LP and QP long before they could reach the
conic fallback. Nothing exercised the fallbacks either: CI installs
gurobipy from pip, which is always current.

Drops the `vectorized` flag and both scalar helpers. add_model_lin_constr
and add_model_soc_constr have no callers left; CPLEX and SCIP define their
own copies rather than inheriting these.

The install docs claimed Gurobi 7.5.2, which the QP interface has not
honoured for some time. Updated to 9.5 (released November 2021).

No behavior change on any supported version, and no performance change.

Co-Authored-By: Claude Opus 5 <[email protected]>

* Address review: promote the multi-cone test to a standard SOCP

Replace the Gurobi-only ``test_gurobi_multiple_socp_cones`` with a new
``socp_4`` standard test: three SOC constraints of dimension 2, 3 and 4
plus an equality, with the third cone inactive at the optimum. Expected
objective, primal and per-cone duals are hardcoded rather than
cross-checked against another solver, and the test is wired up for
GUROBI, CLARABEL and SCS.

Also drop the comment describing the removed membership-test loop and
reword the primal-extraction comment.

Co-Authored-By: Claude Opus 5 <[email protected]>

* Run socp_4 for GUROBI only

Keeps this PR to the interface it is about. The other conic interfaces
that split the SOC block themselves can pick up the standard test
separately.

Co-Authored-By: Claude Opus 5 <[email protected]>

---------

Co-authored-by: Claude Opus 5 <[email protected]>

* Fix CVXOPT inequality presolve: opposite-row deletion and dual recovery (#3390)

* Only merge positive-ratio rows in CVXOPT inequality presolve

compress_matrix, used by the CVXOPT interface to drop redundant
inequality rows when the equality matrix is rank deficient, treated a
row as redundant whenever it was any scalar multiple of a kept row with
a matching right-hand side ratio. For inequalities Gx <= h a negative
multiple describes the opposite half-space (e.g. z >= 1 stored as
-z <= -1 versus z <= 1), so dropping it enlarges the feasible set:
min z s.t. z <= 1, z >= 1 plus duplicated equalities was reported
unbounded instead of optimal at 1.

Require a positive ratio before merging rows. compress_matrix's only
caller compresses the LEQ block of G in cvxopt_conif.py, so no equality
usage relies on merging negative multiples.

Tests: unit test that compress_matrix keeps negative-multiple rows, and
a CVXOPT regression test for the formerly-unbounded problem.

Co-Authored-By: Claude Fable 5 <[email protected]>

* Fix dual recovery for inequalities merged by CVXOPT presolve

When CVXOPT's presolve merged duplicate inequality rows, the duals of
the original constraints were recovered as P_leq @ z, giving every
merged row a scaled copy of the full compressed multiplier. Duplicated
constraints then received duals summing to (1 + ratio^2) times the
correct value, violating KKT stationarity (e.g. z >= 1 stated twice
produced duals 1.0 and 1.0 where the sum must be 1.0).

Recover the duals with an expansion matrix that assigns each compressed
dual entirely to the row that was kept and zero to dropped duplicates,
so G_orig.T @ y_orig == G_compr.T @ z and stationarity holds.
compress_matrix now also returns the kept row indices to build this
expansion.

Tests: CVXOPT regression test checking stationarity of the recovered
duals for exact and scaled duplicate inequalities; existing
compress_matrix unit test extended to cover the kept-row indices.

Co-Authored-By: Claude Fable 5 <[email protected]>

* move CVXOPT-only utility into cvxopt_conif.py

move tests of CVXOPT utility from free-function in test_conic_solvers.py to the class test_conic_solvers.py::TestCVXOPT

* Fixes precommit

---------

Co-authored-by: Claude Fable 5 <[email protected]>
Co-authored-by: Riley Murray <[email protected]>

* Fix trace canonicalization for scalar products (#3492)

* Copy aliased arrays in Constant so cached sign cannot go stale (#3404)

* Copy aliased arrays in Constant so cached sign cannot go stale

Constant stores float64 ndarrays (and same-dtype scipy sparse inputs)
by reference: const_to_matrix returns the caller's array unchanged. But
sign/PSD/symmetry attributes are computed lazily and cached on first
query. Mutating the original array afterwards left `.value` live-updated
while is_nonneg()/sign reported the stale cached answer, e.g.:

    a = np.array([1.0, 2.0]); c = cp.Constant(a)
    c.is_nonneg()              # True (cached)
    a[0] = -5.0
    c.value, c.is_nonneg()     # [-5., 2.], True  -- inconsistent

Fix: defensively copy at construction, but only when the converted
value actually shares memory with the caller's input (checked via
np.may_share_memory); conversions that already produced a fresh array
(lists, scalars, dtype casts, sparse-to-dense) pay nothing. This gives
Constant snapshot semantics, matching its lack of a value setter.
Variable/Parameter are unaffected: their sign/PSD attributes are
declared at construction, not inferred from the value, so aliasing
there cannot desynchronize.

Measured overhead (canonicalization via get_problem_data, CLARABEL,
min of 3): three dense 2000x2000 constants 3.351s -> 3.404s (~1.6%);
2000-term small-constant tree 0.272s -> 0.279s (~2.6%, ~7ms). Micro:
one 2000x2000 Constant costs 2.4ms (a 32MB memcpy). Well under the
acceptable threshold, including internal canonicalization paths that
construct Constants from arrays they own.

Adds a regression test covering dense and sparse aliasing.

Co-Authored-By: Claude Fable 5 <[email protected]>

* Preserve constant CSE after defensive copies

---------

Co-authored-by: Claude Fable 5 <[email protected]>

* Fix DPP AssertionError for cumsum of parameter-only expressions (#3398)

* Fix DPP AssertionError for cumsum of parameter-only expressions

Bug: a DPP problem containing cumsum(p) @ x, with p a Parameter,
passed is_dpp() but crashed during canonicalization with a bare
AssertionError (coeff_extractor.py: assert affine_expr.is_dpp()).
The same crash occurred in constraint position.

Root cause: cumsum_canon unconditionally replaced cumsum(X) with a
fresh auxiliary Variable defined via difference constraints, even when
X contained no variables. That turned the parameter-affine factor
cumsum(p) into a Variable, so cumsum(p) @ x became Variable @ Variable,
which is not DPP-representable.

Fix: when the cumsum argument contains no variables (constant or
parameter-affine), keep the expression parameter-affine by applying an
explicit lower-triangular ones operator along the cumsum axis
(I_post kron L kron I_pre on the F-order vectorization, which handles
1-D, 2-D, and N-D arguments and all axis values including None). The
existing variable path is unchanged, so variable arguments still avoid
dense triangular matrices; parameter-only arguments are typically
small, so the O(n^2) operator is an acceptable trade-off.

Audited the other eliminate_pwl canonicalizers that introduce
auxiliary variables (cummax etc.); they are nonlinear atoms whose
epigraph variables remain valid for parameter-only arguments, and no
concrete failure could be reproduced, so they are left unchanged.

Tests: new regression tests in cvxpy/tests/test_dpp.py cover the
objective and constraint positions, parameter re-solve (DPP cache),
and 2-D parameter cumsum along axis 0, 1, and None.

Co-Authored-By: Claude Fable 5 <[email protected]>

* Fix parametric-lhs matmul dedup in SCIPY and COO backends

Reshaping a parametric constant from column to matrix format
deduplicated entries by param_idx alone, keeping a single entry per
parameter slice. That is only correct for broadcast copies; a
parametric expression such as A @ p has several genuine entries per
slice (a column of A multiplies each parameter entry), and
broadcast_to(p, shape) @ x places the same parameter at several
positions. Both were silently stuffed with wrong coefficients.

Dedup now runs only when a broadcast factor was divided out, and keys
on the (param_idx, position) pair so only true broadcast copies are
dropped. This is what made test_cumsum_of_parameter fail on the
SCIPY/COO CI jobs: cumsum_canon's new tril @ param expression is
exactly the multi-entry-per-slice case.

Co-Authored-By: Claude Fable 5 <[email protected]>

* Address cumsum DPP review feedback

Share the sparse cumsum axis operator with gradient construction, centralize parametric-entry deduplication with an explicit identical-copy invariant, correct backend documentation, and cover the optional RUST backend.

---------

Co-authored-by: Claude Fable 5 <[email protected]>

* Fix zero-add fold dropping numpy broadcasting (#3387)

Adding a zero-valued Constant to an expression folded the zero away
unconditionally (expr + 0 optimization from 9a66a4a), so a
larger-shaped zero like np.zeros(3) added to a scalar Variable was
silently discarded instead of broadcasting. This produced wrong shapes
and wrong solve results, e.g. Minimize(sum(x + np.zeros(3))) returned
1.0 instead of 3.0.

Fix: in __add__/__radd__, only fold the zero when its shape broadcasts
to self.shape without enlarging it, i.e.
np.broadcast_shapes(self.shape, other.shape) == self.shape. The
subtraction paths reuse __add__ and are fixed by the same check.

Adds TestExpressions.test_add_zero covering the fold (same shape and
scalar zero), broadcasting for add/sub in both operand orders,
incompatible shapes raising ValueError, and the solve repro.

Co-authored-by: Claude Fable 5 <[email protected]>

* Fix silently dropped quad objectives for dim-reduced attribute variables (#3384)

Quadratic objectives over variables with dim-reducing attributes
(symmetric, PSD, NSD, diag, sparsity) were silently dropped: an
all-zero P matrix was sent to the solver, producing wrong optimal
values, primals, and duals on the quadratic-objective path (CLARABEL,
OSQP, SCS, ...). E.g. minimizing sum_squares(X) with X symmetric and
trace(X) >= 2 returned 4.5 instead of 2.0.

Root cause: the Dcp2Cone quad canonicalizers (quad_form_canon,
power_canon, quad_over_lin_canon) keep the argument inline in
SymbolicQuadForm whenever it is a Variable. CvxAttr2Constr then runs
after Dcp2Cone and tree_copy-substitutes dim-reduced attribute
variables with affine expressions (e.g. reshape(upper_tri_to_full @ u))
inside the SymbolicQuadForm. CoeffExtractor.extract_quadratic_coeffs
keys each P block on SymbolicQuadForm.args[0].id, which is now the id
of a reshape atom rather than a variable, so the block matches no
variable offset and quad_form() silently emits P = 0.

Fix: in the three quad canonicalizers, route Variables with
dim-reducing attributes (Leaf._has_dim_reducing_attr) through the
existing auxiliary-variable path (t = Variable(shape); arg == t), so
the SymbolicQuadForm argument remains a true Variable after
CvxAttr2Constr. This was preferred over teaching the coefficient
extractor to compose the quadratic through an affine map u -> A u,
which would invalidate the diagonal/block structure assumptions of its
constant-P paths (A'PA is dense even for diagonal P) and require
reworking all four extraction paths. Attributes substituted by a plain
Variable (nonneg, nonpos, pos, neg, bounds) were already correct and
are unaffected.

Tests: differential regression tests in test_qp_solvers.py covering
symmetric/PSD/diag/sparsity variables (optimal value and constraint
dual against the analytic optimum, CLARABEL and OSQP) and a DPP
objective mixing a symmetric variable with a plain variable, solved
and re-solved for two parameter values. Verified quad_form,
quad_over_lin, huber, and NSD variants against the conic reference,
and that nonneg/bounds paths are unchanged.

Review follow-up: added coverage for the NSD attribute (with a
sign-compatible constraint), for power_canon via sum(square(X)), and for
quad_form_canon via an explicit quad_form over a sparsity-constrained
vector. All three fail on master.

Co-authored-by: Claude Fable 5 <[email protected]>

* Support broader axis handling in stats atoms (#3496)

(cherry picked from commit 8ad59df)

* Make geo_mean an AxisAtom, reducing along any axis (#3493)

@PTNobel on #3493: one PR, one implementation. Three things came out of
the review and they all point at AxisAtom, so geo_mean now inherits from
it rather than stacking scalars.

axis and keepdims follow cp.mean's spelling, and tuples of axes come
with them: normalize_axis handles tuples, negative indices, and reducing
over every axis. Shapes and values were checked against
exp(mean(log(x), axis)) for 2-D and 3-D inputs, every single axis and
every pair, with and without keepdims.

The stacking is gone. gm was already vectorised — it takes t.size and
emits one SOC over all entries — and PowConeND takes one cone per column,
so both canonicalizers now cover every output with a single set of
constraints. Measured through geo_mean_approx_canon: 3 SOC constraints
for 2 outputs, and 3 for 200. A test pins that down.

Zero-weight entries used to be indexed out of the argument itself, which
put a stray [True, True, True] into every geo_mean's printed form and
would not have survived an axis. They are now dropped in the aligned view
the canonicalizers work on, so the argument keeps the shape the axis
bookkeeping is based on, and the two label expectations lose the stray
index.

This also fixes a row vector reaching the solver: geo_mean(Variable((1, 4)))
built and evaluated fine on master and then raised IndexError during
canonicalization, because the canonicalizer indexed along axis 0.

Closes the geo_mean half of #3489

(cherry picked from commit 1f673ee)

* Run pyright inside the uv environment (ported from #3339)

The pyright workflow ran `uvx pyright` with no environment configured, so
numpy/scipy imports resolved to Unknown. With #3493's geo_mean rewrite this
tripped reportOptionalIterable: `p = p[self._keep]` assigns an Unknown,
which resets p's narrowed type to its declared Optional parameter type.

Port the workflow half of #3339 from master: add the dev dependency group
(pyright, scipy-stubs) and run pyright via `uv run --no-sync` so it
analyzes real NumPy/SciPy types. The full run is at 0 errors on this branch.

Co-Authored-By: Claude Fable 5 <[email protected]>

* Extract direct cone constraints for Moreau (#3499)

Moreau supports placing cones directly on subvectors of `x`, avoiding redundant identity rows in `A x + s = b`. This PR adds an opt-in `ExtractDirectCones` reduction and advertises Moreau's supported direct cone families through `DIR_CONE_KINDS`.

The reduction recognizes structurally constant NonNeg, SOC, SvecPSD, ExpCone, PowCone3D, and PowConeND blocks after matrix stuffing. Registered extractors emit typed `DirectCone` metadata in `ParamConeProg.dir_cones`, and a shared linear scan removes eligible rows without sorting. Parameter-dependent blocks stay on the slack side, preserving DPP correctness.

The Moreau interface maps unscaled symmetric variables into scaled svec coordinates and recovers direct cone duals from the full quadratic KKT residual. Slack-side PSD constraints remain supported, with mixed PSD/EXP/POW3D ordering covered by an interface regression. CI uses the normal `MOREAU` optional dependency, and the solver support table now lists SDP support.

The interface uses `moreau.DirectConeSpec` and `Cones.dir_cones`. The corresponding API rename is [moreau-project/moreau#3](moreau-project/moreau#3) and is included in the planned Moreau >=0.4.0 release. The cone-ordering fix was merged in [moreau-project/moreau#1](moreau-project/moreau#1).

Moreau integration was verified against a CPU wheel built from the companion branch, whose current version metadata is 0.3.3; this is not validation of a published >=0.4.0 package. The companion branch also passes 111 focused CPU tests. CUDA build and runtime validation are outstanding there.

(cherry picked from commit 606f6b0)

* Add regression tests for Python-scalar values on discrete leaves (#3287)

* fix: handle scalar Python bool/int in leaf.py project()
When val is a plain Python bool (False/True) or int, calling
val.astype() or accessing val.ndim fails with AttributeError.
Convert val to numpy array before accessing .shape and .ndim.
Fixes #3284

* fix: set lp.col_names_ before passModel to preserve variable names in HiGHS

* fix: set lp.col_names_ before passModel to preserve variable names in HiGHS

* style: remove trailing whitespace in highs_conif.py

* fix: move col_names_ before passModel and remove duplicate call

* fix: restore TODO comment in solve_via_data

* Add regression tests for Python-scalar values on discrete leaves (gh-3284)

cp.Parameter(boolean=True, value=False) used to raise AttributeError:
project() called val.astype, which a Python bool/int does not have. Master
no longer reproduces it -- project() was since rewritten to go through
np.array/np.shape/np.ndim -- but nothing in the suite pins that, so the
fix can regress silently.

Cover both the constructor, which is where the issue was actually reported,
and project() directly, asserting a 0-d value survives as 0-d.

Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>

---------

Co-authored-by: William Zijie Zhang <[email protected]>
Co-authored-by: Claude Opus 5 (1M context) <[email protected]>
(cherry picked from commit e7e704c)

* Fix np.matmul semantics for 1-D operands in batched matmul (#3386)

* Fix np.matmul semantics for 1-D operands in batched matmul

Bug: in N-D (batched) matmul, a 1-D operand violated np.matmul
semantics. w (3,) @ A (2,3,4) silently returned shape (2,2,4) with
wrong values, and A (2,3,4) @ v (4,) raised ValueError; numpy gives
(2,4) and (2,3) respectively.

Root cause: MulExpression._broadcast_batch_dims treated a 1-D operand
as a matrix with empty batch dims and broadcast it along the batch
dimensions, instead of applying numpy's documented promotion rule
(prepend/append a 1, matmul, then drop the inserted axis).

Fix: _broadcast_batch_dims now promotes a 1-D operand to 2-D via
reshape ((k,) -> (1,k) for the lhs, (k,) -> (k,1) for the rhs) before
batch broadcasting, and the matmul() entry point (now used by
__matmul__/__rmatmul__) removes the inserted axis from the result.
Shapes, values, and the canonicalization/solve path all match numpy.

Tests: added TestNDRmulEdgeCases::test_1d_operand_batched covering
constant values for (3,)@(2,3,4), (2,3,4)@(4,), and 4-D batches, plus
solve-path checks with 1-D variables on either side, for the SCIPY and
COO backends.

Co-Authored-By: Claude Fable 5 <[email protected]>

* Update cvxpy/atoms/affine/binary_operators.py

* Apply suggestion from @PTNobel

* Address batched 1-D matmul review feedback

Route deprecated star-matmul through the public wrapper, derive result shapes from mul_shapes, preserve user-facing dimension errors, materialize concrete vector operands, and expand variable and parameter coverage.

* Clarify the matmul dim pre-check and drop _supports_cpp from the test

Explain why matmul() checks inner dimensions itself: mul_shapes raises
on a mismatch too, but only after 1-D promotion, so it reports the
promoted (1, k) instead of the (k,) the caller passed.

Replace the test's `assert problem._supports_cpp()` with a direct
assertion on the expression tree: a concrete 1-D operand leaves no
broadcast_to atom, while a symbolic one (which must stay symbolic for
DPP) keeps it. This observes the materialization directly instead of
going through a C++-backend predicate that says nothing about the
optimization once cvxcore is gone.

Co-Authored-By: Claude Opus 5 (1M context) <[email protected]>

---------

Co-authored-by: Claude Fable 5 <[email protected]>
Co-authored-by: William Zijie Zhang <[email protected]>
(cherry picked from commit f31996f)

* Handle Python scalars in the discrete-leaf backport

The regression tests from #3287 require scalar-safe projection, which master already obtained through its later Leaf.project rewrite. Port the np.array, np.shape, and np.ndim conversions to the 1.9.x boolean and integer paths.

* ChatGPT generated: Backport Moreau solver reuse and native duals (#3532)

(cherry picked from commit 50fcc33)

* ChatGPT generated: Backport SCS 3.3 compatibility fixes (#3513)

(cherry picked from commit fb34106)

* ChatGPT generated: Backport scalar sum axis compatibility (#3512)

(cherry picked from commit a4b9e10)

* ChatGPT generated: Backport matrix norm axis validation (#3501)

(cherry picked from commit 2d29f48)

* ChatGPT generated: Backport CI matrix exclusion fix (#3507)

(cherry picked from commit cb7e8b6)

---------

Signed-off-by: Ramakrishna Prabhu <[email protected]>
Co-authored-by: Claude Opus 4.8 (1M context) <[email protected]>
Co-authored-by: Parth Nobel <[email protected]>
Co-authored-by: Ramakrishna Prabhu <[email protected]>
Co-authored-by: Transurgeon <[email protected]>
Co-authored-by: Aditya Karna <[email protected]>
Co-authored-by: Charlie Vanaret <[email protected]>
Co-authored-by: Enzo Boulin <[email protected]>
Co-authored-by: Hugo Delatte <[email protected]>
Co-authored-by: Steven Diamond <[email protected]>
Co-authored-by: Riley Murray <[email protected]>
Co-authored-by: mrlonely <[email protected]>
Co-authored-by: Nicoló Angileri <[email protected]>
Co-authored-by: DHRUV <[email protected]>

* Support tuple axes in std

@PTNobel on #3489, item 2, and following your "happy to do different
behavior on std and norm even if that adds complexity to std": this is
the std half, and it does not touch norm.

A tuple of axes pools the entries of every axis it names into one sample
set per remaining position. norm takes a single axis, so the pooled axes
are moved to the front with moveaxis, flattened into one in Fortran
order, reduced with norm(..., 2, axis=0), and the result folded back into
the shape the pooled axes left behind. That is the same alignment
geo_mean uses since #3493, order included.

std is the one atom here that cannot borrow var's route: var is
sum_squares over the axes and takes tuples already, but sqrt of a convex
expression is not DCP, which is why std goes through norm in the first
place.

Checked against numpy over (0,1), (1,2), (0,2), the one-element tuple
(0,), the full (0,1,2) and the negative (-3,-1), each crossed with
keepdims and ddof in {0, 1}: shapes and values all agree. The int and
None paths are untouched and still match. The tuple form is DCP and
solves.

Note this drops an existing assertion in test_stats that
cp.std(a, axis=(0, 1)) raises ValueError, which was pinning the gap this
closes.

No new backend fallback: std(X, axis=0) on a three-dimensional variable
already falls back to SCIPY because the CPP backend does not take
expressions past two dimensions, so the reshape adds nothing there.

* Reuse the reduction output shape in std