DEMOINTLAB_LARGER Some larger examples with INTLAB

Dense linear systems
Ill-conditioned dense linear systems
Sparse linear systems I
Sparse linear systems II
Larger least squares problems
Sparse least squares problems
An nonlinear optimization problem in 100 unknowns
An nonlinear optimization problem in 10,000 unknowns
Enjoy INTLAB

Following are some larger examples using INTLAB, the Matlab toolbox for Reliable Computing. All computations are on my 3.0 GHz Laptop using Matlab R2020a under Windows.

For more larger sparse linear and nonlinear problems see the demo dsparse.

Dense linear systems

The following generates a dense linear system with n=5000 unknowns randomly with solution approximately all 1's. Since random matrices are generally well-conditioned, this is no real challenge concerning verification of the result.

Here and in the following we display the computing time for the Matlab built-in solver and for our verification routines. Note that this compares apples with oranges: the results of the verification routine are mathematically proved to be correct, including all rounding errors and the proof of non-singularity of the input matrix, whereas approximations are sometimes not correct, even without warning (see e.g. the section "Larger least squares problems").

Following the computing time for the Matlab solver A\b and for the verification INTLAB algorithm verifylss, some components of the solution as well as the minimum and median number of correct digits is displayed.

format short
n = 5000;
A = randn(n);
x = ones(n,1);
b = A*x;

tic
x = A\b;
disp(sprintf('Time for the built-in Matlab solver %5.1f [sec]',toc))

tic
X = verifylss(A,b);
disp(sprintf('Time for the verification algorithm %5.1f [sec]',toc))

format long
disp('First and last three components: approximation and inclusion')
for i=[1:3 n-2:n]
  disp(sprintf('%17.7e %53s',x(i),infsup(X(i))))
  if i==3, disp([blanks(30) '...']), end
end
format short

r = relacc(X);
disp('Minimum and median number of correct digits')
[min(r) median(r)]

Time for the built-in Matlab solver   0.8 [sec]
Time for the verification algorithm  12.9 [sec]
First and last three components: approximation and inclusion
    1.0000000e+00            [   0.99999999999880,   0.99999999999881] 
    1.0000000e+00            [   1.00000000000172,   1.00000000000173] 
    1.0000000e+00            [   1.00000000000158,   1.00000000000160] 
                              ...
    1.0000000e+00            [   0.99999999999882,   0.99999999999883] 
    1.0000000e+00            [   1.00000000000039,   1.00000000000040] 
    1.0000000e+00            [   0.99999999999879,   0.99999999999880] 
Minimum and median number of correct digits
ans =
   1.0e-15 *
    0.2220    0.2220

Since the right hand side b is computed as A*x in floating-point, the true solution is approximately the vector of 1's, but not exactly. To force the solution to include the vector of 1's, the right hand side is computed as an inclusion of A*b. Such methods are often used as tests for verification algorithms.

bb = A*intval(ones(n,1));

tic
X = A\bb;
T = toc

format long
disp('First and last three components: approximation and inclusion')
for i=[1:3 n-2:n]
  disp(sprintf('%17.7e %53s',x(i),infsup(X(i))))
  if i==3, disp([blanks(30) '...']), end
end
format short

r = relacc(X);
disp('Minimum and median number of correct digits')
[min(r) median(r)]

T =
   13.1532
First and last three components: approximation and inclusion
    1.0000000e+00            [   0.99999999595866,   1.00000000404339] 
    1.0000000e+00            [   0.99999999324595,   1.00000000675051] 
    1.0000000e+00            [   0.99999999448354,   1.00000000551352] 
                              ...
    1.0000000e+00            [   0.99999999561357,   1.00000000438854] 
    1.0000000e+00            [   0.99999999819875,   1.00000000180063] 
    1.0000000e+00            [   0.99999999599235,   1.00000000401000] 
Minimum and median number of correct digits
ans =
   1.0e-08 *
    0.0462    0.7116

The computing time is roughly the same, but the inclusion is less accurate. However, the right hand side is now an interval vector, and the solution of all linear systems with a right hand side within bb is included.

For cond(A)~10^k, according to the well-known rule of thumb in numerical analysis, the accuracy of the inclusion should be roughly the number of correct digits in bb minus k. This is indeed true.

accX = median(r)
median(relacc(bb)) - log10(cond(A))

accX =
   7.1156e-09
ans =
   -4.8072

Ill-conditioned dense linear systems

Next an ill-conditioned linear system with n=5000 unknowns is generated with solution again roughly the vector of 1's. The condition number is approximately 10^14.

