Implement BidiagonalConjugationData

Project.toml

@@ -1,6 +1,6 @@
 name = "InfiniteLinearAlgebra"
 uuid = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
-version = "0.8.3"
+version = "0.8.4"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -23,8 +23,9 @@ BandedMatrices = "1.7.2"
 BlockArrays = "1.0"
 BlockBandedMatrices = "0.13"
 FillArrays = "1.0"
-Infinities = "0.1"
 InfiniteArrays = "0.14"
+InfiniteRandomArrays = "0.2"
+Infinities = "0.1"
 LazyArrays = "2.0"
 LazyBandedMatrices = "0.10"
 LinearAlgebra = "1"
@@ -39,9 +40,10 @@ julia = "1.10"
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+InfiniteRandomArrays = "2bc77966-89c7-476d-a40f-269028fac4a9"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "Test", "Random", "SpecialFunctions", "StaticArrays"]
+test = ["Aqua", "Test", "Random", "InfiniteRandomArrays", "SpecialFunctions", "StaticArrays"]

src/InfiniteLinearAlgebra.jl

@@ -21,7 +21,7 @@ import ArrayLayouts: AbstractBandedLayout, AbstractQLayout, AdjQRPackedQLayout,
 
 import BandedMatrices: BandedColumns, BandedMatrix, BandedMatrix, _BandedMatrix, AbstractBandedMatrix,
                        _BandedMatrix, _BandedMatrix, _banded_qr, _banded_qr!, _default_banded_broadcast, banded_chol!,
-                       banded_similar, bandedcolumns, bandeddata, bandwidths, bandwidths
+                       banded_similar, bandedcolumns, bandeddata, bandwidths
 
 import BlockArrays: AbstractBlockLayout, BlockLayout, BlockSlice, BlockSlice1, BlockedOneTo,
                     blockcolsupport, sizes_from_blocks, OneToCumsum, AbstractBlockedUnitRange
@@ -38,9 +38,9 @@ import Infinities: InfiniteCardinal, Infinity
 
 import LazyArrays: AbstractCachedMatrix, AbstractCachedVector, AbstractLazyLayout, ApplyArray, ApplyLayout, ApplyMatrix,
                    CachedArray, CachedLayout, CachedMatrix, CachedVector, LazyArrayStyle, LazyLayout,
-                   LazyLayouts, LazyMatrix, AbstractPaddedLayout, PaddedColumns, _broadcast_sub_arguments,
+                   LazyLayouts, LazyMatrix, LazyVector, AbstractPaddedLayout, PaddedColumns, _broadcast_sub_arguments,
                    applybroadcaststyle, applylayout, arguments, cacheddata, paddeddata, resizedata!, simplifiable,
-                   simplify, islazy, islazy_layout
+                   simplify, islazy, islazy_layout, cache_getindex
 
 import LazyBandedMatrices: AbstractLazyBandedBlockBandedLayout, AbstractLazyBandedLayout, ApplyBandedLayout, BlockVec,
                            BroadcastBandedLayout, KronTravBandedBlockBandedLayout, LazyBandedLayout,
@@ -143,5 +143,6 @@ include("infql.jl")
 include("infqr.jl")
 include("inful.jl")
 include("infcholesky.jl")
+include("banded/bidiagonalconjugation.jl")
 
 end # module

src/banded/bidiagonalconjugation.jl

