pyxu.operator.linop #

Sum(arg_shape, axis=None)[source]#

Multi-dimensional sum reduction.

This operator re-arranges the input array to a multi-dimensional array of shape arg_shape, then reduces it via summation across one or more axis.

For example, assuming the input array \(\mathbf{x} \in \mathbb{R}^{N_1 \times N_2 \times N_3}\) and axis=-1, then

\[\mathbf{y}_{i,j} = \sum_{k}{\mathbf{x}_{i,j,k}}.\]

Parameters:

arg_shape (NDArrayShape) – Shape of the data to be reduced.
axis (NDArrayAxis) – Axis or axes along which a sum is performed. The default, axis=None, will sum all the elements of the input array. If axis is negative it counts from the last to the first axis.

Return type:

OpT

Notes

The matrix operator of a 1D reduction applied to \(\mathbf{x} \in \mathbb{R}^{N}\) is given by

\[\mathbf{A}(x) = \mathbf{1}^{T} \mathbf{x},\]

where \(\sigma_{\max}(\mathbf{A}) = \sqrt{N}\). An ND reduction is a chain of 1D reductions in orthogonal dimensions. Hence the Lipschitz constant of an ND reduction is the product of Lipschitz constants of all 1D reductions involved, i.e.:

\[L = \sqrt{\prod_{i_{k}} N_{i_{k}}},\]

where \(\{i_{k}\}_{k}\) denotes the axes being summed over.

class IdentityOp(dim)[source]#

Bases: OrthProjOp

Identity operator.

Parameters:: dim (Integral) –

class NullOp(shape)[source]#

Bases: LinOp

Null operator.

This operator maps any input vector on the null vector.

Parameters:: shape (tuple[Integral, Integral]) –

NullFunc(dim)[source]#

Null functional.

This functional maps any input vector on the null scalar.

Parameters:: dim (Integral) –
Return type:: OpT

HomothetyOp(dim, cst)[source]#

Constant scaling operator.

Parameters:

dim (Integer) – Dimension of the domain.
cst (Real) – Scaling factor.

Returns:

op – (dim, dim) scaling operator.

Return type:

Note

This operator is not defined in terms of DiagonalOp() since it is array-backend-agnostic.

DiagonalOp(vec, enable_warnings=True)[source]#

Element-wise scaling operator.

Note

DiagonalOp() instances are not arraymodule-agnostic: they will only work with NDArrays belonging to the same array module as vec. Moreover, inner computations may cast input arrays when the precision of vec does not match the user-requested precision. If such a situation occurs, a warning is raised.

Parameters:

vec (NDArray) – (N,) diagonal scale factors.
enable_warnings (bool) – If True, emit a warning in case of precision mis-match issues.

Return type:

OpT

class Pad(arg_shape, pad_width, mode='constant')[source]#

Bases: LinOp

Multi-dimensional padding operator.

This operator pads the input array in each dimension according to specified widths.

Notes

If inputs are D-dimensional, then some of the padding of later axes are calculated from padding of previous axes.
The adjoint of the padding operator performs a cumulative summation over the original positions used to pad. Its effect is clear from its matrix form. For example the matrix-form of Pad(arg_shape=(3,), mode="wrap", pad_width=(1, 1)) is:

\[\begin{split}\mathbf{A} = \left[ \begin{array}{ccc} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{array} \right].\end{split}\]

The adjoint of \(\mathbf{A}\) corresponds to its matrix transpose:

\[\begin{split}\mathbf{A}^{\ast} = \left[ \begin{array}{ccccc} 0 & 1 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 \end{array} \right].\end{split}\]

This operation can be seen as a trimming (\(\mathbf{T}\)) plus a cumulative summation (\(\mathbf{S}\)):

\[\begin{split}\mathbf{A}^{\ast} = \mathbf{T} + \mathbf{S} = \left[ \begin{array}{ccccc} 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \end{array} \right] + \left[ \begin{array}{ccccc} 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \end{array} \right],\end{split}\]

where both \(\mathbf{T}\) and \(\mathbf{S}\) are efficiently implemented in matrix-free form.
The Lipschitz constant of the multi-dimensional padding operator is upper-bounded by the product of Lipschitz constants of the uni-dimensional paddings applied per dimension, i.e.:

\[L \le \prod_{i} L_{i}, \qquad i \in \{0, \dots, N-1\},\]

where \(L_{i}\) depends on the boundary condition at the \(i\)-th axis.

\(L_{i}^{2}\) corresponds to the maximum singular value of the diagonal matrix

\[\mathbf{A}_{i}^{\ast} \mathbf{A}_{i} = \mathbf{T}_{i}^{\ast} \mathbf{T}_{i} + \mathbf{S}_{i}^{\ast} \mathbf{S}_{i} = \mathbf{I}_{N} + \mathbf{S}_{i}^{\ast} \mathbf{S}_{i}.\]
- In mode=”constant”, \(\text{diag}(\mathbf{S}_{i}^{\ast} \mathbf{S}_{i}) = \mathbf{0}\), hence \(L_{i} = 1\).
- In mode=”edge”,
  
  \[\text{diag}(\mathbf{S}_{i}^{\ast} \mathbf{S}_{i}) = \left[p_{lhs}, 0, \ldots, 0, p_{rhs} \right],\]
  
  hence \(L_{i} = \sqrt{1 + \min(p_{lhs}, p_{rhs})}\).
- In mode=”symmetric”, “wrap”, “reflect”, \(\text{diag}(\mathbf{S}_{i}^{\ast} \mathbf{S}_{i})\) equals (up to a mode-dependant permutation)
  
  \[\text{diag}(\mathbf{S}_{i}^{\ast} \mathbf{S}_{i}) = \left[1, \ldots, 1, 0, \ldots, 0\right] + \left[0, \ldots, 0, 1, \ldots, 1\right],\]
  
  hence
  
  \[\begin{split}L^{\text{wrap, symmetric}}_{i} = \sqrt{1 + \lceil\frac{p_{lhs} + p_{rhs}}{N}\rceil}, \\ L^{\text{reflect}}_{i} = \sqrt{1 + \lceil\frac{p_{lhs} + p_{rhs}}{N-2}\rceil}.\end{split}\]

Parameters:

arg_shape (Integral | tuple[Integral, ...]) –
pad_width (Integral | Sequence[Integral] | Sequence[tuple[Integral, Integral]]) –
mode (str | Sequence[str]) –

WidthSpec#: alias of Union[Integral, Sequence[Integral], Sequence[tuple[Integral, Integral]]]

ModeSpec#: alias of Union[str, Sequence[str]]

__init__(arg_shape, pad_width, mode='constant')[source]#

Parameters:

arg_shape (NDArrayShape) – Shape of the input array.
pad_width (WidthSpec) –
Number of values padded to the edges of each axis. Multiple forms are accepted:
- int: pad each dimension’s head/tail by pad_width.
- tuple[int, ...]: pad dimension[k]’s head/tail by pad_width[k].
- tuple[tuple[int, int], ...]: pad dimension[k]’s head/tail by pad_width[k][0] / pad_width[k][1] respectively.
mode (str, list ( str )) –
Padding mode. Multiple forms are accepted:
- str: unique mode shared amongst dimensions. Must be one of:
  - ’constant’ (zero-padding)
  - ’wrap’
  - ’reflect’
  - ’symmetric’
  - ’edge’
- tuple[str, …]: pad dimension[k] using mode[k].
(See numpy.pad() for details.)

Transforms #

class FFT(arg_shape, axes=None, real=False, **kwargs)[source]#

Bases: LinOp

Multi-dimensional Discrete Fourier Transform (DFT) \(A: \mathbb{C}^{N_{1} \times \cdots \times N_{D}} \to \mathbb{C}^{N_{1} \times \cdots \times N_{D}}\).

The FFT is defined as follows:

\[(A \, \mathbf{x})[\mathbf{k}] = \sum_{\mathbf{n}} \mathbf{x}[\mathbf{n}] \exp\left[-j 2 \pi \langle \frac{\mathbf{n}}{\mathbf{N}}, \mathbf{k} \rangle \right],\]

\[(A^{*} \, \hat{\mathbf{x}})[\mathbf{n}] = \sum_{\mathbf{k}} \hat{\mathbf{x}}[\mathbf{k}] \exp\left[j 2 \pi \langle \frac{\mathbf{n}}{\mathbf{N}}, \mathbf{k} \rangle \right],\]

\[(\mathbf{x}, \, \hat{\mathbf{x}}) \in \mathbb{C}^{N_{1} \times \cdots \times N_{D}}, \quad (\mathbf{n}, \, \mathbf{k}) \in \{0, \ldots, N_{1}-1\} \times \cdots \times \{0, \ldots, N_{D}-1\}.\]

The DFT is taken over any number of axes by means of the Fast Fourier Transform algorithm (FFT).

Implementation Notes

The CPU implementation uses SciPy’s FFT implementation.
The GPU implementation uses cuFFT via CuPy.

Examples

1D DFT of a cosine pulse.

from pyxu.operator import FFT
import pyxu.util as pxu

N = 10
op = FFT(
    arg_shape=(N,),
    real=True,
)

x = np.cos(2 * np.pi / N * np.arange(N))
y = pxu.view_as_complex(op.apply(x))
# [0, N/2, 0, 0, 0, 0, 0, 0, 0, N/2]

z = op.adjoint(op.apply(x))
# np.allclose(z, N * x) -> True

1D DFT of a complex exponential pulse.

from pyxu.operator import FFT
import pyxu.util as pxu

N = 10
op = FFT(
    arg_shape=(N,),
    real=False,  # complex-valued inputs
)

x = np.exp(1j * 2 * np.pi / N * np.arange(N))
y = pxu.view_as_complex(
        op.apply(
            pxu.view_as_real(x)
        )
    )
# [0, N, 0, 0, 0, 0, 0, 0, 0, 0]

z = pxu.view_as_complex(
        op.adjoint(
            op.apply(
                pxu.view_as_real(x)
            )
        )
    )
# np.allclose(z, N * x) -> True

2D DFT of an image

from pyxu.operator import FFT
import pyxu.util as pxu

N_h, N_w = 10, 8
op = FFT(
    arg_shape=(N_h, N_w),
    real=True,
)

