# Alternating Direction Method of Multipliers for Separable Convex Optimization of Real Functions in Complex Variables

Authors: Lu Li; Xingyu Wang; Guoqiang Wang
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The alternating direction method of multipliers (ADMM) has been widely explored due to its broad applications, and its convergence has been gotten in the real field. In this paper, an ADMM is presented for separable convex optimization of real functions in complex variables. First, the convergence of the proposed method in the complex domain is established by using the Wirtinger Calculus technique. Second, the basis pursuit (BP) algorithm is given in the form of ADMM in which the projection algorithm and the soft thresholding formula are generalized from the real case. The numerical simulations on the reconstruction of electroencephalogram (EEG) signal are provided to show that our new ADMM has better behavior than the classic ADMM for solving separable convex optimization of real functions in complex variables.


---

## Body

## 1. Introduction

The augmented Lagrangian methods (ALMs) are a certain class of algorithms for solving constrained optimization problems, which were originally known as the method of multipliers in 1969 [1], and were studied much in the 1970s and 1980s as a good alternative to penalty methods. They have similarities to penalty methods in that they replace a constrained optimization problem by a series of unconstrained problems and add a penalty term to the objective. In particular, a variant of the standard ALMs that uses partial updates (similar to the Gauss-Seidel method for solving linear equations) known as the alternating direction method of multipliers (ADMM) gained some attention [2]. The ADMM has been extensively explored in recent years due to broad applications and empirical performance in a wide variety of problems such as image processing [3], applied machine learning and statistics [4], sparse optimizations, and other relevant fields [2]. Specifically, an advantage of the ADMM is that it can handle linear equality constraint of the form x,z∣Ax+Bz=c, which makes distributed optimization by variable splitting in a batch setting straightforward. Recently, the convergence rates of order O(1/k) for the real case are considered under some additional assumptions; see, for example, [5–10]. For a survey on the ALMs and the ADMM, we refer to the references [1, 2, 11–16].Compressed sensing (CS) is a signal processing technique for efficiently acquiring and reconstructing a signal by finding solutions to underdetermined linear systems. In the CS processing, the sparsity of a signal can be exploited to recover it from samples far fewer than required by the Shannon-Nyquist sampling theorem. The idea of CS got a new life in 2006 when Candès et al. [17] and Donoho [18] gave important results on the mathematical foundation of CS. This methodology has attached much attention from applied mathematicians, computer scientists, and engineers for a variety of applications in biology [19], medicine [20], and radar [21], and so forth. Algorithms for signal reconstruction in a CS framework are expressed as sparse signal reconstruction algorithms. One of the most successful algorithms, known as basis pursuit (BP), is on the basis of constrained l1-norm minimization [22]. Most of the work is focused on the optimization in the real number field.Signals in complex variables emerge in many areas of science and engineering and have become the objects of signal processing. There have been many works on the signal processing in complex variable. For example, independent component analysis (ICA) for separating complex-valued signals has found utility in many applications such as face recognition [23], analysis of functional magnetic resonance imaging [24], and electroencephalograph [25]. Taking impropriety and noncircularity of complex-valued signals into consideration, the right type of processing can give significant performance gains [26]. Methods of digital modulation schemes that produce improper complex signals have been studied in [27], such as binary phase shift keying and pulse amplitude modulation. In these researches, most nonlinear optimization methods use the first-order or second-order approximation of the objective function to create a new step or a descent direction, where the approximation is either updated or recomputed in every iteration. Unfortunately, all these functions do not satisfy the Cauchy-Riemann conditions. There exists no Taylor series of f at z0 so that the series converges to f(z) in a neighborhood of z0. A common solution is to convert the optimization problem to the real domain by referring to f as a function of the real and imaginary parts of z. Reformulating an optimization problem which is inherently complex to the real domain would miss important insights on the structure of the problem that might otherwise be exploited [28]. Even so, there are many primal-dual optimization methods for optimization problems with the complex variable in the literatures. The usual method analyzing complex-valued optimization problem is to separate it into the real part and the imaginary part and then to recast it into an equivalent real-valued optimization problem by doubling the size of the constraint conditions; see [29–31] and the references therein. For some other related references on optimization problems in complex variables, we refer to [28, 32].To overcome the above-mentioned difficulties, the purpose of the paper is to generalize ADMM for separable convex optimization in the real number domain to the complex number domain. The concepts of convex function and subgradient are expanded from the real field to the complex field. By introducing the augmented complex variable, the definition of the augmented Lagrange function in complex variables is given. Under some mild assumptions, we establish the convergence of the proposed method. For the applications, we consider the BP algorithm which concludes projection algorithm and the soft thresholding operator in the complex field. Some numerical simulation results are reported to show that the proposed algorithm is indeed more efficient and more robust.The outline of the paper is as follows. In Section2, we recall some elementary theories and methods of the complex analysis and Wirtinger calculus. The ADMM for complex separable convex optimization and its convergence are presented in Section 3. In Section 4, we study the BP algorithm for the equality-constrained l1 minimization problem in the form of ADMM. In Section 5, some numerical simulations are provided. Finally, some conclusions are drawn in Section 6.


## 2. Preliminaries

In this section, we first give some notations used. Vectors are denoted by lower case, for example,z, and matrices are denoted by capital letters, for example, A. The kth entry of a vector z is denoted by zk and element (i,j) of a matrix A by aij. The subscripts ·re and ·im denote the real and imaginary parts, respectively; for example, zre=Re{z} and Aim=Im{A}. The superscripts ·T,·¯,·H and ·-1 are used for the transpose, conjugate, Hermitian conjugate, and matrix inverse. The domf denotes the domain of function f. The identity matrix of order n is denoted by In. The one-norm and two-norm are denoted by ·1 and ·2, respectively. zR denotes the real composite 2n-dimensional vector; for example, zR=(zreT,zimT)T∈R2n, obtained by stacking zre on the top of zim. The notation ∂f denotes the set of all subgradients of f.


### 2.1. Wirtinger Calculus

