This article investigates the methodology and applicability of the statistical linearization (SL) method to incorporating multi-variate non-differentiable nonlinearities, with a focus on floating renewable energy devices. The SL method serves as a highly competitive approach for analyzing floating renewable energy structures, such as wave energy converters (WECs) and floating wind energy turbines, because it inherently combines adequate accuracy and high computational efficiency. The origin of high accuracy comes from its incorporation of nonlinear effects through statistically linearized representations. Yet, the statistically linearized solutions have only been derived and verified for a limited number of nonlinearities of floating renewable energy devices, mostly simply-formed and differentiable in their mathematical expressions. However, floating renewable energy devices usually exhibit a complex dynamic mechanism, in which the relevant nonlinear effects could appear to be highly complex for linearization process to describe. These nonlinear effects could make a significant impact on the system dynamics, exemplified by external machinery force saturation and nonlinear hydrostatics of floaters with a non-uniform geometry. To push forward the boundary of the SL method, it is crucial to demonstrate how it applies to nonlinearities of different features. In this paper, the existing SL method is extended to address the nonlinear effects expressed as multi-variate non-differentiable functions. Several case studies are carried out to exemplify the application of the extended SL approach to the concerned nonlinearities in floating renewable energy devices. The accuracy and computational efficiency of the extended SL approach are evaluated by verifying against the corresponding nonlinear time-domain (TD) and linear frequency-domain (FD) models. Despite the complexity of the given nonlinearities, the relative errors of the SL approach are no more than 6 % while its computational time is comparable to the FD model, being thousands of times faster than the TD model. Comparatively, the FD model leads to a relative error of over 70% in some cases.