# Efficient falsification approach for autonomous vehicle validation using a
parameter optimisation technique based on Reinforcement Learning

Dhanoop Karunakaran1, Stewart Worrall1, Eduardo Nebot1 1D.Karunakaran, S.
Worrall, E. Nebot are with the Australian Centre for Field Robotics (ACFR) at
the University of Sydney (NSW, Australia). E-mails: {d.karunakaran, s.worrall,
<EMAIL_ADDRESS>
###### Abstract

The widescale deployment of Autonomous Vehicles (AV) appears to be imminent
despite many safety challenges that are yet to be resolved. It is well-known
that there are no universally agreed Verification and Validation (VV)
methodologies to guarantee absolute safety, which is crucial for the
acceptance of this technology. The uncertainties in the behaviour of the
traffic participants and of the dynamic world cause stochastic reactions in
advanced autonomous systems. In addition to that, modern autonomous vehicles
will undoubtedly include machine learning and probabilistic techniques. The
addition of ML algorithms adds significant complexity to the process for real-
world testing when compared to traditional methods. The most common approach
for evaluating system performance is based on large scale real-world data
gathering exercises (number of miles travelled), or a test matrix approach.
Recently, the research community and industry have started to move towards
scenario-based testing as an alternative option that is considered to be a
feasible methodology to realise AV use on public roads. Most research in this
area focuses on generating challenging concrete scenarios or test cases to
evaluate the system performance by looking at the frequency distribution of
extracted parameters as collected from the real-world data. These approaches
generally employ Monte-Carlo simulation and importance sampling to generate
critical cases.

This paper presents an efficient falsification method to evaluate the System
Under Test. The approach is based on a parameter optimisation problem to
search for challenging scenarios. The optimisation process aims at finding the
challenging case that has maximum return. Our method is based on Neural
Architecture Search (NAS) to search for the optimal parameter combinations.
NAS applies policy-gradient reinforcement learning algorithm to enable the
learning. The riskiness of the scenario is measured by the well established
RSS safety metric, euclidean distance, and instance of a collision. We
demonstrate that by using the proposed method, we can more efficiently search
for challenging scenarios which could cause the system to fail in order to
satisfy the safety requirements or evaluation metrics.

##  I Introduction

A guarantee of safety is essential for the acceptance of autonomous vehicles,
and society is unlikely to tolerate fatalities caused by failures of an
intelligent system[1][2][3]. The validation of Highly Automated Vehicles (HAV)
is an essential step before deploying them on public roads. Traditionally, the
test matrix approach is considered to be one of the primary validation methods
to evaluate the system. It uses test scenarios extracted from a crash database
[4]. A typical example based on this approach is the Autonomous Emergency
Braking (AEB) test protocol developed as part of the Euro New Car Assessment
Program (Euro-NCAP)[5]. Even though this approach is repeatable and reliable,
these scenarios are fixed and predefined[4]. As a result, the system is likely
to perform well in these tests, but this is not a true test of real-world
performance as the real-world is non-deterministic and challenging[4].

Figure 1: Carla simulation of the experiment: a pedestrian crossing the street
while the vehicle is moving.

One possible method to evaluate the system is to use a Naturalistic Field
Operation Test (N-FOT) that involves driving on public roads to validate the
system statistically. The motivation of this distance-based validation is to
uncover rare events that are not captured during the validation process. Such
occurrences are by definition rare events and are not uniformly distributed,
meaning that the AV will need to drive a large number of miles[2][6][7]. In
addition to that, the collected evidence will be undermined with every
software change [2][6][7]. Thus, this approach is considered to be time-
consuming and expensive[2][8].

Recently, the research community and automotive industry have broadly started
to support scenario-based testing as an alternative strategy
[4][9][10][11][12]. This has been considered as a feasible solution to
realising HAVs on public roads. Most of the driving data collected in N-FOT is
considered non-critical, thus not meaningful for validation purposes. The
scenario-based approach for validating HAVs is promising as it focuses on
critical scenarios without the need of driving millions of kilometres [6][9].
Scenario-based testing has been implemented in real-world and simulated
settings. The detailed definition of the scenario and different
classifications are explained in the methodology section. The aim of this
paper is to find critical, concrete scenarios as defined by this
classification. We use the term ‘scenario’ to refer to a concrete scenario,
and the terms ‘test case’, ‘challenging scenario‘, ‘challenging
case‘,‘critical scenario’ and ‘edge cases’ are all used interchangeably.

From previous works [4][10][11] in the area of scenario-based testing, the
frequency distributions of parameters that are considered influential are
extracted from N-FOT. From these, concrete scenarios are generated based on
observation drawn from the distribution. The system gets evaluated by running
a Monte-Carlo simulation using generated test cases and importance sampling
using the test cases that were identified as critical in the Monte-Carlo run.

In our work, we considered that scenarios are constructed using different
parameters [4][9][10] and that critical scenarios are considered as a
combination of parameter values found to maximise the safety risks. The
parameter optimisation technique is defined as the approach to find the the
combination of parameter values that create challenging cases using an
optimisation engine. The main contribution of this paper is to propose a
innovative falsification method by considering the generation of challenging
scenarios as a parameter optimisation problem. The falsification approach
focuses on finding challenging cases in order to evaluate the System Under
Test(SUT)[13]. The whole approach is considered as an optimisation guided
falsification method, and the optimisation process looks to find the critical
scenario that has the maximum return. The proposed method is based on Neural
Architecture Search(NAS)[14], which is the search for the best neural network
architecture for a specific task. The NAS architecture is based on a policy
gradient reinforcement learning (RL) algorithm. We have designed a reward
function for the RL process in such a way that a positive reward for the
combination of considered parameter values that can create challenging
scenarios. In this experiment, we use the well established Responsibility-
Sensitive Safety (RSS) metric[2][15], the euclidean distance of pedestrian
w.r.t ego-vehicle, and whether there was a collision as the metrics. The
pedestrian crossing illustrated in figure 1 is used as the experimental
context in the simulator to test our method.

We can summarise the contribution of this papers as follows:

  * •

We propose an innovative falsification method to search for critical scenarios
using a parameter optimisation technique based on NAS as the optimisation
engine. The learning of the optimisation engine searches for critical
scenarios that has maximize the measure of safety risk. Due to the RL
strategy, the proposed method has better performance in searching for
challenging scenarios compared to random search as used in many recent
scenarios-based testing approaches[4][10][11][16], or brute-force search.

  * •

We demonstrate that the proposed approach is scalable as it can handle high-
dimension problems compared to our previous work[17]. Also, it can incorporate
more parameters in a computationally efficient manner as opposed to the prior
work.

  * •

The proposed method can create more realistic scenarios compared to similar
work in this area[17][18].

##  II Related work

ISO 26262 is a widely accepted functional safety standard in the automotive
industry, and the focus of this standard is on conventional vehicle
systems[19]. This standard looks at preventing malfunctions(faults and
failures) which are mainly related to internal causes[20]. However, the most
functionally safe vehicle can still have an accident due to the performance
limitations of the system caused by unintended behaviour when the system
interacts with the external environment[2]. The recently released ISO 21448
addresses this problem by aiming to reduce the number of unintended behaviours
that a system can produce[21]. The standard suggests that safety of intended
functionality is ensured by explicitly evaluating the known-unsafe scenario
and evaluating unknown-unsafe scenarios using industry best practice,
systematic analysis, or dedicated experiments[21]. The goal of the standard is
to reduce the known unintended behaviours (known unsafe) and the unknown
potential behaviour (unknown-unsafe) to an acceptable level of residual
risk[22]. These existing standards are essential to provide best practices for
developing safety systems, but they are not comprehensive enough to cover the
safety of Highly Automated Vehicles[23]. These standards assume the ultimate
guarantee of safety is the responsibility of the human driver. This is a
fundamental problem for HAV operating on the road without a human driver.
There are many technologies used in HAV that are inherently incompatible with
existing automobile safety standards[24][22][23].

A common strategy is to use a test matrix for system validation[4]. In this
approach, the system is evaluated using pre-defined test scenarios extracted
from the crash database[4]. One such example is the AEB Autonomous Emergency
Braking test protocol developed as part of the Euro New Car Assessment Program
(Euro-NCAP)[4][25]. Even though this approach is repeatable and reliable,
these scenarios are fixed and predefined[26]. So, the system is likely to
perform well in these tests, but the real-life scenarios are more diverse and
challenging[4].

Evaluating the performance of the system in a real-world environment is an
essential aspect of robotic system design. Several Naturalistic Field
Operation Tests(N-FOT) have been conducted around the world to evaluate the
performance of automated vehicles using real-world data[4]. In a real-world
setting, exposure to rare events is uncommon[2][4]. So this naturalistic data
based approach or distance based approach needs to collect a large number of
miles of data to make statistically significant observations[2]. Although
real-world testing is an important step in the development process, these
approaches are time-consuming and expensive[2].

Recently, a significant body of research focused on scenario-based testing has
been presented as alternative option to validate the SUT compared to
naturalistic data based approach[6][8][27][28]. This type of approach is based
on the assumption that most scenarios in the real world are non-critical and
there is a need to generate critical scenarios intentionally[28][29]. Hence,
it can reduce the required validation effort as a scenario based approach is
focused directly on the relevant critical scenarios[28][29]. In [10],
scenarios are parameterised and stored in a database. Then, Monte-Carlo
simulations are run to generate test cases from the parameterised scenarios.
Similarly, in [4], frequency distributions are extracted from the real-world
data, and Monte-Carlo simulations are used to generate test cases by sampling
from these distributions. The ENABLE-S3 project[11] utilises the work from
[30] for the validation process. A method to extract the scenarios from real-
world data is demonstrated by parameterising these scenarios. Then test cases
can be derived from the frequency distribution generated during the
parameterisation. Similarly, in Pegasus[12], data from different sources are
used to cluster similar scenarios together to generate logical scenarios. Test
cases (concrete scenarios) are created based on the logical scenarios. One of
the challenges of scenario-based testing is parameter space explosion [9]. The
scenario-based testing approach has to consider a high number of influential
parameters. This leads to a high dimensional parameter space for validation,
resulting in parameter space explosion. In other words, this approach requires
a large number of concrete scenarios per logical scenario. In [9], a
functional decomposition method is proposed to reduce the effect of the
parameter space explosion.

Waymo recreates the scenarios from the real-world data in simulation and
generates thousand of variations of the recreated scenarios to validate the
system[31]. They heavily rely on simulation to identify the critical events
that are not encountered during real-world testing by creating variations of
each scenario. This approach enables them to release software updates quickly
rather than requiring millions of additional miles for validation.

[18] presents an interesting approach for validating the autonomous vehicles
at a pedestrian crossing. They have proposed a method to find for most likely
failure scenarios using a process called Adaptive Stress Testing. The primary
strategy of the approach is to search for the most critical scenarios.
Nevertheless, as there are no constraints on the variation of parameters this
approach could generate non realistic scenarios, with non-smooth trajectories
from the simulated participants.

In our work, an efficient falsification approach to evaluate the SUT is
proposed. The method searches for critical scenarios by considering the
generation of scenarios as a parameter optimisation problem. The scenarios are
constructed by varying the parameters and applying a loss function to optimise
towards challenging scenarios. The approach is an optimisation-guided
falsification method, and the optimisation process tends the scenarios towards
the challenging cases with the maximum cost. In this work, we use NAS as an
optimisation engine to create challenging scenarios. Figure 5 is the high-
level overview of the proposed method where the search space represents the
selected parameters and the search strategy is the optimisation approach. A
new scenario is generated in the simulation using a mixture of the values for
each parameter. The reward function is used to determine the level of safety
risk in the generated scenario created from a new set of generated parameters.
The reward function used is based on a combination of the Responsibility
Sensitive Safety (RSS) metric[2], the euclidean distance of pedestrian w.r.t
ego-vehicle, and the instance of a collision.

##  III Methodologies

###  III-A Scenario definition and classification

Figure 2: Example of three main classification of scenarios.

In [6], the scenario is defined as sequence of events. An event is defined as
actions and manoeuvres performed by the host vehicle and other traffic
participants within a scenario. Therefore, a scenario can consist of one or
more events that are typically performed by the host vehicle (turn, lane
change, brake) and the other interacting traffic participants (cutting-in,
cutting-out, braking). In addition, a scenario can contain information about,
for example, the weather conditions (sunny, rainy), type of the road, state of
the driver and the static environment. A scenario is classified into three
categories as shown in figure 2: functional scenario, logic scenario, and
concrete scenario[8]. The functional scenario is concerned about a high-level
aspect of the scenario. For instance, types of road, list of traffic
participants etc. are considered to be components of a functional scenario. We
define a logic scenario as a functional scenario when a parameter range is
specified. This could be, for example, ego-vehicle speed between 50 kph to 100
kph and pedestrian speed between 1 to 2 meters per second. The concrete
scenario is defined when we take concrete values for each parameter defined in
the logic scenario. For instance, 60 kph for ego-vehicle speed and speed of
1.1 m/s for pedestrians. Finally, adding a pass or fail criteria to concrete
scenarios is called a test case.

###  III-B Markov Decision Processes and Reinforcement Learning

A Markov Decision Process (MDP) is a Markov Reward Process (MRP) with decision
capabilities. It represents an environment in which all of the states hold the
Markov property, meaning that the future is independent of the past given the
present[32]. The solution of an MDP is an optimal policy that evaluates the
best action to choose from each state[32]. In MDP, at each time step $t$, the
agent takes an action $a_{t}$ in the environment. Then, it receives the
observation from the environment in the form of a new state $s_{t+1}$ given
the current state $s_{t}$. State $s_{t}$ is the state representation of the
agent constructed by either directly or indirectly observing the environment
state[32]. The reward $r_{t}$ is a scalar quantity to give feedback to the
agent about how appropriate (given the cost function) the action $a_{t}$ was
based on the state $s_{t}$ from the environment. The policy is a mapping from
states to actions[32]. In other words, policy determines how the agent behaves
from a specific state. A value function or state-value is the expected total
reward, starting from state $s$ and acting according to the current
policy[32]. An action-value is the expected return for an agent starting from
state $s$ and taking an arbitrary action before continuing on according to the
current policy[32].

A Reinforcement Learning(RL) algorithm finds an optimal policy that maximizes
the reward by interacting with the environment that is modelled as an MDP
where no prior knowledge about the MDP is available[32][33].

###  III-C REINFORCE

REINFORCE implements Monte-Carlo sampling of a policy gradient method[34].
That means the RL agent samples from the starting state to generate the goal
state directly from the environment, in comparison to other techniques that
use bootstrapping such as Temporal Difference Learning and Dynamic
programming. A policy gradient algorithm is a policy iteration approach where
the policy is directly manipulated to reach the optimal policy that maximises
the expected return. This type of algorithms is known as model-free RL. The
model-free indicates that there is no prior knowledge of the model of the
environment. In other words, we do not need to explicitly know the environment
dynamics or transition probability. Pseudocode for the REINFORCE as
follows[34]:

  1. 1.

Sample N trajectories using policy $\pi_{\theta}$

  2. 2.

Evaluate the gradient of the objective function J using the below expression:

$\nabla J(\theta)\approx\frac{1}{N}\sum_{\tau\in
N}\sum_{t=0}^{T-1}\nabla_{\theta}log\pi_{\theta}(a_{t},s_{t})R(\tau)$ (1)

Where N is the number of trajectories and R is the total return of a
trajectory

  3. 3.

Update the policy parameters

$\theta=\theta+\alpha\nabla J(\theta)$ (2)

  4. 4.

Repeat 1 to 3 until we find the optimal policy $\pi_{\theta}$.

The cross-entropy loss in the Neural Network(NN) can be used as a grad-log
policy for the policy gradient expression. It opens up the use of NN as the
functional approximator to represent the stochastic policy $\pi_{\theta}$.

###  III-D Neural Architecture Search (NAS)

Figure 3: Overview of Neural Architecture Search(NAS)[14]

NAS is a hyperparameter optimisation approach to design a neural network for a
specific task[14]. As shown in figure 3, the controller generates the
architectural hyperparameters of the neural networks. The controller is
implemented as a recurrent neural network. The whole problem is set up as a
reinforcement learning problem to optimise the controller to generate the
combination architectural hyperparameters to design the NN that maximises the
reward. In [34], a policy gradient RL method called REINFORCE is presented to
solve the optimisation process.

###  III-E Responsibility Sensitive Safety (RSS)

The ultimate goal of RSS as a formal method is to guarantee that an agent will
not cause an accident, rather than to guarantee that an agent will not be
involved in an accident[2]. RSS has formalized the set of rules using four
realms: Safe Distance, Dangerous Situation, Proper Response, and
Responsibility[15]. A safe distance is calculated longitudinally and laterally
based on the formula provided by the method. This method considers the worst-
case scenario which eliminates the need for estimating road user intentions.
In this paper, we focus on the longitudinal safe distance. By definition, the
longitudinal safe distance is the minimum distance required for the ego
vehicle to stop in time if a vehicle or object in front brakes abruptly.

$d_{min}=\left[v_{r}\rho+\frac{1}{2}a_{max,a}\rho^{2}+\frac{(v_{r}+\rho
a_{max,a})^{2}}{2a_{min,b}}-\frac{v_{f}^{2}}{2a_{max,b}}\right]$ (3)

where $\ d_{min}$ represents the longitudinal safe distance, $\ v_{r}$ and $\
v_{f}$ are the velocity of the agent vehicle and front vehicle, respectively.
$\ a_{min,b}$ is the minimum reasonable braking force of the agent vehicle, $\
a_{max,b}$ is the maximum braking force of the front vehicle. In terms of
acceleration, $\ a_{max,a}$ is the maximum acceleration of the front vehicle.
$\rho$ is the agent vehicle response time.

Based on the safe distance calculation, RSS can determine whether or not the
ego vehicle is in a dangerous situation. As per the RSS rule, the proper
response should be evaluated and executed once it is determined that the
vehicle is in a dangerous situation. If the ego-vehicle that follows the RSS
involved in an accident, RSS can give an assurance that the responsibility of
the accident does not belong to the ego-vehicle.

###  III-F Carla Simulator

Figure 4: Carla simulator

Carla is an open source simulator for the development, training, and
validation of autonomous urban driving systems[35]. The simulator consists of
realistic participants and objects to represent the scenarios which are
similar to the real world. The provided ROS bridge allows the extraction of
data from the simulation directly in ROS format. As our real world research
vehicle platform uses ROS, any algorithms from the real vehicle can be
connected through the ROS bridge to operate in Carla. This allows us to test
any sub-system in simulation that is used in the real vehicle. The ROS bridge
has been modified for the purpose of building the proposed method.

###  III-G Proposed method

Figure 5: This figure depicts a high-level overview of the proposed
falsification method. The search space represents the range of possible
parameter values. The search strategy for the optimisation process is a policy
gradient RL called REINFORCE to generate scenario in each episode and to learn
towards the challenging case that has maximum return.

The proposed method is a falsification method based on a parameter
optimisation technique to search for the challenging scenarios of the SUT. As
shown in Figure 5, the parameter values are defined in the search space. Our
search strategy is to optimise the search to find the combination of parameter
values that leads to challenging cases. This is done using a policy-gradient
reinforcement learning algorithm. The whole approach comes under optimisation
guided falsification, and the learning process is towards the challenging case
that has the maximum return.

Figure 6: This figure illustrates the detailed architecture of the proposed
method—the overall process considered as Policy gradient RL. The purpose is to
generate challenging scenarios by optimising the controller to produce actions
that leads to a cases with the maximum return, indicating scenarios with
higher safety risk. The reward is used to update the weights of the
controller. Each action is a combination of parameter values that are taken
from the search space.

As shown in figure 6, the architecture of the proposed method is a modified
version of the original Neural Architecture Search(NAS) [14]. The overall idea
is to generate critical scenarios by optimising the controller using the
policy gradient algorithm. The action produced by the controller in each
episode is the combination of the parameter values. The end result is that the
controller will optimise the parameters to tend towards the scenarios with the
highest risk by maximising the cost function. During the initial stages, the
controller has not yet learned suitable parameter combinations, so the actions
initiated by the controller may not create challenging cases. At later stages,
the controller has already been optimised to produce the right mixture of
values as an action. In each episode, the controller outputs the actions as a
set of parameter values. Based on these values, a scenario is constructed in
the Carla simulator. The simulation is carried out to compute the reward,
determining the level of safety risk generated by the action. After the
episode ends, the action and reward is stored. The controller is then updated
in the direction of maximum rewards after N episodes. More specifically, we
are using a Recurrent Neural Network(RNN) as the controller, so the weight of
NN gets updated based on the gradient of the objective function.

The proposed method is an efficient falsification approach to generate
critical scenarios for evaluating the SUT as compared to many random-search
based scenario-based testing and brute-force testing to finding these
scenarios. In this approach, we consider the value range for each parameter
instead of a constant value. The optimisation process is over all of these
parameters. As the proposed method handles this relatively high-dimensional
problem, it is more scalable than our previous work. In addition, the method
is computationally efficient enabling it to incorporate more parameters. The
proposed system can create more realistic scenarios by considering the
existing parameter frequency distribution from previous works to define the
parameter range; for instance, pedestrian velocity and acceleration values can
be obtained from frequency distributions described in earlier works[36][37].
Incorporating these constraints allows for the creation of smooth
trajectories, which results in more realistic scenarios.

##  IV Experiment

Figure 7: Pedestrian crossing intersection

In this work, we conducted experiments at a pedestrian crossing using the
Carla simulator, as shown in Figure 7. We have implemented the proposed
method, as depicted in Figure 6. The goal of this method is to evaluate the
SUT using a falsification approach by searching for challenging cases using a
parameter optimisation technique. The overall process is to teach the
controller to predict the suitable action that can create critical scenarios
by biasing towards the scenario that has the maximum return. At each episode,
the controller predicts an action that is a combination of parameter values
from search space. Based on the list of parameters in the action, a concrete
scenario or test case is constructed. Once the scenario is completed, a reward
function will be used to validate the quality of the action from the
controller. The reward function is then used to update the controller.

We have selected a relatively simple task in order to demonstrate the method,
with more complexity to be added in future work. We expect that this approach
can be used to determine the critical combinations of parameter values for
autonomous vehicle sub-systems such as collision avoidance systems, to create
challenging cases when working under a particular ODD. At this stage, we are
not considering the pedestrian intention and uncertainties from the ego
vehicle’s sensor measurement.

The main components of the experiment are explained below, before discussing
the implementation of the proposed method.

System Under Test (SUT): We are using a Collision Avoidance System (CAS)
provided by the Carla simulator to test our method. As depicted in the Figure
8, the maximum detection range of the system is 10 meters. The systems basic
functionality is to control the speed of the vehicle to avoid a collision when
an object is detected. Importantly, we are treating this system as a blackbox
system.

Figure 8: Collision Avoidance System (CAS)’s detection range

Search space: The search space contains the parameters and their frequency
distribution to create the scenarios. We have considered 5 parameters in this
experiment: Ego-vehicle longitudinal start position (ego-long-pos), pedestrian
acceleration (ped-accel), pedestrian velocity (ped-vel), pedestrian
longitudinal start position (ped-long-pos), and weather. Carla has 15 weather
settings, and values from the weather parameter indicates each of these
settings. The list of parameters and their distribution is shown in the Table
I. No. of samples is an indicator of how many values are sampled from the
distribution for each parameter.

TABLE I: List of parameter in the search space Parameter | Distribution | No. of samples
---|---|---
ego-long-pos | $\mathcal{U}(1,10)m$ | 10
ped-accel | $\mathcal{U}(0,0.1)m/s^{2}$ | 10
ped-vel | $\mathcal{N}(1.46,\,0.24)m/s$ | 25
ped-long-pos | $\mathcal{U}(3,4.5)m$ | 4
weather | $\mathcal{U}(0,14)$ | 10

The samples taken from the distributions in this experiment as follows:

ego-long-pos:$[1,2,3,4,5,6,7,8,9,10]\\\ $

ped-accel: $[0.046,0.051,0.097,0.075,0.099,0.076,0.065,\\\
0.007,0.007,0.013]\\\ $

ped-vel: $[1.803,1.178,1.139,1.476,1.205,1.725,1.142,\\\
1.516,1.201,1.614,1.808,1.303,1.565,1.416,1.247,1.355,\\\
1.755,1.237,1.196,1.303,1.552,1.344,0.937,1.108,0.976]\\\ $

ped-long-pos:$[3,3.5,4,4.5]\\\ $

weather: $[4,1,7,8,6,5,8,12,9,2]\\\ $

State space and Action space: The method is based on NAS[14] and both state
space and action space represents the same item. The controller accept state
$s_{t}$ and reward based on the appropriateness of the action $a_{t-1}$ and
generates the next action $a_{t}$. State represents the previous action
$a_{t-1}$.

$Parameters,\ P=(ego-long-pos,\ ped-accel,\\\ ped-vel,\ ped-long-pos,\
weather)$ (4)

$a_{t}=P,\ s_{t}=a_{t-1}$ (5)

Reward function: This function determines the safety risk for an action given
the state. At each time steps$\ t$, the minimum safe distance is measured
using the RSS formula. This metric compares with the Euclidean distance to the
pedestrian. The expected behaviour of the SUT is to reduce the speed of the
ego vehicle when an object is detected. If the Euclidean distance is less than
the safe distance measured by RSS, it is considered as an improper response
for the SUT. We classify this as a high-risk timestep, $highrisk_{ts}^{RSS}$.
At the end of the episode, a total number of high-risk timesteps get
normalised between 0.01 and -0.01 using the equation 6. It is essential to
normalise the reward to a small value range as a large range can affect the
ability of the controller to learn.

$new_{x}=(b-a)*\frac{x-min}{max-min}+a$ (6)

Where $x$ is the normalised total number of improper responses by the SUT,
minimum (min) value for the no. of high-risk timesteps and maximum(max) value
is the total no. of timesteps in an episode. The new minimum value $a$ is
-0.01 and the new maximum value $b$ is 0.01

A higher number of improper response by the SUT gives a high reward as the
objective is to optimise the parameters towards the challenging cases. If
there is a low number of high risk responses, this indicates that the system
is less likely to have safety risks. If the combination of parameters lead to
a collision, we reward an extra 0.25 to the total reward.

$new_{x}=\frac{x-min}{max-min}$ (7)

The third type of reward is based on the euclidean distance between pedestrian
and ego-vehicle. When the distance gets smaller, it is more likely to be in a
high risk situation due to RSS, or the potential of collision. At the end of
the episode, the final distance is calculated and normalised between 0 and 1
using the equation 7. To provide a higher reward for a smaller distance, we
reverse the normalised distance by subtracting the normalised distance from 1.
Then, using the equation 6, we further normalise the distance between -0.01
and 0.01 for the reward. The final reward gets computed by adding the output
of these three reward functions.

$R=\Bigg{\\{}\begin{matrix}-0.01\ to\ 0.01&\ highrisk_{ts}^{RSS}\\\ -0.01\ to\
0.01&\ euclidean\ distance\\\ 0.25&\ collision\\\ \end{matrix}$ (8)

Episode: Typically, the episode is a sequence of states and actions. In this
architecture, we have only one state and action per episode. A scenario is
constructed from the parameter values in each action. The end of a scenario is
determined either by time elapsed, distance travelled by ego vehicle, or in
the event of a collision. Once the scenario ends, it is considered the end of
the episode and the reward is computed. Each scenario is set to run up to a
maximum of 4000 episodes. We set this number arbitrarily to demonstrate the
proposed method.

Evaluation metric: The evaluation metric determines whether a scenario can be
considered as challenging or non-challenging. The$\ distance_{Euclidean}$ is
the euclidean distance between pedestrian and ego vehicle, and $\
distance_{RSS}$ is the minimum safe distance measured by RSS. We calculate
both values at each timestep of an episode. If $\ distance_{Euclidean}$ < $\
distance_{RSS}$, then it is marked as high-risk timestep. The high-risk
timestep is represented as $\ highrisk_{ts}^{RSS}$. We consider the overall
scenario as a critical scenario if the number of high-risk timesteps is
greater than 50$\%$ of the total timesteps. Otherwise, it will be designated
as a non-challenging scenario. We arbitrarily chose 50$\%$ as a threshold
limit to demonstrate the proposed method, further study into the appropriate
values for this number are part of future work. Another criterion to determine
the challenging scenario is based on the collision. Even if the total
$highrisk_{ts}^{RSS}$ is less than the threshold, a scenario is set to be a
challenging scenario if there is a collision.

$Challenging_{scenario}=\Bigg{\\{}\begin{matrix}\frac{1}{N}\sum_{i=1}^{N}highrisk_{ts}^{RSS}*100\
>=\\\ 50\%\\\ or\\\ Collision\\\ \end{matrix}$ (9)

Controller and training: The learning aspect of the controller as shown in the
figure 6 governs by the policy gradient RL algorithm called REINFORCE. The
goal of the algorithm is to optimise the controller by generating the action
that leads to the most critical events. In RL terms, the controller learns the
policy that maximises the return and the policy is represented in this
experiment by the controller. In REINFORCE, a Neural Network(NN) can be used
as a Universal Function Approximator to represent the policy. We use a
Recurrent Neural Network(RNN) as functional approximator to represent the
policy in this experiment.

Figure 9: Architecture of the controller

As shown in the figure 9, there is a dense layer for each output. In this
experiment, five dense layers represent the five parameters. The number of
output units in each dense layer varies based on the number of samples for
each parameter. Each dense layer has a separate Categorical Cross-entropy loss
function. We choose the Adam optimiser to train the network.

The training process of the NN in REINFORCE is slightly different compared to
a standard NN. Usually, in every epoch the gradient of the loss function is
computed. Then weights get updated according to the direction of the gradient.
In a supervised learning setting, we have the ground truth, and the loss
function can compute by how much the prediction has deviated. In RL, there is
no ground truth, and a reward is an indicator of how well the action generated
the desired response. In this case, it is essential to update the weights in
the direction of the reward.

Instead of updating weights in every epoch, the RL approach waits for N
episodes. Once N trajectory samples are collected (a trajectory contains
state, action, and reward), the gradient of the loss functions is computed
based on the equation 1 of the REINFORCE algorithm. The value of N is set to
25 episodes in this work. In equation 1, the reward is multiplied with
gradient log of the NN output. The NN represents the policy and the gradient
log of the output is the gradient of the standard cross-entropy loss function.
Based on the equation 2, we can update the weights of the NN, $\theta$ with
the gradient expression computed in the equation 1.

In RL, the balance between exploitation vs exploration is an important aspect.
Exploration at the early stages enables the model not to get stuck in a local
optima. Initially, an agent is weighted towards exploration compared to
exploitation as determined by the $\ \epsilon$ value. The initial value of $\
\epsilon$ is set to 1, decaying at each time step by multiplying by 0.995
until $\ \epsilon_{min}$=0.01. During this time, the actions selected for the
agent will be random.

Implementation: The proposed method evaluates the SUT by searching for the
challenging scenarios where the system fails to satisfy the safety
requirements, or the evaluation metric threshold. The search process is
completed by optimising the controller towards the critical case that has
maximum return. Firstly, the search space is created by sampling from the
parameter frequency distribution as shown in the table I. Once the search
space has been defined, the iterative process of learning begins.

At each episode, the controller generates the action given a state. During the
exploration stage, the controller has not yet learned anything, so the output
is towards a randomly generated action. Later, the controller learns the
optimal policy to create an action that leads to a challenging scenario for
the SUT.

The generated action contains five parameter values. The method uses these
values to constructs a scenario in the Carla simulator. In this work, we are
focusing on developing critical cases for the SUT in a pedestrian crossing.
The realistic scenarios get constructed as the velocity and acceleration are
sampled from the frequency distributions based on previous works. More
realistic scenarios are constructed by incorporating the acceleration and
velocity constraints which limits the change in pedestrian speed over
sequences of timesteps.

The reward is based on three criteria: RSS failure, euclidean distance, and
collision. At the end of each episode, the final reward is computed by adding
the individual rewards from these three criteria. A high reward is applied to
the action that leads to failure (high risk) and low reward for the successful
(low risk) scenarios.

As the state is the previous action, we set the last action as the current
state for the controller. The simulator gets reset to generate the next
scenario based on the following action. Action and reward are then stored for
25 episodes. The controller represents the policy, and after every 25
episodes, the policy gets updated using the REINFORCE algorithm to updating
the network weights,$\theta$.

In this work, we have set the controller to operate for 4000 episodes to learn
the optimal policy. Depending on the number of parameters and the values, we
may need to adjust this limit.

##  V Results

The results of our approach are demonstrated using several metrics that show
how the potentially unsafe challenging cases can be learnt using our method.
The goal of this method is to evaluate SUT by searching for challenging
scenarios in which the SUT fails to satisfy the evaluation metrics. As shown
Figure 10, the proposed method begins to converge after the 1500th episode.
The average reward remains almost constant from 2000th episodes onwards. It
means that the controller has learned to predict the correct action that
creates a challenging case with maximised return. In the brute-force search or
enumeration technique, we need to go through $12^{5}=248832$ concrete
scenarios where 12 is the average number of values in a parameter, and 5 is
the total number of parameters. Similarly, in a random search based method, we
do not know how many concrete scenarios needed to be considered. This result
indicates that our method can evaluate the SUT by generating the challenging
scenarios efficiently compared to the random search based scenario-based
testing and brute force search.

Figure 10: Metric used to measure the performance of the method. The plot
shows that after the 1500th episode, the approach learned the action that can
obtain a maximum return. It indicates that the controller has started to
predict the right action that maximises the return just over 1500th episode.
From 2000th episode onwards, the reward is almost constant.

(a) This diagram depicts the challenging scenario that has a maximum
return—this particular scenario ended up in a collision. The figure shows the
speed of the pedestrian and ego-vehicle, and the parameter values that create
the challenging case. The figure illustrates that ego-vehicle reduced the
speed after seeing the pedestrian, but accelerated again which ended up
causing a collision.

(b) This figure depicts the number of high-risk timesteps in a challenging
scenario. It indicates that continuous high-risk timesteps can lead to
collision or challenging scenarios based on the evaluation metric.

Figure 11: This plot illustrates the challenging scenario that finishes in a
crash.

The definition of a challenging case is based on the evaluation metric. A
typical challenging scenario that leads to the collision is depicted in figure
11(a). The diagram illustrates the pedestrian speed(red) and ego vehicle
speed(green) at each time step of the scenario. Also, all of the parameter
values that are required to create the scenario are shown at the top of the
figure. This particular scenario is constructed using the action from the
controller. The controller converged to the optimal policy towards the end of
the learning process. In this scenario, collision occurred as a result of the
speed of the pedestrian being around 1.237 m/s, and acceleration of
0.007$m/s^{2}$ with the starting position of the ego vehicle at 9m and
pedestrian starting position at 3.5m. The contributing factor for this
scenario is that the pedestrian managed to reach the centre of the crosswalk
while ego-vehicle got closer to the pedestrian without reducing the speed.
After reducing the speed, the ego-vehicle accelerated again causing the
collision with the pedestrian.

Figure 12: The figure depicts the comparison between a challenging scenario
and non-challenging scenario based on high-risk and low-risk timesteps. It is
evident that the challenging scenario has more high-risk timesteps than the
non-challenging scenario.

We observed a collision after a number of continuous highrisk timesteps. Based
on 11(b), continuous high-risk timesteps occurred from the 39th timestep until
the end of the episode. As shown in the figure 12, we can see the difference
between the number of high-risk timesteps and low-risk timesteps for a
challenging scenario and non-challenging scenario.

(a) The figure explains the typical non-challenging scenario in the
experiment. Pedestrian reached at the centre of the crossing long before the
ego-vehicle reached the crossing.

(b) The diagram shows that there is a continuous high-risk timesteps for a few
timesteps, but it’s not consistent.

Figure 13: This figure depicts the example of typical non-challenging scenario
based on the evaluation metric.

Figure 13 illustrates the typical example of non-challenging scenario based on
the evaluation metric. In this case, a pedestrian walked almost at the same
speed as in the challenging scenario with comparatively less acceleration. The
contributing factor is that the acceleration and speed of the pedestrian are
higher than the challenging scenario. As a result, the pedestrian reaches the
centre long before the vehicle arrives in the crossing area. There is a
smaller number of high-risk time steps and no collision, which defines this as
a successful (i.e. low risk) scenario.

Figure 14: The figure depicts the approach taken in our previous work[17]. The
changes in the fixed parameters can lead to new different scenarios. As it
considered only one parameter for the optimisation, the whole approach is
solving the low-dimensional problem. This approach is not scalable with more
parameters.

The proposed method is a scalable approach compared to previous work[17]. In
[17], pedestrian speed was the only parameter that was varied for the
optimisation process. The other parameters such as the starting position of
ego vehicle, the pedestrian starting position, and weather were held constant.
The issue here is that we do not know the optimal values for these constant
parameters to generate challenging scenarios. This low dimensional
representation of the problem limited the ability to find the challenging
scenarios. In our proposed method, the architecture enables us to add a set of
parameters with distributions over realistic values. During the optimisation
process, the proposed method searches over the parameter values. In this
experiment, five parameters each defined for a valid range of values are
considered. This approach is more scalable as it can handle the high-
dimensional problem without exponential increase in computational
requirements.

Figure 15: The figure depicts the average reward per episode for the
experiment to evaluate the scalability. We have added two more parameters to
the search space to see whether the learning time has a significant change. It
is evident from the figure that achieving the highest reward did not take a
substantial number of episodes.

To demonstrate the ability of the system to handle higher dimensional
problems, we included two additional parameters which impact on the generation
of the scenarios. The first additional parameter is designed to cause abrupt
changes to the pedestrian speed for five timesteps. The second parameter was
included to indicate a timestep in which the unexpected change was to
commence. The values for two additional parameters added as follows:

ped_speed_change:$[0.50,-0.50,0.75,-0.75]\\\ $

ped_timesteps: $[20,30,40,50,60]\\\ \\\ $

As shown in Figure 15, adding more parameters does not increase the number of
episodes required to converge to an optimal policy that has maximum return.
The average reward per episode is similar to the results in Figure 10, which
was the experiment using five parameters. There is a substantial increases in
the possible combination of parameter values with 5 parameters from
$12^{5}=248832$ to $9^{7}=4782969$ with 7 parameters.This demonstrates that
our approach is scalable as the increase in parameters does not increase the
computational requirements by a significant amount.

One of the contributions of the proposed method is to create more realistic
scenarios by using the parameter distributions measured from existing real-
world data. As shown in the figures 11 and 13, the speed and acceleration of
the pedestrian is more realistic as compared to the previous work [17][18].
The acceleration and velocity is bounded for the pedestrian motion while they
are crossing the road, better reflecting the real world physical constraints.
The pedestrian trajectory is created at every episode based on the selected
value of the velocity parameter, ped-vel and acceleration parameter, ped-
accel. As shown in Figures 11(a) and 13(a), the pedestrian trajectories are
very similar, the main difference being the differences in speed and
acceleration. Using trajectories based on real world physical constraints
makes the validation of SUT more efficient as iterations of the optimiser are
not wasted on unrealistic and infeasible trajectories.

##  VI Discussion and Conclusion

This paper presented a falsification method to evaluate a system under test
with a number of influencing parameters that is more efficient when compared
to existing random search or brute-force search based methods. Overall, the
scenarios generated by the proposed method are more realistic than previous
works as appropriate distributions for the parameters representing the
pedestrian kinematics are incorporated to produce smooth, feasible pedestrian
trajectories.

In future work, we will be increasing the complexity of the experiments by
incorporating more parameters to evaluate the scalability of the algorithms
further. Using our experimental autonomous vehicle platforms, we hope to
assess the on-board safety systems to compare the performance to the simulated
results found in this paper. Finding critical scenarios with a specific design
indicates that we have to redesign the SUT to provide a higher standard of
safety. We aim to demonstrate that we can effectively influence the design of
SUT using this approach.

## Acknowledgment

This work has been funded by the Australian Centre for Field Robotics (ACFR),
University of Sydney and Insurance Australia Group (IAG) and iMOVE CRC and
supported by the Cooperative Research Centres program, an Australian
Government initiative.

## References

  * [1] M. Soledad Elli, I. Alvarez, T. Stolte, F. Wirth, C. Stiller, M. Maurer, and J. Weast, “Workshop on automated vehicle safety: Verification, validation and transparency,” in _2020 IEEE Intelligent Transportation Systems Conference (ITSC)_. IEEE, 2020.
  * [2] S. Shalev-Shwartz, S. Shammah, and A. Shashua, “On a formal model of safe and scalable self-driving cars,” _CoRR_ , vol. abs/1708.06374, 2017. [Online]. Available: http://arxiv.org/abs/1708.06374
  * [3] I. Alvarez, “Workshop on automated vehicle safety: Verification, validation and transparency,” in _2019 IEEE Intelligent Transportation Systems Conference (ITSC)_. IEEE, 2019.
  * [4] D. Zhao, H. Lam, H. Peng, S. Bao, D. J. LeBlanc, K. Nobukawa, and C. S. Pan, “Accelerated evaluation of automated vehicles safety in lane-change scenarios based on importance sampling techniques,” _IEEE transactions on intelligent transportation systems_ , vol. 18, no. 3, pp. 595–607, 2016.
  * [5] A. Sarkar and K. Czarnecki, “A behavior driven approach for sampling rare event situations for autonomous vehicles,” _arXiv preprint arXiv:1903.01539_ , 2019.
  * [6] H. Elrofai, D. Worm, and O. O. den Camp, “Scenario identification for validation of automated driving functions,” in _Advanced Microsystems for Automotive Applications 2016_. Springer, 2016, pp. 153–163.
  * [7] N. Kalra and S. M. Paddock, “Driving to safety: How many miles of driving would it take to demonstrate autonomous vehicle reliability?” _Transportation Research Part A: Policy and Practice_ , vol. 94, pp. 182–193, 2016.
  * [8] T. Menzel, G. Bagschik, and M. Maurer, “Scenarios for development, test and validation of automated vehicles,” in _2018 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2018, pp. 1821–1827.
  * [9] C. Amersbach and H. Winner, “Functional decomposition—a contribution to overcome the parameter space explosion during validation of highly automated driving,” _Traffic injury prevention_ , vol. 20, no. sup1, pp. S52–S57, 2019\.
  * [10] E. de Gelder and J.-P. Paardekooper, “Assessment of automated driving systems using real-life scenarios,” in _2017 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2017, pp. 589–594.
  * [11] E.-S. Consortium _et al._ , “Enable-s3 european project,” _available at: www. enable-s3. eu/(accessed 22 January 2018)_ , 2016.
  * [12] A. Pütz, A. Zlocki, J. Bock, and L. Eckstein, “System validation of highly automated vehicles with a database of relevant traffic scenarios,” _situations_ , vol. 1, pp. 19–22, 2017.
  * [13] C. E. Tuncali and G. Fainekos, “Rapidly-exploring random trees for testing automated vehicles,” in _2019 IEEE Intelligent Transportation Systems Conference (ITSC)_. IEEE, 2019, pp. 661–666.
  * [14] B. Zoph and Q. V. Le, “Neural architecture search with reinforcement learning,” _arXiv preprint arXiv:1611.01578_ , 2016.
  * [15] Mobileye, “Implementing the rss model on nhtsa pre-crash scenarios,” Tech. Rep., 2018.
  * [16] C. E. Tuncali, G. Fainekos, H. Ito, and J. Kapinski, “Simulation-based adversarial test generation for autonomous vehicles with machine learning components,” in _2018 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2018, pp. 1555–1562.
  * [17] D. Karunakaran, S. Worrall, and E. Nebot, “Efficient statistical validation with edge cases to evaluate highly automated vehicles,” _arXiv preprint arXiv:2003.01886_ , 2020.
  * [18] M. Koren, S. Alsaif, R. Lee, and M. J. Kochenderfer, “Adaptive stress testing for autonomous vehicles,” in _2018 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2018, pp. 1–7.
  * [19] International Standardization Organization(ISO), “ISO 26262-1:2018,” 2018\. [Online]. Available: https://www.iso.org/standard/68383.html
  * [20] O. Kirovskii and V. Gorelov, “Driver assistance systems: analysis, tests and the safety case. iso 26262 and iso pas 21448,” in _IOP Conference Series: Materials Science and Engineering_ , vol. 534, no. 1. IOP Publishing, 2019, p. 012019.
  * [21] International Standardization Organization(ISO), “Road vehicles — Safety of the intended functionality,” 2019. [Online]. Available: https://www.iso.org/standard/70939.html
  * [22] Aptiv, “Safety first for automated driving,” Tech. Rep., 2019.
  * [23] P. Koopman, U. Ferrell, F. Fratrik, and M. Wagner, “A safety standard approach for fully autonomous vehicles,” in _International Conference on Computer Safety, Reliability, and Security_. Springer, 2019, pp. 326–332.
  * [24] P. Koopman and M. Wagner, “Challenges in autonomous vehicle testing and validation,” _SAE International Journal of Transportation Safety_ , vol. 4, no. 1, pp. 15–24, 2016.
  * [25] N. Euro, “Euro ncap test protocol–aeb systems,” _no. July_ , 2013.
  * [26] H. Peng and D. Leblanc, “Evaluation of the performance and safety of automated vehicles,” _White Pap. NSF Transp. CPS Work_ , 2012.
  * [27] A. Erdogan, E. Kaplan, A. Leitner, and M. Nager, “Parametrized end-to-end scenario generation architecture for autonomous vehicles,” in _2018 6th International Conference on Control Engineering & Information Technology (CEIT)_. IEEE, 2018, pp. 1–6.
  * [28] N. Weber, D. Frerichs, and U. Eberle, “A simulation-based, statistical approach for the derivation of concrete scenarios for the release of highly automated driving functions,” in _AmE 2020-Automotive meets Electronics; 11th GMM-Symposium_. VDE, 2020, pp. 1–6.
  * [29] C. Amersbach and H. Winner, “Defining required and feasible test coverage for scenario-based validation of highly automated vehicles,” in _2019 IEEE Intelligent Transportation Systems Conference (ITSC)_. IEEE, 2019, pp. 425–430.
  * [30] H. Elrofai, J.-P. Paardekooper, E. de Gelder, S. Kalisvaart, and O. O. den Camp, “Scenario-based safety validation of connected and automated driving,” _Netherlands Organization for Applied Scientific Research, TNO, Tech. Rep_ , 2018.
  * [31] L. Waymo, “On the road to fully self-driving,” _Waymo Safety Report_ , pp. 1–43, 2017.
  * [32] R. S. Sutton and A. G. Barto, _Reinforcement learning: An introduction_. MIT press, 2018.
  * [33] M. Van Otterlo and M. Wiering, “Reinforcement learning and markov decision processes,” in _Reinforcement Learning_. Springer, 2012, pp. 3–42.
  * [34] S. Levine, “Cs 294: Deep reinforcement learning,” 2017.
  * [35] A. Dosovitskiy, G. Ros, F. Codevilla, A. Lopez, and V. Koltun, “CARLA: An open urban driving simulator,” in _Proceedings of the 1st Annual Conference on Robot Learning_ , 2017, pp. 1–16.
  * [36] S. Chandra and A. K. Bharti, “Speed distribution curves for pedestrians during walking and crossing,” _Procedia-Social and Behavioral Sciences_ , vol. 104, pp. 660–667, 2013.
  * [37] M. Liu, W. Zeng, P. Chen, and X. Wu, “A microscopic simulation model for pedestrian-pedestrian and pedestrian-vehicle interactions at crosswalks,” _PLoS one_ , vol. 12, no. 7, p. e0180992, 2017.


# The MatrixSchubert Package for Macaulay2

Ayah Almousa  School of Mathematics, University of Minnesota, Minneapolis, MN,
55455<EMAIL_ADDRESS>https://sites.google.com/view/ayah-almousa ,  Sean
Grate  Department of Mathematics and Statistics, Auburn University, Auburn,
AL, 36849<EMAIL_ADDRESS>https://seangrate.com/ ,  Daoji Huang  School
of Mathematics, University of Minnesota, Minneapolis, MN, 55455
<EMAIL_ADDRESS>https://www.daojihuang.me/ ,  Patricia Klein  Department of
Mathematics, Texas A &M University, College Station, TX, 77840
<EMAIL_ADDRESS>https://patriciajklein.github.io/ ,  Adam LaClair  Department
of Mathematics, Purdue University, West Lafayette, IN, 47907
<EMAIL_ADDRESS>https://sites.google.com/view/adamlaclair/home ,  Yuyuan
Luo  Department of Mathematics, Massachusetts Institute of Technology
<EMAIL_ADDRESS>https://www.mit.edu/~lyuyuan/ and  Joseph McDonough  School of
Mathematics, University of Minnesota, Minneapolis, MN, 55455<EMAIL_ADDRESS>https://jmcdonough98.github.io/

###### Abstract.

We introduce the MatrixSchubert package for the computer algebra system
Macaulay2. This package has tools to construct and study matrix Schubert
varieties and alternating sign matrix (ASM) varieties. The package also
introduces tools for quickly computing homological invariants of such
varieties, finding the components of an ASM variety, and checking if a union
of matrix Schubert varieties is an ASM variety.

###### Key words and phrases:

This work was begun at a Macaulay2 mini-school and workshop funded by NSF
grant DMS-2302476 and hosted by the University of Minnesota. We thank the NSF
for its support and the University of Minnesota for its hospitality.

PK is partially supported by NSF grant DMS-2246962.

AL was partially supported by NSF grant DMS-2100288 and by Simons Foundation
Collaboration Grant for Mathematicians #580839.

## 1\. Introduction

Fulton [Ful92] introduced matrix Schubert varieties in the study of Schubert
varieties in the complete flag variety. In that paper, Fulton showed that
matrix Schubert varieties are Cohen–Macaulay and gave an attractive
description of their codimension and defining equations from a combinatorial
perspective. Since that time, there has been a great deal of interest in their
algebraic, geometric, and combinatorial properties (see, for example, [KM05,
KMY09, Hsi13, EM16, FRS16, RRR+21, PSW22, Por23, RRW23]).

More recently, Weigandt [Wei17] introduced alternating sign matrix (ASM)
varieites, which generalize matrix Schubert varieties, and gave a
combinatorial description of their defining equations that generalizes
Fulton’s description in the case of matrix Schubert varieties (equivalently,
when the ASM happens to be a permutation matrix). There are to date no
combinatorial descriptions of the codimension of an ASM variety nor a
combinatorial criterion to determine if it is Cohen–Macaulay or even unmixed.

The MatrixSchubert package implements many basic functions for permutations
and ASMs (such as extending a partial ASM to an ASM, checking if one
permutation avoids another, and finding the descent set of a permutation) as
well as more complicated ones that rely on one or several theoretical results.
Several such examples are described in greater detail in Sections 5 and 7.

One goal of this package is to implement core results from the past three
decades on commutative algebraic aspects of Schubert calculus (particularly
those found in [Ful92, KM05, Wei17, PSW22]) for the purpose of enjoying the
theory that has already been developed. Another is to facilitate inquiry into
the many open questions that remain, including those surrounding resolutions
and Betti numbers of matrix Schubert varieties and ASM varieties,
Cohen–Macaulayness of ASM varieties, and codimension of ASM varieties.

## 2\. Background and some basic functions

### 2.1 Permutations

Let $S_{n}$ denote the group of permutations of $n$ letters. For a permutation
$w\in S_{n}$, we call the matrix that has $1$’s in the positions $(i,w(i))$
and $0$’s in all other positions the _permutation matrix_ of $w$. (As a
cautionary note, some authors take the convention that the matrix described
here - and used throughout the package MatrixSchubert \- is the permutation
matrix of $w^{-1}$.)

Given $w\in S_{n}$, the _Rothe diagram_ of $w$ is

$D(w)=\\{(i,j):i,j\in[n],w(i)>j,\text{ and }w^{-1}(j)>i\\}.$

The Coxeter length of $w$ satisfies $\ell(w)=|D(w)|$. The _essential set_ of
$w$ is

$\operatorname{Ess}(w)=\\{(i,j)\in D(w):(i+1,j),(i,j+1)\not\in D(w)\\},$

i.e., the maximally southeast elements of the connected components of $D(w)$.

    i1 : needsPackage "MatrixSchubert"
    i2 : w = {2,1,5,4,3};
    i3 : rotheDiagram w
    o3 = {(1, 1), (3, 3), (3, 4), (4, 3)}
    o3 : List
    i4 : essentialSet w
    o4 = {(1, 1), (3, 4), (4, 3)}
    o4 : List

### 2.2 Alternating Sign Matrices

A _partial alternating sign matrix_ (partial ASM) is a matrix with entries in
$\\{-1,0,1\\}$ so that partial sums taken along each row (and column) are all
$0$ or $1$. If the entries of each row (and column) sum to $1$ (in which case
the matrix must be square), we call the partial ASM an ASM. The ASMs whose
entries all lie in $\\{0,1\\}$ are exactly the permutation matrices.

The _rank function_ of the $m\times n$ partial ASM $A=(A_{i,j})$ is defined by
$\operatorname{rk}_{A}(a,b)=\sum_{i=1}^{a}\sum_{j=1}^{b}A_{i,j}$ for $1\leq
a\leq m$, $1\leq b\leq n$. (Note that if $A$ is not a permutation matrix, then
$\operatorname{rk}_{A}(a,b)$ may not be the rank of the submatrix of $A$
consisting of its first $a$ rows and first $b$ columns.)

Given a matrix $M$, let $M_{[i],[j]}$ be the submatrix of $M$ consisting of
the first $i$ rows and $j$ columns. Given an $m\times n$ partial ASM $A$, we
define the _ASM variety_ of $A$ to be

$X_{A}=\\{M\in\mbox{mat}(m,n):\operatorname{rk}(M_{[i],[j]})\leq\operatorname{rk}_{A}(i,j)\text{
for all }(i,j)\leq(m,n)\\}.$

If $A\in S_{n}$, we call $X_{w}$ a _matrix Schubert variety_. For background
on matrix Schubert varieties, including a geometric motivation for their
definition and a description of their connection to Schubert varieties, see
[Ful92, MS05].

We call the defining radical of $X_{A}$ in $R=\kappa[z_{1,1},\ldots,z_{m,n}]$,
$\kappa$ an arbitrary field, an _ASM ideal_ or, when $A\in S_{n}$, a _Schubert
determinantal ideal_.

######  Proposition 2.1 ([Ful92, Proposition 3.3]).

For $w\in S_{n}$, $I_{w}$ is prime, $\mbox{codim}(I_{w})=\ell(w)$, $R/I_{w}$
is Cohen–Macaulay.

The function schubertCodim $w$ computes the codimension of $I_{w}$ (that is,
the codimension of spec$(R/I_{w})$ in spec$(R)$) using [Ful92, Proposition
3.3] and the equality $\ell(w)=|D(w)|$. Returning to the example $w=21543$,
one may compare the outcome below to the size of $D(w)$, computed above.

    i5 : schubertCodim w
    o5 = 4

In [Ful92, Lemma 3.10], Fulton gave a generating set for the Schubert
determinantal ideal $I_{w}$. Fix an $m\times n$ generic matrix $Z=(z_{i,j})$.
We write $I_{k}(Z_{[i],[j]})$ for the ideal of $R$ generated by the $k$-minors
in the submatrix $Z_{[i],[j]}$ of $Z$ consisting of its first $i$ rows and
first $j$ columns. Then

$I_{w}=\sum_{(i,j)\in\operatorname{Ess}(w)}I_{\operatorname{rk}_{w}(i,j)+1}(Z_{[i],[j]}),$

and we call these generators the _Fulton generators_. There is a
generalization of the Fulton generators for ASM ideals (see [Wei17, Lemma
5.9]).

The function schubertDeterminantalIdeal takes in either a permutation (as a
list representing its one-line notation) or a partial ASM matrix and produces
its ASM ideal via its Fulton generators (which typically do not form a minimal
generating set).

    i6 : v = {3,1,4,2};
    i7 : schubertDeterminantalIdeal v
    o7 = ideal (z   , z   , - z   z   + z   z  , - z   z   + z   z    , - z   z   + z   z   )
                 1,1   1,2    1,2 2,1   1,1 2,2     1,2 3,1   1,1  3,2     2,2 3,1   2,1 3,2

A term order on $R$ so that the lead term of the determinant of any submatrix
of $Z$ is the product of terms along its antidiagonal is called an
_antidiagonal_ term order.

We now state a pivotal result, due to Knutson and Miller in the case of
Schubert determinantal ideals. One can extend the result to ASM ideals using
either using Frobenius splitting [Knu09] or a combinatorial argument [Wei17].
For full details on the latter, see [KW].

######  Theorem 2.2.

[KM05, Knu09, Wei17, KW] Fix a partial ASM $A$. The Fulton generators of
$I_{A}$ form a Gröbner basis under any antidiagonal term order $<$.
Consequently, $\operatorname{in}_{<}(I_{A})$ is radical.

Rather than computing an ASM ideal and then afterwards computing a Gröbner
basis, these theorems allow us to get an antidiagonal initial ideal directly
from the permutation or ASM.

    i8 : A = matrix{{0,0,1,0},{1,0,-1,1},{0,0,1,0},{0,1,0,0}};
    i9 : antiDiagInit A
    o9 = monomialIdeal (z   , z   , z   z    , z   z   , z   z   )
                         1,1   1,2   1,3  2,1  1,3  2,2  2,2  3,1

## 3\. Rank tables

The rank function $\operatorname{rk}_{A}(i,j)$ of an $m\times n$ partial ASM
used in the definition ASM ideal can be applied to all $(i,j)\leq(m,n)$. The
entire function can be computed and recorded in an $m\times n$ matrix,
implemented via the method rankTable.

    i10 : M = matrix{{0,1,0},{1,-1,0}};
                   2        3
    o10 : Matrix ZZ  <--- ZZ
    i11 : rankTable M
    o11 = | 0 1 1 |
          | 1 1 1 |
                   2        3
    o11 : Matrix ZZ  <--- ZZ

Given any rank table that could be constructed from a partial ASM, the method
rankTableToASM produces the unique partial ASM of the same size as the input
having that rank table.

    i12 : rankTableToASM matrix{{0,1,1},{0,1,1},{1,2,2}}
    o12 = | 0 1 0 |
          | 0 0 0 |
          | 1 0 0 |
                   3        3
    o12 : Matrix ZZ  <--- ZZ

For the typical ASM, many different rank tables could be used to construct the
same ASM variety $X_{A}$. The method rankTableToASM expects the (unique)
matrix with minimum possible entries, which is the one constructed in Section
2.2 (see [RR86, Lemma 1]). If the user has a non-minimal rank table (as a
matrix of non-negative integers), the function rankTableFromMatrix transforms
the non-minimal rank table into a minimal rank table.

    i13 : rankTableFromMatrix matrix{{0,1,2},{0,4,1},{8,2,4}}
    o13 = | 0 1 1 |
          | 0 1 1 |
          | 1 2 2 |
                   3        3
    o13 : Matrix ZZ  <--- ZZ

Among other uses, rank tables also facilitate efficient addition of ASM
ideals. Every sum of ASM ideals is again an ASM ideal [Wei17, Section 3], and
the rank table of the sum is the entrywise minimum of the rank tables of the
ASMs appearing as summands. If $A_{1},\ldots,A_{k}$ are all $m\times n$
partial ASMs, the method schubertAdd computes the rank tables of the $A_{i}$
(via rankTable), takes entrywise minima (via entrywiseMinRankTable), and
computes an ASM ideal from that rank table. Both the ASM and its rank table
are saved in the cache of the newly computed ASM ideal.

    i14 : N = matrix{{1,0,0},{0,0,1}}
    o14 = | 1 0 0 |
          | 0 0 1 |
                  2        3
    o14 : Matrix ZZ  <--- ZZ
    i15 : idealSum = schubertAdd{M,N}
    o15 = ideal (z   , - z   z    + z    z   )
                  1,1     1,2  2,1   1,1  2,2
    o15 : Ideal of QQ[z   ..z   ]
                       1,1   2,2
    i16 : peek idealSum.cache
    o16 = CacheTable{ASM => | 0 1  0 0 |     }
                            | 1 -1 0 1 |
                            | 0 1  0 0 |
                            | 0 0  1 0 |
                     rankTable => | 0 1 1 1 |
                                  | 1 1 1 2 |
                                  | 1 2 2 3 |
                                  | 1 2 3 4 |
    i17 : getASM idealSum
    o17 = | 0 1  0 0 |
          | 1 -1 0 1 |
          | 0 1  0 0 |
          | 0 0  1 0 |

                   4        4
    o17 : Matrix ZZ  <--- ZZ

## 4\. Pattern Avoidance

The MatrixSchubert package has functions to test pattern avoidance for
permutations.

A permutation is called _vexillary_ if it avoids the permutation $2143$. The
class of one-sided ladder determinantal ideals coincides exactly with the
class of vexillary matrix Schubert varieties. The vexillary condition has a
large number of equivalent definitions. We direct the reader to [KMY09,
Section 3.2] for many of them. Testing if a permutation is vexillary is
implemented via the isVexillary function.

    i18 : w = {7,2,5,8,1,3,6,4};
    i19 : isVexillary w
    o19 = false
    i20 : w = {1,6,9,2,4,7,3,5,8};
    i21 : isVexillary w
    o21 = true

A permutation $w$ is _CDG_ if it avoids all eight of the following patterns:

$13254,21543,214635,215364,215634,241635,315264,4261735.$

The class of CDG permutations was named in [HPW22], where a diagonal Gröbner
basis was conjectured for the class (proved in [Kle]). Every CDG permutation
is vexillary, and the CDG permutations form the largest named class of
permutations for which a diagonal Gröbner basis of their matrix Schubert
varieties is known. Testing if a permutation is CDG is implemented via the
isCDG function.

    i22 : w = {5,7,2,1,6,4,3};
    i23 : isCDG w
    o23 = false
    i24 : w = {1,3,5,7,2,4,6};
    i25 : isCDG w
    o25 = true

We say that $w$ is _Cartwright–Sturmfels_ if it avoids all of the following
twelve patterns:

$12543,13254,13524,13542,21543,125364,125634,215364,215634,315264,315624,315642.$

For background on the Cartwright–Sturmfels property of ideals in general, see
[CDG22]. For a proof that the Cartwright-Sturmfels property is characterized
by the given pattern avoidance condition, that Cartwright-Sturmfels Schubert
determinantal ideals have an universal Gröbner basis, and that any initial
ideal of a Cartwright-Sturmfels Schubert determinantal ideal is
Cohen–Macaulay, see [CDG22, Theorem 4.6]. Every Cartwright–Sturmfels Schubert
determinantal ideal is CDG.

Testing if a permutation is Cartwright-Sturmfels is implemented via the
isCartwrightSturmfels functions.

    i26 : w = {3,1,2,6,5,4};
    i27 : isCartwrightSturmfels w;
    o27 = false
    i28 : w = {6,3,5,2,1,4};
    i29 : isCartwrightSturmfels w
    o29 = true

More generally, avoidsAllPatterns inputs a permutation and a list of patterns
to avoid, and determines if the permutation avoids all of the patterns. This
allows users to test conjectures related to pattern avoidance in a much more
general capacity.

## 5\. Algorithms for Castelnuovo–Mumford regularity

Castelnuovo–Mumford regularity is a fundamental invariant in commutative
algebra and algebraic geometry that in a rough sense gives a measure of the
complexity of a module or sheaf.

In this package we implement (as $\mathtt{schubertRegularity}$) a purely
combinatorial formula developed by Peckenik, Speyer, and Weigandt [PSW22,
Theorem 1.2] for computing the Castelnuovo–Mumford regularity of $S/I_{w}$ for
arbitrary $w\in S_{n}$. We also extend the functionality of
$\mathtt{schubertRegularity}$ so that it can compute the Castelnuovo–Mumford
regularity of the coordinate ring associated to a partial ASM by passing to
the antidiagonal initial ideal, a valid strategy in light of an important
theorem of Conca and Varbaro [CV20] together with Theorem 2.2. See Subsection
5.2 for a fuller explanation.

### 5.1 Matrix Schubert varieties

The theoretical foundation for the $\mathtt{schubertRegularity}$ function is a
result of Pechenik, Speyer, and Weigandt:

######  Theorem 5.1 ([PSW22, Theorem 1.2]).

For $w\in S_{n}$,

$\displaystyle\operatorname{reg}(S/I_{w})=\operatorname{raj}(w)-|D(w)|$

where $\operatorname{raj}(w)$ is the Rajchgot index of a permutation.

We refer the reader to [PSW22] for the definition of the Rajchgot index.

The function $\mathtt{schubertRegularity}$ takes either a permutation in one-
line notation (that is, as a list) or a partial ASM and returns the
Castelnuovo–Mumford regularity of the associated coordinate ring. The
computation of the Rajchot index of a permutation invovles determining longest
subsequences of the permutation subject to certain conditions. We utilize
memoization in determining these longest subsequences which leads to drastic
speed improvements over Macaulay2’s built-in command $\mathtt{regularity}$.
(By that, we mean that there are examples of the type we would be inclined to
compute in the course of research which run reliably faster, not that we have
performed any type of formal efficiency analysis. We include a couple of
illustrating examples here.) The following example run on a computer with an
AMD Ryzen 5 5600U processor demonstrates the extent of this speed-up.

    i30 : w = {1,2,3,9,8,4,5,6,7};
    i31 : I = antiDiagInit(w, CoefficientRing=>ZZ/3001);
                            ZZ
    o31 : MonomialIdeal of ----[z   ..z   ]
                           3001  1,1   9,9
    i32 : M = comodule I;
    i33 : time regularity M
         -- used 1.41504 seconds
    o33 = 6
    i34 : time schubertRegularity w
         -- used 0.000181139 seconds
    o34 = 6
    i35 : time schubertRegularity random toList (1 .. 100)
         -- used 4.03299 seconds
    o35 = 1925

### 5.2 ASM Varieties

The function $\mathtt{schubertRegularity}$ also accepts as input a partial
ASM. First, $\mathtt{schubertRegularity}$ checks whether the matrix is a
permutation matrix; in which case the Castelnuovo–Mumford regularity is
computed via Theorem 5.1. Otherwise, the antidiagonal initial ideal of the ASM
is computed, and the built-in Macaulay2 command $\mathtt{regularity}$ is used.
The Castelnuovo–Mumford regularity of the quotient by the initial ideal will
coincide with the Castelnuovo–Mumford regularity of the coordinate ring
corresponding to the ASM by [CV20, Corollary 2.7] since the antidiagonal
initial ideal is squarefree [KM05, Wei17].

    i36 : A = matrix{{0,0,1,0},{0,1,-1,1},{1,-1,1,0},{0,1,0,0}};
                  4        4
    o36 : Matrix ZZ  <--- ZZ
    i37 : time regularity comodule schubertDeterminantalIdeal A
         -- used 0.00968312 seconds
    o37 = 1
    i38 : time schubertRegularity A
         -- used 0.0100184 seconds

    o38 = 1
    i39 : B = matrix{{1,0,0,0,0,0,0,0},{0,1,0,0,0,0,0,0},{0,0,0,0,1,0,0,0},
             {0,0,0,0,0,1,0,0},{0,0,1,0,0,-1,1,0},{0,0,0,1,-1,1,0,0},
             {0,0,0,0,1,0,0,0},{0,0,0,0,0,0,0,1}};
                  8        8
    o39 : Matrix ZZ  <--- ZZ
    i40 : time regularity comodule schubertDeterminantalIdeal B
         -- used 1.01169 seconds
    o40 = 8
    i41 : time schubertRegularity B
         -- used 0.08511 seconds
    o41 = 8

We are a little bit disappointed that, as we see above with the ASM $A$, the
command schubertRegularity is sometimes slower than the already-available
regularity command applied to an ASM ideal. However, as we see above with the
ASM $B$, schubertRegularity can be a meaningful improvement over regularity in
other cases.111We expect that the issue is that we are for some reason
describing the sets of rows and columns whose minors define the ASM ideal (or
its antidiagonal initial ideal) less efficiently than the minors command does.
This explanation is compatible with our experience that antiDiagInit is
typically slower than schubertDeterminantalIdeal for dominant permutations,
i.e., permutations indexing Schubert determinantal ideals that are already
monomial ideals. As a rule of thumb, we recommend schubertRegularity for ASM
ideals with more generators of higher degrees and regularity for ASM ideals
with fewer generators or generators in lower degrees.

## 6\. Schubert and Grothendieck polynomials

The MatrixSchubert package provides functions to compute Schubert, double
Schubert, and Grothendieck polynomials for permutations. We give a brief
overview of how these families of polynomials are constructed, and refer the
reader to [KM05] for a more detailed (and somewhat more general) treatment.
Let $\mathbf{x}=\\{x_{1},\ldots,x_{n}\\}$, and let $S$ denote the polynomial
ring in $\mathbf{x}$ over the field $\kappa$. Let $\partial_{i}$ be the $i$’th
divided difference operator, which sends $f\in S$ to

$\partial_{i}(f)=\frac{f(x_{1},\ldots,x_{n})-f(x_{1},..,x_{i+1},x_{i},..,x_{n})}{x_{i}-x_{i+1}}.$

It is not obvious that $\partial_{i}(f)\in S$, but it is true. Let
$w_{0}=nn-1\cdots 321$, i.e., the longest word in $S_{n}$. For $w\in S_{n}$,
the _Schubert polynomial_ $\mathfrak{S}_{w}$ is defined recursively as
follows:

$\mathfrak{S}_{w_{0}}(\mathbf{x})=\prod_{i=1}^{n}x_{i}^{n-i}\in\mathbb{Z}[\mathbf{x}]\qquad\text{and}\qquad\mathfrak{S}_{ws_{i}}(\mathbf{x})=\partial_{i}\mathfrak{S}_{w}(\mathbf{x})$

where $s_{i}$ is a right descent of $w$ (i.e., $w(i)>w(i+1)$). The _double
Schubert polynomial_ $\mathfrak{S}_{w}(\mathbf{x},\mathbf{y})$ is defined
using the same recursion, but with initial condition
$\mathfrak{S}_{w_{0}}(\mathbf{x},\mathbf{y})=\prod_{i+j\leq
n}(x_{i}-y_{j})\in\mathbb{Z}[\mathbf{y}][\mathbf{x}]$. Finally, the
_Grothendieck polynomial_ $\mathcal{G}_{w}(\mathbf{x})$ is defined by the
recursion

$\mathcal{G}_{w_{0}}(\mathbf{x})=\prod_{i=1}^{n}x_{i}^{n-i}\qquad\text{and}\qquad\mathcal{G}_{ws_{i}}(\mathbf{x})=\partial_{i}(\mathcal{G}_{w}(\mathbf{x})-x_{i+1}\mathcal{G}_{w}(\mathbf{x}))$

where $s_{i}$ is a right descent of $w$. Computing these polynomials for a
permutation $w\in S_{n}$ given in one-line notation is implemented via the
schubertPolynomial, doubleSchubertPolynomialnomial, and grothendieckPolynomial
functions respectively.

    i42 : w = {2,1,4,3};
    i43 : schubertPolynomial w
           2
    o43 = x  + x x  + x x
           1    1 2    1 3
    o43 : QQ[x ..x ]
              1   4
    i44 : doubleSchubertPolynomialnomial w
           2                                        2
    o44 = x  + x x  + x x  - 2x y  - x y  - x y  + y  - x y  + y y  - x y  + y y
           1    1 2    1 3     1 1    2 1    3 1    1    1 2    1 2    1 3    1 3
    o44 : QQ[x ..x , y ..y ]
              1   4   1   4
    i45 : grothendieckPolynomial w

           2        2      2               2
    o45 = x x x  - x x  - x x  - x x x  + x  + x x  + x x
          1 2 3    1 2    1 3    1 2 3    1    1 2    1 3
    o45 : QQ[x ..x ]
              1   4

The default options for computing Schubert, double Schubert, and Grothendieck
polynomials use the definitions by divided difference operators, where we
deterministically pick one reduced word for each $w$ to apply the divided
difference operators. For Schubert polynomials, we also provide the option
Algorithm=>"Transition" that computes Schubert polynomials via transition
equations; see, e.g., [W21]. For Grothendieck polynomials, we provide three
different implementations: "DividedDifference", "Degree", and "PipeDream". The
"Degree" option computes the twisted $K$-polynomials of the matrix Schubert
variety and should not be used for any practical implementation. The
"PipeDream" option computes Grothendieck polynomials by the pipe dream
formula. We provide these different options for any users who are interested
in comparing efficiency of the different algorithms.

## 7\. Studying ASM varieties via initial ideals

Rank functions induce a lattice structure on the set of $n\times n$ ASMs
defined by $A\geq B$ if and only if
$\operatorname{rk}_{A}(i,j)\leq\operatorname{rk}_{B}(i,j)$ for all
$i,j\in[n]$. The restriction of this partial order to $S_{n}$ recovers
(strong) Bruhat order on $S_{n}$. Define

$\text{perm}(A)=\\{w\in S_{n}:w\geq A,\text{ and, if }w\geq v\geq A\mbox{ for
some $v\in S_{n}$, then }w=v\\}.$

######  Proposition 7.1.

[Wei17, Proposition 5.4],[KW, Lemma 2.6] If $A$ is an ASM and $<$ is an
antidiagonal term order, then

$I_{A}=\bigcap_{w\in\text{perm}(A)}I_{w}\mbox{ \hskip 56.9055pt and \hskip
56.9055pt}\operatorname{in}_{<}(I_{A})=\bigcap_{w\in\text{perm}(A)}\operatorname{in}_{<}(I_{w}).$

By combining Knutson and Miller’s [KM05, Theorem B] with Bergeron and Billey’s
[BB93, Theorem 3.7], one may construct a reduced word for $w\in S_{n}$ from
the indices of the variables generating any minimal prime of
$\operatorname{in}_{<}(I_{w})$.

It is of independent combinatorial interest to understand the lattice of ASMs.
The function permSetOfASM takes in an ASM $A$ and computes $\text{perm}(A)$ by
decomposing the antidiagonal initial ideal of $I_{A}$ and reading a reduced
word from each of the primes appearing in the decomposition. The set of
distinct permutations encountered comprises $\text{perm}(A)$ and therefore
also indexes the components in a prime decomposition of $I_{A}$.

    i46 : A = matrix{{0,1,0},{1,-1,1},{0,1,0}};
    i47 : permSetOfASM A
    o47 = {{3, 1, 2}, {2, 3, 1}}

The function schubertDecompose takes in an ideal, computes its initial ideal
by the default term order in Macaulay2 (which is antidiagonal) and, from the
minimal primes of that ideal, finds and returns the set of permutations with
at least one reduced word given by the set of generators of one of those
minimal primes. The primary use of this function is on an ideal the user knows
to be an ASM ideal $I_{A}$ (such as one arising directly from the matrix $A$
or as a sum of other ASM ideals), in which case the output will be
$\text{perm}(A)$.

    i48 : schubertDecompose schubertDeterminantalIdeal A
    o48 = {{3, 1, 2}, {2, 3, 1}}

If the user has an ideal and is unsure if that ideal is an ASM ideal, they may
use the function isASMIdeal, which takes in an ideal $I$ and first applies
schubertDecompose. It then takes entrywise maxima (using
entrywiseMaxRankTable) among the rank tables of the permutations found from
schubertDecompose. Using rankTableToASM, it constructs the partial ASM $A$
whose ASM ideal is determined by that rank table. Finally, isASMIdeal returns
a boolean indicating if $I_{A}$ is equal to the input ideal. In case it is,
isASMIdeal caches the partial ASM $A$ so that $I=I_{A}$. The ASM $A$ may then
be retrieved via getASM.

    i49 : I1 = schubertDeterminantalIdeal {3,4,1,2};
    o49 : Ideal of QQ[z   ..z   ]
                       1,1   4,4
    i50 : I2 = sub(schubertDeterminantalIdeal {3,2,4,1},ring I1);
    o50 : Ideal of QQ[z   ..z   ]
                       1,1   4,4
    i51 : I = intersect(I1,I2);
    o51 : Ideal of QQ[z   ..z   ]
                       1,1   4,4
    i52 : isASMIdeal I
    o52 = true
    i53 : A = getASM I
    o53 = | 0 0  1 0 |
          | 0 1  0 0 |
          | 1 -1 0 1 |
          | 0 1  0 0 |

The antidiagonal initial ideal of an ASM ideal can be used for additional
computations in light of recent and impactful results of Conca and Varbaro
[CV20]. Fix a partial ASM $A$ and antidiagonal term order $<$. A result from
[CV20] closely related to that discussed in the section on Castelnuovo–Mumford
regulartiy states that $R/I_{A}$ is Cohen–Macaulay if and only if
$R/\operatorname{in}_{<}(I_{A})$ is Cohen–Macaulay (which again uses that
$\operatorname{in}_{<}(I_{A})$ is radical). The function isSchubertCM assesses
Cohen–Macaulayness of $R/I_{A}$ by checking
$\operatorname{pdim}(\operatorname{in}_{<}(I_{A}))==\operatorname{codim}(\operatorname{in}_{<}(I_{A}))$.

We continue from the example above and then consider a familiar non-
Cohen–Macaulay variety.

    i54 : isSchubertCM A
    o54 = true
    i55 : B = matrix{{0,0,1,0,0},{0,0,0,1,0},{1,0,-1,0,1},{0,1,0,0,0},{0,0,1,0,0}};
                     5        5
    o55 : Matrix ZZ  <--- ZZ
    i56 : trim schubertDeterminantalIdeal B
    o56 = ideal (z   , z   , z   , z   , z   z   , z   z   , z   z   , z   z   )
                  2,2   2,1   1,2   1,1   2,3 3,2   1,3 3,2   2,3 3,1   1,3 3,1
    i57 : isSchubertCM B
    o57 = false

The Stanley–Reisner complexes of antidiagonal initial ideals of ASM ideals are
of independent interest. Knutson and Miller [KM05] introduced subword
complexes, of which Stanley–Reisner complexes of antidiagonal initial ideals
of Schubert determinantal ideals are the motivating example. Given a
permutation $w$, the method subwordComplex produces the Stanley–Reisner
complex of the antidiagonal initial ideal of $I_{w}$.

    i58 : w = {2,1,4,3};

    i59 : netList facets subwordComplex w

           +--------------------------------------------------------+
    o59 =  |z   z   z   z   z   z   z   z   z   z   z   z   z   z   |
           | 1,2 1,4 2,1 2,2 2,3 2,4 3,1 3,2 3,3 3,4 4,1 4,2 4,3 4,4|
           +--------------------------------------------------------+
           |z   z   z   z   z   z   z   z   z   z   z   z   z   z   |
           | 1,2 1,3 1,4 2,1 2,3 2,4 3,1 3,2 3,3 3,4 4,1 4,2 4,3 4,4|
           +--------------------------------------------------------+
           |z   z   z   z   z   z   z   z   z   z   z   z   z   z   |
           | 1,2 1,3 1,4 2,1 2,2 2,3 2,4 3,2 3,3 3,4 4,1 4,2 4,3 4,4|
           +--------------------------------------------------------+
    i60 : v = {2,1,6,3,5,4};
    i61 : # facets subwordComplex v
    o61 = 35

Moreover, Knutson and Miller [KM05] showed that the prime components of the
antidiagonal initial ideal of Schubert determinantal ideals are indexed by
combinatorial objects called _pipe dreams_ (by [KM05], now adopted as standard
terminology in the literature) or _RC-graphs_ (by [BB93], to which we refer
the reader for background).

In particular, to read off an associated prime of the antidiagonal initial
ideal from a pipe dream, one simply needs to read off the locations of the $+$
tiles in the pipe dream. This package provides the class PipeDream to display
and manipulate pipe dreams.

    i62 : u = {2,1,4,3,6,5};

    i63 : (pipeDreams u)_0

    o63 = +/+/+/
          //////
          //////
          //////
          //////
          //////

    o63 : PipeDream

    i64 : (decompose antiDiagInit u)_0

    o64 = monomialIdeal (z   , z   , z   )
                         1,1   1,3   1,5

    o64 : MonomialIdeal of QQ[z   ..z   ]
                              1,1   6,6

The generator $z_{1,1}$ corresponds to the $+$ tile in location $(1,1)$ of the
given pipe dream, the generator $z_{1,3}$ to the $+$ tile in location $(1,3)$,
and the generator $z_{1,5}$ to the $+$ tile in location $(1,5)$.

To compare the Macaulay2 drawing of a pipe dream with those in [KM05] and
subsequent literature, + is interpreted as a cross tile, and / is interpreted
as an elbow tile.

Figure 1. One pipe dream for $w=214365$. + is interpreted as a cross tile, and
/ is interpreted as an elbow tile, and entries below the main antidiagonal are
ignored, as is standard.

## 8\. Some diagonal term orders

Knutson and Miller’s [KM05] result that every antidiagonal term order
determines the same initial ideal of a given Schubert determinantal ideal has
been used heavily throughout the `MatrixSchubert` package. By contrast,
different diagonal term orders can yield different initial ideals. The various
diagonal initial ideals of Schubert determinantal ideals are a topic of active
research. Given a permutation in one-line notation or a partial ASM, the
functions `diagLexInitSE`, `diagLexInitNW`, `diagRevLexInit` will each produce
an initial ideal under a distinct diagonal term order.

Precisely, `diagLexInitSE` is the lexicographic order for which $z_{n,n}$ is
largest and the remaining variables are ordered by reading left across the
bottom row, then right to left across row $n-1$, and so on until arriving
finally at $z_{1,1}$, the smallest variable. `diagLexInitNW` is the
lexicographic order for which $z_{1,1}$ is largest and the remaining variables
are ordered by reading right across the top row, then left to right across row
$2$, and so on until arriving finally at $z_{n,n}$, the smallest variable. And
`diagRevLexInit` is the reverse lexicographic order where $z_{n,1}$ is the
smallest (or most penalized) variable followed by the variables encountered
reading left to right along the bottom row, then left to right along row $n-1$
and so on until arriving at $z_{1,n}$.

The example $w=214365$, taken from [KW], is the smallest example of which the
authors are aware of a permutation with different initial ideals for different
diagonal term orders. The authors are unaware of any examples for which
`diagLexInitSE` and `diagRevLexInit` produce different initial ideals. In the
example below, the Macaulay2 output recording the ambient rings of the initial
ideals has been omitted for brevity.

    i65 : w = {2,1,4,3,6,5}
    o65 = {2, 1, 4, 3, 6, 5}
    o65 : List
    i66 : diagLexInitSE w
                                                                     2
    o66 = monomialIdeal (z   z   z   z   z   z   , z   z   z   z   z   ,
                          5,5 4,3 3,4 3,2 2,1 1,3   5,5 4,3 3,4 2,1 1,2
                                            z   z   z   , z   z   z   z   z   z   , z   )
                                            3,3 2,1 1,2   5,5 4,3 3,4 3,1 2,3 1,2   1,1
    i67 : diagLexInitNW w

    o67 = monomialIdeal (z   , z   z   z   , z   z   z   z   z   ,
                          1,1   1,2 2,1 3,3   1,2 2,1 3,4 4,3 5,5
                                      2
                                    z   z   z   z   z   z   , z   z   z   z   z   z   )
                                     1,2 2,3 3,1 3,4 4,3 5,5   1,3 2,1 3,2 3,4 4,3 5,5
    i68 : diagRevLexInit w
                                                                     2
    o68 = monomialIdeal (z   z   z   z   z   z   , z   z   z   z   z   ,
                          5,5 4,3 3,4 3,2 2,1 1,3   5,5 4,3 3,4 2,1 1,2
                                            z   z   z   , z   z   z   z   z   z   , z   )
                                            3,3 2,1 1,2   5,5 4,3 3,4 3,1 2,3 1,2   1,1

## 9\. Available Examples of ASMs

For the convenience of the user, we provide a complete list of ASMs up to
ASM(7). It can be accessed as `ASMFullList n` for any $1\leq n\leq 7$, which
returns a list of objects of type `Matrix`. To access a list of $m$ random
examples of ASMs of size $n$, `ASMRandomList(n,m)` returns a random length-$m$
list of ASMs which are size $n$, presented as `Matrix` objects.

    i69 : ASMRandomList(5,4)

    o69 = {| 1 0 0 0  0 |, | 0 0 0 1  0 |, | 0 0  0 1  0 |, | 0 0 0  0 1 |}
           | 0 0 1 0  0 |  | 1 0 0 -1 1 |  | 0 1  0 -1 1 |  | 0 0 1  0 0 |
           | 0 0 0 1  0 |  | 0 1 0 0  0 |  | 1 -1 1 0  0 |  | 1 0 0  0 0 |
           | 0 1 0 -1 1 |  | 0 0 0 1  0 |  | 0 0  0 1  0 |  | 0 1 -1 1 0 |
           | 0 0 0 1  0 |  | 0 0 1 0  0 |  | 0 1  0 0  0 |  | 0 0 1  0 0 |

Additionally, for $n\leq 6$, lists of non-permutation ASMs that define an
arithmetically Cohen–Macaulay variety, ASMs that do not define an
arithmetically Cohen–Macaulay variety, and antidiagonal initial ideals of ASMs
are provided, and can be accessed with `cohenMacaulayASMsList n`,
`nonCohenMacaulayASMsList n`, and `initialIdealsList n`, respectively.

It is well known that there are $429$ $5\times 5$ ASMs, of which $5!$ are
permutation matrices. Each non-permutation $5\times 5$ ASM is on exactly one
of the lists `cohenMacaulayASMsList 5` or `nonCohenMacaulayASMsList 5`.

    i70 : CM = cohenMacaulayASMsList 5;
    i71 : NCM =  nonCohenMacaulayASMsList 5;
    i72 : #CM+#NCM+5! == 429
    o72 = true
    i73 : #ASMFullList 5
    o73 = 429

## 10\. Acknowledgements

The authors thank Shiliang Gao, Pooja Joshi, and Antsa Tantely Fandresena
Rakotondrafara, contributors to the MatrixSchubert package who provided code,
improved existing code, or enhanced documentation. They also thank Anton
Leykin, Mike Stillman and Gregory Smith, who answered countless questions both
during the Macaulay2 workshop, where this package was begun, and after. They
thank Mahrud Sayrafi for his helpful code review and for implementing the
PipeDream class. The authors thank Anna Weigandt and Jonathan Montaño for
helpful conversations. Finally, the authors would also like to thank Christine
Berkesch, Michael Perlman, and Mahrud Sayrafi for co-organizing the workshop
with the first author where this project began.

## References