Florian Pescher
email: florian.pescher@femto-st.fr

Nils Napp
email: nnapp@cornell.edu

Benoît Piranda
email: benoit.piranda@femto-st.fr

Julien Bourgeois
email: julien.bourgeois@femto-st.fr

GAPCoD: A Generic Assembly Planner by Constrained Disassembly

Keywords:

In the literature we can find many kinds of modular robot that can build a wide variety of structures. In general, finding an assembly order to reach the final configuration, while respecting the insertion constraints of each kind of modular robot is a difficult process that requires system-specific tuning. In this article, we introduce a generic assembly planner by constrained disassembly (GAPCoD) which works with all kinds of modular robots. It outputs a directed acyclic graph where vertices are modules needing to be placed before his child nodes. This graph is obtained through the disassembly of the desired structure submitted to user chosen constraints. We detail the compiler as well as the way to choose constraints and their influence on performance. The robots embed simple path planning algorithm to reach the destination and act as decentralized agents. Examples are provided to show the possibilities that the compiler offers with two very different robot systems and constraints.

INTRODUCTION

Analyzing assemblies and encoding assembly sequences has been explored since the 1980s, and an algorithm for generating the complete set of partial assemblies and assembly sequences, represented as and-or graphs, exists. However, the scaling of the algorithm depends on several problem-specific factors, but is exponential in most cases [START_REF] Luiz | AND / OR REPRESEN-TATION OF ASSEMBLY[END_REF]. In robotic construction and modular robot reconfiguration, instead of assembling larger sub-assemblies, components are typically added one by one and many difficulties arise from the motion of active robots, which is not represented in and-or graphs. As such, robotic reconfiguration systems mostly develop their own planning algorithms and plan representations. In robotic assembly, plans for reconfiguration are often co-developed with a particular set of hardware, and it is difficult to compare between systems [START_REF] Jenett | Material Robot System for Assembly of Discrete Cellular Structures[END_REF][START_REF] Lindsey | Construction of Cubic Structures with Quadrotor Teams[END_REF][START_REF] Seo | Assembly Planning for Planar Structures of a Brick Wall Pattern with Rectangular Modular Robots[END_REF][START_REF] Tucci | A Distributed Self-Assembly Planning Algorithm for Modular Robots[END_REF]. There are a few abstract models for robotic reconfiguration, but while they are typically based on an original robotic system [START_REF] Butler | Distributed motion planning for modular robots with unit-compressible modules[END_REF][START_REF] Murata | A 3-D self-reconfigurable structure[END_REF], it is often not easy to build robots that fit them. However, they allow for the development of generic reconfiguration algorithms, e.g. [START_REF] Lengiewicz | Efficient collective shape shifting and locomotion of massively-modular robotic structures[END_REF] and the references therein describe recent developments. Instead of solving a motion problem, these strategies build an intermediate (porous lattice) structure to facilitate motion [START_REF] Thalamy | Distributed Self-Reconfiguration using a Deterministic Autonomous Scaffolding Structure[END_REF]. In robotic construction [START_REF] Jenett | Material Robot System for Assembly of Discrete Cellular Structures[END_REF][START_REF] Lindsey | Construction of Cubic Structures with Quadrotor Teams[END_REF][START_REF] Werfel | Designing collective behavior in a termite-inspired robot construction team[END_REF], i.e. active robots manipulating passive components, navigating intermediate structures is paramount, and this is often achieved by structuring the travel paths/directions of robots so that they do not collide. Here, we adapt ideas from planners for the TERMES [START_REF] Werfel | Designing collective behavior in a termite-inspired robot construction team[END_REF] systems to lattice-based robotic systems, such as quasi-spherical 3D Catoms [START_REF] Piranda | Designing a quasi-spherical module for a huge modular robot to create programmable matter[END_REF] in the context of the polymer molding application [START_REF] Pescher | Molding a Shape-Memory Polymer with Programmable Matter[END_REF] and cubic Blinky Blocks [START_REF] Kirby | Blinky blocks: a physical ensemble programming platform[END_REF] with balance constraints. The goal of the TERMES compiler is to make the runtime execution distributed and requires essentially no information exchange between robots at runtime (beyond locally avoiding collision). Assembly checks can be fully determined by the local state of the structure. The planning process is termed "compilation" since it prepares a high-level target for execution by the runtime system (i.e. the robots and raw materials) which do not change. Recent work has significantly increased the speed and scalability of the TERMES compiler [START_REF] Deng | A Compiler for Scalable Construction by the TERMES Robot Collective[END_REF]. The key idea in this compiler is to use a partial order reduction, a technique in formal verification to combat state-space explosion [START_REF] Godefroid | Partial-Order Methods for the Verification of Concurrent Systems: An Approach to the State-Explosion Problem[END_REF]. Some states are not explicitly enumerated, since they can be reached by a combination of independent actions, instead of checking the 2 𝑁 possible combinations of 𝑁 actions, only 𝑁 independent actions need to be checked. This idea is exploited implicitly in many assembly planning algorithms in the literature, e.g. [START_REF] Tolley | Stochastic Modular Robotic Systems: A Study of Fluidic Assembly Strategies[END_REF], by asserting that some actions can happen parallel, i.e. their order does not matter. Instead of focusing on a particular hardware platform or type of motion, connection, or stability constraint, we focus on the way the states are enumerated and order dependencies captured, and show that these ideas can be applied to a variety of systems. The resulting assembly orders are directed acyclic graphs (DAGs), which can be interpreted as a partial order on the assembly steps. This representation can only capture a subset of the full states and the compiler commits to specific order choices when it encounters them. Within these restrictions it attempts greedily find orders that allow for as many parallel steps as possible.

THE GAPCOD PLANNER

We are now going to present GAPCoD, a Generic Assembly Planner by Constrained Disassembly. This planner produces an assembly order schedule according to user-defined constraints by disassembling the desired target shape. First, we show the general workflow of GAPCoD, then we explain how GAPCoD Python compiler part works, detail how to implement different constraints for different physical modules or rules of assembly, and finally we discuss the complexity and scalability of GAPCoD and its dependency on the physical constraints.

GAPCoD workflow

GAPCoD is composed of two main parts: a physical part, whose goal is to check the constraints for each modules, and a logical part in charge of producing the output DAG. The key idea is that the logic for exploring and ordering states can be largely separated from the details of why constraints exist and how they are checked for a particular robot system. In our implementation, the physical part is handled by the VisibleSim [START_REF] Piranda | VisibleSim: Your simulator for Programmable Matter[END_REF] simulator allowing us to simulate constraints for different kinds of modules with ease, and the logical part is handled by a Python compiler which is an adaptation of the TERMES compiler [START_REF] Deng | A Compiler for Scalable Construction by the TERMES Robot Collective[END_REF][START_REF] Werfel | Designing collective behavior in a termite-inspired robot construction team[END_REF]. The general workflow of the compiler as well as the interactions between VisibleSim and the Python compiler can be seen in Figure 1. The detailed physical state is only maintained in the simulator, while the state-traversal logic keeps track of candidate locations in the assembly to check next.

In order to plan an assembly order with GAPCoD, we must first determine and implement the constraints that will be applied in our scenario. More details about how to implement constraints can be found in Section 2.3. That part is the specific tuning one must perform in order to use GAPCoD with different kinds of systems. Once the constraints are implemented, we can start the planning for our structures by giving GAPCoD a shape description and seed modules, i.e. the last and first (free) modules to be disassembled. In our current implementation, we limited ourselves to a single initially free module.

Python compiler explanation

The compiler computes a DAG in which the roots are the seed modules and every vertex represent a module that can only be assembled if all of its parents modules have already been assembled. This graph is created and updated by disassembling the complete structure that we are trying to find an assembly order for. All sets and variables used by the compiler are described in Table 1 Algorithm 1 shows the global process of the graph creation until the structure has been fully disassembled. In this section, whenever we refer to a removable or free module, we are talking about a module that can be disassembled without breaking any of the constraints the compiler has been given. Let isRemovable(𝑘,𝐶) be the function that returns true if the module 𝑘 can be removed from 𝐶 without breaking any constraints (and false otherwise).

The compiler is based on two main functions called all along the disassembly of the structure: Disassemble(𝑛,𝐶,𝑅) that disassemble a module 𝑛, updating both the list of removable modules 𝑅 and the structure 𝐶 (detailed in Algorithm 2 ), and Reassemble(𝑛,𝐶,𝑅) that reassemble a disassembled module 𝑛, updating both the list of removable modules 𝑅 and the structure 𝐶 (detailed in Algorithm 3). The disassembly of the currently free modules can be divided into three main steps: Checking blockages (Algorithm 4), finding and updating the dependencies (Algorithm 5) and creating required additional edges (Algorithm 6). One execution of each of those three steps will be called a "level" in the rest of the paper. The algorithms are detailed in the following paragraphs.

VisibleSim Python compiler

Description

2.2.1

Step 1: Checking blockages. Before actually removing a module from the structure, we must ensure that the removal of this module will not prevent the removal of another module trying to disassemble itself this level. The reasons for those blockages depend on the particular constraints. Some examples would be the removal of a module used as a pivot for the movement of another one (as it can be seen in Figure 2), or making another module the only module responsible for maintaining the connectivity in the structure. Detecting and preventing those blockages according to the constraints given to the compiler is the goal of Algorithm 4.

When such a blockage is detected the blocker module is added to a list called BL. The modules in BL will be considered as not free for the rest of the current level. Although considering the blocker modules as not free until the next level solve most problems, we must still be careful about the case of two modules blocking each other, causing a deadlock. When a deadlock happens, the first module tested for blockage is considered free again, thus giving it priority for disassembly. An example showing how this interblocking situation is handled is given in Figure 6. Function CheckBlockages(𝐶,𝑅):

𝐵𝐿 ←-∅ ; ∀𝑖 ∈ 𝐶; 𝐵𝐿 𝑖 ←-∅ ; forall 𝑛 ∈ 𝑅 do {𝐶, 𝑅 ′ } ←-Disassemble(𝑛,𝐶,𝑅); if (∃𝑚 ≠ 𝑛; 𝑚 ∈ 𝑅 ∧ 𝑚 ∉ 𝑅 ′ ) then 𝐵𝐿 ←-𝐵𝐿 ∪ 𝑛; 𝐵𝐿 𝑛 ←-𝐵𝐿 𝑛 ∪ 𝑚; if (𝐵𝐿 𝑚 = {𝑛}) then 𝐵𝐿 ←-𝐵𝐿 -𝑚 ; {𝐶, 𝑅} ←-Reassemble(𝑛,𝐶,𝑅 ′ ); 𝑅 ←-𝑅 -𝐵𝐿 ; return {𝑅, 𝐵𝐿} ; 2.2.2
Step 2: Finding and updating the dependencies. Considering that the list of removable modules has been updated to prevent blockages, we can now start actually disassembling this level of the structure and updating the output graph (cf. Algorithm 5). During this step, we disassemble all the free modules of this level and whenever a new module is freed we create the corresponding edge in the output graph. Some modules might need the removal of several other modules in order to be freed i.e. when we find a newly freed module, we check if the previously removed modules during this level impact the removability of this freed module and create the corresponding dependencies. Figure 4 shows an example of how this combo check is performed on a simple configuration.

This algorithm flags every free module either as "with dependencies" or "without dependencies" with the exception of seed modules that will not free any other module. An example of a module with no dependencies is given in Figure 5.

2.2.3

Step 3: Creating required additional edges. Algorithm 6 is the last step of the disassembly process during which it constructs the graph 𝐺.

𝑅 ′ ←-𝑅 ; forall (𝑛 ∈ 𝑅) do {𝐶, 𝑅 ′ } ←-Disassemble(𝑛,𝐶,𝑅 ′ ); if (∃𝑚; 𝑚 ∈ 𝑅 ′ ∧ 𝑚 ∉ 𝑅) then 𝐸 ←-𝐸 ∪ (𝑚 → 𝑛); 𝐷𝑒𝑝 𝑛 ←-𝑇𝑟𝑢𝑒 ; forall (𝑛 ′ ∈ 𝐷) do {𝐶, 𝑅 ′ } ←-Reassemble(𝑛 ′ ,𝐶,𝑅 ′ ) ; if (𝑚 ∉ 𝑅 ′ ) then 𝐸 ←-𝐸 ∪ (𝑚 → 𝑛 ′ ); 𝐷𝑒𝑝 𝑛 ′ ←-𝑇𝑟𝑢𝑒 ; {𝐶, 𝑅 ′ } ←-Disassemble(𝑛 ′ ,𝐶,𝑅 ′ ); else if (𝑛 ∉ 𝑆) then 𝐷𝑒𝑝 𝑛 ←-𝐹𝑎𝑙𝑠𝑒 ; 𝐷 ←-𝐷 ∪ 𝑛; 𝑅 ←-𝑅 ′ ; return {𝐶, 𝑅 ′ , 𝐷, 𝐸, 𝐷𝑒𝑝} ;
During this step, we create additional edges in 𝐺 between the modules that will still be present at the next level and the ones that were just removed (flagged as "with dependencies"). Those additional edges ensure that when we will try to assemble the structure, all the modules that must be present for this assembly actually will be. Let us show that those edges are necessary to obtain a correct order. For this example we suppose that we have square modules which can only be disassembled/assembled if there is a module or a wall directly connected on their left side and no module above them. We consider the configuration shown in Figure 3 and perform the compiler process on this configuration with (case 𝐴 in grey) and without (case 𝐵) the call to the CreateEdges function described in Algorithm 6.

Considering that all seed modules are placed, if we try to assemble the structure following the graphs, we get 6 order possibilities for case 𝐵, but only one is correct. Case 𝐴 is more efficient as it gives the only correct solution.

{R,BL} ← CheckBlockages(C,R) C={0,1,2,3,4,5}, R={5}, BL={3,4} BL_3={4} BL_4={5} {C,R,E,Dep} ← UpdateDependencies(R,C,E,S) C={0,1,2,3,4}, R={2} E=E U {2->5} Dep_5=True {C,R,E} ← CreateEdges(C,R,E,BL,Dep) C={0,1,2,3,4}, R={2} E ← E U {3->5,4->5} R ← R U BL  R={2,3,4}
Level 2:

{R,BL} ← CheckBlockages(C,R) C={0,1,2,3,4}, R={2,3,4}, BL={3} BL_3={4} {C,R,E,Dep} ← UpdateDependencies(R,C,E,S) C={0,1,3}, R={1} E ← E U {1→4} Dep_4=True (2 frees nothing but is seed) {C,R,E} ← CreateEdges(C,R,E,BL,Dep) C={0,1,3}, R={1} E ← E U {3→4} R ← R U BL  R={1,3} {R,BL} ← CheckBlockages(C,R) C={0,1,3}, R={1,3}, BL={} {C,R,E,Dep} ← UpdateDependencies(R,C,E,S) C={0}, R={0} E ← E U {0->3} Dep_3 ← True (1 frees nothing but is seed) {C,R,E} ← CreateEdges(C,R,E,BL,Dep) C={0}, R={0} E ← E U {}

Implementing new constraints

This section details which kind of constraints may be used with our compiler. A constraint is compatible with GAPCoD if it is associated to a function that returns true if the module can be removed according to the constraint and false otherwise.

All constraints have a few pieces of information that will have an influence on the compiler efficiency:

• The number of influence neighbors. Influence neighbors of a module are the modules that can be affected (meaning modules that can change their removability state) when assembling/disassembling this module according to the considered constraint. • The complexity of the boolean function checking if the constraint (mobility constraint, balance constraint. . . ) is satisfied.

We are now going to detail some examples of constraints: the connectivity constraint as most modular robotic systems needs their modules to be connected to one another at all times, the 3D Catoms mobility constraints and the balance constraint of a set of Blinky Blocks.

Connectivity.

When removing or adding a module according to the connectivity constraint only the physical neighbors are affected. This leads to 6 influence neighbors for modules in a cubic lattice such as the Blinky Blocks or 12 influence neighbors for the 3D Catom that are organized in Face-Centered Cubic (FCC) lattice. Despite having only a few influence neighbors the connectivity constraint is still a global constraint. In order to check if a module 𝑚 can be removed we need to verify that all modules can be reached from any module 𝑠 ∈ 𝑆 even after the removal of 𝑚. This lead to a complexity of 𝑂 (𝑛) where 𝑛 is the number of modules in 𝐶.

3D Catom mobility constraints. The movement of the 3D

Catoms are subject to a lot of constraints. Indeed to perform a single move we must make sure that the required pivot is present and that a bunch of other positions are free. Theoretical motions of 3D Catoms [START_REF] Piranda | Designing a quasi-spherical module for a huge modular robot to create programmable matter[END_REF] give 120 different moves for a 3D Catom. However in the context of the polymer molding application presented in [START_REF] Pescher | Molding a Shape-Memory Polymer with Programmable Matter[END_REF] we can reduce this number to 84.

BL ← {} {C,R'} ← Disassemble(1,C,R) C={0,2,3,4,5},R'={},R={1,4} 4 (R -R')  blockages BL ← BL U 1  BL={1}; BL 1 ={} BL 1 ← BL 1 U 4  BL={1}; BL 1 ={4} {C,R} ← Reassemble(1,C,R') C={0,1,2,3,4,5}, R={1,4} {C,R'} ← Disassemble(4,C,R) C={0,1,2,3,4,5}, R'={}, R={1,4} 1 (R -R')  blockages BL ← BL U 4  BL={1,4}; BL 1 ={4}, BL 4 ={1} BL 4 ← BL 4 U 1  BL={1,4}; BL 1 ={4}; BL 4 ={1} BL 1 ={4}  BL ← BL -{1} {C,R} ← Reassemble(4,C,R') C={0,1,2,3,4,5}, R={1,4} R ← R-BL;  R={1} {R,BL} ← {R,BL} R={1},BL={4},BL 4 ={1}
Those moves can all be checked by knowing if 48 specific neighboring cells of the lattice are filled or empty. We can notice that the 12 actual physical neighbors of the 3D Catom are included in those 48 influence neighbors.

The occupancy requirements (limited to the 12 physical neighbors minus the one above) for all possibles move can be summarized in a 84 × 11 matrix which when multiplied by an occupancy state of those 11 cells where the value is 0 if the cell is empty and 1 if the cell is filled gives us a vector 𝑉 1 of all possible moves. 𝑉 1 is a vector of presence of the neighbors (component equal to 1 if the pivot is present and 0 otherwise).

In the same way we can verify the empty cells requirements by multiplying a 84 × 48 matrix with the opposite of the occupancy state of all 48 influence neighbors (vector 𝑉 2 ). In that way we designed a way to check with a constant complexity if a 3D Catom is able to move.

Blinky Blocks balance of a T-shapes.

The last constraint is the stability of T-shapes like structures limited in the (𝑥, 𝑧) plane for the Blinky Blocks. In this case, the influence neighbors for this global constraint are all the 𝑛 modules of 𝐶.

In order to check if the structure is stable we first compute the gravity center 𝐺 of the shape. We consider that the structure is stable if the vertical line passing by 𝐺 cross the basis of the structure as shown in Figure 7.

The complexity of this check is then the complexity to compute the gravity center of the shape and as such can be done in 𝑂 (𝑛).

Complexity

We detailed how the compiler works in Section 2.2 and how constraints for the compiler can be implemented in Section 2.3, allowing us to express the complexity of the compiler depending on the number of modules in the final assembly, 𝑛. Let 𝑘 be the number of influence neighbor and define the highest complexity to check a constraint as 𝑂 (𝑀).

In the worst case scenario all the 𝑛 modules are free at a given level but only one is definitely removed during this level, leading to 𝑛 levels. However, such a system is both difficult to imagine, and would not be a good candidate for this planner, which exploits the fact that parallel actions reduce the number of states that need to enumerated, and that any any given time only a limited number of actions are potentially parallel.

By definition we have 𝑖𝑠𝑅𝑒𝑚𝑜𝑣𝑎𝑏𝑙𝑒 (𝑘, 𝐶) with a complexity of 𝑂 (𝑀). Knowing that we can now compute that Algorithm 2 and 3 are in 𝑂 (𝑘𝑀). This allows us to find the complexity for Algorithms 4,5 and 6: Algorithm 4 is in 𝑂 (𝑛𝑘𝑀), Algorithm 5 in 𝑂 (𝑛 2 𝑘𝑀) and Algorithm 6 in 𝑂 (𝑛 2 ) or 𝑂 (𝑛𝑘𝑀) whichever is the highest. Finally this leads to a worst case complexity of 𝑂 (𝑛 3 𝑘𝑀) for the complete compiler described in Algorithm 1.

This complexity seems important but let us not forget than in most cases not all modules will be free during a level and that the compiler will usually not perform as many levels as modules in the structure. Plus the main bottleneck in our compiler is the combo check performed in Algorithm 5, however it is this combo check that yields one of the most important feature of our compiler: the ability to find possible parallel operations that will result in a time gain when assembling the structure. In other words we choose to increase compile time in order to reduce assembly time.

EXPERIMENTATION

Finally we are going to use this compiler to find an assembly order in different scenarios with different kind of modules associated to different set of constraints. We study the scalability of those different scenarios and we talk about the complexity of the compiler (using the notations of Section 2.4).

3D Catom assembly

The first scenario we will study is an assembly order with the 3D Catoms module associated to a mobility constraint (can the 3D Catom be inserted at a given position of the grid) and a connectivity constraint (the structure must always be connected).

Figure 9 shows some steps of the assembly of a wall shape made of 220 3D Catoms. The linked video 1 presents the full construction of the wall.

Let us now compute the complexity of the compiler in this particular scenario. The number of influence neighbors is 𝑘 = 48: We have 12 physical neighbors for the connectivity constraint and 48 mobility influence neighbors, including the 12 physical ones. Then the constraint that is checked with the highest complexity is the connectivity which is done in 𝑂 (𝑛) against the mobility constraint which is done in 𝑂 (1), giving us 𝑀 = 𝑂 (𝑛). Given that the complexity of the compiler in any scenario is 𝑂 (𝑛 3 𝑘𝑀), in this case we have a worst case complexity of 𝑂 (𝑛 4 ).

Let us now compute an approximate experimental complexity to compare with this worst case theoretical one. In order to do so we are going to run the compiler on different size of cubes and check both the time it takes the compiler to produce an order as well as the dept of the produced graph to see the efficiency of our method. The results can be found in Table 2.

If we take the disassembly time for the six smaller cubes and compute the polynomial tendencies as a function of the number of module and compare it to the actual disassembly times we get the graph in Figure 10.

Let us now analyze those results. We can see that the experimental disassembly time is best fitted by 𝑛 2.37 which is much better than the theoretical 𝑂 (𝑛 4 ). This means that the compiler performs better than the worst case scenario in this case. This difference is explained by the fact that in the worst case we supposed that all modules were free every level and that we had as many levels (the graph depth) as modules in the structure when in practice that is not the case as shown by the graph depth column in Table 2. The other interesting to notice is that the depth of the output graph does not increase drastically as the number of modules increase, meaning that for bigger structure the compiler finds more possible parallel operations which implies possible time gains when we will try to assemble the structure. In other words we are performing a time expensive computation to produce the assembly graph for a given shape in order to gain time when assembling the structure by maximizing the number of parallel operations. 2. The best fit scaling exponent is 𝑂 (𝑛 2.37 ), which is significantly better than the worst case analysis 𝑂 (𝑛 4 ).

Blinky Blocks T-shape example

We will now look into a T-shape assembly order for the Blinky Blocks submitted to the mechanical stability constraint described in Section 2.3.3 as well as the connectivity constraint. We can see the graph generated by the compiler as well as all the possible assembly orders we can find by following this graph for a small example in Figure 11. Although it should be possible to assemble the last module in any branch first, the way the compiler is currently implemented, the structure is disassembled starting by an initial free module chosen by the user. In this example we chose the last module of the right branch as the only initially free module which is the reason why we end up assembling the last module of the left branch before the last module of the right branch.

Let us now compute the complexity of the compiler in this scenario. The number of influence neighbors is 𝑘 = 𝑛: We have the 6 physical neighbors for the connectivity constraint, and 𝑛 for the stability constraint which affects the whole structure. Then both constraints are checked with a complexity of 𝑂 (𝑛) leading us to 𝑀 = 𝑂 (𝑛). Given that the complexity of the compiler in any scenario is 𝑂 (𝑛 3 𝑘𝑀), in this case we have a complexity of 𝑂 (𝑛 5 ).

Discussion

The experiments show that the compiler can produce assembly orders for different kinds of constraints. However, even though the compiler is adaptable there is a cost in scalability, especially for constraints that affect the whole structure and/or need a global check at each step. We decided to prioritize the correctness of the order and the number of possible parallel operation of the assembly plan over the runtime of the planner. As a result the scaling of the planner depends on the hardest constraint check times the number of influence neighbors. In situations where the influence of a module is limited, the worst-case complexity analysis is very conservative and GAPCoD is ≈ 1.5 orders lower than the worst case polynomial power. Examples for the 3D Catoms and the Blinky Blocks illustrate that global constraints incur an additional 𝑛 scaling factor suggesting some loss of efficient scalability. However global constraints can take into account constraints that can not be checked locally, such as the mechanical constraint of stability for T-shapes. This means that when giving constraints to the compiler a tradeoff between adaptability and scalability must be done. The current implementation does not optimize the removability check in the simulator. However, the way GAPCoD traverses assembly states consecutive checks are usually for similar states, which means they are easy targets for optimization by caching some of the repetitive computation, similar to using a spanning tree to speed up connectivity checking in [START_REF] Deng | A Compiler for Scalable Construction by the TERMES Robot Collective[END_REF].

Another thing to notice in the experiments is that we only used the graph to find a correct assembly order and then sequentially assemble modules. Even though it shows the compiler gives a correct assembly order, the actual assembly does not make use of all the interesting properties of the output graph. Since the graph indicates which modules should be placed before assembling a specific module this allows for some operations to be parallelized once the requirements are met. If we take the graph from Figure 11a as an example, after assembling (0, 5) we can assemble either (-1, 5), [START_REF] Butler | Distributed motion planning for modular robots with unit-compressible modules[END_REF][START_REF] Kirby | Blinky blocks: a physical ensemble programming platform[END_REF] or both at the same time in parallel. By maximizing the number of parallel operations we can reduce the time it takes for the assembly to be performed. If we call a time step the maximum time it takes a single module to be placed for a given structure, by maximizing the number of parallel operations we reduce the number of time steps to the graph depth thus reducing the assembly time. In the example from Figure 11 by maximizing the number of parallel operations we can go from 12 to 10 time steps. This speedup is even more impressive in the example of cubes made of 3D Catoms as can be seen in Table 2 where for a cube made of 1400 modules we go from 1400 time steps to 55 if we perform all possible parallel operations at the same time. This theoretical speedup by maximizing the number of parallel operation performed at the same time seems very promising and will be the object of future works.

CONCLUSION

In this paper we introduced GAPCoD a Generic Assembly Planner by Constrained Disassembly. We explained how this planner produces a DAG representing possible orders of assembly with as much parallel operations as possible. We then tested this method in two different scenarios showing the adaptability of our method to different constraints and physical modules. Future works will focus on choosing between different possible DAG orders to optimize assembly time as well as testing GAPCoD on many other scenarios than the one presented here.

Algorithm 4 :

 4 CheckBlockages(C,R): check blockages and remove interlocks

Figure 2 :

 2 Figure 2: 3D Catom in rotation around a pivot.

Algorithm 6 :

 6 Pushes blocker and no dependencies modules to next level + creates required edges in graph Function CreateEdges(𝐶,𝑅,𝐷,𝐸,𝐵𝐿,𝐷𝑒𝑝): forall (𝑛 ∈ 𝐵𝐿) do forall (𝑚 ∈ 𝐷; 𝐷𝑒𝑝 𝑚 = 𝑇𝑟𝑢𝑒) do 𝐸 ←-𝐸 ∪ (𝑛 → 𝑚); forall (𝑛 ∈ 𝐷; 𝐷𝑒𝑝 𝑛 = 𝐹𝑎𝑙𝑠𝑒) do {𝐶, 𝑅} ←-Reassemble(𝑛,𝐶,𝑅); forall (𝑚; 𝐷𝑒𝑝 𝑚 = 𝑇𝑟𝑢𝑒) do 𝐸 ←-𝐸 ∪ (𝑛 → 𝑚); return {𝐶, 𝑅, 𝐸} ;

Figure 3 :

 3 Figure 3: Creating graphs with (Case A in grey) and without Algorithm 6. The middle column shows the operations performed by the algorithm and their influence of the different sets. Left column is a visual representation of those sets showing the configuration at each step of the compiler and right column shows the graph created at each of those steps.

Figure 4 :

 4 Figure 4: Example 1 illustrating the algorithm. In this example a module can be removed if and only if the structure connectivity is maintained. (Same color scheme as Figure 3.)

Figure 5 :

 5 Figure 5: Example 2 illustrating the algorithm. In this example a module can be removed if and only if the structure connectivity is maintained. (Same color scheme as Figure 3.)

1 4 Stage 2 : 4 BL={4} and Dep 1 =Figure 6 :

 142416 Figure 6: Example 3 illustrating the algorithm. In this example a module can be removed if and only if the structure connectivity is maintained. (Same color scheme as Figure 3.)

Figure 7 :

 7 Figure 7: Example for a stable and an unstable configuration for the T-shape stability test

Figure 8 :

 8 Figure 8: Complexity computation for the algorithms

Figure 9 :

 9 Figure 9: Assembly of a wall of 3D Catoms.

Figure 10 :

 10 Figure 10: Experimental disassembly times starting from a varying number of modules, see Table2. The best fit scaling exponent is 𝑂 (𝑛 2.37 ), which is significantly better than the worst case analysis 𝑂 (𝑛 4 ).

  Graph output for the assembly of the T-shape with Blinky Blocks. Basis of the T has coordinate (0, 0) in the (𝑥, 𝑧) plan. This graph does not represent all possibilities, only the ones found by the compiler according to the constraints and initial free modules are drawn.(b) Different possible orders according to the graph.

Figure 11 :

 11 Figure 11: Compiler output and hardware assembly with Blinky Blocks for a T-shape.

Table 1 :

 1 . Table of sets and variables used by the compiler

	Symbol	Description
	𝐶	Set of connected modules
	𝐺 (𝑉 , 𝐸)	Assembly graph where 𝑉 is the set of vertices and 𝐸 is the set of edges
	𝑆	Set of seed modules
	𝑅, 𝑅 ′	Sets of removable modules
	𝐵𝐿	Set of modules blocking the removal of another module when removes
	𝐵𝐿 𝑖	Set of modules blocked by the removal of module 𝑖 ∈ 𝐵𝐿
	𝐷	Set of disassembled modules
		Boolean: True if there is an edge pointing towards
	𝐷𝑒𝑝 𝑖	the node representing module 𝑖 ∈ 𝐷
		in the graph False otherwise
	𝑛, 𝑛 ′ , 𝑚, 𝑘 Lower case single characters represents a module

1

  YouTube video of the assembly of the wall: https://youtu.be/s_-LyYyVnCU

	Nb of modules Graph Depth Disassembly time (s)
	1	1	0.002279
	12	6	0.366118
	36	11	5.386122
	96	19	67.193171
	175	24	263.753336
	288	29	990.736958
	490	37	3724.907382
	704	41	9992.689347
	1053	48	23334.396296
	1400	55	47485.460037

Table 2 :

 2 Graph depth and disassembly time for different size cubes made of 3D Catoms.

ACKNOWLEDGMENTS

This work was partially supported by the ANR (ANR-16-CE33-0022-02), the French Investissements d'Avenir program, ISITE-BFC project (ANR-15-IDEX-03), Labex ACTION program (ANR-11-LABX-01-01),the Mobilitech project, and the National Science Foundation (NSF #1846340).


Louis Andreoli

Xavier Porte
email: javier.porte@femto-st.fr.

Stéphane Chrétien

Maxime Jacquot

Daniel Brunner

Laurent Larger

Keywords: Boolean learning, neural networks, noise

A high efficiency hardware integration of neural networks benefits from realizing nonlinearity, network connectivity and learning fully in a physical substrate. Multiple systems have recently implemented some or all of these operations, yet the focus was placed on addressing technological challenges. Fundamental questions regarding learning in hardware neural networks remain largely unexplored. Noise in particular is unavoidable in such architectures, and here we experimentally and theoretically investigate its interaction with a learning algorithm using an opto-electronic recurrent neural network. We find that noise strongly modifies the system's path during convergence, and surprisingly fully decorrelates the final readout weight matrices. This highlights the importance of understanding architecture, noise and learning algorithm as interacting players, and therefore identifies the need for mathematical tools for noisy, analogue system optimization.

Introduction

In recent years, neural networks (NNs) take centre-stage in advancing computation [START_REF] Lecun | Deep learning[END_REF]. Optimized by training, such learning machines provide key advantages for solving abstract computational problems and already outperform humans in numerous tasks previously deemed impossible for classically (algorithmically) programmed computers [START_REF] Lecun | Deep learning[END_REF][START_REF] Amos | OpenFace: A General-Purpose Face Recognition Library with Mobile Applications[END_REF][START_REF] Graves | Speech recognition with deep recurrent neural networks[END_REF].

However, NNs are still mostly emulated by traditional Turing/von Neumann computers. The absence of computing hardware supporting fully parallel NNs reduces energy efficiency and overall speed, and new paradigms addressing these problems are desirable. The lack of a parallel network substrate is a fundamental roadblock and is an active area of research since decades, with current analogue hardware either implementing the full network [START_REF] Lin | All-optical machine learning using diffractive deep neural networks[END_REF][START_REF] Shen | Deep learning with coherent nanophotonic circuits[END_REF][START_REF] Tait | Broadcast and weight : an integrated network for scalable photonic spike processing[END_REF] or the neurons [START_REF] Appeltant | Information processing using a single dynamical node as complex system[END_REF][START_REF] Duport | All-optical reservoir computing[END_REF][START_REF] Larger | Photonic information processing beyond Turing: an optoelectronic implementation of reservoir computing[END_REF][START_REF] Brunner | Reconfigurable semiconductor laser networks based on diffractive coupling[END_REF][START_REF] Torrejon | Neuromorphic computing with nanoscale spintronic oscillators[END_REF]. An implementation of nonlinear neurons, fully-parallel information transduction and learning on a substrate level promises a revolution, and photonic NNs [START_REF] Farhat | Optical implementation of the Hopfield model[END_REF][START_REF] Psaltis | Optical information processing based on an associative-memory model of neural nets with thresholding and feedback[END_REF] remain a highly promising avenue [START_REF] Lin | All-optical machine learning using diffractive deep neural networks[END_REF][START_REF] Shen | Deep learning with coherent nanophotonic circuits[END_REF].

Noise is an inseparable companion of analogue hardware [START_REF] Semenova | Fundamental aspects of noise in analog-hardware neural networks[END_REF], yet the fundamental aspects of optimizing a noisy NN [START_REF] Hermans | Embodiment of learning in electro-optical signal processors[END_REF][START_REF] Antonik | Brain-inspired photonic signal processor for generating periodic patterns and emulating chaotic systems[END_REF][START_REF] Bueno | Reinforcement learning in a large scale photonic recurrent neural network[END_REF] have so far hardly been exploredneither in experiments [START_REF] Alata | Phase noise robustness of a coherent spatially parallel optical reservoir[END_REF][START_REF] Soriano | Optoelectronic reservoir computing: tackling noise-induced performance degradation[END_REF] nor in theory [START_REF] Semenova | Fundamental aspects of noise in analog-hardware neural networks[END_REF]. Here, we investigate the interactions between noise, learning rules and the topology of an error landscape for the first time. We experimentally implement a NN with 961 electro-optical neurons via a spatial-light modulator (SLM) [START_REF] Bueno | Reinforcement learning in a large scale photonic recurrent neural network[END_REF], use diffraction [START_REF] Lin | All-optical machine learning using diffractive deep neural networks[END_REF][START_REF] Shen | Deep learning with coherent nanophotonic circuits[END_REF][START_REF] Brunner | Reconfigurable semiconductor laser networks based on diffractive coupling[END_REF][START_REF] Maktoobi | Diffractive coupling for photonic networks: how big can we go?[END_REF] to physically realize the network's internal connections and a digital micro-mirror device (DMD) for programmable Boolean readout weights [START_REF] Bueno | Reinforcement learning in a large scale photonic recurrent neural network[END_REF]. The particular NN task consists in one-step-ahead prediction of the chaotic Mackey-Glass times series. Learning exclusively modifies the readout connections [START_REF] Jaeger | Harnessing nonlinearity: predicting chaotic systems and saving energy in wireless communication[END_REF] via an evolutionary Boolean algorithm based on the error gradient only, and optimization is using either fully random (Markovian) or structured (greedy) exploration.

The statistics of the experimentally obtained learning trajectories prove that noise and exploration strategy strongly interact. Noise induces a kind of random forcing upon the descent algorithm, which strongly modifies the system's path during its convergence towards a local minimum. We find that noise decorrelates the final weight configurations: starting from identical weight configurations and exploring the error landscape's dimensions in identical sequences always leads to clearly differentiated local minima. Quite astonishingly, all minima are spaced at an almost constant distance from each other, which for the generally non-trivial error landscape topologies is unusual at the least. Noise therefore appears to arrange minimizers in periodic positions, much like competitive Brownian walkers with non-local interactions [START_REF] Heinsalu | Competitive brownian and lévy walkers[END_REF]. These fundamental effects highlight the importance of considering hardware architecture, noise and learning algorithm as intimately linked.

Neural network hardware

A recurrent NN inspired by reservoir computing, illustrated in Figure 1a, was our experimentally realized NN test bench. Figure 1b schematically depicts the experiment. An optical plane wave E 0 illuminates the SLM's pixels, and the reflected field is filtered by a polarizing beam splitter (PBS). The SLM combined with the PBS creates a cos(•) non-linearity and the SLM's pixels physically encode the NNs state. A quarter wave plate located between the PBS and the mirror directs the signal towards a camera, and a double pass through the diffractive optical element (DOE) establishes the recurrent connections W DOE [START_REF] Brunner | Reconfigurable semiconductor laser networks based on diffractive coupling[END_REF][START_REF] Bueno | Reinforcement learning in a large scale photonic recurrent neural network[END_REF][START_REF] Maktoobi | Diffractive coupling for photonic networks: how big can we go?[END_REF]. Camera state

x cam i (n) at integer time n x cam i (n) α ∑ N j W DOE i, j E j (n) 2
is combined with external input information u (n+1) and sent to the SLM, creating the network's state according to

x i (n + 1) α|E 0 i | 2 cos 2 β ⋅ x cam i (n) + γW inj i u(n + 1) + θ i . (1)
Here, N = 961 is the recurrent layer's number of nodes, β = 0.8 the feedback gain, γ = 0.4 the input injection gain and α a normalization parameter. The optical electric field and the nonlinearity's bias offset for node i are E 0 i and θ i , respectively. Input information u(n+1) is injected into the system according to random connections W inj . The polarization reflected by the PBS is imaged onto the DMD, whose mirrors are programmed to fixed angles of ±12°from normal incidence. A photodiode only detects optical signals reflected of mirrors with -12°, thus implementing a Boolean readout weight matrix W DMD i 1…N (k). The RC's output is

y out (k, n + 1)∝ ∑ N i W DMD i (k) E 0 i -E i (n + 1) 2 ∝ ∑ N i W DMD i (k)x i (n + 1) 2 . (2)
Here, k is the learning epoch and xi is the optical field of node i arriving at the detector. As in reservoir computing, we restrict learning to the optimization of the readout weights. Finally, the absence of negative weights is partially mitigated by distributing the offset phases θ i | i 1⋯N randomly between θ 0 +δθ i and θ 0 +Δθ+δθ i , where δθ i is a random Gaussian distribution [START_REF] Bueno | Reinforcement learning in a large scale photonic recurrent neural network[END_REF]. Internal W DOE and readout W DMD (k) connections are, however, realized in passive and fully parallel photonic hardware. As the network is constructed of physical neurons it harbours noise, which can either be additive or multiplicative, as well as correlated or uncorrelated [START_REF] Semenova | Fundamental aspects of noise in analog-hardware neural networks[END_REF]. The main sources of noise in our experiment are the SLM and the camera, in relation to which the illumination laser and output detector can be considered as noiseless, and so are the internal coupling and readout matrices implemented by the DOE and DMD, respectively. All relevant noise sources are therefore reservoir-internal, and our following discussion is by no means limited to systems where the readout layer is implement physically. More details about the theoretical treatment and propagation of noise in NNs as well as the individual noise sources and their respective amplitudes and statistics can be found in [START_REF] Semenova | Fundamental aspects of noise in analog-hardware neural networks[END_REF].

Boolean evolutionary learning

Most current learning techniques require complete knowledge of the internal network's state [START_REF] Jaeger | Harnessing nonlinearity: predicting chaotic systems and saving energy in wireless communication[END_REF], all connection weights and potentially all gradients [START_REF] Lecun | Deep learning[END_REF]. In a hardware network this demands probing (and most probably externally storing) the value of each node and connection, which necessitates auxiliary circuitry of a complexity potentially exceeding the actual neural network. This jeopardizes precisely the benefits one targets when mapping a neural network onto hardware. We therefore employ learning that only tracks the computation error's evolution, and hence imposes no constraint on the type of neurons, and more broadly, on hidden layers as a whole. Such an implementation's complexity does not depend on, and hence does not limit the NN's size, which is crucial considering the importance of scalability for computing.

Here, we optimize the DMD's configuration simply by measuring the impact of output mirrors' modifications onto computing error ϵ(k). The objective is to modify W DMD (k) during k = 1, 2, … K learning epochs such that output y out (n + 1) best approximates target τ (n + 1). Our Boolean learning algorithm can be divided into three conceptual sections:

Mutation

W select (k) rand(N) ⋅ W bias (k), (3)
l(k) max W select (k) , (4)
W DMD l(k) (k + 1) ¬ W DMD l(k) (k) , (5)
W bias (k + 1) 1 N + W bias (k), W bias l(k) 0. ( 6
)
We create a vector with N random elements, independently and identically distributed between 0 and 1 (rand (N )). W bias offers the possibility to modifying the otherwise stochastic W select (k) ϵ ℝ N in Eq. ( 3), whose largest entry's position, l(k), determines the Boolean readout weight W DMD l(k) (k) to be mutated via a logical inversion operator (¬(⋅)), see Eqs. ( 4) and [START_REF] Shen | Deep learning with coherent nanophotonic circuits[END_REF].

A fully stochastic Markovian descent is obtained with W bias = 1 and excluding Eq. ( 6). However, here we also investigate exploration which makes mutating a particular connection in near succession unlikely. There, W bias is randomly initialized at k = 1, and at each epoch Eq. [START_REF] Tait | Broadcast and weight : an integrated network for scalable photonic spike processing[END_REF] increases the bias of all connections by 1/N, while the currently modified connection's bias is set to zero. On average, the probability of again probing a particular weight reaches unity only after N learning epochs have passed, and we therefore refer to this biased exploration as greedy learning.

According to these instructions, our algorithm only probes, hence potentially mutates one mirror at a time. We have considered updating more than one mirror at each k, yet simplified numerical simulations indicate that convergence was significantly faster for updating only one weight at a time.

Error and reward signals

ϵ(k) 1 T ∑ n 1 T τ(n + 1) -ỹout (k, n + 1) 2 , (7)
r(k) 1 if Δϵ(k) < 0 0 if Δϵ(k) ≥ 0 (8) ϵ min (k) (1 -r(k))ϵ(k -1) + r(k)ϵ(k), (9)
k min (1 -r(k))k min + r(k)k. (10)
Mean square error ϵ(k) is obtained from a sequence of T data points according to Eq. ( 7), and comparison to the previous error assigns a reward r(k) = 1 only if a modification Δϵ(k) ϵ(k) -ϵ(k -1) was beneficial, see Eq. [START_REF] Duport | All-optical reservoir computing[END_REF]. In that case, the minimum error ϵ min k and the best learning epoch k min are updated following Eqs. ( 9) and [START_REF] Brunner | Reconfigurable semiconductor laser networks based on diffractive coupling[END_REF]. From the system's output y out (k, n + 1) we subtract the mean and normalize by its standard deviation, creating output ỹout (k, n + 1) which is used in Eq. ( 7).

Descent action

W DMD l(k), k r(k)W DMD l(k) (k) + (1 -r(k))W DMD l(k) (k -1) . ( 11
)
Based on reward r(k), the DMD's current configuration either accepts or rejects the previous modification, Eq. ( 11). For a noise-less system, reward r(k) is therefore simply based on the gradient found at position l(k). We will refer to this hypothetical gradient of a noise-less system as the systematic gradient.

Results

While such Boolean learning has been applied to a wide range of computational problems, recurrent neural networks have a particular relevance for dynamical signal processing, and we therefore explore one-step-ahead prediction of the chaotic Mackey-Glass sequence with a Lyapunov exponent of ∼3•10 -3 . This particular input is a commonly employed benchmark test and our results are therefore directly comparable to other works such as Mackey-Glass prediction based on a semiconductor laser delay reservoir with weights optimized and applied in an offline procedure [START_REF] Bueno | Conditions for reservoir computing performance using semiconductor lasers with delayed optical feedback[END_REF], as well as the seminal work on RC [START_REF] Jaeger | Harnessing nonlinearity: predicting chaotic systems and saving energy in wireless communication[END_REF] where however the time step was twice as large.

The chaotic sequence acting as input information u(n+1) has zero mean and is normalized to its standard deviation, making error ϵ(k) the normalized mean square error. Its first two hundred points are used as training signal u(n+1), of which we however removed the first 30 time steps due to their transient nature. The result is a training signal with T = 200 -30 = 170 data points for which the target is τ(n + 1) u(n + 2). Finally, the testing error is determined with an independent set of 9000 data-points unused in the training sequence. Based on mutation of the readout weights, our concept explores an error landscape with height ϵ(k) and position W DMD (k), and reward r(k) drives the configuration from W DMD (1) to a local minima at W DMD (k min ). There our system will remain trapped due to an exploration step size of 1. We will refer to one complete learning process for k:1→k min as a minimizer.

Understanding why generalization is possible for a training set size (T = 170) not orders of magnitude larger than the number of to be optimized weights (N = 961) is an interesting question. Recent results on deep neural networks, triggered by the insightful analysis from Ref. [START_REF] Belkin | Reconciling modern machine learning and the bias-variance trade-off[END_REF], show that overparametrisation may not preclude generalization. See Ref. [START_REF] Hastie | Surprises in high-dimensional ridgeless least squares interpolation[END_REF] for an account to this phenomenon using random matrix theory, starting from simple linear models and generalizing to kernel estimation. In our setting we, however, might additionally postulate that we work below the overparametrization barrier due to the Boolean entries of W DMD (k), which brings substantial rigidity into play. The price one pays is making the problem harder from a computational optimization viewpoint [START_REF] Hadaeghi | Computing optimal discrete readout weights in reservoir computing is NP-hard[END_REF].

Typically, the main metric for evaluating learning are speed of convergence k min and final inference error ϵ min . However, in analogue neural hardware, reproducibility as well as robustness to noise and parameter drifts also play an essential role. We start by collecting statistical information and measure 20 [START_REF] Semenova | Fundamental aspects of noise in analog-hardware neural networks[END_REF] curves for the greedy (Markovian) exploration. All measurements started at an identical position W DMD (1) and we therefore focus on the algorithm's exploration of the error-landscape. Results are shown in Figure 2, with individual learning curves as grey lines and their average as red crosses. Panel (a) shows data for the greedy, panel (b) for Markovian exploration.

Average and local features of convergence and minima

On average, the error landscape topology excellently follows an exponential decay for both exploration strategies, see fit (blue line) to the average error (red crosses) in Figure 2. Comparing individual trajectories, however, reveals strong inter-trial differences significantly exceeding the noise level. This diversity corresponds to the error landscape's topological richness probed by the different random descents, and trajectories range from rather smooth descents to paths including steep drops. No correlation between the starting ϵ(1) and best performance ϵ min was found, and for the many different minimizers our system never got stuck in a local minima with bad performance. The different impact of a greedy or Markovian exploration sequence can be seen when comparing Figure 2a,b, respectively. Greedy exploration arrives at a minimum error of ϵ min (14.2 ± 2.9) ⋅ 10 -3 after k min = 973.6 ± 63.7 learning epochs, while Markovian exploration arrives at a slightly lower error ϵ min (13.4 ± 1.9) ⋅ 10 -3 at the expense of a prolonged learning effort k min = 1856.5 ± 175.1. Crucially, the system's testing error (green line in Figure 2) excellently matches its training error, hence ruling out over fitting. Noteworthy, convergence for both cases scales linear with network size

H(k a , k b ) ∑ i W DMD, a i (k a ) -W DMD, b i (k b ) .
For the 20 minimizers we obtain 20 (20-1)/2 = 190 distances between their respective minima, and their statistical distribution obtained for greedy exploration is shown in Figure 3. The red line is a Gaussian fit centred at H = 419 and with a half width at 1/e of 14. Data shows a very specific and unusual error landscape topology: local minima appear not to be irregularly distributed, nor located in a particular region. Instead, the negligible correlations between the minimas' locations, and the systematic and narrow distribution of inter-minima distances shown in Figure 3 reveals their almost uniform distribution across the error landscape. Again, we find that Markovian exploration results in an identical behaviour.

Noise sensitivity

To further investigate this phenomena, we reduce the number of uncertainties during learning. We measure three minimizer paths starting at the same W DMD 1 . One of the three minimizers acts as a master and defines mutation sequence l(k), k ∈ [1 … K], which the other two minimizers follow as slaves. Crucially, all three compute their own rewards r(k) and hence independently evaluate mutating the same weight. Keeping the potentially systematic error of a slow experimental parameter drift in mind, the three systems are evaluated at each learning epoch k before advancing to k + 1. A single minimizer takes ∼20 h, and sequential evaluation would amplify susceptibility to slow experimental parameter drifts which take place on the scale of hours in our experiment.

Results are shown in Figure 4a. The blue, green and red lines correspond to the different errors ϵ(k), plotted on a semi-logarithmic scale for the master and two slaves, respectively. The different data have a high degree of similarity with an average correlation of 99.4%, yet locally one can identify some significant differences. We computed the temporal evolution of the Hamming distances H(k) of the two slaves to their master (red and green crosses) and between the two slaves (grey crosses); all three grow linearly at essentially the same rate. Without noise, each minimizer's reward r(k) would be identical and they would consequently all follow the same trajectories W DMD (k) and arrive at the same minima W DMD (k min ).

To understand this behaviour, we therefore have to consider the impact of noise and learning upon the system's error ϵ(k). The response of error ϵ(k) to a modification in the system's output Δy out (k) is

Δϵ(k) ε(k) ⋅ Δy out (k), (12)
and Δy out (k) is the mean modification of output y out (k,n+1) within a certain window, which during training contains T sample points. Some general considerations regarding our system are in order. The amplitudes of all network nodes xi are Gaussian-distributed due to the SLM's illumination by a collimated Gaussian beam. Randomly changing one readout weight therefore results in Δy out (k) according to a normalized Gaussian distribution with a width of Δy out, learn (k) σ l (k). We have carefully characterized the noise of all elements in our opto-electronic NN, and their accumulated impact upon y out (k) is excellently approximated by Gaussian white noise with a width of Δy out, noise (k) σ n (k) [START_REF] Semenova | Fundamental aspects of noise in analog-hardware neural networks[END_REF]. Readout weights W DMD remain approximately evenly distributed between zeros and ones for all k. Learning does only modify readout connections and therefore neither modifies xi nor the system's noise, making both independent of learning. We can therefore assume that modifications to y out induced by learning and noise remain constant for all k, hence σ l (k) σ l and σ n (k) σ n . The fact that according to Eq. ( 12) noise (σ n ) and learning (σ l ) both modify the system's error according to the same relationship is of general importance. Convergence during learning is characterized by ϵ(k) and ε(k), and both depend on the particular computational tasks. Yet, the ratio between the error's susceptibility towards σ n and σ l remains constant as both scale with the same constant ε(k). This in turn imposes the same constant relative susceptibility of reward r(k) towards σ n and σ l , meaning that neither ε(k) nor ϵ(k) modify the interplay between noise and learning upon the system's weights. The following discussion, results and observation are fully general and independent of task and learning algorithm as long as these do not change amplitudes σ n and σ l . Noise and weight modifications are therefore independent players, whose action upon learning is somehow competitive.

The objective of modifying a readout weight is to probe the error landscape's systematic gradient. However, this action is contaminated by noise which can potentially exceed the systematic gradient in the opposite direction. The consequence is a change in the sign of Δϵ(k), in which case reward r (k) is inverted. How likely such a modification takes place depends on the relative amplitudes of σ n and σ l , and C is the constant probability of such a modification occurring. The analytical derivation of C is possible, yet beyond the scope of this manuscript.

Probability C is the driving force behind the growing separation between two identical minimizers, and two situations are relevant. The first situation occurs when r(k) for one minimizer is inverted by noise while the other preserves its systematic value, which has a probability of

C(1-C)+(1-C)C = 2C(1-C).
The other situation is if both minimizers have an identical reward r(k), which can either be the consequence of both retaining their systematic result, or for both being inverted by noise, with a combined probability of C(1-C) 2 +C 2 = 1-2C(1-C). The first situation leads to H(k + 1) ≠ H(k), the second to H(k + 1) = H(k), and the Hamming distance's rate equation is

ΔH(k + 1) ρ id (k)2C(1 -C) -ρ op (k)2C(1 -C). (13)
Here, ρ id (k) and ρ op (k) are the probability of finding both minimizers' weights l(k) to be identical or opposite, respectively. Using ρ id (k) 1 -ρ op (k), we arrive at

H(k + 1) H(k) + C 1 -2ρ op (k) , (14)
ΔH(k + 1) C 1 -2ρ op (k) , (15)
where C 2C(1 -C).

The Hamming distance's evolution is therefore governed by noise quantified through constant C, and by how the learning algorithm picks weight W DMD l(k) (k) from a population with a certain ρ op (k), and neither error nor gradient play a role.

For fully random mutation, the probability of a weight to be selected is identical at every k, and hence the Hamming distance at the previous epoch k determines the probability of two weights being opposite in their configuration: ρ op (k) H(k)/N. For greedy learning, the bias term in Eq. ( 3) causes mutations hardly ever to repeat the same weight during an interval specified by k k ' + aN, k ' ∈[1, N], with non-negative integer a. The probability of both minimizers to be configured opposite for all k' and a specific a is therefore their Hamming distance at the end of the previous interval: ρ op (k) H(aN)/N, which results in constant slopes ΔH(a) for each k'.

Figure 4b shows the evolution of Hamming distance H(k/N), and since all minimizers start at the same position W DMD (1) we always have H(1) = 0. Greedy mutations in the experiment (analytics) are the red line (black dashed line), while random mutations in the experiment (analytics) are the blue line (black solid line). For both scenarios, greedy and random descent, the experimental data is the average obtained from 20 minimizers. We then changed the starting conditions and realized two parallel minimizers which started with a separation H(1) > 0, see Figure 4c. In general, the evolution according to Eq. ( 15) perfectly reproduces results of the highly different experimental learning scenarios. In particularly for the averaged data, where we always arrive at H(k)| k→∞ N/2, regardless of the algorithm.

Different minimizers therefore always arrive at final readout configurations which share no common feature. This suggests a closer look into the role and relevance of individual weights: how many induce a systematic contribution to convergence at all, and if their gradients depend on the sequence of previous mutations. We optimized readout weights via two minimizers starting at different random positions W DMD,a ( 1 Weights insensitive (sensitive) to preceding optimizations correspond to linearly independent (linearly dependent) NN dimensions. Linearly independent NN dimensions must always induce the same gradient, regardless of the preceding optimization path, and they therefore have to be located on the red diagonal line in Figure 5b. The Figure's green area indicates the linearly independent criteria when considering the impact of noise σ n , and we find ≈30% of 430 dimensions fall into this category. However, this is only a necessary criterion; a sufficient criteria requires allocating dimensions inside this area for all potential configurations of the remaining m-1 dimensions, which for the 2 429 possibilities is prohibitive to prove experimentally. The NN dimensions whose weight configuration depends on the previously optimized weights lie outside the grey and green areas. This is a sufficient criteria for linear dependent NN dimension, and we find ≈59% of 430 to be contained inside this set. There appears to be no concentration of potentially linearly independent dimensions towards the red diagonal line, and hence we conclude that these also mostly belong to the set of linear dependent dimensions.

Discussion

Our experimental findings and analytical descriptions are the first of their kind and stimulate a fundamental discussion. Equation ( 12) is of interesting consequence for noisy hardware NNs comprising linear readout weights. It links the susceptibility of Δϵ(k) to NN noise to the system's location inside the error landscape ϵ(k). Experimentally we obtained noise induced variations ε(150)σ n 4.4 ⋅ 10 -3 and ε(961)σ n 0.6 ⋅ 10 -3 . Our learning curves excellently agree with exponential convergence ϵ(k) ϵ (0) e -αk , for which ε(k) ∝ ϵ(k). Error and gradient therefore evolve in a linearly proportional manner, and ϵ(150)/ϵ(961) 9.8 in close agreement with the noise sensitivity's evolution ε(150)σ n / ε(961)σ n ≈ 7.3 confirms this fundamental relationship.

We would like to also propose alternative noise-mitigation approaches which are a derivative of our findings. Simply suppressing noise on a hardware level is potentially expensive, and topological requirements can limit mitigation based on connectivity statistics [START_REF] Semenova | Fundamental aspects of noise in analog-hardware neural networks[END_REF]. One might therefore curb the impact of noise by modified learning strategies. Noise will first of all limit the absolute performance, but also cause this performance to fluctuate, which is an effect one could address for example by amending an optimization's cost function by the gradients encountered in the proximity of a neighborhood. According to Eq. ( 12) the local gradients ε, or ε| W DMD (k) probed during learning also determine the system's performance fluctuation. Making them part of the optimization and not only the absolute error would therefore stabilize the system's performance. The influence of the noise on the performance of the algorithm can be mitigated using aggregated information along the trajectory in the spirit of recent work about variance reduction in stochastic gradient schemes and its extensions to non-convex zeroth-order optimization; see e.g. [START_REF] Liu | Zeroth-order stochastic variance reduction for nonconvex optimization[END_REF].

Equation [START_REF] Hermans | Embodiment of learning in electro-optical signal processors[END_REF] shows that for C > 0 the Hamming distance between readout weights of two systems will always tend towards complete decorrelation as H(k)| k→∞ N/2. Even for 100% identical networks one will therefore never obtain similar readout configurations [START_REF] Freiberger | Training passive photonic reservoirs with integrated optical readout[END_REF]. This finding can most likely be extended to the non-Boolean case and to the weights between layers of analogue deep NNs. The field of learning implemented in physical and hence noisy substrates is only in its infancy [START_REF] Shen | Deep learning with coherent nanophotonic circuits[END_REF][START_REF] Antonik | Brain-inspired photonic signal processor for generating periodic patterns and emulating chaotic systems[END_REF][START_REF] Bueno | Reinforcement learning in a large scale photonic recurrent neural network[END_REF], and confirmation of our findings in other hardware systems would prove the generality of our result. Finally, the human brain is a very noisy network indeed [START_REF] Suarez-Perez | Quantification of signal-to-noise ratio in cerebral cortex recordings using flexible MEAs with co-localized platinum black, carbon nanotubes, and gold electrodes[END_REF], which suggests that our findings may also have interesting implications for the field of theoretical neuroscience.

We have shown that the large majority of our NN's dimensions are most likely linear dependent. What this means in pratical terms is that each modification of a weight has to be interpreted in the context of all previous modifications. Each configuration W DMD (k) therefore encodes the history of modifications to the reward due to noise during the previous learning epoch.

One direct consequence for applications is that one cannot simply transfer or swap weight configurations between optimized analogue neural networks, even for potentially available identical twin networks. The reason is that optimized configurations are not only the consequence of error landscape, system properties and noise, but also of the precise history of noise during an exploration path. Even almost perfectly reproducible hardware networks will therefore always have to be individually trained for optimal performance; simply uploading a configuration will potentially not work. A 'school' in which each neural network learns individually might therefore be required. Finally, our findings open a new field where such twin-minimizers could be considered for probing and interrogating unknown hardware neural networks. The average divergences shown in Figure 4b agree exceptionally well with our model, and based on such data one can therefore make accurate inferences about the noise properties of a hardware NN and about its error landscape exploration strategy.

Conclusions

In our work we have investigated the intricate interactions between different learning concepts and the noise inherently present in analogue neural networks. We experimentally showed that trajectories of individual minimizers (i.e. learning trajectories) strongly diverge, and were able to analytically link this divergence to a constant ration between output error and noise susceptibility. Our analytical description only assumes a linear multiplication between a NN's state and its readout weights, and hence should be generally applicable to this wide class of analogue hardware NNs.

Figure 1 :

 1 Figure 1: (a) Schematic illustration of a recurrent neural network. (b) Photonic implementation of a spatio-temporal neural network with 961 nodes. An optical plane wave illuminates the spatial light modulator (SLM), the neural network states are encoded by the SLM pixels. These are imaged on the camera, passing through a polarizing beam splitter (PBS) and the diffractive optical element (DOE) creating the coupling between network states. The information detected by the camera is used to drive the SLM. The network's output weights are realized via a digital micro-mirrors device (DMD) which creates Boolean readout weight matrix W DMD .

Figure 2 :

 2 Figure 2: Learning performance, with individual (averaged) trajectories shown by grey data (red crosses). The green line is the testing error, and the blue is an exponential fit. Panel (a) and (b) were obtained with greedy and Markovian exploration, respectively.

Figure 3 :

 3 Figure 3: Probability distribution of the Hamming distances at k = N, obtained from the different greedy learning curves displayed in Figure 2a. Red curve is a Gaussian fit which is centred at 419 and has a standard deviation of 29.

Figure 4 :

 4 Figure 4: (a) Three minimizers starting from identical position are measured in parallel. The red and green (slaves) optimization paths test the dimension l(k) determined by the blue minimizer (master), reward of each mutation are evaluated and applied individually. The lines show the individual errors on a logarithmic scale, crosses depict the Hamming distance's evolution between the three minimizers. (b) Hamming distance evolution with learning epoch k normalized by the network's Size N. Random mutation (blue line) leads to a smooth saturation function behaviour, while biased (greedy) mutation (red dashed line) results in linear intervals of length N. (c) The same characteristics are found for two minimizers starting already separated by a distance (1)≠0, and the behaviour is therefore generally true for Boolean learning in noisy hardware.

  ) and WDMD,b [START_REF] Lecun | Deep learning[END_REF], which arrived at two distinct local minima M a W DMD, a (k min, a ) andM b W DMD, b (k min, b). Once there, we determined the list of m weights where M b differs from M a . The list is randomly arranged in sequence l ∈ [l(1), l(2), …, l(m -1), l(m)] according to which we invert the corresponding weightsW DMD, a l(k) (k) and W DMD, b l(k) (k) for k ∈ [1, m].Importantly, this mutation is always kept and no optimization based on reward r (k) is taking place. Starting from M a (M b ), this results in a random path P a : M a → M b (P b : M b → M a ). As the weights addressed in sequence l(k) are the ones in an opposite configuration for M a and M b , P a and P b connect both minima along inverted trajectories, see Figure5a. We probe error ϵ(k) along P a and P b and determine error gradientsεa (k) ϵ a (k) -ϵ a (k min, a ) and εb (k) ϵ b (k) -ϵ b (k min, b ).Results are shown in Figure5bin the ( εa , εb )-plane, and for this experiment we obtained m = 430 different dimensions between M a and M b . Only ≈11% of the 430 gradients consistently remained below our system's noise floor εσ n , indicated by the grey circle.

Figure 5 :

 5 Figure 5: (a) Two inverted paths are probed between two local minima. (b) Error gradient for all readout weights encountered along path one and two as x and y axis, respectively. The red diagonal corresponds to linear independent weights, and the uncertainty induced by noise is indicated by the green area. Within, weights are potentially linear independent. Outside the green area weights are linearly dependent, and for data inside the grey circle of diameter σ no classification is possible.

L.Andreoli et al.: Boolean learning under noise-perturbations

Acknowledgement: The authors acknowledge the support of the Region Bourgogne Franche-Comté. This work has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreements No. 860360 (POST DIGITAL) and No. 713694 (MULTIPLY). This work was also supported by the EUR EIPHI program (Contract No. ANR-17-EURE-0002), the BiPhoProc project (Contract No. ANR-14-OHRI-0002-02), and by the Volkwagen Foundation (NeuroQNet). Author contribution: All the authors have accepted responsibility for the entire content of this submitted manuscript and approved submission. Research funding: This research was funded by the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreements No. 860830 (POST DIGITAL) and No. 713694 (MULTIPLY). This work was also supported by the EUR EIPHI program (Contract No. ANR-17-EURE-0002), the BiPhoProc project (Contract No. ANR-14-OHRI-0002-02), and by the Volkwagen Foundation (NeuroQNet).

Conflict of interest statement: The authors declare no conflicts of interest regarding this article.