Nonnegative Tensor Factorization using a proximal algorithm: - - PowerPoint PPT Presentation

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Nonnegative Tensor Factorization using a proximal algorithm: application to 3D fluorescence spectroscopy

Caroline Chaux

Joint work with X. Vu, N. Thirion-Moreau and S. Maire (LSIS, Toulon) Aix-Marseille Univ. I2M

IHP, Feb. 2 2017

1 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Outline

Introduction 3D fluorescence spectroscopy Tensor definition Goal Proximal tools Criterion formulation Proximity operator Proximal algorithm Application to CPD Minimization problem Algorithm Numerical simulations Synthetic case Real case: water monitoring Conclusion and future work

2 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

3D fluorescence spectroscopy

λem Normalized emission spectra

λex Normalized excitation spectra

Experiment Normalized concentrations

3 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Coumpounds characterisation

300 320 340 360 380 400 420 440 460 480 500 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04

λem Normalized emission spectra

200 220 240 260 280 300 320 340 360 380 400 0.005 0.01 0.015 0.02 0.025 0.03 0.035

λex Normalized excitation spectra

λex FEEM

220 240 260 280 300 320 340 360 380 400 350 400 450 500

λex

220 240 260 280 300 320 340 360 380 400 350 400 450 500

λex λem

220 240 260 280 300 320 340 360 380 400 350 400 450 500 1 2 3 4 5 6 7 8

4 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Tensor

An Nth-order tensor is represented by an N-way array in a chosen basis. What is a tensor? Example:

◮ N = 1: a vector. ◮ N = 2: a matrix.

5 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Third-order tensors

◮ A special case: nonnegative third-order tensors (N = 3)

T = (ti1i2i3)i1,i2,i3 ∈ R+I1×I2×I3.

◮ The Canonical Polyadic (CP) decomposition:

Tensor rank T = R

• r=1

¯ a(1)

r

• ¯

a(2)

r

• ¯

a(3)

r

= [ [¯ A(1), ¯ A(2), ¯ A(3)] ] Loading vectors Loading matrices ∀n ∈ {1, 2, 3}, ¯ a(n)

r

∈ R+In and ¯ A(n) ∈ RIn×R

• : the outer product.

◮ Entry-wise form:

ti1i2i3 =

R

• r=1

¯ a(1)

i1r ¯

a(2)

i2r ¯

a(3)

i3r ,

∀(i1, i2, i3)

6 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Standard operations

◮ Outer product: let u ∈ RI, v ∈ RJ,

u◦v = uv⊤ ∈ RI×J

◮ Khatri-Rao product: let U = [u1, u2, . . . , uJ] ∈ RI×J and

V = [v1, v2, . . . , vJ] ∈ RK×J U⊙V = [u1⊗v1, u2⊗v2, . . . , uJ⊗vJ] ∈ RIK×J. where u⊗v = [u1v; . . . ; uIv] ∈ RIK (Kronecker product).

◮ Hadamard division: let U ∈ RI×J, V ∈ RI×J,

U⊘V = (uij/vij)i,j ∈ RI×J

7 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Tensor flattening: example

Objective: to handle matrices instead of tensors.

Tensor Mode 1 Mode 2 Mode 3

8 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

3D fluorescence spectroscopy and tensors

T =

R

• r=1

¯ a(1)

r

• ¯

a(2)

r

• ¯

a(3)

r

λem Normalized emission spectra

λex Normalized excitation spectra

Experiment Normalized concentrations

9 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Objective: tensor decomposition

◮ Input:

◮ Observed tensor T : observation of an original (unknown) tensor T

possibly degraded (noise).

◮ Output:

◮ Estimated loading matrices

A(n) for all n ∈ {1, 2, 3}

◮ Difficulty:

◮ Rank R unknown (i.e.

R = R): generally i) estimated or ii) overestimated.

Formulate the problem under a variational approach. Proposed approach

10 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Minimization problem

◮ Standard problem:

minimize

x∈RL

F(x)

• Fidelity

+ R(x)

Regularization

.

◮ Taking into account several regularizations (J terms):

R(x) =

J

• j=1

Rj(x)

