pymor.algorithms package¶
Submodules¶
adaptivegreedy module¶

class
pymor.algorithms.adaptivegreedy.
AdaptiveSampleSet
(parameter_space)[source]¶ Bases:
pymor.core.base.BasicObject
An adaptive parameter sample set.
Used by
adaptive_weak_greedy
.Methods
map_vertex_to_mu
,refine
,visualize
Attributes
Element

pymor.algorithms.adaptivegreedy.
adaptive_weak_greedy
(surrogate, parameter_space, target_error=None, max_extensions=None, validation_mus=0, rho=1.1, gamma=0.2, theta=0.0, visualize=False, visualize_vertex_size=80, pool=None)[source]¶ Weak greedy basis generation algorithm with adaptively refined training set.
This method extends pyMOR’s default
weak_greedy
greedy basis generation algorithm by adaptive refinement of the parameter training set according to [HDO11] to prevent overfitting of the approximation basis to the training set. This is achieved by estimating the approximation error on an additional validation set of parameters. If the ratio between the estimated errors on the validation set and the validation set is larger thanrho
, the training set is refined using standard grid refinement techniques.Parameters
 surrogate
See
weak_greedy
. parameter_space
The
ParameterSpace
for which to compute the approximation basis. target_error
See
weak_greedy
. max_extensions
See
weak_greedy
. validation_mus
 One of the following:
a list of
Parameters
to use as validation set,a positive number indicating the number of random parameters to use as validation set,
a nonpositive number, indicating the negative number of random parameters to use as validation set in addition to the centers of the elements of the adaptive training set.
 rho
Maximum allowed ratio between maximum estimated error on validation set vs. maximum estimated error on training set. If the ratio is larger, the training set is refined.
 gamma
Weight of the age penalty term in the training set refinement indicators.
 theta
Ratio of training set elements to select for refinement. (One element is always refined.)
 visualize
If
True
, visualize the refinement indicators. (Only available for 2 and 3 dimensional parameter spaces.) visualize_vertex_size
Size of the vertices in the visualization.
 pool
See
weak_greedy
.
Returns
Dict with the following fields
 extensions
Number of greedy iterations.
 max_errs
Sequence of maximum errors during the greedy run.
 max_err_mus
The parameters corresponding to
max_errs
. max_val_errs
Sequence of maximum errors on the validation set.
 max_val_err_mus
The parameters corresponding to
max_val_errs
. refinements
Number of refinements made in each extension step.
 training_set_sizes
The final size of the training set in each extension step.
 time
Duration of the algorithm.

pymor.algorithms.adaptivegreedy.
rb_adaptive_greedy
(fom, reductor, parameter_space, use_error_estimator=True, error_norm=None, target_error=None, max_extensions=None, validation_mus=0, rho=1.1, gamma=0.2, theta=0.0, extension_params=None, visualize=False, visualize_vertex_size=80, pool=None)[source]¶ Reduced basis greedy basis generation with adaptively refined training set.
This method extends pyMOR’s default
rb_greedy
greedy reduced basis generation algorithm by adaptive refinement of the parameter training set [HDO11] to prevent overfitting of the reduced basis to the training set as implemented inadaptive_weak_greedy
.Parameters
 fom
See
rb_greedy
. reductor
See
rb_greedy
. parameter_space
The
ParameterSpace
for which to compute the reduced model. use_error_estimator
See
rb_greedy
. error_norm
See
rb_greedy
. target_error
See
weak_greedy
. max_extensions
See
weak_greedy
. validation_mus
See
adaptive_weak_greedy
. rho
See
adaptive_weak_greedy
. gamma
See
adaptive_weak_greedy
. theta
See
adaptive_weak_greedy
. extension_params
See
rb_greedy
. visualize
See
adaptive_weak_greedy
. visualize_vertex_size
See
adaptive_weak_greedy
. pool
See
weak_greedy
.
Returns
Dict with the following fields
 rom
The reduced
Model
obtained for the computed basis. extensions
Number of greedy iterations.
 max_errs
Sequence of maximum errors during the greedy run.
 max_err_mus
The parameters corresponding to
max_errs
. max_val_errs
Sequence of maximum errors on the validation set.
 max_val_err_mus
The parameters corresponding to
max_val_errs
. refinements
Number of refinements made in each extension step.
 training_set_sizes
The final size of the training set in each extension step.
 time
Duration of the algorithm.
basic module¶
Module containing some basic but generic linear algebra algorithms.

pymor.algorithms.basic.
almost_equal
(U, V, product=None, sup_norm=False, rtol=1e14, atol=1e14)[source]¶ Compare U and V for almost equality.
The vectors of
U
andV
are compared in pairs for almost equality. Two vectorsu
andv
are considered almost equal iffu  v <= atol + v * rtol.
The norm to be used can be specified via the
product
orsup_norm
parameters.If the length of
U
resp.V
is 1, the single specified vector is compared to all vectors of the other array. Otherwise, the lengths of both indexed arrays have to agree.Parameters
 U, V
VectorArrays
to be compared. product
If specified, use this inner product
Operator
to compute the norm. Otherwise, the Euclidean norm is used. sup_norm
If
True
, use the sup norm to compute the error. rtol
The relative tolerance.
 atol
The absolute tolerance.
Defaults
rtol, atol (see
pymor.core.defaults
)

pymor.algorithms.basic.
contains_zero_vector
(vector_array, rtol=None, atol=None)[source]¶ returns
True
iff any vector in the array float_compares to 0s of the same dimParameters
 vector_array
a
VectorArray
implementation rtol
relative tolerance for float_cmp
 atol
absolute tolerance for float_cmp

pymor.algorithms.basic.
project_array
(U, basis, product=None, orthonormal=True)[source]¶ Orthogonal projection of
VectorArray
onto subspace.Parameters
 U
The
VectorArray
to project. basis
VectorArray
of basis vectors for the subspace onto which to project. product
Inner product
Operator
w.r.t. which to project. orthonormal
If
True
, the vectors inbasis
are assumed to be orthonormal w.r.t.product
.
Returns
The projected
VectorArray
.

pymor.algorithms.basic.
relative_error
(U, V, product=None)[source]¶ Compute error between U and V relative to norm of U.
ei module¶
This module contains algorithms for the empirical interpolation of Operators
.
The main work for generating the necessary interpolation data is handled by
the ei_greedy
method. The objects returned by this method can be used
to instantiate an EmpiricalInterpolatedOperator
.
As a convenience, the interpolate_operators
method allows to perform
the empirical interpolation of the Operators
of a given model with
a single function call.

pymor.algorithms.ei.
deim
(U, modes=None, pod=True, atol=None, rtol=None, product=None, pod_options={})[source]¶ Generate data for empirical interpolation using DEIM algorithm.
Given a
VectorArray
U
, this method generates a collateral basis and interpolation DOFs for empirical interpolation of the vectors contained inU
. The returned objects can be used to instantiate anEmpiricalInterpolatedOperator
(withtriangular=False
).The collateral basis is determined by the first
pod
modes ofU
.Parameters
 U
A
VectorArray
of vectors to interpolate. modes
Dimension of the collateral basis i.e. number of POD modes of the vectors in
U
. pod
If
True
, perform a POD ofU
to obtain the collateral basis. IfFalse
,U
is used as collateral basis. atol
Absolute POD tolerance.
 rtol
Relative POD tolerance.
 product
Inner product
Operator
used for the POD. pod_options
Dictionary of additional options to pass to the
pod
algorithm.
Returns
 interpolation_dofs
NumPy array
of the DOFs at which the vectors are interpolated. collateral_basis
VectorArray
containing the generated collateral basis. data
Dict containing the following fields:
 svals
POD singular values.

pymor.algorithms.ei.
ei_greedy
(U, error_norm=None, atol=None, rtol=None, max_interpolation_dofs=None, copy=True, pool=DummyPool('??', '??'))[source]¶ Generate data for empirical interpolation using EIGreedy algorithm.
Given a
VectorArray
U
, this method generates a collateral basis and interpolation DOFs for empirical interpolation of the vectors contained inU
. The returned objects can be used to instantiate anEmpiricalInterpolatedOperator
(withtriangular=True
).The interpolation data is generated by a greedy search algorithm, where in each loop iteration the worst approximated vector in
U
is added to the collateral basis.Parameters
 U
A
VectorArray
of vectors to interpolate. error_norm
Norm w.r.t. which to calculate the interpolation error. If
None
, the Euclidean norm is used. atol
Stop the greedy search if the largest approximation error is below this threshold.
 rtol
Stop the greedy search if the largest relative approximation error is below this threshold.
 max_interpolation_dofs
Stop the greedy search if the number of interpolation DOF (= dimension of the collateral basis) reaches this value.
 copy
If
False
,U
will be modified during executing of the algorithm. pool
If not
None
, theWorkerPool
to use for parallelization.
Returns
 interpolation_dofs
NumPy array
of the DOFs at which the vectors are evaluated. collateral_basis
VectorArray
containing the generated collateral basis. data
Dict containing the following fields:
 errors
Sequence of maximum approximation errors during greedy search.
 triangularity_errors
Sequence of maximum absolute values of interoplation matrix coefficients in the upper triangle (should be near zero).

pymor.algorithms.ei.
interpolate_operators
(fom, operator_names, parameter_sample, error_norm=None, product=None, atol=None, rtol=None, max_interpolation_dofs=None, pod_options={}, alg='ei_greedy', pool=DummyPool('??', '??'))[source]¶ Empirical operator interpolation using the EIGreedy/DEIM algorithm.
This is a convenience method to facilitate the use of
ei_greedy
ordeim
. Given aModel
, names ofOperators
, and a sample ofParameters
, first the operators are evaluated on the solution snapshots of the model for the provided parameters. These evaluations are then used as input forei_greedy
/deim
. Finally the resulting interpolation data is used to createEmpiricalInterpolatedOperators
and a new model with the interpolated operators is returned.Note that this implementation creates one common collateral basis for all specified operators, which might not be what you want.
Parameters
 fom
 operator_names
List of keys in the
operators
dict of the model. The correspondingOperators
will be interpolated. parameter_sample
A list of
Parameters
for which solution snapshots are calculated. error_norm
See
ei_greedy
. Has no effect ifalg == 'deim'
. product
Inner product for POD computation in
deim
. Has no effect ifalg == 'ei_greedy'
. atol
See
ei_greedy
. rtol
See
ei_greedy
. max_interpolation_dofs
See
ei_greedy
. pod_options
Further options for
pod
algorithm. Has no effect ifalg == 'ei_greedy'
. alg
Either
ei_greedy
ordeim
. pool
If not
None
, theWorkerPool
to use for parallelization.
Returns
 eim
Model
withOperators
given byoperator_names
replaced byEmpiricalInterpolatedOperators
. data
Dict containing the following fields:
 dofs
NumPy array
of the DOFs at which theOperators
have to be evaluated. basis
VectorArray
containing the generated collateral basis.
In addition,
data
contains the fields of thedata
dict
returned byei_greedy
/deim
.
eigs module¶

pymor.algorithms.eigs.
_extend_arnoldi
(A, E, V, H, f, p)[source]¶ Extend an existing Arnoldi factorization.

