Multidimensional Maps¶

Module: pycsou.core.map

Abstract classes for multidimensional (potentially nonlinear) maps.

`Map`(shape[, is_linear, is_differentiable])	Base class for multidimensional maps.
`DifferentiableMap`(shape[, is_linear, …])	Base class for multidimensional differentiable maps.
`MapStack`(*maps, axis[, n_jobs, joblib_backend])	Stack maps together.
`DiffMapStack`(*diffmaps, axis[, n_jobs, …])	Stack differentiable maps together.

class Map(shape: Tuple[int, int], is_linear: bool = False, is_differentiable: bool = False)[source]¶

Bases: abc.ABC

Base class for multidimensional maps.

Any instance/subclass of this class must at least implement the abstract method __call__. This class supports the following arithmetic operators +, -, *, @, ** and /, implemented with the class methods __add__/__radd__, __sub__/__neg__, __mul__/__rmul__, __matmul__, __pow__, __truediv__. Such arithmetic operators can be used to add, substract, scale, compose, exponentiate or evaluate Map instances.

Examples

Consider four maps: two nonlinear functionals \(f_1:\mathbb{R}^{10}\to \mathbb{R}\), \(f_2:\mathbb{R}^{10}\to \mathbb{R}\) and two linear operators \(L_1:\mathbb{R}^{10}\to \mathbb{R}^{10}\), \(L_2:\mathbb{R}^{10}\to \mathbb{R}^{10}\).

>>> print(f1.shape, f2.shape, L1.shape, L2.shape)
(1, 10) (1, 10) (10, 10) (10, 10)

We can combine linearly/compose the maps with consistent domains/ranges:

>>> f3 = f1 / 3 + np.pi * f2
>>> np.allclose(f3(x), f1(x) / 3 + np.pi * f2(x))
True
>>> L3 = L1 * 3 - (L2 ** 2) / 6
>>> np.allclose(L3(x), L1(x) * 3 - (L2(L2(x))) / 6)
True
>>> f4 = f3 * L3
>>> np.allclose(f4(x), f3(L3(x)))
True

Note that multiplying a map with an array is the same as evaluating the map at the array.

>>> np.allclose(f3 * x, f3(x))
True

The multiplication operator @ can also be used in place of *, in compliance with Numpy’s interface:

>>> np.allclose(f3 * x, f3 @ x)
True
>>> np.allclose((f1 * L1)(x), (f1 @ L1)(x))
True

Finally, maps can be shifted via the method shifter:

>>> f5=f4.shifter(shift=2 * x)
>>> np.allclose(f5(x), f4(x + 2 * x))
True

__init__(shape: Tuple[int, int], is_linear: bool = False, is_differentiable: bool = False)[source]¶

Parameters

shape (Tuple[int, int]) – Shape of the map.
is_linear (bool) – Whether the map is linear or not.
is_differentiable (bool) – Whether the map is differentiable or not.

abstract __call__(arg: Union[numbers.Number, numpy.ndarray]) → Union[numbers.Number, numpy.ndarray][source]¶

Call self as a function.

Parameters: arg (Union[Number, np.ndarray]) – Argument of the map.
Returns: Value of arg through the map.
Return type: Union[Number, np.ndarray]

apply_along_axis(arr: numpy.ndarray, axis: int = 0) → numpy.ndarray[source]¶

Apply the map to 1-D slices along the given axis.

Parameters

arr (np.ndarray) – Input array.
axis (int) – Axis along which arr is sliced.

Returns

The output array. The shape of the latter is identical to the shape of arr, except along the specified axis dimension. This axis is removed, and replaced with new dimensions equal to self.shape[0]. If self.shape[0]==1 the output array will have one fewer dimensions than arr.

Return type

np.ndarray

Raises

ValueError – If arr.shape[axis] != self.shape[1].

Examples

>>> from pycsou.linop import DownSampling
>>> from pycsou.func import SquaredL2Norm
>>> D=DownSampling(size=20, downsampling_factor=2)
>>> arr=np.arange(200).reshape(20,2,5)
>>> out = D.apply_along_axis(arr, axis=0)
>>> out.shape
(10, 2, 5)
>>> ff = SquaredL2Norm(dim=2)
>>> out = ff.apply_along_axis(arr, axis=1)
>>> out.shape
(20, 5)

