How do I use MATLAB.jl and vpa?

Specific Domains Numerics
1

Hello everyone,

I am want to do something conceptually simple using MATLAB.jl package. I have a BigFloat matrix in Julia, and I want to send to Matlab to compute its inverse/determinant/eigenvalues with arbitrary-precision arithmetic (vpa) and return back to Julia.
Assuming that Matlab is already in local PATH, the following example should work:

using MATLAB
using LinearAlgebra
using GenericLinearAlgebra

m = n = 2
a = rand(m, n)
x = mxarray(a)

# works OK
inv(a)  
mat"inv($x)"

det(a)
mat"det($x)"

eigen(a)
mat"eig($x)"


# NOT OK with big numbers
b = big.(a)
y = mxarray(b)
# Maybe a better way to write it could be
# mat"y = vpa($b,60)" 

det(b) 
mat"det($y)" # returns: Undefined function 'det' for input arguments of type 'cell'.

eigen(b) # ? how to make it work ?
mat"eig($y)" # returns a bigger message complaining about types

By the error messages I guess that MATLAB.jl convert Julia type in some data structure, and send to Matlab to be executed. Any ideas of how to handle my example with BigFloat ?

Thank you

2

Matlab’s variable-precision arithmetic library likely uses a different representation for its floats than GMP, and they probably can’t be translated without some amount of effort. Is there some reason you don’t want to do the calculations in Julia? These all work fine:

using LinearAlgebra, GenericLinearAlgebra

A = big.(rand(3,3) + I)
inv(A)
det(A)
eigvals(A)
1 Like
3

I am doing a mixture codes, some are already written in Matlab and uses vpa. I am just afraid to use BigFloat and don’t have the level of precision that I have with vpa.

By the way, how can I make eigvecs() to work with BigFloats ?

4

BigFloat has arbitrary precision. It looks like ArbitraryLinearAlgebra is limited to Hermitian matrices for eigenvectors:

using LinearAlgebra, GenericLinearAlgebra
setprecision(80)

A = rand(BigFloat,3,3) + I
hA = Hermitian(A*A')
eigvecs(hA)
5

I need to work with non-hermitian matrices and I will use matlab meanwhile to do the job - I hope that someday this operation can be done in Julia, it will help my a lot.

Thank you by the reply

6

Another option is to code it up yourself, using a standard algorithm from a book like Golub and Van Loan: that’s more or less what I did with GenericSVD.jl (you may even be able to reuse some of that code).

7

Using inverse iteration:

using LinearAlgebra, GenericLinearAlgebra
import LinearAlgebra.eigen

function eigen(A::Array{T,2}) where T <: Union{BigFloat,Complex{BigFloat}}
    ni, nj = size(A)
    Λ = eigvals(A)
    if norm(imag(Λ)) < eps(BigFloat)
        Λ = real(Λ)
    end
    
    V = rand(eltype(Λ), nj, nj)
    bkp = similar(V, nj)
    for j = 1:nj
        Δ = BigFloat(Inf)
        bk = view(V,:,j)
        f = qr(A - Λ[j]*I)
        i = 0
        while Δ > eps(BigFloat) && i < 10
            bkp .= f\bk
            ck = norm(bkp)
            bkp ./= ck
            Δ = norm(bkp - bk)
            bk .= bkp
            i += 1
        end
    end
    return Λ, V
end