Approximation of the Determinant of Large Sparse Symmetric Positive Definite Matrices.
This paper is concerned with the  problem of approximating det(A)1/n for a large sparse symmetric positive definite matrix A of order n. It is shown that an efficient solution of this problem is obtained  by using a sparse approximate inverse of A. The method is explained and theoretical properties are discussed. The method is ideal for implementation on a parallel computer. Numerical experiments are described  that illustrate the performance of this new method and provide a comparison with Monte Carlo--type  methods from the literature.
Introduction
. Throughout this paper, A denotes a real symmetric positive
definite matrix of order n with eigenvalues
In a number of applications, for example in lattice Quantum Chromodynamics [12],
certain functions of the determinant of A, such as det(A) 1=n or ln(det(A)) are of
interest. It is well-known (cf. also x2) that for large n the function A ! det(A) has
poor scaling properties and can be very ill-conditioned for certain matrices A. In this
paper we consider the function
A few basic properties of this function are discussed in x2. In this paper we present
a new method for approximating d(A) for large sparse matrices A. The method is
based on replacing A by a matrix which is in a certain sense close to A \Gamma1 and for
which the determinant can be computed with low computational costs. One popular
method for approximating A is based on the construction of an incomplete Cholesky
factorization. This incomplete factorization is often used as a preconditioner when
solving linear systems with matrix A. In this paper we use another preconditioning
technique, namely that of sparse approximate inverses (cf. [1, 7, 9, 11]). In Re-mark
3.10 we comment on the advantages of the use of sparse approximate inverse
preconditoning for approximating d(A). Let A = LL T be the Cholesky decomposition
of A. Then using techniques known from the literature a sparse approximate
inverse GE of L, i.e. a lower triangular matrix GE which has a prescribed sparsity
structure E and which is an approximation of L \Gamma1 , can be constructed. We then
use det(GE
ii as an approximation for d(A). In x3 we explain
the construction of GE and discuss theoretical properties of this sparse approximate
inverse. For example, such a sparse approximate inverse can be shown to exist for
any symmetric positive definite A and has an interesting optimality property related
to d(A). As a direct consequence of this optimality property one obtains that
holds and that the approximation of d(A) by det(GE ) \Gamma2=n becomes
better if a larger sparsity pattern E is used.
Institut fur Geometrie und Praktische Mathematik, RWTH Aachen, Templergraben 55, D-52056
Germany.
A. REUSKEN
In x4 we consider the topic of error estimation. In the paper [2] bounds for the determinant
of symmetric positive definite matrices are derived. These bounds, in which
the Frobenius norm and an estimate of the extreme eigenvalues of the matrix involved
are used, often yield rather poor estimates of the determinant (cf. experiments in [2]).
In x4.1 we apply this technique to the preconditioned matrix GEAG T
E and thus obtain
reliable but rather pessimistic error bounds. It turns out that this error estimation
technique is rather costly. In x4.2 we introduce a simple and cheap Monte Carlo technique
for error estimation. In x5 we apply the new method to a few examples of large
sparse symmetric positive definite matrices.
2. Preliminaries. In this section we discuss a few elementary properties of the
function d. We give a comparision between the conditioning of the function d and
of the fuction A ! We use the notation k \Delta k 2 for the Euclidean
norm and denotes the spectral condition number of A. The trace of
the matrix A is denoted by tr(A).
Lemma 2.1. Let A and ffiA be symmetric positive definite matrices of order n.
The following inequalities hold:
Proof. The result in (2.1a) follows from
Y
The result in (2.1b) follows from the inequality between the geometric and arithmetic
mean:
Y
From the Courant-Fischer characterization of eigenvalues it follows that
for all i. Hence holds. Now note that
Y
Y
APPROXIMATION OF DETERMINANTS 3
Thus the result in (2.1c) is proved.
The result in (2.1c) shows that the function d(A) is well-conditioned for matrices
A which have a not too large condition number (A).
We now briefly discuss the difference in conditioning between the functions A !
det(A). For any symmetric positive definite matrix B of order n we
have
From the Courant-Fischer eigenvalue characterization we obtain
for all i. Hence
B is SPD
B is SPD
with equality for I . Thus for the condition number of the function d we have
=n
Note that for the diagonal matrix
in the inequality in (2.2) one obtains equality for n !1. For this A and with
we have equality in the second inequality in (2.1c), too.
For ~
the condition number is given by
~
i.e. n times larger than the condition number in (2.2). The condition numbers for
d and ~
d give an indication of the sensitivity if the perturbation kffiAk 2 is sufficiently
small. Note that the bound in (2.1c) is valid for arbitrary symmetric positive definite
perturbations ffiA. The bound shows that even for larger perturbations the function
is well-conditioned at A if (A) is not too large. For the function ~
the
effect of relatively large perturbations can be much worse than for the asymptotic
case (ffiA ! 0), which is characterized by the condition number in (2.3). Consider, for
example, for
2 a perturbation
~
~
which is very large if, for example,
The results in this section show that the numerical approximation of the function
is a much easier task than the numerical approximation of A ! det(A).
3. Sparse approximate inverse. In this section we explain and analyze the
construction of a sparse approximate inverse of the matrix A. Let A = LL T be the
Cholesky factorization of A, i.e. L is lower triangular and L \Gamma1 AL I . Note that
ii . We will construct a sparse lower triangular approximation
G of L \Gamma1 and approximate d(A) by d(G)
ii . The construction of a
sparse approximate inverse that we use in this paper was introduced in [9, 10, 11] and
can also be found in [1]. Some of the results derived in this section are presented in
[1], too.
4 A. REUSKEN
We first introduce some notation. Let E ae f(i; ng be a given
sparsity pattern. By #E we denote the number of elements in E. Let SE be the set
of n \Theta n matrices for which all entries are set to zero if the corresponding index is not
in E:
For use the
representation
For we define the projection
Note that the matrix
is symmetric positive definite. Typical choices of the sparsity pattern E (cf. x5) are
such that n i is a very small number compared to n (e.g. In such a case the
projected matrix P i AP T
i has a small dimension.
To facilitate the analysis below, we first discuss the construction of an approximate
sparse inverse ME 2 SE in a general framework. For we use the
representation
Note that if n
For given A; B 2 R n\Thetan with A symmetric positive definite we consider the following
problem:
determine
In (3.3) we have #E equations to determine #E entries in ME . We first give two basic
lemmas which will play an important role in the analysis of the sparse approximate
inverse that will be defined in (3.9).
Lemma 3.1. The problem (3.3) has a unique solution
the ith row of ME is given by m T
i with
i is the ith row of B.
Proof. The equations in (3.3) can be represented as
(b T
for all i with
i is the ith row of ME . Consider an i with
. For the unknown entries in m i we obtain the system of equations
APPROXIMATION OF DETERMINANTS 5
which is equivalent to
The matrix P i AP T
i is symmetric positive definite and thus m i must satisfy
Using
we obtain the result in (3.4). The construction in this proof
shows that the solution is unique.
Below we use the Frobenius norm, denoted by k
Lemma 3.2. Let A = LL T be the Cholesky factorization of A and let ME 2 SE
be the unique solution of (3.3). Then ME is the unique minimizer of the functional
Proof. Let e i be the ith basis vector in R n . Take M 2 SE . The ith rows of M
and B are denoted by m T
The minimum of the functional (3.6) is obtained if in (3.7) we minimize the functionals
for all i with
(3.8) can be rewritten as
The unique minimum of this functional is obtained for "
Using Lemma 3.1 it follows that ME
is the unique minimizer of the functional (3.6).
Sparse approximate inverse. We now introduce the sparse approximate inverse
that will be used as an approximation for L \Gamma1 . For this we choose a lower triangular
ng and we assume that (i; i) 2 E l for all i. The
sparse approximate inverse is constructed in two steps:
l such that ( "
2:
6 A. REUSKEN
The construction of G E l in (3.9) was first introduced in [9]. A theoretical background
for this factorized sparse inverse is given in [11]. The approximate inverse "
(3.9a) is of the form (3.3) with I . From Lemma 3.1 it follows that in (3.9a) there
is a unique solution "
. Note that because E l is lower triangular and (i; i) 2 E l we
have
Hence it follows from Lemma 3.1
that the ith row of "
denoted by g T
i , is given by
The ith entry of g i , i.e. e T
which is strictly positive
because
i is symmetric positive definite. Hence diag( "
contains only strictly
positive entries and the second step (3.9b) is well-defined.
. The sparse
in (3.9a) can be computed by solving the (low dimensional)
symmetric positive definite systems
We now derive some interesting properties of the sparse approximate inverse as in
(3.9). We start with a minimization property of "
Theorem 3.3. Let A = LL T be the Cholesky factorization of A and D :=
l as in (3.9a) is the unique minimizer of the functional
Proof. The construction of "
in (3.9a) is as in (3.3) with I . Hence
Lemma 3.2 is applicable with I . It follows that "
l is the unique minimizer of
Decompose L \GammaT as L \GammaT strictly upper triangular. Then D
and R are lower and strictly upper triangular, respectively, and we obtain:
F
Hence the minimizers in (3.13) and (3.12) are the same.
Remark 3.4. From the result in Theorem 3.3 we see that in a scaled Frobenius
norm (scaling with D
l is the optimal approximation of "
in the set S E l , in
the sense that "
L is closest to the identity. A seemingly more natural minimization
problem is
min
i.e. we directly approximate L \Gamma1 (instead of "
do not use the scaling with
. The minimization problem (3.14) is of the form as in Lemma 3.2 with
APPROXIMATION OF DETERMINANTS 7
. Hence the unique minimizer in (3.14), denoted by ~
must satisfy (3.3)
Because E l contains only indices (i;
~
must satisfy
This is similar to the system of equations in (3.9a), which characterizes "
l . However,
in (3.16) one needs the values L ii , which in general are not available. Hence opposite
to the minimization problem related to the functional (3.12) the minimization problem
(3.14) is in general not solvable with acceptable computational costs. 2
The following lemma will be used in the proof of Theorem 3.7.
Lemma 3.5. Let "
l be as in (3.9a). Decompose "
L), with
D diagonal and "
strictly lower triangular. Define E l
ng.
L is the unique minimizer of the functional
and also of the functional
Furthermore, for "
D we have
Proof. From the construction in (3.9a) it follows that
is such that ( "
. This is of the form (3.3)
From Lemma 3.2 we obtain that "
L is the unique minimizer of
the functional
i.e., of the functional (3.17). From the proof of Lemma 3.2, with
that the minimization problem
min
decouples into seperate minimization problems (cf. (3.8)) for the rows of L:
min
l
Al i g (3.20)
8 A. REUSKEN
for all i with
i and a T
i are the ith rows of L and A, respectively. The
minimization problem corresponding to (3.18) is
min
Y
Y
Al
This decouples into the same minimization problems as in (3.20). Hence the functionals
in (3.17) and (3.18) have the same minimizer.
Using the construction of "
in (3.9a) we
obtain
ii J
D)
(i;k)2E l
Hence "
ii holds for all i, i.e., (3.19) holds.
Corollary 3.6. From (3.19) it follows that diag( "
and thus, using (3.9b) we obtain
for the sparse approximate inverse G E l . 2
The following theorem gives a main result in the theory of approximate inverses.
It was first derived in [11]. A proof can be found in [1], too.
Theorem 3.7. Let G E l be the approximate inverse in (3.9). Then G E l is the
unique minimizer of the functional
Proof. For G 2 S E l we use the decomposition diagonal
and L
. Furthermore, for L
The inequality in (3.23) follows from the inequality between the arithmetic and geometric
in (3.9a) we use the decomposition "
L). For the approximate
L). From Lemma 3.5
it follows that det(J L )  det(J "
. Furthermore from Lemma 3.5
we obtain that for G
L) we have
I and thus equality in
We conclude that G E l is the unique minimizer of the functional
APPROXIMATION OF DETERMINANTS 9
in (3.22).
Remark 3.8. The quantity
can be seen as a nonstandard condition number (cf. [1, 9]). Properties of this quantity
are given in [1] (Theorem 13.5). One elementary property is
Corollary 3.9. For the approximate inverse G E l as in (3.9) we have (cf.(3.21))
i.e.,
Y
Y
Let ~
l be a lower triangular sparsity pattern that is larger than E l , i.e., E l ae ~
ng. From the optimality result in Theorem 3.7 it follows that
~
Motivated
by the theoretical results in Corollary 3.9 we propose to use the sparse approximate
l as in (3.9) for approximating
an estimate for d(A). Some properties of this method are discussed in the following
remark.
Remark 3.10. We consider the method of approximating d(A) by d(G
. The practical realization of this method boils down to chosing a sparsity
and solving the (small) systems in (3.11). We list a few properties of this
1. The sparse approximate inverse exists for every symmetric positive definite
A. Note that such an existence result does not hold for the incomplete Cholesky
factorization. Furthermore, this factorization is obtained by solving low dimensional
symmetric positive definite systems of the form P i AP T
(cf. (3.11)), which can
be realized in a stable way.
2. The systems P i AP T
solved in parallel.
3. For the computation of d(G
need the diagonal
entries of "
(cf. (3.24)). In the systems P i AP T
e i we then only have to
compute the last entry of "
. If these systems are solved using the Cholesky
lower triangular) we only need the
of L i , since ("g i
4. The sparse approximate inverse has an optimality property related to the
determinant: The functional G From
this the inequality (3.24) and the monotonicity result (3.25) follow.
A. REUSKEN
5. From(3.24) we obtain the upper bound 0 for the relative error
1. In x4 we will derive useful lower bounds for this relative error. These are a posteriori
error bounds which use the matrix G E l . 2
4. A posteriori error estimation. In the previous section it has been explained
how an estimate d(G E l ) \Gamma2 of d(A) can be computed. From (3.24) we have
the error bound
In this section we will discuss a posteriori estimators for the error d(A)=d(G
In x4.1 we apply the analysis from [2] to derive an a posteriori lower bound for the
quantity in (4.1). This approach results in safe, but often rather pessimistic bounds
for the error. In x4.2 we propose a very simple stochastic method for error estimation.
This method, although it does not yield guaranteed bounds for the error, turns out
to be very useful in practice.
4.1. estimation based on bounds from [2] . In this section we show
how the analysis from [2] can be used to obtain an error estimator. We first recall a
main result from [2] (Theorem 2). Let A be a symmetric positive matrix of order n,
F and
exp
l
exp
In [2] this result is applied to obtain computable bounds for d(A). Often these
bounds yield rather poor estimates of d(A). In the present paper we approximate
use the result in (4.2) for error estimation. The upper bound
turns out to be satisfactory in numerical experiments (cf. x5). Therefore we
restrict ourselves to the derivation of a lower bound for d(A)=d(G E l ) \Gamma2 , based on the
left inequality in (4.2).
Theorem 4.1. Let G E l be the approximate inverse from (3.9) and 0 ! ff
1. The following holds: ff  1,
and
exp
Proof. The right inequality in (4.3) is already given in (4.1). We introduce the
for the eigenvalues of G E l AG T
l . From (3.21) we obtainn
and from this it follows that ff   1  1 holds. Furthermore,
APPROXIMATION OF DETERMINANTS 11
yields
We now use the left inequality in
(4.2) applied to the matrix G E l AG T
l . Note that
A simple computation yieldsn
l
and
Substitution of (4.5) in (4.4) results inn
l
Using this the left inequality in (4.3) follows from the left inequality in (4.2).
Note that for the lower bound in (4.3) to be computable, we need
F
and a strictly positive lower bound ff for the smallest eigenvalue of G E l AG T
l . We now
discuss methods for computing ff and . These methods are used in the numerical
experiments in x5.
We first discuss two methods for computing ff. The first method, which can be
applied if A is an M-matrix, is based on the following lemma, where we use the
Lemma 4.2. Let A be a symmetric positive definite matrix of order n with A ij  0
for all i 6= j and G E l its sparse approximate inverse (3.9). Furthermore, let z be such
that
Then
holds.
Proof. From the assumptions it follows that A is an M-matrix. In [11] (Theorem
4.1) it is proved that then G E l AG T
l is an M-matrix, too. Let z
Because
nonnegative entries it follows that
A. REUSKEN
Hence
Using  min (G E l AG T
we obtain the result (4.6).
Based on this lemma we obtain the following method for computing ff. Choose
apply the conjugate gradient method to the system G E l AG T
This results in approximations z of z   . One iterates until the stopping criterion
1 . In
view of efficiency one should not take a very small tolerance j. In our experiments
in x5 we use 1. Note that the CG method is applied to a system
with the preconditioned matrix G E l AG T
l . In situations where the preconditioning is
effective one may expect that relatively few CG iterations are needed to compute z j
such that kG E l AG T
numerical experiments are
presented in x5.
As a second method for determining ff, which is applicable to any symmetric positive
definite A, we propose the Lanczos method for approximating eigenvalues applied to
the matrix G E l AG T
l . This method yields a decreasing sequence  (1)
l ) of approximations  (j)
1 of  min (G E l AG T
holds, then
can be used in Theorem 4.1. However, in practice it is usually
not known how to obtain reasonable values for " in (4.7). Therefore, in our experiments
we use a simple heuristic for error estimation (instead of a rigorous bound as
in (4.7)), based on the observed convergence behaviour of  (j)
It is known that for the Lanczos method the convergence to extreme eigenvalues is
relatively fast. Moreover, it often occurs that the small eigenvalues of the preconditioned
are well-separated from the rest of the spectrum, which
has a positive effect on the convergence speed  (j)
In numerical
experiments we indeed observe that often already after a few Lanczos iterations we
have an approximation of  min (G E l AG T
with an estimated relative error of a few
percent. However, for the ff computed in this second method we do not have a rigorous
analysis which guarantees that
from numerical experiments we see that this method is satisfactory. This is partly explained
by the relatively fast convergence of the Lanczos method towards the smallest
eigenvalue. A further explanation follows from the form of the lower bound in (4.3).
For ff  which is typically the case in our experiments in x5, this lower
bound essentially behaves like exp(ffi ln ff) =: g(ff). Note that
holds. Hence the sensitivity of the lower bound with respect to perturbations in ff is
very mild.
We now discuss the computation of the quantity
F , which is needed
in (4.3). Clearly, for computing  one needs the matrices G E l and A. To avoid un-necessary
storage requirements one should not compute the matrix X := G E l AG T
and then determine
F . A with respect to storage more efficient approach
can be based on
where e i is the ith basis vector in R n . For the computation of kG E l AG T
APPROXIMATION OF DETERMINANTS 13
be done in parallel, one needs only sparse matrix-vector
multiplications with the matrices G E l and A. Furthermore, for the computation
of AG T
use that (DG E l
It follows from (3.10) that
holds.
Remark 4.3. Note that for the error estimators discussed in this section the
must be available (and thus stored), whereas for the computation of
the approximation d(G E l ) \Gamma2 of d(A) we do not have to store the matrix G E l (cf.
Remark 3.10 item 3). Furthermore, as we will see in x5, the computation of these
error estimators is relatively expensive. 2
4.2. estimation based on a Monte Carlo approach. In this section
we discuss a simple error estimation method which turns out to be useful in practice.
Opposite to those treated in the previous section this method does not yield (an
approximation of) bounds for the error.
The exact error is given by
l is a sparse symmetric positive definite matrix. Fore ease of
presentation we assume that the pattern E l is sufficiently large such that
holds. In [11] it is proved that if A is an M-matrix or a (block) H-matrix then (4.8)
is satisfied for every lower triangular pattern E l . In the numerical experiments (cf.
with matrices which are not M-matrices or (block) H-matrices (4.8) turns out to
be satisfied for standard choices of E l . We note that if (4.8) does not hold then the
technique discussed below can still be applied if one replaces
a suitable damping factor such that ae(I \Gamma !E
For the exact error we obtain, using a Taylor expansion of ln(I \Gamma B) for B 2 R n\Thetan
with
\Gamman
Hence, an error estimation can be based on estimates for the partial sums Sm :=
The construction of G E l is such that diag(E E l
and thus tr(E
14 A. REUSKEN
For S 3 we obtain
Note that in S 2 and S 3 the quantity tr(E 2
F occurs which
is also used in the error estimator in x4.1. In this section we use a Monte Carlo
method to approximate the trace quantities in Sm . The method we use is based on
the following proposition [8, 3].
Proposition 4.4. Let H be a symmetric matrix of order n with tr(H) 6= 0.
Let V be the discrete random variable which takes the values 1 and \Gamma1 each with
probability 0:5 and let z be a vector of n independent samples from V . Then z T Hz is
an unbiased estimator of tr(H):
E(z
and
For approximating the trace quantity in S 2 we use the following Monte Carlo algorithm

1. Generate z j 2 R n with entries which are uniformly distributed in (0; 1).
2. If (z j
3. y j
Based on Proposition 4.4 and (4.10) we use
as an approximation for S 2 . The corresponding error estimator is
For the approximation of S 3 we replace step 3 in the algorithm above by
3. y j
and we use
as an estimate for S 3 . The corresponding error estimator is
exp(\Gamman
APPROXIMATION OF DETERMINANTS 15
Clearly, this technique can be extended to the partial sums Sm with m ? 3. However,
in our applications we only use "
S 3 for error estimation. It turns out that, at
least in our experiments, the two leading terms in the expansion (4.9) are sufficient
for a reasonable error estimation. Note that due to the truncation of the Taylor ex-
pansion, the estimators E 2 and E 3 for are biased.
It is shown in [3] that based on the so-called Hoeffding inequality (cf. [13]) probabilistic
bounds for
can be derived, where z are
independent random variables as in Proposition 4.4. In this paper we do not use
these bounds. Based on numerical experiments we take a fixed small value for the
parameter M in the Monte Carlo algorithm above (in the experiments in x5:
Remark 4.5. In the setting of this paper Proposition 4.4 is applied with
is a known polynomial of degree 2 or 3. In the Monte Carlo technique
for approximating Proposition 4.4 is applied with
ln(A). The quantity z T ln(A)z, which can be considered as a Riemann-Stieltjes
integral, is approximated using suitable quadrature rules. In [3] this quadrature is
based on a Gauss-Christoffel technique where the unknown nodes and weights in the
quadrature rule are determined using the Lanczos method. For a detailed explanation
of this method we refer to [3].
A further alternative that could be considered for error estimation is the use of this
method from [3]. In the setting here, this method could be used to compute a (rough)
approximation of det(G E l AG T
We did not investigate this possibility. The results
in [2, 3] give an indication that this alternative is probably much more expensive
than the method presented in this section. 2
5. Numerical experiments. In this section we present some results of numerical
experiments with the methods introduced in x3 and x4. All experiments are done
using a MATLAB implementation. We use the MATLAB notation nnz(B) for the
number of nonzero entries in a matrix B.
Experiment 1 (discrete 2D Laplacian). We consider the standard 5-point discrete
Laplacian on a uniform square grid with m mesh points in both directions, i.e.
For this symmetric positive definite matrix the eigenvalues are known:
sin
For the choice of the sparsity pattern E l we use a simple approach based on the
nonzero structure of (powers of) the matrix A:
We first describe some features of the methods for the case
that we will vary m and k. Let A denote the discrete Laplacian for the case
LA its lower triangular part. We then have nnz(LA 2640. For the sparse approximate
inverse we obtain nnz(G E l (2) 6002. The systems P i AP T
that have to be solved to determine G E l (2) (cf. (3.11)) have dimensions between 1 and
7; the mean of these dimensions is 6.7. As an approximation of
obtain
Y
A. REUSKEN
Hence 0:965. For the computation of this approximation along
the lines as described in Remark 3.10, item 3, we have to compute the Cholesky
factorizations
n. For this approximately
are needed (in the MATLAB implementation). If we compare this with the costs of
one matrix-vector multiplication A   x (8760 flops), denoted by MATVEC, it follows
that for computing this approximation of d(A), with an error of 3.5 percent, we need
work comparable to only 5 MATVEC.
We will see that the arithmetic costs for error estimation are significantly higher.
We first consider the methods of x4.1. The arithmetic costs are measured in terms
of MATVEC. For the computation of ff as indicated in Lemma 4.2 with
using the CG method with starting vector need 8 iterations.
In each CG iteration we have to compute a matrix-vector multiplication G E l AG T
which costs approximately 3.7 MATVEC. We obtain ff 0:0155. For the method
based on the Lanczos method for approximating  min (G E l AG T
use the heuristic
stopping criterion
We then need 7 Lanczos iterations, resulting in ff direct computation
results in  min (G E l AG T
For the computation of
F we first computed the lower triangular part
of
l and then computed kXkF (making use of symmetry). The total
costs of this are approximately MATVEC. Application of Lemma 4.1, with ff CG
and ff Lanczos yields the two intervals
which both contain the exact error 0.965. In both cases, the total costs for error
estimation are 40-45 MATVEC, which is approximately 10 times more than the costs
for computing the approximation d(G E l (2) ) \Gamma2 .
We now consider the method of x4.2. We use the estimators E 2 and E 3 from
(4.13), (4.15) with 6. The results are 0:973. Note that
the order of magnitude of the exact error (3:5 percent) is approximated well by both
In step 3 in the Monte Carlo algorithm
for computing "
need one matrix-vector multiplication G E l AG T
MATVEC). The total arithmetic costs for E 2 are approximately 20 MATVEC. For
S 3 we need two matrix-vector multiplications with l in the third step of the Monte
Carlo algorithm. The total costs for E 3 are approximately 40 MATVEC.
In

Table

5.1 we give results for the discrete 2D Laplacian with
We use the sparsity pattern E l (2).
In the third column of this table we give the computed approximation of d(A) and the
corresponding relative error. In the fourth column we give the total arithmetic costs
for the Cholesky factorization of the matrices P i AP T
item 3). In the columns 5-8 we give the results and corresponding arithmetic costs
for the error estimators discussed in x4. The fifth column corresponds to the method
discussed in x4.1 with ff determined using the CG method applied to G E l AG T
with starting vector 1. In the stopping criterion we take
The computed used as input for the lower bound in (4.3). The resulting
bound for the relative error and the arithmetic costs for computing this error bound
are shown in column 5. In column 6 one finds the computed error bounds if ff is
determined using the Lanczos method with stopping criterion (5.2). In the last two
APPROXIMATION OF DETERMINANTS 17

Table
Results for 2D discrete Laplacian with
costs for Thm. 4.1, Thm. 4.1, MC MC

Table
Results for 2D discrete Laplacian with
costs for Thm. 4.1, Thm. 4.1, MC MC
columns the results for the Monte Carlo estimators are given.
In

Table

5.2 we show the results and corresponding arithmetic costs for the method
with sparsity pattern
Concerning the numerical results we note the following. From the third and fourth
column in Table 5.1 we see that using this method we can obtain an approximation
of d(A) with relative error only a few percent and arithmetic costs only a few
MATVEC. Moreover, this efficiency hardly depends on the dimension n. Comparison
of the third and fourth columns of the Tables 5.1 and 5.2 shows that the approximation
significantly improves if we enlarge the pattern from E l (2) to E l (4). The
corresponding arithmetic costs increase by a factor of about 9. This is caused by
the fact that the mean of the dimensions of the systems P i AP T
from approximately 7 (E l (2)) to approximately 20. For
For the other n
values we have similar ratios between the number of nonzeros in the matrices LA and
. Note that the matrix G E l has to be stored for the error estimation but not
for the computation of the approximation d(G E l ) \Gamma2 . The error bounds in the fifth
and sixth column in the Tables 5.1 and 5.2 are rather conservative and expensive.
Furthermore there is some deterioration in the quality and a quite strong increase in
the costs if the dimension n grows. The strong increase in the costs is mainly due to
the fact that the CG and Lanczos method both need significantly more iterations if n
increases. This is a well-known phenomenom (the matrix G E l AG T
E l has a condition
number that is proportional to n). Also note that the costs for these error estimators
are (very) high compared to the costs of the computation of d(G E l ) \Gamma2 . The results
in the last two columns indicate that the Monte Carlo error estimators, although less
reliable, are more favourable.
In

Figure

5.1 we show the eigenvalues of the matrix G E l AG T
l for the case
(computed with the MATLAB function eig). The eigenvalues are in the
A. REUSKEN
interval [0:025; 1:4]. The mean of these eigenvalues is 1
can see that relatively many eigenvalues are close to 1 and only a few eigenvalues are
close to zero.
Fig. 5.1. Eigenvalues of the matrix G E l AG T
l in Experiment 1
Experiment 2 (MATLAB random sparse matrix). The sparsity structure of the
matrices considered in Experiment 1 is very regular. In this experiment we consider
matrices with a pattern of nonzero entries that is very irregular. We used the
MATLAB generator (sprand(n; n; 2=n)) to generate a matrix B of order n with approximately
2n nonzero entries. These are uniformly distributed random entries in
(0; 1). The matrix B T B is then sparse symmetric positive semidefinite. In the generic
case this matrix has many eigenvalues zero. To obtain a positive definite matrix we
generated a random vector d with all entries chosen from a uniform distribution on
the interval (0; 1) (d :=rand(n; 1)). As a testmatrix we used A := B T B+diag(d). We
performed numerical experiments similar to those in Experiment 1 above. We only
consider the case with sparsity pattern (2). The error estimator based on the
CG method is not applicable because the sign condition in Lemma 4.2 is not fulfilled.
For the case 900 the eigenvalues of A and of G E l AG T
are shown in Figure 5.2.
For A the smallest and largest eigenvalues are 0:0099 and 5:70, respectively. The
picture on the right in Figure 5.2 shows that for this matrix A sparse approximate
inverse preconditioning results in a very well-conditioned matrix. Related to this, one
can see in Table 5.3 that for this random matrix A the approximation of d(A) based
on the sparse approximate inverse is much better than for the discrete Laplacian in
Experiment 1. For
and respectively. For the
mean of the dimensions of the systems P i AP T
spectively. In all three cases the costs for a matrix-vector multiplication G E l AG E l x
are approximately 4.3 MV. Furthermore, in all three cases the matrix G E l AG T
l is
well-conditioned and the number of Lanczos iterations needed to satisfy the stopping
criterion (5.2) hardly depends on n. Due to this, for increasing n, the growth in the
costs for the error estimator based on Theorem 4.1 (column 5) is much slower than in
Experiment 1. As in the Tables 5.1 and 5.2, in Table 5.3 the error quantities in the
columns 3, 5,6,7 are bounds or estimates for the relative error
APPROXIMATION OF DETERMINANTS 19
Fig. 5.2. Eigenvalues of the matrices A and G E l AG T
l in Experiment 2
Table
Results for MATLAB random sparse matrices with
costs for Thm. 4.1, MC MC
For the values of d(A) are not given (column 2). This has to do
with the fact that for these matrices with very irregular sparsity patterns the Cholesky
factorization suffers from much more fill-in than for the matrices in the Experiments
1 and 3. For the matrix A in this experiment with
10000 we run into storage problems
if we try to compute the Cholesky factorization using the MATLAB function chol.
Experiment 3 (QCD type matrix). In this experiment we consider a complex
Hermitean positive definite matrix with sparsity structure as in Experiment 1. This
matrix is motivated by applications from the QCD field. In QCD simulations the
determinant of the so-called Wilson fermion matrix is of interest. These matrices
and some of their properties are discussed in [4, 5]. The nonzero entries in a Wilson
fermion matrix are induced by a nearest neighbour coupling in a regular 4-dimensional
grid. These couplings consist of 12 \Theta 12 complex matrices M xy , which have a tensor
product structure M
xy\Omega U xy , where P xy 2 R 4\Theta4 is a projector, U xy 2 C 3\Theta3 is
from SU 3 and x and y denote nearest neighbours in the grid. These coupling matrices
M xy strongly fluctuate as a function of x and y. Here we consider a (toy) problem
with a matrix which has some similarities with these Wilson fermion matrices. We
start with a 2-dimensional regular grid as in Experiment 1 (n grid points). For the
couplings with nearest neighbours we use complex numbers with length 1. These
numbers are chosen as follows. The couplings with south and west neighbours at a
grid point x are exp(2iff S (x)) and exp(2iff W (x)), respectively, where ff S (x) and
ff W (x) are chosen from a uniform distribution on the interval (0; 1). The couplings
with the north and east neighbours are taken such that the matrix is hermitean. To
make the comparison with Experiment 1 easier the matrix is scaled by the factor n,
A. REUSKEN
i.e. the couplings with nearest neighbours have length n. For the diagonal we take flI ,
where fl is chosen such that the smallest eigenvalue of the resulting matrix is approximately
1 (this can be realized by using the MATLAB function eigs for estimating
the smallest eigenvalue). We performed numerical experiments as in Experiment 1
with (2). The number of nonzero entries in LA and G E l are the same as
in Experiment 1. For 900 the eigenvalues of the matrices A and G E l AG T
are
shown in Figure 5.3. These spectra are in the intervals
respectively.
The results of numerical experiments are presented in Table 5.4. Note that the error
Fig. 5.3. Eigenvalues of the matrices A and G E l AG T
l in Experiment 3
estimator from x4.1 in which the CG method is used for computing ff can not be used
for this matrix (assumptions in Lemma 4.2 are not satisfied). We did not consider the
case here because then the application of the eig function for computing
the smallest eigenvalue led to memory problems.
Comparison of the results in Table 5.4 with those in Table 5.1 shows that when the

Table
Results
costs for Thm. 4.1, MC MC
method is applied to the QCD type of problem instead of the discrete Laplacian the
performance of the method does not change very much.
Finally, we note that in all measurements of the arithmetic costs we did not take
into account the costs of determining the sparsity pattern E l (k) and of building the
matrices



--R

Cambridge University Press
Bounds on the trace of the inverse and the determinant of symmetric positive definite matrices
Some large scale matrix computation problems
Progress on lattice QCD algorithms
Exploiting structure in Krylov subspace methods for the Wilson fermion matrix
Matrix Computations
Parallel preconditioning with sparse approximate inverses
A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines
An alternative approach to estimating the convergence rate of the CG method
On a family of two-level preconditionings of the incomplete block factorization type
Factorized sparse approximate inverse precondi- tionings I : Theory
Quantum Fields on a Lattice
Convergence of Stochastic Processes
--TR
sparse approximate inverse;determinant;preconditioning