Kuiper, Lin, Blanchet, and Tarokh

# Generative Learning for Simulation of Vehicle Faults

Patrick Kuiper1 Sirui Lin2 Jose Blanchet2 and Vahid Tarokh1
1Dept. of Electrical and Computer Eng  Duke University  Durham  NC  USA
2Dept. of Management Science and Eng  Stanford University  Stanford  CA  USA
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
## ABSTRACT

We develop a novel generative model to simulate vehicle health and forecast
faults, conditioned on practical operational considerations. The model,
trained on data from the US Army’s Predictive Logistics program, aims to
support predictive maintenance. It forecasts faults far enough in advance to
execute a maintenance intervention before a breakdown occurs. The model
incorporates real-world factors that affect vehicle health. It also allows us
to understand the vehicle’s condition by analyzing operating data, and
characterizing each vehicle into discrete states. Importantly, the model
predicts the time to first fault with high accuracy. We compare its
performance to other models and demonstrate its successful training.

##  1 INTRODUCTION

Many organizations have realized the utility of using machine learning and a
data driven approach to understand and improve the health of their vehicle
fleets; these organizations invest enormous amounts of resources in
transportation systems as they are critical to operations. We focus this
analysis on the United States’ Department of Defense (DoD), where the US Army
alone is projected to spend an estimated $5 billion per year (in 2020 dollar
terms through 2050), developing and acquiring ground vehicles, where ground
vehicles are any vehicles other than aircraft and ships [CBO (2021)].
Maintaining this enormous investment is critical to ensuring combat readiness
across the DoD, where the department spent $90 billion in 2022 on maintaining
vehicles across domains: ground, air, and sea [GAO (2022)].

Predicting requirements is critical to an effective maintenance program. The
application of statistics towards vehicle maintenance prediction is often
referred to as predictive maintenance. Recognizing the importance of
predictive maintenance, in the 2022 National Defense Authorization Act (NDAA)
Congress required the DoD Inspector General Office to review predictive
maintenance practices, originally established by DoD directives in 2002 and
2007 [DoDIG (2023)]. Given the importance of predictive maintenance in
promoting readiness and conserving budgetary resources, the US Army is
employing hardware and software solutions to address Congressional and
Departmental guidance.

We propose a model designed to achieve _predictive maintenance_ , with the
goal of completing corrective actions on a vehicle in order to avoid more
costly or catastrophic damage in the future [Theissler et al. (2021)]. This
predictive maintenance model is _data driven_ , where conclusions are based on
the outcomes observed in a dataset of sensor and fault signals recorded from
an onboard vehicle computer. This data was collected by the US Army under an
effort designated as Condition Based Maintenance (CBM). In this effort
millions of time indexed data points were collected from a number of military
vehicle types, with the goal of informing the development of maintenance
processes, hardware, and analysis software to support predictive maintenance.
We have developed a generative learning approach trained on a publicly
releasable subset of this CBM data to predict future vehicle faults.
Additionally, we simulate outcomes conditioned on practical vehicle
considerations, including vehicle age, location, and mode sub-type.

We begin the discussion of our model with a brief literature review,
introducing general time-series modeling, then focusing on prediction models
for vehicle faults, and relevant technical reports related to the maintenance
of military vehicles. Next, we review an analysis related to the DoD’s efforts
with predictive maintenance and the CBM dataset used to develop our model. We
then introduce the proposed generative fault prediction model. We conclude
with results, providing prediction performance analysis, comparison to
alternate models, and interpretation of the analysis provided by the proposed
model. Additional information regarding model training convergence and
parameter interpretation are provided in the appendices, including a link to a
github repository with all modeling code and processor information in Appendix
B. In summary, our contributions include:

  1. 1.

We develop and train a practical, generative learning model to accurately
predict vehicle faults and time to first fault on a rolling basis in real
time.

  2. 2.

We generate future sensor feature covariates based on the values of previous
covariates, using an attention based transformer.

  3. 3.

We also generate future sensor feature covariates based on the values of
previous data and a selected vehicle use profile, where the profile is based
upon practical considerations, using a Variational Auto-Encoder and K-means
clustering.

  4. 4.

We provide results and conclusions which are interpretable and understandable
among professional practitioners in the application field.

##  2 Related Work

Time series analysis, which leverages time-indexed data to extract meaningful
statistical conclusions, remains an important field in statistics with
applications in finance, medicine, and manufacturing. The wide adoption of
machine learning techniques, specifically generative learning, has led to the
development of a number of innovative time series methods. We discuss these
techniques, along with their application to reliability engineering in
automobiles.

###  2.1 General Time Series Prediction

A number of effective, well established methods exist to perform time series
analysis; these include Auto-Regressive (AR) and Moving Average (MA) models,
in addition to more advanced integration of both AR and MA, with Auto-
Regressive Moving Average (ARMA) and Auto-Regressive Integrated Moving Average
(ARIMA) [Box and Jenkins (1976)]. Approaches employing neural networks
specifically pertaining to the forecasting of time-indexed sequence data have
included Recurrent Neural Networks (RNNs) and Long-Short Term Memory networks
(LSTMs). Of particular interest to this analysis is the Deep-Auto Regressive
(DeepAR) method developed by salinas2019deepar. The DeepAR method employs a
LSTM-based neural network design, with an AR component, that learns from
multiple time series data points. In the proposed model, multiple time series
data includes the fault target and sensor feature covariates.

The DeepAR method produces a probabilistic function via Monte Carlo
simulation, allowing for the model to be queried on quantiles, which then
generates forecasted future data. This enables the generation of multiple
future predictions via the quantile function, improving the utility of the
analysis when compared to point estimates. A critical consideration is that
the DeepAR model requires the incorporation of _future_ covariates when
forecasting a future target variable. A novel approach for generating
multivariate time series predictions is proposed by attn_mech, with the use of
a spatio-temporal attention mechanism (STAM). This model isolates critical
time steps (temporal) and variables (spatial) when forecasting future target
predictions. It does not depend on future covariates during the prediction
step as with DeepAR; however, we found STAM alone does not perform as well on
forecasting vehicle faults when compared to the DeepAR model (see Section 5
for comparison). We propose a hybrid model using both STAM and DeepAR, where
the requirement of the DeepAR model to incorporate future covariates is
addressed by generating these sensor values using STAM. We develop this
approach further in the Section 4 with additional discussion in Section 5.4.

###  2.2 Vehicle Fault Prediction

Using generative learning to anticipate faults in vehicles remains a
relatively new topic. Recently, mult_time_series completed an analysis similar
to our proposed application, using self-supervised learning and a graph based
technique to forecast vehicle faults. The formulation of the technique
developed by mult_time_series differ significantly to our proposed method, and
the dataset is much smaller, with the faults being labeled by experts visually
reviewing sensor signals. This is in contrast to our analysis which leverages
a much larger dataset collected in non-experimental, field conditions, where
the faults are labeled by an automated vehicle onboard computer.

Another similar analysis is provided by fault_class_discrim, where several
discriminative classification methods are evaluated for determining the
occurrence of vehicle faults in several different subsystems. The evaluated
classification model takes real-time vehicle data in the form of sensor data
and Diagnostic Trouble Codes (DTC), with dimensional reduction using Principal
Component Analysis (PCA). The experiments proposed in this analysis are
similar to our proposal in that data was collected from dozens of vehicles of
a single vehicle family (make / model), with the goal of forecasting faults in
real time using machine learning. However, the model only produces a
discriminative classification, not a generative prediction, and no
interpretable information is provided about the causes of vehicle failure from
this analysis.

A general review of machine learning applications for vehicle predictive
maintenance is proposed by pred_maint_survey. Of note from this rigorous
review and categorization of related work is the proposal of several research
challenges in the field. The authors highlight the need for _public real-world
datasets_ , which are _labeled_. An additional relevant challenge referenced
in this work is the use of more complex, _truly deep neural networks_ for
modeling, while maintaining _interpretability_ and _explainability_ given the
need for understanding when addressing maintenance concerns. We work to
address all of these research challenges in this proposal and provide a
practical model for a large family of vehicles in the US Army, which may be
readily applied to additional vehicle types.

###  2.3 Military Maintenance Background

An detailed report of the CBM dataset, along with a number of other predictive
maintenance efforts established by the US Navy and Marine Corps, is provided
by rit_report. This report provides a wealth of technical detail concerning
the vehicle population covered in the CBM dataset, along with the entire DoD
predictive maintenance effort across the uniformed services.

As a check on the CBM dataset, we reference orcen_report who developed an
insightful analysis into historical operational readiness rates of several
common US Army vehicles. This analysis is relevant as it describes with
precision common operational readiness rates observed among organizations in
the US Army. These published readiness rates serve as a validation for
observed rates recorded in the CBM dataset used in the proposed analysis.

##  3 CBM Dataset Overview

We have identified a research demand for a real-world dataset composed of
vehicle sensor readings, with associated vehicle faults. Such a dataset would
allow for the development of a data driven fault prediction model. This
dataset would be of particular use in military and civilian automotive
applications. We now provide further details about the US Army CBM dataset
which we propose meets the challenges presented by the predictive maintenance
research community, and further develop our assumptions that vehicle
electronic fault signals may be used as a proxy for vehicle condition.

###  3.1 US Military Predictive Maintenance Data Collection Efforts

According to the Government Accountability Office (GAO) report titled
“Military Readiness: Actions Needed to Further Predictive Maintenance on
Weapon Sytems”, published in December 2022 [GAO (2022)], the US DOD’s
predictive maintenance efforts may be summarized as

> …any effort that uses condition-monitoring technology or analysis of
> historical data to anticipate maintenance needs in a manner that reduces
> unscheduled reactive maintenance or overly prescriptive preventive
> maintenance. The military services use multiple terms for predictive
> maintenance or predictive maintenance enablers, including Condition-Based
> Maintenance (CBM), Condition-Based Maintenance Plus (CBM+), Prognostic and
> Predictive Maintenance (PPMx), Enhanced Reliability-Centered Maintenance
> (eRCM), and Predictive Maintenance.

The dataset used in this analysis is from the US Army’s CBM effort. As
illustrated in Figure 1, sensor and fault data was recorded over a period of
months from a number of vehicle types. As referenced in this GAO report,
billions of dollars are spent and military readiness is degraded due to
unplanned maintenance associated with US weapons (“weapons” includes
vehicles). This degradation of resources and readiness may be significantly
reduced through the use of effective predictive maintenance, as dictated in
the 2002 DOD predictive maintenance policy.

Figure 1: Description of Condition Based Maintenance dataset. Data from field
vehicle operations, collected over several months from various vehicle types,
are presented. This dataset includes time series of various fault signals,
each with associated recording timestamps, as well as time series from
multiple sensors also timestamped. Our analysis aims to utilize the sensor
data as feature covariates to predict fault occurrences.

###  3.2 Fault Signals and US Army Vehicles

Vehicle fault targets in this analysis are recorded in accordance with the
Society of Automotive Engineers (SAE) J1939 and J1708 data specifications
[Thurston, Michael, G and McConky, Sean, P and Valant, Chris, J and Nenadic,
Nenad, G (2019)]. It is important to note that the analysis we propose relies
on these _fault signals_ , also known as Diagnostic Trouble Codes (DTC), as a
target and not an observation of vehicle failure. This is because failure data
is not available in the analyzed set. The assumption that vehicle fault
signals or DTC data is highly correlated with actual failure is supported by
the detailed work regarding the CBM dataset completed by rit_report as
described directly below:

> DTCs can be very useful for development and enhancement of PHM (Predictive
> Health Monitoring) capabilities, because they contain the ground truth
> information for impending failures. In fact, they represent a level of
> condition monitoring that is already implemented by electronic control units
> (ECUs), employing traditionally conservative thresholds to reduce false
> alarm problems.

While true vehicle failure would be an informative additional covariate, the
goal of predictive maintenance is to take action based on vehicle signals
before failure, where failures have high correlation with fault signals, as
these signals measure the true status of critical vehicle systems (i.e.
engine, transmission, and brakes).

##  4 Generative Fault Prediction Model Definition

In this section we propose a target distribution of vehicle faults conditioned
on vehicle sensor values, which would provide significant utility towards
predictive maintenance. We formulate a model for this target distribution
employing generative learning. We fully specify this generative learning model
and provide details concerning training and experimentation. The proposed
model is designed to predict if a vehicle will experience a fault and time to
first fault, given a fault is predicted to occur.

###  4.1 Target Distribution Definition

The ultimate goal of this effort is to define the most accurate probability
distribution, shown in Equation (1), where $(z_{t},t\geq 0)$ is the fault
_target_ , or fault indicator time series, and $(x_{i,t},1\leq i\leq 38,t\geq
0)$ are the sensor _features_ or covariate time series. Here $t$ is defined as
a time point and $i$ is the label of sensor feature covariates, where there
are $38$ total sensor covariates in the CBM dataset considered. With this
target distribution, a vehicle operator or maintainer would know precisely the
probability of _if_ a vehicle is to experience a fault from time point $t_{0}$
to $t_{0}+T$, or $t=[t_{0}:t_{0}+T]$, and _when_ the fault would occur. We
propose this information will allow for vehicle operators and maintainers to
significantly improve predictive maintenance. Below we define this conditional
probability distribution, which we name the target distribution:

$P(z^{k}_{t_{0}:t_{0}+T}|z^{k}_{0:t_{0}-1},(x^{k}_{i,0:t_{0}-1},1\leq i\leq
38))\>\>\forall k\in K\>\>\forall\>t_{0}\in S,$ (1)

where, $z^{k}_{t_{0}:t_{0}+T}$ represents the time series of fault signals
collected from vehicle $k$ from time point $t_{0}$ to $t_{0}+T$,
$x^{k}_{i,t_{0}:t_{0}-1}$ represents the time series of sensor feature $i$
collected from vehicle $k$ from time point $0$ to $t_{0}-1$, and for model
training and experimental purposes, the set $S$ comprises randomly selected
start times for each vehicle $k$ in $K$, where $K$ is the set of selected
vehicles for analysis.

###  4.2 Generative Model Definition

While ideally we could formulate the probability distribution in Equation (1)
explicitly, in practice this is not feasible so we employ a generative model
to forecast outcomes given the data collected, and subsequently employ these
results to model the target distribution. To accomplish this, we employ a
hybrid of generative models, using DeepAR in conjunction with STAM (see
Section 2.1) and a VAE. The complete prediction process is described below:

  * •

Step 1: We train a DeepAR model $(Q_{\theta})$ on a large (relative to
forecast length) block of vehicle sensor feature and fault target data.

  * •

Step 2: Given time point $t_{0}$ to begin our forecasting period, we train the
STAM model $(g_{\theta})$ and a VAE $(v_{\theta})$ to generate estimates of
future sensor feature covariates.

  * •

Step 3: We generate future targets from times $t_{0}$ to $t_{0}+T$ using the
previously trained DeepAR model ($Q_{\theta}$), and covariates generated from
both the STAM model $(g_{\theta})$ and the VAE model $(v_{\theta})$.

  * •

Step 4: We use classifier ($C$ \- see Equation (3)) and regression model ($R$
\- see Equation (4)) to determine if there is a fault during the forecasted
time period $[t_{0}:t_{0}+T]$, and what the time to first fault is, given a
fault is expected to occur.

The complete process of model training and experimentation on a single vehicle
dataset is shown in Figure 2, with each step identified. We note that each
vehicle $k\in K$ requires the training of a DeepAR model, $Q_{\theta}$, and
each forecast $t_{0}\in S$ requires the training of a STAM model,
$g_{\theta}$, potentially with the evaluation of the VAE model, $v_{\theta}$.
We will further explain this process in the proceeding subsections.

Figure 2: Overview of generative model for conditional vehicle fault
prediction on vehicle $k\in K$. At Step 1, a DeepAR model ($Q_{\theta}$) is
trained using time series of fault target and sensor feature values from time
point $0$ to $q$, and is subsequently evaluated at randomly selected time
periods, with $t_{0}\in S$. At Step 2, a STAM model ($g_{\theta}$) is trained
over a period of $b=$ 3 hours and a VAE model ($v_{\theta}$) is trained on out
of sample data. At Step 3, the DeepAR model ($Q_{\theta}$) generates
representations of the fault indicator $\tilde{h}^{k}_{t_{0}:t_{0}+T}$ (resp.
$\hat{h}^{k}_{t_{0}:t_{0}+T}$) of $T=$ 30 minutes using the sensor feature
covariates generated from the STAM model (resp. the VAE model), where STAM
utilizes a lookback length of 20 minutes. At Step 4, the classifier ($C$) and
the regression model ($R$) use the predictions from both the STAM and VAE
model to forecast if there is a fault and what time the first fault is (given
a fault is predicted) from time point $t_{0}$ to $t_{0}+T$.

###  4.3 Spatio-Temporal Attention Mechanism (STAM) and Variational Auto
Encoder (VAE) for Feature Generation

As discussed in Step 1, in order to develop Equation (1), we reference the
DeepAR model, $Q_{\theta}$ and define a lookback parameter of length $\ell$
and Spatio-Temporal Attention Mechanism (STAM) training length $b$. It is
critical to note that to achieve reasonable performance with the DeepAR model,
_features_ must be available during the prediction time period
$t=[t_{0}:t_{0}+T]$, where $t_{0}$ is the prediction start time ($t_{0}\in S$
for each vehicle $k$ in $K$) and $t_{0}+T$ the end time for the selected
forecast. The requirement of features during the prediction time poses a
problem during practical forecasting because this information is _future
sensor_ data which is not available at prediction time. To meet this
requirement, we produce sensor feature data during the prediction time period
using two generative methods:

  1. 1.

$g_{\theta}\equiv\textit{spatio-temporal attention mechanism (STAM)}$, where
the STAM is trained on recent data:
$g_{\theta}(x^{k}_{i,t_{0}-b:t_{0}-1})\Rightarrow\tilde{x}^{k}_{i,t_{0}:t_{0}+T}$.
This model captures the vehicle’s most recent operational history.

  2. 2.

$v_{\theta}\equiv\textit{variational auto-encoder (VAE)}$, where the VAE is
trained on out of sample data:
$v_{\theta}(\hat{x}^{k^{\prime}}_{i,t_{0}^{\prime}:t_{0}^{\prime}+T})\Rightarrow\hat{x}^{k}_{i,t_{0}:t_{0}+T}$
such that $t_{0}^{\prime}\in S^{\prime}$ and $k^{\prime}\in K^{\prime}$, and
$k^{\prime}\notin K$. This model captures how the vehicle will be operated in
the future, which is either based upon the true data,
$v_{\theta}(x^{k}_{i,t_{0}:t_{0}+T})$, or conditionally selected based on
latent space cluster state,
$v_{\theta}(\hat{x}^{k^{\prime}}_{i,t_{0}^{\prime}:t_{0}^{\prime}+T})$, to
_simulate_ actual vehicle operating conditions.

In the prediction process, with the generated sensor data from both the STAM
and VAE generators, we propose calculating two _representations_ of the fault
indicator future time series:

$Q_{\theta}(z^{k}_{t_{0}:t_{0}+T}|z^{k}_{0:q},x^{k}_{i,0:q},\tilde{x}^{k}_{i,t_{0}:t_{0}+T})\Rightarrow\tilde{h}^{k}_{t_{0}:t_{0}+T}$
(2a)
$Q_{\theta}(z^{k}_{t_{0}:t_{0}+T}|z^{k}_{0:q},x^{k}_{i,0:q},\hat{x}^{k}_{i,t_{0}:t_{0}+T})\Rightarrow\hat{h}^{k}_{t_{0}:t_{0}+T}$
(2b)

where $\tilde{h}^{k}_{t_{0}:t_{0}+T}$ and $\hat{h}^{k}_{t_{0}:t_{0}+T}$, are
hidden representations of the target distribution of $z^{k}_{t_{0}:t_{0}+T}$.
We call these values hidden representations because they are used in the
classification and regression models,
$C(\tilde{h}^{k}_{t_{0}:t_{0}+T},\hat{h}^{k}_{t_{0}:t_{0}+T})=c\in\\{0,1\\}$
and
$R(\tilde{h}^{k}_{t_{0}:t_{0}+T}|C(\tilde{h}^{k}_{t_{0}:T},\hat{h}^{k}_{t_{0}:t_{0}+T})=1)=f\in\mathbb{R}^{+}$,
shown in Equations (3) and (4), where the regression model $R$ is not trained
with the VAE generated hidden representation, $\hat{h}^{k}_{t_{0}:t_{0}+T}$.
Here the predicted value from the regression, $f$, is defined as the time to
_first_ fault observed over the forecasted time period. We focus on predicting
the time to first fault, as it has the most practical implications with
respect to vehicle condition and function: the first fault is most likely to
cause a maintenance event affecting operation. The classification and
regression models take in the hidden representation data given each type of
input and predicts fault or no-fault, and time to first fault. We employ the
random forest classification and random forest regression models for $C$ and
$R$, respectively. These models were selected based upon their performance and
potential for interpretability [Ho (1995)]. Additionally, the parameters of
the DeepAR model, $\theta$, are identical for the evaluation of $Q$ when used
with generated covariates $\tilde{x}^{k}_{i,t_{0}:t_{0}+T}$ and
$\hat{x}^{k}_{i,t_{0}:t_{0}+T}$.

$C(\tilde{h}^{k}_{t_{0}:t_{0}+T},\hat{h}^{k}_{t_{0}:t_{0}+T})=c\in\\{0,1\\}.$
(3)

$R(\tilde{h}^{k}_{t_{0}:t_{0}+T}|C(\tilde{h}^{k}_{t_{0}:t_{0}+T},\hat{h}^{k}_{t_{0}:t_{0}+T})=1)=f\in\mathbb{R}^{+}.$
(4)

