Dave....... I started something again months ago....
Help I'm not sure this is even close to being correct...
I also want to know what the output is if the input is
|ψ₁⟩ ⊗ |ψ₂⟩ ⊗ |ψ₃⟩ ⊗ |ψ₄⟩ = (α₁|00⟩ + β₁|01⟩ + γ₁|10⟩ + δ₁|11⟩) ⊗ (α₂|00⟩ + β₂|01⟩ + γ₂|10⟩ + δ₂|11⟩) ⊗ (α₃|00⟩ + β₃|01⟩ + γ₃|10⟩ + δ₃|11⟩) ⊗ (α₄|00⟩ + β₄|01⟩ + γ₄|10⟩ + δ₄|11⟩)
import org.apache.commons.math3.complex.Complex;
import org.apache.commons.math3.linear.Array2DRowFieldMatrix;
import org.apache.commons.math3.linear.FieldMatrix;
import org.apache.commons.math3.linear.MatrixUtils;
import org.apache.commons.math3.linear.RandomMatrixUtils;
public class QuantumRandomization {
public static void main(String[] args) {
// Assign coefficient values
double alpha1 = 0.5;
double beta1 = 0.3;
double gamma1 = 0.1;
double delta1 = 0.2;
double alpha2 = 0.2;
double beta2 = 0.4;
double gamma2 = 0.6;
double delta2 = 0.8;
double alpha3 = 0.7;
double beta3 = 0.9;
double gamma3 = 0.3;
double delta3 = 0.5;
double alpha4 = 0.4;
double beta4 = 0.2;
double gamma4 = 0.8;
double delta4 = 0.6;
// Create four 2-state quantum systems
FieldMatrix system1 = new Array2DRowFieldMatrix<>(new Complex[][] {
{ new Complex(alpha1, 0), new Complex(beta1, 0) },
{ new Complex(gamma1, 0), new Complex(delta1, 0) }
});
FieldMatrix system2 = new Array2DRowFieldMatrix<>(new Complex[][] {
{ new Complex(alpha2, 0), new Complex(beta2, 0) },
{ new Complex(gamma2, 0), new Complex(delta2, 0) }
});
FieldMatrix system3 = new Array2DRowFieldMatrix<>(new Complex[][] {
{ new Complex(alpha3, 0), new Complex(beta3, 0) },
{ new Complex(gamma3, 0), new Complex(delta3, 0) }
});
FieldMatrix system4 = new Array2DRowFieldMatrix<>(new Complex[][] {
{ new Complex(alpha4, 0), new Complex(beta4, 0) },
{ new Complex(gamma4, 0), new Complex(delta4, 0) }
});
// Take the tensor product of the four systems
FieldMatrix tensorProduct = tensorProduct(tensorProduct(tensorProduct(system1, system2), system3), system4);
// Randomize the tensor product using a unitary matrix
FieldMatrix randomizedTensorProduct = tensorProduct.preMultiply(RandomMatrixUtils.createUnitaryMatrix(tensorProduct.getRowDimension()));
// Measure the randomized tensor product, returning a random state
FieldMatrix measuredTensorProduct = measureAndNormalize(randomizedTensorProduct);
// Print the measured tensor product
System.out.println(measuredTensorProduct);
}
// Compute the tensor product of two matrices
public static FieldMatrix tensorProduct(FieldMatrix matrix1, FieldMatrix matrix2) {
int rows1 = matrix1.getRowDimension();
int cols1 = matrix1.getColumnDimension();
int rows2 = matrix2.getRowDimension();
int cols2 = matrix2.getColumnDimension();
FieldMatrix result = new Array2DRowFieldMatrix<>(rows1 * rows2, cols1 * cols2);
for (int i = 0; i < rows1; i++) {
for (int j = 0; j < cols1; j++) {
for (int k = 0; k < rows2; k++) {
for (int l = 0; l < cols2; l++) {
result.setEntry(i * rows2 + k, j * cols2 + l, matrix1.getEntry(i, j).multiply(matrix2.getEntry(k, l)));
}
}
}
}
return result;
}
// Measure the matrix and normalize
public static FieldMatrix measureAndNormalize(FieldMatrix matrix) {
// Compute the squared absolute
double[] probabilities = new double[matrix.getRowDimension()];
for (int i = 0; i < matrix.getRowDimension(); i++) {
Complex entry = matrix.getEntry(i, 0);
probabilities[i] = entry.real() * entry.real() + entry.imag() * entry.imag();
}
// Sample a random state from matrix
int randomIndex = sampleFromDistribution(probabilities);
// Create a new matrix
FieldMatrix measuredMatrix = new Array2DRowFieldMatrix<>(matrix.getField(), matrix.getRowDimension(), matrix.getColumnDimension());
measuredMatrix.setEntry(randomIndex, 0, Complex.ONE);
// Normalize
double norm = MatrixUtils.createFieldVector(measuredMatrix.getColumn(0)).getNorm();
return measuredMatrix.scalarMultiply(1 / norm);
}
public static int sampleFromDistribution(double[] probabilities) {
double randomValue = Math.random();
double cumulativeProbability = 0;
int index = 0;
while (cumulativeProbability < randomValue && index < probabilities.length) {
cumulativeProbability += probabilities[index];
index++;
}
return index - 1;
}
}