The recent success of - and demand for - compliant mechanisms has increased rapidly within the micro-electromechanical systems-, aircraft-, spacecraft-, surgical-, and precision-instrument industries. Yet, even greater success may be achieved by overcoming the computational cost and instability of the mechanisms' design methodology. This involves large-scale topology optimization, geometrically non-linear structural analysis, and particularly integrating the latter into the former. Therefore, a novel framework is proposed that extends the powerful design freedom of topology optimization with most of the geometrically non-linear qualities, without as much of the computational ramifications. Utilizing a Bayesian-enhanced perturbed analysis, the equilibrium curve is locally approximated by an asymptotic expansion, satisfying the curve's higher-order geometric derivatives at the undeformed state. Each of the latter is recursively and efficiently obtained through a linear solve with the same, Cholesky-factorized stiffness matrix. Furthermore, a tensor reformulation and -decomposition of the four-noded bilinear element's Green-Lagrange strain energy model are successfully exploited. A tight error estimator is derived to govern the Bayesian analysis and balance approximation efficiency and accuracy during optimization. Design-sensitivities of this error-estimator and other design-dependent responses are analytically obtained through the adjoint method, and applied in a few classical density-based topology optimizations; While Bayesian-enhanced perturbed analysis required noticeably less computational effort compared to Newton-Raphson analysis under mildly non-linear conditions, it often resulted in practically identical performance improvements compared to linear analysis.