CODEWITHSUNDEEP
matlabjulialinear-algebraeigenvectormatrix-factorization
Eric Johnson
Eric Johnson

Reputation: 714

Replicate MATLAB's `eig` in Julia

It is known that the Eigen decomposition of a matrix isn't unique.

Yet I wonder if one could replicate MATLAB's eig function in Julia.
Namely have exact output up to numeric errors assuming ascending order of the eigen values.

Assuming we set the options of MATLAB's eig to their default (balanceOption, balance, Using cholesky for symmetric matrices and qz for others as algorithms).

Since MATLAB is using Intel MKL so assuming also MKL.jl is available in Julia (does it have an effect on Julia for the eigen function?).

Upvotes: 1

Views: 372

Answers (1)

Bill
Bill

Reputation: 6086

Here is a recipe for deciding on a standardized eigenvector. There are arbitrary choices here to normalize the vector and to sort complex numbers as tuples of (real, complex) and to choose the signs of the eigenvector components so as to have the largest real value in the normalized vector (see comments):

using LinearAlgebra

""" compare two complex values as if tuples with (real portion, imaginary portion) """
complexlt(x, y) = (real(x), imag(x)) < (real(y), imag(y))

""" get a kind of canonical normalized form of an eigenvector of matrix A """
function sorted_normalized_eigen(A)
    # julia's eigen returns a struct, with eigenvector as values component of the struct
    m = eigen(A).values
    # sort the normalized form of the eigenvector from least to greatest
    normed1 = sort(normalize(m), lt = complexlt)
    # sort the additive inverse of normalized form of the eigenvector from least to greatest
    normed2 = sort(-1 .* normed1, lt = complexlt)
    # choose the normalized vector with the largest real component as our unique return value
    return complexlt(last(normed1), last(normed2)) ? normed2 : normed1
end

B =  [1 2 3; 3 1 2; 2 3 1]

@show sorted_normalized_eigen(B)

sorted_normalized_eigen(B) = ComplexF64[-0.2314550249431379 - 0.13363062095621225im, -0.2314550249431379 + 0.13363062095621225im, 0.9258200997725515 + 0.0im]

Upvotes: 1

Related Questions