pymor.algorithms.eigs.
eigs
(A, E=None, k=3, which='LM', b=None, l=None, maxiter=1000, tol=1e13, imag_tol=1e12, complex_pair_tol=1e12, seed=0)[source]¶ Approximate a few eigenvalues of a linear
Operator
.Computes
k
eigenvaluesw
with corresponding eigenvectorsv
which solve the eigenvalue problem\[A v_i = w_i v_i\]or the generalized eigenvalue problem
\[A v_i = w_i E v_i\]if
E
is notNone
.The implementation is based on Algorithm 4.2 in [RL95].
Parameters
 A
The real linear
Operator
for which the eigenvalues are to be computed. E
The real linear
Operator
which defines the generalized eigenvalue problem. k
The number of eigenvalues and eigenvectors which are to be computed.
 which
A string specifying which
k
eigenvalues and eigenvectors to compute:'LM'
: select eigenvalues with largest magnitude'SM'
: select eigenvalues with smallest magnitude'LR'
: select eigenvalues with largest real part'SR'
: select eigenvalues with smallest real part'LI'
: select eigenvalues with largest imaginary part'SI'
: select eigenvalues with smallest imaginary part
 b
Initial vector for Arnoldi iteration. Default is a random vector.
 l
The size of the Arnoldi factorization. Default is
min(n  1, max(2*k + 1, 20))
. maxiter
The maximum number of iterations.
 tol
The relative error tolerance for the Ritz estimates.
 imag_tol
Relative imaginary parts below this tolerance are set to 0.
 complex_pair_tol
Tolerance for detecting pairs of complex conjugate eigenvalues.
 seed
Random seed which is used for computing the initial vector for the Arnoldi iteration.
Returns
 w
A 1D
NumPy array
which contains the computed eigenvalues. v
A
VectorArray
which contains the computed eigenvectors.
error module¶

pymor.algorithms.error.
reduction_error_analysis
(rom, fom, reductor, test_mus, basis_sizes=0, error_estimator=True, condition=False, error_norms=(), error_norm_names=None, error_estimator_norm_index=0, custom=(), plot=False, plot_custom_logarithmic=True, pool=DummyPool('??', '??'))[source]¶ Analyze the model reduction error.
The maximum model reduction error is estimated by solving the reduced
Model
for given randomParameters
.Parameters
 rom
The reduced
Model
. fom
The highdimensional
Model
. reductor
The reductor which has created
rom
. test_mus
List of
Parameters
to compute the errors for. basis_sizes
Either a list of reduced basis dimensions to consider, or the number of dimensions (which are then selected equidistantly, always including the maximum reduced space dimension). The dimensions are input for the
dim
Parameter ofreductor.reduce()
. error_estimator
If
True
evaluate the error estimator ofrom
on the testParameters
. condition
If
True
, compute the condition of the reduced system matrix for the given testParameters
(can only be specified ifrom
is an instance ofStationaryModel
androm.operator
is linear). error_norms
List of norms in which to compute the model reduction error.
 error_norm_names
Names of the norms given by
error_norms
. IfNone
, thename
attributes of the given norms are used. error_estimator_norm_index
When
error_estimator
isTrue
anderror_norms
are specified, this is the index of the norm inerror_norms
w.r.t. which to compute the effectivity of the error estimator. custom
List of custom functions which are evaluated for each test
parameter values
and basis size. The functions must have the signaturedef custom_value(rom, fom, reductor, mu, dim): pass
 plot
If
True
, generate a plot of the computed quantities w.r.t. the basis size. plot_custom_logarithmic
If
True
, use a logarithmic yaxis to plot the computed custom values. pool
If not
None
, theWorkerPool
to use for parallelization.
Returns
Dict with the following fields
 mus
The test
Parameters
which have been considered. basis_sizes
The reduced basis dimensions which have been considered.
 norms
NumPy array
of the norms of the highdimensional solutions w.r.t. all given testParameters
, reduced basis dimensions and norms inerror_norms
. (Only present whenerror_norms
has been specified.) max_norms
Maxima of
norms
over the given testParameters
. max_norm_mus
Parameters
corresponding tomax_norms
. errors
NumPy array
of the norms of the model reduction errors w.r.t. all given testParameters
, reduced basis dimensions and norms inerror_norms
. (Only present whenerror_norms
has been specified.) max_errors
Maxima of
errors
over the given testParameters
. max_error_mus
Parameters
corresponding tomax_errors
. rel_errors
errors
divided bynorms
. (Only present whenerror_norms
has been specified.) max_rel_errors
Maxima of
rel_errors
over the given testParameters
. max_rel_error_mus
Parameters
corresponding tomax_rel_errors
. error_norm_names
Names of the given
error_norms
. (Only present whenerror_norms
has been specified.) error_estimates
NumPy array
of the model reduction error estimates w.r.t. all given testParameters
and reduced basis dimensions. (Only present whenerror_estimator
isTrue
.) max_error_estimate
Maxima of
error_estimates
over the given testParameters
. max_error_estimate_mus
Parameters
corresponding tomax_error_estimates
. effectivities
errors
divided byerror_estimates
. (Only present whenerror_estimator
isTrue
anderror_norms
has been specified.) min_effectivities
Minima of
effectivities
over the given testParameters
. min_effectivity_mus
Parameters
corresponding tomin_effectivities
. max_effectivities
Maxima of
effectivities
over the given testParameters
. max_effectivity_mus
Parameters
corresponding tomax_effectivities
. errors
NumPy array
of the reduced system matrix conditions w.r.t. all given testParameters
and reduced basis dimensions. (Only present whenconditions
isTrue
.) max_conditions
Maxima of
conditions
over the given testParameters
. max_condition_mus
Parameters
corresponding tomax_conditions
. custom_values
NumPy array
of custom function evaluations w.r.t. all given testParameters
, reduced basis dimensions and functions incustom
. (Only present whencustom
has been specified.) max_custom_values
Maxima of
custom_values
over the given testParameters
. max_custom_values_mus
Parameters
corresponding tomax_custom_values
. time
Time (in seconds) needed for the error analysis.
 summary
String containing a summary of all computed quantities for the largest (last) considered basis size.
 figure
The figure containing the generated plots. (Only present when
plot
isTrue
.)
genericsolvers module¶
This module contains some iterative linear solvers which only use the Operator
interface

pymor.algorithms.genericsolvers.
apply_inverse
(op, V, initial_guess=None, options=None, least_squares=False, check_finite=True, default_solver='generic_lgmres', default_least_squares_solver='generic_least_squares_lsmr')[source]¶ Solve linear equation system.
Applies the inverse of
op
to the vectors inV
using a generic iterative solver.Parameters
 op
The linear, nonparametric
Operator
to invert. V
VectorArray
of righthand sides for the equation system. initial_guess
VectorArray
with the same length asV
containing initial guesses for the solution. Some implementations ofapply_inverse
may ignore this parameter. IfNone
a solverdependent default is used. options
The
solver_options
to use (seesolver_options
). least_squares
If
True
, return least squares solution. check_finite
Test if solution only contains finite values.
 default_solver
Default solver to use (generic_lgmres, generic_least_squares_lsmr, generic_least_squares_lsqr).
 default_least_squares_solver
Default solver to use for least squares problems (generic_least_squares_lsmr, generic_least_squares_lsqr).
Returns
VectorArray
of the solution vectors.Defaults
check_finite, default_solver, default_least_squares_solver (see
pymor.core.defaults
)

pymor.algorithms.genericsolvers.
lgmres
(A, b, x0=None, tol=1e05, maxiter=1000, M=None, callback=None, inner_m=30, outer_k=3, outer_v=None, store_outer_Av=True)[source]¶

pymor.algorithms.genericsolvers.
lsmr
(A, b, damp=0.0, atol=1e06, btol=1e06, conlim=100000000.0, maxiter=None, show=False)[source]¶

pymor.algorithms.genericsolvers.
lsqr
(A, b, damp=0.0, atol=1e08, btol=1e08, conlim=100000000.0, iter_lim=None, show=False)[source]¶

pymor.algorithms.genericsolvers.
solver_options
(lgmres_tol=1e05, lgmres_maxiter=1000, lgmres_inner_m=39, lgmres_outer_k=3, least_squares_lsmr_damp=0.0, least_squares_lsmr_atol=1e06, least_squares_lsmr_btol=1e06, least_squares_lsmr_conlim=100000000.0, least_squares_lsmr_maxiter=None, least_squares_lsmr_show=False, least_squares_lsqr_damp=0.0, least_squares_lsqr_atol=1e06, least_squares_lsqr_btol=1e06, least_squares_lsqr_conlim=100000000.0, least_squares_lsqr_iter_lim=None, least_squares_lsqr_show=False)[source]¶ Returns available solvers with default
solver_options
.Parameters
 lgmres_tol
 lgmres_maxiter
 lgmres_inner_m
 lgmres_outer_k
 least_squares_lsmr_damp
 least_squares_lsmr_atol
 least_squares_lsmr_btol
 least_squares_lsmr_conlim
 least_squares_lsmr_maxiter
 least_squares_lsmr_show
 least_squares_lsqr_damp
 least_squares_lsqr_atol
 least_squares_lsqr_btol
 least_squares_lsqr_conlim
 least_squares_lsqr_iter_lim
 least_squares_lsqr_show
Returns
A dict of available solvers with default
solver_options
.Defaults
lgmres_tol, lgmres_maxiter, lgmres_inner_m, lgmres_outer_k, least_squares_lsmr_damp, least_squares_lsmr_atol, least_squares_lsmr_btol, least_squares_lsmr_conlim, least_squares_lsmr_maxiter, least_squares_lsmr_show, least_squares_lsqr_atol, least_squares_lsqr_btol, least_squares_lsqr_conlim, least_squares_lsqr_iter_lim, least_squares_lsqr_show (see
pymor.core.defaults
)
gram_schmidt module¶

pymor.algorithms.gram_schmidt.
gram_schmidt
(A, product=None, return_R=False, atol=1e13, rtol=1e13, offset=0, reiterate=True, reiteration_threshold=0.9, check=True, check_tol=0.001, copy=True)[source]¶ Orthonormalize a
VectorArray
using the modified GramSchmidt algorithm.Parameters
 A
The
VectorArray
which is to be orthonormalized. product
The inner product
Operator
w.r.t. which to orthonormalize. IfNone
, the Euclidean product is used. return_R
If
True
, the R matrix from QR decomposition is returned. atol
Vectors of norm smaller than
atol
are removed from the array. rtol
Relative tolerance used to detect linear dependent vectors (which are then removed from the array).
 offset
Assume that the first
offset
vectors are already orthonormal and start the algorithm at theoffset + 1
th vector. reiterate
If
True
, orthonormalize again if the norm of the orthogonalized vector is much smaller than the norm of the original vector. reiteration_threshold
If
reiterate
isTrue
, reorthonormalize if the ratio between the norms of the orthogonalized vector and the original vector is smaller than this value. check
If
True
, check if the resultingVectorArray
is really orthonormal. check_tol
Tolerance for the check.
 copy
If
True
, create a copy ofA
instead of modifyingA
inplace.
Returns
 Q
The orthonormalized
VectorArray
. R
The uppertriangular/trapezoidal matrix (if
compute_R
isTrue
).
Defaults
atol, rtol, reiterate, reiteration_threshold, check, check_tol (see
pymor.core.defaults
)

pymor.algorithms.gram_schmidt.
gram_schmidt_biorth
(V, W, product=None, reiterate=True, reiteration_threshold=0.1, check=True, check_tol=0.001, copy=True)[source]¶ Biorthonormalize a pair of
VectorArrays
using the biorthonormal GramSchmidt process.See Algorithm 1 in [BKS11].
Parameters
 V, W
The
VectorArrays
which are to be biorthonormalized. product
The inner product
Operator
w.r.t. which to biorthonormalize. IfNone
, the Euclidean product is used. reiterate
If
True
, orthonormalize again if the norm of the orthogonalized vector is much smaller than the norm of the original vector. reiteration_threshold
If
reiterate
isTrue
, reorthonormalize if the ratio between the norms of the orthogonalized vector and the original vector is smaller than this value. check
If
True
, check if the resultingVectorArray
is really orthonormal. check_tol
Tolerance for the check.
 copy
If
True
, create a copy ofV
andW
instead of modifyingV
andW
inplace.
Returns
The biorthonormalized
VectorArrays
.
greedy module¶

class
pymor.algorithms.greedy.
RBSurrogate
(fom, reductor, use_error_estimator, error_norm, extension_params, pool)[source]¶ Bases:
pymor.algorithms.greedy.WeakGreedySurrogate
Surrogate for the
weak_greedy
error used inrb_greedy
.Not intended to be used directly.
Methods
evaluate
,extend
Attributes

evaluate
(mus, return_all_values=False)[source]¶ Evaluate the surrogate for given parameters.
Parameters
 mus
List of parameters for which to estimate the approximation error. When parallelization is used,
mus
can be aRemoteObject
. return_all_values
See below.
Returns
If
return_all_values
isTrue
, anNumPy array
of the estimated errors. Ifreturn_all_values
isFalse
, the maximum estimated error as first return value and the corresponding parameter as second return value.


class
pymor.algorithms.greedy.
WeakGreedySurrogate
[source]¶ Bases:
pymor.core.base.BasicObject
Surrogate for the approximation error in
weak_greedy
.Methods
evaluate
,extend
Attributes

abstract
evaluate
(mus, return_all_values=False)[source]¶ Evaluate the surrogate for given parameters.
Parameters
 mus
List of parameters for which to estimate the approximation error. When parallelization is used,
mus
can be aRemoteObject
. return_all_values
See below.
Returns
If
return_all_values
isTrue
, anNumPy array
of the estimated errors. Ifreturn_all_values
isFalse
, the maximum estimated error as first return value and the corresponding parameter as second return value.

abstract

pymor.algorithms.greedy.
rb_greedy
(fom, reductor, training_set, use_error_estimator=True, error_norm=None, atol=None, rtol=None, max_extensions=None, extension_params=None, pool=None)[source]¶ Weak Greedy basis generation using the RB approximation error as surrogate.
This algorithm generates a reduced basis using the
weak greedy
algorithm [BCDDPW11], where the approximation error is estimated from computing solutions of the reduced order model for the current reduced basis and then estimating the model reduction error.Parameters
 fom
The
Model
to reduce. reductor
Reductor for reducing the given
Model
. This has to be an object with areduce
method, such thatreductor.reduce()
yields the reduced model, and anexted_basis
method, such thatreductor.extend_basis(U, copy_U=False, **extension_params)
extends the current reduced basis by the vectors contained inU
. For an example seeCoerciveRBReductor
. training_set
The training set of
Parameters
on which to perform the greedy search. use_error_estimator
If
False
, exactly compute the model reduction error by also computing the solution offom
for allparameter values
of the training set. This is mainly useful when no estimator for the model reduction error is available. error_norm
If
use_error_estimator
isFalse
, use this function to calculate the norm of the error. IfNone
, the Euclidean norm is used. atol
See
weak_greedy
. rtol
See
weak_greedy
. max_extensions
See
weak_greedy
. extension_params
dict
of parameters passed to thereductor.extend_basis
method. IfNone
,'gram_schmidt'
basis extension will be used as a default for stationary problems (fom.solve
returnsVectorArrays
of length 1) and'pod'
basis extension (adding a single POD mode) for instationary problems. pool
See
weak_greedy
.
Returns
Dict with the following fields
 rom
The reduced
Model
obtained for the computed basis. max_errs
Sequence of maximum errors during the greedy run.
 max_err_mus
The parameters corresponding to
max_errs
. extensions
Number of performed basis extensions.
 time
Total runtime of the algorithm.

pymor.algorithms.greedy.
weak_greedy
(surrogate, training_set, atol=None, rtol=None, max_extensions=None, pool=None)[source]¶ Weak greedy basis generation algorithm [BCDDPW11].
This algorithm generates an approximation basis for a given set of vectors associated with a training set of parameters by iteratively evaluating a
surrogate
for the approximation error on the training set and adding the worst approximated vector (according to the surrogate) to the basis.The constructed basis is extracted from the surrogate after termination of the algorithm.
Parameters
 surrogate
An instance of
WeakGreedySurrogate
representing the surrogate for the approximation error. training_set
The set of parameter samples on which to perform the greedy search.
 atol
If not
None
, stop the algorithm if the maximum (estimated) error on the training set drops below this value. rtol
If not
None
, stop the algorithm if the maximum (estimated) relative error on the training set drops below this value. max_extensions
If not
None
, stop the algorithm aftermax_extensions
extension steps. pool
If not
None
, aWorkerPool
to use for parallelization. Parallelization needs to be supported bysurrogate
.
Returns
Dict with the following fields
 max_errs
Sequence of maximum estimated errors during the greedy run.
 max_err_mus
The parameters corresponding to
max_errs
. extensions
Number of performed basis extensions.
 time
Total runtime of the algorithm.
hapod module¶

class
pymor.algorithms.hapod.
DistHAPODTree
(slices)[source]¶ Bases:
pymor.algorithms.hapod.Tree
Methods
children
is_leaf
Attributes
depth
,root

class
pymor.algorithms.hapod.
IncHAPODTree
(steps)[source]¶ Bases:
pymor.algorithms.hapod.Tree
Methods
children
is_leaf
Attributes
depth
,root

class
pymor.algorithms.hapod.
LifoExecutor
(executor, max_workers=None)[source]¶ Bases:
object
Methods
done_callback
,run_task
,submit

class
pymor.algorithms.hapod.
Tree
[source]¶ Bases:
pymor.core.base.BasicObject
A rooted tree.
Methods
children
,is_leaf
Attributes
depth
,root

pymor.algorithms.hapod.
dist_hapod
(num_slices, snapshots, eps, omega, product=None, executor=None, eval_snapshots_in_executor=False)[source]¶ Distributed Hierarchical Approximate POD.
This computes the distributed HAPOD from [HLR18].
Parameters
 num_slices
The number of snapshot vector slices.
 snapshots
A mapping
snapshots(slice)
returning for each slice number the associated snapshot vectors. eps
Desired l2mean approximation error.
 omega
Tuning parameter (0 < omega < 1) to balance performance with approximation quality.
 product
Inner product
Operator
w.r.t. which to compute the POD. executor
If not
None
, aconcurrent.futures.Executor
object to use for parallelization. eval_snapshots_in_executor
If
True
also parallelize the evaluation of the snapshot map.
Returns
 modes
The computed POD modes.
 svals
The associated singular values.
 snap_count
The total number of input snapshot vectors.

pymor.algorithms.hapod.
dist_vectorarray_hapod
(num_slices, U, eps, omega, product=None, executor=None)[source]¶ Distributed Hierarchical Approximate POD.
This computes the distributed HAPOD from [HLR18] of a given
VectorArray
.Parameters
 num_slices
The number of snapshot vector slices.
 U
The
VectorArray
of which to compute the HAPOD. eps
Desired l2mean approximation error.
 omega
Tuning parameter (0 < omega < 1) to balance performance with approximation quality.
 product
Inner product
Operator
w.r.t. which to compute the POD. executor
If not
None
, aconcurrent.futures.Executor
object to use for parallelization.
Returns
 modes
The computed POD modes.
 svals
The associated singular values.
 snap_count
The total number of input snapshot vectors.

pymor.algorithms.hapod.
hapod
(tree, snapshots, local_eps, product=None, pod_method=<function default_pod_method>, executor=None, eval_snapshots_in_executor=False)[source]¶ Compute the Hierarchical Approximate POD.
This is an implementation of the HAPOD algorithm from [HLR18].
Parameters
 tree
A
Tree
defining the worker topology. snapshots
A mapping
snapshots(node)
returning for each leaf node the associated snapshot vectors. local_eps
A mapping
local_eps(node, snap_count, num_vecs)
assigning to each tree nodenode
an l2 truncation error tolerance for the local pod based on the number of input vectorsnum_vecs
and the total number of snapshot vectors below the given nodesnap_count
. product
Inner product
Operator
w.r.t. which to compute the POD. pod_method
A function
pod_method(U, eps, root_node, product)
for computing the POD of theVectorArray
U
w.r.t. the given inner productproduct
and the l2 error toleranceeps
.root_node
is set toTrue
when the POD is computed at the root of the tree. executor
If not
None
, aconcurrent.futures.Executor
object to use for parallelization. eval_snapshots_in_executor
If
True
also parallelize the evaluation of the snapshot map.
Returns
 modes
The computed POD modes.
 svals
The associated singular values.
 snap_count
The total number of input snapshot vectors.

pymor.algorithms.hapod.
inc_hapod
(steps, snapshots, eps, omega, product=None, executor=None, eval_snapshots_in_executor=False)[source]¶ Incremental Hierarchical Approximate POD.
This computes the incremental HAPOD from [HLR18].
Parameters
 steps
The number of incremental POD updates.
 snapshots
A mapping
snapshots(step)
returning for each incremental POD step the associated snapshot vectors. eps
Desired l2mean approximation error.
 omega
Tuning parameter (0 < omega < 1) to balance performance with approximation quality.
 product
Inner product
Operator
w.r.t. which to compute the POD. executor
If not
None
, aconcurrent.futures.Executor
object to use for parallelization. eval_snapshots_in_executor
If
True
also parallelize the evaluation of the snapshot map.
Returns
 modes
The computed POD modes.
 svals
The associated singular values.
 snap_count
The total number of input snapshot vectors.

pymor.algorithms.hapod.
inc_vectorarray_hapod
(steps, U, eps, omega, product=None, executor=None)[source]¶ Incremental Hierarchical Approximate POD.
This computes the incremental HAPOD from [HLR18] for a given
VectorArray
.Parameters
 steps
The number of incremental POD updates.
 U
The
VectorArray
of which to compute the HAPOD. eps
Desired l2mean approximation error.
 omega
Tuning parameter (0 < omega < 1) to balance performance with approximation quality.
 product