@@ -0,0 +1,138 @@
+@inline function _to_uplo(char::Symbol)
+    if char == :U
+        'U'
+    elseif char == :L
+        'L'
+    else
+        _throw_uplo()
+    end
+end
+@inline function _to_uplo(char::Char)
+    if char ∈ ('L', 'U')
+        char
+    else
+        _throw_uplo()
+    end
+end
+@noinline _throw_uplo() = throw(ArgumentError("uplo argument must be either :U (upper) or :L (lower)"))
+
+mutable struct BidiagonalConjugationData{T}
+    const U::AbstractMatrix{T} # Typing these concretely prevents the use of Bidiagonal, unless we want LazyBandedMatrices.Bidiagonal
+    const C::AbstractMatrix{T} # Function barriers help to minimise the penalty from this when resizing anyway.
+    const dv::Vector{T}
+    const ev::Vector{T}
+    const uplo::Char
+    datasize::Int # Number of columns
+end
+function BidiagonalConjugationData(U, X, V, uplo::Char)
+    C = X * V
+    T = promote_type(typeof(inv(U[1, 1])), eltype(U), eltype(C)) # include inv so that we can't get Ints
+    dv, ev = T[], T[]
+    return BidiagonalConjugationData(U, C, dv, ev, uplo, 0)
+end
+
+function copy(data::BidiagonalConjugationData)
+    U, C, dv, ev, uplo, datasize = data.U, data.C, data.dv, data.ev, data.uplo, data.datasize
+    return BidiagonalConjugationData(copy(U), copy(C), copy(dv), copy(ev), uplo, datasize)
+end
+
+function _compute_column_up!(data::BidiagonalConjugationData, U, C, i)
+    dv, ev = data.dv, data.ev
+    if i == 1
+        dv[i] = C[1, 1] / U[1, 1]
+    else
+        uᵢ₋₁ᵢ₋₁, uᵢᵢ₋₁, uᵢ₋₁ᵢ, uᵢᵢ = U[i-1, i-1], U[i, i-1], U[i-1, i], U[i, i]
+        cᵢ₋₁ᵢ, cᵢᵢ = C[i-1, i], C[i, i]
+        Uᵢ⁻¹ = inv(uᵢ₋₁ᵢ₋₁ * uᵢᵢ - uᵢ₋₁ᵢ * uᵢᵢ₋₁)
+        dv[i] = Uᵢ⁻¹ * (uᵢ₋₁ᵢ₋₁ * cᵢᵢ - uᵢᵢ₋₁ * cᵢ₋₁ᵢ)
+        ev[i-1] = Uᵢ⁻¹ * (uᵢᵢ * cᵢ₋₁ᵢ - uᵢ₋₁ᵢ * cᵢᵢ)
+    end
+    return data
+end
+
+function _compute_column_lo!(data::BidiagonalConjugationData, U, C, i)
+    dv, ev = data.dv, data.ev
+    uᵢᵢ, uᵢ₊₁ᵢ, uᵢᵢ₊₁, uᵢ₊₁ᵢ₊₁ = U[i, i], U[i+1, i], U[i, i+1], U[i+1, i+1]
+    cᵢᵢ, cᵢ₊₁ᵢ = C[i, i], C[i+1, i]
+    Uᵢ⁻¹ = inv(uᵢᵢ * uᵢ₊₁ᵢ₊₁ - uᵢᵢ₊₁ * uᵢ₊₁ᵢ)
+    dv[i] = Uᵢ⁻¹ * (uᵢ₊₁ᵢ₊₁ * cᵢᵢ - uᵢᵢ₊₁ * cᵢ₊₁ᵢ)
+    ev[i] = Uᵢ⁻¹ * (uᵢᵢ * cᵢ₊₁ᵢ - uᵢ₊₁ᵢ * cᵢᵢ)
+    return data
+end
+
+function _compute_columns!(data::BidiagonalConjugationData, i)
+    U, C = data.U, data.C # Treat _compute_column_(up/lo) as function barriers and take these out early
+    return __compute_columns!(data, U, C, i)
+end
+function __compute_columns!(data::BidiagonalConjugationData, U, C, i)
+    ds = data.datasize
+    up = data.uplo == 'U'
+    for j in (ds+1):i
+        up ? _compute_column_up!(data, U, C, j) : _compute_column_lo!(data, U, C, j)
+    end
+    data.datasize = i
+    return data
+end
+
+function resizedata!(data::BidiagonalConjugationData, n)
+    n ≤ 0 && return data
+    v = data.datasize
+    n = max(v, n)
+    dv, ev = data.dv, data.ev
+    if n > length(ev) # Avoid O(n²) growing. Note min(length(dv), length(ev)) == length(ev)
+        resize!(dv, 2n + 1)
+        resize!(ev, 2n)
+    end
+    n > v && _compute_columns!(data, n)
+    return data
+end
+
+struct BidiagonalConjugationBand{T} <: LazyVector{T}
+    data::BidiagonalConjugationData{T}
+    diag::Bool # true => diagonal, false => offdiagonal 
+end
+@inline size(::BidiagonalConjugationBand) = (ℵ₀,)
+@inline resizedata!(A::BidiagonalConjugationBand, n) = resizedata!(A.data, n)
+
+function _bcb_getindex(band::BidiagonalConjugationBand, I)
+    resizedata!(band, maximum(I) + 1)
+    if band.diag
+        return band.data.dv[I]
+    else
+        return band.data.ev[I]
+    end
+end
+
+@inline getindex(band::BidiagonalConjugationBand, I::Integer) = _bcb_getindex(band, I)
+@inline getindex(band::BidiagonalConjugationBand, I::AbstractVector) = _bcb_getindex(band, I)
+
+copy(band::BidiagonalConjugationBand) = band
+
+const BidiagonalConjugation{T} = Bidiagonal{T,BidiagonalConjugationBand{T}}
+
+"""
+    BidiagonalConjugation(U, X, V, uplo)
+
+Efficiently compute the projection of the matrix product
+`inv(U)XV` onto a bidiagonal matrix. The `uplo` argument 
+specifies whether the projection is upper (`uplo = 'U'`)
+or lower (`uplo = 'L'`) bidiagonal. 
+
+The computation is returned as a `Bidiagonal` matrix whose 
+diagonal and off-diagonal vectors are computed lazily.
+"""
+function BidiagonalConjugation(U, X, V, uplo)
+    _uplo = _to_uplo(uplo)
+    data = BidiagonalConjugationData(U, X, V, _uplo)
+    return _BidiagonalConjugation(data, _uplo)
+end
+
+function _BidiagonalConjugation(data, uplo) # need uplo argument so that we can take transposes
+    dv = BidiagonalConjugationBand(data, true)
+    ev = BidiagonalConjugationBand(data, false)
+    return Bidiagonal(dv, ev, uplo)
+end
+
+copy(A::BidiagonalConjugation) = A # no-op
+
+LazyBandedMatrices.Bidiagonal(A::BidiagonalConjugation) = LazyBandedMatrices.Bidiagonal(A.dv, A.ev, A.uplo)

