Skip to content

Fix: ensure step sizes and tau are passed as scalars - #263

Merged
mrava87 merged 4 commits into
PyLops:devfrom
mrava87:fix_l2_scalar_times_op
Aug 29, 2026
Merged

Fix: ensure step sizes and tau are passed as scalars#263
mrava87 merged 4 commits into
PyLops:devfrom
mrava87:fix_l2_scalar_times_op

Conversation

@mrava87

@mrava87 mrava87 commented Aug 28, 2026

Copy link
Copy Markdown
Contributor

This PR fixes a bug identified when running various solvers in GPU/CuPy mode. Whilst in NumPy mode it is possible to multiply a 1-element np.ndarray (aka scalar array) to an operator, this is not possible in CuPy mode due to the internal CuPy mul failing to cover this scenario:

cupy/_core/core.pyx in cupy._core.core._ndarray_base.__mul__()

TypeError: operand type(s) all returned NotImplemented from __array_ufunc__(<ufunc 'multiply'>, '__call__', array(1.+1.j), <8x8 BlockDiag with dtype=float64>): 'ndarray', 'BlockDiag'

This is fixed both at solver level (ensuring .item() is calling when extracting a value for tau or step size from an iteration-dependent array) and in L2 using math.sqrt instead of np.sqrt as the latter repromotes a 1d array also when the input is a float.

@codacy-production

Copy link
Copy Markdown

Up to standards ✅

🟢 Issues 0 issues

Results:
0 new issues

View in Codacy

🟢 Metrics 0 duplication

Metric Results
Duplication 0

View in Codacy

NEW Get contextual insights on your PRs based on Codacy's metrics, along with PR and Jira context, without leaving GitHub. Enable AI reviewer
TIP This summary will be updated as you push new changes.

@mrava87

mrava87 commented Aug 28, 2026

Copy link
Copy Markdown
Contributor Author

@claude review

@claude

claude Bot commented Aug 28, 2026

Copy link
Copy Markdown

Claude finished @mrava87's task in 3m 5s —— View job


Claude finished reviewing

  • Gather context (PR diff, description, changed files)
  • Review pyproximal/optimization/cls_primal.py
  • Review pyproximal/optimization/cls_primaldual.py
  • Review pyproximal/proximal/L2.py
  • Post review feedback

Summary

The core fix is correct: math.sqrt instead of np.sqrt in L2.prox (pyproximal/proximal/L2.py:271-272) avoids producing a CuPy 0-d array before multiplying with the LinearOperator, which is exactly what triggers 'ndarray', 'BlockDiag' TypeError reported in the PR description. The .item() calls added in the solvers are a good complementary defensive fix, converting per-iteration array lookups (self.epsg[self.iiter], self.tau[self.iiter], self.mu[self.iiter]) into plain Python scalars.

Findings

1. The constant-scalar (ndim == 0) branch is not converted to a Python scalar, so the fix is incomplete for the common case of a fixed (non-per-iteration) step size

In all four places touched by this PR, the fix is only applied in the else branch (per-iteration array), while the ndim == 0 branch just reassigns the stored array directly, leaving it as a 0-d array (CuPy 0-d array in GPU mode) instead of a plain scalar:

  • pyproximal/optimization/cls_primal.py:624-630 (ProximalGradient.step)
  • pyproximal/optimization/cls_primal.py:1078-1084 (AndersonProximalGradient.step)
  • pyproximal/optimization/cls_primal.py:1915-1918 (HQS.step)
  • pyproximal/optimization/cls_primaldual.py:296-306 (PrimalDual.step)

e.g.

if self.tau.ndim == 0:
    tau = self.tau                      # still a 0-d (cupy) array
else:
    tau = self.tau[self.iiter].item()   # python scalar

Since self.tau/self.epsg/self.mu are stored via self.ncp.asarray(...) (cls_primaldual.py:215, cls_primal.py:563/1017/1856), passing a single constant value (the most common usage — a fixed step size) still leaves tau as a 0-d CuPy array in the ndim == 0 branch.

This matters because at least one other proximal operator multiplies tau directly against a LinearOperator without going through sqrt/float():

  • pyproximal/proximal/Quadratic.py:111: Op1 = Identity(self.Op.shape[0], dtype=self.Op.dtype) + tau * self.Op

