From 59ce6bde8d71fc9f0babc5817167c88efa0b8999 Mon Sep 17 00:00:00 2001 From: 40% Date: Fri, 16 Jan 2026 15:21:40 +0800 Subject: [PATCH 01/12] Speed up `np.add.reduce` Implements the __array_ufunc__ method in MatrixExpr to handle numpy ufuncs, specifically enabling correct behavior for reductions like np.add.reduce by delegating to the sum method. This improves compatibility with numpy operations. --- src/pyscipopt/matrix.pxi | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/pyscipopt/matrix.pxi b/src/pyscipopt/matrix.pxi index 1a6a09cf3..d2e9a0ed3 100644 --- a/src/pyscipopt/matrix.pxi +++ b/src/pyscipopt/matrix.pxi @@ -51,6 +51,13 @@ def _matrixexpr_richcmp(self, other, op): class MatrixExpr(np.ndarray): + def __array_ufunc__(self, ufunc, method, *args, **kwargs): + if method == "reduce": + if ufunc is np.add and isinstance(args[0], MatrixExpr): + return args[0].sum(**kwargs) + + return super().__array_ufunc__(ufunc, method, *args, **kwargs) + def sum( self, axis: Optional[Union[int, Tuple[int, ...]]] = None, From 1c66a5c861a23f60938c98e6d18c10b7c55a44ae Mon Sep 17 00:00:00 2001 From: 40% Date: Fri, 16 Jan 2026 16:01:13 +0800 Subject: [PATCH 02/12] BUG: `MatrixExpr.mean(axis=1)` will crush kernel Moved the sum computation logic from MatrixExpr.sum and __array_ufunc__ into a new _core_sum function for better code reuse and maintainability. --- src/pyscipopt/matrix.pxi | 51 +++++++++++++++++++++++----------------- 1 file changed, 30 insertions(+), 21 deletions(-) diff --git a/src/pyscipopt/matrix.pxi b/src/pyscipopt/matrix.pxi index d2e9a0ed3..723072bb1 100644 --- a/src/pyscipopt/matrix.pxi +++ b/src/pyscipopt/matrix.pxi @@ -54,7 +54,7 @@ class MatrixExpr(np.ndarray): def __array_ufunc__(self, ufunc, method, *args, **kwargs): if method == "reduce": if ufunc is np.add and isinstance(args[0], MatrixExpr): - return args[0].sum(**kwargs) + return _core_sum(args[0], **kwargs) return super().__array_ufunc__(ufunc, method, *args, **kwargs) @@ -92,26 +92,7 @@ class MatrixExpr(np.ndarray): a MatrixExpr. """ - axis: Tuple[int, ...] = normalize_axis_tuple( - range(self.ndim) if axis is None else axis, self.ndim - ) - if len(axis) == self.ndim: - res = quicksum(self.flat) - return ( - np.array([res], dtype=object).reshape([1] * self.ndim).view(MatrixExpr) - if keepdims - else res - ) - - keep_axes = tuple(i for i in range(self.ndim) if i not in axis) - shape = ( - tuple(1 if i in axis else self.shape[i] for i in range(self.ndim)) - if keepdims - else tuple(self.shape[i] for i in keep_axes) - ) - return np.apply_along_axis( - quicksum, -1, self.transpose(keep_axes + axis).reshape(shape + (-1,)) - ).view(MatrixExpr) + return _core_sum(self, axis=axis, keepdims=keepdims, **kwargs) def __le__(self, other: Union[float, int, "Expr", np.ndarray, "MatrixExpr"]) -> MatrixExprCons: return _matrixexpr_richcmp(self, other, 1) @@ -168,3 +149,31 @@ class MatrixExprCons(np.ndarray): def __eq__(self, other): raise NotImplementedError("Cannot compare MatrixExprCons with '=='.") + + +def _core_sum( + a: MatrixExpr, + axis: Optional[Union[int, Tuple[int, ...]]] = None, + keepdims: bool = False, + **kwargs, +) -> Union[Expr, MatrixExpr]: + axis: Tuple[int, ...] = normalize_axis_tuple( + range(a.ndim) if axis is None else axis, a.ndim + ) + if len(axis) == a.ndim: + res = quicksum(a.flat) + return ( + np.array([res], dtype=object).reshape([1] * a.ndim).view(MatrixExpr) + if keepdims + else res + ) + + keep_axes = tuple(i for i in range(a.ndim) if i not in axis) + shape = ( + tuple(1 if i in axis else a.shape[i] for i in range(a.ndim)) + if keepdims + else tuple(a.shape[i] for i in keep_axes) + ) + return np.apply_along_axis( + quicksum, -1, a.transpose(keep_axes + axis).reshape(shape + (-1,)) + ).view(MatrixExpr) From a2a57c1f494fee4b815bcc78a9e30218e12d8acc Mon Sep 17 00:00:00 2001 From: 40% Date: Fri, 16 Jan 2026 16:02:10 +0800 Subject: [PATCH 03/12] Add tests for matrix mean performance and type Introduces tests to compare the performance of the mean operation on matrix variables and checks the return types of mean with and without axis argument. --- tests/test_matrix_variable.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/tests/test_matrix_variable.py b/tests/test_matrix_variable.py index e4758f077..071132207 100644 --- a/tests/test_matrix_variable.py +++ b/tests/test_matrix_variable.py @@ -293,6 +293,30 @@ def test_matrix_sum_axis_not_none_performance(n): assert model.isGT(end_orig - start_orig, end_matrix - start_matrix) +@pytest.mark.parametrize("n", [50, 100]) +def test_matrix_mean_performance(n): + model = Model() + x = model.addMatrixVar((n, n)) + + start = time() + np.ndarray.mean(x, axis=0) + orig = time() - start + + start = time() + x.mean(axis=0) + matrix = time() - start + + assert model.isGT(orig, matrix) + + +def test_matrix_mean(): + model = Model() + x = model.addMatrixVar((2, 2)) + + assert isinstance(x.mean(), Expr) + assert isinstance(x.mean(1), MatrixExpr) + + def test_add_cons_matrixVar(): m = Model() matrix_variable = m.addMatrixVar(shape=(3, 3), vtype="B", name="A", obj=1) From 7804c1007ef3eaea146d3d625ba38082228ccf88 Mon Sep 17 00:00:00 2001 From: 40% Date: Fri, 16 Jan 2026 16:05:29 +0800 Subject: [PATCH 04/12] Update changelog for MatrixExpr.add.reduce optimization Documented the speed improvement for MatrixExpr.add.reduce using quicksum in the changelog. --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 52df46228..077d56b12 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ - Fixed lotsizing_lazy example ### Changed - changed default value of enablepricing flag to True +- Speed up MatrixExpr.add.reduce via quicksum ### Removed ## 6.0.0 - 2025.xx.yy From 4a8a23349be168c2fbe7bad0809b5c4c448ac025 Mon Sep 17 00:00:00 2001 From: 40% Date: Fri, 16 Jan 2026 16:08:04 +0800 Subject: [PATCH 05/12] Remove MatrixExpr.sum method and update _core_sum docstring The sum method was removed from the MatrixExpr class, consolidating summation logic in the _core_sum function. The docstring for _core_sum was expanded to include detailed parameter and return value descriptions, improving code clarity and maintainability. --- src/pyscipopt/matrix.pxi | 66 ++++++++++++++++++---------------------- 1 file changed, 30 insertions(+), 36 deletions(-) diff --git a/src/pyscipopt/matrix.pxi b/src/pyscipopt/matrix.pxi index 723072bb1..adb6f10e7 100644 --- a/src/pyscipopt/matrix.pxi +++ b/src/pyscipopt/matrix.pxi @@ -58,42 +58,6 @@ class MatrixExpr(np.ndarray): return super().__array_ufunc__(ufunc, method, *args, **kwargs) - def sum( - self, - axis: Optional[Union[int, Tuple[int, ...]]] = None, - keepdims: bool = False, - **kwargs, - ) -> Union[Expr, MatrixExpr]: - """ - Return the sum of the array elements over the given axis. - - Parameters - ---------- - axis : None or int or tuple of ints, optional - Axis or axes along which a sum is performed. The default, axis=None, will - sum all of the elements of the input array. If axis is negative it counts - from the last to the first axis. If axis is a tuple of ints, a sum is - performed on all of the axes specified in the tuple instead of a single axis - or all the axes as before. - - keepdims : bool, optional - If this is set to True, the axes which are reduced are left in the result as - dimensions with size one. With this option, the result will broadcast - correctly against the input array. - - **kwargs : ignored - Additional keyword arguments are ignored. They exist for compatibility - with `numpy.ndarray.sum`. - - Returns - ------- - Expr or MatrixExpr - If the sum is performed over all axes, return an Expr, otherwise return - a MatrixExpr. - - """ - return _core_sum(self, axis=axis, keepdims=keepdims, **kwargs) - def __le__(self, other: Union[float, int, "Expr", np.ndarray, "MatrixExpr"]) -> MatrixExprCons: return _matrixexpr_richcmp(self, other, 1) @@ -157,6 +121,36 @@ def _core_sum( keepdims: bool = False, **kwargs, ) -> Union[Expr, MatrixExpr]: + """ + Return the sum of the array elements over the given axis. + + Parameters + ---------- + a : MatrixExpr + + axis : None or int or tuple of ints, optional + Axis or axes along which a sum is performed. The default, axis=None, will + sum all of the elements of the input array. If axis is negative it counts + from the last to the first axis. If axis is a tuple of ints, a sum is + performed on all of the axes specified in the tuple instead of a single axis + or all the axes as before. + + keepdims : bool, optional + If this is set to True, the axes which are reduced are left in the result as + dimensions with size one. With this option, the result will broadcast + correctly against the input array. + + **kwargs : ignored + Additional keyword arguments are ignored. They exist for compatibility + with `numpy.ndarray.sum`. + + Returns + ------- + Expr or MatrixExpr + If the sum is performed over all axes, return an Expr, otherwise return + a MatrixExpr. + + """ axis: Tuple[int, ...] = normalize_axis_tuple( range(a.ndim) if axis is None else axis, a.ndim ) From 377ec6a96dcac1a2a91881edd74d36967d3b05f8 Mon Sep 17 00:00:00 2001 From: 40% Date: Fri, 16 Jan 2026 16:17:49 +0800 Subject: [PATCH 06/12] Remove sum method from MatrixExpr stub Deleted the type stub for the sum method in MatrixExpr, likely to reflect changes in the underlying implementation or to correct type information. --- src/pyscipopt/scip.pyi | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/pyscipopt/scip.pyi b/src/pyscipopt/scip.pyi index 831dd02ed..1e703c006 100644 --- a/src/pyscipopt/scip.pyi +++ b/src/pyscipopt/scip.pyi @@ -509,9 +509,6 @@ class MatrixConstraint(numpy.ndarray): def isStickingAtNode(self) -> Incomplete: ... class MatrixExpr(numpy.ndarray): - def sum( # type: ignore[override] - self, axis: Incomplete = ..., keepdims: Incomplete = ..., **kwargs: Incomplete - ) -> Incomplete: ... def __add__(self, other: Incomplete) -> Incomplete: ... def __eq__(self, other: Incomplete)-> Incomplete: ... def __ge__(self, other: Incomplete) -> MatrixExprCons: ... From 1f406a4055ab9ab7b7408789869cfd490621a673 Mon Sep 17 00:00:00 2001 From: 40% Date: Fri, 16 Jan 2026 22:09:30 +0800 Subject: [PATCH 07/12] Improve MatrixExpr ufunc handling and remove unused include Enhanced the __array_ufunc__ method in MatrixExpr to ensure proper array conversion and consistent return types. Added the _ensure_array helper for argument handling. Also removed an unused include of matrix.pxi from expr.pxi. --- src/pyscipopt/expr.pxi | 1 - src/pyscipopt/matrix.pxi | 19 ++++++++++++++++++- 2 files changed, 18 insertions(+), 2 deletions(-) diff --git a/src/pyscipopt/expr.pxi b/src/pyscipopt/expr.pxi index f0c406fcb..49189bc27 100644 --- a/src/pyscipopt/expr.pxi +++ b/src/pyscipopt/expr.pxi @@ -42,7 +42,6 @@ # which should, in princple, modify the expr. However, since we do not implement __isub__, __sub__ # gets called (I guess) and so a copy is returned. # Modifying the expression directly would be a bug, given that the expression might be re-used by the user. -include "matrix.pxi" def _is_number(e): diff --git a/src/pyscipopt/matrix.pxi b/src/pyscipopt/matrix.pxi index adb6f10e7..d2c87aca1 100644 --- a/src/pyscipopt/matrix.pxi +++ b/src/pyscipopt/matrix.pxi @@ -56,7 +56,11 @@ class MatrixExpr(np.ndarray): if ufunc is np.add and isinstance(args[0], MatrixExpr): return _core_sum(args[0], **kwargs) - return super().__array_ufunc__(ufunc, method, *args, **kwargs) + args = _ensure_array(args, convert_scalar=True) + if "out" in kwargs: + kwargs["out"] = _ensure_array(kwargs["out"]) + res = super().__array_ufunc__(ufunc, method, *args, **kwargs) + return res.view(MatrixExpr) if isinstance(res, np.ndarray) else res def __le__(self, other: Union[float, int, "Expr", np.ndarray, "MatrixExpr"]) -> MatrixExprCons: return _matrixexpr_richcmp(self, other, 1) @@ -115,6 +119,18 @@ class MatrixExprCons(np.ndarray): raise NotImplementedError("Cannot compare MatrixExprCons with '=='.") +cdef inline tuple _ensure_array(tuple args, bool convert_scalar = False): + if not convert_scalar: + return tuple( + x.view(np.ndarray) if isinstance(x, np.ndarray) else x + for x in args + ) + return tuple( + x.view(np.ndarray) if isinstance(x, np.ndarray) else np.array(x, dtype=object) + for x in args + ) + + def _core_sum( a: MatrixExpr, axis: Optional[Union[int, Tuple[int, ...]]] = None, @@ -171,3 +187,4 @@ def _core_sum( return np.apply_along_axis( quicksum, -1, a.transpose(keep_axes + axis).reshape(shape + (-1,)) ).view(MatrixExpr) + From a933366e526663cc4fd49b7fdc211d66a51e6fc7 Mon Sep 17 00:00:00 2001 From: 40% Date: Fri, 16 Jan 2026 22:21:48 +0800 Subject: [PATCH 08/12] Fix MatrixExpr matmul return type and update tests Updated MatrixExpr.__matmul__ to return the correct type when the result is not an ndarray. Adjusted tests to reflect the expected return type for 1D matrix multiplication and improved performance test timing logic. --- src/pyscipopt/matrix.pxi | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/pyscipopt/matrix.pxi b/src/pyscipopt/matrix.pxi index d2c87aca1..c083d7034 100644 --- a/src/pyscipopt/matrix.pxi +++ b/src/pyscipopt/matrix.pxi @@ -102,7 +102,8 @@ class MatrixExpr(np.ndarray): return super().__rsub__(other).view(MatrixExpr) def __matmul__(self, other): - return super().__matmul__(other).view(MatrixExpr) + res = super().__matmul__(other) + return res.view(MatrixExpr) if isinstance(res, np.ndarray) else res class MatrixGenExpr(MatrixExpr): pass From 5cc69d2177ff1d3497d7e060c2188b6b5d72dd50 Mon Sep 17 00:00:00 2001 From: 40% Date: Fri, 16 Jan 2026 22:22:23 +0800 Subject: [PATCH 09/12] Refactor matrix sum performance tests timing logic Simplifies timing measurement in matrix sum performance tests by directly calculating elapsed time instead of storing start and end times separately. This improves code readability and reduces variable usage. --- tests/test_matrix_variable.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/tests/test_matrix_variable.py b/tests/test_matrix_variable.py index 071132207..60dd4ce9e 100644 --- a/tests/test_matrix_variable.py +++ b/tests/test_matrix_variable.py @@ -263,16 +263,16 @@ def test_matrix_sum_axis_is_none_performance(n): x = model.addMatrixVar((n, n)) # Original sum via `np.ndarray.sum`, `np.sum` will call subclass method - start_orig = time() + start = time() np.ndarray.sum(x) - end_orig = time() + orig = time() - start # Optimized sum via `quicksum` - start_matrix = time() + start = time() x.sum() - end_matrix = time() + matrix = time() - start - assert model.isGT(end_orig - start_orig, end_matrix - start_matrix) + assert model.isGT(orig, matrix) @pytest.mark.parametrize("n", [50, 100]) @@ -281,16 +281,16 @@ def test_matrix_sum_axis_not_none_performance(n): x = model.addMatrixVar((n, n)) # Original sum via `np.ndarray.sum`, `np.sum` will call subclass method - start_orig = time() + start = time() np.ndarray.sum(x, axis=0) - end_orig = time() + orig = time() - start # Optimized sum via `quicksum` - start_matrix = time() + start = time() x.sum(axis=0) - end_matrix = time() + matrix = time() - start - assert model.isGT(end_orig - start_orig, end_matrix - start_matrix) + assert model.isGT(orig, matrix) @pytest.mark.parametrize("n", [50, 100]) From efc36b04067e7f45ac50a1ac3e9494a20da2f4d1 Mon Sep 17 00:00:00 2001 From: 40% Date: Fri, 16 Jan 2026 22:22:31 +0800 Subject: [PATCH 10/12] Update matmul return type assertion in test Changed the expected type of 1D @ 1D matrix multiplication from MatrixExpr to Expr in test_matrix_matmul_return_type to reflect updated behavior. --- tests/test_matrix_variable.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_matrix_variable.py b/tests/test_matrix_variable.py index 60dd4ce9e..49740d3a4 100644 --- a/tests/test_matrix_variable.py +++ b/tests/test_matrix_variable.py @@ -598,7 +598,7 @@ def test_matrix_matmul_return_type(): # test 1D @ 1D → 0D x = m.addMatrixVar(3) - assert type(x @ x) is MatrixExpr + assert type(x @ x) is Expr # test 1D @ 1D → 2D assert type(x[:, None] @ x[None, :]) is MatrixExpr From b0d190a9704698f71fe76b468091acb80f632968 Mon Sep 17 00:00:00 2001 From: 40% Date: Fri, 16 Jan 2026 22:33:26 +0800 Subject: [PATCH 11/12] Refactor matrix variable tests to use view casting Updated tests to use x.view(MatrixExpr) and x.view(np.ndarray) instead of direct subclass method calls. This clarifies the intent and ensures the correct method resolution for sum and mean operations in performance and result comparison tests. --- tests/test_matrix_variable.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/tests/test_matrix_variable.py b/tests/test_matrix_variable.py index 49740d3a4..7916b4d10 100644 --- a/tests/test_matrix_variable.py +++ b/tests/test_matrix_variable.py @@ -248,11 +248,11 @@ def test_matrix_sum_axis(): ) def test_matrix_sum_result(axis, keepdims): # directly compare the result of np.sum and MatrixExpr.sum - _getVal = np.vectorize(lambda e: e.terms[CONST]) + _getVal = np.vectorize(lambda e: e[CONST]) a = np.arange(6).reshape((1, 2, 3)) np_res = a.sum(axis, keepdims=keepdims) - scip_res = MatrixExpr.sum(a, axis, keepdims=keepdims) + scip_res = a.view(MatrixExpr).sum(axis, keepdims=keepdims) assert (np_res == _getVal(scip_res)).all() assert np_res.shape == _getVal(scip_res).shape @@ -262,9 +262,9 @@ def test_matrix_sum_axis_is_none_performance(n): model = Model() x = model.addMatrixVar((n, n)) - # Original sum via `np.ndarray.sum`, `np.sum` will call subclass method + # Original sum via `np.ndarray.sum` start = time() - np.ndarray.sum(x) + x.view(np.ndarray).sum() orig = time() - start # Optimized sum via `quicksum` @@ -280,9 +280,9 @@ def test_matrix_sum_axis_not_none_performance(n): model = Model() x = model.addMatrixVar((n, n)) - # Original sum via `np.ndarray.sum`, `np.sum` will call subclass method + # Original sum via `np.ndarray.sum` start = time() - np.ndarray.sum(x, axis=0) + x.view(np.ndarray).sum(axis=0) orig = time() - start # Optimized sum via `quicksum` @@ -298,10 +298,12 @@ def test_matrix_mean_performance(n): model = Model() x = model.addMatrixVar((n, n)) + # Original sum via `np.ndarray.sum` start = time() - np.ndarray.mean(x, axis=0) + x.view(np.ndarray).mean(axis=0) orig = time() - start + # Optimized sum via `quicksum` start = time() x.mean(axis=0) matrix = time() - start From 103a96f11985bc9a78c321fa4c75c4c995d1ef35 Mon Sep 17 00:00:00 2001 From: 40% Date: Sat, 17 Jan 2026 00:23:51 +0800 Subject: [PATCH 12/12] define the variable --- src/pyscipopt/matrix.pxi | 1 + 1 file changed, 1 insertion(+) diff --git a/src/pyscipopt/matrix.pxi b/src/pyscipopt/matrix.pxi index c083d7034..5698d53a5 100644 --- a/src/pyscipopt/matrix.pxi +++ b/src/pyscipopt/matrix.pxi @@ -121,6 +121,7 @@ class MatrixExprCons(np.ndarray): cdef inline tuple _ensure_array(tuple args, bool convert_scalar = False): + cdef object x if not convert_scalar: return tuple( x.view(np.ndarray) if isinstance(x, np.ndarray) else x