We note that when training the classification model $C$, $v_{\theta}$ is given
the _true_ future covariates for encoding and decoding, using the trained VAE:
$v_{\theta}(x^{k}_{i,t_{0}:t_{0}+T})\Rightarrow\hat{x}^{k}_{i,t_{0}:t_{0}+T}$.
This is to ensure an accurate fault classification is given when conditioning
on selected latent cluster states. While in practice the _true_ future sensor
covariate timer series, $x^{k}_{i,t_{0}:t_{0}+T}$, would not be available for
the forecast, the information which characterizes the vehicle latent cluster
_state_ would be available. Using the VAE model, we would be able to generate
future fault covariates from the true vehicle state, where covariates from the
common states are known to be close in the latent space.

Next, in Section 4.4, we will discuss the random selection and generation of
sensor data from the trained VAE latent space, and subsequent simulation of
fault probabilities conditioned on this state. Finally, we observe that with
the output of Equations (3) and (4), we may produce a prediction for
${z}^{\star k}_{t_{0}:t_{0}+T}$, or a final time series forecast indicating
exactly when an _initial_ or first fault will occur. We propose that this is
critically important information for predictive maintenance.

###  4.4 Simulation using Generative Modeling Conditioned on Vehicle State

As illustrated in Figure 3, the training of a VAE, $v_{\theta}$, on _out of
sample_ vehicle data
$v_{\theta}(\hat{x}^{k^{\prime}}_{i,t_{0}^{\prime}:t_{0}^{\prime}+T})\Rightarrow\hat{x}^{k}_{i,t_{0}:t_{0}+T}$
such that $t_{0}^{\prime}\in S^{\prime}$ and $k^{\prime}\in K^{\prime}$, and
$k^{\prime}\notin K$ allows us to define latent states, in a reduced
dimension, which represent practical considerations associated with vehicles.
These practical considerations could include: vehicle age, maintenance
condition, location, etc. We enumerate the practical considerations present in
in the CBM dataset and discuss their impact on the vehicle latent states in
Section 5.1. These states are designated by clustering on the reduced latent
dimension produced by the encoding step of the VAE.

Figure 3: Illustration of mapping vehicle operational condition to generate
future sensor readings.

In order to simulate vehicle conditions, which were not observed in the actual
data, we may randomly sample data points from the selected state, decode the
random point to generate sensor data,
$\hat{x}^{k^{\prime}}_{i,t_{0}^{\prime}:t_{0}^{\prime}+T}$. We may then use
this data in the DeepAR model and generate the hidden representation
$\hat{h}^{k^{\prime}}_{t_{0}^{\prime}:t_{0}^{\prime}+T}$, with the simulated
sensor values. This will produce alternate values for the classification of
vehicle faults and time to first fault, using models $C$ and $R$ respectively,
based on practical considerations.

##  5 Results

We develop our results addressing two topics. First, what does the clustering
of vehicle sensor data, discussed in Section 4.4, practically indicate with
respect to vehicle conditions. Second, how accurate is the discriminative
model described in Section 4.1 at predicting if there will be a vehicle fault
observed during the generated time period, and how accurately can we predict
when the first fault will occur. We provide results allowing for model
interpretation, and comparison to alternate generative model formulations. We
discuss the interpretation of model parameters and convergence analysis
regarding the training of the DeepAR model in the appendix.

###  5.1 Simulation using Latent Space Representation and K-means Clustering

In Figure 4 we provide a representation of vehicle out of sample sensor data
in a three dimensional latent space. Here the trained VAE, $v_{\theta}$, maps
historical sensor readings
$\hat{x}^{k^{\prime}}_{i,t_{0}^{\prime}:t_{0}^{\prime}+T}$, such that
$t_{0}^{\prime}\in S^{\prime}$ and $k^{\prime}\in K^{\prime}$, into the latent
representation shown. See Section 4 for definitions of these variables.
Subsequently, the K-means clustering algorithm is applied to the latent
representation, producing five clusters.

(a) K-means clusters with state interpretation.

(b) VAE latent representation with metadata.

Figure 4: K-means cluster representation labeled on VAE latent space with
metadata interpretation.

We referenced metadata with respect to vehicle location, odometer mileage,
engine hours, and sub-family and found that the VAE latent state mapping and
K-means application identified practical clusters or _states_ , numbered 0-4,
in which each vehicle operates. Figure 4(b) provides the latent representation
of each of the metadata categories (location, vehicle odometer mileage, engine
hours, and vehicle sub-family). This figure shows how the logical clusters or
states were interpreted in Figure 4(a). In Figure 5 we provide another view of
the three dimensional latent space, where each point is labeled according to
whether the vehicle truly experienced a fault, or did not experience a fault,
during the generated time period $t=[t_{0}:t_{0}+T]$.

Figure 5: Latent space fault visualization.

As proposed in Section 4, we complete a simulation generating sensor data
using the decoding function of the VAE to condition the trained DeepAR model
$Q_{\theta}$ when producing $\hat{h}^{k}_{t_{0}:t_{0}+T}$, the VAE generated
hidden representation of the true fault indicator function
$z^{k}_{t_{0}:t_{0}+T}$. We employ the trained classifier, $C$, proposed in
Equation (3) to determine if a vehicle experiences a fault during a predicted
time period. This allows practitioners to generate simulated outcomes
conditioned on practical considerations, or “what if” scenarios associated
with their vehicles: i.e. what would happen if they changed the location the
vehicle operates out of, what happens as the vehicle mileage, or engine age
increase, etc. In Table 5.1 we provide these simulated probabilities. In this
simulation, for the same test vehicles and randomly selected time periods used
throughout this analysis, we have generated sensor covariates,
$\hat{x}^{k}_{t_{0}:t_{0}+T}$, by randomly sampling data points in the VAE
latent space from a fixed state, where the members of the state are determined
by the K-means model. For example, if all the vehicles ($k\in K$) are assumed
to be operating in State 0, we may sample hundreds of points in the latent
space only from this fixed state, generate $\hat{x}^{k}_{t_{0}:t_{0}+T}$ using
the trained VAE $v_{\theta}$, and subsequently use this generated data in the
DeepAR model $Q_{\theta}$ to produce the hidden representation
$\hat{h}^{k}_{t_{0}:t_{0}+T}$, which is then evaluated using the trained
classifier $C$. In this example for State 0, we expect a fault rate of
$0.080$.

Table 1: Simulation of fault rates using sensor data generated from vehicle states. State | Conditional Fault Rate
---|---
State 0 | 0.080
State 1 | 0.078
State 2 | 0.069
State 3 | 0.111
State 4 | 0.098

The conditioned fault rates observed in this simulation follow logically given
the practical interpretation of vehicle states given in Figure 4(a). We
observe, relative to other states, that the fault rate of States 3 and 4 are
the highest, where these are vehicles with medium to high engine hours and
lower odometer readings. This may be an indicator of scheduled service
requirements for newer vehicles. Additionally, considerations with locations
$1$ and $3$ may weigh on the fault rates as different locations often have
varying maintenance programs. As noted in the introduction and formulation of
this generative learning model, a great advantage is given by the ability to
visualize and interpret results in order to understand why the model is
generating specific fault forecasts. Furthermore, the values provided in Table
5.1 are validated by US Army published standards and observed operational
readiness rates of vehicles, cited to be between $80\%-100\%$ [Kays, James, L
and Carlton, William, B and Lee, Mark, M and Ratliff, William, L (1998)].

###  5.2 Fault Prediction

In this section, we provide a visual representation of the performance of the
discriminative classifier proposed in Equation (3), which determines if a
fault is going to be forecasted during the generated time period
$t=[t_{0}:t_{0}+T]$. In the models trained and experiments visualized, the
time period of prediction is 30 minutes, where $T=1800\>secs$. To be specific,
we consider two cases: 1. using only the STAM generated covariate data, and 2.
both the STAM and the VAE generated covariate data. In the latter case, for
the VAE we use the actual latent state of the vehicle, by referencing the true
future feature covariates. Using the actual latent state facilitates the
development of the simulation by mapping the correct latent state to the
hidden representation. The latent state is available in practice as the
practitioner knows current attributes of the vehicle (vehicle mileage, engine
age, etc) and may choose the appropriate state.

For each case, we choose 50 pairs of training and test sets through random
permutation of the data. We note that the training and test sets were not
split between distinct vehicles, because in practice the regression model will
be updated as the vehicle is operated. As shown in Figure 6, for the test set
in each of the 50 pairs, we plot the Receiver Operating Characteristics (ROC)
curve using a unique color. The average Area Under the Curve (AUC) for the 50
tests are also computed.

(a) ROC curves using STAM covariates.

(b) ROC curves using STAM and VAE covariates.

Figure 6: ROC curves of our proposed discriminative classifiers. We choose 50
pairs of training and test data by random permutation. We plot the ROC curve
on each test data using different colors after training our classifiers on the
corresponding training data. We compute the average AUC for the 50 ROC curves.
In Figure 6a, we use the representations of the fault indicator generated by
the DeepAR model that takes in only the STAM generated sensor covariates,
while in Figure 6b, we utilize the representations of the fault indicator
generated by the DeepAR model that takes in both the STAM and VAE generated
sensor covariates.

These results articulate the highly accurate and steady performance achieved
when using: 1. only the STAM generated covariate data, resulting in a $0.950$
AUC, and 2. both the STAM and the VAE generated covariate data, resulting in a
$0.978$ AUC. More specifically, Figure 6(a) displays the ROC curves when the
hidden representation ($\tilde{h}^{k}_{t_{0}:t_{0}+T}$) is generated by the
DeepAR model taking in only the STAM generated covariates (see Equation (2a))
and this hidden representation is used in the training and testing of the
classifier: $C(\tilde{h}^{k}_{t_{0}:t_{0}+T})$. In Figure 6(b), the hidden
representations ($\tilde{h}^{k}_{t_{0}:t_{0}+T}$ and
$\hat{h}^{k}_{t_{0}:t_{0}+T}$) generated by the DeepAR models using the
covariates generated from the STAM and the VAE (see both Equations (2a) and
(2b)) are both used in the training and testing of the classifier:
$C(\tilde{h}^{k}_{t_{0}:t_{0}+T},\hat{h}^{k}_{t_{0}:t_{0}+T})$. We observe
that the inclusion of both the STAM and VAE generated covariates improves the
average AUC performance by $2.8\%$. This improvement is expected as additional
information is provided to the DeepAR ($Q_{\theta}$) model via the VAE
($v_{\theta}$), with the actual latent state information. Even without the
latent state information, the model performs with satisfactory accuracy
($0.950$). To further illustrate the advantages of our method in fault
prediction, we compare the AUC performance of our method to that of other
existing methods in Section 5.4.

###  5.3 Time to First Fault Prediction

As discussed in Section 4.2, once a fault is determined to occur during the
predicted time period,
$C(\tilde{h}^{k}_{t_{0}:t_{0}+T},\hat{h}^{k}_{t_{0}:t_{0}+T})=1$, we employ a
random forest regression model,
$R(\tilde{h}^{k}_{t_{0}:t_{0}+T}|C(\tilde{h}^{k}_{t_{0}:t_{0}+T},\hat{h}^{k}_{t_{0}:t_{0}+T})=1)=f\in\mathbb{R}^{+}$,
to predict the time to first fault [Ho (1995)]. The regression model $R$
achieves a coefficient of determination ($r^{2}$) of $0.77$ in predicting time
to first fault, where $R$ was trained on the 50th quantile generated results.
We define the coefficient of determination as
$r^{2}=1-\frac{\sum_{i}(y_{i}-f_{i})^{2}}{\sum_{i}(y_{i}-\bar{y})^{2}}$, where
$y_{i}$ is the true time to first fault, $f_{i}$ is the forecasted time to
first fault for observation $i$, and $\bar{y}$ is the true mean. We note that
all true time to first fault values are no further than 4 minutes
($240\>secs$) from the forecasted time to first fault; however, these
predictions could benefit from further refinement on robustness since the
present prediction indicates potential under-estimation of the time to first
fault. This issue has practical implications for predictive maintenance and we
seek to address this shortcoming in future work.

###  5.4 Comparison to Alternate Models on Fault Prediction

We selected two methods for comparison to the the proposed model. The selected
methods include a LSTM and a STAM model, where both are used to directly
generate hidden representations of the target fault indicator. By contrast,
our proposed method uses the STAM model to generate sensor covariates, and
then inputs the generated sensor covariates into the DeepAR model to generate
hidden representations of the target fault indicator.

The time period is set to be $t=[t_{0}:t_{0}+T]$. As with the proposed
methodology, the directly generated hidden representation, denoted as
$\tilde{h}^{k}_{t_{0}:t_{0}+T}$, serves as the input into a classification
model, $C(\tilde{h}^{k}_{t_{0}:t_{0}+T})=c\in\\{0,1\\}$. The performance of
the comparison models are visualized through plots of ROC curves in Figure 7.
In Table 2, we compare the average AUC for these methods and see clearly that
the performance for both the LSTM and STAM models is significantly poorer when
compared to the proposed model.

(a) ROC curves of a LSTM-based classifier.

(b) ROC curves of a STAM-based classifier.

Figure 7: ROC curves of alternate classification models. We choose 50 pairs of training and test data by random permutation. We plot the ROC curve on each test data using different colors after training our classifiers on the corresponding training data. We compute the average AUC for the 50 ROC curves. In Figure 7a, we we utilize the representations of the fault indicator directly generated by LSTM, while in Figure 7b, we we utilize the representations of the fault indicator directly generated by STAM. Table 2: Model comparison table. Model for Comparison | AUC
---|---
LSTM Baseline | 0.520
STAM Baseline | 0.458
DeepAR with STAM Covariates (ours) | 0.950
DeepAR with STAM and VAE Covariates (ours) | 0.978

###  5.5 Discussion

A clear question arises from the proposed model: why does the use of the STAM
method to generate covariates, and subsequently generating target values using
the DeepAR model produce such an accurate fault representation? This is
clearly true when compared to other baseline methods as presented in Figure 7.
We believe this is due to several attributes of both the DeepAR and STAM
models, namely: (i) DeepAR, being auto-regressive and recurrent, is
particularly designed to learn seasonal behavior well, which is specifically
helpful with larger and longer time series datasets (ii) DeepAR incorporates
the use of covariates focused towards learning a specific target distribution,
(iii) STAM is demonstrated to be particularly scalable and accurate when
forecasting time series values with multiple variables.

Specifically referencing point (i), we know that the vehicle sensor feature
and fault target data exhibits seasonality or cyclic behavior when explored
over significant periods of time. As such, we trained the DeepAR over 36 hours
of data (129600 data points). This allowed for the auto-regressive network to
learn sufficient detail of recurrent structure over the fault targets and
sensor feature covariates, as specified in point (ii). Additionally, while
alone not accurate when predicting target values as shown in Section 5.4, the
STAM model is excellent at generating the covariates when trained over modest
time periods, such as 3 hours of data in the proposed model (10800 data
points). STAM provides sufficient accuracy when forecasting the subsequent 30
minutes (1800 data points) of sensor feature data to produce fault target data
which is sufficiently distinguishable to support highly accurate binary
classification and regression results. Attributes of the STAM and DeepAR
models are combined with appropriate training periods to leverage the
strengths of both methods in order to facilitate forecasts that are well
suited for this vehicle fault prediction problem.

##  6 Conclusion and Further Work

We proposed a modeling approach which allows for the highly accurate
prediction of future vehicle faults, and the simulation of future faults
conditioned on practical considerations. The simulated results are validated
using logical interpretation of vehicle metadata and published documentation.
We demonstrated that the model is interpretable and provides actionable
results for practitioners in the vehicle operation and maintenance fields.
Significant further work remains including:

  1. 1.

Improving robustness of time to first fault forecast (Section 5.3) with
mitigation of underestimation.

  2. 2.

Extending prediction time to the scale of hours.

  3. 3.

Simulating the evolution of states for multiple vehicles to streamline fleet
maintenance scheduling.

We note that a link to a github repository with all modeling code and
processor information is available in Appendix B. The work presented
represents an initial step in this modeling effort, and we will continue to
update our procedures to produce accurate and interpretable results for the
predictive maintenance effort.

##  7 Acknowledgement

The material in this paper is based upon work supported by the Air Force
Office of Scientific Research under award number FA9550-20-1-0397.

##  Appendix A Model Interpretability: Attention Weight Analysis

An important attribute of the proposed analysis, given the methods used, is
the models provide interpretable parameters. These parameters are produced
during the training of the STAM neural network, used to generate future sensor
feature covariates attn_mech. The STAM produces attention weights when
training, where these weights may be interpreted as an importance metric for
each of the covariates or sensors, when generating future values. In Figure 8,
we visualize a mean of recorded attention weights, averaged when executing the
experiments presented in Sections 5.2 and 5.3, over the randomly sampled time
periods for all $k\in K$ and $t_{0}\in S$. Two sets of attention weights are
visualized: one set over the samples where a fault is observed during the time
period, and another where a fault is not observed.

Figure 8: Attention weight analysis

We observe a shift in attention weights when there is a fault observed.
Specifically, the sensors associated with the transmission subsystem become
less important and the engine subsystem becomes more important. These results
follow logically, as the engine subsystem is a more complicated component of
the vehicle with numerous components and potential points of fault. These
results are shown in Table 3.

Table 3: Attention weight comparison table | Attention Weight
---|---
Engine Without Fault | 0.480
Transmission Without Fault | 0.520
Engine With Fault | 0.515
Transmission With Fault | 0.485

##  Appendix B Implementation Code

The code base, written in python, is available at the following github link:
generative predictive maintenance github code. There is only one vehicle
dataset provided for demonstration purposes at this github repository. The
full dataset was not provided due to file size considerations, and maybe
referenced for reproducibility by contacting the following email: author
email. The analysis referenced in the paper was completed with 200 vehicles of
data. Training of the STAM and VAE models were completed on NVIDIA GeForce RTX
3090 GPUs with CUDA Version 12.2. Further details may be referenced in the
linked github repository.

##  Appendix C Model Training Convergence Analysis

We provide a visualization of the loss history recorded during the training of
the DeepAR model. This allows us to confirm the DeepAR training convergence
visually. Given that we trained a single DeepAR model for each vehicle
analyzed in the CBM dataset, we have many training loss profiles. We see
clearly the convergence of the training, for all DeepAR models, over the 100
epoch training period.

Figure 9: Training loss profile

## References

  * Box and Jenkins (1976) Box, G. and G. M. Jenkins. 1976.  Time Series Analysis: Forecasting and Control.  Holden-Day.
  * CBO (2021) CBO 2021.  “Projected Acquisition Costs for the Army’s Ground Combat Vehicles”.
  * DoDIG (2023) DoDIG 2023.  “Management Advisory: Maintenance Concerns for the Army’s Prepositioned Stock–5 Equipment Designated for Ukraine”.  DODIG Report No. DODIG-2023-076.
  * Gangopadhyay et al. (2020) Gangopadhyay, T., S. Y. Tan, Z. Jiang, R. Meng and S. Sarkar. 2020.  “Spatiotemporal Attention for Multivariate Time Series Prediction and Interpretation”.  CoRR abs/2008.04882.
  * GAO (2022) GAO 2022.  “Military Readiness: Actions Needed to Further Predictive Maintenance on Weapon Sytems”.  GAO Report No. GAO-23-105556.
  * Ho (1995) Ho, T. K. 1995.  “Random decision forests”.  In Proceedings of 3rd International Conference on Document Analysis and Recognition, Volume 1, 278–282 vol.1 https://doi.org/10.1109/ICDAR.1995.598994.
  * Hojjati et al. (2023) Hojjati, H., M. Sadeghi, and N. Armanfard. 2023.  “Multivariate Time-Series Anomaly Detection with Temporal Self-supervision and Graphs: Application to Vehicle Failure Prediction”.  In Machine Learning and Knowledge Discovery in Databases: Applied Data Science and Demo Track, edited by G. De Francisci Morales, C. Perlich, N. Ruchansky, N. Kourtellis, E. Baralis, and F. Bonchi, 242–259.  Cham: Springer Nature Switzerland.
  * Kays, James, L and Carlton, William, B and Lee, Mark, M and Ratliff, William, L (1998) Kays, James, L and Carlton, William, B and Lee, Mark, M and Ratliff, William, L 1998.  “Analysis of Operational Readiness Rates”.  Defense Technical Information Center.
  * Salinas et al. (2019) Salinas, D., V. Flunkert, and J. Gasthaus. 2019.  “DeepAR: Probabilistic Forecasting with Autoregressive Recurrent Networks”.
  * Shafi et al. (2018) Shafi, U., A. Safi, A. R. Shahid, S. Ziauddin and M. Q. Saleem. 2018.  “Vehicle remote health monitoring and prognostic maintenance system”.  Journal of advanced transportation 2018:1–10.
  * Theissler et al. (2021) Theissler, A., J. Pérez-Velázquez, M. Kettelgerdes, and G. Elger. 2021.  “Predictive maintenance enabled by machine learning: Use cases and challenges in the automotive industry”.  Reliability Engineering & System Safety 215:107864 https://doi.org/https://doi.org/10.1016/j.ress.2021.107864.
  * Thurston, Michael, G and McConky, Sean, P and Valant, Chris, J and Nenadic, Nenad, G (2019) Thurston, Michael, G and McConky, Sean, P and Valant, Chris, J and Nenadic, Nenad, G 2019.  “Feasibility of Diagnostics, Prognostics and Hybrid Prognostics across Multiple Platforms”.  Defense Technical Information Center.

## AUTHOR BIOGRAPHIES

PATRICK KUIPER received his B.S. in Operations Research from The United States
Military Academy at West Point in 2007, and a M.Eng. in Applied Mathematics
from Harvard University in 2014. He is an active duty officer in the United
States Army and served as an instructor in the West Point Department of
Mathematical Sciences from 2016 to 2019. He is currently pursuing a Ph.D. in
Electrical and Computer Engineering at Duke University. His research interests
include applications of generative modeling and extreme value theory. His
email address is .

SIRUI LIN is a Ph.D. student at the Department of Management Science and
Engineering (MS&E) at Stanford University, with research interest in applied
probability and optimal transport-based statistical inference. His email
address is .

JOSE BLANCHET is a Professor of Management Science and Engineering (MS&E) at
Stanford. Prior to joining Stanford, he was a professor at Columbia
(Industrial Engineering and Operations Research, and Statistics, 2008-2017),
and before that he taught at Harvard (Statistics, 2004-2008). Jose is a
recipient of the 2010 Erlang Prize and several best publication awards in
areas such as applied probability, simulation, operations management, and
revenue management. He also received a Presidential Early Career Award for
Scientists and Engineers in 2010. He is the Area Editor of Stochastic Models
in Mathematics of Operations Research and has served on the editorial board of
Advances in Applied Probability, Bernoulli, Extremes, Insurance: Mathematics
and Economics, Journal of Applied Probability, Queueing Systems: Theory and
Applications, and Stochastic Systems, among others. His email address is .

VAHID TAROKH worked at the AT&T Labs-Research until 2000. From 2000 to 2002 he
was an Associate Professor at Massachusetts Institute of Technology (MIT). In
2002 he joined Harvard University as a Hammond Vinton Hayes Senior Fellow of
Electrical Engineering and Perkins Professor in Applied Mathematics. He joined
Duke University in 2018 as the Rhodes Family Professor in Electrical and
Computer Engineering, Computer Science, and Mathematics and Bass Connections
Endowed Professor. He was also a Gordon Moore Distinguished Research Fellow at
CALTECH in 2018. Since Jan 2019, he has also been named as a Microsoft Data
Science Investigator at Duke University. His email address is .


# Solving systems of linear equations through zero forcing set and application
to lights out

Jianbo Wang111<EMAIL_ADDRESS>University of Electronic Science and
Technology of China, School of Computer Science and Engineering. Chao Xu222
<EMAIL_ADDRESS>University of Electronic Science and Technology of China,
School of Computer Science and Engineering. Siyun Zhou333
<EMAIL_ADDRESS>University of Electronic Science and Technology of
China, School of Mathematical Sciences.

###### Abstract

Let $\mdmathbb{F}$ be any field, we consider solving $Ax=b$ repeatedly for a
matrix $A\in\mdmathbb{F}^{n\times n}$ of $m$ non-zero elements, and multiple
different $b\in\mdmathbb{F}^{n}$. If we are given a zero forcing set of $A$ of
size $k$, we can then build a data structure in $O(mk)$ time, such that each
instance of $Ax=b$ can be solved in $O(k^{2}+m)$ time. As an application, we
show how the lights out game in an $n\times n$ grid is solved in $O(n^{3})$
time, and then improve the running time to $O(n\log n)$ by exploiting the
repeated structure in grids.

##  1 Introduction

Solving the linear system $Ax=b$ has been a fundamental problem in
mathematics. We consider the problem in its full generality. Given a field
$\mdmathbb{F}$, a matrix $A\in\mdmathbb{F}^{n\times n}$ and a vector
$b\in\mdmathbb{F}^{n}$, we are interested to find an $x\in\mdmathbb{F}^{n}$
such that $Ax=b$. Additionally, a zero forcing set associated with $A$ is also
given.

As arguably the most important algorithmic problem in linear algebra, there
are many algorithms designed for solving systems of linear equations, we refer
the reader to [12]. However, most works are for matrices over the real,
complex or rational number fields. As this work is concerned with general
fields, we will describe some known general algorithms below. The Gaussian
elimination is a classic method for solving systems of linear equations, which
requires $O(n^{3})$ time. If the system is of full rank, one can reduce the
problem to a single matrix multiplication, by taking $x=A^{-1}b$. In general
cases, the more recent LSP decomposition, which is a generalization of the LUP
decomposition, allows the running time to match the matrix multiplication [15,
17]. If one views the non-zeros of the matrix as an adjacency matrix of a
graph, the property of the graph can be then used to speed up the algorithm.
For example, if the graph is planar, then the nested-dissection technique can
obtain a $O(n^{3/2})$ time algorithm when the field is real or complex [20].
Later, it was generalized to any non-singular matrix over arbitrary field, and
the running time was also improved to $O(n^{\omega/2})$ [4], where
$\omega<2.3728596$ is the matrix multiplication constant [3]. These algorithms
are fairly complicated, but share one similarity with our work in taking the
advantage of properties of the graph to improve the running time of the
algorithm.

Zero forcing sets were first studied in [2] on graphs relating to the maximum
nullity of a matrix, and the results were later expanded to directed graphs
[5]. The zero forcing set captures some "core" information of the linear
system: $Ax=b$ is uniquely determined by the values of $x$ on a zero forcing
set. Therefore this implies often that one just have to observe part of $x$ to
recover the entirety of $x$. This idea leads to independent discovery of zero
forcing set by physicists for control of quantum systems [8, 21]. Later it was
shown to be applicable to Phasor Measurement Units (PMU) to monitor power
networks [9]. Most studies of zero forcing set are on the algebraic and
combinatorial properties. On the computational front, finding the smallest
zero forcing set of an undirected graph is not known to be solvable in
polynomial time, and it is likely NP-hard. The related problem of finding the
minimum 0-1-weighted zero forcing set in an undirected graph is hard [1]. For
directed graphs, finding the zero forcing set is NP-hard [22]. Exact
algorithms for finding the smallest zero forcing set was considered, see [7]
for a survey.

This work was inspired by an algorithm for the lights out game. In the game,
there is a light on each vertex of the graph, and each light is in a state
either on ($1$) or off ($0$). There is also a button on each vertex. Pressing
the button would flip the state of the light of itself and all its
(out-)neighbors. The goal is to turn off all the lights. The lights out game
is equivalent to solving a system of linear equations $Ax+b=0$ in
$\mdmathbb{F}_{2}$, where $A$ is the adjacency matrix of the graph $G$, and
$b$ is the state of the lights and $\mdmathbb{F}_{2}$ is the finite field of
order $2$. The interpretation is that $x_{v}=1$ means pressing the button at
vertex $v$, and $b_{v}$ is the initial state of the light at vertex $v$. Note
in $\mdmathbb{F}_{2}$, $-b=b$, hence it is equivalent to solving $Ax=b$. There
is a large amount of research in the lights out game, see [10] for a survey.
Finding a solution of the lights out game with the minimum number of button
presses is NP-hard [6].

We focus on the case when $G$ is an $n\times n$ grid, which corresponds to an
$n^{2}\times n^{2}$ matrix. Gaussian elimination would take $O(n^{6})$ time.
An alternative approach is the light-chasing algorithm [19]. Since the lights
in the first row uniquely affect the states of the remaining rows, one can
then look up which action on the first row should be taken according to the
states of the last row. However, the literature does not provide the strategy
for finding the lookup table. Wang claimed there is a $O(n^{3})$ time
algorithm that given the state of the last row, the state of the first row can
be found through solving a linear equation of an $n$-square matrix [23].
Wang’s result is the motivation behind this work. However, there is no proof
of its correctness. If the adjacency matrix has full rank, then by nested
dissection, there is a randomized Las Vegas algorithm with running time of
$O(n)$ [4]. Unfortunately, there are infinite many singular adjacency matrices
[14], the list of such matrices can be found in OEIS A117870 [16]. A SAGE
package for computing the solution of the $n\times n$ grid lights out game
when all the lights are on is provided in [11].

#### Our contribution

We generalize Wang’s method for the lights out game to solving an arbitrary
system of linear equations $Ax=b$, and give a formal proof of its correctness.
We define a structure, the core matrix $B$ of $A$, such that one can solve the
system of linear equations over $B$ instead of $A$, and then lift it to a
solution of $A$ in linear time. As a consequence, we obtain the following
algorithmic result.

######  Theorem 1.1 (Main).

Given a matrix $A\in\mdmathbb{F}^{n\times n}$ of $m$ non-zero elements, and a
zero forcing set of $A$ of size $k$ where $k\leq n$. A data structure of size
$O(k^{2})$ can be computed in $O(mk)$ time, such that for each
$b\in\mdmathbb{F}^{n}$, solving $Ax=b$ takes $O(k^{2}+m)$ time.

We also show that the lights out game on an $n\times n$ grid can be solved in
$O(n\log n)$ time by finding the core matrix using an alternative method.
Additionally, we prove some linear algebraic properties of the zero forcing
set, which might be of independent interest.

##  2 Preliminaries

Define an index set $[n]=\\{1,\ldots,n\\}$. We consider an algebraic model,
where every field operation takes $O(1)$ time, and each element in the field
can be stored in $O(1)$ space.

It is useful to think about a vector $x$ indexed by elements $X$ is a function
$x:X\to\mdmathbb{F}$. Hence, we write $x\in\mdmathbb{F}^{X}$. Define
$\operatorname{supp}(x)$ to be the set of non-zero coordinates of $x$.

For a matrix $A$ and a set of row indices $R$ and a set of column indices $C$,
we define $A_{R,C}$ to be the submatrix of elements indexed by the row and
columns in $R$ and $C$, respectively. We use $*$ to mean the set of all row
indices or column indices, depending on context. For example, $A_{R,*}$ to
denote the submatrix of $A$, which is composed of the rows indexed by $R$. We
define $A_{i,*}$ as the $i$th row vector of $A$. Similarly, $A_{*,j}$ is the
$j$th column vector. $A_{i,j}$ is the element in the $i$th row and $j$th
column. For a vector $x$, $x_{i}$ is the scalar at index $i$, and $x_{I}$ is
the subvector of the elements in the index set $I$.

Solving $Ax=b$ is finding a vector $x$ that satisfies $Ax=b$. For a directed
graph $G=(V,E)$ on $n$ vertices, and a field $\mdmathbb{F}$, we define
$\mathcal{M}(G,\mdmathbb{F})$ to be the set of all matrices
$A\in\mdmathbb{F}^{n\times n}$ such that for all $u\neq v$, $A_{u,v}\neq 0$ if
and only if $(u,v)\in E$. We do not have any restriction on $A_{v,v}$. Let
$N^{+}(u)=\\{v\mid(u,v)\in E\\}$ be the set of all out-neighbors of $u$.

We briefly review some basic concepts and results related to the zero forcing.
Consider the process of coloring a graph. We start with a set of blue colored
vertices, denoted by $Z$. If there exists a blue vertex $v$ with exactly one
non-blue out-neighbor, say $u\in V$, then color $u$ blue. The operation of
turning $u$ blue is called _forcing_ , and we say $u$ is forced by $v$. If
this process ends up with a situation where all vertices are colored blue,
then we say the initial set $Z$ is a _zero forcing set_. The _zero forcing
number_ $Z(G)$ of $G$ is the size of the smallest zero forcing set. We also
call $Z$ a zero forcing set of $A$ if $A\in\mathcal{M}(G,\mdmathbb{F})$ for
some field $\mdmathbb{F}$. The following proposition gives the reason for the
name zero forcing.

######  Proposition 2.1 ([2, 13]).

For a zero forcing set $Z$ of $A$, if $x\in\ker(A)$, then $x_{Z}=0$ implies
$x=0$. Namely, $x$ vanishing at zero forcing set forces $x$ to be $0$.

The converse is not true in general. Indeed, we can simply take $A$ to be a
$2\times 2$ identity matrix. Then, $Z=\\{1,2\\}$ is the unique zero forcing
set of $A$, and $\ker(A)=\\{0\\}$. If $x\in\ker(A)$ and $x_{1}=0$, then we
have $x=0$. But $\\{1\\}$ is clearly not a zero forcing set.

The order of forces for a zero forcing set $Z$ is not unique, although the
final coloring is [2, 13]. For simplicity, we avoid the issue by considering a
particular chronological list of forces , which picks the smallest indexed
vertex that can be forced. is a total ordering of the vertices such that the
elements in $Z$ are ordered arbitrarily, and smaller than all elements in
$V\setminus Z$. For each $v,u\in V\setminus Z$, $v\leq u$ if $v$ is forced no
later than $u$. The forcing graph is a graph where there is an edge between
$u,v$ if $u$ forces $v$. It is well known that such graph is a set of node
disjoint paths that begin in $Z$ [13]. Hence, if $v$ forces $u$, we can define
$u^{\uparrow}=v$ to be the _forcing parent_ , and correspondingly, $u$ is
called the _forcing child_ of $u^{\uparrow}$. A vertex is a _terminal_ , if it
does not have a forcing child. Let $T$ be the set of terminals, then
$|T|=|Z|$.

Consider the algorithm $\textsc{Forcing}(A,Z,b)$ as summarized in Figure 2.1,
which takes a matrix $A$, a vector $b$, and a zero forcing set $Z$ of $A$ as
inputs, and updates $x_{u}$ for all $u\in V$ in each round. We set $x_{u}=0$
as initialization. Then, for each $u\in V\setminus Z$, we update $x_{u}$
according to $\left(b_{u^{\uparrow}}-\sum_{v\in
N^{+}(u^{\uparrow})\setminus\left\\{u\right\\}}A_{u^{\uparrow},v}x_{v}\right)/A_{u^{\uparrow},u}$,
iteratively in the forcing order. Note that it is equivalent to setting
$x_{u}\leftarrow\left(b_{u^{\uparrow}}-A_{u^{\uparrow},*}x\right)/A_{u^{\uparrow},u}$,
since $x_{u}$ is previously $0$.

|  Forcing$(A,Z,b)$
---
$x\leftarrow 0$
for $u\in V\setminus Z$ ordered by some forcing sequence
$x_{u}\leftarrow\left(b_{u^{\uparrow}}-A_{u^{\uparrow},*}x\right)/A_{u^{\uparrow},u}$
return $x$
Figure 2.1: The forcing operation.

The following result formalizes the relation between the solution to a given
linear system and the corresponding zero forcing set.

######  Proposition 2.2.

After the value of $x_{u}$ is updated in $\textsc{Forcing}(A,Z,b)$, we have
$A_{u^{\uparrow},*}x=b_{u^{\uparrow}}$ for all $v\leq u$. In particular, if
$Ax=b$ has a solution where $x_{Z}=0$, then $\textsc{Forcing}(A,Z,b)$ finds
such solution.

  * Proof

One can prove it by induction. The trivial case is when no $x_{u}$ has been
updated, then the conclusion is vacuously true. If some $x_{u}$ is updated by
$\textsc{Forcing}(A,Z,b)$, we let $x^{\prime}$ be the vector before the
update, and $x^{\prime\prime}$ be the one after the update. Then we get
$x^{\prime}_{v}=x^{\prime\prime}_{v}$ for all $v\neq u$, and
$x^{\prime}_{u}=0$. And
$A_{u,*}x^{\prime\prime}=A_{u,*}x^{\prime}+A_{u^{\uparrow},u}x^{\prime\prime}_{u}=A_{u,*}x^{\prime}+A_{u^{\uparrow},u}(b_{u^{\uparrow}}-A_{u^{\uparrow},*}x^{\prime})/A_{u^{\uparrow},u}=b_{u^{\uparrow}}$.
Also, for all $v<u$,
$A_{v^{\uparrow},*}x^{\prime\prime}=A_{v^{\uparrow},*}x^{\prime}=b_{v^{\uparrow}}$.

Observe that the running time of $\textsc{Forcing}(A,Z,b)$ is $O(m)$, where
$m$ is the number of non-zero elements in $A$.

######  Theorem 2.3.

Let $Z$ be a zero forcing set of $A$. The following statements are true.

  1. 1.

The columns in $V\setminus Z$ are linearly independent.

  2. 2.

If $Ax=b$, $Ay=b$, and $x_{Z}=y_{Z}$, then $x=y$.

  3. 3.

Given $x^{\prime}\in F^{V}$ such that
$\operatorname{supp}(x^{\prime})\subseteq Z$. If $Ax=b$ for some $x$ such that
$x_{Z}=x^{\prime}_{Z}$, then
$x=\textsc{Forcing}(A,Z,b-Ax^{\prime})+x^{\prime}$.

  * Proof

See [18] for proof of the first two statements. Note in [18], for the second
statement, only the $b=0$ version was proven, but the same proof still works
for general $b$. For the third statement, if $Ax=b$ where
$x_{Z}=x^{\prime}_{Z}$, then $A(x^{\prime}+(x-x^{\prime}))=b$, or in other
words, $A(x-x^{\prime})=b-Ax^{\prime}$. Hence by section 2,
$\textsc{Forcing}(A,\pi,b-Ax^{\prime})=x-x^{\prime}$, and therefore
$(x-x^{\prime})+x^{\prime}=x$.

The _LSP decomposition_ of $A\in\mdmathbb{F}^{m\times n}$ takes the form
$A=LSP$, where $L\in\mdmathbb{F}^{m\times m}$ is a lower triangular matrix
with value $1$ in the diagonals, $S\in\mdmathbb{F}^{m\times n}$ can be reduced
to a upper triangular matrix if all zero rows are deleted and elements in the
main diagonal after the zero row deletion is non-zero, and
$P\in\mdmathbb{F}^{n\times n}$ is a permutation matrix. The LSP decomposition
can be found in $O(mnr^{\omega-2})$ time [17], where $r$ is the rank of $A$.

######  Theorem 2.4.

One can build a data structure on an $n$-square matrix $A$ in $O(n)$ time,
such that $Ax=b$ can be solved in $O(n^{2})$ time for each $b$.

  * Proof

Computing the LSP decomposition of $A\in\mdmathbb{F}^{n\times n}$ takes $O(n)$
time. Given the LSP decomposition, solving $Ax=b$ can be obtained through
solving $Ly=b$ and $SPx=y$ through back substitution [15], which takes
$O(n^{2})$ time. We just have to store all the involved matrices, which takes
$O(n^{2})$ space.

##  3 Solving $Ax=b$ through zero forcing set

In this section, we first introduce the _core matrix_ that serves as the key
part of our algorithm, followed by its theoretical guarantees. Based on the
core matrix, we then present the detailed algorithm for solving $Ax=b$ with a
given zero forcing set of size $k$, as well as the corresponding computational
analysis.

To greatly simplify the exposition, we set up the instance we are working
with. Throughout this section, we fix a directed graph $G=(V,E)$, a field
$\mdmathbb{F}$, a matrix $A\in\mathcal{M}(G,\mdmathbb{F})$, a zero forcing set
$Z$ of $G$, a forcing order , and the terminals $T$ under the forcing order.
Since the operations are performed on indices, without loss of generality, we
assume that $V=[n]$ and $Z=[k]$.

Let $L(x)=\textsc{Forcing}(A,\pi,x)$, $R(x)=x-AL(x)$. As we will see, both are
linear functions. It can be also observed that
$\operatorname{supp}(R(b))\subseteq T$ and $\operatorname{supp}(L(b))\subseteq
V\setminus Z$ for all $b$. $|V|=n$, $|Z|=k$, and the number of non-zero
elements in $A$ is $m$.

###  3.1 Core matrix

As Theorem 2.3 suggests, the solution to $Ax=b$ can be obtained by knowing
only $x_{Z}$. Hence, it is natural to think of finding the correct $x_{Z}$
instead of acquiring the full information of $x$ for solving $Ax=b$.

Let $a_{v}=A_{*,v}$ be the column of $A$ indexed by $v$. Define a $k\times k$
matrix $B\in\mdmathbb{F}^{T\times Z}$ as $B=(RA)_{T,Z}$. In other words, the
$v$th column equals $R(a_{v})_{T}$. The matrix $B$ is called the _core matrix_
of $A$. We will show that for any $b$, $Ax=b$ if and only if
$Bx_{Z}=R(b)_{T}$. Hence, solving the equation $By=R(b)_{T}$ along with some
post processing is sufficient to give a solution of $Ax=b$.

First, we prove the linearity of functions $L$ and $R$.

######  Lemma 3.1.

$L$ and $R$ are linear functions.

  * Proof

The equation $Rb=(I-AL)b$ indicates that if $L$ is linear, $R$ must be linear.
Hence we only need to show the linearity of $L$.

The proof can be derived by encoding each forcing operation using matrices.
Define $D(u)\in\mdmathbb{F}^{n\times n}$ for each $u$ with
$D(u)_{u,*}=A_{u^{\uparrow},*}/A_{u^{\uparrow},u}$, $D_{v,v}=1$ if $v\neq u$,
and $0$ everywhere else. Also, define $E(u)\in\mdmathbb{F}^{n\times n}$ to be
the matrix that is $0$ everywhere except
$E(u)_{u^{\uparrow},u}=1/A_{u^{\uparrow},u}$. We can then define
$M(v)\in\mdmathbb{F}^{(2n)\times(2n)}$ for each $v\in V$, which is written in
a block form as follows

$\begin{bmatrix}-D(u)&E(u)_{u,u^{\uparrow}}\\\ 0_{n}&I_{n}\end{bmatrix},$

where $I_{n}$ and $0_{n}$ are the identity and the zero matrix of order $n$,
respectively.

Let $s=(x_{1},\ldots,x_{n},b_{1},\ldots,b_{n})^{T}$. Observe that $M(v)s$ and
$s$ differ only in one coordinate, $x_{u}$, and the value of $x_{u}$ is
$b_{u^{\uparrow}}/A_{u^{\uparrow},u}-A_{u^{\uparrow},*}/A_{u^{\uparrow},u}x$,
precisely the relation in forcing.

Now, consider
$\hat{s}=(0,\ldots,0,b_{1},\ldots,b_{n})^{T}\in\mdmathbb{F}^{2n}$, which can
be obtained by a linear transformation from $b$. Let
$\hat{s}^{\prime}=M(v_{n})\ldots M(v_{k+1})\hat{s}$, where
$v_{k+1},\ldots,v_{n}$ are the vertices in $V\setminus Z$ ordered by the
forcing order. Hence, the vector $\hat{s}^{\prime}$ can also be obtained by a
linear transformation from $b$. Assigning the values of the first $n$
coordinates of $\hat{s}^{\prime}$ to $x^{\prime}$, $x^{\prime}$ then coincides
with the value of $L(b)$. This implies that $x^{\prime}$ is obtained from $b$
through a linear transformation, and thus, $L$ is linear.

Based on the linearity above, we will abuse the notation and let $L$ and $R$
to represent the matrices for the corresponding linear transformations.

######  Lemma 3.2.

$Ra_{v}=0$ for $v\not\in Z$.

  * Proof

The forcing algorithm given in Figure 2.1 shows that if $v\neq u$, we will
then have $x_{v}=1$, and $x_{u}=0$. Note that $Ax=a_{v}$, and therefore
$Ra_{v}=a_{v}-Ax=0$.

Before delving into the key theorem, we introduce some useful notations and
definitions that will be repeatedly used in the remaining part of this
section. Let $\rank(A)=r\geq n-k$. Define
$A^{\prime}=A_{*,[k]}\in\mdmathbb{F}^{n\times k}$, and
$A^{\prime\prime}=A_{*,[n]\setminus[k]}\in\mdmathbb{F}^{n\times(n-k)}$. To
facilitate the analysis, we further decompose $A^{\prime}$ into a $1\times 2$
block matrix form $A^{\prime}=[A^{\prime}_{1}\quad A^{\prime}_{2}]$ such that
$A^{\prime}_{1}=[a_{1}\quad\ldots\quad a_{p}]$,
$A^{\prime}_{2}=[a_{p+1}\quad\ldots\quad a_{k}]$, $\rank([A^{\prime}_{1}\quad
A^{\prime\prime}])=r$, and $\rank(A^{\prime}_{1})=p=r-(n-k)$. The matrix $A$
can be then rewritten as

$A=[A^{\prime}_{1}\quad A^{\prime}_{2}\quad A^{\prime\prime}].$ (1)

The following result presents the rank-preserving property of $A^{\prime}_{1}$
under the linear mapping $R$. Based on this property, together with the
assumption of the existence of the solution to $Ax=b$, we can then arrive at
the main theoretical result of this work, which will be detailed in Theorem
3.4 and Theorem 3.5.

######  Lemma 3.3.

$\rank(RA^{\prime}_{1})=\rank(A^{\prime}_{1})$.

  * Proof

For any $s\in\mdmathbb{F}^{n}$, there exist ${}_{k+1},\cdots,{}_{n}$, such
that

$s+\sum_{i=k+1}^{n}{}_{i}a_{i}=Rs,$

which implies that $(R-I)s\in\operatorname{span}(A^{\prime\prime})$. Define
another linear mapping $Q$ as $Q=R-I$, and we have
$\operatorname{span}Q\subseteq\operatorname{span}A^{\prime\prime}$. Then, we
can write $RA_{1}^{\prime}$ as

$RA_{1}^{\prime}=(Q+I)A_{1}^{\prime}=QA_{1}^{\prime}+A_{1}^{\prime},$

where
$QA_{1}^{\prime}\in\operatorname{span}Q\subseteq\operatorname{span}A^{\prime\prime}$.
Hence, the columns of $QA_{1}^{\prime}$ and the ones of $A_{1}^{\prime}$ are
linearly independent, which immediately gives
$\rank(RA^{\prime}_{1})=\rank(A^{\prime}_{1})$.

The following theorem plays a pivotal role in the theoretical guarantees for
our algorithm.

######  Theorem 3.4.

Given $A\in\mdmathbb{F}^{n\times n}$ of the form (1). Let
$M\in\mdmathbb{F}^{k\times n}$ ($k\leq n$) such that
$A^{\prime\prime}\in\ker(M)$, and
$\rank(MA^{\prime}_{1})=\rank(A^{\prime}_{1})$. For $b\in\mdmathbb{F}^{n}$,
suppose that $Ax=b$ has a solution. If $MA^{\prime}y=Mb$ for some
$y\in\mdmathbb{F}^{k}$, then there exists $x\in\mdmathbb{F}^{n}$ such that
$Ax=b$ and $x_{Z}=y$.

  * Proof

Since $A_{2}^{\prime}\in\operatorname{span}(A_{1}^{\prime},A^{\prime\prime})$,
we can rewrite $A_{2}^{\prime}$ as
$A_{2}^{\prime}=A_{1}^{\prime}C_{1}^{\prime}+A^{\prime\prime}C^{\prime\prime}$
where $C_{1}^{\prime}\in\mdmathbb{F}^{p\times(k-p)}$ and
$C^{\prime\prime}\in\mdmathbb{F}^{(n-k)\times(k-p)}$ are coefficient matrices.
The existence of the solution to $Ax=b$ indicates that
$b\in\operatorname{span}(A_{1}^{\prime},A^{\prime\prime})$. Similarly, we can
write $b=A_{1}^{\prime}C_{b}^{\prime}+A^{\prime\prime}C_{b}^{\prime\prime}$
where $C_{b}^{\prime}\in\mdmathbb{F}^{p}$ and
$C^{\prime\prime}_{b}\in\mdmathbb{F}^{n-k}$ are coefficient vectors. Then, we
have

$\displaystyle M[A_{1}^{\prime}\quad A_{2}^{\prime}]y$
$\displaystyle=M[A_{1}^{\prime}\quad
A_{1}^{\prime}C_{1}^{\prime}+A^{\prime\prime}C^{\prime\prime}]y=M[A_{1}^{\prime}\quad
A_{1}^{\prime}C_{1}^{\prime}]y,$ $\displaystyle Mb$
$\displaystyle=M(A^{\prime}_{1}C_{b}^{\prime}+A^{\prime\prime}C_{b}^{\prime\prime})=MA_{1}^{\prime}C_{b}^{\prime},$

where the last equalities of the above two formulas follow from
$A^{\prime\prime}\in\ker(M)$. Then, the equation $MA^{\prime}y=Mb$ can be
written as $[MA_{1}^{\prime}\quad
MA_{1}^{\prime}C_{1}^{\prime}]y=MA_{1}^{\prime}C_{b}^{\prime}$. Decomposing
$y\in\mdmathbb{F}^{k}$ into $y=\begin{bmatrix}y_{1}\\\ y_{2}\end{bmatrix}$
with $y_{1}\in\mdmathbb{F}^{p},y_{2}\in\mdmathbb{F}^{k-p}$ further leads to a
new homogeneous linear system

$MA_{1}^{\prime}(y_{1}+C_{1}^{\prime}y_{2}-C_{b}^{\prime})=0.$

The condition $\rank(MA^{\prime}_{1})=\rank(A^{\prime}_{1})$ gives
$\rank(MA^{\prime}_{1})=p$, or to say, $MA^{\prime}_{1}$ is of full rank.
Thus, we have $y_{1}+C_{1}^{\prime}y_{2}-C_{b}^{\prime}=0$.

In a similar manner, we decompose $x\in\mdmathbb{F}^{n}$ into three subvectors
$x_{1}\in\mdmathbb{F}^{p}$,$x_{2}\in\mdmathbb{F}^{k-p}$, and
$x_{3}\in\mdmathbb{F}^{n-k}$. Then $Ax=b$ gives

$[A_{1}^{\prime}\quad
A_{1}^{\prime}C_{1}^{\prime}+A^{\prime\prime}C^{\prime\prime}\quad
A^{\prime\prime}]\begin{bmatrix}x_{1}\\\ x_{2}\\\
x_{3}\end{bmatrix}=A_{1}^{\prime}C_{b}^{\prime}+A^{\prime\prime}C_{b}^{\prime\prime}.$

Taking $x_{1}=y_{1}$, $x_{2}=y_{2}$ yields

$A_{1}^{\prime}y_{1}+A_{1}^{\prime}C_{1}^{\prime}y_{2}+A^{\prime\prime}C^{\prime\prime}y_{2}+A^{\prime\prime}x_{3}=A_{1}^{\prime}C_{b}^{\prime}+A^{\prime\prime}C_{b}^{\prime\prime},$

which can be rearranged as

$A_{1}^{\prime}(y_{1}+C_{1}^{\prime}y_{2}-C_{b}^{\prime})+A^{\prime\prime}(C^{\prime\prime}y_{2}+x_{3}-C_{b}^{\prime\prime})=0.$

Since $C^{\prime\prime}y_{2}+x_{3}-C_{b}^{\prime\prime}=0$, we get

$A^{\prime\prime}(C^{\prime\prime}y_{2}+x_{3}-C_{b}^{\prime\prime})=0.$

Now, we can set $x_{3}=C_{b}^{\prime\prime}-C^{\prime\prime}y_{2}$, which
immediately gives a solution of $Ax=b$ of the form

$x=\begin{bmatrix}y_{1}\\\ y_{2}\\\
C_{b}^{\prime\prime}-C^{\prime\prime}y_{2}\end{bmatrix}.$

######  Theorem 3.5.

Let $B$ be the core matrix of $A$, and $Ax=b$ has at least one solution. If
$By=R_{T,*}b$, then there exists $x$ such that $Ax=b$ and $x_{Z}=y$.

  * Proof

Since $\operatorname{supp}(Rx),\operatorname{supp}(R_{T,*}x)\in T$ for all
$x$, we have $\rank(RA^{\prime}_{1})=\rank(R_{T,*}A^{\prime}_{1})$. By Theorem
3.4 with $M=R_{T,*}$ and section 3.1, we conclude that such solution always
exists.

###  3.2 The algorithm

We first provide the algorithm FindCore, see Figure 3.1, for computing the
core matrix effectively, with the computational cost given in Theorem 3.6.

|  |  FindCore$(A,Z,b)$
---|---
$\pi,T\leftarrow$ the forcing ordering and terminal set
for $v\in Z$
$z\leftarrow$ Forcing$(A,Z,A_{*,v})$
$B_{*,v}\leftarrow(A_{*,v})_{T}-A_{T,*}z$
Compute the LSP decomposition of $B$
return $B$
Figure 3.1: Find the core matrix.

######  Theorem 3.6.

The algorithm FindCore takes $O(mk)$ time and $O(k^{2})$ space.

  * Proof

