Implement BidiagonalConjugationData #186

DanielVandH · 2024-06-30T21:44:16Z

Closes #185

This PR implements BidiagonalConjugationData which represents products of the form $A = U^{-1}XV$ where $A$ is upper bidiagonal, $U$ is upper Hessenberg, and $X$ and $V$ are banded (doesn't particularly matter what their bandwidth is, but them being banded makes sure that C = X * V in the constructor is efficient).

The equations used in _compute_column! come from noting that, letting $C = XV$,

$$c_{11} = u_{11}a_{11}$$

and then, for $i > 1$,

$$\begin{align*} c_{ii} &= u_{i,i-1}a_{i-1,i} + u_{ii}a_{ii}, \\\ c_{i-1,i} &= u_{i-1,i-1}a_{i-1,i} + u_{i-1,i}a_{ii}, \end{align*}$$

and these two equations get inverted to give $(a_{i-1, i}, a_{ii})$ for each $i>1$.

Testing uses InfRandVectors from InfiniteArrays.jl which is currently only available on InfiniteArrays#master. I probably want to find a fix for JuliaArrays/InfiniteArrays.jl#182 before the new InfiniteArrays release gets tagged so that this PR can also be merged.

dlfivefifty · 2024-07-01T20:47:13Z

src/banded/bidiagonalconjugationdata.jl

+- `V` is banded.
+
+None of these properties are checked internally. It is the user's responsibility
+to ensure these properties hold and that the product `inv(U)XV` is indeed bidiagonal.


Does our algorithm change if we just think of this as a bidiagonal projection of a possibly dense matrix? If not we don't have to include this caveat.

It doesn't change, so yeah I can remove that note.

If this computes bidiagonal projections in the general case, are you saying if I just wrap this in Symmetric(...) then it would work for the orthonormal symmetric tridiagonal Jacobi matrices given connection coefficients in U and V?

@TSGut I'm not exactly familiar with what you mean. Can you give an example?

src/banded/bidiagonalconjugationdata.jl

DanielVandH · 2024-07-02T22:01:33Z

I've redesigned the implementation based on your feedback. It needs JuliaArrays/LazyArrays.jl#326 and, for the tests, it needs InfiniteRandomArrays.jl's registration to complete

DanielVandH · 2024-07-04T14:30:35Z

Tests are passing now, so ready for another review

codecov · 2024-07-04T14:30:38Z

Codecov Report

Attention: Patch coverage is 92.50000% with 6 lines in your changes missing coverage. Please review.

Project coverage is 83.40%. Comparing base (3e81f50) to head (ebc1d2f).

Files	Patch %	Lines
src/banded/bidiagonalconjugation.jl	92.50%	6 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #186      +/-   ##
==========================================
+ Coverage   82.99%   83.40%   +0.41%     
==========================================
  Files          12       13       +1     
  Lines        1511     1591      +80     
==========================================
+ Hits         1254     1327      +73     
- Misses        257      264       +7

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

dlfivefifty · 2024-07-04T14:55:40Z

I still think that the interface of using Bidiagonal{T,<:BidiagonalConjugationBands} where each diagonal points to the same BidiagonalConjugationData is a more versatile and efficient design. For example, in means the bands act like cached vectors so can be accessed very fast, whereas your wrapper interface will be slow.

Perhaps we discuss tomorrow?

DanielVandH · 2024-07-04T14:59:34Z

Yes let's discuss tomorrow. I was thinking to have separate bands since my original attempt to have them be the same failed somewhere, but probably because I was confusing myself.

DanielVandH · 2024-07-06T12:23:15Z

I've redesigned it but it's still a bit broken in some areas... what methods am I missing to allow operations like the ones below to work? I can't see why it's trying to call BidiagonalConjugationData with a single index in some areas.

julia> V = InfRandTridiagonal(); A = InfRandBidiagonal('U'); X = brand(∞, 0, 2); U = X * V * ApplyArray(inv, A);

julia> B = BidiagonalConjugation(U, X, V, 'U');

julia> B * I
ℵ₀×ℵ₀ Bidiagonal{Float64, BroadcastVector{Float64, typeof(*), Tuple{InfiniteLinearAlgebra.BidiagonalConjugationBand{Float64}, Bool}}} with indices OneToInf()×OneToInf():
Error showing value of type Bidiagonal{Float64, BroadcastVector{Float64, typeof(*), Tuple{InfiniteLinearAlgebra.BidiagonalConjugationBand{Float64}, Bool}}}:
ERROR: MethodError: no method matching getindex(::InfiniteLinearAlgebra.BidiagonalConjugationData{Float64}, ::CartesianIndex{1})
Stacktrace:
  [1] getindex(A::InfiniteLinearAlgebra.BidiagonalConjugationBand{Float64}, I::CartesianIndex{1})
    @ LazyArrays C:\Users\djv23\.julia\packages\LazyArrays\VPs0M\src\cache.jl:207
  [2] _broadcast_getindex
    @ .\broadcast.jl:662 [inlined]
  [3] _getindex
    @ .\broadcast.jl:705 [inlined]
  [4] _broadcast_getindex
    @ .\broadcast.jl:681 [inlined]
  [5] getindex
    @ .\broadcast.jl:636 [inlined]
  [6] getindex
    @ C:\Users\djv23\.julia\packages\LazyArrays\VPs0M\src\lazybroadcasting.jl:93 [inlined]
  [7] isassigned(A::BroadcastVector{Float64, typeof(*), Tuple{InfiniteLinearAlgebra.BidiagonalConjugationBand{Float64}, Bool}}, i::Int64)
    @ Base .\multidimensional.jl:1578
  [8] isassigned(A::Bidiagonal{Float64, BroadcastVector{Float64, typeof(*), Tuple{…}}}, i::Int64, j::Int64)
    @ LinearAlgebra C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\bidiag.jl:136
  [9] alignment(io::IOContext{…}, X::AbstractVecOrMat, rows::Vector{…}, cols::Vector{…}, cols_if_complete::Int64, cols_otherwise::Int64, sep::Int64, ncols::Infinities.InfiniteCardinal{…})
    @ Base .\arrayshow.jl:68
 [10] _print_matrix(io::IOContext{…}, X::AbstractVecOrMat, pre::String, sep::String, post::String, hdots::String, vdots::String, ddots::String, hmod::Int64, vmod::Int64, rowsA::InfiniteArrays.InfUnitRange{…}, colsA::InfiniteArrays.InfUnitRange{…})
    @ Base .\arrayshow.jl:207
 [11] print_matrix(io::IOContext{…}, X::Bidiagonal{…}, pre::String, sep::String, post::String, hdots::String, vdots::String, ddots::String, hmod::Int64, vmod::Int64)
    @ Base .\arrayshow.jl:171
 [12] print_matrix
    @ .\arrayshow.jl:171 [inlined]
 [13] print_array
    @ .\arrayshow.jl:358 [inlined]
 [14] show(io::IOContext{Base.TTY}, ::MIME{Symbol("text/plain")}, X::Bidiagonal{Float64, BroadcastVector{Float64, typeof(*), Tuple{…}}})
    @ Base .\arrayshow.jl:399
 [15] (::REPL.var"#55#56"{REPL.REPLDisplay{REPL.LineEditREPL}, MIME{Symbol("text/plain")}, Base.RefValue{Any}})(io::Any)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:273
 [16] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:569
 [17] display(d::REPL.REPLDisplay, mime::MIME{Symbol("text/plain")}, x::Any)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:259
 [18] display(d::REPL.REPLDisplay, x::Any)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:278
 [19] display(x::Any)
    @ Base.Multimedia .\multimedia.jl:340
 [20] #invokelatest#2
    @ .\essentials.jl:892 [inlined]
 [21] invokelatest
    @ .\essentials.jl:889 [inlined]
 [22] print_response(errio::IO, response::Any, show_value::Bool, have_color::Bool, specialdisplay::Union{Nothing, AbstractDisplay})
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:315
 [23] (::REPL.var"#57#58"{REPL.LineEditREPL, Pair{Any, Bool}, Bool, Bool})(io::Any)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:284
 [24] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:569
 [25] print_response(repl::REPL.AbstractREPL, response::Any, show_value::Bool, have_color::Bool)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:282
 [26] (::REPL.var"#do_respond#80"{…})(s::REPL.LineEdit.MIState, buf::Any, ok::Bool)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:911
 [27] (::VSCodeServer.var"#103#106"{REPL.var"#do_respond#80"{…}})(mi::REPL.LineEdit.MIState, buf::IOBuffer, ok::Bool)
    @ VSCodeServer c:\Users\djv23\.vscode\extensions\julialang.language-julia-1.83.2\scripts\packages\VSCodeServer\src\repl.jl:122
 [28] #invokelatest#2
    @ .\essentials.jl:892 [inlined]
 [29] invokelatest
    @ .\essentials.jl:889 [inlined]
 [30] run_interface(terminal::REPL.Terminals.TextTerminal, m::REPL.LineEdit.ModalInterface, s::REPL.LineEdit.MIState)
    @ REPL.LineEdit C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\LineEdit.jl:2656
 [31] run_frontend(repl::REPL.LineEditREPL, backend::REPL.REPLBackendRef)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:1312
 [32] (::REPL.var"#62#68"{REPL.LineEditREPL, REPL.REPLBackendRef})()
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:386
Some type information was truncated. Use `show(err)` to see complete types.

julia> B - A
ℵ₀×ℵ₀ Bidiagonal{Float64, BroadcastVector{Float64, typeof(-), Tuple{InfiniteLinearAlgebra.BidiagonalConjugationBand{Float64}, InfiniteRandomArrays.InfRandVector{Float64, Type{Float64}, Xoshiro}}}} with indices OneToInf()×OneToInf():
Error showing value of type Bidiagonal{Float64, BroadcastVector{Float64, typeof(-), Tuple{InfiniteLinearAlgebra.BidiagonalConjugationBand{Float64}, InfiniteRandomArrays.InfRandVector{Float64, Type{Float64}, Xoshiro}}}}:
ERROR: MethodError: no method matching getindex(::InfiniteLinearAlgebra.BidiagonalConjugationData{Float64}, ::CartesianIndex{1})
Stacktrace:
  [1] getindex(A::InfiniteLinearAlgebra.BidiagonalConjugationBand{Float64}, I::CartesianIndex{1})
    @ LazyArrays C:\Users\djv23\.julia\packages\LazyArrays\VPs0M\src\cache.jl:207
  [2] _broadcast_getindex
    @ .\broadcast.jl:662 [inlined]
  [3] _getindex
    @ .\broadcast.jl:705 [inlined]
  [4] _broadcast_getindex
    @ .\broadcast.jl:681 [inlined]
  [5] getindex
    @ .\broadcast.jl:636 [inlined]
  [6] getindex
    @ C:\Users\djv23\.julia\packages\LazyArrays\VPs0M\src\lazybroadcasting.jl:93 [inlined]
  [7] isassigned(A::BroadcastVector{Float64, typeof(-), Tuple{…}}, i::Int64)
    @ Base .\multidimensional.jl:1578
  [8] isassigned(A::Bidiagonal{Float64, BroadcastVector{Float64, typeof(-), Tuple{…}}}, i::Int64, j::Int64)
    @ LinearAlgebra C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\bidiag.jl:136
  [9] alignment(io::IOContext{…}, X::AbstractVecOrMat, rows::Vector{…}, cols::Vector{…}, cols_if_complete::Int64, cols_otherwise::Int64, sep::Int64, ncols::Infinities.InfiniteCardinal{…})
    @ Base .\arrayshow.jl:68
 [10] _print_matrix(io::IOContext{…}, X::AbstractVecOrMat, pre::String, sep::String, post::String, hdots::String, vdots::String, ddots::String, hmod::Int64, vmod::Int64, rowsA::InfiniteArrays.InfUnitRange{…}, colsA::InfiniteArrays.InfUnitRange{…})
    @ Base .\arrayshow.jl:207
 [11] print_matrix(io::IOContext{…}, X::Bidiagonal{…}, pre::String, sep::String, post::String, hdots::String, vdots::String, ddots::String, hmod::Int64, vmod::Int64)
    @ Base .\arrayshow.jl:171
 [12] print_matrix
    @ .\arrayshow.jl:171 [inlined]
 [13] print_array
    @ .\arrayshow.jl:358 [inlined]
 [14] show(io::IOContext{Base.TTY}, ::MIME{Symbol("text/plain")}, X::Bidiagonal{Float64, BroadcastVector{Float64, typeof(-), Tuple{…}}})
    @ Base .\arrayshow.jl:399
 [15] (::REPL.var"#55#56"{REPL.REPLDisplay{REPL.LineEditREPL}, MIME{Symbol("text/plain")}, Base.RefValue{Any}})(io::Any)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:273
 [16] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:569
 [17] display(d::REPL.REPLDisplay, mime::MIME{Symbol("text/plain")}, x::Any)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:259
 [18] display(d::REPL.REPLDisplay, x::Any)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:278
 [19] display(x::Any)
    @ Base.Multimedia .\multimedia.jl:340
 [20] #invokelatest#2
    @ .\essentials.jl:892 [inlined]
 [21] invokelatest
    @ .\essentials.jl:889 [inlined]
 [22] print_response(errio::IO, response::Any, show_value::Bool, have_color::Bool, specialdisplay::Union{Nothing, AbstractDisplay})
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:315
 [23] (::REPL.var"#57#58"{REPL.LineEditREPL, Pair{Any, Bool}, Bool, Bool})(io::Any)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:284
 [24] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:569
 [25] print_response(repl::REPL.AbstractREPL, response::Any, show_value::Bool, have_color::Bool)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:282
 [26] (::REPL.var"#do_respond#80"{…})(s::REPL.LineEdit.MIState, buf::Any, ok::Bool)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:911
 [27] (::VSCodeServer.var"#103#106"{REPL.var"#do_respond#80"{…}})(mi::REPL.LineEdit.MIState, buf::IOBuffer, ok::Bool)
    @ VSCodeServer c:\Users\djv23\.vscode\extensions\julialang.language-julia-1.83.2\scripts\packages\VSCodeServer\src\repl.jl:122
 [28] #invokelatest#2
    @ .\essentials.jl:892 [inlined]
 [29] invokelatest
    @ .\essentials.jl:889 [inlined]
 [30] run_interface(terminal::REPL.Terminals.TextTerminal, m::REPL.LineEdit.ModalInterface, s::REPL.LineEdit.MIState)
    @ REPL.LineEdit C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\LineEdit.jl:2656
 [31] run_frontend(repl::REPL.LineEditREPL, backend::REPL.REPLBackendRef)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:1312
 [32] (::REPL.var"#62#68"{REPL.LineEditREPL, REPL.REPLBackendRef})()
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:386
Some type information was truncated. Use `show(err)` to see complete types.

julia> B * B
ERROR: InterruptException:
Stacktrace:
 [1] cache_filldata!(::LazyArrays.CachedArray{Float64, 2, Matrix{…}, Zeros{…}}, ::UnitRange{Int64}, ::Base.OneTo{Int64})
   @ LazyArrays C:\Users\djv23\.julia\packages\LazyArrays\VPs0M\src\cache.jl:226
 [2] resizedata!(::DenseColumnMajor, ::ZerosLayout, ::LazyArrays.CachedArray{…}, ::Int64, ::Int64)
   @ LazyArrays C:\Users\djv23\.julia\packages\LazyArrays\VPs0M\src\cache.jl:267
 [3] resizedata!
   @ C:\Users\djv23\.julia\packages\LazyArrays\VPs0M\src\cache.jl:218 [inlined]
 [4] getindex
   @ C:\Users\djv23\.julia\packages\LazyArrays\VPs0M\src\cache.jl:70 [inlined]
 [5] _mul!(C::LazyArrays.CachedArray{…}, A::Bidiagonal{…}, B::Bidiagonal{…}, _add::LinearAlgebra.MulAddMul{…})
   @ LinearAlgebra C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\bidiag.jl:500
 [6] mul!
   @ C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\bidiag.jl:429 [inlined]
 [7] mul!
   @ C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\matmul.jl:237 [inlined]
 [8] *(A::Bidiagonal{…}, B::Bidiagonal{…})
   @ LinearAlgebra C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\matmul.jl:106
 [9] top-level scope
   @ REPL[43]:1
Some type information was truncated. Use `show(err)` to see complete types.

I thought I can set BroadcastStyle(::Type{<:BidiagonalConjugation})= LazyArrayStyle{2}() but that didn't work. I think the last example I gave requires JuliaLinearAlgebra/ArrayLayouts.jl#241 to be merged, but I haven't checked.

dlfivefifty · 2024-07-06T12:29:35Z

I'll have a look

dlfivefifty · 2024-07-06T12:36:17Z

Ah the problem is that AbstractCachedArray assumes the field data stores the entries. So I think we just don't make it <: AbstractCachedVector

dlfivefifty · 2024-07-06T12:41:43Z

I fixed the tests but julia> n = 100_000; @time ArrayLayouts.layout_getindex(B,1:n,1:n) is insanely slow so we probably need to make sure its resizing

dlfivefifty · 2024-07-06T12:42:52Z

Actually there's something wrong in resizedata!: julia> @time LazyArrays.resizedata!(B.dv.data, 100_000); freezes

dlfivefifty · 2024-07-06T12:44:13Z

Oh wait your U is dense

DanielVandH · 2024-07-06T13:17:44Z

I can maybe add a test where all the matrices are banded and check the timing? It's dense because it was the simplest way to get some test matrices where I know what A is exactly.

Are there any other changes you need to commit before I work some more on this?

DanielVandH · 2024-07-06T15:24:39Z

Some benchmarks on the latest commit for sparse U, X, V

julia> using InfiniteLinearAlgebra, ArrayLayouts, BandedMatrices, InfiniteRandomArrays, BenchmarkTools

julia> using InfiniteLinearAlgebra: BidiagonalConjugation

julia> using LazyArrays

julia> @benchmark ArrayLayouts.layout_getindex(A, 1:n, 1:n) setup = begin
           U = brand(∞, ∞, 0, 1)
           X = brand(∞, ∞, -1, 2)
           V = brand(∞, ∞, 1, 0)
           A = BidiagonalConjugation(U, X, V, 'U')
           n = 100_000
       end
BenchmarkTools.Trial: 50 samples with 1 evaluation.
 Range (min … max):   67.689 ms … 170.047 ms  ┊ GC (min … max):  8.29% … 41.76%
 Time  (median):      96.980 ms               ┊ GC (median):     7.30%
 Time  (mean ± σ):   100.439 ms ±  21.801 ms  ┊ GC (mean ± σ):  10.92% ±  8.81%

  ▄     ▁▁▁▁    ▄▄█▁   ▄ ▁▁▄  ▁▁                              ▁
  █▆▁▁▆▆████▁▁▁▁████▁▆▆█▁███▁▆██▁▆▆▆▆▁▁▆▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  67.7 ms          Histogram: frequency by time          170 ms <

 Memory estimate: 77.12 MiB, allocs estimate: 786490.

julia> @benchmark LazyArrays.resizedata!(A.dv.data, n) setup = begin
           U = brand(∞, ∞, 0, 1)
           X = brand(∞, ∞, -1, 2)
           V = brand(∞, ∞, 1, 0)
           A = BidiagonalConjugation(U, X, V, 'U')
           n = 100_000
       end
BenchmarkTools.Trial: 51 samples with 1 evaluation.
 Range (min … max):  82.229 ms … 183.218 ms  ┊ GC (min … max): 6.99% … 37.55%
 Time  (median):     93.903 ms               ┊ GC (median):    7.64%
 Time  (mean ± σ):   98.143 ms ±  16.153 ms  ┊ GC (mean ± σ):  9.14% ±  5.85%

             █  ▁     
  ▄▄▁▆▄▁▄▆▆▇▁█▄▇█▆▇▆▇▄▄▁▁▄▁▁▁▁▄▄▁▁▁▁▁▄▁▁▁▁▁▁▄▁▁▁▄▁▄▁▁▁▁▄▁▁▁▁▁▄ ▁
  82.2 ms         Histogram: frequency by time          132 ms <

 Memory estimate: 73.66 MiB, allocs estimate: 800036.

Quite a lot of allocations... would've thought the function barriers would be handling the non-concretely typed U and C better than that. Unless the resizing of the vectors is not lumped together as a single allocation in the estimate maybe

Nevermind, it's probably because taking slices is resizing for each index rather than doing the maximum

DanielVandH · 2024-07-06T15:45:23Z

It's actually not that bad.

julia> ArrayLayouts.layout_getindex(A, 1:100_000, 1:100_000)
n = 11
n = 23
n = 47
n = 95
n = 191
n = 383
n = 767
n = 1535
n = 3071
n = 6143
n = 12287
n = 24575
n = 49151
n = 98303
100000×100000 BandedMatrix{Float64} with bandwidths (0, 1):
 0.291719  0.157616   ⋅          ⋅          ⋅         ⋅          ⋅       …   ⋅    ⋅         ⋅         ⋅          ⋅         ⋅         ⋅ 
  ⋅        1.11719   0.737074    ⋅          ⋅         ⋅          ⋅           ⋅    ⋅         ⋅         ⋅          ⋅         ⋅         ⋅
  ⋅         ⋅        0.309456  -0.338909    ⋅         ⋅          ⋅           ⋅    ⋅         ⋅         ⋅          ⋅         ⋅         ⋅
  ⋅         ⋅         ⋅         3.58005   -1.29835    ⋅          ⋅           ⋅    ⋅         ⋅         ⋅          ⋅         ⋅         ⋅
  ⋅         ⋅         ⋅          ⋅         0.648818  0.640061    ⋅           ⋅    ⋅         ⋅         ⋅          ⋅         ⋅         ⋅
  ⋅         ⋅         ⋅          ⋅          ⋅        0.0280967  5.15181  …   ⋅    ⋅         ⋅         ⋅          ⋅         ⋅         ⋅
 ⋮                                                   ⋮                   ⋱                 ⋮
  ⋅         ⋅         ⋅          ⋅          ⋅         ⋅          ⋅           ⋅   0.671618  0.474535   ⋅          ⋅         ⋅         ⋅
  ⋅         ⋅         ⋅          ⋅          ⋅         ⋅          ⋅       …   ⋅    ⋅        0.183867  4.7101      ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅          ⋅          ⋅         ⋅          ⋅           ⋅    ⋅         ⋅        0.196432  -1.06407    ⋅         ⋅
  ⋅         ⋅         ⋅          ⋅          ⋅         ⋅          ⋅           ⋅    ⋅         ⋅         ⋅         0.766571  1.58069    ⋅
  ⋅         ⋅         ⋅          ⋅          ⋅         ⋅          ⋅           ⋅    ⋅         ⋅         ⋅          ⋅        0.112282  2.25429       
  ⋅         ⋅         ⋅          ⋅          ⋅         ⋅          ⋅           ⋅    ⋅         ⋅         ⋅          ⋅         ⋅        0.329934

It does resize a bit often but not as much as I thought. I tried getting around it and doing one big resize but couldn't get the details working while making everything else efficient. Maybe this is fine to go? See what you think.

dlfivefifty · 2024-07-07T10:37:17Z

src/banded/bidiagonalconjugation.jl

+    ds = data.datasize
+    n ≤ ds && return data
+    dv, ev = data.dv, data.ev
+    resize!(dv, 2n + 1)


I believe doing this may cause it to be an O(n^2) resize: if it has to copy memory everytime we call resize! then something like

for k=1:n resizedata!(data, k) end

will be $O(n^2)$.

In CachedArray this is avoided by doubling the size whenever we call resize!, see

https://github.com/JuliaArrays/LazyArrays.jl/blob/f2ed39925bad206351cd4fd7c0978310aacd52cc/src/cache.jl#L238

This is why we need datasize instead of just relying on the size of the array.

Maybe I'm not understanding _vec_resizedata! correctly: Aren't I already doubling the size of the array here? That's what I was trying to do with the 2n instead of just n

The main difference is _vec_resizedata! only calls resize! when n is bigger than the length. Your version is always calling resize! so will be $O(n^2)$

Doesn't n ≤ ds && return data make sure that no resizing actually happens when n exceeds the length?

But ds is not the same as length(dv)

I'm a bit confused on how I can change this to address what you mean. Maybe it's obvious but I can't see it sorry. From what I can see, it's resizing as I would expect.

Here's an example. I put a @show n in resizedata! after my check n <= ds to see exactly when a resize is actually being performed. When resizing dv:

V, A, X = InfRandTridiagonal(), InfRandBidiagonal('U'), brand(∞, 0, 2); U = X * V * ApplyArray(inv, A); B = BidiagonalConjugation(U, X, V, 'U'); dv, ev = B.dv, B.ev; for n in 1:24 println("Loop iter dv: $n") LazyArrays.resizedata!(dv, n) end

Loop iter dv: 1 n = 1 Loop iter dv: 2 Loop iter dv: 3 Loop iter dv: 4 n = 4 Loop iter dv: 5 Loop iter dv: 6 Loop iter dv: 7 Loop iter dv: 8 Loop iter dv: 9 Loop iter dv: 10 n = 10 Loop iter dv: 11 Loop iter dv: 12 Loop iter dv: 13 Loop iter dv: 14 Loop iter dv: 15 Loop iter dv: 16 Loop iter dv: 17 Loop iter dv: 18 Loop iter dv: 19 Loop iter dv: 20 Loop iter dv: 21 Loop iter dv: 22 n = 22 Loop iter dv: 23 Loop iter dv: 24

When resizing ev:

B = BidiagonalConjugation(U, X, V, 'U'); dv, ev = B.dv, B.ev; for n in 1:24 println("Loop iter: ev $n") LazyArrays.resizedata!(ev, n) end

Loop iter: ev 1 n = 1 Loop iter: ev 2 Loop iter: ev 3 Loop iter: ev 4 n = 4 Loop iter: ev 5 Loop iter: ev 6 Loop iter: ev 7 Loop iter: ev 8 Loop iter: ev 9 Loop iter: ev 10 n = 10 Loop iter: ev 11 Loop iter: ev 12 Loop iter: ev 13 Loop iter: ev 14 Loop iter: ev 15 Loop iter: ev 16 Loop iter: ev 17 Loop iter: ev 18 Loop iter: ev 19 Loop iter: ev 20 Loop iter: ev 21 Loop iter: ev 22 n = 22 Loop iter: ev 23 Loop iter: ev 24

What that is showing is that it is doubling the populated data everytime it resizes. This is also problematic as it means it will need double the computation necessary.

Ah. Yes good point. I will try and rewrite it to do it efficiently.

What do you think of the new change?

src/banded/bidiagonalconjugation.jl

Implement BidiagonalConjugationData

e643dfa

dlfivefifty requested changes Jul 1, 2024

View reviewed changes

DanielVandH added 2 commits July 2, 2024 22:58

Redesign

3e455e9

Fix import and make infiniterandomarrays an extra

e79939a

DanielVandH added 2 commits July 2, 2024 23:02

typo

e7c37bd

Tests

e2cb231

Redesign

563e0b6

Function barriers

550a03a

BidiagonalConjugationBand should not be <: AbstractCachedVector

c45472f

oops

f9ded91

Changes

833a6a9

dlfivefifty requested changes Jul 7, 2024

View reviewed changes

DanielVandH added 3 commits July 7, 2024 14:17

Make copy a no-op

e9f33d9

Cleanup

a74a47e

Fix resize

ebc1d2f

dlfivefifty approved these changes Jul 12, 2024

View reviewed changes

DanielVandH mentioned this pull request Jul 14, 2024

Implement b=-1 JuliaApproximation/SemiclassicalOrthogonalPolynomials.jl#111

Merged

dlfivefifty merged commit 6e56195 into JuliaLinearAlgebra:master Jul 14, 2024
6 checks passed

DanielVandH deleted the bidiag branch July 14, 2024 09:57

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Implement BidiagonalConjugationData #186

Implement BidiagonalConjugationData #186

DanielVandH commented Jun 30, 2024 •

edited

Loading

dlfivefifty Jul 1, 2024

DanielVandH Jul 1, 2024

TSGut Jul 1, 2024

DanielVandH Jul 4, 2024

DanielVandH commented Jul 2, 2024

DanielVandH commented Jul 4, 2024

codecov bot commented Jul 4, 2024 •

edited

Loading

dlfivefifty commented Jul 4, 2024

DanielVandH commented Jul 4, 2024

DanielVandH commented Jul 6, 2024

dlfivefifty commented Jul 6, 2024

dlfivefifty commented Jul 6, 2024

dlfivefifty commented Jul 6, 2024

dlfivefifty commented Jul 6, 2024

dlfivefifty commented Jul 6, 2024

DanielVandH commented Jul 6, 2024

DanielVandH commented Jul 6, 2024 •

edited

Loading

DanielVandH commented Jul 6, 2024

dlfivefifty Jul 7, 2024

DanielVandH Jul 7, 2024

dlfivefifty Jul 7, 2024

DanielVandH Jul 7, 2024

dlfivefifty Jul 7, 2024

DanielVandH Jul 8, 2024

dlfivefifty Jul 10, 2024

DanielVandH Jul 10, 2024

DanielVandH Jul 12, 2024

Implement BidiagonalConjugationData #186

Implement BidiagonalConjugationData #186

Conversation

DanielVandH commented Jun 30, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

DanielVandH commented Jul 2, 2024

DanielVandH commented Jul 4, 2024

codecov bot commented Jul 4, 2024 • edited Loading

Codecov Report

dlfivefifty commented Jul 4, 2024

DanielVandH commented Jul 4, 2024

DanielVandH commented Jul 6, 2024

dlfivefifty commented Jul 6, 2024

dlfivefifty commented Jul 6, 2024

dlfivefifty commented Jul 6, 2024

dlfivefifty commented Jul 6, 2024

dlfivefifty commented Jul 6, 2024

DanielVandH commented Jul 6, 2024

DanielVandH commented Jul 6, 2024 • edited Loading

DanielVandH commented Jul 6, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

DanielVandH commented Jun 30, 2024 •

edited

Loading

codecov bot commented Jul 4, 2024 •

edited

Loading

DanielVandH commented Jul 6, 2024 •

edited

Loading