How can macroscopic $\bf D$ be related to microscopic $\bf E$ just by a constant $\epsilon$? In electromagnetism, the electric displacement field $\bf D$ and magnetic field strength $\bf H$ are  macroscopic functions defined in matter, and $\bf E$ and $\bf B$ are microscopic functions.
How can the relation $\bf D=\epsilon \bf E$ hold while on one side is a macroscopic quantity and on the other side is a microscopic one?
 A: Hopefully the following will answer what you are getting at. 
In general macroscopic $\vec{D}$ is not related to $\vec{E}$ by a simple scalar proportionality constant $\epsilon$ (nor is $\vec{B}$ related to $\vec{H}$ a simple scalar proportionality constant $\mu$). You likely know that in general $\epsilon$ and $\mu$ are $3\times 3$ matrices because anistropic materials have preferred directions owing to molecular orientation that makes them more polarizable in certain directions than others. 
However, what you may not have heard of and which is a fact that I'm guessing might help you here is that $\epsilon$ and $\mu$ are not even simple matrix multipliers.
In general they can be  nonlocal linear operators defined by spatial convolutions, or even nonlocal, nonlinear spatial integral operators 
This is because the polarizing molecules and atoms in the material affect the nett fields in a sizeable neighbourhood around them, so the effective support of the kernels in such integral operators are of the order of a few molecule widths or, in the case of patterned arrangements of molecules, of the order of the widths of the "unit cells" of the pattern.
But these widths are generally much smaller than the shortest wavelength present in the electromagnetic fields of interest. For example, optical fields vary significantly only over tens of nanometers, whereas the molecular and pattern widths we are talking about are often a fraction of a nanometer across. Therefore the integral operators are well approximated by local constants: a photon spread over a half a micron interacts with billions of atoms at a time, so the medium really does look like a continuum for most optical fields and the kernel width in the integral operator is three orders of magnitude smaller than the photon's region of influence. Therefore, the $\epsilon$ and $\mu$ operators are, to an excellent approximation, local and can be represented by simple proportionality constants, or at least local $3\times3$ matrices. But this depends on the frequency of the field in question. Once we get to gamma frequencies, nonlocal effects will be highly significant and the summations over all the atoms in my expressions below will be highly dependent on position and nonlocal.
Just as an example (you're not meant to grasp this on inspection, it's simply to show you how complicated the general case can be), here is what a general linear $\epsilon$ operator looks like, as taken from my own paper in J. Opt. Soc. Am. B/Vol. 24, No. 4/ April 2007:
$$
{\rm{\varepsilon}}_0
  \left(
   {
    \begin{array}{l}
     1+
      \sum \limits_
      {{\rm{all\;atoms\;m}}}
      {
       \left(
        {
         \frac{
          \sin
          \left(
          {{\rm{\omega}}_{\rm{m}}t}
          \right)
          \ast_t \mathfrak{Re}
          \left(
          {\mathcal{L}_{11,{\rm{m}}}}
          \right)
          -\,\cos
          \left(
          {{\rm{\omega}}_{\rm{m}}t}
          \right)
          \ast_t \mathfrak{Im}
          \left(
          {\mathcal{L}_{11,{\rm{m}}}}
          \right)
         }{{\rm{\omega}}_{\rm{m}}}
        }
       \right)
      }
     + \\ \quad\quad
      \sum \limits_
      {{\rm{all\;atoms\;m}}}
      {
       \left(
        {
         \frac{
          \cos
          \left(
          {{\rm{\omega}}_{\rm{m}}t}
          \right)
          \ast_t \mathfrak{Re}
          \left(
          {\mathcal{L}_{11,{\rm{m}}}}
          \right)
          +\sin
          \left(
          {{\rm{\omega}}_{\rm{m}}t}
          \right)
          \ast_t \mathfrak{Im}
          \left(
          {\mathcal{L}_{11,{\rm{m}}}}
          \right)
         }{{\rm{\omega}}_{\rm{m}}}
        }
       \right)
      }
     \,\mathcal{H}_t
    \end{array}
   }
  \right)
