Using operators" /> 2. Operators - PyOperators

2. Operators

2.1. Using operators

The Operator class and its subclasses are function factories. It means that in order to perform computations using an operator class, we first need to create an instance of it:

A = FFTOperator(1024)

and then use the instance as a function. The object A is a Python callable and takes numpy ’s N-dimensional array object as inputs and outputs:

>>> pyoperators.memory.verbose = True
>>> x = np.arange(1024, dtype=complex)
>>> y = A(x)
Info: Allocating (1024,) complex128 = 0.015625 MiB in FFTOperator.

When using the operator in such a way, a no-side-effect policy is enforced: it is garanteed that the input array will not be modified and as a consequence, a new buffer will be allocated every time the operator is called.

Similarly to ufuncs, such behaviour can be changed by providing the output array:

>>> pyoperators.memory.verbose = True
>>> y = np.empty(1024, dtype=complex)
>>> A(x, out=y)    # or simply A(x, y)

In-place operations can be performed in a similar way by providing the input array as the output.

In this documentation, the same term operator will be used for the Operator class, subclasses and instances, though the context will disambiguate the meaning.

2.2. Manipulating operators

Operators are easy to manipulate: they can be multiplied by a scalar, combined by addition, element-wise multiplication, composition or partition. The synthax uses the usual arithmetic signs, but be careful that the * sign depends on the operands’ linear flag: if both operands are linear, it stands for composition otherwise it stands for element-wise multiplication.

>>> A = DiagonalOperator([1., 2.])
>>> B = DiagonalOperator([1., 1.])
>>> A.flags.linear, B.flags.linear
(True, True)
>>> x = [2, 2]
>>> np.allclose((3 * A)(x), 3 * A(x))
True
>>> np.allclose((A + B)(x), A(x) + B(x))
True
>>> np.allclose(A(B)(x), A(B(x)))
True
>>> np.allclose((A * B)(x), A(B(x)))  # because A and B are linear
True
>>> np.allclose(MultiplicationOperator([A, B])(x), A(x) * B(x))
True

Unless an algebraic simplification is possible, the resulting operator is a composite operator:

addition AdditionOperator \( A(x) + B(x) \)
element-wise multiplication MultiplicationOperator \( A(x) \times B(x) \)
composition CompositionOperator \( A(B(x)) \)
partition BlockColumnOperator \( \begin{bmatrix} A \\ B \end{bmatrix} \)
BlockDiagonalOperator \( \begin{bmatrix} A & 0 \\ 0 & B \end{bmatrix} \)
BlockRowOperator \( \begin{bmatrix} A & B \end{bmatrix} \)

The conjugate, transpose, adjoint and inverse operators, when defined, can be obtained by using the following attributes:

conjugate .C
transpose .T
adjoint .H
inverse .I

For example:

>>> o = Operator()
>>> o.T.C is o.H
True
>>> o = Operator(flags='real,symmetric')
>>> o.T is o.H is o
True

2.3. Creating a new operator

Operators can be created in two ways. The first one is to define a direct function which will replace the usual matrix-vector operation and to instantiate the Operator class by passing the function as the first argument:

>>> def f(x, out):
...     out[...] = 2 * x
>>> P = Operator(f)

Transforming a single-argument ufunc into an Operator instance is straightforward:

>>> sqrt = Operator(np.sqrt)

The alternative way is more flexible, it consists in subclassing the Operator class:

>>> class MyOperator(Operator):
...     def __init__(self, coef):
...         self.coef = coef          
...         Operator.__init__(self)
...     def direct(self, x, out):
...         out[...] = self.coef * x
...
>>> P = MyOperator(2)

2.4. Operator’s data type

The operator’s data type is a data type object (an instance of numpy.dtype class). It can be set using the dtype keyword.

>>> a = Operator(dtype=np.float32)
>>> a.dtype
dtype('float32')

It is used to determine the output data type, as being the common type of the operator and input data types, following the standard coercion rules. For example:

>>> a = DiagonalOperator([1, 2])
>>> a.dtype
dtype('int64')
>>> a([1, 1j]).dtype
dtype('complex128')
>>> a(np.array([2, 2], np.uint8)).dtype
dtype('int64')

The operator’s data type can be equal to None (the default), in which case the output data type is always the input data type.

2.5. Operator’s flags