We next recall some well-known concepts and results on the complex analysis and Wirtinger calculus which will be used in our future analysis. A comprehensive treatment of Wirtinger calculus can be found in [33, 34].Define the complex augmented vector as follows:(1)z_=zT,zHT∈C2n, which is obtained by stacking z on the top of its complex conjugate z¯. The complex augmented vector z_∈C2n is related to the real composite vector zR∈R2n as z_=JnzR and zR=1/2JnHz_, where the real-to-complex transformation(2)Jn=InjInIn-jIn∈C2n×2nis unitary up to a factor of 2; that is, JnJnH=JnHJn=2I2n. The linear map Jn is an isomorphism map from R2n to C2n and its inverse is given by 1/2JnH.Lemma 1.
LetA∈Cp×n, B∈Cp×m, and c∈Cp. Then(3)Ax+Bz=c⟺A_x_+B_z_=c_,where(4)A_=A00A¯∈C2p×2n,B_=B00B¯∈C2p×2m.Proof.
Since(5)A_x_+B_z_-c_=A00A¯xx¯+B00B¯zz¯-cc¯=Ax+Bz-cAx+Bz-c¯,then we have(6)Ax+Bz=c⟺A_x_+B_z_=c_.This completes the proof.Consider a complex-valued function(7)fz=uzre,zim+jvzre,zim, where z=zre+jzim, f:Cn→C, and u,v:Rn×Rn→R. The definition of complex differentiability requires that the derivatives be defined as the limit(8)f′z=limΔz→0⁡fz+Δz-fzΔz is independent of the direction in which Δz approaches zero in the complex plane. This requires that the Cauchy-Riemann equations (9)∂u∂zre=∂v∂zim,∂u∂zim=-∂v∂zre should be satisfied [35]. These conditions are necessary for f(z) to be complex-differentiable. A function which is complex-differentiable on its entire domain is called analytic or holomorphic. Clearly, the Cauchy-Riemann conditions do not hold for real-valued functions which are v(zre,zim)≡0, and thus cost functions are not analytic. These conditions imply complex differentiability which are quite stringent and impose a strong structure on u(zre,zim) and v(zre,zim) and, consequently, on f(z). Obviously, most cost functions do not satisfy the Cauchy-Riemann equations as these functions are typically f:Cn→R with v(zre,zim)=0.To overcome such a difficulty, a sound approach in [33] relaxes this strong requirement for differentiability and defines a less stringent form for the complex domain. More importantly, it describes how this new definition can be used for defining complex differential operators that allow computation of derivatives in a very straightforward manner in the complex number domain, by simply using real differentiation results and procedures. A function is called real differentiable when u(zre,zim) and v(zre,zim) are differentiable as the functions of real-valued variables zre and zim. Then, one can write the two real variables as zre=(z+z¯)/2 and zim=-j(z-z¯)/2 and use the chain rule to derive the operators for differentiation given in the theorem below. The key point in the derivation is regarding the two variables z and z¯ as independent variables, which is also the main approach allowing us to make use of the elegance of Wirtinger calculus.In view of this, we consider the function (7) as f:R2n→C by rewriting it as f(z)=f(u,v) and make use of the underlying R2n structure. The function f(·) can be regarded as either f(zre,zim) with variables zre and zim or f(z,z¯) with variables z and z¯, and it can be simply written as f(z). The functions may take different forms; however, they are equally valued. For convenience, we use the same function f to denote them as follows: (10)fz,z¯=fzre,zim=fz.The main result in this context is stated by Brandwood in [36].Theorem 2.
Letf:R2n→C be a function of real variables zre and zim such that f(z)=f(zre,zim), where z=zre+jzim, and that f is analytic with respect to z and z¯ independently. Then, consider the following:
(1) The partial derivatives(11)∂f∂z=12∂f∂zre-j∂f∂zim,∂f∂z¯=12∂f∂zre+j∂f∂zim can be computed by treating z¯ as a constant in f and z as a constant, respectively.
(2) A necessary and sufficient condition forf to have a stationary point is that(12)∂f∂z=0. Similarly,(13)∂f∂z¯=0 is also a necessary and sufficient condition.As for the applications of Wirtinger derivatives, we consider the following two examples, which will be used in the subsequent analysis.Example 3.
Consider the real function in complex variables as follows:(14)fx,z=2ReyHAx+Bz-c,where x∈Cn, z∈Cm, y∈Cp, c∈Cp, A∈Cp×n, and B∈Cp×m.It follows from Theorem2 that(15)∂f∂x=12∂f∂xre-j∂f∂xim=yreAre+yimAim+jyreAim-yimAre=yHA. Similarly, we have(16)∂f∂z=12∂f∂zre-j∂f∂zim=yreBre+yimBim+jyreBim-yimBre=yHB.Example 4.
Consider the real function in complex variables as follows:(17)fx,z=Ax+Bz-c22, where x∈Cn, z∈Cm, c∈Cp, A∈Cp×n, and B∈Cp×m.We have(18)fx,z=Ax+Bz-cHAx+Bz-c=rre-jrimrre+jrim=rre2+rim2, where r=Ax+Bz-c. Then (19)∂f∂x=12∂f∂xre-j∂f∂xim=12∂rre2+rim2∂xre-j∂rre2+rim2∂xim=rreAre+rimAim+jrreAim-rimAre=rHA,∂f∂z=12∂f∂zre-j∂f∂zim=12∂rre2+rim2∂zre-j∂rre2+rim2∂zim=rreBre+rimBim+jrreBim-rimBre=rHB.


### 2.2. Convex Analysis in the Complex Number Domain

In order to meet the demands of next work, we give some definitions in the complex number domain.Definition 5 (see [37]).
A setA is convex if the line segment between any two points in A lies in A; that is, if for any z1,z2∈A and any θ∈R with 0≤θ≤1, then(20)θz1+1-θz2∈A.Definition 6 (see [34]).
Letz=zre+jzim∈Cn. The complex gradient operator ∂/∂z_ is defined by (21)∂∂z_=∂∂z,∂∂z¯.The linear mapJn also defines a one-to-one correspondence between the real gradient ∂/∂zR and the complex gradient ∂/∂z¯; namely, (22)∂∂zR=JnT∂∂z_.For real function in complex variablef:Cn→R, it has an equivalent form as f(z)=u(zre,zim) according to (10). So we can similarly extend some concepts of the functions in the real number domain [38, 39] to the complex number domain.Definition 7.
A real function in complex variablef:Cn→R is convex if domf is a convex set and if for any x,y∈domf and any θ∈R with 0≤θ≤1, then(23)fθx+1-θy≤θfx+1-θfy.Definition 8.
A real function in complex variablef:Cn→R is proper if its effective domain is nonempty and it never attains -∞.Definition 9.
A real function in complex variablef:Cn→R is closed if, for each α∈R, the sublevel set x∈domf∣fx≤α is a closed set.Definition 10.
Given any real function in complex variablef:Cn→R∪+∞, a vector v∈Cn is said to be a subgradient of f at z0 if(24)fz≥fz0+2RevHz-z0.


## 2.1. Wirtinger Calculus

We next recall some well-known concepts and results on the complex analysis and Wirtinger calculus which will be used in our future analysis. A comprehensive treatment of Wirtinger calculus can be found in [33, 34].Define the complex augmented vector as follows:(1)z_=zT,zHT∈C2n, which is obtained by stacking z on the top of its complex conjugate z¯. The complex augmented vector z_∈C2n is related to the real composite vector zR∈R2n as z_=JnzR and zR=1/2JnHz_, where the real-to-complex transformation(2)Jn=InjInIn-jIn∈C2n×2nis unitary up to a factor of 2; that is, JnJnH=JnHJn=2I2n. The linear map Jn is an isomorphism map from R2n to C2n and its inverse is given by 1/2JnH.Lemma 1.
LetA∈Cp×n, B∈Cp×m, and c∈Cp. Then(3)Ax+Bz=c⟺A_x_+B_z_=c_,where(4)A_=A00A¯∈C2p×2n,B_=B00B¯∈C2p×2m.Proof.
Since(5)A_x_+B_z_-c_=A00A¯xx¯+B00B¯zz¯-cc¯=Ax+Bz-cAx+Bz-c¯,then we have(6)Ax+Bz=c⟺A_x_+B_z_=c_.This completes the proof.Consider a complex-valued function(7)fz=uzre,zim+jvzre,zim, where z=zre+jzim, f:Cn→C, and u,v:Rn×Rn→R. The definition of complex differentiability requires that the derivatives be defined as the limit(8)f′z=limΔz→0⁡fz+Δz-fzΔz is independent of the direction in which Δz approaches zero in the complex plane. This requires that the Cauchy-Riemann equations (9)∂u∂zre=∂v∂zim,∂u∂zim=-∂v∂zre should be satisfied [35]. These conditions are necessary for f(z) to be complex-differentiable. A function which is complex-differentiable on its entire domain is called analytic or holomorphic. Clearly, the Cauchy-Riemann conditions do not hold for real-valued functions which are v(zre,zim)≡0, and thus cost functions are not analytic. These conditions imply complex differentiability which are quite stringent and impose a strong structure on u(zre,zim) and v(zre,zim) and, consequently, on f(z). Obviously, most cost functions do not satisfy the Cauchy-Riemann equations as these functions are typically f:Cn→R with v(zre,zim)=0.To overcome such a difficulty, a sound approach in [33] relaxes this strong requirement for differentiability and defines a less stringent form for the complex domain. More importantly, it describes how this new definition can be used for defining complex differential operators that allow computation of derivatives in a very straightforward manner in the complex number domain, by simply using real differentiation results and procedures. A function is called real differentiable when u(zre,zim) and v(zre,zim) are differentiable as the functions of real-valued variables zre and zim. Then, one can write the two real variables as zre=(z+z¯)/2 and zim=-j(z-z¯)/2 and use the chain rule to derive the operators for differentiation given in the theorem below. The key point in the derivation is regarding the two variables z and z¯ as independent variables, which is also the main approach allowing us to make use of the elegance of Wirtinger calculus.In view of this, we consider the function (7) as f:R2n→C by rewriting it as f(z)=f(u,v) and make use of the underlying R2n structure. The function f(·) can be regarded as either f(zre,zim) with variables zre and zim or f(z,z¯) with variables z and z¯, and it can be simply written as f(z). The functions may take different forms; however, they are equally valued. For convenience, we use the same function f to denote them as follows: (10)fz,z¯=fzre,zim=fz.The main result in this context is stated by Brandwood in [36].Theorem 2.
Letf:R2n→C be a function of real variables zre and zim such that f(z)=f(zre,zim), where z=zre+jzim, and that f is analytic with respect to z and z¯ independently. Then, consider the following:
(1) The partial derivatives(11)∂f∂z=12∂f∂zre-j∂f∂zim,∂f∂z¯=12∂f∂zre+j∂f∂zim can be computed by treating z¯ as a constant in f and z as a constant, respectively.
(2) A necessary and sufficient condition forf to have a stationary point is that(12)∂f∂z=0. Similarly,(13)∂f∂z¯=0 is also a necessary and sufficient condition.As for the applications of Wirtinger derivatives, we consider the following two examples, which will be used in the subsequent analysis.Example 3.
Consider the real function in complex variables as follows:(14)fx,z=2ReyHAx+Bz-c,where x∈Cn, z∈Cm, y∈Cp, c∈Cp, A∈Cp×n, and B∈Cp×m.It follows from Theorem2 that(15)∂f∂x=12∂f∂xre-j∂f∂xim=yreAre+yimAim+jyreAim-yimAre=yHA. Similarly, we have(16)∂f∂z=12∂f∂zre-j∂f∂zim=yreBre+yimBim+jyreBim-yimBre=yHB.Example 4.
Consider the real function in complex variables as follows:(17)fx,z=Ax+Bz-c22, where x∈Cn, z∈Cm, c∈Cp, A∈Cp×n, and B∈Cp×m.We have(18)fx,z=Ax+Bz-cHAx+Bz-c=rre-jrimrre+jrim=rre2+rim2, where r=Ax+Bz-c. Then (19)∂f∂x=12∂f∂xre-j∂f∂xim=12∂rre2+rim2∂xre-j∂rre2+rim2∂xim=rreAre+rimAim+jrreAim-rimAre=rHA,∂f∂z=12∂f∂zre-j∂f∂zim=12∂rre2+rim2∂zre-j∂rre2+rim2∂zim=rreBre+rimBim+jrreBim-rimBre=rHB.


