  Multivariate normal (MVN) probabilities arise in myriad applications, but
they are analytically intractable and need to be evaluated via
Monte-Carlo-based numerical integration. For the state-of-the-art minimax
exponential tilting (MET) method, we show that the complexity of each of its
components can be greatly reduced through an integrand parameterization that
utilizes the sparse inverse Cholesky factor produced by the Vecchia
approximation, whose approximation error is often negligible relative to the
Monte-Carlo error. Based on this idea, we derive algorithms that can estimate
MVN probabilities and sample from truncated MVN distributions in linear time
(and that are easily parallelizable) at the same convergence or acceptance rate
as MET, whose complexity is cubic in the dimension of the MVN probability. We
showcase the advantages of our methods relative to existing approaches using
several simulated examples. We also analyze a groundwater-contamination dataset
with over twenty thousand censored measurements to demonstrate the scalability
of our method for partially censored Gaussian-process models.
