Having the rendering equation one can rewrite it in operator form:
$$L = L_e + TL$$
$$(I-T)L = L_e$$
$$L = (I-T)^{-1}L_e$$
$$L = \sum_{k=0}^{\infty}T^kL_e$$
The last equality holds if $T$ is a contraction ($\|T\| < 1$), and under very specific conditions if $T$ is not a contraction. An example of when the 4th equality does not hold would be picking an enclosed volume bounded by some surface, that is made of a non-energy conserving material, and putting a light source inside. Then the radiance will keep increasing with every bounce so at some point $L$ will be infinity everywhere on the inside of the volume. Meaning that the the series does not converge.
The 4th equality is really the Liouville-Neumann expansion. So you get an infinite sum of increasingly-dimensional integrals, each subsequent term being a bounce away from the previous. Each term gives you the contribution due to a light path of length equal to the index of the term.
More precisely, this solution gives you the steady state $L$ if it exists (we assume that light travels instantaneously, that it is incoherent, linear etc). Note that $L$ is simply the radiance function, so this formalism really gives you the solution at every point in the scene.
In practice you usually need only the radiance measured at the surface of a virtual camera film. So a straightforward optimization is to trace rays starting from the surface of the camera film. This does not guarantee that you get the solution for $L$ at every point of the scene, but is a lot more efficient. You can also trace rays starting at the light source, but you are usually also interested in their contribution to the film, so the measurement function (see Veach's thesis) has to also make it in there, and only a subset of lights may be of interest. Finally you can do photon mapping, or light tracing where you are interested of the radiance over the whole scene, and that would be the closest thing to the original formulation.