If tau arrives here as a 0-d CuPy array (constant step size case, GPU mode), this would hit the same TypeError this PR is fixing. Suggest normalizing both branches consistently, e.g.:

tau = self.tau.item() if self.tau.ndim == 0 else self.tau[self.iiter].item()

Fix this →

2. (Minor, unconfirmed) GeneralizedProximalGradient.step has a similar pattern that wasn't touched

pyproximal/optimization/cls_primal.py:1528: ztmp = proxg.prox(ztmp, self.tau * self.epsg[i] / self.weights[i]). Here self.epsg/self.weights are also stored via self.ncp.asarray(...) (per-operator, not per-iteration), so self.epsg[i] is a 0-d array element on GPU. The resulting tau scalar passed to proxg.prox could be a 0-d CuPy array for the same reason as above. Worth double-checking if this class is exercised in GPU mode; if so it may need the same .item() treatment.

Nitpicks

  • Good catch on the sigma docstring fix (intfloat) in L2.py:41.
  • The .item() calls will trigger a device-to-host sync each call on GPU, but this is unavoidable to get a plain Python scalar and is consistent with how tau is already used elsewhere (e.g. _check_tau compares tau <= 0).

I was not able to run the test suite or a GPU/CuPy environment in this sandbox to confirm the failure mode directly (no cupy tests exist in pytests/), so finding #1 is based on static analysis of how tau propagates into Quadratic.prox — worth confirming with an actual CuPy repro similar to the one described in the PR body.

@mrava87

mrava87 commented Aug 28, 2026

Copy link
Copy Markdown
Contributor Author

@claude review again

@claude

claude Bot commented Aug 28, 2026

Copy link
Copy Markdown

Claude finished @mrava87's task in 2m 38s —— View job


Todo List

  • Gather context (PR diff, previous review findings)
  • Re-check finding Welcome to pyproximal Discussions! #1 (ndim == 0 branches not calling .item())
  • Re-check finding Prepare for v0.0.0 #2 (GeneralizedProximalGradient.step)
  • Check Quadratic.py line 111 concern against the new commit
  • Post updated review

Summary

Commit ab2e7af ("minor: added more .item()") addresses both issues raised in the previous review:

  • The ndim == 0 branches for tau/epsg/mu now call .item() consistently with the else branch, in ProximalGradient.step (pyproximal/optimization/cls_primal.py:625-627), AndersonProximalGradient.step (:1079-1081), HQS.step (:1917-1918), and PrimalDual.step (pyproximal/optimization/cls_primaldual.py:297-300, 303-306).
  • GeneralizedProximalGradient.step now also calls .item() on self.epsg[i] and self.weights[i] (pyproximal/optimization/cls_primal.py:1527-1531).

Good follow-up — thanks for addressing both.

New finding

ProximalGradient/AndersonProximalGradient still pass a non-scalar tau into proxg.prox(), so the original GPU bug can still be triggered through these (very commonly used) solvers

Unlike HQS/PrimalDual, self.tau in these two classes is never reduced to a Python scalar — it's stored as self.ncp.atleast_1d(self.ncp.asarray(tau, ...)) (pyproximal/optimization/cls_primal.py:563, :1017) and then used directly:

  • pyproximal/optimization/cls_primal.py:635, 639-640, 657-658: epsg * self.tau (passed as the tau arg to proxg.prox) is a Python float times a 1-element ncp array → still an array, not a scalar.
  • pyproximal/optimization/cls_primal.py:1022, 1111, 1117, 1128 (AndersonProximalGradient): same pattern, epsg_ * self.tau / epsg * self.tau.
  • The backtracking path is also affected: cast(float, self.tau) at pyproximal/optimization/cls_primal.py:647 is a typing.cast — a no-op at runtime — so tau inside _backtracking (cls_primal.py:47-78) remains the ncp array, and self.tau = tau (line 645) keeps it an array across all subsequent iterations too.

If proxg is e.g. Quadratic (Op1 = Identity(...) + tau * self.Op, pyproximal/proximal/Quadratic.py:111) and the solver is run on CuPy arrays, tau here will still be a 1-element CuPy array multiplying a LinearOperator, which is exactly the TypeError this PR is meant to fix — just reached through ProximalGradient/AndersonProximalGradient instead of HQS/PrimalDual. Since ProximalGradient is the most widely used solver, this is arguably the most important path to cover.

