# The Effect of the Russian-Ukrainian Conflict from the Perspective of
Internet eXchanges

1st Cristian Trusin  University of Twente
The Netherlands
<EMAIL_ADDRESS>2nd Leandro Bertholdo  University of Twente
The Netherlands
<EMAIL_ADDRESS>3rd José Jair Santanna  Northwave and University of
Twente
The Netherlands
<EMAIL_ADDRESS>
###### Abstract

In 2022 the Russian invasion of Ukraine began. It is known that Ukraine faced
outages because of the damage to their infrastructure. It is also known that
Russia was boycotted by the international community. However the impact on the
telecommunications of the two countries remains unknown. In this paper we
quantified the degree to which the Internet was affected in both countries by
analyzing routing tables from five large Internet Exchange Points (IXPs). IXPs
provide a central point of interconnection where internet traffic can be
freely exchanged between Autonomous Systems (ASes). This centrality makes IXPs
a good vantage point for analyzing changes in the Internet infrastructure.
With data collected before and after the start of the conflict we observed
considerable damage to the Ukrainian Internet network with numerous outages
and minimal damage to the Russian network. An average of 11.12% of Ukrainian
ASes were unreachable at each IXP. We identified the biggest outages and the
events responsible for them. This paper highlights resilience issues during
conflicts to the network and management community, and serves as a basis for
future more in-depth research.

###### Index Terms:

Internet Exchange Point, IXP, Russia, Ukraine, War, Conflict, Impact,
Autonomous System, AS, Border Gateway Protocol, BGP.

##  I Introduction

On the 24th of February 2022 Russia invaded Ukraine. Ukraine has been
experiencing blackouts as parts of its infrastructure have been cut off after
bombings [1]. Russia has been boycotted by the international community and
businesses. Sanctions were registered since 24th of February 2022 and are
still ongoing [2]. Major transit providers like Cogent and Lumen, stopped
selling Internet access to Internet Service Providers (ISPs) inside Russia [2,
3]. Furthermore, content providers like TikTok and Netflix have ceased their
services in the country [4]. These events have had an impact on the Internet
in both countries, some services are unreachable or are forced to be routed
through a slower path, but it remains unclear how big this impact is.

In this paper we investigate some of this impact using data from Internet
Exchange Points (IXPs). IXPs are a crucial element of today’s internet [5],
they provide direct peering opportunities between Autonomous Systems (ASes).
Direct peering is a voluntary interconnection of two separate networks with
the purpose of directly exchanging the traffic between the users of those
networks, it improves latency and decreases interconnection costs. Since IXPs
connect many ASes, they are a rich source for networking research [5], and are
suitable to analyze changes in the network, thus they can help us quantify the
effects of the war on the Internet in Ukraine and Russia.

Our methodology to assess this impact is based on Border Gateway Protocol
(BGP) data from the routing tables provided by five large IXPs: AMSIX, LINX,
SIX, AUIX, and SPOIXBR. We only investigated open peering data from these IXPs
and we analyze the period between 19th of February 2022 (before Russian
invasion) and 29th of April 2022 (more than two months after the invasion).
Our goal is to quantify the effect that the war had on the Russian and
Ukrainian networks from the perspective of the IXPs.

To achieve this goal, we have defined the following research questions (RQs)
as a basis for our research:

  * •

RQ1: How many Russian and Ukrainian ASes became unreachable?

  * •

RQ2: For how long were these ASes unreachable?

  * •

RQ3: What were the biggest outages and what events are responsible for them?

This paper is structured as follows: In Section 2 we discuss the necessary
background information and the related work. Section 3 will talk about the
methodology we used to answer the Research Questions. Section 4 will discuss
the results of each research question and their implications. In Section 5 we
conclude our work and discuss possibilities for future work.

Figure 1: Number of Ukrainian Prefix announcements before the start of the
war.

##  II Background

The Internet is a network of networks, it can be broken into smaller networks
called Autonomous Systems (ASes). Each AS is a group of routers that belong to
the same organization. For example, to ISPs, universities, tech companies,
governments etc. Each AS wishing to connect to the Internet and exchange
routing information must have a valid AS Number (ASN) assigned by the Internet
Assigned Numbers Authority.

ASes communicate with each other using the Border Gateway Protocol (BGP). BGP
provides a standardized way of sharing routing information between ASes, it is
the glue that connects all the Autonomous Systems together. At its essence,
BGP is trying to connect the whole Internet together by letting ASes share
their routing data. Then from this routing data, it selects the best path to a
destination based on an attribute called as_path. The as_path is a sequence of
ASNs that a packet will follow to reach the destination. BGP chooses the
shortest as_paths. The details of this process are not important for this
paper, but more information can be found in [6].

For two ASes to interconnect, they need to establish a peering session. One
way of doing this is to connect to an IXP, which is a central point of
interconnection for ASes, and form an open peering agreement where they can
freely exchange traffic with other members of the IXP. This is achieved by a
Route Server. A Route Server is a construct in an IXP which has information
about the routes offered by other member ASes. By connecting to a Route
Server, you can establish connection to multiple ASes by using a single BGP
session, this enables one-to-many peerings instead of traditional one-to-one
peerings [7], which improves scalability and makes routing easier. By
connecting to only 10 IXPs you can reach 56% of all ASes on the Internet [8],
which makes them a good source of insight to answer our research questions. A
good high-level introduction to IXPs can be found in Gerson’s et al. [9], a
more detailed dive into IXPs can be found in [5].

##  III Methodology

In this section we will explain how we performed our research. We will start
by explaining what this dataset is and how it was collected, and then dive
into our implementation. We quantified the outages and changes in the network
by performing Data Analysis on selected BGP attributes from our collected
dataset. The code for this implementation as well as a Database required for
it is provided here [10].

###  III-A Datasets and IXP Selection

From 13th of February 2021 to 29th of April 2022, the routing tables from
AMSIX, LINX, SIX, AUIX, SPOIXBR have been collected daily through a Looking
Glass, as part of the research conducted by Bertholdo et al. [8]. A Looking
Glass allows an IXP to share information from the Route Server. The routing
tables represent only open peering data.

The 5 IXPs listed above were chosen due to their size and representativeness
[8], [11]. AMSIX (Amsterdam Internet Exchange) and LINX (London Internet
Exchange) are in Europe, close to the conflict. LINX is particularly
interesting since they reported sanctions against Russia [12]. AMSIX is the
second biggest IXP in Europe after DE-CIX (Deutscher Commercial Internet
Exchange). We found some gaps in the DE-CIX dataset during war days, so it was
not possible to include it in our analysis. We also included other IXPs
outside Europe to verify if problems were observed on other continents. SIX
(Seattle Internet Exchange) and AUIX (Australia Internet Exchange) were used
to represent North America and Australia in this research. And finally for
South America we had data from SPOIXBR (S$\tilde{a}$o Paulo Internet Exchange
in Brazil).

The gathered routing tables are in a .csv format and represent BGP data. We
have a routing table for every single day and for each IXP. So, for any given
day we have 5 routing tables to analyze since we have 5 IXPs. We will focus on
the as_origin, as_neighbor and as_path attributes of the routing tables.
as_origin denotes the AS that can be reached from this IXP, it is an ASN. If
this AS cannot be reached anymore, its as_origin will disappear from the
routing table. One as_origin can route traffic to multiple prefixes.
as_neighbor from the perspective of an IXP denotes the ASN of the AS that is
directly connected to the IXP. as_path denotes all the ASes the as_origin had
to transit through to reach the IXP. as_path is a list of ASNs, the last entry
in the list is always the as_origin. We will use these attributes to analyze
if and how routes to Ukrainian and Russian ASes have changed.

We know for a fact that certain ASes experienced outages because of the war
[1], [4]. Moreover, the routing paths to some Russian ASes have changed [3].
To quantify these, we needed a baseline. We analyzed several days of data
before the war as can be seen in Figure 1, and concluded that except for 2
outages in September and October 2021, the data showed little variation, so we
chose 19th of February 2022 to serve as our baseline. The last datapoint we
have is on 29th of April (End label in Figure 1).

###  III-B Implementation

The first phase of the research was identifying which ASNs belong to Russia
and Ukraine. We did this by creating a database where we aggregated data from
Regional Internet Registries (RIRs). This data offers information about the
current allocations and assignments of Internet Number Resources (IPs and
ASNs). We used the data from all 5 RIRs: AFRINIC, APNIC, ARIN, LANIC and RIPE
NCC. We filtered it to contain ASN assignments only, and for each ASN we kept
the country to which it belongs. The data from each RIR was later aggregated
together, resulting in a Database that can identify the originating country of
a given ASN. From now on we will refer to this Database as ASNDB. ASNDB allows
us to filter our routing tables and isolate Ukrainian and Russian routes.

####  III-B1 Counting the number of unreachable ASes (RQ1)

To answer RQ1 we needed to count the number of ASes which became unreachable
at any point after the start of the war. We joined ASNDB with our routing
tables, resulting in a table where for each route we also have its originating
country code as a new column. By filtering on this new column, we were able to
find the Ukrainian and Russian as_origins for each IXP. By analyzing our
baseline and looking at our last datapoint, we identified which and how many
as_origins became unreachable after the start of the war for each IXP. This
was achieved by selecting all distinct prewar as_origins and searching for
those as_origins in the last routing table we have collected. If an as_origin
was not found, it was pronounced unreachable. It could be that a temporary
outage was happening during the collection of our last datapoint, making some
as_origins temporarily unreachable. To avoid drawing false conclusions we also
checked if they were unreachable during the 3 days before the collection of
our last datapoint. The results did not show any significant differences.

####  III-B2 Counting the number of days ASes were unreachable (RQ2)

Due to the sheer amount of ASes and variation in the number of days they were
offline for, we could not just count for how long each AS was offline, since
this data would be hard to visualize. We decided to focus on the as_neighbors
instead. An as_neighbor, in the context of an IXP, is the ASN of the AS which
is directly connected to the IXP. This decision was made because if one
as_neighbor were to disconnect, all the prefixes it provided routing
information for will also be disconnected, thus making them unreachable for
this IXP and the member ASes of this IXP. This is indicative of an outage at
the said AS, since the consequences of disconnecting are undesirable. If an
ISP were to disconnect from an IXP, all its clients would lose access to fast
routes connecting them to the global internet. Bearing this in mind, we
counted the number of distinct as_neighbors for each day at every IXP.

####  III-B3 Identifying the outages and the responsible events (RQ3)

To answer RQ3 we needed to find when the outages happened. We could not use
the plot from RQ2 for this purpose, since it only offers a high-level view of
the situation. It could be that the as_neighbor has not disconnected from the
IXP, but it has fewer prefixes to route to because of the damage to the
infrastructure.

Instead, we calculated the number of Ukrainian and Russian prefix
announcements, alongside the number of distinct Ukrainian and Russian
as_origins for every day, for each IXP. This was achieved by joining the ASNDB
with the routing tables and filtering on routes originating from the two
countries. We plotted this data to see where we can find dips in the plot
(outages). Having found a dip we would look at the news for that day alongside
the report from NetBlocks [1] to find events that are responsible for the
outage.

Figure 2: Number of Ukrainian as_neighbors at AMSIX. Figure 3: Number of
Ukrainian Prefix announcements at AMSIX and LINX. Figure 4: Number of
Ukrainian Prefix announcements at SIX, AUIX and SPOIXBR. Figure 5: Number of
Russian Connections at AMSIX. Figure 6: Number of Distinct Russian Prefixes at
LINX. Figure 7: Number of Distinct Russian as_origins at SIX, AUIX and
SPOIXBR.

##  IV Results and Discussion

In this section we show and discuss our findings. Since there is little
overlap between the situation of Ukrainian and Russian Networks, we split it
into 2 subsections: one discussing our results regarding Ukraine and another
regarding Russia. The plots used show data since 19th of February, since we
chose this date as our baseline of the prewar situation. The datapoints before
this day showed little variation as discussed in Section III Part A, so they
were cut off from the resulting plot.

###  IV-A Mapping outages in Ukraine

The Ukrainian network took considerable damage since the start of the war.
Table I provides information about the number of unreachable Ukrainian ASes at
each IXP. On average every single IXP lost 11.12% of Ukrainian ASes since the
war started. The European IXPs were the most affected of all.

TABLE I: Number of unreachable Ukrainian as_origins IXP | Total ASes | Lost ASes | % Lost
---|---|---|---
AUIX | 1016 | 87 | 8.5%
LINX | 1335 | 254 | 19.0%
AMSIX | 1571 | 164 | 10.4%
SPOIXBR | 1021 | 92 | 9.0%
SIX | 1096 | 96 | 8.7%

