Shape issues in pm.MvNormal

[Dpananos]

See here for more.

When attempting to simulate an Nxp random matrix in which the rows are draws from a multivariate normal, a shape issue pops up. The example below should

a) Generate a 10 x 2 array X

b) generate a 2 dimensional array, beta

c) Compute their matrix product to yield a 10 dimensional array y

import pymc3 as pm

N = 10
Sigma = np.eye(2)

with pm.Model() as model:
    
    X = pm.MvNormal('X', mu=np.zeros(2), cov = Sigma, shape = (N,2))
    betas = pm.Normal('betas', 0, 1, shape = 2)
    y = pm.Deterministic('y', pm.math.dot(X,betas))
    
    prior_pred = pm.sample_prior_predictive(1)

Please provide the full traceback.

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-18-0c1a293243c3> in <module>
     10     y = pm.Deterministic('y', pm.math.dot(X,betas))
     11 
---> 12     prior_pred = pm.sample_prior_predictive(1)

/opt/anaconda3/lib/python3.7/site-packages/pymc3/sampling.py in sample_prior_predictive(samples, model, vars, var_names, random_seed)
   1320     names = get_default_varnames(model.named_vars, include_transformed=False)
   1321     # draw_values fails with auto-transformed variables. transform them later!
-> 1322     values = draw_values([model[name] for name in names], size=samples)
   1323 
   1324     data = {k: v for k, v in zip(names, values)}

/opt/anaconda3/lib/python3.7/site-packages/pymc3/distributions/distribution.py in draw_values(params, point, size)
    437                                             point=point,
    438                                             givens=givens.values(),
--> 439                                             size=size)
    440                         evaluated[param_idx] = drawn[(param, size)] = value
    441                         givens[param.name] = (param, value)

/opt/anaconda3/lib/python3.7/site-packages/pymc3/distributions/distribution.py in _draw_value(param, point, givens, size)
    615                 input_vals = []
    616             func = _compile_theano_function(param, input_vars)
--> 617             output = func(*input_vals)
    618             return output
    619     raise ValueError('Unexpected type in draw_value: %s' % type(param))

/opt/anaconda3/lib/python3.7/site-packages/numpy/lib/function_base.py in __call__(self, *args, **kwargs)
   2089             vargs.extend([kwargs[_n] for _n in names])
   2090 
-> 2091         return self._vectorize_call(func=func, args=vargs)
   2092 
   2093     def _get_ufunc_and_otypes(self, func, args):

/opt/anaconda3/lib/python3.7/site-packages/numpy/lib/function_base.py in _vectorize_call(self, func, args)
   2155         """Vectorized call to `func` over positional `args`."""
   2156         if self.signature is not None:
-> 2157             res = self._vectorize_call_with_signature(func, args)
   2158         elif not args:
   2159             res = func()

/opt/anaconda3/lib/python3.7/site-packages/numpy/lib/function_base.py in _vectorize_call_with_signature(self, func, args)
   2185 
   2186         broadcast_shape, dim_sizes = _parse_input_dimensions(
-> 2187             args, input_core_dims)
   2188         input_shapes = _calculate_shapes(broadcast_shape, dim_sizes,
   2189                                          input_core_dims)

/opt/anaconda3/lib/python3.7/site-packages/numpy/lib/function_base.py in _parse_input_dimensions(args, input_core_dims)
   1854     dim_sizes = {}
   1855     for arg, core_dims in zip(args, input_core_dims):
-> 1856         _update_dim_sizes(dim_sizes, arg, core_dims)
   1857         ndim = arg.ndim - len(core_dims)
   1858         dummy_array = np.lib.stride_tricks.as_strided(0, arg.shape[:ndim])

/opt/anaconda3/lib/python3.7/site-packages/numpy/lib/function_base.py in _update_dim_sizes(dim_sizes, arg, core_dims)
   1820             '%d-dimensional argument does not have enough '
   1821             'dimensions for all core dimensions %r'
-> 1822             % (arg.ndim, core_dims))
   1823 
   1824     core_shape = arg.shape[-num_core_dims:]

ValueError: 1-dimensional argument does not have enough dimensions for all core dimensions ('i_0_0', 'i_0_1')

Please provide any additional information below.

Versions and main components

• PyMC3 Version: 3.7
• Theano Version: 1.0.4
• Python Version: 3.7
• Operating system: OSX 10.14.6
• How did you install PyMC3: conda

[lucianopaz]

Related to #2848.
I never had time to address this problem. It involves two big aspects:

1. If chol is given, we have to extend the shape_utils to account for parameters that can have different core ranks, like mu and chol that would have rank 1 and 2 respectively.
2. For other parametrizations, we have to make a batched version of theano's cholesky decomposition.

[lucianopaz]

Ah, sorry @Dpananos, maybe in your example you meant to write X = pm.MvNormal('X', mu=np.zeros(2), cov = Sigma, shape = (N, 2)) instead.
In your example, the distribution's shape doesn't match with the implied shape given the parameters you supplied.

[Dpananos]

@lucianopaz Oops, let me correct that

[Dpananos]

@lucianopaz Making the appropriate change now gives a different traceback.

[lucianopaz]

@Dpananos, could you edit your original code with the new traceback? That will be the relevant error, which is related to #2848

[Dpananos]

@lucianopaz yep, already made the edit. That is the new traceback

[junpenglao]

Also, if you do sample_prior_preditive(k) it samples but silently fail with output shape of 'y' being (k, k)

[canyon289]

Referencing this ticket here #2848

[AlexAndorra]

I'm having this same issue for my PyMCon talk -- sampling from the MvNormal prior predictive but it silently drops the first dimension (the N dimension in Demetri's example).
That'd be nice if I could fix this before submitting my PyMCon talk, but I'm quite time-constrained right now... @lucianopaz, how long and difficult you estimate the fix to be?

[lucianopaz]

@AlexAndorra, I haven't really thought how long it would take to solve the MvNormal problems but it won't be an easy fix.
I never came around to extending shape utils to handle multivariate event shapes because I had trouble thinking of a proper strategy to deal with batching. It would take long, and for now I don't have the bandwidth to do it. On the other hand, implementing a batched cholesky decomposition in theano is completely out of my league.

[AlexAndorra]

Thanks for the quick and detailed answer @lucianopaz ! I'd be interested in working on this too, but I don't have the bandwidth either currently 😩 That would definitely be a very welcome contribution, so I'll mark this issue as "help wanted".
In the meantime, I know how to work around this, so I'll show it in the tutorial 😉

[twiecki]

Closed by #4207.