One way to see this is to consider the fact that the vector space of rank (3,3) completely antisymmetric tensors ($ \Lambda_3^3(R^3) $) has dimension one (it's just a linear algebra exercise). Then define the tensor:
$$
M_{ijk}^{lmn} = \delta_i^{[l} \, \delta_j^m \, \delta_k^{n]} = \frac{1}{3!} \sum_{\sigma \in S_3} sgn(\sigma) \, \delta_i^{\sigma(l)} \, \delta_j^{\sigma(m)} \, \delta_k^{\sigma(n)}
$$
where we are summing over all the permutations of three numbers $\sigma$, and $sgn(\sigma)$ denotes the sign of the permutation. It's worthy to note that
$$
M_{ijk}^{lmn}=\frac{1}{3!}\begin{vmatrix} \delta_i^l & \delta_i^m & \delta_i^n\\
\delta_j^l & \delta_j^m & \delta_j^n \\ \delta_k^l & \delta_k^m & \delta_k^n \end{vmatrix}
$$
by the Leibniz formula of the determinant (http://en.wikipedia.org/wiki/Leibniz_formula_for_determinants). So we have $M \in \Lambda_3^3(R^3) $. Since $M \neq 0$, $M$ is a basis for the space $ \Lambda_3^3(R^3) $. Now consider the tensor
$$
\epsilon_{ijk} \, \epsilon^{lmn} = B_{ijk}^{lmn}
$$
Since $B \in \Lambda_3^3(R^3) $ and $M$ is a basis, there exists a constant $k$ such that
$$
B_{ijk}^{lmn} = k \, M_{ijk}^{lmn} \implies \epsilon_{ijk} \, \epsilon^{lmn} = k \, \delta_i^{[l} \, \delta_j^m \, \delta_k^{n]}
$$
Now to determine $k$ contract $\epsilon_{lmn}$ on both sides and use the fact that $\epsilon_{lmn} \, \epsilon^{lmn} =3!$ (since you sum $3!$ terms equal to one)
$$
\epsilon_{ijk} \, \epsilon^{lmn} \, \epsilon_{lmn} = 3! \, \epsilon_{ijk} = k \, \delta_i^{[l} \, \delta_j^m \, \delta_k^{n]} \epsilon_{lmn} = k \, \delta_i^{l} \, \delta_j^m \, \delta_k^{n} \, \epsilon_{[lmn]} = k \, \delta_i^{l} \, \delta_j^m \, \delta_k^{n} \, \epsilon_{lmn} = k \, \epsilon_{ijk}
$$
So we finally get $k=3!$ and
$$
\epsilon_{ijk} \, \epsilon^{lmn}=\begin{vmatrix} \delta_i^l & \delta_i^m & \delta_i^n\\
\delta_j^l & \delta_j^m & \delta_j^n \\ \delta_k^l & \delta_k^m & \delta_k^n \end{vmatrix}
$$
And you can get the identity you want by contracting. The same argument could be used in any dimension to show that
$$
\epsilon_{i_1,\dots,i_n} \, \epsilon^{j_1,\dots,j_n} = \begin{vmatrix} \delta_{i_1}^{j_1} & \dots & \delta_{i_1}^{j_n}\\
\vdots & & \vdots \\ \delta_{i_n}^{j_1} & \dots & \delta_{i_n}^{j_n} \end{vmatrix}
$$
Don't use indices. Use a sophisticated and efficient notation for tensor manipulations: through clifford algebra.
Let $\gamma^0, \gamma^1, \gamma^2, \gamma^3$ be basis vectors. Under the clifford product operation, they obey the following:
$$\gamma^\mu \gamma^\nu = \begin{cases}\eta^{\mu \nu} & \mu = \nu \\ -\gamma^\mu \gamma^\nu & \mu \neq \nu\end{cases}$$
The product operation is also associative.
Write $E$ and $B$ in terms of $\gamma^1, \gamma^2, \gamma^3$ as you would write any vectors: $E = E_x \gamma^1 + \ldots$ for instance. Then $F$ can be written
$$F = \gamma^0 E - \gamma^1 \gamma^2 \gamma^3 B$$
Define $\epsilon = -\gamma^0 \gamma^1 \gamma^2 \gamma^3$. Then the tensor equation $\epsilon^{\alpha \beta \gamma \delta} F_{\alpha \beta} F_{\gamma \delta}$ can be written
$$\epsilon^{\alpha \beta \gamma \delta} F_{\alpha \beta} F_{\gamma \delta} = \langle \epsilon FF\rangle_0$$
The brackets denote that we're looking for the scalar part of this expression. This is equivalent to $\epsilon \langle FF \rangle_4 = \epsilon (F \wedge F)$ in the common shorthand.
Now we can use the associativity of the product to our advantage. For brevity, let $i = \gamma^1 \gamma^2 \gamma^3$, such that $\gamma^0 i = \epsilon$, and see that
$$\langle \epsilon FF \rangle_0 = \langle \epsilon (\gamma^0 E - iB)(\gamma^0 E - iB) \rangle_0$$
Multiply out the terms:
$$\langle \epsilon FF \rangle_0 = \langle \epsilon (\gamma^0 E \gamma^0 E -\gamma^0 E iB - iB \gamma^0 E + iBiB )\rangle_0$$
Use the properties of the product to deduce that $E\gamma^0 = -\gamma^0 E$, $\gamma^0 B = -B \gamma^0$, $Ei = iE$ and $iB = Bi$, and $ii = 1$. The result reduces to
$$\langle \epsilon FF \rangle_0 = \langle \epsilon (-EE - 2 \epsilon EB + BB) \rangle_0$$
Finally, deduce that $\epsilon \epsilon = -1$, $EE = -|E|^2$ and $BB = -|B|^2$ to get
$$\langle \epsilon FF \rangle_0 = \langle \epsilon (|E|^2 - |B|^2) + 2 E \cdot B \rangle_0$$
The only term that survives the angled brackets is the $2 E \cdot B$ term. Of course, staring us right in the face is the other well-known invariant, $E^2 - B^2$. That is one advantage of this method: both invariants can be derived from the same approach, and at the same time. Moreover, we never have to break into components to evaluate the expression. We merely use the commutation properties of basis elements in the clifford algebra--an algebra already widely used in quantum mechanics.
Best Answer
Following the method in L&L vol. 2 ch. 1, the left-hand side $$\varepsilon^{\alpha \beta \mu \nu} \varepsilon_{\alpha \beta \rho \sigma}$$ is a product of a pseudo-tensor, which is invariant under Lorentz transformations up to a factor of a determinant, and a pseudo-tensor which transforms under the inverse determinant.
This tells us that the right-hand side has to be an invariant tensor and so must be constructed from Kronecker delta's.
The right hand side must also be anti-symmetric in $\mu,\nu$ and anti-symmetric in $\rho,\sigma$ so that on general grounds $$ \varepsilon^{\alpha \beta \mu \nu} \varepsilon_{\alpha \beta \rho \sigma} = - A (\delta^{\mu}_{\rho} \delta^{\nu}_{\sigma} - \delta^{\mu}_{\sigma} \delta^{\nu}_{\rho})$$ must hold, where the minus sign is due to the Minkowski metric, and the factor $A=2$ is fixed by the requirement that $\varepsilon^{\mu \nu \rho \sigma} \varepsilon_{\mu \nu \rho \sigma} = - 4!$ holds, so setting $ - A (\delta^{\mu}_{\mu} \delta^{\nu}_{\nu} - \delta^{\mu}_{\nu} \delta^{\nu}_{\mu}) = - 4!$ gives $A(16 - 4) = 24$ or $A = 2$.
More generally, the same thinking tells us that the product $\varepsilon^{\alpha \beta \mu \nu} \varepsilon_{\alpha \beta \rho \sigma}$ must be an anti-symmetric combination of Kronecker delta's up to an overall normalization factor, where the normalization factor must respect the fact that $\varepsilon^{\mu \nu \rho \sigma} \varepsilon_{\mu \nu \rho \sigma} = - 4!$ should hold, but $\delta^{\mu \nu \rho \sigma}_{\alpha \beta \gamma \delta}$ is an anti-symmetric combination of Kronecker delta's such that $\delta^{\mu \nu \rho \sigma}_{\mu \nu \rho \sigma} = 4!$, and so $$ \varepsilon^{\alpha \beta \mu \nu} \varepsilon_{\alpha \beta \rho \sigma} = - \delta^{\mu \nu \rho \sigma}_{\alpha \beta \rho \sigma}$$ must hold. Thus we can work out contractions like $\varepsilon^{\alpha \beta \mu \nu} \varepsilon_{\alpha \beta \rho \sigma} = - \delta^{\alpha \beta \mu \nu}_{\alpha \beta \rho \sigma}$ directly \begin{align*} \varepsilon^{\alpha \beta \mu \nu} \varepsilon_{\alpha \beta \rho \sigma} &= - \delta^{\alpha \beta \mu \nu}_{\alpha \beta \rho \sigma} \\ &= - (\delta^{\alpha}_{\alpha} \delta^{\beta \mu \nu}_{\beta \rho \sigma} - \delta^{\alpha}_{\beta} \delta^{\beta \mu \nu}_{\alpha \rho \sigma} + \delta^{\alpha}_{\rho} \delta^{\beta \mu \nu}_{\alpha \beta \sigma} - \delta^{\alpha}_{\sigma}\delta^{\beta \mu \nu}_{\alpha \beta \rho}) \\ &= - (4\delta^{\beta \mu \nu}_{\beta \rho \sigma} - \delta^{\alpha \mu \nu}_{\alpha \rho \sigma} + \delta^{\beta \mu \nu}_{\rho \beta \sigma} - \delta^{\beta \mu \nu}_{\sigma \beta \rho}) \\ &= - [4(\delta^{\beta}_{\beta} \delta^{\mu \nu}_{\rho \sigma} - \delta^{\beta}_{\rho} \delta^{\mu \nu}_{\beta \sigma} + \delta^{\beta}_{\sigma}\delta^{\mu \nu}_{\beta \rho}) - ( \delta^{\alpha }_{\alpha} \delta^{\mu \nu}_{\rho \sigma} - \delta^{\alpha}_{\rho} \delta^{\mu \nu}_{\alpha \sigma} + \delta^{\alpha}_{\sigma} \delta^{\mu \nu}_{\alpha \rho}) - \delta^{\beta \mu \nu}_{\beta \rho \sigma} + \delta^{\beta \mu \nu}_{\beta \sigma \rho}] \\ &= - [4(4\delta^{\mu \nu}_{\rho \sigma} - \delta^{\mu \nu}_{\rho \sigma} + \delta^{\mu \nu}_{\sigma \rho}) - ( 4 \delta^{\mu \nu}_{\rho \sigma} - \delta^{\mu \nu}_{\rho \sigma} + \delta^{\mu \nu}_{\sigma \rho}) - \delta^{\beta \mu \nu}_{\beta \rho \sigma} + \delta^{\beta \mu \nu}_{\beta \sigma \rho}] \\ &= - [4(2\delta^{\mu \nu}_{\rho \sigma} ) - 2 \delta^{\mu \nu}_{\rho \sigma} - 2 \delta^{\mu \nu}_{\rho \sigma} + 2 \delta^{\mu \nu}_{\sigma \rho}] \\ &= - [8 \delta^{\mu \nu}_{\rho \sigma} - 6 \delta^{\mu \nu}_{\rho \sigma} ] \\ &= - 2 \delta^{\mu \nu}_{\rho \sigma} \\ &= - 2 (\delta^{\mu}_{\rho} \delta^{\nu}_{\sigma} - \delta^{\mu}_{\sigma} \delta^{\nu}_{\rho}). \end{align*} Similarly, in $d$ dimensions we have $$\varepsilon^{\mu_1 .. \mu_d} \varepsilon_{\mu_1 .. \mu_d} = - d!$$ From this one immediately sees that $$ \varepsilon^{\mu_1 .. \mu_d} \varepsilon_{\nu_1 .. \mu_d} = - \delta^{\mu_1 .. \mu_d}_{\nu_1 .. \nu_d} $$ and so contractions obey identities like $$\varepsilon^{\mu_1 \ldots \mu_r \mu_{r+1} .. \mu_d} \varepsilon_{\mu_1 .. \mu_r \nu_{r+1} \ldots \nu_d} = - A \delta^{\mu_{r+1} .. \mu_d}_{\nu_{r+1} .. \nu_d}$$ must hold in $d$ dimensions, where $A$ can be fixed by expanding $$\varepsilon^{\mu_1 \ldots \mu_r \mu_{r+1} .. \mu_d} \varepsilon_{\mu_1 .. \mu_r \nu_{r+1} \ldots \nu_d} = - \delta^{\mu_1 .. \mu_r \mu_{r+1} .. \mu_d}_{\mu_1 .. \mu_r \nu_{r+1} .. \nu_d}$$ one step at a time as in the example above, and obviously $A = r!$ should hold so that $A = - d!$ when $r = d$, as can be proven by induction.