We introduce a method to infer a variational approximation to the posteriordistribution of solutions in computational imaging inverse problems. Machinelearning methods applied to computational imaging have proven very successful,but have so far largely focused on retrieving a single optimal solution for agiven task. Such retrieval is arguably an incomplete description of thesolution space, as in ill-posed inverse problems there may be many similarlylikely reconstructions. We minimise an upper bound on the divergence betweenour approximate distribution and the true intractable posterior, therebyobtaining a probabilistic description of the solution space in imaging inverseproblems with empirical prior. We demonstrate the advantage of our technique inquantitative simulations with the CelebA dataset and common imagereconstruction tasks. We then apply our method to two of the currently mostchallenging problems in experimental optics: imaging through highly scatteringmedia and imaging through multi-modal optical fibres. In both settings wereport state of the art reconstructions, while providing new capabilities, suchas estimation of error-bars and visualisation of multiple likelyreconstructions.