f = 2/((x1 + Sqrt[3] x2)^2 + (y1 + Sqrt[3] y2)^2)
X = vx1^2*D[D[f, x1], x1] + vx2^2*D[D[f, x2], x2] + 2 vx1*vx2*D[D[f, x1], x2]
Y= vx1*D[f, x1] + vx2*D[f, x2]
g=X+Y^2

H = ImplicitRegion[{3 (-x1^2 - y1^2 + x2^2 + y2^2) +
 2 Sqrt[3] (x1*x2 + y1*y2) == 0}, {x1, y1, x2, y2}]

FindMaximum[{g, {x1, y1, x2, y2} ∈ H}, {x1, y1, x2, y2}]

reg = ImplicitRegion[0 <= x <= 2, x]
f = a*x^2 + b
g = MaxValue[f, Element[{x}, reg]]
g /. {a -> 1, b -> 2}

f[x1_, x2_, y1_, y2_] :=
 2/((x1 + Sqrt[3] x2)^2 + (y1 + Sqrt[3] y2)^2)

X[vx1_, vx2_, x1_, x2_, y1_, y2_] :=
 vx1^2*D[D[f[x1, x2, y1, y2], x1], x1] +
  vx2^2*D[D[f[x1, x2, y1, y2], x2], x2] +
  2 vx1*vx2*D[D[f[x1, x2, y1, y2], x1], x2]

Y[vx1_, vx2_, x1_, x2_, y1_, y2_] :=
 vx1*D[f[x1, x2, y1, y2], x1] + vx2*D[f[x1, x2, y1, y2], x2]

g[vx1_, vx2_, x1_, x2_, y1_, y2_] :=
 X[vx1, vx2, x1, x2, y1, y2] + Y[vx1, vx2, x1, x2, y1, y2]^2

H = ImplicitRegion[{3 (-x1^2 - y1^2 + x2^2 + y2^2) + 2 Sqrt[3] (x1*x2 + y1*y2) == 0}, {x1, y1, x2, y2}]

max[vx1_, vx2_] :=
 FindMaximum[{g[vx1, vx2, x1, x2, y1,
    y2], {x1, y1, x2, y2} ∈ H}, {x1, y1, x2, y2}]

m = max[1000., 1000.]

m[[2, All, 2]] ∈ H