ar
X
iv
:0
80
3.
21
36
v1
[
co
nd
-m
at.
oth
er]
1
4 M
ar
20
08
Typical kernel size and number of sparse random matrices over GF (q) - a statistical
physics approach
Alamino, R.C., Saad D.
Aston University, Neural Computing Research Group, Birmingham B4 7ET, UK
Using methods of statistical physics, we study the average number and kernel size of general sparse
random matrices over GF (q), with a given connectivity profile, in the thermodynamical limit of large
matrices. We introduce a mapping of GF (q) matrices onto spin systems using the representation
of the cyclic group of order q as the q-th complex roots of unity. This representation facilitates
the derivation of the average kernel size of random matrices using the replica approach, under the
replica symmetric ansatz, resulting in saddle point equations for general connectivity distributions.
Numerical solutions are then obtained for particular cases by population dynamics. Similar tech-
niques also allow us to obtain an expression for the exact and average number of random matrices
for any general connectivity profile. We present numerical results for particular distributions.
PACS numbers: 02.10.Yn, 02.70.-c,05.10.-a
Keywords: random matrices, Galois fields, statistical mechanics, replica theory
I. INTRODUCTION
Random matrices over GF (q) are highly important in a number of application areas ranging from biology to
computer science and telecommunication. One of the areas where they play a particularly important role is coding
theory [1]. In particular, linear codes are defined by the kernel of a parity-check matrix, where each kernel vector
is termed a codeword and is associated with an original uncoded message vector by a linear operation defined by a
generator matrix. Well known examples include the Hadamard codes, where properties of the kernel and rank play an
important role [2], and low-density parity-check codes (LDPC) which provide the best performance to date in many
noise regimes. Although the most studied and applied case of LDPC codes is of binary codes over GF (2) there is a
significant body of work, of both practical and theoretical nature [3], on codes over more general finite fields showing
an improvement in performance with respect to the binary version. In particular, statistical physics based analysis of
LDPC codes over GF (q) has been reported in [4].
Low-density parity-check codes are based on random sparse matrices, where the fraction of non-zero elements goes
to zero as the size of the matrix increases. In most studies of LDPC codes, it is assumed that a parity-check matrix
with M rows (parity-checks) and N columns defines a code of rate R = 1−M/N , exactly, which is equivalent to the
assertion that the number of vectors in the kernel (and therefore the number of codewords) is exactly qNR.
In addition to being an interesting applied problem, the properties of these matrices are also of great interest from
the pure mathematical point of view and a number of papers has already tried to answer related questions in different
instances with a mathematical rigorous approach [5, 6, 7].
In this contribution, we address two key properties of sparse random matrices over GF (q), namely the average
dimension of their kernel and the number of matrices for a given connectivity profile, in the case of large matrices.
When the matrices are large, keeping N → ∞ with M/N constant, the problem can be mapped into a system of
interacting “spins” and the powerful machinery developed for the study of disordered spin lattices in condensed matter
physics can then be used, under some assumptions, to obtain the required properties.
In order to keep this paper as self-contained as possible and make it accessible to a broad readership, we provide in
section II a brief introduction to GF (q) matrices and their properties, and to the basic statistical physics methodology
on which we have based our analysis. The usual statistical physics approach to the analysis of LDPC codes over the
binary field GF (2) is generalized in such a way that it can be efficiently applied to any GF (q) for a general connectivity
distribution of non-zero elements and then used to calculate the average kernel dimension of sparse random matrices
(SRM) in section IV. Making use of techniques developed in section IV, the number of matrices for a given distribution
of non-zero elements is then obtained for various connectivity profiles, in section V. Finally, we present a discussion
of the obtained results in section VI.

2II. KEY CONCEPTS
A. GF (q)-Matrices
A Galois field GF (q) is a finite field with q elements, i.e., a set of q elements {0, ..., q − 1}, which we symbolize by
integers for convenience, which is a commutative group under addition ⊕ : GF (q) → GF (q), defined as integer addition
mod q, and with a monoid structure with respect to a commutative multiplication operation ⊗ : GF (q) → GF (q).
The field also includes the zero element ’0’, mapping every other element to itself, and the identity ’1’; an additional
requirement is that the multiplication and addition have the algebraic distributive property. This last requirement
restricts the number of elements to be q = pn, where p is a prime number and n an integer.
Entries in matrices over GF (q) take values of numbers in the field GF (q), where the usual additions and multi-
plications involved in their algebra are defined by the corresponding operations over the Galois field. The kernel, or
null space, of an M ×N matrix A is defined as the set of vectors v ∈ GF (q)N such that Av = 0, with all operations
in the field GF (q). The kernel is a linear vector space and therefore will have qd(A) vectors, where d(A) is the kernel
dimension. The rank r(A) of the matrix is obtained by the rank-nullity theorem as r(A) = N − d(A).
B. Disordered Systems
An interacting spin problem has two main elements: an interaction defined between a number of spin units,
collectively represented by the vector σ = (σ1, ..., σN ), in a lattice and a local field which acts in each variable σi
separately. Disordered spin systems are systems where one or both of these elements (interaction and field) is a
random variable. Usually, we are interested in the properties of very large systems, where the number N of spins
becomes infinite, the so-called thermodynamic limit.
The main properties of the system in the thermodynamic limit can be derived from a key quantity, the free-energy
f , which in probabilistic terms corresponds to the cummulant generating function. For disordered systems, in the
cases where the free-energy is self-averaging with respect to the disorder, we can calculate this quantity as
f = − lim
N→∞
1
βN 〈lnZ〉, (1)
where 〈·〉 indicates the disorder average, Z =∑σ e−βH(σ) is the partition function and H(σ) is the Hamiltonian of
the system. Although the self-averaging property should be rigorously investigated for each system, we will assume
it holds here.
In order to obtain the free-energy, a powerful technique is to make use of the replica method, based on the identity
[ ∂
∂n ln 〈Z
n〉
]
n=0
= 〈lnZ〉 . (2)
Average quantities can then be calculated for integer n and then analytically continued to zero. The replica theory
is commonly used in the area of disordered systems and is known to provide exact results in many regimes, which
include both physical and non-physical systems [8, 9].
Many problems in computing and communication theory can be mapped to spin systems. For instance, error-
correcting codes, in particular LDPC codes [10] and hard computational problems such as K-SAT [11] and graph-
coloring [12, 13], can be mapped to diluted spin systems with random p-spin interactions and local fields. In the
coding example, interactions are defined by the parity-check constraints, while the local fields are induced by the
codeword and received message. In the statistical physics treatment, for mathematical convenience, the message bits
{0, 1} and ’⊕’ operation are mapped onto spin values {+1,−1} and multiplication using the mapping x → (−1)x.
Variables over a general finite field GF (q), q 6= 2 are typically first mapped onto a binary string and then, using the
spin values representation, transformed into a spin system [4].
III. MAPPING GF (q) MATRICES INTO SPIN SYSTEMS
The transformation
σ(v) = (−1)v, (3)
where σ ∈ {+1,−1} and v ∈ {0, 1}, is usually employed to map the GF (2) variables onto the binary representation.
This mapping can be generalized to any GF (q) without an intermediate use of the binary field.

