Skip to content

Commit

Permalink
WIP: Adding logarithm of gate method (#117)
Browse files Browse the repository at this point in the history
* Adding logarithm of gate method

* Adding 2x2 exp/log convenience methods

* Fixing imaginary Exp convention

* Anti-controlled Exp options

* Renaming convenience exp2x2/log2x2/mul2x2

* HamiltonianOp control bit toggles

* Basic HamiltonianOp toggle unit test

* Commentary on HamiltonianOp toggles
  • Loading branch information
WrathfulSpatula authored Dec 26, 2018
1 parent 3d26395 commit fadccd5
Show file tree
Hide file tree
Showing 7 changed files with 232 additions and 90 deletions.
8 changes: 7 additions & 1 deletion include/common/qrack_types.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,8 @@

#pragma once

#include <memory>

#define bitLenInt uint8_t
#define bitCapInt uint64_t
#define bitsInByte 8
Expand Down Expand Up @@ -39,4 +41,8 @@

namespace Qrack {
typedef std::shared_ptr<complex> BitOp;
}

void mul2x2(complex* left, complex* right, complex* out);
void exp2x2(complex* matrix2x2, complex* outMatrix2x2);
void log2x2(complex* matrix2x2, complex* outMatrix2x2);
} // namespace Qrack
32 changes: 30 additions & 2 deletions include/hamiltonian.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -27,37 +27,65 @@ struct HamiltonianOp {
BitOp matrix;
bitLenInt* controls;
bitLenInt controlLen;
bool anti;
bool* toggles;

HamiltonianOp(bitLenInt target, BitOp mtrx)
: targetBit(target)
, matrix(mtrx)
, controls(NULL)
, controlLen(0)
, anti(false)
, toggles(NULL)
{
}

HamiltonianOp(bitLenInt* ctrls, bitLenInt ctrlLen, bitLenInt target, BitOp mtrx)
HamiltonianOp(bitLenInt* ctrls, bitLenInt ctrlLen, bitLenInt target, BitOp mtrx, bool antiCtrled = false,
bool* ctrlToggles = NULL)
: targetBit(target)
, matrix(mtrx)
, controls(new bitLenInt[ctrlLen])
, controlLen(ctrlLen)
, anti(antiCtrled)
, toggles(NULL)
{
std::copy(ctrls, ctrls + ctrlLen, controls);

if (ctrlToggles) {
toggles = new bool[ctrlLen];
std::copy(ctrlToggles, ctrlToggles + ctrlLen, toggles);
}
}

~HamiltonianOp()
{
if (controls) {
delete[] controls;
}

if (toggles) {
delete[] toggles;
}
}
};

/**
* To define a Hamiltonian, give a vector of controlled single bit gates ("HamiltonianOp" instances) that are
* applied by left-multiplication in low-to-high vector index order on the state vector.
*
* \warning Hamiltonian components might not commute, and observe the component factor of 2 * pi.
* To specify Hamiltonians with interaction terms, arbitrary sets of control bits may be specified for each term, in
* which case the term is acted only if the (superposable) control bits are true. The "antiCtrled" bool flips the
* overall convention, such that the term is acted only if all control bits are false. Additionally, for a combination
* of control bits and "anti-control" bits, an array of booleans, "ctrlToggles," of length "ctrlLen" may be specified
* that flips the activation state for each control bit, (relative the global anti- on/off convention,) without altering
* the state of the control bits.
*
* The point of this "toggle" behavior is to allow enumeration of arbitrary local Hamiltonian terms with permutations of
* a set of control bits. For example, a Hamiltonian might represent an array of local electromagnetic potential wells.
* If there are 4 wells, each with independent potentials, control "toggles" could be used on two control bits, to
* enumerate all four permutations of two control bits with four different local Hamiltonian terms.
*
* \warning Hamiltonian components might not commute.
*
* As a general point of linear algebra, where A and B are linear operators, e^{i * (A + B) * t} = e^{i * A * t} *
* e^{i * B * t} might NOT hold, if the operators A and B do not commute. As a rule of thumb, A will commute with B
Expand Down
15 changes: 12 additions & 3 deletions include/qinterface.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -634,11 +634,20 @@ class QInterface {
virtual void Exp(real1 radians, bitLenInt qubitIndex);

/**
* Exponentiation of arbitrary 2x2 gate
* Imaginary exponentiation of arbitrary 2x2 gate
*
* Applies \f$ e^{-i*Op*I} \f$, where "Op" is a 2x2 matrix, (with controls on the application of the gate).
* Applies \f$ e^{-i*Op} \f$, where "Op" is a 2x2 matrix, (with controls on the application of the gate).
*/
virtual void Exp(bitLenInt* controls, bitLenInt controlLen, bitLenInt qubit, complex* matrix2x2);
virtual void Exp(
bitLenInt* controls, bitLenInt controlLen, bitLenInt qubit, complex* matrix2x2, bool antiCtrled = false);

