I think it's important to emphasize a point which was never emphasized to me when I was taking my courses. When a theory such as electromagnetism can be formulated in (roughly) equivalent ways at varying levels of sophistication, the more sophisticated forms should make your life **easier**, not harder.

Now, generally speaking this benefit only comes when you actually *understand* the added sophistication from both a physical and mathematical point of view, so for a time it may seem more complicated. But when that understanding has been achieved, there will be a moment where you suddenly see through the symbols and recognize everything as simpler than you first thought. The higher mathematics is there to serve you, not the other way around.

---

To put this in practice with electromagnetism, start from a relativistic standpoint and imagine that there was a force which is a linear function of the 4-velocity - that is, something of the form

$$\dot{\mathbf p} = qF \mathbf u$$

where $\mathbf p = m\mathbf u$ is the 4-momentum, $\mathbf u$ is the 4-velocity, and $q$ is some coupling constant. What are the properties of this $F$, which must be expressed as a matrix of some kind? Well, the first thing to point out is that because $\mathbf p\cdot \mathbf p = m^2c^4$ is constant, it follows that $\mathbf p \cdot \dot{\mathbf p} = 0$. But this requires that 
$$\mathbf u \cdot F\mathbf u = u^\mu\big(g_{\mu \alpha} F^\alpha_{\ \ \nu} \big)u^\nu \equiv F_{\mu\nu} u^\mu u^\nu =0$$
for any possible (timelike) 4-vector $\mathbf u$. The simplest way to satisfy this is to require that $F_{\mu\nu}$ be antisymmetric.

In any specific frame of reference, such a matrix acting on the 4-velocity can be broken up into a piece which depends on the ordinary 3-velocity $\vec v$ (in that frame) and a piece which does not.  It's a straightforward exercise to show that 

$$qF\mathbf u = q\gamma \pmatrix{\sum_{i=1}^3 F^0_{\ \ i} v^i \\ F^1_{\ \ 0} + \sum_{i=1}^3 F^1_{\ \ i} v^i \\  F^2_{\ \ 0} + \sum_{i=1}^3 F^2_{\ \ i} v^i \\  F^3_{\ \ 0} + \sum_{i=1}^3 F^3_{\ \ i} v^i }$$

This suggests that, for organizational purposes, we might define

$$F^i_{\ \ 0}\equiv E^i \qquad F^i_{\ \ j} = \mathcal B^i_{\ \ j}$$

where $\vec E$ is a 3-vector and $\mathcal B$ is a 3$\times$3 matrix with the property that $\mathcal B^i_{\ \ j} = -\mathcal B^j_{\ \ i}$ for each $i,j$. With this notation being defined, we have

$$F\mathbf u = q\gamma \pmatrix{\vec E \cdot \vec v \\ \vec E + \mathcal B\vec v}$$

Finally, observe that in 3 dimensions, the action of an **antisymmetric** matrix on a vector can be re-expressed as as a cross product:
$$\pmatrix{0 & \mathcal B^1_{\ \ 2} & \mathcal B^1_{\ \ 3} \\ -\mathcal B^1_{\ \ 2} & 0 & \mathcal B^2_{\ \ 3}\\ -\mathcal B^1_{\ \ 3} & -\mathcal B^2_{\ \ 3} & 0}\pmatrix{v^1\\v^2\\v^3} = \pmatrix{\mathcal B^1_{\ \ 2}v^2 + \mathcal B^1_{\ \ 3} v^3 \\ - \mathcal B^2_{\ \ 1} v^1 + \mathcal B^2_{\ \ 3}v^3\\ -\mathcal B^1_{\ \ 3} v^1 - \mathcal B^2_{\ \ 3} v^2 } = \vec v \times \vec B $$
where we have defined the 3-vector $\vec B \equiv (-\mathcal B^2_{\ \ 3}, \mathcal B^1_{\ \ 3}, -\mathcal B^1_{\ \ 2})$. So finally we have obtained

$$F\mathbf u = \gamma \pmatrix{q\vec E \cdot \vec v\\ q\big(\vec E + \vec v\times \vec B\big)}$$

---

The main lessons from this go as follows. 

 1. The fact that the force is linear in $\mathbf u$ and preserves the rest mass of the particle is sufficient to show that in any specific reference frame, it can be decomposed into timelike and spacelike parts - or put differently, into parts proportional to $\vec v$ and parts which are not. 
 2. This splitting is obviously frame-dependent, because $\vec v$ changes from frame to frame.
 3. The $\vec v$-dependent part is the action of an antisymmetric matrix $\mathcal B$ on $\vec v$. The fact that this can be expressed as the cross product of $\vec v$ with a vector $\vec B$ is an accident of 3-dimensions, essentially deriving from the fact that an $n\times n$ antisymmetric matrix has $n(n-1)/2$ independent components.
 4. If we want Lorentz-invariant quantities, we should obtain them at the level of $F$. Examples include $F^2 \equiv F_{\mu\nu}F^{\mu \nu} = 2(B^2-E^2/c^2)$ and $\epsilon^{\mu\nu\alpha\beta}F_{\mu\nu}F_{\alpha\beta} \propto \vec E \cdot \vec B$.

If we wanted to go further and endow this $F$ with dynamics of its own, we could start to think about the simplest, Lorentz-covariant differential equations we could write down. In what may or may not be surprising, the Maxwell equations take the rather straightforward form
$$\nabla_\mu F^{\mu\nu} = \mu_0 J^\nu \qquad \epsilon^{\mu \alpha\beta\gamma}\nabla_\alpha F_{\beta\gamma} = 0$$
where indices have been freely lowered and raised with the metric. The latter can be re-interpreted by defining the dual tensor
$$G^{\alpha\mu} = \frac{1}{2}\epsilon^{\alpha \mu \beta \gamma}F_{\beta \gamma}$$
in which case the second equation becomes $\nabla_\mu G^{\mu\nu} = 0$. These equations can be discussed and motivated at various levels as well, but that's beyond the scope of this answer.