Graphical Newton

arXiv:1508.00952v3 [math.OC] 8 Oct 2017

Akshay Srinivasan1 and Emanuel Todorov2

Abstract— Computing the Newton step for a generic function
f : RN → R takes O(N 3 ) flops. In this paper, we explore
avenues for reducing this bound, when the computational
structure of f is known beforehand. It is shown that the
Newton step can be computed in time, linear in the size of
the computational-graph, and cubic in its tree-width.

Newton’s method is essential to many areas of at the
core of second-order methods in nonlinear-optimization. It’s
applicability to large-scale programming, however, is often
limited due to the run-time complexity in computing the
Newton step.
For a generic function f : RN → R, computing the
Hessian requires atleast O(N 2 ) flops; further inverting the
matrix requires O(N γ ) flops (γ = 3, in practice). This is
computationally infeasible for many problems in practice.
Often, however, one is also given access to the the computational structure of the objective. The computer routine
for calculating the objective f (·) can be represented as a
Directed Acyclic Graph [DAG] mapping inputs to f (·) via
intermediary nodes.
For instance, the objective function for the canonical
optimal-control problem is given by,
"
#
n−1
X
min
J (u0 , . . . , un ) ,
li (xi , ui ) + ln (xn ) ,
i=0

∀i, xi+1 ← f(xi , ui ),
(1)
where the dynamics and local-objectives of the system are
given by f(·, ·), and li (·, ·) respectively. The infix operator
’←’ indicates that the value appearing on the right-hand
side, is given the placeholder symbol present to its left; we
explicitly distinguish this from the ’=’ operator, which is
taken to represent a constraint.

u0

J (u0 , . . . , un ) =l0 (x0 , u0 )+
l1 (f(x0 , u0 ), u1 )+
l2 (f(f(x0 , u0 ), u1 ), u2 ) + . . . .

I. I NTRODUCTION

u0 ,u1 ,...,un−1

the apparent sparsity in (1), is entirely destroyed once all the
placeholders are substituted for,

x1

x2

x3

xn−1

u1

u2

u3

un−1

xn

Fig. 1.
Optimal control problem: The dynamical system states are
represented by {xi }, and the control by nodes {ui }.

The order of computation for the objective (1) can be
represented by a linear-chain (Figure I). Lacking constraints,
1 This work was done when Akshay Srinivasan was at the University of
Washington akshaysrinivasan@gmail.com
2 Emanuel Todorov is with the Depts. of Computer Science & Engineering and Applied Mathematics, University of Washington, Seattle, WA
98195, USA todorov@cs.washington.edu

