# Source code for pymor.algorithms.lyapunov

# This file is part of the pyMOR project (http://www.pymor.org).

import numpy as np
import scipy.linalg as spla

from pymor.core.config import config
from pymor.core.defaults import defaults
from pymor.operators.interface import Operator

_DEFAULT_LYAP_LRCF_SPARSE_SOLVER_BACKEND = ('pymess' if config.HAVE_PYMESS else

_DEFAULT_LYAP_LRCF_DENSE_SOLVER_BACKEND = ('pymess' if config.HAVE_PYMESS else
'slycot' if config.HAVE_SLYCOT else
'scipy')

_DEFAULT_LYAP_DENSE_SOLVER_BACKEND = ('pymess' if config.HAVE_PYMESS else
'slycot' if config.HAVE_SLYCOT else
'scipy')

[docs]@defaults('value')
def mat_eqn_sparse_min_size(value=1000):
"""Returns minimal size for which a sparse solver will be used by default."""
return value

[docs]@defaults('default_sparse_solver_backend', 'default_dense_solver_backend')
def solve_lyap_lrcf(A, E, B, trans=False, options=None,
default_sparse_solver_backend=_DEFAULT_LYAP_LRCF_SPARSE_SOLVER_BACKEND,
default_dense_solver_backend=_DEFAULT_LYAP_LRCF_DENSE_SOLVER_BACKEND):
"""Compute an approximate low-rank solution of a Lyapunov equation.

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

- if trans is False and E is None:

.. math::
A X + X A^T + B B^T = 0,

- if trans is False and E is an |Operator|:

.. math::
A X E^T + E X A^T + B B^T = 0,

- if trans is True and E is None:

.. math::
A^T X + X A + B^T B = 0,

- if trans is True and E is an |Operator|:

.. math::
A^T X E + E^T X A + B^T B = 0.

We assume A and E are real |Operators|, E is invertible, and all the
eigenvalues of (A, E) all lie in the open left half-plane.
Operator B needs to be given as a |VectorArray| from A.source, and
for large-scale problems, we assume len(B) is small.

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

- for sparse problems (minimum size specified by
:func:mat_eqn_sparse_min_size)

1. pymess (see :func:pymor.bindings.pymess.solve_lyap_lrcf),
2. lradi (see :func:pymor.algorithms.lradi.solve_lyap_lrcf),

- for dense problems (smaller than :func:mat_eqn_sparse_min_size)

1. pymess (see :func:pymor.bindings.pymess.solve_lyap_lrcf),
2. slycot (see :func:pymor.bindings.slycot.solve_lyap_lrcf),
3. scipy (see :func:pymor.bindings.scipy.solve_lyap_lrcf).

Parameters
----------
A
The non-parametric |Operator| A.
E
The non-parametric |Operator| E or None.
B
The operator B as a |VectorArray| from A.source.
trans
Whether the first |Operator| in the Lyapunov equation is
transposed.
options
The solver options to use.
See:

- :func:pymor.algorithms.lradi.lyap_lrcf_solver_options,
- :func:pymor.bindings.scipy.lyap_lrcf_solver_options,
- :func:pymor.bindings.slycot.lyap_lrcf_solver_options,
- :func:pymor.bindings.pymess.lyap_lrcf_solver_options.

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

_solve_lyap_lrcf_check_args(A, E, B, trans)
if options:
solver = options if isinstance(options, str) else options['type']
backend = solver.split('_')[0]
else:
if A.source.dim >= mat_eqn_sparse_min_size():
backend = default_sparse_solver_backend
else:
backend = default_dense_solver_backend
if backend == 'scipy':
from pymor.bindings.scipy import solve_lyap_lrcf as solve_lyap_impl
elif backend == 'slycot':
from pymor.bindings.slycot import solve_lyap_lrcf as solve_lyap_impl
elif backend == 'pymess':
from pymor.bindings.pymess import solve_lyap_lrcf as solve_lyap_impl
from pymor.algorithms.lradi import solve_lyap_lrcf as solve_lyap_impl
else:
raise ValueError(f'Unknown solver backend ({backend}).')
return solve_lyap_impl(A, E, B, trans=trans, options=options)

def _solve_lyap_lrcf_check_args(A, E, B, trans):
assert isinstance(A, Operator) and A.linear
assert not A.parametric
assert A.source == A.range
if E is not None:
assert isinstance(E, Operator) and E.linear
assert not E.parametric
assert E.source == E.range
assert E.source == A.source
assert B in A.source

[docs]@defaults('default_solver_backend')
def solve_lyap_dense(A, E, B, trans=False, options=None,
default_solver_backend=_DEFAULT_LYAP_DENSE_SOLVER_BACKEND):
"""Compute the solution of a Lyapunov equation.

Returns the solution :math:X of a (generalized) continuous-time
algebraic Lyapunov equation:

- if trans is False and E is None:

.. math::
A X + X A^T + B B^T = 0,

- if trans is False and E is an |Operator|:

.. math::
A X E^T + E X A^T + B B^T = 0,

- if trans is True and E is None:

.. math::
A^T X + X A + B^T B = 0,

- if trans is True and E is an |Operator|:

.. math::
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 :func:pymor.bindings.pymess.solve_lyap_dense)
2. slycot (see :func:pymor.bindings.slycot.solve_lyap_dense)
3. scipy (see :func: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:

- :func:pymor.bindings.scipy.lyap_dense_solver_options,
- :func:pymor.bindings.slycot.lyap_dense_solver_options,
- :func:pymor.bindings.pymess.lyap_dense_solver_options.

default_solver_backend
Default solver backend to use (pymess, slycot, scipy).

Returns
-------
X
Lyapunov equation solution as a |NumPy array|.
"""

_solve_lyap_dense_check_args(A, E, B, trans)
if options:
solver = options if isinstance(options, str) else options['type']
backend = solver.split('_')[0]
else:
backend = default_solver_backend
if backend == 'scipy':
from pymor.bindings.scipy import solve_lyap_dense as solve_lyap_impl
elif backend == 'slycot':
from pymor.bindings.slycot import solve_lyap_dense as solve_lyap_impl
elif backend == 'pymess':
from pymor.bindings.pymess import solve_lyap_dense as solve_lyap_impl
else:
raise ValueError(f'Unknown solver backend ({backend}).')
return solve_lyap_impl(A, E, B, trans, options=options)

def _solve_lyap_dense_check_args(A, E, B, trans):
assert isinstance(A, np.ndarray) and A.ndim == 2
assert A.shape[0] == A.shape[1]
if E is not None:
assert isinstance(E, np.ndarray) and E.ndim == 2
assert E.shape[0] == E.shape[1]
assert E.shape[0] == A.shape[0]
assert isinstance(B, np.ndarray) and A.ndim == 2
assert not trans and B.shape[0] == A.shape[0] or trans and B.shape[1] == A.shape[0]

[docs]def _chol(A):
"""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).
"""

assert isinstance(A, np.ndarray) and A.ndim == 2
assert A.shape[0] == A.shape[1]

U, s, _ = spla.svd(A, lapack_driver='gesvd')
L = U * np.sqrt(s)
return L