Clear[Evaluate[Context[]<>"*"]];
arg2d[r_]:=-(xx r[[1]])^2-(yy r[[2]])^2;
g2d[x_,y_]:=Exp[e+arg2d[RotationMatrix[α].{x,y}]];
Print["arg = ",Collect[FullSimplify[Log@g2d[x,y],{x>0,y>0,α>0,e>0,xx>0,yy>0}],{x,y}]];
grad=Map[D[(vEx-of-g2d[x-x0,y-y0])^2,#]&,{of,e,xx,yy,x0,y0,α}];
grad=grad/.{g2d[x-x0,y-y0]->vFit};
grad=grad/(vEx-vFit-of)/-2;
Print["d/dof = ",First@grad];
grad=Rest@grad/vFit;
Print[MatrixForm[grad[[1;;3]]]];
grad=Rest@Rest@Rest@grad;
grad=FullSimplify@grad;
grad=grad/.{
xx^2Cos[α]^2+yy^2Sin[α]^2->-xxs,
yy^2Cos[α]^2->-yys-xx^2Sin[α]^2
};
grad=FullSimplify@grad;
grad=grad/.{Sin[2α]->xys/(xx^2-yy^2)};
grad=FullSimplify@grad;
MatrixForm@Collect[grad,{xxs,yys,xys}]
Map[D[(v-v1-v2-v3)^2,#]&,{v1,v2,v3}]

Clear[Evaluate[Context[]<>"*"]];
n=2;
parameter=Join[
{of},
Flatten@Table[
{e[i],xx[i],yy[i],x0[i],y0[i],α[i]},
{i,1,n}
]
];
dxs[i_]:=x-x0[i];
dys[i_]:=y-y0[i];
xxs[i_]:=-(xx[i]Cos[α[i]])^2-(yy[i]Sin[α[i]])^2;
yys[i_]:=-(yy[i]Cos[α[i]])^2-(xx[i]Sin[α[i]])^2;
xys[i_]:=(xx[i]^2-yy[i]^2)2Cos[α[i]]Sin[α[i]];
gausz[i_]:=Exp[e[i]+xxs[i]dxs[i]^2+yys[i]dys[i]^2+xys[i]dxs[i]dys[i]];
tmp=of-val+Sum[gausz[i],{i,1,n}];
abstand=tmp^2;
gradAnalytisch=Map[D[abstand,#]&,parameter]/2;
gradNumerisch=Join[
{tmp},
Flatten@Table[
tmp*gausz[i]*{
1,
-2xx[i](dxs[i]Cos[α[i]]-dys[i]Sin[α[i]])^2,
-2yy[i](dys[i]Cos[α[i]]+dxs[i]Sin[α[i]])^2,
-2dxs[i]xxs[i]-xys[i]dys[i],
-2dys[i]yys[i]-xys[i]dxs[i],
2(xx[i]^2-yy[i]^2)(dys[i]Cos[α[i]]+dxs[i]Sin[α[i]])(dxs[i]Cos[α[i]]-dys[i]Sin[α[i]])
},
{i,1,n}
]
];
Map[#[[1]]/#[[2]]&,Transpose@{gradAnalytisch,gradNumerisch}]//FullSimplify