The problem of numerical computation of a few Lyapunov exponents (LEs) of finite-dimensional dynamical systems is considered from the viewpoint of the differential geometry of Stiefel manifolds. Whether one computes one, many or all LEs of a continuous dynamical system by time integration, discrete or continuous orthonormalization is essential for stable numerical integration. A differential-geometric view of continuous orthogonalization suggests that one restricts the linearized vectorfield to a Stiefel manifold. However, the Stiefel manifold is not in general an attracting submanifold of the ambient Euclidean space: it is a constraint manifold with a weak numerical invariant. New numerical algorithms for this problem are then designed which use the fiber-bundle characterization of these manifolds. This framework leads to a new class of systems for continuous orthogonalization which have strong numerical invariance properties and the strong skew-symmetry property. Numerical integration of these new systems with geometric integrators leads to a new class of numerical methods for computing a few LEs which preserve orthonormality to machine accuracy. This idea is also taken a step further by making the Stiefel manifold an attracting invariant manifold in which case standard explicit Runge–Kutta algorithms can be used. This leads to an algorithm which requires only marginally more computation than a standard explicit integration without orthogonalization. These class of methods should be particularly effective for computing a few LEs for large-dimension dynamical systems. The new schemes are straightforward to implement. A test case is presented for illustration, and an example from dynamical systems is presented where a few LEs are computed for an array of coupled oscillators.