pymor.operators.numpy

Operators based on NumPy arrays.

This module provides the following NumPy-based Operators:

• NumpyMatrixOperator wraps a 2D NumPy array as an Operator.

• NumpyMatrixBasedOperator should be used as base class for all Operators which assemble into a NumpyMatrixOperator.

• NumpyGenericOperator wraps an arbitrary Python function between NumPy arrays as an Operator.

• NumpyCirculantOperator matrix-free operator of a circulant matrix.

• NumpyToeplitzOperator matrix-free operator of a Toeplitz matrix.

• NumpyHankelOperator matrix-free operator of a Hankel matrix.

Module Contents

class pymor.operators.numpy.NumpyCirculantOperator(c, name=None)

Bases: pymor.operators.interface.Operator, pymor.core.cache.CacheableObject

Matrix-free representation of a (block) circulant matrix by a NumPy array.

A (block) circulant matrix is a special kind of (block) Toeplitz matrix which is (block) square and completely determined by its first (matrix-valued) column via

\[\begin{split}C = \begin{bmatrix} c_1 & c_n & c_{n-1} & \cdots & \cdots & c_2 \\ c_2 & c_1 & c_n & & & \vdots \\ c_3 & c_2 & \ddots & & & \vdots \\ \vdots & & & \ddots & c_n & c_{n-1} \\ \vdots & & & c_2 & c_1 & c_n \\ c_n & \cdots & \cdots & c_3 & c_2 & c_1 \end{bmatrix} \in \mathbb{C}^{n*p \times n*m},\end{split}\]

where the so-called circulant vector \(c \in \mathbb{C}^{\times n\times p\times m}\) denotes the first (matrix-valued) column of the matrix. The matrix \(C\) as seen above is not explicitly constructed, only c is stored. Efficient matrix-vector multiplications are realized with DFT in the class’ apply method. See [GVL13] Chapter 4.8.2. for details.

Parameters:

• c – The NumPy array of shape (n) or (n, p, m) that defines the circulant vector.

• name – Name of the operator.

Methods

apply	Apply the operator to a VectorArray.
apply_adjoint	Apply the adjoint operator.

cache_region = 'memory'

linear = True

apply(U, mu=None)

Apply the operator to a VectorArray.

Parameters:

• U – VectorArray of vectors to which the operator is applied.

• mu – The parameter values for which to evaluate the operator.

Returns:

|VectorArray| of the operator evaluations.

apply_adjoint(V, mu=None)

Apply the adjoint operator.

For any given linear Operator op, parameter values mu and VectorArrays U, V in the source resp. range we have:

op.apply_adjoint(V, mu).dot(U) == V.inner(op.apply(U, mu))

Thus, when op is represented by a matrix M, apply_adjoint is given by left-multiplication of (the complex conjugate of) M with V.

Parameters:

• V – VectorArray of vectors to which the adjoint operator is applied.

• mu – The parameter values for which to apply the adjoint operator.

Returns:

|VectorArray| of the adjoint operator evaluations.

class pymor.operators.numpy.NumpyGenericOperator(mapping, adjoint_mapping=None, dim_source=1, dim_range=1, linear=False, parameters={}, solver_options=None, name=None)

Bases: pymor.operators.interface.Operator

Wraps an arbitrary Python function between NumPy arrays as an Operator.

Parameters:

• mapping –

  The function to wrap. If parameters is None, the function is of the form mapping(U) and is expected to be vectorized. In particular:

  mapping(U).shape == (dim_range,) + U.shape[1:].
  

  If parameters is not None, the function has to have the signature mapping(U, mu).

• adjoint_mapping –

  The adjoint function to wrap. If parameters is None, the function is of the form adjoint_mapping(U) and is expected to be vectorized. In particular:

  adjoint_mapping(U).shape == (dim_source,) + U.shape[1:].
  

  If parameters is not None, the function has to have the signature adjoint_mapping(U, mu).

• dim_source – Dimension of the operator’s source.

• dim_range – Dimension of the operator’s range.

