Skip to content

Feat: change internal handling of Zero in hstack/vstack #702

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mrava87 merged 2 commits into from
Oct 14, 2025
Merged

Feat: change internal handling of Zero in hstack/vstack #702

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
2609b33
feat: change internal handling of Zero in hstack/vstack
mrava87 Oct 12, 2025
4db0967
minor: added more info about using Zero in stacking operators
mrava87 Oct 13, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
200 changes: 135 additions & 65 deletions examples/plot_stacking.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,9 +13,10 @@
types of optimization problems with regularization or preconditioning.

Some of this operators naturally lend to embarassingly parallel computations.
Within PyLops we leverage the multiprocessing module to run multiple processes
at the same time evaluating a subset of the operators involved in one of the
stacking operations.
Within PyLops we leverage the
`Multiprocessing <https://docs.python.org/3/library/multiprocessing.html>`_
module to run multiple processes at the same time evaluating a subset of the
operators involved in one of the stacking operations.
"""
import matplotlib.pyplot as plt
import numpy as np
Expand All @@ -35,8 +36,9 @@
# can be performed between two or more linear operators.
# This can be easily achieved using the ``*`` symbol
#
# .. math::
# \mathbf{D_{cat}}= \mathbf{D_v} \mathbf{D_h}
# .. math::
# \mathbf{D_{cat}}= \mathbf{D_v} \mathbf{D_h}

Nv, Nh = 11, 21
X = np.zeros((Nv, Nh))
X[int(Nv / 2), int(Nh / 2)] = 1
Expand All @@ -60,17 +62,18 @@
###############################################################################
# We now want to *vertically stack* three operators
#
# .. math::
# \mathbf{D_{Vstack}} =
# \begin{bmatrix}
# \mathbf{D_v} \\
# \mathbf{D_h}
# \end{bmatrix}, \qquad
# \mathbf{y} =
# \begin{bmatrix}
# \mathbf{D_v}\mathbf{x} \\
# \mathbf{D_h}\mathbf{x}
# \end{bmatrix}
# .. math::
# \mathbf{D_{Vstack}} =
# \begin{bmatrix}
# \mathbf{D_v} \\
# \mathbf{D_h}
# \end{bmatrix}, \qquad
# \mathbf{y} =
# \begin{bmatrix}
# \mathbf{D_v}\mathbf{x} \\
# \mathbf{D_h}\mathbf{x}
# \end{bmatrix}

Nv, Nh = 11, 21
X = np.zeros((Nv, Nh))
X[int(Nv / 2), int(Nh / 2)] = 1
Expand All @@ -91,17 +94,54 @@
plt.tight_layout()
plt.subplots_adjust(top=0.8)

###############################################################################
# If one now wants to have one or more null operators in the stack, it is
# reccomended to use the :py:class:`pylops.Zero` operator (instead of, for example,
# a :py:class:`pylops.MatrixMult` filled with a zero matrix); under the hood, a
# :py:class:`pylops.Zero` operator will be simply by-passed both in the forward and
# adjoint steps.
#
# .. math::
# \mathbf{D_{Vstack}} =
# \begin{bmatrix}
# \mathbf{D_v} \\
# \mathbf{0} \\
# \mathbf{D_h}
# \end{bmatrix}, \qquad
# \mathbf{y} =
# \begin{bmatrix}
# \mathbf{D_v}\mathbf{x} \\
# \mathbf{0} \\
# \mathbf{D_h}\mathbf{x}
# \end{bmatrix}
#
# Note that this feature will become particularly handy when defining
# :py:class:`pylops.Block` operators with large zero blocks.

# With MatrixMult operator filled with zeros
Zop = pylops.MatrixMult(np.zeros((Nv * Nh, Nv * Nh)))
Dstack = pylops.VStack([D2vop, Zop, D2hop])
Y = np.reshape(Dstack * X.ravel(), (Nv * 3, Nh))

# With Zero operator
Zop = pylops.Zero(Nv * Nh)
D1stack = pylops.VStack([D2vop, Zop, D2hop])
Y1 = np.reshape(D1stack * X.ravel(), (Nv * 3, Nh))

print("Y == Y1:", np.allclose(Y, Y1))

###############################################################################
# Similarly we can now *horizontally stack* three operators
#
# .. math::
# \mathbf{D_{Hstack}} =
# \begin{bmatrix}
# \mathbf{D_v} & 0.5*\mathbf{D_v} & -1*\mathbf{D_h}
# \end{bmatrix}, \qquad
# \mathbf{y} =
# \mathbf{D_v}\mathbf{x}_1 + 0.5*\mathbf{D_v}\mathbf{x}_2 -
# \mathbf{D_h}\mathbf{x}_3
# .. math::
# \mathbf{D_{Hstack}} =
# \begin{bmatrix}
# \mathbf{D_v} & 0.5*\mathbf{D_v} & -1*\mathbf{D_h}
# \end{bmatrix}, \qquad
# \mathbf{y} =
# \mathbf{D_v}\mathbf{x}_1 + 0.5*\mathbf{D_v}\mathbf{x}_2 -
# \mathbf{D_h}\mathbf{x}_3

Nv, Nh = 11, 21
X = np.zeros((Nv * 3, Nh))
X[int(Nv / 2), int(Nh / 2)] = 1
Expand All @@ -128,19 +168,20 @@
# We can even stack them both *horizontally* and *vertically* such that we
# create a *block* operator
#
# .. math::
# \mathbf{D_{Block}} =
# \begin{bmatrix}
# \mathbf{D_v} & 0.5*\mathbf{D_v} & -1*\mathbf{D_h} \\
# \mathbf{D_h} & 2*\mathbf{D_h} & \mathbf{D_v} \\
# \end{bmatrix}, \qquad
# \mathbf{y} =
# \begin{bmatrix}
# \mathbf{D_v} \mathbf{x_1} + 0.5*\mathbf{D_v} \mathbf{x_2} -
# \mathbf{D_h} \mathbf{x_3} \\
# \mathbf{D_h} \mathbf{x_1} + 2*\mathbf{D_h} \mathbf{x_2} +
# \mathbf{D_v} \mathbf{x_3}
# \end{bmatrix}
# .. math::
# \mathbf{D_{Block}} =
# \begin{bmatrix}
# \mathbf{D_v} & 0.5*\mathbf{D_v} & -1*\mathbf{D_h} \\
# \mathbf{D_h} & 2*\mathbf{D_h} & \mathbf{D_v} \\
# \end{bmatrix}, \qquad
# \mathbf{y} =
# \begin{bmatrix}
# \mathbf{D_v} \mathbf{x_1} + 0.5*\mathbf{D_v} \mathbf{x_2} -
# \mathbf{D_h} \mathbf{x_3} \\
# \mathbf{D_h} \mathbf{x_1} + 2*\mathbf{D_h} \mathbf{x_2} +
# \mathbf{D_v} \mathbf{x_3}
# \end{bmatrix}

Bop = pylops.Block([[D2vop, 0.5 * D2vop, -1 * D2hop], [D2hop, 2 * D2hop, D2vop]])
Y = np.reshape(Bop * X.ravel(), (2 * Nv, Nh))

Expand All @@ -157,23 +198,51 @@
plt.tight_layout()
plt.subplots_adjust(top=0.8)

###############################################################################
# And here with some Zero blocks
#
# .. math::
# \mathbf{D_{Block}} =
# \begin{bmatrix}
# \mathbf{D_v} & \mathbf{0} & -1*\mathbf{D_h} \\
# \mathbf{D_h} & 2*\mathbf{D_h} & \mathbf{0} \\
# \end{bmatrix}, \qquad
# \mathbf{y} =
# \begin{bmatrix}
# \mathbf{D_v} \mathbf{x_1} - \mathbf{D_h} \mathbf{x_3} \\
# \mathbf{D_h} \mathbf{x_1} + 2*\mathbf{D_h} \mathbf{x_2}
# \end{bmatrix}

# With MatrixMult operator filled with zeros
Zop = pylops.MatrixMult(np.zeros(D2vop.shape))
Bop = pylops.Block([[D2vop, Zop, -1 * D2hop], [D2hop, 2 * D2hop, Zop]])
Y = np.reshape(Bop * X.ravel(), (2 * Nv, Nh))

# With Zero operator
Zop = pylops.Zero(D2vop.shape[0])
B1op = pylops.Block([[D2vop, Zop, -1 * D2hop], [D2hop, 2 * D2hop, Zop]])
Y1 = np.reshape(B1op * X.ravel(), (2 * Nv, Nh))

print("Y == Y1:", np.allclose(Y, Y1))

###############################################################################
# Finally we can use the *block-diagonal operator* to apply three operators
# to three different subset of the model and data
#
# .. math::
# \mathbf{D_{BDiag}} =
# \begin{bmatrix}
# \mathbf{D_v} & \mathbf{0} & \mathbf{0} \\
# \mathbf{0} & 0.5*\mathbf{D_v} & \mathbf{0} \\
# \mathbf{0} & \mathbf{0} & -\mathbf{D_h}
# \end{bmatrix}, \qquad
# \mathbf{y} =
# \begin{bmatrix}
# \mathbf{D_v} \mathbf{x_1} \\
# 0.5*\mathbf{D_v} \mathbf{x_2} \\
# -\mathbf{D_h} \mathbf{x_3}
# \end{bmatrix}
# .. math::
# \mathbf{D_{BDiag}} =
# \begin{bmatrix}
# \mathbf{D_v} & \mathbf{0} & \mathbf{0} \\
# \mathbf{0} & 0.5*\mathbf{D_v} & \mathbf{0} \\
# \mathbf{0} & \mathbf{0} & -\mathbf{D_h}
# \end{bmatrix}, \qquad
# \mathbf{y} =
# \begin{bmatrix}
# \mathbf{D_v} \mathbf{x_1} \\
# 0.5*\mathbf{D_v} \mathbf{x_2} \\
# -\mathbf{D_h} \mathbf{x_3}
# \end{bmatrix}

BD = pylops.BlockDiag([D2vop, 0.5 * D2vop, -1 * D2hop])
Y = np.reshape(BD * X.ravel(), (3 * Nv, Nh))

Expand All @@ -195,6 +264,7 @@
# blockdiagonal structure, it may be convenient to span multiple processes
# handling subset of operators at the same time. This can be easily achieved
# by simply defining the number of processes we want to use via ``nproc``.

X = np.zeros((Nv * 10, Nh))
for iv in range(10):
X[int(Nv / 2) + iv * Nv, int(Nh / 2)] = 1
Expand All @@ -221,21 +291,21 @@
# Finally we use the *Kronecker operator* and replicate this example on
# `wiki <https://en.wikipedia.org/wiki/Kronecker_product>`_.
#
# .. math::
# \begin{bmatrix}
# 1 & 2 \\
# 3 & 4 \\
# \end{bmatrix} \otimes
# \begin{bmatrix}
# 0 & 5 \\
# 6 & 7 \\
# \end{bmatrix} =
# \begin{bmatrix}
# 0 & 5 & 0 & 10 \\
# 6 & 7 & 12 & 14 \\
# 0 & 15 & 0 & 20 \\
# 18 & 21 & 24 & 28 \\
# \end{bmatrix}
# .. math::
# \begin{bmatrix}
# 1 & 2 \\
# 3 & 4 \\
# \end{bmatrix} \otimes
# \begin{bmatrix}
# 0 & 5 \\
# 6 & 7 \\
# \end{bmatrix} =
# \begin{bmatrix}
# 0 & 5 & 0 & 10 \\
# 6 & 7 & 12 & 14 \\
# 0 & 15 & 0 & 20 \\
# 18 & 21 & 24 & 28 \\
# \end{bmatrix}
A = np.array([[1, 2], [3, 4]])
B = np.array([[0, 5], [6, 7]])
AB = np.kron(A, B)
Expand Down
4 changes: 4 additions & 0 deletions pylops/basicoperators/block.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,10 @@ class Block(_Block):
r"""Block operator.

