Might work with numerators to make everything strictly polynomial.
exprs = Subtract @@@ eqns;
exprsnum = Numerator[Together]
(* Out[378]= {-lambda1 theta11 - lambda3 pi1 theta11 +
x11, -lambda1 theta12 - lambda3 pi2 theta12 + x12, -lambda3 p1 +
x13, -lambda2 theta21 - lambda4 pi1 theta21 +
x21, -lambda2 theta22 - lambda4 pi2 theta22 + x22, -lambda4 p2 +
x23, -lambda3 pi1 theta11 - lambda4 pi1 theta21 + x11 +
x21, -lambda3 pi2 theta12 - lambda4 pi2 theta22 + x12 + x22, -1 +
theta11 + theta12, -1 + theta21 + theta22, -1 + p1 + pi1 theta11 +
pi2 theta12, -1 + p2 + pi1 theta21 + pi2 theta22} *)
Form a Groebner basis over the rational function field comprised of the xjk's (this seems to be
what is wanted although the post is not entirely clear on this point).
Timing[
gb = GroebnerBasis[
exprsnum, {p1, p2, pi1, pi2, theta11, theta12, theta21,
theta22}, {lambda1, lambda2, lambda3, lambda4},
CoefficientDomain -> RationalFunctions];]
(* Out[395]= {1.744000, Null} *)
Now solve for the variables of interest in terms of those parameters.
Timing[
solns = Solve[
gb == 0, {p1, p2, pi1, pi2, theta11, theta12, theta21, theta22,
lambda1, lambda2, lambda3, lambda4}]]
(* {0.996000, {{p1 -> x13/(x11 + x12 + x13),
p2 -> x23/(x21 + x22 + x23),
pi1 -> -((-x12 x21 + x11 x22)/(
x12 x21 - x11 x22 - x13 x22 + x12 x23)),
pi2 -> -((-x12 x21 + x11 x22)/(
x12 x21 + x13 x21 - x11 x22 - x11 x23)),
theta11 -> -((
x11 x12 x21 - x11^2 x22 - x11 x13 x22 +
x11 x12 x23)/((x11 + x12 + x13) (-x12 x21 + x11 x22))),
theta12 -> -((
x12^2 x21 + x12 x13 x21 - x11 x12 x22 -
x11 x12 x23)/((x11 + x12 + x13) (-x12 x21 + x11 x22))),
theta21 -> -((-x12 x21^2 + x11 x21 x22 + x13 x21 x22 -
x12 x21 x23)/((x12 x21 - x11 x22) (x21 + x22 + x23))),
theta22 -> -((-x12 x21 x22 - x13 x21 x22 + x11 x22^2 +
x11 x22 x23)/((x12 x21 - x11 x22) (x21 + x22 + x23))),
lambda1 -> 0, lambda2 -> 0, lambda3 -> x11 + x12 + x13,
lambda4 ->
x21 + x22 + x23}, {p1 -> -((-x13 - x23)/(
x11 + x12 + x13 + x21 + x22 + x23)),
p2 -> -((-x13 - x23)/(x11 + x12 + x13 + x21 + x22 + x23)),
pi1 -> ((x11 + x21) (x13 x21 + x13 x22 - x11 x23 -
x12 x23))/((x11 + x12 + x13 + x21 + x22 + x23) (x13 x21 -
x11 x23)),
pi2 -> ((x12 + x22) (x13 x21 + x13 x22 - x11 x23 -
x12 x23))/((x11 + x12 + x13 + x21 + x22 + x23) (x13 x22 -
x12 x23)),
theta11 -> -((
x13 x21 - x11 x23)/(-x13 x21 - x13 x22 + x11 x23 + x12 x23)),
theta12 -> -((
x13 x22 - x12 x23)/(-x13 x21 - x13 x22 + x11 x23 + x12 x23)),
theta21 -> -((
x13 x21 - x11 x23)/(-x13 x21 - x13 x22 + x11 x23 + x12 x23)),
theta22 -> -((
x13 x22 - x12 x23)/(-x13 x21 - x13 x22 + x11 x23 + x12 x23)),
lambda1 -> -((x13 x21 + x13 x22 - x11 x23 - x12 x23)/(x13 + x23)),
lambda2 -> -((-x13 x21 - x13 x22 + x11 x23 + x12 x23)/(x13 + x23)),
lambda3 -> -((-x11 x13 - x12 x13 - x13^2 - x13 x21 - x13 x22 -
x13 x23)/(x13 + x23)),
lambda4 -> -((-x11 x23 - x12 x23 - x13 x23 - x21 x23 - x22 x23 -
x23^2)/(x13 + x23))}}} *)
Since we extracted numeraots, backsubstitute to make sure no denominators get zeroed.
Together[exprs /. solns]
(* Out[402]= {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}} *)
Looks fine.
