Julia is slower than matlab when it comes to matrix diagonalization?

General Usage
, ,
21

The code is below, reproduced from this post:

using ArnoldiMethod, LinearAlgebra, LinearMaps

PQ = P*Q

# Factorizes A and builds a linear map that applies inv(A) to a vector.
function construct_linear_map(A)
    F = factorize(A - ev*I)
    LinearMap{eltype(A)}((y, x) -> ldiv!(y, F, x), size(A,1), ismutating=true)
end

function rerun_ArnoldiMethod(n=1, PQ=PQ)
    Threads.nthreads() == n || error("Restart Julia with -t$n")
    BLAS.set_num_threads(n)
    println("ArnoldiMethod $n Threads:")
    @time begin
        decomp, history = partialschur(construct_linear_map(PQ), nev=NMODES, tol=1e-14, restarts=300, which=LM())
        D2, Exy = partialeigen(decomp)
    end
    println(history)
    D2 = inv.(D2) .+ ev
    @show D2
end

I’m sorry but I don’t have the time to dig up and test the code to recreate the matrix right now. As stated in my original post in that thread, the matrix PQ is a sparse, nonsymmetric, Float64 matrix of order 65274 with 1103319 stored entries.
ev is a scalar with value -12.25.

1 Like
22

There’s this question in the Matlab Answers site which shows how to get output from eigs showing what it is doing. It does look like the default for nonsymmetric problems is probably Krylov-Schur, but I don’t know if there are cases where it might switch over to use ARPACK.

2 Likes
23

Thank you!

My matrix pencil is different (A symmetric, semi-pd, B symmetric, pd). I expected the eigenvectors to be B-orthogonal. They are not. The residual ||Av - Bv*Lam|| is large. So, no, ArnoldiMethod did not succeed.

24

You are right, Arpack is no longer called in Matlab eigs.

2 Likes
26

I am now posting the original matlab code. I am really disappointed if julia is significantly slower than matlab when it comes large sparse matrix diagonalization, which is common problem in computational physics.

clear all;close all; clc;tic
global  dim
global  table

N=12;            % number of atoms
M=N;              % number of sites
% number of basis vectors calculated analytically
dim=prod(1:max(1,N+M-1))/prod(1:max(1,N))/prod(1:max(1,M-1))
% basis is the table to store the information of the basis vectors
basis=zeros(M,dim);
basis(1,1)=N;
s=1;
table=zeros(1,dim);
weight=sqrt(100*(1:M)+3);
table(1)=weight*basis(:,1);
while (basis(M,s)<N)
    s1=M-1;
    while (basis(s1,s)==0)
        s1=s1-1;
    end
    
     for s2=1:s1-1
            basis(s2,s+1)=basis(s2,s);
     end
     basis(s1,s+1)=basis(s1,s)-1;
     basis(s1+1,s+1)=N-sum(basis(1:s1,s+1));
     table(s+1)=weight*basis(:,s+1);
     
     s=s+1;
end
[table,ind]=sort(table);

% to generate the table storing information of the hamiltonian
COO1=zeros(3,dim);
COO2=zeros(3,dim*M*2);
for s1=1:dim
    COO1(1,s1)=s1;
    COO1(2,s1)=s1;
    COO1(3,s1)=0.5*sum(basis(:,s1).^2-basis(:,s1));
end
Int=sparse(COO1(1,:),COO1(2,:),COO1(3,:),dim,dim);
s=0;
target=[(2:M),1];
for s1=1:dim
    if (mod(s1,10000)==0)
        s1
    end
% s1
    for s2=1:M
         if (basis(s2,s1)>0)
             s=s+1;
             final=basis(:,s1);
             final(s2)=final(s2)-1;
             s3=target(s2);
             final(s3)=final(s3)+1;
             value=weight*final;
             index=search(value);
             
             COO2(1,s)=s1;
             COO2(2,s)=ind(index);
             COO2(3,s)=-sqrt(final(s3)*(final(s2)+1));
         end
    end
end
COO2=COO2(:,1:s);
Kin=sparse(COO2(1,:),COO2(2,:),COO2(3,:),dim,dim);
Kin=Kin+Kin';
% save Kin_Int_N10.mat    Kin  Int

toc

