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)
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)``.