qda_linear_solver.cpp

#include <EigenUnsupported/Eigen/KroneckerProduct>
#include <EigenUnsupported/Eigen/MatrixFunctions>
#include "Core/QuantumMachine/OriginQuantumMachine.h"
#include "Core/QuantumCircuit/QGate.h"
#include "Core/QuantumCircuit/QCircuit.h"
#include "Core/QuantumCircuit/QProgram.h"
#include "Core/Utilities/QProgInfo/Visualization/QVisualization.h"
#include "Core/Utilities/Tools/MatrixDecomposition.h"
#include "Core/Utilities/UnitaryDecomposer/QSDecomposition.h"
#include "QAlg/Base_QCircuit/AmplitudeEncode.h"
#include "block_encoding.h"
#include "qda_linear_solver.h"
#include "utils.h"

static VectorXcd v0(2);
static MatrixXcd PauliX(2, 2);
static MatrixXcd PauliY(2, 2);
static MatrixXcd PauliP(2, 2);
static MatrixXcd PauliM(2, 2);
static bool is_init_consts = false;

inline void init_consts() {
  if (is_init_consts) return;
  v0 << 1, 0;
  PauliX << 0, 1,
            1, 0;
  PauliY << 0, dcomplex(0, -1),
            dcomplex(0, 1), 0;
  PauliP = (PauliX + dcomplex(0, 1) * PauliY) / 2;
  PauliM = (PauliX - dcomplex(0, 1) * PauliY) / 2;
  is_init_consts = true;
}

MatrixXcd EF_R_l(MatrixXcd H) {
  const int lmbd = 16;    // 阶数
  const double delta = 0.1;
  const double delta2 = delta * delta;
  const int sign = 1; // := pow(-1, lmbd)

  auto T_l = [=](double x) -> double {
    if (-1 <= x && x <= 1) return cos(lmbd * acos(x));
    else if (x > 1)        return cosh(lmbd * acosh(x));
    else if (x < -1)       return sign * cosh(lmbd * acosh(-x));
  };
  auto t = [=](double x) -> double {
    return (x * x - delta2) / (1 - delta2);
  };
  auto t0 = t(0);
  auto R_l = [=](double x) -> double {
    auto p = T_l(-1 + 2 * t(x));
    auto q = T_l(-1 + 2 * t0);
    return p / q;
  };

  ComplexEigenSolver<MatrixXcd> solver;
  solver.compute(H);
  if (solver.info() != Eigen::Success) abort();
  auto V = solver.eigenvectors();
  auto D = solver.eigenvalues();
  VectorXcd R_l_D(D.size());
  for (int i = 0; i < D.size(); i++)
    R_l_D(i) = R_l(D(i).real());
  return V * R_l_D.asDiagonal() * V.inverse();
}

inline MatrixXcd exp_iHt_approx(MatrixXcd H, float t=1.0, int order=1) {
  // https://en.wikipedia.org/wiki/Matrix_exponential
  // exp(-iHt) ~= I - iHt + (-iHt)^2 / 2 + ... + (-iHt)^k / k!
  int N = H.rows();
  MatrixXcd A = MatrixXcd::Identity(N, N);
  MatrixXcd _iHt = dcomplex(0, -1) * H * t;
  int fac = 1;
  for (int o = 1; o <= order; o++) {
    fac *= o;
    A += _iHt.pow(o) / fac;
  }
  return A;
}

inline QCircuit matrix_decompose(DecompositionMode method, MatrixXcd mat, QVec qv) {
  if (false) {
    auto err = (mat * mat.adjoint() - MatrixXcd::Identity(mat.rows(), mat.cols())).cwiseAbs().mean();
    if (err > 1e-10) {
      cout << mat << endl;
      cout << "[matrix_decompose] err: " << err << endl;
    }
  }

  QCircuit qcir;
  QMatrixXcd qmat = mat;
  switch (method) {
    case DecompositionMode::QR:
      qcir = matrix_decompose_qr(qv, qmat, false);  // NOTE: must keep false here
      break;
    case DecompositionMode::CSD:
      qcir = unitary_decomposer_nq(mat, qv, DecompositionMode::CSD, true);   // NOTE: must keep true here
      break;
    case DecompositionMode::QSD:
      qcir = unitary_decomposer_nq(mat, qv, DecompositionMode::QSD, true);   // NOTE: must keep true here
      break;
    default:
      throw invalid_argument("invalid decompose_method");
  }
  return qcir;
}


// The ideal yet naive QDA implementation
VectorXcd linear_solver_ideal(MatrixXcd A, VectorXcd b, DecompositionMode decompose_method=DecompositionMode::QR) {
  init_consts();

  // time-dependent hamiltonian
  int N = A.rows();
  MatrixXcd Qb = MatrixXcd::Identity(N, N) - b * b.adjoint();
  MatrixXcd H0 = kroneckerProduct(PauliX, Qb);
  MatrixXcd H1 = kroneckerProduct(PauliP, A * Qb) + kroneckerProduct(PauliM, Qb * A);
  auto H_s = [&](float s) -> MatrixXcd { return (1 - s) * H0 + s * H1; };
  VectorXd init_state = kroneckerProduct(v0, b).real();  // |0,b>

  // gated quantum computing
  const int S = 200;    // premise 1: adequate stage count
  const int T = 4;      // premise 2: adequate physical time during each stage
  const size_t n_qubit = ceil(log2(H0.rows()));
  CPUQVM qvm;
  qvm.setConfigure({n_qubit, n_qubit});
  qvm.init();
  QVec qv = qvm.qAllocMany(n_qubit);
  QCircuit qcir;
  // init state
  vector<double> amplitude(init_state.data(), init_state.data() + init_state.size());
  qcir << amplitude_encode(qv, amplitude);  // premise 3: correct initial state preparation
  // adiabatic evolution
  for (int s = 1; s <= S; s++) {
    MatrixXcd H = H_s(float(s) / S);  // premise 4: ideal time-independent hamiltonian schedule
    MatrixXcd iHt = dcomplex(0, -1) * H * T;
    MatrixXcd U_iHt = iHt.exp();      // premise 5: ideal time evolution operator
    QCircuit qc_TE = matrix_decompose(decompose_method, U_iHt, qv);
    qcir << qc_TE << BARRIER(qv);
  }
  // final state
  QProg qprog = createEmptyQProg() << qcir;
  qvm.directlyRun(qprog);
  QStat qs = qvm.getQState();

  // result
  qs = QStat(qs.begin(), qs.begin() + N);  // project only the working qubit
  VectorXcd state = Map<VectorXcd>(qs.data(), qs.size());
  return state /= state.norm();
}


// The basic implementation strictly follows all contest-specified restrictions :(
VectorXcd linear_solver_contest(MatrixXcd A, VectorXcd b, DecompositionMode decompose_method=DecompositionMode::QR) {
  init_consts();

  // time-dependent hamiltonian
  int N = A.rows();
  MatrixXcd Qb = MatrixXcd::Identity(N, N) - b * b.adjoint();
  MatrixXcd H0 = kroneckerProduct(PauliX, Qb);
  MatrixXcd H1 = kroneckerProduct(PauliP, A * Qb) + kroneckerProduct(PauliM, Qb * A);
  VectorXcd init_state = kroneckerProduct(v0, b);
  auto H_s = [&](float s) -> MatrixXcd { return (1 - s) * H0 + s * H1; };

  // gated quantum computing
  const int S = 200;    // restrict 1: fixed as contest required
  const int T = 1;      // restrict 2: fixed as contest required (?)
  const size_t n_qubit = ceil(log2(H0.rows()));   // since we're going to block-encode the hamiltonian H(s), not the matrix A in equation
  const size_t n_ancilla = 1;   // NOTE: modify this according to your block_encode method :)
  const size_t n_qubit_ex = n_ancilla + n_qubit;
  CPUQVM qvm;
  qvm.setConfigure({n_qubit_ex, n_qubit_ex});
  qvm.init();
  QVec qv = qvm.qAllocMany(n_qubit_ex);
  QCircuit qcir;
  // init state
  VectorXcd anc0 = v0;
  for (int i = 1; i < n_ancilla; i++) anc0 = kroneckerProduct(anc0, v0).eval();   // NOTE: avoid aliasing effect
  VectorXd init_state_ex = kroneckerProduct(anc0, init_state).real();
  vector<double> amplitude(init_state_ex.data(), init_state_ex.data() + init_state_ex.size());
  qcir << amplitude_encode(qv, amplitude);   // |anc_BE,0,b>
  // adiabatic evolution
  for (int s = 1; s <= S; s++) {
    MatrixXcd H = H_s(float(s) / S);
    MatrixXcd iHt = exp_iHt_approx(H, T);   // restrict 3: approx as contest required
    iHt = normalize_QSVT(iHt);
    MatrixXcd U_iHt = block_encoding_QSVT(iHt).unitary;  // restrict 4: block_encode as contest required, and as consequence of approx
    QCircuit qc_TE = matrix_decompose(decompose_method, U_iHt, qv);
    qcir << qc_TE << BARRIER(qv);
  }
  // final state
  QProg qprog = createEmptyQProg() << qcir;
  qvm.directlyRun(qprog);
  QStat qs = qvm.getQState();

  // result
  qs = QStat(qs.begin(), qs.begin() + N);  // project only the working qubit
  VectorXcd state = Map<VectorXcd>(qs.data(), qs.size());
  return state /= state.norm();
}