• linear – Set to True if the provided mapping and adjoint_mapping are linear.

• parameters – The Parameters the operator depends on.

• solver_options – The solver_options for the operator.

• name – Name of the operator.

Methods

apply	Apply the operator to a VectorArray.
apply_adjoint	Apply the adjoint operator.

apply(U, mu=None)

Apply the operator to a VectorArray.

Parameters:

• U – VectorArray of vectors to which the operator is applied.

• mu – The parameter values for which to evaluate the operator.

Returns:

|VectorArray| of the operator evaluations.

apply_adjoint(V, mu=None)

Apply the adjoint operator.

For any given linear Operator op, parameter values mu and VectorArrays U, V in the source resp. range we have:

op.apply_adjoint(V, mu).dot(U) == V.inner(op.apply(U, mu))

Thus, when op is represented by a matrix M, apply_adjoint is given by left-multiplication of (the complex conjugate of) M with V.

Parameters:

• V – VectorArray of vectors to which the adjoint operator is applied.

• mu – The parameter values for which to apply the adjoint operator.

Returns:

|VectorArray| of the adjoint operator evaluations.

class pymor.operators.numpy.NumpyHankelOperator(c, r=None, name=None)

Bases: pymor.operators.interface.Operator

Matrix-free representation of a finite dimensional Hankel operator by a NumPy array.

A finite dimensional Hankel operator can be represented by a (block) matrix with constant anti-diagonals (anti-diagonal blocks), i.e.:

\[\begin{split}H = \begin{bmatrix} c_1 & c_2 & c_3 & \cdots & \cdots & r_1 \\ c_2 & c_3 & & & & \vdots \\ c_3 & & & & & \vdots \\ \vdots & & & & & r_{n-2} \\ \vdots & & & & r_{n-2} & r_{n-1} \\ c_m & \cdots & \cdots & r_{n-2} & r_{n-1} & r_n \end{bmatrix} \in \mathbb{C}^{n*p \times k*m},\end{split}\]

where \(c\in\mathbb{C}^{s\times p\times m}\) and \(r\in\mathbb{C}^{k\times p\times m}\) denote the first (matrix-valued) column and last (matrix-valued) row of the (block) Hankel matrix, respectively. The matrix \(H\) as seen above is not explicitly constructed, only the arrays c and r are stored. Efficient matrix-vector multiplications are realized with DFT in the class’ apply method (see [MSKC21] Algorithm 3.1. for details).

Parameters:

• c – The NumPy array of shape (n) or (n, p, m) that defines the first column of the (block) Hankel matrix.

• r – The NumPy array of shape (k,) or (k, p, m) that defines the last row of the (block) Hankel matrix. If supplied, its first entry r[0] has to be equal to c[-1]. Defaults to None. If r is None, the behavior of scipy.linalg.hankel is mimicked which sets r to zero (except for the first entry).

• name – Name of the operator.

Methods

apply	Apply the operator to a VectorArray.
apply_adjoint	Apply the adjoint operator.

linear = True

apply(U, mu=None)

Apply the operator to a VectorArray.

Parameters:

• U – VectorArray of vectors to which the operator is applied.

• mu – The parameter values for which to evaluate the operator.

Returns:

|VectorArray| of the operator evaluations.

apply_adjoint(V, mu=None)

Apply the adjoint operator.

For any given linear Operator op, parameter values mu and VectorArrays U, V in the source resp. range we have:

op.apply_adjoint(V, mu).dot(U) == V.inner(op.apply(U, mu))

Thus, when op is represented by a matrix M, apply_adjoint is given by left-multiplication of (the complex conjugate of) M with V.

Parameters:

• V – VectorArray of vectors to which the adjoint operator is applied.

• mu – The parameter values for which to apply the adjoint operator.

Returns:

|VectorArray| of the adjoint operator evaluations.

class pymor.operators.numpy.NumpyMatrixBasedOperator

Bases: pymor.operators.interface.Operator

Base class for operators which assemble into a NumpyMatrixOperator.

sparse

