BesselJ example

Author: dlfivefifty

Project.toml

@@ -26,6 +26,7 @@ julia = "1"
 [extras]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [targets]
-test = ["Test", "Random"]
+test = ["Test", "Random", "SpecialFunctions"]

README.md

@@ -7,6 +7,63 @@ A Julia repository for linear algebra with infinite banded and block-banded matr
 
 This currently supports the infinite-dimensional QR factorization for banded matrices, also known as the adaptive QR decomposition as the entries of the QR decomposition are determined lazily. 
 
+As a simple example, consider the Bessel recurrence relationship:
+$$
+J_{n-1}(z)  - {2 n \over z} J_n(z) + J_{n+1}(z) = 0
+$$
+This can be recast as an infinite linear system:
+```julia
+julia> using InfiniteLinearAlgebra, InfiniteArrays, BandedMatrices, FillArrays, SpecialFunctions
+
+julia> z = 10_000; # the bigger z the longer before we see convergence
+
+julia> A = BandedMatrix(0 => -2*(0:∞)/z, 1 => Ones(∞), -1 => Ones(∞))
+∞×∞ BandedMatrix{Float64,ApplyArray{Float64,2,typeof(*),Tuple{Array{Float64,2},ApplyArray{Float64,2,typeof(vcat),Tuple{Transpose{Float64,InfiniteArrays.InfStepRange{Float64,Float64}},Ones{Float64,2,Tuple{Base.OneTo{Int64},InfiniteArrays.OneToInf{Int64}}},Ones{Float64,2,Tuple{Base.OneTo{Int64},InfiniteArrays.OneToInf{Int64}}}}}}},InfiniteArrays.OneToInf{Int64}} with indices OneToInf()×OneToInf():
+ 0.0   1.0       ⋅        ⋅        ⋅        ⋅       ⋅        ⋅        ⋅      …  
+ 1.0  -0.0002   1.0       ⋅        ⋅        ⋅       ⋅        ⋅        ⋅         
+  ⋅    1.0     -0.0004   1.0       ⋅        ⋅       ⋅        ⋅        ⋅         
+  ⋅     ⋅       1.0     -0.0006   1.0       ⋅       ⋅        ⋅        ⋅         
+  ⋅     ⋅        ⋅       1.0     -0.0008   1.0      ⋅        ⋅        ⋅         
+  ⋅     ⋅        ⋅        ⋅       1.0     -0.001   1.0       ⋅        ⋅      …  
+  ⋅     ⋅        ⋅        ⋅        ⋅       1.0    -0.0012   1.0       ⋅         
+  ⋅     ⋅        ⋅        ⋅        ⋅        ⋅      1.0     -0.0014   1.0        
+  ⋅     ⋅        ⋅        ⋅        ⋅        ⋅       ⋅       1.0     -0.0016     
+  ⋅     ⋅        ⋅        ⋅        ⋅        ⋅       ⋅        ⋅       1.0        
+  ⋅     ⋅        ⋅        ⋅        ⋅        ⋅       ⋅        ⋅        ⋅      …  
+  ⋅     ⋅        ⋅        ⋅        ⋅        ⋅       ⋅        ⋅        ⋅         
+  ⋅     ⋅        ⋅        ⋅        ⋅        ⋅       ⋅        ⋅        ⋅         
+  ⋅     ⋅        ⋅        ⋅        ⋅        ⋅       ⋅        ⋅        ⋅         
+  ⋅     ⋅        ⋅        ⋅        ⋅        ⋅       ⋅        ⋅        ⋅         
+ ⋮                                         ⋮                                 ⋱  
+```
+The first row corresponds to specifying an initial condition. Thus we can determine $J_n(z)$ via solving the recurrence:
+```julia
+julia> A \ Vcat([besselj(1,z)], Zeros(∞)) 
+∞-element LazyArrays.CachedArray{Float64,1,Array{Float64,1},Zeros{Float64,1,Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf():
+ -0.007096160353406478 
+  0.0036474507555295833
+  0.007096889843557584 
+ -0.0036446119995921654
+ -0.007099076610757337 
+  0.0036389327383035616
+  0.00710271554349564  
+ -0.003630409479651369 
+ -0.007107798116767152 
+  0.0036190370026645442
+  0.007114312383371949 
+ -0.0036048083778978026
+ -0.0071222429618033245
+  0.0035877149947894766
+  0.007131571020789777 
+  ⋮                    
+
+julia> J[1000] - besselj(999,z) # matches besselj to high (relative) accuracy
+-6.8252695162307475e-15
+
+julia> J[11_000] - besselj(11_000-1, z)
+3.3730094946097293e-143
+```
+
 ## Infinite-dimensional QL factorization
 
 

