CODEWITHSUNDEEP
pythonnumpyscipysympy
Nikita
Nikita

Reputation: 61

How to solve equation with two matrices in python?

I need to find matrix A that A*B*A = B*A*B, and A^2 =1 if B =(1,0,0)(0,-1,0)(0,0,1)

I tried that do it in sympy:sympy.solve(A*B*A - B*A*A)

import sympy as sp
a,b,c,d,e,f,g,h,k = sp.symbols('a b c d e f g h k')
B = sp.Matrix([[1,0,1],[0,-1,0],[0,0,1]])
A = sp.Matrix([[a,b,c],[d,e,f],[g,h,k]])
sp.solve([A*B*A - B*A*B,B**2-sp.eye(3)])

and also I tried numpy.linalg, but I don't have any results

Upvotes: 2

Views: 835

Answers (1)

user6655984
user6655984

Reputation:

First of all, linalg will not help, as this is not a linear problem: unknown A gets multiplied by itself. You want to solve a system of 18 quadratic equations with 9 unknowns. For a generic system we'd expect no solutions, but there is a lot of structure here.

In my version of SymPy (1.1.1) direct attempts to solve even one of the matrix equations A*B*A=B*A*B or A*A=I fail to finish in reasonable time. So let's follow the advice of saintsfan342000 and approach the problem numerically, as a minimization problem. This is how I did it:

import numpy as np
from scipy.optimize import minimize    
B = np.array([[1,0,0], [0,-1,0], [0,0,1]])
def func(A, B):
    A = A.reshape((3, 3))
    return np.linalg.norm(A.dot(B).dot(A)-B.dot(A).dot(B))**2 + np.linalg.norm(A.dot(A)-np.eye(3))**2
while True:
    guess = np.random.uniform(-2, 2, size=(9,))
    res = minimize(func, guess, args=(B,))
    if res.fun < 1e-15:
        A = res.x.reshape((3, 3))
        print(A)

The function to minimize is the sum of squares of Frobenius norms of A*B*A-B*A*B and A*A-I. I put minimization in a loop because there are some local minima where minimize will get stuck; so when the minimal value found is not sufficiently close to zero, I ignore the result and start over. After running for a while, the script will print a bunch of matrices like

[[ 0.70386835  0.86117949 -1.40305355]
 [ 0.17193376  0.49999999  0.81461157]
 [-0.25409118  0.73892171 -0.20386834]]

which all share two important features:

  • the central element A[1,1] is 1/2
  • the trace of the matrix (the sum of diagonal elements) is 1.

Let's use this information to help SymPy solve the system. I still don't want to throw both equations at it, so I try to get one at a time.

from sympy import *
var('a:h')                 #  a quick way to declare a bunch of symbols
B = Matrix([[1, 0, 0], [0, -1, 0], [0, 0, 1]])
A = Matrix([[a, b, c], [d, S(1)/2, f], [g, h, S(1)/2-a]])   # S(1)/2 is a way to get rational 1/2 instead of decimal 0.5
print(solve(A*B*A - B*A*B))
print(solve(A*A - eye(3)))

Now solve succeeds and prints the following:

[{b: h*(4*f*h - 3)/(2*g), d: -g/(2*h), a: 2*f*h - 1/2, c: -f*h*(4*f*h - 3)/g}]
[{b: h*(4*f*h - 3)/(2*g), d: -g/(2*h), a: 2*f*h - 1/2, c: -f*h*(4*f*h - 3)/g}]

Whoa! With the two constraints that we found numerically, both matrix equations are equivalent! I did not expect that. So we already have the solution: 

A = Matrix([[2*f*h - S(1)/2, h*(4*f*h - 3)/(2*g), -f*h*(4*f*h - 3)/g], [-g/(2*h), S(1)/2, f], [g, h, 1 - 2*f*h]])

for arbitrary f, g, h.

Note: A=B is a trivial solution, which was excluded above by the requirement of A[1,1]=1/2. I imagine this was the intent; it appears you were looking for faithful 3-dimensional representations of the symmetric group S3.

Upvotes: 3

Related Questions