The numerical computation of Lagrangian invariant subspaces of large scale Hamiltonian matrices is discussed in the context of the solution of Lyapunov and Riccati equations. A new version of the low-rank alternating direction implicit method is introduced, which in order to avoid numerical difficulties with solutions that are of very large norm, uses an inverse-free representation of the subspace and avoids inverses of ill-conditioned matrices. It is shown that this prevents large growth of the elements of the solution which may destroy a low-rank approximation of the solution. A partial error analysis is presented and the behavior of the method is demonstrated via several numerical examples
