Data Unfolding with Wiener-SVD Method arXiv:1705.03568 Tang, a,1 X. - - PowerPoint PPT Presentation

▶
data unfolding with wiener svd method
SMART_READER_LITE
LIVE PREVIEW

Data Unfolding with Wiener-SVD Method arXiv:1705.03568 Tang, a,1 X. - - PowerPoint PPT Presentation

Feb 27, 2023 259 likes •722 views

Data Unfolding with Wiener-SVD Method arXiv:1705.03568 Tang, a,1 X. Li, b,1 X. Qian, a,2 H. Wei, a C. Zhang a W. a Physics Department, Brookhaven National Laboratory, Upton, NY, USA b State University of New Y ork at Stony Brook, Department of

slide-1
SLIDE 1

Data Unfolding with Wiener-SVD Method

W. Tang,a,1 X. Li,b,1 X. Qian,a,2 H. Wei,a C. Zhanga

a Physics Department, Brookhaven National Laboratory, Upton, NY, USA b State University of New

Y

  • rk at Stony Brook, Department of Physics and Astronomy, Stony Brook, NY, USA

Xiaoyue Li DPF2017 - FNAL 1

xiaoyue.li@stonybrook.edu DPF – FNAL, Aug. 2nd, 2017 arXiv:1705.03568

slide-2
SLIDE 2

Outline

§ The unfolding problem § Wiener filter in digital signal processing § SVD unfolding with Wiener filter § Data unfolding example

§ Cross-section § Reactor neutrino flux

§ Summary and recommendation

Xiaoyue Li DPF2017 - FNAL 2

slide-3
SLIDE 3

The data unfolding problem

§ The problem: 𝑔

"#$%&(𝑦) = +𝑆 𝑦 𝑧 𝑔 .&/$ 𝑧 d𝑧

§ 𝑆 𝑦 𝑧 : response function § 𝑛3 = ∑ 𝑆35𝑡5

5

§ 𝑆35 = 𝑄 observed in bin 𝑗 true value in bin 𝑘)

§ Why not do 𝒕 = 𝑆GH I m ? § Due to statistical fluctuation, the unfolded spectrum from direct inversion of response matrix bares no resemblance to the true spectrum. § Introduce additional constraints.

§ E.g. trade bias for smoothness § Bayesian analysis

Xiaoyue Li DPF2017 - FNAL 3

2 4 6 8 10 12 14 16 18 20 200 400 600 800 1000 1200

True distribution

2 4 6 8 10 12 14 16 18 20 200 400 600 800 1000 1200

Smeared distribution Data

slide-4
SLIDE 4

The data unfolding problem

§ The problem: 𝑔

"#$%&(𝑦) = +𝑆 𝑦 𝑧 𝑔 .&/$ 𝑧 d𝑧

§ 𝑆 𝑦 𝑧 : response function § 𝑛3 = ∑ 𝑆35𝑡5

5

§ 𝑆35 = 𝑄 observed in bin 𝑗 true value in bin 𝑘)

§ Why not do 𝒕 = 𝑆GH I m ? § Due to statistical fluctuation, the unfolded spectrum from direct inversion of response matrix bares no resemblance to the true spectrum. § Introduce additional constraints.

§ E.g. trade bias for smoothness § Bayesian analysis

Xiaoyue Li DPF2017 - FNAL 4

2 4 6 8 10 12 14 16 18 20 200 400 600 800 1000 1200

True distribution

2 4 6 8 10 12 14 16 18 20 200 400 600 800 1000 1200

Smeared distribution Data

slide-5
SLIDE 5

The data unfolding problem

§ The problem: 𝑔

"#$%&(𝑦) = +𝑆 𝑦 𝑧 𝑔 .&/$ 𝑧 d𝑧

§ 𝑆 𝑦 𝑧 : response function § 𝑛3 = ∑ 𝑆35𝑡5

5

§ 𝑆35 = 𝑄 observed in bin 𝑗 true value in bin 𝑘)

§ Why not do 𝒕 = 𝑆GH I m ? § Due to statistical fluctuation, the unfolded spectrum from direct inversion of response matrix bares no resemblance to the true spectrum. § Introduce additional constraints.

§ E.g. trade bias for smoothness § Bayesian analysis

Xiaoyue Li DPF2017 - FNAL 5

2 4 6 8 10 12 14 16 18 20 200 400 600 800 1000 1200

True distribution

2 4 6 8 10 12 14 16 18 20 200 400 600 800 1000 1200

Smeared distribution Data

2 4 6 8 10 12 14 16 18 20 500 1000 1500 2000

True distribution Unfolded Data

2 4 6 8 10 12 14 16 18 20 500 1000 1500 2000

True distribution Unfolded Data

slide-6
SLIDE 6

The data unfolding problem

§ The problem: 𝑔

"#$%&(𝑦) = +𝑆 𝑦 𝑧 𝑔 .&/$ 𝑧 d𝑧

§ 𝑆 𝑦 𝑧 : response function § 𝑛3 = ∑ 𝑆35𝑡5

5

§ 𝑆35 = 𝑄 observed in bin 𝑗 true value in bin 𝑘)

§ Why not do 𝒕 = 𝑆GH I m ? § Due to statistical fluctuation, the unfolded spectrum from direct inversion of response matrix bares no resemblance to the true spectrum. § Introduce additional constraints.

§ E.g. trade bias for smoothness § Bayesian analysis

Xiaoyue Li DPF2017 - FNAL 6

2 4 6 8 10 12 14 16 18 20 200 400 600 800 1000 1200

True distribution

2 4 6 8 10 12 14 16 18 20 200 400 600 800 1000 1200

Smeared distribution Data

2 4 6 8 10 12 14 16 18 20 500 1000 1500 2000

True distribution Unfolded Data

2 4 6 8 10 12 14 16 18 20 500 1000 1500 2000

True distribution Unfolded Data

slide-7
SLIDE 7

The data unfolding problem

§ The problem: 𝑔

"#$%&(𝑦) = +𝑆 𝑦 𝑧 𝑔 .&/$ 𝑧 d𝑧

§ 𝑆 𝑦 𝑧 : response function § 𝑛3 = ∑ 𝑆35𝑡5

5

§ 𝑆35 = 𝑄 observed in bin 𝑗 true value in bin 𝑘)

§ Why not do 𝒕 = 𝑆GH I m ? § Due to statistical fluctuation, the unfolded spectrum from direct inversion of response matrix bares no resemblance to the true spectrum. § Introduce additional constraints.

§ E.g. trade bias for smoothness § Bayesian analysis

Xiaoyue Li DPF2017 - FNAL 7

2 4 6 8 10 12 14 16 18 20 200 400 600 800 1000 1200

True distribution

2 4 6 8 10 12 14 16 18 20 200 400 600 800 1000 1200

Smeared distribution Data

2 4 6 8 10 12 14 16 18 20 500 1000 1500 2000

True distribution Unfolded Data

2 4 6 8 10 12 14 16 18 20 500 1000 1500 2000

True distribution Unfolded Data

slide-8
SLIDE 8

Wiener filter in digital signal processing (1)

§ Deconvolution problem:

§ 𝑁 𝑢M = ∫ 𝑆 𝑢, 𝑢M I 𝑇 𝑢 d𝑢

Q GQ

§ True signal 𝑇(𝑢), measured signal 𝑁(𝑢M), response function 𝑆(𝑢, 𝑢′) ≡ 𝑆(𝑢 − 𝑢′)

§ Fourier transform: 𝑁 𝜕 = 𝑆(𝜕) I 𝑇(𝜕) → 𝑇 𝜕 = 𝑁 𝜕 /𝑆(𝜕) → Inverse FFT 𝑇 𝜕 → 𝑇 𝑢 § The response function 𝑆(𝜕) does not address noise contributions to the measured signal. Worse still, 𝑆(𝜕) is generally smaller at higher frequencies due to the shaping features of electronics, resulting in amplification of noises.

Xiaoyue Li DPF2017 - FNAL 8

