The class of matrix optimization problems (MOPs) has been recognized in recent years to be a powerful tool by researchers far beyond the optimization community to model many important applications involving structured low rank matrices. This trend can be credited to some extent to the exciting developments in the emerging field of compressed sensing. The L\"owner operator, which generates a matrix valued function by applying a single-variable function to each of the singular values of a matrix, has played an important role for a long time in solving matrix optimization problems. {However, the classical theory developed for L\"owner operators has become inadequate in these recent applications.} The main objective of this paper is to provide some necessary theoretical foundations for designing numerical methods for solving the MOP. This goal is achieved by introducing and conducting a thorough study on a new class of matrix valued functions, coined as spectral operators of matrices. Several fundamental properties of spectral operators, including the well-definedness, continuity, directional differentiability, Fr\'{e}chet-differentiability, locally Lipschitzian continuity, $\rho$-order B(ouligand)-differentiability ($0<\rho\leq 1$), $\rho$-order G-semismooth ($0<\rho\leq 1$) and the characterization of Clarke's generalized Jacobian, are systematically studied.