Function getRandomUnitary

tests/deprecated/test_utilities.cpp:677–711

}

QMatrix getRandomUnitary(int numQb) {
    DEMAND( numQb >= 1 );

    // create Z ~ random complex matrix (distribution not too important)
    size_t dim = 1 << numQb;
    QMatrix matrZ = getRandomQMatrix(dim);
    QMatrix matrZT = getTranspose(matrZ);

    // create Z = Q R (via QR decomposition) ...
    QMatrix matrQT = getOrthonormalisedRows(matrZ);
    QMatrix matrQ = getTranspose(matrQT);
    QMatrix matrR = getZeroMatrix(dim);

    // ... where R_rc = (columm c of Z) . (column r of Q) = (row c of ZT) . (row r of QT)
    for (size_t r=0; r<dim; r++)
        for (size_t c=r; c<dim; c++)
            matrR[r][c] = matrZT[c] * matrQT[r];

    // create D = normalised diagonal of R
    QMatrix matrD = getZeroMatrix(dim);
    for (size_t i=0; i<dim; i++)
        matrD[i][i] = matrR[i][i] / abs(matrR[i][i]);

    // create U = Q D
    QMatrix matrU = matrQ * matrD;

    // in the rare scenario the result is not sufficiently precisely unitary,
    // replace it with a trivially unitary diagonal matrix
    QMatrix daggerProd = matrU * getConjugateTranspose(matrU);
    QMatrix iden = getIdentityMatrix(dim);
    if( ! areEqual(daggerProd, iden) )
        matrU = getRandomDiagonalUnitary(numQb);

    return matrU;
}

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