I think this part in article is starting point of the confusion :
See how the $I(x+u, y+v)$ changed into a totally different form
$I(x,y)+ uI_x + vI_y)$?
The author should have write a bit more detailed notation as:
$I(x,y)+ uI_x(x,y) + vI_y(x,y)$
i.e. pointing out that these are element-wise spatial derivatives (like Sobel) not matrix derivatives ( ...
Check out Szeliski's book:
There is also an old book on feature detection:
You can read the sections that you care about. Also, I think it is always a good idea to read about scale space theory if you are to ...
I think the first smoothing, by $\sigma_D$, is only done to get more stable derivatives whereas in the second step the convolution by a Gaussian with $\sigma_I$ is done to establish the 'scale-space' in which the operator is applied.
Ignoring $\sigma_D$, this looks like this in Matlab:
dx = [-1 0 1; -1 0 1; -1 0 1]; % Simple mask for derivative
Ix = ...