Testing a Null Hypothesis of the form $H_0\colon d(X, Y)<ε$ Given the issues of normality testing (Is normality testing 'essentially useless'?), in particular for large sample sizes, I was wondering if there are any feasible "approximate" normality tests.
Typically, in Null Hypothesis testing, we test a hypothesis of the form $X = Y$, i.e. one tests if, given all the observed data, is there any credence that two distributions differ.
$$ (1) \qquad \boxed{H_0\colon X=Y} $$
This has the effect that in practice tests for normality almost always fail for large sample sizes. So, instead, I am interested in testing a Null Hypothesis of the form
$$ (2) \qquad \boxed{H_0\colon d(X, Y) < ε}$$
where $d$ is a suitable statistical divergence such as the Kullback–Leibler divergence or Wasserstein metric.
In particular, I am interested in automatically testing whether a large sample $X = (x_i)_{i=1:N}$ in $\mathbb R^d$ comes from a distribution close to a standard normal distribution. Here, the allowed closeness is defined by the design parameter $ε$.For example, we may need this to take into account numerical errors and other noise (PRNGs not producing perfect normal number etc).

Question: What are some numerically feasible test I could use?

Background: Creating a unit test for neural network architectures at initialization.
There is good empirical and theoretical evidence, that, in order to have stable training for deep neural networks, it is necessary that these networks satisfy a normalization condition at initialization; that is, we want that if $x∼(0,1)$ then $f(x, θ_0)∼(0,1)$, and this should hold layer- or block-wise. My goal is to build an automated probabilistic unit-test, that takes as input some neural network architecture and tests whether this condition holds approximately. Towards this goal I want to implement a function
def test_approx_normal(x: list[float], tol: float) -> bool:

that returns False if and only if the Null Hypothesis, i.e. that the sample comes from some distribution that is close to a standard normal distribution, is rejected, and True otherwise.
The standard Kolmogorov–Smirnov test does not seem suitable for this task, at least the following example shows $p$-values all over the place.
import numpy as np
from scipy.stats import kstest, norm as normal
from scipy.stats import norm as normal

for n in (10, 100, 1000, 10_000):
    x = normal.rvs(size=(n,))
    A = normal.rvs(size=(n, n)) / np.sqrt(n)
    y = A @ x
    print(n, kstest(x, normal.cdf), kstest(y, normal.cdf), sep="\n")

This is not surprising because this test is not "numerically robust" as it tests for equality between distributions (1) instead of approximate equality (2). However, we clearly see that the distribution follows a standard normal approximately:

 A: The Kullback-Leibler divergence isn't suitable here because it's designed for mutually absolutely continuous distributions: it doesn't make sense to compute the KL divergence between a discrete empirical distribution and an absolutely continuous reference distribution.
A significant obstacle with a Wasserstein metric would be computing the null distribution of the test statistic, which would be an infimum over an infinite-dimensional $\epsilon$ neighbourhood of the Normal distribution. As far as I know, this distribution isn't known.
If I were going to try this, I'd go for the supremum norm on cumulative distribution functions as the metric, because that's relatively well-studied and computationally simple.  Rather than a test of the null hypothesis that distance $>\epsilon$, I'd still probably want a confidence interval for the distance.
A: I think you can use:
$Z=\displaystyle\frac{\hat{\theta}-\theta}{se(\hat{\theta})} \rightarrow N(0,1)$ as $n \rightarrow \infty$.
where $se$ is the ESTIMATED standard deviation.
So I think you could:
$\theta=d(X,Y)-\epsilon$
$\hat{\theta}=\hat{d}(X,Y)-\epsilon$
where $\hat{d}$ is the estimated KL divergence, you can estimate this by using empiric distributions for both $X$ and $Y$ or using the known distribution (but unknown parameters) and to estimate the parameters for both densities and then evaluate KL divergence.
Finally you got:
$CI=(\hat\theta-z_{\alpha/2}.se(\hat{\theta}),\hat\theta+z_{\alpha/2}.se(\hat{\theta}))$
$CI=(\hat{d}(X,Y)-\epsilon -z_{\alpha/2}.se(\hat{\theta}),\hat{d}(X,Y)-\epsilon+z_{\alpha/2}.se(\hat{\theta}))$
then:
$\hat{d}(X,Y)-\epsilon -z_{\alpha/2}.se(\hat{\theta})\leq d(X,Y)-\epsilon \leq \hat{d}(X,Y)-\epsilon+z_{\alpha/2}.se(\hat{\theta})$
$\hat{d}(X,Y) -z_{\alpha/2}.se(\hat{\theta})\leq d(X,Y) \leq \hat{d}(X,Y)+z_{\alpha/2}.se(\hat{\theta})$
so you look over how $\epsilon$ would lie in it. I'd compute $se(\hat\theta)$ by bootstrap.
DISCLAIMER: I'm not sure if $Z$ distribution will hold for KL divergence.
