As you can see in the following example, the Matrix.eigenvects method returns wrong results when the entries in the matrix are sympy.Float instances:
>>> A = Matrix([[Float('0.1'), Float('0.2')], [Float('0.8'), Float('0.5')]])
>>> A
Matrix([
[0.1, 0.2],
[0.8, 0.5]])
>>> A.eigenvects()
[(-0.147213595499958, 1, [Matrix([
[-0.628960169645094],
[ 0.777437524821136]])]), (0.747213595499958, 1, [Matrix([
[-0.355518216481267],
[ -1.15048111577287]])])]
>>>
Here are the expected results: https://www.wolframalpha.com/input/?i=eigenvectors+%7B%7B0.1%2C+0.2%7D%2C+%7B0.8%2C+0.5%7D%7D.
The output does look pretty messy but it seems to print out the following form:
Note that eigenvectors are not unique.
The following can be used to check that the decomposition is correct:
from sympy import *
A = Matrix([[Float('0.1'), Float('0.2')],
[Float('0.8'), Float('0.5')]])
decomposition = A.eigenvects()
d1, d2 = decomposition
val1, _, vec1 = d1
val2, _, vec2 = d2
vec1, vec2 = vec1[0], vec2[0]
D = diag(val1, val2)
Q = vec1.col_insert(1, vec2)
print(D) # same as Wolfram
print(Q) # you can multiply each col by a constant to get Wolfram's result
print(simplify(Q*D*Q**-1)) # A
Maybe A.diagonalize() might be simpler to work with? I think it provides normalized eigenvectors too.
Edit: I am not sure why Floats makes it a purely numeric computation while other classes such as Rational do not.
You can verify that these are eigenvectors with:
In [17]: [(l1, m1, (v1,)), (l2, m2, (v2,))] = A.eigenvects()
In [18]: A*v1 - l1*v1
Out[18]:
⎡ 0 ⎤
⎢ ⎥
⎣5.55111512312578e-17⎦
In [19]: A*v2 - l2*v2
Out[19]:
⎡ 0 ⎤
⎢ ⎥
⎣-1.11022302462516e-16⎦
The tiny floats are just residual rounding error.
I am not sure why
Floats makes it a purely numeric computation while other classes such asRationaldo not.
This is always the case in sympy. As soon as a float appears the whole expression is considered to be tainted as imprecise.
You can use nsimplify to approximate the Floats as Rational:
In [21]: A.applyfunc(nsimplify)
Out[21]:
⎡1/10 1/5⎤
⎢ ⎥
⎣4/5 1/2⎦
In [22]: A.applyfunc(nsimplify).eigenvals()
Out[22]:
⎧3 √5 3 √5 ⎫
⎨── - ──: 1, ── + ──: 1⎬
⎩10 5 10 5 ⎭
I don't think that the result is wrong if you have given floating point.
But it's wolfram alpha that automatically converts decimals to rationals.
I'd also note that Mathematica gives floating point results even if wolfram alpha may use it as its backend
In[1]:= Eigensystem[{{0.1, 0.2}, {0.8, 0.5}}]
Out[1]= {{0.747214, -0.147214}, {{-0.295242, -0.955423}, {-0.62896, 0.777438}}}
And I'd keep behavior in this way because I have seen many problems with eigenvalues not working with floats, and it's an impractical way to compute eigenvalues for floating points by rationalizing
Note that eigenvectors are not unique.
My bad, forgot about that while trying to fix a failing test case in https://github.com/mathics/Mathics/pull/821 😅️