# Optimization of DSP Applications Using Parameterized Error Models for Low
Power Approximate Adders

Celia Dharmaraj, Vinita Vasudevan, and Nitin Chandrachoodan The authors are
with the Department of Electrical Engineering, Indian Institute of Technology
Madras, India.
E-mail<EMAIL_ADDRESS>
###### Abstract

Approximate circuit design has gained significance in recent years targeting
error tolerant applications. In this paper, we first demonstrate that the
commonly used assumption that the inputs to the adder are uniformly
distributed results in an inaccurate prediction of error statistics for multi-
level circuits. To overcome this problem, we derive parameterized error models
that can be used within any optimization framework in order to optimize the
number of approximate bits. We also show that in order to accurately compute
the MSE, the optimization framework needs to take into account not just the
functionality of the adder, but also its position in the circuit,
functionality of its parents and the number of approximate bits in the parent
blocks. We demonstrate a significant improvement of accuracy in the prediction
of the noise power of DSP systems containing approximate adders.

###### Index Terms:

Approximate computing, error model, static probability, optimization, noise
power.

##  I Introduction

Approximate computing is widely used in signal and image processing
applications to obtain improvements in power and/or speed while maintaining
the required accuracy. Adders are the basic building blocks in these
applications and a typical implementation has a large number of adders. A
variety of approximate adders have been proposed in the literature, with
different levels of trade-offs between accuracy and performance. These adders
can be classified as low-latency [1] (and references therein) and low-power
approximate adders (LPAA) [2, 3, 4, 5, 6, 7]. Low power implementations of DSP
systems using LPAAs optimize the power consumption by maximizing the number of
approximate bits in each adder for a given accuracy.

In the literature, multiple approaches have been proposed to find optimal
approximation levels for adders used in low power implementations. An
approximate Finite Impulse Response (FIR) filter is designed by fixing the
level of approximation of the adders using Monte-Carlo (MC) simulations in
[8]. Approximate mirror adder-5 (AMA-5) [4] modeled assuming uniformly
distributed inputs are used in a 2D Discrete Cosine Transform (DCT) module
constructed using 1D DCT blocks in [9] and the optimization problem is solved
using a mixed integer non-linear problem solver. Cartesian Genetic Programming
(CGP) is used to design various approximate implementations of four point 1D
DCT in [10]. An expression for variance of error of AMA 1-5 adders [4] and
Lower-part-OR adder (LOA) [5] is obtained in [11] empirically by regression
assuming uniform inputs, and heuristics are used to solve the approximation-
level optimization problem. In [12], AMA 1-5 adders and Transmission Gate
based Approximate adders TGA I-II [13] are considered. An expression for mean
square error (MSE) is obtained assuming that the distribution of inputs and
error are uniform. This is then used in a Lagrange multiplier based
optimization approach.

All of the above previous works on optimization use error metrics based on
uniformly distributed inputs. Moreover, the same error model is used for all
adders in the circuit. To verify the validity of these assumptions, we found
the MSE in an approximate $8\times 8$ DCT module [14], that has 288 adders in
a 6-level adder tree. We evaluated the noise power in dB
($10\times\log_{10}MSE$), assuming an error model for LOA based on uniformly
distributed inputs [5] and compared it with values obtained using MC
simulations, with $10^{5}$ uniformly distributed random inputs ($NP_{1}$ and
$NP_{sim}$ respectively in Table I). The numbers after levels ($L_{1}-L_{6}$)
indicate the number of approximate bits. As seen from Table I, the analytical
noise power differs from the simulated value by as much as 11 dB, though the
model worked well for an individual adder.

While the assumption of uniformly distributed lower order bits may be
justified for the primary inputs, neither the output of LPAAs nor the error is
uniformly distributed. A more accurate method of obtaining the probability
mass function (PMF) of error is proposed in [15]. However, including this
method within an optimization routine would require extensive computations.
Moreover, in most applications, an accurate estimate of the mean error and MSE
is sufficient and we do not need the PMF of the error.

TABLE I: Noise power (in dB) of an $8\times 8$ DCT module that uses Lower part
OR adders. $NP_{1}$: Noise power assuming uniform distribution, $NP_{sim}$:
Noise power using simulations, $e=|NP_{sim}-NP_{1}|$.

$NP_{1}$ | $NP_{sim}$ | $e$ | No. of approximate bits
---|---|---|---
-53.45 | -51.20 | 2.25 |  $L_{1}-L_{4}$: 5; $L_{5}-L_{6}$: 6
-51.76 | -44.67 | 7.09 |  $L_{1}-L_{2}$: 5; $L_{3}-L_{6}$: 6
-51.32 | -43.00 | 8.32 |  $L_{1}-L_{2}$: 5; $L_{3}-L_{5}$: 6; $L_{6}$: 7
-50.55 | -41.21 | 9.34 |  $L_{1}-L_{2}$: 5; $L_{3}-L_{4}$: 6; $L_{5}-L_{6}$: 7
-47.59 | -36.50 | 11.09 |  $L_{1}$: 5; $L_{2}-L_{3}$: 6; $L_{4}-L_{5}$: 7; $L_{6}$: 8

In this paper, we show that mean and MSE can be computed accurately without
computing the PMF of the error if the number of approximate bits and the
static probabilities of the inputs is taken into account. To this end, we
derive parameterised error models that can be used within any optimization
framework in order to optimize the number of approximate bits. We also show
that in order to accurately compute the MSE, the optimization framework needs
to take into account not just the functionality of the adder, but also its
position in the circuit, functionality of its parents and the number of
approximate bits in the parent blocks. We demonstrate a significant
improvement of accuracy in the prediction of MSE in applications such as FIR
filter, IIR filter and DCT using a variety of LPAAs.

Each approximate adder requires the static probabilities of its input bits,
which eventually traces back to the primary inputs. Even though the PMF of the
input signal is usually very non-uniform, the PMF of the lower order bits is
often close to uniform, so that simple expressions for the error can be used.
Some justification for this assumption is included in [7], but it is heuristic
and there is no condition that can be checked to test for uniformity. In this
paper, we show that the discrete Fourier transform (DFT) of the input signal
has to satisfy certain constraints if the PMF is uniform. This can be used to
easily find the maximum number of bits that can be considered to be uniformly
distributed.

##  II Parameterized error models for low power approximate adders

As mentioned, the simple error models fail to give accurate estimates of the
mean error and MSE since they typically assume that the inputs or the error is
uniformly distributed. The more complex models that compute the PMF of the
error are more accurate, but unsuitable for use in an optimizer. We propose
parameterized error models with the input static probability as parameters,
that can take into account the input distribution without evaluating the full
PMF. The models are more complex than the expressions in [12], but they are
analytical expressions, suitable for use in an optimizer.

Assume that an $(N,k)$-bit LPAA with inputs $a$ and $b$ has $k$ approximate
bits and $N-k$ accurate bits, with $\hat{s}_{i}=f(a_{i},b_{i},\hat{c}_{i-1})$
and $\hat{c}_{i}=g(a_{i},b_{i},\hat{c}_{i-1})$ denoting the $i^{th}$ bit of
the approximate sum and carry. Let $P_{x_{i}}$ denote the static probability
of a signal $x_{i}$. As with all models in the literature, we continue to
assume that $a$ and $b$ are independent. This assumption is an approximation
when the circuit has reconvergent fanouts. However, it is a reasonable
approximation in many cases as correlations are diluted as the logic depth
increases, as argued in [15].

The error in the output is due to the approximate lower part sum and the
approximate carry to accurate adder. This can be written as

$\displaystyle
e=\sum\limits_{i=0}^{k-1}(a_{i}+b_{i})2^{i}-\sum\limits_{i=0}^{k-1}\hat{s}_{i}2^{i}-\hat{c}_{k-1}2^{k}.$
(1)

The mean error of the approximate adder is given by

$\displaystyle
E\\{e\\}=\sum\limits_{i=0}^{k-1}(P_{a_{i}}+P_{b_{i}})2^{i}-\sum\limits_{i=0}^{k-1}P_{\hat{s}_{i}}2^{i}-P_{\hat{c}_{k-1}}2^{k}.$
(2)

To derive $P_{\hat{s}_{i}}$, we need $P_{\hat{c}_{i-1}}$. One way to do this
is to start with the LSB, for which $P_{\hat{c}_{-1}}=0$ and find the static
probability of all the other carry bits based on the truth table.
Alternatively, we can assume that the static probability of the carry is
independent of the position i.e., $P_{\hat{c}_{i}}=P_{\hat{c}_{i-1}}$. This is
a reasonable assumption as the input static probability is not expected to
vary with location. For example, for AMA-1 if we assume
$P_{\hat{c}_{i}}=P_{\hat{c}_{i-1}}$, we get 0.67 as the static probability,
whereas if we use the first method, we get [0.5,0.625,0.656,0.67,0.67] for the
five LSBs as in [4]. Using 0.67 as the static probability of the carry,
$P_{\hat{s}_{i}}=0.33$ which matches well with actual values for $k>3$. Note
that we have not assumed anything about the PMF of the inputs. If it happens
to be uniform, we can use $P_{a_{i}}=0.5$. However, if the inputs come from
another approximate adder, we just need to use the right value of the static
probability.

The MSE $(E\\{e^{2}\\})$ can be derived in a similar fashion. The expression
for MSE computation requires evaluation of correlated terms such as
$E\\{\hat{s}_{i}\hat{s}_{j}\\}$, $E\\{a_{i}\hat{s}_{j}\\}$,
$E\\{b_{i}\hat{s}_{j}\\}$, $E\\{a_{i}\hat{c}_{k-1}\\}$,
$E\\{b_{i}\hat{c}_{k-1}\\}$ and $E\\{\hat{s}_{i}\hat{c}_{k-1}\\}$ in terms of
the static probabilities. Each of these can be derived using truth table of
the approximate adder. The correlation between bits
$E\\{\hat{s}_{i}\hat{s}_{j}\\}$ can be neglected in most adders i.e., the
individual bits of the sum are independent. However, it is significant in a
few adders such as AMA2 ($\hat{s}_{i}=\overline{\hat{c}_{i}}$), where there is
a close relationship between the sum and carry. An exception to this method
for deriving error models is ETA-I [6], where the lower part sum cannot be
written as a truth table. Its error metrics are derived in our earlier work
[16].

##  III Multi-level systems

In DSP systems, the inputs to the adder are either the primary inputs or they
are output of another approximate adder, subtractor or (in our case, accurate)
multiplier. In each case, the static probabilites of the input needs to be
evaluated correctly. We consider each case in the following subsections.

###  III-A Static Probabilities: Primary inputs

Typical PMF of any primary input such as an image is not uniform. However,
since error expression of a LPAA involves the static probability of the lower
$k$ bits, we would only need to check if the PMF of these $k$ bits is uniform.
A brute-force technique to check would be to compute the PMF of the lower
order bits for each value of $k$ and check for uniformity of the PMF. In this
section, we show that it is possible to check if the distribution is uniform
by computing the DFT of the input. A single DFT is sufficient to find the
values of $k$ for which this assumption is reasonable.

Let $\mathcal{F}_{A}$ be the $2^{N}$-point DFT of the PMF of $N$-bit signal
$A$ and $\mathcal{F}_{A_{L}}$ be the $2^{k}$-point DFT of the PMF of $A_{L}$
($k$ LSBs of $A$). We have,

$\displaystyle\mathcal{F}_{A_{L}}[m]=$
$\displaystyle\sum\limits_{n=0}^{2^{k}-1}P(A_{L}=n)e^{-jmn2\pi/2^{k}}$ (3)
$\displaystyle=$
$\displaystyle\sum\limits_{n=0}^{2^{k}-1}\sum\limits_{l=0}^{2^{N-k}-1}P(A=l2^{k}+n)e^{-jmn2\pi/2^{k}}$
$\displaystyle=$
$\displaystyle\sum\limits_{l=0}^{2^{N-k}-1}\sum\limits_{n^{\prime}=l2^{k}}^{l2^{k}+2^{k}-1}P(A=n^{\prime})e^{-jmn^{\prime}2\pi/2^{k}}e^{jml2^{k}2\pi/2^{k}}$
$\displaystyle=$
$\displaystyle\sum\limits_{n^{\prime}=0}^{2^{N}-1}P(A=n^{\prime})e^{-jmn^{\prime}2\pi/2^{k}}$
$\displaystyle=$ $\displaystyle\mathcal{F}_{A}[m\cdot 2^{N-k}],0\leq m<2^{k}.$
(4)

If $A_{L}$ is uniform, $P(A_{L}=n)=\dfrac{1}{2^{k}},0\leq n<2^{k}$. Hence from
(3), if $A_{L}$ is uniform, we have

$\displaystyle\mathcal{F}_{A_{L}}[m]=$
$\displaystyle\sum\limits_{n=0}^{2^{k}-1}\dfrac{1}{2^{k}}e^{-jmn2\pi/2^{k}}=\begin{cases}1,&\text{if
}m=0\\\ 0,&\text{if }0<m<2^{k}.\end{cases}$ (5)

Since DFT is unique, the converse is also true. Therefore using (4), we have
the following condition to be satisfied for $A_{L}$ to be uniformly
distributed.

$\displaystyle\mathcal{F}_{A}[m\cdot 2^{N-k}]$
$\displaystyle=\begin{cases}1,&\text{if }m=0\\\ 0,&\text{if
}0<m<2^{k}.\end{cases}$ (6)

In [17], they have similar condition for continuous signals that are
quantized, although the derivation is a little more involved.

$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$11$$12$$13$$14$$15$$0$$0.05$$0.1$$0.15$$m$$\mathcal{F}_{A}{[}m\cdot
2^{N-k}{]}$$k=1$$k=2$$k=3$$k=4$$k=5$ Figure 1: Illustration of condition for
$k$ lower-order bits of Cameraman image to be uniform.

$0$$50$$100$$150$$200$$250$ (a) $N=8$

$0$$1$ (b) $k=1$

$0$$1$$2$$3$ (c) $k=2$

$0$$2$$4$$6$ (d) $k=3$

$0$$5$$10$$15$ (e) $k=4$

$0$$5$$10$$15$$20$$25$$30$ (f) $k=5$

Figure 2: (a) PMF of Cameraman image; (b)-(f) PMF of the lower $k$ bits of the
image.

To illustrate this condition (6), let us consider the Cameraman image with
$N=8$. For the image pixel distribution, $\mathcal{F}_{A}[m\cdot 2^{N-k}]$ for
different values of $k$, $m$ varying from $0$ to $2^{k}-1$, is plotted in Fig.
1. It is seen that for lower values of $k$, the value of the transform is very
close to zero for $0<m<2^{k}$. As $k$ increases, the value of transform also
increases and for $k=5$, the values are high. This is confirmed from the
actual PMF of the lower order bits of the image shown in Fig. 2. From the
figure, it is seen that distribution can be considered uniform even if half
the bits are approximated. This turns out to be true for all the standard
images we have looked at (for example, for the Lena image, $k=6$, for Rice,
$k=5$). Hence, we assume that the PMF of the lower order bits of primary
inputs to the approximate adder is uniform i.e., the static probability is
0.5.

###  III-B Static Probabilities: Adders in the higher levels

If the inputs to the adder are the output of another adder as in Fig. 3(a),
the mean and MSE are derived using $P_{\hat{s}_{i}}$ and $P_{\hat{c}_{k-1}}$
as discussed previously.

$\boldsymbol{+}$$\boldsymbol{+}$$\boldsymbol{+}$$N$$N$$N$$N$$m_{1},$$k_{1}$$m_{2},$$k_{2}$$m_{3},$$k_{3}$Adder1Adder2Adder3
(a)

$\boldsymbol{+}$$\boldsymbol{+}$$\boldsymbol{+}$$\boldsymbol{*}$$N$$N$$N$$N$$m_{1},$$k_{1}$$(P_{a_{i}})$$m_{2},$$(P_{i})$$(P_{b_{i}})$$k_{2}$$m_{3},$$k_{3}$cAdder1Adder2Adder3
(b)

Figure 3: Adder tree in a circuit with $N$ bits at primary inputs.

The other possibility is that input is the output of a multiplier. In this
work, we are only optimizing adders and all multipliers are accurate, with the
output truncated to the standard precision used in the circuit. Also, we only
consider linear systems, so that one of the inputs to the multiplier is a
constant coefficient. In Fig. 3(b), consider Adder3 which has an input from
the output of the multiplier. Depending on the value of the constant
coefficient $c$, the probability of the LSBs at the output of the multiplier
($P_{b_{i}}$) will vary. Let $P_{i}$ denote the probability that the $i^{th}$
bit of the $k_{2}$ LSBs of the multiplicand (output of Adder2) is 1. Consider
the following cases.

  1. 1.

When $c=2^{l}$ and $l\geq 0$, the product is the logical left shift of the
multiplicand. So $P_{b_{i}}=0$ for the first $l$ LSBs and $P_{b_{i}}=P_{i-l}$
for the next $k_{2}-l$ LSBs.

  2. 2.

When $c=-2^{l}$ and $l\geq 0$, $P_{b_{i}}=0$ for the first $l$ LSBs and
$P_{b_{i}}=1-P_{i-l}$ for the next $k_{2}-l$ LSBs is a good approximation,
accounting for flipping involved in the two’s complement representation for
negative numbers.

  3. 3.

When $c=2^{l}$ and $l<0$, the product is the right shift of the multiplicand.
So $P_{b_{i}}=P_{i+|l|}$ for $k_{2}-l$ LSBs.

  4. 4.

When $c=-2^{l}$ and $l<0$, $P_{b_{i}}=1-P_{i+|l|}$ for $k_{2}-l$ LSBs.

  5. 5.

If $c$ is an arbitrary constant that is not a power of 2, MC simulations
indicate that the average static probability of the output bits is $0.5\pm
0.03$. This is also intuitively correct, since the product is the sum of
several partial products and the LSBs are truncated to the precision
maintained in the system. Hence, the overall PMF is likely to be symmetric
(moving towards Gaussian), which means that the static probability is 0.5. So
we assume that $P_{b_{i}}=0.5$ for $k_{2}$ LSBs.

###  III-C Truncation and Median adder (MA) in higher levels

Both Truncation and MA [7] have their lower part sum bits fixed to constant
all 0’s and 1’s respectively. In these adders, since the lower part sum is
known, the lower part sum of the adders in higher levels can be fixed more
accurately so that the accuracy of the approximate circuit is improved. In
case of Truncation adder, the approximate sum is obviously zero. In Fig. 3(a)
with MA, Adder1 and Adder2 will have their lower part sum as $2^{k_{1}}-1$ and
$2^{k_{2}}-1$, respectively. For Adder3, we lower the MSE of the circuit by
setting the sum to $2^{k_{3}+1}-1$ instead of $2^{k_{3}}-1$ for the following
cases:

  1. 1.

When $k_{3}\leq k_{1},k_{2}$, the lower part sum is known exactly and is equal
to $2^{k_{3}+1}-2$.

  2. 2.

When $k_{1}\geq k_{3}>k_{2}$, the mean of the sum is $(3\times
2^{k_{3}}+2^{k_{2}}-4)/2$, which is closer to $2^{k_{3}+1}-1$ than
$2^{k_{3}}-1$.

Using this setting, we obtain up to 6 dB improvement in MSE for the adder tree
in Fig. 3(a). Beyond the first level, the input static probability for Median
adders is set to 1.

##  IV Optimization framework

Classical wordlength optimization uses a simple model for the quantization
error and the same model is used for all nodes in the circuit. In the
literature, a similar framework with a single expression for error for all the
adders has been used to find the optimal number of approximate bits. However,
the discussion in the previous section shows that this is inadequate to get
accurate numbers. Hence, we made several modifications to the framework, which
are detailed below.

The input to the optimizer is the circuit implemented using adders,
multipliers and registers and the corresponding signal flow graph. The primary
inputs to the system are normalized to $1.N$ fixed point numbers. For each
functional unit in the system, we use the required number of integer bits
while maintaining the number of fractional bits as $N$. The goal of the
optimizer is to maximize the number of approximate fractional bits of the
adders in the circuit for a given MSE at the output. We use the three step
procedure discussed in [18] and adapt it for approximate computing. It uses
Minimum Width algorithm, Mildest Greedy Ascent algorithm and Tabu search
algorithm. The main steps in our optimization framework are as follows:

  * •

We have a pre-processing step in which adder gets the static probability of
the inputs based on its parent nodes. If the parent is a register, it gets the
static probability of the inputs to the register. Also, the transfer function
from each adder node to the output is computed.

  * •

Next, we run a minimum width algorithm (MWA) that gives the maximum number of
approximate fractional bits at the output of each adder when all the other
adders are accurate and the required MSE constraint is satisfied.

  * •

Starting with the number of approximate bits from the MWA, a greedy descent
procedure is used to decrease the number of approximate bits in the adder that
causes the maximum improvement in MSE. An important difference from the
quantization noise optimization in [18] is that the approximation noise can
worsen even if the number of approximate bits is decreased. For circuits with
multiple outputs, we find the fan-in cone of each output in sequence. The
adders in the fan-in cone are optimized while keeping the adders in the fan-in
cone that increase the MSE of previously targeted outputs untouched.

  * •

Finally we run a tabu search algorithm targeting signals with minimum number
of approximate bits (instead of the most sensitive signal) and keep increasing
the number of approximate bits as long as MSE constraint is met. We found that
this heuristic provides better optimization of approximate bits.

In each of these algorithms, the overall MSE at the output is computed using
the transfer functions from the adder nodes and the parameterized error model
for the adder.

##  V Experimental Results

In this section, we first validate our error model. Each assumption is tested
against simulation and verified for correctness. We then obtain the optimum
number of approximate bits for a given MSE using our optimization framework.
This is done for FIR and IIR filters and an $8\times 8$ DCT module. These
results are used to show that the optimizer requires the parameterized error
models with the correct values of the static probability for accurate
prediction of the MSE.

In order to validate our error model, we used the simple adder tree depicted
in Fig. 3(a).

$2$$3$$4$$5$$6$$7$$8$$0$$1$$2$$3$Number of approximate bits ($k$)Error in
noise power
(dB)$(1)P_{a,i}=P_{b,i}=0.5$$(2)P_{a_{i}a_{j}}=P_{a,i}P_{a,j}$$(3)P_{c,i}=P_{c,i-1}$(4)
Param. model (a) AMA1 in the second level

$2$$3$$4$$5$$6$$7$$8$$0$$2$$4$$6$Number of approximate bits ($k$)Error in
noise power (dB) (b) AMA2 in the second level

Figure 4: Error in noise power computation of Adder3 in Fig.3(a) when the tree
is implemented using (a) AMA-1 and (b) AMA-2 adders.

Some of the LPAAs like AMA-1 and AMA-2 adders involve carry propagation in the
approximate lower part sum. Therefore, in addition to the correct value of the
static probability, evaluation of the MSE also requires correlations between
bits to be taken into account. Fig. 4 shows a comparison of the error in noise
power computation of Adder3 in Fig. 3(a) (1) assuming $P_{a,i}=P_{b,i}=0.5$
(2) assuming that the individual bits of each input are independent (i.e.
$P_{a_{i}a_{j}}=P_{a,i}P_{a_{j}}$) (3) assuming that $P_{c,i}=P_{c,i-1}$ (4)
using the parameterized error model including all the correlations. It can be
seen that, $P_{c,i}=P_{c,i-1}$ is a good approximation for $k>3$. The
correlations between the bits in each input can be ignored in AMA-1, AMA-2
requires all correlations to be taken into account. From the discussion in
Section II, both these results are as expected.

We validated our optimization framework using an 18-tap FIR filter (direct
form I realization), IIR filter (direct form II realization of a 4th order low
pass Butterworth filter) and DCT [14] consisting of 17, 8 and 288 adders
respectively. These are typically the benchmarks that have been used in the
literature. For each of these systems, we obtained the optimum number of
approximate bits for each adder in the system, given an overall MSE
specification. This was done for the Truncation, MA, AMA-5, LOA and ETA-I
adders. These adders have minimal hardware for evaluation of the approximate
sum and are energy-efficient. Using the optimal configuration of approximate
bits for each adder obtained from the optimizer, the circuits are implemented
and simulated with $10^{5}$ uniform random inputs to obtain the actual MSE of
the circuits. The actual MSE is compared with the analytical value obtained
using the parameterized error model.

Figs. 5 and 6 show the comparison of simulated value of noise power and
analytically computed noise power in dB for the FIR filter and IIR filter
respectively using various approximate adders for two cases - (a) assuming
that the input static probabilities are 0.5 for all the adders (b) computing
the input static probabilities for adders as described in Section III-B. It
can be seen that the analytical values obtained by using the parameterized
error model match very well with the actual values, while those obtained using
uniform probabilities deviate significantly (upto 6 dB and 10 dB error in case
of FIR and IIR filters respectively).

$-60$$-50$$-40$$-30$$-20$$-60$$-50$$-40$$-30$$-20$Actual noise power
(dB)Analytical noise power (dB)TruncAMA5LOAETA-IMA (a) Uniform distribution

$-60$$-50$$-40$$-30$$-20$$-60$$-50$$-40$$-30$$-20$Actual noise power
(dB)Analytical noise power (dB)TruncAMA5LOAETA-IMA (b) Parameterized error
model

Figure 5: Comparison of simulated noise power and analytically computed noise
power for FIR filter when input static probability is (a) assumed to be 0.5
(b) computed as in Section III-B for various approximate adders.

$-60$$-55$$-50$$-45$$-40$$-35$$-30$$-60$$-50$$-40$$-30$Actual noise power
(dB)Analytical noise power (dB)TruncAMA5LOAETA-IMA (a) Uniform distribution

$-60$$-55$$-50$$-45$$-40$$-35$$-30$$-60$$-50$$-40$$-30$Actual noise power
(dB)Analytical noise power (dB)TruncAMA5LOAETA-IMA (b) Parameterized error
model

Figure 6: Comparison of simulated noise power and analytically computed noise
power for IIR filter when input static probability is (a) assumed to be 0.5
(b) computed as in Section III-B for various approximate adders.

The DCT module [14] is a multiplierless implementation with adders and
subtractors. Fig. 7(a) and Fig. 7(b) show the comparison of simulated value of
noise power and analytically computed noise power using the parameterized
error model for various approximate adders, for two different images namely
Rice image and Lena image respectively. It can be seen that the analytical
values match reasonably well with the actual values (with a maximum of 2.5 dB
and 3 dB error for Rice and Lena images respectively).

$-40$$-35$$-30$$-25$$-20$$-40$$-35$$-30$$-25$$-20$Simulated noise power
(dB)Analytical noise power (dB)TruncAMA5LOAMA$NP_{anal}=NP_{sim}$ (a) Rice
image input

$-40$$-35$$-30$$-25$$-20$$-40$$-35$$-30$$-25$$-20$Simulated noise power
(dB)Analytical noise power (dB)TruncAMA5LOAMA$NP_{anal}=NP_{sim}$ (b) Lena
image input

Figure 7: Comparison of simulated noise power and analytically computed noise
power for $8\times 8$ DCT using parameterized error model.

##  VI Conclusion

We have proposed parameterized error models for approximate adders using input
static probabilities as parameters and incorporated these error models in an
optimization framework. We have shown that the parameterized error models
provide better noise power prediction than the typical error models that
assume uniform input distribution. The results of FIR and IIR filters and DCT
module show that the use of parameterized error model in the optimization
framework improves the accurate prediction of the overall MSE.

## References

  * [1] M. Shafique, W. Ahmad, R. Hafiz, and J. Henkel, “A low latency generic accuracy configurable adder,” in Proceedings of the 52Nd Annual Design Automation Conference, DAC ’15, (New York, NY, USA), ACM, 2015.
  * [2] H. A. F. Almurib, T. N. Kumar, and F. Lombardi, “Inexact designs for approximate low power addition by cell replacement,” in Design, Automation and Test in Europe (DATE), 2016.
  * [3] Z. Yang, A. Jain, J. Liang, J. Han, and F. Lombardi, “Approximate xor/xnor-based adders for inexact computing,” 2013 13th IEEE International Conference on Nanotechnology (IEEE-NANO 2013), 2013.
  * [4] V. Gupta, D. Mohapatra, A. Raghunathan, and K. Roy, “Low-power digital signal processing using approximate adders,” IEEE Trans. on Comp.-Aided Design of Integrated Circuits and Systems, vol. 32, 1 2013.
  * [5] H. R. Mahdiani, A. Ahmadi, S. M. Fakhraie, and C. Lucas, “Bio-inspired Imprecise computational blocks for efficient VLSI implementation of soft-computing applications,” IEEE Trans. on Circuits and Systems I: Regular Papers, vol. 57, 4 2010.
  * [6] N. Zhu, W. L. Goh, W. Zhang, K. S. Yeo, and Z. H. Kong, “Design of low-power high-speed truncation-error-tolerant adder and its application in digital signal processing,” IEEE TVLSI, vol. 18, no. 8, 2010.
  * [7] D. Celia, V. Vasudevan, and N. Chandrachoodan, “Optimizing power-accuracy trade-off in approximate adders,” in 2018 Design, Automation Test in Europe Conference Exhibition (DATE), March 2018.
  * [8] L. B. Soares, S. Bampi, and E. Costa, “Approximate adder synthesis for area- and energy-efficient fir filters in cmos vlsi,” in 2015 IEEE 13th International New Circuits and Systems Conference, June 2015.
  * [9] F. S. Snigdha, D. Sengupta, J. Hu, and S. S. Sapatnekar, “Optimal design of jpeg hardware under the approximate computing paradigm,” in 2016 53nd ACM/EDAC/IEEE DAC, June 2016.
  * [10] Z. Vasicek, V. Mrazek, and L. S. Brno, “Towards low power approximate dct architecture for hevc standard,” in Design, Automation Test in Europe Conference Exhibition (DATE), 2017, March 2017.
  * [11] D. Sengupta, F. S. Snigdha, Jiang Hu, and S. S. Sapatnekar, “Saber: Selection of approximate bits for the design of error tolerant circuits,” in 2017 54th ACM/EDAC/IEEE DAC, June 2017.
  * [12] M. Pashaeifar, M. Kamal, A. Afzali-Kusha, and M. Pedram, “A theoretical framework for quality estimation and optimization of dsp applications using low-power approximate adders,” IEEE Trans. on Circuits and Systems I: Regular Papers, vol. 66, Jan 2019.
  * [13] Z. Yang, J. Han, and F. Lombardi, “Transmission gate-based approximate adders for inexact computing,” in Proceedings of the 2015 IEEE/ACM NANOARCH´15, July 2015.
  * [14] S. Bouguezel, M. O. Ahmad, and M. N. S. Swamy, “A fast $8\times 8$ transform for image compression,” in 2009 International Conference on Microelectronics - ICM, Dec 2009.
  * [15] D. Sengupta, F. S. Snigdha, J. Hu, and S. S. Sapatnekar, “An analytical approach for error pmf characterization in approximate circuits,” IEEE Trans. on CAD of Integrated Circuits and Systems, vol. 38, Jan 2019\.
  * [16] D. Celia, V. Vasudevan, and N. Chandrachoodan, “Probabilistic error modeling for two-part segmented approximate adders,” in 2018 IEEE International Symposium on Circuits and Systems (ISCAS), May 2018.
  * [17] A. Sripad and D. Snyder, “A necessary and sufficient condition for quantization errors to be uniform and white,” IEEE Trans. on Acoustics, Speech, and Signal Processing, vol. 25, October 1977.
  * [18] D. Menard, N. Herve, O. Sentieys, and H.-N. Nguyen, “High-level synthesis under fixed-point accuracy constraint,” Journal of Electrical and Computer Engineering, vol. 2012, Jan. 2012.


#  GomalizingFlow.jl: A Julia package for Flow-based sampling algorithm for
lattice field theory

Akio Tomiya  Faculty of Technology and Science, International Professional
University of Technology, 3-3-1, Umeda, Kita-ku, Osaka, 530-0001, Osaka, Japan
Satoshi Terasaki  AtelierArith, 980-0004, Miyagi, Japan

###### Abstract

GomalizingFlow.jl: is a package to generate configurations for quantum field
theory on the lattice using the flow-based sampling algorithm in Julia
programming language. This software serves two main purposes: to accelerate
research of lattice QCD with machine learning with easy prototyping, and to
provide an independent implementation to an existing public Jupyter notebook
in Python/PyTorch. GomalizingFlow.jl implements, the flow based sampling
algorithm, namely, RealNVP and Metropolis-Hastings test for two dimension and
three dimensional scalar field, which can be switched by a parameter file. HMC
for that theory also implemented for comparison. This package has a Docker
image, which reduces the effort for environment construction. This code works
on both the CPU and the NVIDIA GPU.

###### keywords:

Lattice QCD , Particle physics , Machine learning , Normalizing flow , Julia

††journal: SoftwareX

## Required Metadata

## Current code version

Nr. |  Code metadata description |  Please fill in this column
---|---|---
C1 |  Current code version |  v1.0.0
C2 |  Permanent link to code/repository used for this code version |  https://github.com/AtelierArith/GomalizingFlow.jl
C3 |  Code Ocean compute capsule |  To be uploaded. Docker image available
C4 |  Legal Code License |  MIT
C5 |  Code versioning system used |  Git
C6 |  Software code languages, tools, and services used |  Julia
C7 |  Compilation requirements, operating environments & dependencies |  Julia 1.8.2, Flux.jl v0.13.5, (Docker image is provided)
C8 |  If available Link to developer documentation/manual |  https://github.com/AtelierArith/GomalizingFlow.jl/blob/main/README.md
C9 |  Support email for questions |<EMAIL_ADDRESS>
Table 1: Code metadata (mandatory)

##  1 Motivation and significance

A lattice gauge theory with Markov chain Monte-Carlo calculations provides
quantitative information on quantum field theory and particle physics [1, 2,
3, 4]. Regularized quantum field theory is defined on discretized spacetime
(lattice), and the discretization makes degrees of freedom finite. This
enables us to perform numerical calculations with Markov chain Monte-Carlo on
supercomputers. Quantum field theory on a lattice has been applied mainly for
QCD (Quantum Chromo-dynamics), which is a fundamental theory for the inside of
nuclei [5]. Lattice calculations of QCD (Lattice QCD) have succeeded in
reproducing the hadron spectrum [6], and enable us to access non-perturbative
and quantitative information which is difficult to obtain with other methods
[7, 8, 9].

The ultraviolet cutoff is necessary to define a finite theory. However, the
ultraviolet cutoff, namely finite lattice spacing, bring unwanted artifacts
into the theory. Finer lattice spacing gives smaller lattice artifacts, but it
causes long autocorrelation among the Monte Carlo samples, which makes Monte
Carlo calculations inefficient because the random samples produced are
correlated, which is called critical slowing down [10]. To make Monte-Calo
samples in good quality in a short time, a new algorithm with shorter
autocorrelation is demanded because the autocorrelation time depends on an
algorithm.

The trivializing map has been investigated to resolve the critical slowing
down in lattice QCD calculations [11]. Trivializing maps is a change of
variables in the (path-)integral between a physical theory to a trivial
theory, which is easy to calculate. A concrete example of the trivializing map
called the Nicoli map has been realized in a supersymmetric field theory in
the first [12]. For non-super symmetric theory, M. Luscher has shown that the
gradient flow with Jacobian is an approximate trivializing map, but his
original idea has not worked well because of the numerical cost.

Machine learning algorithms to generate configurations for lattice field
theory are widely investigated [13]. The first trial has been performed with
the Boltzmann machine [14], and others such as GAN [15, 16] have been tried.
The gauge covariant neural network [17] is applied to perform simulations for
four-dimensional QCD. A neural network parameterized leapfrog integrator
enhances topology change in simulations [18]. The normalizing flow [19, 20,
21, 22, 23, 24, 25, 26, 27] has recently been extensively investigated, and we
discuss this further later. Generative models are machine learning
architecture for image generation, and which has been investigated to generate
good quality images. Images are two-dimensional data with certain structures,
and it seems that it can be applied to generate field configurations.

The flow based sampling algorithm, based on Real-NVP (Non-volume preserving
map), which is a generative model, and is an exact algorithm of Monte-Carlo
configuration generation for quantum field theories associated with the
Metropolis–Hastings test. It has the potential to improve critical slowing
[19, 20], which has been extended to system with SU(N) links [21] and with
fermions [22, 23, 24]. This can be regarded as a realization of a trivializing
map a la M. Luscher using a neural network.

One of the purposes of quantum field theory is to calculate the expectation
value in the following form,

$\displaystyle\left\langle
O\right\rangle\equiv\int\mathcal{D}\phi\;P^{\text{(QFT)}}_{\\{L,m^{2},\lambda\\}}[\phi]\;O[\phi],$
(1)

where

$\displaystyle\mathcal{D}\phi\;=\prod_{n=n_{o}}^{L^{d}}d\phi_{n}$ (2)

and $n_{o}=\underbrace{(1,1,\cdots)}_{d}$, the origin111 In Julia programming
language, index is started from $1$ mostly. We follow this notation. , $d$ is
the dimensionality, and $L$ is a size of the system. And
$P^{\text{(QFT)}}_{\\{L,m^{2},\lambda\\}}[\phi]$ is a probability density of
the system,

$\displaystyle P^{\text{(QFT)}}_{\\{L,m^{2},\lambda\\}}[\phi]$
$\displaystyle=\frac{1}{Z}\text{e}^{-S^{\text{(QFT)}}[\phi]},$ (3)

and $Z$ normalizes the expectation values as $\left\langle 1\right\rangle=1$.
The action is,

$\displaystyle
S^{\text{(QFT)}}[\phi]=-\sum_{n=n_{o}}^{L^{d}}\phi(n_{o})\partial^{2}\phi(n_{o})+\sum_{n=n_{o}}^{L^{d}}V[\phi](n).$
(4)

where
$\partial^{2}\phi(n)=\sum_{\mu}[\phi(n+\hat{\mu})+\phi(n-\hat{\mu})-2\phi(n)]$
is a discretized Laplacian and,

$\displaystyle V[\phi](n)={{m}^{2}}\phi^{2}(n)+\lambda\phi^{4}(n),$ (5)

is a potential term.

To calculate the integral Eq. (1), we use Markov Chain-Monte Carlo. We
generate a sequence of configurations,

$\displaystyle\phi^{(1)}\to\phi^{(2)}\to\phi^{(3)}\to\phi^{(4)}\to\cdots\to\phi^{(N_{\rm
conf})},$ (6)

where $\phi^{(k)}$ is a configuration sampled from Eq. (3) and ${N_{\rm
conf}}$ is the number of samples. Typically, HMC (Hybrid / Hamiltonian Monte
Carlo) has been used [28, 29]. We can calculate the expectation value through

$\displaystyle\left\langle O\right\rangle=\lim_{N_{\rm
conf}\to\infty}\frac{1}{N_{\rm conf}}\sum_{k=1}^{N_{\rm conf}}O[\phi^{(k)}].$
(7)

In practice, we cannot take $N_{\rm conf}$ infinity, the expectation value is
suffered from the statistical error,

$\displaystyle\left\langle O\right\rangle=\frac{1}{N_{\rm
conf}}\sum_{k=1}^{N_{\rm
conf}}O[\phi^{(k)}]+\mathcal{O}\left(\frac{1}{\sqrt{N_{\rm indep}}}\right)$
(8)

where $N_{\rm indep}$ is the number of independent samples defined by

$\displaystyle N_{\rm indep}=\frac{N_{\rm conf}}{2\tau_{ac}}$ (9)

and $\tau_{ac}$ is the autocorrelation time [30]. In some parameter regime,
typically critical regime, $\tau_{ac}$ grows and the number of independent
configurations $N_{\rm indep}$ becomes small, which is called the critical
slowing down. The autocorrelation time $\tau_{ac}$ depends on an algorithm,
and here we implement the flow based sampling algorithm to reduce the
autocorrelation time.

The flow-based sampling algorithm utilizes three probability distributions.
First is, the target distribution. We want to draw a sample from the target as

$\displaystyle\phi\sim P^{\text{(QFT)}}_{\\{L,m^{2},\lambda\\}}[\phi].$ (10)

where “$\sim$” represents “sampling” from a probability distribution.

Second one is, a trivial distribution (a prior distribution in terms of
machine learning),

$\displaystyle\varphi\sim P^{\text{(pri)}}_{\vartheta}[\varphi],$ (11)

which has a set of parameters $\vartheta$, this is a fixed parameter of the
calculations (hyperparameter). This probability distribution should be point-
wise independent distribution and this can be regarded as a physical
distribution at infinite temperature for spin system, or, $\beta=0$
distribution for lattice gauge systems. In terms of quantum field theory, one
can regard a probability distribution for a field without a kinetic term.

We apply a parameterized transformation (like a cooling process),

$\displaystyle\phi=\mathcal{F}^{-1}_{\theta}[\varphi],$ (12)

where $\theta$ is a set of parameters on the map and such that,

$\displaystyle P^{\text{(QFT)}}_{\\{L,m^{2},\lambda\\}}[\phi]\approx
P^{\text{(ML)}}_{\Theta}[\phi]=P^{\text{(pri)}}_{\vartheta}[\varphi]\Bigg{|}\frac{\partial\mathcal{F}^{-1}_{\theta}[\varphi]}{\partial\varphi}\Bigg{|},$
(13)

and $\Theta=\theta\cup\vartheta$ and

$\displaystyle\Bigg{|}\frac{\partial\mathcal{F}^{-1}_{\theta}[\varphi]}{\partial\varphi}\Bigg{|},$
(14)

is the Jacobian for the transformation. The quality of this map is measured by
(shifted and reversed) Kullback–Leibler divergence,

$\displaystyle\int\mathcal{D}\phi\;P^{\text{(ML)}}_{\Theta}[\phi]\ln
P^{\text{(ML)}}_{\Theta}[\phi]/P^{\text{(QFT)}}_{\\{L,m^{2},\lambda\\}}[\phi]-\ln
Z,$ (15)

and this works as a variational energy of approximation. By tuning parameters
$\theta$ to mimize Kullback–Leibler divergence, we can realize approximated
(un) trivializing map.

Fig. 1 is a schematic overview for the flow-based sampling algorithm. First,
it samples from a configuration with point-wise independent distribution (in
terms of gauge theory, it is sampling from $\beta=0$). Next, we apply the flow
transformation $\mathcal{F}^{-1}_{\theta}$, un-trivializing map, written by a
neural network. This makes correlations between sites. Output from the network
can be regarded as a sample from $\beta>0$ approximately. Finally, with the
Metropolis-Hastings test with the Jacobian, we obtain an ensemble which obeys
the correct distribution ${\text{e}}^{-S}/Z$. Training also follows the
process described above.

Figure 1: Schematic picture for the flow based sampling algorithm. Here we use
$\beta$ to indicate triviality of the system, which is a coefficient of the
kinetic term in our current application (please see the main text in details).

Here, we describe how we can realize the map $\mathcal{F}_{\theta}$
(trivializing map) and $\mathcal{F}^{-1}_{\theta}$ (untrivializing map) in
terms of neural networks. The structure of the affine coupling layer is
essential for a tractable Jacobian. We employ checker-board decomposition to
get tractable Jacobian as in the original work [19]. We alternately apply
small and reversible trivializing maps, which are realized using a neural
network, on even and odd sites. By composing these maps, we approximately
realize a map between the target distribution to the trivial distribution and
one for the opposite direction. The affine coupling layer is realized as
follows.

We consider the affine coupling layer on odd sites first. We denote $\phi$ and
$\varphi$ as nontrivial and trivial side variables, respectively. The affine
coupling layer (scaling and shift transformation for field variable) is
defined as,

$\displaystyle\begin{cases}\varphi_{e}&=\text{e}^{s[\phi_{o}]}\phi_{e}+t[\phi_{o}],\\\
\varphi_{o}&=\phi_{o}\end{cases}$ (16)

where $s$ and $t$ are output of a neural network,

$\displaystyle T^{o}_{\theta}:\phi_{o}\mapsto\begin{pmatrix}s\\\
t\end{pmatrix}=\begin{pmatrix}s\\\
t\end{pmatrix}[\phi_{o}]=\begin{pmatrix}s\\\
t\end{pmatrix}(\phi_{1},\phi_{3},\cdots)$ (17)

and the arguments are field values on odd sites for a configuration. We remark
that this transformation is bijective. Transformations for even sites
$T^{e}_{\theta}$ are defined in a similar manner. In this work, this map is
realized using a two- or three-dimensional convolutional neural network. We
can realize an approximate trivializing map by applying this map with
different weights repeatedly,

$\displaystyle\mathcal{F}_{\theta}[\phi]=T^{o}_{\theta_{N_{\rm tri}}}\circ
T^{e}_{\theta_{N_{\rm tri}-1}}\circ\cdots\circ T^{e}_{\theta_{2}}\circ
T^{o}_{\theta_{1}}[\phi].$ (18)

where $\theta=\theta_{1}\cup\theta_{2}\cup\cdots\cup{\theta_{N_{\rm tri}}}$
and $\theta_{i}$ is a set of weights in each neural network, and ${N_{\rm
tri}}$ is the number of affine coupling layers. This map
$\mathcal{F}_{\theta}[\varphi]$ is also bijective, since $T$ is bijective.
Thus, we can obtain an untrivializing map $\mathcal{F}^{-1}_{\theta}$. Due to
the even odd structure of the transformation, the Jacobian calculation is
$O(N)$ [19]. We remark that, (un-)tirivializing procedure is analogous to
integration steps in the gradient flow [17].

This package is implemented by Julia, which is a new scientific programming
language [31]. Typically Julia is as fast as C and Fortran [32], and is fit to
calculations in the lattice QCD community. The flow based sampling algorithm
has been implemented in Python/PyTorch [25], and this package provides another
implementation.

##  2 Software description

###  2.1 Software Architecture

This package works as follows (Fig. 2). First, a parameter configuration file
(e.g. “cfgs/example2d.toml”) is loaded. The parameters in the file are
classified into four groups. $DeviceParams$ ($dp$) specifies what kind of
acceleration device is used (-1 for no accelerator). $TrainingParams$ ($tp$)
controls the training procedure (, e.g. epoch, batch size, optimizer).
$PhysicalParams$ ($pp$) should contain parameters of the target system (mass,
coupling, and system size). $ModelParams$ ($mp$) is a set of parameters for a
neural network in the flow model. These parameters are loaded through
$load\\_hyperparams$ function.

The code constructs objects according to the loaded hyperparameters and
physical parameters in the memory. The objects is for example, action of the
target theory, neural network model, an optimizer, and setting of training
like, the number of epochs, batch size, and scheduling of learning rate
($lr$).

The simulation is performed with parameters with the action defined in
$src/actions.jl$. Jacobian can be calculated from the output of the neural
network.

Figure 2: An illustration of code architecture

###  2.2 Software Functionalities

In the current version, the flow based sampling algorithm for 2 and 3
dimensional scalar field theory can be used in a Docker environment, which is
explained in the following section. HMC for these theories are implemented to
compare results. HMC uses the same struct for parameters.

Jupyter lab environment for analysis, e.g. the condensate, zero-momentum two
point functions, calculation of the autocorrelation time, are available in
$playground/notebook/julia/analysis\\_tool.md$. This file is a markdown file
but will be automatically converted to an ipynb file via jupytext.

##  3 Illustrative Examples

The use of Docker images is the easiest way to use this package222 To use
Docker with NVIDIA GPU, we recommended to follow NVIDIA provided procedure in
https://docs.nvidia.com/datacenter/cloud-native/container-toolkit/install-
guide.html#docker. . In the repository, one can execute the training, with
commands in Listing 1. A file $cfgs/example2d.toml$ is a parameter file, which
is explained below.

⬇

1$ make

2$ docker-compose run \--rm julia julia begin_training.jl cfgs/example2d.toml

Listing 1: Command to start training

Then the training for the two dimensional scalar field will be started. The
results are saved in the $results$ directory.

We provide example files for two-dimensional scalar field theory and three-
dimensional scalar field theory for broken and symmetric phases. Hyper-
parameters for training and parameters for the physical system can be written
as Listing 2. This performs training with 200 epochs and a prior distribution
a Gaussian $N(0,1)$. A neural network has 16 layers with hidden size, $8\times
8$. The training parameter can be understood by comments that started on line
14 of the list.

⬇

1[config]

2version = "0.1.0"

3

4[device]

5# setup DeviceParams for example

6# device_id = -1 # <\--- train with CPU

7# device_id = 0 # <\--- train with GPU its Device ID is 0

8# device_id = 1 # <\--- train with GPU its Device ID is 1

9device_id = -1

10

11

12# setup TrainingParams

13#

14# for _ in 1:epochs

15# for _ in 1:iterations

16# # extract batchsize data

17# z = rand(prior, lattice_shape..., batchsize)

18# gs = Flux.gradient(ps) do

19# do something

20# end

21# Flux.Optimise.update!(opt(base_lr), ps, gs)

22# end

23# end

24[training]

25seed = 12345

26batchsize = 64

27epochs = 200

28iterations = 100

29base_lr = 0.001

30opt = "Adam"

31prior = "Normal{Float32}(0.f0, 1.f0)"

32lr_scheduler = ""

33pretrained = ""

34

35# setup PhysicalParams

36#

37# lattice_shape = (L, L) if Nd = 2

38# lattice_shape = (L, L, L) if Nd = 3

39# action = ScalarPhi4Action(M2, lam)

40[physical]

41L = 8

42Nd = 2

43M2 = -4.0 # m^2

44lam = 8.0 # λ

45

46

47# setup ModelParams

48#

49[model]

50seed = 2021

51n_layers = 16

52hidden_sizes = [8, 8]

53kernel_size = 3

54inC = 1

55outC = 2

56use_final_tanh = true

57use_bn = false

Listing 2: example2d.toml

##  4 Impact

This software enables us to perform extensive research with the flow based
sampling algorithm in two dimension and three dimension on CPU and GPU, which
is preferred since modern computer clusters and supercomputers have GPU
accelerator. More and more work with the flow-based sampling algorithm is
being published, and most of them are based on a Python code [25]. We provide
a new code in the Julia language with Docker environment. Thanks to the Julia
environment, this package can be run on a large-scale supercomputer or a
cluster with GPU if the machine supports Julia language.

This package also has the advantage of being used for prototyping. If one
wants to change the theory, one can achieve it by changing $src/actions.jl$.
Differentiation can be done by Flux.jl/Zygote.jl automatically as is explained
above.

##  5 Conclusions

We implement a flow-based sampling algorithm for two or three-dimensional
lattice scalar field theory in Julia programming language. This package
extensively uses automatic differentiation in Flux.jl/Zygote.jl to calculate
Jacobian in the flow based sampling algorithm. In addition, this package
provides, HMC, analysis tools for autocorrelation and ESS (effective sample
size) calculation with Jupyter notebook. A Docker image is provided, which
enables us to perform simulations without effort on environment construction
and this package works not only on the CPU but also on the NVIDIA GPU.
Hyperparameters and parameters including dimensionality can be specified in
the a parameter file.

This package pave a way of research directions not only the reduction of
autocorrelation time [19, 20, 22], scaling study [26], and investigation with
sign problem [33].

##  6 Conflict of Interest

There are no known conflicts of interest associated with this publication and
there has not been significant financial support for this work that could have
influenced its outcome.

## Acknowledgements

Authors thank to a public notebook [25], and we refer the code for the
implementation. This work of AT was supported by JSPS KAKENHI Grant Number
JP20K14479, JP22H05112, JP22H05111, and JP22K03539. This work was supported by
MEXT as “Program for Promoting Researches on the Supercomputer Fugaku”
(Simulation for basic science: approaching the new quantum era; Grant Number
JPMXP1020230411). This work was supported by MEXT as “Program for Promoting
Researches on the Supercomputer Fugaku” (Search for physics beyond the
standard model using large-scale lattice QCD simulation and development of AI
technology toward next-generation lattice QCD; Grant Number JPMXP1020230409).

## References

  * [1] R. Frezzotti, P. A. Grassi, S. Sint, P. Weisz, Lattice QCD with a chirally twisted mass term, JHEP 08 (2001) 058.  arXiv:hep-lat/0101001, doi:10.1088/1126-6708/2001/08/058.
  * [2] S. Capitani, M. Lüscher, R. Sommer, H. Wittig, Non-perturbative quark mass renormalization in quenched lattice QCD, Nucl. Phys. B 544 (1999) 669–698, [Erratum: Nucl.Phys.B 582, 762–762 (2000)].  arXiv:hep-lat/9810063, doi:10.1016/S0550-3213(98)00857-8.
  * [3] D. E. Groom, et al., Review of particle physics. Particle Data Group, Eur. Phys. J. C 15 (2000) 1–878.
  * [4] Y. Aoki, et al., FLAG Review 2021 (11 2021).  arXiv:2111.09849.
  * [5] C. Gattringer, C. B. Lang, Quantum chromodynamics on the lattice, Vol. 788, Springer, Berlin, 2010.  doi:10.1007/978-3-642-01850-3.
  * [6] S. Durr, et al., Ab-Initio Determination of Light Hadron Masses, Science 322 (2008) 1224–1227.  arXiv:0906.3599, doi:10.1126/science.1163233.
  * [7] A. Bazavov, et al., Equation of state in ( 2+1 )-flavor QCD, Phys. Rev. D 90 (2014) 094503.  arXiv:1407.6387, doi:10.1103/PhysRevD.90.094503.
  * [8] S. Borsanyi, G. Endrodi, Z. Fodor, A. Jakovac, S. D. Katz, S. Krieg, C. Ratti, K. K. Szabo, The QCD equation of state with dynamical quarks, JHEP 11 (2010) 077.  arXiv:1007.2580, doi:10.1007/JHEP11(2010)077.
  * [9] T. Aoyama, et al., The anomalous magnetic moment of the muon in the Standard Model, Phys. Rept. 887 (2020) 1–166.  arXiv:2006.04822, doi:10.1016/j.physrep.2020.07.006.
  * [10] S. Schaefer, R. Sommer, F. Virotta, Critical slowing down and error analysis in lattice QCD simulations, Nucl. Phys. B 845 (2011) 93–119.  arXiv:1009.5228, doi:10.1016/j.nuclphysb.2010.11.020.
  * [11] M. Luscher, Trivializing maps, the Wilson flow and the HMC algorithm, Commun. Math. Phys. 293 (2010) 899–919.  arXiv:0907.5491, doi:10.1007/s00220-009-0953-7.
  * [12] H. Nicolai, On a New Characterization of Scalar Supersymmetric Theories, Phys. Lett. B 89 (1980) 341.  doi:10.1016/0370-2693(80)90138-0.
  * [13] D. Boyda, et al., Applications of Machine Learning to Lattice Quantum Field Theory, in: 2022 Snowmass Summer Study, 2022.  arXiv:2202.05838.
  * [14] A. Tanaka, A. Tomiya, Towards reduction of autocorrelation in HMC by machine learning (12 2017).  arXiv:1712.03893.
  * [15] J. M. Pawlowski, J. M. Urban, Reducing Autocorrelation Times in Lattice Simulations with Generative Adversarial Networks, Mach. Learn. Sci. Tech. 1 (2020) 045011.  arXiv:1811.03533, doi:10.1088/2632-2153/abae73.
  * [16] K. Zhou, G. Endrodi, L.-G. Pang, H. Stöcker, Generative Model Study for 1+1d-Complex Scalar Field Theory, PoS AISIS2019 (2020) 007.  doi:10.22323/1.372.0007.
  * [17] A. Tomiya, Y. Nagai, Gauge covariant neural network for 4 dimensional non-abelian gauge theory (3 2021).  arXiv:2103.11965.
  * [18] S. Foreman, X.-Y. Jin, J. C. Osborn, LeapfrogLayers: A Trainable Framework for Effective Topological Sampling, PoS LATTICE2021 (2022) 508.  arXiv:2112.01582, doi:10.22323/1.396.0508.
  * [19] M. S. Albergo, G. Kanwar, P. E. Shanahan, Flow-based generative models for Markov chain Monte Carlo in lattice field theory, Phys. Rev. D 100 (3) (2019) 034515.  arXiv:1904.12072, doi:10.1103/PhysRevD.100.034515.
  * [20] G. Kanwar, M. S. Albergo, D. Boyda, K. Cranmer, D. C. Hackett, S. Racanière, D. J. Rezende, P. E. Shanahan, Equivariant flow-based sampling for lattice gauge theory, Phys. Rev. Lett. 125 (12) (2020) 121601.  arXiv:2003.06413, doi:10.1103/PhysRevLett.125.121601.
  * [21] D. Boyda, G. Kanwar, S. Racanière, D. J. Rezende, M. S. Albergo, K. Cranmer, D. C. Hackett, P. E. Shanahan, Sampling using $SU(N)$ gauge equivariant flows, Phys. Rev. D 103 (7) (2021) 074504.  arXiv:2008.05456, doi:10.1103/PhysRevD.103.074504.
  * [22] M. S. Albergo, G. Kanwar, S. Racanière, D. J. Rezende, J. M. Urban, D. Boyda, K. Cranmer, D. C. Hackett, P. E. Shanahan, Flow-based sampling for fermionic lattice field theories, Phys. Rev. D 104 (11) (2021) 114507.  arXiv:2106.05934, doi:10.1103/PhysRevD.104.114507.
  * [23] R. Abbott, et al., Gauge-equivariant flow models for sampling in lattice field theories with pseudofermions (7 2022).  arXiv:2207.08945.
  * [24] R. Abbott, et al., Sampling QCD field configurations with gauge-equivariant flow models, in: 39th International Symposium on Lattice Field Theory, 2022\.  arXiv:2208.03832.
  * [25] M. S. Albergo, D. Boyda, D. C. Hackett, G. Kanwar, K. Cranmer, S. Racanière, D. J. Rezende, P. E. Shanahan, Introduction to Normalizing Flows for Lattice Field Theory (1 2021).  arXiv:2101.08176.
  * [26] L. Del Debbio, J. M. Rossney, M. Wilson, Efficient modeling of trivializing maps for lattice $\phi$4 theory using normalizing flows: A first look at scalability, Phys. Rev. D 104 (9) (2021) 094507.  arXiv:2105.12481, doi:10.1103/PhysRevD.104.094507.
  * [27] S. Foreman, T. Izubuchi, L. Jin, X.-Y. Jin, J. C. Osborn, A. Tomiya, HMC with Normalizing Flows, PoS LATTICE2021 (2022) 073.  arXiv:2112.01586, doi:10.22323/1.396.0073.
  * [28] S. Duane, A. D. Kennedy, B. J. Pendleton, D. Roweth, Hybrid Monte Carlo, Phys. Lett. B 195 (1987) 216–222.  doi:10.1016/0370-2693(87)91197-X.
  * [29] M. A. Clark, A. D. Kennedy, Accelerating dynamical fermion computations using the rational hybrid Monte Carlo (RHMC) algorithm with multiple pseudofermion fields, Phys. Rev. Lett. 98 (2007) 051601.  arXiv:hep-lat/0608015, doi:10.1103/PhysRevLett.98.051601.
  * [30] M. Luscher, Schwarz-preconditioned HMC algorithm for two-flavour lattice QCD, Comput. Phys. Commun. 165 (2005) 199–220.  arXiv:hep-lat/0409106, doi:10.1016/j.cpc.2004.10.004.
  * [31] J. Bezanson, S. Karpinski, V. B. Shah, A. Edelman, Julia: A fast dynamic language for technical computing (2012).  doi:10.48550/ARXIV.1209.5145.
URL https://arxiv.org/abs/1209.5145

  * [32] J. contributors, Julia Micro-Benchmarks, https://julialang.org/benchmarks/, [Online; accessed 18-August-2022] (2008).
  * [33] M. Rodekamp, E. Berkowitz, C. Gäntgen, S. Krieg, T. Luu, J. Ostmeyer, Mitigating the Hubbard Sign Problem with Complex-Valued Neural Networks (3 2022).  arXiv:2203.00390.


# Differentiating Through Linear Solvers

Jan Hückelheim  Paul Hovland

###### Abstract

Computer programs containing calls to linear solvers are a known challenge for
automatic differentiation. Previous publications advise against
differentiating through the low-level solver implementation, and instead
advocate for high-level approaches that express the derivative in terms of a
modified linear system that can be solved with a separate solver call. Despite
this ubiquitous advice, we are not aware of prior work comparing the accuracy
of both approaches. With this article we thus empirically study a simple
question: What happens if we ignore common wisdom, and differentiate through
linear solvers?

Automatic Differentiation, Linear Solvers

##  1 Introduction

One of the merits of automatic, or algorithmic, differentiation (AD or
autodiff) is that it is applied directly to algorithms, applying rule-based
differentiation to the low-level functions intrinsic to a particular
programming language and combining those derivatives according to the chain
rule of differential calculus (Griewank & Walther, 2008). Nonetheless, when
applying AD to programs including iterative solvers, conventional wisdom
advises treating the solver as an intrinsic function, rather than “naively”
differentiating through the low-level implementation details. This approach is
particularly common in the case of fixed point iterations, where there is a
wealth of established theory (Christianson, 1994; Griewank et al., 1993) and
many AD tools support high-level differentiation of fixed-point iterative
solvers that have labeled as such.

The situation for linear solvers is more complicated. While stationary
iterative methods adhere to the theory of fixed point iterations, Krylov
methods do not, and the theory in this regime is much less developed. Gratton
et al. (2014) and Christianson (2018) have advanced the theory for conjugate
gradient solvers (see Section 2 for more details). However, it is unclear how
much of this theory carries over to solvers for nonsymmetric systems, such as
GMRES and BiCGStab. Furthermore, even in cases where theory argues for high-
level differentiation, practical considerations such as difficulty in clearly
delineating the solver interface or problems finding an appropriate
preconditioner for the transposed system might argue in favor of low-level
differentiation.

This paper examines empirically the tradeoffs between low-level
differentiation and high-level differentiation of several linear solvers, as
implemented in SPARSKIT (Saad, 1994). The remainder of this paper is organized
as follows. The next section provides additional background on low-level and
high-level differentiation of linear solvers. Section 3 discusses our
experimental design, followed by our experimental results in Section 4. We
conclude with a summary of our findings and our plans for future work.

##  2 Background

Automatic differentiation has two main “modes,” forward (or tangent) and
reverse (or adjoint) mode. In the so-called forward mode, the chain rule is
applied starting with the independent variables and propagating derivatives
forward to the dependent variables. In the so-called reverse mode, the chain
rule is applied starting with the dependent variables and propagating backward
to the independent variables. For this paper, we will focus primarily on the
forward mode.

Conjugate gradient is the best-studied Krylov solver with respect to low-level
differentiation. Gratton et al. (2014) examined the consequences of
differentiating through a conjugate gradient (CG) solver, under the assumption
that all of the iterates $x_{k}$ are well-defined, that is, CG does not reach
the exact solution $A^{-1}b$ before iteration $k$. Christianson (2018) also
considered low-level differentiation of CG, but in the case where the exact
solution is in fact found, proving that in exact arithmetic and with
appropriate modifications, both the solution to the linear system and its
derivatives can be computed in no more than $N$ iterations, for
$A\in\mathcal{R}^{N\times N}$. Low-level differentiation of Krylov solvers for
nonsymmetric systems is at present poorly understood at a theoretical level.

The alternative to low-level differentiation is high-level differentiation,
treating the linear solve as an elementary function and applying the rules of
matrix calculus to compute the derivatives directly. The derivation of such
rules is straightforward for the forward mode. For independent variable $u$
and since $Ax=b$, we have

$\frac{\partial A}{\partial u}x+A\frac{\partial x}{\partial u}=\frac{\partial
b}{\partial u}.$

Consequently, one can compute $\frac{\partial x}{\partial u}$ by solving the
linear system $Ay=\hat{b}$, where $\frac{\partial x}{\partial u}=y$ and

$\hat{b}=\frac{\partial b}{\partial u}-\frac{\partial A}{\partial u}x.$

Because the matrix $A$ is used in both linear solves, we can reuse the
preconditioner used to solve for $x$. The situation is more complicated for
the reverse mode and nonsymmetric matrices, where the linear solve used to
compute derivatives involves the matrix transpose, $A^{T}$.

##  3 Experimental Setup

We implemented low-level differentiation by applying the Tapenade AD tool
(Hascoet & Pascual, 2013) to the SPARSKIT implementation of GMRES, TFQMR and
BiCGStab. We modified the SPARSKIT implementation to replace Fortran language
features that caused compilation errors after differentiation with Tapenade,
such as computed and assigned goto statements, and created a head routine that
implements the choice of linear solver with a branch instead of passing the
name of an external subroutine as an argument. We also replaced the BLAS
implementation in SPARSKIT with a more recent version (Blackford et al.,
2002). We then implemented high-level differentiation of these solvers at the
matrix calculus level, following the derivation in Section 2.

We evaluated both approaches using 65 matrices from the SuiteSparse collection
(Davis & Hu, 2011) (see A for the list), comprising all of the nonsymmetric
matrices from the Bai collection with up to $1,000$ rows and columns. To
facilitate our assessment of the effectiveness of both strategies, we
manufactured the righthand sides $b$ and $\frac{\partial b}{\partial u}$ from
known $x$ and $\frac{\partial x}{\partial u}$, each with elements drawn from
the uniform distribution $[0,1)$. We monitored the progress after $k$
iterations, for $k\in[1,50]$, computing the L2 norm of the difference between
the computed $x$ (or $\frac{\partial x}{\partial u}$) and the reference value.

##  4 Experimental Results

Figures 1 and 2 show the effectiveness of the two differentiation strategies
in comparison with the original, undifferentiated linear solver. For the
original solver, the curve shows the number of problems solved (out of 65)
such that the L2 norm of the difference between the computed $x$ and the
reference value for $x$ is less than some threshold ($10^{-2}$ or $10^{-4}$).
For the two differentiation strategies, the curves show the number of problems
solved such that the L2 norm of the difference between the computed
$\frac{\partial x}{\partial u}$ an the reference value is less than the same
threshold. A higher curve indicates that more problems are solved within a
given number of iterations and a curve further to the left indicates that a
given number of problems are solved in fewer iterations.

A few patterns emerge from our results. First, none of the linear solvers is
able to solve more than two thirds of the problems to within $10^{-4}$ of the
reference solution in 2000 iterations, and restarted GMRES significantly
underperforms TFQMR and BiCGStab. Second, in general high-level
differentiation performs nearly as well as the original solver, but there are
typically a few problems that require more iterations to achieve similar
levels of accuracy. Third, the effectiveness of low-level differentiation is
highly solver-dependent. For TFQMR and restarted GMRES, the low-level
differentiation strategy is nearly as effective as high-level differentiation.
In contrast, for BiCGStab, there is a sizable gap in performance between high-
level and low-level differentiation, with the latter typically only able to
solve half as many problems within a similar number of iterations.

Figure 1: Data profile showing the relative performance of the low-level and
high-level differentiation strategies. The curve for the original linear
solver shows the number of problems solved (out of 65) such that
$\left\|x-x_{\mbox{\scriptsize ref}}\right\|_{2}<10^{-2}$. The curves for the
differentiated linear solvers show the number of problems solved such that
$\left\|\frac{\partial x}{\partial u}-\left(\frac{\partial x}{\partial
u}\right)_{\mbox{\scriptsize ref}}\right\|_{2}<10^{-2}$. Figure 2: Data
profile showing the relative performance of the low-level and high-level
differentiation strategies. The curve for the original linear solver shows the
number of problems solved (out of 65) such that $\left\|x-x_{\mbox{\scriptsize
ref}}\right\|_{2}<10^{-4}$. The curves for the differentiated linear solvers
show the number of problems solved such that $\left\|\frac{\partial
x}{\partial u}-\left(\frac{\partial x}{\partial u}\right)_{\mbox{\scriptsize
ref}}\right\|_{2}<10^{-4}$.

Figures 3 and 4 show detailed results for two matrices, namely BFWA398 and
BFWA62 from the collection. These matrices emerge from computational
electromagnetics problems, and have reasonable condition numbers of
$2.993111e+03$ and $5.530615e+02$, respectively. The original linear solvers
therefore converge well, as do the high-level differentiated solvers. Low-
level differentiation through the solvers leads to very erratic behavior for
BiCGStab, and divergence for GMRES in both cases. For TFQMR both
differentiation strategies appear to work well.

Figure 3: Convergence for the original solver and both differentiation
strategies when applied to the BFWA62 matrix. The curves show the error in the
system solution as a function of the iterations performed, where a lower value
is better. The BICGStab and TFQMR solvers rapidly reach a value close to
machine precision and remain relatively stable at that lavel. GMRES converges
more slowly. High-level differentiation broadly follows this trend for all
solvers. Low-level differentiation performs more erratically, and consistently
underperforms high-level differentiation. Nevertheless, low-level
differentiation appears to work reasonably well for TFQMR. Figure 4:
Convergence for the original solver and both differentiation strategies when
applied to the BFWA398 matrix. Similar to the results in Figure 3, BICGStab
and TFQMR rapidly converge and outperform GMRES. High-level differentiation
performs worse than the original solver, but far better than low-level
differentiation. Once again, TFQMR appears to be better suited for low-level
differentiation and leads to reasonable results for both differentiation
approaches.

##  5 Conclusions

We demonstrate with this article that the common advice is justified, and that
high-level differentiation is indeed usually preferable to low-level
differentiation. However, our work also shows a more nuanced picture: some
linear solvers appear to be better suited for low-level differentiation than
others, and in situations where high-level differentiation is not desirable –
for example due to the manual development effort that is often involved – a
careful solver choice may lead to useful gradients even with low-level
approaches.

In future work, we plan to investigate reverse mode differentiation through
solvers, extend the collection of solvers and preconditioners, and compare
empirical results with the theory presented in previous work (Christianson,
2018). We are also interested in quantifying the influence of roundoff errors
on our results.

##  6 Acknowledgements

This work was partly supported by the Applied Mathematics activity within the
U.S. Department of Energy, Office of Science, Advanced Scientific Computing
Research Program, under contract number DE-AC02-06CH11357.

## References

  * Blackford et al. (2002) Blackford, L. S., Petitet, A., Pozo, R., Remington, K., Whaley, R. C., Demmel, J., Dongarra, J., Duff, I., Hammarling, S., Henry, G., et al.  An updated set of basic linear algebra subprograms (blas).  _ACM Transactions on Mathematical Software_ , 28(2):135–151, 2002.
  * Christianson (1994) Christianson, B.  Reverse accumulation and attractive fixed points.  _Optimization Methods and Software_ , 3(4):311–326, 1994.
  * Christianson (2018) Christianson, B.  Differentiating through conjugate gradient.  _Optimization Methods and Software_ , 33(4-6):988–994, 2018.
  * Davis & Hu (2011) Davis, T. A. and Hu, Y.  The university of florida sparse matrix collection.  _ACM Trans. Math. Softw._ , 38(1), dec 2011.  ISSN 0098-3500.  doi: 10.1145/2049662.2049663.  URL https://doi.org/10.1145/2049662.2049663.
  * Gratton et al. (2014) Gratton, S., Titley-Peloquin, D., Toint, P., and Ilunga, J. T.  Differentiating the method of conjugate gradients.  _SIAM Journal on Matrix Analysis and Applications_ , 35(1):110–126, 2014.
  * Griewank & Walther (2008) Griewank, A. and Walther, A.  _Evaluating derivatives: principles and techniques of algorithmic differentiation_.  SIAM, 2008.
  * Griewank et al. (1993) Griewank, A., Bischof, C., Corliss, G., Carle, A., and Williamson, K.  Derivative convergence for iterative equation solvers.  _Optimization methods and software_ , 2(3-4):321–355, 1993.
  * Hascoet & Pascual (2013) Hascoet, L. and Pascual, V.  The tapenade automatic differentiation tool: Principles, model, and specification.  _ACM Transactions on Mathematical Software (TOMS)_ , 39(3):1–43, 2013.
  * Saad (1994) Saad, Y.  SPARSKIT: a basic tool kit for sparse matrix computations, 1994.

##  Appendix A List of SuiteSparse Matrices

The following are the SuiteSparse (Davis & Hu, 2011) matrices used in our
experiments.

Id | Name | Group | Rows | Cols | Nonzeros | Kind | Date
---|---|---|---|---|---|---|---
292 | bfwa398 | Bai | 398 | 398 | 3,678 | Electromagnetics Problem | 1994
293 | bfwa62 | Bai | 62 | 62 | 450 | Electromagnetics Problem | 1994
294 | bfwa782 | Bai | 782 | 782 | 7,514 | Electromagnetics Problem | 1994
298 | bwm200 | Bai | 200 | 200 | 796 | Chemical Process Simulation Problem | 1992
299 | bwm2000 | Bai | 2,000 | 2,000 | 7,996 | Chemical Process Simulation Problem | 1992
300 | cdde1 | Bai | 961 | 961 | 4,681 | Computational Fluid Dynamics Problem Sequence | 1992
301 | cdde2 | Bai | 961 | 961 | 4,681 | Subsequent Computational Fluid Dynamics Problem | 1992
302 | cdde3 | Bai | 961 | 961 | 4,681 | Subsequent Computational Fluid Dynamics Problem | 1992
303 | cdde4 | Bai | 961 | 961 | 4,681 | Subsequent Computational Fluid Dynamics Problem | 1992
304 | cdde5 | Bai | 961 | 961 | 4,681 | Subsequent Computational Fluid Dynamics Problem | 1992
305 | cdde6 | Bai | 961 | 961 | 4,681 | Subsequent Computational Fluid Dynamics Problem | 1992
306 | ck104 | Bai | 104 | 104 | 992 | 2D/3D Problem | 1986
307 | ck400 | Bai | 400 | 400 | 2,860 | 2D/3D Problem | 1986
308 | ck656 | Bai | 656 | 656 | 3,884 | 2D/3D Problem | 1986
309 | dw1024 | Bai | 2,048 | 2,048 | 10,114 | Electromagnetics Problem | 1993
310 | dw256A | Bai | 512 | 512 | 2,480 | Electromagnetics Problem | 1993
311 | dw256B | Bai | 512 | 512 | 2,500 | Electromagnetics Problem | 1993
312 | dw4096 | Bai | 8,192 | 8,192 | 41,746 | Electromagnetics Problem | 1993
313 | lop163 | Bai | 163 | 163 | 935 | Statistical/Mathematical Problem | 1978
314 | mhda416 | Bai | 416 | 416 | 8,562 | Electromagnetics Problem | 1994
316 | odepa400 | Bai | 400 | 400 | 1,201 | 2D/3D Problem | 1978
318 | olm100 | Bai | 100 | 100 | 396 | Computational Fluid Dynamics Problem | 1994
319 | olm1000 | Bai | 1,000 | 1,000 | 3,996 | Computational Fluid Dynamics Problem | 1994
320 | olm2000 | Bai | 2,000 | 2,000 | 7,996 | Computational Fluid Dynamics Problem | 1994
321 | olm500 | Bai | 500 | 500 | 1,996 | Computational Fluid Dynamics Problem | 1994
322 | olm5000 | Bai | 5,000 | 5,000 | 19,996 | Computational Fluid Dynamics Problem | 1994
323 | pde225 | Bai | 225 | 225 | 1,065 | 2D/3D Problem | 1982
324 | pde2961 | Bai | 2,961 | 2,961 | 14,585 | 2D/3D Problem | 1982
325 | pde900 | Bai | 900 | 900 | 4,380 | 2D/3D Problem | 1982
328 | qh882 | Bai | 882 | 882 | 3,354 | Power Network Problem | 1994
329 | rbsa480 | Bai | 480 | 480 | 17,088 | Robotics Problem | 1993
330 | rbsb480 | Bai | 480 | 480 | 17,088 | Robotics Problem | 1993
331 | rdb2048 | Bai | 2,048 | 2,048 | 12,032 | Computational Fluid Dynamics Problem | 1994
332 | rdb5000 | Bai | 5,000 | 5,000 | 29,600 | Computational Fluid Dynamics Problem | 1994
333 | rdb968 | Bai | 968 | 968 | 5,632 | Computational Fluid Dynamics Problem | 1994
334 | rw136 | Bai | 136 | 136 | 479 | Statistical/Mathematical Problem | 1978
335 | rw496 | Bai | 496 | 496 | 1,859 | Statistical/Mathematical Problem | 1978
336 | rw5151 | Bai | 5,151 | 5,151 | 20,199 | Statistical/Mathematical Problem | 1978
337 | tub100 | Bai | 100 | 100 | 396 | Computational Fluid Dynamics Problem | 1994
338 | tub1000 | Bai | 1,000 | 1,000 | 3,996 | Computational Fluid Dynamics Problem | 1994
1612 | cryg10000 | Bai | 10,000 | 10,000 | 49,699 | Materials Problem | 1996
1613 | cryg2500 | Bai | 2,500 | 2,500 | 12,349 | Materials Problem | 1996
1614 | dw2048 | Bai | 2,048 | 2,048 | 10,114 | Electromagnetics Problem | 1993
1615 | dw8192 | Bai | 8,192 | 8,192 | 41,746 | Electromagnetics Problem | 1993
1616 | dwa512 | Bai | 512 | 512 | 2,480 | Electromagnetics Problem | 1993
1617 | dwb512 | Bai | 512 | 512 | 2,500 | Electromagnetics Problem | 1993
1620 | mhd1280a | Bai | 1,280 | 1,280 | 47,906 | Electromagnetics Problem | 1994
1622 | mhd3200a | Bai | 3,200 | 3,200 | 68,026 | Electromagnetics Problem | 1994
1624 | mhd4800a | Bai | 4,800 | 4,800 | 102,252 | Electromagnetics Problem | 1994
1626 | qh1484 | Bai | 1,484 | 1,484 | 6,110 | Power Network Problem | 1994
1627 | qh768 | Bai | 768 | 768 | 2,934 | Power Network Problem | 1994
1628 | rdb1250 | Bai | 1,250 | 1,250 | 7,300 | Computational Fluid Dynamics Problem | 1994
1629 | rdb1250l | Bai | 1,250 | 1,250 | 7,300 | Computational Fluid Dynamics Problem | 1994
1630 | rdb200 | Bai | 200 | 200 | 1,120 | Computational Fluid Dynamics Problem | 1994
1631 | rdb200l | Bai | 200 | 200 | 1,120 | Computational Fluid Dynamics Problem | 1994
1632 | rdb2048_noL | Bai | 2,048 | 2,048 | 12,032 | Computational Fluid Dynamics Problem | 1994
1633 | rdb3200l | Bai | 3,200 | 3,200 | 18,880 | Computational Fluid Dynamics Problem | 1994
1634 | rdb450 | Bai | 450 | 450 | 2,580 | Computational Fluid Dynamics Problem | 1994
1635 | rdb450l | Bai | 450 | 450 | 2,580 | Computational Fluid Dynamics Problem | 1994
1636 | rdb800l | Bai | 800 | 800 | 4,640 | Computational Fluid Dynamics Problem | 1994
1637 | tols1090 | Bai | 1,090 | 1,090 | 3,546 | Computational Fluid Dynamics Problem | 1991
1638 | tols2000 | Bai | 2,000 | 2,000 | 5,184 | Computational Fluid Dynamics Problem | 1991
1639 | tols340 | Bai | 340 | 340 | 2,196 | Computational Fluid Dynamics Problem | 1991
1640 | tols4000 | Bai | 4,000 | 4,000 | 8,784 | Computational Fluid Dynamics Problem | 1991
1641 | tols90 | Bai | 90 | 90 | 1,746 | Computational Fluid Dynamics Problem | 1991

Table 1: Matrices used in our experiments.