Will quadratic-based estimation (not necessarily MSE) always generate a symmetric residuals after training it? These are error's empirical distribution for XGB, RF and kNN, the last one have taken on another dataset.

Neither of them is normally distribuited but they all are symmetric. None of used algorithms have even used MSE optimum, for example, both XGB and RF make a greedy approach for it due to being decision-tree-based and kNN uses euclidean distance which has nothing to do with MSE because that is not even a error-based estimation, my guess is that happens due to quadratic-based methods ignore error's signal but I can't link that to symmetry in probability density sense.

NO
For a counterexample, fit an OLS regression with an exponential error term.
$$
Y=X\beta+\epsilon\\
\epsilon_i\overset{iid}{\sim}exp(1)
$$
When you do the OLS fit, which uses square loss, you will get asymmetric residuals.
(There’s an annoying issue where $exp(1)$ has a mean of $1$ instead of $0$. You can resolve this if you do a simulation in software like R by subtracting $1$ from every error value. That isn’t an exponential distribution but a shifted exponential, but the asymmetric residuals will be present.)
Here is a simulation demonstrating my counterexample.
set.seed(2021)
N <- 1000
x <- seq(0, 10, 10/(N - 1))
e <- rexp(N, 1) - 1 # subtract 1 so E(e) = 0
y <- x + e
plot(x, y)
L <- lm(y ~ x)
hist(resid(L))
plot(L)

The histogram of the residuals and first two plots of the model L show the marked asymmetry of the residuals.



