Preserve convolution kernel axis order at every filter rank

[antonsoo]

With unordered axes, filters with fewer dimensions than the model are reshaped without reordering their values. When filter and model ranks match, axes is ignored. For example, a kernel applied with axes=(2, 0) can represent a different spatial filter. Forward and adjoint dot-tests alone miss this because both use the same wrong kernel.

Transpose the padded filter into physical axis order before the existing singleton-dimension embedding, at every filter rank. ConvolveND and its SmoothingND wrapper now infer omitted axes from the filter rank, so default 3D-filter and 3D-smoothing calls still work when the transpose is applied consistently. Explicit axes correspond to the filter dimensions in order. Invalid axis counts, repeated normalized axes and filters with more dimensions than the model now raise clear ValueErrors before padding or transposition. Convolve2D keeps its existing default pair of axes.

This example differs by 565 on dev and agrees to floating-point precision with the correction:

import numpy as np
from scipy.signal import convolve2d
from pylops.signalprocessing import Convolve2D

x = np.arange(5 * 3 * 7, dtype=float).reshape(5, 3, 7)
h = np.arange(1, 7, dtype=float).reshape(2, 3)
op = Convolve2D(x.shape, h=h, offset=(0, 1), axes=(2, 0))
expected = np.stack([
    convolve2d(x[:, i, :], h.T, mode="full")[1:6, :7]
    for i in range(3)
], axis=1)
print(np.max(np.abs(op @ x - expected)))

Validation on Linux / RTX 5070:

• The convolution, smoothing and PWD area passes 462 tests on CPU and 454 on CuPy, with eight unsupported PWD cases skipped. The convolution modelling cases use the existing parameter dictionaries, seeded inputs, cyclic/negative axes, dtype checks and explicit dot-test sizes. The SmoothingND cases are at the end of the file, use par1/par2, and check analytic impulse responses for 2D, 3D and 4D models with omitted, None and permuted axes.
• Twelve invalid-input regression cases fail before this follow-up and pass afterward. They cover missing/extra axes, repeated negative aliases and filters with more dimensions than the model, including singleton filters.
• All 288 complete photo workload records remain identical to the earlier independently verified results on CPU and CUDA. References use sparse convolution matrices assembled from physical coordinates, direct sparse solves and separable smoothing averages, checking forward, adjoint and inverse results.
• Project make lint passes. The preceding Sphinx API build passed with gallery execution disabled and one configuration cache warning; this test-only follow-up leaves every library and documentation file unchanged.

AI disclosure: OpenAI Codex was used for investigation, implementation, tests, independent CPU/CUDA validation, documentation and PR prose.

[codacy-production]

Up to standards ✅

🟢 Issues 0 issues

Results:
0 new issues

View in Codacy

🟢 Metrics 0 complexity · 0 duplication

Metric	Results
Complexity	0
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]

@antonsoo thanks for this contribution. It does indeed seem to be a bug!

I will have a closer look in coming days and let you know if I have any further comments (already left a few minor ones)

[antonsoo]

Thanks for checking. I will address the minor points once the review is published; there are no inline comments visible yet.

[mrava87]

@antonsoo thanks for spotting this inconsistency... we should be able to support unordered axes inputs and I think this is the correct way of doing it... I requested a few changes aimed at keeping consistency with the rest of the codebase before I can merge.

[mrava87]

@claude review. Check if any similar pattern should be fixed in other operators using axes

[claude]

Claude finished @mrava87's task in 30s —— View job

---

Review of #818: filter axis order in ConvolveND/Convolve2D

I read the diff and swept the other axes-taking operators. I didn't run the test suite.

The fix is correct

