Why is state-space representation more numerically stable than zeros-poles representation? Matlab documentation says things like:

For high order filters, the state-space form is the most numerically accurate, followed by the zero-pole-gain form. The transfer function coefficient form is the least accurate; numerical problems can arise for filter orders as low as 15.

http://www.mathworks.com/help/signal/ref/besself.html

Conversions between the TF, ZPK, and SS representations involve numerical computations and can incur loss of accuracy when overused. Because the SS and FRD representations are best suited for numerical computations, it is good practice to convert all models to SS or FRD and only use the TF and ZPK representations for construction or display purposes.

http://www.mathworks.com/products/demos/shipping/control/GSModelTypeConversions.html?product=CT#9

You can represent numeric system components using any model type. However, Numeric LTI model types are not equally well-suited for numerical computations. In general, it is recommended that you work with state-space (ss) or frequency response data (frd) models, for the following reasons:
The accuracy of computations using high-order transfer functions (tf or zpk models) is sometimes poor, particularly for MIMO or high-order systems. Conversions to a transfer function representation can incur a loss of accuracy.

http://www.mathworks.com/help/control/ug/conversion-between-model-types.html#f3-1039600
I understand why zpk is better than tf, and know that ss can handle MIMO systems that the others can't, but is ss also numerically better than zpk?  Why?  Matlab generates filter prototypes in zpk and then converts them to ss for doing filter transformations, which seems odd since you can do the filter transformations directly in zpk without converting back and forth.
This is because polynomial root finding algorithms are terrible without any implication that people code it wrongly or else. It is an ill-posed problem whereas eigenvalue computations and hence generalized eigenvalue problem solutions are much more stable. 

Any mistake in any of the pole/zero data permeates to all entries of the polynomial entries and contaminates them which in turn causes big numerical inaccuracies in the polynomial root finding. 

Plus, MIMO system representations in both are simply asking for trouble since common terms are very difficult to retain and after a few manipulations you can easily create dummy roots -which could have been considered to be numerical noise- becoming actual modes of the representations. 

Similarly, canonical representations suffer from the same problem in state space. You can test it via creating 2 identical A matrices and perturb the entries of the row/col that holds the entries. Then create a balanced realization of the same system and perturb the A matrix entry of that representation and you'll see that the difference is usually very big.

In your comment you mention the nonuniqueness of the SS representation and it is true. But the dynamics are identical, hence they form an equivalence class up to transformations with nonzero transformation matrix $T$. That's actually a feature not a bug.