Stochastic Optimal Control and Estimation with Multiplicative and Internal Noise
Authors: Francesco Damiani
Akiyuki Anzai
Jan Drugowitsch
Gregory C DeAngelis
Rubén Moreno-Bote
Venue: NeurIPS 2024 poster (2024)
Keywords: control theory
stochastic optimal control
sensorimotor system
multiplicative and internal noise
motor control
Primary area: neuroscience_and_cognitive_science
TL;DR: We propose a novel iterative algorithm to solve a stochastic optimal control problem under multiplicative and internal noise, outperforming state-of-the-art solutions in the presence of internal noise and before algorithmic convergence.

Abstract: A pivotal brain computation relies on the ability to sustain perception-action loops. Stochastic optimal control theory offers a mathematical framework to explain these processes at the algorithmic level through optimality principles. However, incorporating a realistic noise model of the sensorimotor system — accounting for multiplicative noise in feedback and motor output, as well as internal noise in estimation — makes the problem challenging. Currently, the algorithm that is commonly used is the one proposed in the seminal study in (Todorov, 2005). After discovering some pitfalls in the original derivation, i.e., unbiased estimation does not hold, we improve the algorithm by proposing an efficient gradient descent-based optimization that minimizes the cost-to-go while only imposing linearity of the control law. The optimal solution is obtained by iteratively propagating in closed form the sufficient statistics to compute the expected cost and then minimizing this cost with respect to the filter and control gains. We demonstrate that this approach results in a significantly lower overall cost than current state-of-the-art solutions, particularly in the presence of internal noise, though the improvement is present in other circumstances as well, with theoretical explanations for this enhanced performance. Providing the optimal control law is key for inverse control inference, especially in explaining behavioral data under rationality assumptions.

Stochastic Optimal Control and Estimation with Multiplicative and Internal Noise
Francesco Damiani
Akiyuki Anzai
Center for Brain and Cognition, Department of Engineering
Department of Brain and Cognitive Sciences University of Rochester Rochester, USA
Pompeu Fabra University
Barcelona, ES francesco dami ani @nnfuup.... edu
aanzai@ura rochester....... edu
Jan Drugowitsch Department of Neurobiology
Gregory C. DeAngelis
Department of Brain and Cognitive Sciences University of Rochester Rochester, USA gdeangellssbur.rooco.. rochessse edu
Harvard Medical School
Boston, USA
jan_ www.drrgovittcchhmmms...
harvard.edu hawpeacchhhhnss.aaa.a.eeeeeed
Rubén Moreno-Bote
Center for Brain and Cogmition, Department of Engineering,
Serra Húnter Fellow Programme Pompeu Fabra University Barcelona, ES ruben. moreeonur edu
Abstract
A pivotal brain computation relies on the ability to sustain perception-actio loops. Stochastic optimal control theory offers mathematical framework explain these processes at the algorithmic level through optimality principles. However, incorporating realistic noise model of the sensorimotor system accounting for multiplicative noise feedback and motor output, as well as internal noise in estimation makes the problem challenging. Currently, the algorithm that is commonly used is the one proposed in the seminal study in I. After discovering some pitfalls in the original derivation, i.e., unbiased estimation does not hold, we improve the algorithm by proposing an efficient gradient descent-basee optimization that minimizes the cost-to-go while only imposing linearity of the control law. The optimal solution is obtained by iteratively propagating in closed form the sufficient statistics to compute the expected cost and then minimizing this cost with respect to the filter and control gains. We demonstrate that this approach results significantly lower overall particularly in the presence of internal noise, though the improvement is present other circumstances as well, with theoretical explanations for this enhanced performance. Providing the optimal control law is key for inverse control inference, especially in explaining behavioral data under rationality assumptions.
than current tfate-of-the-ar solutions,
Introduction
The sensorimotor system possesses a remarkable ability reliably execute actions aligned with external and internal goals in spite of the noise sources affecting it 2 3 and the numerous solutions
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

