Scribd
Upload
11 views3 pages

C Implementation of Sparse Signal Algorithms

The document outlines an algorithm for sparse signal representation using basis pursuit and ADMM in C, detailing the steps for optimization while managing memory constraints. It discusses the importance of avoiding overwriting vectors during computations and presents a modified algorithm that introduces an additional variable to preserve necessary data. Additionally, it mentions the potential for efficiency improvements by utilizing a unitary matrix to reduce computational complexity in the update equations.

Uploaded by

scribd
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
11 views3 pages

C Implementation of Sparse Signal Algorithms

The document outlines an algorithm for sparse signal representation using basis pursuit and ADMM in C, detailing the steps for optimization while managing memory constraints. It discusses the importance of avoiding overwriting vectors during computations and presents a modified algorithm that introduces an additional variable to preserve necessary data. Additionally, it mentions the potential for efficiency improvements by utilizing a unitary matrix to reduce computational complexity in the update equations.

Uploaded by

scribd
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd

1 Algorithm implementation in C

1.1 Sparse representation (with equality)


Consider the problem of finding a sparse signal representation of a signal y in a tight frame Φ. One approach
is ‘basis pursuit’ which finds the coefficients a by minimization of the `1 norm:

argmin kλ ak1 (1)


a

such that Φa = y (2)

where the columns of Φ form a tight frame:


Φ Φt = c I.

Using ADMM, an algorithm for solving the optimization problem is:

initialize: µ > 0, d (3)


u ← soft(a + d, 0.5λ/µ) − d (4)
1
d ← Φt (y − Φu) (5)
c
a←d+u (6)
repeat (7)

In our C implementation, when we compute Φu we overwrite (destroy) the contents of u. Therefore, we


can not use u again on the ‘right-hand-side’ of an assignment — for example, we can not use u in (6) because
it was overwritten in (5). The iteration above needs to be modified to account for this. One way to address
this issue is to use save u in a prior to computing Φu in (5). This way we have the following iteration:

initialize: µ > 0, d (8)


u ← soft(a + d, 0.5λ/µ) − d (9)
a←u (10)
1
d ← Φt (y − Φu) (11)
c
a←d+a (12)
repeat (13)

To save memory, we can reuse u in place of d:

initialize: µ > 0, u (14)


u ← soft(a + u, 0.5λ/µ) − u (15)
a←u (16)
1
u ← Φt (y − Φu) (17)
c
a←u+a (18)
repeat (19)

1
In our C implementation, the use of Φt also overwrites the vector to which it is applied. We do not want
to overwrite y in (17) because y is needed on every iteration. To avoid this, note that the step
1 t
u← Φ (y − Φu)
c
can be implemented as

z ← Φu (20)
z←y−z (21)
1
u ← Φt z (22)
c
where we have introduced an additional variable z. This way we avoid overwriting y.
The full algorithm is now:

initialize: µ > 0, u (23)


u ← soft(a + u, 0.5λ/µ) − u (24)
a←u (25)
z ← Φu (26)
z←y−z (27)
1
u ← Φt z (28)
c
a←u+a (29)
repeat (30)

For this version of the iteration it does not matter that Φ and Φt overwrites the vectors to which they are
applied. It is a ‘C-friendly’ algorithm description. Also, in addition to a and y of the original problem
statement, only two additional vectors need to be allocated: u and z. For u, we need to allocate a set of
subbands. For z, we need to allocate a signal.

1.2 Sparse approximation


The same procedure can be used for the case of sparse approximation (for example, basis pursuit denoising):

argmin ky − Axk22 + kλ xk1


x

where the columns of A form a tight frame:

A At = c I.

1.3 MCA
The same procedure can be used for MCA.

2
1.4 Further efficiency
Consider the update equation (5)
1 t
d← Φ (y − Φu).
c
For the some transform Φ we may have
Φ = QΦ̂

where Q is a unitary matrix


Qt Q = QQt = I

and Φ̂ has lower computational complexity that Φ. (For the TQWT, Q is the Fourier transform on the
‘signal side’ of the TQWT.) Then
Φt = Φ̂t Qt

and
Φt Φ = Φ̂t Qt QΦ̂ = Φ̂t I Φ̂ = Φ̂t Φ̂

and equation (5) can be written as


1 t t 1 1
d← Φ̂ Q (y − QΦ̂u) = Φ̂t (Qt y − Qt QΦ̂u) = Φ̂t (y − Φ̂u)
c c c
where
y := Qt y

can be precomputed. So in the iterative algorithm, we can replace Φ with Φ̂ and we can replace Φt with Φ̂t .

You might also like