def METHOD_NAME(X):
    """Compute square of the operator 2-norm (spectral norm) of X.
    This corresponds to the Lipschitz constant of the gradient of the
    squared-loss function:
        w -> .5 * ||y - Xw||^2
    Parameters
    ----------
    X : ndarray, shape (n_samples, n_features)
      Design matrix.
    Returns
    -------
    lipschitz_constant : float
      The square of the spectral norm of X.
    """
    # On big matrices like those that we have in neuroimaging, svdvals
    # is faster than a power iteration (even when using arpack's)
    return linalg.svdvals(X)[0] ** 2