// <copyright file="QuadraticGradientProjectionSearch.cs" company="Math.NET">
|
// Math.NET Numerics, part of the Math.NET Project
|
// http://numerics.mathdotnet.com
|
// http://github.com/mathnet/mathnet-numerics
|
//
|
// Copyright (c) 2009-2017 Math.NET
|
//
|
// Permission is hereby granted, free of charge, to any person
|
// obtaining a copy of this software and associated documentation
|
// files (the "Software"), to deal in the Software without
|
// restriction, including without limitation the rights to use,
|
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
// copies of the Software, and to permit persons to whom the
|
// Software is furnished to do so, subject to the following
|
// conditions:
|
//
|
// The above copyright notice and this permission notice shall be
|
// included in all copies or substantial portions of the Software.
|
//
|
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
// OTHER DEALINGS IN THE SOFTWARE.
|
// </copyright>
|
|
using System;
|
using System.Collections.Generic;
|
using IStation.Numerics.LinearAlgebra;
|
|
namespace IStation.Numerics.Optimization
|
{
|
public static class QuadraticGradientProjectionSearch
|
{
|
public static GradientProjectionResult Search(Vector<double> x0, Vector<double> gradient, Matrix<double> hessian, Vector<double> lowerBound, Vector<double> upperBound)
|
{
|
List<bool> isFixed = new List<bool>(x0.Count);
|
List<double> breakpoint = new List<double>(x0.Count);
|
for (int ii = 0; ii < x0.Count; ++ii)
|
{
|
breakpoint.Add(0.0);
|
isFixed.Add(false);
|
if (gradient[ii] < 0)
|
breakpoint[ii] = (x0[ii] - upperBound[ii]) / gradient[ii];
|
else if (gradient[ii] > 0)
|
breakpoint[ii] = (x0[ii] - lowerBound[ii]) / gradient[ii];
|
else
|
{
|
if (Math.Abs(x0[ii] - upperBound[ii]) < 100 * Double.Epsilon || Math.Abs(x0[ii] - lowerBound[ii]) < 100 * Double.Epsilon)
|
breakpoint[ii] = 0.0;
|
else
|
breakpoint[ii] = Double.PositiveInfinity;
|
}
|
}
|
|
var orderedBreakpoint = new List<double>(x0.Count);
|
orderedBreakpoint.AddRange(breakpoint);
|
orderedBreakpoint.Sort();
|
|
// Compute initial state variables
|
var d = -gradient;
|
for (int ii = 0; ii < d.Count; ++ii)
|
if (breakpoint[ii] <= 0.0)
|
d[ii] *= 0.0;
|
|
|
int jj = -1;
|
var x = x0;
|
var f1 = gradient * d;
|
var f2 = 0.5 * d * hessian * d;
|
var sMin = -f1 / f2;
|
var maxS = orderedBreakpoint[0];
|
|
if (sMin < maxS)
|
return new GradientProjectionResult(x + sMin * d, 0,isFixed);
|
|
// while minimum of the last quadratic piece observed is beyond the interval searched
|
while (true)
|
{
|
// update data to the beginning of the interval we're searching
|
jj += 1;
|
x = x + d * maxS;
|
maxS = orderedBreakpoint[jj+1] - orderedBreakpoint[jj];
|
|
int fixedCount = 0;
|
for (int ii = 0; ii < d.Count; ++ii)
|
if (orderedBreakpoint[jj] >= breakpoint[ii])
|
{
|
d[ii] *= 0.0;
|
isFixed[ii] = true;
|
fixedCount += 1;
|
}
|
|
if (Double.IsPositiveInfinity(orderedBreakpoint[jj + 1]))
|
return new GradientProjectionResult(x, fixedCount, isFixed);
|
|
f1 = gradient * d + (x - x0) * hessian * d;
|
f2 = d * hessian * d;
|
|
sMin = -f1 / f2;
|
|
if (sMin < maxS)
|
return new GradientProjectionResult(x + sMin * d, fixedCount, isFixed);
|
else if (jj + 1 >= orderedBreakpoint.Count - 1)
|
{
|
isFixed[isFixed.Count - 1] = true;
|
return new GradientProjectionResult(x + maxS * d, lowerBound.Count, isFixed);
|
}
|
}
|
}
|
|
public readonly struct GradientProjectionResult
|
{
|
public GradientProjectionResult(Vector<double> cauchyPoint, int fixedCount, List<bool> isFixed)
|
{
|
CauchyPoint = cauchyPoint;
|
FixedCount = fixedCount;
|
IsFixed = isFixed;
|
}
|
public Vector<double> CauchyPoint { get; }
|
public int FixedCount { get; }
|
public List<bool> IsFixed { get; }
|
}
|
}
|
}
|
