  In this paper, we propose, analyze, and test an efficient algorithm for
computing ensemble average of incompressible magnetohydrodynamics (MHD) flows,
where instances/members correspond to varying kinematic viscosity, magnetic
diffusivity, body forces, and initial conditions. The algorithm is decoupled in
Els\"asser variables and permits a shared coefficient matrix for all members at
each time-step. Thus, the algorithm is much more computationally efficient than
separately computing simulations for each member using usual MHD algorithms. We
prove the proposed algorithm is unconditionally stable and convergent. Several
numerical tests are given to support the predicted convergence rates. Finally,
we test the proposed scheme and observe how the physical behavior changes as
the coupling number increases in a lid-driven cavity problem with mean Reynolds
number $Re\approx 15000$, and as the deviation of uncertainties in the initial
and boundary conditions increases in a channel flow past a step problem.
