Gene Cheung
Associate Professor, York University
17th June, 2019

Graph Signal Analysis: Imaging,
Learning, Sampling.
Univ. of Waterloo 6/17/2019

1

Acknowledgement
Collaborators:• X. Liu (HIT, China)• W. Hu, W. Gao (Peking U., China)• L. Fang (Tsinghua, China)• Y. Zhao (BJTU, 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)

Univ. of Waterloo 6/17/2019

2

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

engineering

mathUniv. of Waterloo 6/17/20193

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

engineering

mathUniv. of Waterloo 6/17/20193

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
Univ. of Waterloo 6/17/20194

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
Univ. of Waterloo 6/17/2019

5

Signal Decomposition• Decompose signal into basic components:

x =  X k kk

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

Univ. of Waterloo 6/17/2019

6

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
Univ. of Waterloo 6/17/2019

6

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.
Univ. of Waterloo 6/17/2019

6

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

Univ. of Waterloo 6/17/2019

7

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

1. Data domain is naturally a graph.

signal on graph kernel

2. Embed similarity structure into graph.
•

Eigenvectors have global info aggregated from local info.

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

signal on regular 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.

8

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

1. Data domain is naturally a graph.

signal on graph kernel

2. Embed similarity structure into graph.
•

Eigenvectors have 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.

8

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.

Univ. of Waterloo 6/17/2019*https://en.wikipedia.org/wiki/Second_derivative

9

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

…

8

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.
Univ. of Waterloo 6/17/201910

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

…

8

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.
Univ. of Waterloo 6/17/201910

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

Univ. of Waterloo 6/17/2019*https://en.wikipedia.org/wiki/Delaunay triangulation

11

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

Univ. of Waterloo 6/17/2019*https://en.wikipedia.org/wiki/Delaunay triangulation

11

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

Univ. of Waterloo 6/17/2019*https://en.wikipedia.org/wiki/Delaunay triangulation

11

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

Univ. of Waterloo 6/17/2019*https://en.wikipedia.org/wiki/Delaunay triangulation

11

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
Univ. of Waterloo 6/17/2019

12

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.

13

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.

13

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.

14

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.

15

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.

16

Results: Contrast Enhancement17

Results:
Contrast
Enhancement18

Results: Contrast Enhancement

Univ. of Waterloo 6/17/2019

19

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

20

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

20

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

20

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
21& Total Variation,” IEEE 20th International Workshop on Multimedia Signal Processing, Vancouver, Canada, August 2018.

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

22

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

23

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,” accepted to IEEE International Conference on Image Processing, Taiwan, October 2019.

24

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,” accepted to IEEE International Conference on Image Processing, October 2019.

25

Results: Point Cloud Super-Resolution

26

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
Univ. of Waterloo 6/17/201927

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.

28

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.

28

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.

28

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..

29

Unrolling Graph Laplacian Regularizer

Univ. of Waterloo 6/17/2019[1] J. Zeng et al., “Deep Graph Laplacian Regularization for Robust Denoising of Images,” accepted to NTIRE Workshop, CVPR 2019.

30

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.
Univ. of Waterloo 6/17/2019[1] J. Zeng et al., “Deep Graph Laplacian Regularization for Robust Denoising of Images,” accepted to NTIRE Workshop, CVPR 2019.

31

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.

32

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.

33

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.

DeepGLR34

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.

DeepGLR35

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

Univ. of Waterloo 6/17/2019

36

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.

37

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.

37

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.

38

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.

39

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

Univ. of Waterloo 6/17/201940

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

Univ. of Waterloo 6/17/201940

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

Univ. of Waterloo 6/17/201941

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

Univ. of Waterloo 6/17/201941

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

Univ. of Waterloo 6/17/201942

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.

43

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.

