Source code for pymor.bindings.scipy

# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2019 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)


import numpy as np
from packaging.version import Version
import scipy.version
from scipy.linalg import solve, solve_continuous_lyapunov, solve_continuous_are
from scipy.sparse.linalg import bicgstab, spsolve, splu, spilu, lgmres, lsqr, LinearOperator

from pymor.algorithms.lyapunov import _solve_lyap_lrcf_check_args, _solve_lyap_dense_check_args, _chol
from pymor.algorithms.riccati import _solve_ricc_check_args
from pymor.algorithms.genericsolvers import _parse_options
from pymor.algorithms.to_matrix import to_matrix
from pymor.core.config import config
from pymor.core.defaults import defaults
from pymor.core.exceptions import InversionError
from pymor.core.logger import getLogger
from pymor.operators.numpy import NumpyMatrixOperator


[docs]@defaults('bicgstab_tol', 'bicgstab_maxiter', 'spilu_drop_tol', 'spilu_fill_factor', 'spilu_drop_rule', 'spilu_permc_spec', 'spsolve_permc_spec', 'spsolve_keep_factorization', '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', sid_ignore=('least_squares_lsmr_show', 'least_squares_lsqr_show')) def solver_options(bicgstab_tol=1e-15, bicgstab_maxiter=None, spilu_drop_tol=1e-4, spilu_fill_factor=10, spilu_drop_rule=None, spilu_permc_spec='COLAMD', spsolve_permc_spec='COLAMD', spsolve_keep_factorization=True, lgmres_tol=1e-5, lgmres_maxiter=1000, lgmres_inner_m=39, lgmres_outer_k=3, least_squares_lsmr_damp=0.0, least_squares_lsmr_atol=1e-6, least_squares_lsmr_btol=1e-6, least_squares_lsmr_conlim=1e8, least_squares_lsmr_maxiter=None, least_squares_lsmr_show=False, least_squares_lsqr_damp=0.0, least_squares_lsqr_atol=1e-6, least_squares_lsqr_btol=1e-6, least_squares_lsqr_conlim=1e8, least_squares_lsqr_iter_lim=None, least_squares_lsqr_show=False): """Returns available solvers with default |solver_options| for the SciPy backend. Parameters ---------- bicgstab_tol See :func:`scipy.sparse.linalg.bicgstab`. bicgstab_maxiter See :func:`scipy.sparse.linalg.bicgstab`. spilu_drop_tol See :func:`scipy.sparse.linalg.spilu`. spilu_fill_factor See :func:`scipy.sparse.linalg.spilu`. spilu_drop_rule See :func:`scipy.sparse.linalg.spilu`. spilu_permc_spec See :func:`scipy.sparse.linalg.spilu`. spsolve_permc_spec See :func:`scipy.sparse.linalg.spsolve`. spsolve_keep_factorization See :func:`scipy.sparse.linalg.spsolve`. lgmres_tol See :func:`scipy.sparse.linalg.lgmres`. lgmres_maxiter See :func:`scipy.sparse.linalg.lgmres`. lgmres_inner_m See :func:`scipy.sparse.linalg.lgmres`. lgmres_outer_k See :func:`scipy.sparse.linalg.lgmres`. least_squares_lsmr_damp See :func:`scipy.sparse.linalg.lsmr`. least_squares_lsmr_atol See :func:`scipy.sparse.linalg.lsmr`. least_squares_lsmr_btol See :func:`scipy.sparse.linalg.lsmr`. least_squares_lsmr_conlim See :func:`scipy.sparse.linalg.lsmr`. least_squares_lsmr_maxiter See :func:`scipy.sparse.linalg.lsmr`. least_squares_lsmr_show See :func:`scipy.sparse.linalg.lsmr`. least_squares_lsqr_damp See :func:`scipy.sparse.linalg.lsqr`. least_squares_lsqr_atol See :func:`scipy.sparse.linalg.lsqr`. least_squares_lsqr_btol See :func:`scipy.sparse.linalg.lsqr`. least_squares_lsqr_conlim See :func:`scipy.sparse.linalg.lsqr`. least_squares_lsqr_iter_lim See :func:`scipy.sparse.linalg.lsqr`. least_squares_lsqr_show See :func:`scipy.sparse.linalg.lsqr`. Returns ------- A dict of available solvers with default |solver_options|. """ opts = {'scipy_bicgstab_spilu': {'type': 'scipy_bicgstab_spilu', 'tol': bicgstab_tol, 'maxiter': bicgstab_maxiter, 'spilu_drop_tol': spilu_drop_tol, 'spilu_fill_factor': spilu_fill_factor, 'spilu_drop_rule': spilu_drop_rule, 'spilu_permc_spec': spilu_permc_spec}, 'scipy_bicgstab': {'type': 'scipy_bicgstab', 'tol': bicgstab_tol, 'maxiter': bicgstab_maxiter}, 'scipy_spsolve': {'type': 'scipy_spsolve', 'permc_spec': spsolve_permc_spec, 'keep_factorization': spsolve_keep_factorization}, 'scipy_lgmres': {'type': 'scipy_lgmres', 'tol': lgmres_tol, 'maxiter': lgmres_maxiter, 'inner_m': lgmres_inner_m, 'outer_k': lgmres_outer_k}, 'scipy_least_squares_lsqr': {'type': 'scipy_least_squares_lsqr', 'damp': least_squares_lsqr_damp, 'atol': least_squares_lsqr_atol, 'btol': least_squares_lsqr_btol, 'conlim': least_squares_lsqr_conlim, 'iter_lim': least_squares_lsqr_iter_lim, 'show': least_squares_lsqr_show}} if config.HAVE_SCIPY_LSMR: opts['scipy_least_squares_lsmr'] = {'type': 'scipy_least_squares_lsmr', 'damp': least_squares_lsmr_damp, 'atol': least_squares_lsmr_atol, 'btol': least_squares_lsmr_btol, 'conlim': least_squares_lsmr_conlim, 'maxiter': least_squares_lsmr_maxiter, 'show': least_squares_lsmr_show} return opts
[docs]@defaults('check_finite', 'default_solver', 'default_least_squares_solver') def apply_inverse(op, V, options=None, least_squares=False, check_finite=True, default_solver='scipy_spsolve', default_least_squares_solver='scipy_least_squares_lsmr'): """Solve linear equation system. Applies the inverse of `op` to the vectors in `rhs` using SciPy. 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 :func:`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 (scipy_spsolve, scipy_bicgstab, scipy_bicgstab_spilu, scipy_lgmres, scipy_least_squares_lsmr, scipy_least_squares_lsqr). default_least_squares_solver Default solver to use for least squares problems (scipy_least_squares_lsmr, scipy_least_squares_lsqr). Returns ------- |VectorArray| of the solution vectors. """ assert V in op.range if isinstance(op, NumpyMatrixOperator): matrix = op.matrix else: from pymor.algorithms.to_matrix import to_matrix matrix = to_matrix(op) options = _parse_options(options, solver_options(), default_solver, default_least_squares_solver, least_squares) V = V.to_numpy() promoted_type = np.promote_types(matrix.dtype, V.dtype) R = np.empty((len(V), matrix.shape[1]), dtype=promoted_type) if options['type'] == 'scipy_bicgstab': for i, VV in enumerate(V): R[i], info = bicgstab(matrix, VV, tol=options['tol'], maxiter=options['maxiter']) if info != 0: if info > 0: raise InversionError(f'bicgstab failed to converge after {info} iterations') else: raise InversionError('bicgstab failed with error code {} (illegal input or breakdown)'. format(info)) elif options['type'] == 'scipy_bicgstab_spilu': if Version(scipy.version.version) >= Version('0.19'): ilu = spilu(matrix, drop_tol=options['spilu_drop_tol'], fill_factor=options['spilu_fill_factor'], drop_rule=options['spilu_drop_rule'], permc_spec=options['spilu_permc_spec']) else: if options['spilu_drop_rule']: logger = getLogger('pymor.operators.numpy._apply_inverse') logger.error("ignoring drop_rule in ilu factorization due to old SciPy") ilu = spilu(matrix, drop_tol=options['spilu_drop_tol'], fill_factor=options['spilu_fill_factor'], permc_spec=options['spilu_permc_spec']) precond = LinearOperator(matrix.shape, ilu.solve) for i, VV in enumerate(V): R[i], info = bicgstab(matrix, VV, tol=options['tol'], maxiter=options['maxiter'], M=precond) if info != 0: if info > 0: raise InversionError(f'bicgstab failed to converge after {info} iterations') else: raise InversionError('bicgstab failed with error code {} (illegal input or breakdown)'. format(info)) elif options['type'] == 'scipy_spsolve': try: # maybe remove unusable factorization: if hasattr(matrix, 'factorization'): fdtype = matrix.factorizationdtype if not np.can_cast(V.dtype, fdtype, casting='safe'): del matrix.factorization if Version(scipy.version.version) >= Version('0.14'): if hasattr(matrix, 'factorization'): # we may use a complex factorization of a real matrix to # apply it to a real vector. In that case, we downcast # the result here, removing the imaginary part, # which should be zero. R = matrix.factorization.solve(V.T).T.astype(promoted_type, copy=False) elif options['keep_factorization']: # the matrix is always converted to the promoted type. # if matrix.dtype == promoted_type, this is a no_op matrix.factorization = splu(matrix_astype_nocopy(matrix.tocsc(), promoted_type), permc_spec=options['permc_spec']) matrix.factorizationdtype = promoted_type R = matrix.factorization.solve(V.T).T else: # the matrix is always converted to the promoted type. # if matrix.dtype == promoted_type, this is a no_op R = spsolve(matrix_astype_nocopy(matrix, promoted_type), V.T, permc_spec=options['permc_spec']).T else: # see if-part for documentation if hasattr(matrix, 'factorization'): for i, VV in enumerate(V): R[i] = matrix.factorization.solve(VV).astype(promoted_type, copy=False) elif options['keep_factorization']: matrix.factorization = splu(matrix_astype_nocopy(matrix.tocsc(), promoted_type), permc_spec=options['permc_spec']) matrix.factorizationdtype = promoted_type for i, VV in enumerate(V): R[i] = matrix.factorization.solve(VV) elif len(V) > 1: factorization = splu(matrix_astype_nocopy(matrix.tocsc(), promoted_type), permc_spec=options['permc_spec']) for i, VV in enumerate(V): R[i] = factorization.solve(VV) else: R = spsolve(matrix_astype_nocopy(matrix, promoted_type), V.T, permc_spec=options['permc_spec']).reshape((1, -1)) except RuntimeError as e: raise InversionError(e) elif options['type'] == 'scipy_lgmres': for i, VV in enumerate(V): R[i], info = lgmres(matrix, VV, tol=options['tol'], maxiter=options['maxiter'], inner_m=options['inner_m'], outer_k=options['outer_k']) if info > 0: raise InversionError(f'lgmres failed to converge after {info} iterations') assert info == 0 elif options['type'] == 'scipy_least_squares_lsmr': from scipy.sparse.linalg import lsmr for i, VV in enumerate(V): R[i], info, itn, _, _, _, _, _ = lsmr(matrix, VV, damp=options['damp'], atol=options['atol'], btol=options['btol'], conlim=options['conlim'], maxiter=options['maxiter'], show=options['show']) assert 0 <= info <= 7 if info == 7: raise InversionError(f'lsmr failed to converge after {itn} iterations') elif options['type'] == 'scipy_least_squares_lsqr': for i, VV in enumerate(V): R[i], info, itn, _, _, _, _, _, _, _ = lsqr(matrix, VV, damp=options['damp'], atol=options['atol'], btol=options['btol'], conlim=options['conlim'], iter_lim=options['iter_lim'], show=options['show']) assert 0 <= info <= 7 if info == 7: raise InversionError(f'lsmr failed to converge after {itn} iterations') else: raise ValueError('Unknown solver type') if check_finite: if not np.isfinite(np.sum(R)): raise InversionError('Result contains non-finite values') return op.source.from_numpy(R)
# unfortunately, this is necessary, as scipy does not # forward the copy=False argument in its csc_matrix.astype function
[docs]def matrix_astype_nocopy(matrix, dtype): if matrix.dtype == dtype: return matrix else: return matrix.astype(dtype)
[docs]def lyap_lrcf_solver_options(): """Return available Lyapunov solvers with default options for the SciPy backend. Returns ------- A dict of available solvers with default solver options. """ return {'scipy': {'type': 'scipy'}}
[docs]def solve_lyap_lrcf(A, E, B, trans=False, options=None): """Compute an approximate low-rank solution of a Lyapunov equation. See :func:`pymor.algorithms.lyapunov.solve_lyap_lrcf` for a general description. This function uses `scipy.linalg.solve_continuous_lyapunov`, which is a dense solver for Lyapunov equations with E=I. Therefore, we assume A and E can be converted to |NumPy arrays| using :func:`~pymor.algorithms.to_matrix.to_matrix` and that `B.to_numpy` is implemented. .. note:: If E is not `None`, the problem will be reduced to a standard continuous-time algebraic Lyapunov equation by inverting E. 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 :func:`lyap_lrcf_solver_options`). Returns ------- Z Low-rank Cholesky factor of the Lyapunov equation solution, |VectorArray| from `A.source`. """ _solve_lyap_lrcf_check_args(A, E, B, trans) options = _parse_options(options, lyap_lrcf_solver_options(), 'scipy', None, False) X = solve_lyap_dense(to_matrix(A, format='dense'), to_matrix(E, format='dense') if E else None, B.to_numpy().T if not trans else B.to_numpy(), trans=trans, options=options) return A.source.from_numpy(_chol(X).T)
[docs]def lyap_dense_solver_options(): """Return available dense Lyapunov solvers with default options for the SciPy backend. Returns ------- A dict of available solvers with default solver options. """ return {'scipy': {'type': 'scipy'}}
[docs]def solve_lyap_dense(A, E, B, trans=False, options=None): """Compute the solution of a Lyapunov equation. See :func:`pymor.algorithms.lyapunov.solve_lyap_dense` for a general description. This function uses `scipy.linalg.solve_continuous_lyapunov`, which is a dense solver for Lyapunov equations with E=I. .. note:: If E is not `None`, the problem will be reduced to a standard continuous-time algebraic Lyapunov equation by inverting E. 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 :func:`lyap_dense_solver_options`). Returns ------- X Lyapunov equation solution as a |NumPy array|. """ _solve_lyap_dense_check_args(A, E, B, trans) options = _parse_options(options, lyap_dense_solver_options(), 'scipy', None, False) if options['type'] == 'scipy': if E is not None: A = solve(E, A) if not trans else solve(E.T, A.T).T B = solve(E, B) if not trans else solve(E.T, B.T).T if trans: A = A.T B = B.T X = solve_continuous_lyapunov(A, -B.dot(B.T)) else: raise ValueError(f"Unexpected Lyapunov equation solver ({options['type']}).") return X
[docs]def ricc_lrcf_solver_options(): """Return available Riccati solvers with default options for the SciPy backend. Returns ------- A dict of available solvers with default solver options. """ return {'scipy': {'type': 'scipy'}}
[docs]def solve_ricc_lrcf(A, E, B, C, R=None, S=None, trans=False, options=None): """Compute an approximate low-rank solution of a Riccati equation. See :func:`pymor.algorithms.riccati.solve_ricc_lrcf` for a general description. This function uses `scipy.linalg.solve_continuous_are`, which is a dense solver. Therefore, we assume all |Operators| and |VectorArrays| can be converted to |NumPy arrays| using :func:`~pymor.algorithms.to_matrix.to_matrix` and :func:`~pymor.vectorarrays.interfaces.VectorArrayInterface.to_numpy`. 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 :func:`ricc_lrcf_solver_options`). Returns ------- Z Low-rank Cholesky factor of the Riccati equation solution, |VectorArray| from `A.source`. """ _solve_ricc_check_args(A, E, B, C, R, S, trans) options = _parse_options(options, ricc_lrcf_solver_options(), 'scipy', None, False) if options['type'] != 'scipy': raise ValueError(f"Unexpected Riccati equation solver ({options['type']}).") A_source = A.source A = to_matrix(A, format='dense') E = to_matrix(E, format='dense') if E else None B = B.to_numpy().T C = C.to_numpy() S = S.to_numpy().T if S else None if R is None: R = np.eye(C.shape[0] if not trans else B.shape[1]) if not trans: if E is not None: E = E.T X = solve_continuous_are(A.T, C.T, B.dot(B.T), R, E, S) else: X = solve_continuous_are(A, B, C.T.dot(C), R, E, S) return A_source.from_numpy(_chol(X).T)
[docs]def pos_ricc_lrcf_solver_options(): """Return available positive Riccati solvers with default options for the SciPy backend. Returns ------- A dict of available solvers with default solver options. """ return {'scipy': {'type': 'scipy'}}
[docs]def solve_pos_ricc_lrcf(A, E, B, C, R=None, S=None, trans=False, options=None): """Compute an approximate low-rank solution of a positive Riccati equation. See :func:`pymor.algorithms.riccati.solve_pos_ricc_lrcf` for a general description. This function uses `scipy.linalg.solve_continuous_are`, which is a dense solver. Therefore, we assume all |Operators| and |VectorArrays| can be converted to |NumPy arrays| using :func:`~pymor.algorithms.to_matrix.to_matrix` and :func:`~pymor.vectorarrays.interfaces.VectorArrayInterface.to_numpy`. 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 :func:`pos_ricc_lrcf_solver_options`). Returns ------- Z Low-rank Cholesky factor of the positive Riccati equation solution, |VectorArray| from `A.source`. """ _solve_ricc_check_args(A, E, B, C, R, S, trans) options = _parse_options(options, pos_ricc_lrcf_solver_options(), 'scipy', None, False) if options['type'] != 'scipy': raise ValueError(f"Unexpected positive Riccati equation solver ({options['type']}).") if R is None: R = np.eye(len(C) if not trans else len(B)) return solve_ricc_lrcf(A, E, B, C, -R, S, trans, options)