through which the same goals can be achieved 14. Our nervous system is able to combine the complexity of the mechanical properties of the body with a regulatory control system [3]. How such
control is implemented at a computational and algorithmic level is still an open question in systems neuroscience.
Optimal feedback control provides a valuable framework for understanding how the motor system creates coordinated and adaptable behavior 6 []. However, using optimality principles infer the underlying computation [7] is a powerful yet potentially risky approach. Indeed, multiple independent factors can lead to discrepancies with experimental data. For instance, predictions might be based on unsuitable approximations or conditions. Consequently, although optimal feedback control, and more broadly, stochastic optimal control theory, can be regarded as ideal candidates to understand the principles of motor control [5] 6 7 their effectiveness depends on the mathematical correctness of the derived predictions and the accuracy of their assumptions.
Solving an optimal control problem implies deriving an optimal state-to-action policy, or control-law, to minimize certain cost function, usually embedding control effort and task related goals [6]. The classical framework in which analytical solutions for optimal feedback control of stochastic, partially observable, continuous, non stationary, and high-dimensional systems can be derived, the (LQAG) problem (see Section 2 that assumes linear dynamics, a quadratic cost function and additive Gaussian noise 8 9.. Despite having been used in the past to model motor control [10] [1] 12], these assumptions are too limiting to explain a wide range of observed, relevant behaviors, like smooth velocity profiles [ 13 [4], I speed-accuracy trade-offs
[1 15 sensorimotor system is crucial to fill this gap, even at the cost of decreasing the mathematical tractability of the problem [1. Indeed, accounting for control [1 6 17 20 and -ipeen-de
16 17 and movement corrections [18 [9]. Including realistic noise model for the
[
21
23 [4] noise at the motor output and sensory feedback level, and for internal noise 24 in the estimation process, permits explaining broad range of experimentally observed
22
phenomena
18
25 26), as discussed in Section
The seminal study in [. widely regarded as state-of-the the-art for solving optimal control problems under this extended noise model, offers an iterative algorithm whereby stochastic optimal control problem, incorporating multiplicative motor and sensory noise and additive internal noise, can be efficiently solved. Such an algorithm currently used to explain behavioral data in the context Unfortunately, the derivation used in [ erroneously assumes
of inverse optimal control 27 unbiased estimators. We propose an alternative algorithm that addresses this issue by assuming only linear control. The algorithm leverages the fact that the cost function can be computed from closed-form moment expressions. which can then be minimized numerically. To handle potential high computational costs, we derive an analytical counterpart for the optimization, based on the efficient propagation of cost function derivatives over time. For simplicity, the algorithm is derived for simpler, yet relevant, case outlined in Section 3.3 with extensions to the more general auraorithhm
case also discussed. Our algorithm outperforms the solutions in [[ under internal noise, providing both theoretical and heuristic explanations for the performance differences. In sensorimotor hand- reaching task, it reduces the cost by up to 90% when internal noise constitutes 10% of the total. This reveals qualitatively different behaviors, underscoring the importance of using the actual optimal controller, particularly when explaining behavior in a principled way [27].
In Section 2 to fix notation and ideas, we begin by formalizing the optimal control problem using the classic LQAG framework R addressing partial observability and assuming fully additive noise (Section 2.T. We then introduce the Linear- Miwrraaaiiiiiiiiiiiiiiieeeeeeeee
(LQMG) framework, which extends the noise model to include multiplicative noise in both control and observations, as well as additive internal noise, following the approach of [1] (Section 2.2). We demonstrate that the well-established solution from I produces suboptimal solutions in the presence of internal noise and prior to full algorithmic convergence (Section 2.3 In Section 3 we introduce novel numerical algorithm that achieves optimal solutions and outperforms the approach in 11. as demonstrated empirically in Section 3.2 Finally, in Section 3.3 we present the analytical counterpart to the numerical algorithm.
2

---

2 Control and Estimation with Multiplicative and Internal Noise
2.1 The Classic haapn--Quuudric---ddddtt----sSs (LQAG) Problem
Stochastic optimal control theory formalizes the idea of controlling a dynamical system under partial observability to accomplish a goal [8. In the LQQG problem (typically referred to LQG), linear system with latent state 74 E IIm
(1) is controlled by control signal E IP, with time-independent matrices A R and B E R=xp, and initial condition Ii -onsidering time-dependent matrices is straightforward. The term & E I stands for Gaussian random variable with zero mean and covariance Ze (we always consider i.i.d. random variables, but note that temporally correlated random variables can be generated by filtering the noise with the linear dynamics in Eq. D. In the most relevant case, the controller does not have full access to the latent state Tz: the observation Y E R is a noisy version of øt,
Fe+11 = Az, + Bu + &
y/ H + w/
with observation matrix H E R&xm and W, E R$ being Gaussian random variable with zero mean and covariance Pay Note that all noise sources are additive, that is, state independent, hence we refer to the classic LQG problem as LQAG. The controller = we (e(B)) is constrained to be function of the past observations only, yi (y1...., Vt-1). 1). and must be optimized to minimize the total quadratic cost
T
T
E[J Σ Elid = Σ B Ba + +u, R 1=1
(3)
[=1
where the duration of the task, and jt is the cost per step in a trial, which includes control cost (reflecting the internal goal of minimizing control effort) determined by the symmetric positive definite matrix R3 Rexp, with R3 > 0, and a state cost (modeling potential external goals, such as minimizing the distance to a chosen target), determined by Q, E mm Again, Qe is symmetric and positive definite, Qt > O, and modulates the cost of the state being far from a chosen target. J is the total cost, over whole trial, while B[J] is the total expected cost. Here, the expectation B(SO) denotes an average over all noise random variables with the same initial condition, that is, Bif(-
ʃ dæ.....r f p(z2....T V).....))
The optimal controller can be derived analytically [8 (see Appendix In summary, it is linear function of the state estimate ft, U, = L, where L, E RPxm is the control gain, and ê, is the estimator of the unobserved variable Tt, recursively computed with linear Kalman filter
8 = Aft, + Bu, + Kz H..))
(4)
with filter gains K, E R and initial condition _1 E[xr] (r) is random Gaussian variable with covariance Σm;) - to start with an unbiased estimate of the latent variable. Intuitively, the estimate at time L ++ 1 consists of a next-state prediction term (the first two terms in the r.h.s.) from the current estimate è; plus a correction (third term) that depends on the prediction error, the difference between the new observation y/ and the previous state prediction, weighted by its reliability. For example, if the noise magnitude is very large, K, = O, indicating that an open-loop strategy would be optimal
For the classic LQAG problem, it is well known that the optimal Kalman filter satisfies the orthogonality principle [8], stating that the estimation error is orthogonal to the optimal esti- mator St, i.e. [x[ E, ) = O, where we define from here onwards the expectation as B(S(-)] dæs....T diz.... (i2....TI(() p(22....T, .)) where p is the joint density of latent and estimation variables with initial condition én = B[m)] e.g., BB(x)] = d
Also, it is clear from the analytical expression (Appendix I. the computation for the optimal controller and filter gains are mathematically independent of each other, the so-called separation
principle 28 29 30), which is closely related to the concept of certainty equivalence
2.2 An Extended Noise Model: Optimal Control Beyond the LQAG Framework
Purely additive noise sources alone are insufficient to model the sensorimotor action-perception loop, as multiplicative noise affects both motor control [ 15 [6] and sensory feedback, including
3

---

visual and proprioceptive signals 13 21 22 14 For instance, stronger muscle forces produce greater noise [5 17, and visual sensory noise increases in the periphery relative the fovea Accounting for these characteristics is crucial for explaining and reproducing various behavioral features in reaching movements, such as stereotyped bell-shaped velocity profiles and the speed-accuracy trade-off 15 16 17].
1
23
21
22
1 1
14
These considerations result in the Linear- Munudraaie-Multiplcctive---ee- the following dynamics for state and sensory feedback
(LQMG) model, with
I++1 = Axx + But + E +
(5)
yt H + We + Σ D,
(6)
In comparison to Eqs. account for multiplicative motor noise (Eq. (5) and sensory noise (Eq. 6 Specifically, performing
2 for the classic LQAG model, the LQMG model adds the final terms to
control action de introduces noise proportional to the control signal itself (Eq. (5), while perceiving the state variable 2, induces noise proportional to the observed state (Eq. 6p. Here, C. E Rmx
and D; E R*xmm are constant scaling matrices, and & E IR and E IRd are zero-mean Gaussian noise terms with covariances Ωg I and Q2p I. respectively M. For simplicity, in the expressions derived below, we set d = 1 to improve readability, without loss of generality.
As in the LQAG problem, the objective is find the optimal control signal u, E Rp that depends solely on past observations .t- to minimize the cost function defined in Eq. [] In this case, the optimal state estimate and corresponding filters are state-dependent and in general intractable, but can simplify the problem by assuming, as in (1) that the filter is non-adaptive (i.e. independent of the state estimate, [I), similar to the classic LQAG problem. This leads to the assumption that the state estimate follows the equation
41+1 Ad, + Bu, + (1 H + T/
with the same terminology as in Eq. 4 The initial state 21 and its estimate æ assumed to be independent, are Gaussian variables with the same mean B[x,] = E[ and covariances Σ,, and Es,, respectively. This dynamics of the state estimate only differs from Eq. due to the presence of an additional zero-mean Gaussian noise term 7/ E Rn with covariance S2,, which models the possibility of inefficient filtering of the past observations. This noise term may represent internal
fluctuations in neural activity [2] 22 32
3
or inaccuracies in the filtering process, and is important
for explaining behavioral data D.
Under the new LQMG model (Eqs. 7). the task is to find the optimal control signal Uy (f2).
for 1.... T 1. and the filter gains K1,.. -T-2, (T-2, that minimize the quadratic cost Eq.
2.3 State-of-the-Art Solutions for the LQMG Model: Causes of Suboptimality
is important to first recognize that solving the LQMG problem in Eqs. 5 is significantly more complex than in the classic LQAG problem: while in the latter the separation principle applies, allowing for a direct analytical solution (Appendix I. in the former the principle does not hold, resulting in tightly intertwined controller and filter gains. Notably, the presence of internal noise alone introduces control-estimation interdependencies- factor previously overlooked in earlier approaches.
The algorithm currently used to solve the optimal control problem under the LQMG model (Eqs.
is the one introduced in the seminal work of (1, whose solutions are detailed in Appendix A.2
The impact of this research extends beyond theoretical considerations 227 34 35 36
37 38
We now describe some pitfalls in the original derivation and explain why certain assumptions fail,
leading to suboptimality. In Section we propose an alternative algorithm, and in Section 3.2 we
demonstrate that our solutions outperform the previous ones.
The algorithm in | assumes, throughout the derivation of the optimal control-estimation loop, that the estimator is unbiased, meaning El [fx]] it.. However, this condition is never truly satisfied. To illustrate this conceptually, we can consider one-dimensional toy problem involving
partially observable stochastic process 2, (Fig. [T). Assume that at time - 1. the condition

---

G(mz-1(ě(-11 - holds. Now, suppose that at the same time, a large positive fluctuation, possibly caused by sensory or internal noise, affects the agent's internal estimate. As a result, while the actual state changes only slightly compared to 1;-1, the state estimate ã changes significantly, so that dt > 1. this point, it becomes clear that the expected value of x4 conditioned on f, cannot equal 4t, thus violating the unbiasedness condition (see Fig. Tp). This effect is more pronounced when the state estimate undergoes large fluctuations, but a similar bias, although smaller, would still be present with minor fluctuations.
We demonstrate this issue numerically by considering a one dimensional problem (m = =p=k= 1) with multiplicative, additive, and internal noise (for the details see Appendix In the absence of internal noise, the violation of unbiasedness is still present, though becomes pronounced only for large values of é (Fig. [[..c) However, when internal noise is introduced, the bias increases substantially since the internal fluctuations are not directly attenuated by the gains K2 (see also Appendix A.3)
(a)
7,
2. 0.0
Uxl
-10
Figure 1: The invalidity of the unbiasedness condition. (a) A toy example illustrating estimation bias. The black line represents the dynamics of partially observable process, zz, while the red line shows the state estimate, f, biased by random internal fluctuations from the noise term f) (orange arrow). (b-c) Ef [ for / = 8, as function of ã, for o, 0.0, 0.6, respectively, using the solutions from [. The conditional expectation EE [2, [ is computed through Monte Carlo simulations (dots + error bars: mcan + 1std). The gray straight line represents the identity line, where E[ [ ât.
As said above, in the classic LQAG problem the optimal Kalman filter satisfies the orthogonality principle [8) In contrast, when internal noise non-zero we show that the orthogonality principle does not hold anymore for the optimal filter (see Appendix [A.3II) We further show that the condition of unbiasedness implics the orthogonality principle, but not the other way around. Therefore, applying unbiasedness to solve the optimal control problem introduced in Section 2 leads to suboptimal noise. similar
solutions when internal noise is present, irrespective of control or signal-depende issue arises when the algorithm has not yet converged, even for zero internal noise, because the orthogonality principle only holds for optimized filter gains (and the absence of internal noise). As result, the algorithm in [I] also fails to produce the optimal control when derived with fixed suboptimal estimator, and to produce the optimal filter estimate with a fixed suboptimal controller.
3 A Novel Algorithm for Optimal Control Problems
To address the issues outlined in the previous section, we introduce an alternative method for solving the LQMG problem presented in Section 2.2 We compute the expected total accumulated cost, E[J]. by averaging over all stochastic terms present in Eqs. 5 and as a function of L, and K,- For fixed Le and K,, B[J] serves as the objective function for a standard gradient descent algorithm aimed at minimizing it. In Section 3.1 we detail the computation of this objective function through moment propagation and discuss the minimization process with respect to L1.... T1 and K1... T-2- In Section 3.2 we demonstrate that this approach outperforms state-of-the ofttth---tttttt algorithms. In Section we derive the analytical counterpart to our numerical algorithm.
3.3
3.1 Minimization of Theoretical Expected Cost Through Numerical Gradient Descent (GD)
Our method assumes linear control, where the control signal ue is linear in the internal estimate
L.
(8)
This assumption is not very limiting, as it is correct for the classic LQAG problem and has been used before for the LQMG problem []] Crucially, do not assume unbiasedness to solve the
5

---

optimal control problem. The expected total accumulated cost is computed by propagating the first two moments of x and ğ in closed form. Given that both control and estimation are linear in z and å, and the cost function is quadratic in x: and the first and second moments act as sufficient statistics. As result, no additional approximations (e.g., assuming Gaussianity and () are required find the optimal solutions. By using Eq. [ and the formula for the expected cost of quadratic form, Eq.
3 can be rewritten as
T
T Σ (EE + +
B
1=1 Σ
(9)
where Tr[[]] stands for the trace operation, Mx,, is the covariance matrix of the latent state Iy and Ca, the covariance the state estimate at time t. Note that B[z,]. B| Mz, and Eg, will implicitly depend on L......--1 and K -.. From Eqs. 7 we can derive the update equations propagate the first and second-orde moments E E[za] [[ Σz, and Σg, in a closed-form manner, in order to compute the total expected cost B[J] at fixed L -1.--.T-11 and K1.--T -T.------- (for the derivation see 1 for simplicity (the case c, d 1 follows simply by replacing terms with D or C matrices by C2 and D., respectively, and sum over i).
Appendix A.4.I). Here and the following we set c =
To rewrite our results in a more compact form, we define (similarly to 27 39)
Pize
I0 x4 E
14
(10)
12.1
Exer
23
(11)
Li
A
B
M,
(12)
K,H
A BL
K,H
and
CLi(E + Me H, + Sgg
K.D(Ea. + PZ. D + KR + o.)) (13)
G =
where Dg,,,, = [[ Efa, E and Σe,,,, =ET 22. We have defined as a column vector whose block elements are m dimensional vectors. Similarly, E,, M2 and G2 are block matrices, whose block elements
x matrices.
With these definitions, we see that the first and second moments propagate in a closed manner as
Me+1 = Mepu
(14)
Σ M,E.M + G..
(15)
In other words, if the first and second moments are known at time t, their values can be recursively computed at time & + 1, and no other moments involved in the caleulations.
As result, given the initial conditions for μl and Σ1, we can compute the expected accumulated cost
B[J] at fixed L,,... , T for the algorithm to compute the expected cost [E[] is provided in Appendix A.4.2 Algorithm To find the optimal control and filter gains we would then use [[]] as the objective function of
and Ki.... KK....7-2 -- by using Eqs. 14 15 together with Eq. 9 The pseudo-code
numerical gradient descent procedure. The analytical gradient descent counterpart discussed in
Section 3.3
3.2 Experiments: Enhanced Performance with the GD Algorithm
We apply our algorithm and compare it with the state-of-the-art solutions in two scenarios governed by linear dynamical system (Eqs. 8 Hereafter, GD refers to our numerical algorithm (Section B.I. and TOD refers the algorithm from [] First, in simplified
dieemmmniional reaching task (m p k 1) with all noise sources present, we show that for non-zero internal noise, Ω, > 0, GD outperforms TOD, resulting lower accumulated cost. Second, in a reaching task with cour-dimeniional state and one-dimensional control and sensory feedback (n = 4. p k = 1), GD predicts qualitatively different behavior and shows 90% performance improvement when internal noise contributes 10% of the total.

---

One-Dimensional Case: Understanding the Qualitative Differences We examine the case where m p = L, incorporating multiplicative, additive, and internal noise. The system parameters are provided in Table 2 in Appendix A.5.1 (note that we define the strength of the internal noise as on = =VS,,). With non-zero internal noise, our algorithm achieves lower-cost solutions compared to the method proposed in I (Fig. 2h). This improved performance arises from different modulations
of L and K2 as o, varies (Figs. 2p. c). Crucially, our solution results in control gains that decrease as internal noise increases, while the TOD solution shows little sensitivity to internal noise magnitude (Fig. [p): internal noise increasingly intertwines the optimal solutions for K₁ and L, In Appendix A.5.2 we provide a geometric interpretation of why this modulation is optimal, demonstrating that this optimality enhances adaptability, and showing how internal noise disrupts the orthogonality principle.
As outlined in Section 3 the incorrect unbiasedness condition implies the orthogonality principle. Thus, even if the estimator is biased, the algorithm in ( finds the optimal solution with zero internal noise, as this principle holds for the optimal Kalman filter. However, for non-zero internal noise, on > O, the TOD algorithm underperforms due to the breakdown of the orthogonality principle. In Appendix A.5.3 we also discuss that this suboptimality is observed before the algorithm converges, when the optimal control law is derived from fixed suboptimal filters, and vice versa, regardless of the presence of internal noise.
2.00
(0)
(0)
YOU
EPOXP
-0
1000
10 75
(D)
0.5
10.50
(c)
2.0
16 75
1.5
0.5
0.0
0.5
1.0
1.5
2.0
Internal Dolse dy
Figure 2: Enhanced performance and different solutions with internal noise. (a) Expected accumu- lated cost [E] computed by averaging the quantity from Eq. 3over 50k trials, as a function of the internal noise strength oy for TOD II and GD (Section 3.1 al algorithms (mean ISEM from Monte Carlo simulations, error bars not visible as too small). The expected value aligns with our theoretical estimate of E[J]. as derived in Section 3.1 (b-c) Optimal control and filter gains, Le and K,, for TOD and GD and algorithms. We also present the solutions derived from the analytical counterpart of the numerical GD algorithm to demonstrate that they match the optimal solutions ('Fixed Point Optimization with Moments Propagation' FPOMP algorithm, Section 3.3
multi-Dimensional Case: A Motor Control Application To demonstrate the scalability of our algorithm to more realistic motor control scenarios, we examine problem with a four-dimensional state vector (m = 4) and one-dimensional control and sensory feedback (p k 1). The task is identical to that in [I] except that we include internal noise, which was absent in the original formulation. We model a single-joint reaching movement aimed at minimizing the distance between the hand position pt and a target, while minimizing control effort. The state variable of the problem is I [Pt. P. = dps/ / dt. ft, gt], where fe is the force acting the hand and ge is an auxiliary variable used to filter the control signal un (see [ and 440 for a more detailed discussion). We include control and state-dependent noise, as well internal noise, perturbing the estimate of pt. Note that
now we denote og = =V V. All parameters are listed in Appendix A.5.4
Our results confirm the findings from the previous scenario in this more complex sensorimotor task. The GD algorithm achieves a lower expected accumulated cost, with the performance gap widening as internal noise o, increases (Fig. 3p). This achieved by reducing control gains with increasing
o, (Fig. 3D), resulting in smoother control signal on individual trials (Fig. (3)) and overall reduced control effort (Fig. 3d.. These adjustments lead to two key behavioral outcomes: slower movements
compared to TOD solutions and significantly reduced trialto--trial variability (Fig. 3b). As mentioned
7

---

earlier, GD outperforms the algorithm proposed by [ reducing the cost by up to 90% when internal noise contributes approximately 10% of the total noise (a, = 0.05). To further quantify the impact of internal noise in this scenario and enhance clarity, we calculate the ratio between the average fluctuation amplitude of the state estimate (FA, the standard deviation of the state estimate) and the average range of variation of the state (RV, the range of variation in position Pt. defined as max (pt) min, (P:)). The resulting ratio is FA/RV ~ 0.5 for on = 0.05 (see also Appendix (.555)
We emphasize that the GD algorithm naturally handles arbitrarily high-dimensional problems without requiring any further adjustments. Algorithm I in Appendix A.4.2 which serves as the objective function for the numerical optimization via gradient descent, is designed to accommodate arbitrary dimensions for state, control, and sensory feedback, as well as trial duration. However, the time horizon must remain finite by assumption, similar to [[]. We empirically demonstrate the scalability of our algorithm by applying to a significantly higher-dimensional problem, where the linear dynamical system is governed by random matrices with Gaussian entries. Specifically, we consider m = 10, p = 4, and = 10 for the dimensions of state, control, and observation, respectively (Appendix A.5.6) A.5.6).
GD
(c)
TOD
10
GD
-4
10
1d
10"
10
10
25
50
75
100
25
50
75
100
Internal nalse 0g
(d)
El [u()))d
10¹
0.10
0.05
16
10
25
100
internal moiee Oz
Figure 3: Single-joint reaching task. (a) Expected accumulated cost [[.] as a function of the internal noise (mean + ISEM, error bars not visible as too small), for TOD and GD algorithms.
(b) Magnitude of the control gain vector as function of time for TOD and GD solutions. (c) control signal in a sample trial for the two algorithms for o, = 0.05. (d) Amount of control as o, increases, that is, mean integral of the absolute control signal for the two algorithms.
(e) Mean position over time for the two solutions. All averages are over 50k trials; shadowed areas are SEM.
3.3 An Analytical Approach: FPOMP Algorithm
Although, as discussed in the previous section, our GD algorithm performs well for arbitrarily high-dimensional problems (see also Appendix A.5.6). its application to complex, real-world tasks can become computationally expensive. With Li.....T -T Rpxm and K,,...T- .T-2 E Rm the total number of parameters to optimize via numerical Gradient Descent is mp(TT - 1) + mk(T - 2). which can become quite large for nigh-dimensional state spaces. For example, in the ID problem. the TOD and GD algorithms have comparable and short computation times. However, for the multi- munti-ddimesionall sensorimotor task presented above, the GD algorithm takes significantly longer: while the TOD algorithm completes in just a few minutes on a standard laptop. the GD optimization requires several hours (approximately 4 hours).
To address this, we propose an analytically derived iterative algorithm, where we alternate between finding the optimal (i.e., cost-minimizing) ......7-11 and Ki,...,T T 2, denoted as I and Kü for fixed state and state estimate moments, Mt and E,, and re-computing these moments in light of the updated Lt's 's and K,'s s. We refer to this method as the 'Fixed Point Optimization with Moments
8

---

Propagation' (FPOMP) algorithm. Note that, when optimizing for L, and K₁, we condition also on all the future control and filter gains, therefore not only on μ4 and Σ,, but also on 1t-1... . and Kt+,..., .7..---222 To simplify the notation, we will omit this implicit dependence and explicitly state only the dependence on fe and Σ,. Therefore, at each iteration, we identify the critical points of the total conditional expected cost with respect to L and K,, We can compute the conditional expected
2-
cost per step at time + i, i = O, O..,. T as
T Elje+s (e: E.] = BB I14, Ex + 11/ De+ Me, Ex]+
+ E(ūc++ Me, _ L1 Rr - L (2+ Di]+
(16)
+ Tr Pr(Qu++) E 20 mr S.) + Tr T. Ro- Iuti Spe..
where Bf [Me E.], T6] [8 [Mte En], E I and + My are computed by propagating the + i. Indeed, as discussed Section 3.1 pr and Σ, serve as sufficient statistics for computing the expected cost. To derive L; and K;, we set the derivatives of the expected cost in Eq. 9 to zero. Excluding the constant terms, we obtain
moments f4 and Σ, (Eqs. 14 15) until
0
E Bi Σg] =0
(17)
a
T
Eli++slme 2g = 0.
(18)
0k.
As shown in Appendix 1. and A.6.2 solving Eqs. 17 18
leads to a backward algorithm to compute
L, and Kł,
I. = f(1: E,, L 1.441 T- K .T2.
(19)
K-- = g(vr: E1, T T- K ++1 1+1....TT-2 T
(20)
with 1= 1,.., T- for 1. and 1 == 1,..., T 2 for Kū- From this we can build a recursive relationship that, starting from an initial guess for I. T-1 and Ki.... T iteratively computes all the moments M1,--.TT , and Ei,...,T (Eqs. 15 at fixed Li T-I and Ki.... T- , Given those moments. Li....--11 and K 1.... T - attained. The pseudo-code for the FPOMP algorithm, with its implementation details, can found in Appendix [A.6.3] Algorithm? 2 In such way, we eliminate the numerical optimization procedure, making the algorithm suitable extremely large optimal control problems. The FPOMP algorithm is flexible and works for arbitrary dimensions of state, control, sensory feedback, and trial duration, with computational costs and runtime comparable to those of the approach proposed in [
are updated by using Eqs. 19 20 and so on, until convergence is
T-1
In Appendix I. we explicitly solve Eqs. 17 18 for the one-dimensional case, while in Appendix A.6.2 we extend the approach to a milti-dimensional scenario, considering, for the sake of simplicity, the LQAG problem (but, crucially, including internal noise), to prove the generalizability of Algorithm 2 In Appendix [.662 we provide the solution for Eq. 17] with the same procedure applying Eq.
18 We also discuss the potential extension to the full noise model. Lastly, in Appendix A.7] we
examine the assumption of linear dynamics and extend our approach to switching linear dynamical system to make it less restrictive.
Experiments To empirically validate our iterative algorithm, we apply it the same one- dimensional problem discussed in Section 3.2 The FPOMP algorithm aligns with the optimal solutions found by the GD algorithm, resulting in identical solution and performance (Figs. p..c). In Fig. 2 a, the FPOMP algorithm follows the same cost trend as the GD algorithm. with the curves Appendix AA8.II.
overlapping (the FPOMP curve omitted for clarity, see also Fig. 10
In Appendix A.8.2 we analyze the same colti-dimenssonall task as in Section 3.2, excluding multi- plicative noise but including internal noise. The results show that the FPOMP algorithm matches the GD optimal solutions for the controller and outperforms TOD when o, > 0. Even in its current form, FPOMP surpasses the method from [ when internal noise is considered.
4 Conclusion
In this paper, we provide a novel approach for solving stochastic optimal control problems adapted to the noise characteristics of the human sensorimotor system. Our work builds the seminal

---

study in [ where the classical LQG framework (called here LQAG) is extended to the LQMG framework to include both control and signal-dependent noise, as well as internal noise in the estimation process. This extension provides a more realistic description of the sensorimotor system, enabling the reproduction of larger sample of behavioral data in motor control, albeit at the cost of reduced mathematical tractability.
However, the solution derived in M. which is widely used 227 33
36 37
38
suffers from an ill-conditioned derivation. Specifically, that solution assumes unbiased estimators -a condition that, as we prove numerically and conceptually in this work, does not hold, leading to suboptimal performance when internal noise is considered or before algorithmic convergence. This suboptimality arises from the close relationship between unbiasedness and the orthogonality principle of an optimal estimator, where the former, mathematically, implies the latter. Yet, the orthogonality principle is satisfied by the optimal filter only in the absence of internal noise and under algorithmic convergence.
Assuming only that control is linear in the current state estimate, derive an alternative algorithm that optimizes control and estimation without requiring unbiasedness. The optimal solution is obtained by propagating sufficient statistics compute the expected cost, which is minimized via numerical gradient descent on filter and control gains. For a more constrained, but still relevant, version of the problem, we derive the analytical counterpart, which alternates between forward propagation of moments and backward optimization of control and filter gains until convergence. This makes our approach suitable for high-dimensional problems, significantly reducing computational cost.
We demonstrate superior performance in the presence of internal noise and before convergence is reached (that is, when filter control gains are fixed at suboptimal values, regardless of the level of internal noise), and provide both mathematical and heuristic explanations. Joint modulation of control and filter gains helps filter internal fluctuations, enhancing adaptability and generalization across internal noise levels, as discussed in Appendix A.5.2 By applying our algorithm to a sensorimotor task, make novel behavioral predictions that distinguish our solution from previous ones. Specifically, we find that control gains decrease with increasing internal noise, leading to smoother control signals in individual trials. This results in slower movements with reduced trial-to- trial variability.
In summary, algorithm extends optimal feedhack control broader range of problems in systems neuroscience.
Limitations and Future Work One limitation of our work is the assumption of state- cute-indeppndent filter gains for the optimal estimator: in the presence of multiplicative noise, non-adaptive estimation proves sub-optimal. Additionally, incorporating more realistic cost functions could extend our framework beyond the traditional quadratic dependence. Further investigation into the connections between our optimal control law and biologically plausible learning rules [43] may also be necessary. Morcover, we have not formally demonstrated the convergence properties of our algorithm to a global minimum, although our algorithm is guaranteed to converge at least to a local minima, and did not find any numerical evidence for multiple local minima. The next immediate step is to derive the FPOMP algorithm for the general case with multiplicative noise, as discussed Appendix A.6.2
10

---

Acknowledgments and Disclosure of Funding
This work was supported by: NEI-NIH R01 EY016178; Grant PRE2021-097771 funded by MI- CIU/AEI/IO. 13039/501 10001 and by "ESF+": "Project PID2023-146524NB- financed by MCIN/AEI/1O. (303//50110001 1033/ ERDF, EU, the Spanish State Research Agency and the Euro- pean Union, and ICREA Academia.
References
[1] Emanuel Todorov.
Stochastic optimal control and estimation methods adapted the noise characteristics of the sensorimotor system. Neural computation, 17(5):1084-11 108, 2005.
[2] A Aldo Faisal, Lue PJ Selen, and Daniel M Wolpert. Noise in the nervous system. Nature reviews neuroscience, 944)292-3303 2008.
[3] David W Franklin and Daniel M Wolpert. Computational mechanisms of sensorimotor control. Neuron, 72(3):425 442, 2011.
[4] Emmanuel Guigon, Pierre Baraduc, and Michel Desmurget. Computational motor control: redundancy and invariance. Journal of neurophysiology, 97(1):331-347, 200,
[5] Stephen H Scott. Optimal feedback control and the neural basis of volitional motor control. Nature Reviews Neuroscience, 5(7):532-545, 2004.
[6] Emanuel Todorov and Michael I Jordan. Optimal feedback control as theory of motor coordination. Nature neuroxcience, 1211)): 2002.
[7] Emanuel Todorov. Optimality principles in sensorimotor control. Nature neuroscience, 7(9):907-
915, 2004.
[8] Mark Davis. Stochastic modelling and control. Springer Science & Business Media, 2013.
[9] Robert F Stengel. Optimal control and estimation. Courier Corporation, 1994.
[10] Gerald E Loeb, WS Levine, and Jiping He. Understanding sensorimotor feedback through optimal control. In Cold Spring Harbor symposia on quantitative biology, volume 55, pages 791-803, Cold Spring Harbor Laboratory Press, 1990.
[11] Arthur Kuo. An optimal control model for analyzing human postural balance. transactions on biomedical engineering, 42(1))87-101, 1995.
IEEE
[12] Bruce Richard Hoff. A computational description of the organization of human reaching and prehension. University of Southern California, 1992.
[13] Tamar Flash and Neville Hogan. The coordination of movements: experimentally confirmed mathematical model. Journal of neuroscience, 5(7):168--1703 1985.
[14] Christopher M Harris and Daniel M Wolpert. Signal-dependent noise determines motor planning. Nature, (94(6695)77800-784 1998.
[15] GG Sutton and K Sykes. The variation of hand tremor with force in healthy subjects. The Journal of physiology, 191(3):699-711, 1967.
[16] Emanuel Todorov. Cosine tuning minimizes motor enrors. Neural computation, 14(6):1233- 1260, 2002.
[17] Richard A Schmidt, Howard Zelaznik, Brian Hawkins, James S Frank, and John T Quinn Jr. Motor-output variability: a theory for the accuracy of rapid motor acts. Psychological review, 86(5):415, 1979.
[18] Reza Shadmehr and John W Krakauer. A computational neuroanatomy for motor control. Experimental brain research, 185:359-381, 2008.

---

[19] Joseph Y Nashed, Frédéric Crevecoeur, and Stephen H Scott. Influence of the behavioral goal and environmental obstacles on rapid feedback responses. Journal of neurophysiology. 108(4)999-1100, 2012.
[20] Kelvin E Jones, Antonia F de C Hamilton, and Daniel M Wolpert. Sources of signal-dependent noise during isometric force production. Journal of neurophysiology, 88(3):1533- 1544, 2002.
[21] Chistina A Burbeck and Yen Lee Yap. Two mechanisms for localization? evidence for separation-depeendent and separation-indeeendent processing of position information. Vision research, 30(5):739-1)) 1990.
[22] David Whitaker and Keziah Latham. Disentangling the role of spatial scale, separation and eccentricity weber's law for position. Vision research, 37(5):515-524, 1997.
[23] Emanuel Vassilev Todrov. Studies of goal directed movements. PhD thesis, Massachusetts Institute of Technology, 1998.
[24] Rubén Moreno-Bote, Jeffrey Beck, Ingmar Kanitscheider, Xaq Pitkow, Peter Latham, and Alexandre Pouget. Information-limiting correlations. Nature neumscience, 17(10):1410-1417, 2014.
[25] Philippe Vindras and Paolo Viviani. Frames of reference and control parameters in visuomanual pointing. Journal of Experimental Psychology: Human Perception and Performance, 24(2):569, 1998.
[26] Andrew M Gordon, Goran Westling, Kelly J Cole, and Roland S Johansson. Memory represen- tations underlying motor commands used during manipulation of common and novel objects. Journal of neurophysiology, 69(6):189-1179, 1993.
[27] Matthias Schultheis, Dominik Straub, and Constantin A Rothkopf. Inverse optimal control adapted to the noise characteristics of the human sensorimotor system. Advances in Neural Information Processing Systems, 34:9429-9444 99499,99422244 2021.
[28] Sanjoy K Mitter. Filtering and stochastic control: A historical perspective. IEEE Control Systems Magazine, 16(3):67-76, 1996.
[29] Tryphon T Georgiou and Anders Lindquist. The separation principle in stochastic control, redux. IEEE Transactions on Automatic Control, 58(10))2481-244) 2013.
[30] Peter Joseph and Julius Tou. On linear control theory. Transactions of the American Institute of Electrical Engineers, Part II: Applications and Industry, 80(4): 1961.
[31] Henk Van de Water and Jan Willems. The certainty equivalence property in stochastic control theory. IEEE Transactions on Automatic Control, 26(5):1080-1087, 1981.
[32] Mark M Churchland, Afsheen Afshar, and Krishna V Shenoy. A central source of movement variability. Neuron, 52(6): 1085-166, 2006.
[33] Dominik Straub and Constantin A Rothkopf. Putting perception into action with inverse optimal control for continuous psychophysics. Elife, 11:e76635, 2022.
[34] Jonathon W Sensinger and Strahinja Dosen. A review of sensory feedback in upper-limb prostheses from the perspective of human motor control. Frontiers in neuroscience, 14:345, 2020.
[35] Dan Liu and Emanuel Todorov. Evidence for the flexible sensorimotor strategics predicted by optimal feedback control. Journal of Neuroscience, 27(35):9354 9368, 2007.
[36] Jun Izawa, Tushar Rane, Opher Donchin, and Reza Shadmehr. Motor adaptation as a process of reoptimization. Journal of Neuroscience, 28(1 1):2833-289,, 2008.
[37] Tomohiko Takei, Stephen G Lomber, Douglas J Cook, and Stephen H Scott. Transient deac- tivation of dorsal premotor cortex or parietal area impairs feedback control of the limb in macaques. Current Biology, 3177)1447-1487, 2021.
12

---

[38] Maryam M Shanechi, Ziv M Williams, Gregory W Wornell, Rollin C Hu, Marissa Powers, and Emery N Brown. A real-time brain-maching interface combining motor target and trajectory intent using an optimal feedback control design. PloS one, 8(4))e59049, 2013.
[39] Jur Van Den Berg, Pieter Abbeel, and Ken Goldberg. Lgg-mp: Optimized path planning for robots with motion uncertainty and imperfect state information. 2011.
[40] David Winter. Biomechanics and motor control of human movement. John wiley & sons, 2009.
[41] Philipp Karg, Simon Stoll, Simon Rothfuß, and Sören Hohmann. Inverse stochastic optimal control for linear-quadratic gaussian and linear-quadratic sensorimotor control models. In 2022 IEEE 6lst Conference on Decision and Control (CDC), pages 2801-2808. IEEE, 2022.
[42] Joseph Y Nashed, Frédéric Crevecoeur, and Stephen Scott. Rapid online selection between multiple motor plans. Journal of Neuroscience, 1495)1769-17788)) 2014.
[43] Johannes Friedrich, Siavash Golkar, Shiva Farashahi, Alexander Genkin, Anirvan Sengupta, and Dmitri Chklovskii. Neural optimal feedback control with local learning rules. Advances in Neural Information Processing Systems,
16370, 2021.
[44] Kaare Brandt Petersen, Michael Syskind Pedersen, al. The matrix cookbook. Technical University of Denmark, 7(15):510, 2008.
[45] Philip Becker-Ehmek, Jan Peters, and Patrick Van Der Smagt. Switching linear dynamics for variational bayes filtering. In International conference on machine learning, pages 553-562.
PMLR, 2019.
[46] Steven L Brunton, Bingni W Brunton, Joshua L Proctor, and J Nathan Kutz. Koopman invariant subspaces and finite linear representations of nonlinear dynamical systems for control. PloS one, 112771 2016.
13

---

A Appendix
A.1 Background: Classic LQAG Solutions
The optimal control and filter gains, L1.... T-1 and K1,.. 2- T22, for the classic LQAG problem
defined by Eqs. 2 can be derived analytically and are given by [8
L. - + I BT S
(21)
S. = Q: + A S (A-.))
(22)
K. = ( + Pw)
(23)
1+ S2 (A K,H)E,AT,
(24)
A.2 Background: Solutions from
The algorithm proposed in []] alternates between optimizing control and estimation. As for the classic LQAG problem, the control is optimized iteratively in a backward-in-time fashion, while keeping the filters K, fixed. The solution derived by using the method dynamic programming, writing down the Bellman equation for the optimal cost-to-go assuming the unbiasedness of the estimator [I] From there, the optimal filters K₁ are found at fixed L,, again by minimizing the cost-to-go. Taken together, these two optimization steps lead to iterative algorithm that supposed to provide the optimal solution to the control problem [ For completeness, we present here the optimal solutions for L, and K, from [1]
The optimal control gains are given by the following backward algorithm
L = (R, ST B CT((S + C1+ $
(25)
S Q AT Sz (A ALB
Σ D.K.S + K,D:
(26)
S A Str BL, + (A- (-KKH S4 (A K,H)
(27)
with Sr = Qr and St = Note that in | the optimal control law is defined as u, L,
whereas use Eq. to compare the solutions from 1 with ours, we need to invert the sign of the control gains L,
The optimal filter gains follow instead a forward optimization
=; = AE + Ω + D. + Et + Dfe + Ee D!!
(28)
E1 = 22g + G. + (A KHHEE + CCLEE
(29)
DE = Pr + K.HΣ AT + (A - BL, E (A BL.)
(30)
+ (A BL, HTK. + K,H (A BL.)
Fiee
-++ = (A BI,) Ee (A K2H)T
(31)
D = (
(32)
with 5 = En,, Σ = क and E = 0. Note that in 1 Σg = E[eee]], where e == 3t St Σ = E [ and E = [
A.3 Unbiasedness and Orthogonality Principle: How Internal Noise Affects Optimality
In Section 2.3 we discussed the invalidity of the unbiasedness condition. Here, we provide the details of the dinediimal problem used to numerically validate this assertion, presenting the plots in Figs. A.3b. also for the case where o, 0.3. The system parameters are listed in Table 2 in
14

---

Appendix A.5.1
where on VQa. The trial duration is set to T = 10 time steps, and we vary on across the values 0.0, 0.3, 0.6. At & = 8, we compute [ [ as function of Ft. To do this, we collect the values of za and ēg over 102 trials. bin the data for & with bin size of d& 0.1, and compute the mean of 2 within each bin. The standard deviation is shown as error bars. Note that the choice of is arbitrary.
(e)
= 0.0
(0)
= 0.0
ER 2:] -
[[sla) -
-10
(4)
0.0
03 Elel
[[ F.l
A
Figure 4: The invalidity of the umbiasedness condition. Here we plot B[2; [e), for given value tf as a function of &, for different og, using the solutions from M. The conditional expectation [[xa[(]]] is computed through Monte Carlo (MC) simulations. (a) [[xx[[ee] as function of Gt for On = 0 (dots with error bars given by the std of our MC estimate). The gray dotted line stands for the bisector, where Ef(a [+x] = fz. (b) Same (a), but for on = 0.3. (c) Same as (a) and (b), but for o, = 0.6. (d-f) Absolute value of the distance between [2 [&e] and 2g as a function of [ for
Tn = 0.0, 0.3, 0.6. The gray dotted lines represent E[ [ār] Lz-
A.3.1
Orthogonality Principle and Suboptimality
we outlined the relationship between the unbiasedness condition and the orthogonality principle, highlighting how this connection results in the suboptimality of the approach in | when internal noise is present. For simplicity, and without loss of generality, we consider ID scenario, m p where the orthogonality principle implies
In Section 2.3
Ω E(aa(c, 2,)] E [
E Eira 0
(33)
This results in the estimation error being orthogonal to the estimate itself 2, We also set d 1. We demonstrate here that the condition [R|r [ 4., used
[I] to derive the optimal control law, implies the orthogonality principle and that this principle cannot be satisfied by optimal estimator when internal noise is present.
Assuming
Bir [2x] ft.
(34)
and multiplying by ft on both sides and then taking the expectation over dt we obtain
E E [[r
(35)
corresponding Eq. 33 In the absence of internal noise, the optimal filter gains K₁ can be found by imposing the orthogonality principle, without the need to minimize the cost function. It holds
Q-1 = (K²H² -K,
+ (A2 (AR+K? ABL, - 2AKH - BLK,
(36)
+ (AK,H- - -KH"HH)Z + KKΩ + Pry
If we use Ω = 0, as in due to the initial conditions, we can solve the equation Qt Vt I,.,T, T, obtaining equation for K,,
15

---

K₁ - - D + +2w)
(37)
2(H2T, D2B[=? + S2w)
with
= B[a 01)
(38)
For 22 = 0. Eq. 37 simplifies to
AHT, + AHT;
K2 2(H D2 + S2w)
(39)
Observing that the solution K; would correspond to open-loop strategy. sub-optimal (sensory information would not be integrated) for a stochastic partially observable system as the one we are considering, we get for the optimal filter gains
AHT:
K₁ HJ, D +2w
(40)
It can be shown that the solution in 1 for Q.j 0 aligns with Eq. 40 We observe that Eq. 40 can replace Eq. 20 in Algorithm to optimize the filter gains. For 2n 0, this leads to the optimal solution. Thus, in the absence of internal noise, the optimization of control and estimation can be performed using two separate objective functions: one enforcing the orthogonality principle for the optimal estimator, and the other minimizing the cost function for the optimal controller, regardless of the multiplicative nature of the noise. This could also be relevant for more biologically plausible scenarios
However, when S2, > O, the existence real solution for Eq. 37 depends on the initial conditions and is no longer guaranteed. Moreover, as we demonstrate in Appendix [A5.2) the optimal solutions do not satisfy Q2 = for S2,, > 0. Therefore, the orthogonality principle holds for an optimal Kalman filter only when 22 O. Consequently, the algorithm in [1] assumes unbiased estimation, which should imply the orthogonality principle, even in cases where it no longer applies to the optimal estimator.
Separation Principle, Orthogonality Principle, Unbiasedness: A Brief Digression Unbiased- ness, orthogonality, and the separation principle are related but distinct concepts. Here, we briefly clarify their differences and commonalities.
The separation principle stems from the formulation of the classic LQAG problem, where the optimal solutions for control and estimation are mathematically independent, allowing for their separate optimization. However, with multiplicative noise, this independence is lost M. We have shown that this breakdown occurs even with additive internal noise and zero multiplicative noises.
The orthogonality principle (in ID) states that [[xx B[ meaning that estimation error and estimate are orthogonal. This condition holds for an optimal Kalman filter only in the absence of internal noise (see Fig. 5h.. Internal fluctuations, however, disrupt the mathematical independence between control and estimation, invalidating the orthogonality principle as well. These two concepts distinct: for instance, with internal noise but non-zero multiplicative noise, the orthogonality principle would still hold, yet the mathematical independence between control and estimation would be broken.
The unbiasedness condition (which, as previously discussed, never holds) states that E [z] = at implying the orthogonality principle. This explains the optimality of the solutions in [1] in the absence of internal not due to the validity of the unbiasedness condition, but because the orthogonality condition holds.
A.4 Novel Algorithm for Optimal Control Problems
A.4.1 Derivation of Closed-Foo -for Equations for Moments Propagation
We explicitly derive Eqs. 14 here. Notably, no approximations are required to propagate the first two moments of the joint variable (ar, . in closed form, as both control and estimation are linear in the state and state estimate (see Eqs.
g). Consequently, Eqs. 15 hold regardless of the
16

---

distribution of (a, R). By taking the expected value of Eqs. 7 state estimate and sensory feedback, we obtain
over the joint distribution of state,
+++1] = AB + BLE
(41)
[++1] a + (A BL, - KcH)B[ad)
(42)
which correspond to Eq. 14
Similarly, we compute the second non-central moments of the joint
variable (z, f), resulting in
B 1[1--1 = [ 1L
(43)
+ BL CL +
+4 (A BL - K
+ K/ (A+ ( (44)
+ K,DE| [KDDE| +0,
1( K. ( -
(45)
+ K, + (A+ BL, K,
B[ #
El 12
(46)
From this, we can derive Eq.
15
Since the cost function is quadratic in the state and state estimate, the variables μ4 and E,, defined in Eqs. 10 serve sufficient statistics compute E[J] (Eq. 9). which is all that is needed derive the optimal control and filter gains.
A.4.2 Pseudo-Code
we minimize the expected accumulated cost TE[]], computed through Algorithm [ with respect the filter and control gains Li,... the function "Gradien(Descent(" in the "Optim.jl il" Julia package. The hyper-parameter of the used algorithms are listed in Table in Appendix A.5
For the GD algorithm (Section 3.1
-- Li,...T-1,
and K,,... T -2 using
Algorithm 1 Propagation of the expected cost GD algorithm
1: Input: μ4 Σ1 (initial conditions of the system), TT1, Ki.". K T-2, and the system parame- ters (A, B. H, Ci-I,.., Di=1....d Dgg Swi 2.y).
2: Output: [[J]
3: Algorithm steps:
4: B[J] 0
5: Pold μ1
6: Eoid Σ
7: for each iteration t = 1,2.... 2........ T do
8: E[J] + E[J] E[j.]. (Eq. 9
9: Update M, and G. (Eqs. 12 3
Enem M,Eada + G.
11:
Anew = Mepoad
12: Eold f Enew
13:
Pold f- Pnew
14:
end for
A.5 Experiments: GD Algorithm
The hyper-parameters, of all the used algorithms are provided in Table [ For the GD algorithm (Section 3.1) we minimize the expected accumulated cost [F]]., computed through Algorithm using the function "GradientDescent"" in the "Optim.jl" Julia package.
17

---

Table 1: Hyper-parameten of the used algorithms
Algorithm
Description
value
GD
Number of iterations of the "GradientDeescent(" function
100000
FPOMP
Number of iterations of the control-estimation optimization
1000
TOD
Number of iterations of the control-estimation optimization
1000
A.5.1
One-Dimensional Case: Parameters
We set d = 1.
Table 2: Parameters of the one-dimensional problem
Name
Description
value
Linear map for the system dynamics
1.0
B
Scaling of the control signal
1.0
Scaling matrix for control-dependent noise
0.5
D
Scaling for signal-dependent noise in the sensory feedback
0.5
H
Observation matrix
R
Control-dependen cost at each < T Task-related cost at each time < Task-related cost time | = T
1
Q
Qr
20
time steps
100
E[a]] B[m]]
Initial condition for the mean state and state estimate Initial covariance of the state
1.0
Em
0.0
Es.
Initial covariance the state estimate Covariance matrix of the additive Gaussian noise &
0.0
0.52
Covariance matrix of the additive Gaussian noise We
0.52
On
Standard deviation of the additive internal Gaussian noise M/
{0.0 0.1 2.0)
A.5.2 One-Dimensional Case: Understanding the Qualitative Differences
When the Orthogonality Principle Is No Longer Optimal As discussed in Section 2 and Appendix(A.33. the presence of internal noise causes the optimal estimator to no longer satisfy the orthogonality principle. Here, demonstrate this result numerically, using the same one-dimensional problem presented in Section 3.2
0.9
2.0
6:
0.6
1.5
0.3
1.0
0.0
0.0
0.5
1.0
1.5
2.0
0.0
0.5
1.0
1.5
2.0
0.0
0.5
1.0
1.5
2.0
Internal noisn G7
Intrrnal moksc
O,
Figure 5: Filtering out the internal fluctuations. (a) P2 averaged over time (we indicate the time average with (-)), as function of a, for TOD and GD algorithms. (b) (T) as a function of og, (c) (e) function of on. The error bars (mean 1SEM from Monte Carlo simulations) are not visible as too small.
The optimal solutions do not minimize Ω - [ E [
(Fig. 5). Instead, the optimal strategy uz (Fig. 5p). This allows the system to filter out internal fluctuations affecting the estimation process, reducing their correlation with the latent
appears favor lower values T, E[=]] E
state dynamics, æ. As result, the absolute estimation error, leel = [ J vhl + Ir,
slightly (but significantly) larger for the GD solutions (Fig. 5b), which seems help decorrelate internal noise from the state evolution.
18

---

This "decorrelation mechanism' is achieved through an intertwined modulation of control and filter gains. In the next paragraph, we provide geometric interpretation of this behavior through an eigenvector decomposition of the dynamical system under investigation.
Eigenvector Decomposition and Adaptability of the Solutions In one dimension we can write the update equations for
r3 E[
E [r
(47)
f2 E E[ - E[r
(48)
Te+1 B(+1
2g + C
Mr
+ K2 K2 DZ
(49)
where
(-
A
-BL,
M, = (A-K,H)
(-KHH K,H
(50)
A BL, K.H
The eigenvectors of M, are given by
w = (7) 1 (BLijKiH)
(51)
(52)
Note that the angles 8, between these two eigenvectors are the same as the angles between the eigenvectors of the matrix M,. Indeed, the eigenvectors of M, are given by
0
(53)
0/2 (-BLi/KHH)
(54)
A parity operation (along the araiss) maps ones into the others, preserving the angles.
The optimal solution arises from the adjustment of the angle 0 between the two eigenvectors (in this onc-dimensional case). As 7, increases, 0 also increases, due to the joint modulation of L, and K₁ with o, (Fig. 64). This increase in 0 allows the system filter internal fluctuations more effectively and better generalize to other levels of o, (Fig. (D).
optima
(a)
TOD
CD
(6)
TOD
40
1500
30
Test
1250
1000
20
750
500
100
25
50
75
100
optim
optim
Figure 6: Intertwined modulation of control and filter gains to deal with the internal noise. (a) Angles 0, between the two cigenvectors the matrix M, (Eq, [nd and see next Paragraph), at different levels of internal noise opptim for TOD and GD algorithms. (b) "Adaptability" of the two solutions: the solution found by the GD algorithm (right panel) generalizes better than the one by TOD (left panel) when optimized for certain level of internal noise, o opptim and tested another one, o test, for larger optim the generalization property improves thanks to due modulation of 9-.
By examining the modulation Th as o, changes the (T, - -Q4) plane, we can offer heuristic interpretation of the different solutions found by the TOD and GD algorithms. As on increases,
19

---

the angle between wj and ti grows for both algorithms. However, this modulation is much more pronounced in the GD solution (Fig. (a).
Furthermore, if only additive noise were considered, there would be no modulation 0g with on in the TOD solution (for confirmation, Appendix A.8.2 without multiplicative noise, TOD's derivation does not modulate the control gains with (.))).
The joint modulation of Le and Ke causes W2 move closer to the y-axis in the GD solution (Fig. 7 green line). This configuration results in more effective filtering of internal fluctuations, decoupling them from the latent state dynamics, since these fluctuations occur on S2 (see Eq. 44).. This result aligns with the observed decrease L, as on increases (see Figs. and 3). where lowering the control gain moves dh closer to the y-axis. Thus, this eigenvector analysis qualitatively explains the trends observed in Fig. for (T) as function of in both the TOD and GD solutions.
2
= [pp]] Baah Bpp] 2.
Figure 7: Figenvector decomposition the dynamics. We show here qualitative representation of the cigenvectors of the matrix M, in the plane (2). The black arrow represents the "shared" eigenvector wi while the blue (green) arrow represents W2 for TOD (GD) solution. Note that the optimal L, are negative, while the optimal K₂ are positive (Fig.
A.5.3 One-Dimensional Case: Improving Performance Without Internal Noise
We briefly show that, even in the absence of internal noise, if the algorithm has not yet converged, the solution proposed by [1] does yield the optimal control law. We demonstrate this in one- dimensional example, using the same parameters shown in Appendix A.3. (but the result is valid in general), while only varying the scaling matrix for the multiplicative sensory noise and keeping ", = O. We fix the filter gains at the suboptimal constant value K1,... T = = 1.0, and optimize the vector L L T-1 using TOD and GD algorithms.
TOD
3
GD
2
10
20
0.2
0.4
0.6
0.8
1.0
Figure 8: Enhanced performance when optimizing control at fixed filter gains and zero internal noise. We plot the expected accumulated cost E[J], computed by averaging the quantity from Eq. over 50% trials, as a function of the scaling matrix D. with error bars (mean + 1SEM from Monte Carlo simulations, error bars not visible too small), for the two algorithms TOD and GD.
We find that TOD control law leads to higher expected accumulated cost E[J] (Fig. The improved performance the GD (and FPOMP in Section 3.2 we show that the solutions of the two match) algorithm arises from not assuming unbiasedness when optimizing control. As a result, the
20

---

algorithm adjusts the control gains to account for the bias introduced by the suboptimal estimator. Similar to the case with internal noise, the control gains found by the GD algorithms in this scenario are typically smaller than those found by the algorithm of []] A similar behavior occurs when optimizing the filter gains while keeping the control gains fixed and suboptimal. The difference in performance is due to the fact that the derivation of the optimal estimator in [ is only valid when the
the control gains are optimal. Thus, the algorithm in [I] does not apply when the suboptimality controller needs to be "ballncedd by the estimator. When the algorithm in [ is run in its entirety, optimizing both control and filter gains iteratively, these issues are resolved. Upon convergence, the controller and estimator optimally adjusted to each other, in the absence of internal noise.
A.5.4 Multi-Dimensional Case: Parameters
For the sensorimotor task described Section 3.2 the discrete-time dynamics is the same as in
[
p(l + At) p(() + p())
(55)
p(t At) -p(t) + (tt)2t/mm
(56)
ftt + At) =(()(1 - A1/12) + g(1)
(57)
gt At) +(t At/r) u(t) dtE))t/TT
(58)
We have therefore the following system parameters (with 1)
At
0
0
1
At/m
0
A
0
0
1
At/T2
At/12
(59)
0
0
0
1
A//m.
B
=
(0
0
0
At/T)
(60)
C=
(0
0
0.
(61)
0
0
H
(62)
0
0
0
100
0000
0
D
(63)
1 0
O
O 0
0
Q 111
(64)
0
0
0
Qr = pp
+ OUT
+ fff
(65)
T
R1 T-1
(66)
T
-11
Rr = 0
(67)
T
0
0) 0)
(68)
T = (0
W
0
(69)
f = (0
0
Wh
(70)
0
0
0
0
53
(71)
0
0
0
0
(72)
0
2n =
120
(73)
0
21

---

with the initial conditions given by
B[ri] (z
0)T
(74)
B[s1]
B
E
B((11]
(75)
0
0
Σz1
(76)
0
0
0
0/
0
0
0
Σe.
(77)
0
0
0
The parameters of the problem are listed Table B (std standard deviation).
Table 3: Parameters the sensorimotor task
Name
Description
value
Al
time-step (8)
0.010
m
mass the hand (Kg), modelled as a point mass time constant of the second order low pass filter time constant of the second order low pass filter Control-dependent cost at each < T
T
0.04
0.04
r
1e-5
W
Task-related cost at time t T for the velocity Task-related cost at time
0.2
wy
for the force
0.01
T
time steps
100
2
Target position
0.15
Target position standard deviation std of the additive Gaussian noise & std of the additive Gaussian noise w/ std of the control-dependent noise &
0.0
dg
0.0
0w
0.0
Oe
0.5
0p
std of the signal-dependent noise p
0.5
On
std of the additive internal noise f/e for the position estimate
{0.0, 0.005, 0.5}
On.
std of the additive internal noise 7/ acting on the velocity estimate
Ous
std of the additive internal noise V/e for the force estimate
One
std of the additive internal noise mi for the estimate of g
Note that the initial condition for the state 2 the actual target position: in such a way the control signal aims at minimizing the distance from I, 0.
A.5.5 Multi-Dimensional Case: Impact of Internal Noise
We provide here an extended analysis of the impact of internal noise in the sensorimotor task discussed in Section 3.2. offering additional insights. We compute the posterior variance of (state belief variance) to assess how internal noise affects For on 0, 0.005, 0.05, 0.5 = (an arbitrarily chosen time), we find [²[z]] 2 10
10 8. 10 10-5,2 10 respectively. In the case of O, 0.05, where internal noise accounts for about 10%, there is significant increase in state uncertainty compared to the scenario without internal noise. Indeed,
80, since both the state z; and the state estimate â;j are near zero (with the target position as the reference point in our coordinate system), internal noise becomes the dominant source of fluctuations, as we are only accounting for multiplicative noise in this task.
A.5.6 Scaling to Higher- Dimensional Problems: An Application
We demonstrate how our algorithm scales to high-dimensional problems, building on discussion in the final paragraph of Section 3.2) We implement high-dimensional task to show the generalizability of the GD algorithm. The same results would apply analytical counterpart, the FPOMP
algorithm, discussed in Section 3.3 and Appendices A.8.1 A.8.2 In this scenario, we set the
22

---

dimensions of the state, control, and observation to m 10, = 4, and & = 10, respectively. Note that this significantly increases the dimensionality compared the problem in Section 3.2 (for the multi-dimensional case).
(e)
(
120
TOD
TOD
0.08
110
100
0.06
80
0.02
60
0.00
0.25
0.75 1.00
0.25
0.50
1.00
internal noise B
internal noise By
Figure 9: High-dimensional lask. (a) Expected accumulated cost a function of og for TOD (blue dots) and GD (green dots) algorithms. We see that even in this nigh-dimensional task, GD solutions outperform the ones from (.] To compute the expected cost, we used Algorithm (but the results are confirmed by Monte Carlo simulations). (b) Pscudo-determinan of the control gains L (averaged over time), denoted [[pp as a function oy for TOD (blue dots) and GD (green dots) algorithms.
The system matrices A, B, and are random matrices with elements drawn from a standard normal distribution (mean zero, standard deviation onc), while C' is defined C = ag og The matrix H is the identity matrix, and the time horizon is set to T = 10. All elements the state and state estimate vectors are initialized one. We used og a a o. 0.5 and varied o, across values of 0.0, 0.5, and 1.0. The matrices defining the quadratic cost functions, Q and R, are identity matrices at each time step. All the findings from Section 3.2 are confirmed in this high-dimensional setting (Fig. 9). The GD algorithm continues to outperform the solutions in [1, with performance improving as internal noise increases, and the control gain magnitude decreases internal fluctuations grow. In fact, as internal noise increases, the optimal strategy involves reducing control over the system. To quantify control magnitude, we compute the pseudo-determinant of I,,.. T and average it over time. The pseudo-determinat,, a generalization of the determinant for non-square matrices, provides measure of the volume scaling induced by the control gains.
A.6 FPOMP Algorithm: An Analytical Counterpart the Numerical GD Algorithm
A.6.1 One-Dimensional Case
In the one-dimensional case we have m k = 1. Additionally, to simplify the notation, we set d = We start by defining
F
A (B² + C2
F.
F12
(78)
Fr,
2ABL,
/r,1
KH( + D2)
G,
G1
(A+ BL,)² K² -2AK, 2BL,KH 2BL + 2AK, 22?
(79)
G1,3/
Hr1
AK,HI ABL, BLL BL,KH A - ABL, - -K,H + BL,KHH
f.
H
(80)
Hr,
In one dimension, we can then propagate the non-central moments as
E = Bl, + Fr,
+ F + 22g
(81)
B| = B,, 2 Gr, + K?S + S2
(82)
B| = H., E[
+ H
(83)
23

---

The derivatives of the non-central moments with respect to Lg and K₁ obey the following equations, for
l.... T,
0B[r[+d
at+i-1,124
11ti-1,1
(84)
0h:
@R[12td
Utti-1,244
beti-1.2
(85)
0h.
0141-1344 + 1145-133
(86)
0h.
and
0B[a]+1 a141- 8,41-111
(87)
0K,
0K, a141- + 8141-1.2
(88)
OE [xt+sS++i 0141-133K, + 3145-1.3 kk,
(89)
with ā. U, ā and 3 given by the following recursive equations
F141 σe
01+1,1
d++1 =
01+1,2
G ae
(90)
(1++1,3/
H1+1 02
(b+1,11)
Fi41 6g
b141 =
1.+1.2
=
G by
(91)
b,+1.3/
Hi+1 be
F1+1
[04+1,1
01+1,2
G Cr
(92)
04,1,33/
H44 C2.
(2+1,1
F 3
Be41 =
1++1,
=
G B,
(93)
B,41.3/
A44 B.
The initial conditions for Eqs.
90 93 are
2(B²
C2 EE+))
a
2B
(94)
2
2ABE( DABELT, AB
6,
2B,, (B
(95)
+ E[zx0.1) - BK, E(rr
0
6
2H² + 2F(x++)) + 20 + 22D2
(96)
0
3, 2 B (L
E x+ - E BL, 2 - E(2,2,))
(97)
-
By observing that the expected accumulated cost, Eq. 9 (adapted to the one-dimensional case), will
be a function of E E| and B for t = 1,., T t, and by using Eqs. 84 89
rewrite Eqs.
24

---

17 18888 as
0
E 14, = 2R,
a 1m0
T-t
(98)
E [Quuietti-11 + tet i=1 + (Qitibe4i-111 +
=0
and
0
T-t Eii Mo, Σd]
T-1
Σ
IiQu+a44--111 + 414
0K;
(99)
==1
+ (Q4418e + R
Therefore, from Eqs. 98 99
have the following instantiations of Eqs. 19 20 for the optimal
control and filter gains at time 1, I and K;
L Luum
L;
(100)
L.den
Knum
K2 Kden
(101)
with
Lt Luum
(Qussbei-111 RetiL Detii1)))
(102)
2R, I den
T-t
(103)
Q414++1111 + R4
and
T-1
K
(Q44B44s-111 + (e
(104)
Kden
(Q4+004i-111 Rut.
(105)
We can then use Eqs. 100 101 the one-dimensional problem.
to implement Algorithm and extract Li... T- T-1 and K K - - for
A.6.2 Multi-Dimensional Case
For the i 1,.., and D, for = 1..., d) in the presence of internal noise (S,, 2 0).
case, we derive Eqs. 19 20 for the classic LQAG problem (C; = 0 for
As a title of example, we derive here Eq. 19 for the optimal L (to be used in Algorithm 2 but the approach would be the same for the optimal filter gains K The extension to the more general scenario including the multiplicative sources of noise would follow similar method. As outlined in the expected cost per step is given by
Section 3.1 Eq. 9
E[j++1
E[m++1] [ (14, + E[ J R B(oudt
(106)
+ FrrQQus + Tr[ E
for 0, .,... T
When computing Bijj-a114 Ed] to write down Eq. 17 (with C, == I,..., and D; i 1,... d), and derive Eq. 19 the coefficients mulupplyng E [[,,] coming from the term
25

---

B[ BattsB in Eq. 106 will be the same as the ones multiplying Σs, and coming from the The same holds for the coefficients multiplying respectively E[r,]B[a,]" and
term E.
Similarly, we group together the coefficients coming from the other two factors B[ R2 Bit+ and 14, in Eq. 106
We now note that the terms dependent on L, appearing Eljr- [4,, E,] will show a dependence on the afore-mentioned moments [a [B[ Ea, [B(r [ [[ and Dr..at- More specifically, the quadratic factors in La will only depend on [[ and Ed.- Taken together, these observations lead to the following form for Eq.
17
Te J 19 0 + S P.E O,
(107)
where we have used Dg - [[ = E[ and Er,... + 11 = E|,,,]]
Therefore, find the optimal control gains L; from Eq. 107 only need to compute the coefficients Tr. S₁ and Pt, similar to what we have done for the one-dimensional case in Appendix [A.6.I] As before, we compute the coefficients Te. S, and P; by only looking at the first two terms appearing
[ in Eq. 106 that is B[x TQ++ B[ae+c] and B[ Rg T E [ }. By using ((44)
Oo X
OUT
(108)
0x
Oह XT
(109)
0X
(TT X 2N U0T
(110)
ax
where ü and w are vectors and N is symmetric matrix, we obtain
T-1
[ (Q++ + RetiBee Ve+tr1] J. = 2R, +2
Σ
(111)
i=1
T-t
S = 2 Σ
tti (1 + Rotebet (1 Vit [Rits
thi
(112)
1=1
T-t
P. 2
Σ
Viti-1 Qetti uti (u
LI LI4, Retelete (ui
(113)
with Ve+i given by
V44+ = II (A BL
(114)
j=1
for i = I,...,? T - l, and
V B.
(115)
In Eqs.
112 113
(2 )
c+i
is a vector whose elements are m m matrices:
(+t ( (( (utti (12000 2
(116)
The subscript (:-)) indicates the initial condition 0) for the evolution of 4+i (L,=0.1..) with T denoting the m m identity matrix and 0 being an m x m matrix whose elements all zeros, e.g..
(117)
(L, is updated through the following equations
Miph Miph
for i
11)00 =
s+i
MM+t-11 L.
for i = 2, 2,.., T
(118)
26

---

with Mr given by Eq. 12 and M, having the same form as the block matrix M,, but with L, 0,
M,
(119)
K,H
- A-KH
From Eq. 107 we can then write for Eq. 19
= - (S, +
P.B(O, (AECEE)))) B[G
(120)
where .t denotes the pseudoinverse operation. Note that J. is a symmetric x matrix with del [e. > and therefore invertible. Due to the initial conditions for Es, and E[ the symmetric matrix [[ could have a null determinant: for this reason we use the pseudoinverse operation. This consideration relevant only for an initial transient: after a certain time 0 (depending on the dynamics parameters) would have det TE [ > 0 and [[ E [[ due the properties of the pseudoinverse. With Eq. 120 we can implement Algorithm 2to find the optimal why, mathematically, the control gains decrease when the internal noise level is increased: the factor [[ will get bigger and bigger as S2 gets larger.
control gains IT- Notably, from the form of Eq. 120
The derivation of Eq. 20 for the optimal filter gains K? would follow the same procedure. To extend the presented approach to the case with multiplicative noise, we need to propagate the terms depending on C. 1, I,... and D;, i l.....d coming from the factors (rr(Q++2mm, and similarly what we have done with the other terms in Eq. 118
Tr[L,, Ret:Le+
Eq. 106
but using (144)
a
Tr Tr[AX AT,,
(121)
0x
Tr(AXTB]= BA
(122)
ax
0
Tr = ax
(123)
where Å, B and Ò are matrices. We observe that, even when considering multiplicative noise, Eq. 120 will still be valid: only the matrices Jr, St and Pe will change, including now also the terms depending on C, I,... c, and D,, I..... d.
A.6.3 Pseudo-Code
Algorithm 2 (Section 3.3) extracts the optimal solutions analytically by identifying the critical points of the total expected cost, conditioned the first two moments of the joint variable distribution (t). In the pseudo-code, (k) and k²(*) stand for, respectively, the control and filter gains at time t and at optimization step k. The hyper-parameters of the used algorithms are listed in Table in
Appendix A.5
Algorithm 2 FPOMP algorithm
T-2 (initial guesses for L and K*), and the system parameters (A, B, H, Cl=i....e> D D==1....ds Sgg Zw P,y).
Input: μl, Σn (initial conditions of the system), (1)
1,,- T-1 K (1)
i..., -
2: Output: Li.... 111 Ki...,.- -222 (optimal control and filter gains)
Algorithm steps:
4: for each iteration & 2, ..... optimization steps do
MT,--T-11, E1,, E1..,. <- Eqs. 14 15
using 1,, -1.--T--11 and K (<
(k-1)
1,.-,7-2
6: for each iteration 1,....
do
t- -
L f(m, E., L (k) T- K (3-1) 1+1, 2
8:
1+1.,..
K!a) glue, Ex, 111 1-
K (k)
1+1,-, x++1...,T-2) 2
10: end for
end for
Li,...T T-1 I (k)
12:
"1.....--
2-2 (K)) 1,
Ki..... T1..
27

---

A.7 Extension to Switching Linear Dynamics
We discuss here how to extend our approach to switching linear dynamics. One of the underlying assumptions in this work and in [ is that the agent has complete knowledge of the updating rules of the latent dynamical system. By using the same set of matrices to update the state and the state estimate, we implicitly assume that all uncertainty in the estimation process arises solely from noise sources: the problem of inferring the matrices A and B goes beyond the objectives of this approach. For this reason, to extend our work to the more general and realistic case of Switching Linear Dynamics (SLD), we can consider matrix A depending on the time step t, A, A complete formulation of SLD might require adding another variable, a discrete switch variable 8t regulating the way the matrices A, vary with time and context [45] Given that in our case the agent has access to the updating rules of the dynamical system, we can omit 8, (the agent does not have infer 81 and A₂) and direetly consider the case in which we have predetermined set of matrices A T---- The same applies to the matrix B, that can be replaced by B1... TT-1 Note that to preserve linearity we assume At and B, to be independent
and â. We consider here the multidimensional case to be as general as possible. To extend the GD algorithm we only need to modify the block matrix M, that we use to update the moments Σx, l4 and eventually propagate the expected cost E[J] through Eq. [3] Indeed, once we can compute the expected cost at fixed control and filter gains, Lı.... .1-1, use Algorithm to define the objective function to be minimized through gradient descent with respect to L, and K,. To update the block matrix M2 we have to substitute A and B respectively with A, and B, in Eq. T2 To handle the potentially high computational costs of performing a numerical gradient descent, we introduced the analytical counterpart of the GD algorithm, the FPOMP algorithm. For the one-dimensional case, it supports all the noise sources mentioned in Section 2 (additive, multiplicative and internal). We extended this algorithm to the multi-dimensional case for additive and internal noise for the sake of simplicity, leaving the more general version for future work (A (Appendix AA..2 outlines how this can be done). Here, we extend the afore-ment mentioned FPOMP algorithm (for both mie-dimnnioonl and multi multi-dimensional cases) to switching linear dynamics, following similar procedure to that of the numerical algorithm. For the mne-dimensional case, we replace and B respectively with A, and B, Eqs. [8 and 94 97 For the multi-dimensional case we have to substitute A with Ac+j and B with Besj in Eq. 114 and B with B, in Eq. 115 Finally, previously done, replace A with A, in Eq. 119 for M, With these changes, we can implement Algorithm for the case with switching linear dynamics.
and K,....TT T-
We observe that the extension to switching linear dynamics aims to make the assumption of linearity less restrictive. Additionally, given the flexibility our approach handling high-dimensional systems (Section 3.2), it is reasonable to think that this assumption does not limit the effective description of lower-dimensional nonlinear dynamics, potentially by employing the Koopman operator [46].
A.8
Experiments: FPOMP Algorithm
A.8.1
One-Dimensional Case
We show that the FPOMP and GD algorithms yield the same performance the one-dimensional problem introduced in Section 3.2
confirming their equivalence as discussed in Section 3.2
E|JGD GD J FPOMP!
-10
0.0
0.5
1.5
2.0
internal noise dy
Figure 10: Accumulated cost difference. Difference of [BJ] for GD and FPOMP solutions (computed by averaging the quantity from Eq. over 50k trials), function of og, with error bars (SEM).
28

---

A.8.2 Multi-Dimensional Case
We consider the same problem described in Section 3.2 (see also Appendix A.5.4 but without multiplicative noise sources, to validate the FPOMP algorithm derived in Appendix A.6.2 The parameters of the problem are listed in Table
Table 4: Parameters the problem sensorimotor task without multiplicative noise
Name
Description
value
At
time-step (s)
0.010
mass of the hand (Kg). modelled point mass time constant of the second order low pass filter time constant of the second order low pass filter Control-dependee cost at each / < T Ruskerelated cost time T for the velocity Task-related cost at time = T for the force time steps
1
T1
0.04
72
0.04
le
Wy
0.2
w
0.01
T
50
Target position 0=
2
0.15
Target position standard deviation std of the additive Gaussian noise & of the additive Gaussian noise wir std of the control-dependent noise
0.0
gg
0.5
0w
0.5
0.0
0p
the signal-dependent noise Pr std of the additive internal noise т) for the position estimate
0.0
On
(0.0, 1.0, 2.0}
One
std the additive internal noise T acting on the velocity estimate
Onf
std of the additive internal noise 7, for the force estimate
0
dng
std the additive internal noise Y for the estimate
0
0
(n)
P0D
FPOME
600
TOD FPOMP with edlmmter
10
200
1000
6
10
20
30
40
0
10
20
30
40
internal noise Og
(c)
(d)
FPOMP
10 102
PPOMP with CD catmater
GD
10
1.30
-400
-600
0.0
0.5
1.0
1.5
2.0
0
10
20
30
40
50
10
20
30
40
50
internal noise O
Figure 11: Application of the FPOMP algorithm. (a) R[J], computed by averaging the quantity from Eq. 3 over 50k trials, as a function of a, with error bars (mean + 1SEM from Monte Carlo simulations, error bars not visible too small) and for TOD, GD and FPOMP (at fixed filters given by TOD solution).
(b) Magnitude of the control gain vector as function of time for TOD and FPOMP (at fixed filters given by TOD). (c) Same as (a), but comparing GD and FPOMP (now at fixed filters given by GD). (d) First component of the vector L, for the solution given by GD and FPOMP (now fixed filters given by GD solution).
29

---

We recall that, for readability, we derived solutions only for the optimal controller, the procedure for the optimal estimator follows similarly but was explicitly derived. Extensions the optimal estimator and the general case with multiplicative noise are left for future work, as discussed in Appendix A.6.2
Optimizing the control gains L, using the FPOMP algorithm (with fixed filter gains K₁ from the TOD solution) leads improved performance (Fig. Map. orange dashed line) internal noise increases. However, this solution is not fully optimal, as the estimator is still optimized using the TOD algorithm. When L, and K; are both optimized with the GD method, lower accumulated cost is achieved (green dashed line). An interesting feature of our algorithm is that, being fully analytical, it can enhance numerical solutions. Due to a potentially shallow parameter landscape (vamishing gradient) or limited computation time, the GD optimization may stop near the global optimum without fully reaching We find that re-optimizing the control gains using the FPOMP algorithm after the GI) solution for filter gains yields a small but significant performance boost (Fig. TTc), with minimal changes in the final Lt vector (Fig. (Td). This also confirms that our algorithm finds the optimal solutions. Extensions to estimator optimization and the multiplicative noise case discussed in Appendix A.6.2
The qualitative trends observed in the sensorimotor task (Section 3.2 Fig. are confirmed: control magnitude decreases as internal noise increases (Fig. ITp). Additionally, while the TOD solution does not modulate control with respect to internal noise when only additive noise is present, the FPOMP algorithm introduces such modulation, leading to a lower accumulated cost (Figs. l.b.)). consistent with our discussion in Appendix A.5.2
30

---

NeurIPS Paper Checklist
1. Claims
Question:
Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer:[Yes
Justification: The claims made in the abstract and introduction are accurately confirmed by the results provided in the results and commented in the conclusion.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
The limitations of the work discussed in the conclusion.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
31

---

Justification:
The assumptions and mathematical proofs are explained in the paper and presented in more details in the Appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: All the information is provided in Section 3 and the Appendix, including pseudo-codes and detailed equations for the implementation of the proposed algorithms.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that appropriate to the research performed.
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
32

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification: The codes to generate the data and the figures are provided in the supplemental material.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nnis.cc
aublic//uudees//odeetubiiss caapssionoolyyyy for more details.
While encourage the release of code and data, understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
/ wwtp:/cccppppgguuuiiaeeeae
Stteeee//ooSSSeeissSioica www.ssssiicgcyy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: The details about the implementation of the algorithms and the parameters of the experiments are discussed in the Appendix, in Section A.5
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented the core the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: We present error bars representing the standard error of the mean calculated across multiple trials.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
33

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [No]
Justification: The derived algorithms be efficiently performed standard commercial CPU.
Guidelines:
The answer NA means that the paper does include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental rums as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https::/neurips.. wttpariipp.cpppupc/eeeeeeaaa ctw.eeeeeeeicciiceeeeieenee
Answer: [Yes]
Justification: We present a theoretical study aimed at improving current algorithms for solving optimal control problems. Our research is entirely computational and does not involve human or animal subjects, thus raising no ethical concerns.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: Since study is theoretical, it does not have any clear societal impacts.
Guidelines:
34

---

The answer NA means that there is no societal impact of the work performed.
If the authors answer NA No, they should explain why their work has societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (c.g., gated release models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: The theory presented in this paper carries no potential for misuse.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [NA]
The paper relies solely on original codes and data, without utilizing any existing assets.
Justification:
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
35

---

For scraped data from a particular source (e.g., website). the copyright and terms of service of that souree should be provided.
assets are released, the license, copyright information, and terms use in the package should be provided. For popular datasets, paperssi thood cww datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: The paper does not introduce new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not make any use of human data.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethics, workers involved data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The research presented in this paper does not include any studies involving human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
36

---

Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
37


Reviews:
[{"review_id": "ILnmQyHScI", "reviewer": "Reviewer_5RsR", "summary": "The work introduces a new algorithm that adapts to the noise typical in human sensorimotor systems. It includes different types of noise like those from control movements to improve how accurately we can predict motor behaviors. This joint optimization of control and estimation significantly outperforms traditional models, especially in noisy environments, offering precise predictions for sensorimotor tasks.", "questions": "-\tThe paper does not clearly explain how the proposed gradient descent-based algorithm differs from or improves upon the existing algorithm in Todorov (2005). More discussion is needed on the specific pitfalls discovered in Todorov's derivation and how the proposed approach addresses them\n-\tBackground section lacks detail. Some improvement suggestions: having goal of LQG problem upfront, explaining key assumptions of the LQG problem, better details/introduction to Todorov’s approach. In general, more thorough summary of key concepts & results would make it accessible to reader less familiar with this topic\n-\tMore explanation needed for “realistic noise model”\no\tmultiplicative noise terms $C(u_t)$ and \\sum(\\rho_i D_i x_t_ are introduced, however ranging of c and d are not defined\no\t\\nu_t is referred as \"internal noise” and details regarding how it cannot be directly filtered is discussed. This area is not clear. What does “internal noise” mean here conceptually? Why is it important to include in the model? \no\tIt is unclear why $\\omega_t, \\rho_t$ are modeled separately. What is the practical significance here?\n-\tThe paper claims, η_t in the state estimator dynamics (Eq. 2) leads to challenges for Todorov's approach. Why is that the case? Better explanation will strengthen this work\n-\tThe experiment is done via reaching task (Section 4.2) that uses smaller number of states, control inputs and feedback signals. How do the results generalize to problems with higher-dimensional controls and observations?", "limitations": "-\tThe assumption of linear dynamics (Eq. 1) and linear observations (Eq. 3) is quite restrictive. \n-\tThe assumptions about the nature and impact of noise in the system might not fully capture the complexity of noise in biological systems", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "-\tThe algorithm provides a more realistic noise models. This resonates with unpredictability seen in real-world sensorimotor tasks\n-\tRobust performance even when internal noise is significant.\n-\tEmpirical validation through simulations to demonstrate its effectiveness and cost", "weaknesses": "-\tDoes not provide any theoretical guarantees (e.g. on convergence) for the proposed GD algorithms.\n-\tModel predominantly utilizes quadratic cost functions which might not capture all the nuances.\n-\tIntroduction of more realistic noise models increase mathematical and computational complexity of the algorithm", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their comments. We will respond to each comment point by point.\n- Theoretical guarantees\n\nOur numerical algorithm, being a GD, converges to, at least, a local minima. In [7] convergence to at least local minima is also found. \nConvergence of our algorithm is numerically shown in Fig. 8, 9, 10, as we find that our algorithms match the solutions of [7] for no internal noise. \n- Quadratic cost functions\n\nWe refer to Todorov’s work and the latest applications in inverse optimal control theory [27].  Our approach could be extended to non-quadratic cost functions by approximating state and state estimate distribution using a multivariate Gaussian [27]. For smooth scalar cost functions we can apply Stein’s lemma (Mamis, arXiv, 2022) to compute the expected cost, and generalize Algorithm 1. This solution would significantly increase the mathematical complexity. We will comment on this in the new version of the manuscript.\n\n- more realistic noise models increase complexity \n\nIncluding a realistic noise model is crucial for predicting behavior [7]. While our analytical algorithm becomes more complex in this scenario, it does not prevent us from finding the optimal solution. The algorithm focuses on minimizing the expected cost over control and filter gains, making it easy to interpret. When considering only additive and internal noise, the complexity of our algorithm is similar to that of [7]. Our approach uses a standard minimization procedure to identify critical points, comparable to the derivation of the Bellman Equation in [7]. We provide all necessary codes as a reference.\n\n- how the proposed algorithm differs from or improves upon Todorov’s \n\nOur algorithm outperforms [7] in the presence of internal noise (see Figs. 2a, 4a). \nThe enhanced performance arises from correcting a condition used in the derivation in [7]. Specifically, in [7], when deriving the optimal controller by minimizing the optimal cost-to-go, the estimator is assumed to be unbiased.\nWe prove that this condition fails, especially with internal noise, leading to suboptimal solutions (see Fig 1 in the uploaded pdf and Sections 3.3, A.2 and A.5.2). Our algorithm addresses the issue by deriving the optimal solutions without relying on the incorrect unbiasedness condition. We find a closed-form solution for the expected cost, that can be minimized to find the optimal strategy, without further assumptions, apart from the linearity of the controller. We apologize if  this was unclear. We will further explain the above in the new version of the manuscript.\n\n- Background section \n\nWe will surely include the suggestion in an improved version of the paper. Thanks for noticing this. \n\n- More explanation for “realistic noise model” \n\nThe realistic noise model was not detailed extensively as its formulation is based on [7] and is not novel in our paper. This model is also used in inverse optimal control problems to explain behavioral data [27]. In the next version, we will describe the model and define all terms clearly.\nFor the scaling matrices $C_i$ and $D_i$, the ranges of $c$ and $d$ are not explicitly defined because the general model only implies that the relationships for motor and sensory noise are multiplicative [7,27]. \nIncluding internal noise in the estimation process is crucial, as discussed in [7]. Neurophysiologically neural representations are subject to internal fluctuations [1,2,24,28]. The state estimate, affected by internal noise should model, at the algorithmic level, the computation performed by the agent to infer the latent state of the system.  Behaviorally, adding internal noise makes the model more accurate in explaining data (Section 6 in [7]). \nWe note that the gains $K_t$ of the Kalman filter cannot act directly on $\\eta_t$ (see Eq. 2).\n$\\omega_t$ and $\\rho_t$ allow for both additive and multiplicative noise in the sensory feedback, aiming for a realistic noise model [7, 27].\n\n- $\\eta_t$ challenges Todorov's approach\n \nThe presence of $\\eta_t$ makes the use of the condition of unbiasedness particularly inappropriate to derive the optimal solutions, leading to the suboptimality shown in Figs. 2a and 4a. Indeed in this scenario, the orthogonality principle, which is implied by the unbiasedness, does not hold anymore for the optimal filter (see Sections 3.3, A.2 and A.5.3). \nNoticing, and solving, this pitfall is a central contribution of our paper.\n- generalization to higher-dimensional problems\n\nTo provide a fair comparison with [7], we took the example with no change, apart from including internal noise. This is considered to be a quite realistic model for reaching. Our algorithm does not have any limitation in terms of dimensionality of state, control and observations. To demonstrate this, we have included a higher-dimensional task (see Fig. 2 in the uploaded pdf). \n- assumption of linear dynamics \n\nWe acknowledge that the assumption might seem restrictive. However, we refer to the linear quadratic framework, building on the pivotal work of Todorov [7]. This formulation has proven to be realistic enough to describe behavioral data, as discussed in [7, 27]. Additionally, we have extended our work to linear switching dynamics, as discussed in the general author's rebuttal: this makes the assumption of linear dynamics less restrictive. Lastly, we note that, given the flexibility of our approach to handle high dimensional systems, we might think that this assumption does not preclude the effective description of lower-dimensional nonlinear dynamics (using Koopman operator, see Brunton et al, Plos one, 2016).\n- the noise model and the complexity in biological systems \n\nWe agree with this point. However, to the best of our knowledge, this noise model is the most realistic one that still allows for mathematical tractability, leading to a fully interpretable model. We will further discuss all the above limitations in the new version of our manuscript. Thanks for pointing them out.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "qyOpiD2v6z", "reviewer": "Reviewer_HNsh", "summary": "The paper extends a popular stochastic optimal control framework for explaining the algorithmic function of sensorimotor circuits in the brain to the case where there is multiplicative noise in feedback and motor output, which is a more realistic assumption for noise in the brain. Prior work assumed the estimator was unbiased in the presence of internal noise and prior to convergence. The authors show that this assumption is wrong and propose a means to circumvent it, assuming only linearity of the feeback law. They give two gradient-descent based procedures for optimizing the parameters of the model. In particular, the Iterative-GD algorithm develops a more efficient alternating optimization of the free parameters, improving its scalability to larger problems. The authors evaluate their method on two toy linear dynamical system problems with either internal but additive noise or internal and multiplicative noise, and show that it outperforms the existing prior method.", "questions": "- How well does the developed algorithm scale to larger problems?\n- How much better does it fit to actual experimental data compared to the prior method?", "limitations": "The authors clearly discuss the limitations of their approach in the conclusion, such as state-independent filter gains and the choice of cost function.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "- Modeling sensorimotor circuits in the brain is an important and open problem in systems neuroscience. And the proposed model makes more realistic assumptions about the noise characteristics in the brain.\n- The proposed solution for optimizing the parameters of the model is fairly simple, and the iterative algorithm is especially scalable to larger problems.\n- The example which illustrates how the original formulation of the estimation procedure cannot be unbiased due to internal noise is clear and intuitive.\n- For the toy problems considered, there are clear advantages of the proposed method over prior work.\n- The paper is well organized and overall written well. It provides a thorough related work section and does a good job explaining the novelty and results.", "weaknesses": "- The two systems considered in evaluations are very simple. Although they capture the main issue with the prior method well, it would strengthen the paper to consider larger problems with more degrees of freedom and complex dynamics. For instance, the choice of linear dynamics may not be suitable. A switching linear dynamical system may be a more reasonable assumption here.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their comments. We will respond to each comment point by point.\n\n- Switching linear dynamics \n\nThanks a lot for the interesting suggestion. We discuss here how to extend our approach to switching linear dynamics. \nOne of the underlying assumptions in this work and in [7], as outlined in Section 2, is that the agent has complete knowledge of the updating rules of the latent dynamical system. By using the same set of matrices to update the state and the state estimate, we implicitly assume that all uncertainty in the estimation process arises solely from noise sources: the problem of inferring the matrices A and B goes beyond the objectives of this approach.\nFor this reason, to extend our work to the more general and realistic case of switching linear dynamics (SLD), we can consider a matrix A depending on the time step t, $A_t$. A complete formulation of SLD might require adding another variable, a discrete switch variable $s_t$ regulating the way the matrices $A_t$ vary with time and context (see Becker-Ehmck et al, PMLR, 2019). Given that in our case the agent has access to the updating rules of the dynamical system, we can omit $s_t$ (the agent does not have to infer $s_t$ and $A_t$) and directly consider the case in which we have a predetermined set of matrices $A_t, t=1,...,T-1$. The same applies to the matrix B, that can be replaced by $B_t, t=1,...,T-1$. Note that to preserve linearity we assume $A_t$ and $B_t$ to be independent on $x$ and $\\hat{x}$.\nWe consider here the multidimensional case to be as general as possible. \nTo extend the GD algorithm we only need to modify the block matrix $M_t$ that we use to update the moments $ \\Sigma_t, \\mathbb{S}_t, \\mu_t $ and eventually propagate the expected cost $\\mathbb{E}[J]$ through Eq. 6. Indeed, once we can compute the expected cost at fixed control and filter gains, $L_t, t = 1,...,T-1$, and $K_t, t = 1,...,T-2$, we can use Algorithm 1 to define the objective function to be minimized through Gradient Descent with respect to $L_t$ and $K_t$. To update the block matrix $M_t$ we have to substitute A and B respectively with $A_t$ and $B_t, t=1,...,T$ in Eq. 9.\nTo handle the potentially high computational costs of performing a numerical Gradient Descent, we introduced the analytical counterpart of the GD algorithm, the Iterative-GD. For the one-dimensional case, Section A.4.1, it supports all the noise sources mentioned in Section 2 (additive, multiplicative and internal). We extended this algorithm to the multi-dimensional case for additive and internal noise in Section A.4.2 for the sake of simplicity, leaving the more general version for future work (Section A.4.2 outlines how this can be done). Here, we extend the afore-mentioned Iterative-GD algorithm (for both one-dimensional and multi-dimensional cases) to switching linear dynamics, following a similar procedure to that of the numerical algorithm. \nFor the one-dimensional case, we replace A and B respectively with $A_t$ and $B_t$ in Eqs. 34,35,36,50,51,52 and 53.\n\nFor the multi-dimensional case we have to substitute A with $A_{t+j}$ and B with $B_{t+j}$ in Eq. 70 and B with $B_t$ in Eq. 71.\nFinally, as previously done, we replace A with $A_t$ in Eq. 75 for $\\tilde{M}_t$.\nWith these changes, we can implement Algorithm 2 for the case with switching linear dynamics.\n\nWe will add the optimal controller and estimator for the switching linear dynamics as described above in the improved version of our manuscript. \n\nWe also considered a task with a higher number of degrees of freedom (see Fig.2 in the uploaded pdf).\n\nFinally, we emphasized that the sensorimotor task borrowed from [7] and discussed in Section 4.2, can already provide interesting behavioral predictions that match with realistic scenarios, as also outlined in [7].\n\n- larger problems\n\nThanks again for the comment. The proposed numerical algorithm, referred to as GD in the paper, supports high-dimensional problems without the need for further adjustment. However, as the dimensionality of the problem increases, the computational cost rises significantly. Specifically, this algorithm performs a numerical Gradient Descent over mp(T−1)+mk(T−2) parameters, where m, p, and k are the dimensions of the state, control, and sensory feedback, respectively, and T is the number of time steps. To demonstrate that the algorithm works for higher-dimensional problems, we apply it to a higher-dimensional task, as discussed in the general author rebuttal.\nThe analytical counterpart of the numerical algorithm, called Iterative-GD in the paper, scales effectively to higher-dimensional systems with the same flexibility and computational costs as the approach proposed by Todorov in [7]. \nIt is important to note that this algorithm, which exactly matches the solutions of the numerical one in the described scenario, already outperforms Todorov’s when internal noise is considered, as shown in Fig. 10.\nAdditionally, we emphasize that the problem discussed in Section 4.2, despite its simplicity, already provides a good approximation of a realistic scenario, with the capability of predicting recorded behavioral data (see Section 6 of [7]). \n\n- fit to actual experimental data \n\nWe have not yet used our new approach to model behavioral data, but this is our next planned step. \nWe anticipate that having internal noise to model behavioral data will be very important, as already pointed out by Todorov (2005) (his Figure 1C) who argues that internal noise is needed to model an arm reaching task. Therefore, having a model with the optimal controller with internal noise seems to be a crucial step to not only better fit the data, but also have an interpretable model with the correct rationality assumptions.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "lxyjaLR6rE", "reviewer": "Reviewer_C3X4", "summary": "The authors found that the algorithm developed by Todorov (TOD) for optimal feedback control problems, in the presence of internal and sensory noise, assumes that the state estimation is unbiased at all times. However, they show that this assumption does not hold in the presence of internal noise and, even when internal noise is absent, before convergence of the algorithm. This causes the solution proposed by Todorov to be suboptimal when realistic noise is factored in the system. With this motivation, the paper presents a gradient-descent based algorithm (GD), which works without assuming the state estimation to be unbiased. In systems with internal noise, GD outperforms TOD in terms of cost reduction, producing qualitatively different behavior.", "questions": "* The title seems excessively long, is it possible to make it more concise?\n* I am not familiar with the use of > and >= to indicate positive definite and positive semidefinite matrices, used throughout the paper. Is this a common notation?\n* The authors mention that analytical solutions to the LQG problem can be obtained assuming independence between control and estimation optimizations (separation principle). However, they also show that the separation principle follows from the assumption that the state estimation is unbiased, and not the opposite. They also say that the unbiasedness assumption is what was used in Todorov's paper. It is therefore not clear that the separation principle implies unbiasedness of the state estimate, and therefore the analytical tractability of the LQG problem. Could the authors elaborate on this point?\n* Why are the results of iterative-GD excluded from the main text?\n* What is the difference in computation time between GD, iterative-GD and TOD in the presented experiments?\n* Which experiment supports the statement \"when optimized for high levels of internal noise, the system generalizes well to other noise levels\"?", "limitations": "Although short, the limitations section adequately addresses them.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "The paper is clearly written and describes a well-executed study, where the identified problem is analyzed theoretically, numerically and qualitatively. The motivation is convincing, as the authors find an important problem with a very relevant algorithm used to study behavior, popular in neuroscience research and particularly in motor control. The authors explain the problem in detail, with a theoretical analysis of why the assumptions of TOD are not satisfied in typical scenarios, and also qualitatively showing why it is the case.\nThe proposed gradient descent algorithm seems a valid alternative to solve the optimal feedback control problem avoiding the unbiasedness assumption. The \"iterative-GD\" version of the algorithm seems also a valid alternative to avoid gradient descent updates, replacing them with an alternated optimization scheme.\nThe numerical results are promising, as they show that GD has strong performance compared to TOD.", "weaknesses": "* While the authors identify two reasons why the unbiasedness of the state estimation does not hold, the paper does not discuss the relative importance of the two. From fig. 4, it seems that TOD and GD perform equally well when internal noise is absent, suggesting that the unbiasedness of the state estimation before convergence does not cause TOD to underperform.\n* The moment propagation technique assumes a fixed time horizon, which is a hyperparameter of the algorithm. However, the paper does not discuss how the solution depends on this parameter choice. In general, the choice of the hyperparameters is not discussed. Furthermore, it is not clear how the computational cost of the estimation of E[J] scales with the horizon T.\n* I find the name \"iterative-GD\" improper for the second algorithm, as it rather looks like a fixed point iteration combined with an \"expectation-maximization-style\" update.\n* Due to the large number of results and experiments, a considerable amount of information is provided as supplementary material. This paper would have benefitted from a longer submission format.\n\nMinor:\n* Line 51, \"the estimator will be unbiased\". While the meaning will become clear later in the paper, at this point this statement needs to be better specified.\n* Section 2, the matrices C_i and D_i are only defined in the following subsection\n* Line 112, space before a comma\n* Line 148, \"imposing only the linearity of u\". While I understand what the authors mean, they should specify that they mean that u is a linear function of the state estimate\n* Line 280, \"leads to good a generalizability\"", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their comments. We will respond to each comment point by point.\n\n- unbiasedness\n\nWe first note our new result, not present in the paper: unbiasedness does not hold even in the absence of internal noise  (see Fig 1 of the uploaded pdf). Therefore, for the unbiasedness to break down, one does not need internal noise. The idea is simple: in Fig. 1 of the manuscript, we intuitively explain that a large fluctuation in internal noise can lead to a large $\\hat{x}_t$, while $x_t$ does not change. This causes the the unbiasedness condition  $\\mathbb{E}[x_t|\\hat{x}_t] = \\hat{x}_t$ to break. But the point is that the same would happen, without internal noise, if there is a large fluctuation of sensory noise: in this case again, $\\hat{x}_t$ would be large without changing $x_t$. \nThe reason why Todorov finds the correct solution with zero internal noise is not because the unbiasedness condition is correct (which is wrong in all cases, as described above), but because with zero internal noise the unbiasedness conditions implies orthogonality, and this condition holds true (see Sec. A.2). The orthogonality principle is the one that is satisfied at the optimal solution – but only after convergence of the algorithm, so orthogonality does not hold if the filter is suboptimal. \n\nWe will explain this better in the new version of the paper. Thanks a lot for the question. \n\n- fixed time horizon\n\nThanks a lot for the question. The time horizon T follows [7] for fair comparison. Iterative-GD has similar computational costs to Todorov’s. However, numerical GD significantly increases costs with T, as it optimizes over mp(T−1)+mk(T−2) parameters, where m, p, and k are the dimensions of the state, control, and sensory feedback.\n\n- name \"iterative-GD\" \n\nWe thank the reviewer for the suggestion and propose “fixed point optimization with moments propagation” (FPMP) algorithm. \n\n- Longer submission format\n\nThanks for noting the density of results. In the next version we will try to organize the whole thing better\n\n- Long title\n\nWe agree with the reviewer and propose “Improving optimal control and estimation for realistic noise models of the sensorimotor system”. \n\n- positive matrices\n\nWe followed the notation used in [7], which aligns with standard notation in the field of linear algebra, but please let us know if you think a different notation would be more suited for our paper\n\n- Orthogonality and separation principle, unbiasedness\n\nWe are sorry for being unclear. Unbiasedness, orthogonality and separation principle are related but distinct concepts. \n\nThe separation principle comes from the formulation of the classic LQG problem: the optimal solutions for control and estimator are mathematically independent, allowing for separate optimization of the two. With multiplicative noise, this independence breaks down [7]. We argue that this happens already with additive internal noise, contrary to what said in [7]: if we look at Eq. 4.2 in [7], without multiplicative noise, $L_t$ and $K_t$ do not depend on each others even for non zero internal noise. \n\nThe orthogonality principle states that $\\mathbb{E}[x_t\\hat{x}_t] = \\mathbb{E}[\\hat{x}_t^2]$. This condition holds for an optimal Kalman filter, without internal noise (Section A.2 and Fig. 3a). Then, internal fluctuations break the mathematical independence between control and estimation, invalidating at the same time the orthogonality principle. \n\nHowever, these two concepts are separated: e.g., if we have no internal noise, but non zero multiplicative noise, the orthogonality principle would still hold, but the mathematical independence between control and estimation would break.\n\nIn Section A.2, we show that without internal noise, the orthogonality principle is satisfied when $K_t$ follows Eq. 33. Empirical validation of the optimality of this comes from finding the same $K_t$ with GD, TOD, and Eq. 33 (without internal noise). See Section A.2, see also Fig. 3a.\n\nThe condition of unbiasedness (which is always wrong, as described in the first question) implies that $\\mathbb{E}[x_t|\\hat{x}_t] = \\hat{x}_t$. In Section A.2 we show that if $\\mathbb{E}[x_t|\\hat{x}_t] = \\hat{x}_t$ holds, then also the orthogonality principle does. This explains Todorov’s optimality in the absence of internal noise: not because the unbiasedness condition is true, but because the orthogonality condition holds.\n\nWe intend to add a full section in the Appendix discussing all the above in the new version of our manuscript, and also clarify the central messages in the paper.\n\n- iterative-GD results \n\nWe included those results in the appendix section due to the submission format constraints (Iterative-GD matches the numerical GD solutions). However, we agree that including the Iterative-GD results could be a valuable addition.\n\n- computation time\n\nFor the 1D problem, the three algorithms have comparable and short computation times. For the multi-dimensional problem the numerical GD has a significantly higher computation time compared to the TOD algorithm. This is due to the GD optimization operating over a high-dimensional parameter space. While the TOD algorithm runs in a few minutes on a standard laptop, the GD optimization requires several hours (approximately 4 hours). \nTo address this issue, we derived in the previous paper the Iterative-GD algorithm. Its computation time is comparable to that of TOD, as they rely on relatively similar procedures.\n- adaptability to other noise levels\n\nThanks for noticing the unclarity. Figure 6b in Section A.5.3 supports this statement. We tested the optimal solutions with different internal noise levels, and GD solutions generalize well. This aligns with the eigenvector analysis in Section A.5.3. Todorov’s solutions lack this adaptability (Fig. 6b), likely due to limited modulation of eigenvector angles with noise levels (Fig. 6a).\n\nWe will clarify this point in the revision.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "laGEaHQErP", "reviewer": "Reviewer_EneD", "summary": "This paper targets the stochastic optimal control problem with signal-dependent and internal noise by Todorov (2005). The authors question an assumption made in the original paper and propose an approximate solution to solving the problem without this assumption. The authors empirically show that in the case of high internal noise, the fixed approach can lead to quite different trajectories.", "questions": "How important is the requirement that the internal noise needs to be much greater than one for the assumption to be violated? Based on my intuition, the assumption should be violated even for small noises.\n\nTo me, the assumption of unbiased state estimates seems to be also wrong if no noise is present. As we receive noisy observations, the belief gets also biased by the observations and one could do a similar reasoning as with internal noise: Assume that the received noisy observation in a 1D example is larger than the true state. Then, the prior distribution of the state belief will be distorted upwards and its mean will lie above the true state. This seems to contradict the investigation in the paper where the authors claim that without internal noise, there does not seem to be a problem. Any thoughts on this?\n\nAlgorithm 1 L. 10: How do you update $M_t$ and $G_t$? This seems to be the most important line in the algorithm and but it is not clear what happens here.", "limitations": "Yes", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "There has been some recent work on stochastic optimal control using the considered method proposed by Todorov, so the problem formulation, even if dated back to 2005, could be still relevant these days. Spotting and fixing the inappropriate assumption in the original work is interesting and could potentially benefit other researchers in the field. The paper itself is well-written and well understandable (with minor remarks, see Weaknesses). Even if only two simple simulated examples are considered, the evaluation fulfills its job to provide the insight that, in case of high internal noises, trajectories produced by Todorov's approach and the approach presented in the paper can yield different trajectories.", "weaknesses": "The invalidity of the assumption of unbiased estimation is illustrated by an artificial example (section 3.3.). While this provides an intuition that the assumption is not reasonable, assuming a large positive fluctuation for this example seems to contradict the Gaussian assumption of the noise, limiting the significance of the example. I wonder how important it is that the noise needs to be large (specifically much larger than 1), as the assumption should be violated even for small noises (or in relation to the state representation). I also wonder why internal noise is needed to violate the assumption (see questions for this). I believe the paper would further benefit from a real simulated example that illustrates the problem by computing the expectation values and comparing them to the true state values, thereby indicating the violation of the assumptions. I could imagine that the following point also could be helpful to avoid confusion:\n\nThe variables in the paper are overloaded and confusing to me, making the Math parts hard to follow. It is not clear which variables are random variables or vectors. In section 2, all variables, in particular $x$, $\\hat{x}$, seem to be vectors. In section 3, eq. 7, for example, the authors consider expectations of these values. So to me, it is still difficult to grasp, what $\\mathbb{E}[x | \\hat{x}]$ means (e.g. distribution over of the true state or the belief). I think as the paper has as objective to give a proper treatment of the original formulation by Todorov, the notations should be clear to avoid further issues.\n\nFor me, the structure of section 3 is unintuitive. First, the authors present their method (subsec. 1 and 2), afterwards they explain why the assumption in the original work is not valid (subsec. 3), without even using their presented method. It would have facilitated my understanding of the paper if the intuition was first provided and then the solution was presented.\n\nAs the naive gradient decent approach of sec. 3.1 is computationally heavy, the authors propose a \"gradient decent-based iterative algorithm, derived analytically\" which is claimed as novel. To me it seems not to be based on gradient descent, as in the derivation, the gradients are analytically set to zero and closed-form equations are derived assuming the controller and estimator are alternatingly updated (please correct me if I understood it incorrectly). This approach (and the fundamental equations) seem very close to the original approach by Todorov and while the authors claim that they \"do not use dynamic programming\", it strongly seems to me they do. I think there is nothing wrong with still following the original approach by Todorov but I would have hoped to see an analysis of the equations, which parts differ when the covariances of states and beliefs are taken into account.\n\nThe main proposed algorithm (Iterative-GD) does not seem to be part of the main paper's evaluation. The evaluation plots in Fig. 2-4 could also show a third graph with this method. The term \"TODO\" does not seem to be defined, and I assume it refers to the method by Todorov (2005).\n\nThe approach of formulating block matrices for the state and belief and defining a respective dynamical system is very close to the approach in the work of Schultheis et al. (2021), which seems to stem from the work of Van den Berg et al. (2011) [1] which should be referenced accordingly. They also apply a similar approximation of a Gaussian distribution at every time step, which is not explicitly mentioned.\n\nIn the evaluation of the reaching movement, the internal noise at every time step is quite high (std 0.05, while the state lies in the interval [0, 0.15]). I therefore doubt how \"realistic\" (terminology of the paper) this task with such a high fluctuation of beliefs is. As it is used for illustration purposes and the values are frankly stated this is acceptable but the authors might reconsider their wording. Further, it would be advantageous to indicate the state range of the one-dimensional problem in the evaluation more clearly, or even better, instead of stating the internal noise in absolute form, provide relative quantities to enable a more intuitive understanding of the relevant settings where the different formulation is impactful.\n\nLastly, while fixing the inappropriate assumption is interesting and could help future research, the scope of the contribution is on point and rather limited. For my taste, the paper is a bit overclaimed as the novel approach still shows high similarities to the original one which is not thoroughly investigated. If this work will not be accepted to the conference, the contribution for a future submission could be further increased by additional work, such as: application of the method in the inverse optimal control framework as indicated in the future work section, additional experiments (simulated or real-world),  and/or a more thorough theoretical analysis.\n\n__Minor:__\n\nLine 497 seems to be incomplete.\n\n---\n[1] Van Den Berg, J., Abbeel, P., & Goldberg, K. (2011). LQG-MP: Optimized path planning for robots with motion uncertainty and imperfect state information.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their comments. We will respond to each comment point by point.\n\n- unbiasedness and first question\n\nWe first note our new result, not present in the paper: unbiasedness does not hold even in the absence of internal noise  (see Fig 1 of the uploaded pdf). Therefore, for the unbiasedness to break down, one does not need internal noise. The idea is simple: in Fig. 1 of the manuscript, we intuitively explain that a large fluctuation in internal noise can lead to a large $\\hat{x}_t$, while $x_t$ does not change. This causes the the unbiasedness condition  $\\mathbb{E}[x_t|\\hat{x}_t] = \\hat{x}_t$ to break. But the point is that the same would happen, without internal noise, if there is a large fluctuation of sensory noise: in this case again, $\\hat{x}_t$ would be large without changing $x_t$.\nAs suggested by the reviewer, we numerically computed $\\mathbb{E}[x_t|\\hat{x}_t]$ as a function of $\\hat{x}_t$ (see Fig1 of the uploaded pdf) to test whether the unbiasedness condition holds or not. We find that, without internal noise, the violation of unbiasedness is present, but it is only clearly visible for large $\\hat{x}_t$, because close to $\\hat{x}_t = 0$ the bias will be small (the state estimate experiences small fluctuations if we constrain to small values of $\\hat{x}_t$ – this is also why we chose in our example of Section 3.3 a large fluctuation, but the same bias, albeit small, happens with small sensory or internal noise perturbations). We also find that, when considering internal noise, the bias sensibly increases because the internal fluctuations are not filtered at all by the gains $K_t$.\nTodorov finds the correct solution with zero internal noise because, in this case, the unbiasedness condition implies orthogonality, which holds true (see Sec. A.2). The orthogonality principle is satisfied at the optimal solution but only after the algorithm converges; it does not hold if the filter is suboptimal.\n\nWe will explain this better in the new version of the paper. \n\nThanks a lot for your several comments regarding this crucial point. \n\n- notation (“The variables in the paper…”)\n\nWe apologize for the confusion and we will surely improve our clarity on the revised version. \n\nIn Section 2, $x_t$ $\\hat{x}_t$ $y_t$ are vectors of generic dimensions m,k and p respectively, following [7]. \nThe noise terms $\\xi_t$,$\\epsilon_t$, $\\eta_t$, $\\omega_t$ and $\\rho_t$ are discussed in “A realistic noise model of the sensorimotor system”, Section 2.\n\nAs commented in Section 3, the variable $\\mu_t$ is a vector of vectors. We use block matrices and vectors to provide a more compact description of the problem, similarly to [27] and Van den Berg et al. (2011).\n\nFor unconditional expectations $\\mathbb{E}[ ]$, we mean the expectation over the joint distribution of state, state estimate and sensory feedback, while $\\mathbb{E}[x|\\hat{x}]$ is an expectation taken over the distribution $p(x_t|\\hat{x}_t)$.\n\n- structure of section 3 \n\nWe agree that changing the order of the subsections might improve the readability and understanding of the paper, so we will do it in the revision.\n\n- connections with Todorov (“As the naive gradient…”)\n\nWe agree that the terminology for the Iterative-GD might be misleading, and we propose a new name for that, “fixed point optimization with moments propagation” (FPMP) algorithm. \n\nWe acknowledge that our approach is related to dynamic programming. In the cited sentence, that can be removed, we meant that we do not explicitly write the Bellman equation for the optimal cost-to-go, as in [7]. We compute the expected cost at fixed filter and control gains and derive it with respect to $L_t$ and $K_t$ to find the critical points. \n\nAssuming linear control and propagating the expected cost in closed form allowed us to avoid using unbiasedness during derivation. In contrast, [7] Section 4 requires an expectation conditioned on state, state estimate, and sensory feedback. Unbiasedness is introduced to ensure control depends only on the state estimate. However, unbiasedness is incorrect, leading to suboptimality.\n \n- Iterative-GD in main paper and typo \n\nThe term “TODO” is a typo for “TOD”, as stated by the reviewer.  We apologize for that.\nWe agree that including the Iterative-GD results in the paper would be a valuable addition to the paper. \n\n- block matrices\n\nWe thank the reviewer for the useful suggestions. We will cite the work of Van den Berg et al. (2011) and clarify the connection in our paper. \n\nOur algorithms do not assume Gaussianity. Moment propagation is exact, with the quadratic cost function and linearity of control making the first two moments of state and state estimate sufficient statistics. Deriving optimal solutions requires no Gaussianity approximations. We will highlight this important observation.\n\n- relative quantities for internal noise and contribution of our work (“In the evaluation…”)\n\nWe recomputed the relative importance of internal noise in the sensorimotor task of Section 4.2 to demonstrate the relevance of our results.\n\nUsing Todorov’s solutions, we computed the time-averaged absolute values of the control signal $u_t$ and the state $x_t$ to obtain $\\mathbb{E}[|u|]$ and $\\mathbb{E}[|x|]$. Having multiplicative noise (Table 3), the total average noise magnitude is $NM =C\\mathbb{E}[|u|] + D\\mathbb{E}[|x|] + \\sigma_{\\eta}$. For $\\sigma_{\\eta} = 0.05$ we obtain $\\sigma_{\\eta}/NM  \\approx 0.1$. The performance improvement of the GD algorithm is around 90%. This might be a realistic scenario, so it is important to have good models that produce the actual optimal controller. We thank the reviewer for the relevant suggestion.\n\n- additional work (“Lastly...”)\n\nThis would be good additions for a longer version of the paper.\n\n- Mt and Gt\n\nUpdating $M_t$ and $G_t$ means computing the matrices $M_t$ and $G_t$ with Eqs. (9) and (10), using the values of the needed moments at the time step t, that is $\\mathbb{S}_{old}$.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
- switching linear dynamics

We discuss here how to extend our approach to switching linear dynamics. 
One of the underlying assumptions in this work and in [7], as outlined in Section 2, is that the agent has complete knowledge of the updating rules of the latent dynamical system. By using the same set of matrices to update the state and the state estimate, we implicitly assume that all uncertainty in the estimation process arises solely from noise sources: the problem of inferring the matrices A and B goes beyond the objectives of this approach.
For this reason, to extend our work to the more general and realistic case of switching linear dynamics (SLD), we can consider a matrix A depending on the time step t, $A_t$. A complete formulation of SLD might require adding another variable, a discrete switch variable $s_t$ regulating the way the matrices $A_t$ vary with time and context (see Becker-Ehmck et al, PMLR, 2019). Given that in our case the agent has access to the updating rules of the dynamical system, we can omit $s_t$ (the agent does not have to infer $s_t$ and $A_t$) and directly consider the case in which we have a predetermined set of matrices $A_t, t=1,...,T-1$. The same applies to the matrix B, that can be replaced by $B_t, t=1,...,T-1$. Note that to preserve linearity we assume $A_t$ and $B_t$ to be independent on $x$ and $\hat{x}$.
We consider here the multidimensional case to be as general as possible. 
To extend the GD algorithm we only need to modify the block matrix $M_t$ that we use to update the moments $ \Sigma_t, \mathbb{S}_t, \mu_t $ and eventually propagate the expected cost $\mathbb{E}[J]$ through Eq. 6. Indeed, once we can compute the expected cost at fixed control and filter gains, $L_t, t = 1,...,T-1$, and $K_t, t = 1,...,T-2$, we can use Algorithm 1 to define the objective function to be minimized through Gradient Descent with respect to $L_t$ and $K_t$. To update the block matrix $M_t$ we have to substitute A and B respectively with $A_t$ and $B_t, t=1,...,T$ in Eq. 9.
To handle the potentially high computational costs of performing a numerical Gradient Descent, we introduced the analytical counterpart of the GD algorithm, the Iterative-GD. For the one-dimensional case, Section A.4.1, it supports all the noise sources mentioned in Section 2 (additive, multiplicative and internal). We extended this algorithm to the multi-dimensional case for additive and internal noise in Section A.4.2 for the sake of simplicity, leaving the more general version for future work (Section A.4.2 outlines how this can be done). Here, we extend the afore-mentioned Iterative-GD algorithm (for both one-dimensional and multi-dimensional cases) to switching linear dynamics, following a similar procedure to that of the numerical algorithm. 
For the one-dimensional case, we replace A and B respectively with $A_t$ and $B_t$ in Eqs. 34,35,36,50,51,52 and 53.

For the multi-dimensional case we have to substitute A with $A_{t+j}$ and B with $B_{t+j}$ in Eq. 70 and B with $B_t$ in Eq. 71.
Finally, as previously done, we replace A with $A_t$ in Eq. 75 for $\tilde{M}_t$.
With these changes, we can implement Algorithm 2 for the case with switching linear dynamics.


- Fig.1 in the uploaded pdf: unbiasedness

In Fig.1 we show how the condition of unbiasedness, stating $ \mathbb{E}[x_t|\hat{x}_t] = \hat{x}_t $ does not hold under any condition. Therefore, for the unbiasedness to break down, one does not need internal noise. The idea is simple: in Fig. 1, Section 3.3, in the manuscript, we intuitively explain that a large fluctuation in internal noise can lead to a large $\hat{x}_t$, while $x_t$ does not change. This causes the the unbiasedness condition $\mathbb{E}[x_t|\hat{x}_t] = \hat{x}_t$ to break. But the point is that the same would happen, without internal noise, if there is a large fluctuation of sensory noise: in this case again, $\hat{x}_t$ would be large without changing $x_t$.

To generate the figures, we considered a 1D problem, with the same system dynamics as the one in Section 4.1 (Table 2). We set $\sigma_{\omega} = \sigma_{\xi} = \sigma_{\epsilon} = \sigma_{\rho} = 0.5$ and varied $\sigma_{\eta} = 0.0, 0.3, 0.6$. The initial conditions are $x_0 = \hat{x}_0 = 1$. We considered $ T=10 $ and selected $t=8$ to compute $ \mathbb{E}[x_t|\hat{x}_t] $ as a function of $\hat{x}_t$. Note that the choice of t is arbitrary. We collected the list of states and state estimates at t=8 over 5 *10^7 trials. We then binned the data for $\hat{x}_t$ using $\delta \hat{x} = 0.1$ for the size of the bins. To obtain $ \mathbb{E}[x_t|\hat{x}_t] $  we then computed the mean of all the $x_t$ falling in the same bin, and we used the std for the error bars.

– Fig.2 in the uploaded pdf: high-dimensional task

In Fig.2 we show how our algorithm scales to high-dimensional problems. All the parameters of the dynamical system, that follows Eqs. (1)-(5), are detailed in the caption of Fig. 2. To compute the expected cost, we used Algorithm 1 (but the results are confirmed by Monte Carlo simulations). The pseudo-determinat $|L|_p$ in Fig. 2b is computed by multiplying all the non-zero singular values of the matrices $L_t$ and then averaging this quantity over $t=1,...,T-1$.

Meta-review: ## Recommendation

I recommend "Accept (poster)".

## Summary

The manuscript presents a novel framework on stochastic optimal control where the noise is modeled as multiplicative. The authors present their framework along with 2 gradient-based methods. Interestingly, one of the methods is practical and can be efficiently implemented. The paper includes representative experiments.

## Concerns

I have two main concerns regarding the manuscript:

1) The paper is difficult to read for people without background in stochastic optimal control. The authors should make an effort to enhance readability for the NeurIPS audience (NeurIPS is not a control conference).
2) The evaluation is only performed in "simple" (yet representative) systems/problems. It is important for the NeurIPS community to at least have a detailed discussion (in the text) about the scaling of the proposed method to higher dimensional problems.

## Minor comment

In my opinion, the title of the manuscript is not good. A suggestion is: "Stochastic Optimal Control and Estimation with Multiplicative Noise". Maybe my suggestion is not good, but the current title reads like a sentence from the text.

## Verdict

Despite my concerns, I recommend acceptance (poster) for the following reasons:

1) The authors' responses were clear, concise, and offered valuable discussions. Overall, their comments reflect a strong grasp of the literature.
2) All reviewers had a positive opinion of the work.
Decision: Accept (poster)