Differential Operators¶

Module: pycsou.linop.diff

Discrete differential and integral operators.

This module provides differential operators for discrete signals defined over regular grids or arbitrary meshes (graphs).

Many of the linear operators provided in this module are derived from linear operators from PyLops.

1D Operators

`FirstDerivative`(size[, shape, axis, step, …])	First derivative.
`SecondDerivative`(size[, shape, axis, step, …])	Second derivative.
`GeneralisedDerivative`(size[, shape, axis, …])	Generalised derivative.
`Integration1D`(size[, shape, axis, step, dtype])	1D integral/cumsum operator.

2D Operators

`Gradient`(shape[, step, edge, dtype, kind])	Gradient.
`Laplacian`(shape[, weights, step, edge, dtype])	Laplacian.
`GeneralisedLaplacian`([shape, step, edge, …])	Generalised Laplacian operator.
`FirstDirectionalDerivative`(shape, directions)	First directional derivative.
`SecondDirectionalDerivative`(shape, directions)	Second directional derivative.
`DirectionalGradient`(…)	Directional gradient.
`DirectionalLaplacian`(…[, weights])	Directional Laplacian.

FirstDerivative(size: int, shape: Optional[tuple] = None, axis: int = 0, step: float = 1.0, edge: bool = True, dtype: str = 'float64', kind: str = 'forward') → pycsou.linop.base.PyLopLinearOperator[source]¶

First derivative.

This docstring was adapted from ``pylops.FirstDerivative``.

Approximates the first derivative of a multi-dimensional array along a specific axis using finite-differences.

Parameters

size (int) – Size of the input array.
shape (tuple) – Shape of the input array.
axis (int) – Axis along which to differentiate.
step (float) – Step size.
edge (bool) – For kind='centered', use reduced order derivative at edges (True) or ignore them (False).
dtype (str) – Type of elements in input array.
kind (str) – Derivative kind (forward, centered, or backward).

Returns

First derivative operator.

Return type

PyLopLinearOperator

Raises

ValueError – If shape and size are not compatible.
NotImplementedError – If kind is not one of: forward, centered, or backward.

Examples

>>> x = np.repeat([0,2,1,3,0,2,0], 10)
>>> Dop = FirstDerivative(size=x.size)
>>> y = Dop * x
>>> np.sum(np.abs(y) > 0)
6
>>> np.allclose(y, np.diff(x, append=0))
True

import numpy as np
import matplotlib.pyplot as plt
from pycsou.linop.diff import FirstDerivative

x = np.repeat([0,2,1,3,0,2,0], 10)
Dop_bwd = FirstDerivative(size=x.size, kind='backward')
Dop_fwd = FirstDerivative(size=x.size, kind='forward')
Dop_cent = FirstDerivative(size=x.size, kind='centered')
y_bwd = Dop_bwd * x
y_cent = Dop_cent * x
y_fwd = Dop_fwd * x
plt.figure()
plt.plot(np.arange(x.size), x)
plt.plot(np.arange(x.size), y_bwd)
plt.plot(np.arange(x.size), y_cent)
plt.plot(np.arange(x.size), y_fwd)
plt.legend(['Signal', 'Backward', 'Centered', 'Forward'])
plt.title('First derivative')
plt.show()

(Source code, png, hires.png, pdf)

Notes

The FirstDerivative operator applies a first derivative along a given axis of a multi-dimensional array using either a second-order centered stencil or first-order forward/backward stencils.

For simplicity, given a one dimensional array, the second-order centered first derivative is:

\[y[i] = (0.5x[i+1] - 0.5x[i-1]) / \text{step}\]

while the first-order forward stencil is:

\[y[i] = (x[i+1] - x[i]) / \text{step}\]

and the first-order backward stencil is:

\[y[i] = (x[i] - x[i-1]) / \text{step}.\]

SecondDerivative(size: int, shape: Optional[tuple] = None, axis: int = 0, step: float = 1.0, edge: bool = True, dtype: str = 'float64') → pycsou.linop.base.PyLopLinearOperator[source]¶

Second derivative.

This docstring was adapted from ``pylops.SecondDerivative``.

Approximates the second derivative of a multi-dimensional array along a specific axis using finite-differences.

Parameters

