Form of the Classical EM Lagrangian
The Lagrangian for Electromagnetism follows uniquely from requiring renormalizability and gauge invariance (plus parity time reversal)
U(1) gauge Invariance
if you require your Lagrangian to be locally invariant under symmetry operations of the unitary group U(1) that is under
$$\phi\to e^{i\alpha(x)}\phi$$
all derivatives $\partial_\mu$ have to be replaced by the covariant derivative $D_\mu = \partial_\mu+ieA_\mu$, where, in order to save local invariance the gauge field is introduced. Loosely speaking this is necessary to make fields at different spacetiem points comparable Since two points may have an arbitrary phase difference, due to the fact that we can set $\alpha(x)$ as we wish, something has to compensate this difference, before we can compare fields, which is what differentiation basically does. This is similar to parallel transport in general relativity (the mathematical keyword is connection see wiki: Connection (wiki) The gauge field $A_\mu$ transforms as $A_\mu \to A_\mu - \frac{1}{e}\partial_\mu\alpha(x)$.
Now the question is what kind of Lagrangians we can build with this requirement. For matter (i.e. non-gauge) fields it's easy to construct gauge invariant quantities by just replacing the derivatives with the covariant derivatives, i.e.
$$\bar{\psi}\partial_\mu\gamma^\mu\psi\to \bar{\psi}D_\mu\gamma^\mu\psi$$,
this will yield kinetic terms for the field (the part with the normal derivative), and interactions terms between matter fields and the gauge field.
Gauge-Field only terms
the remaining question is how to construct terms involving only the gauge field and no matter fields (i.e. the 'source-free' terms your question is about). For this we must construct gauge-invariant germs of $A_\mu$.
Once $\alpha(x)$ is chosen we can imagine starting from a point and walking on a loop back to that same point (this is called a wilson loop (wiki)). This must necessarily be gauge invariant since any phase that we pick up on the way we must also loose on the way back. It turns out, that this is exactly the term $F_{\mu\nu}$, i.e. the field strength. (the calculation is a little longer, see Peskin & Schroeder page 484). Actually this is only true for abelian symmetries such as U(1), for non abelian ones such as SU(3) we will get some interaction terms between the gauge fields which is why light does not interact with itself but gluons do.
Bilinear mass terms such as $A_\mu A^\mu$ are not gauge invariant (in the end this is the need for the Higgs meachanism)
Renormalizability
If we wish that our theory is renormalizable, we can only include terms into the lagrangian up to mass dimension 4. Now listing all terms up to mass dimension 4 we arrive at
$$\mathcal{L} = \cdot\bar{\psi}D_\mu\psi - m\bar{\psi}\psi - b\cdot F_{\mu\nu}F^{\mu\nu} + d\cdot \epsilon^{\alpha\beta\gamma\delta}F_{\alpha\beta}F_{\gamma\delta}$$
the last term involves the anti-symmetric tensor $\epsilon^{\alpha\beta\gamma\delta}$ and is therefore not time and parity invariant.
Note that we have not included linear terms here since we will be expanding around a local minimum anyways, so that the linear term will vanish.
Conclusion
if we require U(1) gauge invariance and renormalizability (mass dimension up to 4) and time and parity invariance we only get
$$\mathcal{L} = \cdot\bar{\psi}D_\mu\psi - m\bar{\psi}\psi - b\cdot F_{\mu\nu}F^{\mu\nu}$$
In the source-free case this is
$$\mathcal{L} = - b\cdot F_{\mu\nu}F^{\mu\nu}$$
the overall factor $\frac{1}{4}$ is not important.
The Maxwell term
$$\tag{1} {\cal L}_{\rm Maxwell}~=~-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}$$
emerges naturally for many reasons.
1) Pure EM without matter. There is a very short list of special relativistic renormalizable terms, that one can put in a local Lagrangian density with no higher-order time derivatives, and that is gauge-invariant (up to boundary terms). Among the shortlist of these candidates, the Lagrangian density (1) is the only one (modulo an overall normalization factor and modulo boundary terms) that leads to Maxwell equations (without source terms). The Born-infeld Lagrangian density is an example of a non-local candidate.
2) EM coupled to matter. Additional conditions arise when one tries to couple EM to point charges. One may argue that the relativistic Lagrangian for $n$ point charges $q_1, \ldots, q_n$, at positions ${\bf r}_1, \ldots, {\bf r}_n$, in an EM background is given as
$$\tag{2} L ~=~ -\sum_{i=1}^n \left(\frac{m_{0i}c^2}{\gamma({\bf v}_i)} +q_i\{\phi({\bf r}_i) - {\bf v}_i\cdot {\bf A}({\bf r}_i)\} \right).$$
See also this Phys.SE answer. This Lagrangian (2) e.g. reproduces correctly the Lorentz force. From eq. (2)it is only a small step to conclude that the interaction term ${\cal L}_{\rm int}$ between EM and matter [in the $(-,+,+,+)$ sign convention] is of the form
$$\tag{3} {\cal L}_{\rm int}~=~J^{\mu}A_{\mu}.$$
Also recall that Maxwell's equations with sources are
$$\tag{4} d_{\nu} F^{\nu\mu}~=~-J^{\mu}. $$
If the action
$$\tag{5} S[A]~=~\int\! d^4x~ {\cal L} $$
is supposed to be varied wrt. the $4$-gauge potential $A_{\mu}$, i.e. the $4$-gauge potential $A_{\mu}$ are the fundamental variables of the theory, and if moreover the corresponding Euler-Lagrange equations
$$\tag{6} 0~=~\frac{\delta S}{\delta A_{\mu}} ~\stackrel{?}{=}~d_{\nu} F^{\nu\mu}+J^{\mu}$$
are supposed to reproduce Maxwell's equations (4), it quickly becomes clear that
$$\tag{7} {\cal L}~=~{\cal L}_{\rm Maxwell}+{\cal L}_{\rm int} $$
is the appropriate Lagrangian density for EM with background sources $J^{\mu}$.
3) For a discussion of formulating EM without the $4$-gauge potential $A_{\mu}$, see also this Phys.SE post.
The most economical derivation of the Maxwell equations (i.e., depending on the least number of postulates) is known nowadays by the name of "Feynman's proof of the Maxwell equations". This "proof" was discovered by Feynman in 1948, but was not published because Feynman did not think that it leads to new physics. It was only in 1989 after Feynman's death when his proof was published by Dyson, please see the Dyson's article.
After being published by Dyson, the Feynman's proof was found to contain very deep ideas in Poisson geometry. As well as its generalization to the non-Abelian case leads to the Wong equations describing the motion of a point particle in an external Yang-Mills field, (please see the following exposition by: Montesinos and Abdel Perez-Lorenzana) which in turn is related to the Kaluza-Klein theory. It explains various mechanical ideas, for example the phenomenon of the falling cat. The whole subject is known today by the name of sub-Riemannian geometry.
The postulates of the Feynman's proof are the following:
- The position and velocity of the particle satisfy the canonical Poisson bracket relation with the position.
- The acceleration of a charged particle in an electromagnetic field is a function of the position and velocity only.
Under these assumptions only, Feynman proved (please see the Dyson's article for details) that the electromagnetic force must satisfy the Lorentz law and the electromagnetic field must satisfy the homogeneous Maxwell equations.