Gene Cheung
Associate Professor, York University
16th October, 2019

Graph Signal Analysis: Imaging,
Learning, Sampling.

Acknowledgement
Collaborators:• X. Liu (HIT, China)• W. Hu, W. Gao (Peking U., China)• L. Fang (Tsinghua, China)• C.-W. Lin (National Tsing Hua University, Taiwan)• A. Ortega (USC, USA)• D. Florencio (MSR, USA)• J. Liang, I. Bajic (SFU, Canada)• X. Wu (McMaster U, Canada)• P. Frossard (EPFL, Switzerland)• V. Stankovic (U of Strathclyde, UK)• Y. Nakatsukasa (Oxford, UK)• P. Le Callet (U of Nantes, France)2

Introducing math tools
Students in EECS4452: “This is math, not engineering!”

engineering

math

3

Introducing math tools
Students in EECS4452: “This is math, not engineering!”
Me: “Math is the heart of engineering!”

engineering

math

4

Outline• Defining Graph frequencies• Inverse Imaging• Image denoising• Image contrast enhancement• 3D point cloud denoising / super-resolution• Deep GLR• Semi-Supervised Learning• Graph Sampling• Matrix completion5Outline• Defining Graph frequencies• Inverse Imaging• Image denoising• Image contrast enhancement• 3D point cloud denoising / super-resolution• Deep GLR• Semi-Supervised Learning• Graph Sampling• Matrix completion6

Signal Decomposition• Decompose signal into basic components:

x =  X k kk

• Newton decomposed white light into color components (1730).

7

Signal Decomposition• Decompose signal into basic components:

x =  X k kk

• Newton decomposed white light into color components (1730).

• “Basic” components can be complex exponentials:

x =  X k e j 2ktk

X k =  x(t )e − j 2kt dt
8Signal Decomposition• Decompose signal into basic components:

x =  X k kk

• Newton decomposed white light into color components (1730).

• “Basic” components can be complex exponentials:

x =  X k e j 2ktk

X k =  x(t )e − j 2kt dt• Complex exponentials are eigenfunctions of 2nd derivative operator.
9

Digital Signal Processing• Discrete signals on regular data kernels.
• Ex.1: audio on regularly sampled timeline.
• Ex.2: image on 2D grid.

• Harmonic analysis tools (transforms, wavelets):• Compression, restoration, segmentation, etc.

a =xsparsetransform coeff.

2D DCT basis

desired signaltransform

10Graph Signal Processing• Signals on irregular data kernels described by graphs.
• Graph: nodes and edges.
• Edges reveals node-to-node relationships.

1. Harmonic Analysis of graph signals.

signal on graph kernel

2. Embed pairwise similarity info into graph.
•

Eigenvectors provide global info aggregated from local info.

Graph Signal Processing (GSP) provides spectral
analysis tools for signals residing on graphs.

signal on graph kernel[1] A. Ortega, P. Frossard, J. Kovacevic, J. M. F. Moura, and P. Vandergheynst, “Graph signal processing: Overview, challenges,
and applications,” Proceedings of the IEEE, vol. 106, no. 5, pp. 808–828, 2018.
[2] G. Cheung, E. Magli, Y. Tanaka, and M. K. Ng, “Graph spectral image processing,” Proceedings of the IEEE, vol. 106, no.
5, pp. 907–930, 2018.

11

GSP and Graph-related Research
GSP: SP framework that unifies concepts from multiple fields.

Partial Differential
Eq’ns

MachineLearning

Laplaceequationgraphical model,
manifold learning,classifier learning

LaplaceBeltrami
operator

Computer Graphics

Graph Signal
Processing* (GSP)

DSP

spectralclustering

Computer
VisionMax cut, graphtransformation

eigen-analysis ofgraph Laplacian,adjacency matrices

Combinatorial
Graph Theory

SpectralGraph Theory12

undirected graph

Graph Fourier Transform (GFT)
Graph Laplacian:• Adjacency Matrix A: entry Ai,j has non-negativeedge weight wi,j connecting nodes i and j.
• Degree Matrix D: diagonal matrix w/ entry Di,ibeing sum of column entries in row i of A.
Di ,i =  Ai , jj• Combinatorial Graph Laplacian L: L = D-A• L is related to 2nd derivative.

L3,: x = − x2 + 2 x3 − x4

w1, 2

1

1

2

3

4

 0
w
A =  1, 2
 0

 0

w1, 2

0 0
1 0
0 1

1 0

 w1, 2
 0
D=
 0

 0

0 0
w1, 2 + 1 0 0
0
2 0

0
0 1

1

 w1, 2− w
L =  1, 2
 0

 0

0
1
0
0

− w1, 20
w1, 2 + 1 − 1 0 −1
2 − 1
0−1 1 
0f (x + h ) − 2 f (x ) + f (x − h )h →0h2f ( x ) = lim• L is a differential operator on graph.

13*https://en.wikipedia.org/wiki/Second_derivativeGraph Spectrum from GFT• Graph Fourier Transform (GFT) is eigen-matrix of graph Laplacian L.
eigenvalues along diagonal

L = V  VT

GFT~
x = VT x

1

w1, 2

1
2

1
3

4

…

1st AC eigenvector

eigenvectors in columns
GFT coefficients

1. Eigenvectors aggregates info from weights.
• Constant eigenvector is DC.
• # zero-crossings increases as λ increases.

2. Eigenvalues (≥ 0) as graph frequencies.
• GFT defaults to DCT for un-weighted connected line.
• GFT defaults to DFT for un-weighted connected circle.
148

Graph Spectrum from GFT• Graph Fourier Transform (GFT) is eigen-matrix of graph Laplacian L.
eigenvalues along diagonal

L = V  VT

GFT~
x = VT x

1

w1, 2

1
2

1
3

4

…

1st AC eigenvector

eigenvectors in columns
GFT coefficients

1. Eigenvectors aggregates info from weights.
• Constant eigenvector is DC.
• # zero-crossings increases as λ increases.

2. Eigenvalues (≥ 0) as graph frequencies.
• GFT defaults to DCT for un-weighted connected line.
• GFT defaults to DFT for un-weighted connected circle.
15

8

Graph Frequency Examples (US Temperature)• Weather stations from 100 most populated cities.
• Graph connections from Delaunay Triangulation*.
• Edge weights inverse proportion to distance.

location diff.
 − l −l 2 
 i j 2wi , j = exp
2

 




Edge weights

V1: DC component

16*https://en.wikipedia.org/wiki/Delaunay triangulation

Graph Frequency Examples (US Temperature)• Weather stations from 100 most populated cities.
• Graph connections from Delaunay Triangulation*.
• Edge weights inverse proportion to distance.

location diff.
 − l −l 2 
 i j 2wi , j = exp
2

 




Edge weights

V2: 1st AC component

17*https://en.wikipedia.org/wiki/Delaunay triangulation

Graph Frequency Examples (US Temperature)• Weather stations from 100 most populated cities.
• Graph connections from Delaunay Triangulation*.
• Edge weights inverse proportion to distance.

location diff.
 − l −l 2 
 i j 2wi , j = exp
2

 




Edge weights

V3: 2nd AC component

18*https://en.wikipedia.org/wiki/Delaunay triangulation

Graph Frequency Examples (US Temperature)• Weather stations from 100 most populated cities.
• Graph connections from Delaunay Triangulation*.
• Edge weights inverse proportion to distance.

location diff.
 − l −l 2 
 i j 2wi , j = exp
2

 




Edge weights

V4: 9th AC component

19*https://en.wikipedia.org/wiki/Delaunay triangulation

Outline• Defining Graph frequencies• Inverse Imaging• Image denoising• Image contrast enhancement• 3D point cloud denoising / super-resolution• Deep GLR• Semi-Supervised Learning• Graph Sampling• Matrix completion
20

Graph Laplacian Regularizer• x T L x (graph Laplacian regularizer) [1]) is one smoothness measure.
1
2
x T Lx =  wi , j (xi − x j ) =  k ~xk2
signal contains
2 i, jk• Signal Denoising:observation• MAP Formulation:
fidelity term

signal smooth innodal domain

y = x+ v

min y − x 2 +  x T L x

desired signal
noise

2

x

update edge
weights

mostly low graph freq.

(I+ L ) x * = ysmoothness priorlinear system of eqn’s w/ sparse, symmetric PD matrix[1] P. Milanfar, “A Tour of Modern Image Filtering: New Insights and Methods, Both Practical and Theoretical,” IEEE Signal
Processing Magazine, vol.30, no.1, pp.106-128, January 2013.

21

Graph Laplacian Regularizer• x T L x (graph Laplacian regularizer) [1]) is one smoothness measure.
1
2
x T Lx =  wi , j (xi − x j ) =  k ~xk2
signal contains
2 i, jk• Signal Denoising:• MAP Formulation:

signal smooth innodal domain

y = x+ v
pixel intensity diff.

min y − x 2 +  x T L x
2x

update edge
weights

mostly low graph freq.

(I+ L ) x * = ypixel location diff.

− x −x 2  − l −l 2 ij i j 22 wi , j = exp
exp

 2

 12




2





Bilateral filter weights

linear system of eqn’s w/ sparse, symmetric PD matrix[1] P. Milanfar, “A Tour of Modern Image Filtering: New Insights and Methods, Both Practical and Theoretical,” IEEE Signal
Processing Magazine, vol.30, no.1, pp.106-128, January 2013.