• In pylops/signalprocessing/convolvend.py:128, dimsh[ax] = nh[iax] puts filter dim i at model axis axes[i]. Transposing by argsort(axes) first makes the data match that layout. After the transpose, the filter's dims are in increasing model-axis order, which is what the reshape assumes. Ordered axes give the identity permutation, and negative axes are already normalised into self.axes.
• offset and the padding are applied per filter dimension before the transpose, so offsets stay tied to h's own dims. That is consistent.
• The new test is independent of the implementation. It builds the expected result from shifted slices and checks both forward and adjoint. That matters because a dot-test can't catch this bug.

Open point: the full-rank case (len(dims) == h.ndim)

This is the point in mrava87's inline comment. With a full-rank h the code skips the reshape and never uses axes. So ConvolveND((a, b, c), h3d, axes=(2, 0, 1)) still applies h dim i to model axis i, which contradicts the new docstring wording, "filter dimension i is applied along model axis axes[i]".

For consistency, I suggest one of these:

1. Always transpose, which is what mrava87 hinted at. Drop the if len(self.dims) != len(self.nh) condition around the transpose. Keep the reshape only when the ranks differ. Then h.transpose(argsort(axes)) runs in all cases. With the default axes=(-2, -1) and h.ndim == len(dims) == 2, the permutation is the identity, so default usage is unchanged. For 3D dims with a 3D h and the default axes, the constructor would fail or behave oddly anyway (the dims/axes mismatch, which is not covered here).
2. Keep the full-rank path unchanged, and document that axes is only used when h has fewer dims than the model.

Either way, add a full-rank test with unsorted axes, e.g. dims=(4,5,6), a 3D h, axes=(2,0,1), and a loop-based expected result like the one in the new test.

Other operators using axes

