The ADMM algorithm for audio signal recovery and performance modification with the dual Douglas-Rachford dynamical system

1.   Introduction

Throughout, $\mathcal{H}$ denotes a Hilbert space. Additionally, $\Gamma_{0}(\mathcal{H})$ denotes the set of proper lower semicontinuous convex functions from $\mathcal{H}$ to $[\infty_{-}, \infty_{+}]$. For readers unfamiliar with these notions, it suffices to know that the functions $f, g, f^*, g^*$ from the problems described herein all belong to the class $\Gamma_0(\mathbb{R}^r)$.

The alternating direction method of multipliers (ADMM) solves problems of the form

where $f$ and $g$ are convex functions, $\nu \in \mathbb{R}^r$, $J \in \mathbb{R}^{r \times n}$, and $B \in \mathbb{R}^{r \times m}$. The ADMM algorithm deals with the two functions $f$ and $g$ separately, and so it is said to be an operator splitting method. Depending upon how we choose the functions $f$ and $g$, we can model many different problems in this very flexible form. These include regression and classification problems, image and signal denoising, probabilistic optimization problems like progressive hedging, and basis pursuit ^[13]. We can use the latter to recover an approximate audio signal from a compressed set of sampling data.

This paper is outlined as follows. In Section 2 we introduce the basis pursuit audio problem. In Section 3 we introduce ADMM, motivating with basis pursuit. In Section 4, we describe the Douglas-Rachford method (DR) and use the basis pursuit problem to illustrate DR's duality-like relationship with ADMM. In Section 5, we describe the exemplar Lyapunov surrogate method with which we will work and provide (Theorem 2), a clean way to compute it. The remainder of the paper contains our experimental results, which we summarize in Section 7. This paper is intended to be widely accessible, including to researchers who are unfamiliar with ADMM. As such, we adopt a somewhat tutorial approach to the topic.

2.   Basis pursuit for audio signal recovery

When digital audio is recorded, a microphone's diaphragm vibrates in response to sound waves, and $N$ measurements are taken of the diaphragm's displacement. A typical rate at which these samples are collected is 44,100 per second (44.1 kHz). Recording over a time interval of length $L$, we denote by $v_i$ the measurement value at sample time:

We can define a continuous function

The unique constants $x_k$ can be computed as follows. First, writing down the equality for the $i$th sample yields

Since we have $N$ such equalities (one for each sample $(t_i, v_i)$) in $N$ unknowns ($x_0, \dots, x_{N-1}$), we can write them as the following linear system:

The matrix $M$ is the inverse discrete cosine transform matrix, and it is orthogonal. Its inverse $M^T$—the discrete cosine transform matrix, allows us to compute the $x$ values as $x = M^T v$. The nonzero values in $x$ determine the shifts and dilates of the cosine function (i.e., the frequencies) that are present in the signal $v$. For this reason, $x$ is said to be the spectrum of $v$.

2.1.   Signal compression

To reduce the size of audio files, we can delete some of the samples. The deletion of a sample $(t_{i}, v_{i})$ is reflected in the system by the removal of row $i$ from $M$ and $v$. For example, removing samples $(t_2, v_2)$ and $(t_4, v_4)$, the previous linear system becomes

Of course, the solution $\tilde{x}$ is no longer unique. The set of solutions is a subspace of dimension $2$, since $2$ is the number of samples we have deleted and therefore the number of free variables we have gained. We can still use the Moore-Penrose pseudo-inverse $\tilde{M}^T$ to obtain a spectrum $\tilde{x} = \tilde{M}^T \tilde{v}$. The signal $\tilde{V}$ that corresponds to this spectrum $\tilde{x}$ still matches the samples we have not thrown away. Assuming we have not thrown away many samples, $\tilde{V}$ may still provide a good approximation of $V$. In practice, of course, we want to reduce file size appreciably. When we remove 90% of samples (as is typical), the spectrum $\tilde{M}^T\tilde{v}$ no longer provides a good approximation signal. Note that while we removed entries 2 and 4 for illustrative purposes, there is nothing special about removing even-numbered entries; in our experiments, we remove 90% of entries randomly.