The Hessian of J (·) thus being dense, implies a run-time
that is cubic in the input dimensions for the Newton step
computation; computing the Hessian itself is quadratic.
By contrast, once the problem (1) is written in its constrained form (by replacing ’←’ with ’=’), the sparsity of the
resulting Karush-Kuhn-Tucker [KKT] system, readily allows
for computing the SQP/Lagrange-Newton step in linear time
[1]. Such a transformation, however, comes at the cost of
increasing the size of the optimization problem, abandoning
state feasibility, and increased implementation complexity.
The question which this paper answers, is whether there
exist general techniques, which allow exploiting the sparsity
of the problem, while working solely with the input variables.
Note that these are not questions merely about elimination
orders, but are also verily algebraic in nature.
AUTOMATIC D IFFERENTIATION: Research on Automatic
Differentiation [AD] has produced many techniques for
exploiting the computational structure of generic functions.
They are routinely employed for efficient calculation of gradients and Hessian vector products [2]. The applicability of
AD to second-order optimization is, however, quite limited.
AD is typically used either for computing the entire
Hessian matrix, or for calculating Hessian vector products for
use in Nonlinear Conjugate Gradient Descent [CG]. Hessians
are computed, by accumulating one column at a time via calls
to the Hessian vector product routine [2]. The sparsity of the
Hessian can be exploited in reducing the number of such calls
[3], but structured problems such as (1) will not allow for
any such economy. Compositional chains of functions, such
as those in optimal control (1) (Figure I), not only serve to
make Hessians dense but can also lead to condition numbers,
exponential in their diameters. Large condition numbers are
likely to negate any computational advantages offered by
methods like CG.
The above techniques form the Hessian matrix, directly or
indirectly, before computing the Newton step. This stands in
contrast with the root-finding problem, for which there do
exist methods for directly computing the Newton-Raphson
step [4] [2] [5]. The root finding problem involves the
inversion of the Jacobian of a function - rather than the
Hessian - and these methods reduce this computation to
that of inverting a sparse matrix [4]. The Newton step (for

optimization) can also be computed using this method by
formulating it as a root-finding problem on the gradient.
This, however, results in non-symmetric matrices depending
on the computational graph of the gradient, as opposed to
the function itself. The latter graph is transitively closed,
and hence, analyses for the above Newton-Raphson AD
algorithm apply to the gradient, and are difficult to extend
to the underlying objective [5].
DYNAMIC P ROGRAMMING: The question posed earlier,
has already been answered in the affirmitive, for the optimal
control problem. There exists an algorithm for optimal control, based on Dynamic Programming, that exploits algebraic
dependencies in (1), in order to compute the Newton-step in
only linear time [6] [7] [8].
The run-time of this algorithm is the direct result of
the sparsity of the corresponding constrained problem [6]
[1][9]. The band-structure of the relevant KKT system allows
for solving the system in linear time [1]. The relationship
between computing the Newton-step (Hessian of the objective), and computing the Lagrange-Newton step (Hessian of
the Lagrangian), is established by noting that there exist
multiplier values such that both compute the same result [6].
Such algorithms are routinely employed by practitioners
for updating control policies in real-time, while maintaining
a feasible trajectory. These algorithms have been extended to
Extended Kalman Filtering (EKF), as well as various other
formulations of the control problem [10] [11] [12].
OVERVIEW: We generalize such algorithms, by using
Hessian vector product equations from AD, to relate the
computation of Newton step and Lagrange-Newton step, for
arbitrary structured objectives.
We then extend this framework to structured optimization
problems with equality constraints.
Further, we show that solving the resultant KKT systems
can be accomplished in time Õ(tw3 ), where ’tw’ is the treewidth of the canonical computational graph.
Finally, we show results from numerical experiments.
II. N OTATION
Let G be a Directed Acyclic Graph [DAG], and let each
vertex v ∈ V [G], be associated with state Sv ∈ Uv ⊂ Rnv ,
taking values in an open set. Denote by δ + (v), the parents of
v ∈ V [G], and by δ - (v) its children; let SA be the (labelled)
concatenation of states, associated with vertices in set A ⊂
V [G]. Define the set of input nodes X = {x1 , x2 , . . . , xn } ,
{v | δ + (v) = ∅, v ∈ V [G]}, to be the parentless vertices of
G.
An objective function f : Ux1 × . . . Uxn → R, has the
computational structure given by the tuple (G, {ϕv }, {lv }),
if
Q it can be written as the sum of local objectives lv :
z∈{v}∪δ + (v) Uz → R, on the graph G,
X
f : (Sx1 , . . . , Sxn ) 7→
lv (Sv∪δ+ (v) ),
(2)
v∈V [G]
Sv ← ϕv (Sδ+ (v) ), ∀v ∈ V [G], δ + (v) 6= ∅.
The state of a non-input node v ∈ V [G] in (2), is defined
recursively as Sv ← ϕv (Sδ+ (v) ), for some given function

Q
ϕv :
z∈δ + (v) Uz → Uv . It follows since G is a DAG,
that SV [G] and hence f (·), is uniquely determined from the
input SX , and functions {ϕv }. The order of computation for
the objective is given by the topological ordering of G, and
the DAG G is called the computational graph of f (·). The
computer routine for calculating any objective function, can
be represented by such a structure [2].
In the following sections, the symbolism ∂u v is used as
∂Sv
. The derivatives of functions with
a shorthand for ∂S
u SX
respect to Su are similarly denoted by the operator ∂u ; that
with respect to a (labelled) set A = {v1 , v2 , . . . } ⊂ V [G] by
∂A , [∂a1 , ∂a2 , . . . ].
III. N EWTON STEP
Consider the objective function in (2), defined by the tuple
(G, {ϕv }, {lv }). The optimization problem of interest is the
following,


