Learning with density matrices and random features

A density matrix describes the statistical state of a quantum system. It is a powerful formalism to represent both the quantum and classical uncertainty of quantum systems and to express different statistical operations such as measurement, system combination and expectations as linear algebra operations. This paper explores how density matrices can be used as a building block for machine learning models exploiting their ability to straightforwardly combine linear algebra and probability. One of the main results of the paper is to show that density matrices coupled with random Fourier features could approximate arbitrary probability distributions over Rn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbb {R}^{n}$$\end{document}. Based on this finding the paper builds different models for density estimation, classification and regression. These models are differentiable, so it is possible to integrate them with other differentiable components, such as deep learning architectures and to learn their parameters using gradient-based optimization. In addition, the paper presents optimization-less training strategies based on estimation and model averaging. The models are evaluated in benchmark tasks and the results are reported and discussed.


Introduction
The formalism of density operators and density matrices was developed by von Neumann as a foundation of quantum statistical mechanics (Von Neumann, 1927).From the point of view of machine learning, density matrices have an interesting feature: the fact that they combine linear algebra and probability, two of the pillars of machine learning, in a very particular but powerful way.
The main question addressed by this work is how density matrices can be used in machine learning models.One of the main approaches to machine learning is to address the problem of learning as one of estimating a probability distribution from data: joint probabilities P (x, y) in generative supervised models or conditional probabilities P (y|x) in discriminative models.
The central idea of this work is to use density matrices to represent these probability distributions tackling the important question of how to encode arbitrary probability density functions in R n into density matrices.
The quantum probabilistic formalism of von Neumann is based on linear algebra, in contrast with classical probability which is based on set theory.In the quantum formalism the sample space corresponds to a Hilbert space H and the event space to a set of linear operators in H, the density operators (Wilce, 2021).
The quantum formalism generalizes classical probability.A density matrix in an n-dimensional Hilbert space can be seen as a catalog of categorical distributions on the finite set {1 . . .n}.A direct application of this fact is not very useful as we want to efficiently model continuous probability distributions in R n .One of the main results of this paper is to show that it is possible to model arbitrary probability distributions in R n using density matrices of finite dimension in conjunction with random Fourier features (Rahimi and Recht, 2007).In particular the paper presents a method for non-parametric density estimation that combines density matrices and random Fourier features to efficiently learn a probability density function from data and to efficiently predict the density of new samples.
The fact that the probability density function is represented in matrix form and that the density of a sample is calculated by linear algebra operations makes it easy to implement the model in GPU-accelerated machine learning frameworks.This also facilitates using density matrices as a building block for classification and regression models, which can be trained using gradientbased optimization and can be easily integrated with conventional deep neural networks.The paper presents examples of these models and shows how they can be trained using gradient-based optimization as well as optimization-less learning based on estimation.
The paper is organized as follows: Section 2 covers the background on kernel density estimation, random features, and density matrices; Section 5 presents four different methods for density estimation, classification and regression; Section 6 discusses some relevant works; Section 7 presents the experimental evaluation; finally, Section 8 discusses the conclusions of the work.
2 Background and preliminaries

Kernel density estimation
Kernel Density Estimation (KDE) (Rosenblatt, 1956;Parzen, 1962), also known as Parzen-Rossenblat window method, is a non-parametric density estimation method.This method does not make any particular assumption about the underlying probability density function.Given an iid set of samples X = {x 1 , . . ., x N }, the smooth Parzen's window estimate has the form where k λ (•) is a kernel function, λ is the smoothing bandwidth parameter of the estimate and M λ is a normalizing constant.A small λ-parameter implies a small grade of smoothing.Rosenblatt (1956) and Parzen (1962) showed that eq. ( 1) is an unbiased estimator of the pdf f .If k γ is the Gaussian kernel, eq. ( 1) takes the form where M γ = (π/γ) d 2 .KDE has several applications: to estimate the underlying probability density function, to estimate confidence intervals and confidence bands (Efron, 1992;Chernozhukov et al, 2014), to find local modes for geometric feature estimation (Chazal et al, 2017;Chen et al, 2016), to estimate ridge of the density function (Genovese et al, 2014), to build cluster trees (Balakrishnan et al, 2013), to estimate the cumulative distribution function (Nadaraya, 1964), to estimate receiver operating characteristic (ROC) curves (McNeil and Hanley, 1984), among others.
One of the main drawbacks of KDE is that it is a memory-based method, i.e. it requires the whole training set to do a prediction, which is linear on the training set size.This drawback is typically alleviated by methods that use data structures that support efficient nearest-neighbor queries.This approach still requires to store the whole training dataset.

