Commit 20edde54 authored by Nguyen, Thien Minh's avatar Nguyen, Thien Minh
Browse files

Refactored the ZYZ decomposition inside KAK into a separate circuit



Used it in the optimization routine.
Signed-off-by: Nguyen, Thien Minh's avatarThien Nguyen <nguyentm@ornl.gov>
parent d9e8194d
......@@ -44,6 +44,7 @@ public:
auto aswap = std::make_shared<xacc::circuits::ASWAP>();
auto qfast = std::make_shared<xacc::circuits::QFAST>();
auto kak = std::make_shared<xacc::circuits::KAK>();
auto zyz = std::make_shared<xacc::circuits::ZYZ>();
context.RegisterService<xacc::Instruction>(hwe);
context.RegisterService<xacc::Instruction>(expit);
......@@ -56,6 +57,7 @@ public:
context.RegisterService<xacc::Instruction>(aswap);
context.RegisterService<xacc::Instruction>(qfast);
context.RegisterService<xacc::Instruction>(kak);
context.RegisterService<xacc::Instruction>(zyz);
}
void Stop(BundleContext context) {}
......
......@@ -224,6 +224,167 @@ Eigen::Matrix4cd interactionMatrixExp(double x, double y, double z)
Eigen::MatrixXcd unitary = herm.exp();
return unitary;
}
std::shared_ptr<xacc::CompositeInstruction> singleQubitGateGen(const Eigen::Matrix2cd& in_mat, size_t in_bitIdx)
{
using GateMatrix = Eigen::Matrix2cd;
auto gateRegistry = xacc::getService<xacc::IRProvider>("quantum");
// Use Z-Y decomposition of Nielsen and Chuang (Theorem 4.1).
// An arbitrary one qubit gate matrix can be written as
// U = [ exp(j*(a-b/2-d/2))*cos(c/2), -exp(j*(a-b/2+d/2))*sin(c/2)
// exp(j*(a+b/2-d/2))*sin(c/2), exp(j*(a+b/2+d/2))*cos(c/2)]
// where a,b,c,d are real numbers.
const auto singleQubitGateDecompose = [](const Eigen::Matrix2cd& matrix) -> std::tuple<double, double, double, double> {
if (allClose(matrix, GateMatrix::Identity()))
{
return std::make_tuple(0.0, 0.0, 0.0, 0.0);
}
const auto checkParams = [&matrix](double a, double bHalf, double cHalf, double dHalf) {
GateMatrix U;
U << std::exp(I*(a-bHalf-dHalf))*std::cos(cHalf),
-std::exp(I*(a-bHalf+dHalf))*std::sin(cHalf),
std::exp(I*(a+bHalf-dHalf))*std::sin(cHalf),
std::exp(I*(a+bHalf+dHalf))*std::cos(cHalf);
return allClose(U, matrix);
};
double a, bHalf, cHalf, dHalf;
const double TOLERANCE = 1e-9;
if (std::abs(matrix(0, 1)) < TOLERANCE)
{
auto two_a = fmod(std::arg(matrix(0, 0)*matrix(1, 1)), 2*M_PI);
a = (std::abs(two_a) < TOLERANCE || std::abs(two_a) > 2*M_PI-TOLERANCE) ? 0 : two_a/2.0;
auto dHalf = 0.0;
auto b = std::arg(matrix(1, 1))-std::arg(matrix(0, 0));
std::vector<double> possibleBhalf { fmod(b/2.0, 2 * M_PI), fmod(b/2.0 + M_PI, 2.0 * M_PI) };
std::vector<double> possibleChalf { 0.0, M_PI };
bool found = false;
for (int i = 0; i < possibleBhalf.size(); ++i)
{
for (int j = 0; j < possibleChalf.size(); ++j)
{
bHalf = possibleBhalf[i];
cHalf = possibleChalf[j];
if (checkParams(a, bHalf, cHalf, dHalf))
{
found = true;
break;
}
}
if (found)
{
break;
}
}
assert(found);
}
else if (std::abs(matrix(0, 0)) < TOLERANCE)
{
auto two_a = fmod(std::arg(-matrix(0, 1)*matrix(1, 0)), 2*M_PI);
a = (std::abs(two_a) < TOLERANCE || std::abs(two_a) > 2*M_PI-TOLERANCE) ? 0 : two_a/2.0;
dHalf = 0;
auto b = std::arg(matrix(1, 0))-std::arg(matrix(0, 1)) + M_PI;
std::vector<double> possibleBhalf { fmod(b/2., 2*M_PI), fmod(b/2.+M_PI, 2*M_PI) };
std::vector<double> possibleChalf { M_PI/2., 3./2.*M_PI };
bool found = false;
for (int i = 0; i < possibleBhalf.size(); ++i)
{
for (int j = 0; j < possibleChalf.size(); ++j)
{
bHalf = possibleBhalf[i];
cHalf = possibleChalf[j];
if (checkParams(a, bHalf, cHalf, dHalf))
{
found = true;
break;
}
}
if (found)
{
break;
}
}
assert(found);
}
else
{
auto two_a = fmod(std::arg(matrix(0, 0)*matrix(1, 1)), 2*M_PI);
a = (std::abs(two_a) < TOLERANCE || std::abs(two_a) > 2*M_PI-TOLERANCE) ? 0 : two_a/2.0;
auto two_d = 2.*std::arg(matrix(0, 1))-2.*std::arg(matrix(0, 0));
std::vector<double> possibleDhalf { fmod(two_d/4., 2*M_PI),
fmod(two_d/4.+M_PI/2., 2*M_PI),
fmod(two_d/4.+M_PI, 2*M_PI),
fmod(two_d/4.+3./2.*M_PI, 2*M_PI) };
auto two_b = 2.*std::arg(matrix(1, 0))-2.*std::arg(matrix(0, 0));
std::vector<double> possibleBhalf { fmod(two_b/4., 2*M_PI),
fmod(two_b/4.+M_PI/2., 2*M_PI),
fmod(two_b/4.+M_PI, 2*M_PI),
fmod(two_b/4.+3./2.*M_PI, 2*M_PI) };
auto tmp = std::acos(std::abs(matrix(1, 1)));
std::vector<double> possibleChalf { fmod(tmp, 2*M_PI),
fmod(tmp+M_PI, 2*M_PI),
fmod(-1.*tmp, 2*M_PI),
fmod(-1.*tmp+M_PI, 2*M_PI) };
bool found = false;
for (int i = 0; i < possibleBhalf.size(); ++i)
{
for (int j = 0; j < possibleChalf.size(); ++j)
{
for (int k = 0; k < possibleDhalf.size(); ++k)
{
bHalf = possibleBhalf[i];
cHalf = possibleChalf[j];
dHalf = possibleDhalf[k];
if (checkParams(a, bHalf, cHalf, dHalf))
{
found = true;
break;
}
}
if (found)
{
break;
}
}
if (found)
{
break;
}
}
assert(found);
}
// Final check:
assert(checkParams(a, bHalf, cHalf, dHalf));
return std::make_tuple(a, bHalf, cHalf, dHalf);
};
// Use Z-Y decomposition of Nielsen and Chuang (Theorem 4.1).
// An arbitrary one qubit gate matrix can be writen as
// U = [ exp(j*(a-b/2-d/2))*cos(c/2), -exp(j*(a-b/2+d/2))*sin(c/2)
// exp(j*(a+b/2-d/2))*sin(c/2), exp(j*(a+b/2+d/2))*cos(c/2)]
// where a,b,c,d are real numbers.
// Then U = exp(j*a) Rz(b) Ry(c) Rz(d).
auto [a, bHalf, cHalf, dHalf] = singleQubitGateDecompose(in_mat);
auto composite = gateRegistry->createComposite("__TEMP__COMPOSITE__" + std::to_string(getTempId()));
composite->addInstruction(gateRegistry->createInstruction("Rz", { in_bitIdx }, { 2 * dHalf }));
composite->addInstruction(gateRegistry->createInstruction("Ry", { in_bitIdx }, { 2 * cHalf }));
composite->addInstruction(gateRegistry->createInstruction("Rz", { in_bitIdx }, { 2 * bHalf }));
// Validate U = exp(j*a) Rz(b) Ry(c) Rz(d).
const auto validate = [](const GateMatrix& in_mat, double a, double b, double c, double d) {
GateMatrix Rz_b, Ry_c, Rz_d;
Rz_b << std::exp(-I*b/2.0), 0, 0, std::exp(I*b/2.0);
Rz_d << std::exp(-I*d/2.0), 0, 0, std::exp(I*d/2.0);
Ry_c << std::cos(c/2), -std::sin(c/2), std::sin(c/2), std::cos(c/2);
auto mat = std::exp(I*a)*Rz_b*Ry_c*Rz_d;
return allClose(in_mat, mat);
};
assert(validate(in_mat, a, 2*bHalf, 2*cHalf, 2*dHalf));
return composite;
}
}
using namespace xacc;
......@@ -364,164 +525,6 @@ Eigen::MatrixXcd KAK::KakDecomposition::toMat() const
std::shared_ptr<CompositeInstruction> KAK::KakDecomposition::toGates(size_t in_bit1, size_t in_bit2) const
{
auto gateRegistry = xacc::getService<IRProvider>("quantum");
// Use Z-Y decomposition of Nielsen and Chuang (Theorem 4.1).
// An arbitrary one qubit gate matrix can be written as
// U = [ exp(j*(a-b/2-d/2))*cos(c/2), -exp(j*(a-b/2+d/2))*sin(c/2)
// exp(j*(a+b/2-d/2))*sin(c/2), exp(j*(a+b/2+d/2))*cos(c/2)]
// where a,b,c,d are real numbers.
const auto singleQubitGateDecompose = [](const GateMatrix& matrix) -> std::tuple<double, double, double, double> {
if (allClose(matrix, GateMatrix::Identity()))
{
return std::make_tuple(0.0, 0.0, 0.0, 0.0);
}
const auto checkParams = [&matrix](double a, double bHalf, double cHalf, double dHalf) {
GateMatrix U;
U << std::exp(I*(a-bHalf-dHalf))*std::cos(cHalf),
-std::exp(I*(a-bHalf+dHalf))*std::sin(cHalf),
std::exp(I*(a+bHalf-dHalf))*std::sin(cHalf),
std::exp(I*(a+bHalf+dHalf))*std::cos(cHalf);
return allClose(U, matrix);
};
double a, bHalf, cHalf, dHalf;
const double TOLERANCE = 1e-9;
if (std::abs(matrix(0, 1)) < TOLERANCE)
{
auto two_a = fmod(std::arg(matrix(0, 0)*matrix(1, 1)), 2*M_PI);
a = (std::abs(two_a) < TOLERANCE || std::abs(two_a) > 2*M_PI-TOLERANCE) ? 0 : two_a/2.0;
auto dHalf = 0.0;
auto b = std::arg(matrix(1, 1))-std::arg(matrix(0, 0));
std::vector<double> possibleBhalf { fmod(b/2.0, 2 * M_PI), fmod(b/2.0 + M_PI, 2.0 * M_PI) };
std::vector<double> possibleChalf { 0.0, M_PI };
bool found = false;
for (int i = 0; i < possibleBhalf.size(); ++i)
{
for (int j = 0; j < possibleChalf.size(); ++j)
{
bHalf = possibleBhalf[i];
cHalf = possibleChalf[j];
if (checkParams(a, bHalf, cHalf, dHalf))
{
found = true;
break;
}
}
if (found)
{
break;
}
}
assert(found);
}
else if (std::abs(matrix(0, 0)) < TOLERANCE)
{
auto two_a = fmod(std::arg(-matrix(0, 1)*matrix(1, 0)), 2*M_PI);
a = (std::abs(two_a) < TOLERANCE || std::abs(two_a) > 2*M_PI-TOLERANCE) ? 0 : two_a/2.0;
dHalf = 0;
auto b = std::arg(matrix(1, 0))-std::arg(matrix(0, 1)) + M_PI;
std::vector<double> possibleBhalf { fmod(b/2., 2*M_PI), fmod(b/2.+M_PI, 2*M_PI) };
std::vector<double> possibleChalf { M_PI/2., 3./2.*M_PI };
bool found = false;
for (int i = 0; i < possibleBhalf.size(); ++i)
{
for (int j = 0; j < possibleChalf.size(); ++j)
{
bHalf = possibleBhalf[i];
cHalf = possibleChalf[j];
if (checkParams(a, bHalf, cHalf, dHalf))
{
found = true;
break;
}
}
if (found)
{
break;
}
}
assert(found);
}
else
{
auto two_a = fmod(std::arg(matrix(0, 0)*matrix(1, 1)), 2*M_PI);
a = (std::abs(two_a) < TOLERANCE || std::abs(two_a) > 2*M_PI-TOLERANCE) ? 0 : two_a/2.0;
auto two_d = 2.*std::arg(matrix(0, 1))-2.*std::arg(matrix(0, 0));
std::vector<double> possibleDhalf { fmod(two_d/4., 2*M_PI),
fmod(two_d/4.+M_PI/2., 2*M_PI),
fmod(two_d/4.+M_PI, 2*M_PI),
fmod(two_d/4.+3./2.*M_PI, 2*M_PI) };
auto two_b = 2.*std::arg(matrix(1, 0))-2.*std::arg(matrix(0, 0));
std::vector<double> possibleBhalf { fmod(two_b/4., 2*M_PI),
fmod(two_b/4.+M_PI/2., 2*M_PI),
fmod(two_b/4.+M_PI, 2*M_PI),
fmod(two_b/4.+3./2.*M_PI, 2*M_PI) };
auto tmp = std::acos(std::abs(matrix(1, 1)));
std::vector<double> possibleChalf { fmod(tmp, 2*M_PI),
fmod(tmp+M_PI, 2*M_PI),
fmod(-1.*tmp, 2*M_PI),
fmod(-1.*tmp+M_PI, 2*M_PI) };
bool found = false;
for (int i = 0; i < possibleBhalf.size(); ++i)
{
for (int j = 0; j < possibleChalf.size(); ++j)
{
for (int k = 0; k < possibleDhalf.size(); ++k)
{
bHalf = possibleBhalf[i];
cHalf = possibleChalf[j];
dHalf = possibleDhalf[k];
if (checkParams(a, bHalf, cHalf, dHalf))
{
found = true;
break;
}
}
if (found)
{
break;
}
}
if (found)
{
break;
}
}
assert(found);
}
// Final check:
assert(checkParams(a, bHalf, cHalf, dHalf));
return std::make_tuple(a, bHalf, cHalf, dHalf);
};
const auto singleQubitGateGen = [&](const GateMatrix& in_mat, size_t in_bitIdx) {
// Use Z-Y decomposition of Nielsen and Chuang (Theorem 4.1).
// An arbitrary one qubit gate matrix can be writen as
// U = [ exp(j*(a-b/2-d/2))*cos(c/2), -exp(j*(a-b/2+d/2))*sin(c/2)
// exp(j*(a+b/2-d/2))*sin(c/2), exp(j*(a+b/2+d/2))*cos(c/2)]
// where a,b,c,d are real numbers.
// Then U = exp(j*a) Rz(b) Ry(c) Rz(d).
auto [a, bHalf, cHalf, dHalf] = singleQubitGateDecompose(in_mat);
auto composite = gateRegistry->createComposite("__TEMP__COMPOSITE__" + std::to_string(getTempId()));
composite->addInstruction(gateRegistry->createInstruction("Rz", { in_bitIdx }, { 2 * dHalf }));
composite->addInstruction(gateRegistry->createInstruction("Ry", { in_bitIdx }, { 2 * cHalf }));
composite->addInstruction(gateRegistry->createInstruction("Rz", { in_bitIdx }, { 2 * bHalf }));
// Validate U = exp(j*a) Rz(b) Ry(c) Rz(d).
const auto validate = [](const GateMatrix& in_mat, double a, double b, double c, double d) {
GateMatrix Rz_b, Ry_c, Rz_d;
Rz_b << std::exp(-I*b/2.0), 0, 0, std::exp(I*b/2.0);
Rz_d << std::exp(-I*d/2.0), 0, 0, std::exp(I*d/2.0);
Ry_c << std::cos(c/2), -std::sin(c/2), std::sin(c/2), std::cos(c/2);
auto mat = std::exp(I*a)*Rz_b*Ry_c*Rz_d;
return allClose(in_mat, mat);
};
assert(validate(in_mat, a, 2*bHalf, 2*cHalf, 2*dHalf));
return composite;
};
const auto generateInteractionComposite = [&](size_t bit1, size_t bit2, double x, double y, double z) {
const double xAngle = M_PI * (x * -2 / M_PI + 0.5);
const double yAngle = M_PI * (y * -2 / M_PI + 0.5);
......@@ -976,5 +979,23 @@ KAK::KakDecomposition KAK::canonicalizeInteraction(double x, double y, double z)
assert(allClose(result.toMat(), interactionMatrixExp(x, y, z)));
return result;
}
bool ZYZ::expand(const xacc::HeterogeneousMap& runtimeOptions)
{
Eigen::Matrix2cd unitary;
if (runtimeOptions.keyExists<Eigen::Matrix2cd>("unitary"))
{
unitary = runtimeOptions.get<Eigen::Matrix2cd>("unitary");
}
else
{
xacc::error("unitary matrix is required.");
return false;
}
assert(isUnitary(unitary));
auto decomposed = singleQubitGateGen(unitary, 0);
addInstructions(decomposed->getInstructions());
return true;
}
} // namespace circuits
} // namespace xacc
\ No newline at end of file
......@@ -6,6 +6,18 @@
namespace xacc {
namespace circuits {
// Use Z-Y-Z decomposition of Nielsen and Chuang (Theorem 4.1).
// An arbitrary one qubit gate matrix can be writen as
// U = exp(j*a) Rz(b) Ry(c) Rz(d).
class ZYZ : public xacc::quantum::Circuit
{
public:
ZYZ() : Circuit("z-y-z") {}
bool expand(const xacc::HeterogeneousMap &runtimeOptions) override;
const std::vector<std::string> requiredKeys() override { return { "unitary" }; }
DEFINE_CLONE(ZYZ);
};
// KAK decomposition via *Magic* Bell basis transformation
// Reference:
// https://arxiv.org/pdf/quant-ph/0211002.pdf
......
#include "GateMergeOptimizer.hpp"
#include <cassert>
#include "GateMergeOptimizer.hpp"
#include "GateFusion.hpp"
#include "xacc_service.hpp"
namespace xacc {
namespace quantum {
void MergeSingleQubitGatesOptimizer::apply(std::shared_ptr<CompositeInstruction> program, const std::shared_ptr<Accelerator> accelerator, const HeterogeneousMap &options)
{
auto gateRegistry = xacc::getService<xacc::IRProvider>("quantum");
std::set<size_t> processInstIdx;
for (size_t instIdx = 0; instIdx < program->nInstructions(); ++instIdx)
{
......@@ -16,12 +19,28 @@ void MergeSingleQubitGatesOptimizer::apply(std::shared_ptr<CompositeInstruction>
assert(!xacc::container::contains(processInstIdx, instIdx));
processInstIdx.emplace(instIdx);
}
auto tmpKernel = gateRegistry->createComposite("__TMP__");
for (const auto& instIdx: sequence)
{
auto instrPtr = program->getInstruction(instIdx)->clone();
assert(instrPtr->bits().size() == 1);
// Remap to bit 0 for fusing
instrPtr->setBits({ 0 });
tmpKernel->addInstruction(instrPtr);
}
auto fuser = xacc::getService<xacc::quantum::GateFuser>("default");
fuser->initialize(tmpKernel);
const Eigen::Matrix2cd uMat = fuser->calcFusedGate(1);
std::cout << "Unitary matrix: \n" << uMat << "\n";
auto zyz = std::dynamic_pointer_cast<quantum::Circuit>(xacc::getService<Instruction>("z-y-z"));
const bool expandOk = zyz->expand({
std::make_pair("unitary", uMat)
});
// std::cout << "Found sequence:\n";
// for (const auto& idx : sequence)
// {
// std::cout << program->getInstruction(idx)->toString() << "\n";
// }
assert(zyz->nInstructions() == 3);
std::cout << "Decomposed: \n" << zyz->toString() << "\n";
// TODO: further optimize this sequence: remove trivial rotation (0, 2*pi, etc.)
}
}
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment