Multiscale Graph Signal Transforms on 1D Path

Let us use the unweighted 1D path with a synthetic signal as a simple example to demonstrate the usage of the MultiscaleGraphSignalTransforms.jl.

Set up

We first construct the GraphSig and GraphPart objects of the primal graph $G = P_{64}$.

using MultiscaleGraphSignalTransforms, Plots, LinearAlgebra
import WaveletsExt: wiggle

# construct P64
N = 64
G = gpath(N)

# compute graph Laplacian eigenvectors
W = G.W
L = diagm(sum(W; dims = 1)[:]) - W  # unnormalized graph Laplacian
𝛌, 𝚽 = eigen(L)
𝚽 = 𝚽 .* sign.(𝚽[1,:])'

# perform recursive bipartitioning of G by the Fiedler vectors of Lrw
GP = partition_tree_fiedler(G; swapRegion = false)

# use Chebyshev polynomial T₅(x) (x ∈ [0, 1]) as an example signal
G.f = reshape([16 * x^5 - 20 * x^3 + 5 * x for x in LinRange(0, 1, N)], (N, 1))
plot(G.f; c = :black, lw = 2, legend = false, grid = false, size = (815, 300))

Graph Signal Processing via HGLET/LP-HGLET

## analyze the signal via HGLET
dmatrixH, dmatrixHrw, dmatrixHsym = HGLET_Analysis_All(G, GP)
dvec_hglet, BS_hglet, trans_hglet = HGLET_GHWT_BestBasis(GP, dmatrixH = dmatrixH,
    dmatrixHrw = dmatrixHrw, dmatrixHsym = dmatrixHsym, costfun = 1)

## LP-HGLET
dmatrixsH, dmatrixsHsym = LPHGLET_Analysis_All(G, GP; ϵ = 0.3)
dvec_lphglet, BS_lphglet, trans_lphglet = HGLET_GHWT_BestBasis(GP, dmatrixH = dmatrixsH,
    dmatrixHsym = dmatrixsHsym, costfun = 1)

# find the top 10 HGLET basis vectors
important_idx = sortperm(dvec_hglet[:].^2; rev = true)
hglet_top10 = zeros(N, 10)
for i in 1:10
    w, _ = HGLET_Synthesis(reshape(spike(important_idx[i], N), (N, 1)),
        GP, BS_hglet, G, method = :L)
    hglet_top10[:, i] = w[:]
end
wiggle(hglet_top10; sc = 0.45)
p1 = title!("Top 10 HGLET basis vectors")

# find the top 10 LP-HGLET basis vectors
important_idx = sortperm(dvec_lphglet[:].^2; rev = true)
lphglet_top10 = zeros(N, 10)
for i in 1:10
    w, _ = LPHGLET_Synthesis(reshape(spike(important_idx[i], N), (N, 1)),
        GP, BS_lphglet, G; method = :L, ϵ = 0.3)
    lphglet_top10[:, i] = w[:]
end
wiggle(lphglet_top10; sc = 0.45)
p2 = title!("Top 10 LP-HGLET basis vectors")

plot(p1, p2, layout = Plots.grid(2, 1), size = (815, 600))
xticks!([1; 8:8:64], vcat(string("1"), [string(k) for k in 8:8:64])) # hide
yticks!([0; 1:10], vcat(string(""), [string(k) for k in 1:10])) # hide
plot!(left_margin = 5mm) # hide

Graph Signal Processing via GHWT, eGHWT, etc.

## analyze the signal via GHWT
dmatrix = ghwt_analysis!(G, GP = GP)

## Haar
BS_haar = bs_haar(GP)
dvec_haar = dmatrix2dvec(dmatrix, GP, BS_haar)

## Walsh
BS_walsh = bs_walsh(GP)
dvec_walsh = dmatrix2dvec(dmatrix, GP, BS_walsh)

## GHWT_c2f
dvec_c2f, BS_c2f = ghwt_c2f_bestbasis(dmatrix, GP)

## GHWT_f2c
dvec_f2c, BS_f2c = ghwt_f2c_bestbasis(dmatrix, GP)

## eGHWT
dvec_eghwt, BS_eghwt = ghwt_tf_bestbasis(dmatrix, GP)

We then find the top 10 basis vectors in each case.

## Haar
important_idx = sortperm(dvec_haar[:].^2; rev = true)
haar_top10 = zeros(N, 10)
for i in 1:10
    w = ghwt_synthesis(reshape(spike(important_idx[i], N), (N, 1)), GP, BS_haar)
    haar_top10[:, i] = w[:]
end
wiggle(haar_top10; sc = 0.45)
p1 = title!("Top 10 Haar basis vectors")

## Walsh
important_idx = sortperm(dvec_walsh[:].^2; rev = true)
walsh_top10 = zeros(N, 10)
for i in 1:10
    w = ghwt_synthesis(reshape(spike(important_idx[i], N), (N, 1)), GP, BS_walsh)
    walsh_top10[:, i] = w[:]
end
wiggle(walsh_top10; sc = 0.45)
p2 = title!("Top 10 Walsh basis vectors")

## GHWT_c2f
important_idx = sortperm(dvec_c2f[:].^2; rev = true)
ghwt_c2f_top10 = zeros(N, 10)
for i in 1:10
    w = ghwt_synthesis(reshape(spike(important_idx[i], N), (N, 1)), GP, BS_c2f)
    ghwt_c2f_top10[:, i] = w[:]
end
wiggle(ghwt_c2f_top10; sc = 0.45)
p3 = title!("Top 10 GHWT c2f best basis vectors")

## GHWT_f2c
important_idx = sortperm(dvec_f2c[:].^2; rev = true)
ghwt_f2c_top10 = zeros(N, 10)
for i in 1:10
    w = ghwt_synthesis(reshape(spike(important_idx[i], N), (N, 1)), GP, BS_f2c)
    ghwt_f2c_top10[:, i] = w[:]
end
wiggle(ghwt_f2c_top10; sc = 0.45)
p4 = title!("Top 10 GHWT f2c best basis vectors")

## eGHWT
important_idx = sortperm(dvec_eghwt[:].^2; rev = true)
eghwt_top10 = zeros(N, 10)
for i in 1:10
    w = ghwt_synthesis(reshape(spike(important_idx[i], N), (N, 1)), GP, BS_eghwt)
    eghwt_top10[:, i] = w[:]
end
wiggle(eghwt_top10; sc = 0.45)
p5 = title!("Top 10 eGHWT best basis vectors")

# display the top 10 basis vectors
plot(p1, p2, p3, p4, p5, layout = Plots.grid(5, 1), size = (815, 1500))

Graph Signal Processing via the NGWP dictionaries

To perform the NGWP transforms, we set up the dual graph $G^{\star}$ (which is also $P_{64}$).

# build the dual graph object
Gstar = GraphSig(W)

# perform recursive bipartitioning of Gstar by the Fiedler vectors of Lrw
GstarP = partition_tree_fiedler(Gstar; swapRegion = false)

# perform the pair-clustering algorithm to recursively bipartition G
GP_pc = pairclustering(𝚽, GstarP)  # for PC-NGWP

Now, let us construct the three NGWP dictionaries (i.e., the VM-NGWP, the PC-NGWP, and the LP-NGWP) and use them to analyze the signal, respectively.

VM_NGWP = vm_ngwp(𝚽, GstarP)
PC_NGWP = pc_ngwp(𝚽, GstarP, GP_pc)
LP_NGWP = lp_ngwp(𝚽, W, GstarP; ϵ = 0.3)  # relative action region bandwidth ϵ

# NGWP analysis, i.e., get the expansion coefficient matrix and apply the best
# basis algorithm.
dmatrix_VM = ngwp_analysis(G, VM_NGWP)
dvec_vm_ngwp, BS_vm_ngwp = ngwp_bestbasis(dmatrix_VM, GstarP)
dmatrix_PC = ngwp_analysis(G, PC_NGWP)
dvec_pc_ngwp, BS_pc_ngwp = ngwp_bestbasis(dmatrix_PC, GstarP)
dmatrix_LP = ngwp_analysis(G, LP_NGWP)
dvec_lp_ngwp, BS_lp_ngwp = ngwp_bestbasis(dmatrix_LP, GstarP)

Then, the top 10 NGWP basis vectors selected from each dictionary can be displayed as follows.

important_idx = sortperm(dvec_vm_ngwp[:].^2; rev = true)
wav_vm_top10 = zeros(N, 10)
for i in 1:10
    dr, dc = BS_vm_ngwp.levlist[important_idx[i]]
    wav_vm_top10[:, i] = VM_NGWP[dr, dc, :]
end
wiggle(wav_vm_top10; sc = 0.45)
p1 = title!("Top 10 VM-NGWP basis vectors")

important_idx = sortperm(dvec_pc_ngwp[:].^2; rev = true)
wav_pc_top10 = zeros(N, 10)
for i in 1:10
    dr, dc = BS_pc_ngwp.levlist[important_idx[i]]
    wav_pc_top10[:, i] = PC_NGWP[dr, dc, :]
end
wiggle(wav_pc_top10; sc = 0.45)
p2 = title!("Top 10 PC-NGWP basis vectors")

important_idx = sortperm(dvec_lp_ngwp[:].^2; rev = true)
wav_lp_top10 = zeros(N, 10)
for i in 1:10
    dr, dc = BS_lp_ngwp.levlist[important_idx[i]]
    wav_lp_top10[:, i] = LP_NGWP[dr, dc, :]
end
wiggle(wav_lp_top10; sc = 0.45)
p3 = title!("Top 10 LP-NGWP basis vectors")

plot(p1, p2, p3, layout = Plots.grid(3, 1), size = (815, 900))