// The optimized implementation with various tricks for the best precision :)
VectorXcd linear_solver_ours(MatrixXcd A, VectorXcd b, DecompositionMode decompose_method=DecompositionMode::QR) {
  init_consts();

  // time-dependent hamiltonian
  int N = A.rows();
  MatrixXcd Qb = MatrixXcd::Identity(N, N) - b * b.adjoint();
  MatrixXcd H0 = kroneckerProduct(PauliX, Qb);
  MatrixXcd H1 = kroneckerProduct(PauliP, A * Qb) + kroneckerProduct(PauliM, Qb * A);
  VectorXcd init_state = kroneckerProduct(v0, b);
  // AQC(P) schedule
  float k = 5.82842712474619;   // condition_number of the original A
  float p = 2.0;
  auto f_ = [&](float s) -> float {
    float t = 1 + s * (pow(k, p - 1) - 1);
    return k / (k - 1) * (1 - pow(t, 1 / (1 - p)));
  };
  auto H_s = [&](float s) -> MatrixXcd {
    float f_s = f_(s);
    return (1 - f_s) * H0 + f_s * H1;
  };

  // gated quantum computing
  const int S = 400;    // many logical steps 
  const int T = 2;      // enough long physical time of each step
  const size_t n_qubit = ceil(log2(H0.rows()));   // since we're going to block-encode the hamiltonian H(s), not the matrix A in equation
  const size_t n_ancilla = 1;   // NOTE: modify this according to your block_encode method :)
  const size_t n_qubit_ex = n_ancilla + n_qubit;
  CPUQVM qvm;
  qvm.setConfigure({n_qubit_ex, n_qubit_ex});
  qvm.init();
  QVec qv = qvm.qAllocMany(n_qubit_ex);
  QCircuit qcir;
  // init state
  VectorXcd anc0 = v0;
  for (int i = 1; i < n_ancilla; i++) anc0 = kroneckerProduct(anc0, v0).eval();   // NOTE: avoid aliasing effect
  VectorXd init_state_ex = kroneckerProduct(anc0, init_state).real();
  vector<double> amplitude(init_state_ex.data(), init_state_ex.data() + init_state_ex.size());
  qcir << amplitude_encode(qv, amplitude);   // |anc_BE,0,b>
  // adiabatic evolution
  for (int s = 1; s <= S; s++) {
    MatrixXcd H = H_s(float(s) / S);
    MatrixXcd iHt = exp_iHt_approx(H, T);
    iHt = normalize_QSVT(iHt);
    MatrixXcd U_iHt = block_encoding_QSVT(iHt).unitary;
    QCircuit qc_TE = matrix_decompose(decompose_method, U_iHt, qv);
    qcir << qc_TE << BARRIER(qv);
  }
  // eigen filter
  MatrixXcd EF = EF_R_l(H1);
  //EF = normalize_QSVT(EF);
  MatrixXcd U_EF = block_encoding_QSVT(EF).unitary;
  // err: 1e-9, 精度不够 QR 分解 =_=
  QCircuit qc_EF = matrix_decompose(DecompositionMode::QSD, U_EF, qv);
  qcir << qc_EF << BARRIER(qv);
  // final state
  QProg qprog = createEmptyQProg() << qcir;
  qvm.directlyRun(qprog);
  QStat qs = qvm.getQState();

  // result
  qs = QStat(qs.begin(), qs.begin() + N);  // project only the working qubit
  VectorXcd state = Map<VectorXcd>(qs.data(), qs.size());
  return state /= state.norm();
}


// ↓↓ keep signature for the contest solution
#define DEBUG
#ifdef DEBUG
#include <iostream>
using namespace std;
#endif

qdals_res qda_linear_solver(MatrixXcd A, VectorXcd b) {
  // classical solution |x_r>
  VectorXcd s_r = A.colPivHouseholderQr().solve(b);
  VectorXcd x_r = s_r / s_r.norm();
  // quantum solution |x>
  auto b_norm = b.norm();
  VectorXcd x = linear_solver_ours(A / b_norm, b / b_norm, DecompositionMode::QR);
  x = -x;     // just fix the phase-flip :)
  // fidelity <x_r|x>
  dcomplex fidelity = x_r.adjoint().dot(x);

#ifdef DEBUG
  cout << "Classical solution: " << x_r.real().transpose() << endl;
  cout << "Quantum solution: " << x.real().transpose() << endl;
  cout << "Fidelity: " << abs(fidelity) << endl;
#endif

  // result pack
  return { x, fidelity };
}