From fb971e4e74aeef87655b3c7082449114ef933393 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Thu, 2 Jul 2026 22:29:32 +0100 Subject: [PATCH 1/5] feat: added sparse option to convmtx utility routines --- pylops/avo/poststack.py | 4 ++- pylops/basicoperators/matrixmult.py | 2 +- pylops/utils/backend.py | 52 ++++++++++++++++++++++++++++- pylops/utils/signalprocessing.py | 38 ++++++++++++++++++--- pytests/test_signalutils.py | 49 ++++++++++++++++----------- 5 files changed, 118 insertions(+), 27 deletions(-) diff --git a/pylops/avo/poststack.py b/pylops/avo/poststack.py index 27394ba59..a33390114 100644 --- a/pylops/avo/poststack.py +++ b/pylops/avo/poststack.py @@ -122,7 +122,9 @@ def _PoststackLinearModelling( ) else: Cop = _MatrixMult( - nonstationary_convmtx(wav, nt0, hc=wav.shape[1] // 2, pad=(nt0, nt0)), + nonstationary_convmtx( + wav, nt0, hc=wav.shape[1] // 2, pad=(nt0, nt0), sparse=sparse + ), otherdims=spatdims, dtype=dtype, **args_MatrixMult, diff --git a/pylops/basicoperators/matrixmult.py b/pylops/basicoperators/matrixmult.py index 8b6bd805a..d92c49287 100644 --- a/pylops/basicoperators/matrixmult.py +++ b/pylops/basicoperators/matrixmult.py @@ -23,7 +23,7 @@ class MatrixMult(LinearOperator): Parameters ---------- - A : :obj:`numpy.ndarray` or :obj:`scipy.sparse` matrix + A : :obj:`numpy.ndarray` or :obj:`scipy.sparse.spmatrix` Matrix. otherdims : :obj:`tuple`, optional Number of samples for each other dimension of model diff --git a/pylops/utils/backend.py b/pylops/utils/backend.py index 91493ff38..eb005332c 100644 --- a/pylops/utils/backend.py +++ b/pylops/utils/backend.py @@ -12,6 +12,8 @@ "get_block_diag", "get_toeplitz", "get_csc_matrix", + "get_csr_matrix", + "get_dia_matrix", "get_sparse_eye", "get_lstsq", "get_sp_fft", @@ -35,7 +37,7 @@ import scipy.fft as sp_fft from scipy.linalg import block_diag, lstsq, toeplitz from scipy.signal import convolve, correlate, fftconvolve, oaconvolve -from scipy.sparse import csc_matrix, eye +from scipy.sparse import csc_matrix, csr_matrix, dia_matrix, eye from pylops.utils import deps from pylops.utils.typing import ArrayLike, DTypeLike, NDArray, Tfftengine_ncj @@ -51,6 +53,8 @@ from cupyx.scipy.signal import fftconvolve as cp_fftconvolve from cupyx.scipy.signal import oaconvolve as cp_oaconvolve from cupyx.scipy.sparse import csc_matrix as cp_csc_matrix + from cupyx.scipy.sparse import csr_matrix as cp_csr_matrix + from cupyx.scipy.sparse import dia_matrix as cp_dia_matrix from cupyx.scipy.sparse import eye as cp_eye if deps.jax_enabled: @@ -382,6 +386,52 @@ def get_csc_matrix(x: ArrayLike) -> Callable: return cp_csc_matrix +def get_csr_matrix(x: ArrayLike) -> Callable: + """Returns correct csr_matrix module based on input + + Parameters + ---------- + x : :obj:`numpy.ndarray` + Array + + Returns + ------- + f : :obj:`callable` + Function to be used to process array + + """ + if not deps.cupy_enabled: + return csr_matrix + + if cp.get_array_module(x) == np: + return csr_matrix + else: + return cp_csr_matrix + + +def get_dia_matrix(x: ArrayLike) -> Callable: + """Returns correct dia_matrix module based on input + + Parameters + ---------- + x : :obj:`numpy.ndarray` + Array + + Returns + ------- + f : :obj:`callable` + Function to be used to process array + + """ + if not deps.cupy_enabled: + return dia_matrix + + if cp.get_array_module(x) == np: + return dia_matrix + else: + return cp_dia_matrix + + def get_sparse_eye(x: ArrayLike) -> Callable: """Returns correct sparse eye module based on input diff --git a/pylops/utils/signalprocessing.py b/pylops/utils/signalprocessing.py index 2c615a1c0..ab721c0ea 100644 --- a/pylops/utils/signalprocessing.py +++ b/pylops/utils/signalprocessing.py @@ -18,13 +18,15 @@ from pylops.utils._pwd2d import _conv_allpass, _triangular_smoothing_from_boxcars from pylops.utils.backend import ( get_array_module, + get_csr_matrix, + get_dia_matrix, get_normalize_axis_index, get_toeplitz, ) from pylops.utils.typing import NDArray, Tpwdsmoothing -def convmtx(h: NDArray, n: int, offset: int = 0) -> NDArray: +def convmtx(h: NDArray, n: int, offset: int = 0, sparse: bool = False) -> NDArray: r"""Convolution matrix Makes a dense convolution matrix :math:`\mathbf{C}` @@ -42,12 +44,16 @@ def convmtx(h: NDArray, n: int, offset: int = 0) -> NDArray: Convolution filter (1D array) n : :obj:`int` Number of columns of convolution matrix - offset : :obj:`int` + offset : :obj:`int`, optional Index of the center of the filter + sparse : :obj:`bool`, optional + .. versionadded:: 2.9.0 + + Return dense (``False``) or sparse (``True``) matrix Returns ------- - C : :obj:`numpy.ndarray` + C : :obj:`numpy.ndarray` or :obj:`scipy.sparse.spmatrix` Convolution matrix of size :math:`\text{len}(h)+n-1 \times n` """ @@ -62,12 +68,22 @@ def convmtx(h: NDArray, n: int, offset: int = 0) -> NDArray: ) ncp = get_array_module(h) + + # create Toeplitz matrix nh = len(h) col_1 = ncp.r_[h, ncp.zeros(n + nh - 2, dtype=h.dtype)] row_1 = ncp.r_[h[0], ncp.zeros(n - 1, dtype=h.dtype)] C = get_toeplitz(h)(col_1, row_1) + # apply offset C = C[offset : offset + nh + n - 1] + + # convert to sparse using the following rule-of-thumb: + # - DIA format very short filters (<= 11) + # - CSR format for other filters + if sparse: + C = get_dia_matrix(h)(C) if nh <= 11 else get_csr_matrix(h)(C) + return C @@ -76,6 +92,7 @@ def nonstationary_convmtx( n: int, hc: int = 0, pad: tuple[int] = (0, 0), + sparse: bool = False, ) -> NDArray: r"""Convolution matrix from a bank of filters @@ -97,18 +114,31 @@ def nonstationary_convmtx( Zero-padding to apply to the bank of filters before and after the provided values (use it to avoid wrap-around or pass filters with enough padding) + sparse : :obj:`bool`, optional + .. versionadded:: 2.9.0 + + Return dense (``False``) or sparse (``True``) matrix Returns ------- - C : :obj:`numpy.ndarray` + C : :obj:`numpy.ndarray` or :obj:`scipy.sparse.spmatrix` Convolution matrix """ ncp = get_array_module(H) + # create Toeplitz matrix + nh = H.shape[1] H = ncp.pad(H, ((0, 0), pad), mode="constant") C = ncp.array([ncp.roll(h, ih) for ih, h in enumerate(H)]) C = C[:, pad[0] + hc : pad[0] + hc + n].T # take away edges + + # convert to sparse using the following rule-of-thumb: + # - DIA format very short filters (<= 11) + # - CSR format for other filters + if sparse: + C = get_dia_matrix(H)(C) if nh <= 11 else get_csr_matrix(H)(C) + return C diff --git a/pytests/test_signalutils.py b/pytests/test_signalutils.py index a4b95299a..d1857f8e5 100644 --- a/pytests/test_signalutils.py +++ b/pytests/test_signalutils.py @@ -1,8 +1,16 @@ import os -import numpy as np +if int(os.environ.get("TEST_CUPY_PYLOPS", 0)): + import cupy as np + from cupy.testing import assert_array_almost_equal + + backend = "cupy" +else: + import numpy as np + from numpy.testing import assert_array_almost_equal + + backend = "numpy" import pytest -from numpy.testing import assert_array_almost_equal from pylops.utils.signalprocessing import convmtx, nonstationary_convmtx, slope_estimate @@ -39,29 +47,25 @@ np.random.seed(10) -@pytest.mark.skipif( - int(os.environ.get("TEST_CUPY_PYLOPS", 0)) == 1, reason="Not CuPy enabled" -) @pytest.mark.parametrize("par", [(par1), (par1j), (par2), (par2j)]) -def test_convmtx(par): +@pytest.mark.parametrize("sparse", [False, True]) +def test_convmtx(par, sparse): """Compare convmtx with np.convolve (small filter)""" x = np.random.normal(0, 1, par["nt"]) + par["imag"] * np.random.normal( 0, 1, par["nt"] ) h = np.hanning(par["nh"]) - H = convmtx(h, par["nt"], par["nh"] // 2) + H = convmtx(h, par["nt"], par["nh"] // 2, sparse=sparse) y = np.convolve(x, h, mode="same") - y1 = np.dot(H[: par["nt"]], x) + y1 = (H @ x)[: par["nt"]] assert_array_almost_equal(y, y1, decimal=4) -@pytest.mark.skipif( - int(os.environ.get("TEST_CUPY_PYLOPS", 0)) == 1, reason="Not CuPy enabled" -) @pytest.mark.parametrize("par", [(par1), (par1j), (par2), (par2j)]) -def test_convmtx1(par): +@pytest.mark.parametrize("sparse", [False, True]) +def test_convmtx1(par, sparse): """Compare convmtx with np.convolve (large filter)""" x = np.random.normal(0, 1, par["nt"]) + par["imag"] * np.random.normal( 0, 1, par["nt"] @@ -69,19 +73,20 @@ def test_convmtx1(par): h = np.hanning(par["nh"]) X = convmtx( - x, par["nh"], par["nh"] // 2 - 1 if par["nh"] % 2 == 0 else par["nh"] // 2 + x, + par["nh"], + par["nh"] // 2 - 1 if par["nh"] % 2 == 0 else par["nh"] // 2, + sparse=sparse, ) y = np.convolve(x, h, mode="same") - y1 = np.dot(X[: par["nt"]], h) + y1 = (X @ h)[: par["nt"]] assert_array_almost_equal(y, y1, decimal=4) -@pytest.mark.skipif( - int(os.environ.get("TEST_CUPY_PYLOPS", 0)) == 1, reason="Not CuPy enabled" -) @pytest.mark.parametrize("par", [(par1), (par1j)]) -def test_nonstationary_convmtx(par): +@pytest.mark.parametrize("sparse", [False, True]) +def test_nonstationary_convmtx(par, sparse): """Compare nonstationary_convmtx with convmtx for stationary filter""" x = np.random.normal(0, 1, par["nt"]) + par["imag"] * np.random.normal( 0, 1, par["nt"] @@ -89,7 +94,10 @@ def test_nonstationary_convmtx(par): h = np.hanning(par["nh"]) H = convmtx( - h, par["nt"], par["nh"] // 2 - 1 if par["nh"] % 2 == 0 else par["nh"] // 2 + h, + par["nt"], + par["nh"] // 2 - 1 if par["nh"] % 2 == 0 else par["nh"] // 2, + sparse=sparse, ) H1 = nonstationary_convmtx( @@ -98,7 +106,8 @@ def test_nonstationary_convmtx(par): hc=par["nh"] // 2, pad=(par["nt"], par["nt"]), ) - y = np.dot(H[: par["nt"]], x) + + y = (H @ x)[: par["nt"]] y1 = np.dot(H1, x) assert_array_almost_equal(y, y1, decimal=4) From 03516178c1b5bd638d9b2b5bda9174e2a8fc6b16 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Thu, 2 Jul 2026 22:47:57 +0100 Subject: [PATCH 2/5] fix: avoid DIA format for CuPy --- pylops/utils/signalprocessing.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/pylops/utils/signalprocessing.py b/pylops/utils/signalprocessing.py index ab721c0ea..6a937554f 100644 --- a/pylops/utils/signalprocessing.py +++ b/pylops/utils/signalprocessing.py @@ -137,8 +137,12 @@ def nonstationary_convmtx( # - DIA format very short filters (<= 11) # - CSR format for other filters if sparse: - C = get_dia_matrix(H)(C) if nh <= 11 else get_csr_matrix(H)(C) - + if ncp == np: + C = get_dia_matrix(H)(C) if nh <= 11 else get_csr_matrix(H)(C) + else: + # For CuPy DIA cannot take a dense matrix, so the dense matrix is + # always converted to CSR format + C = get_csr_matrix(H)(C) return C From 7f2692a74af0f77db17438488d692861e75156f7 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Thu, 2 Jul 2026 22:54:16 +0100 Subject: [PATCH 3/5] test: disable CuPy tests for Bilinear with complex numbers --- pytests/test_interpolation.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/pytests/test_interpolation.py b/pytests/test_interpolation.py index d70638ce1..1da7aed5f 100644 --- a/pytests/test_interpolation.py +++ b/pytests/test_interpolation.py @@ -543,6 +543,9 @@ def test_Bilinear_2dsignal(par: InterpolationTestParameters, dtype: np.dtype): @pytest.mark.parametrize("dtype", [np.float32, np.float64]) def test_Bilinear_2dsignal_flatten(par: InterpolationTestParameters, dtype: np.dtype): """Dot-test and forward for Interp operator for 2d signal with forceflat""" + if par.imag == 1j: + pytest.skip("cupy.add.at currently does not support complex numbers") + np.random.seed(1) dtype1 = (np.empty(0, dtype=dtype) + par.imag * np.empty(0, dtype=dtype)).dtype @@ -572,6 +575,9 @@ def test_Bilinear_2dsignal_flatten(par: InterpolationTestParameters, dtype: np.d @pytest.mark.parametrize("dtype", [np.float32, np.float64]) def test_Bilinear_3dsignal(par: InterpolationTestParameters, dtype: np.dtype): """Dot-test and forward for Interp operator for 3d signal""" + if par.imag == 1j: + pytest.skip("cupy.add.at currently does not support complex numbers") + np.random.seed(1) dtype1 = (np.empty(0, dtype=dtype) + par.imag * np.empty(0, dtype=dtype)).dtype From 12a59761ac3d62b156fe9684c1fe375fd39cd645 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Thu, 2 Jul 2026 22:56:57 +0100 Subject: [PATCH 4/5] fix: restore wrong test changed --- pytests/test_interpolation.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/pytests/test_interpolation.py b/pytests/test_interpolation.py index 1da7aed5f..d70638ce1 100644 --- a/pytests/test_interpolation.py +++ b/pytests/test_interpolation.py @@ -543,9 +543,6 @@ def test_Bilinear_2dsignal(par: InterpolationTestParameters, dtype: np.dtype): @pytest.mark.parametrize("dtype", [np.float32, np.float64]) def test_Bilinear_2dsignal_flatten(par: InterpolationTestParameters, dtype: np.dtype): """Dot-test and forward for Interp operator for 2d signal with forceflat""" - if par.imag == 1j: - pytest.skip("cupy.add.at currently does not support complex numbers") - np.random.seed(1) dtype1 = (np.empty(0, dtype=dtype) + par.imag * np.empty(0, dtype=dtype)).dtype @@ -575,9 +572,6 @@ def test_Bilinear_2dsignal_flatten(par: InterpolationTestParameters, dtype: np.d @pytest.mark.parametrize("dtype", [np.float32, np.float64]) def test_Bilinear_3dsignal(par: InterpolationTestParameters, dtype: np.dtype): """Dot-test and forward for Interp operator for 3d signal""" - if par.imag == 1j: - pytest.skip("cupy.add.at currently does not support complex numbers") - np.random.seed(1) dtype1 = (np.empty(0, dtype=dtype) + par.imag * np.empty(0, dtype=dtype)).dtype From 7702af582e84b9edc0be856135a2f3ce761c89b0 Mon Sep 17 00:00:00 2001 From: Matteo Ravasi Date: Thu, 2 Jul 2026 23:20:41 +0100 Subject: [PATCH 5/5] fix: avoid DIA in convmtx for Cupy arrays Refactor sparse matrix conversion logic for clarity. --- pylops/utils/signalprocessing.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/pylops/utils/signalprocessing.py b/pylops/utils/signalprocessing.py index 6a937554f..f5d70bb7e 100644 --- a/pylops/utils/signalprocessing.py +++ b/pylops/utils/signalprocessing.py @@ -82,8 +82,12 @@ def convmtx(h: NDArray, n: int, offset: int = 0, sparse: bool = False) -> NDArra # - DIA format very short filters (<= 11) # - CSR format for other filters if sparse: - C = get_dia_matrix(h)(C) if nh <= 11 else get_csr_matrix(h)(C) - + if ncp == np: + C = get_dia_matrix(h)(C) if nh <= 11 else get_csr_matrix(h)(C) + else: + # For CuPy DIA cannot take a dense matrix, so the dense matrix is + # always converted to CSR format + C = get_csr_matrix(h)(C) return C