Add some functions for leftrightmul and truncatedconvop types

src/ConvolutionOperators.jl

@@ -25,6 +25,7 @@ include("linearcombinations.jl")
 include("leftrightmul.jl")
 include("liftedconvop.jl")
 include("truncatedconvop.jl")
+include("matrixconvop.jl")
 
 
 end

src/companion.jl

@@ -1,65 +1,98 @@
-function companion(Z)
+using LinearAlgebra
 
+```
+    Build the companion matrix for Z that has no tail
+```
+function companion_no_tail(Z)
+    M, N = size(Z)[1:2]
     T = eltype(Z)
-    K = size(Z,3)
+    @assert M == N
+
+    K = size(Z,3) - 1
     @assert K > 1
 
-    M, N = size(Z)[1:2]
-    C = similar(Z, M*(K-1), N*(K-1))
-    fill!(C,0)
+    C = zeros(T, M*(K-1), N*(K-1))   
 
-    @assert M == N
     Id = Matrix{T}(I, M, N)
     for m in 2:K-1
         n = m-1
         C[(m-1)*M+1:m*M, (n-1)*N+1:n*N] = Id
     end
 
-    W = -inv(Z[:,:,1])
+    W = -inv(timeslice(Z,1))
     for n in 1:K-1
         m = 1
-        C[(m-1)*M+1:m*M, (n-1)*N+1:n*N] = W*Z[:,:,n+1]
+        C[(m-1)*M+1:m*M, (n-1)*N+1:n*N] = W*timeslice(Z,n+1)
     end
-
+    
     return C
 end
 
-
-function materializeconvop(Z)
-    M = size(Z,1)
+```
+    Build the companion matrix for Z that has a tail. The `ranktail` specifies the dimension of the range of tail. SVD does not require the tail to be symmetric
+```
+function companion_with_tail(Z; ranktail=size(Z,1))
+    M, N = size(Z)[1:2]
     T = eltype(Z)
-
-    if hastail(Z)
-        kmax = ConvolutionOperators.tailindex(Z)
-        Q = zeros(T, 2M, 2M, kmax+1)
-        for k in 1:kmax
-            Q[1:M,1:M,k] .= timeslice(Z,k)
+    @assert M == N
+    @assert ranktail <= M
+    Mt = min(ranktail, M)
+
+    K = ConvolutionOperators.tailindex(Z) 
+    @assert K > 1
+
+    C = zeros(T, M*(K-2)+Mt, N*(K-2)+Mt)
+
+    Id = Matrix{T}(I, M, N)
+    Ztail = timeslice(Z, K)
+    if Mt < M
+        SVD = svd(Ztail)
+        S = zeros(T, M, Mt)
+        for i in 1:Mt
+            S[i, i] = SVD.S[i]
         end
-        Id = Matrix{T}(LinearAlgebra.I,M,M)
-        Q[M+1:2M,1:M,1] .= -Id
-        Q[M+1:2M,M+1:2M,1] .= Id
-        Q[M+1:2M,M+1:2M,2] .= -Id
-        Q[1:M,M+1:2M,kmax+1] .= timeslice(Z,kmax+1)
-        # return eigvals(companion(Q))
+        LD = SVD.U*S
+        RD = SVD.Vt[1:Mt, :]
     else
-        Q = zeros(T, M, M, size(Z,3))
-        for k in 1:size(Z,3)
-            Q[:,:,k] = timeslice(Z,k)
-        end
+        LD = Ztail
+        RD = Id
+    end
+
+    for m in 2:K-2
+        n = m-1
+        C[(m-1)*M+1:m*M, (n-1)*N+1:n*N] = Id
     end
-    return Q
+    C[(K-2)*M+1:(K-2)*M+Mt, (K-3)*N+1:(K-2)*N] = RD
+
+    W = -inv(timeslice(Z,1))
+    for n in 1:K-2
+        m = 1
+        C[(m-1)*M+1:m*M, (n-1)*N+1:n*N] = W*timeslice(Z,n+1)
+    end
+    C[1:M, (K-2)*N+1:(K-2)*N+Mt] = W*LD
+
+    Idt = Matrix{T}(I, Mt, Mt)
+    C[(K-2)*M+1:(K-2)*M+Mt, (K-2)*N+1:(K-2)*N+Mt] = Idt      
+
+    return C
 end
 
-function polyvals(Z)
-    Q = materializeconvop(Z)
-    C = companion(Q)
+function polyvals(Z; ranktail=size(Z,1))
+    if hastail(Z)
+        C = companion_with_tail(Z; ranktail)
+    else
+        C = companion_no_tail(Z)
+    end
     @show size(C)
     return eigvals(C)
 end
 
-function polyeig(Z)
-    Q = materializeconvop(Z)
-    C = companion(Q)
+function polyeig(Z; ranktail=size(Z,1))
+    if hastail(Z)
+        C = companion_with_tail(Z; ranktail)
+    else
+        C = companion_no_tail(Z)
+    end
     @show size(C)
     w, V = eigen(C)
     M = size(Z,1)

