Natural Graph Wavelet Dictionaries

natural_eigdist(𝚽, 𝛌, Q; α = 1.0, T = :Inf,
                input_format = :zero_measures, distance = :DAG,
                edge_weight = 1, edge_length = 1)

compute natural distances between graph Laplacian eigenvectors.

Input Arguments

• 𝚽::Matrix{Float64}: matrix of (weighted) graph Laplacian eigenvectors.
• 𝛌::Vector{Float64}: vector of eigenvalues.
• Q::Matrix{Float64}: unweighted incidence matrix of the graph.
• α::Float64: ROT parameter. (default: 1.0)
• T::Any: TSD parameter, i.e., the stopping time T in K_functional (default: :Inf)
• input_format::Symbol: options: :zero_measures, :pmf1 and :pmf2 (default: :zero_measures)
• distance::Symbol: options: :ROT, :HAD, :DAG and :TSD (default: :DAG)
• edg_length::Any: vector of edge lengths (default: 1 represents unweighted graphs)
• edge_weight::Any: the weights vector, which stores the affinity weight of each edge (default: 1 represents unweighted graphs).

Output Argument

• dis::Matrix{Float64}: the distance matrix, dis[i,j] = d(𝜙ᵢ₋₁, 𝜙ⱼ₋₁).

ROT_Distance(A, B, Q; edge_length = 1, α = 1.0)

computes the ROT distance matrix from A's column vectors to B's column vectors. If A, B are vector inputs, then it also returns the cost value and the optimal transport plan.

Input Argument

• A::Any: a vector or matrix whose columns are initial probability measures.
• B::Any: a vector or matrix whose columns are terminal probability measures.
• Q::SparseMatrixCSC{Int64,Int64}: the unweighted incidence matrix of the graph.
• edge_length::Any: the length vector (default: 1 represents unweighted graphs)
• α::Float64: default is 1.0. ROT parameter.
• retPlan::Bool: an indicator if return the optimal plan (default: false)

Output Argument

• dis::Matrix{Float64}: distance matrix, dis[i,j] = d_{ROT}(aᵢ, bⱼ; α).

eigROT_Distance(P, Q; edge_length = 1, α = 1.0)

computes the ROT distance matrix of P's column vectors on a graph.

Input Argument

• P::Matrix{Float64}: a matrix whose columns are vector measures with the same total mass.
• Q::SparseMatrixCSC{Int64,Int64}: the unweighted incidence matrix of the graph.
• edge_length::Any: the length vector (default: 1 represents unweighted graphs)
• α::Float64: default is 1.0. ROT parameter.

Output Argument

• dis::Matrix{Float64}: distance matrix, dis[i,j] = d_{ROT}(pᵢ, pⱼ; α).

eigsROT_Distance(P::Matrix{Float64}, W::SparseMatrixCSC{Float64, Int64}, X::Matrix{Float64}; α::Float64 = 1.0)

computes the sROT distance matrix of P's column vectors on a (unweighted) tree.

Input Argument

