summaries: This paper represents a novel approach to learning embeddings for entities and relations from a knowledge base that are capable of directly encoding *reasoning* (i.e. FOL). The principal representation proposed for this representation is a *cone* (i.e. a set which is closed under multiplication by a positive scalar) which is motivated by the fact that cones are closed under intersection and complements. In practice, the authors use *sector cones*, which leads to some incongruity between the original motivation and the model in practice (eg. sector cones are not closed under intersection or union). Even so, the authors apply this embedding method to several standard logic embedding tasks and demonstrate significant performance improvements.

summaries: A recent line of knowledge graph (KG) query systems "compiles" the query into an embedded representation, which is interpreted as a geometric region and then used to find response entities whose embedding representations are contained within the region.  Query2box is an early example, where the regions are multidimensional boxes or rectangles.  To be useful, the regions have to be simple, closed under query operators, and capable of being characterized using a few parameters.  Union (disjunction) and negation have therefore created problem for boxes.  This paper replaces boxes with a Cartesian space of cones in 2d.  By restricting the family of permitted cones and disjunction to the last step of any query, the authors argue that their cone family is closed under query operators.  Experiments with benchmarks that have a good diversity of query graph structures show promising results.

summaries: The paper proposes a new framework for embedding first-order logic (FOL) queries over knowledge bases using region-based embeddings. The embedding space is a Cartesian product of 2D cones. Entities are embedded in the same space with 0-aperture cones (ray starting at origin).  - Authors provide the relevant foundation to cones in 2D space, the final Cartesian product space, and suggest parameterizations of projection, intersection, union, and complement operations. - The choice of using cones facilitates complement operation. - Authors demonstrate empirical improvements over past work in FOL query embedding. - While the proposed cone-based region embeddings allow for clear descriptions for intersection, union, and complement, the implementation does not strictly enforce these geometric concepts.

summaries: The paper presents a novel geometric embeddings called cone embeddings for solving query answering. Query answering requires embeddings of queries that are represented using first order logic. This includes conjunction, disjunction, negation but excluding Universal Quantification. Given this as basis, the related work until now has focused on two different types of embeddings (1) Geometric; and (2) probabilistic query based model. The progression from Geometric based embeddings to probabilistic models was primarily due to the fact that geometric embeddings cannot handle negation.

summaries: The paper proposes a new embedding model for multi-hop reasoning over knowledge graphs. The idea is to embed a complex query as n arc segments each on a 1-sphere. The paper designs several neural logic operators including the intersection, projection and complement operators. For union, the model uses the DNF technique to reorder the operators so that the union operators always appear in the last step. The authors evaluate the model on standard benchmarks and achieve better results than prior methods.   Overall the paper is clean and easy to follow. A general question I have is that the current modeling is to embed the query into several independent arc segments. Have you tried to embed the query into the surface (sphere cap) of a n-sphere?

abstractText: Query embedding (QE)—which aims to embed entities and first-order logical (FOL) queries in low-dimensional spaces—has shown great power in multi-hop reasoning over knowledge graphs. Recently, embedding entities and queries with geometric shapes becomes a promising direction, as geometric shapes can naturally represent answer sets of queries and logical relationships among them. However, existing geometry-based models have difficulty in modeling queries with negation, which significantly limits their applicability. To address this challenge, we propose a novel query embedding model, namely Cone Embeddings (ConE), which is the first geometry-based QE model that can handle all the FOL operations, including conjunction, disjunction, and negation. Specifically, ConE represents entities and queries as Cartesian products of two-dimensional cones, where the intersection and union of cones naturally model the conjunction and disjunction operations. By further noticing that the closure of complement of cones remains cones, we design geometric complement operators in the embedding space for the negation operations. Experiments demonstrate that ConE significantly outperforms existing state-of-the-art methods on benchmark datasets.

authors:
I. Balazevic
C. Allen
T. Hospedales

title: Multi-relational poincaré graph embeddings

venue: In Advances in Neural Information Processing Systems,

authors:
K. Bollacker
C. Evans
P. Paritosh
T. Sturge
J. Taylor

title: Freebase: A collaboratively created graph database for structuring human knowledge

venue: In Proceedings of the 2008 ACM SIGMOD International Conference on Management of Data,

authors:
A. Bordes
N. Usunier
A. Garcia-Duran
J. Weston
O. Yakhnenko

title: Translating embeddings for modeling multi-relational data

venue: In Advances in Neural Information Processing Systems

authors:
S. Boyd
L. Vandenberghe

title: Convex Optimization

authors:
R. Das
A. Neelakantan
D. Belanger
A. McCallum

title: Chains of reasoning over entities, relations, and text using recurrent neural networks

venue: In Proceedings of the 15th Conference of the European Chapter of the Association for Computational Linguistics: Volume 1, Long Papers,

authors:
S. Dasgupta
M. Boratko
D. Zhang
L. Vilnis
X. Li
A. McCallum

title: Improving local identifiability in probabilistic box embeddings

authors:
B.A. Davey
H.A. Priestley

title: Introduction to lattices and order

venue: Cambridge university press,

authors:
O. Ganea
G. Becigneul
T. Hofmann

title: Hyperbolic entailment cones for learning hierarchical embeddings

venue: In Proceedings of the 35th International Conference on Machine Learning,

authors:
D. Garg
S. Ikbal
S.K. Srivastava
H. Vishwakarma
H. Karanam
L.V. Subramaniam

title: Quantum embedding of knowledge for reasoning

authors:
S. Guo
Q. Wang
L. Wang
B. Wang
L. Guo

title: Knowledge graph embedding with iterative guidance from soft rules

venue: In Proceedings of the Thirty-Second AAAI Conference on Artificial Intelligence,

authors:
K. Guu
J. Miller
P. Liang

title: Traversing knowledge graphs in vector space

venue: In Proceedings of the Conference on Empirical Methods in Natural Language Processing,

authors:
W.L. Hamilton
P. Bajaj
M. Zitnik
D. Jurafsky
J. Leskovec

title: Embedding logical queries on knowledge graphs

venue: In Proceedings of the 32nd International Conference on Neural Information Processing Systems,

authors:
D.P. Kingma
J. Ba

title: Adam: A method for stochastic optimization

venue: In Proceedings of the International Conference on Learning Representations (ICLR),

authors:
T. Lacroix
N. Usunier
G. Obozinski

title: Canonical tensor decomposition for knowledge base completion

authors:
A. Lai
J. Hockenmaier

title: Learning to predict denotational probabilities for modeling entailment

authors:
Li
F.-L
M. Qiu
H. Chen
X. Wang
X. Gao
J. Huang
J. Ren
Z. Zhao
W. Zhao
L. Wang
G. Jin
Chu
W. Alime

title: assist: An intelligent assistant for creating an innovative e-commerce experience

venue: In Proceedings of the 2017 ACM on Conference on Information and Knowledge Management,

authors:
X.V. Lin
R. Socher
C. Xiong

title: Multi-hop knowledge graph reasoning with reward shaping

venue: In Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing,

authors:
Y. Lin
Z. Liu
M. Sun
Y. Liu
X. Zhu

title: Learning entity and relation embeddings for knowledge graph completion

venue: In Proceedings of the Twenty-Ninth AAAI Conference on Artificial Intelligence,

authors:
M. Nickel
D. Kiela

title: Poincaré embeddings for learning hierarchical representations

authors:
H. Ren
J. Leskovec

title: Beta embeddings for multi-hop logical reasoning in knowledge graphs

venue: In Neural Information Processing Systems,

authors:
H. Ren
W. Hu
J. Leskovec

title: Query2box: Reasoning over knowledge graphs in vector space using box embeddings

venue: In International Conference on Learning Representations,

authors:
A. Saxena
A. Tripathi
P. Talukdar

title: Improving multi-hop question answering over knowledge graphs using knowledge base embeddings

venue: In Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics,

authors:
H. Sun
T. Bedrax-Weiss
W. Cohen

title: PullNet: Open domain question answering with iterative retrieval on knowledge bases and text

venue: In Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP),

authors:
H. Sun
A.O. Arnold
T. Bedrax-Weiss
F. Pereira
W.W. Cohen

title: Faithful embeddings for knowledge base queries

authors:
Z. Sun
Deng
Z.-H
Nie
J.-Y
J. Tang

title: Rotate: Knowledge graph embedding by relational rotation in complex space

authors:
K. Toutanova
D. Chen

title: Observed versus latent features for knowledge base and text inference

venue: In 3rd Workshop on Continuous Vector Space Models and Their Compositionality,

authors:
T. Trouillon
J. Welbl
S. Riedel
E. Gaussier
G. Bouchard

title: Complex embeddings for simple link prediction

venue: In Proceedings of The 33rd International Conference on Machine Learning,

authors:
I. Vendrov
R. Kiros
S. Fidler
R. Urtasun

title: Order-embeddings of images and language

venue: In 4th International Conference on Learning Representations,

authors:
W. Xiong
T. Hoang
W.Y. Wang

title: DeepPath: A reinforcement learning method for knowledge graph reasoning

venue: In Proceedings of the 2017 Conference on Empirical Methods in Natural Language Processing,

authors:
Z. Zhang
J. Cai
J. Wang

title: Duality-induced regularizer for tensor factorization based knowledge graph completion

authors:
Z. Zhang
J. Cai
Y. Zhang
J. Wang

title: Learning hierarchy-aware knowledge graph embeddings for link prediction

venue: In Proceedings of the Thirty-Fourth AAAI Conference on Artificial Intelligence,

heading: 1 Introduction

