  Given a family of nearly commuting symmetric matrices, we consider the task
of computing an orthogonal matrix that nearly diagonalizes every matrix in the
family. In this paper, we propose and analyze randomized joint diagonalization
(RJD) for performing this task. RJD applies a standard eigenvalue solver to
random linear combinations of the matrices. Unlike existing optimization-based
methods, RJD is simple to implement and leverages existing high-quality linear
algebra software packages. Our main novel contribution is to prove robust
recovery: Given a family that is $\epsilon$-near to a commuting family, RJD
jointly diagonalizes this family, with high probability, up to an error of norm
O($\epsilon$). We also discuss how the algorithm can be further improved by
deflation techniques and demonstrate its state-of-the-art performance by
numerical experiments with synthetic and real-world data.