◮ For large size problem or for other reasons, can be interesting to work

• n data blocks x(j) of size Lj (x = (x(j))1≤j≤J)

R(x) =

J

• j=1

Rj(x(j)) Technical assumptions: F, R and Rj are proper lower semi-continuous

• functions. F is differentiable with a β-Lipschitz gradient. Rj is assumed to

be bounded from below by an affine function, and its restriction to its domain is continuous.

11 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Proximity operator

◮ let ϕ : R → ]−∞, +∞] be a proper lower semi-continuous function.

The proximity operator is defined as proxϕ : R → R: v → arg min

u∈R

1 2 u − v2 + ϕ(u),

◮ let ϕ : RL → ]−∞, +∞] be a proper lower semi-continuous function.

The proximity operator associated with a Symmetric Positive Definite (SPD) matrix P is defined as proxP,ϕ : RL → RL : v → arg min

u∈RL

1 2 u − v2

P + ϕ(u),

where ∀x ∈ RL, x2

P = x, Px and ·, · is the inner product.

Remark : Note that if P reduces to the identity matrix, then the two definitions coincides.

12 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Criterion to be minimized

minimize

x∈RL

F(x) +

J

• j=1

Rj(x(j)) Some solutions (non exhaustive list, CPD oriented):

◮ Proximal Alternating Linearized Minimization (PALM) [Bolte et al., 2014] ◮ A Block Coordinate Descent Method for both CPD and Tucker

decomposition [Xu and Yin, 2013]

◮ An accelerated projection gradient based algorithm [Zhang et al., 2016] ◮ Block-Coordinate Variable Metric Forward-Backward (BC-VMFB)

algorithm [Chouzenoux et al., 2016] Advantages of the BC-VMFB: flexible, stable, integrates preconditionning, relatively fast.

13 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Block coordinate proximal algorithm

