How can I give a function to eigs() in julia like matlab

General Usage
, , , ,
1

I have the following matlab code:

function [ DomRVec, DomEig, Error ] = MyEigs( Tup, Tdn, NumofValue, SpeString, Dbond, Dspin, InitVec )

VecDim = Dbond^2; 
Tup = reshape( permute( Tup, [ 1, 3, 2 ] ), [ Dbond * Dspin, Dbond ] );                      % A(lrm) => A(lm,r), placed at upper
Tdn = reshape( permute( Tdn, [ 3, 2, 1 ] ), [ Dspin * Dbond, Dbond ] );                      % B(lrm) => B(mr,l), placed at lower
ActionFun = @( x ) Action( Tup, Tdn, x, Dbond, Dspin );
[ DomRVec, DomEig ] = eigs( ActionFun, VecDim, NumofValue, SpeString, 'IsFunctionSymmetric', 0, 'Tolerance', 1E-14, 'MaxIterations', 300, 'StartVector', InitVec, ...
    'Display', 0 );

Error = Maxdiff( Action( Tup, Tdn, DomRVec, Dbond, Dspin ), DomRVec * DomEig );

end

function [ y ] = Action( Tup, Tdn, x, Dbond, Dspin )

d = size( x, 2 );
x = reshape( x, [ Dbond, Dbond * d ] );

if( d == 1 )
    y = reshape( Tup * x, [ Dbond, Dspin * Dbond ] );
    y = reshape( y * Tdn, [ Dbond^2, 1 ] );
elseif( d > 1) 
    y = reshape( Tup * x, [ Dbond, Dspin * Dbond, d ] );
    y = reshape( permute( y, [ 1, 3, 2 ] ), [ Dbond * d, Dspin * Dbond ] );
    y = reshape( y * Tdn, [ Dbond, d, Dbond ] );
    y = reshape( permute( y, [ 1, 3, 2 ] ), [ Dbond * Dbond, d ] );    
end

end

After referring to the documentation and source code of Arpack.jl, I wrote the following code. It runs without errors, but the results are significantly different from those obtained with MATLAB. Moreover, the results vary each time I run the code. I can’t figure out what’s causing this inconsistency. Could anyone point out where the mistake might be, or suggest alternative methods to achieve the same functionality?

using LinearAlgebra
using Arpack

struct MyMatrix{T,S} <: AbstractArray{T, 2}
    Tup::S
    Tdn::S
    Dbond::Int
    Dspin::Int
end

function MyMatrix(Tup, Tdn, Dbond, Dspin)
    Tnew = typeof(zero(eltype(Tup))/sqrt(one(eltype(Tup))))
    Tup_new = convert(AbstractMatrix{Tnew}, Tup)
    Tdn_new = convert(AbstractMatrix{Tnew}, Tdn)
    MyMatrix{Tnew,typeof(Tup_new)}(Tup_new, Tdn_new, Dbond, Dspin)
end

function LinearAlgebra.mul!(y::StridedVector{T}, A::MyMatrix{T}, x::StridedVector{T}) where T
    d = 1
    x = reshape(x, (A.Dbond, A.Dbond * d))
    if d == 1
        y = reshape(A.Tup * x, (A.Dbond, A.Dspin * A.Dbond))
        y = reshape(y * A.Tdn, (A.Dbond^2, 1))

    return y
end

Base.size(A::MyMatrix)  = (A.Dbond^2, A.Dbond^2)
LinearAlgebra.ishermitian(A::MyMatrix) = true

Tup, Tdn, Dbond, Dspin = reshape(1:18,3,3,2),reshape(2:19,3,3,2),3,2
Tup = reshape( permutedims( Tup, [ 1, 3, 2 ] ), Dbond * Dspin, Dbond )
Tdn = reshape( permutedims( Tdn, [ 3, 2, 1 ] ), Dspin * Dbond, Dbond )
A = MyMatrix(Tup, Tdn, Dbond, Dspin)
λ, ϕ = eigs(A; nev=1,which=:LR)
2

Did you try LinearMaps?

3

Note that nowadays there are pure-julia implementations of restarted Arnoldi, which may be more convenient than Arpack, such as KrylovKit.jl.

Note that if your linear operator is Hermitian, you should make sure the iterative solver knows this since better methods are available in this case. (In that case, you could also try LOBPCG in IterativeSolvers.jl.)

3 Likes
4

thank u! It’s solved by KrylovKit.jl, but there are still some minor issues, the keyword parameter howany cannot accurately control the number of output eigenvalues/eigenvectors. Take the following code as an example, the number of the output eigenvalues vary each time though I set howany==1:

function Action(Tup, Tdn, x, Dbond, Dspin)
    d = size(x)[2]
    x = reshape(x, Dbond, Dbond * d)
    if d == 1
        y = reshape(Tup * x, Dbond, Dspin * Dbond)
        y = reshape(y * Tdn, Dbond^2,1)
    return y
end
Tup, Tdn, Dbond, Dspin = rand(10,10,2),rand(10,10,2),10,2
Tup = reshape( permutedims( Tup, [ 1, 3, 2 ] ), Dbond * Dspin, Dbond )
Tdn = reshape( permutedims( Tdn, [ 3, 2, 1 ] ), Dspin * Dbond, Dbond );
x0=ones(100,1);

ActionFun = x -> Action(Tup, Tdn, x, Dbond, Dspin);
DomEig, DomRVec = eigsolve(ActionFun,x0,1,:LR);

I don’t know why this is happening :thinking:

5

from the docs (at least):

Compute at least howmany eigenvalues from the linear map encoded in the matrix A

Basically it will stop when at least 1 ev has been found

1 Like
6

Alright, it looks like I asked a stupid question :face_with_open_eyes_and_hand_over_mouth:, thank u.

7

No, but it’s solved by KrylovKit.jl, thanks for your suggestion