Here is an explicit calculation for the Gaussian ensembles of $n\times n$ matrices $W$ with distribution
$$P(W)\propto e^{-\tfrac{1}{2}\beta\,{\rm tr}\,W^2},$$
with $\beta=1$ if $W$ is real symmetric and $\beta=2$ if $W$ is complex Hermitian. The eigenvalues of $W$ are denoted by $w_1,w_2,\ldots w_n$. I consider two quantities:
$$e^{nF(w_i)}=\prod_{j=1, j\ne i}^n|w_i-w_j|\Rightarrow F(w_i)=n^{-1}\sum_{j=1,j\ne i}^n\ln|w_i-w_j|,$$
$$e^{n^2G}=\prod_{i=1}^n\prod_{j=1, j\ne i}^n|w_i-w_j|\Rightarrow G=n^{-2}\sum_{i=1}^n\sum_{j=1,j\ne i}^n\ln|w_i-w_j|.$$
For $n\gg 1$ we may replace sums by integrals, upon the replacements of $w_i\mapsto n^{1/2}\lambda(i/n)$ and of
$$\sum_i u(w_i)\mapsto n\int dx\, u(n^{1/2}\lambda(x))=n\int d\lambda \,\rho(\lambda)u(n^{1/2}\lambda),\;\;\text{with}\;\;\rho(\lambda)=(d\lambda/dx)^{-1}.$$
We thus have the expressions
$$F(n^{1/2}\lambda)=\int \rho(\mu)d\mu\,\ln|n^{1/2}\lambda-n^{1/2}\mu|,$$
$$G=\int\rho(\lambda)d\lambda\int\rho(\mu)d\mu\,\ln |n^{1/2}\lambda-n^{1/2}\mu|.$$
The rescaled eigenvalue density $\rho(\lambda)=(d\lambda/dx)^{-1}$ is given in the large-$n$ limit by the Wigner semicircle),
$$\rho(\lambda)=\frac{1}{\pi}\sqrt{2-\lambda^2},\;\;\lambda\in(-\sqrt 2,\sqrt 2).$$
We thus arrive at
$$F(n^{1/2}\lambda)=\tfrac{1}{2}\ln n+\frac{1}{\pi}\int_{-\sqrt{2}}^{\sqrt{2}}d\mu\,\sqrt{2-\mu^2}\ln|\lambda-\mu|,$$
$$G=\tfrac{1}{2}\ln n-\tfrac{1}{4} -\tfrac{1}{2}\ln 2.$$
For $\lambda=0$ the integral over $\mu$ evaluates to
$$F(0)=\tfrac{1}{2}\ln n-\tfrac{1}{2} -\tfrac{1}{2}\ln 2.$$
Now to make contact with the OP, note that
$$f(\lambda_i)^{1/n}=\prod_{j=1, j\ne i}^n|\lambda_i-\lambda_j|^{1/n}=n^{-1/2}e^{F(n^{1/2}\lambda_i)},$$
hence I conclude that for $n\gg 1$ we have
$$f(\lambda)^{1/n}\rightarrow\exp\left[\frac{1}{\pi}\int_{-\sqrt{2}}^{\sqrt{2}}d\mu\,\sqrt{2-\mu^2}\ln|\lambda-\mu|\right],\;\;\lambda\in(-\sqrt 2,\sqrt 2).$$
For $\lambda=0$ this reduces to
$$f(0)^{1/n}=e^{-(1+\ln 2)/2}=0.428882.$$
I checked this numerically for the GOE, and it agrees (see plot below).
Data points are 40 random realizations of $f(\lambda_{n/2})^{1/n}$ for real symmetric matrices of size $n=800$, red line is the result $e^{-(1+\ln 2)/2}$ from the calculation.
Addendum: I note that the expectation value of the double product in the Gaussian ensembles can be computed exactly for any $n$ using Mehta's integral formula. For the case $\beta=1$ of real symmetric matrices one has
$$\mathbb{E}\left[\prod_{i=1}^n\prod_{j=1, j\ne i}^n|w_i-w_j|\right]=\int_{-\infty}^\infty dw_1\cdots\int_{-\infty}^\infty dw_n\, e^{-\tfrac{1}{2}\sum_i w_i^2}\prod_{i=1}^n\prod_{j=1,j\ne i}^n|w_i-w_j|^{3/2}$$
$$\qquad=\prod_{i=1}^n\frac{\Gamma(1+3i/2)}{\Gamma(1+3/2)}.$$