The computing time for the Matlab solver A\b and for the verification INTLAB algorithm verifylss, some components of the solution as well as the minimum and median number of correct digits is displayed.

The condition number implies that the accuracy of the inclusion should be roughly 16-14 = 2 correct digits. This indeed true.

format short
n = 5000;
A = randmat(n,1e14);
x = ones(n,1);
b = A*x;

tic
x = A\b;
disp(sprintf('Time for the built-in Matlab solver %5.1f [sec]',toc))

tic
X = verifylss(A,b);
disp(sprintf('Time for the verification algorithm %5.1f [sec]',toc))

format long
disp('First and last three components: approximation and inclusion')
for i=[1:3 n-2:n]
  disp(sprintf('%17.7e %53s',x(i),infsup(X(i))))
  if i==3, disp([blanks(30) '...']), end
end
format short

r = relacc(X);
disp('Minimum and median number of correct digits')
[min(r) median(r)]

disp('Median relative error of Matlab''s built-in solver')
median(relerr(x,X))

Time for the built-in Matlab solver   0.9 [sec]
Time for the verification algorithm  21.8 [sec]
First and last three components: approximation and inclusion
    9.9726806e-01            [   1.00032038514138,   1.00032038514140] 
    9.8414449e-01            [   0.99949208517136,   0.99949208517138] 
    1.1167028e+00            [   0.99971686209827,   0.99971686209832] 
                              ...
    1.0034367e+00            [   1.00026905070490,   1.00026905070492] 
    1.0069937e+00            [   0.99992340247084,   0.99992340247088] 
    1.0068903e+00            [   1.00094308064935,   1.00094308064939] 
Minimum and median number of correct digits
ans =
   1.0e-13 *
    0.0666    0.1065
Median relative error of Matlab's built-in solver
ans =
    0.0021

Sparse linear systems I

By the principle of the used method, mainly symmetric positive definite matrices can be treated. The performance for general sparse matrices is not good; alas, basically no better method is known.

Consider for example matrix #356 from the Florida matrix market of dimension 52,329 with 2.6 million nonzero elements. The matrix looks as follows.

Problem = ssget(356);
A = Problem.A;
n = size(A,1)
nnzA = nnz(A)
b = A*ones(n,1);
close
spy(A)

We display the timing the Matlab solver and the verification routine verifylss, and show the minimum and median accuracy of the inclusion. Note that the estimated condition number is 2e14 which is prohibitive to use normal equations.

CndEst = condest(A)

tic
x = A\b;
disp(sprintf('Time for the built-in Matlab solver %5.1f [sec]',toc))

tic
X = verifylss(A,b);
disp(sprintf('Time for the verification algorithm %5.1f [sec]',toc))

format long
disp('First and last three components: approximation and inclusion')
for i=[1:3 n-2:n]
  disp(sprintf('%17.7e %53s',x(i),infsup(X(i))))
  if i==3, disp([blanks(30) '...']), end
end
format short
r = relacc(X);
disp('Minimum and median number of correct digits')
[min(r) median(r)]

CndEst =
   2.2282e+14
Time for the built-in Matlab solver   6.1 [sec]
Time for the verification algorithm   5.9 [sec]
First and last three components: approximation and inclusion
    9.9999998e-01            [   1.00000000105935,   1.00000000105936] 
    9.9999998e-01            [   0.99999999771329,   0.99999999771330] 
    1.0000000e+00            [   0.99999999936704,   0.99999999936705] 
                              ...
    9.9999999e-01            [   1.00000000073318,   1.00000000073319] 
    9.9999998e-01            [   0.99999999832736,   0.99999999832737] 
    1.0000000e+00            [   1.00000000009755,   1.00000000009757] 
Minimum and median number of correct digits
ans =
   1.0e-15 *
    0.2220    0.4441

Note that the verification algorithm and the built-in Matlab solver require about the same computing time, and the result by verifylss is maximally accurate and mathematically verified to be correct.

Sparse linear systems II

Sometimes the verification routine is about as fast or even faster than the built-in Matlab solver. The next test matrix is #938 from the Florida matrix market. This matrix has dimension 36,000 with about 14 million nonzero elements.

Problem = ssget(938);
A = Problem.A;
n = size(A,1)
nnzA = nnz(A)
b = A*ones(n,1);
close
spy(A)

The estimated condition number is about 2.5e7. Again the verification routine is faster than the approximate solver.

CndEst = condest(A)

tic
x = A\b;
disp(sprintf('Time for the built-in Matlab solver %5.1f [sec]',toc))

