numerical linear algebra using interval arithmetic

108 Pages·2012·1.16 MB·English

THÈSE
en vue de l’obtention du grade de
Docteur de l’École Normale Supérieure de Lyon - Université de Lyon
Spécialité : Informatique
Laboratoire de l’Informatique du Parallélisme
École Doctorale de Mathématiques et Informatique Fondamentale
présentée et soutenue publiquement le 18 janvier 2011
par M. Hong Diep NGUYEN
Efficient algorithms for verified scientific computing:
numerical linear algebra using interval arithmetic
Directeur de thèse : M. Gilles VILLARD
Co-directrice de thèse : Mme Nathalie REVOL
Après l’avis de : M. Philippe LANGLOIS
M. Jean-Pierre MERLET
Devant la commission d’examen formée de :
M. Luc JAULIN, Président
M. Philippe LANGLOIS, Rapporteur
M. Jean-Pierre MERLET, Rapporteur
Mme Nathalie REVOL, Co-directrice
M. Siegfried RUMP, Membre
M. Gilles VILLARD, Directeur

Efficient algorithms for verified scientific computing:
numerical linear algebra using interval arithmetic
˜ `
Nguyên Hông Diê.p
18 January 2011

Contents
1 Introduction 1
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Description of the document . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Interval arithmetic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3.1 Operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3.2 Algebraic properties . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3.3 Absolute value properties . . . . . . . . . . . . . . . . . . . . . . . 8
1.3.4 Midpoint and radius properties . . . . . . . . . . . . . . . . . . . . 8
1.3.5 Variable dependency . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.3.6 Wrapping effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.3.7 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3.8 Implementation environment . . . . . . . . . . . . . . . . . . . . . . 15
1.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2 Interval matrix multiplication 17
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2 State of the art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.3 Algorithm 1 using endpoints representation . . . . . . . . . . . . . . . . . . 21
2.3.1 Special cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.3.2 General case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.3.3 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.4 Algorithm 2 using midpoint-radius representation . . . . . . . . . . . . . . 34
2.4.1 Scalar interval product . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.4.2 Matrix and vector operations . . . . . . . . . . . . . . . . . . . . . 38
2.4.3 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3 Verification of the solution of a linear system 45
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.2 State of the art, verifylss function . . . . . . . . . . . . . . . . . . . . . 46
3.3 Floating-point iterative refinement . . . . . . . . . . . . . . . . . . . . . . 48
3.4 Interval iterative refinement . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.4.1 Numerical behavior . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.4.2 Updating the computed solution . . . . . . . . . . . . . . . . . . . . 52
3.4.3 Initial solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.4.4 Interval iterative improvement . . . . . . . . . . . . . . . . . . . . . 54
3.4.5 Detailed algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
i
ii CONTENTS
3.5 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.5.1 Implementation issues . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.5.2 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . 59
3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4 Relaxation method 63
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.2 Interval improvement: Jacobi and Gauss-Seidel iterations . . . . . . . . . . 64
4.2.1 Convergence rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.2.2 Accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.2.3 Initial solution of the relaxed system . . . . . . . . . . . . . . . . . 70
4.3 Implementation and numerical experiments . . . . . . . . . . . . . . . . . . 71
4.3.1 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.3.2 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5 Effects of the computing precision 75
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.2 Implementation using variable computing precisions . . . . . . . . . . . . . 75
5.2.1 Residual computation . . . . . . . . . . . . . . . . . . . . . . . . . 76
5.2.2 Floating-point working precision . . . . . . . . . . . . . . . . . . . . 76
5.2.3 Needed working precision . . . . . . . . . . . . . . . . . . . . . . . 77
5.3 Error analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
5.3.1 Normwise relative error . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.3.2 Componentwise relative error . . . . . . . . . . . . . . . . . . . . . 83
5.4 Doubling the working precision for approximate solution . . . . . . . . . . 84
5.4.1 Implementation issues . . . . . . . . . . . . . . . . . . . . . . . . . 84
5.4.2 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
5.4.3 Other extended precisions . . . . . . . . . . . . . . . . . . . . . . . 85
5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
Conclusions and Perspectives 89
List of Figures
1.1 Wrapping effect. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.1 Performance of igemm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.2 Performance of igemm3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.1 Solution set of the interval linear system (3.3). . . . . . . . . . . . . . . . . 51
3.2 MatLab implementation results for 1000 × 1000 matrices. . . . . . . . . . 60
4.1 Solution set of the relaxed system. . . . . . . . . . . . . . . . . . . . . . . . 65
4.2 Relaxation technique: MatLab implementation results for 1000×1000 ma-
trices. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5.1 Effects of the residual precision. . . . . . . . . . . . . . . . . . . . . . . . . 77
5.2 Effects of the working precision. . . . . . . . . . . . . . . . . . . . . . . . . 78
5.3 Minimal working precision needed. . . . . . . . . . . . . . . . . . . . . . . 78
5.4 Effect of doubling the working precision for the approximate solution. . . . 87
5.5 Effect of other extended precisions. . . . . . . . . . . . . . . . . . . . . . . 88
iii
iv LIST OF FIGURES
List of Algorithms
2.1 Decomposition algorithm for endpoint intervals . . . . . . . . . . . . . . . 29
3.1 Classical iterative refinement . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.2 Interval iterative refinement . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.3 Testing the H-matrix property . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.4 Solving and verifying a linear system . . . . . . . . . . . . . . . . . . . . . 56
4.1 Interval improvement function refine using relaxed Jacobi iterations . . . 72
5.1 Solving and verifying a linear system . . . . . . . . . . . . . . . . . . . . . 76
5.2 Doubling the precision for the approximate solution . . . . . . . . . . . . . 86
v
vi LIST OF ALGORITHMS