Solving for a matrix?

[Matthias Geier]

Dear list.

I have this equation:

a = M * b,

where a and b are column vectors and M is a 4x4 matrix.
a and b consist of quite simple expressions, M is unknown:

>>> import sympy as sp
>>> a0, a1, a2, a3 = sp.symbols('a:4')
>>> a = sp.Matrix([a3, a2, a1, a0])
>>> b = sp.Matrix([a0, a3 + a2 + a1 + a0, a1, 3 * a3 + 2 * a2 + a1])
>>> M = sp.MatrixSymbol('M', 4, 4)
>>> sp.Eq(a, M * b)
Eq(Matrix([[a3], [a2], [a1], [a0]]), M*Matrix([
[ a0],
[a0 + a1 + a2 + a3],
[ a1],
[ a1 + 2*a2 + 3*a3]]))

What's the most straightforward way to solve this for M?

I have actually found a way, but it seems a bit non-obvious to me:

>>> def get_value(i, j):
... return b[i].as_coefficients_dict()[a[j]]
...
>>> M_inv = sp.Matrix(4, 4, get_value)
>>> M_inv.inv()
Matrix([
[ 2, -2, 1, 1],
[-3, 3, -2, -1],
[ 0, 0, 1, 0],
[ 1, 0, 0, 0]])

This is the solution I'm looking for, but I was hoping there is a
better way to get it.
Also, the way I did it I had to know that I could look for the inverse
and then undo the inverse afterwards, but I would prefer if there is a
solution without this manual step.

Thanks in advance!

cheers,
Matthias

[Leonid Kovalev]

You have 4 equations with 16 unknowns, so there are going to be a lot of solutions. 

Many SymPy functions struggle with M[i, j] construction which creates MatrixElement instead of an ordinary Symbol. I would use a Matrix filled with Symbols: 

M = sp.Matrix(sp.symarray('M', (4, 4)))
sp.solve(a - M*b, list(M))

This returns

{M_3_0: (-M_3_1*(a0 + a1 + a2 + a3) - M_3_2*a1 - M_3_3*(a1 + 2*a2 + 3*a3) + a0)/a0, M_2_0: (-M_2_1*(a0 + a1 + a2 + a3) - M_2_2*a1 - M_2_3*(a1 + 2*a2 + 3*a3) + a1)/a0, M_1_0: (-M_1_1*(a0 + a1 + a2 + a3) - M_1_2*a1 - M_1_3*(a1 + 2*a2 + 3*a3) + a2)/a0, M_0_0: (-M_0_1*(a0 + a1 + a2 + a3) - M_0_2*a1 - M_0_3*(a1 + 2*a2 + 3*a3) + a3)/a0}

which, as previously mentioned, is a lot of solutions. You can plug in some arbitrary numbers for the free variables here. 

[Matthias Geier]

Thanks a lot for your answer, Leonid!

I think I'm starting to see the problem ...
When I'm solving the system of equations, I get many solutions, most
of them containing a0, a1, a2 and a3.

However, I'm just interested in exactly one of those solutions where
all occurrences of a0, a1, a2 and a3 cancel out and the resulting
matrix only consists of concrete numbers.

Is there a (simple) way to tell this to SymPy?

In the meantime, I've found a different way to get the desired result.

Let me first repeat the original setup:

>>> import sympy as sp
>>> a0, a1, a2, a3 = sp.symbols('a:4')
>>> a = sp.Matrix([a3, a2, a1, a0])
>>> b = sp.Matrix([a0, a3 + a2 + a1 + a0, a1, 3 * a3 + 2 * a2 + a1])
>>> M = sp.MatrixSymbol('M', 4, 4)
>>> sp.Eq(a, M * b)
Eq(Matrix([[a3], [a2], [a1], [a0]]), M*Matrix([
[ a0],
[a0 + a1 + a2 + a3],
[ a1],
[ a1 + 2*a2 + 3*a3]]))

I stumbled upon
https://groups.google.com/forum/#!searchin/sympy/linear_eq_to_matrix/sympy/Wqs1OhTBexg/eGoKjrWHBwAJ,
which mentions sympy.solvers.solveset.linear_eq_to_matrix().
http://docs.sympy.org/latest/modules/solvers/solveset.html#linear-eq-to-matrix

With this I can get my expected solution:

>>> from sympy.solvers.solveset import linear_eq_to_matrix
>>> linear_eq_to_matrix(list(b), list(a))
(Matrix([
[0, 0, 0, 1],
[1, 1, 1, 1],

[0, 0, 1, 0],

[3, 2, 1, 0]]), Matrix([
[0],
[0],
[0],
[0]]))
>>> _[0].inv()

Matrix([
[ 2, -2, 1, 1],
[-3, 3, -2, -1],
[ 0, 0, 1, 0],
[ 1, 0, 0, 0]])

This still involves the extra step of inverting the result.

Is this already the best way to get this result with SymPy?

cheers,
Matthias

[Jason Moore]

Matrices have a number of solvers, check out: http://docs.sympy.org/latest/modules/matrices/matrices.html#linear-algebra

If you have a simple linear system, I recommend LUsolve().

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+unsubscribe@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/CAFesC-fH%3D8BwFi9ZVLL_4GBzUv3At363X9jhV6mtLOS62n4QMQ%40mail.gmail.com.

For more options, visit https://groups.google.com/d/optout.

[Leonid Kovalev]

I understood what's going on. The variables a0, a1, a2, a3 are a distraction since they are not to appear in the solution. What matters is their coefficients in b, which naturally form a matrix.

>>> sp.Matrix([[eq.coeff(v) for v in a] for eq in b])
Matrix([

[0, 0, 0, 1],
[1, 1, 1, 1],
[0, 0, 1, 0],
[3, 2, 1, 0]])

So that's what is really given. And the desired result is the inverse of that matrix. 

>>> sp.Matrix([[eq.coeff(v) for v in a] for eq in b]).inv()

Matrix([
[ 2, -2,  1,  1],
[-3,  3, -2, -1],
[ 0,  0,  1,  0],
[ 1,  0,  0,  0]])

This is the best way to get the result with SymPy, because the best way to find the inverse of a matrix is to invert that matrix. It's often said (correctly) that solving a linear system is more efficient than inverting its matrix. But here, the task was really to find the inverse of a 4 by 4 matrix.   

[Matthias Geier]

Thanks Leonid, that clears everything up for me!

cheers,
Matthias

> --
> 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/1118d33d-2448-41d5-ab6c-653009b0893a%40googlegroups.com.