# Using Hierarchical Parallelism to Accelerate the Solution of Many Small
Partial Differential Equations

Jacob Merson  Department of Mechanical,
Aerospace
and Nuclear Engineering
Rensselaer Polytechnic Institute
110 Eighth St. Troy, NY 12180
Email<EMAIL_ADDRESS>Mark S. Shephard  Scientific Computation
Research Center
Rensselaer Polytechnic Institute
110 Eighth St. Troy, NY 12180
Email<EMAIL_ADDRESS>
###### Abstract

This paper presents efforts to improve the hierarchical parallelism of a two
scale simulation code. Two methods to improve the GPU parallel performance
were developed and compared. The first used the NVIDIA Multi-Process Service
and the second moved the entire sub-problem loop into a single kernel using
Kokkos hierarchical parallelism and a PackedView data structure. Both
approaches improved parallel performance with the second method providing the
greatest improvements.

##  1 Introduction

Hierarchical multiscale methods are commonly used to model engineering
materials which exhibit complex micromechanical behavior that is not easily
captured with standard constitutive modeling [1, 2, 3, 4]. This behavior is
often caused by a material having a makeup of discrete constituents such as
atoms, molecules, fibers, etc. For a two scale analysis the macroscale, or
engineering scale, partial differential equations are often discretized by
finite element methods. At each material point in the macroscale, a microscale
sub-problem made from discrete components is solved to obtain the material
constitutive properties at that point. This information passing scheme is
often referred to as the FE2 scheme drawing from the fact that there are is a
finite element analysis occurring on two scales [1].

In our multiscale implementation we use parallelism at multiple levels. On the
macroscale, we use domain level parallelism to break up our finite element
mesh. This is implemented using the SCOREC parallel unstructured meshing
infrastructure (PUMI) [5]. Each microscale sub-problem is independent which
leads to an “embarrassingly parallel” algorithm. To couple the scales, we use
the Adaptive Multiscale Simulation Infrastructure which breaks processors on
the target platform into an independent processor set for each scale [6]. In
our previous work, the individual microscale sub-problems were not
parallelized. Due to our increased understanding of the sub-scale physics in
our problem of interest, we have increased the number of degrees of freedom in
the microscale sub-problems by two to three orders of magnitude [7, 8]. This
increase in problem size makes parallelization of the individual sub-problems
essential to performing analyses with physical relevance.

Due to the increased computational cost associated with changes in the
solution method described in section 2 and the increase in microscale problem
size, the solution to the microscale problems became a performance bottleneck.
Therefore, we ported our code to use GPU parallelism for the microscale
problems. For our problem of interest, the microscale problems often have less
than 100,000 degrees of freedom. This is not large enough to saturate the GPU
with our current analysis methods which primarily consists of vector
operations—similar to BLAS level 1 (See figure 1).

To achieve adequate GPU throughput, multiple microscale problems must be
solved on each GPU at a time. This was accomplished using two methods. The
first was using NVIDIA Multi-Process Service (MPS) which unintrusively allows
multiple processes to launch GPU kernels at a time [9]. Although MPS gave a
good speedup, our analysis was still limited by kernel launch overhead and by
the limited number of processors on each node to launch GPU kernels from.
Additionally, great care must be taken to use MPS in an environment with
multiple GPUs per node because improper MPS setup causes a drastic reduction
in the weak scalability. The second method for increasing GPU throughput was
to pack multiple microscale problems into a single kernel launch using the
Kokkos hierarchical parallelism construct [10]. A more thorough discussion of
this method is given in section 3. To aid in understanding the selection of
our parallelization strategy, a more comprehensive discussion of multiscale
modeling techniques is presented below.

##  2 Multiscale Modeling of Biological Tissues

One particular application of these methods is to model biological tissues
which are made of constituent collagen fiber networks. Typically, modeling
biological tissues requires large strain analysis because they are soft and
physiological strains can often exceed 50%. The FE2 method has been extended
to allow for large strains; however, the methods discussed in the literature
utilize implicit finite element methods for both analysis scales [11, 2].
Unfortunately, the deformation of fiber networks is highly nonlinear and the
network can go through bifurcation points, or may not be isostatic, i.e. the
tangent stiffness matrix can be singular during an analysis [12]. As a result,
athermal fiber networks, such as collagen networks, are typically modeled with
explicit finite element methods [13, 14]. The use of a purely explicit
analysis for the sub-scale problems in a multiscale analysis is problematic
because kinetic energy is lost in the microscale-to-macroscale coupling.
Additionally, inertial effects can change the microscale material properties.
To work around these issues, a dynamic relaxation method is used for the
microscale problems. Dynamic relaxation works by mapping a static analysis to
a damped dynamic explicit one where the system residual is monitored for
convergence [15].

The use of the dynamic relaxation method at the microscale greatly improves
the global strains which the multiscale method can achieve for fibrous
materials, however it imposes a significant computational cost. Our previous
studies have shown the multiscale analysis of biological tissues at scale
using homogeneous computing technologies with MPI based parallelism [6, 16].
Due to the increased computational cost of dynamic relaxation, physiologically
relevant problems are no longer accessible through MPI based parallelism
alone. Therefore, the microscale portions of our analysis code have been
ported to use GPU accelerators.

The variability in GPU programming environments across hardware vendors poses
a significant challenge to the maintainability of a GPU accelerated code which
must run on a variety of systems. Therefore, we chose to use the Kokkos C++
library for GPU support. Kokkos is a C++ programming model which is designed
to enable performance portability [10, 17]. The Kokkos team has committed to
maintaining support for all of the vendors who are providing accelerators for
the Department of Energy leadership class computing resources (AMD, Intel,
NVIDIA). This support allows for writing a single version of the analysis code
that will run across most of the easily accessible GPU accelerators.

##  3 Parallel Implementation

Figure 1: Dynamic relaxation algorithm

1:Load mesh and compute edge connectivity

2:Compute the mesh connectivity array

3:Transfer the connectivity array to the device (GPU)

4:Set displacement boundary conditions on fixed nodes (fixed dof vector
operation)

5:Compute mass matrix (finite element integration)

6:getInternalForces($\mathbf{u}$)

7:update accelerations (free dof vector operation)

8:repeat

9: Compute next time step

10: Partial velocity update (free dof vector operation)

11: Update displacements (free dof vector operation)

12: getInternalForces($\mathbf{u}$)

13: Compute damping force (free dof vector operation)

14: Compute force residual (free dof vector reduction)

15: Update Accelerations (free dof vector operation)

16: Partial velocity update (free dof vector operation)

17: Optionally Check Energy Balance (3 vector reductions)

18: Update iteration count

19:until Force residual converged

Figure 1 gives the basic dynamic relaxation algorithm we used for the
microscale sub-problems. This algorithm is identical to a two step central
difference method found in any finite element text book, with the exception
that the convergence criteria is based on a force residual measurement rather
than time. Note that each sub-scale problem converges at a different rate
which can lead to load imbalance.

In the naive approach, this algorithm was carried out using fused kernels for
any subsequent operations with the same loop characteristics. The benefits of
kernel fusion have been discussed extensively in the literature both for the
case of explicit ODEs, and general GPU computations [18, 19, 20]. The
GetInternalForces subroutine accounts for two Kernel launches: the first to
zero the internal force vector, and the second to scatter the elemental
internal forces to the nodes. The current implementation uses atomic
operations to scatter the forces.

In this naive approach, a number of microscale sub-problems were assigned to
each MPI rank, and were executed serially with respect to each other within
each rank. Despite the use of GPU acceleration for the vector operations, this
approach had poor performance for the sub-scale problems with small numbers of
degrees of freedom when compared with a CPU-only implementation with serial
vector operations. To unintrusively improve this naive approach, NVIDIA MPS
was used to allow kernels from multiple MPI ranks to run concurrently. The use
of MPS led to significant performance improvements for small DOF problems
compared with the naive case. Problem size and number of simultaneous MPI
ranks used with MPS can have a drastic effect on performance. All MPS results
presented in section 4 use 32 MPI ranks per GPU which gives the best
performance in the range of problem sizes discussed here.

Since the loop in algorithm 1 executes millions of times per macroscale
simulation step, we observed that this approach had significant kernel launch
overhead. To overcome this, we moved the entire loop into a single kernel.
This was done using Kokkos hierarchical parallelism which uses teams of
threads to enable a 2D map to the hardware. The CUDA reciprocal to this
mechanism is launching a 1D grid of 1D blocks. Since our sub-scale problems
each have less than 10,000 free degrees of freedom, we found that good
performance could be achieved by assigning one thread team to each sub-scale
problem. Here, we juxtapose the free degrees of freedom which are those
without any Dirichlet constraints, to what we call degrees of freedom which
are all potential degrees of freedom. Unlike an implicit FEM method, the
constrained degrees of freedom cannot be completely eliminated as they are
needed for the internal force computation. Reordering the fixed degrees of
freedom to a contiguous block at the end of the displacement array allows most
of the update algorithm to only operate on the smaller proportion of free
degrees of freedom (figure 1).

The choice of number of threads per team had a strong effect on performance.
The ideal number of threads per team is a function of the microscale problem
size. All presented results use 512 threads per team, which provided a good
compromise for the performance of the smallest and largest systems we tested.

A PackedView data structure which has similar semantics to a Kokkos DualView
was used to allow effective access to N-D vector data within each thread team
[21]. This data structure uses a row vector and value vector, similar to those
from compressed row storage (CRS), to store the data associated with all sub-
scale problems on the current MPI rank in a contiguous array in memory. Each
sub-scale problem gains access to the correct portion of memory through a
Kokkos Subview. In some ways, this structure is similar to a Kokkos View of
Views. However, with the current implementation the PackedView can not be
resized after initialization. A comprehensive performance comparison between
the PackedView data structure, and View of Views has not been performed to
date. This differs from the StaticCrsGraph in Kokkos which cannot handle non-
integral datatypes, and does not have DualView semantics.

Although moving the analysis loop inside of a single kernel launch was
effective for our problems of interest, it can easily succumb to low
performance from high register pressure. Significant effort had to be made to
reduce the register pressure and ensure that multiple warps could be
concurrently scheduled. One mechanism we used to reduce register pressure was
to move some of the variables which are carried across loop iterations such as
the pseudo-time and the loop iteration count into shared memory. We found that
performance gain from the reduction in register pressure outweighed the loss
in bandwidth from moving these variables to shared memory. The need to reduce
register usage in this single kernel implementation led us to favor a
stripped-down version of our algorithm which was specific to the physical
system at hand. In other words, flexibility of our code had to be sacrificed
to obtain improved performance characteristics.

##  4 Results

Figure 2: Speedup of the naive loop based analysis normalized by the number of
concurrent sub-problems compared with the loop based single sub-problem. A
flat line corresponds to a linear increase in runtime. Each data point is the
mean of three analysis runs. Figure 3: Speedup of the thread team based
analysis normalized by the number of concurrent sub-problems compared with the
thread team based single sub-problem. Each data point is the mean of three
analysis runs.

The performance results presented here, are all computed on a single Volta
V100 GPU—part of an IBM AC922 node. Each AC922 node contains 2, 20 core IBM
power 9 processors clocked at 3.15GHz, 512 GiB of RAM, and 6 Volta V100 GPUs.
The code is compiled with version 16.1.0 of IBM’s XL compiler for host code,
version 10.1 of Cuda, and version 3.1 of Kokkos. The MPS results make use of
Spectrum MPI version 10.3.

Figures 2 and 3 show the runtime of a single sub-problem divided by runtime
normalized by the number of concurrent sub-problems. This gives a measure of
the speedup of a single sub-problem when computed in a concurrent batch. Since
we are using the analysis technique’s own single sub-problem runtime as a
baseline for the speedup, we call this the “self speedup”. For the naive loop
based case (figure 2), we see the self speedup is very flat which indicates
the expected linear increase in runtime. The smallest problem size sees a
slight self speedup. When thread team based parallelization is used, a
significant self speedup is observed (figure 3). Here we see a initial regime
of linear self speedup and a plateau regime for large numbers of concurrent
sub-problems. In this initial linear scaling regime, the runtime remains flat
since the numerical workload is not large enough to overcome the kernel launch
latency. Interestingly, the initial self speedup is almost identical for each
of the problem sizes we tried. The plateau region show that as the number of
degrees of freedom in the problem increase, the self speedup decreases.

Figure 4: Speedup of the team based (solid line) and MPS based (dashed line)
analysis over the naive approach. The MPS lines correspond to launching GPU
kernels from 32 MPI ranks simultaneously. Each data point is the mean of three
analysis runs.

Figure 4 shows the speedup of the thread team based and MPS based analysis
methods over the naive loop based approach. In this plot, we see an initial
linear scaling and a plateau region. We observe that as the problem size
increases, the speedup obtained from the team thread based method decreases.
This is likely due to a reduction in percentage of the problem which resides
in the cache. We also observe that MPS based parallelism does provide some
speedup over the naive approach, but it is not as effective as the thread team
based approach. Also, for MPS, the speedup plateau does not depend on the
problem size. One way to interpret these results is that for maximum
efficiency, at least 80 concurrent sub-problems should be run on each GPU.
Since the V100 has 80 SMs (streaming multiprocessors), this is consistent with
each Kokkos team (CUDA block) occupying a single SM.

The speedups achieved using thread team parallelism tell a compelling story
that moving an analysis loop inside of a single heavy weight kernel can be an
effective optimization mechanism for problems that need to solve many problems
which cannot saturate the GPU on their own. Although MPS seemed like it might
be a reasonable solution, it suffered from still incurring a high kernel call
latency due to the many kernels which were being called inside of a hot loop.
Additionally, the MPS solution was not able to make as effective use of the
cache since many different sub-scale problems were competing to be scheduled
simultaneously, and each sub-scale problem that was scheduled in an
interleaved fashion would cause cache misses.

## Acknowledgments

This work was supported in part by the National Institutes of Health (NIH)
through Grant No. U01 AT010326-06. Also, this material is based upon work
supported by the National Science Foundation Graduate Research Fellowship
under Grant No. DGE-1744655.

