Learning Q#. Statistical comparison of two sequences of numbers

Learning Q#. Statistical comparison of two sequences of numbers

Welcome to the new world of new computing technologies!

In everyday life, when we look at different objects, we try to understand whether they are similar or not, and how similar they are.

So it is in mathematics – when we look at sequences of numbers, we try to understand – whether they are similar or not, and how similar they are.

One of such criteria of “similarity” is the coincidence of the frequency characteristics of these sequences.

Let’s consider the question of how to sell such a check using quantum computing and write a test program on Q-sharp to check these considerations.

To understand this tutorial, you will need a basic knowledge of

  • probability theory

  • algebra

  • Boolean functions

  • convolution, correlation, scalar product

  • quantum computing (qubits and transformations)

  • programming on Q-sharp

Obtaining frequency characteristics of a numerical sequence

Let be the given sequence x(i) integers with 0..m-1where i=0..l-1

Let’s ask ourselves – is it possible to get a state for a register (array) of qubits SUM SQRT (P (v)) |where P(v) – Symbol frequency v in sequence x(i)?

The answer is why not (but with caveats)!

What does the record mean? x(i)where i=0..l-1 ?

This effectively means that we have a mapping (function) x: 0..l-1 to 0..m-1

Let

class – the number of bits sufficient to store numbers 0..l-1
km – the number of bits sufficient to store numbers 0..m-1

Convert the display x: 0..l-1 to 0..m-1 into a transformation for quantum computing Ux(i,z)=(i,z xor x(i))
where i – register with class qubits, and Z – register with km qubits

It is obvious that if you perform the transformation Ux(H|0>,|0>)then (since i runs the value from 0..2^kl-1) in the register Z a calculation of the statistics of the appearance of various symbols in the sequence will be formed.

But there is a small problem – we do not know the values ​​of the function x on the range of the argument l..2^kl-1.

Are there any solutions? Of course!

The easiest option is to supplement the sequence on the site l..2^kl-1 random equally probable numbers with 0..m-1

And voila…

We have x:0..2^kl-1 to 0..m-1 and Ux(i,z)=(i,z xor x(i))

  1. We set the initial state (i, z) = (H | 0>, | 0>)

  2. We perform the transformation Ux

  3. And in the register Z contains the frequency characteristics of the sequence (well, slightly averaged, the calculation of the generated tail from random, equally probable numbers)

How to implement a Ux transformation?

Let us have a basic set of gates. X and Controlled X

Let x|j:0..2^kl-1 to 0..1 – The reflection obtained by taking Jth bit of the value x:0..2^kl-1 to 0..m-1

That is, x = SUM 2^j*x|j

Like any Boolean function, x|j can be presented in a distributive form, which means Ux|j can be calculated as Ux|j(i,z|j) = SUM x|j(k)(Controlled X(X(i) xor k), z|j)

N.B. We will not now consider the issue of minimizing the number of hets, for example by computing the minimal distribution form

What’s next?

We understood that having an arbitrary sequence of numbers, it is possible to form the state of the qubit register, which will reflect the frequency characteristics of this sequence.

Let us now consider two sequences of numbers from 0..m-1

  • x: 0..lx-1 to 0..m-1

  • y: 0..ly-1 to 0..m-1

Can we, using quantum computing, quickly check whether the frequency characteristics of these two sequences match or not?

Let k – the number of bits sufficient to store numbers from the range 0..max(lx,ly)-1

In the same way as before, we will verify the mapping of x and y on the missing range

And we will build the corresponding transformations

  • Ux(i,z)=(i,z xor x(i))

  • Uy(i,z)=(i,z xor y(i))

Calculation of the scalar product through inversion and reconciliation

Let us have a wedding ring with surgery Op and a couple of functions a and bwhich display a ring in the field, for example, of complex numbers

Reconciliation operation a and b we will call the function c = axbsuch that c(i) = SUB a(ix)b(iy)|i=ix op iy

Let us have a wedding ring with surgery Op and function fwhich display a ring in the field, for example, of complex numbers

Let’s call the function the inversion of the argument inv f such that (inv f)(i) = f(-i)where i op -i == 0

