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	getDiscreteFourierTransform (qvector in, bool inverse)
qvector	getDiscreteFourierTransform (qvector in, vector< int > targs, bool inverse)
qmatrix	getExponentialOfDiagonalMatrix (qmatrix)
qmatrix	getExponentialOfNormalisedPauliVector (qreal arg, qreal x, qreal y, qreal z)
qmatrix	getExponentialOfPauliMatrix (qcomp 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)
qindex	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 48 of file linalg.cpp.

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

getBits()

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

Definition at line 55 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 68 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 283 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 293 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 453 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;
}

getDiscreteFourierTransform() [1/2]

qvector getDiscreteFourierTransform	(	qvector	in,
		bool	inverse )

Definition at line 162 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 = (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;
}

getDiscreteFourierTransform() [2/2]

qvector getDiscreteFourierTransform	(	qvector	in,
		vector< int >	targs,
		bool	inverse )

Definition at line 181 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 = (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;
}

getExponentialOfDiagonalMatrix()

qmatrix getExponentialOfDiagonalMatrix

(

qmatrix

m

)

Definition at line 338 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 359 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	(	qcomp	arg,
		qmatrix	pauli )

Definition at line 350 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 210 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 561 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 525 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 549 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 29 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 465 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 476 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 145 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()

qindex getNumPermutations	(	int	n,
		int	k )

Definition at line 97 of file linalg.cpp.

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

getOrthonormalisedRows()

qmatrix getOrthonormalisedRows

(

qmatrix

matr

)

Definition at line 377 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 222 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 424 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 41 of file linalg.cpp.

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

getPowerOfDiagonalMatrix()

qmatrix getPowerOfDiagonalMatrix	(	qmatrix	diag,
		qcomp	power )

Definition at line 311 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 402 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 411 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 117 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 134 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 486 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 261 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 271 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 579 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 252 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 241 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 506 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 80 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 88 of file linalg.cpp.

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