docs: def bound_pseudo(arnoldifyer, Wt,
                   g_norm=0.,
                   G_norm=0.,
                   GW_norm=0.,
                   WGW_norm=0.,
                   tol=1e-6,
                   pseudo_type='auto',
                   pseudo_kwargs=None,
                   delta_n=20,
                   terminate_factor=1.
                   ):
      r'''Bound residual norms of next deflated system.
  
      :param arnoldifyer: an instance of
        :py:class:`~krypy.deflation.Arnoldifyer`.
      :param Wt: coefficients :math:`\tilde{W}\in\mathbb{C}^{n+d,k}` of the
        considered deflation vectors :math:`W` for the basis :math:`[V,U]`
        where ``V=last_solver.V`` and ``U=last_P.U``, i.e.,
        :math:`W=[V,U]\tilde{W}` and
        :math:`\mathcal{W}=\operatorname{colspan}(W)`. Must fulfill
        :math:`\tilde{W}^*\tilde{W}=I_k`.
      :param g_norm: norm :math:`\|g\|` of difference :math:`g=c-b` of
        right hand sides. Has to fulfill :math:`\|g\|<\|b\|`.
      :param G_norm: norm :math:`\|G\|` of difference
        :math:`G=B-A` of operators.
      :param GW_norm: Norm :math:`\|G|_{\mathcal{W}}\|` of difference
        :math:`G=B-A` of operators restricted to :math:`\mathcal{W}`.
      :param WGW_norm: Norm :math:`\|\langle W,GW\rangle\|_2`.
      :param pseudo_type: One of
  
        * ``'auto'``: determines if :math:`\hat{H}` is non-normal, normal or
          Hermitian and uses the corresponding mode (see other options below).
        * ``'nonnormal'``: the pseudospectrum of the Hessenberg matrix
          :math:`\hat{H}` is used (involves one computation of a pseudospectrum)
        * ``'normal'``: the pseudospectrum of :math:`\hat{H}` is computed
          efficiently by the union of circles around the eigenvalues.
        * ``'hermitian'``: the pseudospectrum of :math:`\hat{H}` is computed
          efficiently by the union of intervals around the eigenvalues.
        * ``'contain'``: the pseudospectrum of the extended Hessenberg matrix
          :math:`\begin{bmatrix}\hat{H}\\S_i\end{bmatrix}` is used
          (pseudospectrum has to be re computed for each iteration).
        * ``'omit'``: do not compute the pseudospectrum at all and just use the
          residual bounds from the approximate Krylov subspace.
      :param pseudo_kwargs: (optional) arguments that are passed to the method
        that computes the pseudospectrum.
      :param terminate_factor: (optional) terminate the computation if the ratio
        of two subsequent residual norms is larger than the provided factor.
        Defaults to 1.
      '''
      if pseudo_kwargs is None:
          pseudo_kwargs = {}
  
      # Arnoldify!
      Hh, Rh, q_norm, vdiff_norm, PWAW_norm = arnoldifyer.get(Wt)
  
      # get original linear system
      ls_orig = arnoldifyer._deflated_solver.linear_system
  
      k = Wt.shape[1]
      if k > 0:
          # smallest singular value of W^*AW
          WAW = Wt.T.conj().dot(arnoldifyer.J.dot(arnoldifyer.L.dot(Wt)))
          sigma_min = numpy.min(scipy.linalg.svdvals(WAW))
  
          if sigma_min <= WGW_norm:
              raise utils.AssumptionError(
                  'sigma_min(W^*AW) > ||W^*GW|| not satisfied.')
          eta = GW_norm/(sigma_min - WGW_norm)
      else:
          eta = 0.
      b_norm = ls_orig.MMlb_norm
      beta = PWAW_norm*(eta*(b_norm + g_norm) + g_norm) + vdiff_norm
  
      # check assumption on g_norm and b_norm
      if g_norm >= b_norm:
          raise utils.AssumptionError(
              '||g_norm|| < ||b_norm|| not satisfied')
  
      # compute residual norms of Hh*z=e_1*b_norm
      ls_small = linsys.LinearSystem(Hh,
                                     numpy.eye(Hh.shape[0], 1) * q_norm,
                                     normal=ls_orig.normal,
                                     self_adjoint=ls_orig.self_adjoint,
                                     positive_definite=ls_orig.positive_definite
                                     )
  
      Solver = type(arnoldifyer._deflated_solver)
      if issubclass(Solver, linsys.Minres) or issubclass(Solver, linsys.Gmres):
          aresnorms = utils.get_residual_norms(Hh,
                                               self_adjoint=ls_orig.self_adjoint)
      else:
          # TODO: compute residuals more efficiently for CG
          try:
              solver = Solver(ls_small, tol=tol, maxiter=Hh.shape[0])
          except utils.ConvergenceError as e:
              # use all residuals that have been computed
              # (useful for short recurrences)
              solver = e.solver
          aresnorms = numpy.array(solver.resnorms)
      # absolute residual norm
      aresnorms = aresnorms * q_norm
  
      if pseudo_type == 'omit':
          return aresnorms / (b_norm - g_norm)
  
      # spectrum of Hh
      evals, evecs = scipy.linalg.eig(Hh)
      if ls_small.self_adjoint:
          evals = numpy.real(evals)
  
      # norm of Hh
      Hh_norm = numpy.linalg.norm(Hh, 2)
  
      def _auto():
          '''determine pseudo automatically'''
          # is Hh Hermitian?
          if numpy.linalg.norm(Hh-Hh.T.conj(), 2) < 1e-14*Hh_norm:
              return 'hermitian'
  
          # is Hh normal?
          if numpy.linalg.cond(evecs, 2) < 1+1e-14:
              return 'normal'
  
          return 'nonnormal'
  
      if pseudo_type == 'auto':
          pseudo_type = _auto()
  
      # for delta >= min(|\lambda|), the pseudospectrum will contain zero and
      # the thus polymax > 1. nevertheless, the bound may provide useful
      # information in early iterations with large values of delta.
      # Therefore, the maximal perturbation is chosen as the maximal
      # eigenvalue of Hh
      delta_max = 1e2*numpy.max(numpy.abs(evals))
  
      # minimal delta is defined via Rh
      # HACK until numpy.linal.svd (and thus numpy.linalg.norm) is fixed
      from scipy.linalg import svd
      _, Rhsvd, _ = svd(Rh[:, :1])
      delta_min = PWAW_norm*(eta*(Hh_norm + G_norm) + G_norm) + numpy.max(Rhsvd)
      if delta_min == 0:
          delta_min = 1e-16
  
      import pseudopy
      if not ls_small.normal:
          # construct pseudospectrum for the expected range
          pseudo = pseudopy.NonnormalAuto(Hh, delta_min*0.99, delta_max*1.01,
                                          **pseudo_kwargs)
      elif not ls_small.self_adjoint:
          pseudo = pseudopy.NormalEvals(evals)
      else:
          pseudo = None
  
      bounds = [aresnorms[0]]
      for i in range(1, len(aresnorms)):
          # compute roots of polynomial
          if issubclass(Solver, linsys.Cg):
              roots = scipy.linalg.eigvalsh(Hh[:i, :i])
          else:
              # TODO: more stable way of computing the roots of the MINRES
              #       poly with exploitation of symmetry?
              HhQ, HhR = scipy.linalg.qr(Hh[:i+1, :i], mode='economic')
              roots_inv = scipy.linalg.eigvals(HhQ[:i, :].T.conj(), HhR)
              roots = 1./roots_inv[numpy.abs(roots_inv) > 1e-14]
  
          if ls_small.self_adjoint:
              roots = numpy.real(roots)
  
          # compute polynomial
          p = utils.NormalizedRootsPolynomial(roots)
          if ls_small.self_adjoint:
              p_minmax_candidates = p.minmax_candidates()
  
          # absolute residual
          aresnorm = aresnorms[i]
  
          # perturbation
          # HACK until numpy.linal.svd (and thus numpy.linalg.norm) is fixed
          from scipy.linalg import svd
          _, Rhsvd, _ = svd(Rh[:, :i])
          Rhnrm = numpy.max(Rhsvd)
          epsilon = PWAW_norm*(eta*(Hh_norm + G_norm) + G_norm) \
              + Rhnrm
              #+ numpy.linalg.norm(Rh[:, :i], 2)
          if epsilon == 0:
              epsilon = 1e-16
  
          if pseudo_type == 'contain':
              raise NotImplementedError('contain not yet implemented')
  
          # exit if epsilon >= delta_max
          if epsilon >= delta_max:
              break
  
          delta_log_range = numpy.linspace(numpy.log10(1.01*epsilon),
                                           numpy.log10(delta_max),
                                           delta_n+2
                                           )[0:-1]
  
          def compute_pseudo(delta_log):
              delta = 10**delta_log
              if ls_small.self_adjoint:
                  # pseudospectrum are intervals
                  pseudo_intervals = utils.Intervals(
                      [utils.Interval(ev-delta, ev+delta) for ev in evals])
  
                  # add roots of first derivative of p
                  candidates = []
                  for candidate in p_minmax_candidates:
                      if pseudo_intervals.contains(candidate):
                          candidates.append(candidate)
                  all_candidates = numpy.r_[pseudo_intervals.get_endpoints(),
                                            numpy.array(candidates)]
  
                  # evaluate polynomial
                  polymax = numpy.max(numpy.abs(p(all_candidates)))
                  pseudolen = 2 * delta
  
              else:
                  # get pseudospectrum paths
                  pseudo_path = pseudo.contour_paths(delta)
  
                  # length of boundary
                  pseudolen = pseudo_path.length()
  
                  # evaluate polynomial on points of path
                  if pseudolen > 0:
                      polymax = numpy.max(numpy.abs(p(pseudo_path.vertices())))
                  else:
                      polymax = numpy.Inf
  
              # compute THE bound
              return pseudolen/(2*numpy.pi*delta) \
                  * (epsilon/(delta-epsilon)*(q_norm + beta) + beta) \
                  * polymax
  
          # minimization
          from scipy.optimize import minimize_scalar
          opt_res = minimize_scalar(compute_pseudo,
                                    bounds=(delta_log_range[0],
                                            delta_log_range[-1]),
                                    method='bounded',
                                    options={'maxiter': delta_n}
                                    )
          # the delta with minimal value is min_delta = 10**opt_res.x
          min_val = opt_res.fun
  
          # minimal bound value
          boundval = aresnorm + min_val
  
          # if not increasing: append to bounds
          if i > 1 and boundval/bounds[-1] > terminate_factor:
              break
          else:
              bounds.append(numpy.min([boundval, bounds[-1]]))
      return numpy.array(bounds) / (b_norm - g_norm)