Solving system of equations from partial derivatives & Lagrange multipliers

Questions
,
1

I’m trying to solve a system of equations using the following code:

Solve[(x11/theta11) == lambda1 + (lambda3*pi1) && (x12/theta12) == 

   lambda1 + (lambda3*pi2) && (x13/p1) == lambda3 && (x21/theta21) == 

   lambda2 + (lambda4*pi1) && (x22/theta22) == 

   lambda2 + (lambda4*pi2) && (x23/p2) == 

   lambda4 && (x11/pi1) + (x21/pi1) == (lambda3*theta11) + (lambda4*

      theta21) && (x12/pi2) + (x22/pi2) == (lambda3*

      theta12) + (lambda4*theta22) && theta11 + theta12 == 1 && 

  theta21 + theta22 == 1 && (theta11*pi1) + (theta12*pi2) + p1 == 

   1 && (theta21*pi1) + (theta22*pi2) + p2 == 1, {theta11, theta12, 

  theta21, theta22, pi1, pi2, p1, p2, lambda1, lambda2, lambda3, 

  lambda4}]

The code has been running for over 24 hrs, so I’m assuming there’s something wrong. The system of equations arises by taking the partial derivatives of each theta, pi, p, and lambda variable using the following code:

D[x11*Log[theta11*pi1] + x12*Log[theta12*pi2] + x13*Log[p1] + 

  x21*Log[theta21*pi1] + x22*Log[theta22*pi2] + x23*Log[p2] - 

  lambda1*(theta11 + theta12 - 1) - 

  lambda2*(theta11*pi1 + theta12*pi2 + p1 - 1) - 

  lambda3*(theta21 + theta22 - 1) - 

  lambda4*(theta21*pi1 + theta22*pi2 + p2 - 1), theta11]

for each theta, pi, p, and lambda variable. Essentially, I’m trying to find parameter values that yield the maximum of the function under certain conditions, so I’m taking the partial derivatives with Lagrange multipliers, then solving the system of equations.

Does anyone see any reason why the above Solve[] code would not run properly? Any suggestions on what else I can do to solve this problem? Any help is appreciated. Thanks!

2

Hello and welcome to the Wolfram Community! Please take a few minutes to read this tutorial about correct posting – especially of Mathematica code:

How to type up a post: editor tutorial & general tips

If you will not follow the tutorial, other members of community may not be able to test your code. To edit your post – click “Edit” in the lower right corner of your post. You can also attach a notebook.

3

Several of your variables are linear combinations of other variables. Eliminating some variables may greatly increase the speed of solutions.

This

Reduce[Eliminate[(x11/theta11) == lambda1 + (lambda3*pi1) &&
 (x12/theta12) == lambda1 + (lambda3*pi2) && (x13/p1) == lambda3 &&
 (x21/theta21) == lambda2 + (lambda4*pi1) &&
 (x22/theta22) == lambda2 + (lambda4*pi2) && (x23/p2) == lambda4 &&
 (x11/pi1) + (x21/pi1) == (lambda3*theta11) + (lambda4*theta21) &&
 (x12/pi2) + (x22/pi2) == (lambda3*theta12) + (lambda4*theta22) &&
 theta11 + theta12 == 1 && theta21 + theta22 == 1 &&
 (theta11*pi1) + (theta12*pi2) + p1 == 1 &&
 (theta21*pi1) + (theta22*pi2) + p2 == 1,
 {p1, p2, theta11, theta21, lambda2}],
 {theta12, theta22, pi1, pi2, lambda1, lambda3, lambda4}]

returns in a few seconds a long list of various alternative solutions, depending on values of some parameters.

Perhaps you can look through those alternatives, possibly dismiss many of those, and select your solution.
From that you can then reconstruct the values of the eliminated variables, if necessary.

Alternately, if you know that some of your variables are not zero, for example, you can provide that information
to Reduce and it will automatically eliminate alternatives that do not satisfy these conditions.

4

Thank you, Bill. I appreciate your thoughtful response; however, I’m not sure this will give me the solutions I’m after. I essentially want to define each parameter (all theta’s, pi’s, and p’s) in terms of the x’s only. I think this is possible. I’m struggling to understand if:

  1. my math calculations are incorrect
  2. there is something wrong with my code (I am a very new Mathematica user), or
  3. I’ve reached the limitations of Mathematica

Since Mathematica gives no indication that something is wrong, I’m inclinded to believe it is still working on a solution. However, it has now been running for nearly a week and there is also no indication that anything has been solved. Does anyone have any other ideas about what may be going on?

5

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.

6

Brilliant Daniel! These solutions look reasonable. Thank you so much for your help!