octave_packages/m/sparse/pcg.m

   1 ## Copyright (C) 2004-2012 Piotr Krzyzanowski
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn  {Function File} {@var{x} =} pcg (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{m1}, @var{m2}, @var{x0}, @dots{})
  21 ## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}, @var{eigest}] =} pcg (@dots{})
  22 ##
  23 ## Solve the linear system of equations @code{@var{A} * @var{x} = @var{b}}
  24 ## by means of the Preconditioned Conjugate Gradient iterative
  25 ## method.  The input arguments are
  26 ##
  27 ## @itemize
  28 ## @item
  29 ## @var{A} can be either a square (preferably sparse) matrix or a
  30 ## function handle, inline function or string containing the name
  31 ## of a function which computes @code{@var{A} * @var{x}}.  In principle
  32 ## @var{A} should be symmetric and positive definite; if @code{pcg}
  33 ## finds @var{A} to not be positive definite, you will get a warning
  34 ## message and the @var{flag} output parameter will be set.
  35 ##
  36 ## @item
  37 ## @var{b} is the right hand side vector.
  38 ##
  39 ## @item
  40 ## @var{tol} is the required relative tolerance for the residual error,
  41 ## @code{@var{b} - @var{A} * @var{x}}.  The iteration stops if
  42 ## @code{norm (@var{b} - @var{A} * @var{x}) <=
  43 ##       @var{tol} * norm (@var{b} - @var{A} * @var{x0})}.
  44 ## If @var{tol} is empty or is omitted, the function sets
  45 ## @code{@var{tol} = 1e-6} by default.
  46 ##
  47 ## @item
  48 ## @var{maxit} is the maximum allowable number of iterations; if
  49 ## @code{[]} is supplied for @code{maxit}, or @code{pcg} has less
  50 ## arguments, a default value equal to 20 is used.
  51 ##
  52 ## @item
  53 ## @var{m} = @var{m1} * @var{m2} is the (left) preconditioning matrix, so that
  54 ## the iteration is (theoretically) equivalent to solving by @code{pcg}
  55 ## @code{@var{P} *
  56 ## @var{x} = @var{m} \ @var{b}}, with @code{@var{P} = @var{m} \ @var{A}}.
  57 ## Note that a proper choice of the preconditioner may dramatically
  58 ## improve the overall performance of the method.  Instead of matrices
  59 ## @var{m1} and @var{m2}, the user may pass two functions which return
  60 ## the results of applying the inverse of @var{m1} and @var{m2} to
  61 ## a vector (usually this is the preferred way of using the preconditioner).
  62 ## If @code{[]} is supplied for @var{m1}, or @var{m1} is omitted, no
  63 ## preconditioning is applied.  If @var{m2} is omitted, @var{m} = @var{m1}
  64 ## will be used as preconditioner.
  65 ##
  66 ## @item
  67 ## @var{x0} is the initial guess.  If @var{x0} is empty or omitted, the
  68 ## function sets @var{x0} to a zero vector by default.
  69 ## @end itemize
  70 ##
  71 ## The arguments which follow @var{x0} are treated as parameters, and
  72 ## passed in a proper way to any of the functions (@var{A} or @var{m})
  73 ## which are passed to @code{pcg}.  See the examples below for further
  74 ## details.  The output arguments are
  75 ##
  76 ## @itemize
  77 ## @item
  78 ## @var{x} is the computed approximation to the solution of
  79 ## @code{@var{A} * @var{x} = @var{b}}.
  80 ##
  81 ## @item
  82 ## @var{flag} reports on the convergence.  @code{@var{flag} = 0} means
  83 ## the solution converged and the tolerance criterion given by @var{tol}
  84 ## is satisfied.  @code{@var{flag} = 1} means that the @var{maxit} limit
  85 ## for the iteration count was reached.  @code{@var{flag} = 3} reports that
  86 ## the (preconditioned) matrix was found not positive definite.
  87 ##
  88 ## @item
  89 ## @var{relres} is the ratio of the final residual to its initial value,
  90 ## measured in the Euclidean norm.
  91 ##
  92 ## @item
  93 ## @var{iter} is the actual number of iterations performed.
  94 ##
  95 ## @item
  96 ## @var{resvec} describes the convergence history of the method.
  97 ## @code{@var{resvec} (i,1)} is the Euclidean norm of the residual, and
  98 ## @code{@var{resvec} (i,2)} is the preconditioned residual norm,
  99 ## after the (@var{i}-1)-th iteration, @code{@var{i} =
 100 ## 1, 2, @dots{}, @var{iter}+1}.  The preconditioned residual norm
 101 ## is defined as
 102 ## @code{norm (@var{r}) ^ 2 = @var{r}' * (@var{m} \ @var{r})} where
 103 ## @code{@var{r} = @var{b} - @var{A} * @var{x}}, see also the
 104 ## description of @var{m}.  If @var{eigest} is not required, only
 105 ## @code{@var{resvec} (:,1)} is returned.
 106 ##
 107 ## @item
 108 ## @var{eigest} returns the estimate for the smallest @code{@var{eigest}
 109 ## (1)} and largest @code{@var{eigest} (2)} eigenvalues of the
 110 ## preconditioned matrix @code{@var{P} = @var{m} \ @var{A}}.  In
 111 ## particular, if no preconditioning is used, the estimates for the
 112 ## extreme eigenvalues of @var{A} are returned.  @code{@var{eigest} (1)}
 113 ## is an overestimate and @code{@var{eigest} (2)} is an underestimate,
 114 ## so that @code{@var{eigest} (2) / @var{eigest} (1)} is a lower bound
 115 ## for @code{cond (@var{P}, 2)}, which nevertheless in the limit should
 116 ## theoretically be equal to the actual value of the condition number.
 117 ## The method which computes @var{eigest} works only for symmetric positive
 118 ## definite @var{A} and @var{m}, and the user is responsible for
 119 ## verifying this assumption.
 120 ## @end itemize
 121 ##
 122 ## Let us consider a trivial problem with a diagonal matrix (we exploit the
 123 ## sparsity of A)
 124 ##
 125 ## @example
 126 ## @group
 127 ## n = 10;
 128 ## A = diag (sparse (1:n));
 129 ## b = rand (n, 1);
 130 ## [l, u, p, q] = luinc (A, 1.e-3);
 131 ## @end group
 132 ## @end example
 133 ##
 134 ## @sc{Example 1:} Simplest use of @code{pcg}
 135 ##
 136 ## @example
 137 ## x = pcg (A,b)
 138 ## @end example
 139 ##
 140 ## @sc{Example 2:} @code{pcg} with a function which computes
 141 ## @code{@var{A} * @var{x}}
 142 ##
 143 ## @example
 144 ## @group
 145 ## function y = apply_a (x)
 146 ##   y = [1:N]' .* x;
 147 ## endfunction
 148 ##
 149 ## x = pcg ("apply_a", b)
 150 ## @end group
 151 ## @end example
 152 ##
 153 ## @sc{Example 3:} @code{pcg} with a preconditioner: @var{l} * @var{u}
 154 ##
 155 ## @example
 156 ## x = pcg (A, b, 1.e-6, 500, l*u)
 157 ## @end example
 158 ##
 159 ## @sc{Example 4:} @code{pcg} with a preconditioner: @var{l} * @var{u}.
 160 ## Faster than @sc{Example 3} since lower and upper triangular matrices
 161 ## are easier to invert
 162 ##
 163 ## @example
 164 ## x = pcg (A, b, 1.e-6, 500, l, u)
 165 ## @end example
 166 ##
 167 ## @sc{Example 5:} Preconditioned iteration, with full diagnostics.  The
 168 ## preconditioner (quite strange, because even the original matrix
 169 ## @var{A} is trivial) is defined as a function
 170 ##
 171 ## @example
 172 ## @group
 173 ## function y = apply_m (x)
 174 ##   k = floor (length (x) - 2);
 175 ##   y = x;
 176 ##   y(1:k) = x(1:k) ./ [1:k]';
 177 ## endfunction
 178 ##
 179 ## [x, flag, relres, iter, resvec, eigest] = ...
 180 ##                    pcg (A, b, [], [], "apply_m");
 181 ## semilogy (1:iter+1, resvec);
 182 ## @end group
 183 ## @end example
 184 ##
 185 ## @sc{Example 6:} Finally, a preconditioner which depends on a
 186 ## parameter @var{k}.
 187 ##
 188 ## @example
 189 ## @group
 190 ## function y = apply_M (x, varargin)
 191 ##   K = varargin@{1@};
 192 ##   y = x;
 193 ##   y(1:K) = x(1:K) ./ [1:K]';
 194 ## endfunction
 195 ##
 196 ## [x, flag, relres, iter, resvec, eigest] = ...
 197 ##      pcg (A, b, [], [], "apply_m", [], [], 3)
 198 ## @end group
 199 ## @end example
 200 ##
 201 ## References:
 202 ##
 203 ## @enumerate
 204 ## @item
 205 ## C.T. Kelley, @cite{Iterative Methods for Linear and Nonlinear Equations},
 206 ## SIAM, 1995. (the base PCG algorithm)
 207 ##
 208 ## @item
 209 ## Y. Saad, @cite{Iterative Methods for Sparse Linear Systems}, PWS 1996.
 210 ## (condition number estimate from PCG) Revised version of this book is
 211 ## available online at @url{http://www-users.cs.umn.edu/~saad/books.html}
 212 ## @end enumerate
 213 ##
 214 ## @seealso{sparse, pcr}
 215 ## @end deftypefn
 216
 217 ## Author: Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
 218 ## Modified by: Vittoria Rezzonico <vittoria.rezzonico@epfl.ch>
 219 ##  - Add the ability to provide the pre-conditioner as two separate matrices
 220
 221 function [x, flag, relres, iter, resvec, eigest] = pcg (A, b, tol, maxit, m1, m2, x0, varargin)
 222
 223   ## M = M1*M2
 224
 225   if (nargin < 7 || isempty (x0))
 226     x = zeros (size (b));
 227   else
 228     x = x0;
 229   endif
 230
 231   if (nargin < 5 || isempty (m1))
 232      exist_m1 = 0;
 233   else
 234      exist_m1 = 1;
 235   endif
 236
 237   if (nargin < 6 || isempty (m2))
 238      exist_m2 = 0;
 239   else
 240      exist_m2 = 1;
 241   endif
 242
 243   if (nargin < 4 || isempty (maxit))
 244     maxit = min (size (b, 1), 20);
 245   endif
 246
 247   maxit += 2;
 248
 249   if (nargin < 3 || isempty (tol))
 250     tol = 1e-6;
 251   endif
 252
 253   preconditioned_residual_out = false;
 254   if (nargout > 5)
 255     T = zeros (maxit, maxit);
 256     preconditioned_residual_out = true;
 257   endif
 258
 259   ## Assume A is positive definite.
 260   matrix_positive_definite = true;
 261
 262   p = zeros (size (b));
 263   oldtau = 1;
 264   if (isnumeric (A))
 265     ## A is a matrix.
 266     r = b - A*x;
 267   else
 268     ## A should be a function.
 269     r = b - feval (A, x, varargin{:});
 270   endif
 271
 272   resvec(1,1) = norm (r);
 273   alpha = 1;
 274   iter = 2;
 275
 276   while (resvec (iter-1,1) > tol * resvec (1,1) && iter < maxit)
 277     if (exist_m1)
 278       if(isnumeric (m1))
 279         y = m1 \ r;
 280       else
 281         y = feval (m1, r, varargin{:});
 282       endif
 283     else
 284       y = r;
 285     endif
 286     if (exist_m2)
 287       if (isnumeric (m2))
 288         z = m2 \ y;
 289       else
 290         z = feval (m2, y, varargin{:});
 291       endif
 292     else
 293       z = y;
 294     endif
 295     tau = z' * r;
 296     resvec (iter-1,2) = sqrt (tau);
 297     beta = tau / oldtau;
 298     oldtau = tau;
 299     p = z + beta * p;
 300     if (isnumeric (A))
 301       ## A is a matrix.
 302       w = A * p;
 303     else
 304       ## A should be a function.
 305       w = feval (A, p, varargin{:});
 306     endif
 307     ## Needed only for eigest.
 308     oldalpha = alpha;
 309     alpha = tau / (p'*w);
 310     if (alpha <= 0.0)
 311       ## Negative matrix.
 312       matrix_positive_definite = false;
 313     endif
 314     x += alpha * p;
 315     r -= alpha * w;
 316     if (nargout > 5 && iter > 2)
 317       T(iter-1:iter, iter-1:iter) = T(iter-1:iter, iter-1:iter) + ...
 318           [1 sqrt(beta); sqrt(beta) beta]./oldalpha;
 319       ## EVS = eig(T(2:iter-1,2:iter-1));
 320       ## fprintf(stderr,"PCG condest: %g (iteration: %d)\n", max(EVS)/min(EVS),iter);
 321     endif
 322     resvec (iter,1) = norm (r);
 323     iter++;
 324   endwhile
 325
 326   if (nargout > 5)
 327     if (matrix_positive_definite)
 328       if (iter > 3)
 329         T = T(2:iter-2,2:iter-2);
 330         l = eig (T);
 331         eigest = [min(l), max(l)];
 332         ## fprintf (stderr, "pcg condest: %g\n", eigest(2)/eigest(1));
 333       else
 334         eigest = [NaN, NaN];
 335         warning ("pcg: eigenvalue estimate failed: iteration converged too fast");
 336       endif
 337     else
 338       eigest = [NaN, NaN];
 339     endif
 340
 341     ## Apply the preconditioner once more and finish with the precond
 342     ## residual.
 343     if (exist_m1)
 344       if (isnumeric (m1))
 345         y = m1 \ r;
 346       else
 347         y = feval (m1, r, varargin{:});
 348       endif
 349     else
 350       y = r;
 351     endif
 352     if (exist_m2)
 353       if (isnumeric (m2))
 354         z = m2 \ y;
 355       else
 356         z = feval (m2, y, varargin{:});
 357       endif
 358     else
 359       z = y;
 360     endif
 361
 362     resvec (iter-1,2) = sqrt (r' * z);
 363   else
 364     resvec = resvec(:,1);
 365   endif
 366
 367   flag = 0;
 368   relres = resvec (iter-1,1) ./ resvec(1,1);
 369   iter -= 2;
 370   if (iter >= maxit - 2)
 371     flag = 1;
 372     if (nargout < 2)
 373       warning ("pcg: maximum number of iterations (%d) reached\n", iter);
 374       warning ("the initial residual norm was reduced %g times.\n", ...
 375                1.0 / relres);
 376     endif
 377   elseif (nargout < 2)
 378     fprintf (stderr, "pcg: converged in %d iterations. ", iter);
 379     fprintf (stderr, "the initial residual norm was reduced %g times.\n",...
 380              1.0/relres);
 381   endif
 382
 383   if (! matrix_positive_definite)
 384     flag = 3;
 385     if (nargout < 2)
 386       warning ("pcg: matrix not positive definite?\n");
 387     endif
 388   endif
 389 endfunction
 390
 391 %!demo
 392 %!
 393 %!      # Simplest usage of pcg (see also 'help pcg')
 394 %!
 395 %!      N = 10;
 396 %!      A = diag ([1:N]); b = rand (N, 1); y =  A \ b; #y is the true solution
 397 %!      x = pcg (A, b);
 398 %!      printf('The solution relative error is %g\n', norm (x - y) / norm (y));
 399 %!
 400 %!      # You shouldn't be afraid if pcg issues some warning messages in this
 401 %!      # example: watch out in the second example, why it takes N iterations
 402 %!      # of pcg to converge to (a very accurate, by the way) solution
 403 %!demo
 404 %!
 405 %!      # Full output from pcg, except for the eigenvalue estimates
 406 %!      # We use this output to plot the convergence history
 407 %!
 408 %!      N = 10;
 409 %!      A = diag ([1:N]); b = rand (N, 1); X =  A \ b; #X is the true solution
 410 %!      [x, flag, relres, iter, resvec] = pcg (A, b);
 411 %!      printf('The solution relative error is %g\n', norm (x - X) / norm (X));
 412 %!      title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||/||b||)');
 413 %!      semilogy([0:iter], resvec / resvec(1),'o-g');
 414 %!      legend('relative residual');
 415 %!demo
 416 %!
 417 %!      # Full output from pcg, including the eigenvalue estimates
 418 %!      # Hilbert matrix is extremely ill conditioned, so pcg WILL have problems
 419 %!
 420 %!      N = 10;
 421 %!      A = hilb (N); b = rand (N, 1); X = A \ b; #X is the true solution
 422 %!      [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], 200);
 423 %!      printf('The solution relative error is %g\n', norm (x - X) / norm (X));
 424 %!      printf('Condition number estimate is %g\n', eigest(2) / eigest (1));
 425 %!      printf('Actual condition number is   %g\n', cond (A));
 426 %!      title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)');
 427 %!      semilogy([0:iter], resvec,['o-g';'+-r']);
 428 %!      legend('absolute residual','absolute preconditioned residual');
 429 %!demo
 430 %!
 431 %!      # Full output from pcg, including the eigenvalue estimates
 432 %!      # We use the 1-D Laplacian matrix for A, and cond(A) = O(N^2)
 433 %!      # and that's the reasone we need some preconditioner; here we take
 434 %!      # a very simple and not powerful Jacobi preconditioner,
 435 %!      # which is the diagonal of A
 436 %!
 437 %!      N = 100;
 438 %!      A = zeros (N, N);
 439 %!      for i=1 : N - 1 # form 1-D Laplacian matrix
 440 %!              A (i:i+1, i:i+1) = [2 -1; -1 2];
 441 %!      endfor
 442 %!      b = rand (N, 1); X = A \ b; #X is the true solution
 443 %!      maxit = 80;
 444 %!      printf('System condition number is %g\n', cond (A));
 445 %!      # No preconditioner: the convergence is very slow!
 446 %!
 447 %!      [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit);
 448 %!      printf('System condition number estimate is %g\n', eigest(2) / eigest(1));
 449 %!      title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)');
 450 %!      semilogy([0:iter], resvec(:,1), 'o-g');
 451 %!      legend('NO preconditioning: absolute residual');
 452 %!
 453 %!      pause(1);
 454 %!      # Test Jacobi preconditioner: it will not help much!!!
 455 %!
 456 %!      M = diag (diag (A)); # Jacobi preconditioner
 457 %!      [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit, M);
 458 %!      printf('JACOBI preconditioned system condition number estimate is %g\n', eigest(2) / eigest(1));
 459 %!      hold on;
 460 %!      semilogy([0:iter], resvec(:,1), 'o-r');
 461 %!      legend('NO preconditioning: absolute residual', ...
 462 %!             'JACOBI preconditioner: absolute residual');
 463 %!
 464 %!      pause(1);
 465 %!      # Test nonoverlapping block Jacobi preconditioner: it will help much!
 466 %!
 467 %!      M = zeros (N, N); k = 4;
 468 %!      for i = 1 : k : N # form 1-D Laplacian matrix
 469 %!              M (i:i+k-1, i:i+k-1) = A (i:i+k-1, i:i+k-1);
 470 %!      endfor
 471 %!      [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit, M);
 472 %!      printf('BLOCK JACOBI preconditioned system condition number estimate is %g\n', eigest(2) / eigest(1));
 473 %!      semilogy ([0:iter], resvec(:,1),'o-b');
 474 %!      legend('NO preconditioning: absolute residual', ...
 475 %!             'JACOBI preconditioner: absolute residual', ...
 476 %!             'BLOCK JACOBI preconditioner: absolute residual');
 477 %!      hold off;
 478 %!test
 479 %!
 480 %!      #solve small diagonal system
 481 %!
 482 %!      N = 10;
 483 %!      A = diag ([1:N]); b = rand (N, 1); X = A \ b; #X is the true solution
 484 %!      [x, flag] = pcg (A, b, [], N+1);
 485 %!      assert(norm (x - X) / norm (X), 0, 1e-10);
 486 %!      assert(flag, 0);
 487 %!
 488 %!test
 489 %!
 490 %!      #solve small indefinite diagonal system
 491 %!      #despite A is indefinite, the iteration continues and converges
 492 %!      #indefiniteness of A is detected
 493 %!
 494 %!      N = 10;
 495 %!      A = diag([1:N] .* (-ones(1, N) .^ 2)); b = rand (N, 1); X = A \ b; #X is the true solution
 496 %!      [x, flag] = pcg (A, b, [], N+1);
 497 %!      assert(norm (x - X) / norm (X), 0, 1e-10);
 498 %!      assert(flag, 3);
 499 %!
 500 %!test
 501 %!
 502 %!      #solve tridiagonal system, do not converge in default 20 iterations
 503 %!
 504 %!      N = 100;
 505 %!      A = zeros (N, N);
 506 %!      for i = 1 : N - 1 # form 1-D Laplacian matrix
 507 %!              A (i:i+1, i:i+1) = [2 -1; -1 2];
 508 %!      endfor
 509 %!      b = ones (N, 1); X = A \ b; #X is the true solution
 510 %!      [x, flag, relres, iter, resvec, eigest] = pcg (A, b, 1e-12);
 511 %!      assert(flag);
 512 %!      assert(relres > 1.0);
 513 %!      assert(iter, 20); #should perform max allowable default number of iterations
 514 %!
 515 %!test
 516 %!
 517 %!      #solve tridiagonal system with 'prefect' preconditioner
 518 %!      #converges in one iteration, so the eigest does not work
 519 %!      #and issues a warning
 520 %!
 521 %!      N = 100;
 522 %!      A = zeros (N, N);
 523 %!      for i = 1 : N - 1 # form 1-D Laplacian matrix
 524 %!              A (i:i+1, i:i+1) = [2 -1; -1 2];
 525 %!      endfor
 526 %!      b = ones (N, 1); X = A \ b; #X is the true solution
 527 %!      [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], [], A, [], b);
 528 %!      assert(norm (x - X) / norm (X), 0, 1e-6);
 529 %!      assert(flag, 0);
 530 %!      assert(iter, 1); #should converge in one iteration
 531 %!      assert(isnan (eigest), isnan ([NaN, NaN]));
 532 %!