# Question regarding numerical stability among Julia versions
**URL:** https://discourse.julialang.org/t/question-regarding-numerical-stability-among-julia-versions/112741
**Category:** Numerics
**Tags:** question
**Created:** [April 9, 2024, 5:24pm UTC](https://discourse.julialang.org/t/question-regarding-numerical-stability-among-julia-versions/112741 "2024-04-09T17:24:06Z")
**Posts on this page:** 17
**Page:** 1
### Author:  [@CodeLenz](https://discourse.julialang.org/u/CodeLenz)
#### Post date: [April 9, 2024, 5:24pm UTC](https://discourse.julialang.org/t/question-regarding-numerical-stability-among-julia-versions/112741/1 "2024-04-09T17:24:06Z")
Dear Julians.
I was testing the following code (M,C and K are sparse, symmetric and positive definite matrices)
> Julia\> M,C,K=THS();
>
> julia\> cholM=cholesky(M);
>
> julia\> Cb = cholM\C;
>
> julia\> Kb = cholM\K;
>
> julia\> F = 0.5\*(Array(Cb) + sqrt(Array(Cb^2 - 4\*Kb)));
>
> julia\> norm(F^2 - F\*Cb + Kb)
Using Julia-dev Version 1.12.0-DEV.317 Commit 0e28cf6abf (2024-04-08 10:46 UTC)
or Julia Version 1.11.0-alpha2 Commit 9dfd28ab751 (2024-03-18 20:35 UTC)
gives a reasonable value
> 1.0530104604285245e-5
but using julia 1.9 (installed using juliaup in a Linux (Fedora 39) machine gives
> 7.606244293550662e52
a very weird result.
More detailed data about the versions
> julia\> versioninfo()
> Julia Version 1.12.0-DEV.317
> Commit 0e28cf6abf (2024-04-08 10:46 UTC)
> Platform Info:
> OS: Linux (x86\_64-redhat-linux)
> CPU: 12 × 12th Gen Intel(R) Core™ i5-1235U
> WORD\_SIZE: 64
> LLVM: libLLVM-16.0.6 (ORCJIT, alderlake)
> Threads: 1 default, 0 interactive, 1 GC (on 12 virtual cores)
> Environment:
> JULIA\_EDITOR = code
and
> julia\> versioninfo()
> Julia Version 1.9.4
> Commit 8e5136fa297 (2023-11-14 08:46 UTC)
> Build Info:
> Official [https://julialang.org/](https://julialang.org/) release
> Platform Info:
> OS: Linux (x86\_64-linux-gnu)
> CPU: 12 × 12th Gen Intel(R) Core™ i5-1235U
> WORD\_SIZE: 64
> LIBM: libopenlibm
> LLVM: libLLVM-14.0.6 (ORCJIT, alderlake)
> Threads: 1 on 12 virtual cores
> Environment:
> JULIA\_EDITOR = code
It does not happen for small, dense, matrices. So I guess is something related to numerical stability of computing the sqrt.
Matrices are 2200 x 2200
K
> [K.txt · GitHub](https://gist.github.com/CodeLenz/fdf69c3b50d90a4a2d15f2bae13a400e)
C = 1E-6\*K
M
> [Matrix M · GitHub](https://gist.github.com/CodeLenz/a667ee800746c015f7a172c06cd96934)
Does anyone can give me some light regarding what may be the cause for this difference?
Thank you!
---
### Author:  [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [April 9, 2024, 11:27pm UTC](https://discourse.julialang.org/t/question-regarding-numerical-stability-among-julia-versions/112741/2 "2024-04-09T23:27:33Z")
Can you identify which operation causes the difference? It may be an improvement (bug fix, or algorithmic) on some operation, or just a case of a very unstable calculation depending on numerical inaccuracies.
What happens in 1.10?
---
### Author:  [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)
#### Post date: [April 10, 2024, 4:00am UTC](https://discourse.julialang.org/t/question-regarding-numerical-stability-among-julia-versions/112741/3 "2024-04-10T04:00:58Z")
@sethaxen revised the matrix `sqrt` code in that time frame; perhaps he eliminated a bug in the process?
[https://github.com/JuliaLang/julia/pull/39973](https://github.com/JuliaLang/julia/pull/39973)
---
### Author:  [@CodeLenz](https://discourse.julialang.org/u/CodeLenz)
#### Post date: [April 10, 2024, 12:05pm UTC](https://discourse.julialang.org/t/question-regarding-numerical-stability-among-julia-versions/112741/4 "2024-04-10T12:05:48Z")
@lmiq and @Ralph_Smith Thank you for your advice. I will take a look at 1.10 and the result of sqrt (as well as the conditioning of both matrices).
---
### Author:  [@CodeLenz](https://discourse.julialang.org/u/CodeLenz)
#### Post date: [April 10, 2024, 6:46pm UTC](https://discourse.julialang.org/t/question-regarding-numerical-stability-among-julia-versions/112741/5 "2024-04-10T18:46:21Z")
It seems to be an overflow problem during the sqrt(A)
In this example, A is non symmetric with norm=1.32E9, smaller eigenvalue is
-4.8e7 and the largest is -4.06 (large condition number). The resulting matrix of sqrt(A) has norm **8.58E39** and does not satisfy (sqrt(A))^2 - A = 0
Actually, the norm of the difference is
> julia\> norm( (sqrt(A)^2 - A)
> 2.853968508045082e63
Indeed, if I scale A by its norm A2 = A/norm(A) it works such that
> julia\> norm( (sqrt(A2)^2 - A2)
> 3.185526215409009e-14
With this example, I observe the same behavior in Julia 1.9, 1.10, 1.11 and dev.
---
### Author:  [@CodeLenz](https://discourse.julialang.org/u/CodeLenz)
#### Post date: [April 11, 2024, 12:37am UTC](https://discourse.julialang.org/t/question-regarding-numerical-stability-among-julia-versions/112741/7 "2024-04-11T00:37:23Z")
Using eigen decomposition works
> julia\> D,U = eigen(A);
>
> julia\> R=sqrt.(Complex.(D));
>
> julia\> norm(U\*diagm(R)\*diagm(R)\*inv(U)-A)
> 9.41423450406504e-5
---
### Author:  [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)
#### Post date: [April 11, 2024, 5:00am UTC](https://discourse.julialang.org/t/question-regarding-numerical-stability-among-julia-versions/112741/8 "2024-04-11T05:00:47Z")
Thanks, but why the overflow in the original post materialize only on certain Julia versions?
---
### Author:  [@sethaxen](https://discourse.julialang.org/u/sethaxen)
#### Post date: [April 11, 2024, 8:40am UTC](https://discourse.julialang.org/t/question-regarding-numerical-stability-among-julia-versions/112741/9 "2024-04-11T08:40:49Z")
> [@Ralph\_Smith](#):
>
> @sethaxen revised the matrix `sqrt` code in that time frame; perhaps he eliminated a bug in the process?
> [https://github.com/JuliaLang/julia/pull/39973](https://github.com/JuliaLang/julia/pull/39973)
No, that wouldn’t be it. That PR was present already in v1.7.0.
---
### Author:  [@sethaxen](https://discourse.julialang.org/u/sethaxen)
#### Post date: [April 11, 2024, 8:42am UTC](https://discourse.julialang.org/t/question-regarding-numerical-stability-among-julia-versions/112741/10 "2024-04-11T08:42:03Z")
@CodeLenz you could try git bisecting a local build of Julia to identify the commit at which the instability goes away.
---
### Author:  [@CodeLenz](https://discourse.julialang.org/u/CodeLenz)
#### Post date: [April 11, 2024, 2:06pm UTC](https://discourse.julialang.org/t/question-regarding-numerical-stability-among-julia-versions/112741/11 "2024-04-11T14:06:52Z")
OK guys. I noticed that I was loading MKL.jl in my startup.jl. I also simplified the problem avoiding the computation of Cb. Thus, matrix A now is just -4\*Kb.
All the eigenvalues of A are negative and inv(A)\*A \approx I, as expected. A is non symmetric with dimension 2200 x 2200.
I am getting very weird results for this “simple” problem. ~~After inspecting the problem, it seems that the culprit is the Schur decomposition used in sqrt AFTER 1.10.2.~~
There is also a problem with MKL.jl
Test script:
> cholM = cholesky(M);
>
> Kb = Array(cholM\K)
>
> A = - 4\*Kb
>
> julia\_root = sqrt(A)
>
> diff = norm(julia\_root\*julia\_root .- A)
>
> if diff\>1
> println("failure ",diff)
> else
> println("Success ",diff)
> end
Results
1. Julia 1.0.4 → Success (norm is 1E-5)
Julia 1.0.4 + MKL → Failure (norm is 3.7E64)
2. Julia 1.10.2 → Success (norm is 4E-5)
Julia 1.10.2 + MKL → Failure (norm is 3.7E64)
3. Julia 1.11.0-beta1 → Failure (norm is 1.4E61)
Julia 1.11.0-beta1 + MKL → Failure (norm is 2.8E63)
4. Julia DEV.322 → Failure (norm is 1E61)
Julia DEV.322 + MKL → Failure (norm is 2.8E63)
Thus, there is some problem with MKL.jl (at least, regarding sqrt of a negative definite, non symmetric matrix like A)
Now, testing with eigen decomposition of A (it should work!)
> # Eigen decomposition
>
> D,U = eigen(A)
>
> R=sqrt.(Complex.(D))
>
> diff\_eigen = norm(U\*diagm(R)\*diagm(R)\*inv(U) .- A)
>
> if diff\_eigen\>1
> println("failure ",diff\_eigen)
> else
> println("Success ",diff\_eigen)
> end
Results
1. Julia 1.0.4 → Success (norm is 7.7E-5)
Julia 1.0.4 + MKL → Success (norm is 0.0002858379704314278)
2. Julia 1.10.2 → Success (norm is 0.0001034042666795483)
Julia 1.10.2 + MKL → Success (norm is 0.0002858379704314278)
3. Julia 1.11.0-beta1 → Success (norm is 9.41423450406504e-5)
Julia 1.11.0-beta1 + MKL → Success (norm is 0.0006612749077237161)
4. Julia DEV.322 → Success (norm is 9.41423450406504e-5)
Julia DEV.322 + MKL → Success (norm is 0.0006612749077237161)
I understand that sqrt(A) is using Schur decomposition. Lets give a try
> # Schur decomposition
>
> S = Schur{Complex}(schur(A))
>
> R = sqrt(Complex.(S.Schur))
>
> diff\_schur = norm( (S.Z)_(R_R)\*adjoint(S.Z) .- A)
>
> if diff\_schur\>1
> println("failure ",diff\_schur)
> else
> println("Success ",diff\_schur)
> end
1. Julia 1.0.4 → Sucess
Julia 1.0.4 + MKL → Failure
2. Julia 1.10.2 → Success
Julia 1.10.2 + MKL → Failure
3. Julia 1.11.0-beta1 → Failure
Julia 1.11.0-beta1 + MKL → Failure
4. Julia DEV.322 → Failure
Julia DEV.322 + MKL → Failure
EDIT
Additional test
> S = Schur{Complex}(schur(A))
>
> R = sqrt(Complex.(S.Schur))
>
> A2 = (S.Z)\*(S.Schur)\*adjoint(S.Z)
>
> @show norm(A2.-A)
and all versions, with and without MKL are giving a correct matrix A2. Thus, the problem is not in the Schur decomposition. I also tested sqrt(Complex.(A)) with the same results as before.
Thus, the problem should be in sqrt(S.Schur) an upper triangular complex matrix.
---
### Author:  [@CodeLenz](https://discourse.julialang.org/u/CodeLenz)
#### Post date: [April 11, 2024, 3:51pm UTC](https://discourse.julialang.org/t/question-regarding-numerical-stability-among-julia-versions/112741/12 "2024-04-11T15:51:18Z")
Dear @sethaxen
I would have to study this procedure, since I never did something similar. I am trying to understand the problem to pinpoit the culprit.
---
### Author:  [@CodeLenz](https://discourse.julialang.org/u/CodeLenz)
#### Post date: [April 12, 2024, 12:07am UTC](https://discourse.julialang.org/t/question-regarding-numerical-stability-among-julia-versions/112741/13 "2024-04-12T00:07:35Z")
Julia version for the algorithm in
[https://mathsfromnothing.au/matrix-square-root/?i=1](https://mathsfromnothing.au/matrix-square-root/?i=1)
works in every version I tested, with and without MKL
```julia
function Sqrt(A)
# Dimension
n = size(A,1)
# Schur decomposition of A
S,U = Schur{Complex}(schur(A))
# Diagonal of S
D = diag(S)
# sqrt of main diagonal of S
X = diagm(sqrt.(D))
if !isdiag(S)
for j=2:n
for i=j-1:-1:1
k = i+1:j-1
somat = transpose(X[i,k])*X[k,j]
X[i,j] = (S[i,j]-somat)/(X[i,i]+X[j,j])
end
end
end
sqrtA = U*X*adjoint(U)
end
```
---
### Author:  [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)
#### Post date: [April 12, 2024, 2:00am UTC](https://discourse.julialang.org/t/question-regarding-numerical-stability-among-julia-versions/112741/14 "2024-04-12T02:00:54Z")
Do your matrices have any repeated eigenvalues? The matrix `sqrt` code uses a blocked scheme for large nonsymmetric matrices which invokes a Sylvester solver that is unreliable if multiple eigenvalues show up in the same block. For some reason the Julia code ignores the error return from LAPACK in such cases.
The simpler approaches, like the one you just listed, do not suffer from this problem.
I say “unreliable” rather than “broken” because whether the solver actually fails depends on how close the eigenvalues appear to be and some floating point details, so that could explain why a given matrix behaves differently with different libraries.
---
### Author:  [@CodeLenz](https://discourse.julialang.org/u/CodeLenz)
#### Post date: [April 12, 2024, 10:11am UTC](https://discourse.julialang.org/t/question-regarding-numerical-stability-among-julia-versions/112741/15 "2024-04-12T10:11:39Z")
@Ralph_Smith Indeed, there are many repeated eigenvalues. Thank your for your feedback.
Nonetheless, sqrt(schur(A).Schur) is the main culprit here and this matrix is UpperTriangular. I was thinking that Julia dispatches to a specific solver for triu matrices, using Lapack. Other interesting fact is the MKL error in all versions I tested.
Thank you.
---
### Author:  [@CodeLenz](https://discourse.julialang.org/u/CodeLenz)
#### Post date: [April 12, 2024, 11:06am UTC](https://discourse.julialang.org/t/question-regarding-numerical-stability-among-julia-versions/112741/16 "2024-04-12T11:06:25Z")
@sethaxen
Something related to sqrt(A) after (not including) 1.10.2. As @Ralph_Smith said, Julia is ignoring the return error flag from LAPACK in this case (A has repeated eigenvalues).
I have opened a ticket with this issue
>
>
> Reference thread
>
> \[(https://discourse.julialang.org/t/question-regarding-nume…rical-stability-among-julia-versions/112741)\]
>
> sqrt(A) is failing when A is non-symmetric with repeated eigenvalues.
---
### Author:  [@CodeLenz](https://discourse.julialang.org/u/CodeLenz)
#### Post date: [April 15, 2024, 9:51am UTC](https://discourse.julialang.org/t/question-regarding-numerical-stability-among-julia-versions/112741/17 "2024-04-15T09:51:07Z")
I made a very small and simple repository to this function (with some checks and obvious optimizations). Just in case someone is interested in this solution.
> **[GitHub - CodeLenz/MySqrt: Hack to bypass julias sqrt(A)](https://github.com/CodeLenz/MySqrt)**
>
> Hack to bypass julias sqrt(A) . Contribute to CodeLenz/MySqrt development by creating an account on GitHub.
---
### Author:  [@CodeLenz](https://discourse.julialang.org/u/CodeLenz)
#### Post date: [April 16, 2024, 12:54pm UTC](https://discourse.julialang.org/t/question-regarding-numerical-stability-among-julia-versions/112741/18 "2024-04-16T12:54:05Z")
Update:
`A[abs.(A).