Image super-resolution with PCA Reduced generalized
Gaussian mixture models
Dang-Phuong-Lan Nguyen, Johannes Hertrich, Jean-Francois Aujol, Yannick
Berthoumieu
To cite this version:
Dang-Phuong-Lan Nguyen, Johannes Hertrich, Jean-Francois Aujol, Yannick Berthoumieu. Image
super-resolution with PCA Reduced generalized Gaussian mixture models. 2022. hal-03664839
HAL Id: hal-03664839
https://hal.science/hal-03664839
Preprint submitted on 11 May 2022
HAL is a multi-disciplinary open access
archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.
L’archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d’enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
Image super-resolution with PCA Reduced generalized
Gaussian mixture models
Dang-Phuong-Lan Nguyen
∗†
Johannes Hertrich
Yannick Berthoumieu
‡
Jean-Francois Aujol
∗
†
May 11, 2022
Abstract
Single Image Super-Resolution algorithms based on patches have been noticed and widely
used over the past decade. Recently, generalized Gaussian mixture models (GGMMs) have been
shown to be a suitable tool for many image processing problems because of the flexible shape
parameter. In this work, we first propose to use a joint GGMM learned from concatenated
vectors of high- and low-resolution training patches. For each low-resolution patch, we compute
the minimum mean square error (MMSE) estimator and generate the high-resolution image
by averaging these estimates. We select the MMSE approach using GGMM as the method
is invariant to affine contrast change and also invariant to a linear super-resolution operator.
Unfortunately, the large dimension of the concatenated high- and low-resolution patches leads
to instabilities and an intractable computational effort when estimating the parameters of the
GGMM. Thus, we propose to combine a GGMM with a principal component analysis and
derive an EM algorithm for estimating the parameters of the arising model. We demonstrate
the performance of our model by numerical examples on synthetic and real images of materials’
microstructures.
Key words: Generalized Gaussian mixture model, image super-resolution, high-dimensional
data, dimensionality reduction.
1
Introduction
Recent developments in imaging techniques deeply modify the way materials science consider their
objects of research. The motivation to develop new super-resolution (SR) methods originates in
these modifications. The material image has moved from micron-scale resolution to true nano scale.
Studying material data is very challenging because high-resolution (HR) and low-resolution (LR)
images have different contrasts. Besides the knowledge of super-resolution operator is not known. In
this article, we introduce a novel SR method that deals with the unknown operator and the problem
of an affine contrast change.
Single-image super-resolution (SISR) aims to reconstruct a high-resolution (HR) image from a
single low-resolution (LR) input image. In general, the LR image is generated by
XL = AXH + ϵ
where XL is the low-resolution image, XH is the high-resolution image, and A is an unknown
super-resolution operator. Image restoration techniques can be grouped into two main approaches:
pixel-based and patch-based. Patch-based methods are a popular and successful class of image
restoration techniques, especially in super-resolution.
∗ Univ.
Bordeaux, Bordeaux INP, CNRS, IMB, UMR 5251, F-33400 Talence, France, {lan.nguyen, jeanfrancois.aujol}@math.u-bordeaux.fr
† Univ. Bordeaux, Bordeaux INP, CNRS, IMS, UMR 5218, F-33400 Talence, France, yannick.berthoumieu@imsbordeaux.fr
‡ TU Berlin, Straße des 17. Juni 136, D-10587 Berlin, Germany, {j.hertrich}@math.tu-berlin.de.
1
Among the patch-based methods for super-resolution images, Zoran and Weiss [1] proposed
a method called expected patch log likelihood (EPLL) which can be applied for several inverse
problems. This method uses the negative log-likelihood function of a Gaussian mixture model
(GMM) as a regularizer of the inverse problem. The high-resolution image XH is estimated by
solving the following:
X
2
log p (XH,i )
argmin ∥AXH − XL ∥ − λ
XH
i∈I
where p is the probability density function of the GMM and (XH,i )i∈I are the patches in the highresolution image XH . The following years, S. Parameswaran et al. extended the EPLL method by
various accelerations in [2] Fast EPLL (called FEPLL).
Before the reconstruction step, these approaches must learn the parameters of the GMM model
from the HR training patches. GMM can be learned by the expectation–maximization algorithm
(EM) [3, 4]. In the past decade various works [5] [6] [7] have shown that generalized Gaussian
mixture model (GGMM) has more flexibility to adapt to the shape of data probability density
function and less sensibility for over fitting the number of classes than the GMM. In [8], C. Deledalle
et al. proposed a method that uses the generalized Gaussian mixture model (GGMM) for the EPLL
algorithm (EPLL-GGMM). They showed that the GGMM gets the distribution of patches better
than a GMM and that it outperforms the GMM when used in the EPLL framework. However,
the EPLL-GGMM model learns the parameters of the mixture model by estimating the covariance
matrix with the following formula:
PN
T
i=1 ξi,k XH,i XH,i
Σk =
.
PN
i=1 ξi,k
N
Here {XH,i }i=1 is a set of all HR patches and ξi,k is the conditional distribution of the component
k given the patch XH,i . Unfortunately, this formula is not true in the GGMM setting.
In [6], Z. Boukouvalas et al. introduced a Riemannian averaged Fixed-point (RA-FP) algorithm
for multivariate generalized Gaussian distribution (MGGD) parameter estimation. The RA-FP
algorithm can effectively estimate the scatter matrix for any value of the shape parameter: for
instance, β = 1 is the Gaussian distribution and β = 0.5 is the Laplacian distribution. In the case
of the generalized Gaussian mixture model, we generalize the algorithm from [9] for the weighted
maximum likelihood estimation using the EM algorithm based on a fixed-point algorithm. This
algorithm estimates the parameters of the mixture model, including the shape parameters of each
component. This algorithm is called FP-EM algorithm for generalized Gaussian mixture model.
Furthermore, in some cases, the HR and LR images of the real data may not have the same
contrast. Both EPLL and FEPLL methods for GMM or GGMM models require knowledge of the
operator A of the super-resolution problem and the contrast change parameters, which we do not
know in practice. To avoid this requirement, P. Sandeep et al. [10] proposed a method for the SR
problem using a joint GMM, which is invariant to the affine change of contrast. It describes the
distribution of the concatenated vectors of high and low resolution patches sampled from a large
database of pairs of high resolution images and the corresponding low resolution one. Each highresolution patch is reconstructed from the low-resolution one by using the Minimum Mean Squared
Error (MMSE) estimator as:
X̂H = argmin MSE (T (XL ) − XL ) .
T
In our work, we provide a method that uses the MMSE estimator for the super-resolution based
on GGMM for joint HR-LR modeling, called MMSE-GGMM. This method adapts the FP-EM
algorithm to learn the GGMM models. Our method not only avoids learning the operator A, but
is also invariant by affine contrast change. This is demonstrated through Theorem 4 as well as the
experimental results in Section 6.
However, the EM algorithm for these models (GMMs or GGMMs) becomes very slow because
the number of data points becomes large and three-dimensional images. To overcome these issues,
we propose a dimensionality reduction method.
2
The standard method for dimensionality reduction is the principal component analysis (PCA)
[11], which was introduced in 1901 by Karl Pearson. Based on this method, [12] introduced the
probabilistic principal component analysis (PPCA) for the Gaussian distribution. It was extended
to the Gaussian Mixture model in [13] with the mixture of probabilistic principal component analysers (MPPCA) model. Following these works, [14, 15] proposed directly a GMM with PCA-based
constrained covariance matrices. J. Hertrich et al. proposed a PCA-GMM model in [16], which
combines the GMM with a PCA by adding the minimization term of the PCA and the negative
log-likelihood function of the GMM on the reduced dimensionality data points. The authors rewrote
this minimization problem again as a negative log-likelihood function of a Gaussian mixture model,
which has additional constraints on the parameters.
To develop the PCA with the GMM model up to the GGMM model, we introduce the PCAGGMM model. The PCA-GGMM model incorporates the PCA model with the negative log likelihood function of the GGMM on the low-dimensional data point. We derive an EM algorithm
with a special M-step for finding a minimizer of our objective function. The M-step requires solutions of maximization problems with constraints on the Stiefel manifold. Due to the complexity
of the GGMMs and the non-convexity of the PCA-GGMM, we use the ”Half-Quadratic splitting”
[17] method to tackle the minimization problem in the M-step of the EM algorithm. Finally, the
PCA-GGMM model can be used for the super-resolution model based on the MMSE method for
GGMM.
The main contribution of the paper is to use the GGMM model for SR. We also show that the
MMSE method for the GGMM model has the affine contrast-invariant property. And finally, this
paper adapts the PCA method of dimensionality reduction to the GGMM model.
The paper is organized as follows. Section 2 briefly recalls the super-resolution problem introduced by [10] using joint mixture models. We revisit the EM algorithm for mixture models and
describe the FP-EM algorithm for generalized Gaussian mixture model in Section 3. In Section
4, the PCA-reduced generalized Gaussian distribution is introduced. We derive an algorithm for
estimating its parameters. In Section 5, the MMSE Estimator for generalized Gaussian distribution
and the affine transformation of contrast and brightness are explained. Finally, Section 6 shows
numerical examples of super-resolution based on our models on 2D and 3D images. The paper
concludes with a discussion in Section 7.
2
Super-Resolution via Joint Mixture Models
In this section, we briefly revisit the super-resolution by [10] using joint mixture models. Let
{Pθ : θ ∈ Θ} be a parametric family of probability distributions with density functions p(·|θ).
Then, a mixture model is a probability distribution defined by the probability density function
F (x|w, Θ) =
K
X
wk p(x|θk ),
w = (w1 , ..., wk ), Θ = (θ1 , ..., θk ),
(1)
k=1
where the weights wk are non-negative and sum up to one. If θ = (µ, Σ) and Pθ = N (µ, Σ), we
speak of a Gaussian mixture model (GMM).
Throughout this paper, we aim to reconstruct the unknown high-resolution image XH based on
a low resolution observation XL . Here, we assume that we have given a reference image, where the
high-resolution image x̃H as well as the low-resolution observation x̃L are both known.
In this setting, Sandeep and Jacob proposed in [10] to reconstruct XH using the following steps.
2
1. Learning a joint mixture model. In a first step, for given low-resolution patches x̃L,i ∈ Rτ
2 2
of an image and their high-resolution counterparts x̃H,i ∈ Rq τ , q ∈ N, q ≥ 2, i = 1, ..., N , we
x̃
H,i
∈ Rn ,
approximate the distribution of the joint high- and low-resolution patches xi =
x̃L,i
n = (q 2 +1)τ 2 by a mixture model. Here, Sandeep and Jacob use GMMs such that the resulting
approximation is characterized by the parameters
w = (wk )k ,
µ = (µk )k ,
3
Σ = (Σk )k
with
µk =
µ
H,k
,
µL,k
Σ=
Σ
H,k
ΣT
HL,k
ΣHL,k
.
ΣL,k
2. Estimation of the high-resolution patches using the MMSE estimator. In the second
2
step, we estimate the high resolution correspondence of a given low resolution patch xL ∈ Rτ .
For this purpose, we first select the component k∗, such that the likelihood that xL belongs
to the k ∗ -th component is maximal, i.e., we compute
k ∗ = argmax wk p(xL |µL,k , ΣL,k )
k=1,...,K
2 2
Afterwards, we estimate the high resolution patch xH ∈ Rq τ as the minimum mean square
estimator (MMSE) of XH given XL = xL for a random variable X = (XH , XL ) ∼ Pθk . For
Gaussian distributions, this corresponds to
x̂H = µH,k + ΣHL,k Σ−1
L (xL − µL,k ).
We give a more detailed explanation on MMSE estimators in Section 5.
3. Reconstruction of the high-resolution image by patch averaging. Finally, we estimate
for any patch within the low resolution image, the high resolution correspondence. Afterwards,
qτ,qτ
we reconstruct the high resolution image as follows: Let xH = (xk,l )qτ
be a twok,l=1 ∈ R
dimensional high-resolution patch. Then, we assign to each pixel xk,l the weight
ρk,l := exp − γ2 (k − qτ2+1 )2 + (l − qτ2+1 )2 .
After that, we add up for each pixel in the high resolution image the corresponding weighted
pixel values and normalize the result by dividing by the sum of the weights.
Contributions of this paper
Sandeep and Jacob considered in [10] Gaussian mixture models. However, it is well-known that
GMMs are very flat-tailed and therefore lack within their ability to approximate patch distributions. Furthermore, for large magnification factors or three-dimensional images, the dimension of
the joint mixture model becomes very high. Unfortunately, this leads to numerical instabilities and
an intractable computational effort for estimating the joint mixture model. Thus, we propose within
this paper two improvements of Sandeeps and Jacobs super-resolution method [10]. First, we replace
the Gaussian distribution by generalized Gaussian distributions. Both contain a shape parameter,
which can adjust to more heavy-tailed or flat-tailed distributions. Second, we incorporate a dimensionality reduction within mixture models of elliptical distributions and derive an algorithm for
estimating its parameters. Additionally to these methodical improvements, we apply our method to
three-dimensional real-world images showing material microstructures.
3
Parameter Estimation of Mixture Models
In this section, we focus on the parameter estimation of mixture models. For this purpose, we aim
to employ the expectation-maximization (EM) algorithm. First, in Subsection 3.1, we revisit the
generic form of the EM algorithm for mixture models. This algorithm contains as a subproblem
the weighted maximum likelihood estimation for the corresponding parametric distribution family.
Thus, we consider this subproblem in Subsection 3.2 for some elliptical distributions, namely, the
Gaussian distribution, the generalized Gaussian distribution.
4
3.1
EM algorithm for Mixture Models
In the following, we consider the EM algorithm to estimate the parameters of the mixture models.
EM algorithms were introduced in [18] and can be used for various probability distributions. We
refer to [3] for an overview. Given samples x1 , ..., xN aims to minimize the negative log-likelihood
function
N
1 X
L(w, Θ) = −
log(F (xi |w, Θ),
N i=1
where F is defined as in (1). Then, the EM algorithm for mixture models reads as Algorithm 1, see
e.g. [19].
Algorithm 1 EM Algorithm for Mixture Models
Input: x = (x1 , ..., xN ) ∈ Rn×N , initial estimate w(0) , Θ(0) .
for r = 0, 1, ... do
E-Step: For k = 1, ..., K and i = 1, . . . , N compute
(r)
(r)
w p(xi |θ )
(r)
αi,k = PK k (r) k (r)
j=1 wj p(xi |θj )
M-Step: For k = 1, ..., K compute
N
1 X (r)
α ,
N i=1 i,k
(r+1)
=
(r+1)
= argmax
wk
θk
θk
N
nX
i=1
o
(r)
αi,k log(p(xi |θk )) .
end for
As outlined in the algorithm, we have to compute the weighted maximum likelihood estimator
as a subproblem of the EM algorithm, i.e., we need to maximize
ℓ(θ) =
N
X
αi log(p(xi |θ))
(2)
i=1
for the underlying parametric distribution family. For readability, in the following article, we denote
αi,k by αi and θk by θ. In the following sections, we show how the EM algorithm can be done for
various distributions.
3.2
Weighted Maximum Likelihood Estimation of generalized Gaussian
Distribution
In order to solve subproblem (2) within the EM algorithm for mixture models, we consider in
the following the weighted maximum likelihood estimation for several elliptical distributions. Let
SPD(n) be the set of positive-definite symmetric matrices. A probability distribution En (µ, Σ, g)
with µ ∈ Rn and Σ ∈ SPD(n) is called an elliptical distribution if it admits a density function
p(x|µ, Σ, g) ∝ Σ−1/2 g((x − µ)T Σ−1 (x − µ)),
where g : R≥0 → R≥0 is Lebesgue measurable with
Z ∞
tn/2−1 g(t)dt < ∞.
0
In this paper, we consider Gaussian, generalized Gaussian as examples of elliptical distributions.
5
Gaussian distribution. A probability distribution is a Gaussian distribution N (µ, Σ) with mean
µ ∈ Rn and covariance matrix Σ ∈ SPD(n) if it has the probability density function
1
p(x|µ, Σ) =
exp − 12 (x − µ)T Σ−1 (x − µ) .
(2π)n/2 |Σ|1/2
The Gaussian distribution is an elliptical distribution with g(t) = exp (− 21 t). For samples x1 , ..., xN
with weights α1 , ..., αN , we can calculate the weighted maximum likelihood estimator of µ and Σ,
i.e., the solution of (2) by setting the gradients of the objective function to zero. This leads to
estimators
1
µ̂ = PN
i=1 αi
N
X
αi xi ,
i=1
1
Σ̂ = PN
i=1 αi
N
X
i=1
N
X
1
αi (xi − µ̂)(xi − µ̂)T = PN
i=1 αi
T
αi xi xT
i − µ̂µ̂ .
i=1
Generalized Gaussian distribution. We call a probability distribution a generalized Gaussian
distribution GG(µ, Σ, β) with location µ ∈ Rn , scatter matrix Σ ∈ SPD(n) and shape parameter
β ∈ R>0 , if it has the probability density function
p(x|µ, Σ, β) =
C(β)
exp −
|Σ|1/2
1
2
(x − µ)T Σ−1 (x − µ)
β
,
where the normalizing constant C is given by
C(β) =
βΓ( n2 )
n .
)
2n/(2β) π n/2 Γ( 2β
Here, Γ denotes the gamma function. For small values of the shape parameter β, we obtain heavytailed distributions, while large values of β mean that the corresponding generalized Gaussian distribution is flat-tailed. The generalized Gaussian distribution includes the Gaussian distribution
for β = 1 and the Laplace distribution for β = 1/2. The generalized Gaussian distribution is an
elliptical distribution with g(t) = exp (− 21 tβ ), β > 0.
In the literature, several algorithms have been proposed to estimate the parameters of a generalized Gaussian distribution of samples x1 , ..., xN . Most of these methods are based on fixed point
(FP) iterations [7, 20] or first-order computations on the Riemannian manifold of parameters [6].
For this work, we generalize the algorithm from [9] for weighted maximum likelihood estimation (2),
i.e., for maximizing
N
X
ℓ(µ, Σ, β) =
αi log(p(xi |µ, Σ, β)),
i=1
for some weights α1 , ..., αN > 0. Here, setting the gradient to zero yields that any minimizer (µ, Σ, β)
of ℓ fulfills
PN
αi δiβ−1 xi
,
µ = Pi=1
N
β−1
i=1 αi δi
PN
αi δiβ−1 (xi − µ)(xi − µ)T
,
Σ = i=1
PN
i=1 αi
β = β + ρ∇β ℓ(µ, Σ, β),
−1
where δi = (xi −µ)Σ (xi −µ) and ρ > 0. This motivates us to maximize ℓ by generating a sequence
(µ(r) , Σ(r) , β (r) ) by fixed-point (FP) iteration.
PN
(r) (r)
αi (δi )β −1 xi
(r+1)
,
µ
= Pi=1
(r) β (r) −1
N
i=1 αi (δi )
PN
(r) (r)
αi (δi )β −1 (xi − µ(r) )(xi − µ(r) )T
,
Σ(r+1) = i=1
PN
i=1 αi
β (r+1) = β (r) + ρ∇β ℓ(µ(r) , Σ(r) , β (r) ),
6
(r)
where δi
= (xi − µ(r) )(Σ(r) )−1 (xi − µ(r) ). Note that ∇β ℓ(µ, Σ, β) can be explicitly derived as
∇β ℓ(µ, Σ, β) =
1
β
+
N
N
X
X
n
n
)
+
log
2
α
+
αi δiβ log(δi ).
ψ(
i
2β
2β 2
i=1
i=1
The EM algorithm with the FP iteration for GGMM is called FP-EM algorithm. For the sake of
completeness, we add the detailed computations in Appendix A.
Remark. In the unweighted case, i.e., α1 = . . . = αN , the authors of [9] show convergence of the
above iteration.
Remark. Note that for the estimation of β, several other approaches could be considered such as
the Newton-Raphson method [21].
4
4.1
Dimensionality-Reduced Generalized Gaussian Distribution
Combining PCA with generalized Gaussian distribution
When considering high-dimensional probability distributions, the EM algorithm 1 for mixture models
can be numerically unstable and computationally costly. Thus, we aim to reduce the dimension of the
corresponding probability distributions. For this purpose, we aim to combine principal component
analysis (PCA) with generalized Gaussian distributions. The PCA finds the affine subspace {U t+b :
t ∈ Rd }, 1 ≤ d ≪ n with the smallest squared distance from the samples x1 , ..., xN by minimizing
P (U, b) =
N
X
∥(U U T − I)(xi − b)∥2 ,
i=1
for b ∈ Rn and U ∈ St(d, n), where St(d, n) := {U ∈ Rn,d : U T U = I} is the Stiefel manifold.
In the following, we consider the samples x1 , ..., xN ∈ Rn , which are approximately located
in a d-dimensional subspace of Rn . Then, our objective is to find simultaneously the subspace
{U t + b : t ∈ Rd }, U ∈ St(d, n), b ∈ Rn that contains the samples xi and a generalized Gaussian
distribution GG(µ, Σ, β) with density function q(·|µ, Σ, β). For this, we propose to minimize the
function
ℓPCA (U, b, µ, Σ, β) :=
N
X
1
∥(U U T − I)(xi − b)∥2 − log(q(U T (xi − b)|µ, Σ, β)),
2
2σ
i=1
(3)
which is a weighted sum of a PCA term and the negative log-likelihood function of q within the
lower dimensional subspace.
The following proposition gives a general property of the function lP CA for the generalized
Gaussian distribution. The proof is given in Appendix B.
Proposition 1. Let GG(µ, Σ, β) be a generalized Gaussian distribution with a density function
q(·|µ, Σ, β). The function ℓPCA (U, b, µ, Σ, β) is up to a constant the negative log-likelihood function
of a probability density
p(x|U, b, µ, Σ, β) :=
1
(2πσ 2 )(n−d)/2
exp(− 2σ1 2 ∥(U U T − I)(x − b)∥2 )q(U T (xi − b)|µ, Σ, β).
The distribution Pθ with θ = (U, b, µ, Σ, β) corresponding to this density is given by
Pθ = r# (GG(µ, Σ, β) ⊗ N (0, σ 2 In−d )),
r(x) = V x + b,
where V = (U |Ũ ) is an orthogonal matrix, ⊗ denotes the product measure and r# P is the pushforward measure of P under r.
7
Due to its relationship to the PCA, we call the probability distribution p(·|U, b, µ, Σ, β) a PCAreduced generalized Gaussian distribution.
Remark (PCA-reduced Gaussian distributions). The authors of [16] used mixture models of PCAreduced Gaussian distributions. In particular, they proved that a PCA-reduced Gaussian distribution is again a Gaussian distribution with a special structured covariance matrix. Using this
result, they derived an EM algorithm for PCA-reduced Gaussian mixture models, which showed
a significant speed-up compared to a high-dimensional GMM. However, for generalized Gaussian
distribution, the PCA-reduced counterpart is not a generalized Gaussian distribution such that the
analysis of [16] is no longer applicable.
In the following subsection, we introduce a method for estimating the parameters of mixture
models for the PCA-reduced generalized Gaussian distribution. We call this mixture model a PCAGGMM model.
4.2
Weighted Maximum Likelihood Estimation
To estimate the parameters of a PCA-reduced generalized Gaussian distribution from samples
x1 , ..., xN with weights α1 , ..., αN , our objective is to minimize the weighted negative log likelihood
function
N
X
αi log(p(xi |U, b, µ, Σ, β)).
i=1
By Proposition 1 this is equivalent to minimizing the function
ℓPCA,α (U, b, µ, Σ, β) :=
N
X
i=1
αi
1
T
2
T
∥(U
U
−
I)(x
−
b)∥
−
log(q(U
(x
−
b)|µ,
Σ,
β))
,
i
i
2σ 2
which is up to the weights α the same function as ℓPCA from (3). Unfortunately, we cannot minimize
ℓPCA,α directly. Instead, we apply a technique called half-quadratic splitting [17]. That is, instead
of considering ℓPCA,α (U, b, µ, Σ, β), we minimize for some large η > 0 the function
Hσ,η (z, U, b, µ, Σ, β) :=
N
X
i=1
αi
1
T
2
T
2
∥(U
U
−
I)(x
−
b)∥
+
η∥U
(x
−
b)
−
z
∥
−
log(q(z
|µ,
Σ,
β))
,
i
i
i
i
2σ 2
i.e., we replace U T (xi − b) in the second summand by zi and penalize the squared distance of
U T (xi − b) and zi . Thus, for η → ∞, the minimizers of Hσ,η converge to the minimizers of ℓPCA,α .
Now, we minimize Hσ,η iteratively with respect to (U, b), z and (µ, Σ, β), i.e., we generate a sequence
(U (r) , b(r) , µ(r) , Σ(r) , β (r) )r by iteration
U
(r)
∈ argmin
N
X
U ∈St(d,n) i=1
b(r) ∈ argmin
b∈Rn
z (r) ∈
N
X
αi
i=1
argmin
1
(r−1) 2
T
(r−1) 2
T
(r−1)
∥(U
U
−
I)(x
−
b
)∥
+
η∥U
(x
−
b
)
−
z
∥
,
i
i
i
2σ 2
(4)
1
(r−1) 2
∥(U (r) (U (r) )T − I)(xi − b)∥2 + η∥(U (r) )T (xi − b) − zi
∥ ,
2
2σ
(5)
αi
N
X
z=(z1 ,...,zN )∈Rd,N i=1
αi η∥(U (r) )T (xi − b(r) ) − zi ∥2 − log(q(zi |µ(r−1) , Σ(r−1) , β (r−1) )) (6)
and finally
(µ(r) , Σ(r) , β (r) ) ∈
argmax
N
X
µ∈Rd ,Σ∈SPD(d),β>0 i=1
(r)
αi log(q(zi |µ, Σ, β))
(7)
The final step (7) is the weighted maximum likelihood estimation of the generalized Gaussian distribution q, which was discussed in the previous section. It remains to solve the first three steps.
8
Solving (4) by Uzawas’ Algorithm.
Using that for any y ∈ Rn and U ∈ St(d, n) it holds that
∥(U U T − I)y∥2 = y T (U U T − I)2 y = y T (U U T U U T − 2U U T + I)y = y T (U U T − 2U U T + I)y
= y T (I − U U T )y = ∥y∥2 − ∥U T y∥2 ,
the optimization problem in (4) reads as
N
X
1
αi − 2 ∥U T (xi − b)∥2 + η∥U T (xi − b) − zi ∥2
2σ
U ∈St(d,n) i=1
argmin
Now, we transform the problem in an unconstrained problem by considering the Lagrangian function
L(U, Λ) =
N
X
1
αi − 2 ∥U T (xi − b)∥2 + η∥U T (xi − b) − zi ∥2 + tr(Λ(U T U − I)).
2σ
i=1
Then, the solution of (4) is equivalent to solving the saddle-point problem
argmin max L(U, Λ).
d,d
U ∈Rn,d Λ∈R
To solve this saddle-point problem, we use the Uzawas’ algorithm [22], which consists of the following
two steps: First, we minimize L with respect to U . Second, we perform a gradient ascent step with
respect to Λ.
The second step is given by
Λ(r+1) = Λ(r) + ρ((U (r+1) )T U (r+1) − I)
In the case of our specific Lagrangian L, the second step is given by the following proposition. The
proof is given in Appendix C.
Proposition 2. Let Λ be fixed. Then, any minimizer of L(U, Λ) solves the Sylvester equation
(η −
1
2σ 2 )
N
X
i=1
N
X
αi (xi − b)(xi − b)T U + U Λ = η
αi (xi − b)ziT .
i=1
Setting the gradient of the objective function to zero leads to
Solving (5).
N
X
i=1
αi
1
T
T
(b − xi ) − ηU zi = 0
I − U U + ηU U
2σ 2
which is equivalent to
b=
PN
i=1
αi xi +
1
2σ 2
−1 PN
η i=1 αi U zi
I − U U T + ηU U T
.
PN
i=1 αi
Solving (6) using a Gradient Descent.
decouples into
(r)
zi
Due to the sum-structure of problem (6), the solution
∈ argmin η∥(U (r) )T (xi − b(r) ) − zi ∥2 − log(q(zi |µ(r−1) , Σ(r−1) , β (r−1) )).
(8)
zi ∈Rd
This is now a d-dimensional optimization problem, which (owing to its low dimension) can be
solved efficiently via a gradient descent scheme.
9
Remark (Differentiability of the objective). In the case of generalized Gaussian distributions, the
function g is given by g(x) = exp(− 12 xβ ), which is differentiable at R>0 . Thus, the objective function
in (8) is differentiable for any zi ̸= µ(r−1) . In the case β ≥ 1, the derivative of g can be continuously
extended to 0. In this case, the objective function in (8) is differentiable.
The reduction in dimensionality of the PCA-GGMM model reduces the execution time of the
E-step, but the M-step is slower than the GGMM due to its higher complexity. However, all
computations in the M-step are on d-dimensional optimization problems, and the FP iteration is
implemented on d-dimensional data. This limits the numerical instability problem to estimate the
generalized Gaussian mixture model of the original data using the FP-EM algorithm. Moreover, the
PCA-reduced generalized Gaussian distribution can be generalized to the PCA-reduced elliptical
distribution by combining PCA with the elliptical distribution En (µ, Σ, g) with the density function
q(·|µ, Σ, g). We denote p(·|U, b, µ, Σ, g) as a probability density function of the PCA-reduced elliptical
distribution. By using the half-quadratic splitting technique, we can apply the iterations (4), (5),
(6) and (7) with q(·|µ, Σ, g) instead of q(·|µ, Σ, β) to estimate the parameters of a PCA-reduced
elliptical distribution p(·|U, b, µ, Σ, g).
5
MMSE Estimator for Generalized Gaussian Distributions
In the following, we consider the computation of a high-resolution patch xH given the corresponding
low-resolution patch xL and a mixture model (1), where the distribution Pθk is a generalized Gaussian
distribution. In particular, we have that the parameters are given by θk = (µk , Σk , βk ) for some
location µk ∈ Rn , a scatter matrix Σk ∈ SPD(n) and β > 0.
As in [10], we first split the parameters
µk =
µH,k
µL,k
,
Σk =
Σ
H,k
ΣT
HL,k
ΣHL,k
ΣL,k
into components belonging to the high- and low-resolution part of the distribution p(·|θk ). Then we
select the component k ∗ , which best matches the low-resolution patch xL by
k ∗ := argmax p(xL |µL,k , ΣL,k , βk ).
k=1,...,K
Afterwards, we compute xH as the minimum mean square estimator (MMSE) of XH given
XL = xL , where X = (XH , XL ) ∼ Pθk .
Remark. Given a random variable XL : Ω → Rd in a probability space (Ω, A, P), we wish to estimate
a random variable XH : Ω → RD , i.e., we seek an estimator T : Rd → RD such that X̂H = T (XL )
approximates XH . A common quality measure for this task is the mean square error E∥XH −
T (XL )∥22 , which gives rise to the definition of the minimum mean square estimator
TMMSE ∈ argmin E∥XH − T (XL )∥22 .
(9)
T
Under weak additional regularity assumptions on the estimator T , the Lehmann-Scheffé theorem [23]
states that the general solution of the minimization problem (9) is given by
TMMSE (XL ) = E(XH |XL ).
In general, it is not possible to give an analytical expression of the MMSE estimator TMMSE . In
the case of elliptical distributions, the following theorem can be found, e.g., in [24, Theorem 8].
Theorem 3. Assume that X = (XH , XL ) : Ω → Rn has an elliptical distribution Pθ with parameters
θ = (µ, Σ, g), where
Σ
ΣHL
µH
H
, Σ=
µ=
T
µL
ΣHL ΣL
10
Then, for each PXL -almost every xL , we have that the conditional distribution PXH |XL =xL is given
by the elliptical distribution Pθ̂ , where the parameters θ̂ = (µ̂, Σ̂, ĝ) are given by
µ̂ = µH + ΣHL Σ−1
L (xL − µL ),
T
Σ̂ = ΣH − ΣHL Σ−1
L ΣHL ,
ĝ(t) = g(t + t0 )
with t0 = (xL − µL )T Σ−1
L (xL − µL ).
Since a generalized Gaussian distribution is also an elliptical distribution, the MMSE estimator
TM M SE for the generalized Gaussian distribution is expressed as
TMMSE (XL ) = E(XH |XL ) = µH + ΣHL Σ−1
L (xL − µL ).
(10)
In some applications, the high-resolution image and its low-resolution correspondence have different brightness and contrast. The following theorem shows that the MMSE estimator is invariant
under affine transformations of the observation. In particular, the MMSE estimator compensates
for the change in contrast and brightness.
Theorem 4. Assume that X = (XH , XL ) : Ω → Rn has a generalized Gaussian distribution
GG(µ, Σ, β), where
Σ
ΣHL
µH
H
, Σ=
µ=
µL
ΣT
ΣL
HL
Further define the random variable X̃L := α1 XL + α2 with α1 > 0, α2 ∈ R. Then it holds
E(XH |X̃L ) = E(XH |XL ) = µH + ΣHL Σ−1
L (XL − µL ).
Proof. The first equality holds true as it holds
σ(X̃L ) = σ({X̃L−1 (] − ∞, x]) : x ∈ Rd }) = σ({X̃L−1 (] − ∞, α1 x + α2 ]) : x ∈ Rd })
= σ({XL−1 (] − ∞, x]) : x ∈ Rd }) = σ(XL ).
The second equality follows directly from Theorem 3.
Remark. (PCA-reduced generalized Gaussian distribution with the MMSE estimator). In the
PCA-GGMM model, the latent data U T (xi − b) have a d-dimensional generalized Gaussian distribution
GG(µ, Σ, β). Due to the following proposition and the definition of the affine subspace
U t + b : t ∈ Rd , the high-dimensional samples x1 , ..., xN ∈ Rn have a generalized Gaussian distribution GG(b + U µ, U ΣU T , β). Finally, we can apply the MMSE estimator in equation (10) to the
generalized Gaussian distribution GG(b + U µ, U ΣU T , β) of the high dimensional data.
Proposition 5. Let Y ∼ Ed (µ, Σ, g), let B be a n × d matrix, and let b ∈ Rn . Then
b + BY ∼ En (b + Bµ, BΣB T , g).
Proof. Lemma 3.1, p. 5 [25].
6
Numerical Results
In this section, we demonstrate the performance of our models: MMSE-GGMM and PCA-GGMM
on some 2D images such as Gold-hill, Barbara, Camera-man, and our material data. All the implementations were done in Matlab. We run all experiments on PlaFRIM 1 with 32 cores (4 Go/core)
AMD EPYC 7452.
1 https://www.plafrim.fr/
11
6.1
Generation of the test examples
In the experiment, we create the low resolution image from the original image by XL,i = AXH,i + ε
where ε is a white Gaussian noise with standard derivation σ = 0.02. For the operator A, we use
the implementation of [2]:
A = SH
(11)
The blur operator H is given by a convolution with a Gaussian kernel with standard deviation 0.5.
The downsampling operator S : Rm1 ,n1 → Rm2 ,n2 (m1 > m2 and n1 > n2 ) is given by
S=
m2 n2 −1
F
DFm1 ,n1 ,
m1 n1 m2 ,n2
where D : Rm1 ,n1 → Rm2 ,n2 and for x ∈ Rm1 ,n1 the (i, j)-th entry of D (x) is given by
if i ≤ m22 and j ≤ n22 ,
xi,j ,
x
if i > m22 and j ≤ n22 ,
i+m1 −m2 ,j ,
if i ≤ m22 and j > n22 ,
xi,j+n1 −n2 ,
xi+m1 −m2 ,j+n1 −n2 , if i > m22 and j > n22 .
Fm1 ,n1 is the two-dimensional discrete Fourier transform which is defined by Fm1 ,n1 := Fn1 ⊗ Fm1 ,
m1 −1
where Fm1 = (exp (−2πikl/m1 ))k,l=0
.
Figure 1: Top: Images for estimating the mixture models. Bottom: Ground truth for reconstruction.
First column: Material ”FS”, second column: Material ”SiC Diamonds”, third column: goldhill
image.
Training data. For estimating the parameters of the mixture model, we used a part of the
original HR image, such as the top left quarter as the top row in Figure 1 and the corresponding
one of the LR image. From the training part of the LR image, we define a set of overlapping patches
2
N
of size τ × τ : {x̃L,i }i=1 with x̃L,i ∈ Rτ . For the HR image XH , we can also have a set of qτ × qτ
12
2
N
2
HR patches {x̃H,i }i=1 where x̃H,i ∈ Rq τ corresponds to the LR patch x̃L,i . Finally, we define a set
2
2
{x } with a vector x ∈ Rτ (q +1) by concatenating the HR and LR patches as
i
i
xi =
6.2
x̃
H,i
∈ Rn .
x̃L,i
Material data preprocessing
In this part, we discuss material data and the pre-processing step for super-resolution. In the frame
of the ITN MUMMERING, a series of multiscale 3D images have been acquired by synchrotron
microtomography at the SLS beamline TOMCAT. The materials of the two samples were selected
to provide 3D images that have diverse levels of complexity:
• The first one is a sample of Fontainebleau sandstone (”FS”), a natural rock rather homogeneous
and commonly used in the oil industry for flow experiments.
• The second one is a composite (”SiC Diamonds”) obtained by microwave sintering of silicon
and diamonds, see [26].
In Figure 2, the first column shows the HR image, and the second column displays the LR image,
which we acquired from synchrotron microtomography. It is easy to see that, the black region in LR
images is even much more than the HR ones. The important task is to find the position corresponding
to HR image of the subject in the LR images.
For more than a decade, the Scale-invariant feature transform SIFT [27] has arguably been the
most popular keypoint detection and matching method. The SIFT algorithm proposed by David
Lowe consists of a keypoint detection which is independent towards image rotation, scale change,
affine transformation, intensity variation, and viewpoint change. Using the SIFT algorithm, we
determine all the match points between the HR and LR images as shown in Figure 3. Then, the
boundary of the LR image will be found by selecting the minimum and maximum coordinates of
the keypoints. Finally, we obtain the LR image which matches with the HR image. These matched
LR images are given in the third column in Figure 2.
Figure 2: Top: Material ”FS”, bottom: Material ”SiC Diamonds”. First column: Ground-truth
HR image, second column: Input LR image, third column: Match LR images. The red lines in the
second column are the boundaries of the LR image that correspond to the HR image.
13
Figure 3: Matched points between the HR and LR image by using SIFT
6.3
MMSE-GGMM
First, to underline the advantages of the proposed method MMSE-GGMM, we compare our method
with the MMSE-GMM [10] and the EPLL-GMM [1] as well as EPLL-GGMM [8]. In this subsection,
we use all models on standard images of size 512 × 512, such as Gold Hill, Camera-Man, Barbara,
and material images for the magnification factor q = 2. Low-resolution images are created using the
operator A in equation (11). To estimate the parameters of the mixture model, we extract the lowresolution image with the size of the overlap patches τ = 4. Thus, the corresponding high-resolution
patch has size qτ = 8. Then we get N ≈ 1500 patches for the training data to learn the mixture
model. Finally, high-resolution images can be reconstructed from the learned mixture model and LR
images based on the MMSE as described in Section 2. To implement the EPLL method, we estimate
the mixture models using the EM algorithm, which is introduced in [8]. Once again, EPLL requires
the knowledge of the operator A in the image estimation step. For the standard image experiments,
we reconstruct the HR image with the operator A defined by equation (11), called (”EPLL with
given A”). Besides that, we implement the EPLL method with the operator A estimated as in [28]
from the reference HR and LR images, called (”EPLL with learned A”). Since the operator A is
learned from the quarter-known image, it is the same size as the reference image. Thus, we upsample
the operator A by a zero padding technique to have the same size as the whole image.
Table 1 shows the performance of the MMSE and EPLL methods on synthetic images. Both
of these methods are compared based on the GMM, LMM, GGMM with fixed shape parameter β
(e.g. β = 0.8, 0.4, 0.25) and the GGMM model which learns the shape parameter based on FP-EM
algorithm. The PSNR values of MMSE-GGMM with the FP algorithm are slightly higher than
those of the MMSE-GMM and LMM models. To compare our EM algorithm for the GGMM model
with the EM algorithm proposed by Deledalle et al. [8], we compare the MMSE method with the
FP algorithm (MMSE-GGMM FP) and Deledalle’s approach (MMSE-GGMM Deledalle) in Table
1. The MMSE-GGMM FP gives better results than the MMSE-GGMM Deledalle. Furthermore,
we observed that our method can obtain results similar to those of EPLL with given A, while we
do not have to learn the super-resolution operator A. Besides the results of our method are 1.5 dB
higher than EPLL with the learned operator A.
In practice, LR and HR images were acquired under different conditions. Therefore, they may
have different contrasts. Thus, we have to estimate the contrast change parameters between the LR
and HR images. However, Theorem 4 proves that MMSE-GGMM is invariant for an affine contrast
change. Therefore, MMSE-GGMM does not require learning these parameters and the knowledge
of the operator A, while the EPLL method does. To demonstrate these arguments, we perform the
MMSE and the EPLL approach for our material data. Table 2 gives the PSNR value of the HR
reconstructions of real material images using the MMSE and EPLL approach. The second column
14
MMSE
EPLL with given A
EPLL with learned A
GMM
LMM
GGMM β = 0.8
GGMM β = 0.4
GGMM β = 0.25
GGMM FP
GGMM Deledalle
GMM
LMM
GGMM β = 0.8
GGMM β = 0.4
GGMM β = 0.25
GGMM
GMM
LMM
GGMM β = 0.8
GGMM β = 0.4
GGMM β = 0.25
GGMM
Hill
31.60
31.69
31.70
31.68
31.60
31.70
31.52
31.62
31.46
31.58
31.43
31.26
31.58
30.64
30.57
30.56
30.42
30.40
30.67
Camera
32.75
32.82
32.81
32.84
32.80
32.86
32.56
32.91
32.85
32.86
32.83
32.65
32.94
30.74
30.88
30.79
30.76
30.67
30.99
Barbara
25.27
25.31
25.29
25.20
25.29
25.33
25.20
25.39
25.30
25.28
25.30
25.21
25.33
24.32
24.29
24.25
24.31
24.19
24.36
Table 1: PSNRs of the reconstructions of synthetically downsampled 2D images using either MMSE
and EPLL approaches for GMM and LMM, GGMM with different shape parameter β. The magnification factor is q = 2.
group gives the PSNR values of the EPLL approach with a learned operator A, which is estimated
from the known part of the HR image. The PSNR values in the third column group are obtained
from the EPLL method with the learning contrast parameters and the operator A. This operator
is learned as in [28] from the observed LR image and the whole HR image of the ground truth. We
call this approach is ”EPLL with given A”. The PSNRs in Table 2 show that our method achieves
significantly better results than the EPLL method with the learning A from the known HR part.
Our method gives PSNRs 1 dB lower than the EPLL method with learned A, while knowledge of
the operator A is not required in our method.
FS
Sic Diamonds
GMM
33.09
28.00
MMSE
LMM GGMM
33.32
33.35
28.07
28.08
EPLL
GMM
32.25
26.75
with learned A
LMM GGMM
32.31
32.39
26.78
26.82
EPLL with given A
GMM LMM GGMM
33.61
33.86
34.22
29.11
29.29
29.45
Table 2: PSNRs of the reconstructions of material 2D images with contrast change problem using
either MMSE estimator and EPLL approach for GMM, LMM, and GGMM. The magnification factor
is q = 2.
6.4
PCA-GGMM
In the second subsection, we focus on the dimensionality reduction task of our PCA-GGMM model.
We compare our model with the PCA-GMM [16] on 2D artificially images as in the previous subsection and the 2-3D material images: FS”, ”SiC Diamonds” with the zooming ratio q = 2, 4.
As in the previous experiments with MMSE-GGMM, we use the patch size τ = 4 for the low
resolution and qτ = 8 (or qτ = 16) for high resolution images. For the material images with size
2500 × 2500 of HR, this leads to N ≈ 400000 patches for q = 2 and N ≈ 100000 patches for q = 4.
These patches will be used for the training step of the PCA-GGMM model. In this step, we use
15
(a) HR
(b) MMSE-GMM
(c) MMSE-GGMM
(d) LR
(e) EPLL-GMM
(f) EPLL-GGMM
Figure 4: Reconstructions of 2D Sic Diamonds image with magnification factor q = 2 by using
MMSE and EPLL method.
2
2
q
+
1
= 80 or
K = 100 components
for
the
mixture
model.
Each
training
patch
V
of
size
D
=
τ
2
2
D = τ q + 1 = 272 is reduced to d low dimensional for d ∈ {4, 8, 12, 16, 20}. The PSNR values
in Tables 3, 4 and 5 demonstrate that our PCA-GGMM model produces results almost as good as
MMSE-GGMM and even slightly better than MMSE-GMM.
For the 3D material data, we crop a 600 × 600 × 600 image from the material images FS and SiC
Diamonds. Training data are taken from the upper left part of the 3D image 300×300×300. We use
the magnification factor q = 2, the low-resolution patch size τ = 4, and the high-resolution patch
size
qτ = 8. Thus, the dimension of the pair high and low resolution patches is D = τ 3 q 3 + 1 = 576.
These high-dimensional data are reduced to d for d ∈ {20, 40, 60} in the PCA-GGMM model. Table
6 shows the PSNR values of the MMSE method for the GMM, LMM, GGMM model and the
PCA-GMM, PCA-GGMM model.
7
Conclusion
This paper proposed a new algorithm to perform image super-resolution. We extended the image
super-resolution using the GMM method provided by Sandeep and Jacob [10] to the GGMM model,
which is learned by the FP algorithm. We also derived a new method based on a Fixed Point
approach to estimate the parameters of GGMM model. As in the previous work by J.Hertrich
et al. [16] (in the case of GMM), we introduced a new Generalized Gaussian Mixture Model for
dimensionality reduction. The new model, called PCA-GGMM, incorporates the PCA model with
the GGMM model on the low-dimensional data. Experiments with 2D, 3D synthetic, and real
material images demonstrated the effectiveness of the PCA-GGMM models. They showed that
16
MMSE-GMM
MMSE-LMM
MMSE-GGMM
PCA-GMM
PCA-LMM
PCA-GGMM
d
80
80
80
20
16
12
8
4
20
16
12
8
4
20
16
12
8
4
Hill
31.60
31.69
31.70
31.59
31.46
31.23
30.83
30.71
31.61
31.60
31.49
31.09
30.80
31.63
31.60
31.54
31.18
31.02
Camera
32.75
32.82
32.86
32.69
32.61
32.58
32.36
32.19
32.72
32.64
32.61
32.46
32.32
32.74
32.70
32.65
32.53
32.40
Barbara
25.27
25.31
25.33
25.26
25.23
25.11
24.82
24.75
25.29
25.28
25.20
25.02
24.93
25.30
25.28
25.25
25.11
24.94
FS
35.48
35.53
35.57
35.45
35.43
35.46
35.29
34.71
35.51
35.49
35.48
35.36
34.85
35.56
35.52
35.50
35.41
34.97
Diamonds
37.23
37.34
37.39
36.22
36.20
36.15
36.67
35.41
37.32
37.29
37.23
36.88
35.57
37.38
37.35
37.30
36.97
35.71
Table 3: PSNRs of the reconstructions of 2D images using MMSE method with GMM, LMM and
GGMM (with learned shape parameter β) models, and the PCA with PCA-GMM, PCA-LMM,
PCA-GGMM. The magnification factor is set to q = 2.
d
MMSE
FS
20
16
12
8
4
PCA
MMSE
Diamonds
20
16
12
8
4
PCA
GMM
33.09
33.03
32.99
32.86
32.41
32.10
28.00
27.99
27.98
27.85
27.67
27.14
LMM
33.32
33.25
33.21
33.17
32.68
32.43
28.07
28.02
28.00
27.91
27.70
27.21
GGMM
33.35
33.29
33.25
33.20
32.91
32.44
28.08
28.05
28.02
27.94
27.70
27.23
Table 4: PSNRs of the reconstructions of real material 2D image for magnification factor is set
to q = 2 using MMSE method with GMM, LMM and GGMM (with learned shape parameter β)
models, and the PCA with PCA-GMM, PCA-LMM, PCA-GGMM.
17
MMSE-GMM
MMSE-LMM
MMSE-GGMM
PCA-GMM
PCA-LMM
PCA-GGMM
d
20
16
12
8
4
20
16
12
8
4
20
16
12
8
4
Hill
27.15
27.16
27.18
26.89
26.81
26.73
26.41
26.22
27.14
27.05
26.84
26.60
26.46
27.15
27.01
26.83
26.71
26.57
Camera
26.19
26.21
26.23
26.15
26.07
25.96
25.69
25.19
26.18
26.12
26.03
25.87
25.34
26.20
26.14
26.03
25.89
25.40
Barbara
23.64
23.67
23.72
23.47
23.32
23.27
23.01
22.86
23.50
23.41
23.29
23.12
22.94
23.65
23.58
23.43
23.21
23.02
FS
30.72
30.76
30.83
30.73
30.71
30.68
30.47
29.80
30.74
30.71
30.70
30.52
29.84
30.81
30.78
30.74
30.59
29.93
Diamonds
30.74
30.79
30.81
30.73
30.61
30.46
30.18
29.23
30.76
30.68
30.53
30.24
29.41
30.78
30.73
30.59
30.30
29.48
Table 5: PSNRs of the reconstructions of 2D images for magnification factor is set to q = 4 using
MMSE method with GMM, LMM and GGMM (with learned shape parameter β) models, and the
PCA with PCA-GMM, PCA-LMM, PCA-GGMM.
d
MMSE
FS
60
40
20
PCA
MMSE
Diamonds
60
40
20
PCA
GMM
33.34
33.28
33.25
33.17
30.68
30.63
30.55
30.36
LMM
33.36
33.31
33.27
33.19
30.70
30.64
30.57
30.42
GGMM
33.41
33.36
33.35
33.23
30.75
30.67
30.62
30.48
Table 6: PSNRs of the reconstructions of 3D images using MMSE method and PCA-GGMM with
q = 2.
18
our model gives a similar quality for the results as MMSE-GGMM and is better than the PCAGMM model of [16]. Furthermore, our method, which estimates the parameters of the PCA-reduced
generalized model, could be extended to the case of elliptical distributions such as the Student-t
distribution or any other one. Finally, we are aware of the effectiveness of deep learning approaches
for super-resolution, e.g. [29, 30, 31]. In future work, we will consider some deep learning approaches
for super-resolution with high magnification factor and high-dimensional data.
Acknowledgment
Funding by the French Agence Nationale de la Recherche (ANR) under reference ANR-18-CE920050 SUPREMATIM as well as by the German Research Foundation (DFG) at Technical University
of Berlin within the project STE 571/16-1 , is gratefully acknowledged. Phuong-Lan NGUYEN expresses her sincere thanks to Dr. Laurent Facq. He taught her how to use PlaFRIM (a computational
platform) of the University of Bordeaux. The experiments presented in thisarticle were carried out
using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), University
of Bordeaux, Bordeaux INP, and Conseil Régional d’Aquitaine (see https://www.plafrim.fr/).
We also would like to express our deep gratitude to Professor Dominique Bernard and Dr Abdellatif
Saadaldin f rom ICMCB, University of Bordeaux, for the generation and discussion of the material
images in Section 6.
A
Computations for Generalized Gaussian Distributions
In the following, we compute the derivatives of
ℓ(µ, Σ, β) :=
N
X
αi log fβ (xi |µ, Σ),
i=1
which is up to a constant equal to
N
X
αi log
i=1
1
β
Cp (β)
exp − (xi − µ)T Σ−1 (xi − µ)
1/2
2
|Σ|
Using the notation δi = (xi − µ)T Σ−1 (xi − µ), the gradient with respect to µ is given by
∇µ ℓ(µ, Σ, β) =
N
X
αi βΣ−1 (xi − µ)δiβ−1
i=1
Setting the gradient to zero yields
PN
µ = Pi=1
N
αi xi δiβ−1
i=1
Plugging in the formulas
∂aT Σ−1 b
= Σ−T abT Σ−T
∂Σ
and
αi δiβ−1
.
∂ log |Σ−1 |
∂ log |Σ−1 | ∂Σ−1
=
= ΣΣ−2 = Σ−1
∂Σ
∂Σ−1
∂Σ
from [32], we obtain that
∇Σ ℓ(µ, Σ, β) =
N
X
i=1
αi
1 −1
2Σ
−
β β−1 −1
δi Σ (xi − µ)(xi − µ)T Σ−1 .
2
19
Setting the gradient to zero and multiplying by Σ from the left and from the right yields that
PN
αi δiβ−1 (xi − µ)(xi − µ)T
.
Σ = i=1
PN
i=1 αi
Finally, the gradient with respect to β is given by
∇β ℓ(µ, Σ, β) =
N
X
αi ∇β log(Cn (β)) − 12 ∇β δiβ
i=1
Here, we have that ∇β δiβ = δiβ log(δi ) and ∇β log(Cn (β)) =
digamma function. Hence, we obtain
∇β ℓ(µ, Σ, β) =
B
1
β
+
1
β
+
n
2β 2
n
ψ( 2β
) + log 2 , where ψ is the
N
N
X
X
n
n
)
+
log
2
α
+
αi δiβ log(δi ).
ψ(
i
2β
2β 2
i=1
i=1
Proof of Proposition 1
Using that U T U = Id and Ũ T Ũ = In−d , we have that the distribution GG(µ, Σ, β) ⊗ N (0, σ 2 In−d )
has at x = (x1 , x2 ) ∈ Rn the density
1
exp(− 2σ1 2 ∥Ũ T Ũ x2 ∥2 )q(U T U x1 |µ, Σ, β)
(2πσ 2 )(n−d)/2
As U T Ũ = 0 and Ũ T U = 0, we obtain, that this is equal to
1
(2πσ 2 )(n−d)/2
=
exp(− 2σ1 2 ∥Ũ T (U x1 + Ũ x2 )∥2 )q(U T (U x1 + Ũ x2 )|µ, Σ, β)
1
exp(− 2σ1 2 ∥Ũ T V x∥2 )q(U T V x|µ, Σ, β)
(2πσ 2 )(n−d)/2
Using the change of variables formula for push-forward measures, we obtain that Pθ has the density
1
(2πσ 2 )(n−d)/2
exp(− 2σ1 2 ∥Ũ T V r−1 (x)∥2 )q(U T V r−1 (x)|µ, Σ, β)|det(∇r−1 (x))|
As it holds r−1 (x) = V T (x − b) and as ∇r−1 (x) = V T is an orthogonal matrix, we obtain that this
is equal to
1
exp(− 2σ1 2 ∥Ũ T (x − b)∥2 )q(U T (x − b)|µ, Σ, β)
(12)
2
(2πσ )(n−d)/2
where we used V V T = In . Finally, the fact that the mapping x 7→ Ũ x is an isometry and that
I − U U T = Ũ Ũ T imply that
∥Ũ T (x − b)∥2 = ∥Ũ Ũ T (x − b)∥2 = ∥(U U T − I)(x − b)∥2
such that the density (12) coincides with p(x|U, b, µ, Σ, β). Taking the negative logarithm and sub2
tracting n−d
2 log(2πσ ) shows that the negative log-likelihood function of p(·|U, b, µ, Σ, β) coincides
up to a constant with ℓPCA (U, b, µ, Σ, β).
□
C
Proof of Proposition 2
We aim to minimize the Lagrangian function
L(U, Λ) =
N
X
1
αi − 2 ∥U T (xi − b)∥2 + η∥U T (xi − b) − zi ∥2 + tr(Λ(U T U − I)).
2σ
i=1
20
for fixed Λ with respect to U by setting the gradient to zero. Since it holds by [32] that
∂aT U U T b
= U T (abT + baT )
∂U
and
∂tr(ΛU T U )
∂tr(ΛU T U ) ∂U T U )
=
= 2ΛT U T ,
∂U
∂(U T U )
∂U
Thus, it holds that
∇U L(U, Λ) = 2(η −
1
2σ 2 )
N
X
αi (xi − b)(xi − b)T U − 2η
N
X
αi (xi − b)ziT + 2U Λ.
i=1
i=1
Setting the gradient to zero, it shows that U is a solution of
(η −
1
2σ 2 )
N
X
i=1
N
X
αi (xi − b)(xi − b)T U + U Λ = η
αi (xi − b)ziT .
i=1
This finishes the proof.
□
References
[1] D. Zoran and Y. Weiss. From learning models of natural image patches to whole image restoration. In Computer Vision (ICCV), 2011 IEEE International Conference on, pages 479–486, 2011.
[2] S. Parameswaran, C. Deledalle, L. Denis, and T. Q. Nguyen. Accelerating gmm-based patch
priors for image restoration: Three ingredients for a 100x speed-up. IEEE Transactions on Image
Processing, 28(2):687-698, 2019.
[3] C. L. Byrne. The EM Algorithm: Theory, Applications and Related Methods. Lecture Notes,
University of Massachusetts, 2017.
[4] A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood from incomplete data via
the EM algorithm. Journal of the Royal Statistical Society. Series B (Methodological), 39(1):1-38,
1977.
[5] M. Saıd Allili. Wavelet modeling using finite mixtures of generalized Gaussian distributions:
Application to texture discrimination and retrieval. IEEE Transactions on Image Processing,
21(4):1452-1464, 2012.
[6] Z. Boukouvalas, S. Said, L. Bombrun, Y. Berthoumieu, and T. Adali. A new Riemannian averaged fixed-point algorithm for MGGD parameter estimation. IEEE Signal Processing Letters,
22(12):2314-2318, Dec 2015.
[7] F. Pascal, L. Bombrun, J.-Y. Tourneret, and Y. Berthoumieu. Parameter estimation for multivariate generalized Gaussian distributions. IEEE Transactions on Image Processing, 61(23):59605971, 2013.
[8] CA Deledalle, S Parameswaran, and TQ Nguyen. Image denoising with generalized Gaussian
mixture model patch priors. SIAM Journal on Imaging Sciences, 11(4):2568-2609.
[9] B. Wang, H. Zhang, Z. Zhao and Y. Sun, Globally Convergent Algorithms for Learning Multivariate Generalized Gaussian Distributions. IEEE Statistical Signal Processing Workshop (SSP),
2021, pp. 336-340, doi: 10.1109/SSP 49050.2021.9513857.
[10] P. Sandeep and T. Jacob. Single image super-resolution using a joint GMM method. IEEE
Transactions on Image Processing, 25(9): 4233-4244, 2016.
[11] K. Pearson. On lines and planes of closest fit to systems of points in space. Philosophical
Magazine, 2(11):559-572, 1901.
21
[12] M. E. Tipping , C. M. Bishop. Probabilistic Principal Component Analysis. Journal of the royal
statistical society, Series B, Vol. 61, No. 3 (1999), 611-622.
[13] M. E. Tipping and C. M. Bishop. Mixture of probabilistic principal component analyzers. Neural
Computation, 11(2):443-482, 1999.
[14] C. Bouveyron, S. Girard, and C. Schmid. High-dimensional data clustering. Computational
Statistics and Data Analysis, 52(1): 502-519, 2007.
[15] A. Houdard, C. Bouveyron, and J. Delon. High-dimensional mixture model for unsupervised
image denoising (HDMI). SIAM Journal on Imaging Sciences, 11(4):2815-2846, 2018.
[16] J. Hertrich, D-P-L. Nguyen, J-F. Aujol, D. Bernard, Y. Berthoumieu, A. Saadaldin, and G.
Steidl. PCA Reduced Gaussian Mixture Models with Applications in Superresolution, Inverse
Problems in Imaging, 2022, 16 (2) : 341-366. doi: 10.3934/ipi.2021053.
[17] D. Geman and C. Yang. Nonlinear image recovery with half-quadratic regularization. IEEE
Transactions on Image Processing, 5(7):932–946, 1995.
[18] A. P. Dempster , N. M. Laird , D. B. Rubin, Maximum likelihood from incomplete data via the
EM algorithm. Journal of the royal statistical society, SERIES B 1977.
[19] G. McLachlan, D. Peel, Finite Mixture Models, Wiley-Interscience, New York, 2000.
[20] F. Najar, S. Bourouis, N. Bouguila and S. Belghith. A Fixed-point Estimation Algorithm for
Learning The Multivariate GGMM: Application to Human Action Recognition. IEEE Canadian
Conference on Electrical & Computer Engineering (CCECE), 2018.
[21] Tjalling J. Ypma, Historical development of the Newton–Raphson method, SIAM Review 37
(4), 531–551, 1995.
[22] K. J. Arrow, L. Hurwicz, and H. Uzawa. Studies in linear and nonlinear programming. In Stanford Mathematical Studies in the Social Sciences, volume II. Stanford University Press, Stanford,
1958.
[23] E. Lehmann and H. Scheffé, Completeness, similar regions, and unbiased estimation: Part I,
Sankhyā, 10 (1950), 305–340.
[24] E. Gomez, M.A Gomez-Villegas and J.M. Marin. A survey on continuous Elliptical vector
distributions.
[25] H. Hult and L. Filip, Multivariate extremes, aggregation and dependence in elliptical distributions. Advances in Applied probability 34.3 (2002): 587-608.
[26] S. Vaucher, P. Unifantowicz, C. Ricard, L. Dubois, M. Kuball, J.-M. Catala-Civera, D. Bernard,
M. Stampanoni and R. Nicula, On-line tools for microscopic and macroscopic monitoring of microwave processing, Physica B: Condensed Matter, 398 (2007), 191-195.
[27] Lowe, David G. (1999). Object recognition from local scale-invariant features. Proceedings of
the International Conference on Computer Vision. 2. pp. 1150–1157.
[28] F. Altekruger, J. Hertrich. WPPNets: Unsupervised CNN Training with Wasserstein Patch
Priors for Image Superresolution. arXiv preprint arXiv:2201.08157, 2022.
[29] B. Lim, S. Son, H. Kim, S. Nah, and K. Mu Lee, Enhanced deep residual networks for single
image super-resolution. Proceedings of the IEEE Conference on Computer Vision and Pattern
Recognition Workshops, 2017, pp. 136–144.
[30] Y. Zhang, Y. Tian, Y. Kong, B. Zhong, and Y. Fu, Residual dense network for image superresolution. Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition,
2018, pp. 2472–2481.
22
[31] C. Ledig, L. Theis, F. Huszar, J. Caballero, A. Cunningham, A. Acosta, A. Aitken, A. Tejani, J.
Totz, Z. Wang et al., Photo-realistic single image super-resolution using a generative adversarial
network. Proceedings of the IEEE Conference on Computer Vision and pattern Recognition, 2017,
pp. 4681–4690.
[32] Petersen, K.B., Pedersen, M.S.: The Matrix Cookbook. Technical University of Denmark, Lecture Notes (2008).
23