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

`AdaptiveSampleSet`	`map_vertex_to_mu`, `refine`, `visualize`
`BasicObject`	`disable_logging`, `enable_logging`

Attributes

`AdaptiveSampleSet`	`Element`
`BasicObject`	`logger`, `logging_disabled`, `name`, `uid`

pymor.algorithms.adaptivegreedy.adaptive_greedy(*args, **kwargs)[source]¶

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 than rho, 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 non-positive 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_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 in adaptive_weak_greedy.

Parameters

fom: See rb_greedy.
reductor: See rb_greedy.
parameter_space: The ParameterSpace for which to compute the reduced model.
use_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, norm=None, rtol=1e-14, atol=1e-14)[source]¶

Compare U and V for almost equality.

The vectors of U and V are compared in pairs for almost equality. Two vectors u and v are considered almost equal iff

||u - v|| <= atol + ||v|| * rtol.

The norm to be used can be specified via the norm or product parameter.

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. product and norm are mutually exclusive.
norm: If specified, must be a callable which is used to compute the norm or, alternatively, one of the strings ‘l1’, ‘l2’, ‘sup’, in which case the respective VectorArray norm methods are used. product and norm are mutually exclusive. If neither is specified, norm='l2' is assumed.
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 dim

Parameters

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 in basis 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 in U. The returned objects can be used to instantiate an EmpiricalInterpolatedOperator (with triangular=False).

The collateral basis is determined by the first pod modes of U.

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 of U to obtain the collateral basis. If False, 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 EI-Greedy algorithm.