x = np.pad(  # an image
    np.ones((N_h//2, N_w//2)),
    pad_width=((0, N_h//2), (0, N_w//2)),
)
y = pxu.view_as_complex(  # sinc
        op.apply(x.reshape(-1))
    ).reshape(N_h, N_w)
z = op.adjoint(
        op.apply(x.reshape(-1))
    ).reshape(N_h, N_w)
# np.allclose(z, (N_h * N_w) * x) -> True

1D DFT of an image’s rows

from pyxu.operator import FFT
import pyxu.util as pxu

N_h, N_w = 10, 8
op = FFT(
    arg_shape=(N_h, N_w),
    axes=-1,
    real=True,
)

x = np.pad(  # an image
    np.ones((N_h//2, N_w//2)),
    pad_width=((0, N_h//2), (0, N_w//2)),
)
y = pxu.view_as_complex(  # sinc
        op.apply(x.reshape(-1))
    ).reshape(N_h, N_w)
z = op.adjoint(
        op.apply(x.reshape(-1))
    ).reshape(N_h, N_w)
# np.allclose(z, N_w * x) -> True

Parameters:

arg_shape (Integral | tuple[Integral, ...]) –
axes (Integral | tuple[Integral, ...]) –
real (bool) –

__init__(arg_shape, axes=None, real=False, **kwargs)[source]#

Parameters:

arg_shape (NDArrayShape) – (N_1, …, N_D) shape of the input array \(\mathbf{x} \in \mathbb{R}^{N_{1} \times \cdots \times N_{D}}\) or \(\mathbb{C}^{N_{1} \times \cdots \times N_{D}}\).
axes (NDArrayAxis) – Axis or axes along which the DFT is performed. The default, axes=None, will transform all dimensions of the input array.
real (bool) –
If True, assumes apply() takes (…, N.prod()) inputs in \(\mathbb{R}^{N}\).

If False, then apply() takes (…, 2N.prod()) inputs, i.e. \(\mathbb{C}^{N}\) vectors viewed as bijections with \(\mathbb{R}^{2N}\).
kwargs (dict) –
Extra kwargs passed to scipy.fft.fftn() or cupyx.scipy.fft.fftn().

Supported parameters for scipy.fft.fftn() are:
- workers: int = 1
Supported parameters for cupyx.scipy.fft.fftn() are:
- NOT SUPPORTED FOR NOW
Default values are chosen if unspecified.

apply(arr)[source]#

Parameters:

arr (NDArray) –

(…, N.prod()) inputs \(\mathbf{x} \in \mathbb{R}^{N_{1} \times \cdots \times N_{D}}\) (real=True.)
(…, 2N.prod()) inputs \(\mathbf{x} \in \mathbb{C}^{N_{1} \times \cdots \times N_{D}}\) viewed as a real array. (See view_as_real().)

Returns:

out – (…, 2N.prod()) outputs \(\hat{\mathbf{x}} \in \mathbb{C}^{N_{1} \times \cdots \times N_{D}}\) viewed as a real array. (See view_as_real().)

Return type:

adjoint(arr)[source]#

Parameters:

arr (NDArray) – (…, 2N.prod()) outputs \(\hat{\mathbf{x}} \in \mathbb{C}^{N_{1} \times \cdots \times N_{D}}\) viewed as a real array. (See view_as_real().)

Returns:

out –

(…, N.prod()) inputs \(\mathbf{x} \in \mathbb{R}^{N_{1} \times \cdots \times N_{D}}\) (real=True.)
(…, 2N.prod()) inputs \(\mathbf{x} \in \mathbb{C}^{N_{1} \times \cdots \times N_{D}}\) viewed as a real array. (See view_as_real().)

Return type:

class NUFFT(shape)[source]#

Bases: LinOp

Non-Uniform Fast Fourier Transform (NUFFT) of Type 1/2/3 (for \(d=\{1,2,3\}\)).

The Non-Uniform Fast Fourier Transform (NUFFT) generalizes the FFT to off-grid data. There are three types of NUFFTs proposed in the literature:

Type 1 (non-uniform to uniform),
Type 2 (uniform to non-uniform),
Type 3 (non-uniform to non-uniform).

See the notes below, including [FINUFFT], for definitions of the various transforms and algorithmic details.

The transforms should be instantiated via type1(), type2(), and type3() respectively. (See each method for usage examples.)

The dimension of each transform is inferred from the dimensions of the input arguments, with support for \(d=\{1,2,3\}\).

Notes

We adopt here the same notational conventions as in [FINUFFT].

Mathematical Definition. Let \(d\in\{1,2,3\}\) and consider the mesh

\[\mathcal{I}_{N_1,\ldots,N_d} = \mathcal{I}_{N_1} \times \cdots \times \mathcal{I}_{N_d} \subset \mathbb{Z}^d,\]

where the mesh indices \(\mathcal{I}_{N_i}\subset\mathbb{Z}\) are given for each dimension \(i=1,\dots, d\) by:

\[\begin{split}\mathcal{I}_{N_i} = \begin{cases} [[-N_i/2, N_i/2-1]], & N_i\in 2\mathbb{N} \text{ (even)}, \\ [[-(N_i-1)/2, (N_i-1)/2]], & N_i\in 2\mathbb{N}+1 \text{ (odd)}. \end{cases}\end{split}\]

Then the NUFFT operators approximate, up to a requested relative accuracy \(\varepsilon \geq 0\), [1] the following exponential sums:

\[\begin{split}\begin{align} (1)\; &u_{\mathbf{n}} = \sum_{j=1}^{M} w_{j} e^{\sigma i\langle \mathbf{n}, \mathbf{x}_{j} \rangle}, \quad &\mathbf{n}\in \mathcal{I}_{N_1,\ldots, N_d},\qquad &\text{Type 1 (non-uniform to uniform)}\\ (2)\; &w_{j} = \sum_{\mathbf{n}\in\mathcal{I}_{N_1,\ldots, N_d}} u_{\mathbf{n}} e^{\sigma i\langle \mathbf{n}, \mathbf{x}_{j} \rangle }, \quad &j=1,\ldots, M,\qquad &\text{Type 2 (uniform to non-uniform)}\\ (3)\; &v_{k} = \sum_{j=1}^{M} w_{j} e^{\sigma i\langle \mathbf{z}_k, \mathbf{x}_{j} \rangle }, \quad &k=1,\ldots, N, \qquad &\text{Type 3 (non-uniform to non-uniform)} \end{align}\end{split}\]

where \(\sigma \in \{+1, -1\}\) defines the sign of the transform and \(u_{\mathbf{n}}, v_{k}, w_{j}\in \mathbb{C}\). For the type-1 and type-2 NUFFTs, the non-uniform samples \(\mathbf{x}_{j}\) are assumed to lie in \([-\pi,\pi)^{d}\). For the type-3 NUFFT, the non-uniform samples \(\mathbf{x}_{j}\) and \(\mathbf{z}_{k}\) are arbitrary points in \(\mathbb{R}^{d}\).
Adjoint NUFFTs. The type-1 and type-2 NUFFTs with opposite signs form an adjoint pair. The adjoint of the type-3 NUFFT is obtained by flipping the transform’s sign and switching the roles of \(\mathbf{z}_k\) and \(\mathbf{x}_{j}\) in (3).
Lipschitz Constants. We bound the Lipschitz constant by the Frobenius norm of the operators, which yields \(L \leq \sqrt{NM}\). Note that these Lipschitz constants are cheap to compute but may be pessimistic. Tighter Lipschitz constants can be computed by calling estimate_lipschitz().
Error Analysis. Let \(\tilde{\mathbf{u}}\in\mathbb{C}^{\mathcal{I}_{N_1,\ldots, N_d}}\) and \(\tilde{\mathbf{w}}\in\mathbb{C}^{M}\) be the outputs of the type-1 and type-2 NUFFT algorithms which approximate the sequences \({\mathbf{u}}\in\mathbb{C}^{\mathcal{I}_{N_1,\ldots, N_d}}\) and \({\mathbf{w}}\in\mathbb{C}^{M}\) defined in (1) and (2) respectively. Then [FINUFFT] shows that the relative errors \(\|\tilde{\mathbf{u}}-{\mathbf{u}}\|_2/\|{\mathbf{u}}\|_2\) and \(\|\tilde{\mathbf{w}}-{\mathbf{w}}\|_2/\|{\mathbf{w}}\|_2\) are almost always similar to the user-requested tolerance \(\varepsilon\), except when round-off error dominates (i.e. very small user-requested tolerances). The same holds approximately for the type-3 NUFFT. Note however that this is a typical error analysis: some degenerate (but rare) worst-case scenarios can result in higher errors.
Complexity. Naive evaluation of the exponential sums (1), (2) and (3) above costs \(O(NM)\), where \(N=N_{1}\ldots N_{d}\) for the type-1 and type-2 NUFFTs. NUFFT algorithms approximate these sums to a user-specified relative tolerance \(\varepsilon\) in log-linear complexity in \(N\) and \(M\). The complexity of the various NUFFTs are given by (see [FINUFFT]):

\[\begin{split}&\mathcal{O}\left(N \log(N) + M|\log(\varepsilon)|^d\right)\qquad &\text{(Type 1 and 2)}\\ &\mathcal{O}\left(\Pi_{i=1}^dX_iZ_i\sum_{i=1}^d\log(X_iZ_i) + (M + N)|\log(\varepsilon)|^d\right)\qquad &\text{(Type 3)}\end{split}\]

where \(X_i = \max_{j=1,\ldots,M}|(\mathbf{x}_{j})_i|\) and \(Z_i = \max_{k=1,\ldots,N}|(\mathbf{z}_k)_i|\) for \(i=1,\ldots,d\). The terms above correspond to the complexities of the FFT and spreading/interpolation steps respectively.
Memory footprint. The complexity and memory footprint of the type-3 NUFFT can be arbitrarily large for poorly-centered data, or for data with a large spread. An easy fix consists in centering the data before/after the NUFFT via pre/post-phasing operations, as described in equation (3.24) of [FINUFFT]. This optimization is automatically carried out by FINUFFT if the compute/memory gains are non-negligible.

Additionally the type-3 summation (eq. 3) can be evaluated block-wise by partitioning the non-uniform samples \((x_{j}, z_{k})\):

\[\begin{align} (4)\;\; &v_{k} = \sum_{p=1}^{P} \sum_{j\in \mathcal{M}_{p}} w_{j} e^{\sigma i\langle \mathbf{z}_{k}, \mathbf{x}_{j} \rangle }, \quad & k\in \mathcal{N}_{q}, \quad q=1,\ldots, Q, \quad & \text{Type 3 (chunked)} \end{align}\]

where \(\{\mathcal{M}_1, \ldots, \mathcal{M}_P\}\) and \(\{\mathcal{N}_1, \ldots, \mathcal{N}_Q\}\) are partitions of the sets \(\{1, \ldots, M\}\) and \(\{1, \ldots, N\}\) respectively. In the chunked setting, the partial sums in (4)

\[\left\{\sum_{j\in \mathcal{M}_p} w_{j} e^{\sigma i\langle \mathbf{z}_k, \mathbf{x}_{j} \rangle }, \, k\in \mathcal{N}_q\right\}_{p,q}\]

are each evaluated via a type-3 NUFFT involving a subset of the non-uniform samples. Sub-problems can be evaluated in parallel. Moreover, for wisely chosen data partitions, the memory budget of each NUFFT can be explicitly controlled and capped to a maximal value. This allows one to perform out-of-core NUFFTs with (theoretically) no complexity overhead w.r.t a single large NUFFT. (See note below however.)

Chunking of the type-3 NUFFT is activated by passing chunked=True to type3() (together with parallel=True for parallel computations). Finally, auto_chunk() can be used to compute a good partition of the X/Z-domains.

Hint

The current implementation of the chunked type-3 NUFFT has computational complexity:

\[\mathcal{O}\left( \sum_{p,q=1}^{P}\Pi_{i=1}^{d} X^{(p)}_{i}Z^{(q)}_{i} \sum_{i=1}^{d} \log(X^{(p)}_{i} Z^{(q)}_{i}) + (M_p + N_q)|\log(\varepsilon)|^{d} \right),\]

where

\[\begin{split}\begin{align*} X^{(p)}_{i} & = \max_{j\in \mathcal{M}_{p}}|(\mathbf{x}_{j})_{i}|, \quad i = \{1,\ldots,d\} \\ Z^{(p)}_{i} & = \max_{k\in \mathcal{N}_{q}}|(\mathbf{z}_k)_{i}|, \quad i = \{1,\ldots,d\} \end{align*}\end{split}\]

and \(M_{p}, N_{q}\) denote the number of elements in the sets \(\mathcal{M}_{p}, \mathcal{N}_{q}\), respectively.

For perfectly balanced and uniform chunks (i.e. \(M_{p}=M/P\), \(N_{q}=N/Q\), \(X^{(p)}_{i} = X_{i}/P\), \(Z^{(q)}_{i} = Z_{i}/Q\) and \(P=Q\)) the complexity reduces to

\[\mathcal{O}\left( \Pi_{i=1}^{d} X_{i}Z_{i}\sum_{i=1}^{d}\log(X_{i}Z_{i}) + P(M + N)|\log(\varepsilon)|^{d} \right).\]

This shows that the spreading/interpolation is \(P\) times more expensive than in the non-chunked case. This overhead is usually acceptable if the number of chunks remains small. With explicit control on the spreading/interpolation steps (which the Python interface to the FINUFFT backend does not currently offer), the spreading/interpolation overhead can be removed so as to obtain a computational complexity on par to that of the non-chunked type-3 NUFFT. This will be implemented in a future release.
Backend. The NUFFT transforms are computed via Python wrappers to FINUFFT and cuFINUFFT. (See also [FINUFFT] and [cuFINUFFT].) These librairies perform the expensive spreading/interpolation between nonuniform points and the fine grid via the “exponential of semicircle” kernel.
Optional Parameters. [cu]FINUFFT exposes many optional parameters to adjust the performance of the algorithms, change the output format, or provide debug/timing information. While the default options are sensible for most setups, advanced users may overwrite them via the kwargs parameter of type1(), type2(), and type3(). See the guru interface from FINUFFT and its companion page for details.

Hint

FINUFFT exposes a dtype keyword to control the precision (single or double) at which transforms are performed. This parameter is ignored by NUFFT: use the context manager Precision to control floating point precision.

Hint

The NUFFT is performed in batches of (n_trans,), where n_trans denotes the number of simultaneous transforms requested. (See the n_trans parameter of finufft.Plan.)

Good performance is obtained when each batch fits easily in memory. This recommendation also applies to Dask inputs which are re-chunked internally to be n_trans-compliant.

This parameter also affects performance of the eps=0 case: increasing n_trans may improve performance when doing several transforms in parallel.

Warning

Since FINUFFT plans cannot be shared among different processes, this class is only compatible with Dask’s thread-based schedulers, i.e.:

scheduler='thread'
scheduler='synchronous'
distributed.Client(processes=False)

Batches are hence processed serially. (Each batch is multi-threaded however.)

See also

FFT

Parameters:: shape (tuple[Integral, Integral]) –

static type1(x, N, isign=1, eps=0.0001, real=False, plan_fw=True, plan_bw=True, enable_warnings=True, **kwargs)[source]#

Type-1 NUFFT (non-uniform to uniform).

Parameters:

x (NDArray) – (M, [d]) d-dimensional sample points \(\mathbf{x}_{j} \in [-\pi,\pi)^{d}\).
N (Integer, tuple ( Integer )) –
([d],) mesh size in each dimension \((N_1, \ldots, N_d)\).

If N is an integer, then the mesh is assumed to have the same size in each dimension.
isign (1, -1) – Sign \(\sigma\) of the transform.
eps (Real) –
Requested relative accuracy \(\varepsilon \geq 0\).

If eps=0, the transform is computed exactly via direct evaluation of the exponential sum using a Numba JIT-compiled kernel.
real (bool) –
If True, assumes apply() takes (…, M) inputs in \(\mathbb{R}^{M}\).

If False, then apply() takes (…, 2M) inputs, i.e. \(\mathbb{C}^{M}\) vectors viewed as bijections with \(\mathbb{R}^{2M}\).
plan_fw/bw (bool) – If True, allocate FINUFFT resources to do the forward (fw) and/or backward (bw) transform. These are advanced options: use them with care. Some public methods in the LinOp interface may not work if fw/bw transforms are disabled. These options only take effect if eps > 0.
enable_warnings (bool) – If True, emit a warning in case of precision mis-match issues.
**kwargs –
Extra kwargs to finufft.Plan. (Illegal keywords are dropped silently.) Most useful are n_trans, nthreads and debug.
plan_fw (bool) –
plan_bw (bool) –

Returns:

op – (2N.prod(), M) or (2N.prod(), 2M) type-1 NUFFT.

Return type:

Examples

import numpy as np
import pyxu.operator as pxo
import pyxu.runtime as pxrt
import pyxu.util as pxu

rng = np.random.default_rng(0)
D, M, N = 2, 200, 5  # D denotes the dimension of the data
x = np.fmod(rng.normal(size=(M, D)), 2 * np.pi)

with pxrt.Precision(pxrt.Width.SINGLE):
    # The NUFFT dimension (1/2/3) is inferred from the trailing dimension of x.
    # Its precision is controlled by the context manager.
    N_trans = 5
    A = pxo.NUFFT.type1(
            x, N,
            n_trans=N_trans,
            isign=-1,
            eps=1e-3,
        )

    # Pyxu operators only support real inputs/outputs, so we use the functions
    # pxu.view_as_[complex/real] to interpret complex arrays as real arrays (and
    # vice-versa).
    arr =        rng.normal(size=(3, N_trans, M)) \
          + 1j * rng.normal(size=(3, N_trans, M))
    A_out_fw = pxu.view_as_complex(A.apply(pxu.view_as_real(arr)))
    A_out_bw = pxu.view_as_complex(A.adjoint(pxu.view_as_real(A_out_fw)))

static type2(x, N, isign=1, eps=0.0001, real=False, plan_fw=True, plan_bw=True, enable_warnings=True, **kwargs)[source]#

Type-2 NUFFT (uniform to non-uniform).

Parameters:

x (NDArray) – (M, [d]) d-dimensional sample points \(\mathbf{x}_{j} \in [-\pi,\pi]^{d}\).
N (Integer, tuple ( Integer )) –
([d],) mesh size in each dimension \((N_1, \ldots, N_d)\).

If N is an integer, then the mesh is assumed to have the same size in each dimension.
isign (1, -1) – Sign \(\sigma\) of the transform.
eps (Real) –
Requested relative accuracy \(\varepsilon \geq 0\).

If eps=0, the transform is computed exactly via direct evaluation of the exponential sum using a Numba JIT-compiled kernel.
real (bool) –
If True, assumes apply() takes (…, N.prod()) inputs in \(\mathbb{R}^{N}\).

If False, then apply() takes (…, 2N.prod()) inputs, i.e. \(\mathbb{C}^{N}\) vectors viewed as bijections with \(\mathbb{R}^{2N}\).
plan_fw/bw (bool) – If True, allocate FINUFFT resources to do the forward (fw) and/or backward (bw) transform. These are advanced options: use them with care. Some public methods in the LinOp interface may not work if fw/bw transforms are disabled. These options only take effect if eps > 0.
enable_warnings (bool) – If True, emit a warning in case of precision mis-match issues.
**kwargs –
Extra kwargs to finufft.Plan. (Illegal keywords are dropped silently.) Most useful are n_trans, nthreads and debug.
plan_fw (bool) –
plan_bw (bool) –

Returns:

op – (2M, N.prod()) or (2M, 2N.prod()) type-2 NUFFT.

Return type:

Examples

import numpy as np
import pyxu.operator as pxo
import pyxu.runtime as pxrt
import pyxu.util as pxu

rng = np.random.default_rng(0)
D, M, N = 2, 200, 5  # D denotes the dimension of the data
N_full = (N,) * D
x = np.fmod(rng.normal(size=(M, D)), 2 * np.pi)

with pxrt.Precision(pxrt.Width.SINGLE):
    # The NUFFT dimension (1/2/3) is inferred from the trailing dimension of x.
    # Its precision is controlled by the context manager.
    N_trans = 5
    A = pxo.NUFFT.type2(
            x, N,
            n_trans=N_trans,
            isign=-1,
            eps=1e-3,
        )

    # Pyxu operators only support real inputs/outputs, so we use the functions
    # pxu.view_as_[complex/real] to interpret complex arrays as real arrays (and
    # vice-versa).
    arr = np.reshape(
               rng.normal(size=(3, N_trans, *N_full))
        + 1j * rng.normal(size=(3, N_trans, *N_full)),
        (3, N_trans, -1),
    )
    A_out_fw = pxu.view_as_complex(A.apply(pxu.view_as_real(arr)))
    A_out_bw = pxu.view_as_complex(A.adjoint(pxu.view_as_real(A_out_fw)))

static type3(

x,

z,

isign=1,

eps=0.0001,

real=False,

plan_fw=True,

plan_bw=True,

enable_warnings=True,

chunked=False,

parallel=False,

**kwargs,

)[source]#

Type-3 NUFFT (non-uniform to non-uniform).

Parameters:

x (NDArray) – (M, [d]) d-dimensional sample points \(\mathbf{x}_{j} \in \mathbb{R}^{d}\).
z (NDArray) – (N, [d]) d-dimensional query points \(\mathbf{z}_{k} \in \mathbb{R}^{d}\).
isign (1, -1) – Sign \(\sigma\) of the transform.
eps (Real) –
Requested relative accuracy \(\varepsilon \geq 0\).

If eps=0, the transform is computed exactly via direct evaluation of the exponential sum using a Numba JIT-compiled kernel.
real (bool) –
If True, assumes apply() takes (…, M) inputs in \(\mathbb{R}^{M}\).

If False, then apply() takes (…, 2M) inputs, i.e. \(\mathbb{C}^{M}\) vectors viewed as bijections with \(\mathbb{R}^{2M}\).
plan_fw/bw (bool) – If True, allocate FINUFFT resources to do the forward (fw) and/or backward (bw) transform. These are advanced options: use them with care. Some public methods in the LinOp interface may not work if fw/bw transforms are disabled. These options only take effect if eps > 0.
enable_warnings (bool) – If True, emit a warning in case of precision mis-match issues.
chunked (bool) – If True, the transform is performed in small chunks. (See Notes for details.)
parallel (bool) – This option only applies to chunked transforms. If True, evaluate chunks in parallel.
**kwargs –
Extra kwargs to finufft.Plan. (Illegal keywords are dropped silently.) Most useful are n_trans, nthreads and debug.
plan_fw (bool) –
plan_bw (bool) –

Returns:

op – (2N, M) or (2N, 2M) type-3 NUFFT.

Return type:

Examples

import numpy as np
import pyxu.operator as pxo
import pyxu.runtime as pxrt
import pyxu.util as pxu

rng = np.random.default_rng(0)
D, M, N = 3, 200, 5  # D denotes the dimension of the data
x = rng.normal(size=(M, D)) + 2000  # Poorly-centered data
z = rng.normal(size=(N, D))
with pxrt.Precision(pxrt.Width.SINGLE):
    # The NUFFT dimension (1/2/3) is inferred from the trailing dimension of x/z.
    # Its precision is controlled by the context manager.
    N_trans = 20
    A = pxo.NUFFT.type3(
            x, z,
            n_trans=N_trans,
            isign=-1,
            eps=1e-6,
         )

    # Pyxu operators only support real inputs/outputs, so we use the functions
    # pxu.view_as_[complex/real] to interpret complex arrays as real arrays (and
    # vice-versa).
    arr =        rng.normal(size=(3, N_trans, M)) \
          + 1j * rng.normal(size=(3, N_trans, M))
    A_out_fw = pxu.view_as_complex(A.apply(pxu.view_as_real(arr)))
    A_out_bw = pxu.view_as_complex(A.adjoint(pxu.view_as_real(A_out_fw)))

Notes (chunked-transform)

An extra initialization step is required before using a chunked-transform:

A = pxl.NUFFT.type3(
        x, z
        isign
        chunked=True,   # with chunk specified
        parallel=True,  # for extra speed (chunked-only)
        **finufft_kwargs,
     )
x_chunks, z_chunks = A.auto_chunk()  # auto-determine a good x/z chunking strategy
A.allocate(x_chunks, z_chunks)  # apply the chunking strategy.

auto_chunk() is a helper method to auto-determine a good chunking strategy.

Its runtime is significant when the number of sub-problems grows large. (1000+) In these contexts, assuming a good-enough x/z-split is known in advance, users can directly supply them to allocate().

apply() / adjoint() runtime is minimized when x/z are well-ordered, i.e. when sub-problems can sub-sample inputs to apply() / adjoint() via slice operators.

To reduce runtime of chunked transforms, allocate() automatically re-orders x/z when appropriate.

The side-effect is the cost of a permutation before/after calls to apply() / adjoint(). This cost becomes significant when the number of non-uniform points x/z is large. (> 1e6)

To avoid paying the re-ordering cost at each transform, it is recommended to supply x/z and apply/adjoint inputs in the “correct” order from the start.

A good re-ordering is computed automatically by allocate() and can be used to initialize a new chunked-transform with better runtime properties as such:

### Initialize a chunked transform (1st try; as above)
A = pxl.NUFFT.type3(
        x, z
        isign
        chunked=True,
        parallel=True,
        **finufft_kwargs,
     )
x_chunks, z_chunks = A.auto_chunk()  # auto-determine a good x/z chunking strategy
A.allocate(x_chunks, z_chunks)  # will raise warning if bad x/z-order detected

### Now initialize a better transform (2nd try)
x_idx, x_chunks = A.order("x")  # get a good x-ordering
z_idx, z_chunks = A.order("z")  # get a good z-ordering
A = pxl.NUFFT.type3(
        x[x_idx], z[z_idx]  # re-order x/z accordingly
        ...                 # same as before
     )
A.allocate(x_chunks, z_chunks)  # skip auto-chunking and apply
                                # optimal x/z_chunks provided.

apply(arr)[source]#

Parameters:

arr (NDArray) –

Type 1 and 3:
- (…, M) input weights \(\mathbf{w} \in \mathbb{R}^{M}\) (real transform).
- (…, 2M) input weights \(\mathbf{w} \in \mathbb{C}^{M}\) viewed as a real array. (See view_as_real().)
Type 2:
- (…, N.prod()) input weights \(\mathbf{u} \in \mathbb{R}^{\mathcal{I}_{N_1, \ldots, N_d}}\) (real transform).
- (…, 2N.prod()) input weights \(\mathbf{u} \in \mathbb{C}^{\mathcal{I}_{N_1, \ldots, N_d}}\) viewed as a real array. (See view_as_real().)

Returns:

out –

Type 1:
(…, 2N.prod()) output weights \(\mathbf{u} \in \mathbb{C}^{\mathcal{I}_{N_1, \ldots, N_d}}\) viewed as a real array. (See view_as_real().)
Type 2:
(…, 2M) output weights \(\mathbf{w} \in \mathbb{C}^{M}\) viewed as a real array. (See view_as_real().)
Type 3:
(…, 2N) output weights \(\mathbf{v} \in \mathbb{C}^{N}\) viewed as a real array. (See view_as_real().)

Return type:

adjoint(arr)[source]#

Parameters:

arr (NDArray) –

Type 1:
(…, 2N.prod()) output weights \(\mathbf{u} \in \mathbb{C}^{\mathcal{I}_{N_1, \ldots, N_d}}\) viewed as a real array. (See view_as_real().)
Type 2:
(…, 2M) output weights \(\mathbf{w} \in \mathbb{C}^{M}\) viewed as a real array. (See view_as_real().)
Type 3:
(…, 2N) output weights \(\mathbf{v} \in \mathbb{C}^{N}\) viewed as a real array. (See view_as_real().)

Returns:

out –

Type 1 and 3:
- (…, M) input weights \(\mathbf{w} \in \mathbb{R}^{M}\) (real transform).
- (…, 2M) input weights \(\mathbf{w} \in \mathbb{C}^{M}\) viewed as a real array. (See view_as_real().)
Type 2:
- (…, N.prod()) input weights \(\mathbf{u} \in \mathbb{R}^{\mathcal{I}_{N_1, \ldots, N_d}}\) (real transform).
- (…, 2N.prod()) input weights \(\mathbf{u} \in \mathbb{C}^{\mathcal{I}_{N_1, \ldots, N_d}}\) viewed as a real array. (See view_as_real().)

Return type:

ascomplexarray(xp=<module 'numpy' from '/root/tools/miniconda3/envs/pyxu/lib/python3.11/site-packages/numpy/__init__.py'>, dtype=None)[source]#

Matrix representation (complex-valued) of the linear operator.

Parameters:

xp (ArrayModule) – Which array module to use to represent the output.
dtype (DType) – Optional (complex-valued) type of the array.

Returns:

A – Array representation of the operator (NUFFT type-dependant).

Type 1: (N.prod(), M)
Type 2: (M, N.prod())
Type 3: (N, M)

Return type:

mesh( xp=<module 'numpy' from '/root/tools/miniconda3/envs/pyxu/lib/python3.11/site-packages/numpy/__init__.py'>, dtype=None, scale='unit', upsampled=False, )[source]#

For type-1/2 NUFFT: compute the transform’s meshgrid \(\mathcal{I}_{N_{1} \times \cdots \times N_{d}} = \mathcal{I}_{N_{1}} \times \cdots \times \mathcal{I}_{N_{d}}\).

For type-3 NUFFT: compute the (shifted) meshgrid used for internal FFT computations.

Parameters:

xp (ArrayModule) – Which array module to use to represent the output.
dtype (DType) – Optional type of the array.
scale (str) –
Grid scale. Options are:
- Type1 and 2:
  - unit, i.e. \(\mathcal{I} = [[-N_{d}//2, \ldots, (N_{d}-1)//2 + 1))^{d}\)
  - source, i.e. \(\mathcal{I} \subset [-\pi, \pi)^{d}\)
- Type 3:
  - unit, i.e. \(\mathcal{I} = [[-N_{d}//2, \ldots, (N_{d}-1)//2 + 1))^{d}\)
  - source, i.e. \(\mathcal{I}_{\text{source}} \subset x^{c} + [-X_{d}, X_{d})^{d}\)
  - target, i.e. \(\mathcal{I}_{\text{target}} \subset z^{c} + [-Z_{d}, Z_{d})^{d}\),
  where \(x^{c}\), \(z^{c}\) denote the source/target centroids, and \(X\), \(Z\) the source/target half-widths.
upsampled (bool) – Use the upsampled meshgrid. (See [FINUFFT] for details.)

Returns:

grid – (N1, …, Nd, d) grid.

Return type:

Examples

import numpy as np
import pyxu.operator as pxo

rng = np.random.default_rng(0)
D, M, N = 1, 2, 3  # D denotes the dimension of the data
x = np.fmod(rng.normal(size=(M, D)), 2 * np.pi)
A = pxo.NUFFT.type1(
    x, N,
    isign=-1,
    eps=1e-3
)
A.mesh()  # array([[-1.],
          #        [ 0.],
          #        [ 1.]])

plot_kernel(ax=None, **kwargs)[source]#

Plot the spreading/interpolation kernel (along each dimension) on its support.

Parameters:

ax (Axes) – Axes to draw on. If None, a new axes is used.
**kwargs – Keyword arguments forwarded to matplotlib.axes.Axes.plot().

Returns:

ax

Return type:

Axes

Examples

import numpy as np
import pyxu.operator as pxo
import matplotlib.pyplot as plt

rng = np.random.default_rng(0)
D, M, N = 1, 2, 3  # D denotes the dimension of the data
x = np.fmod(rng.normal(size=(M, D)), 2 * np.pi)
A = pxo.NUFFT.type1(
    x, N,
    isign=-1,
    eps=1e-9
)
A.plot_kernel()
plt.show()

Notes

Requires Matplotlib to be installed.

params()[source]#

Compute internal parameters of the [FINUFFT] implementation.

Returns:

p – Internal parameters of the FINUFFT implementation, with fields:

upsample_factor: Real
FFT upsampling factor > 1
kernel_width: Integer
Width of the spreading/interpolation kernels (in number of samples).
kernel_beta: Real
Kernel decay parameter \(\beta > 0\).
fft_shape: (d,) [Integer]
Size of the D-dimensional FFT(s) performed internally.
dilation_factor: (d,) [Real]
Dilation factor(s) \(\gamma_{d}\). (Type-3 only)

Return type:

namedtuple()

Notes

When called from a chunked type-3 transform, params() returns parameters of the equivalent monolithic type-3 transform. The monolithic transform is seldom instantiable due to its large memory requirements. This method can hence be used to estimate the memory savings induced by chunking.

auto_chunk(max_mem=10, max_anisotropy=5)[source]#

(Only applies to chunked type-3 transforms.)

Auto-determine chunk indices per domain.

Use this function if you don’t know how to optimally ‘cut’ x/z manually.

Parameters:

max_mem (Real) – Max FFT memory (MiB) allowed per sub-block. (Default = 10 MiB)
max_anisotropy (Real) –
Max tolerated (normalized) anisotropy ratio >= 1.
- Setting close to 1 favors cubeoid-shaped partitions of x/z space.
- Setting large allows x/z-partitions to be highly-rectangular.

Returns:

x_chunks (list[NDArray[int]]) – (x_idx[0], …, x_idx[A-1]) x-coordinate chunk specifier. x_idx[k] contains indices of x which participate in the k-th NUFFT sub-problem.
z_chunks (list[NDArray[int]]) – (z_idx[0], …, z_idx[B-1]) z-coordinate chunk specifier. z_idx[k] contains indices of z which participate in the k-th NUFFT sub-problem.

Return type:

tuple[list[NDArray], list[NDArray]]

Notes

Chunks are identified by a custom hierarchical clustering method, with main steps:

Partition the NUFFT domains. Given a maximum FFT memory budget \(B>0\) and chunk anisotropy \(\alpha\geq 1\), partition the source/target domains into uniform rectangular cells. The (half) widths of the source/target cells \(h_{k}>0\) and \(\eta_{k}>0\) in each dimension \(k=\{1,\ldots d\}\) are chosen so as to:

Minimize the total number of partitions:

\[N_{c} = \underbrace{\prod_{k=1}^{d} \frac{X_k}{h_k}}_{\text{Source partition count}} \times \underbrace{\prod_{k=1}^{d} \frac{Z_k}{\eta_k}}_{\text{Target partition count}}\]

subject to:
1. \[\prod_{k=1}^{d} \eta_k h_k \leq (\pi/2\upsilon)^{d} \frac{B}{\delta \; \texttt{n_trans}},\]
  
  where
  \(\upsilon\) denotes the NUFFT’s grid upsampling factor,
  
  \(\delta\) the number of bytes occupied by a complex number,
  
  \(\texttt{n_trans}\) the number of simultaneous transforms performed.
2. \[\begin{split}\begin{align*} h_{k} & \leq X_{k}, \quad k=\{1,\ldots,d\} \\ \eta_{k} & \leq Z_{k}, \quad k=\{1,\ldots,d\} \end{align*}\end{split}\]
3. \[N_{c} \ge 1.\]
4. \[\begin{split}\begin{align*} \frac{1}{\alpha} & \leq \frac{h_{k}}{h_{j}} \frac{X_{j}}{X_{k}} \leq \alpha, \quad k \ne j, \\ \frac{1}{\alpha} & \leq \frac{\eta_{k}}{\eta_{j}} \frac{Z_{j}}{Z_{k}} \leq \alpha, \quad k \ne j. \end{align*}\end{split}\]
5. \[\begin{align*} \frac{1}{\alpha} & \leq \frac{h_{k}}{\eta_{j}} \frac{Z_{j}}{X_{k}} \leq \alpha, \quad (j, k) = \{1,\ldots,d\}^{2}. \end{align*}\]
Constraint (a) ensures type-3 NUFFTs performed in each sub-problem do not exceed the FFT memory budget. Constraints (b-c) ensure that partitions are non-degenerate/non-trivial respectively. Constraints (d-e) limit the anisotropy of the partitioning cells in each domain and across domains.
Data-independent Chunking. Chunk the data by assigning non-uniform samples to their enclosing cell in the partition. Empty partitions are dropped.
Data-dependent Chunk Fusion. Since (1) is data-independent, data chunks obtained in (2) may split clusters among adjacent chunks, which is undesirable. Clusters whose joint-spread is small enough are hence fused hierarchically.

Warning

This procedure yields a small number of memory-capped and well-separated data chunks in source/target domains. However, it may result in unbalanced chunks, with some chunks containing significantly more data-points than others. FINUFFT mitigates the unbalanced-chunk problem by spawning multiple threads to process dense clusters.

allocate(x_chunks, z_chunks, direct_eval_threshold=10000)[source]#

(Only applies to chunked type-3 transforms.)

Allocate NUFFT sub-problems based on chunk specification.

Parameters:

x_chunks (list[NDArray[int] | slice]) – (x_idx[0], …, x_idx[A-1]) x-coordinate chunk specifier. x_idx[k] contains indices of x which participate in the k-th NUFFT sub-problem.
z_chunks (list[NDArray[int] | slice]) – (z_idx[0], …, z_idx[B-1]) z-coordinate chunk specifier. z_idx[k] contains indices of z which participate in the k-th NUFFT sub-problem.
direct_eval_threshold (Integer) –
If provided: lower bound on len(x) * len(z) below which an NUFFT sub-problem is replaced with direct-evaluation (eps=0) for performance reasons.

(Defaults to 10k as per the FINUFFT guidelines.)

diagnostic_plot(domain)[source]#

(Only applies to chunked type-3 transforms.)

Plot data + tesselation structure for diagnostic purposes.

Parameters:: domain ('x', 'z') – Plot x-domain or z-domain data.
Returns:: fig – Diagnostic plot.
Return type:: Figure

Notes

This method can only be called after allocate().
This method only works for 2D/3D domains.

Examples

import numpy as np
import pyxu.operator as pxo

rng = np.random.default_rng(2)
D, M, N = 2, 500, 200
rnd_points = lambda _: rng.normal(scale=rng.uniform(0.25, 0.5, size=(D,)), size=(_, D))
rnd_offset = lambda: rng.uniform(-1, 1, size=(D,))
scale = 20
x = np.concatenate(
    [
        rnd_points(M) + rnd_offset() * scale,
        rnd_points(M) + rnd_offset() * scale,
        rnd_points(M) + rnd_offset() * scale,
        rnd_points(M) + rnd_offset() * scale,
        rnd_points(M) + rnd_offset() * scale,
    ],
    axis=0,
)
z = np.concatenate(
    [
        rnd_points(N) + rnd_offset() * scale,
        rnd_points(N) + rnd_offset() * scale,
        rnd_points(N) + rnd_offset() * scale,
        rnd_points(N) + rnd_offset() * scale,
        rnd_points(N) + rnd_offset() * scale,
    ],
    axis=0,
)

kwargs = dict(
    x=x,
    z=z,
    isign=-1,
    eps=1e-3,
)
A = pxo.NUFFT.type3(**kwargs, chunked=True)
x_chunks, z_chunks = A.auto_chunk(
    max_mem=.1,
    max_anisotropy=1,
)
A.allocate(x_chunks, z_chunks)
fig = A.diagnostic_plot('x')
fig.show()

Notes

Requires Matplotlib to be installed.

stats()[source]#

(Only applies to chunked type-3 transforms.)

Gather internal statistics about a chunked type-3 NUFFT.

Returns:

p – Statistics on the NUFFT chunks, with fields:

blk_count: Integer
Number of NUFFT chunks.
dEval_count: Integer
Number of chunks directly evaluated via the NUDFT.

Return type:

namedtuple()

Stencils & Convolutions #

class _Stencil(kernel, center)[source]#

Bases: object

Multi-dimensional JIT-compiled stencil. (Low-level function.)

This low-level class creates a gu-vectorized stencil applicable on multiple inputs simultaneously.

Create instances via factory method init().

Example

Correlate a stack of images A with a (3, 3) kernel such that:

\[B[n, m] = A[n-1, m] + A[n, m-1] + A[n, m+1] + A[n+1, m]\]

import numpy as np
from pyxu.operator.linop.stencil._stencil import _Stencil

# create the stencil
kernel = np.array([[0, 1, 0],
                   [1, 0, 1],
                   [0, 1, 0]], dtype=np.float64)
center = (1, 1)
stencil = _Stencil.init(kernel, center)

# apply it to the data
rng = np.random.default_rng()
A = rng.normal(size=(2, 3, 4, 30, 30))  # 24 images of size (30, 30)
B = np.zeros_like(A)
stencil.apply(A, B)  # (2, 3, 4, 30, 30)

Parameters:

kernel (NDArray) –
center (NDArray) –

IndexSpec#: alias of Sequence[Integral]

static init(kernel, center)[source]#

Parameters:

kernel (NDArray) –
(k_1, …, k_D) kernel coefficients.

Only float32/64 kernels are supported.
center (IndexSpec) – (D,) index of the kernel’s center.

Returns:

st – Rank-D stencil.

Return type:

_Stencil

apply(arr, out, **kwargs)[source]#

Evaluate stencil on multiple inputs.

Parameters:

arr (NDArray) – (…, N_1, …, N_D) data to process.
out (NDArray) – (…, N_1, …, N_D) array to which outputs are written.
kwargs (dict) –
Extra kwargs to configure f_jit(), the Dispatcher instance created by Numba.

Only relevant for GPU stencils, with values:
- blockspergrid: int
- threadsperblock: int
Default values are chosen if unspecified.

Returns:

out – (…, N_1, …, N_D) outputs.

Return type:

Notes

arr and out must have the same type/dtype as the kernel used during instantiation.
Index regions in out where the stencil is not fully supported are set to 0.
apply() may raise CUDA_ERROR_LAUNCH_OUT_OF_RESOURCES when the number of GPU registers required exceeds resource limits. There are 2 solutions to this problem:
1. Pass the max_registers kwarg to f_jit()’s decorator; or
2. Limit the number of threads per block.
1. must be set at compile time; it is thus left unbounded.
2. is accessible through .apply(**kwargs).

class Stencil(arg_shape, kernel, center, mode='constant', enable_warnings=True)[source]#

Bases: SquareOp

Multi-dimensional JIT-compiled stencil.

Stencils are a common computational pattern in which array elements are updated according to some fixed pattern called the stencil kernel. Notable examples include multi-dimensional convolutions, correlations and finite differences. (See Notes for a definition.)

Stencils can be evaluated efficiently on CPU/GPU architectures.

Several boundary conditions are supported. Moreover boundary conditions may differ per axis.

Implementation Notes

Numba (and its @stencil decorator) is used behind the scenes to compile efficient machine code from a stencil kernel specification. This has 2 consequences:
- Stencil instances are not arraymodule-agnostic: they will only work with NDArrays
  belonging to the same array module as kernel.
- Compiled stencils are not precision-agnostic: they will only work on NDArrays with the same dtype as kernel. A warning is emitted if inputs must be cast to the kernel dtype.
Stencil kernels can be specified in two forms: (See __init__() for details.)
- A single non-separable \(D\)-dimensional kernel \(k[i_{1},\ldots,i_{D}]\) of shape \((K_{1},\ldots,K_{D})\).
- A sequence of separable \(1\)-dimensional kernel(s) \(k_{d}[i]\) of shapes \((K_{1},),\ldots,(K_{D},)\) such that \(k[i_{1},\ldots,i_{D}] = \Pi_{d=1}^{D} k_{d}[i_{d}]\).

Mathematical Notes

Given a \(D\)-dimensional array \(x\in\mathbb{R}^{N_1\times\cdots\times N_D}\) and kernel \(k\in\mathbb{R}^{K_1\times\cdots\times K_D}\) with center \((c_1, \ldots, c_D)\), the output of the stencil operator is an array \(y\in\mathbb{R}^{N_1\times\cdots\times N_D}\) given by:

\[y[i_{1},\ldots,i_{D}] = \sum_{q_{1},\ldots,q_{D}=0}^{K_{1},\ldots,K_{D}} x[i_{1} - c_{1} + q_{1},\ldots,i_{D} - c_{D} + q_{D}] \,\cdot\, k[q_{1},\ldots,q_{D}].\]

This corresponds to a correlation with a shifted version of the kernel \(k\).

Numba stencils assume summation terms involving out-of-bound indices of \(x\) are set to zero. Stencil lifts this constraint by extending the stencil to boundary values via pre-padding and post-trimming. Concretely, any stencil operator \(S\) instantiated with Stencil can be written as the composition \(S = TS_0P\), where \((T, S_0, P)\) are trimming, stencil with zero-padding conditions, and padding operators respectively. This construct allows Stencil to handle complex boundary conditions under which \(S\) may not be a proper stencil (i.e., varying kernel) but can still be implemented efficiently via a proper stencil upon appropriate trimming/padding.

For example consider the decomposition of the following (improper) stencil operator:

>>> S = Stencil(
...      arg_shape=(5,),
...      kernel=np.r_[1, 2, -3],
...      center=(2,),
...      mode="reflect",
... )

>>> S.asarray()
[[-3   2   1   0   0]
 [ 2  -2   0   0   0]
 [ 1   2  -3   0   0]
 [ 0   1   2  -3   0]
 [ 0   0   1   2  -3]] # Improper stencil (kernel varies across rows)
=
[[0  0  1  0  0  0  0  0  0]
 [0  0  0  1  0  0  0  0  0]
 [0  0  0  0  1  0  0  0  0]
 [0  0  0  0  0  1  0  0  0]
 [0  0  0  0  0  0  1  0  0]]  # Trimming
@
[[-3   0   0   0   0   0   0   0   0]
 [ 2  -3   0   0   0   0   0   0   0]
 [ 1   2  -3   0   0   0   0   0   0]
 [ 0   1   2  -3   0   0   0   0   0]
 [ 0   0   1   2  -3   0   0   0   0]
 [ 0   0   0   1   2  -3   0   0   0]
 [ 0   0   0   0   1   2  -3   0   0]
 [ 0   0   0   0   0   1   2  -3   0]
 [ 0   0   0   0   0   0   1   2  -3]]  # Proper stencil (Toeplitz structure)
@
[[0  0  1  0  0]
 [0  1  0  0  0]
 [1  0  0  0  0]
 [0  1  0  0  0]
 [0  0  1  0  0]
 [0  0  0  1  0]
 [0  0  0  0  1]
 [0  0  0  1  0]
 [0  0  1  0  0]]  # Padding with reflect mode.

Note that the adjoint of a stencil operator may not necessarily be a stencil operator, or the associated center and boundary conditions may be hard to predict. For example, the adjoint of the stencil operator defined above is given by:

>>> S.T.asarray()
[[-3   2   1   0   0],
 [ 2  -2   2   1   0],
 [ 1   0  -3   2   1],
 [ 0   0   0  -3   2],
 [ 0   0   0   0  -3]]

which resembles a stencil with time-reversed kernel, but with weird (if not improper) boundary conditions. This can also be seen from the fact that \(S^\ast = P^\ast S_0^\ast T^\ast = P^\ast S_0^\ast P_0,\) and \(P^\ast\) is in general not a trimming operator. (See Pad.)

The same holds for gram/cogram operators. Consider indeed the following order-1 backward finite-difference operator with zero-padding:

>>> S = Stencil(
...      arg_shape=(5,),
...      kernel=np.r_[-1, 1],
...      center=(0,),
...      mode='constant',
... )

>>> S.gram().asarray()
[[ 1  -1   0   0   0]
 [-1   2  -1   0   0]
 [ 0  -1   2  -1   0]
 [ 0   0  -1   2  -1]
 [ 0   0   0  -1   2]]

We observe that the Gram differs from the order 2 centered finite-difference operator. (Reduced-order derivative on one side.)

Example

Moving average of a 1D signal

Let \(x[n]\) denote a 1D signal. The weighted moving average

\[y[n] = x[n-2] + 2 x[n-1] + 3 x[n]\]

can be viewed as the output of the 3-point stencil of kernel \(k = [1, 2, 3]\).

import numpy as np
from pyxu.operator import Stencil

x = np.arange(10)  # [0  1  2  3  4  5  6  7  8  9]

op = Stencil(
    arg_shape=x.shape,
    kernel=np.array([1, 2, 3]),
    center=(2,),  # k[2] applies on x[n]
)

y = op.apply(x)  # [0  3  8  14  20  26  32  38  44  50]

Non-seperable image filtering

Let \(x[n, m]\) denote a 2D image. The blurred image

\[y[n, m] = 2 x[n-1,m-1] + 3 x[n-1,m+1] + 4 x[n+1,m-1] + 5 x[n+1,m+1]\]

can be viewed as the output of the 9-point stencil

\[\begin{split}k = \left[ \begin{array}{ccc} 2 & 0 & 3 \\ 0 & 0 & 0 \\ 4 & 0 & 5 \end{array} \right].\end{split}\]

import numpy as np
from pyxu.operator import Stencil

x = np.arange(64).reshape(8, 8)  # square image
# [[ 0   1   2   3   4   5   6   7]
#  [ 8   9  10  11  12  13  14  15]
#  [16  17  18  19  20  21  22  23]
#  [24  25  26  27  28  29  30  31]
#  [32  33  34  35  36  37  38  39]
#  [40  41  42  43  44  45  46  47]
#  [48  49  50  51  52  53  54  55]
#  [56  57  58  59  60  61  62  63]]

op = Stencil(
    arg_shape=x.shape,
    kernel=np.array(
        [[2, 0, 3],
         [0, 0, 0],
         [4, 0, 5]]),
    center=(1, 1),  # k[1, 1] applies on x[n, m]
)

y = op.apply(x.reshape(-1)).reshape(8, 8)
# [[ 45   82   91  100  109  118  127   56 ]
#  [ 88  160  174  188  202  216  230  100 ]
#  [152  272  286  300  314  328  342  148 ]
#  [216  384  398  412  426  440  454  196 ]
#  [280  496  510  524  538  552  566  244 ]
#  [344  608  622  636  650  664  678  292 ]
#  [408  720  734  748  762  776  790  340 ]
#  [147  246  251  256  261  266  271  108 ]]

Seperable image filtering

Let \(x[n, m]\) denote a 2D image. The warped image

\[\begin{split}\begin{align*} y[n, m] = & + 4 x[n-1,m-1] + 5 x[n-1,m] + 6 x[n-1,m+1] \\ & + 8 x[n ,m-1] + 10 x[n ,m] + 12 x[n ,m+1] \\ & + 12 x[n+1,m-1] + 15 x[n+1,m] + 18 x[n+1,m+1] \end{align*}\end{split}\]

can be viewed as the output of the 9-point stencil

\[\begin{split}k_{2D} = \left[ \begin{array}{ccc} 4 & 5 & 6 \\ 8 & 10 & 12 \\ 12 & 15 & 18 \\ \end{array} \right].\end{split}\]

Notice however that \(y[n, m]\) can be implemented more efficiently by factoring the 9-point stencil as a cascade of two 3-point stencils:

\[\begin{split}k_{2D} = k_{1} k_{2}^{T} = \left[ \begin{array}{c} 1 \\ 2 \\ 3 \end{array} \right] \left[ \begin{array}{c} 4 & 5 & 6 \end{array} \right].\end{split}\]

Seperable stencils are supported and should be preferred when applicable.

import numpy as np
from pyxu.operator import Stencil

x = np.arange(64).reshape(8, 8)  # square image
# [[ 0   1   2   3   4   5   6   7]
#  [ 8   9  10  11  12  13  14  15]
#  [16  17  18  19  20  21  22  23]
#  [24  25  26  27  28  29  30  31]
#  [32  33  34  35  36  37  38  39]
#  [40  41  42  43  44  45  46  47]
#  [48  49  50  51  52  53  54  55]
#  [56  57  58  59  60  61  62  63]]

op_2D = Stencil(  # using non-seperable kernel
    arg_shape=x.shape,
    kernel=np.array(
        [[ 4,  5,  6],
         [ 8, 10, 12],
         [12, 15, 18]]),
    center=(1, 1),  # k[1, 1] applies on x[n, m]
)
op_sep = Stencil(  # using seperable kernels
    arg_shape=x.shape,
    kernel=[
        np.array([1, 2, 3]),  # k1: stencil along 1st axis
        np.array([4, 5, 6]),  # k2: stencil along 2nd axis
    ],
    center=(1, 1),  # k1[1] * k2[1] applies on x[n, m]
)

y_2D = op_2D.apply(x.reshape(-1)).reshape(8, 8)
y_sep = op_sep.apply(x.reshape(-1)).reshape(8, 8)  # np.allclose(y_2D, y_sep) -> True
# [[ 294   445   520   595   670   745   820   511 ]
#  [ 740  1062  1152  1242  1332  1422  1512   930 ]
#  [1268  1782  1872  1962  2052  2142  2232  1362 ]
#  [1796  2502  2592  2682  2772  2862  2952  1794 ]
#  [2324  3222  3312  3402  3492  3582  3672  2226 ]
#  [2852  3942  4032  4122  4212  4302  4392  2658 ]
#  [3380  4662  4752  4842  4932  5022  5112  3090 ]
#  [1778  2451  2496  2541  2586  2631  2676  1617 ]]

See also

Convolve, _Stencil

Parameters:

arg_shape (Integral | tuple[Integral, ...]) –
kernel (NDArray | Sequence[NDArray]) –
center (Sequence[Integral]) –
mode (str | Sequence[str]) –
enable_warnings (bool) –

KernelSpec#: alias of Union[NDArray, Sequence[NDArray]]

__init__(arg_shape, kernel, center, mode='constant', enable_warnings=True)[source]#

Parameters:

arg_shape (NDArrayShape) – Shape of the rank-\(D\) input array.
kernel (KernelSpec) –
Stencil coefficients. Two forms are accepted:
- NDArray of rank-\(D\): denotes a non-seperable stencil.
- tuple[NDArray_1, …, NDArray_D]: a sequence of 1D stencils such that dimension[k] is filtered by stencil kernel[k], that is:
  
  \[k = k_1\otimes \cdots\otimes k_D,\]
  
  or in Python: k = functools.reduce(numpy.multiply.outer, kernel).
center (IndexSpec) –
(i_1, …, i_D) index of the stencil’s center.

center defines how a kernel is overlaid on inputs to produce outputs.
mode (str, list ( str )) –
Boundary conditions. Multiple forms are accepted:
- str: unique mode shared amongst dimensions. Must be one of:
  - ’constant’ (zero-padding)
  - ’wrap’
  - ’reflect’
  - ’symmetric’
  - ’edge’
- tuple[str, …]: dimension[k] uses mode[k] as boundary condition.
(See numpy.pad() for details.)
enable_warnings (bool) – If True, emit a warning in case of precision mis-match issues.

configure_dispatcher(**kwargs)[source]#

(Only applies if kernel is a CuPy array.)

Configure stencil Dispatcher.

See apply() for accepted options.

Example

import cupy as cp
from pyxu.operator import Stencil

x = cp.arange(10)

op = Stencil(
    arg_shape=x.shape,
    kernel=np.array([1, 2, 3]),
    center=(1,),
)

y = op.apply(x)  # uses default threadsperblock/blockspergrid values

op.configure_dispatcher(
    threadsperblock=50,
    blockspergrid=3,
)
y = op.apply(x)  # supplied values used instead

property kernel: NDArray | Sequence[NDArray]#

Stencil kernel coefficients.

Returns:

kern – Stencil coefficients.

If the kernel is non-seperable, a single array is returned. Otherwise \(D\) arrays are returned, one per axis.

Return type:

KernelSpec

property center: Sequence[Integral]#

Stencil central position.

Returns:: ctr – Stencil central position.
Return type:: IndexSpec

property relative_indices: Sequence[NDArray]#

Relative indices of the stencil.

Returns:: r_idx – Relative indices of the stencil’s kernel in each dimension.
Return type:: list ( NDArray )

Example

S = Stencil(
    arg_shape=(5, 6, 9),
    kernel=[
        np.r_[1, -1],
        np.r_[3, 2, 1],
        np.r_[2, -1, 3, 1],
    ],
    center=(1, 0, 3),
)
S.relative_indices  # [array([-1,  0]), array([0, 1, 2]), array([-3, -2, -1,  0])]

visualize()[source]#

Show the \(D\)-dimensional stencil kernel.

The stencil’s center is identified by surrounding parentheses.

Example

from pyxu.operator import Stencil

S = Stencil(
    arg_shape=(5, 6),
    kernel=[
        np.r_[3, 2, 1],
        np.r_[2, -1, 3, 1],
    ],
    center=(1, 2),
)
print(S.visualize())  # [[6.0 -3.0  9.0  3.0]
                      #  [4.0 -2.0 (6.0) 2.0]
                      #  [2.0 -1.0  3.0 1.0]]

Return type:: str

Correlate#: alias of Stencil

Convolve(arg_shape, kernel, center, mode='constant', enable_warnings=True)[source]#

Multi-dimensional JIT-compiled convolution.

Inputs are convolved with the given kernel.

Notes

Given a \(D\)-dimensional array \(x\in\mathbb{R}^{N_1\times\cdots\times N_D}\) and kernel \(k\in\mathbb{R}^{K_1\times\cdots\times K_D}\) with center \((c_1, \ldots, c_D)\), the output of the convolution operator is an array \(y\in\mathbb{R}^{N_1\times\cdots\times N_D}\) given by:

\[y[i_{1},\ldots,i_{D}] = \sum_{q_{1},\ldots,q_{D}=0}^{K_{1},\ldots,K_{D}} x[i_{1} - q_{1} + c_{1},\ldots,i_{D} - q_{D} + c_{D}] \,\cdot\, k[q_{1},\ldots,q_{D}].\]

The convolution is implemented via Stencil. To do so, the convolution kernel is transformed to the equivalent correlation kernel:

\[\begin{split}y[i_{1},\ldots,i_{D}] = \sum_{q_{1},\ldots,q_{D}=0}^{K_{1},\ldots,K_{D}} &x[i_{1} + q_{1} - (K_{1} - c_{1}),\ldots,i_{D} + q_{D} - (K_{D} - c_{D})] \\ &\cdot\, k[K_{1}-q_{1},\ldots,K_{D}-q_{D}].\end{split}\]

This corresponds to a correlation with a flipped kernel and center.

Examples

import numpy as np
from scipy.ndimage import convolve
from pyxu.operator import Convolve

x = np.array([
     [1, 2, 0, 0],
     [5, 3, 0, 4],
     [0, 0, 0, 7],
     [9, 3, 0, 0],
])
k = np.array([
     [1, 1, 1],
     [1, 1, 0],
     [1, 0, 0],
])
op = Convolve(
    arg_shape=x.shape,
    kernel=k,
    center=(1, 1),
    mode="constant",
)

y_op = op.apply(x.reshape(-1)).reshape(4, 4)
y_sp = convolve(x, k, mode="constant", origin=0)  # np.allclose(y_op, y_sp) -> True
# [[11  10   7   4],
#  [10   3  11  11],
#  [15  12  14   7],
#  [12   3   7   0]]

See also

Stencil

Parameters:

arg_shape (Integral | tuple[Integral, ...]) –
kernel (NDArray | Sequence[NDArray]) –
center (Sequence[Integral]) –
mode (str | Sequence[str]) –
enable_warnings (bool) –

Filters #

MovingAverage(arg_shape, size, center=None, mode='constant', gpu=False, dtype=None)[source]#

Multi-dimensional moving average or uniform filter.

Notes

This operator performs a convolution between the input \(D\)-dimensional NDArray \(\mathbf{x} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}\) and a uniform \(D\)-dimensional filter \(\mathbf{h} \in \mathbb{R}^{\text{size} \times \cdots \times \text{size}}\) that computes the size-point local mean values using separable kernels for improved performance.

\[y_{i} = \frac{1}{|\mathcal{N}_{i}|}\sum_{j \in \mathcal{N}_{i}} x_{j}\]

Where \(\mathcal{N}_{i}\) is the set of elements neighbouring the \(i\)-th element of the input array, and \(\mathcal{N}_{i}\) denotes its cardinality, i.e. the total number of neighbors.

Parameters:

arg_shape (NDArrayShape) – Shape of the input array.
size (int, tuple) –
Size of the moving average kernel.

If a single integer value is provided, then the moving average filter will have as many dimensions as the input array. If a tuple is provided, it should contain as many elements as arg_shape. For example, the size=(1, 3) will convolve the input image with the filter [[1, 1, 1]] / 3.
center (IndexSpec) –
(i_1, …, i_D) index of the kernel’s center.

center defines how a kernel is overlaid on inputs to produce outputs. For odd size, it defaults to the central element (center=size//2). For even size the desired center indices must be provided.
mode (str, list[str]) –
Boundary conditions. Multiple forms are accepted:
- str: unique mode shared amongst dimensions. Must be one of:
  - ’constant’ (default): zero-padding
  - ’wrap’
  - ’reflect’
  - ’symmetric’
  - ’edge’
- tuple[str, …]: the d-th dimension uses mode[d] as boundary condition.
(See numpy.pad() for details.)
gpu (bool) – Input NDArray type (True for GPU, False for CPU). Defaults to False.
dtype (DType) – Working precision of the linear operator.

Returns:

op – MovingAverage

Return type:

Example

import matplotlib.pyplot as plt
from pyxu.operator import MovingAverage

arg_shape = (11, 11)
image = np.zeros(arg_shape)
image[5, 5] = 1.

ma = MovingAverage(arg_shape, size=5)
out = ma(image.ravel())
plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.imshow(image)
plt.colorbar()
plt.subplot(122)
plt.imshow(out.reshape(*arg_shape))
plt.colorbar()

See also

Gaussian

Gaussian(arg_shape, sigma=1.0, truncate=3.0, order=0, mode='constant', sampling=1, gpu=False, dtype=None)[source]#

Multi-dimensional Gaussian filter.

Notes

This operator performs a convolution between the input \(D\)-dimensional NDArray \(\mathbf{x} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}\) and a Gaussian \(D\)-dimensional filter \(\mathbf{h} \in \mathbb{R}^{\text{size} \times \cdots \times \text{size}}\) using separable kernels for improved performance.

\[y_{i} = \sum_{j \in \mathcal{N}_{i}} a_{j} x_{j} \exp(\frac{d_{ij}^{2}}{\sigma^{2}})\]

Where \(\mathcal{N}_{i}\) is the set of elements neighbouring the \(i\)-th element of the input array, and \(a_{j} = \sum_{j \in \mathcal{N}_{i}} a_{j} \exp(\frac{d_{ij}^{2}}{\sigma^{2}})\) normalizes the kernel to sum to one.

Parameters:

arg_shape (NDArrayShape) – Shape of the input array.
sigma (float, tuple) –
Standard deviation of the Gaussian kernel.

If a scalar value is provided, then the Gaussian filter will have as many dimensions as the input array. If a tuple is provided, it should contain as many elements as arg_shape. Use 0 to prevent filtering in a given dimension. For example, the sigma=(0, 3) will convolve the input image in its last dimension.
truncate (float, tuple) – Truncate the filter at this many standard deviations. Defaults to 3.0.
order (int, tuple) –
Gaussian derivative order.

Use 0 for the standard Gaussian kernel.
mode (str, list[str]) –
Boundary conditions. Multiple forms are accepted:
- str: unique mode shared amongst dimensions. Must be one of:
  - ’constant’ (default): zero-padding
  - ’wrap’
  - ’reflect’
  - ’symmetric’
  - ’edge’
- tuple[str, …]: the d-th dimension uses mode[d] as boundary condition.
(See numpy.pad() for details.)
sampling (Real, list[Real]) – Sampling step (i.e. distance between two consecutive elements of an array). Defaults to 1.
gpu (bool) – Input NDArray type (True for GPU, False for CPU). Defaults to False.
dtype (DType) – Working precision of the linear operator.

Returns:

op – Gaussian

Return type:

Example

import matplotlib.pyplot as plt
from pyxu.operator import Gaussian

arg_shape = (11, 11)
image = np.zeros(arg_shape)
image[5, 5] = 1.

gaussian = Gaussian(arg_shape, sigma=3)
out = gaussian(image.ravel())
plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.imshow(image)
plt.colorbar()
plt.subplot(122)
plt.imshow(out.reshape(*arg_shape))
plt.colorbar()

DifferenceOfGaussians( arg_shape, low_sigma=1.0, high_sigma=None, low_truncate=3.0, high_truncate=3.0, mode='constant', sampling=1, gpu=False, dtype=None, )[source]#

Multi-dimensional Difference of Gaussians filter.

Notes

This operator uses the Difference of Gaussians (DoG) method to a \(D\)-dimensional NDArray \(\mathbf{x} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}\) using separable kernels for improved performance. The DoG method blurs the input image with two Gaussian kernels with different sigma, and subtracts the more-blurred signal from the less-blurred image. This creates an output signal containing only the information from the original signal at the spatial scale indicated by the two sigmas.

Parameters:

arg_shape (NDArrayShape) – Shape of the input array.
low_sigma (float, tuple) –
Standard deviation of the Gaussian kernel with smaller sigmas across all axes.

If a scalar value is provided, then the Gaussian filter will have as many dimensions as the input array. If a tuple is provided, it should contain as many elements as arg_shape. Use 0 to prevent filtering in a given dimension. For example, the low_sigma=(0, 3) will convolve the input image in its last dimension.
high_sigma (float, tuple, None) – Standard deviation of the Gaussian kernel with larger sigmas across all axes. If None is given (default), sigmas for all axes are calculated as 1.6 * low_sigma.
low_truncate (float, tuple) – Truncate the filter at this many standard deviations. Defaults to 3.0.
high_truncate (float, tuple) – Truncate the filter at this many standard deviations. Defaults to 3.0.
mode (str, list[str]) –
Boundary conditions. Multiple forms are accepted:
- str: unique mode shared amongst dimensions. Must be one of:
  - ’constant’ (default): zero-padding
  - ’wrap’
  - ’reflect’
  - ’symmetric’
  - ’edge’
- tuple[str, …]: the d-th dimension uses mode[d] as boundary condition.
(See numpy.pad() for details.)
sampling (Real, list[Real]) – Sampling step (i.e. distance between two consecutive elements of an array). Defaults to 1.
gpu (bool) – Input NDArray type (True for GPU, False for CPU). Defaults to False.
dtype (DType) – Working precision of the linear operator.

Returns:

op – DifferenceOfGaussians

Return type:

Example

import matplotlib.pyplot as plt
from pyxu.operator import DoG

arg_shape = (11, 11)
image = np.zeros(arg_shape)
image[5, 5] = 1.

dog = DoG(arg_shape, low_sigma=3)
out = dog(image.ravel())
plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.imshow(image)
plt.colorbar()
plt.subplot(122)
plt.imshow(out.reshape(*arg_shape))
plt.colorbar()

DoG(arg_shape, low_sigma=1.0, high_sigma=None, low_truncate=3.0, high_truncate=3.0, mode='constant', sampling=1, gpu=False, dtype=None)#

Alias of DifferenceOfGaussians().

Parameters:

arg_shape (Integral | tuple[Integral, ...]) –
mode (str | Sequence[str]) –
sampling (Real | Sequence[Real, ...]) –
gpu (bool) –
dtype (dtype[Any] | None | type[Any] | _SupportsDType[dtype[Any]] | str | tuple[Any, int] | tuple[Any, SupportsIndex | Sequence[SupportsIndex]] | list[Any] | _DTypeDict | tuple[Any, Any]) –

Laplace(arg_shape, mode='constant', sampling=1, gpu=False, dtype=None)[source]#

Multi-dimensional Laplace filter.

The implementation is based on second derivatives approximated via finite differences.

Notes

This operator uses the applies the Laplace kernel \([1 -2 1]\) to a \(D\)-dimensional NDArray \(\mathbf{x} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}\) using separable kernels for improved performance. The Laplace filter is commonly used to find high-frequency components in the signal, such as for example, the edges in an image.

Parameters:

arg_shape (NDArrayShape) – Shape of the input array.
mode (str, list[str]) –
Boundary conditions. Multiple forms are accepted:
- str: unique mode shared amongst dimensions. Must be one of:
  - ’constant’ (default): zero-padding
  - ’wrap’
  - ’reflect’
  - ’symmetric’
  - ’edge’
- tuple[str, …]: the d-th dimension uses mode[d] as boundary condition.
(See numpy.pad() for details.)
sampling (Real, list[Real]) – Sampling step (i.e. distance between two consecutive elements of an array). Defaults to 1.
gpu (bool) – Input NDArray type (True for GPU, False for CPU). Defaults to False.
dtype (DType) – Working precision of the linear operator.

Returns:

op – DifferenceOfGaussians

Return type:

Example

import matplotlib.pyplot as plt
from pyxu.operator import Laplace

arg_shape = (11, 11)
image = np.zeros(arg_shape)
image[5, 5] = 1.

laplace = Laplace(arg_shape)
out = laplace(image.ravel())
plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.imshow(image)
plt.colorbar()
plt.subplot(122)
plt.imshow(out.reshape(*arg_shape))
plt.colorbar()

See also

Sobel, Prewitt, Scharr

Sobel(arg_shape, axis=None, mode='constant', sampling=1, gpu=False, dtype=None)[source]#

Multi-dimensional Sobel filter.

Notes

This operator uses the applies the multi-dimensional Sobel filter to a \(D\)-dimensional NDArray \(\mathbf{x} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}\) using separable kernels for improved performance. The Sobel filter applies the following edge filter in the dimensions of interest: [1, 0, -1] and the smoothing filter on the rest of dimensions: [1, 2, 1] / 4. The Sobel filter is commonly used to find high-frequency components in the signal, such as for example, the edges in an image.

Parameters:

arg_shape (NDArrayShape) – Shape of the input array.
axis (int, tuple) –
Compute the edge filter along this axis. If not provided, the edge magnitude is computed.

This is defined as: np.sqrt(sum([sobel(array, axis=i)**2 for i in range(array.ndim)]) / array.ndim) The magnitude is also computed if axis is a sequence.
mode (str, list[str]) –
Boundary conditions. Multiple forms are accepted:
- str: unique mode shared amongst dimensions. Must be one of:
  - ’constant’ (default): zero-padding
  - ’wrap’
  - ’reflect’
  - ’symmetric’
  - ’edge’
- tuple[str, …]: the d-th dimension uses mode[d] as boundary condition.
(See numpy.pad() for details.)
sampling (Real, list[Real]) – Sampling step (i.e. distance between two consecutive elements of an array). Defaults to 1.
gpu (bool) – Input NDArray type (True for GPU, False for CPU). Defaults to False.
dtype (DType) – Working precision of the linear operator.

Returns:

op – Sobel

Return type:

Example

import matplotlib.pyplot as plt
from pyxu.operator import Sobel

arg_shape = (11, 11)
image = np.zeros(arg_shape)
image[5, 5] = 1.

sobel = Sobel(arg_shape)
out = sobel(image.ravel())
plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.imshow(image)
plt.colorbar()
plt.subplot(122)
plt.imshow(out.reshape(*arg_shape))
plt.colorbar()

See also

Prewitt, Scharr

Prewitt(arg_shape, axis=None, mode='constant', sampling=1, gpu=False, dtype=None)[source]#

Multi-dimensional Prewitt filter.

Notes

This operator uses the applies the multi-dimensional Prewitt filter to a \(D\)-dimensional NDArray \(\mathbf{x} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}\) using separable kernels for improved performance. The Prewitt filter applies the following edge filter in the dimensions of interest: [1, 0, -1], and the smoothing filter on the rest of dimensions: [1, 1, 1] / 3. The Prewitt filter is commonly used to find high-frequency components in the signal, such as for example, the edges in an image.

Parameters:

arg_shape (NDArrayShape) – Shape of the input array.
axis (int, tuple) –
Compute the edge filter along this axis. If not provided, the edge magnitude is computed. This is defined as:

np.sqrt(sum([prewitt(array, axis=i)**2 for i in range(array.ndim)]) / array.ndim) The magnitude is also computed if axis is a sequence.
mode (str, list[str]) –
Boundary conditions. Multiple forms are accepted:
- str: unique mode shared amongst dimensions. Must be one of:
  - ’constant’ (default): zero-padding
  - ’wrap’
  - ’reflect’
  - ’symmetric’
  - ’edge’
- tuple[str, …]: the d-th dimension uses mode[d] as boundary condition.
(See numpy.pad() for details.)
sampling (Real, list[Real]) – Sampling step (i.e. distance between two consecutive elements of an array). Defaults to 1.
gpu (bool) – Input NDArray type (True for GPU, False for CPU). Defaults to False.
dtype (DType) – Working precision of the linear operator.

Returns:

op – Prewitt

Return type:

Example

import matplotlib.pyplot as plt
from pyxu.operator import Prewitt

arg_shape = (11, 11)
image = np.zeros(arg_shape)
image[5, 5] = 1.

prewitt = Prewitt(arg_shape)
out = prewitt(image.ravel())
plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.imshow(image)
plt.colorbar()
plt.subplot(122)
plt.imshow(out.reshape(*arg_shape))
plt.colorbar()

See also

Sobel, Scharr

Scharr(arg_shape, axis=None, mode='constant', sampling=1, gpu=False, dtype=None)[source]#

Multi-dimensional Scharr filter.

Notes

This operator uses the applies the multi-dimensional Scharr filter to a \(D\)-dimensional NDArray \(\mathbf{x} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}\) using separable kernels for improved performance. The Scharr filter applies the following edge filter in the dimensions of interest: [1, 0, -1], and the smoothing filter on the rest of dimensions: [3, 10, 3] / 16. The Scharr filter is commonly used to find high-frequency components in the signal, such as for example, the edges in an image.

Parameters:

arg_shape (NDArrayShape) – Shape of the input array.
axis (int, tuple) –
Compute the edge filter along this axis. If not provided, the edge magnitude is computed. This is defined as:

np.sqrt(sum([scharr(array, axis=i)**2 for i in range(array.ndim)]) / array.ndim) The magnitude is also computed if axis is a sequence.
mode (str, list[str]) –
Boundary conditions. Multiple forms are accepted:
- str: unique mode shared amongst dimensions. Must be one of:
  - ’constant’ (default): zero-padding
  - ’wrap’
  - ’reflect’
  - ’symmetric’
  - ’edge’
- tuple[str, …]: the d-th dimension uses mode[d] as boundary condition.
(See numpy.pad() for details.)
sampling (Real, list[Real]) – Sampling step (i.e. distance between two consecutive elements of an array). Defaults to 1.
gpu (bool) – Input NDArray type (True for GPU, False for CPU). Defaults to False.
dtype (DType) – Working precision of the linear operator.

Returns:

op – Scharr

Return type:

Example

import matplotlib.pyplot as plt
from pyxu.operator import Scharr

arg_shape = (11, 11)
image = np.zeros(arg_shape)
image[5, 5] = 1.

scharr = Scharr(arg_shape)
out = scharr(image.ravel())
plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.imshow(image)
plt.colorbar()
plt.subplot(122)
plt.imshow(out.reshape(*arg_shape))
plt.colorbar()

See also

Sobel, Prewitt

class StructureTensor(

arg_shape,

diff_method='fd',

smooth_sigma=1.0,

smooth_truncate=3.0,

mode='constant',

sampling=1,

gpu=False,

dtype=None,

parallel=False,

**diff_kwargs,

)[source]#

Bases: DiffMap

Structure tensor operator.

Notes

The Structure Tensor, also known as the second-order moment tensor or the inertia tensor, is a matrix derived from the gradient of a function. It describes the distribution of the gradient (i.e., its prominent directions) in a specified neighbourhood around a point, and the degree to which those directions are coherent. The structure tensor of a \(D\)-dimensional signal \(\mathbf{f} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}\) can be written as:

\[\begin{split}\mathbf{S}_\sigma \mathbf{f} = \mathbf{g}_{\sigma} * \nabla\mathbf{f} (\nabla\mathbf{f})^{\top} = \mathbf{g}_{\sigma} * \begin{bmatrix} \left( \dfrac{ \partial\mathbf{f} }{ \partial x_{0} } \right)^2 & \dfrac{ \partial^{2}\mathbf{f} }{ \partial x_{0}\,\partial x_{1} } & \cdots & \dfrac{ \partial\mathbf{f} }{ \partial x_{0} } \dfrac{ \partial\mathbf{f} }{ \partial x_{D-1} } \\ \dfrac{ \partial\mathbf{f} }{ \partial x_{1} } \dfrac{ \partial\mathbf{f} }{ \partial x_{0} } & \left( \dfrac{ \partial\mathbf{f} }{ \partial x_{1} }\right)^2 & \cdots & \dfrac{ \partial\mathbf{f} }{ \partial x_{1} } \dfrac{ \partial\mathbf{f} }{ \partial x_{D-1} } \\ \vdots & \vdots & \ddots & \vdots \\ \dfrac{ \partial\mathbf{f} }{ \partial x_{D-1} } \dfrac{ \partial\mathbf{f} }{ \partial x_{0} } & \dfrac{ \partial\mathbf{f} }{ \partial x_{D-1} } \dfrac{ \partial\mathbf{f} }{ \partial x_{1} } & \cdots & \left( \dfrac{ \partial\mathbf{f} }{ \partial x_{D-1}} \right)^2 \end{bmatrix},\end{split}\]

where \(\mathbf{g}_{\sigma} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}\) is a discrete Gaussian filter with standard variation \(\sigma\) with which a convolution is performed elementwise.

However, due to the symmetry of the structure tensor, only the upper triangular part is computed in practice:

\[\begin{split}\mathbf{H}_{\mathbf{v}_1, \ldots ,\mathbf{v}_m} \mathbf{f} = \mathbf{g}_{\sigma} * \begin{bmatrix} \left( \dfrac{ \partial\mathbf{f} }{ \partial x_{0} } \right)^2 \\ \dfrac{ \partial^{2}\mathbf{f} }{ \partial x_{0}\,\partial x_{1} } \\ \vdots \\ \left( \dfrac{ \partial\mathbf{f} }{ \partial x_{D-1}} \right)^2 \end{bmatrix} \mathbf{f} \in \mathbb{R}^{\frac{D (D-1)}{2} \times N_0 \times \cdots \times N_{D-1}}\end{split}\]

Remark#

In case of using the finite differences (diff_type="fd"), the finite difference scheme defaults to central (see PartialDerivative).

Example

import numpy as np
import matplotlib.pyplot as plt
from pyxu.operator import StructureTensor
from pyxu.util.misc import peaks

# Define input image
n = 1000
x = np.linspace(-3, 3, n)
xx, yy = np.meshgrid(x, x)
image = peaks(xx, yy)
nsamples = 2
arg_shape = image.shape  # (1000, 1000)
images = np.tile(image, (nsamples, 1, 1)).reshape(nsamples, -1)
print(images.shape)  # (2, 1000000)
# Instantiate structure tensor operator
structure_tensor = StructureTensor(arg_shape=arg_shape)

outputs = structure_tensor(images)
print(outputs.shape)  # (2, 3000000)
# Plot
outputs = structure_tensor.unravel(outputs)
print(outputs.shape)  # (2, 3, 1000, 1000)
plt.figure()
plt.imshow(images[0].reshape(arg_shape))
plt.colorbar()
plt.title("Image")
plt.axis("off")

plt.figure()
plt.imshow(outputs[0][0].reshape(arg_shape))
plt.colorbar()
plt.title(r"$\hat{S}_{xx}$")
plt.axis("off")

plt.figure()
plt.imshow(outputs[0][1].reshape(arg_shape))
plt.colorbar()
plt.title(r"$\hat{S}_{xy}$")
plt.axis("off")

plt.figure()
plt.imshow(outputs[0][2].reshape(arg_shape))
plt.colorbar()
plt.title(r"$\hat{S}_{yy}$")
plt.axis("off")

../../_images/linop-10_00.png — (`png`, `hires.png`, `pdf`)#

../../_images/linop-10_01.png — (`png`, `hires.png`, `pdf`)#

../../_images/linop-10_02.png — (`png`, `hires.png`, `pdf`)#

../../_images/linop-10_03.png — (`png`, `hires.png`, `pdf`)#

unravel(arr)[source]#

ravel(arr)[source]#

apply(arr)[source]#

Evaluate operator at specified point(s).

Parameters:: arr (NDArray) – (…, M) input points.
Returns:: out – (…, N) output points.
Return type:: NDArray

Parameters:

arg_shape (Integral | tuple[Integral, ...]) –
smooth_sigma (Real | tuple[Real, ...]) –
smooth_truncate (Real | tuple[Real, ...]) –
mode (str | Sequence[str]) –
sampling (Real | tuple[Real, ...]) –
gpu (bool) –
dtype (dtype[Any] | None | type[Any] | _SupportsDType[dtype[Any]] | str | tuple[Any, int] | tuple[Any, SupportsIndex | Sequence[SupportsIndex]] | list[Any] | _DTypeDict | tuple[Any, Any]) –
parallel (bool) –

Derivatives #

class PartialDerivative[source]#

Bases: object

Partial derivative operator based on Numba stencils.

Notes

This operator approximates the partial derivative of a \(D\)-dimensional signal \(\mathbf{f} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}\)

\[\frac{\partial^{n} \mathbf{f}}{\partial x_0^{n_0} \cdots \partial x_{D-1}^{n_{D-1}}} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}\]

where \(\frac{\partial^{n_i}}{\partial x_i^{n_i}}\) is the \(n_i\)-th order partial derivative along dimension \(i\) and \(n = \prod_{i=0}^{D-1} n_{i}\) is the total derivative order.

Partial derivatives can be implemented with finite differences via the finite_difference() constructor, or with the Gaussian derivative via the gaussian_derivative() constructor.
When using the finite_difference() constructor, the adjoint of the resulting linear operator will vary depending on the type of finite differences:
- For forward type, the adjoint corresponds to:
  
  \((\frac{\partial^{\text{fwd}}}{\partial x})^{\ast} = -\frac{\partial^{\text{bwd}}}{\partial x}\)
- For backward type, the adjoint corresponds to:
  
  \((\frac{\partial^{\text{bwd}}}{\partial x})^{\ast} = -\frac{\partial^{\text{fwd}}}{\partial x}\)
- For central type, and for the gaussian_derivative() constructor, the adjoint corresponds to:
  
  \((\frac{\partial}{\partial x})^{\ast} = -\frac{\partial}{\partial x}\)

Warning

When dealing with high-order partial derivatives, the stencils required to compute them can become large, resulting in computationally expensive evaluations. In such scenarios, it can be more efficient to construct the partial derivative through a composition of lower-order partial derivatives.

See also

Gradient, Laplacian, Hessian

static finite_difference(arg_shape, order, scheme='forward', accuracy=1, mode='constant', gpu=False, dtype=None, sampling=1)[source]#

Compute partial derivatives for multi-dimensional signals using finite differences.

Parameters:

arg_shape (NDArrayShape) – Shape of the input array.
order (list[Integer]) – Derivative order for each dimension. The total order of the partial derivative is the sum of the elements of the tuple. Use zeros to indicate dimensions in which the derivative should not be computed.
scheme (str, list[str]) – Type of finite differences: [‘forward, ‘backward, ‘central’]. Defaults to ‘forward’. If a string is provided, the same scheme is assumed for all dimensions. If a tuple is provided, it should have as many elements as order.
accuracy (Integer, list[Integer]) – Determines the number of points used to approximate the derivative with finite differences (see Notes). Defaults to 1. If an int is provided, the same accuracy is assumed for all dimensions. If a tuple is provided, it should have as many elements as arg_shape.
mode (str, list[str]) –
Boundary conditions. Multiple forms are accepted:
- str: unique mode shared amongst dimensions. Must be one of:
  - ’constant’ (default): zero-padding
  - ’wrap’
  - ’reflect’
  - ’symmetric’
  - ’edge’
- tuple[str, …]: the d-th dimension uses mode[d] as boundary condition.
(See numpy.pad() for details.)
gpu (bool) – Input NDArray type (True for GPU, False for CPU). Defaults to False.
dtype (DType) – Working precision of the linear operator.
sampling (Real, list[Real]) – Sampling step (i.e. distance between two consecutive elements of an array). Defaults to 1.

Returns:

op – Partial derivative

Return type:

Notes

We explain here finite differences for one-dimensional signals; this operator performs finite differences for multi-dimensional signals along dimensions specified by order.

This operator approximates derivatives with finite differences. It is inspired by the Finite Difference Coefficients Calculator to construct finite-difference approximations for the desired (i) derivative order, (ii) approximation accuracy, and (iii) finite difference type. Three basic types of finite differences are supported, which lead to the following first-order (order = 1) operators with accuracy = 1 and sampling step sampling = h for one-dimensional signals:

Forward difference: Approximates the continuous operator \(D_{F}f(x) = \frac{f(x+h) - f(x)}{h}\) with the discrete operator

\[\mathbf{D}_{F} f [n] = \frac{f[n+1] - f[n]}{h},\]

whose kernel is \(d = \frac{1}{h}[-1, 1]\) and center is (0, ).
Backward difference: Approximates the continuous operator \(D_{B}f(x) = \frac{f(x) - f(x-h)}{h}\) with the discrete operator

\[\mathbf{D}_{B} f [n] = \frac{f[n] - f[n-1]}{h},\]

whose kernel is \(d = \frac{1}{h}[-1, 1]\) and center is (1, ).
Central difference: Approximates the continuous operator \(D_{C}f(x) = \frac{f(x+h) - f(x-h)}{2h}\) with the discrete operator

\[\mathbf{D}_{C} f [n] = \frac{f[n+1] - f[n-1]}{2h},\]

whose kernel is \(d = \frac{1}{h}[-\frac12, 0, \frac12]\) and center is (1, ).

Warning

For forward and backward differences, higher-order operators correspond to the composition of first-order operators. This is not the case for central differences: the second-order continuous operator is given by \(D^2_{C}f(x) = \frac{f(x+h) - 2 f(x) + f(x-h)}{h}\), hence \(D^2_{C} \neq D_{C} \circ D_{C}\). The corresponding discrete operator is given by \(\mathbf{D}^2_{C} f [n] = \frac{f[n+1] - 2 f[n] + f[n-1]}{h}\), whose kernel is \(d = \frac{1}{h}[1, -2, 1]\) and center is (1, ). We refer to this paper for more details.

For a given derivative order \(N\in\mathbb{Z}^{+}\) and accuracy \(a\in\mathbb{Z}^{+}\), the size \(N_s\) of the stencil kernel \(d\) used for finite differences is given by:

For central differences: \(N_s = 2 \lfloor\frac{N + 1}{2}\rfloor - 1 + a\)
For forward and backward differences: \(N_s = N + a\)

For \(N_s\) given support indices \(\{s_1, \ldots , s_{N_s} \} \subset \mathbb{Z}\) and a derivative order \(N<N_s\), the stencil kernel \(d = [d_1, \ldots, d_{N_s}]\) of the finite-difference approximation of the derivative is obtained by solving the following system of linear equations (see the Finite Difference Coefficients Calculator documentation):

Remark

The number of coefficients of the finite-difference kernel is chosen to guarantee the requested accuracy, but might be larger than requested accuracy. For example, if choosing scheme='central' with accuracy=1, it will create a kernel corresponding to accuracy=2, as it is the minimum accuracy possible for such scheme (see the finite difference coefficient table).

\[\begin{split}\left(\begin{array}{ccc} s_{1}^{0} & \cdots & s_{N_s}^{0} \\ \vdots & \ddots & \vdots \\ s_{1}^{N_s-1} & \cdots & s_{N_s}^{N_s-1} \end{array}\right)\left(\begin{array}{c} d_{1} \\ \vdots \\ d_{N_s} \end{array}\right)= \frac{1}{h^{N}}\left(\begin{array}{c} \delta_{0, N} \\ \vdots \\ \delta_{i, N} \\ \vdots \\ \delta_{N_s-1, N} \end{array}\right),\end{split}\]

where \(\delta_{i, j}\) is the Kronecker delta.

This class inherits its methods from Stencil.

Example

import numpy as np
import matplotlib.pyplot as plt
from pyxu.operator import PartialDerivative
from pyxu.util.misc import peaks

x = np.linspace(-2.5, 2.5, 25)
xx, yy = np.meshgrid(x, x)
image = peaks(xx, yy)
arg_shape = image.shape  # Shape of our image
# Specify derivative order at each direction
df_dx = (1, 0)  # Compute derivative of order 1 in first dimension
d2f_dy2 = (0, 2)  # Compute derivative of order 2 in second dimension
d3f_dxdy2 = (1, 2)  # Compute derivative of order 1 in first dimension and der. of order 2 in second dimension
# Instantiate derivative operators
sigma = 2.0
diff1 = PartialDerivative.gaussian_derivative(order=df_dx, arg_shape=arg_shape, sigma=sigma / np.sqrt(2))
diff2 = PartialDerivative.gaussian_derivative(order=d2f_dy2, arg_shape=arg_shape, sigma=sigma / np.sqrt(2))
diff = PartialDerivative.gaussian_derivative(order=d3f_dxdy2, arg_shape=arg_shape, sigma=sigma)
# Compute derivatives
out1 = (diff1 * diff2)(image.ravel()).reshape(arg_shape)
out2 = diff(image.ravel()).reshape(arg_shape)
# Plot derivatives
fig, axs = plt.subplots(1, 3, figsize=(15, 4))
im = axs[0].imshow(image)
axs[0].axis("off")
axs[0].set_title("f(x,y)")
plt.colorbar(im, ax=axs[0])
axs[1].imshow(out1)
axs[1].axis("off")
axs[1].set_title(r"$\frac{\partial^{3} f(x,y)}{\partial x\partial y^{2}}$")
plt.colorbar(im, ax=axs[1])

axs[2].imshow(out2)
axs[2].axis("off")
axs[2].set_title(r"$\frac{\partial^{3} f(x,y)}{\partial x\partial y^{2}}$")
plt.colorbar(im, ax=axs[2])

# Check approximation error
plt.figure()
plt.imshow(abs(out1 - out2)), plt.colorbar()

../../_images/linop-11_00.png — (`png`, `hires.png`, `pdf`)#

../../_images/linop-11_01.png — (`png`, `hires.png`, `pdf`)#

static gaussian_derivative(arg_shape, order, sigma=1.0, truncate=3.0, mode='constant', gpu=False, dtype=None, sampling=1)[source]#

Compute partial derivatives for multi-dimensional signals using gaussian derivatives.

Parameters:

arg_shape (NDArrayShape) – Shape of the input array.
order (list[Integer]) – Derivative order for each dimension. The total order of the partial derivative is the sum of the elements of the tuple. Use zeros to indicate dimensions in which the derivative should not be computed.
sigma (Real, list[Real]) – Standard deviation for the Gaussian kernel. Defaults to 1.0. If a float is provided, the same sigma is assumed for all dimensions. If a tuple is provided, it should have as many elements as order.
truncate (Real, list[Real]) – Truncate the filter at this many standard deviations (at each side from the origin). Defaults to 3.0. If a float is provided, the same truncate is assumed for all dimensions. If a tuple is provided, it should have as many elements as order.
mode (str, list[str]) –
Boundary conditions. Multiple forms are accepted:
- str: unique mode shared amongst dimensions. Must be one of:
  - ’constant’ (default): zero-padding
  - ’wrap’
  - ’reflect’
  - ’symmetric’
  - ’edge’
- tuple[str, …]: the d-th dimension uses mode[d] as boundary condition.
(See numpy.pad() for details.)
gpu (bool) – Input NDArray type (True for GPU, False for CPU). Defaults to False.
dtype (DType) – Working precision of the linear operator.
sampling (Real, list[Real]) – Sampling step (i.e., the distance between two consecutive elements of an array). Defaults to 1.

Returns:

op – Partial derivative.

Return type:

Notes

We explain here Gaussian derivatives for one-dimensional signals; this operator performs partial Gaussian derivatives for multi-dimensional signals along dimensions specified by order.

A Gaussian derivative is an approximation of a derivative that consists in convolving the input function with a Gaussian function \(g\) before applying a derivative. In the continuous domain, the \(N\)-th order Gaussian derivative \(D^N_G\) amounts to a convolution with the \(N\)-th order derivative of \(g\):

\[D^N_G f (x) = \frac{\mathrm{d}^N (f * g) }{\mathrm{d} x^N} (x) = f(x) * \frac{\mathrm{d}^N g}{\mathrm{d} x^N} (x).\]

For discrete signals \(f[n]\), this operator is approximated by

\[\mathbf{D}^N_G f [n] = f[n] *\frac{\mathrm{d}^N g}{\mathrm{d} x^N} \left(\frac{n}{h}\right),\]

where \(h\) is the spacing between samples and the operator \(*\) is now a discrete convolution.

Warning

The operator \(\mathbf{D}_{G} \circ \mathbf{D}_{G}\) is not directly related to \(\mathbf{D}_{G}^{2}\): Gaussian smoothing is performed twice in the case of the former, whereas it is performed only once in the case of the latter.

Note that in contrast with finite differences (see finite_difference()), Gaussian derivatives compute exact derivatives in the continuous domain, since Gaussians can be differentiated analytically. This derivative is then sampled in order to perform a discrete convolution.

This class inherits its methods from Stencil.

Example

import numpy as np
import matplotlib.pyplot as plt
from pyxu.operator import PartialDerivative
from pyxu.util.misc import peaks

x = np.linspace(-2.5, 2.5, 25)
xx, yy = np.meshgrid(x, x)
image = peaks(xx, yy)
arg_shape = image.shape  # Shape of our image
# Specify derivative order at each direction
df_dx = (1, 0) # Compute derivative of order 1 in first dimension
d2f_dy2 = (0, 2) # Compute derivative of order 2 in second dimension
d3f_dxdy2 = (1, 2) # Compute derivative of order 1 in first dimension and der. of order 2 in second dimension
# Instantiate derivative operators
diff1 = PartialDerivative.gaussian_derivative(order=df_dx, arg_shape=arg_shape, sigma=2.0)
diff2 = PartialDerivative.gaussian_derivative(order=d2f_dy2, arg_shape=arg_shape, sigma=2.0)
diff = PartialDerivative.gaussian_derivative(order=d3f_dxdy2, arg_shape=arg_shape, sigma=2.0)
# Compute derivatives
out1 = (diff1 * diff2)(image.ravel()).reshape(arg_shape)
out2 = diff(image.ravel()).reshape(arg_shape)
plt.figure()
plt.imshow(image),
plt.axis('off')
plt.colorbar()
plt.title('f(x,y)')
plt.figure()
plt.imshow(out1.T)
plt.axis('off')
plt.title(r'$\frac{\partial^{3} f(x,y)}{\partial x\partial y^{2}}$')
plt.figure()
plt.imshow(out2.T)
plt.axis('off')
plt.title(r'$\frac{\partial^{3} f(x,y)}{\partial x\partial y^{2}}$')

../../_images/linop-12_00.png — (`png`, `hires.png`, `pdf`)#

../../_images/linop-12_01.png — (`png`, `hires.png`, `pdf`)#

../../_images/linop-12_02.png — (`png`, `hires.png`, `pdf`)#

Gradient(arg_shape, directions=None, diff_method='fd', mode='constant', gpu=False, dtype=None, parallel=False, **diff_kwargs)[source]#

Gradient operator.

Notes

This operator stacks the first-order partial derivatives of a \(D\)-dimensional signal \(\mathbf{f} \in \mathbb{R}^{N_{0} \times \cdots \times N_{D-1}}\) along each dimension:

\[\begin{split}\boldsymbol{\nabla} \mathbf{f} = \begin{bmatrix} \frac{\partial \mathbf{f}}{\partial x_0} \\ \vdots \\ \frac{\partial \mathbf{f}}{\partial x_{D-1}} \end{bmatrix} \in \mathbb{R}^{D \times N_{0} \times \cdots \times N_{D-1}}\end{split}\]

The partial derivatives can be approximated by finite differences via the finite_difference() constructor or by the Gaussian derivative via gaussian_derivative() constructor. The parametrization of the partial derivatives can be done via the keyword arguments **diff_kwargs, which will default to the same values as the PartialDerivative constructor.

The gradient operator’s method unravel() allows reshaping the vectorized output gradient to [n_dirs, N0, ..., ND] (see the example below).

Parameters:

arg_shape (NDArrayShape) – Shape of the input array.
directions (Integer, list[Integer], None) – Gradient directions. Defaults to None, which computes the gradient for all directions.
diff_method ('gd', 'fd') –
Method used to approximate the derivative. Must be one of:
- ’fd’ (default): finite differences
- ’gd’: Gaussian derivative
mode (str, list[str]) –
Boundary conditions. Multiple forms are accepted:
- str: unique mode shared amongst dimensions. Must be one of:
  - ’constant’ (default): zero-padding
  - ’wrap’
  - ’reflect’
  - ’symmetric’
  - ’edge’
- tuple[str, …]: the d-th dimension uses mode[d] as boundary condition.
(See numpy.pad() for details.)
gpu (bool) – Input NDArray type (True for GPU, False for CPU). Defaults to False.
dtype (DType) – Working precision of the linear operator.
parallel (bool) – If True, use Dask to evaluate the different partial derivatives in parallel.
diff_kwargs (dict) – Keyword arguments to parametrize partial derivatives (see finite_difference() and gaussian_derivative())

Returns:

op – Gradient

Return type:

Example

import numpy as np
import matplotlib.pyplot as plt
from pyxu.operator import Gradient
from pyxu.util.misc import peaks

# Define input image
n = 100
x = np.linspace(-3, 3, n)
xx, yy = np.meshgrid(x, x)
image = peaks(xx, yy)
arg_shape = image.shape  # (1000, 1000)
# Instantiate gradient operator
grad = Gradient(arg_shape=arg_shape)

# Compute gradients
output = grad(image.ravel()) # shape = (2000000, )
df_dx, df_dy = grad.unravel(output) # shape = (2, 1000, 1000)

# Plot image
fig, axs = plt.subplots(1, 3, figsize=(15, 4))
im = axs[0].imshow(image)
axs[0].set_title("Image")
axs[0].axis("off")
plt.colorbar(im, ax=axs[0])

# Plot gradient
im = axs[1].imshow(df_dx)
axs[1].set_title(r"$\partial f/ \partial x$")
axs[1].axis("off")
plt.colorbar(im, ax=axs[1])
im = axs[2].imshow(df_dy)
axs[2].set_title(r"$\partial f/ \partial y$")
axs[2].axis("off")
plt.colorbar(im, ax=axs[2])

PartialDerivative(), Jacobian()

See also

Jacobian(arg_shape, n_channels, directions=None, diff_method='fd', mode='constant', gpu=False, dtype=None, parallel=False, **diff_kwargs)[source]#

Jacobian operator.

Notes

This operator computes the first-order partial derivatives of a \(D\)-dimensional vector-valued signal of \(C\) variables or channels \(\mathbf{f} = [\mathbf{f}_{0}, \ldots, \mathbf{f}_{C-1}]\) with \(\mathbf{f}_{c} \in \mathbb{R}^{N_{0} \times \cdots \times N_{D-1}}\).

The Jacobian of \(\mathbf{f}\) is computed via the gradient as follows:

\[\begin{split}\mathbf{J} \mathbf{f} = \begin{bmatrix} (\boldsymbol{\nabla} \mathbf{f}_{0})^{\top} \\ \vdots \\ (\boldsymbol{\nabla} \mathbf{f}_{C-1})^{\top} \\ \end{bmatrix} \in \mathbb{R}^{C \times D \times N_0 \times \cdots \times N_{D-1}}\end{split}\]

The partial derivatives can be approximated by finite differences via the finite_difference() constructor or by the Gaussian derivative via gaussian_derivative() constructor. The parametrization of the partial derivatives can be done via the keyword arguments **diff_kwargs, which will default to the same values as the PartialDerivative constructor.

The Jacobian operator’s method unravel() allows reshaping the vectorized output Jacobian to [..., n_channels, n_dirs, N0, ..., ND] (see the example below).

Remark

Pyxu’s convention when it comes to field-vectors, is to work with vectorized arrays. However, the memory order of these arrays before raveling should be [S_0, ..., S_D, C, N] shape, with S_0, ..., S_D being stacking dimensions, C being the number of variables or channels, and N being the size of the domain (for example, N = npixels_x * npixels_y for a 2D image).

Parameters:

arg_shape (NDArrayShape) – Shape of the input array.
n_channels (Integer) – Number of channels or variables of the input vector-valued signal. The Jacobian with n_channels==1 yields the gradient.
directions (Integer, list[Integer], None) – Gradient directions. Defaults to None, which computes the gradient for all directions.
diff_method ("gd", "fd") –
Method used to approximate the derivative. Must be one of:
- ’fd’ (default): finite differences
- ’gd’: Gaussian derivative
mode (str, list[str]) –
Boundary conditions. Multiple forms are accepted:
- str: unique mode shared amongst dimensions. Must be one of:
  - ’constant’ (default): zero-padding
  - ’wrap’
  - ’reflect’
  - ’symmetric’
  - ’edge’
- tuple[str, …]: the d-th dimension uses mode[d] as boundary condition.
(See numpy.pad() for details.)
gpu (bool) – Input NDArray type (True for GPU, False for CPU). Defaults to False.
dtype (DType) – Working precision of the linear operator.
parallel (bool) – If True, use Dask to evaluate the different partial derivatives in parallel.
diff_kwargs (dict) – Keyword arguments to parametrize partial derivatives (see finite_difference() and gaussian_derivative())

Returns:

op – Jacobian

Return type:

Example

import numpy as np
import matplotlib.pyplot as plt
from pyxu.operator import Jacobian
from pyxu.util.misc import peaks

x = np.linspace(-2.5, 2.5, 25)
xx, yy = np.meshgrid(x, x)
image = np.tile(peaks(xx, yy), (3, 1, 1))
jac = Jacobian(arg_shape=image.shape[1:], n_channels=image.shape[0])
out = jac.unravel(jac(image.ravel()))
fig, axes = plt.subplots(3, 2, figsize=(10, 15))
for i in range(3):
    for j in range(2):
        axes[i, j].imshow(out[i, j].T, cmap=["Reds", "Greens", "Blues"][j])
        axes[i, j].set_title(f"$\partial I_{{{['R', 'G', 'B'][j]}}}/\partial{{{['x', 'y'][j]}}}$")
plt.suptitle("Jacobian")

Gradient(), PartialDerivative()

See also

Divergence(arg_shape, directions=None, diff_method='fd', mode='constant', gpu=False, dtype=None, parallel=False, **diff_kwargs)[source]#

Divergence operator.

Notes

This operator computes the expansion or outgoingness of a \(D\)-dimensional vector-valued signal of \(C\) variables \(\mathbf{f} = [\mathbf{f}_{0}, \ldots, \mathbf{f}_{C-1}]\) with \(\mathbf{f}_{c} \in \mathbb{R}^{N_{0} \times \cdots \times N_{D-1}}\).

The Divergence of \(\mathbf{f}\) is computed via the adjoint of the gradient as follows:

\[\operatorname{Div} \mathbf{f} = \boldsymbol{\nabla}^{\ast} \mathbf{f} = \begin{bmatrix} \frac{\partial \mathbf{f}}{\partial x_0} + \cdots + \frac{\partial \mathbf{f}}{\partial x_{D-1}} \end{bmatrix} \in \mathbb{R}^{N_{0} \times \cdots \times N_{D-1}}\]

The partial derivatives can be approximated by finite differences via the finite_difference() constructor or by the Gaussian derivative via gaussian_derivative() constructor. The parametrization of the partial derivatives can be done via the keyword arguments **diff_kwargs, which will default to the same values as the PartialDerivative constructor.

When using finite differences to compute the Divergence (i.e., diff_method = "fd"), the divergence returns the adjoint of the gradient in reversed order:

For forward type divergence, the adjoint of the gradient of “backward” type is used.
For backward type divergence, the adjoint of the gradient of “forward” type is used.

For central type divergence, and for the Gaussian derivative method (i.e., diff_method = "gd"), the adjoint of the gradient of “central” type is used (no reversed order).

The divergence operator’s method unravel() allows reshaping the vectorized output divergence to [..., N0, ..., ND].

Parameters:

arg_shape (NDArrayShape) – Shape of the input array.
directions (Integer, list[Integer], None) – Divergence directions. Defaults to None, which computes the divergence for all directions.
diff_method ("gd", "fd") –
Method used to approximate the derivative. Must be one of:
- ’fd’ (default): finite differences
- ’gd’: Gaussian derivative
mode (str, list[str]) –
Boundary conditions. Multiple forms are accepted:
- str: unique mode shared amongst dimensions. Must be one of:
  - ’constant’ (default): zero-padding
  - ’wrap’
  - ’reflect’
  - ’symmetric’
  - ’edge’
- tuple[str, …]: the d-th dimension uses mode[d] as boundary condition.
(See numpy.pad() for details.)
gpu (bool) – Input NDArray type (True for GPU, False for CPU). Defaults to False.
dtype (DType) – Working precision of the linear operator.
parallel (bool) – If True, use Dask to evaluate the different partial derivatives in parallel.
diff_kwargs (dict) – Keyword arguments to parametrize partial derivatives (see finite_difference() and gaussian_derivative())

Returns:

op – Divergence

Return type:

Example

import numpy as np
import matplotlib.pyplot as plt
from pyxu.operator import Gradient, Divergence, Laplacian
from pyxu.util.misc import peaks

n = 100
x = np.linspace(-3, 3, n)
xx, yy = np.meshgrid(x, x)
image = peaks(xx, yy)
arg_shape = image.shape  # (1000, 1000)
grad = Gradient(arg_shape=arg_shape)
div = Divergence(arg_shape=arg_shape)
# Construct Laplacian via composition
laplacian1 = div * grad
# Compare to default Laplacian
laplacian2 = Laplacian(arg_shape=arg_shape)
output1 = laplacian1(image.ravel())
output2 = laplacian2(image.ravel())
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
im = axes[0].imshow(np.log(abs(output1)).reshape(*arg_shape))
axes[0].set_title("Laplacian via composition")
plt.colorbar(im, ax=axes[0])
im = axes[1].imshow(np.log(abs(output1)).reshape(*arg_shape))
axes[1].set_title("Default Laplacian")
plt.colorbar(im, ax=axes[1])

Gradient(), PartialDerivative()

See also

Hessian(arg_shape, directions='all', diff_method='fd', mode='constant', gpu=False, dtype=None, parallel=False, **diff_kwargs)[source]#

Hessian operator.

Notes

The Hessian matrix or Hessian of a \(D\)-dimensional signal \(\mathbf{f} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}\) is the square matrix of second-order partial derivatives:

\[\begin{split}\mathbf{H} \mathbf{f} = \begin{bmatrix} \dfrac{ \partial^{2}\mathbf{f} }{ \partial x_{0}^{2} } & \dfrac{ \partial^{2}\mathbf{f} }{ \partial x_{0}\,\partial x_{1} } & \cdots & \dfrac{ \partial^{2}\mathbf{f} }{ \partial x_{0} \, \partial x_{D-1} } \\ \dfrac{ \partial^{2}\mathbf{f} }{ \partial x_{1} \, \partial x_{0} } & \dfrac{ \partial^{2}\mathbf{f} }{ \partial x_{1}^{2} } & \cdots & \dfrac{ \partial^{2}\mathbf{f} }{\partial x_{1} \,\partial x_{D-1}} \\ \vdots & \vdots & \ddots & \vdots \\ \dfrac{ \partial^{2}\mathbf{f} }{ \partial x_{D-1} \, \partial x_{0} } & \dfrac{ \partial^{2}\mathbf{f} }{ \partial x_{D-1} \, \partial x_{1} } & \cdots & \dfrac{ \partial^{2}\mathbf{f} }{ \partial x_{D-1}^{2}} \end{bmatrix}\end{split}\]

The partial derivatives can be approximated by finite differences via the finite_difference() constructor or by the Gaussian derivative via gaussian_derivative() constructor. The parametrization of the partial derivatives can be done via the keyword arguments **diff_kwargs, which will default to the same values as the PartialDerivative constructor.

The Hessian being symmetric, only the upper triangular part at most needs to be computed. Due to the (possibly) large size of the full Hessian, 4 different options are handled:

[mode 0] directions is an integer, e.g.: directions=0

\[\partial^{2}\mathbf{f}/\partial x_{0}^{2}.\]
[mode 1] directions is a tuple of length 2, e.g.: directions=(0,1)

\[\partial^{2}\mathbf{f}/\partial x_{0}\partial x_{1}.\]
[mode 2] directions is a tuple of tuples, e.g.: directions=((0,0), (0,1))

\[\left(\frac{ \partial^{2}\mathbf{f} }{ \partial x_{0}^{2} }, \frac{ \partial^{2}\mathbf{f} }{ \partial x_{0}\partial x_{1} }\right).\]
[mode 3] directions = 'all' (default), computes the Hessian for all directions (only the upper triangular part of the Hessian matrix), in row order, i.e.:

\[\left(\frac{ \partial^{2}\mathbf{f} }{ \partial x_{0}^{2} }, \frac{ \partial^{2}\mathbf{f} } { \partial x_{0}\partial x_{1} }, \, \ldots , \, \frac{ \partial^{2}\mathbf{f} }{ \partial x_{D-1}^{2} }\right).\]

The shape of the output LinOp depends on the number of computed directions; by default (all directions), we have \(\mathbf{H} \mathbf{f} \in \mathbb{R}^{\frac{D(D-1)}{2} \times N_0 \times \cdots \times N_{D-1}}\).

The Hessian operator’s unravel() method allows reshaping the vectorized output Hessian to [n_dirs, N0, ..., ND] (see the example below).

Parameters:

arg_shape (NDArrayShape) – Shape of the input array.
directions (Integer, (Integer, Integer), ((Integer, Integer), ..., (Integer, Integer)), 'all') – Hessian directions. Defaults to all, which computes the Hessian for all directions. (See Notes.)
diff_method ("gd", "fd") –
Method used to approximate the derivative. Must be one of:
- ’fd’ (default): finite differences
- ’gd’: Gaussian derivative
mode (str, list[str]) –
Boundary conditions. Multiple forms are accepted:
- str: unique mode shared amongst dimensions. Must be one of:
  - ’constant’ (default): zero-padding
  - ’wrap’
  - ’reflect’
  - ’symmetric’
  - ’edge’
- tuple[str, …]: the d-th dimension uses mode[d] as boundary condition.
(See numpy.pad() for details.)
gpu (bool) – Input NDArray type (True for GPU, False for CPU). Defaults to False.
dtype (DType) – Working precision of the linear operator.
parallel (bool) – If True, use Dask to evaluate the different partial derivatives in parallel.
diff_kwargs (dict) – Keyword arguments to parametrize partial derivatives (see finite_difference() and gaussian_derivative())

Returns:

op – Hessian

Return type:

Example

import numpy as np
import matplotlib.pyplot as plt
from pyxu.operator import Hessian, PartialDerivative
from pyxu.util.misc import peaks

n = 100
x = np.linspace(-3, 3, n)
xx, yy = np.meshgrid(x, x)
image = peaks(xx, yy)
arg_shape = image.shape  # (1000, 1000)

# Instantiate Hessian operator
hessian = Hessian(arg_shape=arg_shape, directions="all")
# Compute Hessian
output = hessian(image.ravel()) # shape = (300000,)
d2f_dx2, d2f_dxdy, d2f_dy2 = hessian.unravel(output)

# Plot
fig, axs = plt.subplots(1, 4, figsize=(20, 4))
im = axs[0].imshow(image)
plt.colorbar(im, ax=axs[0])
axs[0].set_title("Image")
axs[0].axis("off")

im = axs[1].imshow(d2f_dx2)
plt.colorbar(im, ax=axs[1])
axs[1].set_title(r"$\partial^{2} f/ \partial x^{2}$")
axs[1].axis("off")

im = axs[2].imshow(d2f_dxdy)
plt.colorbar(im, ax=axs[2])
axs[2].set_title(r"$\partial^{2} f/ \partial x\partial y$")
axs[2].axis("off")

im = axs[3].imshow(d2f_dy2)
plt.colorbar(im, ax=axs[3])
axs[3].set_title(r"$\partial^{2} f/ \partial y^{2}$")
axs[3].axis("off")

Laplacian(arg_shape, directions=None, diff_method='fd', mode='constant', gpu=False, dtype=None, parallel=False, **diff_kwargs)[source]#

Laplacian operator.

Notes

The Laplacian of a \(D\)-dimensional signal \(\mathbf{f} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}\) is the sum of second-order partial derivatives across all input directions:

\[\sum_{d = 0}^{D-1} \dfrac{ \partial^{2}\mathbf{f} }{ \partial x_{d}^{2} }\]

The partial derivatives can be approximated by finite differences via the finite_difference() constructor or by the Gaussian derivative via gaussian_derivative() constructor. The parametrization of the partial derivatives can be done via the keyword arguments **diff_kwargs, which will default to scheme='central' and accuracy=2 for diff_method='fd' (finite difference), and the same values as the PartialDerivative constructor for diff_method='gd' (gaussian derivative).

The Laplacian operator’s method unravel() allows reshaping the vectorized output Laplacian to [..., N0, ..., ND] (see the example below).

Parameters:

arg_shape (NDArrayShape) – Shape of the input array.
directions (Integer, list[Integer], None) – Laplacian directions. Defaults to None, which computes the Laplacian with all directions.
diff_method ("gd", "fd") –
Method used to approximate the derivative. Must be one of:
- ’fd’ (default): finite differences
- ’gd’: Gaussian derivative
mode (str, list[str]) –
Boundary conditions. Multiple forms are accepted:
- str: unique mode shared amongst dimensions. Must be one of:
  - ’constant’ (default): zero-padding
  - ’wrap’
  - ’reflect’
  - ’symmetric’
  - ’edge’
- tuple[str, …]: the d-th dimension uses mode[d] as boundary condition.
(See numpy.pad() for details.)
gpu (bool) – Input NDArray type (True for GPU, False for CPU). Defaults to False.
dtype (DType) – Working precision of the linear operator.
parallel (bool) – If True, use Dask to evaluate the different partial derivatives in parallel.
diff_kwargs (dict) – Keyword arguments to parametrize partial derivatives (see finite_difference() and gaussian_derivative())

Returns:

op – Laplacian

Return type:

Example

import numpy as np
import matplotlib.pyplot as plt
from pyxu.operator import Laplacian
from pyxu.util.misc import peaks

# Define input image
n = 100
x = np.linspace(-3, 3, n)
xx, yy = np.meshgrid(x, x)
image = peaks(xx, yy)

arg_shape = image.shape  # (1000, 1000)
# Compute Laplacian
laplacian = Laplacian(arg_shape=arg_shape)
output = laplacian(image.ravel())
output = laplacian.unravel(output) # shape = (1, 1000, 1000)
# Plot
fig, axs = plt.subplots(1, 2, figsize=(10, 4))
im = axs[0].imshow(image)
plt.colorbar(im, ax=axs[0])
axs[0].set_title("Image")
axs[0].axis("off")

im = axs[1].imshow(output.squeeze())
plt.colorbar(im, ax=axs[1])
axs[1].set_title(r"$\partial^{2} f/ \partial x^{2}+\partial^{2} f/ \partial y^{2}$")
axs[1].axis("off")

fig.show()

DirectionalDerivative(arg_shape, order, directions, diff_method='fd', mode='constant', parallel=False, **diff_kwargs)[source]#

Directional derivative operator.

Notes

The first-order DirectionalDerivative of a \(D\)-dimensional signal \(\mathbf{f} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}\) applies a derivative along the direction specified by a constant unitary vector \(\mathbf{v} \in \mathbb{R}^D\):

\[\boldsymbol{\nabla}_\mathbf{v} \mathbf{f} = \sum_{i=0}^{D-1} v_i \frac{\partial \mathbf{f}}{\partial x_i} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}\]