44

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[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.

45

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
𝐧
46

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.
[1] Y. Bai, G. Cheung, F. Wang, X. Liu, W. Gao, "Reconstruction-Cognizant Graph Sampling Using Gershgorin Disc Alignment,"
IEEE International Conference on Acoustics, Speech and Signal Processing, Brighton, UK, May 2019.

47

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

𝐁 = 𝐇𝑇 𝐇 + 𝜇 𝐋

1

1

1

2

3

4

1 −1 0−1 2 −1
𝐁=0 −1 2
0
0 −1

0
0
−1
1

1

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

𝐂 = 𝐒 𝐁 𝐒 −1

similarity transform

1

2

3

4

Sample set { }Scale factor {1,1,1,1}

diagonal matrix w/ scale factors• Scale row → expand disc radius.
→ shrink neighbors’ disc radius.
[1] Y. Bai, G. Cheung, F. Wang, X. Liu, W. Gao, "Reconstruction-Cognizant Graph Sampling Using Gershgorin Disc Alignment,"
IEEE International Conference on Acoustics, Speech and Signal Processing, Brighton, UK, May 2019.

48

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

𝐁 = 𝐇𝑇 𝐇 + 𝜇 𝐋

1

1

1

2

1

4

3

1 −1 00−1 3 −1 0𝐁=0 −1 2 −1
0
0 −1 1coeff. matrix• Sample node → shift disc.
• Consider similar transform of B:

𝐂 = 𝐒 𝐁 𝐒 −1

similarity transform

1

2

3

4

5

Sample set {2}Scale factor {1,1,1,1}

diagonal matrix w/ scale factors• Scale row → expand disc radius.
→ shrink neighbors’ disc radius.
[1] Y. Bai, G. Cheung, F. Wang, X. Liu, W. Gao, "Reconstruction-Cognizant Graph Sampling Using Gershgorin Disc Alignment,"
IEEE International Conference on Acoustics, Speech and Signal Processing, Brighton, UK, May 2019.

48

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

𝐁 = 𝐇𝑇 𝐇 + 𝜇 𝐋

1

1

1

2

1

4

3

1 −1 00−1 3 −1 0𝐁=0 −1 2 −1
0
0 −1 1coeff. matrix• Sample node → shift disc.
• Consider similar transform of B:

𝐂 = 𝐒 𝐁 𝐒 −1

similarity transform

1

2

3

4

5

Sample set {2}Scale factor {1,s2,1,1}

diagonal matrix w/ scale factors• Scale row → expand disc radius.
→ shrink neighbors’ disc radius.
[1] Y. Bai, G. Cheung, F. Wang, X. Liu, W. Gao, "Reconstruction-Cognizant Graph Sampling Using Gershgorin Disc Alignment,"
IEEE International Conference on Acoustics, Speech and Signal Processing, Brighton, UK, May 2019.

48

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

𝐁 = 𝐇𝑇 𝐇 + 𝜇 𝐋

1

1

1

2

1

4

3

1 −1 00−1 3 −1 0𝐁=0 −1 2 −1
0
0 −1 1coeff. matrix• Sample node → shift disc.
• Consider similar transform of B:

𝐂 = 𝐒 𝐁 𝐒 −1

similarity transform

1

2

3

4

5

Sample set {2}Scale factor {1,s2,1,1}

diagonal matrix w/ scale factors• Scale row → expand disc radius.
→ shrink neighbors’ disc radius.
[1] Y. Bai, G. Cheung, F. Wang, X. Liu, W. Gao, "Reconstruction-Cognizant Graph Sampling Using Gershgorin Disc Alignment,"
IEEE International Conference on Acoustics, Speech and Signal Processing, Brighton, UK, May 2019.

48

Gershgorin Disc Alignment (math)

1

1
1

2

1
3

1
4

• Breadth First Iterative Sampling:• 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.

[1] Y. Bai, G. Cheung, F. Wang, X. Liu, W. Gao, "Reconstruction-Cognizant Graph Sampling Using Gershgorin Disc Alignment,"
IEEE International Conference on Acoustics, Speech and Signal Processing, Brighton, UK, May 2019.

49

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.

[1] Y. Bai, G. Cheung, F. Wang, X. Liu, W. Gao, "Reconstruction-Cognizant Graph Sampling Using Gershgorin Disc Alignment,"
IEEE International Conference on Acoustics, Speech and Signal Processing, Brighton, UK, May 2019.

50

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: Choose samples to maximize the ripple
on the least disturbed node.

[1] Y. Bai, G. Cheung, F. Wang, X. Liu, W. Gao, "Reconstruction-Cognizant Graph Sampling Using Gershgorin Disc Alignment,"
IEEE International Conference on Acoustics, Speech and Signal Processing, Brighton, UK, May 2019.

51

Results: Graph Sampling• GDA has larger λmin and smaller condition number.
• GDA has smaller reconstruction MSE.
• GDA has 10% of computation cost of Spectral proxies.

[1] Y. Bai, G. Cheung, F. Wang, X. Liu, W. Gao, "Reconstruction-Cognizant Graph Sampling Using Gershgorin Disc Alignment,"
IEEE International Conference on Acoustics, Speech and Signal Processing, Brighton, UK, May 2019.

52

Results: Graph Sampling• Journal version [1]:1.
2.

formalize the max λmin GDA optimization problem,
derive linear-time heuristic based on set cover approx. alg.

• Small graphs: GDA has roughly the same reconstruction MSE.
• Large graphs: GDA has smallest reconstruction MSE.
[1] Y. Bai, F. Wang, G. Cheung, F. Wang, Y. Nakatsukasa, W. Gao, "Fast Sampling Set Selection for Graph Signals Using Gershgorin
Disc Alignment,“ in preparation for IEEE Transactions on Signal Processing.

53

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.

54

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.
55[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, C. Yang, G. Cheung, Y. Wang, “Graph Sampling for Matrix Completion Using Recurrent First Eigenvector
Computation,” arXiv, preprint arXiv:1906.01087, 2019.

56

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] F. Wang, C. Yang, G. Cheung, Y. Wang, “Graph Sampling for Matrix Completion Using Recurrent First Eigenvector
Computation,” arXiv preprint arXiv:1906.01087, 2019.

57

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

Univ. of Waterloo 6/17/2019

58

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

Univ. of Waterloo 6/17/2019

59

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 Theory
Univ. of Waterloo 6/17/201960