\section{Introduction}
Modern knowledge-based systems still have difficulties in dealing with information that change over time, known as \emph{dynamic information}. In several real-world scenarios, whatever phenomenon we represent (natural, computational or abstract) is unlikely to be static, therefore representing dynamic information is vital \cite{Fisher2008}. For instance, one can use any knowledge representation model to represent that a given person is the CEO of a company. However, this information no longer holds when another employee assumes that position. This simple example helps to illustrate the need for representing the time span in which a given information is valid.

RDF is a representation language widely used for specifying and querying information in knowledge bases. Although it has proven itself in a widespread range of scenarios, its \emph{subject-predicate-object} (SPO) data model lacks the proper expressiveness for representing temporal information about facts. Different approaches have been proposed over the years for addressing this issue, but their use in current knowledge-based systems is not widespread.

In this work, we approach the problem by using the \emph{hyperknowledge} model \cite{Moreno2017}. Hyperknowledge is a hybrid knowledge representation framework that unifies concepts from hypermedia and knowledge engineering. By using hyperknowledge, one can specify applications linking, for instance, knowledge descriptions and hypermedia anchors (segments of a media object) and  knowledge-aware interactions (e.g., a given content should be presented every time users interact with specific concepts).  Previous work \cite{Moreno2018} extended hyperknowledge with notion of temporal anchors (timed segments in media objects), which allowed the specification of temporal facts in the knowledge base.  In this paper, we revisit that temporal framework and implement it in a query language and engine that is capable of running temporal queries on temporal hyperknowledge models. We illustrate the practical gains by discussing how our approach can be used in a real-world scenario of the Oil \& Gas Industry.

\section{Related Work}

Temporal reasoning in knowledge bases is a well-known problem in AI . There are two key issues that a temporal reasoning framework should address \cite{Vila1994}: (i) an extension to the language/model for representing the temporal aspect of the knowledge, and (ii) a temporal reasoning system. In this paper, we focus on the first issue.

 A proper formalism should provide means to link atemporal assertions (e.g., facts) to temporal references \cite{Pani2001}.  Most of the current AI-based systems rely on the SPO-based RDF data model for representing information. However, representing facts as SPO triples complicates the addition of time information about validity.

 There are some works in literature that address such issue. One of the earliest attempts to add time to RDF was made by Buraga and Ciobanu \cite{Buraga2002}. The authors proposed an extension to the XML representation of RDF triples to add temporal relationships (called \emph{links}) to entities. Their proposal allowed one to express all the Allen operators \cite{Allen1990} among web pages. Despite having a specific scope (expressing timing relationships between websites) and limited expressiveness (it does not promote the specification that a fact is valid over a time period) it helps to illustrate the lack of time awareness in RDF-based languages.

In a well-known work,  \cite{Gutierrez2005} proposes the Temporal RDF model, which extends the SPO data model by labeling triples with temporal values. Such temporal values represent the time in which a triple is valid. This proposal effectively redresses the triple model as a quad model. While being an interesting approach, it does not allow for the explicit specification of multiple validy intervals for a given triple (except for triple reification, which is not ideal in terms of storage and retrieval capacity). A more recent and similar approach to Temporal RDF is presented in \cite{Bereta2013} and but has similar drawbacks.

In \cite{Hoffart2013}, Hoffart \emph{et al.} present YAGO2, an extension of the YAGO knowledge-base that supports the specification of spatiotemporal relationships. Their approach is also based on RDF, using reification to temporaly qualify triples. As mentioned before, such approaches are not ideal in terms of storage and retrieval performance.

There are other works in literature either based on those discussed in this section \cite{Rula2014,Sheth2008} or that extend SPARQL to allow temporal queries in knowledge-bases based on the RDF data model \cite{Lopes2010,Tappolet2009,Perry2011}. RDF lacks enough expressiveness for expressing $n$-ary relationships, which hinders the direct specification of temporal information, leading to several proposals using different techniques to overcome this issue.

Our proposal is built on the hyperknowledge model. Hyperknowledge was proposed in \cite{Moreno2017} as a hybrid knowledge representation model. Hyperknowledge  supports the representationof $n$-ary relationships without reification, which allows one to express using a single relationship all temporal intervals in which a fact holds. This facilitates the implementation of reasoning engines and also avoids the duplication of facts for adding a temporal dimension to them.

\section{Motivational Scenario}
\label{motivational-scenario}

In this section, we present a use case in the Oil \& Gas industry that illustrates the need for representing temporal information in knowledge bases. We will refer to this sceario thoroughout the next sections.

\begin{figure}
	\centering
	\includegraphics[width=0.8\linewidth]{figs/chrono.png}
	\caption{Chronostratigraphic chart of the Scotia Basin, offshore Nova Scotia, Canada \cite{Wade1990} (\emph{Neog} corresponds to Neogene and \emph{E, M, L} abreviate Early, Middle and Late, respectively)}
	\label{fig:chrono}
\end{figure}

A common task in Petroleum Geology is the characterization of basin formation. Geologists must characterize the temporal evolution of  geological formations in a given basin throughout geological time in order to draw conclusions regarding whether it is possible to have a proper condition for hydrocarbon accumulation and/or production. As an illustration, consider Fig.~\ref{fig:chrono}. It depicts a chronostratigraphic chart of the Scotia Basin, located offshore of Nova Scotia, Canada. This type of chart corelates geological periods with litho-stratigraphic units (here called ``facies'') formed on those periods. The first three columns represent geological time in different scales (absolute time in millions of years, periods and epochs, respectively). The fourth column shows information about geological formations/facies and lithology, where each color represents a different formation.

Based on this type of data, geologists can make temporally-situated queries based on temporal relations about the formation of those facies, such as ``\emph{which facies are forming in which geological periods in a given basin}'', or ``\emph{what are the facies of a given type forming after a given epoch}''. These queries require a representation of geological facts that are valid only in a given (geological) timeframe.  They require temporal reasoning in order to calculate temporal relationships between facts.

\section{Preliminaries: The Hyperknowledge Model}

Hyperknowledge was first proposed aiming to better represent relationships among multimedia content (e.g., image, audio, video, text, etc.) as well as conceptual entities, allowing the specification of multimedia-aware knowledge bases. As a first example, consider the model in Fig.~\ref{fig:onto}. It illustrates taxonomic information about geological facies and geolgoical periods present in Fig.~\ref{fig:chrono} as hyperknowledge entities.

\begin{figure}
	\centering
	\includegraphics[width=\linewidth]{figs/onto.png}
	\caption{Ontology using the hyperknowledge model.}
	\label{fig:onto}
\end{figure}

A hyperknowledge data model can be represented as a graph with enhanced vertices. In Fig.~\ref{fig:onto}, there are two types of vertices. Circular labeled vertices are \emph{nodes}, which are entities that represent media content or abstract concepts. For instance, the node labeled \emph{Period} represents the chronological concept of geological periods. Likewise, the node labeled \emph{Geological Facies} represents the homonym abstract concept.

Hyperknowledge nodes are decorated by \emph{anchors}, depicted as small labeled circles on the node boundaries. Anchors are first-class entities that represent (i.e. select) a portion of the node's content. A spatial anchor may represent a subregion of an image. A temporal anchor represents a temporal segment (interval) of a continuous media (e.g., audio, video). In Fig.~\ref{fig:onto}, portions of geological periods that correspond to geological epochs are depicted as labeled anchors. The \emph{lambda} ($\lambda$) anchors represent the whole content of the node.

Darkgrey vertices in Fig.~\ref{fig:onto} are \emph{links}, which define $n$-ary relationships among nodes. One of links' main characteristics is that they can associate nodes only through their anchors.  Links have a \emph{type}, a \emph{predicate}  and some \emph{roles}.  The predicate of a link defines its name. Roles specify the different arguments of the relationship represented by the link. A link type defines its overall role structure. In this paper, we employ \emph{fact} and \emph{hierarchy} links only. Fact links define facts having a subject and an object, similar to SPO triples, but they can also express $n$-ary relationships in hyperknowledge. Hierarchy links are used to define instantiation relationships among nodes. In our scenario, the link named \emph{type} represents a $n$-ary relation between a class and $n-1$ instances of that class.

\section{Representing Temporal Relationships in Hyperknowledge}

In \cite{Moreno2018}, the hyperknowledge model has been extended with a temporal model to represent temporal information. In this section, we revisit that temporal model, providing clearer definitions to some of its main concepts. In the following sections, we employ these notions to propose a temporal query engine.

The main idea of the temporal extension to hyperknowledge is to exploit the association of links and nodes through anchors in order to represent complex temporal relations. A temporal anchor is defined by the hyperknowledge model as being a temporal segment of a node representing a (continuous) media content. For instance, one can define a node that represents a video and create a temporal anchor that represents an interval of that video.

All anchors are assumed to have timing information, including concept nodes. This facilitates the implementation of the reasoning engine (described in next section) because it does not have to handle media nodes differently from others due to their temporal dimension. In fact, any anchor defines a temporal interval and may be used in the temporal algebra implemented by the reasoning system.

We start by defining time. We will consider time to be: (a) \emph{Interval-based} rather then point-based intervals; (b) \emph{unbounded} in the past and future; and (c) \emph{linear}, with no notion of branching, parallelism or circularity. We assume closed time intervals, where $[a, b]$ denotes a closed interval starting at time $a$ and finishing at time $b$. The only exception is with unbounded time intervals $[-\infty,b], [a,+\infty], [-\infty,+\infty]$, in which we assume to be open on the unbounded side.

\begin{definition}
	A \emph{temporal} anchor is an anchor with a ---possibly unbounded --- time interval $[begin, end]$.
\end{definition}

We also refer to a time interval $[begin, end]$ of an anchor $a$ by using the dot notation; i.e.  $a.begin$ and $a.end$.  The interpretation of a temporal anchor is given by its node. For example, a temporal anchor might denote a temporal slice in a piece of media stream or a part of the temporal existence of a given entity.

Temporal anchors can be temporally included in each other:

\begin{definition} Let $x$ and $y$ be two temporal anchors. Anchor $x$ is said to be \emph{included} in $y$ (denoted $in(x, y)$ iff $ y.begin \leq x.begin$ and $x.end \leq y.end$.
\end{definition}

This notion allows us to define temporal lambda anchors. Let $A(n)$ denote the set of anchors of a node $n$. Then::

\begin{definition}A temporal $\lambda$ anchor represents the whole time interval on which a node is defined, such that, given any node $n$, its temporal lambda anchor $\lambda$ and any of its anchors $a \in A(n)$, then $in(n.a,n.\lambda)$.
\end{definition}

An example of temporal lambda anchors are anchors denoting the whole time of a piece of stream, as well as, denoting an entity existence interval.

Reconsider Fig.~\ref{fig:onto}. One can define that the anchors \emph{late}, \emph{middle} and \emph{early} in the \emph{Jurassic} concept node represent time intervals within the Jurassic period. These can be used to represent temporal relationships with facies related to their formation. Fig.~\ref{fig:models1s2} depicts a hyperknowledge specification of the fact that the Scotia Basin has a Carbonate facies that is being formed along the whole of Jurassic Period and early Cretaceus. This fact is modulated by the \emph{when} role, which states that this relationship is true only during the Jurassic Period and the early Cretaceous.

\begin{figure}
\centering
\includegraphics[width=0.9\linewidth]{figs/models1s2.png}
\caption{Hyperknowledge modelling representing S1 and S2.}
\label{fig:models1s2}
\end{figure}

We introduce the role \emph{when} in order to differentiate static and dynamic facts: the former represents aspects of the world that never changes; the latter expresses the dynamic aspect representing an information that may change over the time. Informally, the role \emph{when} defines a set of intervals in which the relationship expressed by a link should be considered valid (or \emph{true}):

\begin{definition} Let $l$ be a link and $W(l)$ be the set of all temporal anchors bound to the \emph{when} role of $l$. The link $l$ is considered \emph{temporally valid} in the interval $[x,y]$ if, for any $t$, such that $x \leq t \leq y$, there is a temporal anchor $a \in W(l)$ such that $a.begin  \leq t \leq a.end$. If $W(l) = \emptyset$, then $l$ is considered atemporally valid (i.e. temporaly valid within $[-\infty,+\infty]$).
\label{def:hklink}
\end{definition}
Note that this definition does not change the semantics of links that do not use the role \emph{when}. For instance, consider the link defining that \emph{Scotia Basin} is an instance of a \emph{Basin.} As this link has no role \emph{when} (i.e., $W(l) = \emptyset$), it is considered valid for any interval. That is, the hierarchy relationship expressed by that link always holds.

\section{Queries}
\label{sec:queries}

In this section, we describe a query language to retrieve information from a knowledge base structured using hyperknowledge.
The queries have the general form ``\emph{\textbf{select} targets \textbf{where} constraints}''. The targets define the type of the output of the query. The possible types are nodes, properties, or anchors.

Query variables refer to \emph{instances} and \emph{classes} of the hyperknowledge model. The constraints are  links, temporal operators, and attribute comparisons. The link constraint has the form ``\emph{predicate(role\_1:a role\_2:b ... role\_n:n)}'', allowing querying $n$-ary relations. Role can be omited if they are not used in constraints. The temporal operators are equivalent to those defined by Allen \cite{Allen1990}, e.g., ``\emph{A before B}''. Finally, constraints may compare the value of a nodes' property with other property or literals.

In the query evaluation process, the reasoning engine must resolve the type hierarchy. In practice, the goal of the engine is to find all instances that have relations and properties that satisfy the query constraints. Therefore, when a node is tested against a constraint, the reasoning engine should check if it is \emph{true} by replacing each instance of that type considering the full type hierarchy.

In our current proposal, there are three types of constraints:

\begin{itemize}
\item
  \emph{Link Constraint:} A link constraint is evaluated to \emph{true}, if a link exists with the given predicate and for each declared role there are nodes that match with the passed identifier, either as an instance or their classes. If the link has a constraint based on a time interval, temporal validty is checked (according with Def. \ref{def:hklink}).
\item
  \emph{Attribute constraint:} It checks if the attribute in an instance (property) satisfies a comparison with the literal of the constraint;
\item
  \emph{Temporal constraints:} Constraint like ``\emph{foo before bar}'', which are evaluated by the reasoning engine by checking if at least one of its anchors satisfies the constraint.
\end{itemize}

Our query language can handle common queries in our domain. For example, fetching all the basins in which carbonates are forming:
\begin{equation}
\begin{aligned}
& \text{\textbf{select }}  \text{{Basin} \textbf{ where }} \\[-4pt]
& \quad \text{\emph{{isFormingIn}}(\emph{subject}:{Carbonate} \emph{object}:Basin)}
\end{aligned}
\end{equation}
This query has only a simple link constraint. It retrieves instances of the node \emph{Basin} that in their links there is a role bound to the instance \emph{Carbonate}. Temporal queries can be specified with the operator \emph{when}, such as retrieving the basin having a carbonate facies forming after the Jurassic period:
\begin{equation}
\label{eqn:2}
\begin{aligned}
& \text{\textbf{select }}  \text{{Basin} \textbf{ where }} \\[-4pt]
	& \quad \text{\emph{{isFormingIn}}(\emph{subject}:{Carbonate} \emph{object}:Basin}\\[-4pt]
	& \quad \quad \text{ \emph{when}:{Period})} \\[-4pt]
	& \quad \text{\textbf{and} {Period} \textbf{after} Jurassic}
\end{aligned}
\end{equation}
In this example, we combine two constraints, the same from previous query and the temporal operation \emph{after}. Here we have used the role \emph{when} and declared the variable \emph{Period} implicitly. Then, we compared \emph{Period} with the interval defined by the \emph{lambda} anchor of the \emph{Jurassic} instance. Anchors can also be refered in the queries, such as when selecting the periods after the Jurassic in wich a carbonate facies has been identified being in formation.
\begin{equation}
\begin{aligned}
&\text{\textbf{select }} \text{\textbf{anchor from} {Period} \textbf{ where }} \\[-4pt]
& \quad \text{\emph{{isFormingIn}}(\emph{subject}:{Carbonate} \emph{object}:Basin}\\[-4pt]
& \quad \quad \text{ \emph{when}:{Period})} \\[-4pt]
& \quad \text{\textbf{and} {Period} \textbf{after} Jurassic}
\end{aligned}
\end{equation}

\section{Implementation}

We have implemented a knowledge-base system (called Hyperknowledge base, or HKBase{} as shorthand) that uses hyperknowledge as data model. Fig.~\ref{fig:architecture} depicts the main components of the HKBase{} architecture.

\begin{figure}
\centering
\includegraphics[width=0.9\linewidth]{figs/architecture.pdf}
\caption{The Hyperknowledge Base architecture overview.}
\label{fig:architecture}
\end{figure}

The HKBase{} implements a CRUD API for knowledge curation and has functionalities for handling multimedia data. The architecture includes two storages, one to store the hyperknowledge content, and a second one to store unstructured content, such as multimedia data that need to be uploaded to the system.  The \emph{Information Extraction Orchestrator} (IEO) is a component that selects an apropriate information extraction service to get concepts from every new media content referenced by the hyperknowledge representation. For instance, if a hyperknowledge node represents an image having an URL property pointing to a valid and accessible file, the HKBase{} is capable of retrieving that image, and to choose an information extraction service (if available) specialized in processing that type of content to extract semantic information from it. The \emph{data structurer} receives the output from IEO and structures it before saving it to the hyperknowledge base. Finally, a reasoning engine processes queries made by the user.

\subsection{Mapping Hyperknowledge to a Knowledge Graph}

The HKBase{} does not rely on a specific storage service, but rather, it defines a driver API (called \emph{IDB -- Interface for DataBase}) which should be implemented by any given storage to be used. The main functionality of an IDB driver is to map the hyperknowledge representation maintained by the HKBase{} to the data model supported by the underlying database. We currently have the Janus Graph\footnote{http://janusgraph.org/} implementation that fully supports the temporal reasoning.

The Janus Graph driver uses the Apache TinkerPop framework\footnote{http://tinkerpop.apache.org/} to convert the hyperknowledge entities to a knowledge graph representation supported. To implement this mapping, it uses different types of vertices and edges. A hyperknowledge node is represented by a vertex named \emph{node} and an anchor of hyperknowledge nodes is represented by an \emph{anchor} vertex. Anchor vertices hold in an internal table the properties of the corresponding hyperknowledge anchor.

Hyperknowledge links are represented by a third type of vertex called \emph{link}, which holds the predicate and roles in its properties table. This design allows representation of $n$-ary relationships, as those discussed above. This addresses the issue of edges having the limitation of connecting exactly two vertices in a graph.

Likewise, different types of edges are used to represent the relationships among vertices in the graph. A \emph{hasAnchor} edge connects a node to its anchors. A \emph{bind} edge connects an \emph{anchor} to a \emph{link} vertice. Bind edges specialize in its properties table to which link role that anchor is bound.

Fig.~\ref{fig:example} depicts the knowledge graph that represents the hyperknowledge model in Fig.~\ref{fig:models1s2}. Arrows represent edges and circles represent vertices in the figure.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{figs/example.png}
\caption{Example of representation of the hyperknowledge model in a Graph Database.}
\label{fig:example}
\end{figure}

The full conversion of all hyperknowledge entities to a knowledge graph requires additional types of vertices (\emph{context, switch, connector, etc.}) and edges (\emph{hasConnector, hasChildren, etc).} However, it is out of the scope of this paper the description of the whole process, but we have focused on the entities involved in the temporal reasoning.

\subsection{Temporal Reasoning}

The reasoning engine receives a query as input, processes it, and returns a set corresponding to the results of the query. For the case of temporal queries, the reasoning process is implemented in a two-stages pipeline composed of the following steps: i) data retrieval; and ii) computation of temporal relations and constraints.

Consider Query \ref{eqn:2} in Section \ref{sec:queries}. In the first stage, the reasoning engine retrieves all links that have the predicate \emph{isFormingIn}, the role \emph{subject} bound to an instance of the class \emph{Carbonate}, the role \emph{object} bound to the node \emph{Basin}, and the role \emph{when} bound to an instance of \emph{Period}. The IDB Driver implements this retrieval by using the Gremlin query language\footnote{https://tinkerpop.apache.org/gremlin.html}.

The output of the first stage becomes the input to the second, which is entirely performed by the reasoning engine. For Query \ref{eqn:2}, it computes which of the returned periods satisfy the temporal constraint \emph{``after Jurassic''}.

The main advantage of this two-stages pipeline for processing queries is that the reasoning process becomes independent of storage services. That is, if all temporal calculus would be performed by the underlying database, the architecture of the HKBase{} would be dependent of that service. Using this two-stages processing, the IDB driver is responsible only for retrieving data from the database, which is an operation supported by any storage service.

\section{Final Remarks}

In this paper, we have revisited an approach for representing temporal information in hyperknowledge representation language. Temporal anchors provide a useful mechanism for expressing richer temporal relationships in hypermedia knowledge bases, being suitable for expressing the intervals in which a given fact holds. We also introduce a query language and engine that can take advantage of these constructs for temporal queries. As future work, we plan to finish the implementation of other IDB drivers and report a quantitative comparison analyses of our reasoning engine working with different storage services. We believe that such quantitative analysis can bring valuable insights to the Multimedia community.  We also intend to conduct user experiments regarding the temporal reasoning and report in a future work the qualitative analysis of users' perspective about the work described in this paper.

\bibliographystyle{IEEEbib}


\section{Introduction}
Orthogonal time frequency space (OTFS) is a kind of modulation technologies which can strongly combat the doppler effect in time-variant channels by spreading each data symbol on the delay-Doppler (DD) plane over the entire time-frequency (TF) plane~\cite{DBLP:journals/corr/abs-1802-02623}. Thus, OTFS is very robust in high-mobility scenarios, which is identified as a main challenge of future wireless communications.

Most existing work assumes that ideal waveforms are adopted in OTFS systems, which ensures the transmitter- and receiver-side waveforms satisfy bi-orthogonality~\cite{7925924,DBLP:journals/corr/abs-1808-00519,8786203}. This assumption can reap the advantages that the inter-symbol interference (ISI) and inter-carrier interference (ICI) are both perfectly eliminated, and constant channel gains are guaranteed for each data symbol in the DD domain~\cite{8424569}. However, the ideal waveform assumption is unavailable in practical scenarios, and the ISI and ICI is inevitably introduced, which causes non-uniform channel gains in the DD domain. As a consequence, channel matrix is no longer doubly circular~\cite{8424569}, which prohibitively increases the complexity of channel equalization, and hence hinders the realization of practical OTFS systems.

Therefore, low-complexity equalization plays a pivotal role in practical OTFS systems. Considering the potential of complexity reduction, linear equalization is usually much preferred, especially when the fractional doppler shift occurs, destroys the sparsity of the equivalent channel in the DD domain, and hence disables the feasibility of nonlinear equalization, e.g., message passing (MP)~\cite{7925924,DBLP:journals/corr/abs-1903-09441}. For the widely-used linear equalization, i.e., zero-forcing (ZF) and minimum mean square error (MMSE), the complexities of the existing direct-matrix-inversion-based schemes are both ${\cal O}( {{{\left( {NM} \right)}^3}} )$, where $M$, $N$ are the numbers of OTFS symbols and subcarriers, respectively. Typically, the value of $NM$ is in the order of thousand or even larger. Thus, even for the linear equalization of practical OTFS systems with non-ideal rectangular waveforms, the complexity is still overly intensive, and heavily limits the realization of practical OTFS systems.

To this end, we concentrate on low-complexity linear equalization algorithms for practical OTFS systems in this paper. Specifically, we first prove the block-circulant property and special Fourier transform structure of the effective channel matrix under rectangular waveforms. Then, low-complexity ZF and MMSE equalization algorithms are proposed by exploiting the proved property and structure. Our proposed algorithms significantly decrease the computational complexities from ${\cal O}( {{{\left( {NM} \right)}^3}} )$ to $\mathcal{O}\left( {{M^2}N\log N} \right)$ for ZF and ${\cal O}\left( {{M^2}{N^2}P} \right)$ for MMSE, where $P$ is the number of channel multi-paths.

In this paper, scalars, matrixes, and vectors are denoted by boldface capital, boldface lowercase, and normal letters, respectively. ${{\bf{I}}_N}$ is the identity matrix of size $N$. ${{\bf{A}}^{ - 1}}$, ${{\bf{A}}^{\rm{T}}}$, ${{\bf{A}}^{\rm{H}}}$ represent the inverse, transpose, and conjugate transpose of ${\bf{A}}$. The operation $\otimes $ represents Kronecker product. The operation ${\rm{circ}}\left\{ {{{\bf{a}}}} \right\}$ forms cyclic square matrix ${{\bf{A}}}$ whose first column is the vector ${{{\bf{a}}}}$, and ${\rm{circ}}\left\{ {{{\bf{A}}_1}, \cdot  \cdot {{\bf{A}}_N}} \right\}$ forms a block-circulant square matrix whose first column block as ${({\bf{A}}_{_1}^{\rm{T}}, \cdot  \cdot {\bf{A}}_{_N}^{\rm{T}})^{\rm{T}}}$. ${{\bf{F}}_M},{\bf{F}}_M^{\rm{H}} \in {\mathbb{C}^{M \times M}}$ are $M$-point DFT and IDFT matrices, respectively. ${[ {{{\bf{A}}}} ]_{i,p}}$ is the ${i_{{\rm{th}}}}$ row, ${p_{{\rm{th}}}}$ column element of ${{{\bf{A}}}}$, ${[ {{{\bf{a}}}} ]_{p}}$ is the ${p_{{\rm{th}}}}$ element of ${{{\bf{a}}}}$, and ${[ {{{\bf{A}}}} ]_{i}}$ is the ${i_{{\rm{th}}}}$ row vector of ${{{\bf{A}}}}$. $a:b$ denotes the indexes increasing from $a$ to $b$ one by one, and ${\left\langle  \cdot  \right\rangle _N}$ denotes the operation of modulus $N$.

\section{OFDM-Based OTFS Systems}
The realization of OTFS is compatible with OFDM by adding a pre-processing and post-processing block~\cite{7925924}, i.e., OFDM technologies can be adopted for multicarrier time-frequency signal transmission in OTFS systems~\cite{8119540,8422467,DBLP:journals/corr/abs-1903-09441}. In this section, a brief review of OFDM-based OTFS systems is given, and then the effective channel matrices of OTFS systems are derived for rectangular and ideal waveforms, respectively.

\subsection{Elements of OFDM-based OTFS Systems}
The time-frequency plane is defined as an $M \times N$ lattice which is sampled by interval $T$ (second) along the time axis and $\Delta f$ (Hz) along the frequency axis, respectively, i.e.,
\[\Lambda {\rm{ = }}\left\{ {\left( {m\Delta f,nT} \right):m = 1,...,M,n = 1,...,N} \right\}.\]
Delay-Doppler plane is defined as an \begin{small} $M \times N$ lattice \[\Gamma {\rm{ = }}\left\{ {\left( {\frac{l}{{M\Delta f}},\frac{k}{{NT}}} \right):l = 1,...,M,k = 1,...,N} \right\}.\vspace{-0.2cm}\]\end{small

As a pre-processing block added to OFDM systems at the transmitter, \emph{inverse symplectic finite Fourier transform} (ISFFT) converts the modulated data in DD domain to TF domain. ISFFT can be written in a matrix form as
\begin{equation}
{{\bf{X}}_{{\rm{TF}}}} = {\rm{ISFFT(}}{{\bf{X}}_{{\rm{DD}}}}{\rm{) = }}{{\bf{F}}_M}{{\bf{X}}_{{\rm{DD}}}}{\bf{F}}_N^{\rm{H}},
\end{equation}
where ${{\bf{X}}_{{\rm{TF}}}}$, \!${{\bf{X}}_{{\rm{DD}}}}\!\in \!{\mathbb{C}^{M \times N}}$ denote the discrete transmit data in a matrix form mapped on $\Lambda $ and $\Gamma $, respectively. At the receiver, the post-processing block additionally executes \emph{symplectic finite Fourier transform} (SFFT), which converts the output signal ${{\bf{R}}_{{\rm{TF}}}}$ on $\Lambda $ plane to DD-domain signal ${{\bf{Y}}_{{\rm{DD}}}}$,
\begin{equation}
{{\bf{Y}}_{{\rm{DD}}}} = {\rm{SFFT(}}{{\bf{R}}_{{\rm{TF}}}}{\rm{) = }}{\bf{F}}_M^{\rm{H}}{{\bf{R}}_{{\rm{TF}}}}{{\bf{F}}_N}.
\end{equation}
For the brevity of subsequent analysis, a DD-domain generalized input-output model is built for vectorized signal as
\begin{equation}\label{endtoend}
{{\bf{y}}_{{\rm{DD}}}} = {{\bf{H}}_{{\rm{e}}{\rm{f}}{\rm{f}}}}{{\bf{x}}_{{\rm{DD}}}}{\rm{ + }}{\bf{w}},
\end{equation}
where ${{\bf{x}}_{{\rm{DD}}}},{{\bf{y}}_{{\rm{DD}}}}\!\! \in\!\! {\mathbb{C}^{NM \times 1}}$ are obtained by column-wise stacking ${{\bf{X}}_{{\rm{DD}}}},{{\bf{Y}}_{{\rm{DD}}}}$ into column vectors, respectively. ${{\bf{H}}_{{\rm{eff}}}} \!\!\in\!\! {\mathbb{C}^{NM \times NM}}$ denotes the DD-domain end-to-end effective channel matrix and ${\bf{w}}$ is the noise vector. Note that we assume the same type of waveform at the transmitter and receiver. Details of ${{\bf{H}}_{{\rm{eff}}}}$ are discussed as follows.

\subsection{Effective Channel Matrix with Rectangular Waveforms}
In~\cite{8422467}, rectangular waveforms are considered in an OFDM-based OTFS system. The effective channel ${{\bf{H}}_{{\rm{eff}}}}$ can be expressed as
\begin{equation}\label{heff}
{{\bf{H}}_{{\rm{eff}}}}\!\! =\!\!(\underbrace {{{\bf{F}}_N} \!\otimes\! {\bf{F}}_M^{\rm{H}}}_{{\text{SFFT}}})\!{{{\bf{\Pi }}_{\rm{r}}}}\!(\!\underbrace{{{{\bf{G}}_{{\rm{r}}}}} \!\otimes \!{{\bf{F}}_M}}_{{{\text{OFDM Demod}}\text{.}}}\!){\bf{\tilde H}}(\underbrace{{{{{\bf{G}}_{{\rm{t}}}}} \!\otimes\! {\bf{F}}_M^{\rm{H}}}}_{{\text{OFDM Mod.}}}){{{\bf{\Pi }}_{\rm{t}}}}(\underbrace{{{\bf{F}}_N^{\rm{H}} \!\otimes \!{{\bf{F}}_M}}}_{\text{ISFFT}\!}),
\end{equation}
where ${{{\bf{\Pi }}_{\rm{r}}}},{{{\bf{\Pi }}_{\rm{t}}}} \!\!\in\!\! {\mathbb{C}^{MN \times MN}}$ and ${{\bf{G}}_{{\rm{r}}}},{{\bf{G}}_{{\rm{t}}}}\!\!\in\!\!{\mathbb{C}^{N \times N}}$ denote the receive, transmit windows and receive, transmit waveform matrices, respectively. Note that ${{{\bf{\Pi }}_{\rm{r}}}},{{{\bf{\Pi }}_{\rm{t}}}}\!\!= \!\!{{\bf{I}}_{MN}}$, and ${{\bf{G}}_{\rm{r}}},{{\bf{G}}_{\rm{t}}}\!\! =\!\! {{\bf{I}}_N}$ when the rectangular windows and waveforms are considered~\cite{8422467}. ${{\bf{\tilde H}}}\!\in \!{\mathbb{C}^{MN \times MN}}$ is the time-domain channel impulse response matrix of the vectorized time-domain signal after cyclic prefix (CP) adding and removing operations~\cite{8422467}. ${{\bf{\tilde H}}}$ is a block-diagonal matrix when the length of CP is larger than the maximum delay of the channel, and ${\bf{\tilde H}}$ is expressed as ${\bf{\tilde H}}\!\! = \!\!{\rm{diag}}({{\bf{\tilde H}}_1},{{\bf{\tilde H}}_2} \cdots {{\bf{\tilde H}}_N})$, where ${{\bf{\tilde H}}_p} \!\!\in\!\! {\mathbb{C}^{M \times M}}$. According to the tap-delay-line model~\cite{8119540}, there are only $P$ non-zero elements in each column and each row of ${{{\bf{\tilde H}}}_p}$ at fixed positions determined by the channel multi-path delay positions on the time axis, which is denoted as ${\bf{d}} \!\!=\! \!\left[ {{D_1},{D_2} \cdots {D_P}} \right]$. Therein, ${{D_P}}$ is the maximum delay position. Finally, for our considered system, the effective channel matrix in (\ref{heff}) is simplified to
\begin{equation}\label{recheff}
 {\bf{H}}_{_{{\rm{eff}}}}^{{\rm{rect}}} = \left( {{{\bf{F}}_N} \otimes {{\bf{I}}_M}} \right){\bf{\tilde H}}\left( {{\bf{F}}_N^{\rm{H}} \otimes {{\bf{I}}_M}} \right).
\end{equation}

\subsection{Channel Mismatch Under Ideal Waveform Assumption}
In practical wireless systems, if the effect of practical waveforms is not considered, i.e., assuming the ideal waveforms which cause that each data symbol experiences the constant channel gain at the receiver~\cite{8786203}, the systems may suffer the channel mismatch issues. We term this case as the ideal waveform assumption. In such a case, the channel gain of the first data symbol of ${{\bf{x}}_{{\rm{DD}}}}$, i.e., the first column of ${\bf{H}}_{{\rm{eff}}}^{{\rm{rect}}}$, is wrongly deemed as the constant channel gain, i.e., $[{\bf{H}}_{{\rm{eff}}}^{{\rm{rect}}}]_1^{\rm{T}}$. Furthermore, $[{\bf{H}}_{{\rm{eff}}}^{{\rm{rect}}}]_1^{\rm{T}}$ can be expressed as ${[ {{\bf{h}}_1^{\rm{T}},{\bf{h}}_2^{\rm{T}},...,{\bf{h}}_N^{\rm{T}}} ]^{\rm{T}}}$, wherein ${\bf{h}}_p \!\!\in \!\! {\mathbb{C}^{M \times 1}} $. Then, we can reshape $[{\bf{H}}_{{\rm{eff}}}^{{\rm{rect}}}]_1^{\rm{T}}$ into a matrix form as ${{\bf{H}}_{{\rm{DD}}}} \!\!=\!\![{{\bf{h}}_{\rm{1}}},{{\bf{h}}_{\rm{1}}},...,{{\bf{h}}_N}]\!\!\in\!\! {\mathbb{C}^{M \times N}}$. Hereby, each received data symbol can be rewritten by a two-dimension circular convolution between ${{\bf{H}}_{{\rm{DD}}}}$ and the transmit data as~\cite{DBLP:journals/corr/abs-1808-00519},
\begin{equation}\label{twocirc}
{\left[ {{{\bf{Y}}_{{\rm{DD}}}}} \right]_{k,l}}\!\! = \!\!\sum\limits_{k' = 1}^M  \sum\limits_{l' = 1}^N {{{\left[ {{{\bf{H}}_{{\rm{DD}}}}} \right]}_{k',l'}}{{\left[ {{{\bf{X}}_{{\rm{DD}}}}} \right]}_{{{\left\langle {k - k'} \right\rangle }_M} + 1,{{\left\langle {l - l'} \right\rangle }_N} + 1}}} .
\end{equation}
Based on (\ref{twocirc}), the mismatched effective channel matrix ${\bf{H}}_{{\rm{eff}}}^{{\rm{ide}}}$ under the ideal waveform assumption is given by
\begin{equation}\label{ideeff}
    {\bf{H}}_{{\rm{eff}}}^{{\rm{ide}}} \!=\! {\rm{circ}}\left\{ {{\rm{circ}}\left\{ {\bf{h}}_1 \right\},{\rm{circ}}\left\{ {\bf{h}}_2 \right\},\! \cdots \! ,{\rm{circ}}\left\{ {\bf{h}}_N \right\}} \right\}.
\end{equation}

\section{Low-Complexity ZF and MMSE Equalizations}
In this section, we first explore a deep insight into ${\bf{H}}_{{\rm{eff}}}^{{\rm{rect}}}$. Then, low-complexity ZF and MMSE equalization methods are proposed for practical OTFS systems with rectangular waveforms based on the structure feature of ${\bf{H}}_{{\rm{eff}}}^{{\rm{rect}}}$.

In order to facilitate the subsequent analysis, we first introduce two new operators. ${\rm{FF}}{{\rm{T}}_{{\rm{Mtx}}}}( \cdot )$ is defined as the Fourier transform between two sets of matrices. Taking ${{{\bf{\tilde H}}}_p},{{\bf{A}}_n} \!\in\! \mathbb{C}^{M \times M}$ as example, ${{{\bf{\tilde H}}}_p}{\rm{ = FF}}{{\rm{T}}_{{\rm{Mtx}}}}({{\bf{A}}_n})$ when the relationship ${{{\bf{\tilde H}}}_p} \!\!=\! \!\sum\nolimits_{n = 1}^N {{e^{-j\frac{{2\pi }}{N}\left( {p - 1} \right)\left( {n - 1} \right)}}}{{\bf{A}}_n}$ is satisfied. Similarly, we can define ${\rm{IFF}}{{\rm{T}}_{{\rm{Mtx}}}}( \cdot )$ by ${{\bf{A}}_n}{\rm{ = IFF}}{{\rm{T}}_{{\rm{Mtx}}}}({{{\bf{\tilde H}}}_p})$. Actually, when observing the fix-position element, e.g., $(m,l)$-element ${{[ {{{\bf{\tilde H}}}_p}]}_{m,l}}$ in the matrix ${{{\bf{\tilde H}}}_p}$, and re-organizing the fix-position elements, e.g., $(m,l)$-elements, from ${{{\bf{A}}_n}}$ as a vector by ${{\bf{a}}^{m,l}}\!\! = \!\!( {{{[{{\bf{A}}_1}]}_{m,l}},...,{{[{{\bf{A}}_N}]}_{m,l}}}) \!\!\in\!\! {\mathbb{C}^{N \times 1}}$, one can find ${{[ {{{\bf{\tilde H}}}_p}]}_{m,l}} =\!\! {\rm{FFT(}}{{{\bf{a}}^{m,l}}})$, where ${\rm{FFT(}}\cdot)$ is the traditional Fourier transform which can be implemented by fast Fourier transform (FFT). Since there are $M \times M$ different ${{\bf{a}}^{m,l}}$, ${{{\bf{\tilde H}}}_p}{\rm{ = FF}}{{\rm{T}}_{{\rm{Mtx}}}}({{\bf{A}}_n})$ is actually $M^2$ times Fourier transform operations. Thus, ${\rm{FF}}{{\rm{T}}_{{\rm{Mtx}}}}( \cdot )$, and ${\rm{IFF}}{{\rm{T}}_{{\rm{Mtx}}}}( \cdot )$ are both linear operators, which will be used for the following theorem.
\begin{theorem}\label{th1}
\emph{${\bf{H}}_{{\rm{eff}}}^{{\rm{rect}}}$ is a block-circulant matrix, i.e., ${{\bf{H}}_{{\rm{eff}}}^{{\rm{rect}}}}\!\! =\!\!{\rm{circ}}\{ {{\bf{A}}_1},{{\bf{A}}_2}\! \cdots,\!{{\bf{A}}_N}\} $, and the submatrices ${{\bf{A}}_n}$ have a Fourier transform relationship with channel impulse response ${{{\bf{\tilde H}}}_p}$ as ${{\bf{A}}_n}\! \!=\! \!{\rm{IFF}}{{\rm{T}}_{{\rm{Mtx}}}}({{{\bf{\tilde H}}}_p})$.}
\end{theorem}
\begin{proof}
${\bf{H}}_{{\rm{eff}}}^{{\rm{rect}}} \!\!\in\!\! {\mathbb{C}^{NM \times NM}}$ can be viewed as a block matrix, including $N \!\times\! N $ small square matrices denoted by ${\bf{H}}_{i,k}^{{\rm{rect}}}\! \in \!{\mathbb{C}^{M \!\times M}}$
, and from (\ref{recheff}), ${\bf{H}}_{i,k}^{{\rm{rect}}}$ is calculated as
\begin{equation}\label{heffgeneral}
{\bf{H}}_{i,k}^{{\rm{rect}}}= \sum\limits_{p = 1}^N {{f_{ip}}{f'_{pk}}} {{\bf{\tilde H}}_{p}}, i,k,p = 1,...,N,
\end{equation}
where ${f_{ip}} \!\!= \!\!{[ {{{\bf{F}}_N}} ]_{i,p}},{f'_{pk}}\!\! = \!\!{[ {{\bf{F}}_N^{\rm{H}}} ]_{p,k}}$. ${f_{ip}}{f'_{pk}}$ is further derived as
\begin{small}
\begin{equation}\label{fourpro}
\!\!{f_{ip}}{f'_{pk}}\!\! = \!\!{e^{ - j\frac{{2\pi }}{N}\left( {i - 1} \right)\left( {p - 1} \right)}}{e^{j\frac{{2\pi }}{N}\left( {p - 1} \right)\left( {k - 1} \right)}}\!\! = \!\!{e^{j\frac{{2\pi }}{N}\left( {p - 1} \right)\left( {k - i} \right)}}.
\end{equation}
\end{small}
The right hand side of (\ref{fourpro}) is only related to the difference of $i$ and $k$. By defining $n = {\left< {k\! - \! i} \right>_N} \!+ \!1$, and substituting $n$ into (\ref{fourpro}), we can obtain ${f_{ip}}{f'_{pk}}\!\! =\!\! {e^{j\frac{{2\pi }}{N}\left( {p - 1} \right)\left( {n - 1} \right)}}\!\!\triangleq\!\! { f_{n,p}}$, and then (\ref{heffgeneral}) can be simplified to
\begin{equation}\label{simgener}
{\bf{H}}_n^{{\rm{rect}}} = \sum\limits_{p = 1}^N {{f_{n,p}}} {{{\bf{\tilde H}}}_p} = {\rm{IFF}}{{\rm{T}}_{{\rm{Mtx}}}}({{{\bf{\tilde H}}}_p}).
\end{equation}
By defining ${{\bf{A}}_n} = {\bf{H}}_n^{{\rm{rect}}} \!\in\! {\mathbb{C}^{M \times M}}$, ${\bf{H}}_{{\rm{eff}}}^{{\rm{rect}}}$ can be rewritten in a block-circulant matrix form as
\begin{equation}\label{hcric}
{\bf{H}}_{{\rm{eff}}}^{{\rm{rect}}} \!\!= \!\!{\rm{circ}}\{ {{\bf{A}}_1},{{\bf{A}}_2} \cdots ,{{\bf{A}}_N}\},{{\bf{A}}_n} \!\!= \!\!{\rm{IFF}}{{\rm{T}}_{{\rm{Mtx}}}}({{{\bf{\tilde H}}}_p}).
\end{equation}
Then, the proof of Theorem~\ref{th1} is finished.
\end{proof}

As discussed before, under the ideal waveform assumption, the submatrices of the effective channel belong to circular matrices. In this case, the low-complexity ZF and MMSE schemes are more easily attained, since ${\bf{H}}_{{\rm{eff}}}^{{\rm{ide}}}$ can be diagonalized by FFT matrix thanks to the doubly circular property~\cite{8786203}.
Oppositely, with rectangular waveforms, the data symbol in DD domain does not experience a constant channel any more, which causes that ${{\bf{A}}_n}$ in~\eqref{hcric} is not a circulant matrix. The appealing diagonalization operation is not available for rectangular waveforms, which calls for new low-complexity methods based on Theorem~\ref{th1}.

\subsection{Low-Complexity ZF Equalization}
According to the vectorized model in (\ref{endtoend}), the vectorized output signal ${\bf{\hat x}}_{_{{\rm{DD}}}}^{\rm{ZF}}$ of ZF equalizer can be written by
\begin{equation}\label{zf}
{\bf{\hat x}}_{_{{\rm{DD}}}}^{{\rm{ZF}}} = {{\bf{W}}_{{\rm{ZF}}}}{{\bf{y}}_{{\rm{DD}}}} = {\bf{H}}_{{\rm{eff}}}^{ - 1}{{\bf{y}}_{{\rm{DD}}}},
\end{equation}
where ${{\bf{W}}_{{\rm{ZF}}}}$ denotes the evaluation matrix for ZF. Since direct inversion requires overwhelmingly high computational complexity, we develop a low-complexity algorithm. By referring to~\cite{1143132}, the inversion of a block-circulant matrix, i.e., ${\bf{H}}_{{\rm{eff}}}^{ - 1}$, is also a block-circulant matrix. Thus, we can represent ${{\bf{W}}_{{\rm{ZF}}}}$ as
\begin{equation}\label{hinv}
{{\bf{W}}_{{\rm{ZF}}}} = {\bf{H}}_{{\rm{eff}}}^{ - 1} = {\rm{circ}}\{ {{\bf{B}}_1},{{\bf{B}}_2} \cdots {{\bf{B}}_N}\},
\end{equation}
where ${{\bf{B}}_q} \!\!\in \!\!{\mathbb{C}^{M \times M}}$ can be obtained by ${{\bf{A}}_n}$ as~\cite{1143132}
\begin{align}\label{bp}
\nonumber {{\bf{B}}_q}\! = \!\frac{1}{N}\sum\limits_{t = 1}^N {{{\left( {{e^{j\frac{{2\pi (p - 1)}}{N}}}} \right)}^{N - t + 1}}{\bf{S}}_t^{ - 1}}& \! = \! \frac{1}{N}{\rm{FF}}{{\rm{T}}_{{\rm{Mtx}}}}({\bf{S}}_t^{ - 1}),\\
{{\bf{S}}_t} = \sum\limits_{n = 1}^N {{e^{j\frac{{\left( {t - 1} \right)\left( {k - 1} \right)}}{N}}}{{\bf{A}}_n}} & = {\rm{IFF}}{{\rm{T}}_{{\rm{Mtx}}}}({{\bf{A}}_n}),
\end{align}
with ${{\bf{S}}_t} \in {\mathbb{C}^{M \times M}}$. According to \eqref{hinv}, (\ref{bp}), to obtain ${{\bf{W}}_{{\rm{ZF}}}}$, the direct inversion of ${\bf{H}}_{{\rm{eff}}}^{{\rm{rect}}} \in {\mathbb{C}^{NM \times NM}}$ can be avoided by calculating the inversion of ${{\bf{S}}_t}$ and executing FFT operations. Moreover, it is observed that ${{\bf{S}}_t}$ has the same structure of ${{\bf{A}}_n}$. By exploiting the sparse property and structure knowledge of ${{\bf{S}}_t}$, we further develop a LU factorization method to obtain ${\bf{S}}_t^{ - 1}$ to proceed with the simplification of ZF equalization algorithm. The entire details of the proposed low-complexity ZF equalization are provided in Algorithm~\ref{alg1}.

\begin{algorithm}[tb]
\caption{Low-Complexity ZF Equalization}
\label{alg1}
\begin{algorithmic}[1]
\begin{small}
\INPUT  Delay position vector ${\bf{d}}$, multi-path number $P$, ${{{\bf{\tilde H}}}_p}$.\\
        \STATE  Define ${\bf{\Phi}} \!\!= \!\!{{\bf{I}}_M}$, vector ${\bf{u}}\!\! =\!\! [{\bf{d}},M \!- \!{D_P}]$, ${\bf{Y}}\!\! =\!\! {{\bf{0}}_M} \in {\mathbb{C}^{M \times M}}$.\\
        \STATE Compute submatrices ${{\bf{A}}_n}{\rm{ = IFF}}{{\rm{T}}_{{\rm{Mtx}}}}({{{\bf{\tilde H}}}_p})$.\\
        \STATE Compute ${{\bf{S}}_t}{\rm{ = IFF}}{{\rm{T}}_{{\rm{Mtx}}}}({{{\bf{A}}}_n})$.
\FOR {$t=1$ to $N$}
        \STATE  From $i \!= \!1$ to $M \!- \!{D_P}$, ${\left[ {\bf{\Phi }} \right]_{i + {\bf{u}},i}} = {\textstyle{{{{\left[ {{{\bf{S}}_t}} \right]}_{i + {\bf{u}},i}}} \over {{{\left[ {{{\bf{S}}_t}} \right]}_{i,i}}}}}
\FOR{$n = 1$ to $P+1$}
      \STATE $i = {[{\bf{u}}]_n} + 1:{[{\bf{u}}]_{n + 1}}$;
       \STATE  ${\left[ {{\bf{LU}}} \right]_i} = {\left[ {{{\bf{S}}_t}} \right]_i} - \sum\limits_{k = 2}^n {\sum\limits_{j = i}^M {{{\left[ {{\bf{\Phi}}} \right]}_{i,i - {D_k}}}{{\left[ {{\bf{LU}}} \right]}_{i - {D_k},j}}} } $.
     \ENDFOR
       \STATE Update the remaining part of ${\bf{\Phi}}$ by regular LU factorization.
\FOR{$n = 1$ to $P$}
      \STATE $k = {[{\bf{u}}]_n} + 1:{[{\bf{u}}]_{n + 1}}$;
      \STATE ${\left[ {\bf{Y}} \right]_k} = {\left[ {{{\bf{I}}_M}} \right]_k} - \sum\limits_{i = 1}^n {{{\left[ {{\bf{\Phi}}} \right]}_{k,k - {\bf{u}}(i)}}{{\left[ {\bf{Y}} \right]}_{k - {\bf{u}}(i)}}} $.
      \ENDFOR
      \STATE From $k \!\!=\!\! {D_P} \!+\! 1$ to $M$, ${\left[ {\bf{Y}} \right]_k} = {\left[ {{{\bf{I}}_M}} \right]_k} - {\left[ {{\bf{\Phi}}} \right]_{k,1:k - 1}}{\left[ {\bf{Y}} \right]_{1:k - 1}}$.
      \STATE From $k \!\!=\!\! N$ to $1$, set $f$ as the maximum between $M - {D_{P}}$ and ${k + 1}$, then ${\left[ {{\bf{S}}_t^{ - 1}} \right]_k} = {\left[ {\bf{Y}} \right]_k} - {\left[ {\bf{\Phi }} \right]_{k,f:M}}{\left[ {\bf{Y}} \right]_{f:M}}$.
    \ENDFOR
    \STATE Compute ${{\bf{B}}_q}{\rm{ = }}\frac{1}{N}{\rm{FF}}{{\rm{T}}_{{\rm{Mtx}}}}({\bf{S}}_{_t}^{ - 1})$.

    \OUTPUT Evaluation matrix ${{\bf{W}}_{{\rm{ZF}}}}\! = \!{\rm{circ}}\{ {{\bf{B}}_1},{{\bf{B}}_2} \cdots {{\bf{B}}_N}\}$.
\end{small}
\end{algorithmic}
\end{algorithm}

\subsection{Low-Complexity MMSE Equalization}
ZF may cause the noise enhancement issue, which can be addressed by the MMSE equalization given by
\begin{equation}\label{mmse}
{\bf{\hat x}}_{_{{\rm{DD}}}}^{{\rm{MMSE}}}\!\! = \!\!{{\bf{W}}_{{\rm{MMSE}}}}{{\bf{y}}_{{\rm{DD}}}} \!\! = \!\! {\left( {{\bf{H}}_{{\rm{eff}}}^{\rm{H}}{{\bf{H}}_{{\rm{eff}}}} + {\sigma ^2}{{\bf{I}}_{MN}}} \right)^{ - 1}}{\bf{H}}_{{\rm{eff}}}^{\rm{H}}{{\bf{y}}_{{\rm{DD}}}},
\end{equation}
where ${{\bf{W}}_{{\rm{MMSE}}}}$ denotes the evaluation matrix for MMSE, and ${{\sigma ^2}}$ is the noise energy.
\begin{proposition}\label{th2}
The evaluation matrix ${{\bf{W}}_{{\rm{MMSE}}}}$ is a block-circulant matrix and the submatrices can be expressed as
\begin{equation}\label{mmseblock}
\begin{small}
{{\bf{W}}_{{\rm{MMSE}}}} = {\rm{circ}}\left\{ {{{{\bf{\tilde W}}}_1},{{{\bf{\tilde W}}}_2} ,\cdots ,{{{\bf{\tilde W}}}_N}} \right\},
\end{small}
\end{equation}
where ${{{\bf{\tilde W}}}_k} \!\!=\!\! \sum\nolimits_{i = 1}^N {{\bf{C}}_i^{ - 1}{\bf{A}}_{{{\left[ {i -k} \right]}_N}\! +\! 1}^{\rm{H}}}\! \!\! \in\! \!\!{\mathbb{C}^{M \times M}}$, and ${\bf{C}}_i^{ - 1}\!\! \in\!\! {\mathbb{C}^{M \times M}}$ is the submatrices of the block-circulant matrix ${{\bf{C}}^{ - 1}} = {\left( {{\bf{H}}_{{\rm{eff}}}^{\rm{H}}{{\bf{H}}_{{\rm{eff}}}} + {\sigma ^2}{\bf{I}}} \right)^{ - 1}} = {\rm{cric}}\{ {\bf{C}}_1^{ - 1},{\bf{C}}_2^{ - 1}, \cdots ,{\bf{C}}_N^{ - 1}\}$.
\end{proposition}
\begin{proof}
The conjugate transpose of effective channel can be expressed as ${\bf{H}}_{{\rm{eff}}}^{\rm{H}}\!\! = \!\!{\rm{circ}}\{ {\bf{A}}_1^{\rm{H}},{\bf{A}}_N^{\rm{H}},{\bf{A}}_{N - 1}^{\rm{H}}, \cdots ,{\bf{A}}_2^{\rm{H}}\}$, and hence ${\bf{C}}\!\! =\!\! \left( {{\bf{H}}_{{\rm{eff}}}^{\rm{H}}{{\bf{H}}_{{\rm{eff}}}} + {\sigma ^2}{\bf{I}}} \right)$ is the product of two block-circulant matrices plus a diagonal matrix. The product of two block-circulant matrices is also a block-circulant matrix, which means both ${\bf{C}}$ and ${{\bf{C}}^{ - 1}}$ are block-circulant matrices. Thus, ${{\bf{W}}_{{\rm{MMSE}}}}$ is a block-circulant matrix, and a general expression of submatrices in ${{\bf{W}}_{{\rm{MMSE}}}}$ can be obtained by using block-circulant matrix multiplication twice, as shown in~(\ref{mmseblock}).
\end{proof}

Rewrite ${\bf{C}}$ as ${\bf{C}} \!\!=\!\! {\rm{cric}}\{ {\bf{C}}_1^{},{\bf{C}}_2^{}, \cdots ,{\bf{C}}_N^{}\}$. According to (\ref{mmseblock}), ${{\bf{C}}^{{\rm{ - 1}}}} \!\!\in\!\! {\mathbb{C}^{NM \times NM}}$ is required when obtaining ${{\bf{W}}_{{\rm{MMSE}}}}$. We can apply (\ref{bp}) to simplify the inversion operation as
\begin{equation}\label{mmseinv}
{\bf{C}}_q^{ - 1} = \frac{1}{N}{\rm{FF}}{{\rm{T}}_{{\rm{Mtx}}}}({\bf{S}}_t^{ - 1}),{{\bf{S}}_t} = {\rm{IFF}}{{\rm{T}}_{{\rm{Mtx}}}}({{\bf{C}}_n}).
\end{equation}
Based on the analysis above, the direct inversion of MMSE equalization can also be avoided by using a similar algorithm as Algorithm~\ref{alg1}.

\subsection{Complexity Analysis}
As there are $M^2$ times $N$-point FFT (IFFT) operations, the computational complexity of ${\rm{FF}}{{\rm{T}}_{{\rm{Mtx}}}}( \cdot )$ (${\rm{IFF}}{{\rm{T}}_{{\rm{Mtx}}}}( \cdot )$) is $\mathcal{O}\left( {{M^2}N\log N} \right)$. Moreover, the LU factorization method reduces the complexity of direct inversion of ${\bf{S}}_t^{}$ from ${\cal O}\left( {{M^3}} \right)$ to ${\cal O}\left( {{M^2}{D_P}} \right)$, usually with ${D_P} \ll M$. As for the multiplication of the matrices in (\ref{mmseblock}), by exploiting the sparsity of the matrices, the complexity is reduced from ${\cal O}( {{{\left( {NM} \right)}^3}} )$ to ${\cal O}\left( {{M^2}{N^2}P} \right)$. Finally, Table~\uppercase\expandafter{\romannumeral1} shows the computational complexity comparison between the direct-matrix-inversion-based equalization and the proposed low-complexity equalization. The effective channel matrix in the DD domain is obtained by (\ref{simgener}), which is also simplified via the operation ${\rm{IFF}}{{\rm{T}}_{{\rm{Mtx}}}}( \cdot )$. For $N=M=32,P=6$, the overall complexity can be decreased by a few thousand times for ZF and a few hundred times for MMSE. More complexity reduction can be achieved for practical systems with larger $N,M,P$.

\section{Simulation Results}
In this section, simulation results are provided. The parameters are set as follows: The carrier frequency is 5~GHz, the subcarrier space $\Delta f$ is 15~KHz, the number of subcarrier $M$ is $64$, the number of OTFS symbols $N$ in a transmission block is $32$, and the channel is modeled as Vehicular Channel B~\cite{ahmadi2007channel}. Hereby, ${D_P}$ is 20 and let ${L_{{\rm{CP}}}}\!\! =\!\! {D_P}\! +\! 1$. The 4-QAM modulated signal is adopted and the maximum doppler spread is ${f_{\max }} \!\!=\!\! 1$~KHz. The simulation results are averaged over 20000 times realizations.
\begin{figure}[!htbp]
  \centering
    \vspace{-0.3cm}

  \includegraphics[width=3in]{fig1.eps}\\
  \setlength{\abovecaptionskip}{0pt}
\setlength{\belowcaptionskip}{-2cm}
  \caption{BER performance comparison of different schemes.}\label{fig1}
    \vspace{-0.3cm}
\end{figure}

Fig.~\ref{fig1} indicates that our proposed low-complexity equalization algorithms result in no performance loss compared with the direct-inversion-based schemes~\cite{8422467}. The significant performance degradation is observed by the inappropriate ideal-waveform assumption~\cite{8786203}, due to that the interference under rectangular waveforms gives rise to the variation of channel gains among different data symbols in the DD domain, which causes a mismatch between the theoretical and practical channels. Note that a two-stage equalization scheme, resembling traditional feedback decision equalization, has been recently proposed in~\cite{DBLP:journals/corr/abs-1709-02505}. From Fig.~\ref{fig1} and Table~\uppercase\expandafter{\romannumeral1} ($\Gamma$ is the complexity induced by the hard-decision in the second stage), the complexity of our proposed schemes is at least a few hundreds times lower than that of~\cite{DBLP:journals/corr/abs-1709-02505} with slight performance loss during the signal-noise-ratio (SNR) range of 5-16~dB. Besides, more advantages of our schemes are as follows: The linear equalization only requires the easily-obtained DD domain channel~\cite{8422467}, while the two-stage equalization requires the intractable frequency-domain response of time-varying channel in practical OTFS systems; the two-stage equalization imposes higher hardware requirements on practical systems.

\begin{table}[!htbp]
  \centering
  \vspace{-0.3cm}
  \caption{Complexities of Different Schemes}\label{time}
  \vspace{-0.1cm}
  \begin{tabular}{|c|c|}
  \hline
  \bfseries  Scheme  &\bfseries Complexity \\
  \hline
  Direct-inversion-based ZF/MMSE~\cite{8422467}  &  $\mathcal{O}( {{{( {NM} )}^3}} )$  \\
  \hline
  Two-stage equalization in~\cite{DBLP:journals/corr/abs-1709-02505}  & $\mathcal{O}( {{{( {NM} )}^3}} ){\rm{ + }}\Gamma $\\
  \hline
  Proposed low-complexity ZF/MMSE    & ${\cal O}\!( {{M^2}N\!\log \!N } \!)$/${\cal O}( {{M^2}{N^2}P} )$ \\
  \hline
  \end{tabular}
  \vspace{-0.5cm}
\end{table}

\section{Conclusion}
In this paper, we identify and prove the effective channel matrix structure of OTFS systems with rectangular waveforms. Based on the proved structure feature, new low-complexity linear equalization methods are proposed for ZF and MMSE, respectively. Analysis and simulation results validate that the proposed linear equalization methods can significantly reduce the computational complexity with no performance loss, which enables the future applications of OTFS technologies to practical communication systems.

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\bibliographystyle{IEEEtran}


\section*{Acknowledgments}
I.P.I. thanks the Institute of Modern Physics, Lanzhou, China, for financial support and hospitality during his stay.
I.P.I. acknowledges funding from the Portuguese
\textit{Fun\-da\-\c{c}\~{a}o para a Ci\^{e}ncia e a Tecnologia} (FCT) through the FCT Investigator
contract IF/00989/2014/CP1214/CT0004 and project PTDC/FIS-PAR/29436/2017,
which are partially funded through POCI, COMPETE, Quadro de Refer\^{e}ncia
Estrat\'{e}gica Nacional (QREN), and the European Union.
I.P.I. also acknowledges the support from National Science Center, Poland,
via the project Harmonia (UMO-2015/18/M/ST2/00518).
P.M.Z. is supported by the National Natural Science Foundation of China (Grant No. 11975320).
A.V.P. and N.K. thank the Chinese Academy of Sciences President's International Fellowship Initiative for the support
via grants No. 2019PM0036 (A.V.P.) and No. 2017PM0043 (N.K.).


\section{Introduction}

Granular materials are conglomerates of discrete, macroscopic, solid particles, such as sand, soils, pills, etc. They are ubiquitous in many industrial and natural processes.
Analytical study is usually difficult due to the complicated nonlinear heterogeneous multi-body interactions in the dense granular system. Numerical methods have remarkable significance in the study of granular materials since most perturbative analytical methods in condense matter physics requires either the low density in the ideal gas limit or a prefect lattice for crystals. On the meanwhile, it is difficult for experimental methods to deal with large number of particles,  measure physical quantities in 3D and investigate phase transitions such as jamming phenomena quantitatively \cite{Kou:2017, Zhang:2014}. Numerical simulations allow access to all detailed information of the particle system, and usually it is straightforward to calculate all the relavant quantities (e.g. fabric, stress, energy) of the system.

However, it remains challenging to develop efficient and robust numerical methods for granular system, especially for large particle systems. Molecular dynamics and energy minimization are two main classes of simulation methods. Discrete element method (soft-particle method) is one of the typical molecular dynamics methods in granular applications \cite{Luding:2008}, where a set of Langevin equations of motion are formulated, with relaxation and random fluctuation terms to quantify the physical relaxation time scale and thermal effect, respectively. Related methods include event-driven method, non-smooth contact dynamics, etc, see references in \cite{Krijgsman:2016}. On the other hand, energy minimization methods can be applied to the situation of quasi-static deformation, when the applied load changes slowly over time with respect to the inertial forces. One thus can significantly speed up the simulation if time scales are not needed to be fully resolved. Also, the high accuracy of the energy minimization methods provides advantage for the quantitative study of the jamming transition and the scaling law behavior close to the jamming point. Although molecular dynamics method have been extensively used in the simulation of granular materials, the application of energy minimization methods is relatively recent. For example, Luding et. al. have applied trust region methods to study the relaxation and shear of granular system \cite{Krijgsman:2016}. In this paper, we will study the performance of several typical minimization methods for the granular simulation, such as conjugate gradient methods and quasi-Newton type methods.

Preconditioning is the main bottleneck for the development of efficient and robust numerical algorithms for large scale molecular simulation \cite{Alrachid:2018, Packwood:2016}, boh in the case of molecular dynamics and energy minimization. General purpose preconditioners are not specifically targeted at large-scale atomistic systems, and are not particularly effective. In this paper, we will also stress on the development of preconditioning techniques for the efficient energy minimization of granular systems.

\begin{Remark}
The molecular details of molecular dynamics methods such as soft-particle methods and energy minimization methods are usually not the same. However, it is important to note that the equilibrium physical quantities such as fabric, stress and energy remain similar \cite{Krijgsman:2016}. Also, there is no temperature in the granular model simulated by energy minimization methods, namely, we are only interested in  athermal granular materials .
\end{Remark}

\subsection{Outline}
The paper is organized as follows. In Section \ref{sec:setup}, we present the physical model for the granular system, as well as the critical jamming transition which is important for the numerical investigation of dense granular system. We present the optimization and preconditioning techniques used for granular simulation in Section \ref{sec:numerics}. The numerical methods are validated and benchmarked in the numerical experiments in Section \ref{sec:numerical experiment}. We conclude in Section \ref{sec:conclusion} and point out some future improvements.

\section{Physical Setup}
\label{sec:setup}

\def\mathrm{tot}{\mathrm{tot}}

In this section, we introduce the basic physical models. The potential energy of the system is the sum of interaction energy of all particles. It is important to note that particles can interact unless they overlap. In other words, there is only repulsive force therefore the system is qualitatively different from the repulsive-attractive forces such as Lennard-Jones \cite{Edward:1924}. For simplicity, we only consider the contact models without tangential forces, and also only the two dimensional case.

\subsection{Granular configuration}
An admissible configuration $\mathcal{C}$ of particles is given by a list of center $\vec{x}_i$ and radius $r_i$ of particles. The volume (or packing) fraction $\phi$ is a key parameter in the simulation of  granular materials,
$$
\phi=\frac{1}{V}\sum_{i}\pi r_i^2
$$
where $V$ is the volume of the containing domain (container).

\begin{figure}[!htbp]
\begin{center}
\begin{tikzpicture}
\draw [black,thick,fill=black!20] (1.5,3.5) circle (1cm) node[left] {$i$};
\fill [black] (1.5,3.5) circle (1.5pt);
\draw [blue,-] (1.5,3.5) -- (2.1,4.3) node[above] at (1.79,3.91) {$r_{i}$};
\draw [black,thick,fill=black!20] (3.5,2) circle (1.5cm) node[below] {$j$};
\fill [black] (3.5,2) circle (1.5pt);
\fill [black] (2.3,2.9) circle (1.5pt) node[right] {$c$};
\draw [blue,-] (3.5,2) -- (5,2) node[above] at (4.25,2) {$r_j$};
\draw [thick, ->] (-0.75,0) -- (5,0) node[above] {$x$};
\draw [thick, ->] (0,-0.75) -- (0,5) node[right] {$y$};
\draw [black,->] (0,0) -- (1.5,3.5) node[above] at (0.7,1.8) {$\vec{x}_{i}$};
\draw [black,->] (0,0) -- (3.5,2) node[below] at (1.75,1) {$\vec{x}_{j}$};
\draw [black,->] (1.5,3.5) -- (3.5,2);
\path[draw] (2.5,2.75) -- (4.25,4.5) -- (4.75,4.5);
\node [above] at (4.5,4.5) {$\vec{l}_c$};
\node [below] at (-0.25,0) {$O$} ;
\end{tikzpicture}
\end{center}
\caption{Two particles in contact}
\label{fig:contact}
\end{figure}

The geometric structure of the granular materials can be characterized by the following fabric tensor,
\begin{equation}
F=\frac{1}{V}\sum_{c \in C}\pi\left(r_i^2 + r_j^2\right)\vec{n}_c\otimes \vec{n}_c
\end{equation}

\subsection{Contact Model}

Here we will introduce a contact potential of the following form,
\begin{equation} \label{eq: pair potential}
E_{ij} = \left\{
    \begin{aligned}
           & k_{ij}\frac{(r_i+r_j - \left|\vec{x}_i-\vec{x}_j\right|)^{\alpha}}{\alpha} &  &\text{for} \quad \left|\vec{x}_i-\vec{x}_j\right|< r_i+r_j, \\
           & 0 &  &\text{for} \quad \left|\vec{x}_i-\vec{x}_j\right|\ge r_i+r_j, ,
    \end{aligned}
\right.   i\ne j
\end{equation}
where $k_{ij}$ denotes the particle stiffness, $r_i$ and $r_j$ are the radii of particle $i$ and $j$, respectively, $\vec{x}_i$ and $\vec{x}_j$ are the positions of the their centers, respectively.

The total potential energy of the system is,
\begin{equation}\label{eq:Etot}
     E_{\mathrm{tot}} = \sum_{i<j}E_{ij}.
\end{equation}
See Figure \ref{fig:contact} for an illustration.

Since the particles interact with each other only if they contact, we can rewrite the potential energy for a single contact $c$ as,
\begin{equation}
E_c = k_c\frac{(\delta_c\vee 0)^{\alpha}}{\alpha},
\end{equation}
where $\delta_c$ is the overlap associated with the contact $c$, and by equation \eqref{eq: pair potential},
\begin{equation}
\delta_c = r_i + r_j - \left|\vec{x}_i-\vec{x}_j\right|.
\end{equation}
In this formula, particle $i$ and $j$ are two particles associated with contact point $c$ (see Figure \ref{fig:contact}).  Therefore, if $\delta_c\geq 0$, the first derivatives, or forces of $E_c$ are,
\begin{equation} \label{eq: first derivative}
  \begin{aligned}
    & \frac{\partial E_c}{\partial \vec{x}_i} = -\vec{f}_c = -k\delta_c^{\alpha-1}\vec{n}_c \\
    & \frac{\partial E_c}{\partial \vec{x}_j} = \vec{f}_c = k\delta_c^{\alpha-1}\vec{n}_c,
  \end{aligned}
\end{equation}
where
\begin{equation}
\vec{n}_c = \frac{\vec{l}_c}{l_{c}} ,\quad \vec{l}_c = \vec{x}_i-\vec{x}_j, \quad l_c = \left| \vec{l}_c \right|,
\end{equation}
and $\vec{f}_c$ is the contact force acting on particle $i$ by particle $j$.
Furthermore, the second derivatives are,
\begin{equation} \label{eq: second derivative}
  \begin{aligned}
    & \frac{\partial^2 E_c}{\partial\vec{x}_i \partial\vec{x}_i} = k\delta_c^{\alpha-2} \left(\left(\alpha-1+\frac{\delta_c}{l_c}\right)\vec{n}_c\otimes \vec{n}_c - I \frac{\delta_c}{l_c} \right) \\
    & \frac{\partial^2 E_c}{\partial\vec{x}_j \partial\vec{x}_j} = \frac{\partial^2 E_c}{\partial\vec{x}_i \partial\vec{x}_i}\\
    & \frac{\partial^2 E_c}{\partial\vec{x}_i \partial\vec{x}_j} = -\frac{\partial^2 E_c}{\partial\vec{x}_i \partial\vec{x}_i},
  \end{aligned}
\end{equation}
where $\otimes$ is the tensor product of two vectors and $I$ is the $2\times 2$ identity matrix.

\begin{remark}
The exponent $\alpha\geq 2$ identifies the nonlinearity of the pair potentials. When $\alpha=2$, the interaction between particles follows the \emph{Hooke's law}, which models a repulsive harmonic spring. The force between two overlapped particles linearly depends on the overlap. Hooke's model is one of the simplest contact model, which is still of great importance. From equation \eqref{eq: first derivative} and \eqref{eq: second derivative}, it is clear that the first derivatives of the Hooke potential are continuous at $\delta_c = 0$, while the second derivatives are discontinuous there.
When $\alpha > 2$, we obtain the so-called \emph{Modlin-Hertzian model} or \emph{Hertzian model}, which is based on the Hertz contact theory \cite{Johnson:1987}. The force between two overlapped particles has a power law nonlinearity.
Both the first and second derivatives of the Hertzian potential are continuous if $\alpha >2$. We take $\alpha = 5/2$ in the numerical examples of this paper.
\end{remark}

\subsection{Macroscopic Variables}

The average stress tensor can be calculated by the following formula,
\begin{equation} \label{eq: def of stress tensor}
\sigma=\frac{1}{V}\sum_{c \in C}\vec{l}_c\otimes \vec{f}_c,
\end{equation}
where $V$ is the volume of the containing domain and $C$ is the collection of all the contact points in the domain. According to equation \eqref{eq: def of stress tensor}, the pressure of the system can be defined as,
\begin{equation}
p=\frac{\sum_{\alpha=1}^d\sigma_{\alpha\alpha}}{d}
\end{equation}
where $d$ is the dimension of the domain. The bulk modulus can then be defined as,
\begin{equation}
B=\phi\frac{dp}{d\phi}.
\end{equation}

For a simple shear \cite{Ogden:1997} such that the deformation matrix is given by
\begin{displaymath}
D=\left(
    \begin{array}{cc}
      1 & \gamma \\
      0 & 1 \\
    \end{array}
  \right),
\end{displaymath}
the shear stress is $-\sigma_{12}$, and the shear modulus can be calculated by:
\begin{equation}\label{eq: cal shear modulus}
G = -\frac{d\sigma_{12}}{d\gamma}
\end{equation}

In addition, the average contact number or coordination number, denoted by $\bar{z}$, is a key geometric parameter in granular simulation. For a granular system consisting of $N$ particles,  $\bar{z}$ can be simply defined by,
$$
\bar{z}=\frac{|C|}{N}
$$
where $|C|$ is the number of all contact points.

We note that there exist particles with zero contacts, namely, the so called "rattlers". In addition, in our 2D frictionless circular disk system, some particle are called "dynamic rattlers" \cite{Fatih:2010} if their numbers of contacts are less than 3, which can lead to mechanical instability. Thus, these rattlers are excluded from the calculation of the coordination number. We denote the number of particles with at least 3 contacts by $N_3$, and the set of contact points for particles with at least 3 contacts by $C_3$. The modified definition of coordination number is,
\begin{equation}
\bar{z}=\frac{|C_3|}{N_3}
\label{eq:z3}
\end{equation}
\subsection{Jamming formation}

The jamming transition is an important physical process in granular system \cite{Majmudar:2005, O'Hern:2003, Zhang:2005}, which refers to the critical slowing down of the system due to "overcrowded" particles. Close to the transition point, the geometric constraints prevent the system to explore the phase space, and the system changes its behavior from gas or liquid like to solid like. During the transition, the system has to adjust its configurations more and more in a collective manner instead of via local movements, in order to achieve energy equilibrium. This corresponds to the essential difficulty for the numerical simulation which will be discussed later.

In this paper, we adopt the approach by O'Hern et al. \cite{O'Hern:2003} to generate a jamming configuration: An initial configuration is sampled with sufficiently low packing fraction $\phi_{\mathrm{init}}$ such that there is no overlapping particles after relaxation (energy minimization). The total energy and pressure for the relaxed system with volume fraction $\phi_{\mathrm{init}}$ is $0$. Starting from $\phi_{\mathrm{init}}$, we repeatedly increase the volume fraction $\phi$ with a small increment $\delta \phi$ by enlarging the size of particles accordingly, and then carry out energy relaxation. At some critical fraction $\phi_{\mathrm{cr}}$, there will be no free space for particles to explore and unavoidably, they will come into contact. As the system is further compressed, some particles overlap, and the total energy and pressure become nonzero and starts to grow. The system will also posses nonzero bulk modulus since the pressure increases upon compression. $\phi_{\mathrm{cr}}$ is called the critical volume fraction or jamming point. In  2D bi-disperse system with number ratio 1:1 and size ratio 1.4:1, $\phi_{\mathrm{cr}}$ is approximately 0.842 \cite{Krijgsman:2016}.

\def\mathrm{app}{\mathrm{app}}
\def\mathrm{pre}{\mathrm{pre}}
We can adopt a bisection algorithm to increase the accuracy of $\phi_{\mathrm{cr}}$. Let $\phi_{\mathrm{cr}}^\mathrm{app}$ be an approximate critical fraction associated with a configuration with nonzero energy and pressure, the previous unjammed configuration with zero energy and pressure has the volume fraction $\phi_{\mathrm{cr}}^{\mathrm{pre}}:=\phi_{\mathrm{cr}}^\mathrm{app} - \delta \phi$. We can then compress the unjammed configuration at $\phi_{\mathrm{cr}}^{\mathrm{pre}}$ by increasing its volume fraction to $\phi_{\mathrm{m}} = \phi_{\mathrm{cr}}^{\mathrm{pre}} + \delta \phi/2$, and see if the jamming transition already happens at $\phi_{\mathrm{m}}$. This procedure can be repeated until for example, $\delta \phi< \varepsilon$, where $\varepsilon$ represents the desired accuracy. We note that this procedure can achieve very high accuracy for the determination of the jamming point, which is not feasible for experimental methods.

The overall procedure is summarized in the following algorithm for the generation of the jamming configuration. Denote by $\min E(\phi, \C_{\mathrm{ref}})$ the energy minimization in the configuration space $\mathcal{C}$ with the packing fraction $\phi$ and the initial configuration $\C_{\mathrm{ref}}$, which can be solved by certain optimization algorithm in Section \ref{sec:numerics}. The minimizer is denoted by
$\mathcal{C}_{\mathrm eq} = \argmin E(\phi, \C_{\mathrm{ref}})$.
Here we use a square shaped container $\Omega$ and the periodic boundary condition for the granular system.
\defp_{\mathrm{ref}}{p_{\mathrm{ref}}}
\def\C_m{\mathcal{C}_m}
\def\C_{\mathrm{ref}}{\mathcal{C}_{\mathrm{ref}}}
\defE_{\mathrm{eq}}{E_{\mathrm{eq}}}
\defp_{\mathrm{eq}}{p_{\mathrm{eq}}}

\renewcommand{\algorithmicrequire}{ \textbf{Input:}}
\renewcommand{\algorithmicensure}{ \textbf{Output:}}

\begin{algorithm}[H]
\caption{Generation of the Jamming Configuration.}
\label{alg:jamming}
\begin{algorithmic}[1]
\REQUIRE ~~\\
particle number $N$;\\
initial packing fraction $\phi_{\mathrm{init}}$;\\
increment $\delta \phi$; \\
prescribed accuracy $\varepsilon$.
\ENSURE ~~\\
critical volume fraction $\phi_{\mathrm{cr}}$;\\
jamming configuration $\C_{\mathrm{jam}}$ with $\phi_{\mathrm{cr}}$.
\STATE Generate an initial configuration $\C_{\mathrm{init}}$ with $N$ particles and volume fraction $\phi_{\mathrm{init}}$;
\STATE let $\C_{\mathrm{ref}} = \argmin (\phi_{\mathrm{init}}, \C_{\mathrm{init}})$, with total energy $E_{\mathrm{ref}}=0$ and pressure $p_{\mathrm{ref}}=0$.
\label{alg:jam:step1}
\STATE // Increment Step
\WHILE{$E_{\mathrm{ref}}=0$ and $p_{\mathrm{ref}}=0$}
\STATE Let $\phi = \phi_{\mathrm{ref}} + \delta \phi$;\\
\STATE Fix the position of the particles, and enlarge their radii uniformly to generate a new configuration $\C_{\mathrm{ref}}$ with volume fraction $\phi$;\\
\STATE Let $\C_{\mathrm{eq}}(\phi) = \argmin E(\phi, \C_{\mathrm{ref}})$,  compute $E(\C_{\mathrm{eq}})$ and $p(\C_{\mathrm{eq}})$.
\IF{$E(\C_{\mathrm{eq}})>0$ \& $p(\C_{\mathrm{eq}})>0$}
\STATE  $\phi_{\mathrm{cr}}^\mathrm{app} = \phi$, $\C_{\mathrm{jam}} = \C_{\mathrm{eq}}$, break;
\ELSE
\STATE  $\phi_{\mathrm{ref}} = \phi$, $\C_{\mathrm{ref}} = \C_{\mathrm{eq}}$;
\ENDIF
\ENDWHILE
\STATE  // Bisection Step
\WHILE{$\phi_{\mathrm{cr}}^\mathrm{app} - \phi_{\mathrm{ref}} > \epsilon$ }
\STATE  $\phi_{\mathrm{m}}= (\phi_{\mathrm{cr}} + \phi_{\mathrm{ref}})/2$;
\STATE  Generate an intermediate configuration $\C_m$ with volume fraction $\phi_{\mathrm{m}}$, starting from $\C_{\mathrm{ref}}$.
\STATE  Let $\C_{\mathrm{eq}} = \argmin (\phi_{\mathrm{m}}, \C_m)$, and calculate $E(\C_{\mathrm{eq}})$ and $p(\C_{\mathrm{eq}})$.
\IF{$E(\C_{\mathrm{eq}})>0$ \& $p(\C_{\mathrm{eq}})>0$}
\STATE   $\phi_{\mathrm{cr}}^\mathrm{app} = \phi_{\mathrm{m}}$;
\ELSE
\STATE $\phi_{\mathrm{ref}} = \phi_{\mathrm{m}}$, $\C_{\mathrm{ref}} = \C_{\mathrm{eq}}$;
\ENDIF
\ENDWHILE

\STATE Denote $\phi_{\mathrm{cr}} = \phi_{\mathrm{cr}}^\mathrm{app}$, $\C_{\mathrm{jam}} = \C_{\mathrm{eq}}$.
\RETURN $\phi_{\mathrm{cr}}$, $\C_{\mathrm{jam}}$;
\end{algorithmic}
\end{algorithm}

\section{Numerical Method}
\label{sec:numerics}

\newcommand{\br}[1]{{(#1)}}

The minimization of the total potential energy \eqref{eq:Etot} relies on the state-of-the-art optimization and precondition techniques which will be presented in this section. The performance of those methods is contingent upon the particular applications. In particular, the construction of a good preconditioner, can be regarded as "an art rather than a science" \cite{Wathen:2015}.

\subsection{Energy minimization}

In this paper, we use the limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) method \cite{Nocedal:1980} and Fletcher-Reeves conjugate gradient (FR-CG) method \cite{Nocedal:2006} for the energy minimization of granular system. We denote the total energy of a system as $f(x)$ where $x$ is the positions of all particles. If $x^\br{k}$ is the iterate at step $k$, we denote the energy at step $k$ by $f^\br{k} = f(x^\br{k})$, and the gradient at step $k$ by $g^\br{k}=\nabla f^\br{k}$.

\subsubsection{L-BFGS method}
L-BFGS is one type of the quasi-Newton methods, namely, it utilizes an approximation of the inverse Hessian matrix to generate the search direction. Instead of the full  dense $n\times n$ approximations to the inverse Hessian ($n$ is the number of variables, which is $Nd$ with $N$ number of particles and $d$ the dimension in our problem), L-BFGS uses a low rank approximation with only a few vectors of length $n$ to represent the approximation implicitly.

We denote the approximate inverse Hessian at step $k$ by $B^{(k)}$. Assuming we have stored the last $m$ updates of the form $s^\br{k}=x^\br{k+1}-x^\br{k}$ and $y^\br{k}=g^\br{k+1}-g^\br{k}$, we define $\displaystyle\rho^\br{k}=\frac{1}{y^{(k),T}s^\br{k}}$.
The L-BFGS search direction $p^{\br{k}} (= B^{(k)}\nabla f^\br{k})$ can be obtained using the following algorithm.

\begin{algorithm}[H]
\caption{L-BFGS two loop recursion.}
\label{alg:L-BFGS}
\begin{algorithmic}[1]

\STATE $q=g^{(k)}$\;

\FOR{ $i=k-1,k-2,...,k-m$}
\STATE $\alpha^\br{i}=\rho^\br{i} s^{(i),T} q$\;
\STATE $q=q-\alpha^\br{i} y^\br{i};$
\ENDFOR
\STATE $\boxed{r=B^{(k)}_0q}$\;

\FOR{$i=k-m,k-m+1,...,k-1$}

\STATE $\beta=\rho^\br{i}y^{(i),T}r$\;
\STATE $r=r+s^\br{i}(\alpha^\br{i}-\beta)$\;
\ENDFOR

\RETURN $p^\br{k} = r (= B^\br{k}\nabla f^\br{k})$\;
\end{algorithmic}
\end{algorithm}

The matrix $B^{(k)}_0$ in the boxed step in the algorithm is an rough approximation to the inverse of the exact Hessian, and an effective choice is:
$$
B^{(k)}_0=\frac{s^{(k-1),T}y^\br{k-1}}{y^{(k-1),T}y^\br{k-1}}I
$$
This choice ensures that the search direction is well scaled and therefore the unit step length is accepted in most iterations. Besides, the diagonal matrix makes it much simpler to compute the multiplication $r=B^{(k)}_0q$.

\subsubsection{FR-CG method}

We use the Fletcher-Reeves variant \cite{Fletcher:1964} of the nonlinear conjugate gradient (CG) method \cite{Pytlak:2008}. As an extension of the linear conjugate gradient method \cite{Hestenes:1952}, the search direction $p^{(k)}$ at $k$th step of FR-CG method is defined by,

\def\mathrm{FR}{\mathrm{FR}}
\begin{algorithm}[H]
\caption{FR-CG.}
\label{alg:FRCG}
\begin{algorithmic}[1]

\STATE $\displaystyle\beta_{\mathrm{FR}}^{\br{k}} = \frac{g^{\br{k}, T}{g^{\br{k}}}}{g^{\br{k-1}, T}{g^{\br{k-1}}}}$;
\RETURN $p^{\br{k}} = -g^{(k)} + \beta_{\mathrm{FR}}^{\br{k}} p^\br{k-1}$.
\end{algorithmic}
\end{algorithm}

\subsubsection{Line Search}

Once we calculate the new search direction using either L-BFGS in Algorithm \ref{alg:L-BFGS} or FR-CG in Algorithm \ref{alg:FRCG}, the next iteration is given by
\begin{equation}
x^\br{k+1}=x^\br{k}+\alpha^\br{k}p^\br{k},
\end{equation}
where $\alpha^\br{k}$ is chosen by a line search method which satisfies Wolfe conditions or strong Wolfe conditions \cite{Nocedal:2006}, such that the update is stable.

\subsection{Preconditioning}
Preconditioning is important for the efficiency of the minimization algorithms, especially for large scale problems. We have to choose the preconidtioner matrix $P^\br{k}$ such that it is similar to the exact Hessian, and at the same time easy to invert.

\subsubsection{Preconditioned minimization algorithm}
For L-BFGS algorithm \ref{alg:L-BFGS}, we replace the boxed step in the algorithm with:
\begin{equation}\label{eq: preconditioned step}
\boxed{r = P^{(k),-1}q}
\end{equation}
in order to obtain a preconditioned L-BFGS algorithm.

For FR-CG algorithm \ref{alg:FRCG}, we have the following preconditioned version for $p^{\br{k}}$.
\def\mathrm{FR}{\mathrm{FR}}
\begin{algorithm}[H]
\caption{Preconditioned FR-CG method.}
\label{alg:PFRCG}
\begin{algorithmic}[1]
\STATE $y^{\br{k}} = P^{\br{k},-1} g^{\br{k}}$;
\STATE $\displaystyle\beta_{*\mathrm{FR}}^{\br{k}} = \frac{y^{\br{k}, T}{g^{\br{k}}}}{y^{\br{k-1}, T}{g^{\br{k-1}}}}$;
\RETURN $p^{\br{k}} = -y^{(k)} + \beta_{*\mathrm{FR}}^{\br{k}} p^\br{k-1}$.
\end{algorithmic}
\end{algorithm}

The equation \eqref{eq: preconditioned step} or Step 1 in Algorithm \ref{alg:PFRCG} is equivalent to solve a linear system:
\begin{equation}\label{eq:solvpsys}
P^\br{k}r=q.
\end{equation}

An alternative point of view for the preconditioning is to make a change of coordinates: $\widetilde{x}= P^{(k),1/2}x$. Considering the function $F(\widetilde{x})=f(P^{(k),-1/2}\widetilde{x})$, we have
\begin{equation}
\nabla F(\widetilde{x})=P^{(k),-1/2}\nabla f(P^{(k),-1/2}\widetilde{x}).
\end{equation}

Applying the L-BFGS algorithm \ref{alg:L-BFGS} or FR-CG algorithm \ref{alg:FRCG} to optimize the transformed function $F(\tilde{x})$, we can obtain the preconditioned version of the corresponding algorithm.

\subsubsection{Construction of preconditioners for large scale molecular system}

\def\mathrm cut{\mathrm cut}
\def\mathrm stab{\mathrm stab}

Preconditioning is well established in numerical linear algebra and numerical PDE problems \cite{Wathen:2015}. However, in many applications general purpose preconditioners do not work particularly well, and specifically designed preconditioned are preferable. For large scale atomic/molecular simulations, Packwood et. al. proposed the so-called "universal preconditioner" $P$ \cite{Packwood:2016}. They have successfully applied this preconditioner to atomic simulation and electronic structure calculation of crystalline system of the size from hundreds to $10^4$ atoms, with a speedup of order 1 or 2 \cite{Mones:2018}. In this paper, we will test the performance of those universal preconditioners for granular systems.

 The universal preconditioner $P$ is defined via the quadratic form
\begin{equation}\label{eq: precond matrix entry}
\begin{aligned}
	u^T P u &= \mu \sum_{0<|r_{ij}|<r_{\mathrm cut}} c_{ij} |u_i - u_j|^2, \\
        c_{ij} &= \exp(-A (r_{ij}/r_{nn} -1)).
        \end{aligned}
\end{equation}
Here $r_{\mathrm cut}$ is the cutoff radius, $\mu$ is the energy scale,
$A$ is a parameter which should be large enough to ensure that the nearest neighbor interactions dominant. The numerical experiments in \cite{Packwood:2016} and in our numerical results indicate that as long as $A$ is of order 1, it does not significantly affect the performance of the preconditioner. In particular, if we take $A=0$, we choose the preconditioner matrix $P$ as the stabilized adjacency matrix of the particles \cite{Packwood:2016}, namely,
\begin{equation}\label{eq: precondition matrix}
\begin{aligned}
P_{ij} &= \left\{
    \begin{aligned}
            -\mu, & &  d_{ij}< r_{\mathrm cut} \\
            0, & &  d_{ij}\ge r_{\mathrm cut}
    \end{aligned},
   \right.\\
P_{ii} &= -\sum_{i \ne j}P_{ij} + \mu C_{\mathrm stab},
\end{aligned}
\end{equation}
where $d_{ij}$ is the distance between particle $i$ and $j$, $C_{\mathrm stab}$ is a stabilization term to make sure the matrix is strictly positive definite, we choose $C_{\mathrm stab}=0.1$ in our application.

In our model problem, the bi-disperse granular system is composed of two types of particles, whose radii are denoted by $r_{\min}$ and $r_{\max}$. We take the cutoff radius $r_{\mathrm cut} > 2r_{\max}$ so that all the possibly overlapping particles are covered. For example, we can take  $r_{\mathrm cut} = 2.2r_{\max}$. The preconditioner matrix will be recalculated when the maximum displacement of particles is beyond $r_{cut}/4$, therefore it will not be recalculated frequently. The energy scale $\mu$ is chosen to make sure that the precondition matrix is of the same order as the real hessian, so we may choose the unit step-length when the L-BFGS method is applied. One way to obtain $\mu$ is,
\begin{equation}
\nu^{T}\left(\nabla E(x_{0}+\nu)-\nabla E(x_{0})\right)=\mu \nu^{T}P_{\mu=1}\nu,
\end{equation}
where $P_{\mu=1}$ is the matrix when $\mu=1$, and $\nu$ is a test perturbation with the following form
\begin{equation}
\nu(x)=M\left(\text{sin}(x/L)\right),
\end{equation}
where $L$ is the lengths of the periodic cell and $M$ is a constant, here we choose $M=10^{-3}r_{cut}$.

Since we need to solve linear system \eqref{eq:solvpsys} with respect to $P$, when the particle number is small, we use direct method and rearrange the index of grains by sparse reverse Cuthill-McKee ordering method \cite{Cuthill:1969}, the band width of $P$ can be reduced so that the equation \eqref{eq:solvpsys} can be easily solved. If the number of particles is large, iterative solvers such as preconditioned CG or AMG \cite{Brandt:1985} can be utilized.

\section{Numerical Results}
\label{sec:numerical experiment}

In our numerical experiments, we use the bidisperse granular system which consists of two types of particles with number ratio 1:1 and  size ratio $r_{\max}/r_{\min} = 1.4$, so the system will not crystallize \cite{Speedy:1999}. The particles are contained in a unit square domain $\Omega = [0,1]^2$. We use the periodic boundary condition in order to reduce the boundary layer effect. Solid wall boundary conditions can also be implemented.

In the numerical experiments, we test the performance of preconditioned L-BFGS method and preconditioned FR-CG method, by running several benchmark simulations: the jamming formation, the scaling law behavior for the jammed configuration, and the shear deformation of the over-jammed states. In those different scenarios, we will show that the proposed preconditioned algorithm has better performance compared to their un-preconditioned counterparts, and preconditioned L-BFGS is currently the method of choice.

\subsection{Jamming formation}

We first run Algorithm \ref{alg:jamming} to generate the jamming configuration, starting from a relaxed configuration with volume fraction $\phi< \phi_{\mathrm{cr}}$, where $\phi_{\mathrm{cr}}$ is the critical volume fraction (or jamming fraction). We let the particles grow with small increment $\delta\phi$, then we use either the preconditioned or un-preconditioned version of Algorithm \ref{alg:L-BFGS} and Algorithms \ref{alg:FRCG}, \ref{alg:PFRCG} as the energy minimizer in Algorithm \ref{alg:jamming}.

 In Figure \ref{fig: gradient vs iteration all method} , the residual norm of each iteration during a relaxation process is plotted. We apply four methods (L-BFGS, preconditioned L-BFGS, FR-CG, Preconditioned FR-CG) to minimize the energy of an unjamming configuration. A high accuracy tolerance is adopted here when an accurate jamming point is pursued. The figure shows that: L-BFGS is better than FR-CG, preconditioned L-BFGS method has the best convergence curve compared to other methods, and preconditioner provides a factor of two speed up for L-BFGS.

\begin{figure}[H]
  \centering

  \includegraphics[width=0.55\linewidth]{gradient_norm_vs_iter_all_method.pdf}\\
  \caption{$\ell^\infty$ norm of gradient vs iteration number for one energy minimization step of a granular system with particle number 4096, $\phi=0.8389$, the radii of particular grow by $r \rightarrow (1+\delta)r, \delta=5\times 10^{-5}$.}\label{fig: gradient vs iteration all method}
\end{figure}

We also notice that convergence curves for all methods are extremely rough and have long asymptotic regimes, that means the energy landscape is very complex and there are many local minimizers to explore.

Figure \ref{fig: jamming progress} shows the evolution of the granular system up to the jamming point. The volume fraction of the initial configuration is small, therefore particles can move around in order to achieve a zero energy configuration. When the volume fraction increases, the particles start to contact with each other, but initially it is still possible to achieve a zero energy configuration through energy relaxation (minimization). As the free space becomes less and less, at a critical volume fraction, the particles are forced to contact and overlap, and the potential energy of the system becomes non-zero.

\begin{figure}[H]
\centering
\subfigure[]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=0.9\linewidth]{p1.png}
\end{minipage}%
}%
\subfigure[]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=0.9\linewidth]{p28.png}
\end{minipage}%
}%
\subfigure[]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=0.9\linewidth]{p50.png}
\end{minipage}
}%
\centering
\caption{Schematic of the evolution from non-jamming to jamming state of bi-disperse granular system using Algorithm \ref{alg:jamming}, \textbf{a} Initial non-jamming configuration. \textbf{b} Intermediate configuration. \textbf{c} Jamming configuration, \emph{Visualization tool: OVITO \cite{Stukowski:2009}}.}\label{fig: jamming progress}
\end{figure}

\subsection{Scaling Law.}
Once the jamming configuration $\C_{\mathrm{jam}}$ is obtained, we can continue increasing the volume fraction in a quasi-static manner. It has been proposed by physicists that the granular system has the so-called scaling law behavior close to the jamming point \cite{O'Hern:2003}. Here, we numerically justify this scaling law using our energy minimization techniques, and more interestingly, we observe that the iteration number of the energy minimization methods also has similar scaling law behavior.

\begin{figure}[H]
  \centering

  \includegraphics[width=0.55\linewidth]{p_vs_Phi.pdf}\\
  \caption{pressure $p$ vs $\phi-\phi_{\mathrm{cr}}$ for a 2D bi-disperse granular system with particle number 4096.}\label{fig:pvsphi}
\end{figure}

In the Figure \ref{fig:pvsphi}, we plot the log-log curve of pressure $p$ vs $\phi-\phi_{\mathrm{cr}}$. It is clear that $p$ follows a power law with respect to $\phi-\phi_{\mathrm{cr}}$ close to the jamming point, and the critical exponent depends on the nature of interactions. For harmonic $(\alpha=2)$ and Hertzian $(\alpha=2.5)$ potentials, the critical exponent are approximately $1.0$ and $1.5$, respectively. In addition, this scaling law between $p$ and $\phi-\phi_{\mathrm{cr}}$ is robust for different initial configurations (they may have different $\phi_{\mathrm{cr}}$).

\begin{figure}[H]
  \centering

  \includegraphics[width=0.55\linewidth]{G_vs_Phi.pdf}\\
  \caption{Shear modulus $G$ vs $\phi-\phi_{\mathrm{cr}}$ for a 2D bi-disperse granular system with particle number 4096.}\label{fig:GvsPhi}
\end{figure}

If we exert small simple shear strain $\gamma$ to these configurations, we can compute their shear modulus using equation \eqref{eq: cal shear modulus}. In Figure \ref{fig:GvsPhi}, we plot shear modulus for different volume fractions. It shows that the shear modulus $G$ of the granular systems also follows a similar power scaling law with respect to $\phi-\phi_{\mathrm{cr}}$ but with different exponents, which is approximately $0.5$ when $\alpha=2.0$, and $1.0$ when $\alpha=2.5$.

\begin{figure}[H]
  \centering

  \includegraphics[width=0.55\linewidth]{z_vs_Phi.pdf}\\
  \caption{Average contact number $\bar{z}$ vs $\phi$ for a 2D bi-disperse granular system with particle number 4096. $z_0$ is the average contact number of the granular system with volume fraction $\phi_{\mathrm{cr}}$. In our 2D simulation, particles are frictionless, $z_0=4.0$.}\label{fig:zvsPhi}
\end{figure}

The average contact number $\bar{z}$ defined in \eqref{eq:z3} also exhibits a scaling law in regards to volume fraction difference $\phi-\phi_{\mathrm{cr}}$. In Figure \ref{fig:zvsPhi}, the lines in the log-log plots have similar slopes close to $0.5$ for both $\alpha=2.0$ and $\alpha=2.5$. The critical average contact number $z_0=4$ is the average contact number when a system reaches exactly the jamming state. The average contact number of a system takes a leap from zero to 4 when jamming transition occurs.

\begin{figure}[H]
  \centering
  \includegraphics[width=0.6\linewidth]{iter_vs_phi_mean_variance.pdf}\\
  \caption{Iteration numbers (average and variance) with respect to $\phi-\phi_{\mathrm{cr}}$. Red circles and blue stars are average iteration numbers of 10 quasi-static simulations using preconditioned L-BFGS method and un-preconditioned L-BFGS method, respectively. In addition, we plot the error bars for those simulations.}
  \label{fig:itervsphi}
\end{figure}

In Figure \ref{fig:itervsphi}, we plot the iteration number vs. volume fraction difference $\phi - \phi_{\mathrm{cr}}$ for the preconditioned L-BFGS and L-BFGS methods. Starting from certain jamming configuration, we compress the granular system by increasing the volume fraction $\phi$. Then for each relaxed over-jammed configuration, we exert small shear strains from $10^{-6}$ to $10^{-5}$ quasi-statically, in 10 steps. The iteration numbers exhibit "random" behavior, especially close to the jamming point. We plot both the average and the variance (as error bars) of the iteration number. We observe that the iteration number of each numerical method also follows a scaling law, which ranges from $10^3-10^4$ close to the jamming point ($\phi-\phi_{\mathrm{cr}}=10^{-6}$), to about hundreds when $\phi-\phi_{\mathrm{cr}} =10^{-2}$. The result shows that the preconditioned algorithm performs better, and also it is more robust to rounding error as the variance of the iteration number for the preconditioned algorithm is also smaller away from the jamming point.

In Table \ref{tab:methods}, we illustrate the performance of four optimization methods: L-BFGS, preconditioned L-BFGS, FR-CG,  and preconditioned FR-CG in three different situations given an over-jammed configuration with 4096 particles and volume fraction $\phi=0.8438$.
\begin{itemize}
\item Case 1, fixing the positions of the particles, and enlarging them slightly with ratio $1+\delta, \delta=5\times 10^{-5}$;
\item Case 2, applying a small shear strain ($\gamma=10^{-4}$) to the configuration;
\item Case 3, exerting a small ($\mathcal{O}(10^{-4})$) random perturbation to the configuration.
\end{itemize}

The comparison of iteration number and computational time in Table \ref{tab:methods} clearly demonstrates that the preconditioner can improve the efficiency of the energy minimization methods. Also, it seems we can simply take the parameter $A=0$ \eqref{eq: precond matrix entry} in the preconditioner.

\newcommand{\tabincell}[2]{\begin{tabular}{@{}#1@{}}#2\end{tabular}}
\begin{table}[H]
  \centering
    \begin{tabular}{|c|c|c|c|c|c|c|}
      \hline
      \multicolumn{2}{|c|}{Method}&L-BFGS&\tabincell{c}{P-L-BFGS\\($A=0$)}&\tabincell{c}{P-L-BFGS\\($A=3$)}&FR-CG&P-FR-CG\\
      \hline
      \multirow{2}*{Case 1}&n\_iter&2003&705&813&4574&1101\\
      \cline{2-7}
      &time/s&24.3&11.8&16.2&58.9&18.5\\
      \hline
      \multirow{2}*{Case 2}&n\_iter&2749&986&1351&7306&2013\\
      \cline{2-7}
      &time/s&29.2&15.8&27.1&90.1&33.5\\
      \hline
      \multirow{2}*{Case 3}&n\_iter&602&380&432&940&528\\
      \cline{2-7}
      &time/s&7.5&6.8&9.1&11.9&9.2\\
      \hline
    \end{tabular}
  \caption{Iteration number and computational time for L-BFGS,  P-L-BFGS, FR-CG, and P-FR-CG method for three different cases ($A$ is the parameter in \eqref{eq: precond matrix entry}). The granular system has 4096 particles and volume fraction $\phi=0.8438$. } \label{tab:methods}
\end{table}

\subsection{Shear Test}

In the shear test, we deform the over-jammed configuration with pure shear (stretching in x direction and compressing in y direction, while keeping the area unchanged). The deformation matrix is given by,
$$
D=\left(
    \begin{array}{cc}
      1+\gamma & 0 \\
      0 & \frac{1}{1+\gamma} \\
    \end{array}
  \right).
$$

\begin{figure}[H]
  \centering

   \includegraphics[width=0.9\linewidth]{pureshear_Fabric_Stress.pdf}\\
  \caption{Stress ($\sigma$) vs. strain and fabric ($F$) vs. strain curve for pure shear, $\gamma$ is the shear strain. The subscripts iso and dev, represent the isotropic and deviatoric parts of these tensors (stress, fabric), respectively. Let $\lambda_1$ and $\lambda_2$ be the eigenvalues of those tensors, the isotropic part is defined as $(\lambda_1+\lambda_2)/2$ and the deviatoric part is $|\lambda_1 - \lambda_2|/2$} \label{fig:StressStrainCurvePureShear}
\end{figure}

The stress vs. strain and fabric vs. strain curve for the pure shear simulation is shown in Figure \ref{fig:StressStrainCurvePureShear}. Notice that we use increments of different sizes, $4\times 10^{-4}$, $2\times 10^{-4}$, and $10^{-4}$ to reach a final 2\% shear strain, we still obtain qualitatively (quantitatively for small strain) similar curves.

\section{Conclusion}\label{sec:conclusion}
In this paper, we introduce the energy minimization techniques for the efficient simulation of dense athermal granular systems in two dimensions. Preconditioners are used to enhance the performance of the simulation. We carry out numerical experiments for some typical scenarios of granular simulation in order to validate the numerical methods, which include the jamming formation, scaling law phenomena close to the jamming point, and deformation of over-jammed states. Speedup of the preconditioned minimization methods are observed.

This work opens avenue for several possible improvements: First of all, we are going to extend the current work to three dimensions, which is physically more relevant to real applications, and still difficult for experimental studies \cite{Kou:2017, Zhang:2014} .

Secondly, from the practical point of view, we will optimize the implementation of preconditioners using, for example, AMG \cite{Brandt:1985} and other types of state-of-the-art techniques, especially for the three dimensional case.

Last but not least, we observe that the iteration number of our current preconditioned algorithms follows a scaling law behavior close to the jamming point, which is a fundamental bottleneck for the energy minimization approach. It remains open whether one can find a preconditioner which is robust close to the jamming point. Usually, an efficient preconditioner takes account of the long wavelength modes of the system. However, close to the phase transition point, more and more high frequency information might be needed. The development of numerical techniques in this direction relies on the understanding of the physical origin of the jamming transition.

\section*{Acknowledgments}
We thank Jie Zhang and Zhaohui Jin for stimulating discussions on granular modeling and experiments. We thank Christoph Ortner and Mingjie Liao for stimulating discussions on numerical methods and their efficient implementations.

\medskip
\bibliographystyle{aims}