## 2.2. Convex Analysis in the Complex Number Domain

In order to meet the demands of next work, we give some definitions in the complex number domain.Definition 5 (see [37]).
A setA is convex if the line segment between any two points in A lies in A; that is, if for any z1,z2∈A and any θ∈R with 0≤θ≤1, then(20)θz1+1-θz2∈A.Definition 6 (see [34]).
Letz=zre+jzim∈Cn. The complex gradient operator ∂/∂z_ is defined by (21)∂∂z_=∂∂z,∂∂z¯.The linear mapJn also defines a one-to-one correspondence between the real gradient ∂/∂zR and the complex gradient ∂/∂z¯; namely, (22)∂∂zR=JnT∂∂z_.For real function in complex variablef:Cn→R, it has an equivalent form as f(z)=u(zre,zim) according to (10). So we can similarly extend some concepts of the functions in the real number domain [38, 39] to the complex number domain.Definition 7.
A real function in complex variablef:Cn→R is convex if domf is a convex set and if for any x,y∈domf and any θ∈R with 0≤θ≤1, then(23)fθx+1-θy≤θfx+1-θfy.Definition 8.
A real function in complex variablef:Cn→R is proper if its effective domain is nonempty and it never attains -∞.Definition 9.
A real function in complex variablef:Cn→R is closed if, for each α∈R, the sublevel set x∈domf∣fx≤α is a closed set.Definition 10.
Given any real function in complex variablef:Cn→R∪+∞, a vector v∈Cn is said to be a subgradient of f at z0 if(24)fz≥fz0+2RevHz-z0.


## 3. ADMM for Convex Separable Optimization

In this section, we will first recall the ADMM for real convex separable optimization. Then we will study the ADMM for convex separable optimization of real functions in complex variables.


### 3.1. ADMM for Real Convex Separable Optimization

The ADMM has been well studied for the following linearly constrained separable convex programming whose objective function is separated into two individual convex functions with nonoverlapping variables as follows:(25)minimizefx+gz:Ax+Bz=c,x∈χ1,z∈χ2,where χ1⊂Rn and χ2⊂Rm are closed convex sets; A∈Rp×n and B∈Rp×m are given matrices; c∈Rp is a given vector; and f:Rn→R∪{+∞} and g:Rm→R∪{+∞} are proper, closed, and convex functions.More specifically, the Lagrangian function and the augmented Lagrangian function of (25) are given by (26)L0x,z,y=fx+gz+yTAx+Bz-c,(27)Lρx,y,z=fx+gz+yTAx+Bz-c+ρ2Ax+Bz-c22, respectively. Then the iterative scheme of the ADMM for solving (25) is given by(28)xk+1=argminx⁡Lρx,zk,yk,zk+1=argminz⁡Lρxk+1,z,yk,yk+1=yk+ρAxk+1+Bzk+1-c.Without loss of generality, we give the following two assumptions.Assumption 11.
The (extended-real-valued) functionsf:Rn→R∪{+∞} and g:Rm→R∪{+∞} are proper, closed, and convex.Assumption 12.
The Lagrangian functionL0 has a saddle point; that is, there exists (x∗,z∗,y∗), not necessarily unique, for which(29)L0x∗,z∗,y≤L0x∗,z∗,y∗≤L0x,z,y∗holds for all x, z, and y.The convergence of the ADMM for real convex separable optimization is established in the following theorem.Theorem 13 (Section  3.2.1 in [2]).
Under Assumptions11 and 12, the ADMM iterates (28) satisfy the following.
(1) Residual Convergence. rk=Axk+Bzk-c→0 as k→∞; that is, the iterates approach feasibility.
(2) Objective Convergence. f(xk)+g(zk)→f(x∗)+g(z∗) as k→∞; that is, the objective function of the iterates approaches the optimal value.
(3) Dual Variable Convergence. yk→y∗ as k→∞, where y∗ is a dual optimal point.


### 3.2. ADMM for Complex Convex Separable Optimization

According to (10), we can consider the real functions in complex variables f:Cn→R∪{+∞} and g:Cm→R∪{+∞}. Then, the convex separable optimization of real functions in complex variables becomes (30)minimizefx+gz:Ax+Bz=c,x∈χ1,z∈χ2, where f and g are proper, closed, and convex functions; χ1⊂Cn and χ2⊂Cm are closed convex sets; A∈Cp×n and B∈Cp×m are given matrices; and c∈Cp is a given vector.From (10) and Lemma 1, we can conclude that the complex convex separable optimization (30) is equivalent to the following convex separable optimization problem: (31)minimizefx,x¯+gz,z¯:A_x_+B_z_=c_,x∈χ1,z∈χ2.The Lagrangian function of (31) is(32)L0x,z,y=fx,x¯+gz,z¯+y_HA_x_+B_z_-c_=fx,x¯+gz,z¯+2ReyHAx+Bz-c, where y∈Cp. Then, the augmented Lagrangian function of (31) is(33)Lρx,z,y=fx,x¯+gz,z¯+y_HA_x_+B_z_-c_+ρ2A_x_+B_z_-c_22=fx,x¯+gz,z¯+2ReyHAx+Bz-c+ρAx+Bz-c22, where ρ>0 is called the penalty parameter. The ADMM for complex convex separable optimization is composed of the iterations (34)xk+1=argminx⁡Lρx,zk,yk,zk+1=argminz⁡Lρxk+1,z,yk,yk+1=yk+ρAxk+1+Bzk+1-c.Let r=Ax+Bz-c. Then we have(35)2ReyHr+ρr22=ρr+u22-ρu22,where u=1/ρy is the scaled dual variable. Using the scaled dual variable, we can express the ADMM iterations (34) as(36)xk+1=argminx⁡fx,x¯+ρAx+Bzk-c+uk22,zk+1=argminz⁡gz,z¯+ρAxk+1+Bz-c+uk22,uk+1=uk+Axk+1+Bzk+1-c.


### 3.3. Optimality Conditions

The necessary and sufficient optimality conditions for the ADMM problem (31) are primal feasibility,(37)Ax∗+Bz∗-c=0,and dual feasibility, (38)0∈∂fx∗+y∗HA,(39)0∈∂gz∗+y∗HB.Because zk+1 minimizes Lρ(xk+1,z,yk), we have(40)0∈∂gzk+1+ykHB+ρrk+1HB=∂gzk+1+yk+1HB. This means that zk+1 and yk+1 always satisfy (39); thus attaining optimality leads to satisfying (37) and (38).Becausexk+1 minimizes Lρ(x,zk,yk), we have (41)0∈∂fxk+1+ykHA+ρAxk+1+Bzk-cHA=∂fxk+1+ykH+ρrk+1HA+ρBzk-zk+1HAor equivalently(42)ρBzk+1-zkHA∈∂fxk+1+yk+1HA.From (38), sk+1=ρ(B(zk+1-zk))HA can be viewed as a residual for the dual feasibility condition. By (37), rk+1=Axk+1+Bzk+1-c can be viewed as a residual for the primal feasibility condition. These two residuals converge to zero as the ADMM proceeds.


### 3.4. Convergence

