Search code examples
juliaeigenvector

Eigenvector in Julia -- result less precise than Mathematica

I'm having a great time learning Julia.

With Wolfram Mathematica, I can get exact Eigenvectors for simple matrices. E.g.: enter image description here

The same example in Julia returns the following:

enter image description here

Is there a way to get Julia's output to be simpler/more exact/more like Mathematica's?

Thank you!

Solution

  • Mathematica tends to use exact, symbolic algorithms in situations like these, even though they don't generalize or scale to larger inputs. Julia uses the same algorithms on both large and small inputs, which tend to be numerical and not exact for integer inputs. The eigenvectors produced by eigvecs are always normalized so that they have unit norm:

    julia> B = eigvecs(A)
3×3 Matrix{Float64}:
 1.0  0.707107  0.557086
 0.0  0.707107  0.742781
 0.0  0.0       0.371391

julia> mapslices(norm, B, dims=1)
1×3 Matrix{Float64}:
 1.0  1.0  1.0

    Mathematica apparently doesn't do that, which lets it give eigenvectors with integral entries. In general, integral eigenvectors don't exist, which raises the question: how does Mathematica decide to use an exact algorithm here? Who knows. What algorithm does it use? Also, who knows. Mathematica is an amazing piece of software, but it's closed source so it's a black box. It's possible they document this, and someone can go look at the docs if they want and post another answer to this question, but unlike open source software, there's no way to actually look at what it's doing.

    If you suspect a matrix has eigenvectors that have integral entries, then you can use the following hack to try to recover them: for each column, append zero to the column and then find smallest positive difference between any two entries in the padded column; then divide the original column by that difference. Here's some code that implements this for your example:

    julia> B ./ mapslices(v -> minimum(filter(>(0), [v;0] .- [v;0]')), B, dims=1)
3×3 Matrix{Float64}:
 1.0  1.0  3.0
 0.0  1.0  4.0
 0.0  0.0  2.0

    As you can see, it recovers the integer-valued eigenvectors that Mathematica reports. I'm a little surprised that Mathematica returns the eigenvectors as rows rather than columns, but that seems to be what it's doing.