Magnetic resonance imaging (MRI) is a key diagnostic technique in healthcare nowadays, and of central importance for experimental research of the brain. However, in the present technique, imaging speed has to be limited in order not to harm patients, and cannot be significantly scaled up by parallelization or additional hardware. To decrease scan time, images have to be reconstructed from undersampled data, which is possible due to recent advances in the understanding of low level image statistics and sparse estimation. In this context, we show how improved MRI sequences can be found through optimization of Bayesian design scores. Our approach is driven by novel algorithms for approximate inference, which run orders of magnitude faster than previous methods, and allow for the estimation of posterior covariances over full-size images. Scalability is achieved by reduction to well-studied primitives of numerical mathematics and signal processing, such as least squares estimation, the Lanczos algorithm, or non-equispaced FFT, which can be scaled up dramatically through the use of commodity graphics hardware. Present applications of our framework lie in offline sequence optimization, but with further advances in achievable, highly parallel digital computation, we can envision dynamic real-time MRI settings, based on semi-automatic Bayesian decision support.