####  IV-A1 Europe

AMSIX experienced the most variation in available Ukrainian connections. Each
as_neighbor provides routes to many as_origins, losing one is significant. At
AMSIX, a Ukrainian as_neighbor provides routes to an average of 400 prefixes
which is around 409.600 IP Addresses. As we can see in the historical graph in
Figure 2, most of these ASes experienced an outage for a day, however some (AS
25133 – McLaut-Invest) were unreachable for 50 days, and some never
reconnected (AS 12963 – Volz). A peculiar finding is Omegatelecom (AS 199995)
which is a new member of AMSIX coming from Ukraine. It connected on April 9th
and is responsible for providing routes to 243 distinct Ukrainian prefixes and
90.264 IP addresses. LINX has not experienced any significant loss of
as_neighbors. The only disconnection happened on 28th until 29th of March
during Event 3 in which LINX lost its only Ukrainian as_neighbor which
provided routes to many Ukrainian prefixes. This is explained by a cyberattack
on Ukrtelecom’s core infrastructure as reported by BBC and NetBlocks [13, 14]
and is also responsible for the decrease in prefix announcements during that
time.

As can be seen in Figure 3, there are significant variations in the number of
prefix announcements at AMSIX and LINX since the start of the war. Event 1
shows a significant loss of connectivity, around 1000 less prefixes were
announced that day. This can be attributed to a major blackout in the region
of Sumy, which is reported as the largest disruption to Telcom services since
the start of the war [15]. As a result of bombing at the local thermal power
plant, the region of Sumy experienced a large blackout which also affected
Internet connectivity. Connectivity started to improve over the course of the
next days, and 10 days later new prefix announcements were observed, even more
than before the war.

Other major disruptions can be seen at Event 2 and Event 3. We could not
identify a concrete event responsible for the former, but the latter can be
attributed to the cyberattack on Ukrtelecom’s core infrastructure mentioned
above [13, 14]. Overall Ukrainian prefixes at AMSIX seem to have rebound to
prewar levels since April.

However, LINX does not show such a rebound. It sustained the biggest Ukrainian
losses, registering a total of 254 out of 1335 ASes. 19% of all Ukrainian
networks were disconnected from this IXP. Figure 3 shows that LINX has been
experiencing a gradual decline in Ukrainian connections since the beginning of
the war. Besides experiencing the same outages as AMSIX, LINX also had a loss
of connectivity during Event 4, an emergency was reported by the operator WNET
which impacted Ukraine’s international connectivity during that period [16].

####  IV-A2 Other continents

Figure 4 shows the data analyzed for SIX (USA), AUIX (Australia) and SPOIXBR
(Brazil). These IXPs follow a similar pattern to Europe, with a significant
loss of connectivity at Event 1, due to the Sumy Blackout [15], unidentified
outage at Event 2 and cyberattack at Event 3 [14]. Another dip in connections
can be seen during Event 5. This can possibly be explained by a global
Internet Disruption in mid-March, when the number of Internet outages
increased by 22% [17]. None of these IXPs have Ukrainian as_neighbors, thus
none of them were lost during the war.

####  IV-A3 Observations and Discussion

Ukraine has sustained damage to its networks since the start of the war, an
average of 11.12% of Ukrainian ASes are now unreachable. We have found 5
outages: Event 1 - The Sumy Oblast Blackout, Event 2 - Unidentified Outage,
Event 3 - Ukrtelecom Cyberattack, Event 4 - WNET International Connectivity
emergency, Event 5 - March Global Internet Disruption. The first three of
these can be observed both in and outside Europe. Event 4 is exclusive to
LINX, and Event 5 is exclusive to IXPs outside Europe. Except for LINX, the
number of Ukrainian connections at IXPs has started recovering and the number
of outages has decreased significantly since April.

LINX had the biggest losses, we observed a decrease of 19% in the number of
reachable Ukrainian as_origins. Judging by the fact that Event 3 marks the
biggest outage for LINX with a loss of 2000 prefix announcements, and that
LINX temporarily lost its only Ukrainian as_neighbor during that event, we can
assume that these losses are attributed to the disconnection of that neighbor.
This as_neighbor provided routes to most Ukrainian prefixes at LINX and
because of the war, it likely lost connection to a lot of ASes, so it could no
longer provide routes to those ASes for LINX.

AMSIX has 12 Ukrainian as_neighbors which provided more redundancy. If one
neighbor lost connection to some ASes, another one might still have a way to
reach them, which explains why AMSIX restored most connections. LINX did not
lose all Ukrainian routes due to the disconnection of the said as_neighbor
because neighboring countries kept providing alternative routes to some
Ukrainian ASes.

These observations invite a discussion about the resiliency of the Ukrainian
network. While damaged, it remains functional. This can be attributed to the
presence of many redundancies in many layers of the network [18]. Ukraine has
19 local IXPs and they share the market almost evenly. This provides
alternative routes in case of outages i.e. redundancy. Moreover, Ukrainian
ISPs are very decentralized, 55% of the ISPs serve less than 1% of the market
[18], meaning that the Ukrainian population is connected through a variety of
ISPs without a single point of failure. However we have observed that energy
and cyberattacks are a weak-point of the Ukrainian network. As seen in Event
1, damaging energy sources has severe effects on the connectivity of the
country since the local ISPs and infrastructure need energy to operate. Event
3 showed the damage a cyberattack can inflict, however judging by the lack of
further outages caused by cyberattacks, it is safe to say that the Ukrainian
network learned to mitigate them. Perhaps it could be attributed to the help
that Ukraine got from companies like Microsoft and ESET [19].

###  IV-B The effects of the boycott on Russia’s connectivity to IXPs

Russia was heavily boycotted for their invasion by most countries in the world
[1]. There have been threats that Russia will be cut off from the global
Internet. Moreover, major Internet Transit Providers like Cogent and Lumen
have stopped offering Russia their services [2], [3]. In this section we
analyze how the Russian network was affected at the IXPs.

Table II provides information about the number of Russian ASes that lost
access to the IXP since the start the of the war. An average of 10.94% ASes
were lost. However, new ASes have connected during that same period, thus
increasing the total number of Russian as_origins and prefixes since the start
of the war.

TABLE II: Number of unreachable Russian as_origins IXP | Total ASes | Lost ASes | % Lost
---|---|---|---
AUIX | 3749 | 117 | 3.1%
LINX | 2886 | 109 | 3.7%
AMSIX | 421 | 62 | 14.7%
SPOIXBR | 415 | 78 | 18.7%
SIX | 419 | 61 | 14.5%

####  IV-B1 Europe

As can be seen in Figure 5, AMSIX shows no decrease in Russian prefixes
pnnouncements and the number of distinct prefixes. A sudden increase in prefix
Announcements was observed for 6 days (Fig. 5 Spike 1). We analyzed the number
of distinct Russian prefixes at AMSIX during this period and found an increase
of about 800 prefixes. Another such spike can be seen (Fig. 5 Spike 2),
however this time it was only followed by an increase of 200 distinct
prefixes, and it is still ongoing. There was no variation in the number of
Russian as_neighbors at AMSIX.

LINX tells a similar story to AMSIX, there was no decrease in the number of
Russian routes and as_neighbors. As can be seen in Figure 6, an increase of
distinct prefixes was observed from 15th of April. These findings are peculiar
since LINX was one of the few IXPs that announced sanctions against Russia,
and one would expect the number of Russian connections to decrease. After
further investigation we found that LINX only disconnected two Russian ASes:
Rostelecom (AS 12389) and Megafon (AS 31133) [12]. These ASes are Russia’s
biggest ISPs, but apparently, they had little to no traffic at LINX and were
planning on abandoning LINX services themselves [20].

####  IV-B2 Other Continents

Figure 7 shows the data for SIX, AUIX and SPOIXBR. They all show that Russia
experienced an outage on 12th until 14th of March (Fig. 7 Outage 1), which is
consistent with the outage observed at these IXPs for Ukraine (Fig 4 Event 5)
in the same period. It is also explained by the mid-March global Internet
Disruptions [17]. Another serious outage was observed on 31st of March (Fig. 7
Outage 2). SIX and AUIX had no Russian as_neighbors prewar, and the situation
remained unchanged. SPOIXBR had one Russian as_neighbor and besides a small
1-day outage on 2nd of March it remained stable. These IXPs also show an
increase in Russian prefixes in April, as well as an increase in distinct
as_origins since the start of the war.

####  IV-B3 Observations and Discussion

Russia has not sustained much damage to their network, contrary to the rumors
circling in press that it will be disconnected from the global Internet. We
found a decrease of 10.94% in the number of distinct as_origins since the
start of the war. However, as we observed in the plots the number of distinct
as_origins has increased after the war, making this loss less significant
since the total number of Russian ASes has grown since the start of the war.
The loss of the old ASes could just signify a change in Russia’s internal
Internet network infrastructure. Moreover, the number of Russian prefixes has
increased, suggesting that the sanctions applied to Russia did not have a
major impact on their network. Our findings align with the ones seen in the
report by ThousandEyes [3].

Major transit providers Cogent and Lumen have stopped their services to
Russia. This did not have any significant effects since we did not register a
drop in overseas Russian prefixes at SIX, AUIX and SPOIXBR. This can be
explained by the numerous alternative transit providers that Russia can still
use [21]. The reason behind the ban was limiting Russia’s capacity to perform
cyberattacks, according to CNN [22], and not disconnecting Russia from the
Internet.

IXPs have not ceased providing services to Russia in any meaningful way. LINX
has stopped providing services to the two biggest Russian ISPs: Rostelecom and
Megafon. However, these ISPs were barely using LINX’s services already, and
LINX continues to advertise prefixes which contain those ISPs in their
as_path. According to the statement provided by LINX to their customers, these
measures were not designed to block Russia from the Internet but rather to
target individuals which are affiliated with the current Russian government
and responsible for the war [12]. Some of them are in direct ownership or
benefit from Rostelecom and Megafon. Europe imposed sanctions against these
persons [23] and LINX was merely complying with them.

##  V Conclusions and Future Work

In this paper we quantified the effect of the conflict between Russia and
Ukraine from the perspective of five large IXPs. We found how many Ukrainian
and Russian ASes were lost, for how long, and we identified the biggest
outages and the events responsible for them.

Our results show that each IXP has lost connection to an average of 11.12% of
Ukrainian ASes. We identified 5 big outages which are responsible for this
loss. LINX lost most Ukrainian connections and has not yet restored all of
them. Other IXPs started recovering since April and the number of outages has
decreased. The Ukrainian network shows resilience, and while clearly damaged
remains functional due to the presence of many redundancies.

We did not find any substantial damage, or loss of connectivity to the Russian
network. Contrary to the sanctions imposed, the Russian network remains
reachable and shows signs of growth.

As this research focused on a limited set of IXPs, future work could analyze
data from other IXPs. In particular DE-CIX, since it is the largest one in
Europe and from the limited data we had, we noticed that it had a lot of
Russian and Ukrainian connections, so it might provide some new insights. It
would also be worthwhile to look at the data we analyzed from a different
angle, by counting the number of unique reachable IP addresses, thus leading
to more fine-grained results. We would also like to analyze changes in the
as_paths, and investigate the impact the war had on other countries, in
particular neighboring ex-soviet republics, since they may rely on routes
provided by Ukraine and Russia.

## Acknowledgment

This work is funded by the Netherlands Organisation for Scientific Research,
PAADDoS project (628.001.029), and CONCORDIA, the Cybersecurity Competence
Network supported by the European Union's Horizon 2020 research and innovation
programme under grant agreement No 830927. We would like to thank Suzan Bayhan
for her continued support for this paper.

