Variance of sum of square differences corrupted by noise I've been trying to extend a proof I have for a univariate function $f(x)$ to a
multivariate function $f(x,y)$ but the result I get looks weird. I'm wondering
whether there is a mistake in any step I'm making.
The background is that $f(x,y)$ is an image that is superimposed with some
gaussian noise $n\sim\mathcal{N}(0,\sigma^2)$. This image is shifted by a true
but unknown offset $d_x$ and $d_y$ and the goal is to find an estimate $(\Delta
x,\Delta y)^\mathsf{T}$ for this offset and to estimate its variance. This leads
to the following expression that is to be minimized:
$$
e(\Delta x,\Delta y) = \iint \left(f(x-d_x+\Delta x+\lambda, y-d_y+\Delta y+\eta) - f(x+\lambda,y+\eta) + n_1(x+\lambda,y+\eta) - n_0(x+\lambda,y+\eta)\right)^2\text{d}\lambda\text{d}\eta
$$
I am computing the second order taylor expansion around $d_x=\Delta x$ and
$d_y=\Delta y$. This yields
$$
a(\Delta x - d_x)^2 + b(\Delta y - d_y)^2 +\\ c(\Delta x - d_x)(\Delta y - d_y) + d(\Delta x - d_x) + e(\Delta y - d_y) + f
$$
where
\begin{eqnarray*}
a \approx \iint f_{xx}(x+\lambda,y+\eta)^2\text{d}\lambda\text{d}\eta\\
b \approx \iint f_{yy}(x+\lambda,y+\eta)^2\text{d}\lambda\text{d}\eta\\
c \approx \iint 2 f_{x}(x+\lambda,y+\eta)f_y(x+\lambda,y+\eta)\text{d}\lambda\text{d}\eta\\
d = \iint 2(n_1(x+\lambda,y+\eta)-n_0(x+\lambda,y+\eta)f_x(x+\lambda,y+\eta)\text{d}\lambda\text{d}\eta\\
e = \iint 2(n_1(x+\lambda,y+\eta)-n_0(x+\lambda,y+\eta)f_y(x+\lambda,y+\eta)\text{d}\lambda\text{d}\eta\\
f = \iint (n_1(x+\lambda,y+\eta)-n_0(x+\lambda,y+\eta)^2 \text{d}\lambda\text{d}\eta
\end{eqnarray*}
Subscripts denote partial derivatives. The second order partial derivative terms  have been omitted in the formula above. This is the only step I don't feel very comfortable with, but it is also done in the original, univariate proof. Note that $a$, $b$ and $c$ are constants and only $d$, $e$ and $f$ are random variables. The vertex of this paraboloid -- if it exists -- is found by the equating the partial derivatives of the paraboloid to 0, which results in
$$
\Delta x - d_x = -\frac{2bd-ce}{4ab-c^2},\;
\Delta y - d_y = -\frac{2ae-cd}{4ab-c^2}\;.
$$
I am interested in the variance of the estimate,
$$
\text{Var}\left(\begin{pmatrix}\Delta x\\\Delta y\end{pmatrix}\right) = \text{Var}\left(\begin{pmatrix}d_x\\d_y\end{pmatrix} - \begin{pmatrix}\frac{2bd-ce}{4ab-c^2}\\\frac{2ae-cd}{4ab-c^2}\end{pmatrix}\right) = 
\text{Var}\left(-\begin{pmatrix}\frac{2bd-ce}{4ab-c^2}\\\frac{2ae-cd}{4ab-c^2}\end{pmatrix}\right)
$$
In the next few steps, I will use $\text{Var}(AX)=A\text{Var}(X)A^\mathsf{T}$ repeatedly. First I pull out the denominator:
$$
\text{Var}\left(\begin{pmatrix}\Delta x\\\Delta y\end{pmatrix}\right) = 
\frac{1}{\left(4ab-c^2\right)^2}\text{Var}\left(\begin{pmatrix}2bd-ce\\2ae-cd\end{pmatrix}\right)
$$
Then I put back the definitions of the polynomial coefficients. I'm using a short-hand notation here where I omit the function arguments and abbreviate $n_1(x+\lambda,y+\eta)-n_0(x+\lambda,y+\eta)$ with $\delta n$.
$$
\text{Var}\left(\begin{pmatrix}\Delta x\\\Delta y\end{pmatrix}\right) = 
\frac{1}{\left(4\iint f_{xx}^2\iint f_{yy}^2 - 4(\iint f_{xy})^2\right)^2}\text{Var}\left(\begin{pmatrix}4\iint f_{yy}^2\iint f_{x}\delta n-4\iint f_{xy}\iint f_{y}\delta_n\\4\iint f_{xx}^2\iint f_{y}\delta n-4\iint f_{xy}\iint f_{x}\delta n\end{pmatrix}\right)
$$
Now I decompose the matrix and pull it out of the variance. 
$$
\text{Var}\left(\begin{pmatrix}\Delta x\\\Delta y\end{pmatrix}\right) = 
\frac{1}{\left(4\iint f_{xx}^2\iint f_{yy}^2 - 4(\iint f_{xy})^2\right)^2}\text{Var}\left(4\begin{pmatrix}\iint f_{yy}^2 & -\iint f_xf_y \\ -\iint f_xf_y & \iint f_{xx}^2\end{pmatrix}\begin{pmatrix}\iint f_x\delta_n\\\iint f_y\delta_n\end{pmatrix}\right)=
\frac{16}{\left(4\iint f_{xx}^2\iint f_{yy}^2 - 4(\iint f_{xy})^2\right)^2}\begin{pmatrix}\iint f_{yy}^2 & -\iint f_xf_y \\ -\iint f_xf_y & \iint f_{xx}^2\end{pmatrix}\text{Var}\left(\begin{pmatrix}\iint f_x\delta_n\\\iint f_y\delta_n\end{pmatrix}\right)\begin{pmatrix}\iint f_{yy}^2 & -\iint f_xf_y \\ -\iint f_xf_y & \iint f_{xx}^2\end{pmatrix}
$$
Then I proceed similarly for the remaining integral:
$$
\text{Var}\left(\begin{pmatrix}\Delta x\\\Delta y\end{pmatrix}\right) = 
\frac{32\sigma^2}{\left(4\iint f_{xx}^2\iint f_{yy}^2 - 4(\iint f_{xy})^2\right)^2}\begin{pmatrix}\iint f_{yy}^2 & -\iint f_xf_y \\ -\iint f_xf_y & \iint f_{xx}^2\end{pmatrix}\begin{pmatrix}\iint f_x\\\iint f_y\end{pmatrix}\begin{pmatrix}\iint f_x\\\iint f_y\end{pmatrix}^\mathsf{T}\begin{pmatrix}\iint f_{yy}^2 & -\iint f_xf_y \\ -\iint f_xf_y & \iint f_{xx}^2\end{pmatrix}
$$
Now here is the problem: I would expect the solution
$$
\text{Var}\left(\begin{pmatrix}\Delta x\\\Delta y\end{pmatrix}\right) = 
\frac{k\sigma^2}{\iint f_{xx}^2\iint f_{yy}^2 - (\iint f_{xy})^2}\begin{pmatrix}\iint f_{yy}^2 & -\iint f_xf_y \\ -\iint f_xf_y & \iint f_{xx}^2\end{pmatrix}
$$
where $k$ is some constant. This could be achieved if the square in the denominator would cancel out with the first three matrices. However it does not, I always get these cross-terms. Hence I am wondering whether there is any mistake in my derivation, or alternatively if I could use some approximation that escapes me to get the desired result.
Thank you very much!
 A: I am not trully following the content of the question, but the three matrices in the last expression remind me of the expression for the asymptotic variance of the maximum likelihood estimator:  
Having a sample $\mathbf x$ from $n$ i.i.d, random variables $X_i$, characterized by an unknown parameter vector $\theta_0$, let $H(x_i;\theta_0)$ be the Hessian of the likelihood (for a typical observation, not the likelihood of the sample), and let $S(x_i;\theta_0)$ be its gradient (both evaluated at the true value of the parameter). Then the quantity $\sqrt n (\hat \theta_{MLE}-\theta_0)$ is asymptotically normal with variance
$$\operatorname{Avar}(\hat \theta_{MLE})=\Big[E(H(x_i;\theta_0)\Big]^{-1}E\left(S(x_i;\theta_0)S(x_i;\theta_0)^T\right)\Big[E(H(x_i;\theta_0)\Big]^{-1}$$
Note that as usual we confuse people by saying that this is the variance of $\hat \theta_{MLE}$ while in reality is the variance of the quantity $\sqrt n (\hat \theta_{MLE}-\theta_0)$.
These expected values are consistently estimated by using the corresponding sample means.
Now the information matrix equality (which holds when the model is correctly specified) states that 
$$-\Big[E(H(x_i;\theta_0)\Big]^{-1} = E\left(S(x_i;\theta_0)S(x_i;\theta_0)^T\right)$$
and so things cancel out and we end up with 
$$\operatorname{Avar}(\hat \theta_{MLE})=-\Big[E(H(x_i;\theta_0)\Big]^{-1} = E\left(S(x_i;\theta_0)S(x_i;\theta_0)^T\right)$$
..which is also analogous to the simplification that you expected/hoped, for your variance expression.
