Lattice kinetic equations for simulating incompressible magnetohydrodynamics in two or three dimensions are constructed. The fluid is simulated via a conventional low Mach number lattice Boltzmann scheme, modified to include the Lorentz force due to the magnetic field. The magnetic field is represented by a separate vector-valued magnetic distribution function which obeys a vector Boltzmann–BGK equation. The two distribution functions are only coupled via the macroscopic density, momentum, and magnetic field evaluated at lattice points. This allows a reduced lattice to be used for the magnetic distribution function, with a corresponding saving in storage, which becomes comparable to that for the scalar hydrodynamic distribution function. The magnetic diffusivity may be adjusted independently of the fluid viscosity, unlike an earlier formulation. Numerical experiments with Hartmann flow, the Orszag–Tang vortex, and the doubly periodic coalescence instability compare favorably with results obtained using a spectral method, and with previously published results. The scheme preserved a consistent approximation to the divergence-free condition ∇·B=0 to round-off error.