Hybrid zero dynamics (HZD) has emerged as a popular framework for dynamic and underactuated bipedal walking, but has significant implementation difficulties when applied to the high degrees of freedom present in humanoid robots. The primary impediment is the process of gait design-it is difficult for optimizers to converge on a viable set of virtual constraints defining a gait. This paper presents a methodology that allows for the fast and reliable generation of efficient multi-contact robotic walking gaits through the framework of HZD, even in the presence of underactuation. To achieve this goal, we unify methods from trajectory optimization with the control framework of multi-domain hybrid zero dynamics. By formulating a novel optimization problem in the context of direct collocation and generating analytic Jacobians for the constraints, solving the resulting nonlinear program becomes tractable for large-scale nonlinear programming solvers, even for systems as high-dimensional as humanoid robots. We experimentally validated our methodology on the spring-legged prototype humanoid, DURUS, showing that the optimization approach yields dynamic and stable 3D walking gaits.