LETTER Communicated by Michael Lewicki
Blind Source Separation by Sparse Decomposition in a Signal
Dictionary
Michael Zibulevsky
Department of Computer Science, University of New Mexico, Albuquerque, NM 87131,
U.S.A.
Barak A. Pearlmutter
Department of Computer Science and Department of Neurosciences, University of New
Mexico, Albuquerque, NM 87131, U.S.A.
The blind source separation problem is to extract the underlying source
signals from a set of linear mixtures, where the mixing matrix is unknown.
This situation is common in acoustics, radio, medical signal and image
processing, hyperspectral imaging, and other areas. We suggest a two-
stage separation process: a priori selection of a possibly overcomplete
signal dictionary (for instance, a wavelet frame or a learned dictionary) in
which the sources are assumed to be sparsely representable, followed by
unmixing the sources by exploiting the their sparse representability. We
consider the general case of more sources than mixtures, but also derive a
more efficient algorithm in the case of a nonovercomplete dictionary and
an equal numbers of sources and mixtures. Experiments with artificial
signals and musical sounds demonstrate significantly better separation
than other known techniques.
1 Introduction
In blind source separation an N-channel sensor signal x(t) arises from M un-
known scalar source signals s
i
(t), linearly mixed together by an unknown
N × M matrix A, and possibly corrupted by additive noise ξ(t),
x(t) = As(t) + ξ(t). (1.1)
We wish to estimate the mixing matrix A and the M-dimensional source
signal s(t). Many natural signals can be sparsely represented in a proper
signal dictionary:
s
i
(t) =
K
X
k=1
C
ik
ϕ
k
(t). (1.2)
The scalar functions ϕ
k
(t) are called atoms or elements of the dictionary.
These elements do not have to be linearly independent and instead may
Neural Computation 13, 863–882 (2001)
c
° 2001 Massachusetts Institute of Technology
864 Michael Zibulevsky and Barak A. Pearlmutter
form an overcomplete dictionary. Important examples are wavelet-related
dictionaries (e.g., wavelet packets, stationary wavelets; see Chen, Donoho,
& Saunders, 1996; Mallat, 1998) and learned dictionaries (Lewicki & Se-
jnowski, in press; Lewicki & Olshausen, 1998; Olshausen & Field, 1996,
1997). Sparsity means that only a small number of the coefficients C
ik
differ
significantly from zero.
We suggest a two-stage separation process: a priori selection of a possibly
overcomplete signal dictionary in which the sources are assumed to be
sparsely representable and then unmixing the sources by exploiting their
sparse representability.
In the discrete-time case t = 1, 2,...,T we use matrix notation. X is an
N × T matrix, with the ith component x
i
(t) of the sensor signal in row i, S is
an M × T matrix with the signal s
j
(t) in row j, and 8 is a K × T matrix with
basis function ϕ
k
(t) in row k. Equations 1.1 and 1.2 then take the following
simple form:
X = AS + ξ (1.3)
S = C8. (1.4)
Combining them, we get the following when the noise is small:
X ≈ AC8.
Our goal therefore can be formulated as follows: Given the sensor signal ma-
trix X and the dictionary 8, find a mixing matrix A and matrix of coefficients
C such that X ≈ AC8 and C is as sparse as possible.
We should mention other problems of sparse representation studied in
the literature. The basic problem is to represent sparsely scalar signal in
given dictionary (see Chen et al., 1996). Another problem is to adapt the
dictionary to the given class of signals
1
(Lewicki & Sejnowski, 1998; Lewicki
& Olshausen, 1998; Olshausen & Field, 1997). This problem is shown to
be equivalent to the problem of blind source separation when the sources
are sparse in time (Lee, Lewicki, Girolami, & Sejnowski, 1999; Lewicki &
Sejnowski, in press). Our problem is different, but we will use and generalize
some techniques presented in these works.
Independent factor analysis (Attias, 1999) and Bayesian blind source sep-
aration (Rowe, 1999) also consider the case of more sources than mixtures.
In our approach, we take an advantage when the sources are sparsely repre-
sentable. In the extreme case, when the decomposition coefficients are very
sparse, the separation becomes practically ideal (see section 3.2 and the six
flutes example in Zibulevsky, Pearlmutter, Bofill, & Kisilev, in press). Never-
theless, detailed comparison of the methods on real-world signals remains
open for future research.
1
Our dictionary 8 may be obtained in this way.
Blind Source Separation by Sparse Decomposition 865
In section 2 we give some motivating examples, which demonstrate how
sparsity helps to separate sources. Section 3 gives the problem formula-
tion in probabilistic framework and presents the maximum a posteriori
approach, which is applicable to the case of more sources than mixtures.
In section 4 we derive another objective function, which provides more ro-
bust computations when there is an equal number of sources and mixtures.
Section 5 presents sequential source extraction using quadratic program-
ming with nonconvex quadratic constraints. Finally, in section 6 we derive
a faster method for nonovercomplete dictionaries and demonstrate high-
quality separation of synthetically mixed musical sounds.
2 Separation of Sparse Signals
In this section we present two examples that demonstrate how sparsity of
source signals in the time domain helps to separate them. Many real-world
signals have sparse representations in a proper signal dictionary but not in
the time domain. The intuition here carries over to that situation, as shown
in section 3.1.
2.1 Example: Two Sources and Two Mixtures. Two synthetic sources
are shown in Figures 1a and 1b. The first source has two nonzero samples,
and the second has three. The mixtures, shown in Figures 1c and 1d, are less
sparse: they have five nonzero samples each. One can use this observation
to recover the sources. For example, we can express one of the sources as
e
s
i
(t) = x
1
(t) + µx
2
(t)
and choose µ so as to minimize the number of nonzero samples k
e
s
i
k
0
, that
is, the l
0
norm of s
i
.
This objective function yields perfect separation. As shown in Figure 2a,
when µ is not optimal, the second source interferes, and the total number
of nonzero samples remains five. Only when the first source is recovered
perfectly, as in Figure 2b, does the number of nonzero samples drop to two
and the objective function achieve its minimum.
Note that the function k
e
s
i
k
0
is discontinuous and may be difficult to
optimize. It is also very sensitive to noise: even a tiny bit of noise would
make all the samples nonzero. Fortunately in many cases, the l
1
norm k
e
s
i
k
1
is a good substitute for this objective function. In this example, it too yields
perfect separation.
2.2 Example: Three Sources and Two Mixtures. The signals are pre-
sented in Figure 3. These sources have about 10% nonzero samples. The
nonzero samples have random positions and are zero-mean unit-variance
gaussian distributed in amplitude. Figure 4 shows a scatter plot of the mix-
tures. The directions of the columns of mixing matrix are clearly visible.
866 Michael Zibulevsky and Barak A. Pearlmutter
(a)
(b)
(c)
(d)
Figure 1: Coefficients of signals, with coefficient identity on the x-axis (10 co-
efficients, arbitrarily ordered) and magnitude on the y-axis (arbitrarily scaled).
Sources (a and b) are sparse. Mixtures (c and d) are less sparse.
Figure 2: Coefficients of signals, with coefficient identity on the x-axis (10 co-
efficients, arbitrarily ordered) and magnitude on the y-axis (arbitrarily scaled).
(a) Imperfect separation. Since the second source is not completely removed,
the total number of nonzero samples remains five. (b) Perfect separation. When
the source is recovered perfectly, the number of nonzero samples drops to two,
and the objective function achieves its minimum.
This phenomenon can be used in clustering approaches to source separa-
tion (Pajunen, Hyvrinen, & Karhunen, 1996; Zibulevsky et al., in press). In
this work we will explore a maximum a posteriori approach.
Blind Source Separation by Sparse Decomposition 867
Figure 3: Coefficients of signals, with coefficient identity on the x-axis (300 co-
efficients, arbitrarily ordered) and magnitude on the y-axis (arbitrarily scaled).
(Top three panels) Sparse sources (sparsity is 10%). (Bottom two panels) Mix-
tures.
3 Probabilistic Framework
In order to derive a maximum a posteriori solution, we consider the blind
source separation problem in a probabilistic framework (Belouchrani & Car-
doso, 1995; Perlmutter & Parra, 1996). Suppose that the coefficients C
ik
in a
source decomposition (see equation 1.4) are independent random variables
with a probability density function (pdf) of an exponential type,
p
i
(C
ik
) ∝ exp −β
i
h(C
ik
). (3.1)
This kind of distribution is widely used for modeling sparsity (Lewicki &
Sejnowski, in press; Olshausen & Field, 1997). A reasonable choice of h(c)
may be
h(c) =|c|
1/γ
γ ≥ 1 (3.2)