2.2.   Basis pursuit for signal recovery

We would expect that the sparsest solution $\tilde{x}$ consists of the frequencies that actually matter, while avoiding the superfluous ones. Finding the sparsest solution to a linear system $Ax = c$ is equivalent to the optimization problem:

This is called the basis pursuit problem. For our signal recovery application, we solve this optimization problem with $A = \tilde{M}$ and $c = \tilde{v}$. Using a Matlab code (P-signal) to generate a signal $c$, we compare in Figure 1^* the solution from (2.2) to a solution obtained using the Moore-Penrose pseudo-inverse. For this problem, $c$ has $15$ frequencies (i.e., Freq = $15$) denoted by $b_1, \dots, b_{15}\in \{1, 2, \dots, 500\}$ with corresponding amplitudes $a_1, \dots, a_{15} \in \left[0, 1\right]$. Here $c$ represents 1/10th of a second of audio, and so is defined over domain $d = [0, \frac{1}{44100}, \frac{2}{44100}, ..., \frac{4410}{44100}]$.

^*The specific signal of Figure 1 can be generated by fixing the random generators with the additional code rng(3).

3.   The alternating direction method of multipliers

We can use ADMM to solve the basis pursuit problem. We first define the following set and associated function:

We can use these to rewrite (2.2) in the following form:

The function $\delta_S$ is called the indicator of the set $S$. We introduce it so that we can take advantage of ADMM's ability to treat the functions $f(x) = \delta_S(x)$ and $g(z) = \|z\|_1$ separately. ADMM solves the problem

by finding a triple $(x^{*}, z^{*}, \lambda^{*})$. This triple is a saddle point for $L_{0}$ in the sense that $(x^*, z^*)$ minimizes $(x, z)\mapsto L_0(x, z, \lambda^*)$ while $\lambda^*$ maximizes $\lambda \mapsto L_0(x^*, z^*, \lambda)$. See, for example, [9, Section 5.5] and [10, Section 4.3]. Here $\rho$ is a chosen constant. The multiplier variable $\lambda$ generalizes the classical notion of a Lagrange multiplier for a general convex optimization problem, and is also called a dual variable. The primal variables $x$ and $z$ have the property that the optimal pair $(x^{*}, z^{*})$ is a solution to (BP-p). The ADMM update steps are

The multiplier update (3.2c) maximizes a concave quadratic approximation to $\lambda \mapsto L_\rho(x^{k+1}, z^{k+1}, \lambda)$. ADMM draws its name from this iterative process of minimizing with respect to the primal variables $x$ and $z$, and then maximizing over the dual variable $\lambda$.

3.1.   ADMM updates for basis pursuit

Let us look at how these updates are computed for the basis pursuit problem. We begin with the $x$ update:

Here $P_S$ is the closest-point projection operator for the set $S$. As $S$ is closed and convex in a Hilbert space, $P_{S}$ has nonempty, single-valued images ^[3]. Moreover, since $S$ is an affine subspace defined by $\tilde{M} x = c$, $P_{S}$ has a closed form ^[11]. We can use that closed form to write (3.3) explicitly:

We compute the $z$ update as follows:

The operator $\mathcal{S}_{\frac{1}{\rho}}(u)$ is known as the soft thresholding function. It is computable as follows:

Both the projection operator $P_S$ and the shrinkage operator $\mathcal{S}_{\frac{1}{\rho}}$ are examples of proximity operators. The proximity operator for a function $h \in \Gamma_{0}(\mathcal{H})$ is defined as follows:

Straight from the definition, one sees that the projection operator $P_S$ is ${\rm prox}_{\delta_S}$, while the soft thresholding operator $\mathcal{S}_{\frac{1}{\rho}}$ is ${\rm prox}_{\| \cdot \|_1}$. Such operators are used to compute the steps for operator splitting methods.

4.   ADMM and DR duality

We will next look at another operator splitting method, one that has a special duality-like relationship with ADMM. Let $I$ be the identity map and $R_{h}$ denote the reflected proximal mapping:

The Douglas-Rachford algorithm (DR) minimizes $f+g$ where $f, g \in \Gamma_{0}(\mathcal{H})$ by generating a sequence $(y^{n})_{n = 1}^{\infty}$ as

In ^[2], Attouch showed^† that if $f$ and $g$ are convex and $0$ belongs to the interior of $\{x: f(x) < \infty \} - \{x: g(x) < \infty \}$, then $f + g$ satisfies the necessary conditions for the classical DR convergence result of Lions and Mercier.

^†Attouch's actual requirements are less restrictive than what we have presented.

Theorem 1. (^[2,26]) If $f, g\in \Gamma_{0}(\mathcal{H})$ and Attouch's criterion holds, then $y^k$ from (4.2) converges weakly to some $y \in \mathcal{H}$ such that ${\rm prox}_{g} y$ minimizes $f + g$.

For our context, we are interested in the setting where the Douglas-Rachford method is applied to the problem

where $A$ is a matrix of appropriate dimension. Here $h^*$ denotes the Fenchel-Moreau-Rockafellar conjugate of $h$ ^[3,24], defined as

We are particularly interested in the case when $A = I$, $f^* = \delta_S^*$, and $g^* = \| \cdot \|_1^*$. The conjugates are

Note that Attouch's condition is satisfied for $\delta_S^*$ and $\|\cdot\|_1^*$, as long as the system $Ax = b$ is consistent. The general problem (d) is said to be dual to the general problem (p), while the specific problem

is dual to the primal basis pursuit problem (BP-p). For fully convex problems, DR and ADMM are closely related, with ADMM for (p) being equivalent to DR for (d) ^{[19,20,24,25]}. The linear relationship between DR and ADMM variables is outlined in Table 1.

From Theorem 1 and the relationships in the table, it is clear that the multiplier sequence $\lambda^k$ converges to the solution of (d). For the basis pursuit problem (BP-p) and its dual problem (BP-d), the relationship in Table 1 is illustrated in Figure 2. Because the DR operator $T_{g^*, f^*}$ has particularly desirable properties, DR convergence results are often cited as particularly elegant proofs of the convergence of ADMM, and the two algorithms are frequently studied together ^[22,23].

5.   A Lyapunov surrogate method for basis pursuit

Benoist, Dao, Tam, Giladi, and Rüffer used Lyapunov functions to describe the dynamics of DR in cases when it exhibited spiraling patterns similar to the one in Figure 2 ^[8,14,21]. All of their Lyapunov functions $U$ shared the property that

Borrowing inspiration from their constructions, the method introduced in ^[24] provides an update candidate that minimizes a spherical surrogate for such a Lyapunov function. Figure 2 shows an example of such a surrogate, and the update obtained from it. When $y^k, y^{k+1}, y^{k+2}$ are collinear, the construction is not possible and so a normal algorithmic update would be accepted. Otherwise, this update, denoted $L_T(y^k)$, is constructed to belong to the 2-dimensional affine subspace containing $y^k$, $y^{k+1}$, and $y^{k+2}$, and to have the property that $\langle L_Ty_k - y^{k+1}, y^{k+1} - y^k \rangle = \langle L_Ty^k - y^{k+2}, y^{k+2} - y^{k+1} \rangle = 0$. A candidate update for the multiplier is

We can propagate the update to the primal variables by $z^{L_{T}} \gets {L_{T}}y^k - \lambda^{L^{T}}$ and $x^{L_T} \gets P_S(z^{L_T}-\lambda^{L_T})$. Then $x^{L_T}$ is a candidate for updating $x$. There are various criteria we can use to decide whether to accept it or reject it in favor of a regular ADMM update. In this paper, the main criterion we use in our experiments is whether $x_{L_T}$ has a smaller 1-norm value than a regular update would. Based on the symbol for the operator $L_T: \mathbb{R}^r \rightarrow \mathbb{R}^r$, we refer to the full algorithm in this case as LT. For comparison, in some experiments we instead accept the $x_{L_T}$ update always. We denote the algorithm in such cases by LTA. The operator $L_T$ may be computed efficiently, as we now show.

Theorem 2. Let $y^0, y^1, y^2 \in \mathcal{H}$. Let $w^1 = y^1-y^0$ and $w^2 = y^2-y^0$. If $\|w^1\|^2\|w^2\|^2-(\langle w^1, w^2\rangle)^2 = 0$, then $y^0, y^1, y^2$ are collinear. Otherwise, let $w = y^0+\mu_1w^1+\mu_2w^2$ where

and $w$ is the unique point in the 2-dimensional affine subspace of $\mathcal{H}$ containing $\{y^0, y^1, y^2\}$ that satisfies

Proof. Letting $u = w-y^0$, the linear system rewrites to

The second equality yields $\langle u, w^1 \rangle = \|w^1\|^2$, whereupon the first equality yields

As $w$ belongs to the 2-dimensional affine subspace containing $\{y^0, y^1, y^2\}$, there must exist $\mu_1, \mu_2 \in \mathbb{R}$ such that $u = \mu_1w^1 + \mu_2 w^2$. Substituting this identity into the two equalities above yields the system:

The left matrix is invertible so long as its determinant $\delta$ is nonzero. Bearing in mind that $\|w^1\| \|w^2\| \cos(\theta) = \langle w^1, w^2 \rangle$ where $\theta$ is the angle between $w^1$ and $w^2$, we have that this determinant is zero exactly when the angle between $w^1$ and $w^2$ is $0$ or $\pi$ (i.e., $y^0, y^1, y^2$ are collinear). When the determinant is nonzero, we take the inverse of the matrix on the left and obtain the claimed identity for $\mu$. □

Remark 1 (Checking collinearity). In computational practice, to decide whether the determinant $\delta$ in Theorem 2 is nonzero, one must choose some $\epsilon > 0$ and replace the condition $\delta \neq 0$ with $|\delta| > \epsilon$. In our experiments, we used the threshold $\epsilon = 10^{-10}$, and never encountered the case when $|\delta| < \epsilon$. If $\epsilon$ is chosen very small and the condition number is very large, the computed $\mu$ may suffer numerical inaccuracy. This is another important reason for the inclusion of the objective function check.

Interestingly, for structured feasibility problems of the kind studied in ^[8,14,21], it was shown in ^[24] that the circumcentered reflections method (see, for example, ^{[1,4,5,6,7,15,17]}) is another example of an algorithm that returns minimizers for surrogates of a Lyapunov function for the DR dynamical system. For reasons described in ^[24], the natural generalization of circumcentered reflections method does not work for basis pursuit problems, and so we do not include it here. Indeed, a key motivation for the algorithm LT is that whenever a Lyapunov function has the property (5.1), LT and LTA return minimizers of spherical surrogates for that Lyapunov function ^[24]. This means that they may be useful for a wider variety of applications than those listed here. Algorithm 1 shows their implementation for basis pursuit specifically. In Table 2, we summarize the criteria that we use in our experiments to decide whether to accept an update candidate derived from $L_{T}(y^k)$ or reject it in favor of a regular update.

6.   Experimental results

Performance profiles are a popular way of comparing algorithm performance; see ^[18], whose notation we closely follow. We use them to compare the performances of LT, LTA, and ADMM. For each problem $p$ and solver $s$, we define

The performance of each algorithm is evaluated relative to the superior algorithm in each instance. This is done according to the following ratio, where $S$ denotes the set of solvers

Our performance profile graphs depict the cumulative distribution functions $\Phi_{s}$, defined as

where $n_{p}$ is the total number of problems solved by $s$. For example, if $\Phi_{s}(1.5) = 0.8$, then algorithm $s$ solved 80% of problems using no more than 150% of the passes through (3.2) required by the algorithm that used the fewest passes.

6.1.   Random problems

In ^[24], LT was found to reliably outperform ADMM for problem (2.2), where the problems were randomly generated according to the specifications in Boyd, Parikh, Chu, Peleato and Eckstein's Matlab scripts ^[12]. The code for these problems is (P-random). Here, $A \in \mathbb{R}^{500 \times 5000}$, with each entry drawn from the standard normal distribution, while $x \in \mathbb{R}^{5000}$ is sparse with approximately $500$ nonzero $x_{i} \sim \mathcal{N}(0, 1)$. The data vector $c$ is given by $c = Ax$. The performance profile for LT, LTA, and ADMM in solving 1000 problems of this kind is Figure 3b. For these problems, the algorithms stopped once locating solutions within $\pm 10^{-4}$ of optimality (i.e., $\epsilon_{\text{abs}} = \epsilon_{\text{rel}} = 10^{-4}$).

Following ^[24], we say that an algorithm shows signs of spiraling when $\|{\rm prox}_{g^*}y^{k+1} - {\rm prox}_{g^*}y^{k}\|$ (i.e., $\|\lambda^{k+1}-\lambda^k\|$) exhibits a distinct "tombstone" pattern (observable in Figure 3a). This phenomenon is often associated with the existence of a Lyapunov function satisfying condition (5.1). Problems generated using (P-random) typically exhibit signs of spiraling.

In Figure 4, we illustrate the behavior of the DR dual sequence $(y^k)_{k \in \mathbb{N}}$ for a problem generated using (P-random)^‡. We visualize $(y^k)$'s behavior by projecting it onto 10 different 2-dimensional subspaces (i.e., plotting coordinate pairs). Namely, we plot ($y_1^n, y_2^n$) (red), ($y_1^n, y_3^n$) (neon green), ($y_1^n, y_4^n$) (blue), ($y_1^n, y_5^n$) (cyan), ($y_2^n, y_3^n$) (magenta), ($y_2^n, y_4^n$) (yellow), ($y_2^n, y_5^n$) (purple), ($y_3^n, y_4^n$) (teal), ($y_3^n, y_5^n$) (orange), and ($y_4^n, y_5^n$) (green). In Figure 4b and 4c, we zoom in on the tail of the sequence ($y_2^n, y_3^n$) and plot the sequence ($(y_{LT}^n)_2, (y_{LT}^n)_{3}$) of coordinate pairs (blue shades) as well as the sequence ($(y_{LTA}^n)_{2}, (y_{LTA}^n)_{3}$) (green shades). The accepted $L_T$ updates are shown as squares.

^‡For the example shown, the random seed state was rng(8).

Because LT and LTA afforded measurable improvement for solving basis pursuit problems with random structure, we next benchmarked their performance on the audio compression and recovery problem.

6.2.   Audio signal recovery

To evaluate the performance of LT and LTA for audio recovery problems, we generated signals using (P-signal). $d = [0, \frac{1}{44100}, \frac{2}{44100}, ..., \frac{4410}{44100}]$ was used as the domain of each $c$, while $c$ is the sum of shifted and dilated sine functions with amplitudes and frequencies specified by $a_{1}, ..., a_{\text{Freq}} \in [0, 1]$ and $b_{1}, ..., b_{\text{Freq}} \in [1, 2, ..., 500]$. The number of frequencies (Freq) corresponds to the number of such functions composing $c$, with a higher number of frequencies resulting in more complex signals. To assess the algorithms on a variety of signals, we ran experiments with $\text{Freq} = 2, 5$ and $10$. For each Freq value, 1000 problems were solved where $A = \tilde{M}^{441\times 4410}$. We include the performance profiles for these problems in Figure 5.

We provide a sense of how the algorithms pursue a sparse solution in Figure 6, where we define the superfluous $x_{i}$ as having absolute values less than $0.7$. This threshold was determined based on plotting the spectra for various solved problems (as in Figure 1). Figures 6a–c depict $Mx^{5}$, $Mx^{30}$ and $Mx^{*}$ ($x^{*} = x^{158}$), as produced by ADMM solving (2.2) with $A = \tilde{M}$ and $c$ specified by (P-signal) with $\text{Freq} = 15$^§. We graph the percentage of $x_{i}^{k}$ that satisfy $|x_i^k| > 0.7$ for $k = 5, 30,158$, in Figure 6d. These signals were reconstructed with $\epsilon_{\text{abs}} = \epsilon_{\text{rel}} = 10^{-3}$ accuracy. Although not perfect, we found that lower $\epsilon_{\text{abs}}$ and $\epsilon_{\text{rel}}$ values did not result in better approximations for $V$, so this accuracy is used throughout our experiments.

^§For the example shown, the random seed state was rng(3).

6.2.1.   Experiments on real audio

We likewise applied LT and LTA to recovering the waveform signal of a performance of Frédéric Chopin's Scherzo No. 2. We decomposed the $2$ minute and $14$ second performance into $1340$ problems, each recovering $0.1$ seconds of audio. Using the file's sampling rate of $44800$ Hz, each problem was formulated with $A = \tilde{M}^{448 \times 4480}$. The performance profile for this experiment is Figure 7.

6.3.   Experiments: Discussion

The meta algorithms exhibited noticeably less improvement relative to ADMM when the random problem structure of (P-random) was replaced by the structure of the simulated audio problem (P-signal), and even less improvement for a real audio problem. For problems we solved using $\tilde{M}$, we typically did not observe signs of spiraling for the dual sequence $y^k$; when we did observe signs of spiraling, we only observed it after the signal was recovered (i.e. when computing more than necessary). Figure 8 shows $y^k$ coordinate pairs for an example problem generated using (P-signal)^¶ with Freq = 10. Enhanced views of sequences ($y^{n}_{1}, y^{n}_{3}$) and ($y^{n}_{1}, y^{n}_{4}$) for this problem are included in Figure 9a. Figure 9b enhances the sequence ($y^{n}_{2}, y^{n}_{3}$), and also shows the sequence of coordinate pairs $((y_{LT}^n)_2, (y_{LT}^n)_3)$ generated by LT, as well as the sequence of coordinate pairs $((y_{LTA}^n)_2, (y_{LTA}^n)_3)$ generated by the variant LTA. For the latter two sequences, regular DR updates are circles while updates accepted from the surrogate are squares.

^¶For the example shown, the random seed state was rng(4)

7.   Conclusions

Our results highlight the importance of testing algorithms on problems with real data. While the Lyapunov surrogate variant LTA is state-of-the-art (much faster than DR and equally reliable) for the real world problem of finding wavelets ^[16], its performance for the real world audio compression and recovery problem is less impressive than experiments on random problems might suggest.

Additionally, our results are consistent with the hypothesis in ^[24] that Lyapunov surrogate methods are more likely to deliver improvement in situations when signs of spiraling are present, as we suspect that LT's inconsistent performance is due to differences in the dynamics of the underlying DR dual algorithm.

Finally, the visualizations in Figures 4 and 9 are interesting in two ways. First, they support the hypothesis that signs of spiraling do, in fact, reflect the dynamic they purport to. The dynamics we observe for the simulated audio problem in Figure 9, when signs of spiraling were not present, appear much less regular than those we observed for the fully random problem in Figure 4b, when signs of spiraling were present. Second, the latter dynamics still exhibit patterns that may lend themselves to extrapolation methods or other dynamics-based schemes. We suggest this as future research.

Use of AI tools declaration

The authors declare that they have not used artificial intelligence tools in the creation of this article.

Conflict of interest

The authors declare no conflict of interest.

Data Availability

The code used to generate all examples in this paper is available at https://github.com/AndrewCalcan/Basis-Pursuit-via-LT-LTA-ADMM.git.