Vlist=0:0.5:20;
Cmax=zeros(1,length(Vlist));
tic
for s1=1:length(Vlist)
    s1
    H=Kin+Vlist(s1)*Int;
    tic
    [V,D]=eigs(H,1,'sa');    % to solve the smallest algebraic eigenvalues
    toc
    d=diag(D);
    [d,ind]=sort(d);
    gstate=V(:,ind(1));
    gstate2=abs(gstate').^2;
    Cmax(s1)=max(gstate2);
end

plot(Vlist,Cmax)
ylim([0 1])
xlim([0 max(Vlist)])
xlabel('V/J','fontsize',14)
ylabel('C^2_{max}','fontsize',14)
str=strcat('N=',num2str(N));
title(str,'fontsize',14)
% legend(str,'fontsize',14)
toc

this script calls a function which is put in another file

function   index=search(value)

global  dim
global  table

L=1;
R=dim;
M=floor((L+R)/2);

while (L<=R)
    if (value>table(M))
        L=M+1;
        M=floor((L+R)/2);
    elseif (value<table(M))
            R=M-1;
            M=floor((L+R)/2);
    elseif (value==table(M))
               index=M;
               return
    end
end
27

I find that if I replace this function with

using KrylovKit

function gstate(Kin, Interaction, Vlist)
    Cmax = zeros(length(Vlist))
    ge = zeros(length(Vlist))
    gstate = randn(size(Kin, 1))
    H  = spzeros(size(Kin)...)
    for s1 = 1:length(Vlist)
        # print(s1, '\n') # commented out because that's obnoxious
        H .= Kin .+ Vlist[s1] .* Interaction
        d, _gstate = eigsolve(H, gstate, 1, :SR)
        gstate  .= _gstate[1]
        Cmax[s1] = maximum(abs2, gstate)
        ge[s1] = d[1]
    end
    return ge, Cmax
end

then the new version of gstate runs about 4 times faster on my machine.

5 Likes
28

I want to give KrylovKit.jl another try, This time not using geneigsolve (that doesn’t work). However, I cannot get eigsolve working either:

convinfo = ConvergenceInfo: 18 converged values after 5000 iterations and 58012 applications of the linear map;                                                          
norms of residuals are given by (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0014937282975521534, 0.026181378600122828).

I am calling it like this:

    Mfactor = cholesky(Symmetric(M))
    z = complex(zero(eltype(M)))
    d, vv, convinfo = eigsolve(
        x -> Mfactor \ (K * x),
        rand(typeof(z), size(K, 1)),
        nev,
        :SR;
        maxiter = maxiter
    )

The problem is (K - lam * M) * v = 0. I convert to (M^-1 * K - lam * 1) * v = 0. Matrix size (624, 624) .

The first six eigenvalues are zero, then come 2, 3, 3, 1, 2, 3 repeated eigenvalues.

I must be missing something. Is there an improvement that someone could suggest?

29

If the problem is definite, maybe you should be using shift-invert? (The Julia libraries won’t switch to shift-invert automatically.)

It’s hard to give more concrete advice without knowing your matrices.

It’s also not clear what is “not working” — the residual norms are small, so what is incorrect about the answers?

2 Likes
30

Only 18 eigenvalues converged out of twenty, using a ridiculous number of iterations. And even with 15000 iterations, I cannot get all of those eigenvalues to converge.

31

Edit: PLEASE DISREGARD.
Perhaps it is the Krylov-Schur method. I cannot make Matlab eigs to converge either. Matlab says it did, but it did not:

=== Simple eigenvalue problem A*x = lambda*x ===

The eigenvalue problem is real non-symmetric.

Computing 20 eigenvalues of type 'smallestabs'.


Parameters passed to Krylov-Schur method:
  Maximum number of iterations: 500
  Tolerance: 2.2204e-16
  Subspace Dimension: 40

Find eigenvalues of A\x = mu*x.
Compute decomposition of A...

--- Start of Krylov-Schur method ---
Iteration   1:  8 of 20 eigenvalues converged. Smallest non-converged residual 5.7e-16 (tolerance 2.2e-16).
Iteration   2:  8 of 20 eigenvalues converged. Smallest non-converged residual 6.5e-16 (tolerance 2.2e-16).
Iteration   3: 18 of 20 eigenvalues converged. Smallest non-converged residual 3.2e-14 (tolerance 2.2e-16).
Iteration   4: 20 of 20 eigenvalues converged.
Elapsed time is 0.033472 seconds.

resid =

  145.8667

The residual resid = norm(M \ (K * v) - v * d) is definitely not zero.

1 Like
32

it is what I usually observe… you need to tune the tolerance, the krylov dim, etc for this to work.

You can also pre-factor Mfactor (lu?)

We cannot know that the spectrum of Mfactor^(-1)K is “compact”. Dont you have additional mathematical properties?

1 Like
33

Yes, M pos def, K pos semi def, both symm. I know accurately all 20 eigenvalues.

https://github.com/PetrKryslUCSD/GEPHelpers.jl:

Arpack.jl, SubSIt.jl converge quickly and accurately.

ArnoldiMethod.jl converges on the eigenvalues, not the eigenvectors.

KrylovKit.jl does not converge.

Edit: ArnoldiMethod.jl and KrylovKit.jl are in fact giving correct solutions. I just need to convert the standard eigenvectors to mass-orthogonal. Next, I will investigate the relative speeds of the above packages. I will post in a separate thread.

34
35

I don’t have a huge amount of practical experience using Krylov methods and this is a shot in the dark that is not based on very much, but seeing an absolute residual like that makes me immediately think of the possibility of not taking into account \|M^{-1} K\|. I would be curious about the relative residual:

resid / ((opnorm(M \ K) + abs(d)) * norm(v))

It seems that you are testing with small enough matrices that you could look at that easily. A small relative residual implies that you have an eigenvalue and eigenvector for a matrix close to M^{-1}K in a norm-wise relative sense and you can’t really expect to do better than that with this approach. Along those lines, poking into the KrylovKit.jl code, it looks like it might default to determining convergence by applying a fixed absolute tolerance of 10^{-12} to a Hessenberg subdiagonal. I could have missed something in the very short time I spent looking at it, but I didn’t notice anywhere in the code that scales that tolerance to take into account the possibility that the Hessenberg matrix in the Krylov factorization might have large norm. If that’s the case and \|M^{-1}K\| is large, then the code might be defaulting to an unrealistic convergence tolerance for your problem.

36

The norm of M\K is ~4e4.

37

My bad. I forgot to shift to account for the singular K.

Matlab ([v, d] = eigs(M \ (K + 0.1 * M), 20, 'SM', tol = eps, maxiter = 500, Display = true);):

=== Simple eigenvalue problem A*x = lambda*x ===

The eigenvalue problem is real non-symmetric.

Computing 20 eigenvalues of type 'smallestabs'.


Parameters passed to Krylov-Schur method:
  Maximum number of iterations: 500
  Tolerance: 2.2204e-16
  Subspace Dimension: 40

Find eigenvalues of A\x = mu*x.
Compute decomposition of A...

--- Start of Krylov-Schur method ---
Iteration   1:  8 of 20 eigenvalues converged. Smallest non-converged residual 1.7e-11 (tolerance 2.2e-16).
Iteration   2:  9 of 20 eigenvalues converged. Smallest non-converged residual 3.3e-12 (tolerance 2.2e-16).
Iteration   3:  9 of 20 eigenvalues converged. Smallest non-converged residual 8.9e-15 (tolerance 2.2e-16).
Iteration   4: 14 of 20 eigenvalues converged. Smallest non-converged residual 3.2e-16 (tolerance 2.2e-16).
Iteration   5: 17 of 20 eigenvalues converged. Smallest non-converged residual 1.7e-11 (tolerance 2.2e-16).
Iteration   6: 17 of 20 eigenvalues converged. Smallest non-converged residual 4.4e-14 (tolerance 2.2e-16).
Iteration   7: 20 of 20 eigenvalues converged.
Elapsed time is 0.147191 seconds.

resid =

   4.1613e-08

KrylovKit.jl:

convinfo = ConvergenceInfo: 23 converged values after 1634 iterations and 26058 applications of the linear map;

Eigenvalues converged, residual a tad bigger than for Matlab.

So now to reduce that ridiculous amount of iterations.

38

From the Matlab eigs documentation:

using 'smallestabs', which employs a Krylov method using the inverse of A.

This is what I suggested above: the Julia libraries won’t do shift/invert automatically (whereas Matlab will). This can drastically change the number of iterations.

So to get a comparable number of iterations to Matlab with M \ (K + 0.1 * M) (which is internally inverted) you’ll want to do something like:

F = cholesky(Symmetric(K + 0.1 * M))
eigsolve(x -> F \ (M * x), ..., :SM)

(:SM is the equivalent of smallestabs in Matlab, not :SR, though :SR should work as well since all the eigenvalues are real and positive). Then invert and shift the resulting eigenvalues to get the eigenvalues of the original problem.

Correction: Once you invert, you should be using :LM to look for largest-magnitude eigenvalues (rather than :SM).

1 Like
39

I think you meant:

    Kfactor = cholesky(Symmetric(K))
    di, vv, convinfo = eigsolve(
        x -> Kfactor \ (M * x),
        rand(typeof(z), size(K, 1)),
        nev,
        :LR;
        maxiter = maxiter,
        krylovdim = 2 * nev
    )
    d = 1 ./ di

Yep. That dramatically reduced the number of iterations.

40

That’s not equivalent to your Matlab code, which added a shift to K (presumably to make it positive definite).

41

Yes, this is inside a function that takes K already shifted. I was just pointing to :LR… :slight_smile:

1 Like