Symmetry of a quantum Hamiltonian. Consider the quantum Heisenberg model:
$$H=-J\sum_{\left< \vec r,\vec r'\right>} \hat S_\vec r\cdot \hat S_{\vec r'}\tag{1}$$
according to David Bar Moshe's answer on a related question this is symmetric under the parity operator:
$$P\circ \sigma_i^a=-\sigma_i^a\tag{2}$$
which manifestly it is in the fact that it does not change the form of (1). But my problem is that it does change the form of the commutation relations of the spins. I.e instead of:
$$ [ S_i, S_j] = i \varepsilon_{ijk} S_k \tag{3}$$
we get
$$ [ S_i, S_j] = -i \varepsilon_{ijk} S_k\tag{4}$$
so can we truly say that parity is a symmetry of the system - do we not care that the spin operators change commutation relation?
I would also be interested to know how (2) acts on the basis vectors of the underlying Hilbert space.
 A: A quantum symmetry does not necessarily preserve the system's operator algebra commutation relations. Let me please refer you to Klaas Landsman's book: Foundations of quantum theory (There is a downloadable version  at researchgate) - chapter 9 page 334, definition 9.2 (Wigner symmetry)
This definition can be rephrased as follows: 
A Wigner symmetry is a continuous bijection of the pure state space which preserves transition probabilities. 
(Landsman gives a total of 6 (almost equivalent) definitions of a symmetry, each emphasizes another aspect of quantum theory).
I'll elaborate the case of the parity symmetry of a spin $ \frac{1}{2}$ system
A pure state of a pin-$\frac{1}{2}$ system can be represented by a density matrix which is also a projector:
$$\rho = \frac{1}{2}(1 + \sum_i x^i \sigma_i)$$
With: $\sum_i x_i^2 = 1$. The last condition ensures that the density matrix is a projector: $\det \rho = 0$.
Since the parity operator reverses the sign of the Pauli matrices $\sigma_i \rightarrow {\sigma}'_i = -\sigma_i$, it acts on the pure states as:
$$\rho \rightarrow   \rho' =  \frac{1}{2}(1 - \sum_i x^i \sigma_i)$$
preserving expectations:
$$\mathrm{tr}( \rho \sigma_i) = \mathrm{tr} ( \rho' {\sigma_i}').$$
Given two pure states $\rho_x = \frac{1}{2}(1 + \sum_i x^i \sigma_i)$ and $\rho_y = \frac{1}{2}(1 + \sum_i y^i \sigma_i)$. Their transition probability can be written as:
$$\tau(\rho_x, \rho_y) = \mathrm{tr}( \rho_x\rho_y)$$
In our case we have 
$$\tau({\rho}'_x, {\rho}'_y) = \mathrm{tr}( \frac{1}{2}(1 - \sum_i x^i \sigma_i) \times \frac{1}{2}(1 - \sum_i y^i \sigma_i) )= \frac{1}{2}(1+\sum_ix^iy^i) = \mathrm{tr}( \rho_x\rho_y) =\tau(\rho_x, \rho_y) $$ 
Thus the parity is a symmetry.
