BigFloat compatibility for qr_jacobimatrix #146

TSGut · 2023-07-26T18:23:04Z

@dlfivefifty @MikaelSlevinsky This PR restores BigFloat compatibility to the QR variants of jacobimatrix computation. The fix in theory also addresses Cholesky but that remains broken because of a similar issue in InfiniteLinearAlgebra.jl exemplified by

julia> P = Normalized(legendre(big(0)..big(1)))
Normalized(Legendre{BigFloat}() affine mapped to 0..1)

julia> X = jacobimatrix(P)
ℵ₀×ℵ₀ LazyBandedMatrices.SymTridiagonal{BigFloat, Fill{BigFloat, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, BroadcastVector{BigFloat, typeof(sqrt), Tuple{BroadcastVector{BigFloat, typeof(*), Tuple{BroadcastVector{BigFloat, typeof(/), Tuple{BroadcastVector{BigFloat, typeof(/), Tuple{InfiniteArrays.InfStepRange{BigFloat, BigFloat}, InfiniteArrays.InfStepRange{Int64, Int64}}}, BigFloat}}, BroadcastVector{BigFloat, typeof(/), Tuple{BroadcastVector{BigFloat, typeof(/), Tuple{InfiniteArrays.InfStepRange{BigFloat, BigFloat}, InfiniteArrays.InfStepRange{Int64, Int64}}}, BigFloat}}}}}}} with indices OneToInf()×OneToInf():
 0.5       0.288675   ⋅         ⋅         ⋅         ⋅         ⋅         ⋅       …  
 0.288675  0.5       0.258199   ⋅         ⋅         ⋅         ⋅         ⋅          
  ⋅        0.258199  0.5       0.253546   ⋅         ⋅         ⋅         ⋅          
  ⋅         ⋅        0.253546  0.5       0.251976   ⋅         ⋅         ⋅          
  ⋅         ⋅         ⋅        0.251976  0.5       0.251259   ⋅         ⋅          
  ⋅         ⋅         ⋅         ⋅        0.251259  0.5       0.250873   ⋅       …  
  ⋅         ⋅         ⋅         ⋅         ⋅        0.250873  0.5       0.25064     
  ⋅         ⋅         ⋅         ⋅         ⋅         ⋅        0.25064   0.5         
  ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅        0.25049     
 ⋮                                                 ⋮                            ⋱  

julia> cholesky(Symmetric(X^2)).U[1,1]
ERROR: UndefRefError: access to undefined reference
Stacktrace:
  [1] getindex
    @ ./essentials.jl:14 [inlined]
  [2] getindex
    @ ./subarray.jl:286 [inlined]
  [3] _getindex
    @ ./abstractarray.jl:1344 [inlined]
  [4] getindex
    @ ./abstractarray.jl:1294 [inlined]
  [5] macro expansion
    @ ~/Documents/Backends/julia-1.9.0/share/julia/stdlib/v1.9/LinearAlgebra/src/generic.jl:221 [inlined]
  [6] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [7] lmul!(s::BigFloat, X::SubArray{BigFloat, 2, Matrix{BigFloat}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false})
    @ LinearAlgebra ~/Documents/Backends/julia-1.9.0/share/julia/stdlib/v1.9/LinearAlgebra/src/generic.jl:220
  [8] _banded_lmul!
    @ ~/.julia/packages/BandedMatrices/qLIMl/src/generic/broadcast.jl:938 [inlined]
  [9] banded_lmul!
    @ ~/.julia/packages/BandedMatrices/qLIMl/src/generic/broadcast.jl:947 [inlined]
 [10] materialize!
    @ ~/.julia/packages/BandedMatrices/qLIMl/src/generic/broadcast.jl:950 [inlined]
 [11] #lmul!#10
    @ ~/.julia/packages/ArrayLayouts/3K92Y/src/lmul.jl:47 [inlined]
 [12] lmul!
    @ ~/.julia/packages/ArrayLayouts/3K92Y/src/lmul.jl:47 [inlined]
 [13] default_blasmul!(α::BigFloat, A::Adjoint{BigFloat, SubArray{BigFloat, 2, BandedMatrix{BigFloat, Matrix{BigFloat}, Base.OneTo{Int64}}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}}, B::SubArray{BigFloat, 2, BandedMatrix{BigFloat, Matrix{BigFloat}, Base.OneTo{Int64}}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}, β::BigFloat, C::SubArray{BigFloat, 2, BandedMatrix{BigFloat, Matrix{BigFloat}, Base.OneTo{Int64}}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/3K92Y/src/muladd.jl:172
 [14] materialize!(M::ArrayLayouts.MulAdd{BandedMatrices.BandedRows{ArrayLayouts.DenseColumnMajor}, BandedMatrices.BandedColumns{ArrayLayouts.DenseColumnMajor}, BandedMatrices.BandedColumns{ArrayLayouts.DenseColumnMajor}, BigFloat, Adjoint{BigFloat, SubArray{BigFloat, 2, BandedMatrix{BigFloat, Matrix{BigFloat}, Base.OneTo{Int64}}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}}, SubArray{BigFloat, 2, BandedMatrix{BigFloat, Matrix{BigFloat}, Base.OneTo{Int64}}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}, SubArray{BigFloat, 2, BandedMatrix{BigFloat, Matrix{BigFloat}, Base.OneTo{Int64}}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/3K92Y/src/muladd.jl:241
 [15] muladd!
    @ ~/.julia/packages/ArrayLayouts/3K92Y/src/muladd.jl:70 [inlined]
 [16] partialcholesky!(F::InfiniteLinearAlgebra.AdaptiveCholeskyFactors{BigFloat, BandedMatrix{BigFloat, Matrix{BigFloat}, Base.OneTo{Int64}}, ApplyArray{BigFloat, 2, typeof(*), Tuple{LazyBandedMatrices.SymTridiagonal{BigFloat, Fill{BigFloat, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, BroadcastVector{BigFloat, typeof(sqrt), Tuple{BroadcastVector{BigFloat, typeof(*), Tuple{BroadcastVector{BigFloat, typeof(/), Tuple{BroadcastVector{BigFloat, typeof(/), Tuple{InfiniteArrays.InfStepRange{BigFloat, BigFloat}, InfiniteArrays.InfStepRange{Int64, Int64}}}, BigFloat}}, BroadcastVector{BigFloat, typeof(/), Tuple{BroadcastVector{BigFloat, typeof(/), Tuple{InfiniteArrays.InfStepRange{BigFloat, BigFloat}, InfiniteArrays.InfStepRange{Int64, Int64}}}, BigFloat}}}}}}}, LazyBandedMatrices.SymTridiagonal{BigFloat, Fill{BigFloat, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, BroadcastVector{BigFloat, typeof(sqrt), Tuple{BroadcastVector{BigFloat, typeof(*), Tuple{BroadcastVector{BigFloat, typeof(/), Tuple{BroadcastVector{BigFloat, typeof(/), Tuple{InfiniteArrays.InfStepRange{BigFloat, BigFloat}, InfiniteArrays.InfStepRange{Int64, Int64}}}, BigFloat}}, BroadcastVector{BigFloat, typeof(/), Tuple{BroadcastVector{BigFloat, typeof(/), Tuple{InfiniteArrays.InfStepRange{BigFloat, BigFloat}, InfiniteArrays.InfStepRange{Int64, Int64}}}, BigFloat}}}}}}}}}}, n::Int64)
    @ InfiniteLinearAlgebra ~/.julia/packages/InfiniteLinearAlgebra/y1GuF/src/infcholesky.jl:35
 [17] getindex(F::InfiniteLinearAlgebra.AdaptiveCholeskyFactors{BigFloat, BandedMatrix{BigFloat, Matrix{BigFloat}, Base.OneTo{Int64}}, ApplyArray{BigFloat, 2, typeof(*), Tuple{LazyBandedMatrices.SymTridiagonal{BigFloat, Fill{BigFloat, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, BroadcastVector{BigFloat, typeof(sqrt), Tuple{BroadcastVector{BigFloat, typeof(*), Tuple{BroadcastVector{BigFloat, typeof(/), Tuple{BroadcastVector{BigFloat, typeof(/), Tuple{InfiniteArrays.InfStepRange{BigFloat, BigFloat}, InfiniteArrays.InfStepRange{Int64, Int64}}}, BigFloat}}, BroadcastVector{BigFloat, typeof(/), Tuple{BroadcastVector{BigFloat, typeof(/), Tuple{InfiniteArrays.InfStepRange{BigFloat, BigFloat}, InfiniteArrays.InfStepRange{Int64, Int64}}}, BigFloat}}}}}}}, LazyBandedMatrices.SymTridiagonal{BigFloat, Fill{BigFloat, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, BroadcastVector{BigFloat, typeof(sqrt), Tuple{BroadcastVector{BigFloat, typeof(*), Tuple{BroadcastVector{BigFloat, typeof(/), Tuple{BroadcastVector{BigFloat, typeof(/), Tuple{InfiniteArrays.InfStepRange{BigFloat, BigFloat}, InfiniteArrays.InfStepRange{Int64, Int64}}}, BigFloat}}, BroadcastVector{BigFloat, typeof(/), Tuple{BroadcastVector{BigFloat, typeof(/), Tuple{InfiniteArrays.InfStepRange{BigFloat, BigFloat}, InfiniteArrays.InfStepRange{Int64, Int64}}}, BigFloat}}}}}}}}}}, k::Int64, j::Int64)
    @ InfiniteLinearAlgebra ~/.julia/packages/InfiniteLinearAlgebra/y1GuF/src/infcholesky.jl:42
 [18] getindex(A::UpperTriangular{BigFloat, InfiniteLinearAlgebra.AdaptiveCholeskyFactors{BigFloat, BandedMatrix{BigFloat, Matrix{BigFloat}, Base.OneTo{Int64}}, ApplyArray{BigFloat, 2, typeof(*), Tuple{LazyBandedMatrices.SymTridiagonal{BigFloat, Fill{BigFloat, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, BroadcastVector{BigFloat, typeof(sqrt), Tuple{BroadcastVector{BigFloat, typeof(*), Tuple{BroadcastVector{BigFloat, typeof(/), Tuple{BroadcastVector{BigFloat, typeof(/), Tuple{InfiniteArrays.InfStepRange{BigFloat, BigFloat}, InfiniteArrays.InfStepRange{Int64, Int64}}}, BigFloat}}, BroadcastVector{BigFloat, typeof(/), Tuple{BroadcastVector{BigFloat, typeof(/), Tuple{InfiniteArrays.InfStepRange{BigFloat, BigFloat}, InfiniteArrays.InfStepRange{Int64, Int64}}}, BigFloat}}}}}}}, LazyBandedMatrices.SymTridiagonal{BigFloat, Fill{BigFloat, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, BroadcastVector{BigFloat, typeof(sqrt), Tuple{BroadcastVector{BigFloat, typeof(*), Tuple{BroadcastVector{BigFloat, typeof(/), Tuple{BroadcastVector{BigFloat, typeof(/), Tuple{InfiniteArrays.InfStepRange{BigFloat, BigFloat}, InfiniteArrays.InfStepRange{Int64, Int64}}}, BigFloat}}, BroadcastVector{BigFloat, typeof(/), Tuple{BroadcastVector{BigFloat, typeof(/), Tuple{InfiniteArrays.InfStepRange{BigFloat, BigFloat}, InfiniteArrays.InfStepRange{Int64, Int64}}}, BigFloat}}}}}}}}}}}, i::Int64, j::Int64)
    @ LinearAlgebra ~/Documents/Backends/julia-1.9.0/share/julia/stdlib/v1.9/LinearAlgebra/src/triangular.jl:232
 [19] top-level scope

Once that is fixed we can remove the _broken from the test added here.

codecov · 2023-07-26T18:39:16Z

Codecov Report

Patch coverage: 100.00% and project coverage change: +0.01% 🎉

Comparison is base (aaa182b) 91.76% compared to head (729870c) 91.78%.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #146      +/-   ##
==========================================
+ Coverage   91.76%   91.78%   +0.01%     
==========================================
  Files          17       17              
  Lines        1822     1825       +3     
==========================================
+ Hits         1672     1675       +3     
  Misses        150      150

Files Changed	Coverage Δ
src/choleskyQR.jl	`100.00% <100.00%> (ø)`

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

TSGut · 2023-07-26T20:18:41Z

Ok @dlfivefifty , after some digging: the error above for Cholesky appears to live in this line https://github.com/JuliaLinearAlgebra/InfiniteLinearAlgebra.jl/blob/0fc35261c6e8a8ce2bb4b647cbf935855e4b81ec/src/infcholesky.jl#L35

ArrayLayouts fails to materialize MulAdd for certain BigFloat matrices. I can avoid this by special casing the line including muladd!() to just to the lazy multiplication α * A * B + β C for T==BigFloat without the fancy BLAS-like machinery. Any alternative ideas I am missing before I do that?

( btw which entry does muladd!() overwrite? is it A? the docs are a bit unclear at https://github.com/JuliaLinearAlgebra/ArrayLayouts.jl/blob/master/src/muladd.jl )

dlfivefifty · 2023-07-27T15:37:52Z

🤷‍♂️ I can't find the code where I had to work around this. But definitely don't specialise to T == BigFloat, just use isprimitivetype

src/choleskyQR.jl

TSGut · 2023-07-27T15:52:14Z

I wasn't aware of isprimitivetype but it's exactly what I was looking for, nice.

TSGut · 2023-07-27T15:55:23Z

We technically don't need Cholesky to work for SemiclassicalOPs but it would feel silly to leave that part of the issue open. I think we just need muladd!() to have a method for bigfloats, even if it isn't as elegant.

dlfivefifty · 2023-07-30T15:19:38Z

src/choleskyQR.jl

-        # we fill 1 entry on the first run
-    dv = zeros(T,2)
+    # we fill 1 entry on the first run and circumvent BigFloat issue
+    dv = zeros(T,2) 


Try not to add extra spaces to ends of lines

TSGut added 4 commits July 26, 2023 19:12

fix QR-X BigFloat issue

954afe3

add test and broken test for Cholesky

41c3095

fix test

ab8e92f

fix isa vs ==

3f8a8f9

dlfivefifty requested changes Jul 27, 2023

View reviewed changes

src/choleskyQR.jl Outdated Show resolved Hide resolved

TSGut added 2 commits July 27, 2023 16:51

use zeros for length 2

df748d6

up version a bit

1d2f29f

use isprimitivetype

729870c

dlfivefifty reviewed Jul 30, 2023

View reviewed changes

dlfivefifty merged commit f193783 into JuliaApproximation:main Jul 30, 2023

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

BigFloat compatibility for qr_jacobimatrix #146

BigFloat compatibility for qr_jacobimatrix #146

Uh oh!

TSGut commented Jul 26, 2023

Uh oh!

codecov bot commented Jul 26, 2023 •

edited

Loading

Uh oh!

TSGut commented Jul 26, 2023

Uh oh!

dlfivefifty commented Jul 27, 2023

Uh oh!

Uh oh!

TSGut commented Jul 27, 2023

Uh oh!

TSGut commented Jul 27, 2023 •

edited

Loading

Uh oh!

dlfivefifty Jul 30, 2023

Uh oh!

Uh oh!

BigFloat compatibility for qr_jacobimatrix #146

BigFloat compatibility for qr_jacobimatrix #146

Uh oh!

Conversation

TSGut commented Jul 26, 2023

Uh oh!

codecov bot commented Jul 26, 2023 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

TSGut commented Jul 26, 2023

Uh oh!

dlfivefifty commented Jul 27, 2023

Uh oh!

Uh oh!

TSGut commented Jul 27, 2023

Uh oh!

TSGut commented Jul 27, 2023 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

dlfivefifty Jul 30, 2023

Choose a reason for hiding this comment

Uh oh!

Uh oh!

codecov bot commented Jul 26, 2023 •

edited

Loading

TSGut commented Jul 27, 2023 •

edited

Loading