Youla Coding and Computation of Gaussian Feedback
Capacity
Chong Li and Nicola Elia

arXiv:1801.04190v1 [cs.IT] 12 Jan 2018

Abstract
In this paper, we propose an approach to numerically compute the feedback capacity of stationary finite dimensional Gaussian channels and construct (arbitrarily close to) capacity-achieving feedback codes. In particular,
we first extend the interpretation of feedback communication over stationary finite dimensional Gaussian channels
as feedback control systems by showing that, the problem of finding stabilizing feedback controllers with maximal
reliable transmission rate over Youla parameters coincides with the problem of finding strictly causal filters to achieve
feedback capacity derived in [2]. This extended interpretation provides an approach to construct deterministic feedback
coding schemes with double exponential decaying error probability. We next propose asymptotic capacity-achieving
upper bounds, which can be numerically evaluated by solving finite dimensional convex optimizations. From the filters
that achieve the upper bounds, we apply the Youla-based interpretation to construct feasible filters, i.e., feedback codes,
leading to a sequence of lower bounds. We prove the sequence of lower bounds is asymptotically capacity-achieving.

I. I NTRODUCTION
To facilitate the reading, we first introduce some notations as follows.
Notation: Uppercase and corresponding lowercase letters (e.g.Y, Z, y, z) denote random variables and realizations,
∞ refers to the set of bounded continuous functions
respectively. log denotes the logarithm base 2 and 0 log 0 = 0. C[a,b]

on [a, b]. We use x0 to denote the transpose of a real vector or matrix x.
RH 2 = RH ∞ denotes the set of real-rational transfer functions corresponding to the z−transform of the impulse
response of linear time invariant finite-dimensional (LTI-FD) causal stable systems [3] 1 . We interchangeably use
RH ∞ and RH 2 throughout the paper. A function f (z) ∈ RH 2 has all its poles strictly inside the unit disc.
In this paper, we consider a discrete-time Gaussian channel with noiseless feedback.
Assumption 1: The additive Gaussian channel is modeled as
Yi = Ui +Wi ,

i = 1, 2, · · ·

(1)

Dr. Chong Li is with Qualcomm Research, Bridgewater, NJ, 08807. chongl@qti.qualcomm.com
Dr. Nicola Elia is with the Department of Electrical and Computer Engineering, Iowa State University, Ames, IA, 50011 nelia@iastate.edu
This paper was presented (in part) in [1].
1 In

general, H∞ ⊂ H2 , for discrete-time systems. When restricted to the real-rational functions, however, RH ∞ = RH 2 . Note that [3] uses

λ -transform instead of the z−tranform, where λ = 1/z.

where the Gaussian noise sequence {Wi }∞
i=1 is assumed to be wide-sense stationary and has power spectral density
Sw (e jθ ) > 0 for ∀θ ∈ [−π, π). Moreover, the power spectral density satisfies the Paley-Wiener condition,
1
2π

Z π
−π

| log Sw (e jθ )|dθ < ∞.

Unless the contrary is explicitly stated, we use the term “stationary” to refer to stationary in wide sense.
For a code with achievable rate R, we specify a sequence of (n, 2nR ) channel codes as follows. M is a uniformly
distributed message index where M = m ∈ {1, 2, 3, · · · , 2nR }. There exists an encoding process Ui (M,Y i−1 ), where
Y i−1 = {Y0 ,Y1 , · · · ,Yi−1 }, for i = 1, 2, · · · , n and U1 (M,Y 0 ) = U1 (M) with average transmit power constraint. That is,
the channel input Ui is determined by the message index M and the previous channel output Y i−1 . Furthermore, there
(n)

exists a decoding function g: Y n → {1, 2, · · · , 2nR } with an error probability satisfying Pe

=

1
2nR

nR

∑2m=1 P(g(Y n ) 6=

m|M = m) ≤ εn , where limn→∞ εn = 0. The objective of communication is to deliver M to the receiver at highest
code rate with arbitrarily small error probability. The feedback capacity C f b is defined as the supremium of all
achievable rates R.
As shown in [2], the feedback capacity of channel (1) with the average power budget P can be characterized by


Z
Sv (e jθ )
1 π 1
jθ 2
+
|1
+
Q(e
)|
dθ ,
log
C f b = max
SV ,Q 2π −π 2
Sw (e jθ )
Z
1 π
s.t.
Sv (e jθ ) + |Q(e jθ )|2 Sw (e jθ )dθ ≤ P,
2π −π
(2)
jθ
Sv (e ) ≥ 0,
Q(e jθ ) =

∞

∑ qk e jkθ ∈ H2 is strictly causal.
k=1

Assumption 2: In this paper, we further assume that noise W is the output of a LTI-FD stable system H ∈ RH 2 ,
not necessarily minimum phase, driven by white Gaussian noise with zero mean and unit variance. The power
spectral density (PSD) of W has a canonical spectral factorization given by Sw (e jθ ) = |H(e jθ )|2 .
Note that any stationary process can be approximated with arbitrary accuracy by this LTI-FD filtering model and
this approximation is very “efficient” as it corresponds to the rational approximation of the spectral density [4].
Given the capacity of white Gaussian channel with feedback has been well known, in this paper we consider
non-white channel noise spectra.
Assumption 3: Noise W has a colored (non-white) power spectral density.
Remark 1: Under Assumption 2 and 3, [2] showed that the optimal solution to (2) must have Sv = 0. Furthermore,
Proposition 5.1 in [2] presented that the capacity is necessarily achieved by a rational filter Q. Therefore, without
loss of optimality, we restrict the search space Q to RH 2 in (2).
Specifically, (2) can be simplified to
C f b = max
Q

s.t.

1
2π

1
2π

Z π
−π

Z π
−π

log |1 + Q(e jθ )|dθ ,

|Q(e jθ )|2 Sw (e jθ )dθ ≤ P,

Q(e jθ ) =

∞

∑ qk e jkθ ∈ RH 2 is strictly causal.
k=1

(3)

Remark 2: Under Assumption 2 and 3, the optimal Q has no zeros on the unit circle (Proposition 5.1 (ii) in [2]).
While the above characterizations are elegant, they are infinite dimensional. As stated in [2],

“· · · except for the first-order ARMA spectrum, it is still a nontrivial task to find analytically (or even numerically)
the optimal filter and corresponding feedback capacity”.

In this paper, we aim to solve the above optimization problem and explicitly construct capacity-achieving feedback
codes. Firstly, we revisit and extend the interpretation of feedback communication over Gaussian channels as
feedback control problems [5]. In particular, we highlight the central role of Youla parameterization over all
stabilizing controllers in connecting these two theories by showing that the characterization of the maximum-rate
over all stabilizing controllers and the feedback capacity over all coding schemes coincide. This new interpretation
provides an approach to construct an explicit (sub-)optimal communication scheme (i.e. encoder and decoder)
directly from the filter Q in (3). Next, to find an optimal Q, we provide an alternative characterization of the
feedback capacity, from which an asymptotic capacity-achieving sequence of upper bounds is derived and can
be numerically evaluated by solving finite dimensional convex optimizations. Furthermore, from the filters Q that
achieve upper bounds, we derive a sequence of lower bounds on the feedback capacity by constructing deterministic
feedback codes with double exponential decaying error probability. Essentially, this constructed coding scheme has
the structure of a generalized Schalkwijk-Kailath scheme which have been extensively studied by Elia [6], Kim [2],
Liu-Elia [7], Shayevitz-Feder [8], Ardestanizadeh-Minero-Franceschetti [9] [10] and others. It is proved that the
sequence of lower bounds converges to the capacity, leading to an asymptotically optimal feedback coding scheme.
It is worth noting that [2] (see Theorem 6.1 and Lemma 6.1) has shown a structural, not computable result that
a k-dimensional generalization of the Schalkwijk-Kailath coding scheme achieves the feedback capacity for any
auto-regression moving average (ARMA) noise spectrum of order k, while we herein provide a numerically explicit
feedback coding scheme from Q by leveraging control-oriented derivations.
A. Related Work
We review the literature along two avenues: information theory and feedback control theory. As a complete
survey is vast and most of them are out of the scope of our discussion, we herein list most relevant results to this
paper. In the field of information theory, the investigation on feedback Gaussian capacity has been experiencing
a decade journey. [11] and its sequel [12] are recognized as the first works on feedback Gaussian channels by
proposing feedback coding schemes. [13] [14] developed an elegant linear feedback coding scheme of achieving
the capacity of additive white Gaussian noise (AWGN) channel with noiseless feedback. Thereafter, several works
by Butman [15], [16], Tiernan [17] [18], Wolfowitz [19] and Ozarow [20] [21] extended this notable result to ARMA
Gaussian channels, with an objective to find channel capacity and optimal feedback codes. As a consequence, many
interesting upper and lower bounds were obtained. Based on the insight/results from aforementioned literature,
[22] made a major breakthrough on characterizing the n-block capacity of arbitrary feedback Gaussian channels
by using asymptotic equipartition property (AEP) theorem. It was also shown that feedback capacity for arbitrary

Gaussian channels cannot be increased by factor two or half bit. This n-block capacity was extended to the case of
feedback Gaussian channels with noisy feedback where capacity bounds and other interesting results were obtained
[23]–[27]. As hinted by this n-block capacity characterization, [2] developed a variational characterization on the
capacity of stationary feedback Gaussian channels, which is an infinite dimensional optimization problem. For the
first-order ARMA noise, this variational characterization yields a closed-form solution on the capacity and shows
the optimality of the Schalkwijk-Kailath scheme.
In the field of feedback control theory, many control-based technical tools have been utilized to attack the problem
of finding feedback Gaussian channel capacity and capacity-achieving codes. [6] proposed the derivation of feedback
communication schemes based on a feedback control method. These results were obtained from considering the
problem of stabilization of a given unstable plant over a Gaussian communication channel. The communication rate
(in the sense of Shannon) over the channel was connected to the degree on instability of the plant. The minimal
transmission power for a given unstable plant was obtained by solving the classical H2 (or Linear Quadratic
Gaussian) problem. However, plants with the same degree of instability may require different transmission power
to be stabilized. [6] provided the plants that can be stabilized most efficiently, i.e. with the least transmission
power for a given degree of instability for special case channels. This approach provides a method of finding
feedback coding scheme for Gaussian channels. The approach has been further extended to Markov channels with
channel state information [28], connected to the classical linear quadratic Gaussian (LQG) control problem [9]
[10]. Specifically, [10] presented a code for an k-receiver additive white Gaussian noise broadcast channel with
feedback and characterized its sum-rate performance by using the tools from LQG control theory. Finally, [7]
extended the convergence of the fundamental limitations of control and communication to include the limitations of
estimation. In light of this unified framework, a set of achievable rates of feedback Gaussian channels were obtained
by constructing specific feedback coding schemes via control-oriented approaches. [29] converted the problem of
finding feedback Gaussian channel capacity into a form of stochastic control and used dynamic programming to
compute the n-block capacity.
II. P RELIMINARIES
In this section we formalize some preliminary results, which directly follow from [2] and are useful for our
derivations. In addition, we briefly review the theory of Youla parameterization, a useful tool to construct capacityachieving feedback codes.
A. Capacity Characterization Revisited
From (2), it is immediate that
1
C f b = sup
2π
SV ,Q
s.t.

1
2π

Z π
−π


Sv (e jθ )
jθ 2
log
+ |1 + Q(e )| dθ ,
2
Sw (e jθ )

Z π
1
−π



Sv (e jθ ) + |Q(e jθ )|2 Sw (e jθ )dθ ≤ P,

jθ

Sv (e ) ≥ 0,
Q ∈ RH 2 is strictly causal,

(4)

Fig. 1. Feedback Gaussian channels with two degree of freedom controllers and exogenous input.