shifter(shift: Union[numbers.Number, numpy.ndarray]) → pycsou.core.map.MapShifted[source]¶

Returns a shifted version of the map.

Parameters: shift (Union[Number, np.ndarray]) – Shift vector.
Returns: Shifted map.
Return type: MapShifted

Notes

Let A be a Map instance and B=A.shifter(y) with y some vector in the domain of A. Then we have: B(x)=A(x+y).

__add__(other: pycsou.core.map.Map) → pycsou.core.map.MapSum[source]¶

Add a map with another map instance.

Parameters: other (Union['Map', np.ndarray]) – The other map or array to be added.
Returns: Sum of the map with another map.
Return type: MapSum
Raises: NotImplementedError – If other is not a map.

__mul__(other: Union[numbers.Number, Map, numpy.ndarray]) → Union[pycsou.core.map.MapComp, numpy.ndarray][source]¶

Multiply a map with another map, a scalar, or an array.

The behaviour of this method depends on the type of other:

If other is another map, then it returns the composition of the map with other.
If other is an array with compatible shape, then it calls the map on other.
If other is a scalar, then it multiplies the map with this scalar.

Parameters: other (Union[Number, 'Map', np.ndarray]) – Scalar, map or array with which to multiply.
Returns: Composition of the map with another map, or product with a scalar or array.
Return type: MapComp, np.ndarray
Raises: NotImplementedError – If other is not a scalar, a map or an array.

See also

__matmul__(), __rmul__()

__matmul__(other: pycsou.core.map.Map) → pycsou.core.map.MapComp[source]¶: Alias for __mul__() offered to comply with Numpy’s interface.

__neg__() → pycsou.core.map.MapComp[source]¶: Negates a map. Alias for self.__mul__(-1).

__sub__(other: Union[Map, numpy.ndarray]) → Union[MapSum, MapBias][source]¶: Substracts a map, scalar, or array to another map. Alias for self.__add__(other.__neg__()).

__pow__(power: int) → pycsou.core.map.MapComp[source]¶: Raise a map to a certain power. Alias for A*A*...*A with power multiplications.

__truediv__(scalar: numbers.Number) → pycsou.core.map.MapComp[source]¶: Divides a map by a scalar. Alias for self.__mul__(1 / scalar).

class DifferentiableMap(shape: Tuple[int, int], is_linear: bool = False, lipschitz_cst: float = inf, diff_lipschitz_cst: float = inf)[source]¶

Bases: pycsou.core.map.Map

Base class for multidimensional differentiable maps.

Any instance/subclass of this class must at least implement the abstract methods __call__ and jacobianT. The attributes lipschitz_cst, diff_lipschitz_cst and the method jacobianT are automatically updated using standard differentiation rules when scaling, composing , summing or shifting differentiable maps.

Examples

Consider three differentiable maps: a nonlinear differentiable functional \(f:\mathbb{R}^{10}\to \mathbb{R}\), and two linear operators \(L_1:\mathbb{R}^{10}\to \mathbb{R}^{10}\), \(L_2:\mathbb{R}^{10}\to \mathbb{R}^{10}\).

>>> print(f.lipschitz_cst, f.diff_lipschitz_cst)
inf 2
>>> print(np.round(L1.lipschitz_cst,1), np.round(L1.diff_lipschitz_cst,1))
2.0 2.0
>>> print(np.round(L2.lipschitz_cst,1), np.round(L2.diff_lipschitz_cst,1))
1.5 1.5
>>> L3 = 2 * L1 + (L2 ** 2) / 3
>>> np.allclose(L3.lipschitz_cst, 2 * L1.lipschitz_cst + (L2.lipschitz_cst ** 2)/3)
True
>>> map_ = f * L1
>>> np.allclose(map_.lipschitz_cst, f.lipschitz_cst * L1.lipschitz_cst)
True
>>> np.allclose(map_.jacobianT(x), L1.jacobianT(x) * f.jacobianT(L1(x)))
True

Notes

This class supports the following arithmetic operators +, -, *, @, ** and /, implemented with the class methods __add__/__radd__, __sub__/__neg__, __mul__/__rmul__, __matmul__, __pow__, __truediv__. Such arithmetic operators can be used to add, substract, scale, compose, exponentiate or evaluate DifferentiableMap instances.

