Function getRandomKrausMap

tests/deprecated/test_utilities.cpp:713–751

}

vector<QMatrix> getRandomKrausMap(int numQb, int numOps) {
    DEMAND( numOps >= 1 );
    DEMAND( numOps <= 4*numQb*numQb );

    // generate random unitaries
    vector<QMatrix> ops;
    for (int i=0; i<numOps; i++)
        ops.push_back(getRandomUnitary(numQb));
        
    // generate random weights
    vector<qreal> weights(numOps);
    for (int i=0; i<numOps; i++)
        weights[i] = getRandomReal(0, 1);
        
    // normalise random weights
    qreal weightSum = 0;
    for (int i=0; i<numOps; i++)
        weightSum += weights[i];
    for (int i=0; i<numOps; i++)
        weights[i] = sqrt((qreal) weights[i]/weightSum);
        
    // normalise ops
    for (int i=0; i<numOps; i++)
        ops[i] *= weights[i];
        
    // check what we produced was a valid Kraus map
    QMatrix iden = getIdentityMatrix(1 << numQb);
    QMatrix prodSum = getZeroMatrix(1 << numQb);
    for (int i=0; i<numOps; i++)
        prodSum += getConjugateTranspose(ops[i]) * ops[i];

    // in the rare scenario it is insufficiently numerically precise,
    // replace the map with trivially precise diagonals
    if( ! areEqual(prodSum, iden) )
        for (int i=0; i<numOps; i++)
            ops[i] = weights[i] * getRandomDiagonalUnitary(numQb);

    return ops;
}

vector<QVector> getRandomOrthonormalVectors(int numQb, int numStates) {
    DEMAND( numQb >= 1 );