## References

  * [1] “Internet disruptions registered as russia moves in on ukraine - netblocks,” accessed: 07-05-2022. [Online]. Available: https://netblocks.org/reports/internet-disruptions-registered-as-russia-moves-in-on-ukraine-W80p4k8K
  * [2] “2022 boycott of russia and belarus - wikipedia,” accessed: 01-07-2022. [Online]. Available: https://en.wikipedia.org/wiki/2022_boycott_of_Russia_and_Belarus
  * [3] “Russia and the global internet - thousandeyes,” accessed: 07-05-2022. [Online]. Available: https://www.thousandeyes.com/blog/russia-global-internet
  * [4] “Russia, blocked from the global internet, plunges into digital isolation - the new york times,” accessed: 07-05-2022. [Online]. Available: https://www.nytimes.com/2022/03/07/technology/russia-ukraine-internet-isolation.html
  * [5] N. Chatzis, G. Smaragdakis, and A. Feldmann, “On the importance of internet exchange points for today’s internet ecosystem,” _arXiv preprint arXiv:1307.5264_ , 2013.
  * [6] S. Halabi, _Internet Routing Architectures_ , 2nd ed. Cisco Press, 8 2000.
  * [7] P. Richter, G. Smaragdakis, A. Feldmann, N. Chatzis, J. Boettger, and W. Willinger, “Peering at peerings: On the role of ixp route servers,” in _Proceedings of the 2014 Conference on Internet Measurement Conference_ , 2014, pp. 31–44.
  * [8] L. M. Bertholdo, J. M. Ceron, L. Z. Granville, and R. van Rijswijk-Deij, “Forecasting the impact of ixp outages using anycast.” in _TMA_ , 2021.
  * [9] P. S. Ryan PhD and J. Gerson, “A primer on internet exchange points for policymakers and non-engineers,” _Available at SSRN 2128103_ , 2012.
  * [10] C. Trusin, “Code repository for this research,” accessed: 04-09-2022. [Online]. Available: https://github.com/Bo0ne/UkraineResearch
  * [11] Euro-IX, “The ixp database,” https://ixpdb.euro-ix.net/en/, 09 2022, (Accessed on 09/15/2022).
  * [12] “Bill woodcock on twitter, march 11 2022,” accessed: 22-06-2022. [Online]. Available: https://t.co/Wytcj098LJ
  * [13] “Ukraine war: Major internet provider suffers cyber-attack - bbc news,” accessed: 02-07-2022. [Online]. Available: https://www.bbc.com/news/60854881
  * [14] “Netblocks on twitter, march 28 2022,” accessed: 23-06-2022. [Online]. Available: https://twitter.com/netblocks/status/1508465391244304389?
s=20&t=lNIGIXV-kAkEttykeSDVCw

  * [15] “Netblocks on twitter, march 3 2022,” accessed: 23-06-2022. [Online]. Available: https://t.co/zNG42jfbIS
  * [16] “Netblocks on twitter, april 8 2022,” accessed: 23-06-2022. [Online]. Available: https://t.co/x6gPCYKyt0
  * [17] “Global internet outages jump 22%, business-hours disruptions fall in mid-march — s&p global market intelligence,” accessed: 23-06-2022. [Online]. Available: https://www.spglobal.com/marketintelligence/en/news-insights/latest-news-headlines/global-internet-outages-jump-22-business-hours-disruptions-fall-in-mid-march-69455283
  * [18] E. Aben, “The internet in ukraine,” https://ripe84.ripe.net/wp-content/uploads/presentations/23-ukraine-internet.emileaben.ripe84.pdf, 05 2022, accessed: 7-9-2022.
  * [19] Sophos, “Six months on: Looking back at the role of cyberattacks in the ukraine war,” https://news.sophos.com/en-us/2022/09/14/six-months-on-looking-back-at-the-role-of-cyberattacks-in-the-ukraine-war/, 09 2022, accessed: 13-9-2022.
  * [20] AK&M Information Agency, “Rostelecom and megafon were disconnected from the linx internet traffic exchange point, which is located in london,” accessed: 25-06-2022. [Online]. Available: https://www.akm.ru/eng/press/rostelecom-and-megafon-were-disconnected-from-the-linx-internet-traffic-exchange-point-which-is-loca/
  * [21] “Cogent and lumen curtail operations in russia — kentik blog,” accessed: 25-06-2022. [Online]. Available: https://www.kentik.com/blog/cogent-disconnects-from-russia
  * [22] “Why internet backbone services cogent and lumen are cutting off russia - cnn,” accessed: 02-07-2022. [Online]. Available: https://edition.cnn.com/2022/03/11/tech/russia-internet-backbone-cogent-lumen/index.html
  * [23] “Eu sanctions against russia following the invasion of ukraine,” accessed: 24-09-2022. [Online]. Available: https://eu-solidarity-ukraine.ec.europa.eu/eu-sanctions-against-russia-following-invasion-ukraine_en


IEEEexample:BSTcontrol

# Minimally Invasive Live Tissue High-fidelity Thermophysical Modeling using
Real-time Thermography

Hamza El-Kebir  Student Member, IEEE  Junren Ran  Yongseok Lee  Leonardo P.
Chamorro  Martin Ostoja-Starzewski  Richard Berlin  Gabriela M. Aguiluz
Cornejo  Enrico Benedetti  Pier C. Giulianotti  and Joseph Bentsman  Senior
Member, IEEE Manuscript submitted for review on July 21, 2022. Research
reported in this publication was supported by the National Institute of
Biomedical Imaging and Bioengineering of the National Institutes of Health
under award number R01EB029766. The content is solely the responsibility of
the authors and does not necessarily represent the official views of the
National Institutes of Health. Copyright (c) 2022 IEEE. Personal use of this
material is permitted. However, permission to use this material for any other
purposes must be obtained from the IEEE by sending an email to pubs-
<EMAIL_ADDRESS>H. El-Kebir is with the Dept. of Aerospace Engineering,
University of Illinois Urbana-Champaign, Urbana, IL 61801 USA (e-mail:
elkebir2@illinois.edu).J. Ran, Y. Lee, L. P. Chamorro, M. Ostoja-Starzweski,
and J. Bentsman are with the Dept. of Mechanical Science and Engineering,
University of Illinois Urbana-Champaign, Urbana, IL 61801 USA (email: {jran2,
yl50, lpchamo, martinos, jbentsma}@illinois.edu). J. Bentsman is the
corresponding author.R. Berlin is with the Department of Trauma Surgery, Carle
Hospital and Department of Biomedical and Translational sciences, Carle
Illinois College of Medicine, Urbana, IL 61801 USA.G. M. Aguiluz Cornejo, E.
Benedetti, and P. C. Giulianotti are with the Dept. of Surgery, University of
Illinois at Chicago, Chicago, IL 60612 USA.

###### Abstract

We present a novel thermodynamic parameter estimation framework for energy-
based surgery on live tissue, with direct applications to tissue
characterization during electrosurgery. This framework addresses the problem
of estimating tissue-specific thermodynamics in real-time, which would enable
accurate prediction of thermal damage impact to the tissue and damage-
conscious planning of electrosurgical procedures. Our approach provides basic
thermodynamic information such as thermal diffusivity, and also allows for
obtaining the thermal relaxation time and a model of the heat source, yielding
in real-time a controlled hyperbolic thermodynamics model. The latter accounts
for the finite thermal propagation time necessary for modeling of the
electrosurgical action, in which the probe motion speed often surpasses the
speed of thermal propagation in the tissue operated on. Our approach relies
solely on thermographer feedback and a knowledge of the power level and
position of the electrosurgical pencil, imposing only very minor adjustments
to normal electrosurgery to obtain a high-fidelity model of the tissue-probe
interaction. Our method is minimally invasive and can be performed _in situ_.
We apply our method first to simulated data based on porcine muscle tissue to
verify its accuracy and then to _in vivo_ liver tissue, and compare the
results with those from the literature. This comparison shows that
parameterizing the Maxwell–Cattaneo model through the framework proposed
yields a noticeably higher fidelity real-time adaptable representation of the
thermodynamic tissue response to the electrosurgical impact than currently
available. A discussion on the differences between the live and the dead
tissue thermodynamics is also provided.

Biomedical infrared imaging, thermography, tissue thermodynamics, real-time
model estimation.

##  1 Introduction

Figure 1: Illustration of the proposed parameter estimation for controlled
tissue thermodynamics (in this case, electrosurgery is used as an example). In
step 1, a probing sequence is determined, consisting of (i) a fixed power
stationary probing task, and (ii) a moving cutting task, possibly with
variable velocity and power. These probing sequences are applied to the tissue
in step 2, during which the tissue thermal response is recorded by a
thermographer. These thermographer images are then processed in real-time in
step 3, where _attention-based noise-robust averaging_ (ANRA) is used to
obtain: (A) the thermal diffusivity, (B) the heat source model, and (C) the
thermal relaxation time. Subtasks (A) and (B) rely on the stationary probe
data, while subtask (C) is based on the moving probe data. This yields an
estimated Maxwell–Cattaneo hyperbolic heat transfer model.

In this work, we evaluate the suitability of the hyperbolic Maxwell–Cattaneo
model of heat propagation to support a high-fidelity thermodynamic
representation of the live tissue response to electrosurgical impact, with
applications to porcine liver tissue. As demonstrated in [1, 2], the movement
speeds encountered in electrosurgical procedures (5–20 mm/s), when compared to
the speed of heat propagation in tissue (around 2 mm/s), shows that wave-like
phenomena due to the finite propagation of thermal information play a dominant
role in the overall thermodynamics. An important effect of the slow thermal
propagation speed of tissue is the fact that the electrosurgical heat source
undergoes distortion from an ideally Gaussian shape to a tear-drop-like shape
[2], causing significant discrepancies between simulations that use the
Maxwell–Cattaneo formulation versus the classical Fourier heat equation in the
case of a moving heat source [3]. Since accurately predicting the thermal
spread in live tissues is of prime importance in enabling damage-conscious
control of the electrosurgical process and preventing the high heat flux from
violating safety margins [4, 5], we focus in this work on developing a method
for estimating the relevant thermophysical parameters for this model based on
thermographer readings.

Previous approaches to parameter estimation in service of modeling the
electrosurgical process have been undertaken by [6, 7], focusing exclusively
on estimating the thermal diffusivity, as well as [1], where the thermal
relaxation time was the subject of study. These methods rely on direct
solution of the underlying PDEs (partial differential equations), to find a
thermal diffusivity that minimizes the error between the simulation result at
specific locations that correspond to locations of embedded thermal sensors,
and the temperatures recorded at the sensor locations. This method is
invasive, requiring embedding sensors into the tissue. In addition, obtaining
parameter estimates requires a significant amount of time because the
governing PDE must be explicitly solved multiple times to obtain an optimal
parameter estimate, precluding the use of such an algorithm for real-time
adaptation and safety region mapping during clinical procedures on live
tissue. More importantly, multivariate optimization often fails to yield
accurate parameter combinations owing to the existence of local minima
encountered by optimization algorithms [8]. These local minima often
necessitate the use of prior knowledge about the parameters to constrain the
search space [6, 7], thereby inadvertently forcing preconceived biases on the
obtained parameter estimates. Our proposed method does not require such
assumptions, mainly by a proposed probing sequence that isolates various
governing phenomena of the underlying heat transfer problem, allowing
parameters to be obtained in series, rather than through the commonly observed
error-prone method of simultaneous parameter fitting.

In contrast to the above methods, in [8] we presented a method that relies
fully on non-invasive remote infrared thermography, in which a thermal camera
is used to capture the surface temperature of the tissue. Our approach does
not rely on explicitly solving the governing PDEs, but instead employs noise-
robust gradient operators to _directly_ obtain the PDE’s governing parameters
by balancing the underlying model differential equations based on the
thermographic data. In this work, we build upon these results to achieve a
method capable of automatic model estimation for a hyperbolic thermal
equation, particularly the Maxwell–Cattaneo heat equation. Central to this
endeavor is the adoption of a simple minimally invasive probing sequence,
which eliminates ambiguity in the parameter estimation step by isolating the
various driving phenomena, namely: diffusion, wave propagation, and external
forcing. Our approach allows for _in situ_ parameter estimation with a
_minimally invasive probing sequence_ (obviating biopsies and laboratory
testing), with results being obtained in _real-time_. These thermophysical
parameters are estimated such that our parsimonious thermodynamics model, the
Maxwell–Cattaneo heat equation, best fits the observed thermodynamics. Any
unmodeled dynamics are subsumed in this limited set of parameters, producing
_effective_ parameter estimates that can directly be used for prediction and
control of the tissue thermodynamics. To the best of the authors’ knowledge,
this combination of probing sequence design and direct parameter estimation
based on thermographer readings has hitherto not been attempted in the
literature. Our method is outlined in Fig. 1.

This paper makes the following contributions to the existing literature: (i)
development of a novel procedure for parameterizing the tissue model in real
time under electrosurgical impact, based on thermographic data; (ii) using the
latter, enabling a higher-fidelity than currently available modeling of the
thermal tissue behavior though parametrizing the Maxwell–Cattaneo tissue
model, directly applicable during regular operation without invasive sensor
insertion; (iii) based on (ii), offering a high fidelity comparison of tissue
thermodynamics under _in vivo_ vs _ex vivo_ electrosurgical procedures.

This paper is organized as follows. In Sec. 2, the Maxwell–Cattaneo heat
transfer model is developed and the parameter estimation problem that will be
considered is introduced. Sec. 3 presents the parameter estimation method,
providing justification for our design choices and illustrating considerations
in composing probing sequences to separate thermodynamic phenomena. In Sec. 4,
simulation results demonstrating the efficacy of our method compared to a
known ground truth are used for verification. In Sec. 5, we apply our method
to real-life thermographic data from electrosurgery on _in vivo_ porcine liver
tissue, discussing differences between the two, and showing improvements in
the simulation results based on the updated parameters. Conclusions and future
research directions are presented in Sec. 6.

##  2 Problem Formulation

As mentioned in the introduction, in this work we seek to obtain a heat
transfer model that captures the salient dynamics observed in energy-based
surgery methods applied to live tissue. In our model, we wish to include the
following phenomena:

  1. i)

Diffusion (heat conduction);

  2. ii)

Presence of a controlled external heat source;

  3. iii)

Finite heat propagation time;

The first of these phenomena, diffusion, is by far the most studied in the
field of heat transfer in live tissue [9, 6, 7, 10] (the former two works
studying porcine aorta, and the latter considering general animal tissues).
The most commonly used formulation is the classical _Fourier heat equation_ ,
which models thermal conductivity. When operating on a one-dimensional domain,
the following expression is obtained based on the thermal diffusivity $a$:

$\dot{x}(t,\eta)=a\frac{\partial^{2}x(t,\eta)}{\partial\eta^{2}},$ (1)

where $x$ denotes the temperature, $t$ the time, and $\eta$ the spatial
location.

After adding the heat source, we obtain the following formulation:

$\dot{x}(t,\eta)=a\frac{\partial^{2}x(t,\eta)}{\partial\eta^{2}}+T_{p(t)}u(t)q(\eta),$
(2)

In this expression, the heat source is position at location $p(t)$, where
$T_{p(t)}$ is an operator that serves to translate the center of the heat
source to $p(t)$. The applied heat source is modeled as a unit heat source
field $q$ scaled by the input power $u(t)$.

A major drawback in the classical Fourier heat equation, is the fact that
changes in the heat source propagate instantaneously throughout the entire
domain [2, 3]. This assumption is unphysical, as heat waves propagate at a
finite speed through any physical medium [1, 2, 3]. An important consequence
of an infinite heat propagation speed is a typical assumption that the tissue
immediately reacts to external thermal excitation, whereas in practice, when
an external heat source is applied there is a finite time before the
temperature starts rising [1]. In addition, in the case of electrosurgery, one
does not observe a rise of temperature in front of the moving needle, even
during regular surgical procedures [1]; this also points at the importance of
a finite heat propagation time. Introduction of the latter is often achieved
by introducing a second time derivative term, so as to induce wave-like
behavior. The formulation thus obtained is commonly referred to as the
_Maxwell–Cattaneo equation_ , or the _telegraph equation_ :

$\tau\ddot{x}(t,\eta)+\dot{x}(t,\eta)=a\frac{\partial^{2}x(t,\eta)}{\partial\eta^{2}}+T_{p(t)}u(t)q(\eta).$
(3)

Here, the second time derivative is scaled by $\tau$, which is known as the
_thermal relaxation time_. A so called _second speed of sound_ is often
obtained by considering both the thermal diffusivity and thermal relaxation
time as follows [2]:

$c=\sqrt{a/\tau}.$ (4)

We wish to obtain parameter estimates on a patient-by-patient basis,
irrespective of underlying conditions. In light of this goal, we seek a method
that is:

  1. 1.

minimally invasive (i.e., no biopsies or unnecessary tissue resection is
required);

  2. 2.

applicable _in situ_ , without the need for external specialized equipment or
laboratory testing;

  3. 3.

runnable in real-time, with results immediately available.

This yields the following general problem statement, which can be classified
as a _parameter estimation problem_ :

###### Problem 1

Given $N$ surface thermographer readings $\\{\hat{x}^{(i)}\\}_{i=1}^{N}$ taken
on subsets of $\hat{\Omega}$ with fixed time interval $\Delta t>0$, as well as
control inputs $(u(t),v(t))$, estimate the following parameters of system (5):

  1. 1.

Thermal diffusivity $a$;

  2. 2.

Unit heat source $q$;

  3. 3.

Thermal relaxation time $\tau$.

We now proceed by briefly introducing a well-posed mathematical formulation of
our model on which our results in this work will be based.

###  2.1 Mathematical Heat Transfer Model

In this work, we consider the following thermodynamic model on a two-
dimensional compact domain $\Omega\subseteq\mathbb{R}^{2}$:

$\begin{bmatrix}\dot{x}(t)\\\ \ddot{x}(t)\\\
\dot{p}(t)\end{bmatrix}=\begin{bmatrix}0&1&0\\\
\frac{a}{\tau}\nabla^{2}&-\frac{1}{\tau}&0\\\
0&0&0\end{bmatrix}\begin{bmatrix}x(t)\\\ \dot{x}(t)\\\
p(t)\end{bmatrix}+\begin{bmatrix}0\\\ \frac{1}{\tau}T_{p(t)}u(t)q\\\
v(t)\end{bmatrix},$ (5)

where $a>0$ denotes the thermal diffusivity, $\tau>0$ denotes the thermal
relaxation time, and $\nabla^{2}$ is the Laplacian operator. $x(t)$ denotes
the temperature field at time $t$, such that $x(t)\in H^{1}(\Omega)$, where
$H^{1}(\Omega)$ refers to the Sobolev space $W^{1,2}(\Omega)$; by $x(t,\eta)$
we denote $x(t)(\eta)$ for $\eta\in\Omega$. The heat source position is
denoted by $p(t)\in\Omega$, which is used as a parameter in the translation
operator $T$, defined as

$\begin{split}T_{p}:H^{1}(\Omega)&\to H^{1}(\Omega),\quad\forall
p\in\Omega,\\\
T_{p}q(\eta)&\mapsto\chi_{\Omega}(\eta)q(\eta-p),\quad\forall\eta\in\Omega,q\in
H^{1}(\Omega),\end{split}$ (6)

where $\chi_{\Omega}$ is the indicator function on the set $\Omega$. In our
model, $q\in H^{1}(\Omega)$ denotes a normalized heat source term, which is
scaled by control input $u(t)\in\mathbb{R}_{+}$, which indicates the
electrical power provided by the electrosurgical power supply. Finally, $v\in
C^{0}(\mathbb{R}^{n})$ is the known controlled heat source velocity, which is
such that $p(t):=p(0)+\int_{0}^{t}v(\tau)\ \mathrm{d}\tau\in\Omega$ for all
$t\in[0,\infty)$. We denote the full state by $z\in H^{1}(\Omega)\times
L^{2}(\Omega)\times C^{1}(\Omega)$, which includes the temperature, the
temperature rate of change, and the probe position.

In this work, we assume having access to thermographer data comprised of
temperatures that are sampled on a finite Cartesian grid with uniform spacing
$\Delta\eta>0$, denoted by $\hat{\Omega}\subset\Omega$. We consider
$\mathrm{dim}\ \hat{\Omega}=M<\infty$, such that
$\hat{\Omega}:=\\{\eta^{(i)}\\}_{i=1}^{M}$.

In the remainder of this work, we make the following assumptions:

  1. (A1)

The initial condition $x(0)\in H^{1}(\Omega)$ has known constant temperature,
i.e., $x(0,\eta)=x_{0}\chi_{\Omega}(\eta)$ for all $\eta\in\Omega$ for some
known $x_{0}\in\mathbb{R}$.

  2. (A2)

$\inf_{t}u(t)\geq 0$ and $\sup_{t}u(t)<\infty$.

  3. (A3)

On the boundary of $\Omega$, $\partial\Omega$, a homogeneous Neumann boundary
condition holds, i.e., $\nabla x(t)|_{\partial\Omega}\equiv 0$.

  4. (A4)

Heat source velocity $v\in C_{\infty}^{0}(\mathbb{R}^{n})$ is continuous and
bounded, such that $p(t):=p(0)+\int_{0}^{t}v(\tau)\ \mathrm{d}\tau\in\Omega$
for all $t\in[0,\infty)$.

As can be seen in (5), we assume that there are no errors in the probe
position, nor the heat source model, i.e., (5) is considered to be a nominal
model of the system’s thermodynamics, which can be used in a subsequent
controller design or for predicting/assessing tissue damage. In the remainder
of this work, we consider control input pair $(u(t),v(t))$ to be given for all
$t\in[0,\infty)$; we also assume that there is no error in the control inputs.

#### Existence and Uniqueness

To prove the existence and uniqueness of (5), we may apply the Picard
iteration method on the system to obtain a unique solution ($x(t,\eta)$) for
known input pair $(u(t),v(t))$ and with any initial condition ($x_{0}(\eta)$)
that satisfies assumption (A1) and (A2) [11, §2.3]; continuity of the ordinary
differential equation (ODE) that defines $p(t)$ follows directly from (A4).
Thus, well-posedness of the heat equation under consideration is guaranteed.

###  2.2 Probing Sequence for Parameter Estimation

Rather than attempting to estimate all three parameters of system (5) at once,
it bears mentioning that, in general, there is no unique solution to Problem 1
for an arbitrary electrosurgical cutting sequence. Indeed, the ratios $a/\tau$
and $q/\tau$ do not necessarily admit unique solutions without imposing
further constraints, especially when considering standard methods that are
prone to finding only locally optimal parameters, not to mention the case of
unmodeled dynamics and noise when dealing with real-life data. To address this
problem, a central idea is to precede the main cutting task by a series of
minimally invasive probing tasks, to limit the number of phenomena
simultaneously being excited to obtain a unique set of parameters.

Let us rewrite the dynamics of $x(t)$ from (5):

$\tau\ddot{x}(t)+\dot{x}(t)=a\nabla^{2}x(t)+T_{p(t)}u(t)q.$ (7)

From observations in [2], it follows that a non-moving source with steady heat
input results in $\tau\ddot{x}(t)$ vanishing after a transient, whose duration
is dictated by thermal relaxation time $\tau$. Therefore, it is possible to
study the following equation in isolation, i.e., the tissue’s _free response_
:

$\dot{x}(t)=a\nabla^{2}x(t).$ (8)

To do so, one would need to sufficiently raise the temperature, e.g., by
sustained application of a small amount of electrosurgical action without
moving the electrosurgical probe’s contact point. Starting from a low or
steady-state temperature will lead to $\dot{x}(t)\approx 0$, which does not
allow for observing $a\nabla^{2}x(t)$. Having raised the temperature, after
disengaging the heat source, one will need to wait for the thermal wave caused
by the sudden lack of external forcing to have subsided as dictated by the
second speed of sound.

It suffices to wait for the final heat wave to propagate outside of the field
of view of the thermographer, after which the free response of the tissue can
be observed. The latter has been the object of study in [8], in which a novel
method known as _attention-based noise-robust averaging_ (ANRA) has been
developed to solve for the tissue’s thermal diffusivity $a$ given the free
response.

We can then turn our attention to the stationary forced response, in which we
study the time during which heat was actively supplied to the tissue. In this
case, provided that a constant non-moving heat source is present, after the
initial thermal wave has subsided, we observe a _forced stationary response_ ,
governed by the following dynamics:

$\dot{x}(t)=a\nabla^{2}x(t)+T_{p_{0}}u_{0}q,$ (9)

where $p_{0}$ is the fixed position of the probe, and $u_{0}$ is the constant
power supplied. In what follows, we shall discuss our model of $q$, and the
method by which we obtain it. It must be noted that $a$ is assumed to be
constant throughout the entire cut. The latter assumption is justified,
considering that during electrosurgery, undamaged tissue is constantly
contacted by the moving probe. In contrast, damaged tissue is often
immediately ablated, or is confined to a small region around the cut [10]. It
also bears mentioning that obtaining $a$ and $q$ from the free and forced
response can be achieved through a single stationary probing event, in which
the tissue is heated (as dictated by (9)) and subsequently allowed to cool
down (as governed by (8)).

Finally, we may obtain the thermal relaxation time $\tau$ by studying any form
of regular cutting action, as governed by (7). The exact method of estimating
$\tau$ will be discussed in the following section.

##  3 Parameter Estimation Method

In this section, we present our parameter estimation method to solve Problem
1. We identify three distinct steps:

  1. A)

Determining the thermal diffusivity $a$;

  2. B)

Determining the heat source $q$;

  3. C)

Determining the thermal relaxation time $\tau$.

It will become clear that these steps must be performed in order using both
the stationary and moving thermal response, as illustrated in Fig. 1.

Our method relies solely on surface thermographer readings, and thus serves to
obtain a model of the apparent thermodynamics of the tissue, as observed on a
2D surface. Extensions to a 3D domain based on 2D observations have been
discussed in related work [8], where an observer for 3D tissue temperature has
been developed; such an extension is, however, beyond the scope of the current
work.

To enable direct estimation of parameters, we leverage a set of noise-robust
gradient operators introduced by P. Holoborodko [12, 13]; these combine the
capabilities of finite difference operators with those of a low-pass filter,
allowing for noise-robust gradient computation [8]. In essence, our method
reduces the parameter estimation problem from an optimization problem over the
solution of (5) to that of balancing differential equations, as will be shown
next. Since it was previously hard to obtain gradient field from thermography
data on account of high-frequency noise, such gradient-based balancing
techniques have not existed before the concept was introduced by the authors
in [8].

In the following, we shall denote noise-robust gradient operators by
$\hat{\nabla}$ and $\hat{\partial}_{t}$. We now proceed with the parameter
estimation method for the thermal diffusivity.

###  3.1 Thermal Diffusivity Estimation

We require that (8) holds, with a sufficiently large magnitude of $\dot{x}$ to
obtain the thermal diffusivity. As mentioned in Sec. 2.2, this is achieved by
raising the temperature of the tissue by application of a small amount of
sustained electrosurgical action, after which the free response is studied.
Focusing on (8), it is apparent that solving for

$\hat{\partial}_{t}\hat{x}(t,\eta)=\hat{a}(\eta)\hat{\nabla}^{2}\hat{x}(t,\eta)$
(10)

at each $\eta\in\hat{\Omega}$, will yield a field $\hat{a}$ as opposed to a
single scalar. Many of these parameter estimates will be erroneous, safe for
those estimates obtained in regions with a high rate of temperature change
[8]. This observation has previously led to the development of _attention-
based noise-robust averaging_ (ANRA) [8], where an _attention layer_ is
constructed to allow for an average parameter to be obtained. We refer the
reader to [8] for a full explanation of this method.

###  3.2 Heat Source Estimation

We assume that the unit heat source is modeled as a Gaussian heat source of
the form [8]:

$q(\eta)=q_{0}\exp(-\|\eta\|^{2}/\rho^{2}),$ (11)

for $\eta\in\Omega$, where $q_{0}$ is such that $\int_{\Omega}q(\eta)\
\mathrm{d}\eta=1$, and $\rho$ controls the spread of heat. Also, the heat
source scales linearly with the power applied, denoted by $u$ in (5). We also
introduce some additional notation: we shall call $\ddot{x}$ the ‘temperature
acceleration,’ and $\dddot{x}$ the ‘temperature jerk.’ When applied to real-
life data, these repeated derivatives are well defined.

Given a diffusivity estimate $\hat{a}$ as obtained from the free response in
Sec. 3.1, we can now consider the forced response of the tissue, i.e., the
time when the tissue was subject to stationary and constant electrosurgical
action. In other words, we stipulate that the forced part of the tissue
response be obtained by subjecting the tissue to a constant power, with a
constant electrosurgical probe position.

To ensure that the forced stationary response model (9) is valid, it is
necessary to wait for the initial thermal waves to subside; the modeling
errors accrued when not doing so will be demonstrated in due detail in Sec. 4.
Since this time period is dependent on the thermal relaxation time $\tau$,
which has not been obtained yet, we elect to study the second time derivative
of the observed temperature field instead. In particular, we consider
$\hat{\ddot{x}}_{\text{max}}[k]:=\max_{i,j}|\hat{\partial}^{2}_{t}\hat{x}_{i,j}[k]|$;
the rationale for studying this value lies in the fact that a small magnitude
of $\ddot{x}$ results in (9) being an accurate approximation of (7). In this
work, we wait until the maximum temperature acceleration
$\hat{\ddot{x}}_{\text{max}}[k]$ is below 10% of the all-time maximum
temperature acceleration $\hat{\ddot{x}}_{\text{max,
total}}=\sup_{k=1,\ldots,N}\hat{\ddot{x}}_{\text{max}}[k]$, before assuming
that (9) is sufficiently valid; this is illustrated in Sec. 4.

Based on (9), we define
$\hat{S}[k]=\hat{\partial}_{t}\hat{x}[k]-\hat{a}\hat{\nabla}^{2}\hat{x}[k]$.

Given that the heat source is located at position $\eta_{0}\in\hat{\Omega}$
and has power $u_{0}>0$, we choose to apply nonlinear least-squares fitting to
obtain $(q_{0},\rho)$ by considering the residual:

$r_{i,j}[k]=S_{i,j}[k]-\hat{Q}\exp\left(-\|\eta_{i,j}-\eta_{0}\|^{2}/\hat{\rho}^{2}\right),$
(12)

where $(i,j)$ are chosen in a small radius around $\eta_{0}$; in this work, a
region around $\eta_{0}$ where $S$ assumes positive values is considered. This
yields estimates $\hat{q}_{0}=\hat{Q}/u_{0}$ and $\hat{\rho}$, after taking an
unweighted time-average over the period following subsidence of the initial
thermal wave, as described previously.

###  3.3 Thermal Relaxation Time Estimation

Having obtained $\hat{a}$, $\hat{q}_{0}$, and $\hat{\rho}$, we may now obtain
an estimate for the thermal relaxation time $\tau$. The latter is particularly
difficult to observe, since its effects manifest only in the case of transient
behavior, such as initial application of electrosurgical power, or when
disabling the electrosurgical probe. In steady-state cutting, as is often
observed in surgical practice, wave-like effects are rarely observed, given
that the process is in a quasi-steady-state where $\ddot{x}$ is of very small
magnitude [2]. Instead, a diffuse shock front can be observed when cutting
speeds are sufficiently high, as is commonly the case [2]; in practice, the
resulting Mach cone is so diffuse that $\ddot{x}$ is not of a significant
magnitude to be useful in estimation. This latter shock front will be the
subject of study in future work.

The equation we wish to solve for is similar to (10), but is instead governed
by (7):

$\begin{split}&-\hat{\partial}_{t}\hat{x}(t,\eta)+\hat{a}\hat{\nabla}^{2}\hat{x}(t,\eta)+T_{p(t)}u(t)\hat{q}_{0}\exp\left(-\frac{\|\eta-p(t)\|^{2}}{\hat{\rho}^{2}}\right)\\\
&=\hat{\tau}(\eta)\hat{\partial}^{2}_{t}\hat{x}(t,\eta),\end{split}$ (13)

where we have used the heat source model obtained in Sec. 3.2.

Our approach is mostly identical to that of Sec. 3.1 (see (10)), with the
exception that the variables are replaced by those appearing in (13). To
obtain a time average, instead of considering the intensity $I[k]$, we elect
to impose an event-driven method. Any time we have a sufficiently large change
in $u(t)$ (i.e., $\partial_{t}|u(t)|>\dot{u}_{\text{threshold}}>0$), we
consider the time starting at that change, up until the time when the third
time-derivative of the temperature,
$\hat{\dddot{x}}_{\text{max}}[k]:=\max_{i,j}|\hat{\partial}^{3}_{t}\hat{x}_{i,j}[k]|$,
is a fraction of the highest recorded value (in this work, 10% is considered).
This rate of acceleration is also known as _jerk_ , a commonly encountered
concept in mechanics. The latter condition based on the maximum absolute jerk
$\hat{\dddot{x}}_{\text{max}}[k]$ provides a clear means of delineating the
time during which the effects of a sudden change are apparent; this is
illustrated in Sec. 4.

##  4 Simulation

To verify the efficacy of the approach outlined in Sec. 3, we simulated (5)
based on the properties of porcine liver tissue. These parameters have been
obtained from the literature [1, 2, 9]. They are as follows: a thermal
diffusivity $a$ of $0.144$ mm2/s, a thermal relaxation time $\tau$ of 0.1 s,
and a heat source spread $\rho$ of 0.5 mm. In all of the simulations, the heat
source has a power of 30 watts when engaged.

The problem domain is a 4 by 6 cm grid, with a grid resolution of 0.2 mm; a
discretization time step of 2 ms is employed. We consider two data sets: one
in which the heat source is stationary (to obtain $a$ and $q$), and the other
in which the heat source is moving (to obtain $\tau$). The movement speed is
set at 20 mm/s. In all cases, the simulation time is 2 seconds, and the probe
is active for 1 second in the stationary case, and 1.5 seconds in the moving
case. The simulation was carried out using a custom simulation code based on a
finite difference scheme, as described in [2].

###  4.1 Thermal Diffusivity

Employing the method introduced in Sec. 3.1, we focus on the last 1 second
portion of the stationary simulation run. A cumulative estimate of the
diffusivity is shown in Fig. 2, expressed as the ratio $\hat{a}/a$, where the
estimate is based on the samples obtained until that instant. The estimated
thermal diffusivity $\hat{a}$ is $0.102$ mm2/s, which deviates 11% from the
ground truth. As one may observe, the estimate quickly converges to the ground
truth within 100 samples (0.2 s), showing little response to perturbations
over time.

Figure 2: Running estimate of the ratio between the ground truth thermal and
the estimated diffusivity, obtained during the unforced part of the stationary
thermal response. The dotted horizontal line corresponds to a perfect
estimate.

###  4.2 Heat Source

As mentioned in Sec. 3.2, the heat source must only be estimated after the
initial heat wave subsides, as illustrated by $\hat{\ddot{x}}_{\text{max}}$.
We find that the sampling index at which magnitudes below 10% of the maximum
acceleration are maintained is 250 (0.5 s).

We apply our heat source fitting method at index 250, after the thermal wave
has subsided, yielding the results shown in Fig. 3. Here, we have shown the
heat source fit with respect to the available data based on (9) (i.e., a
Fourier heat equation assumption); we have compared this data to the actual
case of the Maxwell–Cattaneo model, in which we have considered (7) using the
ground truth value of $\tau=0.1$ seconds. As can be seen in Fig. 3(b), we have
obtained an adequate fit with respect to the unknown Maxwell–Cattaneo-based
heat source. We may assess the accuracy of our model by considering the spread
$\rho$ in (11), which is obtained to be 0.56 mm, as opposed to the ground
truth 0.5 mm.

(a) Estimated Fourier heat source.

(b) Ground truth Maxwell–Cattaneo heat source.

Figure 3: Heat source model fit compared to (a) Fourier heat source assumption
used in fitting the model, (b) ground truth Maxwell–Cattaneo heat source
model, _after_ the initial heat wave has subsided. The solid line labeled
‘fit’ denotes the model, whereas the other marked lines are the heat fluxes in
each of the four cardinal directions. Almost no discrepancies can be observed
among the two models.

To demonstrate what the results would be if we considered the earlier stages
of the forced response when the thermal acceleration is still high (i.e., the
thermal wave has not subsided), we have developed the same estimates as in
Fig. 3 at index 50 (0.1 s), as shown in Fig. 4.

(a) Estimated Fourier heat source.

(b) Ground truth Maxwell–Cattaneo heat source.

Figure 4: Heat source model fit compared to (a) Fourier heat source assumption
used in fitting the model, (b) ground truth Maxwell–Cattaneo heat source
model, _before_ the initial heat wave has subsided. The solid line labeled
‘fit’ denotes the model, whereas the other marked lines are the heat fluxes in
each of the cardinal directions. Significant discrepancies can be observed
owing to the large thermal acceleration.

###  4.3 Thermal Relaxation Time

To obtain the thermal relaxation time $\tau$, we apply the method introduced
in Sec. 3.3 on the moving probe simulation. At each time index, we have given
the estimated thermal relaxation constant, as shown in Fig. 5. Clearly, many
of the values are erroneous, except for the initial few values, and more
accurately, the values that appear shortly after the heat source is disengaged
(as indicated by the dotted vertical line). Following the averaging criteria
given in Sec. 3.3, we obtain an estimate of $\hat{\tau}=0.09$ s, deviating 10%
from the ground truth $\tau=0.1$ s.

Figure 5: Thermal relaxation time estimate for the moving heat source
experiment. The dotted vertical line indicates the time instant at which the
heat source is disengaged. The sampling rate is 28 Hz.

##  5 Application

In this section, we apply our proposed method to the thermal response of _in
vivo_ porcine liver tissue to electrosurgical excitation. The trial was
performed at the University of Illinois at Chicago’s (UIC) Surgical Innovation
Training Laboratory (SITL), in accordance with NIH Vertebrate Animal Section
standards. The specimen was provided by the Biological Resources Laboratory of
the UIC. After exposing the liver tissue, a stationary point probe was
performed using a custom robotic platform based on a Franka Emika Panda
robotic arm, equipped with visual and thermal sensors (four Optris Xi 400
thermographers (Optris GmbH, Berlin, Germany)). The electrosurgical generator
used was a Valleylab Force FX (Valleylab Inc., CO, USA), applying the
monopolar pure cut mode at a power of 5 watts. For the stationary cut, a
Covidien Valleylab 2.8” needle electrode (Medtronic Inc., MN, USA) was
inserted at 5 mm depth.

To illustrate the various segments of the stationary probe experiment, the
maximum temperature recording during the trial over time is shown in Fig. 6.
During the interval 2.8–5.6 seconds, the electrosurgical unit was engaged.
Similar to the procedures in Sec. 4, the period during which the diffusivity
was estimated was chosen to be 6.2–12.4 seconds based on the second derivative
of temperature. Applying the same method as in Sec. 3.1, we obtain at thermal
diffusivity of $0.222$ mm2/s, which is about 53% higher than the non-perfused
ex vivo value reported in [9, Tab. 2.3, p. 16] [14]. This difference likely
stems from the presence of blood perfusion, which causes an increase in
effective diffusivity of 1.4–1.5 [9, Sec. 2.1.5.2] , consistent with our
result of 1.53 here.

Figure 6: Maximum temperature recorded during the stationary cut trial on live
porcine liver tissue. At 2.8 seconds, the power source is engaged, with
disengagement taking place at 5.6 seconds, as indicated by the dotted black
lines.

Having obtained the thermal diffusivity, we may now estimate the heat source
model. We delineate the period during which the heat source model is estimated
based on the maximum absolute temperature acceleration criterion (Sec. 3.2),
obtaining a time period of interest between 3.8–4.8 seconds. In this case, the
heat source estimate field is much more spurious when compared to Sec. 4 as
shown in Fig. 7, but most of the noise has been suppressed by the low-pass
properties of the noise-robust gradient operators. In the line corresponding
to the ‘up’ direction, however, one may notice an increase in heat flux around
2.75 mm; this anomaly can be attributed to thermal sensor noise, as well as
the presence of particulates (smoke), which inevitably introduce
discrepancies. Hence, we have based our heat source model on the average heat
flux across multiple directions.

Figure 7: Live porcine liver heat source model fit compared to Fourier heat
source assumption used in fitting the model after the initial heat wave has
subsided. The solid line labeled ‘fit’ denotes the model, whereas the other
marked lines are the heat fluxes in each of the four cardinal directions.
Almost no discrepancies can be observed among the two models.

To assess the quality of the heat source, we note that our estimated spread is
$\hat{\rho}=2.26$ mm, which is reasonable considering the diameter of the
probe (1 mm). The correctness of this result will further be discussed in Sec.
5.1.

We now consider obtaining the thermal relaxation time. Due to issues of image
alignment as a result of tissue motion caused by breathing and uncertainties
in the thermographer location, we have chosen to estimate the thermal
relaxation time based on the stationary probe response, rather than a moving
probe response. In an ideal situation, the latter method would yield more
conclusive results, since the change in heating may be much more substantial
when compared to what can be achieved by stationary probing, while minimizing
unnecessary tissue damage. While tissue motion rejection will be the subject
of future study, in this work we show that an adequate estimate may also be
obtained based on the stationary response. Applying the method of Sec. 3.3, we
find relaxation times as shown in Fig. 8. In the latter figure, we find that
the initial engagement of the probe results in two spikes in the thermal
relaxation constant estimate. This is a result of the probe temporarily losing
contact with the tissue after initial ablation, after which it reconnects with
the remaining tissue. This effects can also be observed in the maximum
temperature upon close inspection of Fig. 6, where there is a short period of
stagnation at the 3.5 second mark. After averaging the relaxation time as
described in Sec. 3.3, we find a relaxation time of $\hat{\tau}=0.16$ s. This
yields a second speed of sound of 1.2 mm/s by (4), which is close to the
results reported in [1, 2] of around 2 mm/s. The finite heat propagation time
can also be observed in the delay between the times where usable estimates of
$\tau$ are available in Fig. 8 (around $2\tau$ seconds).

Figure 8: Thermal relaxation time estimate for the stationary heat source
experiment on live porcine tissue. The dotted vertical lines indicate the
times at which the heat source is engaged and disengaged, respectively.

We may observe the effect of neglecting the thermal relaxation constant based
on the error between the Fourier-based and Maxwell–Cattaneo-based heat source
approximations with the thermal relaxation time being taken as 0.16 s, as
shown in Fig. 9. For ease of viewing, we have smoothened the raw heat flux
curves using a third-order Savitzky–Golay filter, with a window size of 31
samples. Clearly, _the Fourier heat approximation exhibits a large offset
close to the probe location_ , where an error of 10 W/m2 is observed in Fig.
9(b) at the point of application (0 mm distance); _such an error will yield
poorly agreeing simulation results across the entire domain_. In the next
subsection, we will show that our heat model produces an adequate response
when compared to the real-life thermal measurements.

(a) Estimated Fourier heat source.

(b) Expected Maxwell–Cattaneo heat source.

Figure 9: Heat source model fit compared to (a) Fourier heat source assumption
used in fitting the model, (b) ground truth Maxwell–Cattaneo heat source
model, _before_ the initial heat wave has subsided for the live porcine liver
response. The observed heat flux curves have been smoothened using a
Savitzky–Golay filter for ease of viewing. The solid line labeled ‘fit’
denotes the model, whereas the other marked lines are the heat fluxes in each
of the cardinal directions. Significant discrepancies can be observed owing to
the large thermal acceleration in (a), especially close to the point of
application (0 mm distance).

###  5.1 Comparison with Simulation

We rerun the simulation using the parameter estimates obtained in the
foregoing; we have chosen the initial temperature to be 28 ∘C, which is the
average temperature prior to engaging the electrosurgical unit. To illustrate
the agreement between the observed real-life thermal response, and the
simulated response, the maximum temperature time series plotted in Fig. 10
provides a clear overview. We note that given the parameters obtained by
application of our method, the simulated maximum temperature discrepancy is
under 5 ∘C (6%) with respect to real-life data across the entire simulation
run. To provide further evidence of the parameter fidelity, we compare the
temperature fields obtained at 2.5 seconds, when the probe is disengaged. In
Fig. 11 the observed and simulated temperature fields are shown. Most of the
error can be attributed to misalignment of the camera and the liver tissue’s
curvature, but despite these effects, the maximum temperature error is less
than 12 ∘C (14%). We have previously addressed mismatches between simulation
and sensor measurements in detail, proving that an observer structure with
output feedback can eliminate these discrepancies [8].

Figure 10: Maximum temperature measured over time for both the real-life
observed thermal response, and the simulated response of porcine liver tissue
subjected to a stationary 5 W electrosurgical heat input. The dotted vertical
line indicates the point of probe disengagement. Clear congruence between the
results can be observed, with small overall errors.

(a) Meas. temperature field.

(b) Sim. temperature field.

Figure 11: Comparison between the temperature field at the point of probe
disengagement, for both (a) the measured real-life data, as well as (b) the
simulated data.

##  6 Conclusion

We have presented a novel approach to real-time parameter estimation for a
hyperbolic heat transfer model using thermographic data, as applied to _in
vivo_ porcine liver tissue subjected to electrosurgical action. Unlike past
methods, which require the solution of complex PDE systems, our approach,
based on our previous work on _attention-based noise robust averaging_ (ANRA)
[8], operates directly on thermographic data and requires only the PDE
structure to be given. We have applied our approach to a theoretical test case
involving a Maxwell–Cattaneo model with parameters based on porcine liver
tissue, as well as real-life data from an _in vivo_ porcine liver subject to
electrosurgical action. In all cases, our method produces parameter estimates
close to the ground truth or parameters reported in the literature. Using the
results obtained by our method, we have also presented a comparison between
live and dead tissue thermal behavior. In addition, a simulation was run using
these new parameter estimates, achieving clear correspondence with the real-
life porcine liver thermal response. With this validation, we have shown that
_our method is capable of in situ real-time minimally invasive high-fidelity
thermophysical parameter estimation_.

Some of the limitations of this work, planned to be addressed in our future
research, arise due to the challenging live tissue geometry and breathing
motion, and electrosurgical smoke and particles causing temperature reading
distortion. The future work may extend this approach to other thermodynamic
models, such as the reaction-diffusion equation. In addition, we plan to
develop methods to discern tissue types in real-time based on their thermal
response. Related to this is planned work on comparing tissue properties in
_ex vivo_ and _in vivo_ tissue, with the goal of deriving robust relationships
between the two. We also plan to further investigate thermal shock waves and
their significance to energy-based surgery and the tissue thermal response.

## Animal Use Approval Details

The animal use protocols (ACC NO: 19-151 and 22-097) were reviewed in
accordance with the Animal Care Policies and Procedures of the University of
Illinois at Chicago. Original protocol approval: 9/24/2019 (3 year approval
with annual continuation required). Current approval period: 8/9/2022 to
7/19/2025.

## References

  * [1] A. Madhukar _et al._ , “Heat conduction in porcine muscle and blood: Experiments and time-fractional telegraph equation model,” _Journal of The Royal Society Interface_ , vol. 16, no. 160, p. 20190726, Nov. 2019.
  * [2] J. Ran _et al._ , “Telegraph equation shock waves in biological tissues,” _International Journal for Multiscale Computational Engineering_ , 2022.
  * [3] J. Ran, M. Ostoja-Starzewski, and Y. Povstenko, “Mach fronts in random media with fractal and Hurst effects,” _Fractal and Fractional_ , vol. 5, no. 4, p. 229, 2021.
  * [4] H. El-Kebir and J. Bentsman, “PDE-based modeling and non-collocated feedback control of electrosurgical-probe/tissue interaction,” in _2021 American Control Conference_. New Orleans, Louisiana, USA: IEEE, 2021, pp. 4045–4050.
  * [5] H. El-Kebir _et al._ , “Online hypermodel-based path planning for feedback control of tissue denaturation in electrosurgical cutting,” in _11th IFAC Symposium on Biological and Medical Systems_ , vol. 54. Ghent, Belgium: IFAC, 2021, pp. 448–453.
  * [6] C.-H. Yang, W. Li, and R. K. Chen, “Characterization and modeling of tissue thermal conductivity during an electrosurgical joining process,” _IEEE Transactions on Biomedical Engineering_ , vol. 65, no. 2, pp. 365–370, 2018.
  * [7] ——, “Determination of tissue thermal conductivity as a function of thermal dose and its application in finite element modeling of electrosurgical vessel sealing,” _IEEE Transactions on Biomedical Engineering_ , vol. 67, no. 10, pp. 2862–2869, Oct. 2020.
  * [8] H. El-Kebir _et al._ , “Infinite-dimensional adaptive boundary observer for inner-domain temperature estimation of 3D electrosurgical processes using surface thermography sensing,” in _Proceedings of the 2022 IEEE Control and Decision Conference_. Cancún, Mexico: IEEE, 2022.
  * [9] F. A. Duck, _Physical Properties of Tissues_. San Diego, CA, USA: Academic Press Inc., 1990.
  * [10] S. Singh and R. Melnik, “Thermal ablation of biological tissues in disease treatment: A review of computational models and future directions,” _Electromagnetic Biology and Medicine_ , vol. 39, no. 2, pp. 49–88, Apr. 2020\.
  * [11] L. C. Evans, _Partial Differential Equations_ , 2nd ed., ser. Graduate Studies in Mathematics. Providence, Rhode Island: American Mathematical Society, 2010, no. 19.
  * [12] P. Holoborodko, “Smooth Noise Robust Differentiators,” http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/, 2008\.
  * [13] ——, “Noise Robust Gradient Operators,” http://www.holoborodko.com/pavel/image-processing/edge-detection/, 2009.
  * [14] J. W. Valvano, J. T. Allen, and H. F. Bowman, “The simultaneous measurement of thermal conductivity, thermal diffusivity, and perfusion in small volumes of tissue,” _Journal of Biomechanical Engineering_ , vol. 106, no. 3, pp. 192–197, 1984.


&+\left(\frac{d}{dt}(P)^{H}P\right)RDR^{H}+RDR^{H}\left(P^{H}\frac{d}{dt}(P)\right).\end{split}$
(14)

Due to the fact that $P^{H}P=I_{r}$, we have
$\frac{d}{dt}(P)^{H}P+P^{H}\frac{d}{dt}(P)=0$. Denote
$M=P^{H}\frac{d}{dt}(P)$, then $M$ is an antisymmetric matrix such that
$M+M^{H}=0$. Therefore, we have

