Harris Corners Detector squared derivatives
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 ( those x and y indicate horizontal and vertical spatial derivative).