pymor.operators.numpy
¶
Operators
based on NumPy
arrays.
This module provides the following NumPy
-based Operators
:
NumpyMatrixOperator
wraps a 2DNumPy array
as anOperator
.NumpyMatrixBasedOperator
should be used as base class for allOperators
which assemble into aNumpyMatrixOperator
.NumpyGenericOperator
wraps an arbitrary Python function betweenNumPy arrays
as anOperator
.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, source_id=None, range_id=None, name=None)[source]¶
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.- source_id
The id of the operator’s
source
VectorSpace
.- range_id
The id of the operator’s
range
VectorSpace
.- name
Name of the operator.
Methods
Apply the operator to a
VectorArray
.Apply the adjoint operator.
- apply(U, mu=None)[source]¶
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)[source]¶
Apply the adjoint operator.
For any given linear
Operator
op
,parameter values
mu
andVectorArrays
U
,V
in thesource
resp.range
we have:op.apply_adjoint(V, mu).dot(U) == V.inner(op.apply(U, mu))
Thus, when
op
is represented by a matrixM
,apply_adjoint
is given by left-multiplication of (the complex conjugate of)M
withV
.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={}, source_id=None, range_id=None, solver_options=None, name=None)[source]¶
Bases:
pymor.operators.interface.Operator
Wraps an arbitrary Python function between
NumPy arrays
as anOperator
.Parameters
- mapping
The function to wrap. If
parameters
isNone
, the function is of the formmapping(U)
and is expected to be vectorized. In particular:mapping(U).shape == U.shape[:-1] + (dim_range,).
If
parameters
is notNone
, the function has to have the signaturemapping(U, mu)
.- adjoint_mapping
The adjoint function to wrap. If
parameters
isNone
, the function is of the formadjoint_mapping(U)
and is expected to be vectorized. In particular:adjoint_mapping(U).shape == U.shape[:-1] + (dim_source,).
If
parameters
is notNone
, the function has to have the signatureadjoint_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 providedmapping
andadjoint_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 the operator to a
VectorArray
.Apply the adjoint operator.
- apply(U, mu=None)[source]¶
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)[source]¶
Apply the adjoint operator.
For any given linear
Operator
op
,parameter values
mu
andVectorArrays
U
,V
in thesource
resp.range
we have:op.apply_adjoint(V, mu).dot(U) == V.inner(op.apply(U, mu))
Thus, when
op
is represented by a matrixM
,apply_adjoint
is given by left-multiplication of (the complex conjugate of)M
withV
.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, source_id=None, range_id=None, name=None)[source]¶
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
andr
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 entryr[0]
has to be equal toc[-1]
. Defaults toNone
. Ifr
isNone
, the behavior ofscipy.linalg.hankel
is mimicked which setsr
to zero (except for the first entry).- source_id
The id of the operator’s
source
VectorSpace
.- range_id
The id of the operator’s
range
VectorSpace
.- name
Name of the operator.
Methods
Apply the operator to a
VectorArray
.Apply the adjoint operator.
- apply(U, mu=None)[source]¶
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)[source]¶
Apply the adjoint operator.
For any given linear
Operator
op
,parameter values
mu
andVectorArrays
U
,V
in thesource
resp.range
we have:op.apply_adjoint(V, mu).dot(U) == V.inner(op.apply(U, mu))
Thus, when
op
is represented by a matrixM
,apply_adjoint
is given by left-multiplication of (the complex conjugate of)M
withV
.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[source]¶
Bases:
pymor.operators.interface.Operator
Base class for operators which assemble into a
NumpyMatrixOperator
.- sparse[source]¶
True
if the operator assembles into a sparse matrix,False
if the operator assembles into a dense matrix,None
if unknown.
Methods
Apply the operator to a
VectorArray
.Apply the adjoint operator.
Apply the inverse operator.
Return a
VectorArray
representation of the operator in its range space.Return a
VectorArray
representation of the operator in its source space.Assemble the operator for given
parameter values
.Save the matrix of the operator to a file.
- apply(U, mu=None)[source]¶
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)[source]¶
Apply the adjoint operator.
For any given linear
Operator
op
,parameter values
mu
andVectorArrays
U
,V
in thesource
resp.range
we have:op.apply_adjoint(V, mu).dot(U) == V.inner(op.apply(U, mu))
Thus, when
op
is represented by a matrixM
,apply_adjoint
is given by left-multiplication of (the complex conjugate of)M
withV
.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)[source]¶
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 asV
containing initial guesses for the solution. Some implementations ofapply_inverse
may ignore this parameter. IfNone
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)[source]¶
Return a
VectorArray
representation of the operator in its range space.In the case of a linear operator with
NumpyVectorSpace
assource
, this method returns for givenparameter values
mu
aVectorArray
V
in the operator’srange
, such thatV.lincomb(U.to_numpy()) == self.apply(U, mu)
for all
VectorArrays
U
.Parameters
- mu
The
parameter values
for which to return theVectorArray
representation.
Returns
- V
The
VectorArray
defined above.
- as_source_array(mu=None)[source]¶
Return a
VectorArray
representation of the operator in its source space.In the case of a linear operator with
NumpyVectorSpace
asrange
, this method returns for givenparameter values
mu
aVectorArray
V
in the operator’ssource
, such thatself.range.make_array(V.inner(U).T) == self.apply(U, mu)
for all
VectorArrays
U
.Parameters
- mu
The
parameter values
for which to return theVectorArray
representation.
Returns
- V
The
VectorArray
defined above.
- assemble(mu=None)[source]¶
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 aFixedParameterOperator
. The only assured property of the assembled operator is that it no longer depends on aParameter
.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)[source]¶
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
, theOperator
’sname
is used.- output_format
Output file format. Either
matlab
ormatrixmarket
.- mu
The
parameter values
to assemble the to be exported matrix for.
- class pymor.operators.numpy.NumpyMatrixOperator(matrix, source_id=None, range_id=None, solver_options=None, name=None)[source]¶
Bases:
NumpyMatrixBasedOperator
Wraps a 2D
NumPy array
orSciPy spmatrix
as anOperator
.Note
In the case of a
NumPy array
, theapply_inverse
method by default usescheck_finite=True
andcheck_cond=True
. Setting them toFalse
(e.g., viadefaults
) can significantly speed up the computation, especially for smaller matrices.Parameters
- matrix
The
NumPy array
orSciPy spmatrix
which is to be wrapped.- source_id
The id of the operator’s
source
VectorSpace
.- range_id
The id of the operator’s
range
VectorSpace
.- solver_options
The
solver_options
for the operator.- name
Name of the operator.
Methods
Apply the operator to a
VectorArray
.Apply the adjoint operator.
Apply the inverse operator.
Apply the inverse adjoint operator.
Return a
VectorArray
representation of the operator in its range space.Return a
VectorArray
representation of the operator in its source space.Assemble the operator for given
parameter values
.- apply(U, mu=None)[source]¶
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)[source]¶
Apply the adjoint operator.
For any given linear
Operator
op
,parameter values
mu
andVectorArrays
U
,V
in thesource
resp.range
we have:op.apply_adjoint(V, mu).dot(U) == V.inner(op.apply(U, mu))
Thus, when
op
is represented by a matrixM
,apply_adjoint
is given by left-multiplication of (the complex conjugate of)M
withV
.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')[source]¶
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 asV
containing initial guesses for the solution. Some implementations ofapply_inverse
may ignore this parameter. IfNone
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)[source]¶
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 asU
containing initial guesses for the solution. Some implementations ofapply_inverse_adjoint
may ignore this parameter. IfNone
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)[source]¶
Return a
VectorArray
representation of the operator in its range space.In the case of a linear operator with
NumpyVectorSpace
assource
, this method returns for givenparameter values
mu
aVectorArray
V
in the operator’srange
, such thatV.lincomb(U.to_numpy()) == self.apply(U, mu)
for all
VectorArrays
U
.Parameters
- mu
The
parameter values
for which to return theVectorArray
representation.
Returns
- V
The
VectorArray
defined above.
- as_source_array(mu=None)[source]¶
Return a
VectorArray
representation of the operator in its source space.In the case of a linear operator with
NumpyVectorSpace
asrange
, this method returns for givenparameter values
mu
aVectorArray
V
in the operator’ssource
, such thatself.range.make_array(V.inner(U).T) == self.apply(U, mu)
for all
VectorArrays
U
.Parameters
- mu
The
parameter values
for which to return theVectorArray
representation.
Returns
- V
The
VectorArray
defined above.
- assemble(mu=None)[source]¶
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 aFixedParameterOperator
. The only assured property of the assembled operator is that it no longer depends on aParameter
.Parameters
- mu
The
parameter values
for which to assemble the operator.
Returns
Parameter-independent, assembled
Operator
.
- class pymor.operators.numpy.NumpyToeplitzOperator(c, r=None, source_id=None, range_id=None, name=None)[source]¶
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
andr
are stored. The operator’sapply
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 entryr[0]
has to be equal toc[0]
. Defaults toNone
. Ifr
isNone
, the behavior ofscipy.linalg.toeplitz
is mimicked which setsr = c.conj()
(except for the first entry).- source_id
The id of the operator’s
source
VectorSpace
.- range_id
The id of the operator’s
range
VectorSpace
.- name
Name of the operator.
Methods
Apply the operator to a
VectorArray
.Apply the adjoint operator.
- apply(U, mu=None)[source]¶
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)[source]¶
Apply the adjoint operator.
For any given linear
Operator
op
,parameter values
mu
andVectorArrays
U
,V
in thesource
resp.range
we have:op.apply_adjoint(V, mu).dot(U) == V.inner(op.apply(U, mu))
Thus, when
op
is represented by a matrixM
,apply_adjoint
is given by left-multiplication of (the complex conjugate of)M
withV
.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.