$\displaystyle\frac{d}{dt}\left(RDR^{H}\right)=M^{H}RDR^{H}+RDR^{H}M+\Sigma_{z}^{-1}RDQ^{H}QDR^{H}+RDQ^{H}QDR^{H}\Sigma_{z}^{-1}.$

On the other hand, we have

$\displaystyle\frac{d}{dt}(\Sigma_{z})$
$\displaystyle=\frac{d}{dt}(P^{H}S_{z}P)$
$\displaystyle=\frac{d}{dt}(P)^{H}S_{z}P+P^{H}\frac{d}{dt}(S_{z})P+P^{H}S_{z}\frac{d}{dt}(P)$
$\displaystyle=P^{H}\left(R_{*}DR_{*}^{H}-S_{z}\right)P+\left(\frac{d}{dt}(P)^{H}P\right)\Sigma_{z}+\Sigma_{z}\left(P^{H}\frac{d}{dt}(P)\right)$
$\displaystyle=RDR^{H}-\Sigma_{z}+\left(\frac{d}{dt}(P_{z})^{H}P_{z}\right)\Sigma_{z}+\Sigma_{z}\left(P_{z}^{H}\frac{d}{dt}(P_{z})\right).$

Thus, we arrive at

$\begin{split}\frac{d}{dt}(\Sigma_{z}(s,s))&=(RDR^{H})(s,s)-\Sigma_{z}(s,s).\end{split}$
(15)

Since $\frac{d}{dt}(\text{offdiag}(\Sigma_{z}))=0$, we have

$\displaystyle\text{offdiag}(RDR^{H})=-M^{H}\Sigma_{z}-\Sigma_{z}M=M\Sigma_{z}+\Sigma_{z}M^{H}.$
(16)

By substitution, we get

$\displaystyle\frac{d}{dt}(RDR^{H})(j,j)$
$\displaystyle=e_{j}^{H}(2M\Sigma_{z}M-\Sigma_{z}MM-
MM\Sigma_{z})e_{j}+\frac{2}{\sigma_{j}(Z)}\|QDR^{H}e_{j}\|_{2}^{2}$
$\displaystyle=2\sum_{k=1}^{r}(\sigma_{j}(Z)-\sigma_{k}(Z))M(k,j)^{2}+\frac{2}{\sigma_{j}(Z)}\|QDR^{H}e_{j}\|_{2}^{2}$
(17) $\displaystyle=2\sum_{k\neq
j}\frac{1}{\sigma_{j}-\sigma_{k}}(RDR^{H})(k,j)^{2}+\frac{2}{\sigma_{j}(Z)}\|QDR^{H}e_{j}\|_{2}^{2}.$

Note that

$\displaystyle\sum_{j=1}^{r}\sum_{k=1}^{r}(\sigma_{j}(Z)-\sigma_{k}(Z))M(k,j)^{2}$
$\displaystyle=-\sum_{k=1}^{r}\sum_{j=1}^{r}(\sigma_{k}(Z)-\sigma_{j}(Z))M(k,j)^{2}$
$\displaystyle=-\sum_{k=1}^{r}\sum_{j=1}^{r}(\sigma_{j}(Z)-\sigma_{k}(Z))M(k,j)^{2}.$

Thus, we obtain

$\displaystyle\frac{d}{dt}(\sum_{j=1}^{r}\|D^{\frac{1}{2}}R^{H}e_{j}\|^{2})$
$\displaystyle=\sum_{j=1}^{r}\frac{2}{\sigma_{j}(Z)}\|QDR^{H}e_{j}\|^{2}.$
(18)

  2. 2)

Denote $R=T_{2}\Sigma_{RR}T_{1}^{H}$. Simple calculation gives
$R^{H}R=T_{1}\Sigma_{RR}T_{1}^{H}$ and $Q^{H}Q=T_{1}\Sigma_{QQ}T_{1}^{H}$ with
$\Sigma_{RR}+\Sigma_{QQ}=I_{r}$. We have

$\displaystyle\sum_{j=1}^{r}(RDR^{H})(j,j)$
$\displaystyle=\text{trace}(RDR^{H})=\text{trace}(DR^{H}R)$
$\displaystyle=\text{trace}(T_{1}^{H}DT_{1}\Sigma_{RR})$
$\displaystyle=\sum_{j=1}^{r}t_{j}r_{j}.$

Here, for each $j\in[r]$, $t_{j}=\|D^{\frac{1}{2}}T_{1}e_{j}\|^{2}=\Omega(1)$
and $\Sigma_{RR}(j,j)=r_{j}$. Now, for a bounded sequence $\\{Z_{t}\\}$ we
have

$\displaystyle\sum_{j=1}^{r}\frac{2}{\sigma_{j}}\|QDR^{H}e_{j}\|^{2}$
$\displaystyle\gtrsim\text{trace}(RDQ^{H}QDR^{H})$
$\displaystyle=\text{trace}(T_{1}^{H}DT_{1}\Sigma_{QQ}T_{1}^{H}DT_{1}\Sigma_{RR}).$

Denote $N=T_{1}^{H}DT_{1}$. We obtain the following estimate:

$\displaystyle\sum_{j=1}^{r}\frac{2}{\sigma_{j}}\|QDR^{H}e_{j}\|^{2}$
$\displaystyle\gtrsim\sum_{j=1}^{r}r_{j}(1-r_{j})N(j,j)^{2}$
$\displaystyle\gtrsim\sum_{j=1}^{r}r_{j}(1-r_{j}).$

This gives

$\displaystyle\frac{d}{dt}(\sum_{j=1}^{r}r_{j})\gtrsim\sum_{j=1}^{r}r_{j}(1-r_{j}).$

Further, if $\sigma_{j}=\Omega(1)$, for all $j\in[r]$, by equation (18), we
have

$\displaystyle\frac{d}{dt}(\sum_{j=1}^{r}\|D^{\frac{1}{2}}R^{H}e_{j}\|^{2})$
$\displaystyle=\Omega\left(\sum_{j=1}^{r}\|QDR^{H}e_{j}\|^{2}\right)$
$\displaystyle=\Omega\left(\text{trace}\left(T_{1}^{H}DT_{1}\Sigma_{QQ}T_{1}^{H}DT_{1}\Sigma_{RR}\right)\right)$
$\displaystyle=\Omega\left(\sum_{j=1}^{r}r_{j}(1-r_{j})\right).$

Thus we have
$\frac{d}{dt}\sum_{j=1}^{r}r_{j}=\Omega\left(\sum_{j=1}^{r}(1-r_{j})r_{j}\right)$.

∎

###  B.9 Proof of Lemma 3.18

###### Proof.

For simplicity, the following proof is based on the symmetric case where
$Z=U_{z}\Sigma U_{z}^{H}$ and $X=UDU^{H}$.

  1. 1)

Let $\widetilde{U}_{z}\in\mathbb{F}^{n\times(n-r)}$ be the orthogonal
complement of $U_{z}$, and

$\displaystyle
X=(U_{z},\widetilde{U}_{z})\widetilde{X}(U_{z},\widetilde{U}_{z})^{H}=(U_{z},\widetilde{U}_{z})\begin{pmatrix}X_{11}&X_{12}\\\
X_{21}&X_{22}\end{pmatrix}\begin{pmatrix}U_{z}^{H}\\\
\widetilde{U}_{z}^{H}\end{pmatrix}.$

Then, we get

$\displaystyle\|P_{T_{z}}(Z-X)\|_{F}$
$\displaystyle=\left\|\begin{pmatrix}X_{11}-\Sigma_{z}&X_{12}\\\
X_{21}&0\end{pmatrix}\right\|_{F}=\|X_{11}-\Sigma_{z}\|_{F}^{2}+\|X_{12}\|_{F}^{2}+\|X_{21}\|_{F}^{2},$
$\displaystyle\|Z-X\|_{F}^{2}$
$\displaystyle=\|P_{T_{z}}(Z-X)\|_{F}^{2}+\|X_{22}\|_{F}^{2}.$

We also assume that $U=U_{z}R+\widetilde{U}_{z}S$, where
$R\in\mathbb{F}^{r\times r}$ and $S\in\mathbb{F}^{(n-r)\times r}$. Then, we
obtain

$\displaystyle X_{11}=RDR^{H},\quad X_{12}=RDS^{H},\quad X_{21}=SDR^{H},\quad
X_{22}=SDS^{H}.$

The goal is to find lower and upper bounds for

$\displaystyle s=\frac{\|P_{T_{z}}(Z-X)\|_{F}^{2}}{\|Z-X\|_{F}^{2}}.$

It is obvious that $s\leq 1$. To identify the lower bound we consider

$\displaystyle\phi=\frac{\|Z-X\|_{F}^{2}}{\|X_{22}\|_{F}^{2}},\quad
s=\frac{\phi-1}{\phi},\quad\phi\geq 1.$

For fixed $\Sigma$ and $D$, $\phi$ can be seen as a function of $R$. We
express $\phi(R)$ in terms of $\phi_{1}$ and $\phi_{2}$ defined below:

$\displaystyle\phi(R)=\frac{\phi_{1}(R)}{\phi_{2}(R)},$ $\displaystyle\quad
R\in\mathcal{D}\subset\mathbb{F}^{r\times r},\,\mathcal{D}=\\{R:0\preccurlyeq
R^{H}R\preccurlyeq I_{r}\\},$ $\displaystyle\text{where }\phi_{1}(R)$
$\displaystyle=\|X\|_{F}^{2}+\|\Sigma\|_{F}^{2}-2\langle\Sigma,X_{11}\rangle$
$\displaystyle=\|D\|_{F}^{2}+\|\Sigma\|_{F}^{2}-2\text{tr}(\Sigma RDR^{H}),$
$\displaystyle\phi_{2}(R)$
$\displaystyle=\text{tr}\left((SDS^{H})^{2}\right)=\text{tr}\left((DS^{H}S)^{2}\right)$
$\displaystyle=\text{tr}\left((D(I-R^{H}R))^{2}\right).$

To minimize $t$ for given $D$ and $\Sigma$, the first-order condition is
obtained by taking derivative of $t$ over $R$:

$\displaystyle\frac{\partial t}{\partial R}$
$\displaystyle=\frac{1}{\phi_{2}^{2}}(\frac{\partial\phi_{1}}{\partial
R}\phi_{2}-\frac{\partial\phi_{2}}{\partial R}\phi_{1})$
$\displaystyle=\frac{1}{\phi_{2}^{2}}(-4\phi_{2}\Sigma
RD+4\phi_{1}RD(I-R^{H}R)D)$
$\displaystyle=\frac{4}{\phi_{2}^{2}}(-\phi_{2}\Sigma RD+\phi_{1}RDS^{H}SD).$

Imposing $\partial t/\partial R=0$ gives

$\displaystyle\Sigma RD=\phi\cdot RDS^{H}SD.$ (19)

We now claim that this first-order condition _cannot_ be satisfied in the
interior of the domain $\mathcal{D}$. To see this, note that the above
equation (19) gives

$\displaystyle\Sigma RDR^{H}=\phi\cdot RDS^{H}SDR^{H},\quad\text{ i.e.
}\quad\langle\Sigma,X_{11}\rangle=\phi\|X_{12}\|_{F}^{2}=\phi\|X_{21}\|_{F}^{2}.$

Thus, we have

$\displaystyle\phi$
$\displaystyle=\frac{\|Z-X\|_{F}^{2}}{\|X_{22}\|_{F}^{2}}=\frac{\|\Sigma\|_{F}^{2}+\|X_{11}\|_{F}^{2}-2\langle\Sigma,X_{11}\rangle+\|X_{12}\|_{F}^{2}+\|X_{21}\|_{F}^{2}+\|X_{22}\|_{F}^{2}}{\|X_{22}\|_{F}^{2}}$
$\displaystyle=\frac{\|\Sigma\|_{F}^{2}+\|X_{11}\|_{F}^{2}-(2-2/\phi)\langle\Sigma,X_{11}\rangle}{\|X_{22}\|_{F}^{2}}+1$
$\displaystyle=\frac{\|X_{11}-(1-1/\phi)\Sigma\|_{F}^{2}}{\|X_{22}\|_{F}^{2}}+(\frac{2}{\phi}-\frac{1}{\phi^{2}})\frac{\|\Sigma\|_{F}^{2}}{\|X_{22}\|_{F}^{2}}+1$
$\displaystyle\geq(\frac{2}{\phi}-\frac{1}{\phi^{2}})\frac{\|\Sigma\|_{F}^{2}}{\|X_{22}\|_{F}^{2}}+1.$