22

Results: natural image denoising• Subjective comparisons (  I = 40 )Original

Noisy, 16.48 dB

K-SVD, 26.84 dB

BM3D, 27.99 dB

PLOW, 28.11 dB

OGLR, 28.35 dB[1] J. Pang, G. Cheung, "Graph Laplacian Regularization for Image Denoising: Analysis in the Continuous Domain," IEEE
Transactions on Image Processing, vol. 26, no.4, pp.1770-1785, April 2017.

23

Results: depth image denoising• Subjective comparisons (  I = 30 )Original

Noisy, 18.66 dB

BM3D, 33.26 dB NLGBT, 33.41dB OGLR, 34.32 dB[1] W. Hu et al., "Depth Map Denoising using Graph-based Transform and Group Sparsity," IEEE International Workshop on
Multimedia Signal Processing, Pula (Sardinia), Italy, October, 2013.

24

GLR for Joint Dequantization / Contrast Enhancement• Retinex decomposition model:

reflectancenoise

scalarillumination• Objective: general smoothness for luminance, smoothness
w/ negative edges for reflectance.
generalized smooth

piecewise smooth• Constraints: quantization bin constraints• Solution: Alternating accelerated proximal gradient alg [1].

[1] X. Liu, G. Cheung, X. Ji, D. Zhao, W. Gao, "Graph-based Joint Dequantization and Contrast Enhancement ofPoorly Lit JPEG Images," IEEE Transactions on Image Processing, vol. 28, no.3, pp.1205-1219, March 2019.

25

Results: Contrast Enhancement26

Results:
Contrast
Enhancement27

Results: Contrast Enhancement

28

GTV for Point Cloud Denoising• Acquisition of point cloud introduces noise.
• Point cloud is irregularly sampled 2D manifold in 3D space.
• Not appropriate to apply GTV directly on 3D coordinates [1].
• only a singular 3D point has zero GTV value.

 y −y f (x ) − f (x )
i

i −1

i

i −1

i

i

y

f(x)

geometrysmoothness

functionalsmoothnessxx• Proposal: Apply GTV is to the surface normals of 3D point
cloud—a generalization of TV to 3D geometry.

[1] Y. Schoenenberger, J. Paratte, and P. Vandergheynst, “Graph-based denoising for time-varying point clouds,” inIEEE 3DTV-Conference, 2015, pp. 1–4

29

GTV for Point Cloud Denoising• Acquisition of point cloud introduces noise.
• Point cloud is irregularly sampled 2D manifold in 3D space.
• Not appropriate to apply GTV directly on 3D coordinates [1].
• only a singular 3D point has zero GTV value.

 y −y f (x ) − f (x )
i

i −1

i

i −1

i

i

y

f(x)

geometrysmoothness

functionalsmoothnessxx• Proposal: Apply GTV is to the surface normals of 3D point
cloud—a generalization of TV to 3D geometry.

[1] Y. Schoenenberger, J. Paratte, and P. Vandergheynst, “Graph-based denoising for time-varying point clouds,” inIEEE 3DTV-Conference, 2015, pp. 1–4

30

GTV for Point Cloud Denoising• Acquisition of point cloud introduces noise.
• Point cloud is irregularly sampled 2D manifold in 3D space.
• Not appropriate to apply GTV directly on 3D coordinates [1].
• only a singular 3D point has zero GTV value.

 y −y f (x ) − f (x )
i

i −1

i

i

i

i −1

 n −nii −1 1iyf(x)geometrysmoothness

functionalsmoothnessxx• Proposal: Apply GTV is to the surface normals of 3D point
cloud—a generalization of TV to 3D geometry.

[1] Y. Schoenenberger, J. Paratte, and P. Vandergheynst, “Graph-based denoising for time-varying point clouds,” inIEEE 3DTV-Conference, 2015, pp. 1–4

31

PC Denoising Algorithm• Use GTV of surface normals over the K-NN graph:
n𝒊

n𝒋

𝑗
𝑖
• Denoising problem as l2-norm fidelity plus GTV of surface normals:

min q − p 2 +   wi , j n i − n jsmoothness on surface normals

2p,ni , jE

1• Surface normal estimation of n𝒊 is a nonlinear function of p𝒊 and neighbors.

Proposal:
1. Partition point cloud into two independent classes (say red and blue).
2. When computing surface normal for a red node, use only neighboring blue points.

3. Solve convex optimization for red (blue) nodes alternately.
[1] C. Dinesh, G. Cheung, I. V. Bajic, C. Yang, “Fast 3D Point Cloud Denoising via Bipartite Graph Approximation
32& Total Variation,” IEEE 20th International Workshop on Multimedia Signal Processing, Vancouver, Canada, August 2018.

Results: Point Cloud Denoising
Anchor model (𝜎=0.3)

33

Results: Point Cloud Denoising
Daratech model (𝜎=0.3)

34

PC Super-Res Algorithm• Add new interior points to low-res point cloud.
1. Construct triangular mesh using Delaunay triangulation using known points q.
2. Insert new points at the centroids of triangles.

• Partition point cloud into two independent classes (say red and blue).
• When computing normal for a red node, use only neighboring blue points.
• Use graph total variation (GTV) of surface normals over the K-NN graph:mi, j = n i − n jsmoothness on surface normals

min  wi , j n i − n j
p,ni , jE1

 I − B  m   v 
0 C   p  =  q    • Solved via augmented Lagrangian + ADMM.
[1] C. Dinesh, G. Cheung, I. V. Bajic, C. Yang, “3D Point Cloud Super-Resolution via Graph Total Variation on
Surface Normals,” IEEE International Conference on Image Processing, Taiwan, October 2019.

35

Results: Point Cloud Super-Resolution• APSS and RIMLS schemes generate overly smooth models.
• Existing methods result in distorted surfaces with somedetails lost.

[1] C. Dinesh, G. Cheung, I. V. Bajic, “3D Point Cloud Super-Resolution via Graph Total Variation on Surface
Normals,” IEEE International Conference on Image Processing, October 2019.

36

Results: Point Cloud Super-Resolution37

Outline• Defining Graph frequencies• Inverse Imaging• Image denoising• Image contrast enhancement• 3D point cloud denoising / super-resolution• Deep GLR• Semi-Supervised Learning• Graph Sampling• Matrix completion
38

Unrolling Graph Laplacian Regularizer• Recall MAP formulation of denoising problem with quadratic
graph Laplacian regularizer:

min y − x 2 +  x T L x
2x
fidelity termsmoothness prior• Solution is system of linear equations:(I+ L ) x * = ylinear system of eqn’s w/ sparse, symmetric PD matrix[1] J. Pang, G. Cheung, "Graph Laplacian Regularization for Inverse Imaging: Analysis in the Continuous Domain," IEEE
Transactions on Image Processing, vol. 26, no.4, pp.1770-1785, April 2017.

39

Unrolling Graph Laplacian Regularizer• Recall MAP formulation of denoising problem with quadratic
graph Laplacian regularizer:

min y − x 2 +  x T L x
2x
fidelity termsmoothness prior• Solution is system of linear equations:(I+ L ) x * = ylinear system of eqn’s w/ sparse, symmetric PD matrix

Q: what is the “most appropriate” graph?

[1] J. Pang, G. Cheung, "Graph Laplacian Regularization for Inverse Imaging: Analysis in the Continuous Domain," IEEE
Transactions on Image Processing, vol. 26, no.4, pp.1770-1785, April 2017.

40

Unrolling Graph Laplacian Regularizer• Recall MAP formulation of denoising problem with quadratic
graph Laplacian regularizer:

min y − x 2 +  x T L x
2x
fidelity termsmoothness prior• Solution is system of linear equations:(I+ L ) x * = ylinear system of eqn’s w/ sparse, symmetric PD matrix

Q: what is the “most appropriate” graph?

Bilateral weights:
− x −x 2  − l −l 2 
i
j 2 

 i j 2wi , j = exp
exp

 2

 12




2





[1] J. Pang, G. Cheung, "Graph Laplacian Regularization for Inverse Imaging: Analysis in the Continuous Domain," IEEE
Transactions on Image Processing, vol. 26, no.4, pp.1770-1785, April 2017.

41

Unrolling Graph Laplacian Regularizer• Deep Graph Laplacian Regularization:
1. Learn features f’s using CNN.
2. Compute distance from features.
3. Compute edge weights using Gaussian kernel.
4. Construct graph, solve QP.
featurespre-filter

weight para

[1] M. McCann et al., “Convolutional Neural Networks for Inverse Problems in Imaging,“ IEEE SPM, Nov. 2017.
[2] K. Gregor and Y. LeCun, “Learning fast approximations of sparse coding,” in Proc. 27th Int. Conf. Machine Learning, 2010..

42

Unrolling Graph Laplacian Regularizer

43[1] J. Zeng et al., “Deep Graph Laplacian Regularization for Robust Denoising of Images,” NTIRE Workshop, CVPR 2019.

Unrolling Graph Laplacian Regularizer• Graph Model guarantees numerical stability of solution:(I+ L ) x = y*• Thm 1: condition number κ of matrix satisfies [1]:maximum node degree• Observation: By restricting search space of CNN to degree-bounded graphs,
we achieve robust learning.
44[1] J. Zeng et al., “Deep Graph Laplacian Regularization for Robust Denoising of Images,” NTIRE Workshop, CVPR 2019.

Experimental Results – Numerical Comparison• Trained on AWGN on 5 images, patches of size 26-by-26.
• Batch size is 4, model is trained for 200 epochs.
• Trained for both known and blind noise variance.

[1] Kai Zhang et al, “Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising,” TIP 2017.
[2] Marc Lebrun et al, “The noise clinic: a blind image denoising algorithm,” IPOL 2015.

45

Experimental Results – Numerical Comparison• Cross-domain generalization.
• trained on Gaussian noise, tested on low-light images in (RENOIR).
• Competing methods: DnCNN [1], noise clinic [2].
• outperformed DnCNN by 5.74 dB, and noise clinic by 1.87 dB.

[1] Kai Zhang et al, “Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising,” TIP 2017.
[2] Marc Lebrun et al, “The noise clinic: a blind image denoising algorithm,” IPOL 2015.

46

Experimental Results – Visual Comparison• trained on Gaussian noise, tested on low-light images in (RENOIR).
• Competing methods: DnCNN [1], noise clinic [2].
• outperformed DnCNN by 5.74 dB, and noise clinic by 1.87 dB.

Noise Clinic

CDnCNN[1] Kai Zhang et al, “Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising,” TIP 2017.
[2] Marc Lebrun et al, “The noise clinic: a blind image denoising algorithm,” IPOL 2015.

DeepGLR47

Experimental Results – Visual Comparison• trained on Gaussian noise, tested on low-light images in (RENOIR).
• Competing methods: DnCNN [1], noise clinic [2].
• outperformed DnCNN by 5.74 dB, and noise clinic by 1.87 dB.

Noise Clinic

CDnCNN[1] Kai Zhang et al, “Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising,” TIP 2017.
[2] Marc Lebrun et al, “The noise clinic: a blind image denoising algorithm,” IPOL 2015.

DeepGLR48

Outline• Defining Graph frequencies• Inverse Imaging• Image denoising• Image contrast enhancement• 3D point cloud denoising / super-resolution• Deep GLR• Semi-Supervised Learning• Graph Sampling• Matrix completion

49

Semi-Supervised Graph Classifier Learning• Binary Classifier: given feature vector xi ofdimension K, compute f(xi) ∊ {0,1}.

• Classifier Learning: given partial, noisy labels (xi,yi), train classifier f(xi).

• GSP Approach [1]:example graph-based classifier

1. Construct signed similarity graph with +/- edges.
2. Pose MAP graph-signal restoration problem.
3. Perturb graph Laplacian to ensure PSD.
4. Solve num. stable MAP as sparse lin. system.

[1] Yu Mao, Gene Cheung, Chia-Wen Lin, Yusheng Ji, “Image Classifier Learning from Noisy Labels via Generalized Graph
Smoothness Priors,” IEEE IVMSP Workshop, Bordeaux, France, July 2016. (Best student paper award)[2] G. Cheung, W.-T. Su, Y. Mao, C.-W. Lin, "Robust Semi-Supervised Graph Classifier Learning with Negative Edge Weights," IEEE
Transactions on Signal and Information Processing over Networks, vol. 4, no.4, pp.712-726, December 2018.

50

Semi-Supervised Graph Classifier Learning• Binary Classifier: given feature vector xi ofdimension K, compute f(xi) ∊ {0,1}.

• Classifier Learning: given partial, noisy labels (xi,yi), train classifier f(xi).

• GSP Approach [1]:-1

example graph-based classifier

1. Construct signed similarity graph with +/- edges.
2. Pose MAP graph-signal restoration problem.
3. Perturb graph Laplacian to ensure PSD.
4. Solve num. stable MAP as sparse lin. system.

[1] Yu Mao, Gene Cheung, Chia-Wen Lin, Yusheng Ji, “Image Classifier Learning from Noisy Labels via Generalized Graph
Smoothness Priors,” IEEE IVMSP Workshop, Bordeaux, France, July 2016. (Best student paper award)[2] G. Cheung, W.-T. Su, Y. Mao, C.-W. Lin, "Robust Semi-Supervised Graph Classifier Learning with Negative Edge Weights," IEEE
Transactions on Signal and Information Processing over Networks, vol. 4, no.4, pp.712-726, December 2018.

51

Graph-Signal Smoothness Priorfor signed graphs

w=1• Graph Laplacian Regularizer [1]:GFT coeff

eigenvalues / graph freqs• Promote large / small inter-node differences depending on edge
signs.

Promote large difference

Promote small difference• Sensible, but numerically unstable.