src/InfiniteLinearAlgebra.jl

@@ -12,7 +12,7 @@ import BandedMatrices: BandedMatrix, _BandedMatrix, bandeddata, bandwidths
 import LinearAlgebra: lmul!,  ldiv!, matprod, qr, AbstractTriangular, AbstractQ, adjoint, transpose
 import LazyArrays: CachedArray, CachedMatrix, CachedVector, DenseColumnMajor, FillLayout, ApplyMatrix, check_mul_axes, ApplyStyle, LazyArrayApplyStyle, LazyArrayStyle,
                     CachedMatrix, CachedArray, resizedata!, MemoryLayout, mulapplystyle, LmulStyle, RmulStyle,
-                    colsupport, rowsupport, triangularlayout
+                    colsupport, rowsupport, triangularlayout, factorize
 import MatrixFactorizations: ql, ql!, QLPackedQ, getL, getR, reflector!, reflectorApply!, QL, QR, QRPackedQ
 
 import BlockArrays: BlockSizes, cumulsizes, _find_block, AbstractBlockVecOrMat, sizes_from_blocks

src/banded/infbanded.jl

@@ -8,9 +8,22 @@ const TriPertToeplitz{T} = Tridiagonal{T,Vcat{T,1,Tuple{Vector{T},Fill{T,1,Tuple
 const AdjTriPertToeplitz{T} = Adjoint{T,Tridiagonal{T,Vcat{T,1,Tuple{Vector{T},Fill{T,1,Tuple{OneToInf{Int}}}}}}}
 const InfBandedMatrix{T,V<:AbstractMatrix{T}} = BandedMatrix{T,V,OneToInf{Int}}
 
+function BandedMatrix{T}(kv::Tuple{Vararg{Pair{<:Integer,<:AbstractVector}}},
+                         ::NTuple{2,Infinity},
+                         (l,u)::NTuple{2,Integer}) where T
+    ks = getproperty.(kv, :first)
+    rws = Vcat(permutedims.(getproperty.(kv, :second))...)
+    l,u = -minimum(ks),maximum(ks)
+    c = zeros(T, l+u+1, length(ks))
+    for (k,j) in zip(u .- ks .+ 1,1:length(ks))
+        c[k,j] = one(T)
+    end
+    _BandedMatrix(ApplyArray(*,c,rws), ∞, l, u)
+end
+
 # Construct InfToeplitz
 function BandedMatrix{T}(kv::Tuple{Vararg{Pair{<:Integer,<:Fill{<:Any,1,Tuple{OneToInf{Int}}}}}},
-                         mn::NTuple{2,Integer},
+                         mn::NTuple{2,Infinity},
                          lu::NTuple{2,Integer}) where T
     m,n = mn
     @assert isinf(n)
@@ -25,7 +38,7 @@ function BandedMatrix{T}(kv::Tuple{Vararg{Pair{<:Integer,<:Fill{<:Any,1,Tuple{On
 end
 
 function BandedMatrix{T}(kv::Tuple{Vararg{Pair{<:Integer,<:Vcat{<:Any,1,<:Tuple{<:AbstractVector,Fill{<:Any,1,Tuple{OneToInf{Int}}}}}}}},
-                         mn::NTuple{2,Integer},
+                         mn::NTuple{2,Infinity},
                          lu::NTuple{2,Integer}) where T
     m,n = mn
     @assert isinf(n)

src/infqr.jl

@@ -102,6 +102,7 @@ function lmul!(A::QRPackedQ{<:Any,<:AdaptiveQRFactors}, B::CachedVector{T,Vector
     b = B.data    
 
     sB_2 = length(b)
+    partialqr!(A.factors.data, sB)
 
     @inbounds begin
         for k = sB:-1:1
@@ -131,22 +132,30 @@ function lmul!(adjA::Adjoint{<:Any,<:QRPackedQ{<:Any,<:AdaptiveQRFactors}}, B::C
     end
     Afactors = A.factors
     sB = length(B.data)
+    partialqr!(A.factors.data, min(sB+100,nA))
+    resizedata!(B, min(sB+100,mB))
+    Bd = B.data
     @inbounds begin
         for k = 1:min(mA,nA)
             vBj = B[k]
             allzero = k > sB ? iszero(vBj) : false
             cs = colsupport(Afactors,k)
+            cs_max = maximum(cs)
+            if cs_max > length(Bd) # need to grow the data
+                resizedata!(B, min(cs_max+100,mB)); Bd = B.data
+                partialqr!(A.factors.data, min(cs_max+100,nA))
+            end
             for i = (k+1:mB) ∩ cs
-                Bi = B[i]
+                Bi = Bd[i]
                 if !iszero(Bi)
                     allzero = false
                     vBj += conj(Afactors[i,k])*Bi
                 end
             end
             vBj = conj(A.τ[k])*vBj
-            B[k] -= vBj
+            Bd[k] -= vBj
             for i = (k+1:mB) ∩ cs
-                B[i] -= Afactors[i,k]*vBj
+                Bd[i] -= Afactors[i,k]*vBj
             end
             allzero && break
         end

test/runtests.jl

@@ -2,32 +2,23 @@ using InfiniteLinearAlgebra, BlockBandedMatrices, BlockArrays, BandedMatrices, I
 import InfiniteLinearAlgebra: qltail, toeptail, tailiterate , tailiterate!, tail_de, ql_X!,
                     InfToeplitz, PertToeplitz, TriToeplitz, InfBandedMatrix, 
                     rightasymptotics, QLHessenberg
+
 import BlockArrays: _BlockArray                    
 import BlockBandedMatrices: isblockbanded, _BlockBandedMatrix
 import MatrixFactorizations: QLPackedQ
 import BandedMatrices: bandeddata, _BandedMatrix
 import LazyArrays: colsupport, ApplyStyle, MemoryLayout
 
-@testset "Algebra" begin 
-    @testset "BlockTridiagonal" begin
-        A = BlockTridiagonal(Vcat([fill(1.0,2,1),Matrix(1.0I,2,2),Matrix(1.0I,2,2),Matrix(1.0I,2,2)],Fill(Matrix(1.0I,2,2), ∞)), 
-                            Vcat([zeros(1,1)], Fill(zeros(2,2), ∞)), 
-                            Vcat([fill(1.0,1,2),Matrix(1.0I,2,2)], Fill(Matrix(1.0I,2,2), ∞)))
-                            
-        @test A isa InfiniteLinearAlgebra.BlockTriPertToeplitz                       
-        @test isblockbanded(A)
 
-        @test A[Block.(1:2),Block(1)] == A[1:3,1:1] == reshape([0.,1.,1.],3,1)
 
-        @test BlockBandedMatrix(A)[1:100,1:100] == BlockBandedMatrix(A,(2,1))[1:100,1:100] == BlockBandedMatrix(A,(1,1))[1:100,1:100] == A[1:100,1:100]
 
-        @test (A - I)[1:100,1:100] == A[1:100,1:100]-I
-        @test (A + I)[1:100,1:100] == A[1:100,1:100]+I
-        @test (I + A)[1:100,1:100] == I+A[1:100,1:100]
-        @test (I - A)[1:100,1:100] == I-A[1:100,1:100]
-    end
-    
+
+@testset "Algebra" begin 
     @testset "BandedMatrix" begin
+        A = BandedMatrix(-3 => Fill(7/10,∞), -2 => 1:∞, 1 => Fill(2im,∞))
+        @test A isa BandedMatrix{ComplexF64}
+        @test A[1:10,1:10] == diagm(-3 => Fill(7/10,7), -2 => 1:8, 1 => Fill(2im,9))
+
         A = BandedMatrix(-3 => Fill(7/10,∞), -2 => Fill(1,∞), 1 => Fill(2im,∞))
         Ac = BandedMatrix(A')
         At = BandedMatrix(transpose(A))
@@ -52,6 +43,25 @@ import LazyArrays: colsupport, ApplyStyle, MemoryLayout
         @test (A^3)[1:10,1:10] == (A*A*A)[1:10,1:10] == (A[1:100,1:100]^3)[1:10,1:10]
     end
 
+    @testset "BlockTridiagonal" begin
+        A = BlockTridiagonal(Vcat([fill(1.0,2,1),Matrix(1.0I,2,2),Matrix(1.0I,2,2),Matrix(1.0I,2,2)],Fill(Matrix(1.0I,2,2), ∞)), 
+                            Vcat([zeros(1,1)], Fill(zeros(2,2), ∞)), 
+                            Vcat([fill(1.0,1,2),Matrix(1.0I,2,2)], Fill(Matrix(1.0I,2,2), ∞)))
+                            
+        @test A isa InfiniteLinearAlgebra.BlockTriPertToeplitz                       
+        @test isblockbanded(A)
+
+        @test A[Block.(1:2),Block(1)] == A[1:3,1:1] == reshape([0.,1.,1.],3,1)
+
+        @test BlockBandedMatrix(A)[1:100,1:100] == BlockBandedMatrix(A,(2,1))[1:100,1:100] == BlockBandedMatrix(A,(1,1))[1:100,1:100] == A[1:100,1:100]
+
+        @test (A - I)[1:100,1:100] == A[1:100,1:100]-I
+        @test (A + I)[1:100,1:100] == A[1:100,1:100]+I
+        @test (I + A)[1:100,1:100] == I+A[1:100,1:100]
+        @test (I - A)[1:100,1:100] == I-A[1:100,1:100]
+    end
+    
+
     @testset "Fill" begin
         A = _BandedMatrix(Ones(1,∞),∞,-1,1)
         @test 1.0 .* A isa BandedMatrix{Float64,<:Fill}

test/test_infqr.jl

@@ -1,4 +1,4 @@
-using InfiniteLinearAlgebra, LinearAlgebra, BandedMatrices, InfiniteArrays, MatrixFactorizations, LazyArrays, FillArrays
+using InfiniteLinearAlgebra, LinearAlgebra, BandedMatrices, InfiniteArrays, MatrixFactorizations, LazyArrays, FillArrays, SpecialFunctions
 import LazyArrays: colsupport, rowsupport, MemoryLayout, DenseColumnMajor, TriangularLayout
 import BandedMatrices: _BandedMatrix, _banded_qr!, BandedColumns
 import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData, AdaptiveLayout
@@ -76,9 +76,22 @@ import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData, AdaptiveLayout
 
         @test Q'*b == r
 
+        @test factorize(A) isa typeof(qr(A))
         @test qr(A)\b == A\b
         @test (A*(A\b))[1:100] ≈ [1:3; Zeros(97)]
     end
+
+    @testset "Bessel J" begin
+        z = 1000; # the bigger z the longer before we see convergence
+        A = BandedMatrix(0 => -2*(0:∞)/z, 1 => Ones(∞), -1 => Ones(∞))
+        J = A \ Vcat([besselj(1,z)], Zeros(∞))
+        @test J[1:2000] ≈ [besselj(k,z) for k=0:1999]
+
+        z = 10_000; # the bigger z the longer before we see convergence
+        A = BandedMatrix(0 => -2*(0:∞)/z, 1 => Ones(∞), -1 => Ones(∞))
+        J = A \ Vcat([besselj(1,z)], Zeros(∞))
+        @test J[1:20_000] ≈ [besselj(k,z) for k=0:20_000-1]
+    end
 end