## References

  * [1] F. Feyel and J.-L. Chaboche, “FE2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 183, no. 3, pp. 309–330, Mar. 2000. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0045782599002248
  * [2] R. Smit, W. Brekelmans, and H. Meijer, “Prediction of the mechanical behavior of nonlinear heterogeneous systems by multi-level finite element modeling,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 155, no. 1-2, pp. 181–192, Mar. 1998\. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S0045782597001394
  * [3] C. Miehe, “Computational micro-to-macro transitions for discretized micro-structures of heterogeneous materials at finite strains based on the minimization of averaged incremental energy q,” _Comput. Methods Appl. Mech. Engrg._ , p. 33, 2003\.
  * [4] P. Kanouté, D. P. Boso, J. L. Chaboche, and B. A. Schrefler, “Multiscale Methods for Composites: A Review,” _Archives of Computational Methods in Engineering_ , vol. 16, no. 1, pp. 31–75, Mar. 2009. [Online]. Available: http://link.springer.com/10.1007/s11831-008-9028-8
  * [5] D. A. Ibanez, E. S. Seol, C. W. Smith, and M. S. Shephard, “PUMI: Parallel Unstructured Mesh Infrastructure,” _ACM Transactions on Mathematical Software_ , vol. 42, no. 3, pp. 1–28, May 2016. [Online]. Available: http://dl.acm.org/citation.cfm?doid=2935754.2814935
  * [6] W. R. Tobin, “The Adaptive Multiscale Simulation Infrastructure,” Ph.D. dissertation, Rensselaer Polytechnic Institute, Troy, Ny, Jul. 2018.
  * [7] A. S. Shahsavari and R. C. Picu, “Size effect on mechanical behavior of random fiber networks,” _International Journal of Solids and Structures_ , vol. 50, no. 20, pp. 3332–3338, Oct. 2013. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0020768313002382
  * [8] J. Merson and R. Picu, “Size Effects in Random Fiber Networks Controlled by the Use of Generalized Boundary Conditions,” _Under Review_ , 2020.
  * [9] NVIDIA, “Multi-Process Service,” NVIDIA Corporation, Tech. Rep., Aug. 2019\. [Online]. Available: https://docs.nvidia.com/deploy/pdf/CUDA_Multi_Process_Service_Overview.pdf
  * [10] H. Carter Edwards, C. R. Trott, and D. Sunderland, “Kokkos: Enabling manycore performance portability through polymorphic memory access patterns,” _Journal of Parallel and Distributed Computing_ , vol. 74, no. 12, pp. 3202–3216, Dec. 2014. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S0743731514001257
  * [11] C. Miehe, “Numerical computation of algorithmic (consistent) tangent moduli in large-strain computational inelasticity,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 134, no. 3-4, pp. 223–240, Aug. 1996. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/0045782596010195
  * [12] A. J. Licup, A. Sharma, and F. C. MacKintosh, “Elastic regimes of subisostatic athermal fiber networks,” _Physical Review E_ , vol. 93, no. 1, Jan. 2016. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevE.93.012407
  * [13] M. R. Islam, G. Tudryn, R. Bucinell, L. Schadler, and R. C. Picu, “Stochastic continuum model for mycelium-based bio-foam,” _Materials & Design_, vol. 160, pp. 549–556, Dec. 2018. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0264127518307482
  * [14] S. Deogekar and R. C. Picu, “On the strength of random fiber networks,” _Journal of the Mechanics and Physics of Solids_ , vol. 116, pp. 1–16, Jul. 2018. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0022509618301285
  * [15] P. Underwood, “Dynamic Relaxation,” in _Computational Method for Transient Analysis_ , 1986, vol. 1, pp. 245–263. [Online]. Available: https://ci.nii.ac.jp/naid/10009833325/
  * [16] V. W. L. Chan, W. R. Tobin, S. Zhang, B. A. Winkelstein, V. H. Barocas, M. S. Shephard, and C. R. Picu, “Image-based multi-scale mechanical analysis of strain amplification in neurons embedded in collagen gel,” _Computer Methods in Biomechanics and Biomedical Engineering_ , pp. 1–16, Nov. 2018. [Online]. Available: https://www.tandfonline.com/doi/full/10.1080/10255842.2018.1538414
  * [17] “Kokkos C++ Performance Portability Programming EcoSystem: The Programming Model: Parallel Execution and Memory Abstraction - kokkos/kokkos,” Kokkos. [Online]. Available: https://github.com/kokkos/kokkos
  * [18] M. Korch and T. Werner, “Accelerating explicit ODE methods on GPUs by kernel fusion,” _Concurrency and Computation: Practice and Experience_ , vol. 30, no. 18, p. e4470, 2018. [Online]. Available: http://onlinelibrary.wiley.com/doi/abs/10.1002/cpe.4470
  * [19] G. Wang, Y. Lin, and W. Yi, “Kernel Fusion: An Effective Method for Better Power Efficiency on Multithreaded GPU,” in _2010 IEEE/ACM Int’l Conference on Green Computing and Communications Int’l Conference on Cyber, Physical and Social Computing_ , Dec. 2010, pp. 344–350.
  * [20] M. Wahib and N. Maruyama, “Scalable Kernel Fusion for Memory-Bound GPU Applications,” in _SC ’14: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ , Nov. 2014, pp. 191–202.
  * [21] J. Merson, “Kokkos-packed-data.” [Online]. Available: https://github.com/jacobmerson/kokkos-packed-data


# Sensing-Assisted Receivers for Resilient-By-Design 6G MU-MIMO Uplink
††thanks:  The authors acknowledge the financial support by the Federal
Ministry of Education and Research of Germany in the programme of “Souverän.
Digital. Vernetzt.”. Joint project 6G-life, project identification number:
16KISK002. This work was also supported by the Bavarian Ministry of Economic
Affairs, Regional Development and Energy within the project 6G Future Lab
Bavaria.

Vlad C. Andrei1, Xinyang Li2, Ullrich J. Mönich3 and Holger Boche4  Chair of
Theoretical Information Technology
Technical University of Munich
Munich, Germany
Email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We address the resilience of future 6G MIMO communications by considering an
uplink scenario where multiple legitimate transmitters try to communicate with
a base station in the presence of an adversarial jammer. The jammer possesses
full knowledge about the system and the physical parameters of the legitimate
link, while the base station only knows the UL-channels and the angle-of-
arrival (AoA) of the jamming signals. Furthermore, the legitimate transmitters
are oblivious to the fact that jamming takes place, thus the burden of
guaranteeing resilience falls on the receiver. For this case we derive one
optimal jamming strategy that aims to minimize the rate of the strongest user
and multiple receive strategies, one based on a lower bound on the achievable
signal-to-interference-to-noise-ratio (SINR), one based on a zero-forcing (ZF)
design, and one based on a minimum SINR constraint. Numerical studies show
that the proposed anti-jamming approaches ensure that the sum rate of the
system is much higher than without protection, even when the jammer has
considerably more transmit power and even if the jamming signals come from the
same direction as those of the legitimate users.

###### Index Terms:

Worst-case jammer, Multiple-Input-Multiple-Output (MIMO), 6G, Joint
communication and sensing (JCAS), Beamforming

††publicationid: pubid: 978-1-6654-0579-9/22/$31.00 ©2023 IEEE

##  I Introduction

###  I-A Motivation

While the first generations of mobile communications standards concerned
themselves with increasing data-rates and connectivity, the latter generations
focused on a wide array of use-cases and applications. With the currently
deployed 5G standard, very high data rates, reliability, low latency and
massive connectivity between a very large number of entities are made
possible, enabling applications such as the Tactile Internet for businesses,
virtual reality and advanced, remote-controlled robotics [1]. In the future 6G
standard, data rates and time sensitivity orders of magnitude higher than what
are currently available, as well as connectivity for all things and digital
twins are envisioned [2].

One of the key drivers towards enabling this massive digital transformation is
integrated, or joint communication and sensing (JCAS)[3], which offers radio
sensing services along with communications capabilities, all under a unified
hardware platform [4]. With the emergence of the virtualization not only of
data, but of whole physical objects at once, the question of trustworthiness
becomes crucial [5], which comprises, among other things the resilience of a
system.

While jamming can also be friendly [6], a malicious jammer, who has access to
all physical and system parameters of the environment, such as communication
protocols, signal processing, physical channels, etc. is considered. This is
also known as worst-case or smart jammer. The legitimate entities on the other
hand only have access to their own system and physical parameters, and have
very limited knowledge about the adversarial entity. In [5] it was shown that
in the classical communication setting (without sensing), the question whether
a denial-of-service attack has taken place is undecidable on a Turing machine
[7]. This means that the coding layer can not ensure resilience against such
attacks at the resource allocation level, and that other coordination
mechanisms, also known as common randomness (CR), are required [7]. In this
paper, we explore how to embed sensing information into a system architecture,
which needs to be resilient by design. Our main goal is to optimally leverage
sensing information and the spatial structure of the radio channel in order to
design transmit and receive filters which are resistant against worst-case
jammers. Note that this approach might naturally be interpreted in the CR
context, since the AoAs assure the coordination resources.

###  I-B Related Work and Main Contributions

While anti-jamming is a widely studied subject, to our knowledge there are no
works which address the problem of designing resilient transceivers against a
jammer who employs his optimal strategy by leveraging sensing information and
the spatial structure of the radio channel. In [8] the authors propose a
database approach to nullifying signals coming from the jammer direction using
massive planar arrays, and by only considering the jammer power and not the
its transmit strategy, In [9] and [10], the authors use reinforcement-learning
(RL) to optimize various system parameters under the simplifying assumption
that the jammer transmits uniformly over all streams.

