Joint diagonalization problems of Hermitian or non-Hermitian matrices occur as the final parameter estimation step in several blind source separation problems such as ACMA, JADE, PARAFAC, and SOBI. Previous approaches have been Jacobi iteration schemes and alternating projections. Here we show how the joint diagonalization problem can be formulated as a (weighted) subspace fitting problem so that it can be solved using the efficient Gauss-Newton optimization algorithm proposed for that problem. Since a good initial point is usually available, the algorithm converges very fast.