QL for BBBArrowheadMatrix (WIP)

src/semiseparable_arrowhead.jl

@@ -0,0 +1,41 @@
+struct SemiseparableBBBArrowheadMatrix{T} <: AbstractBlockBandedMatrix{T}
+    # banded parts
+    A::BandedMatrix{T}
+    B::NTuple{2,BandedMatrix{T}} # first row blocks
+    C::NTuple{4,BandedMatrix{T}} # first col blocks
+    D
+
+    # fill parts
+    Asub::NTuple{2,Vector{T}}
+    Asup::Tuple{2,Matrix{T}} # matrices are m × 2
+
+    Bsub::NTuple{2,Vector{T}}
+    Bsup::NTuple{2,NTuple{2,Vector{T}}}
+
+    Csub::NTuple{2,NTuple{2,Vector{T}}}
+    Csup::NTuple{2,Vector{T}}
+
+    A22sub::NTuple{2,Vector{T}}
+    A32sub::NTuple{2,Vector{T}}
+
+    A32extra::Vector{T}
+    A33extra::Vector{T}
+
+    D::DD # these are interlaces
+
+end
+
+axes(::SemiseparableBBBArrowheadMatrix) = ...
+
+function getindex(L::SemiseparableBBBArrowheadMatrix{T}, Kk::BlockIndex{1}, Jj::BlockIndex{1})::T where T
+    K,k = block(Kk),blockindex(Kk)
+    J,j = block(Jj),blockindex(Jj)
+    # TODO: add getindex
+end
+
+
+function getindex(L::SemiseparableBBBArrowheadMatrix, k::Int, j::Int)
+    ax,bx = axes(L)
+    L[findblockindex(ax, k), findblockindex(bx, j)]
+end
+

test/explore_QL.jl

