∎

11institutetext: Mahdi Zarei 22institutetext: University of California San
Francisco,
22email<EMAIL_ADDRESS>
# Feature selection algorithm based on Catastrophe model to improve the
performance of regression analysis

Mahdi Zarei

###### Abstract

In this paper we introduce a new feature selection algorithm to remove the
irrelevant or redundant features in the data sets. In this algorithm the
importance of a feature is based on its fitting to the Catastrophe model.
Akaike information criterion value is used for ranking the features in the
data set. The proposed algorithm is compared with well-known RELIEF feature
selection algorithm. Breast Cancer, Parkinson Telemonitoring data and Slice
locality data sets are used to evaluate the model.

###### Keywords:

Feature selection Catastrophe theory Akaike information criterion RELIEF
feature selection algorithm Regression analysis

##  1 Introduction

Finding the informative features from a data is a complicated process. Many
algorithms have been developed to remove the irrelevant features in the data
set and improve the performance of analysis. For example multivariate feature
selection statistics is used to reduce the complexity of the data analysis
norman2006beyond . Dimension reduction is another method to select informative
features that many researchers applied to the features in the data
ku2008comparison ; o2007theoretical ; mourao2006impact .

In this paper, we introduce a new feature selection algorithm to improve
performance of regression analysis. Akaike information criterion value is used
for ranking the features in the data set. The proposed algorithm is compared
with well-known RELIEF feature selection algorithm. This algorithm is able to
significantly reduce the number of features in this data set improving
regression analysis accuracy.

Since our algorithm is based on the approaches from Catastrophe theory and
Akaike information criterion, we start with a brief description of them.

##  2 Cusp Catastrophe

In this section we give a brief description of cusp model. Consider the
following dynamical system:

$\frac{\partial y}{\partial t}=-\frac{\partial V(y;c)}{\partial t},y\in
R^{k},c\in R^{p},$ (1)

where $V$ is the potential function, $y(t)$ represents the system’s state
variable(s), $c$ shows one or multiple (control) parameter(s) whose value(s)
determine the specific structure of the system. If $y$ is at a point where

$\frac{\partial V(y;c)}{\partial t}=0$ (2)

the system is in equilibrium. The function $V(y;c)$ acquires a minimum with
respect to $y$ at a non-equilibrium point. Equilibrium points that correspond
to minima of $V(y;c)$ are stable equilibrium points because the system will
return to such a point after a small perturbation to the system’s state. The
equilibrium points that correspond to maxima of $V(y;c)$ are unstable
equilibrium points because a perturbation of the system’s state will cause the
system to move away from the equilibrium point towards a stable equilibrium
point. Equilibrium points that correspond neither to maxima nor to minima of
$V(y;c)$, at which the Hessian matrix ($\partial^{2}V(y)/\partial
y_{i}\partial y_{j}$) has eigenvalues equal to zero, are called degenerate
equilibrium points. When the control variables of the system are changed.
System can give rise to unexpected bifurcations in its equilibrium states at
these points when the control variables of the system are changed
saunders1980introduction ; zeeman1977catastrophe ; grasman2009fitting .
Cusp model that is the simplest form of Catastrophe and can be formulated as
follows:

$-V(y;\alpha,\beta)=\alpha y+\frac{1}{2}\beta y^{2}-\frac{1}{4}y^{4},$ (3)

where $V$ is the canonical form of the potential function for the Cusp model
and its equilibrium points is a function of the control parameters $\alpha$
and $\beta$ (see Figure 1). The control parameters are the solution to the
equation

$\alpha+\beta y-y^{3}=0.$ (4)

This equation has one solution if $\delta=27\alpha-4\beta^{3}$ that is greater
than zero, and has three solutions if $\delta<0$ cobb1980estimation ;
grasman2009fitting .

Figure 1: Cusp surface cobb1980statistical

##  3 Akaike information criterion

Akaike information criterion (AIC) is a model quality measure for a given data
akaike1974new ; burnham2004multimodel . For a model AIC measure can be defined
as follow sakamoto1986akaike ; bozdogan2000akaike :

$AIC=-2log\mathrm{L}(\widehat{\theta})+2k,$ (5)

where $\mathrm{L}(\widehat{\theta})$ is the maximized likelihood function and
$k$ is the number of free parameters in the model. The smaller value of AIC
shows that data is the better fit to model. In the proposed algorithm, we used
the reverse value of AIC for ranking the features in our data.

##  4 The feature selection algorithm

In the Catastrophe theory, small change in certain parameters of a system can
cause equilibria to appear or disappear thom1983mathematical ;
zeeman1977catastrophe . We used this characteristic of the Catastrophe model
to find the features that are more affective in regression analysis. In the
proposed algorithm the features that better change the dynamic of outcome
feature or features are considered as informative features. Assume that we are
given a data set $A$ with $N$ features that $z$ is outcome feature. The
algorithm takes each feature $i$ from the data set and considers it as
bifurcation variable in the Cusp Catastrophe model. If this variable affects
the dynamic of the system (outcome feature), it is the informative feature.
The AIC value of the Cusp model is computed for each feature for ranking. The
ranking of a feature $i$ can be formulated as follows:

$AIC_{i}=AIC(-V(y;\alpha,i)),$ (6)

where $V$ is the potential function for the Cusp model (see Equation 3),
$AIC_{i}$ is the AIC value of the Cusp model for the feature $i$ as
bifurcation value ($\beta$) and $\alpha$ is the asymmetric value in the Cusp
model. Figure 2 shows the preparing the input parameters for Cusp model where
the outcome feature is considered as the state variable and the features $i$
and the last features are considered as bifurcation and asymmetric values,
respectively. The state variable and control values can be computed as follows
grasman2009fitting :

$\displaystyle y[t]=w[0]+w[1]*Y[t,1]+...+w[p]*Y[t,p],$ (7)
$\displaystyle\alpha[t]=a[0]+a[1]*X[t,1]+...+a[p]*X[t,p],$ (8)
$\displaystyle\beta[t]=b[0]+b[1]*X[t,1]+...+b[p]*X[t,p],$ (9)

where $X[t,p]$’s are independent and $Y[t,p]$’s are dependent features in the
data set. The vectors $a[j]$’s, $b[j]$’s and $w[j]$’s are estimated by means
of maximum likelihood. The rank of each feature $i$ in the data set can be
calculated as follows:

$rank_{i}\leftarrow\frac{1}{AIC_{i}}.$ (10)

Figure 2: Preparing input features for Cusp Catastrophe model

More details about the model are shown in the Algorithm 1.

Algorithm 1 Feature selection algorithm based on the Cusp Catastrophe model
and AIC ranking

1:(Initialization) $N\leftarrow$ Number of features ,$N_{F}\leftarrow$ Number
of informative features , $\alpha\leftarrow feature_{N}$, $i\leftarrow 1$ and
$\alpha$ is asymmetric variable

2: Let $\beta\leftarrow feature_{i}$ be bifurcation value in the Cusp model

3: (Fitting the Cusp model using $\alpha$ and $\beta$) Let $AIC_{i}$ be the
Akaike information criterion value of the fitting Cusp model using parameters
$\alpha$ and $\beta$

4: (Ranking the feature) $rank_{i}\leftarrow\frac{1}{AIC_{i}}$ is the rank of
feature $i$ in the dataset

5: if $\frac{1}{AIC_{i}}\leq t$ then $feature_{i}$ is not informative and
eliminate it, $i\leftarrow i+1$ and go to 6

6: (Stopping criterion) if $i>N_{F}$ stop. Otherwise go to Step 2

7: (Retraining informative features) Return $N_{F}$ informative features.

Here $N$ is the number of all feature in the data set and $N_{F}$ ($N_{F}<N$)
is the number of informative features. For all features $i$ of the data set
their rank in the data set is computed ($rank_{i}$). The set of informative
features with $N_{F}$ features is the outcome of the algorithm.

##  5 RELIEF feature selection algorithm

Next, we give a brief description of the RELIEF algorithm. More detailed
description can be found in Kira1992 ; Kononenko1994 ; Robnik-Sikonja1997 .
For a given data set with $m$ samples, and threshold of relevancy $\tau$
($0\leq\tau\leq 1$), it detects those features which are statistically
relevant to the target concept ($Y=f(X)$). Differences of feature value
between two instances $X$ and $Y$ are defined by the following function $diff$
kira1992feature .

$diff(x_{k},y_{k})=(x_{k}-y_{k})/nu_{k},$ (11)

where $nu_{k}$ is a normalization unit to normalize the values of $diff$ into
the interval $[0,l]$. RELIEF picks a sample composed of $m$ triplets of an
instance $X$, it’s same-class instance ($nearHit$) and closest different-class
instance ($nearMiss$). RELIEF uses the $p$-dimensional Euclidean distance for
selecting $nearHit$ and $nearMiss$. In every routine the feature weight $W$
vector is updated as follows:

$W_{i}=W_{i-1}-(x_{i}-nearHit_{i})^{2}+(x_{i}-nearMiss_{i})^{2}.$ (12)

Then the average feature weight vector relevance is determined for every
sample triple. Finally, it chooses the features whose average weight is above
the given threshold $\tau$.

##  6 Experimental results

The effectiveness of the proposed algorithm is verified using three different
data sets: Parkinson’s Telemonitoring, Breast Cancer and Slice locality from
UCI machine learning repository blake1998uci . Numerical experiments have been
carried out on a PC with Processor Intel(R) Core(TM) i5-3470S CPU 2.90 GHz and
8 GB RAM running under Windows 7.
In numerical experiments we apply the proposed algorithm to find a ranking
sequence of features in data sets. Then we apply different regression analysis
algorithms from WEKA to compute regression error with subsets of features. The
following regression analysis algorithms from WEKA are used in numerical
experiments:

  * •

Linear regression: Linear regression finds the best curve to fit the data by
computing the relationship between a scalar dependent variable $y$ and one or
more explanatory variables denoted $X$. It applies least squares, which
minimizes the sum of the distance from the line for each of points. The actual
observations, $y_{i}$, may be slightly off the population line because of
variability in the population. The equation is
$y_{i}=\beta_{0}+\beta_{1}x_{i}+\epsilon_{i}$, where $\epsilon_{i}$ is the
deviation from the population line which is called the residual
barlow1993numerical ; neter1983applied .

  * •

K nearest neighbors regressor: The algorithm computes the mean of the function
values of its $K$-nearest neighbours kramer2011unsupervised .

  * •

M5Rulles: It generates rules for numeric prediction by separate-and-conquer
and at each iteration builds a model tree using M5 and makes the ”best” leaf
into a rule Holmes1999 ; Quinlan1992 ; Wang1997

  * •

REPTree: Reptree is a fast tree learner that uses reduced error pruning
witten2005data .

###  6.1 Results for Breast cancer data set

Breast Cancer Wisconsin (Prognostic) Data Set contains 30 features with 569
samples. Each record represents follow-up data for one breast cancer case
mangasarian1995breast ; street1995inductive  . Table 1 presents the error of
analysing the data using for regression analysis algorithms. The second row
shows the number of features before and after feature selection. Results from
this table demonstrate that features selected by the proposed algorithm allow
us to reduce the mean absolute error (MAE) regression. MAE is calculated as
follows:

$MAE=\frac{1}{n}\sum_{i=1}^{n}\left|f_{i}-y_{i}\right|,$ (13)

where $n$ is the number of observation, $f_{i}$ is the predicted and $y_{i}$
is the true values. Although this data set is not noisy the proposed algorithm
is able to significantly reduce the number of features without deteriorating
the regression error. Regression errors with the subsets of features which are
better than that of for all features are presented in bold font.

Table 1: Performance of regression analysis algorithms for breast cancer data set | Original data | After feature selection
---|---|---
Number of features | 30 | 25 | 20 | 15 | 10 | 6 | 5
Linear Regression | 0.003 | 0.003 | 0.003 | 0.003 | 0.003 | 0.003 | 0.004
IBK | 0.008 | 0.008 | 0.008 | 0.007 | 0.007 | 0.006 | 0.007
M5P | 0.003 | 0.003 | 0.003 | 0.003 | 0.003 | 0.003 | 0.004
M5Rules | 0.003 | 0.003 | 0.003 | 0.003 | 0.003 | 0.003 | 0.004

###  6.2 Results for Slice locality data set

Slice locality data set consists of 384 features extracted from 53500 CT
images. The CT images are from 74 different patients (43 male, 31 female). The
class variable of this data set is the location of the CT slice on the axial
axis of the human body graf20112d . This data set is available on UCI Machine
Learning Repository.

Results for 10 subjects of Slice locality data set are presented in Tables
2-5. In these tables regression error obtained by regression algorithms are
given. The second line in all tables contains a number of features of original
data and after feature selection. Table 2 presents results for all subjects
using IBK algorithm. One can see that the IBK algorithm achieved the better
accuracy for all subjects data set except subject number 10 using 380
features. Table 3 presents results for all subjects using Logistic regression
algorithm. The use of the proposed algorithm allows improving the performance
of Logistic regression using 250 features for Subject 1 and 150 features for
Subjects 2 and 3. The best performance for Subject 5 achieved using 100
features. Results are almost the same for other Subjects.

Tables 4 and 5 show results for all patients using M5P and M5Rules algorithms,
respectively. Results for these two algorithms are very similar and one can
see that the proposed algorithm can improve the accuracy of regression
algorithms.

