pymor.algorithms package

Submodules

adaptivegreedy module


class pymor.algorithms.adaptivegreedy.AdaptiveSampleSet(parameter_space)[source]

Bases: pymor.core.interfaces.BasicInterface

An adaptive parameter sample set.

Used by adaptive_weak_greedy.


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=None, 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. If None, the parameter space of fom is used.
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.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.

error module


pymor.algorithms.error.reduction_error_analysis(rom, fom, reductor, test_mus=10, basis_sizes=0, random_seed=None, 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
Either a list of Parameters to compute the errors for, or the number of parameters which are sampled randomly from parameter_space (if given) or rom.parameter_space.
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().
random_seed
If test_mus is a number, use this value as random seed for drawing the Parameters.
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 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, rhs, 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 rhs using a generic iterative solver.

Parameters

op
The linear, non-parametric Operator to invert.
rhs
VectorArray of right-hand sides for the equation system.
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.

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.interfaces.BasicInterface

Surrogate for the approximation error in weak_greedy.

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 each training set Parameter. 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

Attributes

Tree depth, root
BasicInterface 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

Attributes

Tree depth, root
BasicInterface 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.interfaces.BasicInterface

A rooted tree.

Attributes

Tree depth, root
BasicInterface 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.

Attributes

CollectOperatorRangeRules action_apply_operator, action_Concatenation, action_EmpiricalInterpolatedOperator, action_recurse, rules
BasicInterface logger, logging_disabled, name, uid

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

Bases: pymor.algorithms.rules.RuleTable

RuleTable for the estimate_image algorithm.

Attributes

CollectVectorRangeRules action_as_range_array, action_Concatenation, action_recurse, action_VectorArray, rules
BasicInterface 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

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
BasicInterface 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.

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='scipy')[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:

  1. pymess (see pymor.bindings.pymess.solve_lyap_dense)
  2. slycot (see pymor.bindings.slycot.solve_lyap_dense)
  3. 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='lradi', default_dense_solver_backend='scipy')[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:

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, error_norm=None, least_squares=False, miniter=0, maxiter=100, rtol=-1.0, atol=-1.0, relax=1.0, stagnation_window=3, stagnation_threshold=inf, return_stages=False, return_residuals=False)[source]

Basic 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 for which to solve the equation.
error_norm
The norm with which the norm of the residual is computed. If None, the Euclidean norm 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.
rtol
Finish when the residual norm has been reduced by this factor relative to the norm of the initial residual.
atol
Finish when the residual norm is below this threshold.
relax
Relaxation factor for Newton updates.
stagnation_window
Finish when the residual norm has not been reduced by a factor of stagnation_threshold during the last stagnation_window iterations.
stagnation_threshold
See stagnation_window.
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:

error_sequence:NumPy array containing the residual norms after each iteration.
stages:See return_stages.
residuals:See return_residuals.

Raises

NewtonError
Raised if the Netwon 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=4e-08, 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

Attributes

PreAssembleRules action_AdjointOperator, action_assemble, action_identity, action_recurse, action_recurse_and_assemble, rules
BasicInterface 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.

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
BasicInterface 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.

Attributes

ProjectToSubbasisRules action_ConstantOperator, action_IdentityOperator, action_NumpyMatrixOperator, action_ProjectedEmpiciralInterpolatedOperator, action_ProjectedOperator, action_recurse, rules
BasicInterface 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='scipy')[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:

  1. pymess (see pymor.bindings.pymess.solve_pos_ricc_lrcf),
  2. slycot (see pymor.bindings.slycot.solve_pos_ricc_lrcf),
  3. 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='lrradi', default_dense_solver_backend='scipy')[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:

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.interfaces.BasicInterface

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.
rules

list of all defined rules.

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 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:

  1. a is an Operator.
  2. a is a mapping and each of its values is either an Operator or None.
  3. 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, bases, namespace)[source]

Bases: pymor.core.interfaces.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).


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]

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.

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=4e-08, 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 TimeStepperInterface 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(nt)[source]

Bases: pymor.algorithms.timestepping.TimeStepperInterface

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.
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 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(nt, solver_options='operator')[source]

Bases: pymor.algorithms.timestepping.TimeStepperInterface

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.
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 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.TimeStepperInterface[source]

Bases: pymor.core.interfaces.ImmutableInterface

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.

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 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

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
BasicInterface 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 for which to convert op.

Returns

res
The matrix equivalent to op.