# Utility routines

## Contents

- Gaussian elimination using naive interval arithmetic
- Gaussian elimination using affine interval arithmetic
- Gaussian elimination using gradients
- Gaussian elimination using hessians
- Gaussian elimination using k-bit arithmetic
- Gaussian elimination over GF[p]
- Random integers
- Even/odd
- Ill-conditioned matrices
- The parity of a permutation
- Positive definiteness of a matrix
- Prove of regularity of matrix with arbitrary condition number
- Full rank of a matrix
- Random orthogonal matrices
- Relative errors
- Distance in bits and the accuracy of inclusion functions
- The bit pattern
- Minimum precision to store data
- Intelligent help
- Greatest common divisor and least common multiple
- All factors of an integer and perfect numbers
- Symbolic computations with pretty output
- Checking exponential many possibilities
- Traversing all vertices of an interval matrix
- Enjoy INTLAB

## Gaussian elimination using naive interval arithmetic

The generic routine solvewpp (Solve linear system With Partial Pivoting) works for a number of data types. It is based on Gaussian elimination with subsequent forward and backward substitution using corresponding routines luwpp, forward and backward.

The first example is using ordinary interval arithmetic. The entries of the matrix are afflicted with a relative error of 1e-10.

format long n = 5; a = randn(n).*midrad(1,1e-10); b = a*ones(n,1); format long xint = solvewpp(a,b)

intval xint = 1.0000000_______ 1.0000000_______ 1.0000000_______ 1.00000000______ 1.0000000_______

## Gaussian elimination using affine interval arithmetic

Next we check whether the result is narrower using affine arithmetic.

xaff = solvewpp(affari(a),b) [ rad(xint) rad(xaff) ]

affari xaff = 1.00000000______ 1.00000000______ 1.00000000______ 1.000000000_____ 1.0000000_______ ans = 1.0e-07 * 0.170135415755723 0.043519243764223 0.110010376364045 0.064242825503769 0.125490400293415 0.068209489167970 0.031334421723983 0.009750474871950 0.192322044956228 0.125225998459655

## Gaussian elimination using gradients

We may solve the linear system with gradients as well:

xgra = solvewpp(gradientinit(a),b); xgra.x

intval ans = 1.0000000_______ 1.0000000_______ 1.0000000_______ 1.00000000______ 1.0000000_______

The data can be used to analyze the sensitivity of individual solution components with respect to individual matrix components. Suppose we want to know how the 3rd solution component depends on the upper left matrix entry a(1,1). The answer is

sens = xgra(3).dx(1)

intval sens = (1,1) 1.0413362_______

This can be verified by direct computation using the midpoint equation:

xs = a.mid\b.mid; e = 1e-10; aa = a.mid; aa(1,1) = aa(1,1)+e; xss = aa\b.mid; sens_ = ( xss(3)-xs(3) ) / e

sens_ = 1.041333685947166

## Gaussian elimination using hessians

The linear system can be solved with hessians as well, so that information on all partial second derivatives are available. We only display the second derivative of the 3rd solution components with respect to the upper left matrix element a(1,1).

```
format short
xhes = solvewpp(hessianinit(a),b);
struct(xhes)
sens_hx = xhes(3).hx(1)
```

ans = struct with fields: x: [5×1 intval] dx: [25×5 intval] hx: [625×5 intval] intval sens_hx = (1,1) -0.9368

## Gaussian elimination using k-bit arithmetic

INTLAB's fl-package uses k-bit IEEE 754 binary arithmetic. After initialization, in this example to 22 bits and exponent range -99 ... 100, Gaussian elimination is performed in this precision.

```
flinit(22,100)
flinit('DisplayDouble')
xflk = solvewpp(fl(a),b)
```

ans = 'Initialization of fl-format to 22 mantissa bits incl. impl. 1 and exponent range -99 .. 100 for normalized fl-numbers' ===> Display fl-variables as doubles fl-type intval xflk = 1.0000 1.0000 1.0000 1.0000 1.0000

Often it is convenient to display the bit pattern of the answer. The format can be changed by

```
format flbits
xflk
```

===> Display fl-variables by bit representation fl-type intval xflk = [ +1.111111111111100100101 * 2^-1 , +1.000000000000001101100 * 2^0 ] [ +1.111111111111101110101 * 2^-1 , +1.000000000000001001001 * 2^0 ] [ +1.111111111111101100000 * 2^-1 , +1.000000000000001010101 * 2^0 ] [ +1.111111111111111010101 * 2^-1 , +1.000000000000000010100 * 2^0 ] [ +1.111111111111100000010 * 2^-1 , +1.000000000000001110110 * 2^0 ]