__init__(shape: Tuple[int, int], is_linear: bool = False, lipschitz_cst: float = inf, diff_lipschitz_cst: float = inf)[source]¶

Parameters

shape (Tuple[int, int]) – Shape of the differentiable map.
is_linear (bool) – Whether the differentiable map is linear or not.
lipschitz_cst (float) – Lispchitz constant of the differentiable map if it exists/is known. Default to \(+\infty\).
diff_lipschitz_cst (float) – Lispchitz constant of the derivative of the differentiable map if it exists/is known. Default to \(+\infty\).

abstract jacobianT(arg: Union[numbers.Number, numpy.ndarray]) → LinearOperator[source]¶

Transpose of the Jacobian matrix of the differentiable map evaluated at arg.

Parameters: arg (Union[Number, np.ndarray]) – Point at which the transposed Jacobian matrix is evaluated.
Returns: Linear operator associated to the transposed Jacobian matrix.
Return type: LinearOperator

gradient(arg: Union[numbers.Number, numpy.ndarray]) → LinearOperator[source]¶: Alias for self.jacobianT.

compute_lipschitz_cst()[source]¶: User-implemented method to compute the Lipschitz constant of the map.

compute_diff_lipschitz_cst()[source]¶: User-implemented method to compute the Lipschitz constant of the derivative of the map.

shifter(shift: Union[numbers.Number, numpy.ndarray]) → pycsou.core.map.DiffMapShifted[source]¶

Returns a shifted version of the map.

Parameters: shift (Union[Number, np.ndarray]) – Shift vector.
Returns: Shifted map.
Return type: DiffMapShifted

See also

shifter()

class MapStack(*maps, axis: int, n_jobs: int = 1, joblib_backend: str = 'loky')[source]¶

Bases: pycsou.core.map.Map

Stack maps together.

This class constructs a map by stacking multiple maps together, either vertically (axis=0) or horizontally (axis=1):

Vertical stacking: Consider a collection \(\{L_i:\mathbb{R}^{N}\to \mathbb{R}^{M_i}, i=1,\ldots, k\}\) of maps. Their vertical stacking is defined as the operator

\[\begin{split}V:\begin{cases}\mathbb{R}^{N}\to \mathbb{R}^{M_1}\times \cdots \times\mathbb{R}^{M_k}\\ \mathbf{x}\mapsto (L_1\mathbf{x},\ldots, L_k\mathbf{x}). \end{cases}\end{split}\]
Horizontal stacking: Consider a collection \(\{L_i:\mathbb{R}^{N_i}\to \mathbb{R}^{M}, i=1,\ldots, k\}\) of maps. Their horizontal stacking is defined as the operator

\[\begin{split}H:\begin{cases}\mathbb{R}^{N_1}\times \cdots \times\mathbb{R}^{N_k}\to \mathbb{R}^{M}\\ (\mathbf{x}_1,\ldots, \mathbf{x}_k)\mapsto \sum_{i=1}^k L_i \mathbf{x}_i. \end{cases}\end{split}\]

Examples

Consider three maps \(K_1:\mathbb{R}^{10}\to \mathbb{R}\), \(K_2:\mathbb{R}^{20}\to \mathbb{R}\) and \(K_3:\mathbb{R}^{10}\to \mathbb{R}^{10}\)

>>> V = MapStack(K1, K3, axis=0)
>>> V.shape
(11, 10)
>>> np.allclose(V(x), np.concatenate((K1(x).flatten(), K3(x).flatten())))
True
>>> parV = MapStack(K1, K3, axis=0, n_jobs=-1)
>>> np.allclose(V(x), parV(x))
True
>>> H = MapStack(K1,K2, axis=1)
>>> H.shape
(1, 30)
>>> np.allclose(H(np.concatenate((x,y))), K1(x) + K2(y))
True
>>> parH = MapStack(K1,K2, axis=1, n_jobs=-1)
>>> np.allclose(H(np.concatenate((x,y))), parH(np.concatenate((x,y))))
True

__init__(*maps, axis: int, n_jobs: int = 1, joblib_backend: str = 'loky')[source]¶

Parameters

