Linalg

Testing utilities which perform linear algebra routines upon reference qvector and qmatrix. These are slow, serial, un-optimised, defensively- designed routines. More...

Functions
int	getBitAt (qindex num, int ind)
vector< int >	getBits (qindex num, int numBits)
qindex	getBitsAt (qindex num, vector< int > inds)
qmatrix	getConjugate (qmatrix)
qmatrix	getConjugateTranspose (qmatrix)
qmatrix	getControlledMatrix (qmatrix matrix, int numCtrls)
qvector	getDisceteFourierTransform (qvector in, vector< int > targs)
qvector	getDisceteFourierTransform (qvector)
qmatrix	getExponentialOfDiagonalMatrix (qmatrix)
qmatrix	getExponentialOfNormalisedPauliVector (qreal arg, qreal x, qreal y, qreal z)
qmatrix	getExponentialOfPauliMatrix (qreal arg, qmatrix pauli)
qcomp	getInnerProduct (qvector bra, qvector ket)
qmatrix	getKroneckerProduct (qmatrix, int count)
qmatrix	getKroneckerProduct (qmatrix, qmatrix)
qmatrix	getKroneckerProduct (vector< qmatrix >)
int	getLog2 (qindex)
qmatrix	getMixture (vector< qmatrix > densmatrs, vector< qreal > probs)
qmatrix	getMixture (vector< qvector > statevecs, vector< qreal > probs)
qvector	getNormalised (qvector)
int	getNumPermutations (int n, int k)
qmatrix	getOrthonormalisedRows (qmatrix)
qmatrix	getOuterProduct (qvector ket, qvector bra)
qmatrix	getPartialTrace (qmatrix matrix, vector< int > targets)
qindex	getPow2 (int)
qmatrix	getPowerOfDiagonalMatrix (qmatrix diag, qcomp power)
qmatrix	getProjector (int outcome)
qmatrix	getProjector (vector< int > targets, vector< int > outcomes, int numQubits)
qcomp	getSum (qvector)
qreal	getSum (vector< qreal > vec)
qmatrix	getSuperOperator (vector< qmatrix >)
qcomp	getTrace (qmatrix)
qmatrix	getTranspose (qmatrix)
bool	isApproxCPTP (vector< qmatrix >)
bool	isApproxUnitary (qmatrix)
bool	isDiagonal (qmatrix)
qvector	operator* (const qmatrix &, const qvector &)
qindex	setBitAt (qindex num, int ind, int bit)
qindex	setBitsAt (qindex num, vector< int > inds, qindex bits)

Detailed Description

Testing utilities which perform linear algebra routines upon reference qvector and qmatrix. These are slow, serial, un-optimised, defensively- designed routines.

Function Documentation

getBitAt()

int getBitAt	(	qindex	num,
		int	ind )

Definition at line 47 of file linalg.cpp.

                                  {
    DEMAND( num >= 0 );
 
    return (num >> ind) & 1;
}

getBits()

vector< int > getBits	(	qindex	num,
		int	numBits )

Definition at line 54 of file linalg.cpp.

                                             {
    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;
}

getBitsAt()

qindex getBitsAt	(	qindex	num,
		vector< int >	inds )

Definition at line 67 of file linalg.cpp.

                                               {
    DEMAND( num >= 0 );
 
    qindex out = 0;
 
    for (size_t i=0; i<inds.size(); i++)
        out |= getBitAt(num, inds[i]) << i;
    
    return out;
}

getConjugate()

qmatrix getConjugate

(

qmatrix

m

)

Definition at line 281 of file linalg.cpp.

                                {
 
    for (auto& row : m)
        for (auto& elem : row)
            elem = std::conj(elem);
 
    return m;
}

getConjugateTranspose()

qmatrix getConjugateTranspose

(

qmatrix

m

)

Definition at line 291 of file linalg.cpp.

                                         {
    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;
}

getControlledMatrix()