An operator P can have the following algebraic properties, accessible with the flags attribute:

linear
real conj P = P
symmetric P T = P
hermitian P H = P
idempotent P P = P
involutary P P = I
orthogonal P T P = I
unitary P H P = I

They can be set using the flags keyword:

>>> def f(x, out):
...     out[...] = -x
>>> P = Operator(direct=f, flags='symmetric,involutary')

or by using operator decorators:

>>> @pyoperators.flags.symmetric
... @pyoperators.flags.involutary
... class MyOperator(Operator):
...     def direct(self, x, out):
...         out[...] = -x
...
>>> P = MyOperator()

These flags are used to simplify expressions:

>>> C = Operator(flags='idempotent')
>>> C * C is C
True
>>> D = Operator(flags='involutary')
>>> D * D
IdentityOperator()

Other properties are described by the flags attribute:

square input and output have the same shape
separable can be separated if applied by or applying on a block operator
inplace handle inplace operations (more info)
update_output can update the operator output (more info)
alignment_input only handle aligned input
alignment_output only handle aligned output
contiguous_input only handle contiguous input
contiguous_output only handle contiguous output
shape_input explicit, implicit or unconstrained (more info)
shape_output explicit, implicit or unconstrained (more info)

2.6. Operator’s input and output shapes

The Operator class has a flexible way to deal with inputs and outputs. Its output shape can be:

explicit: the operator only returns arrays of a specified shape.

implicit: the output shape is derived from the input shape through the method reshapein.

unconstrained: the output shape does not depend on the input shape.

Likewise, its input shape can be:

explicit: the operator only applies over arrays of a specified shape.

implicit: the operator can only apply over arrays whose shape is derived from the output shape through the method reshapeout.

unconstrained: any input shape can be handled by the operator.

The operator’s flags record this information in the shape_input and shape_output attributes. If the output (input) shape is

Operator with explicit input and output shape.

The explicit shapes are specified using the shapein and shapeout keywords. The following example appends a zero to a fixed-size input.

>>> class Op(Operator):
...     def __init__(self, value):
...         Operator.__init__(self, shapein=3, shapeout=4)
...         self.value = value
...     def direct(self, input, output):
...         output[0:3] = input
...         output[3] = self.value
...
>>> op = Op(0.)
>>> op.shapein, op.shapeout
((3,), (4,))
>>> op([3., 2. ,1.])
array([3., 2., 1., 0.])

Operator with implicit input shape and explicit output shape.

This case does not happen because the input shape can be obtained from the explicit output shape and the reshapeout method.

Operator with unconstrained input shape and explicit output shape.

Such operator is obtained by only setting the shapeout keyword. The following example returns an array of size 2 whose elements are the sum and product of the input elements.

>>> class Op(Operator):
...     def __init__(self):
...         Operator.__init__(self, shapeout=2)
...     def direct(self, input, output):
...         output[0] = np.sum(input)
...         output[1] = np.product(input)
...
>>> op = Op()
>>> op.shapein, op.shapeout
(None, (2,))
>>> op([1., 2., 3., 4.])
array([ 10., 24.])

Operator with explicit input shape and implicit output shape

This case does not happen because the output shape can be obtained from the explicit input shape and the reshapein method.

Operator with implicit input and output shape

This operator defines both reshapein and reshapeout methods. Square-decorated operators are examples of this kind of operators: their reshapein and reshapeout methods simply are lambda x:x by default. The interest for an implicit output shape operator of also handling implicit input shapes arises when a reverse operator such as .T or .H is defined because in which the reshapein and reshapeout methods are swapped. The operator in the following example adds a zero to its input.

>>> @pyoperators.flags.linear
... class Op(Operator):
...     def direct(self, input, output):
...         output[:-1] = input
...         output[-1] = 0
...     def transpose(self, input, output):
...         output[...] = input[:-1]
...     def reshapein(self, shapein):
...         return (shapein[0]+1,)
...     def reshapeout(self, shapeout):
...         return (shapeout[0]-1,)
...
>>> op = Op()
>>> op.shapein, op.shapeout
(None, None)
>>> op([1,2,3])
array([1, 2, 3, 0])
>>> op.T([1, 2, 3, 0])
array([1, 2, 3])
>>> np.array_equal(op.T.todense(4), op.todense(3).T)
True