• SmoothingND and Smoothing2D (pylops/basicoperators/smoothingnd.py:73, smoothing2d.py:74): they build h = ones(nsmooth) and delegate to ConvolveND/Convolve2D, so they get the fix for free. Before it, they were wrong for unsorted axes with different nsmooth per axis. This is worth a small test, e.g. nsmooth=(3, 5) with axes=(2, 0). Smoothing with equal sizes along every axis was never affected.
• FFT family (fft2d.py, fftnd.py): axes goes straight to numpy/scipy fft(s=nffts, axes=axes), so nffts[i] follows axes[i] and unsorted axes are handled natively. No change needed.
• Laplacian (laplacian.py:92): axes, weights and sampling are zipped together, so the order is consistent. No change needed.
• Transpose: the permutation is the point of the operator, so nothing to fix.
• Convolve1D: it takes a single axis, so the issue doesn't arise.
• Not checked line by line: Gradient, Sliding*, Patch*, DWT*/SWT*, DCT, Radon* and the rest. Wherever per-axis parameters (e.g. nwin, level) are paired with axes, it may be worth a quick look for the same reshape-without-transpose pattern. A search for .reshape( right after building a per-axis shape would find it.

Remaining requests from mrava87's review (not yet addressed)

• Remove the CHANGELOG entry.
• Remove the assertion block from examples/plot_convolve.py, since examples shouldn't contain testing commands.
• Revert the docstring additions in convolve2d.py and convolvend.py.
• Rework the test to follow the test_Convolve3D style: np.random.seed() and a "3D on 4D (only modelling)" pattern with shifted axes in place of [0, 1, 2].

I haven't made any changes to the branch. I can apply these if you ask.

[antonsoo]

Updated in 198bfdf: the transpose now applies at every filter rank. ConvolveND infers omitted axes from h.ndim, since retaining the old two-axis default would break the existing 3D-filter calls once the transpose is unconditional. Convolve2D keeps its existing default axes.

The tests now follow the seeded 3D/4D modelling convention with shifted and negative axes. The example plots its results, and the old release entry and extra axis-mapping prose are removed. CPU and CuPy convolution suites pass; independent photograph matrix/inverse checks cover both embedded and full-rank filters, with unchanged controls. CuPy's unsupported multidimensional automatic-method selection is unchanged and documented in the validation notes.

[antonsoo]

Read the Docs exposed a wrapper default I missed: SmoothingND still passed the old two-axis default for a 3D filter. Fixed in 4d81ea4 by letting the wrapper infer omitted axes too.

Both failing examples now finish locally, and the targeted Sphinx gallery executes plot_convolve, plot_smoothing3d and plot_pwdslope successfully. CPU convolution/smoothing/PWD checks pass; CuPy convolution/smoothing checks pass with the unsupported PWD cases skipped. Added default-axis regressions cover 2D, 3D and 4D models. Independent photo averages also match on CPU and CUDA, with the existing controls unchanged. The description is updated; the hosted rebuild still needs to finish.

[mrava87]

A few more comments

[mrava87]

@claude review again

[claude]

Claude finished @mrava87's task in 24s —— View job

---

Review of #818 (second pass)

I read the diff against origin/dev and the final convolvend.py. I didn't run the tests.

Implementation: correct

• convolvend.py:123 transposes by argsort(axes) for every filter rank, after offset padding. The padding is applied per h dimension in h's own order, so offsets stay tied to the right dimensions.
• self.nh is captured before the transpose and dimsh[ax] = nh[iax] uses the original filter order. After the transpose, filter dimension j sits at the j-th smallest axis, so the singleton-embedding reshape is consistent. Full-rank filters skip the reshape and now get the permutation too.
• Defaulting axes=None to the last h.ndim axes is a sensible way to keep 3D-filter calls working. Convolve2D keeps (-2, -1). The SmoothingND default change works the same way. Smoothing2D is unaffected.
• The FFT family, Laplacian and Convolve1D don't need the same fix. Both of those checks come from my earlier pass.

Minor suggestions

• There is no validation that len(axes) == h.ndim or that the axes are unique. With a mismatch, argsort(axes).transpose raises a cryptic error. A short ValueError would be friendlier, but this isn't blocking.
• SmoothingND's docstring says axes is tuple or list, while ConvolveND's still says int. This is pre-existing, so I'd leave it.
• zip(..., strict=True) in test_smoothing.py needs Python ≥ 3.10. Please check this against requires-python in pyproject.toml.

mrava87's latest requests (19:41) still open

The latest commit, 4d81ea4, doesn't address these. They are all in the tests.

pytests/test_convolve.py (test_Convolve3D_axes):

• Use the par1_3d / par2_3d parameter dicts, as test_Convolve3D does, and parametrize only axes and dtype on top. The ndims and dims logic then goes into the dicts.
• Add explanatory comments like those in test_Convolve3D, for example "3D on 3D" and "3D on 4D (only modelling)".
• Rename reference to Cop_ref.
• Pass explicit dimensions to dottest (Cop.shape[0], Cop.shape[1], as already done), with dims/dimsd kept explicit in the same way as the other tests. Please check this matches the convention in test_Convolve3D.

pytests/test_smoothing.py (test_SmoothingND_axes):

• Pass explicit dimensions to dottest, i.e. dottest(op, np.prod(dims), np.prod(dims), ...), instead of dottest(op, rtol=1e-4, backend=backend).
• Drop the .ravel() / .reshape(dims). op @ x and op.H @ x work directly on N-D arrays.

Summary

The core fix is sound and I have no correctness concerns. The review comments that remain are test-style conventions. I made no changes to the branch.

[mrava87]

@antonsoo handle also these two comments:

• There is no validation that len(axes) == h.ndim or that the axes are unique. With a mismatch, argsort(axes).transpose raises a cryptic error. A short ValueError would be friendlier, but this isn't blocking.
• SmoothingND's docstring says axes is tuple or list, while ConvolveND's still says int. This is pre-existing, so I'd leave it.

[antonsoo]

The remaining macOS/Python 3.13t failure is in four existing float32 prestack inversions (LinAlgError from SciPy's SVD). The matching 0fed249 base job passed. Those cases pass on base and head under Linux/Python 3.13.16t with the same NumPy 2.4.3/SciPy 1.17.1, and the least-squares matrix/RHS hashes are identical, including after changing the ambient RNG seed. The prestack code, test and solver path are unchanged. This does not reproduce the macOS environment; could that job be rerun?

[mrava87]

@claude review again and check if comments in previous review are now handled.

[claude]

Claude finished @mrava87's task in 28s —— View job

---

Review of #818 (third pass)

I diffed 0fed249..HEAD and read the final convolvend.py and smoothingnd.py. I didn't run the tests.

Previous comments

Request	Status
Remove CHANGELOG entry	✅ Done. The diff has no CHANGELOG change.
No assertions in examples/plot_convolve.py	✅ Done. The example now plots only.
Revert docstring additions in convolve2d.py and convolvend.py	⚠️ Mostly done. convolve2d.py is untouched. convolvend.py has a one-line addition ("If None, use the last h.ndim axes"), which documents the new default and is reasonable to keep.
Always transpose, at every filter rank	✅ Done (convolvend.py:123).
test_Convolve3D_axes uses par1_3d / par2_3d, with axes and dtype parametrized	✅ Done.
"3D on 3D" and "3D on 4D (only modelling)" comments	✅ Done.
reference renamed to Cop_ref	✅ Done.
Explicit sizes in dottest	✅ Done in both test files.
test_SmoothingND_axes: no .ravel() / .reshape()	✅ Done.
zip(strict=True) and Python version	✅ Fine. requires-python is >=3.11.
Validate len(axes) == h.ndim and that axes are unique	❌ Not addressed. mrava87 asked for this explicitly.
SmoothingND docstring says tuple or list, while ConvolveND says int	⚠️ SmoothingND was updated. ConvolveND still says :obj:int``. I'd fix it, since the same file is already being touched.

Remaining items

1. Add axes validation in ConvolveND.__init__, right after self.axes is built (convolvend.py:90). I'd raise ValueError when len(self.axes) != h.ndim, and another when len(set(self.axes)) != len(self.axes). Without this, a mismatch fails inside transpose with a cryptic message.
2. The new default can fail silently. With axes=None, np.arange(len(self.dims) - h.ndim, len(self.dims)) starts at a negative number if h.ndim > len(dims). The result is wrong axes rather than an error. The same check in item 1 won't catch it, because the lengths still match. Check h.ndim <= len(dims) first. The existing "dims and h of different dimensions" check from an earlier commit may already cover this. Please confirm, and add a pytest.raises test at the top of the test file if it doesn't.
3. Add tests for the new errors. Per CLAUDE.md, pytest.raises tests go above the normal tests.
4. ConvolveND docstring: change axes : :obj:int to `:obj:`tuple` or :obj:`list.

Correctness

No concerns with the core logic:

• self.nh is captured before the transpose.
• dimsh[ax] = nh[iax] uses the original filter order.
• The argsort transpose puts the filter dimensions in increasing model-axis order, which the singleton reshape assumes.
• The SmoothingND default (axes=None) goes through the same inference.

I made no changes to the branch.

[mrava87]

@antonsoo I left one last inline comment, didn't realize the new test in test_smoothing was also not aligned with our test style... plus, look at Remaining items in Claude last comment.

[codacy-production]

Up to standards ✅

🟢 Issues 0 issues

Results:
0 new issues

View in Codacy

🟢 Metrics 0 complexity · 0 duplication

Metric	Results
Complexity	0
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.}

[antonsoo]

Added checks for the number of axes and for duplicates after normalizing negative axes. Filters with more dimensions than the model now fail early, including filters whose dimensions are all 1. The 12 new error cases follow the updated test guide. I left the older axes type wording unchanged as suggested.

CPU: 460 pass. CuPy: 452 pass, eight unsupported PWD skips. All 288 photo records are identical to the earlier independent references. Lint and the API docs build pass; the docs build skips gallery execution.