  In this paper we show that standard implementations of fluctuating Lattice
Boltzmann methods do not obey Galilean invariance at a fundamental level. In
trying to remedy this we are led to a novel kind of multi-relaxation time
lattice Boltzmann methods where the collision matrix depends on the local
velocity. This new method is conceptually elegant but numerically inefficient.
With a small numerical trick, however, this method recovers nearly the original
efficiency and allows the practical implementation of fluctuating lattice
Boltzmann methods with significantly improved Galilean invariance. This will be
important for applications of fluctuating lattice Boltzmann for non-equilibrium
systems involving strong flow fields.
