DEMOINTVAL Interval computations in INTLAB

Contents

A key to interval operations in INTLAB is changing the rounding mode. Following we ensure "rounding to nearest".

setround(0)

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.

u.rad

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

reshape(u,2,2)

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.89999999999999,   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>

or

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

1e4*X

intval ans = 
  1.0e+005 *
[  -3.00001000000001 +  0.19998999999999i,  -2.99998999999999 +  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>

Rigorous output

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

Consider

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]

Newton iteration

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

Bit representation

The bit pattern of floating-point numbers and of intervals as well is displayed as follows.

getbits(X)

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.016191 seconds.
Elapsed time is 0.068021 seconds.
Elapsed time is 0.092340 seconds.
Elapsed time is 11.712556 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.015615 seconds.
Elapsed time is 0.053699 seconds.
Elapsed time is 0.066320 seconds.
Elapsed time is 0.118413 seconds.
ans =
    1.0620

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 =
    3.5986

This simple initialization takes considerable computing time. Compare to

tic
y = intval(x);
y = y.^2 - y;
t2 = toc
ratio = t1/t2

t2 =
    0.0020
ratio =
   1.7740e+03

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
Overestimation of interval operations

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 =
   4.2718e-04
trigorous =
    0.0175
ratio =
   40.8785

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 Euclidian 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 = 
   0.39687953198958
  -0.56740197138354
   1.10082187113049
  -0.35321774537699
  -1.81406407142539
  -0.43453751398342
   2.32650849408573
  -1.2709359097054_
   1.3918810983008_
  -2.0529312207357_
   2.0252875245986_
  -0.809903571216__
  -0.427606984178__
  -1.508668209095__
   3.49383756426___
  -1.55206781379___
   0.97739941047___
  -1.5213393432____
   1.4744584536____
   0.2361192932____
  -0.9590255806____
  -2.878398971_____
   3.150872745_____
  -0.06936361______
  -0.32911216______
  -0.34675168______
   1.65106065______
  -0.5935438_______
  -2.4215111_______
   1.6294674_______
   0.509579________
   0.506362________
  -0.616269________
  -0.965903________
   0.89245_________
   0.13974_________
  -0.62716_________
  -1.0949__________
   3.8834__________
  -4.2718__________
   1.346___________
   0.347___________
   1.388___________
  -3.02____________
   1.58____________
   0.99____________
  -1.27____________
   0.1_____________
  -0.4_____________
   2.2_____________

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:

cond(T)

ans =
  64.270085531579468

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:

verifylss(T,b)

intval ans = 
   0.39687953198958
  -0.56740197138354
   1.10082187113049
  -0.35321774537699
  -1.81406407142539
  -0.43453751398341
   2.32650849408572
  -1.27093590970542
   1.39188109830083
  -2.05293122073570
   2.02528752459860
  -0.80990357121638
  -0.42760698417779
  -1.50866820909509
   3.49383756425849
  -1.55206781379031
   0.97739941047032
  -1.52133934324254
   1.47445845356984
   0.23611929315358
  -0.95902558057949
  -2.87839897051334
   3.15087274452804
  -0.06936360756385
  -0.32911216379048
  -0.34675167739015
   1.65106065205435
  -0.59354382560490
  -2.42151111470462
   1.62946737929466
   0.50957942840624
   0.50636176421430
  -0.61626872221439
  -0.96590273903046
   0.89244969900543
   0.13974069577751
  -0.62716328488942
  -1.09485338033753
   3.88337114287061
  -4.27178901510237
   1.34554963887067
   0.34676468266117
   1.38804080662359
  -3.02446838129360
   1.58368647718130
   0.98666804610288
  -1.27479588612223
   0.02745187799245
  -0.40714027748311
   2.17732447432640

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
   1.00000000000000
   1.00000000000000
   0.99999999999999
   0.99999999999999
ans =
     2.220446049250328e-16

which calculates the maximum relative error of the inclusion radius with respect to the midpoint. The same is done by

max(relerr(X))

ans =
     2.220446049250328e-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.478702657981727e+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+04 *
   7.321109455577261
   3.025159198731216
   1.268943417919181
   0.533418768270654
   0.222293354099652
   0.090898725948031
   0.035978549642143
   0.013486319482052
   0.004729604952212
   0.001663101082006
   0.000341823635279
   0.000049611224054
  -0.000181764854093
  -0.000062232245592
   0.000018163312053
   0.000161586376965
   0.000253435588127
  -0.000295760199754
  -0.000408560294806
  -0.000987709900885

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+004 *
   9.3166__________
   3.962943________
   1.70719503______
   0.73546401752___
   0.313101335106__
   0.1300292967384_
   0.0517511577165_
   0.01919515332241
   0.00644083405422
   0.00200080149058
   0.00033324068084
   0.00002031653199
  -0.00018113850168
  -0.00006028106191
  -0.00000359259829
   0.00012354286209
   0.00024336104730
  -0.00020241879051
  -0.00011392393837
  -0.00038021211346

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

max(relerr(X))

ans =
     3.390173344223075e-06

so that each component is accurately included with at least 7 correct figures.

Verification of regularity of extremely ill-conditioned 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 =
    5.393490108579374e+232
Elapsed time is 3.739011 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.039490 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')
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')
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.198443 seconds.

tic, X = verifylss(A(p,p),b(p)); toc

Elapsed time is 1.711355 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.301890________
   0.44607_________
   0.238157________
   0.080419________
   0.0200421_______
   0.00397193______
   0.00065440______
   0.000092313_____
   0.0000112430____
  -0.00000000000___
  -0.000011243042__
  -0.0000923130036_
  -0.00065439636234
  -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.00065440370822
  -0.0000923571434_
  -0.000011553974__
  -0.00000250882___
  -0.0000115540____
  -0.000092357_____
  -0.00065440______
  -0.0039719_______
  -0.0200421_______
  -0.080419________
  -0.238157________
  -0.44607_________
  -0.301890________
   0.38247_________
intval L = 
  10.74619418290332
intval X = 
   0.55939864230126
   0.41742001280922
   0.16949775589362
   0.04805373844102
   0.01052087952088
   0.00188039874002
   0.0002842567806_
   0.0000372526996_
   0.00000430986___
   0.00000044221___
   0.0000000000____
  -0.000000442_____
  -0.00000431______
  -0.0000372_______
  -0.000284________
  -0.00188_________
  -0.0105__________
  -0.048___________
  -0.169___________
  -0.417___________
  -0.56____________
intval L = 
  10.74619418290339
intval X = 
   0.56____________
   0.42____________
   0.169___________
   0.048___________
   0.0105__________
   0.00188_________
   0.000284________
   0.000037________
   0.00000431______
   0.000000451_____
   0.0000000839____
   0.0000004509____
   0.000004310876__
   0.0000372528328_
   0.0002842568009_
   0.00188039874370
   0.01052087952170
   0.04805373844125
   0.16949775589369
   0.41742001280919
   0.55939864230118

Eigenvalue pairs and invariant subspaces

The largest eigenvalues 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 eigenvalues 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.4173    0.4170    0.3042    0.4043    0.4086    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.

Enjoy INTLAB

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


Published with MATLAB® R2020a