Solution of small system of equations grows past manageability

[janosc...@gmail.com]

I am trying to solve a system of equations with sympy that arises from a constraint of the form:

  (A x B) x C = D

where

* A, B, C and D are 3x3 matrices
* the diagonal of D should be zero
* B is a "rigid motion 2D" transformation, with elements cos(phi), +-sin(phi), x and y
* A and C are fully filled with (supposedly known) values
* I want to solve for phi, x and y

This gives me four equations:

* one for each diagonal element in D
* one additional (quadratic) equation sin^2(phi) + cos^2(phi) = 1

When feeding those to equations directly to sympy, this takes some hours and then breaks with an out of memory message.

My next approach was to help sympy by guiding the solution step by step (*).

* First i took two of the linear equations and let sympy solve for x and y (works great)
* Instead of having cos(phi) and sin(phi) in the B matrix, i introduced new symbols cosphi and sinphi
* Then i took the resulting expressions for x and y, and solve with the third linear equation for the cosphi element (works too)
* Finally i tried to solve the quadratic equation for sinphi by inserting the just gathered cosphi expression
* The last step was not feasible without transforming the expression to a polynom in sinphi and by replacing all coefficient expressions by new symbols, then it worked

The resulting expressions for x, y and phi (written as python expressions) are about 3 MB (!) of text.

This does not seem to be adequate to the problem, and when converting to a theano function i get "maximum recursion depth exceeded".
When i look at the expressions they are very repetitive, so i tried CSE, which brings it down to about 30 KB, but they are still very repetitive and full of patterns.

I suspect that the resulting expressions actually just perform some matrix operations, so probably there would be an efficient way to compute the solution if only one could get back to matrix expressions.
I tried to guess what the appropriate matrix operations are, but without success (**). And this feels of course very wrong and backwards.

Is there some obvious approach to such problems that i missed? Is the problem actually that hard?

I am aiming for a mostly automated solution process without steps like (*) and (**), because i have a hand full of very similar problems ahead ...
Any hint appreciated!

--
Best regards
Janosch

[Jason Moore]

Can you please share the code so we can see what you are doing?

Jason
moorepants.info
+01 530-601-9791

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at https://groups.google.com/group/sympy.
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/8555837f-d87d-484c-b882-2d8f7085d3b2%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[janosc...@gmail.com]

My description was a little compressed, so i had to clean up the code to match my description again ...
The code is available here: http://pastebin.com/MMW3B88h
I hope its readable for you.

[Alan Bromborsky]

Physically what are all the matrices.  Do A and C also describe rotations.  Please give the actual physics problem as well as the resulting math.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/ccde273d-ee0e-46bc-a822-7d1b9ff88f7e%40googlegroups.com.

[janosc...@gmail.com]

Physically, the rows of A are three points fixed on a movable board.

These points run freely in three linear bearings which are placed on a fixed base.

The linear bearings are described in hesse normal form in the rows of matrix C.

The robust motion matrix B is the transformation which transforms points on the board to points in the base.

So together my constraint D = (A * B) * C means
- Transform the points in A from the board to the base: A * B
- Compute the distance from the linear bearings: * C
- Claim that the distances are zero and solve for the motion

I am aware that there are some other approaches to tackle this problem, but i was not able to get a grip on them such that i could formulate them in code.

[janosc...@gmail.com]

I guess its hard to get from my description, so i uploaded a drawing to visualize the physical problem: http://pasteboard.co/1Bvt53hY.png

Thanks for your interest!

[Jason Moore]

Where are phi, x, y on the diagram?

Jason
moorepants.info
+01 530-601-9791

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/f68853ca-2a1f-4ef1-a371-340602cbd63d%40googlegroups.com.

[janosc...@gmail.com]

They describe the location of the board (the blue rectangle) in relation to its "normal" position by a rotation about an angle of phi and a translation of x and y.

[Jason Moore]

If the blue dots are fixed on the board, doesn't the linear bearings remove all degrees of freedom? I don't see how this thing can move.

Jason
moorepants.info
+01 530-601-9791

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/16628e82-8a37-4ce4-8840-e6bcabcb7d20%40googlegroups.com.

[janosc...@gmail.com]

Yes, exactly, its the linear bearings that can be at different locations and force therefore the board to different positions, those are the ones that i am interested in!

[Oscar Benjamin]

On 11 June 2016 at 17:52, <janosc...@gmail.com> wrote:
>
> Yes, exactly, its the linear bearings that can be at different locations and
> force therefore the board to different positions, those are the ones that i
> am interested in!

Rather than thinking about x, y and theta think about the three pin
positions. Give them position vectors r1, r2 and r3. Each is
constrained by a linear bearing and so e.g. r1 = a1 + t1*b1 where a1
and b1 are known vectors and t1 is the unknown line parameter. We have
then three unknown scalars t1, t2, and t3. Although we don't know r1,
r2, or r3 we do know their pairwise distances d12, d13 and d23. This
gives three equations e.g. |r1-r2|**2 = d12**2. Substitute for r1, r2
and r3 into those and we get 3 bivariate quadratic equations for t1,
t2, and t3. I would expect that sympy can solve that quickly for
numeric coefficients (if you can give numeric values for a1, b1, d12,
etc.).

In your example it's clear from symmetry that if phi is a solution
then so is phi+pi so you should expect to get multiple solutions here.
This method will also give upside down solutions that you may need to
prune. Presumably also you'll need to check the solutions for
consistency with the limits on the line bearings.

--
Oscar

[janosc...@gmail.com]

I had the same idea earlier, but i dropped it because my intuition was, that three quadratic equations are worse than three linear and one quadratic equation :-)

Since you brought this approach up again, i tried it now, but sympy does not seem to find a solution.

You can check out my code here: http://pastebin.com/famnqkLC

You wrote you expected sympy to find a solution for numeric coefficients, but i need a symbolic solution because i want to proceed further (by differentiating with respect so some of the parameters for optimization), and don't want a sympy.solve step in each optimization step.

Any idea why my original approach "explodes" in regards of the resulting expressions?

[Oscar Benjamin]

On 14 June 2016 at 20:30, <janosc...@gmail.com> wrote:
> I had the same idea earlier, but i dropped it because my intuition was, that
> three quadratic equations are worse than three linear and one quadratic
> equation :-)
>
> Since you brought this approach up again, i tried it now, but sympy does not
> seem to find a solution.
>
> You can check out my code here: http://pastebin.com/famnqkLC
>
> You wrote you expected sympy to find a solution for numeric coefficients,
> but i need a symbolic solution because i want to proceed further (by
> differentiating with respect so some of the parameters for optimization),
> and don't want a sympy.solve step in each optimization step.

You should perhaps explain more clearly what you do want then. You
want <X> in terms of <Y> but what are X and Y?

> Any idea why my original approach "explodes" in regards of the resulting
> expressions?

Given your three linear equations you can solve them to get x, y, and
cosphi all in terms of sinphi. Then you have that cosphi is:

(-(Cux*Cvy - Cuy*Cvx)*(Awx*Cvz*sinphi - Awy*Cuz*sinphi + Cwz) +
(Cux*Cvz - Cuz*Cvx)*(Avx*Cvy*sinphi - Avy*Cuy*sinphi + Cwy) - (Cuy*Cvz
- Cuz*Cvy)*(Aux*Cvx*sinphi - Auy*Cux*sinphi + Cwx))/(Cux*Cvy*(Awx*Cuz
+ Awy*Cvz) - Cux*Cvz*(Avx*Cuy + Avy*Cvy) - Cuy*Cvx*(Awx*Cuz + Awy*Cvz)
+ Cuy*Cvz*(Aux*Cux + Auy*Cvx) + Cuz*Cvx*(Avx*Cuy + Avy*Cvy) -
Cuz*Cvy*(Aux*Cux + Auy*Cvx))

Which you can substitute for cosphi in cosphi**2+sinphi**2-1 to get a
quadratic in only sinphi.

The explosion comes from expanding the cosphi expression (and then
squaring!). You have 12 different symbols (excluding sinphi) and all
the cross-multiplications gives a combinatoric explosion of terms that
aren't easy to factor.

I guess what you want to do is rearrange the cosphi expression into
the form a*sinphi + b but without expanding a and b.

--
Oscar

[janosc...@gmail.com]

Sorry if i am being imprecise, but actually i cant describe it better than in my opening post. I want X, Y and Phi in terms of everything else.

Maybe i was wrong in assuming that i could simply pass sympy the system of equations and it would deliver the solutions in a form that i could then use for further processing (auto diff), but the intermediate expression you give for cosphi indicates that there is a much more compact representation of the solution in matrix notations. I understand that sympy only got some (non-matrix) equations, and therefore couldn't introduce matrix operations by itself.

But maybe there is some approach i am missing, which would help sympy to use matrix operations in the solution?