or along spatially-varying directions \(\mathbf{v} = [\mathbf{v}_0, \ldots , \mathbf{v}_{D-1}]^\top \in \mathbb{R}^{D \times N_0 \times \cdots \times N_{D-1} }\) where each direction \(\mathbf{v}_{\cdot, i_0, \ldots , i_{D-1}} \in \mathbb{R}^D\) for any \(0 \leq i_d \leq N_d-1\) with \(0 \leq d \leq D-1\) is a unitary vector:

\[\boldsymbol{\nabla}_\mathbf{v} \mathbf{f} = \sum_{i=0}^{D-1} \mathbf{v}_i \odot \frac{\partial \mathbf{f}}{\partial x_i} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}},\]

where \(\odot\) denotes the Hadamard (elementwise) product.

Note that choosing \(\mathbf{v}= \mathbf{e}_d \in \mathbb{R}^D\) (the \(d\)-th canonical basis vector) amounts to the first-order PartialDerivative() operator applied along axis \(d\).

The second-order DirectionalDerivative \(\boldsymbol{\nabla}^2_{\mathbf{v}_{1, 2}}\) is obtained by composing the first-order directional derivatives \(\boldsymbol{\nabla}_{\mathbf{v}_{1}}\) and \(\boldsymbol{\nabla}_{\mathbf{v}_{2}}\):