Suggested fix: extract a scalar tau at the top of ProximalGradient.step/AndersonProximalGradient.step the same way epsg is already extracted (e.g. tau = self.tau.item()), use that scalar in the arithmetic, and store self.tau back as a plain float after backtracking too (self.tau = float(tau) in _backtracking, or convert after the call).

Fix this →

Nitpicks

  • The math.sqrt-in-L2.prox and solver-level .item() fixes look correct and match the reported failure mode for the classes they touch (HQS, PrimalDual, GeneralizedProximalGradient).

I was not able to run the test suite or a GPU/CuPy environment in this sandbox to reproduce the failure directly; the finding above is based on static tracing of how tau flows from ProximalGradient.setup/step into proxg.prox() and then into Quadratic.prox. Worth confirming with a CuPy repro similar to the one in the PR description, using ProximalGradient with a Quadratic proximal operator.

@mrava87

mrava87 commented Aug 28, 2026

Copy link
Copy Markdown
Contributor Author

@claude one more review

@claude

claude Bot commented Aug 28, 2026

Copy link
Copy Markdown

Claude finished @mrava87's task in 2m 57s —— View job


Todo List

  • Gather context (PR diff, previous review findings, latest commit f954050)
  • Re-check finding from last review (ProximalGradient/AndersonProximalGradient tau scalar handling)
  • Review full current diff for any remaining issues
  • Post updated review

Summary

Commit f954050 ("fixed .item() issue also for proximal gradient solvers") correctly resolves the finding from the previous review: ProximalGradient.step and AndersonProximalGradient.step now extract a Python scalar tau up front (pyproximal/optimization/cls_primal.py:625-628, :1086-1089) and use it consistently, including through _backtracking (which now receives/returns a real float since tau is a scalar going in). self.tau is correctly re-stored as a size-1 array after backtracking (cls_primal.py:660). Good follow-up.

New finding

GeneralizedProximalGradient.step calls .item() on self.weights[i], which breaks when weights is passed as a plain Python list/tuple

pyproximal/optimization/cls_primal.py:1539 and :1542:

ztmp = proxg.prox(
    ztmp, self.tau * self.epsg[i].item() / self.weights[i].item()
)
...
x += self.weights[i].item() * self.zs[i]

