When is a Hamiltonian Galilean invariant? The Hamiltonian of a point mass connected to a fixed point by a spring in (1-dimensional) space is
$$H(x,p)=\frac{p^2}{2m}+\frac{1}{2}kx^2\tag{1}\label{eq1}$$
The Hamiltonian equations are
$$\begin{equation}
\frac{dx}{dt}= \frac{\partial H}{\partial p} = \frac{p}{m}\\ 
\frac{dp}{dt}= -\frac{\partial H}{\partial x} = -kx
\tag{2}
\end{equation}$$
If $(x(t),\,p(t))$ is  a solution of this equation then $(x(t)+vt,\,p(t)+mv)$ isn't, so Hamiltonian $(\ref{eq1})$ isn't Galilean-invariant. This isn't surprising if we know the physical background: the "fixed point" is fixed only in one special reference frame. The situation immediately changes when we consider two point-masses connected to each other by a spring. For this system, the Hamiltonian is
$$H(x_1,x_2,p_1,p_2) = \frac{p_1^2}{2m_1}+\frac{p_2^2}{2m_2} + \frac{1}{2}k(x_2-x_1)^2\tag{3}\label{eq3}$$
, and the Hamiltonian equations are:
$$\begin{equation} 
\frac{dx_1}{dt}= \frac{\partial H}{\partial p_1} = \frac{p_1}{m_1}\\ 
\frac{dx_2}{dt}= \frac{\partial H}{\partial p_2} = \frac{p_2}{m_1}\\ 
\frac{dp_1}{dt}= -\frac{\partial H}{\partial x_1} = k(x_2-x_1)\\
\frac{dp_2}{dt}= -\frac{\partial H}{\partial x_2} = k(x_1-x_2)\tag{4}
\end{equation}$$
In this case, together with solution $s_0(t)=(x_1(t),x_2(t),p_1(t),p_2(t))$, $s_v(t) =(x_1(t)+vt,\,x_2(t)+vt,\,p_1(t)+m_1v,\,p_2(t)+m_2v)$ is also a solution, so Hamiltonian $(\ref{eq3})$ is Galilean-invariant as expected.
My question: what is the general, exact mathematical condition for a Hamiltonian $H$ of a system of point masses to be invariant for a Galilean boost?
 A: In order to try and answer your question I'm going to assume that your system has nothing in the way of Dirac's 1st, 2nd, etc.-kind constraints. This is to guarantee that we can shift from the Lagrangian to the Hamiltonian formalism smoothly.
I will also approach it in Cartesian coordinates. This is not because they are more fundamental, but because the action of the Galilean group is far more transparent in those.
The action of the Galilean group is,
$$
\boldsymbol{x}_{i}\mapsto A\boldsymbol{x}_{i}+\boldsymbol{V}t+\boldsymbol{b}
$$
where index $i$ tags the particle, $A$ is an orthogonal matrix (rotation):
$$
AA^{T}=A^{T}A=I
$$
$\boldsymbol{V}$ is a common velocity vector and $\boldsymbol{b}$ is a common translation vector (by "common" I mean neither depends on the particle index $i$).
Suppose your Hamiltonian admits a form,
$$
H=\sum_{i=1}^{n}\frac{\boldsymbol{p}_{i}^{2}}{2m_{i}}+U\left(\boldsymbol{x}_{1},\cdots\boldsymbol{x}_{n}\right)
$$
where,
$$
\boldsymbol{p}_{i}=m_{i}\dot{\boldsymbol{x}}_{i}
$$
The Galilean transformations induce these changes in the other dynamical/kinematical variables:
$$
\dot{\boldsymbol{x}}_{i}\mapsto A\dot{\boldsymbol{x}}_{i}+\boldsymbol{V}
$$
$$
\boldsymbol{p}_{i}=m_{i}\dot{\boldsymbol{x}}_{i}\mapsto A\boldsymbol{p}_{i}+m_{i}\boldsymbol{V}
$$
$$
\boldsymbol{p}_{i}^{2}=\mapsto A\boldsymbol{p}_{i}\cdot A\boldsymbol{p}_{i}+m_{i}^{2}\boldsymbol{V}^{2}+2m_{i}\boldsymbol{V}\cdot\boldsymbol{p}_{i}
$$
$$
A\boldsymbol{p}_{i}\cdot A\boldsymbol{p}_{i}=\boldsymbol{p}_{i}\cdot\boldsymbol{p}_{i}
$$
As for the Hamiltonian:
$$
H\mapsto\sum_{i=1}^{n}\frac{\boldsymbol{p}_{i}^{2}}{2m_{i}}+U\left(A\boldsymbol{x}_{i}+\boldsymbol{V}t+\boldsymbol{b},\cdots,A\boldsymbol{x}_{i}+\boldsymbol{V}t+\boldsymbol{b}\right)+\sum_{i=1}^{n}\frac{1}{2m_{i}}\left(m_{i}^{2}\boldsymbol{V}^{2}+2m_{i}\boldsymbol{V}\cdot\boldsymbol{p}_{i}\right)=
$$
$$
=H'+\sum_{i=1}^{n}\left(\frac{1}{2}m_{i}\boldsymbol{V}^{2}+\boldsymbol{V}\cdot\boldsymbol{p}_{i}\right)
$$
Now, the term $\sum_{i=1}^{n}\left(\frac{1}{2}m_{i}\boldsymbol{V}^{2}+\boldsymbol{V}\cdot\boldsymbol{p}_{i}\right)$ is a total derivative, and thus we can ignore it --going back to Lagrangian formalism.
We need to prove,
$$
H=H'
$$
and we're home-free. For this, it will suffice that,
$$
U\left(A\boldsymbol{x}_{i}+\boldsymbol{V}t+\boldsymbol{b},\cdots,A\boldsymbol{x}_{i}+\boldsymbol{V}t+\boldsymbol{b}\right)=U\left(\boldsymbol{x}_{i},\cdots,\boldsymbol{x}_{i}\right)
$$
It will be enough that $U$ is an $SO\left(3\right)$ scalar in the differences $\boldsymbol{x}_{i}-\boldsymbol{x}_{j}$.
That is:
$$
U\left(\boldsymbol{x}_{i},\cdots,\boldsymbol{x}_{i}\right)=\sum_{i,j=0}^{n}U_{ij}\left(\left\Vert \boldsymbol{x}_{i}-\boldsymbol{x}_{j}\right\Vert \right)
$$
I've proven you this as a sufficient condition. What about the "if and only if"? It's also true. You can find details in Landau-Lifshitz, Course of Theoretical Physics, Vol. 1 (Mechanics). But it's not difficult to convince yourself that this dependence on Cartesian coordinate differences in the potential energy (plus scalar character) is the key.
Your system perfectly illustrates all this, except that it only involves one spatial dimension, and the orthogonal group is $SO\left(1\right)$, instead of $SO\left(3\right)$.
A: Take the position vectors of the particles as
\begin{align*}
\boldsymbol R_i=\begin{bmatrix}
          x_i \\
          y_i \\
          z_i \\
        \end{bmatrix}\qquad 
