Kernel Quantum Probability Library
The KQP library aims at providing tools for working with quantums probabilities
evd_utils.hpp
1 /*
2  This file is part of the Kernel Quantum Probability library (KQP).
3 
4  KQP is free software: you can redistribute it and/or modify
5  it under the terms of the GNU General Public License as published by
6  the Free Software Foundation, either version 3 of the License, or
7  (at your option) any later version.
8 
9  KQP is distributed in the hope that it will be useful,
10  but WITHOUT ANY WARRANTY; without even the implied warranty of
11  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  GNU General Public License for more details.
13 
14  You should have received a copy of the GNU General Public License
15  along with KQP. If not, see <http://www.gnu.org/licenses/>.
16  */
17 
18 #ifndef __KQP_KERNEL_EVD_UTILS_H__
19 #define __KQP_KERNEL_EVD_UTILS_H__
20 
21 #include <cassert>
22 
23 #include <boost/scoped_ptr.hpp>
24 #include <boost/type_traits/is_complex.hpp>
25 
26 #include <kqp/kqp.hpp>
27 #include <Eigen/Eigenvalues>
28 
29 #include <kqp/kernel_evd.hpp>
30 #include <kqp/subset.hpp>
31 
32 namespace kqp {
33 
34  template<typename BaseMatrix>
35  struct ThinEVD {
36  typedef typename BaseMatrix::Scalar Scalar;
37  KQP_SCALAR_TYPEDEFS(Scalar);
38 
39  static void run(const Eigen::SelfAdjointEigenSolver<BaseMatrix> &evd,
40  ScalarMatrix &eigenvectors,
41  RealVector &eigenvalues,
42  ScalarMatrix *nullEigenvectors = nullptr,
43  Real threshold = -1) {
44 
45  // We expect eigenvalues to be sorted by increasing order
46  Index dimension = evd.eigenvectors().rows();
47 
48  const Eigen::Matrix<Real,Dynamic,1> &d = evd.eigenvalues();
49  if (threshold < 0)
50  threshold = Eigen::NumTraits<Scalar>::epsilon() * (Real)d.size() * d.cwiseAbs().maxCoeff();
51 
52  Index n = d.rows();
53  Index negatives = 0, zeros = 0;
54 
55  for(Index i = 0; i < n; i++) {
56  assert(i==0 || d[i-1] <= d[i]);
57  if (-d[i] > threshold) negatives++;
58  else if (d[i] < threshold) zeros++;
59  else break;
60  }
61 
62  Index positives = n - negatives - zeros;
63 
64  eigenvalues.resize(positives+negatives);
65  eigenvalues.head(negatives) = d.head(negatives);
66  eigenvalues.tail(positives) = d.tail(positives);
67 
68  eigenvectors.resize(dimension, positives + negatives);
69  eigenvectors.leftCols(negatives) = evd.eigenvectors().leftCols(negatives);
70  eigenvectors.rightCols(positives) = evd.eigenvectors().rightCols(positives);
71 
72  if (nullEigenvectors) {
73  *nullEigenvectors = evd.eigenvectors().block(0, negatives, dimension, zeros);
74  }
75  }
76 
77 
79  static void run(const Eigen::SelfAdjointEigenSolver<BaseMatrix> &evd,
80  ScalarAltMatrix &eigenvectors,
81  RealAltVector &eigenvalues,
82  ScalarAltMatrix *nullEigenvectors = nullptr,
83  Real threshold = -1) {
84 
85  ScalarMatrix _eigenvectors;
86  RealVector _eigenvalues;
87  boost::scoped_ptr<ScalarMatrix> _nullEigenvectors;
88 
89 
90  if (nullEigenvectors != nullptr)
91  _nullEigenvectors.reset(new ScalarMatrix());
92 
93  ThinEVD<ScalarMatrix>::run(evd, _eigenvectors, _eigenvalues, _nullEigenvectors.get(), threshold);
94 
95  eigenvectors.swap(_eigenvectors);
96  eigenvalues.swap(_eigenvalues);
97  if (nullEigenvectors != nullptr)
98  nullEigenvectors->swap(*_nullEigenvectors);
99 
100  }
101 
102  };
103 
104 
105 
106  // template<typename Scalar, typename Cond = void> struct Orthonormalize;
107 
108  template<typename Scalar>
109  struct Orthonormalize {
110  KQP_SCALAR_TYPEDEFS(Scalar);
111 
113  static void run(const FSpaceCPtr &fs, const FMatrixCPtr &mX,
114  ScalarAltMatrix &mY,
115  RealAltVector &mD) {
116  // FIXME: should swap if dense types
117  ScalarMatrix _mY(mY);
118  RealVector _mD(mD);
119  run(fs, mX, _mY, _mD);
120  mY.swap(_mY);
121  mD.swap(_mD);
122  }
123 
124  static void run(const FSpaceCPtr &fs,
125  const FMatrixCPtr &mX,
126  ScalarMatrix &mY,
127  RealVector &mD) {
128 
129  // Negative case: copy what we need
130  RealVector _mD;
131  ScalarMatrix _mY;
132 
133  Index n = 0;
134  std::vector<bool> selection(mD.rows(),false);
135  for(Index j = 0; j < mD.rows(); j++) {
136  if (mD[j] < 0) {
137  selection[j] = true;
138  n++;
139  }
140  }
141 
142  if (n > 0) {
143  select_rows(selection, mD, _mD);
144  select_columns(selection, mY, _mY);
145  }
146 
147  // Perform the EVD
148  ScalarMatrix m = fs->k(mX, mY, RealAltVector(mD.cwiseAbs().cwiseSqrt()));
149 
150  Eigen::SelfAdjointEigenSolver<decltype(m)> evd(m.template selfadjointView<Eigen::Lower>());
151 
152  // A = mY mD^1/2
153  ScalarMatrix mY2;
154  RealVector mD2;
155  kqp::ThinEVD<ScalarMatrix>::run(evd, mY2, mD2);
156  mD2.array() = mD2.array().cwiseAbs(); // just in case of small rounding errors
157  mY *= mD.cwiseAbs().cwiseSqrt().asDiagonal() * mY2 * mD2.cwiseSqrt().cwiseInverse().asDiagonal();
158  mD = std::move(mD2);
159 
160  // Handles negative eigenvalues
161  if (n > 0) {
162  m = mD.template cast<Scalar>().asDiagonal();
163  m.template selfadjointView<Eigen::Lower>().rankUpdate(mY.adjoint() * fs->k(mX) * _mY * _mD.cwiseAbs().cwiseSqrt().asDiagonal() , -2);
164  ScalarMatrix mU;
165  ThinEVD<ScalarMatrix>::run(Eigen::SelfAdjointEigenSolver<ScalarMatrix>(m.template selfadjointView<Eigen::Lower>()), mU, mD);
166  mY *= mU;
167  }
168  }
169 
170 
171 
172 
173  };
174 }
175 
176 #endif