Create a block operator from N lists of M linear operators each.
Note that in case one or more operators are filled with zeros, it is
recommended to use the :py:class:`pylops.Zero` operator instead of e.g.,
:py:class:`pylops.MatrixMult` with a matrix of zeros, as the former will
be simply by-passed both in the forward and adjoint steps.

Parameters
----------
Expand Down
31 changes: 19 additions & 12 deletions pylops/basicoperators/hstack.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@
from typing import Optional, Sequence

from pylops import LinearOperator
from pylops.basicoperators import MatrixMult
from pylops.basicoperators import MatrixMult, Zero
from pylops.utils.backend import get_array_module, get_module, inplace_add, inplace_set
from pylops.utils.typing import NDArray

Expand All @@ -33,7 +33,12 @@ def _matvec_rmatvec_map(op, x: NDArray) -> NDArray:
class HStack(LinearOperator):
r"""Horizontal stacking.

Stack a set of N linear operators horizontally.
Stack a set of N linear operators horizontally. Note that in case
one or more operators are filled with zeros, it is recommended to use
the :py:class:`pylops.Zero` operator instead of e.g.,
:py:class:`pylops.MatrixMult` with a matrix of zeros, as the former will
be simply by-passed both in the forward and adjoint steps.


Parameters
----------
Expand Down Expand Up @@ -178,11 +183,12 @@ def _matvec_serial(self, x: NDArray) -> NDArray:
)
y = ncp.zeros(self.nops, dtype=self.dtype)
for iop, oper in enumerate(self.ops):
y = inplace_add(
oper.matvec(x[self.mmops[iop] : self.mmops[iop + 1]]).squeeze(),
y,
slice(None, None),
)
if not isinstance(oper, Zero):
y = inplace_add(
oper.matvec(x[self.mmops[iop] : self.mmops[iop + 1]]).squeeze(),
y,
slice(None, None),
)
return y

def _rmatvec_serial(self, x: NDArray) -> NDArray:
Expand All @@ -193,11 +199,12 @@ def _rmatvec_serial(self, x: NDArray) -> NDArray:
)
y = ncp.zeros(self.mops, dtype=self.dtype)
for iop, oper in enumerate(self.ops):
y = inplace_set(
oper.rmatvec(x).squeeze(),
y,
slice(self.mmops[iop], self.mmops[iop + 1]),
)
if not isinstance(oper, Zero):
y = inplace_set(
oper.rmatvec(x).squeeze(),
y,
slice(self.mmops[iop], self.mmops[iop + 1]),
)
return y

def _matvec_multiproc(self, x: NDArray) -> NDArray:
Expand Down
28 changes: 18 additions & 10 deletions pylops/basicoperators/vstack.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@
from typing import Callable, Optional, Sequence

from pylops import LinearOperator
from pylops.basicoperators import MatrixMult
from pylops.basicoperators import MatrixMult, Zero
from pylops.utils.backend import get_array_module, get_module, inplace_add, inplace_set
from pylops.utils.typing import DTypeLike, NDArray

Expand All @@ -33,7 +33,11 @@ def _matvec_rmatvec_map(op: Callable, x: NDArray) -> NDArray:
class VStack(LinearOperator):
r"""Vertical stacking.