[1] J. Pang and G. Cheung, “Graph Laplacian regularization for image denoising: Analysis inn the continuous domain,” in IEEE
Transactions on Image Processing, vol. 26, no.4, April 2017, pp. 1770–1785.

52

Semi-Supervised Learning Formulation• MAP formulation:

l0 fidelity term

perturbation matrix
to ensure PSD!

smoothness prior for
signed graph• One sol’n is △=λmin I, i.e. shift all eigenvalues up by η=λmin.
• Intuition: signal variations + signal energies[1] G. Cheung, W.-T. Su, Y. Mao, C.-W. Lin, "Robust Semi-Supervised Graph Classifier Learning with Negative Edge Weights," IEEE
Transactions on Signal and Information Processing over Networks, vol. 4, no.4, pp.712-726, December 2018.

53

Results: Semi-Supervised Learning• Comparisons w/ other classifiers:

54

Results: Semi-Supervised Learning• Comparisons w/ other classifiers:

55

Results: Semi-Supervised Learning• Comparisons w/ other classifiers:

56

Results: Semi-Supervised Learning• Comparisons w/ other classifiers:57

Outline• Defining Graph frequencies• Inverse Imaging• Image denoising• Image contrast enhancement• 3D point cloud denoising / super-resolution• Deep GLR• Semi-Supervised Learning• Graph Sampling• Matrix completion

58

Graph Sampling (with and without noise)• Q: How to choose best samples for graph-based reconstruction?

• Existing graph sampling strategies extend
Nyquist sampling to graph data kernels:• Assume bandlimited signal.
• Greedily select most “informative” samples bycomputing extreme eigenvectors of sub-matrix.
• Computation-expensive.

[1] A. Anis, A. Gadde, and A. Ortega, “Efficient sampling set selection for bandlimited graph signals using graph spectral proxies,”
IEEE Transactions on Signal Processing, vol. 64, no. 14, pp. 3775–3789, 2016.

59

Related Works
Aggressive sampling [1]A-optimal(M. Tsitsvero.TSP2016; F. Wang. SPL2018)

E-optimal(S. Chen. TSP2015)

Graph sampling

Bandlimitedness

Local measurement [2]Spectral proxies(A. Anis. TSP2016; A. Anis. TIT2018)

Deterministic

Localized coverage(A. Sakiyama. TSP2019)

Selection sampling

Smoothness: GMRF model
Random selection [3](P. Chen. ICASSP2018; Y. Bai. ICASSP 2019)

Eigen-decomposition Free

[1] A. G. Marques, S. Segarra, G. Leus, and A. Ribeiro, “Sampling of graph signals with successive local aggregations.” IEEE
Transactions on Signal Processing, vol. 64, no. 7, pp. 1832–1843, 2016.
[2] X. Wang, J. Chen, and Y. Gu, “Local measurement and reconstruction for noisy bandlimited graph signals,” Signal Processing, vol.
129, pp. 119–129, 2016.
[3] G. Puy, N. Tremblay, R. Gribonval, and P. Vandergheynst, “Random sampling of bandlimited signals on graphs,” Applied andComputational Harmonic Analysis, vol. 44, no. 2, pp. 446–475, 2018.

60

Signal Reconstruction using GLR
sampling matrix
desired signal

observation• Signal Model:

𝐲 = 𝐇𝐱 + 𝐯

noise

1

1

1

1

2

3

4

𝐇=

0
0

1 0
0 0

0
1

Sample set {2, 4}• Signal prior is graph Laplacian regularizer (GLR) [1]:1
2𝐱 Lx = ෍ 𝑤𝑖,𝑗 𝑥𝑖 − 𝑥𝑗 = ෍ 𝜆𝑘 𝑥෤𝑘2
2
𝑇

𝑖,𝑗

𝑘

signal smooth w.r.t. graphsignal contains
mostly low graph freq.

• MAP Formulation:
signal priorfidelity term

min 𝐲 − 𝐇𝐱 22 + 𝜇 𝐱 𝑇 𝐋𝐱
𝐱

𝐇𝑇𝐇 + 𝜇𝐋 𝐱 ∗ = 𝐲linear system of eqn’s solved using conjugate gradient
61

1

1

Stability of Linear System

1• Examine system of linear equations :

1

2

4

3

1 −1 00−1 2 −1 0𝐋=
0 −1 2 −1
0
0 −1 1

𝐇𝑇𝐇 + 𝜇𝐋 𝐱 ∗ = 𝐲coefficient matrix B• Stability depends on the condition number(λmax/ λmin) of coeff. matrix B.

0
0
𝐇T𝐇 =
0
0

0
1
0
0

0
0
0
0

0
0
0
1