1: Let x0 ∈ domR, k ∈ N and γk ∈]0, +∞[

// Initialization step

2: for k = 0, 1, ... do

// k-th iteration of the algorithm

3:

Let jk ∈ {1, ..., J} // Processing of block number jk (chosen, here, according to a quasi cyclic rule)

4:

Let Pjk(xk) be a SPD matrix // Construction of the preconditioner Pjk(xk)

5:

Let ∇jkF(xk) be the Gradient // Calculation of Gradient

6:

˜ x(jk)

k

= x(jk)

k

− γkPjk(xk)−1∇jkF(xk) // Updating of block jk according to a Gradient step

7:

x(jk)

k+1 ∈ proxγ−1

k

Pjk (xk),Rjk

• ˜

x(jk)

k

• // Updating of block jk according to a

Proximal step

8:

x

¯ jk k+1 = x ¯ jk k where¯

j = {1, ..., J} \ {j} // Other blocks are kept unchanged

9: end for

14 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Prox for CP decomposition

CP decomposition: decompose a tensor into a (minimal) sum of rank-1 terms. Order 3: T =

R

• r=1

¯ a(1)

r

• ¯

a(2)

r

• ¯

a(3)

r

= [ [¯ A(1), ¯ A(2), ¯ A(3)] ], (1) Tensor structure: naturally leads to consider 3 blocks corresponding to the loading matrices A(1), A(2) and A(3). minimize

A(n)∈RIn×R, n∈{1,2,3} F(A(1), A(2), A(3))+R1(A(1))+R2(A(2))+R3(A(3)).

Proposed optimization problem Some of the fastest classical approaches: Fast HALS [Phan et al., 2013] and N-Way [Bro, 1997].

15 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Tensor matricization

◮ T (n) In,I−n ∈ RIn×I−n +

the matrix obtained by unfolding the tensor T in the n-th mode where the size I−n is equal to I1I2I3/In

◮ Tensor expressed under matrix form as

T

(n) In,I−n = ¯

A(n)(Z

(−n))⊤,

n ∈ {1, 2, 3} where Z

(−1) = ¯

A(3) ⊙ ¯ A(2) ∈ RI−1×R

+

, Z

(−2) = ¯

A(3) ⊙ ¯ A(1) ∈ RI−2×R

+

, Z

(−3) = ¯

A(2) ⊙ ¯ A(1) ∈ RI−3×R

+

,

16 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Function choice

◮ F(A(1), A(2), A(3)): quadratic data fidelity term

F(A(1), A(2), A(3)) = 1 2T −[ [A(1), A(2), A(3)] ]2

F = 1

2T(n)

In,I−n−A(n)Z(−n)⊤2 F ◮ Rn(A(n)): block dependent penalty terms enforcing sparsity and

nonnegativity Rn(A(n)) =

In

• in=1

R

• r=1

ρn(a(n)

inr )

∀n ∈ {1, 2, 3} where loading matrices are defined element wise as A(n) = (a(n)

inr )(in,r)∈{1,...,In}×{1,...,R} and

ρn(ω) =

• α(n)|ω|π(n)

if η(n)

min ≤ ω ≤ η(n) max

+∞

• therwise

α(n) ∈]0, +∞[, π(n) ∈ N∗, η(n)

min ∈ [−∞, +∞[ and η(n) max ∈ [η(n) min, +∞].

⇒ block dependent but constant within a block regularization parameters.

17 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Preconditionning

Preconditionning similar to the one used in NMF [Lee and Seung, 2001]. The matrix P for the n-th block can be defined as follows ∀n ∈ {1, 2, 3} P(n)(A(1), A(2), A(3)) = A(n)(Z(−n)⊤Z(−n)) ⊘ A(n) , Remark: ∀n ∈ {1, 2, 3}, A(n) must be non zero.

18 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Gradient and proximity operator

◮ Gradient matrices of F with respect to A(n) for all n = 1, . . . , 3, defined

as ∇nF(A(1), A(2), A(3)) = −(T(n)

In,I−n − A(n)Z(−n)⊤)Z(−n). ◮ Proximity operator given by (∀y = (y(i))i∈{1,...,RIn} ∈ RRIn))

proxγ[k]−1P(n)[k],Rn(y) =

• proxγ[k]−1p(n)

i

[k],ρn(y(i))

• i∈{1,...,RIn} .

where ∀i ∈ {1, ..., RIn}, we have (∀υ ∈ R) proxγ[k]−1p(n)

i

,ρn(υ) = min

• η(n)

max, max

• η(n)

min, proxγ[k]α(n)(p(n)

i

[k])−1| . |π(n) (υ)

• (separable structure, diagonal preconditionning matrices,

componentwise calculation)

19 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Proximal algorithm for tensor decomposition

matrices An[0] Inputs: 1) Initial loading 2) Stepsize choice γ 3) Iteration k = 0 Yes Choose randomly a block n ∈ {1, 2, 3} to be updated Preconditioner P(n)[k] Partial gradient ∇n[k] Gradient step Proximal step Is stopping criterion reached ? proxγ−1P(n)[k],Rn

• A(n)[k] = A(n)[k]

−γ∇n[k] ⊘ P(n)[k] ( A(n)[k]) Other blocks unchanged at k matrices A(n) Estimated loading Outputs: No k ← k + 1

Figure: BC-VMFB algorithm for CPD.

20 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Computer simulation: simulated spectroscopy-like data

◮ Simulated tensor: (uni or bimodal type) emission and excitation

spectra, random concentrations ⇒ T ∈ R100×100×100

+

and R = 5.

◮ Simulated observed tensor: T = T + B where B stands for an additive

white Gaussian noise

◮ 2 considered cases :

• 1. Perturbed case (noiseless): no noise added and

R = 6 (overestimation).

• 2. Perturbed case (noisy): B fixed such that SNR = 17.6 dB and

R = 6 (overestimation).

◮ Error measures

• 1. Signal to Noise Ratio defined as SNR = 20 log10

T F T − T F

• 2. Relative Reconstruction Error defined as RRE = 20 log10

T − T 1 T 1

• 3. Estimation error: E1 = 10 log10

3

n=1

A(n)(1 : R) − ¯ A(n)1 3

n=1 ¯

A(n)1

• 4. Over-factoring error: E2 = 10 log10

 

• R
• r=R+1
• a(1)

r

a(2)

r

a(3)

r 1

 

21 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Numerical results

