Nonlinear optics as gauge theory the widely used approach to nonlinear optics is a Taylor expansion of the dielectric displacement field $\mathbf{D} = \epsilon_0\cdot\mathbf{E} + \mathbf{P}$ in a Fourier representation of the polarization $\mathbf{P}$ in terms of the dielectric susceptibility $\mathcal{X}$:
$\mathbf{P} = \epsilon_0\cdot(\mathcal{X}^{(1)}(\mathbf{E}) + \mathcal{X}^{(2)}(\mathbf{E},\mathbf{E}) + \dots)$ .
This expansion does not work anymore if the excitation field has components close to the resonance of the medium. Then, one has to take the whole quantum mechanical situation into account by e.g. describing light/matter interaction by a two-level Hamiltonian.
But this approach is certainly not the most general one.
Intrinsically nonlinear formulations of electrodynamics
So, what kind of nonlinear formulations of electrodynamics given in a Lagrangian formulation are there?
One known ansatz is the Born-Infeld model as pointed out by Raskolnikov. There, the Lagrangian density is given by
$\mathcal{L} = b^2\cdot \left[ \sqrt{-\det (g_{\mu \nu})} - \sqrt{-\det(g_{\mu \nu} + F_{\mu \nu}/b)} \right]$
and the theory has some nice features as for example a maximum energy density and its relation to gauge fields in string theory. But as I see it, this model is an intrinsically nonlinear model for the free-space field itself and not usefull for describing nonlinear matter interaction.
The same holds for an ansatz of the form
$\mathcal{L} = -\frac{1}{4}F^{\mu\nu}F_{\mu\nu} + \lambda\cdot\left( F^{\mu\nu}F_{\mu\nu} \right)^2$
proposed by Mahzoon and Riazi. Of course, describing the system in Quantum Electrodynamics is intrinsically nonlinear and ... to my mind way to complicated for a macroscopical description for nonlinear optics. The question is: Can we still get a nice formulation of the theory, say, as a mean field theory via an effective Lagrangian?
I think a suitable ansatz could be
$\mathcal{L} = -\frac{1}{4}M^{\mu\nu}F_{\mu\nu}$
where $M$ now accounts for the matter reaction and depends in a nonlinear way on $\mathbf{E}$ and $\mathbf{B}$, say
$M^{\mu\nu} = T^{\mu\nu\alpha\beta}F_{\alpha\beta}$
where now $T$ is a nonlinear function of the field strength and might obey certain symmetries. The equation $T = T\left( F \right)$ remains unknown and depends on the material.
Metric vs. $T$ approach
As pointed out by space_cadet, one might ask the question why the nonlinearity is not better suited in the metric itself. I think this is a matter of taste. My point is that explicitly changing the metric might imply a non-stationary spacetime in which a Fourier transformation might not be well defined. It might be totally sufficient to treat spacetime as Lorentzian manifold.
Also, we might need a simple spacetime structure later on to explain the material interaction since the polarization $\mathbf{P}$ depends on the matter response generally in terms of an integration over the past, say
$\mathbf{P}(t) = \int_{-\infty}^{t}R\left[\mathbf{E}\right](\tau )d\tau$
with $R$ beeing some nonlinear response function(al) related to $T^{\mu\nu\alpha\beta}$.
Examples for $T$
To illustrate the idea of $T$, here are some examples.
For free space, $T$ it is given by $T^{\mu\nu\alpha\beta} = g^{\mu\alpha}g^{\nu\beta}$ resulting in the free-space Lagrangian $\mathcal{L} = -\frac{1}{4}T^{\mu\nu\alpha\beta}F_{\alpha\beta}F_{\mu\nu} = -\frac{1}{4}F^{\mu\nu}F_{\mu\nu}$
The Lagrangian of Mahzoon and Riazi can be reconstructed by
$T^{\mu\nu\alpha\beta} = \left( 1 + \lambda F^{\gamma\delta}F_{\gamma\delta} \right)\cdot g^{\mu\alpha}g^{\nu\beta}$.
One might be able to derive a Kerr nonlinearity using this Lagrangian.
So, is anyone familiar in a description of nonlinear optics/electrodynamics in terms of a gauge field theory or something similar to the thoughts outlined here?
Thank you in advance.
Sincerely,
Robert
Comments on the first Bounty
I want to thank everyone actively participating in the discussion, especially Greg Graviton, Marek, Raskolnikov, space_cadet and Willie Wong. I am enjoying the discussion relating to this question and thankfull for all the nice leads you gave. I decided to give the bounty to Willie since he gave the thread a new direction introducing the material manifold to us.
For now, I have to reconsider all the ideas and I hope I can come up with a new revision of the question that should be formulated in a clearer way as it is at the moment.
So, thank you again for your contributions and feel welcome to share new insights.
 A: Just a few random thoughts. 
