-
Notifications
You must be signed in to change notification settings - Fork 7
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
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 <[email protected]>
- Loading branch information
1 parent
3e81f50
commit 6e56195
Showing
5 changed files
with
243 additions
and
8 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 |
6e56195There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@JuliaRegistrator register
6e56195There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Registration pull request updated: JuliaRegistries/General/109177
Tip: Release Notes
Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.
To add them here just re-invoke and the PR will be updated.
Tagging
After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.
This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via: