Kernel Quantum Probability Library
The KQP library aims at providing tools for working with quantums probabilities
incremental.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_INCREMENTAL_BUILDER_H__
19 #define __KQP_INCREMENTAL_BUILDER_H__
20 
21 #include <iostream>
22 #include <limits>
23 
24 #include <kqp/cleanup.hpp>
25 #include <kqp/evd_update.hpp>
26 #include <kqp/kernel_evd.hpp>
27 #include <kqp/alt_matrix.hpp>
28 #include <kqp/evd_utils.hpp>
29 
30 #include <kqp/cleaning/unused.hpp>
31 #include <kqp/cleaning/null_space.hpp>
32 #include <kqp/cleaning/qp_approach.hpp>
33 
34 namespace kqp {
35 
36 # include <kqp/define_header_logger.hpp>
37  DEFINE_KQP_HLOGGER("kqp.kernel-evd.incremental");
38 
39 
40 
45  template <typename Scalar> class IncrementalKernelEVD : public KernelEVD<Scalar> {
46  public:
47  KQP_SCALAR_TYPEDEFS(Scalar);
48 
50 
51 #ifndef SWIG
52  // Indirect comparator (decreasing order)
53  struct Comparator {
54  const RealVector &v;
55  Comparator(const RealVector &v) : v(v) {}
56  bool operator()(Index i, Index j) {
57  return this->v(i) > this->v(j);
58  }
59  };
60 #endif
61 
62  IncrementalKernelEVD(const FSpaceCPtr &fs)
63  : KernelEVD<Scalar>(fs), mX(fs->newMatrix()),
64  preImageRatios(std::numeric_limits<Scalar>::infinity(), std::numeric_limits<Scalar>::infinity()) {}
65  virtual ~IncrementalKernelEVD() {}
66 
67  void reset() {
68  mX = this->getFSpace()->newMatrix();
69  mY = ScalarMatrix();
70  mZ = ScalarMatrix();
71  mD = RealVector();
72  // mutable ScalarMatrix k;
73  }
74 
75  void setSelector(const boost::shared_ptr< const Selector<Real> > &selector) {
76  this->selector = selector;
77  }
78 
80  void setPreImagesPerRank(float minimum, float maximum) {
81  this->preImageRatios = std::make_pair(minimum, maximum);
82  }
83 
84 
85  virtual void _add(Real alpha, const FMatrixCPtr &mU, const ScalarAltMatrix &mA) override {
86  // --- Info
87 
88 // KQP_LOG_DEBUG_F(KQP_HLOGGER, "Dimensions: X [%d], Y [%dx%d], Z [%dx%d], D [%d], U [%d], A [%dx%d]",
89 // %mX->size() %mY.rows() %mY.cols() %mZ.rows() %mZ.cols() %mD.rows() %mU.size() %mA.rows() %mA.cols());
90 
91  // --- Pre-computations
92 
93  // Compute W = Y^T X^T
94 
95  ScalarMatrix mW = getFSpace()->k(mX, mY, mU, mA);
96 
97  // Compute V^T V
98  ScalarMatrix vtv = getFSpace()->k(mU, mA);
99  vtv -= mW.adjoint() * mW;
100 
101 
102  // (thin) eigen-value decomposition of V^T V
103  Eigen::SelfAdjointEigenSolver<ScalarMatrix> evd(vtv);
104  ScalarMatrix mQ;
105  ScalarMatrix mQ0;
106  RealVector mDQ;
107  kqp::ThinEVD<ScalarMatrix>::run(evd, mQ, mDQ, &mQ0);
108 
109  Index rank_Q = mQ.cols();
110  mDQ = mDQ.cwiseAbs().cwiseSqrt();
111  for(Index i = 0; i < rank_Q; i++)
112  mQ.col(i) /= mDQ(i);
113 
114  // --- Update
115 
116  // m is the matrix [WQD^1/2 WQ_0; D^1/2 0 ]
117  ScalarMatrix m(mW.rows() + rank_Q, mW.cols());
118 
119  m.topLeftCorner(mW.rows(), rank_Q) = mW * mQ * mDQ.asDiagonal();
120  m.topRightCorner(mW.rows(), mW.cols() - rank_Q) = mW * mQ0;
121 
122  m.bottomLeftCorner(rank_Q, rank_Q) = mDQ.template cast<Scalar>().asDiagonal();
123  m.bottomRightCorner(rank_Q, mW.cols() - rank_Q).setConstant(0) ;
124 
125 
126  for(Index i = 0; i < m.cols(); i++) {
128  // Update
129  ScalarVector v(m.rows());
130  v.head(mZ.cols()) = mZ.adjoint() * m.col(i).head(mZ.cols());
131  v.tail(m.rows() - mZ.cols()) = m.col(i).tail(m.rows() - mZ.cols());
132 
133 
134  // Rank-1 update:
135  // For better accuracy, we don't decrease the rank yet using the selector,
136  // but this might be an option in the future (so as to improve speed)
137  evdRankOneUpdate.update(mD, alpha, v, false, 0, false, result, &mZ);
138 
139  // Take the new diagonal
140  mD = result.mD;
141  }
142 
143  // Update X and Y if the rank has changed
144  if (rank_Q > 0) {
145  // Add the pre-images from U
146  mX->add(mU);
147 
148  // Update mY
149  Index old_Y_rows = mY.rows();
150  Index old_Y_cols = mY.cols();
151  mY.conservativeResize(mY.rows() + mA.rows(), mY.cols() + mQ.cols());
152 
153  mY.bottomRightCorner(mA.rows(), mQ.cols()) = mA * mQ;
154  mY.bottomLeftCorner(mA.rows(), old_Y_cols).setConstant(0);
155  if (old_Y_rows > 0)
156  mY.topRightCorner(old_Y_rows, mQ.cols()) = - mY.topLeftCorner(old_Y_rows, mW.rows()) * mW * mQ;
157  }
158 
159  // (4) Clean-up
160 
161 
162 
163  // --- Rankselection
164  DecompositionList<Real> list(mD);
165  bool identityZ = false;
166  if (this->selector) {
167  // Selects the eigenvalues
168  this->selector->selection(list);
169 
170  // Remove corresponding entries
171  select_rows(list.getSelected(), mD, mD);
172  mY = mY * mZ;
173  select_columns(list.getSelected(), mY, mY);
174 
175  identityZ = true;
176  }
177 
178 
179  // First, tries to remove unused pre-images images
181 
182  // --- Ensure we have a small enough number of pre-images
183  float maxRank = this->preImageRatios.second * (float)mD.rows();
184  if (mX->size() > maxRank) {
185 
186  // Get rid of Z
187  if (!identityZ) mY = mY * mZ;
188  identityZ = true;
189 
190  // Try again to remove unused pre-images
192  KQP_LOG_DEBUG_F(KQP_HLOGGER, "Rank after unused pre-images algorithm: %d [%d]", %mY.rows() %maxRank);
193 
194  // Try to remove null space pre-images
195  ReducedSetNullSpace<Scalar> nullSpace;
196  auto result = nullSpace.run(getFSpace(), mX, mY);
197  mX = result.mX;
198  mY = std::move(result.mY);
199  KQP_LOG_DEBUG_F(KQP_HLOGGER, "Rank after null space algorithm: %d [%d]", %mY.rows() %maxRank);
200 
201  if (mX->size() > maxRank) {
202  if (getFSpace()->canLinearlyCombine()) {
203  // Easy case: we can linearly combine pre-images
204  mX = getFSpace()->linearCombination(mX, mY);
205  mY = Eigen::Identity<Scalar>(mX->size(), mX->size());
206  } else {
207  // Use QP approach
208  ReducedSetWithQP<Scalar> qp_rs;
209  qp_rs.run(this->preImageRatios.first * (float)mD.rows(), getFSpace(), mX, mY, mD);
210 
211  // Get the decomposition
212  mX = qp_rs.getFeatureMatrix();
213  mY = qp_rs.getMixtureMatrix();
214  mD = qp_rs.getEigenValues();
215 
216  // The decomposition is not orthonormal anymore
217  Orthonormalize<Scalar>::run(getFSpace(), mX, mY, mD);
218  }
219  }
220 
221  }
222 
223  if (identityZ)
224  mZ = Eigen::Identity<Scalar>(mY.cols(), mY.cols());
225 
226  }
227 
228  // Gets the decomposition
229  virtual Decomposition<Scalar> _getDecomposition() const override {
230  Decomposition<Scalar> d(this->getFSpace());
231 
232  d.mX = this->mX;
233 
234  const_cast<ScalarMatrix&>(this->mY) = this->mY * this->mZ;
235  const_cast<ScalarMatrix&>(this->mZ) = Eigen::Identity<Scalar>(mY.rows(), mY.rows());
236 
237  d.mY = ScalarAltMatrix(this->mY);
238 
239  d.mD = this->mD;
240 
241  return d;
242  }
243 
244 
245  private:
247  FMatrixPtr mX;
248 
250  ScalarMatrix mY;
251 
253  ScalarMatrix mZ;
254 
256  RealVector mD;
257 
258  // Rank-one EVD update
259  FastRankOneUpdate<Scalar> evdRankOneUpdate;
260 
261  // Used in computation
262  // mutable ScalarMatrix k;
263 
265  boost::shared_ptr< const Selector<Real> > selector;
266 
268  std::pair<float,float> preImageRatios;
269 
270  };
271 
272 }
273 
274 #endif