The characterization of a worst-case jammer in the context of MIMO
communications has been done in terms of the optimal jammer strategy, albeit
without considering the spatial structure of the radio channel. The authors in
[11] derive the optimal transmit and jammer strategies using knowledge of the
transmitted symbols strategy at the jammer and show that under certain
conditions, this does not affect jammer performance. In [12] the authors prove
that there exists a lower bound on the achievable rate of communication
between the legitimate party, which does not depend on the jammer setup.
Furthermore, they characterize and give a closed form solution of the optimal
jamming strategy in the so-called jammer-dominant regime while in [13] the
authors derive the optimal jamming strategy in closed form, which minimizes
the signal-to-interference-and-noise-ratio (SINR) at the legitimate receiver
in a MIMO system. At last, in [14] the authors use game theoretical tools to
approach the problem of optimal resource allocation in the MIMO multiple
access (MAC) and broadcast channels (BC) under jammer conditions, and show
that the Nash equilibria of the resulting games always exist.

The only link between the concept of jamming and joint communication and
sensing (JCAS) we are aware of, are the works of [15] and [16], where the
authors employ knowledge about an eavesdroppers location in order to send
artificial noise in its direction, therefore increasing the secrecy rate.
Thus, our contributions are as follows:

  * •

We construct a communication model which exploits the spatial structure of the
radio channel and sensing information for a multi-user MIMO MAC scenario.

  * •

We propose an optimal jamming strategy which minimizes the achievable rate of
the strongest user, thus being robust against Successive Interference
Cancellation (SIC) [17].

  * •

Next we propose anti-jamming receive filters at the BS which only use AoA
information. The proposed designs are based on a lower bound on the individual
per-stream rate, a zero-forcing, and a minimum SINR constraint respectively.

  * •

We demonstrate by numerical experiments, that we can ensure constant,
satisfactory performance regardless of the number of antennas at the jammer.
We also demonstrate numerically that the achieved sum-rate mainly depends on
the AoA of the jamming signals and on the jamming power, and that the proposed
designs can ensure a non-zero sum rate even in the most extreme cases.

###  I-C Notation

Throughout this paper we denote the sets of natural, real and complex numbers
by $\mathbb{N}$, $\mathbb{R}$ and $\mathbb{C}$, respectively. We use lower
case letters for scalars $x$, bold lower-case letters for vectors $\bm{x}$ and
bold upper-case letters for matrices $\bm{X}$. We write sets as
$\\{x_{i}\\}_{i=1}^{S}$, where $i\in\mathbb{N}$ is the index and $S$ is the
cardinality of the set. The hermitian transponse and inverse of a matrix
$\bm{X}$ are denoted by $\bm{X}^{\mathsf{T}}$, $\bm{X}^{\mathsf{H}}$, and
$\bm{X}^{-1}$, respectively. The trace, Frobenius norm and rank of a matrix
$\bm{X}$ are written as $\operatorname*{tr}\left(\bm{X}\right)$,
$\lVert\bm{X}\rVert_{F}$ and $\text{rank}(\bm{X})$. The symbol $\succcurlyeq$
denotes the semi-ordering relationship on the cone of positive semi-definite
matrices. We use $\lVert\bm{x}\rVert_{2}$ to denote the Euclidean norm of a
vector and $\bm{e_{k}}$ for the $k$-th canonical basis vector in
$\mathbb{C}^{N}$. Here, $\bm{\mathrm{I}}$ denotes the identity matrix. By
$\bm{x}\sim\mathcal{N}\left(\bm{0},\bm{C}_{\bm{x}}\right)$ we say the random
vector $\bm{x}$ is a zero-mean, proper Gaussian random variable with
covariance matrix $\bm{C}_{\bm{x}}$, and by $\mathcal{U}(a,b)$ we denote the
uniform distribution on the interval $[a,b]\in\mathbb{R}$. Lastly, the mutual
information between the two random variables $\bm{x}$ and $\bm{y}$ is denoted
as $I(\bm{x};\bm{y})$ and the real part of $\bm{z}\in\mathbb{C}^{N}$ as
$\Re\left\\{\bm{z}\right\\}$.

##  II Preliminaries

###  II-A System Model

We assume a setup consisting of $K$ legitimate transmitters, a legitimate
receiver / base station (Bob) and jammer (Jimmy). The transmitters send their
data streams $s_{k}\sim\mathcal{N}(0,1)$ by applying precoding with
$\bm{w}_{k}\in\mathbb{C}^{N_{A_{k}}}$ through $N_{A_{k}}$ antennas. The
precoding vectors $\bm{w}_{k}$ satisfy the power constraint
$\left\lVert\bm{w}_{k}\right\|_{2}^{2}\leq P_{A_{k}}$. The resulting signals
$\bm{x}_{k}=\bm{w}_{k}s_{k}$ propagate through the legitimate channels
$\bm{H}_{k}\in\mathbb{C}^{N_{B}\times N_{A_{k}}}$ and arrive at Bob’s $N_{B}$
antennas, where we assume $N_{B}\geq K+1$. The signal is corrupted by white
noise $\bm{n}\sim\mathcal{N}\left(\bm{0},\sigma^{2}\bm{\mathrm{I}}\right)$, as
well as by the jammer signal
$\bm{z}\sim\mathcal{N}\left(\bm{0},\bm{C}_{\bm{z}}\right)\in\mathbb{C}^{N_{J}}$,
which propagates through the channel $\bm{G}\in\mathbb{C}^{N_{B}\times
N_{J}}$, and satisfies the power constraint
$\operatorname*{tr}\left(\bm{C}_{\bm{z}}\right)\leq P_{J}$. At last, Bob
applies the equalizers $\bm{v}_{k}^{\mathsf{H}}\in\mathbb{C}^{1\times N_{B}}$
to build an estimate $\hat{s}_{k}$ for each the transmitted streams $s_{k}$.
More formally:

$\displaystyle\bm{x}_{k}=\bm{w}_{k}s_{k}$ (1)
$\displaystyle\bm{y}=\sum_{k=1}^{K}\bm{H}_{k}\bm{x}_{k}+\bm{G}\bm{z}+\bm{n}$
(2)
$\displaystyle\hat{s}_{k}=\bm{v}_{k}^{\mathsf{H}}\bm{y}=\bm{v}_{k}^{\mathsf{H}}\left(\sum_{k=1}^{K}\bm{H}_{k}\bm{w}_{k}s_{k}\right)+\bm{v}_{k}^{\mathsf{H}}\bm{G}\bm{z}+\bm{v}_{k}^{\mathsf{H}}\bm{n}$
(3)

The legitimate channels $\bm{H}_{k}$ and the jammer channel $\bm{G}$ are
modelled as spatial beam-space channels, namely

$\displaystyle\bm{H}_{k}$
$\displaystyle=\sum_{l=1}^{L_{H_{k}}}b_{H_{k},l}\bm{a}_{B}(\theta_{H_{k},l})\bm{a}_{A_{k}}(\psi_{H_{k},l})^{\mathsf{H}}$
(4) $\displaystyle\bm{G}$
$\displaystyle=\sum_{l=1}^{L_{G}}b_{G,l}\bm{a}_{B}(\theta_{G,l})\bm{a}_{J}(\psi_{G,l})^{\mathsf{H}},$
(5)

with $L_{H_{k}}$, $L_{G}$ being the number of resolvable paths for each
channel, $\bm{a}_{S},\,S\in\\{A,B,J\\}$ the steering vectors, and
$\theta_{\cdot,l}$, $\psi_{\cdot,l}$, $b_{\cdot,l}$ the angles of arrival,
angles of departure and path gain, corresponding to the $l$-th resolvable path
in each channel, respectively.

For simplicity of analysis, we assume a 2D geometry, and thus all parties
employ uniform linear arrays (ULA), with the steering vector given by

$\displaystyle\bm{a}_{S}(\theta)=\begin{bmatrix}1&\mathrm{e}^{-\mathrm{j}\frac{2\pi
d}{\lambda_{c}}\cdot\sin\theta}&\dots&\mathrm{e}^{-\mathrm{j}\frac{2\pi
d}{\lambda_{c}}\cdot(N-1)\sin\theta}\end{bmatrix}^{\mathsf{T}}$ (6)

where $\lambda_{c}$, $d$ and $N$ denote the wavelength, element spacing, and
the number of elements. Defining the matrices