Elapsed time (s) BC-VMFB without penalty BC-VMFB with penalty N-way fast HALS For 50 iterations 0.2 0.2 11 0.5 Noisy case To reach stopping conditions 102 75 8 8 (actual number of iterations) (48500) (36500) (43) (1856) (SNR,E1, E2) dB (31.3, -12.5, 30.6) (32.7, -11.2, -409) (31.3, -12.5, 30.6) (31.3, -12.5, 30.6) Noiseless case To reach stopping conditions 202 74 80 3.7 (actual number of iterations) (100000) (36500) (838) (308) (RRE,E1, E2) dB (-75.1,-12.4,25.6) (-44.7, -15, -409) (-127.9,-8.7, 31.7) (-63.9, -6.1, 31.7)

Computation time comparison of BC-VMFB in two cases: with or without penalty, with N-way [Bro, 1997] and fast HALS [Phan et al., 2013] using the same initial value in the noiseless and noisy cases.

22 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Influence of the initialization

10 20 30 40 50 60 70 80 90 100 −14 −12 −10 −8 −6 −4 −2 Initializations dB

Error

BC−VMFB with penalty BC−VMFB without penalty Nway fHALS 10 20 30 40 50 60 70 80 90 100 −450 −400 −350 −300 −250 −200 −150 −100 −50 50 Initializations dB

Over−factoring index

BC−VMFB with penalty BC−VMFB without penalty Nway fHALS

Performance versus different initializations (noisy, overestimated case): error index E1, overfactoring error index E2

23 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Visual results: noiseless case

Reference

λex 400 500 λex 400 500 λex 400 500 λex 400 500 λex 400 500 λem λex 350 400 450 500 400 500

BC−VMFB without penalty

λem 350 400 450 500

BC−VMFB with penalty

λem 350 400 450 500 0.05 0.5 5

Figure: FEEM of reference (left) - FEEM reconstructed using BC-VMFB without regularization (middle) and with regularization α = 0.05 (right).

24 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Visual results: noiseless case

0.05 0.1 Excitation spectra 0.02 0.04 Emission spectra 50 100 Concentrations 0.02 0.04 0.02 0.04 100 200 0.05 0.01 0.02 50 100 0.02 0.04 0.01 0.02 100 200 0.02 0.04 0.05 100 200 300 350 400 450 500 0.05 λem 300 350 400 450 500 0.02 0.04 λex 50 100 5 10 Experiments

Figure: R = 6 - reference spectra / BC-VMFB without penalty / BC-VMFB with penalty α = 0.05.

24 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Visual results: noisy case

Reference

λex 400 500 λex 400 500 λex 400 500 λex 400 500 λex 400 500 λem λex 350 400 450 500 400 500

BC−VMFB without penalty

λem 350 400 450 500

BC−VMFB with penalty

λem 350 400 450 500 0.05 0.5 5

Figure: FEEM of reference (left) - FEEM reconstructed using BC-VMFB without regularization (middle) and with regularization α = 0.05 (right).

25 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Visual results: noisy case

0.05 0.1 Excitation spectra 0.02 0.04 Emission spectra 50 100 Concentrations 0.02 0.04 0.02 0.04 100 200 0.05 0.01 0.02 50 100 0.02 0.04 0.01 0.02 100 200 0.02 0.04 0.05 100 200 300 350 400 450 500 0.01 0.02 λem 300 350 400 450 500 0.01 0.02 λex 50 100 10 20 Experiments

Figure: R = 6 - reference spectra / BC-VMFB without penalty / BC-VMFB with penalty α = 0.05.

25 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Computer simulation: real experimental data - water monitoring to detect pollutants

◮ Data were acquired automatically every 3 minutes, during a 10 days

monitoring campaign performed on water extracted from an urban river ⇒ tensor of size 36 × 111 × 2594.

◮ The excitation wavelengths range from 225nm to 400nm with a 5nm

bandwidth, whereas the emission wavelengths range from 280nm to 500nm with a 2nm bandwidth.

◮ The FEEM have been pre-processed using the Zepp’s method (negative

values were set to 0). During this experiment, a contamination with diesel oil appeared 7 days after the beginning of the monitoring. Contamination