True if the operator assembles into a sparse matrix, False if the operator assembles into a dense matrix, None if unknown.

Methods

apply	Apply the operator to a VectorArray.
apply_adjoint	Apply the adjoint operator.
apply_inverse	Apply the inverse operator.
as_range_array	Return a VectorArray representation of the operator in its range space.
as_source_array	Return a VectorArray representation of the operator in its source space.
assemble	Assemble the operator for given parameter values.
export_matrix	Save the matrix of the operator to a file.

linear = True

apply(U, mu=None)

Apply the operator to a VectorArray.

Parameters:

• U – VectorArray of vectors to which the operator is applied.

• mu – The parameter values for which to evaluate the operator.

Returns:

|VectorArray| of the operator evaluations.

apply_adjoint(V, mu=None)

Apply the adjoint operator.

For any given linear Operator op, parameter values mu and VectorArrays U, V in the source resp. range we have:

op.apply_adjoint(V, mu).dot(U) == V.inner(op.apply(U, mu))

Thus, when op is represented by a matrix M, apply_adjoint is given by left-multiplication of (the complex conjugate of) M with V.

Parameters:

• V – VectorArray of vectors to which the adjoint operator is applied.

• mu – The parameter values for which to apply the adjoint operator.

Returns:

|VectorArray| of the adjoint operator evaluations.

apply_inverse(V, mu=None, initial_guess=None, least_squares=False)

Apply the inverse operator.

Parameters:

• V – VectorArray of vectors to which the inverse operator is applied.

• mu – The parameter values for which to evaluate the inverse operator.

• initial_guess – VectorArray with the same length as V containing initial guesses for the solution. Some implementations of apply_inverse may ignore this parameter. If None a solver-dependent default is used.

• least_squares –

  If True, solve the least squares problem:

  u = argmin ||op(u) - v||_2.
  

  Since for an invertible operator the least squares solution agrees with the result of the application of the inverse operator, setting this option should, in general, have no effect on the result for those operators. However, note that when no appropriate solver_options are set for the operator, most implementations will choose a least squares solver by default which may be undesirable.

Returns:

|VectorArray| of the inverse operator evaluations.

Raises:

InversionError – The operator could not be inverted.

as_range_array(mu=None)

Return a VectorArray representation of the operator in its range space.

In the case of a linear operator with NumpyVectorSpace as source, this method returns for given parameter values mu a VectorArray V in the operator’s range, such that

V.lincomb(U.to_numpy()) == self.apply(U, mu)

for all VectorArrays U.

Parameters:

mu – The parameter values for which to return the VectorArray representation.

Returns:

V – The VectorArray defined above.

as_source_array(mu=None)

Return a VectorArray representation of the operator in its source space.

In the case of a linear operator with NumpyVectorSpace as range, this method returns for given parameter values mu a VectorArray V in the operator’s source, such that

self.range.make_array(V.inner(U)) == self.apply(U, mu)

for all VectorArrays U.

Parameters:

mu – The parameter values for which to return the VectorArray representation.

Returns:

V – The VectorArray defined above.

assemble(mu=None)

Assemble the operator for given parameter values.

The result of the method strongly depends on the given operator. For instance, a matrix-based operator will assemble its matrix, a LincombOperator will try to form the linear combination of its operators, whereas an arbitrary operator might simply return a FixedParameterOperator. The only assured property of the assembled operator is that it no longer depends on a Parameter.

Parameters:

mu – The parameter values for which to assemble the operator.

Returns:

Parameter-independent, assembled |Operator|.

export_matrix(filename, matrix_name=None, output_format='matlab', mu=None)

Save the matrix of the operator to a file.

Parameters:

• filename – Name of output file.

• matrix_name – The name, the output matrix is given. (Comment field is used in case of Matrix Market output_format.) If None, the Operator’s name is used.

• output_format – Output file format. Either matlab or matrixmarket.

• mu – The parameter values to assemble the to be exported matrix for.

class pymor.operators.numpy.NumpyMatrixOperator(matrix, solver_options=None, name=None)