$\displaystyle\bm{A}_{B}(\bm{\theta}_{G})$
$\displaystyle=\begin{bmatrix}\bm{a}_{B}(\theta_{G,1})&\dots&\bm{a}_{B}(\theta_{G,L_{G}})\end{bmatrix}$
(7) $\displaystyle\bm{A}_{J}(\bm{\psi}_{G})$
$\displaystyle=\begin{bmatrix}\bm{a}_{J}(\psi_{G,1})&\dots&\bm{a}_{J}(\psi_{G,L_{G}})\end{bmatrix}$
(8) $\displaystyle\bm{B}_{G}$
$\displaystyle=\text{diag}\\{b_{G,l}\\}_{l=1}^{L_{G}}$ (9)

with $\bm{A}_{B}(\bm{\theta}_{G})\in\mathbb{C}^{N_{B}\times L_{G}}$,
$\bm{A}_{J}(\bm{\psi}_{G})\in\mathbb{C}^{N_{J}\times L_{G}}$,
$\bm{B}_{G}\in\mathbb{C}^{L_{G}\times L_{G}}$, the jammer channel $\bm{G}$
becomes

$\bm{G}=\bm{A}_{B}(\bm{\theta}_{G})\bm{B}_{G}\bm{A}_{J}(\bm{\psi}_{G})^{\mathsf{H}}.$
(10)

###  II-B Metrics

The achievable rates $R_{k}^{A}$ for each of the $K$ users under no particular
decoding order at the BS and the per-stream achievable rates $R_{B}^{k}$ are
given by

$\displaystyle R_{k}^{A}$
$\displaystyle=I(\bm{y};s_{k})=\log(1+\gamma_{k}^{A})$ (11) $\displaystyle
R_{k}^{B}$ $\displaystyle=I(\hat{s}_{k};s_{k})=\log(1+\gamma_{k}^{B}),$ (12)

with $\gamma_{k}^{A}$, $\gamma_{k}^{B}$ as in equations 13 and 14.

$\displaystyle\gamma_{k}^{A}$
$\displaystyle=\bm{w}_{k}^{\mathsf{H}}\bm{H}_{k}^{\mathsf{H}}\left(\sum\limits_{k^{\prime}\neq
k}\bm{H}_{k^{\prime}}\bm{w}_{k^{\prime}}\bm{w}_{k^{\prime}}^{\mathsf{H}}\bm{H}_{k^{\prime}}^{\mathsf{H}}+\bm{G}\bm{C}_{\bm{z}}\bm{G}^{\mathsf{H}}+\sigma^{2}\bm{\mathrm{I}}\right)^{-1}\bm{H}_{k}\bm{w}_{k}$
(13) $\displaystyle\gamma_{k}^{B}$
$\displaystyle=\dfrac{\bm{v}_{k}^{\mathsf{H}}\bm{H}_{k}\bm{w}_{k}\bm{w}_{k}^{\mathsf{H}}\bm{H}_{k}^{\mathsf{H}}\bm{v}_{k}}{\bm{v}_{k}^{\mathsf{H}}\left(\sum\limits_{k^{\prime}\neq
k}\bm{H}_{k^{\prime}}\bm{w}_{k^{\prime}}\bm{w}_{k^{\prime}}^{\mathsf{H}}\bm{H}_{k^{\prime}}^{\mathsf{H}}+\bm{G}\bm{C}_{\bm{z}}\bm{G}^{\mathsf{H}}+\sigma^{2}\bm{\mathrm{I}}\right)\bm{v}_{k}}$
(14)

The achievable sum-rates $R^{A}$, $R^{B}$ are then given by summing over all
$R_{k}^{A}$’s and $R_{k}^{B}$’s respectively.

##  III Jammer Model

We now present the smart jammer model used throughout the paper. We assume the
jammer has access to all system and physical parameters of all links, i.e.
$\\{\bm{w}_{k},\bm{H}_{k}\,\bm{v}_{k}\\}_{k=1}^{K}$, $\bm{G}$, $\sigma^{2}$,
but not to the transmitted symbols $s_{k}$, since it has been shown in [11]
that knowledge about $s$ does not increase jammer performance if
$\\{\bm{w}_{k},\bm{H}_{k}\\}_{k=1}^{K}$ are known. The next assumption is that
the jammer has more antennas than all parties, i.e. $N_{J}\geq
N_{A_{k}},N_{B}\,\forall k$, as well as more transmit power than all
transmitters $P_{J}\geq P_{A_{k}}\,\forall k$. Finally, we assume that the
disturbing transmission takes place in the so-called jammer-dominant regime.
This concept was rigorously formalized in [12], and roughly speaking, means
that the jammer power is much higher than the noise power, i.e.
$P_{J}\gg\sigma^{2}$.

The goal of the jammer is to minimize the sum-rate $R^{A}$ of all transmitters
under the power constraint $P_{J}$. In order to achieve this, we consider the
inequality

$R^{A}=\sum_{k=1}^{K}R_{k}^{A}\leq K\max_{k}R_{k}^{A}.$ (15)

Thus, if the rate of the “strongest” user can be brought to $0$, then the sum
rate is also $0$. Since $\log$ is a monotonically increasing function, the
problem can be cast as

$\displaystyle\bm{C}_{\bm{z}}^{*}=\arg\min_{\begin{subarray}{c}\bm{C}_{\bm{z}}\succcurlyeq\bm{0},\,\bm{C}_{\bm{z}}=\bm{C}_{\bm{z}}^{\mathsf{H}}\\\
\operatorname*{tr}\left(\bm{C}_{\bm{z}}\right)\leq
P_{J}\end{subarray}}\max_{k}\gamma_{k}^{A}.$ (16)

This is a convex optimization problem, since both the objective function
111The maximum over a family of convex functions is itself convex, the
function $\bm{x}^{\mathsf{H}}\bm{B}^{-1}\bm{x}$ is convex in symmetric
positive semidefinite $\bm{B}$ [18]. and the constraints are convex. Note that
we do not have to consider the equalizer at the receiver, since $R_{k}^{A}\geq
R_{k}^{B}\,\forall k$ by standard information theoretical arguments.
Furthermore, note that this jamming strategy can be greatly simplified if
Successive Interference Cancellation (SIC) [17] is used at the receiver, since
in this case, only the rate of one user needs to be minimized, rendering the
problem equivalent to that in [13]. In the following sections, we will see how
to ensure protection over a wide range of channel conditions when dealing with
a jammer employing the strategy presented in (16).

##  IV Resilient Receiver Design

###  IV-A General Considerations

We assume the legitimate receiver has perfect knowledge of the channels and
signal processing $\\{\bm{w}_{k},\bm{H}_{k}\\}_{k=1}^{K}$ at the transmitters,
and only has access to the AoAs
$\bm{\theta}_{G}=\\{\theta_{G,l}\\}_{l=1}^{L_{G}}$ of the impinging jammer
signals. Furthermore, the legitimate transmitters do not know the
communication is jammed, and do not cooperate with each other. That being
said, we assume for simplicity of analysis that the transmitters employ
singular value decomposition (SVD), namely

$\bm{w}_{k}=\mathcal{V}_{\text{max}}(\bm{H}_{k}),$ (17)

where $\mathcal{V}_{\text{max}}(\bm{H}_{k})$ denotes the right-singular vector
corresponding to the largest singular value of the matrix $\bm{H}_{k}$.

In the following, we will concentrate on the receiver and derive a bound on
the individual per-stream rate $R_{k}^{B}$, which depends, up to a constant,
solely on its setup. We will then use this lower bound in our designs.

Let $\bm{V}_{k}=\bm{v}_{k}\bm{v}_{k}^{\mathsf{H}}$ and consider

$\displaystyle\bm{v}_{k}^{\mathsf{H}}\left(\bm{G}\bm{C}_{\bm{z}}\bm{G}^{\mathsf{H}}+\sigma^{2}\bm{\mathrm{I}}\right)\bm{v}_{k}$
$\displaystyle=\operatorname*{tr}\left(\left(\bm{G}\bm{C}_{\bm{z}}\bm{G}^{\mathsf{H}}+\sigma^{2}\bm{\mathrm{I}}\right)\bm{V}_{k}\right)$
$\displaystyle=\operatorname*{tr}\left(\bm{G}\bm{C}_{\bm{z}}\bm{G}^{\mathsf{H}}\bm{V}_{k}+\sigma^{2}\bm{V}_{k}\right).$
(18)