tic
X = verifylss(A,b);
disp(sprintf('Time for the verification algorithm %5.1f [sec]',toc))

format long
disp('First and last three components: approximation and inclusion')
for i=[1:3 n-2:n]
  disp(sprintf('%17.7e %53s',x(i),infsup(X(i))))
  if i==3, disp([blanks(30) '...']), end
end
format short

r = relacc(X);
disp('Minimum and median number of correct digits')
[min(r) median(r)]

CndEst =
   2.4397e+07
Time for the built-in Matlab solver  55.0 [sec]
Time for the verification algorithm  53.3 [sec]
First and last three components: approximation and inclusion
    1.0000000e+00            [   1.00000000000328,   1.00000000000329] 
    1.0000000e+00            [   0.99999999999969,   0.99999999999970] 
    1.0000000e+00            [   1.00000000000027,   1.00000000000028] 
                              ...
    1.0000000e+00            [   1.00000000000397,   1.00000000000398] 
    1.0000000e+00            [   0.99999999999869,   0.99999999999870] 
    1.0000000e+00            [   0.99999999999967,   0.99999999999968] 
Minimum and median number of correct digits
ans =
   1.0e-15 *
    0.2220    0.4441

The accuracy of the inclusion is as expected. When applying an a priori minimum degree sorting, the computing time for verification does not change too much, but Matlab's built-in solver becomes much faster:

tic
p = symamd(A);
x = A(p,p)\b(p);
disp(sprintf('Time for the built-in Matlab solver with preordering %5.1f [sec]',toc))

tic
X = verifylss(A(p,p),b(p));
disp(sprintf('Time for the verification algorithm %5.1f [sec]',toc))

Time for the built-in Matlab solver with preordering  10.5 [sec]
Time for the verification algorithm  71.1 [sec]

Larger least squares problems

We first generate a dense 5000x500 matrix with condition number 1e12 to solve the corresponding least squares problem. The right hand side is the vector of 1's. The computing time of the built-in Matlab solver and the verification routine is displayed.

format short
m = 5000; n = 500;
A = randmat([m n],1e12);
b = ones(m,1);

tic
x = A\b;
disp(sprintf('Time for the built-in Matlab solver %5.1f [sec]',toc))

tic
X = verifylss(A,b);
disp(sprintf('Time for the verification algorithm %5.1f [sec]',toc))

Time for the built-in Matlab solver   0.1 [sec]
Time for the verification algorithm   2.0 [sec]

Next we show some components of the approximate solution computed x by Matlab and the verified inclusion X by INTLAB. From the accuracy of the verified inclusion, the accuracy of the Matlab approximation can be judged.

format long
disp('First and last three components: approximation and inclusion')
for i=[1:3 n-2:n]
  disp(sprintf('%17.7e %53s',x(i),infsup(X(i))))
  if i==3, disp([blanks(30) '...']), end
end
format short

r = relacc(X);
disp('Minimum and median number of correct digits')
[min(r) median(r)]

First and last three components: approximation and inclusion
   -1.0755408e+11  [ -1.075553258673534e+011, -1.075553258673533e+011] 
    3.5454138e+10  [  3.545468865163642e+010,  3.545468865163647e+010] 
    1.1586072e+11  [  1.158662070143869e+011,  1.158662070143870e+011] 
                              ...
   -8.5685821e+10  [ -8.568860287185251e+010, -8.568860287185246e+010] 
    3.4676452e+10  [  3.468419051714758e+010,  3.468419051714763e+010] 
    6.5325447e+10  [  6.533178840574877e+010,  6.533178840574879e+010] 
Minimum and median number of correct digits
ans =
   1.0e-15 *
    0.1143    0.3602

Sparse least squares problems

Following we display the timing and accuracy of the built-in Matlab routine and the verification routine verifylss for a larger least squares problem, namely matrix #2201. This is a problem with 37932 for 331 unknowns and about 137 thousand nonzero elements. The right hand side is again the vector of 1's.

Problem = ssget(2201);
A = Problem.A;
[m n] = size(A)
b = ones(m,1);

tic
x = A\b;
disp(sprintf('Time for the built-in Matlab solver %5.1f [sec]',toc))

tic
X = verifylss(A,b);
disp(sprintf('Time for the verification algorithm %5.1f [sec]',toc))

v = [1:3 n-2:n];
format long
disp('Inclusion of the first and last three components')
X(v)
format short

r = relacc(X);
disp('Minimum and median number of correct digits')
[min(r) median(r)]

m =
       37932
n =
   331
