Computing Betti Numbers via Combinatorial Laplacians

Joel Friedman Department of Mathematics 1984 Mathematics Road University of British Columbia Vancouver, BC V6T 1Z2 CANADA

                  jf@math.ubc.ca

June 4, 1997

Abstract

We use the Laplacian and power method to compute Betti numbers of simplicial complexes. This has a number of advantages over other methods, both in theory and in practice. It requires small storage space in many cases. It seems to run quickly in practice, but its running time depends on a ratio, ν, of eigenvalues which we have yet to fully understand.

We numerically verify a conjecture of Björner, Lovász, Vrećica, and Z̆ivaljević on the chessboard complexes C(4, 6), C(5, 7), and C(5, 8). Our verification suffers a technical weakness, which can be overcome in various ways; we do so for C(4, 6) and C(5, 8), giving a completely rigourous (computer) proof of the conjecture in these two cases. This brings up an interesting question in recovering an integral basis from a real basis of vectors.

1 Introduction

A fundamental and important set of invariants of a topological space, X, is its collection of homology groups, H_{i}(X) = H_{i}(X, Z). Computing these groups is of great importance, especially in topology and pure mathematics. However, there are numerous applied areas where computing homology is of interest. These include pattern recognition and classification in biology, chemistry, robotics, and scene analysis, involving low dimensional topological spaces, as well as time series analysis and dynamical systems, involving higher dimensional spaces (see [Cha95]).

A part of the homology groups are the Betti numbers, the i-th Betti number, b_{i} = b_{i}(X), being the rank of H_{i}(X). The Betti numbers often have intuitive meanings. For example, b_{0} is simply the number of connected components of the space. As another example, the oriented two-dimensional manifolds are completely classified by b_{1} = 2g, where g is the genus (i.e. number of “handles”) of the two-manifold. Knowing the Betti numbers is the same as knowing the rational homology (i.e. homology with rational coefficients), which is the same as knowing the homology groups, H_{i}(X), up to torsion.

When one wants to know the homology groups, sometimes it suffices to know the Betti numbers. First, if one is trying to distinguish two objects via their homology, their Betti numbers may already distinguish them; similarly, to prove the nonvanishing of a homology group it suffices (but is not necessary) to have the corresponding Betti number vanish. Second, there are certain situations where one a priori knows that there is no torsion in the homology groups; then homology is completely described by the Betti numbers. Moreover, it is often much simpler to compute the Betti numbers, and so one does this as a first step in determining the homology groups.

A simplicial complex is a collection of subsets (called faces) of a fixed set, S, which is closed under taking subsets. It gives rise to a certain topological space, and the homology of any “reasonable” topological space can be computed from those of an “approximating” simplicial complex (see [Mun84]). In this paper we restrict ourselves to input data being a simplicial complex (actually it suffices to list the maximal faces). In applications one is often given a topological space directly in this form.

Computing homology at present seems hard; currently algorithms require computing the Smith normal form of the matrix. The latter can be done in polynomial time (see [KB79]), but the current polynomial required is still quite large in degree^{1}. Computing Betti numbers, on the other hand, seems much easier. To find b_{i}, one only need compute the ranks of two matrices. This gives what might be called the “standard algorithm,” requiring S = O(T n_{i}(n_{i−1} + n_{i+1})) storage and O(S max(n_{i}, n_{i−1} + n_{i+1}) log T) operations, where n_{i} is the number of i-faces (i.e. faces of dimension i, i.e. subsets of size i + 1), and where T is the number of storage registers needed to store a matrix entry during the rank calculation; unless the model allows storage registers of arbitrary size, T would need to be O(min(n_{i}, n_{i−1} + n_{i+1})). This points to a practical problem: when the n_{j} are large an exact computation of the rank may be unreliable unless a multiple precision package is used; thus the operations become very expensive. It becomes important to have techniques for computing the Betti numbers which are faster and require less storage.

As a step in this direction, [DE93] give an algorithm to compute Betti numbers; it works in almost linear time in the n_{i}, but only works for subtriangulations of R^{3}. Their method is based on a standard type of topological induction, and for more general simplicial complexes it is not clear how to make this method fast.

1In [DC91] it is claimed that this can be done faster for sparse matrices, but the latest version of this paper still has serious flaws.

In this paper we exploit the combinatorial Laplacian to give a simple algorithm to compute the Betti numbers of a simplicial complex. It gives an algorithm which given ϵ, δ > 0 computes b_{i} with probability ≥ 1 − δ, requiring S = O(n_{i}(d_{i} + b_{i})) storage and O(S(b_{i} + 1)r) + O(n_{i}(n_{0} − i)(i + 1) log N) operations; here N = max(n_{i−1}, n_{i}, n_{i+1}), d_{i} is the average degree of the i-th Laplacian (always ≤ 1 + (i + 1)(n_{0} − i − 1)), and r = O(log(n_{i}b_{i}δ^{−1}ϵ^{−1}ν^{−1})/ν), where ν = (3/4)λ_{1}/λ_{max} is a ratio of Laplacian eigenvalues discussed later; here we are assuming that ϵ is chosen so that 2ϵb_{i} < ν. We do not need to know ν or b_{i} before we choose ϵ, but we do need to choose ϵ so that we are confident that 2ϵb_{i} < ν holds; in practice this doesn’t seem to be a problem. It is an important question to determine what ν is, in various situations. In examples of interest here, 1/ν is always ≤ 25 while n_{i} gets as high as around 10, 000. We know that for i = 0, 1, ν can be as small as roughly 1/n_{i}^{2}; we also know that for i = 0 we have ν ≥ 1/2 with high probability in certain random settings.

After running the main algorithm, one can try to verify rigourously (not merely with high probability and not requiring assumptions on ν) that b_{i} is at least the computed value for b_{i} by an “integralizing” process. This requires no more storage and less time than does the main algorithm. However, this does not always work in single or double precision, and brings up an interesting question about the ability to “integralize.”

We use our algorithm to compute the Betti numbers of a number of complexes, especially the “chessboard complexes” (see [BLVZ̆94]); these complexes arise in several contexts, including some combinatorial geometric problems where there homology groups are of interest (see [Z̆V92, ABFK92]). We use our algorithm to verify a conjecture of [BLVZ̆94] involving three of the chessboard complexes. In two cases we are able to “integralize” and give rigourous demonstrations of the conjecture in these cases. We also notice that all eigenvalues of the Laplacians are integers for all the complexes we checked; we conjecture this is always true^{2}.

We also mention that with our algorithm it is much easier to verify (with high probability) that a Betti number vanishes, rather than to show that it doesn’t and to evaluate it. But indeed, for a large class of combinatorial complexes (e.g. those which are shellable) only one Betti number other than b_{0} = 1 does not vanish; for these complexes we can compute all the Betti numbers (using the known Euler characteristic) as soon as we can detect whether or not a Betti number vanishes. This method was used in [Fri92] to help successfully calculate the Betti numbers of a large family of complexes.

We finish by summarizing the rest of this paper. In section 2 we describe the combinatorial Laplacians by way of Hodge theory. In section 3 we analyze some algorithms for determining the largest eigenvalue of a symmetric matrix and the eigenvalue’s multiplicity. In section 4 we prove theorems about our Betti number computing algorithms. In section 5 we give two complexes whose Betti numbers we compute via our algorithm, and we give some computational data about our numerical experiments. In section 6 we describe our findings on the chessboard complexes. In section 7 we describe some provable and some conjectural lower bounds for λ_{1} and λ_{1}/λ_{max} of the i-th Laplacian.

2This has lead to the recent work [FH], where this and part of the [BLVZ̆94] conjecture is proven.

2 Hodge Theory and the Laplacian

Our main tool is to use the combinatorial Laplacians (see [Hod41, Eck45, Dod76, DP76]) to compute the Betti numbers. These Laplacians are most easily described via the Hodge theory of Hodge [Hod41].

Recall that the Betti numbers, b_{i}, are the dimensions of the homology groups, H_{i} = ker(∂_{i})/im(∂_{i+1}) of the chain complex,

⋯ ⟶ C_{i+1} ⟶^{∂_{i+1}} C_{i} ⟶^{∂_{i}} C_{i−1} ⟶ ⋯ ⟶ C_{−1} = 0,  (2.1)

where C_{i} is the space of formal R-linear sums of oriented i-dimensional faces, i.e. oriented subsets of the abstract simplicial complex of size i + 1, and ∂_{i} is the boundary map (see [Mun84]).

Concretely, for each face, A, we fix an ordering of its vertices, A. ∂_{i} is an n_{i−1} × n_{i} matrix indexed on the (i − 1)- and i-faces as follows: if A, B are i- and (i − 1)-faces respectively, with orderings A, B, then ∂_{i} at B, A is zero unless B ⊂ A; if B ⊂ A, then there is a unique permutation π on A such that B is π(A) minus its first element (in order), and ∂_{i} at B, A is the sign of π.

So the “standard algorithm” to compute b_{i} referred to in the introduction is simply to determine the ranks of ∂_{i} and ∂_{i−1}.

Hodge theory works for an arbitrary chain complex over R (or any field of characteristic 0, such as Q or C). Recall that a chain complex is a collection, C_{i}, of vector spaces, with maps ∂_{i}: C_{i} → C_{i−1}, as in equation 2.1, such that ∂_{i−1} ∘ ∂_{i} = 0 for all i. Endowing each C_{i} with an inner product, we get maps ∂_{i}^{∗}: C_{i−1} → C_{i} (i.e. the transpose of ∂_{i}), and thus a Laplacian (as in [Dod76]), Δ_{i}: C_{i} → C_{i}, for each i, defined by

Δ_{i} = ∂_{i+1}∂_{i+1}^{∗} + ∂_{i}^{∗}∂_{i}.

For each i we define the set of harmonic i-forms to be

H_{i} = {c ∈ C_{i}|Δ_{i}c = 0}.

For chain complexes where each C_{i} is a finite dimensional R-vector space, Hodge theory involves only elementary linear algebra, and says:

Proposition 2.1 (Hodge theory) For each i we have H_{i} ≅ H_{i}, in that each member of H_{i} gives rise to a class in H_{i}, and each class in H_{i} contains a unique harmonic form in H_{i}.

In more detail, for each i we have that C_{i} decomposes into orthogonal subspaces

C_{i} = H_{i} ⊕ im(∂_{i+1}) ⊕ im(∂_{i}^{∗}).

The Laplacian is positive definite on the latter two summands and is invariant on each. ∂_{i+1}∂_{i+1}^{∗} is invariant on the middle summand and vanishes on the other two, and similarly for ∂_{i}^{∗}∂_{i}.

Proof Follows easily from the facts that (1) A = ∂_{i}^{∗}∂_{i} and B = ∂_{i+1}∂_{i+1}^{∗} are positive semi-definite and commute, satisfying AB = BA = 0, and (2) imS = im(S ∘ S^{∗}) for any map of finite dimensional inner product spaces, S: V → W.

□

We may therefore calculate the Betti numbers as the dimension of the H_{i}. The algorithm we will give works better the larger the eigenvalues of Δ_{i} are on im(∂_{i+1}) ⊕ im(∂_{i}^{∗}); this space is just H_{i}^{⊥}, by Hodge theory. One easy observation which we will use is:

Proposition 2.2 The set of (non-zero) eigenvalues of Δ_{i} on H_{i}^{⊥} is a subset of the union of those of Δ_{i−1} on H_{i−1}^{⊥} and of Δ_{i+1} on H_{i+1}^{⊥}.

Proof This follows immediately from Hodge theory and the fact that for any matrix, A, we have that AA^{∗} and A^{∗}A have the same set of non-zero eigenvalues.

□

3 The Power Method

In this section we discuss the usual power method for finding the largest eigenvalue. We will give an analysis of certain aspects of it which we will need.

Let A be a symmetric matrix, with eigenvalues λ_{1} ≥ λ_{2} ≥ ⋯ ≥ λ_{n} ≥ 0 and an orthonormal set of eigenvectors, u_{i}, i = 1, …, n, with Au_{i} = λ_{i}u_{i}. (At the end of this section we weaken the restriction λ_{n} ≥ 0.) We wish to find λ_{1}, and perhaps its multiplicity and a basis for its eigenspace.

Our version of the power method finds λ_{1} as follows. We choose a random vector, v, by taking its coordinates to be iid normal with mean 0 and variance 1. We set v_{r} = A^{r}v, and consider

T (v_{r}) = (Av_{r}, Av_{r})/(v_{r}, v_{r}) = (v_{r+1}, v_{r+1})/(v_{r}, v_{r}),

where ( , ) denotes the inner product. T (v_{r}) should converge to λ_{1}^{2}, and accordingly v_{r} should approach (when suitably normalized) an eigenvector. We take r large enough (see below) so that T (v_{r}) is close to λ_{1}^{2}. If we only wish to find λ_{1}, we are done.

If we wish to determine the multiplicity of λ_{1}, we set u_{1} = v_{r}/∥v_{r}∥, and apply the same power method to Ã = P AP where P is the projection onto the orthogonal complement of v_{1}. The power method gives us a sequence v_{i}, each orthogonal to u_{1}, with T (v_{r}) approaching λ_{2}^{2}. Again, for large r we see if T (v_{r}) is close (see below) to our computed value of λ_{1}^{2}. If not we are done; if so we start a third sequence of iterations, on Ã = P AP with P projecting out u_{1}, u_{2}. Continuing in this way we will find orthonormal u_{1}, …, u_{m} with each T (u_{i}) close to λ_{1}^{2}, and the next set of iterations giving v_{i} with T (v_{r}) not close to λ_{1}^{2}; we conclude λ_{1} has multiplicity m.

We now analyze the algorithm.

Let p_{i} = (v, u_{i}). The p_{i} are therefore iid N(0, 1). The quantity

L = (p_{2}^{2} + ⋯ + p_{n}^{2})/p_{1}^{2}

will be crucial to analyzing the convergence of the power method; for example, if λ_{2} ≠ λ_{1} then L < ∞ iff the power method converges. It is not hard to estimate:

Lemma 3.1 For any α > 0, the probability that L ≥ nα is at most \sqrt{2/(πα)}.

Proof See [Fri95].

L’s importance can be seen from the next lemma.

Lemma 3.2 Let T (v_{r}) ≤ λ_{1}^{2}(1−η) for some r and η > 0. Then L ≥ η(1−η)^{−r}. Also, if λ_{1} has multiplicity m and λ_{m+1}/λ_{1} = μ, then L ≥ ημ^{−2r}.

Proof We will easily reduce this to the following lemma:

Lemma 3.3 If q_{2}, …, q_{n} are reals with ∑ q_{i}^{2} ≤ K, then for any λ_{2}, …, λ_{n} ∈ R with |λ_{i}| ≤ 1 for all i we have:

\frac{1 + ∑_{i>1} λ_{i}^{2r+2}q_{i}^{2}}{1 + ∑_{i>1} λ_{i}^{2r}q_{i}^{2}} ≤ 1 − η

implies that K(1 − η)^{r} ≥ η. Also, if λ_{i} ≤ μ whenever λ_{i} ≠ 1, then η ≤ Kμ^{2r}.

Proof Call the left-hand-side of the above equation f, viewed as a function of the λ_{i} and the q_{i}.

First we claim that f is minimized when all the λ_{i} are equal up to sign. Indeed, viewing the q_{i} as fixed, it is clear that f achieves its minimum, f_{min} somewhere with the λ_{i} ranging over the compact set [−1, 1]^{n}. By differentiating g(λ) = (A + q^{2}λ^{2r+2})/(B + q^{2}λ^{2r}) one sees that g is minimized when λ^{2} = (r/(r + 1))g(λ). So at the minimum we have λ_{i}^{2} = (r/(r + 1))f_{min} for each i, and so all λ_{i} are equal there up to sign.

So at its minimum value,

f = \frac{1 + Qλ^{2r+2}}{1 + Qλ^{2r}} = 1 − \frac{λ^{2r} − λ^{2r+2}}{(1/Q) + λ^{2r}}

where Q = ∑ q_{i}^{2}, which is clearly minimized when Q is as large as possible, namely K. So it remains to consider the minimum of:

f(λ) = \frac{1 + Kλ^{2r+2}}{1 + Kλ^{2r}} .

Now if f(λ) ≤ 1 − η then we have λ^{2} ≤ 1 − η, since

f(λ) = (\frac{1}{1 + Kλ^{2r}} ) + (\frac{Kλ^{2r}}{1 + Kλ^{2r}} ) λ^{2} ≥ min(1, λ^{2}).

So f(λ) ≤ 1 − η implies λ^{2} ≤ 1 − η; since

f(λ) = 1 − (1 − λ^{2}) (\frac{Kλ^{2r}}{1 + Kλ^{2r}} ) ≥ 1 − (1 − λ^{2})Kλ^{2r} ≥ 1 − Kλ^{2r},

we can apply λ^{2} ≤ 1 − η to the last term in this string of inequalities and conclude

f(λ) ≥ 1 − K(1 − η)^{r};

but f(λ) ≤ 1 − η applied to this last inequality implies K(1 − η)^{r} ≥ η.

Finally, in the case λ_{i} ≤ μ whenever λ_{i} ≠ 1, a similar argument shows that f’s minimum is achieved when all λ_{i} ≠ 1 are equal, say equal λ, and when Q is as large as possible. Thus f({λ_{i}}, {q_{i}}) is lower bounded by

f(λ) = \frac{1 + Kλ^{2r+2}}{1 + Kλ^{2r}} = 1 − \frac{(1 − λ^{2})Kλ^{2r}}{1 + Kλ^{2r}} ≥ 1 − Kλ^{2r} ≥ 1 − Kμ^{2r}.

So f(λ) ≤ 1 − η implies Kμ^{2r} ≥ η.

□

To finish the proof of the former lemma, observe that

\frac{∥Av_{r}∥^{2}}{∥v_{r}∥^{2}} = \frac{λ_{1}^{2r+2}p_{1}^{2} + ∑_{i>1} λ_{i}^{2r+2}p_{i}^{2}}{λ_{1}^{2r+2}p_{1}^{2} + ∑_{i>1} λ_{i}^{2r}p_{i}^{2}} = \frac{1 + ∑_{i>1}(λ_{i}/λ_{1})^{2r+2}(p_{i}/p_{1})^{2}}{1 + ∑_{i>1}(λ_{i}/λ_{1})^{2r}(p_{i}/p_{1})^{2}} ,

and observe that if q_{i} = p_{i}/p_{1} then ∑_{i>1} q_{i}^{2} = L.

□

Combining the first two lemmas of this section easily yields:

Theorem 3.4 The probability that T (v_{r}) ≤ (1 − η)λ_{1}^{2} is at most

η^{−1/2}(1 − η)^{r/2}\sqrt{2n/π}.

This probability is less than δ provided that r ≥ − log_{1−η}(2nδ^{−2}η^{−1}π^{−1}) (or that r ≥ log(2nδ^{−2}η^{−1}π^{−1})/η since − log(1 − η) > η).

The above theorem tells us how soon the Rayleigh quotient T (v_{r}) approaches λ_{1}^{2}. To find the multiplicity of λ_{1}, we need the following result.

Let u_{1}, …, u_{k} be k be orthonormal vectors such that T (u_{i}) ≥ (λ_{1} − ϵ)^{2} for some ϵ > 0. Let P be the projection onto the subspace orthogonal to the u_{i}, and let Ã = P AP.

Theorem 3.5 If A has eigenvalue λ_{1} with multiplicity m, then Ã has eigenvalue λ_{1} with multiplicity m − k. The second eigenvalue, μ_{2}, of Ã satisfies λ_{m+1} ≤ μ_{2} ≤ λ_{m+1} + ϵk.

Proof This is standard and easy. If w_{1}, …, w_{m−k} are orthonormal, orthogonal to the u_{i}, and in the eigenspace of A with respect to λ_{1}, then each w_{i} is clearly an eigenvector of Ã with eigenvalue λ_{1}. This Ã has eigenvalue λ_{1} with multiplicity at least m − k. If v is a unit vector maximizing R(v) = (Ãv, v)/(v, v) subject to being orthogonal to the w_{i}, then clearly P v = v (or else R(P v) < R(v)); hence R(v) = (Av, v), which is the next largest eigenvalue of Ã. An easy calculation shows that for any symmetric matrix, B, and orthonormal s_{1}, …, s_{r} we have

∑_{i=1}^{r}(Bs_{i}, s_{i}) ≤ ∑_{i=1}^{r} ρ_{i},

where ρ_{i} is the i-th largest eigenvalue of B. In particular we have

mλ_{1} + λ_{m+1} ≥ ∑(Aw_{i}, w_{i}) + ∑(Au_{i}, u_{i}) + (Av, v) ≥ k(λ_{1} − ϵ) + (m − k)λ_{1} + (Ãv, v),

which proves μ_{2} ≤ λ_{m+1} + ϵk. But clearly we can find a v in the span of the λ_{1} and λ_{m+1} eigenvectors of A orthogonal to the w_{i} and u_{i}, and for such a v we have (Av, v) ≥ λ_{m+1}; also for such a v we have P v = v, which shows that μ_{2} ≥ λ_{m+1}.

□

Theorem 3.6 Let A have eigenvalue λ_{1} with multiplicity m. Fix an ϵ > 0 and set μ = mϵ + (λ_{m+1}/λ_{1}). Then in m applications of the power method of r iterations each, we will find orthonormal vectors w_{1}, …, w_{m} vectors with T (w_{i}) ≥ λ_{1}^{2}(1−ϵ)^{2} with probability at least 1 − mμ^{r}\sqrt{2n/(πη)}, where η = 2ϵ − ϵ^{2}. This probability is ≥ 1 − δ if

r ≥ (1/2) log(nm^{2}/(ηδ^{2}))/(1 − μ).

Remark 3.7 Clearly all the theorems and lemmas in this section hold for arbitrary symmetric matrices A, provided that (1) we label the eigenvectors in increasing absolute value, |λ_{1}| ≥ |λ_{2}| ≥ ⋯ ≥ |λ_{n}|, and (2) all statements made about a λ_{i} are replaced with |λ_{i}|.

4 Algorithms and Theorems

Recall that we are computing the Betti number, b_{i}, as the dimension of the kernel of Δ_{i}, the i-th Laplacian, which is a rather sparse, n_{i}×n_{i} positive semi-definite symmetric matrix, where n_{i} is the number of faces of dimension i.

In this section we use the power method to compute the multiplicity of the eigenvalue 0 of Δ_{i}, thereby computing b_{i}. The algorithm actually finds the number of eigenvalues of Δ_{i} less than a prescribed parameter; to verify that all these eigenvalues are 0, we can try an “integralization algorithm” for a rigourous verification (we can also use a lower bound on λ_{1} the first non-zero eigenvalue of Δ_{i}, if such a bound is available). We make some further remarks about computing other eigenvalues of Δ_{i} as well.

First we must construct Δ_{i}; let d_{i} denote its average degree (i.e. average number of nonzero entries per row). Let X be a simplicial complex, and let X_{i} be the i-faces. We assume that each X_{i} is specified as a binary word of length n_{0} in the obvious way (from some ordering of X_{0}).

Proposition 4.1 The matrix Δ_{i} can be obtained from a list of X_{i} and X_{i+1} in space O(n_{i}d_{i}) and time T = O(n_{i}(n_{0} − i)(i + 1) log N) where N = max(n_{i−1}, n_{i}, n_{i+1}).

Proof We begin by sorting the X_{i} and X_{i+1}. Since clearly n_{i+1} ≤ (n_{0} − i)n_{i} (considering the map from X_{i+1} to X_{i} of deleting the smallest vertex), the sorting takes time O(n_{i}(n_{0} − i) log N).

For α, β ∈ X_{i}, consider the Δ_{i} entry corresponding to α, β. If α = β, this entry is equal to the number of γ ∈ X_{i−1} such that γ ⊂ α plus the number of γ ∈ X_{i+1} such that α ⊂ γ; one can count this number in n_{0}(log(max(n_{i−1}, n_{i+1}))). If α ≠ β, then this entry is zero if β does not differ from α by one element (i.e. if |α ∪ β| ≠ i + 1); if β differs from α by one element (of which there are (i + 1)(n_{0} − i − 1) such β’s for a fixed α), then the entry is zero if α ∪ β ∈ X_{i}, and is ±1 otherwise (depending on the position of the elements where the ordered α, β differ). So this entry can be found in time O(log(n_{i+1})) for each α, β (determining the position of the elements where α, β differ comes at O(1) cost if for fixed α we loop through the possible β appropriately).

□

Next we describe two algorithms for Betti numbers, using the power method. The first algorithm just tests if a Betti number is nonzero; we will get a sharp bound on the number of power method iterations needed. The second algorithm (which is slightly more complicated) finds a Betti number under certain assumptions. Both algorithms have two phases, and both begin by picking a δ > 0 and then ensure that each phase fails with probability ≤ δ/2. Both algorithms have the same first phase, which is to bound Δ_{i}’s largest eigenvalue.

Let Δ_{i}’s first (i.e. smallest) non-zero eigenvalue be λ_{1} (notice the change in notation from the previous section), and its largest eigenvalue be λ_{max}. We start by finding with probability ≥ 1 − (δ/2) a good approximation of λ_{max}. We apply the power method to Δ_{i}, simply applying it

r = ⌈log_{9/4}(n_{i}δ^{−2}72/(5π))⌉

times and setting B̃ = \sqrt{T (v_{r})} (notation as in the previous section); applying theorem 3.4 taking η = 5/9 and δ replaced here by δ/2, we conclude that λ_{max} ≤ 3B̃/2 with probability ≥ 1 − (δ/2). Of course, B̃ ≤ λ_{max}.

If X_{1} is empty and X is just a collection of vertices, it is trivial to compute the Betti numbers (b_{0} = n_{0}, b_{i} = 0 for i ≥ 1). If not, then Δ_{i} contains a positive integer somewhere along its diagonal. Hence λ_{max} ≥ 1. It follows that if we round B̃ up to the next integer, B, we have B ≤ 2λ_{max}, and λ_{max} ≤ 3B/2.

Next set A = B · I − Δ_{i} where I is the identity matrix. b_{i} is the multiplicity of B as an eigenvalue of A. A’s smallest eigenvalue, μ_{min}, is at least −B/2; let μ_{2} denote A’s largest eigenvalue < B. Now we apply the power method to A. If T (v_{r}) does not seem to converge to B^{2}, we can stop according to theorem 3.4; i.e. when it tells us that with probability ≥ 1 − (δ/2) A does not have an eigenvalue as large as B.

Specifically, let η_{j} for j = 0, 1, … be determined via T (v_{j}) = (1 − η_{j})B^{2}. Our first algorithm applies the power method to A and stops at the r-th iteration provided that

r ≥ − log_{1−η_{r}}(n_{i}δ^{−2}(1 − η_{r})^{−1}8/π).  (4.1)

At this point we know that A has B as an eigenvalue with probability at most δ/2. Overall, then, we conclude that with probability ≤ δ the Betti number, b_{i}, is positive.

The question, of course, is given that b_{i} = 0, how many iterations will it take until r and η_{r} satisfy equation 4.1.

Theorem 4.2 Assume b_{i} = 0. Set ν = (3/4)(λ_{1}/λ_{max}) with λ_{1}, λ_{max} as before. Then η_{i} ≥ ν for all i. Hence our algorithm detects that b_{i} = 0 with probability at least 1 − δ after at most

⌈log_{9/4}(n_{i}δ^{−2}72/(5π))⌉ + ⌈− log_{1−ν}(n_{i}δ^{−2}(1 − ν)^{−1}8/π)⌉

iterations.

Proof It suffices to prove the first claim. We have

T (v_{r}) ≤ max(μ_{2}^{2}, μ_{min}^{2}),

so that T (v_{r}) ≤ (1 − η)B^{2} where

η = min(1 − (μ_{min}^{2}/B^{2}), 1 − (μ_{2}^{2}/B^{2})) ≥ min(3/4, 1 − (B − λ_{1})^{2}/B^{2})
= min(3/4, 2(λ_{1}/B) − (λ_{1}/B)^{2}).

This last expression is

≥ min(3/4, (3/2)(λ_{1}/B))

since for x ∈ [0, 1] either x ≥ 1/2, in which case 2x − x^{2} ≥ 3/4, or x < 1/2, in which case 2x − x^{2} = x(2 − x) ≥ (3/2)x. Since λ_{max} ≤ 2B we have

min(3/4, (3/2)(λ_{1}/B)) ≥ min(3/4, (3/4)(λ_{1}/λ_{max})) = ν,

and the theorem follows.

□

Our first main theorem easily follows:

Theorem 4.3 Assume b_{i} = 0, and let λ_{1} > 0 and λ_{max} be the smallest and largest eigenvalues of Δ_{i}. Set ν = min(λ_{1}/λ_{max}, 3/4). For any prescribed δ > 0 the above algorithm will determine that b_{i} = 0 with probability ≤ 1 − δ in S = O(n_{i}d_{i}) space and O(Sr) operations, where r = O(log(n_{i}δ^{−1}ν^{−1})/ν).

Proof Each iteration requires a multiplication by A, requiring O(S) operations, and an evaluation of T requiring O(n_{i}) = O(S) operations.

□

Next we describe the second algorithm. When b_{i} ≠ 0, then in performing the power method on A we will (with high probability) see T (v_{i}) approaching B. We will pick a point to stop iterating; namely we will fix a small ϵ > 0 and stop iterating when T (v_{r}) ≥ B^{2}(1−ϵ)^{2}. After m applications of the power method we find u_{1}, …, u_{m} orthonormal with T (u_{i}) ≥ B^{2}(1−ϵ)^{2} for each i. We then perform one more application of the power method (to Ã = P AP with P as in section 3), and have T (v_{r}) ≤ (ν_{2}+mϵ)^{2} no matter how large r is. Assuming that ν_{2} + mϵ is less than B(1 − ϵ) (perhaps much less than, relative to ϵ) we will conclude that b_{i} = m.

Theorem 4.4 Let λ_{1} > 0 and λ_{max} be the smallest non-zero and largest eigenvalues of Δ_{i}. Set ν = (3/4)(λ_{1}/λ_{max}). Let ϵ, δ > 0 be given. Then in O(n_{i}(d_{i} + b_{i})) space and O(Sb_{i}r) operations, where r = O(log(n_{i}δ^{−1}ϵ^{−1}α^{−1})/α) with α = ν − 3b_{i}ϵ/2, the above algorithm will compute b_{i} correctly with probability ≥ 1 − δ under the assumption that α ≥ 3ϵ/2.

We remark that we can simplify the theorem by setting ν = λ_{1}/λ_{max} and insisting that ν ≥ 2b_{i}ϵ; then r becomes O(log(n_{i}δ^{−1}ϵ^{−1}ν^{−1})/ν). Proof Say that u_{1}, …, u_{m} with T (u_{i}) ≥ B^{2}(1 − ϵ)^{2} have been found, and we are iterating Ã = P AP. First of all, n_{i}m storage is required to store the u_{i}. Then applying P requires O(mn_{i}) operations, while applying A requires O(n_{i}d_{i}) operations. If we apply the same number, r, iterations to find each u_{j}, j = 1, …, b_{i}, we require a number of operations proportional to

∑_{j=1}^{b_{i}} n_{i}(d_{i} + j) ≤ n_{i}(d_{i} + b_{i})b_{i},

and O(n_{i}b_{i}) storage (in addition to the O(n_{i}d_{i}) required to store Δ_{i}). Next, by theorem 3.6 with probability ≥ 1 − (δ/2) we need r = O(log(nb_{i}ϵ^{−1}δ^{−1}(1 − μ)^{−1})/(1 − μ)) iterations for each u_{i}, where

μ = b_{i}ϵ + max(μ_{2}, |μ_{min}|)/B ≤ b_{i}ϵ + max(1 − (λ_{1}/B), 1/2).

Hence

1 − μ ≥ min((1/2)(λ_{1}/λ_{max}), 1/2) − b_{i}ϵ ≥ (2/3)ν − b_{i}ϵ.

The final phase of the algorithm, i.e. applying the power method after having found b_{i} approximate eigenvectors, requires (by the same analysis as in theorem 4.3) r = O(log(n_{i}δ^{−1}α^{−1})/α iterations to be able to assert that with probability ≥ 1 − δ we have B has multiplicity ≤ b_{i} in A.

□

Next we describe the “integralizing” procedure to rigourously verify that b_{i} ≥ m for some m. Say that our power method algorithm computes that b_{i} ≥ m, meaning that it has found orthonormal vectors u_{1}, …, u_{m} with T (u_{j}) ≥ B^{2}(1 − ϵ)^{2}. To check that the kernel of Δ_{i} is at least m-dimensional, we proceed as follows. By Gauss-Jordan elimination, and by renumbering the coordinates, we can transform the {u_{j}} into a basis {v_{j}} where for any 1 ≤ k ≤ m we have v_{j}’s k-th coordinate is 1 if j = k and 0 if j ≠ k. If the {u_{j}} were a basis for the kernel of Δ_{i}, then by Cramer’s rule all coordinates of the v_{j}’s would be rational numbers. If each v_{j} has an approximate common denominator, meaning an integer D_{j} > 0 such that D_{j}v_{j} is nearly integral, then we can let w_{j} be the “integralized” (i.e. rounded) D_{j}v_{j}. The w_{j} will still be linearly independent; we can check Δ_{i}w_{j} = 0 by integer arithmetic, getting a definite answer. If we are able to verify this, then we have a rigourous proof that b_{i} ≥ m; of course, our main algorithm tells us that b_{i} ≤ m with probability ≥ 1 − δ.

We we shall see in section 6, this integralizing procedure does not always work. The problem is that the approximate common denominators may be large, requiring extra precision in the computation. This leads us to the following question.

Question 4.5 Let w_{1}, …, w_{m} be vectors in R^{n} whose entries are integers between −D and D for some integer D, and let S denote their span intersected with the unit sphere. Let u_{1}, …, u_{m} be chosen inductively with u_{i} chosen from S_{i} uniformly, with S_{i} the intersection of S with the orthogonal complement of u_{1}, …, u_{i−1}. How much precision (and how much time and space) is required to reconstruct an integral basis for S?

In this question we allow the reconstructed integral basis to have coefficients larger than D in absolute value if need be.

Finally we describe an algorithm to compute other eigenvalues of Δ_{i}; we use it in section 6. Namely, to compute the multiplicity of λ in Δ_{i}, we compute B as before and set C = R − (Δ_{i} − λ)^{2}, where R = max(λ^{2}, ((3/2)B − λ)^{2}). (If we perform more iterations to find B, so that B is with high probability very close to λ_{m}x, we can replace the 3/2 by 1.) Then the multiplicity of R in C is that of λ in Δ_{i}; the next largest eigenvalue of C, which controls the convergence rate and therefore speed of the algorithm, is just R − D^{2}, where D is the distance of λ to the nearest other eigenvalue. We do not analyze this algorithm here, but its analysis is similar to that of the above algorithms.

Also, if we expect that Δ_{i} has, in addition to λ, eigenvalues λ_{1}, …, λ_{r}, then we can multiply the initial vector, v_{1}, in the power method by (Δ_{i} − λ_{1}) ⋯ (Δ_{i} − λ_{r}) (to eliminate the λ_{j} components of v_{1}). This greatly sped up this algorithm for the chessboard complexes of section 6.

5 Experiments on A_{k}^{4,2}

i = n_{i} = λ_{max}  λ_{1} λ_{1}/λ_{max} δ = 10^{−2} δ = 10^{−10}
0 16 16.0000 0.  0.
1 120 16.0000 16.0000 1.0000  32  104
2 560 16.0000 12.0000 .7500  43  133
3 1796 16.0000 8.0000 .5000  62  185
4 4080 16.0000 4.0362 .2522  114  334
5 6520 16.0000 4.0362 .2522  117  337
6 7104 16.0000 1.9632 .1227  223  650
7 4962 16.0000 0.  0.
8 1984 16.0000 1.9175 .1198  213  650
9 376 14.0000 3.3316 .2379  104  337
10 16 11.0000 11.0000 1.0000  28  100

Table 1: Data for A_{0}^{4,2}.

Sometimes we are interested in just verifying that certain Betti numbers vanish. As we mentioned in the introduction, there is a large collection of simplicial complexes for which b_{0} = 1 and only one other Betti number does not vanish. In this case it suffices to check the non-vanishing of all but one of the b_{i}’s with i ≥ 1 to determine all the b_{i}’s. The data we present shows that for the complexes studied here, ν = λ_{1}/λ_{max} is large enough to give a good performance for the algorithm.

We are studying the simplicial complex A_{k}^{4,2} of [Fri92] which is defined as follows. We consider all subsets of the four-dimensional Boolean cube, B^{4}, which do not contain more than k 2-dimensional subcubes. So these are simplicial complexes on n_{0} = 16 vertices.

The data presented is the following. For each i we list n_{i}, then λ_{max}(Δ_{i}), λ_{1}(Δ_{i}), and then the ratio of λ_{1}/λ_{max}. After that we calculate how many iterations^{3} are needed

3Notice that by inspecting Δ_{i} to see if it is a multiple of the identity matrix, we don’t really need

i = n_{i} = λ_{max}  λ_{1} λ_{1}/λ_{max} δ = 10^{−2} δ = 10^{−10}
0  16 16.0000 0.  0.
1 120 16.0000 16.0000 1.0000  32  104
2 560 16.0000 16.0000 1.0000  35  107
3 1820 16.0000 16.0000 1.0000  37  109
4 4368 16.0000 8.0000 .5000  65  189
5 7912 16.0000 3.1547 .1971  146  421
6 10560 16.0000 2.9539 .1846  157  450
7 9762 16.0000 1.5637 .09773  280  809
8 5632 16.0000 0.  0.
9 1672 16.0000 2.3843 .1490  173  529
10 208 15.0000 4.0000 .2666  91  302

Table 2: Data for A_{1}^{4,2}.

before we can be sure with probability ≥ 1 − δ that b_{i} = 0; this number is given in theorem 4.2. For both A_{0}^{4,2} and A_{1}^{4,2} there is only one Betti number other than b_{0} = 1 which does not vanish, and so we can calculate it via the Euler characteristic.

The data for A_{i}^{4,2}, i = 0, 1 shows certain curious features. For example, λ_{max} = 16 for a very large number of the Δ_{i}, and λ_{max} never exceeds 16. The author has no idea why this happens, and is very curious to know if this is a property of the specific complex, or if it reflects some general truth about combinatorial Laplacians.

We also remark that in calculating the λ_{1} and λ_{max}, we kept iterating until these values converged within double precision; this typically took a few hundred iterations, and never took more than a few thousand.

For practical interest, we also include a short table of the time taken to compute 1000 iterations of the power method applied to B − Δ_{i}. There we list i, n_{i}, the sum of the degrees of all vertices of the Laplacian, and the time per 1000 iterations. As expected, for comparable sum of degrees, the Laplacian with n_{i} smaller iterates faster. The times is CPU seconds measured on a DEC 3000/400 workstation.

6 Experiments on Chessboard Complexes

For m, n positive integers, consider an m × n chessboard, and consider those subsets of the squares such that no two squares lie on the same column or row (i.e. if rooks were placed on the squares, no two rooks could take each other); this collection of subsets, C(m, n), is closed under taking subsets. C(m, n) is the chessboard complex of

any iterations when it is, i.e. when λ_{1} = λ_{max}.

i = n_{i} = deg sum time
2 120  120 .13
3 560  560 .65
4 1820  1820 2.25
5 4368  7248 7.41
6 7912 33736 25.68
7 10560 90304 57.83
8 9762 128354 74.69
9 5632 91808 50.96
10 1672 23464 11.45
11 208  2704 1.28

Table 3: CPU Seconds per 1000 iterations.

[Gar79],[BLVZ̆94], however it appears in some combinatorial geometric problems, such as the Colored Tverberg’s Problem, as in [Z̆V92] (see also [ABFK92] for applications). The k-connectivity of C(m, n) was a key fact in [Z̆V92] for certain values of m, n, k; in general, k-connectivity (for k ≥ 1) is equivalent to the triviality of π_{0} and π_{1} and the vanishing of H_{i} for all 1 ≤ i ≤ k. In [BLVZ̆94] the connectivity of C(m, n) was studied; there it was proven that C(m, n) is v −2-connected for v = min(m, n, ⌊(m+n+1)/3⌋). It was conjectured that C(m, n) is not v − 1 connected, i.e. that H_{v−1}(C(m, n)) ≠ 0. This was proven for n ≥ 2m−1 and all m, n with m ≤ n and m ≤ 5, with the exception of (m, n) = (4, 6), (5, 7), (5, 8). The conjecture holds if b_{v−1}(C(m, n)) > 0.

We computed b_{2}(C(4, 6)). It has n_{2} = 480, d_{2} = 13. We found λ_{max}(Δ_{2}) = 16. We took ϵ = 10^{−10} and computed b_{2} = 5 (and λ_{1} = 1). Upon integralizing we verified that b_{2} ≥ 5, finding a basis of 0, ±1 valued vectors in ker(Δ_{2}); of course, we know that b_{2} < 6 with high probability from the algorithm, so that b_{2} = 5 with high probability.

We similarly computed b_{3}(C(5, 7)) = 98 (n_{3} = 4200, d_{3} = 17) and b_{3}(C(5, 8)) = 14 (n_{3} = 8400, d_{3} = 21). Integralizing we found a 0, ±1 basis for ker(Δ_{3}) for C(5, 8), but we failed to find an integral basis for C(5, 7). Thus we rigorously know b_{3} ≥ 14, and know b_{3} = 14 with high probability for C(5, 8); we know b_{3} < 99 with high probability, and b_{3} = 98 under the assumption that λ_{1} > 10^{−8}, for C(5, 7).

We began to suspect that the eigenvalues of any Δ_{i} on C(n, m) are all integers, and we verified this for C(m, n) with m = 3 and n = 3, …, 7, with m = 4 and n = 4, 5, 6, and with m = n = 5. We conjectured this suspicion until its recent proof, in [FH]; there the other calculations made in this section can be verified.

7 Worst Case Eigenvalue Separation

In this section we give a number of arguments which indicate that λ_{1} and/or ν = (3/4)λ_{1}/λ_{max} of Δ_{i}, will be bounded away from zero. Most of the results are worst case estimates, and are not very optimistic; we believe that in many applications the true λ_{1} and/or ν will be much better than given here (see, for example, section 5 and 6).

A major problem which is left open is to bound λ_{1} and ν of Δ_{i} for all i in terms of the n_{i} and d_{i}. This could be done in a worst case, restricted, or probabilistic setting.

7.1 Bounds Based on Random Graphs

Consider Δ_{0}, which just depends on the 0- and 1-faces of the complex, which is just a graph. It is well known that if X is a random graph on n vertices of degree d, then λ_{1}(Δ_{0}) = d − \sqrt{d}(2 + o(1)) with high probability in many situations (e.g. d ≥ O(log^{2} n) and n large or d fixed and even and n large) (see [Fri91]). It follows that ν ≥ 1 − d^{−1/2}(2 + o(1)) for most graphs.

7.2 Bounds Based on Worst Case Graphs

We can get some bounds on Δ_{0} just from graph theory bounds. Indeed, from [Fri94], it follows that for a connected graph on n nodes we have

λ_{1}(Δ) ≥ 2 − 2 cos(π/n) = π^{2}/n^{2} + O(n^{−4})

where Δ is the graph Laplacian. Since Δ_{0} is just a graph Laplacian, it follows that the same bound holds for Δ_{0} (and n = n_{0}). This also gives ν proportional to 1/n^{2}; we don’t know if this is the worst possible ν for a graph on n vertices.

Sometimes a similar bound will hold of Δ_{d} for a d-dimensional complex (meaning that d is the top dimension of a face). Recall that Δ_{d} depends on an ordering of each d-face. We say that two faces, α, β, of size d intersect if |α ∩ β| = d (i.e. Δ_{d} is non-zero at the α, β entry), and that they intersect negatively if α and β have differing signs in the α ∩ β row of ∂_{d} (i.e. Δ_{d} is −1 at the α, β entry).

Definition 7.1 We say that a d-dimensional simplicial complex is orientable if there is an ordering of the d-faces such that any two intersecting d-faces intersect negatively.

In particular, a complex based on an orientable manifold is orientable. Also, in an orientable complex, each (d − 1)-face is incident on at most two d-faces.

For a d-dimensional orientable complex we clearly have that Δ_{d} is a graph Laplacian, where the vertices are the d-faces and the vertices the (d−1)-faces. To summarize we have:

Theorem 7.2 We have that λ_{1}(Δ_{0}) ≥ 2 − 2 cos(π/n_{0}); for an orientable complex we have λ_{1}(Δ_{d}) ≥ 2 − 2 cos(π/n_{d}).

Notice that for d = 2, λ_{1}(Δ_{1}) is bounded below by the min of that of λ_{1}(Δ_{0}) and λ_{1}(Δ_{2}). It follows that:

Corollary 7.3 For an orientable complex of dimension d ≤ 2 we have that for all i we have λ_{1}(Δ_{i}) ≥ 2 − 2 cos(π/ñ), where ñ = max(n_{0}, n_{d}).

7.3 A General Conjecture

Although we can’t in general bound λ(Δ_{i}) for all i when d ≥ 3, we can give some arguments which suggest how λ_{1}(Δ_{i}) changes if we take a fixed simplicial complex and subdivide it more and more finely.

For one thing, it is known (see [Cha84]) that for continuous Laplacians acting on the i-forms of d-dimensional manifolds, we have that the λ_{j} are infinite and grow like λ_{n} ≈ cn^{2/d}. It is known that certain types of refinements of combinatorial Laplacians have limits whose eigenvalues converge to those of the continuous Laplacian in a certain sense (see [DP76]). So it makes sense to conjecture that for regular types of refinements (i.e. those where the aspect ratio of sides of a simplex remains bounded with respect to some fixed, smooth metric) we have λ_{n} ≈ cn^{2/d}.

Conjecture 7.4 For regular refinements of a given simplicial complex of dimension d we have ν = (3/4)λ_{1}/λ_{max} = Ω(n_{i}^{−2/d}) as n_{i} → ∞, i.e. as the refinements get finer and finer.

References

[ABFK92] N. Alon, I. Bárány, Z. Füredi, and D. Kleitman. Point selection and weak ϵ-nets for convex hulls. Combinatorics, Probability, and Computing, 1:189– 200, 1992.

[BLVZ̆94] A. Björner, L. Lovász, S.T. Vrećica, and R.T. Z̆ivaljević. Chessboard complexes and matching complexes. J. London Math. Soc., 49:25–39, 1994.

[Cha84] Isaac Chavel. Eigenvalues in Riemannian Geometry. Academic Press, 1984.

[Cha95] Bernard Chazelle. Computational geometry: A retrospective. In D.-Z. Du and F. Hwang, editors, Computing in Euclidean Geometry (2nd Edition), pages 22–46. World Scientific Press, 1995.

[DC91] Bruce R. Donald and David R. Chang. On computing the homology type of a triangulation. In 32nd Annual Symposium on Foundations of Computer Science, pages 650–661, 1991. also preprint (of the full version).

[DE93] C.J.A. Delfinado and H. Edelsbrunner. An incremental algorithm for Betti numbers of simplicial complexes. In 9th Annual ACM Conference on Computational Geometry, pages 232–239, 1993.

[Dod76] Jozef Dodziuk. Finite-difference approach to the Hodge theory of harmonic forms. American Journal of Mathematics, 98:79–104, 1976.

[DP76] J. Dodziuk and V.K. Patodi. Riemannian structures and triangulations of manifolds. J. of the Indian Math. Soc., 40:1–52, 1976.

[Eck45] B. Eckmann. Harmonische Funktionen und Randvertanfgaben in einem Komplex. Commentarii Math. Helvetici, 17:240–245, 1944–45.

[FH] J. Friedman and P. Hanlon. On the Betti numbers of chessboard complexes. To appear.

[Fri91] J. Friedman. On the second eigenvalue and random walks in random dregular graphs. Combinatorica, 11:331–362, 1991.

[Fri92] J. Friedman. Probabilistic spaces of Boolean functions of a given complexity: Generalities and random k-sat coefficients. July 1992. Preprint.

[Fri94] J. Friedman. Minimum higher eigenvalues of Laplacians on graphs. October 1994. To appear in Duke J. Math.

[Fri95] J. Friedman. Error bounds on the power method for determining the largest eigenvalue. October 1995. Preprint.

[Gar79] P.F. Garst. Some Cohen-Macaulay complexes and group actions. Ph.D. Thesis, 1979. University of Madison, Wisconsin.

[Hod41] W.V.D. Hodge. The Theory and Applications of Harmonic Integrals. Cambridge University Press, 1941.

[KB79] R. Kannan and A. Bachem. Polynomial algorithms for computing the Smith and Hermite normal forms of an integer matrix. SIAM J. Computing, 8:499– 507, 1979.

[Mun84] James R. Munkrees. Elements of Algebraic Topology. Benjamin/Cummings, 1984.

[Z̆V92] R.T. Z̆ivaljević and S.T. Vrećica. The colored Tverberg’s problem and complexes of injective functions. J. Combin. Theory Ser. A, 61:309–318, 1992.