test/runtests.jl

@@ -3,15 +3,18 @@ using InfiniteLinearAlgebra, BlockBandedMatrices, BlockArrays, BandedMatrices, I
 import InfiniteLinearAlgebra: qltail, toeptail, tailiterate, tailiterate!, tail_de, ql_X!,
     InfToeplitz, PertToeplitz, TriToeplitz, InfBandedMatrix, InfBandCartesianIndices,
     rightasymptotics, QLHessenberg, ConstRows, PertConstRows, chop, chop!, pad,
-    BandedToeplitzLayout, PertToeplitzLayout, TridiagonalToeplitzLayout, BidiagonalToeplitzLayout
+    BandedToeplitzLayout, PertToeplitzLayout, TridiagonalToeplitzLayout, BidiagonalToeplitzLayout,
+    BidiagonalConjugation
 import Base: BroadcastStyle, oneto
 import BlockArrays: _BlockArray, blockcolsupport
 import BlockBandedMatrices: isblockbanded, _BlockBandedMatrix
 import MatrixFactorizations: QLPackedQ
 import BandedMatrices: bandeddata, _BandedMatrix, BandedStyle
-import LazyArrays: colsupport, MemoryLayout, ApplyLayout, LazyArrayStyle, arguments, paddeddata, PaddedColumns
+import LazyArrays: colsupport, MemoryLayout, ApplyLayout, LazyArrayStyle, arguments, paddeddata, PaddedColumns, LazyLayout
 import InfiniteArrays: OneToInf, oneto, RealInfinity
 import LazyBandedMatrices: BroadcastBandedBlockBandedLayout, BroadcastBandedLayout, LazyBandedLayout, BlockVec
+import InfiniteRandomArrays: InfRandTridiagonal, InfRandBidiagonal
+import ArrayLayouts: diagonaldata, supdiagonaldata, subdiagonaldata
 
 using Aqua
 @testset "Project quality" begin
@@ -476,3 +479,4 @@ include("test_inful.jl")
 include("test_infcholesky.jl")
 include("test_periodic.jl")
 include("test_infreversecholesky.jl")
+include("test_bidiagonalconjugation.jl")

test/test_bidiagonalconjugation.jl

@@ -0,0 +1,90 @@
+using InfiniteLinearAlgebra, InfiniteRandomArrays, BandedMatrices, LazyArrays, LazyBandedMatrices, InfiniteArrays, ArrayLayouts, Test
+using InfiniteLinearAlgebra: BidiagonalConjugation, OneToInf
+using ArrayLayouts: supdiagonaldata, subdiagonaldata, diagonaldata
+using LinearAlgebra
+using LazyArrays: LazyLayout
+
+@testset "BidiagonalConjugationData" begin
+    @test InfiniteLinearAlgebra._to_uplo('U') == 'U'
+    @test InfiniteLinearAlgebra._to_uplo('L') == 'L'
+    @test_throws ArgumentError InfiniteLinearAlgebra._to_uplo('a')
+    @test InfiniteLinearAlgebra._to_uplo(:U) == 'U'
+    @test InfiniteLinearAlgebra._to_uplo(:L) == 'L'
+    @test_throws ArgumentError InfiniteLinearAlgebra._to_uplo(:a)
+
+    for _ in 1:3
+        V1 = InfRandTridiagonal()
+        A1 = InfRandBidiagonal('U')
+        X1 = brand(∞, 0, 2)
+        U1 = X1 * V1 * ApplyArray(inv, A1)
+        B1 = BidiagonalConjugation(U1, X1, V1, 'U');
+
+        V2 = brand(∞, 0, 1)
+        A2 = LazyBandedMatrices.Bidiagonal(Fill(0.2, ∞), 2.0 ./ (1.0 .+ (1:∞)), 'L') # LinearAlgebra.Bidiagonal not playing nice for this case
+        X2 = InfRandBidiagonal('L')
+        U2 = X2 * V2 * ApplyArray(inv, A2) 
+        B2 = BidiagonalConjugation(U2, X2, V2, :L);
+
+        for (A, B, uplo) in ((A1, B1, 'U'), (A2, B2, 'L'))
+            @test B.dv.data === B.ev.data
+            @test MemoryLayout(B) isa BidiagonalLayout{LazyLayout,LazyLayout}
+            @test diagonaldata(B) === B.dv
+            if uplo == 'U'
+                @test supdiagonaldata(B) === B.ev
+                @test_throws ArgumentError subdiagonaldata(B)
+                @test bandwidths(B) == (0, 1)
+            else
+                @test subdiagonaldata(B) === B.ev
+                @test_throws ArgumentError supdiagonaldata(B)
+                @test bandwidths(B) == (1, 0)
+            end
+            @test size(B) == (ℵ₀, ℵ₀)
+            @test axes(B) == (OneToInf(), OneToInf())
+            @test eltype(B) == Float64
+            for _B in (B, B')
+                BB = copy(_B)
+                @test BB.dv.data === BB.ev.data
+                @test parent(BB).dv.data.datasize == parent(_B).dv.data.datasize
+                # @test !(BB === B) && !(parent(BB).dv.data === parent(B).dv.data) # copy is a no-op
+                @test BB[1:100, 1:100] == _B[1:100, 1:100]
+                @test BB[1:2:50, 1:3:40] == _B[1:2:50, 1:3:40]
+                @test view(BB, [1, 3, 7, 10], 1:10) == _B[[1, 3, 7, 10], 1:10]
+            end
+            @test LazyBandedMatrices.bidiagonaluplo(B) == uplo
+            @test LazyBandedMatrices.Bidiagonal(B) === LazyBandedMatrices.Bidiagonal(B.dv, B.ev, Symbol(uplo))
+            @test B[1:10, 1:10] ≈ A[1:10, 1:10]
+            @test B[230, 230] ≈ A[230, 230]
+            @test B[102, 102] ≈ A[102, 102] # make sure we compute intermediate columns correctly when skipping 
+            @test B[band(0)][1:100] == B.dv[1:100]
+            if uplo == 'U'
+                @test B[band(1)][1:100] == B.ev[1:100]
+                # @test B[band(-1)][1:100] == zeros(100) # This test requires that we define a 
+                # convert(::Type{BidiagonalConjugationBand{T}}, ::Zeros{V, 1, Tuple{OneToInf{Int}}}) where {T, V} method, 
+                # which we probably don't need beyond this test
+            else
+                @test B[band(-1)][1:100] == B.ev[1:100]
+                # @test B[band(1)][1:100] == zeros(100)
+            end
+            @test B.dv[500] == B.dv.data.dv[500]
+            @test B.dv.data.datasize == 501
+            @test B.ev[1005] == B.ev.data.ev[1005]
+            @test B.ev.data.datasize == 1006
+            @test ApplyArray(inv, B)[1:100, 1:100] ≈ ApplyArray(inv, A)[1:100, 1:100] # need to somehow let inv (or even ApplyArray(inv, )) work
+            @test (B+B)[1:100, 1:100] ≈ 2A[1:100, 1:100] ≈ 2B[1:100, 1:100]
+            @test (B*I)[1:100, 1:100] ≈ B[1:100, 1:100]
+            # @test (B*Diagonal(1:∞))[1:100, 1:100] ≈ B[1:100, 1:100] * Diagonal(1:100) # Uncomment once https://github.com/JuliaLinearAlgebra/ArrayLayouts.jl/pull/241 is registered
+
+            # Pointwise tests 
+            for i in 1:10
+                for j in 1:10
+                    @test B[i, j] ≈ A[i, j]
+                    @test B'[i, j] ≈ A[j, i]
+                end
+            end
+            @inferred B[5, 5]
+
+            # Make sure that, when indexing the transpose, B expands correctly
+            @test B'[3000:3005, 2993:3006] ≈ A[2993:3006, 3000:3005]'
+        end
+    end
+end