#
Uniqueness and Complexity of Inverse MDP Models

Marcus Hutter and Steven Hansen
DeepMind
http://www.hutter1.net/

(2 June 2022)

###### Abstract

What action sequence $aa^{\prime}a^{\prime\prime}$ was likely responsible for
reaching state $s^{\prime\prime\prime}$ (from state $s$) in 3 steps?
Addressing such questions is important in causal reasoning and in
reinforcement learning. Inverse “MDP” models
$p(aa^{\prime}a^{\prime\prime}|ss^{\prime\prime\prime})$ can be used to answer
them. In the traditional “forward” view, transition “matrix”
$p(s^{\prime}|sa)$ and policy $π(a|s)$ uniquely determine “everything”: the
whole dynamics
$p(as^{\prime}a^{\prime}s^{\prime\prime}a^{\prime\prime}...|s)$, and with it,
the action-conditional state process
$p(s^{\prime}s^{\prime\prime}...|saa^{\prime}a^{\prime\prime})$, the multi-
step inverse models $p(aa^{\prime}a^{\prime\prime}...|ss^{i})$, etc. If the
latter is our primary concern, a natural question, analogous to the forward
case is to which extent 1-step inverse model $p(a|ss^{\prime})$ plus policy
$π(a|s)$ determine the multi-step inverse models or even the whole dynamics.
In other words, can forward models be inferred from inverse models or even be
side-stepped. This work addresses this question and variations thereof, and
also whether there are efficient decision/inference algorithms for this.

###### Contents

  1. 1 Introduction
  2. 2 Problem Formalization and Preliminaries
  3. 3 Degenerative Cases
  4. 4 (Non)Uniqueness of Inverse MDP Models
  5. 5 Linear Relaxation
  6. 6 Experiments
  7. 7 Computational Complexity
  8. 8 Conclusion
  9. A List of Notation
  10. B Application to Planning
  11. C Counter-Examples in Related Work
  12. D Deterministic Cases
  13. E Characterizing $M$ and $W$ for which EqIM(1) holds
  14. F Characterizing $M$ and $W$ for which EqIM(1) and EqIM(2+) hold
  15. G EqIM(1)$∧$EqIM(2+)$\not\rightarrow$EqIM(3) for full low rank $M$?
  16. H EqIM(1) does not imply EqIM(2) ($⊙$-version)
  17. I $\text{EqIM}(1a)∧...∧\text{EqIM}(ia)$ do not imply EqIM($i+1$) ($⊙$-version)
  18. J Non-Uniqueness of Inverse MDP Models for $i≥2$
  19. K Systems of Quadratic Matrix Equations
  20. L Compact Representation of EqIM(2+)
  21. M Open Problem
  22. N Experimental Details
  23. O How to Deal with 0/0
  24. P Solution Dimensions of $W$ and $B^{aa^{\prime}}$.

Keywords

> inverse models; reinforcement learning; causality; theory; multi-step
> models; planning.

##  1 Introduction

Consider an MDP with actions $a∈\\{0,..,k-1\\}$ and states $s∈\\{1,...,d\\}$.
Rewards play no role in our analysis, so _controlled Markov process_ [DY79] or
_conditional Markov chain_ may be a more apt naming. Transition “matrix”
$p(s^{\prime}|sa)$ (_“Forward model”_) and policy $π(a|s)$ uniquely determine
the whole _dynamics_

$\displaystyle
p(as^{\prime}a^{\prime}s^{\prime\prime}a^{\prime\prime}...|s)~{}=~{}π(a|s)⋅p(s^{\prime}|sa)⋅π(a^{\prime}|s^{\prime})⋅p(s^{\prime\prime}|s^{\prime}a^{\prime})⋅...$
(1)

and also determines the action-conditional state process (“Multi-Step Forward
Model”):

$\displaystyle
p(s^{\prime}s^{\prime\prime}...|saa^{\prime}a^{\prime\prime})~{}=~{}p(as^{\prime}a^{\prime}s^{\prime\prime}a^{\prime\prime}...|s)/∑_{s^{\prime}s^{\prime\prime}...}p(as^{\prime}a^{\prime}s^{\prime\prime}a^{\prime\prime}...|s)$
(2)

Here we consider _Inverse Model_ $p(a|ss^{\prime})$ and _Multi-Step Inverse
Models_
$p(aa^{\prime}a^{\prime\prime}...|ss^{\prime}s^{\prime\prime}s^{\prime\prime\prime}...)$
and $p(a|ss^{i})$ and variations thereof. Inverse MDP models should not be
confused with inverse reinforcement learning [AD21], which infers rewards,
which play no role here.

### Motivation.

One motivation to consider inverse models is causal inference: An inverse
model captures the likelihood that an action $a$ was the cause of the
transition from state $s$ to state $s^{\prime}$. A multi-step inverse model
captures the likelihood that a first action $a$ or action sequence
$aa^{\prime}...a^{i-1}$ was the cause of the state sequence
$ss^{\prime}...s^{i}$ or the cause of the transition from state $s$ to state
$s^{i}$. The latter is the primary goal in (automatic/stochastic) planning
[HSHB99]: to find an action sequence that leads to a desired goal state
$s^{i}=s_{\text{goal}}$. The shortest path, i.e. smallest $i$, that reaches
$s_{\text{goal}}$ (with high probability in the stochastic case) can easily be
found via a trivial search over $i=1,2,3,...$ if the fixed-$i$ planning
problem can be solved efficiently.

Another machine-learning motivation is that inverse models may be
substantially smaller than forward models. For instance, an action-independent
Markov process $p(s^{\prime}|sa)=p(s^{\prime}|s)$ may be very complex for
large $d$, but for a state-independent (known) policy $π(a|s)=π(a)$, the
inverse model $p(aa^{\prime}...|s..s^{\prime\prime}..)=π(a)π(a^{\prime})...$
is trivial (and known). Of course this extreme case is uninteresting, but a
partial similar simplification happens if state $s$ decomposes into
$s=(\dot{s},\ddot{s})$ [EMK+22]. In this case, if the forward model
$p(s^{\prime}|sa)$ factors into a (simple) controlled
$p(\dot{s}^{\prime}|\dot{s}a)$ and (complex) uncontrolled
$p(\ddot{s}^{\prime}|\ddot{s})$, and the policy $π(a|s)=π(a|\dot{s})$ only
depends on (small) $\dot{s}$, then
$p(aa^{\prime}...|s..s^{\prime\prime}..)=p(aa^{\prime}...|\dot{s}..\dot{s}^{\prime\prime}..)$
is independent of (large) $\ddot{s}$. Note that this simplification happens
“automatically”. We do not need to know the factorization structure, say
$(\dot{s},\ddot{s})=f(s)$ for some unknown $f$. Appendix B contains a bit of
practical context/motivation/application.

### Main questions.

The main question we consider here is:
to which extent do inverse model $p(a|ss^{\prime})$ plus policy $π(a|s)$
determine the multi-step inverse model or even the whole dynamics.

For instance, do $p(a|ss^{\prime})$ plus $π(a|s)$ determine

  * (i)

the full dynamics (1),

  * (ii)

the full dynamics, if also $p(aa^{\prime}|ss^{\prime\prime})$ is provided,

  * (iii)

the multi-step inverse model $p(aa^{\prime}...|ss^{i})$ (or
$p(aa^{\prime}...|ss^{\prime}s^{\prime\prime}...)$),

  * (iv)

the multi-step inverse model $p(aa^{\prime}...|ss^{i})$ (or
$p(aa^{\prime}...|ss^{\prime}s^{\prime\prime})$), if also
$p(a|ss^{\prime\prime})$ is provided,

  * (v)

just the initial action $p(a|ss^{\prime\prime})$ from just final state
$s^{\prime\prime}$,

  * (vi)

$p(a|ss^{i})$ if also $p(a|ss^{\prime\prime})$ is provided,

and variations thereof? Also, is there an efficient algorithm that can decide
whether the solution is unique and/or computes any or all of them?

Unlike in the “forward” case (1), the answer to all these questions is
‘complicated’ and ‘sometimes’. For instance, (i) is true iff $k≥d$ and
$p(s^{\prime}|sa)$ has full rank. (ii) seems true for “most” transition
matrices. (iii-vi) can fail, but (iv) and (vi) seem to hold for interesting
cases. In some situations there are efficient algorithms which sometimes work.

### (Un)Related work.

There is of course abundant literature on causal reasoning in general [PGJ16],
and in the modern context of Deep Learning in particular [OKD+21], but to the
best of our knowledge, the setup and questions we are asking are novel, at
least in this generality and rigor.

A special case of our setup is considered in [EMK+22]. The authors consider
Exogenous Block MDPs (EX-BMDPs) which correspond to the motivating
decomposition example above, and formalized in Section 3 as tensor-product
MDPs. Additionally they assume episodic MDPs with near-deterministic dynamics.
On the other hand, they allow for block-observation functions, which we don’t.
This doesn’t increase the model class, but can lead to smaller controlled-
state representations. Their PPE algorithm finds action sequences of high
inverse probability $p(aa^{\prime}...a^{i-1}|ss^{i})$ in polynomial time in
$\dot{s}$ rather than $s$, while our aim is to infer higher- from lower-step
inverse models for general MDPs.

In the context of Deep Learning, there is ample empirical work that would
benefit from a positive answer to our main question: Variational Intrinsic
Control [GRW17] and Diversity is All You Need [EGIL18] are representative of a
broad class of methods that learn diverse options (policies / action
sequences) that are inferrable from their effects on the environment. This
relies on inverse modelling, as their mutual information objective is
decomposed into maximizing skill/policy entropy and minimizing the entropy of
an inverse model:

$I(s^{i};a...a^{i-1}|s)\equiv
H(a...a^{i-1}|s)–H(aa^{\prime}...a^{i-1}|ss^{i})$

This is akin to finding all action sequences of sufficiently high probability
$p(aa^{\prime}...a^{i-1}|ss^{i})$, or all skills when the policy space is
captured by an auxiliary variable $p(z|ss^{i})$. The EDDICT algorithm [HDB+21]
also maximizes this objective, and parameterizes the requisite inverse models
such that they yield forward predictions, but as detailed in Section 4 its
unlikely that such models would yield optimal multi-step inverse predictions
in general. Dynamics-Aware Unsupervised Discovery of Skills [SGL+19]
decomposes the mutual information in the opposite direction, so as to avoid
learning an inverse model and instead relies on a conventional forward model.
Uniting all of the above mentioned methods is that the action sequence/skill
horizon $i$ must be fixed a priori. Inferring long horizon inverse models from
shorter ones (the topic of the present work) would allow all of these methods
to circumvent this constraint.

A second stream of empirical work uses single-step inverse models for
representation learning [BEP+18]. Agent57 is arguably the most prominent of
these methods [BPK+20], and therein the authors note that this choice of
representation limits the generality of their approach, as multi-step effects
can be aliased over. Despite this being a known limitation, multi-step inverse
models are not used as they are too cumbersome to effectively learn online. A
positive result to our questions (iii) or (iv) would allow such methods to
leverage multi-step inverse predictions despite only learning a single-step
model.

These two beneficiaries of improvements to the construction of multi-step
inverse models (filtering action sequences and state abstraction) dovetail
into potential benefits for a broad range of planning algorithms. Exploiting
this relationship between the questions addressed here and planning problems
is left to future work, but we sketch out the motivation more fully in Section
B.

Regarding the computational complexity of our problem: There is abundant
literature on approximately/heuristically solving the NP-complete class of
systems of polynomial/quadratic equations [Stu02], but we could not find an
NP-complete sub-class of ours [HL13]. The (NP-complete) problem of inference
in Bayesian networks [KF09] also is of different nature. We also don’t have a
learning problem to deal with. We do have a latent variable / partial
observability problem ($s^{\prime}...s^{i-1}$ integrated out), but we are not
(directly) interested in inferring them.

### Contents.

In Section 2 we will formalize questions (i)-(vi) in matrix/tensor notation.
Section 3 gives a first probe into these questions by considering various
degenerate cases. In Section 4 we study the solvability and uniqueness
questions (i),(iii),(v), when only $B^{a}$ is given, i.e. the case $i=1$, in
preparation for and showing the necessity of considering $i>1$. In Section 5
we provide a polynomial-time algorithm via linear relaxation that works under
certain conditions. Section 6 provides some validation experiments on toy
domains. Section 7 looks into the possible NP-hardness of this problem, even
setting uniqueness aside. Section 8 concludes, followed by references.

Appendix A contains a list of notation. The other appendices contain further
considerations: a brief discussion of some applications related to the use of
inverse models in planning (B), counter-examples in related work (C),
deterministic cases (D), an alternative derivation of characterizing when (i)
holds (E), an explicit representation as a System of Quadratic Equations (SQE)
(F) with empirical rank analysis (G), counter-examples to (ii),(iv),(vi), and
indeed for all $i>1$, which is surprising given the severe over-determined
nature of the problem (H,I,J), general matrix SQEs (K), a simple and self-
contained instantiation (L) of the main open problem (M), experimental details
and supplementary figures (N), proper handling of 0/0 (O), and formulas for
the dimension of the solution spaces in case they are not unique (P).

##  2 Problem Formalization and Preliminaries

We now formalize our questions (i)-(vi) from the introduction, and for this
purpose introduce some useful matrix notation. We are not aware of prior work
addressing these questions, so quite some ground-work to suitably formalize
the various question is needed, and many little results are derived or
mentioned in passing to give better insight into the structure of the problem.
To avoid clutter, we will not constantly point out edge cases or domain
constraints. For instance quantities that represent probabilities are
obviously non-negative and sum to one. The reader worried about divisions by
$0$ here and there should best assume that all probabilities are strictly
positive, but most considerations and results naturally generalize with some
care, e.g. by adding “almost surely” w.r.t. to the joint distribution (1).
Appendix O contains a proper treatment of 0/0.

### Notation.

Capital letters $B,D,I,M,W,...$ are used for $d×d$ matrices over $[0;1]⊂ℝ$ and
tensors by adding further upper indices, e.g. $M^{⋅}_{⋅⋅}$ is an order-3
tensor, and $M^{a}_{⋅⋅}$ a matrix for each $a∈\\{0:k-1\\}:=\\{0,...,k-1\\}$,
and $A,C,V,...$ are other tensors. We define Id to be the identity (eye)
matrix
$\text{Id}_{ss^{\prime}}:=δ_{ss^{\prime}}:=[\mskip-3.0mu[s=s^{\prime}]\mskip-3.0mu]~{}∀s,s^{\prime}∈\\{1:d\\}$,
and $I$ to be the all-one matrix $I_{ss^{\prime}}=1~{}∀ss^{\prime}$. We drop
all-quantifiers $∀s,s^{\prime},...$ if clear from context. Let $⊙$ denote
element-wise (Hadamard) multiplication
($[A⊙B]_{ss^{\prime}}=A_{ss^{\prime}}B_{ss^{\prime}}$), and similarly $⊘$,
while (no) $⋅$ represents (conventional) matrix multiplication and has
operator preference over $⊙$ and $⊘$. Matrices form a ring under conventional
$(+,⋅)$ and a commutative ring under $(+,⊙)$, but $(A⋅B)⊙C≠A⋅(B⊙C)$. A
diagonal matrix $D$ has the property $D=D⊙\text{Id}$, i.e.
$D_{ss^{\prime}}=D_{ss}[\mskip-3.0mu[s=s^{\prime}]\mskip-3.0mu]$. $V:=I⋅D$ is
a matrix with $D_{ss}$ in the whole of column $s$
($V_{ss^{\prime}}=V_{*s^{\prime}}=D_{s^{\prime}s^{\prime}}$). Note that
$A⋅D=A⊙V$
($[A⋅D]_{ss^{\prime\prime}}=∑_{s^{\prime}}A_{ss^{\prime}}D_{s^{\prime}s^{\prime\prime}}=A_{ss^{\prime\prime}}D_{s^{\prime\prime}s^{\prime\prime}}=A_{ss^{\prime\prime}}V_{*s^{\prime\prime}}=[A⊙V]_{ss^{\prime\prime}}$).
Similar left-right reversed identities hold. $\bot$ denotes ‘undefined’. See
Appendix A for a full List of Notation.

### Matrix/tensor formalization.

We define

$\displaystyle
M^{a}_{ss^{\prime}}~{}:=~{}p(as^{\prime}|s)~{}=~{}π(a|s)p(s^{\prime}|sa)$

Marginalizing out the action, gives

$\displaystyle
p(s^{\prime}|s)~{}=~{}∑_{a}p(as^{\prime}|s)~{}=~{}∑_{a}M^{a}_{ss^{\prime}}~{}=:~{}M^{+}_{ss^{\prime}}$

Marginalizing out the next-state, gives back

$\displaystyle
π(a|s)~{}=~{}∑_{s^{\prime}}p(as^{\prime}|s)~{}=~{}∑_{s^{\prime}}M^{a}_{ss^{\prime}}~{}=:~{}M^{a}_{s+}$

For instance, the multi-step dynamics can be written as

$\displaystyle
p(as^{\prime}a^{\prime}s^{\prime\prime}...|s)~{}=~{}p(as^{\prime}|s)⋅p(a^{\prime}s^{\prime\prime}|a^{\prime})⋅...~{}=~{}M^{a}_{ss^{\prime}}M^{a^{\prime}}_{s^{\prime}s^{\prime\prime}}...$

Marginalizing out the intermediate states gives

$\displaystyle
p(aa^{\prime}...a^{i-1}s^{i}|s)~{}=~{}[M^{a}⋅M^{a^{\prime}}...⋅M^{a^{i-1}}]_{ss^{i}}$

The inverse MDP model can then be expressed as

$\displaystyle
B^{a}_{ss^{\prime}}~{}:=~{}p(a|ss^{\prime})~{}=~{}p(as^{\prime}|s)/p(s^{\prime}|s)~{}=~{}M^{a}_{ss^{\prime}}/M^{+}_{ss^{\prime}}~{}=~{}[M^{a}⊘M^{+}]_{ss^{\prime}}$

The multi-step inverse model given the whole state sequence becomes

$\displaystyle
p(aa^{\prime}...|ss^{\prime}s^{\prime\prime}...)~{}=~{}\frac{p(as^{\prime}|s)p(a^{\prime}s^{\prime\prime}|s^{\prime})...}{p(s^{\prime}|s)p(s^{\prime\prime}|s^{\prime})...}~{}=~{}\frac{M^{a}_{ss^{\prime}}M^{a^{\prime}}_{s^{\prime}s^{\prime\prime}}...}{M^{+}_{ss^{\prime}}M^{+}_{s^{\prime}s^{\prime\prime}}...}~{}=~{}p(a|ss^{\prime})p(a^{\prime}|s^{\prime}s^{\prime\prime})...$
(3)