@@ -0,0 +1,114 @@
+using PiecewiseOrthogonalPolynomials, Plots, BlockArrays, Test
+using MatrixFactorizations, LinearAlgebra, BlockBandedMatrices
+###
+# QL
+####
+function my_ql(A::BBBArrowheadMatrix{T}) where T
+    m,n = size(A.A)
+    l = length(A.D)
+    m2, n2 = size(A.D[1])
+    @assert m == n == l+1
+    @assert m2 == n2
+    #results stored in F and τ
+    F = BlockedArray(Matrix(A), axes(A))
+    τ = zeros(m+l*m2)
+    for j in m2:-1:3
+        for i in l:-1:1
+            upper_entry = F[Block(j-1, j+1)][i, i] #A.D[i][j-2,j]
+            dia_entry = F[Block(j+1, j+1)][i, i] #A.D[i][j,j]
+            #perform Householder transformation
+            dia_entry_new = -sign(dia_entry)*sqrt(dia_entry^2 + upper_entry^2)
+            v = [upper_entry, dia_entry-dia_entry_new]
+            coef = 2/(v[1]^2+v[2]^2)
+            #denote the householder transformation as [c1 s1;c2 s2]
+            c1 = 1 - coef * v[1]^2
+            s1 = - coef * v[1] * v[2]
+            c2 = s1
+            s2 = 1 - coef * v[2]^2
+            print(dia_entry_new)
+            F[m+(j-1)*l+i, m+(j-1)*l+i] = dia_entry_new #update F[Block(j+1, j+1)][i, i]
+            F[m+(j-3)*l+i, m+(j-1)*l+i] = v[1]/v[2] #update F[Block(j-1, j+1)][i, i]
+            τ[m+(j-1)*l+i] = coef*v[2]^2
+            #row recombination(householder transformation) for other columns
+            current_upper_entry = F[Block(j-1, j-1)][i, i] #A.D[i][j-2,j-2]
+            current_lower_entry = F[Block(j+1, j-1)][i, i] #A.D[i][j,j-2]
+            F[m+(j-3)*l+i, m+(j-3)*l+i] = c1 * current_upper_entry + s1 * current_lower_entry #update F[Block(j-1, j-1)][i, i]
+            F[m+(j-1)*l+i, m+(j-3)*l+i] = c2 * current_upper_entry + s2 * current_lower_entry #update F[Block(j+1, j-1)][i, i]
+            if j >= 5
+                #Deal with A.D blocks which do not share common rows with A.C
+                current_entry = F[Block(j-1, j-3)][i, i] #A.D[i][j-2,j-4]
+                F[m+(j-3)*l+i, m+(j-5)*l+i] = c1 * current_entry #update F[Block(j-1, j-3)][i, i]
+                F[m+(j-1)*l+i, m+(j-5)*l+i] = c2 * current_entry #update F[Block(j+1, j-3)][i, i]
+            else
+                #Deal with A.D blocks which share common rows with A.C
+                current_entry = F[Block(j-1, 1)][i, i] #A.C[j-2][i,i]
+                F[m+(j-3)*l+i, i] = c1 * current_entry #update F[Block(j-1, 1)][i, i]
+                F[m+(j-1)*l+i, i] = c2 * current_entry #update F[Block(j+1, 1)][i, i]
+
+                current_entry = F[Block(j-1, 1)][i, i+1] #A.C[j-2][i,i+1]
+                F[m+(j-3)*l+i, i+1] = c1 * current_entry #update F[Block(j-1, 1)][i, i+1]
+                F[m+(j-1)*l+i, i+1] = c2 * current_entry #F[Block(j+1, 1)][i, i+1]
+            end
+        end
+    end
+
+    #Deal with Block(1,3)
+    #vectors x and Λ denote a rank 1 semiseperable matrix
+    λ = 1.0
+    Λ = []
+    x = [F[Block(1,3)][l+1,l]]
+    x_len = abs(x[1])
+    for i in l:-1:2 #consider i=1 later
+        a = F[Block(1,3)][i,i]
+        b = F[Block(1,3)][i,i-1]
+        c = F[Block(3,3)][i,i]
+        v_last = c + sign(c) * sqrt(a^2 + λ^2 * x_len^2 + c^2)
+        v_len = sqrt(a^2 + λ^2 * x_len^2 + v_last^2)
+        F[m+l+i,m+l+i] = -sign(c) * sqrt(a^2 + λ^2 * x_len^2 + c^2)
+        pushfirst!(Λ, λ / v_last)
+        λ = -2/v_len^2 * a * b * λ
+        F[m+l+i, m+l+i-1] = -2/v_len^2 * v_last * a * b
+        x_first = (1 - 2/v_len^2 * a^2) * b / λ
+        pushfirst!(x, x_first)  
+        x_len = sqrt(x_len^2 + x_first^2)  
+        #record information of V
+        F[i+1, m+l+i] = 0
+        F[i, m+l+i] = a / v_last
+        τ[m+l+i] = 2 * v_last^2 / v_len^2
+    end
+    #deal with the last column in Block(1,3)
+    a = F[Block(1,3)][1,1]
+    c = F[Block(3,3)][1,1]
+    v_last = c + sign(c) * sqrt(a^2 + λ^2 * x_len^2 + c^2)
+    v_len = sqrt(a^2 + λ^2 * x_len^2 + v_last^2)
+    pushfirst!(Λ, λ / v_last)
+    F[m+l+1,m+l+1] = -sign(c) * sqrt(a^2 + λ^2 * x_len^2 + c^2)
+    F[2, m+l+1] = 0
+    F[1, m+l+1] = a / v_last
+    τ[m+l+1] = 2 * v_last^2 / v_len^2
+
+    F, τ, x, Λ
+end
+
+
+𝐗 = range(-1,1; length=10)
+C = ContinuousPolynomial{1}(𝐗)
+plot(C[:,Block(2)])
+
+#plot(C[:,Block.(2:3)])
+M = C'C
+#M = grammatrix(C)
+Δ = weaklaplacian(C)
+N = 6
+KR = Block.(Base.OneTo(N))
+Mₙ = M[KR,KR]
+Δₙ = Δ[KR,KR]
+A = Δₙ + 100^2 * Mₙ
+FF,tτ, xx, LΛ = my_ql(A)
+τ = ql(A).τ
+f = ql(A).factors
+
+@test BlockedArray(tτ, (axes(A,2),))[Block.(3:6)] ≈ BlockedArray(τ, (axes(A,2),))[Block.(3:6)]
+@test f[:,Block.(4:6)] ≈ FF[:,Block.(4:6)]
+
+@test tril(xx * LΛ') + FF[Block(1,3)][2:end,:] ≈ f[Block(1,3)][2:end,:]