/**
* Logarithm of arbitrary 2x2 gate
*
* Applies \f$ log(Op) \f$, where "Op" is a 2x2 matrix, (with controls on the application of the gate).
*/
virtual void Log(
bitLenInt* controls, bitLenInt controlLen, bitLenInt qubit, complex* matrix2x2, bool antiCtrled = false);

/**
* Dyadic fraction (identity) exponentiation gate
Expand Down
18 changes: 17 additions & 1 deletion src/qinterface/operators.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -171,7 +171,23 @@ void QInterface::TimeEvolve(Hamiltonian h, real1 timeDiff)
for (int j = 0; j < 4; j++) {
mtrx[j] = opMtrx[j] * (-timeDiff);
}
Exp(op->controls, op->controlLen, op->targetBit, mtrx);
if (op->toggles) {
for (bitLenInt j = 0; j < op->controlLen; j++) {
if (op->toggles[j]) {
X(op->controls[j]);
}
}
}

Exp(op->controls, op->controlLen, op->targetBit, mtrx, op->anti);

if (op->toggles) {
for (bitLenInt j = 0; j < op->controlLen; j++) {
if (op->toggles[j]) {
X(op->controls[j]);
}
}
}
}
}
} // namespace Qrack
103 changes: 101 additions & 2 deletions src/qinterface/protected.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -10,9 +10,9 @@
// See LICENSE.md in the project root or https://www.gnu.org/licenses/lgpl-3.0.en.html
// for details.

#include "common/qrack_types.hpp"
#include <algorithm>

#include "qinterface.hpp"
#include <cmath>

namespace Qrack {

Expand All @@ -36,4 +36,103 @@ void rotate(BidirectionalIterator first, BidirectionalIterator middle, Bidirecti

template void rotate<complex*>(complex* first, complex* middle, complex* last, bitCapInt stride);

void mul2x2(complex* left, complex* right, complex* out)
{
out[0] = (left[0] * right[0]) + (left[1] * right[2]);
out[1] = (left[0] * right[1]) + (left[1] * right[3]);
out[2] = (left[2] * right[0]) + (left[3] * right[2]);
out[3] = (left[2] * right[1]) + (left[3] * right[3]);
}

void _expLog2x2(complex* matrix2x2, complex* outMatrix2x2, bool isExp)
{
// Solve for the eigenvalues and eigenvectors of a 2x2 matrix, diagonalize, exponentiate, return to the original
// basis, and apply.

// Diagonal matrices are a special case.
bool isDiag = true;
if (norm(matrix2x2[1]) > min_norm) {
isDiag = false;
} else if (norm(matrix2x2[2]) > min_norm) {
isDiag = false;
}

complex expOfGate[4];
complex jacobian[4];
complex inverseJacobian[4];
complex tempMatrix2x2[4];

// Diagonalize the matrix, if it is not already diagonal. Otherwise, copy it into the temporary matrix.
if (!isDiag) {
complex trace = matrix2x2[0] + matrix2x2[3];
complex determinant = (matrix2x2[0] * matrix2x2[3]) - (matrix2x2[1] * matrix2x2[2]);
complex quadraticRoot =
sqrt((matrix2x2[0] - matrix2x2[3]) * (matrix2x2[0] - matrix2x2[3]) - (real1)(4.0) * determinant);
complex eigenvalue1 = (trace + quadraticRoot) / (real1)2.0;
complex eigenvalue2 = (trace - quadraticRoot) / (real1)2.0;

if (norm(matrix2x2[1]) > min_norm) {
jacobian[0] = matrix2x2[1];
jacobian[2] = eigenvalue1 - matrix2x2[0];

jacobian[1] = matrix2x2[1];
jacobian[3] = eigenvalue2 - matrix2x2[0];
} else {
jacobian[0] = eigenvalue1 - matrix2x2[3];
jacobian[2] = matrix2x2[2];

jacobian[1] = eigenvalue2 - matrix2x2[3];
jacobian[3] = matrix2x2[2];
}

real1 nrm = std::sqrt(norm(jacobian[0]) + norm(jacobian[2]));
jacobian[0] /= nrm;
jacobian[2] /= nrm;

nrm = std::sqrt(norm(jacobian[1]) + norm(jacobian[3]));
jacobian[1] /= nrm;
jacobian[3] /= nrm;

determinant = (jacobian[0] * jacobian[3]) - (jacobian[1] * jacobian[2]);
inverseJacobian[0] = jacobian[3] / determinant;
inverseJacobian[1] = -jacobian[1] / determinant;
inverseJacobian[2] = -jacobian[2] / determinant;
inverseJacobian[3] = jacobian[0] / determinant;

mul2x2(matrix2x2, jacobian, tempMatrix2x2);
mul2x2(inverseJacobian, tempMatrix2x2, expOfGate);
} else {
std::copy(matrix2x2, matrix2x2 + 4, expOfGate);
}

if (isExp) {
// In this branch, we calculate e^(matrix2x2).

// Note: For a (2x2) hermitian input gate, this theoretically produces a unitary output transformation.
expOfGate[0] = ((real1)std::exp(real(expOfGate[0]))) *
complex((real1)cos(imag(expOfGate[0])), (real1)sin(imag(expOfGate[0])));
expOfGate[1] = complex(ZERO_R1, ZERO_R1);
expOfGate[2] = complex(ZERO_R1, ZERO_R1);
expOfGate[3] = ((real1)std::exp(real(expOfGate[3]))) *
complex((real1)cos(imag(expOfGate[3])), (real1)sin(imag(expOfGate[3])));
} else {
// In this branch, we calculate log(matrix2x2).
expOfGate[0] = complex(std::log(abs(expOfGate[0])), arg(expOfGate[0]));
expOfGate[1] = complex(ZERO_R1, ZERO_R1);
expOfGate[2] = complex(ZERO_R1, ZERO_R1);
expOfGate[3] = complex(std::log(abs(expOfGate[3])), arg(expOfGate[3]));
}

if (!isDiag) {
mul2x2(expOfGate, inverseJacobian, tempMatrix2x2);
mul2x2(jacobian, tempMatrix2x2, expOfGate);
}

std::copy(expOfGate, expOfGate + 4, outMatrix2x2);
}

void exp2x2(complex* matrix2x2, complex* outMatrix2x2) { _expLog2x2(matrix2x2, outMatrix2x2, true); }

void log2x2(complex* matrix2x2, complex* outMatrix2x2) { _expLog2x2(matrix2x2, outMatrix2x2, false); }

} // namespace Qrack
98 changes: 20 additions & 78 deletions src/qinterface/rotational.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -64,90 +64,32 @@ void QInterface::Exp(real1 radians, bitLenInt qubit)
ApplySingleBit(expIdentity, true, qubit);
}

void matrix2x2Mul(complex* left, complex* right, complex* out)
/// Imaginary exponentiate of arbitrary single bit gate
void QInterface::Exp(bitLenInt* controls, bitLenInt controlLen, bitLenInt qubit, complex* matrix2x2, bool antiCtrled)
{
out[0] = (left[0] * right[0]) + (left[1] * right[2]);
out[1] = (left[0] * right[1]) + (left[1] * right[3]);
out[2] = (left[2] * right[0]) + (left[3] * right[2]);
out[3] = (left[2] * right[1]) + (left[3] * right[3]);
}

/// Exponentiate of arbitrary single bit gate
void QInterface::Exp(bitLenInt* controls, bitLenInt controlLen, bitLenInt qubit, complex* matrix2x2)
{
// Solve for the eigenvalues and eigenvectors of a 2x2 matrix, diagonalize, exponentiate, return to the original
// basis, and apply.

// Diagonal matrices are a special case.
bool isDiag = true;
if (norm(matrix2x2[1]) > min_norm) {
isDiag = false;
} else if (norm(matrix2x2[2]) > min_norm) {
isDiag = false;
complex timesI[4];
complex toApply[4];
for (bitLenInt i = 0; i < 4; i++) {
timesI[i] = complex(ZERO_R1, ONE_R1) * matrix2x2[i];
}

complex expOfGate[4];
complex jacobian[4];
complex inverseJacobian[4];
complex tempMatrix2x2[4];

// Diagonalize the matrix, if it is not already diagonal. Otherwise, copy it into the temporary matrix.
if (!isDiag) {
complex trace = matrix2x2[0] + matrix2x2[3];
complex determinant = (matrix2x2[0] * matrix2x2[3]) - (matrix2x2[1] * matrix2x2[2]);
complex quadraticRoot =
sqrt((matrix2x2[0] - matrix2x2[3]) * (matrix2x2[0] - matrix2x2[3]) - (real1)(4.0) * determinant);
complex eigenvalue1 = (trace + quadraticRoot) / (real1)2.0;
complex eigenvalue2 = (trace - quadraticRoot) / (real1)2.0;

if (norm(matrix2x2[1]) > min_norm) {
jacobian[0] = matrix2x2[1];
jacobian[2] = eigenvalue1 - matrix2x2[0];

jacobian[1] = matrix2x2[1];
jacobian[3] = eigenvalue2 - matrix2x2[0];
} else {
jacobian[0] = eigenvalue1 - matrix2x2[3];
jacobian[2] = matrix2x2[2];

jacobian[1] = eigenvalue2 - matrix2x2[3];
jacobian[3] = matrix2x2[2];
}

real1 nrm = std::sqrt(norm(jacobian[0]) + norm(jacobian[2]));
jacobian[0] /= nrm;
jacobian[2] /= nrm;

nrm = std::sqrt(norm(jacobian[1]) + norm(jacobian[3]));
jacobian[1] /= nrm;
jacobian[3] /= nrm;

determinant = (jacobian[0] * jacobian[3]) - (jacobian[1] * jacobian[2]);
inverseJacobian[0] = jacobian[3] / determinant;
inverseJacobian[1] = -jacobian[1] / determinant;
inverseJacobian[2] = -jacobian[2] / determinant;
inverseJacobian[3] = jacobian[0] / determinant;

matrix2x2Mul(matrix2x2, jacobian, tempMatrix2x2);
matrix2x2Mul(inverseJacobian, tempMatrix2x2, expOfGate);
Qrack::exp2x2(timesI, toApply);
if (antiCtrled) {
ApplyAntiControlledSingleBit(controls, controlLen, qubit, toApply);
} else {
std::copy(matrix2x2, matrix2x2 + 4, expOfGate);
ApplyControlledSingleBit(controls, controlLen, qubit, toApply);
}
}

// Note: For a (2x2) hermitian input gate, this theoretically produces a unitary output transformation.
expOfGate[0] =
((real1)exp(-imag(expOfGate[0]))) * complex((real1)cos(real(expOfGate[0])), (real1)sin(real(expOfGate[0])));
expOfGate[1] = complex(ZERO_R1, ZERO_R1);
expOfGate[2] = complex(ZERO_R1, ZERO_R1);
expOfGate[3] =
((real1)exp(-imag(expOfGate[3]))) * complex((real1)cos(real(expOfGate[3])), (real1)sin(real(expOfGate[3])));

if (!isDiag) {
matrix2x2Mul(expOfGate, inverseJacobian, tempMatrix2x2);
matrix2x2Mul(jacobian, tempMatrix2x2, expOfGate);
/// Logarithm of arbitrary single bit gate
void QInterface::Log(bitLenInt* controls, bitLenInt controlLen, bitLenInt qubit, complex* matrix2x2, bool antiCtrled)
{
complex toApply[4];
Qrack::log2x2(matrix2x2, toApply);
if (antiCtrled) {
ApplyAntiControlledSingleBit(controls, controlLen, qubit, toApply);
} else {
ApplyControlledSingleBit(controls, controlLen, qubit, toApply);
}

ApplyControlledSingleBit(controls, controlLen, qubit, expOfGate);
}

/// Exponentiate Pauli X operator
Expand Down
Loading

0 comments on commit fadccd5

Please sign in to comment.