GitHub

Kronecker sum of circulant

This example illustrates efficient computation of the inverse of a Kronecker sum of circulant matrices using the Julia language.

Related to Problem 8.7 in the textbook.

This page comes from a single Julia file: kron-sum-inv.jl.

You can access the source code for such Julia documentation using the 'Edit on GitHub' link in the top right. You can view the corresponding notebook in nbviewer here: kron-sum-inv.ipynb, or open it in binder here: kron-sum-inv.ipynb.

Setup

Add the Julia packages used in this demo. Change false to true in the following code block if you are using any of the following packages for the first time.

if false
    import Pkg
    Pkg.add([
        "FFTW"
        "InteractiveUtils"
        "LinearAlgebra"
        "MIRTjim"
        "Plots"
        "Random"
    ])
end

Tell Julia to use the following packages. Run Pkg.add() in the preceding code block first, if needed.

using FFTW: fft, ifft
using InteractiveUtils: versioninfo
using LinearAlgebra: I
using MIRTjim: jim, prompt
using Plots: savefig
using Random: seed!
seed!(0)
Random.TaskLocalRNG()

The following line helps when running this jl-file as a script; this way it prompts user to hit a key after each image is displayed.

isinteractive() && prompt(:prompt);

Circulant matrix with given first column

function circm(x)
    N = length(x)
    return hcat([circshift(x, k-1) for k in 1:N]...)
end
circm (generic function with 1 method)

Numerical test

First perform a numerical test to verify the formulas.

Test data

M, N = 64, 32
b = randn(N)
c = randn(M)
X = randn(M, N);

Matrix-inverse solution

B = circm(b) # N
C = circm(c) # M
A = kron(B, I(M)) + kron(I(N), C)
y = inv(A) * vec(X);

Denominator

P = fft(c) .+ transpose(fft(b));

DFT solution

Fn = fft(I(N), 1)
Fm = fft(I(M), 1)
Y1 = (Fm' * ((Fm * X * transpose(Fn)) ./ P) * conj(Fn)) / (M * N)
@assert Y1 ≈ real(Y1) # should be real-valued
Y1 = real(Y1); # so discard the imaginary part

FFT solution

Y2 = ifft(fft(X) ./ P)
@assert Y2 ≈ real(Y2)
Y2 = real(Y2);

Verify

@assert y ≈ vec(Y1)
@assert y ≈ vec(Y2)

Image data

Now illustrate visually using additively separable blur with individually invertible blur kernels.

X = ones(60, 64); X[20:50,10:50] .= 2
M, N = size(X)
b = zeros(N); b[1 .+ mod.(-3:3,N)] = (4 .- abs.(-3:3)); b[1] +=  1
b /= sum(b)
c = zeros(M); c[1 .+ mod.(-4:4,M)] = (5 .- abs.(-4:4)); c[1] +=  1
c /= sum(c)
B = circm(b)
C = circm(c)
pb = jim(B, "circulant B", size=(300,300))
pc = jim(C, "circulant C", size=(300,300))
p1 = jim(pb, pc; size=(600,300))
Example block output
prompt()


p2 = jim(X; title="X original", size=(300,300))
Y = C * X + X * transpose(B)
p3 = jim(Y; title="Y blurred", size=(300,300))
p23 = jim(p2, p3; size=(600,300))
Example block output
prompt()

FFT solution

P = fft(c) .+ transpose(fft(b)) # denominator
Xhat = ifft(fft(Y) ./ P)
@assert Xhat ≈ real(Xhat)
Xhat = real(Xhat)
@assert Xhat ≈ X

Reproducibility

This page was generated with the following version of Julia:

using InteractiveUtils: versioninfo
io = IOBuffer(); versioninfo(io); split(String(take!(io)), '\n')
11-element Vector{SubString{String}}:
 "Julia Version 1.11.1"
 "Commit 8f5b7ca12ad (2024-10-16 10:53 UTC)"
 "Build Info:"
 "  Official https://julialang.org/ release"
 "Platform Info:"
 "  OS: Linux (x86_64-linux-gnu)"
 "  CPU: 4 × AMD EPYC 7763 64-Core Processor"
 "  WORD_SIZE: 64"
 "  LLVM: libLLVM-16.0.6 (ORCJIT, znver3)"
 "Threads: 1 default, 0 interactive, 1 GC (on 4 virtual cores)"
 ""

And with the following package versions

import Pkg; Pkg.status()
Status `~/work/book-la-demo/book-la-demo/docs/Project.toml`
  [6e4b80f9] BenchmarkTools v1.5.0
  [aaaa29a8] Clustering v0.15.7
  [35d6a980] ColorSchemes v3.27.1
⌅ [3da002f7] ColorTypes v0.11.5
⌃ [c3611d14] ColorVectorSpace v0.10.0
  [717857b8] DSP v0.7.10
  [72c85766] Demos v0.1.0 `~/work/book-la-demo/book-la-demo`
  [e30172f5] Documenter v1.7.0
  [4f61f5a4] FFTViews v0.3.2
  [7a1cc6ca] FFTW v1.8.0
  [587475ba] Flux v0.14.25
  [a09fc81d] ImageCore v0.10.4
  [71a99df6] ImagePhantoms v0.8.1
  [b964fa9f] LaTeXStrings v1.4.0
  [7031d0ef] LazyGrids v1.0.0
  [599c1a8e] LinearMapsAA v0.12.0
  [98b081ad] Literate v2.20.1
  [7035ae7a] MIRT v0.18.2
  [170b2178] MIRTjim v0.25.0
  [eb30cadb] MLDatasets v0.7.18
  [efe261a4] NFFT v0.13.5
  [6ef6ca0d] NMF v1.0.3
  [15e1cf62] NPZ v0.4.3
  [0b1bfda6] OneHotArrays v0.2.5
  [429524aa] Optim v1.10.0
  [91a5bcdd] Plots v1.40.9
  [f27b6e38] Polynomials v4.0.11
  [2913bbd2] StatsBase v0.34.3
  [d6d074c3] VideoIO v1.1.0
  [b77e0a4c] InteractiveUtils v1.11.0
  [37e2e46d] LinearAlgebra v1.11.0
  [44cfe95a] Pkg v1.11.0
  [9a3f8284] Random v1.11.0
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`

This page was generated using Literate.jl.