docs: def relaxation_decomp(P, p0, obs, times=[1], k=None, ncv=None):
      r"""Relaxation experiment.
  
      The relaxation experiment describes the time-evolution
      of an expectation value starting in a non-equilibrium
      situation.
  
      Parameters
      ----------
      P : (M, M) ndarray
          Transition matrix
      p0 : (M,) ndarray (optional)
          Initial distribution for a relaxation experiment
      obs : (M,) ndarray
          Observable, represented as vector on state space
      times : list of int (optional)
          List of times at which to compute expectation
      k : int (optional)
          Number of eigenvalues and amplitudes to use for computation
      ncv : int (optional)
          The number of Lanczos vectors generated, `ncv` must be greater than k;
          it is recommended that ncv > 2*k
  
      Returns
      -------
      res : ndarray
          Array of expectation value at given times
  
      """
      R, D, L = rdl_decomposition(P, k=k, ncv=ncv)
      """Extract eigenvalues"""
      ev = np.diagonal(D)
      """Amplitudes"""
      amplitudes = np.dot(p0, R) * np.dot(L, obs)
      """Propgate eigenvalues"""
      times = np.asarray(times)
      ev_t = ev[np.newaxis, :] ** times[:, np.newaxis]
      """Compute result"""
      res = np.dot(ev_t, amplitudes)
      """Truncate imgainary part - is zero anyways"""
      res = res.real
      return res