From 85ed6cba3190e988b0a57db07b533fc2cf96ab86 Mon Sep 17 00:00:00 2001 From: lvchien Date: Mon, 12 Dec 2022 10:56:07 +0100 Subject: [PATCH 01/10] truncated convops --- src/truncatedconvop.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/truncatedconvop.jl b/src/truncatedconvop.jl index 6541032..baf1824 100644 --- a/src/truncatedconvop.jl +++ b/src/truncatedconvop.jl @@ -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 \ No newline at end of file From f39dc7178ccf0481e0791be457ab006c4b25e2c4 Mon Sep 17 00:00:00 2001 From: lvchien Date: Mon, 23 Oct 2023 11:32:22 +0200 Subject: [PATCH 02/10] Add eltype of LeftRightMul Convops --- src/leftrightmul.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/leftrightmul.jl b/src/leftrightmul.jl index 2b7a78d..24d59d5 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 From a185adbe36601c36b79f5aad81cb5ef9b26fb816 Mon Sep 17 00:00:00 2001 From: lvchien Date: Mon, 16 Sep 2024 16:34:35 +0200 Subject: [PATCH 03/10] Reduce the size of companion matrix --- src/companion.jl | 56 ++++++++++++++++-------------------------------- 1 file changed, 19 insertions(+), 37 deletions(-) diff --git a/src/companion.jl b/src/companion.jl index 8b9a89c..9884bf3 100644 --- a/src/companion.jl +++ b/src/companion.jl @@ -1,64 +1,46 @@ -function companion(Z) +# TODO: remove the nullspace of the tail +function companion(Z) + M, N = size(Z)[1:2] T = eltype(Z) - K = size(Z,3) - @assert K > 1 + @assert M == N - M, N = size(Z)[1:2] - C = similar(Z, M*(K-1), N*(K-1)) - fill!(C,0) + if hastail(Z) + K = ConvolutionOperators.tailindex(Z) + else + K = size(Z,3) - 1 + end + @assert K > 1 - @assert M == N + C = zeros(T, M*(K-1), N*(K-1)) 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) - 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) - 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)) - else - Q = zeros(T, M, M, size(Z,3)) - for k in 1:size(Z,3) - Q[:,:,k] = timeslice(Z,k) - end + m, n = K-1, K-1 + C[(m-1)*M+1:m*M, (n-1)*N+1:n*N] = Id end - return Q + + return C end function polyvals(Z) - Q = materializeconvop(Z) - C = companion(Q) + C = companion(Z) @show size(C) return eigvals(C) end function polyeig(Z) - Q = materializeconvop(Z) C = companion(Q) @show size(C) w, V = eigen(C) From 35e4db15b4cb611d3a0b26a6cb40c0c566ba9c00 Mon Sep 17 00:00:00 2001 From: lvchien Date: Sat, 28 Sep 2024 23:48:02 +0200 Subject: [PATCH 04/10] Remove the nullspace of tail --- src/companion.jl | 83 ++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 67 insertions(+), 16 deletions(-) diff --git a/src/companion.jl b/src/companion.jl index 9884bf3..9c9fa88 100644 --- a/src/companion.jl +++ b/src/companion.jl @@ -1,20 +1,19 @@ -# TODO: remove the nullspace of the tail +using LinearAlgebra -function companion(Z) +``` + 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 - if hastail(Z) - K = ConvolutionOperators.tailindex(Z) - else - K = size(Z,3) - 1 - end + K = size(Z,3) - 1 @assert K > 1 - C = zeros(T, M*(K-1), N*(K-1)) - Id = Matrix{T}(I, M, N) + C = zeros(T, M*(K-1), N*(K-1)) + 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 @@ -26,22 +25,74 @@ function companion(Z) C[(m-1)*M+1:m*M, (n-1)*N+1:n*N] = W*timeslice(Z,n+1) end - if hastail(Z) - m, n = K-1, K-1 - C[(m-1)*M+1:m*M, (n-1)*N+1:n*N] = Id + return C +end + +``` + 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) + M, N = size(Z)[1:2] + T = eltype(Z) + @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 + LD = SVD.U*S + RD = SVD.Vt[1:Mt, :] + else + 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 + 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) - C = companion(Z) +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) - 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) From 6751c8f112f3f3645848d4ab362abbe53650c285 Mon Sep 17 00:00:00 2001 From: lvchien Date: Wed, 16 Oct 2024 16:29:56 +0200 Subject: [PATCH 05/10] Add matrix of convolution operators --- src/ConvolutionOperators.jl | 1 + src/matrixconvop.jl | 76 +++++++++++++++++++++++++++++++++++++ 2 files changed, 77 insertions(+) create mode 100644 src/matrixconvop.jl 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/matrixconvop.jl b/src/matrixconvop.jl new file mode 100644 index 0000000..a76af54 --- /dev/null +++ b/src/matrixconvop.jl @@ -0,0 +1,76 @@ +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 + blocksize = size(mat[1, 1]) + sz3 = 0 + for i in 1:size(mat, 1) + for j in 1:size(mat, 2) + sz3 = max(sz3, size(mat[i, j], 3)) + end + end + return (blocksize[1]*size(mat, 1), blocksize[2]*size(mat, 2), sz3) +end + +function convolve!(y, Z::MatrixConvOp, x, X, j, k_start=1, k_stop=size(Z,3)) + mat = Z.matconvop + blocksize = size(mat[1, 1]) + + for i1 in 1:size(mat, 1) + for i2 in 1:size(mat, 2) + xi2 = x[(i2-1)*blocksize[2]+1:i2*blocksize[2], :] + Xi2 = X[(i2-1)*blocksize[2]+1:i2*blocksize[2], :] + Y = zeros(size(xi2)[1]) + convolve!(Y, mat[i1, i2], xi2, Xi2, j, k_start, k_stop) + y[(i1-1)*blocksize[1]+1:i1*blocksize[1]] .+= Y + end + end + return y +end + +function timeslice!(Y, Z::MatrixConvOp, k) + fill!(Y, 0) + mat = Z.matconvop + slice = timeslice.(mat, k) + blocksize = size(mat[1, 1]) + for i in 1:size(mat, 1) + for j in 1:size(mat, 2) + Y[(i-1)*blocksize[1]+1:i*blocksize[1], (j-1)*blocksize[2]+1:j*blocksize[2]] .= slice[i, j] + end + 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 From e76fe46a67b7e473f95b285afa00437f14118bfa Mon Sep 17 00:00:00 2001 From: lvchien Date: Fri, 3 Jan 2025 16:27:06 +0700 Subject: [PATCH 06/10] matrixconvop for rectangular matrices --- src/matrixconvop.jl | 36 +++++++++++++++++++++++------------- 1 file changed, 23 insertions(+), 13 deletions(-) diff --git a/src/matrixconvop.jl b/src/matrixconvop.jl index a76af54..660e60f 100644 --- a/src/matrixconvop.jl +++ b/src/matrixconvop.jl @@ -15,28 +15,34 @@ end function Base.size(Z::MatrixConvOp) mat = Z.matconvop - blocksize = size(mat[1, 1]) - sz3 = 0 + 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) - sz3 = max(sz3, size(mat[i, j], 3)) + if i == 1 + N += size(mat[1, j], 2) + end + K = max(K, size(mat[i, j], 3)) end end - return (blocksize[1]*size(mat, 1), blocksize[2]*size(mat, 2), sz3) + return (M, N, K) end function convolve!(y, Z::MatrixConvOp, x, X, j, k_start=1, k_stop=size(Z,3)) mat = Z.matconvop - blocksize = size(mat[1, 1]) - + M0 = 0 for i1 in 1:size(mat, 1) + N0 = 0 for i2 in 1:size(mat, 2) - xi2 = x[(i2-1)*blocksize[2]+1:i2*blocksize[2], :] - Xi2 = X[(i2-1)*blocksize[2]+1:i2*blocksize[2], :] - Y = zeros(size(xi2)[1]) + 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[(i1-1)*blocksize[1]+1:i1*blocksize[1]] .+= Y + y[M0+1:M0+M] .+= Y + N0 += N end + M0 += size(mat[i1, 1], 1) end return y end @@ -44,12 +50,16 @@ end function timeslice!(Y, Z::MatrixConvOp, k) fill!(Y, 0) mat = Z.matconvop - slice = timeslice.(mat, k) - blocksize = size(mat[1, 1]) + M0 = 0 for i in 1:size(mat, 1) + N0 = 0 for j in 1:size(mat, 2) - Y[(i-1)*blocksize[1]+1:i*blocksize[1], (j-1)*blocksize[2]+1:j*blocksize[2]] .= slice[i, j] + 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 From a15cdc800359e5970c625d451b969f7b857e1fce Mon Sep 17 00:00:00 2001 From: lvchien Date: Thu, 20 Feb 2025 13:56:57 +0100 Subject: [PATCH 07/10] update size of truncated convop --- src/truncatedconvop.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/truncatedconvop.jl b/src/truncatedconvop.jl index baf1824..a3d4345 100644 --- a/src/truncatedconvop.jl +++ b/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)) From 0f28e9199bcbb2f40377ed6721c401a8b7238203 Mon Sep 17 00:00:00 2001 From: lvchien Date: Wed, 4 Jun 2025 15:06:21 +0200 Subject: [PATCH 08/10] default full rank tail --- src/companion.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/companion.jl b/src/companion.jl index 9c9fa88..164a09a 100644 --- a/src/companion.jl +++ b/src/companion.jl @@ -31,7 +31,7 @@ end ``` 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) +function companion_with_tail(Z; ranktail=size(Z,1)) M, N = size(Z)[1:2] T = eltype(Z) @assert M == N From 90375dc666702aeb49e8ab7cad20ff418642e7b6 Mon Sep 17 00:00:00 2001 From: lvchien Date: Wed, 3 Dec 2025 14:37:58 +0100 Subject: [PATCH 09/10] tailindex for LeftRightMulConvOp --- src/leftrightmul.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/leftrightmul.jl b/src/leftrightmul.jl index 24d59d5..a237ea9 100644 --- a/src/leftrightmul.jl +++ b/src/leftrightmul.jl @@ -65,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) From 1ce8b863673757ef8ec288b22112069a1f2dec29 Mon Sep 17 00:00:00 2001 From: lvchien Date: Mon, 2 Feb 2026 11:13:47 +0100 Subject: [PATCH 10/10] fix issue with tail index --- src/companion.jl | 6 +++--- src/truncatedconvop.jl | 3 +++ 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/companion.jl b/src/companion.jl index 164a09a..f5c3277 100644 --- a/src/companion.jl +++ b/src/companion.jl @@ -8,7 +8,7 @@ function companion_no_tail(Z) T = eltype(Z) @assert M == N - K = size(Z,3) - 1 + K = size(Z,3) @assert K > 1 C = zeros(T, M*(K-1), N*(K-1)) @@ -38,8 +38,8 @@ function companion_with_tail(Z; ranktail=size(Z,1)) @assert ranktail <= M Mt = min(ranktail, M) - K = ConvolutionOperators.tailindex(Z) - @assert K > 1 + K = ConvolutionOperators.tailindex(Z) + 1 + @assert K > 2 C = zeros(T, M*(K-2)+Mt, N*(K-2)+Mt) diff --git a/src/truncatedconvop.jl b/src/truncatedconvop.jl index a3d4345..c86af74 100644 --- a/src/truncatedconvop.jl +++ b/src/truncatedconvop.jl @@ -18,6 +18,9 @@ function Base.size(x::TruncatedConvOp) 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