From 48a7e4e418ee26465a4d3a24e45e26cf7e90eb71 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Fri, 22 Mar 2013 14:47:58 -0400 Subject: [PATCH] Add a second reference for the PCGM and make it more resistant to accumulated roundoff error. --- ...preconditioned_conjugate_gradient_method.m | 23 +++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) diff --git a/optimization/preconditioned_conjugate_gradient_method.m b/optimization/preconditioned_conjugate_gradient_method.m index b7999d9..4e68ccb 100644 --- a/optimization/preconditioned_conjugate_gradient_method.m +++ b/optimization/preconditioned_conjugate_gradient_method.m @@ -58,8 +58,15 @@ function [x, k] = preconditioned_conjugate_gradient_method(Q, ... % REFERENCES: % % 1. Guler, Osman. Foundations of Optimization. New York, Springer, - % 2010. + % 2010. % + % 2. Shewchuk, Jonathan Richard. An Introduction to the Conjugate + % Gradient Method Without the Agonizing Pain, Edition 1.25. + % August 4, 1994. + % + + % We use this in the inner loop. + sqrt_n = floor(sqrt(length(x0))); % Set k=0 first, that way the references to xk,rk,zk,dk which % immediately follow correspond (semantically) to x0,r0,z0,d0. @@ -85,11 +92,19 @@ function [x, k] = preconditioned_conjugate_gradient_method(Q, ... % do them both, so we precompute the more expensive operation. Qdk = Q * dk; - % After substituting the two previously-created variables, the - % following algorithm occurs verbatim in the reference. alpha_k = rkzk/(dk' * Qdk); x_next = xk + (alpha_k * dk); - r_next = rk + (alpha_k * Qdk); + + % 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. + if (mod(k, sqrt_n) == 0) + r_next = Q*x_next - b; + else + r_next = rk + (alpha_k * Qdk); + end + z_next = M \ r_next; beta_next = (r_next' * z_next)/rkzk; d_next = -z_next + beta_next*dk; -- 2.44.2