Scalar product a and b called (a, b) = SUM a (i) b (i) (=(a,b)(0) where (a, b) (i) = SUM a (ix) b (iy) | i = ix-iy)

It is easy to see for yourself that (a, b) = ax (inv b)

Let for the qubit registers a and b

  • a = SUM A(i)|i>

  • b = SUM B(i)|i>

Then, if Op:

  • xor: (a,b)(i) = SUM A(ix)B(iy)|i>|i=ix xor ~iy

  • add: (a,b)(i) = SUM A(ix)B(iy)|i>|i=ix + (-iy)

Which gives us the calculation of the scalar product function a and b?

The answer is donkeys A(i)~=B(i) to all i and the distributions are different from trivial ones, then the value (a,b)(0) is the maximum (modulo) among other values ​​of this function
So, if we have a pair of qubit registers containing the frequency characteristics of two sequences, then

  1. computing the convolution of these two registers,

  2. checking that the resulting value |0> has the maximum probability

  3. will give us the answer that the two statistics coincide

and if not – then as a bonus – value |j>which has the maximum probability, is most likely to be a “shift” when generating a second distribution from the first (for example, when encrypting a gambit with a repeated key, that is, it will turn out to be the cipher key)

What will the calculation of the scalar product of two qubit registers look like?

Very simple – for example, for xor – it will be a transformation scalar(a,b)=(a,a xor X(b))and for addscalar (a, b) = (a, a add X (b) add | 1>)

Thus, the following scheme of the algorithm was obtained https://drive.google.com/file/d/1KIl4LVzPNoPKOQoURHmkVilz_Czxe4tC/view?usp=sharing

Statistical comparison of two sequences

