[iosys] minor improvements in reductors

- removed experimental algorithms.arnoldi.arnoldi_tangential
- in irka and tf_irka, force_sigma_in_rhp is defaulted to False, as in the
  referenced papers
- changed the convergence criterion names in irka and tf_irka
- added more comments

src/pymor/algorithms/arnoldi.py

@@ -3,7 +3,6 @@
 # License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
 
 from pymor.operators.constructions import LincombOperator, VectorArrayOperator
-from pymor.operators.numpy import NumpyMatrixOperator
 
 
 def arnoldi(A, E, b, sigma, trans=False):
@@ -111,155 +110,3 @@ def arnoldi(A, E, b, sigma, trans=False):
             v = v2
 
     return V
-
-
-def arnoldi_tangential(A, E, B, sigma, directions, trans=False):
-    r"""Tangential Rational Arnoldi algorithm.
-
-    If `trans == False`, using tangential Arnoldi process [DSZ14]_,
-    computes a real orthonormal basis for the rational Krylov subspace
-
-    .. math::
-        \span\{(\sigma_1 E - A)^{-1} B d_1, (\sigma_2 E - A)^{-1} B d_2,
-        \ldots, (\sigma_r E - A)^{-1} B d_r\},
-
-    otherwise, computes the same for
-
-    .. math::
-        \span\{(\sigma_1 E - A)^{-*} B^* d_1, (\sigma_2 E - A)^{-*} B^* d_2,
-        \ldots, (\sigma_r E - A)^{-*} B^* d_r\}.
-
-    Interpolation points in `sigma` are allowed to repeat (in any order).
-    Then, in the above expression,
-
-    .. math::
-        \underbrace{(\sigma_i E - A)^{-1} B^* d_i, \ldots,
-        (\sigma_i E - A)^{-1} B^* d_i}_{m \text{ times}}
-
-    is replaced by
-
-    .. math::
-        (\sigma_i E - A)^{-1} B^* d_i, (\sigma_i E - A)^{-2} B^* d_i,
-        \ldots, (\sigma_i E - A)^{-m} B^* d_i.
-
-    Analogously for the `trans == True` case.
-
-    .. warning::
-        The implementation is still experimental.
-
-    .. [DSZ14] V. Druskin, V. Simoncini, M. Zaslavsky, Adaptive Tangential
-               Interpolation in Rational Krylov Subspaces for MIMO Dynamical
-               Systems,
-               SIAM Journal on Matrix Analysis and Applications, 35(2),
-               476-498, 2014.
-
-    Parameters
-    ----------
-    A
-        Real |Operator| A.
-    E
-        Real |Operator| E.
-    B
-        Real |Operator| B.
-    sigma
-        Interpolation points (closed under conjugation).
-    directions
-        Tangential directions (closed under conjugation), |VectorArray| of
-        length `len(sigma)` from `B.source` or `B.range` (depending on
-        `trans`).
-    trans
-        Boolean, see above.
-
-    Returns
-    -------
-    V
-        Projection matrix.
-    """
-    r = len(sigma)
-    assert len(directions) == r
-    assert (not trans and B.source.dim > 1 and directions in B.source or
-            trans and B.range.dim > 1 and directions in B.range)
-
-    directions.scal(1 / directions.l2_norm())
-
-    V = A.source.empty(reserve=r)
-
-    for i in range(r):
-        if sigma[i].imag == 0:
-            sEmA = LincombOperator((E, A), (sigma[i].real, -1))
-
-            if not trans:
-                if i == 0:
-                    v = sEmA.apply_inverse(B.apply(directions.real[0]))
-                else:
-                    Bd = B.apply(directions.real[i])
-                    VTBd = VectorArrayOperator(V, transposed=True).apply(Bd)
-                    sEmA_proj_inv_VTBd = NumpyMatrixOperator(sEmA.apply2(V, V)).apply_inverse(VTBd)
-                    V_sEmA_proj_inv_VTBd = VectorArrayOperator(V).apply(sEmA_proj_inv_VTBd)
-                    rd = sEmA.apply(V_sEmA_proj_inv_VTBd) - Bd
-                    v = sEmA.apply_inverse(rd)
-            else:
-                if i == 0:
-                    v = sEmA.apply_inverse_adjoint(B.apply_adjoint(directions.real[0]))
-                else:
-                    CTd = B.apply_adjoint(directions.real[i])
-                    VTCTd = VectorArrayOperator(V, transposed=True).apply(CTd)
-                    sEmA_proj_inv_VTCTd = NumpyMatrixOperator(sEmA.apply2(V, V)).apply_inverse_adjoint(VTCTd)
-                    V_sEmA_proj_inv_VTCTd = VectorArrayOperator(V).apply(sEmA_proj_inv_VTCTd)
-                    rd = sEmA.apply_adjoint(V_sEmA_proj_inv_VTCTd) - CTd
-                    v = sEmA.apply_inverse_adjoint(rd)
-
-            if i > 0:
-                v_norm_orig = v.l2_norm()[0]
-                Vop = VectorArrayOperator(V)
-                v -= Vop.apply(Vop.apply_adjoint(v))
-                if v.l2_norm()[0] < v_norm_orig / 10:
-                    v -= Vop.apply(Vop.apply_adjoint(v))
-            v.scal(1 / v.l2_norm()[0])
-            V.append(v)
-        elif sigma[i].imag > 0:
-            sEmA = LincombOperator((E, A), (sigma[i], -1))
-
-            if not trans:
-                if i == 0:
-                    v = sEmA.apply_inverse(B.apply(directions[0]))
-                else:
-                    Bd = B.apply(directions[i])
-                    VTBd = VectorArrayOperator(V, transposed=True).apply(Bd)
-                    sEmA_proj_inv_VTBd = NumpyMatrixOperator(sEmA.apply2(V, V)).apply_inverse(VTBd)
-                    V_sEmA_proj_inv_VTBd = VectorArrayOperator(V).apply(sEmA_proj_inv_VTBd)
-                    rd = sEmA.apply(V_sEmA_proj_inv_VTBd) - Bd
-                    v = sEmA.apply_inverse(rd)
-            else:
-                if i == 0:
-                    v = sEmA.apply_inverse_adjoint(B.apply_adjoint(directions[0]))
-                else:
-                    CTd = B.apply_adjoint(directions[i])
-                    VTCTd = VectorArrayOperator(V, transposed=True).apply(CTd)
-                    sEmA_proj_inv_VTCTd = NumpyMatrixOperator(sEmA.apply2(V, V)).apply_inverse_adjoint(VTCTd)
-                    V_sEmA_proj_inv_VTCTd = VectorArrayOperator(V).apply(sEmA_proj_inv_VTCTd)
-                    rd = sEmA.apply_adjoint(V_sEmA_proj_inv_VTCTd) - CTd
-                    v = sEmA.apply_inverse_adjoint(rd)
-
-            v1 = v.real
-            if i > 0:
-                v1_norm_orig = v1.l2_norm()[0]
-                Vop = VectorArrayOperator(V)
-                v1 -= Vop.apply(Vop.apply_adjoint(v1))
-                if v1.l2_norm()[0] < v1_norm_orig / 10:
-                    v1 -= Vop.apply(Vop.apply_adjoint(v1))
-            v1.scal(1 / v1.l2_norm()[0])
-            V.append(v1)
-
-            v2 = v.imag
-            v2_norm_orig = v2.l2_norm()[0]
-            Vop = VectorArrayOperator(V)
-            v2 -= Vop.apply(Vop.apply_adjoint(v2))
-            if v2.l2_norm()[0] < v2_norm_orig / 10:
-                v2 -= Vop.apply(Vop.apply_adjoint(v2))
-            v2.scal(1 / v2.l2_norm()[0])
-            V.append(v2)
-
-            v = v2
-
-    return V

