---
abstract: 'Using physical layer network coding, compute-and-forward is a promising relaying scheme that effectively exploits the interference between users and thus achieves high rates. In this paper, we consider the problem of finding the optimal integer-valued coefficient vector for a relay in the compute-and-forward scheme to maximize the computation rate at that relay. Although this problem turns out to be a shortest vector problem, which is suspected to be NP-hard, we show that it can be relaxed to a series of equality-constrained quadratic programmings. The solutions of the relaxed problems serve as real-valued approximations of the optimal coefficient vector, and are quantized to a set of integer-valued vectors, from which a coefficient vector is selected. The key to the efficiency of our method is that the closed-form expressions of the real-valued approximations can be derived with the Lagrange multiplier method. Numerical results demonstrate that compared with the existing methods, our method offers comparable rates at an impressively low complexity.'
author:
- 'Baojian Zhou,  Jinming Wen, and Wai Ho Mow, [^1][^2]'
bibliography:
- 'C:/Users/bzhouab/Dropbox/Research/ref/cf\_qpr.bib'
title: |
    A Quadratic Programming Relaxation Approach to Compute-and-Forward\
    Network Coding Design
---

physical layer network coding, AWGN networks, compute-and-forward, quadratic programming, Lagrange multiplier.

Introduction
============

a promising relaying strategy in wireless networks, compute-and-forward (CF) has attracted a lot of research interest since it was proposed in 2008 by Nazer and Gastpar [@Nazer2008]. The advantage of CF is that it achieves higher rates in the medium signal-to-noise ratio (SNR) regime when compared with other relaying strategies, e.g., amplify-and-forward, decode-and-forward. Relays in CF attempt to decode integer linear combinations of the transmitted codewords, rather than the codewords themselves. The integer coefficient vectors corresponding to the linear combinations and the decoded messages are then forwarded to the destination. Under certain conditions (see [@Nazer2011] for details), the destination can recover the original source messages with enough forwarded messages and coefficient vectors from the relay.

The design of the CF scheme lies in selecting the coefficient vectors at the relays. There are many choices of coefficient vectors for one relay, and each may render a different *computation rate* [@Nazer2011] at that relay. Computation rate is defined as the maximum transmission rate from the associated sources to a relay such that the linear combinations at the relay can be successfully decoded. If the coefficient vectors are linearly independent, the *achievable rate* of the network equals the minimum computation rate; otherwise, the achievable rate of the network is zero and none of the messages can be recovered. The objective of designing the CF scheme is to maximize the overall achievable rate of the network. One approximation method is presented in [@Wei2012]. However, for those networks where each relay is allowed to send only one coefficient vector to the destination, and only local channel state information (CSI) is available, i.e., each relay knows merely its own channel vector, one reasonable solution is to select the coefficient vector that maximizes the computation rate at each relay.

In this paper, we consider additive white Gaussian noise (AWGN) networks where only local CSI is available, and focus on the CF network coding design problem with the objective being maximizing the computation rate at a relay by choosing the optimal coefficient vector. It has been shown that this problem reduces to a shortest vector problem (SVP). Different methods have been developed to tackle this problem. The branch-and-bound method proposed in [@Richter2012] finds the optimal solution but its efficiency degrades as the dimension of channel vectors grows according to the simulation results. Although the general SVP is suspected to be NP-hard, Sahraei and Gastpar showed in [@Sahraei2014] that the SVP in the CF design is special, and developed an algorithm (called the “SG” method in this paper) that solves the SVP in polynomial time. For independent and identically distributed (i.i.d.) Gaussian channel entries, the complexity of the SG method is of order 2.5 with respect to the dimension, and is linear with respect to the square root of the signal-to-noise ratio (SNR). A class of methods are those based on lattice reduction (LR) algorithms (e.g., Minkowski, HKZ, LLL, and CLLL LR algorithms; c.f. [@Zhang2012; @Lenstra1982; @Vetter2009; @Ling2013; @Chang2013; @Gan2009]). The method in [@Sakzad2012] based on the LLL LR algorithm [@Lenstra1982] provides close-to-optimal rates and is well-known to be of polynomial time complexity with respect to the vector dimension. In [@Wen2015], we proposed an efficient method based on sphere decoding to find the optimal coefficient vector; however there is no theoretical guarantee on the complexity.

Our goal in this work is to develop a new method that finds a suboptimal coefficient vector for a relay with low complexity compared with the existing methods, while provides a close-to-optimal computation rate at the same time. Taking advantage of some useful properties of the problem, we first show that the original SVP can be approximated by a series of quadratic programmings (QPs). The closed-form solutions of the QPs are derived by use of the Lagrange multiplier method and can be computed with linear complexity with respect to the dimension, which is the key to the efficiency of our method. The solutions of the QPs serve as real-valued approximations of the integer coefficient vector, and are quantized into a set of candidate integer vectors by a successive quantization algorithm. Finally, the integer vector in the candidate set that maximizes the computation rate is selected to be the coefficient vector. The complexity of our method is of order 1.5 with respect to the dimension for i.i.d. Gaussian channel entries, and is lower than the above mentioned methods. Numerical results demonstrate that among existing methods that provide close-to-optimal rates, our method is much more efficient as expected.

As a summary, our contributions in this work include the following:

-   For the real-valued channels, we develop a quadratic programming relaxation approach to find a suboptimal coefficient vector for a relay so that the computation rate at that relay is close-to-optimal. The complexity is $O(L\sqrt{P\norm{\h}^2})$ for a given channel vector $\h\in\Rbb^L$ and signal power constraint $P$, and is of average value $O(P^{0.5}L^{1.5})$ for i.i.d. standard Gaussian channel entries.

-   For the complex-valued channels, we demonstrate how to apply our method in an efficient way to find the complex-valued coefficient vector.

-   Extensive simulation results are presented to compare the effectiveness and efficiency of our method with the existing methods.

Part of this work has been presented in [@Zhou2014]. One main improvement here is the complexity order for i.i.d. Gaussian channel entries is further reduced from 3 to 1.5.

In the following, we will first introduce the system model of AWGN networks as well as the CF network coding design problem in Section \[section:ProblemStatement\]. Then in Section \[section:ProposedMethod\], we will present our proposed method in detail. Numerical results will be shown in Section \[section:NumericalResults\]. Finally, we will conclude our work in Section \[section:Conclusions\].

[*Notation.*]{} Let $\Rbb$ be the real field, $\Cbb$ be the complex field, and $\Zbb$ be the ring of integers. Boldface lowercase letters denote column vectors, and boldface uppercase letters denote matrices, e.g., $\w\in\Rbb^L$ and $\W\in\Rbb^{M\times L}$. $\norm{\w}$ denotes the $\ell^2$-norm of $\w$, and $\w^T$ denotes the transpose of $\w$. For a vector $\w$, let $\w(\ell)$ be the element with index $\ell$, and $\w(i\!:\!j)$ be the vector composed of elements with indices from $i$ to $j$. For a matrix $\W$, let be the submatrix containing elements with row indices from $i$ to $j$ and column indices from $k$ to $\ell$, be the submatrix containing elements with row indices from $i$ to $j$ and column index $k$, $\W(i,k\!:\!\ell)$ be the submatrix containing elements with row index $i$ and column indices from $k$ to $\ell$, and $\W(i,j)$ be the element with row index $i$ and column index $j$. Let $\floor{x}$ and $\ceil{x}$, i.e., the corresponding floor and ceiling functions of $x$, be the maximum integer no greater than $x$ and the minimum integer no less than $x$, respectively. Let $\floor{\w}_\ell$ and $\ceil{\w}_\ell$ be the vectors generated from $\w$ by applying the corresponding operation on the $\ell$-th element only. $\0$ denotes an all-zero vector, and $\I$ denotes an identity matrix. ${\rm sign}(\w)$ returns the vector that contains the signs of the elements in $\w$. ${\rm abs}(\w)$ returns the vector whose elements are the absolute values of the elements in $\w$.

Problem Statement {#section:ProblemStatement}
=================

We consider additive white Gaussian noise (AWGN) networks [@Nazer2011] where sources, relays and destinations are connected by linear channels with AWGN. For the ease of explanation, we first develop our method for real-valued channels, and then demonstrate how to apply our method to complex-valued channels. An AWGN network with real-valued channels is defined as the following.

\[definition:RealChannelModel\] *(Real-Valued Channel Model)* In an AWGN network, each relay (indexed by $m=1,2,\cdots,M$) observes a noisy linear combination of the transmitted signals through the channel, $$\begin{aligned}
\label{equation:RealChannelModel}
\y_m = \sum_{\ell=1}^L \h_m(\ell)\x_\ell + \z_m,\end{aligned}$$ where $\x_\ell \in \Rbb^n$ with the power constraint $\frac{1}{n}\norm{\x_\ell}^2 \leq P$ is the transmitted codeword from source $\ell$ ($\ell = 1,2,\cdots,L$), $\h_m \in \Rbb^L$ is the channel vector to relay $m$, $\h_m(\ell)$ is the $\ell$-th entry of $\h_m$, $\z_m \in \Rbb^n$ is the noise vector with entries being i.i.d. Gaussian, i.e., $\z_m\!\sim\!\bigN\!\left(\0,\I\right)$, and $\y_m$ is the signal received at relay $m$.

In the sequel, we will focus on one relay and thus ignore the subscript “$m$” in $\h_m$, $\a_m$, etc.

In CF, rather than directly decode the received signal $\y$ as a codeword, a relay first applies to $\y$ an amplifying factor $\alpha$ such that $\alpha\h$ is close to an integer *coefficient vector* $\a$, and tries to decode $\alpha\y$ as an integer linear combination, whose coefficients form $\a$, of the original codewords $\{\x_\ell\}$. The *computation rate* [@Nazer2011] is the maximum transmission rate from the associated sources to a relay such that the integer linear combinations at the relay can be decoded with arbitrarily small error probability. Assume the $\log$ function is with respect to base 2, and define $\log^+(w) \triangleq \max \left( \log(w),0 \right)$. The computation rate can be calculated with Theorem \[theorem:ComputationRate\] from [@Nazer2011].

*(Computation Rate in Real-Valued Channel Model)* \[theorem:ComputationRate\] For a relay with coefficient vector $\a$ in the real-valued channel model defined in Definition \[definition:RealChannelModel\], the following computation rate is achievable,

$$\begin{aligned}
\label{equation:ComputationRate}
\bigR \left( \h, \a \right) =
\frac{1}{2} \log^+
\left( \left(
\norm{\a}^2 -
\frac{P (\h^T\a)^2}{1 + P\norm{\h}^2}
\right)^{-1}\right).
\end{aligned}$$

With the computation rate being the metric, we define the optimal coefficient vector as follows.

*(The Optimal Coefficient Vector)* The optimal coefficient vector ${\a^\star}$ for a channel vector $\h$ is the one that maximizes the computation rate, $$\begin{aligned}
\label{equation:aOptimal}
{\a^\star}= \arg \max_{\a \in \Zbb^L \backslash\{\0\}} \bigR \left( \h, \a \right).
\end{aligned}$$

After a few simple manipulations, the optimization problem stated in can be written in the following quadratic form [@Wei2012], $$\label{equation:IntegerQP}
\begin{aligned}
{\a^\star}=
\arg\min_{\a\in\Zbb^L\backslash\{\0\}}
\a^T\bsG\a,
\end{aligned}$$ where $$\begin{aligned}
\label{equation:G}
\bsG \triangleq
\I - \frac{P}{1 + P\norm{\h}^2} \h\h^T.
\end{aligned}$$

If we take $\bsG$, which is positive definite, as the *Gram matrix* of a lattice $\Lambda$, then the problem turns out to be the SVP in the lattice $\Lambda$. In the next section, we will propose an efficient approximation method based on QP relaxation that gives a suboptimal coefficient vector.

Proposed Method {#section:ProposedMethod}
===============

In this section, we will first derive our method for the real-valued channel model, and then extend the method for the complex-valued channel model.

Preliminaries
-------------

We start with investigating some properties of the problem, which is the basis of our new method.

(Signature Matrix) A signature matrix is a diagonal matrix whose diagonal elements are $\pm1$.

(Signed Permutation Matrix) A signed permutation matrix is a generalized permutation matrix whose nonzero entries are $\pm 1$.

After replacing $-1$’s with $1$’s, a signed permutation matrix becomes a permutation matrix. Obviously, signed permutation matrices are unimodular and orthogonal. Every signed permutation matrix can be expressed as $\S=\P\T$, where $\T$ is a signature matrix, and $\P$ is a permutation matrix.

\[theorem:ProblemTransformation\] If ${\a^\star}$ is the optimal coefficient vector for a channel vector $\h$ with power constraint $P$, then for any signed permutation matrix $\S\in\Zbb^{L\times L}$, $\S{\a^\star}$ is optimal for $\S\h$ with the same power constraint $P$, and $\bigR(\h,{\a^\star})=\bigR(\S\h,\S{\a^\star})$.

We first show $\bigR \left( \h, \a \right) = \bigR \left( \S\h, \S\a \right)$ for any $\h$ and $\a$ with the same power constraint $P$. $\S$ is unimodular, then $\S\a$ is an integer vector and can be applied as a coefficient vector. $\S$ is orthogonal, then $\S^T\S = \I$. $\norm{\S\h}^2 = \h^T\S^T\S\h = \h^T\h = \norm{\h}^2$, and similarly $\norm{\S\a}^2 = \norm{\a}^2$. $(\S\h)^T\S\a = \h^T\S^T\S\a = \h^T\a$. According to Theorem \[theorem:ComputationRate\], the computation rate $\bigR \left( \h, \a \right)$ is determined by $P$, $\norm{\h}^2$, $\norm{\a}^2$, and $\h^T\a$. Thus, $\bigR \left( \h, \a \right) = \bigR \left( \S\h, \S\a \right)$.

${\a^\star}$ is optimal for $\h$ means ${\a^\star}$ maximizes $\bigR \left( \h, \a \right)$. Then $\S{\a^\star}$ maximizes $\bigR \left( \S\h, \S\a \right)$ since $\bigR \left( \h, \a \right) = \bigR \left( \S\h, \S\a \right)$ always holds. Therefore, $\S{\a^\star}$ is optimal for $\S\h$ with the same power constraint $P$.

(Nonnegative Ordered Vector) A vector $\h$ is said to be nonnegative ordered if its elements are nonnegative and in nondecreasing order according to their indices.

\[lemma:h2hBar\] For any vector $\h$, there exists a signed permutation matrix $\S$ such that $\S\h$ is nonnegative ordered.

To find such an $\S$ in Lemma \[lemma:h2hBar\], we can simply choose $\S=\P\T$, where $\T$ is a signature matrix that converts all the elements in $\h$ to nonnegative, and $\P$ is a permutation matrix that sorts the elements in $\T\h$ in nondecreasing order.

With Theorem \[theorem:ProblemTransformation\] and Lemma \[lemma:h2hBar\], for any channel vector $\h$, we can first find a signed permuation matrix $\S$ and transform $\h$ to the nonnegative ordered ${\bar{\h}}=\S\h$, then obtain the optimal coefficient vector ${{\bar{\a}}^\star}$ for ${\bar{\h}}$, and finally recover the desired optimal coefficient vector ${\a^\star}=\S^{-1}{{\bar{\a}}^\star}$ for $\h$. In this way, it suffices to focus on solving the problem in for nonnegative ordered channel vectors ${\bar{\h}}$.

\[remark:Transformation\] In implementation, there is no need to use the signed permutation matrix $\S$. It is merely necessary to: 1) record the sign of the elements in $\h$ with a vector $\bst=\fsign{\h}$, and 2) sort $\boldsymbol{\hbar}=\fabs{\h}$ in ascending order as ${\bar{\h}}$ and record the original indices of the elements with a vector $\p$ such that ${\bar{\h}}(\ell)=\boldsymbol{\hbar}(\p(\ell))$, $\ell=1,2,\cdots,L$. After ${{\bar{\a}}^\star}$ for ${\bar{\h}}$ is obtained, ${\a^\star}$ for $\h$ can be recovered with ${\a^\star}(\p(\ell))=\bst(\p(\ell)){{\bar{\a}}^\star}(\ell)$, $\ell=1,2,\cdots,L$.

Given a channel vector as $\h=[-1.9,0.1,1.1]^T$, then $\bst=[-1,1,1]^T$, $\fabs{\h}=[1.9,0.1,1.1]^T$, ${\bar{\h}}=[0.1,1.1,1.9]^T$, and $\p=[2,3,1]^T$. If for certain power $P$, ${{\bar{\a}}^\star}=[0,1,2]^T$, then ${\a^\star}=[-2,0,1]^T$.

According to , if ${\a^\star}$ is optimal for $\h$, then $-{\a^\star}$ is also optimal for $\h$. To reduce redundancy, we restrict the optimal coefficient vector ${\a^\star}$ to be the one such that $\h^T{\a^\star}\geq0$ in the following.

\[lemma:aNonnegative\] If all the elements in a channel vector $\h$ are nonnegative, then all the elements in the optimal coefficient vector ${\a^\star}$ are also nonnegative.

Suppose ${\a^\star}(i) < 0$, and define $\a'$ as: $\a'(i) = 0$, and $\a'(\ell) = {\a^\star}(\ell)$, $\forall \ell \neq i$. Obviously, $\norm{\a'} < \norm{{\a^\star}}$, and $\h^T \a' \geq \h^T {\a^\star}\geq 0$. Then according to , $\bigR \left( \h, \a' \right) > \bigR \left( \h, {\a^\star}\right)$, which implies ${\a^\star}$ is not optimal and leads to a contradiction. Thus, all the elements in ${\a^\star}$ must be nonnegative.

\[lemma:a0\] For a channel vector $\h$ and its corresponding optimal coefficient vector ${\a^\star}$, if $\h(i) = 0$, then ${\a^\star}(i) = 0$.

Suppose $\h(i) = 0$, and ${\a^\star}(i) \neq 0$. Define $\a'$ as: $\a'(i) = 0$, and $\a'(\ell) = {\a^\star}(\ell)$, $\forall \ell \neq i$. Obviously, $\norm{\a'} < \norm{{\a^\star}}$, and $\h^T \a' = \h^T {\a^\star}\geq 0$. Then according to , $\bigR \left( \h, \a' \right) > \bigR \left( \h, {\a^\star}\right)$, which implies ${\a^\star}$ is not optimal. Thus, if $\h(i) = 0$, then ${\a^\star}(i) = 0$.

\[lemma:hijEqual\] For a channel vector $\h$ and its corresponding optimal coefficient vector ${\a^\star}$, if $\h(i) = \h(j)$, $i < j$, then ${\a^\star}(i) = {\a^\star}(j)$ or $\fabs{{\a^\star}(i) - {\a^\star}(j)} = 1$.

Without loss of generality, assume ${\a^\star}(i) - {\a^\star}(j) < -1$. Define $\a'$ as: $\a'(i) = {\a^\star}(i)+1$, $\a'(j) = {\a^\star}(j)-1$, and $\a'(\ell) = {\a^\star}(\ell)$, $\forall \ell \notin \{i,j\}$. Obviously, $\norm{\a'} < \norm{{\a^\star}}$, and $\h^T \a' = \h^T {\a^\star}\geq 0$. Then according to , $\bigR \left( \h, \a' \right) > \bigR \left( \h, {\a^\star}\right)$, which implies ${\a^\star}$ is not optimal. Thus, ${\a^\star}(i) - {\a^\star}(j) \geq -1$. Similarly, ${\a^\star}(j) - {\a^\star}(i) \geq -1$. Therefore, ${\a^\star}(i) = {\a^\star}(j)$ or $\fabs{{\a^\star}(i) - {\a^\star}(j)}=~1$.

\[remark:hijEqual\] In Lemma \[lemma:hijEqual\], for the case where $\h(i) = \h(j)$ with $\fabs{{\a^\star}(i) - {\a^\star}(j)} = 1$, $i < j$, we will always set ${\a^\star}(j) = {\a^\star}(i) + 1$ since setting ${\a^\star}(i) = {\a^\star}(j) + 1$ results the same computation rate. Then, as long as $\h(i) = \h(j)$, $i < j$, it holds that ${\a^\star}(i)\leq{\a^\star}(j)$.

\[theorem:aNonnegativeOrdered\] For a nonnegative ordered channel vector $\h$, the optimal coefficient vector ${\a^\star}$ is also nonnegative ordered.

According to Lemma \[lemma:aNonnegative\], all the elements in ${\a^\star}$ are nonnegative. Suppose ${\a^\star}$ is not nonnegative ordered, then there must exist $i, j $ ($1 \leq i < j \leq L $) such that ${\a^\star}(i) > {\a^\star}(j) \geq 0$. According to Lemma \[lemma:a0\], ${\a^\star}(i) > 0$ implies $\h(i) > 0$. According to Lemma \[lemma:hijEqual\] and Remark \[remark:hijEqual\], ${\a^\star}(i) > {\a^\star}(j)$ implies $\h(i) \neq \h(j)$ and thus $\h(j) > \h(i) > 0$. Then, $\h(i) {\a^\star}(j) + \h(j) {\a^\star}(i) > \h(i) {\a^\star}(i) + \h(j) {\a^\star}(j)$.

Define $\a'$ as: $\a'(i) = {\a^\star}(j)$, $\a'(j) = {\a^\star}(i)$, and $\a'(\ell) = {\a^\star}(\ell)$, $\forall \ell \notin \{i,j\}$. Obviously, $\norm{\a'} = \norm{{\a^\star}}$, and $\h^T \a'
= \sum_{\ell = 1}^L \h(\ell) \a'(\ell)
= \sum_{\ell \notin \{i,j\}} \h(\ell) \a'(\ell) + \h(i) \a'(i) + \h(j) \a'(j)
= \sum_{\ell \notin \{i,j\}} \h(\ell) {\a^\star}(\ell) + \h(i) {\a^\star}(j) + \h(j) {\a^\star}(i)
> \sum_{\ell \notin \{i,j\}} \h(\ell) {\a^\star}(\ell) + \h(i) {\a^\star}(i) + \h(j) {\a^\star}(j)
= \sum_{\ell = 1}^L \h(\ell) {\a^\star}(\ell)
= \h^T {\a^\star}\geq 0$. Then according to , $\mathcal{R} \left( \h, \a' \right) > \mathcal{R} \left( \h, {\a^\star}\right)$, which implies ${\a^\star}$ is not optimal. Therefore, ${\a^\star}$ must be nonnegative ordered.

Relaxation to QPs
-----------------

As stated before, it suffices to obtain the optimal coefficient vector for a nonnegative ordered channel vector. Thus, in the following, we will focus on solving the problem in for a nonnegative ordered channel vector ${\bar{\h}}$. We first relax this problem to a series of QPs.

Denote the optimal coefficient vector for ${\bar{\h}}$ as ${{\bar{\a}}^\star}$. According to Theorem \[theorem:aNonnegativeOrdered\], the maximum element in ${{\bar{\a}}^\star}$ is ${{\bar{\a}}^\star}(L)$. Suppose ${{\bar{\a}}^\star}(L)$ is known to be $\bar{a}^\star_L \in \Zbb \backslash \{0\}$, then the problem in can be relaxed as a QP, $$\begin{aligned}
\label{equation:RelaxedQP}
&\underset{\a}{\text{minimize}}
&&\a^T \bsG \a\\
&\text{subject to}
&&\a \in \Rbb^L\\
&
&&\a(L) = \bar{a}_L^\star
\end{aligned}$$ where $\bsG$ is as defined in . The problem is convex since $\bsG$ is positive definite. Denote the solution of this relaxed problem as ${{\bar{\a}}^\dagger}\in \Rbb^L$. The intuition behind this relaxation is that appropriate quantization of the real-valued optimal ${{\bar{\a}}^\dagger}$ with the constraint $\a(L) = \bar{a}_L^\star$ will lead to the integer-valued optimal ${{\bar{\a}}^\star}$ or at least a close-to-optimal one with a high probability.

However, since ${{\bar{\a}}^\star}(L)$ is unknown, we alternatively approximate the problem in by solving a series of QPs, i.e., solving the following QP multiple times for $k = 1,2,\cdots,K$. $$\begin{aligned}
\label{equation:RelaxedQPs}
&\underset{\a}{\text{minimize}}
&&\a^T \bsG \a\\
&\text{subject to}
&&\a \in \Rbb^L\\
&
&&\a(L) = k
\end{aligned}$$ Denote the solution to the above QP with the constraint $\a(L) = k$ as ${{\bar{\a}}^\dagger}_k$. For simplicity, we use $\{s_k\}$ to denote the set with elements being $s_k, k=1,2,\cdots,K$ in the following. As long as $K \geq {{\bar{\a}}^\star}(L)$, the solution ${{\bar{\a}}^\dagger}$ to the QP in will be included in the set of solutions $\{ {{\bar{\a}}^\dagger}_k \}$ to the QPs in . Fortunately, to obtain the solution set $\{ {{\bar{\a}}^\dagger}_k \}$, it is sufficient to solve merely one QP in with $k = 1$, according to the following theorem.

\[theorem:aLinear\] Denote the solution to the QP in with the constraint $\a(L) = k$ as ${{\bar{\a}}^\dagger}_k$, then ${{\bar{\a}}^\dagger}_k = k {{\bar{\a}}^\dagger}_1$.

$$\begin{aligned}
{{\bar{\a}}^\dagger}_k
&= \arg \min_{ \a \in \Rbb^L, \a(L) = k }
\a^T \bsG \a\\
&= k \left( \arg \min_{ \mathbf{t} \in \Rbb^L, \mathbf{t}(L) = 1 }
(k \mathbf{t}^T) \bsG (k \mathbf{t}) \right) \quad ({\rm Let\ } \a = k \mathbf{t}.)\\
&= k \left( \arg \min_{ \mathbf{t} \in \Rbb^L, \mathbf{t}(L) = 1 }
k^2 \mathbf{t}^T \bsG \mathbf{t} \right)\\
&= k \left( \arg \min_{ \mathbf{t} \in \Rbb^L, \mathbf{t}(L) = 1 }
\mathbf{t}^T \bsG \mathbf{t} \right)\\
&= k {{\bar{\a}}^\dagger}_1. \qedhere\end{aligned}$$

The closed-form expression of ${{\bar{\a}}^\dagger}_1$ can be readily obtained by solving a linear system as stated in the following theorem, which is the key for the low complexity of our method.

\[theorem:aBarDagger1\] Let ${{\bar{\a}}^\dagger}_1$ be the optimal solution to the QP in with the constraint $\a(L) = 1$, then $${{\bar{\a}}^\dagger}_1 =
\begin{bmatrix}
\rr\\
1
\end{bmatrix},$$ where $$\label{equation:r}
\rr = - \Bigl( \bsG(1\!:\!L-1,1\!:\!L-1) \Bigr)^{-1}
\bsG(1\!:\!L-1,L).$$

The QP in has only an equality constraint, and thus is linear and particularly simple [@Luenberger2008]. We now derive the closed-form solution with the Lagrange multiplier method. Let the Lagrange multiplier associated with the constraint $\a(L) = 1$ be $\lambda \geq 0$, then the Lagrangian is $$\mathcal{L}(\a,\lambda) = \a^T \bsG \a + \lambda \left( \a(L) - 1 \right).$$ The optimal solution can be obtained by letting the derivative of the Lagrangian be zero, i.e., $$\begin{aligned}
\frac{\partial}{\partial \a}\mathcal{L}(\a,\lambda)
= (\bsG + \bsG^T) \a +
\begin{bmatrix}
\0\\
\lambda
\end{bmatrix}
= 2 \bsG \a +
\begin{bmatrix}
\0\\
\lambda
\end{bmatrix}
= \0.\end{aligned}$$ Let $\rr = \a(1\!:\!L-1)$, $\lambda = 2 \mu$, and write $\bsG$ and $\a$ as block matrices, then $$\begin{aligned}
\left[
\arraycolsep=1pt\def\arraystretch{1.5}
\begin{array}{c|c}
\bsG(1\!:\!L-1,1\!:\!L-1)
&\bsG(1\!:\!L-1,L)\\
\hline
\bsG(L,1\!:\!L-1)
&\bsG(L,L)
\end{array}
\right]
\begin{bmatrix}
\rr\\
1
\end{bmatrix}
+
\begin{bmatrix}
\0\\
\mu
\end{bmatrix}
= \0.
\end{aligned}$$ In the above equation, observe that $$\begin{aligned}
\left[
\arraycolsep=1pt\def\arraystretch{1.5}
\begin{array}{c|c}
\bsG(1\!:\!L-1,1\!:\!L-1)
&\bsG(1\!:\!L-1,L)
\end{array}
\right]
\begin{bmatrix}
\rr\\
1
\end{bmatrix}
= \0,
\end{aligned}$$ then $$\rr = - \Bigl( \bsG(1\!:\!L-1,1\!:\!L-1) \Bigr)^{-1}
\bsG(1\!:\!L-1,L),$$ and the results follow immediately.

Calculating ${{\bar{\a}}^\dagger}_1$ in Theorem \[theorem:aBarDagger1\] has a complexity order of $O(L^3)$ due to the matrix inversion in the expression of $\bsr$. We note that this complexity order can be reduced to $O(L)$ by the following lemma.

\[lemma:r\] Equation  can be expressed in a simpler form as $$\bsr = \frac{\u(L)}{1-\norm{\u(1\!:\!L-1)}^2}\u(1\!:\!L-1),$$ where the “normalized" channel vector $\u$ is defined as $$\label{equation:u}
\u = \sqrt{\frac{P}{1+P\norm{{\bar{\h}}}^2}}{\bar{\h}}.$$

Express $\bsG$ in and $\bsr$ in in terms of $\u$, $$\begin{gathered}
\bsG = \I - \u\u^T,\\
\bsr = \big(\I^{L-1} - \u(1\!:\!L-1)\u^T(1\!:\!L-1)\big)^{-1}\u(L)\u(1\!:\!L-1),\end{gathered}$$ where $\I^{L-1}$ denotes the identity matrix with dimension $L-1$. Then it is easy to verify that Lemma \[lemma:r\] holds.

With Theorems \[theorem:aLinear\] and \[theorem:aBarDagger1\], the $K$ solutions $\{ {{\bar{\a}}^\dagger}_k \}$ to the $K$ QPs in can be easily obtained.

The next step is to quantize the real-valued approximations $\{ {{\bar{\a}}^\dagger}_k \}$ to integer vectors by applying the floor or the ceiling functions to each of the elements. One issue that still remains is how to determine the value of $K$. Intuitively, the larger $K$, the better. Actually, it is sufficient to set $K$ as $$\begin{aligned}
\label{equation:K}
K = \arg\max_{\norm{\floor{{{\bar{\a}}^\dagger}_k}}^2 < 1 + P \norm{{\bar{\h}}}^2} k
\end{aligned}$$ according to the following lemma from [@Nazer2011].

\[lemma:aRange\] For a given channel vector $\h$, the computation rate $\mathcal{R} \left( \h, \a \right)$ is zero if the coefficient vector $\a$ satisfies $$\begin{aligned}
\norm{\a}^2 \geq 1 + P \norm{\h}^2.\end{aligned}$$

\[remark:KPractical\] For high SNR (i.e., large $P$) and large dimensions of ${\bar{\h}}$, $K$ in can be quite huge. However, as we will show in the next section, for i.i.d. Gaussian channel entries with high SNR, $K$ can be set to a rather small value without degrading the average computation rate.

In practice, for i.i.d. Gaussian channel entries, we set an upper bound for $K$ as $K_u$, which is determined off-line according to the simulation results, such that the simulated average computation rate at 20dB with $K$ being $K_u$ is greater than 99% of that with $K$ being $K_u+1$. We set $K_u$ based on rates at 20dB since the value of $K$ influences more the rates at larger SNR, and 20dB is the the maximum SNR considered in this paper. Then, we set $K$ as the maximum integer that is no greater than $K_u$ while satisfies at the same time, i.e., $$\begin{aligned}
\label{equation:KPractical}
K = \arg\max_{\substack{\norm{\floor{{{\bar{\a}}^\dagger}_k}}^2 < 1 + P \norm{{\bar{\h}}}^2\\
k\leq K_u}} k.
\end{aligned}$$ For implementation, $K$ can be easily determined by using a bi-section search.

Quantization
------------

We propose the *successive quantization* algorithm shown in Algorithm \[agorithm:SuccessiveQuantization\] to quantize the $K$ real-valued approximations $\{ {{\bar{\a}}^\dagger}_k \}$ to integer-valued vectors $\{ {{\bar{\a}}^\diamond}_k \}$ that serve as candidates of a suboptimal coefficient vector ${{\bar{\a}}^\diamond}$. For convenience, define $$\begin{aligned}
\label{equation:f}
f(\w) \triangleq \w^T \bsG \w,\end{aligned}$$ where $\w\in\Rbb^L$, and $\bsG$ is defined in with $\h$ being nonnegative ordered. Also, let $\floor{\w}_\ell$ and $\ceil{\w}_\ell$ be the vectors generated from $\w$ by applying the floor and the ceiling operations on the $\ell$-th element only, respectively.

${{\bar{\a}}^\diamond}_k\leftarrow {{\bar{\a}}^\dagger}_k$

To simplify the inequality condition $f\left(\floor{{{\bar{\a}}^\dagger}_k}_\ell\right) < f\left(\ceil{{{\bar{\a}}^\dagger}_k}_\ell\right)$ in line \[line:FloorCondition\] of Algorithm \[agorithm:SuccessiveQuantization\], we first introduce the following lemma.

\[lemma:FloorCondition\] For the function $f(\w) = \w^T \bsG \w$ defined in where $\bsG^T = \bsG$, the inequality condition $f\left(\floor{\w}_\ell\right) < f\left(\ceil{\w}_\ell\right)$ is equivalent to $$\begin{aligned}
\label{equation:FloorCondition}
2 \left(\floor{\w}_\ell\right)^T \bsG(:,\ell) + \bsG(\ell,\ell) > 0.\end{aligned}$$

$f\left(\floor{\w}_\ell\right) < f\left(\ceil{\w}_\ell\right)$ implies $\floor{\w}_\ell \neq \ceil{\w}_\ell$, i.e., $\w(\ell)$ is not an integer. Let $\e_\ell \in \Rbb^L$ be the vector with only one nonzero element $\e_\ell(\ell) = 1$, then $\ceil{\w}_\ell = \floor{\w}_\ell + \e_\ell$, and $$\begin{aligned}
&f\left(\ceil{\w}_\ell\right)
= \left(\ceil{\w}_\ell\right)^T \bsG \ceil{\w}_\ell\\
&= \left( \floor{\w}_\ell + \e_\ell \right)^T \bsG \left( \floor{\w}_\ell + \e_\ell \right)\\
&= \left(\floor{\w}_\ell\right)^T \bsG \floor{\w}_\ell
+ \left(\floor{\w}_\ell\right)^T \bsG \e_\ell
+ \e_\ell^T \bsG \floor{\w}_\ell
+ \e_\ell^T \bsG \e_\ell\\
&= f\left(\floor{\w}_\ell\right)
+ \left(\floor{\w}_\ell\right)^T \bsG(:,\ell)
+ \bsG(\ell,:) \floor{\w}_\ell
+ \bsG(\ell,\ell)\\
&= f\left(\floor{\w}_\ell\right) + 2 \left(\floor{\w}_\ell\right)^T \bsG(:,\ell) + \bsG(\ell,\ell).
\end{aligned}$$ Obviously, $f\left(\floor{\w}_\ell\right) < f\left(\ceil{\w}_\ell\right)$ is equivalent to $$\begin{aligned}
2 \left(\floor{\w}_\ell\right)^T \bsG(:,\ell) + \bsG(\ell,\ell) > 0.\tag*{\qedhere}\end{aligned}$$

\[lemma:FloorCondition2\] With Lemma \[lemma:FloorCondition\], the inequality condition $f\left(\floor{{{\bar{\a}}^\dagger}_k}_\ell\right) < f\left(\ceil{{{\bar{\a}}^\dagger}_k}_\ell\right)$ in line \[line:FloorCondition\] of Algorithm \[agorithm:SuccessiveQuantization\] can be simplified as $$\begin{aligned}
2\floor{{{\bar{\a}}^\dagger}_k(\ell)} - 2\left(\left(\floor{{{\bar{\a}}^\dagger}_k}_\ell\right)^T\u\right)\u(\ell) + 1 - \u(\ell)^2 < 0,\end{aligned}$$ where $\u$ is the normalized channel vector as defined in .

The proof is straightforward by writing $\bsG$ in terms of $\u$, and thus is omitted here.

After the quantization, a suboptimal coefficient vector ${{\bar{\a}}^\diamond}$ for ${\bar{\h}}$ is obtained with $$\begin{aligned}
\label{equation:aBarDiamond}
{{\bar{\a}}^\diamond}= \arg \min_{ \a \in \{ {{\bar{\a}}^\diamond}_k \} } \a^T \bsG \a.
\end{aligned}$$

Finally, a suboptimal coefficient vector ${\a^\diamond}$ for the original channel vector $\h$ is recovered from ${{\bar{\a}}^\diamond}$ according to Remark \[remark:Transformation\].

We summarize our proposed *QP relaxation* method in Algorithm \[agrm:QPRApproachOutline\]. The pseudocode is shown in Algorithm \[algorithm:QPRApproachCode\], where the function $[\bar{\w},\p]={\rm sort}(\w)$ sorts the elements in $\w$ in ascending order, returns the sorted vector $\bar{\w}$, and stores the original indices of the elements as vector $\p$, the function ${\rm floor}(\w)$ applies the floor operation to each element of $\w$ and returns the resulted integer vector.

1.  \[item:outline:Preprocess\] Preprocess $\h$ to the nonnegative ordered ${\bar{\h}}$ with Remark \[remark:Transformation\].

2.  \[item:outline:CalculateaBarDagger1\] Calculate ${{\bar{\a}}^\dagger}_1$, with Theorem \[theorem:aBarDagger1\] and Lemma \[lemma:r\].

3.  \[item:outline:DetermineK\] Determine $K$ with .

4.  \[item:outline:CalculateaBarDaggerk\] Calculate $\{ {{\bar{\a}}^\dagger}_k \}$, i.e., the real-valued approximations of the optimal coefficient vector for ${\bar{\h}}$, using Theorem \[theorem:aLinear\].

5.  \[item:outline:Quantize\] Quantize $\{ {{\bar{\a}}^\dagger}_k \}$ to integer-valued vectors $\{ {{\bar{\a}}^\diamond}_k \}$ with Algorithm \[agorithm:SuccessiveQuantization\].

6.  \[item:outline:Select\] Select a vector from $\{ {{\bar{\a}}^\diamond}_k \}$ to be a suboptimal coefficient vector ${{\bar{\a}}^\diamond}$ for ${\bar{\h}}$ using .

7.  \[item:outline:Recover\] Recover a suboptimal coefficient vector ${\a^\diamond}$ for $\h$ from ${{\bar{\a}}^\diamond}$ according to Remark \[remark:Transformation\].

$\bst\leftarrow$ $({\bar{\h}},\p)\leftarrow$ $b\leftarrow 1 + P\norm{\h}^2$

$\u\leftarrow (P/b)^{1/2}{\bar{\h}}$ $\bsr\leftarrow \frac{\u(L)}{1-\norm{\u(1:L-1)}^2}\u(1\!:\!L-1)$ ${{\bar{\a}}^\dagger}_1(1\!:\!L-1)\leftarrow \bsr$ ${{\bar{\a}}^\dagger}_1(L)\leftarrow 1$

()[ $K_l\leftarrow 1$ $K\leftarrow K_l$ ]{}

${{\bar{\a}}^\diamond}(1\!:\!L-1)\leftarrow \0$ ${{\bar{\a}}^\diamond}(L)\leftarrow 1$ $f_{\min}\leftarrow \norm{{{\bar{\a}}^\diamond}}^2 - \big(({{\bar{\a}}^\diamond})^T\u\big)^2$

Complexity Analysis
-------------------

Here we analyze the complexity of our algorithm, in terms of the number of flops required. Referring to the outline in Algorithm \[agrm:QPRApproachOutline\], the processing of $\h$ in step \[item:outline:Preprocess\] involves recording the signs of the elements and sorting the elements, and takes $O(L\log(L))$ flops. Calculating ${{\bar{\a}}^\dagger}_1$ in step \[item:outline:CalculateaBarDagger1\] has a complexity of $O(L)$. For the bi-section search applied to determine $K$ in step \[item:outline:DetermineK\], the maximum number of loops required to execute is $O(\log(K_u))$, the number of flops in each loop is $O(L)$, and thus the maximum cost is $O(\log(K_u)L)$. Step \[item:outline:CalculateaBarDaggerk\] takes $O(KL)$ flops. By introducing appropriate temporary variables $b$ and $d$ as shown in Algorithm \[algorithm:QPRApproachCode\], the successive quantization of a real-valued approximation ${{\bar{\a}}^\dagger}_k$ can be implemented in an efficient way in $O(L)$ flops. Thus, the complexity of quantizing all the $K$ real-valued approximations is $O(KL)$. Selecting a coefficient vector from the quantized vector set in step \[item:outline:Select\] has a cost of $O(KL)$. Step \[item:outline:Recover\] takes $O(L)$ flops. In summary, the complexity of the method is $O(L(\log(L)+\log(K_u)+K))$.

However, the above analyzed complexity expression involves the experiment-based $K_u$, and its exact order with respect to the dimension $L$ is intractable. As an alternative, we use an upper bound to approximate the cost. According to , it is easy to see that $K$ and $K_u$ are at most of order $O(\sqrt{P\norm{\h}^2})$. Then, the complexity of our method is $O(L(\log(L)+\sqrt{P\norm{\h}^2}))$. We reserve the power $P$ in the expression since we may also care about how the complexity varies when the SNR gets large.

In the complexity expression above, since the square root function is strictly concave, it follows from Jensen’s inequality that $\mathbb{E}(\sqrt{\norm{\h}^2})\leq\sqrt{\mathbb{E}(\norm{\h}^2)}$. Specifically, for i.i.d. standard Gaussian channel entries, the expectation of $\norm{\h}^2$ is $L$, and thus the corresponding average complexity of the proposed method becomes $O(L\log(L)+P^{0.5}L^{1.5})$. It is easy to see that the complexity is of order 1.5 with respect to the dimension $L$.

Extension to the Complex-Valued Channel Model
---------------------------------------------

We now consider the complex-valued channel model of the AWGN networks, and demonstrate how to apply the proposed QP relaxation method for complex-valued channels. The complex-valued channel model is defined as below.

\[definition:ComplexChannelModel\] *(Complex-Valued Channel Model)* In an AWGN network, each relay (indexed by $m=1,2,\cdots,M$) observes a noisy linear combination of the transmitted signals through the channel, $$\begin{aligned}
\label{equation:ComplexChannelModel}
\y_m = \sum_{\ell=1}^L \h_m(\ell)\x_\ell + \z_m,\end{aligned}$$ where $\x_\ell\in\Cbb^n$ with the power constraint $\frac{1}{n}\norm{\x_\ell}^2 \leq P$ is the transmitted codeword from source $\ell$ ($\ell = 1,2,\cdots,L$), $\h_m\in\Cbb^L$ is the channel vector to relay $m$, $\z_m\in\Cbb^n$ is the noise vector with entries being i.i.d. Gaussian, i.e., $\z_m\sim\mathcal{CN}\!\left(\0,\I\right)$, and $\y_m$ is the signal received at relay $m$.

Similar to what we have done for the real-valued channel model, we will focus on one relay, and ignore the subscript “$m$” for notational convenience.

Writing the summation in in the vector product form, becomes $$\begin{aligned}
\label{equation:ComplexChannelModelVectorForm}
\y = [\x_1,\x_2,\cdots,\x_L]\h + \z.\end{aligned}$$ It is well-known that a complex-valued channel model can be written in its real-valued equivalent form. Let $\Re(\w)$ denote the vector composed of the real part of $\w$, and $\Im(\w)$ denote the vector composed of the imaginary part of $\w$. The complex-valued equation  has the following real-equivalent form $$\small
\begin{aligned}
\label{equation:ComplexChannelModelRealForm}
[\Re(\y),\Im(\y)]=[\Re(\x_1),\Re(\x_2),\cdots,\Re(\x_L),\Im(\x_1),\Im(\x_2),\cdots,\Im(\x_L)]
\times
\begin{bmatrix}
\Re(\h),\Im(\h)\\
-\Im(\h),\Re(\h)
\end{bmatrix}
+ [\Re(\z),\Im(\z)].
\end{aligned}$$ It is obvious that $$\small
\begin{aligned}
\label{equation:ComplexChannelModelRealForm2}
\Re(\y)&=[\Re(\x_1),\Re(\x_2),\cdots,\Re(\x_L), \Im(\x_1),\Im(\x_2),\cdots,\Im(\x_L)]
\times
\begin{bmatrix}
\Re(\h)\\
-\Im(\h)
\end{bmatrix}
+ \Re(\z),\\
\Im(\y)&=[\Re(\x_1),\Re(\x_2),\cdots,\Re(\x_L), \Im(\x_1),\Im(\x_2),\cdots,\Im(\x_L)]
\times
\begin{bmatrix}
\Im(\h)\\
\Re(\h)
\end{bmatrix}
+ \Im(\z).
\end{aligned}$$ Then we can view $\begin{bmatrix}
\Re(\h)\\
-\Im(\h)
\end{bmatrix}$ and $\begin{bmatrix}
\Im(\h)\\
\Re(\h)
\end{bmatrix}$ as two $2L$-dimensional real-valued channels, and view $\Re(\x_\ell)$ and $\Im(\x_\ell)$ as two independent $n$-dimensional real-valued transmitted codewords. Assume equal power allocation on the real part and the imaginary part of each transmitted codeword, i.e., $\norm{\Re(\x_\ell)}^2=\norm{\Im(\x_\ell)}^2$, then the power constraint of each real-valued transmitted codeword is $\frac{1}{n}\norm{\Re(\x_\ell)}^2\leq\frac{1}{2}P$ and $\frac{1}{n}\norm{\Im(\x_\ell)}^2\leq\frac{1}{2}P$.

Based on the above interpretation, we can apply the proposed QP relaxation method to each of the two $2L$-dimensional channels to find the corresponding coefficient vectors. It should be noted that we only need to find the coefficient vector for one of the $2L$-dimensional channels, which saves half of the computation cost. Let $\a$ be a Gaussian integer, and assume the found coefficient vector for $\begin{bmatrix}
\Re(\h)\\
-\Im(\h)
\end{bmatrix}$ is $\begin{bmatrix}
\Re(\a)\\
-\Im(\a)
\end{bmatrix}$. Then, according to Theorem \[theorem:ProblemTransformation\], the coefficient vector for $\begin{bmatrix}
\Im(\h)\\
\Re(\h)
\end{bmatrix}$ is $\begin{bmatrix}
\Im(\a)\\
\Re(\a)
\end{bmatrix}$. In this sense, for each complex-valued channel vector $\h$, we can find a Gaussian integer as the best coefficient vector $\a$ using the QP relaxation method.

Numerical Results {#section:NumericalResults}
=================

![Average computation rate for $L=4$ using our QPR method with different $K$.[]{data-label="figure:K_L4"}](K_L4-eps-converted-to.pdf "fig:"){width="240pt"}\

In this section, we present some numerical results to demonstrate the effectiveness and efficiency of our QP relaxation approach. As explained before, finding the coefficient vector for a complex-valued channel can be transformed to finding the coefficient vector for a real-valued channel. Thus, we focus on the real-valued channels here. We consider the case where the entries of the channel vector $\h$ are i.i.d. standard Gaussian, i.e., $\h\sim\mathcal{N}(\0,\I)$. In our simulations, the dimension $L$ ranges from 2 to 16, and the power $P$ ranges from 0dB to 20dB. For a given dimension and a given power, we randomly generate $10000$ instances of the channel vector, and apply the QP relaxation method to find the coefficient vectors, and calculate the corresponding average computation rate.

We first show that as stated in Remark \[remark:KPractical\], for high dimension and large power, the number of real-valued approximations $K$ can be set to a rather small value without degrading the rate apparently. As shown in Figure \[figure:K\_L4\], for dimension $L=4$ and power $P$ from 0dB to 20dB, the average computation rate quickly converges as $K$ increases from 1 to 4. Further increasing $K$ up to 10 incurs additional computational cost with little improvement in the average computation rate.

With the above observation, it is reasonable to introduce the upper bound $K_u$ for $K$, and adopt the criterion in to determine $K$. $K_u$ can be calculated off-line by simulations prior to applying the method, which incurs no additional processing complexity in real-time. The values of $K_u$ according to the simulation results are listed in Table \[table:Ku\].

[c|\*[15]{}[r]{}]{} $L$ &2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12 & 13 & 14 & 15 & 16\
$K_u$ &2 & 3 & 4 & 5 & 5 & 5 & 6 & 6 & 6 & 6 & 7 & 6 & 6 & 6 & 4\

\

We then show the effectiveness of our method by comparing the average computation rate with those of other existing methods. The methods covered include the following.

-   Our QP relaxation (QPR) method that gives the suboptimal solution.

-   The branch-and-bound (BnB) method proposed by Richter [*et al.*]{} in [@Richter2012] that provides the optimal solution.

-   The method developed by Sahraei and Gastpar in [@Sahraei2014] that finds the optimal solution with an average-case complexity of $O(P^{0.5}L^{2.5})$ for i.i.d. Gaussian channel entries. We refer to this method as the “SG" method for short.

-   The LLL method proposed by Sakzad [*et al.*]{} in [@Sakzad2012], which is based on the LLL lattice reduction (LR) algorithm. The parameter $\delta$ in the LLL LR algorithm is set as 0.75 since further increasing $\delta$ towards 1 achieves little gain in the computation rate but requires more computation labor. Although the LLL LR algorithm has known average complexity for some cases [@Daude1994; @Ling2007], its average complexity for our case is unknown, and the worst-case complexity could be unbounded [@Jalden2008].

-   The quantized search (QS) method developed by Sakzad [*et al.*]{} in [@Sakzad2012]. The search consists of two phases: 1) an integer $\alpha_0$ between 1 and $\floor{P^{1/2}}$ that provides the maximum rate is selected as the initial value of the amplifying factor $\alpha$; 2) the amplifying factor is then refined by searching in $[\alpha_0-1,\alpha_0+1]$ with a step size 0.1. After the amplifying factor $\alpha$ is determined, the coefficient vector $\a$ is set as $\round{\alpha\h}$. An improved version of the QS method is the quantized exhaustive search (QES) method proposed in [@Sakzad2014], which was developed for complex-valued channels.

-   The rounding method that simply sets the coefficient vector by rounding the channel vector to an integer-valued vector.

As shown in Figure \[figure:Rate\], the optimal methods, i.e., the BnB method and the SG method, always provide the highest average computation rates for all dimensions and over the whole SNR regime, as expected. The LLL method provides close-to-optimal average computation rates. Our proposed QPR method also offers close-to-optimal average computation rates for almost all the dimensions and SNR values considered, except that its performance degrades a little bit for high dimensions at high SNR as shown in Figure \[figure:Rate\_L16\]. The performance of our QPR method improves slightly compared with the version we presented in [@Zhou2014; @Wen2015]. The reason is that here we initialize the output coefficient vector as $[0,\cdots,0,1]^T$, which definitely results non-zero computation rate, while in the previous version the output coefficient vector could result zero computation rate.

\
\

Finally, we demonstrate the efficiency of the proposed QPR method by comparing the running time of finding the coefficient vectors for 10000 channel vector samples. The methods considered include those that provide optimal rates and close-to-optimal rates, i.e., the SG method, the BnB method, the LLL method, and our QPR method. The running time varies for different SNR values, and thus we compare the running time with $P$ being 0dB, 10dB and 20dB. As shown in Figure \[figure:Time\], the proposed QPR method is much more efficient than all the other methods, especially for high dimensions. Specifically, the running time of the optimal methods can be one scale larger than that of the QPR method.

In summary, for i.i.d. Gaussian channel entries, our proposed QPR method offers close-to-optimal average computation rates with a much lower complexity than that of the existing optimal and close-to-optimal methods.

Conclusions {#section:Conclusions}
===========

In this paper, we considered the compute-and-forward network coding design problem of finding the optimal coefficient vector that maximizes the computation rate at a relay, and developed the quadratic programming (QP) relaxation method that finds a high quality suboptimal solution. We first revealed some useful properties of the problem, and relaxed the problem to a series of QPs. We then derived the closed-form solutions of the QPs, which is the key to the efficiency of our method, and proposed a successive quantization algorithm to quantize the real-valued solutions to integer vectors that serve as candidates of the coefficient vector. Finally, the candidate that maximizes the computation rate is selected as the best coefficient vector. For $L$-dimensional channel vectors with i.i.d. Gaussian entries, the average-case complexity of the proposed QP relaxation method is of order 1.5 with respect to the dimension $L$. Numerical results demonstrated that our QP relaxation method offers close-to-optimal computation rates, and is much more computationally efficient than the existing methods that provide the optimal computation rates as well as the LLL method that also provides close-to-optimal computation rates.

[^1]: Baojian Zhou and Wai Ho Mow are with the Department of Electronic and Computer Engineering, The Hong Kong University of Science and Technology, Clear Water Bay, NT, Hong Kong (e-mail: {bzhouab, eewhmow}@ust.hk). Baojian Zhou was supported by a grant from University Grants Committee of the Hong Kong Special Administrative Region, China (Project No. AoE/E-02/08).

[^2]: Jinming Wen is with the Laboratoire de l’Informatique du Parallélisme, (CNRS, ENS de Lyon, Inria, UCBL), Université de Lyon, Lyon 69007, France (e-mail: jwen@math.mcgill.ca). Jinming Wen was supported in part by ANR through the HPAC project under Grant ANR 11 BS02 013.