Let’s check the reasoning by writing a program on Q-sharp and running tests

  1. We need arithmetic operations

  • increment of the register (ordered array) of qubits

  • summing the values ​​of individual qubits into a register (ordered array) of qubits

  • measuring the value of the register (ordered array) of qubits in the form of an integer

  • changing the value of the register (ordered array) of qubits according to the binary representation of an integer (xor)

  • copying a qubit register to another qubit register

    /// # Описание
    /// увеличение на единицу числового значения в массиве кубитов (рассматриваемых как регистр)
    /// то есть трансформация вида |k> -> |k+1>
    operation Inc(target: Qubit[]) : Unit is Ctl {
        let n = Length(target);
        for idx in 1..n {
            Controlled X(target[0..n-idx-1], target[n-idx]);
        } 
    }

    /// # Описание
    /// Реализация операции суммирования, то есть трансформации |abcde...> -> |a+b+c+d+e+...>
    /// Это не самая эффективная реализация, но самая простая для кодирования
    operation Sum(qubits: Qubit[], register: Qubit[]) : Unit {
        for q in qubits {
            Controlled Inc([q], register);
        }
    }
    
    /// # Описание
    /// измерение значений (коллапсирование) кубитов в массиве (который рассматриваем как один регистр)
    /// и возврат числа (равного полученной двоичной последовательности)
    operation Measure(qubits: Qubit[]) : Int {
        let results = ForEach(M, qubits);
        let i = ResultArrayAsInt(results);
        return i;
    }

    /// # Описание
    /// побитовое изменение на указанную величину числового значения в массиве кубитов (рассматриваемых как регистр)
    /// то есть трансформация вида |k> -> |k xor value>
    operation Xor(target: Qubit[], value: Int) : Unit {
        let n = Length(target);
        let bools = IntAsBoolArray(value, n);
        for idx in 0..n-1 {
            if(bools[idx]) {
                X(target[idx]);
            }
        }
    }

    /// # Описание
    /// вспомогательный метод для копирования значений массива кубитов
    operation Copy(source: Qubit[], target: Qubit[]) : Unit {
        let n = Length(source);
        for i in 0..(n-1) {
            Controlled X([source[i]], target[i]);
        }
    }
  1. Statistical methods

  • obtaining values ​​with a binomial distribution in the qubit register

  • random value generation with binomial distribution (and xor shift by a given value)

  • generation of a random value with equal probability distribution

    /// # Описание
    /// Реализация операции биноминального распределения, то есть n -> SUM SQRT(C(k,n))|k>
    operation Binom(n: Int, register: Qubit[]) : Unit {
        use qubits = Qubit[n] {
            ApplyToEach(H, qubits);
            Sum(qubits, register);
            ResetAll(qubits);
        }
    }

    /// # Описание
    /// Генерация случайного целого числа из диапазона 0..m-1
    /// Равновероятное распределение
    operation RandomRandom(m: Int) : Int {
        let k = BitSizeI(m-1);
        mutable i = 0;
        repeat {
            use qubits = Qubit[k] {
                ApplyToEach(H, qubits);
                set i = Measure(qubits);
                ResetAll(qubits);
            }
        }
        until (0<=i and i<m);
        return i;
    }

    /// # Описание
    /// Генерация случайного целого числа из диапазона 0..m-1
    /// Биноминальное распределение
    /// Дополнительный параметр a позволяет сместить распределение
    operation RandomBinom(a: Int, m: Int) : Int {
        let k = BitSizeI(m-1);
        use register = Qubit[k] {
            Binom(m-1, register); 
            let i = Measure(register);
            ResetAll(register);
            return Microsoft.Quantum.Bitwise.Xor(i,a)  % m;
        }
    }
  1. Methods of implementation of the described algorithm

  • applying the transformation Ux on the given array x

  • calculation of the scalar product for two registers

  • finding the index of the array element with the largest value

  • Repeatedly measuring the bias of the distribution and returning the best (most likely) value of the bias

    /// # Описание
    /// Выполнение трансформации для получения частотных характеристик последовательности чисел 0..m-1
    operation U(m: Int, l: Int, x: Int[], z: Qubit[]) : Unit {
        let kl = BitSizeI(l-1);
        let km = Length(z);
        let total = 2^kl;

        use (index) = (Qubit[kl]) {
            ApplyToEach(H, index);
            for i in 0..total-1 {
                Xor(index, i);
                mutable value = 0;
                if(i < l) {
                    set value = x[i];
                }
                else {
                    set value = RandomRandom(m);
                }
                let bools = IntAsBoolArray(value, km);
                for j in 0..km-1 {
                    if(bools[j]) {
                        Controlled X(index, z[j]);
                    }
                }
                Xor(index, i);
            }
            ResetAll(index);
        }
    }
    
    /// # Описание
    /// Вычисление скалярного произведения содержимого регистров кубитов
    operation ScalarTransformation(a: Qubit[], b: Qubit[]) : Unit
    {
        let n = Length(a);
        ApplyToEach(X, b);
        for i in 0..n-1 {
            Controlled X([a[i]], b[i]);
        }
    }

    /// # Описание
    /// Нахождение индекса элемента с максимальным значением в массиве чисел
    operation IndexOfMax(arr: Int[]) : Int {
        mutable result = 0;
        let n = Length(arr);
        for i in 1..n-1 {
            if(arr[result]<arr[i]) {
                set result = i;
            }
        }
        return result;
    }


    /// # Описание
    /// Многократное повторение измерения смещения распределения для заданных двух массивов
    /// и возврат лучшего (наиболее верооятного) значения смещения
    operation TestArrays(m: Int, x: Int[], y: Int[], tests: Int): (Int, Int[]) {
        let lx = Length(x);
        let ly = Length(y);
        mutable l = lx;
        if(ly>lx) {
            set l = ly;
        }
        let kl = BitSizeI(l-1);
        let km = BitSizeI(m-1);
        let s = 2^km;

        // Аллокируем массив для подсчёта количества исходов экспериментов
        mutable arr = [0, size = s];
        
        for _ in 1..tests {
            use (a, b) = (Qubit[km], Qubit[km]){
                
                U(m, l, x, a);
                U(m, l, y, b);

                ScalarTransformation(a, b);

                let i = Measure(b);
                
                // Добавим полученное значение в счётчик
                set arr w/= i <- arr[i] + 1;
                
                ResetAll(a);
                ResetAll(b);
            }

        }
        return (IndexOfMax(arr), arr);
    }

    /// # Описание
    /// Многократное повторение измерения смещения распределения для двух массивов
    /// которые создаются с помощью заданных генераторов случайных числел
    /// и возврат лучшего (наиболее верооятного) значения смещения

    operation Test(m: Int, randX: (Int)=>Int, randY: (Int)=>Int, l: Int, tests: Int) : (Int, Int[]) {
        mutable x = [0, size = l];
        mutable y = [0, size = l];

        for i in 0..l-1 {
            set x w/= i <- randX(m);
        }
        for i in 0..l-1 {
            set y w/= i <- randY(m);
        }
        
        Message($"x = {x} / y = {y}");

        return TestArrays(m, x, y, tests);
    }
  1. And, actually, the test itself with launch parameters

    @EntryPoint()
    operation Main(n: Int, l: Int, t: Int) : Unit {
        let m=2^n;

        Message("Hello quantum world!");
        Message($"Binom0/Binom0: {Test(m, RandomBinom(0,_), RandomBinom(0,_), l, t)}");
        Message($"Binom0/Binom1: {Test(m, RandomBinom(0,_), RandomBinom(1,_), l, t)}");
        Message($"Binom3/Binom1: {Test(m, RandomBinom(3,_), RandomBinom(1,_), l, t)}");
        Message($"Binom0/Random: {Test(m, RandomBinom(0,_), RandomRandom, l, t)}");
        Message($"Random/Random: {Test(m, RandomRandom, RandomRandom, l, t)}");
        Message($"Random/Binom0: {Test(m, RandomRandom, RandomBinom(0,_), l, t)}");
    }

