DEMOINTVAL Interval computations in INTLAB
- How to define an interval I
- How to define an interval II
- How to define an interval III
- How to define an interval IV
- How to define an interval V
- Output formats of intervals I
- Output formats of intervals II
- Rigorous output
- What you see is correct
- Output formats of intervals III
- Changing interval output permanently
- Display with uncertainty
- Newton iteration
- Bit representation
- Invoking interval operations
- Interval matrix operations
- Sharp interval multiplication
- Fast interval multiplication
- Acceleration by vector/matrix notation
- Overestimation of interval operations
- Interval standard functions
- Complex interval standard functions
- Standard functions with argument out of range
- A common misuse of interval arithmetic
- Rigorous solution of linear systems
- Accuracy of rigorous linear system solving: Hilbert matrices
- Extremely ill-conditioned linear systems
- Verification of regularity of extremely ill-condtioned matrices
- Structured linear systems
- Sparse linear systems
- Inclusion of eigenvalues and eigenvectors
- Eigenvalue pairs and invariant subspaces
- Eigenvalues of structured matrices
- Nonlinear systems of equations, polynomials, etc.
- Enjoy INTLAB
A key to interval operations in INTLAB is changing the rounding mode. Following we ensure "rounding to nearest".
How to define an interval I
There are four possibilities to generate an interval, the first is a simple typecast of a real or complex quantity, for example a matrix. It uses Matlab conversion, i.e. the first component does not contain "2.3". This is because Matlab first converts "2.3" into binary format, and then the type cast is called.
format compact long infsup u = intval( [ 2.3 -4e1 ; 3 0 ] )
intval u = [ 2.29999999999999, 2.30000000000000] [ -40.00000000000000, -40.00000000000000] [ 3.00000000000000, 3.00000000000000] [ 0.00000000000000, 0.00000000000000]
How to define an interval II
The second possibility is to use INTLAB conversion of constants. In this case the argument is a string and INTLAB verified conversion to binary is called rather than Matlab conversion.
u = intval( '0.1 -3.1 5e2 .3e1' )
intval u = 1.0e+002 * [ 0.00099999999998, 0.00100000000001] [ -0.03100000000001, -0.03099999999998] [ 5.00000000000000, 5.00000000000000] [ 0.03000000000000, 0.03000000000000]
The first component, for example, definitely contains "0.1". Since 0.1 is not finitely representable in binary format, the radius of the first component must be nonzero.
ans = 1.0e-15 * 0.013877787807814 0.444089209850063 0 0
Generating an interval by an input string always produces a column vector. To change "u" into a 2 x 2 matrix, use
intval ans = 1.0e+002 * [ 0.00099999999998, 0.00100000000001] [ 5.00000000000000, 5.00000000000000] [ -0.03100000000001, -0.03099999999998] [ 0.03000000000000, 0.03000000000000]
Note that arrays are stored columnwise in Matlab.
How to define an interval III
The third possibility is by specification of midpoint and radius.
u = midrad( [ -3e1+2i ; .4711 ; 3 ] , 1e-4 )
intval u = [ -30.00010000000001 + 1.99989999999999i, -29.99989999999999 + 2.00010000000001i] [ 0.47099999999998 - 0.00010000000001i, 0.47120000000001 + 0.00010000000001i] [ 2.99989999999999 - 0.00010000000001i, 3.00010000000001 + 0.00010000000001i]
How to define an interval IV
The fourth possibility is by specification of infimum and supremum
u = infsup( [ 1 2 ] , [ 4 7 ] )
intval u = [ 1.00000000000000, 4.00000000000000] [ 2.00000000000000, 7.00000000000000]
How to define an interval V
Yet another possibility is to use [infimum,supremum] or midpoint,radius or giving significant digits in a string:
u = intval('[3,4] <5,.1> 1.23_')
intval u = [ 3.00000000000000, 4.00000000000000] [ 4.89999999999998, 5.10000000000001] [ 1.21999999999999, 1.24000000000001]
Output formats of intervals I
The output format can be changed using the different Matlab formats, for example
format midrad long e X = midrad( [ -3e1+2i ; .4711 ; 3 ] , 1e-4 )
intval X = < -3.000000000000000e+01 + 2.000000000000000e+00i, 1.000000000000001e-004> < 4.711000000000000e-01 + 0.000000000000000e+00i, 1.000000000001111e-004> < 3.000000000000000e+00 + 0.000000000000000e+00i, 1.000000000000001e-004>
format infsup long X
intval X = [ -30.00010000000001 + 1.99989999999999i, -29.99989999999999 + 2.00010000000001i] [ 0.47099999999998 - 0.00010000000001i, 0.47120000000001 + 0.00010000000001i] [ 2.99989999999999 - 0.00010000000001i, 3.00010000000001 + 0.00010000000001i]
or with a common exponent
intval ans = 1.0e+005 * [ -3.00001000000001 + 0.19998999999999i, -2.99998999999998 + 0.20001000000001i] [ 0.04709999999999 - 0.00001000000001i, 0.04712000000001 + 0.00001000000001i] [ 0.29998999999999 - 0.00001000000001i, 0.30001000000001 + 0.00001000000001i]
Output formats of intervals II
With two arguments the functions infsup and midrad define an interval, with one argument they control the output of an interval:
format short u = infsup( [ 1 2 ] , [ 4 7 ] ); infsup(u), midrad(u)
intval u = [ 1.0000, 4.0000] [ 2.0000, 7.0000] intval u = < 2.5000, 1.5000> < 4.5000, 2.5000>
Note that output in INTLAB is rigorous. That means,
left <= ans <= right for inf/sup notation ans in mid+/-rad for mid/rad notation
where ans is the true (real or complex) answer, and left,right, mid,rad are the numbers corresponding to the displayed decimal figures.
What you see is correct
format short midrad x = midrad(-5,1e-20)
intval x = < -5.0000, 0.0001>
One might get the impression that this is a very crude display. However, the radius of the interval is nonzero. Hence, in four decimal places, the displayed radius 0.0001 is the best possible. Changing the display format gives more information:
format long x
intval x = < -5.00000000000000, 0.00000000000001>
Output formats of intervals III
A convenient way to display narrow intervals is the following:
x=midrad(pi,1e-8); format short, infsup(x), midrad(x), disp_(x) format long, infsup(x), midrad(x), disp_(x) format short
intval x = [ 3.1415, 3.1416] intval x = < 3.1416, 0.0001> intval x = 3.1415 intval x = [ 3.14159264358979, 3.14159266358980] intval x = < 3.14159265358979, 0.00000001000001> intval x = 3.1415927_______
Mathematically the following is true: Form an interval of the displayed midpoint and a radius of 1 unit of the last displayed decimal figure, then this is a correct inclusion of the stored interval.
Changing interval output permanently
The interval output format can be changed permanently, for example, to infimum/supremum notation:
u = midrad( [ -3e1+2i ; .4711 ; 3 ] , 1e-4 ); format infsup u
intval u = [ -30.0002 + 1.9998i, -29.9998 + 2.0002i] [ 0.4709 - 0.0002i, 0.4713 + 0.0002i] [ 2.9998 - 0.0002i, 3.0002 + 0.0002i]
or to midpoint/radius notation:
format midrad u
intval u = < -30.0000 + 2.0000i, 0.0002> < 0.4711 + 0.0000i, 0.0002> < 3.0000 + 0.0000i, 0.0002>
or to display with uncertainties depicted by "_":
format _ u
intval u = -30.0000 + 2.0000i 0.4711 + 0.000_i 3.0000 + 0.000_i
Display with uncertainty
Display with uncertainty makes only sense for sufficiently narrow intervals. If the accuracy becomes too poor, INTLAB changes automatically to inf-sup or mid-rad display for real or complex intervals, respectively:
for k=-5:-1 disp_(midrad(pi,10^k)) end
intval = 3.1416 intval = 3.142_ intval = 3.14__ intval = 3.1___ [ 3.0415, 3.2416]
The following code is an interval Newton iteration to include sqrt(2).
format long _ f = @(x)(x^2-2); % Function f(x) = x^2-2 X = infsup(1.4,1.7); % Starting interval for i=1:4 xs = X.mid; % Midpoint of current interval Y = f(gradientinit(X)); % Gradient evaluation of f(X) X = intersect( X , xs-f(intval(xs))/Y.dx ) % Interval Newton step with intersection end
intval X = 1.4_____________ intval X = 1.4142__________ intval X = 1.414213562_____ intval X = 1.41421356237309
The "format _" output format shows nicely the quadratic convergence.
The last displayed result (which is in fact an interval) proves that the true value of sqrt(2) is between 1.41421356237308 and 1.41421356237310. Indeed, sqrt(2)=1.41421356237309504...
The format "long e" in Matlab displays the most figures. With this we see that the internal accuracy of the final X is in fact even better, the width is only 2 units in the last place.
format long e X format short
intval X = 1.414213562373095e+000
The bit pattern of floating-point numbers and of intervals as well is displayed as follows.
ans = 4×62 char array 'infimum ' ' +1.0110101000001001111001100110011111110011101111001011 * 2^0' 'supremum ' ' +1.0110101000001001111001100110011111110011101111001101 * 2^0'
That shows indeed that the left and right bound of X coincide up to 2 bits.
Invoking interval operations
An operation uses interval arithmetic if at least one of the operands is of type intval. For example, in
u = intval(5); y = 3*u-exp(u)
intval y = -133.4131
the result y is an inclusion of 15-exp(5). However, in
u = intval(5); z = 3*pi*u-exp(u)
intval z = -101.2892
the first multiplication "3*pi" is a floating point multiplication. Thus it is not guaranteed that the result z is an inclusion of 15pi-exp(5).
The following ensures a correct inclusion of the mathematical 3 \pi u - exp(u):
Z = 3*intval('pi')*u - exp(u) format long infsup(Z) format short
intval Z = -101.2892 intval Z = 1.0e+002 * [ -1.01289269298730, -1.01289269298729]
Interval matrix operations
INTLAB is designed to be fast. Case distinctions in interval multiplication can slow down computations significantly due to interpretation overhead. Therefore, there is a choice between 'fast' and 'sharp' evaluation of interval matrix products. This applies only to 'thick' intervals, i.e. intervals with nonzero diameter.
Sharp interval multiplication
In the following example, c is a real random matrix, C is an interval matrix with diameter zero (a thin interval matrix), and CC is an interval matrix with nonzero diameter (a thick interval matrix), all of dimension nxn for n=1000. First we measure the computing time with option 'SharpIVmult'.
n = 1000; c = randn(n); C = intval(c); C_ = midrad(c,.1); intvalinit('SharpIVmult') tic, scc = c*c; toc tic, sCC = C*C; toc tic, sCC = C*C_; toc tic, sCC__ = C_*C_; toc
===> Slow but sharp interval matrix multiplication in use Elapsed time is 0.053658 seconds. Elapsed time is 0.116197 seconds. Elapsed time is 0.177183 seconds. Elapsed time is 14.014845 seconds.
As can be seen, there is not much penalty if not both matrices are thick interval matrices; for both thick intervals, however, computation is slowed down significantly.
Fast interval multiplication
intvalinit('FastIVmult') tic, fcc = c*c; toc tic, fCC = C*C; toc tic, fCC = C*C_; toc tic, fCC__ = C_*C_; toc max(max(diam(fCC__)./diam(sCC__)))
===> Fast interval matrix multiplication in use (maximum overestimation factor 1.5 in radius, see FAQ)Elapsed time is 0.042636 seconds. Elapsed time is 0.106939 seconds. Elapsed time is 0.162284 seconds. Elapsed time is 0.212915 seconds. ans = 1.0615
As can be seen there is again not much penalty if not both matrices are thick. However, the 'fast' implementation is much faster than the 'sharp' at the cost of a little wider output. If intervals are very wide and any overestimation cannot be afforded (as in global optimization), the option 'SharpIVmult' may be used. It is shown in
S.M. Rump: Fast and parallel interval arithmetic. BIT Numerical Mathematics, 39(3):539-560, 1999
that the maximum (componentwise) overestimation by the option 'FastIVmult' compared to 'SharpIVmult' is a factor 1.5, for real and complex intervals. For not too thick intervals, however, it is shown that the overestimation is negligible.
Acceleration by vector/matrix notation
It is advisable to use vector/matrix notation when using interval operations. Consider
n = 10000; x = 1:n; y = intval(x); tic for i=1:n y(i) = y(i)^2 - y(i); end t1 = toc
t1 = 4.7558
This simple initialization takes considerable computing time. Compare to
tic y = intval(x); y = y.^2 - y; t2 = toc ratio = t1/t2
t2 = 0.0048 ratio = 991.7781
Sometimes code looks more complicated, a comment may help. It is worth it.
Overestimation of interval operations
Note that the main principle of interval arithmetic is that for given intervals A,B and an operation o, the result a o b is included in the interval result A o B for all a in A and all b in B. Since the result must be an interval, overestimations cannot be avoided in many situations. For example, in
close, kmax = 40; i = sqrt(-1); a=midrad(2,.7); b=midrad(1-i,1); plotintval(3*a-exp(b)); hold on phi = linspace(0,2*pi,kmax); [A,B] = meshgrid( mid(a)+rad(a)*exp(i*phi) , mid(b)+rad(b)*exp(i*phi) ); plot(3*A-exp(B)) hold off
the blue circle is the result of the interval operations, whereas the many circles approximate the power set operation (see also the DEMOINTLAB). Another reason for overestimation are dependencies, see below.
Interval standard functions
Interval standard functions in INTLAB are rigorous. For a given interval X and a function X let Y be the computed value of f(X). Then f(x) is in Y for all x in X. For example
x = intval(1e10); format long sin(x)
intval ans = -0.48750602508751
Note that the result is rigorous (try sin(2^1000) or similar). For timing comparison consider
format short n=10000; x=random(n,1); X=intval(x); tic, exp(x); tapprox = toc tic, exp(X); trigorous = toc ratio = trigorous/tapprox
tapprox = 0.0025 trigorous = 0.0061 ratio = 2.4648
Complex interval standard functions
Complex interval standard functions are rigorous as well, for example
format long Z = midrad(3+4i,1e-7); R = sin(Z)
intval R = 3.85374_________ - 27.01681_________i
It is mathematically correct, that sin(z) is an element of R for every complex z with Eucledian distance less than or equal to 1e-7 from 3+4i.
Standard functions with argument out of range
When entering a real argument leading to a complex value of a standard function, there are three possibilities to be specified by intvalinit:')
intvalinit('RealStdFctsExcptnNan'); sqrt(intval(-2)) intvalinit('RealStdFctsExcptnWarn'); sqrt(intval(-2)) intvalinit('RealStdFctsExcptnAuto'); sqrt(intval(-2))
===> Result NaN for real interval input out of range intval ans = NaN ===> Complex interval stdfct used automatically for real interval input out of range, but with warningintval ans = -0.00000000000000 + 1.41421356237309i ===> Complex interval stdfct used automatically for real interval input out of range (without warning)intval ans = -0.00000000000000 + 1.41421356237309i
A common misuse of interval arithmetic
The dependency problem is the most serious problem of (naive) interval arithmetic. The following procedure:
" Take some numerical algorithm and replace every operation by its corresponding interval operation. Then the computed interval result(s) definitely contain the true result which would be obtained without the presence of rounding errors. "
will most certainly fail in practice. Although a true statement (if no exception like divide by a zero interval occurs), the computed result interval(s) will, for very modest problem size, most certainly be of huge diameter and thus useless.
Consider, for example, the triangular matrix T where all elements on and below the diagonal are equal to 1, and take a randomly generated right hand side. The following lines do this for dimension n=50:
n = 50; T = tril(ones(n)); b = randn(n,1);
Then perform a standard forward substitution to compute an inclusion T\b. Note that X is defined to be an interval vector, so all operations are indeed interval operations (see above section "Invoking interval operations").
X = intval(zeros(n,1)); for i=1:n X(i) = b(i) - T(i,1:i-1)*X(1:i-1); end X
intval X = 1.54496886565095 -1.73641760729346 0.92018854401184 0.65197328202597 -2.33827755592769 0.16629614976347 2.43421718014301 -1.1757752311920_ -0.3244045981527_ -0.3109015016353_ 0.8243081711104_ 0.062040815522__ -0.257822709060__ -1.33661157686___ 1.18458604457___ -0.33469668576___ 0.11636807926___ -0.1864900273____ 1.0265758952____ -1.1168262990____ 0.3330179902____ 0.397103103_____ -2.959729555_____ 3.799835771_____ -0.87738401______ -2.15995104______ 2.00125210______ -0.7093691_______ 0.1979432_______ 0.1816328_______ -0.793923________ 1.278688________ -1.115279________ 0.108663________ 2.59745_________ -1.95887_________ 1.58669_________ 0.7210__________ -1.4453__________ -1.4330__________ 0.951___________ -1.588___________ 1.736___________ -2.084___________ 1.10____________ 0.06____________ -1.63____________ 1.2_____________ 0.6_____________ 0.4_____________
The result is displayed with uncertainty, perfectly making visible the exponential loss of accuracy. This is due to one of the most common misuses of interval arithmetic, also called "naive interval arithmetic". For more details and examples cf.
S.M. Rump: Verification methods: Rigorous results using floating-point arithmetic. Acta Numerica, 19:287-449, 2010.
to be downloaded from "www.ti3.tuhh.de/rump". See, in particular,
Chapter 10.1: The failure of the naive approach: interval Gaussian elimination (IGA)
Note that the linear system above is very well-conditioned:
ans = 64.270085531579554
By the well-known rule of thumb of numerical analysis we expect at least 14 correct digits in a floating-point approximation T\b. Using a proper (non-naive) method, an inclusion of this quality is indeed achieved:
intval ans = 1.54496886565095 -1.73641760729346 0.92018854401184 0.65197328202597 -2.33827755592769 0.16629614976346 2.43421718014301 -1.17577523119201 -0.32440459815271 -0.31090150163526 0.82430817111040 0.06204081552205 -0.25782270906023 -1.33661157686160 1.18458604456945 -0.33469668576355 0.11636807926340 -0.18649002731532 1.02657589515373 -1.11682629899407 0.33301799020482 0.39710310294512 -2.95972955511010 3.79983577121854 -0.87738400810208 -2.15995104006627 2.00125209821602 -0.70936905983630 0.19794316597518 0.18163284792087 -0.79392266302791 1.27868838309226 -1.11527942977767 0.10866280278142 2.59745200007108 -1.95886795922010 1.58668718772608 0.72102363349527 -1.44532990095592 -1.43303475689456 0.95128027871488 -1.58813522029171 1.73628335123737 -2.08378329283865 1.10008919761228 0.05734129838531 -1.62759197715764 1.24373722971505 0.59753040062128 0.38015759066768
Such methods are called "self-validating methods" or "verification methods". For some details see the reference above or
S.M. Rump: Self-validating methods. Linear Algebra and its Applications (LAA), 324:3-13, 2001.
Due to an improved evaluation of the residual (default option "intvalinit('ImprovedResidual')" , see also function "lssresidual.m") 15 correct decimal digits of the result are computed.
Rigorous solution of linear systems
The INTLAB linear system solver can be called with "\" or "verifylss".' For example, [bare with me, I am often in Japan where the backslash appears like japanese Yen.]
n = 100; A = randn(n); b = A*ones(n,1); X = verifylss(A,b);
generates and solves a randomly generated 100x100 linear system. The inclusion and its quality is checked by
X([1:3 98:100]) max( X.rad ./ abs(X.mid) )
intval ans = 1.00000000000000 0.99999999999999 0.99999999999999 0.99999999999999 1.00000000000000 0.99999999999999 ans = 5.551115123125781e-16
which calculates the maximum relative error of the inclusion radius with respect to the midpoint. The same is done by
ans = 4.440892098500623e-16
Accuracy of rigorous linear system solving: Hilbert matrices
For estimating accuracy, try
format long e n = 10; H = hilb(n); b = ones(n,1); X = verifylss(H,b)
intval X = -9.99830__________e+000 9.89853__________e+002 -2.375688_________e+004 2.402116_________e+005 -1.261125_________e+006 3.783408_________e+006 -6.726110_________e+006 7.000691_________e+006 -3.937911_________e+006 9.237120_________e+005
The notoriously ill-conditioned Hilbert matrix is given by H_ij := 1/(i+j-1). To estimate the accuracy, we use the symbolic toolbox to compute the perturbation of the solution when perturbing only the (7,7)-element of the input matrix by 2^(-52):
Hs = sym(H,'f'); Hs(7,7) = Hs(7,7)*(1+sym(2^(-52))); double( Hs \ b )
ans = -9.997198050336207e+00 9.897576960772726e+02 -2.375483680969619e+04 2.401930526008937e+05 -1.261036061017173e+06 3.783164425266260e+06 -6.725710140930751e+06 7.000304229698808e+06 -3.937707813532462e+06 9.236673822756851e+05
The statement "sym(H,'f')" makes sure that no conversion error occurs when changing H into symbolic format.
Extremely ill-conditioned linear systems
By default, all computations in INTLAB are, like in Matlab, performed in double precision. This restricts treatable linear systems to a maximum condition number of roughly below 10^16.
Starting with INTLAB Version 7, I rewrote my linear system solver completely. Now, although only double precision is used, linear systems with larger condition numbers are solvable. Consider
format long _ n = 20; A = invhilb(n); condA = norm(double(inv(sym(A))))*norm(A)
condA = 5.478702657981728e+27
The common rule of thumb tells that for a condition number 10^k, an algorithm in double precision should produce 16-k correct digits. In our case this means roughly 16-27=-11 "correct" digits, namely none. For a random right hand side Matlab computes
b = A*randn(n,1); x = A\b
x = 1.0e+05 * -3.526561759042908 -1.454567045706571 -0.592520417296815 -0.235078095024219 -0.089144109711128 -0.031474527357279 -0.009882856050908 -0.002503614144822 -0.000398423032250 0.000031205126401 0.000038665477879 0.000006113262176 -0.000029375253182 0.000009822814340 0.000032598385605 0.000010223756347 0.000008182354189 -0.000003888579149 0.000057657456205 0.000191873032237
A corresponding warning indicates the difficulty of the problem. Note that in this case the Matlab guess of the condition number is pretty good.
A verified inclusion of the solution is computed by
X = verifylss(A,b,'illco')
intval X = 1.0e+005 * -3.79143_________ -1.65610894______ -0.709759661_____ -0.2950853244____ -0.1169928675649_ -0.04311353750180 -0.01409374776675 -0.00368058769965 -0.00054223078157 0.00011892565814 0.00010397698909 0.00001270523760 -0.00005403043476 -0.00001774136504 0.00001222576888 -0.00001114755428 -0.00003401120403 -0.00009069689194 -0.00009515901791 -0.00004152613564
As expected the Matlab approximation differs significantly from the true values, for some components the sign is incorrect. The maximum relative error of the components of the computed inclusion is
ans = 3.511080212154039e-07
so that each component is accurately included with at least 7 correct figures.
Verification of regularity of extremely ill-condtioned matrices
For accurate dot product computations some reference implementations are used which suffer from interpretation overhead. This makes it time consuming to solve linear systems with extremely ill-conditioned matrix.
However, checking regularity of very extremely ill-conditioned matrices can be done very fast using mod-p arithmetic. Consider
n = 50; cnd = 1e200; digits(400) A = randmat(n,cnd); tic, Cnd = norm(A,inf)*norm(double(inv(sym(A,'f'))),inf), toc
Cnd = 7.086324698036915e+226 Elapsed time is 4.833116 seconds.
As can be seen the matrix A is indeed extremely ill-conditioned. However, checking regularity in mod-p arithmetic is rigorous and very fast:
tic, isregular(A), toc
ans = 1 Elapsed time is 0.010447 seconds.
With a chance of about 1/p regularity cannot be checked. Since p is of the order 1e8, this chance is slim. It can be further diminished by using several primes, see "isregular".
Note that verification of regularity, i.e. a sufficient criterion, is simple, whereas a necessary and sufficient check of singularity is much more involved.
Structured linear systems
In general, intervals are treated as independent quantities. If, however, there are dependencies, then taking them into account may shrink the solution set significantly. An example is
format short n = 4; e = 1e-3; intvalinit('displayinfsup'); A = midrad( toeplitz([0 1+e e 1+e]),1e-4); b = A.mid*ones(n,1); Amid = A.mid X1 = verifylss(A,b)
===> Default display of intervals by infimum/supremum (e.g. [ 3.14 , 3.15 ]) Amid = 0 1.0010 0.0010 1.0010 1.0010 0 1.0010 0.0010 0.0010 1.0010 0 1.0010 1.0010 0.0010 1.0010 0 intval X1 = [ 0.3327, 1.6673] [ 0.3327, 1.6673] [ 0.3327, 1.6673] [ 0.3327, 1.6673]
First the matrix has been treated as a general interval matrix without dependencies. Recall that only the midpoint is displayed above; all entries of the interval matrix have a uniform tolerance of 1e-4.
Several structural information may be applied to the input matrix, for example,
X2 = verifystructlss(structure(A,'symmetric'),b); X3 = verifystructlss(structure(A,'symmetricToeplitz'),b); X4 = verifystructlss(structure(A,'generalToeplitz'),b); X5 = verifystructlss(structure(A,'persymmetric'),b); X6 = verifystructlss(structure(A,'circulant'),b); res = [ X1 X2 X3 X4 X5 X6 ]; rad(res)
ans = 0.6672 0.5062 0.1689 0.3376 0.5062 0.0003 0.6672 0.5062 0.1688 0.3375 0.5062 0.0003 0.6672 0.5062 0.1688 0.3375 0.5062 0.0003 0.6672 0.5062 0.1689 0.3376 0.5062 0.0003
Here only the radii of the inclusions are displayed. Note that the inclusion may become much narrower, in particular treating the input data as a circulant matrix.
Sparse linear systems
The following generates a random sparse system with about 9 nonzero elements per row.
format short n = 10000; A = sprand(n,n,2/n)+speye(n); A = A*A'; b = ones(n,1);
The linear system is generated to be symmetric positive definite. Before calling the verified linear system solver, the fill-in should reduced. The original matrix looks like
spy(A) title('sparsity pattern of A')
whereas after minimum degree reordering the matrix looks like
p = symamd(A); spy(A(p,p)) title('sparsity pattern of renumbered A')
The timing for the built-in (floating point) solver compared to the verified solver is as follows:
tic, x = A(p,p)\b(p); toc
Elapsed time is 0.525967 seconds.
tic, X = verifylss(A(p,p),b(p)); toc
Elapsed time is 2.297256 seconds.
Inclusion of eigenvalues and eigenvectors
To compute verified inclusions of eigenvalue/eigenvector pairs of simple or multiple eigenvalues, consider, for example, the famous Wilkinson(21) matrix. Following inclusions of the last four eigenvalue/eigenvector pairs are displayed. Those are the most ill-conditioned.
format long _ W = wilkinson(21); % generation of the matrix [V,D] = eig(W); % eigenvalue/eigenvector approximations for k=18:21 [L,X] = verifyeig(W,D(k,k),V(:,k)) % inclusions for the small eigenvalues end
intval L = 9.21067864730491 intval X = -0.38247_________ 0.30189_________ 0.44607_________ 0.23816_________ 0.080419________ 0.020042________ 0.0039719_______ 0.00065440______ 0.000092313_____ 0.0000112430____ -0.00000000000___ -0.000011243042__ -0.000092313004__ -0.0006543963623_ -0.00397193251082 -0.02004206756034 -0.08041877341334 -0.23815677108033 -0.44606931512504 -0.30188982395948 0.38246757537813 intval L = 9.21067864736133 intval X = 0.38246757533800 -0.30188982390622 -0.44606931509072 -0.23815677111720 -0.08041877354260 -0.02004206794300 -0.00397193399399 -0.0006544037082_ -0.000092357143__ -0.00001155397___ -0.00000250882___ -0.0000115540____ -0.000092357_____ -0.00065440______ -0.0039719_______ -0.020042________ -0.080419________ -0.23816_________ -0.44607_________ -0.30189_________ 0.38247_________ intval L = 10.74619418290332 intval X = 0.55939864230126 0.41742001280922 0.16949775589362 0.04805373844102 0.01052087952088 0.00188039874002 0.0002842567806_ 0.000037252700__ 0.00000430986___ 0.0000004422____ -0.000000000_____ -0.00000044______ -0.0000043_______ -0.000037________ -0.00028_________ -0.00188_________ -0.0105__________ -0.048___________ -0.170___________ -0.42____________ -0.56____________ intval L = 10.74619418290339 intval X = 0.56____________ 0.42____________ 0.170___________ 0.048___________ 0.0105__________ 0.00188_________ 0.00028_________ 0.000037________ 0.0000043_______ 0.00000045______ 0.000000084_____ 0.0000004509____ 0.00000431088___ 0.000037252833__ 0.0002842568009_ 0.00188039874370 0.01052087952171 0.04805373844126 0.16949775589369 0.41742001280919 0.55939864230118
Eigenvalue pairs and invariant subspaces
The largest eigenvlues are 10.74619418290332 and 10.74619418290339, where all displayed digits are verified to be correct. Invariant subspaces of nearby eigenvalues are in general ill-conditioned. Nearby eigenvalues can also be regarded as clusters. From the inclusions above we can judge how narrow the eigenvalues are. So one of the approximations can be used as an approximation of the pair.
[L,X] = verifyeig(W,D(18,18),V(:,18:19)) % inclusion of the 18/19 eigenvalue pair [L,X] = verifyeig(W,D(20,20),V(:,20:21)) % inclusion of the 20/21 eigenvalue pair
intval L = 9.2106786473____ + 0.0000000000____i intval X = -0.38246721566493 0.38246757533800 0.30188954003016 -0.30188982390622 0.44606889559399 -0.44606931509072 0.23815654709237 -0.23815677111720 0.08041869777897 -0.08041877354260 0.02004204871064 -0.02004206794300 0.00397192877519 -0.00397193399399 0.00065439574687 -0.00065440370821 0.00009231291679 -0.00009235714338 0.00001124303112 -0.00001155397350 -0.00000000000118 -0.00000250882017 -0.00001124304199 -0.00001155396293 -0.00009231300365 -0.00009235705656 -0.00065439636234 -0.00065440309275 -0.00397193251082 -0.00397193025836 -0.02004206756034 -0.02004204909331 -0.08041877341334 -0.08041869790822 -0.23815677108033 -0.23815654712923 -0.44606931512504 -0.44606889555967 -0.30188982395948 -0.30188953997691 0.38246757537813 0.38246721562480 intval L = 10.7461941829034_ + 0.0000000000000_i intval X = 0.55939864230126 0.53926516857120 0.41742001280922 0.40239653183024 0.16949775589362 0.16339731453128 0.04805373844102 0.04632422283758 0.01052087952089 0.01014221959041 0.00188039874009 0.00181272078417 0.00028425678094 0.00027402603484 0.00003725270198 0.00003591205802 0.00000430988244 0.00000415574012 0.00000044236679 0.00000043485205 0.00000000151024 0.00000008241241 -0.00000042613744 0.00000045076775 -0.00000415472850 0.00000431085765 -0.00003591192480 0.00003725283040 -0.00027402601454 0.00028425680050 -0.00181272078050 0.00188039874363 -0.01014221958959 0.01052087952169 -0.04632422283735 0.04805373844125 -0.16339731453120 0.16949775589369 -0.40239653183027 0.41742001280919 -0.53926516857129 0.55939864230118
Note that interval output with uncertainty ("_") is used, so all displayed decimal places of the bases of the invariant subspaces are verified to be correct. As explained in section "Output formats of intervals III", the inclusion 10.7461941829034_ of the two smallest eigenvlues reads [10.7461941829033,10.7461941829035], thus including the true eigenvalues as displayed above.
The mathematical statement is that the displayed intervals for the cluster contain (at least) two eigenvalues of the Wilkinson matrix W. The size of the cluster is determined by the number of columns of the invariant subspace approximation.
Eigenvalues of structured matrices
As for linear systems, the interval input matrix may be structured. Taking into account such structure information may shrink the inclusion. As an example consider
format short e = 1e-3; A = midrad( toeplitz([0 1+e -e/2 1+e]),1e-4); [v,d] = eig(A.mid); xs = v(:,2:3); lambda = d(2,2); X1 = verifyeig(A,lambda,xs); X2 = verifystructeig(structure(A,'symmetric'),lambda,xs); X3 = verifystructeig(structure(A,'symmetricToeplitz'),lambda,xs); X4 = verifystructeig(structure(A,'generalToeplitz'),lambda,xs); X5 = verifystructeig(structure(A,'persymmetric'),lambda,xs); X6 = verifystructeig(structure(A,'circulant'),lambda,xs); res = [ X1 X2 X3 X4 X5 X6 ]; rad(res)
ans = 1.0e-03 * 0.7401 0.6455 0.3378 0.4850 0.5700 0.4000
As for linear systems, only the radii of the inclusions are displayed.
Nonlinear systems of equations, polynomials, etc.
For inclusions of systems of nonlinear equations, of roots of polynomials etc. cf. the corresponding demos.
INTLAB was designed and written by S.M. Rump, head of the Institute for Reliable Computing, Hamburg University of Technology. Suggestions are always welcome to rump (at) tuhh.de