solving large nonlinear systems in parallel
The main point of Krylov-Newton is that it does not require computing and storing the whole Jacobian matrix. It only requires the ability to apply the matrix to a vector. 

The stardard approach is to do a little bit of pencil and paper work to figure out how to express the action of the derivative without actually computing and storing the whole matrix. If this is not possible, one can try to figure out a compressed representation of the matrix.

For problems where the matrices are dense, but the long range interactions are low rank (I think many radiative transfer problems would fall in this category..?), fast multipole or H-matrix type methods tend to be very effective.

It's difficult to give further advice without knowing more specifics of the problem.