Similar to the ADMM for separable convex optimization in the real number domain, we can establish the convergence of the ADMM for complex separable convex optimization.In this paper, we make the following two assumptions on the separable convex optimization of the real functions in complex variables.Assumption 14.
The (extended-real-valued) functionsf:Cn→R∪{+∞} and g:Cm→R∪{+∞} are proper, closed, and convex.Assumption 15.
The Lagrangian functionL0 (32) has a saddle point; that is, there exists (x∗,z∗,y∗), not necessarily unique, for which(43)L0x∗,z∗,y≤L0x∗,z∗,y∗≤L0x,z,y∗ holds for all x, z, and y.Theorem 16.
Under Assumptions14 and 15, the ADMM iterations (36) have the following conclusions.
(1) Residual Convergence. rk=Axk+Bzk-c→0 as k→∞; that is, the iterates approach feasibility.
(2) Objective Convergence. f(xk)+g(zk)→f(x∗)+g(z∗) as k→∞; that is, the objective function of the iterates approaches the optimal value.
(3) Dual Variable Convergence. yk→y∗ as k→∞, where y∗ is a dual optimal point.Proof.
Let(x∗,z∗,y∗) be the saddle point for L0 and q∗=f(x∗)+g(x∗). Then we have(44)L0x∗,z∗,y∗≤L0xk+1,zk+1,y∗. Since Ax∗+Bz∗-c=0 and L0(x∗,z∗,y∗) is equivalent to q∗, then we have (45)q∗≤qk+1+2Rey∗HAxk+1+Bzk+1-c. From Theorem 2 and Examples 3 and 4, we get(46)∂Lρ∂x=∂fxk+1,x¯k+1∂x+ykHA+ρAxk+1+Bzk-cHA. Note that Lρ is a real-valued function; then we get(47)∂Lρ∂x¯=∂Lρ∂x¯.By (36), xk+1 minimizes Lρ(x,zk,yk) for f is convex, which is subdifferentiable, and so is Lρ(x,z,y). Based on Theorem 2, the optimality condition is(48)0∈∂Lρxk+1,zk,yk=∂fxk+1,x¯k+1+ykHA+ρAxk+1+Bzk-cHA. Since(49)yk+1=yk+ρAxk+1+Bzk+1-c, we have(50)0∈∂fxk+1,x¯k+1+yk+1-ρBzk+1-zkHA, which implies that xk+1 minimizes (51)fx,x¯+2Reyk+1-ρBzk+1-zkHAx. Similarly, we may have that zk+1 minimizes (52)gz,z¯+2Reyk+1HBz. From (51) and (52), we have(53)fxk+1,x¯k+1+2Reyk+1-ρBzk+1-zkHAxk+1≤fx∗,x¯∗+2Reyk+1-ρBzk+1-zkHAx∗,gzk+1,z¯k+1+2Reyk+1HBzk+1≤gz∗,z¯∗+2Reyk+1HBz∗. From (53) and Ax∗+Bz∗-c=0, we can make the conclusion that (54)pk+1-p∗≤2Re-yk+1Hrk+1-ρBzk+1-zkH-rk+1+Bzk+1-z∗, where rk+1=Axk+1+Bzk+1-c.
Adding (45) and (54), we get (55)4Reyk+1-y∗Hrk+1-ρBzk+1-zkHrk+1+ρBzk+1-zkHBzk+1-z∗≤0. Let(56)wk=1ρyk-y∗22+ρBzk-z∗22. By following manipulation and rewriting of (55), we have (57)wk+1≤wk-ρrk+122-ρBzk+1-zk22.Rewriting the first term of (55) as (58)4Reyk+1-y∗Hrk+1-yk+ρrk+1-y∗Hrk+1+yk-y∗Hrk+1+ρrk+122+ρrk+122and then substituting(59)rk+1=1ρyk+1-yk into the first two terms in (58) give(60)2Re2ρyk-y∗Hyk+1-yk+1ρyk+1-yk22+ρrk+122. Since(61)yk+1-yk=yk+1-y∗-yk-y∗, this can be expressed as (62)1ρyk+1-y∗22-yk-y∗22+ρrk+122. Now let us regroup the remaining terms, that is, (63)ρrk+122-2Re2ρBzk+1-zkHrk+1+2ρBzk+1-zkHBzk+1-z∗, where ρrk+122 is taken from (62). Substituting(64)zk+1-z∗=zk+1-zk+zk-z∗ into the last term in (63) yields(65)ρrk+1-Bzk+1-zk22+ρBzk+1-zk22+4ρReBzk+1-zkHBzk-z∗ and substituting(66)zk+1-zk=zk+1-z∗-zk-z∗ into the last two terms in (63), we get(67)ρrk+1-Bzk+1-zk22+ρBzk+1-zk22-Bzk-z∗22. It implies that (55) can be expressed as (68)wk-wk+1≥ρrk+1-Bzk+1-zk22. To obtain (57), it suffices to show that the middle term(69)-4Reρrk+1HBzk+1-zk of the expanded right-hand side of (68) is positive. To understand this, by reviving that zk+1 minimizes g(z,z¯)+2Re{(yk+1)HBz} and zk minimizes g(z,z¯)+2Re{(yk)HBz}, we can add(70)gzk+1,z¯k+1+2Reyk+1HBzk+1≤gzk,z¯k+2Reyk+1HBzk,gzk,z¯k+2ReykHBzk≤gzk+1,z¯k+1+2ReykHBzk+1 to obtain that(71)2Reyk+1-ykHBzk+1-zk≤0. Since ρ>0, if we substitute(72)yk+1-yk=ρrk+1, we can get (57).
This means thatwk decreases in each iteration by an amount depending on the norm of the residual and on the change in z over one iteration. Since wk≤w0, it follows that yk and Bzk are bounded. Iterating the inequality above gives(73)ρ∑k=0∞rk+122+Bzk+1-zk22≤w0, implying that rk→0 and B(zk+1-zk)→0 as k→∞. From (45), we have(74)fxk,x¯k+gzk,z¯k⟶fx∗,x¯∗+gz∗,z¯∗ as k→∞. Furthermore, since wk→0 as k→∞, we have yk→y∗ as k→∞. This completes the proof.


### 3.5. Stopping Criterion

We can find that(75)-rk+1+Bzk+1-zk=-Axk+1-x∗. Substituting this into (54), we get (76)pk+1-p∗≤2Re-yk+1Hrk+1+xk+1-x∗Hsk+1.This means that when the two residuals are small, the error must be small. Thus an appropriate termination criterion is that the primal residuals rk+1 and dual residuals sk+1 are small simultaneously; that is, rk+12≤εpri and sk+12≤εdual, where εpri and εdual are tolerances for the primal and dual feasibility, respectively.


## 3.1. ADMM for Real Convex Separable Optimization

The ADMM has been well studied for the following linearly constrained separable convex programming whose objective function is separated into two individual convex functions with nonoverlapping variables as follows:(25)minimizefx+gz:Ax+Bz=c,x∈χ1,z∈χ2,where χ1⊂Rn and χ2⊂Rm are closed convex sets; A∈Rp×n and B∈Rp×m are given matrices; c∈Rp is a given vector; and f:Rn→R∪{+∞} and g:Rm→R∪{+∞} are proper, closed, and convex functions.More specifically, the Lagrangian function and the augmented Lagrangian function of (25) are given by (26)L0x,z,y=fx+gz+yTAx+Bz-c,(27)Lρx,y,z=fx+gz+yTAx+Bz-c+ρ2Ax+Bz-c22, respectively. Then the iterative scheme of the ADMM for solving (25) is given by(28)xk+1=argminx⁡Lρx,zk,yk,zk+1=argminz⁡Lρxk+1,z,yk,yk+1=yk+ρAxk+1+Bzk+1-c.Without loss of generality, we give the following two assumptions.Assumption 11.
The (extended-real-valued) functionsf:Rn→R∪{+∞} and g:Rm→R∪{+∞} are proper, closed, and convex.Assumption 12.
The Lagrangian functionL0 has a saddle point; that is, there exists (x∗,z∗,y∗), not necessarily unique, for which(29)L0x∗,z∗,y≤L0x∗,z∗,y∗≤L0x,z,y∗holds for all x, z, and y.The convergence of the ADMM for real convex separable optimization is established in the following theorem.Theorem 13 (Section  3.2.1 in [2]).
Under Assumptions11 and 12, the ADMM iterates (28) satisfy the following.
(1) Residual Convergence. rk=Axk+Bzk-c→0 as k→∞; that is, the iterates approach feasibility.
(2) Objective Convergence. f(xk)+g(zk)→f(x∗)+g(z∗) as k→∞; that is, the objective function of the iterates approaches the optimal value.
(3) Dual Variable Convergence. yk→y∗ as k→∞, where y∗ is a dual optimal point.


## 3.2. ADMM for Complex Convex Separable Optimization

According to (10), we can consider the real functions in complex variables f:Cn→R∪{+∞} and g:Cm→R∪{+∞}. Then, the convex separable optimization of real functions in complex variables becomes (30)minimizefx+gz:Ax+Bz=c,x∈χ1,z∈χ2, where f and g are proper, closed, and convex functions; χ1⊂Cn and χ2⊂Cm are closed convex sets; A∈Cp×n and B∈Cp×m are given matrices; and c∈Cp is a given vector.From (10) and Lemma 1, we can conclude that the complex convex separable optimization (30) is equivalent to the following convex separable optimization problem: (31)minimizefx,x¯+gz,z¯:A_x_+B_z_=c_,x∈χ1,z∈χ2.The Lagrangian function of (31) is(32)L0x,z,y=fx,x¯+gz,z¯+y_HA_x_+B_z_-c_=fx,x¯+gz,z¯+2ReyHAx+Bz-c, where y∈Cp. Then, the augmented Lagrangian function of (31) is(33)Lρx,z,y=fx,x¯+gz,z¯+y_HA_x_+B_z_-c_+ρ2A_x_+B_z_-c_22=fx,x¯+gz,z¯+2ReyHAx+Bz-c+ρAx+Bz-c22, where ρ>0 is called the penalty parameter. The ADMM for complex convex separable optimization is composed of the iterations (34)xk+1=argminx⁡Lρx,zk,yk,zk+1=argminz⁡Lρxk+1,z,yk,yk+1=yk+ρAxk+1+Bzk+1-c.Let r=Ax+Bz-c. Then we have(35)2ReyHr+ρr22=ρr+u22-ρu22,where u=1/ρy is the scaled dual variable. Using the scaled dual variable, we can express the ADMM iterations (34) as(36)xk+1=argminx⁡fx,x¯+ρAx+Bzk-c+uk22,zk+1=argminz⁡gz,z¯+ρAxk+1+Bz-c+uk22,uk+1=uk+Axk+1+Bzk+1-c.


## 3.3. Optimality Conditions

The necessary and sufficient optimality conditions for the ADMM problem (31) are primal feasibility,(37)Ax∗+Bz∗-c=0,and dual feasibility, (38)0∈∂fx∗+y∗HA,(39)0∈∂gz∗+y∗HB.Because zk+1 minimizes Lρ(xk+1,z,yk), we have(40)0∈∂gzk+1+ykHB+ρrk+1HB=∂gzk+1+yk+1HB. This means that zk+1 and yk+1 always satisfy (39); thus attaining optimality leads to satisfying (37) and (38).Becausexk+1 minimizes Lρ(x,zk,yk), we have (41)0∈∂fxk+1+ykHA+ρAxk+1+Bzk-cHA=∂fxk+1+ykH+ρrk+1HA+ρBzk-zk+1HAor equivalently(42)ρBzk+1-zkHA∈∂fxk+1+yk+1HA.From (38), sk+1=ρ(B(zk+1-zk))HA can be viewed as a residual for the dual feasibility condition. By (37), rk+1=Axk+1+Bzk+1-c can be viewed as a residual for the primal feasibility condition. These two residuals converge to zero as the ADMM proceeds.


## 3.4. Convergence

Similar to the ADMM for separable convex optimization in the real number domain, we can establish the convergence of the ADMM for complex separable convex optimization.In this paper, we make the following two assumptions on the separable convex optimization of the real functions in complex variables.Assumption 14.
The (extended-real-valued) functionsf:Cn→R∪{+∞} and g:Cm→R∪{+∞} are proper, closed, and convex.Assumption 15.
The Lagrangian functionL0 (32) has a saddle point; that is, there exists (x∗,z∗,y∗), not necessarily unique, for which(43)L0x∗,z∗,y≤L0x∗,z∗,y∗≤L0x,z,y∗ holds for all x, z, and y.Theorem 16.
Under Assumptions14 and 15, the ADMM iterations (36) have the following conclusions.
(1) Residual Convergence. rk=Axk+Bzk-c→0 as k→∞; that is, the iterates approach feasibility.
(2) Objective Convergence. f(xk)+g(zk)→f(x∗)+g(z∗) as k→∞; that is, the objective function of the iterates approaches the optimal value.
(3) Dual Variable Convergence. yk→y∗ as k→∞, where y∗ is a dual optimal point.Proof.
Let(x∗,z∗,y∗) be the saddle point for L0 and q∗=f(x∗)+g(x∗). Then we have(44)L0x∗,z∗,y∗≤L0xk+1,zk+1,y∗. Since Ax∗+Bz∗-c=0 and L0(x∗,z∗,y∗) is equivalent to q∗, then we have (45)q∗≤qk+1+2Rey∗HAxk+1+Bzk+1-c. From Theorem 2 and Examples 3 and 4, we get(46)∂Lρ∂x=∂fxk+1,x¯k+1∂x+ykHA+ρAxk+1+Bzk-cHA. Note that Lρ is a real-valued function; then we get(47)∂Lρ∂x¯=∂Lρ∂x¯.By (36), xk+1 minimizes Lρ(x,zk,yk) for f is convex, which is subdifferentiable, and so is Lρ(x,z,y). Based on Theorem 2, the optimality condition is(48)0∈∂Lρxk+1,zk,yk=∂fxk+1,x¯k+1+ykHA+ρAxk+1+Bzk-cHA. Since(49)yk+1=yk+ρAxk+1+Bzk+1-c, we have(50)0∈∂fxk+1,x¯k+1+yk+1-ρBzk+1-zkHA, which implies that xk+1 minimizes (51)fx,x¯+2Reyk+1-ρBzk+1-zkHAx. Similarly, we may have that zk+1 minimizes (52)gz,z¯+2Reyk+1HBz. From (51) and (52), we have(53)fxk+1,x¯k+1+2Reyk+1-ρBzk+1-zkHAxk+1≤fx∗,x¯∗+2Reyk+1-ρBzk+1-zkHAx∗,gzk+1,z¯k+1+2Reyk+1HBzk+1≤gz∗,z¯∗+2Reyk+1HBz∗. From (53) and Ax∗+Bz∗-c=0, we can make the conclusion that (54)pk+1-p∗≤2Re-yk+1Hrk+1-ρBzk+1-zkH-rk+1+Bzk+1-z∗, where rk+1=Axk+1+Bzk+1-c.
Adding (45) and (54), we get (55)4Reyk+1-y∗Hrk+1-ρBzk+1-zkHrk+1+ρBzk+1-zkHBzk+1-z∗≤0. Let(56)wk=1ρyk-y∗22+ρBzk-z∗22. By following manipulation and rewriting of (55), we have (57)wk+1≤wk-ρrk+122-ρBzk+1-zk22.Rewriting the first term of (55) as (58)4Reyk+1-y∗Hrk+1-yk+ρrk+1-y∗Hrk+1+yk-y∗Hrk+1+ρrk+122+ρrk+122and then substituting(59)rk+1=1ρyk+1-yk into the first two terms in (58) give(60)2Re2ρyk-y∗Hyk+1-yk+1ρyk+1-yk22+ρrk+122. Since(61)yk+1-yk=yk+1-y∗-yk-y∗, this can be expressed as (62)1ρyk+1-y∗22-yk-y∗22+ρrk+122. Now let us regroup the remaining terms, that is, (63)ρrk+122-2Re2ρBzk+1-zkHrk+1+2ρBzk+1-zkHBzk+1-z∗, where ρrk+122 is taken from (62). Substituting(64)zk+1-z∗=zk+1-zk+zk-z∗ into the last term in (63) yields(65)ρrk+1-Bzk+1-zk22+ρBzk+1-zk22+4ρReBzk+1-zkHBzk-z∗ and substituting(66)zk+1-zk=zk+1-z∗-zk-z∗ into the last two terms in (63), we get(67)ρrk+1-Bzk+1-zk22+ρBzk+1-zk22-Bzk-z∗22. It implies that (55) can be expressed as (68)wk-wk+1≥ρrk+1-Bzk+1-zk22. To obtain (57), it suffices to show that the middle term(69)-4Reρrk+1HBzk+1-zk of the expanded right-hand side of (68) is positive. To understand this, by reviving that zk+1 minimizes g(z,z¯)+2Re{(yk+1)HBz} and zk minimizes g(z,z¯)+2Re{(yk)HBz}, we can add(70)gzk+1,z¯k+1+2Reyk+1HBzk+1≤gzk,z¯k+2Reyk+1HBzk,gzk,z¯k+2ReykHBzk≤gzk+1,z¯k+1+2ReykHBzk+1 to obtain that(71)2Reyk+1-ykHBzk+1-zk≤0. Since ρ>0, if we substitute(72)yk+1-yk=ρrk+1, we can get (57).
This means thatwk decreases in each iteration by an amount depending on the norm of the residual and on the change in z over one iteration. Since wk≤w0, it follows that yk and Bzk are bounded. Iterating the inequality above gives(73)ρ∑k=0∞rk+122+Bzk+1-zk22≤w0, implying that rk→0 and B(zk+1-zk)→0 as k→∞. From (45), we have(74)fxk,x¯k+gzk,z¯k⟶fx∗,x¯∗+gz∗,z¯∗ as k→∞. Furthermore, since wk→0 as k→∞, we have yk→y∗ as k→∞. This completes the proof.


## 3.5. Stopping Criterion

We can find that(75)-rk+1+Bzk+1-zk=-Axk+1-x∗. Substituting this into (54), we get (76)pk+1-p∗≤2Re-yk+1Hrk+1+xk+1-x∗Hsk+1.This means that when the two residuals are small, the error must be small. Thus an appropriate termination criterion is that the primal residuals rk+1 and dual residuals sk+1 are small simultaneously; that is, rk+12≤εpri and sk+12≤εdual, where εpri and εdual are tolerances for the primal and dual feasibility, respectively.


## 4. Basis Pursuit with Complex ADMM

Consider the equality-constrainedl1 minimization problem in the complex number domain (77)minimizex1:Ax=b,x∈Cn, where A∈Cp×n is a given matrix, Rank(A)=p, and b∈Cp is a given vector.Recall that(78)x1=∑k=1nxk=∑k=1nxkre2+xkim2. Then(79)x_1=2x1. In the form of the ADMM, the BP method can be expressed as (80)minimizefx+z1:x=z,x,z∈Cn, where f is the indicator function of X=x∈Cn∣Ax=b; that is, f(x)=0 for x∈X and f(x)=+∞ otherwise. Then, with the idea in [40], the ADMM iterations are provided as follows:(81)xk+1=argminx⁡fx+ρx-zk+uk22,zk+1=argminz⁡z1+ρxk+1-z+uk22,uk+1=uk+xk+1-zk+1.The x-update, which involves solving a linearly constrained minimum Euclidean norm problem, can be written as (82)minimizexk+112xk+1-zk-uk22:Ax=b.Let x~=xk+1-(zk-uk). Then, (82) is equivalent to (83)minimizex~12x~22:Ax~=b~,where b~=b-A(zk-uk).Lemma 17 (see [41]).
The minimum-norm least-squares solution ofAx~=b~ is x~=A†b~, where A† is the Moore-Penrose inverse of matrix A.Theorem 18.
Thex-update of (81) is (84)xk+1=I-AHAAH-1Azk-uk+AHAAH-1b≔Πzk-uk.Proof.
AsA is of full row rank, its full-rank factorization is Am×n=Bm×mCm×n. Then it yields that [41](85)A†=CHCCH-1BHB-1BH. From Lemma 17, we have(86)x~=CHCCH-1BHB-1BHb~; that is, (87)xk+1-zk-uk=CHCCH-1BHB-1BHb-Azk-uk. Since B is of full rank, we can rearrange (87) to obtain(88)xk+1=I-AHAAH-1Azk-uk+AHAAH-1b. This completes the proof.If problem (82) is in the real number domain, we have(89)xk+1=I-ATAAT-1Azk-uk+ATAAT-1b, which is the same one obtained in Section 6.2 [2].Thez-update can be solved by the soft thresholding operator S in the following theorem, which is a generalization of the soft thresholding in [2].Theorem 19.
Lett=1/(2ρ). Then one has the following.
(1) Ifxk+1+uk is real-valued, that is, xk+1+uk=a, the soft thresholding operator is(90)Sta=max⁡0,a-t-max⁡0,-a-t=a-t,a>t;0,a≤t;a+t,a<-t.
(2) Ifxk+1+uk is purely imaginary, that is, xk+1+uk=bj, the soft thresholding operator is(91)Stbj=jmax⁡0,b-t-max⁡0,-b-t=b-tj,b>t;0,b≤t;b+tj,b<-t.
(3) Ifxk+1+uk=a+bj, the soft thresholding operator is (92)Sta+bj=max⁡0,a-tre-max⁡0,-a-tre+jmax⁡0,b-tim-max⁡0,-b-tim=0,a2+b2≤t;a-tre+jb-tim,a2+b2>t,a>0,b>0;a+tre+jb-tim,a2+b2>t,a<0,b>0;a-tre+jb+tim,a2+b2>t,a>0,b<0;a+tre+jb+tim,a2+b2>t,a<0,b<0,where tre=ta2/a2+b2 and tim=tb2/a2+b2.Proof.
(1) Assume thatxk+1+uk=a, the updating of z becomes minimizing the following function: (93)Fz=z1+ρz-a22=∑k=1nzrek2+zimk2+ρ∑k=1nzrek-ak2+zimk2.From Theorem 2, we have(94)∂F∂zimk=zimkzrek2+zimk2+ρzimk=0. This implies that (zim)k=0. Then (93) can be rewritten as(95)Fz=zre1+ρzre-a22=∑k=1nzrek+ρ∑k=1nzrek-ak2. When ak≥0, (zre)k should be positive. Then, we have(96)Fz=∑k=1nzrek+ρ∑k=1nzrek-ak2.It is clear that this is a simple parabola, and its results are(97)Ifak>t=12ρ,zrek=zrek-tk;if0≤ak≤t,zrek=0.When ak<0, we can get the similar results as follows:(98)Ifak<-t,zrek=zrek+tk;if-t≤ak<0,zrek=0.From the above discussion, we can complete the proof.
(2) Assume thatxk+1+uk=bj; we have(99)Fz=z1+ρz-bj22=∑k=1nzrek2+zimk2+ρ∑k=1nzrek2+zimk-bk2, and by adopting the same approach in the above (1), we can get the results.
(3) Assume thatxk+1+uk=a+bj and satisfy a2+b2>t,a>0,b>0; then(100)Fz=z1+ρz-a+bj22=∑k=1nzrek2+zimk2+ρ∑k=1nzrek-a2+zimk-bk2.It follows from Theorem 2 that (101)∂f∂zrek=zrekzrek2+zimk2+2ρzrek-ak=0,∂f∂zimk=zimkzrek2+zimk2+2ρzimk-bk=0.By resolving it, we may get (zre)k=ak-(tre)k, (zim)k=bk-(tim)k, where tre=ta2/a2+b2 and tim=tb2/a2+b2. Other cases can be discussed similarly. Thus we omit the proof here. This completes the proof.From what has been discussed above onx-update and z-update, the iteration of the BP algorithm is (102)xk+1=Πzk-uk,zk+1=S1/2ρxk+1+uk,uk+1=uk+xk+1-zk+1, where Π is projection onto x∈Cn∣Ax=b and S is the soft thresholding operator in the complex number domain.


## 5. Numerical Simulation

We give two numerical simulations with random data and EEG data. All our numerical experiments are carried out on a PC with Intel Core i7-4710MQ CPU at 2.50 GHz and 8 GB of physical memory. The PC runs MATLAB Version: R2013a on Window 7 Enterprise 64-bit operating system.


### 5.1. Numerical Simulation of the BP Algorithm with Random Data

Assume thatx∈Cn is a discrete complex signal interested. x itself is r-sparse, which contains (at most) r nonzero entries with r≪n. Select p(p<n) measurements uniformly at random matrix Ap×n via Ap×nx=b. Hence reconstructing signal x from measurement b is generally an ill-posed problem which is an undetermined system of linear equations. However, the sparsest solution can be obtained by solving the constrained optimization problem: (103)minimizex0:Ax=b,x∈Cn,where x0 is the l0-norm of x. Unfortunately, (103) is a combinatorial optimization problem of which the computational complexity grows exponentially with the signal size n. A key result in [17, 18] is that if x is sparse, the sparsest solution of (103) can be obtained with overwhelming probability by solving the convex optimization problem (77). Next, we consider CS which is actually a kind of application of the BP method in complex variables.


#### 5.1.1. The Effects of the Parameterρ

We demonstrate the complex signal sampling and recovery techniques with a discrete-time complex signalx of length n=300 with sparsity r=30 which is generated randomly. Ap×n is a random sensing matrix with p=120. The variables u0 and z0 are initialized to be zero. We set the two tolerances of primal and dual residuals equal to 10-6. In order to understand the effects of the parameter ρ on the convergence, we set the penalty parameter ρ from 0.1 to 20 with the step 0.1. We have repeated the same experiment 100 times with the same parameter ρ. The average runtime, the average numbers of iterations, and the average primal and dual errors of the ADMM for the different choices of the parameter ρ are presented in Figure 1.Figure 1
Comparison with different parameterρ.It is clear from Figure1 (top) that when 2≤ρ≤4, the average runtime and the average iterations are reasonable. From Figure 1 (bottom), we can observe that, with the parameter ρ increasing, the primal error decreases while the dual error becomes bigger. Numerical simulations suggest that choosing ρ∈[2,4] could accelerate the convergence of the ADMM.


#### 5.1.2. The Effects of Tolerances for the Primal and Dual Residuals

Now, we take the different tolerances for the primal residualsϵpri and the dual residuals ϵdual to analyze the performance of the ADMM, where sparse r=0.1n, measurements p=4r, the penalty parameter ρ=2, and ϵpri=ϵdual=ϵ. We take two different signal lengths n=400 and n=600. We have repeated the same experiment 100 times by a set of randomly generated data. For different choices of ϵ, the average numbers of iterations and the executing time of the above ADMM algorithm (102) are presented in Table 1. It is shown that, with the increasing of precision, the number of iterations increases accordingly while increasing of executing time is not obvious.Table 1
Numerical results on the different tolerances for residuals.
n
ϵ
Iteration
Time
n
ϵ
Iteration
Time
400
10
-
1
21.56
0.0459
600
10
-
1
23.34
0.1182
400
10
-
2
37.80
0.0482
600
10
-
2
40.52
0.1474
400
10
-
3
55.68
0.0571
600
10
-
2
59.20
0.1720
400
10
-
4
78.04
0.0603
600
10
-
4
92.46
0.2238
400
10
-
5
119.80
0.0951
600
10
-
5
114.98
0.2590
400
10
-
6
139.00
0.0956
600
10
-
6
149.26
0.3063
400
10
-
7
161.64
0.1056
600
10
-
7
230.24
0.4346In Figure2(a), the full line describes the changes of the primal residuals rk. In Figure 2(b), the full line describes the changes of the dual residuals sk. The dotted lines in Figures 2(a) and 2 (b) represent the original residual tolerance ϵpri=10-6 and dual residuals tolerance ϵdual=10-6, respectively. From Figure 2, we can see that the two residuals descend monotonously.Figure 2
Norms of primal residual (a) and dual residual (b) versus iteration.
(a)
(b)


### 5.2. Reconstruction of Electroencephalogram Signal by Using Complex ADMM

Electroencephalogram (EEG) signal is a weak bioelectricity of brain cells group, which can be recorded by placing the available electrodes on the scalp or intracranial detect. The EEG signal could reflect the brain bioelectricity rhythmic activity regularity of random nonstationary signal. In this area, much is known for clinical diagnosis and brain function [42]. Because EEG data is large, a very meaningful work is to compress EEG data. It is to effectively reduce the amount of data at the same time and to guarantee that the main features basically remain unchanged [43, 44].In this paper, EEG signals are recorded with a g.USBamp and a g.EEGcap (Guger Technologies, Graz, Austria) with a sensitivity of 100 V, band pass filtered between 0.1 and 30 Hz, and sampled at 256 Hz. Data are recorded and analyzed using the ECUST BCI platform software package developed through East China University of Science and Technology [45, 46].We get the complex signalxo by performing the discrete Fourier transform (DFT) on EEG signal s; that is, xo=F(s), with signal s with length n=2000. Original EEG signal s is not sparse, but its DFT signal xo becomes the approximate sparse signal; see Figure 3. The hard threshold of xo is properly set, which leads to x with 90 percent zero valued entries, and x is the approximation of xo.Figure 3
Original EEG signal and its DFT signal.
(a)
Original EEG
(b)
Real part of DFT
(c)
Imaginary part of DFTTo get the compression of sparse signalx, we can first calculate b from b=Ax. A is a random complex matrix of size p×n, where p=800, n=2000, and the sampling rate is 40%. Substituting A and b into the optimization model as presented in (77), that is,(104)minimizex~1:Ax~=b,x~∈Cn, we can obtain the sparse optimal solution x~ by employing the ADMM algorithm (102) in Section 4, in which x~ is a good approximation of x. By applying the inverse discrete Fourier transform (IDFT) on x~, we can get the approximation s~ of original signal s; that is,(105)s~=F-1x~.The original signal s and its reconstruction signal s~ can be seen in Figure 4, in which (a) is the original signal s, (b) is the reconstruction signal s~, and (c) is the comparison of them. With the comparison of (a), (b), and (c) in Figure 4, we can observe that the reconstructed signal is in good agreement with the original signal and retains the leading characteristic. The relative error δ=norm(s-s~)/norm(s)=0.1858.Figure 4
Original EEG signal and its reconstructed signal in the complex number domain.
(a)
Original EEG signal
(b)
Reconstructed EEG signal in complex number domain
(c)
Comparison of original EEG signal and its reconstructed signalNow we separate the complex signalx into the real part xre and the imaginary part xim and then recast it into an equivalent real-valued optimization problem. We can calculate b1 from b1=Axre and b2 from b2=Axim, respectively. Here A is a random real matrix of size p×n, where p=800, n=2000, and the sampling rate is 40% which is similar to the one used in the complex number domain. Substitute A and b1,b2 into the optimization model as presented in (77); that is,(106)minimizex~re1:Ax~re=b1,x~re∈Rn,minimizex~im1:Ax~im=b2,x~im∈Rn. Although x~re and x~im can approach xre and xim, respectively, the reconstructed signal s~=F-1(x~re+ix~im) is not consistent with the original signal s; see Figure 5. The relative error δ=norm(s-s~)/norm(s)=0.5879. It can be seen that our new ADMM proposed in this paper performs better than the classic ADMM.Figure 5
Original EEG signal and its reconstructed signal in the real number domain.
(a)
Original EEG signal
(b)
Reconstructed EEG signal in the real number domain
(c)
Comparison of original EEG signal and its reconstructed signal


## 5.1. Numerical Simulation of the BP Algorithm with Random Data

Assume thatx∈Cn is a discrete complex signal interested. x itself is r-sparse, which contains (at most) r nonzero entries with r≪n. Select p(p<n) measurements uniformly at random matrix Ap×n via Ap×nx=b. Hence reconstructing signal x from measurement b is generally an ill-posed problem which is an undetermined system of linear equations. However, the sparsest solution can be obtained by solving the constrained optimization problem: (103)minimizex0:Ax=b,x∈Cn,where x0 is the l0-norm of x. Unfortunately, (103) is a combinatorial optimization problem of which the computational complexity grows exponentially with the signal size n. A key result in [17, 18] is that if x is sparse, the sparsest solution of (103) can be obtained with overwhelming probability by solving the convex optimization problem (77). Next, we consider CS which is actually a kind of application of the BP method in complex variables.


### 5.1.1. The Effects of the Parameterρ

We demonstrate the complex signal sampling and recovery techniques with a discrete-time complex signalx of length n=300 with sparsity r=30 which is generated randomly. Ap×n is a random sensing matrix with p=120. The variables u0 and z0 are initialized to be zero. We set the two tolerances of primal and dual residuals equal to 10-6. In order to understand the effects of the parameter ρ on the convergence, we set the penalty parameter ρ from 0.1 to 20 with the step 0.1. We have repeated the same experiment 100 times with the same parameter ρ. The average runtime, the average numbers of iterations, and the average primal and dual errors of the ADMM for the different choices of the parameter ρ are presented in Figure 1.Figure 1
Comparison with different parameterρ.It is clear from Figure1 (top) that when 2≤ρ≤4, the average runtime and the average iterations are reasonable. From Figure 1 (bottom), we can observe that, with the parameter ρ increasing, the primal error decreases while the dual error becomes bigger. Numerical simulations suggest that choosing ρ∈[2,4] could accelerate the convergence of the ADMM.


### 5.1.2. The Effects of Tolerances for the Primal and Dual Residuals

Now, we take the different tolerances for the primal residualsϵpri and the dual residuals ϵdual to analyze the performance of the ADMM, where sparse r=0.1n, measurements p=4r, the penalty parameter ρ=2, and ϵpri=ϵdual=ϵ. We take two different signal lengths n=400 and n=600. We have repeated the same experiment 100 times by a set of randomly generated data. For different choices of ϵ, the average numbers of iterations and the executing time of the above ADMM algorithm (102) are presented in Table 1. It is shown that, with the increasing of precision, the number of iterations increases accordingly while increasing of executing time is not obvious.Table 1
Numerical results on the different tolerances for residuals.
n
ϵ
Iteration
Time
n
ϵ
Iteration
Time
400
10
-
1
21.56
0.0459
600
10
-
1
23.34
0.1182
400
10
-
2
37.80
0.0482
600
10
-
2
40.52
0.1474
400
10
-
3
55.68
0.0571
600
10
-
2
59.20
0.1720
400
10
-
4
78.04
0.0603
600
10
-
4
92.46
0.2238
400
10
-
5
119.80
0.0951
600
10
-
5
114.98
0.2590
400
10
-
6
139.00
0.0956
600
10
-
6
149.26
0.3063
400
10
-
7
161.64
0.1056
600
10
-
7
230.24
0.4346In Figure2(a), the full line describes the changes of the primal residuals rk. In Figure 2(b), the full line describes the changes of the dual residuals sk. The dotted lines in Figures 2(a) and 2 (b) represent the original residual tolerance ϵpri=10-6 and dual residuals tolerance ϵdual=10-6, respectively. From Figure 2, we can see that the two residuals descend monotonously.Figure 2
Norms of primal residual (a) and dual residual (b) versus iteration.
(a)
(b)


## 5.1.1. The Effects of the Parameterρ

We demonstrate the complex signal sampling and recovery techniques with a discrete-time complex signalx of length n=300 with sparsity r=30 which is generated randomly. Ap×n is a random sensing matrix with p=120. The variables u0 and z0 are initialized to be zero. We set the two tolerances of primal and dual residuals equal to 10-6. In order to understand the effects of the parameter ρ on the convergence, we set the penalty parameter ρ from 0.1 to 20 with the step 0.1. We have repeated the same experiment 100 times with the same parameter ρ. The average runtime, the average numbers of iterations, and the average primal and dual errors of the ADMM for the different choices of the parameter ρ are presented in Figure 1.Figure 1
Comparison with different parameterρ.It is clear from Figure1 (top) that when 2≤ρ≤4, the average runtime and the average iterations are reasonable. From Figure 1 (bottom), we can observe that, with the parameter ρ increasing, the primal error decreases while the dual error becomes bigger. Numerical simulations suggest that choosing ρ∈[2,4] could accelerate the convergence of the ADMM.


## 5.1.2. The Effects of Tolerances for the Primal and Dual Residuals

Now, we take the different tolerances for the primal residualsϵpri and the dual residuals ϵdual to analyze the performance of the ADMM, where sparse r=0.1n, measurements p=4r, the penalty parameter ρ=2, and ϵpri=ϵdual=ϵ. We take two different signal lengths n=400 and n=600. We have repeated the same experiment 100 times by a set of randomly generated data. For different choices of ϵ, the average numbers of iterations and the executing time of the above ADMM algorithm (102) are presented in Table 1. It is shown that, with the increasing of precision, the number of iterations increases accordingly while increasing of executing time is not obvious.Table 1
Numerical results on the different tolerances for residuals.
n
ϵ
Iteration
Time
n
ϵ
Iteration
Time
400
10
-
1
21.56
0.0459
600
10
-
1
23.34
0.1182
400
10
-
2
37.80
0.0482
600
10
-
2
40.52
0.1474
400
10
-
3
55.68
0.0571
600
10
-
2
59.20
0.1720
400
10
-
4
78.04
0.0603
600
10
-
4
92.46
0.2238
400
10
-
5
119.80
0.0951
600
10
-
5
114.98
0.2590
400
10
-
6
139.00
0.0956
600
10
-
6
149.26
0.3063
400
10
-
7
161.64
0.1056
600
10
-
7
230.24
0.4346In Figure2(a), the full line describes the changes of the primal residuals rk. In Figure 2(b), the full line describes the changes of the dual residuals sk. The dotted lines in Figures 2(a) and 2 (b) represent the original residual tolerance ϵpri=10-6 and dual residuals tolerance ϵdual=10-6, respectively. From Figure 2, we can see that the two residuals descend monotonously.Figure 2
Norms of primal residual (a) and dual residual (b) versus iteration.
(a)
(b)


## 5.2. Reconstruction of Electroencephalogram Signal by Using Complex ADMM

Electroencephalogram (EEG) signal is a weak bioelectricity of brain cells group, which can be recorded by placing the available electrodes on the scalp or intracranial detect. The EEG signal could reflect the brain bioelectricity rhythmic activity regularity of random nonstationary signal. In this area, much is known for clinical diagnosis and brain function [42]. Because EEG data is large, a very meaningful work is to compress EEG data. It is to effectively reduce the amount of data at the same time and to guarantee that the main features basically remain unchanged [43, 44].In this paper, EEG signals are recorded with a g.USBamp and a g.EEGcap (Guger Technologies, Graz, Austria) with a sensitivity of 100 V, band pass filtered between 0.1 and 30 Hz, and sampled at 256 Hz. Data are recorded and analyzed using the ECUST BCI platform software package developed through East China University of Science and Technology [45, 46].We get the complex signalxo by performing the discrete Fourier transform (DFT) on EEG signal s; that is, xo=F(s), with signal s with length n=2000. Original EEG signal s is not sparse, but its DFT signal xo becomes the approximate sparse signal; see Figure 3. The hard threshold of xo is properly set, which leads to x with 90 percent zero valued entries, and x is the approximation of xo.Figure 3
Original EEG signal and its DFT signal.
(a)
Original EEG
(b)
Real part of DFT
(c)
Imaginary part of DFTTo get the compression of sparse signalx, we can first calculate b from b=Ax. A is a random complex matrix of size p×n, where p=800, n=2000, and the sampling rate is 40%. Substituting A and b into the optimization model as presented in (77), that is,(104)minimizex~1:Ax~=b,x~∈Cn, we can obtain the sparse optimal solution x~ by employing the ADMM algorithm (102) in Section 4, in which x~ is a good approximation of x. By applying the inverse discrete Fourier transform (IDFT) on x~, we can get the approximation s~ of original signal s; that is,(105)s~=F-1x~.The original signal s and its reconstruction signal s~ can be seen in Figure 4, in which (a) is the original signal s, (b) is the reconstruction signal s~, and (c) is the comparison of them. With the comparison of (a), (b), and (c) in Figure 4, we can observe that the reconstructed signal is in good agreement with the original signal and retains the leading characteristic. The relative error δ=norm(s-s~)/norm(s)=0.1858.Figure 4
Original EEG signal and its reconstructed signal in the complex number domain.
(a)
Original EEG signal
(b)
Reconstructed EEG signal in complex number domain
(c)
Comparison of original EEG signal and its reconstructed signalNow we separate the complex signalx into the real part xre and the imaginary part xim and then recast it into an equivalent real-valued optimization problem. We can calculate b1 from b1=Axre and b2 from b2=Axim, respectively. Here A is a random real matrix of size p×n, where p=800, n=2000, and the sampling rate is 40% which is similar to the one used in the complex number domain. Substitute A and b1,b2 into the optimization model as presented in (77); that is,(106)minimizex~re1:Ax~re=b1,x~re∈Rn,minimizex~im1:Ax~im=b2,x~im∈Rn. Although x~re and x~im can approach xre and xim, respectively, the reconstructed signal s~=F-1(x~re+ix~im) is not consistent with the original signal s; see Figure 5. The relative error δ=norm(s-s~)/norm(s)=0.5879. It can be seen that our new ADMM proposed in this paper performs better than the classic ADMM.Figure 5
Original EEG signal and its reconstructed signal in the real number domain.
(a)
Original EEG signal
(b)
Reconstructed EEG signal in the real number domain
(c)
Comparison of original EEG signal and its reconstructed signal


## 6. Conclusions

In this paper, the ADMM for separable convex optimization of real functions in complex variables has been studied. By using Wirtinger calculus, we have established the convergence of the algorithm, which is the generalization of the one obtained in real variables. Furthermore, the BP algorithm is given in the form of the ADMM, in which projection algorithm and the soft thresholding formula are generalized from the real number domain to the complex case. The simulation results demonstrate that the ADMM can quickly solve convex optimization problems in complex variables within the scopes of the signal compression and reconstruction, which is better than the results in the real number domain.



---

*Source: 104531-2015-12-31.xml*