Kernel Quantum Probability Library
The KQP library aims at providing tools for working with quantums probabilities
evd_update.inc.hpp
1 #include <algorithm>
2 #include <bitset>
3 #include <vector>
4 #include <sstream>
5 #include <cmath>
6 #include <limits.h>
7 #include <exception>
8 #include <complex>
9 
10 #include <kqp/kqp.hpp>
11 #include <kqp/logging.hpp>
12 #include <kqp/evd_update.hpp>
13 
14 #include <kqp/define_header_logger.hpp>
15 DEFINE_KQP_HLOGGER("kqp.evd-update");
16 
17 namespace kqp {
18 
19  // Useful Scalar dependant functions
20 
21  inline bool is_real(double) { return true; }
22  inline bool is_real(float) { return true; }
23  template <typename Scalar> inline bool is_real(const std::complex<Scalar>& f) { return std::imag(f) == 0.0; }
24 
25  double real(double f) { return f; }
26  float real(float f) { return f; }
27  template <typename Scalar> inline Scalar real(const std::complex<Scalar>& f) { return std::real(f); }
28 
29 
30 
31  template<typename Scalar> EigenList<Scalar>::~EigenList() {}
32 
39  template <typename Scalar> class IndexedValue {
40  public:
41  typedef typename Eigen::NumTraits<Scalar>::Real Real;
42 
46  size_t newPosition;
47 
51  Real lambda;
52 
56  Real d;
57 
61  Scalar z;
62 
66  std::bitset<2> status;
67 
71  long position;
72 
78  void setSelected(bool value) {
79  status.set(0, value);
80  }
81 
86  bool isRemoved() const {
87  return status[1];
88  }
89 
93  void setRemoved(bool value) {
94  status.set(1, value);
95  }
96 
101  bool isSelected() const {
102  return status[0];
103  }
104 
105 
106  IndexedValue() : position(-1) {
107  }
108 
109  IndexedValue(size_t position, Real d, Scalar z) {
110  this->position = position;
111  this->d = d;
112  this->z = z;
113  this->lambda = d;
114  }
115 
116  std::string toString() {
117  std::stringstream ss;
118  ss << "rank" << newPosition << ", position" << position << ", lambda" << lambda
119  << ", d=" << d << ", z=" << z << ", s=" << isSelected() << ", r=" << isRemoved();
120  return ss.str();
121  }
122  };
123 
124 
125  template<typename Scalar>
127  public:
128  LambdaComparator() {}
129  bool operator() (const IndexedValue<Scalar>* i, const IndexedValue<Scalar>* j) {
130  // Special rule to order after the removed parts
131  int z = (i->isRemoved() ? 1 : 0) - (j->isRemoved() ? 1 : 0);
132  if (z != 0)
133  return z;
134 
135  // We have to invert the order since we want biggest
136  // values
137  // first - that is, returns -1 if o1 > o2
138  return i->lambda > j->lambda;
139  }
140  };
141 
142  template<typename Scalar>
144  DiagonalComparator() {}
145  bool operator() (const IndexedValue<Scalar>* i, const IndexedValue<Scalar>* j) {
146  return i->d > j->d;
147  }
148  };
149 
150 
151 
157  template<typename Scalar> class EigenValues : public EigenList<typename Eigen::NumTraits<Scalar>::Real> {
158  public:
159  typedef typename Eigen::NumTraits<Scalar>::Real Real;
160  std::vector<IndexedValue<Scalar>*>& values;
161  Index minRemoved;
162  Index rank;
163 
164  EigenValues(std::vector<IndexedValue<Scalar>*>& _values) : values(_values) {
165  this->rank = values.size();
166  this->minRemoved = values.size();
167  }
168 
169 
170  Real get(Index index) const {
171  check(index);
172  return values[index]->lambda;
173  }
174 
175 
176  void remove(Index index) {
177  check(index);
178  if (!values[index]->isRemoved()) {
179  minRemoved = std::min(index, minRemoved);
180  values[index]->setRemoved(true);
181  rank--;
182  }
183  }
184 
185 
186  Index size() const {
187  return values.size();
188  }
189 
190 
191  bool isSelected(std::size_t i) const {
192  check(i);
193  return !values[i]->isRemoved();
194  }
195 
196 
197  Index getRank() const {
198  return rank;
199  }
200 
201  private:
202 #ifdef NDEBUG
203  inline void check(std::size_t ) const {}
204 #else
205  inline void check(std::size_t i) const {
206  assert(i < (std::size_t)size());
207  }
208 #endif
209  };
210 
211 
212 
216  template <class Compare, typename Scalar> void sortValues(std::vector<IndexedValue<Scalar>*>& v, size_t from,
217  Compare &compare) {
218 
219  std::sort(v.begin() + from, v.end(), compare);
220  }
221 
222  inline double norm(double x) { return x*x; }
223  inline float norm(float x) { return x*x; }
224  template <typename Scalar> inline Scalar norm(const std::complex<Scalar> &z) { return std::norm(z); }
225 
226 
239  template <typename Scalar> Scalar computeZ(const std::vector<IndexedValue<Scalar>*>& v, int M,
240  double lambda0, int i, const IndexedValue<Scalar>& vi, double di) {
241  double normZ = -(di - lambda0);
242 
243  // lambda_j < di
244  for (int j = i + 1; j < M; j++) {
245  IndexedValue<Scalar> &vj = *v[j];
246  normZ *= (di - vj.lambda) / (di - vj.d);
247  }
248 
249  for (int j = 1; j <= i; j++) {
250  normZ *= (di - v[j]->lambda) / (di - v[j - 1]->d);
251  }
252 
253  // ensures the norm of zi is the same as newz
254  Scalar new_z = vi.z * (Scalar)( sqrt(normZ) / sqrt(kqp::norm(vi.z)));
255  KQP_HLOG_DEBUG( "New z" << convert(i) << " = " << convert(new_z) << " / old was " << convert(vi.z));
256  // newz = kqp::real(vi.z) >= 0 ? sqrt(newz) : -sqrt(newz);
257 
258  return new_z;
259  }
260 
261 
262 
263  template<typename Scalar> class Rotation {
264  public:
265  // The rotation matrix [ c, s; -s, c ]
266  Scalar c, s;
267 
268  // Which is the singular value column we rotated with
269  const IndexedValue<Scalar> *vi, *vj;
270 
271  Rotation(Scalar c, Scalar s, const IndexedValue<Scalar>* _vi, const IndexedValue<Scalar>* _vj): vi(_vi), vj(_vj) {
272  this->c = c;
273  this->s = s;
274  }
275  };
276 
277 
278  template<typename Scalar>
280 
281 
282 
283  template <typename Scalar>
284  void FastRankOneUpdate<Scalar>::update(const Eigen::Matrix<Real,Dynamic,1> & D,
285  double rho, const Eigen::Matrix<Scalar,Dynamic,1> & z,
286  bool computeEigenvectors, const Selector<typename FastRankOneUpdate<Scalar>::Real> *selector, bool keep,
287  EvdUpdateResult<Scalar> &result,
288  Eigen::Matrix<Scalar,Dynamic,Dynamic> * Z) {
289 
290  typedef Eigen::Matrix<Scalar,Dynamic,Dynamic> Matrix;
291  typedef Eigen::Matrix<Scalar,Dynamic,1> Vector;
292 
293 
294  // ---
295  // --- Deflate the matrix (see. G.W. Steward, p. 176)
296  // ---
297 
298  // The deflated diagonal and corresponding z
299  // The matrix we are working on is a (N + 1, N)
300  Index N = z.size();
301  Index rankD = D.rows();
302 
303  KQP_HLOG_DEBUG_F( "EVD rank-one update in dimension %d", %std::max(rankD,N));
304  if (rankD > N)
305  KQP_THROW_EXCEPTION_F(illegal_argument_exception, "D and z are of incompatible dimensions (%d and %d)", %rankD %N);
306 
307 
308  // The algorithm assumes that rho > 0, so that we keep track of this
309  // here, and modify rho so that it is strictly positive
310  bool negativeUpdate = rho < 0;
311  if (negativeUpdate)
312  rho = -rho;
313 
314  // The norm of ||D||
315  double normD = 0;
316 
317  // Our store for singular values, D, and z
318 
319  std::vector<IndexedValue<Scalar> > indexedValues(N); // this one is just for memory allocation
320 
321  typedef std::vector<IndexedValue<Scalar>*> ScalarPtrVector;
322  ScalarPtrVector v(N);
323 
324  // Copy the diagonal entries (possibly inserting zeros if needed)
325  Index offset = N - rankD;
326 
327  // i' is the position on the D diagonal
328  long iprime = rankD - 1;
329  bool foundNonNegative = false;
330  bool toSort = false;
331 
332  for (Index i = N; --i >= 0;) {
333  Real di = 0;
334 
335  if (iprime >= 0 && !is_real(D(iprime)))
336  BOOST_THROW_EXCEPTION(illegal_argument_exception() << errinfo_message("Diagonal value is not real"));
337 
338 
339  // Check the current D[i', i'] if necessary
340  if (iprime >= 0 && offset > 0 && !foundNonNegative)
341  foundNonNegative = kqp::real(D(iprime)) >= 0;
342 
343  size_t position = 0;
344  size_t index = negativeUpdate ? N - 1 - i : i;
345 
346  // If i' points on the first non negative
347  Scalar zpos = 0;
348  if (iprime < 0 || (offset > 0 && foundNonNegative)) {
349  // di is zero in that part
350  // and zi comes from the end of the z vector
351  offset--;
352  position = z.size() + offset - (N - rankD);
353  zpos = z(position);
354  } else {
355  di = kqp::real(D(iprime));
356  di = negativeUpdate ? -di : di;
357  position = i - offset;
358  zpos = z(position);
359  iprime--;
360  }
361 
362  normD += di * di;
363 
364  // Just checking
365  long previousIndex = negativeUpdate ? N - 1 - (i + 1) : i + 1;
366  if (i != N - 1 && (negativeUpdate ^ (v[previousIndex]->d > di))) {
367  toSort = true;
368  }
369 
370  // If the update is negative, we have to reverse the order since
371  // we take the opposite of diagonal entries
372  v[index] = &(indexedValues[index] = IndexedValue<Scalar>(position, di, ((Scalar)sqrt(rho)) * zpos));
373  }
374  normD = sqrt(normD);
375 
376  // Sort if needed
377  if (toSort) {
378  static const DiagonalComparator<Scalar> comparator;
379  sortValues(v, 0, comparator);
380  }
381 
382  // M is the dimension of the deflated diagonal matrix
383  int M = 0;
384 
385  double tau = gamma * epsilon();
386  double tauM2 = tau * normD;
387  double mzNorm = 0;
388 
389  // The list of rotations
390  std::vector<Rotation<Scalar> > rotations;
391 
392  // Deflate the matrix and order the singular values
393  IndexedValue<Scalar> *last = 0;
394  for (int i = 0; i < N; i++) {
395  IndexedValue<Scalar> &vi = *v[i];
396  Scalar zi = vi.z;
397 
398  if (std::abs(zi) <= tauM2) {
399  } else if (M > 0 && (last->d - vi.d <= tauM2)) {
400  double r = sqrt(kqp::norm(last->z) + kqp::norm(zi));
401  rotations.push_back(Rotation<Scalar>(last->z / (Scalar)r, zi / (Scalar)r, last, &vi));
402  last->z = r;
403  vi.z = 0;
404  } else {
405  // Else just copy the values
406  last = &vi;
407  mzNorm += kqp::norm(zi);
408  M++;
409  vi.setSelected(true);
410  }
411  }
412 
413  // Order the array v
414  // so that the first M ones are within the
415  // deflated matrix and the N-M last are outside
416  // and preserving the order
417  // last is the last non selected here
418  int lastFree = -1;
419  if (N != M)
420  for (int i = 0; i < N; i++) {
421  IndexedValue<Scalar> &vi = *v[i];
422  if (!vi.isSelected()) {
423  if (lastFree < 0)
424  lastFree = i;
425  } else if (lastFree >= 0) {
426  // We have some room here
427  v[i] = v[lastFree];
428  v[lastFree] = &vi;
429  if (lastFree + 1 < i)
430  lastFree++;
431  else
432  lastFree = i;
433  }
434  }
435 
436  // ---
437  // --- Search for singular values (solving the secular equation)
438  // ---
439 
440  // For the stopping criterion
441  double e = gamma * epsilon() * M;
442  KQP_HLOG_DEBUG( "Computing " << convert(M) << " eigenvalues");
443  for (int j = 0; j < M; j++) {
444  IndexedValue<Scalar> &svj = *v[j];
445  double diagj = svj.d;
446 
447  double interval = (j == 0 ? mzNorm : v[j - 1]->d - diagj) / 2;
448  double middle = diagj + interval;
449 
450  // Stopping criteria from Gu & Eisenstat
451  double psi = 0;
452  double phi = 0;
453 
454  // "Searching for singular value between %e and %e (interval %e)",
455  // diagj, diagj + interval * 2, interval);
456 
457  double nu = -interval;
458  double f = -1;
459  do {
460  // Update nu
461  // TODO enhance the root finder by using a better
462  // approximation
463  double oldnu = nu;
464  if (f < 0)
465  nu += interval;
466  else
467  nu -= interval;
468 
469  if (nu == oldnu) {
470  // Stopping since we don't change f anymore
471  break;
472  }
473 
474  // Compute the new phi, psi and f
475  psi = phi = 0;
476 
477  // lambda is between diagj and (diagj1 + diagj)/2
478  for (int i = j; i < M; i++) {
479  IndexedValue<Scalar> &vi = *v[i];
480  psi += kqp::norm(vi.z) / (vi.d - middle - nu);
481  }
482 
483  for (int i = 0; i < j; i++) {
484  IndexedValue<Scalar> &vi = *v[i];
485  phi += kqp::norm(vi.z) / (vi.d - middle - nu);
486  }
487 
488  f = 1 + psi + phi;
489 
490  interval /= 2;
491  } while (std::abs(f) == std::numeric_limits<double>::infinity()
492  || std::abs(f) > (1 + std::abs(psi) + std::abs(phi)) * e);
493 
494  // Done, store the eigen value
495  KQP_HLOG_DEBUG( "Eigenvalue: old = " << convert(svj.lambda) << ", computed = " << convert(middle + nu));
496  svj.lambda = middle + nu;
497 
498  // Because of rounding errors, that can happen
499  if (svj.lambda < diagj) {
500  svj.lambda = diagj;
501  } else {
502  double max = j == 0 ? mzNorm + diagj : v[j - 1]->d;
503  if (svj.lambda > max) {
504  svj.lambda = max;
505  }
506  }
507  }
508 
509  // }
510 
511  // ---
512  // --- Compute the singular vectors
513  // ---
514 
515  // First, recompute z to match the singular values we have
516 
517  double lambda0 = v.empty() ? 0 : v[0]->lambda;
518 
519  for (int i = 0; i < M; i++) {
520  IndexedValue<Scalar> vi = *v[i];
521  double di = vi.d;
522  Scalar newz = computeZ(v, M, lambda0, i, vi, di);
523 
524  // Remove z too close to 0
525  if (std::abs(newz) < tauM2) {
526  v[i]->setSelected(false);
527  }
528  vi.z = newz;
529  }
530 
531  // --- Let's construct the result ---
532 
533  // --- First, take the opposite of eigenvalues if
534  // --- we are doing a negative update
535  if (negativeUpdate)
536  for(typename ScalarPtrVector::iterator iv = v.begin(); iv != v.end(); iv++) {
537  (*iv)->d = -(*iv)->d;
538  (*iv)->lambda = -(*iv)->lambda;
539  }
540 
541  // --- Set eigen values (and the rank)
542 
543  // Select the eigenvalues if needed
544  static const LambdaComparator<Scalar> lambdaComparator;
545  std::sort(v.begin(), v.end(), lambdaComparator);
546 
547  Index rank = v.size();
548  if (selector) {
549  EigenValues<Scalar> list(v);
550  selector->selection(list);
551  rank = list.rank;
552 
553  // Reorder if needed
554  if (rank < N && (list.minRemoved != rank)) {
555  std::sort(v.begin() + list.minRemoved, v.end(), lambdaComparator);
556  }
557 
558  }
559 
560  // then store them,
561  int nbSelected = 0;
562  int nbNaN = 0;
563  KQP_HLOG_DEBUG_F( "Final rank is %d", %rank);
564  result.mD = RealVector(rank);
565  for (Index i = 0; i < rank; i++) {
566  v[i]->newPosition = i;
567  result.mD(i) = v[i]->lambda;
568  if (v[i]->isSelected())
569  nbSelected++;
570  if (isNaN(v[i]->lambda))
571  nbNaN++;
572  }
573 
574  if (nbNaN > 0)
575  KQP_THROW_EXCEPTION_F(arithmetic_exception, "We had %d eigen value(s) that is/are NaN", %nbNaN);
576 
577  // --- Compute the eigenvectors
578 
579  if (computeEigenvectors || Z) {
580  // Creates the matrix
581  result.mQ.setZero(N, rank);
582  Matrix &Q = result.mQ;
583 
584  // Set the new values: work eigenvector by eigenvector (indexed by j)
585  for (int j = 0; j < rank; j++) {
586  IndexedValue<Scalar> &vj = *v[j];
587  if (!vj.isSelected()) {
588  Q(vj.position, j) = 1;
589  } else {
590  // Compute the new vector
591  double columnNorm = 0;
592  int iM = 0;
593  for (Index i = 0; i < N && iM < M; i++) {
594  IndexedValue<Scalar> &vi = *v[i];
595  if (vi.isSelected()) {
596  double di = vi.d;
597  Scalar x = vi.z / (Scalar)(di - vj.lambda);
598  columnNorm += kqp::norm(x);
599  Q(vi.position, j) = x;
600 
601  if (isNaN(x))
602  KQP_THROW_EXCEPTION(arithmetic_exception, "NaN value in matrix Q");
603  iM++;
604  }
605  }
606 
607  // Normalize
608  Q.col(j) /= std::sqrt(columnNorm);
609  }
610  }
611 
612  // --- Rotate the vectors that need to be rotated
613  for (size_t r = 0; r < rotations.size(); r++) {
614  Rotation<Scalar> &rot = rotations[r];
615  size_t i = rot.vi->position;
616  size_t j = rot.vj->position;
617  // TODO: use Eigen rotation
618 
619  // Rotation only affect the two rows i and j
620  for (int col = 0; col < rank; col++) {
621  Scalar x = Q(i,col);
622  Scalar y = Q(j,col);
623  Q(i, col) = x * rot.c - y * rot.s;
624  Q(j, col) = x * rot.s + y * rot.c;
625  }
626  }
627 
628  if (!keep && rank < N)
629  Q = Q.block(0,0,rank,rank);
630 
631  if (Z) {
632  if (Z->cols() < rank) {
633  Index z_cols = Z->cols();
634  Index z_rows = Z->rows();
635 
636  Index diff = rank - z_cols;
637  Z->conservativeResize(Z->rows() + diff, Z->cols() + diff);
638 
639  Z->bottomLeftCorner(diff, z_cols).setConstant(0);
640  Z->bottomRightCorner(diff, diff).setIdentity(diff, diff);
641  Z->topRightCorner(z_rows, diff).setConstant(0);
642  }
643  *Z = ((*Z) * Q).eval();
644 
645  }
646 
647  }
648 
649 
650  }
651 
652  //explicit instantiation of
653 #define RANK_ONE_UPDATE(Scalar) template class FastRankOneUpdate<Scalar>; template class EigenList<Scalar>;
654 
655  RANK_ONE_UPDATE(double);
656  RANK_ONE_UPDATE(float);
657  RANK_ONE_UPDATE(std::complex<double>);
658  RANK_ONE_UPDATE(std::complex<float>);
659 }
660 
661