Smail Ait Bouhsain

Rachid Alami

Thierry Siméon

Simultaneous Action and Grasp Feasibility Prediction for Task and Motion Planning through Multi-Task Learning

In this paper, we address task and motion planning (TAMP) which is an important yet challenging robotics problem. It is known to suffer from the high combinatorial complexity of discrete search, often requiring a large number of geometric planning calls. We build upon recent works in TAMP by taking advantage of learning methods to provide action feasibility information as a heuristic to the symbolic planner, thus guiding it to a geometrically feasible solution and reducing geometric planning time. We propose AGFP-Net, a multi-task neural network predicting not only action feasibility, but also the feasibility of a set of grasp types. We also propose an improved feasibility-informed TAMP algorithm capable of solving more complex problems, and handling goals which are not fully specified. Comparative results obtained on different problems of varying complexity show that our method is able to greatly reduce task and motion planning time.

I. INTRODUCTION

Task and motion planning (TAMP) [START_REF] Alami | A geometrical approach to planning manipulation tasks. the case of discrete placements and grasps[END_REF]- [START_REF] Ait Bouhsain | Learning to predict action feasibility for task and motion planning in 3d environments[END_REF] is a robotics problem that involves determining a sequence of actions that a robot must take to achieve a desired goal, along with the corresponding motions of the robot. It combines discrete symbolic planning with continuous geometric planning. TAMP problems are challenging due to the high dimensionality of the state and action spaces, but also to the combinatorial complexity of discrete search.

Geometric planning is a major bottleneck to TAMP. Since symbolic planners do not have any geometric reasoning capabilities, they generate task plans without a guarantee of feasibility. A geometric planner has to be used in order to verify the feasibility of each generated task plan and construct the corresponding motions. However, in complex problems, symbolic planners generate a high number of geometrically infeasible solutions before finding a feasible one, each plan requiring a call to the geometric planner which results in long planning time. Moreover, geometric planning might be time consuming even in the case of feasible tasks. Since not all grasps are feasible or lead to a feasible motion, finding the right one might require a lot of time.

This work aims at providing the symbolic planner with a feasibility prediction capability, which increases the probability that the generated task plans are feasible, and thus reduces the number of calls to the geometric planner. First, we improve upon our previous approach [START_REF] Ait Bouhsain | Learning to predict action feasibility for task and motion planning in 3d environments[END_REF] in order to enrich the information provided to the task planner, and reduce the dimensionality of the action space. Indeed, instead of predicting the feasibility of a P ick or P lace action with a single specific grasp type, we propose a multi-task learning neural network capable of predicting the feasibility of an action as well as the feasibility of each grasp type. We also propose a new training strategy allowing a better generalization to new environments. In addition, we propose a more powerful TAMP algorithm capable of solving a wider range of problems with higher complexity, in particular ones with semi-specified goals. We develop a method for combining probabilities of feasibility allowing the task planner to make better use of the predictions obtained by the neural network, and reduce the branching factor of the search. We also introduce a new cost function that incorporates both the action feasibility as well as the grasp feasibility predictions. Furthermore, we leverage the latter to accelerate geometric planning, thus reducing the time spent on finding a grasp that leads to a feasible motion.

II. RELATED WORK

Task and motion planning combines discrete symbolic planning and continuous geometric planning. Early works [START_REF] Alami | A geometrical approach to planning manipulation tasks. the case of discrete placements and grasps[END_REF]- [START_REF] Lamiraux | Prehensile manipulation planning: Modeling, algorithms and implementation[END_REF] view the problem from the geometrical perspective, by formalizing it as a multi-modal motion planning problem.

Most existing TAMP methods [START_REF] Cambon | A hybrid approach to intricate motion, manipulation and task planning[END_REF]- [START_REF] Garrett | Integrated task and motion planning[END_REF] combine a task planner with a geometric planner, while using geometric backtracking as an interface between the two. These methods are known to suffer from the high combinatorial complexity of discrete search and the slow geometric planning time, since a call to the geometric planner is still needed to verify the feasibility of symbolic plans.

Recent works in TAMP [START_REF] Carpentier | Learning feasibility constraints for multi-contact locomotion of legged robots[END_REF]- [START_REF] Ait Bouhsain | Learning to predict action feasibility for task and motion planning in 3d environments[END_REF] propose to accelerate the planning process by leveraging learning methods. They either aim at guiding the planning process towards a solution, or reducing the number of calls to the geometric planner by providing a learned heuristic to the task planner. Learning approaches have also been proposed for grasp planning [START_REF] Mahler | Dex-net 1.0: A cloud-based network of 3d objects for robust grasp planning using a multi-armed bandit model with correlated rewards[END_REF]- [START_REF] Wu | Grasp proposal networks: An end-to-end solution for visual learning of robotic grasps[END_REF], these works aim at finding possible grasps of complex objects, which is a step of TAMP often overlooked. Works such as [START_REF] Wells | Learning feasibility for task and motion planning in tabletop environments[END_REF] [24] [START_REF] Driess | Deep visual reasoning: Learning to predict action sequences for task and motion planning from an initial scene image[END_REF] propose to train a learning model to predict the geometrical feasibility of P ick and P lace actions in tabletop environments. Bouhsain et al. [START_REF] Ait Bouhsain | Learning to predict action feasibility for task and motion planning in 3d environments[END_REF] extend this idea to 3D environments and proposes a feasibilityinformed task and motion planner which takes advantage of feasibility prediction. This work improves upon [START_REF] Ait Bouhsain | Learning to predict action feasibility for task and motion planning in 3d environments[END_REF] by proposing a multi-task learning neural network providing richer geometric feedback to the task and motion planner. We also develop a TAMP algorithm which uses a new method for handling feasibility predictions, and is capable of solving more complex problems, for instance ones with semi-specified goals.

III. PROBLEM DESCRIPTION

We focus on manipulation problems involving the rearrangement of a set of objects in a three-dimensional manipulation environment E. The latter is composed of a single robot arm with a parallel gripper, n O box-shaped movable objects, n ss stable support surfaces and n obs fixed obstacles. A state s of the planning scene is defined by the configuration of the robot and of all movable objects, i.e the support surface and the pose of each object. We denote the configuration of a movable object O in the state s as s(O). An action a is defined as the transition between two distinct states s and s ′ . It involves picking a single movable object O from its configuration in s, and placing it at a new configuration q. The action a can then be explicitly formulated as:

a = M ove O (s(O) → q) (1)
Finding a solution to a TAMP problem consists in finding a sequence of actions τ = {a 0 , a 1 , ..., a K }, as well as the corresponding sequence of motions Π = {π 0 , π 1 , ..., π K }, bringing the environment from an initial state s 0 to a goal state s goal . Each movable object might have either a fullyspecified goal placement (i.e a support surface and a pose), or a semi-specified goal placement (i.e a support surface only), or no specified goal at all.

We also define G as the set of 6 grasp types described in [START_REF] Ait Bouhsain | Learning to predict action feasibility for task and motion planning in 3d environments[END_REF]. They correspond to the side from which the object is grasped in the perspective of the robot. Hence, each grasp type is a subspace of grasps from the corresponding side of the object.

IV. ACTION AND GRASP FEASIBILITY PREDICTION

As in our previous work [START_REF] Ait Bouhsain | Learning to predict action feasibility for task and motion planning in 3d environments[END_REF], we aim at reducing the number of calls to the geometric planner by providing a feasibility prediction capability to the symbolic planner, allowing it to generate task plans with a high probability of feasibility first. In this paper, rather than predicting the feasibility of an action with a specific grasp type, we train a neural network to simultaneously predict the general feasibility of a P ick or P lace action without considering the grasp, as well as the feasibility of each grasp type. This means that, contrarily to the previous method, the grasp is no longer an input to the neural network but rather an output.

A. Neural Network

The proposed neural network takes as input the 3D scene and the action to be tested. We use the same representation of the 3D environment as in [START_REF] Ait Bouhsain | Learning to predict action feasibility for task and motion planning in 3d environments[END_REF], using 5 depth images which correspond to different views of the scene (i.e top, front, rear, right, left). Each view corresponds to a specific plane in the 3D space. The top view shows a projection of the environment on the (x, y) plane, the front and rear views correspond to the (y, z) plane, and the right and left views correspond to the (x, z) plane.

Regarding action representation, we consider P ick and P lace tasks as equivalent. Indeed, if we construct these actions as starting and ending at the same home configuration of the robot, the motion for placing an object at a configuration is the same as the reversed motion of picking the same object from the same configuration. Using this assumption 1 , we always consider the action to be tested as a P ick task, which reduces the number of inputs to the neural network, and allows to speedup data generation and annotation (Section IV-B).

As in [START_REF] Ait Bouhsain | Learning to predict action feasibility for task and motion planning in 3d environments[END_REF], the object of interest is represented using 5 masks over the scene views, showing only the object at the configuration it should be picked or placed at. Instead of combining the scene views and the object masks into a single 10-channel image, we combine them into 3 separate images depending on the plane the views correspond to. The first is a 2-channel image I xy consisting of the top scene view and mask. The second and third are two 4-channel images I xz and I yz , consisting of the left-right and the frontrear scene views and masks respectively 2 . These images are encoded using three ResNet-based convolutional neural networks [START_REF] He | Deep residual learning for image recognition[END_REF]. The resulting embeddings are concatenated, then fed to a fully connected layer to obtain an encoding of the 3D scene and the action. A multi-layer perceptron followed by a Sigmoid activation function are then used to obtain 2 predictions. The first one is the probability of feasibility of the action p F . The second is a vector p G composed of 6 values, each one being the probability of feasibility of one of the grasp types in G. Figure 2 shows the complete architecture of our proposed neural network.

B. Data Generation and Annotation

In order to train the proposed neural network, we develop an automatic data generation method which can be divided into three steps. The first is scene generation, which follows the same procedure detailed in [START_REF] Ait Bouhsain | Learning to predict action feasibility for task and motion planning in 3d environments[END_REF]. It generates scenes containing a large table, up to 4 elevated support surfaces, and 2 objects. The second is the construction of the dataset. For each object O in a generated scene, we generate one P ick action from its configuration q. Using the assumption that P ick and P lace actions are equivalent, these generated P ick actions also count as P lace actions at the same configuration q.

The third step of data generation is annotation. In this step, we need to label each generated action but also each grasp type as feasible or infeasible. We consider a grasp type as feasible if at least one grasp corresponding to this type is feasible. All the generated actions causing a collision are guaranteed to be infeasible, and are thus automatically labeled as such. We annotate the rest of the actions using a geometric planner composed of a grasp sampler, an inverse kinematics solver and a sampling based motion planner.

V. TAMP ALGORITHM Compared to our prior work [START_REF] Ait Bouhsain | Learning to predict action feasibility for task and motion planning in 3d environments[END_REF], we propose an improved TAMP algorithm, capable of solving a wider range of more complex problems. In addition to leveraging feasibility predictions to accelerate the planning process, the branching factor of our task planner is reduced by using a different definition of actions at the symbolic level. While in our previous algorithm, an action was defined as either a P ick or a P lace task with a specific grasp type, in this work, we define actions as combined P ick -P lace tasks without any grasp type specified.

The main algorithm remains similar to our previous planner. It combines a symbolic task planner which finds a sequence of actions bringing the environment to a goal state, and a geometric planner which checks the feasibility and finds the motions corresponding to each action. The task planner is a best-first tree search in which nodes are states of the environment and edges are actions, as defined in Section III. At each iteration, the node s with the lowest cost is retrieved from the set of open nodes. If s is a goal state, we reconstruct the complete action sequence τ leading to s. The geometric planner is then called to find the corresponding sequence of motions Π. If geometric planning succeeds, then a geometrically feasible solution was found and the problem is solved. On the other hand, if one of the actions is infeasible, we prune all children of the first infeasible node in the task plan from the tree, and continue the search. child ← nextState(s, a)

6:
child.cost ← computeCost(child, s goal , p F (a))

a ← M ove O (s(O) → s goal (O)) 6: A ← A ∪ a 7: else if s goal (O) is semi-specified then 8:
R ← sampleP lacements(E, s goal (O))

9:
for each q in R do 10:

a ← M ove O (s(O) → q)

11:

A ← A ∪ a 12:

end for R ← sampleRandomP lacements(E)

16:
for each q in R do 17:

a ← M ove O (s(O) → q) 18: A ← A ∪ a 19:
end for 20: end for 21: return A In case s is not a goal node, we call the f indChildren function described in Algorithm 1. We first proceed to finding the actions applicable at s by calling the function f indP ossibleActions detailed in Algorithm 2. For each movable object O in the scene, we generate two types of actions. The first action type tries to move the object O from its placement at s denoted s(O), directly to its goal placement s goal (O). Within this type of actions, there are two different scenarios. If O has a fully-specified goal (i.e support surface and pose), then only one action is generated that moves the object to that goal (Alg. 2, lines 4-6). However, if O has a semi-specified goal (i.e it has to be placed on a specific support surface whatever the pose), then the number of goal placements is infinite. In this case, we sample a fixed number of placements on the goal support surface. We then generate multiple actions, each one moving object O to one of the sampled placements (Alg. 2, lines 7-12).

The second action type aims at moving the object O to a random placement without considering the goal. We sample a fixed number of random placements, then generate multiple actions moving object O to each one of them (Alg. 2, lines [START_REF] Silva | A new approach to combined symbolic-geometric backtracking in the context of human-robot interaction[END_REF][START_REF] Gharbi | Combining symbolic and geometric planning to synthesize human-aware plans: toward more efficient combined search[END_REF][START_REF] Bidot | Geometric backtracking for combined task and motion planning in robotic systems[END_REF][START_REF] Garrett | Ffrob: Leveraging symbolic planning for efficient task and motion planning[END_REF][START_REF] Garrett | Sampling-based methods for factored task and motion planning[END_REF]. In both action types, the number of placements sampled is a parameter chosen by the user 3 For each possible action a found, we compute its probability of feasibility p F (a) as well as the probability of feasibility of each grasp type p G (a) (Alg. 1, line 5). Since the neural network predicts the feasibility of P ick and P lace tasks separately, we need a strategy for combining these predictions in order to obtain the probability of feasibility of the complete P ick -P lace action. This is done in the predictF easibility function described in Algorithm 3. First the action a is decomposed into a P ick task and a P lace task. Then AGFP-Net is queried twice to obtain the probability of feasibility and the grasp types probabilities for the P ick and P lace actions separately. For a single grasp type g, we define the probability p g (a) that the complete action a is feasible using g as follows:

p g (a) = p g (P ick) × p g (P lace)

with p g (P ick) and p g (P lace) being the probability that the P ick and the P lace tasks, respectively, are feasible using g.

Using this definition, we can obtain the complete vector of grasp types probabilities p G (a):

p G (a) = p G (P ick) ⊗ p G (P lace) (3)
where ⊗ represents the element-wise product.

Regarding the probability of feasibility of the action p F (a), multiplying the probabilities of the P ick and P lace tasks is not sufficient. The reason is that even if both these tasks are feasible, the grasps leading to a feasible P ick might not be the same as the ones leading to a feasible P lace. Indeed, the complete action a is feasible iff. both the P ick and the P lace tasks are feasible, and have at least one feasible grasp type in common. Therefore, we define p F (a) as:

p F (a) = p F (P ick) × p F (P lace) × max(p G (a)) (4)
where p F (P ick) and p F (P lace) are the probabilities of feasibility of the P ick and the P lace tasks respectively, and max(p G (a)) is the highest probability in the vector p G (a).

After the feasibility prediction step, we construct a new child as the result of applying action a to state s (Alg. 1, line 5). We then compute its cost which is defined as:

C T otal = C SoF ar + C T oGoal + C F easibility (5)
where C SoF ar is the number of actions in the branch leading to the child node, C T oGoal is the number of objects that are not at their goal placement, and C F easibility is a feasibility cost calculated using the expression:

C F easibility = 1 p F (a) -1 ∈]0, ∞[ (6)
This definition allows to incorporate the feasibility predictions in the state costs. Thus, the higher the probability of feasibility of an action leading to a state is, the lower its cost will be. After all the children are generated, they are added to the set of open nodes which is sorted by increasing cost.

One important novelty of this work is the use of grasp feasibility predictions to accelerate geometric planning. When the geometric planner is called to find the motion of an action a, we use the previously computed grasp type probabilities p G (a) in order to prioritize feasible grasp types during grasp sampling. Thanks to this heuristic, the geometric planning time on feasible actions can be reduced.

VI. EXPERIMENTS

In order to test our method, we first train and test the proposed neural network on a generated dataset. Then, we define a set of TAMP problems to evaluate the performance of our algorithm.

A. Neural Network Training and Testing

Using the data generation method described in IV-B, we construct a training dataset based on 20 ′ 000 generated scenes. We also generate a validation set and a testing set, based on 5 ′ 000 generated scenes each, in order to evaluate the performance of the trained neural network. For annotation, we use an adapted version of Moveit! Task Constructor [START_REF] Görner | Moveit! task constructor for task-level motion planning[END_REF] with a BiTRRT motion planner [START_REF] Devaurs | Enhancing the transition-based rrt to deal with complex cost spaces[END_REF]. The model is implemented on Pytorch [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF] and trained using a weighted binary cross entropy loss to account for the imbalance between feasible and infeasible datapoints. We use the ADAM optimizer with a learning rate of 0.0001 and a batch size of 128. To prevent overfitting, we use dropout with a probability of 0.1, and weight decay of 0.0001.

Moreover, in order to increase the size of our training set, we use online data augmentation. Since the robot can turn 360 • around its base, rotating the scene around the z axis does not affect the feasibility of the action or the grasp types. This allows to obtain new labeled datapoints from already annotated ones. During training, scenes are either kept as they are, or rotated by a multiple of π 2 around the z axis. This data augmentation method improves training by allowing the neural network to learn more scenarios, and ensuring that the predictions are consistent between equivalent scenes. It also allows for a better generalization to new environments. The neural network is trained for 200 epochs, which takes approximately 24 hours.

B. Test Problems for the TAMP Algorithm

In order to compare the performance of our proposed approach with our prior work, we reuse the 5-object versions of the three TAMP problems described in [START_REF] Ait Bouhsain | Learning to predict action feasibility for task and motion planning in 3d environments[END_REF]. The first is the Reorder domain in which the goal is to move a set of movable objects from a shelf to another while changing the order. The second is the U npack problem, in which a set of objects on a tray have to be organized on a shelf. The third is the Swap domain in which the goal is to swap the placements of different objects. In these problems, the challenge comes either from grasp choice, the proximity between objects, or the occupancy of goal placements.

We also define two more complex TAMP problems that demonstrate the additional capabilities of our new method. The Access domain, shown in Figure 1a, involves moving a single object (red) to its goal placement while ensuring that the other movable objects are back at their initial placement. In this problem, the challenge comes from the fact that each object blocks access to the object next to it. The robot has to move all blocking objects one by one to a temporary placement, move the red object to its goal placement and finally put back the previously moved objects at their initial placements in the same order.

The Sort problem, shown in Figure 1b aims at demonstrating the ability of our algorithm to handle semi-specified and unspecified goals. The goal is to move a set of objects from a shelf and to one of two small tables depending on their color. Here, the challenge comes from placement sampling, meaning that the planner must find the set of placements allowing all objects of the same color to be on a narrow support surface. The initial proximity of objects is also a challenge, since the objects have to be moved in a specific order. Moreover, we make the problem more complex by adding blocking objects (orange) with unspecified goals on the tables.

The Access and Sort problems also aim at testing the generalizability of AGFP-Net to more realistic scenes, with a higher number of objects and vertical obstacles.

VII. RESULTS

We first analyze the performance of AGFP-Net on the generated test data. We then evaluate the performance of the proposed TAMP algorithm on the test TAMP problems.

A. Neural Network Performance

Table I shows the performance of the proposed neural network. Comparing the results obtained by AGFP-Net to the one yielded by [START_REF] Ait Bouhsain | Learning to predict action feasibility for task and motion planning in 3d environments[END_REF] shows an overall improvement in performance with an F1-Score of 92.4%. Also, the yielded true positive rate (TPR) of 93.7% and true negative rate (TNR) of 97% show that AGFP-Net is as good at predicting feasibility as it is at predicting infeasibility. Regarding grasp feasibility prediction, the results show that the model is able to accurately predict the feasibility of each grasp type with an F1-Score of 86.6%, a TPR of 93.3% and a TNR of 97.1%.

B. Performance of the feasibility-informed TAMP algorithm

We ran our algorithm on all five TAMP problems defined in Section VI-B with 10 runs each. In order to measure the performance gain of action and grasp feasibility prediction, we do the same with a non-informed version of our TAMP algorithm. This is done by disabling feasibility prediction and setting all probabilities of feasibility to 1. We first compare our approach to the one proposed in [START_REF] Ait Bouhsain | Learning to predict action feasibility for task and motion planning in 3d environments[END_REF]. Figure 3 and Table II show a comparison of the two methods on the Reorder, U npack and Swap problems. Results obtained without feasibility prediction show that the new TAMP algorithm generates less geometrically infeasible task plans, which translates into less geometric planning time spent on infeasible actions. This is due to the fact that grasp types are no longer considered by the task planner, but rather by the geometric planner. As a result, in problems where the grasp choice is important such as the Reorder and U npack domains, the task planner does not generate task plans that are infeasible due to the choice of grasp. In the Swap problem, however, this gain in performance is not For each method, we report the mean and standard deviation of the planning time with and without feasibility prediction, using AFP-Net for [START_REF] Ait Bouhsain | Learning to predict action feasibility for task and motion planning in 3d environments[END_REF] and AGFP-Net for the proposed approach.

as pronounced since in this domain, infeasibility is mostly due to collisions. On the other hand, comparing the performances yielded using the feasibility prediction heuristic shows that thanks to the improved accuracy of AGFP-Net (producing less misclassifications), and the new probability combination strategy, the proposed TAMP algorithm is able to reduce the number of infeasible task plans generated to 0 in most runs, even in problems where [START_REF] Ait Bouhsain | Learning to predict action feasibility for task and motion planning in 3d environments[END_REF] has some failures such as the U npack and the Swap domains. Figures 3b and3c show that geometric planning time on infeasible actions is eliminated using our feasibility-informed planner.

Moreover, the results demonstrate the improvement in planning time yielded using grasp feasibility prediction. Indeed, comparing geometric planning time on feasible actions with and without AGFP-Net shows that using the grasp feasibility heuristic allows the geometric planner to find grasps leading to a feasible motion faster. This decrease in geometric planning time not only compensates for the time spent on feasibility prediction, but also reduces the total planning time as shown particularly in Figure 3a.

On the Access problem, the method proposed in [START_REF] Ait Bouhsain | Learning to predict action feasibility for task and motion planning in 3d environments[END_REF] fails to find a geometrically feasible solution in a reasonable amount of time, whether feasibility prediction is used or not. This is due to the high combinatorial complexity of the problem. Thanks to the reduced branching factor, our algorithm is able to solve this problem in approximately 25 minutes when no feasibility heuristic is used. However, Table II shows that the number of infeasible task plans generated before finding a feasible solution is high with an average of 339.3. Using AGFP-Net, this number is reduced to an average of 1.9, yielding a 93% reduction in total planning time. For the Sort domain, the number of infeasible task plans generated is 108, which is reduced to 1 using AGFP-Net. Hence, the total planning time is 12 times faster.

Furthermore, these results show that, although AGFP-Net has been trained on scenes containing two objects only and no vertical obstacles, it is able to generalize to scenes with a higher number of objects and a limited number of obstacles while keeping feasibility predictions accurate.

VIII. CONCLUSION

In this work, we propose a feasibility-informed task and motion planner which improves total planning times on three fronts. First, the branching factor of the task planner is reduced allowing better handling of the combinatorial complexity of the search. Second, we take advantage of feasibility predictions to generate geometrically feasible task plans. Finally, we leverage grasp feasibility predictions to accelerate geometric planning time on feasible actions. We also demonstrate the performance gain yielded by our method on 5 different TAMP problems.

As future work, we consider the extension of this approach to multi-robot TAMP problems such as ones involving multiple robot arms, or mobile manipulators. Also, It is also important to test and generalize the proposed method to reallife scenarios where the shapes of objects, support surfaces and obstacles might be more complex.

  Fig. 1: A visualization of the initial and goal states of two complex TAMP problems solved by the proposed planner, and on which Bouhsain et al. 2023 [30] fails.

Fig. 2 :

 2 Fig. 2: The architecture of the proposed neural network AGFPNet.

Algorithm 1 findChildren

 1 Input: s, E, s goal 1: children ← ∅ 2: A ← f indP ossibleActions(s, E, s goal ) 3: for each a in A do 4: [p F (a), p G (a)] ← predictF easibility(s, E, a) 5:

Fig. 3 :

 3 Fig.3: Planning time comparison between the proposed approach and Bouhsain et al. 2023[START_REF] Ait Bouhsain | Learning to predict action feasibility for task and motion planning in 3d environments[END_REF] averaged over 10 runs. For each method, we report the mean and standard deviation of the planning time with and without feasibility prediction, using AFP-Net for[START_REF] Ait Bouhsain | Learning to predict action feasibility for task and motion planning in 3d environments[END_REF] and AGFP-Net for the proposed approach.

Fig. 4 :

 4 Fig. 4: Detailed planning time of the proposed approach on the Sort and Access problems with/without AGFP-Net.

  .

	Algorithm 3 predictFeasibility
	Input: s, E, a
	1: O ← getObjectT oM ove(a)
	2: q start ← s(O)
	3: q end ← getEndP lacement(a)
	4: [p

F (P ick), p G (P ick)] ← predict(s, E, O, q start ) 5: [p F (P lace), p G (P lace)] ← predict(s, E, O, q end ) 6: p G (a) ← p G (P ick) ⊗ p G (P lace) 7: p F (a) = p F (P ick) × p F (P lace) × max(p G (a)) 8: return [p F (a), p G (a)]

TABLE I :

 I Comparison of the performance of AGFP-Net and our prior model AFP-Net.

	Model	Prediction Task	F1-Score	TPR	TNR
	AFP-Net [30]	Action Feasibility	90%	91.5% 84.3%
	AGFP-Net	Action Feasibility Grasp Feasibility	92.4% 86.6%	93.7% 93.3% 98.1% 97%

TABLE II :

 II Planning performances with and without using a feasibility heuristic, averaged over 10 runs.

				Infeasible	Total
	Domain	Method	Heuristic	Task	Planning
				Plans	Time
	Reorder	[30] Proposed	None AFP-Net None AGFP-Net	5.0 0.4 0.0 0.0	61.0 20.5 15.6 6.4
	Unpack	[30] Proposed	None AFP-Net None AGFP-Net	11.8 3.3 3.0 0.1	71.7 21.9 28.5 7.6
	Swap	[30] Proposed	None AFP-Net None AGFP-Net	30.6 3.9 17.9 0.2	125.8 22.2 117.2 7.4
		[30]		Not handled
	Sort	Proposed	None AGFP-Net	108.3 1.1	599.1 50.0
		[30]	Planning Failure
	Access	Proposed	None AGFP-Net	339.3 1.9	1500.1 95.2

LAAS-CNRS, Toulouse, France, {saitbouhsa, alami, simeon}@laas.fr

Note that this assumption is used only at the feasibility prediction level. The geometric planner used in our TAMP approach considers P ick and P lace as different tasks that might start and end at different configurations of the robot.

This separation of the views allows a more sane scene representation and prevents the neural network from confusing the views

In case none of the sampled placements lead to a geometrically feasible solution, a node can be expanded more than once in order to resample more placements.


Martin J Gander

Michal Outrata

Spectral analysis of implicit 2 stage block Runge-Kutta preconditioners

Keywords: implicit Runge-Kutta methods, two stages, preconditioners, GMRES, bounds Classication: 65F08, 65F10

Introduction

Runge-Kutta methods are a well-established family of one-step solvers for systems of ordinary dierential equations (ODEs; see [START_REF] Wanner | Solving Ordinary Dierential Equations I : Non-Sti Problems[END_REF][START_REF] Wanner | Solving Ordinary Dierential Equations II : Sti and Dierential-Algebraic Problems[END_REF] for an overview and further references). For implicit methods (IRK), their eciency relies on having an ecient solver for the so-called stage equations { in general a system of ms non-linear equations, where m is the number of scalar ODEs in the system and s is the number of stages of the Runge-Kutta method. An important application arises from the space discretization of time-dependent partial dierential equations (PDEs), resulting in a system of ODEs with very large m. If the spatial operator is linear, then the stage equations also become a system of linear algebraic equations, which are often solved by an iterative solver, e.g., a Krylov method. In [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF], the authors introduced a family of preconditioners for GMRES for the stage equations, numerically showing that these preconditioners give an outstanding performance, especially under renement of the spatial mesh, i.e., as m grows. Recently, there has been also other contributions in the direction of preconditioning the fully implicit Runge-Kutta stage equations for PDEs, see [START_REF] Southworth | Fast solution of fully implicit Runge{Kutta and discontinuous Galerkin in time for numerical PDEs, Part I: The linear setting[END_REF][START_REF] Southworth | Fast solution of fully implicit Runge{Kutta and discontinuous Galerkin in time for numerical PDEs, Part II: nonlinearities and DAEs[END_REF] but also [START_REF] Neytcheva | Numerical Solution Methods for Implicit Runge-Kutta Methods of Arbitrarily High Order[END_REF] and [START_REF] Clines | Ecient order-optimal preconditioners for implicit Runge-Kutta and Runge-Kutta-Nyström methods applicable to a large class of parabolic and hyperbolic PDEs[END_REF], introducing new ideas and testing these numerically on a variety of test problems.

We focus on the setting considered by Howle et al. in [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF] and analyze the convergence of the preconditioned GMRES method, giving a theoretical background for the performance observed for two stage methods. The general s-stage case will be treated in a follow-up manuscript. In Section 2 we summarize some preliminary knowledge and introduce the problem and the preconditioners followed by some general results for Kronecker-like matrices in Section 3. We then give the analysis of each of the block type of the preconditioners in Section 4.1 { 4.3 and also comment on the possibility of numerical optimization to improve the solution process in Section 4.4.

Model problem and preliminaries

As our model problem we consider the heat equation on the unit square and a time interval (0, T end ), i.e.,

∂ ∂t u = ∆u + f in Ω × (0, T end ), u = g on ∂Ω × (0, T end ) and u = u 0 in Ω × {0}, (1)
where ∆ is the Laplace operator, f, g, u 0 are given functions and Ω is the unit square Ω := (0, 1) × (0, 1). We discretize in space using nite dierence scheme on an equidistant grid with N + 1 rows and columns and with the mesh size h = 1/N as in Figure 1. The values at the interior grid points become unknown functions of time, which are governed by the system of ODEs,

∂ ∂t u i (t) = u i-N (t) + u i-1 (t) -4u i (t) + u i+1 (t) + u i+N (t) h 2 + b (ST ) i (t), (2)
for i = N + 1, . . . , N (N -1) -1, where b (ST) i (t) collects the known values from the source terms, given by g and f , at the given point. Combining the unknowns in each grid column into one vector denoted by u k (t), i.e.,

u k (t) := u N k+2 u N k+3 • • • u N (k+1)-1 T (t), u(t) := u 1 (t) • • • u N -1 (t) T ,
and also analogously for b k (t) and b(t), we rewrite (2) as

∂ ∂t u(t) = 1 h 2 Lu(t) + b (ST) (t), (3)
with

L =       T I I . . . . . . . . . . . . I I T       , T =       -4 1 1 . . . . . . . . . . . . 1 1 -4       , I =       1 . . . . . . 1       , ( 4
)
where L is of dimension n := (N -1) 2 , the blocks T, I are of dimension N -1. We discretize [0, T end ] with M T end + 1 equidistant time points with time step τ = T end /M T end , i.e.,

{0 = t 0 < t 1 • • • < t M T end -1 < t M T end = T end }, τ = T end M T end and t m = τ • m, m = 0, . . . , M T end . Ω = (0, 1) × (0, 1) × x0,0 × x1,0 × x2,0 × x3,0 × x0,1 × x1,1 × x2,1 × x3,1 × × × × × x4,1 × × × x4,4 h h Ω = (0, 1) × (0, 1) × u1 × u2 × u3 × u4 × u5 × u6 × u7 × u8 × × × × × u N (N -1)+1 × × × u N 2 h h
the corresponding IRK method applied to (3) at the m-th time step gives the approximation

u m ≈ u(t m ) as u m = u m-1 + τ s i=1 b i k m i , (6)
where the vectors k m 1 , . . . , k m s ∈ R n are the solutions of the linear system

      I . . . I    - τ h 2    a 1,1 L . . . a 1,s L . . . . . . . . . a s,1 L . . . a s,s L       ≡Is⊗In-τ h 2 (A⊗L)=:M k m =    1 h 2 Lu m-1 + b (ST) (t m-1 + c 1 τ ) . . . 1 h 2 Lu m-1 + b (ST) (t m-1 + c s τ )    , (7)
with

k m := k m 1 • • • k m s T ∈ R ns .
The symbol ⊗ stands for the Kronecker product (see [START_REF] Van Loan | The ubiquitous Kronecker product[END_REF] and references therein) and we would like to note here that (7) can be reformulated into a matrix equation, which is in general better suited for using a Krylov solver (see [START_REF] Palitta | Optimality properties of Galerkin and Petrov{Galerkin methods for linear matrix equations[END_REF]). Here we focus on the analysis of the results in [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF] and thus we do not address this any further but a study of the preconditioners from [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF] in the matrix equations setting seems worthwhile. Having p ≤ 2s as the order of convergence of the IRK method we assume that it is balanced with the spatial discretization error, i.e., that h 2 = C e τ p for some

C e > 0.
The problem [START_REF] Liesen | GMRES convergence analysis for a convection-diusion model problem[END_REF] with the system matrix M sparse can be very large for h (and τ ) small, suggesting an iterative solver such as GMRES, BiCG or GCR should be used which in turn usually requires a preconditioner to become truly ecient. In [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF], the authors introduce the block preconditioners

P d = I s ⊗ I n - τ h 2 diag(A) ⊗ L, P u = I s ⊗ I n - τ h 2 D A U A ⊗ L and P l = I s ⊗ I n - τ h 2 L A D A ⊗ L, (8)
where L A , D A , U A are the LDU factors of the Butcher tableau matrix A. In addition, the authors also consider the block triangular preconditioners

P GSL = I s ⊗ I n - τ h 2 A L ⊗ L and P GSU = I s ⊗ I n - τ h 2 A U ⊗ L, (9)
where GSL/GSU stands for Gauss-Seidel lower/upper, and A L,U is the lower/upper triangular part of A, i.e.,

(A L ) ij = a ij if i ≥ j 0 otherwise , (A U ) ij = a ij if i ≤ j 0 otherwise .
Some of these were considered already in [START_REF] Sta | Preconditioning of fully implicit Runge-Kutta schemes for parabolic PDEs[END_REF]. Notice that if a ii > 0 for all i = 1, . . . , s, then the preconditioners are invertible as L is symmetric negative-denite. More general conditions for non-singularity of the preconditioners can be also derived analogously to [16, Lemma 1].

Using GMRES for a linear system Cx = f with C being diagonalizable, i.e., C = SΛS -1 and Λ = diag(λ 1 , . . . , λ d ), a standard convergence bound for the residuals r reads r r 0 ≤ κ(S) min

ϕ(0)=1 deg(ϕ)≤ max 1≤i≤d |ϕ(λ i )|, (10)
where κ(S) is the 2-norm condition number of the matrix S, see, e.g., [START_REF] Liesen | Krylov Subspace Methods: Principles and Analysis[END_REF]Section 5.7.2]. We would like to highlight some aspects of this bound that is often used to study GMRES convergence behavior.

Remark 1. As indicated above, the spectral information of the system matrix in GMRES (in our case of the preconditioned system) does not generally govern the convergence (see [START_REF] Greenbaum | Matrices that generate the same Krylov residual spaces[END_REF], [START_REF] Greenbaum | Any nonincreasing convergence curve is possible for GMRES[END_REF] and [START_REF] Arioli | Krylov sequences of maximal length and convergence of GMRES[END_REF] and also [8, Chapter 2 and 5.7] and the references therein). If the system matrix is normal, i.e., it is diagonalizable with S unitary, then the spectral information is enough to evaluate the ideal GMRES bound. However, if C is non-normal, then a convincing argument needs to be put forward to validate linking spectral information with the convergence behavior of GMRES as the authors in [8, p. 303, Remark 1] point out. Moreover, particular knowledge of the interaction of S and the initial residual r 0 can lead to a qualitative and quantitative improvement on [START_REF] Neytcheva | Numerical Solution Methods for Implicit Runge-Kutta Methods of Arbitrarily High Order[END_REF], see, e.g., [START_REF] Liesen | GMRES convergence analysis for a convection-diusion model problem[END_REF]. However, studying GMRES bahavior with the bound [START_REF] Neytcheva | Numerical Solution Methods for Implicit Runge-Kutta Methods of Arbitrarily High Order[END_REF], this interaction is completely lost.

Analysis of the block preconditioners

We start by transforming the calculations into the eigenbasis of the spatial operator. Denoting the eigenpairs of L by (λ k , v k ), we organize the eigenvectors into an n-by-n matrix V and dene the block transformation matrix Q,

V := [v 1 , . . . , v n ] , and Q :=    V . . . V    ∈ R sn×sn . ( 11
)
Transforming M blockwise into the V basis gives

M := QM Q T , M =    I . . . I    - τ h 2    a 1,1 Λ . . . a 1,s Λ . . . . . . . . . a s,1 Λ . . . a s,s Λ    , (12)
with Λ = diag(λ 1 , . . . , λ n ). With the preconditioners proposed in [START_REF] Liesen | Krylov Subspace Methods: Principles and Analysis[END_REF][START_REF] Mardal | Order-optimal preconditioners for implicit Runge{Kutta schemes applied to parabolic PDEs[END_REF] we write the spectrum of the preconditioned system as

sp(M P -1 ) = sp(Q T M P -1 Q) = sp(Q T M QQ T P -1 Q) = sp M P -1 ,
where P := Q T P Q stands for one of the right-preconditioners P d,GSU,u and an analogous formulation follows also for the left-preconditioners P l,GSL . As the preconditioners are dened blockwise as scalar multiplications of L and I, their blockwise transformation into the eigenbasis of L is a straight-forward calculation -replacing L with Λ (and keeping I). Next, such matrices { block matrices with each block being a square, diagonal matrix { can be permuted into classical blockdiagonal matrices as the following lemma shows.

Lemma 1. Let C ∈ R ns×ns be a real matrix with block structure such that every block is a square diagonal matrix, i.e.,

C =    Λ 11 . . . Λ 1s . . . . . . . . . Λ s1 . . . Λ ss    , with Λ ij = diag λ (ij) 1 , . . . , λ (ij) n ∀ij. (13)
Then there exists a permutation matrix Π ∈ R ns×ns such that 11) . . . λ

Π T CΠ =    C 1 . . . C n    with C =    λ (
. . . . . . . . .

λ (s1) . . . λ (ss)    ∈ R s×s , (14)
for any = 1, . . . , n.

Hence, C is diagonalizable if and only if C is diagonalizable for all = 1, . . . , n and if

C = V -1 D V is the eigendecomposition of C with D = diag(µ (1) , . . . , µ (s) ), then sp(C) = n =1 s i=1 µ (i) , and if (µ, v) is an eigenpair of some C , then µ, Π T (v ⊗ e ) is an eigenpair of C. As a result, if C is diagonalizable with C = V -1 DV , then κ(V ) = max =1,...,s κ(V ), (15)
where κ(•) is the 2-norm condition number.

Proof. Setting E = diag(0, . . . 0, 1, 0, . . . , 0) ∈ R n×n as the matrix with the only non-zero entry being at the position ( , ) with value one, we observe that

C = n =1 C ⊗ E .
Using the Kronecker product permutation property from [START_REF] Van Loan | The ubiquitous Kronecker product[END_REF]Eqn. (1) and below], we take Π such that

Π T CΠ = n =1 E ⊗ C =    C 1 . . . C n    ,
proving the rst part of the statement. The rest follows by a direct calculation and the properties of block-diagonal matrices.

Using Lemma 1 we can analyse the eigenproperties of M P -1 directly and thereby evaluate the GMRES bound.

Remark 2. We note that an analogous lemma to Lemma 1 can also be formulated for non-normal matrices (replacing Q T by Q -1 ). Considering the Jordan canonical (or the Schur decomposition form) of C , Lemma 1 can be reformulated to obtain a block upper bidiagonal (or block uppertriangular) matrix.

To shorten the notation we set

θ k := τ h 2 λ k and Θ := τ h 2 Λ, (16)
as these quantities appear always together in the computations. By a direct calculation (see [START_REF] Outrata | Schwarz methods, Schur complements, preconditioning and numerical linear algebra[END_REF]Appendix B.8]) we get the limit behavior of θ k as τ, h → 0,

(θ n , θ 1 ) → (- 8 C e , 0), (θ -1 1 , θ -1 n ) → -∞, - C e 8 , (LIM) p=1 (θ n , θ 1 ) → (-∞, 0), (θ -1 1 , θ -1 n ) → (-∞, 0), (LIM) p>1 (17)
and continue by explicit calculations for the two stage methods.

Two stage methods

For the case s = 2 we have

P d = I -a 11 Θ 0 0 I -a 22 Θ , P GSL = I -a 11 Θ 0 -a 21 Θ I -a 22 Θ , P GSU = I -a 11 Θ -a 12 Θ 0 I -a 22 Θ , P u = I -a 11 Θ -a 12 Θ 0 I -a 22 -a 21 a 12 a 11 Θ , P l = I -a 11 Θ 0 -a 21 Θ I -a 22 -a 21 a 12 a 11 Θ , (18)
and we start with a useful lemma summarizing some direct calculations.

Lemma 2 ([11, Appendix B.8]). Let

C ∈ R 2×2 . (i) The eigenvalues µ 1,2 are given by µ 1,2 = c 11 +c 22 2 ± 1 2 √ D with D = (c 11 -c 22 ) 2 + 4c 12 c 21 . In particular, if D = 0, then C is diagonalizable.
(ii-a) If c 12 = 0 and D = 0, then C is diagonalizable and the eigenvectors v 1,2 are given by v

1,2 = 1 ṽ1,2 ṽ1,2 with ṽ1,2 = 1 α 1,2 and α i = 0 if c 11 = λ i , -c 11 -c 22 ∓ √ D 2c 12 if c 11 = λ 1,2 .
(ii-b) If c 21 = 0 and D = 0, then C is diagonalizable and the eigenvectors v 1,2 are given by v

1,2 = 1 ṽ1,2 ṽ1,2 with ṽ1,2 = α 1,2 1
and

α i = 0 if c 22 = λ i , --c 11 +c 22 ∓ √ D 2c 21 if c 22 = λ 1,2 .
(iii) If c 12 , c 21 = 0 and D = 0, then the condition number κ(V ) of the matrix of eigenvectors of C is given by κ

(V ) = ṽ1 ṽ2 + √ (1+α 1 α 2 )(1+α 1 α 2 ) ṽ1 ṽ2 - √ (1+α 1 α 2 )(1+α 1 α 2 )
, where α i is given as in (ii -a).

(iv-a) If c 21 = 0, c 12 = 0 and c 11 = c 22 , then C is diagonalizable with real eigenvalues c 11 , c 22 and eigenvectors e 1 , v 2 , and the formula from (iii) simplies to κ(V ) = ṽ2 +1 ṽ2 -1 , where v 1 is the vector with non-zero entries from (ii-a), and α 1 is its rst component. We analyze rst the block diagonal preconditioners, and then continue with the block triangular ones. The calculations below give insight into the results presented in [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF], e.g., give explicit formulas for the results in Figure 4.1 and 4.3, Table 4.3 and Table 5.1, 5.2. and 5.3 from [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF] for s = 2.

Block diagonal preconditioner

A direct calculation gives

M ( P diag ) -1 = I -a 12 Θ (I -a 22 Θ) -1 -a 21 Θ (I -a 11 Θ) -1 I , (19)
and using Lemma 1 the eigen-information of the preconditioned system can be obtained from the 2-by-2 matrices

X d k := 1 -a 12 θ k 1-a 22 θ k -a 21 θ k 1-a 11 θ k 1 .
We immediately notice that X d k is diagonalizable if and only if 1

a 12 = 0 ⇐⇒ a 21 = 0. ( 20
)
1 If a12 = a21 = 0, then A is in fact diagonal and hence M = P d , making this case uninteresting.

7

Assuming a 12 , a 21 = 0, we calculate the characteristic polynomial

p X d k (λ) of X d k , p X d k (λ) = λ 2 -2λ + 1 - a 12 a 21 θ 2 k (1 -a 11 θ k )(1 -a 22 θ k ) ,
and therefore the eigenvalues

ξ (k) 1,2 of X d k are given by ξ (k) 1,2 = 1 ± D k with D k = a 12 a 21 (|θ -1 k | + a 11 )(|θ -1 k | + a 22 ) . ( 21
)
We write ξ (k) 1,2 as functions of

|θ k | -1 , ξ (k) 1 = 1 ± φ(|θ k | -1 ) with φ(α) = a 12 a 21 (α + a 11 )(α + a 22 )
,

and α ∈ (|θ 1 | -1 , |θ n | -1
) { an interval converging towards the limit interval in [START_REF] Sta | Preconditioning of fully implicit Runge-Kutta schemes for parabolic PDEs[END_REF]. If a 11 , a 22 ≥ 0 (e.g., Gauss, Radau or Lobatto methods), then ξ (k) 1,2 lie on a line segment in C going through the point 1, which is either a part of the real axis (if sign(a 12 a 21 ) ≥ 0) or on the line 1 + βi, β ∈ R (otherwise). We have

|ξ (k) 1,2 -1| = φ(|θ -1 k |) = a 12 a 12 (|θ -1 k | + a 11 )(|θ -1 k | + a 22 )
, and hence2 the maximum of

|ξ (k) 1,2 -1| as a function of |θ -1 k |
is attained either at one of the endpoints of the interval in [START_REF] Sta | Preconditioning of fully implicit Runge-Kutta schemes for parabolic PDEs[END_REF] or at an interior stationary point. Calculating the derivative, we get

(|φ|) (α) = -sign (φ(α)) a 12 a 21 2α + a 11 + a 22 (α + a 11 ) 2 (α + a 22 ) 2 ,
and thus the only candidate for a stationary point is -(a 11 + a 22 )/2 assuming it belongs to the domain of φ (as mentioned above, this is not the case for the commonly used Gauss, Radau or Lobatto methods). Assuming it does not, e.g., because a 11 , a 22 > 0, the maximum is attained at the left endpoint of the interval in [START_REF] Sta | Preconditioning of fully implicit Runge-Kutta schemes for parabolic PDEs[END_REF], bounded from above by the value at α = 0 which gives

|ξ (k) 1,2 -1| ≤ a 12 a 21 a 11 a 22 . (22)
Thus making this quantity small will make the eigenvalues cluster tightly around 1. Notice that the bound above is suggesting to make the diagonal entries large compared to the o-diagonal ones, making the matrix diagonal in the limit (and hence making the preconditioner exact). Assuming a ij = 0 for i, j = 1, 2 and a 11 , a 22 ≥ 0 we use Lemma 2, and the condition number of the matrix of eigenvectors Recalling [START_REF] Sta | Preconditioning of fully implicit Runge-Kutta schemes for parabolic PDEs[END_REF], we see that κ S d k has no singularities in it and the limits at 0 and +∞ are also real and bounded. Therefore the conditioning of the eigenbasis of the preconditioned system will be uniformly bounded with respect to mesh renement in space and time { just as the clustering diameter.

S d k of X d k is given by κ S d k = 1 + a 21 (|θ k | -1 +a 22 ) a 12 (|θ k | -1 +a 11 ) + 1 -a 21 (|θ k | -1 +a 22 ) a 12 (|θ k | -1 +a 11 ) 1 + a 21 (|θ k | -1 +a 22 ) a 12 (|θ k | -1 +a 11 ) -1 -a 21 (|θ k | -1 +a 22 ) a 12 (|θ k | -1 +a 11 ) . ( 23
Remark 3. Some authors call such preconditioners, i.e., preconditioners such that the eigenproperties of the preconditioned system can be bounded independently of h and τ , order-optimal and the general way to show order-optimality for these kind of preconditioners has been laid out in [START_REF] Mardal | Order-optimal preconditioners for implicit Runge{Kutta schemes applied to parabolic PDEs[END_REF]. However, we want to emphasize that this independence does not mean, practically speaking, that the given preconditioner is in some sense optimal or even well-performing (and this is even more pronounced if the order-optimality is considered only wrt the spectrum, omitting the conditioning of the eigenbasis of the preconditioned system) { but in particular settings with additional reasoning this might be a useful property as we will demonstrate.

Given a Butcher tableau A, the bound (10) can be further approximated using the Joukowsky bound with the Chebyshev polynomials in the complex plane, see [14, Section 6.11 and Corollary 6.33] and also [8, Section 5.7.2] and references therein. A relatively direct and elementary calculation then allows us to evaluate the bound [START_REF] Neytcheva | Numerical Solution Methods for Implicit Runge-Kutta Methods of Arbitrarily High Order[END_REF], and in Figure 2 we show the GMRES convergence behavior together with these bounds for dierent Butcher tableaus for P d on the left (for validation of the above formulas and detailed calculations we refer the interested reader to [START_REF] Outrata | Schwarz methods, Schur complements, preconditioning and numerical linear algebra[END_REF]Section 7.3.1]). Notice that as h is xed for all methods we obtain dierent τ for each of the IRK methods and therefore the scaling factor τ /h 2 becomes h 2(1/p-1) . Hence, the dierence is not only in the choice of A but also in the interval spanned by θ k and how close (or far) these are to the limit in [START_REF] Sta | Preconditioning of fully implicit Runge-Kutta schemes for parabolic PDEs[END_REF].

We also want to address the notable stair-like behavior. Since GMRES is invariant (in exact arithmetic) to an orthogonal transformation we can focus on GMRES applied to a problem with the block-diagonal matrix

X d := diag(X d 1 , X d 2 , . . . , X d n ).
As noted in [START_REF] Faber | On Chebyshev polynomials of matrices[END_REF], the optimal polynomial3 that solves the min-max problem

r r 0 ≤ κ(S d ) min ϕ(0)=1 deg(ϕ)≤ max 1≤k≤n ϕ(X d k ) ,
equioscillates over the blocks { but only for even degrees, i.e., for = 2l for some l ∈ N. Considering our situation, the optimal polynomial is trying to minimize its value over these 2-by-2 blocks, which in our case have eigenvalues

ξ (k) 1 = ξ (k)
2 . Hence, the optimal polynomial aims to be small (in the maximum norm) on the line segment connecting 1 + i |(a 12 a 21 )/(a 11 a 22 )| and 1 -i |(a 12 a 21 )/(a 11 a 22 )|, symmetrically about 1. However, polynomials with real coecients of odd degree have at least one real root. This is why we see notably less improvement at odd iterations { the added degree of freedom for the optimization (placement of one additional root) is constrained in a quite inconvenient way for our setting. In contrast to even iterations, on odd iterations we can make the polynomial small only in the vicinity of the real axis, which after a couple of odd iterations decreases the overall maximum norm only very little, e.g., for LobattoIIIC in Figure 4 we see that initially the real root can contribute almost equally to the decrease of the maximum norm of the polynomial, but after seven or nine iterations this is no longer true. Remark 4. This reasoning is similar to the one in [8, Section 5.7.2, p.291], where the GMRES convergence bound [START_REF] Neytcheva | Numerical Solution Methods for Implicit Runge-Kutta Methods of Arbitrarily High Order[END_REF] is adapted to a symmetric indenite problem with spectrum in two intervals

I -∪ I + ≡ [-1, ν] ∪ [ν, 1]
. Then [START_REF] Neytcheva | Numerical Solution Methods for Implicit Runge-Kutta Methods of Arbitrarily High Order[END_REF] can be simplied to

r r 0 ≤ 2 κ -1 κ + 1 [ /2] , (24)
where κ = 1/ν is the condition number of the given system matrix and [ /2] is the integer part of /2. Notice that this bound gives linear convergence but only over every two iterations of GMRES { precisely the observed behavior in Figure 4.

However adapting this bound to our setting becomes very quickly very complicated and goes beyond the scope of this text. The main diculty, in our eyes, is in the treatment of the optimal polynomials with real coecients on two separated line-segments in C. For more details on the real case and the derivation of (24) we refer also to [4, Chapter 3].

Block upper-triangular preconditioner

We consider the preconditioners P GSU , P u and where necessary we join the quantities corresponding to the preconditioners by writing, e.g., X GSU,u k instead of X GSU k and X u k . By analogous calculations to Section 4.1, we obtain the formulas

X GSU k = 1 0 -a 21 θ k 1-a 11 θ k ξ GSU k and X u k = 1 0 -a 21 θ k 1-a 11 θ k ξ u k ,
hence obtaining the spectrum of the preconditioned systems as the union of (25)

Moreover, Lemma 2 shows that for

X GSU,d k = (S GSU,d k ) -1 diag(1, ξ GSU,u k )S GSU,d k we have κ(S GSU,u k ) = 1 + (β GSU,u ) 2 + 1 1 + (β GSU,u ) 2 -1 ,
where the scalars β GSU,u are given as functions of θ k , .

|β GSU | = |a 12 | ||θ k | -1 +
The conditioning of the eigenbasis S GSU,u k of X GSU,u k can be treated similarly, rst observing that

κ(S GSU,u k
) is a decreasing function of |β GSU,u | ∈ (0, +∞) with a singularity at 0. Hence, as h, τ → 0, the matrix X GSU 1 (and X u n ) becomes non-diagonalizable 4 as |θ 1 | → 0 (and |θ n | → +∞). We show the preconditioned GMRES convergence behavior and the evaluated bounds in Figure 2; for detailed calculations and validation of the above formulas we refer the interested reader to [START_REF] Outrata | Schwarz methods, Schur complements, preconditioning and numerical linear algebra[END_REF]Section 7.3.2]. Also, notice that the stair-like behavior of P d is not present as the entire spectrum is real and covers reasonably uniformly an interval on a real line and a single point 1. As the authors point out in [8, Section 5.6.2 and 5.7.2], as long as the condition number κ(S GSU,u ) is not too large, the classical linear bound based on the condition number of the preconditioned system matrix (i.e., the condition number of M P GSU,u -1 ) can be quite descriptive for the worst-case GMRES behavior. If the preconditioned system is not far away from being normal, then this explains why we only see the linear convergence { the spectrum is populating the interval considered densely enough so that the superlinear convergence argument used in exact arithmetic is not applicable (i.e., we cannot single out any outliers if there are none, see [8, Section 5.6.4 (Figure 5.7 in particular) and also Section 5.6.4]). Hence, linear convergence is, in principle, to be expected. Also, this explains the observation that the number of GMRES iterations does not grow as h → 0 as observed in [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF] (an analogous argument applies also to the block-diagonal preconditioner). Notice that the accuracy of the bound (24) is usually supported by the same kind of arguments. By analogy, this also qualitatively explains the linear in the linear-over-two-iterations convergence behavior of the block-diagonal preconditioned system.

4.3

Block lower-triangular preconditioner

The results for P GSL , P l are completely analogous to the ones from Section 4.2 and hence we just present these without much comment; more details can be found in [START_REF] Outrata | Schwarz methods, Schur complements, preconditioning and numerical linear algebra[END_REF]Section 7.3.3]. We have

X GSL k = 1 -a 12 θ k 1-a 11 θ k 0 ξ GSL k and X l k = 1 -a 12 θ k 1-a 11 θ k 0 ξ l k ,
hence obtaining the spectrum of the preconditioned systems as the union of 1, ξ GSL,l

k n k=1 with ξ GSL k = ξ GSU k = |θ| 2 det(A) + |θ|(a 11 + a 22 ) + 1 |θ| 2 a 11 a 22 + |θ|(a 11 + a 22 ) + 1 , ξ l k = ξ u k = |θ| 2 det(A) + |θ|(a 11 + a 22 ) + 1 |θ| 2 det(A) + |θ|(a 11 + a 22 -a 21 a 12 a 11 ) + 1 . (26)
Moreover, Lemma 2 shows that for

X GSL,l k = (S GSL,l k ) -1 diag(1, ξ GSL,l k )S GSL,l k we have κ(S GSL,l k ) = 1 + (β GSL,l ) 2 + 1 1 + (β GSL,l ) 2 -1 ,
where the scalars β GSL,l are given as functions of θ k ,

|β GSL | = |a 21 | ||θ k | -1 + a 22 |
, and

|β l | = |a 21 /a 11 | |θ k | det(A) a 11 + a 11 .
Notice that the results are either identical or very similar to the ones obtained with P GSU , P u and a comparison of the convergence behavior as well as of the bounds in Figure 2 and Figure 3 reects this.

Optimized Butcher tableaus

In [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF], the authors consider taking a dierent A, let us denote it by Ã, for the construction of the preconditioner, motivated by [START_REF] Sta | Preconditioning of fully implicit Runge-Kutta schemes for parabolic PDEs[END_REF]. Ã is chosen as a result of an optimization routine seeking to minimize κ( Ã-1 A) (or κ(A Ã-1 ) for a right preconditioner) subject to a particular non-zero pattern and diag( Ã) = diag(A). This optimization thus evaluates only quantities corresponding to the s-bys Butcher tableaus, which should be negligible in cost compared to the solution process of the stage equations. Numerical results in [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF]Table 4.3] show that the resulting preconditioner lowers the condition number of the preconditioned system. Having explicit formulas from Section 4.1{4.3 we can move from minimizing the quantity κ( Ã-1 A) (or κ(A Ã-1 )) to minimizing quantities in the bound [START_REF] Neytcheva | Numerical Solution Methods for Implicit Runge-Kutta Methods of Arbitrarily High Order[END_REF] so that a provable bound can be obtained based on the result.

First, adapting the eigenpair formulas in Section 4.1 { 4.3 for the case A = Ã we could optimize the eigenproperties of the preconditioned system and thus obtain a theoretical bound on the resulting preconditioner. Using Lemma 1, the matrices X k are easily accessible with any A xed.

Second, having A = Ã as in Section 4.1 { 4.3, we can still consider the same optimization problem { but naturally with extra constraints so that the resulting IRK method has some desired order of convergence and stability properties.

We numerically test these approaches, dening the objective function to be minimized as a weighted sum of the cluster diameter5 |ξ k -1| and the conditioning of the eigenbasis of the preconditioned system matrix κ(S),

f obj (A, Ã) := max k |ξ k -1| + ωκ(S),
with some positive small weight ω > 0. Minimizing f obj then aims to minimize the min-max part of (10) while keeping the eigenbasis conditioning under control. In case that the clustering is very tight we expect to get good GMRES convergence. Before showing the results we comment on the usefulness of the optimization approach for the case A = Ã.

Remark 5. In order to optimize over A = Ã we have to keep in mind that we are at the same time changing the IRK method we are using to solve the system of ODEs. This introduces some optimization constraints to ensure good order of convergence (equality constraints on A, b, c) and stability (equality and inequality constraints on A, b, c). In our experience, in order to obtain a notable improvement in minimization of f obj over choosing A = A Gauss,RadauIIA,... , we need to relax some of the qualities of the IRK methods, e.g., decrease the order of convergence by one or give up some stability. This seems quite unappealing as the maximal order of convergence p Gauss = 4 is already not too large for practical use.

In our eyes, this makes the optimization for A = Ã better suited for larger s (e.g., s ≥ 5), where decreasing the order of convergence by one or two might not be a signicant drawback as the order of convergence of the standard IRK methods is much higher; see [START_REF] Outrata | Schwarz methods, Schur complements, preconditioning and numerical linear algebra[END_REF]Section 7.5].

Note that even though Remark 5 suggests that the optimization with A = Ã is not really of interest here, it is still reasonable to start studying it for s = 2, to get some insight in the case where the formulas are explicitly available. We show the numerical results in Figure 4 and the resulting matrices in Table 1.

We see that for the setting A = Ã in the rst row of Figure 4 the optimized A are such that we converge after two to ve iterations, basically turning GMRES into a direct solver 6 . This seems natural looking at Table 1 { we observe that the Butcher tableau matrix A adapts the non-zero structure of the preconditioner, e.g., A becomes close to diagonal for P d , so that the preconditioner then becomes almost identical to M , reinforcing the point of Remark 5. Notably, this is achieved while not exploding the condition number of the eigenbasis of the preconditioned system.

For the second row, i.e., A = à and optimizing à that is used to construct the preconditioner, we see that for the block-diagonal preconditioner we can still obtain a considerable speed-up, as well as for the preconditioners P GSL,GSU . This is not the case, however, for the preconditioners P l,u . The reason is that with A = Ã, the spectrum of M (P l,u ) -1 was real but this property is lost in the general case A = Ã. Hence, even though we have tightened the clustering of the eigenvalues by optimizing à this wasn't signicant enough to o-set the introduction of the complex eigenvalues of the preconditioned system.

We would like to note that we obtained similar results to the rst row when considering p = 3 and A-stability (i.e., giving up the more restrictive property of L-stability) but for p = 3 and Lstable methods there seemed to be next to no gains from the extra optimization (see [START_REF] Outrata | Schwarz methods, Schur complements, preconditioning and numerical linear algebra[END_REF]Section 7.5]). block diagonal block lower-triangular block upper-triangular

P d P GSL P l P GSU P u A = Ã 1.8 • 10 -1 -3.8 • 10 -6 2.3 • 10 -6 3.9 • 10 -1 3.5 • 10 -1 -4.5 • 10 -9 2.5 • 10 -6 2.3 • 10 -1 3.2 • 10 -1 -3.0 • 10 -8 9.6 • 10 -1 2.7 • 10 -1 3.5 • 10 -1 4.2 • 10 -1 -5.8 • 10 -9 2.3 • 10 -1 4.9 • 10 -1 -5.2 • 10 -2 6.3 • 10 -8 1.7 • 10 -2 A = Ã 9.3 • 10 -1 4.8 • 10 -2 4.0 • 10 -1 1.3 • 10 -1 5.1 • 10 -1 -6.6 • 10 -2 8.0 • 10 -1 3.0 • 10 -1 4.1 • 10 -1 -8.4 • 10 -2 7.4 • 10 -1 2.2 • 10 -1 4.1 • 10 -1 -8.2 • 10 -2 7.0 • 10 -1 3.7 • 10 -1 4.1 • 10 -1 -8.1 • 10 -2 7.2 • 10 -1 2.3 • 10 -1
Table 1: The resulting matrices A = à (rst row) or à (second row) corresponding to the results presented in Figure 4.

Last but not least we comment on the computational costs. In order to evaluate the complete eigenproperties of the preconditioned system M (P ) -1 (or (P ) -1 M ) we rst calculate the eigendecomposition of L (which is prohibitively costly) and then we calculate eigendecompositions of n : GMRES convergence (and bounds) for block preconditioners P d (left column), P GSL , P l (middle column) and P GSU , P u (right column) for the two stage RadauIIA Butcher tableau (lled markers, same results as in Figure 2 and 3) and also the optimized setting considered above (halflled markers). In the second row the optimization process changes only the matrix à used to construct the preconditioner (and not M ) in contrast to the rst row where the optimized is used for construction of both the system matrix and the preconditioner (in this case we imposed A to be such that the resulting IRK method is at least of order two and is L-stable, the constraints can be found in [START_REF] Outrata | Schwarz methods, Schur complements, preconditioning and numerical linear algebra[END_REF]Section 7.3.4]). We used ω = 10 -5 , N = 100 and τ = h 2/p as above.

matrices X k , each of dimension s, which can be done in parallel. Even without any parallelization the cost is linear in n with the constant corresponding to s 3 . In practice, we often have (or can reasonably cheaply obtain) some estimate of the eigeninformation of L, e.g., estimates µ 1 , µ n of the true eigenvalues λ 1 , λ n or on the conditioning of V . The eigendecomposition of L can then be replaced by, e.g., considering only q \fake" λ k sampled from (µ min , µ max ) (assuming we know the spectrum of L is real). We show a numerical illustration, considering the same setting as for Figure 4 and sample only q distinct λ k (and thus θ k ) in the interval µ min := θ min , µ max := θ max . We show the evolution of the number of preconditioned GMRES with the optimized preconditioners depending on q in Figure 5. We see that for the spatial dimension n = 99 2 we already get comparable results to Figure 4 by sampling only very few \fake" θ k compared to optimizing over the actual roughly ten thousands of them. Also, the bounds stay still quite descriptive { only for q = 40 and P GSU have we found that running the optimization as a black-box does not give a useful descriptive bound 7 . This observation is key for either of the optimization approaches to be viable and further study of q for more involved settings (and its dependency on n) and/or sampling strategies for complex spectra is necessary. However, as the main interest is in s large, this will be presented in the follow-up work. 5

Generalizations and conclusion

We have shown that for two stage IRK methods the preconditioners from [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF] can be analyzed explicitly using spectral techniques and the GMRES convergence behavior can be reasonably predicted using the worst-case GMRES bound. This bound can be evaluated directly before the computations and gives a theoretical background to the results observed in [START_REF] Rana | A New Block Preconditioner for Implicit Runge-Kutta Methods for Parabolic PDE Problems[END_REF]. These results are given for a simple test problem but the analysis clearly extends to any diagonalizable spatial operator L { not just the Laplacian. The same is true also for the discretization scheme used { a nite elements scheme analysis can be done completely analogously (see [START_REF] Outrata | Schwarz methods, Schur complements, preconditioning and numerical linear algebra[END_REF]Section 7.7]). Some of the above can be generalized to s-stage IRK methods but due to the space restrictions these results will be presented in an upcoming manuscript. We would also like to mention that both in [START_REF] Neytcheva | Numerical Solution Methods for Implicit Runge-Kutta Methods of Arbitrarily High Order[END_REF] and [START_REF] Southworth | Fast solution of fully implicit Runge{Kutta and discontinuous Galerkin in time for numerical PDEs, Part I: The linear setting[END_REF], the authors multiplied M with A -1 ⊗I n from the left. The analysis above for this case is analogous and reveals that for some preconditioners we obtain much better properties and for others the performance deteriorates. This is naturally very dependent on the choice of A as well; for more details see [START_REF] Outrata | Schwarz methods, Schur complements, preconditioning and numerical linear algebra[END_REF]Section 7.6].

Figure 1 :

 1 Figure 1: Left: grid points for N + 1 = 4; right: lexicographical ordering of the unknowns for N + 1 = 4.

  (iv-b) If c 12 = 0, c 21 = 0 and c 11 = c 22 , then C is diagonalizable with real eigenvalues c 11 , c 22 and eigenvectors v 1 , e 2 , and the formula from (iii) simplies to κ(V ) = ṽ1 +1 ṽ1 -1 , where v 2 is the vector with non-zero entries from (ii-b).

Figure 2 :

 2 Figure2: We show the GMRES performance together with the evaluated bound[START_REF] Neytcheva | Numerical Solution Methods for Implicit Runge-Kutta Methods of Arbitrarily High Order[END_REF] for dierent preconditioners and Butcher tableaus. We set N = 100, τ = h 2/p (p is the order of convergence of the IRK method; p Gauss = 4, p RadauIIA = 3, p LobattoIIIC = 2) and use a random right-hand side.

Figure 3 :

 3 Figure3: We show the GMRES performance together with the evaluated bound[START_REF] Neytcheva | Numerical Solution Methods for Implicit Runge-Kutta Methods of Arbitrarily High Order[END_REF] for dierent preconditioners and Butcher tableaus. We set N = 100, τ = h 2/p (p is the order of convergence of the IRK method; p Gauss = 4, p RadauIIA = 3, p LobattoIIIC = 2) and use a random right-hand side.

Figure 4

 4 Figure4: GMRES convergence (and bounds) for block preconditioners P d (left column), P GSL , P l (middle column) and P GSU , P u (right column) for the two stage RadauIIA Butcher tableau (lled markers, same results as in Figure2 and 3) and also the optimized setting considered above (halflled markers). In the second row the optimization process changes only the matrix à used to

  Figure 5

  1, ξ GSU,u det(A) + |θ|(a 11 + a 22 ) + 1 |θ| 2 a 11 a 22 + |θ|(a 11 + a 22 ) + 1 , ξ u k = |θ| 2 det(A) + |θ|(a 11 + a 22 ) + 1 |θ| 2 det(A) + |θ|(a 11 + a 22 -a 21 a 12 a 11 ) + 1 .

			k	n k=1	with
	ξ GSU k	=	|θ| 2

  a 22 | , and |β u | = |a 12 /a 11 | |θ k | det(A) a 11 + a 11 . |θ| 2 a 12 a 21 |θ| 2 a 11 a 22 + |θ|(a 11 + a 22 ) + 1 , |θ|a 12 a 21 /a 11 |θ| 2 det(A) + |θ|(a 11 + a 22 -a 21 a 12 a 11 ) + 1 , and assuming a 11 , a 22 ≥ 0, elementary calculus reveals that the cluster diameter for P GSU is bounded from above by |a 12 a 21 /(a 11 a 22 )| { the limit as |θ k | → +∞. For P u analogous calculations show that the cluster diameter is maximized at |θ

	Further calculations give
	|ξ GSU k	-1| =
	|ξ u k -1| =
		|ξ u k -1| ≤	a 21 a 12 11 + det(A) + 2a 11 a 2 √ det(A)

k | = det(A) -1/2 (assuming det(A) ∈ (|θ n |, |θ 1 |))

with the value

We assumed at the beginning of the section that the inverse ( P d ) -1 exists and hence the denominator of φ is non-zero. Hence |φ(α)| is a smooth function.

In[START_REF] Faber | On Chebyshev polynomials of matrices[END_REF], the authors call these the Chebyshev polynomials of the given matrix { of X d in this case.

We can see this already in (25) since in the limits considered we get ξ GSU,u k = 1 thus obtaining a 2-by-2 Jordan block.

The clustering point is for all of our preconditioners naturally equal to 1.

Notice that in practice the solve accuracy also needs to be balanced with the discretization errors, hence even ve iterations to obtain a relative residual smaller than machine precision might eb considered a direct solver (i.e., converging after one or two steps) in some applications.

In fact for that particular setting the optimization routine sets o in a direction of inadmissible A (singular) and does not nd it's way back. Fine tuning the parameters of the optimization routine does, however, x this issue.

this work was partially supported by the SNF grant number 178752.

† Section de Mathématiques, Université de Gen eve; this work was partially supported by the FCS Swiss Excellence PhD Fellowship program of the Swiss Federation.