There is something important in your observation that the Born-Infeld model is essentially a free-space model. It is known to Boillat and Plebanski (separately in 1970) that the Born-Infeld model is the only model of electromagnetism (as a connection on a $U(1)$ vector bundle) that satisfies the following conditions


*

*Covariance under Lorentz transformations

*Reduces to Maxwell's equation in the small-field strength limit

*$U(1)$ gauge symmetry

*Integrable energy density for a point-charge

*No birefringence (speed of light independent of polarization). 


(the linear Maxwell system fails condition 4.) (See Michael Kiessling, "Electromagnetic field theory without divergence problems", J. Stat. Phys. (2004) doi:10.1023/B:JOSS.0000037250.72634.2a for an exposition on this and related issues.)
Now, since you are interested in nonlinear optics inside a material, instead of in vacuum, I think conditions 1 and 5 can safely be dropped. (Though you may want to keep 5 as a matter of course.) Condition 4 is intuitively pleasing, but maybe not too important, at least not until you have some candidate theories in mind that you want to distinguish. Condition 3 you must keep. Condition 2, on the other hand, really depends on what kind of material you have in mind. 
In any case, a small suggestion: personally I think it is better to, from the get-go, write your proposed Lagrangian as
$$ L = T^{abcd} F_{ab}F_{cd} $$
instead of $M^{ab}F_{cd}$. I think it is generally preferable to consider Lagrangian field theories of at least quadratic dependence on the field variables. A pure linear term suggests to me an external potential which I don't think should be built into the theory. 
If you want something like condition 2, but with a dielectric constant or such, then you must have that $T^{abcd}$ admit a Taylor expansion looking something like
$$ T^{abcd} = \tilde{g}^{ac}\tilde{g}^{bd} + O(|F|) $$
where $\tilde{g}$ is some effective metric for the material. Birefringence, however, you don't have to insert in explicitly: most likely a generic (linear or nonlinear) $T^{abcd}$ you write down will have birefringence; it is only when you try to rule it out that you will bring in some constraints. 
An interesting thing is to consider what it means to have an analogous notion to condition 1. In the free-space case, condition 1 implies that the Lagrangian should only be a function of the Lorentz invariant $B^2 - E^2$ (in natural units) and of the pseudo-scalar invariant $B\cdot E$. In terms of the Faraday tensor these two invariants are $F^{ab}F_{ab}$ and $F^{ab}{}^*F_{ab}$ respectively, where ${}^*$ denote the Hodge dual. The determination of the linear part of your theory (of electromagnetic waves in a material) is essentially by what you will use to replace condition 1. If you assume your material is isotropic and homogeneous, then some similar sort of scalar + pseudo-scalar invariants is probably a good bet. 
A: Nonlinear is a buzzword used to cover anything that is not linear. Depending on what kind of nonlinearity is involved, and thus what kind of material, there could be one symmetry or another, or there could be no symmetry at all. For instance, in superconductors, gauge symmetry is broken and photons behave as if they have acquired a mass. The result is that magnetic fields have limited penetration in the superconductor. And I think this is still described by linear equations.
I know of one gauge-invariant theory that is non-linear, this model is called the Born-Infeld model.
A: You have been asking some seriously interesting questions! Here's my take on this one ...
You say this about the Born-Infeld action:

But as I see it, this model is an intrinsically nonlinear model for the free-space field itself and not useful for describing nonlinear matter interaction.

