1. Instantons
1.1 Instantons as a classical solution
An instanton, is pretty much exactly what you say: An (anti-)self-dual configuration of the curvature of a principal bundle. The curvature $F_A$ of a principal bundle with connection $A$ is the field strength tensor of the physical gauge theory, while the connection is called the gauge potential. In electromagnetism, for example, the non-zero components of $F$ are exactly the electric and magnetic fields, so it has direct physical relevance. In general, the action1 of a Yang-Mills gauge theory is given by the integral
$$ S_\text{YM}[A] = \int_M\mathrm{Tr}_\mathfrak{g}(F\wedge \ast F)$$
whose Euler-Lagrange equations are the classical equations of motion, i.e. the classical solutions are stationary points of this functional. Now, we can decompose the field strength into a self-dual part $F^+$ and an anti-self-dual part $F^-$ which are orthogonal to each other w.r.t. the inner product
$$ (G,H) = \int G\wedge\ast H$$
Plugging this in gives
$$ S_\text{YM}[A] = \int \mathrm{Tr}_\mathfrak{g}(F^+\wedge \ast F^+) + \int \mathrm{Tr}_\mathfrak{g}(F^-\wedge \ast F^-)$$
and comparing this to the second Chern class $C_2(A) := \int\mathrm{Tr}_\mathfrak{g}(F\wedge F)$ we can see that $S_\text{YM}[A] \geq \lvert C_2(A)\rvert $, i.e. is locally minimized when equality holds. But that equality holds exactly when either $F^+ = 0$ or $F^- = 0$, i.e. when the full field strength $F$ is itself either self-dual or anti-self-dual (or vanishing). Thus, the classical equations of motion are equivalent to $\ast F = \pm F$, i.e. instantons are just classical solutions to the equation of motion.
Already classically, the size of the moduli space of instantons is somewhat interesting, since it shows whether or not there are possible different solutions to the classical equation of motion. However, instantons may be related by large gauge transformations which one may quotient out classically, so counting the instantons as such is not enough information to be interesting.
1.2 Instantons as "vacua" of the quantum theory
Standard perturbation theory in quantum field theory proceeds by expanding the exponential of the action function around a classical solution. The existence of multiple classical solutions means we should expand around every single one of them, and sum up the solutions with some weight $f$.2 If we denote by $\mathcal{D}A_k$ a Feynman path integral measure (with gauge equivalence classes quotiented out) that is over perturbations of the instanton labeled by $k = \frac{1}{8\pi}C_2(A)$3, we get for the Euclidean path integral
$$ Z = \sum_k f(k)\int \mathrm{e}^{-S_\text{YM}[A_k]}\mathcal{D}A_k$$
By a standard decomposition argument that this integral should factor into local parts when we consider things like $A_k = A_{k_1} + A_{k_2}$ where $A_{k_1},A_{k_2}$ live on different regions, one heuristically finds that $f(k_1+k_2) = f(k_1)f(k_2)$ should hold, and since everything in a physicist's world is smooth except when it isn't, this means that $f$ is itself an exponential $f(k) = \mathrm{e}^{\mathrm{i}\theta k}$ for some $\theta\in\mathbb{R}$. Remember how $k$ was defined, we get a modified action
$$ S_\theta[A] = \int\mathrm{Tr}_\mathfrak{g}(F\wedge\ast F) + \frac{\mathrm{i}\theta}{8\pi^2}\int\mathrm{Tr}_\mathfrak{g}(F\wedge F)$$
with which we can drop the sum and write
$$ Z = \int \mathrm{e}^{-S_\theta[A]}\mathcal{D}A$$
where $\mathcal{D}A$ now ranges over the whole space of connections.
From this a development of other theories may follow, e.g. one may now try to promote the $\theta$ itself to a dynamical field, then called the axion in Peccei-Quinn theory.
However, there is an issue when one tries to do this sum in theories coupled to fermions, which is the appearence of a quantum anomaly.
1.3 Instantons and quantum anomalies
We consider a theory of quantum fermion fields coupled to a classical gauge field background (i.e. an instanton)
$$ S_\text{Dirac}[\psi,A] = \bar\psi (\mathrm{i}D - m)\psi - \frac{1}{4}F\wedge\ast F$$
where $D$ is the Dirac operator $D = D_\mu\gamma^\mu$ with $D_\mu$ the covariant derivative belonging to $A$ and $\gamma^\mu$ the usual gamma matrices acting on the spinors $\psi$. We would like to write
$$ Z[A] = \int\mathcal{D}\psi\mathcal{D}\bar\psi \mathrm{e}^{-S_\text{Dirac}[\psi,A]}$$
but there is an issue with the usual definition of $\mathcal{D}\psi\mathcal{D}\bar\psi$ in terms of a limit of integrations over the modes of the field. See this answer of mine for the definition and regularization of the measure and this answer of mine for how, in the end, the anomaly - i.e. the problem in defining the measure in an invariant way under large gauge transformations - is related to the index of $D$, which by the Atiyah-Singer index theorem is closely related to the instanton number $k$ as the Chern class.
1Actually, we consider the so-called Euclidean action, to which the actual action is related by Wick rotation
2Perturbatively, the states belonging to such different sectors are unrelated, but with the theory of instantons, one can see that there are non-perturbative amplitudes between them, see this answer of mine
3There may be multiple instantons with the same number, the sum is then over each on of them (since in general it is not guaranteed those are related by a gauge transformation)
2. Donaldson invariants
This section attempts to be a brief retelling of Witten's "Topological Quantum Field Theory"
2.1 Twisted $\mathcal{N}=2$ supersymmetric Yang-Mills theory
The physical setting in which the Donaldson invariants will appear is a Yang-Mills theory living on the four-dimensional manifold $M$ coupled to certain fields such that the total action enjoys a supersymmetry. The action of this theory looks admittedly horrible:
$$ S_\text{SYM} = \int_M \mathrm{Tr}_\mathfrak{g}\left(\frac{1}{4}F_{\mu\nu}F^{\mu\nu} + \frac{1}{4}F_{\mu\nu}(\ast F)^{\mu\nu} + \frac{1}{2}\phi D_\mu D^\mu \lambda - \mathrm{i}\eta D_\mu \psi^\mu + \mathrm{i}D_\mu\psi_\nu \chi^{\mu\nu} - \frac{\mathrm{i}}{8}\phi[\chi_{\mu\nu},\chi^{\mu\nu}] - \frac{\mathrm{i}}{2}\lambda[\psi_\mu,\psi^\mu] - \frac{\mathrm{i}}{2}\phi[\eta,\eta] - \frac{1}{8}[\phi,\lambda]^2 \right)\sqrt{g}\mathrm{d}^4 x$$
Here $D_\mu = \nabla_\mu + [A_\mu,\dot{}]$ is the gauge covariant derivative (with $\nabla$ the ordinary Riemannian covariant derivative) and all fields are $\mathfrak{g}$-valued. There is a $\mathbb{Z}_2$-grading on the space of fields, where we call the different classes "bosonic" (or even) and "fermionic" (or odd). The bosonic fields are $\phi,\lambda$, the fermionic ones are $\eta,\psi,\chi$, and $\chi$ is additionally constrained to be self-dual. This action is invariant under the infinitesimal symmetry
\begin{align}
\delta_\epsilon A = \mathrm{i}\epsilon\psi \quad \delta_\epsilon\phi = 0 \quad \delta_\epsilon\lambda = 2\mathrm{i}\epsilon\eta \\
\delta_\epsilon\psi = -\epsilon D\phi \quad \delta_\epsilon\eta = \frac{1}{2}\epsilon[\phi,\lambda] \quad \delta_\epsilon \chi = \epsilon(F + \ast F) \tag{1}
\end{align}
with $\epsilon$ a fermionic infinitesimal parameter. As with all transformations, we think of this one as having a generator - its supercharge - $Q$, which gives all transformations as $\delta \alpha = -\mathrm{i}\epsilon Q(\alpha)$ where $\alpha$ is any field. In a Hamiltonian formulation, $Q(\alpha)$ would be the Poisson bracket $\{Q,\alpha\}$, but on a general manifold, we don't have that option. By explicit computation, one finds that $$\delta_\epsilon\delta_\zeta X - \delta_\zeta\delta_\epsilon X = -2\mathrm{i}\epsilon\zeta\phi$$
for every field $X$ except $A$, where it is $-D \phi$. This holds only on-shell for $\chi$, but off-shell for all others. Therefore, the commutator of two such transformations is a gauge transformation, and hence has no physical impact. In the Hamiltonian formulation, we would write this $\{Q,Q\} = 0$ (modulo gauge transformations), that is, $Q$ is akin to a BRST charge. The conserved current associated to this symmetry is
$$ J = \mathrm{Tr}_\mathfrak{g}\left((F_{\mu\nu} + (\ast F)_{\mu\nu})\psi^\nu - \eta D_\mu \phi - D^\nu\phi \chi_{\mu\nu} - \frac{1}{2}\psi_\mu[\lambda,\phi]\right)\mathrm{d}x^\mu$$
where conservation means that $\ast J$ is closed, so for any homology 3-cycle $\gamma$, the integral
$$ Q[\gamma] = \int_\gamma \ast J $$
depends only on the homology class of $\gamma$. Furthermore, one may show that the energy-momentum tensor $T_{\mu\nu} = 2 \frac{\delta S_\text{SYM}}{\delta g^{\mu\nu}}$ of this theory is an infinitesimal transform $T_{\mu\nu} = \{Q,\lambda_{\mu\nu}\}$ for another ugly expression $\lambda$ (see Witten's eq. (2.34)).
2.2 Donaldson invariants as path integrals
In the following, the path integral measure $\mathcal{D}X$ includes all fields, and also intends to have gauge equivalence classes quotiented out. The generic object we consider is the (unnormalized) expectation value of any observable $\mathcal{O}$, where $\mathcal{O}$ is any nice functional in the fields:
$$ Z(\mathcal{O}) = \int \mathcal{O}\mathrm{e}^{-S_\text{SYM}/e^2}\mathcal{D}X $$
If the supersymmetry transformation is non-anomalous, we have $Z(\{Q,\mathcal{O}\}) = 0$ for every observable. We now claim that $Z = Z(1)$ is a smooth invariant, and in particular will turn out to be a Donaldson invariant. For $Z$ to be a smooth invariant, it must be invariant under changes in the metric. The change of the action under a change of metric is by definition $\delta S = \frac{1}{2}\int_M \sqrt{g} T_{\mu\nu}\delta g^{\mu\nu}$ and this leads to
$$\delta Z(1) = -\frac{1}{e^2}Z(\{Q,\int_M \sqrt{g}\delta g^{\mu\nu}\lambda_{\mu\nu}\}) = 0 $$
so $Z(1)$ is invariant under changes of the metric. Similarily, it is invariant under changes of the gauge coupling constant $e$, as long as it stays non-zero. But in the limit of small coupling, the path integral is strongly dominated by classical minima of the free theories...and the classical minima of the free gauge theory are the anti-self-dual instantons! The self-dual ones are not minima because we added the topological $F\wedge F$ to the Lagrangian. So we may evaluate $Z$ by looking at the instanton contributions. However, as in 1.3 above, the path integral measure is not invariant if the fermionic zero modes are mismatched. The equations for the $\psi$ zero modes turn out to be the exact same equation as that for an infinitesimal perturbation $\delta A$ to an instanton configuration $A$ to be an instanton:
$$ D_\mu Y_\nu - D_\nu Y_\mu + \epsilon_{\mu\nu\sigma\rho}D^\sigma Y^\rho \quad\land\quad D_\mu Y^\mu = 0$$
for $Y$ either $\delta A$ or $\psi$. But the number of possible independent perturbations of an instanton that are again an instanton is exactly what one would call the dimension $\dim(\mathcal{M})$ of the moduli space at the point $A$ if it was a proper smooth non-singular space, so the number of $\psi$ zero modes is the same as the dimension of the moduli space of instantons. In general, there seems to be an index theorem that says the total number of zero modes is equal to the formal dimension of the moduli space, but Witten does not give its name or application here.
Specializing to the situation in the question, where the instanton moduli space is discrete and has dimension zero, we thus have that $Z$ is a smooth invariant. For a fixed instanton background and in the limit of weak coupling, it is enough to look at the lowest order terms in the fields. Those are quadratic terms of the form $\Phi \Delta \Phi$ and $\Psi \mathscr{D} \Psi$ where $\Delta$ is a second order elliptic operator on the bosonic fields $\Phi = (A,\phi,\lambda)^T$ and $\mathscr{D}$ is first order real skew-symmetric on the fermionic fields $\Psi = (\eta,\psi,\chi)^T$. This means the path integral over $\Phi$ and $\psi$ degenerates into a Gaussian integral.
Furthermore, their eigenvalues are related by supersymmetry: Looking at the supersymmetry transformation $(1)$, we see that the classical solutions $F = -\ast F$ and $\phi,\lambda,\eta,\psi,\chi = 0$ are invariant under supersymmetry, so the quantum excitations when expanding about them are related by supersymmetry, too. For each non-zero $\lambda$ that is an eigenvalue of $\mathscr{D}$ (which come in pairs since it is skew-symmetric), there is an eigenvalue $\lambda^2$ of $\Delta$.This is shown in D'Adda, DiVecchia, "Supersymmetry and instantons"
Now, the Gaussian path integral over $\int_M \left(\Phi \Delta \Phi + \mathrm{i}\Psi \mathscr{D} \Psi \right)\sqrt{g}$ yields $\mathrm{Pf}(\mathscr{D}) / \sqrt{\det(\Delta)}$, and the Pfaffian and the determinant differ only by a sign, so this becomes $\pm \prod_i \mathrm{sgn}(\lambda_i)$, where the product runs over all non-zero eigenvalues. The $\pm$ is due to us having to pick an orientation that determines the sign of the Pfaffian. So we pick any one instanton $A_0$ and declare this product is $+1$ for it. Now, one picks any other instanton $A_i$ and determines how often $\mathscr{D}$ has a zero eigenvalue along the homotopy $A_t = t A_0 + (1-t) A_i$. Everytime it obtains a zero eigenvalue, the sign of $\mathrm{Pf}(\mathscr{D})$ is defined to change. (I think this is basically transporting $\mathscr{D}$ along the curve $A_t$ in the space $\mathcal{A}/\mathcal{G}$, where $\mathcal{A}$ is the space of connections.) This gives a well-defined way to define the sign of the Pfaffian for $A_i$ if it gives the same sign regardless of which homotopy $A_t$ is chosen - and since concatenating two of them gives a transformation of $A_0$ onto itself, this is equivalent to the requirement that the sign must not change along any loop in $\mathcal{A}/\mathcal{G}$ based at $A_0$.
If this is given, then we obtain that $Z = \sum_i(-1)^{n_i}$, where the sum is over all instantons and the $n_i$ are determined in the above fashion. This is now, finally, exactly the same sketchily defined-invariant as in the question.