You can get the full program code here https://github.com/dprotopopov/qcompare

We run and get the result …

PS C:\Projects\qcompare> dotnet run -n 3 -l 32 -t 100
Hello quantum world!
x = [7,6,5,4,2,3,3,3,2,4,4,5,6,3,4,4,5,4,4,4,5,4,3,4,5,4,5,4,6,3,2,3] / y = [2,5,1,2,2,3,6,6,3,2,1,5,5,5,3,3,6,5,2,4,4,6,6,2,4,5,5,4,4,4,3,3]
Binom0/Binom0: (7, [15,16,16,6,4,10,15,18])
x = [4,5,6,3,5,3,2,5,7,6,3,4,5,1,2,1,4,2,4,2,1,4,4,2,2,6,2,3,3,5,4,3] / y = [4,5,4,2,3,2,4,5,0,3,1,2,4,2,4,5,3,2,2,5,2,3,5,2,2,2,2,2,5,7,2,3]
Binom0/Binom1: (1, [14,21,9,18,3,6,16,13])
x = [3,5,7,5,6,1,1,7,1,7,6,1,1,7,7,1,1,0,7,0,6,1,1,6,6,4,6,1,0,0,6,7] / y = [2,5,5,2,2,4,0,5,3,2,5,7,4,5,2,0,5,5,5,0,2,5,2,0,2,3,2,4,4,0,5,4]
Binom3/Binom1: (2, [9,6,24,20,18,10,10,3])
x = [6,4,4,2,2,3,3,4,3,4,2,3,3,4,3,4,4,3,4,2,3,5,4,2,5,4,4,3,3,3,5,5] / y = [6,6,0,2,4,1,2,5,3,1,0,2,4,4,3,3,5,0,1,6,3,2,1,3,2,1,0,0,5,7,4,0]
Binom0/Random: (1, [13,16,11,15,7,9,15,14])
x = [4,4,3,4,4,1,3,7,5,6,1,3,1,4,0,1,5,4,0,0,1,2,4,5,5,0,0,4,2,6,7,3] / y = [4,6,7,5,2,4,5,2,6,6,5,7,0,4,3,0,5,3,0,1,7,3,6,7,6,7,0,7,4,1,5,1]
Random/Random: (3, [11,10,11,21,9,10,18,10])
x = [7,2,3,3,5,7,4,4,3,0,2,0,5,3,0,4,7,6,6,7,7,6,2,4,2,6,7,4,0,0,0,4] / y = [5,3,5,5,4,3,3,4,3,3,2,5,3,3,3,2,4,3,5,3,4,6,2,4,4,3,4,3,4,3,3,5]
Random/Binom0: (7, [10,15,10,15,14,7,12,17])

It seems that everything is correct …

Note. The fact that Binom0/Binom0 coincided with the displacement 7 is also the correct answer, since the binomial distribution is symmetrical, i.e. C(k, 7) = C(n xor 7, 7)

Link

Before

Related posts