s) µ Time ( 20 40 60 80 100 V (mV) 0.05 0.1 0.15

  • 3

10 ×

Single Electron Response

s) µ Time ( 20 40 60 80 100

  • e

1000 2000 3000 4000

Signal: 200k electrons

Response function 𝑆(𝑢, 𝑢′) True signal 𝑇(𝑢) Measured signal, with statistical fluctuation

s) µ Time ( 20 40 60 80 100 ADC 10 20 30

Simulated Measured Signal

slide-9
SLIDE 9

Wiener filter in digital signal processing (1)

§ Deconvolution problem:

§ 𝑁 𝑢M = ∫ 𝑆 𝑢, 𝑢M I 𝑇 𝑢 d𝑢

Q GQ

§ True signal 𝑇(𝑢), measured signal 𝑁(𝑢M), response function 𝑆(𝑢, 𝑢′) ≡ 𝑆(𝑢 − 𝑢′)

§ Fourier transform: 𝑁 𝜕 = 𝑆(𝜕) I 𝑇(𝜕) → 𝑇 𝜕 = 𝑁 𝜕 /𝑆(𝜕) → Inverse FFT 𝑇 𝜕 → 𝑇 𝑢 § The response function 𝑆(𝜕) does not address noise contributions to the measured signal. Worse still, 𝑆(𝜕) is generally smaller at higher frequencies due to the shaping features of electronics, resulting in amplification of noises.

Xiaoyue Li DPF2017 - FNAL 9

s) µ Time ( 20 40 60 80 100 V (mV) 0.05 0.1 0.15

  • 3

10 ×

Single Electron Response

s) µ Time ( 20 40 60 80 100

  • e

1000 2000 3000 4000

Signal: 200k electrons

Response function 𝑆(𝑢, 𝑢′) True signal 𝑇(𝑢) Measured signal, with statistical fluctuation

s) µ Time ( 20 40 60 80 100 ADC 10 20 30

Simulated Measured Signal (MHz) ω 0.2 0.4 0.6 0.8 1 ) ω R( 10 20 30

Response function in Frequency domain

Response function 𝑆(𝜕)

(MHz) ω 0.2 0.4 0.6 0.8 1 ) ω M( 500 1000 1500 2000

Data in Frequency domain

Measured signal in frequency domain 𝑁 𝜕

slide-10
SLIDE 10

Wiener filter in digital signal processing (1)

§ Deconvolution problem:

§ 𝑁 𝑢M = ∫ 𝑆 𝑢, 𝑢M I 𝑇 𝑢 d𝑢

Q GQ

§ True signal 𝑇(𝑢), measured signal 𝑁(𝑢M), response function 𝑆(𝑢, 𝑢′) ≡ 𝑆(𝑢 − 𝑢′)

§ Fourier transform: 𝑁 𝜕 = 𝑆(𝜕) I 𝑇(𝜕) → 𝑇 𝜕 = 𝑁 𝜕 /𝑆(𝜕) → Inverse FFT 𝑇 𝜕 → 𝑇 𝑢 § The response function 𝑆(𝜕) does not address noise contributions to the measured signal. Worse still, 𝑆(𝜕) is generally smaller at higher frequencies due to the shaping features of electronics, resulting in amplification of noises.

Xiaoyue Li DPF2017 - FNAL 10

s) µ Time ( 20 40 60 80 100 V (mV) 0.05 0.1 0.15

  • 3

10 ×

Single Electron Response

s) µ Time ( 20 40 60 80 100

  • e

1000 2000 3000 4000

Signal: 200k electrons

Response function 𝑆(𝑢, 𝑢′) True signal 𝑇(𝑢) Measured signal, with statistical fluctuation

s) µ Time ( 20 40 60 80 100 ADC 10 20 30

Simulated Measured Signal (MHz) ω 0.2 0.4 0.6 0.8 1 ) ω R( 10 20 30

Response function in Frequency domain

Response function 𝑆(𝜕)

(MHz) ω 0.2 0.4 0.6 0.8 1 ) ω M( 500 1000 1500 2000

Data in Frequency domain

Measured signal in frequency domain 𝑁 𝜕

s) µ Time ( 20 40 60 80 100

  • e
  • 10000

10000

6

10 ×

Deconvoluton without filter

Deconvoluted signal 𝑇(𝑢)

slide-11
SLIDE 11

Wiener filter in digital signal processing (1)

§ Deconvolution problem:

§ 𝑁 𝑢M = ∫ 𝑆 𝑢, 𝑢M I 𝑇 𝑢 d𝑢

Q GQ

§ True signal 𝑇(𝑢), measured signal 𝑁(𝑢M), response function 𝑆(𝑢, 𝑢′) ≡ 𝑆(𝑢 − 𝑢′)

§ Fourier transform: 𝑁 𝜕 = 𝑆(𝜕) I 𝑇(𝜕) → 𝑇 𝜕 = 𝑁 𝜕 /𝑆(𝜕) → Inverse FFT 𝑇 𝜕 → 𝑇 𝑢 § The response function 𝑆(𝜕) does not address noise contributions to the measured signal. Worse still, 𝑆(𝜕) is generally smaller at higher frequencies due to the shaping features of electronics, resulting in amplification of noises.

Xiaoyue Li DPF2017 - FNAL 11

s) µ Time ( 20 40 60 80 100 V (mV) 0.05 0.1 0.15

  • 3

10 ×

Single Electron Response

s) µ Time ( 20 40 60 80 100

  • e

1000 2000 3000 4000

Signal: 200k electrons

Response function 𝑆(𝑢, 𝑢′) True signal 𝑇(𝑢) Measured signal, with statistical fluctuation

s) µ Time ( 20 40 60 80 100 ADC 10 20 30

Simulated Measured Signal (MHz) ω 0.2 0.4 0.6 0.8 1 ) ω R( 10 20 30

Response function in Frequency domain

Response function 𝑆(𝜕)

(MHz) ω 0.2 0.4 0.6 0.8 1 ) ω M( 500 1000 1500 2000

Data in Frequency domain

Measured signal in frequency domain 𝑁 𝜕

s) µ Time ( 20 40 60 80 100

  • e
  • 10000

10000

6

10 ×

Deconvoluton without filter

Deconvoluted signal 𝑇(𝑢)

slide-12
SLIDE 12

Wiener filter in digital signal processing (2)

§ To address the issue with noise:

§ 𝑇 X 𝜕 =

Y Z [(Z) I 𝐺(𝜕)

§ Wiener filter:

§ 𝐺 𝜕 =

]^(Z) ]^(Z)_`^(Z)

§ The functional form of Wiener filter is obtained by minimizing 𝐹 𝐺 𝜕 I 𝑁 𝜕 − 𝑇(𝜕) b = 𝐹 𝐺 𝜕 I (𝑇 𝜕 + 𝑂 𝜕 ) − 𝑇(𝜕) b § Wiener filter is designed to achieve the best signal/noise.

Xiaoyue Li DPF2017 - FNAL 12

Noise filter

(MHz) ω 0.2 0.4 0.6 0.8 1 ) ω F( 0.2 0.4 0.6 0.8 1

Wiener Filter @ Frequency Domain

s) µ Time (

  • 40
  • 20

20 40 F(t) 0.01 0.02 0.03

Wiener Filter @ Time Domain

Wiener filter in frequency domain Wiener filter in time domain

Noise

slide-13
SLIDE 13

Wiener filter in digital signal processing (3)

Xiaoyue Li DPF2017 - FNAL 13

s) µ Time ( 20 40 60 80 100

  • e
  • 10000

10000

6

10 ×

Deconvoluton without filter

s) µ Time ( 20 40 60 80 100

  • e

1000 2000 3000

Deconvoluted result Truth

Deconvoluton with Wiener filter

Deconvoluted result

80

Truth

Deconvolution without filter Deconvolution with Wiener filter

slide-14
SLIDE 14

Reminder of SVD Method

§ Minimize 𝜓b: 𝜓b(𝑡) = 𝐧 − 𝐬 I 𝑡 h I 𝐷𝑝𝑤GH I 𝐧 − 𝐬 I 𝑡

§ Cholesky decomposition: 𝐷𝑝𝑤GH = 𝑅h I 𝑅 § Pre-scaling: 𝑁 ≔ 𝑅 I 𝐧, 𝑆 ≔ 𝑅 I 𝐬 § ML estimator: 𝑡̂ = 𝑆h𝑆 GH I 𝑆h I 𝑁

§ Equivalently, the estimator can be written as 𝑡̂ = 𝑆h𝑆 GH I 𝑆h I (𝑆 I 𝑡opqr + 𝑂) where 𝑂, the “noise” coming from uncertainties, follows a normal distribution.

Xiaoyue Li DPF2017 - FNAL 14

slide-15
SLIDE 15

Reminder of SVD Method

§ Minimize 𝜓b: 𝜓b(𝑡) = 𝐧 − 𝐬 I 𝑡 h I 𝐷𝑝𝑤GH I 𝐧 − 𝐬 I 𝑡

§ Cholesky decomposition: 𝐷𝑝𝑤GH = 𝑅h I 𝑅 § Pre-scaling: 𝑁 ≔ 𝑅 I 𝐧, 𝑆 ≔ 𝑅 I 𝐬 § ML estimator: 𝑡̂ = 𝑆h𝑆 GH I 𝑆h I 𝑁

§ Equivalently, the estimator can be written as 𝑡̂ = 𝑆h𝑆 GH I 𝑆h I (𝑆 I 𝑡opqr + 𝑂) where 𝑂, the “noise” coming from uncertainties, follows a normal distribution.

Xiaoyue Li DPF2017 - FNAL 15

Singular value decomposition (SVD) 𝑆 = 𝑉 I 𝐸 I 𝑊h 𝑉, 𝑊: orthogonal 𝐸: diagonal (diagonal elements ranked by magnitude)

slide-16
SLIDE 16

Reminder of SVD Method

§ Minimize 𝜓b: 𝜓b(𝑡) = 𝐧 − 𝐬 I 𝑡 h I 𝐷𝑝𝑤GH I 𝐧 − 𝐬 I 𝑡

§ Cholesky decomposition: 𝐷𝑝𝑤GH = 𝑅h I 𝑅 § Pre-scaling: 𝑁 ≔ 𝑅 I 𝐧, 𝑆 ≔ 𝑅 I 𝐬 § ML estimator: 𝑡̂ = 𝑆h𝑆 GH I 𝑆h I 𝑁

§ Equivalently, the estimator can be written as 𝑡̂ = 𝑆h𝑆 GH I 𝑆h I (𝑆 I 𝑡opqr + 𝑂) where 𝑂, the “noise” coming from uncertainties, follows a normal distribution.

Xiaoyue Li DPF2017 - FNAL 16

𝑡̂ = 𝑊 I 𝐸GH I 𝑉h I 𝑆 I 𝑡.&/$ + 𝑂 = 𝑊 I 𝐸GH I (𝑆v I 𝑡.&/$ + 𝑂v) = 𝑊 I 𝐸GH I 𝑁v

§ Because 𝑉 is an orthogonal matrix, 𝑂v are still uncorrelated and follows a normal distribution.

§ 𝑁v = 𝑉h I 𝑁 is the measurement in the effective frequency domain

𝑂v=𝑉h I 𝑂 𝑁v=𝑉h I 𝑁

Singular value decomposition (SVD) 𝑆 = 𝑉 I 𝐸 I 𝑊h 𝑉, 𝑊: orthogonal 𝐸: diagonal (diagonal elements ranked by magnitude)

slide-17
SLIDE 17

Reminder of SVD Method

§ Minimize 𝜓b: 𝜓b(𝑡) = 𝐧 − 𝐬 I 𝑡 h I 𝐷𝑝𝑤GH I 𝐧 − 𝐬 I 𝑡

§ Cholesky decomposition: 𝐷𝑝𝑤GH = 𝑅h I 𝑅 § Pre-scaling: 𝑁 ≔ 𝑅 I 𝐧, 𝑆 ≔ 𝑅 I 𝐬 § ML estimator: 𝑡̂ = 𝑆h𝑆 GH I 𝑆h I 𝑁

§ Equivalently, the estimator can be written as 𝑡̂ = 𝑆h𝑆 GH I 𝑆h I (𝑆 I 𝑡opqr + 𝑂) where 𝑂, the “noise” coming from uncertainties, follows a normal distribution.

Xiaoyue Li DPF2017 - FNAL 17

𝑡̂ = 𝑊 I 𝐸GH I 𝑉h I 𝑆 I 𝑡.&/$ + 𝑂 = 𝑊 I 𝐸GH I (𝑆v I 𝑡.&/$ + 𝑂v) = 𝑊 I 𝐸GH I 𝑁v

§ Because 𝑉 is an orthogonal matrix, 𝑂v are still uncorrelated and follows a normal distribution.

§ 𝑁v = 𝑉h I 𝑁 is the measurement in the effective frequency domain

Recall in signal processing 𝑇 𝜕 = [𝑆(𝜕)]GH𝑁 𝜕

𝑂v=𝑉h I 𝑂 𝑁v=𝑉h I 𝑁

Singular value decomposition (SVD) 𝑆 = 𝑉 I 𝐸 I 𝑊h 𝑉, 𝑊: orthogonal 𝐸: diagonal (diagonal elements ranked by magnitude)

slide-18
SLIDE 18

Reminder of Regularization unfolding

§ The estimator can be obtained by finding the maximum of 𝜚 𝑡 = log 𝑀 𝑡 + 𝜐Σ(𝑡) § E.g. Tikhonov regularization Σ~ 𝑡 𝐹 = − + d~𝑡(𝐹) d~𝐹

b

d𝐹,𝑙 = 0, 1, 2, … § When 𝑙 = 0, the unfolded result can be written as 𝑡̂ = 𝑆h𝑆 + 𝜐𝐽 GH I 𝑆h I 𝑁 where 𝐽 is an identity matrix and 𝜐 is the regularization strength

Xiaoyue Li DPF2017 - FNAL 18

Likelihood function Regularization strength Regularization function

slide-19
SLIDE 19

Reminder of Regularization unfolding

§ The estimator can be obtained by finding the maximum of 𝜚 𝑡 = log 𝑀 𝑡 + 𝜐Σ(𝑡) § E.g. Tikhonov regularization Σ~ 𝑡 𝐹 = − + d~𝑡(𝐹) d~𝐹

b

d𝐹,𝑙 = 0, 1, 2, … § When 𝑙 = 0, the unfolded result can be written as 𝑡̂ = 𝑆h𝑆 + 𝜐𝐽 GH I 𝑆h I 𝑁 where 𝐽 is an identity matrix and 𝜐 is the regularization strength

Xiaoyue Li DPF2017 - FNAL 19

Regularization strength Regularization function

𝑆 = 𝑉 I 𝐸 I 𝑊h

Likelihood function

slide-20
SLIDE 20

Reminder of Regularization unfolding

§ The estimator can be obtained by finding the maximum of 𝜚 𝑡 = log 𝑀 𝑡 + 𝜐Σ(𝑡) § E.g. Tikhonov regularization Σ~ 𝑡 𝐹 = − + d~𝑡(𝐹) d~𝐹

b

d𝐹,𝑙 = 0, 1, 2, … § When 𝑙 = 0, the unfolded result can be written as 𝑡̂ = 𝑆h𝑆 + 𝜐𝐽 GH I 𝑆h I 𝑁 where 𝐽 is a unit matrix and 𝜐 is the regularization strength

Xiaoyue Li DPF2017 - FNAL 20

Regularization strength Regularization function

𝑆 = 𝑉 I 𝐸 I 𝑊h 𝑡̂ = 𝐵 I 𝑆h𝑆 GH I 𝑆h I 𝑁 where 𝐵 = 𝑊 I 𝐺 I 𝑊h, 𝐺35 =

†‡

^

†‡

^_ˆ𝜀35 (𝑒3 is the diagonal terms of 𝐸)

Likelihood function

slide-21
SLIDE 21

Reminder of Regularization unfolding

§ The estimator can be obtained by finding the maximum of 𝜚 𝑡 = log 𝑀 𝑡 + 𝜐Σ(𝑡) § E.g. Tikhonov regularization Σ~ 𝑡 𝐹 = − + d~𝑡(𝐹) d~𝐹

b

d𝐹,𝑙 = 0, 1, 2, … § When 𝑙 = 0, the unfolded result can be written as 𝑡̂ = 𝑆h𝑆 + 𝜐𝐽 GH I 𝑆h I 𝑁 where 𝐽 is a unit matrix and 𝜐 is the regularization strength

Xiaoyue Li DPF2017 - FNAL 21

Regularization strength Regularization function

𝑆 = 𝑉 I 𝐸 I 𝑊h 𝑡̂ = 𝐵 I 𝑆h𝑆 GH I 𝑆h I 𝑁 where 𝐵 = 𝑊 I 𝐺 I 𝑊h, 𝐺35 =

†‡

^

†‡

^_ˆ𝜀35 (𝑒3 is the diagonal terms of 𝐸)

Regularization is effectively introducing an additional smearing to the unfolded results

Likelihood function

slide-22
SLIDE 22

Reminder of Regularization unfolding

§ The estimator can be obtained by finding the maximum of 𝜚 𝑡 = log 𝑀 𝑡 + 𝜐Σ(𝑡) § E.g. Tikhonov regularization Σ~ 𝑡 𝐹 = − + d~𝑡(𝐹) d~𝐹

b

d𝐹,𝑙 = 0, 1, 2, … § When 𝑙 = 0, the unfolded result can be written as 𝑡̂ = 𝑆h𝑆 + 𝜐𝐽 GH I 𝑆h I 𝑁 where 𝐽 is a unit matrix and 𝜐 is the regularization strength

Xiaoyue Li DPF2017 - FNAL 22

Regularization strength Regularization function

𝑆 = 𝑉 I 𝐸 I 𝑊h 𝑡̂ = 𝐵 I 𝑆h𝑆 GH I 𝑆h I 𝑁 where 𝐵 = 𝑊 I 𝐺 I 𝑊h, 𝐺35 =

†‡

^

†‡

^_ˆ𝜀35 (𝑒3 is the diagonal terms of 𝐸)

𝑡̂ = 𝑊 I 𝐺 I 𝐸GH I 𝑆v I 𝑡opqr + 𝑂v , where 𝐺35 = 𝑒3

b

𝑒3

b + 𝜐 𝜀35

𝑂v are suppressed by 𝑒3/(𝑒3

b + 𝜐), regardless of signal strength

Likelihood function

slide-23
SLIDE 23

§ Recall: 𝑡̂ = 𝑊 I 𝐺 I 𝐸GH I (𝑆v I 𝑡opqr + 𝑂v),

Xiaoyue Li DPF2017 - FNAL 23

𝐺

35 =

𝑒3

b

𝑒3

b + 𝜐 𝜀35

slide-24
SLIDE 24

Wiener-SVD

§ Recall: 𝑡̂ = 𝑊 I 𝐺 I 𝐸GH I (𝑆v I 𝑡opqr + 𝑂v),

Xiaoyue Li DPF2017 - FNAL 24

𝐺

35 =

𝑒3

b

𝑒3

b + 𝜐 𝜀35

W W: Wiener filter

slide-25
SLIDE 25

Wiener-SVD

§ Recall: 𝑡̂ = 𝑊 I 𝐺 I 𝐸GH I (𝑆v I 𝑡opqr + 𝑂v), § In Wiener-SVD method, the functional form

  • f W (replace 𝐺 as filter) also factors in the

signal expectation in the effective frequency domain 𝑁v ≔ 𝑉h I 𝑁 Ž = 𝑉h I 𝑆 I 𝑡̅ = 𝐸 I 𝑊h I 𝑡̅

Xiaoyue Li DPF2017 - FNAL 25

𝐺

35 =

𝑒3

b

𝑒3

b + 𝜐 𝜀35

W W: Wiener filter 𝑋

35 =

𝑁v3

b

𝑁v3

b + 𝑂v3 b 𝜀35

slide-26
SLIDE 26

Wiener-SVD

§ Recall: 𝑡̂ = 𝑊 I 𝐺 I 𝐸GH I (𝑆v I 𝑡opqr + 𝑂v), § In Wiener-SVD method, the functional form

  • f W (replace 𝐺 as filter) also factors in the

signal expectation in the effective frequency domain 𝑁v ≔ 𝑉h I 𝑁 Ž = 𝑉h I 𝑆 I 𝑡̅ = 𝐸 I 𝑊h I 𝑡̅ § Signal: 𝑁v3

b = 𝑒3 b I ∑ 𝑊 35 h I 𝑡5

Ž

5 b

§ Noise: 1 (normal distribution)

Xiaoyue Li DPF2017 - FNAL 26

𝐺

35 =

𝑒3

b

𝑒3

b + 𝜐 𝜀35

W W: Wiener filter 𝑋

35 =

𝑁v3

b

𝑁v3

b + 𝑂v3 b 𝜀35

slide-27
SLIDE 27

Wiener-SVD

§ Recall: 𝑡̂ = 𝑊 I 𝐺 I 𝐸GH I (𝑆v I 𝑡opqr + 𝑂v), § In Wiener-SVD method, the functional form

  • f W (replace 𝐺 as filter) also factors in the

signal expectation in the effective frequency domain 𝑁v ≔ 𝑉h I 𝑁 Ž = 𝑉h I 𝑆 I 𝑡̅ = 𝐸 I 𝑊h I 𝑡̅ § Signal: 𝑁v3

b = 𝑒3 b I ∑ 𝑊 35 h I 𝑡5

Ž

5 b

§ Noise: 1 (normal distribution)

Xiaoyue Li DPF2017 - FNAL 27

𝐺

35 =

𝑒3

b

𝑒3

b + 𝜐 𝜀35

W W: Wiener filter 𝑋

35 =

𝑁v3

b

𝑁v3

b + 𝑂v3 b 𝜀35

𝑋

35 =

𝑒3

b I ∑ 𝑊 35 h I 𝑡5

Ž

5 b

𝑒3

b I ∑ 𝑊 35 h I 𝑡5

Ž

5 b + 1

𝜀35

slide-28
SLIDE 28

§ In Tikhonov regularization, when 𝑙 ≠ 0, the solution is 𝑡̂ = 𝑆h𝑆 + 𝜐 𝐷h𝐷 GH I 𝑆h I 𝑁 Ž where in the case of 𝑙 = 2, 𝐷 = −1 1 1 −2 1 1 −2 ⋯ ⋯ ⋯ ⋮ ⋮ ⋮ ⋱ ⋯ ⋯ −2 1 1 −2 1 −1 1

Xiaoyue Li DPF2017 - FNAL 28

slide-29
SLIDE 29

Generalized Wiener-SVD

§ In Tikhonov regularization, when 𝑙 ≠ 0, the solution is 𝑡̂ = 𝑆h𝑆 + 𝜐 𝐷h𝐷 GH I 𝑆h I 𝑁 Ž where in the case of 𝑙 = 2, 𝐷 = −1 1 1 −2 1 1 −2 ⋯ ⋯ ⋯ ⋮ ⋮ ⋮ ⋱ ⋯ ⋯ −2 1 1 −2 1 −1 1 § Similarly to Tikhonov regularization, we can add an additional matrix C to the equation 𝑁 Ž = 𝑆 I 𝐷GH I (𝐷 I 𝑡̅ )

Xiaoyue Li DPF2017 - FNAL 29

slide-30
SLIDE 30

Generalized Wiener-SVD

§ In Tikhonov regularization, when 𝑙 ≠ 0, the solution is 𝑡̂ = 𝑆h𝑆 + 𝜐 𝐷h𝐷 GH I 𝑆h I 𝑁 Ž where in the case of 𝑙 = 2, 𝐷 = −1 1 1 −2 1 1 −2 ⋯ ⋯ ⋯ ⋮ ⋮ ⋮ ⋱ ⋯ ⋯ −2 1 1 −2 1 −1 1 § Similarly to Tikhonov regularization, we can add an additional matrix C to the equation 𝑁 Ž = 𝑆 I 𝐷GH I (𝐷 I 𝑡̅ )

Xiaoyue Li DPF2017 - FNAL 30

Since the effective frequency domain is determined by the smearing matrix 𝑆, the inclusion of 𝐷 would alter the basis of the effective frequency domain.

slide-31
SLIDE 31

Generalized Wiener-SVD

§ In Tikhonov regularization, when 𝑙 ≠ 0, the solution is 𝑡̂ = 𝑆h𝑆 + 𝜐 𝐷h𝐷 GH I 𝑆h I 𝑁 Ž where in the case of 𝑙 = 2, 𝐷 = −1 1 1 −2 1 1 −2 ⋯ ⋯ ⋯ ⋮ ⋮ ⋮ ⋱ ⋯ ⋯ −2 1 1 −2 1 −1 1 § Similarly to Tikhonov regularization, we can add an additional matrix C to the equation 𝑁 Ž = 𝑆 I 𝐷GH I (𝐷 I 𝑡̅ )

Xiaoyue Li DPF2017 - FNAL 31

Since the effective frequency domain is determined by the smearing matrix 𝑆, the inclusion of 𝐷 would alter the basis of the effective frequency domain.

SVD 𝑆𝐷GH = 𝑉• I 𝐸• I 𝑊

  • h
slide-32
SLIDE 32

Generalized Wiener-SVD

§ In Tikhonov regularization, when 𝑙 ≠ 0, the solution is 𝑡̂ = 𝑆h𝑆 + 𝜐 𝐷h𝐷 GH I 𝑆h I 𝑁 Ž where in the case of 𝑙 = 2, 𝐷 = −1 1 1 −2 1 1 −2 ⋯ ⋯ ⋯ ⋮ ⋮ ⋮ ⋱ ⋯ ⋯ −2 1 1 −2 1 −1 1 § Similarly to Tikhonov regularization, we can add an additional matrix C to the equation 𝑁 Ž = 𝑆 I 𝐷GH I (𝐷 I 𝑡̅ )

Xiaoyue Li DPF2017 - FNAL 32

Since the effective frequency domain is determined by the smearing matrix 𝑆, the inclusion of 𝐷 would alter the basis of the effective frequency domain.

SVD 𝑆𝐷GH = 𝑉• I 𝐸• I 𝑊

  • h

The Wiener-SVD solution becomes 𝑡̂ = 𝐷GH I 𝑊

  • I 𝑋
  • I 𝑊
  • h I 𝐷 I 𝑆h𝑆 GH I 𝑆h I 𝑁

Ž

(𝑋•)35= 𝑒3

b I ∑ 𝑊

  • h

35 I ∑ 𝐷5– I 𝑡̅– – 5 b

𝑒3

b I ∑

𝑊

  • h

35 I ∑ 𝐷5– I 𝑡̅– – 5 b + 1

slide-33
SLIDE 33

Regularization interpretation of Wiener-SVD method

§ The Wiener-SVD unfolding method is equivalent to regularization w.r.t. minimizing the signal to noise ratio in the effective frequency domain 𝑁v 3 = ∑ 𝑊

35 h𝑡5 5

: 𝜚 𝑡 = log𝑀(𝑡) + 1 2 — log 𝑁v 3

b

𝑂b

3

= log𝑀 𝑡 + 1 2 — log ∑ 𝑊

35 h𝑡5 5 b

1

3

Xiaoyue Li DPF2017 - FNAL 33

slide-34
SLIDE 34

Application on neutrino cross-section measurements (1)

Xiaoyue Li DPF2017 - FNAL 34

Energy 2 4 6 8 10 12 Number of events 50 100 150 200 250

true spectrum Asimov spectrum measured spectrum

2 4 6 8 10 12 14 2 4 6 8 10 12 14 0.1 0.2 0.3 0.4 0.5 0.6 True energy bin Reconstructed energy bin 2 4 6 8 10 12 14 2 4 6 8 10 12 14 20 40 60 80 100 120 Reconstructed energy bin Reconstructed energy bin Detector response Statistics only uncertainties

Toy experiment

slide-35
SLIDE 35

Application on neutrino cross-section measurements (2)

Xiaoyue Li DPF2017 - FNAL 35

Effective frequency bin 5 10 15 0.5 1 1.5 2

Wiener-SVD w/ C Regularization w/ C

Effective frequency bin 5 10 15 0.5 1 1.5 2

2

Wiener-SVD w/ C

2

Regularization w/ C

Filter shape

Regularization 𝑡̂ = 𝑊 I 𝐺 I 𝐸GH I 𝑁 𝐺35 = 𝑒3

b

𝑒3

b + 𝜐 𝜀35

Wiener-SVD 𝑡̂ = 𝑊 I 𝑋 I 𝐸GH I 𝑁

𝑋

33 =

𝑒3

b I ∑ 𝑊 35 h I 𝑡̅5 5 b

𝑒3

b I ∑ 𝑊 35 h I 𝑡̅ 5 5 b + 1

slide-36
SLIDE 36

Application on neutrino cross-section measurements (3)

Xiaoyue Li DPF2017 - FNAL 36

Energy 2 4 6 8 10 12 Number of events 50 100 150 200 250

true spectrum unfolded spectrum unfolded data

0.2 0.4 0.6 0.8 2 4 6 8 10 12 14 2 4 6 8 10 12 14 True energy bin True energy bin Energy 2 4 6 8 10 12 Number of events 50 100 150 200 250

true spectrum unfolded spectrum unfolded data

0.2 0.4 0.6 0.8 2 4 6 8 10 12 14 2 4 6 8 10 12 14 True energy bin True energy bin

Additional smearing matrix 𝐵• Unfolded covariance Wiener-SVD w/ 𝐷b Tikhonov regularization w/ 𝐷b

𝑡̂ = 𝐵• I 𝑆h𝑆 GH I 𝑆h I 𝑁 Ž

10 − 10 20 30 40 2 4 6 8 10 12 14 2 4 6 8 10 12 14 True energy bin True energy bin 10 − 10 20 30 40 50 2 4 6 8 10 12 14 2 4 6 8 10 12 14 True energy bin True energy bin

slide-37
SLIDE 37

Application on neutrino cross-section measurements (4)

Xiaoyue Li DPF2017 - FNAL 37

τ 0.2 0.4 0.6 0.8 1 MSE

2

10

3

10

4

10

5

10

6

10

Regularization w/ C

1

Regularization w/ C

2

Regularization w/ C Wiener-SVD w/ C

1

Wiener-SVD w/ C

2

Wiener-SVD w/ C

variance

3

10

4

10

2

b

2

10

3

10

4

10

5

10

6

10

7

10

Regularization w/ C

1

Regularization w/ C

2

Regularization w/ C Wiener-SVD w/ C

1

Wiener-SVD w/ C

2

Wiener-SVD w/ C

Regularization strength

slide-38
SLIDE 38

Application on reactor neutrino measurement (1)

Xiaoyue Li DPF2017 - FNAL 38

Energy (MeV) 2 4 6 8 10 12 Counts 1000 2000 3000 4000

Neutrino spectrum Prompt spectrum

Neutrino Energy (MeV)

2 4 6 8 10 12

Prompt Energy (MeV)

2 4 6 8 10 12 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Prompt Energy (MeV) 2 4 6 8 10 12 Prompt Energy (MeV) 2 4 6 8 10 12 1000 2000 3000 4000 5000 6000 7000 Prompt Energy Bin after Pre-scaling 5 10 15 20 25 Counts 5 10 15 20 Neutrino Energy Bin after Pre-scaling 5 10 15 20 25 Prompt Energy Bin after Pre-scaling 5 10 15 20 25 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04

Neutrino spectrum Prompt spectrum

Covariance matrix includes both statistical and systematic uncertainties

Neutrino Energy (MeV)

2 4 6 8 10 12

Prompt Energy (MeV)

2 4 6 8 10 12

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Detector energy response matrix 𝐬 as reported by the Daya Bay collaboration Pre-scaled prompt energy distribution

𝑁 ≔ 𝑅 I 𝐧 Pre-scaled detector energy response matrix 𝑆 ≔ 𝑅 I 𝐬 𝐷𝑝𝑤GH = 𝑅h I 𝑅

slide-39
SLIDE 39

Application on reactor neutrino measurement (2)

Xiaoyue Li DPF2017 - FNAL 39

Neutrino Energy (MeV) 2 4 6 8 10 12 Counts 1000 2000 3000 4000

True

2

Wiener-SVD with C

2

Regularization C

Neutrino Energy (MeV) 2 4 6 8 10 12 Neutrino Energy (MeV) 2 4 6 8 10 12

1000 2000 3000 4000 5000 6000

Neutrino Energy (MeV) 2 4 6 8 10 12 Neutrino Energy (MeV) 2 4 6 8 10 12 0.1 − 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Effective Frequency Bin 5 10 15 20 25 0.2 0.4 0.6 0.8 1

Regularization Wiener-SVD

C

Effective Frequency Bin 5 10 15 20 25 0.2 0.4 0.6 0.8 1

Regularization Wiener-SVD 2

C

Additional smearing matrix

“Filter” shape in the effective frequency domain w/ 𝐷˜ w/ 𝐷b § In practice, the “true” energy spectrum is unknown; all we have are various models à take average § The 𝜓~

b for model k is calculated by comparing the

predicted spectrum 𝒏~ ≔ 𝒔 I 𝒕~ and measurement 𝒏

Neutrino Energy (MeV) 2 4 6 8 10 12 Neutrino Energy (MeV) 2 4 6 8 10 12

1000 2000 3000 4000 5000 6000

Unfolded covariance matrix

Unfolded distributions

slide-40
SLIDE 40

Application on reactor neutrino measurement (3)

Xiaoyue Li DPF2017 - FNAL 40

2

σ

5

10

6

10

2

b

3

10

4

10

5

10

6

10

7

10

Regularization Wiener-SVD Wiener-SVD w/incorrect prediction

2

C τ

7 −

10

6 −

10

5 −

10

4 −

10

3 −

10

2 −

10 MSE

4

10

5

10

6

10

2

C

τ 0.010.020.030.040.05

3 −

10 × MSE 4300 4400 4500 4600

Using 𝑡̅ to represent 𝑡.&/$ Using an incorrect 𝑡 to represent 𝑡.&/$

slide-41
SLIDE 41

Discussions and recommendations (1)

§ Pros & cons of Wiener-SVD

§ It does not require the evaluation of regularization strength. § It is model-dependent. However, the dependence can be mitigated by taking into account a number of different models, and by reporting the additional smearing matrix.

§ Both the traditional regularization or Wiener-SVD are equivalent to applying a new smearing matrix, which suppresses the large oscillation (high variance) of the direct matrix inversion unfolded results, but also introduces bias.

Xiaoyue Li DPF2017 - FNAL 41

𝑡̂ = 𝐵• I 𝑆h𝑆 GH I 𝑆h I 𝑁 Ž Where 𝐵• = 𝐷GH I 𝑊

  • I 𝑋
  • I 𝑊
  • h I 𝐷

(𝑋

  • )35=

𝑒3

b I ∑ 𝑊

  • h

35 I ∑ 𝐷5– I 𝑡̅– – 5 b

𝑒3

b I ∑

𝑊

  • h

35 I ∑ 𝐷5– I 𝑡̅– – 5 b + 1

slide-42
SLIDE 42

Discussions and recommendations (2)

§ The Wiener filter should be constructed in ways which takes into account a range

  • f prior expectations.

§ C2 matrix generally yields better results than C0 and C1. § We recommend reporting this additional smearing matrix 𝐵• in the publication together with the unfolded results to enable a more direct comparison of expectations (e.g. from new theoretical calculations) with the unfolded results. In practice, the new smearing matrix should be applied to the theoretical calculation before comparing to the unfolded results.

Xiaoyue Li DPF2017 - FNAL 42

https://github.com/BNLIF/Wiener-SVD-Unfolding arXiv:1705.03568

slide-43
SLIDE 43

Supplemental

Xiaoyue Li DPF2017 - FNAL 43

slide-44
SLIDE 44

Regularization interpretation of Wiener-SVD method

§ The Wiener-SVD unfolding method is equivalent to regularization w.r.t. minimizing the signal to noise ratio in the effective frequency domain 𝑁v

3 =

∑ 𝑊

35 h𝑡5 5

: 𝜚 𝑡 = log𝑀(𝑡) + 1 2— log 𝑁v

3 b

𝑂 b

3

= log𝑀 𝑡 + 1 2 —log ∑ 𝑊

35 h𝑡5 5 b

1

3

Minimizing Eq. (7) yields the following estimator 𝑡̂ = −𝑌GH I 𝑍 I 𝑁v Where 𝑌35 = 𝜖b𝜚b 𝜖𝑡3𝜖𝑡5 = − 𝑆h𝑆 35 − — 𝑊

3~ I

1 𝑁v ~

b I 𝑊 ~5 h ~

𝑍

35 =

𝜖b𝜚b 𝜖𝑡3𝜖𝑁v5 = 𝑆35

h

§ With 𝑌 and 𝑍 evaluated at the expectation of 𝑡̅ and 𝑁 Ž, one recovers Eq.(6) 𝑡̂ = 𝑊 I 𝑋 I 𝑊h I 𝑆h I 𝑁 Ž With 𝑋

35 =

𝑒3

b

𝑒3

b+ 1

𝑁v

b 3

𝜀35 = 𝑒3

b I ∑ 𝑊 35 h I 𝑡̅5 5 b

𝑒3

b I ∑ 𝑊 35 h I 𝑡̅5 5 b + 1

𝜀35

Xiaoyue Li DPF2017 - FNAL 44