Inner product
Operator
w.r.t. which to compute the POD. executor
If not
None
, aconcurrent.futures.Executor
object to use for parallelization. eval_snapshots_in_executor
If
True
also parallelize the evaluation of the snapshot map.
Returns
 modes
The computed POD modes.
 svals
The associated singular values.
 snap_count
The total number of input snapshot vectors.
image module¶

class
pymor.algorithms.image.
CollectOperatorRangeRules
(source, image, extends)[source]¶ Bases:
pymor.algorithms.rules.RuleTable
RuleTable
for theestimate_image
algorithm.Methods
append_rule
,apply
,apply_children
,breakpoint_for_name
,breakpoint_for_obj
,get_children
,insert_rule
,replace_children
Attributes
action_apply_operator
,action_ConcatenationOperator
,action_EmpiricalInterpolatedOperator
,action_recurse
,rules

class
pymor.algorithms.image.
CollectVectorRangeRules
(image)[source]¶ Bases:
pymor.algorithms.rules.RuleTable
RuleTable
for theestimate_image
algorithm.Methods
append_rule
,apply
,apply_children
,breakpoint_for_name
,breakpoint_for_obj
,get_children
,insert_rule
,replace_children
Attributes
action_as_range_array
,action_ConcatenationOperator
,action_recurse
,action_VectorArray
,rules

pymor.algorithms.image.
estimate_image
(operators=(), vectors=(), domain=None, extends=False, orthonormalize=True, product=None, riesz_representatives=False)[source]¶ Estimate the image of given
Operators
for all mu.Let
operators
be a list ofOperators
with common source and range, and letvectors
be a list ofVectorArrays
or vectorlikeOperators
in the range of these operators. Given aVectorArray
domain
of vectors in the source of the operators, this algorithms determines aVectorArray
image
of range vectors such that the linear span ofimage
contains:op.apply(U, mu=mu)
for all operatorsop
inoperators
, for all possibleParameters
mu
and for allVectorArrays
U
contained in the linear span ofdomain
,U
for allVectorArrays
invectors
,v.as_range_array(mu)
for allOperators
invectors
and all possibleParameters
mu
.
The algorithm will try to choose
image
as small as possible. However, no optimality is guaranteed. The image estimation algorithm is specified byCollectOperatorRangeRules
andCollectVectorRangeRules
.Parameters
 operators
See above.
 vectors
See above.
 domain
See above. If
None
, an emptydomain
VectorArray
is assumed. extends
For some operators, e.g.
EmpiricalInterpolatedOperator
, as well as for all elements ofvectors
,image
is estimated independently from the choice ofdomain
. Ifextends
isTrue
, such operators are ignored. (This is useful in case these vectors have already been obtained by earlier calls to this function.) orthonormalize
Compute an orthonormal basis for the linear span of
image
using thegram_schmidt
algorithm. product
Inner product
Operator
w.r.t. which to orthonormalize. riesz_representatives
If
True
, compute Riesz representatives of the vectors inimage
before orthonormalizing (useful for dual norm computation when the range of theoperators
is a dual space).
Returns
The
VectorArray
image
.Raises
 ImageCollectionError
Is raised when for a given
Operator
no image estimate is possible.

pymor.algorithms.image.
estimate_image_hierarchical
(operators=(), vectors=(), domain=None, extends=None, orthonormalize=True, product=None, riesz_representatives=False)[source]¶ Estimate the image of given
Operators
for all mu.This is an extended version of
estimate_image
, which callsestimate_image
individually for each vector ofdomain
.As a result, the vectors in the returned
image
VectorArray
will be ordered by thedomain
vector they correspond to (starting with vectors which correspond to the elements ofvectors
and toOperators
for which the image is estimated independently fromdomain
).This function also returns an
image_dims
list, such that the firstimage_dims[i+1]
vectors ofimage
correspond to the firsti
vectors ofdomain
(the firstimage_dims[0]
vectors correspond tovectors
and toOperators
with fixed image estimate).Parameters
 operators
See
estimate_image
. vectors
See
estimate_image
. domain
See
estimate_image
. extends
When additional vectors have been appended to the
domain
VectorArray
afterestimate_image_hierarchical
has been called, andestimate_image_hierarchical
shall be called again for the extendeddomain
array,extends
can be set to(image, image_dims)
, whereimage
,image_dims
are the return values of the lastestimate_image_hierarchical
call. The olddomain
vectors will then be skipped during computation andimage
,image_dims
will be modified inplace. orthonormalize
See
estimate_image
. product
See
estimate_image
. riesz_representatives
See
estimate_image
.
Returns
 image
See above.
 image_dims
See above.
Raises
 ImageCollectionError
Is raised when for a given
Operator
no image estimate is possible.
krylov module¶
Module for computing (rational) Krylov subspaces’ bases.

pymor.algorithms.krylov.
rational_arnoldi
(A, E, b, sigma, trans=False)[source]¶ Rational Arnoldi algorithm.
If
trans == False
, using Arnoldi process, computes a real orthonormal basis for the rational Krylov subspace\[\mathrm{span}\{ (\sigma_1 E  A)^{1} b, (\sigma_2 E  A)^{1} b, \ldots, (\sigma_r E  A)^{1} b \},\]otherwise, computes the same for
\[\mathrm{span}\{ (\sigma_1 E  A)^{T} b^T, (\sigma_2 E  A)^{T} b^T, \ldots, (\sigma_r E  A)^{T} b^T \}.\]Interpolation points in
sigma
are allowed to repeat (in any order). Then, in the above expression,\[\underbrace{ (\sigma_i E  A)^{1} b, \ldots, (\sigma_i E  A)^{1} b }_{m \text{ times}}\]is replaced by
\[(\sigma_i E  A)^{1} b, (\sigma_i E  A)^{1} E (\sigma_i E  A)^{1} b, \ldots, \left((\sigma_i E  A)^{1} E\right)^{m  1} (\sigma_i E  A)^{1} b.\]Analogously for the
trans == True
case.Parameters
Returns
 V
Orthonormal basis for the Krylov subspace
VectorArray
.

pymor.algorithms.krylov.
tangential_rational_krylov
(A, E, B, b, sigma, trans=False, orth=True)[source]¶ Tangential Rational Krylov subspace.
If
trans == False
, computes a real basis for the rational Krylov subspace\[\mathrm{span}\{ (\sigma_1 E  A)^{1} B b_1, (\sigma_2 E  A)^{1} B b_2, \ldots, (\sigma_r E  A)^{1} B b_r \},\]otherwise, computes the same for
\[\mathrm{span}\{ (\sigma_1 E  A)^{T} B^T b_1, (\sigma_2 E  A)^{T} B^T b_2, \ldots, (\sigma_r E  A)^{T} B^T b_r \}.\]Interpolation points in
sigma
are assumed to be pairwise distinct.Parameters
 A
Real
Operator
A. E
Real
Operator
E. B
Real
Operator
B. b
VectorArray
fromB.source
, iftrans == False
, orB.range
, iftrans == True
.
 sigma
Sequence of interpolation points (closed under conjugation), of the same length as
b
. trans
Boolean, see above.
 orth
If
True
, orthonormalizes the basis usingpymor.algorithms.gram_schmidt.gram_schmidt
.
Returns
 V
Optionally orthonormal basis for the Krylov subspace
VectorArray
.
lincomb module¶

class
pymor.algorithms.lincomb.
AssembleLincombRules
(coefficients, solver_options, name)[source]¶ Bases:
pymor.algorithms.rules.RuleTable
Methods
append_rule
,apply
,apply_children
,breakpoint_for_name
,breakpoint_for_obj
,get_children
,insert_rule
,replace_children
Attributes
action_BlockDiagonalOperator
,action_BlockOperatorBase
,action_BlockSpaceIdentityOperator
,action_call_assemble_lincomb_method
,action_failed
,action_IdentityAndSecondOrderModelOperator
,action_IdentityOperator
,action_merge_into_low_rank_updated_operator
,action_merge_low_rank_operators
,action_only_one_operator
,action_VectorArrayOperator
,action_zero_coeff
,action_ZeroOperator
,rules

pymor.algorithms.lincomb.
assemble_lincomb
(operators, coefficients, solver_options=None, name=None)[source]¶ Try to assemble a linear combination of the given operators.
Returns a new
Operator
which represents the sumc_1*O_1 + ... + c_N*O_N
where
O_i
areOperators
andc_i
scalar coefficients.This function is called in the
assemble
method ofLincombOperator
and is not intended to be used directly. The assembledOperator
is expected to no longer be aLincombOperator
nor should it contain anyLincombOperators
. If an assembly of the given linear combination is not possible,None
is returned. The special case of aLincombOperator
with a single operator (i.e. a scaledOperator
) is allowed as assemble_lincomb implementsapply_inverse
for this special case.To form the linear combination of backend
Operators
(containing actual matrix data),_assemble_lincomb
will be called on the firstOperator
in the linear combination.Parameters
 operators
List of
Operators
O_i
whose linear combination is formed. coefficients
List of the corresponding linear coefficients
c_i
. solver_options
solver_options
for the assembled operator. name
Name of the assembled operator.
Returns
The assembled
Operator
if assembly is possible, otherwiseNone
.
line_search module¶

pymor.algorithms.line_search.
armijo
(f, starting_point, direction, grad=None, initial_value=None, alpha_init=1.0, tau=0.5, beta=0.0001, maxiter=10)[source]¶ Armijo line search algorithm.
This method computes a step size such that the Armijo condition (see [NW06], p. 33) is fulfilled.
Parameters
 f
Realvalued function that can be evaluated for its value.
 starting_point
A
VectorArray
of length 1 containing the starting point of the line search. direction
Descent direction along which the line search is performed.
 grad
Gradient of
f
in the pointstarting_point
. initial_value
Value of
f
in the pointstarting_point
. alpha_init
Initial step size that is gradually reduced.
 tau
The fraction by which the step size is reduced in each iteration.
 beta
Control parameter to adjust the required decrease of the function value of
f
. maxiter
Use
alpha_init
as default if the iteration count reaches this value without finding a point fulfilling the Armijo condition.
Returns
 alpha
Step size computed according to the Armijo condition.
Defaults
alpha_init, tau, beta, maxiter (see
pymor.core.defaults
)
lradi module¶

pymor.algorithms.lradi.
lyap_lrcf_solver_options
(lradi_tol=1e10, lradi_maxiter=500, lradi_shifts='projection_shifts', projection_shifts_init_maxiter=20, projection_shifts_init_seed=None)[source]¶ Return available Lyapunov solvers with default options.
Parameters
 lradi_tol
See
solve_lyap_lrcf
. lradi_maxiter
See
solve_lyap_lrcf
. lradi_shifts
See
solve_lyap_lrcf
. projection_shifts_init_maxiter
 projection_shifts_init_seed
Returns
A dict of available solvers with default solver options.
Defaults
lradi_tol, lradi_maxiter, lradi_shifts, projection_shifts_init_maxiter, projection_shifts_init_seed (see
pymor.core.defaults
)

pymor.algorithms.lradi.
projection_shifts
(A, E, V, prev_shifts)[source]¶ Find further shift parameters for lowrank ADI iteration using Galerkin projection on spaces spanned by LRADI iterates.
See [PK16], pp. 9295.
Parameters
 A