## Gaussian elimination over GF[p]

INTLAB's gfp-package implements a finite field GF[p] for some prime p. Floating-point numbers are rational and can thus be regarded as elements of a finite field. An example is the inverse over GF[17]:

gfpinit(17) n = 5; A = gfp(hilbert(n)) R = solvewpp(A,eye(n))

Galois field prime set to 17 ans = 17 gfp A = 4 2 7 1 11 2 7 1 11 12 7 1 11 12 3 1 11 12 3 9 11 12 3 9 8 gfp R = 2 10 16 7 13 10 10 1 5 12 16 1 6 10 14 7 5 10 5 16 13 12 14 16 9

Here hilbert(n) denotes the Hilbert matrix multiplied by the least common multiple of all denominators. Correctness of the computed inverse is verified by

R*A

gfp ans = 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1

## Random integers

Sometimes random integers are needed. Here is how to generate an m x n matrix of integers in the range [1 ... k] or [ -k ... k] .

m = 2; n = 7; k = 10; ipos = randint(k,m,n) igen = randint(-k,m,n)

ipos = 10 2 1 1 4 9 3 6 2 2 10 7 8 2 igen = -7 9 8 -10 3 2 10 -5 -6 -5 -10 -4 3 -3

## Even/odd

Sometimes it may be useful to check the parity of an integer. As the name suggests, odd(n) is 1 for odd n. It may be applied to a vector as well.

N = randint(-10,1,5) Nodd = odd(N) Neven = even(N) Nodd | Neven

N = 6 -8 6 -2 9 Nodd = 1×5 logical array 0 0 0 0 1 Neven = 1×5 logical array 1 1 1 1 0 ans = 1×5 logical array 1 1 1 1 1

## Ill-conditioned matrices

The routine randmat generates extremely ill-conditioned matrices.

n = 50; cnd = 1e25; A = randmat(n,cnd); c = cond(A)

c = 4.8279e+18

Of course, the Matlab routine cond cannot approximate the true condition number of about cnd=1e25. The true condition number can be computed by INTLAB's especially designed inclusion algorithm for extreme condition numbers (up to about 1e30).

```
C = cond(intval(A),'illco')
```

intval C = 1.0e+024 * 9.6694

Note that C is an inclusion of the true 2-norm condition number.

## The parity of a permutation

Mathematically, independent of the matrix B, the matrix A computed by

n = 5; B = boothroyd(n)/10; A = B'*B

A = 1.0e+03 * 0.2225 0.7292 0.9277 0.5373 0.1187 0.7292 2.3930 3.0465 1.7650 0.3902 0.9277 3.0465 3.8800 2.2485 0.4973 0.5373 1.7650 2.2485 1.3033 0.2883 0.1187 0.3902 0.4973 0.2883 0.0638

must be positive semi-definite. However, due to rounding errors it might not. The method of choice to factor A would be Cholesky decomposition. However, Gaussian elemination may produce a permutation although mathematically not necessary:

```
[L,U,p] = lu(A,'vector')
```

L = 1.0000 0 0 0 0 0.7860 1.0000 0 0 0 0.2398 0.8511 1.0000 0 0 0.5791 -0.4832 0.3932 1.0000 0 0.1280 -0.1829 0.1851 0.5220 1.0000 U = 1.0e+03 * 0.9277 3.0465 3.8800 2.2485 0.4973 0 -0.0017 -0.0033 -0.0024 -0.0006 0 0 0.0000 0.0000 0.0000 0 0 0 0.0000 0.0000 0 0 0 0 0.0000 p = 3 2 1 4 5

To verify that, at least numerically, the determinant of A is not negative the parity is used:

d = prod(diag(U)) par = parity(p)

d = -1.0000e-10 par = -1

## Positive definiteness of a matrix

The following verifies rigorously that the matrix A of the previous example is indeed positive definite:

isspd(A)

ans = logical 1

## Prove of regularity of matrix with arbitrary condition number

Based on number theory and mod(p) calculations it is possible to verify the regularity of arbitrarily ill-conditioned matrices in short computing time. If, by chance, the test should not work, a second parameter specifying a kind of probability is possible.

An example with extreme condition number of about 1e100 is as follows.