size (int) – Size of the input array.
shape (tuple) – Shape of the input array.
axis (int) – Axis along which to differentiate.
step (float) – Step size.
edge (bool) – Use reduced order derivative at edges (True) or ignore them (False).
dtype (str) – Type of elements in input array.

Returns

Second derivative operator.

Return type

PyLopLinearOperator

Raises

ValueError – If shape and size are not compatible.

Examples

>>> x = np.linspace(-2.5, 2.5, 100)
>>> z = np.piecewise(x, [x < -1, (x >= - 1) * (x<0), x>=0], [lambda x: -x, lambda x: 3 * x + 4, lambda x: -0.5 * x + 4])
>>> Dop = SecondDerivative(size=x.size)
>>> y = Dop * z

import numpy as np
import matplotlib.pyplot as plt
from pycsou.linop.diff import SecondDerivative

x = np.linspace(-2.5, 2.5, 200)
z = np.piecewise(x, [x < -1, (x >= - 1) * (x<0), x>=0], [lambda x: -x, lambda x: 3 * x + 4, lambda x: -0.5 * x + 4])
Dop = SecondDerivative(size=x.size)
y = Dop * z
plt.figure()
plt.plot(np.arange(x.size), z)
plt.title('Signal')
plt.show()

(Source code, png, hires.png, pdf)

../../_images/pycsou-linop-diff-2_00_00.png

plt.figure()
plt.plot(np.arange(x.size), y)
plt.title('Second Derivative')
plt.show()

(png, hires.png, pdf)

../../_images/pycsou-linop-diff-2_01_00.png

Notes

The SecondDerivative operator applies a second derivative to any chosen direction of a multi-dimensional array.

For simplicity, given a one dimensional array, the second-order centered second derivative is given by:

\[y[i] = (x[i+1] - 2x[i] + x[i-1]) / \text{step}^2.\]

GeneralisedDerivative(size: int, shape: Optional[tuple] = None, axis: int = 0, step: float = 1.0, edge: bool = True, dtype: str = 'float64', kind_op='iterated', kind_diff='centered', **kwargs) → pycsou.core.linop.LinearOperator[source]¶

Generalised derivative.

Approximates the generalised derivative of a multi-dimensional array along a specific axis using finite-differences.

Parameters

size (int) – Size of the input array.
shape (tuple) – Shape of the input array.
axis (int) – Axis along which to differentiate.
step (float) – Step size.
edge (bool) – For kind_diff = 'centered', use reduced order derivative at edges (True) or ignore them (False).
dtype (str) – Type of elements in input array.
kind_diff (str) – Derivative kind (forward, centered, or backward).
kind_op (str) –
Type of generalised derivative ('iterated', 'sobolev', 'exponential', 'polynomial'). Depending on the cases, the GeneralisedDerivative operator is defined as follows:
- 'iterated': \(\mathscr{D}=D^N\),
- 'sobolev': \(\mathscr{D}=(\alpha^2 \mathrm{Id}-D^2)^N\), with \(\alpha\in\mathbb{R}\),
- 'exponential': \(\mathscr{D}=(\alpha \mathrm{Id} + D)^N\), with \(\alpha\in\mathbb{R}\),
- 'polynomial': \(\mathscr{D}=\sum_{n=0}^N \alpha_n D^n\), with \(\{\alpha_0,\ldots,\alpha_N\} \subset\mathbb{R}\),
where \(D\) is the FirstDerivative() operator.
kwargs (Any) –
Additional arguments depending on the value of kind_op:
- 'iterated': kwargs={order: int} where order defines the exponent \(N\).
- 'sobolev', 'exponential': kwargs={order: int, constant: float} where order defines the exponent \(N\) and constant the scalar \(\alpha\in\mathbb{R}\).
- kind_op='polynomial': kwargs={coeffs: Union[np.ndarray, list, tuple]} where coeffs is an array containing the coefficients \(\{\alpha_0,\ldots,\alpha_N\} \subset\mathbb{R}\).

Returns

A generalised derivative operator.

Return type

pycsou.core.linop.LinearOperator

Raises

NotImplementedError – If kind_op is not one of: 'iterated', 'sobolev', 'exponential', 'polynomial'.

Examples

