Kernel Quantum Probability Library
The KQP library aims at providing tools for working with quantums probabilities
qp_approach.inc.hpp
1 #include <kqp/coneprog.hpp>
2 
3 #include <kqp/cleaning/qp_approach.hpp>
4 
5 #include <Eigen/Cholesky>
6 
7 #include <kqp/define_header_logger.hpp>
8 DEFINE_KQP_HLOGGER("kqp.qp-approach");
9 
10 namespace kqp {
11 
12  template<typename Derived>
13  bool isnan(const Eigen::MatrixBase<Derived>& x)
14  {
15  return !(x.array() == x.array()).all();
16  }
17 
18  // --- QP solver
19 
20  template<typename Scalar>
21  class KQP_KKTSolver : public cvxopt::KKTSolver<Scalar> {
22  const Eigen::LLT<KQP_MATRIX(Scalar)> &cholK;
23  const KQP_MATRIX(Scalar) &B;
24  const KQP_MATRIX(Scalar) &BBT;
25  const KQP_VECTOR(Scalar) &nu;
26 
27  // Number of pre-images
28  int n;
29  // Rank (number of basis vectors)
30  int r;
31 
32  // Wdi
33  KQP_VECTOR(Scalar) Wd;
34 
35  std::vector<Eigen::LLT<KQP_MATRIX(Scalar)> > L22, L33;
36  std::vector<KQP_MATRIX(Scalar)> L32, L42, L43;
37  Eigen::LLT<KQP_MATRIX(Scalar)> L44;
38  public:
39  KQP_KKTSolver(const Eigen::LLT<KQP_MATRIX(Scalar)> &cholK, const KQP_MATRIX(Scalar) &B, const KQP_MATRIX(Scalar) &BBT, const cvxopt::ScalingMatrix<Scalar> &w, const KQP_VECTOR(Scalar) &nu)
40  : cholK(cholK), B(B), BBT(BBT), nu(nu), n(B.cols()), r(w.d.size() / n / 2) , Wd(w.d)
41  {
42  // The scaling matrix has dimension r * n * 2
43  KQP_HLOG_DEBUG( "Preparing the KKT solver of dimension r=" << convert(r) << " and n=" << convert(n));
44 
45  // Gets U and V
46  KQP_VECTOR(Scalar) U = w.d.topRows(r*n);
47  U.array() = U.array() * U.array();
48 
49  KQP_VECTOR(Scalar) V = w.d.bottomRows(r*n);
50  V.array() = V.array() * V.array();
51 
52  // Computes L22[i]
53  L22.resize(r);
54  for(int i = 0; i < r; i++) {
55  KQP_MATRIX(Scalar) mX = nu[i] * nu[i] * BBT;
56  mX.diagonal() += U.segment(i*n, n);
57  L22[i].compute(mX);
58  if (L22[i].info() != Eigen::ComputationInfo::Success)
59  KQP_THROW_EXCEPTION_F(arithmetic_exception, "Error [%d] in computing L22[%d]", %L22[i].info() %i);
60  }
61 
62 
63  // Computes L32[i]
64  L32.resize(r);
65  for(int i = 0; i < r; i++) {
66  // Solves L32 . L22' = - B . B^\top
67  L32[i] = -nu[i] * nu[i] * BBT;
68  L22[i].matrixU().template solveInPlace<Eigen::OnTheRight>(L32[i]);
69  if (isnan(L32[i])) KQP_THROW_EXCEPTION_F(arithmetic_exception, "NaN in L32[%d]", %i);
70  }
71 
72  // Computes L33
73  L33.resize(r);
74  for(int i = 0; i < r; i++) {
75  // Cholesky decomposition of V + BBT - L32 . L32.T
76  KQP_MATRIX(Scalar) mX = nu[i] * nu[i] * BBT - L32[i] * L32[i].adjoint();
77  mX.diagonal() += V.segment(i*n, n);
78  L33[i].compute(mX);
79  if (L33[i].info() != Eigen::ComputationInfo::Success)
80  KQP_THROW_EXCEPTION_F(arithmetic_exception, "Error [%d] in computing L33[%d]", %L33[i].info() %i);
81  }
82 
83  // Computes L42
84  L42.resize(r);
85  for(int i = 0; i < r; i++) {
86  // Solves L42 L22' = Id
87  L42[i].setIdentity(n,n);
88  L22[i].matrixU().template solveInPlace<Eigen::OnTheRight>(L42[i]);
89  if (isnan(L42[i])) KQP_THROW_EXCEPTION_F(arithmetic_exception, "NaN in L42[%d]", %i);
90  }
91 
92 
93  // Computes L43
94  L43.resize(r);
95  for(int i = 0; i < r; i++) {
96  // Solves L43 L33' = Id - L42 L32'
97  L43[i].noalias() = - L42[i] * L32[i].adjoint();
98  L43[i].diagonal().array() += 1;
99  L33[i].matrixU().template solveInPlace<Eigen::OnTheRight>(L43[i]);
100  if (isnan(L43[i])) KQP_THROW_EXCEPTION_F(arithmetic_exception, "NaN in L43[%d]", %i);
101  }
102 
103  // Computes L44: Cholesky of
104  KQP_MATRIX(Scalar) _L44(n,n);
105  _L44.setConstant(0.);
106  for(int i = 0; i < r; i++) {
107  _L44 += L43[i] * L43[i].adjoint() + L42[i] * L42[i].adjoint();
108  }
109  L44.compute(_L44);
110  if (L44.info() != Eigen::ComputationInfo::Success)
111  KQP_THROW_EXCEPTION_F(arithmetic_exception, "Error [%d] in computing L44", %L44.info());
112 
113  }
114 
116  void reconstruct(KQP_MATRIX(Scalar) &L) const {
117  L.resize(n*(r+1) + 2*n*r, n*(r+1) + 2*n*r);
118 
119  for(int i = 0; i < r; i++) {
120  L.block(i*n, i*n, n, n) = cholK.matrixL();
121  L.block((i+r)*n, (i+r)*n, n,n) = L22[i].matrixL();
122  L.block((i+r)*n, i*n, n,n) = - nu[i] * B;
123  L.block((i+2*r)*n, i*n, n,n) = nu[i] * B;
124 
125  L.block((i+2*r)*n, (i+r)*n, n,n) = L32[i];
126  L.block((i+2*r)*n, (i+2*r)*n, n,n) = L33[i].matrixL();
127 
128  L.block((3*r)*n, (i+r)*n, n,n) = L42[i];
129  L.block((3*r)*n, (i+2*r)*n, n,n) = L43[i];
130  }
131  L.block((3*r)*n, (3*r)*n, n,n) = L44.matrixL();
132 
133 
134  Eigen::DiagonalMatrix<double, Dynamic> D(L.rows());
135  D.diagonal().setConstant(1.);
136  D.diagonal().segment(r*n,2*r*n).setConstant(-1);
137 
138  L = (L * D * L.adjoint()).eval();
139  }
140 
141 
142  virtual void solve(KQP_VECTOR(Scalar) &x, KQP_VECTOR(Scalar) &/*y*/, KQP_VECTOR(Scalar) & z) const {
143  if (isnan(x)) KQP_THROW_EXCEPTION(arithmetic_exception, "NaN in x");
144  if (isnan(z)) KQP_THROW_EXCEPTION(arithmetic_exception, "NaN in z");
145 
146  // Prepares the access to sub-parts
147  Eigen::VectorBlock<KQP_VECTOR(Scalar)> b = x.segment(r*n,n);
148 
149  // First phase
150  for(int i = 0; i < r; i++) {
151  Eigen::VectorBlock<KQP_VECTOR(Scalar)> ai = x.segment(i*n, n);
152  Eigen::VectorBlock<KQP_VECTOR(Scalar)> ci = z.segment(i*n, n);
153  Eigen::VectorBlock<KQP_VECTOR(Scalar)> di = z.segment((i+r)*n, n);
154 
155  cholK.matrixL().solveInPlace(ai);
156 
157  // Solves L22 x = - B * ai - ci
158  KQP_MATRIX(Scalar) Bai = nu[i] * B * ai;
159 
160  ci *= -1.;
161  ci -= Bai;
162  L22[i].matrixL().solveInPlace(ci);
163 
164  di *= -1.;
165  di += Bai - L32[i] * ci;
166  L33[i].matrixL().solveInPlace(di);
167 
168  b += L42[i] * ci + L43[i] * di;
169  }
170  L44.matrixL().solveInPlace(b);
171 
172  // Second phase
173  L44.matrixU().solveInPlace(b);
174  for(int i = 0; i < r; i++) {
175  Eigen::VectorBlock<KQP_VECTOR(Scalar)> ai = x.segment(i*n, n);
176  Eigen::VectorBlock<KQP_VECTOR(Scalar)> ci = z.segment(i*n, n);
177  Eigen::VectorBlock<KQP_VECTOR(Scalar)> di = z.segment((i+r)*n, n);
178 
179  di -= L43[i].adjoint() * b;
180  L33[i].matrixU().solveInPlace(di);
181 
182 
183  ci -= L32[i].adjoint() * di + L42[i].adjoint() * b;
184  L22[i].matrixU().solveInPlace(ci);
185 
186  ai += B.adjoint() * nu[i] * (ci - di);
187  cholK.matrixU().solveInPlace(ai);
188  }
189 
190 
191  // Scale z
192  z.array() *= Wd.array();
193 
194  // Detect NaN
195 
196 
197  if (isnan(x)) KQP_THROW_EXCEPTION(arithmetic_exception, "NaN in solved x");
198  if (isnan(z)) KQP_THROW_EXCEPTION(arithmetic_exception, "NaN in solved z");
199  }
200 
201 
202 
203  };
204 
206  template<typename Scalar>
207  KQP_KKTPreSolver<Scalar>::KQP_KKTPreSolver(const KQP_MATRIX(Scalar)& gramMatrix, const KQP_VECTOR(KQP_REAL_OF(Scalar)) &nu) :
208  // Compute the Cholesky decomposition of the gram matrix
209  nu(nu)
210  {
211  Index n =gramMatrix.rows();
212 
213  // Computes B in B A' = G0 (L21 and L31)
214  // i.e. computes A B' = G0
215  if (boost::is_complex<Scalar>::value) {
216  // FIXME: optimise using 2 Cholesky decompositions (see TR)
217  KQP_MATRIX(Real) dGram;
218  dGram.bottomRightCorner(n,n) = dGram.topLeftCorner(n,n) = gramMatrix.real();
219  dGram.bottomLeftCorner(n,n) = dGram.topRightCorner(n,n) = -gramMatrix.imag();
220  lltOfK.compute(dGram);
221 
222  // Complex case, G0 = [Id, Id; Id, -Id]
223  B.resize(2*n, 2*n);
224  auto Idn = Eigen::Matrix<Real,Dynamic,Dynamic>::Identity(n,n);
225  B.topLeftCorner(n,n) = Idn;
226  B.topRightCorner(n,n) = Idn;
227  B.bottomLeftCorner(n,n) = Idn;
228  B.bottomRightCorner(n,n) = -Idn;
229  } else {
230  lltOfK.compute(gramMatrix.real());
231  B.setIdentity(gramMatrix.rows(), gramMatrix.cols());
232  }
233 
234  lltOfK.matrixU().template solveInPlace<Eigen::OnTheRight>(B);
235 
236  // Computing B * B.T
237  BBT.noalias() = B * B.adjoint();
238  }
239 
240  template<typename Scalar>
241  cvxopt::KKTSolver<typename KQP_KKTPreSolver<Scalar>::Real> *KQP_KKTPreSolver<Scalar>::get(const cvxopt::ScalingMatrix<typename KQP_KKTPreSolver<Scalar>::Real> &w) {
242  KQP_HLOG_DEBUG( "Creating a new KKT solver");
243  return new KQP_KKTSolver<Real>(lltOfK, B, BBT, w, nu);
244  }
245 
246 
267  template<typename Scalar>
268  class KMult : public cvxopt::QPMatrix<KQP_REAL_OF(Scalar)> {
269  Index n, r;
270  const KQP_MATRIX(Scalar) &g;
271  public:
272  typedef KQP_REAL_OF(Scalar) Real;
273 
274  KMult(Index n, Index r, const KQP_MATRIX(Scalar) &g) : n(n), r(r), g(g) {}
275 
276  virtual void mult(const KQP_VECTOR(Real) &x, KQP_VECTOR(Real) &y, bool /*trans*/ = false) const {
277  y.resize(x.rows());
278  y.tail(n).setZero();
279 
280  if (boost::is_complex<Scalar>::value) {
281  for(int i = 0; i < r; i++) {
282  const int i_re = 2 * i, i_im = 2 * i + 1;
283  y.segment(i_re*n, n).noalias() = g.real() * x.segment(i_re*n,n) - x.segment(i_im*n,n);
284  y.segment(i_im*n, n).noalias() = - g.imag() * x.segment(i_re*n,n) + g.real() * x.segment(i_im*n,n);
285  }
286  } else {
287  for(int i = 0; i < r; i++)
288  y.segment(i*n, n).noalias() = g.real() * x.segment(i*n,n);
289  }
290  }
291  virtual Index rows() const { return g.rows(); }
292  virtual Index cols() const { return g.cols(); }
293 
294  };
295 
313  template<typename Scalar>
314  class QPConstraints : public cvxopt::QPMatrix<KQP_REAL_OF(Scalar)> {
315  public:
316  typedef KQP_REAL_OF(Scalar) Real;
317  private:
318  Index n, r;
319  const KQP_VECTOR(Real) &nu;
320 
321  public:
322  QPConstraints(Index n, Index r, const KQP_VECTOR(Real) &nu) : n(n), r(r), nu(nu) {}
323 
325  inline auto get(const KQP_VECTOR(Real) &x, Index i) const -> decltype((x.segment(i*n, n))) { return x.segment(i*n,n); }
326  inline auto get(KQP_VECTOR(Real) &x, Index i) const -> decltype((x.segment(i*n, n))) { return x.segment(i*n,n); }
327 
328  virtual void mult(const KQP_VECTOR(Real) &x, KQP_VECTOR(Real) &y, bool trans = false) const {
329  // assumes y and x are two different things
330  assert(&x != &y);
331 
332  if (boost::is_complex<Scalar>::value) {
333  // Complex case
334  if (!trans) {
335  // y = G x
336  y.resize(2*n*2*r);
337 
338  for(Index i = 0; i < 2 * r; i++) {
339  // Index of the "real part"
340  int di = i % 2 == 0 ? 0 : -1;
341  Index ip = i + di;
342 
343  get(y, i+2*r) = nu[i] * (get(x,ip) + (nu[i] * (Real)di) * get(x,ip+1));
344  get(y, i) = - get(y, i+2*r) - get(x, 2*r);
345  get(y, i+2*r) -= get(x, 2*r);
346  }
347  } else {
348  // y = G' x
349  y.resize(n * (2*r+1));
350  get(y, 2*r).array().setConstant(0);
351 
352  for(int i = 0; i < 2*r; i++) {
353  // Index of the "real part"
354  int di = i % 2 == 0 ? 0 : -1;
355  Index ip = i + di;
356 
357  get(y,i) = nu[i] * (- get(x,ip) - di * get(x,ip+1) + get(x,ip+r) + di * get(x,ip+1+r) );
358  get(y,i) -= get(x,i) + get(x,i+r);
359  }
360 
361  }
362  } else {
363  // Real case
364  if (!trans) {
365  // y = G x
366  y.resize(2*n*r);
367 
368  for(int i = 0; i < r; i++) {
369  y.segment(i*n, n) = - nu[i] * x.segment(i*n, n) - x.segment(n*r, n);
370  y.segment((i+r)*n, n) = nu[i] * x.segment(i*n, n) - x.segment(n*r, n);
371  }
372  } else {
373  // y = G' x
374  y.resize(n * (r+1));
375  y.segment(n * r, n).array().setConstant(0);
376 
377  for(int i = 0; i < r; i++) {
378  y.segment(i*n, n) = nu[i] * (- x.segment(i*n, n) + x.segment((i+r)*n, n));
379  y.segment(n * r, n) -= x.segment(i*n, n) + x.segment((i+r)*n, n);
380  }
381 
382  }
383  }
384 
385  }
386  virtual Index rows() const { return 2*n*r; }
387  virtual Index cols() const { return n * (r+1); }
388 
389  };
390 
391 
392  template<typename Scalar>
393  void solve_qp(int r,
394  KQP_REAL_OF(Scalar) lambda,
395  const KQP_MATRIX(Scalar) &gramMatrix,
396  const KQP_MATRIX(Scalar) &alpha,
397  const KQP_VECTOR(KQP_REAL_OF(Scalar)) &nu,
398  kqp::cvxopt::ConeQPReturn<KQP_REAL_OF(Scalar)> &result,
399  const cvxopt::ConeQPOptions<KQP_REAL_OF(Scalar)>& options) {
400  typedef typename Eigen::NumTraits<Scalar>::Real Real;
401  const bool isComplex = boost::is_complex<Scalar>::value;
402 
403  KQP_HLOG_DEBUG( "Gram matrix:\n" << gramMatrix);
404  KQP_HLOG_DEBUG( "Alpha:\n" << alpha);
405  KQP_HLOG_DEBUG( "nu:\n" << nu);
406 
407  Index rp = isComplex ? 2 * r : r;
408 
409  Index n = gramMatrix.rows();
410  KQP_VECTOR(Real) c(n*rp + n);
411 
412  for(int i = 0; i < rp; i++)
413  if (isComplex) {
414  Index j = i / 2;
415  if (i % 2 == 0)
416  c.segment(i*n, n) = - (gramMatrix.real() * alpha.col(j).real() + gramMatrix.imag() * alpha.col(j).imag());
417  else
418  c.segment(i*n, n) = - (- gramMatrix.imag() * alpha.col(j).real() + gramMatrix.real() * alpha.col(j).imag());
419  } else
420  c.segment(i*n, n) = - gramMatrix.real() * alpha.real().block(0, i, n, 1);
421 
422  c.segment(rp*n,n).setConstant(lambda / Real(2));
423 
424  QPConstraints<Scalar> G(n, r, nu);
425  KQP_KKTPreSolver<Scalar> kkt_presolver(gramMatrix, nu);
426 
427  KQP_HLOG_DEBUG( "c:\n" << c.adjoint());
428 
429  cvxopt::coneqp<Real>(KMult<Scalar>(n,r, gramMatrix), c, result,
430  false /* No initial value */,
431  cvxopt::Dimensions(),
432  &G, NULL, NULL, NULL,
433  &kkt_presolver,
434  options
435  );
436 
437  }
438 }