Galilean invariance is a law of physics. If a numerical model is not galilean invariant, then it is non-physical.
Imagine a fluid simulation with periodic boundary conditions and no external force. The fluid has an imposed initial space-dependent velocity profile u(t=0,x) = u0(x), which could for example consist of a single vortex, or any other initial condition which is divergence-free and has an average velocity of 0. As the simulation goes on, the fluid develops a certain flow pattern and ultimately tends towards a rest state, because of energy dissipation. Let’s call this flow “FLOW 1”, and say that, by definition, it is simulated in a reference frame at rest.
Next, consider another flow, “FLOW 2”, with the same initial condition, except that the velocity is shifted in each point by a constant offset v: u(t=0,x) = u0(x)+v. What you expect from the time evolution of FLOW 2, is that it is the same as FLOW 1, with exception of a constant offset v which persists in time. In particular, the asymptotic solution at large time is not a rest state, but a flow in constant movement u(t=infinity,x) = v.
The reason why you know that this must be the case (physically speaking) is Galilean invariance. An observer moving at constant speed measures the same laws of physics as an observer at rest. Imagine an observer moving with constant speed v and looking at FLOW 2. This observer will have the impression to be at rest and to be observing FLOW 1. In conclusion, FLOW 2 has no choice but doing the same as FLOW 1, except for the speed offset.
If your numerical model yields a different solution for FLOW 2 than for FLOW 1, then it disrespects Galilean invariance. This means essentially that it has a preferred frame, or a “preferred velocity offset”. You can compare this to a numerical model that violates rotational invariance. Such a model has preferred directions, typically the lattice directions. The simulation results may be reasonably good if the main flow direction is aligned with this preferred direction, but quite bad along a diagonal direction.
When you use BGK with one of the common lattices (D2Q9, D3Q13, D3Q15, D3Q19, D3Q27), there is actually a preferred velocity: the model is not fully Galilean invariant. This reference frame is defined by a zero velocity. You can easily test this in any simulation of your choice. If you add a large velocity offset (say, 10 in lattice units) to the velocity, including to the velocity boundary conditions, the result is not simply shifted, as expected, but plain wrong. You may argue that this is the expected behaviour because the flow is in a supersonic regime (the Mach number is larger than 1), which is not supported by BGK. But this is an invalid argument. The Mach number is not defined from an absolute velocity, but rather from a velocity difference, for example the flow velocity as compared to the no-slip boundary of an obstacle.
All of this is usually not a problem, because simulations are always executed at small velocity, and I am not sure in which circumstances the lack of (exact) Galilean invariance becomes an issue. My guess would be than when you develop a LB model for high Mach number flows it becomes important to respect Galilean invariance with higher precision than in BGK.
As a matter of fact, you can recover full Galilean invariance by taking the continuum Boltzmann equation and discretizing systematically the space of velocity on a base of Hermite polynomials. This approach is for example developed in X. Shan e.a., 2006. If I understand this paper right, you need to go to a lattice with 39 neighbors (D3Q39) to obtain what you need. Now this may seem quite an expensive price to pay for a single law of physics. But on the other hand, consider that Galilei almost paid with his life to defend his point of view on physics.