Saleh Jafari
Saleh Jafari

Reputation: 1

Why does finding the eigenvalues of a 4*4 matrix by z3py take so much time and do not give any solutions?

I'm trying to calculate the eigenvalues of a 4*4 matrix called A in my code (I know that the eigenvalues are real values). All the elements of A are z3 expressions and need to be calculated from the previous constraints. The code below is the last part of a long code that tries to calculate matrix A, then its eigenvalues. The code is written as an entire but I've split it into two separate parts in order to debug it: part 1 in which the code tries to find the matrix A and part 2 which is eigenvalues' calculation. In part 1, the code works very fast and calculates A in less than a sec, but when I add part 2 to the code, it doesn't give me any solutions after.

I was wondering what could be the reason? Is it because of the order of the polynomial (which is 4) or what? I would appreciate it if anyone can help me find an alternative way to calculate the eigenvalues or give me some hints on how to rewrite the code so it can solve the problem.

(Note that A2 in the actusl code is a matrix with all of its elements as z3 expressions defined by previous constraints in the code. But, here I've defined the elements as real values just to make the code executable. In this way, the code gives a solution so fast but in the real situation it takes so long, like days. for example, one of the elements of A is almost like this:

0 +
 1*Vq0__1 +
 2 * -Vd0__1 +
 0 +
 ((5.5 * Iq0__1 - 0)/64/5) * 
 (0 +
  0 * (Vq0__1 - 0) +
  -521702838063439/62500000000000 * (-Vd0__1 - 0)) +
 ((.10 * Id0__1 - Etr_q0__1)/64/5) * 
 (0 +
  521702838063439/62500000000000 * (Vq0__1 - 0) +
  0.001 * (-Vd0__1 - 0)) +
 0 +
 0 + 0 +
 0 +
 ((100 * Iq0__1 - 0)/64/5) * 0 +
 ((20 * Id0__1 - Etr_q0__1)/64/5) * 0 +
 0 +
 -5/64

All the variables in this example are z3 variables.)

from z3 import *
import numpy as np


def sub(*arg):
    counter = 0
    for matrix in arg:
        if counter == 0: 
            counter += 1
            Sub = [] 
            for i in range(len(matrix)):
                Sub1 = []
                for j in range(len(matrix[0])):
                    Sub1 += [matrix[i][j]]
                Sub += [Sub1]
        else:
            row = len(matrix)
            colmn = len(matrix[0])
            for i in range(row):
                for j in range(colmn):
                    Sub[i][j] = Sub[i][j] - matrix[i][j]  
    return Sub

Landa = RealVector('Landa', 2) # Eigenvalues considered as real values
LandaI0 = np.diag(  [ Landa[0] for i in range(4)]  ).tolist()

ALandaz3 = RealVector('ALandaz3', 4 * 4 )

############# Building ( A - \lambda * I ) to find the eigenvalues ############

A2 = [[1,2,3,4],
      [5,6,7,8],
      [3,7,4,1],
      [4,9,7,1]]

s = Solver()

for i in range(4):
    for j in range(4):
        s.add( ALandaz3[ 4 * i + j ] == sub(A2, LandaI0)[i][j] )
ALanda = [[ALandaz3[0], ALandaz3[1], ALandaz3[2], ALandaz3[3] ],
          [ALandaz3[4], ALandaz3[5], ALandaz3[6], ALandaz3[7] ],
          [ALandaz3[8], ALandaz3[9], ALandaz3[10], ALandaz3[11]],
          [ALandaz3[12], ALandaz3[13], ALandaz3[14], ALandaz3[15] ]]
Determinant = (
 ALandaz3[0] * ALandaz3[5] * (ALandaz3[10] * ALandaz3[15] - ALandaz3[14] * ALandaz3[11]) -
 ALandaz3[1] * ALandaz3[4] * (ALandaz3[10] * ALandaz3[15] - ALandaz3[14] * ALandaz3[11]) +
 ALandaz3[2] * ALandaz3[4] * (ALandaz3[9]  * ALandaz3[15] - ALandaz3[13] * ALandaz3[11]) -
 ALandaz3[3] * ALandaz3[4] * (ALandaz3[9]  * ALandaz3[14] - ALandaz3[13] * ALandaz3[10]) )

tol = 0.001 

s.add( And( Determinant >= -tol, Determinant <= tol ) )   # giving some flexibility instead of equalling to zero

print(s.check())
print(s.model())

Upvotes: 0

Views: 147

Answers (1)

JohanC
JohanC

Reputation: 80459

Note that you seem to be using Z3 for a type of equations it absolutely isn't meant for. Z is a sat/smt solver. Such a solver works internally with a huge number of boolean equations. Integers and fractions can be converted to boolean expressions, but with general floats Z3 quickly reaches its limits. See here and here for a lot of typical examples, and note how floats are avoided.

Z3 can work in a limited way with floats, converting them to fractions, but doesn't work with approximations and accuracies as in needed in numerical algorithms. Therefore, the results are usually not what you are hoping for.

Finding eigenvalues is a typical numerical problem, where accuracy issues are very tricky. Python has libraries such as numpy and scipy to efficiently deal with those. See e.g. numpy.linalg.eig.

If, however your A2 matrix contains some symbolic expressions (and uses fractions instead of floats), sympy's matrix functions could be an interesting alternative.

Upvotes: 0

Related Questions