# Compression of volume-surface integral equation matrices via Tucker
decomposition for magnetic resonance applications

Ilias I. Giannakopoulos, , Georgy D. Guryev, José E. C. Serrallés, , Ioannis
P. Georgakis, Luca Daniel, , Jacob K. White, , Riccardo Lattanzi This work was
supported by NIH R01 EB024536 and by NSF 1453675. It was performed under the
rubric of the Center for Advanced Imaging Innovation and Research (CAI2R,
www.cai2r.net), a NIBIB Biomedical Technology Resource Center (NIH P41
EB017183).Ilias I. Giannakopoulos, Ioannis P. Georgakis and Riccardo Lattanzi
are with Center for Advanced Imaging Innovation and Research (CAI2R),
Department of Radiology, New York University Grossman School of Medicine, NY,
USA.Georgy D. Guryev, José E.C. Serrallés, Luca Daniel, and Jacob K. White are
with the Research Laboratory of Electronics, Department of Electrical
Engineering and Computer Science, Massachusetts Institute of Technology,
Cambridge, MA, USA.Riccardo Lattanzi is also with the Bernard and Irene
Schwartz Center for Biomedical Imaging, Department of Radiology, New York
University Grossman School of Medicine, NY, USA and the Vilcek Institute of
Graduate Biomedical Sciences, New York University Grossman School of Medicine,
NY, USA.

###### Abstract

In this work, we propose a method for the compression of the coupling matrix
in volume-surface integral equation (VSIE) formulations. VSIE methods are used
for electromagnetic analysis in magnetic resonance imaging (MRI) applications,
for which the coupling matrix models the interactions between the coil and the
body. We showed that these effects can be represented as independent
interactions between remote elements in 3D tensor formats, and subsequently
decomposed with the Tucker model. Our method can work in tandem with the
adaptive cross approximation technique to provide fast solutions of VSIE
problems. We demonstrated that our compression approaches can enable the use
of VSIE matrices of prohibitive memory requirements, by allowing the effective
use of modern graphical processing units (GPUs) to accelerate the arising
matrix-vector products. This is critical to enable numerical MRI simulations
at clinical voxel resolutions in a feasible computation time. In this paper,
we demonstrate that the VSIE matrix-vector products needed to calculate the
electromagnetic field produced by an MRI coil inside a numerical body model
with $1$ mm3 voxel resolution, could be performed in $\sim 33$ seconds in a
GPU, after compressing the associated coupling matrix from $\sim 80$ TB to
$\sim 43$ MB.

###### Index Terms:

Cross approximation, Global Maxwell Tomography, graphical processing unit,
magnetic resonance imaging, matrix-vector product, Tucker decomposition,
volume-surface integral equation.

##  I Introduction

Magnetic resonance (MR) imaging (MRI) provides high-resolution images of the
interior anatomical and physiological structure of the human body, with
exquisite soft-tissue contrast. The quality of MR images, as well as the
achievable spatial and temporal resolution, depend on the available signal-to-
noise ratio (SNR). SNR increases with the main magnetic field strength. This
fact motivated the recent development of $7$ Tesla (T) clinical MR scanners
and research-only scanners with field strengths as high as $11.7$ T [1]. At
ultra-high-field (UHF) MRI ($\geq 7$ T), the radio frequency (RF) wavelength
is short. This results in strong interactions between biological tissues and
the electromagnetic (EM) field generated by the RF coils [2, 3, 4, 5]. Such
interactions could compromise image quality and patient safety. To address
these issues, EM modeling is often used to predict and manipulate the EM field
distribution during RF coil design.

Integral equation (IE) methods are suitable options for EM analysis in MRI.
First, they do not suffer from grid dispersion errors [6, 7], in contrast with
the finite-difference-time-domain (FDTD) and finite-element-methods (FEM),
because the Green’s functions in the underlying IE formulation act as an exact
EM field propagator from a source to an observation point. Second, for the
case of single-frequency problems, IE algorithms can be extensively customized
with the use of numerical linear algebra techniques for fast and accurate
simulations, tailored to specific applications [8, 9, 10, 11].

For example, the MAgnetic-Resonance Integral Equation (MARIE) suite [11, 12]
was developed to numerically compute the EM field distribution generated by RF
coils in the human body during MRI. MARIE combines surface and volume integral
equations (SIE,VIE), employing a triangular tessellation for the RF coils’
conductors and a uniform voxelized grid discretization for the body models.
RWG basis functions [13] and polynomial basis functions [12, 14] are used to
compute the unknowns of the surface and volume IE, respectively. Matrix-vector
products are accelerated using the fast Fourier transform (FFT).

The VIE computational engine of MARIE has been recently employed for the
forward problem in Global Maxwell Tomography (GMT) [15], a technique that
iteratively solves an ill-conditioned inverse problem to extract electrical
properties from volumetric MR measurements. In the first experimental
demonstration of GMT with a uniform phantom, constant incident fields were
used for all iterations [15]. More recently, it was shown in simulation that
GMT could accurately reconstruct brain electrical properties at $7$ T using a
tailored RF coil array [16]. However, in order to confirm this with in-vivo
experiments, the currents on the coil conductors cannot be simply considered
constant as in the initial experiment with a uniform phantom. Instead, the
incident fields must be updated at each iteration of GMT to account for
changes in the sample electrical properties distribution. Therefore, GMT must
be implemented with a volume-surface IE (VSIE) framework, in which a coupling
matrix represents the coil-body interactions in the IE system of equations
[11]. Such approach requires a large amount of memory, which could prevent
using clinically-relevant voxel resolutions and fine coil meshes.

The aim of this work is to use Tucker decomposition [17] to perform a column-
wise compression of the VSIE coupling matrix, in order to limit the associated
memory demand and enable the computation of the relevant matrix-vector
products in GPUs. Our approach was motivated by previous work [10] on the
reduction of the memory footprint of FFT-based VIE Green’s function tensors
and the acceleration of matrix-vector products in VIE using GPU. Tucker
decomposition belongs to a larger family of tensor decompositions and have
been used successfully in the past for matrix compression inside IE-based
simulations for EM applications. Examples include EM simulations of realistic
body models simulations for UHF MRI [18, 10] and capacitance extraction [19,
20, 21]. Other tensor decompositions could be used [22, 23, 24], but for the
intrinsic 3D nature of the problem at hand, Tucker optimizes operations and
memory complexity. In cases where the coil is placed far from the scatterer,
the coupling matrix can be first compressed with a 2D cross approximation
method [25, 26, 27] and then further compressed by applying our proposed
technique to the resulting matrices. Towards this direction, we developed an
algorithm based on the adaptive cross approximation (ACA) [28, 29] to
efficiently perform our compression approach within the iterative loop of ACA.

The remainder of this paper is organized as follows. In Section II, we
summarize the relevant technical background and equations related to the
standard VSIE method. We also show, as an example application, the use of the
VSIE coupling matrix in the GMT inverse problem formulation. In addition, we
outline the Tucker decomposition and the ACA method. In Section III, we
introduce the Tucker-based assembly technique of the coupling matrix along
with the novel algorithms for the matrix-vector products implementation.
Moreover, we present a new algorithm for a memory friendly assembly of the
coupling matrix, based on the combination of ACA with the Tucker
decomposition. Section IV describes a series of numerical experiments aimed at
investigating the multilinear rank and the memory requirements of the
compressed matrix in different scenarios, along with time footprint of the
matrix-vector product for various mesh discretizations. Section V discusses
the results. Section VI summarizes the work and provides a few take home
points. The following TABLE I lists the notations used in this work.

TABLE I: Notation Notation | Description
---|---
$a$ | Scalar
$\bm{a}$ | Vector in $\mathbb{C}^{3}$
$\mathbf{a}$ | Vector in $\mathbb{C}^{n}$
$A$ | Matrix in $\mathbb{C}^{n_{1}\times n_{2}}$
$A^{T}$ | Transpose of matrix
$A^{*}$ | Conjugate transpose of matrix
${\mathcal{A}}$ | Tensor in $\mathbb{C}^{n_{1}\times n_{2}\times n_{3}}$
$\mathscr{A}$ | Struct containing Tucker compressed tensors
${\cal{A}}$ | Operator acting on vectors in $\mathbb{C}^{3}$
$\times_{i}$ | n-mode products
${\mathrm{i}}$ | Imaginary unit ${\mathrm{i}}^{2}=-1$

##  II Technical Background

###  II-A Volume-Surface Integral Equation

####  II-A1 Coupled Linear System

Let us consider an IE-based solver for the EM wave analysis in MRI
applications. The body domain $\Omega$ can be modeled with the current-based
VIE (JVIE), while the conducting surfaces (receive and transmit RF coil
arrays, coil shields, and gradient coils) with the SIE. The IEs are solved
with the Method of Moments (MoM) [30]. The resulting system of equations can
be written in the following block matrix form.

$\begin{bmatrix}Z_{\rm cc}&T^{T}_{\rm bc}\\\ T_{\rm bc}&Z_{\rm
bb}\end{bmatrix}\begin{bmatrix}J_{\rm s}\\\ J_{\rm
p}\end{bmatrix}=\begin{bmatrix}V\\\ 0\end{bmatrix}.$ (1)

Here, $Z_{\rm cc}\in\mathbb{C}^{m\times m}$ is the Galerkin matrix that models
interactions between the coil’s discretization elements (triangles’ common
edges), with the aid of the free-space Green’s function that appears in the
SIE formulation. It associates the equivalent surface currents on the coil
($J_{\rm s}\in\mathbb{C}^{m\times p}$), with the voltage excitation matrix
$V\in\mathbb{C}^{m\times p}$. The coil conductors are modeled with triangular
meshes, and the unknown surface equivalent currents are approximated with RWG
basis functions [13]. $m$ is the number of RWG, or non-boundary, triangle
edges that appear on the mesh, whereas $p$ is the number of excitation ports
(i.e., number of coil’s channels) of the conducting surfaces.

The matrix $Z_{\rm bb}\in\mathbb{C}^{qn_{\rm v}\times qn_{\rm v}}$ is another
Galerkin matrix which models the EM interactions of an external volumetric EM
field, produced by the coil, with the body. Specifically, the matrix relates
the polarization currents in $\Omega$ to an external incident EM field
produced by the conducting surfaces. $n_{\rm v}$ is the number of voxels and
$q$ the number of basis functions per voxel. Differently than $Z_{\rm cc}$,
$Z_{\rm bb}$ requires a large amount of memory, even for coarse resolutions.
To handle that, $\Omega$ can be discretized over a voxelized uniform grid,
giving $Z_{\rm bb}$ a three-level Block-Toeplitz with Toeplitz Blocks (BTTB)
structure. As a result, only the defining columns of the BTTB matrix need to
be stored and the matrix-vector product can be accelerated using the FFT, as
in [31, 32, 33, 34, 35, 36, 37, 12]. The unknown polarization currents
($J_{\rm p}\in\mathbb{C}^{qn_{\rm v}\times p}$) can be discretized with
polynomial basis functions, either piecewise constant [12] (PWC, $3$ unknowns
per voxel) or piecewise linear [14] (PWL, $12$ unknowns per voxel), and a
single-voxel support.

The presence of conductive tissue near the coil conductors perturbs $J_{\rm
s}$ from their values in free-space. In fact, the voltage excitations at the
coil’s ports create incident EM fields that scatter from the dielectric body
back to the coil conductors, changing their current distribution. The coupling
matrix $T_{\rm bc}\in\mathbb{C}^{qn_{\rm v}\times m}$ is used to account for
this effect, by modeling the coupling interactions between the dyadic Green’s
function [38] of the SIE and the VIE formulations. Specifically, in equation
(1), $T_{\rm bc}$ maps electric surface equivalent currents to electric fields
through the ${\cal{N}}$ Green’s function operator of VIE:

${\cal{N}}\left(\bm{s}\right)\triangleq\nabla\times\nabla\times\int\limits_{\Omega}g\left(\bm{r}-\bm{r}^{\prime}\right)\bm{s}\left(\bm{r}^{\prime}\right)d^{3}\bm{r}^{\prime}.$
(2)

$g$ is the free-space Green’s function, or fundamental Helmholtz solution, and
it is equal to

$g\left(\bm{r}-\bm{r}^{\prime}\right)=\frac{e^{-{\mathrm{i}}k_{0}\lvert\bm{r}-\bm{r}^{\prime}\rvert}}{4\pi\lvert\bm{r}-\bm{r}^{\prime}\rvert},$
(3)

where $k_{0}$ is the free-space wavenumber, $\bm{r}$ the source point, and
$\bm{r}^{\prime}$ the observation point. Each element of the VSIE coupling
matrix is a 5D integral formed from the inner product between the discretized
${\cal{N}}$ operator applied on a VIE basis function, and an RWG basis
function.

####  II-A2 VSIE implementation of GMT

GMT estimates tissue electrical properties from MR measurements by solving an
inverse problem [15]. In GMT, the cost function compares actual measurements
against simulated measurements of the relative $b_{1}^{+}$ fields generated by
multiple sources (e.g., multiport transmit coils) inside a sample and
iteratively updates the estimate of the sample’s electrical properties. GMT
was initially demonstrated using the JVIE formulation for the solutions of the
forward problem, therefore ignoring the effect of the dielectric sample on the
incident fields. However, these interactions must be taken into account for
accurate in-vivo experiments with close-fitting RF coils. In other words, the
GMT framework must be ported from a VIE to a VSIE formulation, in which the
incident fields are not constant but calculated at each GMT iteration as

