diff --git a/src/ConvolutionOperators.jl b/src/ConvolutionOperators.jl index 6d0f5da..867510f 100644 --- a/src/ConvolutionOperators.jl +++ b/src/ConvolutionOperators.jl @@ -25,6 +25,7 @@ include("linearcombinations.jl") include("leftrightmul.jl") include("liftedconvop.jl") include("truncatedconvop.jl") +include("matrixconvop.jl") end diff --git a/src/companion.jl b/src/companion.jl index 8b9a89c..f5c3277 100644 --- a/src/companion.jl +++ b/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) + @assert M == N + K = size(Z,3) @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) + 1 + @assert K > 2 + + 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) diff --git a/src/leftrightmul.jl b/src/leftrightmul.jl index 2b7a78d..a237ea9 100644 --- a/src/leftrightmul.jl +++ b/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) diff --git a/src/matrixconvop.jl b/src/matrixconvop.jl new file mode 100644 index 0000000..660e60f --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/truncatedconvop.jl b/src/truncatedconvop.jl index 6541032..c86af74 100644 --- a/src/truncatedconvop.jl +++ b/src/truncatedconvop.jl @@ -12,12 +12,15 @@ 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)) + if k_start > Z.kmax + return y + end k_stop = min(k_stop, Z.kmax) convolve!(y, Z.convop, x, X, j, k_start, k_stop) end @@ -27,6 +30,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 \ No newline at end of file