From 4c5f15b69f6795a22c662a9e22ddd649c65872cb Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Fri, 22 Mar 2013 22:30:26 -0400 Subject: [PATCH] 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. --- ...preconditioned_conjugate_gradient_method.m | 27 ++++++++++--------- 1 file changed, 14 insertions(+), 13 deletions(-) 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 -- 2.43.2