\boldsymbol v_i=\boldsymbol{\dot{R}}_i =   \begin{bmatrix}
          \dot x_i \\
          \dot y_i \\
          \dot z_i \\
        \end{bmatrix}.
\end{align*}
It follows that the kinetic energy can written as
$$T=\frac{1}{2}\sum_{i=1}^{n_p} m_i\boldsymbol{{v}}^T_i\,\boldsymbol{{v}}_i,$$
and the potentail energy as
$$U=\frac{1}{2}\sum_{l=1}^{n_k}\,k_l\,U_l$$
where the spring potential energy  $~U_l~$ can be
$$U_l=U_l\left(\boldsymbol R_i^T\, \boldsymbol R_i\right)$$
or
$$U_l=U_l\left[(\boldsymbol R_j-\boldsymbol R_i)^T\,(\boldsymbol R_j-\boldsymbol R_i)\right]$$
The equations of motion
$$\boldsymbol M\boldsymbol{\ddot{q}}+\boldsymbol C\,\boldsymbol q=0\tag{1}$$
with
\begin{align*}
\boldsymbol M &=\frac{\partial }{\partial \boldsymbol{\dot{q}}} \left(\frac{\partial T}{\partial \boldsymbol{\dot{q}}}\right)\\
\boldsymbol C &=\frac{\partial }{\partial \boldsymbol{q}}\left(\frac{\partial U}{\partial \boldsymbol{q}}\right)\\
\boldsymbol q &=\left[x_1~,y_1~,z_1~,x_2~,y_2~,z_2~,
\ldots~,x_{n_p}~,y_{n_p}~,z_{n_p}\right]^T
\end{align*}
Galilean transformation
\begin{align*}
\boldsymbol R_i \mapsto \boldsymbol R_i+w\,t\,\underbrace{\begin{bmatrix}
                               1 \\
                               1 \\
                               1 \\
                             \end{bmatrix}}_{\boldsymbol e}\qquad 
  &\boldsymbol v_i\mapsto \boldsymbol v_i+w\,\boldsymbol e
\end{align*}
$$ T_G=\frac{1}{2}\sum_{i=1}^{n_p} m_i\left(\boldsymbol{{v}}_i+w\boldsymbol e\right)^T\,\left(\boldsymbol{{v}}_i+w\boldsymbol e\right)$$
$$U_G=\underbrace{\frac{1}{2}\sum_{l}\,k_l\left( (\boldsymbol R_i+w\,t \boldsymbol e )^T\,(\boldsymbol R_i+w\,t\boldsymbol e )\right)^T}_{U_t}+\frac 12\sum_{m}\,k_m\left[(\boldsymbol R_j-\boldsymbol R_i)^T\,(\boldsymbol R_j-\boldsymbol R_i)\right]$$
Therefore
$$\boldsymbol M_g=\boldsymbol M,$$
the equations of motion are Galilean invariant if
$$ \frac{\partial^2 U_G}{\partial t^2}=0~,\Rightarrow~ U_t=0 $$
the Hamiltonian is
$$\boldsymbol p^T\,\boldsymbol{\dot{q}}-\mathcal L(\boldsymbol q~,\boldsymbol p)$$
with
\begin{align*}
\mathcal{L} &=T_G-U_G\\
T &=\frac{1}{2}\,\boldsymbol{\dot{q}}^T\boldsymbol M\,\boldsymbol{\dot{q}}+\boldsymbol c\,\cdot \boldsymbol{\dot{q}}+b\\
\boldsymbol c &=w\,[m_1~,m_1~,m_1~,\ldots~,m_{n_p}~,m_{n_p}~,m_{n_p}]^T\\
b &=\sum_{i}^{n_p}\frac{3\,m_i}{2}\,w^2\\
\frac{\partial T}{\partial\boldsymbol{\dot{q}}} &=\boldsymbol M\,\boldsymbol{\dot{q}}+\boldsymbol c=\boldsymbol p
 ~\Rightarrow~,\boldsymbol{\dot{q}}=M^{-1}(\boldsymbol p-\boldsymbol c)
\end{align*}