$\displaystyle E^{\rm inc}\left(\mathbf{\epsilon_{r}}\right)$
$\displaystyle=T_{\rm bc}J_{\rm s}\left(\mathbf{\epsilon_{r}}\right),$ (4)
$\displaystyle H^{\rm inc}\left(\mathbf{\epsilon_{r}}\right)$
$\displaystyle=Z_{\rm bc}J_{\rm s}\left(\mathbf{\epsilon_{r}}\right).$

where $E^{\rm inc}$ and $H^{\rm inc}$ are the discretized incident electric
and magnetic fields respectively. $\epsilon_{r}$ is the complex permittivity.
$Z_{\rm bc}$ maps the equivalent surface electric currents to the magnetic
fields with the aid of the ${\cal{K}}$ operator:

${\cal{K}}\left(\bm{s}\right)\triangleq\nabla\times\int\limits_{\Omega}g\left(\bm{r}-\bm{r}^{\prime}\right)\bm{s}\left(\bm{r}^{\prime}\right)d^{3}\bm{r}^{\prime}.$
(5)

In addition, in the new implementation, the gradient of the GMT’s cost
function will require to solve a Hermitian adjoint system of equations that
includes multiplications with the conjugate transpose of $Z_{\rm bc}$.

Matrix-vector products involving the coupling matrix are typically performed
without storing the full matrix, due to its intractably large size. In the
case of iterative inverse problem solutions, such as in GMT, this approach
could considerably increase the computation time, because it requires the re-
assembly of the full matrix at each iteration. In the next sections, we
propose a compression algorithm that reduces the computational burden, by
enabling one to assembly the full coupling matrix only once and then perform
just the matrix-vector multiplications in each of GMT’s iterations.

###  II-B Numerical Linear Algebra Methods

####  II-B1 Tucker Decomposition

A 3D tensor ${\mathcal{A}}\in\mathbb{C}^{n_{1}\times n_{2}\times n_{3}}$ can
be decomposed with the Tucker model [17] in the following form:

$\displaystyle{\mathcal{A}}$
$\displaystyle={\mathcal{G}}\times_{1}U^{1}\times_{2}U^{2}\times_{3}U^{3},\>\text{or}$
(6) $\displaystyle{\mathcal{A}}_{ijk}$
$\displaystyle=\sum_{\chi=1}^{r_{1}}\sum_{\psi=1}^{r_{2}}\sum_{\zeta=1}^{r_{3}}{\mathcal{G}}_{\chi\psi\zeta}U^{1}_{i\chi}U^{2}_{j\psi}U^{3}_{k\zeta}.$

Here $U^{\gamma}\in\mathbb{C}^{n_{\gamma}\times r_{\gamma}},\gamma=1,2,3$, are
unitary matrices, dubbed as Tucker factors, while
${\mathcal{G}}\in\mathbb{C}^{r_{1}\times r_{2}\times r_{3}}$ is the Tucker
core. The dimensions of the Tucker core indicate the multilinear (or Tucker)
ranks of ${\mathcal{A}}$. The symbols $\times_{\gamma}$ are called $n$-mode
products and perform a convolution over the $\times_{\gamma}$ axis, for
example, the $\times_{1}$ product performs the following operation:

${\mathcal{P}}={\mathcal{G}}\times_{1}U^{1}\Leftrightarrow{\mathcal{P}}_{i\psi\zeta}=\sum\limits_{\chi=1}^{r_{1}}{\mathcal{G}}_{\chi\psi\zeta}U^{1}_{i\chi}.$
(7)

Here, ${\mathcal{P}}\in\mathbb{C}^{n_{1}\times r_{2}\times r_{3}}$. The
expansion of ${\mathcal{A}}$ in equation (6) can be truncated to a desired
tolerance, and return an approximation of ${\mathcal{A}}$. A visual
representation of Tucker decomposition can be seen in Fig. 1.

Figure 1: Visual representation of Tucker decomposition.

To compute the above-mentioned Tucker components, one has to choose a suitable
compression algorithm. The higher order singular value decomposition (HOSVD)
is an orthogonal Tucker decomposition, widely used because it has a proven
upper error bound [39]. Moreover, the algorithm is based entirely on singular
value decomposition (SVD), which provides a robust and stable approximation of
the initial tensor. Note that SVD requires the assembly of the initial array,
which could be challenging for large tensors. In such cases, one could
implement the Tucker decomposition using a 3D cross approximation algorithm
[40].

####  II-B2 Cross Approximation

A matrix $A\in\mathbb{C}^{n_{1}\times n_{2}}$ can be approximated with the so-
called 2D cross approximation method [25, 26] as follows:

$A\approx UV^{*}.$ (8)

Here, $U\in\mathbb{C}^{n_{1}\times r_{c}}$ and $V\in\mathbb{C}^{n_{2}\times
r_{c}}$. $r_{c}$ represents the column rank of matrix $A$. Cross approximation
algorithms construct the decomposition of $A$ by using only some rows and
columns of it, differently than SVD, which depends on the availability of the
full matrix. Several algorithms have been developed over the previous decades
for the implementation of cross approximation, with two being the most used
ones: the ACA [28, 29], and the maximum volume-based cross algorithm [27, 41].
The latter requires the implementation of LU, QR and SVD for its efficient
implementation. Therefore, the memory demand of the algorithm increases
drastically for large tall matrices, such as the coupling matrices $T_{\rm
bc}$ and $Z_{\rm bc}$ in the case fine voxel resolutions. On the other hand,
the memory demand of ACA is only dictated by the size of matrices $U$ and $V$.

##  III Tucker-Based Compression Algorithm

In VSIE, the columns of the coupling matrix describe interactions between coil
and body basis functions through the Green’s functions of equations (2) and
(5); therefore, they represent well-separated geometrical blocks. Due to the
3D nature of the problem, the key idea of our proposed compression algorithm
is to reshape these columns as tensors and approximate them with the low
multilinear Tucker model. This compression strategy enabled us to develop a
new method to efficiently perform the matrix-vector product and an extension
to ACA, which are described later in this section.

###  III-A Matrix Assembly

Each of the $m$ columns of the coupling matrices $Z_{\rm bc}$ and $T_{\rm bc}$
can be seen as the concatenation of $q$ vectors, where each vector represents
the component-wise interaction between one RWG element on the coil and the
basis functions of all the voxels in the body domain. For PWC, $q=3$, whereas
for PWL, $q=12$. Since these vectors model interactions between remote
discretization elements, they have low-rank properties [42]. To exploit the
low-rank, each column of the coupling matrix can be reshaped as $q$ 3D tensors
${\mathcal{Z}}^{kj}\in\mathbb{C}^{n_{1}\times n_{2}\times n_{3}}$, $k=1:q$,
$n_{\rm v}=n_{1}\times n_{2}\times n_{3}$, which are compressible with the
Tucker decomposition [43]. A graphical description of the algorithm is shown
in Fig. 2 for $Z_{\rm bc}$ and PWC basis functions,.

Figure 2: Visual representation of the Tucker-based algorithm for the
compression of the $Z_{\rm bc}$ matrix, in the case of PWC basis functions.
Each vector can be reshaped into a 3D tensor that is then compressed via
Tucker decomposition.

If the coupling matrix is first approximated with the ACA as $UV^{*}$, then
our approach can still be used to compress the $r_{c}$ columns of $U$. In
fact, cross approximation is a well-conditioned operation, therefore the
Tucker ranks of the reshaped columns of $U$ will be similar to the ones of the
reshaped columns of the coupling matrix. The $V$ matrix here is usually much
smaller than $U$ and does not require any additional compression.

TABLE II shows the total memory footprint associated with the assembly of the
coupling matrix: Full assembly, assembly with ACA, and assembly with our
proposed method by compressing either the columns of the coupling matrix
(Tucker) or the columns of $U$ (ACA+Tucker). The memory required after
compressing the coupling matrix with Tucker is $n_{\rm
v}/\left(r^{3}+3nr\right)$ times smaller than the memory required by the full
matrix, where $n$ and $r$ refer to the tensor’s linear dimension and Tucker
rank, respectively. If our Tucker-based compression method is instead applied
after the ACA assembly, then the total compression improves by a factor of
$\sim m/r_{c}$, given that $r_{c}$ is small. TABLE II also shows the
computational complexity of the assembly operations. The multiplicative
constant factor $c_{1}$, which is present in all cases, represents the cost to
compute the elements of the coupling matrix and is usually large. In fact,
each element requires a 5D integration, whose computational cost depends on
the number of both surface and volume quadrature integration points. As a
result, the assembly of the matrix is extremely inefficient and should be
implemented in parallel for multiple voxel-RWG basis function interactions.

Note that in certain cases, for example when the coil is close to the body,
ACA may not achieve a meaningful compression and would not be advantageous to
combine it with Tucker decomposition. In such cases, the preferable approach
would be to divide the coupling matrix in $q$ blocks of size $n_{\rm v}\times
m$, assembled them in parallel, and then compress their tensor components with
a Tucker-based method like HOSVD. Alternatively, if the coupling matrix is
sufficiently small, one could assemble it in its full form and then apply
Tucker directly to compress it.

TABLE II: Complexity for Constructing the Coupling Matrix Assembly Method | Operations | Memory
---|---|---
Full | ${\mathcal{O}}\left(c_{1}qn_{\rm v}m\right)$ | $qn_{\rm v}m$
ACA | ${\mathcal{O}}\left(c_{1}r_{c}^{2}\left(qn_{\rm v}+m\right)\right)$ | $qn_{\rm v}r_{c}+mr_{c}$
Tucker | Full + ${\mathcal{O}}\left(rqn_{\rm v}m\right)$ | $q\left(r^{3}+3nr\right)m$
ACA+Tucker | ACA + ${\mathcal{O}}\left(rqn_{\rm v}r_{c}\right)$ | $q\left(r^{3}+3nr\right)r_{c}+mr_{c}$

###  III-B Matrix-Vector Product

Decompressing the coupling matrix to compute the matrix-vector product
$\mathbf{y}=Z_{\rm bc}\mathbf{x}$, like in equations (4), may not be possible
due to computer or GPU memory limitations. To address this, we propose a novel
approach to efficiently compute the matrix-vector product without fully
decompressing the coupling matrix. We initiate $\mathbf{y}$ as a vector of
zeros. Inside a loop that cycles over the RWG basis functions, we decompress
the $q$ tensors of a column $j\in[1,m]$, reshape them as vectors, and
concatenate them to form the $j$-th column of the original $Z_{\rm bc}$
matrix. The vector-scalar product between the $j$-th column and
$\mathbf{x}_{j}$ is then computed, and the result increments the elements of
$\mathbf{y}$. The same algorithm can be followed for the matrix-matrix product
$Y=Z_{\rm bc}X$.

The conjugate transpose matrix-vector product $\mathbf{y}=Z^{*}_{\rm
bc}\mathbf{x}$ is required for the computation of the gradient of the cost
function in the VSIE-based GMT implementation. This case is slightly different
than the standard matrix-vector product: inside the loop cycling through the
RWG functions, a row-column vector product must be computed between the
conjugate transpose of the decompressed $j$-th column of $Z_{\rm bc}$ and
$\mathbf{x}$, which yields the scalar $\mathbf{y}_{j}$. The algorithm remains
instead the same for the conjugate transpose matrix-matrix products. Both
algorithms (for $p$-columned matrices $X$ and $Y$) are summarized in the
pseudocode below:

Algorithm 1: $\>Y=Z_{\rm bc}X$

1:for k=1:$q$ do

2: $Y^{k}=0$

3:for j=1:$m$ do

4: for k=1:$q$ do

5: $\text{Decompress}\>{\mathcal{Z}}^{kj}$

6: $Y^{k}+={\mathcal{Z}}^{kj}(:)\mathbf{X_{j:}}$

7:$Y=\begin{bmatrix}Y^{1}\\\ \cdots\\\ Y^{q}\end{bmatrix}$

Algorithm 2: $\>Y=Z^{*}_{\rm bc}X$

1:$Y=0$

2:for j=1:$m$ do

3: for k=1:$q$ do

4: $\text{Decompress}\>{\mathcal{Z}}^{kj}$

5: $Y_{j:}=\begin{bmatrix}{\mathcal{Z}}^{1j}(:)\\\ \cdots\\\
{\mathcal{Z}}^{qj}(:)\end{bmatrix}^{*}X$

In Algorithm 1, $X\in\mathbb{C}^{m\times p}$, $Y\in\mathbb{C}^{qn_{\rm
v}\times p}$ (vice-versa for Algorithm 2), and ${\mathcal{Z}}^{kj}(:)$ is the
reshaped column vector of the tensor component ${\mathcal{Z}}^{kj}$. The
algorithms remain the same if $Z_{\rm bc}$ is compressed with ACA first
($Z_{\rm bc}=UV^{*}$). One has to replace $Z_{\rm bc}$ with $U$, $m$ with
$r_{c}$, and assign $X=V^{*}X$ for Algorithm 1, and $Y=VY$ for Algorithm 2.
Both the standard and the conjugate transpose matrix-vector products have the
same complexity, shown in TABLE III, for the full, ACA, Tucker, and ACA+Tucker
compressed cases. The full matrix-vector product is faster than the Tucker-
compressed approach by a factor of $(r+p)/p$, which depends on the number of
columns of $X$ and $Y$ and the Tucker rank. ACA can be faster than the full
case for small values of $r_{c}$. Although the approach based on Tucker
decomposition is slower because it requires additional flops compared to the
other methods, it is more likely to fit in GPUs, thanks to its small memory
footprint.