Computing $L(A_{*,v})$ for all $v\in Z$ takes $O(m(1+|Z|))=O(mk)$ time. The
computation of $B_{*,v}$ for a $v\in Z$ consists of a vector-vector difference
and a matrix-vector product, which can be implemented in $O(m)$ linear time.
Thus, computing the core matrix takes $O(mk)$ time in total.

Next, the algorithm SolveLinearSystemGivenCore$(A,Z,B,b)$ in Figure 3.2 shows
how the solution of $Ax=b$ is obtained using the computed core matrix and all
the information about the zero forcing set. The computational cost is provided
in Theorem 3.7.

|  |  SolveLinearSystemGivenCore$(A,Z,B,b)$
---|---
$\pi,T\leftarrow$ forcing sequence and terminal set $T$
$z\leftarrow\textsc{Forcing}(A,Z,b)$
$b^{\prime}\leftarrow b_{T}-A_{T,*}z$
$y\leftarrow$ solution to $By=b^{\prime}$
if $y$ does not exists
return "NO SOLUTION"
$x^{\prime}\leftarrow$ the vector where $x^{\prime}_{Z}=y$ and 0 everywhere
else
$x\leftarrow x^{\prime}+\textsc{Forcing}(A,Z,b-Ax^{\prime})$
if $Ax=b$
return $x$
else
return "NO SOLUTION"
Figure 3.2: Solve a linear system $Ax=b$ given a core matrix.

By Theorem 3.5, we know there exists a solution that matches $y$ at the zero
forcing set. By Theorem 2.3, we obtain the remaining part of the solution
through forcing.

######  Theorem 3.7.

Given a matrix $A\in\mdmathbb{F}^{n\times n}$ of $m$ non-zero elements, its
zero forcing set of size $k$, and a core matrix $B$ represented by its LSP
decomposition. The system of linear equations $Ax=b$ can be solved in
$O(k^{2}+m)$ time.

  * Proof

Following the algorithm in Figure 3.2, the computation of
$b^{\prime}=R_{T,*}b$ by forcing takes $O(m)$ time. By Theorem 2.4, if the LSP
decomposition of the matrix is given, then solving a system of linear
equations on $k$ variables and $k$ equations takes $O(k^{2})$ time. Once the
solution of the linear system for the core matrix is obtained, we can find the
solution to the original problem through forcing in $O(m)$ time. The total
running time is then $O(k^{2}+m)$.

Combining Theorem 3.6 and Theorem 3.7, we finally arrive at our main theorem.
See 1.1

##  4 Lights out game on a grid

In this section, we show that the lights out game in an $n\times n$ grid can
be solved in $O(n\log n)$ time. Consider the lights out game on an $n\times n$
grid graph. We number the vertices in position $(i,j)$ with index $(i-1)n+j$.
Hence, the vertices are in $[n^{2}]$. Let $Z=[n]$, we see that $Z$ is a zero
forcing set. If we apply Theorem 1.1 directly to the adjacency matrix of the
grid graph, then in this special case, we will obtain precisely the algorithm
of [23]. Since $m=n^{2}$ and $k=n$, the algorithm takes $O(n^{3})$ time. The
computational bottleneck is the calculation of the core matrix. Fortunately,
by exploiting the repeated structure in grids, we can significantly improve
the running time of computing the core matrix to $O(n\log n)$.

The forcing operation for the lights out game is greatly simplified, because
the field is $\mdmathbb{F}_{2}$. In this case, $x_{u}=1$ or $x_{u}=0$ can be
interpreted as pressing the button at vertex $u$ or not, respectively. The
$b_{u}$ can be understood as the state of light at vertex $u$, where the value
$1$ means on, and $0$ means off. When operating on the $u$th vertex, the
forcing operation sets $x_{u}=1$ if and only if $b^{\prime}_{u^{\uparrow}}=1$.
Here $b^{\prime}$ is the state of the board after applying all previous button
presses. In other words, the forcing operation is iteratively setting
$x_{u}=b^{\prime}_{u^{\uparrow}}$.

We want to encode the operation of forcing, where we are given the state of
first and second rows and want to compute the force operations on the second
row ensuring the states of the first row can be all $0$. It should then output
the state of the second row and the state of the third row after all the
button presses. To this end, we define such matrix to be
$N(n)\in\mdmathbb{F}_{2}^{2n\times 2n}$. The vertices of first row are indexed
from $1$ to $n$. The vertex of the second row below the vertex $i$ is $n+i$.
One can easily verify that $N$ could be written in the following block form

$N(n)=\begin{bmatrix}N^{\prime}(n)&I_{n}\\\ I_{n}&0_{n}\end{bmatrix},$

where $N^{\prime}(n)$ is the matrix satisfying that $N^{\prime}(n)_{i,j}=1$ if
and only if $|i-j|\leq 1$.

$\left[\begin{array}[]{c:c}\begin{array}[]{c c c}1&1&0\\\ 1&1&1\\\ 0&1&1\\\
\end{array}&\begin{array}[]{c c c}1&0&0\\\ 0&1&0\\\ 0&0&1\\\ \end{array}\\\
\hdashline\begin{array}[]{c c c}1&0&0\\\ 0&1&0\\\ 0&0&1\\\
\end{array}&\begin{array}[]{c c c}0&0&0\\\ 0&0&0\\\ 0&0&0\\\
\end{array}\end{array}\right]$ Figure 4.1: The matrix $N(3)$.

Now, define $M=N^{n-1}(n)$, which requires $O(n\log n)$ time to compute using
exponentiation by squaring. Let $a_{j}$ to be the $j$th column of $A$. Since
$\operatorname{supp}(a_{j})$ is a subset of the first two rows if $j\in[n]$,
we first let $a_{j}^{\prime}=(a_{j})_{[2n]}$. Next, we compute
$y=Ma_{j}^{\prime}$. Let $t$ be the last $n$ coordinates of $y$. This then
gives us the desired $R_{T,*}a_{j}$. Note that we can batch all the
multiplications together, that is, we can compute
$M[a_{j}^{\prime}|j\in[n]]=MA_{[n],[2n]}$ in one go, which takes $O(n)$
running time, and recover the desired core matrix from it. The procedure
described above is summarized in Figure 4.2.

|  |  FindGridCore$(n)$
---|---
$A\leftarrow$ adjacency matrix of an $n\times n$ grid
$M\leftarrow(N(n))^{n-1}$
return $(MA_{[n],[2n]})_{[n],[n]}$
Figure 4.2: Find core matrix for a grid graph.

At this point, we have showed a $O(n\log n)$ time algorithm to solve the
lights out problem on an $n\times n$ grid. Moreover, the algorithm can also be
applied to $n\times m$ grids with $n\leq m$, and the running time becomes
$O(n\log m+nm)$ accordingly.

## References

  * [1] Aazami, Ashkan.  Hardness results and approximation algorithms for some problems on graphs.  PhD thesis, 2008.
  * [2] AIM Minimum Rank – Special Graphs Work Group.  Zero forcing sets and the minimum rank of graphs.  Linear Algebra and its Applications, 428(7):1628–1648, April 2008\.
  * [3] Josh Alman and Virginia Vassilevska Williams.  A Refined Laser Method and Faster Matrix Multiplication, pages 522–539.
  * [4] Noga Alon and Raphael Yuster.  Matrix sparsification and nested dissection over arbitrary fields.  Journal of the ACM, 60(4):1–18, August 2013.
  * [5] Francesco Barioli, Shaun Fallat, H. Hall, Daniel Hershkowitz, Leslie Hogben, Hein Van der Holst, Bryan Shader, Bryan Shader, Bryan Shader, Bryan Shader, and Bryan Shader.  On the minimum rank of not necessarily symmetric matrices: A preliminary study.  The Electronic Journal of Linear Algebra, 18, January 2009.
  * [6] Abraham Berman, Franziska Borer, and Norbert Hungerbühler.  Lights Out on graphs.  Mathematische Semesterberichte, 68(2):237–255, October 2021.
  * [7] Boris Brimkov, Caleb C. Fast, and Illya V. Hicks.  Computational approaches for zero forcing and related problems.  European Journal of Operational Research, 273(3):889–903, 2019\.
  * [8] Daniel Burgarth and Vittorio Giovannetti.  Full control by locally induced relaxation.  Phys. Rev. Lett., 99:100501, Sep 2007.
  * [9] Nathaniel Dean, Alexandra Ilic, Ignacio Ramirez, Jian Shen, and Kevin Tian.  On the power dominating sets of hypercubes.  In 2011 14th IEEE International Conference on Computational Science and Engineering, pages 488–491, 2011.
  * [10] Rudolf Fleischer and Jiajin Yu.  A Survey of the Game “Lights Out!”, pages 176–198.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2013.
  * [11] Jim Fowler.  Scanning algorithms for _Lights Out_.  https://github.com/kisonecat/lights-out.  Accessed: 2022-07-28.
  * [12] Gene H Golub and Charles F Van Loan.  Matrix Computations.  Johns Hopkins Studies in the Mathematical Sciences. Johns Hopkins University Press, Baltimore, MD, 4 edition, February 2013.
  * [13] Leslie Hogben.  Minimum rank problems.  Linear Algebra and its Applications, 432(8):1961–1974, April 2010\.
  * [14] Markus Hunziker, António Machiavelo, and Jihun Park.  Chebyshev polynomials over finite fields and reversibility of -automata on square grids.  Theoretical Computer Science, 320(2):465–483, 2004.
  * [15] Oscar H. Ibarra, Shlomo Moran, and Roger Hui.  A generalization of the fast LUP matrix decomposition algorithm and applications.  Journal of Algorithms, 3(1):45–56, 1982.
  * [16] OEIS Foundation Inc.  Entry A117870 in the on-line encyclopedia of integer sequences.  http://oeis.org/A117870, 2022.
  * [17] Claude-Pierre Jeannerod.  LSP matrix decomposition revisited.  Technical Report 2006-28, École Normale Supérieure de Lyon, September 2006.
  * [18] Franklin H.J. Kenter and Jephian C.-H. Lin.  On the error of a priori sampling: Zero forcing sets and propagation time.  Linear Algebra and its Applications, 576:124–141, September 2019\.
  * [19] C. David Leach.  Chasing the Lights in Lights Out.  Mathematics Magazine, 90(2):126–133, April 2017.
  * [20] Richard J. Lipton, Donald J. Rose, and Robert Endre Tarjan.  Generalized Nested Dissection.  SIAM Journal on Numerical Analysis, 16(2):346–358, April 1979.
  * [21] Simone Severini.  Nondiscriminatory propagation on trees.  Journal of Physics A: Mathematical and Theoretical, 41(48):482002, oct 2008.
  * [22] Maguy Trefois and Jean-Charles Delvenne.  Zero forcing number, constrained matchings and strong structural controllability.  Linear Algebra and its Applications, 484:199–218, 2015.
  * [23] Zheng Wang (axpokl).  点灯游戏Flip game的$O(n^{3})$算法.  https://zhuanlan.zhihu.com/p/53646257, 2018.