since RH 2 is a dense subspace of H2 .
Proposition 1: Under Assumption 2 the “sup” in (4) is achieved.
Proof: If the noise spectrum is non-white, then the result follows from (3) and Sv must be equal to 0. If the
noise spectrum is white, it is well known that feedback does not change the capacity. Thus, we can let Sv be given
by the classical water-filling solution and Q = 02
One will see in Section III that the above capacity characterization with Q ∈ RH 2 (equivalently, RH ∞ ) is
useful to connect feedback communications and feedback control, a crucial step to derive the capacity-achieving
feedback codes.
B. Youla Parameterization
We consider the closed loop set-up shown in Figure 1 where F is a single-input-single-output (SISO) LTI-FD
plant represented by a rational transfer function, the additive disturbance w is generated according to Assumption 2,
K̃ is a two-degree of freedom (2dof) stabilizing controller, and v is an exogenous input with power spectral density
Sv . In what follows, we review a fundamental result of linear control theory: Youla-Kuc̆era parameterization (also
simply known as Youla parameterization or Q-parameterization) [30]. In particular, we consider the set of two
degree of freedom stabilizing controllers for a given plant, and the resulting achievable closed loop maps.
Lemma 1: (Youla Parameterization of Controllers, Chapter 5.2 in [30]) Given a SISO plant F, there exist
N, D, X, U ∈ RH ∞ , such that
F=

N
,
D

with
NU + DX = 1.
2 Note

that there may exist infinite number of solutions with Q 6= 0 ∈ RH 2 in this case.

(5)

The above factorization
(5) of plant F is called coprime factorization. Moreover, the set of all 2dof rational
 
v
controllers K̃ :   → u, that stabilize the closed loop of plant F is given by (Chapter 5.4, Theorem 2 in [30])3
y
K̃ =

1
[Qv , (U + DQ)] , for Qv ∈ RH ∞ , Q ∈ RH ∞ .
X + NQ

Note that the order of the controller is not restricted in any way. The main benefit of Youla paramterization will
be seen later in the next section. All omitted technical proofs in the following sections are provided in Appendix.
III. F EEDBACK C ONTROL I NTERPRETATION OF F EEDBACK C APACITY FOR G AUSSIAN C HANNELS
In this section, we propose a control-theoretic approach to derive the feedback capacity formula for the finite
dimensional LTI Gaussian channels. The proposed approach based on [6] reveals the essential role of Youlaparameter in connecting the theories of feedback communications and feedback control, and provides feasible
feedback communication schemes with guaranteed transmission rate. As will be proved in the paper, this constructed
coding scheme achieves (arbitrarily close to) the feedback capacity.
Consider a channel
Yi = Ui +Wi ,

(6)

where Wi satisfies Assumption 2. As shown in Fig. 1, we are interested in the closed loop stabilization problem
over the given channel. Following [6], we consider the map from Ui to Yi as the plant, F, where F = 1 is stable.
Now, according to Youla parametrization of stabilizing controllers, let D = X = 1, N = 1 and U = 0. Then all
the 2dof LTI-FD stabilizing controllers for the plant F = 1 have the following expression, represented as transfer
functions:
K̃ = [Kv , K] = (1 + Q)−1 [Qv , Q]

(7)

where Q ∈ RH ∞ , Qv ∈ RH ∞ .
Working with Q (and Qv ) instead of K̃ is more convenient. The main advantage comes from the fact that the
above transformation convexifies the set of achievable closed loop maps by a stabilizing controller. In particular,
Y

= W + KvV + KY
= (1 − K)−1 (W + KvV )

(8)

(a)

= (1 + Q)W + QvV,

where step (a) follows from (7). Similarly,
U = Y −W = QW + QvV.

(9)

Remark 3: In what follows, we restrict our attention to strictly causal feedback operations. This is equivalent to
restricting Q to be strictly proper. From (9), we can see that this way u is not instantaneous function of w.
3 Note

that the controller parameterizaiton in [30] is given by X − NQ with a negative sign on NQ. This should not be confusing as the

closed-loop systems considered in [30] has negative feedback while it is positive in our model.

In information theory of feedback communication systems, it is known that directed information

4

from channel

inputs to outputs measures the causal information transfer for channels with noiseless feedback and can be used
to characterize the feedback channel capacity. In our model, the directed information from channel inputs U T to
channel outputs Y T is defined by
T

I(U T → Y T ) = ∑ I(U i ;Yi |Y i−1 ),
i=1

where

I(U i ;Yi |Y i−1 )

denotes the conditional mutual information. We next characterize the average directed infor-

mation and the average power of channel input U in terms of the Youla parameters.
Lemma 2: Given the Youla parametrization in (7) with Q strictly proper,


Z
|Qv (e jθ )|2 Sv (e jθ )
1 π
1
T
T
jθ 2
log
+ |1 + Q(e )| dθ ,
lim I(U → Y ) =
T →∞ T
4π −π
Sw (e jθ )
1
T →∞ T
lim

T

∑ Ut2 =

t=1

1
2π

Z π
−π

|Q(e jθ )|2 Sw (e jθ ) + |Qv (e jθ )|2 Sv (e jθ )dθ .

Note that the derivation of the characterization of the above averaged directed information mainly follows from
the proof of Theorem 1 in [6]. The detailed proof of this lemma can be found in Appendix. With this lemma in
hand, we are now interested in finding the largest directed information rate for the closed-loop stabilization under
constrained average power, by searching over Q ∈ RH ∞ strictly proper, Qv ∈ RH ∞ and PSD Sv . We see, however,
that Qv can simply be chosen to be equal to 1. Putting above together, we can directly obtain the following theorem.
Theorem 1: Consider an additive Gaussian channel in (6) under Assumption 2. Given the average channel input
power budget P, the largest directed information rate (in the sense of Shannon) over all strictly causal LTI stabilizing
controllers can be characterized by (4).
This theorem indicates that the feedback capacity characterization (4) derived from information theory can be
equivalently derived from a control approach based on Youla parameterization. According to Remark 1 and
Assumption 3, one must have Sv (e jθ ) = 0. As a consequence, we obtain (3).
In summary, the above derivation
1) extends the feedback control interpretation of feedback communication system over Gaussian channels with
access to feedback and shows how the Youla parameter Q is central to the feedback capacity problem;
2) motivates, in the next section, an explicit construction of capacity-achieving feedback coding schemes from
Youla parameter Q, resolving an issue left open in [2].
A. Construction of Capacity-achieving Feedback Codes from Youla Parameter Q
Once a feasible Q is found for the above optimization (3), which is possibly optimal or arbitrarily close to
optimal, we construct a controller K = Q(1 + Q)−1 as defined in (7) to stabilize the channel within the prescribed
input average power limit.
We next show how to construct a feasible feedback coding scheme from K, a scheme which is deterministic
(time-invariant) and has double exponential decaying decoding error probability. We follow [6].
4 Directed

information, first defined by Massey [31], has been widely used in characterizing the capacity of channels with feedback [32]–[36].

Moreover, it has interpretation on portfolio theory, data compression and hypothesis testing [37].

Fig. 2. Decomposition of controller K into feedback encoder and decoder.

K can always be realized with the following observable and controllable state-space realization (Chapter 3 in
[38]):

A

= s
Xu (k + 1)
0

K:

Xs (k + 1)



h
U(k) = Cs

0



Xs (k)





Bs




 +   Y (k)
Xu (k)
Bu


i Xs (k)
.
Cu 
Xu (k)
Au

(10)

Based on Remark 2, we assume that the eigenvalues of Au are strictly outside the unit disc while the eigenvalues
of As are strictly inside the unit disc. Without loss of generality we assume that As and Au are in Jordan form.
Assume Au has m eigenvalues, denoted by λi (Au ), i = 1, 2, · · · , m.
Starting with the decoder, we decompose K as follows. We present the simplest solution here, others are possible.
This coding scheme is illustrated in Fig.2.

Linear Decoder
The decoder runs K driven by Y .
Xs (k + 1)

=

As Xs (k) + BsY (k),

Xs (0) = 0.

X̂u (k + 1) = Au X̂u (k) + BuY (k),

X̂u (0) = 0.

It produces two signals: an estimate of the negative initial condition of the encoder
X̂u 0 (k) = A−k−1
X̂u (k + 1).
u
and a feedback signal


h
i Xs (k)
.
Û(k) = Cs Cu 
X̂u (k)
Linear Encoder
The encoder runs the following dynamics
X̃u (k + 1) = Au X̃u (k), X̃u (0) = Xu,0 ,
Ũu (k) = Cu X̃u (k),

where Xu,0 represents the message index. It receives Û and produces the channel input
U(k) = Ũu (k) + Û(k)
= Cu (X̃u (k) + X̂u (k)) +Cs Xs (k).
Since the closed loop is stable, U(k) goes to zero with time if the noise is not present. Given the system is
observable, this implies that X̂u (k) → −X̃u (k). Thus, −X̂u 0 (k) is an estimate at time k of X̃u (0) = Xu,0 .
0
−k 0
In the presence of noise, the above coding scheme leads to X̂u 0 (k) v N (−Xu,0 , A−k
u E[X̂u (k)X̂u (k) ](Au ) ) for

large k, where E[X̂u (k)X̂u (k)0 ] represents the state covariance matrix. Note that, since the system is observable and
0
−k 0
controllable, the matrix E[X̂u (k)X̂u (k)0 ] is positive definite and so is A−k
u E[X̂u (k)X̂u (k) ](Au ) (Theorem 4.3 in [6]).

Remark 4: The above proposed coding scheme can be viewed as a generalized Schalkwijk-Kailath scheme. Based
on the Schalkwijk’s scheme in [14], the channel input (encoder) and the message estimate (decoder) for k ≥ 2 can
be summarized as follows by using the notations in this paper.
q
U(k) = A2u − 1Ak−1
u (X̂u 0 (k − 1) + Xu 0 ),
(11)
q
2 − 1Y (k),
X̂u 0 (k) =X̂u 0 (k − 1) − A−k−1
A
u
u
r
P+σw2
where Au =
and σw is the variance of the additive white Gaussian noise in the forward channel 5 . Now, if
σ2
w

in our scheme we let

s
Au =

P + σw2
,
σw2

p
A2u − 1
Bu = −
,
Au

Cu =

q
A2u − 1,

and As = Bs = Cs = 0, then the channel input and the message estimate in our scheme become identical to the
Schalkwijk’s scheme. Specifically, based on the proposed linear coding scheme, we have
U(k) =Ũu (k) + Û(k)
=Cu (X̃u (k) + X̂u (k))
=Cu (Aku Xu 0 + Aku X̂u 0 (k − 1))

(12)

=Cu Aku (Xu 0 + X̂u 0 (k − 1))
q
= A2u − 1Aku (X̂u 0 (k − 1) + Xu 0 ).
X̂u 0 (k) =A−k−1
X̂u (k + 1)
u
=A−k−1
(Au X̂u (k) + BuY (k))
u
−k−1
=A−k
BuY (k)
u X̂u (k) + Au

(13)

=X̂u 0 (k − 1) + A−k−1
BuY (k)
u
q
=X̂u 0 (k − 1) − A−k−2
A2u − 1Y (k).
u
By scaling the message Xu 0 and the corresponding estimate X̂u 0 by factor Au , we recover the dynamics of the
Schalkwijk’s scheme. Note that this constant scaling on the message index Xu 0 have no effect on the reliable
5 We

clarify that in our scheme X̂u 0 (k) is the estimate of −Xu 0 , or −θ in [14], at time instance k.

transmission rate and the power cost at channel input. Furthermore, the first transmission instance (k = 1) of
Schalkwijk’s scheme in [14] is differentiated from the above dynamics such that the scheme is optimal for all
transmission instances. Our scheme, however, is optimal in the steady state which is determined by the above
dynamics.
The next theorem describes how fast messages associated with each Xu,0 are transferred to −X̂u 0 (k) in the presence
of the channel noise.
Theorem 2: Consider an additive Gaussian channel in (6) with Assumption 2 and 3. Given a strictly causal
stable Youla parameter Q(e jθ ) ∈ RH ∞ , the coding scheme described above based on the decomposition of K =
Q(1 + Q)−1 achieves a reliable transmission rate (in the sense of Shannon) at
1
2π

Z π
−π

m

log |1 + Q(e jθ )|dθ = ∑ log |λi (Au )| bits/channel use
i=1

and has double exponential decaying error probability.
The proof is omitted as it directly follows Theorem 4.3 in [6]. Notice that the above achievable rate has the same
form as the objective function in (3), implying that a capacity-achieving feedback code can be constructed from K
if Q is an optimal solution.
In summary, the above discussion provides an approach to construct a feasible feedback coding scheme with
rate

1 Rπ
jθ
2π −π log |1 + Q(e )|dθ

over a stationary finite dimensional Gaussian channel, by leveraging Youla parameter

Q. However, we need first to obtain an optimal Q by solving optimization (3), which is an infinite dimensional
non-convex optimization problem. In the next section, we provide an approach, by solving finite dimensional convex
problems, to find an asymptotic capacity-achieving upper bounds on the capacity. From the filter Q to achieve the
upper bounds, we can construct asymptotically capacity-achieving feedback codes as described in this section.
IV. U PPER B OUNDS ON C APACITY
In this section, we first present an alternative characterization of Gaussian feedback capacity by leveraging the
inverse Fourier transform. Based on this characterization, a sequence of asymptotic capacity-achieving upper bounds
is proposed and evaluated by solving finite dimensional convex optimization problems.
A. Alternative characterization of Gaussian Feedback Capacity
We focus on the optimization problem (3). In what follows, we characterize the Gaussian feedback capacity by
imposing the causality constraints in terms of the inverse Fourier transform.
Lemma 3: Under Assumption 2 and 3, let Q(e jθ ) = a(θ ) + jb(θ ), the feedback capacity can be characterized by
C f b = max
Γ

s.t.

1
2π

1
4π

Z π
−π

Z π

log((1 + a(θ ))2 + b(θ )2 )dθ

−π


a2 (θ ) + b2 (θ ) Sw (θ )dθ ≤ P,

(strict causality constraints in frequency domain)
Z π

Z π

a(θ ) cos(nθ )dθ +
−π

n = 0, 1, 2, · · · , ∞

b(θ ) sin(nθ )dθ = 0
−π

(14)

where the maximum is taken over a functional set Γ defined as
Γ ={a(θ ), b(θ ) : [−π, π] → R

| a(θ ), b(θ ) ∈ L2 }.

(15)

With a bit abuse of notation, Sw (θ ) refers to Sw (e jθ ) for simplicity. The basic idea of this characterization is
that the strict causality can be imposed on the non-positive index coefficients of the inverse Fourier transform of
Q(e jθ ) by setting them to zeros. See Appendix for the detailed proof.
B. Upper bounds
We next obtain upper bounds to C f b by taking into account only a finite number of causality constraints. The
h-upper-bound, denoted by C f b (h), is defined as follows:
Lemma 4: Under Assumption 2 and 3, define C f b (h) with h ∈ Z+ by
C f b (h) = sup
Γ

s.t.

1
2π

1
4π

Z π
−π

Z π

log((1 + a(θ ))2 + b(θ )2 )dθ

−π


a2 (θ ) + b2 (θ ) Sw (θ )dθ ≤ P,
(16)

Z π

Z π

a(θ ) cos(nθ )dθ +
−π

b(θ ) sin(nθ )dθ = 0
−π

n = 0, 1, 2, · · · , h.
Then, C f b (h) ≥ C f b (h + 1),C f b (h) ≥ C f b for any h ≥ 0, and
C f b = lim C f b (h).
h→∞

Proof: C f b (h) is a monotonically non-increasing sequence bounded below, therefore it has a limit. The limit
value must be C f b , otherwise we can construct a strictly causal filter that achieves a rate strictly greater than the
feedback capacity, which is impossible.
The achievability of “supremum” will be proved in the next theorem.
Notice that C f b (h) is still an infinite dimensional problem. To solve C f b (h), we next provide a theorem that
characterizes the Lagrangian dual of C f b (h) and show that there is no duality gap between the infinite dimensional
primal and the finite dimensional dual problems. This result provides a convex optimization approach to compute
C f b (h), a sequence of asymptotically capacity-achieving upper bounds.
Theorem 3: (Main result) Under Assumption 2 and 3, let
A(θ ) = [cos(θ ), cos(2θ ), · · · , cos(hθ )]0 ,
B(θ ) = [sin(θ ), sin(2θ ), · · · , sin(hθ )]0 .
For λ > 0, η ∈ Rh , and η0 ∈ R, define
r2 (θ ) = (2λ Sw (θ ) + η 0 A(θ ) + η0 )2 + (η 0 B(θ ))2 .

Then, the following statements are true.
a) The Lagrangian dual of (16) is given by
(D) : µh = −

max

λ >0,η∈Rh ,η0 ∈R

g(λ , η, η0 )

(17)

where
1
g(λ , η, η0 ) =
2π


r2 (θ )
1
log(2λ Sw (θ ) − ν(θ )) −
+ λ Sw (θ ) dθ − λ P + η0 + .
2
2ν(θ )
2

Z π 
1
−π

with
ν(θ ) =

−r2 (θ ) +

p

r4 (θ ) + 8λ Sw (θ )r2 (θ )
.
2

(18)

(19)

b) The dual problem (D) in (17) is equivalent to the following convex optimization problem
µh = −

max
g̃(λ , η, η0 , ν(θ ))
h
λ > 0, η ∈ R , η0 ∈ R

(20)

∞
ν(θ ) ≥ 0 ∈ C[−π,π]

where
1
g̃(λ , η, η0 , ν(θ )) =
2π


r2 (θ )
1
log(2λ Sw (θ ) − ν(θ )) −
+ λ Sw (θ ) dθ − λ P + η0 + ,
2
2ν(θ )
2

Z π 
1
−π

(21)

and the optimal ν(θ ) is characterized by (19).
c) Furthermore, C f b (h) = µh , and an optimal filter Qh (e jθ ) = a(θ ) + jb(θ ) for C f b (h) exists and is characterized
by
2λ Sw (θ ) + η 0 A(θ ) + η0
−1
ν(θ )
η 0 B(θ )
b(θ ) =
a.e.
ν(θ )
a(θ ) =

a.e.
(22)

where (λ , η 0 , η0 , ν(θ )) are obtained by solving (17) (or (20)).
C. Computing C f b (h)
Although the dual problem of C f b (h) can be cast into a convex optimization (17) with finite number of variables.
The problem is not easily computable since the cost is an integral, not explicitly computable in terms of the variables.
A natural practical approach would be to approximate the integral with a finite sum by discretizing θ . We apply
such discretization to (20) (with spacing

π
m)

and introduce the following finite dimensional convex problem. Given

m, consider
µh (m) = −

max

λ >0,η∈Rh ,η0 ∈R,νi ≥0

g̃m (λ , η, η0 , νi )

where
g̃m (λ , η, η0 , νi ) =



1 2m 1
r2 (θi )
1
log(2λ
S
(θ
)
−
ν
)
+
λ
S
(θ
)
−
−
λ
P
+
η
+
,
w i
i
w i
0
∑ 2
2m i=1
2νi
2

(23)

π
(i − 1).
and θi = −π + m

Under Assumption 1 and 2 on Sw (θ ), we know that for any given feasible (λ , η, η0 , ν(θ )) in (21) the integrand
function, defined in a compact set [−π, π], is bounded (from Assumption 1) and continuous (from Assumption 2)
almost everywhere6 . This implies the Lebesgues criterion for integrablility holds (or equivalently, Riemann integral
holds). Thus, µh (m) is an approximation of µh and limm→∞ µh (m) = µh .
Notice that the optimization (23) is in a simple convex form. In particular, the log of an affine function is
concave.

r2 (θi )
νi

is a quadratic (composed with an affine function of the variables) over linear function, therefore

convex. Thus, (23) can be efficiently solved with standard convex optimization tools, e.g. CVX, a package for
specifying and solving convex programs [39], [40].
Based on the solution to (23), we can actually obtain a guaranteed upper bound on C f b (h) for each m using the
(m)

(m)

upper bound property of dual feasible solutions. Let λ (m) , η (m) , η0 , νi

be the optimal solution to (23). Let

(m)

C f b (m, h) = −g(λ (m) , η (m) , η0 )

(24)

where g(·) is defined in (18). Clearly, C f b (m, h) is computable to arbitrary accuracy.
Corollary 1: Given h ≥ 0, C f b (m, h) ≥ C f b (h) ≥ C f b for ∀m > 0 and
C f b = lim C f b (h) = lim lim C f b (m, h).
h→∞ m→∞

h→∞

V. L OWER B OUNDS ON C APACITY
In the previous section we have introduced a finite dimensional convex optimization (23). From its optimal cost
we were able to obtain a sequence of convergent upper bounds on C f b . In this section, we show that from the
solution to (23) we can obtain lower bounds on C f b by explicitly constructing feedback codes. We show that these
codes provide lower bounds arbitrarily close to the capacity, providing a capacity-achieving feedback code.
The results of this section are summarized as follows.

Constructing Feedback Codes
1) Filter Construction:
Given h ∈ Z+ and m ∈ Z+ with 2m > h. Solve (23) to obtain solution (λh,m , ηh,m , η0,h,m , νi,h,m ).
π
For i = 1, . . . , 2m and θi = −π + m
(i − 1), if νi,h,m > 0 compute ai and bi , from

ai =
bi =

6 Discontinuity

0 A(θ ) + η
2λh,m Sw (θi ) + ηh,m
i
0,h,m

νi,h,m
0 B(θ )
ηh,m
i

νi,h,m

−1

.

may exist when ν(θ ) = 0, but has zero measure. a fact that has been proved in the arguments before (48)

(25)

If νi,h,m = 07 for some i’s, ai , bi can be obtained by solving (70), the dual problem of (23).8
N

Then construct a strictly causal filter Qm
N (z) =

∑ qmn z−n where N = 2m − h − 1 and

n=1
2m

qm
n =

1
∑ ai cos(nθi ) − bi sin(nθi ).
2m i=1

(26)

2) Power Scale:
Let
p :=
Scale Qm
N by αm,N =

1
2π

Z π
−π

jθ 2
jθ
|Qm
N (e )| Sw (e )dθ .

p
m
P/p, i.e., Q̄m
N = αm,N QN , to satisfy the power budget P.

3) Coding Scheme Construction: Construct a feedback coding scheme as described in Section III.A by transm −1 into the state-space representation.9
forming K = Q̄m
N (1 + Q̄N )

Lemma 5: Under Assumptions 2 and 3, and given m, h ∈ Z+ with 2m > h and N = 2m − h − 1, the above coding
scheme achieves a rate (in the sense of Shannon)
RN (m) =

1
2π

Z π
−π

jθ
log |1 + Q̄m
N (e )|dθ .

This result directly follows from Theorem 2. The next theorem shows that the above constructed feedback code
achieves capacity for sufficiently large N and m.
Theorem 4: Under Assumptions 2 and 3, and given h, m ∈ Z+ with 2m > h,
lim lim RN (m) = C f b .

N→∞ m→∞

The proof is rather technical and thus presented in VIII-F. The result of the theorem says that for m and N large
enough there are codes whose rates are arbitrarily close to the feedback capacity. Equivalently, this implies that we
can find m and h = 2m − 1 − N so that the resulting code achieves a rate arbitrarily close to the feedback capacity.
VI. N UMERICAL E XAMPLES
In this section, we provide examples to verify our results. We first summarize the procedure of computing the
feedback capacity and constructing the capacity-achieving codes as follows. Given a finite dimensional Gaussian
channel, we can compute a capacity upper bound C f b (m, h) by solving the finite dimensional convex optimization
(23). Then we construct a feedback coding scheme as described in Section V. One will see that for the examples
the capacity can be evaluated with arbitrary accuracy and the capacity-achieving code can be constructed explicitly.
7 This
8 We

case is not expected to occur for large enough m.

can also compute a j , b j by completion as they together with (25) need to satisfy the causality constraints.
1 2m
1 2m
∑ ai cos(nθi ) + 2m ∑ bi sin(nθi ) = 0,
2m i=1
i=1

9 Technically,

n = 0, 1, · · · , h.

Hankel Singular Value decomposition (H-SVD) procedure [41] can be applied to arbitrarily well approximate an exact state-

space realization of an FIR filter. As widely used to reduce the order of a system realization by such an approximation, one can use H-SVD to
construct a low-order state-space representation of Q̄m
N.

A. First-order moving average Gaussian process
Consider a first-order moving average (i.e. MA(1)) Gaussian process Wi = Ui + 0.4Ui−1 for i ∈ Z, where Ui is a
white Gaussian process with zero mean and unit variance. The power spectral density is Sw (e jθ ) = |1 + 0.4e− jθ |2 .
Given power constraint P = 10, our proposed upper and lower bounds converge to 1.8819 bits/channel use. This
is consistent with that computed from the closed form solution (44) in [2]. In particular, for a stationary noise
process given by
Wi + βWi−1 = Ui + αUi−1 ,
the feedback capacity is
C f b = − log x0 ,
where x0 is the unique positive root of the fourth-order polynomial in x,
Px2 =

(1 − x2 )(1 + σ αx)2
,
(1 + σ β x)

and
σ = sign(β − α).
To obtain a quantitative flavor of the capacity gain of using feedback, we compute the non-feedback capacity as
a comparison, according to the formula (2), (3) in [2]. It turns out that the non-feedback capacity for this channel
is 1.7402 bits/channel use. That is, the capacity gain is 0.1417 bits/channel use by using feedback.
B. Second-order moving average Gaussian process
Consider the following second order moving average Gaussian process with
H(z) = 1 + 0.1z−1 + 0.5z−2
with associated Sw (e jθ ) = |H(e jθ )|2 . While neither the value of capacity or the optimal code are known for this
generalized Gaussian noise, both of them can be efficiently obtained from our approach. With power constraint P =
10, the capacity is evaluated to be 1.9194 bits/channel use (rounded to 4 decimals). Fig. 3 shows the exponentially
fast convergence of the upper and lower bounds in this example. Furthermore, we note that the non-feedback
capacity is 1.7466 bits/channel use, from which we see the noticeable capacity gain by using feedback.
The optimal feedback coding scheme K after order reduction via H-SVD on the finite impulse response is given
by
K=

0.22026(z + 13.84)z2
(z2 + 0.01755z + 0.03498)(z2 + 0.4115z + 3.783)

.

Then, we can obtain the corresponding state-space representation in the real block diagonal form and construct the
explicit coding scheme as shown in Fig. 2. Note that K has two complex conjugate unstable poles at
p1,2 = −0.2057 ± i1.9340,
which would not be easy to find using the approach of [6]. According to Theorem 2, we know the reliable
transmission rate is only determined by the unstable poles of the constructed controller K. Therefore, the fact that

100

Gap [bit/ch.use]

10-2

10-4

10-6

10-8

10-10

10-12

1

2

3

4

5

6

7

h (number of causality constraints)

Fig. 3. Gap between upper and lower bounds exponentially vanishes as the number of causality constraints increases.

probability of error

100

10-1

10-2

10-3
10

15

20

25

30

35

40

code block length k

Fig. 4. Decoding error probability at rate R = 0.95C f b .

the above capacity-achieving coding scheme K is unstable is expected. Also it can be verified from Theorem 2 that
the achievable rate of K is log(|p1 ||p2 |) = 1.9194 bits/channel use.
Moveover, the corresponding optimal closed loop Sensitivity function is
1+Q =

(z2 + 0.01755z + 0.03498)(z2 + 0.4115z + 3.783)
(z2 + 0.1088z + 0.2644)(z2 + 0.1z + 0.5)

As expected, the Sensitivity has the corresponding non-minimum phase zeros at the location of the unstable poles
of K. Note also that the optimal closed loop system includes dynamics that (partially) cancel the noise dynamics.
The term (z2 + 0.1z + 0.5) is the numerator of H(z). This feature is to be consistent with the optimal filter structure
identified in [2].

Fig. 5. Feedback capacity curve for the noise given by (27).

Fig. 4

10

shows that the error probability of the constructed code decays to zero in a double exponential manner.

This implies the coding latency can be very small to achieve certain performance.
C. Third-order ARMA Gaussian process
Finally, Fig. 5 shows the feedback capacity curve as a function of power budget P for an ARMA(3) noise with
H(z) =

z3 − 0.3z2 + 0.5z + 0.2
.
z3 + 0.1z2 + 0.6z + 0.5

(27)

It is noteworthy that the curve is achievable (within any epsilon) by the codes we have proposed and such curve
was unknown until now.
VII. C ONCLUSION
This paper studied the problem of computing the feedback capacity of stationary finite dimensional Gaussian channels and found the asymptotically capacity-achieving codes. Firstly, the interpretation of feedback communication as
feedback control over Gaussian channels was extended by leveraging Youla parameterization. This new interpretation
provides an approach to construct feasible feedback coding schemes with double exponentially decaying error
probability. We next derived an asymptotic capacity-achieving upper bounds, which can be numerically computed
by solving finite dimensional convex optimizations. From the resulting filters that achieve upper bounds, feasible
10 The

simulation procedure, including how to define the message index and determine the detection rule/threshold, can be found in the proof

of Theorem 4.3 in [6].

feedback coding schemes were constructed, leading to a sequence of lower bounds. We proved that the sequence
of lower bounds is asymptotically capacity-achieving.
In summary, the paper provides a computational approach to compute the feedback capacity of Gaussian channels
and construct the capacity-achieving feedback codes, both of which have been open problems in the literature.
Furthermore, we verified our results by numerical examples. In particular, we computed the capacity and constructed
the capacity-achieving feedback codes for a MA(2) channel, and also presented the capacity curve of an ARMA(3)
channel. These results were not known/reported in the literature. Moreover, the resulting structure of the capacityachieving feedback codes, obtained from our control-based numerical approach, could also provide insight for
further investigation on the analytical/closed-form solutions on feedback capacity and capacity-achieving coding
scheme.
VIII. A PPENDIX
A. Proof of Lemma 2
Denote h(·) as the differential entropy, then
→

I (U T → Y T )
T

= ∑ I(U t ;Yt |Y t−1 )
t=0
T

= ∑ h(Yt |Y t−1 ) − h(Yt |Y t−1 ,U t )
t=0
T

= ∑ h(Yt |Y t−1 ) − h(Wt +Ut |Y t−1 ,U t )
t=0
T

= ∑ h(Yt |Y t−1 ) − h(Wt |W t−1 ,U t )
t=0
T

= ∑ h(Yt |Y t−1 ) − h(Wt |W t−1 )
t=0

= h(Y T ) − h(W T )
1
Det[ΣY T ]
,
= log
2
Det[ΣW T ]
where the last step follows from the fact that Y T ,W T are Gaussian vectors and ΣY T denotes the covariance matrix
of sequence Y T . Now, according to (8) (similar to Theorem 4.6 in [6]), it yields
1 → T
I (U → Y T )
T →∞ T
1
Det[ΣY T ]
= lim
log
T →∞ 2T
Det[ΣW T ]
Z
1 π 1
|Y(e jθ )|2
=
log
dθ
2π −π 2
|W(e jθ )|2
Z π
jθ
1
SY (e )
1
=
log
dθ
2π −π 2
Sw (e jθ )
Z π
jθ
1
|1 + Q(e )|2 Sw (e jθ ) + |Qv (e jθ )|2 Sv (e jθ ))
=
log
dθ
4π −π
Sw (e jθ )


Z π
1
|Qv (e jθ )|2 Sv (e jθ )
=
log |1 + Q(e jθ )|2 +
dθ .
4π −π
Sw (e jθ )
lim

The average power formula is derived as follows,
1
T →∞ T
lim

T

∑ Ut2

t=1

1 π
(28)
|Su (e jθ )|2 dθ
=
2π −π
Z π
(b) 1
=
|Q(e jθ )|2 Sw (e jθ ) + |Qv (e jθ )|2 Sv (e jθ )dθ ,
2π −π
where step (a) follows from Parseval’s theorem and step (b) follows from (9) and the independence of W and V .
Z

(a)

B. Proof of Lemma 3
Without loss of generality, we assume a feasible Q(e jθ ) satisfying Q(e jθ ) = Q∗ (e− jθ ). Let Q(e jθ ) = a(θ ) +
− jnθ where a(θ ) = ∞
jb(θ ) = ∑∞
∑n=−∞ cn cos(nθ ) and b(θ ) = ∑∞
n=−∞ cn e
n=−∞ −cn sin(nθ ) represent the real and

imaginary part of Q(e jθ ), respectively. The average power constraint implies a(θ ), b(θ ) ∈ L2 . It is worth noting
that the condition Q(e jθ ) = Q∗ (e− jθ ) indicates Q(e jθ ) is the Fourier transform of a real signal (impulse response
of filter Q(e jθ )) [42]. This implies
a(−θ ) = a(θ ),

b(−θ ) = −b(θ ).

(29)

Next, a strictly causal Q(e jθ ) implies the Fourier coefficients
Z π

cn =

−π

Q(e jθ )e jnθ

dθ
2π

(30)

satisfy cn = 0 for all n ≤ 0. Specifically,
Z π

cn =

(a(θ ) + jb(θ ))(cos(nθ ) + j sin(nθ ))
−π

dθ
=0
2π

(31)

for n ≤ 0. After some elementary algebra and based on the fact that function a(θ ) is even and function b(θ ) is
odd, we have
cn = 0 ⇔

Z π

a(θ ) cos(nθ )dθ −

Z π

−π

b(θ ) sin(nθ )dθ = 0.

(32)

−π

Now we change the sign of index n and have
Z π

c−n =

Z π

a(θ ) cos(nθ )θ +
−π

b(θ ) sin(nθ )dθ = 0

(33)

−π

for n ≥ 0. The proof is complete.
C. Proof of Theorem 3
To facilitate the technical reading, we first provide the high-level discussion on the proof as follows. To derive
the characterization of the dual function, we leverage the optimality conditions when minimizing the Lagrangian
function. These conditions are presented in the form of equations linking the primal and dual variables, based on
which we successfully characterize the minimized Lagrangian function, or equivalently, the dual function. Next, to
show the strong duality, we explicitly construct a primal solution to the primal problem based on the selected dual
variables. Then we show that the dual objective value based on the selected dual variables and the primal objective
value based on the constructed primal solution coincide.

For the ease of derivation, we first change some notations in C f b (h). Let x(θ ) = a(θ ) + 1, y(θ ) = b(θ ). Denote
A(θ ) = [cos(θ ), cos(2θ ), · · · , cos(hθ )]0 , and B(θ ) = [sin(θ ), sin(2θ ), · · · , sin(hθ )]0 . It is true that A(θ ) and B(θ )
have full row rank and moreover the row of A(θ ) and B(θ ) are orthonormal. Consider
(P) :

µhP = infx,y,c

1
4π

Z π

− log(x2 (θ ) + y2 (θ ) + c2 (θ ))dθ

−π

s.t.
Z
1 π
Sw (θ )(x2 (θ ) + y2 (θ ) + c2 (θ ) − 2x(θ ) + 1)dθ ≤ P,
2π Z−π
1 π
(A(θ )x(θ ) + B(θ )y(θ ))dθ = 0,
2π Z−π
π
1
x(θ )dθ = 1, x(θ ), y(θ ), c(θ ) ∈ L2 .
2π −π

(34)

We see that C f b (h) = −µhP if c(θ ) = 0, a.e. Alternatively, µhP can be directly derived from (4). That is, letting
c2 (θ ) =

Sv (e jθ )
,
Sw (e jθ )

Q(e jθ ) = a(θ ) + jb(θ ) = x(θ ) − 1 + jy(θ ) in (4), and following the derivation of Lemma 3 by

respectively replacing a(θ ), b(θ ) with x(θ ) − 1 and y(θ ), we can directly obtain the causality constraints (the last
two equalities) in µhP . Then plugging Q(e jθ ) = x(θ ) − 1 + jy(θ ) into the objective function and the power constraint
in (4) leads to the characterization on µhP . As will be seen, we always have c(θ ) = 0, a.e., for non-flat Gaussian noise.
1) Derivation of dual problem:

Notice that the problem of µhP is neither convex nor quasi-convex. We next construct the dual and analyze
optimality conditions. Consider the Lagrangian
L(x, y, c, λ , η, η0 )
Z 
1 π
1
=
− log(x2 (θ ) + y2 (θ ) + c2 (θ ))
2π −π
2
+λ Sw

(35)

(θ )(x2 (θ ) + y2 (θ ) + c2 (θ ) − 2x(θ ) + 1)

−η 0 (A(θ )x(θ ) + B(θ )y(θ )) − η0 x(θ )] dθ − λ P + η0 ,
where λ ≥ 0, η ∈ Rh , η0 ∈ R. The Lagrangian dual function is
g(λ , η, η0 ) = inf L(x, y, c, λ , η, η0 ).
x,y,c

(36)

Note that for any λ ≥ 0, η ∈ Rh , η0 , g(λ , η, η0 ) provides a lower bound to µhP . The optimal dual objective value is
given as follows,
µhD =

sup g(λ , η, η0 ).

(37)

λ ≥0,η,η0

In sequel, we use “L” to represent “L(x, y, c, λ , η, η0 )” unless specified otherwise. In what follows, we derive the
dual function g(λ , η, η0 ) by leveraging the optimality conditions.

First of all, the optimality conditions for g are obtained by taking functional derivative11 of L(x, y, c, λ , η, η0 )
with respect to function x, y and c. In particular,
(x) :
(y) :
(c) :

x(θ )
x2 (θ )+y2 (θ )+c2 (θ )
y(θ )
x2 (θ )+y2 (θ )+c2 (θ )
c(θ )
x2 (θ )+y2 (θ )+c2 (θ )

= 2λ Sw (θ )x(θ ) − 2λ Sw (θ ) − η 0 A(θ ) − η0
= 2λ Sw (θ )y(θ ) − η 0 B(θ )

(38)

= 2λ Sw (θ )c(θ ).

Note that substituting these conditions back into (35) (hint: taking (x) · x + (y) · y + (c) · c),
Z 
1 π
1
L=
− log(x2 (θ ) + y2 (θ ) + c2 (θ )) + λ Sw (θ )
2π −π
2

−λ Sw (θ )(x2 (θ ) + y2 (θ ) + c2 (θ )) dθ − λ P + η0 + 1,

(39)

which is monotonically decreasing with (x2 (θ ) + y2 (θ ) + c2 (θ )).
Let
ν(θ ) = −

1
+ 2λ Sw (θ ).
x2 (θ ) + y2 (θ ) + c2 (θ )

Then
(x2 (θ ) + y2 (θ ) + c2 (θ )) =

1
,
2λ Sw (θ ) − ν(θ )

2λ Sw (θ ) > ν(θ ).

(40)

Furthermore, the optimality conditions (38) imply that the following equalities need to be satisfied
ν(θ )x(θ ) = 2λ Sw (θ ) + η 0 A(θ ) + η0 ,
ν(θ )y(θ ) = η 0 B(θ ),

(41)

ν(θ )c(θ ) = 0.
For further use, the above assignments yield that
−(2λ Sw (θ ) + η 0 A(θ ) + η0 )x(θ ) + ν(θ )x2
and
− η 0 B(θ )y(θ ) +

2 (θ )

0

= − 2λ Sw (θ )+η2 A(θ )+η0 x(θ ),

ν(θ )y2 (θ )
η 0 B(θ )
=−
y(θ ).
2
2

(42)

(43)

Squaring and summing both sides of (41), we have
ν(θ )2 (x2 (θ ) + y2 (θ ) + c2 (θ )) = r2 (θ ),

(44)

r2 (θ ) = (2λ Sw (θ ) + η 0 A(θ ) + η0 )2 + (η 0 B(θ ))2 , ∀θ ∈ [−π, π].

(45)

where

Substituting (x(θ )2 + y(θ )2 + c2 (θ )) from (40) and rearranging the equation, we have
r2 (θ )(2λ Sw (θ ) − ν(θ )) = ν 2 (θ ),

(46)

which is a second order equation in ν(θ ) (for each θ ) that depends on the dual variables λ , η, η0 . Solving for ν(θ )
we have
ν(θ )1,2 =
11 Note

−r2 (θ ) ±

p

r4 (θ ) + 8λ Sw (θ )r2 (θ )
.
2

that when the space of functions is defined in a Banach space, the functional derivative becomes known as the Fréchet derivative;

when the space of functions is defined in a more general locally convex space, one uses the Gâteaux derivative [43], [44].

From (40), we notice that a larger ν(θ ) leads to a greater (x2 (θ ) + y2 (θ ) + c2 (θ )) that decreases L(x, y, c, λ , η, η0 )
in (39). Therefore, since λ ≥ 0, ν(θ ) ≥ 0 can always be found and we can concentrate on the positive solution
p
−r2 (θ ) + r4 (θ ) + 8λ Sw (θ )r2 (θ )
ν(θ ) =
.
(47)
2
The above derivations show that the desired ν(θ ) with 0 ≤ ν(θ ) < 2λ Sw (θ ) can always be found. Here it is true
that λ > 0. This is because based on (39) the case of λ = 0 can be trivially ruled out as L(x, y, c, λ , η, η0 ) → −∞
by increasing x2 (θ ) + y2 (θ ) + c2 (θ ).
We next remark that, ν(θ ) = 0 for only finite number of θ . The reason is given as follows. Given λ > 0, it is
true that ν(θ ) = 0 for θ = θ̄ if and only if r(θ̄ ) = 0. If r2 (θ ) = 0 for an infinite number of θi , according to (44),
since B(θ ) is full finite column rank, we have η = 0, then 2λ Sw (θi ) + η0 = 0. However, since Sw (θ ) = |H(e jθ )|2
is the spectral density of the output of a finite dimensional LTI stable system driven by white noise, it is a rational
function of cos(θ ). Therefore there must exist a finite number of roots for the equation 2λ Sw (θi ) + η0 = 0, unless
Sw (θ ) is constant, which is excluded by Assumption 3.
Based on this discussion, x(θ ), y(θ ) and c(θ ) satisfying the optimality conditions can be obtained from (41) for
almost all θ , namely,
x(θ ) =

2λ Sw (θ ) + η 0 A(θ ) + η0
,
ν(θ )
y(θ ) =

η 0 B(θ )
,
ν(θ )

c(θ ) = 0,

a.e.

a.e.

(48)
(49)
(50)

a.e.

At those finite number of θ ’s, if any, where ν(θ ) = 0, x and y may not be well-defined and can be discontinuous,
however, these points have zero measure. Furthermore, one can check that, under the above assignments of x(θ ),
y(θ ) and c(θ ) with ν(θ ) in (47), the optimality conditions (38) for the dual function g(λ , η, η0 ) are satisfied.
Note that the above analysis shows that for non-flat channels, the optimal solution must have the feedforward
component c(θ ) = 0 (i.e. Sv (θ ) = 0 a.e.) and therefore the only contribution to the communication rate is obtained
by feedback. From (44) and (46), we have
λ Sw (x2 (θ ) + y2 (θ ) + c2 (θ )) =

1 r2 (θ )
+
.
2 2ν(θ )

(51)

Plugging (44), (46) and (51) into the Lagrangian, the Lagrangian Dual function g(λ , η, η0 ) is characterized by

Z

1
g(λ , η, η0 ) =
2π

π

−π




1
r2 (θ )
1
log(2λ Sw (θ )−ν(θ ))−
+λ Sw (θ ) dθ−λ P+η0 + ,
2
2ν(θ )
2

(52)

and
µhD =max g(λ , η, η0 ),
λ >0,η,η0

(53)

where r2 (θ ) is given by (45) and ν(θ ) given by (47) is a continuous bounded function on a closed set θ ∈ [−π, π],
∞
i.e., ν(θ ) ∈ C[−π,π]
. Notice that “sup” is replaced by “max” due to the existence of solutions in g(λ , η, η0 ), which

is proved in the lemma below.



1
−1
r2 (θ )
∂ g̃
=
+ 2
= 0,
∂ ν(θ ) 2π 2(2λ Sw (θ ) − ν(θ )) 2ν (θ )
1
r2 (θ )
(2λ Sw (θ ) + η 0 A(θ ) + η0 )2 (η 0 B(θ ))2
= 2
=
+
,
2λ Sw (θ ) − ν(θ ) ν (θ )
ν 2 (θ )
ν 2 (θ )

(55)

where r2 (θ ) = (2λ Sw (θ ) + η 0 A(θ ) + η0 )2 + (η 0 B(θ ))2 in (45).


1
2λ Sw (θ ) + η 0 A(θ ) + η0
Sw (θ )
(λ ) :
−2
+ 1 dθ = P
2λ Sw (θ ) − ν(θ )
ν(θ )
−π

Z π 
0
0
2λ Sw (θ ) + A (θ )η + η0
1
B (θ )η
A(θ )
dθ = 0
(η) :
−
+ B(θ )
2π Z−π
ν(θ )
ν(θ )
0
π
1
2λ Sw (θ ) + η A(θ ) + η0
(η0 ) :
−
dθ + 1 = 0
2π −π
ν(θ )
1
(2λ Sw (θ ) + η 0 A(θ ) + η0 )2 (η 0 B(θ ))2
r2 (θ )
(ν(θ )) :
=
+
=
.
2λ Sw (θ ) − ν(θ )
ν 2 (θ )
ν 2 (θ )
ν 2 (θ )
1
2π



Z π

(56)

Lemma 6: Under Assumption 2 and Assumption 3 (i.e. c(θ ) = 0, a.e.), an optimal bounded solution (λ , η, η0 )
to µhD exists.
The proof is given in Appendix VIII-D. As a result, the characterization of dual problem (17) completes. Notice
the sign change µh = −µhD due to the assumption of C f b (h) = −µhP at the beginning of this proof.
The above optimization problem has finite number of variables, however, the cost is an integral function of the
dual variables. To simplify the computation of the dual optimization, we introduce a relaxation of (52) by treating
function ν(θ ) as a free variable, namely,

Z

π

1 log(2λ Sw (θ )−ν(θ )) r2 (θ )
1
−
+λ Sw (θ )dθ−λ P+η0 + .
µ̃hD = max
2π
2
2ν(θ
)
2
λ > 0, η, η
0

(54)

−π

ν(θ ) ≥ 0
∞
ν(θ ) ∈ C[−π,π]

We next show that this relaxation does not lose optimality, i.e., µ̃hD = µhD . Firstly, it is true that µ̃hD ≥ µhD as the
problem of µ̃hD is less constrained. To optimize µ̃hD , we take gradients over λ , η, η0 and functional derivative w.r.t.
function ν(θ ) and equate them to zero. This leads to the equations in (56). Using Leibniz’s rule we can obtain the
gradient w.r.t. to λ , η, η0 . For the functional derivative of the objective function of µ̃hD w.r.t. to ν(θ ), note that the
integral is minimized w.r.t. ν if the integrand (denoted by g̃) is minimized w.r.t. ν(θ ) for all θ , leading to (55) (56).
Note that the last condition in (56) corresponds to (46) and thus (47). This guarantees that ν(θ ) is a continuous
function. In addition, this indicates that given feasible λ > 0, η and η0 , the cost functions of the two optimization
problems coincide (under optimized ν(θ )). Therefore, µ̃hD = µhD .
2) recover the optimal primal value:
Now, we are ready to show the strong duality. The approach is to find a dual solution (λ ∗ , η ∗ , η0∗ , ν ∗ (θ )) based
on which we construct a primal solution (x∗ (θ ), y∗ (θ ), c∗ (θ )). Then we show the dual value and primal value

1
L=
2π

Z π
−π

#
0
 1 (2λ ∗ Sw (θ ) + η ∗ A(θ ) + η0∗ )2 (η ∗0 B(θ ))2
1
∗2
∗2
∗
−
+ λ Sw (θ ) dθ − λ ∗ P + η0∗ .
− log x (θ ) + y (θ ) + −
2
2
2ν ∗ (θ )
2ν ∗ (θ )

"

(61)

coincide.
Assume (λ ∗ , η ∗ , η0∗ , ν ∗ (θ )) is an optimal solution of (52) (or (54)). Then we construct
0

x (θ ) =

2λ ∗ Sw (θ ) + η ∗ A(θ ) + η0∗
,
ν ∗ (θ )

y∗ (θ ) =

η ∗ B(θ )
,
ν ∗ (θ )

∗

(57)

a.e

0

c∗ (θ ) = 0,

(58)

a.e.

(59)

a.e.

Notice that x∗ (θ ), y∗ (θ ), c∗ (θ ) have the same form of (48), (49) and (50), respectively, and, with this association,
we see that (56) coincide with the primal constraints of (P). Moreover, from the fourth expression in (56), we
obtain
x∗2 (θ ) + y∗2 (θ ) =
Since
Rπ

−π x

Rπ

−π log(2λ

∗ S (θ ) − ν ∗ (θ ))dθ
w

∗2 (θ ) + y∗2 (θ )dθ .

1
2λ ∗ Sw (θ ) − ν ∗ (θ )

is bounded, we have bounded

(60)

a.e.

Rπ

−π 2λ

∗ S (θ ) − ν ∗ (θ )dθ
w

and thus bounded

This further indicates x(θ ), y(θ ) ∈ L2 . Therefore x∗ (θ ), y∗ (θ ) are feasible for (P). The

primal cost achieved by this primal feasible solution is
1 π 1
− log(x∗2 (θ ) + y∗2 (θ ))dθ
2π −π 2
Z
1 π 1
=
log(2λ ∗ Sw (θ ) − ν ∗ (θ ))dθ .
2π −π 2
Z

Next we show that
1
2π

Z π
1
−π

2

log(2λ ∗ Sw (θ ) − ν ∗ (θ ))dθ

is equal to the optimal dual cost.
In particular, we use the expressions for x∗ (θ ) and y∗ (θ ) in (57), (58) and take substitutions to obtain (61) with
constraints under the maximum conditions

1 Rπ
∗2
∗2
∗
2π −π Sw (θ ) x (θ ) + y (θ ) − 2x (θ ) + 1 dθ
1 Rπ
∗
∗
2π −π (A(θ )x (θ ) + B(θ )y (θ )) = 0,
1 Rπ ∗
2π −π x (θ )dθ = 1.

= P,

That is,
L=

1
2π

1
− log(x∗2 (θ ) + y∗2 (θ ))dθ + X,
2
−π

Z π

where X is defined in (63). Next we show that X = 0. From equations (42), (43) and (62),


 R

ν ∗ (θ ) ∗2
1 Rπ
1
∗2 (θ )) dθ + λ ∗ 1 π S (θ )(1 − 2x∗ (θ ))dθ − P .
+
(x
(θ
)
+
y
X = 2π
−π 2
2
2π −π w

(62)

#
0
1 (2λ ∗ Sw (θ ) + η ∗ A(θ ) + η0∗ )2 (η∗0 B(θ ))2
∗
−
−
+ λ Sw (θ ) dθ − λ ∗ P + η0∗
2ν(θ )
2ν ∗ (θ )
−π 2
#
"
0
Z
1 π 1 (2λ ∗ Sw (θ ) + η ∗ A(θ ) + η0∗ ) ∗
(η∗0 B(θ )) ∗
∗
=
−
x (θ ) −
y (θ ) + λ Sw (θ ) dθ − λ ∗ P + η0∗ .
2π −π 2
2
2

1
X=
2π

Z π

"

Adding and subtracting λ Sw (θ )(x∗2 (θ ) + y∗2 (θ )) we have (64).

Z 
1 π 1 ν ∗ (θ ) ∗2
∗2
∗
∗2
∗2
X=
+
(x (θ ) + y (θ )) − λ Sw (θ )(x (θ ) + y (θ )) dθ
2π −π 2
2

 Zπ

1
∗2
∗2
∗
∗
Sw (θ ) x (θ ) + y (θ ) − 2x (θ ) + 1 dθ − P .
+λ
2π −π

(63)

(64)

From equations (44), (51) and (62), one can check that X = 0. Therefore the dual cost is equal to the primal cost.
D. Proof of Lemma 6
Recall that the dual problem (53) for the primal problem (34) is given by,
µhD

= supλ ,η,η0 infx,y L(x, y, c, λ , η, η0 )

(65)

= supλ ,η,η0 g(λ , η, η0 ),
where λ > 0, η = [η1 , η2 , · · · , ηh ]0 ∈ Rh , η0 ∈ R and the Lagrangian function L is given by (35) with c(θ ) = 0.
Clearly, µhD is a lower bound on the primal problem (34). Also, it can be verified that the optimal µhD > −∞
by plugging a specific (λ , η, η0 ) into g(λ , η, η0 ), e.g., λ = 1, η = [0, 0, · · · , 0]0 ∈ Rh , η0 = 0. Next, we show by
contradiction that the optimal dual variables are bounded (or exist).
First of all, recall a filter defined by
+∞

Q(e jθ ) = a(θ ) + jb(θ ) =

∑

cn e− jθ .

n=−∞

Let x(θ ) = a(θ ) + 1, y(θ ) = b(θ ). From (34) we know the equality constraints
1
2π

Z π

(A(θ )x(θ ) + B(θ )y(θ ))dθ = 0,
−π

1
2π

Z π

x(θ )dθ = 1
−π

impose zeros on the Fourier coefficients c0 , c−1 , c−2 , · · · , c−h of the filter Q(e jθ ). Specifically,
c0 =
c−n =

1
2π

1
2π

Z π

x(θ )dθ − 1,

−π

Z π

(cos(nθ )x(θ ) + sin(nθ )y(θ ))dθ ,
−π

for n = 1, 2, · · · , h. Also, the constraint
1
P̃(x, y) :=
2π

Z π
−π

Sw (θ )(x2 (θ ) + y2 (θ ) − 2x(θ ) + 1)dθ − P

imposes the power cost of the filter less than P. Equivalently, the Lagrangian function (35) can be rewritten as
L=

1
2π

h
1
− log(x2 (θ ) + y2 (θ ))dθ + λ P̃(x, y) + ∑ ηn c−n + η0 c0 .
2
−π
n=1

Z π

(66)

Now, assume at least one of the optimal dual variables is unbounded, i.e., there exists a sequence ξm =
(λm , ηm , η0,m ) with limm→∞ ||ξm ||2 = ∞ that satisfies
µhD = lim g(ξm ) = lim min L(x, y, ξm ) > −∞.
m→∞

m→∞ x,y

It is worth noting that the case of unbounded λ can be trivially ruled out since P̃(x, y) can be always negative
with a feasible (x, y). Next, for each ξm , we now construct a filter Qm (e jθ ) = x̃m (θ ) − 1 + jỹm (θ ) by assigning a
sequence of Fourier coefficients cn for n ∈ Z where c−n for n = 0, 1, · · · , h has an opposite sign to its corresponding
dual variable as shown in (66). This implies the term ∑hn=1 ηn c−n + η0 c0 → −∞ as m increases. Note that, since the
dual variable λ > 0, we construct the filter satisfying the power constraint by scaling down cn if needed. Therefore,
if limm→∞ ||ξm ||2 = ∞, we have
µhD = lim min L(x, y, ξm ) ≤ lim L(x̃m , ỹm , ξm ) = −∞.
m→∞ x,y

m→∞

This contradicts µhD > −∞. The proof is complete.
E. Proof of Corollary 1
(m)

(m)

From the solution of (23), λ (m) , η (m) , η0 , νi
any feasible dual solution provides a cost

(m)

, we know λ (m) , η (m) , η0

(m)
−g(λ (m) , η (m) , η0 ),

are feasible for (D) in (17), and

which is an upper bound on C f b (h). Namely,

C f b (m, h) ≥ C f b (h).
(m)

(m) ∞
}m=1

Next, consider a sequence {λ (m) , η (m) , η0 , νi

of optimal solutions of (23) for each m. Then, the fact that
(m)

limm→∞ µh (m) = µh implies that limm→∞ maxν(θ )≥0 g̃(λ (m) , η (m) , η0 , ν(θ )) = −µh . Based on the equivalence of
(m)

a) and b) in Theorem 3, we further have limm→∞ g(λ (m) , η (m) , η0 ) = −µh .
Therefore, we have limm→∞ C f b (m, h) = µh . According to c) in Theorem 3, i.e., C f b (h) = µh , and Lemma 4, the
result directly follows.
F. Proof of Theorem 4
The main idea of the proof is the following: we first approximate the feedback capacity characterization (4), C f b ,
by invoking a FIR filter Q of order N and sampling the filter in frequency domain with equal spacing

π
m.

We denote

this two-degree approximation by CNfb (m) and prove that CNfb (m) converges to C f b as N and m increase. Then we
show the constructed code is essentially an optimal solution to CNfb (m), implying that the proposed feedback code
is C f b -achieving as N and m are sufficiently large. We next prove Theorem 4 by providing detailed derivations in
each aforementioned step.

1) A Primal Optimization Problem Equivalent to µh (m):

First of all, we present a finite dimensional optimization and we show that the dual of the optimization is
equivalent to (23). This is a necessary result in order to show that the proposed code construction is essentially an
optimal solution to CNfb (m)

We start from (4), rewritten as follows

C f b = max
U,Q

1
2π

s.t.

1
2π

Z π
1
−π

2



log U(e jθ ) + |1 + Q(e jθ )|2 dθ ,

Z π 
−π


U(e jθ ) + |Q(e jθ )|2 Sw (e jθ )dθ ≤ P,

(67)

U(e jθ ) ≥ 0, Q(e jθ ) ∈ RH 2 , strictly causal,
where U(e jθ ) =

Sv (e jθ )
Sw (e jθ )

since by assumption Sw (e jθ ) > 0 for all θ ∈ [−π, π]. Letting c2 (θ ) =

Sv (e jθ )
Sw (e jθ )

and following

the same steps as in (14), we obtain its generalization as
C f b = max
Γ

1
2π

s.t.

1
4π

Z π
−π

Z π

log(c2 (θ ) + (1 + a(θ ))2 + b(θ )2 )dθ

−π


a2 (θ ) + b2 (θ ) + c2 (θ ) Sw (θ )dθ ≤ P,

Z π

(68)

Z π

a(θ ) cos(nθ )dθ +
−π

b(θ ) sin(nθ )dθ = 0,

n≥0

−π

Γ ={a(θ ), b(θ ), c(θ ) : [−π, π] → R|a(θ ), b(θ ), c(θ ) ∈ L2 }.
Note that we recover (14) if c(θ ) = 0. We are now ready to consider a finite-dimensional approximation of the
above problem. Define
1 2m
∑ log((1 + ai )2 + b2i + c2i )
ai ,bi ,ci ∈R 4m i=1

C f b (m, h) = max
s.t.


1 2m 2
ai + b2i + c2i Sw (θi ) ≤ P,
∑
2m i=1
2m

2m

(69)

1
1
∑ ai cos(nθi ) + 2m ∑ bi sin(nθi ) = 0,
2m i=1
i=1
π
n = 0, 1, 2, · · · , h. θi = −π + (i − 1).
m
As a result of the power constraint, the feasible set of {ai , bi , ci ∈ R} is compact. Therefore, the existence of an
optimal solution is guaranteed.

Lemma 7: (69) is equivalent to the following convex optimization problem:
1 2m
∑ log(Wi )
Wi ,xi ,bi ∈R 4m i=1

C f b (m, h) = max
s.t.

1 2m
∑ (Wi − 2xi + 1) Sw (θi ) ≤ P,
2m i=1
1 2m
∑ xi = 1
2m i=1
1 2m
1 2m
x
cos(nθ
)
+
i
i
∑
∑ bi sin(nθi ) = 0
2m i=1
2m i=1

(70)

n = 1, 2, · · · , h.
π
θi = −π + (i − 1)
m


Wi xi bi




 xi 1 0  ≥ 0.


bi 0 1
Proof: Let xi = 1 + ai and Wi = xi2 + b2i + c2i . Since c2i is free, Wi must equivalently satisfy Wi ≥ xi2 + b2i , which
can be rewritten as an LMI by the Schur complement. The feasible sets of the two problems are equivalent by
simple transformations.
Problem (70) can be solved efficiently using standard interior point methods. It is interesting to note that if the
optimal solution satisfies the Linear Matrix Inequality constraint with equality for all i = 1, . . . , 2m, then the solution
exclusively uses feedback, as the ci ’s, which represent the feedforward component, will be all equal to zero. On
the other hand, it is possible that the optimal solution of these finite-dimensional optimization may require ci 6= 0
for some i. More discussion can be found in the proof of Corollary 2 .
Lemma 8: Given m > 0 and h > 0, µh (m) in (23) is a dual optimization of C f b (m, h) in (69), and µh (m) = C f b (m, h).

Proof: See Appendix VIII-G.
Corollary 2: Given h and m, let ai , bi and ci , i = 1, . . . 2m be the solution to (69), obtained from Wi , xi , bi , the
solution of (70) with ai = xi − 1, and ci = Wi − xi2 − b2i . Let λ , η, η0 , and νi , i = 1, . . . 2m, be the optimal solution to
µh (m) in (23). Then for all i’s such that νi > 0,
2λ Sw (θi ) + η 0 A(θi ) + η0
− 1,
νi
η 0 B(θi )
bi =
.
νi

ai =

(71)

Proof: See Appendix VIII-H.
Remark 5: The reader should notice the similarity with (22) obtained in the semi-infinite dimensional case.
Note however, that the primal solution to C f b (m, h) may require contribution from both feedforward and feedback
components while the solution to C f b (h) is always guaranteed to be exclusively obtained through feedback when
the channel is not flat. Based on Lemma 8, therefore, for fixed h and large enough m the solution to C f b (m, h) will

tend to be exclusively generated from feedback.

2) FIR approximations on C f b :
Since FIR solutions are dense in RH 2 , C f b in (67) can be arbitrarily (uniformly) well approximated by an FIR
Q of order N, for N large enough. This motivates us to introduce
Z


1 π 1
log UN (e jθ ) + |1 + QN (e jθ )|2 dθ ,
CNfb =maxUN ,QN
2π −π 2
Z

1 π
s.t.
UN (e jθ ) + |QN (e jθ )|2 Sw (e jθ )dθ ≤ P,
2π −π

(72)

UN (e jθ ) ≥ 0,
N

QN (e jθ ) =

∑ qk e−k jθ ; FIR order N.
k=1

Remark 6: An optimal solution to the above problem exists since the solution of (4) exists and the set of FIR of
order N is closed in H2 .
Therefore, we have
lim CNfb = C f b .

N→∞

(73)

Now, from the solution to (72) we define
C̃Nfb =

1
2π

Z π
−π

log |1 + QN (e jθ )|dθ .

Note the above integral is well-defined as QN is FIR. Clearly, we have C f b ≥ CNfb ≥ C̃Nfb . Moreover, from (3) and
Remark 1, it follows that for N large enough the contribution to the capacity of the feedforward component UN
must be zero or going to zero (with arbitrarily small power allocation). From the above property, Remark 2 and
the fact that the Fourier series uniformly converge to RL2 functions, we expect 1 + QN to have no zeros on the
unit circle for large enough N. For simplicity we make this a standing assumption. Therefore, we have
lim C̃Nfb −CNfb = 0.

N→∞

(74)

3) Finite-dimensional approximations on C f b :
Firstly, we note that the discretize-then-optimize methodology has been widely investigated and used to solve
optimization problems involving partial differential equations (PDEs) with integration. An overview can be found
in [45], [46] and the reference therein. In what follows, we utilize this methodology in our proof. In particular, a
finite dimensional approximation of C f b in (67) is obtained for large N and by fine frequency discretization of (72),
CNfb (m) =
s.t.



1 2m 1
(e jθi )|2
log Um
(e jθi ) + |1 + Qm
∑
N
N
jθi
m 2m
Um
i=1 2
N (e )≥0,qk
max


1 2m  m jθi
jθi 2
Sw (e jθi ) ≤ P,
(UN (e ) + |Qm
∑
N (e )|
2m i=1
N

jθi
Qm
N (e ) =

∑ qmk e−k jθi ,
k=1

θi = −π +

π
(i − 1).
m

qm
k ∈ R,

(75)

It follows that12
lim CNfb (m) = CNfb ,

m→∞

(77)

According to (73), we conclude that
lim lim CNfb (m) = C f b .

N→∞ m→∞

Similarly to C̃ f b , we introduce13
C̃Nfb (m) =

1 2m
∑ log |1 + QmN (e jθi )|.
2m i=1

(78)

Note that C̃Nfb (m) ≤ CNfb (m). We next show that
lim lim C̃Nfb (m) = C f b .

N→∞ m→∞

(79)

We present the proof as follows. Since we know that limN→∞ limm→∞ CNfb (m) = C f b , we have
lim sup lim sup C̃Nfb (m) ≤ C f b .
N→∞

m→∞

We claim that equality holds. Assume by contradiction that we can always find N and m >> N such that
C f b − C̃Nfb (m) ≥ γ > 0,
while |C f b −CNfb (m)| ≤ ε. C f b − C̃Nfb (m) ≥ γ implies that
PfmN =

1 2m m jθi
∑ UN (e )Sw (e jθi ) ≥ η > 0.
2m i=1

However, omitting some of the tedious steps, which follow similar previous arguments, this would imply that nonzero
feedforward power is required to achieve capacity. Definitely this contradicts the fact that, under our assumption
of non-white noise spectrum, such power allocation cannot be capacity achieving.
Putting (73), (77) and (79) together, we have that, for a given ε > 0, there exist large enough N and m such that
|C f b − C̃Nfb (m)| = |C f b −CNfb +CNfb −CNfb (m) +CNfb (m) − C̃Nfb (m)| < 3ε.
12

(80)

Assume the contrary, i.e., there is always an m for which
either CNfb −CNfb (m) ≥ γ,
or CNfb (m) −CNfb ≥ γ

(76)

for some fixed γ > 0. Take the optimal solution to (72). Since QN (e jθ ) is continuous and bounded, 1 + QN has no zeros on the unit circle,
and UN (e jθ )Sw (e jθ ), (if nonzero), water-fills the modified noise spectrum |1 + QN (e jθ )|2 Sw (e jθ ). Both integrals in (72) can be arbitrarily
well approximated by the corresponding sums since the integrands are Riemann integrable. But this means that within a small perturbation,
(diminishing to zero with m increasing to infinity) we can always find, for all m large enough, a feasible solution to (75) with cost within
any ε > 0, from CNfb , a fact that contradicts the first inequality in (76). To invalidate the second inequality, take the sequence of such m’s
ms
ms jθi
jθi
s N
and call it ms . Assume {qm
k }k=1 leads to 1 + QN with no zeros on the unit circle (otherwise perturb it). UN (e )Sw (e ) will water-fill
ms jθ
jθi 2
jθi
s
|1 + Qm
N (e )| Sw (e ) Now take any δ > 0 and an large enough ms , and consider QN (e ) and the corresponding water-filling solution
jθ
s
Um
N (e ) to a budget power P + δ . This is always possible. Now since the Riemann sums are arbitrarily close to the integrals, the cost will be

arbitrarily close to CNfb + γ. But this is a contradiction, as δ was arbitrary, and we know that with power P the largest rate cannot be strictly
greater than CNfb .
13 We

jθi
m N
assume 1 + Qm
N (e ) 6= 0 for all i. Otherwise, one can apply an arbitrarily small perturbation to {qk }k=1 (resulting in a nearly optimal

jθ
solution) to exclude unit-circle zeros Qm
N (e ). In fact, as a consequence of Remark 2, this assumption always holds for sufficiently large m, N

4) Proof of Theorem 4:

First of all, we need the following lemma to prove our theorem.
Lemma 9: Given N ∈ Z+ with N < 2m, the optimization CNfb (m) in (75) is equivalent to the optimization C f b (m, h)
in (69) with h = 2m − N − 1, i.e.
CNfb (m) = C f b (m, 2m − N − 1).
jθi
2
m jθi
Proof: We start with (75). Let Qm
N (e ) = ai + jbi and ci = UN (e ). The proof follows from the properties

of the DFT applied to the periodic spectrum (using the synthesis and analysis equations in Section 4.21 [47]).
Specifically, we follow the proof of Lemma 3 (in particular, follow the discretized (32)) by imposing zeros on
1 2m
Fourier coefficients qm
∑ ai cos(kθi ) − bi sin(kθi ) for k = N + 1, · · · , 2m.
k =
2m i=1
In what follows, we follow the three-step procedure of feedback codes construction. Based on Lemma 8 and
Lemma 9, we have
CNfb (m) = C f b (m, 2m − N − 1) = µ2m−N−1 (m).
This implies that the feedback code constructed from the optimal solution to µ2m−N−1 (m) is an optimal solution
to CNfb (m). In particular, let a, b be the optimal solution of CNfb (m). a, b can be computed from (69), or directly
from (25) for νi,h,m > 0, as described in the first step of feedback codes construction.
jθ
Then, Qm
N (e ) can be computed from (26) as described in the third step, and requires power

pm
N =

1
2π

Z π
−π

jθ 2
jθ
|Qm
N (e )| Sw (e )dθ .

jθi
In addition, Qm
N (e ) will satisfy the power constraint of (75) and thus,

1 2m m jθi 2
∑ |QN (e )| Sw (e jθi ) = P − PfmN ,
2m i=1
jθi
since there may exist few ci = Um
N (e ) ≥ 0.
m
Then, ηNm = pm
N − (P − Pf N ) must go to zero with m given that N is fixed, as the spectrum of an FIR of order N
jθ
can be arbitrarily well approximated by m samples with m >> N. In addition, the rate for such Qm
N (e ) is given

by
ĈNfb (m) =

1
2π

Z π
−π

jθ
log |1 + Qm
N (e )|dθ ,

and for m > m2 ,
|ĈNfb (m) − C̃Nfb (m)| < ε.
Combined with (80), it yields that for m > max{m0 , m1 , m2 },
|C f b − ĈNfb (m)| < 4ε.
m
N
Now let δNm = pm
N − P. If δN ≤ 0 then Ĉ f b (m) is an achievable rate based on Theorem 2, within 4ε from C f b . If

δNm > 0 then the power constraint is violated, albeit by a negligible amount for large enough N and m.

jθ
We can rescale Qm
N (e ) by αN,m =

q

P
pm
N

=

q

m

δ
m jθ
m
1 − P+δ
m so that αN,m QN (e ) uses power P. When δN ≤ 0,
N

αN,m ≥ 1 and we improve the lower bound, ĈNfb (m), by the rescaling. When δNm > 0, ĈNfb (m) is not a lower bound
jθ
m jθ
on the rate, however, αN,m < 1, and we obtain a lower bound using αN,m Qm
N (e ) in place of QN (e ) .

Based on Theorem 2, the last step of feedback code construction leads to a rate
RN (m) =

1
2π

Z π
−π

jθ
log |1 + αN,m Qm
N (e )|dθ .

We next note that
m
m
m
m
m
δNm = pm
N − P = pN − P + Pf N − Pf N = ηN − Pf N .

Furthermore,
δNm = ηNm − PfmN + Pf N − Pf N ,
|δNm | ≤ |ηNm | + |Pf N | + |PfmN − Pf N |,
where Pf N =

Rπ

−π UN (e

jθ )S (e jθ )dθ .
w

As m → ∞, |ηNm | → 0 and |PfmN − Pf N | → 0. Based on Remark 1, as N → ∞,

|Pf N | → 0. Thus
lim lim δNm = 0, and lim lim αN,m = 1.

N→∞ m→∞

N→∞ m→∞

Therefore, for large enough N and m,
|ĈNfb (m) − RN (m)| < ε,
and
C f b − RN (m) < 5ε.
The proof is complete.
G. Proof of Lemma 8
We consider the dual of (70). Since the Slater’s condition holds and we can always have a feasible primal solution,
e.g.,Wi = xi = 1 and bi = 0, the dual problem has no gap. Moreover, the optimal dual solution exists.
From standard derivations of dual problems involving LMI constraints [48], the dual problem is given by the
following optimization,
C f b (m, h) = −

max

λ ≥0,η0 ,η∈Rh ,ν̃i ,ξ4i ,ξ6i ,ξ5i

gD (λ , η0 , η, νi , ξ4i , ξ6i )

s.t. Ξi ≥ 0, i = 1, . . . , 2m,
r1i = 2λ Sw (θi ) + η 0 A(θi ) + η0 ,

(81)

r2i = η 0 B(θi ),
θi = −π +
where

π
(i − 1),
m


ν̃i


Ξi =  − r21i

− r22i

− r21i
ξ4i
ξ5i

− r22i





ξ5i  ,

ξ6i

(82)

gD (λ , η0 , η, ν̃i , ξ4i , ξ6i ) =

1
1 2m 1
∑ 2 log (2λ Sw (θi ) − 2ν̃i ) + λ Sw (θi ) − ξ4i − ξ6i − λ P + η0 + 2 ,
2m i=1

and A(θ ) and B(θ ) are defined in Theorem 3.
Further analysis leads to the following consequences. If optimal dual λ , η0 , η lead to r1i = r2i = 0 for some i, in
order to maximize gD , we must have the corresponding optimal ν̃i = 0, and the optimal ξ4i = ξ5i = ξ6i = 0. Thus
Ξi = 0. The resulting optimal i-th component of the cost, denoted by gD
i , is then given by


1
1 1
D
log(2λ Sw (θi )) + λ Sw (θi ) − λ P + η0 +
.
gi =
2m 2
2
On the other hand, if the optimal ν̃i is ν̃i > 0 for some i then, using Schur complement on Ξi ≥ 0, the LMI Ξi ≥ 0
is equivalent to



ξ4i

ξ5i

ξ5i

ξ6i





≥

In this case the optimal Ξi is rank one and can
 √
ν̃i

r1i

− √
Ξi = 
 2 ν̃i

r2i
− √
2 ν̃i

− r21i
− r22i


 ν̃i−1

h

− r21i

− r22i .

i

.

be obtained by substituting the above inequality into (82). Namely,


 √

ν̃i



r1i
− √
2 ν̃i

r2i
− √
2 ν̃i


≥ 0.

(83)

The optimal i-th component of the cost is then given by the following expression


r2 (θi )
1
1 1
log(2λ
S
(θ
)
−
2
ν̃
)
+
λ
S
(θ
)
−
−
λ
P
+
η
+
,
gD
=
w i
i
w i
0
i
2m 2
4ν̃i
2
2 + r 2 . Therefore, u (m) in (23) coincides with the cost of (81) and of (70) or (69). Note that ν in
where r2 (θi ) = r1i
n
i
2i

g̃m of (23) is equal to 2ν̃i in gD , while the other dual variables η0 and η are the same between the two problems.
The result then follows from Lemma 7.
H. Proof of Corollary 2
From the proof of Lemma 8, (82) and the complementary slackness condition require that



ν̃i
− r21i − r22i
Wi xi bi



 r1i


 −2
ξ4i
ξ5i   xi 1 0  = 0.



− r22i
ξ5i
ξ6i
bi 0 1
When the optimal dual solution is such that ν̃i > 0, then from (83)


  Wi
√
r1i
r2i

ν̃i − √
− √
 x
2 ν̃i
2 ν̃i  i
bi

this condition reduces to

xi bi


1 0  = 0.

0 1

r1i
r2i
− 1 and bi =
, where νi = 2ν̃i . Since ci = Wi − xi2 − b2i in
νi
νi
Corollary 2, Wi = xi2 + b2i indicates ci = 0. Namely, when the optimal ν̃i > 0, then ci = 0. In other words, the optimal
which when solved leads to Wi = xi2 + b2i , ai =

primal solution is exclusively obtained by feedback and the feedforward contribution must be zero. Moreover, in
this case, the primal optimal solution (25), the components of Q, can be computed from the dual solution. Namely
2λ Sw (θi ) + η 0 A(θi ) + η0
− 1,
νi
η 0 B(θi )
bi =
,
νi

ai =

(84)

where νi = 2ν̃i .
R EFERENCES
[1] C. Li and N. Elia, “Control approach to computing the feedback capacity for stationary finite dimensional Gaussian channels,” 53rd Annual
Allerton Conference on Communication, Control, and Computing (Allerton), pp. 1038–1045, 2015.
[2] Y. H. Kim, “Feedback capacity of stationary Gaussian channels,” IEEE Transactions on Information Theory, vol. 56, no. 1, pp. 57–85,
2010.
[3] B. Francis, “Linear systems and control,” [online] Available: http:// individual.utoronto.ca/ brucefrancis1/ second level/ ECE1638 Francis.
pdf , 2006.
[4] E. J. Hannan and M. Deistler, The Statistical Theory of Linear Systems (Classics in Applied Mathematics). SIAM-Society for Industrial
and Applied Mathematics, 2012.
[5] S.Mitter, “Towards a unified theory of communication and control,” Seminar, LIDS-MIT, 2001.
[6] N. Elia, “When Bode meets Shannon: control-oriented feedback communication schemes,” IEEE Transactions on Automatic Control,
vol. 49, no. 9, pp. 1477–1488, 2004.
[7] J. L. Liu and N. Elia, “Convergence of fundamental limitations in feedback communication, estimation, and feedback control over Gaussian
channels,” Communications in Information and Systems (CIS), vol. 14, no. 3, pp. 161–211, 2014.
[8] O. Shayevitz and M. Feder, “Optimal feedback communication via posterior matching,” IEEE Transactions on Information Theory, vol. 57,
no. 3, pp. 1186 – 1222, 2011.
[9] E. Ardestanizadeh and M. Franceschetti, “Control-theoretic approach to communication with feedback,” IEEE Transactions on Automatic
Control, vol. 57, no. 10, pp. 2576–2587, 2012.
[10] E. Ardestanizadeh, P. Minero, and M. Franceschetti, “Lqg control approach to Gaussian broadcast channels with feedback,” IEEE
Transactions on Information Theory, vol. 58, no. 8, pp. 5265–5278, 2012.
[11] P. Elias, “Channel capacity without coding,” MIT Research Lab. of Electronics, Cambridge, MA, Quarterly Progr. Rep., 1956.
[12] ——, “Networks of Gaussian channels with applications to feedback systems,” IEEE Transactions on Information Theory, vol. 13, pp.
493 – 501, 1967.
[13] J. P. M. Schalkwijk and T. Kailath, “A coding scheme for additive noise channels with feedback, part i: no bandwidth constraint,” IEEE
Transactions on Information Theory, vol. IT-12, no. 2, pp. 172–182, 1966.
[14] J. P. M. Schalkwijk, “A coding scheme for additive noise channels with feedback, part ii: band-limited signals,” IEEE Transactions on
Information Theory, vol. IT-12, no. 2, pp. 183–189, 1966.
[15] S. Butman, “A general formulation of linear feedback communication systems with solutions,” IEEE Transactions on Information Theory,
vol. 15, no. 3, pp. 392 – 400, 1969.
[16] ——, “Linear feedback rate bounds for regressive channels,” IEEE Transactions on Information Theory, vol. 22, no. 3, pp. 363 – 366,
1976.
[17] J. C. Tiernan and J. P. M. Schalkwijk, “An upper bound to the capacity of the band-limited Gaussian autoregressive channel with noiseless
feedback,” IEEE Transactions on Information Theory, vol. 20, pp. 311 – 316, 1974.
[18] J. C. Tiernan, “Analysis of the optimum linear system for the autoregressive forward channel with noiseless feedback,” IEEE Transactions
on Information Theory, vol. 22, pp. 359 – 363, 1976.
[19] J. Wolfowitz, “Signalling over a Gaussian channel with feedback and autoregressive noise,” J. Appl. Probability, vol. 12, no. 4, pp. 713 –
723, 1975.
[20] L. H. Ozarow, “Random coding for additive Gaussian channels with feedback,” IEEE Transactions on Information Theory, vol. 36, no. 1,
pp. 17 – 22, 1990.

[21] ——, “Upper bounds on the capacity of Gaussian channels with feedback,” IEEE Transactions on Information Theory, vol. 36, no. 1, pp.
151 – 161, 1990.
[22] T. M. Cover and S. Pombra, “Gaussian feedback capacity,” IEEE Transactions on Information Theory, vol. 35, no. 1, pp. 37–43, 1989.
[23] J. Ham and D. J. Love, “On the achievable rate of the additive Gaussian noise channel with noisy feedback,” IEEE International Symposium
on Information Theory, pp. 226–230, 2010.
[24] C. Li and N. Elia, “Bounds on the achievable rate of noisy feedback Gaussian channels under linear feedback coding scheme,” IEEE
International Symposium on Information Theory, pp. 169–173, 2011.
[25] ——, “Upper bound on the capacity of Gaussian channels with noisy feedback,” 49th Annual Allerton Conference, pp. 84–89, 2011.
[26] ——, “Noisy feedback communications with side information at the decoder,” 50th Annual Allerton Conference, pp. 1856–1863, 2012.
[27] C. Li, “Fundamental limitations on communication channels with noisy feedback: information flow, capacity and bounds,” Iowa State
University Ph.D. dissertation, 2013.
[28] J. L. Liu, N. Elia, and S. Tatikonda, “Capacity-achieving feedback scheme for markov channels with channel state information,” IEEE
International Symposium on Information Theory, pp. 71–75, 2004.
[29] S. Yang, “On the feedback capacity of power-constrained Gaussian noise channels with memory,” IEEE Transactions on Information
Theory, vol. 53, no. 3, pp. 929 – 954, 2007.
[30] J. C. Doyle, B. A. Francis, and A. R. Tannenbaum, Feedback Control Theory.

New York: Macmillan, 1992.

[31] J. L. Massey, “Causality, feedback and directed information,” In Proc. Int. Symp. Inf. Theory Applic., pp. 303–305, 1990.
[32] S. Yang, A. Kavcić, and S. Tatikonda, “Feedback capacity of finite state machine channels,” IEEE Transactions on Information Theory,
vol. 51, no. 3, pp. 799 – 810, 2005.
[33] Y. H. Kim, “A coding theorem for a class of stationary channels with feedback,” IEEE Transactions on Information Theory, vol. 54, no. 4,
pp. 1488–1499, 2011.
[34] S. Tatikonda and S. Mitter, “The capacity of channels with feedback,” IEEE Transactions on Information Theory, vol. 55, no. 1, pp.
323–349, 2009.
[35] C. Li and N. Elia, “The information theoretic characterization of the capacity of channels with noisy feedback,” IEEE International
Symposium on Information Theory, pp. 174–178, 2011.
[36] G. Kramer, “Directed information for channels with feedback,” Ph.D. dissertation, Swiss Federal Institute of Technology (ETH), Zurich,
Switzeland, 1998.
[37] H. H. Permuter, Y. H. Kim, and T. Weissman, “Interpretations of directed information in portfolio theory, data compression, and hypothesis
testing,” IEEE Transactions on Information Theory, vol. 57, no. 6, pp. 3248–3259, 2011.
[38] W. L. Brogan, Modern Control Theory.

Prentice Hall, 1991.

[39] M. Grant and S. Boyd, “Cvx: Matlab software for disciplined convex programming, version 2.0 beta,” 2013.
[40] ——, “Graph implementations for nonsmooth convex programs, recent advances in learning and control,” Lecture Notes in Control and
Information Sciences, Springer, pp. 95–110, 2008.
[41] S. Kung, “A new identification and model reduction algorithm via singular value decompositions,” Conf. on Circuits, Systems and Computers,
no. 6-8, pp. 705–714, 1978.
[42] D. Rowell, “Determining a system’s causality from its frequency response,” Lecture notes, MIT Signal Processing Course 2.161, 2009.
[43] D. Behmardi and E. D. Nayeri, “Introduction of Fréchet and Gâteaux derivative,” Applied Mathematical Sciences, vol. 2, no. 20, pp.
975–980, 2008.
[44] D. G. Luenberger, Optimization by Vector Space Methods.

John Wiley & Sons, Inc, 1969.

[45] K. Ghobadi, N. Nedialkov, and T. Terlaky, “On the discretize then optimizae approach,” Industrial and Systems Engineering, [online],
https://ise.lehigh.edu/sites/ise.lehigh.edu/files/09t 005.pdf, 2014.
[46] J. Liu and Z. Wang, “Non-commutative discretize-then-optimize algorithms for elliptic pde-constrained optimal control problems,” [online].
https://arxiv.org/abs/1706.07652, 2017.
[47] J. G. Proakis and D. G. Manolakis, Digital Signal Processing, 4th Edition.
[48] S. Boyd and L. Vandenberghe, Convex Optimization.

Pearson Prentice Hall, 2007.

Cambridge: Cambridge university, 2004.

