  We derive a family of efficient constrained dynamics algorithms by
formulating an equivalent linear quadratic regulator (LQR) problem using Gauss
principle of least constraint and solving it using dynamic programming. Our
approach builds upon the pioneering (but largely unknown) O(n + m^2d + m^3)
solver by Popov and Vereshchagin (PV), where n, m and d are the number of
joints, number of constraints and the kinematic tree depth respectively. We
provide an expository derivation for the original PV solver and extend it to
floating-base kinematic trees with constraints allowed on any link. We make new
connections between the LQR's dual Hessian and the inverse operational space
inertia matrix (OSIM), permitting efficient OSIM computation, which we further
accelerate using matrix inversion lemma. By generalizing the elimination
ordering and accounting for MUJOCO-type soft constraints, we derive two
original O(n + m) complexity solvers. Our numerical results indicate that
significant simulation speed-up can be achieved for high dimensional robots
like quadrupeds and humanoids using our algorithms as they scale better than
the widely used O(nd^2 + m^2d + d^2m) LTL algorithm of Featherstone. The
derivation through the LQR-constrained dynamics connection can make our
algorithm accessible to a wider audience and enable cross-fertilization of
software and research results between the fields
