Use StochasticRounding.jl ?

[milankl]

Just came across https://arxiv.org/abs/2403.12278 🎉 And I see that you're actually using Julia. Just to point you to StochasticRounding.jl which implements stochastic rounding for floating-point arithmetics using binary operations (and not floating-point arithmetic itself what seems what you're doing here) -- it's also a registered package. We don't implement stochastic for arbitrary precision (yet?) but only for Float16, BFloat16, Float32, Float64 but new types can easily be added. Feel free to reach out!

[cboutsikas]

Hi Milan, thank you for getting in touch and for sharing your package! I will definetely take a look at StochasticRounding.jl.

We will be in touch, thanks!

[milankl]

We implement stochastic rounding as new types <:AbstractFloat so that you can use many of the existing functionality in Julia out of the box, e.g.

julia> using StochasticRounding
julia> A = randn(Float32sr,3,3)
3×3 Matrix{Float32sr}:
 -1.09967  -0.049254  -0.87082
 -0.58354  -1.52668    0.123254
  1.78912   0.539212  -0.903346

julia> b = randn(Float32sr,3)
3-element Vector{Float32sr}:
  0.6314426f0
  0.29993263f0
 -0.7961567f0

julia> A\b
3-element Vector{Float32sr}:
 -0.49193418f0
 -0.016741086f0
 -0.10295086f0

The LU-decomposition isn't implemented by us, we only define the arithmetic operations but Julia compiles this to our new type (often called package composability). On that note, I realised that LinearAlgebra.svd escapes the type (start with Float32sr but ending up with Float32 so the stochastic rounding is lost somewhere on the way

julia> using LinearAlgebra
julia> svd(rand(Float32sr,2,2))
SVD{Float32, Float32, Matrix{Float32}, Vector{Float32}}
U factor:
2×2 Matrix{Float32}:
 -0.626842  -0.779147
 -0.779147   0.626842
singular values:
2-element Vector{Float32}:
 1.4279796
 0.35908842
Vt factor:
2×2 Matrix{Float32}:
 -0.806354  -0.591433
 -0.591433   0.806354

Using GenericLinearAlgebra.jl I just managed to get a fully (i.e. on every operation) stochastically rounded SVD to work

julia> using GenericLinearAlgebra
julia> A
3×3 Matrix{Float32sr}:
  1.8358      -1.09203      -0.000478083
 -8.40294f-6  -1.89813       0.968928
  1.5081f-7   -0.000401923  -0.668636

julia> GenericLinearAlgebra.svd!(A)
SVD{Float32sr, Float32sr, Matrix{Float32sr}, Vector{Float32sr}}
U factor:
3×3 Matrix{Float32sr}:
  0.701216   -0.706439  -0.0961228
  0.70912     0.677127   0.19659
 -0.0737917  -0.206015   0.975762
singular values:
3-element Vector{Float32sr}:
 2.5803685f0
 1.605163f0
 0.56269586f0
Vt factor:
3×3 Matrix{Float32sr}:
  0.498876  -0.81838    0.285266
 -0.807945  -0.320054   0.494762
 -0.313603  -0.477304  -0.820874

julia> GenericLinearAlgebra.svd!(A)
SVD{Float32sr, Float32sr, Matrix{Float32sr}, Vector{Float32sr}}
U factor:
3×3 Matrix{Float32sr}:
 -0.701203   -0.706458  -0.0960811
 -0.709097    0.677029   0.197008
  0.0741282  -0.206274   0.975683
singular values:
3-element Vector{Float32sr}:
 2.5804255f0
 1.6050962f0
 0.56262755f0
Vt factor:
3×3 Matrix{Float32sr}:
 0.498859   0.818239  -0.285701
 0.807995  -0.319839   0.49482
 0.313502  -0.47769   -0.820688

For some reason I had to define

julia> LinearAlgebra.floatmin2(::Type{Float32sr}) = reinterpret(Float32sr, 0x26000000)

👀

[cboutsikas]

Wow. The stochastically rounded SVD might be very helpful. I will definetely run some tests on that. Thanks!