# NoPropaganda at SemEval-2020 Task 11: A Borrowed Approach to Sequence
Tagging and Text Classification

Ilya Dimov
ABBYY ARD NLP
Lomonosov Moscow State University
<EMAIL_ADDRESS>
&Vladislav Korzun
ABBYY ARD NLP
MIPT ABBYY Lab
<EMAIL_ADDRESS>
Ivan Smurov
ABBYY ARD NLP
MIPT ABBYY Lab
<EMAIL_ADDRESS>

###### Abstract

This paper describes our contribution to SemEval-2020 Task 11: Detection Of
Propaganda Techniques In News Articles. We start with simple LSTM baselines
and move to an autoregressive transformer decoder to predict long continuous
propaganda spans for the first subtask. We also adopt an approach from
relation extraction by enveloping spans mentioned above with special tokens
for the second subtask of propaganda technique classification. Our models
report an F-score of 44.6% and a micro-averaged F-score of 58.2% for those
tasks accordingly.

##  1 Introduction

In recent years natural language processing has experienced rapid development.
In particular, application of deep learning to NLP [Collobert et al., 2011],
the introduction of pre-trained context-independent word embeddings such as
word2vec [Mikolov et al., 2013a, Mikolov et al., 2013b], GloVe [Pennington et
al., 2014] and fastText [Bojanowski et al., 2017], and utilization of RNN-
based [Hochreiter and Schmidhuber, 1997] pipelines such as CharCNN-BLSTM-CRF
[Lample et al., 2016, Ma and Hovy, 2016] allowed to improve SOTA for an
overwhelming majority of NLP tasks. In 2018 a monumental breakthrough happened
when context-dependent embeddings based on pre-trained language models
emerged, such as ELMo [Peters et al., 2018], ULMFit [Howard and Ruder, 2018],
and BERT [Devlin et al., 2018]. The resulting boost in the models’ performance
cannot be underestimated, amounting to a 30% relative error reduction rate.

Researchers compare this latest development to the introduction of ImageNet to
Computer Vision. Furthermore, as breakthroughs in CV allowed such technologies
as DeepFake to arise, similar processes are happening to NLP: the authors of
GPT2 [Radford et al., 2019] anticipated that their model could be used for
malicious purposes such as fake news generation or even impersonation. Given
that at present, the problem of media swaying public opinion to the benefit of
certain parties is more acute than ever, we must keep our news fair and
unbiased; thus, techniques for fake news and propaganda detection need to be
developed. SemEval-2020 task 11 is an essential step in that direction.

Another noteworthy consequence of contextualized embeddings introduction is
the generalization of models. While many models are developed for a particular
task, they can be effectively utilized in order to solve other tasks. Our main
contribution is the utilization of two such models for new tasks, namely
LaserTagger [Malmi et al., 2019], initially developed for summarization, which
is used for span identification, while R-BERT [Wu and He, 2019], originally
developed for relation extraction, is used for span classification.

Section 2 contains the description of our systems111The code is released at
https://github.com/hawkeoni/Semeval2020_task11., section 3 contains evaluation
results and error analysis, while part 4 contains the concluding remark.
11footnotetext:  This work is licensed under a Creative Commons Attribution
4.0 International Licence. Licence details:
http://creativecommons.org/licenses/by/4.0/.

##  2 System description

All of the following models based on BERT [Devlin et al., 2018] use BERTBASE
pretrained model weights. Hyperparameter tuning and cross-validation was
performed on 5 folds split by articles. The data used for training and
evaluation consists of various news articles and is described in detail in [Da
San Martino et al., 2019]. The tasks will be briefly explained in the
following sections; their full description can be found in the [Da San Martino
et al., 2020].

###  2.1 Span Identification

The goal of this task is to find continuous spans of text which contain any
propaganda technique. Consider the following example, which has loaded
language marked up in bold and double squared brackets:

> It’s essentially an admission of guilt, given that it is [[absolutely
> ludicrous]] to think that “national security” would be threatened by the
> release of the CIA’s long-secret JFK-assassination-related records.

We treat this task as a binary sequence tagging where each token is assigned a
label: 0 - for normal text, 1 - for propaganda text. Although there are spans
which overlap multiple sentences, we decided to split such spans and solve the
problem on a sentence level. This approach has a distinct disadvantage, as the
connection between sentences is lost.

At first we started with several baseline bidirectional LSTM models with
linear layer for token class prediction:

  * •

biLSTMGLOVE + charLSTM \- a bidirectional LSTM over GloVe embeddings and a
character-level LSTM. This is a simple sequence tagging model proposed
originally in [Lample et al., 2016].

  * •

biLSTMELMO \- a bidirectional LSTM over ELMo embeddings from [Peters et al.,
2018].

Preliminary results of those models achieved an F-score of 24-26% on the
development set. Postprocessing model predictions by labeling all tokens
between the first sentence-wise positive label and the last sentence-wise
positive label as propaganda pushed F-score to 34.5%. A CRF instead of a
linear layer provided no significant benefit.

Our next model BERTLINEAR, which consists of a linear layer over BERTBASE,
achieved a stronger baseline with F-score of 40.8% without any postprocessing.
However, unlike in previous models postprocessing only made things worse,
degrading prediction quality.

Generating token tags with a linear layer has one big flaw: predictions for
each token are made independently. CRF mitigates this problem and improves the
metrics, but we decided to employ LaserTagger \- an autoregressive transformer
decoder from [Malmi et al., 2019]. One of the critical ideas of LaserTagger is
directly consuming corresponding encoder activation without learning encoder-
decoder attention weights. LaserTagger achieved an even stronger baseline with
an F-score of 42% on development set, which at the time was only several
points away from top-performing teams. It was decided to keep this model and
extensively tune it during further experiments.

Figure 1: LaserTagger model architecture. Ti are token embeddings for BERT, Ei
are BERT outputs, Hi represent LaserTagger hidden states and Li are predicted
labels.

We have employed several techniques to boost model performance:

  * •

We fine-tune BERT as a masked language model on task corpus for 3 epochs
without next sentence prediction loss. This provides no noticeable increase in
F-score.

  * •

Teacher forcing (scheduled sampling) [Bengio et al., 2015] is a method of
improving performance for autoregressive models. As LaserTagger predictions
explicitly depend on predictions from the previous step, errors accumulated
earlier may lead to a completely wrong result. Teacher forcing aims to close
the gap between training and inference by supplying a model with its own
(possibly wrong) predictions at a specific rate during training. This allows
the model to explore more of its state space during training, which in turn
increases its robustness during inference. We linearly decay teacher forcing
rate from 1 (only correct labels) to 0 (only model predictions) during
training. The downside of this approach is a significant slowdown in training
as we do not utilize the transformer’s capability to process the whole
sequence at once.

  * •

Label smoothing was used to lower the number of false-positive predictions and
mitigating minor data discrepancies like inconsistent inclusion of dots and
quotations in propaganda spans, which were caused by limiting spans to
sentences and other preprocessing imperfections.

Our final submission was a single LaserTagger model with BERTBASE encoder,
single-layered decoder with hidden dimension of 128 and 4 attention heads. The
model was trained on an effective batch size of 32 achieved by gradient
accumulation over every 2 steps; we used Adam optimizer with learning rate of
2e-5 and learning rate warmup for 10% of the steps. Teacher forcing was set to
linearly decay from 1 to 0 for the duration of the whole training. This model
achieved F-score of 46.1% on development set and 44.6% on test set.

###  2.2 Technique Classification

The second subtask was to classify spans from previous subtasks into 14
classes of propaganda techniques. As previously mentioned, spans can overlap
several sentences, so for this subtask, we had to use all overlapping
sentences for correct classification. As spans with different techniques may
overlap, or a single span may be associated with several techniques, we treat
this subtask as a multilabel classification problem.

Our main idea was taken from R-BERT [Wu and He, 2019], where special tokens
were inserted around entities to highlight them in the task of relation
extraction. As in R-BERT, we surround propaganda spans with special tokens,
namely the ^ token. BERT representations of propaganda span tokens are then
averaged or taken as a weighted sum based on results of the linear layer
applied to BERT token activations in the following way: let
$e_{i}...e_{i+k}\in R^{n}$ be BERT representations of propaganda tokens in the
current span and $w\in R^{n},b\in R$ \- learnable parameters. Then the
weighted sum is calculated as $\sum_{j=i}^{i+k}\alpha_{j}e_{j}$, where
$\alpha_{j}$ is the result of softmax function over scalar product of $e_{j}$
and $w$:

$\alpha_{j}=\frac{\exp{((w,e_{j})+b)}}{\sum_{j=i}^{i+k}{\exp{((w,e_{j})+b)}}}$
(1)

The resulting vector is concatenated with BERT output of the [CLS] token for
final multilabel prediction. We also combine these approaches with previously
mentioned BERT intask finetuning [Sun et al., 2019], which proves to be
fruitful for this task as it is not too dissimilar to text classification. We
additionally train a simple BERT model, which predicts technique class using
only [CLS] token representation. All models were trained with linear warmup
for 10% of the training steps using the Adam optimizer with learning rate of
2e-5.

Our final submission for task TC consists of predictions averaged from several
models, which are described in Section 3.2.

Figure 2: Task TC model architecture based on R-BERT. Ti are token embeddings
for BERT; ST0 and ST1 represent special tokens used to highlight the
propaganda span; Ei and ESTi are BERT outputs

.

##  3 Experimental Results

###  3.1 Span Identification Results

Model | F1 | Precision | Recall
---|---|---|---
biLSTMGLOVE + charLSTM | 0.249 | 0.302 | 0.212
biLSTMGLOVE + charLSTM \+ postprocessing | 0.349 | 0.304 | 0.411
biLSTMELMO \+ postprocessing | 0.345 | 0.327 | 0.367
BERTLINEAR | 0.408 | 0.357 | 0.477
BERT + CRF | 0.362 | 0.301 | 0.454
BERT + LSTM + CRF | 0.414 | 0.363 | 0.482
LaserTagger | 0.420 | 0.384 | 0.464
LaserTagger + TF | 0.451 | 0.423 | 0.482
LaserTagger + TF + LS | 0.461 | 0.406 | 0.533
LaserTagger + TF + LS (test) | 0.446 | 0.556 | 0.373
Table 1:  Task SI model results on development and test sets. TF and LS stand
for teacher forcing and label smoothing accordingly.

As Table 1 shows BERT-based approaches significantly outperform suggested LSTM
baselines. Main improvements over simple BERTLINEAR model is LaserTagger
decoder, which eliminates the flaw of independent tag prediction. It was
further improved by linearly decreasing teacher forcing rate over the training
and label smoothing. Our model scored 7th place out of 35 in the final table
of the competition.

Boundary errors
---
Gold: This is a [[smear]], plain and simple.
Prediction: This is [[a smear]], plain and simple.
Gold: Luke Harding and the Guardian Publish Still More [[Blatant MI6 Lies]]
Prediction: Luke Harding and the Guardian Publish Still More [[Blatant]] MI6
Lies
Quotation and punctuation errors
Gold: While trying to appear as a moderate, ”Gillum is a Progressive” and ”He
is a part of the [[crazy, crazy, crazies]].”
Prediction: While trying to appear as a moderate, ”Gillum is a Progressive”
and ”[[He is a part of the crazy, crazy, crazies]].”
Gold: ”[[I don’t care who you are]], determine the scope and the parameters of
an investigation.
Prediction: ”[[I don’t care who you are, determine the scope and the
parameters of an investigation]].
Gold: “[[This is devastating]] for our Democratic Latino community.
Prediction: “[[This is devastating for our Democratic Latino community.]]
Table 2:  LaserTagger + TF + LS errors.

Table 2 illustrates the errors of our final model. Besides sometimes
incorrectly adding or losing boundary propaganda span tokens, LaserTagger + TF
+ LS seems to have a problem with quotation and punctuation. When it
encounters quoted text or a part of a sentence separated by commas, it either
labels the whole highlighted area as propaganda or misses it completely. This
may be partly attributed to the nature of quotes - speech often conveys more
emotional words than written text. Another problem is that when the model
finds comma, it often closes the propaganda span and does not reopen it.

###  3.2 Technique Classification Results

Model | micro-averaged F-measure
---|---
BERTCLS | 0.573
R-BERT | 0.592
R-BERTFT | 0.584
R-BERTW | 0.590
R-BERTW + FT | 0.590
All models | 0.606
All models (test) | 0.582
Table 3:  Task TI model results on development and test sets. W and FT stand
for weighted sum instead of averaging BERT token representation and intask
finetuning accordingly.

All models from the table use the ^ special symbol to highlight propaganda
spans because more than 20% of the sentences contain multiple propaganda
spans. For example, if the model has 5 tokens $e_{1}...e_{5}$ and two
propaganda spans $e_{1},e_{2}$ and $e_{4},e_{5}$ we would create two training
samples with the use token ^, which we will write for readability as st:
$st,e_{1},e_{2},st,e_{3},e_{4},e_{5}$ and
$e_{1},e_{2},e_{3},st,e_{4},e_{5},st.$ This allows us to point our model to a
specific propaganda span inside a sentence.

Models from Table 3 are the following:

  * •

BERTCLS \- a model that predicts technique solely from [CLS] token embedding.

  * •

R-BERT \- a model that predicts technique from [CLS] token embedding and
averaged propaganda span representation;

  * •

R-BERTW \- a model that predicts technique from [CLS] token embedding and
weighted propaganda span representation;

It is important to understand that due to class imbalance models with the same
micro-averaged F-score correspond to different F-scores over specific
techniques. While ensembling the models increases overall quality, two
classes, namely ”Bandwagon, Reduction Ad Hitlerum” and ”Whataboutism, Straw
Men, Red Herring”, are entirely missed by all our models. This is not
surprising: those techniques have a very limited amount of positive samples,
which is the reason they were merged into classes. Class wise F-score seems to
respond to the quantity of corresponding class in training dataset: our model
performs best on Loaded Language and Name Calling, Labeling followed by Flag-
Waving and Doubt.

We did not implement anything to boost performance on underrepresented
classes, while some of the other teams found a way to get an F-score of 20%
and higher. Our model scored 6th place out of 31 in the final table of the
competition.

###  3.3 Dataset Exploration

During the competition, most teams had a substantial shift between recall and
precision: precision was 10-20% lower than recall. On the contrary, test
dataset evaluation resulted in a different picture: top-performing teams’
precision is 5-15% higher than recall. Additionally, on task TC all teams
suffered a downgrade in micro-averaged F-score of approximately 5%.

This can be attributed to the fact that the datasets have imbalanced topic
distributions. While the training dataset is not very large, it is well-
rounded and covers a wide range of topics. The development set is focused on
the US news, while the test set seems to have almost an even split between the
UK and the US news.

##  4 Conclusion

In this work, we propose to use two models: autoregressive transformer decoder
for sequence tagging and R-BERT-like approach for propaganda spans technique
classification. The former is taken from such tasks as sentence fusion,
sentence splitting, and abstractive summarization, while the latter was
originally used for relation extraction.

Instead of exploring various pretrained transformers like BERTLARGE, RoBERTa,
ALBERT and XLNet or engineering domain-specific features and approaches, we
focused on improving architectures and trying to push them as far as possible.

In the future the results of this task may be used for propaganda detection in
various media as well as unbiasing texts and making them fair.

## References

  * [Bengio et al., 2015]  Samy Bengio, Oriol Vinyals, Navdeep Jaitly, and Noam Shazeer.  2015\.  Scheduled sampling for sequence prediction with recurrent neural networks.  In Advances in Neural Information Processing Systems, pages 1171–1179.
  * [Bojanowski et al., 2017]  Piotr Bojanowski, Edouard Grave, Armand Joulin, and Tomas Mikolov.  2017\.  Enriching word vectors with subword information.  Transactions of the Association for Computational Linguistics, 5:135–146.
  * [Collobert et al., 2011]  Ronan Collobert, Jason Weston, Léon Bottou, Michael Karlen, Koray Kavukcuoglu, and Pavel Kuksa.  2011\.  Natural language processing (almost) from scratch.  Journal of machine learning research, 12(Aug):2493–2537.
  * [Da San Martino et al., 2019]  Giovanni Da San Martino, Seunghak Yu, Alberto Barrón-Cedeño, Rostislav Petrov, and Preslav Nakov.  2019\.  Fine-grained analysis of propaganda in news articles.  In Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing, EMNLP-IJCNLP 2019, EMNLP-IJCNLP 2019, Hong Kong, China, November.
  * [Da San Martino et al., 2020]  Giovanni Da San Martino, Alberto Barrón-Cedeño, Henning Wachsmuth, Rostislav Petrov, and Preslav Nakov.  2020\.  SemEval-2020 task 11: Detection of propaganda techniques in news articles.  In Proceedings of the 14th International Workshop on Semantic Evaluation, SemEval 2020, Barcelona, Spain, September.
  * [Devlin et al., 2018]  Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  2018\.  BERT: pre-training of deep bidirectional transformers for language understanding.  CoRR, abs/1810.04805.
  * [Hochreiter and Schmidhuber, 1997]  Sepp Hochreiter and Jürgen Schmidhuber.  1997\.  Long short-term memory.  Neural computation, 9(8):1735–1780.
  * [Howard and Ruder, 2018]  Jeremy Howard and Sebastian Ruder.  2018\.  Universal language model fine-tuning for text classification.  arXiv preprint arXiv:1801.06146.
  * [Lample et al., 2016]  Guillaume Lample, Miguel Ballesteros, Sandeep Subramanian, Kazuya Kawakami, and Chris Dyer.  2016\.  Neural architectures for named entity recognition.  arXiv preprint arXiv:1603.01360.
  * [Ma and Hovy, 2016]  Xuezhe Ma and Eduard Hovy.  2016\.  End-to-end sequence labeling via bi-directional lstm-cnns-crf.
  * [Malmi et al., 2019]  Eric Malmi, Sebastian Krause, Sascha Rothe, Daniil Mirylenka, and Aliaksei Severyn.  2019\.  Encode, tag, realize: High-precision text editing.  arXiv preprint arXiv:1909.01187.
  * [Mikolov et al., 2013a]  Tomas Mikolov, Kai Chen, Greg Corrado, and Jeffrey Dean.  2013a.  Efficient estimation of word representations in vector space.  CoRR, abs/1301.3781.
  * [Mikolov et al., 2013b]  Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg S Corrado, and Jeff Dean.  2013b.  Distributed representations of words and phrases and their compositionality.  In C. J. C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Q. Weinberger, editors, Advances in Neural Information Processing Systems 26, pages 3111–3119. Curran Associates, Inc.
  * [Pennington et al., 2014]  Jeffrey Pennington, Richard Socher, and Christopher D. Manning.  2014\.  Glove: Global vectors for word representation.  In Empirical Methods in Natural Language Processing (EMNLP), pages 1532–1543.
  * [Peters et al., 2018]  Matthew E Peters, Mark Neumann, Mohit Iyyer, Matt Gardner, Christopher Clark, Kenton Lee, and Luke Zettlemoyer.  2018\.  Deep contextualized word representations.  arXiv preprint arXiv:1802.05365.
  * [Radford et al., 2019]  Alec Radford, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei, and Ilya Sutskever.  2019\.  Language models are unsupervised multitask learners.  OpenAI Blog, 1(8):9.
  * [Sun et al., 2019]  Chi Sun, Xipeng Qiu, Yige Xu, and Xuanjing Huang.  2019\.  How to fine-tune bert for text classification?  In China National Conference on Chinese Computational Linguistics, pages 194–206. Springer.
  * [Wu and He, 2019]  Shanchan Wu and Yifan He.  2019\.  Enriching pre-trained language model with entity information for relation classification.  In Proceedings of the 28th ACM International Conference on Information and Knowledge Management, pages 2361–2364.


# Tractable Identification of Electric Distribution Networks

Ognjen Stanojev, Lucien Werner, Steven Low and Gabriela Hug This research was
supported by the Swiss National Science Foundation under NCCR Automation,
grant agreement 51NF40_180545.Ognjen Stanojev and Gabriela Hug are with the
Power Systems Laboratory, ETH Zürich, Zürich, Switzerland, emails: {ognjens,
<EMAIL_ADDRESS>Werner and Steven Low are with the Department of
Computing and Mathematical Sciences at California Institute of Techonology,
Pasadena, California, USA, emails: {lwerner<EMAIL_ADDRESS>
###### Abstract

The identification of distribution network topology and parameters is a
critical problem that lays the foundation for improving network efficiency,
enhancing reliability, and increasing its capacity to host distributed energy
resources. Network identification problems often involve estimating a large
number of parameters based on highly correlated measurements, resulting in an
ill-conditioned and computationally demanding estimation process. We address
these challenges by proposing two admittance matrix estimation methods. In the
first method, we use the eigendecomposition of the admittance matrix to
generalize the notion of stationarity to electrical signals and demonstrate
how the stationarity property can be used to facilitate a maximum a posteriori
estimation procedure. We relax the stationarity assumption in the second
proposed method by employing Linear Minimum Mean Square Error (LMMSE)
estimation. Since LMMSE estimation is often ill-conditioned, we introduce an
approximate well-conditioned solution. Our quantitative results demonstrate
the improvement in computational efficiency compared to the state-of-the-art
methods while preserving the estimation accuracy.

##  I Introduction

Electric distribution networks are a vital component of the energy
infrastructure, serving as the final layer in power delivery to residential
and commercial users. The increasing integration of renewable energy sources
and the implementation of decarbonization policies require modernization of
the present control and monitoring practices in power distribution systems.
The admittance matrix is at the heart of numerous power system analysis
techniques, including (optimal) power flow, state estimation, and short
circuit analysis [1]. It bears the structure of graph Laplacian matrices [2],
thus unambiguously explaining the network topology and related line
parameters. However, distribution utilities often lack accurate topology and
parameter information, hindering the construction of the admittance matrix and
the use of the available analysis tools [3].

The recent installation of a significant number of micro Phasor Measurement
Units ($\mu$PMUs) [4] and smart meters [5] in distribution grids provides
network operators with high-precision and high-sampling-rate measurements.
These data streams enhance the observability of distribution grids and enable
network identification. In general, network identification problems involve
determining network connectivity (i.e., topology) [3, 6, 7, 8], line
parameters [9], or both [10, 11, 12, 13, 14, 15, 16, 17], using bus voltage,
current injection, or branch flow measurements. In this work, we address the
problem of estimating the admittance matrix using $\mu$PMU measurements of bus
voltage and current injection phasors. Although prior research has examined
various network identification problems, we only provide a focused overview of
the literature on estimating the admittance matrix in the following paragraph.
Please refer to [3] for a more detailed review.

In [10], matrix least squares estimation is applied on phasor measurements to
approximate the admittance matrix. A constrained least squares approach is
developed in [11] to enforce the Laplacian matrix structure in the least
squares estimate. Instead of batch processing, a recursive least squares
method is developed in [12] to enable frequent online updates. In [13], a
sparsity promoting $\ell_{1}$-norm regularizer is introduced to enhance the
least squares estimation when the admittance matrix is known to be sparse. An
alternative approach to promote sparsity is proposed in [14], where lines with
small conductance values are progressively removed after performing the least
squares estimation. The works above assume noise-free measurements of the
independent variables, which leads to biased estimates when using realistic
data with errors in all measurements (variables). This limitation of the least
squares approaches can be successfully tackled by error-in-variables methods,
such as total least squares [15]. A weighted total least squares method is
introduced in [16] and then extended in [17] to a Bayesian framework that
allows exploiting different forms of prior knowledge of the admittance matrix,
thus creating a flexible framework that can achieve high estimation accuracy.

Despite previous methods laying a solid foundation for admittance matrix
estimation, challenges involving ill-conditioning and high computational and
memory requirements in the estimation process remain unaddressed. Poor
conditioning is common to least squares approaches [10, 11, 13] and arises
even in the estimation of small-size networks due to high correlations in
voltage or current measurements. Significant computational burden and memory
requirements arise when solving the weighted total least squares [16, 17]
since a substantial number of measurements are required for accuracy, and a
large number of parameters contained within the admittance matrix need to be
estimated.

In this paper, we address the aforementioned challenges by proposing two
admittance matrix estimation methods that are computationally efficient and
numerically stable. The first method is motivated by the recent developments
in the graph signal processing community [18] on the identification of graph
filters [19, 20]. Expanding on these works, we use the eigendecomposition of
the admittance matrix to generalize the notion of stationarity to electrical
signals in power networks with a constant reactance-resistance ratio.
Subsequently, we demonstrate how a Maximum a Posteriori (MAP) estimation
method (resembling [16, 17]) can be simplified when the current injections are
stationary.

In the second proposed method, we relax the adopted assumptions and consider a
Linear Minimum Mean Square Error (LMMSE) estimation method which is applicable
to general power networks and generic current statistics. The solution to
LMMSE is known as the Wiener filter, which may suffer from poor conditioning
and might not respect the Laplacian structure of the admittance matrix. To
address the ill-conditioning issue, we introduce an approximate solution based
on eigenvalue truncation. Furthermore, we demonstrate that the Laplacian
structure can be enforced via a postfiltering procedure without significant
additional computational effort. Previous works which examined the voltage and
current injection statistics in a similar way focused on topology
identification rather than admittance matrix estimation [6, 7].

Notation. We denote the sets of real and complex numbers by $\mathbb{R}$ and
$\mathbb{C}$. Given a matrix $A$, $A^{\top}$ denotes its transpose, $A^{*}$
denotes the entrywise conjugate, and $A^{\mathsf{H}}=(A^{*})^{\top}$ denotes
the conjugate transpose. For column vectors $x\in\mathbb{C}^{n}$ and
$y\in\mathbb{C}^{m}$, we use
$(x,y)\coloneqq[x^{\top},y^{\top}]^{\top}\in\mathbb{C}^{n+m}$ to denote a
stacked vector. For a random vector $X$, we use $\mathbb{E}[X]$ to denote its
mean and $\Sigma_{X}$ to denote its covariance matrix. Finally,
$\mathcal{I}_{n}$ denotes the $n\times n$ identity matrix, $\mathbbl{1}_{n}$
and $\mathbbl{0}_{n}$ are $n$-dimensional vectors of all ones and zeros,
respectively.

##  II Admittance Matrix Model of Power Grids

Consider a static (steady-state), single-phase equivalent distribution network
composed of $n$ nodes $\mathcal{N}=\\{1,\dots,n\\}$ and $m$ undirected
branches $\mathcal{E}\subseteq\mathcal{N}\times\mathcal{N}$. The network is
modeled as a connected and undirected graph
$\mathcal{G}\coloneqq(\mathcal{N},\mathcal{E},\mathcal{W})$, with complex-
valued edge weights
$\mathcal{W}\coloneqq\\{y_{ij}\in\mathbb{C}:y_{ij}=g_{ij}+\mathrm{j}b_{ij},g_{ij}>0,b_{ij}\leq
0,\forall\\{i,j\\}\in\mathcal{E}\\}$ representing the series admittances in
the standard lumped $\pi$-model of a transmission line. For the purpose of
defining the incidence matrix of $\mathcal{G}$, let us assign to each edge a
unique identifier $e=\\{1,\dots,m\\}$ and an arbitrary orientation.
Admittances connected to the ground are defined by $y_{i0}\coloneqq
g_{i0}+\mathrm{j}b_{i0}$, with $g_{i0}\geq 0,b_{i0}\geq 0,\forall
i\in\mathcal{N}$, and referred to as shunt admittances. Each node
$i\in\mathcal{N}$ in the network is associated with a nodal current injection
$\bar{I}_{i}\in\mathbb{C}$ and a nodal voltage $\bar{V}_{i}\in\mathbb{C}$.
Kirchhoff’s and Ohm’s laws lead to the following model of the considered
electric network [2]:

$\bar{I}={B}\operatorname{diag}(\\{y_{e}\\}_{e=1}^{m}){B}^{\top}\bar{V}+\operatorname{diag}(\\{y_{i0}\\}_{i=1}^{n})\bar{V}={Y}\bar{V},$
(1)

where $\bar{V}\coloneqq(\bar{V}_{1},\dots,\bar{V}_{n})$ and
$\bar{I}\coloneqq(\bar{I}_{1},\dots,\bar{I}_{n})$ collect the bus voltages and
the current injections, respectively, and $B\in\\{-1,0,1\\}^{n\times m}$ is
the node-edge incidence matrix of $\mathcal{G}$.

######  Definition 1.

The admittance matrix $Y$ is a complex symmetric matrix, with diagonal
elements given by $Y_{ii}=\sum_{j=1,j\neq i}^{n}y_{ij}+y_{i0},\forall
i\in\mathcal{N}$, and off-diagonal elements defined by
$Y_{ij}=-y_{ij},\forall\\{i,j\\}\in\mathcal{E}$ and $Y_{ij}=0$ otherwise.

Under the adopted assumptions on the real and imaginary parts of the series
and shunt admittances, the necessary and sufficient condition for the
invertibility of $Y$ is the existence of at least one shunt admittance. The
assumptions made are reasonable for distribution networks, and we refer the
reader to [21, 22] for a broader discussion on the invertibility of $Y$.

Therefore, singular admittance matrices under the adopted assumptions have
zero row sums, i.e., $Y\mathbbl{1}_{n}=\mathbbl{0}_{n}$ iff $Y$ is singular.
The linear map defined by the singular admittance matrix
$Y:\mathbb{C}^{n}\rightarrow\mathbb{C}^{n}$ has the nullspace of dimension one
consisting of vectors in
$\mathrm{span}(\mathbbl{1}_{n})=\\{\alpha\mathbbl{1}_{n},\alpha\in\mathbb{C}\\}$.
Hence, the Moore-Penrose pseudoinverse of $Y$, denoted by $Y^{\dagger}$, can
be used to form a subspace of solutions to (1), given by

$\bar{V}=Y^{\dagger}\bar{I}+\alpha\mathbbl{1}_{n},$ (2)

where $\alpha\in\mathbb{C}$. The preceding relationship holds if and only if
the current injections are balanced, that is,
$\mathbbl{1}_{n}^{\top}\bar{I}=0$. The same relationship can be used for
invertible $Y$ in which case $\bar{I}\in\mathbb{C}^{n}$ is unrestricted and
$\alpha=0$. Equation (2) represents the so-called impedance matrix model of
the network.

It is evident from Definition 1 that $Y$ is not necessarily a normal matrix
since it is non-Hermitian complex symmetric. Hence, it is not always unitarily
diagonalizable. Nevertheless, if the conductance-susceptance ratio111The
conductance-susceptance ratio is more commonly referred to as the reactance-
resistance ratio or the “$x/r$ ratio” in the power systems literature. These
terms are used interchangeably in this work. is identical for all lines across
the network, then $Y$ is normal [21]. This assumption is typically valid for
lines at the same voltage level.

######  Theorem 1 (Pseudoinverse of $Y$).

Suppose the conductance-susceptance ratio $\frac{g_{ij}}{b_{ij}}$ is the same
for all $(i,j)\in\mathcal{E}$. Then

  1. 1.

$Y$ is a normal matrix and has a spectral decomposition
$Y=\mathbbl{W}\mathbbl{\Lambda}\mathbbl{W}^{\mathsf{H}}$, where
$\mathbbl{\Lambda}$ is the diagonal matrix with the eigenvalues of $Y$ on its
diagonal, and the columns of $\mathbbl{W}$ are the corresponding eigenvectors.

  2. 2.

The Moore-Penrose pseudoinverse of $Y$ is
$Y^{\dagger}=\mathbbl{W}\mathbbl{\Lambda}^{\dagger}\mathbbl{W}^{\mathsf{H}}$,
where $\mathbbl{\Lambda}^{\dagger}$ is the diagonal matrix obtained from
$\mathbbl{\Lambda}$ by replacing nonzero eigenvalues of $Y$ by their
reciprocals.

##  III Network Identification Problem

###  III-A Problem Formulation

We take a statistical perspective and consider zero-mean random vectors
$I\coloneqq\bar{I}-\mathbb{E}[\bar{I}]$ and
$V\coloneqq\bar{V}-\mathbb{E}[\bar{V}]$, defined such that their means satisfy
$\mathbb{E}[I]=Y\mathbb{E}[V]$. Let us further define the voltage covariance
matrix $\Sigma_{V}=\mathbb{E}[{V}{V}^{\mathsf{H}}]$, the current injection
covariance matrix $\Sigma_{I}=\mathbb{E}[{I}{I}^{\mathsf{H}}]$, and their
cross-covariance $\Sigma_{{I}{V}}=\mathbb{E}[{I}{V}^{\mathsf{H}}]$. The
considered network identification problem is formally stated in the following.

######  Problem 1.

Given a set
$\mathcal{S}\coloneqq\\{(\tilde{{V}}^{1},\tilde{{I}}^{1}),\dots,(\tilde{{V}}^{N},\tilde{{I}}^{N})\\}$
of $N$ pairs of noisy, zero-centered bus voltage and current injection
measurements pertaining to different steady-state operating points, the
objective is to infer the underlying distribution network – its edges and the
associated admittance values – or equivalently, the admittance matrix $Y$
(Def. 1).

The $N$ measurements of voltage and current phasors are further collected in
matrices
$\tilde{\mathbbl{V}}=\begin{bmatrix}\tilde{V}^{1}&\tilde{V}^{2}&\dots&\tilde{V}^{N}\end{bmatrix}$
and
$\tilde{\mathbbl{I}}=\begin{bmatrix}\tilde{{I}}^{1}&\tilde{{I}}^{2}&\dots&\tilde{{I}}^{N}\end{bmatrix}$
that will be used to characterize solutions to Problem 1. We are putting forth
the following two assumptions to limit the scope of our analysis.

######  Assumption 1.

The measurement matrices are assumed to be full rank, i.e.,
$\operatorname{rank}(\tilde{\mathbbl{V}})=n$ and
$\operatorname{rank}(\tilde{\mathbbl{I}})=n$.

######  Assumption 2.

The nodes $j$ where the current is neither injected nor extracted ($I_{j}=0$)
have been removed a priori by applying the Kron [23] or subKron [10]
reduction.

The first assumption implies that $N\geq n$. Additionally, experiment design
[24] might be required to guarantee that the measurement matrices are full
rank. The second assumption implies that all the nodes $i$ at which $I_{i}\neq
0$ are collected in $\mathcal{N}$ and are assumed to be observed. Note that
Kron reduction of the network may lead to a non-sparse admittance matrix, and
the Kron-reduced graph $\mathcal{G}$ may not be a tree graph, which are common
assumptions in distribution network studies [8].

###  III-B Measurement Model

The available $\mu$PMU or smart meter measurements are corrupted by
measurement noise. We adopt a generic linear statistical model to represent
the individual bus voltage $\tilde{V}\in\mathbb{C}^{n}$ and current injection
$\tilde{I}\in\mathbb{C}^{n}$ observations as follows:

$\tilde{V}=V+\varepsilon_{v},\qquad\tilde{I}=I+\varepsilon_{i},$ (3)

where $\varepsilon_{v}$ and $\varepsilon_{i}$ are complex random vectors
describing the measurement noise. We assume that the noise vectors follow
uncorrelated complex multivariate Gaussian distributions with zero mean:
$\varepsilon_{v}\sim\mathcal{NC}(\mathbbl{0}_{n},\sigma_{v}^{2}\mathcal{I}_{n}),\varepsilon_{i}\sim\mathcal{NC}(\mathbbl{0}_{n},\sigma_{i}^{2}\mathcal{I}_{n})$.
The covariances $\sigma_{v}$ and $\sigma_{i}$ might be time-varying, but there
is no temporal or spatial correlation in the measurement noise. Note that
other noise models may be applicable [16, 25].

##  IV Network Identification under Stationary Current Injections

This section extends the standard notions of wide-sense stationarity in
discrete time to define stationarity with respect to $Y$ for networks that
satisfy the following assumption.

######  Assumption 3.

The conductance-susceptance ratio $\frac{g_{ij}}{b_{ij}}$ is the same for all
lines $(i,j)\in\cal E$ across the network.

It is first demonstrated that the eigenvectors of $Y$ can be identified from
the voltage covariance matrix when the current injections satisfy the
stationarity property. Subsequently, we employ a MAP procedure to estimate the
corresponding eigenvalues, hence identifying $Y$ according to Theorem 1.

###  IV-A Recovering the Eigenvectors of $Y$

######  Definition 2.

Given an admittance matrix $Y$ with a spectral decomposition
$Y=\mathbbl{W}\mathbbl{\Lambda}\mathbbl{W}^{\mathsf{H}}$, a zero-mean random
variable $X$ is said to be Wide-Sense Stationary (WSS) with respect to $Y$ if
its covariance matrix $\Sigma_{X}$ also has a spectral decomposition with the
unitary $\mathbbl{W}$, i.e.,
$\Sigma_{X}=\mathbbl{W}\mathbbl{\Lambda}_{X}\mathbbl{W}^{\mathsf{H}}$ where
$\mathbbl{\Lambda}_{X}$ is a diagonal matrix of eigenvalues of $\Sigma_{X}$.

A practically relevant example of a variable that is WSS with respect to $Y$
is white noise $W$, characterized by $\mathbb{E}[W]=0$ and
$\mathbb{E}[WW^{\mathsf{H}}]=\sigma^{2}\mathcal{I}_{n}$. Current injections in
distribution grids are predominantly determined by loads that can reliably be
modeled as white noise over short time intervals (on the order of seconds).
Hence, the assumption of white noise current injection statistics has been
common in the distribution network identification literature [26]. The
following proposition establishes a connection between the voltage statistics
and the admittance matrix when the network is subjected to stationary current
injections.

######  Proposition 1.

Let $I$ be WSS with respect to $Y$ and $I=YV$. Then $V$ is also WSS with
respect to $Y$.

###### Proof.

Since $V$ is already defined to be zero-mean, we only need to show that the
covariance matrix $\Sigma_{V}$ of $V$ is unitarily diagonalized by
$\mathbbl{W}$, which is derived by considering that

$\displaystyle\Sigma_{V}$
$\displaystyle=\mathbb{E}[VV^{\mathsf{H}}]=\mathbb{E}[Y^{\dagger}I(Y^{\dagger}I)^{\mathsf{H}}]=Y^{\dagger}\Sigma_{I}(Y^{\dagger})^{\mathsf{H}}$
$\displaystyle=\mathbbl{W}\mathbbl{\Lambda}^{\dagger}\mathbbl{W}^{\mathsf{H}}\mathbbl{W}\mathbbl{\Lambda}_{I}\mathbbl{W}^{\mathsf{H}}(\mathbbl{W}\mathbbl{\Lambda}^{\dagger}\mathbbl{W}^{\mathsf{H}})^{\mathsf{H}}$
$\displaystyle=\mathbbl{W}(|\mathbbl{\Lambda}^{\dagger}|^{2}\mathbbl{\Lambda}_{I})\mathbbl{W}^{\mathsf{H}},$

where $\Sigma_{I}=\mathbbl{W}\mathbbl{\Lambda}_{I}\mathbbl{W}^{\mathsf{H}}$ is
the covariance of $I$. Hence, $\Sigma_{V}$ is unitarily diagonalized by
$\mathbbl{W}$, which concludes the proof. ∎

Remarkably, the eigenvectors of the voltage covariance matrix $\Sigma_{V}$ are
the eigenvectors of the admittance matrix $Y$ given WSS current injections. In
general, we cannot verify if $I$ is WSS since the admittance matrix is
unknown. However, in a practically relevant case when $I$ is white noise,
stationarity holds trivially, and the eigendecomposition of $\Sigma_{V}$ can
be performed to identify the eigenvectors of $Y$. On the other hand, the
eigenvalues cannot be recovered similarly since only their magnitude can be
computed from the above decomposition, i.e., from
$\mathbbl{\Lambda}_{V}=|\mathbbl{\Lambda}^{\dagger}|^{2}\mathbbl{\Lambda}_{I}$,
but not the phase.

###  IV-B Maximum a Posteriori Estimation

Upon recovering the eigenvectors of $Y$, maximum a posteriori estimation can
be leveraged to determine the eigenvalues in
$\operatorname{diag}{(\mathbbl{\Lambda})}$. The MAP estimate gives the most
likely choice of the latent variables $(V,I,Y)$ given the observations
$(\tilde{{V}},\tilde{{I}})$. To this end, the posterior distribution can be
formulated using Bayes’ rule and the conditional independence axioms:

$\displaystyle p(V,I,Y|\tilde{{V}},\tilde{{I}})$ $\displaystyle\propto
p(\tilde{{V}}|{V},{Y})p(\tilde{{I}}|{I},{Y})\frac{p({V},{I})}{p(\tilde{{V}},\tilde{{I}})}p({Y})$
$\displaystyle\,\mathrm{s.t.}\quad{I}={YV},$ (4)

where the admittance matrix is assumed to be independent of the electric
variables and their measurements. According to the measurement model in (3),
the distributions $p(\tilde{{V}}|{V},{Y})$ and $p(\tilde{{I}}|{I},{Y})$ are
Gaussian and can be expressed using the change of variables formula. For
simplicity, priors on voltages and currents are considered noninformative,
thus represented as uniform distributions over their respective domains. Under
this assumption, the quotient of priors
${p({V},{I})}/{p(\tilde{V},\tilde{I})}$ can be neglected. Finally, a prior
commonly imposed on $Y$ assumes a unit variance Gaussian distribution on all
entries of $Y$. Such prior can be represented by a matrix Gaussian
distribution $p(Y)=\exp{({-\operatorname{trace}{(YY^{\mathsf{H}}})})}$ and
leads to ridge regularization. An elaborate discussion on other practically
relevant prior distributions $p(Y)$ is given in [17]. The negative $\log$
minimization of the posterior distribution given in (IV-B) is constructed,
resulting in

$\displaystyle\min_{I,V,Y}$
$\displaystyle\,\,\,\|\tilde{{V}}-{V}\|_{2}^{2}+\|\tilde{{I}}-{I}\|_{2}^{2}+\beta\|Y\|_{\mathrm{F}}^{2}$
(5) $\displaystyle\mathrm{s.t.}$ $\displaystyle\quad{I}={YV},$

where $\beta>0$ is a constant regularization parameter proportional to the
measurement noise variance. The problem at hand is nonconvex and is
characterized by a large number of decision variables. Furthermore, previous
works [17] apply vectorization of the admittance matrix as a part of the
solution approach, which further increases the scale of the problem. We next
demonstrate how the formulation can be simplified by leveraging the obtained
spectral template $\mathbbl{W}$.

Changing the coordinates to the orthonormal basis consisting of the columns of
$\mathbbl W$, the bus voltage and current injection vectors are defined as
${\nu}\coloneqq\mathbbl{W}^{\mathsf{H}}{V}$ and
${\varphi}\coloneqq\mathbbl{W}^{\mathsf{H}}{I}$. Then, (i) the network model
(1) reduces to $\varphi=\mathbbl{\Lambda}\nu$ and the model (2) to
$\nu=\mathbbl{\Lambda}^{\dagger}\varphi$; (ii) given that the Frobenius norm
is unitarily invariant, the regularization term can be reformulated as
$\|Y\|_{\mathrm{F}}=\|\mathbbl{W}^{\mathsf{H}}Y\mathbbl{W}\|_{\mathrm{F}}=\|{\lambda}\|_{2}$,
where $\lambda=\operatorname{diag}(\mathbbl{\Lambda})$; (iii) the complex
power loss222An interesting physical interpretation can be given to the
voltage and current representations based on this property. By taking the
complex magnitude of the loss, we obtain
$|S_{\mathrm{loss}}|=\sum_{i=1}^{n}|\lambda_{i}^{*}||\nu_{i}|^{2}$. The
components in $\nu$ weighted by higher magnitude eigenvalues contribute more
to the power loss magnitude. Thus, the magnitudes of entries in $\nu$ inform
the complex power loss-efficiency of a steady-state operating point. in the
network is given by

$\displaystyle S_{\mathrm{loss}}$ $\displaystyle\ =\
\operatorname{trace}\left(VI^{\mathsf{H}}\right)\ =\
\operatorname{trace}\left(\mathbbl{W}\left(\nu\nu^{\sf
H}\right)\mathbbl{\Lambda}^{\sf H}\mathbbl{W}^{\sf H}\right)$ $\displaystyle\
=\ \operatorname{trace}\left(\left(\nu\nu^{\sf H}\right)\mathbbl{\Lambda}^{\sf
H}\right)\ =\
\operatorname{trace}\left(\mathbbl{\Lambda}^{{\dagger}}\left(\varphi\varphi^{\sf
H}\right)\right),$

that is, $S_{\mathrm{loss}}=\ \sum_{i=1}^{n}\lambda_{i}^{*}|\nu_{i}|^{2}=\
\sum_{i=1}^{n}\lambda_{i}^{\dagger}|\varphi_{i}|^{2}$.

The optimization in (5) is equivalent to the following problem in voltages
$\nu$, currents $\varphi$, and eigenvalues $\lambda$ of $Y$:

$\displaystyle\min_{{\varphi},{\nu},{\lambda}}$
$\displaystyle\,\,\,\|\tilde{{V}}-{\mathbbl{W}{\nu}}\|_{2}^{2}+\|\tilde{{I}}-{\mathbbl{W}{\varphi}}\|_{2}^{2}+\beta\|\lambda\|_{2}^{2}$
(6a) $\displaystyle\mathrm{s.t.}$
$\displaystyle\quad\varphi=\mathbbl{\Lambda}\nu,$ (6b)

with a convex objective function and bilinear constraints
$\mathbbl{\Lambda}\nu$ enforcing the network model. The classical algorithm
for solving this problem is the alternate block coordinate descent [27], which
alternates between setting $({\varphi},{\nu})$ constant and solving for
${\lambda}$, and vice versa, until convergence. Therefore, the two
optimization problems that need to be iteratively solved are given by

$\displaystyle\hat{\lambda}$
$\displaystyle=\operatorname*{arg\,min}\,\|\varphi-\operatorname{diag}{(\nu)}\lambda\|_{2}^{2}+\beta\|\lambda\|_{2}^{2},$
(7) $\displaystyle\hat{\nu}$
$\displaystyle=\operatorname*{arg\,min}\|\tilde{{V}}-{\mathbbl{W}{\nu}}\|_{2}^{2}+\|\tilde{{I}}-{\mathbbl{W}\mathbbl{\Lambda}{\nu}}\|_{2}^{2},$
(8)

where (8) is a quadratic program and (7) is a regularized least-squares
problem, both solvable in closed-form as

$\displaystyle\hat{\lambda}$
$\displaystyle=\big{(}\beta\mathcal{I}_{n}+\operatorname{diag}{(|\nu|)}\big{)}^{-1}\operatorname{diag}{(\varphi\nu^{\mathsf{H}})},$
(9) $\displaystyle\hat{\nu}$
$\displaystyle=\big{(}\mathcal{I}_{n}+\mathbbl{\Lambda}^{2}\big{)}^{-1}(\mathbbl{W}^{\mathsf{H}}\tilde{V}+\mathbbl{\Lambda}\mathbbl{W}^{\mathsf{H}}\tilde{I}).$
(10)

The matrices subject to inversion are guaranteed to be invertible since both
$\lambda^{2}$ and $|\nu|$ are nonnegative and $\beta>0$.

##  V Wiener Filter-based Network Identification

In this section, we consider a more general approach when the assumptions of
current injection stationarity (Def. 2) and constant $x/r$ ratio (Assumption
3) used in the previous section may not hold. To this end, let us define
${Z}\coloneqq({{I}},{{V}})$, with the corresponding joint covariance matrix
given by

${\Sigma_{Z}}=\begin{bmatrix}{\Sigma_{{I}}}&{\Sigma_{{I}{V}}}\\\
{\Sigma_{{I}{V}}}^{\mathsf{H}}&{\Sigma_{{V}}}\end{bmatrix}.$ (11)

In linear minimum mean square error estimation, the aim is to estimate ${I}$
from ${V}$ using a filter ${Y}$ such that the estimate ${Y}{V}$ minimizes the
mean square error $\mathbb{E}[\|{V}{Y}-{I}\|_{2}^{2}]$. Assuming that
${\Sigma}_{{{V}}}$ is full rank, the Wiener-Hopf equation admits a simple
closed-form solution, namely

${Y_{\mathrm{W}}}\coloneqq{\Sigma_{IV}}{\Sigma_{{V}}}^{-1},$ (12)

called the Wiener filter. The corresponding minimum mean square error matrix
is the Shur complement of $\Sigma_{V}$ in the joint covariance matrix, that is
$\Sigma_{I}-\Sigma_{IV}\Sigma_{V}^{-1}\Sigma_{IV}^{\mathsf{H}}$. The quality
of the Wiener filter estimate might be degraded by the effect of additive
noise in the current injection and bus voltage measurements (3). Furthermore,
large condition numbers of ${\Sigma_{{V}}}$ are commonly encountered and might
hinder the numerical computation of ${\Sigma_{{V}}}^{-1}$. These issues are
addressed in the next subsection.

###  V-A A Well-Conditioned Wiener Filter Approximation

A square matrix is ill-conditioned if it is invertible but becomes singular
for a small perturbation of some of its entries. More formally, given a normal
matrix $A\in\mathbb{C}^{n\times n}$ the condition number
$\kappa(A)=|\lambda_{\mathrm{max}}(A)|/|\lambda_{\mathrm{min}}(A)|$ is the
ratio of its largest eigenvalue $\lambda_{\mathrm{max}}(A)$ to its smallest
eigenvalue $\lambda_{\mathrm{min}}(A)$ by moduli. If $\kappa(A)$ is high, $A$
is said to be ill-conditioned. We begin our analysis by illustrating the
physical nature of the conditioning issue in network identification.

Figure 1: A 3-bus example demonstrating the physical origin of poor
conditioning of the covariance matrix $\Sigma_{V}$.

######  Example 1.

Let us consider a simple 3-bus example in Fig. 1. Without loss of generality,
we set $V_{1}$ to $1\,\mathrm{p.u.}$ The covariance matrix $\Sigma_{V}$ is
obtained by averaging the outer products of the form below over a large number
of samples:

$VV^{\mathsf{H}}=\mathbbl{1}_{3}\mathbbl{1}_{3}^{\top}+\begin{bmatrix}0&\Delta
v_{12}^{*}&\Delta v_{13}^{*}\\\ \Delta v_{12}&\Delta v_{12}^{*}+\Delta
v_{12}&\Delta v_{13}^{*}+\Delta v_{12}\\\ \Delta v_{13}&\Delta
v_{12}^{*}+\Delta v_{13}&\Delta v_{13}^{*}+\Delta v_{13}\end{bmatrix},$

where $\Delta v_{1j}=z_{1j}\Delta i_{1j},\forall j\in\\{2,3\\}$. The cross-
product terms of voltage drops are neglected. The matrix is close to
singularity in two cases: (i) if $\Delta v_{12}\approx\Delta v_{13}$ or (ii)
if the voltage drops $\Delta v_{12},\Delta v_{13}\ll 1$ are close to the
machine precision. The first condition occurs for similarly loaded lines, and
the second in light loading conditions.

Following [28], we say that a matrix is $L$-well-conditioned if it can be
computed without any inverse larger than $L\times L$. A well-conditioned
Wiener filter solution can be established by truncating the smallest
eigenvalues of the joint covariance matrix. Furthermore, it is well known that
discarding the smallest eigenvalues and the corresponding eigenvectors leads
to denoising. The eigendecomposition of the joint covariance matrix can be
performed to obtain

${\Sigma_{Z}}=\begin{bmatrix}{X_{I}}\\\
{X_{V}}\end{bmatrix}{S_{Z}}\begin{bmatrix}{X_{I}}\\\
{X_{V}}\end{bmatrix}^{\mathsf{H}},$ (13)

with ${X_{Z}}\coloneqq({X_{I}},{X_{V}})$ and eigenvalues ordered from largest
to smallest. Now let us partition the eigenvector matrices into an $n\times L$
and an $n\times M$ matrix such that

${X_{I}}=\begin{bmatrix}{X}_{{I},L}&{X}_{{I},M}\end{bmatrix},\,\,{X_{V}}=\begin{bmatrix}{X}_{{V},L}&{X}_{{V},M}\end{bmatrix},$

with $L+M=2n$ and $L\leq n$. Similarly,
$S_{Z}=\operatorname{blkdiag}(S_{Z,L},S_{Z,M})$ is partitioned into two square
diagonal matrices of sizes $L\times L$ and $M\times M$. Furthermore, we define
the Karhunen-Loève transform of ${Z}$ by
${K_{Z}}\coloneqq{X_{Z}}^{\mathsf{H}}{Z}$, from where ${Z}={X_{Z}}{K_{Z}}$. In
terms of subvectors we have ${{I}}={X_{I}}{K_{Z}}$ and ${{V}}={X_{V}}{K_{Z}}$.
Now let ${K}_{{Z},L}$ be the top $L$ submatrix of ${K_{Z}}$ so that
${{I}}\approx{X}_{{I},L}{K}_{{Z},L}$ and ${{V}}\approx{X}_{{V},L}{K}_{{Z},L}$.
A least squares approximation of the transform is given by
${K}_{{Z},L}\approx({X}_{{V},L}^{\mathsf{H}}{X}_{{V},L})^{-1}{X}_{{V},L}^{\mathsf{H}}{{V}}$.
Using this estimate, we can obtain a simple approximate filter:

${Y_{\mathrm{WCWF}}}={X}_{{I},L}({X}_{{V},L}^{\mathsf{H}}{X}_{{V},L})^{-1}{X}_{{V},L}^{\mathsf{H}},$
(14)

which is well-conditioned, i.e., the matrix inverses are $L\times L$.

######  Lemma 1.

Let
$\rho_{L}\coloneqq\operatorname{trace}{(S_{Z})-\operatorname{trace}{(S_{Z,L})}}$
define the truncation power loss. The derived $Y_{\mathrm{WCWF}}$ filter
converges to the Wiener filter $Y_{\mathrm{W}}$ as $\rho_{L}\rightarrow 0$.

The proof of the lemma above is available in [28]. The preceding lemma
demonstrates that the obtained filter corresponds to the Wiener filter in the
limit, despite being well-conditioned while the Wiener filter may not be.

###  V-B Embedding the Laplacian Structure via Postfiltering

As discussed in Sec. II, the admittance matrix is symmetric, and additionally,
if the shunt admittances are neglected or nonexistent, $Y$ has zero row-sums.
These properties are not guaranteed to hold for the Wiener filter estimate
(12) or its well-conditioned counterpart (14). In this section, we derive a
simple postfiltering procedure that can be used to enforce the Laplacian
matrix structure.

Since $Y$ has a known structure, some entries are redundant in the sense that
they can be deduced from this structure. Firstly, the admittance matrix is
symmetric, thus requiring solely $n_{d}\coloneqq\frac{1}{2}n(n+1)$ elements to
be stored in a vector $\operatorname{vech}{(Y)}$ such that
$\operatorname{vec}{(Y)}=D\operatorname{vech}{(Y)}$, where
$D\in\\{0,1\\}^{n^{2}\times n_{d}}$ is a full rank matrix called the
duplication matrix and $\operatorname{vec}{(Y)}$ is the column vector stacking
the columns of $Y$. Furthermore, the diagonal elements are redundant as they
can be expressed as a negative sum of the off-diagonal elements in each row.
That is, $\operatorname{vech}{(Y)}=R\operatorname{vechrs}{(Y)}$, where
$R\in\\{-1,0,1\\}^{n_{d}\times n_{r}}$, with
$n_{r}\coloneqq\frac{1}{2}n(n-1)$, is also full rank and
$\operatorname{vechrs}{(Y)}$ collects the off-diagonal elements.

######  Problem 2.

Let us assume that $\bar{Y}$ is an admittance matrix estimate obtained via (6)
or (14). We consider a problem of determining $\hat{{Y}}$ that is symmetric,
has zero row-sums, and is closest to $\bar{Y}$ in the Frobenius norm sense,
i.e.,

$\displaystyle\hat{{Y}}=$
$\displaystyle\operatorname*{arg\,min}_{{Y}\in\mathbb{C}^{n\times
n}}\|{\bar{Y}}-{Y}\|_{\mathrm{F}}^{2}$ (15)
$\displaystyle\quad\mathrm{s.t.}\quad{Y}={Y}^{\top},{Y}_{ii}=-\sum_{j\neq
i}{Y}_{ij},\forall i.$

Interestingly, a closed-form solution based on pseudoinverses of $D$ and $R$
can be obtained.

######  Proposition 2.

The offdiagonal entries of $\hat{Y}$, the solution to (15), are given by
$\operatorname{vechrs}{(\hat{{Y}})}=R^{\dagger}{D}^{\dagger}\operatorname{vec}{(\bar{Y})}$.

###### Proof.

By applying the $\operatorname{vec}{(\cdot)}$ operator on the objective
function we obtain
$\|\bar{Y}-Y\|_{\mathrm{F}}^{2}=\|\operatorname{vec}{(\bar{Y})}-\operatorname{vec}{(Y)}\|_{2}^{2}$.
The constraints can be included by considering that
$\operatorname{vec}{(Y)}=DR\operatorname{vechrs}{(Y)}$. Upon substituting the
previous equality in the objective function, the least squares solution is
given by
$\operatorname{vechrs}{(\hat{Y})}=(DR)^{\dagger}\operatorname{vec}{(\bar{Y})}=R^{\dagger}D^{\dagger}\operatorname{vec}{(\bar{Y})}$.
The last equality holds as $R$ and $D$ are full rank. ∎

Therefore, applying $R^{\dagger}D^{\dagger}$ as a postfilter to an obtained
$Y$ estimate enforces the Laplacian structure. We note that the pseudoinverses
can be constructed efficiently as their structure is generic [29], and depends
only on the number of buses in the network. To further motivate and justify
use of the postfilter, let us consider the following problem.

######  Problem 3.

Using the available measurements $\mathbbl{\tilde{I}}$ and
$\mathbbl{\tilde{V}}$ a constrained least squares network identification
problem respecting the Laplacian structure of $Y$ is formulated as

$\displaystyle\hat{{Y}}=$
$\displaystyle\operatorname*{arg\,min}_{{Y}\in\mathbb{C}^{n\times
n}}\|\mathbbl{\tilde{I}}-Y\mathbbl{\tilde{V}}\|_{\mathrm{F}}^{2}$ (16)
$\displaystyle\quad\mathrm{s.t.}\,\,{Y}={Y}^{\top},{Y}_{ii}=-\sum_{j\neq
i}{Y}_{ij},\forall i.$

Solution to the constrained least squares problem is given in the proposition
below. The proof is omitted as it resembles the proof of Proposition 2, with
the additional identity
$\operatorname{vec}{(Y\mathbbl{\tilde{V}})}=(\mathbbl{\tilde{V}}^{\top}\otimes\mathcal{I}_{n})\operatorname{vec}{(Y)}$
required.

######  Proposition 3.

The off-diagonal entries of the solution to (16) are given by

$\displaystyle\operatorname{vechrs}{(\hat{{Y}})}$
$\displaystyle={R}^{\dagger}{D}^{\dagger}\,(\mathbbl{\tilde{V}}^{\top}\otimes\mathcal{I}_{n})^{\dagger}\operatorname{vec}{(\mathbbl{\tilde{I}})}$
(17)
$\displaystyle={R}^{\dagger}{D}^{\dagger}\operatorname{vec}{(\hat{{Y}}_{\mathrm{LS}})},$
(18)

where
$\hat{{Y}}_{LS}\coloneqq(\mathbbl{\tilde{V}}^{\top}\otimes\mathcal{I}_{n})^{\dagger}\operatorname{vec}{(\mathbbl{\tilde{I}})}$
denotes the unconstrained least squares solution.

The result above shows that the optimal solution to the constrained least
squares problem (16) can be obtained by applying postfiltering (15) to the
unconstrained least squares solution. The equivalence between the
unconstrained least squares solution and the Wiener filter solution (12) when
the same measurements $\mathbbl{\tilde{I}}$ and $\mathbbl{\tilde{V}}$ are used
to compute the sample covariances motivates the use of the postfiltering for
ensuring the Laplacian structure in our proposed solution (14). Note that the
least squares problem in (16) was previously considered in [11]. However, the
connection between the unconstrained and constrained least squares problems
via the postfilter in (15) was not recognized.

Figure 2: Single line diagram of the IEEE 33-bus test system. Figure 3:
Comparison of the relative error of existing and proposed methods for various
noise levels.

##  VI Results

The proposed identification methods are evaluated on the IEEE 33-bus network
presented in Fig. 2. We assume that a $\mu$PMU device is placed on each node
in the network, measuring both voltage and current phasors. The procedure to
generate the estimation data follows [17], where synthetic load profiles are
created using the GENETX generator, and the power flow procedure is run using
the PandaPower library. To realistically represent $\mu$PMU measurements, the
voltage and current phasors are corrupted with $0.01\%$ standard deviation
Gaussian noise. Measurements collected at 50Hz frequency are averaged over a
minute and 7 days of thus constructed data ($10080$ tuples of voltage and
current phasors) are used in the considered estimation procedures.

###  VI-A Estimation Performance under Varying Noise Levels

In this section, we compare the estimation accuracy of the proposed methods:
(i) identification procedure in (6), labeled by $\mathrm{MAP}_{\lambda}$
henceforth, and (ii) the well-conditioned Wiener filter (14), denoted as
$\mathrm{WCWF}$ hereafter, to the state-of-the-art approaches from the
literature under varying noise levels. More precisely, we use Ordinary Least
Squares (OLS) [11], Lasso [13], and MAP with $\ell_{2}$ regularization [17]
for benchmarking. The metric used to evaluate the accuracy of an estimation
procedure is the relative Frobenius norm
$\varepsilon_{\mathrm{F}}=\|\hat{Y}-Y\|_{\mathrm{F}}/\|Y\|_{\mathrm{F}}$,
where $Y$ is the true admittance matrix and $\hat{Y}$ denotes the estimate.
The results in Fig. 3 show significant estimation bias when performing (6) due
to the violation of stationarity and constant $x/r$ ratio assumptions.
Furthermore, the non-errors in variables models, OLS and Lasso, demonstrate
high sensitivity to the measurement noise. Finally, the $\ell_{2}$ regularized
MAP from [17] and the proposed well-conditioned Wiener filter (14) demonstrate
satisfactory performance over a large range of noise levels.

TABLE I: Accuracy and computation time of the considered estimation methods on
three distribution grid test cases.

Network | Lasso height 2ex |  $\mathrm{MAP}_{\ell_{2}}$ height 2ex | WCWF
---|---|---|---
$\varepsilon_{\mathrm{F}}[\%]$ | $\tau[s]$ | $\varepsilon_{\mathrm{F}}[\%]$ | $\tau[s]$ | $\varepsilon_{\mathrm{F}}[\%]$ | $\tau[s]$
CIGRE10 | $51.8$ | $0.5$ | $10.8$ | $33$ | $1.9$ | $0.002$
IEEE33 | $59.5$ | $14.3$ | $2.9$ | $780.5$ | $2.1$ | $0.06$
IEEE123 | $60.1$ | $125.3$ | $6.36$ | $3527$ | $4.2$ | $0.71$

TABLE II: The condition number of matrices subjected to inversion in the least
squares and Wiener filter methods.

Network | $\kappa(\Sigma_{V})$ | $\kappa({\Sigma_{I}})$ | $\kappa({X}_{{V},L}^{\mathsf{H}}{X}_{{V},L})$
---|---|---|---
CIGRE10 | $2\times 10^{12}$ | $3\times 10^{10}$ | $1\times 10^{7}$
IEEE33 | $7\times 10^{12}$ | $4\times 10^{11}$ | $1\times 10^{7}$
IEEE123 | $1\times 10^{13}$ | $3\times 10^{12}$ | $7\times 10^{6}$

###  VI-B Computational Efficiency and Conditioning Analysis

To analyze the computational efficiency, we perform parameter estimation on
three benchmark distribution grids: the 10-bus CIGRE MV feeder, the previously
considered IEEE 33-bus system, and the three-phase part of the IEEE 123-bus
system consisting of 56 buses. Table I summarizes the estimation results in
the form of accuracy $\varepsilon_{\mathrm{F}}$ and computation time $\tau$.
Only Lasso, MAP with $\ell_{2}$ regularization, and the well-conditioned
Wiener filter are considered for brevity. The nominal Gaussian noise of
$0.01\%$ is used. The table shows that WCWF outperforms the other methods in
terms of accuracy and computation time across all tested scenarios.

Table II presents the condition numbers of three matrices, namely
$\Sigma_{V}$, $\Sigma_{I}$, and ${X}_{{V},L}^{\mathsf{H}}{X}_{{V},L}$, which
are subjected to inversion in different identification methods. Specifically,
$\Sigma_{V}$ is inverted in both the OLS approach and the Wiener filter method
in (12), $\Sigma_{I}$ is commonly inverted in the impedance matrix estimation
process [10], and ${X}_{{V},L}^{\mathsf{H}}{X}_{{V},L}$ with $L=n$ is inverted
in the proposed well-conditioned Wiener filter approach (14). The results
indicate that the proposed approach offers significant improvement in
conditioning compared to the other methods across all three test cases.

###  VI-C MAP Estimation under Stationary Current Injections

The MAP estimation of $\lambda$ has thus far demonstrated unsatisfactory
performance – see Fig. 3. However, in the previous simulation setting, neither
the stationarity nor the constant $x/r$ ratio assumptions were valid. We now
enforce the two assumptions and analyze how the performance degrades when
deviations from these assumptions are imposed. To this end, a modified version
of the IEEE 33-bus system is created by setting
$b_{ij}=-r_{ij},\forall(i,j)\in\mathcal{E}$, thus achieving a constant
conductance-susceptance ratio of one throughout the network. We introduce
operator $\mathrm{ndiag}(\cdot):\mathbb{C}^{n\times
n}\rightarrow\mathbb{C}^{n\times n}$, which converts the diagonal entries of a
matrix to zeros and keeps the off-diagonal elements. To measure the deviation
of a matrix $A\in\mathbb{C}^{n\times n}$ from being unitarily diagonalizable
by $\mathbbl{W}$ we define the relative distance
$\mathrm{dist}{W}(A)=\|\mathrm{ndiag}(\mathbbl{W}^{\mathsf{H}}A\mathbbl{W})\|_{\mathrm{F}}/\|\mathbbl{W}^{\mathsf{H}}A\mathbbl{W}\|_{\mathrm{F}}$.

The results in Fig. 4 show the increase in estimation error with the increase
in $\mathrm{dist}{W}(\Sigma_{I})$ which quantifies the violation of the
stationarity property. The estimator is demonstrated to be unbiased when
applied to a constant $x/r$ ratio network. However, a significant bias of
approximately $10\%$ error is present when estimating the original IEEE 33-bus
network. A constant estimation error is shown for the original network
corresponding to $\Sigma_{I}=\mathcal{I}_{n}$ since $\mathbbl{W}$ is not well-
defined for non-normal admittance matrices.

Figure 4: Dependency of the estimation error to the deviation of the current
injection covariance matrix to a matrix diagonalizable by $\mathbbl{W}$.

##  VII Conclusion

This paper shows how a maximum a posteriori admittance matrix estimation can
be simplified when the current injections are stationary. Nevertheless, the
approach only performed well for estimating networks with a constant
conductance-susceptance ratio. We have subsequently derived a more general and
practical admittance matrix estimation approach based on linear minimum mean
square error estimation. Our results demonstrate that the proposed approach is
more accurate and computationally efficient than the state-of-the-art when
applied to standard test networks.

## References

  * [1] A. M. Kettner and M. Paolone, “On the properties of the compound nodal admittance matrix of polyphase power systems,” _IEEE Transactions on Power Systems_ , vol. 34, no. 1, pp. 444–453, 2019.
  * [2] F. Dörfler, J. W. Simpson-Porco, and F. Bullo, “Electrical networks and algebraic graph theory: Models, properties, and applications,” _Proceedings of the IEEE_ , vol. 106, no. 5, pp. 977–1005, 2018.
  * [3] D. Deka _et al._ , “Learning distribution grid topologies: A tutorial,” _IEEE Transactions on Smart Grid_ , pp. 1–1, 2023.
  * [4] A. von Meier, E. Stewart, A. McEachern, M. Andersen, and L. Mehrmanesh, “Precision micro-synchrophasors for distribution systems: A summary of applications,” _IEEE Transactions on Smart Grid_ , vol. 8, no. 6, pp. 2926–2936, 2017.
  * [5] P. Schirmer and I. Mporas, “Non-intrusive load monitoring: A review,” _IEEE Transactions on Smart Grid_ , vol. 14, no. 1, pp. 769–784, 2023.
  * [6] S. Bolognani, N. Bof, D. Michelotti, R. Muraro, and L. Schenato, “Identification of power distribution network topology via voltage correlation analysis,” in _52nd IEEE Conference on Decision and Control_ , 2013, pp. 1659–1664.
  * [7] D. Deka, S. Backhaus, and M. Chertkov, “Structure learning in power distribution networks,” _IEEE Transactions on Control of Network Systems_ , vol. 5, no. 3, pp. 1061–1074, 2018.
  * [8] G. Cavraro and V. Kekatos, “Inverter probing for power distribution network topology processing,” _IEEE Transactions on Control of Network Systems_ , vol. 6, no. 3, pp. 980–992, 2019.
  * [9] A. Wehenkel, A. Mukhopadhyay, J.-Y. L. Boudec, and M. Paolone, “Parameter estimation of three-phase untransposed short transmission lines from synchrophasor measurements,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 69, no. 9, pp. 6143–6154, 2020.
  * [10] K. Moffat _et al._ , “Unsupervised impedance and topology estimation of distribution networks—limitations and tools,” _IEEE Transactions on Smart Grid_ , vol. 11, no. 1, pp. 846–856, 2020.
  * [11] Y. Yuan, S. H. Low, O. Ardakanian, and C. J. Tomlin, “Inverse power flow problem,” _IEEE Transactions on Control of Network Systems_ , vol. 10, no. 1, pp. 261–273, 2023.
  * [12] E. Fabbiani, P. Nahata, G. De Nicolao, and G. Ferrari-Trecate, “Identification of AC distribution networks with recursive least squares and optimal design of experiment,” _IEEE Transactions on Control Systems Technology_ , vol. 30, no. 4, pp. 1750–1757, 2022.
  * [13] O. Ardakanian, V. W. S. Wong, R. Dobbe, S. H. Low, A. von Meier, C. J. Tomlin, and Y. Yuan, “On identification of distribution grids,” _IEEE Transactions on Control of Network Systems_ , vol. 6, no. 3, pp. 950–960, 2019\.
  * [14] J. Yu, Y. Weng, and R. Rajagopal, “PaToPa: A data-driven parameter and topology joint estimation framework in distribution grids,” _IEEE Transactions on Power Systems_ , vol. 33, no. 4, pp. 4335–4347, 2018.
  * [15] R. K. Gupta, F. Sossan, J. Y. Le Boudec, and M. Paolone, “Compound admittance matrix estimation of three-phase untransposed power distribution grids using synchrophasor measurements,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 70, pp. 1–13, 2021.
  * [16] J.-S. Brouillon _et al._ , “Bayesian methods for the identification of distribution networks,” in _Proceedings of the IEEE Conference on Decision and Control_ , 2021, pp. 3646–3651.
  * [17] ——, “Bayesian error-in-variables models for the identification of distribution grids,” _IEEE Transactions on Smart Grid_ , vol. 14, no. 2, pp. 1289–1299, 2023.
  * [18] A. Sandryhaila and J. M. F. Moura, “Discrete signal processing on graphs,” _IEEE Transactions on Signal Processing_ , vol. 61, no. 7, pp. 1644–1656, 2013.
  * [19] N. Perraudin and P. Vandergheynst, “Stationary signal processing on graphs,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 13, pp. 3462–3477, 2017.
  * [20] A. G. Marques, S. Segarra, G. Leus, and A. Ribeiro, “Stationary graph processes and spectral estimation,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 22, pp. 5911–5926, 2017.
  * [21] S. H. Low, _Power System Analysis: A Mathematical Approach_. Lecture Notes, Caltech, 2023.
  * [22] D. Turizo and D. K. Molzahn, “Invertibility conditions for the admittance matrices of balanced power systems,” _IEEE Transactions on Power Systems_ , vol. 38, no. 4, pp. 3841–3853, 2023.
  * [23] F. Dörfler and F. Bullo, “Kron reduction of graphs with applications to electrical networks,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , vol. 60, no. 1, pp. 150–163, 2013.
  * [24] X. Du, A. Engelmann, Y. Jiang, T. Faulwasser, and B. Houska, “Optimal experiment design for AC power systems admittance estimation,” _IFAC-PapersOnLine_ , vol. 53, no. 2, pp. 13 311–13 316, 2020.
  * [25] A. C. Varghese, A. Pal, and G. Dasarathy, “Transmission line parameter estimation under non-gaussian measurement noise,” _IEEE Transactions on Power Systems_ , vol. 38, no. 4, pp. 3147–3162, 2023.
  * [26] M. Bariya, D. Deka, and A. Von Meier, “Guaranteed phase topology identification in three phase distribution grids,” _IEEE Transactions on Smart Grid_ , vol. 12, no. 4, pp. 3605–3612, 2021.
  * [27] H. Zhu, G. Leus, and G. B. Giannakis, “Sparsity-cognizant total least-squares for perturbed compressive sampling,” _IEEE Transactions on Signal Processing_ , vol. 59, no. 5, pp. 2002–2016, 2011.
  * [28] E. K. P. Chong, “Well-conditioned linear minimum mean square error estimation,” _IEEE Control Systems Letters_ , vol. 6, pp. 2431–2436, 2022\.
  * [29] J. R. Magnus and H. Neudecker, “The elimination matrix: Some lemmas and applications,” _SIAM Journal on Algebraic Discrete Methods_ , vol. 1, no. 4, pp. 422–449, 1980.