src/pymor/algorithms/cholp.py

@@ -11,7 +11,6 @@ def cholp(A, copy=True):
     """Low-rank approximation using pivoted Cholesky decomposition
 
     .. note::
-
         Should be replaced with LAPACK routine DPSTRF (when it becomes available in NumPy).
 
     .. [H02] N. J. Higham, Accuracy and Stability of Numerical Algorithms,

src/pymor/reductors/bt.py

@@ -61,14 +61,18 @@ def bt(discretization, r=None, tol=None, typ='lyap', me_solver=None, method='bfs
     assert r is None or 0 < r < discretization.n
     assert method in ('sr', 'bfsr')
 
-    sv, U, V = discretization.sv_U_V(typ, me_solver=me_solver)
+    # compute gramian factors
     cf = discretization.gramian(typ, 'cf', me_solver=me_solver)
     of = discretization.gramian(typ, 'of', me_solver=me_solver)
 
     if r is not None and r > min([len(cf), len(of)]):
         raise ValueError('r needs to be smaller than the sizes of Gramian factors.' +
                          ' Try reducing the tolerance in the low-rank Lyapunov equation solver.')
 
+    # compute "Hankel" singular values and vectors
+    sv, U, V = discretization.sv_U_V(typ, me_solver=me_solver)
+
+    # find reduced order if tol is specified
     if tol is not None:
         bounds = np.zeros((discretization.n,))
         sv_reverse = np.zeros((discretization.n,))
@@ -84,9 +88,9 @@ def bt(discretization, r=None, tol=None, typ='lyap', me_solver=None, method='bfs
         r_tol = np.argmax(bounds <= tol) + 1
         r = r_tol if r is None else min([r, r_tol])
 
+    # compute projection matrices and find the reduced model
     V = cf.lincomb(V[:r])
     W = of.lincomb(U[:r])
-
     if method == 'sr':
         alpha = 1 / np.sqrt(sv[:r])
         V.scal(alpha)

src/pymor/reductors/lti.py

@@ -5,9 +5,10 @@
 import numpy as np
 import scipy.linalg as spla
 
-from pymor.algorithms.arnoldi import arnoldi, arnoldi_tangential
+from pymor.algorithms.arnoldi import arnoldi
 from pymor.algorithms.gram_schmidt import gram_schmidt
-from pymor.operators.constructions import LincombOperator
+from pymor.algorithms.to_matrix import to_matrix
+from pymor.operators.constructions import IdentityOperator, LincombOperator
 from pymor.reductors.basic import reduce_generic_pg
 
 
@@ -28,6 +29,7 @@ def interpolation(discretization, sigma, b, c, use_arnoldi=False):
         `discretization.C.range`.
     use_arnoldi
         Should the Arnoldi process be used for rational interpolation.
+        Available only for SISO systems. Otherwise, it is ignored.
 
     Returns
     -------
@@ -43,22 +45,17 @@ def interpolation(discretization, sigma, b, c, use_arnoldi=False):
     assert b in discretization.B.source and len(b) == r
     assert c in discretization.C.range and len(c) == r
 
-    if use_arnoldi:
-        if discretization.m == 1:
-            V = arnoldi(discretization.A, discretization.E, discretization.B, sigma)
-        else:
-            V = arnoldi_tangential(discretization.A, discretization.E, discretization.B, sigma, b)
-        if discretization.p == 1:
-            W = arnoldi(discretization.A, discretization.E, discretization.C, sigma, trans=True)
-        else:
-            W = arnoldi_tangential(discretization.A, discretization.E, discretization.C, sigma, c, trans=True)
+    if use_arnoldi and discretization.m == 1 and discretization.p == 1:
+        V = arnoldi(discretization.A, discretization.E, discretization.B, sigma)
+        W = arnoldi(discretization.A, discretization.E, discretization.C, sigma, trans=True)
     else:
+        # rescale tangential directions (could avoid overflow or underflow)
         b.scal(1 / b.l2_norm())
         c.scal(1 / c.l2_norm())
 
+        # compute projection matrices
         V = discretization.A.source.empty(reserve=r)
         W = discretization.A.source.empty(reserve=r)
-
         for i in range(r):
             if sigma[i].imag == 0:
                 sEmA = LincombOperator((discretization.E, discretization.A), (sigma[i].real, -1))
@@ -90,8 +87,8 @@ def interpolation(discretization, sigma, b, c, use_arnoldi=False):
     return rom, rc, reduction_data
 
 
-def irka(discretization, r, sigma=None, b=None, c=None, tol=1e-4, maxit=100, verbose=False, force_sigma_in_rhp=True,
-         use_arnoldi=False, conv_crit='rel_sigma', compute_errors=False):
+def irka(discretization, r, sigma=None, b=None, c=None, tol=1e-4, maxit=100, verbose=False, force_sigma_in_rhp=False,
+         use_arnoldi=False, conv_crit='rel_sigma_change', compute_errors=False):
     r"""Reduce using IRKA.
 
     .. [GAB08] S. Gugercin, A. C. Antoulas, C. A. Beattie,
