Statistics and Numerics

Statistics and Numerics

Lecture Prof. Dr. Jens Timmer

Exercises Helge Hass, Mirjam Fehling-Kaschek

Exercise Sheet Nr. 6

Exercise 1: Bias and variance of ill-posed inverse problems

Create a function to compute theN×NHilbert-matrix A_{i j}= 1

i+j−1 given the dimensionNas input.

a) Generate simulated data~bfor normal distributed noise:

• choose:

x_j=sin(2π(j−1)/(N−1)), j=1, . . . ,N

• compute:

b˜_i=

N

∑

A_{i j}x_j (1)

• add noiseε_i∼N(0,σ²):

b_i=b˜_i+ε_i b) Estimate~xfrom the data~b:

• Calculate the singular-value decomposition[U,s,V]of A(numpy.linalg.svd). How do the singular valuess_irelate to the eigenvalues in case of a symmetric matrix?

• Compute the condition numberκ(A)vianumpy.aminandnumpy.amax.

• Write a function to estimatex_ifor the givenb_ifrom Eq. (1), using the inverted matrixA⁻¹based on its singular value decomposition:

A⁻¹=V^T[1/s_i]U^T . (2)

Use the built-in functionsnumpy.diag,numpy.divide,numpy.matmulandnumpy.transpose.

c) Test several setups regarding the effect of the regularization, 1/s_i=0i f max(s)/s_i>κ_max, for different choices of the maximal condition numberκ_max, ofNand of the noiseσ.

• N=4,σ=0.001 for either no regularization or forκmax=1000

• N=7,σ=10⁻⁵,κmax∈ {100,10⁵,10¹⁰}

• N=10,σ=0,κ_max=in f (float(“inf”) in python). Compare the result with the one you obtain from using an in-build matrix inversion function (numpy.linalg.inv) instead of Eq. (2). How areA⁻¹and thex_iaffected?

• N=42,σ=0,κ_max∈ {100,10⁵,10¹⁰}

helge.hass@fdm.uni-freiburg.de mirjam.fehling@physik.uni-freiburg.de

http://jeti.uni-freiburg.de/vorles_stat_num/vorles_stat_num.html