Stack a set of N linear operators vertically.
Stack a set of N linear operators vertically. Note that in case
one or more operators are filled with zeros, it is recommended to use
the :py:class:`pylops.Zero` operator instead of e.g.,
:py:class:`pylops.MatrixMult` with a matrix of zeros, as the former will
be simply by-passed both in the forward and adjoint steps.

Parameters
----------
Expand Down Expand Up @@ -178,9 +182,12 @@ def _matvec_serial(self, x: NDArray) -> NDArray:
)
y = ncp.zeros(self.nops, dtype=self.dtype)
for iop, oper in enumerate(self.ops):
y = inplace_set(
oper.matvec(x).squeeze(), y, slice(self.nnops[iop], self.nnops[iop + 1])
)
if not isinstance(oper, Zero):
y = inplace_set(
oper.matvec(x).squeeze(),
y,
slice(self.nnops[iop], self.nnops[iop + 1]),
)
return y

def _rmatvec_serial(self, x: NDArray) -> NDArray:
Expand All @@ -191,11 +198,12 @@ def _rmatvec_serial(self, x: NDArray) -> NDArray:
)
y = ncp.zeros(self.mops, dtype=self.dtype)
for iop, oper in enumerate(self.ops):
y = inplace_add(
oper.rmatvec(x[self.nnops[iop] : self.nnops[iop + 1]]).squeeze(),
y,
slice(None, None),
)
if not isinstance(oper, Zero):
y = inplace_add(
oper.rmatvec(x[self.nnops[iop] : self.nnops[iop + 1]]).squeeze(),
y,
slice(None, None),
)
return y

def _matvec_multiproc(self, x: NDArray) -> NDArray:
Expand Down
Loading
Loading