>>> x = np.linspace(-5, 5, 100)
>>> z = np.sinc(x)
>>> Dop = GeneralisedDerivative(size=x.size, kind_op='iterated', order=1, kind_diff='forward')
>>> D = FirstDerivative(size=x.size, kind='forward')
>>> np.allclose(Dop * z, D * z)
True

import numpy as np
import matplotlib.pyplot as plt
from pycsou.linop.diff import GeneralisedDerivative

x  = np.linspace(-2.5, 2.5, 500)
z = np.exp(-x ** 2)
Dop_it = GeneralisedDerivative(size=x.size, kind_op='iterated', order=4)
Dop_sob = GeneralisedDerivative(size=x.size, kind_op='sobolev', order=2, constant=1e-2)
Dop_exp = GeneralisedDerivative(size=x.size, kind_op='exponential', order=4, constant=-1e-2)
Dop_pol = GeneralisedDerivative(size=x.size, kind_op='polynomial',
                                coeffs= 1/2 * np.array([1e-8, 0, -2 * 1e-4, 0, 1]))
y_it = Dop_it * z
y_sob = Dop_sob * z
y_exp = Dop_exp * z
y_pol = Dop_pol * z
plt.figure()
plt.plot(x, z)
plt.title('Signal')
plt.show()

(Source code, png, hires.png, pdf)

../../_images/pycsou-linop-diff-3_00_00.png

plt.figure()
plt.plot(x, y_it)
plt.plot(x, y_sob)
plt.plot(x, y_exp)
plt.plot(x, y_pol)
plt.legend(['Iterated', 'Sobolev', 'Exponential', 'Polynomial'])
plt.title('Generalised derivatives')
plt.show()

(png, hires.png, pdf)

../../_images/pycsou-linop-diff-3_01_00.png

Notes

Problematic values at edges are set to zero.

FirstDirectionalDerivative(shape: tuple, directions: numpy.ndarray, step: Union[float, Tuple[float, ...]] = 1.0, edge: bool = True, dtype: str = 'float64', kind: str = 'centered') → pycsou.linop.base.PyLopLinearOperator[source]¶

First directional derivative.

Computes the directional derivative of a multi-dimensional array (at least two dimensions are required) along either a single common direction or different directions for each entry of the array.

Parameters

shape (tuple) – Shape of the input array.
directions (np.ndarray) – Single direction (array of size \(n_{dims}\)) or different directions for each entry (array of size \([n_{dims} \times (n_{d_0} \times ... \times n_{d_{n_{dims}}})]\)). Each column should be normalised.
step (Union[float, Tuple[float, ..]]) – Step size in each direction.
edge (bool) – For kind = 'centered', use reduced order derivative at edges (True) or ignore them (False).
dtype (str) – Type of elements in input vector.
kind (str) – Derivative kind (forward, centered, or backward).

Returns

Directional derivative operator.

Return type

pycsou.linop.base.PyLopLinearOperator

Examples

>>> x = np.linspace(-2.5, 2.5, 100)
>>> X,Y = np.meshgrid(x,x)
>>> Z = peaks(X, Y)
>>> direction = np.array([1,0])
>>> Dop = FirstDirectionalDerivative(shape=Z.shape, directions=direction, kind='forward')
>>> D = FirstDerivative(size=Z.size, shape=Z.shape, kind='forward')
>>> np.allclose(Dop * Z.flatten(), D * Z.flatten())
True

import numpy as np
import matplotlib.pyplot as plt
from pycsou.linop.diff import FirstDirectionalDerivative, FirstDerivative
from pycsou.util.misc import peaks

x  = np.linspace(-2.5, 2.5, 25)
X,Y = np.meshgrid(x,x)
Z = peaks(X, Y)
directions = np.zeros(shape=(2,Z.size))
directions[0, :Z.size//2] = 1
directions[1, Z.size//2:] = 1
Dop = FirstDirectionalDerivative(shape=Z.shape, directions=directions)
y = Dop * Z.flatten()

plt.figure()
h = plt.pcolormesh(X,Y,Z, shading='auto')
plt.quiver(x, x, directions[1].reshape(X.shape), directions[0].reshape(X.shape))
plt.colorbar(h)
plt.title('Signal and directions of derivatives')
plt.figure()
h = plt.pcolormesh(X,Y,y.reshape(X.shape), shading='auto')
plt.colorbar(h)
plt.title('Directional derivatives')
plt.show()