qmatrix getControlledMatrix	(	qmatrix	matrix,
		int	numCtrls )

Definition at line 451 of file linalg.cpp.

                                                          {
 
    size_t dim = getPow2(numCtrls) * matrix.size();
    size_t off = dim - matrix.size();
    
    qmatrix out = getIdentityMatrix(dim);
    setSubMatrix(out, matrix, off, off);
 
    return out;
}

getDisceteFourierTransform() [1/2]

qvector getDisceteFourierTransform	(	qvector	in,
		vector< int >	targs )

Definition at line 179 of file linalg.cpp.

                                                                  {
    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 = 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;
}

getDisceteFourierTransform() [2/2]

qvector getDisceteFourierTransform

(

qvector

in

)

Definition at line 160 of file linalg.cpp.

                                               {
    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 = 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;
}

getExponentialOfDiagonalMatrix()

qmatrix getExponentialOfDiagonalMatrix

(

qmatrix

m

)

Definition at line 336 of file linalg.cpp.

                                                  {
    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;
}

getExponentialOfNormalisedPauliVector()

qmatrix getExponentialOfNormalisedPauliVector	(	qreal	arg,
		qreal	x,
		qreal	y,
		qreal	z )

Definition at line 357 of file linalg.cpp.

                                                                                    {
 
    // 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;
}

getExponentialOfPauliMatrix()

qmatrix getExponentialOfPauliMatrix	(	qreal	arg,
		qmatrix	pauli )

Definition at line 348 of file linalg.cpp.

                                                          {
    
    // 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;
}

getInnerProduct()

qcomp getInnerProduct	(	qvector	bra,
		qvector	ket )

Definition at line 208 of file linalg.cpp.

                                                {
    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;
}

getKroneckerProduct() [1/3]

qmatrix getKroneckerProduct	(	qmatrix	m,
		int	count )

Definition at line 559 of file linalg.cpp.

                                                  {
    DEMAND( count >= 1 );
 
    qmatrix out = getIdentityMatrix(1);
 
    for (int n=0; n<count; n++)
        out = getKroneckerProduct(out, m);
 
    return out;
}

getKroneckerProduct() [2/3]

qmatrix getKroneckerProduct	(	qmatrix	a,
		qmatrix	b )

Definition at line 523 of file linalg.cpp.

                                                  {
 
    // 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;
}

getKroneckerProduct() [3/3]

qmatrix getKroneckerProduct

(

vector< qmatrix >

matrices

)

Definition at line 547 of file linalg.cpp.

                                                      {
 
    qmatrix out = getIdentityMatrix(1);
 
    // matrices[n-1] (x) ... (x) matrices[0]
    for (auto& m : matrices)
        out = getKroneckerProduct(m, out);
 
    return out;
}

getLog2()

int getLog2

(

qindex

a

)

Definition at line 28 of file linalg.cpp.

                      {
    DEMAND( a >= 0 );
    DEMAND( (a & (a - 1)) == 0 );  // is pow2
 
    int n = 0;
    while (a >>= 1)
        n++;
 
    return n;
}

getMixture() [1/2]

qmatrix getMixture	(	vector< qmatrix >	densmatrs,
		vector< qreal >	probs )

Definition at line 463 of file linalg.cpp.

                                                                   {
    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;
}

getMixture() [2/2]

qmatrix getMixture	(	vector< qvector >	statevecs,
		vector< qreal >	probs )

Definition at line 474 of file linalg.cpp.

                                                                   {
 
    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);
}

getNormalised()

qvector getNormalised

(

qvector

vec

)

Definition at line 143 of file linalg.cpp.

                                   {
 
    // 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;
}

getNumPermutations()

int getNumPermutations	(	int	n,
		int	k )

Definition at line 96 of file linalg.cpp.

                                     {
    DEMAND( n >= k );
    DEMAND( n <= 11 ); // else int overflow
 
    // P(n, k) = n! / (n-k)!
    qindex p = 1;
    for (int t=n-k+1; t<=n; t++)
        p *= t;
 
    return p;
}

