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

**URL:** <https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658>\
**Category:** General Usage\
**Tags:** question, eigs, eigenvalue-problem\
**Created:** [January 11, 2024, 3:42am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658 "2024-01-11T03:42:27Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![jiang\_ming\_zhang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jiang_ming_zhang/32/204063_2.png) [@jiang\_ming\_zhang](https://discourse.julialang.org/u/jiang_ming_zhang)\
**Post date:** [January 11, 2024, 3:42am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/1 "2024-01-11T03:42:27Z")

</div>

The following julia script is a literal translation of the original matlab code,

```julia
using LinearAlgebra, SparseArrays
using Arpack, Plots

function bhm_basis(N, M, dim, weight)
    basis = zeros(UInt8, M, dim)
    basis[1, 1] = N
    s = 1
    while basis[M, s] < N
        s1 = M - 1
        while basis[s1, s] == 0
            s1 = s1 - 1
        end

        basis[1:s1-1, s+1] = basis[1:s1-1, s]
        basis[s1, s+1] = basis[s1, s] - 1
        basis[s1+1, s+1] = N - sum(basis[1:s1, s+1])

        s = s + 1
    end
    table = (weight * basis)'
    ind = sortperm(table)
    table = table[ind]
    mini = minimum(table[2:end] - table[1:end-1]) / 4.0
    return basis, table, ind, mini
end

function hint(dim, basis, M)
    COO1 = zeros(3, dim)
    for s1 = 1:dim
        COO1[1, s1] = s1
        COO1[2, s1] = s1
        COO1[3, s1] = 0.5 * sum(basis[:, s1] .^ 2 - basis[:, s1])
    end
    return Interaction = sparse(COO1[1, :], COO1[2, :], COO1[3, :], dim, dim)
end

function hkin(dim, basis, M, table, ind, mini)
    COO2 = zeros(3, dim * M * 2)
    s = 0
    target = [Array(2:M); 1]
    for s1 = 1:dim
        if mod(s1, 10000) == 0
            println(s1)
        end
        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 = searchsortedfirst(table, value - mini)

                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)
    return Kin = Kin + Kin'
end

function gstate(Kin, Interaction, Vlist)
    Cmax = zeros(length(Vlist))
    ge = zeros(length(Vlist))
    for s1 = 1:length(Vlist)
          print(s1, '\n')
        H = Kin + Vlist[s1] * Interaction
        d, gstate = eigs(H, nev=1, which=:SR) # to solve the smallest algebraic eigenvalues
        gstate2 = abs2.(gstate)
        Cmax[s1] = maximum(gstate2)
        ge[s1] = d[1]
    end
    return ge, Cmax
end

@time begin
    N = 12
    M = 12
    dim = binomial(N + M - 1, M - 1)
    weight = sqrt.(100 * collect(1:M) .+ 3.0)'
    basis, table, ind, mini = bhm_basis(N, M, dim, weight)
    Interaction = hint(dim, basis, M)
    Kin = hkin(dim, basis, M, table, ind, mini)

    Vlist = 0:0.5:20
    ge, Cmax = gstate(Kin, Interaction, Vlist)
end

plot(Vlist, Cmax)
plot(Vlist, ge)

```

It is slower than matlab. Actually, it is much faster than matlab for all functions except for the function ‘gstate’, which calls the ‘eigs’ function from Arpack.

It seems that matlab is really good at matrix computation. I checked it. On my desktop, with matlab, each eigs call costs 9 seconds, while for julia, it is 14-18 seconds.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [January 11, 2024, 4:05am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/2 "2024-01-11T04:05:24Z")

</div>

Does `using MKL` at the top fix it?

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [January 11, 2024, 4:58am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/3 "2024-01-11T04:58:35Z")

</div>

Note that Arpack.jl is a wrapper around the FORTRAN Arpack code, so what you seem to observe is FORTRAN running slower than Matlab, which is a bit surprising. Deep inside, Matlab is probably calling the same library. Worth also trying out ArnoldiMethod.jl, which is a Julia translation of Arpack, and provides `eigs` as well.

Is it possible that the tolerance that Matlab uses to assess convergence is different, or that it runs fewer iterations?

---

<div class="post-metadata">

**Author:** ![mkitti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkitti/32/12459_2.png) [@mkitti](https://discourse.julialang.org/u/mkitti)\
**Post date:** [January 11, 2024, 5:05am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/4 "2024-01-11T05:05:32Z")

</div>

There are some opportunities to create view here instead of allocating new arrays.

> [@jiang\_ming\_zhang](#):
>
> ```julia
> basis[1:s1-1, s+1] = @view basis[1:s1-1, s]
> basis[s1, s+1] = basis[s1, s] - 1
> basis[s1+1, s+1] = N - sum(@view basis[1:s1, s+1])
> 
> ```

We can avoid the creation an intermediate array.

> [@jiang\_ming\_zhang](#):
>
> ```julia
> Cmax[s1] = maximum(abs2, gstate)
> 
> ```

---

<div class="post-metadata">

**Author:** ![jiang\_ming\_zhang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jiang_ming_zhang/32/204063_2.png) [@jiang\_ming\_zhang](https://discourse.julialang.org/u/jiang_ming_zhang)\
**Post date:** [January 11, 2024, 5:08am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/5 "2024-01-11T05:08:50Z")

</div>

actually, even without ‘view’, julia is much faster than matlab here. The problem is the eigs calls.

---

<div class="post-metadata">

**Author:** ![jiang\_ming\_zhang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jiang_ming_zhang/32/204063_2.png) [@jiang\_ming\_zhang](https://discourse.julialang.org/u/jiang_ming_zhang)\
**Post date:** [January 11, 2024, 5:14am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/6 "2024-01-11T05:14:45Z")

</div>

no.

no improvement with using MKL

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [January 11, 2024, 5:21am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/7 "2024-01-11T05:21:15Z")

</div>

Could you try `using MKLSparse` instead of MKL?

---

<div class="post-metadata">

**Author:** ![mkitti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkitti/32/12459_2.png) [@mkitti](https://discourse.julialang.org/u/mkitti)\
**Post date:** [January 11, 2024, 5:45am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/8 "2024-01-11T05:45:03Z")

</div>

> [@jiang\_ming\_zhang](#):
>
> `eigs(H`

What is the size of `H`? I’m not entirely sure that `eigs` is the best function to use.

---

<div class="post-metadata">

**Author:** ![jiang\_ming\_zhang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jiang_ming_zhang/32/204063_2.png) [@jiang\_ming\_zhang](https://discourse.julialang.org/u/jiang_ming_zhang)\
**Post date:** [January 11, 2024, 6:11am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/9 "2024-01-11T06:11:29Z")

</div>

the size of the matrix is about 1.35 million.

---

<div class="post-metadata">

**Author:** ![mkitti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkitti/32/12459_2.png) [@mkitti](https://discourse.julialang.org/u/mkitti)\
**Post date:** [January 11, 2024, 6:48am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/10 "2024-01-11T06:48:39Z")

</div>

1.35 million bytes? Elements? On one side?

---

<div class="post-metadata">

**Author:** ![jiang\_ming\_zhang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jiang_ming_zhang/32/204063_2.png) [@jiang\_ming\_zhang](https://discourse.julialang.org/u/jiang_ming_zhang)\
**Post date:** [January 11, 2024, 7:13am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/11 "2024-01-11T07:13:14Z")

</div>

1.35 million by 1.35 million

---

<div class="post-metadata">

**Author:** ![JM\_Beckers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jm_beckers/32/22482_2.png) [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)\
**Post date:** [January 11, 2024, 7:41am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/12 "2024-01-11T07:41:34Z")

</div>

> [@jishnub](#):
>
> tolerance that Matlab uses to assess convergence is different

This might be an explanation, in matlab the default tolerance seems to be 1E-14

> **[Subset of eigenvalues and eigenvectors - MATLAB eigs
- MathWorks Benelux](https://nl.mathworks.com/help/matlab/ref/eigs.html#bu2_q3e-2)**
>
> This MATLAB function returns a vector of the six largest magnitude eigenvalues of matrix A.

whereas in Julia it is set to 0.0

```julia
?eigs

# eigs(A; nev=6, ncv=max(20,2*nev+1), which=:LM, tol=0.0, maxiter=300, sigma=nothing, ritzvec=true, explicittransform=:auto, v0=zeros((0,)), check=0) -> (d,[v,],nconv,niter,nmult,resid)

```

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [January 11, 2024, 2:45pm UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/13 "2024-01-11T14:45:47Z")

</div>

> [@JM\_Beckers](#):
>
> [Subset of eigenvalues and eigenvectors - MATLAB eigs - MathWorks Benelux](https://nl.mathworks.com/help/matlab/ref/eigs.html#bu2_q3e-2)

I believe that means that in Arpack (eigs) the tolerance is reset to machine precision. Perhaps the OP could try setting tol to the same number?

---

<div class="post-metadata">

**Author:** ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)\
**Post date:** [January 11, 2024, 3:38pm UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/14 "2024-01-11T15:38:24Z")

</div>

Recent versions of Matlab don’t use Arpack. Perhaps the algorithms in KrylovKit.jl or ArnoldiMethod.jl would be more efficient for this problem (IIRC they are similar to the newer scheme in Matlab).

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [January 11, 2024, 10:45pm UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/15 "2024-01-11T22:45:47Z")

</div>

I’ve had mixed success with [ArnoldiMethod](https://github.com/JuliaLinearAlgebra/ArnoldiMethod.jl/issues/114). The modes were not at all mass-orthogonal. Still not sure why.

[KrylovKit](https://github.com/Jutho/KrylovKit.jl/issues/73) does not converge at all.

As far as I know Arpack is still the only package in julia that actually works.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [January 12, 2024, 2:12am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/16 "2024-01-12T02:12:15Z")

</div>

That isn’t clear to me. Arpack is still cited in the Matlab manual today.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/1/4/1401320c81879589d770d1129989563cf349bcca.png)

---

<div class="post-metadata">

**Author:** ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)\
**Post date:** [January 12, 2024, 2:44am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/17 "2024-01-12T02:44:40Z")

</div>

Here is my source:

> [@Some eigenpairs from a large, sparse, nonsymmetric matrix: Julia vs Matlab](https://discourse.julialang.org/t/some-eigenpairs-from-a-large-sparse-nonsymmetric-matrix-julia-vs-matlab/93742/18):
>
> Actually, Matlab does not use Arpack anymore. Instead, they have implemented their own version of the Krylov-Schur method accessible from the Matlab prompt with edit eigs.m. I am not sure whether the following tweak will resolve your issue, however, I remember from my own work on a modesolver that SuiteSparse conducts an iterative refinement of the solutions by default. Disabling this refinement using SuiteSparse.UMFPACK.umf\_ctrl[8] = 0 improved the performance of eigs significantly for my us…

Further down that thread is an updated way to get the faster sparse solution, which may help @jiang_ming_zhang as well.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [January 12, 2024, 3:21am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/18 "2024-01-12T03:21:12Z")

</div>

Interesting. Did you get the eigenvectors to converge with ArnoldiMethod? I did not.

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petersimon/32/25193_2.png) [@PeterSimon](https://discourse.julialang.org/u/PeterSimon)\
**Post date:** [January 12, 2024, 3:32am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/19 "2024-01-12T03:32:05Z")

</div>

> [@PetrKryslUCSD](#):
>
> Did you get the eigenvectors to converge with ArnoldiMethod?

Yes

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [January 12, 2024, 4:29am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/20 "2024-01-12T04:29:21Z")

</div>

Could you possibly provide the code and the data?

[Next page](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658.md?page=2)