Plugging in the expression in (10) and introducing the matrix
$\bm{E}=\bm{B}_{G}\bm{A}_{J}(\bm{\psi}_{G})^{\mathsf{H}}\bm{C}_{\bm{z}}^{1/2}$,
we obtain for the first term inside the trace

$\displaystyle\operatorname*{tr}\left(\bm{G}\bm{C}_{\bm{z}}\bm{G}^{\mathsf{H}}\bm{V}_{k}\right)=\operatorname*{tr}\left(\bm{A}_{B}(\bm{\theta}_{G})\bm{E}\bm{E}^{\mathsf{H}}\bm{A}_{B}(\bm{\theta}_{G})^{\mathsf{H}}\bm{V}_{k}\right)$
(19)
$\displaystyle=\operatorname*{tr}\left(\bm{E}\bm{E}^{\mathsf{H}}\bm{A}_{B}(\bm{\theta}_{G})^{\mathsf{H}}\bm{V}_{k}\bm{A}_{B}(\bm{\theta}_{G})\right)$
(20)
$\displaystyle\leq\left\lVert\bm{E}\right\|_{F}^{2}\operatorname*{tr}\left(\bm{A}_{B}(\bm{\theta}_{G})^{\mathsf{H}}\bm{V}_{k}\bm{A}_{B}(\bm{\theta}_{G})\right)$
(21)
$\displaystyle\leq\eta\operatorname*{tr}\left(\bm{A}_{B}(\bm{\theta}_{G})^{\mathsf{H}}\bm{V}_{k}\bm{A}_{B}(\bm{\theta}_{G})\right),$
(22)

with $\eta=P_{J}N_{J}L_{G}\left\lVert\bm{B}_{G}\right\|_{F}^{2}$. (20) follows
from the commutativity of trace, (21) follows from the fact that
$\operatorname*{tr}\left(\bm{AB}\right)\leq\operatorname*{tr}\left(\bm{A}\right)\operatorname*{tr}\left(\bm{B}\right)$
for $\bm{A}$, $\bm{B}$ symmetric positive-semidefinite [19, Section 1], and
lastly, (22) follows from applying the sub-multiplicativity of the Frobenius
norm twice, and from equations (4) and (7), respectively. With the shorthands,

$\displaystyle\bm{A}_{k}$
$\displaystyle=\bm{H}_{k}\bm{w}_{k}\bm{w}_{k}^{\mathsf{H}}\bm{H}_{k}^{\mathsf{H}}$
(23) $\displaystyle\bm{B}_{k}$ $\displaystyle=\sum\limits_{k^{\prime}\neq
k}\bm{H}_{k^{\prime}}\bm{w}_{k^{\prime}}\bm{w}_{k^{\prime}}^{\mathsf{H}}\bm{H}_{k^{\prime}}^{\mathsf{H}}+\sigma^{2}\bm{\mathrm{I}},$
(24)

we have

$\displaystyle\gamma_{k}^{B}\geq\dfrac{\bm{v}_{k}\bm{A}_{k}\bm{v}_{k}}{\bm{v}_{k}^{\mathsf{H}}\left(\bm{B}_{k}+\eta\bm{A}_{B}(\bm{\theta}_{G})\bm{A}_{B}(\bm{\theta}_{G})^{\mathsf{H}}\right)\bm{v}_{k}}\overset{\Delta}{=}\tilde{\gamma}_{k}^{B}.$
(25)

We note that equality in (25) is achieved if $P_{J}=0$ (no jammer) or if
$\bm{A}_{B}(\bm{\theta}_{G})^{\mathsf{H}}\bm{V}_{k}=\bm{0}$, i.e. equalizer
$\bm{v}_{k}$ lies in a subspace orthogonal to the array manifold
$\bm{A}_{B}(\bm{\theta}_{G})$.

We note that aside from $\eta$, all other terms in (25) are known at the
receiver. In the jammer-dominant regime the term
$\eta\operatorname*{tr}\left(\bm{A}_{B}(\bm{\theta}_{G})^{\mathsf{H}}\bm{V}_{k}\bm{A}_{B}(\bm{\theta}_{G})\right)$
is much greater that $\sigma^{2}$ (since it scales linearly in $P_{J}$ and
$N_{J}$), effectively dominating the influence of the noise.

###  IV-B Closed-Form Design

The receivers’ goal is to ensure reliable communication for all uplink
participants, i.e. we are interested in maximizing the individual per-stream
SINRs $\gamma_{k}^{B}$ or at least ensure a certain quality of service (QoS),
e.g. the $\gamma_{k}^{B}$’s lie above a given threshold $\gamma_{0}$.

The maximization of the lower bound in Eq. (25) is a Rayleigh quotient
maximization problem with the standard solution being given by:

$\displaystyle\bm{v}_{k}=\left(\bm{B}_{k}+\eta\bm{A}_{B}(\bm{\theta}_{G})\bm{A}_{B}(\bm{\theta}_{G})^{\mathsf{H}}\right)^{-1}\bm{H}_{k}\bm{w}_{k}$
(26)

for all $k=1,\dots,K$. This is the most simple way to leverage AoA information
at the receiver, since the solution is given in closed form and
computationally tractable, since
$\bm{A}_{B}(\bm{\theta}_{G})\bm{A}_{B}(\bm{\theta}_{G})^{\mathsf{H}}$ can be
easily computed without matrix multiplications by only considering the phase
shifts between the impinging jamming signals.

Note, that the resulting filter still depends on $\eta$, which in general is
not available at the receiver. The most easy way to address this problem is to
consider $\eta$ as a hyperparameter which controls how much the jamming
signals are suppressed. Indeed, this is justified by the fact that in the
jammer-dominant regime, the contribution of the noise variance $\sigma^{2}$ is
small compared to the one of the jammer, i.e
$\bm{v}_{k}^{\mathsf{H}}\left(\eta\bm{A}_{B}(\bm{\theta}_{G})\bm{A}_{B}(\bm{\theta}_{G})^{\mathsf{H}}+\sigma^{2}\bm{\mathrm{I}}\right)\bm{v}_{k}\approx\eta\bm{v}_{k}^{\mathsf{H}}\bm{A}_{B}(\bm{\theta}_{G})\bm{A}_{B}(\bm{\theta}_{G})^{\mathsf{H}}\bm{v}_{k}$.
In the following, we will propose other methods to overcome this dependency.

###  IV-C Zero Forcing Design

As it can be seen in Section IV-B, the main contribution to the jamming
signals is given by the term

$\displaystyle\beta_{k}=\left\lVert\bm{A}_{B}(\bm{\theta}_{G})^{\mathsf{H}}\bm{v}_{k}\right\|_{2}^{2}$
(27)

which can be naturally interpreted as the receive beampattern of the receiver
antenna array [20]. This indicates how much the signals coming from the jammer
direction are emphasized, thus the anti-jamming criterion can be formalized as
a beampattern minimization problem in the same spirit as [21, 15, 22].
Furthermore, a natural way to ensure good per-stream SINRs for all users is to
remove the interference, which can be realized with a zero-forcing (ZF)
design. Thus, the beampattern minimization problem under ZF constraint can be
formulated analogously to [21] as

$\displaystyle\min_{\begin{subarray}{c}\bm{v}_{k}\\\
k=1,\dots,K\end{subarray}}\sum_{k=1}^{K}\beta_{k}\quad\text{s.t.}\quad\bm{P}^{\mathsf{H}}\bm{v}_{k}=\bm{e}_{k}\,\forall
k$ (28)

where we defined

$\bm{P}=\begin{bmatrix}\bm{H}_{1}\bm{w}_{1}&\dots&\bm{H}_{K}\bm{w}_{K}\end{bmatrix}\in\mathbb{C}^{N_{B}\times
K}$ (29)

and $\bm{e}_{k}$ is $k$-th standard basis vector in $\mathbb{R}^{K}$. Note
that this is a convex optimization problem easily solvable in polynomial time.