Given a VectorArray U, this method generates a collateral basis and interpolation DOFs for empirical interpolation of the vectors contained in U. The returned objects can be used to instantiate an EmpiricalInterpolatedOperator (with triangular=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, the WorkerPool 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 EI-Greedy/DEIM algorithm.

This is a convenience method to facilitate the use of ei_greedy or deim. Given a Model, names of Operators, and a sample of Parameters, first the operators are evaluated on the solution snapshots of the model for the provided parameters. These evaluations are then used as input for ei_greedy/deim. Finally the resulting interpolation data is used to create EmpiricalInterpolatedOperators 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: The Model whose Operators will be interpolated.
operator_names: List of keys in the operators dict of the model. The corresponding Operators will be interpolated.
parameter_sample: A list of Parameters for which solution snapshots are calculated.
error_norm: See ei_greedy. Has no effect if alg == 'deim'.
product: Inner product for POD computation in deim. Has no effect if alg == '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 if alg == 'ei_greedy'.
alg: Either ei_greedy or deim.
pool: If not None, the WorkerPool to use for parallelization.

Returns

eim

Model with Operators given by operator_names replaced by EmpiricalInterpolatedOperators.

data

Dict containing the following fields:

dofs

NumPy array of the DOFs at which the Operators have to be evaluated.

basis

VectorArray containing the generated collateral basis.

In addition, data contains the fields of the data dict returned by ei_greedy/deim.

eigs module¶

pymor.algorithms.eigs._arnoldi(A, E, l, b)[source]¶: Compute an Arnoldi factorization.

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

pymor.algorithms.eigs._qr_iteration(H, shifts)[source]¶: Perform the QR iteration.

pymor.algorithms.eigs.eigs(A, E=None, k=3, which='LM', b=None, l=None, maxiter=1000, tol=1e-13, imag_tol=1e-12, complex_pair_tol=1e-12, seed=0)[source]¶

Approximate a few eigenvalues of a linear Operator.

Computes k eigenvalues w with corresponding eigenvectors v 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 not None.

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, estimator=True, condition=False, error_norms=(), error_norm_names=None, 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 random Parameters.

Parameters

rom

The reduced Model.

fom

The high-dimensional 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 of reductor.reduce().

estimator

If True evaluate the error estimator of rom on the test Parameters.

condition

If True, compute the condition of the reduced system matrix for the given test Parameters (can only be specified if rom is an instance of StationaryModel and rom.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. If None, the name attributes of the given norms are used.

estimator_norm_index

When estimator is True and error_norms are specified, this is the index of the norm in error_norms w.r.t. which to compute the effectivity of the estimator.

custom

List of custom functions which are evaluated for each test parameter values and basis size. The functions must have the signature

def 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 y-axis to plot the computed custom values.

pool

If not None, the WorkerPool 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 high-dimensional solutions w.r.t. all given test Parameters, reduced basis dimensions and norms in error_norms. (Only present when error_norms has been specified.)

max_norms

Maxima of norms over the given test Parameters.

max_norm_mus

Parameters corresponding to max_norms.

errors

NumPy array of the norms of the model reduction errors w.r.t. all given test Parameters, reduced basis dimensions and norms in error_norms. (Only present when error_norms has been specified.)

max_errors

Maxima of errors over the given test Parameters.

max_error_mus

Parameters corresponding to max_errors.

rel_errors

errors divided by norms. (Only present when error_norms has been specified.)

max_rel_errors

Maxima of rel_errors over the given test Parameters.

max_rel_error_mus

Parameters corresponding to max_rel_errors.

error_norm_names

Names of the given error_norms. (Only present when error_norms has been specified.)

estimates

NumPy array of the model reduction error estimates w.r.t. all given test Parameters and reduced basis dimensions. (Only present when estimator is True.)

max_estimate

Maxima of estimates over the given test Parameters.

max_estimate_mus

Parameters corresponding to max_estimates.

effectivities

errors divided by estimates. (Only present when estimator is True and error_norms has been specified.)

min_effectivities

Minima of effectivities over the given test Parameters.

min_effectivity_mus

Parameters corresponding to min_effectivities.

max_effectivities

Maxima of effectivities over the given test Parameters.

max_effectivity_mus

Parameters corresponding to max_effectivities.

errors

NumPy array of the reduced system matrix conditions w.r.t. all given test Parameters and reduced basis dimensions. (Only present when conditions is True.)

max_conditions

Maxima of conditions over the given test Parameters.

max_condition_mus

Parameters corresponding to max_conditions.

custom_values

NumPy array of custom function evaluations w.r.t. all given test Parameters, reduced basis dimensions and functions in custom. (Only present when custom has been specified.)

max_custom_values

Maxima of custom_values over the given test Parameters.

max_custom_values_mus

Parameters corresponding to max_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 is True.)

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 in V using a generic iterative solver.

Parameters

op: The linear, non-parametric Operator to invert.
V: VectorArray of right-hand sides for the equation system.
initial_guess: VectorArray with the same length as V containing initial guesses for the solution. Some implementations of apply_inverse may ignore this parameter. If None a solver-dependent default is used.
options: The solver_options to use (see solver_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=1e-05, 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=1e-06, btol=1e-06, conlim=100000000.0, maxiter=None, show=False)[source]¶

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

pymor.algorithms.genericsolvers.solver_options(lgmres_tol=1e-05, lgmres_maxiter=1000, lgmres_inner_m=39, lgmres_outer_k=3, least_squares_lsmr_damp=0.0, least_squares_lsmr_atol=1e-06, least_squares_lsmr_btol=1e-06, 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=1e-06, least_squares_lsqr_btol=1e-06, 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: See scipy.sparse.linalg.lgmres.
lgmres_maxiter: See scipy.sparse.linalg.lgmres.
lgmres_inner_m: See scipy.sparse.linalg.lgmres.
lgmres_outer_k: See scipy.sparse.linalg.lgmres.
least_squares_lsmr_damp: See scipy.sparse.linalg.lsmr.
least_squares_lsmr_atol: See scipy.sparse.linalg.lsmr.
least_squares_lsmr_btol: See scipy.sparse.linalg.lsmr.
least_squares_lsmr_conlim: See scipy.sparse.linalg.lsmr.
least_squares_lsmr_maxiter: See scipy.sparse.linalg.lsmr.
least_squares_lsmr_show: See scipy.sparse.linalg.lsmr.
least_squares_lsqr_damp: See scipy.sparse.linalg.lsqr.
least_squares_lsqr_atol: See scipy.sparse.linalg.lsqr.
least_squares_lsqr_btol: See scipy.sparse.linalg.lsqr.
least_squares_lsqr_conlim: See scipy.sparse.linalg.lsqr.
least_squares_lsqr_iter_lim: See scipy.sparse.linalg.lsqr.
least_squares_lsqr_show: See scipy.sparse.linalg.lsqr.

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=1e-13, rtol=1e-13, offset=0, reiterate=True, reiteration_threshold=0.1, check=True, check_tol=0.001, copy=True)[source]¶

Orthonormalize a VectorArray using the modified Gram-Schmidt algorithm.

Parameters

A: The VectorArray which is to be orthonormalized.
product: The inner product Operator w.r.t. which to orthonormalize. If None, 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 the offset + 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 is True, re-orthonormalize 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 resulting VectorArray is really orthonormal.
check_tol: Tolerance for the check.
copy: If True, create a copy of A instead of modifying A in-place.

Returns

Q: The orthonormalized VectorArray.
R: The upper-triangular/trapezoidal matrix (if compute_R is True).

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 Gram-Schmidt 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. If None, 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 is True, re-orthonormalize 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 resulting VectorArray is really orthonormal.
check_tol: Tolerance for the check.
copy: If True, create a copy of V and W instead of modifying V and W in-place.

Returns

The biorthonormalized VectorArrays.

greedy module¶

class pymor.algorithms.greedy.RBSurrogate(fom, reductor, use_estimator, error_norm, extension_params, pool)[source]¶

Bases: pymor.algorithms.greedy.WeakGreedySurrogate

Surrogate for the weak_greedy error used in rb_greedy.

Not intended to be used directly.

Methods

`RBSurrogate`	`evaluate`, `extend`
`BasicObject`	`disable_logging`, `enable_logging`

Attributes

BasicObject

logger, logging_disabled, name, uid

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 a RemoteObject.
return_all_values: See below.

Returns

If return_all_values is True, an NumPy array of the estimated errors. If return_all_values is False, 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

`WeakGreedySurrogate`	`evaluate`, `extend`
`BasicObject`	`disable_logging`, `enable_logging`

Attributes

BasicObject

logger, logging_disabled, name, uid

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 a RemoteObject.
return_all_values: See below.

Returns

If return_all_values is True, an NumPy array of the estimated errors. If return_all_values is False, the maximum estimated error as first return value and the corresponding parameter as second return value.

pymor.algorithms.greedy.greedy(*args, **kwargs)[source]¶

pymor.algorithms.greedy.rb_greedy(fom, reductor, training_set, use_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 a reduce method, such that reductor.reduce() yields the reduced model, and an exted_basis method, such that reductor.extend_basis(U, copy_U=False, **extension_params) extends the current reduced basis by the vectors contained in U. For an example see CoerciveRBReductor.
training_set: The training set of Parameters on which to perform the greedy search.
use_estimator: If False, exactly compute the model reduction error by also computing the solution of fom for all parameter values of the training set. This is mainly useful when no estimator for the model reduction error is available.
error_norm: If use_estimator is False, use this function to calculate the norm of the error. If None, 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 the reductor.extend_basis method. If None, 'gram_schmidt' basis extension will be used as a default for stationary problems (fom.solve returns VectorArrays 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 after max_extensions extension steps.
pool: If not None, a WorkerPool to use for parallelization. Parallelization needs to be supported by surrogate.

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

`DistHAPODTree`	`children`
`Tree`	`is_leaf`
`BasicObject`	`disable_logging`, `enable_logging`

Attributes

`Tree`	`depth`, `root`
`BasicObject`	`logger`, `logging_disabled`, `name`, `uid`

class pymor.algorithms.hapod.FakeExecutor[source]¶

Bases: object

Methods

FakeExecutor

submit

class pymor.algorithms.hapod.IncHAPODTree(steps)[source]¶

Bases: pymor.algorithms.hapod.Tree

Methods

`IncHAPODTree`	`children`
`Tree`	`is_leaf`
`BasicObject`	`disable_logging`, `enable_logging`

Attributes

`Tree`	`depth`, `root`
`BasicObject`	`logger`, `logging_disabled`, `name`, `uid`

class pymor.algorithms.hapod.LifoExecutor(executor, max_workers=None)[source]¶

Bases: object

Methods

LifoExecutor

done_callback, run_task, submit

class pymor.algorithms.hapod.Tree[source]¶

Bases: pymor.core.base.BasicObject

A rooted tree.

Methods

`Tree`	`children`, `is_leaf`
`BasicObject`	`disable_logging`, `enable_logging`

Attributes

`Tree`	`depth`, `root`
`BasicObject`	`logger`, `logging_disabled`, `name`, `uid`

pymor.algorithms.hapod.default_pod_method(U, eps, is_root_node, product)[source]¶

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 l2-mean 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, a concurrent.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 l2-mean 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, a concurrent.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 node node an l2 truncation error tolerance for the local pod based on the number of input vectors num_vecs and the total number of snapshot vectors below the given node snap_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 the VectorArray U w.r.t. the given inner product product and the l2 error tolerance eps. root_node is set to True when the POD is computed at the root of the tree.
executor: If not None, a concurrent.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 l2-mean 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, a concurrent.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 l2-mean 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, a concurrent.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.std_local_eps(tree, eps, omega, pod_on_leafs=True)[source]¶

image module¶

class pymor.algorithms.image.CollectOperatorRangeRules(source, image, extends)[source]¶

Bases: pymor.algorithms.rules.RuleTable

RuleTable for the estimate_image algorithm.

Methods

`RuleTable`	`append_rule`, `apply`, `apply_children`, `breakpoint_for_name`, `breakpoint_for_obj`, `get_children`, `insert_rule`, `replace_children`
`BasicObject`	`disable_logging`, `enable_logging`

Attributes

`CollectOperatorRangeRules`	`action_apply_operator`, `action_Concatenation`, `action_EmpiricalInterpolatedOperator`, `action_recurse`, `rules`
`BasicObject`	`logger`, `logging_disabled`, `name`, `uid`

class pymor.algorithms.image.CollectVectorRangeRules(image)[source]¶

Bases: pymor.algorithms.rules.RuleTable

RuleTable for the estimate_image algorithm.

Methods

`RuleTable`	`append_rule`, `apply`, `apply_children`, `breakpoint_for_name`, `breakpoint_for_obj`, `get_children`, `insert_rule`, `replace_children`
`BasicObject`	`disable_logging`, `enable_logging`

Attributes

`CollectVectorRangeRules`	`action_as_range_array`, `action_Concatenation`, `action_recurse`, `action_VectorArray`, `rules`
`BasicObject`	`logger`, `logging_disabled`, `name`, `uid`

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 of Operators with common source and range, and let vectors be a list of VectorArrays or vector-like Operators in the range of these operators. Given a VectorArray domain of vectors in the source of the operators, this algorithms determines a VectorArray image of range vectors such that the linear span of image contains:

op.apply(U, mu=mu) for all operators op in operators, for all possible Parameters mu and for all VectorArrays U contained in the linear span of domain,
U for all VectorArrays in vectors,
v.as_range_array(mu) for all Operators in vectors and all possible Parameters mu.

The algorithm will try to choose image as small as possible. However, no optimality is guaranteed. The image estimation algorithm is specified by CollectOperatorRangeRules and CollectVectorRangeRules.

Parameters

operators: See above.
vectors: See above.
domain: See above. If None, an empty domain VectorArray is assumed.
extends: For some operators, e.g. EmpiricalInterpolatedOperator, as well as for all elements of vectors, image is estimated independently from the choice of domain. If extends is True, 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 the gram_schmidt algorithm.
product: Inner product Operator w.r.t. which to orthonormalize.
riesz_representatives: If True, compute Riesz representatives of the vectors in image before orthonormalizing (useful for dual norm computation when the range of the operators 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 calls estimate_image individually for each vector of domain.

As a result, the vectors in the returned image VectorArray will be ordered by the domain vector they correspond to (starting with vectors which correspond to the elements of vectors and to Operators for which the image is estimated independently from domain).

This function also returns an image_dims list, such that the first image_dims[i+1] vectors of image correspond to the first i vectors of domain (the first image_dims[0] vectors correspond to vectors and to Operators 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 after estimate_image_hierarchical has been called, and estimate_image_hierarchical shall be called again for the extended domain array, extends can be set to (image, image_dims), where image, image_dims are the return values of the last estimate_image_hierarchical call. The old domain vectors will then be skipped during computation and image, image_dims will be modified in-place.
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

A: Real Operator A.
E: Real Operator E.
b: Real vector-like operator (if trans is False) or functional (if trans is True).
sigma: Sequence of interpolation points (closed under conjugation).
trans: Boolean, see above.

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 from B.source, if trans == False, or: B.range, if trans == 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 using pymor.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

`RuleTable`	`append_rule`, `apply`, `apply_children`, `breakpoint_for_name`, `breakpoint_for_obj`, `get_children`, `insert_rule`, `replace_children`
`BasicObject`	`disable_logging`, `enable_logging`

Attributes

AssembleLincombRules

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

BasicObject

logger, logging_disabled, name, uid

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 sum

c_1*O_1 + ... + c_N*O_N

where O_i are Operators and c_i scalar coefficients.

This function is called in the assemble method of LincombOperator and is not intended to be used directly. The assembled Operator is expected to no longer be a LincombOperator nor should it contain any LincombOperators. If an assembly of the given linear combination is not possible, None is returned. The special case of a LincombOperator with a single operator (i.e. a scaled Operator) is allowed as assemble_lincomb implements apply_inverse for this special case.

To form the linear combination of backend Operators (containing actual matrix data), _assemble_lincomb will be called on the first Operator 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, otherwise None.

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: Real-valued 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 point starting_point.
initial_value: Value of f in the point starting_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=1e-10, 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: See projection_shifts_init.
projection_shifts_init_seed: See projection_shifts_init.

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 low-rank ADI iteration using Galerkin projection on spaces spanned by LR-ADI iterates.

See [PK16], pp. 92-95.

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 low-rank ADI iteration using Galerkin projection on spaces spanned by LR-ADI iterates.

See [PK16], pp. 92-95.

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 low-rank solution of a Lyapunov equation.

See pymor.algorithms.lyapunov.solve_lyap_lrcf for a general description.

This function uses the low-rank ADI iteration as described in Algorithm 4.3 in [PK16].

Parameters

A: The non-parametric Operator A.
E: The non-parametric Operator E or None.
B: The operator B as a VectorArray from A.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: Low-rank Cholesky factor of the Lyapunov equation solution, VectorArray from A.source.

lrradi module¶

pymor.algorithms.lrradi.hamiltonian_shifts(A, E, B, R, K, Z, shift_options)[source]¶

Compute further shift parameters for low-rank 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. 318-321.

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 low-rank 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. 318-321.

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=1e-10, 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: See hamiltonian_shifts_init.
hamiltonian_shifts_init_seed: See hamiltonian_shifts_init.
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, S=None, trans=False, options=None)[source]¶

Compute an approximate low-rank 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 or None.
B: The operator B as a VectorArray from A.source.
C: The operator C as a VectorArray from A.source.
R: The operator R as a 2D NumPy array or None.
S: The operator S as a VectorArray from A.source or None.
trans: Whether the first Operator in the Riccati equation is transposed.
options: The solver options to use. (see ricc_lrcf_solver_options)

Returns

Z: Low-rank Cholesky factor of the Riccati equation solution, VectorArray from A.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) continuous-time algebraic Lyapunov equation:

if trans is False and E is None:

\[A X + X A^T + B B^T = 0,\]
if trans is False and E is an Operator:

\[A X E^T + E X A^T + B B^T = 0,\]
if trans is True and E is None:

\[A^T X + X A + B^T B = 0,\]
if trans is True and E is an Operator:

\[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 (see pymor.bindings.pymess.solve_lyap_dense)
slycot (see pymor.bindings.slycot.solve_lyap_dense)
scipy (see pymor.bindings.scipy.solve_lyap_dense)

Parameters

A

The operator A as a 2D NumPy array.

E

The operator E as a 2D NumPy array or None.

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 low-rank solution of a Lyapunov equation.

Returns a low-rank Cholesky factor \(Z\) such that \(Z Z^T\) approximates the solution \(X\) of a (generalized) continuous-time algebraic Lyapunov equation:

if trans is False and E is None:

\[A X + X A^T + B B^T = 0,\]
if trans is False and E is an Operator:

\[A X E^T + E X A^T + B B^T = 0,\]
if trans is True and E is None:

\[A^T X + X A + B^T B = 0,\]
if trans is True and E is an Operator:

\[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 half-plane. Operator B needs to be given as a VectorArray from A.source, and for large-scale problems, we assume len(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)
1. pymess (see pymor.bindings.pymess.solve_lyap_lrcf),
2. lradi (see pymor.algorithms.lradi.solve_lyap_lrcf),
for dense problems (smaller than mat_eqn_sparse_min_size)
1. pymess (see pymor.bindings.pymess.solve_lyap_lrcf),
2. slycot (see pymor.bindings.slycot.solve_lyap_lrcf),
3. scipy (see pymor.bindings.scipy.solve_lyap_lrcf).

Parameters

A

The non-parametric Operator A.

E

The non-parametric Operator E or None.

B

The operator B as a VectorArray from A.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: Low-rank Cholesky factor of the Lyapunov equation solution, VectorArray from A.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=1e-07, 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 the jacobian interface method.
rhs: VectorArray of length 1 containing the vector V.
initial_guess: If not None, a VectorArray of length 1 containing an initial guess for the solution U.
mu: The parameter values for which to solve the equation.
range_product: The inner product Operator on operator.range with which the norm of the resiudal is computed. If None, the Euclidean inner product is used.
source_product: The inner product Operator on operator.source with which the norm of the solution and update vectors is computed. If None, 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 the armijo 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 last stagnation_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 a VectorArray of the intermediate approximations of U after each iteration.
return_residuals: If True, return a VectorArray 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=1e-07, atol=0.0, l2_err=0.0, method='method_of_snapshots', orth_tol=1e-10)[source]¶

Proper orthogonal decomposition of A.

Viewing the VectorArray A as a A.dim x len(A) matrix, the return values of this method are the VectorArray of left singular vectors and a NumPy array of singular values of the singular value decomposition of A, where the inner product on R^(dim(A)) is given by product 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 first modes 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 l2-approximation 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 n-th 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: One-dimensional 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

`RuleTable`	`append_rule`, `apply`, `apply_children`, `breakpoint_for_name`, `breakpoint_for_obj`, `get_children`, `insert_rule`, `replace_children`
`BasicObject`	`disable_logging`, `enable_logging`

Attributes

`PreAssembleRules`	`action_AdjointOperator`, `action_assemble`, `action_identity`, `action_recurse`, `action_recurse_and_assemble`, `rules`
`BasicObject`	`logger`, `logging_disabled`, `name`, `uid`

pymor.algorithms.preassemble.preassemble(obj)[source]¶

Preassemble non-parametric operators.

If obj is a non-parametric Operator, return obj.assemble() otherwise return obj. Recursively replaces children of obj.

projection module¶

class pymor.algorithms.projection.ProjectRules(range_basis, source_basis)[source]¶

Bases: pymor.algorithms.rules.RuleTable

RuleTable for the project algorithm.

Methods

`RuleTable`	`append_rule`, `apply`, `apply_children`, `breakpoint_for_name`, `breakpoint_for_obj`, `get_children`, `insert_rule`, `replace_children`
`BasicObject`	`disable_logging`, `enable_logging`

Attributes

`ProjectRules`	`action_AdjointOperator`, `action_AffineOperator`, `action_apply_basis`, `action_BlockOperatorBase`, `action_Concatenation`, `action_ConstantOperator`, `action_EmpiricalInterpolatedOperator`, `action_LincombOperator`, `action_no_bases`, `action_SelectionOperator`, `action_ZeroOperator`, `rules`
`BasicObject`	`logger`, `logging_disabled`, `name`, `uid`

class pymor.algorithms.projection.ProjectToSubbasisRules(dim_range, dim_source)[source]¶

Bases: pymor.algorithms.rules.RuleTable

RuleTable for the project_to_subbasis algorithm.

Methods

`RuleTable`	`append_rule`, `apply`, `apply_children`, `breakpoint_for_name`, `breakpoint_for_obj`, `get_children`, `insert_rule`, `replace_children`
`BasicObject`	`disable_logging`, `enable_logging`

Attributes

`ProjectToSubbasisRules`	`action_ConstantOperator`, `action_IdentityOperator`, `action_NumpyMatrixOperator`, `action_ProjectedEmpiciralInterpolatedOperator`, `action_ProjectedOperator`, `action_recurse`, `rules`
`BasicObject`	`logger`, `logging_disabled`, `name`, `uid`

pymor.algorithms.projection.project(op, range_basis, source_basis, product=None)[source]¶

Petrov-Galerkin projection of a given Operator.

Given an inner product ( ⋅, ⋅), source vectors b_1, ..., b_N and range vectors c_1, ..., c_M, the projection op_proj of op is defined by

[ op_proj(e_j) ]_i = ( c_i, op(b_j) )

for all i,j, where e_j denotes the j-th canonical basis vector of R^N.

In particular, if the c_i are orthonormal w.r.t. the given product, then op_proj is the coordinate representation w.r.t. the b_i/c_i bases of the restriction of op to span(b_i) concatenated with the orthogonal projection onto span(c_i).

From another point of view, if op is viewed as a bilinear form (see apply2) and ( ⋅, ⋅ ) is the Euclidean inner product, then op_proj represents the matrix of the bilinear form restricted to span(b_i) / span(c_i) (w.r.t. the b_i/c_i bases).

How the projection is realized will depend on the given Operator. While a projected NumpyMatrixOperator will again be a NumpyMatrixOperator, only a generic ProjectedOperator can be returned in general. The exact algorithm is specified in ProjectRules.

Parameters

range_basis: The vectors c_1, ..., c_M as a VectorArray. If None, no projection in the range space is performed.
source_basis: The vectors b_1, ..., b_N as a VectorArray or None. If None, no restriction of the source space is performed.
product: An Operator representing the inner product. If None, 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 upper-left (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=1e-15, 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 the VectorArray 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 of A is given by range_product and the inner product of the source of A is given by source_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 the VectorArray Q whose vectors form an approximate orthonomal basis for the range of A.

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.
q: The number of power iterations.
l: The block size of the normalized power iterations.
iscomplex: If True, the random vectors are chosen complex.

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, S=None, trans=False, options=None, default_dense_solver_backend='pymess')[source]¶

Compute an approximate low-rank solution of a positive Riccati equation.

Returns a low-rank Cholesky factor \(Z\) such that \(Z Z^T\) approximates the solution \(X\) of a (generalized) positive continuous-time algebraic Riccati equation:

if trans is False

\[A X E^T + E X A^T + (E X C^T + S) R^{-1} (E X C^T + S)^T + B B^T = 0.\]
if trans is True

\[A^T X E + E^T X A + (E^T X B + S) R^{-1} (E^T X B + S)^T + C^T C = 0.\]

If E is None, it is taken to be identity, and similarly for R. If S is None, it is taken to be zero.

If the solver is not specified using the options argument, a solver backend is chosen based on availability in the following order:

pymess (see pymor.bindings.pymess.solve_pos_ricc_lrcf),
slycot (see pymor.bindings.slycot.solve_pos_ricc_lrcf),
scipy (see pymor.bindings.scipy.solve_pos_ricc_lrcf).

Parameters

A

The non-parametric Operator A.

E

The non-parametric Operator E or None.

B

The operator B as a VectorArray from A.source.

C

The operator C as a VectorArray from A.source.

R

The operator R as a 2D NumPy array or None.

S

The operator S as a VectorArray from A.source or None.

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: Low-rank Cholesky factor of the positive Riccati equation solution, VectorArray from A.source.

Defaults

default_dense_solver_backend (see pymor.core.defaults)

pymor.algorithms.riccati.solve_ricc_lrcf(A, E, B, C, R=None, S=None, trans=False, options=None, default_sparse_solver_backend='pymess', default_dense_solver_backend='pymess')[source]¶

Compute an approximate low-rank solution of a Riccati equation.

Returns a low-rank Cholesky factor \(Z\) such that \(Z Z^T\) approximates the solution \(X\) of a (generalized) continuous-time algebraic Riccati equation:

if trans is False

\[A X E^T + E X A^T - (E X C^T + S) R^{-1} (E X C^T + S)^T + B B^T = 0.\]
if trans is True

\[A^T X E + E^T X A - (E^T X B + S) R^{-1} (E^T X B + S)^T + C^T C = 0.\]

If E is None, it is taken to be identity, and similarly for R. If S is None, it is taken to be zero.

We assume:

A and E are real Operators,
B, C and S are real VectorArrays from A.source,
R is a real NumPy array,
(E, A, B, C) is stabilizable and detectable,
R is symmetric positive definite, and
\(B B^T - S R^{-1} S^T\) (\(C^T C - S R^{-1} S^T\)) is positive semi-definite trans is False (True).

For large-scale problems, we additionally assume that len(B) and len(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)
1. pymess (see pymor.bindings.pymess.solve_ricc_lrcf),
2. lrradi (see pymor.algorithms.lrradi.solve_ricc_lrcf),
for dense problems (smaller than mat_eqn_sparse_min_size)
1. pymess (see pymor.bindings.pymess.solve_ricc_lrcf),
2. slycot (see pymor.bindings.slycot.solve_ricc_lrcf),
3. scipy (see pymor.bindings.scipy.solve_ricc_lrcf).

Parameters

A

The non-parametric Operator A.

E

The non-parametric Operator E or None.

B

The operator B as a VectorArray from A.source.

C

The operator C as a VectorArray from A.source.

R

The operator R as a 2D NumPy array or None.

S

The operator S as a VectorArray from A.source or None.

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: Low-rank Cholesky factor of the Riccati equation solution, VectorArray from A.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 of rules, stored in the rules attributes, which can be applied to given objects.

A new table is created by subclassing RuleTable and defining new methods which are decorated with match_class, match_generic or another rule subclass. The order of the method definitions determines the order in which the defined rules are applied.

Parameters

use_caching: If True, cache results of apply.

Methods

`RuleTable`	`append_rule`, `apply`, `apply_children`, `breakpoint_for_name`, `breakpoint_for_obj`, `get_children`, `insert_rule`, `replace_children`
`BasicObject`	`disable_logging`, `enable_logging`

Attributes

`RuleTable`	`rules`
`BasicObject`	`logger`, `logging_disabled`, `name`, `uid`

rules¶: list of all defined rules.

__repr__()[source]¶: Return repr(self).

apply(obj)[source]¶

Sequentially apply rules to given object.

This method iterates over all rules of the given RuleTable. For each rule, it is checked if it matches the given object. If False, the next rule in the table is considered. If True the corresponding action is executed with obj as parameter. If execution of action raises RuleNotMatchingError, the rule is considered as not matching, and execution continues with evaluation of the next rule. Otherwise, execution is stopped and the return value of rule.action is returned to the caller.

If no rule matches, a NoMatchingRuleError 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 the children parameter or automatically inferred by the get_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:

a is an Operator.

a is a mapping and each of its values is either an Operator or None.

a is an iterable and each of its elements is either an Operator or None.

replace_children(obj, children=None)[source]¶

Replace children of object according to rule table.

Same as apply_children, but additionally calls obj.with_ to replace the children of obj with the result of the corresponding apply call.

class pymor.algorithms.rules.RuleTableMeta(name, parents, dct)[source]¶

Bases: pymor.core.base.UberMeta

Meta class for RuleTable.

Methods

`ABCMeta`	`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 to 1.

__repr__()[source]¶: Return repr(self).

__str__()¶: Return repr(self).

pymor.algorithms.rules.format_rules(rules)[source]¶

class pymor.algorithms.rules.match_always(action)[source]¶

Bases: pymor.algorithms.rules.rule

rule that always matches.

Methods

rule

matches

Attributes

`match_always`	`condition_type`
`rule`	`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.

Methods

rule

matches

Attributes

`match_class`	`condition_type`
`rule`	`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.

Methods

rule

matches

Attributes

`match_class_all`	`condition_type`
`rule`	`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.

Methods

rule

matches

Attributes

`match_class_any`	`condition_type`
`rule`	`action`, `action_description`, `condition_description`, `source`

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

Bases: pymor.algorithms.rules.rule

Methods

rule

matches

Attributes

rule

action, action_description, condition_description, condition_type, source

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.

Methods

rule

matches

Attributes

`match_generic`	`condition_type`
`rule`	`action`, `action_description`, `condition_description`, `source`

pymor.algorithms.rules.print_children(obj)[source]¶

pymor.algorithms.rules.print_rules(rules)[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 rule matches. Matching conditions are specified by subclassing and overriding the matches method.

If an action is decorated by multiple rules, all these rules must match for the action to apply.

Methods

rule

matches

Attributes

rule

action, action_description, condition_description, condition_type, source

action¶: Method to call in case the rule matches.

__call__(action)[source]¶: Call self as a function.

__repr__()[source]¶: Return repr(self).

matches(obj)[source]¶: Returns True if given object matches the condition.

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=1e-10, imagtol=1e-06, conjtol=1e-08, dorqitol=0.0001, rqitol=1e-10, 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 or None.

B

The operator B as a VectorArray from A.source.

C

The operator C as a VectorArray from A.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 two-sided Rayleigh quotient iteration is executed.

rqitol

Tolerance for the relative change of a pole in the two-sided 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 two-sided 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)

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=1e-07, atol=0.0, l2_err=0.0)[source]¶

SVD of a VectorArray using the method of snapshots.

Viewing the VectorArray A as a A.dim x len(A) matrix, the return value of this method is the singular value decomposition of A, where the inner product on R^(dim(A)) is given by product and the inner product on R^(len(A)) is the Euclidean inner product.

Warning

The left singular vectors may not be numerically orthonormal for ill-conditioned 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 first modes 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 l2-approximation 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 n-th singular value.

Returns

U: VectorArray of left singular vectors.
s: One-dimensional 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=4e-08, atol=0.0, l2_err=0.0)[source]¶

SVD of a VectorArray using Gram-Schmidt orthogonalization.

Viewing the VectorArray A as a A.dim x len(A) matrix, the return value of this method is the singular value decomposition of A, where the inner product on R^(dim(A)) is given by product 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 first modes 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 l2-approximation 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 n-th singular value.

Returns

U: VectorArray of left singular vectors.
s: One-dimensional 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 a NumPy array using to_matrix.
E: Real Operator or None (then assumed to be the identity).
Er: Real Operator or None (then assumed to be the identity). It is converted into a NumPy array using to_matrix.
B: Real Operator or None.
Br: Real Operator or None. It is assumed that Br.range.from_numpy is implemented.
C: Real Operator or None.
Cr: Real Operator or None. It is assumed that Cr.source.from_numpy is implemented.

Returns

V: Returned if B and Br are given, VectorArray from A.source.
W: Returned if C and Cr are given, VectorArray from A.source.

Raises

ValueError: If V and W cannot be returned.

timestepping module¶

This module provides generic time-stepping 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 time-stepping are provided. The TimeStepper defines a common interface that has to be fulfilled by the time-steppers used by InstationaryModel. The classes ExplicitEulerTimeStepper and ImplicitEulerTimeStepper encapsulate explicit_euler and implicit_euler to provide this interface.

class pymor.algorithms.timestepping.ExplicitEulerTimeStepper(*args, **kwargs)[source]¶

Bases: pymor.algorithms.timestepping.TimeStepper

Explicit Euler time-stepper.

Solves equations of the form

M * d_t u + A(u, mu, t) = F(mu, t).

Parameters

nt: The number of time-steps the time-stepper will perform.

Methods

`ExplicitEulerTimeStepper`	`solve`
`ImmutableObject`	`with_`, `__setattr__`
`BasicObject`	`disable_logging`, `enable_logging`

Attributes

BasicObject

logger, logging_disabled, name, uid

solve(initial_time, end_time, initial_data, operator, rhs=None, mass=None, mu=None, num_values=None)[source]¶

Apply time-stepper to the equation

M * d_t u + A(u, mu, t) = F(mu, t).

Parameters

initial_time: The time at which to begin time-stepping.
end_time: The time until which to perform time-stepping.
initial_data: The solution vector at initial_time.
operator: The Operator A.
rhs: The right-hand side F (either VectorArray of length 1 or Operator with source.dim == 1). If None, zero right-hand side is assumed.
mass: The Operator M. If None, the identity operator is assumed.
mu: Parameter values for which operator and rhs are evaluated. The current time is added to mu with key t.
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(*args, **kwargs)[source]¶

Bases: pymor.algorithms.timestepping.TimeStepper

Implict Euler time-stepper.

Solves equations of the form

M * d_t u + A(u, mu, t) = F(mu, t).

Parameters

nt: The number of time-steps the time-stepper will perform.
solver_options: The solver_options used to invert M + dt*A. The special values 'mass' and 'operator' are recognized, in which case the solver_options of M (resp. A) are used.

Methods

`ImplicitEulerTimeStepper`	`solve`
`ImmutableObject`	`with_`, `__setattr__`
`BasicObject`	`disable_logging`, `enable_logging`

Attributes

BasicObject

logger, logging_disabled, name, uid

solve(initial_time, end_time, initial_data, operator, rhs=None, mass=None, mu=None, num_values=None)[source]¶

Apply time-stepper to the equation

M * d_t u + A(u, mu, t) = F(mu, t).

Parameters

initial_time: The time at which to begin time-stepping.
end_time: The time until which to perform time-stepping.
initial_data: The solution vector at initial_time.
operator: The Operator A.
rhs: The right-hand side F (either VectorArray of length 1 or Operator with source.dim == 1). If None, zero right-hand side is assumed.
mass: The Operator M. If None, the identity operator is assumed.
mu: Parameter values for which operator and rhs are evaluated. The current time is added to mu with key t.
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(*args, **kwargs)[source]¶

Bases: pymor.core.base.ImmutableObject

Interface for time-stepping algorithms.

Algorithms implementing this interface solve time-dependent problems of the form

M * d_t u + A(u, mu, t) = F(mu, t).

Time-steppers used by InstationaryModel have to fulfill this interface.

Methods

`TimeStepper`	`solve`
`ImmutableObject`	`with_`, `__setattr__`
`BasicObject`	`disable_logging`, `enable_logging`

Attributes

BasicObject

logger, logging_disabled, name, uid

abstract solve(initial_time, end_time, initial_data, operator, rhs=None, mass=None, mu=None, num_values=None)[source]¶

Apply time-stepper to the equation

M * d_t u + A(u, mu, t) = F(mu, t).

Parameters

initial_time: The time at which to begin time-stepping.
end_time: The time until which to perform time-stepping.
initial_data: The solution vector at initial_time.
operator: The Operator A.
rhs: The right-hand side F (either VectorArray of length 1 or Operator with source.dim == 1). If None, zero right-hand side is assumed.
mass: The Operator M. If None, the identity operator is assumed.
mu: Parameter values for which operator and rhs are evaluated. The current time is added to mu with key t.
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.

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

`RuleTable`	`append_rule`, `apply`, `apply_children`, `breakpoint_for_name`, `breakpoint_for_obj`, `get_children`, `insert_rule`, `replace_children`
`BasicObject`	`disable_logging`, `enable_logging`

Attributes

`ToMatrixRules`	`action_AdjointOperator`, `action_BlockOperator`, `action_ComponentProjection`, `action_Concatenation`, `action_IdentityOperator`, `action_LincombOperator`, `action_LowRankOperator`, `action_LowRankUpdatedOperator`, `action_NumpyMatrixOperator`, `action_VectorArrayOperator`, `action_ZeroOperator`, `rules`
`BasicObject`	`logger`, `logging_disabled`, `name`, `uid`

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 appropriate SciPy spmatrix. If None, a choice between dense and sparse format is automatically made.
mu: The parameter values for which to convert op.

Returns

res: The matrix equivalent to op.