X
lv (Sv∪δ+ (v) ) ,
min f ,
Sx1 ,...,Sxn
(3)
v∈V [G]
Sv ← ϕv (Sδ+ (v) ), ∀v ∈ V [G], δ + (v) 6= ∅,
and the corresponding constrained problem is obtained by
replacing the operator ’←’ by ’=’ in (3).
In the following, we consider first the constrained formulation of (3), and define the KKT system involved in computing
the Lagrange-Newton step; we then relate these to computing
the Newton step.
A. Lagrange-Newton
The Lagrangian for the constrained form of (3), is given
by,
L(SV [G] , λ) ,
X
lv (Sv∪δ+ (v) ) +
v∈V [G]

X

λTv hv (Sv∪δ+ (v) ),

v∈V [G],
δ + (v)6=∅

where,
∀v ∈ V [G], δ + (v) 6= ∅,

hv (Sv∪δ+ (v) ) , ϕv (Sδ+ (v) ) − Sv ,
(4)
and the vector λ is the labelled concatenation of all λv ’s.
The necessary first order conditions for optimality of this
problem are given by [13],
∂V L(S∗V , λ∗ ) = 0,

h(S∗V ) = 0.

(5)

The Lagrange-Newton step for solving this system of equations, around a nominal (S˜V , λ), entails solving the following
KKT system [13],
 2

 

∂V L ∂V hT
δSV
−∂V L
=
.
(6)
∂V h
0
δλ
−h
Sequential Quadratic Programming [SQP], involves taking
a step along (δSV , δλ) and iteratively solving for the first
order conditions. In the following section, it will be shown
that there exist values for Lagrange multipliers, depending
only on the inputs, such that the solution to (6), yields the
Newton step for the unconstrained objective.

B. Unconstrained Newton
We recollect certain defintions from AD, and then continue
to present one of the central results of the paper.
R EVERSE AD: The first derivatives of the objective f (·)
can be calculated by applying the chain rule over G,
X
X
∀v, ∂v f =
∂v ls +
∂d f T ∂v d;
s∈v∪δ - (v)

δ[∂v f ] =

X

s∈v∪δ - (v)


X

2
∂av
ls · δSa  +

a∈v∪δ + (s)


X



δ[∂d f ]T ∂v d +

d∈δ - (v)

∀a,

X

2
(∂d f T ∂av
d) · δSa ) ;

a∈δ + (d)

X

δSa =

(8)
These equations can be solved, for a given δSX , by a
forward-backward recursion similar to the one used for
solving (7) [2]. Computing the Hessian-vector product in this
manner takes time Õ(ω(Ĝ)2 ) 1 [2], where ω(Ĝ) is the clique
number of the moralization of G.
N EWTON STEP: The problem of interest is, however, the
exact inverse: find a δSX , such that δ[∂X f ] = −∂X f . This
question is answered by the following theorem.
Theorem 1: (Newton step) The Newton step for the objective (2) is given by the Lagrange-Newton step (6), when
SV is feasible and when ∀v, λv = ∂v f as defined in (7).
Proof. The second equation in (8) is equivalent to ∂V h ·
δSV = −h, in (6). Rearranging the first equation from (8),
and setting δ[∂v f ] = −∂v f for all inputs, we obtain ∀v,
0=
s∈v∪δ - (v),
a∈v∪δ + (s)

−

1 We

s∈v∪δ - (v),
a∈v∪δ + (s)

−

X

2
∂va
ls δSa +

