A Monte Carlo algorithm for efficient large matrix inversion

L. A. Garcı́a-Cortés ∗†and C. Cabrillo ‡

arXiv:cs/0412107v2 [[Link]] 10 Jan 2005

3rd April 2018

Abstract
This paper introduces a new Monte Carlo algorithm to invert large matrices. It is based on
simultaneous coupled draws from two random vectors whose covariance is the required inverse.
It can be considered a generalization of a previously reported algorithm for hermitian matrices
inversion based in only one draw. The use of two draws allows the inversion on non-hermitian
matrices. Both the conditions for convergence and the rate of convergence are similar to the
Gauss-Seidel algorithm. Results on two examples are presented, a real non-symmetric matrix
related to quantitative genetics and a complex non-hermitian matrix relevant for physicists.
Compared with other Monte Carlo algorithms it reveals a large reduction of the processing
time showing eight times faster processing in the examples studied.
Keywords: Monte Carlo, inverse matrix, Gibbs sampler, sparse matrices.

1 Introduction
The computation of the inverse of a matrix or some of its elements is one of the main topics in
numerical analysis. Large sparse matrices can be usually inverted from its factors obtained by
using sparse matrix techniques [11]. Although it can be alleviated by reordering rows and columns,
these techniques suffer from the fill-in required by the LU factorization and they usually result in
algorithms demanding large random access memory (RAM) values.
Monte Carlo algorithms are used to estimate empirically the expectation of random variables.
Using them to invert matrices consist basically in obtaining and averaging several realized values
of a variable whose expectation is the required inverse. Although Monte Carlo algorithms present
a stochastic error, they can be used to obtain estimates of the inverse when the matrix is too large
or complicate to be inverted by using conventional or sparse algorithms. These algorithms have
high parallel efficiency, they do not demand a great amount of RAM and the processing time is
proportional to the number of nonzero elements of the matrix to be inverted. In general, they can
be very efficient when rough estimates of the inverse matrix are required.
The most common Monte Carlo algorithms to invert matrices are based on finite discrete
Markov chains [5]. In these algorithms several random trajectories are calculated to solve a linear
system of equations for each matrix column. Alternatively, there are Monte Carlo methods based
on random variables whose expectation is the whole inverse instead one of its columns. Typically
the rational is to consider the sought for inverse matrix C−1 as the covariance matrix of a normally
distributed random vector,

z ∼ N 0, C−1 . (1)
After generating by any Monte Carlo algorithm sampled values of z (something for which C−1
is not needed but just C), the inverse matrix is obtained by an estimation of the variances and
covariances of the sample.
∗ Dpto. de Mejora Genética, CIT-INIA, Carretera de La Corua, Km. 7,5 Madrid, E-28040, Spain. e-mail:

garcia@[Link]
† Dpto. de Produccin Animal. ETSI Agrnomos. Universidad Politcnica de Madrid
‡ Instituto de Estructura de la Materia, Serrano 123, Madrid, E-28006, Spain. e-mail: ccabrilo@[Link]

1
 Recently the Gibbs sampler [9, 10] algorithm has been proposed to draw z values from (1)
using the elements of C [8, 14]. The procedure consists in updating successively the elements of z
for several cycles by using the formula
i−1 n
(k) (k) 1 1 X (k) 1 X (k−1)
zi = φi √ − zj cij − z cij (2)
cii cii j=1 cii j=i+1 j

where the subscript i means the ith element of z, the superscript (k) means the k th cycle and φ is
an independent random noise drawn from a standard normal distribution. After a given period
of
(k) (k) †(k) −1
convergence
 usually
 called the burn-in period, z will satisfy at each cycle E z z = C
or E z†(k) Qz(k) = trQC−1 (the † denotes hermitian conjugate, i.e., transposition and complex
conjugation). These expectations will provide the elements of the inverse or linear functions of
them respectively. Hereafter, we will denote this algorithm by GS. Although GS does not provide
an independent but a serially correlated set of samples, the convergence is fast and, contrary to
Metropolis-like algorithms, it profits of all the draws.
The obvious disadvantage of the algorithms based in Eq. (1) is a rather narrow range of
applicability because of the necessary condition of C being positive definite. For many statistical
applications this is not a problem [14] but it is a serious drawback in Physics applications (for
instance in Lattice QCD [17]). The problem can be bypassed by inverting C† C and applying
the result to C but this procedure seriously reduces the efficiency. Rather surprisingly this is not
the case for the GS algorithm. Once recourse is made to the Gibbs sampler to implement the
covariance matrix rational, the convergence of the algorithm, at least if it is explicitly written as
in Eq. (2), is governed by the sampling method and not by the implications of (1). As we will
see, the convergence conditions are then more flexible but stillCmust be symmetric (or hermitian
if defined over the complex numbers).
A different approach usually known as stochastic estimation (SE from now on), avoid this
problem [13] by recourse to the repeated solution of linear systems like in the discrete Markov
Chain methods mentioned above but this time the number of systems to be solved does not
depend on the rank of the matrix to be inverted. The rational is to solve a series of linear

systems
†
Cv = φ where φ is drawn from a random variable satisfying
E (φ) = 0 and E φφ = I. The
inverse is then given by taking the ensemble average E vφ† = C−1 . Rather surprisingly it was
not recognized until recently that this method do not rely on the gaussianity of φ. Thus, Dong
and Liu [6] and, independently, Garcı́a-Cortés [7] have shown that the efficiency of the method
improves substantially when is drawn as φi = 2B − 1, where B is a Bernoulli random variable (a
Z 2 noise).
When applicable, however, it is expected that GS will be more efficient than SE. This comes
from the fact that GS mimics Gauss-Seidel algorithm but including noise at each cycle [18]. The
processing time to obtain each realized value of zz† or z† Qz is then slightly greater than the
processing time of a Gauss-Seidel iteration, while the processing time to obtain a SE realized
value is equivalent to the processing time needed to solve a linear system.
In this paper we introduce a new ”noisy” Gauss-Seidel method which include as a particular
case the GS algorithm. Our algorithm, that we will call it the correlated chains sampling algorithm
(CC from now on), is applicable to nondefinite positive matrices and profits of the larger efficiency
of the Z 2 noise. The power of the algorithm is assessed with two different non-hermitian matrix,
a real one used in genetic improvement in animal breeding and a complex one typical of particle
Physics. Numerical tests show how our proposal can be around an order of magnitude faster than
SE. The sufficient and necessary conditions for the convergence of the method are also given.

2 The correlated chains sampling algorithm

The CC algorithm introduced in this section is based on the simultaneous update

of a couple of
random vectors z and w such as, after convergence of the algorithm, E zw† = C−1 being C the

2
matrix to be inverted. The algorithm updates simultaneously each element of z and w using
i−1 n
(k) (k) 1 1 X (k) 1 X (k−1)
zi = φi √ − zj cij − z cij (3)
cii cii j=1 cii j=i+1 j

i−1 n
(k) (k) 1 1 X (k) ∗ 1 X (k−1) ∗
wi = φi p ∗ − ∗ wj cji − ∗ w cji (4)
cii cii j=1 cii j=i+1 j

where Φ(k) is a set of independent noise vectors

 such  as
(k)
E Φ =0 (5)

and  
E Φ(k) Φ(l)† = Iδk,l . (6)

Note that the noise term in (3) and in (4) are both the same.
After discarding N cycles of a total of M as a period of convergence or burn-in, the Monte
Carlo estimation of the inverse can be obtained from,
M
X
1
C−1 ≃ z(k) w(k)† . (7)
M −N
k=N +1


Expression E z† Qz = trQC−1 also holds, and its Monte Carlo estimation via the CC algo-
rithm is given by
XM
1
−1
trQC ≃ z(k)† Qw(k) (8)
M −N
k=N +1

2.1 Demonstration
Let us start by rewriting Eqs. (3) and (4) in matrix form,
1 1 1
z(k) = √ Φ(k) − Lz(k) − Uz(k−1) , (9)
D D D
1 1 1
w(k) = √ Φ(k) − † U† w(k) − † L† w(k−1) , (10)
D † D D
where L, D and U are the lower triangle, the diagonal and the upper triangle of C respectively,
so that L + D + U = C. Solving for z(k) and w(k)† from Eqs.(9) and (10) we have,
1 √ 1
z(k) = D Φ(k) − Uz(k−1) , (11)
(D + L) D+L
√ 1 1
w(k)† = Φ(k)† D − w(k−1)† L . (12)
(D + U) D+U
After running the algorithm during n updating cycles (11) and (12) leads to,
n
X
z(n) = (−T)k Θ(n−k) , (13)
k=0
Xn
w(n)† = Γ(n−k)† (−S)k , (14)
k=0

3
where for the ease of the notation we have defined,
1
T = U, (15)
D+L
1
S = L , (16)
D+U
1 √
Θ(k) = D Φ(k) , (17)
(D + L)
1 √
Γ(k) = D† Φ(k) , (18)
(D + U)†

for all k except for k = 0 for which, again for the ease of notation, Θ(0) and Γ(0) denote just
the initial vectors z(0) and w(0) respectively. For a bounded noise such us the Z 2 noise, z(n) and
w(n) are bounded by
n
X n
X
(n) k 1 √ 1 √
z ≤ (−T) D B≤ D B Tk ,
(D + L) (D + L)
k=0 k=0
Xn √ √ n
X
1 1
w(n) ≤ (−S)k D B≤ D B Sk ,
(D + U) (D + U)
k=0 k=0

where B is the upper bound of the noise Φ and therefore the absolute convergence of the alternating
power series T and S warrants the convergence of the algorithm. In fact, this is a necessary and
sufficient condition for the convergence of the algorithm for any possible value of the noise chain.
To see this notice that the worst case corresponds to a noise such as the components of Tk and Sk
with the maximum absolute values attains positive signs. But the maximum of the absolute values
of the components of a matrix is a norm. The convergence of these components, which warrants
the convergence of any other of the components, corresponds therefore to the absolute convergence
of the series under such a norm and all the norms are equivalent as far as the convergence of matrix
series is concerned.
Given the convergence of the alternating power series in T and S we have finite z(∞) and w(∞)
and for a given integer N , the z(n) and w(n) corresponding to n ≥ N can be split in series of N
terms plus and arbitrarily small remainder for a sufficiently large N ,
N
X −1
(n) N +1
z = (−T)k Θ(n−k) + O(kTk ),
k=0
N
X −1
w(n)† = Γ(n−k)† (−S)k + O(kSkN +1 ) .
k=0

The sample average of the zw† product discarding N burn-in cycles of a total of M ≫ N , again
retaining explicitly only terms up to order N , gives
M N
X X −1 N
X −1
1
hzw† i = (−1)k+j Tk Θ(n−k) Γ(n−j)† Sj +
M −N j=0
n=N k=0

O(kTkN +1 , kSkN +1 ) =
N
X −1 N
X −1 M
X
1
(−1)k+j Tk Θ(n−k) Γ(n−j)† Sj +
M −N
k=0 j=0 n=N

O(kTkN +1 , kSkN +1 ) , (19)

Taking the limit of large M , Eq. (19) yields,

4
 lim hzw† i = E(zw† ) =
M→∞
N
X −1 N
X −1  
N +1 N +1
(−1)k+j Tk E Θ(n−k) Γ(n−j)† Sj + O(kTk , kSk ) =
k=0 j=0
N
X −1 N
X −1
1 √ √ 1 N +1 N +1
(−1)k+j Tk D Iδk,j D Sj + O(kTk , kSk ) =
D+L D+U
k=0 j=0
N
X −1
1 1 N +1 N +1
Tk D Sk + O(kTk , kSk ). (20)
D+L D+U
k=0

In the limit of large N , Eq. (20) reduces to

∞
X 1 1
lim E(zw† ) = Tk D Sk . (21)
N →∞ D+L D+U
k=0

Equation (21) represents the stationary value of the zw† average. That the series in (21)
approach C−1 , when they converge, is easily verified by iteration of the following recursive formula,
1 1 1 1
= D + T S, (22)
C D+L D+U C
which can be derived from the identity,
1 1
(C − U) (C − L) = C − (U + L) + U L ,
C C
since for the particular case in which C = D + U + L it reduces to,
1 1
(D + L) (D + U) = D + U L .
C C
Making use of the definitions (15) and (16) formula (22) follows trivially.

2.2 Convergence analysis

As shown in the previous section, the convergence of the algorithm is determined by the absolute
convergence of the series,
∞
X
(−T)k , (23)
k=0
∞
X
(−S)k . (24)
k=0

We therefore center in searching for the necessary and sufficient conditions for the convergence of
the series (23) and (24). Restricting the analysis to the series (23) as the same rational applies to
(24), first notice that the absolute convergence of (23) implies that the spectral radius of T, i.e.,
1/k
limk→∞ Tk ≡ sp(T) [2], is strictly below 1. This comes as a consequence of the Cauchy root
1/m
convergence test [1] which implies kTm k < 1 for some sufficiently large m if (23) converges.
1/m
Conversely if sp(T) < 1, then kTm k < 1 for some large enough m which in turn implies
kTm k < 1. We now split the series of the absolute values of (23) in a finite part with the terms
up to m − 1 and the rest, so that,
∞
X m−1
X ∞
X m−1
X ∞
X k
Tk = Tk + Tm+k ≤ Tk + kTm k .
k=0 k=0 k=0 k=0 k=0

5
Given that kTm k < 1the last term converges to 1/(1−kTk) and therefore (23) converges absolutely.
Applying the same rational to (24) we have then, that the algorithm converges in the sense that
both z(∞) and w(∞) are finite for any possible drawn values of the noise Φ, if and only if sp(T) < 1
and sp(S) < 1 . In strict sense, these are the necessary and sufficient conditions for the convergence
of Eqs. (11) and (12). The final algorithm is implemented using (9) and (10) and therefore also
the non-singularity of Dis needed.

2.3 The burn-in period and Monte Carlo error

As usual in this kind of algorithms, the efficiency of the code is increased by discarding a number of
iterations, N , called burn-in period, which are too influenced by the initial values. Given that the
CC path does not converge to a deterministic value but to a random variable an estimation of N
is not as easy as in a deterministic iterative procedures (for instance, the Gauss-Seidel algorithm).
To overcome the difficulties associated to the random nature of the algorithm we will use the so
called coupling method [16]. It consists in running a couple of paths for both z and w† with
different initial values but with the same random numbers for each couple. When the difference
between paths of a couple reach a given tolerance, the burn-in is assumed to be finished. From
(9) the difference between two paths z1 and z2 with different initial values but with the same Φ(k)
results in,
(k) (k) 1  (k) (k)
 1  (k−1) (k−1)

z2 − z1 = − L z2 − z1 − U z2 − z2 . (25)
D D
But Eq. (25) is nothing else that the formula corresponding to a Gauss-Seidel iterative algorithm
[20] for the solution of the linear system,

C (z2 − z1 ) = 0 ,

whose solution is z2 − z1 = 0, for a non-singular C. Not surprisingly, the condition for the
−1
convergence of the Gauss-Seidel algorithm represented by (25) is, sp((D + L) U) < 1 , that is,
sp(T) < 1. As expected, the same rational applied to Eqs. (10) leads to a Gauss-Seidel algorithm

−1 †
for C† (w2 − w1 ) = 0 with a convergence condition sp( D† + U† L ) < 1, that is, sp(S) < 1
. These two Gauss-Seidel procedures can be useful to analyze the convergence rate of the CC
algorithm in a given instance.
It is also relevant to calculate the Monte Carlo error of the expectations (7). The Monte Carlo
error is Gaussian because of the central limit theorem, but the realized values provided by the
CC algorithm are serially correlated and its Monte Carlo variance is bigger than the Monte Carlo
variance expected from an independent set of M − N samples. To calculate the Monte Carlo error
instead of M − N samples, we used a higher effective number of samples as proposed in [12].

2.4 Algorithm
Algorithm 1 shows the correlated chains algorithm including the determination of the burn-in
period (N ). During the burn-in period, four chains are computed until the paths converge to a
given tolerance, typically a number small enough to avoid the dependency of the chains on the
starting values. In order to simplify the algorithm outlined here, we consider a fixed chain length
(M ), i.e., we are assuming that the M − N cycles after burn-in are enough to obtain accurate
estimates of the elements of the inverse. Nevertheless, in the numerical tests presented in this
paper, we have included the determination of M − N by inserting the calculation of the effective
number of cycles during the iteration process.

Algorithm 1 Correlated chains algorithm including the determination of the burn-

in period
Given n, Q, C, B, tol
1. Set arbitrary starting values for z, z∗ , w′ and w′∗ , for instance zi =0, zi∗ =i, wi =0 and
wi∗ =i.

6
 2. Sample φ as a vector containing independent draws according with definition (5) and (6)
3. Update z, z∗ , w′ and w′∗ by using equations (3) and (4)
4. p ← p + 1
5. Go to step 2 until z′ z∗ < tol and w′ w∗ < tol
6. Sample φ as a vector containing independent draws according with definition (5) and (6)
7. Update z and w by using equations (3) and (4)
8. Accumulate z′ Qw in s
9. Go to step 6 to compute the next round of iteration (B − p times)
10. Set the final estimate: trQC−1 ← s/ (B − p)

3 Numerical tests
To check the computational efficiency of the proposed algorithm we will compare its performance
with respect to the efficiency of the SE algorithm which has shown itself as a efficient inver-
sion method [6]. We will make use of two different examples from two disparate fields: genetic
improvement in animal breeding and quantum field theory in the lattice.

3.1 Example 1. Wu and Schaeffer’s real asymmetric matrix

The Henderson’s mixed model equations [15] are routinely used in animal breeding to evaluate
the candidates to the artificial selection in livestock populations. This method takes into account
both the performance records of the animals in a given population and the pedigree relationships
between animals. In animal breeding, the quality of the estimations provided by a given model
corresponds to the diagonal elements of the inverse of a coefficient matrix. For instance, in the
case of the simplest model, these diagonal elements, known as the prediction error variances, take
the form  !−1 
X† X X†
diag  σ2 σe2  (26)
X I + A−1 σe2
a

where X is an incidence matrix mapping animals into herds. Its nonzero elements are xij = 1,
which indicates that animal i was recorded in herd j. I is the identity matrix with order equal
to the number of animals in the population. Variances σe2 and σa2 are assumed to be known and
correspond to the additive genetic variance and the residual variance. A is known as the numerator
relationship matrix and it maps the genetic relationships between animals. For instance, aij = 0.25
if animals i and j are half-sibs, aij = 0.5 if animals i and j are parent-progeny related, aii = 1
if the animal i is not inbred, etc. The inverse of the numerator relationship matrix A−1 can be
easily obtained by using the Henderson’s rules [15], but the whole coefficient matrix in equation
(26) has to be explicitly inverted.
The amount of animals used in a typical analysis is around hundreds of thousands and at least
a matrix row per animal is needed so that the rank of the coefficient matrix in Eq. (26) is of that
order. This coefficient matrix is positive definite and the Gibbs sampler based algorithm [14] can
be easily implemented. Nevertheless, the artificial selection based on the Henderson’s mixed model
equations tend to select animals coming from a small number of families and the percentage of
inbreeding increases significantly after a few generations of artificial selection. Wu and Schaeffer
[19] proposed an original method to select artificially the populations keeping the genetic gain
very close to the Henderson’s optimum but reducing significantly the rate of inbreeding. The
numerator relationship matrix has to be replaced by a customary asymmetric relationship matrix
Ae , and they provide simple rules to obtain A e −1 . For each animal i in the pedigree, with sire s
and dam d,
e −1 ,
1. Add ((1 − λ) δi + λ) to the (i, i) position of A
e −1 ,
2. Add − (1 − λ) δi /2 to the (i, s) and (i, d) positions of A

7
 Table 1: Topics of the Wu and Schaeffer’s coefficient matrix in two cases of different size.
Case 1 2
Number of animals 50000 100000
Number of herds 5000 10000
Rank of the coefficient matrix 55000 110000
Nonzero elements 424978 848982
σe2 /σa2 3.0 3.0
λ 0.2 0.2

e −1 ,
3. Add −δi /2 to the (s, i) and (d, i) positions of A
e −1 ,
4. Add −δi /4 to the (s, s), (d, d), (s, d) and (d, s) positions of A

where δi = 2 when both parents of the ith animal are known, δi = 4/3 when one parent is known
and δi = 1 when both parents are unknown. λ ∈ [0, 1] is a known coefficient which determine
the weight of the family records in each animal genetic evaluation, for instance, setting λ = 0
correspond to the conventional Henderson’s mixed model equations. These rules result in an
asymmetric coefficient matrix for λ 6= 0, non suitable for GS. Both SE and CC can be used to
compute the inverse required in expression (26) even in cases where A is replaced by its non-
symmetric counterpart A. e We will now compare the performance of our algorithm against the SE
estimator using two cases of Wu and Schaeffer’s coefficient matrices as described in Table 1.
The correlated chains sampling was implemented by setting the tolerance for burn-in as 5 10−5 .
After the burn-in period was finished, realized values of tr(z† Qw) were obtained for each cycle
and averaged. Convergence  after burn-in
 was checked every 100th iteration. The standard error
†
of the final estimate of E tr(z Qw) was obtained from the variance between realized values of
tr(z† Qw) and an effective chain length as described in [12]. Making recourse of this effective length
comes as a consequence of the correlated nature of the realized values in the CC algorithm. At
any rate, we tested the Geyer’s method against an empirical standard error estimate obtained by
replicating one hundred times the whole analysis varying the random seed number. The algorithm
was assumed to reach the convergence when the standard error yielded a relative error smaller
than the required tolerance of 5 10−5 .
SE method was implemented by solving the linear systems by both the iterative Gauss-Seidel
algorithm and the bi-conjugate gradient method. The latter resulted in a more efficient SE algo-
rithm in terms of CPU time and it is the only one presented in this manuscript.
Table 2 shows the results of both algorithms in cases 1 and 2. Relevant topics concerning
the computing efficiency of the CC algorithm are the number of rounds computed to reach the
burn-in tolerance, the number of rounds required to reach the required standard error of the final
estimate, and the CPU time required to compute both each burn-in cycle (four chains) and each
after burn-in cycle (two chains). Relevant topics for the SE algorithm are the average number
of rounds required to solve each linear system via the bi-conjugate gradient method (maximum
change between successive iterations was set to 5 10−5), the number of linear systems to be solved
in order to reach the required standard error of the final estimate and the CPU time per bi-
conjugate gradient round of iteration. All CPU times in Table 2 are referred to the time taken to
complete a cycle in the CC algorithm in the smallest case.
Table 2 shows also the agreement between the results provided by both the CC and the SE
algorithms. The CC algorithm is vastly more efficient in terms of the final CPU time required for
the computation of the inverse being around eight times faster then its SE counterpart.

3.2 Example 2. The Dirac’s free fermion complex non-hermitian matrix

We now proceed further in assessing the performance of the CC algorithm by including complex
matrix components as well as by strengthen the number of non-zero elements (in the millions
range). The tolerance now set to 10−5 with respect to the absolute value rather than for real and

8
 Table 2: Results of CC and SE on the two Wu and Schaeffer’s matrices
CC method SE method
Case 1 Case 2 Case 1 Case 2
N 103 104 Rounds per system 34.63 38.04
M −N 39200 23800 Total rounds 285732 202448
Eff. length size 15788 8163 Number of systems 8251 5322
E [tr(z′′ Qw)] 10371 20738 E [tr(z′′ Qw)] 10371 20738
V ar [tr(z′′ Qw)] 3932 8123 V ar [tr(z′′ Qw)] 2391 5301
MC St. error 0.499 0.998 MC St. error 0.499 0.998
Empirical St. error 0.502 1.047
CPU time per burn-in cycle 1.97 3.99
CPU time per cycle 1 2.03 CPU time per round 1.03 1.99
CPU time per eff. cycle 2.48 5.92 CPU time per system 35.69 75.80
Total CPU time 39403 48729 Total CPU time 294303 402872

imaginary part independently. The chosen example is very well representative of the computational
demanding tasks typical of the physical sciences. It belongs to the realm of elementary particle
physics known as Lattice Quantum Chromodynamics briefly described in the following.
In elementary particle physics the evolution of the particles is described in terms of a field,
ψ (x), over the space-time (x, is a four components vector, one representing time, the other three
the spatial coordinates) which governs the annihilation and creation of particles. More specifically
ψ (x)is an operator whose action on the elements of its domain (a Hilbert space) represents the
annihilation of a particle at time x0 and spatial coordinates (x1 , x2 , x3 ) . Conversely, the action
of the adjoin field ψ † (x) (the hermitian conjugate in a matrix representation of ψ (x)), describes
the creation of a particle at space-time coordinates x . For a class of particles called 1/2-spin
Fermions, such as the electron, ψis in turn a four components object (each component being an
operator) called a spinor, an object without a classical analog (it is not a vector since under a
2πrotation changes its sing). Free evolution of such particles, i.e., without interaction with other
particles, is governed by the Dirac’s equation which in the so called Euclidean representation reads

L ψ ≡ (∂µ γ µ + m)ψ = 0, (27)

where m is the mass of the particle, µ runs over the four spatio-temporal coordinates, ∂µ denotes
derivative along the µ coordinate, summation over repeated indexes is assumed and γ µ is a set of
four non-hermitian matrices of dimension four satisfying

γ µ γ ν + γ ν γ µ = 2δµν
(4)
,
(4)
where δµν denotes a four by four unit matrix multiplied by the standard Kronecker delta. Alter-
natively the evolution can be described through the associated Green function of Eq. (27), that
is, the solution of

L G(x, x′ ) = δ(x − x′ ), (28)

where δ(x − x′ )represents the Dirac’s delta. G(x, x′ )is known also as a propagator of the field
since given an initial configuration for the field ψ(x0 ), the field at any other space-time event is
obtained by Z
ψ(x) = G(x, x0 )ψ(x0 )dx0 ,

that is, G(x, x0 )represents the propagation of a Dirac’s delta signal from the event x0 to the event
x. In probabilistic terms, it gives the probability of finding a particle at the space-time event x
given that it was at the event x0 . Propagators are of paramount importance in quantum field
theory.

9
 In general the evolution is far more complicated than that described by (27) since interactions
among several fields will be effective. Thus, electrons, having electric charge will interact between
each other through the electromagnetic field (a gauge field in the jargon of quantum field theory)
or, in terms of particles, by exchange of photons. Eq. (27) must then be supplemented with both
a free term for the evolution of photons and an interaction term between photons and electrons.
Within this frame (Quantum Electrodynamics) the interaction is weak enough for a perturvative
approach in many interesting situations. However, there are other cases in which the coupling
is so strong that perturvative techniques are useless. The paradigmatic case is the interaction
of quarks (also fermions, like the electrons) inside nucleons (protons and neutrons), this time by
exchange of particles called gluons describing an interaction known as Quantum Chromodynamics
(QCD). In such a case, recourse is made to a discretization of the space-time (a space-time lattice)
suitable for a numerical solution of the problem. Under such discretization Eq. (28) leads to a
matrix equation of the form L G = I, so that the propagator is given by L−1 . Interactions make
L depend on the value of the gluon field at each spacetime event but still the discretized quark
propagator is given by a matrix inversion. However, an extra average over an ensemble of different
realizations of the fields is necessary. Altogether leads to an extremely demanding computational
task and Lattice QCD is responsible of a mayor part of the CPU time consumed in Science (for a
series of reviews in the computational aspects of Lattice QCD see [3] and in particular [13]). Here
we are only interested in testing the CC algorithm so that we restrict ourselves to the free case
described by (27 ). Although trivial from the physical point of view, this simple case can be solved
analytically [4] and therefore the solution can be checked. The final result after certain technical
details reads (in a compact notation) [4],
X (4) (4)

(4) µ µ
L m n = δm n + K (1 + γ ) δ mν +δµν , nν + (1 − γ ) δ mν −δµν , nν , (29)
µ

where Latin indexes run over lattice points, Greek indexes over space-time dimensions, K is a
constant approaching 1/8 in the continuum limit and δµν is the standard Kronecker delta. In this
notation a Latin index, say m, labels a lattice point while the ν spatio-temporal coordinate of
such a point is denoted by mν . Thus, in each lattice point a four dimensional matrix acting on
the spinor components is defined.
Explicitly in terms of both the spatio-temporal coordinates (Latin indexes) and spinor com-
ponents (Greek indexes) L reads

Lµ ν a b c d m n l k = δµν δa m δb n δc l δd k +
  1

1


K δb n δc l δd k δa+1 m δµν + γµν + δa−1 m δµν − γµν + ··· . (30)

In the next term in the sum γ 1 changes to γ 2 , a interchanges with band m with n and so for until
the indexes are exhausted. The γ’s matrices are defined as
 
0 σi
γi = , i = 1, 2, 3
σi 0

and  
4 0 1
γ =
1 0
where the σ’s, the 1 and the 0 must be understood as two by two matrices. The σ’s are the show
called Pauli matrices given by
 
0 1
σ1 = ,
1 0

 
0 −i
σ2 = ,
−i 0

10
  
3 1 0
σ = ,
0 −1

with i representing the imaginary unit. Eq. (30) is useful to implement L in a program but to
apply the CC algorithm as described in in section 2.4, a one to one mapping of the ten indexes of
L to only two is also needed. We use the following one usual in QCD

m′ = 1 + a + N1 (b + N2 (c + N3 (d + N0 µ)))
′
n = 1 + m + N1 (n + N2 (l + N3 (k + N0 ν))) ,

Here N1 , N2 , N3 and N0 are the number of discrete points in each spatio-temporal coordinate of
the lattice (as above, 0 corresponds to the time coordinate). The rank of the corresponding matrix
is then 4 N0 N1 N2 N3 while expression (30) yields fourteen non zero elements per row so that the
total number of non zero elements is 56 N0 N1 N2 N3 .
Again two cases were tested. The parameters are given in table 3 and corresponds to two
different spatio-temporal lattice sizes in both cases with the same number of discrete values for
the four coordinates. The value of K was arbitrarily chosen to 0.1.

Table 3: Topics of the Dirac’s free fermion matrix in two cases with different lattice sizes
Case 1 2
N1 , N2 , N3 , N0 18 20
Rank of C 419904 640000
Nonzero elements 5878656 8960000
K 0.1 0.1

Table 4: Results of CC and SE on Dirac’s matrices (first case).

CC method SE method
N 18 Rounds per system 55.03
M −N 10832 Total rounds 219167
Eff.
 length size
 10805.0 Number of systems 3983
E tr(z† Qw) exact value 413007.84+0i  
† †
E tr(z
 Qw)  413004.47-1.87i E tr(z
 Qw)  413005.08-1.98i
V ar tr(z† Qw) 184130.43 V ar tr(z† Qw) 67827.29
MC St. error 4.128 MC St. error 4.130
CPU time per burn-in cycle 1.43
CPU time per cycle 0.4746 CPU time per round 0.40
CPU time per eff. cycle 1.00 CPU time per system 21.89
Total CPU time 10859 Total CPU time 87167

11
 Table 5: Results of CC and SE on Dirac’s matrices (second case).
CC method SE method
N 18 Rounds per system 56.08
M −N 6782 Total rounds 146379
Eff.
 length size
 6848.8 Number of systems 2610
†
E tr(z Qw) exact value 629489.14+0i  
† †
E tr(z
 Qw)  629480.89-0.53i E tr(z
 Qw)  629482.78-0.1i
V ar tr(z† Qw) 270373.73 V ar tr(z† Qw) 103395.82
MC St. error 6.283 MC St. error 6.294
CPU time per burn-in cycle 2.26
CPU time per cycle 1.54 CPU time per round 0.61
CPU time per eff. cycle 1.53 CPU time per system 33.99
Total CPU time 10503 Total CPU time 88721

From inspection of tables 4 and 5 it is clear the advantage of the CC algorithm whose perfor-
mance again is around eight times higher than that of the SE. As explained above, this time it is
possible to check the results against exact deterministic calculations. As expected they coincide
with the estimated values within the imposed tolerance.

4 Discussion and Conclusion

Once the far superior efficiency of the CC algorithm has been demonstrated some comments about
its applicability seems in order but first a note about the GS . Notice that the GS as implemented
in (2) is just the particular case of CC corresponding to hermitian matrices since in this case
z† = w. Therefore its convergence is determined by sp(T) < 1 and the no singularity of D. These
conditions do not implies C being positive definite. For instance, in the trivial case of C being
(k) (k) √
diagonal, sp(T) = 0, and Eq. (2) reduces to zi = φi / cii which, although implies imaginary
(k)
values of zi works perfectly.
In the general case when the convergence criteria are not satisfied at least two alternatives ex-
its. One is to rewrite the algorithm with a different partition than that shown here. For instance,
suppose that D is singular. Then, a possibility is to use a partition with relatively small non singu-
lar diagonal blocks surrounding the zeros of D, amenable of been inverted deterministically with
the memory recourses available. In fact, in terms of time efficiency, there will be partitions more
efficient than the one used here but, obviously, without the appealing implementation simplicity
of Eqs. (3) and (4). As a rule of the thumb, the convergence will be guaranteed for a sufficiently
“heavy” D. A new partition in the way just described could solve problems with small elements
in the original D not only with zeros.
Another simpler possibility is just to reorder the rows and the columns of C trying to locate
in the diagonal large enough elements. In general, reordering which implies a low time penalty,
could be advantageous to improve the convergence of the method.
In summary, we have presented an efficient stochastic algorithm based in correlated Markov
chains suitable for the inversion of very large matrix whenever the memory recourses are not
enough for the application of the standard deterministic methods. The efficiency of the algorithm
has been tested in a couple of numerical examples rendering a dramatic improving of eight times
faster runs with respect to the best stochastic method known by the authors. The necessary and
sufficient conditions for the convergence of the algorithm have been also given.

References
[1] G. Arfen. Mathematical methods for physicists, chapter 4, page 281. Academic Press, London,
3th edition, 1985.

12
 [2] H. Baumgartel. Analytic Perturbation Theory for Matrices and Operators, volume 15 of
Operator Theory: Advances and Applications, page 41. Birkhauser Verlag, Basel, 1th edition,
1985.

[3] N. Cabibbo, Y. Iwasaki, and K. schilling (Guess Editors). High performance computing in
lattice QCD. Parallel Computing, 25:1197–1407, 1999.
[4] M. Creutz. Quarks, gluons and lattices, chapter 5, page 21. Cambridge University Press,
Cambridge, 1th edition, 1983.
[5] I. T. Dimov, T. T. Gurov, and T. V. Gurov. A new iterative Monte Carlo approach for inverse
matrix problem. Journal of Computational and Applied Mathematics, 92:15–35, 1998.

[6] S. Dong and K. Liu. Stochastic estimation with Z 2 noise. Physics Letters B, 328:130–136,
1994.
[7] L. A. García-Cortés. Multiple trait estimation of variance components in animal models with
different design matrices. In 6th World Congress on Genetics Applied to Livestock Production,
volume 18, pages 370–373. University of Guelph, Ontario, Canada, 1994.
[8] L. A. García-Cortés and D. Sorensen. Cálculo de la inversa de la matriz de coeficientes de las
ecuaciones de modelo mixto utilizando muestreo de gibbs y aplicaciones. ITEA, Vol. Extra
16:263–265, 1995.
[9] A. E. Gelfand and A. F. M. Smith. Sampling-based approaches to calculating marginal
densities. J. Amer. Statist. Assoc., 85:398–409, 1990.
[10] S. Geman and D. Geman. Stochastic relaxation, Gibbs distributions and the bayesian restora-
tion of images. IEEE Trans. Pattn. Anal. Mach. Iintell., 6:721–741, 1984.
[11] A. George and J. Liu. Computer solution of large sparse positive definite systems. Prencice
Hall, Englewood Cliffs, NJ, 1981.
[12] C. J. Geyer. Practical markov chain Monte Carlo. Statistical Science, 4:473–483, 1992.
[13] S. Gusken. Stochastic estimator techniques and their implementation on distributed parallel
computers. Parallel Computing, 25:1371–1381, 1999.

[14] D. A. Harville. Use of the gibbs sampler to invert large, possibly sparse, positive definite
matrices. Linear Albegra and its Applications, 289:203–224, 1999.
[15] C. R. Henderson. Applications of Linear Models in Animal Breeding. University of Guelph,
1984.
[16] V. E. Johnson. Studying convergence of markov chain monte carlo algorithms using coupled
sample paths. J. Am. Stat. Assoc., 91:154–166, 1996.

[17] J. B. Kogut. The lattice gauge-theory approach to quantum chromodynamics. Reviews of

Modern Physics, 55:775–836, 1983.
[18] R. Thompson. Integrating best linear unbiased prediction and maximum likelihood estima-
tion. In 5th World Congr. Genet. Appl. Livest. Prod., volume 18, pages 337–340, 1995.
[19] L. Wu and L. R. Schaeffer. Reducing the effect of parent averages from animal solutions in
mixed model equations. Journal of Animal Breeding and Genetics, 117:361–374, 2000.

[20] D. M. Young. Iterative solutions of large linear systems. Academic Press, Orlando, 1971.

13