Sparsity-Promoting Methods inRemote Control Systems
Masaaki Nagahara
Graduate School of Informatics, Kyoto University,
Yoshida Honmachi Sakyo-ku, Kyoto, 606-8501, JAPANnagahara@ieee.org

Abstract—In remote control systems, efficient representation ofcontrol signals is one of the crucial issues because of bandwidth
limitation of the communication channel, such as a wirelesscommunication link between the controller and the controlled
object. Recently, a new method based on the technique ofcompressive sampling has been proposed, in which control
signals are sparsely representation based on convex relaxation.
There exist however so many sparsity-promoting methods forcompressive sampling. In this study, we perform a comparative
study of sparsity-promoting methods and reveal their advantages
and disadvantages by simulation in view of remote control overrate-limited networks.
Index Terms—remote control, networked control, compressive
sampling, compressed sensing, sparse representation

I. I NTRODUCTION
In remote control [1], the controlled objects are locatedaway from controllers, and control signals are to be transmitted
through rate-limited channels such as wireless channels [2] orthe Internet [3]. To send control signals through such communication channels, efficient signal compression or representation is essential. For this, an approach has been proposed in[4], [5] for sparsifying control signals using the compressive
sampling technique [6], [7], [8] for remote control systems.
The aim of compressive sampling (aka compressed sensing) is to reconstruct signals from much less information by
assuming that the original signal is sparse [9], [10], [11]. The
core idea used in this area is to introduce a sparsity index in
the optimization. To be more specific, the sparsity index ofa vector v is defined by the amount of nonzero elements inv and is usually denoted by ∥v∥0 , called the “ℓ0 norm.” Thecompressed sensing problem is then formulated by ℓ0 -norm
optimization. By ℓ0 -norm optimization, one obtain a sparse
vector contains many 0-valued elements, and can be highly
compressed by only coding a few nonzeros and their locations.
A well-known example of this sparsity-inducing compressionis JPEG [12].
The associated optimization problem is however hard tosolve, since it is a combinatorial one [13]. One approach tothe combinatorial optimization is an iterative greedy algorithm
called Orthogonal Matching Pursuit (OMP) [14] and its extension, Compressive Sampling Matching Pursuit (CoSaMP)[15]. Another tractable approach for ℓ0 -optimization is convex
This research was supported in part under the JSPS Grant-in-Aid for
Scientific Research (C) No. 24560543.

relaxation to use the ℓ1 norm instead of the highly nonconvexℓ0 norm. This is called Basis Pursuit Denoising (BPDN) in
signal processing, and can be also effectively solved via fastalgorithms such as NESTA (Nestrov’s Algorithm) based on
Nesterov’s method. [19], and FISTA (Fast IST Algorithm) [16]based on iterative shrinkage-thresholding [17], [18].
In this study, we perform a comparative study of thealgorithms, OMP, CoSaMP, NESTA, and FISTA, for remotecontrol systems. We reveal their advantages and disadvantages
by simulation in view of remote control over rate-limited
networks.
Notation
For vector v = [v1 , . . . , vn ]⊤ ∈ Rn , the ℓ1 and ℓ2 norms
are respectively defined by
( n
)1/2
n
∑
∑
2
∥v∥1 :=|vi |, ∥v∥2 :=
vi
.
i=1

i=1

The set of indices of nonzero elements in v is denoted bysupp(v) := {j : vj ̸= 0}, and the “ℓ0 ” norm is definedby ∥v∥0 := |supp(v)|, where |supp(v)| is the number ofmembers in the set supp(v). For a matrix Φ, its norm isdefined by∥Φv∥2∥Φ∥ := max= σmax (Φ),∥v∥2 ̸=0 ∥v∥2
where σmax (Φ) is the maximum singular value of Φ. We
denote by L2 [0, T ] the set of all square integrable functions
on [0, T ] (T > 0) , endowed with the inner product∫ T⟨f, g⟩ :=f (t)g(t)dt, f, g ∈ L2 [0, T ]0and the L2 norm ∥f ∥2 := ⟨f, f ⟩1/2 .
II. C ONTROL P ROBLEM
A. Plant model
In this study, we consider a control problem of a linear
system P on a finite time interval (or horizon) [0, T ], modeledby the following state-space equations:{ẋ(t) = Ax(t) + bu(t),P :(1)y(t) = c⊤ x(t), x(0) = x0 ∈ Rν , t ∈ [0, T ],where A ∈ Rν×ν and b, c ∈ Rν×1 . The initial state x0 ∈ Rνis assumed to be given. We also assume that the system isstable.

r(t)

θ

Ψ

K

u(t)

y(t)
P

By this, the sampled tracking error
Ed (u) = hN∑|y(tn ) − r(tn )|2

n=1

x0

is rewritten as
Fig. 1.

2

Remote control system

Ed (θ) = h ∥Gθ − Hx0 − r∥2 ,where

B. Tracking problem
For the system (1), we here consider the tracking problem;
given a reference signal r(t) over the time interval [0, T ] andthe state-space observation x(0) at the initial time t = 0, finda control signal u(t) over [0, T ] that minimizes the tracking
error
E(u) := ∥y − r∥22 = ∥P u − r∥22 .
In addition, we assume remote control where the controller,
denoted by K, that generates the control signal is located awayfrom the controlled plant P . Fig. 1 shows a block diagram of
the remote control system where the dotted lines indicate ratelimited networks such as wireless networks. In this system, thecontrol signal u(t) over [0, T ] should be represented in a finite
dimensional vector, denoted by θ, which is transmitted throughthe rate-limited network. The received vector is transformed
to the control signal u(t), t ∈ [0, T ] by the decoder Ψ definedbyM∑
u = Ψθ =θ m ψm ,m=−M

where ψm is the m-th Fourier basis function on [0, T ], that is,
1
2πm
ψm (t) := √ exp(jωm t), t ∈ [0, T ], ωm :=.
T
T
Note that the control signal u is band-limited up to thefrequency ωM = 2πM/T . Let VM denote the set of suchsignals, that is,
VM := span{ψm : m = −M, . . . , M } ⊂ L2 [0, T ].
Then we have u ∈ VM . We also assume that the reference
signal r is in VM .
Under the above assumptions, the output y(t) can be rewritten in terms of the Fourier coefficients θi as
y(t) = c⊤ exp(tA)x0 +M∑θm ⟨κ(t, ·), ψm ⟩,

m=−M

where

{

κ(t, τ ) :=c⊤ exp [(t − τ )A] b,0,

if 0 ≤ τ < t ≤ T,otherwise.

Sampling the continuous-time signal y(t) with sampling time
h := T /(N − 1), where N := 2M + 1, the dimension of thesignal subspace VM , we havey(tn ) = c⊤ exp(tn A)x0 +M∑θm ⟨ϕn , ψm ⟩,m=−M

tn := (n − 1)h, ϕn := κ(tn , ·), n = 1, 2, . . . , N.



⟨ϕ1 , ψ−M ⟩ ⟨ϕ2 , ψ−M ⟩
G := 
..

.

...
...
..
.


⟨ϕ1 , ψM ⟩⟨ϕ2 , ψM ⟩  ∈ CN ×N ,
..

.

(2)

⟨ϕN , ψ−M ⟩ . . . ⟨ϕN , ψM ⟩
 ⊤c exp (t1 A)
r(t1 )
 c⊤ exp (t2 A) 
 r(t2 ) r :=  .  ∈ RN , H :=  ∈ RN ×ν .
..
 .. .
c⊤ exp (tN A)

r(tN )

C. Sparse control vector
Since the network is rate-limited, the vector θ should becompressed (or represented) in as small data size as possible.
For this purpose, we have proposed in [4] to use the methodof compressive sampling.
Let U be a random “decimation” matrix of the form


ei(1)
 ei(2) 
U =  .  ∈ {0, 1}K×N ,
.
 . 
ei(K)
where i(1) < i(2) < · · · < i(K) are the random variables of
the uniform distribution on {1, 2, . . . , N }, andi∨ei := [0, . . . , 0, 1, 0, . . . , 0],

i = 1, 2, . . . , N.

This is a model of low rate random sampling of a signal on[0, T ]. We define the random sampling instants byT, k = 1, 2, . . . , K < N.
N −1
By using the matrix U , random sampling of y(t) on [0, T ] isgiven by:) (y (ti(1) ) y ti(2) 
y= = U Gθ + U Hx0 .
..


.
(
)
y ti(K)ti(k) = i(k) · h = i(k) ·The associated optimization problem is formulated bymin ∥θ∥0 subject to ∥Φθ − α∥22 ≤ ϵ,

θ∈RN

(3)

where Φ := U G and α = U (r − Hx0 ). The objective
here is to minimize ∥θ∥0 to obtain a sparse control vector
with an ℓ2 norm constraint for tracking performance. The
optimization (3) is a combinatorial one, which is hard tosolve in a reasonable computational time. For this, we adopt
computationally tractable algorithms in the next section.

Algorithm 1 OMP for sparse control vector θ
Require: α ∈ RK {observed vector}
Ensure: θ {sparse control vector}
θ[0] := 0, r[0] := α − Φθ[0], S[0] := supp{θ[0]} = ∅.
k := 0.
while ∥r[k]∥2 > ϵ do
for j = 1 to N doϕ⊤j r[k]= arg min ∥ϕj z − r[k]∥22 .
zj :=∥ϕj ∥22z∈R
ej := ∥ϕj zj − r[k]∥22 .
end forFind a minimizer j0 ̸∈ S[k] such that ej0 ≤ ej , for allj ̸∈ S[k].
S[k + 1] := S[k] ∪ {j0 }θ[k + 1] :=arg min∥Φθ − α∥22 .
supp{θ}=S[k+1]r[k + 1] := α − Φθ[k + 1].
k := k + 1.
end while
return θ = θ[k].

III. S PARSITY-P ROMOTING M ETHODS
There are many algorithms to obtain a numerical solutionof (3): greedy pursuit, bayesian framework, nonconvex optimization, and brute force [20]. In this article, we adopt greedy
pursuit and convex optimization, in particular.
A. Greedy pursuit
Among many methods of greedy pursuit, the orthogonalmatching pursuit (OMP) [14] is the most popular. Algorithm
1 shows the algorithm of OMP. If the sparsity of the solution, ∥θ∥0 is previously known, then one can adopt a fastermethod called CoSaMP (Compressive Sampling Matching
Pursuit) [15]. The CoSaMP algorithm takes advantage ofRIP (Restricted Isometry Property) of sparse signals. For thealgorithm, see [15].
B. Convex optimizationmin ∥θ∥1 subject to ∥Φθ − α∥22 ≤ ϵ.

(I + δλk Φ⊤ Φ)η[k] = δλk Φ⊤ α + q[k].
Solve the following linear equation for ζ[k]:(I + δλk Φ⊤ Φ)ζ[k] = δλk Φ⊤ α + θ[0] + q[k].
2θ[k + 1] := k+4
ζ[k] + k+2
k+4 η[k].
w[k + 1] := satδ (θ[k + 1]).
1v[k + 1] := 2k+4w[k + 1] + v[k].
k := k + 1.
until |fδ (θ[k−1])−fδ (θ[k−2])| ≤ EPS or k ≥ MAXITER.
return θ = θ[k − 1].

By using Huber’s function, the optimization (4) is approximated as
min fδ (θ) subject to ∥Φθ − α∥2 ≤ ϵ.

θ∈RN

(5)

Note that for any θ ∈ RN , fδ (θ) → ∥θ∥1 as δ → 0. Applying
the Nesterov’s method to the optimization problem (5), we
obtain the NESTA algorithm shown in Algorithm 2. In thisalgorithm, satδ (θ) is a saturation function defined by{1
θi ,if |θi | < µ,[satδ (θ)]i := µsgn(θi ), otherwise,
where [·]i is the i-th element of a vector, and{1, if x ≥ 0,
sgn(x) :=−1, if x < 0.
Note that the function satδ (θ) is the gradient of fδ (θ), i.e.,

Another method to solve the optimization (3) is convex
relaxation. For the relaxation, we replace the term ∥θ∥0 by theℓ1 norm ∥θ∥1 , and solve the following convex optimization:
θ∈RN

Algorithm 2 NESTA for sparse control vector θ
Require: α ∈ RK {observed vector}
Ensure: θ {sparse control vector}
θ[0] := 0, w[0] := satδ (θ[0]), v[0] := 21 w[0].
k := 1.
repeatq[k] := θ[k]{− δw[k].
}1λk := max 0, δϵ∥α − Φq[k]∥2 − 1δ .
Solve the following linear equation for η[k]:(4)The solution of this optimization can be obtained by a standard
method of convex programming, such as interior point method.
However, there exist much faster method called NESTA(Nesterov’s Algorithm) [19]. The main idea of the NESTA
algorithm is to replace the non-smooth ∥θ∥1 term by a smoothfunction called Huber’s function defined as
{
N
1 2
∑
θ ,
if |θi | < µ,fδ (θ) =|θi |δ , |θi |δ := 2µ i µ|θi | − 2 , otherwise.
i=1

satδ (θ) = ∇fδ (θ).
It is known that the sequence {θ[0], θ[1], . . . } converges to thesolution of (5) for any initial vector θ[0] at the convergence
rate O(1/k 2 ).
We can also adopt another method for the optimization (4)by considering the following Lagrange form:
min F (θ) = min ∥Φθ − α∥22 + µ∥θ∥1 .

θ∈RN

θ

(6)

For this optimization problem, FISTA (Fast Iterative
Shrinkage-Thresholding Algorithm) is known to be a fastalgorithm [16]. Algorithm 3 shows the FISTA algorithm. In
this algorithm, shrinkµ (θ) is a shrinkage function defined by[shrinkµ (θ)]i := sgn(ηi )(|ηi | − λ/c)+ ,Algorithm 3 FISTA for sparse control vector θ
Require: α ∈ RK {observed vector}
Ensure: θ {sparse control vector}
θ[0] := 0, θ̃[1] := 0, β[1] := 0.
k := 1.
repeat(1 ⊤)θ[k] := shrinkµ √c Φ (α − Φθ̃[k]) + θ̃[k] .

1
2
4 + β[k] .
()θ̃[k + 1] := θ[k] + β[k]−1
β[k+1] θ[k] − θ[k − 1] .

β[k + 1] := 12 +k := k + 1.
until |F (θ[k −1])−F (θ[k −2])| ≤ EPS or k ≥ MAXITER.
return θ = θ[k − 1].

where (x)+ := max{x, 0} for x ∈ R. If the parameter c
is chosen such that c ≥ ∥Φ∥, the sequence {θ[0], θ[1], . . . }converges the optimal solution of (6) for any initial vector
θ[0]. The convergence rate is known to be O(1/k 2 ).
IV. S IMULATION R ESULTS
We here perform a comparative study of existing methods
for the optimization (3) based on simulation. The state-space
matrices of the controlled plant P are given by
[
]
[ ]
[
]
0
1
0
−0.5
A=
, b=
, c=
.
−0.5 −1.5
11
The transfer function of this system is given by
P (s) =

s − 0.5.
(s + 0.5)(s + 1)

The period T is 2π. The number of basis {ψm } is N = 2M +1 = 51 (M = 25). The reference signal r(t) is given byr(t) = sin(5t) + cos(12.5t).
For compressive sampling, we take K = 51/3 = 17 random
samples among N = 51 sampled data, that is the compression
ratio is 1/3.
With these parameters, we run 100 simulations with randomly generated initial state x0 and matrix U . We adopt
four existing methods: OMP, CoSaMP, NESTA, and FISTA
introduced in the previous section. Table I shows the average values of the control performance (or tracking error)∥Φθ − α∥2 , sparsity ∥θ∥0 , computational (CPU) time (sec),
and the number of iterations. The result indicates that OMP
shows the best control performance and needs the fewest
iterations, while FISTA gives the sparsest control vector andneeds the shortest computational time. Note that NESTA gives
dense vectors since the algorithm uses Huber’s function fδ ,which leads to vectors many elements of which are almost zero(not exactly zero). To obtain a sparse vector with NESTA, one
needs truncation. FISTA is quite fast and leads to very sparse
vectors, but one should choose an appropriate parameter µ > 0
in the algorithm. This requires some trial-and-error process. As
a result, OMP is the best choice for the remote control systemin the simulation.

TABLE IC OMPARISON
methodOMP
CoSaMP
NESTA
FISTA∥Φθ − α∥20.330681.454
0.62325
1.1084

∥θ∥0
6.66
5
51
3.79

CPU time (sec)0.026208
0.0562191.183
0.0016118# of iteration
9.81
345.581000
17.59

V. C ONCLUSION
In this article, we have shown a comparative study ofsparsity-promoting methods in remote control systems bysimulation. The simulation result shows OMP (orthogonalmatching pursuit) is the best for the remote control system
considered in the simulation. Since many other algorithms
have been also proposed for sparsity, it is important to seek
for more suitable algorithm for remote control systems.
R EFERENCES[1] C. Sayers, Remote Control Robotics. Springer, 1999.
[2] A. F. T. Winfield and O. E. Holland, “The application of wireless local
area network technology to the control of mobile robots,” Microprocessors and Microsystems, vol. 23, pp. 597–607, Mar. 2000.
[3] R. C. Luo and T. M. Chen, “Development of a multibehavior-based
mobile robot for remote supervisory control through the internet,”
IEEE/ASME Trans. Mechatron., vol. 5, no. 4, Dec. 2000.
[4] M. Nagahara, T. Matsuda, and K. Hayashi, “Compressive sampling for
remote control systems,” IEICE Trans. on Fundamentals, vol. E95-A,
no. 4, pp. 713–722, Apr. 2012.
[5] M. Nagahara, D. E. Quevedo, T. Matsuda, and K. Hayashi, “Compressive sampling for networked feedback control,” Proc. IEEE ICASSP, pp.
2733–2736, Mar. 2012.
[6] E. J. Candes, “Compressive sampling,” Proc. International Congress ofMathematicians, vol. 3, pp. 1433–1452, Aug. 2006.
[7] E. J. Candes and M. B. Wakin, “An introduction to compressive
sampling,” IEEE Signal Processing Magazine, vol. 25, no. 2, pp. 21–30,
Mar. 2008.
[8] D. L. Donoho, “Compressed sensing,” IEEE Trans. Information Theory,
vol. 52, no. 4, pp. 1289–1306, Apr. 2006.
[9] M. Elad, Sparse and Redundant Representations. Springer, 2010.
[10] J.-L. Starck, F. Murtagh, and J. M. Fadili, Sparse Image and Signal
Processing. Cambridge University Press, 2010.
[11] Y. C. Eldar and G. Kutyniok, Compressed Sensing: Theory and Applications. Cambridge University Press, 2012.
[12] G. K. Wallace, “The JPEG still picture compression standard,” Commun.
ACM, vol. 34, no. 4, pp. 30–44, Apr. 1991.
[13] B. K. Natarajan, “Sparse approximate solutions to linear systems,” SIAM
J. Comput., vol. 24, no. 2, pp. 227–234, 1995.
[14] S. G. Mallat and Z. Zhang, “Matching pursuits with time-frequency
dictionaries,” IEEE Trans. Signal Processing, vol. 41, no. 12, pp. 3397–
3415, Nov. 1993.
[15] D. Needell and J. Tropp, “CoSaMP: Iterative signal recovery fromincomplete and inaccurate samples,” Applied and Computational Harmonic Analysis, vol. 26, no. 3, pp. 301–321, 2009.
[16] A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding
algorithm for linear inverse problems,” SIAM J. Imaging Sciences, vol. 2,
no. 1, pp. 183–202, Jan. 2009.
[17] M. A. Figueiredo and R. Nowak, “An EM algorithm for wavelet-based
image restoration,” IEEE Trans. Image Processing, vol. 12, no. 8, pp.
906–916, Aug. 2003.
[18] M. Zibulevsky and M. Elad, “L1-L2 optimization in signal and image
processing,” IEEE Signal Processing Magazine, vol. 27, pp. 76–88, May
2010.
[19] S. Becker, J. Bobin, and E. J. Candes, “NESTA: A fast and accurate
first-order method for sparse recovery,” SIAM J. Imaging Sci., vol. 4,
no. 1, pp. 1–39, 2011.
[20] J. A. Tropp and S. J. Wright, “Computational methods for sparse solutionof linear inverse problems,” Proc. IEEE, vol. 98, no. 6, pp. 948–958,
June 2010.