System of Equations Question

[Ben Goodrich]

Hi,

I maintain an R package, but there is one place where a symbolic
solution is needed to verify a result. I would like to write an R
function that prints proper SAGE input so that users can easily feed
it to SAGE. However, I have not yet been able to get SAGE to produce
the expected behavior for a relatively simple published example and
have run out of ideas.

The example (for anyone interested) is from

http://scholar.google.com/scholar?num=100&hl=en&lr=&q=allintitle%3A%22Uniqueness+does+not+imply+identification%22&btnG=Search

and involves a system of 10 equations and 10 unknowns. But some of the
equations are redundant, and there are multiple solutions. I was
expecting SAGE to give me a symbolic solution where at least one
unknown was expressed in terms of another unknown, but instead it says
empty set. If there were a unique solution, I would want SAGE to
express lambda_** and Theta2_** in terms of Sigma_**

Here is the input that causes [] to be outputted:

# Define variables
lambda_11, lambda_21, lambda_31, lambda_41, lambda_12, lambda_22,
lambda_32, lambda_42, Phi_11, Phi_21, Phi_22, Theta2_11, Theta2_22,
Theta2_33, Theta2_44, Sigma_11, Sigma_21, Sigma_31, Sigma_41,
Sigma_22, Sigma_32, Sigma_42, Sigma_33, Sigma_43, Sigma_44 =
var('lambda_11 lambda_21 lambda_31 lambda_41 lambda_12 lambda_22
lambda_32 lambda_42 Phi_11 Phi_21 Phi_22 Theta2_11 Theta2_22 Theta2_33
Theta2_44 Sigma_11 Sigma_21 Sigma_31 Sigma_41 Sigma_22 Sigma_32
Sigma_42 Sigma_33 Sigma_43 Sigma_44')

# Specify constraints
lambda_12 = 0/1
lambda_31 = 0/1
lambda_41 = 0/1
Phi_11 = 1/1
Phi_22 = 1/1

# Write out 10 equations from the matrix algebra
eq1 = Sigma_11 == ( (lambda_11 * Phi_11 + lambda_12 * Phi_21) *
lambda_11 + (lambda_11 * Phi_21 + lambda_12 * Phi_22) * lambda_12 ) +
Theta2_11

eq2 = Sigma_21 == ( (lambda_21 * Phi_11 + lambda_22 * Phi_21) *
lambda_11 + (lambda_21 * Phi_21 + lambda_22 * Phi_22) * lambda_12 )

eq3 = Sigma_31 == ( (lambda_31 * Phi_11 + lambda_32 * Phi_21) *
lambda_11 + (lambda_31 * Phi_21 + lambda_32 * Phi_22) * lambda_12 )

eq4 = Sigma_41 == ( (lambda_41 * Phi_11 + lambda_42 * Phi_21) *
lambda_11 + (lambda_41 * Phi_21 + lambda_42 * Phi_22) * lambda_12 )

eq5 = Sigma_22 == ( (lambda_21 * Phi_11 + lambda_22 * Phi_21) *
lambda_21 + (lambda_21 * Phi_21 + lambda_22 * Phi_22) * lambda_22 ) +
Theta2_22

eq6 = Sigma_32 == ( (lambda_31 * Phi_11 + lambda_32 * Phi_21) *
lambda_21 + (lambda_31 * Phi_21 + lambda_32 * Phi_22) * lambda_22 )

eq7 = Sigma_42 == ( (lambda_41 * Phi_11 + lambda_42 * Phi_21) *
lambda_21 + (lambda_41 * Phi_21 + lambda_42 * Phi_22) * lambda_22 )

eq8 = Sigma_33 == ( (lambda_31 * Phi_11 + lambda_32 * Phi_21) *
lambda_31 + (lambda_31 * Phi_21 + lambda_32 * Phi_22) * lambda_32 ) +
Theta2_33

eq9 = Sigma_43 == ( (lambda_41 * Phi_11 + lambda_42 * Phi_21) *
lambda_31 + (lambda_41 * Phi_21 + lambda_42 * Phi_22) * lambda_32 )

eq10 = Sigma_44 == ( (lambda_41 * Phi_11 + lambda_42 * Phi_21) *
lambda_41 + (lambda_41 * Phi_21 + lambda_42 * Phi_22) * lambda_42 ) +
Theta2_44

# Try to solve for 10 unknowns
solutions = solve([eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10],
lambda_11, lambda_21, lambda_22, lambda_32, lambda_42, Phi_21,
Theta2_11, Theta2_22, Theta2_33, Theta2_44)

solutions # empty set


But if you fill in the values from a known solution, then eq1 through
eq10 all evaluate to True

# Solution from Table 1 in paper
lambda_11 = 1/2
lambda_21 = 1/1
lambda_22 = 1/5
lambda_32 = 2/1
lambda_42 = 3/2

Phi_21 = 4/5

Theta2_11 = 1/1
Theta2_22 = 1/2
Theta2_33 = 3/2
Theta2_44 = 2/1

Sigma_11 = 5/4
Sigma_21 = 58/100
Sigma_31 = 4/5
Sigma_41 = 3/5
Sigma_22 = 186/100
Sigma_32 = 2/1
Sigma_42 = 3/2
Sigma_33 = 55/10
Sigma_43 = 3/1
Sigma_44 = 425/100

print(eq1); print(eq2); print(eq3); print(eq4); print(eq5);
print(eq6); print(eq7); print(eq8); print(eq9); print(eq10) # all True

Thus, I do not understand why SAGE says there is no solution to this
system of equations. What should I be doing differently?

Thanks in advance,
Ben

[David Joyner]

I believe Sage simply calls Maxima for the solution. Since you
obviously know the
most about the problem, perhaps the easiest thing to do would be to determine
that it is Sage and not Maxima that is at fault. Perhaps you could see if
the solution is obtained in Maxima? (On the command line of Sage, you
can simply type "maxima_colsole()"; in the notebook I believe there is
a drop down
menu for maxima mode.)

[Ben Goodrich]

On Feb 18, 3:14 pm, Ben Goodrich <goodrich....@gmail.com> wrote:
> I maintain an R package, but there is one place where a symbolic
> solution is needed to verify a result. I would like to write an R
> function that prints proper SAGE input so that users can easily feed
> it to SAGE. However, I have not yet been able to get SAGE to produce
> the expected behavior for a relatively simple published example and
> have run out of ideas.

To (partially) answer my own question, it seems that it is necessary
to specify that the equations be solved for Sigma_** as well (but I
don't know why). E.g.

# Substitute this line

solutions = solve([eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10],
lambda_11, lambda_21, lambda_22, lambda_32, lambda_42, Phi_21,

Theta2_11, Theta2_22, Theta2_33, Theta2_44,
Sigma_11, Sigma_21, Sigma_31, Sigma_41, Sigma_22, Sigma_32, Sigma_42,

Sigma_33, Sigma_43, Sigma_44, solution_dict = True, explicit_solution
= True); solutions

which yields

[{Sigma_41: r31*r35*r36, Sigma_42: r32*r35*r36 + r33*r35, lambda_42:
r35, Theta2_22: r38, Phi_21: r36, Theta2_33: r39, Sigma_44: r40 +
r35^2,
Sigma_43: r34*r35, Sigma_22: r38 + 2*r32*r33*r36 + r33^2 + r32^2,
lambda_22: r33, lambda_21: r32, Sigma_21: r31*r33*r36 + r31*r32,
Sigma_31: r31*r34*r36, Sigma_33: r39 + r34^2, lambda_32: r34,
Sigma_32:
r32*r34*r36 + r33*r34, Theta2_44: r40, Theta2_11: r37, Sigma_11: r37 +
r31^2, lambda_11: r31}]

But is there a way to get SAGE to express the answer in terms of the
original variables instead of these terms involving r? Is this SAGE's
way of telling me (correctly) that there is no solution where
lambda_** and Theta2_** can be expressed solely in terms of Sigma_** ?

Thanks,
Ben

[Ben Goodrich]

On Feb 18, 3:58 pm, "David Joyner" <wdjoy...@gmail.com> wrote:
> I believe Sage simply calls Maxima for the solution. Since you
> obviously know the
> most about the problem, perhaps the easiest thing to do would be to determine
> that it is Sage and not Maxima that is at fault. Perhaps you could see if
> the solution is obtained in Maxima? (On the command line of Sage, you
> can simply type "maxima_colsole()"; in the notebook I believe there is
> a drop down
> menu for maxima mode.)
>

I tried running it via maxima mode under the Online Notebook, but it
ran for about an hour without producing any output. In SAGE mode, it
(now) produces something in a few seconds. So, I think I would be okay
if I knew how to interpret the output. For example, what does the r35
mean in "lambda_42: r35"? This syntax does not seem to appear in the
documentation.

[David Joyner]

On Feb 18, 2008 8:29 PM, Ben Goodrich <goodri...@gmail.com> wrote:
>
> On Feb 18, 3:58 pm, "David Joyner" <wdjoy...@gmail.com> wrote:
> > I believe Sage simply calls Maxima for the solution. Since you
> > obviously know the
> > most about the problem, perhaps the easiest thing to do would be to determine
> > that it is Sage and not Maxima that is at fault. Perhaps you could see if
> > the solution is obtained in Maxima? (On the command line of Sage, you
> > can simply type "maxima_colsole()"; in the notebook I believe there is
> > a drop down
> > menu for maxima mode.)
> >
>
> I tried running it via maxima mode under the Online Notebook, but it
> ran for about an hour without producing any output. In SAGE mode, it

To me that suggests a Maxima bug.

> (now) produces something in a few seconds. So, I think I would be okay

Something *correct*?

> if I knew how to interpret the output. For example, what does the r35
> mean in "lambda_42: r35"? This syntax does not seem to appear in the

I don't see how that contains debugging info. I thin it means it is the
35th R-expression that Sage has had to keep track of to run that
particular session.

> documentation.
>
> >
>

[William Stein]

No, Maxima creates variables r[number] for each free
parameter in the solution set it finds to a problem, e.g.,:

sage: var('x,y')
(x, y)
sage: solve(3*x+5*y+7==0,[x,y])
[[x == (-5*r1 - 7)/3, y == r1]]
sage: var('x,y,z')
(x, y, z)
sage: solve(3*x+5*y+7*z+13==0,[x,y])
[[x == (-7*z - 5*r2 - 13)/3, y == r2]]

[David Joyner]

On Feb 18, 2008 9:09 PM, William Stein <wst...@gmail.com> wrote:
>
> On Feb 18, 2008 6:04 PM, David Joyner <wdjo...@gmail.com> wrote:
> >
> > On Feb 18, 2008 8:29 PM, Ben Goodrich <goodri...@gmail.com> wrote:
> > >
> > > On Feb 18, 3:58 pm, "David Joyner" <wdjoy...@gmail.com> wrote:
> > > > I believe Sage simply calls Maxima for the solution. Since you
> > > > obviously know the
> > > > most about the problem, perhaps the easiest thing to do would be to determine
> > > > that it is Sage and not Maxima that is at fault. Perhaps you could see if
> > > > the solution is obtained in Maxima? (On the command line of Sage, you
> > > > can simply type "maxima_colsole()"; in the notebook I believe there is
> > > > a drop down
> > > > menu for maxima mode.)
> > > >
> > >
> > > I tried running it via maxima mode under the Online Notebook, but it
> > > ran for about an hour without producing any output. In SAGE mode, it
> >
> > To me that suggests a Maxima bug.
> >
> > > (now) produces something in a few seconds. So, I think I would be okay
> >
> > Something *correct*?
> >
> > > if I knew how to interpret the output. For example, what does the r35
> > > mean in "lambda_42: r35"? This syntax does not seem to appear in the
> >
> > I don't see how that contains debugging info. I thin it means it is the
> > 35th R-expression that Sage has had to keep track of to run that
> > particular session.
> >
>
> No, Maxima creates variables r[number] for each free

What I said makes no sense - R isn't even running. Time to quit
for the night!

[Ben Goodrich]

On Feb 18, 9:09 pm, "William Stein" <wst...@gmail.com> wrote:
> No,  Maxima creates variables r[number] for each free
> parameter in the solution set it finds to a problem, e.g.,:
>
> sage: var('x,y')
> (x, y)
> sage: solve(3*x+5*y+7==0,[x,y])
> [[x == (-5*r1 - 7)/3, y == r1]]
> sage: var('x,y,z')
> (x, y, z)
> sage: solve(3*x+5*y+7*z+13==0,[x,y])
> [[x == (-7*z - 5*r2 - 13)/3, y == r2]]

Thank you for this clarification. Maybe my questions pertain to the
underlying Maxima functions, but I am still a bit confused on some
things.

For example
sage: var('x,y,z')
sage: solve([3*x+5*y+7*z+13==0],x)
[x == (-7*z - 5*y - 13)/3]

expresses x in terms of z and y, instead of r[numbers]. That was
basically the behavior I was trying to achieve in the first post, to
express lambda_** and Theta2_** in terms of Sigma_** . But that
resulted in a purported empty set, unless I also requested that it
solve for Sigma_** (third post), which resulted in a solution with
r[numbers]. Is it the case that the r[numbers] will only appear if
there is no unique solution to the system of equations?

In the same vein,

sage: var('x,y,z')
sage: solve([3*x+5*y+7*z+13==0, x == y, z == x + y + 1], x,y) #
error message
sage: solve([3*x+5*y+7*z+13==0, x == y, z == x + y + 1], x,y,z) #
correct solution

sage: solve([3*x+5*y+7*z+13==0, x == y], x,y) # solves for x and y
in terms of z
sage: solve([3*x+5*y+7*z+13==0, x == y], x,y,z) # solution with
r[numbers]

The r[numbers] are not so bad in a simple example where you can just
eyeball the substitutions, but it would be nice to have the solution
expressed in the original variables in big systems of equations. That
is what I was trying to achieve by not solving for all the variables,
but it seems that the behavior of the function is shaky unless all
unknowns are solved for.

Thanks,
Ben