$$
where $\mathcal{H}_t$ is the time Hilbert transform, $*_t$ stands for time convolution and $\mathcal{L}_{11}$ is the spatial convolution:
$$
\mathcal{L}_{11,{\rm{m}}}
  {
   \rm
   {\bf \overset{\smile}{{E}}}
  }
  \left({{\rm{\bf r}},t}\right)
  \,= \\ \frac{1}{2}
   \sum \limits_
   {{\rm{all\;atoms\;m}}}
   {
    \left(
     {
      \left(
       {
        {\rm{\bf C}}_{+{\rm{,m}}}^
        {\rm{\bf \dagger}}
        \left(
         {
          {\rm{\bf r}}_{\rm{m}}-
          {\rm{\bf r}}
         }
        \right)
        +{\rm{\bf C}}_{-{\rm{,m}}}^
        {\rm{\bf \dagger}}
        \left(
         {
          {\rm{\bf r}}_{\rm{m}}-
          {\rm{\bf r}}
         }
        \right)
       }
      \right)
      \,
      \left.
       {
        \left(
         {
          \,
          \left(
           {
            \mathfrak{L}_{BB}
            \left(
             {
              {\rm{\bf C}}_{+,{\rm{m}}}
              \left({\rm{\bf u}}\right)+
              {\rm{\bf C}}_{-,{\rm{m}}}
              \left({\rm{\bf u}}\right)
             }
            \right)
           }
          \right)
          \;\ast_{\rm{\bf r}}\,
          {
           \rm
           {\bf \overset{\smile}{{E}}}
          }
          \left({{\rm{\bf u}},t}\right)
          \,
         }
        \right)
       }
      \right|
      _
      {
       {\rm{\bf u}}={\rm{\bf r}}_
       {\rm{m}}
      }
     }
    \right)
   }
  $$
and $*_r$ stands for spatial convolution and the $\mathbf{C}$ matrices are spatial Fourier transforms of the coupling between an instance of the particular atom species in question and the free, plane wave photons.
A: In general, it doesn't. There are two approximations being made here:


*

*The medium is considered isotropic, i.e. it behaves the same in any direction. This is often very much not the case in crystals where there is long range order but holds for amorphous solids, liquids or gases. In the anisotropic case, the number $\epsilon$ is replaced by a $3\times 3$ matrix of the same name. Then, $\mathbf D$ and $\mathbf E$ will not be parallel in general. This gives rise to e.g. birefringence.

*Also, the linear relation is an approximation, too. You can think of it as the first order approximation in terms of a Taylor series in $|\mathbf E|$. For very large external excitations, we also have to consider the second order correction which goes with $|\mathbf E|^2$. This is used in second-harmonic generation where the frequency of some portion of a highly intense laser beam is doubled by sending it though what is called a non-linear crystal.

A: Even if the linear, isotropic, spatially homogeneous, approximation can be made, there are many situations (already in lenses), where $\varepsilon$ and $\mu$ are dependent on the frequency $\omega$, i.e., for dispersive and absorptive media (often $\mu$ can be set equal to its vacuum value)  This happens if
$${D}(\mathbf{x},t)={E}(\mathbf{x},t)+{P}(\mathbf{x},t),$$
where the polarisation ${P}(\mathbf{x},t)$ is given by
$${P}(\mathbf{x},t)=\int_{t_{0}}^{t}ds{\chi}(\mathbf{x},t-s)\cdot
{E}(\mathbf{x},s),$$
and similar for the magnetisation.
Fourier transforming then gives$$
\tilde{D}(\mathbf{x},\omega )=\varepsilon (\mathbf{x},\omega )\tilde{E}(\mathbf{%
x},\omega ).$$We see that the frequency dependence is a direct consequence of the time delay in ${\chi}(\mathbf{x},t-s)$.
Since dispersion and absorption are always present, the assumption of a frequency-independent $\varepsilon$ involves a further approximation.
As to your question, the relation between $D$ and $E$ in a macroscopic medium is usually derived in textbooks, for instance in my (old) edition of Jackson's "Classical Electrodynamics". Note further that the macroscopic field $E$ does not necessarily coincide with the microscopic electric field.   
