  We develop a new homotopy method for solving multiparameter eigenvalue
problems (MEPs) called the fiber product homotopy method. For a $k$-parameter
eigenvalue problem with matrices of sizes $n_1,\dots ,n_k = O(n)$, fiber
product homotopy method requires deformation of $O(1)$ linear equations, while
existing homotopy methods for MEPs require $O(n)$ nonlinear equations. We show
that the fiber product homotopy method theoretically finds all eigenpairs of an
MEP with probability one. It is especially well-suited for dimension-deficient
singular MEPs, a weakness of all other existing methods, as the fiber product
homotopy method is provably convergent with probability one for such problems
as well, a fact borne out by numerical experiments. More generally, our
numerical experiments indicate that the fiber product homotopy method
significantly outperforms the standard Delta method in terms of accuracy, with
consistent backward errors on the order of $10^{-16}$, even for
dimension-deficient singular problems, and without any use of extended
precision. In terms of speed, it significantly outperforms previous
homotopy-based methods on all problems and outperforms the Delta method on
larger problems, and is also highly parallelizable. We show that the fiber
product MEP that we solve in the fiber product homotopy method, although
mathematically equivalent to a standard MEP, is typically a much better
conditioned problem.