and can easily be computed from the 1-step inverse models. To answer the
primary question: which action sequence can lead to (desired) state $s^{i}$
from state $s$, we need to marginalize out $s^{\prime}...s^{i-1}$. For
instance, the two-step inverse model from $s$ to $s^{\prime\prime}$ with
$s^{\prime}$ marginalized out becomes

$\displaystyle
B^{aa^{\prime}}_{ss^{\prime\prime}}~{}:=~{}p(aa^{\prime}|ss^{\prime\prime})~{}=~{}\frac{∑_{s^{\prime}}M^{a}_{ss^{\prime}}M^{a^{\prime}}_{s^{\prime}s^{\prime\prime}}}{∑_{s^{\prime}}M^{+}_{ss^{\prime}}M^{+}_{s^{\prime}s^{\prime\prime}}}~{}=~{}[M^{a}⋅M^{a^{\prime}}⊘(M^{+})^{2}]_{ss^{\prime\prime}}$
(4)

Note that unlike the forward case, $B^{aa^{\prime}}≠B^{a}⋅B^{a^{\prime}}$,
which is responsible for all the problems we will face. Also $B^{a+}≠B^{a}$
but $B^{+}=1=B^{++}$. We always use brackets to denote and disambiguate
(matrix) powers $()^{2}$ from upper indices $M^{a}$. The initial-action 2-step
(and similarly $i$-step) inverse models follow from further marginalizing
$a^{\prime}a^{\prime\prime}...$:

$\displaystyle B^{a+}_{ss^{\prime\prime}}~{}$
$\displaystyle=~{}p(a|ss^{\prime\prime})~{}=~{}[M^{a}M^{+}⊘(M^{+})^{2}]_{ss^{\prime\prime}},$
$\displaystyle B^{a+^{i-1}}_{ss^{i}}~{}$
$\displaystyle=~{}p(a|ss^{i})~{}=~{}[M^{a}(M^{+})^{i-1}⊘(M^{+})^{i}]_{ss^{i}}$
(5)

With this notation, questions (i-vi) in the introduction can formally be
written as

  * (i)

Can $M$ be inferred from $B^{a}:=M^{a}⊘M^{+}$?

  * (ii)

Can $M$ be inferred from $B^{a}$ and
$B^{aa^{\prime}}:=M^{a}M^{a^{\prime}}⊘(M^{+})^{2}$?

  * (iii)

Can $B^{aa^{\prime}...a^{i}}:=M^{a}M^{a^{\prime}}...M^{a^{i}}⊘(M^{+})^{i}$ be
inferred from $B^{a}$?

  * (iv)

Can $B^{aa^{\prime}...a^{i}}$ be inferred from $B^{a}$ and $B^{aa^{\prime}}$?

  * (v)

Can $B^{a+}:=M^{a}M^{+}⊘(M^{+})^{2}$ be inferred from $B^{a}$?

  * (vi)

Can $B^{a++}:=M^{a}(M^{+})^{2}⊘(M^{+})^{3}$ be inferred from $B^{a}$ and
$B^{a+}$?

Each question comes in two versions, given also $π$, or not knowing $π$. We
mainly consider the former version, i.e. knowing $M^{a}_{s+}$:

$\displaystyle\text{Constraint on $M$ for known
$π$:}~{}~{}~{}M^{a}_{s+}~{}=~{}π(a|s)~{}~{}~{}\text{and in
particular}~{}~{}~{}M^{+}_{s+}~{}=~{}1$ (6)

Questions (i)-(vi) also have multiple variations:

  * (I)

Assume some arbitrary $B^{a}$ (and $B^{aa^{\prime}}$) is given, but not
defined via $M$.
Is there no, exactly one, or multiple $M$ consistent with these $B$?

  * (II)

Is there an efficient algorithm that can decide the previous question?

  * (III)

Is there an efficient algorithm that can compute any/all solutions if one/many
exist, and halts/loops if not (4 non-trivial combinations of $/$).

  * (IV)

Can we efficiently determine the “number” of solutions,
e.g. the dimension of the variety formed by the set of all solutions.

### Formulation of the uniqueness questions.

Abstractly, these questions ask whether $M$ (in case of (i-ii)) or $g(M)$ for
some function $g$ (in case of (iii-vi)) can be inferred from some other
function $f(M)$. Let us define another MDP $q(s^{\prime}|sa)$ with same policy
$π(s|a)$ and shorthand

$\displaystyle W_{ss^{\prime}}^{a}~{}:=~{}π(a|s)q(s^{\prime}|sa)$

(In applications, $B^{a}$ would be learned from data, and $W$ or
$B^{aa^{\prime}...}$ inferred from $B^{a}$ in the hope that $W≈M$.) One way to
rephrase the questions is whether $f(M)=f(W)$ implies $M=W$ or $g(M)=g(W)$ for
all (or most or some) $M$ and $W$. The condition that $π$ is the same for $p$
and $q$, translates to

$\displaystyle\text{Constraint on $M$ and
$W$:}~{}~{}~{}M^{a}_{s+}=π(a|s)=W^{a}_{s+}~{}~{}~{}\text{and in
particular}~{}~{}~{}M^{+}_{s+}~{}=~{}1~{}=~{}W^{+}_{s+}$ (7)

We name the two most interesting equation versions as follows:

$\displaystyle\text{EqIM}(ia):~{}~{}$ $\displaystyle
B^{aa^{\prime}...a^{i}}~{}$ $\displaystyle:=$ $\displaystyle
M^{a}M^{a^{\prime}}...M^{a^{i}}⊘(M^{+})^{i}~{}$
$\displaystyle\stackrel{{\scriptstyle?}}{{=}}~{}W^{a}W^{a^{\prime}}...W^{a^{i}}⊘(W^{+})^{i}$
(8) $\displaystyle\text{EqIM}(i+):~{}~{}$ $\displaystyle B^{a+...+}~{}$
$\displaystyle:=$ $\displaystyle M^{a}(M^{+})^{i-1}⊘(M^{+})^{i}~{}$
$\displaystyle\stackrel{{\scriptstyle?}}{{=}}~{}W^{a}(W^{+})^{i-1}⊘(W^{+})^{i}$
(9)

We allow $M^{+}_{ss^{\prime}}=0$ and keep probabilistic convention that
$p(a|ss^{\prime})=π(a|s)p(s^{\prime}|sa)/p(s^{\prime}|s)$ is undefined iff
$p(s^{\prime}|s)=0$ (see end of Appendix I and Appendix O for more
discussion). Formally, $B^{a}_{ss^{\prime}}=\bot=0/0$ iff
$M^{+}_{ss^{\prime}}=0$, also $W^{+}_{ss^{\prime}}=0$ iff
$M^{+}_{ss^{\prime}}=0$, and similarly for larger $i$.

##  3 Degenerative Cases

To get some feeling about why these questions are so more intricate than
analogous ones in forward models, we consider some simple examples and special
cases first Some further special cases (deterministic planning, deterministic
reachability, and deterministic inverse models) are considered in Appendix D.
There is a strong relationship between the examples violating (i,iii,v) and
counter-examples to seemingly different conjectures found in related work. See
Section C for details.

### Example violating (i,iii,v).

A specific example for $M$ and $W$ which satisfy EqIM$(1)$ but violate
EqIM$(2+)$ and hence EqIM(2$a$) is as follows:

$\displaystyle\textstyle M^{0}={\textstyle\frac{1}{4}}{0~{}2\choose
1~{}1},~{}~{}~{}M^{1}={\textstyle\frac{1}{4}}{2~{}0\choose
1~{}1},~{}~{}~{}W^{0}={\textstyle\frac{1}{2}}{0~{}1\choose
1~{}0},~{}~{}~{}W^{1}={\textstyle\frac{1}{2}}{1~{}0\choose 1~{}0}$

which satisfies (7) ($M^{a}_{s+}={\textstyle\frac{1}{2}}=W^{a}_{s+}$). In this
example, $M^{+}={\textstyle\frac{1}{2}}{1~{}1\choose 1~{}1}$ and
$W^{+}={\textstyle\frac{1}{2}}{1~{}1\choose 2~{}0}$, which shows
$M^{a}⊘M^{+}=W^{a}⊘W^{+}$, except that $W^{+}_{22}=0≠1=M^{+}_{22}$, hence
there is one “dubious” $1\smash{\stackrel{{\scriptstyle?}}{{=}}}0/0$ case. A
simple calculation shows that EqIM(2+) is violated (w/o any division by 0).
The division by 0 can easily avoided by mixing
$U^{a}_{ss^{\prime}}≡{\textstyle\frac{1}{4}}$ into $M$ and $W$, e.g.
$M↝{\textstyle\frac{1}{2}}(M+U)$ and $W↝{\textstyle\frac{1}{2}}(W+U)$. This
means that the 1-step inverse model $B^{a}$ does not always uniquely determine
the 2-step inverse model $B^{aa^{\prime}}$, i.e. (i,iii,v) can fail.

###  $M=W$.

This trivially implies $g(M)=g(W)$. This means if (i) is true, then trivially
also (iii) and (v), and if (ii) is true, then trivially also (iv) and (vi).

###  $M$ and $W$ are independent $a$.

Note that $M^{a}_{ss^{\prime}}≡p(s^{\prime}|sa)$ independent $a$ implies
$M^{a}_{s+}$ independent $a$, hence $π(a|s)=M^{a}_{s+}=1/k$ independent $a$ as
well, hence $M^{a}={\textstyle\frac{1}{k}}M^{+}$. The latter implies
$B^{aa^{\prime}...a^{i}}=k^{-i}$ is independent of $M$ hence is the same as
for $W$. Since we can choose $M≠W$, this shows that (i) and (ii) and higher
order analogues fail for these degenerate $M$ and $W$.

###  $M$ and $W$ are nearly independent $a$.

The above degeneracy generalizes to $M^{a}_{ss^{\prime}}=M_{ss^{\prime}}π_{a}$
and $W^{a}_{ss^{\prime}}=W_{ss^{\prime}}π_{a}$, i.e. action-independent
dynamics, and state-independent actions, which in turn is a special case of
the tensor product below (with $s=\ddot{s}$ and $\dot{s}≡0$).

###  $M$ and $W$ are independent $s^{\prime}$.

In this case,
$M^{a}_{ss^{\prime}}={\textstyle\frac{1}{d}}M^{a}_{s+}={\textstyle\frac{1}{d}}π(a|s)=W^{a}_{ss^{\prime}}$,
hence is a special case of case $M=W$ above.

###  $M$ and $W$ are independent $s$.

In this case,
$[M^{a}M^{a^{\prime}}]_{ss^{\prime\prime}}=∑_{s^{\prime}}M^{a}_{*s^{\prime}}M^{a^{\prime}}_{*s^{\prime\prime}}=π(a|*)M^{a^{\prime}}_{*s^{\prime\prime}}$.
Also the policy $π(a|s)=M^{a}_{s+}$ is independent $s$. If we assume EqIM(1),
this implies

$\displaystyle[M^{a}M^{a^{\prime}}⊘(M^{+})^{2}]_{ss^{\prime\prime}}~{}=~{}\frac{π(a|*)M^{a^{\prime}}_{*s^{\prime\prime}}}{π(+|*)M^{+}_{*s^{\prime\prime}}}~{}=~{}\frac{π(a|*)W^{a^{\prime}}_{*s^{\prime\prime}}}{π(+|*)W^{+}_{*s^{\prime\prime}}}~{}=~{}[W^{a}W^{a^{\prime}}⊘(W^{+})^{2}]_{ss^{\prime\prime}}$

hence EqIM(2) holds and similarly $\text{EqIM}(i)∀i$. As an example, consider

$\displaystyle\textstyle M^{0}~{}:=~{}{\textstyle\frac{1}{2}}{0~{}1\choose
0~{}1},~{}~{}~{}M^{1}~{}:=~{}{\textstyle\frac{1}{2}}{1~{}0\choose
1~{}0},~{}~{}~{}W^{0}~{}:=~{}{\textstyle\frac{1}{3}}{0~{}1\choose
0~{}1},~{}~{}~{}W^{1}~{}:=~{}{\textstyle\frac{2}{3}}{1~{}0\choose 1~{}0}$

These $M≠W$ satisfy EqIM(1) ($M^{a}⊘M^{+}=2M^{a}=W^{a}⊘W^{+}$), hence
constitute another failure case of (i) and (ii).

### Block-diagonal $M$ and $W$.

For $M={\dot{M}~{}0~{}\choose~{}0~{}\ddot{M}}$ and
$W={\dot{W}~{}0~{}\choose~{}0~{}\ddot{W}}$, all operations ($+-×/⊙⊘$) preserve
the block structure, so the above degenerative cases can be combined, one for
the upper-left block and another for the lower-right block.

### Tensor-product $M$ and $W$.

Let
$[\dot{M}⊗\ddot{M}]_{ss^{\prime}}:=\dot{M}_{\dot{s}\dot{s}^{\prime}}\ddot{M}_{\ddot{s}\ddot{s}^{\prime}}$
with $s:=(\dot{s},\ddot{s})$ and
$s^{\prime}:=(\dot{s}^{\prime},\ddot{s}^{\prime})$ be the tensor product of
$\dot{M}$ and $\ddot{M}$ (not to be confused with the element-wise product
$⊙$). Assume $M^{a}=\dot{M}^{a}⊗\ddot{M}$, where the second factor is action-
independent. In this case,
$M^{a}M^{a^{\prime}}...=(\dot{M}^{a}\dot{M}^{a^{\prime}}...)⊗(\ddot{M}\ddot{M}...)$,
and similarly if $a,a^{\prime},...$ is replaced by $+$, hence
$M^{a}M^{a^{\prime}}...M^{a^{i}}⊘(M^{+})^{i}=\dot{M}^{a}\dot{M}^{a^{\prime}}...\dot{M}^{a^{i}}⊘(\dot{M}^{+})^{i}$
is independent of $\ddot{M}$, and similarly for $W^{a}=\dot{W}^{a}⊗\ddot{W}$.
That means, EqIM($i$) hold if $\dot{M}^{a}=\dot{W}^{a}$, whatever $\ddot{M}$
and $\ddot{W}$ are. This formalizes our motivating example that if some part
of the state ($\ddot{s}$) is not controlled (by $a$) and the dynamics
factorizes
($p(s^{\prime}|sa)=p(\dot{s}^{\prime}|\dot{s}a)p(\ddot{s}^{\prime}|\ddot{s})$)
and the policy is independent $\ddot{s}$ ($π(a|s)=π(a|\dot{s})$), then the
multi-step inverse models (3-5) become much simpler than the forward model
(2), namely independent $\ddot{s}$. This case has been studied in [EMK+22] for
episodic near-deterministic $M$.

##  4 (Non)Uniqueness of Inverse MDP Models

We will now consider EqIM(1) and EqIM(2). We first provide a dimensional
analysis which gives some insight and tentative answers about the solution
space for $W$ (given $B$ or $M$): No, one, finitely many, or a polynomial
variety (of some dimension) of solutions. We then consider EqIM(1) only and
characterize $M$ and $W$ for which it holds. This will be used to provide an
algorithm that can determine a (and in some sense all) solution for $W$ and
hence $B^{aa^{\prime}...}$, given only $B^{a}$. EqIM(1) is quite simple, since
it is effectively linear, but EqIM(2) is quadratic in $W$, which is where the
difficulties start.

### Dimensional analysis / counting solutions.

Assume $k≤d$ and $B^{⋅}$ or $M^{⋅}$ are given. The $kd^{2}$ equations EqIM(1)
in $W$ constitute $(k-1)d^{2}$ (linear) constraints on (the $kd^{2}$ real
entries in) $W$. It’s only $(k-1)d^{2}$, since summing over $a$ gives $d^{2}$
vacuous equations $B^{+}=1=W^{+}⊘W^{+}$. There are $kd$ further (linear)
constraints $W^{a}_{s+}=π(a|s)$. Assuming no further (missed/accidental)
redundancies, this leads to a $kd^{2}-(k-1)d^{2}-kd=d(d-k)$ dimensional
(linear) solution space for $W$. This is consistent with the algorithm below
inferring $B^{aa^{\prime}}$ from $B^{a}$ if all $B^{a}$ have full rank. Hence
the set of solutions for $B^{aa^{\prime}}$ forms a polynomial variety of
dimension at least $d(d-k)$.

If also $B^{a+}$ is given, EqIM($2+$) provides $(k-1)d^{2}$ further
(quadratic) constraints (EqIM($ia$) even provides $(k^{i}-1)d^{2}$
constraints). Since $d(d-k)<(k-1)d^{2}$, this now gives an over-determined
system which generally has no solution. But by assumption, $M$ is a solution,
which gives hope that there may be only one or a finite number of solutions.

We can use the $kd+(k-1)d^{2}$ linear equations to eliminate this number of
variables in $W$, which leaves $(k-1)d^{2}$ quadratic equations, now in only
$d(d-k)$ variables, and no further equality constraints. By Bézout’s bound
[FW89], such a System of Quadratic Equations (SQE), either has a continuum
number of solutions (as in the counter-example of Appendix J) or at most
$2^{d(d-k)}$ solutions (as possibly in the counter-example in Appendix I).
Multiple discrete solutions are often caused by symmetries, so for random
$B^{a}$ and $B^{a+}$ consistent with $M$, the solution may indeed be unique.

### Inferring _some_ $B^{aa^{\prime}}$ from $B^{a}$.

Even if $B^{a}$ does not uniquely determine $B^{aa^{\prime}}$, we can ask for
an _algorithm_ inferring _some_ consistent $B^{aa^{\prime}}$ from $B^{a}$.
Indeed this was our primary goal before realizing that the answer is not
always unique. We know that $B^{a}=W^{a}⊘W^{+}$ for _some_ $W$. This implies
$W^{a}=B^{a}⊙W^{+}$. So $W^{a}=B^{a}⊙J$ for some $J$ independent $a$. We need
to ensure proper normalization $W^{a}_{s+}=π(a|s)$, i.e.
$[B^{a}⊙J]_{s+}=π(a|s)$. This leads to the following algorithm to produce some
(and indeed all) $B^{aa^{\prime}}$:

  * •

Given inverse 1-step model $B^{a}_{ss^{\prime}}:=p(a|ss^{\prime})$ and policy
$π(a|s)$

  * •

For each $s$, choose _some_ $d$-vector $J_{s⋅}$
satisfying the $k$ linear equations
$∑_{s^{\prime}}B^{a}_{ss^{\prime}}J_{ss^{\prime}}=π(a|s)$

  * •

Compute forward model $W^{a}:=B^{a}⊙J$

  * •

Compute 2-step inverse model
$B^{aa^{\prime}}:=W^{a}W^{a^{\prime}}⊘(W^{+})^{2}$

  * •

Then $p(aa^{\prime}|ss^{\prime\prime})≡B^{aa^{\prime}}_{ss^{\prime\prime}}$ is
_some_ solution.

If for every $s$, matrix $B^{⋅}_{s⋅}$ has rank $d$, then $B^{aa^{\prime}}$ is
unique. The equations have no solution _iff_ $B$ is invalid in the sense that
no underlying MDP $M$ could have produced such $B$. This can only happen for
$k>d$, i.e. $B$ based on $M$ have some intrinsic constraints beyond $B^{+}=1$
for $k>d$. For instance $B^{0}={\textstyle\frac{1}{2}}{1~{}0\choose 1~{}0}$,
$B^{1}={\textstyle\frac{1}{2}}{0~{}1\choose 0~{}1}$,
$B^{2}={\textstyle\frac{1}{2}}{1~{}1\choose 1~{}1}$ is inconsistent with
$π(a|s)={\textstyle\frac{1}{3}}$. For unknown $π$, any $J$ with $J_{s+}=1$
will do. In general, the valid $J$ span a linear subspace, but the set of all
consistent $B^{aa^{\prime}}$ forms an algebraic variety of equal or lower
dimension. $B^{aa^{\prime}}$ may even be unique even if $J$ and $W$ are not
(see Section 3). Noting that the ranks of $M^{⋅}_{s⋅}$ and $W^{⋅}_{s⋅}$ are
the same, this gives the precise conditions under which (i) is true:

######  Proposition 1 (Conditions under which (i) is true)

$M^{a}⊘M^{+}=W^{a}⊘W^{+}$ implies $M=W$ iff $M^{⋅}_{s⋅}$ has rank $≥d$ for
every $s$.

For this to be possible at all, we need $k≥d$, i.e. more actions than states.
This is typically not the most interesting regime. See Appendix E for an
alternative derivation of this result without an intermediary algorithm.

We will next show that EqIM(2) removes this limitation, but we do not know of
a general and efficient _algorithm_ for inferring (some)
$B^{aa^{\prime}a^{\prime\prime}}$ from $B^{a}$ and $B^{aa^{\prime}}$. We
cannot even rule out that finding approximate solutions is NP-hard.

### (Non)Uniqueness of Inverse MDP Models for $i≥2$.

Above we have established that $B^{a}$ does not uniquely determine
$B^{aa^{\prime}}$ for the interesting regime of $k<d$. From the dimensional
analysis, providing 2-step inverse model $B^{aa^{\prime}}$ in addition, has
the potential of uniquely determining forward model $W$ and/or multi-step
inverse models $B^{aa^{\prime}a^{\prime\prime}...}$. We have numerically
verified that this is indeed the case for $B^{a}$ and $B^{aa^{\prime}}$ based
on random $M^{a}$. A more detailed analysis of the linear/quadratic structure
of the problem is provided in Appendix F and a rank analyses in Appendices G
and P. Unfortunately, even providing $B^{a}$ _and_ $B^{aa^{\prime}}$ does not
always uniquely determine $M^{a}$, nor higher $B$, and (ii,iv,vi) fail for
some $M^{a}$. Furthermore this remains true for higher $i$-versions, i.e. even
EqIM(1)…EqIM($i$) do not always uniquely determine EqIM($i+1$). We provide
(potential) counter-examples in Appendices H and I, but they involve “bad”
0/0. We discuss what this means at the end of Appendix I. We provide a fully
satisfactory counter-example in Appendix J. If the solution is not unique, the
set of solutions forms a polynomial variety. Its (local) dimension measures
the “number” of other solutions (in a neighborhood). In Appendix P we provide
explicit expressions for the tangent spaces from which these dimension can
efficiently be calculated.

##  5 Linear Relaxation

In Section 4 we provided an algorithm if only $B^{a}$ is given. Here we
consider the $i>1$ case, and derive an algorithm for $k^{i}≥d$, provided the
solution is unique and further conditions on $B$ are met. That is, we require
$i≥\log_{k}(d)$, which is greater than the minimum necessary in theory $i=2$
from the dimensional analysis. E.g. for $i=1$ we recover $k≥d$, and $i=2$
improves this to $k≥\sqrt{d}$, and $i=\lceil\log_{2}(d)\rceil$ works for all
$k$.

### Recursive formulation.

From EqIM(1) we know that $W^{a}=B^{a}⊙W^{+}$. Plugging this into EqIM($ia$)
and abbreviating $a^{:i}:=aa^{\prime}...a^{i}$ and
$a^{<i}:=aa^{\prime}...a^{i-1}$ and $j:=i+1$, this gives

$\displaystyle
B^{a^{:i}}⊙(W^{+})^{i}~{}=~{}(B^{a}⊙W^{+})⋅...⋅(B^{a^{i}}⊙W^{+})$ (10)

If we plug EqIM($(i-1)a$) into EqIM($ia$) and abbreviate $V:=(W^{+})^{i-1}$
this simplifies to

$\displaystyle B^{a^{:i}}⊙(V⋅W^{+})~{}=~{}(B^{a^{<i}}⊙V)⋅(B^{a^{i}}⊙W^{+})$

which written out becomes

$\displaystyle∑_{s^{i}}B^{a^{:i}}_{ss^{j}}V_{ss^{i}}W^{+}_{s^{i}s^{j}}~{}=~{}∑_{s^{i}}B^{a^{<i}}_{ss^{i}}V_{ss^{i}}B^{a^{i}}_{s^{i}s^{j}}W^{+}_{s^{i}s^{j}}$
(11)

### Linear relaxation.

We can consider a linear relaxation of this System of Polynomial Equations
(SPE) by introducing new variables $U_{ss^{i}s^{j}}$ (aiming at
$U_{ss^{i}s^{j}}=V_{ss^{i}}W^{+}_{s^{i}s^{j}}$):

$\displaystyle∑_{s^{i}}A^{a^{:i}}_{ss^{i}s^{j}}U_{ss^{i}s^{j}}~{}=~{}0~{}~{}~{}~{}~{}$
$\displaystyle\text{with}~{}~{}~{}~{}~{}A^{a^{:i}}_{ss^{i}s^{j}}~{}:=~{}B^{a^{:i}}_{ss^{j}}-B^{a^{<i}}_{ss^{i}}B^{a^{i}}_{s^{i}s^{j}}$
(12)

These are $k^{i}d^{2}$ potentially independent linear equations in $d^{3}$
unknowns $U$. The solution can only be unique if $k^{i}≥d$. For random $B$,
for each fixed $(s,s^{j})$, the $k^{i}×d$ matrix $A^{⋅⋅⋅}_{s⋅s^{j}}$ has
indeed full rank $\min\\{k^{i},d\\}≥d$, hence $U_{ss^{i}s^{j}}≡0$ is the only
solution. This is inconsistent with the constraints (7), and hence shows that
(unrestricted random) $B$ do not come from some $M$. This makes the validity
of the $B$’s sometimes semi-decidable in time $O(d^{4}(d+k^{i}))$ or
typically/randomized time $O(d^{5})$. For the $B$’s originating from some $M$,
$\hat{U}_{ss^{i}s^{j}}=(M^{+})^{i-1}_{ss^{i}}M^{+}_{s^{i}s^{j}}$ solves (12).
Since for different $ss^{j}$ the equations in (12) are independent,
$U_{ss^{i}s^{j}}:=\hat{U}_{ss^{i}s^{j}}K_{ss^{j}}$ also solves (12) for any
$K$. In other words, the rank of $A^{⋅⋅⋅}_{s⋅s^{j}}$ is bounded by
$\min\\{k^{i},d-1\\}$, and achieved e.g. for random matrices $B$ _consistent_
with $M$. Since the solution is not unique, for many solutions $U$ there will
be no $W^{+}$ satisfying
$U_{ss^{i}s^{j}}=(W^{+})^{i-1}_{ss^{i}}W^{+}_{s^{i}s^{j}}$, not to speak of
$M^{+}$, even if the original problem (10)+(7) has a unique solution.

### Unique solution by lifted constraints.

So we must (and at least for random $M$ can) make the solution unique by
taking into account the linear constraints (7). Applying them to
$s↝s^{i},s^{\prime}↝s^{j},a↝a^{i}$ and multiplying from the left with
$V_{ss^{i}}$ and using $V_{ss^{i}}=U_{ss^{i}+}$ we lift them to

$\displaystyle∑_{s^{j}}B^{a^{i}}_{s^{i}s^{j}}U_{ss^{i}s^{j}}~{}=~{}U_{ss^{i}+}π(a^{i}|s^{i})~{}~{}~{}~{}~{}\text{and}~{}~{}~{}~{}~{}U_{s++}~{}=~{}1$
(13)

These $kd^{2}+d$ further linear constraints have the potential to make the
solution of (12) unique, i.e. resolve the $d^{2}$ degeneracy $K_{ss^{i}}$. If
so, we can recover
$M^{+}_{s^{i}s^{j}}=W^{+}_{s^{i}s^{j}}=U_{ss^{i}s^{j}}/V_{ss^{i}}$ (and
finally $M^{a}=W^{a}=B^{a}⊙W^{+}$) in polynomial time. It actually suffices to
solve (12) and (13) for one fixed $s$, e.g. $s=1$, which with some care can be
done in time $O(d^{4})$. In practice, for approximate $B$ one would solve a
least-squares problem using all equations or a random projection for speed.

### Algorithm.

Putting pieces together, we have the following algorithm for computing $W^{a}$
and hence $B^{a^{:j}}$ for all $j$ via EqIM($ja$) from $B^{a}$ and
$B^{a^{<i}}$ and $B^{a^{:i}}$

  * •

Given: Policy $π(a|s)$ and for $j-1:=i≥2$, inverse $1,i-1,i$-step models
$B^{a}_{ss^{\prime}}=p(a|ss^{\prime})$ and
$B^{a^{<i}}_{ss^{j}}=p(a^{<i}|ss^{i})$ and
$B^{a^{:i}}_{ss^{j}}=p(a^{:i}|ss^{j})$

  * •

Do the following calculations for one $s$ (e.g. $s=1$),
or a few or all $s$ or some random linear combinations of $s$:

  * •

For each $s^{j}$, let $\hat{U}_{ss^{i}s^{j}}$ be a solution of (12) with
$\hat{U}_{s+s^{j}}=1$

  * •

If a non-zero solution does not exist, set $\hat{U}_{ss^{i}s^{j}}=0~{}∀s^{i}$.

  * •

Optional: If multiple solutions exist, return “$W$ _may_ not be unique”

  * •

If $\hat{U}_{s++}=0$, return “$B$ is not consistent with any $M$”

  * •

Solve $∑_{s^{j}}C^{a^{i}}_{ss^{i}s^{j}}K_{ss^{j}}=0$ and $K_{s+}=1$ for
$K_{s*}$, where
$C^{a^{i}}_{ss^{i}s^{j}}:=(B^{a^{i}}_{s^{i}s^{j}}-π(a^{i}|s^{i}))\hat{U}_{ss^{i}s^{j}}$

  * •

If no solution, return “$B$ is not consistent with any $M$”

  * •

Optional: If multiple solutions exist, return “$W$ _may_ not be unique”

  * •

$\tilde{U}_{ss^{i}s^{j}}:=\hat{U}_{ss^{i}s^{j}}K_{ss^{j}}$,
$U_{ss^{i}s^{j}}:=\tilde{U}_{ss^{i}s^{j}}/\tilde{U}_{s++}$,
$V_{ss^{i}}:=U_{ss^{i}+}$, $W^{+}_{s^{i}s^{j}}:=U_{ss^{i}s^{j}}/V_{ss^{i}}$

  * •

Optional: If different $s$ lead to different $W^{+}$ or $V≠(W^{+})^{i-1}$,
return “$W$ _may_ not be unique”

  * •

Return forward model $W^{a}:=B^{a}⊙W^{+}$ and other inverse $B^{⋅⋅⋅}$ computed
via (8)

### Variations that don’t work.

For unknown $π$, we only have $d$ lifted constraints $U_{s++}=1$, which are
not sufficient to make the solution unique, also resulting in too many
solutions for the relinearization trick [CKPS00] to work. The same is true if
we had relaxed $U_{ss^{\prime}s^{j}}=W^{+}_{ss^{\prime}}V_{s^{\prime}s^{j}}$.
If we had applied linear relaxation directly to EqIM($ia$), this would have
led to order-$i+1$ tensors and require $k≥d^{1-1/i}$, which is much worse than
$k≥d^{1/i}$ for $i>2$. Including $B^{a^{:j}}$ and EqIM($ja$) for some or all
$j<i-1$ is not only unhelpful but even counter-productive.

##  6 Experiments

The algorithm described in Section 5 was motivated by the dimensional analysis
and properties of random matrices. Namely, that $A^{⋅⋅⋅}_{s⋅s^{j}}$ is likel
“full” rank, and thus yielding a unique solution. In order to explore the
plausibility of this assumption in practice, we have evaluated the algorithm
on a set of toy (but structured) environments. This includes the canonical
‘four-rooms’ grid-world and samples from the distribution over all grid-worlds
of that size. All environments have $k=5$ (local movement on the grid) and
$d=24$, thus satisfying the $k\geq d^{2}$ constraint which permits solving
EqIM($2$).

### Experiments on naturalistic environments.

As detailed in Appendix N, for all environments tested the algorithm yielded a
unique solution (recovering $M^{a}$) up to a reasonable level of numerical
precision. This remained true even after injecting noise (across several
orders of magnitude) into the environmental transition dynamics. This is in
contrast to related methods which rely on near-deterministic environments
[EMK+22].

This result is non-trivial, as the statistics of these environments differ
significantly from those produced by random matrices. For example, grid-world
dynamics are both local and sparse, unlike random matrix dynamics which almost
always have non-zero probability for all transitions. It remains to be seen
whether or not larger-scale environments yield similar results, but it is at
least non-obvious what additional environmental properties would break the
constraints of the algorithm.

Figure 1: Environments, their transition matrices (i.e. $M^{+}$) and the
matrices inferred by the algorithm (i.e. $W^{+}$). Results shown on the most
and least noisiest variants of each environment. _Top_ ‘four-rooms’ grid-
world. _Bottom_ One of the randomly generated grid-worlds.

### Experiments illustrating robustness to noise.

The propositions (and previous experimental result) assume that we know the
one and two step inverse models ($B1:=B^{a}$, $B2:=B^{a+}$) exactly, but in
practice these distributions must be estimated from data. Here we investigate
the extent to which our algorithm is robust to noise arising from learning.

Rather than committing to a specific learning algorithm, we instead directly
inject noise into the true inverse distributions. Figure 3 shows that noise
doesn’t substantially degrade performance across several orders of magnitude
(see Appendix N for details). Additionally, the effect of this noise is
substantially diminished as the horizon of the inverse model is increased
(from $B1:=B^{a}$ to $B3:=B^{a++}$). While the is perhaps not surprising, as
the entropy of such inverse distributions increases monotonically with the
horizon, it still shows that noise is not compounding in a way that renders
long-horizon predictions meaningless.

### Experiments on the Tensor-product special case.

As detailed in Section 3, if $M$ factors into two processes
$\dot{M}^{a}⊗\ddot{M}$, where $\ddot{M}$ is action-independent, then only the
complexity of the action-dependent process $\dot{M}^{a}$ matters for all of
our questions. The significance of this special case, as well as the details
of environments construction, can be found in Appendix N.

The linear algorithm of Section 4 can (implicitly) output all $W$ and $B2$
consistent with $B1$, and the formulas derived in Appendix P allow to
(explicitly) calculate the dimensions of the solution spaces.

In the experiments shown in Figure 3, the environments complexity is
systematically varied. The results show that the space of forward dynamics $W$
is always larger than the space of the 2-step inverse models ($B2$). This
confirms that inverse models can be simpler than forward models.

Figure 2: Noise-induced reconstruction error: In practice $W$ must be inferred
from learned estimates of $B1$ and $B2$. We investigate the effect of the
resulting error on the inverse models ($B1,B2,B3$) recovered from the inferred
$W$ in terms of their proximity to the ground truth distributions. At each
noise level the algorithm was run on 10 randomly generated grids, with the
shaded region representing $±2σ$.

Figure 3: Solution dimensions of $W$ and $B2$ given $B1$: When the solution to
an inverse model ($B2$) given only $B1$ is not unique, we can characterize the
solution space in terms of its manifold dimension. By comparing this to the
dimension of that of the inferred forward model ($W$), we can see that our
algorithm has narrowed down the space of inverse models further. If also $B2$
is given, the solution dimension of $W$ reduces from $d_{W}$ (blue curve) to
$d_{W}-d_{B}$ (blue minus orange curve).

##  7 Computational Complexity

Maybe even just characterizing all $M$ for which EqIM(1) and EqIM(2) uniquely
determine $W$ is hopeless, not to speak of finding some or all $W$ in case
not. More formally, we can ask the question of whether there exists an
efficient algorithm that can decide whether $\text{EqIM}(i)$ has a unique
solution. We provide some weak preliminary evidence, why this problem may be
NP-hard. Appendix M contains fully self-contained a few versions of this open
problem in their simplest instantiation and most elegant form.

### Decidability and computability.

EqIM(2) converted to (27) and (7), or (14) or (15) below form a System of
Quadratic Equations (SQE). The constraint $W≠M$ can also be expressed as a
quadratic equation (see below). As such, the existence and uniqueness of
solutions is formally decidable by computing a Gröbner basis [Stu02], and
(some) solutions can be found by cylindrical algebraic decomposition in
(double) exponential time. $ε$-approximate solutions can of course be found by
exponential brute-force search through all $W$ on a finite $ε^{\prime}$-grid,
and verified in polynomial time.

### Complexity considerations.

3SAT is NP complete. A CNF formula in $n$ boolean variables can easily be
converted to a multivariate (cubic) polynomial equation $p(x_{1},...,x_{n})=0$
over binary variables $x_{i}$. The binary variables can be relaxed to real-
valued variables by considering
$p(x_{1}...x_{n})^{2}+∑_{i}x_{i}^{2}(1-x_{i}^{2})=0$. By introducing dummy
variables, e.g. $z_{i}=x_{i}^{2}$ and others, this can further be reduced to a
System of Quadratic Equations (SQE). Therefore SQE is also NP hard. EqIM(2+)
explicitly written in quadratic form

$\displaystyle M^{a}M^{+}⊙(W^{+})^{2}-W^{a}W^{+}⊙(M^{+})^{2}~{}=~{}0$ (14)

constitutes an SQE in $W$ given $M$, also if we include linear EqIM(1) and
probability constraints (7). Non-negativity of $W$ can be enforced with
(slack) variables $(Y^{a}_{ss^{\prime}})^{2}=W^{a}_{ss^{\prime}}$. (Similarly
(20) plus constraints (19) constitute an SQE in $Δ$.) To reduce the uniqueness
question to a solvability problem we need to avoid the trivial solution $W≡M$,
e.g. by introducing further (slack) variables $t∈ℝ$ and
$Γ^{a}_{ss^{\prime}}:=(W^{a}_{ss^{\prime}}-M^{a}_{ss^{\prime}})^{2}$ and
constraint $t⋅Γ^{+}_{++}=1$. Due to the minus sign in (14), this can _not_ be
converted to a convex (optimization) problem. The choice of $M$ gives
significant freedom in creating SQE problems, even if only considering
permutation matrices $M^{a}∈\\{0,1\\}^{d×d}$. If one could show that every SQE
can be represented as (14) [plus $W≠M$ constraint] for a suitable choice of
$M$, this would imply that proving the existence of $W≠M$ satisfying (14) is
NP hard. This in turn would imply that computing (any)
$p(a|ss^{\prime\prime\prime})$ from $p(a|ss^{\prime})$ and
$p(a|ss^{\prime\prime})$ is NP hard.On the other hand, matrix multiplication
$W^{a}W^{b}$ is a very specific quadratic form, which may not be flexible
enough to incorporate every SQE within (14).

We could not find any work on NP-hardness of Systems of Polynomial Matrix
Equations (SPME). There is work on the NP-hardness of tensor problems [HL13],
but this refers to the design tensors, e.g.
$∑_{jk}A^{jk}_{i}x_{j}x_{k}+∑_{j}B^{j}_{i}x_{j}+C_{i}=0~{}∀i$, but the
unknowns are always treated as scalars or vectors. Of course
$[X⋅Y]_{ik}=∑_{abcd}A^{abcd}_{ik}X_{ab}Y_{cd}$, but
$A^{abcd}_{ik}=δ_{ai}δ_{dk}δ_{bc}$ is a very special fixed tensor (actually of
low tensor rank $d$) with no flexibility of encoding NP-hard problems therein.

That inference in Bayesian networks is NP-complete [KF09] does not help us
either for two reasons: First, in our problem the probability distribution
over states and actions is only partially given. More importantly, our network
for $i=2$ has only 5 nodes ($s,a,s^{\prime},a^{\prime},s^{\prime\prime}$),
while the NP-hardness proofs we are aware of require large networks. Even for
fixed $i>2$, it is not obvious how to encode NP-hard problems into EqIM(i),
due to the severe structural constraints in EqIM(i) compared to a general
network with $2i+3$ nodes. It is not clear how to exploit the fact that our
(few) state nodes are large.

SQE are polynomially equivalent to Systems of Quadratic Matrix Equations
(SQME), which may be the reason complexity theorists have ignored the latter.
We suspect but do not know whether SQME of _bounded_ structural complexity
(only the definitions of the constant matrices scale with $d×d$) is NP-hard
(Open Problem 7). If we allow sparse encoding of SQE variables in $W$, i.e. we
allow one equation involving $⊙$ of the form $B⊙W=0$ with boolean matrix $B$,
then bounded SQME becomes NP-hard. See Appendix K for details. Below we
directly reduce 1in3SAT to a Bounded-SQME with $⊙$ that resembles our problem
as close as we were able to make it.

### An NP-complete matrix problem.

From EqIM(1) we know that $W^{a}=B^{a}⊙W^{+}$. Plugging this into EqIM($2a$)
gives

$\displaystyle
B^{aa^{\prime}}⊙(W^{+}⋅W^{+})~{}=~{}(B^{a}⊙W^{+})(B^{a^{\prime}}⊙W^{+})~{}~{}~{}\text{with
constraints}~{}~{}~{}[B^{a}⊙W^{+}]_{s+}=π(a|s)$ (15)

This set of equations is purely in terms of what is given ($B^{a}$ and
$B^{aa^{\prime}}$) and only involves unknowns $W^{+}$ without reference to
$W^{a}$. See Appendix L for some further simplification and discussion. We
will show:

######  Proposition 2 (An NP-complete matrix problem)

Given $A,B,C,Π$, deciding whether the following quadratic matrix problem has a
solution in $W$ is NP-complete:

$\displaystyle
A⊙(W⋅W)=(C⊙W)(C⊙W),~{}~{}~{}~{}~{}[B⊙W]_{s+}=1,~{}~{}~{}~{}~{}Π⋅W=W$ (16)

This has some resemblance to (15). Since the boundary between P and NP is very
fractal/subtle, this in-itself may not imply much, but is more meant as a
demonstration of how one may approach proving NP-hardness of (15).

Proof. We reduce 1in3SAT, which is an NP-complete variant of 3SAT, where each
clause must have exactly one satisfying assignment, to (16). A 3CNF($n,m,g$)
formula is a boolean conjunction of $m$ clauses in $n$ variables, where each
clause $c_{i}=\ell_{i1}\dot{∨}\ell_{i2}\dot{∨}\ell_{i3}$ for $i∈\\{1:m\\}$ is
a 1-in-3 disjunction of 3 literals, and each literal is $\ell_{ia}=x_{j}$ or
it’s complement $\ell_{ia}=¬x_{j}≡\bar{x}_{j}$, where $j=g(i,a)$ is the
variable index of clause $i$ in position $a$.

We arithmetize the 3CNF expression in the standard way by replacing True$↝1$,
False$↝0$, and $\dot{∨}↝+$, i.e. we ask whether the system of linear equations
$\ell_{i1}+\ell_{i2}+\ell_{i3}=1~{}∀i$ has a solution in $x_{j}∈\\{0,1\\}$. We
need to encode the $x$’s into $W$ somehow: We aim at the following embedding:

$\displaystyle
W~{}=~{}\left(\begin{array}[]{cccccccc}x_{1}&\bar{x}_{1}&...&x_{n}&\bar{x}_{n}&y_{0}&...&y_{k}\\\
⋮&⋮&\ddots&⋮&⋮&⋮&\ddots&⋮\\\
x_{1}&\bar{x}_{1}&...&x_{n}&\bar{x}_{n}&y_{0}&...&y_{k}\end{array}\right)$

The $y$ are $k+1:=\max\\{1,m-n+2\\}$ extra dummy variables to make the matrix
a square $d×d$ matrix with $d:=\max\\{m+n+2,2n+1\\}$.

Choosing a cyclic permutation matrix $Π=[234...d1]$ ensures that all rows of
$W$ are indeed the same via $Π⋅W=W$. The standard way of achieving
$x_{j},y_{j}∈\\{0,1\\}$ is via $x_{j}^{2}=x_{j}$ and $y_{j}^{2}=y_{j}$. This
can be achieved via $(\text{Id}⊙W)^{2}=\text{Id}⊙W$, were
$\text{Id}_{ss^{\prime}}=δ_{ss^{\prime}}$ is the identity matrix.

We use $[B⊙W]_{s+}=1$ to ensure $\bar{x}_{j}=1-x_{j}$, $y_{0}=1$, and
$y_{1}=...=y_{k}=0$ and $\ell_{i1}+\ell_{i2}+\ell_{i3}=1$ by setting
$B_{s,2s-1}=B_{s,2s}=1$ for $s∈\\{1:n\\}$, and $B_{i+n,2j-1}=1$ if
$\ell_{ia}=x_{j}$ and $B_{i+n,2j}=1$ if $\ell_{ia}=¬x_{j}$ for $i∈\\{1:m\\}$
and $a∈\\{1,2,3\\}$, and $B_{d-1,2n+1}=...=B_{d-1,2n+m}=1$, and
$B_{d,2n+1}=1$, and $B_{ss^{\prime}}=0$ for all other $ss^{\prime}$. This also
ensures that all rows of $W$ sum to $n+1$, hence $W⋅W=(n+1)W$, so
$x_{j}∈\\{0,1\\}$ can be achieved via $C=\text{Id}$ and
$A={\textstyle\frac{1}{n+1}}\text{Id}$ in $A⊙(W⋅W)=(C⊙W)(C⊙W)$.

The construction implies that the 3CNF($n,m,g$) formula is satisfiable _iff_
(16) has a solution in $W$ with the $A,B,C,Π$ as constructed above. This shows
NP-hardness of deciding whether (16) has a solution. A solution can trivially
be verified (in the rationals or to $ε$-precision over the reals) in time
$O(d^{3})$, hence the problem is in NP, hence NP-complete.

######  Conjecture 3 (NP-hardness)

Deciding (ii), (iv), (vi) is NP-hard. Deciding whether $B^{a}$ and
$B^{aa^{\prime}}$ are consistent with some $M$ is also NP-hard. Computing some
solution is FNP-hard.

##  8 Conclusion

### Summary.

We have shown that the 1-step inverse model $p(a|ss^{\prime})$ does not
uniquely determine the 2-step probabilities $p(a|ss^{\prime\prime})$ if there
are less actions than states ($k<d$). Even for $k≥d$, the implication can
fail, e.g. if the extra actions are ineffective, but if
$p(s^{\prime}|sa)=M^{a}_{ss^{\prime}}$ considered as matrices in $a$ and
$s^{\prime}$ for each $s$ have full rank, the implication holds. Even
providing $p(aa^{\prime}...a^{j-1}|ss^{j})$ for all $j<i$ not necessarily
determines $p(a|ss^{i})$. Since the involved SPE is (heavily) over-determined,
we expect the failure cases to be sparse/rare in some sense. For ($B$ based
on) random $M$, we provided evidence that $a=2$ suffices to determine $M$ and
hence $p(aa^{\prime}...|ss^{\prime}s^{\prime\prime}...)$ from
$p(a|ss^{\prime})$ and $p(a|ss^{\prime\prime})$. For low-rank $M$ the
implication may fail. Finally we investigated whether determining (uniqueness
of) $M$ given $B^{a}$ and $B^{aa^{\prime}}$ could be an NP-hard problem, but
the question remains open.

### Open Problems.

Even if our problem is NP-hard, we can still come up with uniqueness
characterizations, but their evaluation would be NP-hard, so their usefulness
probably limited. On the other hand, we have derived efficient algorithms
(essentially solving some linear equations) which cannot only efficiently
determine whether the solution is locally unique (Appendix F) but also the
solution dimension in a neighborhood of a solution (Appendix P). Note that it
is always possible to find all $ε$-approximate solutions of (14) by an
exponential brute-force search through all $W$ on a finite $ε^{\prime}$-grid,
or verify an $ε$-approximate solution in polynomial time.

The counter-examples in Appendices H,I,J involve cyclic permutations and
larger $i$ required larger cycles and hence larger $d$. It would be
interesting to know whether solutions become unique for a-periodic MDPs and/or
fixed $d$ for $i\textrightarrow∞$. The low-rank analysis weakly indicates the
existence of counter-examples for mixing MDPs. If counter-examples have to
scale exponentially with $i$ as they currently do, then they may not be
relevant in practice. The state space in practice is infinite, but the
relevant part may behave like finite.

### Discussion.

Given our analysis, we would expect that in practice, $B^{a}$ and
$B^{aa^{\prime}}$ determines $B^{aa^{\prime}a^{\prime\prime}...}$ and $W$
sufficiently well. Sufficiently well in case of $W$ means all and only those
aspects of the forward model relevant for the inverse model. Then of course
the question remains how to compute the/an answer. While the linear relaxation
developed in Section 5 fails for $k<d^{1/i}$ as an exact method, it might
still lead to useful approximate solutions [Stu02] without formal guarantees.
Indeed, EqIM($ia$) is heavily over-determined for $i≥2$, and heuristic solvers
often work well in this regime.

Handling non-uniqueness: In practice, the state space is very often infinite,
and no finite amount of data will determine even $B^{a}$ uniquely without
further structural assumptions. Neural networks intrinsically restrict the
solution space, but this may not suffice for modern over-parametrized deep
networks. Aiming for the maximum-entropy distribution consistent with the
(constraints from) data is popular, and could make the solution unique, as
well as any other optimization constraint.

## References

  * [AD21] Saurabh Arora and Prashant Doshi.  A survey of inverse reinforcement learning: Challenges, methods and progress.  Artificial Intelligence, 297:103500, August 2021.

  * [BEP+18] Yuri Burda, Harri Edwards, Deepak Pathak, Amos Storkey, Trevor Darrell, and Alexei A Efros.  Large-scale study of curiosity-driven learning.  arXiv preprint arXiv:1808.04355, 2018.
  * [BPK+20] Adrià Puigdomènech Badia, Bilal Piot, Steven Kapturowski, Pablo Sprechmann, Alex Vitvitskyi, Zhaohan Daniel Guo, and Charles Blundell.  Agent57: Outperforming the atari human benchmark.  In International Conference on Machine Learning, pages 507–517. PMLR, 2020.
  * [CKPS00] Nicolas Courtois, Alexander Klimov, Jacques Patarin, and Adi Shamir.  Efficient Algorithms for Solving Overdefined Systems of Multivariate Polynomial Equations.  In Gerhard Goos, Juris Hartmanis, Jan van Leeuwen, and Bart Preneel, editors, Advances in Cryptology — EUROCRYPT 2000, volume 1807, pages 392–407. Springer Berlin Heidelberg, Berlin, Heidelberg, 2000.
  * [DY79] E. B. Dynkin and A. A. Yushkevich.  Controlled Markov processes.  Number 235 in Grundlehren der mathematischen Wissenschaften. Springer-Verlag, Berlin ; New York, 1979.
  * [EGIL18] Benjamin Eysenbach, Abhishek Gupta, Julian Ibarz, and Sergey Levine.  Diversity is all you need: Learning skills without a reward function.  arXiv preprint arXiv:1802.06070, 2018.
  * [EMK+22] Yonathan Efroni, Dipendra Misra, Akshay Krishnamurthy, Alekh Agarwal, and John Langford.  Provable RL with Exogenous Distractors via Multistep Inverse Dynamics.  arXiv:2110.08847 [cs], March 2022.
  * [FW89] William Fulton and Richard Weiss.  Algebraic Curves: An Introduction to Algebraic Geometry.  Addison-Wesley, 1989.
  * [GRW17] Karol Gregor, Danilo Jimenez Rezende, and Daan Wierstra.  Variational Intrinsic Control.  In Workshop, February 2017.
  * [HDB+21] Steven Stenberg Hansen, Guillaume Desjardins, Kate Baumli, David Warde-Farley, Nicolas Heess, Simon Osindero, and Volodymyr Mnih.  Entropic Desired Dynamics for Intrinsic Control.  In Advances in Neural Information Processing Systems, May 2021\.
  * [HL13] Christopher J. Hillar and Lek-Heng Lim.  Most Tensor Problems Are NP-Hard.  Journal of the ACM, 60(6):1–39, November 2013.
  * [HSHB99] Jesse Hoey, Robert St-Aubin, Alan Hu, and Craig Boutilier.  SPUDD: Stochastic planning using decision diagrams.  In Proceedings of the Fifteenth Conference on Uncertainty in Artificial Intelligence, pages 279–288, 1999.
  * [KF09] Daphne Koller and Nir Friedman.  Probabilistic Graphical Models: Principles and Techniques.  Adaptive Computation and Machine Learning. MIT Press, Cambridge, MA, 2009.
  * [LFLDP21] Bonnie Li, Vincent François-Lavet, Thang Doan, and Joelle Pineau.  Domain adversarial reinforcement learning.  arXiv preprint arXiv:2102.07097, 2021.
  * [MHKL20] Dipendra Misra, Mikael Henaff, Akshay Krishnamurthy, and John Langford.  Kinematic state abstraction and provably efficient rich-observation reinforcement learning.  In International conference on machine learning, pages 6961–6971. PMLR, 2020.
  * [MJR15] Shakir Mohamed and Danilo Jimenez Rezende.  Variational information maximisation for intrinsically motivated reinforcement learning.  Advances in neural information processing systems, 28, 2015.
  * [OKD+21] Pedro A. Ortega, Markus Kunesch, Grégoire Delétang, Tim Genewein, Jordi Grau-Moya, Joel Veness, Jonas Buchli, Jonas Degrave, Bilal Piot, Julien Perolat, Tom Everitt, Corentin Tallec, Emilio Parisotto, Tom Erez, Yutian Chen, Scott Reed, Marcus Hutter, Nando de Freitas, and Shane Legg.  Shaking the foundations: Delusions in sequence models for interaction and control.  arXiv:2110.10819 [cs], October 2021.
  * [PGJ16] Judea Pearl, Madelyn Glymour, and Nicholas P. Jewell.  Causal Inference in Statistics: A Primer.  Wiley, Chichester, West Sussex, 2016.
  * [Pre00] Doina Precup.  Temporal Abstraction in Reinforcement Learning, 2000.
  * [SGL+19] Archit Sharma, Shixiang Gu, Sergey Levine, Vikash Kumar, and Karol Hausman.  Dynamics-Aware Unsupervised Discovery of Skills.  In International Conference on Learning Representations, September 2019.
  * [SP02] Martin Stolle and Doina Precup.  Learning Options in Reinforcement Learning.  In Sven Koenig and Robert C. Holte, editors, Abstraction, Reformulation, and Approximation, Lecture Notes in Computer Science, pages 212–223, Berlin, Heidelberg, 2002. Springer.
  * [Stu02] Bernd Sturmfels.  Solving Systems of Polynomial Equations.  Number 97 in Regional Conference Series in Mathematics. American Mathematical Society, Providence, RI, 2002.
  * [WDG+16] Grady Williams, Paul Drews, Brian Goldfain, James M Rehg, and Evangelos A Theodorou.  Aggressive driving with model predictive path integral control.  In 2016 IEEE International Conference on Robotics and Automation (ICRA), pages 1433–1440. IEEE, 2016.

##  Appendix A List of Notation

Symbol  |  Type |  Explanation
---|---|---
$\bot$ |  | undefined
$[\mskip-3.0mu[\text{bool}]\mskip-3.0mu]$ |  $∈\\{0,1\\}$ | =1 if bool=True, =0 if bool=False
$δ_{ss^{\prime}}$ |  $:=[\mskip-3.0mu[s=s^{\prime}]\mskip-3.0mu]$ | Kronecker delta
$d$ |  $∈ℕ$ | number of states
$k$ |  $∈ℕ$ | number of actions
$i,j$ |  $∈ℕ$ | time index/step
$\\{i:j\\}$ |  $⊂ℤ$ | Set of integers from $i$ to $j$ (empty if $j<i$)
$s,s^{\prime},...,s^{i}$ |  $∈\\{1:d\\}$ | state at time step 1,2,…,$i$
$a,a^{\prime},...,a^{i}$ |  $∈\\{0:k-1\\}$ | action at time step 1,2,…,$i$
$b,b^{\prime},...,b^{i}$ |  $∈\\{0:k-1\\}$ | alternative action at time step 1,2,…,$i$
$a^{:i}$ |  $:=~{}aa^{\prime}...a^{i}$ | sequence of $i$ actions
$a^{<i}$ |  $:=~{}aa^{\prime}...a^{i-1}$ |  sequence of $i-1$ actions
$\dot{s},\ddot{s}$ |  $∈\\{1:\dot{d}\\}$ | parts of state, usually $s=(\dot{s},\ddot{s})$
$ε$ |  $>0$ | small number $>0$
$p(...)$ |  $∈[0;1]$ | (conditional) probability distribution over states and actions
$π(a|s)$ |  $∈[0;1]$ | policy. Probability of action $a$ in state $s$
$M^{a},W^{a}$ |  $∈[0;1]^{d×d}$ | transition-policy tensor $M^{a}_{ss^{\prime}}=p(s^{\prime}|sa)⋅π(a|s)$, similarly $W=q$
$B^{a}$ |  $∈[0;1]^{d×d}$ | inverse 1-step model $B^{a}_{ss^{\prime}}=p(a|ss^{\prime})$ for each action $a$
$B^{a++}_{ss^{\prime\prime\prime}}$ |  $∈[0;1]$ | 3-step first-action inverse model $p(a|ss^{\prime\prime\prime})$
$J,K,Δ$ |  $∈ℝ^{d×d}$ | action-independent $d×d$ “transition” matrices
${}^{+}~{}_{+}$ |  $⋅^{n}\textrightarrow⋅$ | index summation, e.g. $M^{+}_{s+}=∑_{as^{\prime}}M^{a}_{ss^{\prime}}$
$~{}⋅$ |  $(⋅,⋅)\textrightarrow⋅$ | matrix multiplication: $[AB]_{ss^{\prime\prime}}=∑_{s^{\prime}}A_{ss^{\prime}}B_{s^{\prime}s^{\prime\prime}}$
$⊙$ |  $(⋅,⋅)\textrightarrow⋅$ | element-wise multiplication of matrix elements: $[A⊙B]_{ss^{\prime}}=A_{ss^{\prime}}B_{ss^{\prime}}$
$⊘$ |  $(⋅,⋅)\textrightarrow⋅$ | element-wise division of matrix elements: $[A⊘B]_{ss^{\prime}}=A_{ss^{\prime}}/B_{ss^{\prime}}$
$⊗$ |  $(⋅,⋅)\textrightarrow⋅$ | tensor product: $[\dot{M}⊗\ddot{M}]_{ss^{\prime}}:=\dot{M}_{\dot{s}\dot{s}^{\prime}}\ddot{M}_{\ddot{s}\ddot{s}^{\prime}}$ with $s=(\dot{s},\ddot{s})$ and $s^{\prime}=(\dot{s}^{\prime},\ddot{s}^{\prime})$

##  Appendix B Application to Planning

In Section 1, various streams of applied work were highlighted; here we focus
on spelling out the overarching impact that compositional inverse models (an
affirmative answer to question (iv)) would have for planning problems.

Many forms of planning involve the evaluation of candidate $i$-step action
sequences (e.g. model predictive path integral control [WDG+16]). Ideally, all
possible action sequences would be evaluated, but as the space of $i$-step
action sequences grows exponentially in $i$, this is often intractable.

Access to the $i$-step inverse distribution $p(a...a^{i}|s...s^{i+1})$ allows
determining the subset of action sequences that likely reach state $s^{i+1}$
post-execution (e.g. those whose probability is above some threshold). It is
often the case that only action sequences that are distinguished in this way
are of interest (e.g. goal-reach tasks), thus access to an inverse model of
the appropriate horizon allows for filtering candidates. This filtering method
is a particularly appealing approach when the cost/reward function is
initially unknown and frequently changes, as in [MJR15].

### Motivating Example.

Consider an agent who has control over $\dot{s}$ but not over $\ddot{s}$. For
instance a robot equipped with a camera can control its position and
orientation, but not the shape and color of objects in its path. The forward
model $p(s^{\prime}|as)$ essentially involves modelling the whole observable
world. The inverse model $p(a|ss^{\prime})$ on the other hand can ignore
inputs that the agent has no control over. Of course in practice, $s$ does not
come neatly separated into $\dot{s}$ and $\ddot{s}$, so a (say) deep neural
network still has to learn the controllable features, but neither needs to
learn nor predict the uncontrollable features (under the factorization
assumptions described in Section 3, now in feature space).

If the goal is to navigate from $s$ to $s^{i}$ in $i$ time steps, and open-
loop control suffices as e.g. in (near)-deterministic problems [EMK+22], then
action sequences for which $p(aa^{\prime}...a^{i-1}|ss^{i})$ is large are the
most likely that caused the transition to $s^{i}$, hence these sequences are
promising candidates for macro actions (temporally extended actions, options)
in Reinforcement Learning [SP02, Pre00].

Since the action space is typically much smaller than the state space (the
former often finite, the latter often even infinite-dimensional), even
learning $p(aa^{\prime}...a^{i-1}|s...s^{i})$ directly for all small $i$ can
be feasible and may be more efficient than learning the one-step forward
model. A closed-loop alternative would be to learn only $p(a|s...s^{i})$, find
the likely first action $a$ that caused the ultimate transition to $s^{i}$,
then take action $a$, iterate, and store the resulting sequence as an option.

The required sample complexity to learn inverse MDP models for larger $i$
directly from data may grow exponentially in $i$, which is why inferring
$i$-step inverse models from 1-step and 2-step inverse models would be useful.
The fact that this problem borders NP-hardness probably prevents even powerful
transformer models to finding the structure in
$p(aa^{\prime}...a^{i-1}|s...s^{i})$ by themselves.

##  Appendix C Counter-Examples in Related Work

In Section 3 we presented a counter-example to questions (i,iii,v). Question
(i) (i.e. Can $M$ be inferred from $B^{a}:=M^{a}⊘M^{+}$?) has been implicitly
addressed in previous work. In [EMK+22, App.A.3] the authors present a
counter-example to the claim that a state representation constructed via an
inverse model (i.e. two states have the same representation iff they yield the
same inverse distribution for all of their possible successor states) is
sufficient for representing a set of policies that differentially visit all
states. This fails whenever two states are aliased by the inverse model.
Technically, as per their Definition 2, this ‘policy cover’ need only account
for all ‘endogenous’ states. But omit the ‘exogenous’ states from their
counter-example and it can be seen to address our question (i).

Note that this failure of state representation learning implies a negative
answer to our question (i), as $W$ would differ from $M$ on these aliased
states. Unlike our counter-example, theirs involves deterministic forward
dynamics, and therefor buttresses our claims by showing that $M$ cannot always
be inferred even in this simpler case. Similar to our counter-example in
Section 3, [MHKL20] proposes a stochastic counter-example to inverse modeling
for state representation learning.

In general, the transferability of these counter-examples suggests a strong
relationship between the literature on using single-step inverse models for
state representation learning and using them for inferring the forward model.
It is an interesting open question whether or not algorithms for
representation learning on the basis of multi-step inverse models (like those
put forward in [EMK+22]) might be used to shed light on the questions put
forward here and vice versa.

##  Appendix D Deterministic Cases

### Deterministic planning / reachability problem.

If we are only interested in finding _some_ action sequence
$aa^{\prime}...a^{i}$ that leads to $s^{i}$, the problem becomes easy: The
only thing that matters is the support of the various matrices, not the
numerical values themselves. Since $B^{a}_{ss^{\prime}}>0$ iff
$M^{a}_{ss^{\prime}}>0$ (either assuming $M^{+}_{ss^{\prime}}>0$ or regarding
$\bot>0$ as False), and similarly for higher orders, we can replace $M^{a}$ by
$B^{a}$ in (iii), and get $B^{aa^{\prime}...a^{i}}_{ss^{i+1}}>0$ iff
$[B^{a}B^{a^{\prime}}...B^{a^{i}}]_{ss^{i+1}}>0$. We could also replace
$M^{a}$ by
$G^{a}_{ss^{\prime}}:=[\mskip-3.0mu[B^{a}_{ss^{\prime}}>0]\mskip-3.0mu]$, then
$[G^{a}G^{a^{\prime}}...G^{a^{i}}]_{ss^{i+1}}>0$ counts the number of paths of
length $i$ from $s$ to $s^{i+1}$ via action sequence $aa^{\prime}...a^{i}$,
and hence determines whether $s^{i+1}$ can be reached. Similarly
$(G^{+})^{i}>0$ iff there is _some_ action sequence that can reach $s^{i+1}$
from $s$. An action $a$ such that $G^{a}(G^{+})^{i}>0$ can be chosen as the
first action of such a sequence if it exists, and
$a^{\prime},a^{\prime\prime}...$ can be found the same way by recursion. So
this deterministic planning/reachability problem has a “unique” solution,
which can be found in time $O(i⋅d⋅(d+k))$ (for fixed $s$ and $s^{i+1}$).

###  $B$ is deterministic.

Assume
$M^{a}_{ss^{\prime}}/M^{+}_{ss^{\prime}}=:B^{a}_{ss^{\prime}}∈\\{0,1,\bot\\}$.
This is true if and only if $M^{a}$ has disjoint support for different $a$,
i.e. iff $M^{a}⊙M^{b}=0~{}∀a≠b$. This in turn means that
$B^{a}_{ss^{\prime}}=[\mskip-3.0mu[W^{a}_{ss^{\prime}}>0]\mskip-3.0mu]$ for
any and only those $W$ with same support as $M$, and hence also
$W^{a}⊙W^{b}=0~{}∀a≠b$, which is another failure case of (i). Here we have
included the case where _no_ action leads from $s$ to $s^{\prime}$, in which
case $W^{+}_{ss^{\prime}}=0$ and $B^{a}$ is undefined ($\bot$). This readily
extends to higher orders: If $B^{aa^{\prime}...}∈\\{0,1,\bot\\}$, then
$B^{aa^{\prime}...}=[\mskip-3.0mu[W^{a}W^{a^{\prime}...}⊘(W^{+})^{i}>0]\mskip-3.0mu]$
iff $W^{a}W^{a^{\prime}}...$ has the same support as $M^{a}M^{a^{\prime}}...$
and

$\displaystyle
W^{a}W^{a^{\prime}}...W^{a^{i}}⊙W^{b}W^{b^{\prime}}...W^{b^{i}}=0~{}~{}~{}∀aa^{\prime}...a^{i}≠bb^{\prime}...b^{i}$
(17)

Note that $W^{a}⊙W^{b}=0$ does not necessarily imply (17), e.g. for
$W^{0}={\textstyle\frac{1}{2}}{1~{}0\choose 0~{}1}$ and
$W^{1}={\textstyle\frac{1}{2}}{0~{}1\choose 1~{}0}$,
$(W^{0})^{2}=(W^{1})^{2}$. In Appendices H&I&J we construct $W$ such that (17)
holds for larger $i$.

##  Appendix E Characterizing $M$ and $W$ for which EqIM(1) holds

$\displaystyle
M^{a}⊘M^{+}=W^{a}⊘W^{+}~{}~{}~{}~{}\Longleftrightarrow~{}~{}~{}~{}W^{a}=M^{a}⊙J~{}~{}~{}~{}\text{with}~{}~{}~{}~{}J:=W^{+}⊘M^{+}$

That is, $J$ is independent of $a$. Phrased differently

For any $M$ and $W$, EqIM(1) is satisfied iff $W^{a}⊘M^{a}$ is independent
$a$. (18)

For a given $M$, this allows to determine all $W$ consistent with EqIM(1), by
just multiplying with any $a$-independent $J≥0$. Not all $J$ though lead to
$W$ consistent with (7). In order to also satisfy (7), $J$ needs to be
restricted as follows: With $Δ_{ss^{\prime}}:=J_{ss^{\prime}}-1$, (7) becomes

$\displaystyle
0~{}\stackrel{{\scriptstyle!}}{{=}}~{}W^{a}_{s+}-M^{a}_{s+}~{}=~{}∑_{s^{\prime}}M^{a}_{ss^{\prime}}(Δ_{ss^{\prime}}+1)-M^{a}_{s+}~{}=~{}∑_{s^{\prime}}M^{a}_{ss^{\prime}}Δ_{ss^{\prime}}$
(19)

For each fixed $s$, these are $k$ homogenous linear equations (one for each
$a$) in $d$ variables. Given $M$, all and only the $W$ consistent with EqIM(1)
_and_ (7) can be obtained via $W^{a}=M^{a}⊙(1+Δ)$ with $Δ$ satisfying
$M^{⋅}_{s⋅}Δ_{s⋅}=0$.

As a special case, $Δ=0$ necessarily if and only if the rank of $M^{⋅}_{s⋅}$
is $≥d$ for every $s$. This gives the precise conditions as stated in
Proposition 1 under which $(i)$ is true. We will next show that EqIM(2)
removes this limitation.

##  Appendix F Characterizing $M$ and $W$ for which EqIM(1) and EqIM(2+) hold

From Appendix E we know that the most general Ansatz for $W^{a}$ satisfying
EqIM(1) is $M^{a}⊙(1+Δ)$. Plugging this into (14) and expanding in Δ, we get

$\displaystyle 0~{}$
$\displaystyle=~{}M^{a}M^{+}⊙(M^{+})^{2}-M^{a}M^{+}⊙(M^{+})^{2}$
$\displaystyle~{}~{}+~{}~{}M^{a}M^{+}⊙[M^{+}(M^{+}⊙Δ)+(M^{+}⊙Δ)⊙M^{+}]-[(M^{a}⊙Δ)M^{+}M^{a}(M^{+}⊙Δ)]⊙(M^{+})^{2}]$
$\displaystyle~{}~{}+~{}~{}M^{a}M^{+}⊙(M^{+}⊙Δ)^{2}-(M^{a}⊙Δ)(M^{+}⊙Δ)⊙(M^{+})^{2}$

This is a collection of quadratic equations in $Δ$. The $Δ$-independent first
line is $0$. We can write this in canonical form:

$\displaystyle
Σ_{kl}A^{a}_{ss^{\prime\prime},kl}Δ_{kl}~{}=~{}R^{a}_{kl}(Δ)~{}~{}~{}\text{with}~{}~{}~{}$
(20) $\displaystyle
A^{a}_{ss^{\prime\prime},kl}~{}:=~{}(Σ_{s^{\prime}}M^{a}_{ss^{\prime}}M^{+}_{s^{\prime}s^{\prime\prime}})(M^{+}_{sk}M^{+}_{ks^{\prime\prime}}δ_{ls^{\prime\prime}}+M^{+}_{sl}M^{+}_{ls^{\prime\prime}}δ_{sk}-M^{a}_{sk}M^{+}_{ks^{\prime\prime}}δ_{ls^{\prime\prime}}-M^{a}_{sl}M^{+}_{ls^{\prime\prime}}δ_{sk})$
$\displaystyle
R^{a}(Δ)~{}:=~{}(M^{a}⊙Δ)(M^{+}⊙Δ)⊙(M^{+})^{2}-M^{a}M^{+}⊙(M^{+}⊙Δ)^{2}$

Let us consider $A^{a}$ as a $d^{2}×d^{2}$ matrix for each $a$, $Δ$ as a
vector of length $d^{2}$, and (wrongly) presume $R^{a}≡0$ at first. $A^{a}$ is
a sum of 4 terms. The second and fourth terms are block-diagonal matrices ($d$
blocks of size $d×d$ in the diagonal) due to the $δ_{sk}$. The first and third
terms are scrambled block-diagonal matrices due to the
$δ_{ls^{\prime\prime}}$, or more precisely, consist of $d×d$ blocks, each bock
being a $d×d$ diagonal matrix. If $M^{a}$ has full rank, each of the four
terms has full rank $d^{2}$, but $A^{a}$ itself can have lower rank,
$0$-eigenvalues due to some cancellations. Random $M$ apparently achieves the
highest rank, but even then, $A^{a}$ itself has only rank $d(d-1)$.

Actually, $A^{a}Δ=0$ is required to hold for all $a$, so the rank of $A$ as a
$kd^{2}×d^{2}$ matrix may still be $d^{2}$. But $A^{+}≡0$ for $k=2$ implies
$A^{0}=-A^{1}$, hence the rank is still at most $d(d-1)$. $k>2$ may rectify
this, but there is an alternative, which works for all $a$: $Δ$ also needs to
satisfy $\eqref{eq:Dconstr}$, which can be rewritten as

$\displaystyle∑_{kl}C^{a}_{s,kl}Δ_{kl}~{}=~{}0~{}~{}~{}\text{with}~{}~{}~{}C^{a}_{s,kl}~{}:=~{}M^{a}_{sl}δ_{sk}$
(21)

These give another $kd$ constraints, and apparently often $d$ new ones from
random $M$. If we combine $A^{\prime}:={A^{⋅}\choose C^{⋅}}$, this implies
that $A^{\prime}$ has often rank $d^{2}$, so $A^{\prime}Δ=0$ can only be
satisfied for $Δ=0$. For $k=2$, $A^{+}=0$, so inclusion of either $A^{0}$ or
$A^{1}$ in $A^{\prime}$ would suffice, but $C^{0}$ and $C^{1}$ are potentially
independent, so both have to be included.

Let us now return to the real case of $R^{a}≠0$ for full random $M$, hence
full-rank $A^{\prime}$. With $R^{\prime}:={R^{⋅}\choose 0}$, we need to solve
$A^{\prime}Δ=R^{\prime}$. Note that $R^{\prime}=R^{\prime}(Δ)$ is not a
constant, but a (homogenous) quadratic function of $Δ$ itself. Consider any
$Δ=Θ(ε)$, then $A^{\prime}Δ=Θ(ε)$ while $R^{\prime}(Δ)=Θ(ε^{2})$, which is a
contradiction for sufficiently small $ε$ (this argument can be made rigorous).
This implies that no $Δ$ with $0<||Δ||<ε$ can satisfy
$A^{\prime}Δ=R^{\prime}(Δ)$. In conclusion,

######  Proposition 4 (Random $M$ and full-rank $A^{\prime}$)

If $A^{\prime}$ has full rank and $W$ is close to $M$, then EqIM(1) and
EqIM(2) imply $W=M$.
Empirically $A^{\prime}$ has full rank for random $M$.

This of course implies $\text{EqIM}(i)∀i$ and also (iv). Globally, i.e. if $W$
is not close to $M$, these implications may not hold.

We have yet to establish sufficient conditions which $M^{a}$ lead to full-rank
$A^{\prime}$. Empirically, this has been true for random $M^{a}$, so should
hold almost surely if $M$ are sampled uniformly. One might conjecture that
full-rank $M^{a}$ are sufficient, but this is not the case. For instance, if
$M^{a}$ is independent $a$, then $A^{\prime}≡0$.

### Zero $A$ and $R$ for full-rank $\dot{M}^{a}$.

We finally we note that $A$ and $R$ can have low rank, indeed $A≡0≡R$ even for
$a$-dependent full-rank $M^{a}$: Consider the example $\dot{M}^{a}$ from (J)
or its generalization (J): First, if for two matrices $M^{a}$ and
$M^{a^{\prime}}$ only one $s^{\prime}$ (depending on $s$ and
$s^{\prime\prime}$) contributes to the sum in $M^{a}M^{a^{\prime}}$ then
$(M^{a}⊙J)(M^{a^{\prime}}⊙J)=M^{a}M^{b}⊙K$ for some $K$. This makes (22) valid
for $M^{a}:=\dot{M}^{a}$ and $W^{a}:=\dot{M}^{a}⊙J$ for any $J$, since for
$aa^{\prime}≠bb^{\prime}$ both sides are 0 by construction of $\dot{M}^{a}$
(the $⊙K$ does nothing to it), and are trivially equal for
$aa^{\prime}=bb^{\prime}$. By summing over $a^{\prime}bb^{\prime}$, also (14)
is valid for any $J$, hence of course also for $J=1+Δ$ for any $Δ$. Since (20)
is equivalent to (14), (20) holds for any $Δ$. This can only be true for $A≡0$
and $R≡0$. This degeneracy in itself does not violate (ii), since the
probability constraints require $W=M$, as established earlier.

##  Appendix G EqIM(1)$∧$EqIM(2+)$\not\rightarrow$EqIM(3) for full low rank
$M$?

The following numerical approach may lead to counter-examples with full
support to (v) without any divisions by $0$ ($M^{+}_{ss^{\prime}}>0$ and
$W^{+}_{ss^{\prime}}>0$ $∀ss^{\prime}$). We now consider full $M^{a}$ but of
rank $r<d$. The most interesting case is where all $M^{a}$ span the same row-
space, i.e. $M^{a}=L^{a}⋅R$, where $L^{a}$ are $d×r$ matrices and $R$ is a
$r×d$ matrix. Recall $A^{\prime}:={A^{⋅}\choose C^{⋅}}$ with $A^{a}$ and
$C^{a}$ defined in (20) and (21). Empirically, for $k=2$, the rank of
$A^{\prime}$ typically is $\min\\{d^{2},(3r-1)d-r(r-1)\\}$, never more, and
only in degenerate cases less. Hence for $r=2$, $A^{\prime}$ is singular for
$d≥5$. Hence for $d≥5$, there exist $Δ≠0$ with $A^{\prime}Δ=0$,

For $Δ_{0}:=Δ=Θ(ε)$, this is an approximate $Θ(ε^{2})$ solution of
$A^{\prime}Δ=R^{\prime}(Δ)$. By iterating $Δ\leftarrow
Δ_{0}+A^{{}^{\prime}+}R^{\prime}(Δ)$, where $A^{{}^{\prime}+}$ is the pseudo-
inverse of $A^{\prime}$, we get an $Θ(ε^{i})$-approximation after $i-2$
iterations. This should rapidly converge to an “exact” non-zero(!) solution
$A^{\prime}Δ=R^{\prime}(Δ)$. This would show that (ii) can fail for full $M$.
Generically, this solution also violates EqIM(3), i.e. also (vi) can fail. By
this we mean, for randomly sampled $L^{a}$ and $R$ (for $a=r=2$ and $d≥5$) and
performing the procedure above, EqIM(3) does not hold. There is a caveat with
this argument, namely if $R^{\prime}$ is not in the range of $A^{\prime}$,
then this construction fails.

##  Appendix H EqIM(1) does not imply EqIM(2) ($⊙$-version)

We have already given a simple example that violates (v) in Section 3, but the
example and methodology provided here generalizes to (vi) and even larger $i$.
We consider deterministic reversible forward dynamics for any policy
$π(a|s)>0~{}∀as$. For simplicity we assume $k=2$ and uniform policy
$π(a|s)={\textstyle\frac{1}{2}}$. We defer a discussion of $0/0$ to the end of
the next Appendix.

We consider $M^{a}$ and $W^{a}$ that permute states. That is,
$M^{⋅}_{ss^{\prime}}:=[\mskip-3.0mu[s^{\prime}=π^{⋅}(s)]\mskip-3.0mu]$ and
$W^{⋅}_{ss^{\prime}}:=[\mskip-3.0mu[s^{\prime}=σ^{⋅}(s)]\mskip-3.0mu]$ for
some permutations $π^{⋅},σ^{⋅}:\\{1,...,d\\}\textrightarrow\\{1,...,d\\}$.
Strictly speaking, we should multiply this by
$π(a|s)={\textstyle\frac{1}{k}}$, but this global factor plays no role here,
so will be dropped everywhere. Matrix multiplication corresponds to
permutation composition:
$[M^{⋅}W^{⋅}]_{ss^{\prime\prime}}=[\mskip-3.0mu[s^{\prime\prime}=σ^{⋅}(π^{.}(s)]\mskip-3.0mu]$.
We denote example permutation (matrices) by $[π]=[π(1)...π(d)]$.

We now construct a counter-example for (v): For $d=4$, let
$M^{0}=W^{0}=\text{Id}=[1234]$ be the identity matrix/permutation. Let
$W^{1}=[2341]$ be the cyclic permutation $1\textrightarrow 2\textrightarrow
3\textrightarrow 4\textrightarrow 1$, and $M^{1}=[2143]$ the cycle pair
$1\leftrightarrow 2$ and $3\leftrightarrow 4$. We know from (18) that EqIM(1)
holds _iff_ $W^{a}⊘M^{a}$ is independent $a$ $(=J)$ _iff_
$W^{a}⊘M^{a}=W^{b}⊘M^{b}~{}∀a,b∈\\{0,1\\}$ _iff_ $W^{a}⊙M^{b}=M^{a}⊙W^{b}$.
Case $a=b$ is trivial, so only $W^{0}⊙M^{1}=M^{0}⊙W^{1}$ needs to be verified.
Now $M^{⋅}⊙W^{⋅}$ of two permutations matrices is _not_ a permutation matrix
(unless $M^{⋅}=W^{⋅}$). It still a 0-1 matrix with at most one non-zero entry
in each row and column. We can generalize the permutation notation to “sub-
permutations” by defining $π(s)=∅$ if row $s$ is empty. For instance
$M^{1}⊙W^{1}=[2∅4∅]$. EqIM(1) holds, since $W^{0}⊙M^{1}=[∅∅∅∅]=M^{0}⊙W^{1}$.

Similarly EqIM($2a$) holds _iff_ $W^{a}W^{a^{\prime}}⊘M^{a}M^{a^{\prime}}$ is
independent $a,a^{\prime}$ _iff_

$\displaystyle
W^{a}W^{a^{\prime}}⊙M^{b}M^{b^{\prime}}~{}=~{}M^{a}M^{a^{\prime}}⊙W^{b}W^{b^{\prime}}~{}~{}~{}∀a,a^{\prime},b,b^{\prime}.$
(22)

But for $a=a^{\prime}=0$ and $b=b^{\prime}=1$ we have

$\displaystyle(W^{0})^{2}⊙(M^{1})^{2}~{}=~{}[1234]⊙[1234]~{}=~{}[1234]~{}≠~{}[∅∅∅∅]~{}=~{}[1234]⊙[3412]~{}=~{}(M^{0})^{2}⊙(W^{1})^{2}$

hence EqIM(1) does not necessarily imply EqIM(2). The advantage of formulation
(22) over (8) is that matrix sums $M^{+}$ and $W^{+}$ are more complicated
objects than the sub-permutation matrices (22). Like random matrices,
permutation matrices, have full rank, but unlike random matrices they can
violate (ii), (iv), and (vi).

##  Appendix I $\text{EqIM}(1a)∧...∧\text{EqIM}(ia)$ do not imply EqIM($i+1$)
($⊙$-version)

Counting variables and equations made the possibility of violating (v) for
$k<d$ plausible (cf. positive result for $k≥d$). A similar counting argument
indicates that (vi) and higher $i$ analogues might actually hold.
Unfortunately this is not the case. I.e. even providing inverse models for all
action sequences up to length $i$ is not sufficient to always uniquely
determine the probability of longer action sequences. This is true even for
deterministic reversible forward dynamics for any policy $π(a|s)>0~{}∀as$. As
for $i=1$, we assume $k=2$, $π(a|s)={\textstyle\frac{1}{2}}$, gloss over
$0/0$, and don’t normalize $M$ and $W$.

For $i=2$, $M^{0}:=W^{0}:=\text{Id}=[123456]$ and $W^{1}:=[234561]=:σ$ ($σ$
for ‘cycle’) and $M^{1}:=[231564]=:π$ can be shown to satisfy EqIM(1) and
EqIM($2a$) but violate EqIM(3). The calculations are not to onerous, but lets
consider directly the general $i$ case: Consider even $d=:2d^{\prime}$ and
identity and cycle (pair)

$\displaystyle M^{0}~{}$
$\displaystyle=~{}W^{0}~{}=~{}\text{Id}~{}=~{}[1,2,...,d-1,d],~{}~{}~{}$
$\displaystyle W^{1}~{}$
$\displaystyle=~{}[2,3,...,d,1],~{}~{}~{}M^{1}~{}=~{}[2,3,...,d^{\prime},1,d^{\prime}+2,...d-1,d,d^{\prime}+1]$

EqIM($ia$) holds _iff_
$W^{a}W^{a^{\prime}}...⊘M^{a}M^{a^{\prime}}...=W^{+}W^{+}...⊘M^{+}M^{+}...$ is
independent $aa^{\prime}...$ _iff_

$\displaystyle
W^{a}W^{a^{\prime}}...W^{a^{i}}⊙M^{b}M^{b^{\prime}}...M^{b^{i}}~{}=~{}M^{a}M^{a^{\prime}}...M^{a^{i}}⊙W^{b}W^{b^{\prime}}...W^{b^{i}}~{}~{}~{}∀aa^{\prime}...a^{i},bb^{\prime}...b^{i}$
(23)

(While this looks like $k^{2i}$ matrix equations, by chaining, checking
$k^{i}$ pairs suffices, which is the same number as in EqIM($ia$)). Now
$W^{a}W^{a^{\prime}}...W^{a^{i}}$ consists of only two types of matrices, a
cycle for $W^{1}=σ$ and identity $W^{0}$. The $W^{0}=\text{Id}$ can be
eliminated, leading to $(W^{1})^{a^{+}}$, where
$a^{+}:=a+a^{\prime}+...+a^{i}$. Similarly
$M^{b}M^{b^{\prime}}...M^{b^{i}}=(M^{1})^{b^{+}}$, etc. Hence we only need to
verify

$\displaystyle(W^{1})^{a^{+}}⊙(M^{1})^{b^{+}}~{}=~{}(M^{1})^{a^{+}}⊙(W^{1})^{b^{+}}~{}~{}~{}\text{for}~{}~{}~{}0≤a^{+},b^{+}≤i$
(24) $\displaystyle(W^{1})^{a^{+}}~{}$
$\displaystyle=~{}[{a^{+}}+1,{a^{+}}+2,...,d,1,2,...,{a^{+}}],~{}~{}~{}\text{while}$
$\displaystyle(M^{1})^{b^{+}}~{}$
$\displaystyle=~{}[{b^{+}}+1,...,d^{\prime},1,...,{b^{+}},d^{\prime}+1+{b^{+}},...,d,d^{\prime}+1,...,d^{\prime}+{b^{+}}]$

hence $(W^{1})^{a^{+}}⊙(M^{1})^{b^{+}}=[∅...∅]=0$ for
$0≤{a^{+}}≠{b^{+}}<d^{\prime}$. For $a^{+}=b^{+}$ both sides of (24) are equal
too. Hence if we choose $d^{\prime}=i+1$, (24) and hence EqIM(1)…EqIM($ia$)
are all satisfied. If we choose $d^{\prime}=i$, $a^{+}=d^{\prime}$, $b^{+}=0$,
(24) reduces to

$\displaystyle(W^{1})^{d^{\prime}}⊙(M^{1})^{0}~{}$
$\displaystyle=~{}[d^{\prime}+1,...,d,1,...,d^{\prime}]⊙\text{Id}~{}=~{}0,~{}~{}~{}\text{and}$
$\displaystyle(M^{1})^{d^{\prime}}⊙(W^{1})^{0}~{}$
$\displaystyle=~{}\text{Id}⊙\text{Id}~{}=~{}\text{Id}$

which are of course not equal. Hence EqIM($i$) fails for $d^{\prime}=i$.
Summing over all $a^{\prime}...a^{d^{\prime}}$ and
$b^{\prime}...b^{d^{\prime}}$, and noting that all other terms are $0$ or
cancel, shows that EqIM($i+$) fails too. Together this shows for
$d^{\prime}=i+1$ that EqIM(1)…EqIM($ia$) do not imply any version of
EqIM($i+1)$.

Despite $M^{a}$ having full rank, $A$ and $A^{\prime}$ defined in Appendix F
have very low rank, indicating potentially many more consistent $W$.

A downside of this example is that it strictly only applies to the $⊙$-version
(23). Many entries of $M^{+}$ and $W^{+}$ and powers thereof are $0$, so (8)
contains many divisions by zero. We were not able to extend this example by
mixing in e.g. a uniform matrix as done in the first counter-example to (v).

Many real-world MDPs are sparse. Only a subset $G⊆S×S$ of transitions
$s\textrightarrow s^{\prime}$ is possible. For $(s,s^{\prime})\not∈G$,
$p(s^{\prime}|sa)=0~{}∀a$, or formally
$M^{a}_{ss^{\prime}}=M^{+}_{ss^{\prime}}=0$. In this case, no action causes
$s\textrightarrow s^{\prime}$ and
$p(a|ss^{\prime})=M^{a}_{ss^{\prime}}/M^{+}_{ss^{\prime}}$ being undefined is
actually appropriate. So we could restrict $(s,s^{\prime})$ to $G$ (and
analogously $(s,...,s^{i})$ and $(ss^{i})$ by chaining $G$) in the conditions
and conclusions of the various conjectures. It is then also natural to
restrict the model class to ${\cal
M}:=\\{M^{⋅}:M^{+}_{ss^{\prime}}>0~{}\Leftrightarrow~{}(s,s^{\prime})∈G\\}$.
For unknown $G$, the condition $M,W∈{\cal M}$ then becomes
$M^{+}_{ss^{\prime}}>0~{}\Leftrightarrow~{}W^{+}_{ss^{\prime}}>0$.
Unfortunately the above counter-example does not even satisfy this weaker
condition, but the more complicated example of Appendix J does. See Appendix O
for how to treat 0/0 in practice.

##  Appendix J Non-Uniqueness of Inverse MDP Models for $i≥2$

In Appendices H/I we provided conjectured/unsatisfactory counter-examples to
$\text{EqIM}(1:i)⇒\text{EqIM}(i+1)$. Here we provide a fully satisfactory
counter-example that avoids the “bad” 0/0.

###  $\text{EqIM}(1)$ and $\text{EqIM}(2a)$ do not imply EqIM(3).

Consider two matrices $\dot{M}^{0}$ and $\dot{M}^{1}$ with disjoint support,
i.e. $\dot{M}^{0}⊙\dot{M}^{1}=0$. In this case
$\smash{\dot{M}^{a}⊘\dot{M}^{+}∈\\{0,1,\bot\\}^{\dot{d}×\dot{d}}}$ is a
partial binary matrix with entry undefined ($\bot$) wherever $\dot{M}^{+}=0$
but otherwise 0 wherever $\dot{M}^{a}=0$ and 1 wherever $\dot{M}^{a}>0$. That
is, it is insensitive to the actual (non-zero) values of $\dot{M}^{a}$. A
simple such $\dot{M}$ is $\dot{M}^{0}={1~{}0\choose 0~{}1}$ and
$\dot{M}^{1}={0~{}1\choose 1~{}0}$, ignoring normalization. For now we ignore
$ss^{\prime}$ for which $\dot{M}^{+}_{ss^{\prime}}=0$ and return to this issue
later.

We consider $M^{a}$ and $W^{a}$ that permute states. That is,
$M^{⋅}_{ss^{\prime}}:=[\mskip-3.0mu[s^{\prime}=π^{⋅}(s)]\mskip-3.0mu]$ and
$W^{⋅}_{ss^{\prime}}:=[\mskip-3.0mu[s^{\prime}=σ^{⋅}(s)]\mskip-3.0mu]$ for
some permutations $π^{⋅},σ^{⋅}:\\{1,...,d\\}\textrightarrow\\{1,...,d\\}$.
Strictly speaking, we should multiply this by e.g.
$π(a|s)={\textstyle\frac{1}{k}}$, but this global factor plays no role here,
so will be dropped everywhere. Matrix multiplication corresponds to
permutation composition:
$[M^{⋅}W^{⋅}]_{ss^{\prime\prime}}=[\mskip-3.0mu[s^{\prime\prime}=σ^{⋅}(π^{.}(s)]\mskip-3.0mu]$.
We denote example permutation (matrices) by $[π]=[π(1)...π(d)]$. Consider now

$\displaystyle\dot{M}^{0}\dot{M}^{0}~{}=~{}[123456]$
$\displaystyle\dot{M}^{0}~{}:=~{}[456123]~{}=:~{}[π_{0}]~{}~{}~{}~{}\smash{\raisebox{-4.30554pt}{$\Longrightarrow$}}~{}~{}~{}~{}~{}~{}$
$\displaystyle\dot{M}^{0}\dot{M}^{1}~{}=~{}[564312]$ (25)
$\displaystyle\dot{M}^{1}~{}:=~{}[231645]~{}=:~{}[π_{1}]~{}~{}~{}~{}\phantom{\Longrightarrow}~{}~{}~{}~{}~{}~{}$
$\displaystyle\dot{M}^{1}\dot{M}^{0}~{}=~{}[645231]$
$\displaystyle\dot{M}^{1}\dot{M}^{1}~{}=~{}[312564]$

No column contains the same number twice, hence this not only satisfies
$\dot{M}^{0}⊙\dot{M}^{1}=0$ but also

$\displaystyle\dot{M}^{a}\dot{M}^{a^{\prime}}⊙\dot{M}^{b}\dot{M}^{b^{\prime}}~{}=~{}0~{}~{}~{}\text{unless
$a=b$ and $a^{\prime}=b^{\prime}$}$ (26)

That $6\textrightarrow 5\textrightarrow 4\textrightarrow 6$ is in reverse oder
to $1\textrightarrow 2\textrightarrow 3\textrightarrow 1$ is crucial for
making $\dot{M}^{0}$ and $\dot{M}^{1}$ not commute. Note that (26) remains
valid if each $1$-entry of $\dot{M}^{a}$ is replaced by a different non-zero
scalar, since (26) is purely multiplicative. So if
$\dot{W}^{a}=\dot{M}^{a}⊙\dot{J}$ for some $J>0$, then
$\dot{W}^{a}\dot{W}^{a^{\prime}}=\dot{M}^{a}\dot{M}^{a^{\prime}}⊙K$ for some
$K>0$. Let $\dot{W}^{a}$ be such a matrix. Then
$[\dot{W}^{a}\dot{W}^{a^{\prime}}⊘\dot{W}^{+}\dot{W}^{+}]_{\dot{s}\dot{s}^{\prime\prime}}=1$
if $[\dot{M}^{a}\dot{M}^{a^{\prime}}]_{\dot{s}\dot{s}^{\prime\prime}}>0$ and 0
(or undefined) otherwise, i.e. is independent of the choice of $J$. So such
$\dot{W}≠\dot{M}$ satisfies EqIM($2a$). Unfortunately the probability
constraints $W^{a}_{s+}=1$ require $J^{a}_{ss^{\prime}}=1$ when
$M^{+}_{ss^{\prime}}>0$, and hence $W=M$. But the general idea is sound and
can be made work as follows:

We split one state, e.g. $s=6$ into two states $s=6a$ and $s=6b$. We leave the
permutation structure intact, except that all deterministic transitions into
$s=6$ are split into stochastic transitions to $s=6a$ and $s=6b$, and
transitions from $6a$ and $6b$ will be to the same state as from original $6$.
Condition (26) is still satisfied, so the above argument still goes through,
but now we can choose different stochastic transitions to $s=6a$ and $s=6b$ in
$W$ and $M$.

Finally, we have to show violation of EqIM(3). EqIM($ia$) holds _iff_
$W^{a}W^{a^{\prime}}...⊘M^{a}M^{a^{\prime}}...=W^{+}W^{+}...⊘M^{+}M^{+}...$ is
independent $aa^{\prime}...$ _iff_

$\displaystyle
W^{a}W^{a^{\prime}}...W^{a^{i}}⊙M^{b}M^{b^{\prime}}...M^{b^{i}}~{}=~{}M^{a}M^{a^{\prime}}...M^{a^{i}}⊙W^{b}W^{b^{\prime}}...W^{b^{i}}~{}~{}~{}∀aa^{\prime}...a^{i},bb^{\prime}...b^{i}$
(27)

(While this looks like $k^{2i}$ matrix equations, by chaining, checking
$k^{i}$ pairs suffices, which is the same number of equations as in
EqIM($ia$)).

It is easier to split _every_ state into two states: $s:=(\dot{s},\ddot{s})$
with $\dot{s}∈\\{1,...,6\\}$ as before and splitter $\ddot{s}∈\\{0,1\\}$.
$M^{a}_{ss^{\prime}}:=\dot{M}^{a}_{\dot{s}\dot{s}^{\prime}}\ddot{M}^{a\dot{s}}_{\ddot{s}\ddot{s}^{\prime}}$.
Note that $\ddot{M}$ is flexible enough to expand each 1-entry in
$\dot{M}^{a}$ to a different $2×2$ (stochastic) matrix, while the 0-entries
become $0~{}0\choose 0~{}0$. This flexibility is important: $\ddot{M}$
independent $a$ or independent $\dot{s}$ would not work. Now let us write out

$\displaystyle[M^{a}M^{a^{\prime}}M^{a^{\prime\prime}}]_{ss^{\prime\prime\prime}}~{}=~{}∑_{\dot{s}^{\prime}\dot{s}^{\prime\prime}}\dot{M}^{a}_{\dot{s}\dot{s}^{\prime}}\ddot{M}^{a^{\prime}}_{\dot{s}^{\prime}\dot{s}^{\prime\prime}}\dot{M}^{a^{\prime\prime}}_{\dot{s}^{\prime\prime}\dot{s}^{\prime\prime\prime}}∑_{\ddot{s}^{\prime}\ddot{s}^{\prime\prime}}\dot{M}^{a\dot{s}}_{\ddot{s}\ddot{s}^{\prime}}\ddot{M}^{a^{\prime}\dot{s}^{\prime}}_{\ddot{s}^{\prime}\ddot{s}^{\prime\prime}}\ddot{M}^{a^{\prime\prime}\dot{s}^{\prime\prime}}_{\ddot{s}^{\prime\prime}\ddot{s}^{\prime\prime\prime}}$
(28)

The crucial difference to the $i=2$ case (26) is that now there are difference
permutation sequences leading to the same permutation, for instance
$\dot{M}^{0}\dot{M}^{0}\dot{M}^{1}=\dot{M}^{1}=\dot{M}^{1}\dot{M}^{0}\dot{M}^{0}$.
Let us choose $aa^{\prime}a^{\prime\prime}=001$ and $\dot{s}=1$, then only
$\dot{s}^{\prime}=π_{0}(\dot{s})=4$ and
$\dot{s}^{\prime\prime}=π_{0}(\dot{s}^{\prime})=1$ contribute to the sum and
$\dot{s}^{\prime\prime\prime}=π_{1}(\dot{s}^{\prime\prime})=2$. For this
choice, (28) becomes
$1⋅1⋅1⋅[\ddot{M}^{01}\ddot{M}^{04}\ddot{M}^{11}]_{\ddot{s}\ddot{s}^{\prime\prime\prime}}$.
If we replace $aa^{\prime}a^{\prime\prime}$ in (28) by
$bb^{\prime}b^{\prime\prime}$ and then choose
$bb^{\prime}b^{\prime\prime}=100$ and again $\dot{s}=1$, then only
$\dot{s}^{\prime}=π_{1}(\dot{s})=2$ and
$\dot{s}^{\prime\prime}=π_{0}(\dot{s}^{\prime})=5$ contribute and
$\dot{s}^{\prime\prime\prime}=π_{0}(\dot{s}^{\prime\prime})=2$. For this
choice, (28) becomes
$1⋅1⋅1⋅[\ddot{M}^{11}\ddot{M}^{02}\ddot{M}^{05}]_{\ddot{s}\ddot{s}^{\prime\prime\prime}}$.
We now define
$W^{a}_{ss^{\prime}}:=\dot{M}^{a}_{\dot{s}\dot{s}^{\prime}}\ddot{W}^{a\dot{s}}_{\ddot{s}\ddot{s}^{\prime}}$.
Since $\dot{M}$ remains the same, the same action and state sequences above
lead to the same result for $W$, just with $\ddot{M}$ replaced by $\ddot{W}$.
If we plug the four expressions into (27) (for $i=3$) we get

$\displaystyle\ddot{W}^{01}\ddot{W}^{04}\ddot{W}^{11}⊙\ddot{M}^{11}\ddot{M}^{02}\ddot{M}^{05}~{}=~{}\ddot{M}^{01}\ddot{M}^{04}\ddot{M}^{11}⊙\ddot{W}^{11}\ddot{W}^{02}\ddot{W}^{05}$

Since this expressions involves 10 different $2×2$ stochastic matrices, there
are plenty of choices to make both sides different. If we choose all $2×2$
matrices to have full support, then by construction, $W$ and $M$ have the same
support, hence constitute a proper counter-example to EqIM(3). We now extend
this construction to $i>2$.

###  $\text{EqIM}(1a)∧...∧\text{EqIM}(ia)$ do not imply EqIM($i+1$).

The construction in the previous paragraph generalizes to $i>2$: We need to
find two permutations $\dot{M}^{0}=π_{0}$ and $\dot{M}^{1}=π_{1}$ such that
for each fixed $j≤i$ all possible $2^{j}$ concatenations (products) of these
permutation (matrices) differ in the sense that no $s$ is mapped to the same
$s^{j}$ (they have disjoint support). Since all
$\dot{M}^{a}\dot{M}^{a^{\prime}}...\dot{M}^{a^{j}}∈\\{0,1\\}$, we can write
this condition compactly as

$\displaystyle∑_{aa^{\prime}...a^{j}}\dot{M}^{a}\dot{M}^{a^{\prime}}...\dot{M}^{a^{j}}~{}∈~{}\\{0,1\\}^{d×d}$

By factoring the sum, this is equivalent to
$(\dot{M}^{+})^{j}∈\\{0,1\\}^{d×d}$. Note that $[(\dot{M}^{+})^{j}]_{ss^{i}}$
counts the number of action sequences $aa^{\prime}...a^{j}$ of length $j$ that
lead from $s$ to $s^{i}$. For $j=i+1$, we want this condition to be violated.
So in order to disprove the implication we need to find two permutations
$M^{0}$ and $M^{1}$ such that

$\displaystyle(\dot{M}^{+})^{j}∈\\{0,1\\}^{d×d}~{}~{}~{}∀j≤i~{}~{}~{}\text{but}~{}~{}~{}(\dot{M}^{+})^{i+1}\not∈\\{0,1\\}^{d×d}$
(29)

The rest of the argument is the same as for the $i=2$ case above: creating two
versions $M^{a}$ and $W^{a}$ of $\dot{M}^{a}$ by spitting one or all states
into two, and replacing the 1s by $2×2$ different stochastic matrices. As for
the choice of $\dot{M}^{a}$, for $i=3$ we can choose 3-cycle and 5-cycle

$\displaystyle\dot{M}^{0}~{}$
$\displaystyle=~{}[6,7,8,9,10,11,12,13,14,15,1,2,3,4,5]$
$\displaystyle=~{}(1,6,11)(2,7,12)(3,8,13)(4,9,14)(5,10,15)$ (30)
$\displaystyle\dot{M}^{1}~{}$
$\displaystyle=~{}[2,3,4,5,1,8,9,10,6,7,14,15,11,12,13]$
$\displaystyle=~{}(1,2,3,4,5)(6,8,10,7,9)(11,14,12,15,13)$

where we also provide the more conventional cycle notation in round brackets.
Crucially the 5-cycles have been chosen to not commute with the 3-cycles
($M^{0}M^{1}≠M^{1}M^{0}$). Conditions (29) can easily be verified numerically.
For higher $i$ we need $p$ cycles and $q$ cycles, where $p$ and $q$ are
relative prime and sufficiently large. We need at least $d=p⋅q≥2^{i}$,
otherwise $\dot{M}^{+}\not∈\\{0,1\\}^{d×d}$ by a simple pigeon-hole argument.
To prove $\text{EqIM}(1a)∧...∧\text{EqIM}(ia)\not\Rightarrow\text{EqIM}(i+1)$
in general for arbitrarily large $i$, we need to invoke some group theory.
All-together we have shown that

######  Proposition 5 ((i)-(vi) can fail)

$\text{EqIM}(1a)∧...∧\text{EqIM}(ia)$ do not necessarily imply EqIM($i+1$) for
any $i$. This in turn implies that (i)-(vi) each can fail for some $M^{⋅}$.

##  Appendix K Systems of Quadratic Matrix Equations

A System of Polynomial Equations (SPE) is a set of multivariate polynomial
equations $\text{Poly}_{j}(x,y,z,...)=0$ over $ℝ$ in $n$ variables
$x,y,z,u,v,w,...∈ℝ$ for $j∈\\{1:m\\}$. This class is NP-hard (via a simple
reduction from 1in3SAT, see Section 7). We can recursively replace each
product $xy$ (sum $bu+cv$) in the polynomials by a new variable $z$ ($w$) and
add “polynomial” equation $z=xy$ ($w=bu+cv$). This results in SPEs consisting
of only linear equations with a single $+$ ($bu+cv=w$) and quadratic equations
without any $+$ ($xy=z$), which are still (even with all $a=b=1$ and $x=y=z$)
NP-hard. We call them Simple Systems of Quadratic Equations (Simple SQE). For
the reduction process to actually work we need one further dummy variable and
equation $q=1$ (to reduce $bu+c=w$). Alternatively, with some extra work, we
can reduce any SPE into a Simple SQE asking for a _non-zero_ solution. We will
pursue the latter, since this is closer to our interest (SQE (20) with
solution $Δ\not≡0$). We can even merge the linear and quadratic equations into
a single form $xy=bu+cv$ by choosing $b=1$ and $c=0$ (replacing $xy$ by $w$
and adding $xy=0⋅u+1⋅w$).

We define a System of Polynomial/Quadratic Matrix Equations (SPME/SQME) as a
set of $m$ multivariate (quadratic) polynomials
$\text{Poly}_{j}(Δ,Γ,...|A,B,C,...)=0$ in the (unknown) matrix variables
$Δ,Γ,...$ _and_ the (given) matrix constants (“coefficients”)
$A,B,C,...$.Alternatively, $\text{Poly}_{j}$ might be viewed as generalized
polynomials over a _non_ -commutative matrix ring in the unknowns only. In any
case, note that

$\displaystyle
A⋅Δ⋅A^{\prime}⋅Δ⋅A^{\prime\prime}+B⋅Δ⋅B^{\prime}+C~{}~{}≠~{}~{}(A⋅A^{\prime}⋅A^{\prime\prime})⋅Δ^{2}+(B⋅B^{\prime})⋅Δ+C$

By writing out all matrix operations in terms of their scalar operations, SPME
is of course a sub-class of SPE. SPE is also a sub-class of SPME (choose all
matrices to be $1×1$ matrices), which implies SPME is NP-hard. But we are
interested in NP-hard _small_ subclasses of SPME, so will construct a more
economical embedding: Assume we have a Simple SQE with $n$ variables
$x,y,z,u,v,...$. We place them into $d×d$ matrix $Δ$ ($d≥\sqrt{n}$)
introducing dummy variables for the remaining entries. We can extract variable
$w=Δ_{ss^{\prime}}$ via $w={\rm{e}}^{s\top}⋅Δ⋅{\rm{e}}^{s^{\prime}}$, where
${\rm{e}}^{s}$ is basis vector ($d×1$-matrix)
$({\rm{e}}^{s})_{s^{\prime}1}=δ_{ss^{\prime}}$. If we replace all variables in
the Simple SQE expressions $xy=au+bv$ by such expressions, we get a Simple
SQME with $\text{Poly}_{j}$ equations of the form (dropping $⋅$ as usual)

$\displaystyle a^{j}ΔA^{\prime j}Δa^{\prime\prime j}~{}=~{}b^{j}Δb^{\prime
j}+c^{j}Δc^{\prime j}~{}~{}~{}∀j$ (31)

While these are scalar equations, since the outer matrices are $1×d$ on the
left and $d×1$ on the right, technically they are matrix equations. We could
pad all involved matrices, including the outer ones, with zeros to square
$ℝ^{d×d}$ matrices of the same size (for sufficiently large $d$, and only
polynomial overhead).

We can reduce (31) to just one equation at the cost of making the equations
more complicated as follows: Write each equation $\text{Poly}_{j}=0$ in the
form ${\rm{e}}^{s}⋅\text{Poly}_{j}⋅{\rm{e}}^{s^{\prime}\top}=0$, with a
different $(s,s^{\prime})$-pair for each $j$. These are now “proper” matrix
equations, but with all entries identically $0$ except entry $(s,s^{\prime})$
being $\text{Poly}_{j}$. This allows us to sum all equations without
conflating them into one (complex) matrix equations

$\displaystyle\textstyle∑_{j}A^{j}ΔA^{\prime j}ΔA^{\prime\prime
j}~{}=~{}∑_{j}B^{j}ΔB^{\prime j}+C^{j}ΔC^{\prime j}$ (32)

Another way to combine (31) into one equation is by putting all $M^{j}$ for
all $j$ into one block-diagonal matrix
$\tilde{M}:=\text{Diag}(M^{1},...,M^{m})$ for
$M∈\\{a,A^{\prime},a^{\prime\prime},b,b^{\prime},c,c^{\prime},Δ\\}$. For
$\tilde{Δ}$ we need to ensure that indeed all blocks $Δ^{j}=Δ$ are equal. This
can be done via $\tilde{Π}^{\top}\tilde{Δ}\tilde{Π}=\tilde{Δ}$ for some cyclic
block permutation $\tilde{Π}$. We further need to ensure that the off-diagonal
blocks of $\tilde{Δ}$ are zero. We can zero each block with one equation, but
it seems impossible to zero all with a bounded number of Simple QMEs. We can
modify the decision problem to decide whether specific sparse solutions
$\tilde{Δ}$ exist. Formally, we can introduce element-wise multiplication $⊙$
and allow one equation of the form $\tilde{B}⊙\tilde{Δ}=0$ with $\tilde{B}$
being $0/1$ on the on/off-diagonal blocks. This leads to a Simple SQME with
$⊙$ in 3 equations (dropping the $\sim$)

$\displaystyle
AΔA^{\prime}ΔA^{\prime\prime}~{}=~{}BΔB^{\prime}+CΔC^{\prime},~{}~{}~{}~{}~{}Π^{\top}\mskip-3.0muΔΠ=Δ,~{}~{}~{}~{}~{}B⊙Δ=0$
(33)

######  Proposition 6 (NP-hardness of Simple SQME)

Systems of Polynomial Equations (SPE) can be polynomially reduced to Simple
Systems of Quadratic Matrix Equations (Simple SQME) (31). The number of
equations can be reduced to 1 at the expense of making the equations complex
(32), or to 2 by asking for sparse solutions or by enforcing sparsity via
$B⊙Δ=0$ (33). Since SPE are NP-hard, deciding the existence of non-zero
solutions for all three SQME versions is also NP-hard.

An NP-hardness proof for a Simple SQME with $⊙$ with 3 equations via reduction
from 1in3SAT that looks much closer to the desired form (15) or (34) is given
in Section 7. By a similar reduction, encoding all $n$ variables and their
complement in the diagonal of $Δ=\text{Diag}(x,\bar{x},y,\bar{y},...,)$, one
can also show that solvability of

$\displaystyle
Δ^{2}=Δ,~{}~{}~{}~{}~{}AΔ1=1,~{}~{}~{}~{}~{}\text{Id}⊙Δ=Δ,~{}~{}~{}~{}~{}\text{with}~{}~{}~{}A∈\\{0,1\\}^{m×2n}$

is NP-complete (1 is the all-1 vector, sparse $A$ with 2 or 3 ones in each row
suffice), but not all SPE can be reduced to this form.

######  Open Problem 7 (Are Bounded SPME NP-hard?)

Are Systems of Polynomial Matrix Equations (without $⊙$) of bounded structural
complexity NP-hard? Bounded means, only the definitions of the constant
matrices scale with $d×d$, but the polynomial degrees, number of equations,
and number of matrix operations are bounded.

##  Appendix L Compact Representation of EqIM(2+)

If only $B^{a+}$ (EqIM(2+)) is given, we can sum (15) over $a^{\prime}$. If we
further assume $a=2$ and define $B=B^{0}$ and $A=B^{0+}$ and $W=W^{+}$ and
exploit $B^{+}=B^{++}=1$, this reduces to the elegant quadratic matrix
equation

$\displaystyle A⊙(W⋅W)~{}=~{}(B⊙W)⋅W~{}~{}~{}~{}~{}~{}~{}~{}~{}$ (34)

with constraints as in (15), or even simpler $W_{s+}=1$ if $π$ is unknown.
This is the most pure formulation of the problem we are trying but are unable
to solve we could come up with. For $A$ and $B$ defined via $M$, we know that
(34) has a solution (namely $W=M^{+}$).

We neither know whether there exists an efficient algorithm to find _some_
solution (34), nor to find _the_ solution in case it is unique, nor to decide
whether there exist solutions in case $A$ and $B$ are chosen arbitrarily.

The condition $W_{s+}=1$ can be relaxed to $W_{s+}>0$. If $W_{ss^{\prime}}$ is
a solution of (34), then also $v_{s}^{-1}W_{ss^{\prime}}v_{s^{\prime}}$ for
any $v_{⋅}>0$ (most easily checked via (11)). Every non-negative matrix has a
real non-negative Eigenvector $v$, and $W_{s+}>0$ implies $v_{s}>0$ and
Eigenvalue $λ>0$, hence for
$W^{\text{norm}}_{ss^{\prime}}:=(λv_{s})^{-1}W_{ss^{\prime}}v_{s^{\prime}}$,
we have $W^{\text{norm}}_{s+}=1$.

$B^{a}≥0$ and $B^{+}=1$ iff $B∈[0;1]$ (and $B^{1}=1-B$). $B^{a+}≥0$ and
$B^{++}=1$ iff $A∈[0;1]$ (and $B^{1+}=1-A$). But we can scale back any $A$ and
$B$ by the same $0<λ<1$ to satisfy these without changing (34), i.e. these
extra conditions ($A$ and $B$ bounded by 1) do not make the problem any
simpler.

##  Appendix M Open Problem

We present the most important open problem(s) in their simplest instantiation
and most elegant form, fully self-contained here: Consider matrices
$A,B,W∈[0;1]^{d×d}$ with $d∈ℕ$, tied by the quadratic matrix equation

$\displaystyle
A⊙(W⋅W)~{}=~{}(B⊙W)⋅W~{}~{}~{}~{}~{}\text{and}~{}~{}~{}~{}~{}W_{s+}=1~{}∀s$
(35)

where $⊙$ is element-wise (Hadamard) multiplication and $⋅$ is standard matrix
multiplication. The open problems are as follows: Given $A$ and $B$, are there
efficient algorithms which

  * (a)

decide whether there exists a $W$ satisfying (35)?

  * (b)

decide whether the solution is unique, assuming (35) has a solution?

  * (c)

compute _a_ solution, assuming (35) has a solution?

  * (d)

compute _the_ solution, assuming (35) has a unique solution?

Computing a real number means, given any $ε>0$, computing an
$ε$-approximation. Efficient means running time is polynomial in $d$, ideally
with a degree independent of $1/ε$. General systems of quadratic equations are
known to be NP-hard, but we do not know the complexity of this particular
matrix sub-class.

The upper bounds $A,B,W≤1$ can always be satisfied by scaling, hence are
irrelevant. $W_{s+}=1$ can be relaxed to $W_{s+}>0$ except in the uniqueness
questions. If helpful: One may assume $A,B,W$ strictly positive. Also, any
finite ($d$-independent) number of equations of the form
$A^{\prime}⊙(W⋅W)~{}=~{}(B^{\prime}⊙W)⋅W$ with other _general_ matrices
$A^{\prime},B^{\prime}∈[0;1]^{d×d}$ may be added, which further constrain the
solution space.

##  Appendix N Experimental Details

Here we provide further experiments supporting and illustrating the theory. In
Appendix O we show how we numerically dealt with $B=0/0=\bot$. Appendix P
derives the formulas for the plotted solution dimensions.

Figure 4: Reconstructing inverse and forward models from inverse models with
noise injected. Rows, from top to bottom, show reconstructions of
$B^{a}$,$B^{a+}$,$B^{a++}$, and $M^{a}$. Noise increases exponentially across
columns, from left to right, $[0,10^{-6},10^{-5},10^{-4},10^{-3}]$. The
subplot titles show the average KL divergence of the recovered distribution
from the ground truth. Figure 5: Reproduced from [LFLDP21], this ‘half-
cheetah’ environment has been augmented with videos of complex scenes. This
highlights how non-controllable aspects of the environment can be made more
complex without changing the underlying control problem. The fact that such
environments are of interest motivates our focus on the Tensor-product special
case.

### Experiments illustrating robustness to noise.

As mentioned in the main text, rather than committing to a specific learning
algorithm, we instead directly inject noise into the true inverse
distributions. This is done by adding $ε×10^{c}$ to the true distribution and
renormalizing $B$, where $ε$ is drawn from the unit uniform distribution:
$ε\sim\mathcal{U}[0,1]$). In Figure 3, this noise is evaluated across several
orders of magnitude ($c$ varied $-7$ to $0$).

The main text also mentions that the effect of this noise is substantially
diminished as the horizon of the inverse model is increased (from $B1:=B^{a}$
to $B3:=B^{a++}$). Figure 4 buttresses this interpretation by showing that the
recovered $B^{a++}$ is qualitatively similar to the ground truth even with
substantial noise.

### Experiments on the Tensor-product special case.

As mentioned in the main text, if $M$ factors into two processes
$\dot{M}^{a}⊗\ddot{M}$, where $\ddot{M}$ is action-independent, then only the
complexity of the action-dependent process $\dot{M}^{a}$ matters for all of
our questions.

This particular special case is important because of its frequency in applied
work. Many environments have most of their complexity in sub-spaces that the
agent has no control over. This is illustrated by Figure 5, reproduced from
[LFLDP21], wherein naturalistic videos are superimposed on relatively simple
continuous control environments. Clearly, the background dynamics can be
arbitrarily complex without impacting the underlying control problem.

We can construct small environments of this form via a simple procedure. We
construct $\dot{M}^{a}$ with $\dot{d}$ states and $k$ actions by sampling each
element of the appropriately sized matrices from $\mathcal{U}[0,1]$ and then
normalizing. $\ddot{M}$ has $\ddot{d}=2$ states that transition uniformly
regardless of the action. For the results shown in Figure 3, $k=5$ as in the
main text, and $d=2\dot{d}$ is varied from 2 to 32.

Note that in Figure 3, the solution dimension is non-zero even when
$\dot{d}≤k<d=2\dot{d}$ (here $k=5$, hence for $d=6|8|10$) despite there
necessarily being a unique solution as per Section 4. This is due to the fact
that the algorithm does not exploit knowledge of the fact that $M$ is a tensor
product, resulting in the solution dimension being correct for the more
general case where $W$ is not confined to being tensor product.

##  Appendix O How to Deal with 0/0

If for some pair of states $(s,s^{\prime})$, no action $a$ of positive
$π$-probability leads from state $s$ to $s^{\prime}$, i.e. if
$M^{+}_{ss^{\prime}}=0$, then $B^{+}_{ss^{\prime}}$ and
$B^{a}_{ss^{\prime}}∀a$ are $0/0=\bot=$ undefined. To also handle
$B^{⋅}_{ss^{\prime}}=\bot$, we need to adapt the linear algorithm in Section
4. We provide 2 different ways of doing so, with a couple of variations, all
leading to the same correct result.

We have to restrict the sum in
$∑_{s^{\prime}}B^{a}_{ss^{\prime}}J_{ss^{\prime}}=π(a|s)$ to those
$s^{\prime}$ for which $B^{a}_{ss^{\prime}}$ is defined. We then solve for
$J_{ss^{\prime}}$, again for $s^{\prime}$ for which $B^{a}_{ss^{\prime}}$ is
defined, and set $J_{ss^{\prime}}=0$ for those $s^{\prime}$ for which
$B^{a}_{ss^{\prime}}=\bot$. Technically this can be achieved by removing the
$s^{\prime}$ columns from matrix $B^{⋅}_{s⋅}$ and $J_{⋅⋅}$ for which
$B^{⋅}_{ss^{\prime}}=\bot$, solve the reduced linear equation system, and
finally reinsert $J_{ss^{\prime}}=0$ for the removed $s^{\prime}$. Simpler is
to replace $B^{a}_{ss^{\prime}}=\bot$ by $B^{a}_{ss^{\prime}}=0$, solve the
equation for $J$, and then set $J_{ss^{\prime}}=0$ for the $s^{\prime}$ for
which the original $B^{a}_{ss^{\prime}}$ was $\bot$. Some solvers
automatically result in $J_{ss^{\prime}}=0$, since this is the minimum norm
solution, but it is better not to reply on this. Instead of setting
$J_{ss^{\prime}}=0$ after solving the linear system, one could also augment
$B^{⋅}_{s⋅}$ with extra rows that enforce $J_{ss^{\prime}}=0$.

Alternatively, we could replace $B^{⋅}_{ss^{\prime}}=\bot$ by a random vector
which sums to 1, e.g. $B^{a}_{ss^{\prime}}=r_{a}/r_{+}$, where $r_{a}=-\log
u_{a}$ with $u_{a}\sim$Uniform$[0;1]$. Provided that the solution is unique,
this also leads to the correct solution (almost surely), and in this way
$J_{ss^{\prime}}=0$ automatically. If the solution is not unique, $W^{⋅}$ will
still satisfy $B^{a}=W^{a}⊘W^{+}$ when for $B^{a}_{ss^{\prime}}≠\bot$, but
$W^{⋅}_{ss^{\prime}}$ may not be $0$.

The adaptation of the Linear Relaxation Algorithm in Section 5 follows the
same pattern: $A^{⋅⋅}_{ss^{i}s^{j}}=\bot$ in (12), whenever one of the three
involved $B$’s is undefined. For such $ss^{i}s^{j}$, we need to ensure that
$\hat{U}_{ss^{i}s^{j}}=0$, which can be done with any of the variations
described above. Once we have $\hat{U}_{ss^{i}s^{j}}$, we set
$C^{a^{i}}_{ss^{i}s^{j}}=0$ if $B^{a^{i}}_{s^{i}s^{j}}=\bot$. No further
intervention is needed, since $\hat{U}_{ss^{i}s^{j}}=0$ already.

##  Appendix P Solution Dimensions of $W$ and $B^{aa^{\prime}}$.

In Section 4 we presented an algorithm for inferring $W$ and $B^{aa^{\prime}}$
from $B^{a}$. Even if $M$ cannot uniquely be reconstructed $\neg$(i),
$B^{aa^{\prime}}$ may still be unique (iii). More generally, the solutions $J$
and $W^{a}$ form linear spaces of dimension $d_{J}=d_{W}≤d(d-1)$
($d_{J}≥d_{W}$ since $W^{⋅}$ is a linear function of $J$ and $d_{J}≤d_{W}$
since $W^{+}=J$). $B^{aa^{\prime}}$ is a (non-linear, polynomial) variety of
dimension $d_{B}≤d_{W}$ at regular points (it is a smooth function of $W$).

### Parameterizing the solutions for $J$ and $W$ and $B$.

We can determine the solution dimensions $d_{J}$, $d_{W}$, and $d_{B}$ as
follows: Let $Γ_{ss^{\prime}}$ be a solution of $[B^{a}⊙Γ]_{s+}=0$. If
$J_{ss^{\prime}}$ is a solution of $[B^{a}⊙J]_{s+}=π(a|s)$, then so is
$\bar{J}:=J+Γ$, hence $W^{a}:=M^{a}+Λ^{a}$ is a solution of
$B^{a}=W^{a}⊘W^{+}$ and $W^{a}_{s+}=π(a|s)$, where $M^{a}:=B^{a}⊙J$ and
$Λ^{a}:=B^{a}⊙Γ$.

If we plug in $W^{a}≡M^{a}+Λ^{a}$ into
$\bar{B}^{aa^{\prime}}:=W^{a}W^{a^{\prime}}⊘(MW+)^{2}$, we get the variety of
$\bar{B}^{aa^{\prime}}$ parameterized in terms $Λ^{a}$. If we expand this non-
linear expression up to linear order in $Λ^{a}$, we get after some algebra

$\displaystyle
B^{aa^{\prime}}=[M^{a}M^{a^{\prime}}+M^{a}Λ^{a^{\prime}}\mskip-3.0mu+\mskip-3.0muΛ^{a}M^{a^{\prime}}-(M^{a}M^{a^{\prime}})⊘(M^{+})^{2}⊙(M^{+}Λ^{+}\mskip-3.0mu+\mskip-3.0muΛ^{+}M^{+})]⊘(M^{+})^{2}+O(Λ^{2})$
(36)

The linear part forms a tangent direction on the $\bar{B}^{aa^{\prime}}$
variety at $B^{aa^{\prime}}:=M^{a}M^{a^{\prime}}⊘(M^{+})^{2}$.

### Determining the solution dimensions for $J$ and $W$ and $B$.

Now, for each $s$, let $Γ_{ss^{\prime}}^{r}$ for $r∈\\{1:d_{Js}\\}$ span all
solutions of $[B^{a}⊙Γ]_{s+}=0$, which can easily be determined by SVD:
$d_{Js}$ is the number zero singular values of matrix $B^{⋅}_{s⋅}$:, and
$Γ_{s⋅}^{r}$ the corresponding singular vectors. Then,
$\bar{J}_{ss^{\prime}}=J_{ss^{\prime}}+∑_{r}Γ^{r}_{ss^{\prime}}z_{sr}$ for any
$z∈ℝ^{d_{J}}$ with $d_{J}=∑_{s}d_{Js}$ is a solution of
$[B^{a}⊙J]_{s+}=π(a|s)$.

Similarly,
$W^{a}_{ss^{\prime}}:=M^{a}_{ss^{\prime}}+∑_{r}Λ^{ar}_{ss^{\prime}}z_{sr}$
with $Λ^{ar}:=B^{a}⊙Γ^{r}$ span all solutions consistent with $B^{a}$ and $π$.
The solution dimension is $d_{W}=∑_{s}d_{Ws}$, where for each $s$, $d_{Ws}$ is
the rank of $Λ^{⋅⋅}_{s⋅}$ if interpreted as a $kd×d_{Js}$ matrix in
$as^{\prime}×r$. $d_{Ws}$ may be smaller than $d_{Js}$, since unlike
$Γ_{s⋅}^{r}$, $Λ^{⋅⋅}_{s⋅}$ may not be full rank.

If we plug $Λ^{a}_{ss^{\prime}}=∑_{r}Λ^{ar}_{ss^{\prime}}z_{sr}$ into (36),
after some index manipulation we get

$\displaystyle\bar{B}^{aa^{\prime}}~{}$
$\displaystyle=~{}B^{aa^{\prime}}~{}+~{}∑_{t=1}^{d}∑_{r=1}^{d_{Jt}}C^{aa^{\prime}rt}z_{tr}⊘(M^{+})^{2}⊘(M^{+})^{2}~{}+~{}O(z^{2})~{}~{}~{}\text{with}~{}~{}~{}$
$\displaystyle C^{aa^{\prime}rt}_{ss^{\prime\prime}}~{}$
$\displaystyle:=~{}(M^{a}_{st}Λ^{a^{\prime}r}_{ts^{\prime\prime}}+[Λ^{ar}M^{a^{\prime}}]_{ss^{\prime\prime}}δ_{ts})[(M^{+})^{2}]_{ss^{\prime\prime}}~{}-~{}[M^{a}M^{a^{\prime}}]_{ss^{\prime\prime}}(M^{+}_{st}Λ^{+r}_{ts^{\prime\prime}}+[Λ^{+r}M^{+}]_{ss^{\prime\prime}}δ_{ts})$

$\bar{B}^{aa^{\prime}}(z)$ is a local parametrization of $B$, and if we drop
the $+O(z^{2})$, it parameterizes its tangential hyperplane at
$B^{aa^{\prime}}$. Its dimension $d_{B}$ is the rank of $C$ interpreted as a
$k^{2}d^{2}×d_{J}$ matrix in $aa^{\prime}ss^{\prime\prime}×rt$. Again, $d_{B}$
may be smaller than $d_{W}$, since $C$ may not be full rank.

### Remarks.

For $r∈\\{1:d_{J}\\}$, the columns of matrix $C$ span the tangential space of
“rescaled” variety $\bar{B}^{aa^{\prime}}$ at $B^{aa^{\prime}}$. Again, the
columns may not be linearly independent. If
$[(M^{+})^{2}]_{ss^{\prime\prime}}=0$, then
$B^{aa^{\prime}}_{ss^{\prime\prime}}=\bot~{}∀aa^{\prime}$, hence all such
$ss^{\prime\prime}$ should be ignored in
$C^{aa^{\prime}r}_{ss^{\prime\prime}}$, but since the corresponding rows in
$C$ are $0$, they don’t contribute to the rank anyway. Numerically, we need to
regard all singular values below some threshold as 0. For (to numerical
precision) exact $B$, the threshold can be fairly small ($10^{-13}$ in all our
experiments). For approximate/learned $B$, the threshold needs to be of the
order of the accuracy of $B$.

### Sampling estimate of $d_{B}$.

A simpler, but less elegant, and more fragile method to estimate $d_{B}$ is as
follows: Fix one solution $J$. Add random noise in direction of the null-space
spanned by $Γ^{r}$ so that it stays a solution, i.e. compute
$\bar{J}=J+∑_{r}Γ^{r}z_{⋅r}$ for random $z$, and from this, $W$ and
$\bar{B}^{aa^{\prime}}$ for many such random $J$. The resulting point cloud
spans covers the solution variety $\bar{B}^{aa^{\prime}}$. Various tools could
be used to analyze this point cloud, e.g. determine its dimension. If $z$ is
chosen small, the point cloud concentrates around $B^{aa^{\prime}}$ and forms
a near-linear space, whose dimension $d_{B}$ can easily be determined by PCA.

### Higher-order $B$ and higher $i$.

In the same way we can derive the solution dimensions $d_{B^{...}}$ for
higher-order $B^{...}$. Also, even though we don’t have (yet) an efficient
algorithm for solving EqIM($i$) for $i>1$ if the solution is not unique, we
still can determine the dimension of the solutions (at a particular point
$M$). Algorithmically already covered is the case of $W$ satisfying
EqIM(1)$∧$EqIM(2), whose solution dimension turns out to be $d_{W}-d_{B}$. The
general procedure is to plug $W=M+Λ$ into and linearly expand EqIM($i$) for
$i$ we to hold. Together they form a system of linear equations whose solution
dimension can be determined by SVD as above.