Inverse problems, spring 2014

== Inversio-ongelmat, kevät 2014 ==

== Inverse problems, spring 2014 ==

=== (% style="color: rgb(0, 0, 0); font-size: 10pt; line-height: 13pt; color: rgb(255, 0, 0)" %)**The course is lectured in English.**(%%) ===

6

7

[[image:attach:Nut3web.jpg]]

8

9

== (% style="font-size: 10.0pt;line-height: 13.0pt;" %)Introduction(%%) ==

10

11

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)Inverse problems research is an active area of mathematics.

12

13

At the Department of Mathematics and Statistics the field is represented by three research groups belonging to(% style="font-size: 10.0pt;line-height: 13.0pt;" %)

14

a [[Centre of Excellence of Academy of Finland>>url:http://wiki.helsinki.fi/display/inverse/Home;jsessionid=BCF26178580CA3F00E25F4F6A441879A||shape="rect"]]:

15

16

* (% style="font-size: 10.0pt;line-height: 13.0pt;" %)[[The team of Samuli Siltanen>>url:http://wiki.helsinki.fi/display/inverseproblems/Computational+Inverse+Problems||shape="rect"]]

17

* (% style="font-size: 10.0pt;line-height: 13.0pt;" %)[[The team of Matti Lassas>>url:http://wiki.helsinki.fi/display/inverseproblems/Geometric+Inverse+problems+and+applications||shape="rect"]]

18

* (% style="font-size: 10.0pt;line-height: 13.0pt;" %)[[The team of Aapo Hyvärinen>>url:http://www.hiit.fi/neuro||shape="rect"]]

19

20

Inverse problems are about interpreting indirect measurements, and they have

21

applications to many areas of science and technology. Examples include

22

23

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)-Sharpening a misfocused photograph

24

-Three-dimensional X-ray imaging ([[more information>>url:http://www.siltanen-research.net/project_Xray.html||shape="rect"]])

25

-Recovering the inner structure of the Earth based on earthquake measurements

26

-Reconstructing electric conductivity from current-to-voltage boundary measurements (see [[this page>>url:https://wiki.helsinki.fi/display/inverse/Electrical+Impedance+Tomography||shape="rect"]] and [[this page>>url:http://www.siltanen-research.net/project_EIT.html||shape="rect"]])

27

-Finding cracks inside solid structures

28

-Prospecting for oil and minerals

29

-Monitoring underground contaminants

30

-Finding the shape of asteroids based on light-curve data (see [[this page>>url:http://www.rni.helsinki.fi/~~mjk/asteroids.html||shape="rect"]])

31

32

33

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)The common features of all this problems are the need to understand indirect measurements and to overcome extreme sensitivity to noise and modelling inaccuracies.

34

35

=== What does the course contain? ===

36

37

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)The goals of the course are

38

39

* (% style="font-size: 10.0pt;line-height: 13.0pt;" %)introduce discrete matrix models of some widely used measurements, such as tomography and convolution

40

* (% style="font-size: 10.0pt;line-height: 13.0pt;" %)show how to detect ill-posedness (sensitivity to measurement noise) in matrix models using Singular Value Decomposition

41

* (% style="font-size: 10.0pt;line-height: 13.0pt;" %)compute noise-robust reconstructions using //regularization//

42

* (% style="font-size: 10.0pt;line-height: 13.0pt;" %)write Matlab algorithms for sharpening photographs and computing tomographic reconstructions

43

44

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)The lectures make up 10 credit units. In addition to lectures the course involves a **project work**. It is done in teams of two and gives 5 credit units to each student.

45

46

The course is in total **15 credit units**.

----

(% style="color: rgb(255,0,0);" %)Left: original image. Middle: misfocused image. Right: sharpened image.

51

52

[[image:attach:hammer.jpg]]

----

=== Lecturer ===

[[Samuli Siltanen>>doc:mathstatHenkilokunta.Siltanen, Samuli]]

=== Assistants ===

[[Esa Niemi>>url:http://wiki.helsinki.fi/display/mathstatHenkilokunta/Niemi%2C+Esa||shape="rect"]], [[Teemu Saksala>>url:http://wiki.helsinki.fi/display/mathstatHenkilokunta/Teemu+Saksal||shape="rect"]]

=== Credit units ===

15 op.

=== Type ===

Advanced studies (syventävä opintojakso)

73

74

=== Prerequisites ===

75

76

Recommended courses to take before this course: Linear algebra 1 and 2, Applications of matrix computations.

77

78

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)Some previous experience with Matlab programming is very helpful. Tutorials for learning basics of Matlab are available for example at the following locations:

79

80

* [[MathWorks tutorials>>url:http://www.mathworks.com/academia/student_center/tutorials/launchpad.html||shape="rect"]]

81

* [[University of Utah tutorials>>url:http://www.math.utah.edu/lab/ms/matlab/matlab.html||shape="rect"]]

82

* [[Clarkson University tutorials>>url:http://www.cyclismo.org/tutorial/matlab/||shape="rect"]]

----

== Project work ==

There are two project work topics: image deblurring and real-data X-ray tomography.

89

90

The idea is to study an inverse problem both theoretically and computationally in **teams of two students**. The end product is a scientific poster that the team will present in a poster session on **Thursday, May 8, 2014 at 14-16**in the Industrial Mathematics Laboratory (Exactum C131). The poster can be printed using the laboratory's large scale printer. (% style="font-size: 10.0pt;line-height: 13.0pt;" %)The classical table of contents is recommended for structuring the poster:

91

92

1 Introduction

93

2 Materials and methods

3 Results

4 Discussion

Section 2 is for describing the data and the inversion methods used. In section 3 those methods are applied to the data and the results are reported with no interpretation; just facts and outcomes of computations are described. Section 4 is the place for discussing the results and drawing conclusions.

98

99

Here are more specific instructions concerning the two project work topics.

100

101

//Image deblurring~:// think about an experiment where you can take three photos of the same object: sharp, a little blurred (misfocused) and a bit more blurred. The object should be essentially two-dimensional, such as text or image on paper or a coin. You can either take three separate photos using three different focus settings in the camera, or you can place three identical objects at three different distances from the lens. It is advisable that each object contains a small black dot on a white surface; this will give you an empirical point spread function. The aim of the project work is to implement and test a deblurring algorithm for the photographs and compare the final result to sharply photographed ground truth. You can choose the large-scale inversion method you want to use in the project. Preferably, take one of the following:

102

103

* (% style="font-size: 10.0pt;line-height: 13.0pt;" %)matrix-free Tikhonov with the conjugate gradient method,

104

* (% style="font-size: 10.0pt;line-height: 13.0pt;" %)matrix-free approximate total variation with the Barzilai-Borwein method.

105

106

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)If you plan to use some other method, discuss this plan with the lecturer in the First Goal meeting. Optimally, you should have an automatic method for choosing the regularization parameter and an automatic stopping criteria for the iteration. These are both difficult requirements, so have a simple approach as plan B if a more complicated approach does not work. Here is a [[test sheet>>attach:Testsheet.pdf]](%%) (% style="font-size: 10.0pt;line-height: 13.0pt;" %)used in previous years' courses, you can use it if you want. Also, here are extra instructions which you need not follow in detail: [[attach:project_deblur.pdf]](%%). **Photographic d(% style="font-size: 10.0pt;line-height: 13.0pt;" %)ata collection will be done on Friday, 4.4.2014, at 12:15-14:00(%%)**(% style="font-size: 10.0pt;line-height: 13.0pt;" %) in the Industrial Mathematics Laboratory (Exactum C131).

107

108

//X-ray tomography~:// Download the set of measured walnut [[data>>attach:CTdataIP2014.zip]]. The dataset contains

109

110

* X-ray projection data in sinogram form,

111

* measurement matrix, and

112

* a filtered back-projection reconstruction from complete data for comparison.

113

114

In the project you are supposed to take a subset of the data with only few projections, such as 20 projection directions. Use one of these methods to recover the walnut slice from sparse data:

115

116

* Tikhonov regularization based on conjugate gradient method,

117

* approximate total variation regularization implemented iteratively with the Barzilai-Borwein method.

118

119

Optimally, you should have an automatic method for choosing the regularization parameter and an automatic stopping criteria for the iteration. These are both difficult requirements, so have a simple approach as plan B if a more complicated approach does not work.\\

120

121

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)**First goal consists of two things:** (a) two first sections should be preliminary written in LaTeX (not necessarily in poster format yet) and (b) the Matlab codes at the following webpages should be run and studied (only the page related to your topic):

122

123

* [[Image deblurring>>url:http://wiki.helsinki.fi/display/mathstatHenkilokunta/2D+deconvolution||shape="rect"]]

124

* [[X-ray tomography with sparse data>>url:http://wiki.helsinki.fi/display/mathstatHenkilokunta/Matrix-free+X-ray+tomography+with+sparse+data||shape="rect"]]

125

126

Two things will be graded in the meeting about the first goal: (a) the draft of project work and (b) your understanding of the Matlab codes at the above webpage relevant to your topic. The grade represents 30% of the final grade of the project work. Please agree on a meeting time (in the period April 1-4) with the lecturer for reviewing and grading the first goal.

127

128

**Second and final goal:** poster is presented in the session on May 8. The poster will be printed in size A1. You may create your own poster (from scratch), or you can use e.g. [[this template>>attach:posterA1_templ_IP2014.zip]] as a starting point and edit its layout, colors, fonts, etc. as much as you like.

129

Please send your poster via email as a pdf attachment to Esa Niemi by Monday 5th May, 12 pm. Then your poster will be printed by the poster session on May 8.

130

**The poster session will be held on Thursday, May 8, 14:15-15:30** in the Industrial Mathematics Laboratory (Exactum C131).

----

== Lectures ==

(% style="color: rgb(128,0,128);" %)**Period III:**(%%) Lectures as follows:

139

140

Tuesday 10-12 in room D123

141

Wednesday 12-14 in room D123

142

Friday 12-14 in room C123.

143

144

Two hours of exercise classes per week.

145

146

(% style="color: rgb(128,0,128);" %)**Period IV:**(%%) Lectures and exercises in the beginning of the period. Later project work, which is reported as a poster in a poster session.

147

148

**14.1.2014 Tuesday:** Introduction to inverse problems. Specific problems discussed were [[image deblurring>>attach:Deblurring.pdf]], [[X-ray tomography>>url:http://www.siltanen-research.net/project_Xray.html||shape="rect"]], [[electrical impedance tomography>>url:http://www.siltanen-research.net/project_EIT.html||shape="rect"]] and [[glottal inverse filtering>>attach:GIF.pdf]].

149

Lecture slides: [[attach:Inverse_introduction2_web.pdf]], [[attach:GIF.pdf]], [[attach:Deblurring.pdf]].

150

Tomography videos: [[construction of the sinogram>>url:http://www.youtube.com/edit?video_id=8_QU2k4TnSk&video_referrer=watch||shape="rect"]], [[reconstruction using Filtered Back-Projection>>url:http://www.youtube.com/edit?video_id=-SiVPCh92TA&video_referrer=watch||shape="rect"]].

151

Book sections (Mueller-Siltanen 2012): Chapter 1.

152

153

**15.1.2014 Wednesday:** Basics of one-dimensional convolution, including the continuous model and its matrix approximation.

154

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)Book sections (Mueller-Siltanen 2012): Sections 2.1.1. and 2.1.2.

155

156

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)Matlab codes: [[attach:DC_convmtx.m]], [[attach:DC_PSF.m]], [[attach:DC_PSF_plot.m]], [[attach:DC_target.m]], [[attach:DC1_cont_data_comp.m]], [[attach:DC1_cont_data_plot.m]], [[attach:DC2_discretedata_comp.m]], [[attach:DC2_discretedata_plot.m]].

157

158

(% style="color: rgb(0,0,0);" %)**21.1.2014 Tuesday:**(%%) Analysis of 1-dimensional discrete deconvolution continues. Assessment of ill-posedness of convolution matrices using Singular Value Decomposition (SVD) and observing the decay of singular values to zero. Numerical experiments showing the failure of naive inversion. A first look at robust reconstructions based on truncated SVD.

159

160

Book sections: 2.1.3, 2.1.4, 3.5, 3.6.

161

162

Matlab codes: [[attach:DC3_naive_plot.m]], [[attach:DC4_show_SVD.m]], [[attach:DC4_truncSVD_comp.m]].

163

164

(% style="color: rgb(0,0,0);" %)**22.1.2014 Wednesday:**(%%) Overview of the topics that will be covered in the course:

165

-Motivation using 1D deconvolution and 2D tomography

166

-Theory of inverse problems: existence, uniqueness, stability, and regularization

167

-Regularization methods: truncated SVD, Tikhonov regularization, total variation regularization

168

-Large-scale computation for practical problems using iterative techniques

169

170

After the overview we moved on to study the basics of X-ray attenuation and two-dimensional tomography. We studied the ill-posedness of tomography by looking at the singular values of tomography matrices. Also, naive inversion (based on least-squares solution) fails in the case of 2D tomography as it failed for 1D deconvolution. Finally, we had a second look at robust reconstructions based on truncated SVD.

171

172

Lecture slides: [[attach:Xray_tomo.pdf]]

173

174

Book sections: 2.3.1, 2.3.2, 2.3.4.

175

176

Matlab codes: [[attach:phantom.m]], [[attach:XRM1_matrix_comp.m]],

177

[[attach:XRM2_naive_comp.m]], [[attach:XRM2_naive_plot.m]],

178

[[attach:XRM3_NoCrimeData_comp.m]], [[attach:XRM3_NoCrimeData_plot.m]], [[attach:XRM4_naive_comp.m]], [[attach:XRM4_naive_plot.m]],

179

[[attach:XRM5_SVD_comp.m]], [[attach:XRM5_SVD_plot.m]], [[attach:XRM6_truncSVD_comp.m]]

180

Remarks about the codes:

181

182

(% style="list-style-type: square;" %)

183

* The routines need to be run in the above order for each value of //n// separately.

184

* You will need the image processing toolbox; otherwise the routine //radon.m// is not available for you.

185

* (% style="font-size: 10.0pt;line-height: 13.0pt;" %)You need to have a subfolder called "data" in your working folder. That's where the data files of type ".mat" will be saved. The routine //XRM1_matrix_comp.m// actually creates the folder "data" automatically for you; if you run it again, Matlab may give you a warning that the directory already exists. You can ignore this warning.

186

* (% style="font-size: 10.0pt;line-height: 13.0pt;" %)If you run the routine XRM3_NoCrimeData_comp.m with one value of //n// first and then with another value of //n//, there is an error message. I have no idea why. If you find the reason please let me know. Anyway, type "clear all" in Matlab after seeing the error message and you're good to go.

187

* (% style="font-size: 10.0pt;line-height: 13.0pt;" %)Computing XRM5_SVD_comp.m can take a really long time when N is large (32 or bigger). For your convenience, I upload here a couple of SVD results: [[attach:XRME_SVD16.mat]], [[attach:XRME_SVD32.mat]].

188

189

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)**24.1.2014 Friday: **Introduction to the general framework of inverse problems. Least squares solution and minimum norm solution for matrix equations.

190

191

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)Lecture slides: [[attach:GeneralTheory.pdf]]

192

193

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)Book sections: 3.1, 3.4 and 4.1.

194

195

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)**28.1.2014 Tuesday:** Continuous model of one-dimensional convolution and deconvolution, including Fourier transform analysis.

196

197

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)Radon transform in X-ray tomography. Fourier slice theorem.

198

199

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)Lecture notes: [[Convolution and Fourier transform>>attach:FourierConv.pdf]], [[attach:Somersalo_Inversio_Liite2.pdf]] (the latter is in Finnish).

200

201

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)Book section: 2.3.3.

202

203

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)**29.1.2014 Wednesday:** Two different derivations of the Filtered Back-Projection (FBP) algorithm for X-ray tomography.

204

205

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)Complementary lecture notes are available [[in this file>>attach:FBP.pdf]].

206

207

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)This is the Matlab file created in the end of the lecture: [[attach:Radon_example.m]]

208

209

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)Book section: 2.3.3.

210

211

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)**31.1.2014 Friday:** Invisible structures in X-ray tomography. Interplay between continuous and discrete tomography. Introduction to Tikhonov regularization.

212

213

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)Lecture notes: [[attach:Xray_tomo_web1.pdf]] [[attach:SmithSolmonWagner1977.pdf]]

214

215

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)Book section 5.1.

216

217

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)**4.2.2014 Tuesday:** Tikhonov regularization continued. Derivation of normal equations. Writing the Tikhonov regularization task in stacked form.

218

219

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)Matlab file: [[attach:DC5_Tikhonov_comp.m]]

220

221

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)Book section 5.2.

222

223

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)**5.2.2014 Tuesday:** Generalized Tikhonov regularization and the corresponding normal equations. A few words and numerical examples about iterative solution of linear equations. The L-curve method for choosing the regularization parameter in Tikhonov regularization.

224

225

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)Matlab routines: [[attach:DC5_Tikhonov_Lcurve.m]], [[attach:DC6_TikhonovD_comp.m]], [[attach:iterfun.m]], [[attach:itersoltest.m]]. For matrix-free iterative regularization of tomography, see [[this page>>url:http://wiki.helsinki.fi/display/mathstatHenkilokunta/Matrix-free+X-ray+tomography||shape="rect"]] and [[this page>>url:http://wiki.helsinki.fi/display/mathstatHenkilokunta/Matrix-free+X-ray+tomography+with+sparse+data||shape="rect"]].

226

227

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)Book sections: 5.3, 5.4.2.

228

229

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)**7.2.2014 Friday:** Conjugate gradient method for large-scale Tikhonov regularization.

230

231

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)Book section: 5.5.

232

233

(% style="color: rgb(0,0,0);" %)**18.2.2014 Tuesday:** Total variation regularization.

234

235

(% style="color: rgb(0,0,0);" %)Book sections: 6.1 and 6.2.

236

237

(% style="color: rgb(0,0,0);" %)Matlab resources: [[attach:DC7_TotalVariation_comp.m]] (note that Matlab Optimization Toolbox is needed)

238

239

(% style="color: rgb(0,0,0);" %)**19.2.2014 Wednesday:** The S-curve method for sparsity-based determination of TV regularization parameter. Large-scale implementation of TV regularization via smoothing and Barzilai-Borwein minimization.

240

241

(% style="color: rgb(0,0,0);" %)Book sections: 6.3 and 6.4.

242

243

(% style="color: rgb(0,0,0);" %)Matlab resources (convolution, Matlab Optimization Toolbox needed):

244

[[attach:DC8_TV_sparsitychoice_comp.m]](% style="font-size: 10.0pt;line-height: 13.0pt;" %), [[attach:DC9_TV_sparsitychoice_plot.m]]

245

246

(% style="color: rgb(0,0,0);" %)Matlab resources (matrix-based tomography):

247

[[attach:XRMH_aTV_comp.m]], [[attach:XRMH_aTV_plot.m]], [[attach:XRMH_misfit.m]], [[attach:XRMH_misfit_grad.m]],

248

[[attach:XRMH_aTV.m]], [[attach:XRMH_aTV_grad.m]], [[XRMH_aTV_feval,m>>attach:XRMH_aTV_feval.m]], [[attach:XRMH_aTV_fgrad.m]]

249

250

(% style="color: rgb(0,0,0);" %)Matlab resources (matrix-free large-scale tomography):

251

[[see this page>>url:http://wiki.helsinki.fi/display/mathstatHenkilokunta/Matrix-free+X-ray+tomography+with+sparse+data||shape="rect"]].

252

253

(% style="color: rgb(255,0,0);" %)**14.3.2014 Friday:** Home exam return session & project work assignment.

----

== Exercises ==

Tuesday 16-18 in room D210 (Physicum)

260

Friday 14-16 in room C128 (Exactum)

261

262

You may choose which one of the exercise sessions you participate. The theoretical exercises (T1, T2, ...) should be done in advance and they will be checked in the beginning of the exercise session. The Matlab exercises (M1, M2, ...) as well as the LaTeX exercises (L1, L2, ...) can be done in advance or during the exercise session.

263

\\

264

265

[[Exercise 1>>attach:Laskari01.pdf]] (January 21-24, 2014), [[Solutions 1>>attach:Ex1 answers.pdf]]

266

267

[[Exercise 2>>attach:Laskari02.pdf]] (January 28-31, 2014), [[Solutions 2>>attach:Ex2 answers.pdf]]

268

269

[[Exercise 3>>attach:Laskari03.pdf]] (February 4-7, 2014), [[attach:DC_PSF.m]], [[attach:Somersalo_Inversio_Liite2.pdf]],[[ Solutions 3>>attach:Laskari03.pdf]]

270

271

[[Exercise 4>>attach:Laskari04.pdf]] (February 11-14, 2014), [[attach:reportdraft1.tex]], [[Solutions 4>>attach:Laskari04.pdf]]

272

273

[[Exercise 5>>attach:Laskari05.pdf]] (February 18-21, 2014), feel free to modify this file: [[attach:DC08_Tikhonov_Morozov1.m]] (the routine will give an error message without appropriate modification). [[Solutions 5>>attach:Laskari05.pdf]]

(((

----

)))

== Course Materials ==

280

281

The course follows the book //Linear and Nonlinear Inverse Problems and Practical Applications //by J.L. Mueller and S. Siltanen (SIAM 2012).

282

Part I of the book will be covered.

283

284

The book is available at [[Google books>>url:http://books.google.fi/books?id=nDa83ldfqf8C&lpg=PP1&pg=PP1#v=onepage&q&f=false||shape="rect"]] and at the [[SIAM bookstore>>url:http://www.ec-securehost.com/SIAM/CS10.html||shape="rect"]].

285

286

Also, the relevant chapters of the book are available for copying in Exactum room C326. **You can now also find chapters 6 and 9 from the copy room**.

287

288

The Matlab codes used in the book are collected to [[this page>>url:http://wiki.helsinki.fi/display/mathstatHenkilokunta/Inverse+Problems+Book+Page||shape="rect"]].

289

290

[[image:attach:FrontCover_small.png]]\\

----

=== Exam ===

The answers to the exam questions should be returned in the (% style="font-size: 10.0pt;line-height: 13.0pt;" %)examination session on Friday, March 14, 2014, at 12:15. The hall is Exactum C123.

297

298

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)More instructions are given in the exam itself: [[attach:2014_home_exam.pdf]].

299

300

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)Here is the draft report file for Problem 3: [[attach:Laplacereport.tex]].

301

302

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)Matlab resources:

303

304

(% style="font-size: 10.0pt;line-height: 13.0pt;" %)[[attach:XR_aTV_feval.m]]