ʻOhana
Population structure, admixture history, and selection using learning methods.
jade.qpas.hpp
1 /* -------------------------------------------------------------------------
2  Ohana
3  Copyright (c) 2015-2020 Jade Cheng (\___/)
4  Jade Cheng <info@jade-cheng.com> (='.'=)
5  ------------------------------------------------------------------------- */
6 
7 #ifndef JADE_QPAS_HPP__
8 #define JADE_QPAS_HPP__
9 
10 #include "jade.matrix.hpp"
11 
12 namespace jade
13 {
14  ///
15  /// A template for a class that implements the Quadratic Programming for
16  /// Active Set algorithm.
17  ///
18  template <typename TValue>
19  class basic_qpas
20  {
21  public:
22  /// The value type.
23  typedef TValue value_type;
24 
25  /// The matrix type.
26  typedef basic_matrix<value_type> matrix_type;
27 
28  ///
29  /// Loops over the active set and computes a delta vector and a new
30  /// active set.
31  ///
32  /// \param b_vec The B vector.
33  /// \param coefficients_mat The coefficients matrix.
34  /// \param hessian_mat The Hessian matrix.
35  /// \param derivative_vec The derivative vector.
36  /// \param fixed_active_set The fixed active set.
37  /// \param active_set The active set.
38  /// \param delta_vec The delta vector.
39  ///
40  static void loop_over_active_set(
41  const matrix_type & b_vec,
42  const matrix_type & coefficients_mat,
43  const matrix_type & hessian_mat,
44  const matrix_type & derivative_vec,
45  const std::vector<size_t> & fixed_active_set,
46  std::vector<size_t> & active_set,
47  matrix_type & delta_vec)
48  {
49  assert(!fixed_active_set.empty() || !active_set.empty());
50 
51  const auto K = hessian_mat.get_height();
52  assert(K > 0);
53 
54  const auto inequality_constraint_count =
55  b_vec.get_length() - fixed_active_set.size();
56 
57  assert(b_vec.is_vector());
58  assert(coefficients_mat.get_height() == b_vec.get_length());
59  assert(fixed_active_set.size() <= K);
60  assert(delta_vec.is_vector());
61  assert(delta_vec.get_height() == K);
62 
63  std::set<unsigned long> visited_sets;
64 
65  const auto insert_key = [&visited_sets, &active_set]() -> bool
66  {
67  auto key = 0UL;
68  for (const auto i : active_set)
69  key |= 1UL << i;
70  if (visited_sets.find(key) != visited_sets.end())
71  return false;
72  visited_sets.insert(key);
73  return true;
74  };
75 
76  while (insert_key())
77  {
78  std::vector<size_t> merged_active_set;
79  merged_active_set.insert(
80  merged_active_set.end(),
81  active_set.begin(),
82  active_set.end());
83  merged_active_set.insert(
84  merged_active_set.end(),
85  fixed_active_set.begin(),
86  fixed_active_set.end());
87 
88  matrix_type lagrangian_vec (
89  merged_active_set.size(),
90  merged_active_set.empty() ? 0 : 1);
91 
92  matrix_type try_delta_vec (K, 1);
93 
94  _kkt(b_vec,
95  coefficients_mat,
96  hessian_mat,
97  derivative_vec,
98  merged_active_set,
99  try_delta_vec,
100  lagrangian_vec);
101 
102  std::vector<size_t> violated_indices;
103 
104  if (active_set.size() < K - fixed_active_set.size())
105  {
106  for (size_t i = 0; i < inequality_constraint_count; i++)
107  {
108  if (merged_active_set.end() != std::find(
109  merged_active_set.begin(),
110  merged_active_set.end(),
111  i))
112  continue;
113 
114  const auto lhs = coefficients_mat
115  .multiply_row(i, try_delta_vec);
116 
117  if (lhs > b_vec[i])
118  violated_indices.push_back(i);
119  }
120  }
121 
122  if (violated_indices.empty())
123  {
124  delta_vec = try_delta_vec;
125 
126  auto lagrangian_index = index_not_found;
127  for (size_t i = 0; i < active_set.size(); i++)
128  {
129  const auto lms_i = lagrangian_vec[i];
130  if (lms_i < value_type(0))
131  continue;
132  if (lagrangian_index == index_not_found ||
133  lms_i > lagrangian_vec[lagrangian_index])
134  lagrangian_index = i;
135  }
136 
137  if (lagrangian_index == index_not_found)
138  return;
139 
140  active_set.erase(
141  active_set.begin() +
142  int(lagrangian_index));
143  }
144  else
145  {
146  const auto k_violated = _backtrack(
147  b_vec,
148  coefficients_mat,
149  delta_vec,
150  try_delta_vec,
151  violated_indices,
152  delta_vec);
153 
154  if (k_violated == index_not_found)
155  continue;
156 
157  assert(std::none_of(
158  active_set.begin(),
159  active_set.end(),
160  [k_violated](const size_t index) -> bool
161  {
162  return index == k_violated;
163  }));
164 
165  assert(k_violated != index_not_found);
166  assert(k_violated < b_vec.get_length());
167  active_set.push_back(k_violated);
168  }
169  }
170  }
171 
172  private:
173  // --------------------------------------------------------------------
174  static constexpr auto index_not_found =
175  std::numeric_limits<size_t>::max();
176 
177  static constexpr auto epsilon =
178  value_type(0.000001);
179 
180  // --------------------------------------------------------------------
181  static size_t _backtrack(
182  const matrix_type & b_vec,
183  const matrix_type & coefficients_mat,
184  const matrix_type & old_delta_vec,
185  const matrix_type & violated_delta_vec,
186  const std::vector<size_t> & violated_indices,
187  matrix_type & new_delta_vec)
188  {
189  assert(violated_delta_vec.is_size(old_delta_vec));
190  assert(b_vec.get_length() == coefficients_mat.get_height());
191 
192  assert(std::all_of(
193  violated_indices.begin(),
194  violated_indices.end(),
195  [&b_vec](const size_t i) -> bool
196  {
197  return i < b_vec.get_length();
198  }));
199 
200  auto min_i = index_not_found;
201  auto min_t = value_type(0);
202 
203  const auto diff_vec = violated_delta_vec - old_delta_vec;
204 
205  for (size_t i = 0; i < violated_indices.size(); i++)
206  {
207  const auto violated_index = violated_indices[i];
208  const auto b_scalar = b_vec[violated_index];
209 
210  const auto coefficient_vec =
211  coefficients_mat.copy_row(violated_index);
212 
213  const auto denominator =
214  (coefficient_vec * diff_vec).get_value(0, 0);
215 
216  const auto numerator = b_scalar -
217  (coefficient_vec * old_delta_vec).get_value(0, 0);
218 
219  if (std::fabs(denominator) < epsilon)
220  continue;
221 
222  const auto t = numerator / denominator;
223  if (min_i == index_not_found || t < min_t)
224  {
225  min_i = i;
226  min_t = t;
227  }
228  }
229 
230  if (min_i == index_not_found)
231  return index_not_found;
232 
233  new_delta_vec = (min_t * diff_vec) + old_delta_vec;
234  assert(!new_delta_vec.contains_inf());
235  assert(!new_delta_vec.contains_nan());
236  return violated_indices[min_i];
237  }
238 
239  // --------------------------------------------------------------------
240  static void _create_kkt_mat(
241  const matrix_type & b_vec,
242  const matrix_type & coefficients_mat,
243  const matrix_type & hessian_mat,
244  const matrix_type & derivative_vec,
245  const std::vector<size_t> & active_set,
246  matrix_type & kkt_mat)
247  {
248  const auto K = hessian_mat.get_height();
249  const auto active_count = active_set.size();
250 
251  assert(b_vec.get_length() == coefficients_mat.get_height());
252  assert(hessian_mat.is_square());
253  assert(kkt_mat.get_height() == K + active_count);
254  assert(kkt_mat.get_width() == K + active_count + 1);
255  assert(coefficients_mat.get_width() == K);
256  assert(derivative_vec.is_vector());
257  assert(derivative_vec.get_length() == K);
258 
259  assert(std::all_of(
260  active_set.begin(),
261  active_set.end(),
262  [&b_vec](size_t i) -> bool {
263  return i < b_vec.get_length();
264  }));
265 
266  assert(std::all_of(
267  active_set.begin(),
268  active_set.end(),
269  [&coefficients_mat](size_t i) -> bool {
270  return i < coefficients_mat.get_height();
271  }));
272 
273  kkt_mat.set_values(value_type(0));
274  for (size_t k1 = 0; k1 < K; k1++)
275  for (size_t k2 = 0; k2 < K; k2++)
276  kkt_mat(k1, k2) = hessian_mat(k1, k2);
277 
278  for (size_t i = 0; i < active_count; i++)
279  for (size_t k = 0; k < K; k++)
280  kkt_mat(i + K, k) = kkt_mat(k, i + K) =
281  coefficients_mat(active_set[i], k);
282 
283  const auto right_column_index = kkt_mat.get_width() - 1;
284 
285  for (size_t k = 0; k < K; k++)
286  kkt_mat(k, right_column_index) = -derivative_vec[k];
287 
288  for (size_t i = 0; i < active_count; i++)
289  kkt_mat(K + i, right_column_index) = b_vec[active_set[i]];
290  }
291 
292  // --------------------------------------------------------------------
293  static void _kkt(
294  const matrix_type & b_vec,
295  const matrix_type & coefficients_mat,
296  const matrix_type & hessian_mat,
297  const matrix_type & derivative_vec,
298  const std::vector<size_t> & active_set,
299  matrix_type & delta_vec,
300  matrix_type & lagrangian_vec)
301  {
302  const auto active_count = active_set.size();
303  const auto K = derivative_vec.get_length();
304 
305  assert(delta_vec.is_vector());
306  assert(delta_vec.get_length() == K);
307  assert(lagrangian_vec.is_empty() || lagrangian_vec.is_vector());
308  assert(lagrangian_vec.get_length() == active_count);
309 
310  matrix_type kkt_mat (K + active_count, K + active_count + 1);
311  _create_kkt_mat(
312  b_vec,
313  coefficients_mat,
314  hessian_mat,
315  derivative_vec,
316  active_set,
317  kkt_mat);
318 
319  assert(kkt_mat.all_of([](const value_type v)
320  { return std::isfinite(v); }));
321 
322  kkt_mat.gesv();
323 
324  assert(kkt_mat.all_of([](const value_type v)
325  { return std::isfinite(v); }));
326 
327  const auto last_column = kkt_mat.get_width() - 1;
328 
329  for (size_t k = 0; k < K; k++)
330  delta_vec[k] = kkt_mat(k, last_column);
331 
332  for (size_t i = 0; i < active_count; i++)
333  lagrangian_vec[i] = kkt_mat(K + i, last_column);
334  }
335  };
336 }
337 
338 #endif // JADE_QPAS_HPP__
jade::basic_matrix::get_length
size_t get_length() const
Definition: jade.matrix.hpp:624
jade::basic_matrix::is_vector
bool is_vector() const
Definition: jade.matrix.hpp:925
jade::basic_matrix::get_height
size_t get_height() const
Definition: jade.matrix.hpp:603
jade::basic_matrix::multiply_row
void multiply_row(const size_t row, const value_type value)
Multiplies a row by a specified value.
Definition: jade.matrix.hpp:976
jade::basic_qpas::value_type
TValue value_type
The value type.
Definition: jade.qpas.hpp:23
jade::basic_qpas::loop_over_active_set
static void loop_over_active_set(const matrix_type &b_vec, const matrix_type &coefficients_mat, const matrix_type &hessian_mat, const matrix_type &derivative_vec, const std::vector< size_t > &fixed_active_set, std::vector< size_t > &active_set, matrix_type &delta_vec)
Loops over the active set and computes a delta vector and a new active set.
Definition: jade.qpas.hpp:40
jade::basic_qpas::matrix_type
basic_matrix< value_type > matrix_type
The matrix type.
Definition: jade.qpas.hpp:26
jade::basic_matrix< value_type >
jade::basic_qpas
A template for a class that implements the Quadratic Programming for Active Set algorithm.
Definition: jade.qpas.hpp:20
Copyright © 2015-2020 Jade Cheng