title: PyMC3 (or Theano) allocating too much virtual memory while sampling
body: <p>I'm running a PyMC3 sampler on a remote cluster with a total of 3000 tuning steps, 3000 draws, and 2 chains. The data that I'm trying to fit consists of ~6000 data points and I'm setting up 7 variables (I think) to be optimized by the sampler. Maybe around 20 steps into the sampling, my program crashes because it hits the virtual memory allotment of <strong>15G</strong>. I can increase the allotment but then my program sits in the job queue on the cluster for much longer before running, which is not ideal.</p>
  
  <p>I'm relatively unfamiliar with the intricacies of PyMC3 and Theano, but I suspect one of these two is responsible. I've included the code for my PyMC3 model (which incorporates a separate astrophysics program, <em>exoplanet</em>) and the code where I set up the sampling. I have the Theano automatic garbage collection already turned on, but it doesn't seem to help much. Thanks for any help you can provide!</p>
  
  <p>This is part of my current script, utilizing PyMC3 + Theano with another astrophysics program <em>exoplanet</em>:</p>
  
  <pre><code>import numpy as np
  import pymc3 as pm
  import exoplanet as xo
  
  if mask is None:
      mask = np.ones(len(x), dtype=bool)
  with pm.Model() as model:
  
      mean = pm.Normal("mean", mu=0.0, sd=1.0)
  
      BoundedNormal_t0 = pm.Bound(pm.Normal, lower=t0_true-0.5, upper=t0_true+0.5)
      t0 = BoundedNormal_t0("t0", mu=t0_true, sd=1.0, shape=1)
  
      BoundedNormal_per = pm.Bound(pm.Normal, lower=(per_true*0.9), upper=(per_true*1.1))
      period = BoundedNormal_per("period", mu=per_true, sd=1.0, shape=1)
  
      u = pm.Uniform("u", lower=0.0, upper=1.0, shape=2)
      rho = pm.Uniform("rho", lower=0.0, upper=1000, shape=1)
  
      BoundedNormal_logr = pm.Bound(pm.Normal, lower=-10.0, upper=0.0)
      logr = BoundedNormal_logr("logr", mu=np.log(rp_rs_true), sd=10.0, shape=1)
      r = pm.Deterministic("r", tt.exp(logr))
  
      b = xo.distributions.ImpactParameter("b", ror=r, shape=1)
  
      # Set up a Keplerian orbit for the planets
      orbit = xo.orbits.KeplerianOrbit(period=period, t0=t0, b=b, rho_star=rho)#, r_star=rstar)#, ecc=ecc, omega=omega)
  
      # Compute the model light curve using starry
      light_curves = xo.LimbDarkLightCurve(u).get_light_curve(
          orbit=orbit, r=r, t=x[mask], texp=texp)
      light_curve = pm.math.sum(light_curves, axis=-1) + mean
  
      # The likelihood function assuming known Gaussian uncertainty
      pm.Normal("obs", mu=light_curve, sd=yerr[mask], observed=y[mask], dtype='float32')
  
      map_soln = xo.optimize(start=start)
      map_soln = xo.optimize(start=map_soln, vars=[logr], verbose=False)
      map_soln = xo.optimize(start=map_soln, vars=[b], verbose=False)
      map_soln = xo.optimize(start=map_soln, vars=[period, t0], verbose=False)
      map_soln = xo.optimize(start=map_soln, vars=[u], verbose=False)
      map_soln = xo.optimize(start=map_soln, vars=[logr], verbose=False)
      map_soln = xo.optimize(start=map_soln, vars=[b], verbose=False)
      map_soln = xo.optimize(start=map_soln, vars=[rho], verbose=False)
      map_soln = xo.optimize(start=map_soln, vars=[mean], verbose=False)
      map_soln = xo.optimize(start=map_soln, verbose=False)
  
  np.random.seed(42)
  with model:
      tr = pm.sample(chains=2
                     start=map_soln,
                     step=xo.get_dense_nuts_step(target_accept=0.95,start=map_soln),
                     tune=3000,
                     draws=3000)
  </code></pre>