All right, let's work out _the_ quantum toy model: a good old spin $1/2$...

From Hilbert spaces to the algebraic formalism...
=================================================

The Hilbert space formulation of a spin $1/2$ is as follow:

 1. the Hilbert space is $\mathcal{H} := \mathbb{C}^2$

 2. the elementary observables are the [Pauli matrices](https://en.wikipedia.org/wiki/Pauli_matrices) $\hat{\sigma}_1, \hat{\sigma}_2, \hat{\sigma}_3$, which satisfy $$\hat{\sigma}_i \hat{\sigma}_j = \delta_{ij} \mathbf{1} + i \varepsilon_{ijk} \hat{\sigma}_k$$

From there we can build the _abstract_ *-algebra $\mathcal{A}$ generated by 3 elements $\sigma_1,\sigma_2,\sigma_3$ satisfying $\sigma_i \sigma_j = \delta_{ij} 1 + i \varepsilon_{ijk} \sigma_k$ and $\sigma _i^* = \sigma _i$. Turns out it is the 4-dimensional algbra $M_2(\mathbb{C})$ mentioned by @Phoenix87, a basis of which is $$1,\sigma_1,\sigma_2,\sigma_3$$

For any $a = \lambda 1 + \mu_i \sigma_i \in \mathcal{A}$, we define $$\hat{a} := \lambda \mathbf{1} + \mu_i \hat{\sigma}_i$$ Note that $\forall a,b \in \mathcal{A}, \widehat{ab} = \hat{a} \hat{b}$ by construction, and that $\left\{ \hat{a} \middle| a \in \mathcal{A} \right\}$ spans the full space of operators on $\mathcal{H}$ (namely, the space of $2\times 2$ complex matrices, which constitutes the defining representation of the abstract algebra $M_2(\mathbb{C})$).

Now, given a density matrix $\hat{\rho}$ on $\mathcal{H}$, we define a linear form $\rho: \mathcal{A} \rightarrow \mathbb{C}$ by $$\forall a \in \mathcal{A}, \rho(a) := \text{Tr} \hat{\rho} \hat{a}$$ We can check that $\rho$ is an algebraic state on $\mathcal{A}$ ie:

 1. $\rho(1) = 1$ ($\text{Tr} \hat{\rho} = 1$)

 2. $\forall a \in \mathcal{A}, \rho(a^*) = {\rho(a)}^*$ (for $\widehat{a^*} = \hat{a}^{\dagger}$)

 3. $\forall a \in \mathcal{A}, \rho(a^* a) \geq 0$ (we can write $\hat{\rho} = \sum_{k} p_{k} \left| \psi_{k} \middle\rangle\middle\langle \psi_{k} \right|$, so $\text{Tr} \hat{\rho} \hat{a}^{\dagger} \hat{a} = \sum_{k} p_{k} \left\langle \hat{a} \psi_{k} \middle| \hat{a} \psi_{k} \right\rangle \geq 0$)

So we can associate to each density matrix an algebraic state on $\mathcal{A}$. Moreover the map $\hat{\rho} \mapsto \rho$ has the following properties:

 1. It is _injective_: $\rho$ completely encodes the quantum state $\hat{\rho}$. Mathematically, this holds because for any $\left|e\right\rangle, \left|f\right\rangle \in \mathcal{H}$, there exists $a \in \mathcal{A}$ such that $\hat{a} = \left|f\middle\rangle\middle\langle e\right|$ (as stressed above, _any_ operator on $\mathcal{H}$ can be written as $\hat{a}$ for some $a \in \mathcal{A}$), so $\left\langle e \middle| \hat{\rho} \middle| f \right\rangle$ can be read out from $\text{Tr} \hat{\rho} \hat{a} = \rho(a)$. Physically, all we use $\hat{\rho}$ for is to compute expectation values, and these are precisely what the linear form $\rho$ supplies: to each observable $a$, it maps its expectation value.

 2. It is also _surjective_: this can be seen by noting that any algebraic state $\tau$, being a linear form over $\mathcal{A}$, is entirely characterized by the 4 complex numbers $\tau(1), \tau(\sigma_1), \tau(\sigma_2), \tau(\sigma_3)$. Imposing normalization (condition 1. above), we get $\tau(1) = 1$. Condition 2. implies that the numbers $\tau_i := \tau(\sigma_i)$ are actually _real_. Finally, the positivity of $\tau$ (condition 3.) requires that, for any $a = \lambda\,1 + \mu_i \sigma_i \in \mathcal{A}$,
$$\tau(a^* a) = \lambda_R^2 + \lambda_I^2 + \left\|\vec{\mu}_R\right\|^2 + \left\|\vec{\mu}_I\right\|^2 + 2 (\lambda_R \vec{\mu}_R + \lambda_I \vec{\mu}_I + \vec{\mu}_I \times \vec{\mu}_R) . \vec{\tau} \geq 0$$
(where the subscripts ${}_R$ resp. ${}_I$ denote real resp. imaginary part).
Minimizing over $\lambda_R$ and $\lambda_I$, this is equivalent to the condition
$$\forall \vec{u}, \vec{v} \in \mathbb{R}^3, (\vec{u}.\vec{\tau})^2 + (\vec{v}.\vec{\tau})^2 + 2 (\vec{u}\times\vec{v}).\vec{\tau} \leq \left\|\vec{u}\right\|^2 + \left\|\vec{v}\right\|^2$$
which is fulfilled iff $\vec{\tau}$ belongs to the _unit ball_ of $\mathbb{R}^3$. We can then check that $$\hat{\tau} := \frac{1}{2} \left[ \mathbf{1} + \tau_i \hat{\sigma}_i \right]$$ is a well-defined density matrix on $\mathcal{H}$, which maps to the algebraic state $\tau$ (since $\text{Tr} \mathbf{1} = 2$ and $\text{Tr} \hat{\sigma}_i = 0$).

Thus, the space of density matrices over $\mathcal{H}$ can be identified _one-to-one_ with the space of algebraic states over $\mathcal{A}$, so we can do physics without ever referring to $\mathcal{H}$: knowing $\mathcal{A}$, we can _parametrize_ the space of states of our quantum theory by algebraic states, and, given one such algebraic state, we can compute the corresponding expectation values for measurements.

...and back (aka. GNS construction)
===================================

Reciprocally, if we start from a quantum theory formulated the algebraic way, the GNS construction allows to recover an Hilbert space formulation for it. This is important as this guarantees that, by using the algebraic formalism, we are not going to start constructing theories that have no counterpart in the Hilbert space formulation, ie. theories that wouldn't be admissible quantum theories in the usual sense.

For the algebra $\mathcal{A} = M_2(\mathbb{C})$, this works as follow. We need to first choose a state $\omega$ (the "vacuum"), characterized by $\vec{\omega}$ in the unit ball of $\mathbb{R}^3$ as described above. Then, we compute its "kernel" (I'm not completely sure it's called like that, but can't find a better term):
$$\mathcal{K}_{\omega} := \left\{ a \in \mathcal{A} \middle| \omega(a^* a)=0 \right\}$$
Writing $a = \lambda\,1 + \mu_i \sigma_i \in \mathcal{A}$, we have $a \in \mathcal{K}_{\omega}$ iff $$\lambda_R^2 + \lambda_I^2 + \left\|\vec{\mu}_R\right\|^2 + \left\|\vec{\mu}_I\right\|^2 + 2 (\lambda_R \vec{\mu}_R + \lambda_I \vec{\mu}_I + \vec{\mu}_I \times \vec{\mu}_R) . \vec{\omega} = 0$$
Since we know that $\lambda_{R/I}^2 + 2 \lambda_{R/I} (\vec{\mu}_{R/I} . \vec{\omega}) + (\vec{\mu}_{R/I} . \vec{\omega})^2 \geq 0$ and $\left\|\vec{\mu}_R\right\|^2 + \left\|\vec{\mu}_I\right\|^2 - (\vec{\mu}_R .\vec{\omega})^2 - (\vec{\mu}_I .\vec{\omega})^2 - 2 (\vec{\mu}_R \times\vec{\mu}_I).\vec{\omega} \geq 0$ (see above), this becomes
$$a \in \mathcal{K}_{\omega} \;\Leftrightarrow\; \lambda_{R/I} = - \vec{\mu}_{R/I} . \vec{\omega} \;\;\&\;\; (\vec{\mu}_R .\vec{\omega})^2 + (\vec{\mu}_I .\vec{\omega})^2 + 2 (\vec{\mu}_R \times\vec{\mu}_I).\vec{\omega} = \left\|\vec{\mu}_R\right\|^2 + \left\|\vec{\mu}_I\right\|^2$$

At this point, we need to distinguish cases. In the case where $\vec{\omega}$ is actually on the unit _sphere_ (ie $\|\vec{\omega}\| = 1$), we find
$$\mathcal{K}_{\omega} = \text{Span} \left\{ 1 - \omega_i \sigma_i, (u_j + i v_j) \sigma_j \right\}$$
where $\vec{u}, \vec{v}$ are some unit vectors such that $\vec{\omega}, \vec{u}, \vec{v}$ form an _orthonormal_ basis.
So, in this case, $\mathcal{K}_{\omega}$ is 2-dimensional, and therefore $\mathcal{H}_{\omega} := \mathcal{A} / \mathcal{K}_{\omega}$ is $4-2=2$-dimensional.
$\mathcal{H}_{\omega}$ is spanned eg. by $\left [1\right ]$ and $\left [u_j \sigma_j\right ]$ (where $\left [\,\cdot\,\right ]$ denotes equivalence classes modulo $\mathcal{K}_{\omega}$).

On $\mathcal{H}_{\omega}$ we can define a scalar product $$\left\langle \left[a\right] \middle| \left[b\right] \right\rangle_{\omega} := \omega(a^*b)$$ and a representation $a \mapsto \hat{a}^{\omega}$ of the algebra $\mathcal{A}$ $$\hat{a}^{\omega} \left[b\right] := \left[ab\right]$$ (check that these are well-defined, ie that the result of each map does not depend on the representatives used to express its arguments, and that each map has the required properties to be a scalar product, resp. a representation of $\mathcal{A}$).

Finally, we want to confirm that the thus obtained Hilbert space/representation reproduces the $\mathcal{H}$ we started from at the very beginning. Indeed, denoting by $\left|\phi_{\omega}\right\rangle$ the eigenvector of $\omega_i \hat{\sigma}_i$ with eigenvalue $+1$ (ie spin $+1$ in the direction of $\vec{\omega}$), $\mathcal{H}_{\omega}$ can be identified with $\mathcal{H}$, via the map $$\begin{array}{ccc}\mathcal{H}_{\omega} & \rightarrow & \mathcal{H}\\ \left[a\right] & \mapsto & \hat{a} \left|\phi_{\omega}\right\rangle \end{array}$$ which is well-defined, bijective and matches the scalar products and representations on both sides.

This can be checked directly, or by observing that $\omega(a) = \left\langle \phi_{\omega} \middle| \hat{a} \middle| \phi_{\omega} \right\rangle$ (as the algebraic state $\omega$ corresponds to the density matrix $\left|\phi_{\omega}\middle\rangle\middle\langle \phi_{\omega}\right|$ via the identification discussed in the previous section). For example, this implies that $\mathcal{K}_\omega$ consists of the $a \in \mathcal{A}$ such that $\left\langle \hat{a} \phi_{\omega} \middle| \hat{a} \phi_{\omega} \right\rangle = 0$, ie $\hat{a} \left|\phi_{\omega}\right\rangle = 0$, and therefore $\left[a\right] = \left[b\right] \Leftrightarrow a-b \in \mathcal{K}_{\omega} \Leftrightarrow \hat{a} \left|\phi_{\omega}\right\rangle = \hat{b} \left|\phi_{\omega}\right\rangle$.

In other words, we have been able to _entirely reconstruct_ $\mathcal{H}$ just by knowing the algebra $\mathcal{A}$ (choosing some algebraic state $\omega$ on $\mathcal{A}$ to "seed" $\mathcal{H}_{\omega}$). Crucially, this reconstruction didn't assume any _a priori_ knowledge of $\mathcal{H}$ (the definitions of $\mathcal{K}_{\omega}$, $\mathcal{H}_{\omega}$, $\left\langle\cdot\middle|\cdot\right\rangle_{\omega}$ and $a \mapsto \hat{a}^{\omega}$ were purely in terms of 1/ the operations in $\mathcal{A}$, ie addition, multiplication and adjoint, and 2/ the state $\omega$).

What about mixed states?
========================

But what would have happened if we had chosen a state characterized by a $\vec{\omega}$ in the _interior_ of the unit ball (ie $\|\vec{\omega}\| < 1$)? In this case, $\mathcal{K}_{\omega}$ would be just $\{0\}$ so $\mathcal{H}_{\omega}$ would be 4-dimensional. One can check that the representation created by the GNS construction can then be identified with $\mathcal{H} \oplus \mathcal{H}$ via:
$$\begin{array}{ccc} \mathcal{H}_{\omega } \approx  \mathcal{A} & \rightarrow  & \mathcal{H}^{(1)} \oplus  \mathcal{H}^{(2)}\\ \left[a\right] = a & \mapsto  & \hat{a} \left( \cos \eta  \left|\phi ^+_{\omega }\right\rangle ^{(1)} + \sin \eta  \left|\phi ^-_{\omega }\right\rangle ^{(2)} \right) \end{array}$$
where $\eta  \in  \left[0,\frac{\pi }{2}\right[$ is defined by $\cos^2 \eta  := \frac{1 + \|\omega \|}{2}$ and $\left|\phi ^{\pm }_{\omega }\right\rangle$ denote the eigenvectors of $\frac{\omega _i}{\|\omega \|} \hat{\sigma }_i$ with eigenvalue $\pm 1$ respectively.

It is a general property of the GNS construction that the representation arising from a _mixed state_ will typically be a _direct sum_ of _independent_ representations (aka. [superselection sectors](http://en.wikipedia.org/wiki/Superselection)), while _pure states_ lead to irreducible representations (ie. representations that _cannot_ be split as a direct sum of simpler representations).

Since density matrices are now $4\times 4$ matrices, does it means that we suddenly have more quantum states that we had on just $\mathcal{H}$? Well, actually _no_. Density matrices on $\mathcal{H}^{(1)} \oplus  \mathcal{H}^{(2)}$ can be written as block matrices: $$\hat{\rho }^{\omega } = \left( \begin{array}{cc} \hat{\rho }^{(1)} & \hat{\gamma } \\ \hat{\gamma }^{\dagger} & \hat{\rho }^{(2)} \end{array} \right)$$ with $\hat{\rho }^{(1)}, \hat{\rho }^{(2)}$ positive semi-definite and $\text{Tr} \hat{\rho }^{(1)} + \text{Tr} \hat{\rho }^{(2)} = 1$ (as well as suitable bounds on $\hat{\gamma }$ for $\hat{\rho }^{\omega }$ to be positive semi-definite). But since observables have the form: $$\hat{a}^{\omega } = \left( \begin{array}{cc} \hat{a} & 0 \\ 0 & \hat{a} \end{array} \right)$$ expectation values (which are the only physically measurable quantities) are simply given by $$\text{Tr} \hat{\rho }^{\omega } \hat{a}^{\omega } = \text{Tr} \left[ \left(\hat{\rho }^{(1)} + \hat{\rho }^{(2)}\right) \hat{a} \right]$$
So we don't have any new quantum state, we just have many different density matrices corresponding to the _same_ quantum state, namely the state that on $\mathcal{H}$ was represented by the density matrix $\hat{\rho } = \hat{\rho }^{(1)} + \hat{\rho }^{(2)}$. In this sense, $\mathcal{H}_{\omega }$ describes just the _same_ quantum theory as $\mathcal{H}$, but with some extra _redondancy_ in the representation of quantum states.

A word of caution: foliums
==========================

To be honest, this toy model works a bit _too well_:

 1. all algebraic states can be represented as density matrices on $\mathcal{H}$

 2. all algebraic states give rise, via the GNS construction, to one or more copies of $\mathcal{H}$

Propositions 1. and 2. are in fact equivalent: they express the fact that there is only one _"folium"_ of states on $M_2(\mathbb{C})$. The folium of an algebraic state $\omega$ is defined as the space of algebraic states that can be represented as density matrices of the GNS representation $\mathcal{H}_{\omega }$ arising from $\omega$. If $\omega$ is pure, all states in its folium will have one or more copies of $\mathcal{H}_{\omega }$ as their GNS representation. But in general an algebra of observables can admit many disjoint foliums, or, equivalently, many _inequivalent_ representations.

For example, if we were to redo the analysis for the *-algebra generated again by 3 self-adjoint elements $\sigma_1,\sigma_2,\sigma_3$ but now only assuming the weaker relations $\left[\sigma_i, \sigma_j\right] = 2i \varepsilon_{ijk} \sigma_k$, aka the [universal enveloping algebra](https://en.wikipedia.org/wiki/Universal_enveloping_algebra) of the Lie algebra $\mathfrak{so}(3) = \mathfrak{su}(2)$, we would find infinitely many inequivalent representations, corresponding to all possible spins $j$.

In the case of ordinary quantum mechanics with an algebra of observables generated by _exponentiated_ position and momentum operators (to avoid issues with unbounded operators), the [Stone–von Neumann theorem](http://en.wikipedia.org/wiki/Stone%E2%80%93von_Neumann_theorem) implies that there is only one folium of _"sufficiently regular"_ states, although if we drop the regularity condition more exotic states and associated GNS representations can be constructed. In QFT, on the other hand, there is a plethora of distinct foliums. Moreover, on _curved spacetimes_, there may be no way to single out a _preferred_ vacuum state/folium and associated representation: that's where the algebraic formalism truly shines, because it allows to discuss all foliums at once, without having to specialize to an arbitrarily chosen one as we would have to do to get an Hilbert space formulation.