self.weights is set at cls_primal.py:1445-1447 as np.ones(len(proxgs)) / len(proxgs) if weights is None else weights — i.e. when the user supplies their own weights, whatever type they pass through (list, tuple, ndarray) is stored as-is, with no np.asarray/self.ncp.asarray conversion. If a user passes weights=[0.5, 0.5] (a plain Python list — this is exactly the kind of value used in pytests/test_solver.py:149 for test_GPG_weights, just with a sum that doesn't pass the ==1.0 check so it never reaches .step()), then self.weights[i] is a native Python float, and float has no .item() method, so this now raises:

AttributeError: 'float' object has no attribute 'item'

This is a regression that would surface on plain NumPy/CPU as well as GPU, whenever weights is user-supplied as a list/tuple with valid values (e.g. weights=[0.5, 0.5] summing to 1, matching the number of proxgs) — not just the CuPy 0-d-array scenario this PR targets. It isn't currently caught by the test suite because the only test exercising list-typed weights is the error-path test at pytests/test_solver.org:141 that fails validation before .step() is ever called.

Suggested fix: normalize self.weights to an ncp/numpy array once at setup time (similar to how self.epsg is handled at cls_primal.py:1453), e.g. self.weights = self.ncp.asarray(self.weights, dtype=np.float32), rather than calling .item() defensively at each use inside step().

Fix this →

Nitpicks

  • pyproximal/optimization/cls_primal.py:660: self.tau = self.ncp.atleast_1d(self.ncp.asarray(tau, dtype=np.float32)) after backtracking hardcodes dtype=np.float32, but this matches the existing convention used in setup (cls_primal.py:563), so not a new issue.
  • The rest of the diff (cls_primaldual.py, L2.py's math.sqrt fix) is unchanged from the previous review and still looks correct.

Note: I don't have a Python execution environment available in this sandbox (all python3 invocations required approval that wasn't grantable in this non-interactive run), so I could not actually execute a reproduction script — this finding is based on static tracing of self.weights' type through setup()/step(), cross-referenced against the list literal already used for weights in pytests/test_solver.py:149.

@mrava87

mrava87 commented Aug 28, 2026

Copy link
Copy Markdown
Contributor Author

@claude review again

@claude

claude Bot commented Aug 28, 2026

Copy link
Copy Markdown

Claude finished @mrava87's task in 1m 48s —— View job


Todo List

  • Gather context (PR diff, previous review findings, latest commits ab2e7af, f954050, d634983)
  • Re-check finding from last review (GeneralizedProximalGradient.weights .item() on plain list/tuple)
  • Review full current diff for any remaining issues
  • Post updated review

Summary

All previously-reported tau/epsg/mu scalar-extraction findings are now resolved:

  • ProximalGradient.step and AndersonProximalGradient.step extract a scalar tau up front (pyproximal/optimization/cls_primal.py:625-628, :1086-1089) and use it consistently throughout, including through _backtracking (:651-660), which now genuinely receives/returns a Python float (the old cast(float, self.tau) no-op is gone).
  • HQS.step (:1927-1930) and PrimalDual.step (pyproximal/optimization/cls_primaldual.py:297-303) call .item() consistently in both the constant and per-iteration branches.
  • GeneralizedProximalGradient.step calls .item() on self.epsg[i] (cls_primal.py:1539).
  • L2.prox uses math.sqrt instead of np.sqrt (pyproximal/proximal/L2.py:271-272), avoiding the 0-d array re-promotion that triggered the original CuPy TypeError.

Remaining finding

GeneralizedProximalGradient: the weights[i].item() fix (cls_primal.py:1539, 1542) is still only papered over by documentation/tests, not by the code

The latest commit (d634983) only changes the docstrings (weights : :obj:float → `:obj:`numpy.ndarray, in both cls_primal.py:1411 and pyproximal/optimization/primal.py:464) and updates the tests (pytests/test_solver.py:150, :704) to pass np.array(...) instead of a plain list/tuple. This avoids tripping the new AttributeError: 'float' object has no attribute 'item' bug in the test suite, but the underlying code at cls_primal.py:1445-1450 still stores self.weights as-is with no np.asarray/self.ncp.asarray conversion:

self.weights = (
    np.ones(len(proxgs)) / len(proxgs) if weights is None else weights
)
if len(self.weights) != len(self.proxgs) or np.sum(self.weights) != 1.0:
    ...

Any caller who follows the (pre-existing) type hint weights: NDArray | None = None loosely — e.g. passing a plain Python list/tuple of floats, which is idiomatic and duck-type-compatible with numpy.ndarray everywhere else weights is validated (len(), np.sum()) — will still hit AttributeError: 'float' object has no attribute 'item' at cls_primal.py:1539/1542, since a list/tuple element is a plain Python float with no .item() method. This is a regression introduced by this PR (previously self.weights[i] was used directly with no .item() call, so lists worked fine), and it now affects plain NumPy/CPU usage too, not just the CuPy scenario this PR targets.

Suggested fix: normalize once at setup time, matching how epsg is already handled two lines below (cls_primal.py:1453):

self.weights = self.ncp.asarray(
    np.ones(len(proxgs)) / len(proxgs) if weights is None else weights,
    dtype=np.float32,
)

This removes the need for defensive .item() calls in step() and makes list/tuple/ndarray input all work uniformly, on both CPU and GPU.

Fix this →

Nitpicks

  • The rest of the diff (cls_primaldual.py, primal.py docstring, L2.py's math.sqrt fix) is correct and unchanged in substance from the previous review.

I was not able to run the test suite or a GPU/CuPy environment in this sandbox (no interactive Python execution available), so the finding above is based on static tracing of weights through setup()/step(), matching the exact scenario flagged in the previous review round.

@mrava87
mrava87 merged commit f552a8a into PyLops:dev Aug 29, 2026
16 checks passed
@mrava87
mrava87 deleted the fix_l2_scalar_times_op branch August 29, 2026 12:57
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant