Nicholas Nell

Nathalie Mitton
email: nathalie.mitton@inria.fr

Thomas Niesler

Riaan Wolhuter
email: wolhuter@sun.ac.za2

An Ad-Hoc Sensor Network for Vineyard Monitoring

Keywords: Ad-hoc Networking, Wireless Sensor Networks, Short-Range Communication, Low-Power Operation, Vineyard Monitoring, Precision Agriculture

This paper explores the use of an ad-hoc wireless sensor network in a vineyard-based environment to increase the effective range of sensor nodes without the use of high power communication techniques. Physical sensor and gateway nodes were developed in order to test the use of an adapted AODV routing protocol in a wireless sensor network. The network was deployed to a vineyard environment where the effectiveness of the strategy was evaluated based on the ability of the network to self-configure. First tests indicate that the protocol performs sufficiently well to be used in vineyard environments for data collection. We conclude that an ad-hoc sensor network is a feasible approach to low-power sensing over wide areas, as required for precision viticulture.

I. INTRODUCTION

The collection and analysis of vineyard telemetry data can significantly contribute to viticultural efficiency and productivity. Such telemetry data can, for example, be used to increase the quality and amount of grapes produced in a vineyard, while at the same time minimising the risk of wasteful operations. This can contribute to increased efficiency by reducing the cost of production and increasing environmental sustainability [START_REF] Wolhuter | The Softweb Solutions Website[END_REF]. Thus, an improvement in the effectiveness and ease of collecting such telemetry data will have a direct impact on the efficiency of the vineyard and allow methods of precision agriculture to be implemented.

Precision agriculture refers to the use of technology to make farming as efficient as possible by means of exact measurement. Generally, this is achieved by breaking farmland into smaller manageable sections, and then obtaining precise and continuous measurement data for each section. This allows for the determination of, for example, the variability in inter and intra-field condition of the plants and soil. If the exact conditions of every square meter of soil and every plant were known, the farmer could locally apply precisely the right quantity of water or fertiliser, thereby reducing waste and improving the quality of the produce. By using technology to virtually divide large areas of land into much smaller sections, the farmer is given far greater control over the crop. In the long term, this technology will also allow the collection and subsequent analysis of large sets of data, revealing trends that can be exploited for even more efficient farming and improved decision making [START_REF] Marios | Precision agriculture: Challenges in sensors and electronics for real-time soil and plant monitoring[END_REF].

Since these advances depend on the collection of telemetry data from the vineyards, a means must be developed to allow such data collection to be set up quickly, easily, at manageable cost and with reliable results. Currently, farmers must operate either without such data collection methods (especially in the case of small scale farms) or must use very labour-intensive approaches, such as on-site manual measurement and data collection. Where farmers have tried to deploy more advanced data collection methods, they have frequently been frustrated by the failure of traditional communication techniques. This could, for example, be due to a lack of cellular signal for GSM enabled nodes or power usage demands which cannot be met due to a lack of grid-tied power, especially when using high power communication techniques such as satellite or GSM. Other difficulties arise in vineyards which are situated in mountainous or hilly landscapes, where direct line-of-sight between sensor nodes in a network is difficult.

On-the-go data collection is employed in situations where wireless sensor networks are not easily implemented. A relatively inexpensive data collection technique, on-the-go data collection is based on a sensor attached to pre-existing mobile farming equipment, and where the data is collected while the equipment is in use. Examples include sensors attached to tractors or irrigation systems [START_REF] Adamchuk | On-the-go soil sensors for precision agriculture[END_REF]. A limitation associated with on-the-go sensors is that the equipment it is attached to has to be in use for data to be collected and that the data usually has to be manually retrieved by the farmer.

Many of these limitations can be overcome by using adhoc wireless sensor networks that use low power, shortrange communication techniques in conjunction with in-situ environmental, soil and plant sensors. Ad-hoc networks are especially useful because it is not necessary to configure each sensor node individually. Low-power designs allow the nodes to be battery-powered, and can thus be placed anywhere in a vineyard. The use of short-range communication requires sensor nodes to communicate with their nearest neighbours, and automatically discover routes to a gateway node without manual intervention. This approach offers a practically feasible means of setting up a data collection system that can provide autonomously gathered real-time data to the farmer. This paper will explore the design and implementation of such a low-power ad-hoc vineyard based wireless sensor network. A description of the network setup will be followed by the hardware design and implementation of the sensors. This is followed by a description of the network protocol used.

Finally, the results of the first practical test will be presented.

II. HARDWARE IMPLEMENTATION

Our sensor nodes consist of an MSP432p401r microprocessor (MCU) attached to several sensors and peripherals [START_REF]MSP432p401r Product Page[END_REF]. The use of an ultra-low-power processor allows the nodes to be battery-powered for extended periods.

Peripherals include a low-power RFM95W LoRa radio module, operating in the 868M Hz frequency band [START_REF]RFM95W Datasheet[END_REF]. This radio module is responsible for all inter-node communication.

Every node also includes a low-power Adafruit PA101D GPS module [START_REF]Adafruit Mini GPS Module[END_REF]. This module supports multiple positioning systems, including GPS, GLONASS, GALILEO and QZSS. Besides providing each node with positional data, the most important function of the GPS module is to provide a stable, shared time base to all nodes. The accurate pulse-per-second signal, which commences when a positional lock is established, allows the MCU to synchronise its onboard real-time clock (RTC). This is critical to the operation of the network, as will be discussed in Section V.

To enable data transfer to cloud storage from the sensor network, the root node includes a GSM module based on the GL865-Quad [START_REF]GSM-Click Product Page[END_REF]. The GSM module allows the gateway sensor node to periodically establish an HTTP connection to a web server and upload all received telemetry data.

Each node also includes a low-power BME280 temperature, pressure and humidity sensor. This integrated device simplifies the hardware design by avoiding the need for multiple separate sensors. For light intensity measurement, the MAX44009 light sensor is included, and to measure the volumetric water content of the soil, the EC-5 soil probe is used.

Lithium polymer (LiPo) batteries provide power to the nodes. This type of battery offers a high energy density and reliability. To regulate the voltage provided by the battery, a DC-DC regulator circuit was designed around the TPS63001 buck-boost converter [START_REF]TPS63001 Datasheet[END_REF].

All components were housed in an IP65 rated ABS enclosure to provide environmental protection and allow continuous operation in a vineyard, shown in Figure 1.

In the current prototype, the data collected by the network is uploaded to a web server running on a Raspberry Pi using the REST API.

III. NETWORK SETUP

Our sensor network consists of a mixture of edge nodes (nodes containing telemetry sensors) and gateway nodes (nodes that receive telemetry data from edge nodes). Edge nodes and gateway nodes have identical hardware, with the exception of the inclusion of the GSM module on the gateway nodes. The advantage of using the same underlying architecture for both node types is that it simplifies the maintenance of the hardware and makes it easy to increase the number of gateway nodes by simply adding a GSM module to an edge node.

Edge nodes primarily operate in a sleep state to conserve power. They are programmed to wake up periodically to collect telemetry data and to initiate a data transfer to a gateway node. Edge nodes also act as routing nodes when the need arises, for instance when there is no direct connection between an edge node and a gateway node. The network protocol will be discussed in greater detail in Section V.

In the current prototype, only a single gateway node is used. However, the network allows for the addition of multiple gateway nodes with no modification to the MAC or network layer protocols.

IV. MAC LAYER IMPLEMENTATION

A. S-MAC

The MAC layer protocol used by our network is an adaptation of the S-MAC protocol. S-MAC is specifically designed to reduce the power consumption of a wireless sensor network. It has been shown to perform well in situations where nodes have single tasks, are deployed in an ad-hoc fashion and where nodes spend most of their time in an inactive or sleep state [START_REF] Ye | An energy-efficient mac protocol for wireless sensor networks[END_REF]. The implemented protocol uses the GPS pulse-persecond signal to synchronise all network node in time. The nodes can then use time slots to synchronise communication. The use of time slots reduces the chance of message collisions on a shared communication medium. However, the main reason for the use of a time-slotted approach is to reduce node energy consumption. Careful synchronisation allows nodes to bring the onboard radio module out of its sleep state only within a specified time slot.

The implemented MAC protocol borrows a few key techniques from the S-MAC protocol. Neighbouring nodes share their transmission time slots, eliminating the need for nodes to listen to the shared communication channel in every time slot. Thus nodes can remain in a sleep state while not in their own, or a neighbouring node's, data transmission time slot. Another borrowed feature is that nodes will place their radio modules in a sleep state when overhearing a data message for which they are not the intended recipient. This happens when a node listens to the shared communication channel, expecting a message, but overhears a message with a different destination. It can then conserve power by entering a sleep state for the rest of the transaction. Synchronised global transmission time slots are used to allow the nodes of the network to perform the initial routing setup when there are no known routes or neighbours.

B. Datagram Structure

Messages sent using the MAC protocol consist of a header, followed by a data section and ended with a network statistics section as seen in Table I. The message header contains information about the route the packet is on, how many hops it has gone through, information about the sending node's transmission slot and current slot as well as a message type flag. The message type flag indicates to the receiving node how to interpret the data contained in the data section of the message. Finally, the network statistics section of the message is used only to analyse the performance of the network and is not necessary for network operation.

The message header structure is illustrated in Table II. The Next Hop byte contains the ID of the next node along the route to the destination node, whose ID is contained in the Network Destination byte. The Local Source byte contains the ID of the node from which the message was most recently passed, while the Network Source byte contains the ID of the node from which the message originated. The two source ID's will be the same at the start of the lifecycle of a message, when it is generated. The Tx Slot contains two bytes that indicate the transmission slot number of the local source node, while the Current Slot contains two bytes which indicate to the receiving node in which slot the message was sent by the local source node. This information is used to ensure that the network remains synchronised.

The data section of the message can be interpreted in one of three ways, depending on the value of the Message Type byte in the message header. First, it can be interpreted as telemetry data, thus containing all the sensor data dispatched by the source node. Second, the data can be interpreted as a route request and thirdly a route reply. The route reply and route request messages will be further discussed in Section V.

V. NETWORK PROTOCOL

A. AODV

AODV (Ad-hoc On-Demand Distance Vector Routing) was used to achieve routing within our network [START_REF] Perkins | Ad-hoc on-demand distance vector routing[END_REF]. The advantages of AODV include loop-free route discovery, quick repair of broken links, and avoiding the need for individual nodes to store information about the entire network. AODV makes use of route request (RREQ) and route reply (RREP) messages to establish paths within the network. A node wanting to find a route to another node (in this implementation a gateway node), will broadcast an RREQ packet to all its neighbouring nodes, (A neighbouring node is classified as any node that can receive a radio message from the specified node). The nodes that receive the RREQ will check to see whether they have a valid and up-to-date route to the destination node. If they do, they unicast an RREP packet back to the source node. If they do not, they rebroadcast the RREQ to all of their neighbours, and store reverse-path information indicating the way to the source node if an RREP packet is later received. Once the RREQ packet reaches its destination, or reaches a node that knows a valid route to the destination, an RREP packet is unicast back to the source node along the reverse path stored by each of the intermediate nodes along which the RREQ had travelled. Each of these nodes in turn set up forward path information to the destination node. Each node along the path from the source node to the destination node will now have a route to the destination node.

B. AODV Modifications

Some simplification to the AODV protocol is possible based on the needs of the overall system. Gateway nodes will never have to initiate a route discovery, and will only receive data from edge nodes. Edge nodes only need to store routes to gateway nodes, and not to other edge nodes, as an edge node will never be the final destination of a data message. Nodes must however store information about all of their neighbours, including their TX slots and their ID. This is because a node will not know whether it is an intermediate node along a route from another edge node to a gateway node until it receives a data message that needs to be forwarded along the route. When an edge node initiates a data transfer to a gateway node, it will look up the route to the gateway node in its routing table, and then send the message in the appropriate transmission slot. The source node will wait for an acknowledgement message from the node on the next hop along the route, and will then assume the message has successfully been sent. This eliminates the need for acknowledgement messages to be routed back from the destination gateway node. If an acknowledgement is not received from the next-hop node, the transaction will be reattempted in the next available global transmission slot. This is repeated until either an acknowledgement is received or the maximum number of permitted retries per message has been reached. Once this maximum is reached, a dead route is assumed and deleted from the sending nodes routing table. Then, a new route discovery sequence is initiated.

Each route recorded by an edge node has an associated expiration timer. The route will have to be checked by the edge node when the timer expires to ensure that the route is still active. The use of an expiration timer on a route allows the quick discovery of broken routes. The edge node updates its routing information by sending out a new RREQ with which to discover a new route to the gateway node.

VI. EXPERIMENTAL SETUP

For initial testing, several statistics regarding the performance of the sensor network were collected. These include the number of times a node had to initiate a route discovery, the message overhead that is required for successful data transfer from an edge node to a gateway node, and the time it takes a node to discover a route to a gateway node.

In our experiments the nodes are configured to collect and transmit sensor data once every twenty minutes, as the environmental data does not change rapidly. The gateway node will thus attempt to upload the data collected by all nodes at the start of the twenty-minute window. This is to reduce the risk of a node sending a data message during the GSM upload. If the GSM upload fails, the gateway will store the data messages that have not been uploaded, and retry the upload halfway through the twenty-minute window. Each of the GSM upload slots is limited to a maximum of one minute, to reduce the chance of messages arriving during the upload. When the gateway is in the process of uploading GSM data, it is unable to receive any messages over the radio channel due to the processor not being able to service both modules at the same time. This can be overcome in future by implementing a dedicated co-processor for the GSM communication.

A total of five nodes were installed in a hillside vineyard in Stellenbosch, South Africa. They were arranged in a star configuration so that each node can see the gateway node. The node placement is shown in Figure 2. This was done to test the loop-free nature of the routing protocol, as it is still possible for neighbouring nodes to respond to RREQ messages even if the node doing the route request has direct contact with the gateway node.

The gateway node was installed and activated first, after which the remaining nodes were enabled one by one. After a node is powered on, it waits for a GPS lock to synchronise the onboard RTC. Once the RTC is initialised, nodes will start the route discovery process. In each global transmission slot (global transmission slots are open every 10 seconds), each node is given a 20% chance to send an RREQ, if no routes to the gateway node are known, otherwise, if at least one route to the gateway is known, the edge nodes will listen in the global transmission slot. This is to reduce congestion at the start of the network's life, or when multiple nodes are activated and start the route discovery process at the same time.

After the nodes were installed this way, they were left for 48 hours to both record vineyard telemetry and to monitor the stability of the network by means of the data that was uploaded to the webserver.

VII. EXPERIMENTAL RESULTS

A. Initial Observations

As expected, it was found that the nodes that were activated directly after the gateway node found routes to the gateway node quickly and with a small number of RREQ messages. For example, The first node to be activated after the gateway (Node AB) found a route in 21 seconds with only 1 RREQ message. The second node (node DE) to be activated found a route to the gateway node after 31 seconds, also using a single RREQ. The remaining two nodes experienced the situation described below.

In situations where two edge nodes can both see the gateway, but one has a much stronger link to the gateway than the other, a direct route to the gateway is not always found. For example, for nodes AB and BC in Figure 2, AB sends an RREQ that is serviced with an RREP by both the gateway node and BC. Due to the capture effect present in LoRa demodulation, the source node AB would often only hear the stronger RREP from BC. This would cause AB to set up a route to the gateway node through node BC. This situation causes unnecessary traffic on the network, since a direct route to the gateway is in fact possible. To address this, a priority system was introduced, where all nodes other than gateway nodes would perform a brief carrier sense on the channel before sending an RREP. This allows the gateway node to always send an RREP first if it is in range.

Another addition that improved the route discovery rate of the network, is the introduction of periodic HELLO messages sent out by the gateway node. When an edge node receives a HELLO message and does not yet know a route to the gateway, it can add a route to its routing table. The edge node that received the HELLO message then does not have to initiate a route discovery. It also further prevented the situation described above, because when a node hears a HELLO message from the root, it can replace any other route to the gateway node with a direct link, if the node had stored a previously inefficient route.

Another situation that led to incorrect or inefficient routing data, was the hidden node problem. In this situation, two nodes hear an RREQ and both respond to it with an RREP. However, the two nodes are not able to hear each other. This can cause the source node to receive incorrect routing data.

B. Network Results

After the changes mentioned above, all nodes were reset. The network performed noticeably better, with all edge nodes now finding the most direct route to the gateway node. A summary of the initial route discovery statistics is shown in Table III. Each edge node required only a single RREQ to discover a route to the gateway node, except for node BC. This node heard a HELLO message before it began route discovery and therefore did not have to perform this sequence. The gateway node sent out three RREQ messages, further confirming that it received each of the RREQ messages.

After the initial discovery, the edge nodes sent telemetry data to the gateway node once in every twenty-minute time window. The expiration time on routes was set to sixty minutes and did not reset even when a data message was successfully sent. This was done to test the situation where a node has to rediscover a route to the gateway node. This was found to occur successfully.

After 10 hours of operation, during which edge nodes performed route discovery every hour, the edge nodes had sent 131 data messages to the gateway, and 148 messages in total. This demonstrates that even with the artificially high incidence of broken links precipitated by not resetting the expiration timers the network protocol overhead was only 12.9%. The 148 messages were sent within a single one second time slot each, with 3540 global transmission slots available in the same 10 hour period. It is clear that with only 4 edge nodes and one gateway node, the network is far from capacity.

C. GSM Upload

During the testing period, we observed that it took the gateway node required up to 46 seconds to upload the five stored messages to the webserver. This delay occurs because the GSM module has to perform one push request to the web server for each message, and response time from the webserver is dependent on multiple uncontrollable factors. Thus, the blocking nature of the GSM uploads is a key limitation in the number of edge nodes a gateway node can accommodate in the current implementation. We have presented the design of an ad-hoc sensor network intended for real-time viticultural monitoring that uses an adapted AODV routing protocol to allow easy installation without prior configuration. Once placed in the vineyard and enabled, nodes independently discover routes to the gateway node quickly and without excessive routing overhead. The telemetry data that was gathered by each node was successfully sent to the gateway node and uploaded to a web server, where it could be viewed on a dashboard [START_REF]ThingsBoard -Dashboard[END_REF].

The upload of data using the GSM module could be improved by better scheduling to avoid the transmission slots of neighbouring nodes, or by offloading the task to a dedicated co-processor to avoid blocking any inter-node data traffic.

The network can be improved by adding additional priority levels to nodes during the route discovery process, where nodes that are closer to the gateway node (and therefore have routes with fewer hops), have a higher priority in sending RREP messages back along a path to a source node that initiated an RREQ. Nodes with a higher priority would then have to perform carrier sensing on the channel for a shorter duration of time compared to nodes with lower priority levels, allowing them to transmit RREP messages sooner.

Fig. 1 .

 1 Fig. 1. Node hardware in enclosure.

Fig. 2 .

 2 Fig. 2. Node placement in the vineyard.

ACKNOWLEDGEMENT

We thank Albert Strever and Talitha Venter for arranging access to the test vineyard. This research was supported by the Wine Industry Network of Expertise and Technology (Winetech). We thank Telkom South Africa for additional support.


R Priam

Probabilistic Elastic Embedding Model: Comparison of Alternative Models *

In data visualization, Elastic Embedding adds an exponential penalty to an Euclidean criterion. It is able to separate the natural classes but its lacks a probabilistic generative setting which brings more flexibility to the modeling and the inference. Hence, it is proposed a new generative interpretation of Elastic Embedding which is closely related to LargeVis. Numerical experiments compare the proposed model and several alternative ones via two new visual indicators among different approaches.

Introduction

In data visualization, the observations in the available sample are vectorial: the rows or the columns of numerical data matrices or tables. A family of methods is dedicated to the symmetric numerical matrices which contain the distances or similarities between high-dimensional data vectors. An extensive literature exists and diverse approaches have been developped until today in this domain of research but they can appear in very different domains of research even if they are closely related: visualization, word embedding, statistical analysis and clustering analysis. Currently, the method t-Distributed Stochastic Neighbor Embedding (t-SNE or tSNE) [START_REF] Van Der Maaten | Visualizing high-dimensional data using t-sne[END_REF] is the state of art for large datasets, it is non generative and scaled well for large datasets. For a generative/probabilistic basis, it can be cited for instance Probabilistic Principal Component * A nearly similar version of this document was sent to review previously from 2018, and after a shorter version without section 3.3 was accepted at a conference in 2020 but canceled.

† rpriam@gmail.com

Analysis [START_REF] Tipping | Probabilistic principal component analysis[END_REF], Glove [START_REF] Pennington | Glove: Global vectors for word representation[END_REF], and LargeVis [START_REF] Tang | Visualizing large-scale and high-dimensional data[END_REF] with a very good scaling. Nextafter, this section presents the notation, the purpose and the plan of the paper.

Data notation and reduction: Let have the available high-dimensional data as a set of data vectors in a space with M dimensions as follows:

X = {x i ∈ R M ; 1 ≤ i ≤ N} .
Let define W = (w i j ) from the distance between x i and x j , for instance d i j = x i -x j , with w i j = d i j , or more generally a function of d i j . Typically, a weighted nearest neighboors graph such as a heat matrix with for instance w i j = e -0.5d 2 i j /τ for τ > 0 when d i j is enough small, and w i j = 0 otherwise. A weighted graph G = (V, E,W ) is defined from V , E, and W which stands respectively for the set of vertices, edges and weights: an edge e i j comes from from a pair of vertices (v i , v j ) in the graph of nearest neighboors. It is also denoted Ē for the set of pairs of vertices that are not neighboors.

The purpose of a reduction is to summarize X by finding relevant lower dimensional representations:

Y = {y i ∈ R S ; 1 ≤ i ≤ N} .
Generally S = 2 as the visualization appears in the two dimensional plane, even if three dimensions or even more remain possible. Existing models consider an embedding of low dimensional positions via their pairwise distances/similarities plus bias/intercept terms. The parameterization approximate the true distances d i j between data pairs (x i ,x j ) via the distances between pairs of S-dimensional lower position vectors (y i ,y j ) denoted δ i j . Let also 1 denote the Poisson mass functions as follows, ϕ PSN (., δ ) = 1 .! (.) e δ e -e δ .

Purpose of the paper: Among the existing non parametric methods, the Elastic embedding [START_REF] Carreira-Perpiñan | The elastic embedding algorithm for dimensionality reduction[END_REF] has been shown to perform well on several datasets such that it may be extended for even better results. EE is defined via a penalized Eudlidean criterion which aims at retrieving the true pairwise distances between each couple of data through new variables which represent the data in a lower dimensional space. Parametric probabilistic models for visualization are generally more flexible than ad hoc criterion and the probabilities have an intuitive interpretation, hence introducing such framework for EE is of main interest.

Plan: Section 2 presents EE, its probabilistic version, links with the current literature and its illustration with real data. Sections 3 presents an analytical evaluation of the model under some hypothesis with a new framework for the comparisons. Section 4 concludes the paper with future perspectives.

Generative Elastic Embedding

In this section, Elastic Embedding is reviewed before introducing a new parameterization in order to deduce a probabilistic model for visualization.

Elastic Embedding

The model EE was introduced as an alternative of SNE [START_REF] Hinton | Stochastic neighbor embedding[END_REF] for dimensionality reduction. The author Miguel A. Carreira-Perpinan notices a link between SNE [START_REF] Hinton | Stochastic neighbor embedding[END_REF] and Laplacian eigenmaps (LE) [START_REF] Belkin | Laplacian eigenmaps for dimensionality reduction and data representation[END_REF]. The SNE criterion is rewritten into two parts: a) the first term is exactly equal to the LE objective except that normalised affinities would be considered, and b) the second term encourages latent points to be far away with a log of a sum of exponential quantities. The author explains that the SNE criterion induces that the data vectors in a same neighboor are mapped into a same area of the low dimensional view and also separates every projections points instead of focusing on the separation of the clusters themself.

In the more formal setting, when w + i j and w - i j stands for given weights, and δ i j = y i -y j 2 , the objective function of EE is as follows,

C EE (Y , γ) = ∑ i, j w + i j δ i j + γ ∑ i, j
w - i j exp(-δ i j ) .

For choosing the input values for the two sets of weights w + i j and w - i j , many approaches are possible as explained by the author of the method.

-In the original paper of EE it is choosen the following one:

w + i j = e -0.5d 2 i j /τ w - i j = d 2 i j if i = j 0 otherwise.
.

It is underlined in the research paper on EE that other weighting are possible such as sparser graphs.

-For sparser weights, it can be noticed that for x i and x j enough near the weight w - i j vanishes, hence summing for only non neighboors pairs for the original space leads to a very similar criterion. In the same idea for x i and x j enough far away the weight w + i j vanishes hence summing for only neighboors pairs leads also to a similar criterion. These approximations have for corresponding criterion CEE with new weights w+ i j and wi j as follows,

w+ i j = w + i j if j ∈ N i 0 otherwise w- i j = w - i j if j ∈ N i 0 otherwise .
This difference of modeling for each part of the graph E and Ē recalls recents models in the literature, in particular LargeVis. A limit of the approximation is to not separate well projection in a given cluster but for very large dataset, this problem may be limited because viewing separately each point appears not relevant for a human and only zooming at a given area of the map might suffer from this approximation. Some authors add an additional pernalty for the local projection but this is not studied herein and only this criterion is considered hereafter, with same name EE.

Generative version

For a generative setting of EE, it is proposed to rewrite the approximated criterion denoted CEE with the constant 1 C w . The weights are supposed w + i j integer and w - i j binary, such that for α > 0 small,

-CEE (Y , γ) -C w = + ∑ (i, j)∈E w + i j log e -δ i j +log α -γ ∑ (i, j)∈ Ēw - i j e -δ i j -C w ≈ + ∑ (i, j)∈E log{e -δ i j +log α } w + i j -αw + i j e -δ i j -log w i j + ! +γ ∑ (i, j)∈ Ē log{e -e -δ i j } = log ∏ (i, j)∈E ( δi j ) w i j e -δi j w i j ! ∏ (i, j)∈ Ēe -δi j γ .
Hence, it is obtained,

argmin Y CEE (Y , γ) ≈ argmax Y LEE (Y , γ) , when, δi j = e -δ i j +log w - i j denoted δ - i j for Ē e -δ i j +log α denoted δ + i j for E
More formally, this is the behavior of the Poisson mass distribution associated to Euclidean distance which is in stake here, Let denote δ a real and w a positive integer. Then the logarithm of the mass distribution ϕ PSN (w, e -δ +log α ) behaves like the limit value -wδ for a real enough small α. This is a consequence of the logarithm of ϕ PSN with the introduced parameterization is equal to -wδ + w log ααe -δlog w!. This induces that only the quantity -wδ is in stake in the maximization and is minimized as wanted for EE. The proposed likelihood for EE is then,

Lα EE (Y , γ) = ∏ (i, j)∈E ϕ PSN (w + i j , δ + i j ) ∏ (i, j)∈ Ēϕ PSN (0, γ δ - i j ) .
Note that this is enough general to apply to some other methods for visualization in order to define a generative version. This model can be seen as a mixture between the observed nodes and the non observed nodes with a classifying likelihood. It is named genEE and the inference for the parameters Y is discussed next subsection.

1 C w = -∑ (i, j)∈E log w + i jlog α∑ (i, j)∈E w + i j

Links with the literature

Several models are closely related to genEE as latent position models have a long history in the statistical and computational literature. The generative model for correspondence analysis [START_REF] Beh | Simple correspondence analysis: A bibliographic review[END_REF] is based on a Poisson with a parameterization e y T i y j +b i +b j . The scalar quantities aggregated in B = (b i ) i are fixed and equal to the logarithm of the margin sums normalized to one: the similarity is different to the Euclidean distance and less flexible than in Glove [START_REF] Mikolov | Distributed representations of words and phrases and their compositionality[END_REF]. Here δi j is written with a scalar product instead of the less flexible Euclidean distance. This model is also defined for a rectangular matrix of positive integer and was shown to lead to almost perfectly equal estimation of the latent position than the matricial correspondence analysis [START_REF] Benzecri | Correspondence Analysis Handbook[END_REF]. In genEE, a criterion very similar to LargeVis (LV) is recognized except that a Poisson is the foundation for the underlying distributional hypothesis instead of a Bernouilli distribution, as suggested in [START_REF] Priam | Symmetric generative methods and tsne: A short survey[END_REF]. The probabilistic interpretation with w i j an integer remains the product between the likelihood of the observed edges with the likelihood of the non observed edges with a weighting γ for regulating the importance of each distribution. As a remark, alternative distributions such as multinomial, negative poisson, binomial or even normal may be also possible. The main difference with genEE is that the weights are included in the model hence a fully probabilistic approach may be possible if the fitting is statistically relevant. The model genEE may be also able to exactly behave like LV when a suitable parameterization is chosen for the expression of the expectations.

Empirical illustration

In this section, the models are compared with a real dataset in order to discuss the quality of the clusters obtained on the nonlinear embedding by several methods. Several indicators are considered for the evaluation of the compactness and the distance between clusters.

-Data: The MNIST dataset is a subset of a larger 20x20 arrays with grey level and then centered in 28x28 arrays.

-Experimental settings: For the learning of the quantities Y in EE, a sequential algorithm is processed, this approach is recently widely studied. In particular, it is implemented the probabilistic edge sampling [START_REF] Tang | Visualizing large-scale and high-dimensional data[END_REF] and negative sampling [START_REF] Mikolov | Distributed representations of words and phrases and their compositionality[END_REF] for an approximate stochastic gradient descent algorithm without normalization. For tSNE, LargeVis and UMAP [START_REF] Mcinnes | Umap: Uniform manifold approximation and projection for dimension reduction[END_REF], the available packages in language python have been used.

-Empirical results: On the Figure 2, the 10 classes are almost perfectly separated on the map from the recents methods tSNE, LargeVis and UMAP with large frontiers for the two last ones. EE performs less well for this dataset with smaller frontiers around the clusters.

Next section is dedicated to explaining the differences between the results of the different methods.

Numerical comparison of the criteria

A study of the behavior of the criteria under simple hypothesis are presented in this subsection in order to help the understanding of their empirical difference in practice. A trade-off between the clusters compactness, clusters separation plus the quality of the global and local projection needs to be found in order to bring the best map. According to the indicators in this section, two proporties of the criteria a checked.

For the compactness of a cluster, the coefficient of variation (CV) has to be enough small. For the separation between the clusters, either the clusters are degenerated into their center either their compactness has an influence with the frontiers.

Distance between two clusters

Following [START_REF] Noack | An energy model for visual graph clustering[END_REF], the distance between two subgraphes can be modeled via the approximation that every nodes in each subgraph have same final positions.

Let denote the number of nodes in first subgraph n 1 and n 2 the number for the second subgraph, and the numbe of edges n e , with n e n 1 n 2 and f = n 1 n 2 /n e . The criterion has the following form:

u(δ ) = n e g(δ ) + n 1 n 2 h(δ ) .
The function g(.) and h(.) depends on the criterion in stake, several functions are presented in [START_REF]An energy model for visual graph clustering[END_REF] for older methods. Note that the function g(.) and h(.) have the following proporties: g is increasing and h is decreasing from a value δ > 0 hence h (δ ) > 0 and h (δ ) < 0 while the second order derivatives are supposed positive for insuring a solution. The functions g(.) and h(.) lead to competitive results when they lead to small clusters and large frontiers between the clusters. For the distances, this framework results into the solutions a) for [START_REF] Eades | A heuristic for graph drawing[END_REF] such as δ [START_REF] Davidson | Drawing graphs nicely using simulated annealing[END_REF] and δ = f for LinLog in [START_REF] Noack | An energy model for visual graph clustering[END_REF][START_REF]An energy model for visual graph clustering[END_REF]. This leads to justify the better results for the method LinLog with an empirical illustration with perfectly symmetric clusters where the frontiers was shown to be more apparent. For Elastic Embedding with full weigths, it may be found δ = log(γ f ). A limit of this criterion is when only the nearest neighboors are involved for the local projection, then the function g(.) may not anymore be involved too. A complement of this study is the influence of the components of the criteria to the compactness of the clusters which is underlined hereafter.

2 log x = f , b) δ = 3 √ f in [16], c) δ = 4 √ f in

Compactness of the clusters: indicator

For a better idea of how behave a criterion from the point of view of the visualization of clusters it may be observed the first moments of the part related to the local projection and the global projection. To our knowledge, this question has not been studied in the literature for the quality of a visualization in a general situation.

-Hypothesis: it is supposed the Gaussianity of the clusters with ellipsoidal shapes. The random variables in stake are,

δ i j = y i -y j 2
With G (µ, Σ) is for a Gaussian random variable with mean µ and variance Σ, the distribution is as follows,

y i ∼ G (µ, 1 2 Σ) .
When y i and y j are independent, y i -y j ∼ G (0, Σ) hence their difference has a similar distribution but with a zero expectation and a double variance. For diagonal matrices Σ = diag k (σ 2 k ) with non null elements σ 2 k for the variances per dimension, according to [START_REF] Kettani | On the distribution of the distance between two multivariate normally distributed points[END_REF], the distribution of δ i j is the sum of S random variables, say informally,

δ i j ∼ ∑ k Γ(0.5, 2σ 2 k ) .
In the simple case when σ k = 1, it is known that the Euclidean distance follows a χ 2 hence the variance is 2S and the mean is S. For the more general case, under the hypothesis of diagonal Gaussianity for y i , the expectation and variance of δ i j are respectively,

µ δ = E[δ i j ] = ∑ k σ 2 k σ 2 δ = Var[δ i j ] = 2 ∑ k σ 4 k .
It is retrieved the case of the χ 2 when it is replaced σ k by 1 with corresponding mean and variance.

-Coefficient of variation (CV): this is a measure of the dispersion of a probability distribution. It is written as the ratio of the standard deviation to the mean. A function with smaller coefficient of variation may lead to more compact clusters because the values of the function has less variations during the projection where the positions y i are constructed. The expression of the indicator for the measure of the scattering is as follows,

I CV = √ Var[g(δ i j )] E[g(δ i j )] ≈ (g (µ δ )) 2 σ 2 δ + (g (µ δ )) 2 σ 4 δ 4 g(µ δ )+ g (µ δ ) σ 2 δ 2 .
For a given function g(δ i j ), a functional expan- sion suggests an approximate value. This indicator is for instance estimated via simulation of Gaussian distributions for y i , and the computations of the resulting distance. It is computed via simulation of (1000 here) vectors for S = 2 for some Gaussian distribution, in the Table 1. Note that the function e atan( √ δ i j ) is also evaluated because of its derivative proportional to (1 + δ i j ) -1 which is very similar to the logarithmic case for log(1 + δ i j ). The table 1 presents the indicators for different shapes of projection: it may confirm the compacity of the clusters for LargeVis in comparison to other methods with a smaller indicator. The alternative functions with an even smaller indicator I CV are expected candidates to behave well at least for small datasets and well separated classes. Next a more general model is proposed for evaluating a criterion when the clusters have spherical shapes.

Compactness of the clusters: model

The generalized criterion of EE is approximated for binary weigths, with E = ∪ E , as follows,

C EE (Y , γ) = ∑ (i, j)∈E g(δ i j ) + γ h (i, j)∈ Ē(δ i j ) ≈ ∑ k ∑ (i, j)∈E g(δ i j ) + γ ∑ (i, j)∈ Ēh(δ i j ) ≈ ∑ n g(δ ) f + (δ )dδ + γ n e h(δ ) f -(δ )dδ .
Here it is clear that under an hypothesis of clustering, the distances δ i j are considered separately in each cluster for the positive interactions (left part) while an unique cluster is considered for the negative interactions (right part). As seen in the paragraph above, under Gaussianity, the densities f + and f -are the Gamma depending on the variances of δ i j , respectively in E and Ē.

The hypothesis of indendence between the distances may be strong but is taken nextafter. For comparison purpose, it is supposed only two clusters and all the variances per dimension equal to σ 2 hence in each cluster the distribution of δ i j is chosen:

f + ∼ Γ(1, 2σ 2 ) .
For modeling the distribution outside the clusters it is taken an ad hoc value equal to 2.1 σ 2 which is the double of one cluster plus some variability from the frontier such that the distribution involved is chosen:

f -∼ Γ(1, 4.2σ 2 ) .
Even if this choice is not optimal because the variance for the interactions may vary from a model to another, it leads to a first informative result. The optimization for a projection is reduced at finding the value of σ 2 which makes minimal the simplified criterion for several choice of function g(.) and h(.).

The first results with such model is very encouraging. The obtained curves of the criteria against the variable σ 2 are shown just after for Elastic Embedding (EE), LinLog (LL) and LargeVis (LV), when n 1 = n 2 = 50 and n 2 = 350 while γ = 1.

The proposed modeling of a criterion for visualization is able to show without any simulation of the related method how would behave the method in practice. In particular, LinLog seems more prone to construct projections with larger clusters. The result for LargeVis leads to an optimal variance higher than Elastic Embedding for this choice of parameters. The corresponding optimal values for σ 2 are approximatively after graph reading,equal to 0.4 for EE, 1.8 for LL and 1.2 for LV. Note that for LV, the optimum exists and for a larger interval of values the concavity of the curve leads to a clear maximization of the log-likelihood. The fourth methods is with the first part of EE and the second part of LV with its optimum around 0.8. In future, a comparison with different settings of the parameters and an improvement of the modeling for the interactions, attractions and repulsions is required for further informations. For instance, the shape for the macro cluster can be chosen ellipsoisal as its accounts for the two spherical clusters.

Conclusion and perspectives

In this article, it is shown that Elastic Embedding is directly related to LargeVis with a particular choice of sparse weighting and with a Poisson mass function instead of the Bernoulli one. The method is compared from different points of view with alternative ones in order to understand further their differences. Two new approaches are proposed as a complement of the estimation of the distance in a particular sit-uation when two clusters collapse into their centers with zero variances. In the two proposed approaches, the variances of the clusters are in stake in order to discuss the reasons why some methods lead to smaller clusters during the projection than other ones.

Figure 1 :

 1 Figure 1: Nonlinear embeddings for MNIST.

Function

  g or h I CV I CV I CV a = 1.0 a = 0.5 a = 2.0 b = 1.0 b = 1.0 b = 1.0 δ i j

Figure 2 :

 2 Figure 2: Values of the criteria against the variable σ 2 .

Table 1 :

 1 Indicators of compactness with Σ = diag(a, b).


Etienne Polack

Geneviève Dusson

Benjamin Stamm

Filippo Lipparini
email: filippo.lipparini@unipi.it

Étienne Polack

Grassmann extrapolation of

in performing an extrapolation of the density matrix stems from the non-linearity of the problem. In other words, a linear combination of idempotent density matrices is not an idempotent density matrix, as density matrices are elements of a manifold and not of a vector space. To circumvent this problem, strategies that extrapolate the Fock or ks matrix [START_REF] Pulay | Fock Matrix Dynamics[END_REF][START_REF] Herbert | Energy-Conserving Born-Oppenheimer Molecular Dynamics via Fock Matrix Extrapolation[END_REF] or that use orbital transformation methods [START_REF] Hutter | Exponential Transformation of Molecular Orbitals[END_REF][START_REF] Vandevondele | An Efficient Orbital Transformation Method for Electronic Structure Calculations[END_REF][START_REF] Vandevondele | Quickstep: Fast and Accurate Density Functional Calculations Using a Mixed Gaussian and Plane Waves Approach[END_REF] have been proposed.

A completely different strategy has been proposed by Niklasson and coworkers. [START_REF] Niklasson | Time-Reversible Born-Oppenheimer Molecular Dynamics[END_REF][START_REF] Niklasson | Extended Born-Oppenheimer Molecular Dynamics[END_REF][START_REF] Niklasson | Extended Lagrangian Born-Oppenheimer Molecular Dynamics with Dissipation[END_REF] In the extended Lagrangian Born-Oppenheimer (xlbo) method, an auxiliary density is propagated in a time-reversible fashion and then used as a guess for the scf procedure. The strategy is particularly successful, as it combines an accurate guess with excellent stability properties.

In particular, the xlbo method allows one to perform accurate simulations converging the scf to average values (for instance, 10 -5 in the root-mean-square (rms) norm of the density increment), which are usually insufficient to compute accurate forces. An xlbo-based bomd strategy has been recently developed by some of us in the context of polarizable multiscale bomd simulations of both ground and excited states. [START_REF] Loco | Hybrid QM/MM Molecular Dynamics with AMOEBA Polarizable Embedding[END_REF][START_REF] Loco | Towards large scale hybrid QM/MM dynamics of complex systems with advanced point dipole polarizable embeddings[END_REF][START_REF] Nottoli | Excited State Born-Oppenheimer Molecular Dynamics through a coupling between Time Dependent DFT and AMOEBA[END_REF][START_REF] Bondanza | Polarizable embedding QM/MM: the future gold standard for complex (bio)systems?[END_REF] Multiscale strategies can be efficiently combined, in a focused model spirit, to bomd simulations to extend the size of treatable systems. Using a polarizable embedding allows one to achieve good accuracy in the description of environmental effects, especially if excited states or molecular properties are to be computed. In such a context, the xlbo guessing strategy allows one to perform stable simulation even using the modes 10 -5 rms convergence threshold, which, thanks to the quality of the xlbo guess, typically requires only about 4 scf iterations. Recently, scf-less formulations of the xlbo schemes have also been proposed. [START_REF] Niklasson | Next generation extended Lagrangian first principles molecular dynamics[END_REF][START_REF] Niklasson | Density-Matrix Based Extended Lagrangian Born-Oppenheimer Molecular Dynamics[END_REF] Unfortunately, the performances of the xlbo-based bomd scheme are not so good when a tighter scf convergence is required, which can be the case when one wants to perform md simulations using post Hartree-Fock (hf) methods or for excited states described in a timedependent dft framework. [START_REF] Nottoli | Excited State Born-Oppenheimer Molecular Dynamics through a coupling between Time Dependent DFT and AMOEBA[END_REF][START_REF] Nottoli | An enhanced sampling QM/AMOEBA approach: The case of the excited state intramolecular proton transfer in solvated 3-hydroxyflavone[END_REF] In fact, such methods require the solution to a second set of qm equations which are typically non-variational, making them more susceptible to numerical errors and instabilities. Computing the forces for non-scf energies requires therefore a more accurate scf solution.

The present work builds on all previous methods for density matrix extrapolation and aims at proposing a simple framework to overcome the difficulties associated with the nonlinearity of the problem. The strategy that we propose is based on a differential geometry approach and is particularly simple. First, we introduce a molecular descriptor, i.e., a function of the molecular geometry and other molecular parameters that represents the molecular structure in a natural way that respects the invariance properties of the molecule within a vector space. At the (n + 1)-th step of an md trajectory, we fit the new descriptor in a least-square fashion using the descriptors available at a number of previous steps and obtain a new set of coefficients. However, we do not use them to directly extrapolate the density. Instead, we first map the unknown density matrix, that we aim to approximate, from the manifold where it is defined to its tangent space. We then perform the extrapolation to approximate the representative density matrix in the tangent space, before mapping this approximation back to the manifold in order to obtain an extrapolated density matrix that satisfies the required physical constraints. This geometrical strategy, that has recently been introduced in the context of density matrix approximation by us, [START_REF] Polack | An Approximation Strategy to Compute Accurate Initial Density Matrices for Repeated Self-Consistent Field Calculations at Different Geometries[END_REF] allows one to use standard linear extrapolation machinery without worrying about the non-linear physical constraints on the density matrix, since both the space of descriptors and the tangent space are vector spaces. As the mapping between the manifold and the tangent space is locally bijective, no concerns about redundant degrees of freedom (such as rotations that mix occupied orbitals) arise. The map and its inverse, which are known as Grassmann Logarithm and Exponential, are easily computed and the implementation of the strategy is straightforward. We shall denote this approach as Grassmann extrapolation (g-ext).

In this contribution, we choose a simple, yet effective molecular descriptor and, for the extrapolation, a least square strategy. These are not the only choices. As our strategy allows one to use any linear extrapolation technique between two vector spaces, which can be in turn coupled with any choice of molecular descriptor, more advanced strategies can be proposed, including machine learning. Our approach ensures that the extrapolated density, independent of how it is obtained, satisfies all the physical requirements of a density stemming from a single Slater determinant.

The paper is organized as follows. In the upcoming Section 2, we present all necessary theoretical foundations required for the development and implementation of the presented g-ext approach. Section 3 then presents detailed numerical tests illustrating the performance of the extrapolation scheme, including realistic applications of bomd within a qm/molecular mechanics (mm)-context before we draw the conclusion in Section 4.

Theory

We consider Born-Oppenheimer ab-initio bomd simulations where the position vector R ∈ R 3M evolves in time according to classical mechanics as

M i Ri (t) = F i (t, R(t)), (1)
where R i (t), F i (t) ∈ R 3 denote the position of the i-th atom with mass M i respectively the force acting on it at time t. We consider a general qm/mm-method but the setting also trivially applies to pure qm-models. The forces at a given time t and position R of the nuclei arise from different interactions, namely qm-qm, qm-mm and mm-mm interactions.

The computationally expensive part is to determine the state of the electronic structure, which is modelled here at the dft level with a given basis set of dimension N . Note that considering hf instead of dft would not change much in the presentation of the method. It consists of computing the instantaneous non-linear eigenvalue problem

             F R (D R )C R = S R C R E R C T R S R C R = Id N D R = C R C T R , (2)
where C R ∈ R N ×N and D R ∈ R N ×N denote the coefficients respectively of the occupied orbitals and density matrix and E R ∈ R N ×N the diagonal matrix containing the energy levels. Further, F R denotes the dft-operator acting on the density matrix and S R the customary overlap matrix.

At this point it is useful to note that the slightly modified coefficient matrix C R := S 1 /2 R C R belongs to the so-called Stiefel manifold defined as follows

St(N, N ) := V ∈ R N ×N V T V = Id N , (3)
due to the second equation in Equation (2). In consequence the normalized density matrix

D R = C R C T R = S 1 /2 R D R S 1 /2
R belongs to the following set

Gr(N, N ) := D ∈ R N ×N D 2 = D, D T = D, Tr D = N , (4)
which can be identified with the Grassmann manifold of N -dimensional subspaces of R N by means of the spectral projectors. In the following, we always assume that the density matrix has been orthonormalized, and therefore drop the ∼ from the notation. For any D ∈ Gr(N, N ), one can associate the tangent space T D which has the structure of a vector space.

The evolution of the electronic structure can therefore be seen as a trajectory t → D R(t) on Gr(N, N ) where t → R(t) denotes the trajectory of the nuclei.

The goal of the present work is to find a good approximation for the electronic density matrix at the next step of md trajectory by extrapolating the densities at previous steps.

More precisely, based on the knowledge of the density matrices Thus, the problem formulation can be seen as an extrapolation problem of the following form: given the set of couples (R(t i ), D i ) and a new position vector R(t n ), provide a guess for the solution D n . Here and in the remaining part of the article, we restrict ourselves on the positions of the qm-atoms, i.e., with slight abuse of notation we denote from now on by R the set of qm-positions only, even within a qm/mm-context.

D i := D R(t i ) , i = n - N t , . . . ,
In order to approximate the mapping R → D R , we split this mapping in several sub-maps that will be composed as follows:

R 3M → M → T D 0 → Gr(N, N ) R → d R → Γ R → D R = Exp D 0 (Γ R ), (5)
where the first line shows the concatenation of maps in terms of spaces and the second in terms of variables. The different mappings will be presented and motivated in the following.

The first map is a mapping of the nuclear coordinates R ∈ R 3M to a (possibly highdimensional) molecular descriptor d R ∈ M that accounts for certain symmetries and invariances of the molecule. The last map, known as the Grassmann exponential, is introduced in order to obtain a resulting density matrix belonging to Gr(N, N ) and thus to guarantee that the guess fulfils all properties of a density matrix. As Gr(N, N ) is a manifold this is not straightforward. The second mapping is the one that we aim to approximate but before we do that, let us first introduce those two special mappings, i.e., the molecular descriptor and the Grassmann exponential, in more details.

Molecular descriptors

The map R → d R is a map from atomic positions to molecular descriptors. These descriptors are used as fingerprints for the considered molecular configurations. Such molecular descriptors have been widely used in the past decades e.g., to learn potential energy surfaces (pes), [START_REF] Behler | Neural network potential-energy surfaces in chemistry: a tool for large-scale simulations[END_REF][START_REF] Rupp | Fast and Accurate Modeling of Molecular Atomization Energies with Machine Learning[END_REF][START_REF] Bartók | Machine-learning approach for one-and two-body corrections to density functional theory: Applications to molecular and condensed water[END_REF][START_REF] Behler | Constructing high-dimensional neural network potentials: A tutorial review[END_REF][START_REF] Manzhos | Neural network-based approaches for building high dimensional and quantum dynamics-friendly potential energy surfaces[END_REF][START_REF] Chmiela | Machine learning of accurate energy-conserving molecular force fields[END_REF][START_REF] Chmiela | Towards exact molecular dynamics simulations with machine-learned force fields[END_REF] or to predict other quantities of interest. Among widely used descriptors, one can find Behler-Parinello symmetry functions, [START_REF] Behler | Generalized neural-network representation of highdimensional potential-energy surfaces[END_REF] Coulomb matrix, [START_REF] Rupp | Fast and accurate modeling of molecular atomization energies with machine learning[END_REF] smooth overlap of atomic positions (soap), [START_REF] Goscinski | The role of feature space in atomistic learning[END_REF] permutationally invariant polynomials, [START_REF] Braams | Permutationally invariant potential energy surfaces in high dimensionality[END_REF] or the atomic cluster expansion (ace). [START_REF] Drautz | Atomic cluster expansion for accurate and transferable interatomic potentials[END_REF][START_REF] Bachmayr | Atomic Cluster Expansion: Completeness, Efficiency and Stability[END_REF] These molecular descriptors are usually designed to retain similar symmetries as the targeted quantities of interest.

In this work, the quantity we are approximating is the density matrix, which is invariant with respect to translations as well as permutations of like particles. The transformation of the density matrix with respect to a global rotation of the system depends on the implementation, as it is possible to consider either a fixed Cartesian frame or one that moves with respect to the molecular system. In the former case, there is an equivariance with respect to rotations of the molecular system, while in the latter, the density matrix is invariant. We should therefore in principle use a molecular descriptor satisfying those properties.

However, the symmetry properties we will rely on are mostly translation and rotation invariance. Therefore, we will use a simple descriptor in form of the Coulomb-matrix denoted by d R , given by

(d R ) ij =      0.5z 2.4 i if i = j z i z j R(t i ) -R(t j ) otherwise . ( 6
)
Note that such a descriptor is not invariant (nor equivariant) with respect to permutations of identical particles. However, we have found this descriptor to offer a good trade-off between simplicity and efficiency. Note that since we aim to extrapolate the density matrix from previous time-steps, permutations of identical particles never occur from one time-step to another and we do not need to rely on this property. Nevertheless, we expect that a better description could be achieved by using more flexible descriptors, such as ace polynomials or the soap descriptors, where the descriptors themselves can be tuned.

The Grassmann exponential

We only give a brief overview as the technical details have already been reported elsewhere. [START_REF] Polack | An Approximation Strategy to Compute Accurate Initial Density Matrices for Repeated Self-Consistent Field Calculations at Different Geometries[END_REF][START_REF] Edelman | The Geometry of Algorithms with Orthogonality Constraints[END_REF][START_REF] Zimmermann | Manifold Interpolation and Model Reduction[END_REF] The set Gr(N, N ) is a smooth manifold and thus, at any point, say D 0 ∈ Gr(N, N ) in our application, there exists the tangent space T D 0 such that one can asso- Exp D 0 (0) = D 0 . These mappings are not only abstract tools from differential geometry but can be computed by means of performing an singular value decomposition (svd). [START_REF] Polack | An Approximation Strategy to Compute Accurate Initial Density Matrices for Repeated Self-Consistent Field Calculations at Different Geometries[END_REF][START_REF] Edelman | The Geometry of Algorithms with Orthogonality Constraints[END_REF][START_REF] Zimmermann | Manifold Interpolation and Model Reduction[END_REF] In our application we use the same reference point D 0 in all cases which brings some computational advantages as will be discussed in more detail the the upcoming Section 2.3.

The approximation problem

Since the tangent space T D 0 is a (linear) vector space, we can now aim to approximate the mapped density matrix on the tangent space T D 0 . To simplify the presentation, we shift the indices in the following and describe the extrapolation method for the first N t time steps.

In the general setting, we should consider the positions R(t i ) for i = n -N t , . . . , n -1, to extrapolate the density matrix at position R(t n ), where n is the current time step of the md.

We look for parameter functions c i , such that, given previous snapshots Γ i = Log D 0 (D i ) for i from 1 to N t , corresponding to some R(t i )'s, the approximation of any density matrix on the tangent space is written as

R → Γ app (R) = Nt i=1 c R,i Γ i ∈ T D 0 , (7)
with

Γ i = Γ R(t i ) .
The question is then how to find these coefficient functions c R,i and we propose to find those via the resolution of a (standard) least-square minimization problem. For a given position R, we look for coefficients that minimise the 2 -error between the descriptor d R and a linear combination of the previous ones d R(t i )

min c R ∈R N t d R - Nt i=1 c R,i d R(t i ) 2 . (8)
In matrix form, this simply reads

min c R ∈R N t d R -P T c R 2 , ( 9
)
where P is the matrix of size N t × N d containing the descriptors P i,j := (d R(t i ) ) j . Note that we only fit on the level of the descriptor, i.e., the mapping from the position vector R to the descriptor d R , and that this method is similar to the ones used by Alfè 2 , Arias et al. [START_REF] Arias | Ab Initio Molecular-Dynamics Techniques Extended to Large-Length-Scale Systems[END_REF] , where the descriptors they used were the positions of the atoms and only considered the previous three time-steps of the md.

If the system is underdetermined, we select the vector c R that has the smallest norm.

However, in general, the system is overdetermined as we have more descriptors than snapshots. This implies that this formulation verifies the interpolation principle: for every i and j from 1 to N t , the solution of Problem ( 8) at the positions R(t j ) satisfies c R(t j ),i = δ ji .

In principle, should we consider a large amount of previous descriptors, then the system may become undetermined and violates the interpolation principle. To mitigate this, we can use a stabilization scheme, as explained in the upcoming subsection.

Note that once we have computed the coefficients c R by solving Problem (12), one computes the initial guess for the density by using the same coefficients in the linear combination on the tangent space as in Equation ( 7) and finally take the exponential (see Equation ( 5)).

The rational for this step is that, if the second mapping in Equation ( 5), that we denote here by F : M → T D 0 , was linear, then there would hold

F Nt i=1 c R,i d R i = Nt i=1 c R,i F(d R i ) = Nt i=1 c R,i Γ i . (10)
In practice, the mapping is however not linear and this approach works well in the test cases we considered. A possible explanation for this is the unfolding of the nuclear coordinates into a high-dimensional descriptor-space M. Indeed, the high-dimensionality of M seems to allow an accurate approximation of F by a linear map. Further, if the system is overdetermined, the scheme satisfies the interpolation property Γ j = Γ(R(t j )), and hence we recover the expected density matrix D R(t j ) = Exp D 0 (Γ j ).

Stabilization

To stabilize the extrapolation by limiting high oscillations of the coefficients, we apply a Tikhonov regularization

min c R ∈R N t   d R - Nt i=1 c R,i d R i 2 + ε c R 2   , (11)
for some choice of ε. This problem is always well-posed, and corresponds to solving the following problem

min c R ∈R N t d R -P T • c R 2 , (12)
where d R ∈ R N d +Nt is the vector d R padded with N t zeros and P ∈ R Nt × R N d +Nt is the P matrix padded with the square diagonal matrix ε Id Nt . We observe in practice that using such a stabilization makes possible to use more previous points without degradation of the initial guess.

The final algorithm

Given previous density matrices D R(t j ) for j = 1, . . . , N t , the initial guess is computed following Algorithm 1. That is, we start by computing the logarithms of the density matrices D R(t j ) , from the coefficients C R(t j ) that are first orthonormalized by performing

C R = S 1 /2 R C R .
Here, we remark that we assume that the density matrices D R(t j ) have been previously Löwdin orthonormalized.

We then compute the descriptors needed to build the P matrix and solve Problem (12). This provides the coefficients in the linear combination of the Γ i s on the tangent space.

Finally, we compute the exponential of the linear combination in order to obtain the predicted

); c ← LeastSquares(desc, p n ); Γ app ← n-1 i=n-1-k c(i ) • gmat(:, :, i ); C app ← Exp(cref, Γ app ); return 2 • C app • C T app ;
density matrix.

Note that the reference point D 0 is chosen once and for all, which makes the computations of these logarithms lighter, even though there is no theoretical justification for keeping a single point D 0 as a reference. Indeed, it is known that the formulae are only correct locally (around D 0 ) on the manifold. However, in practice we have never observed the need to change the reference point. This enables us to compute only one logarithmic map per time step; and hence, only two svd in total per time step. To have a robust algorithm that will work even in this edge case, it will be sufficient to check that the exponential and logarithmic maps are still inverse of one another.

Finally, to be on the safe-side with respect to the computations of the exponential, we have added a check on the orthogonality of the matrix that is obtained: If the residue is higher than a certain threshold, we then perform an orthogonalization of the result.

Numerical tests

In this section we present a series of numerical tests of the newly developed strategy. We test our method on four different systems. All the systems have been object of a previous or current study by some of us, and can therefore be considered representative of real-life applications.

The first system is 3-hydroxyflavone (3hf) in acetonitrile. [START_REF] Nottoli | An enhanced sampling QM/AMOEBA approach: The case of the excited state intramolecular proton transfer in solvated 3-hydroxyflavone[END_REF] Two systems (ocp and appa) are chromophores embedded in a biological matrix -namely, a carotenoid in the orange carotenoid protein (ocp) and flavine in acid phosphatase (appa), a blue light-using flavine photoreceptor. [START_REF] Bondanza | The Multiple Roles of the Protein in the Photoactivation of Orange Carotenoid Protein[END_REF][START_REF] Bondanza | Molecular Mechanisms of Activation in the Orange Carotenoid Protein Revealed by Molecular Dynamics[END_REF][START_REF] Hashem | From crystallographic data to the solution structure of photoreceptors: the case of the AppA BLUF domain[END_REF] The fourth system is dimethylaminobenzonitrile (dmabn) in methanol. [START_REF] Nottoli | Excited State Born-Oppenheimer Molecular Dynamics through a coupling between Time Dependent DFT and AMOEBA[END_REF] The main characteristics of the systems used for testing are recapitulated in Table 1. The systems used for testing include a quite large qm chromophore, the ocp and three medium-sized systems, embedded in large (appa, 3hf) and medium-sized environments (dmabn) and are representative of different possible scenarios.

To test the performances of the new g-ext strategy, we performed three sets of short (1 ps) multiscale bomd simulations on ocp, appa, 3hf, and dmabn. ks density functional theory was used to model the qm subsystem, using the B3LYP 38 hybrid functional and Pople's 6-31G(d) basis set. [START_REF] Hehre | Self -Consistent Molecular Orbital Methods . XII . Further Extensions of Gaussian -Type Basis Sets for Use in Molecular Orbital Studies of Organic Molecules Publishing Articles you may be interested in Self -consistent molecular or[END_REF] For the stability and energy conservation of the method, we did a longer and more realistic simulations of 10 ps on 3hf, where the flavone moiety was described using the ωB97X hybrid functional 40 and Pople's 6-31G(d) basis set. In all cases, the environment was modeled using the AMOEBA polarizable force field. [START_REF] Ponder | Current Status of the AMOEBA Polarizable Force Field[END_REF] All the simulations have been performed using the Gaussian-Tinker interface previously developed by some of us. [START_REF] Loco | Hybrid QM/MM Molecular Dynamics with AMOEBA Polarizable Embedding[END_REF][START_REF] Loco | Towards large scale hybrid QM/MM dynamics of complex systems with advanced point dipole polarizable embeddings[END_REF] In particular, we use a locally modified development version of Gaussian 42 to compute the qm, electrostatic and polarization energy and forces, and Tinker [START_REF] Rackers | Tinker 8: Software Tools for Molecular Design[END_REF] to compute all others contributions to the qm/mm energy. We implemented the g-ext extrapolation scheme in Tinker, that acts as the main driver for the md simulation, being responsible of summing together all the various contributions to the forces and propagating the trajectory. At each md step, using the GauOpen interface, 44 the density matrix, molecular orbital (mo) coefficients, and overlap matrix produced by Gaussian are retrieved.

These are used to compute the extrapolated density as described in Section 2. The density is then passed back to Gaussian to be used for the next md step. All the simulations were carried out in the NVE ensemble, using the velocity Verlet integrator and a 0.5 fs time step. Concerning stabilization, we found that good overall results were obtained using a parameter ε := 10 3 • r scf , where r scf is the tolerance of the scf algorithm.

Numerical results

To assess the performance of the g-ext guess we perform 1 ps md simulations on the four systems described in Section 3 starting from the same exact conditions (positions and initial velocities) and using various strategies to compute the guess density for the scf solver. We compare various flavors of the g-ext method with the the xlbo extrapolation scheme. [START_REF] Niklasson | Extended Born-Oppenheimer Molecular Dynamics[END_REF] Here, we note that the original xlbo method performs a propagation of an auxiliary density matrix, which is then used as a guess. The latter is not idempotent: to restore such a property, we perform a purification step at the beginning of the SCF procedure using McWeeny's algorithm. [START_REF] Mcweeny | Some Recent Advances in Density Matrix Theory[END_REF] In the following, we therefore compare our method, where we use 3 to 6 previous points for the fitting and extrapolation, to both the standard xlbo and to xlbo followed by purification (xlbo/mw). We use an scf convergence threshold of 10 -5 with respect to the rms variation of the density.

We report in Table 2, for each method, the average number of scf iterations performed along the md simulation together with the associated standard deviation. As the xlbo strategy requires 8 previous points, during which a standard scf is performed, we discard the first points from the evaluation of the aforementioned quantities to have a fairer comparison.

We do not report the total time required to compute the guess, as it is in all cases very small (up to 0.1 s wall clock time for the largest system using the g-ext(6) guess). This is an important consideration, as the g-ext method requires one to perform various linearalgebra operations (in particular, thin svd) that can in principle be expensive. Thanks to the availability of optimized LAPACK libraries, this is in practice not a problem.

Table 2: Performances of the g-ext method for different number of extrapolation points, compared with the xlbo algorithm with and without McWeeny purification. All the results were obtained using a 10 -5 convergence threshold for the root-mean-square increment of the density matrix and are derived from a 1 ps long md simulation, using a 0.5 fs time step. We report the average number of iterations required to converge the scf, together with the associated standard deviation. Note that the first 8 steps were discarded. effect on the performances of the xlbo method only for the largest system, ocp, where it results in the gain of almost one scf iteration on average. On the other systems, the purification step has a smaller effect.

ocp
In all the systems we tested, the performances of the g-ext method are systematically better than in xlbo, including with McWeeny purification. The effectiveness of the g-ext extrapolation increases when going from 3 to 6 points, but quickly stagnates. We have performed further tests with more than 6 (up to 20) extrapolation points, but never noted any further gain.

We observe a reduction in the number of iterations that goes from 0.5 in dmabn to 0.75 in ocp (1.62 when compared to xlbo without McWeeny purification). We remark that these gains, while apparently not so large, are greatly amplified during the md simulation, due to the large number of steps that need to be performed.

The tests performed with a 10 -5 convergence threshold are representative of a standard, dft ground state bomd simulation. When performing a more sophisticated quantum mechanical calculation, such as a bomd on an excited state pes, 18 such a convergence threshold may not be sufficient to guarantee the stability of the simulation, as the scf solution is used to set up the linear response equations and the numerical error can be amplified, resulting in poorly accurate forces.

We tested the g-ext algorithm in its best-performing version, the one that uses six extrapolation points, with a tighter, 10 -7 threshold, again for the rms variation of the density.

The results are reported in Table 3, where we compare the g-ext(6) scheme with the xlbo method with McWeeny purification.

The xlbo method is based on the propagation of an auxiliary density and therefore the accuracy of the guess it generates depends little on the accuracy of the previous scf densities.

As a consequence, its performances are reduced if a tighter convergence is required. The g-ext guess, on the other hand, uses previously computed densities as its building blocks and one can expect the accuracy of the resulting guess to be linked to the convergence threshold used during the simulation.

This is exactly what we observe. Using a threshold of 10 -7 , the g-ext(6) guess exhibits significantly better performances than xlbo, gaining, on average, from about 0.7 to about 3 scf iterations on the tested systems.

Stability

The good performances of the g-ext guess come, however, at a price, namely, the lack of time reversibility. We can thus expect the total energy in a NVE simulation to exhibit a long-time drift (ltd). Time reversibility and long-time energy conservation are, on the other hand, one of the biggest strengths of the xlbo method.

To investigate the stability of bomd simulations using the g-ext guess, we build a chal- lenging case, where we start a bomd simulation far from well-equilibrated conditions. We use the 3hf system as a test case and achieve the noisy starting conditions by starting from a well-equilibrated structure and changing the dft functional from B3LYP to ωB97XD. This way, we have a physically acceptable structure, with no close atoms or other problematic structural situations, but obtain starting conditions that are far from equilibrium.

We report in Figure 1 the total energy along a 10 ps bomd simulation of 3hf in acetonitrile using either a 10 -5 scf convergence threshold (left panel) or a 10 -7 one (right panel). The same results for a 10 -6 threshold are reported in the supporting information. We compare the g-ext(3) and g-ext (6) methods to the xlbo one including McWeeny purification. As already noted, while in principle the purification may spoil the time reversibility, this has no noticeable effect in practice.

The very noisy starting conditions are apparent from the energy profiles, that exhibits large oscillations in the first couple hundreds femtoseconds.

To better estimate the short-and long-time energy stability, we report in Table 4 the average short-time fluctuation (stf) and ltd of the energy. The former is computed by taking the rms of the energy fluctuation every 50 fs and averaging the results over the trajectory, discarding the first 500 fs, the latter by fitting the energy with a linear function and taking the slope.

All methods show comparable short-term stability, which is to be mainly ascribed to the chosen integration time-step. On the other hand, from both the results in Table 4 and Figure 1, we observe a clear drift of the energy when the g-ext method is used. In particular, the system cools of about 10 kcal/mol with either g-ext(3) or g-ext (6). The xlbo trajectory, despite the McWeeny purification, exhibits an almost perfect energy conservation.

These results are not surprising, but should be taken into account when choosing to use the g-ext guess, which, if coupled to a 10 -5 scf convergence threshold, cannot guarantee long-term energy conservation. The drift is overall not too large and can be handled by using a thermostat. Whether or not the trade between performances and energy conservation is acceptable for a production simulation is a decision that ultimately lies with the user.

Increasing the accuracy of the scf computation improves the overall stability for g-ext,

which is already good at 10 -6 and becomes virtually identical to the one offered by the xlbo method at 10 -7 .

Conclusion

In this contribution, we presented an extrapolation scheme to predict initial guesses of the density matrix for the scf-iterations within bomd. What makes our approach new is that we enforce the idempotency of the density matrix by extrapolating not the densities themselves, but their map onto a vector space, which is the tangent plane to the manifold of the physically acceptable densities. Such a map is locally bijective, so that after performing the As such, the technique presented in this paper can be seen as a simple case of a more general framework. Such a framework allows one to recast the problem of predicting a guess density by extrapolating information available from previous md steps as a mapping between two vector spaces, i.e., the space of molecular descriptors and the tangent plane. This geometric approach can be seen as an alternative to extrapolating quantities derived from the density, such as the Fock or Kohn-Sham matrix, as proposed by Pulay and Fogarasi 4 and by Herbert and Head-Gordon. [START_REF] Herbert | Energy-Conserving Born-Oppenheimer Molecular Dynamics via Fock Matrix Extrapolation[END_REF] However, the framework we developed, using molecular descriptors and a general linear extrapolation technique, can in principle be easily extended to such approaches.

That being said, our choices of both the molecular descriptor and of the extrapolation strategies are far from being unique. In recent years, molecular descriptors gained attraction within the rise of machine-learning (ml) techniques. Our choice, namely, using the Coulomb matrix, is only one of the many possibilities, and while being simple and effective, more ad-vanced descriptors may be used and possibly improve the overall performances of the method.

We also used a straightforward (stabilized) least-square interpolation of the descriptors at previous point to compute the extrapolation coefficients for the densities. This strategy is, again, simple yet effective. However, many other approaches can be used. In particular, ml techniques may not only provide a very accurate approximated map, but also benefit of a larger amount of information (i.e., use the densities computed at a large number of previous steps), further improving the accuracy of the guess. Improvements on the descriptors and extrapolation strategies are not the only possible extensions of the proposed method. A natural extension that is under active investigation is the application to the g-ext guess to geometry optimization, for which the xlbo scheme cannot be used.

Overall, even the simple choices made in this contribution produced an algorithm that exhibits promising performances. In all our tests, the g-ext method outperformed the wellestablished xlbo technique, especially for tighter scf accuracies which may be relevant for post-scf bomd computations, including computations on excited-state pes. While we tested the method only for ks dft, it can also be used for Hartree-Fock or semiempirical calculations. The main disadvantage of the proposed strategy with respect to the xlbo method is, however, the lack of time reversibility, which manifests itself as a lack of long-term energy conservation. In particular, for longer md simulations, the total energy may exhibit a visible drift, which is something that the user must be aware of. In our test, the observed drift was relatively small and the use of a thermostat should be enough to avoid problems in practical cases, however, this is a clear, and expected, limitation of the proposed approach.

We note that, using a tighter scf convergence, which is also the case where the proposed method shows its best performances, produces an energy conserving trajectory, even starting from very noisy conditions. A time-reversible generalization of the g-ext method is anyways particularly attractive, and is at the moment under active investigation.

  ciate nearby points D ∈ Gr(N, N ) to tangent vectors Γ(D) ∈ T D 0 . The mapping D → Log D 0 (D) = Γ(D) is known as the Grassmann logarithm and its inverse mapping as the Grassmann exponential Γ → Exp D 0 (Γ) = D. There also holds that Log D 0 (D 0 ) = 0 and

Algorithm 1 :

 1 Density extrapolation framework g-ext Data: Array desc containing the descriptors for k previous time-steps, p n the descriptor for the current position, C n-1 and S n-1 respectively the molecular orbitals and overlap matrices of the previous time-step, and cref the reference point on the Grassmannian Result: Guess density matrix for time-step n > 1 begin cmat(:, :, n -1) ← Orthonormalization(C n-1 , S n-1 ); gmat(:, :, n -1) ← Log(cref, cmat(:, :, n -1)); desc, p n ← Stabilization(desc, p n

Figure 1 :

 1 Figure 1: Total energy (kcal/mol) as a function of simulation time (fs) for 3hf comparing g-ext(3), g-ext(6) and xlbo with McWeeny purification, using a convergence threshold for the scf algorithm of 10 -5 (left panel) and 10 -7 (right panel). The total energy was shifted of +505 000 kcal/mol for readability.

  n -1, at N t previous times t i , one aims to compute an accurate guess of the density matrix D n at time t n .

Table 1 :

 1 Overview of the system size in terms of number of qm-atoms (N QM ), number of mm-atoms (N M M ) and the total number of (qm) basis functions (N ).

	System N QM N M M ocp 129 4915 1038 N
	appa	31 16 449 309
	dmabn	21	6843 185
	3hf	28 15 018 290

Table 2

 2

				dmabn	appa		3hf
	Method	Average	σ	Average	σ	Average	σ	Average	σ
	xlbo	3.82	0.66	3.98	0.16	3.00	0.03	4.00	0.14
	xlbo/mw	2.95	0.31	3.76	0.56	3.00	0.34	3.96	0.31
	g-ext(3)	2.57	0.84	3.54	0.78	2.95	0.50	3.09	0.41
	g-ext(4)	2.48	0.88	3.14	0.62	2.51	0.50	3.25	0.68
	g-ext(5)	2.25	0.96	3.23	0.75	2.51	0.50	3.30	0.72
	g-ext(6)	2.20	0.96	2.99	0.02	2.51	0.50	3.14	0.56
	From the results in

, we see that the g-ext algorithms systematically outperforms the xlbo method. It is interesting to note that the McWeeny purification step has a sizeable

Table 3 :

 3 Performances of the g-ext(6) method compared with the xlbo algorithm with McWeeny purification. All the results were obtained using a 10 -7 convergence threshold for the root-mean-square increment of the density matrix and are derived from a 1 ps long md simulation, using a 0.5 fs time step. We report the average number of iterations required to converge the scf, together with the associated standard deviation. Note that the first 8 steps were discarded.

		ocp		dmabn	appa		3hf
	Method	Average	σ	Average	σ	Average	σ	Average	σ
	xlbo/mw	5.02	0.17	7.30	0.64	7.49	0.84	7.47	0.63
	g-ext(6)	3.58	0.79	4.23	0.50	4.39	0.57	6.81	0.78

Table 4 :

 4 Short and long-term stability analysis of the g-ext(3) and g-ext(6) methods, compared to the xlbo algorithm with McWeeny purification, for the 3hf system. For each method we report the stf and the ltd and the average number of scf iterations, for three convergence thresholds of the scf algorithm. Conv. 10 -5 Conv. 10 -6 Conv. 10 -7

	Method	stf	ltd	stf	ltd	stf	ltd
	xlbo/mw 0.55 -0.04 0.55 -0.03 0.57 -0.03 g-ext(3) 0.55 -0.42 0.57 -0.15 0.53 -0.04 g-ext(6) 0.56 -0.53 0.52 -0.13 0.57 -0.04

Funding

Part of this work was supported by the French "Investissements d'Avenir" program, project ISITE-BFC (contract ANR-15-IDEX-0003).

Supporting Information Available

A Julia template of the g-ext algorithm is available at https://github.com/epolack/ GExt.jl. The figure representing the total energy computation with an scf convergence threshold of 10 -6 for the molecule 3hf and formulas for the exponential and logarithm functions are available in the supplementary information.

Grassmann Exponential and Logarithm maps

The Grassmann manifold is a differential manifold and, for any given D 0 = C 0 C T 0 ∈ Gr(N, N ) with D 0 := D R 0 and C 0 := C R 0 for fixed R 0 , the tangent space is characterized by

Note that the tangent space is a linear space. One can then introduce the Grassmann exponential which maps tangent vectors on T D 0 to the manifold Gr(N, N ) in a locally bijective manner around D 0 . Indeed, it is not only an abstract tool from differential geometry, but it can be computed in practice involving the matrix exponential. By complementing C 0 with orthonormal columns to obtain (C 0 , C ⊥ ) ∈ O(N ), where O(N ) denotes the group of orthogonal matrices of dimension N × N , and Γ ∈ T D 0 we have

Here, exp denotes the matrix exponential function, the matrix B ∈ R (N -N )×N contains expansion coefficients of columns of Γ in a span of columns of C ⊥ such that Γ = C ⊥ B and I N ,N = (I N , 0) T ∈ R N ×N are the first N columns of the N × N identity matrix. As described in [1, 2], the Grassmann exponential can then be equivalently expressed by

by means of a singular value decomposition (SVD) of the matrix Γ = U e Σ e V T e . The inverse function is the so-called Grassmann logarithm Log D 0 (see, e.g., [1, 2]) which maps any D = CC T ∈ Gr(N, N ) in a neighborhood of D 0 to the tangent space T D 0 by Log D 0 (D) = U arctan(Σ )V T , (4) using the following SVD decomposition