Here's a helpful start. Notice that
\begin{align}
(M(x-u))^T M(y-u)
&= (x-u)^T M^T M(y-u).
\end{align}
Let
- $a = (x-u).$ So $a \in \mathbf{R}^d$.
- $b= (y-u).$ So $b \in \mathbf{R}^d$.
- $W = M^T M.$ So $W \in \mathbf{R}^{d\times d}$.
- $X = ab^T.$ So $X \in \mathbf{R}^{d\times d}$.
Then
\begin{align}
(x-u)^T M^T M(y-u)
&= a^T W b \\
&= \sum_{i,j} W_{ij} \, a_i \, b_j \\
&= \sum_{i,j} W_{ij} (a b^T)_{ij} \\
&= \sum_{i,j} W_{ij} \, X_{ij} \\
&= \langle W, X \rangle_{\text{matrix}}
\end{align}
where $\langle \cdot\,, \cdot \rangle_{\text{matrix}}$ is an inner product on $d\times d$ matrices, which here can be thought of as vectors with $d^2$ entries. (See the wikipedia entry on the Frobenius inner product.)
In other words, your original function $I(M,x,y) = \text{sign}[(M(x-u))^T M(y-u)]$ looks just like a perceptron, $f(\mathbf{x};\mathbf{w}) = \text{sign}(\mathbf{w}^T \mathbf{x})$ but where the input vectors $X$ have dimension $d^2$ rather than $d$.
One thing that's unclear in your question is whether "$u$" and "$M$" are part of the inputs to shatter, or whether they are the model parameters.
Assuming they are the model parameters, then the VC-dimension of the function "$I$" is going to be related to the VC-dimension of the perceptron on $d^2$-dimensional inputs. I say "related to" rather than "equal to" because the inputs $X = ab^T$ are constrained to be rank-1 matrices, and therefore not all $d^2$-input patterns are possible. And I'm not sure whether this changes the answer to the original question...