maps (Map) – List of maps to stack.
axis – Stacking direction: 0 for vertical and 1 for horizontal stacking.
n_jobs (int) – Number of cores to be used for parallel evaluation of the map stack. If n_jobs==1, the map stack is evaluated sequentially, otherwise it is evaluated in parallel. Setting n_jobs=-1 uses all available cores.
joblib_backend (str) – Joblib backend (more details here).

class MapVStack(*maps, n_jobs: int = 1, joblib_backend: str = 'loky')[source]¶

Bases: pycsou.core.map.MapStack

Alias for vertical stacking, equivalent to MapStack(*maps, axis=0).

Examples

>>> V1 = MapStack(K1, K3, axis=0)
>>> V2 = MapVStack(K1, K3)
>>> np.allclose(V1(x), V2(x))
True

__init__(*maps, n_jobs: int = 1, joblib_backend: str = 'loky')[source]¶

Parameters

maps (Map) – List of maps to stack.
n_jobs (int) – Number of cores to be used for parallel evaluation of the map stack. If n_jobs==1, the map stack is evaluated sequentially, otherwise it is evaluated in parallel. Setting n_jobs=-1 uses all available cores.
joblib_backend (str) –
Joblib backend (more details here).

class MapHStack(*maps, n_jobs: int = 1, joblib_backend: str = 'loky')[source]¶

Bases: pycsou.core.map.MapStack

Alias for horizontal stacking, equivalent to MapStack(*maps, axis=1).

Examples

>>> V1 = MapStack(K1, K2, axis=1)
>>> V2 = MapHStack(K1, K2)
>>> np.allclose(V1(np.arange(V1.shape[1])), V2(np.arange(V1.shape[1])))
True

__init__(*maps, n_jobs: int = 1, joblib_backend: str = 'loky')[source]¶

Parameters

maps (Map) – List of maps to stack.
n_jobs (int) – Number of cores to be used for parallel evaluation of the map stack. If n_jobs==1, the map stack is evaluated sequentially, otherwise it is evaluated in parallel. Setting n_jobs=-1 uses all available cores.
joblib_backend (str) –
Joblib backend (more details here).

class DiffMapStack(*diffmaps, axis: int, n_jobs: int = 1, joblib_backend: str = 'loky')[source]¶

Bases: pycsou.core.map.MapStack, pycsou.core.map.DifferentiableMap

Stack differentiable maps together.

This class constructs a differentiable map by stacking multiple maps together, either vertically (axis=0) or horizontally (axis=1). The attributes lipschitz_cst, diff_lipschitz_cst and the method jacobianT are inferred from those of the stacked maps.

Vertical stacking: Consider a collection \(\{L_i:\mathbb{R}^{N}\to \mathbb{R}^{M_i}, i=1,\ldots, k\}\) of maps, with Lipschitz constants \(\{\beta_i\in\mathbb{R}_+\cup\{+\infty\}, i=1,\ldots, k\}\), Jacobian matrices \(\{\mathbf{J}_i(\mathbf{x})\in\mathbb{R}^{M_i\times N}, i=1,\ldots, k\}\) for some \(\mathbf{x}\in\mathbb{R}^N\). Then the vertically stacked operator

\[\begin{split}V:\begin{cases}\mathbb{R}^{N}\to \mathbb{R}^{M_1}\times \cdots \times\mathbb{R}^{M_k}\\ \mathbf{x}\mapsto (L_1\mathbf{x},\ldots, L_k\mathbf{x}). \end{cases}\end{split}\]

has a Lipschitz constant bounded by \(\sqrt{\sum_{i=1}^k \beta_i^2}\). Moreover the Jacobian matrix of \(V\) is obtained by stacking the individual Jacobian matrices vertically:

\[\begin{split}\mathbf{J}(\mathbf{x})=\left[\begin{array}{c}\mathbf{J}_1(\mathbf{x})\\\vdots\\ \mathbf{J}_k(\mathbf{x}) \end{array}\right]\in\mathbb{R}^{(\sum_{i=1}^k M_i)\times N}.\end{split}\]
Horizontal stacking: Consider a collection \(\{L_i:\mathbb{R}^{N_i}\to \mathbb{R}^{M}, i=1,\ldots, k\}\) of maps, with Lipschitz constants \(\{\beta_i\in\mathbb{R}_+\cup\{+\infty\}, i=1,\ldots, k\}\), Jacobian matrices \(\{\mathbf{J}_i(\mathbf{x}_i)\in\mathbb{R}^{M\times N_i}, i=1,\ldots, k\}\) for some \(\mathbf{x}_i\in\mathbb{R}^{N_i}\). Then the horizontally stacked operator