Furthermore, if the matrix
$\bm{A}_{B}(\bm{\theta}_{G})\bm{A}_{B}(\bm{\theta}_{G})^{\mathsf{H}}$ has
full-rank, we can also give an analytic solution. Letting
$\bm{X}=\bm{A}_{B}(\bm{\theta}_{G})\bm{A}_{B}(\bm{\theta}_{G})^{\mathsf{H}}$,
we construct the Lagrangian

$\displaystyle\mathcal{L}\left(\bm{v}_{k},\bm{\mu}_{k}\right)=\sum_{k=1}^{K}\bm{v}_{k}^{\mathsf{H}}\bm{X}\bm{v}_{k}+2\Re\left\\{\bm{\mu}_{k}^{\mathsf{H}}\left(\bm{P}^{\mathsf{H}}\bm{v}_{k}-\bm{e}_{k}\right)\right\\}$
(30)

where $\bm{\mu}_{k}\in\mathbb{C}^{K}$ are the Lagrange multipliers. Since the
ZF-constraint is an equality constraint, the Karush-Kuhn-Tucker conditions are
both necessary and sufficient. Computing the gradients wrt. $\bm{v}_{k}$ and
$\bm{\mu}_{k}$, setting them to $\bm{0}$, and solving for $\bm{v}_{k}$, we
obtain

$\displaystyle\bm{v}_{k}=\bm{X}^{-1}\bm{P}\left(\bm{P}^{\mathsf{H}}\bm{X}^{-1}\bm{P}\right)^{-1}\bm{e}_{k}$
(31)

By standard linear algebra $\bm{X}$ has full rank iff $L_{G}\geq N_{B}$, which
might pose a problem if the propagation channel of the jammer shows
$L_{G}<N_{B}$ dominant paths. In that case, one can easily “pad” the resolved
AoAs $\bm{\theta_{G}}$ with a set of different angles
$\\{\phi_{G,l}\\}_{l=1}^{N_{B}-L_{G}}$ thus enlarging the grid over which the
beampattern should be minimized. This idea can be naturally extended to the
designs considered in Sections IV-B and IV-D.

###  IV-D Minimum SINR Design

Another way to approach this problem is to minimize the influence of the
jamming signals, while ensuring a minimum SINR $\gamma_{0}$ in the jammer-free
case. Using this method one can address the disadvantages ZF-filters generally
have in the low SNR regime, i.e. $\sigma^{2}\gg 0$. Rewriting the beampattern
$\beta_{k}$ as

$\displaystyle\beta_{k}=\left\lVert\bm{A}_{B}(\bm{\theta}_{G})^{\mathsf{H}}\bm{v}_{k}\right\|_{2}^{2}=\operatorname*{tr}\left(\bm{A}_{B}(\bm{\theta}_{G})^{\mathsf{H}}\bm{V}_{k}\bm{A}_{B}(\bm{\theta}_{G})\right)$
(32)

we cast the resulting problem as

$\displaystyle\min_{\begin{subarray}{c}\bm{V}_{k}\succcurlyeq\bm{0}\,\bm{V}_{k}=\bm{V}_{k}^{\mathsf{H}}\\\
\text{rank}(\bm{V}_{k})=1\\\
k=1,\dots,K\end{subarray}}\,\sum_{k=1}^{K}\beta_{k}\quad\text{s.t.}\quad\dfrac{\operatorname*{tr}\left(\bm{A}_{k}\bm{V}_{k}\right)}{\operatorname*{tr}\left(\bm{B}_{k}\bm{V}_{k}\right)}\geq\gamma_{0}$
(33)

with $\bm{V}_{k}=\bm{v}_{k}\bm{v}_{k}^{\mathsf{H}}$ as before and $\bm{A_{k}}$
and $\bm{B_{k}}$ as in equations (23) and (24). Note that this problem is non-
convex due to the rank-$1$ constraint on $\bm{V}_{k}$. In order to efficiently
solve this problem, we apply semidefinite relaxation [23], i.e. we relax the
rank-$1$ constraint, solve the resulting convex problem and recover the
solution via eigenvalue decomposition (EVD).

##  V Numerical Results

###  V-A Methodology

In this section we present the performance of the proposed designs with
respect to the jammer setup and channel conditions. To this end, we consider
three parameters, namely the number of antennas at the jammer $N_{J}$, the
power budget of the jammer $P_{J}$, and the AoAs $\bm{\theta}_{G}$ of the
impinging signals. We shall compare the sum-rates obtained by the proposed
designs with the sum-rates in the jammer-free case, as well as in the case
when the parties do not ensure any jammer protection. We consider a setup of
$K=3$ users, each with $N_{A_{k}}=8$ and $P_{A_{k}}=5$\mathrm{dBm}$$. The
number of antennas at the legitimate transmitter $N_{B}$ is kept fixed at
$16$, and the noise variance $\sigma^{2}=-10$\mathrm{dB}$$. In order to
generate the channel matrices, we draw all path gains and AoDs in (4) randomly
from $\mathcal{N}(0,1)$ and $\mathcal{U}(-5^{\circ},5^{\circ})$. The AoAs are
generated as $\theta_{G,l}=\theta_{J}+\phi_{l}$ and
$\theta_{H_{k},l}=\theta_{A_{k}}+\omega_{l}$, with $\phi_{l}$ and $\omega_{l}$
drawn i.i.d from $\mathcal{U}(-5^{\circ},5^{\circ})$, and
$\theta_{A_{k}}=\\{-10^{\circ},0^{\circ},10^{\circ}\\}$ from each other. The
direction dependency of the channels are then characterized by the central AoA
$\theta_{J}$ of the impinging jamming signal.

###  V-B Simulation Studies

We first start by assessing the influence of the number of antennas at the
jammer on the communication quality. To this end, we vary the number of
antennas at the jammer $N_{J}$ and average the computed sum-rates over
$\theta_{J}$ and $P_{J}$. We first observe that the obtained sum rates are
independent of the number of antennas at the jammer, since the performance
curves remain relatively flat along the whole range. All of the proposed
designs offer a much higher rate than in the cases without protection, with a
minimum gain of $10\text{bits}/$\mathrm{s}$/$\mathrm{Hz}$$. In terms of
performance ranking, the ZF-based design from equation (28) (denoted as “ZF”
in Figures 1-3) performs best, followed by the optimal analytic filter based
on the surrogate objective in equation (26) (“Analytic” in Figures 1-3). The
filter based on the minimum QoS (Equation (33)), denoted “MinSINR” in the
plots, comes last.

We attribute this ranking to the fact that the filters “Analytic” and
“MinSINR” are much more sensitive to the jamming power, since the former is
coupled to it by the constant $\eta$, and the latter is based on an objective,
that should ensure a minimum SINR at all times. This might not be fully-
possible when jamming powers are high. In order to test this conjecture, we
vary the jammer power $P_{J}$ and average the computed sum-rates over
$\theta_{J}$ and $N_{J}$. Indeed, while the filters based on the minimum SINR
criterion and analytic formula decline sharply over the studied range, the
filter based on the ZF criterion stays approximately constant over the whole
range, only declining for high jammer powers.

Finally, we plot the performance of the proposed designs with respects to the
AoA $\theta_{J}$ of the impinging jamming signals in Figure 3. Note, that the
AoAs of the users are marked with a vertical, red line. We observe that the
sum-rate decreases as the AoAs of the jamming signals become similar to those
of the legitimate parties. We furthermore observe, that the jammer can not
drive the sum-rate to $0$ in the studied configuration. In terms of
performance, the filter based on the analytic formula shows slightly better
performance, than the ZF approach, in the case that the jamming and legitimate
AoAs are similar.