Random features
Random Fourier features (RFF) (Rahimi and Recht, 2007) is a method that builds an embedding . One of the main applications of RFF is to speedup kernel methods, being data independence one of its advantages.
The RFF method is based on the Bochner's theorem.In layman's terms, Bochner's theorem shows that a shift invariant positive-definite kernel k(•) is the Fourier transform of a probability measure p(w).Rahimi and Recht (2007) use this result to approximate the kernel function by designing a sample procedure that estimates the integral of the Fourier transform.The first step is to draw D iid samples {w 1 , . . .w D } from p and D iid samples {b 1 , . . .b D } from a uniform distribution in [0, 2π].Then, define: (3) Rahimi and Recht (2007) showed that the expected value of ϕ T rff (x)ϕ rff (y) uniformly converges to k(x, y): Theorem 1 (Rahimi and Recht, 2007) Let M be a compact subset of R d with a diameter diam(M).Then for the mapping ϕ rff defined above, we have

Pr sup
x,y∈M where, σ 2 p is the second momentum of the Fourier transform of k.In particular, for the Gaussian kernel σ 2 p = 2dγ, where γ is the spread parameter (see Eq. 2).
Different approaches to compute random features for kernel approximation have been proposed based on different strategies: Monte Carlo sampling (Le et al, 2013;Yu et al, 2016), quasi-Monte-Carlo sampling (Avron et al, 2016;Shen et al, 2017), and quadrature rules (Dao et al, 2017).
RFF may be used to formulate a non-memory based version of KDE.For the Gaussian kernel we have: Φ train in eq. ( 5) can be efficiently calculated during training time, since is just an average of the RFF embeddings of the training samples.The time complexity of prediction, eq. ( 5), is constant on the size of the training dataset.The price of this efficiency improvement is a loss in precision, since we are using an approximation of the Gaussian kernel.

Density estimation with density matrices
The Gaussian kernel satisfy ∀x, y ∈ R d , k γ (x, y) > 0, however the RFF estimation may be negative.To alleviate this we could estimate the square of the kernel and use the fact that k γ (x, y) = k 2 γ/2 (x, y).In this case we have: In eq. ( 6) it is important to take into account that the parameters of the RFF embedding, ϕ rff , are sampled using a parameter γ/2 for the Gaussian kernel.
Proposition 2 Let M be a compact subset of R d with a diameter diam(M), let X = {x i } i=1...N ⊂ M a set of iid samples, then fρtrain (eq.( 6)) and fγ satisfy: The Parzen's estimator is an unbiased estimator of the true density function from which the samples were generated and Proposition 2 shows that fρtrain (x) can approximate this estimator.
A further improvement to the fρtrain (x) estimator is to normalize the RFF embedding as follows: Here we use the Dirac notation to emphasize the fact that φrff is a quantum feature map.This has the effect that the estimation k γ (x, x) = φrff (x) φrff (x) = 1 will be exact and ∀x, y ∈ R d , φrff (x) φrff (y) ≤ 1.
During the training phase ρ train is estimated as the average of the cross product of the normalized RFF embeddings of the training samples: The time complexity of calculating ρ train is O(D 2 N ), i.e. linear on the size of the training dataset.One advantage over conventional KDE is that we do not need to store the whole training dataset, but a more compact representation.
During the prediction phase the density of a new sample is calculated as: The fρtrain estimator has an important advantage over the Parzen's estimator, its computational complexity.The time to calculate the Parzen's estimator (eq.( 2)) is O(dN ) while the time to estimate the density based on the density matrix ρ train (eq.( 10)) is O(D 2 ), which is constant on the size of the training dataset.
The ρ train matrix in eq. ( 9) is a well known mathematical object in quantum mechanics, a density matrix, and eq.( 10) is an instance of the Born rule which calculates the probability that a measurement of a quantum system produces a particular result.This connection and the basic ideas behind density matrices are discussed in the next section.