@@ -132,16 +129,17 @@ def irka(discretization, r, sigma=None, b=None, c=None, tol=1e-4, maxit=100, ver
     verbose
         Should consecutive distances be printed.
     force_sigma_in_rhp
-        If `True`, new interpolation points are always in the right
-        half-plane. Otherwise, they are reflections of reduced order model's
-        poles.
+        If 'False`, new interpolation are reflections of reduced order
+        model's poles. Otherwise, they are always in the right
+        half-plane.
     use_arnoldi
         Should the Arnoldi process be used for rational interpolation.
+        Available only for SISO systems. Otherwise, it is ignored.
     conv_crit
         Convergence criterion:
-            - `'rel_sigma'`: relative change in interpolation points
-            - `'max_sin_PG'`: maximum of sines in Petrov-Galerkin subspaces
-            - `'rel_H2'`: relative H_2 distance of reduced order models
+            - `'rel_sigma_change'`: relative change in interpolation points
+            - `'subspace_sin'`: maximum of sines of Petrov-Galerkin subspaces
+            - `'rel_H2_dist'`: relative H_2 distance of reduced order models
     compute_errors
         Should the relative :math:`\mathcal{H}_2`-errors of intermediate
         reduced order models be computed.
@@ -168,12 +166,16 @@ def irka(discretization, r, sigma=None, b=None, c=None, tol=1e-4, maxit=100, ver
         - relative :math:`\mathcal{H}_2`-errors `errors` (if
           `compute_errors` is `True`).
     """
+    if not discretization.cont_time:
+        raise NotImplementedError
     assert 0 < r < discretization.n
     assert sigma is None or len(sigma) == r
     assert b is None or b in discretization.B.source and len(b) == r
     assert c is None or c in discretization.C.range and len(c) == r
-    assert conv_crit in ('rel_sigma', 'max_sin_PG', 'rel_H2')
+    assert conv_crit in ('rel_sigma_change', 'subspace_sin', 'rel_H2_dist')
 
+    # basic choice for initial interpolation points and tangential
+    # directions
     if sigma is None:
         sigma = np.logspace(-1, 1, r)
     if b is None:
@@ -199,7 +201,9 @@ def irka(discretization, r, sigma=None, b=None, c=None, tol=1e-4, maxit=100, ver
     L = [c]
     if compute_errors:
         errors = []
+    # main loop
     for it in range(maxit):
+        # interpolatory reduced order model
         rom, rc, reduction_data = interpolation(discretization, sigma, b, c, use_arnoldi=use_arnoldi)
 
         if compute_errors:
@@ -210,16 +214,21 @@ def irka(discretization, r, sigma=None, b=None, c=None, tol=1e-4, maxit=100, ver
                 rel_H2_err = np.inf
             errors.append(rel_H2_err)
 
-        sigma, Y, X = spla.eig(rom.A._matrix, rom.E._matrix, left=True, right=True)
+        # new interpolation points
+        if isinstance(rom.E, IdentityOperator):
+            sigma, Y, X = spla.eig(to_matrix(rom.A), left=True, right=True)
+        else:
+            sigma, Y, X = spla.eig(to_matrix(rom.A), to_matrix(rom.E), left=True, right=True)
         if force_sigma_in_rhp:
             sigma = np.array([np.abs(s.real) + s.imag * 1j for s in sigma])
         else:
             sigma *= -1
         Sigma.append(sigma)
 
-        if conv_crit == 'rel_sigma':
+        # compute convergence criterion
+        if conv_crit == 'rel_sigma_change':
             dist.append(spla.norm((Sigma[-2] - Sigma[-1]) / Sigma[-2], ord=np.inf))
-        elif conv_crit == 'max_sin_PG':
+        elif conv_crit == 'subspace_sin':
             if it == 0:
                 V_new = reduction_data['V'].data.T
                 W_new = reduction_data['W'].data.T
@@ -232,7 +241,7 @@ def irka(discretization, r, sigma=None, b=None, c=None, tol=1e-4, maxit=100, ver
                 sinV = spla.norm(V_new - V_old.dot(V_old.T.dot(V_new)), ord=2)
                 sinW = spla.norm(W_new - W_old.dot(W_old.T.dot(W_new)), ord=2)
                 dist.append(np.max([sinV, sinW]))
-        elif conv_crit == 'rel_H2':
+        elif conv_crit == 'rel_H2_dist':
             if it == 0:
                 rom_new = rom
                 dist.append(np.inf)
@@ -252,16 +261,19 @@ def irka(discretization, r, sigma=None, b=None, c=None, tol=1e-4, maxit=100, ver
             else:
                 print('{:4d} | {:15.9e}'.format(it + 1, dist[-1]))
 
+        # new tangential directions
         Y = rom.B.range.make_array(Y.conj().T)
         X = rom.C.source.make_array(X.T)
         b = rom.B.apply_adjoint(Y)
         c = rom.C.apply(X)
         R.append(b)
         L.append(c)
 
+        # check if convergence criterion is satisfied
         if dist[-1] < tol:
             break
 
+    # final reduced order model
     rom, rc, reduction_data = interpolation(discretization, sigma, b, c, use_arnoldi=use_arnoldi)
 
     reduction_data.update({'dist': dist, 'Sigma': Sigma, 'R': R, 'L': L})