Kernel Quantum Probability Library
The KQP library aims at providing tools for working with quantums probabilities
null_space.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_REDUCED_SET_NULL_SPACE_H__
19 #define __KQP_REDUCED_SET_NULL_SPACE_H__
20 
21 #include <algorithm> // sort
22 
23 #include <kqp/kqp.hpp>
24 #include <Eigen/Eigenvalues>
25 
26 #include <kqp/feature_matrix.hpp>
27 #include <kqp/subset.hpp>
28 #include <kqp/cleaning/unused.hpp>
29 
30 namespace kqp {
31 # include <kqp/define_header_logger.hpp>
32  DEFINE_KQP_HLOGGER("kqp.cleaning.null_space");
33 
34 #ifndef SWIG
35  template <class ComparableArray, typename Index = int>
36  struct IndirectSort {
37  const ComparableArray &array;
38  IndirectSort(const ComparableArray &array) : array(array) {}
39  bool operator() (int i,int j) { return (array[i] < array[j]);}
40  };
41 
42 
43  template <class ComparableArray, typename Index = int>
44  struct AbsIndirectSort {
45  const ComparableArray &array;
46  AbsIndirectSort(const ComparableArray &array) : array(array) {}
47  bool operator() (int i,int j) { return std::abs(array[i]) < std::abs(array[j]);}
48  };
49 
50 
51  template<typename Scalar>
53  KQP_SCALAR_TYPEDEFS(Scalar);
54  FMatrixPtr mX;
55  ScalarAltMatrix mY;
56  };
57 
58  template<typename Scalar>
60  public:
61  KQP_SCALAR_TYPEDEFS(Scalar);
62  private:
63  double m_epsilon;
64  std::pair<float,float> m_preImageRatios;
65  std::pair<Index,Index> m_preImagesRange;
66 
67  public:
69  m_epsilon(Eigen::NumTraits<Real>::epsilon()),
70  m_preImageRatios(std::make_pair(0, std::numeric_limits<Real>::infinity())),
71  m_preImagesRange(std::make_pair(0, std::numeric_limits<Index>::max()))
72  {}
73 
74  ReducedSetNullSpace & epsilon(double _epsilon) { this->m_epsilon = _epsilon; return *this; }
75 
77  void setPreImagesPerRank(float reset, float maximum) {
78  this->m_preImageRatios = std::make_pair(reset, maximum);
79  }
80 
81  void setRankRange(Index reset, Index maximum) {
82  this->m_preImagesRange = std::make_pair(reset, maximum);
83  }
97  FMatrixPtr remove(const FMatrixCPtr &mF, ScalarMatrix &kernel, Eigen::PermutationMatrix<Dynamic, Dynamic, Index>& mP, const RealVector &weights) const {
98  typedef typename Eigen::PermutationMatrix<Dynamic, Dynamic, Index> Permutation;
99 
100  // FIXME: should be wiser
101  double delta = 1e-4;
102 
103  // --- Check if we have something to do
104  if (mF->size() == 0)
105  return mF->copy();
106 
107  // --- Look up at the indices to remove
108 
109  // Get the pre-images available for removal (one per vector of the null space, i.e. by columns in kernel)
110  // Summarize the result per pre-image (number of vectors where it appears)
111  // Eigen::Matrix<int,Dynamic,1> num_vectors = (kernel.array().abs() > kernel.array().abs().colwise().maxCoeff().replicate(kernel.rows(), 1) * delta).template cast<int>().rowwise().sum();
112 
113  std::vector<Index> to_remove;
114  for(Index i = 0; i < /*num_vectors*/kernel.rows(); i++)
115  to_remove.push_back(i);
116 
117  // Sort
118  const Index pre_images_count = kernel.rows();
119  const std::size_t remove_size = kernel.cols();
120  const Index keep_size = mF->size() - remove_size;
121 
122  mP = Permutation(mF->size());
123  mP.setIdentity();
124 
125  std::sort(to_remove.begin(), to_remove.end(), IndirectSort<RealVector>(weights));
126  std::vector<bool> selection(mF->size(), true);
127  std::vector<bool> used(to_remove.size(), false);
128 
129  ScalarVector v;
130 
131  // --- Remove the vectors one by one (Pivoted Gauss)
132  Index last = mF->size() - 1;
133 
134  Index remaining = remove_size;
135 
136  while (remaining > 0) {
137  // Select the pre-image to remove along with the null space vector
138  // that will be used
139  // TODO: find a better way
140 
141  Eigen::Matrix<Real,Dynamic,1> colnorms = kernel.colwise().norm();
142 
143  size_t i = 0;
144  Index j = -1;
145  Real max = 0;
146  Real threshold;
147  for(size_t index = 0; index < to_remove.size(); index++) {
148  // remove the ith pre-image
149  i = to_remove[index];
150  if (!selection[i]) continue;
151 
152  // Searching for the highest magnitude
153  j = -1;
154  // ... but we want it above a given threshold
155  max = 0;
156  for(Index k = 0; k < kernel.cols(); k++) {
157  Real x = kernel.array().abs()(i, k);
158  if (!used[k] && x > max && x > delta * colnorms(k)) {
159  threshold = delta * colnorms(k); // just for debug
160  max = x;
161  j = k;
162  }
163  }
164  if (j >= 0) break;
165  }
166  if (j == -1)
167  KQP_THROW_EXCEPTION_F(assertion_exception,
168  "Could not find a way to remove the a pre-image with null space (%d/%d pre-images). %d remaining.",
169  %kernel.cols() %kernel.rows() %remaining);
170  KQP_HLOG_DEBUG_F("Selected pre-image %d with basis vector %d [%g > %g; norm=%g]", %i %j %max %threshold %colnorms(j));
171 
172  remaining--;
173  used[j] = true;
174 
175  // Update permutation by putting this vector at the end
176  selection[i] = false;
177  mP.indices()(i) = j + keep_size;
178  last -= 1;
179 
180  // Update the matrix
181  Scalar kij = kernel(i,j);
182  KQP_HLOG_DEBUG_F("Normalizing column %d [norm %g] with the inverse of %g", %j %kernel.col(j).norm() %kij)
183  v = kernel.col(j) / kij;
184  kernel.col(j) /= -kij;
185 
186  assert(!kqp::isNaN(kernel.col(j).squaredNorm()));
187 
188  kernel(i,j) = 0;
189 
190  kernel = ((Eigen::Identity<Scalar>(pre_images_count, pre_images_count)
191  - v * ScalarVector::Unit(pre_images_count, i).adjoint()) * kernel).eval();
192  }
193 
194 
195  // --- Remove the vectors from mF and set the permutation matrix
196 
197  // Remove unuseful vectors
198 
199  select_rows(selection.begin(), selection.end(), kernel, kernel);
200 
201  // Complete the permutation matrix
202  Index count = 0;
203  for(size_t index = 0; index < selection.size(); index++)
204  if (selection[index]) {
205  mP.indices()(index) = count++;
206  }
207 
208  return mF->subset(selection.begin(), selection.end());
209  }
210 
218  ReducedSetNullSpaceResult<Scalar> run(const FSpaceCPtr &fs, const FMatrixCPtr &mF, const ScalarAltMatrix &mY) const {
219  ReducedSetNullSpaceResult<Scalar> result;
220  result.mX = mF->copy();
221  result.mY = mY;
222 
223  auto & _mX = result.mX;
224  auto & _mY = result.mY;
225 
226  // Dimension of the problem
227  const Index N = _mY.rows();
228  const Index rank = _mY.cols();
229  assert(N == _mX->size());
230 
231  // Select with max rank and threshold
232  Index target = N;
233  if (N > this->m_preImageRatios.second * (Real)rank)
234  target = std::min(target, (Index)(this->m_preImageRatios.first * (Real)rank));
235 
236  if (N > this->m_preImagesRange.second)
237  target = std::min(target, this->m_preImagesRange.first);
238 
239  // Don't run the EVD if possible
240  if (N == target && m_epsilon == 0)
241  return std::move(result);
242 
243  // Removes unused pre-images
244  CleanerUnused<Scalar>::run(_mX, _mY);
245 
246  // EVD
247  Eigen::SelfAdjointEigenSolver<ScalarMatrix> evd(fs->k(_mX));
248  const Eigen::Matrix<Real,Dynamic,1> &d = evd.eigenvalues();
249 
250  std::vector<Index> list(N);
251  for(Index i = 0; i < N; i++)
252  list[i] = i;
253  std::sort(list.begin(), list.end(), AbsIndirectSort< Eigen::Matrix<Real,Dynamic,1>, Index >(d));
254 
255 
256  Real threshold = m_epsilon * (Real)d.size() * std::abs(d[list.back()]);
257  Index nullSize = N - target;
258  while (nullSize < N && std::abs(d[nullSize]) < threshold) {
259  nullSize++;
260  }
261 
262  KQP_HLOG_DEBUG_F("Rank of used null space is %d (image %d)", %nullSize %(N-nullSize));
263 
264  assert(nullSize >= 0);
265  if (nullSize <= 0)
266  return std::move(result);
267 
268  // True if the null space is not really null
269  bool forced = std::abs(d[nullSize-1]) > threshold;
270  KQP_HLOG_DEBUG_F("Last eigenvalue is %g (threshold %g): forced %d", %std::abs(d[nullSize-1]) %threshold %forced);
271 
272  // Compute the kernel
273  ScalarMatrix kernel(_mX->size(), nullSize);
274  for(Index i = 0; i < nullSize; i++)
275  kernel.col(i) = evd.eigenvectors().col(list[i]);
276 
277  // Remove pre-images using the kernel
278  RealVector weights = _mY.rowwise().squaredNorm().array() * fs->k(_mX).diagonal().array().abs();
279  Eigen::PermutationMatrix<Dynamic, Dynamic, Index> mP;
280  _mX = remove(_mX, kernel, mP, weights);
281 
282  if (forced) {
283  // Compute _mY so that _mX * _mY is orthonormal, ie _mY' _mX' _mX _mY is the identity
284  Eigen::SelfAdjointEigenSolver<ScalarMatrix> evd(fs->k(_mX).template selfadjointView<Eigen::Lower>());
285  ScalarMatrix __mY = std::move(evd.eigenvectors());
286  __mY *= evd.eigenvalues().cwiseAbs().cwiseSqrt().cwiseInverse().asDiagonal();
287  __mY *= fs->k(_mX, __mY, mF, mY);
288  _mY = __mY;
289  } else {
290  // Y <- (Id A) P Y
291  ScalarMatrix mY2(_mY);
292  mY2= (mP * mY2).topRows(_mX->size()) + kernel * (mP * mY2).bottomRows(_mY.rows() - _mX->size());
293 
294  _mY.swap(mY2);
295  }
296  // Removes unused pre-images
297 
298  CleanerUnused<Scalar>::run(_mX, _mY);
299  KQP_HLOG_DEBUG_F("Finished null space cleaning: number of pre-images is %d", %_mX->size());
300  return std::move(result);
301  }
302 
303  };
304 
305 #endif
306 
307  template<typename Scalar> class CleanerNullSpace: public Cleaner<Scalar> {
308  ReducedSetNullSpace<Scalar> m_cleaner;
309  public:
310  typedef typename Eigen::NumTraits<Scalar>::Real Real;
311 
312  CleanerNullSpace() {
313  }
314 
315  void epsilon(double _epsilon) {
316  m_cleaner.epsilon(_epsilon);
317  }
318 
319  void setPreImagesPerRank(float reset, float maximum) {
320  if (reset < 0)
321  KQP_THROW_EXCEPTION_F(kqp::illegal_argument_exception, "Reset rank (%g) is below 0", %reset);
322  if (maximum < reset)
323  KQP_THROW_EXCEPTION_F(kqp::illegal_argument_exception, "Maximum rank (%g) is out of bounds (should be greater than reset %g)", %maximum %reset);
324  m_cleaner.setPreImagesPerRank(reset, maximum);
325  }
326 
327  void setRankRange(Index reset, Index maximum) {
328  if (reset < 0)
329  KQP_THROW_EXCEPTION_F(kqp::illegal_argument_exception, "Reset rank (%d) is below 0", %reset);
330  if (maximum < reset)
331  KQP_THROW_EXCEPTION_F(kqp::illegal_argument_exception, "Maximum rank (%d) is out of bounds (should be greater than reset %d)", %maximum %reset);
332  m_cleaner.setRankRange(reset, maximum);
333  }
334 
335  virtual void cleanup(Decomposition<Scalar> &d) const {
336  auto result = m_cleaner.run(d.fs, d.mX, d.mY);
337  d.mX = std::move(result.mX);
338  d.mY = std::move(result.mY);
339  }
340 
341  virtual picojson::value save() const override {
342  picojson::object json;
343  json["name"] = picojson::value("null-space");
344  json["scalar"] = picojson::value(ScalarInfo<Scalar>::name());
345  json["epsilon"] = picojson::value(m_epsilon);
346  return picojson::value(json);
347  }
348 
349  private:
350  Real m_epsilon;
351  };
352 
353 }
354 
355 #endif