TABLE III: Matrix-Vector Product Complexity Matrix Form | Operations Complexity
---|---
Full | ${\mathcal{O}}\left(qn_{\rm v}mp\right)$
ACA |  ${\mathcal{O}}\left(qn_{\rm v}r_{c}p\right)$ \+ ${\mathcal{O}}\left(r_{c}mp\right)$
Tucker |  ${\mathcal{O}}\left(rqn_{\rm v}m\right)$ \+ ${\mathcal{O}}\left(qn_{\rm v}mp\right)$
ACA+Tucker |  ${\mathcal{O}}\left(rqn_{\rm v}r_{c}\right)$ \+ ${\mathcal{O}}\left(qn_{\rm v}r_{c}p\right)$ \+ ${\mathcal{O}}\left(r_{c}mp\right)$

###  III-C Tucker-based ACA

If the coupling matrix is first compressed with ACA, the previous methods for
matrix assembly and matrix-vector product could still be applied to the matrix
$U$ of the cross approximation. However, for the case of realistic body models
discretized with fine voxel resolutions, the traditional implementation of ACA
(a detailed description can be found in [44]) might fail due to random access
memory (RAM) overflow because of the size of $U$ (see section IV.B.2). To
address this, we propose an extension of ACA in which the matrix $U$ is
assembled directly in a compressed form, based on our proposed Tucker
decomposition technique. The algorithm is summarized in pseudocode bellow:

Algorithm 3: ACA of $Z\in\mathbb{C}^{m_{1}\times m_{2}}$. Assembly with
compressed $U$

1:$\epsilon=\text{tolerance}$, $i=1$, $\mathbf{s}_{1}=0$

2:$\mathscr{U}=[]$, $V=[]$

3:for $k=1:\text{min}(m_{1},m_{2})$ do

4: $\mathbf{r}\leftarrow Z_{i:}$

5: if $k>1$ then

6: $[f_{1},f_{2},f_{3},p]\leftarrow i$

7: for $l=1:size(\mathscr{U},2)$ do

8: $[{\mathcal{G}},U^{1},U^{2},U^{3}]\leftarrow U_{pl}$

9:
$\mathbf{t}_{l}={\mathcal{G}}\times_{1}U_{f_{1},:}^{1}\times_{2}U_{f_{2},:}^{2}\times_{3}U_{f_{3},:}^{3}$

10: $\mathbf{r}\mathrel{+}=-\mathbf{t}V^{*}$

11: $j\leftarrow\text{index of max element of}\>\mathbf{r}$

12: $\mathbf{y}\leftarrow\left(\mathbf{r}/\mathbf{r}_{j}\right)^{*}$

13: $\mathbf{x}\leftarrow Z_{:j}$

14: if $k>1$ then

15: $\mathbf{x}\mathrel{+}=-\textbf{Alg1}(\mathscr{U},V_{j:}^{*})$

16:
$\mathbf{s}_{k+1}\leftarrow\mathbf{s}_{k}+(\left\lVert\mathbf{x}\right\rVert\left\lVert\mathbf{y}\right\rVert)^{2}$

17: if $k>1$ then

18:
$\mathbf{s}_{k+1}\mathrel{+}=2\sum\left[\mathrm{Re}\,\\{\left(\textbf{Alg2}(\mathscr{U},\mathbf{x})\right)\odot\left(V^{*}\mathbf{y}\right)\\}\right]$

19: Reshape $\mathbf{x}$ to $q$ ${\mathcal{X}}^{q}$ tensors.

20:
$\mathscr{U}=[\mathscr{U}\>\text{HOSVD}({\mathcal{X}}^{1},3\epsilon)\>\cdots\>\text{HOSVD}({\mathcal{X}}^{q},3\epsilon)]$

21: $V=[V\>y]$

22: if
$\left\lVert\mathbf{x}\right\rVert\left\lVert\mathbf{y}\right\rVert\leq\epsilon\sqrt{\mathbf{s}_{k+1}}$
then break

23: $\mathbf{x}=\lvert\mathbf{x}\rvert$, $\mathbf{x}_{i}=0$

24: $i\leftarrow\text{index of max element of}\>\mathbf{x}$

Here $\mathscr{U}$ is a struct of size $q\times r_{c}$ ($q=3$ for PWC, $q=12$
for PWL, $r_{c}$ is the rank of $Z$) than contains tensors. Each time a new
column of $U$ is computed, it is reshaped to $q$ tensors, which are then
compressed with a truncated HOSVD of tolerance $3\epsilon$ (line 19,20). The
HOSVD tolerance has to be higher than the ACA tolerance since the irrelevant
numerical digits ($<1e-3$) appearing in $U$ are incompressible. We found that
a $3$ times higher tolerance is a good choice for our numerical examples. To
perform matrix- and conjugate transpose matrix-vector products with the
compressed $U$ we followed Algorithms 1 (line 15) and 2 (line 18). Finally,
when a row of $U$ is requested in Algorithm 3, we first calculate the voxel
and basis function component corresponding to that row (line 6) and then
decompress, using equation (6), only the required elements from the Tucker
compressed components of $U$ (line 8,9). The proposed algorithm avoids RAM
overflowing, but it is slower than the traditional ACA due to the multiple
tensor decompressions. Nevertheless, it could always be accelerated via GPU,
since its memory demand is as low as for one column of the coupling matrix.

##  IV Numerical Experiments

###  IV-A Tucker Rank Behavior

In this section, we study the low-Tucker rank properties of the $Z_{\rm bc}$
coupling matrix. We considered multiple geometrical scenarios and altered the
distance between the conductive surface (coil) and the VIE domain, the
operating frequency and the conductive surface’s discretization. The tensor
components of the columns of the coupling matrix were compressed with the
HOSVD algorithm and a tolerance of $1e-8$, which yielded a relative error
similar to the tolerance for all cases, due to the robustness of the SVD
itself. Such error can be considered negligible, because the tolerance of the
iterative solver used in FFT-based VIE systems is usually orders of magnitude
higher.

####  IV-A1 Tucker Rank vs. Distance

It is well established that the Green’s function integro-differential
operators between well-separated geometries present low-rank properties [42].
Here we studied the relation between the low multilinear rank of the
compressed coupling matrix $Z_{\rm bc}$ and the distance between the body’s
domain and the coil. We set the frequency to $298.06$ MHz, the operating
frequency of $7$ Tesla MRI. We modeled a single perfectly electric conducting
(PEC) loop coil of radius $\rho=0.50$ m and discretized it with $125$
triangular elements. The coil was centered in $\left(0,d,0\right)$, where $d$
was varied as $0.55$, $0.6$, $\dots$, $1$ m. The domain was a cuboid with edge
length of $1$ m, centered at $(0,0,0)$ and discretized with voxels of $1$ cm
isotropic resolution and PWC basis functions (Fig. 3). As a result, the
tensor’s size was $101\times 101\times 101$ and the memory required by the
fully assembled $Z_{\rm bc}$ was $5848$ MBs.

Figure 3: Loop-cubic domain geometry. The loop coil was shifted on the
$\hat{y}$ direction, for $10$ discrete distances between $0.55$ to $1$ m from
the center of the cube.

Fig. 4 illustrates the reduction of the maximum rank (maximum of all Tucker
ranks for all components) of the coupling matrix (right axis), and the total
memory of the compressed matrix using the algorithm described in Section IV
(left axis). It is evident that the greater the distance between the domain
and the coil, the lower the multilinear ranks and the required memory. The
compression factor varied between $\sim 50$ and $190$, depending on the
distance.

Figure 4: Memory footprint (left) and maximum rank (rank) of the compressed
$Z_{\rm bc}$ matrix, for different distances between the loop and the cubic
domain.

####  IV-A2 Tucker Rank vs. Frequency

The work presented in [42] showed that the rank of the integral operators for
3D problems increases linearly with respect to the operating frequency. This
was confirmed in [10], for the BTTB defining tensors of the FFT-based VIE
systems (discretized integral operators). These tensors were columns of the
corresponding Galerkin MoM matrices and modeled the interactions between one
voxel’s basis function and the whole domain via the ${\cal{N}}$ or ${\cal{K}}$
operators. In the present study, the tensors are columns of the coupling
matrix and model the interactions between one RWG basis function and the whole
body domain via the same operators. Since in both cases the interactions
between separated geometry blocks are modeled, one can expect a similar
behavior for the Tucker ranks.

To confirm this, we performed a frequency sweep ($300$, $600$, $\dots$,
$2700$MHz) for the setup in Fig. 3. The coil was discretized with $125$
elements, whereas the voxel’s isotropic resolution was set to $\lambda/20$,
with $\lambda$ being the wavelength. We repeated the calculations for three
positions of the coil ($d=0.55$, $0.65$, and $0.8$ m). The memory footprint
(left) of the compressed matrix, along with the maximum rank (right), are
shown in the dual axis chart of Fig. 5. The memory footprint increased
linearly with frequency, whereas the maximum rank grew at different rates for
the three investigated cases. This is expected because the maximum rank
represents the worst-case scenario among all ranks, whereas the memory
footprint summarizes the overall effect of all ranks.

Figure 5: Memory footprint (left) and maximum rank (rank) of the compressed
$Z_{\rm bc}$ matrix, for different operating frequencies. Results are shown
for three different distances between the loop and the domain.

####  IV-A3 Tucker Rank vs. Surface Mesh

Let us consider a fixed mesh for the domain and refine only the surface mesh
of the coil. As the coil mesh is refined, the edges of the triangular elements
become smaller and the number of columns of the coupling matrix increases,
with each column representing more remote element interactions between an edge
of the triangular mesh and the voxelized domain. As a result, we should expect
a drop in the Tucker ranks. To verify this, we used the same domain of the
previous sections, and a PEC equilateral triangle with centroid at
$(0,0.55\text{m},0)$ and one vertex at $(0,0.55\text{m},0.5\text{m})$. The
triangle was discretized with $10$ different meshes, starting from triangular
element’s edge of $0.5$m and reducing it by a factor of $\sqrt{2}$, which
resulted in $4$, $6$, $11$, $30$, $48$, $102$, $184$, $358$, $727$, and $1480$
elements. Fig. 6 reports the maximum rank as a function of the length of the
triangular element’s edge, confirming that the rank is smaller when the PEC
triangle’s mesh is finer.

Figure 6: Maximum rank of the compressed $Z_{\rm bc}$ matrix, for various PEC
triangle’s meshes. The rank drops as we refine the mesh.

###  IV-B Application to VSIE-based MRI Simulations

Here we aim to validate the performance of the proposed algorithms for the
assembly of the coupling matrix $Z_{\rm bc}$, and the matrix-vector
implementation for two VSIE-based MRI applications. Both numerical experiments
were implemented in Matlab, except for the matrix assembly part which was
written in C++. For the GPU computations, we used an NVIDIA Quadro Volta GV100
32GB HBM2 PCIe. For the CPU computations, in the first experiment we used a
server with CentOS 6.9 operating system and an Intel(R) Xeon(R) CPU E5-2699 v3
at 2.30GHz, while for the second experiment we used a server with Ubuntu
18.04.5 LTS operating system and an Intel(R) Xeon(R) Gold 6248 CPU at 2.50GHz.
We parallelized on $12$ workers where needed.

####  IV-B1 Head Coil Experiments

We first demonstrated the proposed compression method for an 8-ports close-
fitting head coil, previously designed for GMT [16], which we loaded with the
“Billie” realistic head model from the virtual family population [45] (Fig.
7). The operating frequency was set to $298$ MHz.

Figure 7: Coil-head geometry. The RF coil (discretized with $2380$ triangular
element edges) was loaded with the voxelized realistic human head model
“Billie” (discretized with voxels of $1$ mm isotropic resolution).

The VSIE-based implementation of GMT requires performing operations on the
coupling matrix $Z_{\rm bc}$ and its conjugate transpose. We analyzed the
memory footprint reduction for the compressed coupling matrix and measured the
computation time for both the matrix- and conjugate transpose matrix-vector
products using the algorithms presented in section III. The coil was
discretized with both a coarse ($516$ RWG) and a fine ($2380$ RWG) mesh
resolution. For the VIE domain enclosing the head, we tested three different
voxel resolutions, namely $5$, $2$, and $1$ mm3, which resulted in $34\times
38\times 45$, $84\times 96\times 116$, and $168\times 188\times 222$ voxels,
respectively. Both PWC ($3$ unknowns per voxel) and PWL ($12$ unknown per
voxel) VIE basis functions were considered.

We used a tolerance of $1e-8$ for HOSVD, which would ensure accurate
estimation of electrical properties in an actual GMT experiment. Since the
coil closely fits the head, ACA (or SVD-based methods in general) are expected
to provide negligible compression with a tolerance of $1e-8$. We confirmed
this for the case of PWC basis functions, $5$ mm3 voxel resolution, and fine
coil discretization, for which, in fact, we found that $2238$ of the $2380$
singular values would be needed to accurately represent $Z_{\rm bc}$,
compressing the matrix from $6.18$ GB to $6.07$ GB. Consequently, for the head
coil experiments we did not use the Tucker-based ACA algorithm, but instead we
compressed the columns of the coupling matrix only with the HOSVD-based
method.

##### Memory Compression

The memory footprint for the assembly of the coupling matrix $Z_{\rm bc}$ is
shown in TABLE IV. The memory required to assemble the full matrix was
considerably larger than for the HOSVD-compressed matrix. For example, for PWC
basis functions, voxel resolution of $1$ mm3, and fine coil mesh, the required
memory in the full matrix case was $>740$ GBs, whereas the compressed matrix
required only $2.6$ GBs. Note that in the challenging case of PWL basis
functions, $1$ mm3 voxel resolution, and fine coil mesh, it was not feasible
to apply our compression method. In fact, the memory requirements even for
just one of the $q$ blocks of the matrix (see Section IV. A.) were
prohibitively large for our server. While we could have still implemented the
HOSVD compression by further dividing the matrix in smaller blocks, that would
have required $\sim 1$ month of computations. An alternative method for such
costly cases is mentioned in the discussion and will be pursued in future
work.

TABLE IV: Memory Requirements (GBs) of $Z_{\rm bc}$ Voxel Res. | Assembly | PWC-coarse | PWC-fine | PWL-coarse | PWL-fine
---|---|---|---|---|---
$5$ mm3 | Full | 1.34 | 6.18 | 5.36 | 24.74
HOSVD | 0.20 | 0.88 | 0.86 | 3.77
$2$ mm3 | Full | 21.57 | 99.52 | 86.30 | 398.09
HOSVD | 0.40 | 1.74 | 1.79 | 7.49
$1$ mm3 | Full | 161.73 | 745.99 | 646.95 | 2983.99
HOSVD | 0.63 | 2.62 | 2.85 | N/A

Fig. 8 shows that the compression factor, defined as the memory of the full
matrix over the memory of the compressed one, decreased as the voxel
resolutions ($h$) of the VIE domain’s grid became coarser. The behavior of the
compression factor was similar for PWC or PWL basis functions, either with
fine or coarse mesh discretization. This confirms the excellent stability of
our compression method.

Figure 8: Compression factor of the compressed matrix $Z_{\rm bc}$. Results
are shown for all investigated head and coil discretizations.

Fig. 9 shows the maximum Tucker rank, obtained with HOSVD, for all tensor
components of the coupling matrix. The rank decreased slowlier than the
compression factor (Fig. 8) for coarser discretizations of the VIE domain. For
example, when PWC basis functions and fine coil resolution were used (PWC-
fine), the maximum rank decreased by only $1.5$ times (from $42$ to $28$) when
the isotropic voxel size increased from $1$ to $5$ mm3, which corresponds to a
$5$ times smaller grid in all directions. For all cases, the maximum rank was
smaller for finer coil meshes, which is in agreement with the results shown in
section V.A.3.

Figure 9: Maximum rank as a function of voxel resolution of the VIE domain.
The maximum rank was calculated among all Tucker ranks of the decomposed
tensors, which were obtained with HOSVD. Results are shown for all
investigated head and coil discretizations. For each voxel resolution, the
size of the corresponding VIE discretization grid is indicated.

##### Computation Time

TABLE V reports the computation time for the assembly of the full and the
HOSVD-compressed coupling matrix (rounded to the nearest higher second). For
PWC, we used one quadrature integration point per triangle and voxel, while
for PWL, two for each triangle and eight for each voxel. For the low
resolution matrices, the assembly time for the compressed matrix was larger
than the one for the full matrix by a factor $\leq q$, since in the HOSVD case
the compressed matrix was assembled as $q$ sequential compressed matrix
blocks. For the larger matrices, our server could not perform the assembly of
the full matrix, due to the prohibitively large memory requirements, but it
was able to assemble the compressed matrix using our proposed method.

TABLE V: Time Footprint (hh:mm:ss) of $Z_{\rm bc}$ Assembly Voxel Res. | Assembly | PWC-coarse | PWC-fine | PWL-coarse | PWL-fine
---|---|---|---|---|---
$5$ mm3 | Full | 00:00:12 | 00:00:45 | 00:03:31 | 00:14:17
HOSVD | 00:00:34 | 00:02:33 | 00:32:27 | 02:21:21
$2$ mm3 | Full | 00:02:48 | N/A | 01:09:49 | N/A
HOSVD | 00:06:16 | 00:33:28 | 07:15:45 | 30:32:42
$1$ mm3 | Full | N/A | N/A | N/A | N/A
HOSVD | 00:27:59 | 02:31:52 | 52:17:11 | N/A

TABLE VI and VII summarize the computation times for the matrix- and conjugate
transpose matrix-vector products. Compared to the full form case (Full-CPU),
the compressed matrix-vector product requires additional operations for the
decompression of the tensors. While Algorithms 1 and 2 can reduce the memory
requirements of the matrix-vector products, the time footprint varies based on
how these algorithms are implemented. In particular, the nested loops over the
$m$ RWG functions can be either parallelized on a CPU, if the RAM can support
multiple tensor decompressions in parallel (HOSVD-CPU), or performed
sequentially using a GPU (HOSVD-GPU). In our tests, we multiplied $Z_{\rm bc}$
with $X\in\mathbb{C}^{m\times 8}$ (TABLE VI) and $Z^{*}_{\rm bc}$ with
$\Phi\in\mathbb{C}^{qn_{\rm v}\times 8}$ (TABLE VII), where both $X$ and
$\Phi$ were random matrices to keep the results general. The eight columns of
$X$ could correspond, for example, to the currents associated with the eight
channels of the coil in Fig. 7.

TABLE VI: Time Footprint (hh:mm:ss) of $Y=Z_{\rm bc}X$ Voxel Res. | Form | PWC-coarse | PWC-fine | PWL-coarse | PWL-fine
---|---|---|---|---|---
$5$ mm3 | Full-CPU | 00:00:01 | 00:00:02 | 00:00:02 | 00:00:06
HOSVD-CPU | 00:00:03 | 00:00:07 | 00:00:07 | 00:00:24
HOSVD-GPU | 00:00:03 | 00:00:11 | 00:00:10 | 00:00:44
$2$ mm3 | Full-CPU | 00:00:06 | N/A | 00:00:25 | N/A
HOSVD-CPU | 00:00:25 | 00:01:46 | 00:01:37 | 00:06:24
HOSVD-GPU | 00:00:04 | 00:00:14 | 00:00:13 | 00:00:56
$1$ mm3 | Full-CPU | N/A | N/A | N/A | N/A
HOSVD-CPU | 00:02:54 | 00:11:25 | 00:11:30 | N/A
HOSVD-GPU | 00:00:13 | 00:00:53 | 00:00:52 | N/A
TABLE VII: Time Footprint (hh:mm:ss) of $\Psi=Z^{*}_{\rm bc}\Phi$ Voxel Res. | Form | PWC-coarse | PWC-fine | PWL-coarse | PWL-fine
---|---|---|---|---|---
$5$ mm3 | Full-CPU | 00:00:01 | 00:00:02 | 00:00:02 | 00:00:06
HOSVD-CPU | 00:00:02 | 00:00:04 | 00:00:04 | 00:00:26
HOSVD-GPU | 00:00:02 | 00:00:07 | 00:00:05 | 00:00:21
$2$ mm3 | Full-CPU | 00:00:06 | N/A | 00:00:22 | N/A
HOSVD-CPU | 00:00:13 | 00:00:45 | 00:00:49 | 00:03:25
HOSVD-GPU | 00:00:03 | 00:00:11 | 00:00:10 | 00:00:42
$1$ mm3 | Full-CPU | N/A | N/A | N/A | N/A
HOSVD-CPU | 00:01:20 | 00:05:19 | 00:05:23 | N/A
HOSVD-GPU | 00:00:11 | 00:00:45 | 00:00:41 | N/A

For $5$ mm3 isotropic voxel resolution, the Full-CPU matrix-vector product was
the fastest for all cases, because the coupling matrix is small. For $2$ and
$1$ mm3 voxel resolution, the HOSVD-GPU implementation was the fastest. Note
that the Full-CPU case could not be performed for high voxel and coil mesh
resolutions, due to the excessive memory requirements. The HOSVD-CPU was
slower than HOSVD-GPU, except for the $5$ mm3 voxel resolution.

####  IV-B2 Body Coil Experiments

For the second MRI experiment, we simulated the volume bodycoil of a
commercial 3T MRI scanner [46, 47] and we loaded it with “Billie”, from the
virtual family population [45] (Fig. 10). The frequency was set to $123$ MHz,
corresponding to $3$ Tesla MRI. The coil has $32$ legs, a radius of $35.5$ cm,
length of $45$ cm, and is centered at $\left(0,0,0\right)$. We also modeled
the system conductive shield, which has a radius of $37.2$ cm, a length of
$1.5$ m and is centered at $\left(0,0,0\right)$. The distance between the coil
and the cuboid VIE domain enclosing “Billie” was $15.5$ cm and $24.5$ cm in
the $x$ and $y$, respectively. In contrast with the previous case where the
coil tightly fitted the head, here the coil is remote enough to allow a good
compression of the coupling matrix $Z_{\rm bc}$ with ACA. For this experiment,
we used PWC and PWL basis functions and three voxel resolutions ($5$, $2$, and
$1$ mm3), which corresponded to $81\times 44\times 108$, $205\times 109\times
270$, and $409\times 219\times 541$ voxels for the VIE domain. For the coil
and the shield we used $9450$ RWG basis functions. Two quadrature integration
points were used for each triangle and eight for each voxel, both for PWC and
PWL basis functions.

Figure 10: Coil-body geometry. The RF coil and shield (discretized with $9450$
triangular element edges) was loaded with part of the voxelized realistic
human body model “Billie” (discretized with voxels of $2$ mm isotropic
resolution).

##### Matrix Assembly

TABLE VIII, summarizes the memory requirements and the assembly time for the
coupling matrix $Z_{\rm bc}$. The ACA tolerance was set to $1e-3$ to achieve
good compression. The ACA rank of $Z_{\rm bc}$ was $250$ for the $5$ mm3 cases
and $287$ for the $2$ and $1$ mm3 cases. The maximum Tucker rank of $U$ was
between $15$ and $18$ for all cases. In the $5$ mm3 case, our results show
that ACA could offer an excellent compression of the coupling matrix and the
assembly could be rapidly performed in CPU. For $2$ mm3, ACA’s memory
requirements were large and ACA was outperformed in speed by our proposed
Algorithm 3 (Tucker-based ACA), for which the low memory footprint allowed
using a GPU. For $1$ mm3 resolution, the standard ACA algorithm could not be
performed even on a server equipped with hundreds of GB’s of RAM, due to
overwhelming memory requirements. On the other hand, our proposed ACA
extension in Algorithm 3 kept the memory demand small, enabling for fast
matrix assembly in GPU. Note that the full matrix assembly was only possible
for $5$ mm3 voxel resolution and PWC basis functions.

TABLE VIII: Memory Demand (GB) and Time Footprint (hh:mm:ss) of $Z_{\rm bc}$ Assembly Voxel Res. | Form | PWC | PWL
---|---|---|---
memory | time | memory | time
$5$ mm3 | Full-CPU | 162.60 | 00:15:49 | 650.42 | N/A
ACA-CPU | 4.33 | 00:01:53 | 17.24 | 00:06:39
Algorithm 3-GPU | 0.036 | 00:03:26 | 0.036 | 00:09:58
$2$ mm3 | Full-CPU | 2548 | N/A | 10195 | N/A
ACA-CPU | 77.44 | 00:37:16 | 309.65 | 04:30:36
Algorithm 3-GPU | 0.041 | 00:26:57 | 0.042 | 01:28:53
$1$ mm3 | Full-CPU | 20471 | N/A | 81884 | N/A
ACA-CPU | 621.75 | N/A | 2486 | N/A
Algorithm 3-GPU | 0.041 | 03:35:36 | 0.042 | 15:20:41

For the coarser case of $5$ mm3 voxel resolution and PWC basis functions, the
time footprint of Algorithm 3 for CPU (not shown in the table) was 00:17:50,
which is $\sim 5$ times slower than for the GPU execution.

##### Matrix-Vector Product Performance

The time footprints for the matrix-vector product between the compressed
coupling matrix $Z_{\rm bc}$ and a random vector
$\mathbf{x}\in\mathbb{C}^{m\times 1}$ are shown in TABLE IX. ACA-CPU
corresponds to performing the product $U(V^{*}\mathbf{x})$ in CPU. For
ACA+HOSVD-GPU, $Z_{\rm bc}$ was compressed with Algorithm 3, and the matrix-
vector product was performed with Algorithm 1 in GPU. For $5$ mm3 voxel
resolution, the efficiency is similar for both approaches. For $2$ mm3 voxel
resolution, ACA+HOSVD-GPU outperformed ACA by a factor of $3$, because, due to
its low memory demand, it could be executed on a GPU, whereas ACA could not.

TABLE IX: Time Footprint (hh:mm:ss) of $\mathbf{y}=Z_{\rm bc}\mathbf{x}$ Voxel Res. | Form | PWC | PWL
---|---|---|---
$5$ mm3 | Full-CPU | 00:00:11 | N/A
ACA-CPU | 00:00:01 | 00:00:01
ACA+HOSVD-GPU | 00:00:01 | 00:00:02
$2$ mm3 | Full-CPU | N/A | N/A
ACA-CPU | 00:00:05 | 00:00:18
ACA+HOSVD-GPU | 00:00:02 | 00:00:06
$1$ mm3 | Full-CPU | N/A | N/A
ACA-CPU | N/A | N/A
ACA+HOSVD-GPU | 00:00:09 | 00:00:33

The relative error of $\mathbf{y}$ obtained with ACA+HOSVD-GPU relative to
Full-CPU (ground truth) is shown on the right axis of Fig. 11 for the case of
$5$ mm3 voxel resolution and PWC basis functions. The plot shows how the error
changes as a function of the tolerance ($1e-3$, $\dots$, $1e-8$) used for ACA.
In particular, the relative error remained approximately an order of magnitude
higher than ACA’s tolerance. Fig. 11 also shows plots for the ACA rank and the
maximum Tucker rank of $Z_{\rm bc}$ (values on the left axis). Both ranks
increased as the tolerance of ACA was decreased. We expect similar results for
the other cases, but we were unable to assemble the full coupling matrix, due
to its vast memory footprint.

Figure 11: (left axis) ACA rank and maximum Tucker rank (obtained with HOSVD)
of $Z_{\rm bc}$. (right axis) Error in $\mathbf{y}=Z_{\rm bc}\mathbf{x}$
calculated with ACA+HOSVD-GPU relative to Full-CPU.

##  V Discussion

We showed that our new Tucker-based algorithms could effectively compress the
coupling VSIE matrices in the case of fine voxel resolutions. Thanks to the
achieved compression factors, the matrix-vector products can be performed in
GPU’s, yielding faster execution. The proposed approach will be key for a
VSIE-based in vivo implementation of GMT, for which reducing memory
requirements and computation time are both critical factors.

For cases in which the coil is placed at some distance from the imaging
object, the coupling matrix is low-rank, thus can be compressed using ACA as
$UV^{*}$. For coarse voxel resolutions, such compression strategy alone is
effective and allows for a rapid CPU implementation of matrix-vector products.
However, as the voxel grid of the VIE domain is refined, the memory demand of
$U$ increases until the standard ACA becomes impractical for applications that
need high accuracy, such as GMT. In fact, in order to keep using ACA, one
would have to relax the tolerance, sacrificing accuracy. To avoid this, in
this work we introduced an extension of ACA (Algorithm 3), for which the
memory demand remains as a low as the memory required by one column of the
coupling matrix for any tolerance. Furthermore, Algorithm 3 can be executed in
GPU for rapid computations also in the case of fine voxel resolutions.

Other memory-friendly techniques are available for the fast implementation of
matrix-vector products in VSIE simulations: the magnetic resonance Green’s
function (MRGF), the fast multipole method (FMM), and the precorrected FFT
method (pFFT). The MRGF [11] is a model order reduction technique that can
considerably accelerate the solutions of the VSIE system. However, the
required computational time can be overwhelming when fine voxel resolutions
and PWL basis functions are used. The FMM has been extensively used for the
compression of the MoM matrix [48, 49, 50]. However, in our case, the FMM
could lead to unnecessary overheads, because we are only interested in the
compression of an off-diagonal block of the full MoM matrix (i.e., the
coupling matrix), since the remaining blocks can be handled efficiently with
other methods [13, 12, 10]. Finally, the pFFT method [8] could be used to
project the discretized coil’s elements onto an extended VIE domain, where the
Green’s function tensors are compressible with the Tucker decomposition
(pFFT+Tucker) [10]. However, this approach would be effective only when the
coil is close to the scatterer, like for the close-fitting coil in section
IV.B.1 [51]. In fact, in such situation the extended VIE domain would be
larger than the original one by only a few voxels in each direction, allowing
the matrix-vector products to fit in a GPU. As a result, and also because it
relies on the FFT, such pFFT+Tucker could be more efficient than our proposed
method, although more complex to implement.

An important aspect of the approach presented in this manuscript is that it
can effectively compress the coupling matrix both when the coil is close to or
far from the scatterer. Because of this, our method allows using GPUs to
accelerate the matrix-vector products for most applications. For example, if
the close-fitting coil geometry in Fig. 7 were integrated with a head gradient
insert with a surrounding conductive shield at a certain distance, with our
approach the coupling matrix would still be compressible and fit in the memory
of a GPU. In fact, the interactions between the shield and the head would have
lower Tucker ranks than the ones between the coil’s conductors and the head
(see Section IV.A.1). On the other hand, the pFFT+Tucker method would no
longer be efficient because it would require extending the VIE domain to fully
include the shield. Even though the Green’s function tensors of the extended
domain would still be compressible with Tucker decomposition, the unknowns
that multiply such tensors element-wise would not be. In fact, their
dimensions would have to significantly increase in order to perform the
relevant FFT-based matrix-vector products.

For the previous example, in order to fully exploit the highly parallel
architecture of the GPU with a minimal number of matrix-vector products
operations, one could use a hybrid method that combines pFFT and Algorithm 3.
To do that, the near and far interactions between the VIE domain and the
conducting surfaces would need to be separated. Then, the pFFT could be used
to model the near interactions between the VIE domain and the coil as in [51],
while the arising Green’s function operators could be compressed with the
Tucker decomposition as in [10]. Finally, the remaining far interactions
between the VIE domain and the coil could be modeled with a coupling matrix,
which would be vastly compressed with Algorithm 3. Such hybrid method could
enable us to rapidly execute the matrix-vector products in GPU for most cases,
including complex coil-shield geometries, fine voxel resolutions, and PWL
basis functions. The described hybrid method will be investigated in future
work.

The main limitation of our proposed method is that it requires a considerable
amount of time for the assembly of the coupling matrix, when the matrix is not
compressible with ACA. It is especially slow because each element of the
coupling matrix is a 5D integral. We showed that this could be addressed by
implementing matrix assembly and compression in parallel, but such approach is
not always possible due to memory limitations (IV.B.1). For such cases, one
could alternatively employ 3D cross-Tucker approximation algorithms [40, 18],
which are less efficient than HOSVD for the tensor dimensions in this work,
but do not suffer from memory constraints in case of large tensors. In fact,
3D cross-Tucker methods require only a small number of rows, columns, and
fibers of the tensor they approximate, and they can be implemented with linear
complexity (with respect to tensor’s linear size). In future work, we will
explore the execution of multiple 3D cross-Tucker steps in parallel to avoid
memory overflows when assembling a compressed coupling matrix in the case of
extremely fine resolutions.

##  VI Conclusion

We presented a memory compression technique for the coupling matrix in VSIE
systems. Our method enables one to form and store the coupling matrix even
when its full form size is prohibitively large ($\sim 80$ TB). Specifically,
in this work we were able to achieve a compression between $\sim 0.5$ (PWC)
and $\sim 2$ (PWL) million times when simulating interactions between MRI
coils and realistic body models with some distance between them, in the case
of fine voxel resolutions of $1$ mm3. The error was around one order of
magnitude higher than the tolerance employed for the algorithm. The stored,
compressed matrices could be used multiple times without the need to repeat
the assembly. For example, this would allow one to rapidly perform EM
simulations for the same coil with different body geometries, as far as they
are contained in the original computational domain. For most cases, our
compression method enables fitting large coupling matrices in GPUs, resulting
in rapid execution of the VSIE matrix-vector product (from $1$ to $56$ seconds
for the studied scenarios). Finally, the proposed method could facilitate the
implementation of VSIE-based GMT for in vivo mapping of tissue electrical
properties at clinically-relevant voxel resolutions.

## References

  * [1] S. W. Anderson _et al._ , “Effect of disease progression on liver apparent diffusion coefficient values in a murine model of NASH at 11.7 Tesla MRI,” _Journal of Magnetic Resonance Imaging_ , vol. 33, no. 4, pp. 882–888, 2011.
  * [2] J. Jin and J. Chen, “On the SAR and field inhomogeneity of birdcage coils loaded with the human head,” _Magnetic resonance in medicine_ , vol. 38, no. 6, pp. 953–963, 1997.
  * [3] R. Lattanzi _et al._ , “Electrodynamic constraints on homogeneity and radiofrequency power deposition in multiple coil excitations,” _Magnetic resonance in medicine_ , vol. 61, no. 2, pp. 315–334, 2009.
  * [4] X. Zhang _et al._ , “From complex B1 mapping to local SAR estimation for human brain MR imaging using multi-channel transceiver coil at 7T,” _IEEE transactions on medical imaging_ , vol. 32, no. 6, pp. 1058–1067, 2013\.
  * [5] M. Cosottini _et al._ , “Short-term side-effects of brain MR examination at 7 T: a single-centre experience,” _European radiology_ , vol. 24, no. 8, pp. 1923–1928, 2014.
  * [6] A. Taflove and K. R. Umashankar, “Review of FD-TD numerical modeling of electromagnetic wave scattering and radar cross section,” _Proceedings of the IEEE_ , vol. 77, no. 5, pp. 682–699, 1989.
  * [7] R. Lee and A. C. Cangellaris, “A study of discretization error in the finite element approximation of wave solutions,” _IEEE Transactions on Antennas and Propagation_ , vol. 40, no. 5, pp. 542–549, 1992.
  * [8] J. R. Phillips and J. K. White, “A precorrected-FFT method for electrostatic analysis of complicated 3-D structures,” _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_ , vol. 16, no. 10, pp. 1059–1072, 1997.
  * [9] A. A. Tambova _et al._ , “On the generalization of directfn for singular integrals over quadrilateral patches,” _IEEE Transactions on Antennas and Propagation_ , vol. 66, no. 1, pp. 304–314, 2017.
  * [10] I. I. Giannakopoulos, M. S. Litsarev, and A. G. Polimeridis, “Memory footprint reduction for the fft-based volume integral equation method via tensor decompositions,” _IEEE Transactions on Antennas and Propagation_ , vol. 67, no. 12, pp. 7476–7486, 2019.
  * [11] J. F. Villena _et al._ , “Fast electromagnetic analysis of MRI transmit RF coils based on accelerated integral equation methods,” _IEEE Transactions on Biomedical Engineering_ , vol. 63, no. 11, pp. 2250–2261, 2016.
  * [12] A. G. Polimeridis _et al._ , “Stable FFT-JVIE solvers for fast analysis of highly inhomogeneous dielectric objects,” _Journal of Computational Physics_ , vol. 269, pp. 280–296, 2014.
  * [13] S. Rao, D. Wilton, and A. Glisson, “Electromagnetic scattering by surfaces of arbitrary shape,” _IEEE Transactions on antennas and propagation_ , vol. 30, no. 3, pp. 409–418, 1982.
  * [14] I. P. Georgakis _et al._ , “A fast volume integral equation solver with linear basis functions for the accurate computation of electromagnetic fields in MRI,” _IEEE Transactions on Antennas and Propagation, Early Access_ , 2020.
  * [15] J. E. Serrallés _et al._ , “Noninvasive Estimation of Electrical Properties from Magnetic Resonance Measurements via Global Maxwell Tomography and Match Regularization,” _IEEE Transactions on Biomedical Engineering_ , vol. 67, no. 1, pp. 3–15, 2019.
  * [16] I. Giannakopoulos _et al._ , “Magnetic-resonance-based electrical property mapping using Global Maxwell Tomography with an 8-channel head coil at 7 Tesla: a simulation study,” _IEEE Transactions on Biomedical Engineering_ , vol. 68, no. 1, pp. 236–246, 2021.
  * [17] L. R. Tucker, “Some mathematical notes on three-mode factor analysis,” _Psychometrika_ , vol. 31, no. 3, pp. 279–311, 1966.
  * [18] I. I. Giannakopoulos, M. S. Litsarev, and A. G. Polimeridis, “3D cross-Tucker approximation in FFT-based volume integral equation methods,” in _2018 IEEE International Symposium on Antennas and Propagation & USNC/URSI National Radio Science Meeting_. IEEE, 2018, pp. 2507–2508.
  * [19] J. Zhang, Y. Han, and J. Jiang, “Tucker decomposition-based tensor learning for human action recognition,” _Multimedia Systems_ , vol. 22, no. 3, pp. 343–353, 2016.
  * [20] M. Wang _et al._ , “VoxCap: FFT-Accelerated and Tucker-Enhanced Capacitance Extraction Simulator for Voxelized Structures,” _arXiv preprint arXiv:2004.02609_ , 2020.
  * [21] C. Qian and A. C. Yucel, “On the Compression of Translation Operator Tensors in FMM-FFT-Accelerated SIE Simulators via Tensor Decompositions,” _arXiv preprint arXiv:2010.00520_ , 2020.
  * [22] I. V. Oseledets, “Tensor-train decomposition,” _SIAM Journal on Scientific Computing_ , vol. 33, no. 5, pp. 2295–2317, 2011.
  * [23] B. N. Khoromskij and I. Oseledets, “Quantics-TT collocation approximation of parameter-dependent and stochastic elliptic PDEs,” _Computational Methods in Applied Mathematics Comput. Methods Appl. Math._ , vol. 10, no. 4, pp. 376–394, 2010.
  * [24] L. Grasedyck, “Hierarchical singular value decomposition of tensors,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 31, no. 4, pp. 2029–2054, 2010.
  * [25] E. Tyrtyshnikov, “Mosaic-skeleton approximations,” _Calcolo_ , vol. 33, no. 1-2, pp. 47–57, 1996.
  * [26] E. Tyrtyshnikov, “Mosaic ranks and skeletons,” in _International Workshop on Numerical Analysis and Its Applications_. Springer, 1996, pp. 505–516.
  * [27] S. A. Goreinov, E. E. Tyrtyshnikov, and N. L. Zamarashkin, “A theory of pseudoskeleton approximations,” _Linear algebra and its applications_ , vol. 261, no. 1-3, pp. 1–21, 1997.
  * [28] S. Kurz, O. Rain, and S. Rjasanow, “The adaptive cross-approximation technique for the 3D boundary-element method,” _IEEE transactions on Magnetics_ , vol. 38, no. 2, pp. 421–424, 2002.
  * [29] M. Bebendorf and S. Rjasanow, “Adaptive low-rank approximation of collocation matrices,” _Computing_ , vol. 70, no. 1, pp. 1–24, 2003.
  * [30] R. F. Harrington, _Field computation by moment methods_. Wiley-IEEE Press, 1993.
  * [31] M. F. Catedra, E. Gago, and L. Nuno, “A numerical scheme to obtain the RCS of three-dimensional bodies of resonant size using the conjugate gradient method and the fast Fourier transform,” _IEEE transactions on antennas and propagation_ , vol. 37, no. 5, pp. 528–537, 1989.
  * [32] P. Zwamborn and P. M. Van Den Berg, “The three dimensional weak form of the conjugate gradient FFT method for solving scattering problems,” _IEEE Transactions on Microwave Theory and Techniques_ , vol. 40, no. 9, pp. 1757–1766, 1992.
  * [33] H. Gan and W. C. Chew, “A discrete BCG-FFT algorithm for solving 3D inhomogeneous scatterer problems,” _Journal of Electromagnetic Waves and Applications_ , vol. 9, no. 10, pp. 1339–1357, 1995.
  * [34] J. Jin _et al._ , “Computation of electromagnetic fields for high-frequency magnetic resonance imaging applications,” _Physics in Medicine & Biology_, vol. 41, no. 12, p. 2719, 1996.
  * [35] M. Van Beurden and S. Van Eijndhoven, “Well-posedness of domain integral equations for a dielectric object in homogeneous background,” _Journal of Engineering Mathematics_ , vol. 62, no. 3, pp. 289–302, 2008.
  * [36] J. Markkanen _et al._ , “Analysis of volume integral equation formulations for scattering by high-contrast penetrable objects,” _IEEE Transactions on Antennas and Propagation_ , vol. 60, no. 5, pp. 2367–2374, 2012.
  * [37] P. Yla-Oijala _et al._ , “Surface and volume integral equation methods for time-harmonic solutions of Maxwell’s equations,” _Progress In Electromagnetics Research_ , vol. 149, pp. 15–44, 2014.
  * [38] C.-T. Tai, _Dyadic Green functions in electromagnetic theory_. Institute of Electrical & Electronics Engineers (IEEE), 1994.
  * [39] L. De Lathauwer, B. De Moor, and J. Vandewalle, “A multilinear singular value decomposition,” _SIAM journal on Matrix Analysis and Applications_ , vol. 21, no. 4, pp. 1253–1278, 2000.
  * [40] I. V. Oseledets, D. Savostianov, and E. E. Tyrtyshnikov, “Tucker dimensionality reduction of three-dimensional arrays in linear time,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 30, no. 3, pp. 939–956, 2008.
  * [41] S. A. Goreinov and E. E. Tyrtyshnikov, “The maximal-volume concept in approximation by low-rank matrices,” _Contemporary Mathematics_ , vol. 280, pp. 47–52, 2001.
  * [42] W. Chai and D. Jiao, “Theoretical study on the rank of integral operators for broadband electromagnetic modeling from static to electrodynamic frequencies,” _IEEE Transactions on Components, Packaging and Manufacturing Technology_ , vol. 3, no. 12, pp. 2113–2126, 2013.
  * [43] M. A. Francavilla _et al._ , “Maxwell parallel imaging,” _arXiv preprint arXiv:2008.09042_ , 2020.
  * [44] K. Zhao, M. N. Vouvakis, and J.-F. Lee, “The adaptive cross approximation algorithm for accelerated method of moments computations of EMC problems,” _IEEE transactions on electromagnetic compatibility_ , vol. 47, no. 4, pp. 763–773, 2005.
  * [45] A. Christ _et al._ , “The Virtual Family-development of surface-based anatomical models of two adults and two children for dosimetric simulations,” _Physics in Medicine & Biology_, vol. 55, no. 2, p. N23, 2009\.
  * [46] SIEMENS Healthineers, “Magnetom skyra.” [Online]. Available: https://www.siemens-healthineers.com/magnetic-resonance-imaging/3t-mri-scanner/magnetom-skyra
  * [47] E. Milshteyn _et al._ , “Individualized sar calculations using computer vision-based mr segmentation and a fast electromagnetic solver,” _Magnetic Resonance in Medicine_ , vol. 85, no. 1, pp. 429–443, 2021.
  * [48] L. Greengard and V. Rokhlin, “A fast algorithm for particle simulations,” _Journal of computational physics_ , vol. 73, no. 2, pp. 325–348, 1987.
  * [49] R. Coifman, V. Rokhlin, and S. Wandzura, “The fast multipole method for the wave equation: A pedestrian prescription,” _IEEE Antennas and Propagation Magazine_ , vol. 35, no. 3, pp. 7–12, 1993.
  * [50] B. Shanker and H. Huang, “Accelerated Cartesian expansions–a fast method for computing of potentials of the form R- $\nu$ for all real $\nu$,” _Journal of Computational Physics_ , vol. 226, no. 1, pp. 732–753, 2007.
  * [51] G. G. Guryev _et al._ , “Fast field analysis for complex coils and metal implants in MARIE 2.0.” in _Proc. ISMRM_ , 2019, p. 1035.


# Gröbner-Shirshov bases for some Lie algebras111Supported by the NNSF of
China (11171118), the Research Fund for the Doctoral Program of Higher
Education of China (20114407110007), the NSF of Guangdong Province
(S2011010003374) and the Program on International Cooperation and Innovation,
Department of Education, Guangdong Province (2012gjhz0007).

Yuqun Chen, Yu Li and Qingyan Tang

School of Mathematical Sciences, South China Normal University
Guangzhou 510631, P. R. China
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
Abstract: We give Gröbner-Shirshov bases for Drinfeld-Kohno Lie algebra
$\textbf{L}_{n}$ in [27] and Kukin Lie algebra $A_{P}$ in [32], where $P$ is a
semigroup. As applications, we show that as $\mathbb{Z}$-module
$\textbf{L}_{n}$ is free and a $\mathbb{Z}$-basis of $\textbf{L}_{n}$ is
given. We give another proof of Kukin Theorem: if semigroup $P$ has the
undecidable word problem then the Lie algebra $A_{P}$ has the same property.

Key words:  Gröbner-Shirshov basis, Lie algebra, Drinfeld-Kohno Lie algebra,
word problem, semigroup.

AMS 2000 Subject Classification: 17B01, 16S15, 3P10, 20M05, 03D15

##  1 Introduction

Gröbner bases and Gröbner-Shirshov bases were invented independently by A.I.
Shirshov for ideals of free (commutative, anti-commutative) non-associative
algebras [40, 41], free Lie algebras [39, 41] and implicitly free associative
algebras [39, 41] (see also [2, 4]), by H. Hironaka [29] for ideals of the
power series algebras (both formal and convergent), and by B. Buchberger [20]
for ideals of the polynomial algebras.

Gröbner bases and Gröbner-Shirshov bases theories have been proved to be very
useful in different branches of mathematics, including commutative algebra and
combinatorial algebra, see, for example, the books [1, 18, 21, 22, 24, 25],
the papers [2, 3, 4, 7, 9, 10, 11, 12, 13, 20, 23, 30, 35], and the surveys
[5, 6, 14, 15, 16, 17].

A.A. Markov [34], E. Post [37], A. Turing [42], P.S. Novikov [36] and W.W.
Boone [19] constructed finitely presented semigroups and groups with the
undecidable word problem. This result also follows from the Higman theorem
[28] that any recursive presented group is embeddable into finitely presented
group. A weak analogy of Higman theorem for Lie algebras was proved in [3]
that was enough for existence of a finitely presented Lie algebra with the
undecidable word problem.

In [32], Kukin constructed the Lie algebra $A_{P}$ for a semigroup $P$ such
that if $P$ has the undecidable word problem then $A_{P}$ has the same
property. In this paper we find a Gröbner-Shirshov basis for $A_{P}$ and then
give another proof of the above result.

The Drinfeld-Kohno Lie algebra $\textbf{L}_{n}$ appears in [26, 31] as the
holonomy Lie algebra of the complement of the union of the diagonals
$z_{i}=z_{j},\ i<j$. The universal Knizhnik-Zamolodchikov connection takes
values in this Lie algebra. In this paper, we give a Gröbner-Shirshov basis
for the Drinfeld-Kohno Lie algebra $\textbf{L}_{n}$ over $\mathbb{Z}$ and a
$\mathbb{Z}$-basis of $\textbf{L}_{n}$. As an application we get a simple
proof that $\textbf{L}_{n}$ is an iterated semidirect product of free Lie
algebras.

We are very grateful to Professor L.A. Bokut for his guidance and useful
discussions.

##  2 Composition-Diamond lemma for Lie algebras over a field

For the completeness of the paper, we formulate the Composition-Diamond lemma
for Lie algebras over a field in this section, see [6, 33, 38, 41] for
details.

Let $k$ be a filed, $I$ a well ordered index set, $X=\\{x_{i}|i\in I\\}$ a
set, $X^{*}$ the free monoid generated by $X$ and $Lie(X)$ the free Lie
algebra over $k$ generated by $X$.

We order $X=\\{x_{i}|i\in I\\}$ by $x_{i}>x_{t}$ if $i>t$ for any $i,t\in I$.

We use two linear orders on $X^{*}$: for any $u,v\in X^{*}$,

(i) (lex order) $1\succ t$ if $t\neq 1$ and, by induction, if $u=x_{i}u_{i}$
and $v=x_{j}v_{j}$ then $u\succ v$ if and only if $x_{i}>x_{j}$, or
$x_{i}=x_{j}$ and $u_{i}\succ v_{j}$;

(ii) (deg-lex order) $u>v$ if and only if $deg(u)>deg(v)$, or $deg(u)=deg(v)$
and $u\succ v$, where $deg(u)$ is the length of $u$.

We regard $Lie(X)$ as the Lie subalgebra of the free associative algebra
$k\langle X\rangle$, which is generated by $X$ under the Lie bracket
$[u,v]=uv-vu$. Given $f\in k\langle X\rangle$, denote $\bar{f}$ the leading
word of $f$ with respect to the deg-lex order. $f$ is monic if the coefficient
of $\bar{f}$ is 1.

###### Definition 2.1

An associative word $w\in X^{*}\setminus\\{1\\}$ is an associative Lyndon-
Shirshov word (ALSW for short) if

$(\forall u,v\in X^{*},u,v\neq 1)\ w=uv\Rightarrow w>vu.$

A non-associative word $(u)$ in $X$ is a non-associative Lyndon-Shirshov word
(NLSW for short), denoted by $[u]$, if

(i) $u$ is an ALSW;

(ii) if $[u]=[(u_{1})(u_{2})]$ then both $(u_{1})$ and $(u_{2})$ are NLSW’s;

(iii) if $[u]=[[[u_{11}][u_{12}]][u_{2}]]$ then $u_{12}\preceq u_{2}$.

We denote the set of all ALSW’s and NLSW’s in $X$ by $ALSW(X)$ and $NLSW(X)$
respectively.

For an ALSW $w$, there is a unique bracketing $[w]$ such that $[w]$ is NLSW:
$[w]=[[u][v]]$ if $deg(w)>1$, where $v$ is the longest proper associative
Lyndon-Shirshov end of $w$.

Shirshov Lemma Suppose that $w=aub$, where $w,u\in ALSW(X)$. Then

  1. (i)

$[w]=[a[uc]d],$ where $b=cd$ and possibly $c=1$.

  2. (ii)

Represent $c$ in the form $c=c_{1}c_{2}\ldots c_{n},$ where
$c_{1},\ldots,c_{n}\in ALSW(X)$ and $c_{1}\leq c_{2}\leq\ldots\leq c_{n}$.
Replacing $[uc]$ by $[\ldots[[u][c_{1}]]\ldots[c_{n}]]$ we obtain the word

$[w]_{u}=[a[\ldots[[[u][c_{1}]][c_{2}]]\ldots[c_{n}]]d]$

which is called the Shirshov special bracketing of $w$ relative to $u$.

  3. (iii)

$\overline{[w]}_{u}=w.$

###### Definition 2.2

Let $S\subset Lie(X)$ with each $s\in S$ monic, $a,b\in{X^{*}}$ and $s\in S$.
If $a\bar{s}b$ is an ALSW, then we call
$[asb]_{\bar{s}}=[a\bar{s}b]_{\bar{s}}|_{[\bar{s}]\mapsto{s}}$ a special
normal $S$-word, where $[a\bar{s}b]_{\bar{s}}$ is defined as in Shirshov
Lemma. An $S$-word $(asb)$ is called a normal $S$-word if
$\overline{(asb)}=a\overline{s}b$.

Suppose that $f,\ g\in S$. Then, there are two kinds of compositions:

  1. (i)

If $w=\bar{f}=a\bar{g}b$ for some $a,b\in X^{*}$, then the polynomial
$(f,g)_{w}=f-[agb]_{\bar{g}}$ is called the composition of inclusion of $f$
and $g$ with respect to $w$.

  2. (ii)

If $w$ is a word such that $w=\bar{f}b=a\bar{g}$ for some $a,b\in X^{*}$ with
$deg(\bar{f})+deg(\bar{g})>deg(w)$, then the polynomial
$(f,g)_{w}=[fb]_{\bar{f}}-[ag]_{\bar{g}}$ is called the composition of
intersection of $f$ and $g$ with respect to $w$.

In (i) and (ii), $w$ is called an ambiguity.

Let $h$ be a Lie polynomial and $w\in X^{*}$. We shall say that $h$ is trivial
modulo $(S,w)$, denoted by $h\equiv_{Lie}0\ mod(S,w)$, if
$h=\sum_{i}\alpha_{i}(a_{i}s_{i}b_{i})$, where each $(a_{i}s_{i}b_{i})$ is a
normal $S$-word and $a_{i}\bar{s_{i}}b_{i}<w$.

The set $S$ is called a Gröbner-Shirshov basis in $Lie(X)$ if any composition
in $S$ is trivial modulo $S$ and corresponding $w$.

###### Theorem 2.3

(Composition-Diamond lemma for Lie algebras over a field) Let
$S\subset{Lie(X)}\subset{k\langle X\rangle}$ be nonempty set of monic Lie
polynomials. Let $Id(S)$ be the ideal of $Lie(X)$ generated by $S$. Then the
following statements are equivalent.

  1. (i)

$S$ is a Gröbner-Shirshov basis in $Lie(X)$.

  2. (ii)

