I am looking at the treatment of the chiral anomaly in Fujikawa and Suzuki's "Path Integrals and Quantum Anomalies." To illustrate the quantum breaking of chiral symmetry (section 4.3), they start with the QED Lagrangian and impose a chiral transformation on the fields, as in
\begin{align} \psi(x) &\mapsto \psi^{\prime}(x)=e^{i \alpha(x) \gamma_{5}} \psi(x), \notag\\\label{EQawdwd22ee2} \bar\psi(x) &\mapsto\bar{\psi}^{\prime}(x)=\bar{\psi}(x) e^{i \alpha(x) \gamma_{5}}. \end{align}
Under the assumption of invariance of the path integral and the measure under this transformation, one obtains the ''naive'' chiral identity
\begin{align}\label{EQ23ddww} \partial_{\mu}\left\langle\bar{\psi}(x) \gamma^{\mu} \gamma_{5} \psi(x)\right\rangle= 2 m i\left\langle\bar{\psi}(x) \gamma_{5} \psi(x)\right\rangle. \end{align}
Now is where my confusion begins. Fujikawa claims that this identity does not hold in Lorentz invariant perturbation theory. He claims that one can use covariant regularization of the current to show the deviation from the naive chiral identity, avoiding all perturbative calculations. In particular, he goes on to calculate
\begin{align}\label{EQawd2f3f} \left\langle j_5^{\mu}\right\rangle_{\mathrm{cov}} = \left\langle\bar{\psi}(x) \gamma^{\mu} \gamma_{5} \psi(x)\right\rangle_{\mathrm{cov}} :=\lim _{y \rightarrow x} \operatorname{Tr} \gamma^{\mu} \gamma_{5} \frac{1}{i \not D-m} f(\not D^{2} / \Lambda^{2}) \delta^{(4)}(y-x), \end{align} where the trace is a sum over Dirac indices, and we use the (hermitian) Dirac operator in the Euclidean theory \begin{align} \not D &:=\gamma ^\mu D_\mu\\ % % & =\gamma^{\mu}\left(\partial_{\mu}-i e A_{\mu}\right), \end{align} where $D_{\mu}:=\partial_{\mu}-i e A_{\mu}$ is the gauge covariant derivative. The regulator function $f(x)$ is taken to be an arbitrary smooth function $f(x)$ which vanishes sufficiently rapidly at large values of $x$;
\begin{align}\label{EQd2ddddd11d1} f(0)=1, \quad \lim _{x \rightarrow \infty} f(x)=0. \end{align}
I have no intuition as to where this procedure comes from and why it might work, or what it does. Furthermore, why exactly is there a need to regularize? My understanding is that regularization is applied to remedy unphysical divergences in our theory, and so its not presently clear to me why we do this.
I also don't see how this regularization procedure may implement gauge invariance / covariance as implied by the name.
On a more general note, my understanding is that anomalies are defined as symmetries in the classical action that do not carry over to the full quantum theory. Therefore, I am curious as to how the regularization process outlined above implements this upgrade from the classical to the quantum theory, so as to bring the anomalous non-conservation of chiral current to surface. For reference, the final result is
\begin{align} \partial_{\mu}\left\langle\bar{\psi}(x) \gamma^{\mu} \gamma_{5} \psi(x)\right\rangle_{\mathrm{cov}}&=2 i m\left\langle\bar{\psi}(x) \gamma_{5} \psi(x)\right\rangle_{\mathrm{cov}}\\ % % &\quad +2 i\left(\frac{e^{2}}{32 \pi^{2}}\right) \epsilon^{\mu \nu \alpha \beta} F_{\mu \nu} F_{\alpha \beta}. \end{align}
with the anomalous term in the second line.