I'm not sure exactly what you mean by "free-space" field. I take it that you're referring to $ F_{\mu\nu} $. Well there is no reason why one cannot define an $ F_{\mu\nu} $ for waves propagating non-linearly, within a medium or in a vacuum.
The matter-light interaction can be specified (at least in part if not wholly) by the form of $ g_{\mu\nu} $. Now bear with me for a minute. I'm not referring to the metric generated by some kind of matter. The metric in question does not, a priori, satisfy the Einstein equations. It is instead the effective metric experienced by the light-rays propagating within the given material. See these excellent papers by Ulf Leonhardt and Thomas Philbin [1],[2] for more details on this notion. In brief the off-diagonal components $ g_{ij}$ (where $ (i,j \in \{1,2,3\}\,\, i \neq j) $ encode the susceptibility tensor and the diagonal components $ g_{0i} $ determine the mixing between the electric and magnetic components of the wave.
As for the lagrangian density for the matter-light interaction you posit:
$$ \mathcal{L}_{int} \propto M^{\mu\nu} F_{\mu\nu} = T^{\mu\nu\alpha\beta} F_{\alpha\beta} F_{\mu\nu} $$
for flat space (or no-medium) $ T^{\mu\nu\alpha\beta} = g^{\mu\nu}g^{\alpha\beta} $, this term reduces to $ F^{\mu\nu} F_{\mu\nu}$ which is nothing more than the Maxwell term ! On the face of it this gives us nothing new, unless we adopt the route outlined above and use the metric $g_{\mu\nu}$ to encode the optical properties of the medium.
Another line of thought which exploits this notion of the metric to allow one to speak of an analogy between optical processes and the big-bang is the phenomenal work of Igor Smolyaninov [3]. This paper was accepted by PRL btw, so its nothing to sneeze at.
Assuming that the above line of reasoning is not fatally flawed, and that one can encode the effects of the medium in the metric, it seems that either the Maxwell or the Born-Infeld action are  perfectly good candidates of gauge-invariant actions for your purposes.
                                Cheers,


Edit: Non-linearity redux
As @Raskolnikov pointed out, the identification of the components $g_{ab}$ with the optical susceptibilities of a material, does not give us a nonlinear material. For that, you have to have a dependence of the susceptibilities on the field strengths themselves. So you have a feedback mechanism $ \mathbf{g} \rightarrow \mathbf{F} \rightarrow \mathbf{g} $ and therefore the non-linearity ! Therefore in general, as @robert has been trying to convey to me without success, $\mathbf{g}$ should in general be a function of $\mathbf{F}$.
But then you start treading dangerously close to the speculation that somehow the eventual picture (for the fully non-linear case) might be somehow general relativistic. That is a very tempting idea, but I leave that for another time.
A: In a condensed matter field theory course, I learned the following: microscopically, the Lagrangian for the electromagnetic field looks like it is supposed to, coupling minimally to the particle coordinates.
$$ L = \sum_i\left( \frac m2 (p_i-\frac ec \mathbf A(r_i))^2 - e\Phi(r_i) + \dots \right) .$$
On a macroscopic level, however, after getting rid of all the individual particle degrees of freedom via the grand canonical ensemble, new behavior may emerge. Namely, the effective Lagrangian for the electromagnetic field in the body may look very different from a linear one. For example, the effective action for the e.m. field in a superconductor is
$$ S_{\text{eff}}[\mathbf A] = \frac\beta2 \int d^3r \mathbf A^\perp(r) \left(-\frac 1{\mu_0}\nabla^2 + \frac {n_s}m \right)\mathbf A^\perp(r)$$
where $\mu_0$ is the vacuum permeability, $n_s$ the superfluid density, $m$ the electron mass and $\mathbf A^\perp$ is the perpendicular component of the gauge field, defined in Fourier space as $\mathbf A^\perp(q) = \mathbf A(q) - q(q\cdot \mathbf A(q))/q^2$. The difference to the vacuum action is the additional "mass term" $n_s/m$, which causes the Meissner effect.

I suppose that you are asking for the most general form that such effective actions may have? I don't have an answer, but I don't see why a most general form should actually exist in the first place.