Sample set {2, 4}• λmax is upper-bounded by 1+μ2*dmax.
• Goal: select samples to maximize λmin(without computing eigen-pairs)!

• Also minimizes worst-case MSE:𝐱ො − 𝐱 2 ≤ 𝜇

1𝜆𝑚𝑖𝑛 𝐁

෥
𝐋 𝐱+𝐧
2

2+

෥ 2
𝐧
62

Gershgorin Circle Theorem• Gershgorin Circle Theorem:• Row i of L maps to a Gershgorin
disc w/ centre Lii and radius Ri

1

1

1

2

3

4

1 −1 0−1 2 −1
𝐋=
0 −1 2
0
0 −1

0
0
−1
1

1

𝑅𝑖 = ෍ 𝐿𝑖𝑗𝑗≠𝑖• λmin is lower-bounded by smallest
left-ends of Gershgorin discs:

1

2

3

4

min 𝐿𝑖,𝑖 − 𝑅𝑖 ≤ 𝜆min𝑖

• Graph Laplacian L has all Gershgorin discleft-ends at 0 → L is psd.
63

Gershgorin Disc Alignment• Main Idea: Select samples to max smallest
disc left-end of coefficient matrix B:

𝐁 = 𝐇𝑇 𝐇 + 𝜇 𝐋

1

1

1

2

1

11
−1
−1
𝐁
𝐁=
=
00
00

4

3

−1
00
−1 00
33
−1 00
2 −1
−1
−1 22 −1
−1
00 −1
−1 11

coeff. matrix• Sample node → shift disc.
• Consider similar transform of B:

𝐂 = 𝐒 𝐁 𝐒 −1

similarity transform

11

22

33

44

55

SampleSample set
set {2}
{2}
{}
Scale
Scale factor
factor {1,s{1,1,1,1}
2,1,1}

diagonal matrix w/ scale factors• Scale row → expand disc radius.
→ shrink neighbors’ disc radius.
64

Aligning discs at threshold T

1

1
1

2

1
3

1
4

• Breadth First Iterative Sampling (BFIS):• Given initial node set, threshold T.
1.

Sample chosen node i(shift disc)2.
3.

Scale row i(expand disc radius i to T )If disc left-end of connected node j >T,
Scale row j

Else,

4.
5.

(expand disc radius j to T )Add node j to node set.

Goto 1 if node set not empty.
Output sample set and count K.
65

5

Gershgorin Disc Alignment (math)• Binary Search with BFIS:• Sample count K inverse proportional to threshold T.
• Binary search on T to drive count K to budget.

• Example: line graph with equal edge weight.
• Uniform sampling.

66

Disc-based Sampling (intuition)• Analogy: throw pebbles into a pond.
• Disc Shifting: throw pebble at sample node i.
• Disc Scaling: ripple to neighbors of node i.
• Goal: Select min # of samples so ripple ateach node is at least T.

67

Results: Graph Sampling• GDA is 100x to 1000x faster than state-of-art methods computing e-vectors.
• GDA is “comparable” in complexity to Random [23] and Ed-free [8].

[1] Yuanchao Bai, Fen Wang, Gene Cheung, Yuji Nakatsukasa, Wen Gao, "Fast Graph Sampling Set Selection Using Gershgorin Disc
Alignment," submitted to IEEE Transactions on Signal Processing, July 2019.

68

Results: Graph Sampling• Small graphs: GDA has roughly the same reconstruction MSE.
• Random sensor graph of size 500 for two signal types.

[1] Yuanchao Bai, Fen Wang, Gene Cheung, Yuji Nakatsukasa, Wen Gao, "Fast Graph Sampling Set Selection Using Gershgorin Disc
Alignment," submitted to IEEE Transactions on Signal Processing, July 2019.

69

Results: Graph Sampling• Large graphs: GDA has smallest reconstruction MSE.
• Minnesota road graph of size 2642 and for two signal types.

[1] Yuanchao Bai, Fen Wang, Gene Cheung, Yuji Nakatsukasa, Wen Gao, "Fast Graph Sampling Set Selection Using Gershgorin Disc
Alignment," submitted to IEEE Transactions on Signal Processing, July 2019.

70

Matrix Completion• Fill in missing entries in a matrix:

1

2
4(Low-rank matrix recovery problem)

3
2

min
rank (X )
M N

3

XR

s.t. X i , j = M i , j , i, j  S• Examples of applications:• Recommendation system—making rating prediction.
• Remote sensing—infer full covariance matrix from partial correlations.
• Structure-from-motion in computer vision.

[1] E. Candes and Y. Plan, “Matrix completion with noise,” Proceedings of the IEEE, vol. 98, no.6, April 2010, pp. 925–936.

71

Synthetic Netflix matrix

Matrix Completion• Convex relaxation to nuclear norm:

min
X*
M N

XR

s.t. X i , j = M i , j , i, j  SGraph Fourier transform on L r Graph Fourier transform on L c• Proximal Gradient: SVD plus singular value soft-thresholding.

• Use dual graph-signal smoothness prior to promote low rank [1]:

(

)

(

)

min
tr X L r X +  tr X L c X +  S M − S X F
M N

XR

T

T

2

• Unconstrained convex objective solvable via ADMM, conjugate gradient.
72[1] V. Kalofolias, X. Bresson, M. Bronstein, and P. Vandergheynst, “Matrix completion on graphs,”arXiv preprint arXiv:1408.1717, 2014.

Results: Sampling for matrix completion(1) Noiseless(2) Noisy

Noiseless synthetic rating matrix

Figure: Reconstruction MSE of different sampling methods on synthetic dataset. The
reconstruction method for matrix completion is dual graph smoothness based method.
•••

Noisy synthetic rating matrix

Comparison methods: PG [1]; GWC-random [2]; LOC [3][1] Guillermo Ortiz-Jiménez, Mario Coutino, Sundeep Prabhakar Chepuri, and Geert Leus. “Sampling andreconstruction of signals on product graphs”. arXiv preprint arXiv:1807.00145, 2018.
[2] G. Puy, N. Tremblay, R. Gribonval, and P. Vandergheynst, “Random sampling of bandlimited signals on graphs,”
Applied and Computational Harmonic Analysis, vol. 44, no. 2, pp. 446–475, 2018.
[3] A. Sakiyama, Y. Tanaka, T. Tanaka, and A. Ortega, “Eigendecomposition-free sampling set selection for graphsignals,”IEEE Transactions on Signal Processing, 2019.

[1] F. Wang, Y. Wang, G. Cheung, C. Yang, "Graph Sampling for Matrix Completion Using Recurrent Gershgorin Disc Shift,"
submitted to IEEE Transactions on Signal Processing, October 2019.

73

Results: Sampling for matrix completion
SVT [1]GRALS [2]

GMC [3]

NMC [4]G1

1.021 | 1.031

0.947 | 0.931

1.036 | 1.037 0.890 | 0.888

G2

1.021 | 0.983

0.945 | 0.893

1.118 | 1.054 0.890 | 0.858

Table: RMSE of different matrix completion methods on Mocielens_100k dataset withdifferent sampling strategies on Feature-based graph (G1) and Content-based graph(G2). In each grid, the value on left side belongs to random sampling; the right side value
is of our proposed IGCS sampling. The best performance in each row is marked in bold
and red. In our experiments, the sampling budget is 80k out of 100k available ratings; We
first use random 60k samples as given, and then proceed to sample the next 20k samples
base on random sampling or the proposed IGCS sampling.

•
•
•
•

[1] J. Cai, E. J. Candes, and Z. Shen. “A singular value thresholding algorithm for matrix completion”. preprint, 2008.
[2] N. Rao, H.-F. Yu, P. K. Ravikumar, and I. S. Dhillon. “Collaborative filtering with graph information: Consistency
and scalable methods”. In Proc. NIPS, 2015.
[3] V. Kalofolias, X. Bresson, M. M. Bronstein, and P. Vandergheynst.” Matrix completion on graphs. ”2014.
[4] D. M. Nguyen, E. Tsiligianni, and N. Deligiannis, “Extendable neural matrix completion,” in Proc. IEEE Int. Conf.
Acoust., Speech Signal Process., 2018, pp. 1–5.

[1] [1] F. Wang, Y. Wang, G. Cheung, C. Yang, "Graph Sampling for Matrix Completion Using Recurrent Gershgorin Disc Shift,"
submitted to IEEE Transactions on Signal Processing, October 2019.

74Summary• Graph Spectral Analyis Tools• Similarity graph, graph frequencies.
Image Processing
Denoising, contrast
enhancement, PC
denoising / SR

Maching Learning
Semi-supervised graphclassifier learning

GSP

Data Mining
Fast sampling for matrix
completion

75

Q&A
• Email: genec@yorku.ca• Homepage: https://www.eecs.yorku.ca/~genec/index.html

76Primal Sample Selection Problem• Optimization: Select sample vector a and scalars s:
smallest disc left-end of C

C is similar transform of coeff. matrix
sample vector a is binary
and within budget K
scalars s are positive• Difficulty: max-min objective is hard to optimize.

77

Dual Sample Selection Problem• Dual Formulation: Select sample vector a and scalars s:
total number of samples

all disc left-ends are at least T• Proposition: If there exists threshold T s.t. optimal sol’n (a,s) to dual
satisfies ∑ai = K, one dual sol’n is also optimal to primal.

78