Table 2: IBK algorithm performance for 10 subjects from Slice locality data | Original data | After feature selection
---|---|---
Number of features | 385 | 380 | 350 | 300 | 250 | 200 | 150 | 100
Patient1 | 0.059 | 0.059 | 0.059 | 0.060 | 0.061 | 0.065 | 0.063 | 0.083
Patient2 | 0.080 | 0.080 | 0.081 | 0.081 | 0.082 | 0.083 | 0.085 | 0.103
Patient3 | 0.076 | 0.076 | 0.076 | 0.075 | 0.076 | 0.077 | 0.086 | 0.115
Patient4 | 0.060 | 0.060 | 0.060 | 0.061 | 0.062 | 0.063 | 0.066 | 0.081
Patient5 | 0.078 | 0.078 | 0.078 | 0.079 | 0.080 | 0.088 | 0.086 | 0.090
Patient6 | 0.349 | 0.349 | 0.349 | 0.349 | 0.336 | 0.346 | 0.456 | 0.466
Patient7 | 0.081 | 0.081 | 0.081 | 0.081 | 0.081 | 0.087 | 0.091 | 0.099
Patient8 | 0.087 | 0.087 | 0.087 | 0.087 | 0.086 | 0.086 | 0.093 | 0.099
Patient9 | 0.364 | 0.364 | 0.370 | 0.370 | 0.364 | 0.380 | 0.494 | 0.516
Patient10 | 0.098 | 0.098 | 0.100 | 0.104 | 0.103 | 0.105 | 0.110 | 0.139
Table 3: Logistic regression algorithm performance for 10 subjects from Slice locality data | Original data | After feature selection
---|---|---
Number of features | 385 | 380 | 350 | 300 | 250 | 200 | 150 | 100
Patient1 | 0.354 | 0.392 | 0.250 | 0.267 | 0.284 | 0.326 | 0.411 | 0.570
Patient2 | 0.496 | 0.435 | 0.398 | 0.367 | 0.332 | 0.309 | 0.376 | 0.621
Patient3 | 0.258 | 0.256 | 0.266 | 0.228 | 0.226 | 0.226 | 0.247 | 0.361
Patient4 | 0.282 | 0.294 | 0.305 | 0.281 | 0.294 | 0.269 | 0.373 | 0.476
Patient5 | 0.928 | 1.742 | 2.413 | 0.512 | 0.440 | 0.469 | 0.572 | 0.529
Patient6 | 0.435 | 0.439 | 0.456 | 0.440 | 0.456 | 0.572 | 2.232 | 1.514
Patient7 | 0.515 | 0.500 | 0.460 | 0.426 | 0.420 | 0.414 | 0.443 | 0.756
Patient8 | 1.306 | 1.272 | 1.275 | 1.275 | 1.449 | 1.234 | 1.457 | 2.025
Patient9 | 0.549 | 0.539 | 0.567 | 0.532 | 0.497 | 0.860 | 1.857 | 7.839
Patient10 | 0.570 | 0.565 | 0.513 | 0.522 | 0.508 | 0.492 | 0.506 | 0.681
Table 4: M5P algorithm performance for 10 subjects from Slice localization data | Original data | After feature selection
---|---|---
Number of features | 385 | 380 | 350 | 300 | 250 | 200 | 150 | 100
Patient1 | 0.299 | 0.299 | 0.301 | 0.297 | 0.294 | 0.293 | 0.298 | 0.338
Patient2 | 0.455 | 0.455 | 0.440 | 0.443 | 0.441 | 0.471 | 0.451 | 0.452
Patient3 | 0.352 | 0.352 | 0.352 | 0.349 | 0.358 | 0.343 | 0.342 | 0.337
Patient4 | 0.341 | 0.347 | 0.348 | 0.350 | 0.339 | 0.310 | 0.319 | 0.325
Patient5 | 0.458 | 0.458 | 0.427 | 0.404 | 0.395 | 0.375 | 0.385 | 0.396
Patient6 | 1.334 | 1.297 | 1.289 | 1.326 | 1.357 | 1.136 | 1.229 | 1.291
Patient7 | 0.472 | 0.467 | 0.472 | 0.472 | 0.469 | 0.476 | 0.475 | 0.490
Patient8 | 0.782 | 0.797 | 0.801 | 0.801 | 0.720 | 0.744 | 0.728 | 0.728
Patient9 | 1.214 | 1.214 | 1.175 | 1.189 | 1.152 | 1.020 | 1.683 | 1.754
Patient10 | 0.561 | 0.546 | 0.542 | 0.513 | 0.513 | 0.519 | 0.509 | 0.519
Table 5: M5Rules algorithm performance for 10 subjects for Slice localization data | Original data | After feature selection
---|---|---
Number of features | 385 | 380 | 350 | 300 | 250 | 200 | 150 | 100
Patient1 | 0.331 | 0.319 | 0.313 | 0.368 | 0.370 | 0.322 | 0.272 | 2.217
Patient2 | 0.455 | 0.455 | 0.360 | 0.339 | 0.347 | 0.557 | 0.445 | 0.490
Patient3 | 0.508 | 0.508 | 0.508 | 0.477 | 0.432 | 0.413 | 0.388 | 0.420
Patient4 | 0.328 | 0.307 | 0.311 | 0.328 | 0.333 | 0.294 | 0.309 | 0.317
Patient5 | 0.481 | 0.479 | 0.410 | 0.507 | 0.508 | 0.458 | 0.492 | 0.412
Patient6 | 1.562 | 1.320 | 1.231 | 1.313 | 1.338 | 1.030 | 1.480 | 1.242
Patient7 | 0.783 | 0.783 | 0.784 | 0.783 | 0.559 | 0.500 | 0.412 | 0.611
Patient8 | 0.686 | 0.687 | 0.696 | 0.696 | 0.853 | 0.822 | 0.755 | 2.506
Patient9 | 1.476 | 1.476 | 1.220 | 1.249 | 1.162 | 1.260 | 0.968 | 1.952
Patient10 | 0.815 | 0.693 | 0.727 | 0.714 | 0.688 | - | 1.926 | 0.586

###  6.3 Results for Parkinsons Telemonitoring data set

In this paper, we present the results for Parkinsons Telemonitoring data set.
This data set composed of a range of biomedical voice measurements from 42
people with early-stage Parkinson’s disease. Here we analyzed 15 subjects from
this data set. Results for subjects of Parkinsons Telemonitoring data set are
presented in Tables 6-9. This is illustration of a number of features in
original data and after feature selection. The number of features in original
data is 18.
Table 6 shows the results for the error of the data using IBK regressor
algorithm. The use of a very small subset of features can provide better
performance for almost all subjects. Table 7 presents the results for Logistic
regression algorithm. The proposed algorithm can reduce the error of more than
70% of cases. The situation is almost the same for the M5P algorithm 8, but
M5Rulles algorithm provides better performance and the accuracy is increased
for all subjects except Subjects 14 and 15.

Table 6: IBK algorithm performance for Parkinson’s disease data | Original data | After feature selection
---|---|---
Number of features | 18 | 11 | 10 | 9 | 8 | 7 | 6 | 5 | 4
Subject1 | 0.037 | 0.038 | 0.037 | 0.038 | 0.038 | 0.040 | 0.044 | 0.041 | 0.042
Subject2 | 0.039 | 0.037 | 0.039 | 0.038 | 0.039 | 0.040 | 0.036 | 0.040 | 0.042
Subject3 | 0.030 | 0.027 | 0.027 | 0.027 | 0.027 | 0.026 | 0.029 | 0.027 | 0.027
Subject4 | 0.039 | 0.034 | 0.034 | 0.035 | 0.035 | 0.036 | 0.034 | 0.035 | 0.037
Subject5 | 0.037 | 0.033 | 0.032 | 0.032 | 0.030 | 0.029 | 0.029 | 0.030 | 0.031
Subject6 | 0.034 | 0.037 | 0.035 | 0.033 | 0.033 | 0.033 | 0.031 | 0.031 | 0.031
Subject7 | 0.040 | 0.033 | 0.033 | 0.034 | 0.030 | 0.034 | 0.035 | 0.036 | 0.035
Subject8 | 0.032 | 0.031 | 0.033 | 0.034 | 0.032 | 0.033 | 0.036 | 0.036 | 0.036
Subject9 | 0.041 | 0.038 | 0.038 | 0.039 | 0.039 | 0.037 | 0.036 | 0.036 | 0.039
Subject10 | 0.044 | 0.037 | 0.039 | 0.039 | 0.044 | 0.046 | 0.044 | 0.042 | 0.040
Subject11 | 0.022 | 0.022 | 0.022 | 0.021 | 0.021 | 0.021 | 0.020 | 0.021 | 0.023
Subject12 | 0.030 | 0.024 | 0.024 | 0.028 | 0.029 | 0.030 | 0.032 | 0.032 | 0.030
Subject13 | 0.040 | 0.042 | 0.044 | 0.042 | 0.047 | 0.039 | 0.051 | 0.049 | 0.049
Subject14 | 0.032 | 0.030 | 0.030 | 0.031 | 0.031 | 0.031 | 0.031 | 0.033 | 0.033
Subject15 | 0.032 | 0.031 | 0.030 | 0.030 | 0.032 | 0.031 | 0.030 | 0.032 | 0.032
Table 7: Linear regression algorithm performance for Parkinson’s disease data | Original data | After feature selection
---|---|---
Number of features | 18 | 11 | 10 | 9 | 8 | 7 | 6 | 5 | 4
Subject1 | 0.030 | 0.028 | 0.028 | 0.028 | 0.028 | 0.028 | 0.028 | 0.029 | 0.029
Subject2 | 0.028 | 0.028 | 0.030 | 0.031 | 0.030 | 0.030 | 0.030 | 0.030 | 0.030
Subject3 | 0.018 | 0.019 | 0.018 | 0.020 | 0.020 | 0.020 | 0.022 | 0.021 | 0.021
Subject4 | 0.029 | 0.028 | 0.027 | 0.027 | 0.027 | 0.027 | 0.027 | 0.027 | 0.028
Subject5 | 0.024 | 0.025 | 0.025 | 0.025 | 0.025 | 0.026 | 0.026 | 0.026 | 0.027
Subject6 | 0.024 | 0.025 | 0.024 | 0.025 | 0.025 | 0.025 | 0.025 | 0.025 | 0.025
Subject7 | 0.024 | 0.024 | 0.024 | 0.023 | 0.023 | 0.024 | 0.024 | 0.024 | 0.025
Subject8 | 0.027 | 0.031 | 0.035 | 0.034 | 0.034 | 0.034 | 0.033 | 0.031 | 0.034
Subject9 | 0.029 | 0.029 | 0.030 | 0.030 | 0.030 | 0.030 | 0.037 | 0.037 | 0.038
Subject10 | 0.033 | 0.033 | 0.033 | 0.033 | 0.032 | 0.032 | 0.032 | 0.032 | 0.034
Subject11 | 0.017 | 0.017 | 0.016 | 0.016 | 0.016 | 0.016 | 0.016 | 0.017 | 0.017
Subject12 | 0.019 | 0.018 | 0.017 | 0.021 | 0.021 | 0.021 | 0.020 | 0.020 | 0.021
Subject13 | 0.031 | 0.030 | 0.031 | 0.030 | 0.032 | 0.033 | 0.033 | 0.033 | 0.035
Subject14 | 0.024 | 0.020 | 0.019 | 0.019 | 0.019 | 0.020 | 0.020 | 0.020 | 0.027
Subject15 | 0.019 | 0.020 | 0.018 | 0.018 | 0.018 | 0.018 | 0.021 | 0.021 | 0.022
Table 8: M5P algorithm performance for Parkinson’s disease data  | Original data | After feature selection
---|---|---
Number of features | 18 | 11 | 10 | 9 | 8 | 7 | 6 | 5 | 4
Subject1 | 0.030 | 0.028 | 0.028 | 0.028 | 0.028 | 0.028 | 0.029 | 0.029 | 0.029
Subject2 | 0.027 | 0.028 | 0.029 | 0.030 | 0.030 | 0.030 | 0.030 | 0.030 | 0.030
Subject3 | 0.018 | 0.019 | 0.018 | 0.020 | 0.020 | 0.020 | 0.022 | 0.021 | 0.021
Subject4 | 0.028 | 0.027 | 0.027 | 0.027 | 0.025 | 0.027 | 0.027 | 0.027 | 0.028
Subject5 | 0.024 | 0.023 | 0.023 | 0.023 | 0.023 | 0.022 | 0.022 | 0.023 | 0.024
Subject6 | 0.025 | 0.025 | 0.024 | 0.025 | 0.025 | 0.025 | 0.025 | 0.025 | 0.025
Subject7 | 0.024 | 0.024 | 0.024 | 0.023 | 0.023 | 0.024 | 0.023 | 0.024 | 0.025
Subject8 | 0.024 | 0.026 | 0.029 | 0.029 | 0.030 | 0.030 | 0.029 | 0.030 | 0.030
Subject9 | 0.029 | 0.029 | 0.028 | 0.028 | 0.029 | 0.029 | 0.031 | 0.031 | 0.031
Subject10 | 0.034 | 0.034 | 0.034 | 0.034 | 0.032 | 0.032 | 0.032 | 0.032 | 0.032
Subject11 | 0.017 | 0.017 | 0.016 | 0.016 | 0.016 | 0.016 | 0.016 | 0.017 | 0.017
Subject12 | 0.020 | 0.018 | 0.017 | 0.021 | 0.021 | 0.021 | 0.020 | 0.020 | 0.021
Subject13 | 0.033 | 0.031 | 0.031 | 0.031 | 0.032 | 0.033 | 0.033 | 0.033 | 0.035
Subject14 | 0.019 | 0.019 | 0.019 | 0.019 | 0.019 | 0.021 | 0.020 | 0.020 | 0.023
Subject15 | 0.019 | 0.021 | 0.019 | 0.019 | 0.019 | 0.019 | 0.023 | 0.023 | 0.022
Table 9: M5Rules algorithm performance for Parkinson’s disease data | Original data | After feature selection
---|---|---
Number of features | 18 | 11 | 10 | 9 | 8 | 7 | 6 | 5 | 4
Subject1 | 0.030 | 0.028 | 0.029 | 0.029 | 0.028 | 0.028 | 0.029 | 0.029 | 0.029
Subject2 | 0.027 | 0.029 | 0.029 | 0.031 | 0.030 | 0.030 | 0.030 | 0.030 | 0.030
Subject3 | 0.019 | 0.019 | 0.018 | 0.020 | 0.020 | 0.021 | 0.022 | 0.021 | 0.021
Subject4 | 0.029 | 0.028 | 0.028 | 0.027 | 0.026 | 0.028 | 0.030 | 0.030 | 0.028
Subject5 | 0.025 | 0.023 | 0.023 | 0.023 | 0.023 | 0.022 | 0.022 | 0.023 | 0.024
Subject6 | 0.024 | 0.025 | 0.024 | 0.025 | 0.025 | 0.025 | 0.025 | 0.025 | 0.025
Subject7 | 0.024 | 0.024 | 0.024 | 0.023 | 0.023 | 0.025 | 0.024 | 0.024 | 0.025
Subject8 | 0.029 | 0.038 | 0.043 | 0.043 | 0.031 | 0.031 | 0.044 | 0.045 | 0.031
Subject9 | 0.031 | 0.030 | 0.029 | 0.029 | 0.030 | 0.031 | 0.032 | 0.032 | 0.031
Subject10 | 0.035 | 0.035 | 0.035 | 0.034 | 0.033 | 0.033 | 0.033 | 0.033 | 0.033
Subject11 | 0.017 | 0.017 | 0.016 | 0.016 | 0.016 | 0.016 | 0.016 | 0.017 | 0.017
Subject12 | 0.019 | 0.018 | 0.017 | 0.021 | 0.021 | 0.021 | 0.020 | 0.020 | 0.021
Subject13 | 0.033 | 0.031 | 0.031 | 0.031 | 0.032 | 0.033 | 0.033 | 0.033 | 0.035
Subject14 | 0.019 | 0.020 | 0.020 | 0.019 | 0.019 | 0.021 | 0.020 | 0.020 | 0.023
Subject15 | 0.020 | 0.021 | 0.020 | 0.020 | 0.020 | 0.020 | 0.023 | 0.024 | 0.022