26 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Results: what about the rank ?

(1)

λex 300 350 400 450 500 250 300 350 400

(2)

300 350 400 450 500 250 300 350 400

(3)

λem λex 300 350 400 450 500 250 300 350 400

Estimated FEEM

(4)

λem 300 350 400 450 500 250 300 350 400 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10

(1)

λex 300 350 400 450 500 250 300 350 400

(2)

300 350 400 450 500 250 300 350 400

(3)

λem λex 300 350 400 450 500 250 300 350 400

Estimated FEEM

(4)

λem 300 350 400 450 500 250 300 350 400 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10

penalized BC-VMFB algorithm Bro’s N-way algorithm Case R = 4

27 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Results: what about the rank ?

(1)

λex 300 350 400 450 500 250 300 350 400

(2)

300 350 400 450 500 250 300 350 400

(3)

λex 300 350 400 450 500 250 300 350 400

(4)

300 350 400 450 500 250 300 350 400

(5)

λem λex 300 350 400 450 500 250 300 350 400

Estimated FEEM

(6)

λem 300 350 400 450 500 250 300 350 400 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10

(1)

λex 300 350 400 450 500 250 300 350 400

(2)

300 350 400 450 500 250 300 350 400

(3)

λex 300 350 400 450 500 250 300 350 400

(4)

300 350 400 450 500 250 300 350 400

(5)

λem λex 300 350 400 450 500 250 300 350 400

Estimated FEEM

(6)

λem 300 350 400 450 500 250 300 350 400 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 10

penalized BC-VMFB algorithm Bro’s N-way algorithm Case R = 6

27 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Results: concentrations

200 400 600 800 1000 1200 1400 2 4 6 8 10 12 14 x 10

4

Experiment

Normalized concentrations

(1) (2) (3) (4) 200 400 600 800 1000 1200 1400 2 4 6 8 10 12 14 x 10

4

Experiment

Normalized concentrations

(1) (2) (3) (4)

penalized BC-VMFB algorithm Bro’s N-way algorithm Case R = 4

28 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Concentrations estimated by BC-VMFB

200 400 600 800 1000 1200 1400 2 4 6 8 10 12 14 x 10

4

Experiment

Normalized concentrations

(1) (2) (3) (4)

Case R = 4

29 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Concentrations estimated by BC-VMFB

200 400 2

29 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

200 400 600 800 1000 1200 1400 2 4 6 8 10 12 14 x 10

4

Experiment

Normalized concentrations

(1) (2) (3) (4) (5) (6) 200 400 600 800 1000 1200 1400 2 4 6 8 10 12 14 x 10

4

Experiment

Normalized concentrations

(1) (2) (3) (4) (5) (6)

penalized BC-VMFB algorithm Bro’s N-way algorithm Case R = 6

30 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Concentrations estimated by BC-VMFB

200 400 600 800 1000 1200 1400 2 4 6 8 10 12 14 x 10

4

Experiment

Normalized concentrations

(1) (2) (3) (4) (5) (6)

Case R = 6

31 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Concentrations estimated by BC-VMFB

200 400 2

31 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Conclusion

◮ clear theoretical and mathematical framework for CPD decomposition; ◮ interesting properties of the proposed approach: reliability, robustness

versus noise and overestimation of the rank, good performance despite model errors and relative quickness;

◮ promising results on simulated and real data.

Perspectives:

◮ extension to higher order tensor (order N; LVA-ICA Grenoble 21-23

• Feb. 2017);

◮ possibility of considering missing data; ◮ study other preconditionning stategies.

32 / 32

Introduction Proximal tools Application to CPD Numerical simulations Conclusion and future work

Conclusion

◮ clear theoretical and mathematical framework for CPD decomposition; ◮ interesting properties of the proposed approach: reliability, robustness

versus noise and overestimation of the rank, good performance despite model errors and relative quickness;

◮ promising results on simulated and real data.

Perspectives:

◮ extension to higher order tensor (order N; LVA-ICA Grenoble 21-23

• Feb. 2017);

◮ possibility of considering missing data; ◮ study other preconditionning stategies.

Questions ? Thank you ! ?

32 / 32