# KrylovKit eigsolve of LinearMap giving different values in CPU and GPU

**URL:** <https://discourse.julialang.org/t/krylovkit-eigsolve-of-linearmap-giving-different-values-in-cpu-and-gpu/124032>\
**Category:** GPU\
**Created:** [December 20, 2024, 5:52pm UTC](https://discourse.julialang.org/t/krylovkit-eigsolve-of-linearmap-giving-different-values-in-cpu-and-gpu/124032 "2024-12-20T17:52:23Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![victor\_vhrn](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/victor_vhrn/32/214016_2.png) [@victor\_vhrn](https://discourse.julialang.org/u/victor_vhrn)\
**Post date:** [December 20, 2024, 5:52pm UTC](https://discourse.julialang.org/t/krylovkit-eigsolve-of-linearmap-giving-different-values-in-cpu-and-gpu/124032/1 "2024-12-20T17:52:23Z")

</div>

I am trying to compare the eigsolve function using CPU and GPU, using the following code:

```julia
using CUDA
using CUDA.CUSPARSE
using LinearAlgebra
using LinearMaps
using SparseArrays
using Random
using KrylovKit

function eign_cpu(M::Union{AbstractMatrix,LinearMap{T}}, v0::Vector{T}) where {T<:Number}
    tol = 1e-10
    (E, V, info) = eigsolve(M, v0, 1, :SR; ishermitian=true, tol=tol)
    return (E[1], V[1], info)
end

function eign_gpu(M::Union{CuArray,LinearMap{T}}, v0::CuVector{T}) where {T<:Number}
    tol = 1e-10
    (E, V, info) = eigsolve(M, v0, 1, :SR; ishermitian=true, tol=tol)
    return (E[1], V[1], info)
end

Random.seed!(0)

n = 4
m = 4
Es = [1, 1, 2, 3]
Ed = [2, 3, 4, 4]

w = ones(4)
sw = sqrt.(w)
sw_d = CuArray(sw)

C = LinearMap{Float64}((out::AbstractArray, y::AbstractArray) -> (out .= sw; out .*= -dot(sw, y)), n, ismutating=true, issymmetric=true)
C_d = LinearMap{Float64}((out::CuArray, y::CuArray) -> (out .= sw_d; out .*= -dot(sw_d, y)), n, ismutating=true, issymmetric=true)

M = sparse(Es, Ed, ones(Float64, m), n, n)
M_d = CuSparseMatrixCSC(M)

adj_cpu = p::AbstractArray -> (M.nzval .= p ./ 2; Symmetric(M))
adj_gpu = p::CuArray -> (M_d.nzVal .= p ./ 2; Symmetric(M_d))

p = rand(m, 1)
p_d = CuArray(p)

aux = Vector{Float64}(undef, n)
eig_grad = randn(Float64, n)
v0 = ((randn!(aux) .*= 1e-4) .+= eig_grad)
v0_d = CuArray(v0)

if isapprox(v0, Array(v0_d))
    println("v and v0 are approximate")
end
if isapprox(sw, Array(sw_d))
    println("C and C_d are approximate")
end
if isapprox(M, SparseMatrixCSC(M_d))
    println("M and M_d are approximate")
end
S = C + adj_cpu(p)
(E_cpu, V_cpu, info_cpu) = eign_cpu(S, v0)
println("info cpu = ", info_cpu)
S_d = C_d + adj_gpu(p_d)
(E_gpu, V_gpu, info_gpu) = eign_gpu(S_d, v0_d)
println("info gpu = ", info_gpu)

if !isapprox(V_cpu, Array(V_gpu))
    println("WARN: diverging eigenvalues between cpu and gpu:\n",
        "E_cpu = ", E_cpu, " V_cpu = ", V_cpu, "\nE_gpu = ", E_gpu, " V_gpu = ", V_gpu)
    readline()
end

```

The isapprox checks for the input data all return true. However, the output is:

```julia
info cpu = ConvergenceInfo: 4 converged values after 1 iterations and 4 applications of the linear map;
norms of residuals are given by (1.3683491030815672e-36, 6.545806401840613e-35, 6.344902430327517e-34, 8.827885896605726e-33).

info gpu = ConvergenceInfo: 4 converged values after 1 iterations and 4 applications of the linear map;
norms of residuals are given by (7.776628566574466e-20, 1.087879714604506e-17, 4.976482876540383e-17, 3.382716406082749e-17).

WARN: different eigenvalues between cpu and gpu:
E_cpu = -3.6566857184706016 V_cpu = [-0.5154937434198816, -0.46021338805521755, -0.5355038728938648, -0.485495046386023]
E_gpu = -3.7700145739756743 V_gpu = [0.48275625184457943, 0.47793647243109605, 0.52265534023327, 0.5151257370300364]

```

I expected the results to be approximately equal.  
Am I missing something in this experiment? Some wrong assumption maybe?

---

<div class="post-metadata">

**Author:** ![dubosipsl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dubosipsl/32/52998_2.png) [@dubosipsl](https://discourse.julialang.org/u/dubosipsl)\
**Post date:** [December 20, 2024, 11:31pm UTC](https://discourse.julialang.org/t/krylovkit-eigsolve-of-linearmap-giving-different-values-in-cpu-and-gpu/124032/2 "2024-12-20T23:31:50Z")

</div>

Just to be sure, did you check that `S` and `S_d` are the same ?  
Is there a guarantee that `M.nzval` and `M_d.nzval` refer to the same non-zero elements, in the same order ?

---

<div class="post-metadata">

**Author:** ![victor\_vhrn](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/victor_vhrn/32/214016_2.png) [@victor\_vhrn](https://discourse.julialang.org/u/victor_vhrn)\
**Post date:** [December 21, 2024, 1:22am UTC](https://discourse.julialang.org/t/krylovkit-eigsolve-of-linearmap-giving-different-values-in-cpu-and-gpu/124032/3 "2024-12-21T01:22:53Z")

</div>

When running:

```julia
println("M = ", M)
println("M_d = ", SparseMatrixCSC(M_d))

```

I get

```julia
M = sparse([1, 1, 2, 3], [2, 3, 4, 4], [1.0, 1.0, 1.0, 1.0], 4, 4)
M_d = sparse(Int32[1, 1, 2, 3], Int32[2, 3, 4, 4], [1.0, 1.0, 1.0, 1.0], 4, 4)

```

So the only difference is the type of the indexes, do you think it has some influence?

---

<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:** [December 21, 2024, 1:43am UTC](https://discourse.julialang.org/t/krylovkit-eigsolve-of-linearmap-giving-different-values-in-cpu-and-gpu/124032/4 "2024-12-21T01:43:34Z")

</div>

I think you have encountered a bug in the CUDA sparse array implementation: multiplication by `adj_gpu` is incomplete.  
There is [an open issue](https://github.com/JuliaGPU/CUDA.jl/issues/2572) for CUDA.jl

---

<div class="post-metadata">

**Author:** ![victor\_vhrn](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/victor_vhrn/32/214016_2.png) [@victor\_vhrn](https://discourse.julialang.org/u/victor_vhrn)\
**Post date:** [December 23, 2024, 3:03pm UTC](https://discourse.julialang.org/t/krylovkit-eigsolve-of-linearmap-giving-different-values-in-cpu-and-gpu/124032/5 "2024-12-23T15:03:38Z")

</div>

Thank you. That’s unfortunate.

I tried removing the Symmetric call from the adjs:

```julia
adj_cpu = p::AbstractArray -> (M.nzval .= p ./ 2; M)
adj_gpu = p::CuArray -> (M_d.nzVal .= p ./ 2; M_d)

```

And now the GPU and CPU returns the same minimum eigenvalue but with different eigenvectors:

```julia
WARN: diverging eigenvalues between cpu and gpu:
E_cpu = -3.770014573975674 V_cpu = [-0.4827562518445795, -0.47793647243109577, -0.5226553402332699, -0.5151257370300364]
E_gpu = -3.7700145739756743 V_gpu = [0.48275625184457943, 0.47793647243109605, 0.52265534023327, 0.5151257370300364]

```

However, The eigenvalue-eigenvector equality is not holding in both GPU and CPU:

```julia
if !isapprox(V_cpu, Array(V_gpu))
    println("WARN: diverging eigenvalues between cpu and gpu:\n",
        "E_cpu = ", E_cpu, " V_cpu = ", V_cpu, "\nE_gpu = ", E_gpu, " V_gpu = ", V_gpu)
    lhs_gpu = S_d * V_gpu
    rhs_gpu = E_gpu * V_gpu
    println("lhs gpu = ", lhs_gpu)
    println("rhs gpu = ", rhs_gpu)

    lhs_cpu = S * V_cpu
    rhs_cpu = E_cpu * V_cpu
    println("lhs cpu = ", lhs_cpu)
    println("rhs cpu = ", rhs_cpu)
    readline()
end

```

Output:

```julia
WARN: diverging eigenvalues between cpu and gpu:
E_cpu = -3.770014573975674 V_cpu = [-0.4827562518445795, -0.47793647243109577, -0.5226553402332699, -0.5151257370300364]
E_gpu = -3.7700145739756743 V_gpu = [0.48275625184457943, 0.47793647243109605, 0.52265534023327, 0.5151257370300364]
lhs gpu = [-1.8836115939631806, -1.7764183292957163, -1.976330898719767, -1.998473801538982]
rhs gpu = [-1.8199981051319354, -1.8018274664997551, -1.9704182498456424, -1.942031536033198]
lhs cpu = [1.8836115939631803, 1.7764183292957159, 1.9763308987197665, 1.9984738015389816]
rhs cpu = [1.8199981051319356, 1.8018274664997542, 1.970418249845642, 1.942031536033198]

```

Do you think this is also a bug or I am doing something wrong in my test?

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [December 23, 2024, 5:59pm UTC](https://discourse.julialang.org/t/krylovkit-eigsolve-of-linearmap-giving-different-values-in-cpu-and-gpu/124032/6 "2024-12-23T17:59:43Z")

</div>

> [@victor\_vhrn](#):
>
> but with different eigenvectors:

It is the same eigenvector. Just note that `E_cpu ≈ -E_gpu`, which is totally allowed (eigenvectors are generally only unique upto a factor and since KrylovKit normalizes the vector they are unique upto a global sign).

I cannot tell you why the eigenvalue/eigenvector pair seems to be off though. But again you can see that the values are almost identical between GPU and CPU (neglecting the global sign).

---

<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:** [December 23, 2024, 7:55pm UTC](https://discourse.julialang.org/t/krylovkit-eigsolve-of-linearmap-giving-different-values-in-cpu-and-gpu/124032/7 "2024-12-23T19:55:59Z")

</div>

Now the operator is nonsymmetric, but it seems that you told `eigsolve` otherwise.

---

<div class="post-metadata">

**Author:** ![victor\_vhrn](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/victor_vhrn/32/214016_2.png) [@victor\_vhrn](https://discourse.julialang.org/u/victor_vhrn)\
**Post date:** [December 23, 2024, 9:31pm UTC](https://discourse.julialang.org/t/krylovkit-eigsolve-of-linearmap-giving-different-values-in-cpu-and-gpu/124032/8 "2024-12-23T21:31:07Z")

</div>

That’s true, I forgot about that! By passing `issymmetric=false` to eigsolve, both eigenvalues and eigenvectors returned are equal.