$f\in{Id(S)}\Rightarrow{\bar{f}=a\bar{s}b}$ for some $s\in{S}$ and
$a,b\in{X^{*}}$.

  3. (iii)

$Irr(S)=\\{[u]\in NLSW(X)\ |\ u\neq{a\bar{s}b},\ s\in{S},\ a,b\in{X^{*}}\\}$
is a linear basis for $Lie(X|S)=Lie(X)/Id(S)$.

Remark: The above Composition-Diamond lemma is also valid if we replace the
base field $k$ by an arbitrary commutative ring $K$ with identity. If this is
the case, as $K$-module, $Lie(X|S)$ is free with a $K$-basis $Irr(S)$.

##  3 Kukin’s construction of a Lie algebra with unsolvable word problem

Let $P=sgp\langle x,y|u_{i}=v_{i},\ i\in I\rangle$ be a semigroup. Consider
the Lie algebra

$A_{P}=Lie(x,\hat{x},y,\hat{y},z|S),$

where $S$ consists of the following relations:

  1. (1)

$[\hat{x}x]=0,\ [\hat{x}y]=0,\ [\hat{y}x]=0,\ [\hat{y}y]=0$,

  2. (2)

$[\hat{x}z]=-[zx],\ [\hat{y}z]=-[zy]$,

  3. (3)

$\lfloor zu_{i}\rfloor=\lfloor zv_{i}\rfloor,\ i\in I$.

Here, $\lfloor zu\rfloor$ means the left normed bracketing.

In this section, we give a Gröbner-Shirshov basis for Lie algebra $A_{P}$ and
by using this result we give another proof for Kukin’s theorem, see Corollary
3.2.

Let the order $\hat{x}>\hat{y}>z>x>y$ and $>$ the deg-lex order on
$\\{\hat{x},\hat{y},x,y,z\\}^{*}$. Let $\rho$ be the congruence on
$\\{x,y\\}^{*}$ generated by $\\{(u_{i},v_{i}),\ i\in I\\}$. Let

  1. $(3)^{\prime}$

$\lfloor zu\rfloor=\lfloor zv\rfloor,\ (u,v)\in\rho$ with $u>v$.

###### Theorem 3.1

With the above notation, the set $S_{1}=\\{(1),(2),(3)^{\prime}\\}$ is a
Gröbner-Shirshov basis in $Lie(\hat{x},\hat{y},x,y,z)$.

Proof: For any $u\in\\{x,y\\}^{*}$, by induction on $|u|$, $\overline{\lfloor
zu\rfloor}=zu$. All possible compositions in $S_{1}$ are intersection of (2)
and $(3)^{\prime}$, and inclusion of $(3)^{\prime}$ and $(3)^{\prime}$.

For $(2)\wedge(3)^{\prime},\ w=\hat{x}zu,\ (u,v)\in\rho,\ u>v,\
f=[\hat{x}z]+[zx],\ g=\lfloor zu\rfloor-\lfloor zv\rfloor$. We have

$\displaystyle([\hat{x}z]+[zx],\lfloor zu\rfloor-\lfloor
zv\rfloor)_{w}=[fu]_{\bar{f}}-[\hat{x}g]_{\bar{g}}$ $\displaystyle\equiv$
$\displaystyle\lfloor([\hat{x}z]+[zx])u\rfloor-[\hat{x}(\lfloor
zu\rfloor-\lfloor zv\rfloor)]$ $\displaystyle\equiv$ $\displaystyle\lfloor
zxu\rfloor+\lfloor\hat{x}zv\rfloor\equiv\lfloor zxu\rfloor-\lfloor
zxv\rfloor\equiv 0\ \ mod(S_{1},w).$

For $(3)^{\prime}\wedge(3)^{\prime},\ w=zu_{1}=zu_{2}e,\ e\in\\{x,y\\}^{*},\
(u_{i},v_{i})\in\rho,\ u_{i}>v_{i},\ i=1,2$. We have

$\displaystyle(\lfloor zu_{1}\rfloor-\lfloor zv_{1}\rfloor,\lfloor
zu_{2}\rfloor-\lfloor zv_{2}\rfloor)_{w}\equiv(\lfloor zu_{1}\rfloor-\lfloor
zv_{1}\rfloor)-\lfloor(\lfloor zu_{2}\rfloor-\lfloor zv_{2}\rfloor)e\rfloor$
$\displaystyle\equiv$ $\displaystyle\lfloor\lfloor zv_{2}\rfloor
e\rfloor-\lfloor zv_{1}\rfloor\equiv\lfloor zv_{2}e\rfloor-\lfloor
zv_{1}\rfloor\equiv 0\ \ mod(S_{1},w).$

Thus, the set $S_{1}=\\{(1),(2),(3)^{\prime}\\}$ is a Gröbner-Shirshov basis
in $Lie(\hat{x},\hat{y},x,y,z)$. $\blacksquare$

###### Corollary 3.2

(Kukin [32]) Let $u,v\in\\{x,y\\}^{*}$. Then

$u=v\ \mbox{ in the semigroup}\ P\Leftrightarrow\lfloor zu\rfloor=\lfloor
zv\rfloor\ \mbox{ in the Lie algebra }\ A_{P}.$

Proof: Suppose that $u=v\ \mbox{ in the semigroup}\ P$. Without loss of
generality, we may assume that $u=au_{1}b,\ v=av_{1}b$ for some
$a,b\in\\{x,y\\}^{*}$ and $(u_{1},v_{1})\in\rho$. For any $r\in\\{x,y\\}$, by
the relations (1), we have $[\hat{x}r]=0$ and so $\lfloor
zxc\rfloor=\lfloor[z\hat{x}]c\rfloor=[\lfloor zc\rfloor\hat{x}],\ \lfloor
zyc\rfloor=[\lfloor zc\rfloor\hat{y}]$ for any $c\in\\{x,y\\}^{*}$. From this
it follows that in $A_{P}$, $\lfloor zu\rfloor=\lfloor
zau_{1}b\rfloor=\lfloor\lfloor zau_{1}\rfloor b\rfloor=\lfloor\lfloor
zu_{1}\widehat{\overleftarrow{a}}\rfloor b\rfloor=\lfloor
zu_{1}\widehat{\overleftarrow{a}}b\rfloor=\lfloor
zv_{1}\widehat{\overleftarrow{a}}b\rfloor=\lfloor zav_{1}b\rfloor=\lfloor
zv\rfloor$, where $\overleftarrow{x_{i_{1}}x_{i_{2}}\cdots
x_{i_{n}}}=x_{i_{n}}x_{i_{n-1}}\cdots x_{i_{1}}$ and
$\widehat{x_{i_{1}}x_{i_{2}}\cdots
x_{i_{n}}}=\widehat{x_{i_{1}}}\widehat{x_{i_{2}}}\cdots\widehat{x_{i_{n}}}$,
$x_{i_{j}}\in\\{x,y\\}$. Moreover, $(3)^{\prime}$ holds in $A_{P}$.

Suppose that $\lfloor zu\rfloor=\lfloor zv\rfloor\ \mbox{ in the Lie algebra
}\ A_{P}$. Then both $\lfloor zu\rfloor$ and $\lfloor zv\rfloor$ have the same
normal form in $A_{P}$. Since $S_{1}$ is a Gröbner-Shirshov basis in $A_{P}$
by Theorem 3.1, both $\lfloor zu\rfloor$ and $\lfloor zv\rfloor$ can be
reduced to the same normal form of the form $\lfloor zc\rfloor$ for some
$c\in\\{x,y\\}^{*}$ only by the relations $(3)^{\prime}$. This implies that in
$P$, $u=c=v$. $\blacksquare$

By the above corollary, if the semigroup $P$ has the undecidable word problem
then so does the Lie algebra $A_{P}$.

##  4 Gröbner-Shirshov basis for the Drinfeld-Kohno Lie algebra
$\textbf{L}_{n}$

In this section we give a Gröbner-Shirshov basis for the Drinfeld-Kohno Lie
algebra $\textbf{L}_{n}$.

###### Definition 4.1

([27]) Let $n>2$ be an integer. The Drinfeld-Kohno Lie algebra
$\textbf{L}_{n}$ over $\mathbb{Z}$ is defined by generators $t_{ij}=t_{ji}$
for distinct indices $1\leq i,\ j\leq n-1$, and relations

$\displaystyle t_{ij}t_{kl}=0,$ $\displaystyle t_{ij}(t_{ik}+t_{jk})=0,$

where $i,\ j,\ k,\ l$ are distinct.

Clearly, $\textbf{L}_{n}$ has a presentation $Lie_{\mathbb{Z}}(T|S)$, where
$T=\\{t_{ij}|\ 1\leq i<j\leq n-1\\}$ and $S$ consists of the following
relations

$\displaystyle t_{ij}t_{kl}=\ 0,\ \ \ \ k<i<j,\ k<l,\ l\neq\ i,\ j$ (1)
$\displaystyle t_{jk}t_{ij}+t_{ik}t_{ij}=\ 0,\ \ \ \ i<j<k$ (2) $\displaystyle
t_{jk}t_{ik}-t_{ik}t_{ij}=\ 0,\ \ \ i<j<k$ (3)

Now we order $T$: $t_{ij}<t_{kl}$ if either $i<k$ or $i=k$ and $j<l$. Let $<$
be the deg-lex order on $T^{*}$.

###### Theorem 4.2

Let $S=\\{(\ref{a}),(\ref{b}),(\ref{c})\\}$ be as before, $<$ the deg-lex
order on $T^{*}$. Then $S$ is a Gröbner-Shirshov basis for $\textbf{L}_{n}$.

Proof. We list all the possible ambiguities. Denote $(i)\wedge(j)$ the
composition of the type $(i)$ and type $(j)$.

For $(1)\wedge(n)$, $1\leq n\leq 3$, the possible ambiguities $w$’s are:

$\displaystyle(1)\wedge(1)$ $\displaystyle t_{ij}t_{kl}t_{mr},\ (k<i<j,\ k<l,\
l\neq\ i,j,\ m<k<l,\ m<r,\ r\neq\ k,l),$ $\displaystyle(1)\wedge(2)$
$\displaystyle t_{ij}t_{kl}t_{mk},\ (\ k<i<j,\ k<l,\ l\neq\ i,j,\ m<k<l),$
$\displaystyle(1)\wedge(3)$ $\displaystyle t_{ij}t_{kl}t_{ml},\ (k<i<j,\ k<l,\
l\neq\ i,j,\ m<k<l).$

For $(2)\wedge(n)$, $1\leq n\leq 3$, the possible ambiguities $w$’s are:

$\displaystyle(2)\wedge(1)$ $\displaystyle t_{jk}t_{ij}t_{mr},\ (\ m<i<j<k,\
m<r,\ r\neq\ i,j),$ $\displaystyle(2)\wedge(2)$ $\displaystyle
t_{jk}t_{ij}t_{mi},\ (\ m<i<j<k),$ $\displaystyle(2)\wedge(3)$ $\displaystyle
t_{jk}t_{ij}t_{mj},\ (\ m<i<j<k).$

For $(3)\wedge(n)$, $1\leq n\leq 3$, the possible ambiguities $w$’s are:

$\displaystyle(3)\wedge(1)$ $\displaystyle t_{jk}t_{ik}t_{mr},\ (\ m<i<j<k,\
m<r,\ r\neq\ i,k),$ $\displaystyle(3)\wedge(2)$ $\displaystyle
t_{jk}t_{ik}t_{mi},\ (\ m<i<j<k),$ $\displaystyle(3)\wedge(3)$ $\displaystyle
t_{jk}t_{ik}t_{mk},\ (\ m<i<j<k).$

We claim that all compositions are trivial relative to $S$.

Here, we only prove cases $(1)\wedge(1)$, $(1)\wedge(2)$, $(2)\wedge(1)$,
$(2)\wedge(2)$, and the other cases can be proved similarly.

For $(1)\wedge(1)$, let $f=t_{ij}t_{kl},\ g=t_{kl}t_{mr},\ k<i<j,\ k<l,\
l\neq\ i,j,\ m<k<l,\ m<r,\ r\neq\ k,l.$ Then $w=t_{ij}t_{kl}t_{mr}$ and

$\displaystyle(f,g)_{w}$ $\displaystyle=$
$\displaystyle(t_{ij}t_{kl})t_{mr}-t_{ij}(t_{kl}t_{mr})$ $\displaystyle=$
$\displaystyle\ (t_{ij}t_{mr})t_{kl}\ mod(S,w).$

There are three subcases to consider: $r\neq\ i,j$, $r=i$, $r=j$.

Subcase 1. If $r\neq\ i,j$, then

$\displaystyle(t_{ij}t_{mr})t_{kl}\equiv\ 0\ mod(S,w).$

Subcase 2. If $r=i$, then

$\displaystyle(t_{ij}t_{mr})t_{kl}$ $\displaystyle=$
$\displaystyle(t_{ij}t_{mi})t_{kl}$ $\displaystyle\equiv$ $\displaystyle\
t_{kl}(t_{mj}t_{mi})$ $\displaystyle\equiv$ $\displaystyle\
(t_{kl}t_{mj})t_{mi}+t_{mj}(t_{kl}t_{mi})$ $\displaystyle\equiv$
$\displaystyle\ 0\ mod(S,w).$

Subcase 3. If $r=j$, then

$\displaystyle(t_{ij}t_{mr})t_{kl}$ $\displaystyle=$
$\displaystyle(t_{ij}t_{mj})t_{kl}$ $\displaystyle\equiv$ $\displaystyle\
-t_{kl}(t_{mj}t_{mi})$ $\displaystyle\equiv$ $\displaystyle\
-(t_{kl}t_{mj})t_{mi}-t_{mj}(t_{kl}t_{mi})$ $\displaystyle\equiv$
$\displaystyle\ 0\ mod(S,w).$

