def transform_eigenspace(eigvecs, eigvals, transform_function):
    """Return a function that multiplies a vector by a matrix with
    transformed eigenvalues.

    Let ``eigvecs`` and ``eigvals`` be selected eigenvectors and
    eigenvalues.  The columns of ``eigvecs``
    are the eigenvectors associated with the corresponding entry of
    ``eigvals``.  A function, ``a_mult``, is returned.  This function
    is the identity for all directions orthogonal to ``eigvecs``,
    and has eigenvalues ``transform_function(eigvals)`` in the
    space spanned by ``eigvecs``.

    Parameters
    ------------
    eigvecs: `numpy.ndarray` (N, K)
        The eigenvectors.
    eigvals: `numpy.ndarray` (K,)
        The eigenvalues
    transform_function: callable
        A function from ``eigvals`` to a vector of the same length.
        The output of ``transform_function(eigvals)`` will be the new
        eigenvalues.

    Returns
    -----------
    a_mult: callable
        A linear function from a length-``K`` numpy vector to another
        length-``K`` vector with the above-described eigendecomposition.
    """

    if eigvecs.ndim != 2:
        raise ValueError('``eigvecs`` must be 2d.')
    if eigvals.ndim != 1:
        raise ValueError('``eigvals`` must be 1d.')
    if eigvecs.shape[1] != len(eigvals):
        raise ValueError(
            'The columns of ``eigvecs`` and length of ``eigvals`` must match.')

    new_eigvals = transform_function(eigvals)

    def a_mult(vec):
        vec_loadings = eigvecs.T @ vec
        # Equivalent to the more transparent:
        # vec_perp = vec - eigvecs @ vec_loadings
        # return vec_perp + eigvecs @ (new_eigvals * vec_loadings)
        return vec + eigvecs @ ((new_eigvals - 1) * vec_loadings)

    return a_mult