From 6e5619552c93ae91124cd31c4fb1d3d1b8495dac Mon Sep 17 00:00:00 2001 From: Daniel VandenHeuvel <95613936+DanielVandH@users.noreply.github.com> Date: Sun, 14 Jul 2024 10:56:51 +0100 Subject: [PATCH] Implement BidiagonalConjugationData (#186) * Implement BidiagonalConjugationData * Redesign * Fix import and make infiniterandomarrays an extra * typo * Tests * Redesign * Function barriers * BidiagonalConjugationBand should not be <: AbstractCachedVector * oops * Changes * Make copy a no-op * Cleanup * Fix resize --------- Co-authored-by: Sheehan Olver --- Project.toml | 8 +- src/InfiniteLinearAlgebra.jl | 7 +- src/banded/bidiagonalconjugation.jl | 138 ++++++++++++++++++++++++++++ test/runtests.jl | 8 +- test/test_bidiagonalconjugation.jl | 90 ++++++++++++++++++ 5 files changed, 243 insertions(+), 8 deletions(-) create mode 100644 src/banded/bidiagonalconjugation.jl create mode 100644 test/test_bidiagonalconjugation.jl diff --git a/Project.toml b/Project.toml index 030395a..7d88418 100644 --- a/Project.toml +++ b/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"] diff --git a/src/InfiniteLinearAlgebra.jl b/src/InfiniteLinearAlgebra.jl index 19203f3..175f1c2 100644 --- a/src/InfiniteLinearAlgebra.jl +++ b/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 diff --git a/src/banded/bidiagonalconjugation.jl b/src/banded/bidiagonalconjugation.jl new file mode 100644 index 0000000..57eaa06 --- /dev/null +++ b/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) \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 2fbaece..4481369 100644 --- a/test/runtests.jl +++ b/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") \ No newline at end of file diff --git a/test/test_bidiagonalconjugation.jl b/test/test_bidiagonalconjugation.jl new file mode 100644 index 0000000..510d180 --- /dev/null +++ b/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