Density matrices
This section introduces some basic mathematical concepts that are part of the mathematical framework that supports quantum mechanics and discusses their connection with the ideas introduced in the previous subsection.The contents of this section are not necessary for understanding the rest of the paper and are included to better explain the connection of the ideas presented in this paper with the quantum mechanics mathematical framework.
The state of a quantum system is represented by a vector ψ ∈ H, where H is the Hilbert space of the possible states of the system.Usually1 H = C d .
As an example, consider a system that could be in two possible states, e.g. the spin of an electron that could be up (↑) or down (↓) with respect to some axis z.The state of this system is, in general, represented by a regular column vector |ψ⟩ = (α, β), with |α| 2 + |β| 2 = 1.This state represents a system that is in a superposition of the two basis states |ψ⟩ = α ↑ +β ↓.The outcome of a measurement of this system, along the z axis, is determined by the Born rule: the spin is up with probability |α| 2 and down with probability |β| 2 .Notice that α and β could be negative or complex numbers, but the Born rule guarantees that we get valid probabilities.
The normalized RFF mapping (eq.( 8)) can be seen as a function that maps a sample to the state of a quantum system.In quantum machine learning literature, there are different approaches to encode data in quantum states (Schuld, 2018).The use of RFF as a data quantum encoding strategy was first proposed by (González et al, 2020;González et al, 2021).
The probabilities that arise from the superposition of states in the previous example is a manifestation of the uncertainty that is inherent to the nature of quantum physical systems.We call this kind of uncertainty quantum uncertainty.Other kind of uncertainty comes, for instance, from errors in the measurement or state-preparation processes, we call this uncertainty classical uncertainty.A density matrix is a formalism that allows us to represent both types of uncertainty.To illustrate it, let's go back to our previous example.The density matrix representing the state ψ is: As a concrete example, consider 2 the corresponding density matrix is: which represents a superposition state where we have a 1 2 probability of measuring any of the two states.Notice that the probabilities for each state are in the diagonal of the density matrix.ρ 1 is a rank-1 density matrix, and this means that it represents a pure state.A mixed state, i.e. a state with classical uncertainty, is represented by a density matrix with the form: where , and {ψ i } i=1...N are the states of a an ensemble of N quantum systems, where each system has an associated probability p i .The matrix ρ train in eq. ( 9) is in fact a density matrix that represents the state of an ensemble of quantum systems where each system corresponds to a training data sample.The probability is the same for all the N elements of the ensemble, 1 N .As a concrete example of a mixed state consider two pure states ψ 2 = (1, 0) and ψ ′ 2 = (0, 1), and consider a system that is prepared in state ψ 2 with probability 1 2 and in state ψ ′ 2 with probability 1 2 as well.The state of this system is represented by the following density matrix: At first sight, states ρ 1 and ρ 2 may be seen as representing the same quantum system, one where the probability of measuring an up state (or down state) in the z axis is 1 2 .However they are different systems, ρ 1 represents a system with only quantum uncertainty, while ρ 2 corresponds a system with classical uncertainty.To better observe the differences of the two systems we have to perform a measurement along a particular axis.To do so, we use the following version of the Born rule for density matrices: which calculates the probability of measuring the state φ in a system in state ρ.

Density matrix kernel density estimation (DMKDE)
In this subsection we present a model for non-parametric density estimation based on the ideas discussed in subsection 3. The model receives an input x ∈ R d , represents it using a RFF quantum feature map (eq.( 3)) and estimates the density of it using eq.( 10).The model can be trained by averaging the density matrices corresponding to the training samples or by using stochastic gradient descent.The second approach requires a re-parametrization of the model that we discuss next.
The main parameter of the model is ρ train , which is a Hermitian matrix.To ensure this property, we can represent it using a factorization as follows: 16) where V ∈ R r×D , Λ ∈ R r×r is a diagonal matrix and r < D is the reduced rank of the factorization.With this new representation, eq. ( 10) can be re-expressed as: This reduces the time to calculate the density of a new sample to O(Dr).
x The model is depicted in Fig. 1 and its function is governed by the following equations: The hyperparameters of the model are the dimension of the RFF representation D, the spread parameter γ of the Gaussian kernel and the rank r of the density matrix factorization.The parameters are the weights and biases of the RFF, W rff ∈ R D×d and b rff ∈ R d (corresponding to the w i and b i parameters in Eq. 3), and the components of the factorization, V ∈ R r×D and λ ∈ R r , the vector with the elements in the diagonal of Λ.
The training process of the model is as follows: using the random Fourier features method described in Section 2.2 for approximating a Gaussian kernel with parameters γ/2 and D. 3. Apply φrff (eq.( 8)): 4. Estimate ρ train : 5. Make a spectral decomposition of rank r of ρ train : Notice that this training procedure does not require any kind of iterative optimization.The training samples are only used once and the time complexity of the algorithm is linear on the number of training samples.The complexity of step 4 is O(D 2 N ) and of step 5 is O(D 3 ).
Since the operations defined in eq. ( 18) are differentiable, it is possible to use gradient-descent to minimize an appropriate loss function.For instance, we can minimize the negative log-likelihood: In contrast with the learning procedure based on density matrix estimation, using SGD does not guarantee that we will approximate the real density function.If we train all the parameters, maximizing the likelihood becomes an ill-posed problem because of singularities (a Gaussian with arbitrary small variance centered in one training point) (Bishop, 2006).Keeping fixed the RFF parameters and optimizing the parameters of the density matrix, V and λ has shown a good experimental performance.The version of the model trained with gradient descent is called DMKDE-SGD.Something interesting to notice is that the process described by eqs.( 19) and ( 20) generalizes density estimation for variables with a categorical distribution, i.e. x ∈ {1, . . ., K}.To see this, we replace φrff in eq. ( 19) by the well-known one-hot-encoding feature map: where E i is the unit vector with a 1 in position i and 0 in the other positions.It is not difficult to see that in this case

Density matrix kernel density classification (DMKDC)
The extension of kernel density estimation to classification is called kernel density classification (Hastie et al, 2009).The posterior probability is calculated as where π j and fj are respectively the class prior and the density estimator of class j.We follow this approach to define a classification model that uses the density estimation strategy based on RFF and density matrices described in the previous section.The input to the model is a vector x ∈ R d .The model is depicted in Fig. 2 and defined by the following equations: The hyperparameters of the model are the dimension of the RFF representation D, the spread parameter γ of the Gaussian kernel, the class priors π i and the rank r of the density matrix factorization.The parameters are the weights and biases of the RFF, W rff ∈ R D×d and b rff ∈ R d , and the components of the factorization, V i ∈ R r×D and λ i ∈ R for i = 1 . . .K.
The model can be trained using two different strategies: one, using DMKDE to estimate the density matrices of each class; two, use stochastic gradient descent over the parameters to minimize an appropriate loss function.
The training process based on density matrix estimation is as follows: 1. Use the RFF method to calculate W rff and b rff .2. For each class i: (a) Estimate π i as the relative frequency of the class i in the dataset.(b) Estimate ρ i using eq.( 20) and the training samples from class i.
(c) Find a factorization of rank r of ρ i : Notice that this training procedure does not require any kind of iterative optimization.The training samples are only used once and the time complexity of the algorithm is linear on the number of training samples.The complexity of step 2(b) is O(D 2 N ) and of 2(c) is O(D 3 ).
Since the operations defined in eqs.(25a) to (25d) are differentiable, it is possible to use gradient-descent to minimize an appropriate loss function.For instance, we can use categorical cross entropy: where y = (y 1 , . . ., y K ) corresponds to the one-hot-encoding of the real label of the sample x.The version of the model trained with gradient descent is called DMKDC-SGD.An advantage of this approach is that the model can be jointly trained with other differentiable architecture such as a deep learning feature extractor.

Quantum measurement classification (QMC)
In DMKDC we assume a categorical distribution for the output variable.If we want a more general probability distribution we need to define a more general classification model.The idea is to model the joint probability of inputs and outputs using a density matrix.This density matrix represents the state of a bipartite system whose representation space is H X ⊗ H Y where H X is the representation space of the inputs, H Y is the representation space of the outputs and ⊗ is the tensor product.A prediction is made by performing a measurement with an operator specifically prepared from a new input sample.This model is based on the one described by González et al (2020) and is depicted in Figure 3 and works as follows: • Measurement operator.The effect of this measurement operator is to collapse, using a projector to z, the part H X of the bipartite system while keeping the H Y part unmodified.This is done by defining the following operator: where Id H Y is the identity operator in H Y .• Apply the measurement operator to the training density matrix: • Calculate the partial trace of ρ with respect to X to obtain a density matrix that encodes the prediction: The parameter of the model, without taking into account the parameters of the feature maps, is the ρ train ∈ R D X D Y ×D X D Y density matrix, where D X and D Y are the dimensions of H X and H Y respectively.As discussed in Section 5.1, the density matrix ρ train can be factorized as: where V ∈ R r×D X D Y , Λ ∈ R r×r is a diagonal matrix and r < D X D Y is the reduced rank of the factorization.This factorization not only helps to reduce the space necessary to store the parameters, but learning V and Λ, instead of ρ train , helps to guarantee that ρ train is a valid density matrix.
As in Subsection 5.2, we described two different approaches to train the system: one based on estimation of the ρ train and one based on learning ρ train using gradient descent.QMC can be also trained using these two strategies.
In the estimation strategy, given a training data set {(x i , y i )} i=1...N the training density matrix is calculated by: The computational cost is O(N D2 X D 2 Y ).For the gradient-descent-based strategy (QMC-SGD) we can minimize the following loss function: where ρ Yii is the i-th diagonal element of ρ Y .
As in DMKDC-SGD, this model can be combined with a deep learning architecture and the parameters can be jointly learned using gradient descent.
QMC can be used with different feature maps for inputs and outputs.For instance, if ϕ X = ϕ rff (eq.( 3)) and ϕ Y = ϕ ohe (eq.( 22)), QMC corresponds to DMKDC.However, in this case DMKDC is preferred because of its reduced computational cost.

Quantum measurement regression (QMR)
In this section we show how to use QMC to perform regression.For this we will use a feature map that allows us to encode continuous values.The feature map is defined with the help of D equally distributed landmarks in the [0, 1] interval 2 : The following function (which is equivalent to a softmax) defines a set of unimodal probability density functions centered at each landmark: where β controls the shape of the density functions.
The feature map is defined as: This feature map is used in QMC as the feature map of the output variable (ϕ Y ).To calculate the prediction for a new sample x we apply the process described in Subsection 5.3 to obtain ρ Y .Then the prediction is given by: Note that this framework also allows to easily compute confidence intervals for the prediction.The model can be trained using the strategies discussed in Subsection 5.3.For gradient-based optimization we use a mean squared error loss function: where the second term correspond to the variance of the prediction and α controls the trade-off between error and variance.

Related Work
The ability of density matrices to represent probability distributions has been used in previous works.The early work by Wolf (2006) uses the density matrix formalism to perform spectral clustering, and shows that this formalism not only is able to predict cluster labels for the objects being classified, but also provides the probability that the object belongs to each of the clusters.Similarly, Tiwari and Melucci (2019) proposed a quantum-inspired binary classifier using density matrices, where samples are encoded into pure quantum states.In a similar fashion, Sergioli et al (2018) proposed a quantum nearest mean classifier based on the trace distance between the quantum state of a sample, and a quantum centroid that is a mixed state of the pure quantum states of all samples belonging to a single class.Another class of proposals directly combine these quantum ideas with customary machine learning techniques, such as frameworks for multi-modal learning for sentiment analysis (Li et al, 2021;Li et al, 2020;Zhang et al, 2018).Since its inception, random features have been used to improve the performance of several kernel methods: kernel ridge regression (Avron et al, 2017), support vector machines (SVM) (Sun et al, 2018), and nonlinear component analysis (Xie et al, 2015).Besides, random features have been used in conjunction with deep learning architectures in different works (Arora et al, 2019;Ji and Telgarsky, 2019;Li et al, 2019).
The combination of RFF and density matrices was initially proposed by González et al (2020).In that work, RFF are used as a quantum feature map, among others, and the QMC method (Subsection 5.3) was presented.In González et al (2020) the coherent state kernel showed better performance than the Gaussian kernel.It is important to notice that the coherent state kernel was calculated exactly while the Gaussian kernel was approximated using RFF.It is possible to apply RFF to approximate the coherent state kernel and use it as the quantum feature map in the models presented in this paper.The emphasis of González et al (2020) is to show that quantum measurement can be used to do supervised learning.In contrast, the present paper addresses a wider problem with several new contributions: a new method for density estimation based on density matrices and RFF, the proof of the connection between this method and kernel density estimation, and new differentiable models for density estimation, classification and regression.
The present work can be seen as a type of quantum machine learning (QML), which is generally referred as the field in the intersection of quantum computing and machine learning (Schuld et al, 2015;Schuld, 2018).In particular, the methods in this paper are in the subcategory of QML called quantum inspired classical machine learning, where theory and methods from quantum physics are borrowed and adapted to machine learning methods intended to run in classical computers.Works in this category include: quantum-inspired recommendation systems (Tang, 2019a), quantum-inspired kernel-based classification methods (Tiwari et al, 2020;González et al, 2020), conversational sentiment analysis based on density matrix-like convolutional neural networks (Zhang et al, 2019), dequantised principal component analysis (Tang, 2019b), among others.
Being a memory-based strategy, KDE suffers from large-scale, high dimensional data.Due to this issue, fast approximate evaluation of non-parametric density estimation is an active research topic.Different approaches are proposed in the literature: higher-order divide-and-conquer method (Gray and Moore, 2003), separation of near and far-field (pruning) (March et al, 2015), and hashing based estimators (HBE) (Charikar and Siminelakis, 2017).Even though the purpose of the present work was not to design methods for fast approximation of KDE, the use of RFF to speed KDE seems to be a promising research direction.Comparing DMKDE against fast KDE approximation methods is part of our future work.

Experimental Evaluation
In this section we perform some experiments to evaluate the performance of the proposed methods in different benchmark tasks.The experiments are organized in three subsections: density estimation evaluation, classification evaluation and ordinal regression evaluation.The source code of the methods and the scripts of the experiments are available at https://drive.google.com/drive/folders/16pHMLjIvr6v1zY6cMvo11EqMAMqjn3Xa as Jupyter notebooks.

Density estimation evaluation
The goal of these experiments is to evaluate the efficacy and efficiency of DMKDE to approximate a pdf.We compare it against conventional Gaussian KDE.

Data sets and experimental setup
We used three datasets: • 1-D synthetic.The first synthetic dataset corresponds to a mixture of univariate Gaussians as shown in Figure 4.The mixture weights are 0.3 and 0.7 respectively and the parameters are (µ 1 = 0, σ = 1) and (µ 1 = 5, σ = 1).We generated 10,000 samples for training and use as test dataset 1,000 samples equally spaced in the interval [−5, 10].• 2-D synthetic.This dataset corresponds to three spirals as depicted in Figure 6.The training and test datasets have 10,0000 and 1,000 points respectively, all of them generated with the same stochastic procedure.• MNIST dataset.We used PCA to reduce the original 784 dimension to 40.
The resulting vectors were scaled to [0, 1].We used stratified sampling to choose 10,000 and 1,000 samples for training and testing respectively.
We performed two types of experiments over the three datasets.In the first, we wanted to evaluate the accuracy of DMKDE.In the second, we evaluated the time to predict the density on the test set.
In the first experiment, DMKDE was run with different number of RFF to see how the dimension of the RFF representation affected the accuracy of the estimation.For the 1-D dataset, both the DMKDE prediction and the KDE prediction were compared against the true pdf using root mean squared error (RMSE).For the 2-D dataset the RMSE between the DMKDE prediction and the KDE prediction was evaluated.In the case of MNIST, and because of the small values for the density, we calculated the RMSE between the log density predicted by DMKDE and KDE.The number of eigencomponents (r) was chosen by sorting the eigenvalues in descending order and plotting them to look for the curve elbow.For the 1-D and 2-D datasets, the γ value was chosen to get a good approximation of the data density, this was visually verified.For the MNIST dataset, the γ value was chosen by looking at a histogram of pairwise distances of the data.The value of the parameters were: (γ = 16, r = 30) for the 1-D dataset, (γ = 256, r = 100) for the 2-D dataset, (γ = 1, r = 150) for the MNIST dataset.
For the second experiment, we measured the time taken to predict 1,000 test samples for both KDE and DMKDE using different number of train samples.KDE was implemented in Python using liner algebra operations accelerated by numpy.At least for the experiments reported, our implementation was faster than other KDE implementations available such as the one provided by scikit learn (https://scikit-learn.org/stable/modules/density.html), which is probably optimized for other use cases.DMKDE was implemented in Python using Tensorflow.The main reason for using Tensorflow was its ability to automatically calculate the gradient of computational graphs.KDE could not benefit from this feature, on the contrary, its performance could be hurt by Tensorflow's larger memory footprint.Another advantage of Tensorflow is its ability to generate code optimized for a GPU, so both methods were run on a 2.20 GHz dual-core Intel(R) Xeon(R) CPU without a GPU to avoid any unfair advantage.

Results and discussion
Figure 5 shows how the accuracy of DMKDE increases with an increasing number of RFF.For each configuration 30 experiments were run and the blue solid line represents the mean RMSE of the experiments and the blue region represents the 95% confidence interval.In all the three datasets, 2 10 RFF achieved a low RMSE.The variance also decreases when the number of RFF is increased.For the 1-D dataset both KDE and DMKDE are compared against the true density.For the two other datasets the difference between KDE and DMKDE is calculated.In all the cases the RMSE is calculated.The blue shaded zone represents the 95% confidence interval.
Figure 6 shows the 2-D spirals dataset (left) and the density estimation of both KDE (center) and DMKDE (right).The density calculated by DMKDE is very close to the one calculated with KDE.
Figure 7 shows a comparison of the log density predicted by KDE and DMKDE.Both models were applied to test samples and samples generated randomly from a uniform distribution.As expected points are clustered around the diagonal.The DMKDE log density of test samples (left) seems to be more accurately predicted than the one of random samples.The reason is that the density of random samples is smaller than the density of test samples and the difference is amplified by the logarithm.

Classification evaluation
In this set of experiments, we evaluated DMKDC over different well known benchmark classification datasets.

Data sets and experimental setup
Six benchmark data sets were used.The details of these datasets are shown in Table 1.In the case of Gisette and Cifar, we applied a principal component analysis algorithm using 400 principal components in order to reduce the dimension.DMKDC was trained using the estimation strategy (DMKDC) and an ADAM stochastic gradient descent strategy (DMKDC-SGD).As baseline we compared against a linear support vector machine (SVM) trained using the  same RFF as DMKDC.The SVM was trained using the LinearSVC model from scikit-learn, which is based in an efficient C implementation tailored to linear SVMs.In the case of MNIST and Cifar, we additionally built a union of a LeNet architecture (LeCun et al, 1989), as a feature extraction block, and DMKDC-SGD as the classifier layer.The LeNet block is composed of two To make the comparison with baseline models fair, in all the cases the RFF layer of DMKDC-SGD is frozen, so its weights are not modified by the stochastic gradient descent learning process.
For each data set, we made a hyper parameter search using a five-fold crossvalidation with 25 randomly generated configurations.The number of RFF was set to 1000 for all the methods.For each dataset we calculated the inter-sample median distance µ and defined an interval around γ = 1 2σ 2 .The C parameter of the SVM was explored in an exponential scale from 2 −5 to 2 10 .For the ADAM optimizer in DMKDC-SGD with and without LeNet we choose the

Results and discussion
Table 3 shows the results of the classification experiments.DMKDC is a shallow method that uses RFF, so a SVM using the same RFF is fair and strong baseline.In all the cases, except one, DMKDC-SGD outperforms the SVM, which shows that it is a very competitive shallow classification method.DMKDC trained using estimation shows less competitive results, but they are still remarkable taking into account that this is an optimization-less training strategy that only passes once over the training dataset.For MNIST and Cifar the use of a deep learning feature extractor is mandatory to obtain competitive results.The results show that DMKDC-SGD can be integrated with deep neural network architectures to obtain competitive results.The improvement on classification performance of DMKC-SGD comes at the cost of increased training time.The training of DMKDC is very efficient since it corresponds to do an average of the training density matrices.Linear SVM training is also very efficient.In contrast, DMKDC-SGD requires an iterative training process that has to be tuned to get it to converge to a good local optimum, as is the case for current deep learning models.

Ordinal regression evaluation
Many multi-class classification problems can be seen as ordinal regression problems.That is, problems where labels not only indicate class membership, but also an order.Ordinal regression problems are halfway between a classification problem and a regression problem, and given the discrete probability distribution representation used in QMR, ordinal regression seems to be a suitable problem to test it.

Data sets and experimental setup
Nine standard benchmark data sets for ordinal regression were used.The details of each data set are reported in Table 2.These data sets are originally used in metric regression tasks.To convert the task into an ordinal regression one, the target values were discretized by taking five intervals of equal length over the target range.For each set, 20 different train and test partitions are made.These partitions are the same used by Chu and Ghahramani (2005) and several posterior works, and are publicly available at http://www.gatsby.ucl.ac.uk/ ∼ chuwei/ordinalregression.html.The models were evaluated using the mean absolute error (MAE), which is a popular and widely used measure in ordinal regression (Gutiérrez et al, 2016;Garg and Manwani, 2020).QMR was trained using the estimation strategy (QMR) and an ADAM stochastic gradient descent strategy (QMR-SGD).For each data set, and for each one of the 20 partitions, we made a hyper parameter search using a fivefold cross-validation procedure.The search was done generating 25 different random configuration.The range for γ was computed in the same way as for the classification experiments, β ∈ (0, 25), the number of RFF randomly chosen between the number of attributes and 1024, and the number of eigencomponents of the factorization was chosen from {0.1, 0.2, 0.5, 1} where each number represents a percentage of the RFF.For the ADAM optimizer in QMR-SGD we choose the learning rate in the interval (0, 0.001] and α ∈ (0, 1).The RFF layer was always set to trainable, and the criteria for selecting the best parameter configuration was the MAE performance.

Results and discussion
For each data set, the means and standard deviations of the test MAE for the 20 partitions are reported in Table 4, together with the results of previous state-of-the-art works on ordinal regression: Gaussian Processes (GP) and support vector machines (SVM) (Chu and Ghahramani, 2005), Neural Network Rank (NNRank) (Cheng et al, 2008), Ordinal Extreme Learning Machines (ORELM) (Deng et al, 2010) and Ordinal Regression Neural Network (ORNN) (Fernandez-Navarro et al, 2014).
QMR-SGD shows a very competitive performance.It outperforms the baseline methods in six out of the nine data sets.The training strategy based on estimation, QMR, did not performed as well.This evidences that for this problem a fine tuning of the representation is required and it is successfully accomplished by the gradient descent optimization.

Conclusions
The mathematical framework underlying quantum mechanics is a powerful formalism that harmoniously combine linear algebra and probability in the form of density matrices.This paper has shown how to use these density matrices as a building block for designing different machine learning models.The main contribution of this work is to show a novel perspective to learning that combines two very different and seemingly unrelated tools, random features and density matrices.The, somehow surprising, connection of this combination with kernel density estimation provides a new way of representing and learning probability density functions from data.The experimental results showed some evidence that this building block can be used to build competitive models for some particular tasks.However, the full potential of this new perspective is still to be explored.Examples of directions of future inquire include using complex valued density matrices, exploring the role of entanglement and exploiting the battery of practical and theoretical tools provided by quantum information theory.
Fig. 4 1-D synthetic dataset.The gray zone is the area of the true density.The estimated pdf by DMKDE (γ = 2) and KDE (γ = 4) is shown.

Fig. 5
Fig.5Accuracy of the density estimation of DMKDE for different number of RFF for the 1-D dataset (top left), 2-D dataset (top right) and MNIST dataset (bottom).For the 1-D dataset both KDE and DMKDE are compared against the true density.For the two other datasets the difference between KDE and DMKDE is calculated.In all the cases the RMSE is calculated.The blue shaded zone represents the 95% confidence interval.

Fig
Fig. 6 2-D spirals dataset (top left) and the density estimation of both KDE (top right) and DMKDE (bottom).

Figure 8
Figure 8 shows the time of both methods for different sizes of the training dataset.The prediction time of KDE depends on the size of the training dataset, while the time of DMKDE does not depend on it.The advantage of DMKDE in terms of computation time is clear for training datasets above 10 4 data samples.

Fig. 7
Fig. 7 Scatter-plots comparing the log density predicted by KDE and DMKDE: test samples (top left), uniformly random generated samples (top right), both test and random samples (bottom).

Table 1
Data sets used for classification evaluation.

Table 2
Specifications of the data sets used for ordinal regression evaluation.Train and Test indicate the number of samples, which is the same for all the twenty partitions.

Table 3
Accuracy test results for DMKDC and DMKDC-SGD compared against a linear support vector machine over RFF (SVM-RFF).Two deep learning models are also evaluated on the two image datasets: a convolutional neural network (LeNet) and its combination with DMKDC-SGD (LeNet DMKDC).

Table 4
MAE test results for QMR, QMR-SGD and different baseline methods: support vector machines (SVM), Gaussian Processes (GP), Neural Network Rank (NNRank), Ordinal Extreme Learning Machines (ORELM) and Ordinal Regression Neural Network (ORNN).The results are the mean and standard deviation of the MAE for the twenty partitions of each dataset.The best result for each data set is in bold face.