The
Operator
A from the corresponding Lyapunov equation. E
The
Operator
E from the corresponding Lyapunov equation. V
A
VectorArray
representing the currently computed iterate. prev_shifts
A
NumPy array
containing the set of all previously used shift parameters.
Returns
 shifts
A
NumPy array
containing a set of stable shift parameters.

pymor.algorithms.lradi.
projection_shifts_init
(A, E, B, shift_options)[source]¶ Find starting shift parameters for lowrank ADI iteration using Galerkin projection on spaces spanned by LRADI iterates.
See [PK16], pp. 9295.
Parameters
 A
The
Operator
A from the corresponding Lyapunov equation. E
The
Operator
E from the corresponding Lyapunov equation. B
The
VectorArray
B from the corresponding Lyapunov equation. shift_options
The shift options to use (see
lyap_lrcf_solver_options
).
Returns
 shifts
A
NumPy array
containing a set of stable shift parameters.

pymor.algorithms.lradi.
solve_lyap_lrcf
(A, E, B, trans=False, options=None)[source]¶ Compute an approximate lowrank solution of a Lyapunov equation.
See
pymor.algorithms.lyapunov.solve_lyap_lrcf
for a general description.This function uses the lowrank ADI iteration as described in Algorithm 4.3 in [PK16].
Parameters
 A
The nonparametric
Operator
A. E
The nonparametric
Operator
E orNone
. B
The operator B as a
VectorArray
fromA.source
. trans
Whether the first
Operator
in the Lyapunov equation is transposed. options
The solver options to use (see
lyap_lrcf_solver_options
).
Returns
 Z
Lowrank Cholesky factor of the Lyapunov equation solution,
VectorArray
fromA.source
.
lrradi module¶

pymor.algorithms.lrradi.
hamiltonian_shifts
(A, E, B, R, K, Z, shift_options)[source]¶ Compute further shift parameters for lowrank RADI iteration.
Compute Galerkin projection of Hamiltonian matrix on space spanned by last few columns of \(Z\) and return the eigenvalue of the projected Hamiltonian with the most impact on convergence as the next shift parameter.
See [BBKS18], pp. 318321.
Parameters
 A
The
Operator
A from the corresponding Riccati equation. E
The
Operator
E from the corresponding Riccati equation. B
The
VectorArray
B from the corresponding Riccati equation. R
A
VectorArray
representing the currently computed residual factor. K
A
VectorArray
representing the currently computed iterate. Z
A
VectorArray
representing the currently computed solution factor. shift_options
The shift options to use (see
ricc_lrcf_solver_options
).
Returns
 shifts
A
NumPy array
containing a set of stable shift parameters.

pymor.algorithms.lrradi.
hamiltonian_shifts_init
(A, E, B, C, shift_options)[source]¶ Compute initial shift parameters for lowrank RADI iteration.
Compute Galerkin projection of Hamiltonian matrix on space spanned by \(C\) and return the eigenvalue of the projected Hamiltonian with the most impact on convergence as the next shift parameter.
See [BBKS18], pp. 318321.
Parameters
 A
The
Operator
A from the corresponding Riccati equation. E
The
Operator
E from the corresponding Riccati equation. B
The
VectorArray
B from the corresponding Riccati equation. C
The
VectorArray
C from the corresponding Riccati equation. shift_options
The shift options to use (see
ricc_lrcf_solver_options
).
Returns
 shifts
A
NumPy array
containing a set of stable shift parameters.

pymor.algorithms.lrradi.
ricc_lrcf_solver_options
(lrradi_tol=1e10, lrradi_maxiter=500, lrradi_shifts='hamiltonian_shifts', hamiltonian_shifts_init_maxiter=20, hamiltonian_shifts_init_seed=None, hamiltonian_shifts_subspace_columns=6)[source]¶ Returns available Riccati equation solvers with default solver options.
Parameters
 lrradi_tol
See
solve_ricc_lrcf
. lrradi_maxiter
See
solve_ricc_lrcf
. lrradi_shifts
See
solve_ricc_lrcf
. hamiltonian_shifts_init_maxiter
 hamiltonian_shifts_init_seed
 hamiltonian_shifts_subspace_columns
See
hamiltonian_shifts
.
Returns
A dict of available solvers with default solver options.
Defaults
lrradi_tol, lrradi_maxiter, lrradi_shifts, hamiltonian_shifts_init_maxiter, hamiltonian_shifts_init_seed, hamiltonian_shifts_subspace_columns (see
pymor.core.defaults
)

pymor.algorithms.lrradi.
solve_ricc_lrcf
(A, E, B, C, R=None, trans=False, options=None)[source]¶ Compute an approximate lowrank solution of a Riccati equation.
See
pymor.algorithms.riccati.solve_ricc_lrcf
for a general description.This function is an implementation of Algorithm 2 in [BBKS18].
Parameters
 A
The
Operator
A. E
The
Operator
E orNone
. B
The operator B as a
VectorArray
fromA.source
. C
The operator C as a
VectorArray
fromA.source
. R
The operator R as a 2D
NumPy array
orNone
. trans
Whether the first
Operator
in the Riccati equation is transposed. options
The solver options to use. (see
ricc_lrcf_solver_options
)
Returns
 Z
Lowrank Cholesky factor of the Riccati equation solution,
VectorArray
fromA.source
.
lyapunov module¶

pymor.algorithms.lyapunov.
_chol
(A)[source]¶ Cholesky decomposition.
This implementation uses SVD to compute the Cholesky factor (can be used for singular matrices).
Parameters
 A
Symmetric positive semidefinite matrix as a
NumPy array
.
Returns
 L
Cholesky factor of A (in the sense that L * L^T approximates A).

pymor.algorithms.lyapunov.
mat_eqn_sparse_min_size
(value=1000)[source]¶ Returns minimal size for which a sparse solver will be used by default.
Defaults
value (see
pymor.core.defaults
)

pymor.algorithms.lyapunov.
solve_lyap_dense
(A, E, B, trans=False, options=None, default_solver_backend='pymess')[source]¶ Compute the solution of a Lyapunov equation.
Returns the solution \(X\) of a (generalized) continuoustime algebraic Lyapunov equation:
if trans is
False
and E isNone
:\[A X + X A^T + B B^T = 0,\]if trans is
False
and E is anOperator
:\[A X E^T + E X A^T + B B^T = 0,\]if trans is
True
and E isNone
:\[A^T X + X A + B^T B = 0,\]if trans is
True
and E is anOperator
:\[A^T X E + E^T X A + B^T B = 0.\]
We assume A and E are real
NumPy arrays
, E is invertible, and that no two eigenvalues of (A, E) sum to zero (i.e., there exists a unique solution X).If the solver is not specified using the options argument, a solver backend is chosen based on availability in the following order:
pymess
(seepymor.bindings.pymess.solve_lyap_dense
)slycot
(seepymor.bindings.slycot.solve_lyap_dense
)scipy
(seepymor.bindings.scipy.solve_lyap_dense
)
Parameters
 A
The operator A as a 2D
NumPy array
. E
The operator E as a 2D
NumPy array
orNone
. B
The operator B as a 2D
NumPy array
. trans
Whether the first operator in the Lyapunov equation is transposed.
 options
The solver options to use. See:
 default_solver_backend
Default solver backend to use (pymess, slycot, scipy).
Returns
 X
Lyapunov equation solution as a
NumPy array
.
Defaults
default_solver_backend (see
pymor.core.defaults
)

pymor.algorithms.lyapunov.
solve_lyap_lrcf
(A, E, B, trans=False, options=None, default_sparse_solver_backend='pymess', default_dense_solver_backend='pymess')[source]¶ Compute an approximate lowrank solution of a Lyapunov equation.
Returns a lowrank Cholesky factor \(Z\) such that \(Z Z^T\) approximates the solution \(X\) of a (generalized) continuoustime algebraic Lyapunov equation:
if trans is
False
and E isNone
:\[A X + X A^T + B B^T = 0,\]if trans is
False
and E is anOperator
:\[A X E^T + E X A^T + B B^T = 0,\]if trans is
True
and E isNone
:\[A^T X + X A + B^T B = 0,\]if trans is
True
and E is anOperator
:\[A^T X E + E^T X A + B^T B = 0.\]
We assume A and E are real
Operators
, E is invertible, and all the eigenvalues of (A, E) all lie in the open left halfplane. Operator B needs to be given as aVectorArray
fromA.source
, and for largescale problems, we assumelen(B)
is small.If the solver is not specified using the options argument, a solver backend is chosen based on availability in the following order:
for sparse problems (minimum size specified by
mat_eqn_sparse_min_size
)pymess
(seepymor.bindings.pymess.solve_lyap_lrcf
),lradi
(seepymor.algorithms.lradi.solve_lyap_lrcf
),
for dense problems (smaller than
mat_eqn_sparse_min_size
)pymess
(seepymor.bindings.pymess.solve_lyap_lrcf
),slycot
(seepymor.bindings.slycot.solve_lyap_lrcf
),scipy
(seepymor.bindings.scipy.solve_lyap_lrcf
).
Parameters
 A
The nonparametric
Operator
A. E
The nonparametric
Operator
E orNone
. B
The operator B as a
VectorArray
fromA.source
. trans
Whether the first
Operator
in the Lyapunov equation is transposed. options
The solver options to use. See:
 default_sparse_solver_backend
Default sparse solver backend to use (pymess, lradi).
 default_dense_solver_backend
Default dense solver backend to use (pymess, slycot, scipy).
Returns
 Z
Lowrank Cholesky factor of the Lyapunov equation solution,
VectorArray
fromA.source
.
Defaults
default_sparse_solver_backend, default_dense_solver_backend (see
pymor.core.defaults
)
newton module¶

pymor.algorithms.newton.
newton
(operator, rhs, initial_guess=None, mu=None, range_product=None, source_product=None, least_squares=False, miniter=0, maxiter=100, atol=0.0, rtol=1e07, relax='armijo', line_search_params=None, stagnation_window=3, stagnation_threshold=inf, error_measure='update', return_stages=False, return_residuals=False)[source]¶ Newton algorithm.
This method solves the nonlinear equation
A(U, mu) = V
for
U
using the Newton method.Parameters
 operator
The
Operator
A
.A
must implement thejacobian
interface method. rhs
VectorArray
of length 1 containing the vectorV
. initial_guess
If not
None
, aVectorArray
of length 1 containing an initial guess for the solutionU
. mu
The
parameter values
for which to solve the equation. range_product
The inner product
Operator
onoperator.range
with which the norm of the resiudal is computed. IfNone
, the Euclidean inner product is used. source_product
The inner product
Operator
onoperator.source
with which the norm of the solution and update vectors is computed. IfNone
, the Euclidean inner product is used. least_squares
If
True
, use a least squares linear solver (e.g. for residual minimization). miniter
Minimum amount of iterations to perform.
 maxiter
Fail if the iteration count reaches this value without converging.
 atol
Finish when the error measure is below this threshold.
 rtol
Finish when the error measure has been reduced by this factor relative to the norm of the initial residual resp. the norm of the current solution.
 relax