Bases: NumpyMatrixBasedOperator

Wraps a 2D NumPy array or SciPy spmatrix as an Operator.

Note

In the case of a NumPy array, the apply_inverse method by default uses check_finite=True and check_cond=True. Setting them to False (e.g., via defaults) can significantly speed up the computation, especially for smaller matrices.

Parameters:

• matrix – The NumPy array or SciPy spmatrix which is to be wrapped.

• solver_options – The solver_options for the operator.

• name – Name of the operator.

Methods

apply	Apply the operator to a VectorArray.
apply_adjoint	Apply the adjoint operator.
apply_inverse	Apply the inverse operator.
apply_inverse_adjoint	Apply the inverse adjoint operator.
as_range_array	Return a VectorArray representation of the operator in its range space.
as_source_array	Return a VectorArray representation of the operator in its source space.
assemble	Assemble the operator for given parameter values.
from_file

apply(U, mu=None)

Apply the operator to a VectorArray.

Parameters:

• U – VectorArray of vectors to which the operator is applied.

• mu – The parameter values for which to evaluate the operator.

Returns:

|VectorArray| of the operator evaluations.

apply_adjoint(V, mu=None)

Apply the adjoint operator.

For any given linear Operator op, parameter values mu and VectorArrays U, V in the source resp. range we have:

op.apply_adjoint(V, mu).dot(U) == V.inner(op.apply(U, mu))

Thus, when op is represented by a matrix M, apply_adjoint is given by left-multiplication of (the complex conjugate of) M with V.

Parameters:

• V – VectorArray of vectors to which the adjoint operator is applied.

• mu – The parameter values for which to apply the adjoint operator.

Returns:

|VectorArray| of the adjoint operator evaluations.

apply_inverse(V, mu=None, initial_guess=None, least_squares=False, check_finite=True, check_cond=True, default_sparse_solver_backend='scipy')

Apply the inverse operator.

Parameters:

• V – VectorArray of vectors to which the inverse operator is applied.

• mu – The parameter values for which to evaluate the inverse operator.

• initial_guess – VectorArray with the same length as V containing initial guesses for the solution. Some implementations of apply_inverse may ignore this parameter. If None a solver-dependent default is used.

• least_squares –

  If True, solve the least squares problem:

  u = argmin ||op(u) - v||_2.
  

  Since for an invertible operator the least squares solution agrees with the result of the application of the inverse operator, setting this option should, in general, have no effect on the result for those operators. However, note that when no appropriate solver_options are set for the operator, most implementations will choose a least squares solver by default which may be undesirable.

• check_finite – Test if solution only contains finite values.

• check_cond – Check condition number in case the matrix is a NumPy array.

• default_sparse_solver_backend – Default sparse solver backend to use (scipy, generic).

Returns:

|VectorArray| of the inverse operator evaluations.

Raises:

InversionError – The operator could not be inverted.

apply_inverse_adjoint(U, mu=None, initial_guess=None, least_squares=False)

Apply the inverse adjoint operator.

Parameters:

• U – VectorArray of vectors to which the inverse adjoint operator is applied.

• mu – The parameter values for which to evaluate the inverse adjoint operator.

• initial_guess – VectorArray with the same length as U containing initial guesses for the solution. Some implementations of apply_inverse_adjoint may ignore this parameter. If None a solver-dependent default is used.

• least_squares –

  If True, solve the least squares problem:

  v = argmin ||op^*(v) - u||_2.
  

  Since for an invertible operator the least squares solution agrees with the result of the application of the inverse operator, setting this option should, in general, have no effect on the result for those operators. However, note that when no appropriate solver_options are set for the operator, most operator implementations will choose a least squares solver by default which may be undesirable.

Returns:

|VectorArray| of the inverse adjoint operator evaluations.

Raises:

InversionError – The operator could not be inverted.

as_range_array(mu=None)

Return a VectorArray representation of the operator in its range space.

In the case of a linear operator with NumpyVectorSpace as source, this method returns for given parameter values mu a VectorArray V in the operator’s range, such that

