Understanding the (Partial) Converse to Cauchy-Riemann We have that for a function $f$ defined on some open subset $U  \subset \mathbb{C}$ then the following if true:

Suppose $u=\mathrm{Re}(f), v=\mathrm{Im}(f)$ and that all partial derivatives $u_x,u_y,v_x,v_y$ exists and are continuous on $U$. Suppose further that they satisfy the Cauchy-Riemann equations. Then $f$ is holomorphic on $U$.

The proof for this is readily available, though there is a subtlety that I can't understand.
We essentially want to compute $\lim_{h \rightarrow 0} \dfrac{f(z+h)-f(z)}{h}$ where $h=p+qi \in \mathbb{C}$.
We need a relationship like $u(x+p,y+q)-u(x,y)=pu_x(x,y)+qu_y(x,y)+o(|p|+|q|)$. Why is this relationship true?
 A: One way to look at this is by interpreting $f$ as a function from $\mathbb{R}^2 \to \mathbb{R}^2$. This allows you to apply whatever you know about differentiability in $\mathbb{R}^2$ to $\mathbb{C}$. In other words, you interpret $f$ as  function $$
  f \,:\, (u,v) \to \left(\textrm{Re}(f(u+iv)), \textrm{Im}(f(u+iv)\right) =: (f_1(u,v),f_2(u,v)) \text{.}
$$
If all the partial derivatives $\frac{\partial f_1}{\partial u}$, $\frac{\partial f_1}{\partial v}$, $\frac{\partial f_2}{\partial u}$, $\frac{\partial f_2}{\partial v}$ exist and are continuous, $f$ is differentiable in the $\mathbb{R}^2$-sense, i.e. for every $\mathbf{x}=(u,v)$ there's a linear function $df$ which approximates $f$ around $\mathbf{x}$. In other words, $$\begin{aligned}
  f(\mathbf{x} + \Delta) &= f(\mathbf{x}) + df_{\mathbf{x}}\cdot \Delta + R_{\mathbf{x}}(\Delta)
 \quad\text{where}\quad
  \lim_{\Delta\to 0} \frac{\left|\left| R_{\mathbf{x}}(\Delta) \right|\right|}{||\Delta||} = 0
  \quad\text{and}\quad \\
  df_{\mathbf{x}} &= \left(\begin{matrix}
  \frac{\partial f_1}{\partial u} &\frac{\partial f_1}{\partial v} \\
  \frac{\partial f_2}{\partial u} &\frac{\partial f_2}{\partial v}
  \end{matrix}\right) \text{.}
\end{aligned}$$
If the cauchy-riemann differential equations hold, $df_\mathbf{x}$ is then a scaled rotation, i.e. there are $A_\mathbf{x}$ and $B_\mathbf{x}$ such that $$
  df_\mathbf{x} = \left(\begin{matrix}
  A_\mathbf{x}& B_\mathbf{x} \\
  -B_\mathbf{x}& A_\mathbf{x}
  \end{matrix}\right) \text{.}
$$
This form of $df_\mathbf{x}$ allows you to translate back into $\mathbb{C}$, because its exactly these matrices which correspond to complex multiplication. Simply set $$
  f'(z) = f'(u+iv) = A_{(u,v)} + iB_{(u,v)} \text{.}
$$
and translate the approximation back to $\mathbb{C}$, i.e. $$
  f(z+\Delta) = f(z) + f'(z)\cdot\Delta + R_z(\Delta)
  \quad\text{where}\quad
  \lim_{\Delta\to 0} \frac{\left| R_{\mathbf{x}}(\Delta) \right|}{|\Delta|} = 0 \text{,}
$$
which is exactly what you want to show.
A: Denote $h=\pmatrix{p\\q}.$ 
Since $u(x,\ y)$ is differentiable at the point $(x,\ y),$   increment of $u$ can be represented as
 $$u(x+p,\ y+q)-u(x,\ y)=Du\;h+o(\Vert h\Vert)=pu_x(x,\ y)+qu_y(x,\ y)+o(|p|+|q|),$$
where $Du=\left(\dfrac{\partial{u}}{\partial{x}}, \ \dfrac{\partial{u}}{\partial{y}} \right).$
