Support A[:,3:end] when A isa InfToeplitz (#2)

* Support A[:,3:end] when A isa InfToeplitz

* General pert Toeplitz QL

* fix submaterialize

* fix banded Toeplitz

* resave

* improve q

* tests pass again

Author: dlfivefifty

Project.toml

@@ -14,7 +14,7 @@ MatrixFactorizations = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
 
 
 [compat]
-BandedMatrices = "0.12.2"
+BandedMatrices = "0.12.3"
 BlockArrays = "0.10"
 BlockBandedMatrices = "0.5.1"
 FillArrays = "0.7.1"

examples/toeplitz.jl

@@ -1,7 +1,5 @@
-using Revise, InfiniteLinearAlgebra, BlockBandedMatrices, BlockArrays, BandedMatrices, LazyArrays, FillArrays, MatrixFactorizations, Plots
-import MatrixFactorizations: reflectorApply!, QLPackedQ
-import InfiniteLinearAlgebra: blocktailiterate, _ql, qltail, rightasymptotics
-import BandedMatrices: bandeddata,_BandedMatrix
+using InfiniteLinearAlgebra, BandedMatrices, PyPlot
+
 
 function ℓ11(A,λ; kwds...) 
     try 
@@ -16,24 +14,23 @@ function findsecond(λ)
     λ[j], j
 end
 
-function qlplot(A; branch=findmax, x=range(-4,4; length=200), y=range(-4,4;length=200), kwds...)
+qlplot(A::AbstractMatrix; kwds...) = qlplot(BandedMatrix(A); kwds...)
+function qlplot(A::BandedMatrix; branch=findmax, x=range(-4,4; length=200), y=range(-4,4;length=200), kwds...)
     z = ℓ11.(Ref(A), x' .+ y.*im; branch=branch)
     contourf(x,y,z; kwds...)
 end
 
-function symbolplot!(A::BandedMatrix; kwds...)
+
+
+toepcoeffs(A::BandedMatrix) = InfiniteLinearAlgebra.rightasymptotics(A.data).args[1]
+toepcoeffs(A::Adjoint) = reverse(toepcoeffs(A'))
+
+function symbolplot(A::BandedMatrix; kwds...)
     l,u = bandwidths(A)
-    a = rightasymptotics(A.data).args[1]
+    a = toepcoeffs(A)
     θ = range(0,2π; length=1000)
-    i = Vector{ComplexF64}()
-    for t in θ
-        r = 0.0+0.0im
-        for k=1:length(a)
-            r += exp(im*(k-l-1)*t) * a[k]
-        end
-        push!(i,r)
-    end
-    plot!(i; kwds...)
+    i = map(t -> dot(exp.(im.*(u:-1:-l).*t),a),θ)
+    plot(real.(i), imag.(i); kwds...)
 end
 
 
@@ -42,212 +39,30 @@ end
 ###
 
 A = BandedMatrix(1 => Fill(2im,∞), 2 => Fill(-1,∞), 3 => Fill(2,∞), -2 => Fill(-4,∞), -3 => Fill(-2im,∞))
-p =plot(); symbolplot!(A)
+clf();qlplot(A; x=range(-10,7; length=100), y=range(-7,8;length=100)); symbolplot(A; color=:black); title("Trefethen & Embree")
+clf(); qlplot(transpose(A); x=range(-10,7; length=200), y=range(-7,8;length=200)); symbolplot(A; color=:black); title("Trefethen & Embree, transpose")
 
 ###
-# Non-normal
+# limaçon
 ###
 
-A = BandedMatrix(-1 => Vcat(Float64[], Fill(1/4,∞)), 0 => Vcat([1im],Fill(0,∞)), 1 => Vcat(Float64[], Fill(1,∞)))
-qlplot(A; title="A", linewidth=0, x=range(-2,2; length=100), y=range(-2,2;length=100))
-symbolplot!(A; linewidth=2.0, linecolor=:blue, legend=false)
-qlplot(BandedMatrix(A'); title="A', 1st", linewidth=0, nlevels=100, x=range(-2,2; length=100), y=range(-2,2;length=100))
-symbolplot!(A; linewidth=2.0, linecolor=:blue, legend=false)
-qlplot(BandedMatrix(A'); branch=findsecond, title="A'", linewidth=0, nlevels=100, x=range(-2,2; length=100), y=range(-2,2;length=100))
-
-
-heatmap!(x,y,fill(-100,length(y),length(x)); legend=false, color=:grays, fillalpha=(z -> isnan(z) ? 0.0 : 1.0).(z))
-
-Matrix((A)[1:1000,1:1000]) |> eigvals |> scatter
-
-ql(A-2I)
-
-
-θ = range(0,2π; length=1000)
-a = z -> z + 0.25/z
-plot!(a.(exp.(im.*θ)); linewidth=2.0, linecolor=:blue, legend=false)
-
-
-
-
-
-A = BandedMatrix(-10 => Fill(4.0,∞), 1 => Fill(1,∞))
-ql(A+0im*I)
-
-
-
-
-θ = range(0,2π-0.5; length=1000)
-a = z -> 4z^10 + 1/z
-θ = range(0,2π; length=1000)
-plot!(a.(exp.(im.*θ)); linewidth=2.0, linecolor=:blue, legend=false)
-
-import InfiniteLinearAlgebra: tail_de
-a = reverse(A.data.args[1]) .+ 0im
-
-a = randn(10) .+ im*randn(10); a[1] += 10; a[end] = 1;
-de = tail_de(a)
-
-@which tail_de(a)
-
-ql([transpose(a); 0 transpose(de)])
-
-ql([transpose(a); 0 transpose(de2)])
-
-de
-
-contourf(x,y,angle.(a.(x' .+ y.*im)))
-
-A'
-Fun(a, Laurent())
-
-A'
-
-
-a(1.0exp(0.5im))
+A = BandedMatrix(-2 => Fill(1.0,∞), -1 => Fill(1.0,∞), 1 => Fill(eps(),∞))
+qlplot(A; x=range(-2,3; length=100), y=range(-2.5,2.5;length=100)); symbolplot(A; color=:black); title("Limacon")
+clf(); qlplot(A; branch=findsecond, x=range(-2,3; length=100), y=range(-2.5,2.5;length=100)); symbolplot(A; color=:black); title("Limacon, second branch")
+clf(); qlplot(transpose(A); x=range(-2,3; length=100), y=range(-2.5,2.5;length=100)); symbolplot(A; color=:black); title("Limacon, transpose")
 
-#### To be cleaned
 
+ql(A-(0.5+0.000001im)*I; branch=findsecond)
 
-function reduceband(A)
-        l,u = bandwidths(A)
-        H = _BandedMatrix(A.data, ∞, l+u-1, 1)
-        Q1,L1 = ql(H)
-        D = Q1[1:l+u+1,1:1]'A[1:l+u+1,1:u-1]
-        D, Q1, L1
-end
-_Lrightasymptotics(D::Vcat) = D.args[2]
-_Lrightasymptotics(D::ApplyArray) = D.args[1][2:end] * Ones{ComplexF64}(1,∞)
-Lrightasymptotics(L) = _Lrightasymptotics(rightasymptotics(parent(L).data))
-
-function qdL(A)
-    l,u = bandwidths(A)
-    H = _BandedMatrix(A.data, ∞, l+u-1, 1)
-    Q1,L1 = ql(H)
-    D1, Q1, L1 = reduceband(A)
-    T2 = _BandedMatrix(Lrightasymptotics(L1), ∞, l, u)
-    l1 = L1[1,1]
-    A2 = [[D1 l1 zeros(1,10-size(D1,2)-1)]; T2[1:10-1,1:10]] # TODO: remove
-    B2 = _BandedMatrix(T2.data, ∞, l+u-2, 2)
-    B2 = _BandedMatrix(T2.data, ∞, l+u-2, 2)
-    D2, Q2, L2 = reduceband(B2)
-    l2 = L2[1,1]
-    # peroidic tail
-    T3 = _BandedMatrix(Lrightasymptotics(L2), ∞, l+1, u-1)
-    A3 = [[D2 l2 zeros(1,10-size(D2,2)-1)]; T3[1:10-1,1:10]] # TODO: remove
-
-    Q3,L3 = ql( [A2[1,1] A2[1:1,2:3]; [Q2[1:3,1:1]' * T2[1:3,1]  A3[1:1,1:2] ]])
-
-    fd_data = hcat([0; L3[:,1]; Q2[1:3,2:3]' * T2[1:3,1]], [L3[:,2]; T3[1:3,1]], [L3[2,3]; T3[1:4,2]])
-    B3 = _BandedMatrix(Hcat(fd_data, T3.data), ∞, l+u-1, 1)
-
-    ql(B3).L
-end
-# Bull's head
-A = BandedMatrix(-3 => Fill(7/10,10), -2 => Fill(1,11), 1 => Fill(2im,9))
+### 
+# bull-head
+###
 
 A = BandedMatrix(-3 => Fill(7/10,∞), -2 => Fill(1,∞), 1 => Fill(2im,∞))
-qlplot(A;  title="largest", linewidth=0)
-
-
-qlplot(A; branch=findsecond, title="second largest", linewidth=0)
-
-
-
-sortp
-
-ql(A-(1+2im)*I)
-
-θ = range(0,2π; length=1000)
-a = z -> 2im/z + z^2 + 7/10 * z^3
-plot!(a.(exp.(im.*θ)); linewidth=2.0, linecolor=:blue, legend=false)
-qlplot(A; branch=2, title="branch=2", linewidth=0)
-
-
-
+clf();qlplot(A; x=range(-10,7; length=100), y=range(-7,8;length=100)); symbolplot(A; color=:black); title("Bull-head")
+clf(); qlplot(transpose(A); x=range(-10,7; length=100), y=range(-7,8;length=100)); symbolplot(A; color=:black); title("Bull-head, transpose")
 
-θ = range(0,2π; length=1000)
-a = z -> 2im/z + z^2 + 7/10 * z^3
-plot!(a.(exp.(im.*θ)); linewidth=2.0, linecolor=:blue, legend=false)
-
-
-ql(A - (5+2im)*I; branch=1)
-
-ℓ = λ -> try abs(ql(A-λ*I).L[1,1]) catch DomainError 
-            (-1) end
-x,y = range(-4,4; length=200),range(-4,4;length=200)
-z = ℓ.(x' .+ y.*im)
-contourf(x,y,z; nlevels=100, title="A", linewidth=0)
-θ = range(0,2π; length=1000)
-a = z -> 2im/z + z^2 + 7/10 * z^3
-plot!(a.(exp.(im.*θ)); linewidth=2.0, linecolor=:blue, legend=false)
-
-Ac = BandedMatrix(A')
-
-ℓ = λ -> try abs(qdL(Ac-conj(λ)*I)[1,1]) catch DomainError 
-            (-1.0) end
-x,y = range(-4,4; length=100),range(-4,4;length=100)
-@time z = ℓ.(x' .+ y.*im)
-contourf(x,y,z; nlevels=100, title="A'", linewidth=0.0)
-θ = range(0,2π; length=1000)
-a = z -> 2im/z + z^2 + 7/10 * z^3
-plot!(a.(exp.(im.*θ)); linewidth=2.0, linecolor=:blue, legend=false)
 
+###
 # Grcar
-A = BandedMatrix(-3 => Fill(1,∞), -2 => Fill(1,∞), -1 => Fill(1,∞), 0 => Fill(1,∞), 1 => Fill(-1,∞))
-ℓ = λ -> try abs(ql(A-λ*I).L[1,1]) catch DomainError 
-        -1 end
-x,y = range(-4,4; length=200),range(-4,4;length=200)
-z = ℓ.(x' .+ y.*im)
-
-θ = range(0,2π; length=1000)
-a = z -> -1/z + 1 + z + z^2 + z^3
-contourf(x,y,z; nlevels=50)
-plot!(a.(exp.(im.*θ)); linewidth=2.0, linecolor=:blue, legend=false)
-
-
-# Triangle
-A = BandedMatrix(-2 => Fill(1/4,∞), 1 => Fill(1,∞))
-ℓ = λ -> try abs(ql(A-λ*I).L[1,1]) catch DomainError 
-            (-1) end
-x,y = range(-2,2; length=200),range(-2,2;length=200)
-z = ℓ.(x' .+ y.*im)
-θ = range(0,2π; length=1000)
-a = z -> z^2/4 + 1/z
-
-contourf(x,y,z; nlevels=50)
-plot!(a.(exp.(im.*θ)); linewidth=2.0, linecolor=:blue, legend=false)
-
-
-function Toep_L11(T)
-        l,u = bandwidths(T)
-        @assert u == 2
-        # shift by one
-        H = _BandedMatrix(T.data, ∞, l+1, 1)
-        Q1,L1 = ql(H)
-
-        d = Q1[1:3,1]'T[1:1+l,1]
-        ℓ = Q1.factors.data.args[2].args[1][2:end] # new L
-        T2 = _BandedMatrix(Hcat([[zero(d); d; ℓ[3:end]] L1[1:5,1]], ℓ*Ones{eltype(T)}(1,∞)), ∞, 3, 1)
-        Q2,L2 = ql(T2)
-        D = (Q2')[1:5,1:4] * (Q1')[1:4,1:3] * T[3:5,1:3]
-        X = [Matrix(T[3:4,1:6]); [zeros(2,2) [-1 0; 0 1]*D[1:2,:] [-1 0; 0 1]*L2[1:2,2]]]
-        ql(X).L[1,3]
-end
-
-
-
-A = BandedMatrix(2 => Fill(1/4,∞), -1 => Fill(1,∞), -2 => Fill(0.0, ∞))
-@time Toep_L11(A-(0.1+0im)I)
-
-T = A-(0.1+0im)I
-
-
-ℓ = λ -> abs(Toep_L11(A-λ*I)) 
-x,y = range(-2,2; length=51),range(-2,2;length=50)
-z = ℓ.(x' .+ y.*im)
-θ = range(0,2π; length=1000)
-a = z -> z^2/4 + 1/z
-
-contour(x,y,z; nlevels=50)
-plot!(a.(exp.(im.*θ)); linewidth=2.0, linecolor=:blue, legend=false)
+###

src/InfiniteLinearAlgebra.jl

@@ -2,22 +2,23 @@ module InfiniteLinearAlgebra
 using BlockArrays, BlockBandedMatrices, BandedMatrices, LazyArrays, FillArrays, InfiniteArrays, MatrixFactorizations, LinearAlgebra
 
 import Base: +, -, *, /, \, ^, OneTo, getindex, promote_op, _unsafe_getindex, print_matrix_row, size, axes,
-            AbstractMatrix, AbstractArray, Matrix, Array, Vector, AbstractVector,
+            AbstractMatrix, AbstractArray, Matrix, Array, Vector, AbstractVector, Slice,
             show, getproperty
 import Base.Broadcast: BroadcastStyle
 
-import InfiniteArrays: OneToInf, InfUnitRange, Infinity, InfStepRange
+import InfiniteArrays: OneToInf, InfUnitRange, Infinity, InfStepRange, AbstractInfUnitRange
 import FillArrays: AbstractFill, getindex_value
-import BandedMatrices: BandedMatrix, _BandedMatrix, bandeddata, bandwidths
+import BandedMatrices: BandedMatrix, _BandedMatrix, AbstractBandedMatrix, bandeddata, bandwidths, BandedColumns, bandedcolumns
 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, factorize
+                    colsupport, rowsupport, triangularlayout, factorize, subarraylayout, sub_materialize,
+                    @lazymul
 import MatrixFactorizations: ql, ql!, QLPackedQ, getL, getR, reflector!, reflectorApply!, QL, QR, QRPackedQ
 
 import BlockArrays: BlockSizes, cumulsizes, _find_block, AbstractBlockVecOrMat, sizes_from_blocks
 
-import BandedMatrices: BandedMatrix, bandwidths, AbstractBandedLayout, _banded_qr!, _banded_qr
+import BandedMatrices: BandedMatrix, bandwidths, AbstractBandedLayout, _banded_qr!, _banded_qr, _BandedMatrix
 
 import BlockBandedMatrices: _BlockSkylineMatrix, _BandedMatrix, AbstractBlockSizes, cumulsizes, _BlockSkylineMatrix, BlockSizes, blockstart, blockstride,
         BlockSkylineSizes, BlockSkylineMatrix, BlockBandedMatrix, _BlockBandedMatrix, BlockTridiagonal

src/banded/hessenbergq.jl

@@ -90,14 +90,24 @@ adjoint(Q::LowerHessenbergQ) = UpperHessenbergQ(adjoint.(Q.q))
 check_mul_axes(A::AbstractHessenbergQ, B, C...) =
     axes(A,2) == axes(B,1) || throw(DimensionMismatch("Second axis of A, $(axes(A,2)), and first axis of B, $(axes(B,1)) must match"))
 
+@lazymul AbstractHessenbergQ
+
+# ambiguities
+for Arr in (:AbstractBandedMatrix, :StridedVector, :StridedMatrix)
+    @eval begin
+        *(A::AbstractHessenbergQ, B::$Arr) = apply(*, A, B)
+        *(A::$Arr, B::AbstractHessenbergQ) = apply(*, A, B)
+    end
+end
 
 function lmul!(Q::LowerHessenbergQ{T}, x::AbstractVector) where T
     t = Array{T}(undef, 2)
+    nz = nzzeros(x,1)
     for n = 1:length(Q.q)
         v = view(x, n:n+1)
         mul!(t, Q.q[n], v)
         v .= t
-        all(iszero,t) && return x
+        n > nz && norm(t) ≤ 10floatmin(real(T)) && return x
     end
     x
 end
@@ -112,6 +122,13 @@ function lmul!(Q::UpperHessenbergQ{T}, x::AbstractVector) where T
     x
 end
 
+function lmul!(Q::AbstractHessenbergQ, X::AbstractMatrix)
+    for j in axes(X,2)
+        lmul!(Q, view(X,:,j))
+    end
+    X
+end
+
 
 ###
 # getindex