The Gradient Operator of a Vectorized Image in Matrix Form
It is pretty simple to create those Matrices.
The real issue with them is their size which is enormous for real world images.

For small kernels they are sparse which saves the day.
Indeed for the Derivative Operator, which has only 2 elements, they are highly sparse.

I built them in MATLAB using:

mI = im2double(imread(imageFileName));
mI = mI(11:410, 201:600, 1);

% mI = mI(11:20, 201:210, 1);

numRows     = size(mI, 1);
numCols     = size(mI, 2);
numPixels   = numRows * numCols;

mO = zeros([numRows, numCols, length(vParamLambda)]); %&lt;! Output

vDx = [1, -1]; %&lt;! Matrix is doing Correlation, this for Convolution
vDy = [1; -1];

% Sanity Check
mIxRef = conv2(mI, vDx, 'valid');
mIyRef = conv2(mI, vDy, 'valid');

% mIx = reshape(mDx * mI(:), numRows, numCols - 1);
% mIy = reshape(mDy * mI(:), numRows - 1, numCols);

mDh = sparse(numPixels - numCols, numPixels);
mDv = sparse(numPixels - numRows, numPixels);

% mDx = zeros(numPixels - numCols, numPixels);
% mDy = zeros(numPixels - numRows, numPixels);

tic();
colShift = 0;
for ii = 1:(numPixels - numRows)
    if(mod(ii + colShift, numRows) == 0)
        colShift = colShift + 1;
    end
    mDv(ii, ii + colShift) = -1;
    mDv(ii, ii + colShift + 1) = 1;
end

for ii = 1:(numPixels - numCols)
    mDh(ii, ii) = -1;
    mDh(ii, ii + numCols) = 1;
end
toc();


mIx = reshape(mDh * mI(:), numRows, numCols - 1);
mIy = reshape(mDv * mI(:), numRows - 1, numCols);

mE = mIy - mIyRef;
disp(['Maximum Absolute Error Between Matrix Form to Convolution (Vertical Derivative) - ', num2str(max(abs(mE(:))))]);

mE = mIx - mIxRef;
disp(['Maximum Absolute Error Between Matrix Form to Convolution (Horizontal Derivative) - ', num2str(max(abs(mE(:))))]);


If you load image into mI you see the result is identical to the convolution operator.  

The full MATLAB code which implements the solution is at my Signal Processing StackExchange Question 50329 - GitHub Repository.  

Remark
In practice usually those kind of equation are solved using Preconditioned Conjugate Gradient Method.
The nice trick is that method requires only the result of the operators applied on a vectorized image. This can be done using the convolution operator itself instead of using the matrix form.