=80pt

# The Role of Intent-Based Networking in ICT Supply Chains

Mounir Bensalem, Jasenka Dizdarević, Francisco Carpio, and Admela Jukan
Technische Universität Braunschweig, Germany  {mounir.bensalem, j.dizdarevic,
f.carpio<EMAIL_ADDRESS>
###### Abstract

The evolution towards Industry 4.0 is driving the need for innovative
solutions in the area of network management, considering the complex, dynamic
and heterogeneous nature of ICT supply chains. To this end, Intent-Based
networking (IBN) which is already proven to evolve how network management is
driven today, can be implemented as a solution to facilitate the management of
large ICT supply chains. In this paper, we first present a comparison of the
main architectural components of typical IBN systems and, then, we study the
key engineering requirements when integrating IBN with ICT supply chain
network systems while considering AI methods. We also propose a general
architecture design that enables intent translation of ICT supply chain
specifications into lower level policies, to finally show an example of how
the access control is performed in a modeled ICT supply chain system.

##  I Introduction

In the recent years intent-based networking (IBN) has been the emerging
networking paradigm, levering the integration of AI/ML techniques with network
orchestration to improve upon the network automation, complexity and control
functionalities. Compared to traditional network management solutions it has
shown great potential in numerous aspects, reducing the time for network
configuration deployment while improving on its scalability[1], increasing
both the efficiency and reliability, as well as reducing costs in particularly
complex industrial networks [2] and integrating security through its intent
functions [3]. The key behind IBN’s increasing abilities to solving complex
network tasks lies is in their high-level policies specification capabilities
that allow network administrator to define goals and demands on how the
network should behave.

One of the particularly significant IBN potential application domains are ICT
supply chains, since the fast-growing technological evolution and the digital
transformation behind Industry 4.0 and ICT supply chains have created a gap
between traditional networking solutions and these emerging complex networks.
As now ICT supply chains can span across multiple vendor infrastructures and
devices, assuming different configuration requirements and capabilities with
continuous addition of new functionalities, the interest in novel network
management and orchestration approaches has been growing rapidly. The
utilization of IBN in this domain could facilitate automation, control and
orchestration processes, while improving on the network robustness and
operation dynamicity.

This paper introduces an IBN approach to automating the configuration
management and control of ICT supply chain networks through high-level
intents. The objective is to first analyse the key engineering requirements in
integrating IBN approach with these complex networks, considering the
questions of centralized vs decentralized orchestration, influence of
different AI techniques, and manner of specifying high level intents to ICT
supply chain policies. To this end we propose an architectural design of
intent based translator for ICT supply chains, as well as an example in the
context of supply chain assets access control.

The rest of the paper is organized as follows. Section 2 describes the
background and motivation behind using intent-based approaches. Section 3
gives an overview of engineering requirements, the general architectural
design for utilization of IBN in ICT supply chains, and a use case scenario
for assets access control. Finally, Section 4 concludes the paper.

##  II Background and motivation

In this section we analyse the common architectural building blocks of any
intent-based system, according to the already existing solutions. While the
design of intent-based translators can follow guidelines depending on the
scenario to be applied, this analysis was used as the initial step in
proposing an IBN architectural design of the ICT supply chain system. Table I
shows a comparative analysis of the intent abstraction layer design for ICT
systems. Several proposals from the literature are investigated to understand
the main software components responsible on interfacing with the user,
storage, management, policies, networking, monitoring, and enabling
intelligence. It is important to mention that few proposals highlighted the
role of AI and ML in translating user intents [4] and improving policies using
insight from telemetry [5, 6, 7].

Work |  Interface |  Management |  Storage |  Policies |  Network |  AI/ML |  Monitoring |  Application
---|---|---|---|---|---|---|---|---
[8] |  GUI |  intent parser & compiler |  service library |  policy builder |  switch controller |  no |  no |  switches programmability
[9] |  extended NBI |  intent manager & compiler |  vocabulary store |  conflict checking |  SDN controller |  no |  yes |  —
[10] |  GUI, Alexa |  intent engine & compiler |  MariaDB |  |  floodlight controller |  — |  no |  network visualization
[11] |  UI |  intent manager & compiler |  service registry |  service integration element |  SDN controller |  no |  no |  SDN applications
[12] |  GUI, Alexa |  intent engine |  database |  |  Ryu controller |  no |  — |  self-Healing SDN
[7] |  NBI |  intent manager |  policy store |  policy configurator |  resource controller |  LSTM, GAN |  yes |  orchestration of network slicing
[5] |  NBI |  policy parser |  policy store |  conflicts management |  SDN controller |  heuristics |  yes |  SDN-packet processing
[13] |  REST API |  intent engine |  no |  PCE |  Ryu Controller |  no |  no |  mobile backhauling
[6] |  GUI |  intent manager |  policy store |  policy configurators |  FlexRAN, open MANO |  DL, GAN |  yes |  network slicing
[4] |  slack API |  intent manager & compiler |  ontology and RDF |  conflict checking |  — |  LSTM |  — |  high-speed networks
[14] |  NETCONF |  extractor & compiler |  — |  convertor |  — |  no |  no |  IoT device configuration
our work |  GUI |  intent manager & compiler |  policy store, intent store, telemetry |  policy configurator |  resource, controller |  RL, LSTM |  yes |  ICT supply chain
TABLE I: Comparative analysis of intent-based architectural components

###  II-A Intent-Based Translator Design

####  II-A1 High-Level Language

The definition of a high level intent-based language is essential to allow the
user to express its requirements to the system. Multiple works in the
literature try to define this type of language, for instance, Riftadi et al
[8] employed Nile, a network intent language that provides an intermediate
layer between natural language and lower-level policies. In [15] a high level
declarative policy language is proposed to control the architecture of
software-defined networking (SDN) based on the notion of functional reactive
programming where policies are used to state the change in network conditions.
In [5], it is defined an Open Software Defined Framework (OSDF) policy
language used to provide a high level API in order to express network
requirements by managers and network administrators. [16] uses a Controlled
Natural Languages (CNL), which sets restrictions on inputs by the definition
of grammar model. In this case, CLN was employed to control blockchain
selection queries, as mapping requirements to a particular blockchain needs
expertise and knowledge about the existing blockchain implementations. Recent
development of natural language processing (NLP) algorithms [10, 12] show how
voice-assisted technologies have been used to manage the network as well.

####  II-A2 Interfaces

The initial works on intent-based interfaces support imply their abilities of
managing internal and external interactions independently from the underlying
network technologies. In the approach by Szyrkowiec et al [17], the intents
issued by a client application, are submitted through a representational state
transfer (REST) northbound interface (NBI). Similarly, in the work by Zeydan
et al [18], intent based NBI of NIC allows descriptive way to obtain what is
desired from the network using any kind of protocols, including OpenFlow, BGP,
and Netconf. In the work by Pham et al [11], the controller’s NBI and core
services are used as the atomic layer in the service orchestration because the
application wants to retrieve abstract network topology provided by the
controller. The intent layer in Han et al [9] provides extended North Bound
Interface (NBI) that is used by various applications.

Mahtout et al [4], deployed EVIAN, an approach for using ML for intent-based
provisioning, the EVIAN client provides an intelligent chatbot interface, via
Slack, to users in order to understand their network requirements and science
needs.

####  II-A3 Manager/Parser

In the [2], authors have recognized a core component of IBNs, which interacts
with most of the other components, naming it an Intent Processor. This core
component is responsible for accepting high level intents expressed by users
(applications) using a semantic language or through other appropriate user
interfaces. For instance, it can receive input in plain English and, in turn,
it can translate that into say, Algebraic Modelling Language, which can be
easily understood by a mathematical solver. In [5], the similar functionality
has been been identified under the name Policy Parser, responsible for
analyzing application-based policies.

####  II-A4 Policy Builder

In the work by Abbas et al [6] on creating core network slices through IBN,
the policy configurator extracts the information related to resources from the
graph, provided by the intent manager and converts the specifications into a
slice template. That slice template is according to the network orchestrator
acceptable format; for example, for the OSM core, the policy configurator can
give the slice template in form of a JSON file which contains the information
related to the deployment of network functions and their mapping.

####  II-A5 Compiler

The intent compiler is one of the key functionalities related to IBN system
design, with the task of converting the translated specifications or policies
into technical configurations expected to realise the intent. As such it has
been a focus of numerous research works. In [14] it is used to generate a low-
level configuration XML format from the low-level configuration data using
context-free grammar CFG. In [19] extracted data is converted to network
security function required data via Data Converter. In [8] an element called
P4 code generator is proposed to generate a working P4 program with the input
of a policy graph introduced by the user as an intent, while in [20] generic
constraints are converted into technology-specific ones, to be later applied
by network controller.

####  II-A6 Knowledge Base (KB)

The KB is a database that stores all the necessary information for the intent
based system, such as predefined policies, intent vocabularies and telemetry
used for improving policy building and intent translation process. In the
recent IBN works, KB has been mostly recognized under the name policy store
[7, 6, 5]. Its main tasks include storing and retrieving application-based
network policies that a user enters to the system. As such, it has been used
to access the current active policies in the system, allowing network
administrators to create, read, update, delete policies at runtime. In systems
like [8], the policy store represents an action library containing P4 code
templates, that can be parsed by the intent engine into a policy graph.

###  II-B Application in ICT systems

With an increased interest in the automation of network orchestration, the
emerging IBN approaches are being extensively analysed for its utilization in
numerous application fields, from vehicular, IoT and 5G scenarios, to its
potential in provisioning intent-based security policies and services.

One of the most interesting application assumed to benefit of IBN has been the
Vehicle Network Management with its distributed and self-organizing networks
composed of diverse and enormous number of high-speed vehicles; Here, IBN-
based policies are expected to manage vehicular networks towards an
intelligent network behavior according to the underlying applications or IoT
workflows, as well as help in minimizing the energy consumption and reducing
of network latency. An example can be seen in [21], where the authors designed
an IBN control framework over the SDN architecture in the vehicular edge
computing ecosystem, achieving the set goals of improving on network energy
efficient and latency reduction. network latency. High-level programming
abstractions and declarative languages [13] are posed to significantly benefit
mobile networks and applications as well, due to the increase in their
complexity over the years. This can be seen in the design and prototype
implementation of an Intent-based mobile backhauling interface for 5G networks
presented in [13] and in [22], where a wireless home automation system with
automatic spoken word into text commands translation was proposed.

The possibilities of smart network management and control with IBNs have also
been appealing with industrial solutions, with Huawei [23] developing an IBN
product called Intent-Driven Network (IDN) for its cloud solution. Likewise
OpenDayLight Controller Platform has been involved in numerous strategies on
enabling the utilization of intents with Network Intent Composition (NIC)
[24], Group-based Policy (GBP) intent abstraction [25], and NEMO [26]
language.

The full potential of IBN approach is probably most evident in network
security, where the security can be incorporated as an integral part of
intent-based networking across “intent,” “automation,” and “assurance”
functions. With their continuous loop frameworks, IBNs have proven to provide
ongoing protection and alignment to security policy and compliance
requirements [3]. Moreover, IBN implementation can provide traffic
segmentation based on business roles for connections to remote or cloud
computing sites to simplify threat protection and reduce the attack surface
[3]. A policy-driven security orchestration framework were presented in has
also been presented in [27] in SDN/NFV-aware IoT scenarios.

##  III System Design for ICT Supply Chains

Figure 1: Intent-based translator design for ICT supply chain networks

###  III-A Requirements Engineering

Deployment of IBN network orchestration and policy-based automation solutions
in ICT supply chain management carries a great potential. But, in order to
leverage this potential the future system design needs to adequately address
some of the main requirements which affect both of these domains, and can be
assumed to be critical in their integration.

Centralized/decentralized orchestrating approach: When opting for a
straightforward automation and deployment of intents, with the total control
being handed over to a single node with the complete knowledge of the network,
the choice would be a centralized network orchestrator. This approach comes at
a price of that the central node machine running the orchestrator would have
to have significant storage and processing power in order to maintain
information for potentially very large network. With a decentralized approach,
distributed orchestrators maintain an incomplete information of the network,
allowing higher levels of scalability but introducing its own challenge
regarding the efficient processing of intents across multiple orchestrating
nodes [1]. A similar discussion is led in assessing the advantages and
disadvantages of centralized and decentralized control and orchestration
strategies in todays complex ICT supply chains environments has to be
considered during the design phase of the management system. The more
traditional centralized approaches which often assume complete information
sharing among different SC organizational system components with the a top
hierarchical decision making authority have been noted to offer more stable
environments, faster decision making and easier coordination. On the other
side decentralized approach, with orchestration and decision making left to
individual SC organizational system components, benefits the SC efficiency.

Intent specification: One of the crucial issues in IBN has been on how to
actually specify high-level intents, and here there is still a lack of
standardized solutions. And while there are ongoing efforts towards improving
the ways intents can be described and approaching the naturally spoken
languages (e.g NEMO and INDIRA), how this will relate to translating intent-
based policies is to be discussed, but a general language model will be
required in order to provide a scalable solution for different ICT supply
chain systems.

AI technologies: Both IBN and ICT supply chains have been greatly influenced
by the advances in the domains of AI and ML. In case of AI-powered IBNs, the
performances have been improved by using AI in network analysis and intent
implementation as well as using AI-based heuristics for solving constrained
optimization problems. Here the challenge lies in a critical issue for any
AI/ML algorithm - the lack of data. To tackle this we can assume that IBN will
have to operate for certain amount of time to have enough data available, so
it will be necessary to assess how to use data collected from the pre-IBN
deployment. Moreover, when leveraging the potential of AI techniques for IBN
applications in supply chain management (SCM) there are two additional
questions that need to be addressed. Firstly, the question of leveraging the
potential of natural language processing (NLP) [28] in improving human–machine
interactions in SCM should be addressed. Secondly, as supply chains are
particularly important from industrial networks perspective, in case of an
incident it is a requirement to obtain the explanation behind it, which can be
difficult in case it is a consequence of an IBN automated algorithm reaction
[2].

Closed-loop verification: With closed-loop verification method, IBNs maintain
the enforced network state and configuration through intents, which can be
expressed as series of smaller network element configurations [29]. The
problem lies in the complexity of this process, which assumes constant system
and network configuration adjustments in order to satisfy desired system
requirements and can result in autonomous remediation errors or unwanted
delays. The complex system architectures as the ones behind SC scenarios will
need to improve on this method, as a seemingly small configuration adjustments
can affect the entire network across.

Single vs Multi-Domain: Majority of the IBN deployment and research work
focuses on single-domain network approach, with the orchestrator
centralization. However, this kind of approach is insufficient in addressing
the SCM scenarios, as they normally extend across multiple companies,
decision-makers and network domains. To address this particular issue there is
a need for shifting towards a more multi-domain oriented IBN approach, where
an intent framework will be implemented in a manner that would ensure end-to-
end network orchestration across multiple potential SC domains.

###  III-B Architecture Design

The intent based translator aims to automate the configuration management of
ICT supply chain networks. Thus, we propose a general software architecture
for the intent abstraction layer, as shown in Fig.1. Supply chains include
various organizational partners that own different parts of the supply chain
and need to cooperate in order to successfully deliver a final product or
service. The management solution of such network will allow different network
administrators to easily configure and control the network via high-level
intents. A plain text, that is restricted by a predefined high-level language,
is input by the user, using a dashboard that communicates the information via
an intent-based interface. This interface is developed as to facilitate the
human machine interactions, considering the ability to control the
infrastructure using high-level intents for security assessment, network
service management and control. From here, as shown in Fig.1, high-level
intents are translated and parsed by the intent manager into a structured
format, using natural language processing techniques. A knowledge base is then
used to store the intent vocabularies in a json format, where it is used by
the intent manager to determine the intent specifications. The vocabularies
inside the knowledge base can be enriched dynamically, improving the system
scalability and capability to manage supply chains of different contexts.
After producing a structured format of an intent with user requirement
specifications, the intent is sent to a so-called policy configurator
component, responsible on matching the user requirements with existing
policies that can be applied in the infrastructure. This component will then
check the existing policies stored in the knowledge base in a defined policy
store that can be updated and accessed by developers and experts in order to
define new policies. The created policies are mapped by the policy
configurator into a lower level format, after checking for possible conflicts
among them. The policy configurator is also responsible for triggering the
network controller, using defined APIs. In proposed system we also consider
the possibility of smartly reconfiguring the network based on insights from
the system. Metrics are collected by monitoring tools from all the network
elements, and used to periodically update the knowledge base, which includes
telemetry as well as storing policies and intent vocabularies. A component
called AI Engine is responsible on analyzing the data coming from the system,
predicting the network status, and automatically configuring the policies.

###  III-C Example of Supply Chain Assets Access Control

An overview of the SC abstraction proposed in our example is illustrated in
Fig. 2. Organizations are in the top of the hierarchy and represent companies,
consortiums or institutions that belongs to the supply chain. Realms consist
of an environment within an organization that shares common security policies
and rules. Domains are a set of assets with certain relationship, for
instance, they share the same network, location or infrastructure.

Figure 2: Supply chain abstraction building block

####  III-C1 Intent Language Definition

We define a simple grammar for SC asset access control as an intent language,
based on a context-free grammar standard: the Extended Backus-Naur Form (EBNF)
[30]. The defined grammar in Grammar 1 is a CNL that allows the system
administrator to grant access to users for a certain asset in the system,
restricted to specific organization, realm, or domain at a certain working
time. This proposal is a general definition of the grammar, and can be
enriched with more details related to the nature of an asset, the composition
of a supply chain, the structure of a company, and finally the type of users.

¡intent¿ ::= ¡users¿ is ¡permission¿ to access to ¡asset¿ ¡spot¿ at
¡timeframes¿

¡users¿ ::= ¡user¿ (, and ) ¡user¿

¡user¿ ::= [a-z0-9]+

¡permission¿ ::= allowed — blocked

¡asset¿ ::= domain — realm — organization

¡spot¿ ::= [a-z0-9]+

¡timeframes¿ ::= ¡timeframe¿ (, and ) ¡timeframe¿

¡timeframe¿ ::= ¡morning shift¿ — ¡late shift¿ — ¡night shift¿

¡morning shift¿ ::= [6:00-13:59]

¡late shift¿ ::= [14:00-21:59]

¡night shift¿ ::= [22:00-5:59]

List of Grammars 1 SC Assets Access Control Intent

####  III-C2 Intent Manager

It parses the input plain text into a json file using the defined high-level
language and regular expressions (regex approach). The parsing process can be
improved using ML techniques to allow the user to employ multiple vocabulary
options. The intent manager reads the intent store from the knowledge base and
try to match the existing vocabularies and structures with the input text.
Finally a json file with user specifications is generated and sent to the
policy configurator.

####  III-C3 Policy Configurator

It matches the requirements specifications received from the intent manager
with the existing policies stored in the policy store. It solves conflict
between intents and generates valid policies. Conflicts can happen between
policies while compiling intent specifications. A user can be restricted from
access to assets in realm _o1-r1_ using a following intent : _user-x is not
allowed to access to realm o1-r1_. The policy configurator will map the intent
into high level policies such as:
_check user-x in database of Users
block user-x to access assets in o1-r1
alert admin in o1_
When the administrator wants to give the user access to assets in organization
_o1_ , the following intent can be written: _user-x is allowed to access to
organization o1_. This intent will create a conflict with the previous intent
during the policy translation, as it will allow the user to access to realm
_o1-r1_ , which is inside the organization _o1_. In order to solve this issue,
conflict checking needs to be considered and implemented. A solution to the
mentioned conflict can be translated to the following high-level policies:
_check user-x in database of Users
allow user-x to access assets in o1 except o1-r1
alert admin in o1_
Other conflicts should be analysed considering the organizational structure,
the existing policies and the possible human intervention.

####  III-C4 Knowledge Base

It consists of three main components: intent store, policy store, and
telemetry. The intent store is used to store a dictionary of vocabularies and
options that are required to parse intents. The policy store represents the
policy graphs defined by experts, in our case the policy graph represent the
access to assets based on the organizational hierarchy. And finally the
telemetry is used to store the needed metrics to detect anomalies in the
system.

####  III-C5 AI Engine

NLP techniques, such as long-short time memory (LSTM), can be implemented to
enhance the proposed CNL and provide more flexibility to the user while
writing the intents. The parsing of the input plain text considers a rich
dictionary of vocabularies and expressions that can mean the same input
specifications. AI techniques can also be used to enhance the policy
configuration by solving policy conflicts. Metrics collected from system and
stored in the knowledge base can be used to predict if a user has to be
blocked or not based on his monitored behavior, and ML techniques are a great
candidate for such task of pattern recognition and security assessment.

##  IV Conclusion

This paper proposes a general design of an intent-based translation system
that manages ICT supply chain networks, considering AI enabling technologies
to ensure a valid translation of user requirements, network security and
scalability. With our system design we intend to move the complexity of
network management, which requires very high level of expertise from the
system administrator by automating policy building process through AI and ML
techniques. We analysed an use case example of an ICT supply chain management
task - access control policy configuration using high level intents. Also, we
discussed the engineering requirements and challenges in integrating IBN
approach with ICT supply chains, in terms of network orchestration complexity,
the application of AI techniques, and the ambiguity of defining high level
intents for ICT supply chain policies.

## Acknowledgment

This work is partially funded by European Commission under the H2020-952644
contract for project FISHY: A coordinated framework for cyber resilient supply
chain systems over complex ICT infrastructures.

## References

  * [1] E. Shry, “The state of intent-based networking,” 2019. [Online]. Available: http://www.cse.wustl.edu/~jain/cse570-19/ftp/intent/index.html
  * [2] B. K. Saha _et al._ , “Intent-based networks: An industrial perspective,” in _Proceedings of the 1st International Workshop on Future Industrial Communication Networks_ , 2018, pp. 35–40.
  * [3] B. Laliberte, “The journey to intent-based networking: Ten key principles for accelerating adoption,” 2018.
  * [4] H. Mahtout _et al._ , “Using machine learning for intent-based provisioning in high-speed science networks,” in _Proceedings of the 3rd International Workshop on Systems and Network Telemetry and Analytics_ , 2020, pp. 27–30.
  * [5] D. Comer and A. Rastegatnia, “Osdf: An intent-based software defined network programming framework,” in _2018 IEEE 43rd Conference on Local Computer Networks (LCN)_ , 2018, pp. 527–535.
  * [6] K. Abbas _et al._ , “Slicing the core network and radio access network domains through intent-based networking for 5g networks,” _Electronics_ , vol. 9, no. 10, p. 1710, 2020.
  * [7] T. A. Khan _et al._ , “Intent-based orchestration of network slices and resource assurance using machine learning,” in _NOMS 2020-2020 IEEE/IFIP Network Operations and Management Symposium_. IEEE, 2020, pp. 1–2.
  * [8] M. Riftadi and F. Kuipers, “P4i/o: Intent-based networking with p4,” in _2019 IEEE Conference on Network Softwarization (NetSoft)_. IEEE, 2019, pp. 438–443.
  * [9] Y. Han _et al._ , “An intent-based network virtualization platform for sdn,” in _2016 12th International Conference on Network and Service Management (CNSM)_. IEEE, 2016, pp. 353–358.
  * [10] A. Chaudhari _et al._ , “Vivonet: Visually-represented, intent-based, voice-assisted networking,” _arXiv preprint arXiv:1904.03228_ , 2019.
  * [11] M. Pham and D. B. Hoang, “Sdn applications-the intent-based northbound interface realisation for extended applications,” in _2016 IEEE NetSoft Conference and Workshops (NetSoft)_. IEEE, 2016, pp. 372–377.
  * [12] M. Jain _et al._ , “Intent-based, voice-assisted, self-healing sdn framework,” _Journal of Network Communications and Emerging Technologies (JNCET)_ , vol. 10, no. 2, 2020.
  * [13] T. Subramanya _et al._ , “Intent-based mobile backhauling for 5g networks,” in _2016 12th International Conference on Network and Service Management (CNSM)_. IEEE, 2016, pp. 348–352.
  * [14] C. Chung and J. P. Jeong, “A design of iot device configuration translator for intent-based iot-cloud services,” in _2020 22nd International Conference on Advanced Communication Technology (ICACT)_. IEEE, 2020, pp. 52–56.
  * [15] A. Voellmy _et al._ , “Procera: A language for high-level reactive network control,” in _Proceedings of the First Workshop on Hot Topics in Software Defined Networks_ , ser. HotSDN ’12. New York, NY, USA: Association for Computing Machinery, 2012, p. 43–48. [Online]. Available: https://doi.org/10.1145/2342441.2342451
  * [16] E. J. Scheid _et al._ , “A controlled natural language to support intent-based blockchain selection,” in _2020 IEEE International Conference on Blockchain and Cryptocurrency (ICBC)_. IEEE, 2020, pp. 1–9.
  * [17] T. Szyrkowiec _et al._ , “Automatic intent-based secure service creation through a multilayer sdn network orchestration,” _Journal of Optical Communications and Networking_ , vol. 10, no. 4, pp. 289–297, 2018.
  * [18] E. Zeydan and Y. Turk, “Recent advances in intent-based networking: A survey,” in _2020 IEEE 91st Vehicular Technology Conference (VTC2020-Spring)_. IEEE, 2020, pp. 1–5.
  * [19] J. Yang and J. P. Jeong, “An automata-based security policy translation for network security functions,” in _2018 International Conference on Information and Communication Technology Convergence (ICTC)_. IEEE, 2018, pp. 268–272.
  * [20] P. Sköldström _et al._ , “Dismi-an intent interface for application-centric transport network services,” in _2017 19th International Conference on Transparent Optical Networks (ICTON)_. IEEE, 2017, pp. 1–4.
  * [21] A. Singh _et al._ , “Intent-based network for data dissemination in software-defined vehicular edge computing,” _IEEE Transactions on Intelligent Transportation Systems_ , 2020.
  * [22] S. Amrutha _et al._ , “Voice controlled smart home,” _International Journal of Emerging Technology and Advanced Engineering_ , vol. 5, no. 1, pp. 272–275, 2015.
  * [23] Huawei, “Huawei intent-driven network white paper,” 2019, accessed: 04.05.2021. [Online]. Available: https://bit.ly/3nOKPue
  * [24] OpenDayLight, “Network intent composition (nic) user guide,” 2016, accessed: 04.05.2021. [Online]. Available: https://bit.ly/3tiN36b
  * [25] ——, “Group based policy user guide,” 2016, accessed: 04.05.2021. [Online]. Available: https://bit.ly/3b26dXv
  * [26] onosproject, “Nemo language,” 2016, accessed: 04.05.2021. [Online]. Available: https://bit.ly/3eShzOX
  * [27] A. M. Zarca _et al._ , “Semantic-aware security orchestration in sdn/nfv-enabled iot systems,” _Sensors_ , vol. 20, no. 13, p. 3622, 2020\.
  * [28] R. Toorajipour _et al._ , “Artificial intelligence in supply chain management: A systematic literature review,” _Journal of Business Research_ , vol. 122, pp. 502–517, 2021. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S014829632030583X
  * [29] L. Pang _et al._ , “A survey on intent-driven networks,” _IEEE Access_ , vol. 8, pp. 22 862–22 873, 2020.
  * [30] R. S. Scowen, “Generic base standards,” in _Proceedings 1993 Software Engineering Standards Symposium_. IEEE, 1993, pp. 25–34.


# Hybrid explicit-implicit learning for multiscale problems with time
dependent source

Yalchin Efendiev111Department of Mathematics, Texas A&M University, College
Station, TX 77843, USA, Wing Tat Leung222Department of Mathematics, City
University of Hong Kong, Hong Kong SAR, Wenyuan Li333Department of
Mathematics, Texas A&M University, College Station, TX 77843, USA, Zecheng
Zhang444Department of Mathematical Sciences, Carnegie Mellon University,
Pittsburgh, PA 15213, USA

## Abstract

The splitting method is a powerful method for solving partial differential
equations. Various splitting methods have been designed to separate different
physics, nonlinearities, and so on. Recently, a new splitting approach has
been proposed where some degrees of freedom are handled implicitly while other
degrees of freedom are handled explicitly. As a result, the scheme contains
two equations, one implicit and the other explicit. The stability of this
approach has been studied. It was shown that the time step scales as the
coarse spatial mesh size, which can provide a significant computational
advantage. However, the implicit solution part can be expensive, especially
for nonlinear problems. In this paper, we introduce modified partial machine
learning algorithms to replace the implicit solution part of the algorithm.
These algorithms are first introduced in ‘HEI: Hybrid Explicit-Implicit
Learning For Multiscale Problems’, where a homogeneous source term is
considered along with the Transformer, which is a neural network that can
predict future dynamics. In this paper, we consider time-dependent source
terms which is a generalization of the previous work. Moreover, we use the
whole history of the solution to train the network. As the implicit part of
the equations is more complicated to solve, we design a neural network to
predict it based on training. Furthermore, we compute the explicit part of the
solution using our splitting strategy. In addition, we use Proper Orthogonal
Decomposition based model reduction in machine learning. The machine learning
algorithms provide computational saving without sacrificing accuracy. We
present three numerical examples which show that our machine learning scheme
is stable and accurate.

##  1 Introduction

Many physical problems vary at multiple space and time scales. The solutions
of these problems are difficult to compute as we need to resolve the time and
space scales. For example, porous media flow and transport occur over many
space and time scales. Many multiscale and splitting algorithms have been
proposed to solve these problems [28, 18, 41]. In this paper, we would like to
combine machine learning and splitting/multiscale approaches to predict the
dynamics of complex forward problems. Our goal is to perform machine learning
on a small part of the solution space, while computing “the rest” of the
solution fast. Our proposed splitting algorithms allow doing it, which we will
describe below.

Our approaches use multiscale methods in their construction to achieve coarser
time stepping. Multiscale methods have extensively been studied in the
literature. For linear problems, homogenization-based approaches [23, 39, 40],
multiscale finite element methods [23, 31, 35], generalized multiscale finite
element methods (GMsFEM) [12, 13, 14, 17, 21], Constraint Energy Minimizing
GMsFEM (CEM-GMsFEM) [15, 16, 7], nonlocal multi-continua (NLMC) approaches
[19], metric-based upscaling [44], heterogeneous multiscale method [20],
localized orthogonal decomposition (LOD) [30], equation-free approaches [45,
46], multiscale stochastic approaches [33, 34, 32], and hierarchical
multiscale methods [5] have been studied. Approaches such as GMsFEM and NLMC
are designed to handle high-contrast problems. For nonlinear problems, we can
replace linear multiscale basis functions with nonlinear maps [26, 27, 22].

Splitting methods are introduced when solving dynamic problems [3, 47, 29].
There is a wide range of applications for splitting methods, e.g. splitting
physics, scales, and so on. Splitting methods allow simplifying the large-
scale problem into smaller-scale problems. We couple them and obtain the final
solution after we solve the smaller-scale problems. Recently, a partially
explicit splitting method [28, 18, 24] is proposed for linear and nonlinear
multiscale problems. In this approach, we treat some solution parts implicitly
while treating the other parts explicitly. The stability analysis shows that
the time step only depends on the explicit part of the solution. With the
careful splitting of the space, we are able to make the time step scale as the
coarse mesh size and independent of contrast.

To illustrate the main idea of our methods, we consider

$u_{t}+f(u)=g(x,t).$

In this paper, we combine machine learning and splitting approaches (such as
the partially explicit method). In the splitting algorithm, we write the
solution as $u=u_{1}+u_{2}$, where $u_{1}$ is computed implicitly and $u_{2}$
is computed explicitly. In our previous work, the Transformer has been used to
solve this kind of problem in [25, 11]. In [25], we consider zero source term
and use the Transformer to predict the future dynamics of the solution. In our
current paper, we consider full dynamics and time-dependent source terms. The
latter occurs in porous media applications and represents the rates of wells,
which vary in time. The goal in these simulations is to compute the solutions
at all time steps for a given time-dependent source term in an efficient way,
which is not addressed before.

Many machine learning and multiscale methods are designed for solving
multiscale parametric PDEs [11, 48, 49, 9, 43, 42]. The advantage of machine
learning is that it treats PDEs as an input-to-output map. By training, and
possibly, using the observed data, machine learning avoids computing solutions
for each input (such as parameters in boundary conditions, permeability fields
or sources). In this paper, we use machine learning to learn the difficult-to-
compute part of the solution that is, in general, computed implicitly. This is
a coarse-grid part of the solution, which can also be observed given some
dynamic data. The splitting algorithm gives two coupled equations. The first
equation for $u_{1}$ is implicit while the second equation for $u_{2}$ is
explicit. We design a neural network to predict the first component of the
solution. Then, we calculate the second component using the discretized
equations. The inputs of our neural network are the parameters that
characterize the source term of the equation. The outputs are the whole
history of $u_{1}$, i.e. $u_{1}$ at all time steps. The output dimension of
our neural network is the dimension of $u_{1}$ times the number of time steps,
which is considerably large. Therefore, we use Proper Orthogonal Decomposition
(POD) [48, 10, 38] to reduce the dimension of the solution space. We
concatenate all solution vectors of every training input together to generate
a large-scale matrix. After doing the singular value decomposition, we select
several dominant eigenvectors to form the POD basis vectors. Thus, the output
becomes the coordinate vectors based on the POD basis vectors and its
dimension is greatly reduced, so is the number of parameters inside the neural
networks.

In this paper, we show several numerical examples to illustrate our learning
framework which is based on the partially explicit scheme. In our numerical
examples, the diffusion term $f(u)$ is $-div(\kappa(x,u)\nabla u)$ and the
reaction term is $g(x,t)$. $g(x,t)$ is a source term that is nonzero on
several coarse blocks and is zero elsewhere. We choose such $g(x,t)$ as it
simulates the wells in petroleum applications. By comparing with the solution
which is obtained using fine grid basis functions, we find that the relative
$L^{2}$ error is small and similar to the scheme where we calculate both
components. Thus, the learning scheme can achieve similar accuracy while
saving computation.

This paper is organized as follows. We present preliminaries and review the
partially explicit splitting scheme in Section 2. In Section 3, we show
details about our learning framework. We present numerical experiments in
Section 4.

##  2 Problem Setting

We consider the following parabolic partial differential equation

$\displaystyle u_{t}+f(u)=g(x,t).$ (2.1)

Let $V$ be a Hilbert space. In this equation, $f=\cfrac{\delta F}{\delta u}\in
V^{\ast}$ with $F$ being the variational derivative of energy
$F(u):=\int_{\Omega}E(u)$. In this paper, we use $(\cdot,\cdot)$ to denote
$(\cdot,\cdot)_{V^{*},V}$ for simplicity. We assume that $f$ is contrast
dependent (linear or nonlinear). Here, $g(x,t)$ is a time dependent source
term with parameters. The weak formulation of Equation (2.1) is the following.
We want to find $u\in V$ such that

$\displaystyle(u_{t},v)+(f(u),v)=(g(x,t),v)\quad\forall v\in V.$ (2.2)

The $(\cdot,\cdot)$ is the $L^{2}$ inner product. We assume that

$\displaystyle(f(v),v)>0,\quad\forall v\in V.$ (2.3)

######  Example 1.

In the heat equation, $F(u^{\prime})=\cfrac{1}{2}\int_{\Omega}\kappa|\nabla
u^{\prime}|^{2}$. Let $V=H^{1}(\Omega)$. We have $f(u)=\cfrac{\delta F}{\delta
u^{\prime}}(u)\in(H^{1}(\Omega))^{*}=H^{1}(\Omega)$. Then,

$(\cfrac{\delta F}{\delta u^{\prime}}(u),v)=\int_{\Omega}\kappa\nabla
u\cdot\nabla v.$

Therefore, we have

$f(u)=-\nabla\cdot(\kappa\nabla u).$

The weak formulation (2.2) becomes

$(u_{t},v)-(\nabla\cdot(\kappa\nabla u),v)=(g(x,t),v).$

The standard approach to solving this problem is the finite element method.
Let $V_{H}$ be the finite element space and $u_{H}\in V_{H}$ is the numerical
solution satisfying

$\displaystyle(u_{H,t},v)+(f(u_{H}),v)=(g(x,t),v)\quad\forall v\in V_{H}.$
(2.4)

The Backward Euler (implicit) time discretization for this Equation is

$\displaystyle(\frac{u^{n+1}_{H}-u^{n}_{H}}{\Delta
t},v)+(f(u^{n+1}_{H}),v)=(g(x,t^{n+1}),v)\quad\forall v\in V_{H}.$ (2.5)

The Forward Euler (explicit) time discretization for Equation (2.4) is

$\displaystyle(\frac{u^{n+1}_{H}-u^{n}_{H}}{\Delta
t},v)+(f(u^{n}_{H}),v)=(g(x,t^{n}),v)\quad\forall v\in V_{H}.$ (2.6)

One can find details about the stability analysis of the Backward and Forward
Euler method in [18].

###  2.1 Partially Explicit Scheme

In this paper, we use the partially explicit splitting scheme introduced in
[18]. We split the space $V_{H}$ to be the direct sum of two subspaces
$V_{H,1}$ and $V_{H,2}$, i.e. $V_{H}=V_{H,1}\oplus V_{H,2}$. Then the solution
$u_{H}=u_{H,1}+u_{H,2}$ with $u_{H,1}\in V_{H,1}$ and $u_{H,2}\in V_{H,2}$
satisfies

$\displaystyle((u_{H,1}+u_{H,2})_{t},v_{1})+(f(u_{H,1}+u_{H,2}),v_{1})=(g(x,t),v_{1})\quad\forall
v_{1}\in V_{H,1},$
$\displaystyle((u_{H,1}+u_{H,2})_{t},v_{2})+(f(u_{H,1}+u_{H,2}),v_{2})=(g(x,t),v_{2})\quad\forall
v_{2}\in V_{H,2}.$

We use the partially explicit time discretization for this problem. Namely, we
consider finding $\\{u_{H}^{n}\\}_{n=0}^{K}$ (with
$u_{H}^{n}=u_{H,1}^{n}+u_{H,2}^{n},\;\forall n$) such that

$\displaystyle(\frac{u_{H,1}^{n+1}-u_{H,1}^{n}}{\Delta
t}+\frac{u_{H,2}^{n}-u_{H,2}^{n-1}}{\Delta
t},v_{1})+(f(u_{H,1}^{n+1}+u_{H,2}^{n}),v_{1})=(g(x,t^{n}),v_{1})\quad\forall
v_{1}\in V_{H,1},$ (2.7)
$\displaystyle(\frac{u_{H,1}^{n}-u_{H,1}^{n-1}}{\Delta
t}+\frac{u_{H,2}^{n+1}-u_{H,2}^{n}}{\Delta
t},v_{2})+(f(u_{H,1}^{n+1}+u_{H,2}^{n}),v_{2})=(g(x,t^{n}),v_{2})\quad\forall
v_{2}\in V_{H,2}.$ (2.8)

We refer to [28, 18] for the details of convergence assumptions and stability
analysis of the above scheme.

####  2.1.1 Construction of $V_{H,1}$ and $V_{H,2}$

In this section, we present how to construct the two subspaces for our
partially explicit splitting scheme. We will first present the CEM method
which is used to construct $V_{H,1}$. Next, with the application of eigenvalue
problems, we build $V_{H,2}$. In the following discussion, let
$S\subset\Omega$ and then we define $V(S)=H_{0}^{1}(S)$.

CEM method

We introduce the CEM finite element method in this section. Let
$\mathcal{T}_{H}$ be a coarse grid partition of $\Omega$. For
$K_{i}\in\mathcal{T}_{H}$, we need to construct a collection of auxiliary
basis functions in $V(K_{i})$. We assume that $\\{\chi_{i}\\}$ are basis
functions that form a partition of unity, e.g., piecewise linear functions,
see [4]. We first find $\lambda_{j}^{(i)}$ and $\psi_{j}^{(i)}$ such that the
following equation holds:

$\displaystyle\int_{K_{i}}\kappa\nabla\psi_{j}^{(i)}\cdot\nabla
v=\lambda_{j}^{(i)}s_{i}(\psi_{j}^{(i)},v),\quad\forall v\in V(K_{i}),$

where

$\displaystyle
s_{i}(u,v)=\int_{K_{i}}\tilde{\kappa}uv,\quad\tilde{\kappa}=\kappa
H^{-2}\;\text{or}\;\tilde{\kappa}=\kappa\sum_{i}\left|\nabla\chi_{i}\right|^{2}.$

With rearrangement if necessary, we select the $L_{i}$ eigenfunctions of the
corresponding first $L_{i}$ smallest eigenvalues. We define

$V_{aux}^{(i)}:=\text{span}\\{\psi_{j}^{(i)}:\;1\leq j\leq L_{i}\\}.$

Define a projection operator $\Pi:L^{2}(\Omega)\mapsto V_{aux}\subset
L^{2}(\Omega)$

$s(\Pi u,v)=s(u,v),\quad\forall v\in
V_{aux}:=\sum_{i=1}^{N_{e}}V_{aux}^{(i)},$

with $s(u,v):=\sum_{i=1}^{N_{e}}s_{i}(u|_{K_{i}},v|_{K_{i}})$ and $N_{e}$
being the number of coarse elements. We denote an oversampling domain which is
several coarse blocks larger than $K_{i}$ as $K_{i}^{+}$[15]. For every
auxiliary basis $\psi_{j}^{(i)}$, we solve the following two equations and
find $\phi_{j}^{(i)}\in V(K_{i}^{+})$ and $\mu_{j}^{(i)}\in V_{aux}$

$\displaystyle a(\phi_{j}^{(i)},v)+s(\mu_{j}^{(i)},v)$
$\displaystyle=0,\quad\forall v\in V(K_{i}^{+}),$ $\displaystyle
s(\phi_{j}^{(i)},\nu)$ $\displaystyle=s(\psi_{j}^{(i)},\nu),\quad\forall\nu\in
V_{aux}(K_{i}^{+}).$

Then we can define the CEM finite element space as

$\displaystyle V_{cem}$ $\displaystyle:=\text{span}\\{\phi_{j}^{(i)}:\;1\leq
i\leq N_{e},1\leq j\leq L_{i}\\}.$

We define $\tilde{V}:=\\{v\in V:\;\Pi(v)=0\\}$ and we will need it when we
construct $V_{H,2}$ in the next subsection.

Construction of $V_{H,2}$

We build the second subspace $V_{H,2}$ based on $V_{H,1}$ and $\tilde{V}$.
Similarly, we first have to construct the auxiliary basis functions by solving
the eigenvalue equations. For every coarse element $K_{i}$, we search for
$(\xi_{j}^{(i)},\gamma_{j}^{(i)})\in(V(K_{i})\cap\tilde{V})\times\mathbb{R}$
such that

$\displaystyle\int_{K_{i}}\kappa\nabla\xi_{j}^{(i)}\cdot\nabla v$
$\displaystyle=\gamma_{j}^{(i)}\int_{K_{i}}\xi_{j}^{(i)}v,\;\ \forall v\in
V(K_{i})\cap\tilde{V}.$ (2.9)

We choose the $J_{i}$ eigenfunctions corresponding to the smallest $J_{i}$
eigenvalues. The second type of auxiliary space is defined as
$V_{aux,2}:=\text{span}\\{\xi_{j}^{(i)}:1\leq i\leq N_{e},1\leq j\leq
J_{i}\\}$. For every auxiliary basis $\xi_{j}^{(i)}\in V_{aux,2}$, we look for
a basis function $\zeta_{j}^{(i)}\in V(K_{i}^{+})$ such that for some
$\mu_{j}^{(i),1}\in V_{aux}$, $\mu_{j}^{(i),2}\in V_{aux,2}$, we have

$\displaystyle a(\zeta_{j}^{(i)},v)+s(\mu_{j}^{(i),1},v)+(\mu_{j}^{(i),2},v)$
$\displaystyle=0,\quad\forall v\in V(K_{i}^{+}),$ (2.10) $\displaystyle
s(\zeta_{j}^{(i)},\nu)$ $\displaystyle=0,\quad\forall\nu\in V_{aux},$ (2.11)
$\displaystyle(\zeta_{j}^{(i)},\nu)$
$\displaystyle=(\xi_{j}^{(i)},\nu),\quad\forall\nu\in V_{aux,2}.$ (2.12)

We define

$V_{H,2}=\text{span}\\{\zeta_{j}^{(i)}|\;1\leq i\leq N_{e},\;1\leq j\leq
J_{i}\\}.$

##  3 Machine Learning Approach

In Equation (2.1), we consider $g(x,t)$ as a source term with varying
parameters representing temporal changes. For every parameter, we want to find
the corresponding solution $u(x,t)$. For the traditional finite element
methods, we need to construct meshes and basis functions, assemble the
stiffness and mass matrices and solve the matrices equations, which makes the
computation cost large, especially for the nonlinear PDEs as we have to
assemble the stiffness matrix at every time step. Designing machine learning
algorithms for the partially explicit scheme can alleviate this problem. Once
the neural network is trained, the computation of solution becomes the
multiplication of matrices and applying activation functions, which is more
efficient in contrast to the finite element methods.

The Equation (2.7) is implicit for $u^{n+1}_{H,1}$, which makes the
computations more complicated. We design a machine learning algorithm to
predict $u^{n+1}_{H,1}$ and then we compute $u^{n+1}_{H,2}$ using Equation
(2.8).

Figure 3.1: Flowchart.

We present a flowchart in Figure 3.1 to show the workflow of our method. We
first need to generate the dataset which will be used to train the neural
network. Let $w_{m}=(w_{m1},w_{m2},\cdots,w_{mk})$ be the parameters in
$g_{w}(x,t)$. We set $w_{mi}\in[L,U]$ for $i=1,2,\cdots,k$ with $L$ and $U$
being the lower bound and upper bound, respectively. We generate the training
space $\\{w_{m}\\}_{m=1}^{M}$ by drawing samples from a uniform distribution
over the interval $[L,U]$. For every $w_{m}$ in the training space, we
calculate its corresponding partially explicit solution using Equations (2.7)
and (2.8). Let $N$ be the number of time steps and $M$ be the size of our
training sample space. We denote this solution space as
$\\{\\{u_{H,m}^{n}\\}_{n=2}^{N}\\}_{m=1}^{M}$, where
$u_{H,m}^{n}=u_{H,m,1}^{n}+u_{H,m,2}^{n}$. We remark that the reason why $n$
starts from $2$ is that we need $u^{0}$ and $u^{1}$ to compute $u^{2}$ in the
partially explicit scheme. The $u^{1}$ is computed by the Backward Euler
scheme using basis functions from both $V_{H,1}$ and $V_{H,2}$. We only need
$\\{\\{u_{H,m,1}^{n}\\}_{n=2}^{N}\\}_{m=1}^{M}$ for our machine learning
algorithm. Before we move on, we discuss a useful tool called Proper
Orthogonal Decomposition [6].

For a given set of vectors $\\{x_{i}\\}_{i=1}^{c}$ with
$x_{i}\in\mathbb{R}^{d}$, we want to construct a smaller dimensional vector
space that can represent the space spanned by $\\{x_{i}\\}_{i=1}^{c}$
accurately, i.e., we want to find orthonormal vectors $\\{y_{j}\\}_{j=1}^{l}$
that solve

$\displaystyle\min\limits_{y_{1},\cdots,y_{l}}\sum_{j=1}^{c}\|x_{j}-\sum_{i=1}^{l}\langle
x_{j},y_{l}\rangle_{l^{2}}y_{l}\|^{2}_{l^{2}},$ (3.1)

where $l\leq r$ and $r$ is the rank of $(x_{1},x_{2},\cdots,x_{c})$. We apply
the singular value decomposition for this. Next, we briefly remind the POD
procedure [6]. We first define a matrix $X=(x_{1},x_{2},\cdots,x_{c})$ whose
dimension is $d\times c$. Then $X^{T}X$ and $XX^{T}$ are $c\times c$ and
$d\times d$ dimensional matrices. Consider the eigenvalues and eigenvectors of
$X^{T}X$ as well as $XX^{T}$,

$X^{T}Xz_{i}=\sigma_{i}^{2}z_{i},\quad XX^{T}y_{i}=\sigma_{i}^{2}y_{i}.$

Next, we define

$Y=(y_{1},\cdots,y_{d}),\quad Z=(z_{1},\cdots,z_{c}).$

We now have

$Y^{T}XZ=\Sigma,$

where $\Sigma$ is the diagonal matrix whose main diagonal entries are
$\sigma_{i}^{2}$. After doing rearrangement if necessary, we may assume

$\sigma_{1}^{2}\geq\sigma_{2}^{2}\geq\cdots>0.$

The solution of (3.1) is given by the $l$ eigenvectors of $XX^{T}$
corresponding to the largest $l$ eigenvalues. Thus, the desired $l$
dimensional subspace is the subspace spanned by eigenvectors corresponding to
the largest $l$ eigenvalues.

There are $N\times M$ column vectors ($\in\mathbb{R}^{\dim(V_{H,1})}$) in the
solution space $\\{\\{u_{H,m}^{n}\\}_{n=2}^{N}\\}_{m=1}^{M}$. We would like to
output the solution $\\{u_{H}^{n}\\}_{n=2}^{N}$ at the same time so that we
don’t have to do the iteration for every time step recursively. In this case,
the output dimension is $\dim(V_{H,1})\times(N-1)$, which is too large
($\dim(V_{H,1})$ is $300$ in our examples). As a result, the number of
parameters in the neural network will also be large, which increases the
computation cost. We define

$\displaystyle
S=(u_{H,1,1}^{2},u_{H,1,1}^{3},\cdots,u_{H,1,1}^{N},u_{H,2,1}^{2},u_{H,2,1}^{3},\cdots,u_{H,2,1}^{N},\cdots,u_{H,M,1}^{2},u_{H,M,1}^{3},\cdots,u_{H,M,1}^{N}).$

Hence, $\dim(S)=\dim(V_{H,1})\times((N-1)\cdot M)$. After the singular value
decomposition to $S$, we define

$P=(p_{1},p_{2},\cdots,p_{l}),$

where $p_{i}\in\mathbb{R}^{\dim(V_{H,1})}$ is the eigenvector corresponding to
the $i$-th largest eigenvalue. We use $P$ to get a reduced dimensional
coordinate vectors for every $u_{H,m,1}^{n}$,i.e.

$\displaystyle
u_{H,m,1}^{n,\prime}=(P^{T}P)^{-1}P^{T}u_{H,m,1}^{n}\in\mathbb{R}^{l}$ (3.2)

Thus, the output dimension is reduced to $l\times(N-1)$. In the testing stage,
after we get the output $\\{u_{H,1}^{n,\prime}\\}_{n=2}^{N}$, we use $P$ to
get the coarse grid solution
$\\{u_{H,1}^{n}=Pu_{H,1}^{n,\prime}\\}_{n=2}^{N}$.

As discussed above, in our neural networks, the inputs are the parameters $w$
in $g_{w}(x,t)$ and outputs are $\\{u_{H,m,1}^{n,\prime}\\}_{n}$. That is

$\displaystyle\\{u_{H,m,1}^{n,\prime}\\}_{n}=\mathcal{N}(w),$

where $\mathcal{N}$ is the neural network. To better illustrate,

$\displaystyle Y_{1}=\phi_{1}(wW_{1}+b_{1}),$

where $Y_{1}$, $\phi_{1}$, $w$, $W_{1}$ and $b_{1}$ are the output of the
first layer, the activation function, inputs, weight matrix and the biased
term. Here, $Y_{1}$, $w$ and $b_{1}$ are row vectors. Similarly,

$\displaystyle Y_{i+1}=\phi_{i+1}(Y_{i}W_{i+1}+b_{i+1})$

and

$\displaystyle Y_{I}=Y_{I-1}W_{I}+b_{I},$

where $I$ is the number of layers and $i=1,2,\cdots,I-1$. Then, we can obtain
$\\{u_{H,m,1}^{n,\prime}\\}_{n}$ by reshaping $Y_{I}$.

Next, we would like to prove two theoretical results. Under the setting in
Section 3, we first define the map from the input to the output of our neural
network as

$\displaystyle h:W$ $\displaystyle\rightarrow\mathbb{R}^{l\times(N-1)},$
$\displaystyle w$
$\displaystyle\mapsto(u_{H,w,1}^{2,\prime,1},u_{H,w,1}^{2,\prime,2},\dots,u_{H,w,1}^{2,\prime,l},u_{H,w,1}^{3,\prime,1},u_{H,w,1}^{3,\prime,2},\dots,u_{H,w,1}^{3,\prime,l},\dots,u_{H,w,1}^{N,\prime,1},u_{H,w,1}^{N,\prime,2},\dots,u_{H,w,1}^{N,\prime,l}),$

where $W$ is a compact subset of $\mathbb{R}^{k}$ with norm $\|\cdot\|_{W}$.
We will show that the neural network is capable of approximating $h$
precisely. We then show that when the loss function (mean squared error) goes
to zero, the output of the neural network converges to the exact solution
(which is computed via the partially explicit scheme). We state the following
well-known theorem [36].

######  Theorem 3.1.

Universal Approximation Theorem. Let $\sigma:\mathbb{R}\rightarrow\mathbb{R}$
be any non-affine continuous function. Also, $\sigma$ is continuously
differentiable at at least one point and with nonzero derivative at that
point. Let $\mathcal{N}_{k,D}^{\sigma}$ be the space of neural networks with
$k$ inputs, $D$ outputs, $\sigma$ as the activation function for every hidden
layers and the identity function as the activation function for the output
layers. Then given any $\epsilon>0$ and any function $\mathcal{H}\in
C\left(W,\mathbb{R}^{D}\right)$, there exists
$\hat{\mathcal{H}}\in\mathcal{N}_{k,D}^{\sigma}$ such that

$\sup_{w\in
W}\|\hat{\mathcal{H}}(w)-\mathcal{H}(w)\|_{l_{\infty}(\mathbb{R{{}^{D}})}}<\epsilon.$

The universal approximation theorem relies on the continuity assumption of the
mapping $\mathcal{H}$. More specifically, to apply the universal approximation
theorem, we now need to show that $h$ is continuous.

######  Lemma 3.1.

We define $\|u\|_{a}^{2}=a(u,u)=(f(u),u)$ and assume $f$ is linear. We know
$V_{H}=V_{H,1}\oplus V_{H,2}$. The two subspaces $V_{H,1}$ and $V_{H,2}$ are
finite dimensional subspaces with trivial intersection. By the strengthened
Cauchy Schwarz inequality [2], there exists a constant $\gamma$ depending on
$V_{H,1}$ and $V_{H,2}$ such that

$\displaystyle 0<\gamma:=\sup_{v_{1}\in V_{H,1},v_{2}\in
V_{H,2}}\cfrac{(v_{1},v_{2})}{\|v_{1}\|\|v_{2}\|}<1.$ (3.3)

If

$\displaystyle\Delta t\leq(1-\gamma)\left(\sup_{v_{2}\in
V_{H,2}}\cfrac{\|v_{2}\|_{a}^{2}}{\|v_{2}\|^{2}}\right)^{-1},$ (3.4)

then $h$ is continuous with respect to $w\in W$.

The proof is given in Appendix A. With the continuity of $h$, we have the
following theorem.

######  Theorem 3.2.

For any $\epsilon>0$, there exists $\hat{h}\in\mathcal{N}_{k,l(N-1)}^{\sigma}$
such that

$\sup_{w\in
W}\|\hat{h}(w)-h(w)\|_{l_{\infty}(\mathbb{R}^{l(N-1)})}<\epsilon/3.$

The theorem follows from the continuity of $h$ and the Universal Approximation
Theorem. Finally, we present the estimation of the generalization error.

######  Theorem 3.3.

Let $W_{T}\subset W$ be the set containing all training samples. Fix
$\epsilon>0$. The functions $h$ and $\hat{h}$ are uniformly continuous since
they are continuous and $W$ is compact. It follows that there exists
$\delta>0$ such that

$\displaystyle\|\hat{h}(w^{\prime})-\hat{h}(w^{\prime\prime})\|_{l_{\infty}}<\epsilon/3,\quad\|h(w^{\prime})-h(w^{\prime\prime})\|_{l_{\infty}}<\epsilon/3,$

for any $w^{\prime},w^{\prime\prime}\in W$ and
$\|w^{\prime}-w^{\prime\prime}\|<\delta$. We assume that for any $w\in W$,
there exists $w_{i}\in W_{T}$ such that $\|w-w_{i}\|_{W}<\delta$. Let
$\hat{h}\in\mathcal{N}_{k,l(N-1)}^{\sigma}$ be the mapping in Theorem 3.2.
Then we have

$\sup_{w\in W}\|\hat{h}(w)-h(w)\|_{l_{\infty}(\mathbb{R}^{l(N-1)})}<\epsilon.$

###### Proof.

For any $w\in W$, we can find $w_{i}\in W_{T}$ such that
$\|w-w_{i}\|_{W}<\delta$, it follows that:

$\displaystyle\|\hat{h}(w)-h(w)\|_{l_{\infty}}$
$\displaystyle=\|\hat{h}(w)-\hat{h}(w_{i})+\hat{h}(w_{i})-h(w_{i})+h(w_{i})-h(w)\|_{l_{\infty}}$
$\displaystyle\leq\|\hat{h}(w)-\hat{h}(w_{i})\|_{l_{\infty}}+\|\hat{h}(w_{i})-h(w_{i})\|_{l_{\infty}}+\|h(w_{i})-h(w)\|_{l_{\infty}}$
$\displaystyle\leq\epsilon/3+\epsilon/3+\epsilon/3=\epsilon.$

As $\epsilon$ and $w$ are arbitrary, we obtain the conclusion. ∎

##  4 Numerical Results

In this section, we present numerical examples. In all cases, the domain
$\Omega$ is $[0,1]^{2}$. The permeability field $\kappa$ we use is shown in
Figure 4.1. Note that there are high contrast channels in the permeability
field. The fine mesh size is $\sqrt{2}/100$ and the coarse mesh size is
$\sqrt{2}/10$.

Figure 4.1: The permeability field $\kappa$.

We define the following notation.

  * •

$u_{H,l}$ is obtained by predicting $u_{H,1}$ and computing $u_{H,2}$ by
Equation (2.8).

  * •

$u_{H,c}$ is obtained by computing $u_{H,1}$ and $u_{H,2}$ by Equation (2.7)
and (2.8).

  * •

$u_{H,f}$ is obtained using the fine grid basis functions and Equation (2.5).

We also need to define the following relative error for better illustration.
Let $\|\cdot\|$ be the $L^{2}$ norm. We know that
$u_{H,l}=u_{H,l,1}+u_{H,l,2}$ and $u_{H,c}=u_{H,c,1}+u_{H,c,2}$.

  * •

$e_{1}=\cfrac{\|u_{H,l}-u_{H,f}\|}{\|u_{H,f}\|}\times 100\%$,
$e_{2}=\cfrac{\|u_{H,c}-u_{H,f}\|}{\|u_{H,f}\|}\times 100\%$.

  * •

$e_{3}=\cfrac{\|u_{H,l,1}-u_{H,c,1}\|}{\|u_{H,c,1}\|}\times 100\%$,
$e_{4}=\cfrac{\|u_{H,l}-u_{H,c}\|}{\|u_{H,c}\|}\times 100\%$.

$e_{1}$ and $e_{2}$ are used to measure the difference of accuracy between our
learning scheme and the scheme in which we compute both $u_{H,1}$ and
$u_{H,2}$. $e_{3}$ and $e_{4}$ are defined to measure how close the learning
output is to their target. The neural network is fully connected and is
activated by the scaled exponential linear unit (SELU). To obtain the solution
of the minimization problem, we use the gradient-based optimizer Adam [37].
The loss function we use is the mean square error (MSE). Our machine learning
algorithm is implemented using the open-source library Keras [8] with the
Tensorflow [1] backend. The neural network we use is a simple neural network
with $5$ fully connected layers. We remark that with this neural network, we
can obtain high accuracy.

In the first two examples, we consider the following linear parabolic partial
differential equations

$\displaystyle u_{t}-\nabla\cdot(\kappa\nabla u)$
$\displaystyle=g_{w}(x,t),\;(x,t)\in\Omega\times[0,T],$ $\displaystyle\nabla
u(x)\cdot\vec{n}$ $\displaystyle=0,\;x\in\partial\Omega,$ $\displaystyle
u(x,0)$ $\displaystyle=u_{0},$

where $\vec{n}$ is the outward normal vector and $u_{0}$ is the initial
condition which is presented it in Figure 4.2. In our first example, the
source term is defined as follows,

$\displaystyle g_{w}(x,t)=\begin{cases}100(w_{1}\sin(\cfrac{2\pi
t}{T})+w_{2}\sin(\cfrac{5.2\pi t}{T}))&\text{ for }x\in[0.2,0.3]^{2}\\\
100(w_{3}\sin(\cfrac{2.4\pi t}{T})+w_{4}\sin(\cfrac{4\pi t}{T}))&\text{ for
}x\in[0.8,0.9]^{2}\\\ 0&\text{ otherwise, }\end{cases}$

where $(w_{1},w_{2},w_{3},w_{4})$ are parameters. In this case, we let $1\leq
w_{i}\leq 10$ for $i=1,2,3,4$. We show a simple sketch of $g_{w}(x,t)$ in
Figure 4.2 to illustrate. This sketch only shows the position of the coarse
blocks where the source term is nonzero. We choose this source, which has some
similarities to the wells in petroleum engineering. The final time is
$T=0.01$, there are $N=100$ time steps and the time step size is $\Delta
t=T/N=10^{-4}$.

Figure 4.2: Left: Initial condition $u_{0}$. Right: Source term $g(x,t)$. The
square boxes indicate the positions of the wells.

The training space contains $1000$ samples and the testing space contains
$500$ samples. Every sample’s parameter $w$ is drawn from the uniform
distribution over $[1,10]$. The number of POD basis vectors is $l=15$. We
choose this $l$ based on the eigenvalue behavior in the process of singular
value decomposition. We present $e_{1}$ and $e_{2}$ defined before in Figure
4.3. The error is calculated at each time step as the average over all samples
in the testing space. From the plot in Figure 4.3, we find that the curves for
$e_{1}$ and $e_{2}$ coincide, which means our learning algorithms can obtain
similar accuracy as the scheme in which we compute $u_{H}$ by Equation (2.7)
and (2.8). The average of $e_{3}$ and $e_{4}$ over time are $0.151\%$ and
$0.155\%$, respectively, which means the output of our learning scheme nearly
resembles the corresponding target (the solution we obtain by computing both
equations).

Figure 4.3: $e_{1}$ and $e_{2}$ of the linear case with $4$ parameters.

In the next case, we test our algorithms using a more complicated source term
with more parameters. The source term we use is,

$\displaystyle g_{w}(x,t)=\begin{cases}100(w_{1}\sin(\cfrac{\pi
t}{T})+w_{2}\cos(\cfrac{3.2\pi t}{T}))&\text{ for }x\in[0.2,0.3]^{2}\\\
100(w_{3}\sin(\cfrac{2.2\pi t}{T})+w_{4}\sin(\cfrac{1.6\pi t}{T}))&\text{ for
}x\in[0.8,0.9]^{2}\\\ 100(w_{5}\cos(\cfrac{3\pi t}{T})+w_{6}\cos(\cfrac{4.6\pi
t}{T}))&\text{ for }x\in[0.2,0.3]\times[0.8,0.9]\\\
100(w_{7}\cos(\cfrac{1.4\pi t}{T})+w_{8}\sin(\cfrac{5\pi t}{T}))&\text{ for
}x\in[0.8,0.9]\times[0.2,0.3]\\\ 100(w_{9}\sin(\cfrac{2.8\pi
t}{T})+w_{10}\sin(\cfrac{4\pi t}{T}))&\text{ for }x\in[0.5,0.6]^{2}\\\
0&\text{ otherwise, }\end{cases}$

where $(w_{1},w_{2},w_{3},\cdots,w_{10})$ are parameters. In this case, the
parameters $w_{i}$ of every sample are drawn from the uniform distribution
over $[1,20]$. We present an easy sketch of $g_{w}(x,t)$ in Figure 4.4 to
illustrate the positions of the sources. The initial condition $u_{0}$ is also
shown in Figure 4.4. The final time $T$, the number of time steps $N$ and the
time step size are $0.01$, $100$ and $10^{-4}$, respectively.

Figure 4.4: Left: Initial condition $u_{0}$. Right: Source term $g_{w}(x,t)$.
The square boxes in the plot indicate the positions of the wells.

The model is trained using 2000 samples and we test it with 1000 samples. The
number of POD basis vectors is $l=25$. In Figure 4.5, the $e_{1}$ and $e_{2}$
are presented. By observing the error plot, we find that $e_{1}$ and $e_{2}$
are close to each other, which implies that our learning algorithms can obtain
similar accuracy as the scheme where we calculate both $u_{H,1}$ and $u_{H,2}$
using Equation (2.7) and (2.8). The average of $e_{3}$ and $e_{4}$ over time
are $0.673\%$ and $0.677\%$, respectively. Such small $e_{3}$ and $e_{4}$
indicate that the prediction is accurate.

Figure 4.5: $e_{1}$ and $e_{2}$ of the linear case with $10$ parameters.

In the next case, we consider the following equation

$\displaystyle u_{t}-\nabla\cdot(\kappa e^{u}\nabla u)$
$\displaystyle=g(x,t),\;(x,t)\in\Omega\times[0,T],$ $\displaystyle\nabla
u(x)\cdot\vec{n}=0,\;x\in\partial\Omega,$ $\displaystyle u(x,0)=u_{0},$

where $\vec{n}$ is the outward normal vector. This equation is nonlinear and
its computation cost is much larger than the linear case because we have to
assemble the stiffness matrix at every time step. We use the same initial
condition $u_{0}$ and source term $g(x,t)$ as the second example. The final
time is $T=0.001$, the number of time steps is $N=40$ and the time step size
is $\Delta t=T/N=2.5\times 10^{-5}$. We train our neural network using $1000$
samples and the network is tested using $500$ samples. In Figure 4.6, we
present $e_{1}$ and $e_{2}$ defined before. From the first plot, we find that
the curves for $e_{1}$ and $e_{2}$ nearly coincide, which implies that our
machine learning algorithms can have the same accuracy as we compute the
solution by Equations (2.7) and (2.8). The average of $e_{3}$ and $e_{4}$ over
time are $0.460\%$ and $0.450\%$, respectively, which means our neural network
can approximate the target well and achieve a fantastic accuracy.

Figure 4.6: $e_{1}$ and $e_{2}$ of the nonlinear case with 10 parameters.

##  5 Conclusion

In this paper, we construct neural networks based on the partially explicit
splitting scheme. We predict $u_{1}$, the implicit part of the solution, while
solving for $u_{2}$, the explicit part of the solution. The output dimension
resulting from $u_{1}$ is large and we apply the POD to reduce the dimension.
The POD transformation greatly reduces the output dimension and at the same
time the number of parameters in the neural networks. The POD error is small
compared to the overall error. We discuss three numerical examples. We present
a theoretical justification. We conclude that our machine learning algorithms
can obtain similar accuracy as the partially explicit scheme.

##  Appendix A Proof of Lemma 3.1

In this section, we prove the Lemma 3.1, i.e., the function

$\displaystyle h:W$ $\displaystyle\rightarrow\mathbb{R}^{l\times(N-1)},$
$\displaystyle w$
$\displaystyle\mapsto(u_{H,w,1}^{2,\prime,1},u_{H,w,1}^{2,\prime,2},\dots,u_{H,w,1}^{2,\prime,l},u_{H,w,1}^{3,\prime,1},u_{H,w,1}^{3,\prime,2},\dots,u_{H,w,1}^{3,\prime,l},\dots,u_{H,w,1}^{N,\prime,1},u_{H,w,1}^{N,\prime,2},\dots,u_{H,w,1}^{N,\prime,l})$

is continuous. We equip $\mathbb{R}^{l\times(N-1)}$ with the $l_{\infty}$
norm. Let us decompose the mapping as $h=\mathcal{P}\circ h_{2}\circ h_{1}$,
where

$\displaystyle
h_{1}:(W,\|\cdot\|_{\infty})\to(L^{2}(\Omega\times[0,T],\mathbb{R}),\|\cdot\|_{L^{2}})\quad
w\mapsto g_{w}(x,t),$ $\displaystyle
h_{2}:(L^{2}(\Omega\times[0,T],\mathbb{R}),\|\cdot\|_{L^{2}})\to(H^{1}(\Omega),\|\cdot\|_{H^{1}}),\quad
g_{w}(x,t)\mapsto\\{u_{H}(x,t^{n})\\}_{n=1}^{N},$

and $\mathcal{P}$ is the POD transformation defined by (3.2). The map
$\mathcal{P}$ is a projection, thus continuous. We assume that $g_{w}$ is
continuous with respect to $w$ and therefore, $h_{1}$ is continuous. To show
that $h$ is continuous, we only need to prove $h_{2}$ is continuous. The proof
is given in the next lemma.

######  Lemma A.1.

Let us assume $f$ is linear, and we define a bilinear form $a(u,v):=(f(u),v)$
and its co-responding norm as $\|u\|_{a}:=\sqrt{a(u,u)}$. If

$\displaystyle\Delta t\leq(1-\gamma)\left(\sup_{v_{2}\in
V_{H,2}}\cfrac{\|v_{2}\|_{a}^{2}}{\|v_{2}\|^{2}}\right)^{-1},$ (A.1)

the function $h_{2}$ is continuous with respect to the source term
$g_{w}(x,t)$.

###### Proof.

Case I. The solution $u_{i}(x,\Delta t)$ for $i=1,2$ is continuous with
respect to the source term $g_{w}(x,t)$.
The initial condition $u(x,0)$ is given and the solution at time step
$1\cdot\Delta t$ is calculated via the fully implicit scheme. Thus, we want to
show the continuity from the source term to the solution $u_{i}^{1}$ of the
following scheme

$\displaystyle(\frac{u^{1}_{H}-u^{0}_{H}}{\Delta
t},v)+(f(u^{1}_{H}),v)=(g_{w}(x,t^{1}),v)\quad\forall v\in V_{H},$

where $V_{H}=V_{H,1}\oplus V_{H,2}$.

In the following proof, we abbreviate the subscript $H$ for simplicity.

Let $w_{1}$ and $w_{2}$ be different, and let us denote the corresponding
sources as $g_{w_{1}}$ and $g_{w_{2}}$. The discretized solutions
$u^{1}_{w_{1}}$ and $u^{1}_{w_{2}}$ evolved for one time step then satisfy:

$\displaystyle(\frac{u^{1}_{w_{1}}-u^{0}}{\Delta
t},v)+(f(u^{1}_{w_{1}}),v)=(g_{w_{1}}(x,t^{1}),v)\quad\forall v\in V_{H},$
$\displaystyle(\frac{u^{1}_{w_{2}}-u^{0}}{\Delta
t},v)+(f(u^{1}_{w_{2}}),v)=(g_{w_{2}}(x,t^{1}),v)\quad\forall v\in V_{H}.$

It follows immediately that,

$\displaystyle(\frac{u^{1}_{w_{1}}-u^{1}_{w_{2}}}{\Delta
t},v)+(f(u^{1}_{w_{1}})-f(u^{1}_{w_{2}}),v)=(g_{w_{1}}(x,t^{1})-g_{w_{2}}(x,t^{1}),v)\quad\forall
v\in V_{H}.$

Test with $v=u^{1}_{w_{1}}-u^{1}_{w_{2}}$, the linearity of $f$ then yields,

$\displaystyle\begin{split}\frac{1}{\Delta
t}\|u^{1}_{w_{1}}-u^{1}_{w_{2}}\|^{2}+(f(u^{1}_{w_{1}}-u^{1}_{w_{2}}),u^{1}_{w_{1}}-u^{1}_{w_{2}})&\leq\|g_{w_{1}}(x,t^{1})-g_{w_{2}}(x,t^{1})\|\|u^{1}_{w_{1}}-u^{1}_{w_{2}}\|\\\
&\leq\frac{\Delta
t}{2}\|g_{w_{1}}(x,t^{1})-g_{w_{2}}(x,t^{1})\|^{2}+\frac{1}{2\Delta
t}\|u^{1}_{w_{1}}-u^{1}_{w_{2}}\|^{2}\end{split}$

After rearrangement, we have

$\displaystyle\begin{split}\frac{1}{2\Delta
t}\|u^{1}_{w_{1}}-u^{1}_{w_{2}}\|^{2}+\|u^{1}_{w_{1}}-u^{1}_{w_{2}}\|_{a}^{2}\leq\frac{\Delta
t}{2}\|g_{w_{1}}(x,t^{1})-g_{w_{2}}(x,t^{1})\|^{2}.\end{split}$

Apply the strengthened Cauchy Schwartz inequality (3.3), let
$u_{w_{m}}^{n}=u_{1,w_{m}}^{n}+u_{2,w_{m}}^{n}$ with $u_{1,w_{m}}^{n}\in
V_{H,1}$ and $u_{2,w_{m}}^{n}\in V_{H,2}$, we have,

$\displaystyle\|u_{w_{1}}^{n}-u_{w_{2}}^{n}\|^{2}$
$\displaystyle=\|u_{1,w_{1}}^{n}-u_{1,w_{2}}^{n}\|^{2}+\|u_{2,w_{1}}^{n}-u_{2,w_{2}}^{n}\|^{2}+2(u_{1,w_{1}}^{n}-u_{1,w_{2}}^{n},u_{2,w_{1}}^{n}-u_{2,w_{2}}^{n})$
$\displaystyle\geq\|u_{1,w_{1}}^{n}-u_{1,w_{2}}^{n}\|^{2}+\|u_{2,w_{1}}^{n}-u_{2,w_{2}}^{n}\|^{2}-2\gamma(u_{1,w_{1}}^{n}-u_{1,w_{2}}^{n},u_{2,w_{1}}^{n}-u_{2,w_{2}}^{n})$
$\displaystyle\geq(1-\gamma)(\|u_{1,w_{1}}^{n}-u_{1,w_{2}}^{n}\|^{2}+\|u_{2,w_{1}}^{n}-u_{2,w_{2}}^{n}\|^{2}),$

for all $n>0$.

Thus,

$\displaystyle\frac{1-\gamma}{2\Delta
t}\sum_{i}\|u^{1}_{i,w_{1}}-u^{1}_{i,w_{2}}\|^{2}+\|u^{1}_{w_{1}}-u^{1}_{w_{2}}\|_{a}^{2}\leq\frac{\Delta
t}{2}\|g_{w_{1}}(x,t^{1})-g_{w_{2}}(x,t^{1})\|^{2}.$ (A.2)

Case II. The solution $u_{1}(x,n\Delta t)$ is continuous with respect to the
source term $g_{w}(x,t)$, for $n=2,3,\dots,N$.

In the partially explicit scheme, for the source term $g_{w_{1}}$, the
solution $u_{i,w_{1}}^{n+1}$ satisfies

$\displaystyle(\frac{u_{1,w_{1}}^{n+1}-u_{1,w_{1}}^{n}}{\Delta
t}+\frac{u_{2,w_{1}}^{n}-u_{2,w_{1}}^{n-1}}{\Delta
t},v_{1})+(f(u_{1,w_{1}}^{n+1}+u_{2,w_{1}}^{n}),v_{1})=(g_{w_{1}}(x,t^{n}),v_{1})\quad\forall
v_{1}\in V_{H,1},$ (A.3)
$\displaystyle(\frac{u_{2,w_{1}}^{n+1}-u_{2,w_{1}}^{n}}{\Delta
t}+\frac{u_{1,w_{1}}^{n}-u_{1,w_{1}}^{n-1}}{\Delta
t},v_{2})+(f(u_{1,w_{1}}^{n+1}+u_{2,w_{1}}^{n}),v_{2})=(g_{w_{1}}(x,t^{n}),v_{2})\quad\forall
v_{2}\in V_{H,2}.$ (A.4)

Similarly, for $g_{w_{2}}$, the solution $u_{i,w_{2}}^{n+1}$ satisfies

$\displaystyle(\frac{u_{1,w_{2}}^{n+1}-u_{1,w_{2}}^{n}}{\Delta
t}+\frac{u_{2,w_{2}}^{n}-u_{2,w_{2}}^{n-1}}{\Delta
t},v_{1})+(f(u_{1,w_{2}}^{n+1}+u_{2,w_{2}}^{n}),v_{1})=(g_{w_{2}}(x,t^{n}),v_{1})\quad\forall
v_{1}\in V_{H,1},$ (A.5)
$\displaystyle(\frac{u_{2,w_{2}}^{n+1}-u_{2,w_{2}}^{n}}{\Delta
t}+\frac{u_{1,w_{2}}^{n}-u_{1,w_{2}}^{n-1}}{\Delta
t},v_{2})+(f(u_{1,w_{2}}^{n+1}+u_{2,w_{2}}^{n}),v_{2})=(g_{w_{2}}(x,t^{n}),v_{2})\quad\forall
v_{2}\in V_{H,2}.$ (A.6)

Let us denote $r^{n}_{i}=u_{i,w_{1}}^{n}-u_{i,w_{2}}^{n}$ for $i=1,2$, and
subtract (A.5) from (A.3) and (A.6) from (A.4), it follows that,

$\displaystyle(\frac{r_{1}^{n+1}-r_{1}^{n}}{\Delta
t}+\frac{r_{2}^{n}-r_{2}^{n-1}}{\Delta
t},v_{1})+(f(u_{1,w_{1}}^{n+1}+u_{2,w_{1}}^{n})-f(u_{1,w_{2}}^{n+1}+u_{2,w_{2}}^{n}),v_{1})$
$\displaystyle=(g_{w_{1}}(x,t^{n})-g_{w_{2}}(x,t^{n}),v_{1})\quad\forall
v_{1}\in V_{H,1},$ $\displaystyle(\frac{r_{2}^{n+1}-r_{2}^{n}}{\Delta
t}+\frac{r_{1}^{n}-r_{1}^{n-1}}{\Delta
t},v_{2})+(f(u_{1,w_{1}}^{n+1}+u_{2,w_{1}}^{n})-f(u_{1,w_{2}}^{n+1}+u_{2,w_{2}}^{n}),v_{2})$
$\displaystyle=(g_{w_{1}}(x,t^{n})-g_{w_{2}}(x,t^{n}),v_{2})\quad\forall
v_{2}\in V_{H,2}.$

The linearity of $f$ then yields,

$\displaystyle(\frac{r_{1}^{n+1}-r_{1}^{n}}{\Delta
t}+\frac{r_{2}^{n}-r_{2}^{n-1}}{\Delta
t},v_{1})+(f(r_{1}^{n+1}+r_{2}^{n}),v_{1})=(g_{w_{1}}(x,t^{n})-g_{w_{2}}(x,t^{n}),v_{1})\quad\forall
v_{1}\in V_{H,1},$ $\displaystyle(\frac{r_{2}^{n+1}-r_{2}^{n}}{\Delta
t}+\frac{r_{1}^{n}-r_{1}^{n-1}}{\Delta
t},v_{2})+(f(r_{1}^{n+1}+r_{2}^{n}),v_{2})=(g_{w_{1}}(x,t^{n})-g_{w_{2}}(x,t^{n}),v_{2})\quad\forall
v_{2}\in V_{H,2}.$

Let $v_{1}=r_{1}^{n+1}-r_{1}^{n}$ and $v_{2}=r_{2}^{n+1}-r_{2}^{n}$. We get

$\displaystyle(\frac{r_{1}^{n+1}-r_{1}^{n}}{\Delta
t}+\frac{r_{2}^{n}-r_{2}^{n-1}}{\Delta
t},r_{1}^{n+1}-r_{1}^{n})+(f(r_{1}^{n+1}+r_{2}^{n}),r_{1}^{n+1}-r_{1}^{n})=(g_{w_{1}}(x,t^{n})-g_{w_{2}}(x,t^{n}),r_{1}^{n+1}-r_{1}^{n}),$
$\displaystyle(\frac{r_{2}^{n+1}-r_{2}^{n}}{\Delta
t}+\frac{r_{1}^{n}-r_{1}^{n-1}}{\Delta
t},r_{2}^{n+1}-r_{2}^{n})+(f(r_{1}^{n+1}+r_{2}^{n}),r_{2}^{n+1}-r_{2}^{n})=(g_{w_{1}}(x,t^{n})-g_{w_{2}}(x,t^{n}),r_{2}^{n+1}-r_{2}^{n}).$

Sum up these two equations, we have

$\displaystyle\frac{1}{\Delta
t}\sum_{i}\|r_{i}^{n+1}-r_{i}^{n}\|^{2}+\frac{1}{\Delta t}\sum_{i\neq
j}(r_{i}^{n+1}-r_{i}^{n},r_{j}^{n}-r_{j}^{n-1})+(f(r_{1}^{n+1}+r_{2}^{n}),r^{n+1}-r^{n})$
$\displaystyle=(g_{w_{1}}(x,t^{n})-g_{w_{2}}(x,t^{n}),r^{n+1}-r^{n}),$ (A.7)

where $r^{n}=r^{n}_{1}+r^{n}_{2}$. The second term of (A.7) on the left hand
side can be further estimated as,

$\displaystyle\frac{1}{\Delta t}\sum_{i\neq
j}|(r_{i}^{n+1}-r_{i}^{n},r_{j}^{n}-r_{j}^{n-1})|$
$\displaystyle\leq\frac{\gamma}{\Delta t}\sum_{i\neq
j}\|r_{i}^{n+1}-r_{i}^{n}\|\|r_{j}^{n}-r_{j}^{n-1}\|$
$\displaystyle\leq\frac{\gamma}{2\Delta t}\sum_{i\neq
j}(\|r_{i}^{n+1}-r_{i}^{n}\|^{2}+\|r_{j}^{n}-r_{j}^{n-1}\|^{2}),$

where we apply the strenthened Cauchy Schwartz inequality. It follows that,

$\displaystyle\frac{1}{\Delta
t}\sum_{i}\|r_{i}^{n+1}-r_{i}^{n}\|^{2}+\frac{1}{\Delta t}\sum_{i\neq
j}(r_{i}^{n+1}-r_{i}^{n},r_{j}^{n}-r_{j}^{n-1})$ $\displaystyle\geq$
$\displaystyle\frac{2-\gamma}{2\Delta
t}\sum_{i}\|r_{i}^{n+1}-r_{i}^{n}\|^{2}-\frac{\gamma}{2\Delta
t}\sum_{i}\|r_{i}^{n}-r_{i}^{n-1}\|^{2}$ $\displaystyle=$
$\displaystyle(\frac{\gamma}{2\Delta t}+\frac{1-\gamma}{\Delta
t})\sum_{i}\|r_{i}^{n+1}-r_{i}^{n}\|^{2}-\frac{\gamma}{2\Delta
t}\sum_{i}\|r_{i}^{n}-r_{i}^{n-1}\|^{2}.$

Note that
$(f(r_{1}^{n+1}+r_{2}^{n}),r^{n+1}-r^{n})=a(r_{1}^{n+1}+r_{2}^{n},r^{n+1}-r^{n})$,
we then can estimate the third term of (A.7),

$\displaystyle(f(r_{1}^{n+1}+r_{2}^{n}),r^{n+1}-r^{n})$ $\displaystyle=$
$\displaystyle
a(r^{n+1},r^{n+1}-r^{n})-a(r_{2}^{n+1}-r_{2}^{n},r^{n+1}-r^{n})$
$\displaystyle=$
$\displaystyle\frac{1}{2}(\|r^{n+1}\|_{a}^{2}+\|r^{n+1}-r^{n}\|^{2}_{a}-\|r^{n}\|_{a}^{2})-a(r_{2}^{n+1}-r_{2}^{n},r^{n+1}-r^{n})$
$\displaystyle\geq$
$\displaystyle\frac{1}{2}(\|r^{n+1}\|_{a}^{2}+\|r^{n+1}-r^{n}\|^{2}_{a}-\|r^{n}\|_{a}^{2})-\frac{1}{2}\|r_{2}^{n+1}-r_{2}^{n}\|_{a}^{2}-\frac{1}{2}\|r^{n+1}-r^{n}\|_{a}^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{2}(\|r^{n+1}\|_{a}^{2}-\|r^{n}\|_{a}^{2})-\frac{1}{2}\|r_{2}^{n+1}-r_{2}^{n}\|_{a}^{2}.$

Now (A.7) can be estimated as:

$\displaystyle(\frac{\gamma}{2\Delta t}+\frac{1-\gamma}{\Delta
t})\sum_{i}\|r_{i}^{n+1}-r_{i}^{n}\|^{2}-\frac{\gamma}{2\Delta
t}\sum_{i}\|r_{i}^{n}-r_{i}^{n-1}\|^{2}+\frac{1}{2}(\|r^{n+1}\|_{a}^{2}-\|r^{n}\|_{a}^{2})-\frac{1}{2}\|r_{2}^{n+1}-r_{2}^{n}\|_{a}^{2}$
$\displaystyle\leq(g_{w_{1}}(x,t^{n})-g_{w_{2}}(x,t^{n}),r^{n+1}-r^{n}).$

After rearranging, we have

$\displaystyle\frac{\gamma}{2\Delta
t}\sum_{i}\|r_{i}^{n+1}-r_{i}^{n}\|^{2}+\frac{1}{2}\|r^{n+1}\|_{a}^{2}+\frac{1-\gamma}{\Delta
t}\sum_{i}\|r_{i}^{n+1}-r_{i}^{n}\|^{2}-\frac{1}{2}\|r_{2}^{n+1}-r_{2}^{n}\|_{a}^{2}$
$\displaystyle\leq$ $\displaystyle\frac{\gamma}{2\Delta
t}\sum_{i}\|r_{i}^{n}-r_{i}^{n-1}\|^{2}+\frac{1}{2}\|r^{n}\|_{a}^{2}+(g_{w_{1}}(x,t^{n})-g_{w_{2}}(x,t^{n}),r^{n+1}-r^{n})$
$\displaystyle\leq$ $\displaystyle\frac{\gamma}{2\Delta
t}\sum_{i}\|r_{i}^{n}-r_{i}^{n-1}\|^{2}+\frac{1}{2}\|r^{n}\|_{a}^{2}+\gamma\|g_{w_{1}}(x,t^{n})-g_{w_{2}}(x,t^{n})\|\|r^{n+1}-r^{n}\|$
$\displaystyle\leq$ $\displaystyle\frac{\gamma}{2\Delta
t}\sum_{i}\|r_{i}^{n}-r_{i}^{n-1}\|^{2}+\frac{1}{2}\|r^{n}\|_{a}^{2}+\frac{\gamma^{2}\Delta
t}{1-\gamma}\|g_{w_{1}}(x,t^{n})-g_{w_{2}}(x,t^{n})\|^{2}+\frac{1-\gamma}{4\Delta
t}\|r^{n+1}-r^{n}\|^{2}$ $\displaystyle\leq$
$\displaystyle\frac{\gamma}{2\Delta
t}\sum_{i}\|r_{i}^{n}-r_{i}^{n-1}\|^{2}+\frac{1}{2}\|r^{n}\|_{a}^{2}+\frac{\gamma^{2}\Delta
t}{1-\gamma}\|g_{w_{1}}(x,t^{n})-g_{w_{2}}(x,t^{n})\|^{2}+\frac{1-\gamma}{2\Delta
t}\sum_{i}\|r_{i}^{n+1}-r_{i}^{n}\|^{2}.$

We assume the CFL condition (A.1), specifically,

$\displaystyle(\frac{1-\gamma}{\Delta t}-\frac{1-\gamma}{2\Delta
t})\sum_{i}\|r_{i}^{n+1}-r_{i}^{n}\|^{2}-\frac{1}{2}\|r_{2}^{n+1}-r_{2}^{n}\|_{a}^{2}\geq
0.$

This yields,

$\displaystyle\begin{split}&\frac{\gamma}{2\Delta
t}\sum_{i}\|r_{i}^{n+1}-r_{i}^{n}\|^{2}+\frac{1}{2}\|r^{n+1}\|_{a}^{2}\\\
\leq&\frac{\gamma}{2\Delta
t}\sum_{i}\|r_{i}^{n}-r_{i}^{n-1}\|^{2}+\frac{1}{2}\|r^{n}\|_{a}^{2}+\frac{\gamma^{2}\Delta
t}{1-\gamma}\|g_{w_{1}}(x,t^{n})-g_{w_{2}}(x,t^{n})\|^{2}\\\
\leq&\frac{\gamma}{2\Delta
t}\sum_{i}\|u_{i,w_{1}}^{1}-u_{i,w_{2}}^{1}-(u_{i,w_{1}}^{0}-u_{i,w_{2}}^{0})\|^{2}+\frac{1}{2}\|u_{w_{1}}^{1}-u_{w_{2}}^{1}\|_{a}^{2}+\frac{\gamma^{2}\Delta
t}{1-\gamma}\sum_{k=1}^{n}\|g_{w_{1}}(x,t^{k})-g_{w_{2}}(x,t^{k})\|^{2}\\\
=&\frac{\gamma}{2\Delta
t}\sum_{i}\|u_{i,w_{1}}^{1}-u_{i,w_{2}}^{1}\|^{2}+\frac{1}{2}\|u_{w_{1}}^{1}-u_{w_{2}}^{1}\|_{a}^{2}+\frac{\gamma^{2}\Delta
t}{1-\gamma}\sum_{k=1}^{n}\|g_{w_{1}}(x,t^{k})-g_{w_{2}}(x,t^{k})\|^{2}\\\
\leq&\frac{\gamma\Delta
t}{2(1-\gamma)}\|g_{w_{1}}(x,t^{1})-g_{w_{2}}(x,t^{1})\|^{2}+\frac{\gamma^{2}\Delta
t}{1-\gamma}\sum_{k=1}^{n}\|g_{w_{1}}(x,t^{k})-g_{w_{2}}(x,t^{k})\|^{2},\end{split}$
(A.8)

where we use (A.2) in the last inequality. Then,

$\displaystyle\frac{\gamma}{2\Delta
t}\sum_{i}\|r_{i}^{n+1}\|^{2}\leq\frac{\gamma}{\Delta
t}\sum_{i}(\|r_{i}^{n+1}-r_{i}^{n}\|^{2}+\|r_{i}^{n}\|^{2})$
$\displaystyle\leq$ $\displaystyle\frac{\gamma\Delta
t}{1-\gamma}\|g_{w_{1}}(x,t^{1})-g_{w_{2}}(x,t^{1})\|^{2}+\frac{2\gamma^{2}\Delta
t}{1-\gamma}\sum_{k=1}^{n}\|g_{w_{1}}(x,t^{k})-g_{w_{2}}(x,t^{k})\|^{2}+\frac{\gamma}{\Delta
t}\sum_{i}\|r_{i}^{n}\|^{2}$ $\displaystyle\leq$
$\displaystyle\frac{n\gamma\Delta
t}{1-\gamma}\|g_{w_{1}}(x,t^{1})-g_{w_{2}}(x,t^{1})\|^{2}+\frac{2\gamma^{2}\Delta
t}{1-\gamma}\sum_{ind=1}^{n}\sum_{k=1}^{ind}\|g_{w_{1}}(x,t^{k})-g_{w_{2}}(x,t^{k})\|^{2}+\frac{\gamma}{\Delta
t}\sum_{i}\|r_{i}^{1}\|^{2}$ $\displaystyle\leq$
$\displaystyle\frac{(n+1)\gamma\Delta
t}{1-\gamma}\|g_{w_{1}}(x,t^{1})-g_{w_{2}}(x,t^{1})\|^{2}+\frac{2\gamma^{2}\Delta
t}{1-\gamma}\sum_{ind=1}^{n}\sum_{k=1}^{ind}\|g_{w_{1}}(x,t^{k})-g_{w_{2}}(x,t^{k})\|^{2},$

where we apply (A.2) in the last inequality.

Therefore,

$\displaystyle\sum_{i}\|r_{i}^{n+1}\|^{2}\leq$
$\displaystyle\frac{2(n+1)(\Delta
t)^{2}}{1-\gamma}\|g_{w_{1}}(x,t^{1})-g_{w_{2}}(x,t^{1})\|^{2}+\frac{4\gamma(\Delta
t)^{2}}{1-\gamma}\sum_{ind=1}^{n}\sum_{k=1}^{ind}\|g_{w_{1}}(x,t^{k})-g_{w_{2}}(x,t^{k})\|^{2}$
$\displaystyle\leq$ $\displaystyle\frac{2T(\Delta
t)}{1-\gamma}\|g_{w_{1}}(x,t^{1})-g_{w_{2}}(x,t^{1})\|^{2}+\frac{4\gamma
T^{2}}{1-\gamma}\max_{1\leq k\leq
N-1}\|g_{w_{1}}(x,t^{k})-g_{w_{2}}(x,t^{k})\|^{2},$

where $T=N\Delta t$ is the final time. Also, from (A.8), we have

$\displaystyle\|r^{n+1}\|_{a}^{2}\leq$ $\displaystyle\frac{\gamma\Delta
t}{1-\gamma}\|g_{w_{1}}(x,t^{1})-g_{w_{2}}(x,t^{1})\|^{2}+\frac{2\gamma^{2}\Delta
t}{1-\gamma}\sum_{k=1}^{n}\|g_{w_{1}}(x,t^{k})-g_{w_{2}}(x,t^{k})\|^{2}$
$\displaystyle\leq$ $\displaystyle\frac{\gamma\Delta
t}{1-\gamma}\|g_{w_{1}}(x,t^{1})-g_{w_{2}}(x,t^{1})\|^{2}+\frac{2\gamma^{2}T}{1-\gamma}\max_{1\leq
k\leq N-1}\|g_{w_{1}}(x,t^{k})-g_{w_{2}}(x,t^{k})\|^{2}.$

With these two inequalities, we can obtain the conclusion.

∎

## References

  * [1] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  Software available from tensorflow.org.
  * [2] J. Aldaz.  Strengthened Cauchy-Schwarz and Hölder inequalities.  arXiv preprint arXiv:1302.2254, 2013.
  * [3] U. M. Ascher, S. J. Ruuth, and R. J. Spiteri.  Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations.  Applied Numerical Mathematics, 25(2-3):151–167, 1997.
  * [4] I. Babuška and J. M. Melenk.  The partition of unity method.  Int. J. Numer. Meth. Engrg., 40:727–758, 1997.
  * [5] D. L. Brown, Y. Efendiev, and V. H. Hoang.  An efficient hierarchical multiscale finite element method for stokes equations in slowly varying media.  Multiscale Modeling & Simulation, 11(1):30–58, 2013.
  * [6] A. Chatterjee.  An introduction to the proper orthogonal decomposition.  Current Science, 78(7):808–817, 2000.
  * [7] B. Chetverushkin, E. Chung, Y. Efendiev, S.-M. Pun, and Z. Zhang.  Computational multiscale methods for quasi-gas dynamic equations.  Journal of Computational Physics, 440:110352, 2021.
  * [8] F. Chollet et al.  Keras, 2015.
  * [9] E. Chung, Y. Efendiev, W. T. Leung, S.-M. Pun, and Z. Zhang.  Multi-agent reinforcement learning accelerated mcmc on multiscale inversion problem.  arXiv preprint arXiv:2011.08954, 2020.
  * [10] E. Chung, Y. Efendiev, S.-M. Pun, and Z. Zhang.  Computational multiscale methods for parabolic wave approximations in heterogeneous media.  arXiv preprint arXiv:2104.02283, 2021.
  * [11] E. Chung, W. T. Leung, S.-M. Pun, and Z. Zhang.  A multi-stage deep learning based algorithm for multiscale model reduction.  Journal of Computational and Applied Mathematics, 394:113506, 2021\.
  * [12] E. T. Chung, Y. Efendiev, and T. Hou.  Adaptive multiscale model reduction with generalized multiscale finite element methods.  Journal of Computational Physics, 320:69–95, 2016.
  * [13] E. T. Chung, Y. Efendiev, and C. Lee.  Mixed generalized multiscale finite element methods and applications.  SIAM Multiscale Model. Simul., 13:338–366, 2014.
  * [14] E. T. Chung, Y. Efendiev, and W. T. Leung.  Generalized multiscale finite element methods for wave propagation in heterogeneous media.  Multiscale Modeling & Simulation, 12(4):1691–1721, 2014.
  * [15] E. T. Chung, Y. Efendiev, and W. T. Leung.  Constraint energy minimizing generalized multiscale finite element method.  Computer Methods in Applied Mechanics and Engineering, 339:298–319, 2018.
  * [16] E. T. Chung, Y. Efendiev, and W. T. Leung.  Constraint energy minimizing generalized multiscale finite element method in the mixed formulation.  Computational Geosciences, 22(3):677–693, 2018.
  * [17] E. T. Chung, Y. Efendiev, and W. T. Leung.  Fast online generalized multiscale finite element method using constraint energy minimization.  Journal of Computational Physics, 355:450–463, 2018.
  * [18] E. T. Chung, Y. Efendiev, W. T. Leung, and W. Li.  Contrast-independent partially explicit time discretizations for nonlinear multiscale problems, 2021.
  * [19] E. T. Chung, Y. Efendiev, W. T. Leung, M. Vasilyeva, and Y. Wang.  Non-local multi-continua upscaling for flows in heterogeneous fractured media.  Journal of Computational Physics, 372:22–34, 2018.
  * [20] W. E and B. Engquist.  Heterogeneous multiscale methods.  Comm. Math. Sci., 1(1):87–132, 2003.
  * [21] Y. Efendiev, J. Galvis, and T. Hou.  Generalized multiscale finite element methods (GMsFEM).  Journal of Computational Physics, 251:116–135, 2013.
  * [22] Y. Efendiev, J. Galvis, G. Li, and M. Presho.  Generalized multiscale finite element methods. nonlinear elliptic equations.  Communications in Computational Physics, 15(3):733–755, 2014.
  * [23] Y. Efendiev and T. Hou.  Multiscale Finite Element Methods: Theory and Applications, volume 4 of Surveys and Tutorials in the Applied Mathematical Sciences.  Springer, New York, 2009.
  * [24] Y. Efendiev, W. T. Leung, W. Li, S.-M. Pun, and P. N. Vabishchevich.  Nonlocal transport equations in multiscale media. modeling, dememorization, and discretizations, 2022.
  * [25] Y. Efendiev, W. T. Leung, G. Lin, and Z. Zhang.  Hei: hybrid explicit-implicit learning for multiscale problems, 2021.
  * [26] Y. Efendiev and A. Pankov.  Numerical homogenization of monotone elliptic operators.  SIAM J. Multiscale Modeling and Simulation, 2(1):62–79, 2003.
  * [27] Y. Efendiev and A. Pankov.  Homogenization of nonlinear random parabolic operators.  Advances in Differential Equations, 10(11):1235–1260, 2005.
  * [28] W. T. L. Eric Chung, Yalchin Efendiev and P. N. Vabishchevich.  Contrast-independent partially explicit time discretizations for multiscale flow problems.  arXiv:2101.04863.
  * [29] J. Frank, W. Hundsdorfer, and J. Verwer.  On the stability of implicit-explicit linear multistep methods.  Applied Numerical Mathematics, 25(2-3):193–205, 1997.
  * [30] P. Henning, A. Målqvist, and D. Peterseim.  A localized orthogonal decomposition method for semi-linear elliptic problems.  ESAIM: Mathematical Modelling and Numerical Analysis, 48(5):1331–1349, 2014.
  * [31] T. Hou and X. Wu.  A multiscale finite element method for elliptic problems in composite materials and porous media.  J. Comput. Phys., 134:169–189, 1997.
  * [32] T. Y. Hou, D. Huang, K. C. Lam, and P. Zhang.  An adaptive fast solver for a general class of positive definite matrices via energy decomposition.  Multiscale Modeling & Simulation, 16(2):615–678, 2018.
  * [33] T. Y. Hou, Q. Li, and P. Zhang.  Exploring the locally low dimensional structure in solving random elliptic pdes.  Multiscale Modeling & Simulation, 15(2):661–695, 2017.
  * [34] T. Y. Hou, D. Ma, and Z. Zhang.  A model reduction method for multiscale elliptic pdes with random coefficients using an optimization approach.  Multiscale Modeling & Simulation, 17(2):826–853, 2019.
  * [35] P. Jenny, S. Lee, and H. Tchelepi.  Multi-scale finite volume method for elliptic problems in subsurface flow simulation.  J. Comput. Phys., 187:47–67, 2003.
  * [36] P. Kidger and T. Lyons.  Universal Approximation with Deep Narrow Networks.  In J. Abernethy and S. Agarwal, editors, Proceedings of Thirty Third Conference on Learning Theory, volume 125 of Proceedings of Machine Learning Research, pages 2306–2327. PMLR, 09–12 Jul 2020.
  * [37] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization, 2017.
  * [38] K. Kunisch and S. Volkwein.  Galerkin proper orthogonal decomposition methods for parabolic problems.  Numerische Mathematik, 90:117–148, 2001.
  * [39] C. Le Bris, F. Legoll, and A. Lozinski.  An MsFEM type approach for perforated domains.  Multiscale Modeling & Simulation, 12(3):1046–1077, 2014.
  * [40] W. T. Leung, G. Lin, and Z. Zhang.  Nh-pinn: Neural homogenization based physics-informed neural network for multiscale problems.  arXiv preprint arXiv:2108.12942, 2021.
  * [41] W. Li, A. Alikhanov, Y. Efendiev, and W. T. Leung.  Partially explicit time discretization for nonlinear time fractional diffusion equations.  arXiv preprint arXiv:2110.13248, 2021.
  * [42] G. Lin, C. Moya, and Z. Zhang.  Accelerated replica exchange stochastic gradient langevin diffusion enhanced bayesian deeponet for solving noisy parametric pdes.  arXiv preprint arXiv:2111.02484, 2021.
  * [43] G. Lin, Y. Wang, and Z. Zhang.  Multi-variance replica exchange stochastic gradient mcmc for inverse and forward bayesian physics-informed neural network.  arXiv preprint arXiv:2107.06330, 2021.
  * [44] H. Owhadi and L. Zhang.  Metric-based upscaling.  Comm. Pure. Appl. Math., 60:675–723, 2007.
  * [45] A. Roberts and I. Kevrekidis.  General tooth boundary conditions for equation free modeling.  SIAM J. Sci. Comput., 29(4):1495–1510, 2007.
  * [46] G. Samaey, I. Kevrekidis, and D. Roose.  Patch dynamics with buffers for homogenization problems.  J. Comput. Phys., 213(1):264–287, 2006.
  * [47] H. Shi and Y. Li.  Local discontinuous galerkin methods with implicit-explicit multistep time-marching for solving the nonlinear cahn-hilliard equation.  Journal of Computational Physics, 394:719–731, 2019.
  * [48] Y. Yang, M. Ghasemi, E. Gildin, Y. Efendiev, and V. Calo.  Fast Multiscale Reservoir Simulations With POD-DEIM Model Reduction.  SPE Journal, 21(06):2141–2154, 06 2016.
  * [49] Z. Zhang, E. T. Chung, Y. Efendiev, and W. T. Leung.  Learning algorithms for coarsening uncertainty space and applications to multiscale simulations.  Mathematics, 8(5):720, 2020.