Note that when $R\in int(\mathcal{D})$, $R$ is full-rank. Since $D,\Sigma$ are
also full-rank, (19) gives $\phi^{2}=\|\Sigma\|_{F}^{2}/\|X_{22}\|_{F}^{2}$ by
some simple matrix manipulation:

$\displaystyle\Sigma RD=\phi\cdot RDS^{H}SD$ $\displaystyle\Rightarrow\quad$
$\displaystyle R^{-1}\Sigma R=\phi\cdot DS^{H}S$
$\displaystyle\Rightarrow\quad$ $\displaystyle\text{tr}\left((R^{-1}\Sigma
R)^{2}\right)=\phi^{2}\text{tr}\left((DS^{H}S)^{2}\right)$
$\displaystyle\Rightarrow\quad$
$\displaystyle\text{tr}\left(\Sigma^{2}\right)=\phi^{2}\text{tr}\left((SDS^{H})^{2}\right)$
$\displaystyle\Rightarrow\quad$
$\displaystyle\|\Sigma\|_{F}^{2}=\phi^{2}\|X_{22}\|_{F}^{2}.$

Hence, we obtain

$\displaystyle\phi\geq(\frac{2}{\phi}-\frac{1}{\phi^{2}})\phi^{2}+1=2\phi-1+1=2\phi,$

which contradicts the fact that $\phi\geq 1$.

Therefore, either the extreme values of $t$ are only achieved on
$\partial\mathcal{D}$, or the full-rankness of $D$ is violated. In either
case, the RHS of (19) is rank-deficient. So the LHS of (19) (and thus $R$ and
$X_{11}$) is also rank-deficient. Therefore, we conclude that

$\displaystyle\phi-1\geq\frac{\|X_{11}-\Sigma\|_{F}^{2}}{\|X_{22}\|_{F}^{2}}\geq\frac{\|X_{11}-\Sigma\|_{F}^{2}}{\|X\|_{F}^{2}}\geq\frac{\sigma_{r}^{2}}{\|X\|_{F}^{2}},$

where $\sigma_{r}$ is the smallest nonzero eigenvalue of $Z$. Finally, we
obtain

$\displaystyle s\geq\frac{\sigma_{r}^{2}}{\sigma_{r}^{2}+\|X\|_{F}^{2}}.$

Now we have $\|P_{T_{z}}(Z-X)\|_{F}^{2}\geq s\|Z-X\|_{F}^{2}$ for $s>0$ as
long as $Z$ is of rank $r$ and the smallest nonzero singular value of $Z$ is
bounded away from 0.

  2. 2)

By Lemma C.2, when $\|Z-X\|_{F}<\frac{d_{r}}{2}$, we have
$\sigma_{r}>\frac{d_{r}}{2}$. Further, by (1) we have
$\frac{\|P_{T_{z}}(Z-X\|_{F}}{\|Z-X\|_{F}}\gtrsim d_{r}$. By Lemma 3.13, we
have the locally linear convergence towards $X$ with
$c=-\log\left(1-\Omega(\alpha d_{r}^{2})\right)$. To prove the result for
single-step convergence, we consider

$\displaystyle f_{k+1}$ $\displaystyle\leq
f_{k}-C_{2}\|P_{T_{z}}(Z_{k}-X)\|_{F}\|Z_{k+1}-Z_{k}\|_{F}$ $\displaystyle\leq
f_{k}-C_{2}(\alpha+o(\alpha))\cdot\|P_{T_{z}}(Z_{k}-X)\|_{F}^{2}$
$\displaystyle\leq(1-\Omega(\alpha)C_{1}^{2}C_{2})f_{k}.$

Here, $f_{k}:=\frac{1}{2}\|Z_{k}-X\|_{F}^{2}$. The first inequality is from
condition (D). The second inequality is by the first-order retraction, and the
third inequality is from condition (L). Therefore, we have

$\displaystyle\|Z_{k+1}-X\|_{F}\leq\sqrt{1-\Omega(\alpha)C_{1}^{2}C_{2}}\|Z_{k}-X\|_{F}.$

That is, $\|Z_{k+1}-X\|_{F}\leq e^{-c}\|Z_{k}-X\|_{F}$, with
$c=-\log(1-\Omega\left(\alpha d_{r}^{2})\right)$.

∎

###  B.10 Proof of Theorem 4.1

###### Proof.

The proof of Theorem 4.1 can be reduced to deducing Conditions (L) and (D)
from Conditions (1) and (2) in the assumptions. Note that in the setting of
$T:\mathcal{M}\rightarrow\mathbb{R}^{m}$,
$T(Z)=\frac{1}{\sqrt{m}}\begin{pmatrix}\langle A_{1},Z\rangle\\\ \langle
A_{2},Z\rangle\\\ ...\\\ \langle A_{m},Z\rangle\end{pmatrix}$,
$f(Z)=\frac{1}{2}\|T(Z)-T(X)\|_{2}^{2}$, Conditions (L) and (D) can be
formulated as follows.

  1. 1)

Condition (L): Łojasiewicz gradient inequality

$\left(\frac{1}{2m}\sum_{j}\langle A_{j},Z-X\rangle^{2}\right)^{1-\omega}\leq
C_{l}\|P_{T_{z}}(\frac{1}{m}\sum_{j}\langle A_{j},Z-X\rangle A_{j})\|_{F}$

holds with $\omega=\frac{1}{2}$.

  2. 2)

Condition (D):

$f_{k}-f_{k+1}\geq C_{d}\|P_{T_{Z_{k}}}(\frac{1}{m}\sum_{j}\langle
A_{j},Z_{k}-X\rangle A_{j})\|_{F}\|Z_{k+1}-Z_{k}\|_{F}.$

The proof now goes as follows.

  1. 1)

To prove condition (L): By (1) and (2) we have

$\displaystyle\|P_{T_{z}}(\nabla f(Z))\|_{F}$
$\displaystyle\geq\frac{1}{C_{3}}\|Z-X\|_{F}$
$\displaystyle\geq\frac{1}{C_{2}C_{3}}\|T(Z)-T(X)\|$
$\displaystyle=\frac{\sqrt{2}}{C_{2}C_{3}}\left|f(Z)-f(X)\right|^{\frac{1}{2}},$

which implies that (L) holds with $\omega=\frac{1}{2}$ and
$C_{l}=\frac{\sqrt{2}C_{2}C_{3}}{2}>0$ is an abosolute constant.

  2. 2)

To prove condition (D): we first consider

$\displaystyle f_{k}-f_{k+1}$
$\displaystyle=\frac{1}{2}\|T(Z_{k})-T(X)\|_{2}^{2}-\frac{1}{2}\|T(Z_{k+1})-T(X)\|_{2}^{2}$
$\displaystyle=\frac{1}{2}\langle T(Z_{k}+Z_{k+1}-2X),T(Z_{k}-Z_{k+1})\rangle$
$\displaystyle=\langle
T^{*}T(X-Z_{k}),Z_{k+1}-Z_{k}\rangle-\frac{1}{2}\|T(Z_{k+1}-Z_{k})\|_{F}^{2}.$

Assume that $Z_{k+1}=Z_{k}+\alpha\widetilde{\xi}_{k}$, and $-\alpha
P_{T_{Z_{k}}}(\nabla f(Z_{k}))=\alpha\xi_{k}$. By first-order retraction
property and condition (1), we get

$\displaystyle f_{k}-f_{k+1}$ $\displaystyle\geq\langle-\nabla
f(Z_{k}),\alpha\widetilde{\xi}_{k}\rangle-\frac{C_{2}^{2}\alpha^{2}}{2}\|\widetilde{\xi}_{k}\|^{2}$
$\displaystyle=\langle-\nabla
f(Z_{k}),\alpha\xi_{k}\rangle+o(\alpha\|\xi_{k}\|_{F}^{2})$
$\displaystyle=\langle-\nabla f(Z_{k}),-P_{T_{Z_{k}}}(\nabla
f(Z_{k}))\rangle+o(\alpha\|\xi_{k}\|_{F}^{2})$
$\displaystyle=\alpha\|P_{T_{Z_{k}}}(\nabla
f(Z_{k}))\|_{F}^{2}+o(\alpha\|\xi_{k}\|_{F}^{2}).$

On the other hand, we also obtain

$\displaystyle C_{d}\|P_{T_{Z_{k}}}(\nabla
f(Z_{k}))\|_{F}\|Z_{k+1}-Z_{k}\|_{F}$
$\displaystyle=C_{d}\|P_{T_{Z_{k}}}(\nabla
f(Z_{k}))\|_{F}\|\alpha\widetilde{\xi}_{k}\|_{F}$
$\displaystyle=C_{d}\alpha\|P_{T_{Z_{k}}}(\nabla
f(Z_{k}))\|_{F}^{2}+o(\alpha\|\xi_{k}\|_{F}^{2}).$

By choosing $C_{d}>0$ small enough, we have

$f_{k}-f_{k+1}\geq C_{d}\|P_{T_{Z_{k}}}(\frac{1}{m}\sum_{j}\langle
A_{j},Z_{k}-X\rangle A_{j})\|_{F}\|Z_{k+1}-Z_{k}\|_{F},$

i.e. Condition (D) holds.

From Theorem 3.4, we conclude that projected gradient descent for the least
squares loss function $f(Z)=\frac{1}{2}\|T(Z)-T(X)\|_{2}^{2}$ converges to its
global minimum linearly. ∎

##  Appendix C Auxiliary Lemmas

######  Lemma C.1 (The $\text{sin}(\theta)$ Theorem, [19]).

Let $A$ be a Hermitian operator. Assume that

$\displaystyle
A=\begin{pmatrix}E_{0}&E_{1}\end{pmatrix}\begin{pmatrix}A_{0}&0\\\
0&A_{1}\end{pmatrix}\begin{pmatrix}E_{0}^{H}\\\ E_{1}^{H}\end{pmatrix}$

is an invariant subspace decomposition (i.e. generalized eigenvalue
decomposition) of $A$. Let

$\displaystyle B=A+\Delta,\quad
B=\begin{pmatrix}F_{0}&F_{1}\end{pmatrix}\begin{pmatrix}B_{0}&0\\\
0&B_{1}\end{pmatrix}\begin{pmatrix}F_{0}^{H}\\\ F_{1}^{H}\end{pmatrix}.$

Let $\Theta_{0}$ be the angle matrix between subspaces $E_{0}$ and $F_{0}$.
Define the residual as

$\displaystyle R=BE_{0}-E_{0}A_{0}.$

If there is an interval $[\beta,\alpha]$ and $\delta>0$, such that the
spectrum of $A_{0}$ lies entirely in $[\beta,\alpha]$, while that of $B_{1}$
lies entirely in $(-\infty,\beta-\delta]\cup[\alpha+\delta,+\infty)$, then for
every unitary-invariant norm, we have

$\displaystyle\delta\|\text{sin}\Theta_{0}\|\leq\|R\|.$

In particular, this holds true for the matrix 2-norm and the Frobenius norm.

######  Lemma C.2 (Hoffman-Wielandt Theorem).

  1. 1)

Assume $Z,Z^{\prime}\in\mathbb{R}^{n\times n}$ or $\mathbb{C}^{n\times n}$ are
normal matrices, and their corresponding ordered spectra are
$\\{\lambda_{j}\\}$ and $\\{\tilde{\lambda_{j}}\\}$. Then, we have

$\displaystyle\sqrt{\sum_{j=1}^{n}|\tilde{\lambda_{j}}-\lambda_{j}|^{2}}\leq\|Z^{\prime}-Z\|_{F}.$

  2. 2)

Assume $Z$, $Z^{\prime}\in\mathbb{R}^{m\times n}$ or $\mathbb{C}^{m\times n}$,
and denote their singular values in descending order as $\\{\sigma_{j}\\}$ and
$\\{\tilde{\sigma_{j}}\\}$. Then, we have

$\displaystyle\sqrt{\sum_{j=1}^{n}|\tilde{\sigma_{j}}-\sigma_{j}|^{2}}\leq\|Z^{\prime}-Z\|_{F}.$

######  Lemma C.3 (Dynmical low-rank approximation, [30]).

Assume $Z=USV^{H}$ with
$M:=P_{T_{Z}}(\frac{d}{dt}g(Z))=\nabla_{\mathcal{M}_{r}}g(Z)$, then gradient
flow satisfies:

$\displaystyle\frac{d}{dt}S=U^{H}MV,$
$\displaystyle\frac{d}{dt}U=P_{U}^{\perp}MVS^{-1},$
$\displaystyle\frac{d}{dt}V=P_{V}^{\perp}M^{H}US^{-\top}.$

Here, $P_{U}^{\perp}=I-UU^{H}$ and $P_{V}^{\perp}=I-VV^{H}$.