3Under the operation ⊕, GF (q) is homeomorphic to the cyclic group of order q and therefore has a representation
as the complex q-th roots of unity with the group homeomorphism σ : GF (q) → C given by
σ(v) = exp
(
2πi
q v
)
, (4)
such that for every v1, v2 ∈ GF (q)
σ(v1 ⊕ v2) = exp
[
2πi
q (v1 ⊕ v2)
]
= exp
[
2πi
q (v1 + v2)
]
= exp
(
2πi
q v1
)
exp
(
2πi
q v2
)
= σ(v1)σ(v2).
(5)
This mapping has a clear geometric interpretation: 2πv/q is an angle in the unit circle, such that each element
of the Galois field is being mapped onto a spin variable “pointing” in one of q possible angles. Using this mapping
allows one to write the null-space constraint for a general vector v =
(
v1, ..., vN
)
∈ GF (q)N as
δ(Av, 0) =
M
∏
i=1
δ


N
⊕
j=1
(
Aij ⊗ vj
)
, 0

, (6)
with
δ


N
⊕
j=1
(
Aij ⊗ vj
)
, 0

 =
1
∆(q)
q−1
∏
m=1



1− exp
(
−2πiq m
) N
∏
j=1
exp
[
2πi
q
(
Aij ⊗ vj
)
]



, (7)
and
∆(q) =
q−1
∏
m=1
[
1− exp
(
−2πiq m
)]
. (8)
Using the properties of the complex roots of unity, the above quantity ∆(q) can be shown (see appendix A) to be
real and equal to the order q of the field.
Based on this representation, we can now define the “magnetization” of the original system in analogy with the
spin system as
m = 1N
N
∑
j=1
σj , (9)
and the overlap between two configurations σ and σ′ as
ρ = 1N
N
∑
j=1
σjσ′j , (10)
where we are now working with the spin variables already mapped to the the complex field C and therefore the
operations of multiplication and addition correspond to the usual ones in C.
It turns out that this kind of representation allows a factorization of the terms simplifying the equations and making
the replica calculations simpler, as we will see in the following.
IV. AVERAGE PROPERTIES OF THE KERNEL
The dimension of the kernel of an M ×N matrix A over GF (q) can be written as d(A) = logq Ω where
Ω =
∑
v
δ(Av, 0), (11)

4is the number of vectors in the kernel, δ is the Kroenecker delta and v ∈ GF (q)N . Direct calculation of Ω from
equation (11) by straightforwardly substituting the Kroenecker delta by its integral representation trivially reproduces
the rank-nullity theorem. This calculation is not presented here.
The quantity we are interested in here is the average kernel dimension, more specifically, its density in the limit of
large matrices, defined as Ts where
s ≡ 1T limN→∞
〈d(A)〉A
N = limN→∞
1
N 〈lnΩ〉A, (12)
where 1/T = ln q and M/N ≡ λ, with λ a finite positive constant. Using the replica identity (2), we can write
s = lim
N→∞
[ ∂
∂n ln 〈Ω
n〉A
]
n=0
. (13)
The randomly chosen sparse matrices A have exactly Ki non-zero elements in the i-th row with probability P(K),
K ≡ (K1, ...,KM ), and Cj elements in the j-th column with probability P(C), C ≡ (C1, ..., CN ), obeying the
constraint Λ ≡∑i Ki =
∑
j Cj , where Λ is the total number of non-zero elements of the matrix. The elements of A
are sampled from the finite field GF (q) with independent equal probabilities P(Aij).
Let us define, for brevity of notation, Zn ≡ 〈Ωn〉A. Although the calculations, presented in appendix B, are similar
to related calculations in [14, 15], we will use a different approach which is conceptually clearer and has the advantage
of allowing later generalizations. In this approach, we sum directly over all entries of the matrix instead of defining a
connectivity tensor as used elsewhere [14, 15],
Zn =
〈
1
N
∑
{Aij}


∏
i,j
P(Aij)




M
∏
i=1
δ


N
∑
j=1
χ(Aij),Ki






N
∏
j=1
δ
(
M
∑
i=1
χ(Aij), Cj
)


×
n
∏
a=1
[
∑
va
δ(Ava, 0)
]〉
K,C,Λ
,
(14)
where the average is over the probability distribution P(K,C,Λ) with χ(Aij) = 0 if Aij = 0 and 1 otherwise, and
the normalization N gives the number of matrices which obey the constraints averaged over the distributions of the
entries. In this way, any type of constraint on the matrix can be readily included in the calculation, which could
be rather cumbersome in other approaches, based on the introduction of a connectivity tensor as the corresponding
constraints have to be written in terms of the tensor elements, which can be extremely complicated.
We refer the reader to appendix B for details of the calculations. Using the replica symmetric ansatz, which is
shown to be exact for this problem (see appendix D) we arrive at the following self-consistent saddle point equations
πˆ(xˆ) = 1αǫ(α)
M
∑
i=1
〈
αΛ
Λ!
Kiδ
(
xˆ−
Ki−1
∏
l=1
xl
)〉
x,K,C,Λ
, (15)
π(x) = 1αǫ(α)
N
∑
j=1
〈
αΛ
Λ!
Cjδ
(
x−
∏Cj−1
l=1 [1+(q−1)xˆl]−
∏Cj−1
l=1 (1−xˆl)
∏Cj−1
l=1 [1+(q−1)xˆl]+(q−1)
∏Cj−1
l=1 (1−xˆl)
)〉
xˆ,K,C,Λ
, (16)
0 =
〈αΛ
Λ!
(
1− Λα
)〉
K,C,Λ
, (17)
with
ǫ(α) =
〈αΛ
Λ!
〉
K,C,Λ
, (18)
and to the corresponding expression for s
s = −λ ln q − αN 〈ln [1 + (q − 1)xxˆ]〉x,xˆ
+
1
Nǫ(α)
∑
i
〈
αΛ
Λ!
〈
ln
[
1 + (q − 1)
Ki
∏
l=1
xl
]〉
x
〉
K,C,Λ
+
1
Nǫ(α)
∑
j
〈
αΛ
Λ!
〈
ln



Cj
∏
l=1
[1 + (q − 1)xˆl] + (q − 1)
Cj
∏
l=1
(1− xˆl)



〉
xˆ
〉
K,C,Λ
.
(19)

5It must be noted that the above equations are only meaningful if Λ ∝ N . A striking property of the above equations
is that they are completely independent of the specific distribution of the individual elements of the matrix, depending
only on the distribution of K and C (and, obviously, of Λ).
There exists two straightforward analytical solutions of the above equations, namely, the paramagnetic one given
by
πˆ(xˆ) = δ(xˆ), π(x) = δ(x), (20)
and the ferromagnetic solution
πˆ(xˆ) = δ(xˆ− 1), π(x) = δ(x− 1). (21)
When substituted in the above equations, the paramagnetic solution gives the average kernel density as Ts =
1 − λ = 1 −M/N independently of the order q of the finite field used. In the case of LDPC codes defined by such
matrices, this corresponds to random parity-check matrices that defines a code of rate R = 1− λ. The average rank
density in this case is λ. The ferromagnetic solution gives Ts = 0 and the matrix is full rank; which incidentally
means that such matrices cannot be used to define a parity-check code due to the lack of redundancy.
These quantities can be associated to analogous quantities in the statistical mechanics framework. We start by
associating the average rank density with the free-energy f and writing
f ≡ 〈r(A)〉AN = 1− Ts, (22)
which allows one to associate s with the entropy and the internal energy density being constrained to be u = 1.
Defining β = 1/T , equation (22) becomes
βf = 1− 1N
〈
ln
∑
v
δ(Av, 0)
〉
A
=
1
N
〈
N ln eβ − ln
∑
v
δ(Av, 0)
〉
A
= − 1N
〈
ln
∑
v
δ(Av, 0) e−βN
〉
A
= − 1N
〈
ln
∑
v
e−βH(v)
〉
A
,
(23)
where the Hamiltonian of the corresponding statistical mechanical system is formally
H(v) ≡ N − ln δ(Av, 0). (24)
We solved the saddle point equations by means of population dynamics for three different cases, in all of which we
keep K fixed
1. Regular matrices - C and K fixed;
2. Fixed K and C drawn from a multinomial uniform probability
P(C) = (MK)!∏
j Cj !
1
NMK ; (25)
3. Fixed K while C values are drawn from a Poisson integer distribution of mean Λ/N = λK, for each column
separately, until the limit of MK non-zero elements is reached.
Results for the various cases are presented in Fig. 1. The top left plot shows the theoretical thermodynamically
dominant solutions (paramagnetic in the range 0 ≤ λ ≤ 1 and ferromagnetic for λ > 1) having the lower free energy.
The top right plot shows the results for the regular case (i). Solutions were obtained numerically by iterating
equations (15) and (16) for the case of q = 4 and K = 200; C was varied from 2 to 250. Repeating the calculations
for different values of q and K have produced similar results. We see that the stable solution is always paramagnetic,
but becomes unphysical at λ = 1 once the entropy, and consequently the dimension of the kernel, become negative.

60.0 0.5 1.0
λ
0.0
0.5
1.0
Kernel Dimension
Rank
0.0 0.5 1.0
λ
0.0
0.5
1.0
Kernel Dimension
Rank
0.0 0.5 1.0
λ
0.0
0.5
1.0
Kernel Dimension
Rank
0.0 0.5 1.0
0.0
1.0
Kernel Dimension
Rank
FIG. 1: Average kernel dimension density (continuous lines) and average rank density (dashed lines) calculated as solutions to
the replica symmetric saddle point equations. The top left plot shows the thermodynamically favored solution (paramagnetic
for 0 ≤ λ ≤ 1 and ferromagnetic for λ > 1). The top right shows the regular case (i) for fixed K and C. Cases (ii) and (iii) are
presented at the bottom left and right, respectively. Note that numerical instabilities occur for specific λ values.
In the case of parity-check codes, this result means that the typical parity-check matrix defines a code of rate exactly
(N −M)/N . This is assumed for any parity-check matrix in most calculations in the literature and is confirmed by
our results to be true on average; however, it is important to point out that the result is true in the limit of large
matrices and is likely to have finite size corrections which may affect practical applications.
Cases (ii) and (iii) are presented, respectively, at the bottom left and right of Fig. 1. Although these cases do not
rigorously obey the constraint that each Cj must be at most M , for large matrices and small values of K (which is
what happens in practice) Cj is unlikely to exceed this value. However, instabilities can and indeed occur for specific
λ values, presumably due to instances where Cj takes higher values.
The bottom left plot shows results for the case (ii), with q = 3, K = 4, N = 1000 and 1 ≤ M ≤ 1250. Also in this
case, the stable dominant solution is paramagnetic. Numerical instabilities, which disappear slowly with the increase
in the number of fields and steps in the population dynamics, emerge in the unphysical region and are shown in the
figure.
The behavior for case (iii) is a little more complex due to the nature of the distribution chosen. Using the average
value λK for the variables Cj implies that, as λ varies, their average value also changes. The plot shown was obtained
for q = 2, K = 4, N = 250 and 1 ≤ M ≤ 300. There are clearly special points in this plot, which distinguish it from
the previous cases. The first point separates λ values which give rise to average connectivity values lower/higher than
1 (left and right, respectively). Up to this point, the matrix has too many zero columns, pushing the kernel size to
cover the full space of vectors. The other two points are where numerical instabilities emerge. Further calculations
with different K values indicate that these points appear around the extremes of the interval 2/K ≤ λ ≤ 3/K.
Inside this interval, the average value of the Cj ’s equal to 2 (once we take it to be an integer). This value marks the
percolation transition for binary matrices. Apart from these differences, the resulting curve seems to coincide with
those obtained for the previous cases.
The solution of kernel size problem is mathematically equivalent to the solution of LDPC in channels with infinite
noise. As the solution in the latter is paramagnetic, we are led to speculate that it is the dominant solution also here
up to the point where the quantity s, analogous to the entropy, becomes negative. From this point and on the solution
becomes ferromagnetic. The numerical results seem to support this conjecture, although more careful calculations,
varying all the parameters involved must be carried out to confirm this hypothesis more generally.

7V. NUMBER OF MATRICES
The number of GF (q) matrices given a connectivity profile is of significant interest within the discrete mathematics
community. Exact results have been obtained for the case of finite binary matrices [16] in the form of a formula that
facilitates the calculation of their precise number. In this paper we will analyze the case of large GF (q) matrices and
provide an expression for both their exact and average number. Given the precise number of non-zero elements per
row K = (K1, ...,KM ) and per column C = (C1, ..., CN ), one can write the number of matrices as
NA =
∑
{Aij}


M
∏
i=1
δ


N
∑
j=1
χ(Aij),Ki






N
∏
j=1
δ
(
M
∑
i=1
χ(Aij), Cj
)

. (26)
Note that we are using the summation directly over the entries of the matrix instead of the introduction of a
connectivity tensor. In this way, the calculations are similar to the ones for obtaining the kernel dimension with the
details given in C. The final result is
NA = (q − 1)Λ
Λ!
∏
i Ki!
∏
j Cj !
. (27)
Note that the component on the right represents the number of binary matrices with the given non-zero elements
profile. The factor (q− 1)Λ is the multiplicity of the non-zero entries which can have any non-zero value in the Galois
field.
If we consider a distribution P(K,C,Λ), we can look at the average number of matrices
N¯A =
〈
(q − 1)Λ Λ!∏
i Ki!
∏
j Cj !
〉
K,C,Λ
. (28)
Note that we can write the joint probability distribution as
P(K,C,Λ) = P(K|Λ,C)P(Λ|C)P(C), (29)
and that P(Λ|C) = δ
(
Λ,∑j Cj
)
. Therefore, we have obtained for the average number of matrices
N¯A =
∑
K
∑
C
P(K|C)P(C)(q − 1)
P
j Cj
(
∑
j Cj
)
!
∏
i Ki!
∏
j Cj !
, (30)
where the distribution P(K|C) includes the constraint δ
(
∑
i Ki,
∑
j Cj
)
.
A simple calculation shows that for the regular case, where all Cj ’s and Ki’s are fixed (to C and K, respectively),
and q = 2, the number of matrices scales as NCN . Therefore, a more appropriate quantity to calculate instead of the
average number of matrices would be the quenched entropy
Ξ ≡
〈
1
N lnNA
〉
=
1
N
∑
K
∑
C
P(K|C)P(C) ln

(q − 1)
P
j Cj
(
∑
j Cj
)
!
∏
i Ki!
∏
j Cj !

, (31)
which scales as lnN .
We analyze the behavior of this quantity for three different cases. We choose each Cj to be i.i.d. and K to be
chosen from a multinomial distribution
P(K) = (
∑
i Ki)!
∏
i Ki!
1
N
P
i Ki
δ


∑
i
Ki,
∑
j
Cj

, (32)
for each realization of C. The three probability distributions for the variables Cj to be analyzed are
1. uniform in the interval [0, 2C¯]
P(Cj) = 1/(2C¯ + 1); (33)

8TABLE I: Asymptotic values of Ξ∗ for large λ
C¯ As. Value
5 29.66
10 60.73
20 123.20
2. binomial in the interval [0,M ]
P(Cj) =
(M
Cj
)( C¯
M
)Cj(
1− C¯M
)M−Cj
; (34)
3. Zipf distribution for Cj = 1, ...,M
P(Cj) =
C−sj
∑K
n=1 n−s
, (35)
where C¯ is the mean of the distributions. The motivation for choosing these connectivity profiles is that they appear
to be the most commonly analyzed and feature (especially the latter) in recent analysis and modeling of networks.
Results for the binomial (dashed line) and uniform (dotted line) distributions with means C¯ = 5.0, 10.0, 20.0, q = 2
and N = 300 are plotted in Fig. 2, together with the value of Ξ with constant Cj = C¯ and Kj = C¯/λ values for all i
and j. This function is explicitly given by
Ξ∗ = C¯ ln(q − 1)− ln C¯! + 1N ln (NC)!− λ ln
(
C¯/λ
)
!, (36)
and we can obtain its asymptotic behavior for small and large λ as
λ ≪ 1 ⇒Ξ∗ = C¯ ln(q − 1)− ln C¯! + C¯ lnλN, (37)
λ ≫ 1 ⇒Ξ∗ = C¯ ln(q − 1)− ln C¯! + C¯ ln C¯N + (γ − 1)C¯, (38)
where γ ≈ 0.577216 is the Euler-Mascheroni constant. Asymptotic limits for large λ are given in table I.
For large λ values the result for constant C and K upper-bounds the other two distributions. Additional calculations
seem to indicate that it is always the case for any distribution, although a proof for this conjecture is still sought.
This implies that if we keep the number of columns constant and increase the ratio λ by adding rows, whenever the
number of rows is much larger than the number of columns, the average number of matrices becomes independent of
both the ratio and number of rows. The plots also suggest that the average number of matrices in these cases are
basically defined by the average value of the C distributions.
For small values of λ, the uniform distribution continues to be upper-bounded by the constant distribution. The
binomial distribution, however, is higher for a small interval around zero. This behavior is shown in the inset where
lower C values give rise to higher Ξ as λ becomes smaller.
Figure 3 shows the results for the Zipf distribution with different values for the power s compared with a uniform
distribution in the range [0,M ]. In this case, the mean of the distributions vary with λ. We see that, although the
average value of the Zipf distributions increasingly differs from the uniform value M/2 as s increases, the average
number of matrices actually becomes highly similar.
VI. CONCLUSIONS
We have introduced a new mapping of Galois matrices to spin systems based on the group homeomorphism between
GF (q) under addition mod q (denoted by ⊕) and the complex q-th roots of unity. In addition, we have introduced a
different way for summing over random matrices that can be generalized to include any kind of connectivity constraint
and is conceptually cleaner and simpler than the existing approaches. The new mapping and alternative summation
over random matrices allows for a factorization of the constraints, which simplifies calculations of the kernel and the
number of matrices under various connectivity profiles.
Using the replica approach and these new introduced techniques, we calculated the average dimension of the kernel
for a general distribution of non-zero entries and solved the resulting equations numerically, finding that the average

90.0 0.5 1.0 1.5 2.0
λ
-40
-20
0
20
40
60
80
Ξ 5.0
10.0
20.0
Constant C and K
Binomial Distribution
Uniform Distribution
0 0.07
0
20
FIG. 2: Values of the quenched entropy Ξ versus λ for the different distributions and various C values (C = 5, 10, 20), with
multinomial K: constant (continuous line), binomial (dashed line) and uniform (dotted line). The inset shows in detail the
small λ regime, where just the binomial and constant distributions are represented. The higher lines on the right correspond
to the higher C values.
0.0 0.5 1.0 1.5 2.0
λ
0
100
200
300
400
500
Ξ
FIG. 3: Values of Ξ versus λ for the uniform distribution (dashed line) and the Zipf distribution (continuous lines) for
s = 1, 3, 4, 10, respectively, from bottom to top.
kernel density is 1 −M/N in all cases studied. We conjecture that this result is always valid. Based on the analogy
with thermodynamical quantities corresponding to free energy, internal energy and Hamiltonian, we showed that
the replica symmetric ansatz in this case must be exact. With the same techniques, we were also able to find the
total number of large matrices for fixed K and C and their average number, which was then computed for different
distributions of theoretical and practical relevance.
The results presented have practical relevance in a number of areas, including coding network modeling and some
biological models. With respect to LDPC codes, the average kernels density result implies that randomly generated

10
LDPC codes typically define codes of rate exactly 1−M/N , an assumption which is generally made but lacks rigorous
derivations. Also, as the parity-parity check matrix can represent the connectivities in graphs (see [17]), the results
obtained for the average number of matrices provide a principled approach to determine the average number of possible
graphs with a given connectivity distributions of a more general nature than the connectivity profiles examined in
this paper.
Acknowledgements
Support from EPSRC grant EP/E049516/1 is gratefully acknowledged. R.C.A. would also like to thank Dr. Juan
P. Neirotti for useful discussions.
APPENDIX A: PROOF OF ∆(q) = q
In this appendix we prove the statement made in section IV that ∆(q) = q where
∆(q) =
q−1
∏
m=1
[
1− exp
(
−2πiq m
)]
. (A1)
Let us use the notation
z(m) ≡ e 2piiq m, (A2)
and noting that unit complex roots appear in complex conjugate pairs, we write
∆(q) =





∏(q−1)/2
m=1 [1− z(m)][1− z¯(m)], q odd,
2
∏(q−2)/2
m=1 [1− z(m)][1− z¯(m)], q even,
(A3)
where the bar indicates a complex conjugate. Using
[1− z(m)][1− z¯(m)] = 2− 2Re[z(m)]
= 2
[
1− cos
(
m2πq
)]
= 4 sin2
(
mπq
)
,
(A4)
equation (A3) becomes
∆(q) =









2q−1
[
∏(q−1)/2
m=1 sin
(
mpiq
)]2
, q odd,
2q−1
[
∏(q−2)/2
m=1 sin
(
mpiq
)]2
, q even.
(A5)
As the sin function is positive in the interval (0, π) and sin(π/2) = 1 we can write, for any q,
∆(q) = 2q−1
q−1
∏
m=1
sin
(
mπq
)
. (A6)
Using the known identity [18]
sin(qx) = 2q−1
q−1
∏
m=0
sin
(
x+mπq
)
, (A7)
divided by sinx and taking x → 0, one obtains
q−1
∏
m=1
sin
(
mπq
)
=
q
2q−1
, (A8)
which by substituting into equation (A6) gives the desired result.

11
APPENDIX B: REPLICA SYMMETRIC SADDLE POINT EQUATIONS
Using integral representations for the first two sets of Kroenecker delta functions, we can write the averaged
replicated kernel size defined in equation (14) as
Zn =
〈
1
N
∑
{va}
∮
DWDZ
∑
{Aij}


∏
i,j
P(Aij)(WiZj)χ(Aij)


×
M
∏
i=1
∏
a
δ


N
⊕
j=1
(
Aij ⊗ vja
)
, 0


〉
K,C,Λ
,
(B1)
where ⊗ and ⊕ indicate multiplication and summation on GF (q), respectively, and
DWDZ =
[
M
∏
i=1
dWi
WKi+1i
]


N
∏
j=1
dZj
ZCj+1j

. (B2)
Using the representation of the parity-check constraint given in equation (6), the product over replica indices of the
delta function can be written as
∏
a
δ


N
⊕
j=1
(
Aij ⊗ vja
)
, 0

 =
∏
a
1
q
q−1
∏
m=1



1− exp
(
−2πiq m
) N
∏
j=1
exp
[
2πi
q
(
Aij ⊗ vja
)
]



=
1
qn
∏
a
[
1 +
q−1
∑
s=1
Fi(s, a)G(s)
]
=
1
qn
n
∑
r=0
∑
〈a1···ar〉
∑
s1,...,sr
G(s1) · · ·G(sr)Fi(s1, a1) · · ·Fi(sr, ar),
(B3)
with
G(s) ≡
∑
〈m1···ms〉
(−1)s exp
(
−2πiq m1
)
· · · exp
(
−2πiq ms
)
, (B4)
and
Fi(s, a) ≡ exp
[
2πi
q
(
Ai1 ⊗ v1a
)
]
· · · exp
[
2πi
q
(
AiN ⊗ vNa
)
]
=
N
∏
j=1
γj(s, a, Aij),
(B5)
where we defined, for simplicity,
γj(s, a, Aij) ≡ exp
[
2πi
q s
(
Aij ⊗ vja
)
]
. (B6)
We can now write the partition function as
Zn =
〈
1
N
∑
{va}
∮
DZ
M
∏
i=1
1
qn
n
∑
r=0
∑
〈a1···ar〉
∑
s1,...,sr
G(s1) · · ·G(sr)
∮ dWi
2πi
1
WKi+1i
Γi
〉
K,C,Λ
, (B7)

12
where
Γi =
∑
Ai1,...,AiN


∏
j
P(Aij)(WiZj)χ(Aij)


∏
j
γj(s1, a1, Aij) · · · γj(sr, ar, Aij)
=
∏
j
∑
Aij
P(Aij)(WiZj)χ(Aij)γj(s1, a1, Aij) · · · γj(sr, ar, Aij)
= pN
∏
j
[
1 +
1
p
q−1
∑
h=1
P(Aij = h)WiZjγj(s1, a1, h) · · · γj(sr, ar, h)
]
,
(B8)
where we define, for convenience, p ≡ P(Aij = 0). Let us define a probability distribution over the values of h as
P(h) = P(Aij = h)
1− p , (B9)
in such a way that h varies from 1 to q − 1 and the probability over this range is correctly normalized. Then
Γi = pN
∏
j
[
1 +
(
1− p
p
)
WiZj〈γj(s1, a1, h) · · · γj(sr, ar, h)〉h
]
= pN
N
∑
l=0
∑
〈j1···jl〉
(
1− p
p
)l
W liZj1 · · ·Zjl
× 〈γj1(s1, a1, h) · · · γj1(sr, ar, h)〉h · · · 〈γjl(s1, a1, h) · · ·γjl(sr, ar, h)〉h.
(B10)
The integrals over the Wi’s, acting on the Γi’s, select the power of Wi to be Ki and we therefore obtain
Zn =
〈
κ
∑
{va}
∮
DZ
M
∏
i=1





n
∑
r=0
∑
〈a1···ar〉
∑
s1,...,sr
G(s1) · · ·G(sr)
∑
〈j1···jKi〉
Zj1 · · ·ZjKi
× 〈γj1(s1, a1, h) · · · γj1(sr, ar, h)〉h · · ·
〈
γjKi (s1, a1, h) · · · γjKi (sr, ar, h)
〉
h
}〉
K,C,Λ
≈
〈
κ
∑
{va}
∮
DZ
M
∏
i=1



n
∑
r=0
∑
〈a1···ar〉
∑
s1,...,sr
G(s1) · · ·G(sr)
× N
Ki
Ki!


1
N
N
∑
j=1
Zj〈γj(s1, a1, h) · · · γj(sr, ar, h)〉h


Ki





〉
K,C,Λ
(B11)
where
κ = pNM
(
1− p
p
)
P
i Ki
N−1q−nM . (B12)
The calculation of N is similar to the calculation of the number of matrices shown in appendix C and we end up
with
κ = 1
qnMN (2)A
, (B13)
where N (2)A is exactly the number of binary matrices (q = 2) as calculated in appendix C. Introducing the replica
overlaps
Qs1,...,sr〈a1···ar〉 ≡
1
N
N
∑
j=1
Zj〈γj(s1, a1, h) · · · γj(sr, ar, h)〉h, (B14)

13
and the corresponding auxiliary variables Qˆs1,...,sr〈a1···ar〉 by means of Dirac delta functions, we can express the partition
function as
Zn =
∫
DQDQˆ exp
(
−N
∑
Qs1,...,sr〈a1···ar〉Qˆ
s1,...,sr
〈a1···ar〉
)
×
〈
κN
P
i Ki
∏
i Ki!
∏
i
[
∑
G(s1) · · ·G(sr)
(
Qs1,...,sr〈a1···ar〉
)Ki
]
×
∏
j





∑
{vja}
∮
DZj exp
[
Zj
∑
Qˆs1,...,sr〈a1···ar〉〈γj(s1, a1, h) · · ·γj(sr, ar, h)〉h
]





〉
K,C,Λ
=
∫
DQDQˆ exp
(
−N
∑
Qs1,...,sr〈a1···ar〉Qˆ
s1,...,sr
〈a1···ar〉
)
×
〈
q−nM N
P
i Ki
(
∑
i Ki)!
∏
i
[
∑
G(s1) · · ·G(sr)
(
Qs1,...,sr〈a1···ar〉
)Ki
]
×
∏
j





∑
{vja}
[
∑
Qˆs1,...,sr〈a1···ar〉〈γj(s1, a1, h) · · · γj(sr, ar, h)〉h
]Cj





〉
K,C,Λ
(B15)
where
DQDQˆ ≡
(
∏ dQdQˆ
2πi/N
)
, (B16)
and the summations run over all the allowed values of r, 〈a1 · · · ar〉 and s1, . . . sr.
Under the assumption of replica symmetry in the form
Qs1,...,sr〈a1···ar〉 = Q0〈x
r〉x, (B17)
Qˆs1,...,sr〈a1···ar〉 = Qˆ0〈xˆ
r〉xˆ, (B18)
where the averages over x and xˆ are taken with respect to the field distributions π(x) and πˆ(xˆ) respectively, we can
show by straightforward algebraic manipulations that
∑
Qs1,...,sr〈a1···ar〉Qˆ
s1,...,sr
〈a1···ar〉 = Q0Qˆ0〈[1 + (q − 1)xxˆ]
n〉x,xˆ, (B19)
∑
G(s1) · · ·G(sr)
(
Qs1,...,sr〈a1···ar〉
)Ki
= QKi0
〈{
1 +
[
∑
s
G(s)
]
Ki
∏
l=1
xl
}n〉
x
, (B20)
where it is easy to see that
∑
s
G(s) = ∆(q)− 1 = q − 1, (B21)
and
∑
{vja}
[
∑
Qˆs1,...,sr〈a1···ar〉〈γj(s1, a1, h) · · ·γj(sr, ar, h)〉h
]Cj
=
QˆCj0
〈



q−1
∑
v=0
Cj
∏
l=1
[1 + ω(v, hl)xˆl]



n
〉
xˆ,h
,
(B22)
with
ω(v, hl) ≡
q−1
∑
s=1
exp
[
i2πsq (hl ⊗ v)
]
=
{
q − 1, ifhl ⊗ v = 0,
−1, otherwise. (B23)

14
We can simplify the last equation by noting that
q−1
∑
v=0
Cj
∏
l=1
[1 + ω(v, hl)xˆl] =
Cj
∏
l=1
[1 + (q − 1)xˆl] + (q − 1)
Cj
∏
l=1
(1− xˆl). (B24)
Let us write
Zn =
∫
DQDQˆ eNs˜, (B25)
with
s˜ = − 1N lnN
(2)
A − nλ ln q −Q0Qˆ0〈[1 + (q − 1)xxˆ]
n〉x,xˆ +
1
N lnΦ, (B26)
where
Φ =
〈
NΛ
Λ!
QΛ0 QˆΛ0
∏
i
〈[
1 + (q − 1)
Ki
∏
l=1
xl
]n〉
x
×
∏
j
〈



Cj
∏
l=1
[1 + (q − 1)xˆl] + (q − 1)
Cj
∏
l=1
(1− xˆl)



n
〉
xˆ
〉
K,C,Λ
(B27)
Let us define α ≡ NQ0Qˆ0. For n ≪ 1, we can consider only the leading contributions in the number of replicas,
which gives
lnΦ = ln ǫ(α) + nǫ(α)
∑
i
〈
αΛ
Λ!
〈
ln
[
1 + (q − 1)
Ki
∏
l=1
xl
]〉
x
〉
K,C,Λ
n
ǫ(α)
∑
j
〈
αΛ
Λ!
〈
ln



Cj
∏
l=1
[1 + (q − 1)xˆl] + (q − 1)
Cj
∏
l=1
(1− xˆl)



〉
xˆ
〉
K,C,Λ
,
(B28)
with
ǫ(α) =
〈αΛ
Λ!
〉
K,C,Λ
. (B29)
Substituting the above formulas in s˜ for n → 0, the extremization with respect to Q0, Qˆ0, π(x) and πˆ(xˆ) leads to
the saddle point equations (15), (16) and (17).
APPENDIX C: NUMBER OF MATRICES
Here we give the detailed calculation of the average number of GF (q) (M) × N matrices for large N and N .
Repeating the formula given in section V, we have
NA =
∑
{Aij}


M
∏
i=1
δ


N
∑
j=1
χ(Aij),Ki






N
∏
j=1
δ
(
M
∑
i=1
χ(Aij), Cj
)

. (C1)

15
with χ(Aij) = 0 if Aij = 0 and 1 otherwise. Following a similar procedure as in B, we use the integral representations
of the Kroenecker delta functions to write it as
NA =
∮
DWDZ
∏
i,j
∑
Aij
(WiZj)χ(Aij)
=
∮
DWDZ
∏
i,j
[1 + (q − 1)WiZj]
=
∮
DWDZ
∏
i

1 +
N
∑
r=1
(q − 1)rW ri
∑
〈j1···jr〉
Zj1 · · ·Zjr


=
∮
DWDZ

1 +
M
∑
s=1
∑
〈i1···is〉
∑
r1,...,rs
(q − 1)r1+···+rsW r1i1 · · ·W
rs
is F (r1, Z) · · ·F (rs, Z)

,
(C2)
where
F (r, Z) ≡
∑
〈j1···jr〉
Zj1 · · ·Zjr . (C3)
The integrals over the W ’s can pass through the summations and will factorize to give the corresponding Kroenecker
delta functions resulting in
NA = (q − 1)
P
i Ki
∮
DZF (K1, Z) · · ·F (KM , Z)
= (q − 1)Λ
∮
DZF (K1, Z) · · ·F (KM , Z)
= (q − 1)Λ
∮
DZ
∏
i
∑
〈j1···jKi〉
Zj1 · · ·ZjKi
≈ (q − 1)Λ
∮
DZ
∏
i
1
Ki!


N
∑
j=1
Zj


Ki
= (q − 1)Λ
∮
DZ 1∏
i Ki!


N
∑
j=1
Zj


P
i Ki
=
(q − 1)Λ
∏
i Ki!
∮
DZ
∑
j1,...,jΛ
Zj1 · · ·ZjΛ
=
(q − 1)Λ
∏
i Ki!
(
Λ
C1
)(
Λ− C1
C2
)
· · ·
(
Λ− C1 − · · · − CN−1
CN
)
,
(C4)
which gives the final result
NA =
(q − 1)ΛΛ!
∏
i Ki!
∏
j Cj !
. (C5)
APPENDIX D: PROOF OF REPLICA SYMMETRY
Using the fact that the random matrices can be seen as statistical physics systems with Hamiltonian H(v) ≡
N − ln δ(Av, 0) we now prove that this implies that the replica symmetric solution is the exact one. In fact, the form
of the Hamiltonian implies that
P(v) =
[
∑
v
δ(Av, 0)
]−1
= q−d(A). (D1)

16
The distribution of the overlaps of the spins is given by
P(ρ) =
〈
δ

ρ− 1N
N
∑
j=1
σjσ′j


〉
σ,σ′
= q−2d(A)
∑
v,v′
δ(Av, 0)δ(Av′, 0)δ

ρ− 1N
N
∑
j=1
exp
(
2πi
q
(
vj + v′j
)
)

.
(D2)
Let us call
g(v,v′) ≡ δ

ρ− 1N
N
∑
j=1
exp
(
2πi
q
(
vj + v′j
)
)

, (D3)
and note that g(v,v′) = g(0,v ⊕ v′). Therefore we can write
P(ρ) = q−2d(A)
∑
v,v′
δ(Av, 0)δ(Av′, 0)g(0,v⊕ v′)
= q−2d(A)
∑
v,v′
δ(Av, 0)δ(Av′, 0)
∑
u
δ(u,v ⊕ v′)g(0,u)
= q−2d(A)
∑
u
g(0,u)
[
∑
v
δ(Av, 0)
∑
v′
δ(Av′, 0)δ(u,v ⊕ v′)
]
= q−2d(A)
∑
u
g(0,u)
[
∑
v
δ(Av, 0)δ(A(u⊕ (−v)), 0)
]
= q−d(A)
∑
u
δ(Au, 0)g(0,u)
=
〈
δ

ρ− 1N
N
∑
j=1
σj


〉
σ
.
(D4)
Therefore, the distribution of the overlaps is the same as the distribution of the magnetization in the spin systems.
This implies that there is no spin glass phase in the system and, therefore, no replica symmetry breaking [9]. The
above calculation can also be viewed as a consequence of the gauge invariance of the Hamiltonian with respect to the
transformation v → v ⊕ v′, where Av′ = 0, which leads basically to the same calculation above.
[1] R. McEliece, Theory of Information & Coding (Cambridge University Press, Cambridge, MA, 2002 2nd edition).
[2] K. T. Phelps, J. Rifa`, and M. Villanueva, IEEE Trans. Inf. Theory 51, 3931 (2005).
[3] M. Davey and D. MacKay, IEEE Communications Letters 2, 165 (1998).
[4] K. Nakamura, Y. Kabashima, and D. Saad, Eurphys. Lett. 56, 610 (2001).
[5] C. Cooper, Random Structures and Algorithms 16, 209 (2000).
[6] J. Blo¨mer, R. Karp, and E. Weiz, Random Structures and Algorithms 10, 407 (1998).
[7] X. Feng and Z. Zhang, Applied Mathematics and Computation 185, 689 (2007).
[8] M. Me´zard, G. Parisi, and M. Virasoro, Spin Glass Theory and Beyond (World Scientific Publishing Co., Singapore, 1987).
[9] H. Nishimori, Statistical Physics of Spin Glasses and Information Processing (Oxford University Press, Oxford, UK, 2001).
[10] Y. Kabashima and D. Saad, J. Phys. A. 37, R1 (2004).
[11] R. Monasson and R. Zecchina, Phys. Rev. Lett. 76, 3881 (1996).
[12] J. van Mourik and D. Saad, Phys. Rev. E 66, 056120 (2002).
[13] R. Mulet, A. Pagnani, M. Weigt, and R. Zecchina, Phys. Rev. Lett. 89, 268701 (2002).
[14] R. C. Alamino and D. Saad, J. Phys A: Math. Theor. 40, 12259 (2007).
[15] T. Tanaka and D. Saad, Technical report (unpublished).
[16] B.-Y. Wang and F. Zhang, Discrete Mathematics 187, 211 (1998).
[17] R. Vicente, D. Saad, and Y. Kabashima, Europhys. Lett. 51, 698 (2000).
[18] I. S. Gradshteyn and I. M. Ryzhik, in Table of Integrals, Series, and Products, edited by A. Jeffrey and D. Zwillinger
(Academic Press, USA, 1993).