If real valued, relaxation factor for Newton updates; otherwise
'armijo'
to indicate that thearmijo
line search algorithm shall be used. line_search_params
Dictionary of additional parameters passed to the line search method.
 stagnation_window
Finish when the error measure has not been reduced by a factor of
stagnation_threshold
during the laststagnation_window
iterations. stagnation_threshold
See
stagnation_window
. error_measure
If
'residual'
, convergence depends on the norm of the residual. If'update'
, convergence depends on the norm of the update vector. return_stages
If
True
, return aVectorArray
of the intermediate approximations ofU
after each iteration. return_residuals
If
True
, return aVectorArray
of all residual vectors which have been computed during the Newton iterations.
Returns
 U
VectorArray
of length 1 containing the computed solution data
Dict containing the following fields:
 solution_norms
NumPy array
of the solution norms after each iteration. update_norms
NumPy array
of the norms of the update vectors for each iteration. residual_norms
NumPy array
of the residual norms after each iteration. stages
See
return_stages
. residuals
See
return_residuals
.
Raises
 NewtonError
Raised if the Newton algorithm failed to converge.
Defaults
miniter, maxiter, rtol, atol, relax, stagnation_window, stagnation_threshold (see
pymor.core.defaults
)
pod module¶

pymor.algorithms.pod.
pod
(A, product=None, modes=None, rtol=1e07, atol=0.0, l2_err=0.0, method='method_of_snapshots', orth_tol=1e10)[source]¶ Proper orthogonal decomposition of
A
.Viewing the
VectorArray
A
as aA.dim
xlen(A)
matrix, the return values of this method are theVectorArray
of left singular vectors and aNumPy array
of singular values of the singular value decomposition ofA
, where the inner product on R^(dim(A)
) is given byproduct
and the inner product on R^(len(A)
) is the Euclidean inner product.Parameters
 A
The
VectorArray
for which the POD is to be computed. product
Inner product
Operator
w.r.t. which the POD is computed. modes
If not
None
, at most the firstmodes
POD modes (singular vectors) are returned. rtol
Singular values smaller than this value multiplied by the largest singular value are ignored.
 atol
Singular values smaller than this value are ignored.
 l2_err
Do not return more modes than needed to bound the l2approximation error by this value. I.e. the number of returned modes is at most
argmin_N { sum_{n=N+1}^{infty} s_n^2 <= l2_err^2 }
where
s_n
denotes the nth singular value. method
Which SVD method from
svd_va
to use ('method_of_snapshots'
or'qr_svd'
). orth_tol
POD modes are reorthogonalized if the orthogonality error is above this value.
Returns
 POD
VectorArray
of POD modes. SVALS
Onedimensional
NumPy array
of singular values.
Defaults
rtol, atol, l2_err, method, orth_tol (see
pymor.core.defaults
)
preassemble module¶

class
pymor.algorithms.preassemble.
PreAssembleRules
[source]¶ Bases:
pymor.algorithms.rules.RuleTable
Methods
append_rule
,apply
,apply_children
,breakpoint_for_name
,breakpoint_for_obj
,get_children
,insert_rule
,replace_children
Attributes
action_AdjointOperator
,action_assemble
,action_identity
,action_recurse
,action_recurse_and_assemble
,rules

pymor.algorithms.preassemble.
preassemble
(obj)[source]¶ Preassemble nonparametric operators.
If
obj
is a nonparametricOperator
, returnobj.assemble()
otherwise returnobj
. Recursively replaces children ofobj
.
projection module¶

class
pymor.algorithms.projection.
ProjectRules
(range_basis, source_basis)[source]¶ Bases:
pymor.algorithms.rules.RuleTable
RuleTable
for theproject
algorithm.Methods
append_rule
,apply
,apply_children
,breakpoint_for_name
,breakpoint_for_obj
,get_children
,insert_rule
,replace_children
Attributes
action_AdjointOperator
,action_AffineOperator
,action_apply_basis
,action_BlockOperatorBase
,action_ConcatenationOperator
,action_ConstantOperator
,action_EmpiricalInterpolatedOperator
,action_LincombOperator
,action_no_bases
,action_SelectionOperator
,action_ZeroOperator
,rules

class
pymor.algorithms.projection.
ProjectToSubbasisRules
(dim_range, dim_source)[source]¶ Bases:
pymor.algorithms.rules.RuleTable
RuleTable
for theproject_to_subbasis
algorithm.Methods
append_rule
,apply
,apply_children
,breakpoint_for_name
,breakpoint_for_obj
,get_children
,insert_rule
,replace_children
Attributes
action_BlockColumnOperator
,action_BlockRowOperator
,action_ConstantOperator
,action_IdentityOperator
,action_NumpyMatrixOperator
,action_ProjectedEmpiciralInterpolatedOperator
,action_ProjectedOperator
,action_recurse
,action_VectorArrayOperator
,action_ZeroOperator
,rules

pymor.algorithms.projection.
project
(op, range_basis, source_basis, product=None)[source]¶ PetrovGalerkin projection of a given
Operator
.Given an inner product
( ⋅, ⋅)
, source vectorsb_1, ..., b_N
and range vectorsc_1, ..., c_M
, the projectionop_proj
ofop
is defined by[ op_proj(e_j) ]_i = ( c_i, op(b_j) )
for all i,j, where
e_j
denotes the jth canonical basis vector of R^N.In particular, if the
c_i
are orthonormal w.r.t. the given product, thenop_proj
is the coordinate representation w.r.t. theb_i/c_i
bases of the restriction ofop
tospan(b_i)
concatenated with the orthogonal projection ontospan(c_i)
.From another point of view, if
op
is viewed as a bilinear form (seeapply2
) and( ⋅, ⋅ )
is the Euclidean inner product, thenop_proj
represents the matrix of the bilinear form restricted tospan(b_i) / span(c_i)
(w.r.t. theb_i/c_i
bases).How the projection is realized will depend on the given
Operator
. While a projectedNumpyMatrixOperator
will again be aNumpyMatrixOperator
, only a genericProjectedOperator
can be returned in general. The exact algorithm is specified inProjectRules
.Parameters
 range_basis
The vectors
c_1, ..., c_M
as aVectorArray
. IfNone
, no projection in the range space is performed. source_basis
The vectors
b_1, ..., b_N
as aVectorArray
orNone
. IfNone
, no restriction of the source space is performed. product
An
Operator
representing the inner product. IfNone
, the Euclidean inner product is chosen.
Returns
The projected
Operator
op_proj
.

pymor.algorithms.projection.
project_to_subbasis
(op, dim_range=None, dim_source=None)[source]¶ Project already projected
Operator
to a subbasis.The purpose of this method is to further project an operator that has been obtained through
project
to subbases of the original projection bases, i.e.project_to_subbasis(project(op, r_basis, s_basis, prod), dim_range, dim_source)
should be the same as
project(op, r_basis[:dim_range], s_basis[:dim_source], prod)
For a
NumpyMatrixOperator
this amounts to extracting the upperleft (dim_range, dim_source) corner of its matrix.The subbasis projection algorithm is specified in
ProjectToSubbasisRules
.Parameters
 dim_range
Dimension of the range subbasis.
 dim_source
Dimension of the source subbasis.
Returns
The projected
Operator
.
randrangefinder module¶

pymor.algorithms.randrangefinder.
adaptive_rrf
(A, source_product=None, range_product=None, tol=0.0001, failure_tolerance=1e15, num_testvecs=20, lambda_min=None, iscomplex=False)[source]¶ Adaptive randomized range approximation of
A
.This is an implementation of Algorithm 1 in [BS18].
Given the
Operator
A
, the return value of this method is theVectorArray
B
with the property\[\Vert A  P_{span(B)} A \Vert \leq tol\]with a failure probability smaller than
failure_tolerance
, where the norm denotes the operator norm. The inner product of the range ofA
is given byrange_product
and the inner product of the source ofA
is given bysource_product
.Parameters
 A
The
Operator
A. source_product
Inner product
Operator
of the source of A. range_product
Inner product
Operator
of the range of A. tol
Error tolerance for the algorithm.
 failure_tolerance
Maximum failure probability.
 num_testvecs
Number of test vectors.
 lambda_min
The smallest eigenvalue of source_product. If
None
, the smallest eigenvalue is computed using scipy. iscomplex
If
True
, the random vectors are chosen complex.
Returns
 B
VectorArray
which contains the basis, whose span approximates the range of A.
Defaults
tol, failure_tolerance, num_testvecs (see
pymor.core.defaults
)

pymor.algorithms.randrangefinder.
rrf
(A, source_product=None, range_product=None, q=2, l=8, iscomplex=False)[source]¶ Randomized range approximation of
A
.This is an implementation of Algorithm 4.4 in [HMT11].
Given the
Operator
A
, the return value of this method is theVectorArray
Q
whose vectors form an approximate orthonomal basis for the range ofA
.Parameters
Returns
 Q
VectorArray
which contains the basis, whose span approximates the range of A.
Defaults
q, l (see
pymor.core.defaults
)
riccati module¶

pymor.algorithms.riccati.
solve_pos_ricc_lrcf
(A, E, B, C, R=None, trans=False, options=None, default_dense_solver_backend='pymess')[source]¶ Compute an approximate lowrank solution of a positive Riccati equation.
Returns a lowrank Cholesky factor \(Z\) such that \(Z Z^T\) approximates the solution \(X\) of a (generalized) positive continuoustime algebraic Riccati equation:
if trans is
False
\[A X E^T + E X A^T + E X C^T R^{1} C X E^T + B B^T = 0.\]if trans is
True
\[A^T X E + E^T X A + E^T X B R^{1} B^T X E + C^T C = 0.\]
If E is None, it is taken to be identity, and similarly for R.
If the solver is not specified using the options argument, a solver backend is chosen based on availability in the following order:
pymess
(seepymor.bindings.pymess.solve_pos_ricc_lrcf
),slycot
(seepymor.bindings.slycot.solve_pos_ricc_lrcf
),scipy
(seepymor.bindings.scipy.solve_pos_ricc_lrcf
).
Parameters
 A
The nonparametric
Operator
A. E
The nonparametric
Operator
E orNone
. B
The operator B as a
VectorArray
fromA.source
. C
The operator C as a
VectorArray
fromA.source
. R
The operator R as a 2D
NumPy array
orNone
. trans
Whether the first
Operator
in the positive Riccati equation is transposed. options
The solver options to use. See:
 default_dense_solver_backend
Default dense solver backend to use (pymess, slycot, scipy).
Returns
 Z
Lowrank Cholesky factor of the positive Riccati equation solution,
VectorArray
fromA.source
.
Defaults
default_dense_solver_backend (see
pymor.core.defaults
)