(
δλv
0

2
(λTd ∂av
d) δSa +

d∈δ - (v),
a∈δ + (d)

δ + (v) 6= ∅
otherwise

(10)

!
X

+

(∂v d)T δλd ,

d∈δ - (v)

where,
∂v L =
P
P

 P
s∈v∪δ - (v) ∂v ls +
d∈δ - (v) λd · ∂v d,
P
s∈v∪δ - (v) ∂v ls +
d∈δ - (v) λd · ∂v d

,

−λv

δ + (v) = ∅
otherwise
(11)

The result follows from equations (7), (9), (10) & (11).

G RAPHICAL N EWTON: The above theorem immediately
yields the following optimization algorithm,
Algorithm 1 Graphical Newton
1: Input: initial SX , tuple (G, {ϕv }, {lv })
2: repeat
3:
Compute f , {∂v f }, {∂ 2 ϕv } from (2), (7).
4:
Compute the SQP step from (6), with λv = ∂v f , ∀v.
5:
Compute step-length η via linesearch on inputs SX .
6:
Update inputs: SX ← SX + ηδSX .
7: until ||∂X f || ≤ 

∂d a · δSd .

d∈δ + (a)

X

− ∂v L =
X

d∈δ - (v)

(7)
∂ϕd (Sδ+ (d) )
v ∈ δ + (d) ⇒ ∂v d ,
.
∂Sv
Since G is a DAG, there exist child-less nodes (i.e δ - (v) =
∅), from which the above recursion can be initialized. The
recursion then proceeds backward in the depth first search
order on G. This algorithm is known as reverse-mode AD
[2].
H ESSIAN VECTOR AD: A change in the inputs δSX ,
results in the first-order change in the derivative, δ[∂v f ] ,
2
∂Xv
f · δSX , which is given by the Hessian vector product.
Computing the Newton step is thus, equivalent to finding a
δSX such that, δ[∂X f ] = −∂X f .
Applying chain-rule over the DAG G, for all terms in (7),
we obtain,
∀v,

Similarly, expanding the top block in (6) using the definitions
in (3) & (6), we obtain ∀v,

2
∂va
ls δSa +

X

2
(∂d f T ∂av
d) δSa +

d∈δ - (v),
a∈δ + (d)

(
!
X
δ[∂v f ] δ + (v) 6= ∅
+
(∂v d)T δ[∂d f ].
−∂v f otherwise
d∈δ - (v)
(9)

use Õ(.) to hide factors linear in |E| + |V |.

The run-time of every iteration in Algorithm 1 depends
crucially upon the time required to solve (6). The run-time
bounds for solving such KKT systems is taken up later in
the paper.
C. Extension to equality constraints
Consider optimization problems, which have equality constraints in addition to the structured objective from before,

min

Sx1 ,...,Sxn

f ,


X

lv (Sv∪δ+ (v) ) ,

v∈V [G]

(12)

Sv ← ϕv (Sδ+ (v) ), ∀v ∈ V [G], δ + (v) 6= ∅,
c(SC ) = 0,
where c(·) = 0 is an additional equality constraint, which
depends on the variables C ⊂ V [G]. The Lagrangian for this
problem is given by,
L̂(SV [G] , λ) = L(SV [G] , λV \X ) + λTh c(SC ),

(13)

where L is as defined in (13), and λV \X is the corresponding
set of multipliers; the variable λ, being the concatenation of
λc and all multipliers, λV \X , appearing in (13).
Theorem 1 can be applied to this problem by treating
λTc c(SC ) as another cost function in the objective, while also

including the constraint in the KKT system (6). The iteration
can then proceed by solving the KKT system with λv =
∂v (f +λTc c), ∀v, and using a merit function for the linesearch
procedure; the variables (SX , λc ) are updated accordingly.
We omit the proof for the validity of this method.
IV. M ESSAGE PASSING
The classical run-time bound for Cholesky factorization
(i.e Gaussian Belief Propagation 2 ) [14] [15], cannot be extended to problems such as (6), because of the appearance of
linear constraints. Such bounds for structured KKT systems,
do not appear to be known within the sparse linear algebra
community [16].
In this section, we provide a Message Passing algorithm
for solving such KKT systems, and show that it has a runtime bound of Õ(tw3 ) 3 , given the tree-decomposition.
A. Hypergraph structured QPs
For a hypergraph H, denote the adjacency and incidence
matrices by A[H] & B[H] respectively,
A[H] ∈ R|V [H]|×|V [H]| ,
B[H] ∈ R|E[H]|×|V [H]| ,
(
1 ∃e ∈ E[H], u, v ∈ e
A[H]uv =
0 otherwise
(14)
(
1 u∈e
B[H]eu =
0 otherwise
Given such a hypergraph H, the family of QPs we’re
interested in solving is the following,
X 1
ST Qe Se − bTe Se ,
min
x
2 e
(15)
e∈E[H]
∀e ∈ E[H],

Ge Se = he .

Assuming that the QP has a bounded solution and that the
constraints are full rank, the minimizer to (15) is given by
the solution to the following KKT system,


 

Q GT
x
b
=
,
G 0
λ
h
(16)
x, b ∈ R|V | ,

λ, h ∈ RM ,

where Q, G, λ, x, b are concatenation of terms defined in (15)
respectively. The sparsity/support of (16) is closely related
to H, since,
supp(Q) ⊆ supp(A[H]),
∀i, ∃e, supp(Gi,: ) ⊆ supp(B[H]e,: ).
Every row of the constraint, Gi,: , has the same sparsity as
some edge e ∈ E[H].
T REE DECOMPOSITION: Extending the notion of Dynamic
Programming to non-trees (including Hypergraphs) requires
a partitioning of the graph so as to satisfy a lifted notion of
being a tree [17]. Tree decomposition captures the essence
of such graph partitions,
2 Gaussian-BP,
3 The

computes the LU decomposition of a matrix
tilde hides factors linear in |V [H]|, |E[H]|.

Definition 1: (Tree decomposition) A tree-decomposition
of a hypergraph H consists of a tree T and a map χ :
V [T ] → 2V [H] , such that,
i (Vertex cover) ∪i∈V [T ] χ(i) = V [H].
ii (Edge cover) ∀e ∈ E[H], ∃i ∈ V [T ], e ⊂ χ(i).
iii (Induced sub-tree) ∀u ∈ V [H], Tu , T [{i ∈ V [T ]|u ∈
χ(i)}] is a non-empty subtree
The tree-width of a tree-decomposition T is defined to be
tw(T ) = maxv∈V [T ] |χ(v)| − 1. The tree-width of a graph
H is defined to be the minimal tree-width attained by any
tree-decomposition of H.
We define the vertex-induced subgraph in what follows to
be H[S] , (V [H], {e∩S, e ∈ E[H]}). The following lemma
ensures that such a decomposition ensures local dependence
[17].
Lemma 1: (Edge separation) Deleting the edge xy ∈
E[T ], renders H[V \(χ(x) ∩ χ(y))] disconnected.
H YPERTREE STRUCTURED QP: The tree-decomposition
itself can be considered a Hypergraph, (V [H], {χ(u), ∀u ∈
V [T ]}). Such a Hypertree4 can also be thought of as a
Chordal graph [15]. We assume henceforth that the given
graph H is a hypertree, and that T is its tree-decomposition.
Algorithm 2 Graphical QP
1: Given: T , H, {Qe }, {be }, {Ge }, {he }.
2:
3:
4:
5:
6:

function GatherMessage(l, p, T )
(Q̃l , b̃l , G̃l , h̃l ) ← (Ql , bl , Gl , hl )
for c ∈ δ T (l)\p do
(Qc→l , Gc→l , bc→l , hc→l ) ← GatherMessage(c, p, T )

(Q̃l , b˜l ) ← (Q̃l , b˜l ) + (Qc→l , bc→l )
G̃l ← [G̃l ; Gc→l ], h̃l ← [h̃l ; hc→l ]
9: end for
10: return Factorize(χ(l), χ(p), Q̃l , b̃l , G̃l , h̃l )
7:
8:

11:
12:
13:
14:

function Factorize(χ(l), χ(p), Q̃, b̃, G̃, h̃)
(ξ, ι) ← (χ(l)\χ(p), χ(l) ∩ χ(p))
r ← rank(Q̃ι,ι )
15: return Gaussian-BP messages from (17).
16: return
The gather stage of the Message Passing algorithm, is
illustrated in Algorithm 2. 5
The function, Factorize, computes the partial LU decomposition of its arguments; we describe below, its operation.
Denote the vertices that are interior to l by ι = χ(l) ∩ χ(p),
and those on the boundary (i.e common to p, l) by ξ =
χ(l)\χ(p), and let r = rank(Q̃ι,ι ). The function computes
Gaussian-BP messages from block pivots {2, 3} to {1, 4} in
(17). Note that, unlike Gaussian-BP, the matrices in (17) are
not necessarily positive definite, but are however invertible.
4 There are multiple definitions of a Hypertree; we use the term to mean a
maximal Hypergraph, whose tree-decomposition can be expressed in terms
of its edges.
5 Note that the addition is performed vertex label-wise in Line 6 of
Algorithm 2.

Fig. 2. LEFT: An optimal limit cycle for the system ẍ = −(x3 + ẋ3 )/6 + u. MIDDLE: Convergence of the objective function for the limit-cycle problem
(20). RIGHT: Convergence in norm, of the Lagrangian gradient, and constraint deviation.


T
T
T
 Q̃ξξ Q̃ιξ G̃:r,ξ G̃r:,ι


 Q̃ιξ Q̃ιι G̃T:r,ι G̃Tr:,ι


 G̃:r,ξ G̃:r,ι
0
0


G̃r:,ξ G̃r:,ι 0
0













Sξ
Sι
λ:r
λr:





 b̃ξ
 
  b̃
= ι
 
  h̃
  :r

h̃r:












solver [19]. For ensuring convergence in constrained problems, an augmented Lagragian merit function was used
[20]. The implementation was tested on the following nonstandard control problems.

(17)
Gaussian Belief-Propagation is essentially a restatement of
LU decomposition [14]. Gaussian-BP consists of messages
of the form [18] [15],
X
−1
µi→j := [Ji→j , hi→j ] = [Jii , hi ] −
Jik Jk→i
[Jki , hk→i ],
k∈δ(i)\j

µi =

−1
Ji→j
(hi→j

− Jij µj ),

(18)
where Jµ = h is the equation that is to be solved. These
can be replaced by appropriate square-root forms to obtain
instead, an LDL decomposition.
Theorem 2: The linear equation (16) can be solved in
time Õ(tw(H)3 ), given the minimal tree-decomposition via
Algorithm 2.
Proof. The correctness of the algorithm follows from
Lemma 1. The bound holds trivially if, rank G̃ ≤ rank Q̃ι,ι ,
at every step of the algorithm. Otherwise, by realizing that
G̃l→p , can’t have rank more than |χ(p)|, the proof follows.

It follows from Theorem 2, that the KKT system in
Algorithm 1 can be solved in time Õ(tw(Ĝ)3 ), where Ĝ is
the moralization of the computational graph G.
The above proof also ensures that the equivalent sparse
LU/LDL decomposition [14], with the same pivot order, also
has the same run-time. Since decompositions of indefinite
systems are subject to instability, use of specialized solvers
is generally preferable.
V. N UMERICAL E XPERIMENTS
In this section, we present preliminary numerical results
with an implementation of Algorithm 1, using the MA57

Fig. 3. The computational graph for the limit cycle problem, with N = 20

S PRING - DAMPER LIMIT CYCLE: Consider the following
spring-damper limit-cycle problem [11] (Figure V),
Z
min
`(x, ẋ, u)dt,
x0 ,u[0,T ]

ẍ = −(x3 + ẋ3 )/6 + u,

(19)

x(0) = x(T ) = x0 , ẋ(0) = ẋ(T ) = ẋ0 ,
where,
2
2
1
`(x, ẋ, u) = (1 − e(ẋi −2) − e−(ẋi +2) ) + ||ui ||22 .
2
Discretising the derivatives by finite differences, ẋ ≈
∆xi /∆t = (xi − xi−1 )/∆t, this can be written as the
following structured optimization problem,

min m

x0 ,{ui }0

N
X

`(xi , ∆xi /∆t, ui ),

1



xi+1 ← xi + ∆xi + (∆t)2 −(x3 + (∆xi /∆t)3 )/6 + ui .
x0 = xN −2 , x1 = xN −1 .
(20)
For N = 100, ∆t = 0.1, with random initializations, the
problem showed robust convergence; often taking no more

than ten SQP iterations. The optimal limit cycle, and the
convergence curves for one run of the algorithm are shown
in (Figure IV-A).
ACROBOT: We also considered the Acrobot problem,
using a discretized second-order dynamics as before, and
`2 and total-variation penaltes over the control sequence;
the end-condition being enforced using a constraint over the
final cartesian position. For N = 100, ∆t = 0.04d0, with
random initializations, the optimizer converged often in about
20 iterations (Figure V).

the Hessian vector product takes time Õ(ω(Ĝ)2 ), where Ĝ is
a moralization of the computational graph [5]. By contrast,
if the computational graph were chordal, then computing
the Newton-step via Algorithm 1 is only Õ(ω(Ĝ)3 ). The
latter is more economical when the cliques of a graph are
small in comparison to the order of the graph. The illconditioned nature of structured objectives may also lead to
bad convergence properties for CG algorithms.
For problems whose tree-widths are large, the iterative
method is obviously more viable. However, following the
rapid advances in approximate inference in the past two
decades [15], we hope that the explicit algebraic connection
to graphical models made in this paper, can be exploited in
coming up with less-agnostic iterative methods.
R EFERENCES

Fig. 4. TOP: Convergence of the objective function for the Acrobot problem.
BOTTOM : Convergence in norm, of the Lagrangian gradient, and constraint
deviation.

VI. D ISCUSSION
We have shown that the Newton step can be computed in
time Õ(tw3 ), where ’tw’ is the tree-width of the computational graph. We have also derived extensions to constrained
problems, and provided numerical examples. The technique
presented herein, also generalizes many specialized algorithms in control.
In certain control problems, the solution to the KKT
system, itself can be written in feedback form. Given a
LU decomposition of the KKT system, one can replace
the backsubstitution phase by U , with a function evaluation
that uses L as a control feedback [7]. It is unclear if such
techniques can be generalized, and whether they can be made
independent of the pivot-order used for solving the system.
A competing method for exploiting the structure of objectives such as (3), is by the use Hessian vector product AD
routines in conjugation with CG-like methods. Computing

[1] S. J. Wright, “Solution of discrete-time optimal control problems on
parallel computers,” Parallel Computing, vol. 16, no. 2, pp. 221–237,
1990.
[2] A. Griewank and A. Walther, Evaluating Derivatives: Principles and
Techniques of Algorithmic Differentiation. SIAM, 2008.
[3] T. F. Coleman and J. J. More, “Estimation of sparse Hessian matrices
and graph coloring problems,” Mathematical Programming, vol. 28,
no. 3, pp. 243–270, 1984.
[4] J. Utke, “Efficient Newton steps without Jacobians,” in Computational
Differentiation: Techniques, Applications, and Tools, M. Berz, C. H.
Bischof, G. F. Corliss, and A. Griewank, Eds. Philadelphia, PA:
SIAM, 1996, pp. 253–264.
[5] L. Dixon, “Automatic Differentiation: Calculation of Newton Steps,”
in Encyclopedia of Optimization. Springer, 2009, pp. 137–142.
[6] J. De O. Pantoja, “Differential Dynamic Programming and Newton’s
method,” International Journal of Control, vol. 47, no. 5, pp. 1539–
1553, 1988.
[7] D. H. Jacobson and D. Q. Mayne, Differential Dynamic Programming.
North-Holland, 1970.
[8] E. Todorov and W. Li, “A generalized iterative LQG method for
locally-optimal feedback control of constrained nonlinear stochastic
systems,” in American Control Conference, 2005. Proceedings of the
2005. IEEE, 2005, pp. 300–306.
[9] D. Ralph, “A parallel method for unconstrained discrete-time optimal
control problems,” SIAM Journal on Optimization, vol. 6, no. 2, pp.
488–512, 1996.
[10] M. Toussaint and C. Goerick, “A Bayesian view on motor control and
planning,” in From Motor Learning to Interaction Learning in Robots.
Springer, 2010, pp. 227–252.
[11] Y. Tassa, T. Erez, and E. Todorov, “Optimal limit-cycle control recast
as Bayesian inference,” in Proceedings of thh IFAC world congress.
Citeseer, 2011.
[12] S. J. Wright, “Interior point methods for optimal control of discrete
time systems,” Journal of Optimization Theory and Applications,
vol. 77, no. 1, pp. 161–187, 1993.
[13] J. Nocedal and S. Wright, Numerical Optimization. Springer Science
& Business Media, 2006.
[14] T. A. Davis, Direct Methods for Sparse Linear Systems. SIAM, 2006,
vol. 2.
[15] M. J. Wainwright and M. I. Jordan, “Graphical Models, Exponential
Families, and Variational Inference,” Foundations and Trends R in
Machine Learning, vol. 1, no. 1-2, pp. 1–305, 2008.
[16] R. Bridson, “An ordering method for the direct solution of saddle-point
matrices,” Preprint, 2007.
[17] J. Kleinberg and É. Tardos, Algorithm Design. Pearson, AddisonWesley, 2006.
[18] J. S. Yedidia, W. T. Freeman, and Y. Weiss, “Generalized Belief
Propagation,” in NIPS, vol. 13, 2000, pp. 689–695.
[19] I. S. Duff, “MA57—a code for the solution of sparse symmetric
definite and indefinite systems,” ACM Transactions on Mathematical
Software (TOMS), vol. 30, no. 2, pp. 118–144, 2004.
[20] P. E. Gill, W. Murray, M. A. Saunders, and M. H. Wright, “Some
theoretical properties of an augmented lagrangian merit function.”
DTIC Document, Tech. Rep., 1986.

