function x_star = conjugate_gradient_method(A, b, x0, tolerance) % % Solve, % % Ax = b % % or equivalently, % % min [phi(x) = (1/2)* + ] % % using Algorithm 5.2 in Nocedal and Wright. % % INPUT: % % - ``A`` -- The coefficient matrix of the system to solve. Must % be positive definite. % % - ``b`` -- The right-hand-side of the system to solve. % % - ``x0`` -- The starting point for the search. % % - ``tolerance`` -- How close ``Ax`` has to be to ``b`` (in % magnitude) before we stop. % % OUTPUT: % % - ``x_star`` - The solution to Ax=b. % % NOTES: % % All vectors are assumed to be *column* vectors. % zero_vector = zeros(length(x0), 1); k = 0; xk = x0; rk = A*xk - b; % The first residual must be computed the hard way. pk = -rk; while (norm(rk) > tolerance) alpha_k = step_length_cgm(rk, A, pk); x_next = xk + alpha_k*pk; r_next = rk + alpha_k*A*pk; beta_next = (r_next' * r_next)/(rk' * rk); p_next = -r_next + beta_next*pk; k = k + 1; xk = x_next; rk = r_next; pk = p_next; end x_star = xk; end