n = 200; cnd = 1e100; A = randmat(n,cnd); tic reg = isregular(A) T = toc

reg = 1 T = 0.1649

Note that the proof of regularity is rigorous. As is well-known, a proof of singularity in the same manner is an ill-posed problem.

## Full rank of a matrix

Similarly it can be verified that a rectangular matrix has full rank. An answer 1 proves that the matrix has full rank. If the following

A = randn(5,3)*rand(3,5)

A = 1.2583 0.8856 0.2296 1.1701 0.5128 -3.0938 -1.3562 -0.2933 -2.4345 -1.8984 0.1640 -0.6927 0.1683 -0.5619 0.6026 0.0197 -0.0578 -0.4114 0.2490 0.1520 1.5667 -0.2286 0.1080 0.5594 1.6126

is computed in real arithmetic, the result matrix A has at most rank 3. However, due to rounding errors the result has full rank:

isfullrank(A)

ans = 1

For integer input no rounding error occur and the resulting matrix is indeed rank-deficient:

A = randint(-10,5,3)*randint(-10,3,5) isfullrank(A)

A = 64 -31 -76 30 10 -24 53 62 -32 -14 28 21 -28 -2 14 55 23 -16 0 -11 -13 18 15 -11 3 ans = 0

Here randint(-10,5,3) produces a 5x3 matrix with random integers in [-10,10]. The rank-deficiency of A in some Galois field can be verified:

gfpinit(67108859) [L,U,p,q,growth,r] = luwtp(gfp(A))

Galois field prime set to 67108859 ans = 67108859 gfp L = 1 0 0 0 0 0 1 0 0 0 33554435 54907255 1 0 0 16 6100827 62984854 1 0 67108844 67108839 41240078 0 1 gfp U = 67108857 14 28 67108831 21 0 67108848 55 67108843 23 0 0 67108322 6101082 45755788 0 0 0 0 0 0 0 0 0 0 p = 3 4 5 2 1 q = 4 5 1 3 2 growth = 1 r = 3

Mathematically it is shown that A has rank 3 over GF[p], it might have a higher rank over R.

## Random orthogonal matrices

One way to generate matrices with specified singular values uses random orthogonal matrices.