\[\boldsymbol{\nabla}_{\mathbf{v}_{1}} (\boldsymbol{\nabla}_{\mathbf{v}_{2}} \mathbf{f}) = \boldsymbol{\nabla}_{\mathbf{v}_{1}} (\sum_{i=0}^{D-1} {v_{2}}_{i} \frac{\partial \mathbf{f}}{\partial x_i}) = \sum_{j=0}^{D-1} {v_{1}}_{j} \frac{\partial}{\partial x_j} (\sum_{i=0}^{D-1} {v_{2}}_{i} \frac{\partial \mathbf{f}}{\partial x_i}) = \sum_{i, j=0}^{D-1} {v_{1}}_{j} {v_{2}}_{i} \frac{\partial}{\partial x_j} \frac{\partial}{\partial x_i}\mathbf{f} = \mathbf{v}_{1}^{\top}\mathbf{H}\mathbf{f}\mathbf{v}_{2},\]

where \(\mathbf{H}\) is the discrete Hessian operator, implemented via Hessian.

Higher-order DirectionalDerivative \(\boldsymbol{\nabla}^{N}_\mathbf{v}\) can be obtained by composing the first-order directional derivative \(\boldsymbol{\nabla}_\mathbf{v}\) \(N\) times.

The directional derivative operator’s method unravel() allows reshaping the vectorized output directional derivative to [..., N0, ..., ND] (see the example below).

