linalg.cpp

/** @file
 * Testing utilities which perform linear algebra
 * routines upon reference qvector and qmatrix. 
 * These are slow, serial, un-optimised, defensively-
 * designed routines.
 *
 * @author Tyson Jones
 */
 
#include "qvector.hpp"
#include "qmatrix.hpp"
#include "linalg.hpp"
#include "macros.hpp"
#include "compare.hpp"
 
#include <algorithm>
#include <vector>
#include <limits>
 
using std::vector;
 
 
 
/*
 * SCALAR OPERATIONS
 */
 
 
int getLog2(qindex a) {
    DEMAND( a >= 0 );
    DEMAND( (a & (a - 1)) == 0 );  // is pow2
 
    int n = 0;
    while (a >>= 1)
        n++;
 
    return n;
}
 
 
qindex getPow2(int a) {
    DEMAND( a >= 0 );
 
    return ((qindex) 1) << a;
}
 
 
int getBitAt(qindex num, int ind) {
    DEMAND( num >= 0 );
 
    return (num >> ind) & 1;
}
 
 
vector<int> getBits(qindex num, int numBits) {
    DEMAND( numBits > 0 );
 
    // out ordered least to most significant
    vector<int> out(numBits);
 
    for (int i=0; i<numBits; i++)
        out[i] = getBitAt(num, i);
 
    return out;
}
 
 
qindex getBitsAt(qindex num, vector<int> inds) {
    DEMAND( num >= 0 );
 
    qindex out = 0;
 
    for (size_t i=0; i<inds.size(); i++)
        out |= getBitAt(num, inds[i]) << i;
    
    return out;
}
 
 
qindex setBitAt(qindex num, int ind, int bit) {
    DEMAND( num >= 0 );
 
    qindex one = 1;
    return (num & ~(one << ind)) | (bit << ind);
}
 
 
qindex setBitsAt(qindex num, vector<int> inds, qindex bits) {
 
    for (size_t i=0; i<inds.size(); i++)
        num = setBitAt(num, inds[i], getBitAt(bits, i));
 
    return num;
}
 
 
qindex getNumPermutations(int n, int k) {
    DEMAND( n >= k );
 
    constexpr auto max = std::numeric_limits<qindex>::max();
 
    // P(n, k) = n! / (n-k)!
    qindex p = 1;
    for (int t=n-k+1; t<=n; t++)
        p *= (p < max / t)? t : 0; // set to 0 on overflow
 
    return p;
}
 
 
 
/*
 * VECTOR OPERATIONS
 */
 
 
qcomp getSum(qvector vec) {
 
    qcomp sum = qcomp(0,0);
    qcomp y, t, c=sum;
    
    // complex Kahan summation
    for (auto& x : vec) {
        y = x - c;
        t = sum + y;
        c = ( t - sum ) - y;
        sum = t;
    }
 
    return sum;
}
 
 
qreal getSum(vector<qreal> vec) {
 
    // in = real(vec)
    qvector in = getZeroVector(vec.size());
    for (size_t i=0; i<in.size(); i++)
        in[i] = qcomp(vec[i],0);
 
    return std::real(getSum(in));
}
 
 
qvector getNormalised(qvector vec) {
 
    // prob[i] = abs(vec[i])^2
    vector<qreal> probs(vec.size());
    for (size_t i=0; i<vec.size(); i++)
        probs[i] = std::norm(vec[i]);
 
    // normalise vector
    qreal norm = getSum(probs);
    qreal fac = 1 / std::sqrt(norm);
    for (auto& x : vec)
        x *= fac;
 
    return vec;
}
 
 
qvector getDiscreteFourierTransform(qvector in, bool inverse) {
    DEMAND( in.size() > 0 );
    
    size_t dim = in.size();
    qvector out = getZeroVector(dim);
 
    // PI must be accurate here
    qreal pi = 3.14159265358979323846;
    qreal a = 1 / std::sqrt(dim);
    qreal b = (inverse ? -1 : 1) * 2 * pi / dim;
    
    for (size_t x=0; x<dim; x++)
        for (size_t y=0; y<dim; y++)
            out[x] += a * std::exp(b * x * y * 1_i) * in[y];
 
    return out;
}
 
 
qvector getDiscreteFourierTransform(qvector in, vector<int> targs, bool inverse) {
    DEMAND( in.size() > 0 );
 
    size_t dim = in.size();
    qvector out = getZeroVector(dim);
    
    qindex len = getPow2(targs.size());
    qreal pi = 3.14159265358979323846;
    qreal a = 1 / std::sqrt(len);
    qreal b = (inverse ? -1 : 1) * 2 * pi / len;
 
    for (size_t i=0; i<dim; i++) {
        size_t x = getBitsAt(i, targs);
        for (size_t y=0; y<len; y++) {
            qindex j = setBitsAt(i, targs, y);
            out[j] += a * std::exp(b * x * y * 1_i) * in[i];
        }
    }
 
    return out;
}
 
 
 
/*
 * VECTOR & VECTOR OPERATIONS
 */
 
 
qcomp getInnerProduct(qvector bra, qvector ket) {
    DEMAND( bra.size() == ket.size() );
 
    qcomp out = 0;
 
    for (size_t i=0; i<bra.size(); i++)
        out += std::conj(bra[i]) * ket[i];
 
    return out;
}
 
 
qmatrix getOuterProduct(qvector ket, qvector bra) {
    DEMAND( bra.size() == ket.size() );
 
    qmatrix out = getZeroMatrix(bra.size());
 
    for (size_t i=0; i<ket.size(); i++)
        for (size_t j=0; j<ket.size(); j++)
            out[i][j] = ket[i] * std::conj(bra[j]);
 
    return out;
}
 
 
 
/*
 * MATRIX OPERATIONS
 */
 
 
bool isDiagonal(qmatrix m) {
 
    for (size_t r=0; r<m.size(); r++)
        for (size_t c=0; c<m.size(); c++)
            if (r!=c && m[r][c] != 0_i)
                return false;
 
    return true;
}
 
 
bool isApproxUnitary(qmatrix m) {
 
    // should be identity
    qmatrix md = m * getConjugateTranspose(m);
    qmatrix id = getIdentityMatrix(m.size());
    return doMatricesAgree(md, id);
}
 
 
qcomp getTrace(qmatrix m) {
 
    qcomp out = 0;
    for (size_t r=0; r<m.size(); r++)
        out += m[r][r];
 
    return out;
}
 
 
qmatrix getTranspose(qmatrix m) {
 
    qmatrix out = getZeroMatrix(m.size());
 
    for (size_t r=0; r<m.size(); r++)
        for (size_t c=0; c<m.size(); c++)
            out[r][c] = m[c][r];
 
    return out;
}
 
 
qmatrix getConjugate(qmatrix m) {
 
    for (auto& row : m)
        for (auto& elem : row)
            elem = std::conj(elem);
 
    return m;
}
 
 
qmatrix getConjugateTranspose(qmatrix m) {
    DEMAND( m.size() > 0 );
 
    // unlike most functions which assume qmatrix
    // is square, this one cheekily handles when
    // 'm' is non-square, since necessary for
    // computing partial traces
 
    qmatrix out(m[0].size(), qvector(m.size()));
 
    for (size_t r=0; r<out.size(); r++)
        for (size_t c=0; c<out[0].size(); c++)
            out[r][c] = std::conj(m[c][r]);
 
    return out;
}
 
 
qmatrix getPowerOfDiagonalMatrix(qmatrix m, qcomp p) {
    DEMAND( isDiagonal(m) );
 
    qmatrix out = getZeroMatrix(m.size());
 
    // pow(qcomp,qcomp) introduces wildly erroneous
    // imaginary components when both base is real
    // and negative, and exponent is real and integer
    // (so ergo does not produce complex numbers).
    // We divert to real-pow in that scenario!
 
    for (size_t i=0; i<m.size(); i++) {
        bool mIsRe  = std::imag(m[i][i]) == 0;
        bool mIsNeg = std::real(m[i][i]) < 0;
        bool pIsRe  = std::imag(p) == 0;
        bool pIsInt = std::trunc(std::real(p)) == std::real(p);
 
        // use pow(qreal,qreal) or pow(qcomp,qcomp)
        out[i][i] = (mIsRe && mIsNeg && pIsRe && pIsInt)?
            qcomp(std::pow(std::real(m[i][i]), std::real(p)),0):
            std::pow(m[i][i], p);
    }
 
    return out;
}
 
 
qmatrix getExponentialOfDiagonalMatrix(qmatrix m) {
    DEMAND( isDiagonal(m) );
 
    qmatrix out = getZeroMatrix(m.size());
 
    for (size_t i=0; i<m.size(); i++)
        out[i][i] = std::exp(m[i][i]);
 
    return out;
}
 
 
qmatrix getExponentialOfPauliMatrix(qcomp arg, qmatrix m) {
    
    // exp(-i arg/2 m) where m = prod(paulis)
    qmatrix id = getIdentityMatrix(m.size());
    qmatrix out = std::cos(arg/2)*id - 1_i*std::sin(arg/2)*m;
    return out;
}
 
 
qmatrix getExponentialOfNormalisedPauliVector(qreal arg, qreal x, qreal y, qreal z) {
 
    // exp(-arg/2 i [x^ X + y^ Y + z^ Z])
    qreal n = std::sqrt(x*x + y*y + z*z);
    x /= n;
    y /= n;
    z /= n;
 
    qmatrix id = getIdentityMatrix(2);
    qmatrix out = std::cos(arg/2)*id - 1_i*std::sin(arg/2)*(
        x * getPauliMatrix(1) +
        y * getPauliMatrix(2) +
        z * getPauliMatrix(3));
 
    return out;
}
 
 
qmatrix getOrthonormalisedRows(qmatrix matr) {
 
    // perform the Gram-Schmidt process, processing each row of matr in-turn
    for (size_t i=0; i<matr.size(); i++) {
        qvector row = matr[i];
        
        // compute new orthogonal row by subtracting proj row onto prevs
        for (int k=i-1; k>=0; k--) {
 
            // compute inner_product(row, prev) = row . conj(prev)
            qcomp prod = getInnerProduct(matr[k], row);
                            
            // subtract (proj row onto prev) = (prod * prev) from final row
            matr[i] -= prod * matr[k];
        }
    
        // normalise the row 
        matr[i] = getNormalised(matr[i]);
    }
    
    // return the new orthonormal matrix
    return matr;
}
 
 
qmatrix getProjector(int outcome) {
    DEMAND( outcome == 0 || outcome == 1 );
 
    qmatrix out = getZeroMatrix(2);
    out[outcome][outcome] = 1.;
    return out;
}
 
 
qmatrix getProjector(vector<int> targets, vector<int> outcomes, int numQubits) {
    DEMAND( targets.size() == outcomes.size() );
    DEMAND( numQubits > *std::max_element(targets.begin(), targets.end()) );
 
    // prepare { |0><0|, I, I, |1><1|, ... }
    vector<qmatrix> matrices(numQubits, getIdentityMatrix(2));
    for (size_t i=0; i<targets.size(); i++)
        matrices[targets[i]] = getProjector(outcomes[i]);
 
    return getKroneckerProduct(matrices);
}
 
 
qmatrix getPartialTrace(qmatrix in, vector<int> targets) {
    DEMAND( in.size() > getPow2(targets.size()) );
 
    auto numTargs = targets.size();
    auto numQubits = getLog2(in.size());
    auto numTargVals = getPow2(numTargs);
 
    qmatrix out = getZeroMatrix(getPow2(numQubits - numTargs));
 
    for (qindex v=0; v<numTargVals; v++) {
 
        // prepare { |0>, I, I, |1>, ... }
        vector<qmatrix> matrices(numQubits, getIdentityMatrix(2));
        for (size_t t=0; t<numTargs; t++) {
            int bit = getBitAt(v, t);
            matrices[targets[t]] = {
                {bit? 0.:1.}, 
                {bit? 1.:0.}};
        }
 
        qmatrix ket = getKroneckerProduct(matrices);
        qmatrix bra = getConjugateTranspose(ket);
        out += bra * in * ket;
    }
 
    return out;
}
 
 
qmatrix getControlledMatrix(qmatrix matrix, int numCtrls) {
 
    size_t dim = getPow2(numCtrls) * matrix.size();
    size_t off = dim - matrix.size();
    
    qmatrix out = getIdentityMatrix(dim);
    setSubMatrix(out, matrix, off, off);
 
    return out;
}
 
 
qmatrix getMixture(vector<qmatrix> densmatrs, vector<qreal> probs) {
    DEMAND( densmatrs.size() > 0 );
 
    qmatrix out = getZeroMatrix(densmatrs[0].size());
    for (size_t i=0; i<densmatrs.size(); i++)
        out += probs[i] * densmatrs[i];
 
    return out;
}
 
 
qmatrix getMixture(vector<qvector> statevecs, vector<qreal> probs) {
 
    vector<qmatrix> densmatrs(statevecs.size());
    for (size_t i=0; i<statevecs.size(); i++)
        densmatrs[i] = getOuterProduct(statevecs[i], statevecs[i]);
 
    return getMixture(densmatrs, probs);
}
 
 
qmatrix getSuperOperator(vector<qmatrix> matrices) {
    DEMAND( matrices.size() > 0 );
 
    size_t dim = matrices[0].size();
 
    // out = sum_m conj(m) (x) m
    qmatrix out = getZeroMatrix(dim * dim);
    for (auto& matr : matrices)
        out += getKroneckerProduct(getConjugate(matr), matr);
 
    return out;
}
 
 
 