n = 100; Q = randorth(n); norm(Q'*Q-eye(n),inf)

ans = 1.6708e-14

## Relative errors

Displaying relative errors has well-known drawbacks. First, there is the problem of numbers near zero: For x=1e-20, is x+1e-25 nearby? Second, the relative error depends on how close a number is to the next power of 2, and whether it is slightly smaller or larger.

The predecessor and successor of 1 and the relative error to 1 is computed by

succ1 = succ(1) pred1 = pred(1) relsucc = relerr(1,succ1) relpred = relerr(1,pred1)

succ1 = 1.0000 pred1 = 1.0000 relsucc = 2.2204e-16 relpred = 1.1102e-16

## Distance in bits and the accuracy of inclusion functions

For very narrow intervals it might sometimes be simpler to use the distance in bits:

distsucc = distbits(1,succ1) distpred = distbits(1,pred1)

distsucc = 1 distpred = -1

For performance reasons the maximum distance in bits is restricted to 100. The function distbits can be used to track the accuracy of computed inclusions:

x = linspace(0,200,1e6); Y = gamma(intval(x)); close, plot(x,distbits(Y))

## The bit pattern

Also the bit pattern is sometimes useful:

getbits([pred1 1 succ1])

ans = 3×63 char array ' +1.1111111111111111111111111111111111111111111111111111 * 2^-1' ' +1.0000000000000000000000000000000000000000000000000000 * 2^0 ' ' +1.0000000000000000000000000000000000000000000000000001 * 2^0 '

The function applies to intervals as well. It can be used to verify that the bounds are indeed as close as possible:

```
Pi = intval('pi')
getbits(Pi)
distbits(Pi)
```

intval Pi = 3.1415 ans = 4×62 char array 'infimum ' ' +1.1001001000011111101101010100010001000010110100011000 * 2^1' 'supremum ' ' +1.1001001000011111101101010100010001000010110100011001 * 2^1' ans = 1

Displaying bits applies to single precision numbers or the INTLAB fl-toolbox as well:

getbits(single(pi)) flinit(7,50); getbits(fl(5)) getbits(fl(infsup(11,13)))

ans = ' +1.10010010000111111011011 * 2^1' ans = ' +1.010000 * 2^2' ans = 4×16 char array 'infimum ' ' +1.011000 * 2^3' 'supremum ' ' +1.101000 * 2^3'

## Minimum precision to store data

Sometimes it is convenient to know the minimum necessary precision to store a floating-number. For example

x = single(3.375) getbits(x) p = minprec(x)

x = single 3.3750 ans = ' +1.10110000000000000000000 * 2^1' p = 5

The best approximation for x=0.3 requires 53 bits:

x = 0.3 getbits(x) minprec(x)

x = 0.3000 ans = ' +1.0011001100110011001100110011001100110011001100110011 * 2^-2' ans = 53

Here p=5 shows that x requires at least 5 bits to be stored without rounding error. That includes the implicit 1.

## Intelligent help

Sometimes it is not clear how to spell a name and/or we know only parts of a function name, and even that is uncertain. In that case sometimes helpp, an intelligent help function, may be useful. After initialization by just calling "helpp" the following may be obtained:

```
helpp tschebischeff
```

ans = 4×85 char array 'level 6 C:\Intlab Versions\intlab\intval\checkrounding.m ' 'level 6 C:\Intlab Versions\intlab\utility\Fletcher.m ' 'level 9 C:\Intlab Versions\intlab\polynom\chebyshev.m ' 'level 9 C:\Intlab Versions\intlab\polynom\chebyshev2.m '

The answer comes with certain levels representing some kind of distance. Sometimes even gross misspelling still retrieves the correct name. Needless to say that it works the better the longer the string to look for.

## Greatest common divisor and least common multiple

The smallest common denominator of rational number is the least common multiple of the denominators. For example

lcmvec([12 15 36])

ans = 180

Correspondingly, greatest common divisor is computed:

v = [12 15 36] g = gcdvec(v)

v = 12 15 36 g = 3

## All factors of an integer and perfect numbers

The list of all factors of an integer is computed by

factors(123456789)

ans = Columns 1 through 6 1 3 9 3607 3803 10821 Columns 7 through 11 11409 32463 34227 13717421 41152263

It can be used to verify that a number is perfect:

n = 28 f = factors(n) sum(f)

n = 28 f = 1 2 4 7 14 ans = 28

## Symbolic computations with pretty output

The pretty output of symbolic expressions is often useful:

syms x P = 1; for i=1:5 P = P * (x^2-i); end pretty(collect(P,x))

10 8 6 4 2 x - 15 x + 85 x - 225 x + 274 x - 120

However, the statement leading to the pretty output is not visible. That can be done as follows:

```
pretty('collect(P,x)')
```

collect(P,x) = 10 8 6 4 2 x - 15 x + 85 x - 225 x + 274 x - 120

## Checking exponential many possibilities

Given an interval matrix A it is known to be an NP-hard problem to verify non-singularity. The best criterion known still needs 2^n steps for a given n x n matrix. The result is as follows:

Theorem (Rohn). A given interval matrix mA +/- r*E for E denoting the matrix of all 1's, is non-singular if and only if r < 1/sum(abs(inv(A)*x)) for all vectors x with components +/-1 .

The maximum quantity r is the so-called radius of non-singularity. Compared to checking all vertices the theorem reduces the number of determinants from 2^(n^2) to 2^n.

The function bin2vec is useful to generate all signature matrices. The call v = bin2vec(i,n) generates the vector of the n least significant bits of i.

n = 4; A = hadamard(n) Ainv = inv(A); R = 0; for i=1:2^n x = (2*bin2vec(i,n)-1)'; r = sum(abs( Ainv * x )); if r>R, R=r; end end singRad = 1/R

A = 1 1 1 1 1 -1 1 -1 1 1 -1 -1 1 -1 -1 1 singRad = 0.5000

## Traversing all vertices of an interval matrix

The function bin2vec can also be used to visit all vertices of an interval matrix. These are, however, 2^(n^2) vertices. With the method before we checked that the radius of non-singularity of the 8x8 Hadamard matrix is 0.5. Thus A=midrad(hadamard,0.499) must be non-singular. This is checked as follows.

tic n = 4; A = hadamard(n); s = sign(det(A)); for i=1:2^(n^2) D = det( A + 0.499*reshape( 2*bin2vec(i,n^2)-1 , n,n ) ); if D*s<0 disp('The interval matrix is singular') break end end if D*s>0 disp('The interval matrix is non-singular') end T = toc

The interval matrix is non-singular T = 0.5413

The function v = base2vec(i,b,n) does the same for base b.

## 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