Operator with unconstrained input shape and implicit output shape

This operator can handle arbitrary inputs and its output dimensions are derived from those of the input. To obtain such operator, one simply has to define a reshapein method, which takes the input shape as input and returns the output shape. The following example stretches the input by a factor 2 along the first dimension.

>>> class Op(Operator):
...     def direct(self, input, output):
...         output[::2,...] = input
...         output[1::2,...] = input
...     def reshapein(self, shapein):
...         return (2*shapein[0],) + shapein[1:]
...
>>> op = Op()
>>> op.shapein, op.shapeout
(None, None)
>>> op.reshapein((2,3))
(4,3)
>>> op([[1, 2, 3],
...     [2, 3, 4]])
array([[1, 2, 3],
       [1, 2, 3],
       [2, 3, 4],
       [2, 3, 4]])

Operator with explicit input shape and unconstrained output shape

Such operator is obtained by only setting the shapein keyword. The following example fills the output with an arithmetic progression whose coefficients are given by the two elements of the input.

>>> class Op(Operator):
...     def __init__(self):
...         Operator.__init__(self, shapein=2)
...     def direct(self, input, output):
...         output[...] = input[0] + input[1] * np.arange(output.size).reshape(output.shape)
...
>>> op = Op()
>>> op.shapein, op.shapeout
((2,), None)
>>> y = np.empty((2,2), int)
>>> op([1, 2], out=y)
array([[1, 3],
       [5, 7]])

Operator with implicit input shape and unconstrained output shape

This kind of operator is of limited interest and is shown here for completeness. Since there is no way to compute the output shape from the input, the output argument must be provided. The operator in this example returns the input from which the last element has been removed.

>>> class Op(Operator):
...     def direct(self, input, output):
...         output[...] = input[:-1]
...     def reshapeout(self, shapeout):
...         return (shapeout[0] + 1,)
...
>>> op = Op()
>>> op.shapein, op.shapeout
(None, None)
>>> op([1,2,3])
ValueError: The output shape of an implicit input shape and unconstrained output shape operator cannot be inferred.
>>> y = np.empty(4)
>>> op([1.,2.,3.,4.,5.], out=y)
array([ 1.,  2.,  3.,  4.])
>>> op(np.arange(4.), out=y)
ValueError: The input has an invalid shape '(4,)'. Expected shape is '(5,)'.

Operator with unconstrained input and output shape

This kind of operator is the default: this behaviour is obtained by not setting the shapein and shapeout keywords and not defining the reshapein and reshapeout methods. ConstantOperator is an example of such an operator. If the output is not provided as argument, it is assumed that the output shape is that of the input. The following example fills the output with the value of the sum of the elements of the input.

>>> class Op(Operator):
...     def __init__(self):
...         Operator.__init__(self)
...     def direct(self, input, output):
...         output[...] = np.sum(input)
...
>>> op = Op()
>>> op.shapein, op.shapeout
(None, None)
>>> y = np.empty((2,2), int)
>>> op([2, 1], out=y)
array([[3, 3],
       [3, 3]])

2.7. Operator’s input and output shape validation.

It is possible to validate the shapes of the input or output of an operator by defining the validatein or validateout methods, which take the shapes of the input or output and raise a ValueError exception to signal a validation failure. There are two advantages of using these methods instead of performing the validation in the direct method:

The following example checks the number of axes and the dimension along the axes.

>>> class Op(Operator):
...     def direct(self, input, output):
...         output[...] = 1
...     def validatein(self, shapein):
...         if len(shapein) != 2:
...             raise ValueError('The number of axes of the input is not 2.')
...         if shapein[0] % 16 != 0 or shapein[1] % 16 != 0:
...             raise ValueError('The input dimensions are not a multiple of 16.')
...
>>> op = Op()
>>> op(np.ones(10))
ValueError: The number of axes of the input is not 2.
>>> op(np.ones(16,17))
ValueError: The input dimensions are not a multiple of 16.

2.8. Class propagation

An operator propagates the class of the input array, if the latter subclasses numpy.ndarray:

>>> class ndarray2(np.ndarray):
...     pass
>>> type(I(ndarray2((2,2))))
__main__.ndarray2

It is possible to change the subclass too, by using the classout keyword:

>>> class ndarray3(np.ndarray):
...     pass
>>> I2 = IdentityOperator(classout=ndarray3)
>>> type(I2(ndarray((2,2))))
__main__.ndarray3
>>> type(I2(ndarray2((2,2))))
__main__.ndarray3

2.9. Attribute propagation

An operator also propagates attributes:

>>> class ndarray2(np.ndarray):
...     pass
...
>>> x = ndarray2((2,2))
>>> x.foo = 'bar'
>>> I(x).foo
'bar'

It is also possible to add or change an attribute, by setting the attrout keyword to a dictionary whose keys are the attribute names and whose values are the attribute values:

>>> I2 = IdentityOperator(attrout={'foo':'new_bar'})
>>> I2(x).foo
'new_bar'
>>> x.foo
'bar'
>>> I2(np.ones((2,2))).foo
'new_bar'

More flexibility is possible by passing a function, instead of a dictionary, to the attrout keyword. This function expects a dictionary whose keys are the attribute names and whose values are the attribute values (it is the output’s __dict__):

>>> class I2Operator(IdentityOperator):
...     def __init__(self):
...         IdentityOperator.__init__(self, attrout=self.add_history)
...     def add_history(self, attr):
...         from time import ctime
...         if 'history' not in attr:
...             attr['history'] = []
...         attr['history'] += [ctime() + ' : ' + self.__class__.__name__]
...
>>> I2 = I2Operator()
>>> I2(I2(0)).history

2.10. Operator’s unary rules

The set of operators that can be obtained by using any number of times the '.C', '.T', '.H' and '.I' attributes over a given operator is finite. The set is made of the operator itself, its conjugate, transpose, adjoint, inverse, inverse-conjugate, inverse-transpose and inverse-adjoint.

It is possible to define such operators by using a unary rule of the form ‘subject’ is ‘predicate’. The subject being

'C' operator’s conjugate
'T' operator’s transpose
'H' operator’s adjoint
'I' operator’s inverse
'IC' operator’s inverse-conjugate
'IT' operator’s inverse-transpose
'IH' operator’s inverse-adjoint

and the predicate being

'.' the operator itself
function callback function

For instance, the rule 'C''.' means that the operator’s conjugate is the operator itself. This is how the real decorator is translated internally. Likewise, the involutary decorator is represented by the rule 'I''.'.

Unary rules are attached to an operator by using the set_rule method:

>>> @pyoperators.flags.linear
... class Op1(Operator):
...     def __init__(self, **keywords):
...         Operator.__init__(self, **keywords)
...         self.set_rule('T', lambda s: ReverseOperatorFactory(Op2, s))
...     def direct(self, input, output):
...         output[:-1] = input
...         output[-1] = 0
...     def reshapein(self, shapein):
...         return (shapein[0]+1,)
...
>>> @pyoperators.flags.linear
... class Op2(Operator):
...     def __init__(self, **keywords):
...         Operator.__init__(self, **keywords)
...         self.set_rule('T', lambda s: ReverseOperatorFactory(Op1, s))
...     def direct(self, input, output):
...         output[...] = input[:-1]
...     def reshapein(self, shapein):
...         return (shapein[0]-1,)
...
>>> op1 = Op1()
>>> op2 = op1.T
>>> type(op2)
__main__.Op2
>>> array_equal(op1.todense(4), op2.todense(5).T)
True

The ReverseOperatorFactory is a convenience factory which returns an operator of the type its first argument and by swapping the attributes attrin(out), classin(out), shapein(out), reshapein(out), toshapein(out), validatein(out) from the second argument.

2.11. Operator’s binary rules

PyOperators has the ability to reduce arithmetic expressions involving operators. This is achieved by means of binary rules which are applied to pairs of operators. Such a rule is specific to an operation (addition, multiplication or composition) and has the form ‘subject1 \& subject2 → predicate’, in which subject1 and subject2 are properties that will be matched against the pair of actual operators. One of them must be ‘.’, and stands for the reference operator to which the rule is attached. For instance and for a given operator P, the subject ('.', '.') will match the pair (P, P), the subject ('T', '.') the pair (P.T, P) and ('.', 'H') the pair (P, P.H).
In addition to the subjects that can be used for unary rules, it is possible to specify an Operator subclass for subclass matching. For example the subject (‘.’, Operator) will match any Operator instance on the right of the reference operator.

