optimization/simple_preconditioned_cgm.m

   1 function [x, k] = simple_preconditioned_cgm(Q, ...
   2                                             M, ...
   3                                             b, ...
   4                                             x0, ...
   5                                             tolerance, ...
   6                                             max_iterations)
   7   %
   8   % Solve,
   9   %
  10   %   Qx = b
  11   %
  12   % or equivalently,
  13   %
  14   %   min [phi(x) = (1/2)*<Qx,x> + <b,x>]
  15   %
  16   % using the preconditioned conjugate gradient method (14.54 in
  17   % Guler).
  18   %
  19   % INPUT:
  20   %
  21   %   - ``Q`` -- The coefficient matrix of the system to solve. Must
  22   %     be positive definite.
  23   %
  24   %   - ``M`` -- The preconditioning matrix. If the actual matrix used
  25   %     to precondition ``Q`` is called ``C``, i.e. ``C^(-1) * Q *
  26   %     C^(-T) == \bar{Q}``, then M=CC^T. Must be symmetric positive-
  27   %     definite. See for example Golub and Van Loan.
  28   %
  29   %   - ``b`` -- The right-hand-side of the system to solve.
  30   %
  31   %   - ``x0`` -- The starting point for the search.
  32   %
  33   %   - ``tolerance`` -- How close ``Qx`` has to be to ``b`` (in
  34   %     magnitude) before we stop.
  35   %
  36   %   - ``max_iterations`` -- The maximum number of iterations to
  37   %     perform.
  38   %
  39   % OUTPUT:
  40   %
  41   %   - ``x`` - The solution to Qx=b.
  42   %
  43   %   - ``k`` - The ending value of k; that is, the number of
  44   %   iterations that were performed.
  45   %
  46   % NOTES:
  47   %
  48   % All vectors are assumed to be *column* vectors.
  49   %
  50   % REFERENCES:
  51   %
  52   %   1. Guler, Osman. Foundations of Optimization. New York, Springer,
  53   %   2010.
  54   %
  55
  56   % This isn't great in practice, since the CGM is usually used on
  57   % huge sparse systems.
  58   Ct = chol(M);
  59   C = Ct';
  60   C_inv = inv(C);
  61   Ct_inv = inv(Ct);
  62
  63   Q_bar = C_inv * Q * Ct_inv;
  64   b_bar = C_inv * b;
  65
  66   % But it sure is easy.
  67   [x_bar, k] = vanilla_cgm(Q_bar, b_bar, x0, tolerance, max_iterations);
  68
  69   % The solution to Q_bar*x_bar == b_bar is x_bar = Ct*x.
  70   x = Ct_inv * x_bar;
  71 end