/*
 * MATRIX & VECTOR OPERATIONS
 */
 
 
qvector operator * (const qmatrix& m, const qvector& v) {
    DEMAND( m.size() == v.size() );
 
    qvector out = getZeroVector(v.size());
 
    for (size_t r=0; r<v.size(); r++)
        for (size_t c=0; c<v.size(); c++)
            out[r] += m[r][c] * v[c];
 
    return out;
}
 
 
 
/*
 * MATRIX & MATRIX OPERATIONS
 */
 
 
qmatrix getKroneckerProduct(qmatrix a, qmatrix b) {
 
    // we permit the matrices to be non-square which is
    // pretty cheeky (since qmatrix is assumed square with
    // a 2^N dimension by most other functions), but is
    // necessary for us to compute partial traces
 
    size_t aRows = a.size();
    size_t bRows = b.size();
    size_t aCols = a[0].size();
    size_t bCols = b[0].size();
 
    qmatrix out(aRows * bRows, qvector(aCols * bCols));
    
    for (size_t r=0; r<bRows; r++)
        for (size_t c=0; c<bCols; c++)
            for (size_t i=0; i<aRows; i++)
                for (size_t j=0; j<aCols; j++)
                    out[r+bRows*i][c+bCols*j] = a[i][j] * b[r][c];
 
    return out;
}
 
 
qmatrix getKroneckerProduct(vector<qmatrix> matrices) {
 
    qmatrix out = getIdentityMatrix(1);
 
    // matrices[n-1] (x) ... (x) matrices[0]
    for (auto& m : matrices)
        out = getKroneckerProduct(m, out);
 
    return out;
}
 
 
qmatrix getKroneckerProduct(qmatrix m, int count) {
    DEMAND( count >= 1 );
 
    qmatrix out = getIdentityMatrix(1);
 
    for (int n=0; n<count; n++)
        out = getKroneckerProduct(out, m);
 
    return out;
}
 
 
 
/*
 * MATRIX COLLECTIONS
 */
 
 
bool isApproxCPTP(vector<qmatrix> matrices) {
    DEMAND( matrices.size() >= 1 );
 
    size_t dim = matrices[0].size();
    qmatrix id = getIdentityMatrix(dim);
    qmatrix sum = getZeroMatrix(dim);
 
    for (auto& m : matrices)
        sum += getConjugateTranspose(m) * m;
 
    return doMatricesAgree(sum, id);
}