For $(1)\wedge(2)$, let $f=t_{ij}t_{kl},\ g=t_{kl}t_{mk}+t_{ml}t_{mk},\
k<i<j,\ k<l,\ l\neq\ i,j,\ m<k<l.\ $Then $w=t_{ij}t_{kl}t_{mk}$ and

$\displaystyle(f,g)_{w}$ $\displaystyle=$
$\displaystyle(t_{ij}t_{kl})t_{mk}-t_{ij}(t_{kl}t_{mk}+t_{ml}t_{mk})$
$\displaystyle=$ $\displaystyle\ (t_{ij}t_{mk})t_{kl}-t_{ij}(t_{ml}t_{mk})$
$\displaystyle\equiv$ $\displaystyle\ -t_{ij}(t_{ml}t_{mk})$
$\displaystyle\equiv$ $\displaystyle\
-(t_{ij}t_{ml})t_{mk}-t_{ml}(t_{ij}t_{mk})$ $\displaystyle\equiv$
$\displaystyle\ 0\ mod(S,w).$

For $(2)\wedge(1)$, let $f=t_{jk}t_{ij}+t_{ik}t_{ij},\ g=t_{ij}t_{mr},\
m<i<j<k,\ \ m<r,\ r\neq\ i,j.\ $Then $w=t_{jk}t_{ij}t_{mr}$ and

$\displaystyle(f,g)_{w}$ $\displaystyle=$
$\displaystyle(t_{jk}t_{ij}+t_{ik}t_{ij})t_{mr}-t_{jk}(t_{ij}t_{mr})$
$\displaystyle\equiv$ $\displaystyle\
(t_{jk}t_{mr})t_{ij}+(t_{ik}t_{mr})t_{ij}\ mod(S,w).$

There are two subcases to consider: $r\neq\ k$, $r=k$.

Subcase 1. If $r\neq\ k$, then

$\displaystyle(t_{jk}t_{mr})t_{ij}+(t_{ik}t_{mr})t_{ij}\equiv\ 0\ mod(S,w).$

Subcase 2. If $r=k$, then

$\displaystyle(t_{jk}t_{mr})t_{ij}+(t_{ik}t_{mr})t_{ij}$ $\displaystyle=$
$\displaystyle(t_{jk}t_{mk})t_{ij}+(t_{ik}t_{mk})t_{ij}$ $\displaystyle\equiv$
$\displaystyle\ -t_{ij}(t_{mk}t_{mj})-t_{ij}(t_{mk}t_{mi})$
$\displaystyle\equiv$ $\displaystyle\
(t_{ij}t_{mj})t_{mk}+(t_{ij}t_{mi})t_{mk}$ $\displaystyle\equiv$
$\displaystyle\ -t_{mk}(t_{mj}t_{mi})+t_{mk}(t_{mj}t_{mi})$
$\displaystyle\equiv$ $\displaystyle\ 0\ mod(S,w).$

For $(2)\wedge(2)$, let $f=t_{jk}t_{ij}+t_{ik}t_{ij},\
g=t_{ij}t_{mi}+t_{mj}t_{mi},\ m<i<j<k.\ $Then $w=t_{jk}t_{ij}t_{mi}$ and

$\displaystyle(f,g)_{w}$ $\displaystyle=$
$\displaystyle(t_{jk}t_{ij}+t_{ik}t_{ij})t_{mi}-t_{jk}(t_{ij}t_{mi}+t_{mj}t_{mi})$
$\displaystyle=$ $\displaystyle\
(t_{jk}t_{mi})t_{ij}+(t_{ik}t_{mi})t_{ij}+t_{ik}(t_{ij}t_{mi})-t_{jk}(t_{mj}t_{mi})$
$\displaystyle\equiv$ $\displaystyle\
t_{ij}(t_{mk}t_{mi})-t_{ik}(t_{mj}t_{mi})-t_{jk}(t_{mj}t_{mi})$
$\displaystyle\equiv$ $\displaystyle\
-(t_{ij}t_{mi})t_{mk}+(t_{ik}t_{mi})t_{mj}-(t_{jk}t_{mj})t_{mi}$
$\displaystyle\equiv$ $\displaystyle\
-t_{mk}(t_{mj}t_{mi})-(t_{mk}t_{mi})t_{mj}+(t_{mk}t_{mj})t_{mi}$
$\displaystyle\equiv$ $\displaystyle\ 0\ mod(S,w).$

So $S$ is a Gröbner-Shirshov basis for $\textbf{L}_{n}$. $\blacksquare$

Let $L$ be a Lie algebra over a commutative ring $K$, $L_{1}$ an ideal of $L$
and $L_{2}$ a subalgebra of $L$. We call $L$ a semidirect product of $L_{1}$
and $L_{2}$ if $L=L_{1}\oplus L_{2}$ as $K$-modules.

By Theorems 2.3 and 4.2, we have immediately the following corollaries.

###### Corollary 4.3

The Drinfeld-Kohno Lie algebra $\textbf{L}_{n}$ is a free $\mathbb{Z}$-module
with a $\mathbb{Z}$-basis

$Irr(S)=\\{[t_{ik_{1}}t_{ik_{2}}\cdots\ t_{ik_{m}}]\ |\
t_{ik_{1}}t_{ik_{2}}\cdots\ t_{ik_{m}}\ is\ an\mbox{ ALSW }in\ T^{*},\
m\in\mathbb{N}\\}.$

###### Corollary 4.4

([27]) $\textbf{L}_{n}$ is an iterated semidirect product of free Lie
algebras.

Proof. Let $A_{i}$ be the free Lie algebra generated by $\\{t_{ij}\ |\
i<j\leq\ n-1\\}$. Clearly,

$\textbf{L}_{n}=A_{1}\oplus\ A_{2}\oplus\ \cdots\ \oplus\ A_{n-2}$

as $\mathbb{Z}$-modules, and from the relations $(1),(2),(3)$, we have

$A_{i}\triangleleft\ A_{i}+A_{i+1}+\ \cdots\ +A_{n-2}.$

$\blacksquare$

Remark: In this section, if we replace the base ring $\mathbb{Z}$ by an
arbitrary commutative ring $K$ with identity, then all results hold.

## References

  * [1] W.W. Adams, P. Loustaunau, An introduction to Gröbner bases, Graduate Studies in Mathematics, Vol. 3, American Mathematical Society (AMS), 1994.
  * [2] G.M. Bergman, The diamond lemma for ring theory, Adv. Math. 29 (1978) 178-218.
  * [3] L.A. Bokut, Insolvability of the word problem for Lie algebras and subalgebras of finitely presented Lie algebras, Izvestija AN USSR (mathem.) 36 (1972) 1173-1219.
  * [4] L.A. Bokut, Imbeddings into simple associative algebras, Algebra i Logika 15 (1976) 117-142.
  * [5] L.A. Bokut, Y.Q. Chen, Gröbner-Shirshov bases: Some new results, Proceedings of the Second International Congress in Algebra and Combinatorics, World Scientific, (2008) 35-56.
  * [6] L.A. Bokut, Y.Q. Chen, Gröbner-Shirshov bases and their calculation, arxiv.org/abs/1303.5366.
  * [7] L.A. Bokut, Y.Q. Chen, Y.S. Chen, Composition-Diamond lemma for tensor product of free algebras, _Journal of Algebra_ 323 (2010) 2520-2537.
  * [8] L.A. Bokut, Y.Q. Chen, Y.S. Chen, Groebner-Shirshov bases for Lie algebras over a commutative algebra, _Journal of Algebra_ 337 (2011) 82-102.
  * [9] L.A. Bokut, Y.Q. Chen, X.M. Deng, Gröbner-Shirshov bases for Rota-Baxter algebras, _Siberian Math. J._ 51 (6) (2010) 978-988.
  * [10] L.A. Bokut, Y.Q. Chen, Y. Li, Lyndon-Shirshov basis and anti-commutative algebras, _Journal of Algebra_ 378 (2013) 173-183.
  * [11] L.A. Bokut, Y.Q. Chen, C.H. Liu, Gröbner-Shirshov bases for dialgebras,_International Journal of Algebra and Computation_ 20 (3) (2010) 391-415.
  * [12] L.A. Bokut, Y.Q. Chen, Q.H. Mo, Gröbner-Shirshov bases and embeddings of algebras, _International Journal of Algebra and Computation_ 20 (2010) 875-900.
  * [13] L.A. Bokut, Y.Q. Chen, Q.H. Mo, Gröbner-Shirshov bases for semirings, _Journal of Algebra_ 385 (2013) 47-63.
  * [14] L.A. Bokut, Y.Q. Chen, K.P. Shum, Some new results on Groebner-Shirshov bases, in: Proceedings of International Conference on Algebra 2010, Advances in Algebraic Structures, (2012) 53-102.
  * [15] L.A. Bokut, Y. Fong, W.-F. Ke, P.S. Kolesnikov, Gröbner and Gröbner-Shirshov bases in algebra and conformal algebras, Fundamental and Applied Mathematics 6 (3) (2000) 669-706.
  * [16] L.A. Bokut, P.S. Kolesnikov, Gröbner-Shirshov bases: from their incipiency to the present, J. Math. Sci. 116 (1) (2003) 2894-2916.
  * [17] L.A. Bokut, P.S. Kolesnikov, Gröbner-Shirshov bases, conformal algebras and pseudo-algebras, J. Math. Sci. 131 (5) (2005) 5962-6003.
  * [18] L.A. Bokut, G. Kukin, Algorithmic and Combinatorial algebra, Kluwer Academic Publ., Dordrecht, 1994
  * [19] W.W. Boone, The word problem, Ann. Math. 70 (1959) 207-265.
  * [20] B. Buchberger, An algorithmical criteria for the solvability of algebraic systems of equations, Aequationes Math. 4 (1970) 374-383.
  * [21] B. Buchberger, G.E. Collins, R. Loos, R. Albrecht, Computer algebra, symbolic and algebraic computation, Computing Supplementum, Vol.4, New York: Springer-Verlag, 1982.
  * [22] B. Buchberger, Franz Winkler, Gröbner bases and applications, London Mathematical Society Lecture Note Series, Vol.251, Cambridge: Cambridge University Press, 1998.
  * [23] Y.S. Chen, Y.Q. Chen, Groebner-Shirshov bases for matabelian Lie algebras, _Journal of Algebra_ 358 (2012) 143-161.
  * [24] D.A. Cox, J. Little, D. O’Shea, Ideals, varieties and algorithms: An introduction to computational algebraic geometry and commutative algebra, Undergraduate Texts in Mathematics, New York: Springer-Verlag, 1992.
  * [25] D. Eisenbud, Commutative algebra with a view toward algebraic geometry, Graduate Texts in Math., Vol.150, Berlin and New York: Springer-Verlag, 1995.
  * [26] V.G. Drinfeld, Quasi-Hopf algebras, Algebra i Analiz 1 (1989) 114-148.
  * [27] P. Etingof, A. Henriques, J. Kamnitzer, E.M. Rains, The cohomology ring of the real locus of the moduli space of stable curves of genus 0 with marked points, Ann. Math. 171 (2010) 731-777.
  * [28] G. Higman, Subgroups of finitely presented groups. Proc. Royal Soc. London (Series A) 262 (1961) 455-475.
  * [29] H. Hironaka, Resolution of singularities of an algebraic variety over a field if characteristic zero, I, II, Ann. Math. 79 (1964) 109-203, 205-326.
  * [30] S.-J. Kang, K.-H. Lee, Gröbner-Shirshov bases for irreducible $sl_{n+1}$-modules, Journal of Algebra 232 (2000) 1-20.
  * [31] T. Kohno, Serie de Poincare Koszul associee aux groupes de tresses pures, Invent. Math. 82 (1985) 57-75.
  * [32] G.P. Kukin, On the word problem for Lie algebras, Sibirsk. Mat. Zh. 18 (1977) 1194-1197.
  * [33] Lyndon, R.C.: On Burnside’s problem I, Trans. Amer. Math. Soc. 77 (1954) 202-215.
  * [34] A.A. Markov, Impossibility of some algorithms in the theory of some associative system, Dokl. Akad. Nauk SSSR 55 (1947) 587-590.
  * [35] A.A. Mikhalev, A.A. Zolotykh, Standard Gröbner-Shirshov bases of free algebras over rings, I. Free associative algebras, _International Journal of Algebra and Computation_ 8 (6) (1998) 689-726.
  * [36] P.S. Novikov, On algorithmic undecidability of the word problem in the theory of groups, Trudy MKat. Inst. Steklov. 44 (1955) 1-144.
  * [37] E. Post, A variant of a recursively unsolvable problem, Bull. Amer. Math. Soc. 52 (1946) 264-268.
  * [38] A.I. Shirshov, On free Lie rings, Mat. Sb. 45 (2) (1958) 113-122.
  * [39] A.I. Shirshov, Some algorithmic problem for Lie algebras, Sibirsk. Mat. Zh. 3 (2) (1962) 292-296; English translation in SIGSAM Bull. 33 (1999) 3-6.
  * [40] A.I. Shirshov, Some algorithmic problem for $\varepsilon$-algebras, Sibirsk. Mat. Z. 3 (1962) 132-137.
  * [41] Selected works of A.I. Shirshov, Eds L.A. Bokut, V. Latyshev, I. Shestakov, E. Zelmanov, Trs M. Bremner, M. Kochetov, Birkhäuser, Basel, Boston, Berlin, 2009.
  * [42] A.M. Turing, The word problem in semi-groups with cancellation, Ann. Math. 52 (1950) 191-505.