ENH: Add the idea of iterator axis to SumMultiply

The idea is that one plate axis could be used as an iterator. This
allows SumMultiply to compute the result in smaller blocks instead of a
possibly huge resulting array. This iterator could be used, for
instance, in GaussianMarkovChain for the state dynamics matrix if it
is different at each time instance. Because the state dynamics matrix is
needed only for one time instance at a time.
ENH: Implement rotation speed-ups for Gaussian arrays.

A new transformation/rotation class RotateGaussianArrayARD computes the
rotation for GaussianArrayARD nodes. These rotations can be used for
speeding up the inference of linear state-space models.
DEMO: Update PCA demo to use rotations.

PCA demo (demos/demo_pca.py) uses GaussianArrayARD which are now
rotated in order to speed up inference.
ENH: Implement plate rotation speed-up for Gaussian array node.

Rotation transformations can speed up the inference of linear
state-space models significantly. The rotation requires the rotation of
the state dynamics matrix which leads to also rotating the plate axes in
addition to the variable dimension axes. This plate rotation is now
implemented in GaussianArrayARD for performing the rotation and
RotateGaussianArrayARD for optimizing the rotation. Previous
implementation in RotateGaussianARD is deprecated because the new
implementation is more general.

See: Fast Variational Bayesian Linear State-Space Model (Luttinen, ECML
2013)
DEMO: PCA and LSSM demos use the improved rotation speed ups.
TST: Add a few tests for the plate rotation of Gaussian array nodes.
BUG: Fix the shape of scalar Gaussian array random samples

The shape of random samples from GaussianArrayARD were wrong in a
particular case. If GaussianArrayARD was used for scalars which had axes
variable dimensions (of length 1), they were incorrectly handled and the
result had too many unit axes.
MAINT: Remove a debugging message in GaussianArrayARD
MAINT: Improve error message for incorrectly shaped intialization
MAINT: Add more verbosity to optimize.check_gradient function
ENH: Add random sampling method to gamma node.
TST: Add test for rotating Gaussian ARD arrays.
BUG: Fix CGF in the rotation method of Gaussian ARD arrays

Rotate method was incorrect for Gaussian arrays with more than one
variable axes. The update of CGF in the rotation of Gaussian ARD arrays
did not take into account other variable axes. Now, CGF is multiplied by
the product of variable axes that were not rotated.
MAINT: Remove old un-used code from Gaussian ARD array

New implementation of rotate_plates has replaced the old one which has
been lying there un-used. Cleaning the code by removing it.
ENH: Add atol and rtol keyword parameters to TestCase.assertAllClose
ENH: In optimize.check_gradient, handle zero gradients properly
BUG: Fix speed-up rotations for Gaussian ARD arrays.

Rotation transformations for speeding up VB learning of Gaussian ARD
arrays contained several bugs when using the transformations for
variables with multiple variable axes or plate axes. In addition, alpha
hyperparameter was not handled properly unless it had plates only for
the rotated variable axis. This collection of bugs should now be fixed
and RotateGaussianArrayARD can be used for large range of different
model blocks.
TST: Add unit tests for speed-up rotation of Gaussian ARD arrays.

These unit tests check that the rotation bound is computed correctly and
that the gradient is correct. These checks are performed also for the
plate rotation. Tests are in test_transformations module.
DEMO: Add command line options to PCA demo

PCA demo (bayespy/demos/demo_pca.py) can be run as
    python demo_pca.py <options>
with several command line arguments, such as dimensionality of the data
or the latent space and whether to use speed-up rotations or not.
DEMO: Add command line options to LSSM demo

Linear state-space model (LSSM) demo (bayespy/demos/demo_lssm.py) can be
run as
    python demo_lssm.py <options>
with several command line arguments, such as dimensionality of the data
vectors or the latent space of the model and whether to use speed-up
rotations or not. It is highly recommended to use rotations, otherwise
it may take a very long time to reach convergence.