docs: def arnoldi_projected(H, P, k, ortho='mgs'):
      """Compute (perturbed) Arnoldi relation for projected operator.
  
      Assume that you have computed an Arnoldi relation
  
      .. math ::
  
          A V_n = V_{n+1} \\underline{H}_n
  
      where :math:`V_{n+1}\\in\\mathbb{C}^{N,n+1}` has orthogonal columns
      (with respect to an inner product :math:`\\langle\\cdot,\\cdot\\rangle`)
      and :math:`\\underline{H}_n\\in\\mathbb{C}^{n+1,n}` is an extended
      upper Hessenberg matrix.
  
      For :math:`k<n` you choose full rank matrices
      :math:`X\\in\\mathbb{C}^{n-1,k}` and :math:`Y\\in\\mathbb{C}^{n,k}` and
      define :math:`\\tilde{X}:=A V_{n_1}X = V_n \\underline{H}_{n-1} X` and
      :math:`\\tilde{Y}:=V_n Y` such that
      :math:`\\langle \\tilde{Y}, \\tilde{X} \\rangle = Y^*\\underline{H}_{n-1} X`
      is invertible. Then the projections :math:`P` and :math:`\\tilde{P}`
      characterized by
  
      * :math:`\\tilde{P}x = x -
        \\tilde{X} \\langle \\tilde{Y},\\tilde{X} \\rangle^{-1}
        \\langle\\tilde{Y},x\\rangle`
      * :math:`P = I - \\underline{H}_{n-1}X (Y^*\\underline{H}_{n-1}X)^{-1}Y^*`
  
      are well defined and :math:`\\tilde{P}V_{n+1} = [V_n P, v_{n+1}]` holds.
  
      This method computes for :math:`i<n-k` the Arnoldi relation
  
      .. math ::
  
          (\\tilde{P}A + E_i) W_i
          = W_{i+1} \\underline{G}_i
  
      where :math:`W_{i+1}=V_n U_{i+1}` has orthogonal columns with respect
      to :math:`\\langle\\cdot,\\cdot\\rangle`,
      :math:`\\underline{G}_i` is an extended upper Hessenberg matrix
      and :math:`E_i x = v_{n+1} F_i \\langle W_i,x\\rangle` with
      :math:`F_i=[f_1,\ldots,f_i]\\in\\mathbb{C}^{1,i}`.
  
      The perturbed Arnoldi relation can also be generated with the operator
      :math:`P_{V_n} \\tilde{P} A`:
  
      .. math ::
  
          P_{V_n} \\tilde{P} A W_i
          = W_{i+1} \\underline{G}_i.
  
      In a sense the perturbed Arnoldi relation is the best prediction for the
      behavior of the Krylov subspace :math:`K_i(\\tilde{P}A,\\tilde{P}v_1)`
      that can be generated only with the data from :math:`K_{n+1}(A,v_1)` and
      without carrying out further matrix-vector multiplications with A.
  
      :param H: the extended upper Hessenberg matrix
        :math:`\\underline{H}_n` with ``shape==(n+1,n)``.
      :param P: the projection
        :math:`P:\\mathbb{C}^n\\longrightarrow\\mathbb{C}^n` (has to be
        compatible with :py:meth:`get_linearoperator`).
      :param k: the dimension of the null space of P.
      :returns: U, G, F where
  
        * U is the coefficient matrix :math:`U_{i+1}` with ``shape==(n,i+1)``,
        * G is the extended upper Hessenberg matrix :math:`\\underline{G}_i`
          with ``shape==(i+1,i)``,
        * F is the error matrix :math:`F_i` with ``shape==(1,i)``.
      """
      n = H.shape[1]
      dtype = find_common_dtype(H, P)
      invariant = H.shape[0] == n
      hlast = 0 if invariant else H[-1, -1]
      H = get_linearoperator((n, n), H if invariant else H[:-1, :])
      P = get_linearoperator((n, n), P)
      v = P * numpy.eye(n, 1)
      maxiter = n-k+1
      F = numpy.zeros((1, maxiter), dtype=dtype)
      PH = lambda x: P*(H*x)
      PH = LinearOperator((n, n), dtype, PH)
      _arnoldi = Arnoldi(PH, v, maxiter=maxiter, ortho=ortho)
      while _arnoldi.iter < _arnoldi.maxiter and not _arnoldi.invariant:
          u, _ = _arnoldi.get_last()
          F[0, _arnoldi.iter] = hlast * u[-1, 0]
          _arnoldi.advance()
      U, G = _arnoldi.get()
      return U, G, F[[0], :_arnoldi.iter]