V.lincomb(U.to_numpy()) == self.apply(U, mu)

for all VectorArrays U.

Parameters:

mu – The parameter values for which to return the VectorArray representation.

Returns:

V – The VectorArray defined above.

as_source_array(mu=None)

Return a VectorArray representation of the operator in its source space.

In the case of a linear operator with NumpyVectorSpace as range, this method returns for given parameter values mu a VectorArray V in the operator’s source, such that

self.range.make_array(V.inner(U)) == self.apply(U, mu)

for all VectorArrays U.

Parameters:

mu – The parameter values for which to return the VectorArray representation.

Returns:

V – The VectorArray defined above.

assemble(mu=None)

Assemble the operator for given parameter values.

The result of the method strongly depends on the given operator. For instance, a matrix-based operator will assemble its matrix, a LincombOperator will try to form the linear combination of its operators, whereas an arbitrary operator might simply return a FixedParameterOperator. The only assured property of the assembled operator is that it no longer depends on a Parameter.

Parameters:

mu – The parameter values for which to assemble the operator.

Returns:

Parameter-independent, assembled |Operator|.

classmethod from_file(path, key=None, solver_options=None, name=None)

class pymor.operators.numpy.NumpyToeplitzOperator(c, r=None, name=None)

Bases: pymor.operators.interface.Operator

Matrix-free representation of a finite dimensional Toeplitz matrix by a NumPy array.

A finite dimensional Toeplitz operator can be represented by a (block) matrix with constant diagonals (diagonal blocks), i.e.:

\[\begin{split}T = \begin{bmatrix} c_1 & r_2 & r_3 & \cdots & \cdots & r_n \\ c_2 & c_1 & r_2 & & & \vdots \\ c_3 & c_2 & \ddots & & & \vdots \\ \vdots & & & \ddots & r_2 & r_3 \\ \vdots & & & c_2 & c_1 & r_2 \\ c_m & \cdots & \cdots & c_3 & c_2 & c_1 \end{bmatrix} \in \mathbb{C}^{n*p \times k*m},\end{split}\]

where \(c\in\mathbb{C}^{n\times p\times m}\) and \(r\in\mathbb{C}^{k\times p\times m}\) denote the first (matrix-valued) column and first (matrix-valued) row of the (block) Toeplitz matrix, respectively. The matrix \(T\) as seen above is not explicitly constructed, only the arrays c and r are stored. The operator’s apply method takes advantage of the fact that any (block) Toeplitz matrix can be embedded in a larger (block) circulant matrix to leverage efficient matrix-vector multiplications with DFT.

Parameters:

• c – The NumPy array of shape either (n) or (n, p, m) that defines the first column of the (block) Toeplitz matrix.

• r – The NumPy array of shape (k,) or (k, p, m) that defines the first row of the Toeplitz matrix. If supplied, its first entry r[0] has to be equal to c[0]. Defaults to None. If r is None, the behavior of scipy.linalg.toeplitz is mimicked which sets r = c.conj() (except for the first entry).

• name – Name of the operator.

Methods

apply	Apply the operator to a VectorArray.
apply_adjoint	Apply the adjoint operator.

linear = True

apply(U, mu=None)

Apply the operator to a VectorArray.

Parameters:

• U – VectorArray of vectors to which the operator is applied.

• mu – The parameter values for which to evaluate the operator.

Returns:

|VectorArray| of the operator evaluations.

apply_adjoint(V, mu=None)

Apply the adjoint operator.

For any given linear Operator op, parameter values mu and VectorArrays U, V in the source resp. range we have:

op.apply_adjoint(V, mu).dot(U) == V.inner(op.apply(U, mu))

Thus, when op is represented by a matrix M, apply_adjoint is given by left-multiplication of (the complex conjugate of) M with V.

Parameters:

• V – VectorArray of vectors to which the adjoint operator is applied.

• mu – The parameter values for which to apply the adjoint operator.

Returns:

|VectorArray| of the adjoint operator evaluations.