This thesis develops and analyzes novel randomized algorithms for solving various classes of joint eigenvalue problems. Unlike the relatively well-studied standard or generalized eigenvalue problems, which involve only one or two matrices, joint eigenvalue problems deal with families of matrices and aim to extract shared structures such as common eigenvectors and joint eigenvalues. The central idea is to employ random linear combinations to reduce them to instances of the standard or generalized eigenvalue problems. This reduction enables the use of well-established eigenvalue solvers---built on decades of advances in numerical linear algebra---as reliable black-box tools for tackling these less-explored joint eigenvalue problems.

The first class of joint eigenvalue problems considered in this thesis is joint diagonalization (JD) of (nearly) commuting families of symmetric matrices, a fundamental task that frequently arises in signal processing and machine learning. We propose to diagonalize the input family by diagonalizing one random linear combination of the family. As a by-product, this approach yields an efficient randomized algorithm for diagonalizing normal matrices, achieved by diagonalizing one random linear combination of their Hermitian and skew-Hermitian parts. Moreover, the approach extends to multimodal spectral clustering via a new selection criterion that directly quantifies cross-modality smoothness of embeddings obtained from random linear combinations.

The second class of joint eigenvalue problems considered is simultaneous diagonalization by congruence (SDC) of symmetric matrices, a natural generalization of joint diagonalization (JD) that also arises frequently in signal processing. We propose to diagonalize the input family via congruence by solving a generalized eigenvalue problem formed from two random linear combinations of the family. Numerical evidence indicates that optimization-based algorithms, such as quasi-Newton methods, are promising candidates for efficiently refining this randomized approach.

The third class concerns the computation of joint eigenvalues of general commuting matrices, a problem central to multivariate polynomial system solving and multidimensional parameter identification in signal processing. We propose to approximate the joint eigenvalues by forming Rayleigh quotients of eigenvectors derived from one random linear combination.

Throughout this thesis, we establish theoretical guarantees for the exactness and robustness achieved by using random linear combinations across the different classes of joint eigenvalue problems considered. We further emphasize the practical side by demonstrating the efficiency and accuracy of the proposed algorithms through extensive numerical experiments on both synthetic and real-world data.