text:
Multi-hop reasoning over knowledge graphs (KGs)—which aims to find answer entities of given queries using knowledge from KGs—has attracted great attention from both academia and industry recently [24, 23, 16]. In general, it involves answering first-order logic (FOL) queries over KGs using operators including existential quantification (∃), conjunction (∧), disjunction (∨), and negation (¬). A popular approach to multi-hop reasoning over KGs is to first transform a FOL query to its corresponding computation graph—where each node represents a set of entities and each edge represents a logical operation—and then traverse the KG according to the computation graph to identify the answer set. However, this approach confronts two major challenges. First, when some links are missing in KGs, it has difficulties in identifying the correct answers. Second, it needs to deal with all the intermediate entities on reasoning paths, which may lead to exponential computation cost.
To address these challenges, researchers have paid increasing attention to the query embedding (QE) technique, which embeds entities and FOL queries in low-dimensional spaces [12, 22, 21, 25]. QE models associate each logical operator in computation graphs with a logical operation in embedding
∗Corresponding author.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
spaces. Given a query, QE models generate query embeddings following the corresponding computation graph. Then, they determine whether an entity is a correct answer based on similarities between the query embeddings and entity embeddings.
Among the existing QE models, geometry-based models that embed entities and queries into geometric shapes have shown promising performance [12, 11, 5, 22]. Geometry-based models usually represent entity sets as "regions" (e.g., points and boxes) in Euclidean spaces and then design set operations upon them. For example, Query2Box [22] represents entities as points and queries as boxes. If a point is inside a box, then the corresponding entity is the answer to the query. Compared with non-geometric methods, geometric shapes provide a natural and easily interpretable way to represent sets and logical relationships among them.
However, existing geometry-based models have difficulty in modeling queries with negations, which significantly limits their applicability. For example, GQE [12] and Query2Box [22]—which embed queries to points and boxes, respectively—cannot handle queries with negation, as the complement of a point/box is no longer a point/box. To tackle this problem, Ren & Leskovec [21] propose a probabilistic QE model using Beta distributions. However, it does not have some advantages of geometric models. For example, using Beta distributions, it is unclear how to determine whether an entity is an answer to a query as that in the box case [22]. Therefore, proposing a geometric QE model that can model all the FOL queries is still challenging but promising.
In this paper, we propose a novel geometry-based query embedding model—namely, Cone Embeddings (ConE)—which represents entities and queries as Cartesian products of two-dimensional cones. Specifically, if the cones representing entities are subsets of the cones representing queries, then these entities are the answers to the query. To perform multi-hop reasoning in the embedding space, we define the conjunction and disjunction operations that correspond to the intersection and union of cones. Further, by noticing that the closure of complement of cones are still cones, we correspondingly design geometric complement operators in the embedding space for the negation operations. To the best of our knowledge, ConE is the first geometry-based QE model that can handle all the FOL operations, including conjunction, disjunction, and negation. Experiments demonstrate that ConE significantly outperforms existing state-of-the-art methods on benchmark datasets.

heading: 2 Related Work

text:
Our work is related to answering multi-hop logical queries over KGs and geometric embeddings.
Answering multi-hop logical queries over KGs. To answer multi-hop FOL queries, path-based methods [30, 17, 10] start from anchor entities and require traversing the intermediate entities on the path, which leads to exponential computation cost. Embedding-based models are another line of works, which embed FOL queries into low-dimensional spaces. For example, existing works embed queries to geometric shapes [12, 22, 11, 5], probability distributions [21], and complex objects [9, 25]. Our work also embeds queries to geometric shapes. The main difference is that our work can handle all the FOL operations, while existing works cannot.
Other Geometric embeddings. Geometric embeddings are popular in recent years. For example, geometric operations including translation [3, 18], rotation [28, 26, 14, 32], and complex geometric operations [1, 33] have been widely used in knowledge graph embeddings. Other geometric embedding methods also manage to use boxes [6], convex cones [15, 29, 19], etc. For example, Lütfü Özçep et al. [19] use axis-aligned cones to embed ontologies expressed in the ALC description logic, and use polars of cones to model negation operators. Recent years have also witnessed the development of embeddings in non-Euclidean geometry, such as Poincaré embeddings [1, 20] and hyperbolic entailment cones [8]. Notably, although there exist works that also use cone embeddings [15, 29, 8, 19], they are not designed for the multi-hop reasoning task and their definition of cones are different from that in our work.

heading: 3 Preliminaries

text: In this section, we review the background of query embeddings in Section 3.1 and introduce some basic concepts of two-dimensional cones in Section 3.2.

heading: 3.1 Backgrounds

text:
Knowledge Graphs (KGs). Given a set V of entities (vertices) and a set E of relations (edges), a knowledge graph G = {(si, pj , ok)} ⊂ V × E × V is a set of factual triple, where pj ∈ E is a predicate, and si, ok ∈ V are subject and object, respectively. Suppose that rj(·, ·) ∈ R is a binary function rj : V × V → {True,False} corresponding to pj , where rj(si, ok) = True if and only if (si, pj , ok) is a factual triples. Then, for all (si, pj , ok) ∈ G, we have rj(si, ok) = True. Note that both E and R are involved with relations, while E is a set of relation instances and R is a set of relational functions.
First-Order Logic (FOL). FOL queries in the query embedding literature involve logical operations including existential quantification (∃), conjunction (∧), disjunction (∨), and negation (¬). Universal quantification (∀) is not included, as no entity connects with all other entities in real-world KGs [21]. We use FOL queries in its Disjunctive Normal Form (DNF) [7], which represents FOL queries as a disjunction of conjunctions. To formulate FOL queries, we assume that Va ⊂ V is the non-variable anchor entity set, V1, . . . , Vk are existentially quantified bound variables, and V? is the target variable, i.e., the answers to a certain query. Then, a FOL query q in the disjunctive normal form is:
q[V?] = V?. ∃V1, . . . , Vk : c1 ∨ c2 ∨ · · · ∨ cn.
Specifically, ci are conjunctions, i.e., ci = ei1 ∧ · · · ∧ eim, where eij = r(va, V ) or ¬r(va, V ) or r(V ′, V ) or ¬r(V ′, V ), va ∈ Va, V ∈ {V?, V1, . . . , Vk}, V ′ ∈ {V1, . . . , Vk}, and V ̸= V ′. Using the aforementioned notations, answering a query q is equivalent to finding the set of entities JqK ⊂ V , where v ∈ JqK if and only if q[v] is True. Computation Graphs. Given a query, we represent the reasoning procedure as a computation graph (see Figure 1a for an example), of which nodes represent entity sets and edges represent logical operations over entity sets. We map edges to logical operators according to the following rules.
• Relation Traversal→Projection Operator P . Given a set of entities S ⊂ V and a relational function r ∈ R, the projection operator P outputs all the adjacent entities ∪v∈SN(v, r), where N(v, r) is the set of entities such that r(v, v′) = True for all v′ ∈ N(v, r).
• Conjunction→Intersection Operator I. Given n sets of entities {S1, S2, . . . , Sn}, the intersection operator I performs set intersection to obtain ∩ni=1Sn.
• Disjunction→Union Operator U . Given n sets of entities {S1, S2, . . . , Sn}, the union operator U performs set union to obtain ∪ni=1Sn.
• Negation→Complement Operator C. Given an entity set S ⊂ V , C gives S̄ = V\S.
Query Embeddings (QE). QE models generate low-dimensional continuous embeddings for queries and entities, and associate each logical operator for entity sets with an operation in embedding spaces. Since an entity is equivalent to a set with a single element and each query q is corresponding to a unique answer set JqK, the aim of QE models is equivalent to embedding entity sets that can be answers to some queries.

heading: 3.2 Cones in Two-Dimensional Spaces

text:
To represent FOL queries as Cartesian products of two-dimensional cones, we introduce some definitions about cones and the parameterization method of a special class of cones. Definition 1 (Boyd & Vandenberghe [4]). A set C ⊂ R2 is called a cone, if for every x ∈ C and λ ≥ 0, we have λx ∈ C . A set is a convex cone if it is convex and a cone, which means that for any x1, x2 ∈ C and λ1, λ2 ≥ 0, we have λ1x1 + λ2x2 ∈ C.
By letting λ = 0 in Definition 1, we know that a cone must contain the origin. In view of this property, we define a new operation called closure-complement for cones. Definition 2. Suppose that C ⊂ R2 is a cone. Then, the closure-complement of C is defined by C̃ = cl(R2\C), where cl(·) is the closure of a set.
Next, we introduce a class of cones that can be parameterized in a scalable way. Definition 3. A 2D closed cone is called a sector-cone, if its closure-complement or itself is convex.
The set of sector-cones is closed under closure-complement and their union and intersection are still cones. Besides, we have the following proposition, whose proof is provided in Appendix A. Proposition 1. A sector-cone is always axially symmetric.
Parameterization of 2D Sector-Cones. Proposition 1 suggests that we can use a pair of parameters to represent a two-dimensional sector-cone:
(θax, θap),where θax ∈ [−π, π), θap ∈ [0, 2π].
Specifically, θax represents the angle between the symmetry axis of the sector-cone and the positive x axis. θap represents the aperture of the sector-cone. For any points in the cone, its phase will be in [θax − θap/2, θax + θap/2]. Figure 1b gives examples of several (sector-)cones. One may notice that sector-cones share some similarities with boxes defined in Query2Box, which also involves region representations. However, we argue that sector-cones are more expressive than boxes, of which the details are provided in Appendix F.
Let K be the space consisting of all (θax, θap). We can represent an arbitrary sector-cone C0 as C0 = (θax, θap) ∈ K. Then, for a d-ary Cartesian product of sector-cones
C = C1 × C2 × . . .× Cd, (1)
we represent it via a d-dimensional vector in Kd: C = ( (θ1ax, θ 1 ap), . . . , (θ d ax, θ d ap) ) ⊂ Kd. (2)
where θiax ∈ [−π, π), θiap ∈ [0, 2π], for i = 1, . . . , d. Or equivalently, C = (θax,θap), where θax = (θ 1 ax, . . . , θ d ax) ∈ [−π, π)d and θap = (θ1ap, . . . , θdap) ∈ [0, 2π]d.

heading: 4 Cone Embeddings

text: In this section, we propose Cone Embeddings (ConE) for multi-hop reasoning over KGs. We first introduce cone embeddings for conjunctive queries and entities in Section 4.1. Afterwards, we introduce the logical operators and the methods to learn ConE in Sections 4.2 and 4.3.

heading: 4.1 Cone Embeddings for Conjunctive Queries and Entities

text:
As introduced in Section 3.1, conjunctive queries constitute the basis of all queries in the DNF form. Embeddings of all queries can be generated by applying logical operators to conjunctive queries’ embeddings. Thus, we design embeddings for conjunctive queries in this section. We model queries with disjunction using the Union Operator U in Section 4.2. In general, the answer entities to a conjunctive query q have similar semantics. For example, answers to the query that "List all the directors of American movies" should all be persons; answers to the query that "List all the Asian cities that ever held Olympic Games" should all be places. If we embed an entity set JqK into an embedding space, we expect entities in JqK to have similar embeddings. Thus,
we expect their embeddings to form a "region" in the embedding space. If the embedding of an entity is inside the region, then the entity is likely to be an answer. Further, we can find a semantic center and a boundary for the region, where the semantic center represents the semantics of JqK and the boundary designates how many entities are in JqK. To model the embedding region of JqK, we propose to embed it to a Cartesian product of sector-cones. Specifically, we use the parameter θiax to represent the semantic center, and the parameter θ i ap to determine the boundary of JqK. If we use a d-ary Cartesian product, i.e., the embedding dimension is d, we define the embedding of JqK as
Vcq = (θax,θap),
where θax ∈ [−π, π)d are axes and θap ∈ [0, 2π]d are apertures. An entity v ∈ V is equivalent to an entity set with a single element, i.e., {v}. We propose to represent an entity as a Cartesian product of cones with apertures 0, where the axes indicates the semantics of the entity. Formally, if the embedding dimension is d, the cone embedding of v is v = (θax, 0), where θax ∈ [−π, π)d is the axis embedding and 0 is a d-dimensional vector with all elements being 0.

heading: 4.2 Logical Operators for Cone Embeddings