\[\begin{split}H:\begin{cases}\mathbb{R}^{N_1}\times \cdots \times\mathbb{R}^{N_k}\to \mathbb{R}^{M}\\ (\mathbf{x}_1,\ldots, \mathbf{x}_k)\mapsto \sum_{i=1}^k L_i \mathbf{x}_i. \end{cases}\end{split}\]

has a Lipschitz constant bounded by \(\max_{i=1}^k \beta_i\). Moreover the Jacobian matrix of \(H\) is obtained by stacking the individual Jacobian matrices horizontally:

\[\mathbf{J}(\mathbf{x}_1,\ldots,\mathbf{x}_k)=\left[\begin{array}{c}\mathbf{J}_1(\mathbf{x}_1)&\cdots& \mathbf{J}_k(\mathbf{x}_k) \end{array}\right]\in\mathbb{R}^{M\times (\sum_{i=1}^k N_i)}.\]

Examples

>>> V = DiffMapStack(K1, K3, axis=0)
>>> np.allclose(V.lipschitz_cst, np.sqrt(K1.lipschitz_cst ** 2 + K3.lipschitz_cst ** 2))
True
>>> parV = DiffMapStack(K1, K3, axis=0, n_jobs=-1)
>>> np.allclose(V.jacobianT(x) * np.ones(shape=(V.jacobianT(x).shape[1])), parV.jacobianT(x) * np.ones(shape=(V.jacobianT(x).shape[1])))
True

__init__(*diffmaps, axis: int, n_jobs: int = 1, joblib_backend: str = 'loky')[source]¶

Parameters

diffmaps (DifferentiableMap) – List of differentiable maps to be stacked
axis (int) – Stacking direction: 0 for vertical and 1 for horizontal stacking.
n_jobs (int) – Number of cores to be used for parallel evaluation of the map stack and its Jacobian matrix. If n_jobs==1, the map stack and its Jacobian matrix are evaluated sequentially, otherwise they are evaluated in parallel. Setting n_jobs=-1 uses all available cores.
joblib_backend (str) –
Joblib backend (more details here).

jacobianT(arg: Union[numbers.Number, numpy.ndarray]) → Union[LinOpHStack, LinOpVStack][source]¶

Transpose of the Jacobian matrix of the differentiable map evaluated at arg.

Parameters: arg (Union[Number, np.ndarray]) – Point at which the transposed Jacobian matrix is evaluated.
Returns: Linear operator associated to the transposed Jacobian matrix.
Return type: LinearOperator

class DiffMapVStack(*diffmaps, n_jobs: int = 1, joblib_backend: str = 'loky')[source]¶

Bases: pycsou.core.map.DiffMapStack

Alias for vertical stacking of differentiable maps, equivalent to DiffMapStack(*maps, axis=0).

__init__(*diffmaps, n_jobs: int = 1, joblib_backend: str = 'loky')[source]¶

Parameters

diffmaps (DifferentiableMap) – List of differentiable maps to be stacked
n_jobs (int) – Number of cores to be used for parallel evaluation of the map stack and its Jacobian matrix. If n_jobs==1, the map stack and its Jacobian matrix are evaluated sequentially, otherwise they are evaluated in parallel. Setting n_jobs=-1 uses all available cores.
joblib_backend (str) –
Joblib backend (more details here).

class DiffMapHStack(*diffmaps, n_jobs: int = 1, joblib_backend: str = 'loky')[source]¶

Bases: pycsou.core.map.DiffMapStack

Alias for horizontal stacking of differentiable maps, equivalent to DiffMapStack(*maps, axis=1).

__init__(*diffmaps, n_jobs: int = 1, joblib_backend: str = 'loky')[source]¶

Parameters

diffmaps (DifferentiableMap) – List of differentiable maps to be stacked.
n_jobs (int) – Number of cores to be used for parallel evaluation of the map stack and its Jacobian matrix. If n_jobs==1, the map stack and its Jacobian matrix are evaluated sequentially, otherwise they are evaluated in parallel. Setting n_jobs=-1 uses all available cores.
joblib_backend (str) –
Joblib backend (more details here).