src/leftrightmul.jl

@@ -32,6 +32,8 @@ function Base.:*(A::LinearMap, C::AbstractConvOp) LeftRightMulCVO(C,A,I) end
 function Base.:*(C::AbstractConvOp, B::AbstractArray) LeftRightMulCVO(C,I,B) end
 function Base.:*(A::AbstractArray, C::AbstractConvOp) LeftRightMulCVO(C,A,I) end
 
+function Base.eltype(x::LeftRightMulCVO) eltype(x.convop) end
+
 function convolve!(y, Z::LeftRightMulCVO, x, X, j, k_start=1, k_stop=size(Z,3))
 
     CVO = Z.convop
@@ -63,6 +65,10 @@ function timeslice!(Y, Z::LeftRightMulCVO, k)
     Y .= Matrix(A*C*B)
 end
 
+function tailindex(Z::LeftRightMulCVO)
+    tailindex(Z.convop)
+end
+
 
 function hastail(Z::LeftRightMulCVO)
     hastail(Z.convop)

src/matrixconvop.jl

@@ -0,0 +1,86 @@
+struct MatrixConvOp <: AbstractConvOp
+    matconvop::Matrix{AbstractConvOp}
+end
+
+function Base.eltype(Z::MatrixConvOp)
+    mat = Z.matconvop
+    type = eltype(mat[1, 1])
+    for i in 1:size(mat, 1)
+        for j in 1:size(mat, 2)
+            type = promote_type(type, eltype(mat[i, j]))
+        end
+    end
+    return type
+end
+
+function Base.size(Z::MatrixConvOp)
+    mat = Z.matconvop
+    M, N, K = 0, 0, 0
+    for i in 1:size(mat, 1)
+        M += size(mat[i, 1], 1)
+        for j in 1:size(mat, 2)
+            if i == 1 
+                N += size(mat[1, j], 2) 
+            end
+            K = max(K, size(mat[i, j], 3))
+        end
+    end
+    return (M, N, K)
+end
+
+function convolve!(y, Z::MatrixConvOp, x, X, j, k_start=1, k_stop=size(Z,3))
+    mat = Z.matconvop
+    M0 = 0
+    for i1 in 1:size(mat, 1)
+        N0 = 0
+        for i2 in 1:size(mat, 2)
+            M, N = size(mat[i1, i2])[1:2]
+            xi2 = x[N0+1:N0+N, :]
+            Xi2 = X[N0+1:N0+N, :]
+            Y = zeros(M)
+            convolve!(Y, mat[i1, i2], xi2, Xi2, j, k_start, k_stop)
+            y[M0+1:M0+M] .+= Y
+            N0 += N
+        end 
+        M0 += size(mat[i1, 1], 1)
+    end
+    return y
+end
+
+function timeslice!(Y, Z::MatrixConvOp, k)
+    fill!(Y, 0)
+    mat = Z.matconvop
+    M0 = 0
+    for i in 1:size(mat, 1)
+        N0 = 0
+        for j in 1:size(mat, 2)
+            slice = timeslice(mat[i, j], k)
+            M, N = size(slice)
+            Y[M0+1:M0+M, N0+1:N0+N] .= slice
+            N0 += N
+        end
+        M0 += size(mat[i, 1], 1)
+    end
+    return Y
+end
+
+function tailindex(Z::MatrixConvOp) 
+    mat = Z.matconvop
+    ti = 0
+    for i in 1:size(mat, 1)
+        for j in 1:size(mat, 2)
+            ti = max(ti, tailindex(mat[i, j]))
+        end
+    end
+    return ti
+end
+
+function hastail(Z::MatrixConvOp)
+    mat = Z.matconvop
+    for i in 1:size(mat, 1)
+        for j in 1:size(mat, 2)
+            if hastail(mat[i, j]) return true end
+        end
+    end
+    false
+end

src/truncatedconvop.jl

@@ -12,9 +12,9 @@ function Base.eltype(x::TruncatedConvOp)
 end
 
 function Base.size(x::TruncatedConvOp)
-    size(x.convop)
-    # ln = minimum(size(x.convop)[3], x.kmax)
-    # return (size(x.convop)[1:2]..., ln)
+    # size(x.convop)
+    ln = min(size(x.convop)[3], x.kmax)
+    return (size(x.convop)[1:2]..., ln)
 end
 
 function convolve!(y, Z::TruncatedConvOp, x, X, j, k_start=1, k_stop=size(Z,3))
@@ -27,6 +27,8 @@ function timeslice!(Y, Z::TruncatedConvOp, k)
     timeslice!(Y, Z.convop, k)
 end
 
+tailindex(Z::TruncatedConvOp) = tailindex(Z.convop)
+
 function hastail(Z::TruncatedConvOp)
     false
 end