pymor.algorithms.riccati.
solve_ricc_lrcf
(A, E, B, C, R=None, trans=False, options=None, default_sparse_solver_backend='pymess', default_dense_solver_backend='pymess')[source]¶ Compute an approximate lowrank solution of a Riccati equation.
Returns a lowrank Cholesky factor \(Z\) such that \(Z Z^T\) approximates the solution \(X\) of a (generalized) continuoustime algebraic Riccati equation:
if trans is
False
\[A X E^T + E X A^T  E X C^T R^{1} C X E^T + B B^T = 0.\]if trans is
True
\[A^T X E + E^T X A  E^T X B R^{1} B^T X E + C^T C = 0.\]
If E is None, it is taken to be identity, and similarly for R.
We assume:
A and E are real
Operators
,B and C are real
VectorArrays
fromA.source
,R is a real
NumPy array
,(E, A, B, C) is stabilizable and detectable, and
R is symmetric positive definite.
For largescale problems, we additionally assume that
len(B)
andlen(C)
are small.If the solver is not specified using the options argument, a solver backend is chosen based on availability in the following order:
for sparse problems (minimum size specified by
mat_eqn_sparse_min_size
)pymess
(seepymor.bindings.pymess.solve_ricc_lrcf
),lrradi
(seepymor.algorithms.lrradi.solve_ricc_lrcf
),
for dense problems (smaller than
mat_eqn_sparse_min_size
)pymess
(seepymor.bindings.pymess.solve_ricc_lrcf
),slycot
(seepymor.bindings.slycot.solve_ricc_lrcf
),scipy
(seepymor.bindings.scipy.solve_ricc_lrcf
).
Parameters
 A
The nonparametric
Operator
A. E
The nonparametric
Operator
E orNone
. B
The operator B as a
VectorArray
fromA.source
. C
The operator C as a
VectorArray
fromA.source
. R
The operator R as a 2D
NumPy array
orNone
. trans
Whether the first
Operator
in the Riccati equation is transposed. options
The solver options to use. See:
 default_sparse_solver_backend
Default sparse solver backend to use (pymess, lrradi).
 default_dense_solver_backend
Default dense solver backend to use (pymess, slycot, scipy).
Returns
 Z
Lowrank Cholesky factor of the Riccati equation solution,
VectorArray
fromA.source
.
Defaults
default_sparse_solver_backend, default_dense_solver_backend (see
pymor.core.defaults
)
rules module¶

class
pymor.algorithms.rules.
RuleTable
(use_caching=False)[source]¶ Bases:
pymor.core.base.BasicObject
Define algorithm by a table of match conditions and corresponding actions.
RuleTable
manages a table ofrules
, stored in therules
attributes, which can beapplied
to given objects.A new table is created by subclassing
RuleTable
and defining new methods which are decorated withmatch_class
,match_generic
or anotherrule
subclass. The order of the method definitions determines the order in which the definedrules
are applied.Parameters
 use_caching
If
True
, cache results ofapply
.
Methods
append_rule
,apply
,apply_children
,breakpoint_for_name
,breakpoint_for_obj
,get_children
,insert_rule
,replace_children
Attributes

apply
(obj)[source]¶ Sequentially apply rules to given object.
This method iterates over all rules of the given
RuleTable
. For eachrule
, it is checked if itmatches
the given object. IfFalse
, the nextrule
in the table is considered. IfTrue
the correspondingaction
is executed withobj
as parameter. If execution ofaction
raisesRuleNotMatchingError
, the rule is considered as not matching, and execution continues with evaluation of the next rule. Otherwise, execution is stopped and the return value ofrule.action
is returned to the caller.If no
rule
matches, aNoMatchingRuleError
is raised.Parameters
 obj
The object to apply the
RuleTable
to.
Returns
Return value of the action of the first matching
rule
in the table.Raises
 NoMatchingRuleError
No
rule
could be applied to the given object.

apply_children
(obj, children=None)[source]¶ Apply rules to all children of the given object.
This method calls
apply
to each child of the given object. The children of the object are either provided by thechildren
parameter or automatically inferred by theget_children
method.Parameters
 obj
The object to apply the
RuleTable
to. children
None
or a list of attribute names defining the children to consider.
Returns
Result of : meth:
apply
for all given children.

classmethod
breakpoint_for_name
(name)[source]¶ Add a conditional breakpoint for objects of given name.
Break execution in
apply
, when being applied to an object with a certain name.Parameters
 name
name
of the object for which to add the conditional breakpoint.

classmethod
breakpoint_for_obj
(obj)[source]¶ Add a conditional breakpoint for given object.
Break execution in
apply
, when being applied to a certain object.Parameters
 obj
Object for which to add the conditional breakpoint.

classmethod
get_children
(obj)[source]¶ Determine children of given object.
This method returns a list of the names of all attributes
a
, for which one of the folling is true:

replace_children
(obj, children=None)[source]¶ Replace children of object according to rule table.
Same as
apply_children
, but additionally callsobj.with_
to replace the children ofobj
with the result of the correspondingapply
call.

class
pymor.algorithms.rules.
RuleTableMeta
(name, parents, dct)[source]¶ Bases:
pymor.core.base.UberMeta
Meta class for
RuleTable
.Methods
register
,__instancecheck__
,__subclasscheck__
type
mro
,__dir__
,__prepare__
,__sizeof__
,__subclasses__

static
__new__
(cls, name, parents, dct)[source]¶ I copy docstrings from base class methods to deriving classes.
Copying of docstrings is disabled when the
PYMOR_WITH_SPHINX
environment variable is set to1
.

__str__
()¶ Return repr(self).

static

class
pymor.algorithms.rules.
match_always
(action)[source]¶ Bases:
pymor.algorithms.rules.rule
rule
that always matches.Attributes
condition_type
action
,action_description
,condition_description
,source

class
pymor.algorithms.rules.
match_class
(*classes)[source]¶ Bases:
pymor.algorithms.rules.match_class_base
rule
that matches when obj is instance of one of the given classes.Attributes
condition_type
action
,action_description
,condition_description
,source

class
pymor.algorithms.rules.
match_class_all
(*classes)[source]¶ Bases:
pymor.algorithms.rules.match_class_base
rule
that matches when each item of obj is instance of one of the given classes.Attributes
condition_type
action
,action_description
,condition_description
,source

class
pymor.algorithms.rules.
match_class_any
(*classes)[source]¶ Bases:
pymor.algorithms.rules.match_class_base
rule
that matches when any item of obj is instance of one of the given classes.Attributes
condition_type
action
,action_description
,condition_description
,source

class
pymor.algorithms.rules.
match_class_base
(*classes)[source]¶ Bases:
pymor.algorithms.rules.rule

class
pymor.algorithms.rules.
match_generic
(condition, condition_description=None)[source]¶ Bases:
pymor.algorithms.rules.rule
rule
with matching condition given by an arbitrary function.Parameters
 condition
Function of one argument which checks if given object matches condition.
 condition_description
Optional string describing the condition implemented by
condition
.
Attributes
condition_type
action
,action_description
,condition_description
,source

class
pymor.algorithms.rules.
rule
[source]¶ Bases:
object
Decorator to make a method a rule in a given
RuleTable
.The decorated function will become the
action
to perform in case the rulematches
. Matching conditions are specified by subclassing and overriding thematches
method.If an action is decorated by multiple rules, all these rules must match for the action to apply.

action
¶ Method to call in case the rule matches.

samdp module¶

pymor.algorithms.samdp.
_select_max_eig
(H, G, X, V, B, C, which)[source]¶ Compute poles sorted from largest to smallest residual.
Parameters
 H
The
NumPy array
H from the SAMDP algorithm. G
The
NumPy array
G from the SAMDP algorithm. X
A
VectorArray
describing the orthogonal search space used in the SAMDP algorithm. V
A
VectorArray
describing the orthogonal search space used in the SAMDP algorithm. B
The
VectorArray
B from the corresponding LTI system modified by deflation. C
The
VectorArray
C from the corresponding LTI system modified by deflation. which
A string that indicates which poles to select. See
samdp
.
Returns
 poles
A
NumPy array
containing poles sorted according to the chosen strategy. rightevs
A
NumPy array
containing the right eigenvectors of the computed poles. leftevs
A
NumPy array
containing the left eigenvectors of the computed poles. residue
A 1D
NumPy array
containing the norms of the residues.

pymor.algorithms.samdp.
_twosided_rqi
(A, E, x, y, theta, init_res, imagtol, rqitol, maxiter)[source]¶ Refine an initial guess for an eigentriplet of the matrix pair (A, E).
Parameters
 A
The
Operator
A from the LTI system. E
The
Operator
E from the LTI system. x
Initial guess for right eigenvector of matrix pair (A, E).
 y
Initial guess for left eigenvector of matrix pair (A, E).
 theta
Initial guess for eigenvalue of matrix pair (A, E).
 init_res
Residual of initial guess.
 imagtol
Relative tolerance for imaginary parts of pairs of complex conjugate eigenvalues.
 rqitol
Convergence tolerance for the relative change of the pole.
 maxiter
Maximum number of iteration.
Returns
 x
Refined right eigenvector of matrix pair (A, E).
 y
Refined left eigenvector of matrix pair (A, E).
 theta
Refined eigenvalue of matrix pair (A, E).
 residual
Residual of the computed triplet.

pymor.algorithms.samdp.
samdp
(A, E, B, C, nwanted, init_shifts=None, which='LR', tol=1e10, imagtol=1e06, conjtol=1e08, dorqitol=0.0001, rqitol=1e10, maxrestart=100, krestart=20, rqi_maxiter=10, seed=0)[source]¶ Compute the dominant pole triplets and residues of the transfer function of an LTI system.
This function uses the subspace accelerated dominant pole (SAMDP) algorithm as described in [RM06] in Algorithm 2 in order to compute dominant pole triplets and residues of the transfer function
\[H(s) = C (s E  A)^{1} B\]of an LTI system. It is possible to take advantage of prior knowledge about the poles by specifying shift parameters, which are injected after a new pole has been found.
Parameters
 A
The
Operator
A. E
The
Operator
E orNone
. B
The operator B as a
VectorArray
fromA.source
. C
The operator C as a
VectorArray
fromA.source
. nwanted
The number of dominant poles that should be computed.
 init_shifts
A
NumPy array
containing shifts which are injected after a new pole has been found. which
A string specifying the strategy by which the dominant poles and residues are selected. Possible values are:
'LR'
: select poles with largest norm(residual) / abs(Re(pole))'LS'
: select poles with largest norm(residual) / abs(pole)'LM'
: select poles with largest norm(residual)
 tol
Tolerance for the residual of the poles.
 imagtol
Relative tolerance for imaginary parts of pairs of complex conjugate eigenvalues.
 conjtol
Tolerance for the residual of the complex conjugate of a pole.
 dorqitol
If the residual is smaller than dorqitol the twosided Rayleigh quotient iteration is executed.
 rqitol
Tolerance for the relative change of a pole in the twosided Rayleigh quotient iteration.
 maxrestart
The maximum number of restarts.
 krestart
Maximum dimension of search space before performing a restart.
 rqi_maxiter
Maximum number of iterations for the twosided Rayleigh quotient iteration.
 seed
Random seed which is used for computing the initial shift and random restarts.
Returns
 poles
A 1D
NumPy array
containing the computed dominant poles. residues
A 3D
NumPy array
of shape(len(poles), len(C), len(B))
containing the computed residues. rightev
A
VectorArray
containing the right eigenvectors of the computed poles. leftev
A
VectorArray
containing the left eigenvectors of the computed poles.
Defaults
which, tol, imagtol, conjtol, dorqitol, rqitol, maxrestart, krestart, init_shifts, rqi_maxiter, seed (see
pymor.core.defaults
)
simplify module¶

