Journal of Machine Learning Research 6 (2005) 695–709

Submitted 11/04; Revised 3/05; Published 4/05

Estimation of Non-Normalized Statistical Models
by Score Matching
Aapo Hyvärinen

aapo.hyvarinen@helsinki.fi

Helsinki Institute for Information Technology (BRU)
Department of Computer Science
FIN-00014 University of Helsinki, Finland

Editor: Peter Dayan

Abstract
One often wants to estimate statistical models where the probability density function is
known only up to a multiplicative normalization constant. Typically, one then has to resort
to Markov Chain Monte Carlo methods, or approximations of the normalization constant.
Here, we propose that such models can be estimated by minimizing the expected squared
distance between the gradient of the log-density given by the model and the gradient of
the log-density of the observed data. While the estimation of the gradient of log-density
function is, in principle, a very difficult non-parametric problem, we prove a surprising
result that gives a simple formula for this objective function. The density function of theobserved data does not appear in this formula, which simplifies to a sample average of asum of some derivatives of the log-density given by the model. The validity of the method
is demonstrated on multivariate Gaussian and independent component analysis models,
and by estimating an overcomplete filter set for natural image data.
Keywords: statistical estimation, non-normalized densities, pseudo-likelihood, Markov
chain Monte Carlo, contrastive divergence

1. Introduction
In many cases, probabilistic models in machine learning, statistics, or signal processing are
given in the form of non-normalized probability densities. That is, the model contains an
unknown normalization constant whose computation is too difficult for practical purposes.
Assume we observe a random vector x ∈ Rn which has a probability density function(pdf) denoted by px (.). We have a parametrized density model p(.; θ), where θ is an mdimensional vector of parameters. We want to estimate the parameter θ from x, i.e. wewant to approximate px (.) by p(.; θ̂) for the estimated parameter value θ̂. (We shall here
consider the case of continuous-valued variables only.)
The problem we consider here is that we only are able to compute the pdf given by themodel up to a multiplicative constant Z(θ):
p(ξ; θ) =1
q(ξ; θ).
Z(θ)

That is, we do know the functional form of q as an analytical expression (or any form thatcan be easily computed), but we do not know how to easily compute Z which is given by
c 2005 Aapo Hyvärinen.

Hyvärinen

an integral that is often analytically intractable:Zq(ξ; θ) dξ.
Z(θ) =
ξ∈Rn

In higher dimensions (in fact, for almost any n > 2), the numerical computation of thisintegral is practically impossible as well.
Usually, estimation of non-normalized models is approached by Markov Chain Monte
Carlo (MCMC) methods, which are very slow, or by making some approximations, which
may be quite poor (Mackay, 2003).
Non-normalized models are often encountered in continous-valued Markov random fields,
which are widely used in image modelling, see e.g. (Bouman and Sauer, 1993; Li, 2001).
In general, undirected graphical models cannot be normalized except in the Gaussian case.
Other recent work in image modelling also includes non-normalized models (Hyvärinen andHoyer, 2001; Teh et al., 2003). Presumably, the number of useful applications for nonnormalized models is much larger than the present literature suggests. Non-normalized
models have been avoided because their estimation has been considered too difficult; theadvent of efficient estimation methods may significantly increase their utility.
In this paper, we propose a simple method for estimating such non-normalized models.
This is based on minimizing the expected squared distance of the score function of x andthe score function given by the model. (By score function, we mean here the gradientof log-density.) We show that this distance can be estimated by a very simple formula
involving only sample averages of some derivatives of the logarithm of the pdf given by themodel. Thus, the computations involved are essentially not more complicated than in thecase where we know an analytical expression for the normalization constant. The proposed
formula is exact and does not involve any approximations, which is why we are able toprove the local consistency of the resulting method. Minimization of the proposed objectivefunction thus provides an estimation method that is computationally simple yet statistically
locally consistent.

2. Estimation by Score Matching
In the following, we use extensively the gradient of the log-density with respect to the data
vector. For simplicity, we call this the score function, although according the conventional
definition, it is actually the score function with respect to a hypothetical location parameter(Schervish, 1995). For the model density, we denote the score function by ψ(ξ; θ): ∂ log p(ξ;θ)ψ1 (ξ; θ)
∂ξ1

 

..
..
=
ψ(ξ; θ) = 
 = ∇ξ log p(ξ; θ).
.
.
∂ log p(ξ;θ)ψn (ξ; θ)∂ξnThe point in using the score function is that it does not depend on Z(θ). In fact we
obviously haveψ(ξ; θ) = ∇ξ log q(ξ; θ).
(1)
Likewise, we denote by ψ x (.) = ∇ξ log px (.) the score function of the distribution of observed
data x. This could in principle be estimated by computing the gradient of the logarithm of696

Estimation by Score Matchinga non-parametric estimate of the pdf—but we will see below that no such computation isnecessary. Note that score functions are mappings from Rn to Rn .
We now propose that the model is estimated by minimizing the expected squared distance between the model score function ψ(.; θ) and the data score function ψ x (.). We define
this squared distance asZ
1px (ξ)kψ(ξ; θ) − ψ x (ξ)k2 dξ.
(2)J(θ) =
2 ξ∈RnThus, our score matching estimator of θ is given byθ̂ = arg min J(θ).
θ

The motivation for this estimator is that the score function can be directly computedfrom q as in (1), and we do not need to compute Z. However, this may still seem to be avery difficult way of estimating θ, since we might have to compute an estimator of the datascore function ψ x from the observed sample, which is basically a non-parametric estimation
problem. However, no such non-parametric estimation is needed. This is because we canuse a simple trick of partial integration to compute the objective function very easily, as
shown by the following theorem:
Theorem 1 Assume that the model score function ψ(ξ; θ) is differentiable, as well as someweak regularity conditions.1
Then, the objective function J in (2) can be expressed as

n 
X
1
2
px (ξ)∂i ψi (ξ; θ) + ψi (ξ; θ) dξ + const.
J(θ) =
2
ξ∈Rn
Z

(3)

i=1

where the constant does not depend on θ,
ψi (ξ; θ) =∂ log q(ξ; θ)∂ξiis the i-th element of the model score function, and∂i ψi (ξ; θ) =∂ψi (ξ; θ)∂ 2 log q(ξ; θ)
=
∂ξi
∂ξi2

is the partial derivative of the i-th element of the model score function with respect to thei-th variable.
The proof, given in the Appendix, is based a simple trick of partial integration that has
previously been used in the theory of independent component analysis for modelling thedensities of the independent components (Pham and Garrat, 1997).
We have thus proven the remarkable fact that the squared distance of the model scorefunction from the data score function can be computed as a simple expectation of certain
1. Namely: the data pdf px (ξ) is differentiable, the expectations Ex {kψ(x; θ)k2 } and Ex {kψ x (x)k2 } are
finite for any θ, and px (ξ)ψ(ξ; θ) goes to zero for any θ when kξk → ∞.

697

Hyvärinenfunctions of the non-normalized model pdf. If we have an analytical expression for the
non-normalized density function q, these functions are readily obtained by derivation using(1) and taking further derivatives.
In practice, we have T observations of the random vector x, denoted by x(1), . . . , x(T ).
The sample version of J is obviously obtained from (3) as
1
˜
J(θ)
=
T

T X
n 
X
t=1 i=1
1
2∂i ψi (x(t); θ) + ψi (x(t); θ) + const.
2(4)

which is asymptotically equivalent to J due to the law of large numbers. We propose to
estimate the model by minimization of J˜ in the case of a real, finite sample.
One may wonder whether it is enough to minimize J to estimate the model, or whether
the distance of the score functions can be zero for different parameter values. Obviously, if
the model is degenerate in the sense that two different values of θ give the same pdf, wecannot estimate θ. If we assume that the model is not degenerate, and that q > 0 always,
we have local consistency as shown by the following theorem and the corollary:
Theorem 2 Assume the pdf of x follows the model: px (.) = p(.; θ ∗ ) for some θ ∗ . Assume
further that no other parameter value gives a pdf that is equal2 to p(.; θ ∗ ), and that q(ξ; θ) >0 for all ξ, θ. Then
J(θ) = 0 ⇔ θ = θ ∗ .
For a proof, see the Appendix.
Corollary 3 Under the assumptions of the preceding Theorems, the score matching estimator obtained by minimization of J˜ is consistent, i.e. it converges in probability towards
the true value of θ when sample size approaches infinity, assuming that the optimization
algorithm is able to find the global minimum.
The corollary is proven by applying the law of large numbers. As sample size approaches
infinity, J˜ converges to J (in probability). Thus, the estimator converges to a point where
J is globally minimized. By Theorem 2, the global minimum is unique and found at thetrue parameter value (obviously, J cannot be negative).
This result of consistency assumes that the global minimum of J˜ is found by the optimization algorithm used in the estimation. In practice, this may not be true, in particular
because there may be several local minima. Then, the consistency is of local nature, i.e.,
the estimator is consistent if the optimization iteration is started sufficiently close to thetrue value. Note that consistency implies asymptotic unbiasedness.

3. Examples
Here, we provide three simulations to illustrate how score matching works, as well as toconfirm its consistency and applicability to real data.
2. In this theorem and its proof, equalities of pdf’s are to be taken in the sense of equal almost everywherewith respect to the Lebesgue measure.

698

Estimation by Score Matching

3.1 Multivariate Gaussian Density
As a very simple illustrative example, we consider estimation of the parameters of themultivariate Gaussian density.
3.1.1 Estimation
The probability density function is given byp(x; M, µ) =11exp(− (x − µ)T M(x − µ)),Z(M, µ)2where M is a symmetric positive-definite matrix (the inverse of the covariance matrix). Of
course, the expression for Z is well-known in this case, but this serves as an illustrationof the method. As long as there is no chance of confusion, we use x here as the generaln-dimensional vector. Thus, here we have1
q(x) = exp(− (x − µ)T M(x − µ)),2

(5)and we obtain
ψ(x; M, µ) = −M(x − µ),and∂i ψ(x; M, µ) = −mii .
Thus, we obtain
1
˜
J(M,
µ) =
T

T X
X
1
[
−mii + (x(t) − µ)T MM(x(t) − µ)].
2
t=1

(6)

i

To minimize this with respect to µ, it is enough to compute the gradient1∇µ J˜ = MMµ − MM
T

T
X

x(t),

t=1

P
which is obviously zero if and only if µ is the sample average T1 Tt=1 x(t). This is truly aminimum because the matrix MM that defines the quadratic form is positive-definite.
Next, we compute the gradient with respect to M, which gives∇M J˜ = −I + M

T

T

t=1

t=1

1 X
1 X
[ (x(t) − µ)(x(t) − µ)T ]M,(x(t) − µ)(x(t) − µ)T +2T
2TPwhich is zero if and only if M is the inverse of the sample covariance matrix T1 Tt=1 (x(t) −µ)(x(t) − µ)T , which thus gives the score matching estimate.
Interestingly, we see that score matching gives exactly the same estimator as maximum
likelihood estimation. In fact, the estimators are identical for any sample (and not just
asymptotically). The maximum likelihood estimator is known to be consistent, so the scorematching estimator is consistent as well.
699

Hyvärinen

3.1.2 Intuitive Interpretation
This example also gives some intuitive insight into the principle of score matching. Let us
consider what happened if we just maximized the non-normalized log-likelihood, i.e., log ofq in (5). It is maximized when the scale parameters in M are zero, i.e., the model variances
are infinite and the pdf is completely flat. This is because then the model assigns the sameprobability to all possible values of x(t), which is equal to 1. In fact, the same applies tothe second term in (6), which thus seems to be closely connected to maximization of the
non-normalized log-likelihood.
Therefore, the first term in (3) and (6), involving second derivatives of thePlogarithmof q, seems to act as a kind of a normalization term. Here it is equal to − i mii . To
minimize this, the mii should be made as large (and positive) as possible. Thus, this term
has the opposite effect to the second term. Since the first term is linear and the second
term polynomial in M, the minimum of the sum is different from zero.
A similar interpretation applies to the general non-Gaussian case. The second term in(3), expectation of the norm of score function, is closely related to maximization of nonnormalized likelihood: if the norm of this gradient is zero, then in fact the data point is ina local extremum of the non-normalized log-likelihood. The first term then measures whatkind of an extremum this is. If it is a minimum, the first term is positive and the value of J
is increased. To minimize J, the first term should be negative, in which case the extremum
is a maximum. In fact, the extremum should be as steep a maximum (as opposed to aflat maximum) as possible to minimize J. This counteracts, again, the tendency to assign
the same probability to all data points that is often inherent in the maximization of the
non-normalized likelihood.
3.2 Estimation of Basic Independent Component Analysis Model
Next, we show the validity of score matching in estimating the following modellog p(x) =n
X

G(wkT x) + Z(w1 , . . . , wn ),

(7)

k=1

which is the basic form of the independent component analysis (ICA) model. Again, thenormalization constant is well-known and equal to − log | det W| where the matrix W has
the vectors wi as rows, but this serves as an illustration of our method.
The nice thing about this model is that we can easily generate data that follows thismodel. In fact, if latent variables si , i = 1 . . . , n are independently distributed and have the
pdf given by exp(G(si )), the linear transformationx = As(8)

with A = W−1 follows the pdf’s given in (7), see e.g. (Hyvärinen et al., 2001). Thus, wewill be estimating the generative model in (8) using the non-normalized likelihood in (7).
Here, we choose the distribution of the components si to be so-called logistic withπG(s) = −2 log cosh( √ s) − log 4.
2 3700

Estimation by Score MatchingThis distribution is normalized to unit variance as typical in the theory of ICA. The scorefunction of the model in (7 is given byψ(x; W) =nX

wk g(wkT x),

(9)

k=1

where the scalar nonlinear function g is given byππg(s) = − tanh( √ s).
3
2 3
The relevant derivatives of the score function are given by:∂i ψi (x) =

n
X

2 0
wki
g (wkT x),

k=1

and the sample version of the objective function J˜ is given by


T X
n
n
n
n
X
X
X
X
1
1
2 0

J˜ =
wji g(wjT x(t))
wki g(wkT x(t))
wki
g (wkT x(t)) +T2
t=1 i=1

=

n
X
k=1

j=1

k=1

kwk k2

k=1

T

n

T

t=1

j,k=1

t=1

1 X T
1X
1X 0 Tg (wk x(t)) +wj wk
g(wkT x(t))g(wjT x(t)). (10)T
2T

We performed simulations to validate the consistency of score matching estimation, andto compare its efficiency with respect to maximum likelihood estimation. We generated
data following the model as described above, where the dimension was chosen to be n = 4.
Score matching estimation consisted of minimizing J˜ in (10) by a simple gradient descent;likelihood was maximized using a natural gradient method (Amari et al., 1996; Hyvärinen
et al., 2001), using the true value of Z. We repeated the estimation for several different
sample sizes: 500, 1000, 2000, 4000, 8000, and 16000. For each sample size, the estimation
was repeated 11 times using different random initial points in the optimization, and different
random data sets. For each estimate, a measure of asymptotic variance was computed asfollows. The matrix ŴA, where Ŵ is the estimate was normalized row-by-row so that thelargest value on each row had an absolute value of 1. Then, the sum of squares of all theelements was computed, and 4 (i.e. the sum of the squares of the four elements equal toone) was subtracted. This gives a measure of the squared error of the estimate (we cannot
simply compare ŴA with identity because the order of the components is not well-defined).
For each sample size and estimator type (score matching vs. maximum likelihood) we thencomputed the median error.
Figure 1 shows the results. The error of score matching seems to go to zero, which
validates the theoretical consistency result of Theorem 2. Score matching gives slightly
larger errors than maximum likelihood, which is to be expected because of the efficiency
results of maximum likelihood estimation (Pham and Garrat, 1997).
In the preceding simulation, we knew exactly the proper function g to be used in thescore function. To investigate the robustness of the method to misspecification of the score
701

Hyvärinen

-0.4
-0.6
-0.8
-1
-1.2
-1.4
-1.6
-1.8
-2
-2.2
-2.4
2.6

2.8

3

3.2

3.4

3.6

3.8

4

4.2

4.4

Figure 1: The estimation errors of score matching (solid line) compared with errors of maximum likelihood estimation (dashed line) for the basic ICA model. Horizontal
axis: log10 of sample size. Vertical axis: log10 of estimation error.

function (a well-known problem in ICA estimation), we ran the same estimation methods,score matching and maximum likelihood, for data that was generated by a slightly different
distribution. Specifically, we generated the data so that the independent components s i
had Laplacian distributions of unit variance (Hyvärinen et al., 2001). We then estimatedthe model using exactly the same g as before, which was not theoretically correct. The
estimation errors are shown in Figure 2. We see that score matching still seems consistent.
Interestingly, it now performs slightly better than maximum likelihood estimation (which
would more properly be called quasi-maximum likelihood estimation due to the misspecification (Pham and Garrat, 1997)).
3.3 Estimation of an Overcomplete Model for Image Data
Finally, we show image analysis results using an overcomplete version of the ICA model.
The likelihood is defined almost as in (7), but the number of components m is larger than thedimension of the data n, see e.g. (Teh et al., 2003), and we introduce some extra parameters.
The likelihood is given bylog p(x) =mX

αk G(wkT x) + Z(w1 , . . . , wn , α1 , . . . , αn ),

(11)

k=1

where the vectors wk = (wk1 , . . . , wkn ) are constrained to unit norm (unlike in the preceding
example), and the αk are scaling parameters. We introduce here the extra parameters α k toaccount for different distributions for different projections. Constraining α k = 1 and m = n
702

Estimation by Score Matching

-1.4
-1.6
-1.8
-2
-2.2
-2.4
-2.6
-2.8
-3
-3.2
-3.4
2.6

2.8

3

3.2

3.4

3.6

3.8

4

4.2

4.4

Figure 2: The estimation errors of score matching compared with errors of maximum likelihood estimation for the basic ICA model. This time, the pdf of the independent
components was slightly misspecified. Legend as in Fig. 1.

and allowing the wk to have any norm, this becomes the basic ICA model of the preceding
subsection.
The model is related to ICA with overcomplete bases (Hyvärinen et al., 2001; Hyvärinen
and Inki, 2002; Olshausen and Field, 1997), i.e. the case where there are more independent
components and basis vectors than observed variables. In contrast to most ICA models,the overcompleteness is expressed as overcompleteness of filters wk which seems to makethe problem a bit simpler because no latent variables need to be inferred. However, thenormalization constant Z is not known when G is non-quadratic, i.e. when the model is
non-Gaussian, which is why previous research had to resort to MCMC methods (Teh et al.,2003) or some approximations (Hyvärinen and Inki, 2002).
We have the score function
ψ(x; W, α1 , . . . , αm ) =mX

αk wk g(wkT x),

k=1

where g is the first derivative of G. Going through similar developments as in the case of
the basic ICA model, the sample version of the objective function J˜ can be shown to equal

J˜ =

m
X
k=1

T

m

T

t=1

j,k=1

t=1

1X 0 T
1 X
1X
αkg (wk x(t)) +αj αk wjT wk
g(wkT x(t))g(wjT x(t)).
T
2
T
703

(12)

Hyvärinen

Figure 3: The overcomplete set of filters wi estimated from natural image data. Notethat no dimension reduction was performed, and we show filters instead of basis
vectors, which is why the results are much less smooth and “beautiful” than somepublished ICA results (Hyvärinen et al., 2001).

We estimated the model for image patches of 8 × 8 pixels taken from natural images,
see P.O. Hoyer’s imageica package.3 As preprocessing, the DC component (i.e. the meangray-scale value) was removed from each image patch, reducing the effective dimensionalityof the data to n = 63. The data was also whitened, i.e. the model was used in a linearlytransformed space (the exact method of whitening has no significance). We set m = 200.
We also took the tanh function as g, which corresponds to G(u) = log cosh(u) (we did not
bother to find the right scaling as in the basic ICA case). The objective function J˜ in (12)was optimized by gradient descent. The wi were set to random initial values, and the αiwere all set to the initial value 1.5 that was found to be close to the optimal value in pilot
experiments.
The obtained vectors wi are shown in Figure 3. For the purposes of visualization, thevectors were converted back to the original space from the whitened space. The optimal α iwere in the range 0.5 . . . 2.
To show that the method correctly found different vectors and not duplicates of a smaller
set of vectors, we computed the dot-products between the vectors, and for each w i , weselected the largest absolute value of dot-product |wiT wj |, j 6= i. The dot-products werecomputed in the whitened space. The histogram of these maximal dot-products is shownin Figure 4. They are all much smaller than 1 (in absolute value), in fact all are smaller
than 0.5. Since the vectors wi were normalized to unit norm, this shows that no two wi
were close to equal, and we did find m different vectors.

4. Discussion
Here we discuss the connections of our method to two well-known methods before concluding
the paper.
3. The package can be downloaded at http://www.cs.helsinki.fi/patrik.hoyer/.

704

Estimation by Score Matching

25

20

15

10

5

0
0.25

0.3

0.35

0.4

0.45

0.5

Figure 4: The distribution of maximal dot-products of a filter wi with all other filters,
computed in the whitened space.

4.1 Comparison with Pseudo-Likelihood Estimation
A related method for estimating non-normalized models is maximization of pseudo-likelihood(Besag, 1974). The idea is to maximize the product of marginal conditional likelihoods. The
pdf is approximated bylog ppseudo (x) =

n
X
i=1

p(xi |x1 , . . . , xi−1 , xi+1 , . . . , xn ),(13)

and the likelihood is computed using this approximation. The idea was originally developedin connection with Markov random fields, in which context it is quite natural because
the conditional probabilities are often given as part of the model specification. The idea
can still be used in the general case considered in this article. However, the conditional
probabilities in (13) are not necessarily readily available and need to be computed. In
particular, these conditional densities need to be normalized. The computational burden
needed in the normalization is reduced from the original problem since we only need tonumerically compute n one-dimensional integrals which is far more feasible than a single
n-dimensional integral. However, compared to score matching, this is a computationallyexpensive method since score matching avoids the need for numerical integration altogether.
The question of consistency of pseudo-likelihood estimation seems to be unclear. Someconsistency proofs were provided by Besag (1974, 1977), but these only apply to special
cases such as Gaussian or binary random fields. Sufficiently general consistency results on
pseudo-likelihood estimation seem to be lacking. This is another disadvantage with respect
to score matching, which was shown above to be (locally) consistent.
705

Hyvärinen

4.2 Comparison with Contrastive Divergence
An interesting approximative MCMC method called contrastive divergence was recently
proposed by Hinton (2002). The basic principle is to use an MCMC method for computing
the derivative of the logarithm of the normalization factor Z, but the MCMC is allowed to
run for only a single iteration (or a few iterations) before doing the gradient step.
The method is generally biased, even asymptotically (Carreira-Perpiñán and Hinton,
2005b), except in some special cases such as the multivariate Gaussian distribution (CarreiraPerpiñán and Hinton, 2005a). Score matching is thus preferable if a consistent estimator iswanted.
The computational efficiency of contrastive divergence is difficult to evaluate since it
is not really a single method but a family of methods, depending on the MCMC method
used. For the case of continuous-valued variables that we consider here, a Metropolis-type
algorithm would probably be the method of choice, but there is a large number of different
variants whose performances are likely to be quite different.
Nevertheless, contrastive divergence is a much more general method than score matching
since it is applicable to intractable latent variable models. It can also handle binary/discretevariables—in fact, it is probably much easier to implement, using Gibbs sampling, for binary
variables than for continous-valued variables. Extension of score matching to these two cases
is an important problem for future research.
4.3 Conclusion
We have proposed a new method, score matching, to estimate statistical models in the casewhere the normalization constant is unknown. Although the estimation of the score function
is computationally difficult, we showed that the distance of data and model score functions
is very easy to compute. The main assumptions in the method are: 1) all the variables
are continuous-valued and defined over Rn , 2) the model pdf is smooth enough. Score
matching provides a computationally simple yet locally consistent alternative to existing
methods, such as MCMC and various approximative methods.

Acknowledgments
I am grateful to Patrik Hoyer, Jarmo Hurri, and Shohei Shimizu for comments on the
manuscript, to Sam Roweis for interesting discussions, and to Miguel Carreira-Perpiñán
and Geoffrey Hinton for providing access to unpublished results. The work was supportedby the Academy of Finland, Academy Research Fellow position and project #48593.

Appendix A. Proof of Theorem 1
Definition (2) gives
J(θ) =

Z


1
1
2
2
T
px (ξ) kψ x (ξ)k + kψ(ξ; θ)k − ψ x (ξ) ψ(ξ; θ) dξ.
2
2


(14)

(For simplicity, we omit the integration domain in here.) The first term in brackets does
not depend on θ, and can be ignored. The integral of the second term is simply the integral
706

Estimation by Score Matchingof the sum of the second terms in brackets in (3). Thus, the difficult thing to prove is that
integral of the third term in brackets in (14) equals the integral of the sum of the first terms
in brackets in (3). This term equals
XZ−
px (ξ)ψx,i (ξ)ψi (ξ; θ)dξ,iwhere ψx,i (ξ) denotes the i-th element of the vector ψ x (ξ). We can consider the integral
for a single i separately, which equalsZZZpx (ξ) ∂px (ξ)∂px (ξ)∂ log px (ξ)ψi (ξ; θ)dξ = −ψi (ξ; θ)dξ = −ψi (ξ; θ)dξ.
− px (ξ)∂ξipx (ξ) ∂ξi∂ξi
The basic trick of partial integration needed the proof is simple: for any one-dimensional
pdf p and any function f , we haveZZZ
Z
p0 (x)
p(x)(log p)0 (x)f (x)dx = p(x)f (x)dx = p0 (x)f (x)dx = − p(x)f 0 (x)dxp(x)
under some regularity assumptions that will be dealt with below.
To proceed with the proof, we need to use a multivariate version of such partial integration:
Lemma 4lim

a→∞,b→−∞

f (a, ξ2 , . . . , ξn )g(a, ξ2 , . . . , ξn ) − f (b, ξ2 , . . . , ξn )g(b, ξ2 , . . . , ξn )Z ∞Z ∞∂f (ξ)
∂g(ξ)
dξ1 +
g(ξ)
dξ1 ,
=
f (ξ)
∂ξ1
∂ξ1
−∞
−∞

assuming that f and g are differentiable. The same applies for all indices of ξ i , but fornotational simplicity we only write the case i = 1 here.
Proof of lemma:∂g(ξ)∂f (ξ)∂f (ξ)g(ξ)= f (ξ)+ g(ξ)
.
∂ξ1
∂ξ1
∂ξ1

We can now consider this as a function of ξ1 alone, all other variables being fixed. Then,integrating over ξ1 ∈ R, we have proven the lemma.
Now, we can apply this lemma on px and ψ1 (ξ; θ) which were both assumed to bedifferentiable in the theorem, and we obtain:
−

Z


Z Z
∂px (ξ)∂px (ξ)ψ1 (ξ; θ)dξ = −ψ1 (ξ; θ)dξ1 d(ξ2 , . . . , ξn )
∂ξ1
∂ξ1
Z 
lim
[px (a, ξ2 , . . . , ξn )ψ1 (a, ξ2 , . . . , ξn ; θ)=−a→∞,b→−∞− px (b, ξ2 , . . . , ξn )ψ1 (b, ξ2 , . . . , ξn ; θ)]Z∂ψ1 (ξ; θ)−px (ξ)dξ1 d(ξ2 , . . . , ξn ).
∂ξ1707

Hyvärinen

For notational simplicity, we consider the case of i = 1 only, but this is true for any i.
The limit in the above expression is zero for any ξ2 , . . . , ξn , θ because we assumed thatpx (ξ)ψ(ξ; θ) goes to zero at infinity. Thus, we have proven thatZZ
∂px (ξ)∂ψi (ξ; θ)−
ψi (ξ; θ)dξ =px (ξ)dξ,∂ξi∂ξithat is, integral of the the third term in brackets in (14) equals the integral of the sum of
the first terms in brackets in (3), and the proof of the theorem is complete.

Appendix B. Proof of Theorem 2
Assume J(θ) = 0. Then, the assumption q > 0 implies px (ξ) > 0 for all ξ, which implies
that ψ x (.) and ψ(.; θ) are equal. This implies log px (.) = log p(.; θ) + c for some constantc. But c is necessarily 0 because both px and p(.; θ) are pdf’s. Thus, px = p(.; θ). By
assumption, only θ = θ ∗ fulfills this equality, so necessarily θ = θ ∗ , and we have proven theimplication from left to right. The converse is trivial.

References
S.-I. Amari, A. Cichocki, and H. H. Yang. A new learning algorithm for blind source
separation. In Advances in Neural Information Processing Systems 8, pages 757–763.
MIT Press, 1996.
J. Besag. Spatial interaction and the statistical analysis of lattice systems. Journal of theRoyal Statistical Society, Series B, 36(2):192–236, 1974.
J. Besag. Efficiency of pseudolikelihood estimation for simple gaussian fields. Biometrika,
64(3):616–618, 1977.
C. Bouman and K. Sauer. A generalized gaussian image model for edge-preserving MAP
estimation. IEEE Transactions on Image Processing, 2(3):296–310, 1993.
M. Á. Carreira-Perpiñán and G. E. Hinton. On contrastive divergence (CD) learning.
Technical report, Dept of Computer Science, University of Toronto, 2005a. In preparation.
M. Á. Carreira-Perpiñán and G. E. Hinton. On contrastive divergence learning. In Proceedings of the Workshop on Artificial Intelligence and Statistics (AISTATS2005), Barbados,
2005b.
G. E. Hinton. Training products of experts by minimizing contrastive divergence. NeuralComputation, 14(8):1771–1800, 2002.
A. Hyvärinen and P. O. Hoyer. A two-layer sparse coding model learns simple and complex
cell receptive fields and topography from natural images. Vision Research, 41(18):2413–
2423, 2001.
A. Hyvärinen and M. Inki. Estimating overcomplete independent component bases fromimage windows. Journal of Mathematical Imaging and Vision, 17:139–152, 2002.
708

Estimation by Score Matching

A. Hyvärinen, J. Karhunen, and E. Oja. Independent Component Analysis. Wiley Interscience, 2001.
S. Z. Li. Markov Random Field Modeling in Image Analysis. Springer, 2nd edition, 2001.
D. J. C. Mackay. Information Theory, Inference and Learning Algorithms. Cambridge
University Press, 2003.
B. A. Olshausen and D. J. Field. Sparse coding with an overcomplete basis set: A strategy
employed by V1? Vision Research, 37:3311–3325, 1997.
D.-T. Pham and P. Garrat. Blind separation of mixture of independent sources througha quasi-maximum likelihood approach. IEEE Transactions on Signal Processing, 45(7):1712–1725, 1997.
M. Schervish. Theory of Statistics. Springer, 1995.
Y. W. Teh, M. Welling, S. Osindero, and G. E. Hinton. Energy-based models for sparse
overcomplete representations. Journal of Machine Learning Research, 4:1235–1260, 2003.

709