TL;DR
For an image of size $m\times n$ you can solve this problem in $O(nm(\log(n) + \log(m)))$. In fact, there is nothing to "solve", the solution can be written down analytically. The differencing operator is linear space invariant, i.e. a convolution, and you are asking how to do a "de-convolution". This can be done efficiently in the frequency domain: transform the image with a 2D FFT and divide by the transfer function of the differencing filter, then transform back. There are some details that need to be addressed: the convolution is not cyclic, and the solution you ask for is a least-squares one. This is explained below, where I show the relation to finite-difference approximations of the Poisson equation and Sylvester equations.
Details
To make the notations simpler, let's assume $m = n$. Using the finite-difference matrix
$$D = \begin{bmatrix}1 &-1& & & \\ & 1 & -1 & & \\ & & \ddots & \ddots & \\ & & & 1 & -1\end{bmatrix},$$
we can write your equation as
$$\begin{bmatrix}I \otimes D\\ D^T \otimes I\end{bmatrix} x\approx \begin{bmatrix}b_r\\b_c\end{bmatrix},$$
which is to be solved in the least-squares sense, with $b_r, b_c$ the differences in the row and column directions. The normal equations for the least-squares problem are obtained by multiplying both sides by $\left[I\otimes D^T, \ D\otimes I\right]$, which gives
$$(I\otimes D^TD + DD^T\otimes I)x=(I\otimes D^T) b_r + (D\otimes I)b_c.$$
This can be written as a Sylvester equation (see https://en.wikipedia.org/wiki/Kronecker_product),
$$A X + XB = C,$$
where $A = D^TD$, $B = DD^T$ and $C = DB_r + B_cD$, and the upper-case letters are matrix versions of the lower-case letters representing column vectors. In particular, $X$ is just the original $n\times n$ image. We note the $A$ and $B$ are finite-difference 1D Laplacian matrices, so your problem is reduced to a discrete version of the Poisson equation.
In your question you did not mention any boundary conditions, without which the problem is not well defined. If you only specify the differences between neighboring pixels, you can add an arbitrary constant to your solution without affecting the least-squares error. There are several possible boundary-conditions to choose from, for example, you could use Dirichlet b.c. for all edges. Choosing Neumann b.c. for all edges will not yield a unique solution, so that should be avoided. Whichever b.c. you choose, it will determine the exact form of the $A, \ B$ (and $C$) matrices, which will be slightly different from $DD^T$ and $D^TD$. For the solution to be unique the spectrum of $A$ and $-B$ must be disjoint (https://en.wikipedia.org/wiki/Sylvester_equation).
For all common choices of b.c., the resulting $A$ and $B$ have known eigenvalues and eigenvectors (https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors_of_the_second_derivative), and we may write,
$$A = U_A\Lambda_A U_A^T, \ B = U_B\Lambda_B U_B^T,$$
where $\Lambda_A = \mathrm{diag}(\lambda_A), \Lambda_B = \mathrm{diag}(\lambda_B)$, and $U_A^T U_A = I, U_B^T U_B = I$. Moreover, every choice of b. c. is associated with some version of the discrete cosine or sine transforms for which fast algorithms exist, so that computing the product of $U_A, \ U_B$ times a vector can be done in $O(n \log n)$ (https://en.wikipedia.org/wiki/Discrete_sine_transform).
Inserting the diagonalizations into the Sylvester equation we have,
$$U_A\Lambda_A U_A^T X + X U_B\Lambda_B U_B^T = C \implies \Lambda_A (U_A^T X U_B) + (U_A^TX U_B)\Lambda_B = U_A^T C U_B,$$
which is a Sylvester equation in the transformed image $\hat{X} = U_A^T X U_B$ with diagonal coefficient matrices. Such equations have a known solution,
$$\hat{X}_{ij} = \frac{(U_A^T C U_B)_{ij}}{{(\lambda_A)}_i + {(\lambda_B)}_j}$$
(note the relation to the condition for the existence of the solution to the Sylvester equation). Lastly, transform back to the image to obtain,
$$X = U_A \hat{X} U_B^T.$$