class
pymor.algorithms.simplify.
ExpandRules
[source]¶ Bases:
pymor.algorithms.rules.RuleTable
Methods
append_rule
,apply
,apply_children
,breakpoint_for_name
,breakpoint_for_obj
,get_children
,insert_rule
,replace_children
Attributes
action_ConcatenationOperator
,action_LincombOperator
,action_recurse
,rules

pymor.algorithms.simplify.
expand
(obj)[source]¶ Expand concatenations of LincombOperators.
To any given
Operator
orModel
, the following transformations are applied recursively:Concatenations
ofLincombOperators
are expanded. E.g.(O1 + O2) @ (O3 + O4)
becomes:
O1 @ O3 + O1 @ O4 + O2 @ O3 + O2 @ O4
LincombOperators
insideLincombOperators
are merged into a singleLincombOperator
ConcatenationOperators
insideConcatenationOperators
are merged into a singleConcatenationOperator
.
Parameters
Returns
The transformed object.
svd_va module¶
Module for SVD method of operators represented by VectorArrays
.

pymor.algorithms.svd_va.
method_of_snapshots
(A, product=None, modes=None, rtol=1e07, atol=0.0, l2_err=0.0)[source]¶ SVD of a
VectorArray
using the method of snapshots.Viewing the
VectorArray
A
as aA.dim
xlen(A)
matrix, the return value of this method is the singular value decomposition ofA
, where the inner product on R^(dim(A)
) is given byproduct
and the inner product on R^(len(A)
) is the Euclidean inner product.Warning
The left singular vectors may not be numerically orthonormal for illconditioned
A
.Parameters
 A
The
VectorArray
for which the SVD is to be computed. product
Inner product
Operator
w.r.t. which the SVD is computed. modes
If not
None
, at most the firstmodes
singular values and vectors are returned. rtol
Singular values smaller than this value multiplied by the largest singular value are ignored.
 atol
Singular values smaller than this value are ignored.
 l2_err
Do not return more modes than needed to bound the l2approximation error by this value. I.e. the number of returned modes is at most
argmin_N { sum_{n=N+1}^{infty} s_n^2 <= l2_err^2 }
where
s_n
denotes the nth singular value.
Returns
 U
VectorArray
of left singular vectors. s
Onedimensional
NumPy array
of singular values. Vh
NumPy array
of right singular vectors.
Defaults
rtol, atol, l2_err (see
pymor.core.defaults
)

pymor.algorithms.svd_va.
qr_svd
(A, product=None, modes=None, rtol=4e08, atol=0.0, l2_err=0.0)[source]¶ SVD of a
VectorArray
using GramSchmidt orthogonalization.Viewing the
VectorArray
A
as aA.dim
xlen(A)
matrix, the return value of this method is the singular value decomposition ofA
, where the inner product on R^(dim(A)
) is given byproduct
and the inner product on R^(len(A)
) is the Euclidean inner product.Parameters
 A
The
VectorArray
for which the SVD is to be computed. product
Inner product
Operator
w.r.t. which the left singular vectors are computed. modes
If not
None
, at most the firstmodes
singular values and vectors are returned. rtol
Singular values smaller than this value multiplied by the largest singular value are ignored.
 atol
Singular values smaller than this value are ignored.
 l2_err
Do not return more modes than needed to bound the l2approximation error by this value. I.e. the number of returned modes is at most
argmin_N { sum_{n=N+1}^{infty} s_n^2 <= l2_err^2 }
where
s_n
denotes the nth singular value.
Returns
 U
VectorArray
of left singular vectors. s
Onedimensional
NumPy array
of singular values. Vh
NumPy array
of right singular vectors.
Defaults
rtol, atol, l2_err (see
pymor.core.defaults
)
sylvester module¶

pymor.algorithms.sylvester.
solve_sylv_schur
(A, Ar, E=None, Er=None, B=None, Br=None, C=None, Cr=None)[source]¶ Solve Sylvester equation by Schur decomposition.
Solves Sylvester equation
\[A V E_r^T + E V A_r^T + B B_r^T = 0\]or
\[A^T W E_r + E^T W A_r + C^T C_r = 0\]or both using (generalized) Schur decomposition (Algorithms 3 and 4 in [BKS11]), if the necessary parameters are given.
Parameters
 A
Real
Operator
. Ar
Real
Operator
. It is converted into aNumPy array
usingto_matrix
. E
Real
Operator
orNone
(then assumed to be the identity). Er
Real
Operator
orNone
(then assumed to be the identity). It is converted into aNumPy array
usingto_matrix
. B
Real
Operator
orNone
. Br
Real
Operator
orNone
. It is assumed thatBr.range.from_numpy
is implemented. C
Real
Operator
orNone
. Cr
Real
Operator
orNone
. It is assumed thatCr.source.from_numpy
is implemented.
Returns
 V
Returned if
B
andBr
are given,VectorArray
fromA.source
. W
Returned if
C
andCr
are given,VectorArray
fromA.source
.
Raises
 ValueError
If
V
andW
cannot be returned.
timestepping module¶
This module provides generic timestepping algorithms for the solution of instationary problems.
The algorithms are generic in the sense that each algorithms operates exclusively
on Operators
and VectorArrays
. In particular, the algorithms
can also be used to turn an arbitrary stationary Model
provided
by an external library into an instationary Model
.
Currently, implementations of explicit_euler
and implicit_euler
timestepping are provided. The TimeStepper
defines a
common interface that has to be fulfilled by the timesteppers used
by InstationaryModel
. The classes ExplicitEulerTimeStepper
and ImplicitEulerTimeStepper
encapsulate explicit_euler
and
implicit_euler
to provide this interface.

class
pymor.algorithms.timestepping.
ExplicitEulerTimeStepper
(nt)[source]¶ Bases:
pymor.algorithms.timestepping.TimeStepper
Explicit Euler timestepper.
Solves equations of the form
M * d_t u + A(u, mu, t) = F(mu, t).
Parameters
 nt
The number of timesteps the timestepper will perform.
Methods
Attributes

solve
(initial_time, end_time, initial_data, operator, rhs=None, mass=None, mu=None, num_values=None)[source]¶ Apply timestepper to the equation
M * d_t u + A(u, mu, t) = F(mu, t).
Parameters
 initial_time
The time at which to begin timestepping.
 end_time
The time until which to perform timestepping.
 initial_data
The solution vector at
initial_time
. operator
The
Operator
A. rhs
The righthand side F (either
VectorArray
of length 1 orOperator
withsource.dim == 1
). IfNone
, zero righthand side is assumed. mass
The
Operator
M. IfNone
, the identity operator is assumed. mu
Parameter values
for whichoperator
andrhs
are evaluated. The current time is added tomu
with keyt
. num_values
The number of returned vectors of the solution trajectory. If
None
, each intermediate vector that is calculated is returned.
Returns
VectorArray
containing the solution trajectory.

class
pymor.algorithms.timestepping.
ImplicitEulerTimeStepper
(nt, solver_options='operator')[source]¶ Bases:
pymor.algorithms.timestepping.TimeStepper
Implict Euler timestepper.
Solves equations of the form
M * d_t u + A(u, mu, t) = F(mu, t).
Parameters
 nt
The number of timesteps the timestepper will perform.
 solver_options
The
solver_options
used to invertM + dt*A
. The special values'mass'
and'operator'
are recognized, in which case the solver_options of M (resp. A) are used.
Methods
Attributes

solve
(initial_time, end_time, initial_data, operator, rhs=None, mass=None, mu=None, num_values=None)[source]¶ Apply timestepper to the equation
M * d_t u + A(u, mu, t) = F(mu, t).
Parameters
 initial_time
The time at which to begin timestepping.
 end_time
The time until which to perform timestepping.
 initial_data
The solution vector at
initial_time
. operator
The
Operator
A. rhs
The righthand side F (either
VectorArray
of length 1 orOperator
withsource.dim == 1
). IfNone
, zero righthand side is assumed. mass
The
Operator
M. IfNone
, the identity operator is assumed. mu
Parameter values
for whichoperator
andrhs
are evaluated. The current time is added tomu
with keyt
. num_values
The number of returned vectors of the solution trajectory. If
None
, each intermediate vector that is calculated is returned.
Returns
VectorArray
containing the solution trajectory.

class
pymor.algorithms.timestepping.
TimeStepper
[source]¶ Bases:
pymor.core.base.ImmutableObject
Interface for timestepping algorithms.
Algorithms implementing this interface solve timedependent problems of the form
M * d_t u + A(u, mu, t) = F(mu, t).
Timesteppers used by
InstationaryModel
have to fulfill this interface.Methods
Attributes

abstract
solve
(initial_time, end_time, initial_data, operator, rhs=None, mass=None, mu=None, num_values=None)[source]¶ Apply timestepper to the equation
M * d_t u + A(u, mu, t) = F(mu, t).
Parameters
 initial_time
The time at which to begin timestepping.
 end_time
The time until which to perform timestepping.
 initial_data
The solution vector at
initial_time
. operator
The
Operator
A. rhs
The righthand side F (either
VectorArray
of length 1 orOperator
withsource.dim == 1
). IfNone
, zero righthand side is assumed. mass
The
Operator
M. IfNone
, the identity operator is assumed. mu
Parameter values
for whichoperator
andrhs
are evaluated. The current time is added tomu
with keyt
. num_values
The number of returned vectors of the solution trajectory. If
None
, each intermediate vector that is calculated is returned.
Returns
VectorArray
containing the solution trajectory.

abstract

pymor.algorithms.timestepping.
explicit_euler
(A, F, U0, t0, t1, nt, mu=None, num_values=None)[source]¶

pymor.algorithms.timestepping.
implicit_euler
(A, F, M, U0, t0, t1, nt, mu=None, num_values=None, solver_options='operator')[source]¶
to_matrix module¶

class
pymor.algorithms.to_matrix.
ToMatrixRules
(format, mu)[source]¶ Bases:
pymor.algorithms.rules.RuleTable
Methods
append_rule
,apply
,apply_children
,breakpoint_for_name
,breakpoint_for_obj
,get_children
,insert_rule
,replace_children
Attributes
action_AdjointOperator
,action_BlockOperator
,action_ComponentProjectionOperator
,action_ConcatenationOperator
,action_IdentityOperator
,action_LincombOperator
,action_LowRankOperator
,action_LowRankUpdatedOperator
,action_NumpyMatrixOperator
,action_VectorArrayOperator
,action_ZeroOperator
,rules

pymor.algorithms.to_matrix.
to_matrix
(op, format=None, mu=None)[source]¶ Convert a linear
Operator
to a matrix.Parameters
 op
The
Operator
to convert. format
Format of the resulting matrix:
NumPy array
if ‘dense’, otherwise the appropriateSciPy spmatrix
. IfNone
, a choice between dense and sparse format is automatically made. mu
The
parameter values
for which to convertop
.
Returns
 res
The matrix equivalent to
op
.