I think it might be conceptually easier to start backwards. Let us assume we have $N$ decoupled harmonic oscillators:
$$
H = \sum_{a=1}^N \left[ \frac{(\pi_a)^2}{2m} + \frac{1}{2} m \omega_i^2 \phi_a^2 \right]
$$
Now imagine applying the unitary transformation $x_j = {U_j}^a \phi_a$ and $p_j = {U_j}^a \pi_a$ where $U^T U = 1$ in order that both $[x,p]$ and $[\phi,\pi]$ have the standard commutation relation. Such a transformation will act simply on the $\pi_i^2$ term, but the $\omega_i$'s in the second contribution to $H$ will lead to interesting quadratic couplings between the $x_j$. (If you further want unequal masses (I'm not sure you really do in the context of entanglement entropy of a free scalar field), you can perform the scaling $x_i \to \sqrt{m_i} x_i$ and $p_i \to p_i / \sqrt{m_i}$ without changing the canonical commutation relation between the $x$'s and $p$'s.)
Now the question is about a very particular nearest neighbor coupling between the $x_i$. The usual way to think about this transformation involves Fourier series. If I'm not mistaken, something like
$$
x_j = \sum_{a=1}^N e^{2 \pi i j a /N} \phi_a
$$
should act to diagonalize your Hamiltonian if I further make the assumption that $x_{N+1} = x_1$, i.e. I'm living on a circle.
In the context of the Srednicki paper and entanglement entropy, there are some more modern approaches involving two-point functions. You might want to look at section 2.2 of the review http://arxiv.org/abs/0905.2562 or section III of http://arxiv.org/abs/0906.1663.
Both hamiltonian pieces operate on both mutually equivalent bases--it's just that their action on each is different, as in "non diagonal".
Given both diagonal pieces in some basis, however, you may easily reconstruct the non-diagonal pieces you skipped in the respective bases. No further work is required. (I suppose the time hiatus ensures I am not
doing your homework for you.)
Your coupled basis vectors (call them v), are the eigenstates of $S$ and $S_z$, whence $H_1$ as well. While your uncoupled-basis vectors, w, consist of the eigenstates of $s_{1z}\otimes s_{2z}$, whence $H_2$ as well,
$$
v=\begin{bmatrix}
\frac{\uparrow \downarrow -\downarrow \uparrow}{\sqrt 2} \\
\uparrow \uparrow \\
\frac{\uparrow \downarrow +\downarrow \uparrow}{\sqrt 2} \\
\downarrow \downarrow \\
\end{bmatrix} ,\qquad
w=\begin{bmatrix}
\uparrow \uparrow \\
\uparrow \downarrow \\
\downarrow \uparrow \\
\downarrow \downarrow
\end{bmatrix}.
$$
The two bases are, of course, interconvertible through the orthogonal Clebsch-Gordan matrix, which I trust your instructor has introduced you to,
$$
Uw=v, \qquad \begin{bmatrix}
0& \frac{1}{\sqrt 2} & \frac{-1}{\sqrt 2} &0 \\
1&0& 0 &0 \\
0& \frac{1}{\sqrt 2} & \frac{1}{\sqrt 2} &0 \\
0 &0 &0 &1
\end{bmatrix}, \qquad U^T U=\mathbb{1} ~.
$$
So, then, setting $\hbar=1$ and $\mu_B B=1$ as you all but did, you may write your Hermitean hamiltonian in either basis, always with only a piece of it diagonal, but never both, as you implicitly observe:
$$ H_v= \frac{D}{4} \begin{bmatrix}
-3& 0 & 0 &0 \\
0& 1 & 0 &0 \\
0& 0 & 1 &0 \\
0 &0 &0 & 1
\end{bmatrix}
+\frac{1}{2}\begin{bmatrix}
0&0& g_{1}-g_{2} & 0 \\
0& g_{1}+g_{2} & 0 &0 \\
g_{1}-g_{2} &0&0&0 \\
0 &0 &0 & -g_{1}-g_{2}
\end{bmatrix},$$
$$H_w=U^T H_v U=\frac{D}{4} \begin{bmatrix}
1 & 0 & 0 &0 \\
0& -1 & 2 &0 \\
0& 2 & -1 &0 \\
0 &0 &0 & 1
\end{bmatrix} +\frac{1}{2}\begin{bmatrix}
g_{1}+g_{2} & 0 & 0 &0 \\
0& g_{1}-g_{2} & 0 &0 \\
0& 0 & -g_{1}+g_{2} &0 \\
0 &0 &0 & -g_{1}-g_{2}
\end{bmatrix}.$$
Both Hermitean, of course, and only diagonal in their leading and trailing pieces, respectively, as you evidently encountered. $H_{2v}$, naturally, was found from the $UH_{2w}U^T$ you computed.
In the coupled basis, you see the non-diagonal $H_2$ leaves the $S_z=\pm 1$ components alone, up to multiplication by these eigenvalues. But permutes the $S=1$ with the $S=0$ components that share a common $S_z=0$.
Analogously, in the uncoupled basis $H_1$, again the $|\uparrow \uparrow\rangle$ and $|\downarrow \downarrow\rangle$ components are left alone, but the singlet and triplet $S_z=0$ components are mixed thoroughly.
Best Answer
The Hilbert space of the system is the direct product of the Hilbert spaces corresponding to separate oscillators. Operators for every specific oscillator are transformations in its own Hilbert space, which do not affect other Hilbert spaces.
However, the above is true for the operators taken at the same time instant. If we now work in Heisenberg picture, write Heisenberg equations of motion for operators and solve them, the operators at time $t$ will be linear combinations of the operators at some previous time $t'$ and vice versa. These operators taken at different time moments will not commute, generally speaking, even for the same oscilator.