When the pair of operators is matched by the subject, it is replaced by the predicate:

'.' reference operator
'C' reference operator’s conjugate
'T' reference operator’s transpose
'H' reference operator’s ajoint
'I' reference operator’s inverse
'1' identity operator
function callback function

For a non-commutative operation (such as composition), the rules attached to the left-hand side operator applying on a right-hand-side operator are combined to the rules of the right-hand side operators applying on a left-hand-side operator:

Rules attached to an operator are also prioritised:

Let’s examine some examples for the composition operation:

The rules attached to an operator are stored in the rules attribute. It is a dictionary whose keys can be AdditionOperator, CompositionOperator or MultiplicationOperator. For the non-commutative composition operation, rules are further split according to whether the operator is matched on the left or right hand side of the pair.
Rules are added by using the set_rule method, which has 3 arguments:

The predicate can be a function, in which case it must be a static method with two arguments for the pair of operators that matches the subject. Rules, such as the ones set in the superclass, can be deleted with the del_rule method. In the following example, we construct an operator Power that raises its input to a specified power. We add a rule to reduce the composition of two instances of the Power class.

>>> @pyoperators.flags.square
... class Power(Operator):
...     def __init__(self, exponent):
...         Operator.__init__(self)
...         self.exponent = exponent
...         self.set_rule(('.', Power), self.rule_power, CompositionOperator)
...     def direct(input, output):
...         output[...] = input ** self.exponent
...     @staticmethod
...     def rule_power(p1, p2):
...         return Power(p1.exponent * p2.exponent)
...
>>> p1 = Power(2)
>>> p2 = Power(3)
>>> p = p1(p2)
>>> type(p)
__main__.Power
>>> p.exponent
6

2.12. In-place and out-of-place operators

in-place operator: operator whose direct method can handle input and output arguments pointing to the same memory location.

out-of-place operator: operator whose direct method cannot.

Such a property is set by the inplace decorator. For instance, the operator DiagonalOperator is an in-place operator and the following computation

>>> d = DiagonalOperator([1, 2])
>>> x = np.array([0, 1])
>>> d(x, x)
array([0, 2])

is done in-place and is equivalent to:

>>> x = np.array([0, 1])
>>> x *= [1, 2]

An in-place operator also has to handle out-of-place operations. In the following example, the operation is performed out-of-place to satisfy the no-side-effect policy:

>>> d = DiagonalOperator([1, 2])
>>> x = np.array([0, 1])
>>> y = d(x)
>>> y
array([0, 2])

which is equivalent to:

>>> x = np.array([0, 1])
>>> y = x * [1, 2]

The in-place property is useful to avoid intermediate variables. To minimise the memory-cache transfers during a composition, an algorithm has been put in place to determine the intermediate variables to be extracted from the memory manager, by maximising the temporal locality (i.e. by reusing as much as possible the same memory area). This algorithm depends on:

As an example, let’s consider the composition of an in-place operator IN by an out-of-place operator OUT. The out-of-place composition

>>> (IN * OUT)(x, out=y)

requires an intermediate variable only if the size of OUT ‘s output is larger than that of the y variable. Otherwise, y’s buffer is used as output for the OUT operator and the application of the operator IN is performed in-place on y.

Concerning the in-place composition:

>>> (IN * OUT)(x, out=x)

it is not possible to avoid the use of a temporary variable, since a buffer different from the x variable is required for the OUT operator. The application of the operator IN is performed out-of-place because its output is the x variable.

2.13. The flag update_output

Another possibility to avoid temporaries in the addition or element-wise multiplication of operators is to give an operator the possibility to directly update the running output argument in the direct method. It is particularly interesting for sums of matrices with few non-zero rows, since it gives a mechanism to skip the zero ones.
Let’s first see the default behaviour, in which the operands do not update their output argument:

>>> (o1 + o2)(x, out=y)

The following steps are performed:

If we now assume that the operator o2 can update its output argument, the intermediate variable is not required anymore:

One can enable such property by adding the flag update_output and the keyword operation to the operator’s direct method.