$20$$40$$60$$80$$100$$120$$10$$20$$30$$40$$\displaystyle N_{J}$Rate (bits / s
/ Hz)No JammerNo ProtectionAnalyticMinSINRZF Figure 1: Averaged sum-rates for
different number of antennas at the jammer, $\gamma_{0}=20$\mathrm{dB}$$,
$\eta=1$. $10$$20$$30$$40$$50$$60$$0$$10$$20$$30$$40$$\displaystyle
P_{J}($\mathrm{dBm}$)$Rate (bits / s / Hz)No JammerNo
ProtectionAnalyticMinSINRZF Figure 2: Averaged sum-rates for different jammer
powers, $\gamma_{0}=20$\mathrm{dB}$$, $\eta=1$.
$-30$$-20$$-10$$0$$10$$20$$30$$0$$10$$20$$30$$40$User 1User 2User
3$\displaystyle\theta_{J}$ (deg.)Rate (bits / s / Hz)No JammerNo
ProtectionAnalyticMinSINRZF Figure 3: Averaged sum-rates for different
directions of arrival, $\gamma_{0}=20$\mathrm{dB}$$, $\eta=1$.

##  VI Conclusions

In this work, we derived a communication model which exploits the physical
structure of the MIMO-MAC in the context of anti-jamming We derived one
optimal strategy at the jammer, which seeks to minimize the rate of the
strongest user. We furthermore proposed an array of methods for protection
against this worst-case jammer, possessing a far more powerful setup than the
legitimate parties. We have shown experimentally, that we can achieve a non-
zero sum-rate even if the jamming and legitimate signals come from roughly the
same directions, and even if the jammer transmits with significantly more
power. We conclude, that resilience can not be fully guaranteed by only using
the spatial dimension of communication. Thus, our future work will be
dedicated to the design of optimal signaling and frame structures, as well as
to the study of optimal jamming strategies in these cases.

## References

  * [1] G. P. Fettweis, “The Tactile Internet: Applications and Challenges,” _IEEE Vehicular Technology Magazine_ , vol. 9, no. 1, pp. 64–70, 2014.
  * [2] H. Tataria, M. Shafi, A. F. Molisch, M. Dohler, H. Sjöland, and F. Tufvesson, “6G Wireless Systems: Vision, Requirements, Challenges, Insights, and Opportunities,” _Proceedings of the IEEE_ , vol. 109, no. 7, pp. 1166–1199, 2021.
  * [3] J. A. Zhang, M. L. Rahman, K. Wu, X. Huang, Y. J. Guo, S. Chen, and J. Yuan, “Enabling Joint Communication and Radar Sensing in Mobile Networks—A Survey,” _IEEE Communications Surveys & Tutorials_, vol. 24, pp. 306–345, 2022.
  * [4] J. A. Zhang, F. Liu, C. Masouros, R. W. Heath, Z. Feng, L. Zheng, and A. Petropulu, “An overview of signal processing techniques for joint communication and radar sensing,” _IEEE Journal of Selected Topics in Signal Processing_ , 2021.
  * [5] G. P. Fettweis and H. Boche, “On 6G and Trustworthiness,” _Commun. ACM_ , vol. 65, no. 4, p. 48–49, Mar 2022.
  * [6] W. Shen, P. Ning, X. He, and H. Dai, “Ally Friendly Jamming: How to Jam Your Enemy and Maintain Your Own Wireless Connectivity at the Same Time,” in _IEEE Symposium on Security and Privacy_ , 2013, pp. 174–188.
  * [7] H. Boche, R. Schaefer, and H. Poor, “Denial-of-Service Attacks on Communication Systems: Detectability and Jammer Knowledge,” _IEEE Transactions on Signal Processing_ , vol. 68, pp. 3754–3768, May 2020.
  * [8] M. Chehimi, E. Yaacoub, A. Chehab, and M. Al-Husseini, “Physical Layer Anti-jamming Technique Using Massive Planar Antenna Arrays,” in _International Wireless Communications and Mobile Computing (IWCMC)_ , 2020, pp. 1740–1745.
  * [9] L. Xiao, Y. Li, C. Dai, H. Dai, and H. V. Poor, “Reinforcement Learning-Based NOMA Power Allocation in the Presence of Smart Jamming,” _IEEE Transactions on Vehicular Technology_ , vol. 67, no. 4, pp. 3377–3389, 2018.
  * [10] H. Yang, Z. Xiong, J. Zhao, D. Niyato, Q. Wu, H. V. Poor, and M. Tornatore, “Intelligent Reflecting Surface Assisted Anti-Jamming Communications: A Fast Reinforcement Learning Approach,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 3, pp. 1963–1974, 2021.
  * [11] A. Kashyap, T. Basar, and R. Srikant, “Correlated jamming on MIMO Gaussian fading channels,” _IEEE Transactions on Information Theory_ , vol. 50, no. 9, pp. 2119–2123, 2004.
  * [12] J. Gao, S. A. Vorobyov, H. Jiang, and H. V. Poor, “Worst-Case Jamming on MIMO Gaussian Channels,” _IEEE Transactions on Signal Processing_ , vol. 63, no. 21, pp. 5821–5836, 2015.
  * [13] Q. Liu, M. Li, X. Kong, and N. Zhao, “Disrupting MIMO Communications With Optimal Jamming Signal Design,” _IEEE Transactions on Wireless Communications_ , vol. 14, no. 10, pp. 5313–5325, 2015.
  * [14] M. H. Brady, M. Mohseni, and J. M. Cioffi, “Spatially-correlated Jamming in Gaussian Multiple Access and Broadcast Channels,” in _40th Annual Conference on Information Sciences and Systems_. IEEE, 2006, pp. 1635–1639.
  * [15] N. Su, F. Liu, and C. Masouros, “Secure Radar-Communication Systems With Malicious Targets: Integrating Radar, Communications and Jamming Functionalities,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 1, pp. 83–95, 2021.
  * [16] N. Su, F. Liu, Z. Wei, Y.-F. Liu, and C. Masouros, “Secure Dual-Functional Radar-Communication Transmission: Exploiting Interference for Resilience Against Target Eavesdropping,” _IEEE Transactions on Wireless Communications_ , vol. 21, no. 9, pp. 7238–7252, 2022.
  * [17] X. Zhang and M. Haenggi, “The Performance of Successive Interference Cancellation in Random Wireless Networks,” _IEEE Transactions on Information Theory_ , vol. 60, no. 10, pp. 6368–6388, 2014.
  * [18] S. Boyd and L. Vandenberghe, _Convex Optimization_. Cambridge University Press, March 2004.
  * [19] K. Shebrawi and H. Albadawi, “Trace Inequalities for Matrices,” _Bulletin of the Australian Mathematical Society_ , vol. 87, no. 1, p. 139–148, 2013.
  * [20] H. Van Trees, _Optimum Array Processing: Part IV of Detection, Estimation, and Modulation Theory_ , ser. Detection, Estimation, and Modulation Theory. Wiley, 2004.
  * [21] F. Liu, C. Masouros, A. Li, H. Sun, and L. Hanzo, “MU-MIMO Communications With MIMO Radar: From Co-Existence to Joint Transmission,” _IEEE Transactions on Wireless Communications_ , vol. 17, no. 4, pp. 2755–2770, 2018\.
  * [22] J. A. Zhang, F. Liu, C. Masouros, R. W. Heath, Z. Feng, L. Zheng, and A. Petropulu, “An Overview of Signal Processing Techniques for Joint Communication and Radar Sensing,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 15, no. 6, pp. 1295–1315, 2021.
  * [23] Z.-q. Luo, W.-k. Ma, A. M.-c. So, Y. Ye, and S. Zhang, “Semidefinite Relaxation of Quadratic Optimization Problems,” _IEEE Signal Processing Magazine_ , vol. 27, no. 3, pp. 20–34, 2010.


* [45] R. Stanley, _Enumerative combinatorics_ , 2nd ed., vol. 1, Cambridge University Press, 2011.
  * [46] R. Trombetti and F. Zullo, _On maximum additive hermitian rank-metric codes_ , Journal of Algebraic Combinatorics 54 (2021), no. 1, 151–171.
  * [47] B. Wu and Z. Liu, _Linearized polynomials over finite fields revisited_ , Finite Fields and Their Applications 22 (2013), 79–100.
  * [48] Y. Zhou, _On equivalence of maximum additive symmetric rank-distance codes_ , Designs, Codes and Cryptography 88 (2020), no. 5, 841–850.