# 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:** 2

<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, 5:08am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/21 "2024-01-12T05:08:49Z")

</div>

The code is below, reproduced from [this post](https://discourse.julialang.org/t/some-eigenpairs-from-a-large-sparse-nonsymmetric-matrix-julia-vs-matlab/93742/14):

```julia
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`.

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [January 12, 2024, 10:42am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/22 "2024-01-12T10:42:08Z")

</div>

There’s this [question](https://www.mathworks.com/matlabcentral/answers/1795480-what-is-matlab-doing-under-the-hood-when-i-solve-this-generalized-eigenvalue-problem) 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.

---

<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, 1:29pm UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/23 "2024-01-12T13:29:45Z")

</div>

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 ||A_v - B_v\*Lam|| is large. So, no, ArnoldiMethod did not succeed.

---

<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:12pm UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/24 "2024-01-12T14:12:01Z")

</div>

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

---

<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 13, 2024, 12:23pm UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/26 "2024-01-13T12:23:08Z")

</div>

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.

```julia
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

```julia
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

```

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [January 13, 2024, 12:59pm UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/27 "2024-01-13T12:59:04Z")

</div>

> [@jiang\_ming\_zhang](#):
>
> ```julia
> 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
> 
> ```

I find that if I replace this function with

```julia
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.

---

<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 13, 2024, 5:04pm UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/28 "2024-01-13T17:04:11Z")

</div>

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

```julia
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:

```julia
    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?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 13, 2024, 6:21pm UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/29 "2024-01-13T18:21:56Z")

</div>

> [@PetrKryslUCSD](#):
>
> 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?

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?

---

<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 13, 2024, 6:35pm UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/30 "2024-01-13T18:35:05Z")

</div>

> [@stevengj](#):
>
> It’s also not clear what is “not working” — the residual norms are small, so what is incorrect about the answers?

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.

---

<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 13, 2024, 6:46pm UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/31 "2024-01-13T18:46:28Z")

</div>

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:

```julia
=== 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.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [January 13, 2024, 6:50pm UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/32 "2024-01-13T18:50:32Z")

</div>

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?

---

<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 13, 2024, 6:51pm UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/33 "2024-01-13T18:51:42Z")

</div>

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

[https://github.com/PetrKryslUCSD/GEPHelpers.jl:](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._

---

<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 13, 2024, 6:58pm UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/34 "2024-01-13T18:58:19Z")

</div>

> [@stevengj](#):
>
> It’s hard to give more concrete advice without knowing your matrices.

> <https://github.com/PetrKryslUCSD/GEPHelpers.jl/tree/main/test>
>
> //github.com/PetrKryslUCSD/GEPHelpers.jl/tree/main/test

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [January 13, 2024, 7:50pm UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/35 "2024-01-13T19:50:09Z")

</div>

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:

```julia
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.

---

<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 13, 2024, 7:55pm UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/36 "2024-01-13T19:55:14Z")

</div>

The norm of M\K is ~4e4.

---

<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 13, 2024, 8:20pm UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/37 "2024-01-13T20:20:42Z")

</div>

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);`):

```julia
=== 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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 13, 2024, 9:51pm UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/38 "2024-01-13T21:51:57Z")

</div>

> [@PetrKryslUCSD](#):
>
> `Computing 20 eigenvalues of type 'smallestabs'.`

From the Matlab `eigs` documentation:

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

This is what I [suggested above](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/29): 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:

```julia
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`).

---

<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 13, 2024, 11:38pm UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/39 "2024-01-13T23:38:51Z")

</div>

I think you meant:

```julia
    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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 14, 2024, 1:47am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/40 "2024-01-14T01:47:33Z")

</div>

> [@PetrKryslUCSD](#):
>
> I think you meant:

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

---

<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 14, 2024, 1:50am UTC](https://discourse.julialang.org/t/julia-is-slower-than-matlab-when-it-comes-to-matrix-diagonalization/108658/41 "2024-01-14T01:50:08Z")

</div>

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

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