>>> import operator
>>> import pyoperators
>>> from pyoperators import Operator, operation_assignment
>>> pyoperators.memory.verbose=True
>>> @pyoperators.flags.linear
... @pyoperators.flags.update_output
... class P(Operator):
...     def __init__(self, index):
...          self.index = index
...          Operator.__init__(self, shapein=1, shapeout=11)
...     def direct(self, input, output, operation=operation_assignment):
...         if operation is operation_assignment:
...             output[...] = 0
...         elif operation is not operator.iadd:
...             raise NotImplementedError()
...         output[self.index] += input
...
>>> Q = sum(P(i) for i in [0, 5, 10])
>>> y = np.empty(11)
>>> Q([2.], out=y)
array([ 2.,  0.,  0.,  0.,  0.,  2.,  0.,  0.,  0.,  0.,  2.])

No intermediate variable is needed and all operations are performed in-place. The AdditionOperator Q calls the direct method of its first operand P with the operation keyword set to the function operation_assignment, which means that P ‘s output should be assigned to the output argument. In the subsequent calls, the operation keyword is set to the iadd function from Python’s standard operator module (If Q were a MultiplicationOperator, it would be imul), which results in adding the operand’s output to the output argument in-place. This example shows the interest of the flag update_output when only few rows are non-zero: the output only need to be updated for these rows, which avoids thrashing the cache with zeros everytime an operand is called.

2.14. Operator partition

Currently, there are two ways to partition operators, depending on whether the partitioning is done in an already existing dimension or not.

Stack partition: outputs are stacked along a new dimension
Chunk partition: outputs are concatenated along an existing dimension

Stack partition

A stack partition operator is an instance of BlockRowOperator initialised with the new_axisin keyword, BlockDiagonalOperator with the keywords new_axisin or new_axisout or BlockColumnOperator initialised with the new_axisout keyword. There is a strong constraint on the input and output shapes of the blocks, since they must be the same for all blocks. The partition along the specified axis is always explicit and is equal to a tuple of ones with as many elements as the number of blocks. In the following example is shown a stack partition block column operator:

>>> C = BlockColumnOperator([I, 2*I, 3*I], new_axisout=0)
>>> C(np.ones(2))
array([[ 1.,  1.],
       [ 2.,  2.],
       [ 3.,  3.]])

The outputs of the blocks are stacked along the dimension specified by new_axisout keyword, i.e the first one. The following example shows a stack partition block diagonal operator:

>>> D = BlockDiagonalOperator([I, 2*I, 3*I], new_axisin=-1)
>>> x = np.arange(2*3).reshape((2,3))
>>> x
array([[0, 1, 2],
       [3, 4, 5]])
>>> D(x)
array([[ 0,  2,  6],
       [ 3,  8, 15]])

It can be seen that the first block I is applied over x[:,0], the second block 2*I over x[:,1], and the third block 3*I over x[:,2].

Chunk partition

A chunk partition operator can be obtained as an instance of BlockRowOperator initialised with the axisin keyword, BlockDiagonalOperator with the keywords axisin or axisout or BlockColumnOperator initialised with the axisout keyword. There is a lesser constraint on the shapes of the partitioned input and output of the blocks: all dimensions must be same except along the partitioned dimension(s). The partition along the specified axis can be implicit or explicit, which is reflected in the keywords partitionin and partitionout. The following example shows a block diagonal chunk partition operator with an explicit partition

>>> D = BlockDiagonalOperator([I, 2*I, 3*I], axisin=-1, partitionin=(2,3,2))
>>> D(np.ones(7))
array([ 1.,  1.,  2.,  2.,  2.,  3.,  3.])

It can be seen that the first block I is applied over the first two elements, the second block 2*I over the three following elements, and the third block 3*I over the last two elements. Block row and block column chunk partition operators can handle implicit partitions:

>>> C = BlockColumnOperator([I, 2*I, 3*I], axisout=-1)
>>> C(np.ones((2,2)))
array([[ 1.,  1.,  2.,  2.,  3.,  3.],
       [ 1.,  1.,  2.,  2.,  3.,  3.]])

Likewise, the outputs of the blocks are concatenated along the dimension specified by the axisout keyword, i.e the last one in this case.

2.15. List of available operators

Linear operators

Non-linear operators

Composite operators

MPI operators

More operators can be found in the project PySimulators, projection, discrete differences, compression, downsampling, MPI operators, etc.
There is an on-going work consisting in migrating the most generic ones into the pyoperators package.