Time for the built-in Matlab solver   0.2 [sec]
Time for the verification algorithm   1.1 [sec]
Inclusion of the first and last three components
intval ans = 
   0.86939173774958
   0.91891829446876
   0.93916745272555
  -0.99656428275776
  -0.79264335897383
   0.20386414427256
Minimum and median number of correct digits
ans =
   1.0e-15 *
    0.0009    0.2220

In this case we can judge from the inclusion that about 16 digits of the approximation are correct. With that information we ca judge that indeed the Matlab approximate solution is accurate to 13 digits as well.

An nonlinear optimization problem in 100 unknowns

Consider the problem "bdqrtic" in

  Source: Problem 61 in
     A.R. Conn, N.I.M. Gould, M. Lescrenier and Ph.L. Toint,
     "Performance of a multifrontal scheme for partially separable optimization",
      Report 88/4, Dept of Mathematics, FUNDP (Namur, B), 1988.
      Copyright (C) 2001 Princeton University
      All Rights Reserved
  see bottom of file test_h.m

The model problem is

   N = length(x);      % model problem: initial approximation x=ones(N,1);
   I = 1:N-4;
   y = sum( (-4*x(I)+3.0).^2 ) + sum( ( x(I).^2 + 2*x(I+1).^2 + ...
             3*x(I+2).^2 + 4*x(I+3).^2 + 5*x(N).^2 ).^2 );

This function is evaluated by

   index = 2;
   y = test_h(x,index);

We first solve the corresponding nonlinear system in only 100 unknowns to compare with Matlab's built-in fminsearch.

n = 100;
index = 2;

disp('Floating-point approximation by fminsearch with restart')
optimset.Display = 'off';
x = ones(n,1);
tic
for i=1:5
  x = fminsearch(@(x)test_h(x,index),x,optimset);
  y = test_h(x,index);
  disp(sprintf('iteration %1d and current minimal value %7.1f',i,y))
end
disp(sprintf('Time for fminsearch with 5 restarts %5.1f [sec]',toc))
disp(' ')

xs = ones(n,1);
tic
X = verifylocalmin('test_h',xs,0,index);
disp(sprintf('Time for the verification algorithm %5.1f [sec]',toc))
Y = test_h(X,index);
disp(sprintf('Minimal value for stationary point %7.1f',Y.mid))

r = relacc(X);
disp('Minimum and median number of correct digits of stationary point')
[min(r) median(r)]

Floating-point approximation by fminsearch with restart
iteration 1 and current minimal value  6733.3
iteration 2 and current minimal value  2529.4
iteration 3 and current minimal value   587.9
iteration 4 and current minimal value   406.9
iteration 5 and current minimal value   378.9
Time for fminsearch with 5 restarts  25.2 [sec]
 
Time for the verification algorithm   0.5 [sec]
Minimal value for stationary point   378.8
Minimum and median number of correct digits of stationary point
ans =
   1.0e-13 *
    0.1015    0.4064

After 5 restarts, the approximation by fminsearch is of reasonable accuracy. However, that we know only by the verification method. The built-in Matlab routine fminsearch uses the Nelder-Mead algorithm without derivative, thus it is slow even for few unknowns.

An nonlinear optimization problem in 10,000 unknowns

Next we solve the previous nonlinear system in 10,000 unknowns with verification. The given starting vector is again x = ones(n,1). Note that during the computation x will be a vector of Hessians, each carrying a Hessian matrix, in total 10000^3 = 1e12 elements or 8 TeraByte - if not stored sparse.

sparsehessian(0)    % store Hessians always sparse
n = 10000;
index = 2;
tic
X = verifylocalmin('test_h',ones(n,1),0,index);
disp(sprintf('Time for the verification algorithm %5.1f [sec]',toc))

r = relacc(X);
disp('Minimum and median number of correct digits')
[min(r) median(r)]

===> Hessian derivatives always stored sparse
ans =
     0
Time for the verification algorithm  46.6 [sec]
Minimum and median number of correct digits
ans =
   1.0e-11 *
    0.0071    0.4517

Notice the high accuracy of the result. Mathematically, the interval vector X is proved to contain not only a stationary point, but a true (local) minimum.

The proof for positive definiteness is included in the verification routine verifylocalmin. That proof may be performed separately as follows.

tic
y = test_h(hessianinit(X),index);
isLocalMinimum = isspd(y.hx)
disp(sprintf('Time for the verification algorithm %5.1f [sec]',toc))

isLocalMinimum =
  logical
   1
Time for the verification algorithm   0.8 [sec]

The latter command verified that the Hessian at all points in X is s.p.d., among them at the stationary point xx.

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