Research
Matters
How
and Why
to Estimate
Condition February
Numbers
Matrix Functions
25,for
2009
Nick
Higham
Nick
Higham
Director
SchoolofofResearch
Mathematics
The
University
of Manchester
School
of Mathematics
http://www.maths.manchester.ac.uk/~higham
@nhigham, nickhigham.wordpress.com
Joint work with Sam Relton.
Computational Linear Algebra and Optimization
for the Digital Economy
1/6
Edinburgh, October 2013
Sources of Error
f : Cn×n → Cn×n .
Want to compute f (A) but given A + ∆A not A.
∆A may come from
Rounding errors in storing A.
Measurement errors.
Errors from an earlier computation.
Nick Higham
Matrix Function Condition Numbers
2 / 33
Sources of Error
f : Cn×n → Cn×n .
Want to compute f (A) but given A + ∆A not A.
∆A may come from
Rounding errors in storing A.
Measurement errors.
Errors from an earlier computation.
Photocopying errors!
Nick Higham
Matrix Function Condition Numbers
2 / 33
Photocopier
http://en.wikipedia.org/wiki/Photocopier
“Most current photocopiers use a technology called
xerography, a dry process that uses electrostatic
charges on a light sensitive photoreceptor to first attract
and then transfer toner particles (a powder) onto paper
in the form of an image”
“There is an increasing trend for new photocopiers to
adopt digital technology, thus replacing the older
analog technology. With digital copying, the copier
effectively consists of an integrated scanner and laser
printer.”
Nick Higham
Matrix Function Condition Numbers
3 / 33
Xerox WorkCentre 7535, 7556
August 2013: German computer scientist David Kriesel
discovered the machines often change “6” to “8”.
Jbig2 compression algorithm implicated.
Nick Higham
Matrix Function Condition Numbers
4 / 33
Maxim
Every code for solving a numerical problem should return
an error estimate or bound with the computed result.
Nick Higham
Matrix Function Condition Numbers
5 / 33
Conditioning of f (x)
Scalar function f ∈ C 2 , and y = f (x), y + ∆y = f (x + ∆x).
Then
0 ∆y
xf (x) ∆x
=
+ O ∆x 2 .
y
f (x)
x
The relative condition number
0 xf (x) .
c(x) = f (x) Nick Higham
Matrix Function Condition Numbers
6 / 33
Conditioning of f (x)
Scalar function f ∈ C 2 , and y = f (x), y + ∆y = f (x + ∆x).
Then
0 ∆y
xf (x) ∆x
=
+ O ∆x 2 .
y
f (x)
x
The relative condition number
0 xf (x) .
c(x) = f (x) Condition number of the condition number: assuming
xf 0 (x) > 0 and f (x) > 0,
0 00
xc (x) f (x) f 0 (x) [2]
c (x) := = 1+x
−
.
c(x) f 0 (x)
f (x) Nick Higham
Matrix Function Condition Numbers
6 / 33
Fréchet Derivative
Fréchet derivative of f : Cn×n → Cn×n at X ∈ Cn×n
A linear mapping Lf : Cn×n → Cn×n s.t. for all E ∈ Cn×n
f (X + E) − f (X ) − Lf (X , E) = o(kEk).
Nick Higham
Matrix Function Condition Numbers
7 / 33
Fréchet Derivative
Fréchet derivative of f : Cn×n → Cn×n at X ∈ Cn×n
A linear mapping Lf : Cn×n → Cn×n s.t. for all E ∈ Cn×n
f (X + E) − f (X ) − Lf (X , E) = o(kEk).
Examples:
(X + E)2 − X 2 = XE + EX + E 2 ⇒ Lx 2 (X , E) = XE + EX .
(X + E)−1 − X −1 = −X −1 EX −1 + O(E 2 )
⇒ Lx −1 (X , E) = −X −1 EX −1 .
Nick Higham
Matrix Function Condition Numbers
7 / 33
Fréchet Derivative
Fréchet derivative of f : Cn×n → Cn×n at X ∈ Cn×n
A linear mapping Lf : Cn×n → Cn×n s.t. for all E ∈ Cn×n
f (X + E) − f (X ) − Lf (X , E) = o(kEk).
Examples:
(X + E)2 − X 2 = XE + EX + E 2 ⇒ Lx 2 (X , E) = XE + EX .
(X + E)−1 − X −1 = −X −1 EX −1 + O(E 2 )
⇒ Lx −1 (X , E) = −X −1 EX −1 .
Scalar case: Lf (x, e) = f (0 x)e.
Nick Higham
Matrix Function Condition Numbers
7 / 33
Gâteaux Derivative
Gâteaux derivative of matrix function f at A in direction E:
f (A + E) − f (A)
.
→0
Gf (A, E) = lim
Weaker notion of differentiability.
Nick Higham
Matrix Function Condition Numbers
8 / 33
Applications of Fréchet Derivative
Computation of correlated choice probabilities (2013),
registration of MRI images (2013),
Markov models applied to cancer data (2013),
matrix geometric mean computation (2012),
model reduction (2012),
first and second Fréchet derivative in Halley’s method
on Banach space (2003).
Nick Higham
Matrix Function Condition Numbers
9 / 33
Componentwise Sensitivity (1)
Let E = ei ejT . Then
f (X + ei ejT ) − f (X )
lim
= Lf (X , ei ejT ).
→0
(Lf (X , ei ejT ))rs measures the sensitivity of f (X )rs to
perturbations in aij .
Nick Higham
Matrix Function Condition Numbers
10 / 33
Componentwise Sensitivity (2)
Since Lf is a linear operator,
vec(Lf (A, E)) = Kf (A)vec(E)
where Kf (A) ∈ Cn
derivative.
2 ×n2
is the Kronecker form of the Fréchet
a11 a21 a31 . . . ann
f11
f21
f31
..
.
Kf (A) ≡
Lf (A, ei ejT )
fnn
Nick Higham
Matrix Function Condition Numbers
11 / 33
Condition Number
kf (A + E) − f (A)k
.
→0 kEk≤
condabs (f , A) = lim sup
kLf (A)k := max
E6=0
kLf (A, E)k
.
kEk
Lemma
condabs (f , A) = kLf (A)k.
Nick Higham
Matrix Function Condition Numbers
12 / 33
Condition Number
kf (A + E) − f (A)k
.
→0 kEk≤
condabs (f , A) = lim sup
kLf (A)k := max
E6=0
kLf (A, E)k
.
kEk
Lemma
condabs (f , A) = kLf (A)k.
Scalar case: condabs (f , x) = |f 0 (x)|.
Nick Higham
Matrix Function Condition Numbers
12 / 33
Relative Condition Number
kf (A + E) − f (A)k
→0 kEk≤kAk
kf (A)k
condrel (f , A) := lim
sup
= condabs (f , A)
Nick Higham
kAk
.
kf (A)k
Matrix Function Condition Numbers
13 / 33
Computing Lf : Methods for Specific f
exponential Kenney & Laub (1998),
Al-Mohy & H (2009)
logarithm Al-Mohy, H & Relton (2013)
fractional power H & Lin (2013)
Latter three methods obtained by differentiating
alg for f .
Nick Higham
Matrix Function Condition Numbers
14 / 33
Computing Lf : Via 2n × 2n Matrix
Theorem (Mathias, 1996)
If f is 2n − 1 times ctsly diffble,
A E
f (A) Lf (A, E)
f
=
.
0 A
0
f (A)
Note that Lf (A, αE) = αLf (A, E), but α may effect alg
used for the evaluation.
Nick Higham
Matrix Function Condition Numbers
15 / 33
Computing Lf : Complex Step
Assume that f : Rn×n → Rn×n and A, E ∈ Rn×n . Then
f (A + i hE) − f (A) − i hLf (A, E) = o(h).
Thus (Al-Mohy & H, 2010)
f (A) ≈ Re f (A + i hE),
Lf (A, E) ≈ Im
f (A + i hE)
.
h
Errors O(h2 ).
h not restricted by fl pt arith considerations. Can take
h = 10−100 .
f alg must not employ complex arith.
Nick Higham
Matrix Function Condition Numbers
16 / 33
Condition Estimation (1)
Since Lf is a linear operator,
vec(Lf (A, E)) = Kf (A)vec(E)
where Kf (A) ∈ Cn
derivative.
2 ×n2
is the Kronecker form of the Fréchet
Can show
kLf (A)kF = kKf (A)k2 ,
kLf (A)k1
≤ kKf (A)k1 ≤ nkLf (A)k1 .
n
So problem reduces to matrix norm estimation.
Nick Higham
Matrix Function Condition Numbers
17 / 33
Condition Estimation (2)
Use the block 1-norm estimator of H & Tisseur (2000).
For kBk1 it needs Bx and B ∗ y for several x and y .
Let vec(X ) = x. Then Kf (A)x = vec(Lf (A, X )).
Theorem (H & Lin, 2013)
Let f ∈ C 2n−1 and and ef (z) := f (z).
If ef (A)∗ = ef (A∗ ) for all A ∈ Cn×n then
Kf (A)∗ x = vec(Lef (A, X ∗ )∗ ), where vec(X ) = x.
ef = f for most functions of interest.
ef = f ⇒ ef (A∗ )∗ = ef (A∗ ).
Nick Higham
Matrix Function Condition Numbers
18 / 33
Software for 1-Norm Estimator
MATLAB: normest1.
NAG library: F04YD, F04ZD, Mark 24.
Python:
Nick Higham
Matrix Function Condition Numbers
19 / 33
Python: SciPy 0.13.0
linalg/_expm_frechet.py
linalg/_expm_multiply.py
linalg/_sqrtm.py (blocked)
Nick Higham
Matrix Function Condition Numbers
20 / 33
Higher Derivatives
Needed for
level-2 condition number,
understanding accuracy of algorithms for computing
Lf (A, E).
(2)
Second Fréchet derivative Lf (A, E1 , E2 ) is unique
mutlilinear function of E1 , E2 satisfying
(2)
Lf (A+E2 , E1 ) − Lf (A, E1 ) − Lf
Nick Higham
(A, E1 , E2 ) = o(kE2 k).
Matrix Function Condition Numbers
21 / 33
Higher Derivatives
Needed for
level-2 condition number,
understanding accuracy of algorithms for computing
Lf (A, E).
(2)
Second Fréchet derivative Lf (A, E1 , E2 ) is unique
mutlilinear function of E1 , E2 satisfying
(2)
Lf (A+E2 , E1 ) − Lf (A, E1 ) − Lf
(A, E1 , E2 ) = o(kE2 k).
k th Fréchet derivative defined by
(k −1)
Lf
(k −1)
(A+Ek , E1 , . . . , Ek −1 ) − Lf
(k )
− Lf
Nick Higham
(A, E1 , . . . , Ek −1 )
(A, E1 , . . . , Ek ) = o(kEk k).
Matrix Function Condition Numbers
21 / 33
Existing Literature
Large literature on Fréchet derivatives in Banach
space.
Need specialized results for matrix functions.
Nick Higham
Matrix Function Condition Numbers
22 / 33
Existence & Continuity of Fréchet Derivatives
D = open subset of C.
Cn×n (D, p) = matrices with spectrum in D and largest
Jordan block of size p.
Theorem (H & Relton, 2013)
Let f be 2k p − 1 times continuously differentiable on D.
(k )
Then for A ∈ Cn×n (D, p) the kth Fréchet derivative Lf (A)
(k )
exists and Lf (A, E1 , . . . , Ek ) is continuous in A and
E1 , . . . , Ek ∈ Cn×n .
Proof uses Gâteaux derivative.
k = 1: Mathias (1996).
Nick Higham
Matrix Function Condition Numbers
23 / 33
Properties
Assume from now on conditions of theorem satisfied.
Then the Ei are interchangeable:
(2)
(2)
Lf (A, E1 , E2 ) = Lf (A, E2 , E1 ).
Indeed
(k )
Lf (A, E1 , . . . , Ek )
∂
f (A + s1 E1 + · · · + sk Ek ).
=
∂s1 · · · ∂sk s=0
Nick Higham
Matrix Function Condition Numbers
24 / 33
The Setup
Wish to estimate norms and condition numbers.
Knowing how to evaluate exact quantities helps
develop and test the estimates.
In practice we only ever work with n × n matrices.
Nick Higham
Matrix Function Condition Numbers
25 / 33
How to Compute Lf(2)
X1 =
Let
h
i
A E f (A) Lf (A,E)
1
.
Know
f
(X
)
=
.
1
0
f (A)
0 A