Warning

DirectionalDerivative() instances are not array module-agnostic: they will only work with NDArrays belonging to the same array module as directions. Inner computations may recast input arrays when the precision of directions does not match the user-requested precision.
directions are always normalized to be unit vectors.

Parameters:

arg_shape (NDArrayShape) – Shape of the input array.
order (Integer) – Which directional derivative (restricted to 1: First or 2: Second, see Notes).
directions (NDArray, list) –
For order=1, it can be a single direction (array of size \((D,)\), where \(D\) is the number of axes) or spatially-varying directions:
- array of size \((D, N_0 \times \ldots \times N_{D-1})\) for order=1, i.e., one direction per element in the gradient, and,
- array of size \((D * (D + 1) / 2, N_0 \times \ldots \times N_{D-1})\) for order=2, i.e., one direction per element in the Hessian.
For order=2, it can be a tuple/list of two single directions or spatially-varying dimensions of the same shape, which will compute \(\mathbf{v}_{1}^{\top}\mathbf{H}\mathbf{f}\mathbf{v}_{2}\) or a single direction, as for order=1, which will compute \(\mathbf{v}_{1}^{\top}\mathbf{H}\mathbf{f}\mathbf{v}_{1}\). Note that for order=2, even though directions are spatially-varying, no differentiation is performed for this parameter.
diff_method ("gd", "fd") –
Method used to approximate the derivative. Must be one of:
- ’fd’ (default): finite differences
- ’gd’: Gaussian derivative
mode (str, list[str]) –
Boundary conditions. Multiple forms are accepted:
- str: unique mode shared amongst dimensions. Must be one of:
  - ’constant’ (default): zero-padding
  - ’wrap’
  - ’reflect’
  - ’symmetric’
  - ’edge’
