A frequency domain methodology is developed for stochastic response determination of multi-degree-of-freedom (MDOF) linear and nonlinear structural systems with singular matrices. This system modeling can arise when a greater than the minimum number of coordinates/DOFs is utilized, and can be advantageous, for instance, in cases of complex multibody systems where the explicit formulation of the equations of motion can be a nontrivial task. In such cases, the introduction of additional/redundant DOFs can facilitate the formulation of the equations of motion in a less labor intensive manner. Specifically, relying on the generalized matrix inverse theory, a Moore-Penrose (M-P) based frequency response function (FRF) is determined for a linear structural system with singular matrices. Next, relying on the M-P FRF a spectral input-output (excitation-response) relationship is derived in the frequency domain for determining the linear system response power spectrum. Further, the above methodology is extended via statistical linearization to account for nonlinear systems. This leads to an iterative determination of the system response mean vector and covariance matrix. Furthermore, to account for singular matrices, the generalization of a widely utilized formula that facilitates the application of statistical linearization is proved as well. The formula relates to the expectation of the derivatives of the system nonlinear function and is based on a Gaussian response assumption. Several linear and nonlinear MDOF structural systems with singular matrices are considered as numerical examples for demonstrating the validity and applicability of the developed frequency domain methodology.