A
E
E
0
1
2
 0 A 0 E2 
0 1

X2 = I2 ⊗ X1 +
⊗ E2 = 
 0 0 A E1  .
0 0
0 0 0 A
Then

(2)
f (A) Lf (A, E1 ) Lf (A, E2 ) Lf (A, E1 , E2 )
 0
f (A)
0
Lf (A, E2 ) 
.
f (X2 ) = 
 0
0
f (A)
Lf (A, E1 ) 
0
0
0
f (A)

Nick Higham
Matrix Function Condition Numbers
26 / 33
)
How to Compute L(k
f
i
i
Define Xi ∈ C2 n×2 n by
0 1
Xi = I2 ⊗ Xi−1 +
⊗ I2i−1 ⊗ Ei ,
0 0
X0 = A.
Theorem (H & Relton, 2013)
(k )
The (1, n) block of f (Xk ) is Lf (A, E1 , . . . , Ek ).
Nick Higham
Matrix Function Condition Numbers
27 / 33
Level-2 Condition Number
“Condition number of the condition number”.
Demmel (1987) showed that for
matrix inversion (and D. J. Higham, 1995),
the eigenproblem,
polynomial zero-finding,
pole assignment in linear control problems
(relative) level-1 and level-2 cond no’s are equivalent.
Cheung & Cucker (2005) show same holds when
“condition number = 1 / distance to nearest ill-posed
problem”.
Nick Higham
Matrix Function Condition Numbers
28 / 33
Level-2 Condition Number
[2]
|condabs (f , A + Z ) − condabs (f , A)|
.
→0 kZ k≤
condabs (f , A) = lim sup
Theorem (H & Relton, 2013)
[2]
(2)
condabs (f , A) ≤ kKf (A)k2 ,
(2)
for a Kronecker matrix Kf (A) ∈ Cn
Nick Higham
4 ×n2
.
Matrix Function Condition Numbers
29 / 33
Level-2 Condition Number: Exponential
Theorem
Let A ∈ Cn×n be normal. Then in the 2-norm
[2]
condabs (exp, A) = condabs (exp, A).
Nick Higham
Matrix Function Condition Numbers
30 / 33
Level-2 Condition Number: Exponential
Theorem
Let A ∈ Cn×n be normal. Then in the 2-norm
[2]
condabs (exp, A) = condabs (exp, A).
Theorem
Let A ∈ Cn×n be normal. Then in the 2-norm
[2]
1 ≤ condrel (exp, A) ≤ 2condrel (exp, A) + 1.
Nick Higham
Matrix Function Condition Numbers
30 / 33
Level-2 Condition Number: Matrix Inverse
Theorem
For nonsingular A ∈ Cn×n ,
[2]
condabs (x −1 , A) = 2condabs (x −1 , A)3/2 .
Latter is what is expected from the scalar case:
f (x) = x −1 ⇒ |f 00 | = |2(f 0 )3/2 |.
Have experimental evidence that matrix case can be
similar to scalar case.
Nick Higham
Matrix Function Condition Numbers
31 / 33
Condition Number of the Fréchet Derivative
kLf (A + ∆A, E + ∆E) − Lf (A, E)k
.
→0 k∆Ak≤
condabs (Lf , A, E) = lim sup
k∆Ek≤
Theorem (H & Relton, 2013)
We have
condabs (f , A) ≤ condabs (Lf , A, E)
(2)
≤ max kLf (A, E, Z )k + condabs (f , A).
kZ k=1
We have developed method to estimate the bounds.
Nick Higham
Matrix Function Condition Numbers
32 / 33
Summary and Open Questions
Showed that given an O(n3 ) flop method for f :
(k )
Computing Lf (A, E) costs O(8k n3 ) flops.
(k )
Computing Kf (A) costs O(8k n3+2k ) flops.
Cost of estimates is O(n3 ) flops, given an O(n3 ) flop
method for Lf .
Open questions about relation between level-1 and
level-2 condition numbers.
How can we exploit the symmetry of
(k )
Lf (E1 , E2 , . . . , Ek ) in the Ei ?
Looking at componentwise condition numbers.
Nick Higham
Matrix Function Condition Numbers
33 / 33
References I
S. D. Ahipasaoglu, X. Li, and K. Natarajan.
A convex optimization approach for computing
correlated choice probabilities with many alternatives.
Preprint 4034, Optimization Online, 2013.
26 pp.
A. H. Al-Mohy and N. J. Higham.
Computing the Fréchet derivative of the matrix
exponential, with an application to condition number
estimation.
SIAM J. Matrix Anal. Appl., 30(4):1639–1657, 2009.
Nick Higham
Matrix Function Condition Numbers
1/8
References II
A. H. Al-Mohy and N. J. Higham.
The complex step approximation to the Fréchet
derivative of a matrix function.
Numer. Algorithms, 53(1):133–148, 2010.
A. H. Al-Mohy, N. J. Higham, and S. D. Relton.
Computing the Fréchet derivative of the matrix logarithm
and estimating the condition number.
SIAM J. Sci. Comput., 35(4):C394–C410, 2013.
S. Amat, S. Busquier, and J. M. Gutiérrez.
Geometric constructions of iterative functions to solve
nonlinear equations.
J. Comput. Appl. Math., 157(1):197–205, 2003.
Nick Higham
Matrix Function Condition Numbers
2/8
References III
J. Ashburner and G. R. Ridgway.
Symmetric diffeomorphic modelling of longitudinal
structural MRI.
Frontiers in Neuroscience, 6, 2013.
Article 197.
J. W. Demmel.
On condition numbers and the distance to the nearest
ill-posed problem.
Numer. Math., 51:251–289, 1987.
Nick Higham
Matrix Function Condition Numbers
3/8
References IV
B. García-Mora, C. Santamaría, G. Rubio, and J. L.
Pontones.
Computing survival functions of the sum of two
independent Markov processes. An application to
bladder carcinoma treatment.
Int. Journal of Computer Mathematics, 2013.
D. J. Higham.
Condition numbers and their condition numbers.
Linear Algebra Appl., 214:193–213, 1995.
Nick Higham
Matrix Function Condition Numbers
4/8
References V
N. J. Higham.
Functions of Matrices: Theory and Computation.
Society for Industrial and Applied Mathematics,
Philadelphia, PA, USA, 2008.
ISBN 978-0-898716-46-7.
xx+425 pp.
N. J. Higham and A. H. Al-Mohy.
Computing matrix functions.
Acta Numerica, 19:159–208, 2010.
N. J. Higham and L. Lin.
An improved Schur–Padé algorithm for fractional
powers of a matrix and their Fréchet derivatives.
SIAM J. Matrix Anal. Appl., 34(3):1341–1360, 2013.
Nick Higham
Matrix Function Condition Numbers
5/8
References VI
N. J. Higham and S. D. Relton.
The condition number of the Fréchet derivative of a
matrix function.
MIMS EPrint, Manchester Institute for Mathematical
Sciences, The University of Manchester, UK, 2013.
In preparation.
N. J. Higham and S. D. Relton.
Higher order Fréchet derivatives of matrix functions and
the level-2 condition number.
MIMS EPrint 2013.58, Manchester Institute for
Mathematical Sciences, The University of Manchester,
UK, Nov. 2013.
19 pp.
Nick Higham
Matrix Function Condition Numbers
6/8
References VII
B. Jeuris, R. Vandebril, and B. Vandereycken.
A survey and comparison of contemporary algorithms
for computing the matrix geometric mean.
Electron. Trans. Numer. Anal., 39:379–402, 2012.
C. S. Kenney and A. J. Laub.
A Schur–Fréchet algorithm for computing the logarithm
and exponential of a matrix.
SIAM J. Matrix Anal. Appl., 19(3):640–663, 1998.
R. Mathias.
A chain rule for matrix functions and applications.
SIAM J. Matrix Anal. Appl., 17(3):610–620, 1996.
Nick Higham
Matrix Function Condition Numbers
7/8
References VIII
D. Petersson.
A Nonlinear Optimization Approach to H2 -Optimal
Modeling and Control.
PhD thesis, Department of Electrical Engineering,
Linköping University, Sweden, SE-581 83 Linköping,
Sweden, 2013.
Dissertation No. 1528.
D. Petersson and J. Löfberg.
Model reduction using a frequency-limited H2 -cost.
Technical report, Department of Electrical Engineering,
Linköpings Universitet, SE-581 83, Linköping, Sweden,
Mar. 2012.
13 pp.
Nick Higham
Matrix Function Condition Numbers
8/8