System identification and system inversion are two closely related problems. First considering classical, i.e. non-quantum, signals and systems, the basic version of system identification concerns single-input single-output (SISO) systems. It consists of estimating the unknown parameter values of such a system (i.e. of the transform that it performs) belonging to a known class, by using known values of its input (source signal ) and output (signal ). This version  is stated to be “non-blind” by the signal and image processing community 
or “supervised” by the machine learning and data analysis community. The more challenging version of that problem is the blind  or unsupervised one, where the input values are unknown and uncontrolled, but it may be known that the input signal belongs to a given class (due to this partial knowledge, these methods are sometimes stated to be semi-blind). Both versions may then be extended to multiple-input multiple-output (MIMO) systems.
Besides, in various applications, what is needed is not the direct transform performed by the above system, but the inverse of that transform (assuming it is invertible). For SISO non-blind and blind configurations, this is motivated by the fact that one eventually only accesses the ouput of the above direct system, and one aims at deriving a signal which ideally restores the original source signal . To this end, one may first use the above-mentioned system identification methods in order to estimate the direct system, then derive its inverse and eventually transfer the output of the direct system through the inverse system. Alternatively, one may develop methods for initially identifying the inverse system itself. Extended versions of this “(unknown) system inversion” task concern MIMO configurations, where a set of original source signals to are to be respectively restored on the outputs to y of the inverse system.
: as in system inversion, BSS aims at canceling the contributions of all sources but one in each output signal of the separating system; however, in BSS, one often allows each output signal to be equal to a source signal only up to an acceptable residual transform. These transforms, called indeterminacies, cannot be avoided because only limited constraints are set on the source signals and on the direct system which combines (i.e., “mixes”, in BSS terms) these signals. In particular, the first class of BSS methods that was developed and that is still of major importance is Independent Component Analysis, or ICA,,
, which may be seen as an extension of more conventional Principal Component Analysis, or PCA (PCA alone cannot achieve BSS ). ICA is a statistical approach, which essentially requires statistically independent random source signals. Thus, for the simplest class of mixtures, ICA is guaranteed to restore the source signals up to limited indeterminacies ,,.
We now consider quantum information processing (QIP)  and quantum machine learning ,,, i.e. processing of quantum data and/or processing with quantum means, and we still focus on system identification and inversion problems. Among these problems, the one which was first studied is the quantum version of non-blind system identification, especially  introduced in 1997 in  and called “quantum process tomography” or QPT by the QIP community: see e.g. , , , , , , , , , . The quantum version of the above-mentioned classical source separation, called Quantum Source Separation, or QSS, and especially its blind version, or BQSS, were then introduced in 2007 in 
. Two main classes of BQSS methods were developed since then. The first one may be seen as a quantum extension of the above-mentioned classical ICA methods, since it takes advantage of the statistical independence of the parameters that define random source quantum states (qubit states). It is called Quantum Independent Component Analysis (or QICA, see e.g.,) or, more precisely, Quantum-Source Independent Component Analysis (or QSICA, see e.g. ) to insist on the quantum nature of the considered source data, whereas it uses classical processing means (after quantum/classical data conversion). The second main class of BQSS methods was introduced in 2013-2014 in , and then especially detailed in . It is based on the unentanglement of the considered source quantum states and it typically uses quantum processing means to restore these unknown states from their coupled version. Independently from the above quantum extensions of ICA, a quantum version of PCA was introduced in 2014 in .
In the present paper, our first contribution (see Section IV) concerns yet another type of unsupervised quantum machine learning methods , namely Blind Quantum Process Tomography, or BQPT. Here again, the term “blind”, or “unsupervised”, refers to the fact that we consider situations where the input values of the process to be identified are unknown, but they are requested to meet some (hereafter statistical) properties. We briefly introduced that BQPT concept in 2015 in  and we only outlined some resulting BQPT methods in that and some subsequent short conference papers, but only as spin-offs of corresponding BQSS methods. On the contrary, the present paper is the first one where we provide a detailed description of a method which combines the following features:
This method is primarily intended for BQPT, not for BQSS. To this end, it only uses classical processing means: on the contrary, using quantum processing means requires one to precisely characterize them beforehand, which is a significant drawback here, since BQPT, as QPT, is especially developed as a tool for characterizing quantum gates, as discussed in Sections V.4 and VII.
This BQPT method therefore first performs measurements at the output of the system, i.e. quantum process, to be identified (see Section III), in order to convert quantum states into classical-form data before they are processed with classical means.
Moreover, for the Heisenberg coupling process considered below as an example (see Section II), we aim at minimizing the number of types of measurements performed to fully characterize that process.
As detailed further in this paper, usual, i.e. non-blind, QPT, as well as the above first form of BQPT, use sample frequencies of the above-mentioned measurement outcomes at the output of the process (i.e. normalized cumulative values associated with these outcomes), derived from many copies of each considered state value. To apply such methods, one should therefore be able to prepare many copies of the same input state, which is cumbersome. Our second contribution in this paper (see Section V) then consists of extensions of the above BQPT methods, which also allow one to use few copies or even one instance (i.e. preparation) of each quantum state. The numerical performance of this second type of methods is reported in Section VI, whereas its applications are presented in Section VII, together with conclusions drawn from this investigation.
Ii Considered quantum process and state properties
In QPT and BQPT problems, the system or process to be identified receives a quantum state , e.g. associated with a set of qubits considered at time . This system outputs a quantum state associated with the same set of qubits at a later time . The behavior of the system is defined by an operator which may be rather general (e.g. an arbitrary unitary operator) or which may be restricted to a given class of operators corresponding to a specific type of physical devices (see e.g. ). In this paper, we address the second case and we consider a device composed of two distinguishable qubits  implemented as electron spins 1/2, that are coupled according to the cylindrical-symmetry Heisenberg model, which is e.g. relevant for spintronics applications. We stress that this type of coupling is only used as a concrete example, to show how to fully implement the proposed concepts in a relevant case, but that these concepts and resulting practical BQPT algorithms may then be transposed by the reader to other classes of quantum processes and associated applications.
The symmetry axis of the Heisenberg model is here denoted as . The considered spins are supposed to be placed in a magnetic field (also oriented along and with a magnitude ) and thus coupled to it. Moreover, we assume an isotropic tensor, with principal value . The time interval when these spins are considered is supposed to be short enough for their coupling with their environment to be negligible. In these conditions, the temporal evolution of the state of the device composed of these two spins is governed by the following Hamiltonian:
, where is the Bohr magneton, i.e. and is the reduced Planck constant,
are the three components of the vector operatorassociated with spin in a cartesian frame,
and are the principal values of the exchange tensor.
Among the above parameters, the value of may be experimentally determined, and can be measured. The values of and are here assumed to be unknown.
We here suppose that each spin , with is prepared, i.e. initialized, at a given time , in the pure state
where and are eigenkets of
, for the eigenvaluesand respectively. We will further use the polar representation of the qubit parameters and , which reads
where is the imaginary unit, and with and
because each spin state has unit norm. Moreover, for each couple of phase parameters and , only their difference has a physical meaning. After they have been prepared, these spins are coupled according to the above-defined model for .
Hereafter, we consider the state of the overall system composed of these two distinguishable spins. At time , this state is equal to the tensor product of the states of both spins defined in (2). It therefore reads
in the four-dimensional basis .
The state of this two-spin system then evolves with time. Its value at any subsequent time may be derived from its above-defined Hamiltonian. It is defined  by
where and are the column vectors of components of and , respectively, in basis . For instance, as shown by (6),
where stands for transpose. Moreover, the matrix of (7), which defines the transform applied to , reads
and equal to
The four real (angular) frequencies to in (LABEL:eq-opmixdiagdef) depend on the physical setup. In , it was shown that they read
In this paper, we address the (B)QPT problem, i.e. we aim at estimating the matrix involved in (7), which defines the considered quantum process. Moreover, we estimate it in a blind, i.e. unsupervised, way, that is:
by using values of the output state of this process,
without using nor knowing values of its input state ,
but by knowing and exploiting some properties of these states . In this paper, these requested properties are as follows. The states are required to be unentangled (as shown by (6)). Besides, the proposed BQPT methods are statistical approaches and the six parameters , and , with , defined in (3) are constrained to have properties that are similar to those requested in the above-mentioned QSICA methods: (i) these parameters are random valued, so that we here consider random pure quantum states (see  for more details) and (ii) some combinations of the random variables (RVs) , and are statistically independent and have a few known statistical features, as detailed further in this paper.
As explained in Section I, the considered BQPT task is performed by using only classical-form processing means. To this end, the available quantum-form data, namely the output states , are first converted into classical-form data, by means of measurements, as described hereafter.
Iii Measurements for process outputs
The first type of proposed BQPT approaches uses a set of copies of each output state . For each copy, it measures the components of the considered two spins along the above-defined direction. The result of each such measurement has four possible values, that is , , or in normalized units (see Appendix A.1
). Their probabilities are respectively denoted asto hereafter. Using the polar representation (3), these probabilities read ,
Probability is not considered hereafter because the sum of to is equal to 1.
In practice, for each value of state , estimates of probabilities to are derived, typically as the sample frequencies of the associated measurement outcomes obtained for all copies of (see e.g. ,,,).
Similarly, these BQPT approaches use another set of copies of each output state , by measuring the two spin components along an axis which is orthogonal to . These measurements yield the same four possible outcomes as above, but with different probabilities, which are denoted as to hereafter. As shown in , these probabilities have the following properties
The value of in (25) is known, since it can be derived from the above-defined known quantities. Moreover
The BQPT methods proposed in this paper therefore consist of two major steps. The first step aims at estimating the unknown values of the parameters , and of the mappings from the parameters , and of the initial qubit states to the probabilities of measurement outcomes, namely and , with to 4, or their combinations. The second step then uses the estimated values of , and to derive an estimate of matrix of (9) and hence of the complete matrix of (9), which defines the considered process. We now proceed to the description of these methods.
Iv Multiple-preparation BQPT methods
We first consider the estimation of parameter . This is achieved by exploiting (15). If we were developing a conventional, i.e. non-blind, QPT method, we would use one or several instances of Eq. (15), and each of these instances would involve (i) known values of the input parameters , , and hence of and (ii) an estimate of the set of output probabilities , derived from a set of copies of . This approach is constraining because it requires one to precisely prepare each state value (for state preparation and associated errors, see e.g. ), otherwise the errors in and yield errors in the estimated value of .
The above drawback is avoided by our BQPT methods. Following the above-defined terminology, these methods are blind in the sense that they estimate by using only a set of estimated values of output probabilities , without knowing the values of the input parameters , , and hence , but requesting them to have some known properties. More precisely, we here consider statistical methods, which operate with a set of random states and which thus only set constraints on some of the statistical parameters of (combinations of) , , , not on their individual values for each state . In particular, the versions of these BQPT methods considered in this paper use only the first-order mean statistics of the available quantities , i.e. their expectations . When assuming , and to be statistically independent RVs, (15) yields
In this equation, is known: in practice, it is estimated as the sample mean of the estimates of all values of , themselves typically estimated with sample frequencies, as explained above. Besides, as detailed e.g. in , for BQSS methods intended for the Heisenberg coupling model, setting the constraint
with and . Taking the sample mean of any function of defined by (LABEL:eq-qubitindexstdstateplusmodulus-sol-both-from-mlsp) then yields estimates of all statistics of involved in (28). Finally, we only set the following constraint on one statistical parameter of the used values of and , again without having to know their individual values. We request the states to be prepared with a procedure which is such that the value of is known. With these constraints on input state statistics, the only unknown in (28) is . By solving this type of equations, this BQPT method then yields the desired estimate of . In particular, a simple case consists of using (which may e.g. be achieved by preparing the two spins with states such that and are statisticially independent and have the same statistics): then, (28) straightforwardly yields
In some configurations the sign of is known ,, so that the value of may be derived from (31). Otherwise, it may be derived from another instance of (28), using data that yield another value of : details about how to solve this sign indeterminacy and how to also estimate parameters and are provided below for an improved version of our methods. Indeed, the above version of BQPT is attractive because it does not require each value of to be known, but it still yields a limitation: it requires one to be able to prepare the same value a large number of times, to derive an associated frequency-based estimate of each set of probabilities . This still requires some control of the input states of the process, that we would like to avoid, in order to simplify the practical operation of BQPT methods and to make them “blinder”. We hereafter show how to avoid this preparation of many copies of each state .
V Single-preparation BQPT methods
v.1 Single-preparation QIP
As a second contribution in this paper, we now extend (B)QPT methods so that they can operate with a few copies or even a single instance of each considered input state . For non-blind methods as defined above, this does not seem to be possible, because they need many copies of each state and associated outcomes of measurements performed for each state , in order to derive a frequency-based estimate of each set of probabilities . On the contrary, our blind versions of QPT can be extended so as to reach this goal, because they only need one to estimate expectations of these (now random) probabilities , i.e. , not each of their individual values for each state . In the short conference paper , we very recently introduced a general QIP framework (i.e., not restricted to BQPT) for estimating expectations of probabilities of outcomes of general types of quantum measurements. Its principle is summarized hereafter, whereas its detailed description and properties are provided in Appendix A.
For each expectation of a random probability to be estimated, as discussed above, in practice the expectation operator is replaced by a sample mean, i.e. by a sum (of values, moreover normalized). Similarly, each probability is replaced by a sample frequency, i.e. by a sum (of 1 and 0, depending whether the considered event occurs or not for each trial defined by a preparation of the initial quantum states (2) and by an associated measurement of a couple of spin components; this summation is here again followed by a normalization, by the total number of trials). is therefore estimated by a (normalized) “sum of sums”, which may then be reinterpreted as a single global sum, and what primarily matters is the total number of preparations of initial quantum states (2) involved in that global sum, whereas the number of preparations for each state value (2) may be decreased, down to 1, as confirmed by simulations in Section VI (which also justifies why even better performance is obtained when decreasing the number of preparations per state for a given total number of preparations, i.e. while increasing accordingly the considered number of different states). The corresponding BQPT methods are therefore called single-preparation BQPT methods. It should be clear that they can be freely used with either one instance or several (e.g. many) copies per state, i.e. the above terminology means that these methods allow one to use a single instance of each state. On the contrary, our so-called multiple-preparation BQPT methods force one to use many state copies to achieve good performance.
v.2 Estimating the parameter of measurements
We here aim at using the single-preparation approach of Section V.1 to estimate the parameter involved in the probabilities . We hereafter again take advantage of (28), and especially of its version (31), derived in Section IV for our multiple-preparation BQPT method. However, these expressions involve and which, unlike in Section IV, cannot here be estimated by using the expectation of the square of (LABEL:eq-qubitindexstdstateplusmodulus-sol-both-from-mlsp) because this involves expectations of nonlinear combinations of the above probabilities to , whereas we here aim at developing a single-preparation algorithm, for which Section V.1 only defined how to estimate the expectations of to themselves. We here solve this problem by using a modified approach, where we first take the expectation of (14) and (16), again for statistically independent RVs and . This yields
These equations involve only the unknown of interest, and . Again setting the constraint (29), they yield the unique solution
again with and . If the sign of is known, the value of may thus be derived from (31), therefore using data such that . Otherwise, (31) is first used to estimate , which yields , and the sign of is then derived from another set of spin state preparations, now considering the case when . Eq. (28) then yields
Taking the sign of this equation, where a factor is guaranteed to be positive, results in
For this second set of spin state preparations, (i) we do not request the value of but only its sign to be known, (ii) the values of , and may again be estimated as explained above. Also using the above estimate of , Eq. (36) then allows one to estimate the sign of .
v.3 Estimating the parameters of measurements
We then show how to estimate the parameters and of (20), using measurements of spin components along the axis, in addition to the axis, and the corresponding expectations and . Here again, we only constrain the statistical parameters of the RVs and , not their individual deterministic values, in order to be able to solve (20) with respect to and . More precisely, the RVs , , and are here statistically independent. Besides, and have the same statistics. Finally, and have the same statistics, moreover with
which is e.g. obtained with RVs
whose probability density functions are even and non-zero on. In that case, (14), (20)-(22) and (LABEL:eq-twoqubitsprobaplusplusdirxx-plus-twoqubitsprobaminusminusdirxx-explicit) yield  (with the same statistics for and 2):
Once , and have been estimated as explained above, Eq. (39), with due to , yields
Using (25), (44) and (45), Eq. (41) and (42) then yield estimates of and . The only unknowns of (40) are then and . One could try and solve a single equation (40), by taking into account that and are the cosine and sine of the same angle (see (23)-(24)). However, the solutions of such an equation yield a problematic indeterminacy. This problem is avoided by creating two linearly independent equations (40), by using two sets of statistics for , , and . Solving these two equations yields and .
v.4 Estimating the quantum process
We finally show how the estimates of the parameters , and obtained above may be used to estimate the matrix of (LABEL:eq-opmixdiagdef) and hence the complete matrix of (9), which defines the considered process in the standard basis.
where is one determination associated with the actual value , i.e. is equal to up to the additive constant , where is an integer.
Eq. (46)-(49) define the expressions of the scaled actual principal values and with respect to the actual values of , and . The latter values are unknown but, in practice, the procedure defined in Section V.2 yields an estimate of the value of (for the considered value ). From this, one may derive an estimate of by using in (47). One would then like to derive an estimate of from (46). But one does not know the actual value involved in (46) in the fully blind case considered here, i.e. when no prior information is available about the value of (as opposed to the case when one at least knows in which range of values is situated, which defines the minimum and maximum possible values of ). In this blind method, one can then only select an arbitrary integer and derive the corresponding scaled “shifted estimate” of by using
Similarly, the procedure defined in Section V.3 yields estimates and of the values of and (for the considered value ). From this, one first derives an estimate of by using and in (49). Then, based on (48), one derives a scaled shifted estimate of by using
The shifted estimates and provided by this method are therefore equal to the quantities of interest, that is and , only up to (the above neglected estimation errors and) additive constants which are integer multiples of . These constants are the “undeterminacies” of this method in the classical BSS sense, i.e. the undesired remaining differences between the above estimated and actual quantities, from the point of view of the quantities and . They then yield the following indeterminacies from the point of view of the matrix of the considered quantum process, which is eventually to be estimated. Using the above estimates and , one derives the associated estimate of the matrix (i) by inserting these estimates, which may be expressed as (52)-(53), into (LABEL:eq-opmixdiagdef)-(13), which yields the corresponding estimate of , and (ii) finally by using (9) and (10) to derive the associated estimate of . These calculations especially yield (taking into account that and )
The estimate provided by this first method is therefore equal to the actual matrix up to the phase factor . More specifically, this factor is equal to 1 and thus diseappears for part of the possible values of the integers and , e.g. when is a multiple of 2 and is a multiple of 4. This yields the same phenomenon for . The general phase factor cannot be avoided with this method if no additional information is available. It is the only and quite weak indeterminacy of this BQPT method from the point of view of and . Moreover, we hereafter introduce an extended version of that method, which completely removes this indeterminacy by taking the typical applications of (B)QPT methods into account.
As discussed e.g. in , , , , , , QPT (and hence our blind extension) may especially be used as a tool for characterizing quantum gates, which are the building blocks of quantum computers. This characterization is typically performed before using the considered gates for quantum computation, thus leading to a two-phase approach, composed of an “identification phase” and then of a “computation phase”, for these quantum processes/gates. Moreover, one may consider scenarios where these processes/gates are used in coherent but somewhat different conditions during the identification and computation phases. We hereafter propose such an approach for extending the above BQPT method so as to remove its indeterminacy. We do not claim that the Heisenberg coupling model considered in this paper could be used as a suitable process/gate for quantum computers: it is just used as an example hereafter, to illustrate a possible procedure for removing BQPT indeterminacies, thus then allowing the reader to extend this procedure to other processes/gates that could be of interest in other configurations.
The approach that we propose uses three values of the time interval involved in (LABEL:eq-opmixdiagdef), that we hereafter denote as , and . The first step of the identification phase uses the time interval , essentially to obtain an estimate of associated with this value , that we hereafter denote as for the sake of clarity. More precisely, this first step of the identification phase derives the shifted estimate in the same way as in the above first BQPT method, i.e. using (50), with here replaced by . Therefore, when neglecting estimation errors, this again yields (52), but with our modified notations, that is
The second step of the identification phase then uses the time interval , with ( may instead be set to any other even multiple of , but we keep the values of , and as close as possible to one another, in order to minimize the differences in the conditions of operation in the two steps of the identification phase and in the computation phase). This second step of the identification phase essentially aims at obtaining an estimate of associated with the value , that we therefore hereafter denote as . More precisely, this second step derives the shifted estimate in the same way as in the above first BQPT method, except that this step is here performed with , so that it uses (51) with replaced by , moreover taking into account that the term of this modified version of (51) is here obtained as being equal to the value of this second BQPT method multiplied by 2. When neglecting estimation errors, taking the difference between the modified versions of (48) and (51), and using (57), Eq. (53) is thus replaced by
The computation phase then involves the same type of quantum process (9)-(LABEL:eq-opmixdiagdef), but with a time interval , between input state preparation at time and output state use at time , which is set to , with (again, may instead be set to any other even multiple of , but we keep the values of , and as close as possible to one another). This computation phase should then be analyzed as follows. During that phase, the considered actual process is defined by (9)-(13), but with replaced by . From the point of view of that computation phase, the estimate of that actual process is obtained by replacing and by and in (9)-(13), the latter estimates being derived as explained above by our extended BQPT method (up to the factors and ). These estimates have the properties defined by (57) and (58). Since and , this yields