getOrthonormalisedRows()

qmatrix getOrthonormalisedRows

(

qmatrix

matr

)

Definition at line 375 of file linalg.cpp.

                                             {
 
    // 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;
}

getOuterProduct()

qmatrix getOuterProduct	(	qvector	ket,
		qvector	bra )

Definition at line 220 of file linalg.cpp.

                                                  {
    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;
}

getPartialTrace()

qmatrix getPartialTrace	(	qmatrix	matrix,
		vector< int >	targets )

Definition at line 422 of file linalg.cpp.

                                                         {
    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;
}

getPow2()

qindex getPow2

(

int

a

)

Definition at line 40 of file linalg.cpp.

                      {
    DEMAND( a >= 0 );
 
    return ((qindex) 1) << a;
}

getPowerOfDiagonalMatrix()

qmatrix getPowerOfDiagonalMatrix	(	qmatrix	diag,
		qcomp	power )

Definition at line 309 of file linalg.cpp.

                                                     {
    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;
}

getProjector() [1/2]

qmatrix getProjector

(

int

outcome

)

Definition at line 400 of file linalg.cpp.

                                  {
    DEMAND( outcome == 0 || outcome == 1 );
 
    qmatrix out = getZeroMatrix(2);
    out[outcome][outcome] = 1.;
    return out;
}

getProjector() [2/2]

qmatrix getProjector	(	vector< int >	targets,
		vector< int >	outcomes,
		int	numQubits )

Definition at line 409 of file linalg.cpp.

                                                                               {
    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);
}

getSum() [1/2]

qcomp getSum

(

qvector

vec

)

Definition at line 115 of file linalg.cpp.

                          {
 
    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;
}

getSum() [2/2]

qreal getSum

(

vector< qreal >

vec

)

Definition at line 132 of file linalg.cpp.

                                {
 
    // 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));
}

getSuperOperator()

qmatrix getSuperOperator

(

vector< qmatrix >

matrices

)

Definition at line 484 of file linalg.cpp.

                                                   {
    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;
}

getTrace()

qcomp getTrace

(

qmatrix

m

)

Definition at line 259 of file linalg.cpp.

                          {
 
    qcomp out = 0;
    for (size_t r=0; r<m.size(); r++)
        out += m[r][r];
 
    return out;
}

getTranspose()

qmatrix getTranspose

(

qmatrix

m

)

Definition at line 269 of file linalg.cpp.

                                {
 
    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;
}

isApproxCPTP()

bool isApproxCPTP

(

vector< qmatrix >

matrices

)

Definition at line 577 of file linalg.cpp.

                                            {
    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);
}

isApproxUnitary()

bool isApproxUnitary

(

qmatrix

m

)

Definition at line 250 of file linalg.cpp.

                                {
 
    // should be identity
    qmatrix md = m * getConjugateTranspose(m);
    qmatrix id = getIdentityMatrix(m.size());
    return doMatricesAgree(md, id);
}

isDiagonal()

bool isDiagonal

(

qmatrix

m

)

Definition at line 239 of file linalg.cpp.

                           {
 
    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;
}

operator*()

qvector operator*	(	const qmatrix &	m,
		const qvector &	v )

Definition at line 504 of file linalg.cpp.

                                                        {
    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;
}

setBitAt()

qindex setBitAt	(	qindex	num,
		int	ind,
		int	bit )

Definition at line 79 of file linalg.cpp.

                                              {
    DEMAND( num >= 0 );
 
    qindex one = 1;
    return (num & ~(one << ind)) | (bit << ind);
}

setBitsAt()

qindex setBitsAt	(	qindex	num,
		vector< int >	inds,
		qindex	bits )

Definition at line 87 of file linalg.cpp.

                                                            {
 
    for (size_t i=0; i<inds.size(); i++)
        num = setBitAt(num, inds[i], getBitAt(bits, i));
 
    return num;
}