text:
In this section, we introduce our designed logical operators of ConE in the embedding space, including projection, intersection, union, and complement.
It is worth noting that, the composition of logical operators may lead to non-sense queries. For example, the queries "List the intersection/union of American movies and non-American movies" and "List the intersection of American movies and Asian movies" make no sense in real-world applications. However, the main aim of a query embedding model is to represent all entity sets that can be answer to some real-world query. Therefore, we do not need to model the entity sets that only correspond to theoretically possible queries [21].
Projection Operator P . The goal of P is to represent an entity’s adjacent entities that are linked by a given relation. It maps an entity set to another entity set (see Figure 2a). Thus, we define a relation-dependent function in the embedding space for P:
fr : Kd → Kd, Vcq 7→ V c′ q .
We implement fr in a neural way. First, we represent relations as relational translations of query embeddings and assign each relation with an embedding r = (θax,r,θap,r). Then, we define fr as
fr(Vq) = g(MLP([θax + θax,r;θap + θap,r])) (3)
where MLP : R2d → R2d is a multi-layer perceptron network, [·; ·] is the concatenation of two vectors, and g is a function that generates θ′ax ∈ [−π, π)d and θ′ap ∈ [0, 2π]d. We define g as:
[g(x)]i = { θ′iax = π tanh(λ1xi), if i ≤ d, θ′i−dap = π tanh(λ2xi) + π, if i > d.
where [g(x)]i denotes the i-th element of g(x), λ1 and λ2 are two fixed parameters to control the scale. Note that the range of the hyperbolic tangent function (tanh) are open sets. Thus, we cannot indeed get the boundary value θ′iax = −π and θ′iap = 0, 2π. However, when we implement g in experiments, the value of g can be very close to 0 and 2π, which is equivalent to the closed set numerically.
Intersection Operator I. Given a query q that is the conjunction of conjunctive queries qi, the goal of I is to represent JqK = ∩ni=1JqiK. Since the conjunction of conjunctive queries are still conjunctive
queries, the entities in JqK should have similar semantics. Recall that we only need to model entity sets that can be answers. We still use a Cartesian product of sector-cones to represent JqK (see Figure 2b). Suppose that Vcq = (θax,θap) and V c qi = (θi,ax,θi,ap) are cone embeddings for JqK and JqiK, respectively. We define the intersection operator as follows:
θax = SemanticAverage(Vcq1 , . . . ,V c qn), θap = CardMin(Vcq1 , . . . ,V c qn),
where SemanticAverage(·) and CardMin(·) generates semantic centers and apertures, respectively. In the following, we introduce these two functions in detail.
SemanticAverage. As the semantic center of Vcq, θax should be close to all the semantic centers θi,ax. Thus, we propose to represent θax as a semantic average of θi,ax. We note that the ordinary weighted average may lead to inconsistent semantics. For example, when d = 1, if θ1,ax = π − ϵ and θ2,ax = −π + ϵ (0 < ϵ < π/4), then we expect θax to be around π. However, if we use the ordinary weighted sum, θax will be around 0 with a high probability. To tackle this issue, we propose a semantic average scheme, which takes periodicity of axes into account. For a figure illustration of the difference between the ordinary and semantic average, please refer to Appendix D.
Specifically, we first map [θi,ax]j to points on the unit circle. Then, compute the weighted average of the points using an attention mechanism. Finally, map the points back to angles that represent axes. Formally, the computation process is
[x;y] = n∑ i=1 [ai ◦ cos(θi,ax); ai ◦ sin(θi,ax)],
θax = Arg(x,y),
where cos and sin are element-wise cosine and sine functions; ai ∈ Rd are positive weights vectors that satisfy ∑n i=1[ai]j = 1 for all j = 1, . . . , d; ◦ is the element-wise multiplication; Arg(·) is the function that computes arguments of (a Cartesian of) 2D points. Noticing that the weights ai are relevant to both axes and apertures, we compute [ai]j via the following attention mechanism:
[ai]j = exp([MLP([θi,ax − θi,ap/2;θi,ax + θi,ap/2])]j)∑n
k=1 exp([MLP([θk,ax − θk,ap/2;θk,ax + θk,ap/2])]j) ,
where MLP : R2d → Rd is a multi-layer perceptron network, [·; ·] is the concatenation of two vectors. We can see θi,ax − θi,ap/2 and θi,ax + θi,ap/2 as the lower and upper bound of sector-cones. We use Arg(·) to recover angles of 2D points. Suppose that βi = arctan([y]i/[x]i), then
[θax]i =  βi + π, if [x]i < 0, [y]i > 0, βi − π, if [x]i < 0, [y]i < 0,
βi, otherwise.
Note that [x]i = 0 will lead to an illegal division. In experiments, we manually set [x]i to be a small number (e.g., 10−3) when [x]i = 0.
CardMin. Since JqK is the subset of all JqiK, θiap should be no larger than any apertures θij,ap. Therefore, we implement CardMin by a minimum mechanism with cardinality constraints:
θiap = min{θi1,ap, . . . , θin,ap} · σ([DeepSets({Vqj}nj=1)]i),
where σ(·) is the element-wise sigmoid function, θij,ap is the i-th element of θj,ap, DeepSets(·) is a permutation-invariant function [31]. Specifically, DeepSets({Vqj}nj=1) is computed by
MLP ( 1
n ∑n j=1 MLP ([θj,ax − θj,ap/2;θj,ax + θj,ap/2]) ) .
Union Operator U . Given a query q that is the disjunction of conjunctive queries qi, the goal of the union operator U is to represent JqK = ∪ni=1JqiK. As noted by Ren et al. [22], directly modeling the disjunction leads to unscalable models. Thus, we adopt the DNF technique [22], in which the union operation only appears in the last step in computation graphs.
Suppose that Vcqi = (θi,ax,θi,ap) are cone embeddings for JqiK. To represent the union of several cones (see Figure 2c), we represent JqK as a set of Vcqi :
Vdq = {V c q1 , . . . ,V c qn},
where n may be various in different queries. Equivalently, Vdq can be written as Vdq = ( {(θ11,ax, θ11,ap), . . . , (θ1n,ax, θ1n,ap)}, . . . , {(θd1,ax, θd1,ap), . . . , (θdn,ax, θdn,ap)} ) .
As {(θi1,ax, θi1,ap), . . . , (θin,ax, θin,ap)} are the union of d sector-cones, it is also a cone. Thus, the cone embedding of q is also a Cartesian product of two-dimensional cones.
Complement Operator C. Given an conjunctive query q and the corresponding entity set JqK, the aim of C is to identify the set J¬qK, which is the complementary of JqK, i.e., V\JqK. Since the set of sector-cones is closed under closure-complement, we define C using the closure-complement. Thus, the apertures of Vq plus the apertures of V¬q should be a vector with all elements being 2π. Moreover, to represent the semantic difference between JqK and J¬qK, we assume that their semantic centers to be opposite. Please refer to Figure 2d for a figure illustration.
Suppose that Vq = (θax,θap) and V¬q = (θ′ax,θ′ap). We define the complement operator C as:
[θ′ax]i = { [θax]i − π, if [θax]i ≥ 0, [θax]i + π, if [θax]i < 0,
[θ′ap]i = 2π − [θap]i.

heading: 4.3 Learning Cone Embeddings

text:
To learn cone embeddings, we expect that the cone embeddings of entities v ∈ JqK are inside the cone embeddings of q, and the cone embeddings of entities v′ /∈ JqK are far from the cone embedding of q. This motivates us to define a distance function to measure the distance between a given query embedding and an entity embedding, and a training objective with negative sampling.
Distance Function. We first define the distance function for conjunctive queries. Inspired by Ren et al. [22], we divide the distance d into two parts—the outside distance do and the inside distance di. Figure 2e gives an illustration of the distance function d. Suppose that v = (θvax, 0), V c q = (θax,θap), θL = θax − θap/2 and θU = θax + θap/2. We define the distance as dcon(v;V c q) = do(v;V c q) + λdi(v;V c q).
The outside distance and the inside distance are do = ∥min {|sin (θvax − θL) /2| , |sin (θvax − θU ) /2|}∥1 , di = ∥min {|sin (θvax − θax) /2| , |sin (θap) /2|}∥1 ,
where ∥ · ∥1is the L1 norm, sin(·) and min(·) are element-wise sine and minimization functions. Note that as axes and apertures are periodic, we use the sine function to enforce two close angles have small distance. The parameter λ ∈ (0, 1) is fixed during training, so that v is encouraged to be inside the cones represented by Vcq , but not necessarily be equal to the semantic center of V c q .
Since we represent the disjunctive queries as a set of embeddings, we cannot use dcon to directly compute the distance. Nonetheless, the distance between a point and the union of several sets is equal to the minimum distance between the point and each of those sets. Therefore, for a query q = q1 ∨ · · · ∨ qn in the Disjunctive Normal Form, the distance between q and an entity is
ddis(v;Vdq) = min{dcon(v;V c q1), . . . , dcon(v;V c qn)}.
If we use Vq to represent embeddings of both kinds of queries, the unified distance function d is
d(v;Vq) = { dcon(v;Vq), if q is conjunctive queries, ddis(v;Vq), if q is disjunctive queries.
Training Objective. Given a training set of queries, we optimize a negative sampling loss
L = − log σ(γ − d(v;Vq))− 1
k ∑k i=1 log σ(d(v′i;Vq)− γ),
where γ > 0 is a fixed margin, v ∈ JqK is a positive entity, v′i /∈ JqK is the i-th negative entity, k is the number of negative entities, and σ(·) is the sigmoid function.

heading: 5 Experiments

text: In this section, we conduct experiments to demonstrate that: 1) ConE is a powerful model for the multi-hop reasoning over knowledge graphs; 2) the aperture embeddings of ConE are effective in modeling cardinality (i.e., the number of elements) of answer sets. We first introduce experimental settings in Section 5.1 and then present the experimental results in Sections 5.2 and 5.3. The code of ConE is available on GitHub at https://github.com/MIRALab-USTC/QE-ConE.

heading: 5.1 Experimental Settings

text:
We adopt the commonly used experimental settings for query embeddings [12, 22, 21].
Datasets and Queries. We use three datasets: FB15k [2], FB15k-237 (FB237) [27], and NELL995 (NELL) [30]. QE models focus on answering queries involved with incomplete KGs. Thus, we aim to find non-trivial answers to FOL queries that cannot be discovered by traversing KGs. For a fair comparison, we use the same query structures as those in Ren & Leskovec [21]. The training and validation queries consist of five conjunctive structures (1p/2p/3p/2i/3i) and five structures with negation (2in/3in/inp/pni/pin). We also evaluate models’ generalization ability, i.e., answering queries with structures that models have never seen during training. The extra query structures include ip/pi/2u/up. Please refer to Appendix B.1 for more details about datasets and query structures.
Training Protocol. We use Adam [13] as the optimizer, and use grid search to find the best hyperparameters based on the performance on the validation datasets. For the search range and best hyperparameters, please refer to Appendix B.2.
Evaluation Protocol. We use the same evaluation protocol as that in Ren & Leskovec [21]. We first build three KGs: the training KG Gtrain, the validation KG Gvalid, and the test KG Gtest using training edges, training+validation edges, training+validation+test edges, respectively. Given a test (validation) query q, we aim to discover non-trivial answers JqKtest\JqKvalid (JqKvalid\JqKtrain). In other words, to answer an entity, we need to impute at least one edge to create an answer path to it. For each non-trivial answer v of a test query q, we rank it against non-answer entities V\JqKtest. We denote the rank as r and calculate the Mean Reciprocal Rank (MRR), of which the definition is provided in Appendix B.3. Higher MRR indicates better performance.
Baselines. We compare ConE against three state-of-the-art models, including GQE [12], Query2Box (Q2B) [22], and BETAE [21]. GQE and Q2B are trained only on five conjunctive structures as they cannot model the queries with negation. Since the best embedding dimension d for ConE is 800, we retrain all the baselines with d = 800. The results of GQE and Q2B are better than those reported in Ren & Leskovec [21], while the results of BETAE become slightly worse. Therefore, we reported the results of GQE and Q2B with d = 1600 and BETAE with d = 400. For the results of BETAE with d = 800, please refer to Appendix C.1.

heading: 5.2 Main Results

text:
We compare ConE against baseline models on queries with and without negation. We run our model five times with different random seeds and report the average performance. For the error bars of the performance, please refer to Appendix C.5.
Queries without Negation. Table 1 shows the experimental results on queries without negation, i.e., existentially positive first-order (EPFO) queries, where AVG denotes average performance. Overall, ConE significantly outperforms compared models. ConE achieves on average 19.7%, 12.0%, and 10.6% relative improvement MRR over previous state-of-the-art BETAE on the three datasets, which demonstrates the superiority of geometry-based models. Compared with Q2B, which uses Query2Box to embed queries, ConE gains up to 24.2% relative improvements. ConE also gains an impressive improvement on queries ip/pi/2u/up, which are not in the training graph. For example, ConE outperforms BETAE by 38.9% for 2u query on FB15k. The results show the superior generality ability of ConE. Since ConE is capable of modeling complement, we can also implement disjunctive queries using De Morgan’s law. However, using De Morgan’s law always results in sector-cones, which may be inconsistent with the real set union. Thus, the models with DNF outperforms those with De Morgan’s law. We include the detailed results in Appendix C.2 due to the space limit.
Queries with Negation. Table 2 shows the results of ConE against BETAE on modeling FOL queries with negation. Since GQE and Q2B are not capable of handling the negation operator, we do not include their results in the experiments. Overall, ConE outperforms BETAE by a large margin. Specifically, ConE achieves on average 25.4%, 9.3%, and 8.5% relative improvement MRR over BETAE on FB15k, FB237, and NELL, respectively.

heading: 5.3 Modeling the Cardinality of Answer Sets

text:
As introduced in Section 4.1, the aperture embeddings can designate the cardinality (i.e., the number of elements) of JqK. In this experiment, we demonstrate that although we do not explicitly enforce ConE to learn cardinality during training, the learned aperture embeddings are effective in modeling the cardinality of answer sets. The property partly accounts for the empirical improvements of ConE.
We compute the correlations between learned aperture embeddings and the cardinality of answer sets. Specifically, for the cone embedding Vq = (θax,θap) of a given query q, we use the L1 norm of θap to represent the learned cardinality of JqK. Then, we compute the Spearman’s rank correlation (SRC) between the learned cardinality and the real cardinality, which measures the statistical dependence between the ranking of two variables. Higher correlation indicates that the embeddings can better model the cardinality of answer sets. As we model queries with disjunction using the DNF technique, we do not include the results of disjunctive queries following Ren & Leskovec [21].
Table 3 shows the results of SRC for ConE, Query2Box (Q2B), and BETAE on FB15k. For the results on FB237 and NELL, please refer to Appendix C.3. As Query2Box cannot handle queries with negation, we do not include its results on these queries. On all query structures, ConE outperforms the previous state-of-the-art method BETAE. Note that BETAE is a probabilistic model, of which the authors claim that it can well handle the uncertainty of queries, i.e., the cardinality of answer set. Nonetheless, ConE still outperforms BETAE by a large margin, which demonstrates the expressiveness of cone embeddings. We also conduct experiments using Pearson’s correlation, which measures the linear correlation between two variables. Please refer to Appendix C.3 for the results.

heading: 5.4 The Designed Operators and the Real Set Operations

text:
As introduced in Section 4.2, the designed union (using DNF technique) and complement operators for ConE are non-parametric. They correspond to exact set union and complement. Meanwhile, we define neural operators to approximate the projection and intersection operators to achieve a tractable training process and better performance. Notably, the designed neural operators may not exactly match the real set operations. However, experiments on some example cases demonstrate that these neural operators provide good approximations for real set operations. In the following, we show the experimental results for the operators including projection and intersection. In all experiments, the ConE embeddings are trained on FB15k.
Projection. Suppose that a set A is included by a set B, then we expect the projection of A is also included by the projection of B. We randomly generate 8000 pairs of sector-cones (Ai, Bi), where Ai ⊂ Bi. Then, for each i, we randomly select a relation ri and calculate the projections Pri(Ai) and Pri(Bi). Ideally, the projected cones should satisfy Pri(Ai) ⊂ Pri(Bi). We calculate the ratio ri = |Pri(Ai)∩Pri(Bi)|/|Pri(Ai)| to measure how many elements in Pri(Ai) are included in Pri(Ai)∩Pri(Bi). Finally, we get an average ratio r = 0.8113. That is to say, the learned Pri(Ai) are included in Pri(Ai) with a high probability. The learned projection operators approximate the real set projection well.
Intersection. To validate that the learned intersection can well approximate real set intersection, we randomly generate 8000 pairs of sector-cones (Ci, Di), where Ci ∩ Di is not guaranteed to be a sector-cone. Then, we generate embeddings I(Ci, Di) for the intersection Ci ∩ Di. Ideally, the learned cones I(Ci, Di) should be the same as the real cones Ci ∩Di. We calculate the ratio ri = |I(Ci, Di) ∩ (Ci ∩Di)|/|I(Ci, Di) ∪ (Ci ∩Di)| to measure the overlap between I(Ci, Di) and Ci ∩Di, and obtain an average ratio of r = 0.6134. Note that the experiments are conducted on the test set, and we did not explicitly train our model on these queries. The relatively high overlap ratio r demonstrates that the learned intersection is a good approximation of the real set intersection.
We further conduct experiments to demonstrate that the learned intersection operators can well handle empty intersections. Following Ren et al. [22], on FB15k, we randomly generate 10k queries of two types: (a) intersection queries with more than five answers, and (b) intersection queries with empty answer sets. We found that the average aperture is 2.495 for type (a) queries, while 1.737 for type (b) queries. The results demonstrate that although we have never trained ConE on the type (b) queries, the empty intersection sets are much more likely to have smaller apertures than queries with non-zero answers (with a 0.9632 ROC-AUC score). In other words, though we did not train ConE on datasets with empty intersection sets, we can distinguish empty answer sets by the learned apertures.
We also conduct experiments to demonstrate the difference between the learned union operator with De Morgan’s law and the real set union. Please refer to Appendix C.6 for details.

heading: 6 Conclusion

text:
In this paper, we propose a novel query embedding model, namely Cone Embeddings (ConE), to answer multi-hop first-order logical (FOL) queries over knowledge graphs. We represent entity sets as Cartesian products of cones and design corresponding logical operations. To the best of our knowledge, ConE is the first geometric query embedding models that can model all the FOL operations. Experiments demonstrate that ConE significantly outperforms previous state-of-the-art models on benchmark datasets. One future direction is to adapt ConE to queries in the natural language, which will further improve ConE’s applicability.
Acknowledgments and Disclosure of Funding
We would like to thank all the anonymous reviewers for their insightful comments. This work was supported in part by National Science Foundations of China grants 61822604, U19B2026, 61836006, and 62021001, and the Fundamental Research Funds for the Central Universities grant WK3490000004.

title: ConE: Cone Embeddings for Multi-Hop Reasoning over Knowledge Graphs


summaries: This paper studies the problem of dynamic programming with continuous state and action spaces.  The authors first propose the VI which naively discretizes the state and action space and performs the classical value iteration. The authors use this method as a baseline for comparison with CVI developed in this paper.  Regarding the CVI, this algorithm performs VI in the the conjugate domain. In particular, by taking the dual of the objective function, we take the problem to the dual space. By doing so, instead of taking minimum over all the actions, the problem is reduced to a summation. Hence, the complexity of the CVI is better than the VI in the order of the cardinality of the action space.

summaries: The authors consider value iteration [TJ](x) = min_u [c(x,u) + gamma E J(x’)] for problems where the transition function and cost functions are additively separable: c(x,u)=c1(x)+c2(u) and x’=f1(x)+f2(u) + w (Assumptions 3.1). Furthermore, f1(x) is Lipschitz continuous, f2(u) is linear, c1(x) is convex and Lipschitz continuous, and c2(u) is linear in u. The disturbance w has a finite support (Assumption 2.1). The state space for x and decision space for u are continuous. Thus, the authors consider problems such that [TJ](x) = c1(x) + min_u [c2(u) + gamma sum_i  p_i J(f1(x)+Bu+w_i)]. A problem that fits this description is the stabilization of an inverted pendulum with discretized disturbances, studied in Section 4. With this choice of assumptions, and assuming the value function is convex (which holds true if f1(x) is linear or linearized, that is, if x’=Ax+Bu+w), convex duality can be applied to the minimization problem over (u,x”) given x and subject to x”=Ax+Bu, which gives another representation of the optimization problem, based on the convex conjugate of c1 and the convex conjugate of e(x”)=gamma E J(x”+w).  The authors propose and study an algorithm that relies on computing these convex conjugate functions numerically in order to provide an alternative way to the computation of the value function J. They establish error bounds (Theorem 3.13) in terms of theoretical quantities of the problem data.

summaries: This paper presents a fast method for solving certain classes of optimal control problems, by implementing a value iteration (VI) algorithm that takes advantage of convex duality.  The basic idea is to take advantage of the fact that for two functions $f_1, f_2$, we have $(f_1 \square f_2)^* = f_1^* + f_2^*$, where $f^*$ denotes the Legendre-Fenchel transform of $f$, and $(f_1 \square f_2)(x) = \inf \\{f_1(x_1) + f_2(x_2) : x_1 + x_2 = x \\} $ denotes the infimal convolution of $f_1, f_2$. Thus, the slow $\inf$ operation can be replaced with a fast addition operation. The paper provides analyses of the convergence, time complexity, and error of the algorithm. When compared to a "naive" VI algorithm in which each iteration takes $O(X U)$ (where $X$, $U$ are the grid sizes of the state and input spaces, respectively), the algorithm presented in this paper can achieve $O(X+U)$. Finally, the paper provides a numerical example implemented in MATLAB.

summaries: This paper proposed a new value iteration algorithm (CVI) for optimal control of stochastic system with continuous state space. The proposed algorithm discretizes the state space and performs the value iteration step in the "conjugate domain" to reduce one-step computational complexity. It was shown the proposed algorithm converges under conditions on the stochastic system, discretization, and etc. The proposed algorithm CVI, its variant CVI-d are compared to classical value iteration (VI) over a noisy inverted pendulum experiment.

summaries: The paper develops an approximate value iteration method for an infinite-horizon, discounted-cost Markov Decision Processes (MDPs) that satisfy a given set of regulatory assumptions. The idea is to work on a dual space that replaces the DP value function by a reformulation written in terms of biconjugate operations. The authors present an algorithm that applies the conjugate operators iteratively, and present convergence and other structural results associated with the resulting errors. Numerical results compare the proposed approach against a traditional value iteration and a variant that generated dynamic discretization grids.

abstractText: In this study, we consider the infinite-horizon, discounted cost, optimal control of stochastic nonlinear systems with separable cost and constraints in the state and input variables. Using the linear-time Legendre transform, we propose a novel numerical scheme for implementation of the corresponding value iteration (VI) algorithm in the conjugate domain. Detailed analyses of the convergence, time complexity, and error of the proposed algorithm are provided. In particular, with a discretization of size X and U for the state and input spaces, respectively, the proposed approach reduces the time complexity of each iteration in the VI algorithm from O(XU) to O(X + U), by replacing the minimization operation in the primal domain with a simple addition in the conjugate domain.

authors:
Y. Achdou
F. Camilli
L. Corrias

title: On numerical approximation of the Hamilton-Jacobi-transport system arising in high frequency approximations

venue: Discrete & Continuous Dynamical Systems-Series B,

authors:
M. Akian
S. Gaubert
A. Lakhoua

title: The max-plus finite element method for solving deterministic optimal control problems: Basic properties and convergence analysis

venue: SIAM Journal on Control and Optimization,

title: Max-plus matching pursuit for deterministic Markov decision processes

venue: arXiv preprint arXiv:1906.08524,

authors:
R. Bellman
W. Karush

title: Mathematical programming and the maximum transform

venue: Journal of the Society for Industrial and Applied Mathematics,

authors:
E. Berthier
F. Bach

title: Max-plus linear approximations for deterministic continuous-state markov decision processes

venue: IEEE Control Systems Letters,

title: Linear convex stochastic control problems over an infinite horizon

venue: IEEE Transactions on Automatic Control,

title: Dynamic Programming and Optimal Control, Vol. II

venue: Athena Scientific,

title: Convex Optimization Theory

title: Reinforcement Learning and Optimal Control

authors:
L. Busoniu
R. Babuska
B. De Schutter
D. Ernst

title: Reinforcement learning and dynamic programming using function approximators

venue: CRC press,

authors:
R. Carpio
T. Kamihigashi

title: Fast value iteration: an application of Legendre-Fenchel duality to a class of deterministic dynamic programming problems in discrete time

venue: Journal of Difference Equations and Applications,

authors:
L. Contento
A. Ern
R. Vermiglio

title: A linear-time approximate convex envelope algorithm using the double Legendre–Fenchel transform with application to phase separation

venue: Computational Optimization and Applications,

title: Fast Legendre-Fenchel transform and applications to Hamilton-Jacobi equations and conservation laws

venue: SIAM Journal on Numerical Analysis,

title: Lebacque. A variational formulation for higher order macroscopic traffic flow models: Numerical investigation

venue: Transportation Research Part B: Methodological,

authors:
A.O. Esogbue
C.W. Ahn

title: Computational experiments with a class of dynamic programming algorithms of higher dimensions

venue: Computers & Mathematics with Applications,

authors:
P.F. Felzenszwalb
D.P. Huttenlocher

title: Distance transforms of sampled functions

venue: Theory of computing,

authors:
M. Jacobs
F. Léger

title: A fast approach to optimal transport: The back-and-forth method

venue: Numerische Mathematik,

authors:
C.M. Klein
T.L. Morin

title: Conjugate duality and the curse of dimensionality

venue: European Journal of Operational Research,

authors:
A.S. Kolarijani
S.C. Bregman
P. Mohajerin Esfahani
T. Keviczky

title: A decentralized event-based approach for robust model predictive control

authors:
M.A.S. Kolarijani
P. Mohajerin Esfahani

title: Fast approximate dynamic programming for input-affine dynamics

title: Conjugate value iteration (ConjVI) MATLAB package. Licensed under the MIT License, available online at https://github.com/AminKolarijani/ ConjVI, 2021

title: Discrete conjugate dynamic programming (d-CDP) MATLAB package. Licensed under the MIT License, available online at https://github.com/ AminKolarijani/d-CDP, 2021

title: Faster than the fast Legendre transform, the linear-time Legendre transform

venue: Numerical Algorithms,

title: New sequential exact Euclidean distance transform algorithms based on convex analysis

venue: Image and Vision Computing,

title: Max-plus eigenvector representations for solution of nonlinear H∞ problems: basic concepts

title: Approximate Dynamic Programming: Solving the Curses of Dimensionality

authors:
A. Sidford
M. Wang
X. Wu
Y. Ye

title: Variance reduced value iteration and faster algorithms for solving Markov decision processes

venue: In Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms,

authors:
R.S. Sutton
A.G. Barto

title: Reinforcement Learning: An Introduction

heading: 1 Introduction

text:
Value iteration (VI) is one of the most basic and wide-spread algorithms employed for tackling problems in reinforcement learning (RL) and optimal control [9, 28] formulated as Markov decision processes (MDPs). The VI algorithm simply involves the consecutive applications of the dynamic programming (DP) operator
T J(xt) = min ut
{ C(xt, ut) + γEJ(xt+1) } ,
where C(xt, ut) is the cost of taking the control action ut at the state xt. This fixed point iteration is known to converge to the optimal value function for discount factors γ ∈ (0, 1). However, this algorithm suffers from a high computational cost for large-scale finite state spaces. For problems with a continuous state space, the DP operation becomes an infinite-dimensional optimization problem, rendering the exact implementation of VI impossible in most cases. A common solution is to incorporate function approximation techniques and compute the output of the DP operator for a finite sample (i.e., a discretization) of the underlying continuous state space. This approximation again suffers from a high computational cost for fine discretizations of the state space, particularly in high-dimensional problems. We refer the reader to [9, 26] for various approximations of VI.
For some problems, however, it is possible to partially address this issue by using duality theory, i.e., approaching the minimization problem in the conjugate domain. In particular, as we will see
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
in Section 3, the minimization in the primal domain in DP can be transformed to a simple addition in the dual domain, at the expense of three conjugate transforms. However, proper application of this transformation relies on efficient numerical algorithms for conjugation. Fortunately, such an algorithm, known as linear-time Legendre transform (LLT), has been developed in late 90s [23]. Other than the classical application of LLT (and other fast algorithms for conjugate transform) in solving Hamilton-Jacobi equation [1, 13, 14], these algorithms are used in image processing [24], thermodynamics [12], and optimal transport [17].
The application of conjugate duality for the DP problem is not new and actually goes back to Bellman [4]. Further applications of this idea for reducing the computational complexity were later explored in [15, 18]. However, surprisingly, the application of LLT for solving discrete-time optimal control problems, has been limited. In particular, in [11], the authors propose the “fast value iteration” algorithm (without a rigorous analysis of the complexity and error of the proposed algorithm) for a particular class of infinite-horizon optimal control problems with state-independent stage cost C(x, u) = C(u) and deterministic linear dynamics xt+1 = Axt + But, where A is a non-negative, monotone, invertible matrix. More recently, in [20], we also considered the application of LLT for solving the DP operation in finite-horizon, optimal control of input-affine dynamics xt+1 = fs(xt) + But with separable cost C(x, u) = Cs(x) + Ci(u). In particular, we introduced the “discrete conjugate DP” (d-CDP) operator, and provided a detailed analysis of its complexity and error. As we will discuss shortly, the current study is an extension of the corresponding d-CDP algorithm that, among other things, considers infinite horizon, discounted cost problems. We note that the algorithms developed in [16, 24] for distance transform can also potentially tackle the optimal control problems similar to the ones of interest in the current study. In particular, these algorithms require the stage cost to be reformulated as a convex function of the “distance” between the current and next states. While this property might arise naturally, it can generally be restrictive, as it is in the problem class considered in this study. Another line of work that is closely related to ours invloves utilizing max-plus algebra in solving deterministic, continuous-state, continuous-time, optimal control problems; see, e.g., [2, 25]. These works exploit the compatibility of the DP operation with maxplus operations, and approximate the value function as a max-plus linear combination. Recently, in [3, 5], the authors used this idea to propose an approximate VI algorithm for continuous-state, deterministic MDPs. In this regard, we note that the proposed approach in the current study also involves approximating the value function as a max-plus linear combination, namely, the maximum of affine functions. The key difference is however that by choosing a grid-like (factorized) set of slopes for the linear terms (i.e., the basis of the max-plus linear combination), we take advantage of linear time complexity of LLT in computing the constant terms (i.e., the coefficients of the max-plus linear combination).
Main contribution. In this study, we focus on an approximate implementation of VI involving discretization of the state and input spaces for solving the optimal control problem of discrete-time systems, with continuous state-input space. Building upon our earlier work [20], we employ conjugate duality to speed-up VI for problems with separable stage cost (in state and input) and input-affine dynamics. We propose the conjugate VI (ConjVI) algorithm based on a modified version of the d-CDP operator introduced in [20], and extend the existing results in three directions: We consider infinite-horizon, discounted cost problems with stochastic dynamics, while incorporating a numerical scheme for approximation of the conjugate of input cost. The main contributions of this paper are then as follows: (i) we provide sufficient conditions for the convergence of ConjVI (Theorem 3.9); (ii) we show that ConjVI can achieve a linear time complexity of O(X+U) in each iteration (Theorem 3.10), compared to the quadratic time complexity of O(XU) of the standard VI, where X and U are the cardinalities of the discrete state and input spaces, respectively; (iii) we analyze the error of ConjVI (Theorem 3.11), and use that result to provide specific guidelines on the construction of the discrete dual domain (Section 3.4); (iv) we provide a MATLAB package for the implementation of the proposed ConjVI algorithm.
Notations. The standard inner product in Rn and the corresponding induced 2-norm are denoted by ⟨·, ·⟩ and ∥·∥2, respectively. ∥·∥∞ denotes the infinity norm. We use the superscript d to denote finite (discrete) sets (as in Xd) and discrete functions (as in hd : Xd → R). We use the superscript g to denote grid-like finite sets (as in Xg = Πni=1X g i where X g i ⊂ R). We also use X g sub to denote the sub-grid of Xg derived by omitting the smallest and the largest elements of Xg in each dimension. The cardinality of the finite set Xd or Xg is denoted by X . We use h̃d : Rn → R = R ∪ {∞} to denote a generic extension of a discrete function hd, and hd : Rn → R to denote multilinear interpolation
and extrapolation (LERP) extension of a discrete function with a grid-like domain. Let X,Y be two arbitrary sets in Rn. co(X) is the convex hull of X. We use d(X,Y) := infx∈X,y∈Y ∥x− y∥2 to denote the distance between X and Y. The one-sided Hausdorff distance from X to Y is defined as dH(X,Y) := supx∈X infy∈Y ∥x− y∥2. Let h : Rn → R be an extended real-valued function with a non-empty effective domain dom(h) = X := {x ∈ Rn : h(x) < ∞}. The range of h is denoted by rng(h) := maxx∈X h(x)−minx∈X h(x), and the subdifferential of h at a point x ∈ X is defined as ∂h(x) := { y ∈ Rn : h(x̃) ≥ h(x) + ⟨y, x̃− x⟩ ,∀x̃ ∈ X } . We define L(h) :=
Πni=1 [ L−i (h),L + i (h) ] , where L+i (h) (resp. L − i (h)) is the maximum (resp. minimum) slope of the function h along the i-th dimension. Note that ∂h(x) ⊆ L(h) for all x ∈ X. We report the complexities using the standard big-O notations O and Õ, where the latter hides the logarithmic factors. In this study, we are mainly concerned with the dependence of the computational complexities on the size of the finite sets involved (discretization of the primal and dual domains). In particular, we ignore the possible dependence of the computational complexities on the dimension of the variables, unless they appear in the power of the size of those discrete sets.
We note that the extended version of this article, including the technical proofs, is available in the supplementary material.

heading: 2 VI in primal domain

text:
We are concerned with the infinite-horizon, discounted cost, optimal control problems of the form
J⋆(x) =min Ewt [ ∞∑ t=0 γtC(xt, ut) ∣∣∣∣x0 = x ]
s.t. xt+1 = g(xt, ut, wt), xt ∈ X, ut ∈ U, wt ∼ P(W), ∀t ∈ {0, 1, . . .}, where xt ∈ Rn, ut ∈ Rm, and wt ∈ Rl are the state, input and disturbance variables at time t, respectively; γ ∈ (0, 1) is the discount factor; C : X × U → R is the stage cost; g : Rn × Rm × Rl → Rn describes the dynamics; X ⊂ Rn and U ⊂ Rm describe the state and input constraints, respectively; and, P(·) is the distribution of the disturbance over the support W ⊂ Rl. Assuming the stage cost C is bounded, the optimal value function solves the equation J⋆ = T J⋆, where T is the DP operator (C and J are extended to infinity outside their effective domains) [7, Prop. 1.2.2]
T J(x) := min u
{ C(x, u) + γ · EwJ ( g(x, u, w) )} , ∀x ∈ X. (1)
Indeed, T is γ-contractive in the infinity-norm [7, Prop. 1.2.4]. This property then gives rise to the VI algorithm Jk+1 = T Jk, which converges to J⋆ as k → ∞, for arbitrary initialization J0. Moreover, assuming that the composition J ◦ g (for each w) and the cost C are jointly convex in the state and input variables, T also preserves convexity [8, Prop. 3.3.1]. For numerical implementation of VI, we need to address three issues. First, we need to compute the expectation in (1). In order to simplify the exposition and include the computational cost of this operation explicitly, we consider disturbances with finite support in this study: Assumption 2.1 (Disturbance with finite support). The disturbance w has a finite support Wd ⊂ Rl with a given probability mass function (p.m.f.) p : Wd → [0, 1].
Under the preceding assumption, we have EwJ ( g(x, u, w) ) = ∑ w∈Wd p(w) · J ( g(x, u, w) ) . The second and more important issue is that the optimization problem (1) is infinite-dimensional for the continuous state space X. This renders the exact implementation of VI impossible, except for a few cases with available closed-form solutions. A common solution to this problem is to deploy a sample-based approach, accompanied by a function approximation scheme. To be precise, for a finite subset Xd of X, at each iteration k = 0, 1, . . ., we take the discrete function Jdk : Xd → R as the input, and compute the discrete function Jdk+1 = [ T J̃dk ]d : Xd → R, where J̃dk : X → R is an extension of Jdk . Finally, for each x ∈ Xd, we have to solve the minimization problem in (1) over the control input. Since this minimization problem is often a difficult, non-convex problem, a common approximation involves enumeration over a discretization Ud ⊂ U. Incorporating these approximations, we end up with the approximate VI algorithm Jdk+1 = T dJdk , characterized by the discrete DP (d-DP) operator
T dJd(x) := min u∈Ud
{ C(x, u) + γ · ∑ w∈Wd p(w) · J̃d ( g(x, u, w) )} , ∀x ∈ Xd. (2)
The convergence of approximate VI described above depends on the properties of the extension operation [̃·]. In particular, if the extension operation is non-expansive (in the infinity-norm), then T d is also γ-contractive. The error of this approximation also depends on the extension operation and its representative power. We refer the interested reader to [7, 10, 26] for detailed discussions on the convergence and error of different approximation schemes for VI.
The d-DP operator and the corresponding VI algorithm will be our benchmark for evaluating the performance of the alternative algorithm developed in this study. To this end, we finish this section with some remarks on the time complexity of the d-DP operation. Let the time complexity of a single evaluation of the extension operator [̃·] in (2) be of O(E). Then, the time complexity of the d-DP operation (2) is of O ( XUWE ) . In this regard, note that the scheme described above essentially involves approximating a continuous-state/action MDP with a finite-state/action MDP, and then applying the VI algorithm. This, in turn, implies the lower bound Ω(XU) for the time complexity (corresponding to enumeration over u ∈ Ud for each x ∈ Xd). This lower bound is also compatible with the best existing time complexities in the literature for VI for finite MDPs; see, e.g., [3, 27]. However, as we will see in the next section, for a particular class of problems, it is possible to exploit the structure of the underlying continuous system in order to achieve a better time complexity in the corresponding discretized problem.

heading: 3 Reducing complexity via conjugate duality

text: In this section, we present the class of problems that allows us to employ conjugate duality and propose an alternative path for solving the corresponding DP operator. We also present the numerical scheme for implementing the proposed alternative path, and analyze its convergence, complexity, and error. We note that the proposed algorithm and its analysis are based on the d-CDP algorithm presented in [20, Sec. 5] for finite-horizon, optimal control of deterministic systems. Here, we extend those results for infinite-horizon, discounted cost, optimal control of stochastic systems. Moreover, unlike [20], our analysis includes the case where the conjugate of input cost is not analytically available and has to be computed numerically; see [20, Assump. 5.1] for more details.

heading: 3.1 VI in conjugate domain

text:
Throughout this section, we assume that the problem data satisfy the following conditions. Assumption 3.1 (Problem class). The problem data has the following properties: (i) The dynamics is of the form g(x, u, w) = f(x, u) + w = fs(x) + Bu + w, with additive disturbance, where fs : Rn → Rn is a Lipschitz continuous, possibly nonlinear map, and B ∈ Rn×m. (ii) The stage cost C is separable in state and input; that is, C(x, u) = Cs(x) + Ci(u), where the state cost Cs : X → R and the input cost Ci : U → R are Lipschitz continuous. (iii) The constraint sets X ⊂ Rn and U ⊂ Rm are compact. Moreover, for each x ∈ X, the set of admissible inputs U(x) := {u ∈ U : g(x, u, w) ∈ X, ∀w ∈ Wd} is nonempty.
Some remarks are in order regarding the preceding assumptions. We first note that the setting of Assumption 3.1 goes beyond the classical LQR. In particular, it includes nonlinear dynamics, state and input constraints, and non-quadratic stage costs. Second, the properties laid out in Assumption 3.1 imply that the set of admissible inputs U(x) is a compact set for each x ∈ X. This, in turn, implies that the optimal value in (1) is achieved if J : X → R is also assumed to be lower semi-continuous. For the problem class of Assumption (3.1), we can use duality theory to present an alternative path for computing the output of the DP operator. This path forms the basis for the algorithm proposed in this study. Fix x ∈ X and consider the following reformulation of the optimization problem (1)
T J(x) =Cs(x) + min u,z {Ci(u) + γ · EwJ(z + w) : z = f(x, u)} ,
where we used additivity of disturbance and separability of stage cost. The corresponding dual problem then reads as
T̂ J(x) := Cs(x) + max y min u,z {Ci(u) + γ · EwJ(z + w) + ⟨y, f(x, u)− z⟩} , (3)
where y ∈ Rn is the dual variable corresponding to the equality constraint. For the deterministic dynamics of Assumption 3.1-(i), we can obtain the following representation for the dual problem.
Proposition 3.2 (CDP operator). The dual problem (3) equivalently reads as
ϵ(x) := γ · EwJ(x+ w), x ∈ X, (4a) ϕ(y) := C∗i (−B⊤y) + ϵ∗(y), y ∈ Rn, (4b) T̂ J(x) = Cs(x) + ϕ∗ ( fs(x) ) , x ∈ X, (4c)
where [·]∗ denotes the conjugate operation.
Following [20], we call the operator T̂ in (4) the conjugate DP (CDP) operator. We next provide an alternative representation of the CDP operator that captures the essence of this operation.
Proposition 3.3 (CDP reformulation). The CDP operator T̂ equivalently reads as
T̂ J(x) = Cs(x) + min u
{ C∗∗i (u) + γ · [EwJ(·+ w)]∗∗ ( f(x, u) )} , (5)
where [·]∗∗ denotes the biconjugate operation.
The preceding result implies that the indirect path through the conjugate domain essentially involves substituting the input cost and (expectation of the) value function by their biconjugates. In particular, it points to a sufficient condition for zero duality gap.
Corollary 3.4 (Equivalence of T and T̂ ). T̂ J = T J if Ci : U → R and J : X → R are convex.
Hence, T̂ has the same properties as T if Ci and J are convex. More importantly, if T and T̂ preserve convexity, then the conjugate VI (ConjVI) algorithm Jk+1 = T̂ Jk, also converges to the optimal value function J⋆, with arbitrary convex initialization J0. For convexity to be preserved, however, we need two more additional assumptions. First, the state cost Cs : X → R needs to be also convex. Then, for T̂ J to be convex, a sufficient condition is the convexity of J ◦ f (jointly in x and u), given that J is convex. The following assumption summarizes the sufficient conditions for equivalence of VI and ConjVI algorithms. Assumption 3.5 (Convexity). Consider the following properties for the constraints, costs, and dynamics: (i) The sets X ⊂ Rn and U ⊂ Rm are convex. (ii) The costs Cs : X → R and Ci : U → R are convex. (iii) The deterministic dynamics f : Rn × Rm → Rn is such that given a convex function J : X → R, the composition J ◦ f is jointly convex in the state and input variables.
We note that the last condition in the preceding assumption usually does not hold for nonlinear dynamics; however, for fs(x) = Ax with A ∈ Rn×n, this is indeed the case for problems satisfying Assumptions 3.1 and 3.5 [6]. Note that, if convexity is not preserved, then the alternative path suffers from duality gap in the sense that in each iteration it uses the convex envelop of (the expectation of) the output of the previous iteration.

heading: 3.2 ConjVI algorithm

text:
The approximate ConjVI algorithm involves consecutive applications of an approximate implementation of the CDP operator (4) until some termination condition is satisfied. Algorithm 1 provides the pseudo-code of this procedure. In particular, we consider solving (4) for a finite set Xd ⊂ X, and terminate the iterations when the difference between two consecutive discrete value functions (in the infinity-norm) is less than a given constant et > 0. Since we are working with a finite subset of the state space, we can restrict the feasibility condition of Assumption 3.1-(iii) to all x ∈ Xd: Assumption 3.6 (Feasibile discretization). We have U(x) ̸= ∅ for all x ∈ Xd.
In what follows, we describe the main steps within the initialization and iterations of Algorithm 1. In particular, the conjugate operations in (4) are handled numerically via the linear-time Legendre transform (LLT) algorithm [23]. LLT is an efficient algorithm for computing the discrete conjugate function over a finite grid-like dual domain. Precisely, to compute the conjugate of the function h : X → R, LLT takes its discretization hd : Xd → R as an input, and outputs hd∗d : Yg → R, for the grid-like dual domain Yg. We refer the reader to [23] for a detailed description of LLT. The main steps of the proposed approximate implementation of the CDP operator (4) are as follows: (i) For the expectation operation in (4a), by Assumption 2.1, we again have EwJ(· + w) =
Algorithm 1 ConjVI: Approximate VI in conjugate domain
Input: dynamics fs : Rn → Rn, B ∈ Rn×m; finite state space Xd ⊂ X; finite input space Ud ⊂ U; state cost function Cds : Xd → R; input cost function Cdi : Ud → R; finite disturbance space Wd and its p.m.f. p : Wd → [0, 1]; discount factor γ; termination bound et. Output: discrete value function Ĵd : Xd → R. initialization:
1: construct the grids Vg, Zg, and Yg; 2: use LLT to compute Cd∗di : Vg → R from Cdi : Ud → R; 3: Jd(x)← 0 and Jd+(x)← Cds (x)−minCdi for x ∈ Xd;
iteration: 4: while ∥∥Jd+ − Jd∥∥∞ ≥ et do 5: Jd ← Jd+; 6: Jd+ ← T̂ dJd : Xd → R according to (6) [Cd∗di in (6c) is already computed in line 2]; 7: end while 8: output Ĵd ← Jd+.
∑ w∈Wd p(w) · J(· + w). Hence, we need to pass the value function Jd : Xd → R through the “scaled expection filter” to obtain εd : Xd → R in (6a) as an approximation of ϵ in (4a). Notice that here we are using an extension J̃d : X → R of Jd (recall that we only have access to the discrete value function Jd). (ii) In order to compute ϕ in (4b), we need access to two conjugate functions. First, for ϵ∗, we use the approximation εd∗d : Yg → R in (6b), by applying LLT to the data points εd : Xd → R for a properly chosen state dual grid Yg ⊂ Rn. We also need the conjugate C∗i of the input cost. If this function is not analytically available, we approximate it as follows: For a properly chosen input dual grid Vg ⊂ Rm, we employ LLT to compute Cd∗di : Vg → R in (6c), using the data points Cdi : Ud → R, where Ud is a finite subset of U. With these conjugate functions at hand, we can now compute φd : Yg → R in (6d), as an approximation of ϕ in (4b). In particular, notice that we use the LERP extension Cd∗di of C d∗d i to approximate C d∗ i at the required point −B⊤y for each y ∈ Yg. (iii) To be able to compute the output according to (4c), we need to perform another conjugate transform. In particular, we need the value of ϕ∗ at fs(x) for x ∈ Xd. Here, we use the approximation φd∗d : Zg → R in (6e), by applying LLT to the data points φd : Yg → R for a properly chosen grid Zg ⊂ Rn. Finally, we use the LERP extension φd∗d of φd∗d to approximate φd∗ at the required point fs(x) for each x ∈ Xd, and compute T̂ dJd in (6f) as an approximation of T̂ J in (4c). With these approximations, we introduce the discrete CDP (d-CDP) operator as follows
εd(x) := γ · ∑
w∈Wd p(w) · J̃d(x+ w), x ∈ Xd, (6a)
εd∗d(y) = max x∈Xd
{ ⟨x, y⟩ − εd(x) } , y ∈ Yg, (6b)
Cd∗di (v) = max u∈Ud
{ ⟨u, v⟩ − Cdi (u) } , v ∈ Vg, (6c)
φd(y) := Cd∗di (−B ⊤y) + εd∗d(y), y ∈ Yg, (6d) φd∗d(z) = max y∈Yg { ⟨y, z⟩ − φd(y) } , z ∈ Zg, (6e)
T̂ dJd(x) := Cs(x) + φd∗d ( fs(x) ) , x ∈ Xd. (6f)
The proper construction of the grids Yg, Vg, and Zg will be discussed in Section 3.4.

heading: 3.3 Analysis of ConjVI algorithm

text:
We now provide our main theoretical results concerning the convergence, complexity, and error of the proposed algorithm. Let us begin with presenting the assumptions to be called in this subsection.
Assumption 3.7 (Grids). Consider the following properties for the grids in Algorithm 1 (consult the Notations in Section 1): (i) The grid Vg is constructed such that co(Vgsub) ⊇ L(Cdi ). (ii) The grid Zg is constructed such that co(Zg) ⊇ fs ( Xd ) . (iii) The construction of Yg, Vg, and Zg requires at most
O(X + U) operations. The cardinality of the grids Yg and Zg (resp. Vg) in each dimension is the same as that of Xd (resp. Ud) in that dimension so that Y,Z = X and V = U .
Assumption 3.8 (Extension operator). Consider the following properties for the operator [̃·] in (6a): (i) The extension operator is non-expansive in the infinity norm; (ii) Given a function J : X → R and its discretization Jd : Xd → R, we have ∥J − J̃d∥∞ ≤ ee for some constant ee ≥ 0.
Our first result concerns the contractiveness of the d-CDP operator.
Theorem 3.9 (Convergence). Let Assumptions 3.7-(ii) and 3.8-(i) hold. Then, the d-CDP operator (6) is γ-contractive in the infinity-norm.
The preceding theorem implies that the approximate ConjVI Algorithm 1 is indeed convergent given that the required conditions are satisfied. In particular, for deterministic dynamics, co(Zg) ⊇ fs ( Xd ) is sufficient for Algorithm 1 to be convergent. We next consider the complexity of our algorithm.
Theorem 3.10 (Complexity). Let Assumption 3.7-(iii) hold. Also assume that each evaluation of the extension operator [̃·] in (6a) requires O(E) operations. Then, the time complexities of initialization and each iteration in Algorithm 1 are of O(X + U) and Õ(XWE), respectively.
The requirements of Assumption 3.7-(iii) will be discussed in Section 3.4. Recall that each iteration of VI (in primal domain) has a complexity of O(XUWE), where E denotes the complexity of the extension operation used in (2). This observation points to a basic characteristic of the proposed approach: ConjVI reduces the quadratic complexity of VI to a linear one by replacing the minimization operation in the primal domain with a simple addition in the conjugate domain. Hence, for problem class of Assumption 3.1, ConjVI is expected to lead to a reduction in the computational cost. We note that ConjVI, like VI and other approximation schemes that utilize discretization/abstraction of the continuous state and input spaces, still suffers from the so-called “curse of dimensionality.” This is because the sizes X and U of the discretizations increase exponentially with the dimensions n and m of the corresponding spaces. However, for ConjVI, this exponential increase is of rate max{m,n}, compared to the rate m+ n for VI.
Theorem 3.11 (Error). Let Assumptions 3.5, 3.7-(i,ii), and 3.8-(i) hold. Consider the true optimal value function J⋆ = T J⋆ : X → R and its discretization Jd⋆ : Xd → R, and let Assumption 3.8-(ii) hold for J⋆. Also, let Ĵd : Xd → R be the output of Algorithm 1. Then, ∥Ĵd − Jd⋆ ∥∞ ≤ γ(ee+et)+ed 1−γ , where ed = eu + ev + ex + ey + ez, and
eu = cu · dH(U,Ud), ev = cv · dH ( co(Vg),Vg ) , ex = cx · dH ( X,Xd ) ,
ey = cy ·maxx∈Xd d ( ∂(J⋆ − Cs)(x),Yg ) , ez = cz · dH ( fs(Xd),Zg ) ,
(7)
with constants cu, cv, cx, cy, cz > 0 depending on the problem data.
Let us first note that Assumption 3.5 implies that the DP and CDP operators preserve convexity, and they both have the true optimal value function J⋆ as their fixed point (i.e., the duality gap is zero). Otherwise, the proposed scheme can suffer from large errors due to dualization. Moreover, Assumptions 3.7-(i,ii) on the grids Vg and Zg are required for bounding the error of approximate discrete conjugations using LERP in (6d) and (6f); see the proof of Lemmas A.5 and A.7 of the supplementary material. The remaining sources of error in the proposed approximate implementation of ConjVI are captured by the three main error terms: (i) ee is due to the approximation of the value function using the extension operator [̃·]; (ii) et corresponds to the termination of the algorithm after a finite number of iterations; (iii) ed captures the error due to the discretization of the primal and dual state and input domains.

heading: 3.4 Construction of the grids

text:
In this subsection, we provide specific guidelines for the construction of the grids Yg, Vg and Zg. We note that these discrete sets must be grid-like since they form the dual grid for the three conjugate transforms that are handled using LLT. The presented guidelines aim to minimize the error terms in (7) while taking into account the properties laid out in Assumption 3.7. In particular, the schemes described below satisfy the requirements of Assumption 3.7-(iii).
Construction of Vg. Assumption 3.7-(i) and the error term ev in (7) suggest that we find the smallest input dual grid Vg such that co(Vgsub) ⊇ L(Cdi ). This latter condition essentially means that Vg must “more than cover the range of slope” of the function Cdi ; recall that L(Cdi ) = Πmj=1 [ L−j (C d i ),L − j (C d i ) ] , where L−j (C d i ) (resp. L + j (C d i )) is the minimum (resp. maximum) slope of Cdi along the j-th dimension. Hence, we need to compute/approximate L ± j (C d i ) for j = 1, . . . ,m. A conservative approximation is L−j (Ci) = min ∂Ci/∂uj and L + j (Ci) = max ∂Ci/∂uj , assuming Ci is differentiable. Alternatively, we can directly use the discrete input cost Cdi for computing L±j (C d i ). In particular, if the domain Ud = Ug = Πmj=1U g j of C d i is grid-like and Ci is convex, we can take L−j (C d i ) (resp. L + j (C d i )) to be the minimum first forward difference (resp. maximum last backward difference) of Cdi along the j-th dimension (this scheme requires O(U) operations). Having L±j (C d i ) at our disposal, we can then construct V g sub = Π m j=1V g subj
such that, in each dimension j, Vgsubj is uniform and has the same cardinality as U g j , and co(V g subj ) = [ L−j (C d i ),L + j (C d i ) ] . Finally, we construct Vg by extending Vgsub uniformly in each dimension (by adding a smaller and a larger element to Vgsub in each dimension). Construction of Zg. According to Assumption 3.7-(ii), the grid Zg must be constructed such that co(Zg) ⊇ fs ( Xd ) . This can be simply done by finding the vertices of the smallest box that contains
the set fs ( Xd ) . Those vertices give the diameter of Zg in each dimension. We can then, for example, take Zg to be the uniform grid with the same cardinality as Yg in each dimension (so that Z = Y ). This way, dH ( fs(Xd),Zg ) ≤ dH ( co(Zg),Zg ) , and hence ez in (7) reduces by using finer grids Zg. This construction has a time complexity of O(X). Construction of Yg. Construction of the state dual grid Yg is more involved. According to Theorem 3.11, we need to choose a grid that minimizes ey in (7). This can be done by choosing Yg such that Yg ∩ ∂(J⋆ − Cs) ̸= ∅ for all x ∈ Xd so that ey = 0. Even if we had access to the optimal value function J⋆, satisfying such a condition could lead to a dual grid Yg ⊂ Rn of size O(Xn). Such a large size violates Assumption 3.7-(iii) on the size of Yg, and essentially renders the proposed algorithm impractical for dimensions n ≥ 2. A more practical condition is co(Yg)∩ ∂(J⋆ −Cs) ̸= ∅ for all x ∈ Xd so that maxx∈Xd d ( ∂(J⋆ − Cs)(x),Yg ) ≤ dH ( co(Yg),Yg ) , and hence ey reduces by using a finer grid Yg. The latter condition is satisfied if co(Yg) ⊇ L(J⋆ − Cs), i.e., if Yg “covers the range of slops” of (J⋆−Cs). Hence, we need to approximate the range of slopes of (J⋆−Cs). To this end, we first use the fact that J⋆ is the fixed point of DP operator (1) to approximate rng(J⋆−Cs) by R = rng(C d i )+γ·rng(C d s )
1−γ . We then construct the gird Y g = Πni=1Y g i such that, for each dimension i,
we have ±αR/∆iXd ∈ co(Y g i ), where ∆ i Xd denotes the diameter of the projection of X d on the i-th dimension, and α > 0 is a scaling factor mainly depending on the dimension of the state space. This construction has a one-time computational cost of O(X + U). Dynamic construction of Yg. Alternatively, we can construct Yg dynamically at each iteration in order to minimize the corresponding error in each application of the d-CDP operator given by ey = cy ·maxx∈Xd d ( ∂(T J−Cs)(x),Yg ) ; see Lemma A.6 and Proposition A.8 of the supplementary material. This means that the construction of Yg in Algorithm 1 is moved from line 1 to inside the iterations, after line 5. Similar to the static scheme described above, the aim here is to construct Yg such that co(Yg) ⊇ L(T J − Cs). Since we do not have access to T J (it is the output of the current iteration), we can again use the definition of the DP operator (1) to approximate rng(T J − Cs) by R = rng(Cdi ) + γ · rng(Jd) where Jd is the output of the previous iteration. We then construct the gird Yg = Πni=1Y g i such that, for each dimension i, we again have ±αR/∆iXd ∈ co(Y g i ). This construction has a one-time computational cost of O(U) for computing rng(Cdi ) and a per iteration computational cost of O(X) for computing rng(Jd). Notice, however, that under this dynamic construction, the error bound of Theorem 3.11 does not hold true. More importantly, with a dynamic grid Yg, there is no guarantee for ConjVI to converge.

heading: 4 Numerical simulations

text:
We now showcase the application of the ConjVI algorithm through two numerical examples, and compare its performance with the VI algorithm. The details of these simulations (and another numerical example) are provided in Section 4 of the supplementary material. We also provide the
ConjVI MATLAB package [21] for the implementation of the proposed algorithm. The package also includes the numerical simulations of this section. We note that multiple routines in the developed package are borrowed from the d-CDP MATLAB package [22]. Also, for the discrete conjugation (LLT), we used the MATLAB package (in particular, the LLTd routine) provided in [23].
Example 1. We use the setup of [20, App. C.2.2] for the optimal control of a noisy inverted pendulum with two state variables and one input channel. In particular, we use nearest neighbor extension (which is non-expansive) for the extension operators in (2) for VI and in (6a) for ConjVI. We also consider the dynamic scheme for the construction of Yg in ConjVI (hereafter, referred to as ConjVI-d). Moreover, in each implementation of VI and ConjVI(-d), all of the involved discrete sets are uniform grids with the same size N in each dimension, i.e., X,Y, Z = N2 and U, V = N . We are particularly interested in the performance of these algorithms, as the size of the discretizations increases. The results of our numerical simulations are shown in the top panels of Figure 1. As shown in Figures 1a, both VI and ConjVI algorithms are convergent with a rate ≤ γ = 0.95. In particular, ConjVI terminates in 57 iterations, compared to 101 iterations required for VI to terminate. As expected, this faster convergence, combined with the lower time complexity of ConjVI in each iteration, leads to a significant reduction in the running time of this algorithm compared to VI; see Figure 1b. Since we do not have access to the true optimal value function, in order to evaluate the performance of the outputs of the these algorithms, we consider the performance of the greedy policy with respect to the discrete value function Jd computed using these algorithms. Figure 1c reports the average cost of 100 random instances over T = 100 steps. As shown, the reduction in the running time in ConjVI comes with an increase in the cost of the controlled trajectories, particularly for coarse discretizations of the state-action space. On the other hand, using the dynamic scheme for construction of Yg, we see that ConjVI-d, with a small increase in the running time compared to ConjVI, achieves almost the same performance as VI concerning the quality of the greedy actions.
Example 2. We now consider the optimal control of an unstable deterministic batch reactor with four states and two inputs. The setup is borrowed from [19, Sec. 6]. Once again, all of the involved discrete sets are uniform grids with the same size N in each dimension, i.e., X,Y, Z = N4 and U, V = N2. We note that we use multi-linear interpolation and extrapolation for the extension operator in (2) for VI. Due to the extrapolation, the extension operator is no longer non-expansive, and hence the convergence of VI is not guaranteed. On the other hand, since the dynamics is deterministic, there is no need for extension in ConjVI (the scaled expectation in (6a) in ConjVI reduces to the simple scaling εd = γ · Jd for deterministic dynamics), and hence the convergence of ConjVI only requires
co(Zg) ⊇ fs ( Xg ) . The results of simulations are shown in the bottom panels of Figure 1. Once again, we see the trade-off between the time complexity and the greedy control performance in VI and ConjVI. On the other hand, ConjVI-d has the same control performance as VI with an insignificant increase in running time compared to ConjVI. Moreover, in Figure 1d, we observe the non-monotone behavior of ConjVI-d. In this regard, recall that when the grid Yg is constructed dynamically and varies at each iteration, the d-CDP operator is not necessarily contractive, which is here the case for the first six iterations. The VI algorithm is also showing a non-monotone behavior, where for the first nine iterations the d-DP operation is actually expansive. As we noted earlier, this is because the multi-linear extrapolation operation used for extension is expansive.

heading: 5 Final remarks

text:
In this paper, we proposed the ConjVI algorithm which reduces the time complexity of the VI algorithm from O(XU) to O(X +U). This better time complexity however comes at the expense of restricting the class of problem. In particular, there are two main conditions that must be satisfied in order to be able to apply the ConjVI algorithm: The dynamics must be of the form x+ = fs(x) +Bu+ w, and the stage cost C(x, u) = Cs(x) + Ci(u) must be separable. Moreover, since ConjVI essentially solves the dual problem, for non-convex problems, it suffers from a non-zero duality gap. Based on our simulation results, we also notice a trade-off between computational complexity and control action quality: While ConjVI has a lower computational cost, VI generates better control actions. However, the dynamic scheme for the construction of state dual grid Yg allows us to achieve almost the same performance as VI when it comes to the quality of control actions, with a small extra computational burden. In what follows, we provide our final remarks on the limitations of the proposed ConjVI algorithm and its relation to existing approximate VI algorithms.
Relation to existing approximate VI algorithms. The basic idea for complexity reduction introduced in this study can be potentially combined with and further improve the existing sample-based VI algorithms. These sample-based algorithms solely focus on transforming the infinite-dimensional optimization in DP problems into computationally tractable ones, and in general, they have a time complexity of O(XU), depending on the product of the cardinalities of the discrete state and action spaces. The proposed ConjVI algorithm, on the other hand, focuses on reducing this time complexity to O(X + U), by avoiding the minimization over input in each iteration. Take, for example, the aggregation technique in [26, Sec. 8.1] that leads to a piece-wise constant approximation of the value function. It is straightforward to combine ConjVI with this type of state space aggregation. Indeed, the first numerical example of Section 4 (optimal control of inverted pendulum) essentially uses such an aggregation by approximating the value function via nearest neighbor extension.
Cost functions with a large Lipschitz constant. Recall that for the proposed ConjVI algorithm to be computationally efficient, the size Y of the state dual grid Yg must be controlled by the size X of the discrete state space Xd (Assumption 3.7-(iii)). Then, as the range of slope of the value function J⋆ increases, the corresponding error ey in (7) due to discretization of the dual state space increases. The proposed dynamic approach for construction of Yg partially addresses this issue by focusing on the range of slope of Jdk in each iteration in order to minimize the discretization error of the same iteration k. However, when the cost function has a large Lipschitz constant, even this latter approach can fail to provide a good approximation of the value function. (See the supplementary material for a numerical example).
Gradient-based algorithms for solving the minimization over input. Let us first note that the minimization over u in sample-based VI algorithms usually involves solving a difficult non-convex problem. This is particularly due to that fact that the extension operation employed in these algorithms for approximating the value function using the sample points does not lead to a convex function in u (e.g., take kernel-based approximations or neural networks). This is why in MDP and RL literature, it is actually quite common to consider a finite action space in the first place [10, 26]. Moreover, the minimization over u again must be solved for each sample point in each iteration, while application of ConjVI avoids solving this minimization in each iteration. In this regard, let us note that ConjVI actually uses a convex approximation of the value function, which allows for application of a gradient-based algorithm for minimization over u within the ConjVI algorithm. However, such an algorithm has a per iteration complexity of O(XY ) = O(X2) , which is practically inefficient.

heading: Acknowledgments

text: This research is part of a project that has received funding from the European Research Council (ERC) under the grant TRUST-949796. The authors are also grateful to anonymous reviewers for their comments concerning the three remarks in Section 5.

title: Fast Approximate Dynamic Programming for Infinite-Horizon Markov Decision Processes