• P::Matrix{Float64}: a matrix whose columns are vector measures with the same total mass.
• W::SparseMatrixCSC{Float64, Int64}: the weight matrix of the tree.
• X::Matrix{Float64}: the node positions (i-th row represents node i's location)
• α::Float64: default is 1.0. ROT parameter.

Output Argument

• dist_sROT::Matrix{Float64}: distance matrix, distsROT[i,j] = d{sROT}(pᵢ, pⱼ; α).
• Ws::SparseMatrixCSC{Float64, Int64}: the weight matrix of the simplified tree.
• Xs::Matrix{Float64}: the node locations of the simplified tree
• 𝚯::Matrix{Float64}: the shortened pmfs from P.

find_subgraph_inds(Wc::SparseMatrixCSC{Float64, Int64})

find all subgraph indices of a tree. (subgraph includes: branches and junctions)

Input Argument

• Wc::SparseMatrixCSC{Float64, Int64}: the weight matrix of the tree chopped by the junctions.

Output Argument

• Ind::Matrix{Float64}: a matrix whose columns represent the subgraph node indices in binary form.

eigHAD_Affinity(𝚽, 𝛌; indexEigs = 1:size(𝚽,2))

EIGHAD_AFFINITY compute Hadamard (HAD) affinity between pairwise graph Laplacian eigenvectors.

Input Arguments

• 𝚽::Matrix{Float64}: matrix of graph Laplacian eigenvectors, 𝜙ⱼ₋₁ (j = 1,...,size(𝚽,1)).
• 𝛌::Array{Float64}: array of eigenvalues. (ascending order)
• indexEigs::Int: default is all eigenvectors, indices of eigenvectors considered.

Output Argument

• A::Matrix{Float64}: a numEigs x numEigs affinity matrix, A[i,j] = a_HAD(𝜙ᵢ₋₁, 𝜙ⱼ₋₁).

eigHAD_Distance(𝚽, 𝛌; indexEigs = 1:size(𝚽,2))

compute HAD "distance" (not really a distance) between pairwise graph Laplacian eigenvectors, i.e., dHAD(𝜙ᵢ₋₁, 𝜙ⱼ₋₁) = √(1 - aHAD(𝜙ᵢ₋₁, 𝜙ⱼ₋₁)²).

Input Arguments

• 𝚽::Matrix{Float64}: matrix of graph Laplacian eigenvectors, 𝜙ⱼ₋₁ (j = 1,...,size(𝚽,1)).
• 𝛌::Array{Float64}: array of eigenvalues. (ascending order)
• indexEigs::Int: default is all eigenvectors, indices of eigenvectors considered.

Output Argument

• dis::Matrix{Float64}: the HAD distance matrix, dis[i,j] = d_HAD(𝜙ᵢ₋₁, 𝜙ⱼ₋₁).

eigDAG_Distance(𝚽, Q, numEigs; edge_weights = 1)

compute DAG distances between pairwise graph Laplacian eigenvectors.

Input Arguments

• 𝚽::Matrix{Float64}: matrix of graph Laplacian eigenvectors, 𝜙ⱼ₋₁ (j = 1,...,size(𝚽,1)).
• Q::Matrix{Float64}: incidence matrix of the graph.
• numEigs::Int64: number of eigenvectors considered.
• edge_weight::Any: default value is 1, stands for unweighted graph (i.e., all edge weights equal to 1). For weighted graph, edge_weight is the weights vector, which stores the affinity weight of each edge.

Output Argument

• dis::Matrix{Float64}: a numEigs x numEigs distance matrix, dis[i,j] = d_DAG(𝜙ᵢ₋₁, 𝜙ⱼ₋₁).

K_functional(𝐩::Vector{Float64}, 𝐪::Vector{Float64}, 𝚽::Matrix{Float64},
             ∇𝚽::Matrix{Float64}, 𝛌::Vector{Float64}; length::Any = 1,
             T::Any = :Inf, dt::Float64 = 0.5/maximum(𝛌), tol::Float64 = 1e-5)

computes the K_functional between two vector meassures 𝐩 and 𝐪 on a graph.

Input Argument

• 𝐩::Vector{Float64}: the source vector measure.
• 𝐪::Vector{Float64}: the destination vector measure.
• 𝚽::Matrix{Float64}: matrix of the unweighted graph Laplacian eigenvectors.
• ∇𝚽::Matrix{Float64}: gradient of unweighted graph Laplacian eigenvectors.
• 𝛌::Vector{Float64}: vector of eigenvalues.
• length::Any: vector of edge lengths (default: 1 represents unweighted graphs)
• T::Any: the stopping time T in K_functional (default: :Inf)
• tol::Float64: tolerance for convergence (default: 1e-5)

Output Argument

• K::Float64: TSD distance d_{TSD}(p, q; T).
• E::Float64: an estimated upper bound on the absolute error. In general, E <= tol * norm(K).

eigTSD_Distance(P::Matrix{Float64}, 𝚽::Matrix{Float64}, 𝛌::Vector{Float64},
                Q::SparseMatrixCSC{Int64,Int64}; length::Any = 1,
                T::Any = :Inf, tol::Float64 = 1e-5)

computes the TSD distance matrix of P's column vectors on a graph.

Input Argument

• P::Matrix{Float64}: vector measures with the same total mass 0.
• 𝚽::Matrix{Float64}: matrix of the unweighted graph Laplacian eigenvectors.
• 𝛌::Vector{Float64}: vector of eigenvalues.
• Q::SparseMatrixCSC{Int64,Int64}: the unweighted incidence matrix.
• length::Any: vector of edge lengths (default: 1 represents unweighted graphs)
• T::Any: the stopping time T in K_functional (default: :Inf)
• tol::Float64: tolerance for integral convergence (default: 1e-5)

Output Argument

• dis::Matrix{Float64}: distance matrix, d_{TSD}(φᵢ, φⱼ; T).

dualgraph(dist::Matrix{Float64}; method::Symbol = :inverse, σ::Float64 = 1.0)

build the dual graph's weight matrix based on the given non-trivial eigenvector metric.

Input Arguments

• dist::Matrix{Float64}: eigenvector distance matrix
• method::Symbol: default is by taking inverse of the distance between eigenvectors. Ways to build the dual graph edge weights. Option: :inverse, :gaussian.
• σ::Float64: default is 1.0. Gaussian variance parameter.

Output Argument

• G_star::GraphSig: A GraphSig object containing the weight matrix of the dual graph.

varimax(A; gamma = 1.0, minit = 20, maxit = 1000, reltol = 1e-12)

VARIMAX perform varimax (or quartimax, equamax, parsimax) rotation to the column vectors of the input matrix.

Input Arguments

• A::Matrix{Float64}: input matrix, whose column vectors are to be rotated. d, m = size(A).
• gamma::Float64: default is 1. gamma = 0, 1, m/2, and d(m - 1)/(d + m - 2), corresponding to quartimax, varimax, equamax, and parsimax.
• minit::Int: default is 20. Minimum number of iterations, in case of the stopping criteria fails initially.
• maxit::Int: default is 1000. Maximum number of iterations.
• reltol::Float64: default is 1e-12. Relative tolerance for stopping criteria.

Output Argument

• B::Matrix{Float64}: output matrix, whose columns are already been rotated.

Implemented by Haotian Li, Aug. 20, 2019

vm_ngwp(𝚽::Matrix{Float64}, GP_star::GraphPart)

construct varimax NGWP and GP_star.tag in place.

Input Arguments

• 𝚽::Matrix{Float64}: graph Laplacian eigenvectors 𝚽
• GP_star::GraphPart: GraphPart object of the dual graph

Output Argument

• wavelet_packet::Array{Float64,3}: the varimax NGWP. The first index is for selecting wavelets at a fixed level; the second index is for selecting the level j; the third index is for selecting elements in the wavelet vector.

pairclustering(𝚽::Matrix{Float64}, GP_star::GraphPart)

construct the GraphPart object of the primal graph via pair-clustering algorithm.

Input Arguments

• 𝚽::Matrix{Float64}: graph Laplacian eigenvectors 𝚽
• GP_star::GraphPart: GraphPart object of the dual graph

Output Argument

• GP::GraphPart: GraphPart object of the primal graph via pair-clustering

mgslp(A::Matrix{Float64}; tol::Float64 = 1e-12, p::Float64 = 1.0)

Modified Gram-Schmidt Process Orthogonalization with ℓᵖ pivoting algorithm (MGSLp)

Input Arguments

• A::Matrix{Float64}: whose column vectors are to be orthogonalized.

Output Argument

• A::Matrix{Float64}: orthogonalization matrix of A's column vectors based on ℓᵖ pivoting.

pc_ngwp(𝚽::Matrix{Float64}, GP_star::GraphPart, GP::GraphPart)

construct pair-clustering NGWP and GP.tag in place.

Input Arguments

• 𝚽::Matrix{Float64}: graph Laplacian eigenvectors 𝚽
• GP_star::GraphPart: GraphPart object of the dual graph
• GP::GraphPart: GraphPart object of the primal graph

Output Argument

• wavelet_packet::Array{Float64,3}: the pair-clustering NGWP. The first index is for selecting wavelets at a fixed level; the second index is for selecting the level j; the third index is for selecting elements in the wavelet vector.

lp_ngwp(𝚽::Matrix{Float64}, W_dual::SparseMatrixCSC{Float64, Int64},
        GP_dual::GraphPart; ϵ::Float64 = 0.3)

construct the lapped NGWP and GP.tag in place.

Input Arguments

• 𝚽::Matrix{Float64}: graph Laplacian eigenvectors
• W_dual::SparseMatrixCSC{Float64, Int64}: weight matrix of the dual graph
• GP_dual::GraphPart: GraphPart object of the dual graph

Output Argument

• wavelet_packet::Array{Float64,3}: the lapped NGWP dictionary. The first index is for selecting wavelets at a fixed level; the second index is for selecting the level j; the third index is for selecting elements in the wavelet vector.

ngwp_analysis(G::GraphSig, wavelet_packet::Array{Float64,3})

For a GraphSig object G, generate the matrix of NGWP expansion coefficients.

Input Arguments

• G::GraphSig: an input GraphSig object
• wavelet_packet::Array{Float64,3}: the varimax wavelets packet.

Output Argument

• dmatrix::Array{Float64,3}: the expansion coefficients matrix.

(dvec_ngwp, BS_ngwp) = ngwp_bestbasis(dmatrix::Array{Float64,3}, GP_star::GraphPart;
                                      cfspec::Any = 1.0, flatten::Any = 1.0,
                                      j_start::Int = 1, j_end::Int = size(dmatrix, 2),
                                      useParent::Bool = true)

Select the best basis from the matrix of NGWP expansion coefficients.

Input Arguments

• dmatrix::Array{Float64,3}: the matrix of expansion coefficients
• GP_star::GraphPart: an input GraphPart object of the dual graph
• cfspec::Any: the specification of cost functional to be used (default = 1.0, i.e., 1-norm)
• flatten::Any: the method for flattening vector-valued data to scalar-valued data (default = 1.0, i.e, 1-norm)
• useParent::Bool: the flag to indicate if we update the selected best basis subspace to the parent when parent and child have the same cost (default = false)

Output Arguments

• dvec_ngwp::Matrix{Float64}: the vector of expansion coefficients corresponding to the NGWP best basis
• BS_ngwp::BasisSpec: a BasisSpec object which specifies the NGWP best basis

function NGWP_jkl(GP_star::GraphPart, drow::Int, dcol::Int)

Generate the (j,k,l) indices for the NGWP basis vector corresponding to the coefficient dmatrix(drow,dcol)

Input Arguments

• GP_star::GraphPart: a GraphPart object of the dual grpah
• drow::Int: the row of the expansion coefficient
• dcol::Int: the column of the expansion coefficient

Output Argument

• j: the level index of the expansion coefficient
• k: the subregion in dual graph's index of the expansion coefficient
• l: the tag of the expansion coefficient

frame_approx(f, U, V; num_kept = length(f))

approximate signal f by the frame U.

Input Arguments

• f::Vector{Float64}: input graph signal
• U::Matrix{Float64}: a frame operator (matrix or dictionary)
• V::Matrix{Float64}: the dual frame operator of U
• num_kept::Int64: number of kept coefficients (NCR)

Output Argument

• rel_error::Vector{Float64}: the relative errors
• f_approx::Vector{Float64}: the approximated signal

nat_spec_filter(l, D; σ = 0.25 * maximum(D), method = :regular, thres = 0.2)

assemble the natural spectral graph filter centered at the l-th eigenvector via the distance matrix D.

Input Arguments

• l::Int64: index of the centered eigenvector
• D::Matrix{Float64}: non-trivial distance matrix of the eigenvectors
• σ::Float64: Gaussian window width parameter (default: 0.25 * maximum(D))
• method::Symbol: :regular or :reduced (default: :regular)
• thres::Float64: cutoff threshold ∈ (0, 1).

Output Argument

• 𝛍::Vector{Float64}: the natural spectral graph filter

ngwf_all_vectors(D, 𝚽; σ = 0.2 * maximum(D))

assemble the whole NGWF dictionary.

Input Arguments

• D::Matrix{Float64}: non-trivial distance matrix of the eigenvectors
• 𝚽::Matrix{Float64}: graph Laplacian eigenvectors
• σ::Float64: Gaussian window width parameter (default: 0.25 * maximum(D))

Output Argument

• 𝓤::Matrix{Float64}: the NGWF dictionary

rngwf_all_vectors(D, 𝚽; σ = 0.2 * maximum(D), thres = 0.2)

assemble the reduced NGWF (rNGWF) dictionary.

Input Arguments

• D::Matrix{Float64}: non-trivial distance matrix of the eigenvectors
• 𝚽::Matrix{Float64}: graph Laplacian eigenvectors
• σ::Float64: Gaussian window width parameter (default: 0.25 * maximum(D))
• thres::Float64: cutoff threshold ∈ (0, 1).

Output Argument

• 𝓤::Matrix{Float64}: the rNGWF dictionary
• dic_l2x::Dict: a dictionary to store the filtered locations by QR at the l-th centered eigenvector

rngwf_lx(dic_l2x)

find the sequential subindices of rNGWF vectors.

Input Arguments

• dic_l2x::Dict: a dictionary to store the filtered locations by QR at the l-th centered eigenvector

Output Argument

• Γ::Vector{Tuple{Int64,Int64}}: the sequential subindices of rNGWF vectors.