docs: def economic_qs(K, epsilon=sqrt(finfo(float).eps)):
      r"""Economic eigen decomposition for symmetric matrices.
  
      A symmetric matrix ``K`` can be decomposed in
      :math:`\mathrm Q_0 \mathrm S_0 \mathrm Q_0^\intercal + \mathrm Q_1\
      \mathrm S_1 \mathrm Q_1^ \intercal`, where :math:`\mathrm S_1` is a zero
      matrix with size determined by ``K``'s rank deficiency.
  
      Args:
          K (array_like): Symmetric matrix.
          epsilon (float): Eigen value threshold. Default is
                           ``sqrt(finfo(float).eps)``.
  
      Returns:
          tuple: ``((Q0, Q1), S0)``.
      """
  
      (S, Q) = eigh(K)
  
      nok = abs(max(Q[0].min(), Q[0].max(), key=abs)) < epsilon
      nok = nok and abs(max(K.min(), K.max(), key=abs)) >= epsilon
      if nok:
          from scipy.linalg import eigh as sp_eigh
  
          (S, Q) = sp_eigh(K)
  
      ok = S >= epsilon
      nok = logical_not(ok)
      S0 = S[ok]
      Q0 = Q[:, ok]
      Q1 = Q[:, nok]
      return ((Q0, Q1), S0)