From: Michael Orlitzky <michael@orlitzky.com>
Date: Sat, 23 Mar 2013 02:30:26 +0000 (-0400)
Subject: Avoid divide-by-zero in preconditioned_conjugate_gradient_method().
X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=4c5f15b69f6795a22c662a9e22ddd649c65872cb;p=octave.git

Avoid divide-by-zero in preconditioned_conjugate_gradient_method().
Use the infinity norm (instead of the 2-norm) in the PCGM.
Simplify the stopping conditions in the PCGM.
---

diff --git a/optimization/preconditioned_conjugate_gradient_method.m b/optimization/preconditioned_conjugate_gradient_method.m
index 35fecaa..cd88608 100644
--- a/optimization/preconditioned_conjugate_gradient_method.m
+++ b/optimization/preconditioned_conjugate_gradient_method.m
@@ -77,14 +77,7 @@ function [x, k] = preconditioned_conjugate_gradient_method(Q, ...
   zk = M \ rk;
   dk = -zk;
 
-  while (k <= max_iterations)
-
-    if (norm(rk) < tolerance)
-      % Check our stopping condition. This should catch the k=0 case.
-      x = xk;
-      return;
-    end
-
+  while (k <= max_iterations && norm(rk, 'inf') > tolerance)
     % Used twice, avoid recomputation.
     rkzk = rk' * zk;
 
@@ -92,13 +85,22 @@ function [x, k] = preconditioned_conjugate_gradient_method(Q, ...
     % do them both, so we precompute the more expensive operation.
     Qdk = Q * dk;
 
-    alpha_k = rkzk/(dk' * Qdk);
+    % We're going to divide by this quantity...
+    dkQdk = dk' * Qdk;
+
+    % So if it's too close to zero, we replace it with something
+    % comparable but non-zero.
+    if (abs(dkQdk) < eps)
+      dkQdk = sign(dkQdk)*eps;
+    end
+
+    alpha_k = rkzk/dkQdk;
     x_next = xk + (alpha_k * dk);
 
     % The recursive definition of r_next is prone to accumulate
     % roundoff error. When sqrt(n) divides k, we recompute the
-    % residual to minimize this error. This modification is due to the
-    % second reference.
+    % residual to minimize this error. This modification was suggested
+    % by the second reference.
     if (mod(k, sqrt_n) == 0)
       r_next = Q*x_next - b;
     else
@@ -120,7 +122,6 @@ function [x, k] = preconditioned_conjugate_gradient_method(Q, ...
     dk = d_next;
   end
 
-  % The algorithm didn't converge, but we still want to return the
-  % terminal value of xk.
+  % If we make it here, one of the two stopping conditions was met.
   x = xk;
 end