conjugate_gradient.py
4b6da740ede9master

conjugate_gradient.py
View Options

			#!/usr/bin/env python3

			import numpy as np


			def solve(func, x0, jac=None, hess=None,
			max_iter=50, tol=1e-8,
			callback=None):
			'Minimize function with home-made conjugate gradient'

			if jac is None:
			raise RuntimeError('this method needs a jacobian to be provided')

			if hess(0)[1, 0] != hess(0)[0, 1]:
			raise RuntimeError('this method requires a symmetric matrix A')

			if hess is None:
			raise RuntimeError('this method needs a hessian to be provided')

			x_old = x0
			r_old = -jac(x0)
			p_old = r_old

			# conjugate iteration loop
			for i in range(0, max_iter):

			# fetch the hessian
			A = hess(x_old)

			Ap_old = A@p_old

			alpha = (r_old@r_old) / (p_old@Ap_old)

			x_new = x_old + alpha * p_old
			r_new = r_old - alpha * Ap_old

			# call the callback over iterations
			if callback:
			callback(x_new)

			# evaluate stop criterion
			if np.linalg.norm(x_new-x_old) < tol:
			return x_new

			beta = (r_new@r_new) / (r_old@r_old)

			p_old = r_new + beta * p_old
			r_old = r_new
			x_old = x_new

			return False