Error in Sqrt for julia > 1.10.2

CodeLenz commented 5 months ago

Reference thread

https://discourse.julialang.org/t/question-regarding-numerical-stability-among-julia-versions/112741

sqrt(A) is failing when A is non-symmetric with repeated eigenvalues.

CodeLenz commented 4 months ago

The culprit seems to be schur(A), computed inside sqrt(A).

dkarrasch commented 4 months ago

@aravindh-krishnamoorthy You may be interested in this. I saw you're working on efficient square roots of matrices.

aravindh-krishnamoorthy commented 4 months ago

Interesting problem! I see from the Discourse thread that @CodeLenz has studied it quite a bit. I'll take a closer look at the issue.

CodeLenz commented 4 months ago

@aravindh-krishnamoorthy

I tried the code in your repository, but the same pattern remains. I managed to reduce the test Matrix to

function Test_sqrt(N=1000)

    # Initialize matrix A
    A = zeros(N,N)

    # A is tridiagonal
    E = 210E9
    p = 7.3E-4
    l = 20.0/N

    # Main diagonal
    cte = E/(p*l^2)
    A[diagind(A,0)].= 2.0*cte

    # Upper and lower first diagonals
    cte = -E/(p*l^2)
    A[diagind(A,1)].= cte
    A[diagind(A,-1)].= cte

    # Sqrt
    sA = sqrt(A)

    # Test sqrt
    return A, norm(sA*sA .- A) 

end

aravindh-krishnamoorthy commented 4 months ago

@aravindh-krishnamoorthy

I tried the code in your repository, but the same pattern remains. I managed to reduce the test Matrix to

[snip|

Hello @CodeLenz. I just started with your function Test_sqrt provided above. But, as seen from below, there's no difference in the error values produced between versions 1.9.4 and 1.11.0-alpha2 (which is > 1.10.2) on my PC. Could you please help me identify the exact problem to be solved? Do you think it's better for me to look at the original problem reported in the Discourse thread with the square root of A = -4*Kb having an issue between these versions?

julia> VERSION
v"1.9.4"
julia> using LinearAlgebra
julia> include("test_sqrt.jl")
Test_sqrt (generic function with 2 methods)
julia> A, e = Test_sqrt()
([1.4383561643835615e18 -7.191780821917807e17 … 0.0 0.0; -7.191780821917807e17 1.4383561643835615e18 … 0.0 0.0; … ; 0.0 0.0 … 1.4383561643835615e18 -7.191780821917807e17; 0.0 0.0 … -7.191780821917807e17 1.4383561643835615e18], 5.647486380061708e6)

julia> VERSION
v"1.11.0-alpha2"
julia> using LinearAlgebra
julia> include("test_sqrt.jl")
Test_sqrt (generic function with 2 methods)
julia> A, e = Test_sqrt()
([1.4383561643835615e18 -7.191780821917807e17 … 0.0 0.0; -7.191780821917807e17 1.4383561643835615e18 … 0.0 0.0; … ; 0.0 0.0 … 1.4383561643835615e18 -7.191780821917807e17; 0.0 0.0 … -7.191780821917807e17 1.4383561643835615e18], 5.647486380061708e6)

JuliaLang / julia

Error in Sqrt for julia > 1.10.2 #54062