- tuple[str, …]: the d-th dimension uses mode[d] as boundary condition.
(See numpy.pad() for details.)
parallel (bool) – If True, use Dask to evaluate the different partial derivatives in parallel.
diff_kwargs (dict) – Keyword arguments to parametrize partial derivatives (see finite_difference() and gaussian_derivative())

Returns:

op – Directional derivative

Return type:

Example

import numpy as np
import matplotlib.pyplot as plt
from pyxu.operator import DirectionalDerivative
from pyxu.util.misc import peaks

x = np.linspace(-2.5, 2.5, 25)
xx, yy = np.meshgrid(x, x)
z = peaks(xx, yy)
directions = np.zeros(shape=(2, z.size))
directions[0, : z.size // 2] = 1
directions[1, z.size // 2:] = 1
dop = DirectionalDerivative(arg_shape=z.shape, order=1, directions=directions)
out = dop.unravel(dop(z.ravel()))
dop2 = DirectionalDerivative(arg_shape=z.shape, order=2, directions=directions)
out2 = dop2.unravel(dop2(z.ravel()))
fig, axs = plt.subplots(1, 3, figsize=(15, 5))
axs = np.ravel(axs)
h = axs[0].pcolormesh(xx, yy, z, shading="auto")
axs[0].quiver(x, x, directions[1].reshape(xx.shape), directions[0].reshape(xx.shape))
plt.colorbar(h, ax=axs[0])
axs[0].set_title("Signal and directions of first derivatives")

h = axs[1].pcolormesh(xx, yy, out.squeeze(), shading="auto")
plt.colorbar(h, ax=axs[1])
axs[1].set_title("First-order directional derivatives")

h = axs[2].pcolormesh(xx, yy, out2.squeeze(), shading="auto")
plt.colorbar(h, ax=axs[2])
axs[2].set_title("Second-order directional derivative")

Gradient(), DirectionalGradient()

See also

DirectionalGradient(arg_shape, directions, diff_method='fd', mode='constant', parallel=False, **diff_kwargs)[source]#

Directional gradient operator.

Notes

The DirectionalGradient of a \(D\)-dimensional signal \(\mathbf{f} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}\) stacks the directional derivatives of \(\mathbf{f}\) along a list of \(m\) directions \(\mathbf{v}_i\) for \(1 \leq i \leq m\):

\[\begin{split}\boldsymbol{\nabla}_{\mathbf{v}_1, \ldots ,\mathbf{v}_m} \mathbf{f} = \begin{bmatrix} \boldsymbol{\nabla}_{\mathbf{v}_1} \\ \vdots\\ \boldsymbol{\nabla}_{\mathbf{v}_m}\\ \end{bmatrix} \mathbf{f} \in \mathbb{R}^{m \times N_0 \times \cdots \times N_{D-1}},\end{split}\]

where \(\boldsymbol{\nabla}_{\mathbf{v}_i}\) is the first-order directional derivative along \(\mathbf{v}_i\) implemented with DirectionalDerivative(), with \(\mathbf{v}_i \in \mathbb{R}^D\) or \(\mathbf{v}_i \in \mathbb{R}^{D \times N_0 \times \cdots \times N_{D-1}}\).

Note that choosing \(m=D\) and \(\mathbf{v}_i = \mathbf{e}_i \in \mathbb{R}^D\) (the \(i\)-th canonical basis vector) amounts to the Gradient() operator.

The directional gradient operator’s method unravel() allows reshaping the vectorized output directional derivative to [..., n_dirs, N0, ..., ND] (see the example below).

Parameters:

arg_shape (NDArrayShape) – Shape of the input array.
directions (list[NDArray]) – List of directions, either constant (array of size \((D,)\)) or spatially-varying (array of size \((D, N_0, \ldots, N_{D-1})\))
diff_method ("gd", "fd") –
Method used to approximate the derivative. Must be one of:
- ’fd’ (default): finite differences
- ’gd’: Gaussian derivative
mode (str, list[str]) –
Boundary conditions. Multiple forms are accepted:
- str: unique mode shared amongst dimensions. Must be one of:
  - ’constant’ (default): zero-padding
  - ’wrap’
  - ’reflect’
  - ’symmetric’
  - ’edge’
- tuple[str, …]: the d-th dimension uses mode[d] as boundary condition.
(See numpy.pad() for details.)
parallel (bool) – If True, use Dask to evaluate the different partial derivatives in parallel.
diff_kwargs (dict) – Keyword arguments to parametrize partial derivatives (see finite_difference() and gaussian_derivative())

Returns:

op – Directional gradient

Return type:

Example

import numpy as np
import matplotlib.pyplot as plt
from pyxu.operator import DirectionalGradient
from pyxu.util.misc import peaks

x = np.linspace(-2.5, 2.5, 25)
xx, yy = np.meshgrid(x, x)
z = peaks(xx, yy)
directions1 = np.zeros(shape=(2, z.size))
directions1[0, :z.size // 2] = 1
directions1[1, z.size // 2:] = 1
directions2 = np.zeros(shape=(2, z.size))
directions2[1, :z.size // 2] = -1
directions2[0, z.size // 2:] = -1
arg_shape = z.shape
dop = DirectionalGradient(arg_shape=arg_shape, directions=[directions1, directions2])
out = dop.unravel(dop(z.ravel()))
plt.figure()
h = plt.pcolormesh(xx, yy, z, shading='auto')
plt.quiver(x, x, directions1[1].reshape(arg_shape), directions1[0].reshape(xx.shape))
plt.quiver(x, x, directions2[1].reshape(arg_shape), directions2[0].reshape(xx.shape), color='red')
plt.colorbar(h)
plt.title(r'Signal $\mathbf{f}$ and directions of derivatives')
plt.figure()
h = plt.pcolormesh(xx, yy, out[0], shading='auto')
plt.colorbar(h)
plt.title(r'$\nabla_{\mathbf{v}_0} \mathbf{f}$')
plt.figure()
h = plt.pcolormesh(xx, yy, out[1], shading='auto')
plt.colorbar(h)
plt.title(r'$\nabla_{\mathbf{v}_1} \mathbf{f}$')

Gradient(), DirectionalDerivative()

../../_images/linop-19_00.png — (`png`, `hires.png`, `pdf`)#

../../_images/linop-19_01.png — (`png`, `hires.png`, `pdf`)#

../../_images/linop-19_02.png — (`png`, `hires.png`, `pdf`)#

See also

DirectionalLaplacian(arg_shape, directions, weights=None, diff_method='fd', mode='constant', parallel=False, **diff_kwargs)[source]#

Directional Laplacian operator.

Sum of the second directional derivatives of a multi-dimensional array (at least two dimensions are required) along multiple directions for each entry of the array.

Notes

The DirectionalLaplacian of a \(D\)-dimensional signal \(\mathbf{f} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}\) sums the second-order directional derivatives of \(\mathbf{f}\) along a list of \(m\) directions \(\mathbf{v}_i\) for \(1 \leq i \leq m\):

\[\boldsymbol{\Delta}_{\mathbf{v}_1, \ldots ,\mathbf{v}_m} \mathbf{f} = \sum_{i=1}^m \boldsymbol{\nabla}^2_{\mathbf{v}_i} \mathbf{f} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}},\]

where \(\boldsymbol{\nabla}^2_{\mathbf{v}_i}\) is the second-order directional derivative along \(\mathbf{v}_i\) implemented with DirectionalDerivative().

Note that choosing \(m=D\) and \(\mathbf{v}_i = \mathbf{e}_i \in \mathbb{R}^D\) (the \(i\)-th canonical basis vector) amounts to the Laplacian() operator.

The directional Laplacian operator’s method unravel() allows reshaping the vectorized output directional Laplacian to [..., N0, ..., ND] (see the example below).

Parameters:

arg_shape (NDArrayShape) – Shape of the input array.
directions (list[NDArray]) – List of directions, either constant (array of size \((D,)\)) or spatially-varying (array of size \((D, N_0, \ldots, N_{D-1})\))
weights (list[Real], None) – List of optional positive weights with which each second directional derivative operator is multiplied.
diff_method ("gd", "fd") –
Method used to approximate the derivative. Must be one of:
- ’fd’ (default): finite differences
- ’gd’: Gaussian derivative
mode (str, list[str]) –
Boundary conditions. Multiple forms are accepted:
- str: unique mode shared amongst dimensions. Must be one of:
  - ’constant’ (default): zero-padding
  - ’wrap’
  - ’reflect’
  - ’symmetric’
  - ’edge’
- tuple[str, …]: the d-th dimension uses mode[d] as boundary condition.
(See numpy.pad() for details.)
parallel (bool) – If True, use Dask to evaluate the different partial derivatives in parallel.
diff_kwargs (dict) – Keyword arguments to parametrize partial derivatives (see finite_difference() and gaussian_derivative())

Returns:

op – Directional Laplacian

Return type:

Example

import numpy as np
import matplotlib.pyplot as plt
from pyxu.operator import DirectionalLaplacian
from pyxu.util.misc import peaks

x = np.linspace(-2.5, 2.5, 25)
xx, yy = np.meshgrid(x, x)
z = peaks(xx, yy)
directions1 = np.zeros(shape=(2, z.size))
directions1[0, :z.size // 2] = 1
directions1[1, z.size // 2:] = 1
directions2 = np.zeros(shape=(2, z.size))
directions2[1, :z.size // 2] = -1
directions2[0, z.size // 2:] = -1
arg_shape = z.shape
dop = DirectionalLaplacian(arg_shape=arg_shape, directions=[directions1, directions2])
out = dop.unravel(dop(z.ravel()))
plt.figure()
h = plt.pcolormesh(xx, yy, z, shading='auto')
plt.quiver(x, x, directions1[1].reshape(arg_shape), directions1[0].reshape(xx.shape))
plt.quiver(x, x, directions2[1].reshape(arg_shape), directions2[0].reshape(xx.shape), color='red')
plt.colorbar(h)
plt.title('Signal and directions of derivatives')
plt.figure()
h = plt.pcolormesh(xx, yy, out.squeeze(), shading='auto')
plt.colorbar(h)
plt.title('Directional Laplacian')

../../_images/linop-20_00.png — (`png`, `hires.png`, `pdf`)#

../../_images/linop-20_01.png — (`png`, `hires.png`, `pdf`)#

DirectionalHessian(arg_shape, directions, diff_method='gd', mode='constant', parallel=False, **diff_kwargs)[source]#

Directional Hessian operator.

Notes

The DirectionalHessian of a \(D\)-dimensional signal \(\mathbf{f} \in \mathbb{R}^{N_0 \times \cdots \times N_{D-1}}\) stacks the second-order directional derivatives of \(\mathbf{f}\) along a list of \(m\) directions \(\mathbf{v}_i\) for \(1 \leq i \leq m\):

\[\begin{split}\mathbf{H}_{\mathbf{v}_1, \ldots ,\mathbf{v}_m} \mathbf{f} = \begin{bmatrix} \boldsymbol{\nabla}^2_{\mathbf{v}_0} & \cdots & \boldsymbol{\nabla}_{\mathbf{v}_0} \boldsymbol{\nabla}_{\mathbf{v}_{m-1}} \\ \vdots & \ddots & \vdots \\ \boldsymbol{\nabla}_{\mathbf{v}_{m-1}} \boldsymbol{\nabla}_{\mathbf{v}_0} & \cdots & \boldsymbol{\nabla}^2_{\mathbf{v}_{m-1}} \end{bmatrix} \mathbf{f},\end{split}\]

where \(\boldsymbol{\nabla}_{\mathbf{v}_i}\) is the first-order directional derivative along \(\mathbf{v}_i\) implemented with DirectionalDerivative().

However, due to the symmetry of the Hessian, only the upper triangular part is computed in practice:

\[\begin{split}\mathbf{H}_{\mathbf{v}_1, \ldots ,\mathbf{v}_m} \mathbf{f} = \begin{bmatrix} \boldsymbol{\nabla}^2_{\mathbf{v}_0}\\ \boldsymbol{\nabla}_{\mathbf{v}_0} \boldsymbol{\nabla}_{\mathbf{v}_{1}} \\ \vdots \\ \boldsymbol{\nabla}^2_{\mathbf{v}_{m-1}} \end{bmatrix} \mathbf{f} \in \mathbb{R}^{\frac{m (m-1)}{2} \times N_0 \times \cdots \times N_{D-1}}\end{split}\]

Note that choosing \(m=D\) and \(\mathbf{v}_i = \mathbf{e}_i \in \mathbb{R}^D\) (the \(i\)-th canonical basis vector) amounts to the Hessian() operator.

The directional Hessian operator’s method unravel() allows reshaping the vectorized output directional Hessian to [..., n_dirs, N0, ..., ND] (see the example below).

Parameters:

arg_shape (NDArrayShape) – Shape of the input array.
directions (list[NDArray]) – List of directions, either constant (array of size \((D,)\)) or spatially-varying (array of size \((D, N_0, \ldots, N_{D-1})\))
diff_method ("gd", "fd") –
Method used to approximate the derivative. Must be one of:
- ’fd’ (default): finite differences
- ’gd’: Gaussian derivative
mode (str, list[str]) –
Boundary conditions. Multiple forms are accepted:
- str: unique mode shared amongst dimensions. Must be one of:
  - ’constant’ (default): zero-padding
  - ’wrap’
  - ’reflect’
  - ’symmetric’
  - ’edge’
- tuple[str, …]: the d-th dimension uses mode[d] as boundary condition.
(See numpy.pad() for details.)
parallel (bool) – If True, use Dask to evaluate the different partial derivatives in parallel.
diff_kwargs (dict) – Keyword arguments to parametrize partial derivatives (see finite_difference() and gaussian_derivative())

Returns:

op – Directional Hessian

Return type:

Example

import numpy as np
import matplotlib.pyplot as plt
from pyxu.operator import DirectionalHessian
from pyxu.util.misc import peaks

x = np.linspace(-2.5, 2.5, 25)
xx, yy = np.meshgrid(x, x)
z = peaks(xx, yy)
directions1 = np.zeros(shape=(2, z.size))
directions1[0, :z.size // 2] = 1
directions1[1, z.size // 2:] = 1
directions2 = np.zeros(shape=(2, z.size))
directions2[1, :z.size // 2] = -1
directions2[0, z.size // 2:] = -1
arg_shape = z.shape
d_hess = DirectionalHessian(arg_shape=arg_shape, directions=[directions1, directions2])
out = d_hess.unravel(d_hess(z.ravel()))
plt.figure()
h = plt.pcolormesh(xx, yy, z, shading='auto')
plt.quiver(x, x, directions1[1].reshape(arg_shape), directions1[0].reshape(xx.shape))
plt.quiver(x, x, directions2[1].reshape(arg_shape), directions2[0].reshape(xx.shape), color='red')
plt.colorbar(h)
plt.title(r'Signal $\mathbf{f}$ and directions of derivatives')
plt.figure()
h = plt.pcolormesh(xx, yy, out[0], shading='auto')
plt.colorbar(h)
plt.title(r'$\nabla^2_{\mathbf{v}_0} \mathbf{f}$')
plt.figure()
h = plt.pcolormesh(xx, yy, out[1], shading='auto')
plt.colorbar(h)
plt.title(r'$\nabla_{\mathbf{v}_0} \nabla_{\mathbf{v}_{1}} \mathbf{f}$')
plt.figure()
h = plt.pcolormesh(xx, yy, out[2], shading='auto')
plt.colorbar(h)
plt.title(r'$\nabla^2_{\mathbf{v}_1} \mathbf{f}$')

Hessian(), DirectionalDerivative()

../../_images/linop-21_00.png — (`png`, `hires.png`, `pdf`)#

../../_images/linop-21_01.png — (`png`, `hires.png`, `pdf`)#

../../_images/linop-21_02.png — (`png`, `hires.png`, `pdf`)#

../../_images/linop-21_03.png — (`png`, `hires.png`, `pdf`)#

See also

Tensor Products #

kron(A, B)[source]#

Kronecker product \(A \otimes B\) between two linear operators.

The Kronecker product \(A \otimes B\) is defined as

\[\begin{split}A \otimes B = \left[ \begin{array}{ccc} A_{11} B & \cdots & A_{1N_{A}} B \\ \vdots & \ddots & \vdots \\ A_{M_{A}1} B & \cdots & A_{M_{A}N_{A}} B \\ \end{array} \right],\end{split}\]

where \(A : \mathbb{R}^{N_{A}} \to \mathbb{R}^{M_{A}}\), and \(B : \mathbb{R}^{N_{B}} \to \mathbb{R}^{M_{B}}\).

Parameters:

A (OpT) – (mA, nA) linear operator
B (OpT) – (mB, nB) linear operator

Returns:

op – (mA*mB, nA*nB) linear operator.

Return type:

Notes

This implementation is matrix-free by leveraging properties of the Kronecker product, i.e. \(A\) and \(B\) need not be known explicitly. In particular \((A \otimes B) x\) and \((A \otimes B)^{*} x\) are computed implicitly via the relation:

\[\text{vec}(\mathbf{A}\mathbf{B}\mathbf{C}) = (\mathbf{C}^{T} \otimes \mathbf{A}) \text{vec}(\mathbf{B}),\]

where \(\mathbf{A}\), \(\mathbf{B}\), and \(\mathbf{C}\) are matrices.

khatri_rao(A, B)[source]#

Column-wise Khatri-Rao product \(A \circ B\) between two linear operators.

The Khatri-Rao product \(A \circ B\) is defined as

\[A \circ B = \left[ \begin{array}{ccc} \mathbf{a}_{1} \otimes \mathbf{b}_{1} & \cdots & \mathbf{a}_{N} \otimes \mathbf{b}_{N} \end{array} \right],\]

where \(A : \mathbb{R}^{N} \to \mathbb{R}^{M_{A}}\), \(B : \mathbb{R}^{N} \to \mathbb{R}^{M_{B}}\), and \(\mathbf{a}_{k}\) (repectively \(\mathbf{b}_{k}\)) denotes the \(k\)-th column of \(A\) (respectively \(B\)).

Parameters:

A (OpT) – (mA, n) linear operator
B (OpT) – (mB, n) linear operator

Returns:

op – (mA*mB, n) linear operator.

Return type: