L1-NORM PENALIZED LEAST SQUARES WITH SALSA

IVAN SELESNICK

Abstract. This lecture note describes an iterative optimization algorithm, ‘SALSA’, for solving L1-norm
penalized least squares problems. We describe the use of SALSA for sparse signal representation and
approximation, especially with overcomplete Parseval transforms. We also illustrate the use of SALSA to
perform basis pursuit (BP), basis pursuit denoising (BPD), and morphological component analysis (MCA).
The algorithm, ‘SALSA’, was developed by Afonso, Bioucas-Dias, and Figueiredo.

1. Introduction

Numerous sparsity-based signal processing methods are based on `1 -norm penalized least squares. This
approach has been used for denoising, deconvolution, missing data estimation, signal separation, and other
problems. It has been demonstrated that combining the augmented Lagrangian approach and the variable
splitting technique is an effective algorithmic approach for solving linear inverse problems with sparse regu-
larization [1]. An algorithm, called SALSA, developed in Ref. [1], is based on this approach. This algorithm
is notable due to (i ) its flexibility in handling various problems, and (ii ) its fast convergence in practice.
More generally, the alternating direction method of multipliers (ADMM) has been shown lately to be highly
effective for large scale non-smooth optimization [3].
This note is intended to complement the tutorial [7] which intentionally omitted detailed descriptions of
algorithms for solving the `1 -norm optimization problems described therein. In particular, this note describes
the derivation of SALSA to solve two problems. The first problem is `1 -norm penalized least squares; i.e.,
1
xopt = arg min ky − Axk22 + λkxk1 . (BPD)
x 2
The second problem is that of finding the solution to a system of linear equations with minimal `1 -norm;
i.e.,
xopt = arg min kxk1 (BP)
x

such that Ax = y.
These problems are sometimes referred to as basis pursuit denoising (BPD) and basis pursuit (BP), respec-
tively [4].
For a vector x ∈ CN , the `1 and `2 norms are defined by
N
X −1 N
X −1
kxk1 := |xn |, kxk22 := |xn |2 . (1)
n=0 n=0
Sections 2 and 3 derive iterative algorithms to solve BPD and BP, respectively. Based on these algorithms,
Section 4 derives iterative algorithms for ‘dual BPD’ and ‘dual BP’. These algorithms can be used to
implement morphological component analysis (MCA) for nonlinear signal decomposition. Section 5 describes

Date: November, 2011. Last edit: January 20, 2014.

Matlab software available at [Link]
This tutorial is a Connexions module (m48933 at [Link]
Cite as: I. Selesnick. L1-Norm Penalized Least Squares with SALSA. Connexions, 2014. [Link]
Support from NSF under Grant CCF-1018020 is gratefully acknowledged.
1
2 IVAN SELESNICK

transforms, A, useful for sparse signal representation and approximation. Section 6 presents examples of
BP, BPD, dual BP, and dual BPD applied to simple signals.

1.1. The Augmented Lagrangian. For the constrained optimization problem,

arg min E(z) (2)

z

such that Cz − b = 0,

the augmented Lagrangian is defined as

LA (z, α, µ) = E(z) + αT (Cz − b) + µkCz − bk22 . (3)

The vector, α, are Lagrange multipliers. Version-2 of the augmented Lagrangian method (ALM) [1], to solve
the constrained problem is given by

initialize: µ > 0, d
repeat
z ← arg min E(z) + µkCz − dk22 (4a)
z

d ← d − (Cz − b) (4b)
end

The indented assignment operations are iterated until convergence. This method is also known as the method
of multipliers (MM); so this iterative algorithm is referred to as ALM/MM in [1].
The ALM/MM algorithm calls for a positive scalar, µ, which is like a step-size parameter. Its value can
affect the convergence speed of the algorithm. But it does not affect the solution to which it converges.

2. L1 norm regularized least squares (BPD)

Given an observed vector y and matrix A, consider the problem of finding a sparse vector x such that
y ≈ Ax. Using the `1 norm as a measure of sparsity, the problem can be formulated as
1
xopt = arg min ky − Axk22 + kλ xk1 (5)
x 2
The notation λ x denotes element-wise multiplication of the equal-size vectors λ and x; i.e., [λ x]i = λi xi .
When all elements of vector λ are the same value (i.e., λi = λ ∈ R+ ), then (5) can be written as
1
arg min ky − Axk22 + λkxk1 (6)
x 2
which is the more common form. However, it will sometimes be useful to allow non-uniform regularization
of x, so we will use the form (5).
Applying variable splitting to (5) yields
1
arg min ky − Axk22 + kλ uk1 (7)
2
x,u

such that u − x = 0

Variable splitting introduces an auxiliary variable, but it also decouples the terms of the objective function.
(Actually, it moves the coupling into the constraint, which is handled subsequently through alternating
minimization.)
Problem (7) can be put written in the form of (2) by setting
" #
z1 h i
z1 = x, z2 = u, z = , C = I, −I , b = 0, (8)
z2
 L1-NORM PENALIZED LEAST SQUARES WITH SALSA 3

and
1
E(z) = ky − Az1 k22 + kλ z2 k1 . (9)
2
Now that the problem is expressed in the form of (2), the ALM/MM algorithm (4) can be applied. Using
ALM/MM to solve problem (7), we obtain the iterative algorithm:

initialize: µ > 0, d
repeat
1 µ
x, u ← arg min ky − Axk22 + kλ uk1 + ku − x − dk22 (10a)
2
x,u 2
d ← d − (u − x) (10b)
end

The vector d must be initialized prior to the iteration. We usually initialize d to the zero vector.
As proven by Eckstein and Bertsekas, in a more general setting, if the minimization in (10a) is performed
alternately between x and u, the algorithm will still converge to the global minimum [5]. This technique
is known as alternating direction method of multipliers (ADMM). Alternating between minimization with
respect to each of x and u, we obtain the algorithm:

initialize: µ > 0, d
repeat
µ
u ← arg min kλ ku − x − dk22
uk1 + (11a)
u 2
1 µ
x ← arg min ky − Axk22 + ku − x − dk22 (11b)
x 2 2
d ← d − (u − x) (11c)
end

This algorithm is called SALSA (split augmented Lagrangian shrinkage algorithm) in Ref. [1]. In fact,
SALSA is more general, as it allows a general regularizer, φ(x), not just the `1 norm.
The minimizations (11a) and (11b) can be performed in explicit form. The minimization problem in (11a)
is separable in ui . Its solution is expressed explicitly in terms of the soft-thresholding rule (see Appendix A).
The minimization problem in (11b) is a constrained least squares problem; hence, its solution is available in
explicit form (in terms of a matrix inverse). Utilizing the explicit forms for the two minimization problems,
we obtain the following algorithm.

initialize: µ > 0, d
repeat
u ← soft(x + d, λ/µ) (12a)

−1
x ← AH A + µ I AH y + µ (u − d)


(12b)
d←d−u+x (12c)
end

The operator AH is the complex conjugate (Hermitian) transpose of A.

With a change of variables, v = u − d, the arithmetic operations can be slightly reduced.
4 IVAN SELESNICK

Algorithm 1: Algorithm for basis pursuit denoising (5).

initialize: µ > 0, d
repeat
v ← soft(x + d, λ/µ) − d (13a)

−1
x ← AH A + µ I AH y + µ v


(13b)
d←x−v (13c)
end

Sometimes, in (13b) it can be useful to use the matrix inverse lemma (see Appendix B) to write

−1 1 1
−1
µ I + AH A = I − AH µ I + AAH A, (14)
µ µ
because in certain cases, µ I + AAH is easier to invert than µ I + AH A.

2.1. When A is a Tight Frame. In some signal processing applications, A is a ‘wide’ matrix satisfying

A AH = p I, p > 0. (15)

In this case, it is sometimes said that the columns of A form a tight frame. The matrix AH can also be
considered an overcomplete Parseval transform. For example, the columns of A may be an overcomplete set
of a complex sinusoids with closely spaced frequencies.
In many cases, we will have p = 1 in (15). However, for some problems (e.g., dual BPD in Sect. 4) we
will have p 6= 1.
Using (15) in (14), we obtain:

−1 1 1
µ I + AH A = I− AH A. (16)
µ µ(µ + p)
Then the update equation for x in (13b) becomes:
1 1
AH y + µ v − AH A AH y + µ v



x← (17)
µ µ(µ + p)
which simplifies to
1 H p 1
x← A y+v− AH y − AH A v (18a)
µ µ(µ + p) µ+p
1 1
= AH y + v − AH A v (18b)
µ+p µ+p
1
=v+ AH (y − A v) (18c)
µ+p
Therefore, Algorithm 1 can be written as follows.

initialize: µ > 0, d
repeat
v ← soft(x + d, λ/µ) − d (19a)
1
x←v+ AH (y − A v) (19b)
µ+p
d←x−v (19c)
end

The algorithm can be simplified by a slight rearrangement of operations, as follows.

 L1-NORM PENALIZED LEAST SQUARES WITH SALSA 5

Algorithm 2: Algorithm for basis pursuit denoising (5) with AAH = p I.

initialize: µ > 0, d
repeat
v ← soft(x + d, λ/µ) − d (20a)
1
d← AH (y − A v) (20b)
µ+p
x←d+v (20c)
end
Note that this algorithm does not involve any matrix inverse. If fast implementations are available for A
and AH , then each iteration of the algorithm is fast.

3. L1 Norm Regularized Solutions to Linear Systems (BP)

Given an observed signal y, consider the problem of finding a sparse vector x that solves Ax = y. Using
the `1 norm as a measure of sparsity, the problem can be formulated as:
arg min kλ xk1 (21)
x

such that A x = y
This problem is known as basis pursuit [4]. By applying the variable splitting technique, we obtain an
equivalent optimization problem:
arg min kλ uk1 (22)
x,u

such that A x = y
u−x=0
We will use the ‘partly’ augmented Lagrangian:
LA (x, u, λ, µ) = kλ uk1 + λT (u − x) + 0.5µku − xk22 + λ2 (A x − y) (23)
Using ALM/MM so solve the problem, we obtain the algorithm:
initialize: µ > 0, d
repeat

arg min kλ uk1 + 0.5µku − x − dk22
x, u ← x,u (24a)
such that A x = y

d ← d − (u − x) (24b)
end
By alternately minimizing with respect to x and u (as in Sec. 2), we obtain the algorithm:
initialize: µ > 0, d
repeat
u ← arg min kλ uk1 + 0.5µku − x − dk22 (25a)
u

arg min ku − x − dk2
2
x← x (25b)
such that A x = y
6 IVAN SELESNICK

d ← d − (u − x) (25c)
end
The minimization with respect to u in (25a) can be expressed explicitly in terms of soft-thresholding. The
minimization with respect to x in (25b) is a constrained least squares problem which admits an explicit
solution in terms of matrix inverses. Using the explicit solution to each of the two minimization problems,
we obtain the algorithm:
initialize: µ > 0, d
repeat
u ← soft(x + d, λ/µ) (26a)
x ← (u − d) + AH (AAH )−1 (y − A(u − d)) (26b)
d ← d − (u − x) (26c)
end
With a change of variables, v = u − d, the arithmetic operations can be slightly reduced, as follows.
initialize: µ > 0, d
repeat
v ← soft(x + d, λ/µ) − d (27a)
x ← v + AH (AAH )−1 (y − Av) (27b)
d←x−v (27c)
end
The algorithm can be further simplified by a slight rearrangement of operations, as follows.
Algorithm 3: Algorithm for basis pursuit (21).
initialize: µ > 0, d
repeat
v ← soft(x + d, λ/µ) − d (28a)
H H −1
d ← A (AA ) (y − Av) (28b)
x←d+v (28c)
end
Note that at every iteration, x satisfies Ax = y. This is because
A(d + v) = A AH (AAH )−1 (y − Av) + v
 
(29a)
= AAH (AAH )−1 (y − Av) + Av (29b)
= (y − Av) + Av (29c)
=y (29d)

3.1. When A is a Tight Frame. Consider the BP problem (21) when the columns of A form a tight
frame; i.e., when A satisfies (15). Then Algorithm 3 can be written as follows.
Algorithm 4: Algorithm for basis pursuit (21) with AAH = p I.
initialize: µ > 0, d
 L1-NORM PENALIZED LEAST SQUARES WITH SALSA 7

repeat
v ← soft(x + d, λ/µ) − d (30a)
1 H
d← A (y − Av) (30b)
p
x←d+v (30c)
end
Note that this is very similar to Algorithm 2; only the constant in (20b) is different. Likewise, if A and AH
are fast, then then the algorithm as a whole is fast.

4. Dual BP and dual BPD

In several signal processing applications, it is useful to to model a signal y as

y ≈ A1 x 1 + A2 x 2 . (31)
In particular, this model is used in morphological component analysis (MCA) for the nonlinear separation
of signal components [11, 10]. There are several ways to formulate the MCA problem. Two approaches for
MCA are based on forms of BP and BPD.
In the following, we assume that Ai are tight frames with frame constant p = 1; i.e.,
A1 AH
1 = I, A2 AH
2 = I. (32)

4.1. Dual BPD. If the signal, y, is noisy, then it is appropriate to allow a residual. In this case, MCA may
be formulated as,
1
arg min ky − A1 x1 − A2 x2 k22 + kλ1 x1 k1 + kλ2 x2 k1 . (33)
x1 ,x2 2
This is a special case of BPD (5) with
" # " #
h i x1 λ1
A = A1 A2 , x = , λ= . (34)
x2 λ2
Since Ai are tight frames (32), we have:
" #
h i AH
1
AAH = A1 A2 = A1 AH H
1 + A2 A2 = 2 I. (35)
AH2

Therefore, we can use Algorithm 2 in Sect. 2.1 with p = 2. Hence, we obtain the following algorithm for
dual BPD.
initialize: µ > 0, d
repeat
" # " # " # " #! " #
v1 x1 d1 λ1 /µ d1
← soft + , − (36a)
v2 x2 d2 λ2 /µ d2
" # " # " #!
d1 1 AH1
h i v
1
← y − A 1 A 2 (36b)
d2 µ + 2 AH 2 v 2
" # " # " #
x1 d1 v1
← + (36c)
x2 d2 v2
end
This algorithm can be expressed as follows.
8 IVAN SELESNICK

Algorithm 5: Algorithm for dual BPD (33) with Ai AH

i = I.

initialize: µ > 0, di
repeat
vi ← soft(xi + di , λi /µ) − di , i = 1, 2 (37a)
c ← y − A1 v 1 − A2 v 2 (37b)
1
di ← AH c, i = 1, 2 (37c)
µ+2 i
xi ← di + vi , i = 1, 2 (37d)
end

4.2. Dual Basis Pursuit. If the signal, y, is noise-free, then it is appropriate to use an equality constraint.
In this case MCA may be formulated as,
arg min kλ1 x1 k1 + kλ2 x2 k1 (38)
x1 ,x2

such that y = A1 x1 + A2 x2 .
This is a special case of BP (21) with
" # # "
h i x1 λ1
A = A1 A2 , x = , λ= . (39)
x2 λ2
Since we assume that Ai are tight frames (32), we have (35). Therefore, we can use Algorithm 4 with p = 2.
Hence, we obtain the following algorithm for dual BP.
initialize: µ > 0, d
repeat
" # " # " # " #! " #
v1 x1 d1 λ1 /µ d1
← soft + , − (40a)
v2 x2 d2 λ2 /µ d2
" # " # " #!
d1 1 AH 1
h i v
1
← H
y − A 1 A 2 (40b)
d2 2 A2 v2
" # " # " #
x1 d1 v1
← + (40c)
x2 d2 v2
end
This algorithm can be expressed as follows.
Algorithm 6: Algorithm for dual BP (38) with Ai AH
i = I.

initialize: µ > 0, di
repeat
vi ← soft(xi + di , λi /µ) − di , i = 1, 2 (41a)
c ← y − A1 v1 − A2 v2 (41b)
1 H
di ← A c, i = 1, 2 (41c)
2 i
xi ← di + vi , i = 1, 2 (41d)
end
 L1-NORM PENALIZED LEAST SQUARES WITH SALSA 9

Note that this is the same as Algorithm 5 except for a constant in (37c).

Note that Algorithms 5 and 6 involve no matrix inverses. If Ai and AH

i are fast, then these algorithms as
a whole are fast. For example, if the Ai are FFTs and/or short-time Fourier transforms, then they have low
implementation complexity and can be admit high parallelism. Such a combination of transforms is useful
for decomposing a signal into narrow-band and wide-band signal components, even when the components
overlap in both time and frequency [9]. Alternately, by taking the Ai as wavelet transforms with different
Q-factors, a signal can be decomposed into low and high resonance components [8].

5. Transforms for Sparse Signal Representation

In order to apply BP and BPD, a transform, A, is need for the sparse representation of the signal of
interest. In dual BP and dual BPD, two transforms, A1 and A2 , are needed. The transforms should be
chosen such that they enables a sparse representation (or approximation) of the signals of interest.
To emphasize that the representation of a signal y is in terms of transform coefficients, we use the letter
‘c’ for coefficients. For example, we write y = Ac as a representation of signal y with respect to transform
A where c is the vector of transform coefficients.

5.1. Zero-padded DFT. To sparsely represent a real or complex set of sinusoids, we take A to be the
normalized inverse of an K-point DFT with K ≥ N . Specifically, A : CK → CN is defined by
√
[Ac]n = K [DFT−1
K {c}]n , c ∈ CK , n ∈ ZN . (42)

where DFTK is the discrete Fourier transform (DFT) with zero-padding up to a total length of K. The
√
multiplication by K normalizes A so that AAH = IN . Note that, when K > N , the matrix A is ‘wide’
rather than square. Accordingly, in the definition of A, the inverse K-point DFT is truncated down to N
samples. The matrix A is a sub-matrix of the the inverse K-point DFT matrix (the first N rows of the
K × K inverse DFT matrix). Consequently, AH : CN → CK is defined by
1
[AH x]k = √ [DFTK {x}]k , x ∈ CN , k ∈ ZK (43)
K
where the N -point vector x is zero-padded to length K prior to the DFT computation.
The `2 norm of all the columns of A are equal, specifically,
r
N
kak2 = . (44)
K
This value can be used for setting regularization parameters, λ, in BPD and dual BPD. When the DFT is
critically sampled (i.e., K = N ), then A is simply the conventional DFT, normalized so as to be unitary, in
which case (44) gives unity, as expected.
The DFT operator A can be implemented in MATLAB as
truncate = @(c, N) c(1:N);
A = @(c) sqrt(K) * truncate(ifft(c), N);
and AH as
AH = @(x) fft(x, K)/sqrt(K);
These fast matrix-free implementations of A and AH can be used for dual BP and dual BPD only if the
utilized optimization algorithms are also ‘matrix-free’, as are the SALSA algorithms.
10 IVAN SELESNICK

5.2. STFT. To sparsely represent a signal composed of oscillatory pulses, we take A to be the normalized
inverse of a short-time Fourier transform (STFT). The STFT has several parameters: the frame length,
overlapping factor, and DFT length. We typically use 50% overlapping and a DFT length equal to at least
the frame length. Consequently, the STFT is at least two-times over-sampled. If the time-frequency array
of STFT coefficients is of size M × K, for a signal of length N , then A : CM ×K → CN is defined as

[Ac]n = [STFT−1 {c}]n , n ∈ ZN (45)

and AH : CN → CM ×K is defined as

[AH x](m,k) = [STFT{x}](m,k) , m ∈ ZM , k ∈ ZK . (46)

With a suitably implemented STFT, we have AAH = IN .

The `2 norm of all the columns of A are equal. Specifically, when 50% overlapping is used, the norm is
given by
r
R
kak2 = . (47)
2K
This value can be used for setting regularization parameters, λ, in BPD and dual BPD.
We implement the STFT operators A and AH in MATLAB as
A = @(c) ipSTFT(c, R, N);
AH = @(x) pSTFT(x, R, K);
where pSTFT and ipSTFT are our implementations of the STFT, and its inverse, designed to satisfy AAH =
IN . The ‘p’ stands for ‘Parseval’. The parameter R is the frame length. The parameter K is the DFT
length, with K ≥ R. (The function pSTFT is not the built-in MATLAB spectrogram function, which is not
designed to ensure invertibility.) The implementation of the STFT, so as to satisfy AAH = I, is described
in [6].

6. Examples

Example 1 (BP). This example illustrates the sparse representation of complex sinusoids in white complex
noise using BP with a zero-padded DFT. We assume that the signal y admits a sparse representation of the
form

y = Ac, y ∈ CN , c ∈ CK , A ∈ CN ×K (48)

where A is a zero-padded DFT and c is a sparse set of DFT coefficients. Given y, we find a sparse c by
solving the basis pursuit problem,

copt = arg min kck1 (49)

c

such that Ac = y.

Figure 1 shows the (real part of the) complex signal, y, in (a); the DFT coefficients, AH y in (b); and the
BPD-optimized coefficients, copt , in (c). In this example, the signal y is of length N = 100 samples. In Fig. 1
(and subsequent figures), the sampling rate is taken to be one sample/sec. in the axis labeling.

Example 2 (BPD). This example illustrates the estimation of real-valued sinusoids in white noise using
basis pursuit denoising (BPD). We assume the noisy data, y, is given by

y = Ac + w, y, w ∈ CN , c ∈ CK , A ∈ CN ×K (50)
 L1-NORM PENALIZED LEAST SQUARES WITH SALSA 11

2
Signal
1

real( y )
0

−1

−2
0 20 40 60 80 100
Time (sec)

8
DFT of signal
6
abs( AH y )

0
0 0.2 0.4 0.6 0.8 1
Frequency (Hz)

15
BP−optimized DFT coefficients

10
abs( c )

0
0 0.2 0.4 0.6 0.8 1
Frequency (Hz)

Figure 1. Example 1. Sparse representation of complex sinusoids in white complex noise

using BP with a zero-padded DFT. (a) Complex sinusoids, y. (b) DFT of data, AH y. (c)
Sparse DFT coefficients, copt , obtained by solving the BP problem.

where c is a sparse set of DFT coefficients and w is a white Gaussian vector. This example uses real-valued
sinusoids and noise, but we may still assume complex-valued signals to solve the problem. Due to conjugate
symmetry properties of the DFT, the denoised signal will be real-valued.
Given y, we estimate the sinusoids by solving the BPD problem,
1
copt = arg min ky − Ack22 + λkck1 . (51)
c 2
For the transform, A, we use the zero-padded DFT (100 signal samples in the time domain, 256 DFT
coefficients in the frequency domain). The optimal coefficients, copt , are obtained using Algorithm 2. Figure 2
shows the noisy data, y, in (a); the noisy DFT coefficients, AH y in (b); the BPD-optimized coefficients, copt ,
in (c); and the denoised signal, Acopt , in (d). In the figure, only the positive frequency axis (0 ≤ f ≤ 0.5) is
shown because the DFT coefficients are conjugate symmetric.

Example 3 (BPD). This example illustrates the estimation of a real-valued pulse in white noise using basis
pursuit denoising (BPD). We assume the noisy data, y, is given by

y = Ac + w, y, w ∈ CN , (52)

where c is a sparse set of STFT coefficients and w is a white Gaussian vector.

12 IVAN SELESNICK

3
Noisy data
2
1
0

y
−1
−2
−3
0 20 40 60 80 100
Time (sec)

3
DFT of noisy data
abs( AH y )

0
0 0.1 0.2 0.3 0.4 0.5
Frequency (Hz)

5
BPD DFT (beta = 2.50)
4
abs( c )

0
0 0.1 0.2 0.3 0.4 0.5
Frequency (Hz)

3
Denoised data
2
1
real( Ac )

0
−1
−2
−3
0 20 40 60 80 100
Time (sec)

Figure 2. Example 2. Estimation of sinusoids in white noise using BPD with a zero-
padded DFT. (a) Sinusoids in white noise, y. (b) DFT of noisy data, AH y. (c) Sparse DFT
coefficients, copt , obtained by solving the BPD problem. (d) Denoised signal, Acopt .

Given y, we estimate the pulse by solving the dual BP problem,

1
copt = arg min ky − Ack22 + λkck1 . (53)
c 2
For the transform, A, we use the STFT. The optimal coefficients, copt , are obtained using Algorithm 2.
Figure 3 shows the noisy data, y, in (a); the noisy STFT coefficients, AH y in (b); the BPD-optimized
coefficients, copt , in (c); and the denoised signal, ŝ = Acopt , in (d).

Example 4 (dual-BP). We illustrate the separation of a sinusoid and a pulse in the noise-free case. The
data is given by

y = s1 + s2 , y, s1 , s2 ∈ RN . (54)
 L1-NORM PENALIZED LEAST SQUARES WITH SALSA 13

3
Noisy data
2

y
−1

−2

−3
0 20 40 60 80 100
Time (sec)

0.5
STFT of noisy data
0.4
Frequency (Hz)

0.3

0.2

0.1

0
0 20 40 60 80 100
Time (sec)

0.5
BPD STFT coefficients (beta = 2.50)
0.4
Frequency (Hz)

0.3

0.2

0.1

0
0 20 40 60 80 100
Time (sec)

3
Denoised signal
2

1
real( Ac )

−1

−2

−3
0 20 40 60 80 100
Time (sec)

Figure 3. Example 3. Estimation of a pulse in white noise using BPD with a STFT. (a)
Pulse in white noise, y. (b) STFT of noisy data, AH y. (c) Sparse STFT coefficients, copt ,
obtained by solving the BPD problem. (d) Denoised signal, Acopt .

Given y, we estimate s1 and s2 by solving the dual BP problem,

{copt opt
1 , c2 } = arg min λ1 kc1 k1 + λ2 kc2 k1 (55)
c1 ,c2

such that y = A 1 c1 + A 2 c2 .

The optimal coefficients, ci , are found using Algorithm 6. We then set ŝi = Ai copt
i for i = 1, 2.
The example is illustrated in Fig. 4. For transform A1 , we use a zero-padded DFT with two-times
oversampling; i.e., N = 100 and K = 200 in (42) and (43). For transform A2 , we use the STFT with frame
14 IVAN SELESNICK

length R = 16, DFT length K = 16, and 50% overlapping. Hence, both A1 and A2 are oversampled by two,
√
and ka1 k2 = ka1 k2 = 1/ 2.
In the dual BP problem, we set λ1 = λ2 = 0.5. To solve the dual BP problem, we run Algorithm 6. As a
result, we obtain sparse coefficients, copt opt
i , as illustrated in Fig. 4. The coefficient vector, c1 , is a sparse set
opt
of DFT coefficients. The coefficient array, c2 , is a sparse set of STFT coefficients, displayed as an image
in a time-frequency plane. From the coefficients, we construct the two signals, x1 and x2 as xi = Ai ci , for
i = 1, 2.

Example 5 (dual-BPD). We illustrate the separation and estimation of a sinusoid and a pulse in the case
of additive white Gaussian noise. The data is given by
y = s1 + s2 + w, y, s1 , s2 , w ∈ RN (56)
where w is a zero-mean white Gaussian vector with variance σ 2 . Figure 5 shows the noisy signal, y.
Given y, we estimate s1 and s2 by solving the dual BPD problem,
1
{copt opt
1 , c2 } = arg cmin ky − A1 c1 − A2 c2 k22 + λ1 kc1 k1 + λ2 kc2 k1 . (57)
1 2 2
,c

The components si are then estimated as,

ŝi = Ai copt
i , i = 1, 2. (58)
For transform A1 , we use a zero-padded DFT; i.e., N = 100 and K = 256 in (42) and (43) (100 signal
samples in the time domain, 256 DFT coefficients in the frequency domain). Hence, the dual BPD algorithm
p
can be run using radix-2 FFTs exclusively. For this A1 , we have ka1 k2 = 100/256 = 5/8. For transform
A2 , we use the STFT with frame length R = 16, DFT length K = 16, and 50% overlapping. For this A2 ,
√
we have ka2 k2 = 1/ 2.
In the dual BPD problem, we set λ1 = βka1 k2 σ, and λ2 = βka2 k2 σ where β = 2.5. This choice of λi is
discussed in [other notes]. Generically, one may set β ∈ [2.5, 3]. To solve the dual BPD problem, we run
Algorithm 5. As a result, we obtain sparse coefficients, copt opt
1 and c2 , as illustrated in Fig. 5. (c1 is a sparse
vector of DFT coefficients, c2 is a sparse two-dimensional array of STFT coefficients). From the coefficients,
we obtain the two signals, x1 and x2 as xi = Ai ci , for i = 1, 2.

Example 6 (dual-BP). This example illustrates dual BPD with a speech waveform.
To be completed . . . .

7. Conclusion

This note has described the algorithm, SALSA, for standard `1 norm minimization problems arising in
sparse signal processing. SALSA can also be used for more general problems (not only the quadratic data
fidelity) and more general regularization terms (not only the `1 norm penalty); see Ref. [1] for more details.
An extension of SALSA to the constrained formulation of the sparsity-penalized least squares problem, called
CSALSA, is developed in Ref. [2].
In this note, we emphasize that when Ai are tight frames, then the presented algorithms for BP, BPD,
dual BP, and dual BPD are:
(1) Matrix-free: The Ai and AH i appear only at operators. No elements of Ai need to be individually
accessed. Hence, Ai do not need to be stored as matrices. It is sufficient to implement the operators
as algorithms. Fast algorithm for Ai and Ai can be exploited.
(2) Low complexity: The main computation is Ai and AH i .
(3) Globally convergent: Any initialization leads to an optimal solution (the objective functions are
convex).
 L1-NORM PENALIZED LEAST SQUARES WITH SALSA 15

4
y = s1 + s2

−2

−4
0 10 20 30 40 50 60 70 80 90 100
Time (sec)

5
DFT coefficients (c1) [Output of dual−BP]
4
abs( c1 )

0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
Frequency (Hz)

0.5
STFT coefficients (c2) [Output of dual−BP]
0.4
Frequency (Hz)

0.3

0.2

0.1

0
0 10 20 30 40 50 60 70 80 90 100
Time (sec)

4
Sinusoid component (reconstructed from c1)

2
real(A1 c)1)

−2

−4
0 10 20 30 40 50 60 70 80 90 100
Time (sec)

4
Pulse component (reconstructed from c2)

2
real(A2 c2)

−2

−4
0 10 20 30 40 50 60 70 80 90 100
Time (sec)

Figure 4. Example 4. Separation of a sinusoid and pulse using dual-BP with a zero-padded
DFT and STFT. (a) Data, x. (b) Sparse DFT coefficients, copt 1 . (c) Sparse STFT coefficients,
copt
2 . (b) Estimated sinusoid, A c
1 1
opt
. (c) Estimated pulse, A opt opt
2 c2 . The coefficients, c1 and
copt
2 , are obtained by solving the dual-BP problem.
16 IVAN SELESNICK

4 4
Noise−free signal (s1 + s2) y = s1 + s2 + white noise

2 2

0 0

−2 −2

−4 −4
0 20 40 60 80 100 0 20 40 60 80 100
Time (sec) Time (sec)

7 4
DFT coefficients (c1) [Output of dual−BPD] Sinusoid component (reconstructed from c1)
6
5 2

real(A1 c1)
abs( c1 )

4
0
3
2 −2
1
0 −4
0 0.1 0.2 0.3 0.4 0.5 0 20 40 60 80 100
Frequency (Hz) Time (sec)

0.5 4
STFT coefficients (c2) [Output of dual−BPD] Pulse component (reconstructed from c2)

0.4
2
Frequency (Hz)

real(A2 c2)

0.3
0
0.2

−2
0.1

0 −4
0 20 40 60 80 100 0 20 40 60 80 100
Time (sec) Time (sec)

Figure 5. Example 5. Estimation of a sinusoid and pulse in white noise using dual-BPD
with a zero-padded DFT and STFT. (a) Nois-free data. (b) Noisy data, x. (c) Sparse
DFT coefficients, copt opt opt
1 . (d) Estimated sinusoid, A1 c1 . (e) Sparse STFT coefficients, c2 .
(f) Estimated pulse, A2 copt opt
2 . The coefficients, c1 and copt
2 , are obtained by solving the
dual-BPD problem.

If Ai and AHi are fast, then these algorithms as a whole are fast. For example, if the Ai are FFTs
and/or short-time Fourier transforms, then they have low implementation complexity and can admit high
parallelism. Such a combination of transforms is useful for decomposing a signal into narrow-band and wide-
band signal components, even when the components overlap in both time and frequency. An application
of this method to radar signal processing is described in Ref. [9]. Alternately, by taking the Ai as wavelet
transforms with different Q-factors, a signal can be decomposed into low and high resonance components
[8].

Appendix A. Soft Threshold Function

The soft-thresholding function, soft : C × R+ → C, is defined as

soft(x, T ) = max 1 − T /|x|, 0 · x. (59)
 L1-NORM PENALIZED LEAST SQUARES WITH SALSA 17

soft(x, T)
10

0
T=2
−2

−4

−6

−8

−10
−10 −8 −6 −4 −2 0 2 4 6 8 10
x

Figure 6. The soft threshold function on the real line.

The soft-threshold function on the real-line is illustrated in Fig. 6.

When we apply soft-thresholding to a vector, x ∈ CN , we apply it component-wise; i.e.
[soft(x, T )]i = soft(xi , T ). (60)
If both x and T are vectors of equal length, then
[soft(x, T)]i = soft(xi , Ti ). (61)

Appendix B. Matrix Inverse Lemma

The matrix inverse lemma is given by

−1
−1
(A + BCdD) = A−1 − A−1 B C−1 + DA−1 B DA−1 . (62)
From (62), we obtain

−1 1 1
−1
µ I + AH A = I − AH µ I + AAH A. (63)
µ µ
References
[1] M. V. Afonso, J. M. Bioucas-Dias, and M. A. T. Figueiredo. Fast image recovery using variable splitting and constrained
optimization. IEEE Trans. Image Process., 19(9):2345–2356, September 2010.
[2] M. V. Afonso, J. M. Bioucas-Dias, and M. A. T. Figueiredo. An augmented Lagrangian approach to the constrained
optimization formulation of imaging inverse problems. IEEE Trans. Image Process., 20(3):681–695, March 2011.
[3] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein. Distributed optimization and statistical learning via the alternating
direction method of multipliers. Foundations and Trends in Machine Learning, 3(1):1–122, 2011.
[4] S. Chen, D. L. Donoho, and M. A. Saunders. Atomic decomposition by basis pursuit. SIAM J. Sci. Comput., 20(1):33–61,
1998.
[5] J. Eckstein and D. Bertsekas. On the Douglas-Rachford splitting method and the proximal point algorithm for maximal
monotone operators. Math. Program., 5:293–318, 1992.
[6] I. Selesnick. The short-time Fourier transform and speech denoising. Connexions, 2009. [Link]
[7] I. Selesnick. Introduction to sparsity in signal processing. Connexions, 2012. [Link]
[8] I. W. Selesnick. Resonance-based signal decomposition: A new sparsity-enabled signal analysis method. Signal Processing,
91(12):2793 – 2809, 2011.
18 IVAN SELESNICK

[9] I. W. Selesnick, K. Y. Li, S. U. Pillai, and B. Himed. Doppler-streak attenuation via oscillatory-plus-transient decomposition
of IQ data. In IET Int. Conf. Radar Systems, 2012.
[10] J.-L. Starck, M. Elad, and D. Donoho. Redundant multiscale transforms and their application for morphological component
analysis. Advances in Imaging and Electron Physics, 132:287–348, 2004.
[11] J.-L. Starck, M. Elad, and D. Donoho. Image decomposition via the combination of sparse representation and a variational
approach. IEEE Trans. Image Process., 14(10):1570–1582, October 2005.