I found an excellent survey paper that answers your question is full detail: http://arxiv.org/abs/0811.0882 . Thus I'll be brief in the following overview.
Every orbit of an autonomous dynamical system has an associated sequence of Lyapunov exponents. The first (or maximal) Lyapunov exponent $\lambda_1$ tells how fast the distance of points in an infinitesimal neighborhood of a point of the orbit increases or decreases with time in the limit of very long times. (''Very long'' may be a second for a microscopic system, or many hundred millions of years for a planetary system.)
The asymptotic law is exponential, with dominating term $ce^{\lambda_1 t}$. For example, a dissipative system that contracts to a point has $\lambda_1<0$, a distance-preserving system has $\lambda_1<0$, and a chaotic system has $\lambda_1>0$. The latter condition says that there is a very sensitive dependence on initial conditions, and is a bit more general than chaoticity. The other Lyapunov exponents express similar properties related to areas, 3-volumes, etc. in place of distances.
The Lyapunov exponents are independent of the initial condition on the orbit, but may be different for initial conditions defining different orbits. They are a property of the system only if the latter is topologically transitive (ergodic), so that all orbits behave essentially the same. (An example of the distribution of $\lambda_1$ for a large sample of initial conditions in an application in meteorology is given in Figure 9 of a paper by Shepherd et al..
For their computation of Lyapunov exponents in case of a conservative system with a given Hamiltonian one indeed needs to do a stability analysis, but in a globalized form. One linearizes the system in a neighborhood of the orbit under study, and represents the resulting linear, time-dependent system in a product form, writing (for systems with finitely many degrees of freedom) the evolution matrices $A(t)$ (that tell how an initial deviation vector $d(0)$ is mapped to $d(t)=A(t)d(0)$) in the form $A(t)=Q(t)R(t)$ with orthogonal $Q(t)$ and upper triangular $R(t)$. (This is a continuous version of the orthogonal factorization $A=QR$ of a matrix $A$, heavily used in numerical analysis.) The Lyapunov exponents can then be read off from the time-dependence of the diagonal elements of $R(t)$. To find $Q(t)$ and $R(t)$ one must solve a coupled system of ordinary differential equations for the components.