Zero magnetization of spin model without external magnetic field For a given Hamiltonian with spin interaction, say Ising model
$$H=-J\sum_{i,j} s_i s_j$$
in which there are no external magnetic field. The Hamiltonian is invariant under transformation $s_i \rightarrow -s_i$, so there are always two spin states with exactly same energy.
For the magnetization $M = \sum_i s_i$, we can take the ensemble average
$$\left\langle M\right\rangle =\sum_{\{s_{i}\}}M\exp(-\beta E)$$
and the result should be $\left\langle M\right\rangle  = 0$ because the two states with all spin flipped will exactly cancel each other. There is an argument for this in the wiki.
So the question: how is this situation handled for finite and infinite lattice? How can they obtain the non-zero magnetization for the 2D Ising model?
$$M=\left(1-\left[\sinh\left(\log(1+\sqrt{2})\frac{T_{c}}{T}\right)\right]^{-4}\right)^{1/8}$$

Some information: for a 1D Ising model with external magnetic field, one can solve the Hamiltonian
$$H=-J\sum_{i,j} s_i s_j - \mu B \sum_i s_i$$
and obtains the magnetization as
$$\left\langle M\right\rangle =\frac{N\mu\sinh(\beta\mu B)}{\left[\exp(-4\beta J)+\sinh^{2}(\beta\mu B)\right]^{1/2}}$$
It gives the result $\left\langle M\right\rangle  \rightarrow 0$ when $B \rightarrow 0$ for any temperature and this match the definition of the magnetization above. However, it gives $\left\langle M\right\rangle  \rightarrow N\mu$ when we take the limit of $T \rightarrow 0$. It suggests the ordering of limit is important, but we still get $\left\langle M\right\rangle  = 0$ when there is no external magnetic field.

Reminder: There are some precautions needed to care when you run computer simulation using the definition of magnetization above directly, otherwise, you will always get 0. These methods are similar to create a spontaneous symmetry breaking manually, for the Ising model, the following may be used:


*

*Use the $\left\langle |M|\right\rangle$ instead

*Fix the state of one spin so that the system will close to either $M=+1$ or $M=-1$ at low temperature.


In general, the finite size scaling should be used because we are likely interested in the thermodynamic limit of the system.

Update:
Visualization should explain the problem better. Here is the figure of the canonical distribution as a function of magnetization $M$ and temperature $T$ for 3D Ising model ($L = 10$).

At a fixed temperature $T \lesssim T_c$, there are two symmetric peaks with opposite magnetization. If we blindly use the $\left\langle M\right\rangle$ defined before, we will get $\left\langle M\right\rangle = 0$. So how do we deal with this situation?
For a finite system, there is a finite probability that the transition between two peaks can occur. However, the "valley" between two peaks will become deeper and deeper when the size of the system $L$ increase. When $L \rightarrow \infty$, the transition probability tends to zero and those two configuration space should be separated. Note that the "flat mountain" in the figure at $L = 10$ will also become a very very sharp peak when $L \rightarrow \infty$.
One method discussed in the answers below is to consider the average of $M$ for separate configuration spaces. This seems reasonable for infinite system, but becomes a problem for finite system. Another problem raised here is that how to find each separated configuration spaces?

Thanks people try to give the answers to this question. In the following discussion, Kostya gives the typical treatment of spontaneous symmetry breaking. Marek discusses the ensemble average below and above the critical temperature. Greg Graviton gives an analogue for the real space spontaneous symmetry breaking.
If anyone can explain better to the problem of configuration space and how to take the average for Ising model, other spin models or general case, you are welcome to leave answer here.
 A: If I understand you correctly, you seem to be asking the fundamental question of how on earth broken symmetries are possible. After all, there is an apparent contradiction: on one hand, the thermal average $\langle M \rangle= \sum_{\{s_i\}} M \exp(-\beta E)$ must vanish due to symmetry, while on the other hand, physicists claim that the 2D Ising model gives rise to a nonzero magnetization. How can that be?
The following thought experiment will provide an instructive answer:
Imagine a table with an upside-down pendulum attached to it. We assume that the head of the pendulum can only move from left to right or the other way round until it hits the table, i.e. it is confined to one dimension and its position can be described by a single angle $\phi\in[-\pi,\pi]$.

The potential energy of the pendulum is given by $E(\phi)=-A\cos \phi$. Clearly, the whole situation is symmetric with respect to a reflection $\phi \mapsto -\phi$; the potential energy is invariant under this reflection.
But this symmetry of the energy does not imply that the stable states are symmetric themselves. Namely, the energy has two, asymmetric minima at $\phi=\pi$ and $\phi=-\pi$ respectively (the pendulum is stopped by the table); the stable states are no longer symmetric themselves. Of course, the symmetry is not lost; if you have calculated one stable state, you can apply the reflection and get the other, but the point is that there are several stable states now. The energy is symmetric, but its minima are not.

Now, imagine the following experiment: the table with the pendulum is located in one room and you are located in a different room. A friend of yours puts the pendulum exactly in the $\phi=0$ position (which is unstable) and then lets go. Since you are in another room, you can't see what is going to happen next in the thermodynamic system consisting of the pendulum and your friend.
Question: before you go into the other room, what is the average position, the "ensemble average" of the pendulum? Well, it will fall to either the left or the right with equal probability, so after a short amount of time in which the pendulum hits the table and thermal equilibrium is restored, the average will be…
$$\langle \phi \rangle = \frac12 \pi + \frac12 (-\pi) = 0 \quad!$$
Of course, if you now run into the other room expecting the pendulum to be at position $\phi=0$, dangerously balancing on top of the table, you will be disappointed to learn that it has fallen to either the left $\phi=\pi$ or to the right $\phi=-\pi$.
What has happened? The point is that the thermal average is indeed zero, but the situation is so unstable that as you soon as you acquire a teensy bit of information, like a short glimpse into the other room, or knowledge of the position of just one single molecule of the pendulum, the situation will split in two and you will observe one of two different outcomes. The outcomes have equal probability, that's why you get an ensemble average of zero, either one could have happened. But they occupy quite different parts of the configuration space, and you suddenly get a nonzero angle or nonzero magnetization. Compare this to an ideal gas where getting information about the position of a single molecule has a vanishing impact on the situation at large. For the pendulum, knowing the position of just one molecule changes everything, and that's why the thermal average does not give an adequate description of what you will measure.