Figure 3 demonstrates applying different classifiers for Parkinson’s disease
data set. Figure 3 indicates that cusp model is reduced the error of
classifiers for almost all subjects from Parkinson’s disease data set.

Figure 3: Classification algorithms performance for Parkinson’s disease data
using all features and after feature selection using cusp catastrophe feature
selection algorithm

Figures 4 show the Equilibrium surface (3 dimensional) and control surface (2
dimensional) of fitting the most irrelevant (left) and the most significant
features in different data sets using the Cusp Catastrophe model. The
informative features have more affect on the system and put the system closer
to the bifurcation situation.

Figure 4: Cusp plot: The most least informative features (left) and the most
least informative features (right) base on proposed algorithm for subject 1 to
subject 6

Tables 10\- 15 show the ranking of the features using the proposed and RELIEF
algorithms. The ranking values are not exactly the same, but the for almost
all cases the informative features’ levels are similar in both ranking
results. For example, for the first subject, the informative features of 3,
14, 4 and 6 are in the top of the table in both algorithms and less-
significant features 2 and 17 are at the bottom.

Table 10: Ranking of the features using the proposed and RELIEF algorithms for subject 1 from Parkinsons disease data Feature selection algorithm based on the Cusp model | RELIEF algorithms
---|---
Attribute ID | Rank | Attribute ID | Rank
3 | 0.003144 | 14 | 0.030901
14 | 0.003096 | 3 | 0.014302
4 | 0.003052 | 6 | 0.014158
6 | 0.002947 | 4 | 0.011554
15 | 0.002923 | 5 | 0.009576
5 | 0.002732 | 7 | 0.009572
7 | 0.002731 | 15 | 0.006487
9 | 0.002685 | 12 | 0.004949
12 | 0.002586 | 16 | 0.004764
16 | 0.002569 | 9 | 0.004034
8 | 0.002565 | 11 | 0.001722
11 | 0.002564 | 13 | 0.0016
10 | 0.0025 | 10 | 0.001595
13 | 0.0025 | 2 | 0.000525
17 | 0.002358 | 8 | -0.00004
2 | 0.002351 | 1 | -0.00254
1 | 0.002351 | 17 | -0.00378
Table 11: Ranking of the features using the proposed and RELIEF algorithms for subject 2 from Parkinsons disease data Feature selection algorithm based on the Cusp model | RELIEF algorithms
---|---
Attribute ID | Rank | Attribute ID | Rank
16 | 0.003336 | 14 | 0.0195
6 | 0.003299 | 6 | 0.00943
14 | 0.003241 | 16 | 0.00907
9 | 0.003212 | 12 | 0.00784
12 | 0.0032 | 9 | 0.00706
4 | 0.003173 | 3 | 0.00582
3 | 0.003167 | 4 | 0.00309
15 | 0.003154 | 15 | 0.003
8 | 0.003134 | 11 | 0.0026
13 | 0.003069 | 13 | 0.00244
10 | 0.003069 | 10 | 0.00244
11 | 0.003057 | 7 | 0.00243
7 | 0.00304 | 5 | 0.00242
5 | 0.00304 | 8 | 0.0023
2 | 0.002727 | 2 | 0.00133
1 | 0.002727 | 1 | 0.00107
17 | 0.002718 | 17 | -0.00237
Table 12: Ranking of the features using the proposed and RELIEF algorithms for subject 3 from Parkinsons disease data Feature selection algorithm based on the Cusp model | RELIEF algorithms
---|---
Attribute ID | Rank | Attribute ID | Rank
15 | 0.003585 | 15 | 0.024669
6 | 0.003473 | 14 | 0.018446
3 | 0.003261 | 6 | 0.016579
4 | 0.003031 | 3 | 0.013203
14 | 0.002994 | 4 | 0.010286
7 | 0.002946 | 5 | 0.007498
9 | 0.002946 | 7 | 0.00748
5 | 0.002945 | 11 | 0.005778
12 | 0.002937 | 12 | 0.003904
8 | 0.002934 | 9 | 0.00329
11 | 0.002933 | 1 | 0.003219
10 | 0.00287 | 8 | 0.002655
13 | 0.002869 | 10 | 0.002304
16 | 0.002627 | 13 | 0.002297
1 | 0.002595 | 17 | 0.002161
2 | 0.002589 | 2 | 0.000729
17 | 0.002565 | 16 | -0.00093
Table 13: Ranking of the features using the proposed and RELIEF algorithms for subject 4 from Parkinsons disease data Feature selection algorithm based on the Cusp model | RELIEF algorithms
---|---
Attribute ID | Rank | Attribute ID | Rank
3 | 0.004621 | 6 | 0.02566
4 | 0.00456 | 3 | 0.02124
6 | 0.004473 | 17 | 0.01921
5 | 0.003827 | 4 | 0.01823
7 | 0.003826 | 14 | 0.01734
14 | 0.003417 | 5 | 0.01714
15 | 0.003254 | 7 | 0.01711
9 | 0.002984 | 2 | 0.00843
8 | 0.002968 | 15 | 0.00774
13 | 0.002935 | 13 | 0.00711
10 | 0.002935 | 10 | 0.00711
12 | 0.00293 | 11 | 0.00695
11 | 0.00291 | 12 | 0.00676
17 | 0.002904 | 8 | 0.00671
16 | 0.002793 | 9 | 0.00613
1 | 0.002771 | 1 | 0.00519
2 | 0.00277 | 16 | 0.00168
Table 14: Ranking of the features using the proposed and RELIEF algorithms for subject 5 from Parkinsons disease data Feature selection algorithm based on the Cusp model | RELIEF algorithms
---|---
Attribute ID | Rank | Attribute ID | Rank
14 | 0.003896 | 14 | 0.02979
3 | 0.003671 | 6 | 0.02661
4 | 0.003533 | 4 | 0.02327
6 | 0.003529 | 3 | 0.01819
7 | 0.003189 | 7 | 0.01289
5 | 0.003185 | 5 | 0.01287
15 | 0.003059 | 9 | 0.01101
16 | 0.00253 | 15 | 0.0101
9 | 0.00248 | 12 | 0.00659
12 | 0.002401 | 11 | 0.00414
8 | 0.002372 | 10 | 0.00354
11 | 0.002363 | 13 | 0.00354
10 | 0.002343 | 8 | 0.00331
13 | 0.002343 | 16 | 0.00278
2 | 0.002339 | 2 | 0.00244
1 | 0.002324 | 17 | 0.00116
17 | 0.002314 | 1 | -0.00406
Table 15: Ranking of the features using the proposed and RELIEF algorithms for subject 6 from Parkinsons disease data Feature selection algorithm based on the Cusp model | RELIEF algorithms
---|---
Attribute ID | Rank | Attribute ID | Rank
15 | 0.003297 | 14 | 0.014851
4 | 0.003173 | 6 | 0.014336
3 | 0.003093 | 15 | 0.014142
6 | 0.003076 | 17 | 0.01388
14 | 0.003062 | 4 | 0.012454
7 | 0.002854 | 3 | 0.010648
5 | 0.002854 | 7 | 0.008541
9 | 0.002691 | 5 | 0.008525
12 | 0.002649 | 2 | 0.003976
8 | 0.002644 | 12 | 0.002502
11 | 0.002619 | 1 | 0.001973
16 | 0.002597 | 11 | 0.001794
10 | 0.002565 | 9 | -6.8E-05
13 | 0.002565 | 8 | -0.0006
1 | 0.002418 | 13 | -0.00153
2 | 0.00234 | 10 | -0.00153
17 | 0.002315 | 16 | -0.00244

Tables 16-23 show the mean absolute error and root mean square error for
Regression analysis before and after feature selection for 15 subjects. We
separated the results of different algorithms from each other. Tables 16 and
17 shows the results of Linear regression algorithm. The accuracy of analyzing
all subjects except subject 2, 9 and 14 using the proposed algorithm compared
with original data is improved. The RELIEF algorithm has improvement for
almost all subjects, but our algorithm has better performance than RELIEF
algorithm.

Tables 18-19 are the related results for K-nearest neighbors algorithm and
they show that both algorithms have better accuracy only for 60% of subjects
and the same situation happened for M5Rulles (see the tables 20-21) and
REPTree (22-23) algorithms, but for some subjects the RELIEF algorithm has
better performance.

Table 16: Mean absolute error of Linear regression algorithm after feature selection using the proposed and RELIEF algorithms for Slice locality data set | MAE of Linear Regression
---|---
Subject | Original data |  | After Feature selection
---
using feature selection algorithm
based on the Cusp model
| After feature
---
selection using RELIEF algorithms
1 | 0.0295 | 0.0291 | 0.0282
2 | 0.0276 | 0.028 | 0.028
3 | 0.0183 | 0.0183 | 0.0182
4 | 0.0292 | 0.029 | 0.0292
5 | 0.0235 | 0.0235 | 0.0235
6 | 0.0239 | 0.0239 | 0.024
7 | 0.0243 | 0.0242 | 0.0244
8 | 0.0266 | 0.0266 | 0.028
9 | 0.0286 | 0.0288 | 0.0288
10 | 0.0333 | 0.0333 | 0.0333
11 | 0.0169 | 0.0167 | 0.017
12 | 0.0193 | 0.0187 | 0.0194
13 | 0.0305 | 0.0297 | 0.0315
14 | 0.019 | 0.0193 | 0.0188
15 | 0.0266 | 0.0261 | 0.0266
Table 17: Root mean square error of Linear regression algorithm after feature selection using the proposed and RELIEF algorithms for Slice locality data set | RMSE of Linear Regression
---|---
Subject | Original data |  | After Feature selection
---
using feature selection algorithm
based on the Cusp model
| After feature
---
selection using RELIEF algorithm
1 | 0.0386 | 0.0381 | 0.0384
2 | 0.0372 | 0.0377 | 0.0377
3 | 0.0249 | 0.0249 | 0.0248
4 | 0.042 | 0.0418 | 0.042
5 | 0.0336 | 0.0336 | 0.0336
6 | 0.0325 | 0.0325 | 0.0322
7 | 0.0338 | 0.0338 | 0.0335
8 | 0.0401 | 0.0401 | 0.0424
9 | 0.0376 | 0.0377 | 0.0375
10 | 0.0472 | 0.0472 | 0.0461
11 | 0.0239 | 0.0237 | 0.024
12 | 0.025 | 0.0245 | 0.0256
13 | 0.0404 | 0.0392 | 0.0425
14 | 0.0248 | 0.0253 | 0.0246
15 | 0.0322 | 0.0317 | 0.0319
Table 18: Mean absolute error of IBK algorithm after feature selection using the proposed and RELIEF algorithms for Slice locality data set | MAE of IBK
---|---
Subject | Original data |  | After Feature selection
---
using feature selection algorithm
based on the Cusp model
| After feature
---
selection using RELIEF algorithm
1 | 0.037 | 0.038 | 0.042
2 | 0.0389 | 0.0411 | 0.0411
3 | 0.0304 | 0.0297 | 0.0311
4 | 0.0394 | 0.0387 | 0.0372
5 | 0.0369 | 0.0344 | 0.0356
6 | 0.034 | 0.0334 | 0.0355
7 | 0.0404 | 0.0385 | 0.0389
8 | 0.0321 | 0.032 | 0.032
9 | 0.0405 | 0.0399 | 0.0399
10 | 0.0439 | 0.044 | 0.0433
11 | 0.0218 | 0.0231 | 0.0224
12 | 0.0297 | 0.0295 | 0.0308
13 | 0.0402 | 0.0411 | 0.0402
14 | 0.0317 | 0.03 | 0.0307
15 | 0.0338 | 0.0352 | 0.0335
Table 19: Root mean square error of IBK algorithm after feature selection using the proposed and RELIEF algorithms for Slice locality data set | RMSE of IBK
---|---
Subject | Original data |  | After Feature selection
---
using feature selection algorithm
based on the Cusp model
| After feature
---
selection using RELIEF algorithm
1 | 0.0493 | 0.0506 | 0.0548
2 | 0.0526 | 0.0537 | 0.0537
3 | 0.0379 | 0.0379 | 0.0401
4 | 0.0569 | 0.0565 | 0.0567
5 | 0.0499 | 0.047 | 0.0477
6 | 0.0453 | 0.0447 | 0.0457
7 | 0.0527 | 0.0504 | 0.0507
8 | 0.0462 | 0.0458 | 0.0466
9 | 0.0531 | 0.0528 | 0.0536
10 | 0.056 | 0.0562 | 0.0538
11 | 0.0285 | 0.0316 | 0.029
12 | 0.0385 | 0.0399 | 0.0402
13 | 0.0533 | 0.0539 | 0.0532
14 | 0.038 | 0.0359 | 0.0378
15 | 0.0426 | 0.0435 | 0.0427
Table 20: Mean absolute error of M5Rules algorithm after feature selection using the proposed and RELIEF algorithms for Slice locality data set | MAE of M5Rules
---|---
Subject | Original data |  | After Feature selection
---
using feature selection algorithm
based on the Cusp model
| After feature
---
selection using RELIEF algorithm
1 | 0.0299 | 0.0299 | 0.0292
2 | 0.0273 | 0.0285 | 0.0276
3 | 0.0188 | 0.0181 | 0.0203
4 | 0.0291 | 0.0291 | 0.0278
5 | 0.0246 | 0.0248 | 0.0233
6 | 0.0241 | 0.024 | 0.024
7 | 0.0237 | 0.0233 | 0.0235
8 | 0.0286 | 0.0275 | 0.0262
9 | 0.0306 | 0.0319 | 0.0306
10 | 0.0349 | 0.0349 | 0.034
11 | 0.0167 | 0.0169 | 0.0175
12 | 0.019 | 0.0188 | 0.0197
13 | 0.0333 | 0.0313 | 0.032
14 | 0.0196 | 0.021 | 0.0209
15 | 0.0246 | 0.0249 | 0.0256
Table 21: Root mean square error of M5Rules algorithm after feature selection using the proposed and RELIEF algorithms for Slice locality data set | RMSE of M5Rules
---|---
Subject | Original data |  | After Feature selection
---
using feature selection algorithm
based on the Cusp model
| After feature
---
selection using RELIEF algorithm
1 | 0.0393 | 0.0393 | 0.0393
2 | 0.0366 | 0.0377 | 0.0366
3 | 0.0259 | 0.0252 | 0.0284
4 | 0.0423 | 0.0423 | 0.0415
5 | 0.0343 | 0.0345 | 0.0318
6 | 0.0327 | 0.0327 | 0.0327
7 | 0.033 | 0.0324 | 0.0328
8 | 0.0457 | 0.0444 | 0.0388
9 | 0.0403 | 0.0424 | 0.0401
10 | 0.0488 | 0.0488 | 0.0478
11 | 0.0219 | 0.0225 | 0.0235
12 | 0.0244 | 0.0244 | 0.026
13 | 0.044 | 0.041 | 0.0424
14 | 0.0258 | 0.0286 | 0.0275
15 | 0.0305 | 0.0305 | 0.0314
Table 22: Mean absolute error of REPTree algorithm after feature selection using the proposed and RELIEF algorithms for Slice locality data set | MAE of REPTree
---|---
Subject | Original data |  | After Feature selection
---
using feature selection algorithm
based on the Cusp model
| After feature
---
selection using RELIEF algorithm
1 | 0.0357 | 0.0357 | 0.0353
2 | 0.0344 | 0.0347 | 0.0347
3 | 0.0223 | 0.0228 | 0.0226
4 | 0.0312 | 0.0308 | 0.0304
5 | 0.0272 | 0.0273 | 0.0276
6 | 0.0278 | 0.028 | 0.0278
7 | 0.0273 | 0.0276 | 0.0276
8 | 0.03 | 0.0311 | 0.03
9 | 0.0387 | 0.0381 | 0.0387
10 | 0.0358 | 0.0358 | 0.0349
11 | 0.0183 | 0.018 | 0.0184
12 | 0.0261 | 0.0267 | 0.0261
13 | 0.043 | 0.043 | 0.0428
14 | 0.0263 | 0.0263 | 0.0263
15 | 0.0288 | 0.0289 | 0.0293
Table 23: Root mean square error of REPTree algorithm after feature selection using the proposed and RELIEF algorithms for Slice locality data set | RMSE of REPTree
---|---
Subject | Original data |  | After Feature selection
---
using feature selection algorithm
based on the Cusp model
| After feature
---
selection using RELIEF algorithm
1 | 0.0458 | 0.0458 | 0.0453
2 | 0.0449 | 0.0448 | 0.0448
3 | 0.0284 | 0.0288 | 0.0288
4 | 0.0449 | 0.0446 | 0.0437
5 | 0.0363 | 0.0363 | 0.0367
6 | 0.0371 | 0.038 | 0.0371
7 | 0.0379 | 0.0387 | 0.0383
8 | 0.0506 | 0.0547 | 0.0506
9 | 0.0519 | 0.0513 | 0.052
10 | 0.0458 | 0.0458 | 0.0454
11 | 0.0251 | 0.025 | 0.0252
12 | 0.0336 | 0.0356 | 0.0335
13 | 0.0538 | 0.0538 | 0.0538
14 | 0.0339 | 0.0338 | 0.0339
15 | 0.0362 | 0.0363 | 0.0366

Figure 5 provides a comparison between proposed algorithm and the well known
RELIEF algorithm for Slice locality data set. Mean absolute error and root
mean square error of four classifiers of original data and after feature
selection are shown in the figures. The graph show that the proposed algorithm
is improved the accuracy of classification algorithms for almost all subjects
using different classifiers.

Figure 5: Mean square error and root mean square error of classifiers after
feature selection using the proposed and RELIEF algorithms for Slice locality
data set

##  7 Conclusions

In this paper, we introduced a new feature selection algorithms to remove the
irrelevant or redundant features in the data sets. This algorithm removes the
irrelevant or redundant features of a regression data sets. This algorithm
selects significant features based on their fitting to the Catastrophe model
and the features that better change the dynamics of the outcome feature or
features are considered as informative features. The Akaike information
criterion value of the Cusp model is computed for ranking of each feature. We
applied this algorithm to three different data sets: Parkinson’s
Telemonitoring, Breast Cancer and Slice locality from UCI machine learning
repository. Results show that the proposed algorithm is efficient in finding
the significant subset of features in a data set.

## References

  * (1) Akaike, H. A new look at the statistical model identification.  Automatic Control, IEEE Transactions on 19, 6 (1974), 716–723.
  * (2) Barlow, J. L. Numerical aspects of solving linear least squares problems.  Handbook of Statistics 9 (1993), 303–376.
  * (3) Blake, C., and Merz, C. J. $\\{$UCI$\\}$ repository of machine learning databases.
  * (4) Bozdogan, H. Akaike’s information criterion and recent developments in information complexity.  Journal of mathematical psychology 44, 1 (2000), 62–91.
  * (5) Burnham, K. P., and Anderson, D. R. Multimodel inference understanding aic and bic in model selection.  Sociological methods & research 33, 2 (2004), 261–304.
  * (6) Cobb, L. Estimation theory for the cusp catastrophe model.  In Proceedings of the Section on Survey Research Methods (1980), pp. 772–776.
  * (7) Cobb, L., and Watson, B. Statistical catastrophe theory: An overview.  Mathematical Modelling 1, 4 (1980), 311–317.
  * (8) Graf, F., Kriegel, H.-P., Schubert, M., Pölsterl, S., and Cavallaro, A. 2d image registration in ct images using radial image descriptors.  In Medical Image Computing and Computer-Assisted Intervention–MICCAI 2011. Springer, 2011, pp. 607–614.
  * (9) Grasman, R. P., van der Maas, H. L., and Wagenmakers, E.-J. Fitting the cusp catastrophe in r: A cusp-package primer.  Journal of Statistical Software 32, 8 (2009), 1–28.
  * (10) Holmes, G., Hall, M., and Frank, E. Generating rule sets from model trees.  In Twelfth Australian Joint Conference on Artificial Intelligence (1999), Springer, pp. 1–12.
  * (11) Kira, K., and Rendell, L. A. The feature selection problem: Traditional methods and a new algorithm.  In AAAI (1992), pp. 129–134.
  * (12) Kira, K., and Rendell, L. A. A practical approach to feature selection.  In Ninth International Workshop on Machine Learning (1992), D. H. Sleeman and P. Edwards, Eds., Morgan Kaufmann, pp. 249–256.
  * (13) Kononenko, I. Estimating attributes: Analysis and extensions of relief.  In European Conference on Machine Learning (1994), F. Bergadano and L. D. Raedt, Eds., Springer, pp. 171–182.
  * (14) Kramer, O. Unsupervised k-nearest neighbor regression.  arXiv preprint arXiv:1107.3600 (2011).
  * (15) Ku, S.-p., Gretton, A., Macke, J., and Logothetis, N. K. Comparison of pattern recognition methods in classifying high-resolution bold signals obtained at high magnetic field in monkeys.  Magnetic resonance imaging 26, 7 (2008), 1007–1014.
  * (16) Mangasarian, O. L., Street, W. N., and Wolberg, W. H. Breast cancer diagnosis and prognosis via linear programming.  Operations Research 43, 4 (1995), 570–577.
  * (17) Mourão-Miranda, J., Reynaud, E., McGlone, F., Calvert, G., and Brammer, M. The impact of temporal compression and space selection on svm analysis of single-subject and multi-subject fmri data.  NeuroImage 33, 4 (2006), 1055–1065.
  * (18) Neter, J., Wasserman, W., and Kutner, M. H. Applied linear regression models.  Irwin Homewood, IL, 1983.
  * (19) Norman, K., Polyn, S., Detre, G., and Haxby, J. Beyond mind-reading: multi-voxel pattern analysis of fmri data.  Trends in cognitive sciences 10, 9 (2006), 424–430.
  * (20) O’Toole, A. J., Jiang, F., Abdi, H., Pénard, N., Dunlop, J. P., and Parent, M. A. Theoretical, statistical, and practical perspectives on pattern-based classification approaches to the analysis of functional neuroimaging data.  Journal of cognitive neuroscience 19, 11 (2007), 1735–1752.
  * (21) Quinlan, R. J. Learning with continuous classes.  In 5th Australian Joint Conference on Artificial Intelligence (Singapore, 1992), World Scientific, pp. 343–348.
  * (22) Robnik-Sikonja, M., and Kononenko, I. An adaptation of relief for attribute estimation in regression.  In Fourteenth International Conference on Machine Learning (1997), D. H. Fisher, Ed., Morgan Kaufmann, pp. 296–304.
  * (23) Sakamoto, Y., Ishiguro, M., and Kitagawa, G. Akaike information criterion statistics.  Dordrecht, The Netherlands: D. Reidel (1986).
  * (24) Saunders, P. T. An introduction to catastrophe theory.  Cambridge University Press, 1980.
  * (25) Street, W. N., Mangasarian, O. L., and Wolberg, W. H. An inductive learning approach to prognostic prediction.  In ICML (1995), Citeseer, pp. 522–530.
  * (26) Thom, R. Mathematical models of morphogenesis.  New York (1983).
  * (27) Wang, Y., and Witten, I. H. Induction of model trees for predicting continuous classes.  In Poster papers of the 9th European Conference on Machine Learning (1997), Springer.
  * (28) Witten, I. H., and Frank, E. Data Mining: Practical machine learning tools and techniques.  Morgan Kaufmann, 2005.
  * (29) Zeeman, E. C. Catastrophe theory: Selected papers, 1972–1977. Addison-Wesley, 1977.


Tensor Networks for Simulating Quantum Circuits on FPGAs
Maksim Levental

Most research in quantum computing today is performed against simulations
of quantum computers rather than true quantum computers. Simulating
a quantum computer entails implementing all of the unitary operators
corresponding to the quantum gates as tensors. For high numbers of
qubits, performing tensor multiplications for these simulations becomes
quite expensive, since $N$-qubit gates correspond to $2^{N}$-dimensional
tensors. One way to accelerate such a simulation is to use field programmable
gate array (FPGA) hardware to efficiently compute the matrix multiplications.
Though FPGAs can efficiently perform tensor multiplications, they
are memory bound, having relatively small block random access memory.
One way to potentially reduce the memory footprint of a quantum computing
system is to represent it as a tensor network; tensor networks are
a formalism for representing compositions of tensors wherein economical
tensor contractions are readily identified. Thus we explore tensor
networks as a means to reducing the memory footprint of quantum computing
systems and broadly accelerating simulations of such systems.

§ INTRODUCTION

Quantum computing (QC) refers to the manipulation and exploitation
of properties of quantum mechanical (QM) systems to perform computation.
QM systems exhibit properties such as superposition and entanglement
and clever quantum algorithms operate on these systems to
perform general computation. Unsurprisingly, the technique was intially
conceived of as a way to simulate physical systems themselves:

“… [N]ature isn't classical, dammit, and if you want to
make a simulation of nature, you'd better make it quantum mechanical,
and by golly it's a wonderful problem, because it doesn't look so

This closing remark from the keynote at the 1st
Physics of Computation Conference in 1981, delivered by the late Richard
Feynman [10], succinctly, but accurately,
expresses that initial goal of quantum computing. Although modeling
and simulating physical systems on quantum computers remains a thriving
area of research we narrow our focus here to QC as it pertains to
solving general computational problems. Such problems include unstructured
search [16], integer factorization [31],
combinatorial optimization [9], and many others.
It is conjectured that some quantum algorithms enable quantum computers
to exceed the computational power of classical computers [39].

QC systems are composed of so-called quantum bits, or qubits,
that encode initial and intermediate states of computations. Transformations
between states are effected by time-reversible transforms, called
unitary operators. A formalism for representing
quantum computation is the quantum circuit formalism, where
semenatically related collections of $N$ qubits are represented as
registers and transformations are represented as gates,
connected to those registers by wires, and applied in sequence.
As already mentioned, in hardware, quantum circuits correspond to
physical systems that readily exhibit quantum mechnical properties;
examples of physical qubits include transmons, ion traps and topological
quantum computers [26]. Current state of the art QC systems
are termed Noisy Intermediate-Scale Quantum (NISQ) systems. Such systems
are characterized by moderate quantities of physical qubits (50-100)
but relatively few logical qubits (i.e. qubits robust to inteference
and noise). Due to these limitations (and, more broadly, the relative
scarcity of functioning QC systems), most research on quantum algorithms
is performed with the assistance of simulators of QC systems. Such
simulators perform simulations by representing $N$-qubit circuits
as $2^{N}$-dimensional complex vectors and transformations on those
vectors as $2^{N}$-dimensional complex matrix-vector multiplication.
Naturally, due to this exponential growth, naively executing such
simulations quickly become infeasible for $N>50$ qubits [27],
both due to memory constraints and compute time.

It's the case that matrices are a subclass of a more general mathematical
object called a tensor and composition of matrices can be
expressed as tensor contraction. Tensor networks
(TNs) are decompositions (i.e. factorizations) of very high-dimensional
tensors into networks (i.e. graphs) of low-dimensional tensors. TNs
have been successfullly employed in reducing the memory requirements
of simulations of QC systems [27].
The critical feature of tensor networks that make them useful for
QC is the potential to perform tensor contractions on the low-dimensional
tensors in an order such that, ultimately, the memory and compute
time requirements are lower than for the traditional representation.
Existing applications of TNs to quantum circuits focus primarily on
memory constraints on general purpose computers [11]
and distributed environments [22].

FPGAs are known to be performant for matrix multiplication uses cases
[25]. Though FPGAs typically run at lower
clock speeds (100-300MHz) than either conventional processors or even
graphics processors they, nonetheless, excel at latency constrained
computations owing to their fully “synchronous” nature (all modules
in the same clock domain execute simultaneously). At first
glance FPGAs seem like a suitable platform for performant simulation
of quantum systems when runtime is of the essence. Unfortunately,
RAM is one of the more limited resources on an FPGA and therefore
it becomes necessary to explore memory reduction strategies for simulations
(as well as runtime reduction strategies). Hence, we explore tensor
networks as a means of reducing the memory footprint of quantum circuits
with particular attention to dimensions of the designs space as they
pertain to deployment to FPGAs.

The remainder of this report is organized as follows: section <ref>
covers the necessary background wherein subsection <ref>
very briefly reviews quantum computation and quantum circuits (with
particular focus on aspects that will be relevant for tensor networks
and FPGAs), section <ref> defines tensors and
tensor networks fairly rigorously and discusses algorithms for identifying
optimal contraction orders, section <ref> discusses the
constraints imposed by virtue of deploying to FPGA, section <ref>
describes our implementation of TNs on FPGAs, section <ref>
reports our results on various random circuits, and section <ref>
concludes with future research directions.

§ BACKGROUND

§.§ Quantum Computing

We very (very) quickly review quantum computing and quantum circuits
as they pertain to our project. For a much more pedagogically sound
introduction consult [17]. As already alluded to,
quantum computing exploits properties of quantum mechanical systems
in order to perform arbitrary computation. The fundamental unit of
quantum computation is a qubit, defined as two-dimensional quantum
system with state vector $\psi$ an element of a Hilbert space[A Hilbert space $H$ is a vector space augmented with an inner product
such that, with respect to the metric induced by that inner product,
all Cauchy sequences converge.] $H$:

\[
\psi:=\alpha\twovec 10+\beta\twovec 01\equiv\twovec{\alpha}{\beta}
\]

where $\alpha,\beta\in\mathbb{\C}$ and $\left|\alpha\right|^{2}+\left|\beta\right|^{2}=1$.
This exhibits the superposition property of the qubit[We say that the qubit is in a superposition of the basis vectors/states.]
in that the squares of the coefficients are the probabilities of measuring
the system in the corresponding basis state. Collections of qubits
have state vectors that represented by the tensor product
of the individual states of each qubit; for example, two qubits $\psi,\phi$
have state vector

\[
\psi\otimes\phi:=\twovec{\alpha}{\beta}\otimes\twovec{\alpha'}{\beta'}\equiv\begin{pmatrix}\alpha\alpha'\\
\alpha\beta'\\
\beta\alpha'\\
\beta\beta'
\end{pmatrix}
\]

where the second $\otimes$ is the Kronecker product and $\alpha\alpha'$
indicates conventional complex multiplication. Note that the basis
relative to which $\psi\otimes\phi$ is represented is the standard
basis for $\C^{4}$ and thus we observe exponential growth in the
size of the representation of an $N$-qubit system. An alternative
notation for state vectors is Dirac notation; for example, for a single

\[
\ket{\psi}\equiv\alpha\kket 0+\beta\kket 1
\]

and a 2-qubit system

\begin{align*}
\ket{\psi}\otimes\ket{\phi} & \equiv\left(\alpha\kket 0+\beta\kket 1\right)\otimes\left(\alpha'\kket 0+\beta'\kket 1\right)\\
& \equiv\alpha\alpha'\kket 0\otimes\kket 0+\alpha\beta'\kket 0\otimes\kket 1+\beta\alpha'\kket 1\otimes\kket 0+\beta\beta'\kket 1\otimes\kket 1\\
& \equiv\alpha\alpha'\kket 0\kket 0+\alpha\beta'\kket 0\kket 1+\beta\alpha'\kket 1\kket 0+\beta\beta'\kket 1\kket 1\\
& \equiv\alpha\alpha'\kket{00}+\alpha\beta'\kket{01}+\beta\alpha'\kket{10}+\beta\beta'\kket{11}\\
& \equiv\alpha\alpha'\kket 0+\alpha\beta'\kket 1+\beta\alpha'\kket 2+\beta\beta'\kket 3
\end{align*}

where in the last line we've used the decimal representation for the
bit strings identifying the basis states. Of particular import for
QC are the entangled or bell states; they correspond
to multi-qubit states, such as

\[
\ket{\xi}=\frac{1}{\sqrt{2}}\ket{00}+\frac{1}{\sqrt{2}}\ket{11}
\]

that cannot be “factored” into component states[$\xi$ cannot be factored because there is no solution to the set
of equations (for $\alpha,\alpha',\beta,\beta'$)

\[
\alpha\alpha'=\frac{1}{\sqrt{2}},\quad\alpha\beta'=0,\quad\beta\alpha'=0,\quad\beta\beta'=\frac{1}{\sqrt{2}}
\]

]. Then, naturally, changes in qubit states are represented as unitary[A matrix $U$ is unitary iff $UU^{\dagger}=U^{\dagger}U=I$, i.e.
it is its own Hermitian conjugate or more simply if it is “self-inverse”.] matrices $U$; for example

\[
\psi'=U\psi=\begin{pmatrix}U_{00} & U_{01}\\
U_{10} & U_{11}
\end{pmatrix}\twovec{\alpha}{\beta}=\twovec{U_{00}\alpha+U_{01}\beta}{U_{10}\alpha+U_{11}\beta}
\]

Matrix representations of transformations on multi-qubit states are
constructed using the Kronecker product on the individual matrix representations;
for example

\[
U\otimes V:=\begin{pmatrix}U_{00}V & U_{01}V\\ U_{10}V & U_{11}V \end{pmatrix}:=\begin{pmatrix}
U_{00} V_{00} & U_{00} V_{01} & U_{01} V_{00} & U_{01} V_{01} \\
U_{00} V_{10} & U_{00} V_{11} & U_{01} V_{10} & U_{01} V_{11} \\
U_{10} V_{00} & U_{10} V_{01} & U_{11} V_{00} & U_{11} V_{01} \\
U_{10} V_{10} & U_{10} V_{11} & U_{11} V_{10} & U_{11} V_{11}
\end{pmatrix}
\]
Here we see again an exponential growth in representation size as
a function of number of qubits.

\begin{equation*}
\psi\begin{cases}\qquad\qquad\quad
\Qcircuit @C=1.0em @R=0.7em {
\lstick{ q_0 := \ket{0} } & \gate{H} & \ctrl{1} & \ctrl{2} & \qw & \qw & \qw
\\
\lstick{ q_1 := \ket{0}  } & \gate{H} & \targ & \qw & \gate{H} & \qw & \qw
\\
\lstick{ q_2 := \ket{0}  } & \gate{H} & \qw & \targ & \gate{H} & \qw & \qw
\\
%\lstick{c:} & {/_{_{1}}} \cw & \cw & \cw & \cw & \cw & \cw\\
\end{cases}
\Rightarrow\; \frac{\ket{000} + \ket{111}}{\sqrt{2}}
\end{equation*}

Quantum Circuit representing 3-qubit $q_{0},q_{1},q_{2}$ entanglement
effected by application of successive Hadamard gates.

As already alluded to, quantum circuits are a formalism for representing
quantum computation in general and algorithms designed for quantum
computers in particular. In the quantum circuit formalism qubit states
are represented by wires and unitary transformations are represented
by gates (see figure <ref>), much
like classical combinational logic circuits might be, though, whereas
combinational logic is “memoryless”[The output of a combinational logic circuit at any time is only a
function of the elements of the circuit and its inputs.], sequences of quantum gates specified by a quantum circuit do in
fact connote the evolution (in time) of the qubits. In addition quantum
gates, as opposed to classical gates, are necessarily reversible and
hence there are no quantum analogs to some classical gates such as
NOT and OR.

§.§ Tensors and Tensor Networks

We quickly define tensors and tensor networks and then move on to
tensor network methods for simulating quantum circuits.

§.§.§ Tensors

One definition of a tensor[There are several more at varying levels of mathematical sophistication.
Chapter 14 of [28] is the standard reference.
Ironically, it is this author's opinion that one should shy away from
physics oriented expositions on tensors.] $T$ is as an element of a tensor product space[The collection of tensor products of elements of the component spaces
quotiented by an equivalence relation.]:

\[
T\in\underbrace{V\otimes\cdots\otimes V}_{p{\text{ copies}}}\otimes\underbrace{V^{*}\otimes\cdots\otimes V^{*}}_{q{\text{ copies}}}
\]

where $V^{*}$ is dual[The dual space to a vector space $V$ is the vector $V^{*}$ consisting
of linear maps $f:V\rightarrow\R$. The dual basis of the dual space
consists of $f_{i}$ such that $f_{i}\left(\mathbf{e}_{i}\right)=\delta_{ij}$.
It is convention to write $f_{i}$ as $\mathbf{e}^{i}$ (note the
superscript index).] to $V$. Then $T$, in effect, acts a multilinear map

\[
{\displaystyle T:\underbrace{V^{*}\times\dots\times V^{*}}_{p{\text{ copies}}}\times\underbrace{V\times\dots\times V}_{q{\text{ copies}}}\rightarrow\R}
\]

by “applying” $p$ elements from $V$ to $p$ elements of $V^{*}$
and $q$ elements from $V^{*}$ to $q$ elements of $V$. Note the
swapping of the orders of $V,V^{*}$ in both the definitions and the
description. $T$'s coordinate basis representation

\begin{equation}
{\displaystyle T\equiv T_{j_{1}\dots j_{q}}^{i_{1}\dots i_{p}}\;\mathbf{e}_{i_{1}}\otimes\cdots\otimes\mathbf{e}_{i_{p}}\otimes\mathbf{e}^{j_{1}}\otimes\cdots\otimes\mathbf{e}^{j_{q}}}\label{eq:tensor_coord}
\end{equation}

is determined by its evaluation on each set of bases

\[
{\displaystyle T_{j_{1}\dots j_{q}}^{i_{1}\dots i_{p}}:=T\left(\mathbf{e}^{i_{1}},\ldots,\mathbf{e}^{i_{p}},\mathbf{e}_{j_{1}},\ldots,\mathbf{e}_{j_{q}}\right)}
\]

The pair $\left(p,q\right)$ is called the type or
valence of $T$ while $\left(p+q\right)$ is the order
of the tensor. Note that we do not use rank to mean either
of these things[The rank of a tensor is the minimum number of distinct basis
tensors necessary to define it; the tensor in eqn. (<ref>)
is in fact a rank 1 tensor. The definition is a generalization of
the rank of a matrix (which, recalling, is the dimension of its column
space, i.e. number of basis elements). Despite this obvious, reasonable
definition for rank, one should be aware that almost all literature
in this area of research uses rank to mean order.]. Furthermore, eqn. (<ref>) in fact represents a
linear sum of basis elements, as it employs Einstein summation convention[Repeated indices in juxtapose position indicate summation $a_{i}b^{i}:=\sum_{i}a[i]b[i]$.].
Note we make liberal use of summation convention in the following
but occasionally use explicit sums when it improves presentation (i.e.
when we would like to emphasize a particular contraction).

There are two important operations on tensors we need to define. Firstly,
we can form the tensor product $Z$ of two tensors $T,W$, of types
$\left(p,q\right),\left(r,s\right)$ respectively, to obtain a tensor
of type $\left(p+r,q+s\right)$:

\begin{align*}
Z & :=T\otimes W\\
& \;=\left(T_{j_{1}\dots j_{q}}^{i_{1}\dots i_{p}}\;\mathbf{e}_{i_{1}}\otimes\cdots\otimes\mathbf{e}_{i_{p}}\otimes{\mathbf{e}}^{j_{1}}\otimes\cdots\otimes{\mathbf{e}}^{j_{q}}\right)\otimes\left(W_{l_{1}\dots l_{s}}^{k_{1}\dots k_{r}}\;\mathbf{e}_{k_{1}}\otimes\cdots\otimes\mathbf{e}_{k_{r}}\otimes{\mathbf{e}}^{l_{1}}\otimes\cdots\otimes{\mathbf{e}}^{l_{s}}\right)\\
& \;=\left(T_{j_{1}\dots j_{q}}^{i_{1}\dots i_{p}}W_{l_{1}\dots l_{s}}^{k_{1}\dots k_{r}}\;\mathbf{e}_{i_{1}}\otimes\cdots\otimes\mathbf{e}_{i_{p}}\otimes\mathbf{e}_{k_{1}}\otimes\cdots\otimes\mathbf{e}_{k_{r}}\otimes{\mathbf{e}}^{j_{1}}\otimes\cdots\otimes{\mathbf{e}}^{j_{q}}\otimes{\mathbf{e}}^{l_{1}}\otimes\cdots\otimes{\mathbf{e}}^{l_{s}}\right)\\
& :=Z_{j_{1}\dots j_{q+s}}^{i_{1}\dots i_{p+r}}\;\mathbf{e}_{i_{1}}\otimes\cdots\otimes\mathbf{e}_{i_{p+r}}\otimes{\mathbf{e}}^{j_{1}}\otimes\cdots\otimes{\mathbf{e}}^{j_{q+s}}
\end{align*}

The attentive reader should notice that the coordinate representation
of two tensors is

Despite it being obvious, its important to note that the tensor product
$Z$ produces a tensor of order $\left(p+r+q+s\right)$, i.e. higher
than either of the operands. On the contrary, tensor contraction
reduces the order of a tensor. We define the contraction $Y$ of type
$\left(a,b\right)$ of a tensor $T$ to be the “pairing” of the
$a$th and $b$th bases:

\begin{align*}
Y & :=T_{j_{1}\dots j_{q}}^{i_{1}\dots i_{p}}\;\mathbf{e}_{i_{1}}\otimes\cdots\otimes\mathbf{e}_{i_{a-1}}\otimes\mathbf{e}_{i_{a+1}}\otimes\cdots\otimes\mathbf{e}_{i_{p}}\otimes\left(\mathbf{e}_{i_{a}}\cdot\mathbf{e}^{j_{b}}\right)\otimes\mathbf{e}^{j_{1}}\otimes\cdots\otimes\mathbf{e}^{j_{b-1}}\otimes\mathbf{e}^{j_{b+1}}\otimes\cdots\otimes\mathbf{e}^{j_{q}}\\
& \;=T_{j_{1}\dots j_{q}}^{i_{1}\dots i_{p}}\delta_{i_{a}}^{j_{b}}\;\mathbf{e}_{i_{1}}\otimes\cdots\otimes\mathbf{e}_{i_{a-1}}\otimes\mathbf{e}_{i_{a+1}}\otimes\cdots\otimes\mathbf{e}_{i_{p}}\otimes\mathbf{e}^{j_{1}}\otimes\cdots\otimes\mathbf{e}^{j_{b-1}}\otimes\mathbf{e}^{j_{b+1}}\otimes\cdots\otimes\mathbf{e}^{j_{q}}\qquad\left(\text{since }\mathbf{e}_{i}\cdot\mathbf{e}^{j}=\delta_{i}^{j}\right)\\
& \;=\sum_{j_{b}}T_{j_{1}\dots j_{b}\dots j_{q}}^{i_{1}\dots j_{b}\dots i_{p}}\;\mathbf{e}_{i_{1}}\otimes\cdots\otimes\mathbf{e}_{i_{a-1}}\otimes\mathbf{e}_{i_{a+1}}\otimes\cdots\otimes\mathbf{e}_{i_{p}}\otimes\mathbf{e}^{j_{1}}\otimes\cdots\otimes\mathbf{e}^{j_{b-1}}\otimes\mathbf{e}^{j_{b+1}}\otimes\cdots\otimes\mathbf{e}^{j_{q}}\\
& :=Y_{j_{1}\dots i_{b-1}i_{b+1}\dots j_{q}}^{i_{1}\dots i_{a-1}i_{a+1}\dots i_{p}}\;\mathbf{e}_{i_{1}}\otimes\cdots\otimes\mathbf{e}_{i_{a-1}}\otimes\mathbf{e}_{i_{a+1}}\otimes\cdots\otimes\mathbf{e}_{i_{p}}\otimes\mathbf{e}^{j_{1}}\otimes\cdots\otimes\mathbf{e}^{j_{b-1}}\otimes\mathbf{e}^{j_{b+1}}\otimes\cdots\otimes\mathbf{e}^{j_{q}}
\end{align*}

where $\left(\cdot\right)$ means inner product. Notice that the order
of $Y$ is $\left(p-1,q-1\right)$. Finally notice that we can omit
writing out bases and just manipulate coordinates. We shall do as
such when it simplifies presentation.

As mentioned in the introduction, matrices can be represented as tensors;
for example, the two dimensional $N\times N$ matrix $M$ is taken
to be a tensor of type $\left(1,1\right)$ with basis representation

\[
{\displaystyle M\equiv M_{j}^{i}\;\mathbf{e}_{i}\otimes\mathbf{e}^{j}}
\]

where upper indices correspond to the row index and lower indices
correspond to the column index of the conventional matrix representation
and both range from $1$ to $N$. The attentive reader will notice
that the coordinate representation of the tensor product for type
$\left(1,1\right)$ tensors is exactly the Kronecker product for matrices.
Similarly, tensor contraction for type $\left(1,1\right)$ tensors
is the familiar matrix trace:

\[
\]

More usefully, we can express matrix-vector multiplication in terms
of tensor contraction; let

\[
\mathbf{x}:=\twovec{x^{1}}{x^{2}}\equiv x^{1}\mathbf{e}_{1}+x^{2}\mathbf{e}_{2}\equiv x^{i}\mathbf{e}_{i}
\]

where we switch to valence index notation in the column vector for
closer affinity with tensor notation. Then it must be the case that

\[
\mathbf{y}=M\mathbf{x}=\left(M_{j}^{i}\mathbf{e}_{i}\otimes\mathbf{e}^{j}\right)\left(x^{k}\mathbf{e}_{k}\right)=\left(M_{j}^{i}x^{k}\mathbf{e}_{i}\right)\left(\mathbf{e}^{j}\cdot\mathbf{e}_{k}\right)=M_{j}^{i}x^{k}\delta_{k}^{j}\mathbf{e}_{i}=M_{j}^{i}x^{j}\mathbf{e}_{i}
\]

Letting $y^{i}:=M_{j}^{i}x^{j}$ we recognize conventional matrix-vector
multiplication. Employing tensor contraction in this way extends to
matrix-matrix multiplication (and tensor composition more broadly);
for two type $\left(1,1\right)$ tensors $M,L$ we can form the type
$\left(1,1\right)$ tensor $Z$ corresponding to matrix product $M\cdot L$
of $N\times N$ by first taking the tensor product

\[
\]

The attentive reader will notice that the coordinate representation
of two tensors is exactly the Kronecker product of two matrices. Then
contracting along the off diagonal

\begin{equation}
\end{equation}

One can confirm that this is indeed conventional matrix multipliation
of two $N\times N$ matrices. In general, stated simply, when contracting
indices of a tensor product, contraction can be understood to be a
sum over shared indices.

§.§.§ Tensor Networks

[dashed]  picgrid
(0.3,4.7) node$\mathbf{S}\boldsymbol{'}$;
(0.35,0.35) rectangle (4.65,4.65);
[rotate around=45:(45:2.5),scale=sqrt(2),
shift=(0,-0.75),transform shape]  picgrid=nx=3,ny=3
Projected Entangled Pair State (PEPS) for a 3 × 3 lattice
with open boundary conditions.

[Contraction ]

[State vector representation]

[State vector factorization]

Tensor network diagrammatic contraction. [8]

Tensor networks (TNs) are a way to factor tensors with large orders
into networks of tensors with lower orders; since the number of parameters
a tensor consists of is exponential in the order of the tensor, smaller
order tensors are much preferrable computationally. They were first
used to study ground states of one dimensional quantum many-body systems
[36] but have since been applied in other areas
(such as machine learning [13]). TNs lend
themselves to a diagrammatic representation which can be used to reason
about such factorizations (figure <ref>).
We will primarily be interested in TNs as a means to factoring the
state-vector of an $N$-qubit system (see figure <ref>)

\begin{equation}
\ket{\psi}:=\sum_{i_{1}i_{2}\dots i_{N}}C^{i_{1}i_{2}\dots i_{N}}\ket{i_{1}}\ket{i_{2}}\cdots\ket{i_{N}}\label{eq:state_vector}
\end{equation}

for which its common to propose an ansatz factorizations:

* Matrix Product States (MPS) [19], which yields

\[
C^{i_{1}i_{2}\dots i_{N}}\equiv A_{j_{1}}^{i_{1}}A_{j_{2}}^{i_{1}j_{1}}\cdots A_{j_{N-1}}^{i_{N-1}j_{N-2}}A^{i_{N}j_{N-1}}
\]

where $j$ are called bond indices. If each index $i$ has
dimension $d$ (i.e. takes on values 1 to $d$) then $C$ is specified
by $d^{N}$ parameters and can always be represented by an MPS factorization
$Ndm^{2}$ parameters, where $m:=d^{N/2}$ is the bond dimension.
While for this naive representation $d^{N}<Ndm^{2}$, in practice
$m$ is fixed to some moderate size such that $d^{N}>Ndm^{2}$ and
the MPS factorization functions as an approximation.

* Projected Entangled Pair States (PEPS) [32],
which is a generalization of MPS to higher spatial dimensions, i.e.
TNs that correspond to lattices of contractions of tensors, which
themselves represent pairwise entangled quantum systems. Naturally,
such a series of contractions doesn't lend itself to being expressed
in traditional notation and thus we observe the power of tensor network
diagrams (see PEPS in figure <ref>).

* Tree Tensor Networks (TTN) [30], a
further generalization where tensors are entangled (and therefore
contracted) hierarchically. In fact TTNs bear the closest
resemblance to quantum circuits.

* Multi-scale Entanglement Renormalization Ansatz (MERA) [33],
a specific type of TTN where the tensors are alternatingly unitaries
and isometries[A tensor, seen as a multlinear map, that preserves distances under
the ambient distance metric.].

§.§.§ TNs for Simulating Quantum Circuits

Factoring eqn. (<ref>) is only the first step to
successfully simulating a quantum circuit. By representing some final
state as a tensor as well, and contracting across all indices (called
contracting the network), we can calculate the amplitude
for that particular state. Since tensor contraction is associative[This can be observed by noting that summing is an associative operation
(or by analog with matrix-matrix multiplication).], the order in which tensors are actually contracted is a “hyperparameter”
of TN methods; finding the optimal contraction order, with respect
to accuracy (assuming some approximation has been made in constructing
the factorization), compute time, and memory requirements is critical.

In particular we focus on contraction orders for TTNs as they most
closely resemble quantum circuits. For a TTN consisting of $N$ tensors
(corresponding to $N$ gates) with maximum order $p$, worst case,
we can see that contraction time takes $O\left(N\exp\left(O\left(p\right)\right)\right)$
since, in general, contracting across all indices of a pair of tensors
is exponential in their orders[Consider contracting two $\left(1,1\right)$ tensors (as in eqn. (<ref>)),
i.e. two order 2 tensors, which effectively is matrix multplication
followed by trace. The complexity of this contraction is then $O\left(N^{2+1}+N\right)\equiv O\left(\exp\left(2\log N\right)\left(1+N\right)\right)$
(where $N$ here is the characteristic dimension of the matrix). Assuming
the ranges of all tensor indices is the same (i.e. $N$ is constant
across all tensors), for example $N=2$ as in the case of matrices
derived of unitary transformations operating on single qubits, we
recover the stated complexity. ]. Markov et al. [21] showed that there in fact
exists a contraction ordering which results in a contraction time
of $O\left(N^{O\left(1\right)}\exp\left(O\left(\operatorname{tw}\left(G^{L}\right)\right)\right)\right)$
where $G^{L}$ is the line graph[A line graph captures edge adjacency; given a graph $G$,
$G^{L}$ is defined such that each edge of $G$ corresponds to a vertex
of $G^{L}$ and two vertices are are connected in $G^{L}$ if the
edges in $G$ that they correspond to are adjacent on the same vertex
(in $G$). ] of the tensor network and $\operatorname{tw}\left(G^{L}\right)$
is the tree-width[A tree decomposition of a graph $G$ is a tree $T$ and a
mapping from the vertices of $G$ into “bags” that satisfy the
following properties

* Each vertex must appear in at least one bag.

* For each edge in $G$, at least one bag must contain both of the vertices
it is adjacent on.

* All bags containing a given vertex in $G$ must be connected in $T$.

The width $w$ of a tree decomposition is the cardinality of the largest
bag (minus one). Finally the tree-width of $G$ is the minimum
width over all possible tree decompositions. Intuitively, a graph
has low tree-width if it can be constructed by joining small graphs
together into a tree.] of $G^{L}$.

Thus, in general, we can simulate arbitrary quantum circuits by

* Constructing the TTN (with graph $G$) that corresponds to the quantum

* Computing a tree-decomposition of $G^{L}$ with width $\operatorname{tw}\left(G^{L}\right)$.

* Finding a contraction ordering[Markov et al. also show that one can recover a contraction order from
the tree-decomposition of $G^{L}$ in polynomial time.] for the TTN and fully contracting to a single, real, number (i.e.
order 0 tensor).

Note that, in general, computing the tree-width of a graph is NP-hard
[3] and Markov et al.'s results rely on
generating any tree-decomposition that has a tree-width that
suits their needs (i.e. overall runtime complexity).

For quantum circuits consisting of many few qubit gates this technique
produces a much more (runtime) efficient evaluation of the circuit;
indeed Markov et al. further show that any TTN corresponding to a
quantum circuit with $N$ gates, where the number of gates that act
on any pair of qubits is bounded by $r$, has contraction time $O\left(N^{O\left(1\right)}\exp\left(O\left(r\right)\right)\right)$.

Markov et al.'s results are not tight; their construction finds some
tree-decomposition with the correct corresponding tensor contraction
order that suits their aim (overall runtime complexity of the translation
from quantum circuit to TTN and the ultimate contraction). In reality
there are often contraction orders that are much more space and runtime
efficient. Though, in general problem is NP-hard [3],
for particular TTNs (corresponding to circuits) there are heuristics,
such as non-adjacent contractions [27],
that produce more efficient orders. Alternatively, randomized search
and Bayesian optimization can be used to identify efficient contraction
orders [15, 11].

§.§ FPGAs

FPGA floorplan diagram [24].

A field-programmable gate array (FPGA) is a device designed to be
configured by a user into various arrangements of (classical) gates
and memory. FPGAs consist of arrays (hence the name) of configurable
logic blocks (CLBs), static ram (SRAM), and programmable busses that
connect CLBs and SRAM into various topologies (see figure <ref>).
The CLBs typically contain arithmetic units (such as adders, multipliers,
and accumulators) and lookup tables (LUTs), that can be programmed
to represent truth tables for many boolean functions. Using hardware
description languages (such as VHDL or Verilog) designers specify
modules and compose them into circuits (also known as a dataflows)
that perform arbitrary computation. These circuits then go through
a place and route procedure before ultimately being instantiated
on the FPGA as processing elements (PEs) and connections
between PEs.

While modules consisting purely of combinational logic compute their
outputs at the stated clock speed of the FPGA, inevitably I/O (i.e.
fetching data from memory) interleaved with such modules (otherwise
arranged into a pipeline architecture) creates pipeline stalls. Thus,
it's essential that FPGA designs are as compute bound as possible
(rather than I/O bound). In particular, we explore I/O minimal generalized
matrix multiplication (GEMM) [7] and other
systolic array architectures [34, 12].
A systolic architecture [20] is a gridded, pipelined, array
of PEs that processes data as the data flows[The relationship to cardiovascular “systolic” is in association
with the flow of data into the array, akin to how blood flows through
the veins into the human heart.] through the array. Crucially, a systolic architecture propagates
partial results as well as input data through the pipeline (see figure
<ref>). Systolic arrays are particularly suited
for I/O efficient matrix multiplication owing to the pipelining of
inputs (see figure <ref>).

[Gemmini systolic array architecture. The processing elements (PEs)
are either of type Weight Stationary (WS) or Output Stationary (OS).
[Systolic array architecture implementing matrix multiplication. Input
matrices $A$ and $B$ stream by to produce output matrix $C$ via
successive mutiply-accumulate (MAC) operations. Note that $C$ remains
in the processing elements (i.e. this is a diagram of an OS architecture).
[38].]Systolic arrays

One remaining hurdle to simulating quantum compuations (i.e. carrying
out tensor contractions) on FPGAs is SRAM. The standard remediation
is to perform arithmetic with reduced precision[Germaine to this issue is the fact that arithmetic on FPGAs is typically
performed in fixed precision (via an integer representation), owing
to higher compute cost incurred for floating point arithmetic.]. There is evidence that suggests that simulations of quantum circuits,
of varying depths [6], are robust to reduced
precision computation as long as that loss of precision is uncorrelated
[37] i.e. insofar as it can be treated
as uncorrelated noise.

§ IMPLEMENTATION

\begin{equation*}
\Qcircuit @C=0.5em @R=1.0em @!R { 	 	\lstick{ {q}_{0} :  } & \gate{H} & \ctrl{1} & \gate{R_z} & \gate{R_x} & \qw & \ctrl{1} & \gate{R_z} & \gate{R_x} & \qw & \ctrl{1} & \gate{R_z} & \gate{R_x} & \qw & \ctrl{1} & \gate{R_z} & \gate{R_x} & \gate{H} & \qw & \qw & \qw && \cdots\\ 	 	\lstick{ {q}_{1} :  } & \gate{H} & \targ & \gate{R_z} & \ctrl{1} & \gate{R_x} & \targ & \gate{R_z} & \ctrl{1} & \gate{R_x} & \targ & \gate{R_z} & \ctrl{1} & \gate{R_x} & \targ & \gate{R_z} & \ctrl{1} & \gate{R_x} & \gate{H} & \qw & \qw && \cdots\\ 	 	\lstick{ {q}_{2} :  } & \gate{H} & \ctrl{1} & \gate{R_z} & \control\qw & \gate{R_x} & \ctrl{1} & \gate{R_z} & \control\qw & \gate{R_x} & \ctrl{1} & \gate{R_z} & \control\qw & \gate{R_x} & \ctrl{1} & \gate{R_z} & \control\qw & \gate{R_x} & \gate{H} & \qw & \qw && \cdots\\ 	 	\lstick{ {q}_{3} :  } & \gate{H} & \targ & \gate{R_z} & \ctrl{1} & \gate{R_x} & \targ & \gate{R_z} & \ctrl{1} & \gate{R_x} & \targ & \gate{R_z} & \ctrl{1} & \gate{R_x} & \targ & \gate{R_z} & \ctrl{1} & \gate{R_x} & \gate{H} & \qw & \qw && \cdots\\ 	 	\lstick{ {q}_{4} :  } & \gate{H} & \ctrl{1} & \gate{R_z} & \control\qw & \gate{R_x} & \ctrl{1} & \gate{R_z} & \control\qw & \gate{R_x} & \ctrl{1} & \gate{R_z} & \control\qw & \gate{R_x} & \ctrl{1} & \gate{R_z} & \control\qw & \gate{R_x} & \gate{H} & \qw & \qw && \cdots\\ 	 	\lstick{ {q}_{5} :  } & \gate{H} & \targ & \gate{R_z} & \gate{R_x} & \qw & \targ & \gate{R_z} & \gate{R_x} & \qw & \targ & \gate{R_z} & \gate{R_x} & \qw & \targ & \gate{R_z} & \gate{R_x} & \gate{H} & \qw & \qw & \qw && \cdots\\ 	 	\lstick{\vdots} & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots && \vdots \\ \\	 }
\end{cases}}_{k}
\end{equation*}

Test quantum circuit for $n$ qubits and rounds $k$.

We use  [14] to specify quantum circuits
and generate TNs therefrom. In particular we simulate circuits for
various $n$ qubits and rounds $k$ where each consists of
alternating qubit couplings according of the form in figure <ref>.
We also use Bayesian parameter optimization (BPO) [15]
to find tensor contraction orders and compare against naive greedy
search. Note that for both strategy we set a timeout of 600 seconds.
We then deploy the contraction strategy that produces the fewest number
of tensor contractions balanced against the orders of intermediate
tensors[This choice was purely due to platform constraints in that large intermediate
tensors could not be effectively simulated.]. In order to expedite the process of deploying we precompute some
first few tensor contraction such that all tensors deployed to the
FPGA are square and congruent (i.e. all of the same dimensions). For
tensors of order greater than $\left(1,1\right)$ (i.e. tensors that
are not matrices) we transform them into $\left(1,1\right)$ tensors
by taking the Kronecker product of all component $\left(1,1\right)$
tensors; to be precise we perform the following operation on the $\left(p,q\right)$
mats = [t[idx] for idx in np.ndindex(t.shape[:-2])]
block_mat = block_diag(*qubit_mats)

where builds a block diagonal matrix
of its arguments. All of our code has been made available on GitHub[<https://github.com/makslevental/fpga_stuff/> on the

For deploying circuits to FPGAs we use Chisel [4] as a
HDL, by way of an adaptation of the Gemmini systolic array generator
[12]. Notably we experiment with using Gemmini
as an accelerator (i.e. fully parameterizing the weights/entries of
the matrices) and “hardcoding” certain gates/tensors. One possible
advantage of the latter approach over the former is a reduction in
loads from memory for the weights. The success of the chosen approach
depends heavily on whether certain sequences of fixed gates can actually
be pipelined or alternatively deployed in toto to the FPGA. We hypothesize
that this might depend on the depth and gate count of the circuits/TNs.
See figure <ref> for the netlists corresponding
to our systolic array and matmul implementations. Note that (complex)
arithmetic was done in 32 bit fixed precision for both implementations,
with 28 bits allocated behind the binary point.

One challenge we faced was in deploying to real hardware[A challenge not unfamiliar to the seasoned QC researcher.];
unfortunately time and administrative challenges[We were not able to get allocations on CHI@TACC in a timely fashion
(the issue is ongoing...).] prevented us from actually deploying to real FPGAs. As a substitute
we used the well-known and trusted Verilog simulator[It is simulations all the way down.]
Verilator[<https://www.veripool.org/verilator/>], which
transpiles Verilog (which Chisel generates) to a cycle-accurate model
in C++. We then executed this model to collect proxy measurements.
Note that for certain configurations (generally those with high qubit
and round count) we could successfully simulate due to memory constraints
on the workstation running the Verilator produced model.

[Output stationary systolic array processing element. ]

[$4\times4$ matrix multiplication systolic array implementation. ]

[Naive $4\times4$ matrix multiplication implementation. ]

FPGA implementations

§ EVALUATION

[Runtimes for computing BPO contraction strategy.].5![Runtimes for computing greedy contraction strategy.].5!Runtimes for computing contraction strategies for circuits
for various $n$ qubits and rounds $k$. Color scale represents time
and marker size represents, qualitatively, the number of tensors in
the tensor network corresponding to the quantum circuit. Note that
absent markers correspond to searches that didn't converge.

[Runtimes for contracting with systolic array.].5![Runtimes for contracting with matmul.].5!Runtimes for test circuits for various $n$ qubits and rounds $k$.
Note that absent markers correspond to contractions that could not
be simulated due to memory constraints imposed by using the Verilator
model rather than an actual FPGA.

We perform two sets of evaluations. Even though it was not the central
goal of our exploration we first compare the time required to compute
a tensor contraction strategy across $n$ qubits and rounds $k$ for
the greedy search strategy and the BPO search strategy. We then address
our central goal in comparing the actual runtime for performing the
discovered tensor contraction on both the systolic array implementation
(see fig. <ref>) and the “hardcoded”
naive matmul (see fig. <ref>).

Some interesting things to note searching for contractions: computing
(not evaluating) the optimal contraction strategy (i.e. using BPO)
is generally more performant that greedy search (see figures <ref>,
<ref>). The likely reason for this is that
BPO converges more quickly and more efficiently searches the space
of possible contraction orders than greedy search (which greedily
optimizes some surrogate objective). Also note that, in fact, for
certain configurations greedy didn't converge within the timeout.

Regarding the differences in the evaluation times of the contraction
orders (figures <ref>, <ref>)
it's clear that the systolic array implementation is more performant
in terms of both memory requirements and runtime. This is paradoxically
both obvious and suprising. As already mentioned, one expects systolic
arrays to have improved performance relative to naive matrix multiplication
for streaming data (and indeed, in general, they do) but for this
use case (where all matrix elements are known at deploy time) one
also expects that latency incurred by pipelining would offset that
performance improvement. One hypothesis for this is that the difference
is an artifact of simulating the FPGA implementations insofar as simulating
a more densely connected FPGA implementation (see the differences
between <ref> and <ref>)
is more compute intensive, especially with respect to heap allocations
(since systolic arrays incur more loads from memory). To corroborate
this hypothesis we used Intel's Quartus EDA[Electronic design automation.]
tool to synthesize and place and route (for an Arria II GX). Indeed
the naive implementation an order of magnitude (sometimes several)
of each type of resource (see table <ref>).

Systolic     NaiveLUTs     512     110,074Registers     896     2,048Pins     642     3,074DSPs     32     232

Synthesis and place and route for an Arria II GX for both systolic
arrays and naive matmul.

§ CONCLUSION

We explored tensor networks deployed to FPGAs as a mean of accelerating
simulations of quantum circuits. In order to accomplish this goal
we expressed tensor contraction as sequences of matrix multiplication
and implemented two different matrix multiplication FPGA designs:
systolic arrays, which operate on streaming matrix elements and naive
matrix multiplication, which wholesale instantiates all the necessary
MAC operations. In order to choose the contraction orders we used
an “off the shelf” library which searches for a suitable contraction
by either performing greedy search or Bayesian optimization. We compared
the performance of both the contraction search strategy and each contraction
evaluation implementation. Unfortunately we were unable to obtain
acecss to FPGA devices and thus we made due with cycle-accurate simulations.
Results for both comparison were generally in agreement with intuition:
BPO converged to a contraction order more effectively (more quickly
and more robustly) than greedy search and systolic arrays evaluated
the contraction more efficiently than naive matrix multiplication.

Possible future work includes actually deploying to real FPGAs and
then further comparing performance to the simulations performed here.
Another particularly interesting research direction is the tangential
problem of discovering optimal tensor contraction orders. Finding
such tensor contraction orders is ultimately a combinatorial optimization
problem. It occurs to us that possibly a deep learning approach could
be effective. Recently there has been work on RL for combinatorial
optimization[5] and MCTS
for combinatorial optimization[2] that could, possibly,
be adapted to this problem in a straightforward fashion.

[1]
K. Abdelouahab, M. Pelcat, J. Sérot, C. Bourrasset, and F. Berry.
Tactics to directly map cnn graphs on embedded fpgas.
IEEE Embedded Systems Letters, 9(4):113–116, 2017.

[2]
Kenshin Abe, Zijian Xu, Issei Sato, and Masashi Sugiyama.
Solving np-hard problems on graphs with extended alphago zero, 2020.

[3]
Stefan Arnborg, Derek G. Corneil, and Andrzej Proskurowski.
Complexity of finding embeddings in a k-tree.
SIAM JOURNAL OF DISCRETE MATHEMATICS, 8(2):277–284, 1987.

[4]
Jonathan Bachrach, Huy Vo, Brian Richards, Yunsup Lee, Andrew Waterman, Rimas
Avižienis, John Wawrzynek, and Krste Asanović.
Chisel: Constructing hardware in a scala embedded language.
In DAC Design Automation Conference 2012, pages 1212–1221,

[5]
Thomas Barrett, William Clements, Jakob Foerster, and Alex Lvovsky.
Exploratory combinatorial optimization with reinforcement learning.
Proceedings of the AAAI Conference on Artificial Intelligence,
34(04):3243–3250, Apr. 2020.

[6]
Santiago I. Betelu.
The limits of quantum circuit simulation with low precision
arithmetic, 2020.

[7]
Johannes de Fine Licht, Grzegorz Kwasniewski, and Torsten Hoefler.
Flexible communication avoiding matrix multiplication on fpga with
high-level synthesis.
In Proceedings of the 2020 ACM/SIGDA International Symposium on
Field-Programmable Gate Arrays, FPGA '20, page 244–254, New York, NY, USA,
2020. Association for Computing Machinery.

[8]
Glen Evenbly.
Tensortrace: an application to contract tensor networks, 2019.

[9]
Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.
A quantum approximate optimization algorithm, 2014.

[10]
Richard P Feynman.
Simulating physics with computers.
International journal of theoretical physics, 21(6/7):467–488,

[11]
E. Schuyler Fried, Nicolas P. D. Sawaya, Yudong Cao, Ian D. Kivlichan,
Jhonathan Romero, and Alán Aspuru-Guzik.
qtorch: The quantum tensor contraction handler.
PLOS ONE, 13(12):e0208510, Dec 2018.

[12]
Hasan Genc, Ameer Haj-Ali, Vighnesh Iyer, Alon Amid, Howard Mao, John Wright,
Colin Schmidt, Jerry Zhao, Albert Ou, Max Banister, Yakun Sophia Shao,
Borivoje Nikolic, Ion Stoica, and Krste Asanovic.
Gemmini: An agile systolic array generator enabling systematic
evaluations of deep-learning architectures, 2019.

[13]
Ivan Glasser, Ryan Sweke, Nicola Pancotti, Jens Eisert, and J. Ignacio Cirac.
Expressive power of tensor-network factorizations for probabilistic
modeling, with applications from hidden markov models to quantum machine
learning, 2019.

[14]
Johnnie Gray.
quimb: A python package for quantum information and many-body
Journal of Open Source Software, 3(29):819, 2018.

[15]
Johnnie Gray and Stefanos Kourtis.
Hyper-optimized tensor network contraction.
Quantum, 5:410, Mar 2021.

[16]
Lov K. Grover.
A fast quantum mechanical algorithm for database search.
In Proceedings of the Twenty-Eighth Annual ACM Symposium on
Theory of Computing, STOC '96, page 212–219, New York, NY, USA, 1996.
Association for Computing Machinery.

[17]
Abhijith J., Adetokunbo Adedoyin, John Ambrosiano, Petr Anisimov, Andreas
Bärtschi, William Casper, Gopinath Chennupati, Carleton Coffrin, Hristo
Djidjev, David Gunter, Satish Karra, Nathan Lemons, Shizeng Lin, Alexander
Malyzhenkov, David Mascarenas, Susan Mniszewski, Balu Nadiga, Daniel
O'Malley, Diane Oyen, Scott Pakin, Lakshman Prasad, Randy Roberts, Phillip
Romero, Nandakishore Santhi, Nikolai Sinitsyn, Pieter J. Swart, James G.
Wendelberger, Boram Yoon, Richard Zamora, Wei Zhu, Stephan Eidenbenz,
Patrick J. Coles, Marc Vuffray, and Andrey Y. Lokhov.
Quantum algorithm implementations for beginners, 2020.

[18]
Liancheng Jia, Liqiang Lu, Xuechao Wei, and Yun Liang.
Generating systolic array accelerators with reusable blocks.
IEEE Micro, 40(4):85–92, 2020.

[19]
A Klümper, A Schadschneider, and J Zittartz.
Matrix product ground states for one-dimensional spin-1 quantum
Europhysics Letters (EPL), 24(4):293–297, Nov 1993.

[20]
Why systolic architectures?
Computer, 15(1):37–46, 1982.

[21]
Igor L. Markov and Yaoyun Shi.
Simulating quantum computation by contracting tensor networks.
SIAM J. Comput., 38(3):963–981, June 2008.

[22]
Alexander McCaskey, Eugene Dumitrescu, Mengsu Chen, Dmitry Lyakh, and Travis
Validating quantum-classical programming models with tensor network
PLOS ONE, 13(12):e0206704, Dec 2018.

[23]
Thierry Moreau, Tianqi Chen, Luis Vega, Jared Roesch, Eddie Yan, Lianmin Zheng,
Josh Fromm, Ziheng Jiang, Luis Ceze, Carlos Guestrin, and Arvind
A hardware–software blueprint for flexible deep learning
IEEE Micro, 39(5):8–16, 2019.

[24]
K. E. Murray, M. A. Elgammal, V. Betz, T. Ansell, K. Rothman, and A. Comodi.
Symbiflow and vpr: An open-source design flow for commercial and
novel fpgas.
IEEE Micro, 40(04):49–57, jul 2020.

[25]
Eriko Nurvitadhi, Ganesh Venkatesh, Jaewoong Sim, Debbie Marr, Randy Huang,
Jason Ong Gee Hock, Yeong Tat Liew, Krishnan Srivatsan, Duncan Moss, Suchit
Subhaschandra, and Guy Boudoukh.
Can fpgas beat gpus in accelerating next-generation deep neural
In Proceedings of the 2017 ACM/SIGDA International Symposium on
Field-Programmable Gate Arrays, FPGA '17, page 5–14, New York, NY, USA,
2017. Association for Computing Machinery.

[26]
National Academies of Sciences Engineering and Medicine.
Quantum Computing: Progress and Prospects.
The National Academies Press, Washington, DC, 2019.

[27]
Edwin Pednault, John A. Gunnels, Giacomo Nannicini, Lior Horesh, Thomas
Magerlein, Edgar Solomonik, Erik W. Draeger, Eric T. Holland, and Robert
Pareto-efficient quantum circuit simulation using tensor contraction
deferral, 2020.

[28]
S. Roman.
Advanced Linear Algebra.
Graduate Texts in Mathematics. Springer New York, 2007.

[29]
Justin Sanchez, Nasim Soltani, Pratik Kulkarni, Ramachandra Vikas Chamarthi,
and Hamed Tabkhi.
A reconfigurable streaming processor for real-time low-power
execution of convolutional neural networks at the edge.
In Shijun Liu, Bedir Tekinerdogan, Mikio Aoyama, and Liang-Jie Zhang,
editors, Edge Computing – EDGE 2018, pages 49–64, Cham, 2018.
Springer International Publishing.

[30]
Y.-Y. Shi, L.-M. Duan, and G. Vidal.
Classical simulation of quantum many-body systems with a tree tensor
Phys. Rev. A, 74:022320, Aug 2006.

[31]
P.W. Shor.
Algorithms for quantum computation: discrete logarithms and
In Proceedings 35th Annual Symposium on Foundations of Computer
Science, pages 124–134, 1994.

[32]
F. Verstraete and J. I. Cirac.
Renormalization algorithms for quantum-many body systems in two and
higher dimensions.
7 2004.

[33]
G. Vidal.
Entanglement renormalization.
Phys. Rev. Lett., 99:220405, Nov 2007.

[34]
Jie Wang, Licheng Guo, and Jason Cong.
Autosa: A polyhedral compiler for high-performance systolic arrays on
In The 2021 ACM/SIGDA International Symposium on
Field-Programmable Gate Arrays, FPGA '21, page 93–104, New York, NY, USA,
2021. Association for Computing Machinery.

[35]
Pete Warden.
Why gemm is at the heart of deep learning, Apr 2015.

[36]
Steven R. White.
Density matrix formulation for quantum renormalization groups.
Phys. Rev. Lett., 69:2863–2866, Nov 1992.

[37]
Xin-Chuan Wu, Sheng Di, Emma Maitreyee Dasgupta, Franck Cappello, Hal Finkel,
Yuri Alexeev, and Frederic T. Chong.
Full-state quantum circuit simulation by using data compression.
In Proceedings of the International Conference for High
Performance Computing, Networking, Storage and Analysis, SC '19, New York,
NY, USA, 2019. Association for Computing Machinery.

[38]
Zhijie Yang, Lei Wang, Dong Ding, Xiangyu Zhang, Yu Deng, Shiming Li, and Qiang
Systolic array based accelerator and algorithm mapping for deep
learning algorithms.
In Feng Zhang, Jidong Zhai, Marc Snir, Hai Jin, Hironori Kasahara,
and Mateo Valero, editors, Network and Parallel Computing, pages
153–158, Cham, 2018. Springer International Publishing.

[39]
Han-Sen Zhong, Hui Wang, Yu-Hao Deng, Ming-Cheng Chen, Li-Chao Peng, Yi-Han
Luo, Jian Qin, Dian Wu, Xing Ding, Yi Hu, Peng Hu, Xiao-Yan Yang, Wei-Jun
Zhang, Hao Li, Yuxuan Li, Xiao Jiang, Lin Gan, Guangwen Yang, Lixing You,
Zhen Wang, Li Li, Nai-Le Liu, Chao-Yang Lu, and Jian-Wei Pan.
Quantum computational advantage using photons.
Science, 370(6523):1460–1463, 2020.