SimplexDownhill.h
Go to the documentation of this file.
1 //===========================================================================
2 /*!
3  *
4  *
5  * \brief Nelder-Mead Simplex Downhill Method
6  *
7  *
8  *
9  * \author T. Glasmachers
10  * \date 2015
11  *
12  *
13  * \par Copyright 1995-2015 Shark Development Team
14  *
15  * <BR><HR>
16  * This file is part of Shark.
17  * <http://image.diku.dk/shark/>
18  *
19  * Shark is free software: you can redistribute it and/or modify
20  * it under the terms of the GNU Lesser General Public License as published
21  * by the Free Software Foundation, either version 3 of the License, or
22  * (at your option) any later version.
23  *
24  * Shark is distributed in the hope that it will be useful,
25  * but WITHOUT ANY WARRANTY; without even the implied warranty of
26  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
27  * GNU Lesser General Public License for more details.
28  *
29  * You should have received a copy of the GNU Lesser General Public License
30  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
31  *
32  */
33 //===========================================================================
34 
35 #ifndef SHARK_ALGORITHMS_SIMPLEXDOWNHILL_H
36 #define SHARK_ALGORITHMS_SIMPLEXDOWNHILL_H
37 
38 
40 #include <boost/serialization/vector.hpp>
41 #include <vector>
42 
43 
44 namespace shark {
45 
46 ///
47 /// \brief Simplex Downhill Method
48 ///
49 /// \par
50 /// The Nelder-Mead Simplex Downhill Method is a deterministic direct
51 /// search method. It is known to perform quite well in low dimensions,
52 /// at least for local search.
53 ///
54 /// \par
55 /// The implementation of the algorithm is along the lines of the
56 /// Wikipedia article
57 /// https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method
58 ///
60 {
61 public:
62  /// \brief Default Constructor.
64  { }
65 
66  /// \brief From INameable: return the class name.
67  std::string name() const
68  { return "SimplexDownhill"; }
69 
70  //from ISerializable
71  virtual void read( InArchive & archive )
72  {
73  archive >> m_simplex;
74  archive >> m_best.point;
75  archive >> m_best.value;
76  }
77 
78  virtual void write( OutArchive & archive ) const
79  {
80  archive << m_simplex;
81  archive << m_best.point;
82  archive << m_best.value;
83  }
84 
85  /// \brief Initialization of the optimizer.
86  ///
87  /// The initial simplex is created is a distance of about one around the proposed starting point.
88  ///
89  virtual void init(ObjectiveFunctionType& objectiveFunction, SearchPointType const& startingPoint)
90  {
91  objectiveFunction.init();
92  checkFeatures(objectiveFunction);
93  size_t dim = startingPoint.size();
94 
95  // create the initial simplex
96  m_best.value = 1e100;
97  m_simplex = std::vector<SolutionType>(dim + 1);
98  for (size_t j=0; j<=dim; j++)
99  {
100  RealVector p(dim);
101  for (size_t i=0; i<dim; i++) p(i) = startingPoint(i) + ((i == j) ? 1.0 : -0.5);
102  m_simplex[j].point = p;
103  m_simplex[j].value = objectiveFunction.eval(p);
104  if (m_simplex[j].value < m_best.value) m_best = m_simplex[j];
105  }
106  }
108 
109  /// \brief Step of the simplex algorithm.
110  void step(ObjectiveFunctionType const& objectiveFunction)
111  {
112  size_t dim = m_simplex.size() - 1;
113 
114  // step of the simplex algorithm
115  sort(m_simplex.begin(), m_simplex.end());
116  SolutionType& best = m_simplex[0];
117  SolutionType& worst = m_simplex[dim];
118 
119  // compute centroid
120  RealVector x0(dim, 0.0);
121  for (size_t j=0; j<dim; j++) x0 += m_simplex[j].point;
122  x0 /= (double)dim;
123 
124  // reflection
125  SolutionType xr;
126  xr.point = 2.0 * x0 - worst.point;
127  xr.value = objectiveFunction(xr.point);
128  if (xr.value < m_best.value) m_best = xr; // keep track of best point
129  if (best.value <= xr.value && xr.value < m_simplex[dim-1].value)
130  {
131  // replace worst point with reflected point
132  worst = xr;
133  }
134  else if (xr.value < best.value)
135  {
136  // expansion
137  SolutionType xe;
138  xe.point = 3.0 * x0 - 2.0 * worst.point;
139  xe.value = objectiveFunction(xe.point);
140  if (xe.value < m_best.value) m_best = xe; // keep track of best point
141  if (xe.value < xr.value)
142  {
143  // replace worst point with expanded point
144  worst = xe;
145  }
146  else
147  {
148  // replace worst point with reflected point
149  worst = xr;
150  }
151  }
152  else
153  {
154  // contraction
155  SolutionType xc;
156  xc.point = 0.5 * x0 + 0.5 * worst.point;
157  xc.value = objectiveFunction(xc.point);
158  if (xc.value < m_best.value) m_best = xc; // keep track of best point
159  if (xc.value < worst.value)
160  {
161  // replace worst point with contracted point
162  worst = xc;
163  }
164  else
165  {
166  // reduction
167  for (size_t j=1; j<=dim; j++)
168  {
169  m_simplex[j].point = 0.5 * best.point + 0.5 * m_simplex[j].point;
170  m_simplex[j].value = objectiveFunction(m_simplex[j].point);
171  if (m_simplex[j].value < m_best.value) m_best = m_simplex[j]; // keep track of best point
172  }
173  }
174  }
175  }
176 
177  /// \brief Read access to the current simplex.
178  std::vector<SolutionType> const& simplex()
179  { return m_simplex; }
180 
181 protected:
182  std::vector<SolutionType> m_simplex; ///< \brief Current simplex (algorithm state).
183 };
184 
185 
186 }
187 #endif