Update
To elaborate on the last paragraph:
It is important to understand what the thermal average $\langle \phi \rangle$ does. A very good explanation is provided in E.T. Jaynes' article The Evolution of Carnot's Principle. Namely, the thermal average is

the average over the probability distribution that maximizes the "number of ways in which the system can be realized" (= entropy), subject to the constraint that the average energy has a fixed value $\langle E \rangle = E_0$.

In other words, if we only know the average energy of the system and are ignorant about everything else, the maximum entropy principle gives a probability distribution with which we can calculate other quantities, like the average magnetization $\langle M \rangle$. For our pendulum, this reasoning makes sense: if we are in the other room, the best prediction about the position of the pendulum that we can make is indeed $\langle \phi \rangle = 0$.
If we go into the other room and look at the pendulum, we acquire slightly more information than just the average energy, for instance the position of a single molecule of the pendulum. The point is that the probability distribution is peaked around $\phi=\pi$ and $\phi=-\pi$, and the little information we acquire is enough to see that the system is located at one of the peaks. In other words, lifting our ignorance just a little bit gives rise to a marked change in prediction.
Compare this to an ideal gas, where learning of the position of a single molecule does not change anything about the probability distribution at large and won't change our prediction of macroscopic quantities.
A: Maybe I didn't get the question, but the whole point of the discussion of Ising model is not getting zero magnetization. This is a topic of the spontaneous symmetry breaking (further "SSB"). The symmetry you are talking about is broken spontaneously leading to nonzero magnetization. I must admit that I never got into details of Onsager solution, but I know that there are three approaches to this very important concept of modern physics.


*

*Limit of small "magnetic field".
You can consider a situation, when you got your system in a magnetic field $H$ first -- in that case there is no mentioned symmetry, and the magnetization gets along the magnetic field. Then you turn the magnetic field off.
The point is that the limits $H\to\infty$ and $N\to\infty$ do not commute, so if you first take the $N\to\infty$ limit and then $H\to\infty$ you will end up with spontaneous magnetization.
As far as I know -- this is the oldest approach to SSB and it is repeated in many statistical physics textbooks.  

*Boundary condition.
That so trivial that it is boring -- you just "put in" the preferred magnetization by setting appropriate boundary conditions for your system. Just fix the direction of your spins at the boundary or "at infinity".
As far as I know -- this is the approach adopted in mathematical physics analysis of SSB.

*No transition in your lifetime
That is the most "useless" and the most beautiful in my opinion. There is always a probability for your system to go to the state with opposite direction of magnetization due to some fluctuation with really small probability. So the magnetization is still zero if you average over infinite time. But those transitions are so improbable that you will definitely not see one even if you will wait for "billions of zillions" of years.
A: Let me ask a seemingly naive question first: what is actually a statistical ensemble? Mathematically, it's a prescription of an assignment of probability to every microstate. So this is a probability measure. But we have to be careful with measures defined on infinite systems. You certainly can't blindly work with $\exp(-\beta E)$ factors as energy will generally be infinite for an infinite system. But fortunately in the Ising model case (and lots of other spin models), it turns out that there exists a very nice Gibbs measure. I am not going to derive all the results on Gibbs measures here (see e.g. Georgii's bible) but basically one consistently prescribes how the measure should behave in all finite regions when subjected to some boundary conditions. The logic being that in finite regions everything behaves nicely (and energy is finite) and we can recover infinite-volume thermodynamic behavior by some limiting procedures.
Now, if you are talking about the ensemble average, you are supposing there is only one such ensemble and therefore only one Gibbs measure for a given model. It turns out that this assumption fails. There can exist more than one measure and this is actually equivalent with an appearance of a first-order phase transition. In other words, you have only one measure (or one phase) for $\beta \leq \beta_c$. For $\beta > \beta_c$ you will always have two (equally good) measures to choose from and actually you can also take their convex combination but it turns out only extreme points of the convex set of all Gibbs measures are what we would call a phase. The reason you obtained zero in your average is precisely because you've picked such a combination $1/2 (\mu^{\beta}_+ + \mu^{\beta}_-)$. But this is not a macroscopic phase but rather a description of your (complete) uncertainty that you don't yet know whether system is in $+$ state or the $-$ state.
The second-order phase transition (the one that corresponds to SSB here) is associated with "merging" of the Gibbs measures $$\lim_{\beta \to \beta_c} \mu^{\beta}_+ = \lim_{\beta \to \beta_c} \mu^{\beta}_- = \mu^{\beta_c}$$
Or perhaps more natural point of view is the "splitting" of the measure as you decrease the temperature, so that system has to "decide" which of the low-temperature phases it will pick.
Note: what I implicitily assumed above was that the measures are shift-invariant (homogeneous, if you will). If this condition is dropped there exist many more measures (associated with phases) which describe an interfaces. E.g. consider an Ising system in a square with $+$ boundary condition on the bottom of the box and $-$ on the top. In a thermodynamic limit one obtains a system which is mostly $+$ on the bottom, mostly $-$ on the top and with a fluctuating interface in the middle.
