BFGS Broyden Fletcher Goldfarb Shanno multivariable optimization

In this article I develop the BFGS / Broyden Fletcher Goldfarb Shanno class which is a quasi Newton numerical method that can be used for unconstrained optimization problems.

Most of the functionality of this class resides in the previously created quasi Newton base class, Quasi Newton Base, from which the BFGS class derives. The primary thing that needs to be done for this class after inheriting from Quasi Newton Base is to override the Update Inverse Hessian virtual function and use it to compute the approximate inverse Hessian using the BFGS update formula.

The BFGS inverse Hessian update formula is shown below.

This equation is implemented in the overridden virtual function Update Inverse Hessian.

The code for the BFGS class is given below.

public class BFGS : QuasiNewtonBase
    {
        public BFGS(Func<double> function, Parameter[] optimizationParameters, double inverseHessianResetMultiple, int numberOfDerivativePoints) :
            base(function, optimizationParameters, inverseHessianResetMultiple, numberOfDerivativePoints)
        {
        }
 
        public BFGS(Func<double> function, Parameter[] optimizationParameters)
            : this(function, optimizationParameters, 2.0, 3)
        {
        }
 
        protected override void UpdateInverseHessian()
        {
            int numberOfParameters = _solvedForParameters.Length;
 
            double[] positionDelta = new double[numberOfParameters];
            double[] gradientDelta = new double[numberOfParameters];
            for (int i = 0; i < numberOfParameters; i++)
            {
                gradientDelta[i] = _currentGradient[i] - _lastGradient[i];
                positionDelta[i] = _currentPosition[i] - _lastPosition[i];
            }
 
            Matrix YXt = ComputeQNProduct(gradientDelta, positionDelta);
            double YtX = ComputeInnerProduct(gradientDelta, positionDelta);
            Matrix YXtOverYtX = YXt / YtX;
            Matrix IdentMinusYXtOverYtX = Matrix.Identity(numberOfParameters) - YXtOverYtX;
 
            Matrix T1 = IdentMinusYXtOverYtX.Transpose() * _inverseHessian * IdentMinusYXtOverYtX;
 
            Matrix XXt = ComputeQNProduct(positionDelta, positionDelta);
            Matrix T2 = XXt / YtX;
            Matrix matInverseHessianNew = T1 + T2;
            _inverseHessian = matInverseHessianNew;
        }
    }