I used different tools to do differently expressed genes analysis. From the results I see a general difference in logFC values for genes from limma-voom, edgeR and DESeq2. Although some people saw difference in logFC between edgeR and EDSeq2 and got good answers (seeBig differences in logfoldchange for some genes between edgeR and DESeq2 ), in my case edgeR and DESeq2 get very simmilar results but quite different from limma-voom.
There are a number of reasons why the different packages will give different fold change estimates. It's pretty clear what the culprit is in your case, but I'll list all the possible reasons before I explain that:
Normalization. You have to apply the same normalization method and fit the same model in each package, otherwise the different packages will be estimating different things.
Arithmetic. There is an arithmetic difference between the linear model fold changes computed by limma vs the negative binomial generalized linear model fold changes computed by edgeR and DESeq2. This is the effect that you asked about and which Aaron has discussed. As Aaron explained, the differences are usually either small or not important.
Shrinkage. None of the three packages report "raw" fold changes (unless you ask them to). All three packages instead shrink the log fold changes towards zero and the default amount of shrinkage is different in the three packages. The defaults in limma and edgeR are similar, with limma doing a bit more shrinkage than edgeR. The shrinkage that edgeR does is user-settable although I assume you have used the default. DESeq2 might do a lot more shrinkage if you ask it to estimate the amount of shrinkage to do. This is the effect that was discussed in the previous question & answer that you gave the link to.
Dispersion. All the packages estimate the dispersions somewhat differently and this also has an influence.
The most important factors, the ones that can potentially give big differences, are (1) and (3). In your case it appears that the culprit is (1). Your log fold changes from limma are not shrunk (closer to zero) compared to edgeR and DESeq2, but rather are substantially shifted (more negative, with smaller positive values and larger negative values). Such a shift in the logFCs can only occur because you have normalized differently or fitted a different model in limma compared to the other two packages.
In summary, the differences you see are not intrinsic to the packages but have almost certainly resulted from the way that you used the packages.
I don't think the differences are due to the two formulas you give. The computations done by the three packages are actually more complicated than these formulas and none of the three packages uses the formula that you give as the "actual" logFC.
People often make the mistake of posting on this support site results that they got from different packages and asking why the numbers are different. This overlooks the fact that the packages are not canned analyses but rather flexible pipelines with lots of options. So it isn't meaningful to give results from a package without explaining the details of how the package was used. So you have to give all the steps and options chosen (e.g., by posting detailed code) leading up to the results that are presented.
Your question boils down to whether the arithmetic or geometric means (of expression within each group) are better for computing log-fold changes. I dare say that this doesn't really have a clear answer - both means are valid measures of location, and the choice generally depends on what is mathematically convenient for your particular model. limma operates on the log-expression values, so it makes sense to set up the null hypothesis on the mean log-value (i.e., the log-geometric mean). edgeR and DESeq2 operate on the raw counts, so it is more convenient to set up the null hypothesis on the mean count (i.e., the arithmetic mean).
In general, I would not be particularly worried about the discrepancies between the two definitions. These should only occur when the variances are large and different between groups (a second-order Taylor approximation shows that the difference between E(log(x)) and log(E(x)) is a term dependent on the variance). If the differences between groups are strong enough to overcome the large variances, the discrepancy should be minor compared to the magnitude of the log-FCs and should not change your conclusions. If your DE is weak, then the affected gene wouldn't show up as significant in the first place, and so discrepancies can be ignored.
I guess that, with very large sample sizes, you can get significant hits for small log-fold changes that change in sign depending on whether you use one log-FC definition or the other. In such cases, you'll just have to keep in mind that each model is doing exactly what it's told; it's entirely possible for the arithmetic mean to increase in one group compared to another, while the geometric mean decreases. I suspect that this would require violation of the distributional assumptions, though, or at least some pathological scenarios.